Simulates data from a Gaussian Graphical Model (GGM) with a stochastic
block model (SBM) structure. Nodes are partitioned into communities, with
edge density controlled separately within and between communities via a
blocks x blocks density matrix. Absolute edge weights are drawn from
a Weibull distribution whose parameters are predicted from the network size,
sample size, and signal-to-noise ratio using a Seemingly Unrelated
Regression (SUR) model fitted to 194 empirical psychometric networks
(Huth et al., 2025). The resulting population partial correlation matrix is
used to generate multivariate normal (or skewed) data.
A diffusion parameter controls the minimum proportion of the
strongest edges that are reserved for within-community positions. Because
the remaining edges are shuffled randomly, the effective within-community
advantage will generally exceed 1 - diffusion; see Details.
Parameters do not have default values (except negative_proportion,
diffusion, target_condition, max_correlation, and
max_iterations) and must each be set. See Details and Examples to
get started.
simulate_sbm(
nodes,
blocks,
density_matrix,
snr = 1,
diffusion = 0.3,
diffusion_range = NULL,
negative_proportion,
sample_size,
skew = 0,
skew_range = NULL,
target_condition = 30,
max_correlation = 0.8,
max_iterations = 100
)A named list with four elements:
Numeric matrix of dimension sample_size x sum(nodes)
containing the simulated observations drawn from the population GGM.
Rows are cases; columns are variables named V01, V02,
etc. Values are continuous (or skewed continuous when skew != 0).
To produce ordinal data, pass the columns through
categorize.
A list of input, derived, and estimated parameters:
nodes --- Integer vector of length blocks giving the
number of nodes per block (scalar input is expanded to this length)
blocks --- Number of community blocks
sample_size --- Number of simulated observations
skew --- Named numeric vector of per-variable skew values
actually applied (after rounding and possible resampling)
density_matrix --- The blocks x blocks density matrix
as supplied
negative_proportion --- The proportion of between-community
edges assigned a negative sign, either as supplied or as sampled from the
empirical distribution
weibull --- Named numeric vector of length 2 giving the
Weibull shape and scale parameters of the absolute edge
weight distribution actually used. If ridge conditioning was applied,
these are re-estimated from the conditioned network via MLE.
diffusion --- Numeric vector of length 2 giving the
within-block reservation range as a proportion of within-community edges,
on the internal scale used by the sampler (i.e., range(1 -
diffusion) or range(1 - diffusion_range)). Both values are equal
when diffusion is scalar. Note this is the complement of the
user-supplied diffusion value and represents the fraction of
within-community slots filled from the top-ranked draws.
omega --- Smallworldness omega statistic of the generated
network (Telesford et al., 2011). Values near zero indicate small-world
structure; negative values indicate lattice-like structure; positive
values indicate random-like structure (see smallworldness)
Q --- Newman-Girvan modularity of the population network
(Omega) with respect to the block membership, computed via
igraph::modularity on absolute edge weights. Provides a summary
of the community contrast actually achieved after weight assignment and
any ridge conditioning.
Population-level network parameters:
R --- Population correlation matrix derived from the GGM
via pcor2cor
Omega --- Population partial correlation matrix (the GGM
edge weight matrix), with zeros for absent edges
membership --- Named integer vector of length
sum(nodes) giving the community block assignment (1 to
blocks) for each node
Iteration and conditioning diagnostics:
iterations --- Number of sampling attempts needed to find a
valid network (including graph structure and edge weight draws)
rejections --- Character vector of length
max_iterations + 1 recording the rejection reason for each failed
attempt; entries for successful or unused iterations are empty strings.
Common reasons include disconnected graph structure, condition number
exceeding target_condition, maximum correlation exceeding
max_correlation, and Weibull parameters falling outside empirical
bounds after ridge conditioning.
lambda --- Ridge regularization parameter \(\lambda\)
added to the diagonal of the precision matrix to ensure positive
definiteness; NA if no conditioning was required
condition --- Condition number of the final population
correlation matrix R, computed via kappa with
exact = TRUE
Numeric (length = 1 or blocks).
Number of nodes per community block. Can be a single value applied to all
blocks, or a vector of length blocks specifying each block's size
individually. Minimum of three nodes per block. The total number of nodes
(sum(nodes)) should be between 8 and 54 to remain within the range
of the empirical networks used to fit the Weibull parameter model; values
outside this range are accepted but will trigger extrapolation and may
produce a warning from weibull_parameters.
Numeric (length = 1). Number of community blocks.
Matrix (dim = blocks x blocks).
A symmetric numeric matrix specifying edge probabilities within and between
community blocks. Diagonal entry \([i,i]\) gives the within-block edge
probability for community \(i\); off-diagonal entry \([i,j]\)
(\(i \neq j\)) gives the between-block edge probability for the pair of
communities \(i\) and \(j\). All entries must be in \([0, 1]\) and
the matrix must be symmetric (i.e., density_matrix[i,j] ==
density_matrix[j,i]). See Details for construction guidance.
Numeric (length = 1).
Signal-to-noise ratio of the absolute partial correlation weights,
defined as \(\bar{|w|} / \mathrm{SD}(|w|)\). This ratio governs the
shape and scale of the Weibull distribution used to generate edge weights:
values below 1 produce wider, more heterogeneous weight distributions;
values above 1 produce narrower, more homogeneous distributions. Note that
the same SNR can arise from different combinations of mean and standard
deviation (e.g., mean = 0.05, SD = 0.05 and mean = 0.15, SD = 0.15 both
give SNR = 1), so SNR alone does not determine the absolute magnitude of
edge weights, which is additionally governed by nodes and
sample_size via weibull_parameters. Empirically
observed SNR values ranged from 0.648 to 1.712; values outside this range
are accepted but will trigger a warning. Defaults to 1.
Numeric (length = 1).
Controls the minimum proportion of the top-ranked edges (by absolute
weight) that are guaranteed to be placed within communities rather than
distributed freely across the network. Specifically, 1 - diffusion
of the within-community edges are reserved for the highest-weight draws;
the remaining edges (both within- and between-community) are filled from
a randomly shuffled pool. Because shuffled edges can land within communities
by chance, the actual proportion of top edges that end up within communities
will on average exceed 1 - diffusion. Lower values of
diffusion (e.g., 0.10) therefore produce stronger community
contrast than higher values (e.g., 0.90), but diffusion
should be interpreted as a floor on within-community weight concentration,
not an exact control. Defaults to 0.30. Must be between 0 and 1.
Numeric (length = 2).
If provided, overrides diffusion by drawing the diffusion proportion
uniformly from this interval on each call. Useful for introducing
replication-to-replication variability in community contrast. For example,
diffusion_range = c(0.05, 0.20) samples a value between 5% and
20% on each draw. The same floor interpretation applies as for
diffusion. Both values must be between 0 and 1.
Numeric (length = 1). Proportion of between-community edges assigned a negative sign (inhibitory connections). Each between-community edge is independently signed negative with this probability (i.e., a Bernoulli draw per edge), so the realized proportion will vary around the specified value. Must be between 0 and 1. Applies only to between-community edges; all within-community edges are positive. If not provided, a value is sampled from a truncated normal distribution reflecting the empirical distribution of true negative partial correlations across 194 psychometric networks: mean = 0.34, SD = 0.086, bounded to \([0.083, 0.55]\).
Numeric (length = 1).
Number of observations to generate from the population multivariate
distribution. Also influences the predicted Weibull scale parameter via
weibull_parameters: larger samples are associated with
smaller, more precisely estimated edge weights.
Numeric (length = 1 or sum(nodes)).
Skew applied to each variable after generation from the multivariate
normal. Can be a single value applied to all variables, or one value per
variable. Values are rounded to the nearest 0.05 increment and must be
in \([-2, 2]\). Defaults to 0 (no skew).
Numeric (length = 2).
If provided, overrides skew by drawing a skew value independently
and uniformly from this interval for each variable. Both values must be
in \([-2, 2]\).
Numeric (length = 1).
Target condition number (computed via kappa with
exact = TRUE) used when ridge regularization is required to recover
a positive definite precision matrix. The smallest ridge penalty
\(\lambda\) that brings the condition number to this target is found via
root-finding (uniroot), subject to a maximum shrinkage of
approximately 23% following Peeters et al. (2020). After conditioning,
the Weibull bounds are re-verified on the updated edge weights; draws that
fall outside the empirical bounds after conditioning are rejected. Lower
values produce better-conditioned (more stable) matrices. Defaults to
30. Values up to 100 are accepted but not recommended.
Numeric (length = 1).
Maximum allowed absolute pairwise correlation in the population correlation
matrix R. Any draw where max(abs(R[lower.tri(R)])) >
max_correlation is rejected and a new attempt is made. Must be between
0 and 1. Defaults to 0.80.
Numeric (length = 1).
Maximum number of attempts to find a connected network with valid edge
weights before stopping with an error. The error message reports a
frequency table of rejection reasons to assist with diagnosing convergence
failures. Defaults to 100.
Alexander P. Christensen <alexpaulchristensen@gmail.com>
Constructing density_matrix
The density_matrix is a blocks x blocks symmetric matrix
where entry \([i,i]\) gives the within-block edge probability for
community \(i\), and entry \([i,j]\) (for \(i \neq j\)) gives the
between-block edge probability for communities \(i\) and \(j\). The
simplest construction uses a uniform off-diagonal density with
block-specific diagonals:
# Uniform within (0.90) and between (0.20) density for 3 blocks
dm <- matrix(0.20, nrow = 3, ncol = 3)
diag(dm) <- 0.90
For asymmetric community structure, each diagonal entry can differ:
# Varying within-block density per community
dm <- matrix(0.20, nrow = 3, ncol = 3)
diag(dm) <- c(0.85, 0.90, 0.95)
For full pairwise control over between-block densities, specify the complete symmetric matrix directly:
dm <- matrix(c(
0.90, 0.20, 0.10,
0.20, 0.85, 0.30,
0.10, 0.30, 0.95
), nrow = 3, ncol = 3)
Diffusion and within-community weight concentration
The diffusion parameter does not exactly fix the proportion of
top-ranked edges placed within communities. Instead, 1 - diffusion
of the within-community edge slots are filled deterministically from the
highest-weight draws. The remaining edge slots (within- and
between-community alike) are filled from a randomly shuffled pool of
lower-ranked weights, meaning some additional high-weight edges will land
within communities by chance. The realized within-community weight
advantage will therefore always be at least as large as implied by
1 - diffusion, and typically larger. The Q field in the
returned parameters list (Newman-Girvan modularity) provides a
post-hoc summary of the actual community contrast achieved.
Seminal introduction to Stochastic Block Models
Holland, P. W., Laskey, K. B., & Leinhardt, S. (1983).
Stochastic blockmodels: First steps.
Social Networks, 5(2), 109--137.
Empirical network data used to fit the Weibull SUR model
Huth, K. B. S., Haslbeck, J. M. B., Keetelaar, S., Van Holst, R. J., &
Marsman, M. (2025). Statistical evidence in psychological networks.
Nature Human Behaviour.
Maximum ridge shrinkage bound
Peeters, C. F., van de Wiel, M. A., & van Wieringen, W. N. (2020).
The spectral condition number plot for regularization parameter evaluation.
Computational Statistics, 35(2), 629--646.
# Construct density matrix for 3 blocks with uniform densities
dm <- matrix(0.20, nrow = 3, ncol = 3)
diag(dm) <- 0.90
# Basic 3-block simulation with equal-sized communities
result <- simulate_sbm(
nodes = 6, # 6 nodes per block = 18 total
blocks = 3,
sample_size = 500,
density_matrix = dm
)
# Unequal block sizes
result <- simulate_sbm(
nodes = c(4, 6, 8), # 18 total nodes
blocks = 3,
sample_size = 500,
density_matrix = dm
)
# Varying within-block density per community
dm_varying <- matrix(0.20, nrow = 3, ncol = 3)
diag(dm_varying) <- c(0.85, 0.90, 0.95)
result <- simulate_sbm(
nodes = 6,
blocks = 3,
sample_size = 500,
density_matrix = dm_varying
)
# Full pairwise between-block density control
dm_pairwise <- matrix(c(
0.90, 0.20, 0.10,
0.20, 0.85, 0.30,
0.10, 0.30, 0.95
), nrow = 3, ncol = 3)
result <- simulate_sbm(
nodes = 6,
blocks = 3,
sample_size = 500,
density_matrix = dm_pairwise
)
# Fix the proportion of negative between-community edges
result <- simulate_sbm(
nodes = 6,
blocks = 3,
sample_size = 500,
density_matrix = dm,
negative_proportion = 0.20
)
# Introduce variability in diffusion across replications
result <- simulate_sbm(
nodes = 6,
blocks = 3,
sample_size = 500,
density_matrix = dm,
diffusion_range = c(0.30, 0.70)
)
Run the code above in your browser using DataLab