Package 'psychonetrics'

Description

Multi-group (dynamical) structural equation models in combination with confirmatory network models from cross-sectional, time-series and panel data <doi:10.31234/osf.io/8ha93>. Allows for confirmatory testing and fit as well as exploratory model search.

Details

The DESCRIPTION file:

Index of help topics:

addMIs                  Model updating functions
aggregate_bootstraps    Aggregate Bootstrapped Models
allequal                Constrain all elements of a matrix equal within
                        groups
bifactor                Bi-factor models
BlumeCapel              Blume-Capel model
bootstrap               Bootstrap a psychonetrics model
changedata              Change the data of a psychonetrics object
checkJacobian           Diagnostic functions
CIplot                  Plot Analytic Confidence Intervals
compare                 Model comparison
covML                   Maximum likelihood covariance estimate
dlvm1                   Lag-1 dynamic latent variable model family of
                        psychonetrics models for panel data
duplicationMatrix       Model matrices used in derivatives
emergencystart          Reset starting values to simple defaults
equalityScoreTest       Score (Lagrange Multiplier) test for
                        cross-group equality constraints
esa                     Ergodic Subspace Analysis
factorscores            Compute factor scores
find_penalized_lambda   Automated Lambda Selection for Penalized
                        Estimation
fit                     Print fit indices
fixpar                  Parameters modification
fixstart                Attempt to Fix Starting Values
generate                Generate data from a fitted psychonetrics
                        object
getmatrix               Extract an estimated matrix
getVCOV                 Obtain the asymptotic covariance matrix
groupequal              Group equality constrains
Ising                   Ising model
Jonas                   Jonas dataset
latentgrowth            Latnet growth curve model
logbook                 Retrieve the psychonetrics logbook
loop_psychonetrics      Parallel loop for psychonetrics
lvm                     Continuous latent variable family of
                        psychonetrics models
meta_lvm                Meta-analytic latent variable model
meta_var1               Meta-analytic VAR(1) model
meta_varcov             Variance-covariance and GGM meta analysis
MIs                     Print modification indices
ml_lvm                  Multi-level latent variable model family
ml_tsdlvm1              Multi-level Lag-1 dynamic latent variable model
                        family of psychonetrics models for time-series
                        data
modelsearch             Stepwise model search
NA2020                  Network Analysis 2020 Self-Report Data
parameters              Print parameter estimates
parequal                Set equality constrains across parameters
partialprune            Partial pruning of multi-group models
penalize                Penalized Maximum Likelihood Functions
prune                   Stepdown model search by pruning
                        non-significant parameters.
psychonetrics_bootstrap-class
                        Class '"psychonetrics_bootstrap"'
psychonetrics_log-class
                        Class '"psychonetrics"'
psychonetrics-class     Class '"psychonetrics"'
psychonetrics-package   Structural Equation Modeling and Confirmatory
                        Network Analysis
ri_clpm                 Random intercept cross-lagged panel models
runmodel                Run a psychonetrics model
setestimator            Convenience functions
setverbose              Should messages of computation progress be
                        printed?
simplestructure         Generate factor loadings matrix with simple
                        structure
StarWars                Star Wars dataset
stepup                  Stepup model search along modification indices
transmod                Transform between model types
tsdlvm1                 Lag-1 dynamic latent variable model family of
                        psychonetrics models for time-series data
unionmodel              Unify models across groups
update_baseline         Set a custom baseline or saturated reference
                        model
var1                    Lag-1 vector autoregression family of
                        psychonetrics models
varcov                  Variance-covariance family of psychonetrics
                        models
write_psychonetrics     Write comprehensive model output to a text file

This package can be used to perform Structural Equation Modeling and confirmatory network modeling. Current implemented families of models are (1) the variance–covariance matrix (varcov), (2) the latent variable model (lvm), (3) the lag-1 vector autoregression model (var1), and (4) the dynamical lag-1 latent variable model for panel data (dlvm1) and for time-series data (tsdlvm1).

Author(s)

Sacha Epskamp [aut, cre]

Maintainer: Sacha Epskamp <[email protected]>

References

More information: psychonetrics.org

Description

Aggregates bootstrap results into a psychonetrics_bootstrap object

Usage

aggregate_bootstraps(sample, bootstraps, remove_problematic = TRUE)

Arguments

Details

After running this function, the helper functions parameters, fit, and CIplot can be used to investigate bootstrap results.

Value

An object of the class psychonetrics_bootstrap

Author(s)

Sacha Epskamp

Description

The allequal function constrains all free elements of a model matrix to be equal to one another within each group (a single shared parameter per group). In contrast to groupequal, which constrains corresponding elements equal across groups, allequal shares one parameter across the elements of the matrix and does so separately for each group, so the constraint is not imposed across groups. This is convenient for, for example, constraining all thresholds (tau) or all Blume-Capel quadratic potentials (delta) to a common value, or fitting an equal-edge-weight network (omega).

Usage

allequal(x, matrix, row, col, group, verbose, log = TRUE,
         runmodel = FALSE, identify = TRUE, ...)

Arguments

Details

Only free parameters are affected; elements that are fixed (e.g. fixed to zero) are left unchanged. For symmetric matrices the off-diagonal free elements are constrained equal. Within each group the selected free elements are given the same parameter number and a common starting value (the mean of their current values).

Value

An object of the class psychonetrics (psychonetrics-class)

Author(s)

Sacha Epskamp

See Also

groupequal, parequal, fixpar

Examples

## Not run: 
# Blume-Capel model with all quadratic (delta) potentials constrained equal,
# and all thresholds (tau) constrained equal, within the group:
mod <- BlumeCapel(data, responses = c(-1, 0, 1)) %>%
  allequal("delta") %>%
  allequal("tau") %>%
  runmodel

## End(Not run)

Description

Wrapper to lvm to specify a bi-factor model.

Usage

bifactor(data, lambda, latents, bifactor = "g", ...)

Arguments

Value

An object of the class psychonetrics (psychonetrics-class)

Author(s)

Sacha Epskamp

Description

The Blume-Capel model is the full Spin distribution model. In addition to the linear node potentials (tau) and pairwise interactions (omega) of the Ising model, it estimates a quadratic node potential delta (the Blume-Capel single-ion / crystal-field term). The Ising model is the special case obtained by fixing all delta parameters to zero, so the two models are nested. By convention the Blume-Capel model is defined for the ternary c(-1, 0, 1) response coding (a warning is raised for any other response set), but any number of ordered numeric response options is allowed (the same set for every variable, and not necessarily integers).

Usage

BlumeCapel(data, omega = "full", tau, delta, beta, beta_model =
                 c("beta", "log_beta"), vars, groups, covs, means,
                 nobs, covtype = c("choose", "ML", "UB"), responses,
                 missing = "listwise", equal = "none",
                 baseline_saturated = TRUE, estimator = "default",
                 optimizer, storedata = FALSE, WLS.W, sampleStats,
                 identify = TRUE, verbose = FALSE, maxNodes = 20,
                 maxStates = 2^maxNodes, min_sum = -Inf,
                 bootstrap = FALSE, boot_sub,
                 boot_resample, penalty_lambda = NA,
                penalty_alpha = 1, penalize_matrices)

Arguments

Details

The Blume-Capel model is built on the Spin distribution:

$\Pr(\boldsymbol{Y} = \boldsymbol{y}) = \frac{\exp\left( -\beta H\left(\boldsymbol{y}; \boldsymbol{\tau}, \boldsymbol{\delta}, \boldsymbol{\Omega}\right)\right)}{Z(\boldsymbol{\tau}, \boldsymbol{\delta}, \boldsymbol{\Omega})}$

With Hamiltonian:

$H\left(\boldsymbol{y}; \boldsymbol{\tau}, \boldsymbol{\delta}, \boldsymbol{\Omega}\right) = -\sum_{i=1}^{m} \tau_i y_{i} + \sum_{i=1}^{m} \delta_i y_{i}^2 - \sum_{i=2}^{m} \sum_{j=1}^{i-1} \omega_{ij} y_i y_j.$

and Z representing the partition function or normalizing constant. Equivalently,

$\Pr(\boldsymbol{Y} = \boldsymbol{y}) \propto \exp\left( \sum_i \tau_i y_i - \sum_i \delta_i y_i^2 + \sum_{i<j} \omega_{ij} y_i y_j \right).$

For the ternary states c(-1, 0, 1) this is the usual Blume-Capel parameterization: tau controls the tendency toward the positive versus the negative state, delta controls the tendency toward the zero/middle state versus the active/extreme states (positive delta favours the middle category), and omega controls pairwise alignment. For more general ordered numeric states the same parameterization defines an ordinal spin model.

The Ising model is nested within the Blume-Capel model by fixing all delta parameters to zero. Exact ML estimation enumerates every response pattern, so the cost grows as length(responses)^nNode (see maxStates); the model is therefore limited to relatively small networks.

Value

An object of the class psychonetrics

Author(s)

Sacha Epskamp <[email protected]>

References

Epskamp, S., Maris, G., Waldorp, L. J., & Borsboom, D. (2018). Network Psychometrics. In: Irwing, P., Hughes, D., & Booth, T. (Eds.), The Wiley Handbook of Psychometric Testing, 2 Volume Set: A Multidisciplinary Reference on Survey, Scale and Test Development. New York: Wiley.

See Also

Ising

Examples

library("dplyr")

# Simulate a small Blume-Capel data set (requires the development version of
# IsingSampler that supports the 'delta' argument):
## Not run: 
library("IsingSampler")
nNode  <- 5
graph  <- 0.4 * (matrix(1, nNode, nNode) - diag(nNode)) *
            (matrix(c(0,1,0,0,1, 1,0,1,0,0, 0,1,0,1,0, 0,0,1,0,1, 1,0,0,1,0), 5, 5))
data <- IsingSampler(1000, graph = graph, thresholds = rep(0, nNode),
                     beta = 1, delta = 1, responses = c(-1, 0, 1))

# Fit the Blume-Capel model:
mod <- BlumeCapel(data, responses = c(-1, 0, 1)) %>% runmodel

# Estimated network:
getmatrix(mod, "omega")

# Estimated quadratic node potentials:
getmatrix(mod, "delta")

## End(Not run)

Description

This function will bootstrap the data (once) and return a new unevaluated psychonetrics object. It requres storedata = TRUE to be used when forming a model.

Usage

bootstrap(x, replacement = TRUE, proportion = 1, verbose = TRUE, storedata = FALSE, 
          baseline_saturated = TRUE)

Arguments

Value

An object of the class psychonetrics (psychonetrics-class)

Author(s)

Sacha Epskamp

Description

This function can be used to change the data in a psychonetrics object.

Usage

changedata(x, data, covs, nobs, means, groups, missing = "listwise")

Arguments

Value

An object of the class psychonetrics (psychonetrics-class)

Author(s)

Sacha Epskamp

Description

Function to plot analytic confidence intervals (CI) of matrix elements estimated in psychonetrics.

Usage

CIplot(x, matrices, alpha_ci = 0.05, alpha_color = c(0.05,
                   0.01, 0.001, 1e-04), labels, labels2, labelstart,
                   print = TRUE, major_break = 0.2, minor_break = 0.1,
                   split0, prop0, prop0_cex = 1, prop0_alpha = 0.95,
                   prop0_minAlpha = 0.25)

Arguments

Value

A single ggplot2 object, or a list of ggplot2 objects for each matrix requested.

Author(s)

Sacha Epskamp

Examples

### Example from ?ggm ###
# Load bfi data from psych package:
library("psychTools")
data(bfi)

# Also load dplyr for the pipe operator:
library("dplyr")

# Let's take the agreeableness items, and gender:
ConsData <- bfi %>% 
  select(A1:A5, gender) %>% 
  na.omit # Let's remove missingness (otherwise use Estimator = "FIML)

# Define variables:
vars <- names(ConsData)[1:5]

# Let's fit an empty GGM:
mod0 <- ggm(ConsData, vars = vars)

# Run the model:
mod0 <- mod0 %>% runmodel

# Labels:
labels <- c(
  "indifferent to the feelings of others",
  "inquire about others' well-being",
  "comfort others",
  "love children",
  "make people feel at ease")

# Plot the CIs:
CIplot(mod0,  "omega", labels = labels, labelstart = 0.2)

### Example from ?gvar ###
library("dplyr")
library("graphicalVAR")

beta <- matrix(c(
  0,0.5,
  0.5,0
),2,2,byrow=TRUE)
kappa <- diag(2)
simData <- graphicalVARsim(50, beta, kappa)

# Form model:
model <- gvar(simData)

# Evaluate model:
model <- model %>% runmodel

# Plot the CIs:
CIplot(model,  "beta")

Description

This function will print a table comparing multiple models on chi-square, AIC and BIC.

Usage

compare(...)

## S3 method for class 'psychonetrics_compare'
print(x, ...)

Arguments

Value

A data frame with chi-square values, degrees of freedoms, RMSEAs, AICs, and BICs.

Author(s)

Sacha Epskamp

Description

These functions complement the base R cov function by simplifying obtaining maximum likelihood (ML) covariance estimates (denominator n) instead of unbiased (UB) covariance estimates (denominator n-1). The function covML can be used to obtain ML estimates, the function covUBtoML transforms from UB to ML estimates, and the function covMLtoUB transforms from UB to ML estimates.

Usage

covML(x, ...)
covUBtoML(x, n, ...)
covMLtoUB(x, n, ...)

Arguments

Author(s)

Sacha Epskamp <[email protected]>

Examples

data("StarWars")
Y <- StarWars[,1:10]

# Unbiased estimate:
UB <- cov(Y)

# ML Estimate:
ML <- covML(Y)

# Check:
all(abs(UB - covMLtoUB(ML, nrow(Y))) < sqrt(.Machine$double.eps))
all(abs(ML - covUBtoML(UB, nrow(Y))) < sqrt(.Machine$double.eps))

Description

The 'checkJacobian' function can be used to check if the analytic gradient / Jacobian is aligned with the numerically approximated gradient / Jacobian, and the 'checkFisher' function can be used to check if the analytic Hessian is aligned with the numerically approximated Hessian.

Usage

checkJacobian(x, f = "default", jac = "default", transpose = FALSE, 
          plot = TRUE, perturbStart = FALSE, method = "Richardson")

checkFisher(x, f = "default", fis = "default", transpose = FALSE, 
          plot = TRUE,  perturbStart = FALSE)

Arguments

Author(s)

Sacha Epskamp

Description

This is the family of models that models a dynamic factor model on panel data. All functions accept both wide-format data (with a design matrix for vars) and long-format data (with a character vector for vars plus idvar and optionally beepvar). The format is auto-detected from the type of vars. There are four covariance structures that can be modeled in different ways: within_latent, between_latent for the within-person and between-person latent (contemporaneous) models respectively, and within_residual, between_residual for the within-person and between-person residual models respectively. The panelgvar wrapper function sets the lambda to an identity matrix, all residual variances to zero, and models within-person and between-person latent (contemporaneous) models as GGMs. The panelvar wrapper does the same but models contemporaneous relations as a variance-covariance matrix. Finally, the panellvgvar wrapper automatically models all latent networks as GGMs. The old name panel_lvgvar is deprecated in favor of panellvgvar. For long-format data, pass the idvar (and optionally beepvar) arguments to panelgvar or panelvar via ...; see the examples below.

Usage

dlvm1(data, vars, datatype = c("auto", "wide", "long"),
                   idvar, beepvar,
                   standardize = c("none", "z", "quantile", "z_per_wave"),
                   lambda, within_latent = c("cov", "chol",
                   "prec", "ggm"), within_residual = c("cov", "chol",
                   "prec", "ggm"), between_latent = c("cov", "chol",
                   "prec", "ggm"), between_residual = c("cov", "chol",
                   "prec", "ggm"), beta = "full", omega_zeta_within =
                   "full", delta_zeta_within = "diag", kappa_zeta_within
                   = "full", sigma_zeta_within = "full",
                   lowertri_zeta_within = "full", omega_epsilon_within =
                   "zero", delta_epsilon_within = "diag",
                   kappa_epsilon_within = "diag", sigma_epsilon_within =
                   "diag", lowertri_epsilon_within = "diag",
                   omega_zeta_between = "full", delta_zeta_between =
                   "diag", kappa_zeta_between = "full",
                   sigma_zeta_between = "full", lowertri_zeta_between =
                   "full", omega_epsilon_between = "zero",
                   delta_epsilon_between = "diag", kappa_epsilon_between
                   = "diag", sigma_epsilon_between = "diag",
                   lowertri_epsilon_between = "diag", nu, mu_eta,
                   identify = TRUE, identification = c("loadings",
                   "variance"), latents, groups, groupvar, covs,
                   cors, means, nobs, corinput, start
                   = "version2", covtype = c("choose", "ML", "UB"),
                   missing = "auto", equal = "none",
                   baseline_saturated = TRUE, estimator = "ML",
                   optimizer, storedata = FALSE, verbose = FALSE,
                   sampleStats, baseline =
                   c("stationary_random_intercept", "stationary",
                   "independence", "none"), bootstrap = FALSE, boot_sub,
                   boot_resample, within, between,
                   penalty_lambda = NA, penalty_alpha = 1,
                   penalize_matrices)

panelgvar(data, vars, within_latent = c("ggm","chol","cov","prec"), 
          between_latent = c("ggm","chol","cov","prec"), ...)

panelvar(data, vars, within_latent = c("cov","chol","prec","ggm"), 
          between_latent = c("cov","chol","prec","ggm"), ...)

panellvgvar(...)

panel_lvgvar(...)

Arguments

Value

An object of the class psychonetrics (psychonetrics-class)

Author(s)

Sacha Epskamp

Examples

library("dplyr")

# Smoke data cov matrix, based on LISS data panel https://www.dataarchive.lissdata.nl
smoke <- structure(c(47.2361758611759, 43.5366809116809, 41.0057465682466, 
                     43.5366809116809, 57.9789886039886, 47.6992521367521, 
                     41.0057465682466, 
                     47.6992521367521, 53.0669434731935), .Dim = c(3L, 3L), 
                   .Dimnames = list(
                       c("smoke2008", "smoke2009", "smoke2010"), c("smoke2008", 
                   "smoke2009", "smoke2010")))

# Design matrix:
design <- matrix(rownames(smoke),1,3)

# Form model:
mod <- panelvar(vars = design, 
                covs = smoke, nobs = 352
)


# Run model:
mod <- mod %>% runmodel

# Evaluate fit:
mod %>% fit


# panelgvar and panelvar also accept long-format data.
# For long-format data, pass vars as a character vector
# and provide idvar (and optionally beepvar) via ...:

# Simulate some long-format panel data:
longdata <- data.frame(
  id = rep(1:50, each = 3),
  time = rep(1:3, 50),
  x = rnorm(150),
  y = rnorm(150)
)

# Use panelgvar with long-format data:
mod_long <- panelgvar(data = longdata,
                      vars = c("x", "y"),
                      idvar = "id",
                      beepvar = "time")

Description

These matrices are used in the analytic gradients

Usage

duplicationMatrix(n, diag = TRUE)

eliminationMatrix(n, diag = TRUE)

diagonalizationMatrix(n)

Arguments

Value

A sparse matrix

Author(s)

Sacha Epskamp

Examples

# Duplication matrix for 10 variables:
duplicationMatrix(10)

# Elimination matrix for 10 variables:
eliminationMatrix(10)

# Diagonailzation matrix for 10 variables:
diagonalizationMatrix(10)

Description

This function overwrites the starting values to simple defaults. This can help in cases where optimization fails.

Usage

emergencystart(x)

Arguments

Value

A psychonetrics model.

Author(s)

Sacha Epskamp

Description

Computes the score test (also known as the Lagrange Multiplier test) for releasing cross-group equality constraints in a multi-group psychonetrics model. Both the per-group univariate score statistics and the joint multivariate score statistic (one statistic per equality-constrained parameter, with df = G - 1) are returned.

The algorithm used by method = "jacobian" (the default) is a direct port of the score-test logic in lavaan::lavTestScore (Rosseel, 2012) and reproduces lavaan's output to numerical precision. The derivation closely follows the implementation in R/lav_test_score.R of the lavaan source code (https://github.com/yrosseel/lavaan), reused here with thanks to Yves Rosseel under the GPL ($\geq$ 2) license shared by both packages. The underlying gradient and Fisher-information machinery is provided by psychonetrics, so no refitting and no run-time dependence on lavaan is required.

With the default method = "jacobian" the algorithm reproduces lavaan::lavTestScore exactly: the model is augmented with G-1 extra parameters per equality-constrained tuple (one per non-reference group), the gradient and expected information are computed at the constrained MLE, and for each constraint the test statistic

$T = N \cdot s' \cdot \tilde{J}^{-1} \cdot s$

is formed, where $s$ is the gradient of the log-likelihood and $\tilde{J}^{-1}$ is the upper-left block of the generalised inverse of the bordered information matrix $[J\ R'^T;\ R\ 0]$, with $R$ the rows of the constraint Jacobian for the constraints that are NOT being tested. This is the same construction used by lavaan::lavTestScore.

method = "schur" uses a faster Schur-complement construction (the same used by psychonetrics for ordinary modification indices). It is much faster for complex models because it inverts only the (smaller) block of the information matrix corresponding to the currently free parameters, but it is not exactly equivalent to lavaan in unbalanced multi-group samples (it agrees in balanced samples). Both methods follow the same asymptotic chi-square reference distribution.

Usage

equalityScoreTest(x, matrices, top = 10, verbose = TRUE,
                  method = c("jacobian","schur"))

Arguments

Details

For two groups ($G = 2$) the joint test reduces to the univariate test (both have df = 1). The two tests differ for $G \ge 3$, where the joint test correctly accounts for the covariance between the released parameters via the (Schur-complemented) Fisher information, instead of summing or otherwise aggregating univariate statistics.

Value

Invisibly returns a list with two data frames:

uni

Univariate score tests, one row per released parameter (i.e. one row per group beyond the reference group), sorted from largest to smallest $X^2$.

total

Joint multivariate score tests, one row per equality-constrained parameter (i.e. one row per (matrix, row, column) triple), sorted from largest to smallest $X^2$. The joint test has df = G - 1 when the parameter is constrained equal across all $G$ groups.

Author(s)

Sacha Epskamp <[email protected]>. The method = "jacobian" implementation is a port of the score-test machinery in lavaan::lavTestScore by Yves Rosseel.

References

Rosseel, Y. (2012). lavaan: An R Package for Structural Equation Modeling. Journal of Statistical Software, 48(2), 1-36. doi:10.18637/jss.v048.i02

lavaan source code: https://github.com/yrosseel/lavaan (in particular R/lav_test_score.R).

See Also

MIs, partialprune

Examples

library("dplyr")
library("psychTools")
data(bfi)

ExData <- bfi %>% select(E1:E5, gender) %>% na.omit
vars <- names(ExData)[1:5]

mod <- ggm(ExData, vars = vars, groups = "gender") %>%
  runmodel %>%
  groupequal("omega") %>%
  runmodel

equalityScoreTest(mod)

Description

These functions implement Ergodic Subspace Analysis by von Oertzen, Schmiedek and Voelkle (2020). The functions can be used on the output of a dlvm1 model, or manually by supplying a within persons and between persons variance-covariance matrix.

Usage

esa(x, cutoff = 0.1,
    between = c("crosssection", "between"))
esa_manual(sigma_wp, sigma_bp, cutoff = 0.1)
## S3 method for class 'esa'
print(x, printref = TRUE, ...)
## S3 method for class 'esa_manual'
print(x, printref = TRUE, ...)
## S3 method for class 'esa'
plot(x, plot = c("observed", "latent"), ...)
## S3 method for class 'esa_manual'
plot(x,  ...)

Arguments

Value

For each group a esa_manual object with the following elements:

Author(s)

Sacha Epskamp <[email protected]>

References

von Oertzen, T., Schmiedek, F., and Voelkle, M. C. (2020). Ergodic Subspace Analysis. Journal of Intelligence, 8(1), 3.

Description

Currently, only the lvm framework with single group and no missing data is supported.

Usage

factorscores(data, model, method = c("bartlett", "regression"))

Arguments

Author(s)

Sacha Epskamp <[email protected]>

Description

Selects the optimal penalty strength (lambda) for penalized maximum likelihood (PML) or penalized full-information maximum likelihood (PFIML) models via EBIC-based grid search. A log-spaced grid of lambda values is searched from lambda_max (the smallest lambda for which all penalized parameters are exactly zero, derived from KKT conditions) down to lambda_max * lambda_min_ratio. Each lambda is optimized using ISTA (Iterative Shrinkage-Thresholding Algorithm) with warm starts from the previous solution, and the model with the best EBIC is selected. After selection, a beta_min threshold is applied to set near-zero penalized parameters exactly to zero.

This function is called automatically by runmodel when the model contains auto-select penalty parameters (i.e., penalty_lambda = NA in the parameter table, which is the default when using estimator = "PML" or estimator = "PFIML"). It can also be called directly for more control over the search.

Usage

find_penalized_lambda(model, criterion = "ebic.5", nLambda = 50,
                      lambda_min_ratio = 1e-4,
                      beta_min = c("theoretical", "numerical"),
                      verbose)

Arguments

Details

The search proceeds as follows:

1. Compute lambda_max: The analytical upper bound is derived from KKT (Karush-Kuhn-Tucker) conditions at the null-penalized model (all penalized parameters fixed at zero). Specifically, lambda_max = max(|grad_j|) / alpha where grad_j are the gradient elements for penalized parameters at the null solution and alpha is the elastic net mixing parameter.

2. Build lambda grid: A log-spaced sequence of nLambda values from lambda_max down to lambda_max * lambda_min_ratio.

3. Grid search with warm starts: For each lambda (from largest to smallest), the model is optimized using ISTA. The solution from the previous lambda serves as the starting point (warm start), which substantially reduces computation time. The EBIC (Extended Bayesian Information Criterion) is computed using the unpenalized log-likelihood (Convention A), counting only parameters with |estimate| >= beta_min as non-zero.

4. Beta-min thresholding: Penalized parameters with |estimate| < beta_min are set exactly to zero in the final model.

The EBIC is computed as:

$EBIC = -2 \cdot LL + k \cdot \log(n) + 4 \cdot k \cdot \gamma \cdot \log(p_v)$

where $LL$ is the unpenalized log-likelihood, $k$ is the effective number of parameters (unpenalized free parameters plus non-zero penalized parameters), $n$ is the sample size, $\gamma$ is the EBIC tuning parameter, and $p_v$ is the number of observed variables.

The returned model stores metadata about the search in model@optim$lambda_search, a list containing: criterion, gamma, best_lambda, best_criterion, beta_min, lambda_max, nLambda, and n_thresholded.

Value

An object of class psychonetrics (psychonetrics-class) with:

• The selected lambda assigned to all auto-select penalty parameters

• Parameters optimized at the best lambda

• Near-zero penalized parameters set to exactly zero (via beta_min)

• Search metadata stored in model@optim$lambda_search

After running find_penalized_lambda, use refit to obtain valid standard errors and fit indices via post-selection inference.

Author(s)

Sacha Epskamp

References

Foygel, R., & Drton, M. (2010). Extended Bayesian Information Criteria for Gaussian Graphical Models. In Advances in Neural Information Processing Systems (Vol. 23).

See Also

penalize for manually setting penalty parameters, runmodel for running models (calls find_penalized_lambda automatically when auto-select penalties are present), refit for post-selection inference after penalized estimation.

Examples

# Load bfi data from psych package:
library("psychTools")
library("dplyr")
data(bfi)

ConsData <- bfi %>%
  select(A1:A5) %>%
  na.omit

vars <- names(ConsData)

# Automatic lambda selection (called implicitly by runmodel):
mod <- ggm(ConsData, vars = vars, estimator = "PML")
mod <- mod %>% runmodel  # automatically searches for best lambda

# Check lambda search results:
mod@optim$lambda_search

# Post-selection refit for standard errors:
mod_refit <- mod %>% refit
mod_refit %>% parameters

# Direct call with custom criterion:
mod2 <- ggm(ConsData, vars = vars, estimator = "PML")
mod2 <- find_penalized_lambda(mod2, criterion = "ebic1",
                              nLambda = 100)

Description

This function will print all fit indices of the model/

Usage

fit(x)

Arguments

Value

Invisibly returns a data frame with fit measure estimates.

Author(s)

Sacha Epskamp

Examples

# Load bfi data from psych package:
library("psychTools")
data(bfi)

# Also load dplyr for the pipe operator:
library("dplyr")

# Let's take the agreeableness items, and gender:
ConsData <- bfi %>% 
  select(A1:A5, gender) %>% 
  na.omit # Let's remove missingness (otherwise use Estimator = "FIML)

# Define variables:
vars <- names(ConsData)[1:5]

# Let's fit an empty GGM:
mod0 <- ggm(ConsData, vars = vars, omega = "zero")

# Run model:
mod0 <- mod0 %>% runmodel

# Inspect fit:
mod0 %>% fit # Pretty bad fit...

Description

The fixpar function can be used to fix a parameter to some value (Typically zero), and the freepar function can be used to free a parameter from being fixed to a value.

Usage

fixpar(x, matrix, row, col, value = 0, group, verbose, 
        log = TRUE, runmodel = FALSE, ...)

freepar(x, matrix, row, col, start, group, verbose, log =
        TRUE, runmodel = FALSE, startEPC = TRUE, ...)

Arguments

Value

An object of the class psychonetrics (psychonetrics-class)

Author(s)

Sacha Epskamp

Description

This function attempts to fix starting values by comparing the analytic gradient to a numerically approximated gradient. Parameters with a difference between the analytic and numeric gradient that exceeds 'maxdiff' will be reduced by a factor of 'reduce' in each iteration until the average absolute difference between analytic and numeric gradients is lower than 'tol'. Only off-diagonal elements in omega, sigma, kappa, lowertri or rho matrices or any element in beta matrices are adjusted.

Usage

fixstart(x, reduce = 0.5, maxdiff = 0.1, tol = 0.01, maxtry = 25)

Arguments

Author(s)

Sacha Epskamp

Description

This function will generate new data from the estimated mean and variance-covariance structure of a psychonetrics model.

Usage

generate(x, n = 500)

Arguments

Value

A data frame with simulated data

Author(s)

Sacha Epskamp

Description

This function will extract an estimated matrix, and will either return a single matrix for single group models or a list of such matrices for multiple group models.

Usage

getmatrix(x, matrix, group, threshold = FALSE, alpha = 0.01,
           adjust = c("none", "holm", "hochberg", "hommel",
           "bonferroni", "BH", "BY", "fdr"), mode = c("tested",
           "all"), diag = TRUE)

Arguments

Value

A matrix of parameter estimates, of a list of such matrices for multiple group models.

Author(s)

Sacha Epskamp

Examples

# Load bfi data from psych package:
library("psychTools")
data(bfi)

# Also load dplyr for the pipe operator:
library("dplyr")

# Let's take the agreeableness items, and gender:
ConsData <- bfi %>% 
  select(A1:A5, gender) %>% 
  na.omit # Let's remove missingness (otherwise use Estimator = "FIML)

# Define variables:
vars <- names(ConsData)[1:5]

# Let's fit a full GGM:
mod <- ggm(ConsData, vars = vars, omega = "full")

# Run model:
mod <- mod %>% runmodel

# Obtain network:
mod %>% getmatrix("omega")

Description

This function can be used to obtain the estimated asymptotic covariance matrix from a psychonetrics object.

Usage

getVCOV(model, approximate_SEs = FALSE)

Arguments

Value

This function returns a matrix.

Author(s)

Sacha Epskamp

Description

The groupequal function constrains parameters equal across groups, and the groupfree function frees equality constrains across groups.

Usage

groupequal(x, matrix, row, col, verbose, log = TRUE, runmodel =
                    FALSE, identify = TRUE, ...)

groupfree(x, matrix, row, col, verbose, log = TRUE, runmodel =
                    FALSE, identify = TRUE, ...)

Arguments

Value

An object of the class psychonetrics (psychonetrics-class)

Author(s)

Sacha Epskamp

Description

This is the family of Ising models fit to datasets with two or more ordered integer response options. The classical case is a dichotomous dataset (e.g., encoded with -1 and 1 or with 0 and 1), but the same parameterization—one threshold (tau) per variable and a pairwise network (omega)—is also used for any number of ordered response options (e.g., c(-1,0,1) or seq(-5,5); the values need not be integers), encoded identically across all variables. Only the partition function changes: it sums over all length(responses)^nNode response patterns rather than only 2^nNode. Note that the input matters (see also https://arxiv.org/abs/1811.02916) in this model! Models based on a dataset that is encoded with -1 and 1 are not entirely equivalent to models based on datasets encoded with 0 and 1 (non-equivalences occur in multi-group settings with equality constrains).

Usage

Ising(data, omega = "full", tau, beta, beta_model =
                 c("beta", "log_beta"), vars, groups, covs, means,
                 nobs, covtype = c("choose", "ML", "UB"), responses,
                 missing = "listwise", equal = "none",
                 baseline_saturated = TRUE, estimator = "default",
                 optimizer, storedata = FALSE, WLS.W, sampleStats,
                 identify = TRUE, verbose = FALSE, maxNodes = 20,
                 maxStates = 2^maxNodes, min_sum = -Inf,
                 bootstrap = FALSE, boot_sub,
                 boot_resample, penalty_lambda = NA,
                penalty_alpha = 1, penalize_matrices)

Arguments

Details

The Ising Model takes the following form:

$\Pr(\boldsymbol{Y} = \boldsymbol{y}) = \frac{\exp\left( -\beta H\left(\boldsymbol{y}; \boldsymbol{\tau}, \boldsymbol{\Omega}\right)\right)}{Z(\boldsymbol{\tau}, \boldsymbol{\Omega})}$

With Hamiltonian:

$H\left(\boldsymbol{y}; \boldsymbol{\tau}, \boldsymbol{\Omega}\right) = -\sum_{i=1}^{m} \tau_i y_{i} - \sum_{i=2}^{m} \sum_{j=1}^{i-1} \omega_{ij} y_i y_j.$

And Z representing the partition function or normalizing constant.

The responses $y_i$ need not be dichotomous: they may take any number of ordered values (the same set for every variable, and not necessarily integers), supplied through responses or detected automatically from the data. The Hamiltonian above is unchanged; only the partition function $Z$ (and the expected values and Hessian derived from it) sums over all length(responses)^nNode response patterns instead of only 2^nNode. With two response options the model is the usual Ising model and reproduces the same estimates as before. With more than two response options the model remains a maximum-entropy distribution matching the means (via tau) and the pairwise products (via omega); it does not model higher-order interactions or the within-variable second moments, and the enumeration cost grows steeply with the number of nodes and response options (see maxStates).

Value

An object of the class psychonetrics

Author(s)

Sacha Epskamp <[email protected]>

References

Epskamp, S., Maris, G., Waldorp, L. J., & Borsboom, D. (2018). Network Psychometrics. In: Irwing, P., Hughes, D., & Booth, T. (Eds.), The Wiley Handbook of Psychometric Testing, 2 Volume Set: A Multidisciplinary Reference on Survey, Scale and Test Development. New York: Wiley.

See Also

BlumeCapel for the more general Blume-Capel model, which adds a quadratic node potential (delta) and reduces to the Ising model when all delta are fixed to zero.

Examples

library("dplyr")
data("Jonas")

# Variables to use:
vars <- names(Jonas)[1:10]

# Arranged groups to put unfamiliar group first (beta constrained to 1):
Jonas <- Jonas[order(Jonas$group),]

# Form saturated model:
model1 <- Ising(Jonas, vars = vars, groups = "group")

# Run model:
model1 <- model1 %>% runmodel
# Note: SEs are approximated automatically because there are zeroes
# in the crosstables of the "Knows Jonas" group, which makes the
# Fisher information matrix singular. A warning is shown when
# this fallback occurs.

# Prune-stepup to find a sparse model:
model1b <- model1 %>% prune(alpha = 0.05) %>%  stepup(alpha = 0.05)

# Equal networks:
model2 <- model1 %>% groupequal("omega") %>% runmodel

# Prune-stepup to find a sparse model:
model2b <- model2 %>% prune(alpha = 0.05) %>% stepup(mi = "mi_equal", alpha = 0.05)

# Equal thresholds:
model3 <- model2 %>% groupequal("tau") %>% runmodel

# Prune-stepup to find a sparse model:
model3b <- model3 %>% prune(alpha = 0.05) %>% stepup(mi = "mi_equal", alpha = 0.05)

# Equal beta:
model4 <- model3 %>% groupequal("beta") %>% runmodel

# Prune-stepup to find a sparse model:
model4b <- model4 %>% prune(alpha = 0.05) %>% stepup(mi = "mi_equal", alpha = 0.05)

# Compare all models:
compare(
  `1. all parameters free (dense)` = model1,
  `2. all parameters free (sparse)` = model1b,
  `3. equal networks (dense)` = model2,
  `4. equal networks (sparse)` = model2b,
  `5. equal networks and thresholds (dense)` = model3,
  `6. equal networks and thresholds (sparse)` = model3b,
  `7. all parameters equal (dense)` = model4,
  `8. all parameters equal (sparse)` = model4b
) %>% arrange(BIC)

Description

Responses of 10 attitude items towards a researcher named Jonas. Participants were shown three photos of Jonas with the text: "This is Jonas, a researcher from Germany who is now becoming a PhD in Psychology". Subsequently, the participants had to answer 10 yes / no questions starting with "I believe that Jonas...", as well as rate their familliarity with Jonas. The sample consists of people familiar with Jonas and not familiar with Jonas, and allows for testing Attitudinal Entropy Framework <doi:10.1080/1047840X.2018.1537246>.

Usage

data("Jonas")

Format

A data frame with 215 observations on the following 12 variables.

scientist

... is a good scientist

jeans

... Is a person that wears beautiful jeans

cares

... really cares about people like you

economics

... would solve our economic problems

hardworking

... is hardworking

honest

... is honest

intouch

... is in touch with ordinary people

knowledgeable

... is knowledgeable

makeupmind

... can't make up his mind

getsthingsdone

... gets things done

familiar

Answers to the question "How familiar are you with Jonas?" (three responses possible)

group

The question 'familiar' categorized in two groups ("Knows Jonas" and "Doesn't Know Jonas")

Examples

data(Jonas)

Description

Wrapper to lvm to specify a latent growth curve model.

Usage

latentgrowth(vars, time = seq_len(ncol(vars)) - 1, covariates =
                   character(0), covariates_as = c("regression",
                   "covariance"), ...)

Arguments

Details

See https://github.com/SachaEpskamp/SEM-code-examples/tree/master/Latent_growth_examples/psychonetrics for examples

Value

An object of the class psychonetrics (psychonetrics-class). See for an example https://github.com/SachaEpskamp/SEM-code-examples/tree/master/Latent_growth_examples/psychonetrics.

Author(s)

Sacha Epskamp

Examples

library("dplyr")

# Smoke data cov matrix, based on LISS data panel https://www.dataarchive.lissdata.nl
smoke <- structure(c(47.2361758611759, 43.5366809116809, 41.0057465682466, 
                     43.5366809116809, 57.9789886039886, 47.6992521367521, 
                     41.0057465682466, 
                     47.6992521367521, 53.0669434731935), .Dim = c(3L, 3L), 
                   .Dimnames = list(
                       c("smoke2008", "smoke2009", "smoke2010"), c("smoke2008", 
                   "smoke2009", "smoke2010")))

# Design matrix:
design <- matrix(rownames(smoke),1,3)

# Form model:
mod <- latentgrowth(vars = design, 
                covs = smoke, nobs = 352
)

## Not run: 
# Run model:
mod <- mod %>% runmodel

# Evaluate fit:
mod %>% fit

# Look at parameters:
mod %>% parameters

## End(Not run)

Description

This function can be used to retrieve the logbook of a 'psychonetrics' object.

Usage

logbook(x, log = TRUE)

Arguments

Author(s)

Sacha Epskamp

Description

A convenience wrapper for running repeated psychonetrics analyses in parallel. Replaces the typical boilerplate of makeCluster / clusterExport / parLapply / stopCluster with a single function call. Particularly useful for bootstrapping, simulation studies, and cross-validation.

Variables referenced in the expression that exist in the calling environment are automatically exported to parallel workers, so explicit clusterExport calls are not needed in most cases.

Usage

loop_psychonetrics(expr, reps = 1000, nCores = 1,
                   export = character(0),
                   packages = c("psychonetrics", "dplyr"),
                   verbose = TRUE, envir = parent.frame())

Arguments

Details

When nCores > 1, a PSOCK cluster is created and automatically cleaned up when the function exits (even on error). The specified packages are loaded on each worker, and variables referenced in expr are automatically detected and exported.

For bootstrapping, use this function together with bootstrap = "nonparametric" in the model constructor and aggregate_bootstraps to summarize the results.

Value

A list of length reps, where each element is the result of evaluating expr once.

See Also

aggregate_bootstraps, bootstrap, CIplot

Examples

library(dplyr)
data(NA2020)

# Fit the original model:
model <- ggm(NA2020, estimator = "FIML") %>% runmodel

# Bootstrap with pruning (parallel, 100 reps):
bootstraps <- loop_psychonetrics({
  ggm(NA2020, bootstrap = "nonparametric",
      estimator = "FIML") %>%
    runmodel %>%
    prune(alpha = 0.05)
}, reps = 100, nCores = 2)

# Aggregate and inspect:
boot_agg <- aggregate_bootstraps(
  sample = model %>% prune(alpha = 0.05),
  bootstraps = bootstraps
)

# View results:
parameters(boot_agg)
CIplot(boot_agg, "omega", split0 = TRUE)

Description

This is the family of models that models the data as a structural equation model (SEM), allowing the latent and residual variance-covariance matrices to be further modeled as networks. The latent and residual arguments can be used to define what latent and residual models are used respectively: "cov" (default) models a variance-covariance matrix directly, "chol" models a Cholesky decomposition, "prec" models a precision matrix, and "ggm" models a Gaussian graphical model (Epskamp, Rhemtulla and Borsboom, 2017). The wrapper lnm() sets latent = "ggm" for the latent network model (LNM), the wrapper rnm() sets residual = "ggm" for the residual network model (RNM), and the wrapper lrnm() combines the LNM and RNM.

Usage

lvm(data, lambda, latent = c("cov", "chol", "prec",
                   "ggm"), residual = c("cov", "chol", "prec", "ggm"),
                   sigma_zeta = "full", kappa_zeta = "full", omega_zeta =
                   "full", lowertri_zeta = "full", delta_zeta = "full",
                   sigma_epsilon = "diag", kappa_epsilon = "diag",
                   omega_epsilon = "zero", lowertri_epsilon = "diag",
                   delta_epsilon = "diag", beta = "zero", nu, nu_eta, tau,
                   identify = TRUE, identification = c("loadings",
                   "variance"), vars, ordered = character(0), latents,
                   groups, groupvar, covs, cors, means, nobs,
                   missing = "auto", equal = "none",
                   baseline_saturated = TRUE, estimator = "default",
                   optimizer, storedata = FALSE, WLS.W, covtype =
                   c("choose", "ML", "UB"), corinput, standardize = c("none", "z",
                   "quantile"), sampleStats, verbose = FALSE,
                   simplelambdastart = FALSE, start = "default",
                   bootstrap = FALSE,
                   boot_sub, boot_resample, penalty_lambda = NA,
                   penalty_alpha = 1, penalize_matrices)

lnm(...)
rnm(...)
lrnm(...)

Arguments

Details

The model used in this family is:

$\mathrm{var}( \boldsymbol{y} ) = \boldsymbol{\Lambda} (\boldsymbol{I} - \boldsymbol{B})^{-1} \boldsymbol{\Sigma}_{\zeta} (\boldsymbol{I} - \boldsymbol{B})^{-1\top} \boldsymbol{\Lambda}^{\top} + \boldsymbol{\Sigma}_{\varepsilon}$

$\mathcal{E}( \boldsymbol{y} ) = \boldsymbol{\nu} + \boldsymbol{\Lambda} (\boldsymbol{I} - \boldsymbol{B})^{-1} \boldsymbol{\nu}_eta$

in which the latent covariance matrix can further be modeled in three ways. With latent = "chol" as Cholesky decomposition:

$\boldsymbol{\Sigma}_{\zeta} = \boldsymbol{L}_{\zeta}\boldsymbol{L}_{\zeta}$,

with latent = "prec" as Precision matrix:

$\boldsymbol{\Sigma}_{\zeta} = \boldsymbol{K}_{\zeta}^{-1}$,

and finally with latent = "ggm" as Gaussian graphical model:

$\boldsymbol{\Sigma}_{\zeta} = \boldsymbol{\Delta}_{\zeta}(\boldsymbol{I} - \boldsymbol{\Omega}_{\zeta})^(-1) \boldsymbol{\Delta}_{\zeta}$.

Likewise, the residual covariance matrix can also further be modeled in three ways. With residual = "chol" as Cholesky decomposition:

$\boldsymbol{\Sigma}_{\varepsilon} = \boldsymbol{L}_{\varepsilon}\boldsymbol{L}_{\varepsilon}$,

with latent = "prec" as Precision matrix:

$\boldsymbol{\Sigma}_{\varepsilon} = \boldsymbol{K}_{\varepsilon}^{-1}$,

and finally with latent = "ggm" as Gaussian graphical model:

$\boldsymbol{\Sigma}_{\varepsilon} = \boldsymbol{\Delta}_{\varepsilon}(\boldsymbol{I} - \boldsymbol{\Omega}_{\varepsilon})^(-1) \boldsymbol{\Delta}_{\varepsilon}$.

Value

An object of the class psychonetrics (psychonetrics-class)

Author(s)

Sacha Epskamp

References

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

Muthen, B. (1984). A general structural equation model with dichotomous, ordered categorical, and continuous latent variable indicators. Psychometrika, 49(1), 115-132.

Examples

library("dplyr")

### Confirmatory Factor Analysis ###

# Example also shown in https://youtu.be/Hdu5z-fwuk8

# Load data:
data(StarWars)

# Originals only:
Lambda <- matrix(1,4)

# Model:
mod0 <- lvm(StarWars, lambda = Lambda, vars = c("Q1","Q5","Q6","Q7"), 
            identification = "variance", latents = "Originals")
            
# Run model:
mod0 <- mod0 %>% runmodel

# Evaluate fit:
mod0 %>% fit


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

# Observed variables:
obsvars <- paste0("Q",1:10)

# Latents:
latents <- c("Prequels","Original","Sequels")

# Make model:
mod1 <- lvm(StarWars, lambda = Lambda, vars = obsvars, 
            identification = "variance", latents = latents)

# Run model:
mod1 <- mod1 %>% runmodel

# Look at fit:
mod1

# Look at parameter estimates:
mod1 %>% parameters

# Look at modification indices:
mod1 %>% MIs

# Add and refit:
mod2 <- mod1 %>% freepar("sigma_epsilon","Q10","Q4") %>% runmodel

# Compare:
compare(original = mod1, adjusted = mod2)

# Fit measures:
mod2 %>% fit

### Path diagrams ###
# semPlot is not (yet) supported by default, but can be used as follows:
# Load packages:
if (requireNamespace("semPlot", quietly = TRUE)){
library("semPlot")

# Estimates:
lambdaEst <- getmatrix(mod2, "lambda")
psiEst <- getmatrix(mod2, "sigma_zeta")
thetaEst <- getmatrix(mod2, "sigma_epsilon")

# LISREL Model: LY = Lambda (lambda-y), TE = Theta (theta-epsilon), PS = Psi
mod <- lisrelModel(LY =  lambdaEst, PS = psiEst, TE = thetaEst)

# Plot with semPlot:
semPaths(mod, "std", "est", as.expression = "nodes")


# We can make this nicer (set whatLabels = "none" to hide labels):
semPaths(mod,

# this argument controls what the color of edges represent. In this case, 
# standardized parameters:
    what = "std", 
    
# This argument controls what the edge labels represent. In this case, parameter 
# estimates:
    whatLabels = "est", 
    
# This argument draws the node and edge labels as mathematical exprssions:    
    as.expression = "nodes", 
  
# This will plot residuals as arrows, closer to what we use in class:
    style = "lisrel",
    
# This makes the residuals larger:
    residScale = 10, 
    
# qgraph colorblind friendly theme:
    theme = "colorblind",
    
# tree layout options are "tree", "tree2", and "tree3":
    layout = "tree2", 

# This makes the latent covariances connect at a cardinal center point:
    cardinal = "lat cov",

# Changes curve into rounded straight lines:
    curvePivot = TRUE, 
    
# Size of manifest variables:
    sizeMan = 4, 
    
# Size of latent varibales:
    sizeLat = 10, 
    
# Size of edge labels:
    edge.label.cex = 1,
    
# Sets the margins:
    mar = c(9,1,8,1), 
    
# Prevents re-ordering of ovbserved variables:
    reorder = FALSE, 
    
# Width of the plot:
    width = 8, 
    
# Height of plot:
    height = 5, 

# Colors according to latents:
    groups = "latents",
    
# Pastel colors:
    pastel = TRUE, 
    
# Disable borders:
    borders = FALSE 
    )
    
# Use arguments filetype = "pdf" and filename = "semPlotExample1" to store PDF
} # end semPlot block

### Latent Network Modeling ###

# Latent network model:
lnm <- lvm(StarWars, lambda = Lambda, vars = obsvars,
           latents = latents, identification = "variance",
           latent = "ggm")

# Run model:
lnm <- lnm %>% runmodel

# Look at parameters:
lnm %>% parameters

# Remove non-sig latent edge:
lnm <- lnm %>% prune(alpha = 0.05)

# Compare to the original CFA model:
compare(cfa = mod1, lnm = lnm)

# Plot network:
library("qgraph")
qgraph(lnm@modelmatrices[[1]]$omega_zeta, labels = latents,
       theme = "colorblind", vsize = 10)

# A wrapper for the latent network model is the lnm function:
lnm2 <- lnm(StarWars, lambda = Lambda, vars = obsvars,
            latents = latents, identification = "variance")
lnm2 <- lnm2 %>% runmodel %>% prune(alpha = 0.05)
compare(lnm, lnm2) # Is the same as the model before.

# I could also estimate a "residual network model", which adds partial correlations to 
# the residual level:
# This can be done using lvm(..., residal = "ggm") or with rnm(...)
rnm <- rnm(StarWars, lambda = Lambda, vars = obsvars,
           latents = latents, identification = "variance")
# Stepup search:
rnm <- rnm %>% stepup

# It will estimate the same model (with link Q10 - Q4) as above. In the case of only one 
# partial correlation, There is no difference between residual covariances (SEM) or 
# residual partial correlations (RNM).


# For more information on latent and residual network models, see:
# Epskamp, S., Rhemtulla, M.T., & Borsboom, D. Generalized Network Psychometrics: 
# Combining Network and Latent Variable Models 
# (2017). Psychometrika. doi:10.1007/s11336-017-9557-x

### Gaussian graphical models ###

# All psychonetrics functions (e.g., lvm, lnm, rnm...) allow input via a covariance 
# matrix, with the "covs" and "nobs" arguments.
# The following fits a baseline GGM network with no edges:
S <- (nrow(StarWars) - 1)/ (nrow(StarWars)) * cov(StarWars[,1:10])
ggmmod <- ggm(covs = S, nobs = nrow(StarWars))

# Run model with stepup search and pruning:
ggmmod <- ggmmod%>% prune  %>% modelsearch

# Fit measures:
ggmmod %>% fit

# Plot network:
nodeNames <- c(
"I am a huge Star Wars\nfan! (star what?)",
"I would trust this person\nwith my democracy.",
"I enjoyed the story of\nAnakin's early life.",
"The special effects in\nthis scene are awful (Battle of\nGeonosis).",
"I would trust this person\nwith my life.",
"I found Darth Vader's big\nreveal in 'Empire' one of the greatest
moments in movie history.",
"The special effects in\nthis scene are amazing (Death Star\nExplosion).",
"If possible, I would\ndefinitely buy this\ndroid.",
"The story in the Star\nWars sequels is an improvement to\nthe previous movies.",
"The special effects in\nthis scene are marvellous (Starkiller\nBase Firing)."
)
library("qgraph")
qgraph(as.matrix(ggmmod@modelmatrices[[1]]$omega), nodeNames = nodeNames, 
    legend.cex = 0.25,  theme = "colorblind", layout = "spring")

# We can actually compare this model statistically (note they are not nested) to the 
# latent variable model:
compare(original_cfa = mod1, adjusted_cfa = mod2, exploratory_ggm = ggmmod)


### Meausrement invariance ###
# Let's say we are interested in seeing if people >= 30 like the original trilogy better 
# than people < 30.
# First we can make a grouping variable:
StarWars$agegroup <- ifelse(StarWars$Q12 < 30, "young", "less young")

# Let's look at the distribution:
table(StarWars$agegroup) # Pretty even...

# Observed variables:
obsvars <- paste0("Q",1:10)

# Let's look at the mean scores:
StarWars %>% group_by(agegroup) %>% summarize(across(all_of(obsvars), mean))
# Less young people seem to score higher on prequel questions and lower on other 
# questions

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

# Residual covariances:
Theta <- diag(1, 10)
Theta[4,10] <- Theta[10,4] <- 1

# Latents:
latents <- c("Prequels","Original","Sequels")

# Make model:
mod_configural <- lvm(StarWars, lambda = Lambda, vars = obsvars,
            latents = latents, sigma_epsilon = Theta,
            identification = "variance",
            groups =  "agegroup")

# Run model:
mod_configural <- mod_configural %>% runmodel

# Look at fit:
mod_configural
mod_configural %>% fit

# Looks good, let's try weak invariance:
mod_weak <- mod_configural %>% groupequal("lambda") %>% runmodel

# Compare models:
compare(configural = mod_configural, weak = mod_weak)

# weak invariance can be accepted, let's try strong:
mod_strong <- mod_weak %>% groupequal("nu") %>% runmodel
# Means are automatically identified

# Compare models:
compare(configural = mod_configural, weak = mod_weak, strong = mod_strong)

# Questionable p-value and AIC difference, but ok BIC difference. This is quite good, but 
# let's take a look. I have not yet implemented LM tests for equality constrains, but we 
# can look at something called "equality-free" MIs:
mod_strong %>% MIs(matrices = "nu", type = "free")

# Indicates that Q10 would improve fit. We can also look at residuals:
residuals(mod_strong)

# Let's try freeing intercept 10:
mod_strong_partial <- mod_strong %>% groupfree("nu",10) %>% runmodel

# Compare all models:
compare(configural = mod_configural,
        weak = mod_weak,
        strong = mod_strong,
        strong_partial = mod_strong_partial)

# This seems worth it and lead to an acceptable model! It seems that older people find 
# the latest special effects more marvellous!
mod_strong_partial %>% getmatrix("nu")

# Now let's investigate strict invariance:
mod_strict <- mod_strong_partial %>% groupequal("sigma_epsilon") %>% runmodel

# Compare all models:
compare(configural = mod_configural,
        weak = mod_weak,
        strong_partial = mod_strong_partial,
        strict = mod_strict)
# Strict invariance can be accepted!

#  Now we can test for homogeneity!
# Are the latent variances equal?
mod_eqvar <- mod_strict %>% groupequal("sigma_zeta") %>% runmodel

# Compare:
compare(strict = mod_strict, eqvar = mod_eqvar) 

# This is acceptable. What about the means? (alpha = nu_eta)
mod_eqmeans <- mod_eqvar %>% groupequal("nu_eta") %>% runmodel

# Compare:
compare(eqvar = mod_eqvar, eqmeans = mod_eqmeans)

# Rejected! We could look at MIs again:
mod_eqmeans %>% MIs(matrices = "nu_eta", type = "free")

# Indicates the strongest effect for prequels. Let's see what happens:
eqmeans2 <- mod_eqvar %>% 
  groupequal("nu_eta",row = c("Original","Sequels")) %>% runmodel

# Compare:
compare(eqvar = mod_eqvar, eqmeans = eqmeans2)
# Questionable, what about the sequels as well?

eqmeans3 <- mod_eqvar %>% groupequal("nu_eta", row = "Original") %>% runmodel

# Compare:
compare(eqvar = mod_eqvar, eqmeans = eqmeans3)

# Still questionable.. Let's look at the mean differences:
mod_eqvar %>% getmatrix("nu_eta")

# Looks like people over 30 like the prequels better and the other two trilogies less!

Description

Single-stage meta-analytic latent variable model (LVM) for covariance or correlation matrices. Combines a latent variable model (CFA/SEM) for the mean structure with a random-effects model for between-study heterogeneity. Based on the MAGNA framework (Epskamp et al., 2024).

Usage

meta_lvm(covs, nobs, data, cors, studyvar, groups, groupvar,
                   corinput, Vmats, Vmethod = c("individual", "pooled",
                   "metaSEM_individual", "metaSEM_weighted"), Vestimation
                   = c("averaged", "per_study"),
                   lambda, beta = "zero",
                   latent = c("cov", "chol", "prec", "ggm"),
                   sigma_zeta = "full", kappa_zeta = "full",
                   omega_zeta = "full", lowertri_zeta = "full",
                   delta_zeta = "full",
                   residual = c("cov", "chol", "prec", "ggm"),
                   sigma_epsilon = "diag", kappa_epsilon = "diag",
                   omega_epsilon = "zero", lowertri_epsilon = "diag",
                   delta_epsilon = "diag",
                   identify = TRUE,
                   identification = c("loadings", "variance"),
                   randomEffects = c("chol", "cov", "prec", "ggm", "cor"),
                   sigma_randomEffects = "full",
                   kappa_randomEffects = "full",
                   omega_randomEffects = "full",
                   lowertri_randomEffects = "full",
                   delta_randomEffects = "full",
                   rho_randomEffects = "full",
                   SD_randomEffects = "full",
                   vars, latents,
                   baseline_saturated = TRUE, optimizer,
                   estimator = c("FIML", "ML"),
                   sampleStats, verbose = FALSE,
                   bootstrap = FALSE, boot_sub, boot_resample)

Arguments

Details

This function implements a single-stage meta-analytic latent variable model. The model specifies that the mean of the vectorized correlation matrices follows an LVM-implied correlation structure:

$\mu = \textrm{vechs}(\Lambda (I-B)^{-1} \Sigma_\zeta (I-B)^{-\top} \Lambda' + \Sigma_\varepsilon)$

where diagonal elements of $\Sigma_\varepsilon$ are derived from the constraint that $\mu$ represents correlations (diagonal of implied covariance = 1).

Between-study heterogeneity is modeled as:

$\Sigma = \Sigma^{(\textrm{ran})} + V$

where $V$ is the known sampling error covariance and $\Sigma^{(\textrm{ran})}$ captures true between-study variation.

Value

An object of the class psychonetrics (psychonetrics-class)

Author(s)

Sacha Epskamp <[email protected]>

References

Epskamp, S. et al. (2024). Meta-Analytic Gaussian Network Aggregation. Psychometrika.

Jak, S., and Cheung, M. W. L. (2019). Meta-analytic structural equation modeling with moderating effects on SEM parameters. Psychological Methods.

See Also

meta_varcov, lvm

Examples

## Not run: 
# Generate simulated data
set.seed(42)
nStudies <- 10
nvar <- 6
nlat <- 2

# True factor loading matrix
lambda_true <- matrix(0, nvar, nlat)
lambda_true[1:3, 1] <- c(0.7, 0.8, 0.6)
lambda_true[4:6, 2] <- c(0.7, 0.8, 0.6)

# Generate correlation matrices per study
cors <- list()
nobs <- sample(100:400, nStudies)
for (i in 1:nStudies) {
  sigma <- lambda_true %*% t(lambda_true) + diag(1 - rowSums(lambda_true^2))
  dat <- MASS::mvrnorm(nobs[i], rep(0, nvar), sigma)
  cors[[i]] <- cor(dat)
  colnames(cors[[i]]) <- rownames(cors[[i]]) <- paste0("V", 1:nvar)
}

# Fit meta-analytic CFA
mod <- meta_lvm(covs = cors, nobs = nobs, lambda = lambda_true != 0)
mod <- mod %>% runmodel

# Inspect results
print(mod)
parameters(mod)
getmatrix(mod, "lambda")

## End(Not run)

Description

Single-stage meta-analytic vector autoregressive model (VAR(1)) for time-series data collected across multiple studies. Pools temporal (lag-1) and contemporaneous network structures across studies using a random-effects framework. The meta_gvar wrapper sets contemporaneous = "ggm" for graphical VAR estimation.

Usage

meta_var1(data, covs, nobs,
          vars, studyvar, idvar, dayvar, beepvar,
          contemporaneous = c("cov", "chol", "prec", "ggm"),
          beta = "full",
          omega_zeta = "full", delta_zeta = "full",
          kappa_zeta = "full", sigma_zeta = "full",
          lowertri_zeta = "full",
          randomEffects = c("chol", "cov", "prec", "ggm", "cor"),
          sigma_randomEffects = "full",
          kappa_randomEffects = "full",
          omega_randomEffects = "full",
          lowertri_randomEffects = "full",
          delta_randomEffects = "full",
          rho_randomEffects = "full",
          SD_randomEffects = "full",
          Vmats,
          Vmethod = c("individual", "pooled"),
          Vestimation = c("averaged", "per_study"),
          baseline_saturated = TRUE, optimizer,
          estimator = c("FIML", "ML"),
          sampleStats, verbose = FALSE,
          bootstrap = FALSE, boot_sub, boot_resample)

meta_gvar(...)

Arguments

Details

This function implements a single-stage meta-analytic VAR(1) model. For p observed variables, each study's time-series data is transformed into a Toeplitz covariance structure from which two blocks are extracted:

• Sigma_0: the p x p stationary covariance (symmetric, p(p+1)/2 unique elements)

• Sigma_1: the p x p lag-1 cross-covariance (non-symmetric, p^2 elements)

The VAR(1) structural model implies:

$\textrm{vec}(\Sigma_0) = (I - \beta \otimes \beta)^{-1} \textrm{vec}(\Sigma_\zeta)$

$\Sigma_1 = \beta \Sigma_0$

The meta-level mean is $\mu = [\textrm{vech}(\Sigma_0), \textrm{vec}(\Sigma_1)]$ and heterogeneity is modeled as $\Sigma = \Sigma^{(\textrm{ran})} + V$.

Note: the exogenous (lagged) covariance block from the full Toeplitz matrix is not modeled, as it is a nuisance parameter at the meta-analytic level.

Value

An object of the class psychonetrics (psychonetrics-class)

Author(s)

Sacha Epskamp <[email protected]>

See Also