cusna gives R a native, bit-reproducible engine for stochastic
actor-oriented models (SAOM — the model family of
RSiena) and exponential random
graph models. The engine is a C ABI over CUDA kernels and C++ host logic,
compiled with the package — no Python runtime. Following the design of
the torch package, the CRAN
build is CPU-only from source; the GPU path is compiled from source when a
CUDA 12.x toolkit is detected at configure time.
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).
Nothing to set up — the engine is already compiled in:
library(cusna)
cusna_abi_version()
#> [1] 1
cusna_has_cuda() # FALSE on the CPU-only build
#> [1] FALSE
# 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)
Model terms are plain R objects:
effects <- list(
cusna_effect("density"),
cusna_effect("recip"),
cusna_effect("transTrip")
)
effects[[1]]
#> <cusna effect> density
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):
alcohol <- c(3, 2, 2, 1, 3) # one value per actor (toy)
cusna_effect("egoX", covariate = alcohol)
#> <cusna effect> egoX + covariate
User interactions multiply the change statistics of two or three component effects, which must appear earlier in the list. Component positions are 1-based:
list(
cusna_effect("egoX", covariate = alcohol),
cusna_effect("altX", covariate = alcohol),
cusna_interaction(c(1, 2)) # egoX x altX
)
#> [[1]]
#> <cusna effect> egoX + covariate
#>
#> [[2]]
#> <cusna effect> altX + covariate
#>
#> [[3]]
#> <cusna interaction> components 1 x 2
#>
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:
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
#> <cusna SAOM data> 20 actors, 3 waves (2 periods)
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.)
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
#> cusna MoM fit (cpu backend, conditional)
#> effect estimate se t_conv
#> density -0.2864 0.1519 -1.037
#> recip 0.0535 0.3201 -0.171
#> overall max. convergence ratio: 1.1289
#> 1800 simulations, 0.30s
The result is an ordinary R object:
summary(fit) # RSiena-style table (rates, estimates, s.e., t-conv)
#> rate period 1 2.4456 (sd of durations 0.4700)
#> rate period 2 2.3418 (sd of durations 0.3770)
#> parameter estimate s.e. t-conv
#> density -0.2864 0.1519 -1.037
#> recip 0.0535 0.3201 -0.171
#> overall max. convergence ratio: 1.1289
#> total simulations: 1800, wall time: 0.30s, backend: cpu, conditional: TRUE
coef(fit) # named parameter vector
#> density recip
#> -0.2864204 0.0535384
round(vcov(fit), 4) # delta-method covariance
| density | recip | |
|---|---|---|
| density | 0.023 | -0.023 |
| recip | -0.023 | 0.102 |
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):
dat <- saom_data(list(w1, w2, w3), behavior = drink) # drink: n-by-waves
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)
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.
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:
str(mom_control())
#> List of 9
#> $ firstg : num 0.4
#> $ nsub : int 4
#> $ n2 : num [1:4] 40 40 60 80
#> $ batch2 : int 256
#> $ n1 : int 1000
#> $ n3 : int 4000
#> $ nD : int 1000
#> $ diagonalize: num 0.2
#> $ seed : int 1234
#> - attr(*, "class")= chr "cusna_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").
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:
x <- w1 # any 0/1 adjacency
terms <- list(ergm_term("edges"), ergm_term("mutual"))
ergm_stats(x, terms, directed = TRUE)
#> edges mutual
#> 37 5
ergm_simulate(x, coef = c(-2, 1), terms = terms, nsim = 5,
directed = TRUE, seed = 1)
| edges | mutual |
|---|---|
| 48 | 6 |
| 55 | 12 |
| 47 | 3 |
| 57 | 7 |
| 60 | 6 |
Estimators (each validated against its R reference):
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
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.
An NVIDIA GPU is optional throughout: the CPU build compiles from source on
any C++17 toolchain, and cusna_has_cuda() reports which build you have.