mlGGM: Multi-Level Gaussian Graphical Model Estimation

Description

Estimates within-cluster and between-cluster partial correlation networks from cross-sectional multilevel data (e.g., individuals nested in classrooms, patients nested in clinics) using a single-step nodewise regression approach. Unlike mlVAR, which uses a two-step procedure (temporal model first, then contemporaneous), mlGGM fits all effects in a single multilevel model per variable. This makes it suitable for cross-sectional data where no temporal ordering exists.

For each variable $y_i$, a single mixed-effects model is fitted with within-cluster centered predictors and cluster-mean predictors. Two undirected networks (GGMs) are extracted: one capturing partial correlations at the within-cluster level (individual deviations from cluster means) and one at the between-cluster level (cluster means).

NOTE: This function is currently experimental and may change in future versions.

Usage

mlGGM(data, vars, idvar = "id",
      estimator = c("lmer"),
      randomeffects = c("default", "correlated", "orthogonal", "fixed"),
      scale = TRUE, na.rm = TRUE, verbose = TRUE,
      full_detrend = FALSE)
# S3 method for mlGGM
print(x, ...)
# S3 method for mlGGM
summary(object, show = c("fit", "within", "between"),
        round = 3, ...)
# S3 method for mlGGM
plot(x, type = c("within", "between"),
     partial = TRUE, SD = FALSE, subject, order,
     nonsig = c("default", "show", "hide", "dashed"),
     rule = c("or", "and"), alpha = 0.05,
     layout = "spring", verbose = TRUE, ...)

Value

An object of class "mlGGM" containing:

results

A list with:

within: Within-cluster network estimates (covariance, correlation, partial correlation matrices, each with $mean and $subject components)
between: Between-cluster network estimates (covariance, correlation, partial correlation matrices, with $mean component)
Gamma_within: Within-cluster regression coefficients with $mean, $SE, $P, and $subject components
Gamma_between: Between-cluster regression coefficients with $mean, $SE, and $P components
mu: Intercept estimates with $mean, $SE, and $SD components

output

A list of fitted lmer model objects, one per variable.

fit

A data frame with AIC and BIC for each variable's model.

data

The augmented data frame used in estimation.

model

The predictor model specification data frame.

input

A list of input arguments: vars, estimator, randomeffects, idvar.

Arguments

data: A data frame containing the observed variables and a cluster identifier variable.
vars: A character vector of variable names to include in the GGM estimation. Must contain at least 2 variables.
idvar: A string indicating the name of the cluster/group identifier variable in data. Defaults to "id".
estimator: The estimator to be used. Currently only "lmer" (sequential univariate multilevel estimation via lmer) is supported.
randomeffects: How should random effects be estimated? "correlated" estimates correlated random slopes for within-cluster predictors (using (1 + predictors | cluster)), "orthogonal" estimates uncorrelated random slopes (using (1 + predictors || cluster)), and "fixed" estimates only a random intercept with no random slopes. "default" selects "correlated" when the number of variables is 6 or fewer and "orthogonal" otherwise, as correlated random effects become computationally expensive with many variables.
scale: Logical. Should variables be grand-mean standardized before estimation? Defaults to TRUE. Standardization ensures that edge weights are comparable across variables.
na.rm: Logical. Should rows with missing values be removed? Defaults to TRUE.
verbose: Logical indicating if console messages and progress bars should be shown.
full_detrend: Logical. If TRUE, standardize each variable per observation position across all clusters before estimation. For example, all first observations within each cluster are standardized together, all second observations together, etc. Requires a balanced design (equal number of observations per cluster). Applied before grand-mean standardization (scale). Defaults to FALSE.
x: An mlGGM object.
object: An mlGGM object.
show: Character vector indicating which sections to show in the summary. Options: "fit" for information criteria (AIC/BIC), "within" for within-cluster partial correlations and p-values, "between" for between-cluster partial correlations and p-values.
round: Number of decimal places for rounding in the summary output. Defaults to 3.
type: The type of network to plot or extract: "within" for the within-cluster partial correlation network, "between" for the between-cluster partial correlation network.
partial: Logical. If TRUE (default), plots partial correlations. If FALSE, plots zero-order correlations.
SD: Logical. If TRUE, plots random effect standard deviations instead of fixed effects. Only available for the within-cluster network.
subject: Integer specifying which cluster's network to plot. Only available for the within-cluster network. When specified, the cluster-specific network (fixed effects + random effects for that cluster) is plotted.
order: Character vector or numeric vector specifying the order of nodes in the plot.
nonsig: How to handle non-significant edges in the plot: "show" shows all edges regardless of significance, "hide" removes non-significant edges, "dashed" shows non-significant edges as dashed lines. "default" hides non-significant edges for partial correlation networks and shows them otherwise.
rule: Significance rule for the undirected network: "or" requires at least one regression direction to be significant, "and" requires both directions to be significant. The "and" rule is more conservative and tends to improve specificity, especially for the between-cluster network.
alpha: Significance level for thresholding edges. Defaults to 0.05.
layout: Layout algorithm for the network plot, passed to qgraph. Defaults to "spring".
...: Additional arguments passed to qgraph (for plot) or ignored.

Author

Sacha Epskamp

Details

mlGGM estimates multi-level Gaussian graphical models by fitting one mixed-effects regression model per variable. For each variable $y_i$, the model is:

$$y_i = \beta_0 + \sum_{j \neq i} \gamma^{(w)}_{ij} (y_j - \bar{y}_j) + \sum_{j \neq i} \gamma^{(b)}_{ij} \bar{y}_j + u_{0} + \sum_{j \neq i} u_j (y_j - \bar{y}_j) + \varepsilon$$

where $(y_j - \bar{y}_j)$ are within-cluster centered predictors and $\bar{y}_j$ are cluster-mean predictors. The random effects $u_0, u_j$ are placed on the intercept and within-cluster predictors only (not on between-cluster means).

Two GGMs are extracted using the relationship between nodewise regression coefficients and the precision matrix:

Within-cluster network: The precision matrix is computed as $K_w = D_w (I - \Gamma_w)$, where $D_w = \mathrm{diag}(1/\sigma^2_i)$ contains the inverse residual variances and $\Gamma_w$ contains the fixed effects of the within-cluster centered predictors. The matrix is symmetrized and forced to be positive definite. Partial correlations are derived from the precision matrix.

Between-cluster network: The precision matrix is computed as $K_b = D_b (I - \Gamma_b)$, where $D_b = \mathrm{diag}(1/\mu^2_{SD,i})$ contains the inverse squared random intercept standard deviations and $\Gamma_b$ contains the fixed effects of the cluster-mean predictors.

Cluster-specific networks: When randomeffects is "correlated" or "orthogonal", cluster-specific within-cluster networks are also computed by adding the estimated random effects to the fixed effects for each cluster. These can be accessed via the subject argument in the plot method.

Network matrices can be extracted using getNet with type = "within" or type = "between". The rule argument controls the significance thresholding rule.

References

Epskamp, S., Waldorp, L. J., Mottus, R., & Borsboom, D. (2018). The Gaussian graphical model in cross-sectional and time-series data. Multivariate Behavioral Research, 53(4), 453-480.

Epskamp, S., Rhemtulla, M. T., & Borsboom, D. (2017). Generalized network psychometrics: Combining network and latent variable models. Psychometrika, 82(4), 904-927.

Examples

Run this code

if (FALSE) {
# Generate simulated multilevel data with known GGM structure:
library(bootnet)

nNode <- 4
nCluster <- 30
clusterSize <- 50

# True within-cluster and between-cluster networks:
set.seed(1)
net_within <- genGGM(nNode, propPositive = 0.8)
net_between <- genGGM(nNode, p = 1, propPositive = 0.5)

# Generate within-cluster deviations:
within_data <- ggmGenerator()(nCluster * clusterSize, net_within)
within_data <- as.data.frame(within_data)
vars <- names(within_data)
within_data$cluster <- rep(1:nCluster, each = clusterSize)

# Generate between-cluster means:
between_data <- ggmGenerator()(nCluster, net_between)
between_data <- as.data.frame(between_data)
between_data$cluster <- 1:nCluster

# Combine: observed = within deviation + cluster mean
merged <- dplyr::left_join(within_data, between_data, by = "cluster",
                           suffix = c("_w", "_b"))
data <- merged[, grep("_w$", names(merged))] +
        merged[, grep("_b$", names(merged))]
names(data) <- vars
data$cluster <- merged$cluster

# Fit the multi-level GGM:
fit <- mlGGM(data, vars = vars, idvar = "cluster")

# Inspect results:
summary(fit)

# Plot within-cluster and between-cluster networks:
plot(fit, type = "within")
plot(fit, type = "between")

# Extract network matrices (with AND rule for better specificity):
within_net <- getNet(fit, "within", rule = "and")
between_net <- getNet(fit, "between", rule = "and")
}

Run the code above in your browser using DataLab