var1: Lag-1 vector autoregression family of psychonetrics models

Description

This is the family of models that models time-series data using a lag-1 vector autoregressive model (VAR; Epskamp,Waldorp, Mottus, Borsboom, 2018). The model is fitted to the Toeplitz matrix, but unlike typical SEM software the block of covariances of the lagged variables is not used in estimating the temporal and contemporaneous relationships (the block is modeled completely separately using a cholesky decomposition, and does not enter the model otherwise). The contemporaneous argument can be used to define what contemporaneous model is used: contemporaneous = "cov" (default) models a variance-covariance matrix, contemporaneous = "chol" models a Cholesky decomposition, contemporaneous = "prec" models a precision matrix, contemporaneous = "ggm" (alias: gvar()) models a Gaussian graphical model, also then known as a graphical VAR model, and contemporaneous = "cor" models a correlation matrix with separate standard deviations.

Usage

var1(data, contemporaneous = c("cov", "chol", "prec",
                   "ggm", "cor"), beta = "full", omega_zeta = "full", delta_zeta
                   = "full", kappa_zeta = "full", sigma_zeta = "full",
                   lowertri_zeta = "full", mu, beepvar, dayvar, idvar,
                   vars, groups, groupvar, covs, cors, means,
                   nobs, corinput, missing = "auto",
                   centerWithin = FALSE,
                   equal = "none", baseline_saturated = TRUE, estimator =
                   "ML", optimizer, storedata = FALSE, covtype =
                   c("choose", "ML", "UB"), standardize = c("none", "z",
                   "quantile"), sampleStats, verbose = FALSE, bootstrap =
                   FALSE, boot_sub, boot_resample,
                   penalty_lambda = NA, penalty_alpha = 1,
                   penalize_matrices, rho_zeta = "full", SD_zeta = "full",
                   temporal = c("raw","PDC"), PDC = "full")
gvar(...)

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.
contemporaneous: The type of contemporaneous model used. See description.
beta: A model matrix encoding the temporal relationships (transpose of temporal network). A 0 encodes a fixed to zero element, a 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. Can also be "full" for a full temporal network or "zero" for an empty temporal network.
omega_zeta: Only used when contemporaneous = "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.
delta_zeta: Only used when contemporaneous = "ggm". The scaling (standard deviation) matrix, which is diagonal. Either "full" (the default) or "diag" to freely estimate all diagonal scalings, "zero" (not recommended) 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.
kappa_zeta: Only used when contemporaneous = "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.
sigma_zeta: Only used when contemporaneous = "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.
lowertri_zeta: Only used when contemporaneous = "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.
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.
beepvar: Optional string indicating assessment beep per day. Adding this argument will cause non-consecutive beeps to be treated as missing!
dayvar: Optional string indicating assessment day. Adding this argument makes sure that the first measurement of a day is not regressed on the last measurement of the previous day. IMPORTANT: only add this if the data has multiple observations per day.
idvar: Optional string indicating the subject ID
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, the matrix is treated as a covariance matrix with a warning (appropriate for standardized data). Requires nobs. Note: corinput = TRUE is not supported in var1.
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.
corinput: Logical. Not supported in var1() and will produce an error if set to TRUE.
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.
centerWithin: Logical, should the data be centered within persons (per idvar) before estimation? Defaults to FALSE (the historical behavior). Only used when raw data is supplied.
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. When missing = "auto" (default), "ML" is automatically switched to "FIML" and "PML" to "PFIML" if missing data is detected.
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.
sampleStats: An optional sample statistics object. Mostly used internally.
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.
verbose: Logical, should messages be printed?
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). NA (default) triggers automatic selection via EBIC-based grid search when a penalized estimator is used; set to a specific numeric value to use a fixed penalty strength (0 = no penalty). See find_penalized_lambda and penalize.
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.
rho_zeta: Only used when contemporaneous = "cor". Either "full" or "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. Encodes the contemporaneous (innovation) correlations. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
SD_zeta: Only used when contemporaneous = "cor". Either "full" or "diag", 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. Encodes the diagonal matrix of contemporaneous (innovation) standard deviations. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
temporal: Parameterization of the temporal effects: "raw" (default) models beta directly; "PDC" models the partial directed correlations directly, with the PDC matrix encoding from = row and to = column (the transposed orientation relative to beta, matching getmatrix(..., "PDC")). A (pattern-constrained) PDC model is fit-equivalent to the corresponding raw model, but standard errors, confidence intervals, equality constraints and modification indices then apply directly on the PDC scale.
PDC: Only used when temporal = "PDC". Model matrix for the partial directed correlations (from = row, to = column): "full", "zero", "diag", or a pattern matrix with 0s indicating elements fixed to zero, 1s indicating free elements, and higher integers indicating equality constraints.
...: Arguments sent to var1

Author

Sacha Epskamp

Details

The lag-1 vector autoregression (VAR(1)) model describes a multivariate time series through a lag-0 (contemporaneous) and a lag-1 (temporal) structure. The temporal structure is captured by the regression matrix beta (whose transpose is the temporal network), relating each observation to the previous measurement occasion. The contemporaneous structure captures the associations among the residuals (innovations) within the same measurement occasion, and can be parameterized in several ways through the contemporaneous argument: as a covariance matrix ("cov"), a Cholesky decomposition ("chol"), a precision matrix ("prec"), a Gaussian graphical model ("ggm"; the partial correlation network of the residuals), or a correlation matrix with separate standard deviations ("cor"). The gvar wrapper fixes contemporaneous = "ggm" to directly estimate a graphical VAR model. Combined, the temporal and contemporaneous networks imply the lag-0 and lag-1 (block-Toeplitz) covariance structure that is fitted to the data.

References

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

Examples

Run this code

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

# Parameter estimates:
model %>% parameters

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

# Note: this example is wrapped in a dontrun environment because the data is not 
# available locally.
if (FALSE) {
# Longer example:
#
# Obtain the data from:
#
# Epskamp, S., van Borkulo, C. D., van der Veen, D. C., Servaas, M. N., Isvoranu, A. M., 
# Riese, H., & Cramer, A. O. (2018). Personalized network modeling in psychopathology: 
# The importance of contemporaneous and temporal connections. Clinical Psychological 
# Science, 6(3), 416-427.
# 
# Available here: https://osf.io/c8wjz/

tsdata <- read.csv("Supplementary2_data.csv")

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

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

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

# Estimate, prune with FDR, and perform stepup search:
model_FDRprune <- gvar(
  tsdata, 
  vars = vars, 
  dayvar = "Day",
  estimator = "FIML"
  ) %>% 
  runmodel %>% 
  prune(adjust = "fdr", recursive = FALSE) %>% 
  stepup(criterion = "bic")

# Estimate with greedy stepup search:
model_stepup <- gvar(
  tsdata, 
  vars = vars, 
  dayvar = "Day",
  estimator = "FIML",
  omega_zeta = "zero",
  beta = "zero"
) %>% 
  runmodel %>% 
  stepup(greedy = TRUE, greedyadjust = "bonferroni", criterion = "bic")

# Compare models:
compare(
  FDRprune = model_FDRprune,
  stepup = model_stepup
)
# Very similar but not identical. Stepup is prefered here according to AIC and BIC

# Stepup results:
temporal <- getmatrix(model_stepup, "PDC") # PDC = Partial Directed Correlations
contemporaneous <- getmatrix(model_stepup, "omega_zeta")

# Average layout:
library("qgraph")
L <- averageLayout(temporal, contemporaneous)

# Labels:
labs <- gsub("\\.","\n",vars)

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

Run the code above in your browser using DataLab