varcov: Variance-covariance family of psychonetrics models

Description

This is the family of models that models only a variance-covariance matrix with mean structure. The type argument can be used to define what model is used: type = "cov" (default) models a variance-covariance matrix directly, type = "chol" (alias: cholesky()) models a Cholesky decomposition, type = "prec" (alias: precision()) models a precision matrix, type = "ggm" (alias: ggm()) models a Gaussian graphical model (Epskamp, Rhemtulla and Borsboom, 2017), and type = "cor" (alias: corr()) models a correlation matrix.

Usage

varcov(data, type = c("cov", "chol", "prec", "ggm", "cor"),
                   sigma = "full", kappa = "full", omega = "full",
                   lowertri = "full", delta = "diag", rho = "full", SD =
                   "full", mu, tau, vars, ordered = character(0), groups,
                   groupvar, covs, cors, means, nobs, missing = "auto", equal =
                   "none", baseline_saturated = TRUE, estimator =
                   "default", likelihood = c("normal", "wishart"),
                   fixed_x = character(0),
                   optimizer, storedata = FALSE, WLS.W,
                   sampleStats, meanstructure, corinput, verbose = FALSE,
                   covtype = c("choose", "ML", "UB"), standardize =
                   c("none", "z", "quantile"), fullFIML = FALSE,
                   bootstrap = FALSE, boot_sub, boot_resample,
                   penalty_lambda = NA, penalty_alpha = 1,
                   penalize_matrices, sampling_weights)
cholesky(...)
precision(...)
prec(...)
ggm(...)
corr(...)

Value

An object of the class psychonetrics

Arguments

data: A data frame encoding the data used in the analysis. Can be missing if covs and nobs are supplied.
type: The type of model used. See description.
sigma: Only used when type = "cov". Either "full" to estimate every element freely, "diag" to only include diagonal elements, or a matrix of the dimensions node x node with 0 encoding a fixed to zero element, 1 encoding a free to estimate element, and higher integers encoding equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
kappa: Only used when type = "prec". Either "full" to estimate every element freely, "diag" to only include diagonal elements, or a matrix of the dimensions node x node with 0 encoding a fixed to zero element, 1 encoding a free to estimate element, and higher integers encoding equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
omega: Only used when type = "ggm". Either "full" to estimate every element freely, "zero" to set all elements to zero, or a matrix of the dimensions node x node with 0 encoding a fixed to zero element, 1 encoding a free to estimate element, and higher integers encoding equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
lowertri: Only used when type = "chol". Either "full" to estimate every element freely, "diag" to only include diagonal elements, or a matrix of the dimensions node x node with 0 encoding a fixed to zero element, 1 encoding a free to estimate element, and higher integers encoding equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
delta: Only used when type = "ggm". Either "diag" or "zero" (not recommended), or a matrix of the dimensions node x node with 0 encoding a fixed to zero element, 1 encoding a free to estimate element, and higher integers encoding equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
rho: Only used when type = "cor". Either "full" to estimate every element freely, "zero" to set all elements to zero, or a matrix of the dimensions node x node with 0 encoding a fixed to zero element, 1 encoding a free to estimate element, and higher integers encoding equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
SD: Only used when type = "cor". Either "diag" or "zero" (not recommended), or a matrix of the dimensions node x node with 0 encoding a fixed to zero element, 1 encoding a free to estimate element, and higher integers encoding equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
mu: Optional vector encoding the mean structure. Set elements to 0 to indicate fixed to zero constrains, 1 to indicate free means, and higher integers to indicate equality constraints. For multiple groups, this argument can be a list or array with each element/column encoding such a vector.
tau: Optional list encoding the thresholds per variable.
vars: An optional character vector encoding the variables used in the analysis. Must equal names of the dataset in data.
groups: Deprecated. Use groupvar instead. An optional string indicating the name of the group variable in data.
groupvar: An optional string indicating the name of the group variable in data. Replaces the deprecated groups argument; if both are supplied, groupvar takes precedence with a warning.
covs: A sample variance--covariance matrix, or a list/array of such matrices for multiple groups. Make sure covtype argument is set correctly to the type of covariances used.
cors: A sample correlation matrix, or a list/array of such matrices for multiple groups. When supplied, corinput defaults to TRUE and the matrix is used in place of covs. Requires nobs.
means: A vector of sample means, or a list/matrix containing such vectors for multiple groups.
nobs: The number of observations used in covs and means, or a vector of such numbers of observations for multiple groups.
covtype: If 'covs' is used, this is the type of covariance (maximum likelihood or unbiased) the input covariance matrix represents. Set to "ML" for maximum likelihood estimates (denominator n) and "UB" to unbiased estimates (denominator n-1). The default will try to find the type used, by investigating which is most likely to result from integer valued datasets.
missing: How should missingness be handled in computing the sample covariances and number of observations when data is used. Can be "auto" (default) for automatic detection, "listwise" for listwise deletion, or "pairwise" for pairwise deletion. When "auto", the function checks for missing data and switches ML to FIML, PML to PFIML, or defaults to listwise for LS estimators.
equal: A character vector indicating which matrices should be constrained equal across groups.
baseline_saturated: A logical indicating if the baseline and saturated model should be included. Mostly used internally and NOT Recommended to be used manually.
estimator: The estimator to be used. Currently implemented are "ML" for maximum likelihood estimation, "FIML" for full-information maximum likelihood estimation, "PML" for penalized maximum likelihood estimation, "PFIML" for penalized full-information maximum likelihood estimation, "ULS" for unweighted least squares estimation, "WLS" for weighted least squares estimation, and "DWLS" for diagonally weighted least squares estimation. "WLSMV" is accepted as a synonym for "DWLS" (both use DWLS estimation with a mean-and-variance adjusted scaled test statistic). When missing = "auto" (default), "ML" is automatically switched to "FIML" and "PML" to "PFIML" if missing data is detected. Defaults to "ML" for continuous data and "DWLS" when ordinal variables are specified via ordered.
likelihood: The Gaussian likelihood scaling for maximum-likelihood estimation. "normal" (the default) uses the \(n\) (biased) sample-covariance denominator, so the chi-square is \(N \hat{F}\) and the parameter covariance is \(\mathrm{Info}^{-1}/N\). "wishart" uses the \(n-1\) (unbiased) sample covariance per group, so the chi-square uses \((n_g - 1)\) multipliers and the standard errors are inflated by \(\sqrt{n_g/(n_g-1)}\), matching lavaan's likelihood = "wishart" (Rosseel, 2012). Only available for complete-data maximum likelihood (estimator = "ML"); it errors for FIML, least-squares and ordinal estimators and for raw time-series input. Default "normal" reproduces the behaviour of previous versions exactly.
fixed_x: Character vector of exogenous variable names whose means and mutual (co)variances are fixed to their sample values and excluded from the free-parameter count and the degrees-of-freedom statistic count, matching lavaan's sem(..., fixed.x = TRUE) (Rosseel, 2012). The model is thereby conditioned on these variables, which relaxes the distributional (multivariate-normality) assumption on them; the cross-covariances between the fixed-x and the endogenous variables remain free. Only supported for type = "cov" and complete-data estimator = "ML". The reported log-likelihood is the conditional log-likelihood (of the endogenous variables given x). Default character(0) (no fixed.x; behaviour unchanged). The baseline and saturated reference models are conditioned on the exogenous block as well, so the incremental fit indices (CFI/TLI/NFI/...) and the absolute fit measures (chi-square, df, RMSEA, AIC, BIC, log-likelihood) and all estimates and standard errors match lavaan.
sampling_weights: Optional single column name in data giving (survey) SAMPLING weights, enabling pseudo-maximum-likelihood estimation. The moments are computed as weighted means/covariances (weights normalized so that they sum, over all groups, to the total sample size -- lavaan's sampling.weights.normalization = "total"), and the model is fit with robust (MLR) Huber-White sandwich standard errors and a Yuan-Bentler-Mplus scaled test statistic, exactly matching lavaan's sampling.weights=. Point estimates, standard errors, the scaled chi-square and the robust fit indices agree with lavaan to numerical precision (single- and multi-group). Requires complete-data continuous input (raw data); not supported for ordinal data, missing data (FIML), or likelihood = "wishart". References: Asparouhov (2005); Rosseel (2012).
optimizer: The optimizer to be used. Can be one of "nlminb" (the default R nlminb function), "ucminf" (from the optimr package), "nloptr_TNEWTON" (preconditioned truncated Newton via nloptr), and "LBFGS++" (pure C++ L-BFGS-B). Defaults to "nlminb".
storedata: Logical, should the raw data be stored? Needed for bootstrapping (see bootstrap).
standardize: Which standardization method should be used? "none" (default) for no standardization, "z" for z-scores, and "quantile" for a non-parametric transformation to the quantiles of the marginal standard normal distribution.
WLS.W: Optional WLS weights matrix.
sampleStats: An optional sample statistics object. Mostly used internally.
verbose: Logical, should progress be printed to the console?
ordered: A vector with strings indicating the variables that are ordered categorical, or set to TRUE to model all variables as ordered categorical.
meanstructure: Logical, should the meanstructure be modeled explicitly?
corinput: Logical, is the input a correlation matrix?
fullFIML: Logical, should row-wise FIML be used? Not recommended!
bootstrap: Should the data be bootstrapped? If TRUE the data are resampled and a bootstrap sample is created. These must be aggregated using aggregate_bootstraps! Can be TRUE or FALSE. Can also be "nonparametric" (which sets boot_sub = 1 and boot_resample = TRUE) or "case" (which sets boot_sub = 0.75 and boot_resample = FALSE).
boot_sub: Proportion of cases to be subsampled (round(boot_sub * N)).
boot_resample: Logical, should the bootstrap be with replacement (TRUE) or without replacement (FALSE)
penalty_lambda: Numeric penalty strength for penalized ML estimation (PML/PFIML). Default is NA, which triggers automatic lambda selection via EBIC grid search when estimator = "PML" or "PFIML" (see find_penalized_lambda). Set to a specific numeric value (e.g., 0.1) for manual lambda, or 0 for no penalty.
penalty_alpha: Elastic net mixing parameter: 1 = LASSO (default), 0 = ridge.
penalize_matrices: Character vector of matrix names to penalize. If missing, defaults are selected based on the model type.
...: Arguments sent to varcov

Author

Sacha Epskamp

Details

The model used in this family is:

\(\mathrm{var}(\boldsymbol{y} ) = \boldsymbol{\Sigma}\)

\(\mathcal{E}( \boldsymbol{y} ) = \boldsymbol{\mu}\)

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

\(\boldsymbol{\Sigma} = \boldsymbol{L}\boldsymbol{L}^{\top}\),

with type = "prec" as Precision matrix:

\(\boldsymbol{\Sigma} = \boldsymbol{K}^{-1}\),

with type = "ggm" as Gaussian graphical model:

\(\boldsymbol{\Sigma} = \boldsymbol{\Delta}(\boldsymbol{I} - \boldsymbol{\Omega})^(-1) \boldsymbol{\Delta}\),

and finally with type = "cor" as a correlation matrix scaled by standard deviations:

\(\boldsymbol{\Sigma} = \boldsymbol{D}\boldsymbol{R}\boldsymbol{D}\),

in which \(\boldsymbol{D}\) is a diagonal matrix of standard deviations and \(\boldsymbol{R}\) a correlation matrix.

References

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

Rosseel, Y. (2012). lavaan: An R package for structural equation modeling. Journal of Statistical Software, 48(2), 1--36. tools:::Rd_expr_doi("10.18637/jss.v048.i02")

Examples

Run this code

# 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 (or use missing = "auto" default which auto-selects FIML)

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

# Saturated estimation:
mod_saturated <- ggm(ConsData, vars = vars)

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

# We can look at the parameters:
mod_saturated %>% parameters

# Labels:
labels <- c(
  "indifferent to the feelings of others",
  "inquire about others' well-being",
  "comfort others",
  "love children",
  "make people feel at ease")
  
# Plot CIs:
CIplot(mod_saturated,  "omega", labels = labels, labelstart = 0.2)


# \donttest{
# We can also fit an empty network:
mod0 <- ggm(ConsData, vars = vars, omega = "zero")

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

# We can look at the modification indices:
mod0 %>% MIs

# To automatically add along modification indices, we can use stepup:
mod1 <- mod0 %>% stepup

# Let's also prune all non-significant edges to finish:
mod1 <- mod1 %>% prune

# Look at the fit:
mod1 %>% fit

# Compare to original (baseline) model:
compare(baseline = mod0, adjusted = mod1)

# We can also look at the parameters:
mod1 %>% parameters

# Or obtain the network as follows:
getmatrix(mod1, "omega")

# Penalized GGM estimation with automatic lambda selection:
mod_pml <- ggm(ConsData, vars = vars, estimator = "PML")
mod_pml <- mod_pml %>% runmodel

# Check selected lambda:
mod_pml@optim$lambda_search

# Obtain the sparse network:
getmatrix(mod_pml, "omega")

# Wishart likelihood (n-1 denominator), matching lavaan likelihood = "wishart":
mod_wis <- varcov(ConsData, vars = vars, type = "cov", likelihood = "wishart")
mod_wis <- mod_wis %>% runmodel

# fixed.x: condition on A1 and A2 as exogenous covariates (their means and
# mutual covariances are fixed to the sample values, excluded from npar/df):
mod_fx <- varcov(ConsData, vars = vars, type = "cov",
                 fixed_x = c("A1", "A2"))
mod_fx <- mod_fx %>% runmodel
# }

Run the code above in your browser using DataLab