--- title: "Introduction to cusna" output: litedown::html_format vignette: > %\VignetteEngine{litedown::vignette} %\VignetteIndexEntry{Introduction to cusna} %\VignetteEncoding{UTF-8} --- `cusna` gives R a native, bit-reproducible engine for stochastic actor-oriented models (SAOM --- the model family of [RSiena](https://cran.r-project.org/package=RSiena)) and exponential random graph models. The engine is a small C interface over C++ host logic, compiled with the package --- **no Python runtime**. The package builds the engine's CPU backend (C++17, multi-threaded with OpenMP); a CUDA backend implements the same interface and is built separately from source. Because the engine ships inside the package, the estimation chunks in this vignette actually run (on small synthetic data, with reduced simulation counts to keep the build fast). ## Setup Nothing to set up --- the engine is already compiled in: ```{r} library(cusna) cusna_abi_version() cusna_has_cuda() # FALSE on the CPU-only build # Keep the vignette build well within CRAN's two-core policy: the native CPU # backend is OpenMP-parallel over simulation chains, so run it single-threaded # here. cusna_set_threads() calls the OpenMP runtime directly, so it takes # effect immediately (unlike Sys.setenv(OMP_NUM_THREADS = ...), whose effect # on an already-initialised OpenMP runtime is platform-dependent). The demo # data is tiny, so a single thread is still fast. old_threads <- cusna_set_threads(1L) ``` ## Specifying a model Model terms are plain R objects: ```{r} effects <- list( cusna_effect("density"), cusna_effect("recip"), cusna_effect("transTrip") ) effects[[1]] ``` Covariate effects take a covariate; the engine centers it for you. A covariate is a numeric vector of length *n* (constant) or an *n*-by-waves matrix (changing): ```{r} alcohol <- c(3, 2, 2, 1, 3) # one value per actor (toy) cusna_effect("egoX", covariate = alcohol) ``` User interactions multiply the change statistics of two or three component effects, which must appear earlier in the list. Component positions are 1-based: ```{r} list( cusna_effect("egoX", covariate = alcohol), cusna_effect("altX", covariate = alcohol), cusna_interaction(c(1, 2)) # egoX x altX ) ``` ## Data and estimation Build a panel from a list of adjacency matrices (missing ties are `NA` or the code 9; structural zeros/ones are 10/11). `saom_data()` applies RSiena's imputation and per-period masks --- a construction validated bit-for-bit against the reference implementation: ```{r} set.seed(7) n <- 20 w1 <- matrix(as.integer(runif(n * n) < 0.12), n, n); diag(w1) <- 0L w2 <- w1; flip <- sample(n * n, 40); w2[flip] <- 1L - w2[flip]; diag(w2) <- 0L w3 <- w2; flip <- sample(n * n, 40); w3[flip] <- 1L - w3[flip]; diag(w3) <- 0L dat <- saom_data(list(w1, w2, w3)) dat ``` Estimation runs the RSiena-style Robbins--Monro method of moments on the native simulator. (The tiny `mom_control()` counts below keep this vignette quick; the defaults follow RSiena and are the ones to use in practice.) ```{r} fit <- mom_estimate( dat, effects = list(cusna_effect("density"), cusna_effect("recip")), conditional = TRUE, # RSiena's default for one network control = mom_control(n1 = 100, nsub = 1, n2 = 10, batch2 = 50, n3 = 200)) fit ``` The result is an ordinary R object: ```{r} summary(fit) # RSiena-style table (rates, estimates, s.e., t-conv) coef(fit) # named parameter vector round(vcov(fit), 4) # delta-method covariance ``` ### Covariates and behavior co-evolution Actor covariates enter through the covariate effects; network--behavior co-evolution adds a `behavior` matrix and behavior effects, and uses unconditional estimation (as in RSiena for more than one dependent variable): ```{r, eval=FALSE} dat <- saom_data(s50$waves, behavior = s50$alcohol) # the bundled s50 panel fit <- mom_estimate( dat, effects = list(cusna_effect("density"), cusna_effect("recip"), cusna_effect("transTrip"), cusna_effect("altX", dyn = TRUE), # covariate = the behavior cusna_effect("egoX", dyn = TRUE)), beh_effects = list(cusna_beh_effect("linear"), cusna_beh_effect("quad"), cusna_beh_effect("avSim")), conditional = FALSE) ``` ### Simulated networks and goodness of fit With `return_nets = TRUE` the fit keeps the end networks of its phase-3 simulations (RSiena's `returnDeps`): one chains-by-n-by-n array per period, the layout the native goodness-of-fit primitive `cusna_gof_distribution()` takes. Comparing the simulated and observed indegree distributions of the last wave: ```{r} fit_gof <- mom_estimate( dat, effects = list(cusna_effect("density"), cusna_effect("recip")), conditional = TRUE, return_nets = TRUE, control = mom_control(n1 = 100, nsub = 1, n2 = 10, batch2 = 50, n3 = 200)) sims <- fit_gof$sims[[2]] # period 2 ends at wave 3 dim(sims) sim_indeg <- cusna_gof_distribution(sims, kind = 0L, levls = 0:4) obs_indeg <- cusna_gof_distribution(array(w3, c(1, n, n)), kind = 0L, levls = 0:4) rbind(observed = obs_indeg[1, ], simulated = colMeans(sim_indeg)) ``` Several dependent networks over one actor set co-evolve through `saom_multinet_data()` + `mom_estimate_multinet()`, with cross-network effects (`"crprod"`, `"from"`, ...) pointing at their source network via `net_ref`. ## Tuning the estimator `mom_control()` exposes the Robbins--Monro schedule. The defaults follow RSiena; raise `diagonalize` towards 1 and add phase-3 simulations for harder, nearly collinear specifications: ```{r} str(mom_control()) ``` | Argument | Meaning | |---|---| | `firstg` | initial gain | | `nsub`, `n2` | phase-2 subphases and their iteration counts | | `batch2` | simulation batch per phase-2 iteration | | `n1`, `nD`, `n3` | phase-1 / derivative / phase-3 simulation counts | | `diagonalize` | shrink the derivative towards its diagonal (0--1) | | `seed` | random seed | On nearly collinear pairs (e.g. `transTrip` + `gwesp`) plain Robbins--Monro may not converge from a cold start: increase `diagonalize`, or drive RSiena's own estimator through the native simulator with `cusna_fran()` --- see `vignette("siena07-backend")`. ## Beyond SAOM: ERGM, TERGM, STERGM, ALAAM The same engine serves cross-sectional and temporal exponential random graph models, and autologistic actor attribute models. Sufficient statistics and a TNT sampler are direct calls: ```{r} x <- w1 # any 0/1 adjacency terms <- list(ergm_term("edges"), ergm_term("mutual")) ergm_stats(x, terms, directed = TRUE) ergm_simulate(x, coef = c(-2, 1), terms = terms, nsim = 5, directed = TRUE, seed = 1) ``` Estimators (each validated against its R reference): ```{r, eval=FALSE} ergm_mcmle(x, terms, directed = TRUE) # MCMC-MLE (matches ergm::ergm) ergm_mple(x) # pseudo-likelihood demo tergm_mple(list(w1, w2, w3)) # TERGM (matches btergm) tergm_simulate(w2, lag = w1, coef = c(-4, 3, 2, 0.5)) stergm_cmle(list(w1, w2, w3), # separable TERGM (matches tergm CMLE) formation = list(ergm_term("edges"), ergm_term("mutual"))) alaam_mple(y, net = x) # ALAAM MPLE (matches glm) alaam_mcmle(y, net = x) # ALAAM MCMC-MLE ``` ## Validation The native stack is validated in layers: the CUDA kernels are bit-identical to the reference GPU implementation; the host statistics reproduce RSiena moment targets on public datasets (s50, Knecht, the COW alliance panel) to machine precision; the R-side data preparation, effect preprocessing, and moment targets are bit-identical to the reference on seven datasets; and the full Robbins--Monro estimator agrees with the reference estimator within simulation standard errors on baseline, covariate, co-evolution, and multi-network models. The ERGM samplers match `ergm::simulate` distributionally, and the ERGM/TERGM/STERGM/ALAAM estimators match `ergm`/`btergm`/`tergm`/`glm` on their benchmark models. No GPU is needed: the package compiles from source on any C++17 toolchain, and `cusna_has_cuda()` returns `FALSE` for it.