Example Scripts

This example fits a saturated Gaussian Graphical Model (GGM) to the Big Five Inventory personality data. The saturated model estimates all partial correlations and provides valid confidence intervals and p-values for each edge.

# Load required packages
library(psychonetrics)
library(psych)
library(qgraph)

# Load Big Five Inventory data and remove missing values
data(bfi)
bfi_clean <- na.omit(bfi[, 1:25])  # 25 personality items

# Fit saturated GGM
mod_saturated <- ggm(bfi_clean) %>% runmodel()

# Inspect results
mod_saturated
mod_saturated %>% fit
mod_saturated %>% parameters

# Confidence interval plot for all partial correlations
CIplot(mod_saturated, "omega")

# Extract thresholded network (only significant edges)
omega <- getmatrix(mod_saturated, "omega", threshold = TRUE, alpha = 0.01)

# Visualize the network
qgraph(omega,
       layout = "spring",
       labels = colnames(bfi_clean),
       theme = "colorblind",
       cut = 0,
       title = "Big Five Personality Network (thresholded)")

This example demonstrates exploratory model search on a GGM using prune and stepup. Starting from a saturated model, we prune non-significant edges and then use stepup to add back edges that improve model fit.

# Load required packages
library(psychonetrics)
library(psych)
library(qgraph)

# Load Big Five Inventory data and remove missing values
data(bfi)
bfi_clean <- na.omit(bfi[, 1:25])  # 25 personality items

# Fit saturated GGM
mod_saturated <- ggm(bfi_clean) %>% runmodel()

# Prune non-significant edges at alpha = 0.01
mod_pruned <- mod_saturated %>% prune(alpha = 0.01, adjust = "fdr")

# Use stepup to add back significant edges
mod_stepup <- mod_pruned %>% stepup()

# Compare models
compare(saturated = mod_saturated,
        pruned = mod_pruned,
        stepup = mod_stepup)

# Extract partial correlation network (omega)
omega <- getmatrix(mod_stepup, "omega")

# Visualize the network
qgraph(omega,
       layout = "spring",
       labels = colnames(bfi_clean),
       theme = "colorblind",
       cut = 0,
       title = "Big Five Personality Network (model search)")

# View model fit
mod_stepup %>% fit

This example demonstrates how to discover a network structure on a training set and then confirm it on a held-out test set. This approach helps validate that the discovered structure generalizes to new data.

# Load required packages
library(psychonetrics)
library(psych)

# Load data and remove missing values
data(bfi)
bfi_clean <- na.omit(bfi[, 1:25])

# Split into training (first 1000) and test sets
set.seed(123)
indices <- sample(nrow(bfi_clean))
train_data <- bfi_clean[indices[1:1000], ]
test_data <- bfi_clean[indices[1001:nrow(bfi_clean)], ]

# TRAINING PHASE: Discover structure
mod_train_saturated <- ggm(train_data) %>% runmodel()
mod_train_pruned <- mod_train_saturated %>%
  prune(alpha = 0.01, adjust = "fdr")
mod_train_final <- mod_train_pruned %>% stepup()

# Extract adjacency matrix (network structure)
adj <- 1*(getmatrix(mod_train_final, "omega") != 0)

# TEST PHASE: Confirm structure on test data
mod_test_confirmatory <- ggm(test_data,
                              omega = adj) %>%
  runmodel()

# Evaluate confirmatory fit on test data
mod_test_confirmatory

# Compare with saturated model on test data
mod_test_saturated <- ggm(test_data) %>% runmodel()

compare(confirmatory = mod_test_confirmatory,
        saturated = mod_test_saturated)

# Check if confirmatory model fits adequately
fit(mod_test_confirmatory)

This example demonstrates multi-group network analysis, comparing personality networks between males and females. After pruning each group separately, the group structures generally differ, so we first take the union of the two structures with unionmodel() — this gives both groups the same set of edges (a prerequisite for meaningfully constraining them equal) — and then test whether the edge weights are equal across groups with groupequal().

# Load required packages
library(psychonetrics)
library(dplyr)
library(psych)

# Load data and prepare
data(bfi)
bfi_clean <- na.omit(bfi[, c(1:25, 26)])  # Items + gender
bfi_clean <- bfi_clean[bfi_clean$gender %in% c(1, 2), ]  # 1=male, 2=female

# Fit multi-group GGM with different networks per group
mod_configural <- ggm(bfi_clean,
                      vars = colnames(bfi_clean)[1:25],
                      groups = "gender") %>%
  runmodel()

# Prune each group separately
mod_pruned <- mod_configural %>%
  prune(alpha = 0.01)

# Take the union of the two structures, so both groups
# have the same edges (required before equality constraints):
mod_union <- mod_pruned %>%
  unionmodel() %>%
  runmodel()

# Constrain edge weights to be equal across groups
mod_equal <- mod_union %>%
  groupequal("omega") %>%
  runmodel()

# Compare models
compare(configural = mod_union,
        equal_networks = mod_equal)

# Extract networks for each group from the pruned model
omega_male   <- getmatrix(mod_pruned, "omega", group = 1)
omega_female <- getmatrix(mod_pruned, "omega", group = 2)

# Compare network densities
mean(omega_male[upper.tri(omega_male)] != 0)
mean(omega_female[upper.tri(omega_female)] != 0)

The chi-square difference test rejects exact equality of the edge weights (p < 0.0001), although BIC slightly favors the more parsimonious equal model. Use partialprune() to identify which edges differ between groups.

Note: Penalized estimation is experimental and has not yet been fully validated.

This example demonstrates penalized maximum likelihood estimation for sparse network estimation. The PML estimator automatically searches for an optimal regularization parameter and produces a sparse network.

# Load required packages
library(psychonetrics)
library(psych)
library(qgraph)

# Load data
data(bfi)
bfi_clean <- na.omit(bfi[, 1:25])

# Fit with penalized maximum likelihood
# Automatic lambda selection via EBIC
mod_pml <- ggm(bfi_clean,
               estimator = "PML") %>%
  runmodel()

# View the lambda search results
mod_pml@optim$lambda_search

# Extract the sparse network
omega_sparse <- getmatrix(mod_pml, "omega")

# Count number of edges
n_edges <- sum(omega_sparse[upper.tri(omega_sparse)] != 0)
cat("Number of edges:", n_edges, "\n")

# Refit without penalty for inference
mod_refit <- mod_pml %>% refit

# View refitted model with standard errors
mod_refit

# Visualize sparse network
qgraph(omega_sparse,
       layout = "spring",
       labels = colnames(bfi_clean),
       theme = "colorblind",
       cut = 0,
       title = "Penalized GGM (PML)")

# Compare sparsity
omega_ml <- ggm(bfi_clean) %>% runmodel() %>% getmatrix("omega")
cat("ML edges:", sum(omega_ml[upper.tri(omega_ml)] != 0), "\n")
cat("PML edges:", sum(omega_sparse[upper.tri(omega_sparse)] != 0), "\n")

This example demonstrates how to handle missing data using Full Information Maximum Likelihood (FIML). This approach uses all available information and is generally superior to listwise deletion.

# Load required packages
library(psychonetrics)
library(psych)

# Load data WITH missing values (do not remove NAs)
data(bfi)
bfi_missing <- bfi[, 1:25]

# Check missing data pattern
cat("Proportion missing:", mean(is.na(bfi_missing)), "\n")
cat("Complete cases:", sum(complete.cases(bfi_missing)), "/", nrow(bfi_missing), "\n")

# Fit with FIML (automatically handles missing data)
mod_fiml <- ggm(bfi_missing,
                missing = "auto") %>%  # Automatically switches to FIML
  runmodel()

# Prune and stepup
mod_fiml_pruned <- mod_fiml %>%
  prune(alpha = 0.01, adjust = "fdr") %>%
  stepup()

# Compare with listwise deletion approach
bfi_complete <- na.omit(bfi_missing)
mod_listwise <- ggm(bfi_complete) %>%
  runmodel() %>%
  prune(alpha = 0.01, adjust = "fdr") %>%
  stepup()

# Extract networks
omega_fiml <- getmatrix(mod_fiml_pruned, "omega")
omega_listwise <- getmatrix(mod_listwise, "omega")

# Compare network densities
cat("FIML edges:", sum(omega_fiml[upper.tri(omega_fiml)] != 0), "\n")
cat("Listwise edges:", sum(omega_listwise[upper.tri(omega_listwise)] != 0), "\n")

# Compare sample sizes used
cat("FIML uses all n =", nrow(bfi_missing), "\n")
cat("Listwise uses n =", nrow(bfi_complete), "\n")

# View FIML model
mod_fiml_pruned

This example demonstrates how to fit a GGM when you only have summary statistics (correlation matrix and sample size) rather than raw data. This is useful when working with published correlation matrices.

# Load required packages
library(psychonetrics)
library(psych)
library(qgraph)

# Load Thurstone correlation matrix (built into psych package)
data(Thurstone)

# The Thurstone object contains correlation matrices
# Use the first correlation matrix
cor_matrix <- Thurstone

# Sample size (reported in Thurstone dataset)
n <- 213

# Fit GGM from correlation matrix
mod_summary <- ggm(covs = cor_matrix,
                   nobs = n) %>%
  runmodel()

# Prune and stepup
mod_summary_pruned <- mod_summary %>%
  prune(alpha = 0.01, adjust = "fdr") %>%
  stepup()

# View model
mod_summary_pruned

# Extract network
omega <- getmatrix(mod_summary_pruned, "omega")

# Visualize
qgraph(omega,
       layout = "spring",
       labels = rownames(cor_matrix),
       theme = "colorblind",
       cut = 0,
       title = "Thurstone Abilities Network")

# Note: Some functionality requiring raw data will not be available
# (e.g., bootstrapping, FIML with missing data)
# But model estimation, pruning, and comparison work fine

# Can also do multi-group from multiple correlation matrices
# Example structure (not run):
# mod_mg <- ggm(covs = list(cor1, cor2),
#               nobs = c(n1, n2),
#               groups = c(rep(1, n1), rep(2, n2)))

This example demonstrates a classic confirmatory factor analysis with three correlated factors. We use the Holzinger and Swineford dataset with visual, textual, and speed factors.

# Load required packages
library(psychonetrics)
library(lavaan)

# Load classic CFA dataset
data(HolzingerSwineford1939)
hs_data <- HolzingerSwineford1939[, c(paste0("x", 1:9))]

# Define lambda (factor loading) matrix
# 3 factors: visual (x1-x3), textual (x4-x6), speed (x7-x9)
lambda <- matrix(0, 9, 3)
lambda[1:3, 1] <- 1  # Visual factor
lambda[4:6, 2] <- 1  # Textual factor
lambda[7:9, 3] <- 1  # Speed factor

# Add labels for factor names
colnames(lambda) <- c("visual", "textual", "speed")

# Fit CFA model
mod_cfa <- lvm(hs_data,
               lambda = lambda,
               identification = "variance") %>%  # Fix latent variances to 1
  runmodel()

# View results
mod_cfa

# Check fit indices
fit(mod_cfa)

# Extract factor loadings
params <- parameters(mod_cfa)
loadings <- params[params$matrix == "lambda", ]
print(loadings)

# Extract latent covariance matrix
sigma_zeta <- getmatrix(mod_cfa, "sigma_zeta")
print(sigma_zeta)

# Modification indices
mod_cfa %>% MIs

# R-squared for each indicator
residual_var <- diag(getmatrix(mod_cfa, "sigma_epsilon"))
total_var <- diag(var(hs_data))
r_squared <- 1 - (residual_var / total_var)
names(r_squared) <- colnames(hs_data)
print(r_squared)

This example demonstrates a Latent Network Model where the relationships between latent factors are modeled as a network (GGM) rather than a covariance matrix. This allows testing which latent factors are conditionally independent. Note that lnm() already sets latent = "ggm", so this argument should not be passed again.

# Load required packages
library(psychonetrics)
library(dplyr)
library(qgraph)

# Load Star Wars dataset (10 items, 3 trilogy factors)
data(StarWars)

lambda <- matrix(0, 10, 3)
lambda[1:4, 1]       <- 1  # Prequels:  Q1-Q4
lambda[c(1, 5:7), 2] <- 1  # Original:  Q1, Q5-Q7
lambda[c(1, 8:10), 3] <- 1 # Sequels:   Q1, Q8-Q10

obsvars <- paste0("Q", 1:10)
latents <- c("Prequels", "Original", "Sequels")

# Fit Latent Network Model (lnm() sets latent = "ggm" internally)
mod_lnm <- lnm(StarWars,
               lambda = lambda,
               vars = obsvars,
               latents = latents,
               identification = "variance") %>%
  runmodel()

# Prune non-significant edges in the latent network
mod_lnm_pruned <- mod_lnm %>%
  prune(alpha = 0.05, matrices = "omega_zeta")

# Compare models
compare(saturated = mod_lnm,
        pruned = mod_lnm_pruned)

# Extract and visualize the latent network
omega_latent <- getmatrix(mod_lnm_pruned, "omega_zeta")
qgraph(omega_latent,
       layout = "spring",
       labels = latents,
       theme = "colorblind",
       cut = 0, vsize = 14,
       title = "Star Wars Latent Factor Network")

The edge between Prequels and Original ratings can be removed without significant loss of fit (p = 0.58): liking the prequels and liking the original trilogy are conditionally independent given liking of the sequels.

This example demonstrates a Residual Network Model where residual associations between indicators are modeled as a network of partial correlations. This identifies direct associations between items beyond the latent factors. Note that rnm() already sets residual = "ggm"; with ML estimation a saturated residual network is not identified, so the standard workflow starts from an empty residual network (omega_epsilon = "zero", the default) and adds edges with stepup().

# Load required packages
library(psychonetrics)
library(dplyr)
library(qgraph)

data(StarWars)

# Same factor structure as the LNM example
lambda <- matrix(0, 10, 3)
lambda[1:4, 1]       <- 1
lambda[c(1, 5:7), 2] <- 1
lambda[c(1, 8:10), 3] <- 1

obsvars <- paste0("Q", 1:10)
latents <- c("Prequels", "Original", "Sequels")

# Reference CFA (no residual network):
mod_cfa <- lvm(StarWars, lambda = lambda, vars = obsvars,
               latents = latents, identification = "variance") %>%
  runmodel()

# RNM: start with an empty residual network (default)
mod_rnm <- rnm(StarWars, lambda = lambda, vars = obsvars,
               latents = latents, identification = "variance",
               omega_epsilon = "zero") %>%
  runmodel()

# Add significant residual edges via stepup
mod_rnm_stepup <- mod_rnm %>%
  stepup(matrices = "omega_epsilon", alpha = 0.01)

# Compare with the standard CFA
compare(cfa = mod_cfa,
        rnm_stepup = mod_rnm_stepup)

# Extract and visualize the residual network
omega_residual <- getmatrix(mod_rnm_stepup, "omega_epsilon")
qgraph(omega_residual,
       layout = "spring",
       labels = obsvars,
       theme = "colorblind",
       cut = 0,
       title = "Star Wars Residual Network")

This example reproduces the classic Bollen (1989) industrialization-and-democracy structural equation model: industrialization in 1960 (ind60) predicts democracy in 1960 (dem60) and 1965 (dem65). The textbook model includes six correlated residuals between democracy indicators measured with the same method or at adjacent time points. Remember that in psychonetrics specification matrices 0 means fixed to zero and 1 means free — so the diagonal of sigma_epsilon must be 1 (free residual variances).

# Load required packages
library(psychonetrics)
library(dplyr)
library(lavaan)

# Load political democracy dataset
data(PoliticalDemocracy)

# y1-y4: democracy indicators 1960; y5-y8: democracy indicators 1965
# x1-x3: industrialization indicators 1960
vars <- c(paste0("y", 1:8), paste0("x", 1:3))
data_sem <- PoliticalDemocracy[, vars]

# Measurement model (lambda): 3 latent factors
lambda <- matrix(0, 11, 3)
lambda[1:4, 2]  <- 1     # dem60 measured by y1-y4
lambda[5:8, 3]  <- 1     # dem65 measured by y5-y8
lambda[9:11, 1] <- 1     # ind60 measured by x1-x3
colnames(lambda) <- c("ind60", "dem60", "dem65")
rownames(lambda) <- vars

# Structural model (beta): rows = outcomes, columns = predictors
beta <- matrix(0, 3, 3)
beta[2, 1] <- 1  # dem60 ~ ind60
beta[3, 1] <- 1  # dem65 ~ ind60
beta[3, 2] <- 1  # dem65 ~ dem60
rownames(beta) <- colnames(beta) <- c("ind60", "dem60", "dem65")

# Residual covariance structure:
# free residual variances (diag = 1) plus the six classic
# correlated residuals: y1~~y5, y2~~y4, y2~~y6, y3~~y7, y4~~y8, y6~~y8
sigma_epsilon <- diag(11)
resid_pairs <- rbind(c(1,5), c(2,4), c(2,6), c(3,7), c(4,8), c(6,8))
for (i in seq_len(nrow(resid_pairs))) {
  sigma_epsilon[resid_pairs[i,1], resid_pairs[i,2]] <- 1
  sigma_epsilon[resid_pairs[i,2], resid_pairs[i,1]] <- 1
}

# Fit SEM (default identification: first loading fixed to 1)
mod_sem <- lvm(data_sem,
               lambda = lambda,
               beta = beta,
               sigma_epsilon = sigma_epsilon,
               vars = vars) %>%
  runmodel

# View results
mod_sem
fit(mod_sem)

# Structural coefficients
params <- parameters(mod_sem)
params[params$matrix == "beta" & !params$fixed,
       c("var1", "var2", "est", "se", "p")]

Industrialization predicts democracy at both time points, and democracy is highly stable over time. The model fits well: χ²(35) = 38.13, p = 0.33.

This example tests measurement invariance of the classic three-factor Holzinger–Swineford model across the two schools in the data. Note that with raw-data input, groups must be the name of a column in the data (not a separate vector). We progressively constrain parameters to test configural, weak, strong, and strict invariance.

# Load required packages
library(psychonetrics)
library(dplyr)
library(lavaan)

data(HolzingerSwineford1939)

# Keep the grouping column IN the data frame:
hs_data <- HolzingerSwineford1939[, c("school", paste0("x", 1:9))]

# Define factor structure
lambda <- matrix(0, 9, 3)
lambda[1:3, 1] <- 1   # visual
lambda[4:6, 2] <- 1   # textual
lambda[7:9, 3] <- 1   # speed
colnames(lambda) <- c("visual", "textual", "speed")

# CONFIGURAL INVARIANCE (same structure, all parameters free)
mod_configural <- lvm(hs_data,
                      lambda = lambda,
                      vars = paste0("x", 1:9),
                      groups = "school",     # column name in hs_data
                      identification = "variance") %>%
  runmodel()

# WEAK INVARIANCE (equal factor loadings)
mod_weak <- mod_configural %>%
  groupequal("lambda") %>%
  runmodel()

# STRONG INVARIANCE (equal loadings + equal intercepts)
mod_strong <- mod_weak %>%
  groupequal("nu") %>%
  runmodel()

# STRICT INVARIANCE (+ equal residual variances)
mod_strict <- mod_strong %>%
  groupequal("sigma_epsilon") %>%
  runmodel()

# Compare all models
compare(configural = mod_configural,
        weak = mod_weak,
        strong = mod_strong,
        strict = mod_strict)

# Latent means in group 2 (group 1 is fixed to zero for identification;
# meaningful under at least strong invariance). The latent mean
# matrix in lvm models is called "nu_eta":
params <- parameters(mod_strong)
params[params$matrix == "nu_eta" & params$group_id == 2,
       c("var1", "est", "se", "p")]

Weak (metric) invariance holds (p = 0.22), but strong (scalar) invariance is clearly rejected (p < 0.0001) — some intercepts differ between the schools. Use MIs(matrices = "nu", type = "free") to locate the non-invariant intercepts and free them with groupfree() for partial invariance (see the Star Wars invariance example).

This example demonstrates how to obtain bootstrap confidence intervals for model parameters. The loop_psychonetrics() function handles parallel cluster setup and variable export automatically. We use a small number of bootstrap samples for demonstration purposes, but real analyses should use 1000 or more bootstraps.

# Load required packages
library(psychonetrics)
library(dplyr)
library(lavaan)

# Load data
data(HolzingerSwineford1939)
hs_data <- HolzingerSwineford1939[, c(paste0("x", 1:9))]

# Define CFA model
lambda <- matrix(0, 9, 3)
lambda[1:3, 1] <- 1
lambda[4:6, 2] <- 1
lambda[7:9, 3] <- 1
colnames(lambda) <- c("visual", "textual", "speed")

# First, fit the original model:
mod <- lvm(hs_data,
           lambda = lambda,
           identification = "variance") %>%
  runmodel()

# Run 100 bootstrap replicates in parallel (use 1000+ in real analysis):
bootstraps <- loop_psychonetrics({
  lvm(hs_data, lambda = lambda,
      identification = "variance",
      bootstrap = "nonparametric") %>%
    runmodel()
}, reps = 100, nCores = 4)

# Aggregate bootstrap results:
boot_agg <- aggregate_bootstraps(
  sample = mod,
  bootstraps = bootstraps
)

# View bootstrap parameters (includes se_boot, CIs, prop_non0):
boot_agg %>% parameters

# Plot bootstrap confidence intervals:
CIplot(boot_agg)

This example demonstrates fitting an Ising model to binary data. The Ising model is specifically designed for binary variables and models the log-odds of conditional relationships between variables.

# Load required packages
library(psychonetrics)
library(qgraph)

# Load Jonas dataset (binary symptom data)
data(Jonas)

# Inspect data structure
head(Jonas)
dim(Jonas)  # 10 binary items

# Select only the binary items (columns 1-10)
binary_items <- Jonas[, 1:10]

# Check that data is binary
apply(binary_items, 2, table)

# Fit saturated Ising model
mod_ising_saturated <- Ising(binary_items) %>%
  runmodel()

# Prune non-significant edges
mod_ising_pruned <- mod_ising_saturated %>%
  prune(alpha = 0.01, adjust = "fdr")

# Use stepup to refine
mod_ising_stepup <- mod_ising_pruned %>%
  stepup()

# Compare models
compare(saturated = mod_ising_saturated,
        pruned = mod_ising_pruned,
        stepup = mod_ising_stepup)

# Extract network (omega matrix contains edge weights)
omega_ising <- getmatrix(mod_ising_stepup, "omega")

# Extract thresholds
thresholds <- getmatrix(mod_ising_stepup, "tau")
print(thresholds)

# Visualize Ising network
qgraph(omega_ising,
       layout = "spring",
       labels = colnames(binary_items),
       theme = "colorblind",
       cut = 0,
       title = "Ising Network",
       edge.labels = FALSE)

# View model
mod_ising_stepup

# Centrality analysis
library(qgraph)
cent <- centrality_auto(omega_ising)
print(cent$node.centrality)

# Note: In Ising models, edges represent log-odds ratios
# Positive edge: both variables tend to be in same state
# Negative edge: variables tend to be in opposite states

This example demonstrates multi-group Ising analysis with progressive equality constraints, comparing attitude networks between people who know researcher Jonas personally and those who do not. With raw-data input, groups must be a column name in the data. The equal-network model converges to a local optimum from the default start values, so we use emergencystart() to compute fresh starting values for that step.

# Load required packages
library(psychonetrics)
library(dplyr)

# Load Jonas dataset (10 binary items + a "group" column)
data(Jonas)
ising_vars <- colnames(Jonas)[1:10]

# DENSE MODELS (start with saturated)

# 1. Configural: different networks per group
mod_configural <- Ising(Jonas,
                        vars = ising_vars,
                        groups = "group") %>%   # column name in Jonas
  runmodel()

# 2. Equal networks (omega). The default start values lead to a
#    local optimum here; emergencystart() computes new start values.
mod_equal_omega <- mod_configural %>%
  groupequal("omega") %>%
  emergencystart() %>%
  runmodel()

# 3. Equal networks + equal thresholds
mod_equal_omega_tau <- mod_equal_omega %>%
  groupequal("tau") %>%
  runmodel()

# 4. Fully equal (also the inverse temperature beta)
mod_fully_equal <- mod_equal_omega_tau %>%
  groupequal("beta") %>%
  runmodel()

# Compare the dense models:
compare(configural = mod_configural,
        equal_omega = mod_equal_omega,
        equal_omega_tau = mod_equal_omega_tau,
        fully_equal = mod_fully_equal)

# SPARSE MODELS (prune first, then constrain)
mod_sparse <- mod_configural %>%
  prune(alpha = 0.01)

# Union of group structures, then equality:
mod_sparse_equal <- mod_sparse %>%
  unionmodel() %>%
  groupequal("omega") %>%
  runmodel()

compare(sparse_configural = mod_sparse,
        sparse_equal = mod_sparse_equal)

# Compare thresholds across groups (from the configural model):
tau_group1 <- getmatrix(mod_configural, "tau", group = 1)
tau_group2 <- getmatrix(mod_configural, "tau", group = 2)
data.frame(item = ising_vars,
           group1 = round(tau_group1, 2),
           group2 = round(tau_group2, 2))

The network structure (omega) can be constrained equal across groups without significant loss of fit (p = 0.45), but the thresholds clearly differ (p = 0.0029): people who know Jonas endorse the items at different base rates, while the conditional association structure is the same.

This example uses open survey data on sleep and wellbeing to estimate a GGM with Full Information Maximum Likelihood (FIML) for handling missing data. We demonstrate network pruning, exploratory stepup search from an empty model, and both bootnet-based and psychonetrics-native bootstrapping for edge stability.

# Load packages:
library("psychonetrics")
library("dplyr")
library("qgraph")

# Load data (included in psychonetrics):
data(NA2020)

# Variables (8 items on sleep and wellbeing):
#   regular_sleep     - I try to keep a regular sleep pattern
#   worried_sleep     - I am worried about my current sleeping behavior
#   sleep_interfere   - My sleep interferes with my daily functioning
#   happy_health      - I am happy with my physical health
#   optimistic_future - I feel optimistic about the future
#   very_happy        - I am very happy
#   feel_alone        - I often feel alone
#   happy_love_life   - I am happy with my love life

# Fit a saturated GGM using FIML (handles missing data):
mod <- ggm(NA2020, estimator = "FIML") %>% runmodel

# Confidence interval plot for all partial correlations:
CIplot(mod, "omega")

# Extract thresholded network (alpha = 0.05):
network <- mod %>% getmatrix("omega", threshold = TRUE, alpha = 0.05)

# Custom layout and labels:
Layout <- matrix(c(
  1, 0.28, 0.77, 0.01, -0.93, -0.75, 0.04, -1,
  0.59, -1, -0.71, 1, 0.87, 0.27, -0.04, -0.61
), ncol = 2)

Labels <- c("regular\nsleep", "worried\nsleep", "sleep\ninterfere",
            "happy\nhealth", "optimistic\nfuture", "very\nhappy",
            "feel\nalone", "happy\nlove life")

# Visualize:
qgraph(network, layout = Layout, labels = Labels,
       vsize = 15, theme = "colorblind", cut = 0, label.scale.equal = TRUE)

# --- Model search strategies ---

# Strategy 1: Prune non-significant edges
mod_prune <- mod %>% prune(alpha = 0.05)

# Strategy 2: Prune then stepup
mod_stepup <- mod_prune %>% stepup

# Strategy 3: Start empty, then stepup
mod_empty_stepup <- ggm(NA2020, estimator = "FIML", omega = "zero") %>%
  runmodel %>% stepup

# Compare all strategies:
compare(
  saturated = mod,
  pruned    = mod_prune,
  stepup    = mod_stepup,
  from_empty = mod_empty_stepup
)

# --- Bootstrapping with bootnet ---
library("bootnet")

# Custom estimation function for bootnet:
f <- function(x) {
  library("psychonetrics")
  library("dplyr")
  ggm(x, estimator = "FIML") %>%
    runmodel %>%
    prune(alpha = 0.05) %>%
    getmatrix("omega")
}

# Estimate and bootstrap
# (100 bootstraps for illustration; use 1000+ in real analyses):
net <- estimateNetwork(NA2020, f = f)
boots <- bootnet(net, nBoots = 100, nCores = 4, type = "nonparametric")
plot(boots, split0 = TRUE, order = "sample", plot = "interval")

# --- Native psychonetrics bootstrapping ---
# loop_psychonetrics() handles cluster setup and variable export automatically:
# (100 replicates for illustration; use 1000+ in real analyses)
bootstraps <- loop_psychonetrics({
  ggm(NA2020, estimator = "FIML",
      bootstrap = "nonparametric") %>%
    runmodel %>%
    prune(matrices = "omega", alpha = 0.05)
}, reps = 100, nCores = 4)

# Aggregate and inspect:
mod_agg <- aggregate_bootstraps(sample = mod_prune, bootstraps = bootstraps)
mod_agg %>% parameters
mod_agg %>% fit
CIplot(mod_agg, prop0_cex = 2)

This example uses a Star Wars fandom questionnaire to fit a three-factor CFA (Prequels, Original Trilogy, Sequels), then extends it to a latent network model where latent covariances are replaced by a GGM. We also demonstrate residual network models, modification indices, and model comparison via χ² difference tests and information criteria.

# Load packages:
library("psychonetrics")
library("dplyr")
library("qgraph")

# Load data (included in psychonetrics):
data(StarWars)
Data <- StarWars

# Variables:
# Q1:  I am a huge Star Wars fan! (star what?)
# Q2:  I would trust this person with my democracy
# Q3:  I enjoyed the story of Anakin's early life
# Q4:  The special effects in this scene are awful (Battle of Geonosis)
# Q5:  I would trust this person with my life
# Q6:  I found Darth Vader's big reveal one of the greatest movie moments
# Q7:  The special effects in this scene are amazing (Death Star Explosion)
# Q8:  If possible, I would definitely buy this droid
# Q9:  The story in the sequels is an improvement to previous movies
# Q10: The special effects in this scene are marvellous (Starkiller Base)

# --- Confirmatory Factor Analysis ---

# Define factor loading structure:
Lambda <- matrix(0, 10, 3)
Lambda[1:4, 1]     <- 1   # Prequels factor
Lambda[c(1,5:7), 2] <- 1  # Original Trilogy factor
Lambda[c(1,8:10), 3] <- 1 # Sequels factor

obsvars <- paste0("Q", 1:10)
latents <- c("Prequels", "Original", "Sequels")

# Fit CFA:
mod1 <- lvm(Data, lambda = Lambda, vars = obsvars,
            identification = "variance", latents = latents)
mod1 <- mod1 %>% runmodel

# Inspect:
mod1 %>% parameters
mod1 %>% MIs

# Add residual covariance suggested by MIs:
mod2 <- mod1 %>%
  freepar("sigma_epsilon", "Q10", "Q4") %>%
  runmodel

# Compare models:
compare(original = mod1, adjusted = mod2)
mod2 %>% fit

# Stepup search on loadings and residuals:
mod3 <- mod2 %>% stepup(matrices = c("lambda", "sigma_epsilon"))
compare(stepup = mod3, adjusted = mod2)

# --- Latent Network Model ---
# Replace latent covariance matrix with a GGM:
lnm_mod <- lvm(Data, lambda = Lambda, vars = obsvars,
               latents = latents, identification = "variance",
               latent = "ggm")
lnm_mod <- lnm_mod %>% runmodel

# Prune non-significant latent edges:
lnm_mod <- lnm_mod %>% prune(alpha = 0.05)

# Compare CFA vs latent network model:
compare(cfa = mod1, lnm = lnm_mod)

# Plot latent network:
qgraph(lnm_mod@modelmatrices[[1]]$omega_zeta,
       labels = latents, theme = "colorblind", cut = 0, vsize = 10)

# --- Residual Network Model ---
# Add partial correlations at the residual level:
rnm_mod <- rnm(Data, lambda = Lambda, vars = obsvars,
               latents = latents, identification = "variance")
rnm_mod <- rnm_mod %>% stepup

# --- GGM on observed variables ---
# For comparison, fit an exploratory GGM directly:
S <- (nrow(Data) - 1) / nrow(Data) * cov(Data[, 1:10])
ggm_mod <- ggm(covs = S, nobs = nrow(Data), omega = "zero")
ggm_mod <- ggm_mod %>% stepup %>% prune

# Compare all approaches:
compare(
  cfa_original   = mod1,
  cfa_adjusted   = mod2,
  latent_network = lnm_mod,
  exploratory_ggm = ggm_mod
)

This example demonstrates multi-group measurement invariance testing using the Star Wars questionnaire. We split respondents by age (above/below 30) and test configural, weak, strong, and strict invariance, followed by latent mean comparisons. This showcases the groupequal() and groupfree() functions for constraining and freeing parameters across groups.

# Load packages:
library("psychonetrics")
library("dplyr")

# Load data (included in psychonetrics):
data(StarWars)
Data <- StarWars

# Create age grouping variable:
Data$agegroup <- ifelse(Data$Q12 < 30, "young", "less young")
table(Data$agegroup)

# Setup:
obsvars <- paste0("Q", 1:10)
latents <- c("Prequels", "Original", "Sequels")

Lambda <- matrix(0, 10, 3)
Lambda[1:4, 1]      <- 1
Lambda[c(1,5:7), 2] <- 1
Lambda[c(1,8:10), 3] <- 1

# Allow residual covariance between Q4 and Q10:
Theta <- diag(1, 10)
Theta[4, 10] <- Theta[10, 4] <- 1

# --- Configural invariance ---
mod_configural <- lvm(Data, lambda = Lambda, vars = obsvars,
                      latents = latents, sigma_epsilon = Theta,
                      identification = "variance",
                      groups = "agegroup")
mod_configural <- mod_configural %>% runmodel
mod_configural %>% fit

# --- Weak invariance (equal loadings) ---
mod_weak <- mod_configural %>% groupequal("lambda") %>% runmodel
compare(configural = mod_configural, weak = mod_weak)

# --- Strong invariance (equal intercepts) ---
mod_strong <- mod_weak %>% groupequal("nu") %>% runmodel
compare(configural = mod_configural, weak = mod_weak,
        strong = mod_strong)

# Check which intercept is problematic:
mod_strong %>% MIs(matrices = "nu", type = "free")

# Free intercept for Q10 (partial strong invariance):
mod_strong_partial <- mod_strong %>%
  groupfree("nu", 10) %>% runmodel

compare(configural = mod_configural, weak = mod_weak,
        strong = mod_strong, strong_partial = mod_strong_partial)

# --- Strict invariance (equal residual variances) ---
mod_strict <- mod_strong_partial %>%
  groupequal("sigma_epsilon") %>% runmodel

compare(configural = mod_configural, weak = mod_weak,
        strong_partial = mod_strong_partial, strict = mod_strict)

# --- Latent mean comparisons ---
# Equal latent variances:
mod_eqvar <- mod_strict %>% groupequal("sigma_zeta") %>% runmodel
compare(strict = mod_strict, eqvar = mod_eqvar)

# Equal latent means:
mod_eqmeans <- mod_eqvar %>% groupequal("nu_eta") %>% runmodel
compare(eqvar = mod_eqvar, eqmeans = mod_eqmeans)

# Rejected! Free prequel means (largest MI):
mod_eqmeans_partial <- mod_eqvar %>%
  groupequal("nu_eta", row = c("Original", "Sequels")) %>%
  runmodel
compare(eqvar = mod_eqvar, eqmeans_partial = mod_eqmeans_partial)

# Look at mean differences:
mod_eqvar %>% getmatrix("nu_eta")
# People over 30 like the prequels better!

This example fits a graphical VAR model to ESM data from a single clinical patient (Epskamp, van Borkulo, van der Veen, Servaas, Isvoranu, Riese, & Cramer, 2018). The model estimates both a temporal network (directed lag-1 effects) and a contemporaneous network (partial correlations among innovations). Data available at https://osf.io/c8wjz/.

# Load required packages
library("psychonetrics")
library("dplyr")
library("qgraph")

# Download the data from https://osf.io/c8wjz/
tsdata <- read.csv("Supplementary2_data.csv")

# Encode time variable:
tsdata$time <- as.POSIXct(tsdata$time, tz = "Europe/Amsterdam")

# Extract day variable:
tsdata$Day <- as.Date(tsdata$time, tz = "Europe/Amsterdam")

# Variables to use (7 ESM items):
vars <- c("relaxed", "sad", "nervous", "concentration",
          "tired", "rumination", "bodily.discomfort")

# Fit graphical VAR with FIML for missing data:
model <- gvar(tsdata, vars = vars, dayvar = "Day",
              estimator = "FIML")
model <- model %>% runmodel

# Inspect model:
model
model %>% fit
model %>% parameters

# Confidence interval plot for temporal effects:
CIplot(model, "beta")

# Model search:
model_pruned <- model %>% prune(alpha = 0.01)
model_final <- model_pruned %>% stepup(criterion = "bic")

# Compare models:
compare(saturated = model, pruned = model_pruned, final = model_final)

# Extract networks (use PDC, not beta!):
temporal <- getmatrix(model_final, "PDC")
contemporaneous <- getmatrix(model_final, "omega_zeta")

# Visualize both networks side by side:
labs <- gsub("\\.", "\n", vars)
L <- averageLayout(temporal, contemporaneous)

layout(t(1:2))
qgraph(temporal, layout = L, theme = "colorblind", cut = 0,
       directed = TRUE, diag = TRUE,
       title = "Temporal (PDC)", vsize = 12, mar = rep(6, 4),
       asize = 5, labels = labs)
qgraph(contemporaneous, layout = L, theme = "colorblind", cut = 0,
       title = "Contemporaneous (PCC)", vsize = 12, mar = rep(6, 4),
       asize = 5, labels = labs)

This example fits a random-intercept cross-lagged panel model (RI-CLPM) to simulated 4-wave panel data on three variables. The data-generating process combines correlated stable trait levels (random intercepts) with a stationary within-person VAR(1) process. We then use ri_clpm_search() to test which stationarity constraints the data support. See the panel-data page for background.

library(psychonetrics)
library(dplyr)

# --- Simulate 4-wave panel data (n = 500, 3 variables) ---
set.seed(42)
n <- 500; nVar <- 3; nWave <- 4

# Correlated random intercepts (stable trait levels):
Sigma_RI <- matrix(0.4, nVar, nVar); diag(Sigma_RI) <- 1
RI <- matrix(rnorm(n * nVar), n, nVar) %*% chol(Sigma_RI)

# Within-person VAR(1): autoregression 0.4, cross-lag 0.15
Beta <- matrix(c(0.4, 0.15, 0,
                 0,   0.4,  0.15,
                 0,   0,    0.4), nVar, nVar, byrow = TRUE)
Psi <- matrix(0.3, nVar, nVar); diag(Psi) <- 1

within <- matrix(rnorm(n * nVar), n, nVar) %*% chol(Psi)
Data <- matrix(NA, n, nVar * nWave)
for (w in seq_len(nWave)) {
  Data[, (w - 1) * nVar + seq_len(nVar)] <- RI + within
  if (w < nWave)
    within <- within %*% t(Beta) +
              matrix(rnorm(n * nVar), n, nVar) %*% chol(Psi)
}
colnames(Data) <- paste0("V", rep(1:nVar, nWave),
                         "_W", rep(1:nWave, each = nVar))
Data <- as.data.frame(Data)

# --- Fit the RI-CLPM ---
# Design matrix: rows = variables, columns = waves
design <- matrix(colnames(Data), nVar, nWave)

mod_riclpm <- ri_clpm(Data, vars = design) %>% runmodel

mod_riclpm %>% fit

# Wave-2-on-wave-1 temporal effects (recovers the true Beta):
round(getmatrix(mod_riclpm, "beta")[4:6, 1:3], 2)

# Random-intercept covariances (last 3 latents):
round(getmatrix(mod_riclpm, "sigma_zeta")[13:15, 13:15], 2)

# --- Stationarity search ---
search <- ri_clpm_search(mod_riclpm, criterion = "BIC")
search

The model recovers the simulated temporal effects and random-intercept correlations. Because the data were simulated as fully stationary, the search accepts every stationarity constraint and ends at the panel-VAR model (wave 1 endogenous). The selected fitted model is available as search@selected.

This example pools the correlation matrices of 12 simulated studies into a single meta-analytic Gaussian Graphical Model (MAGNA; Epskamp, Isvoranu & Cheung, 2022). The studies share a chain network (V1–V2–V3–V4, partial correlations 0.3) with mild between-study heterogeneity. We use the Olkin–Siotani sampling-covariance method added in 0.16. See the meta-analysis page for details.

library(psychonetrics)
library(dplyr)

# --- Simulate per-study correlation matrices ---
set.seed(123)
nStudy <- 12; nv <- 4

# True partial-correlation network (chain):
Omega <- matrix(0, nv, nv)
Omega[1, 2] <- Omega[2, 3] <- Omega[3, 4] <- 0.3
Omega <- Omega + t(Omega)
Delta <- diag(1 / sqrt(diag(solve(diag(nv) - Omega))))
trueSigma <- Delta %*% solve(diag(nv) - Omega) %*% Delta

cors <- list(); nobs <- numeric(nStudy)
for (s in seq_len(nStudy)) {
  nobs[s] <- round(runif(1, 200, 600))
  # study-specific population correlations (heterogeneity):
  P <- cov2cor(trueSigma +
               crossprod(matrix(rnorm(nv * nv, 0, 0.08), nv, nv)))
  X <- matrix(rnorm(nobs[s] * nv), nobs[s], nv) %*% chol(P)
  cors[[s]] <- cor(X)
  colnames(cors[[s]]) <- rownames(cors[[s]]) <- paste0("V", 1:nv)
}

# --- Fit and prune the meta-analytic GGM ---
mod_meta <- meta_ggm(cors = cors, nobs = nobs,
                     Vmethod = "OS_pooled") %>% runmodel

# Pooled network:
round(getmatrix(mod_meta, "omega_y"), 2)

# Between-study heterogeneity (SD per correlation):
round(sqrt(diag(getmatrix(mod_meta, "sigma_randomEffects"))), 3)

# Prune null edges:
mod_meta_pruned <- mod_meta %>% prune(alpha = 0.05, matrices = "omega_y")
compare(saturated = mod_meta, pruned = mod_meta_pruned)

The pooled network recovers the data-generating chain, and the three null edges are pruned without loss of fit (p = 0.78).

Since 0.16, Ising() accepts more than two response options, and since 0.16 the Blume-Capel model adds a per-node quadratic potential δ controlling the preference for neutral versus extreme responses. Here we simulate 5 ternary (−1/0/+1) variables from a true Blume-Capel chain network by enumerating all 3⁵ = 243 response patterns exactly, then compare the Ising and Blume-Capel fits. See the Ising & Blume-Capel page.

library(psychonetrics)
library(dplyr)

# --- Simulate from a true Blume-Capel model ---
set.seed(1)
nNode <- 5
resp  <- c(-1L, 0L, 1L)

tau   <- rep(0.1, nNode)    # thresholds
delta <- rep(0.5, nNode)    # quadratic potential (favours x = 0)
Omega <- matrix(0, nNode, nNode)               # chain network
Omega[1,2] <- Omega[2,3] <- Omega[3,4] <- Omega[4,5] <- 0.4
Omega <- Omega + t(Omega)

# Exact distribution over all 3^5 = 243 patterns:
states <- as.matrix(expand.grid(rep(list(resp), nNode)))
H <- states %*% tau - (states^2) %*% delta +
     rowSums((states %*% Omega) * states) / 2
P <- exp(H - max(H)); P <- P / sum(P)

Data <- as.data.frame(states[sample(nrow(states), 800,
                                    replace = TRUE, prob = P), ])
colnames(Data) <- paste0("X", 1:nNode)

# --- Fit multi-state Ising (delta fixed to 0) and Blume-Capel ---
mod_ising <- Ising(Data, responses = resp) %>% runmodel
mod_bc    <- BlumeCapel(Data, responses = resp) %>% runmodel

# The Blume-Capel model fits much better:
compare(ising = mod_ising, blume_capel = mod_bc)

# Recovered parameters (true: delta = 0.5, chain edges = 0.4):
round(as.vector(getmatrix(mod_bc, "delta")), 2)
round(getmatrix(mod_bc, "omega")[1, ], 2)

The Blume-Capel model clearly outperforms the plain multi-state Ising model on these data (ΔBIC ≈ 225) and recovers the simulated δ ≈ 0.5 and edge weights. Note that the exact ML estimator enumerates length(responses)^nNode patterns, so multi-state models are limited to small networks (see the maxStates argument).