This is the family of models that models a lag-1 vector autoregression (VAR) process with random intercepts on panel data of observed variables: a within-person stationary VAR(1) process (temporal matrix beta and contemporaneous structure sigma_zeta_within) combined with a between-person (random intercept) structure (sigma_zeta_between). The panelvar function models the contemporaneous and between-person structures as variance--covariance matrices by default, and the panelgvar wrapper models them as Gaussian graphical models (GGMs), yielding the panel graphical VAR model.
This model family is the observed-variable special case of dlvm1 (a factor loading matrix fixed to identity, zero residual variances and observed means). Since psychonetrics 0.16.2 it is implemented as its own model framework: results are identical to routing through dlvm1 with dummy matrices, but estimation skips all factor-loading algebra and is faster. The mean structure is a single mu matrix of observed stationary means (in older versions these appeared as nu); with stationary_means = FALSE mu instead holds wave-specific means, one per variable per wave (the trend-stationary case, see the details below). Latent panel models (lambda not identity) are still handled by dlvm1 and its wrapper panellvgvar.
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.
panelvar(data, vars, within_latent = c("cov", "chol", "prec",
"ggm", "cor"), between_latent = c("cov", "chol",
"prec", "ggm", "cor"), datatype = c("auto", "wide",
"long"), idvar, beepvar, standardize = c("none",
"z", "quantile", "z_per_wave"), beta = "full",
omega_zeta_within = "full", delta_zeta_within =
"diag", kappa_zeta_within = "full",
sigma_zeta_within = "full", lowertri_zeta_within =
"full", omega_zeta_between = "full",
delta_zeta_between = "diag", kappa_zeta_between =
"full", sigma_zeta_between = "full",
lowertri_zeta_between = "full", mu, 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,
rho_zeta_within = "full", SD_zeta_within = "diag",
rho_zeta_between = "full", SD_zeta_between = "diag",
temporal = c("raw","PDC"), PDC = "full",
stationary_means = TRUE)panelgvar(data, vars, within_latent = c("ggm","chol","cov","prec","cor"),
between_latent = c("ggm","chol","cov","prec","cor"), ...)
An object of the class psychonetrics (psychonetrics-class)
A data frame encoding the data used in the analysis. Can be missing if covs and nobs are supplied.
Required argument. For wide-format data, this must be a *matrix* with each row indicating a variable and each column indicating a measurement. The matrix must be filled with names of the variables in the dataset corresponding to variable i at wave j. NAs can be used to indicate missing waves. The rownames of this matrix will be used as variable names. For long-format data, this should be a character vector of variable names in the dataset.
The type of within-person contemporaneous model to be used. Since there are no latent variables in panelvar, the shorter alias within is preferred.
The type of between-person model to be used. Since there are no latent variables in panelvar, the shorter alias between is preferred.
Data format: "auto" (default) auto-detects based on whether vars is a matrix (wide) or character vector (long). "wide" forces wide-format interpretation. "long" forces long-format interpretation with automatic reshape to wide.
Optional string indicating the subject/cluster ID variable in data. Required when datatype = "long" or when vars is a character vector.
Optional string indicating the time point / measurement occasion variable. If missing when datatype = "long", a sequential measurement variable is created per subject.
Standardization method: "none" (default), "z" for global z-scores per variable across all waves (uses the overall mean and standard deviation for each variable computed across all waves of data), "quantile" for quantile normalization across all waves, "z_per_wave" for independent z-scores per wave column (note: this imposes stationarity as means become 0 and variances become 1 at each wave).
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.
Only used when within_latent = "ggm". Can be "full", "zero", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
Only used when within_latent = "ggm". Can be "diag", "zero" (not recommended), or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
Only used when within_latent = "prec". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
Only used when within_latent = "cov". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
Only used when within_latent = "chol". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
Only used when between_latent = "ggm". Can be "full", "zero", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
Only used when between_latent = "ggm". Can be "diag", "zero" (not recommended), or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
Only used when between_latent = "prec". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
Only used when between_latent = "cov". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
Only used when between_latent = "chol". Can be "full", "diag", or a typical model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this argument can be a list or array with each element/slice encoding such a matrix.
Optional vector encoding the observed stationary means. 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. When stationary_means = FALSE this is instead a nVar by nTime matrix (variables in the rows, waves in the columns), or an array/list of such matrices for multiple groups; elements at variable--wave combinations that are not in the vars design matrix are always fixed to zero.
Deprecated. Use groupvar instead. An optional string indicating the name of the group variable in data.
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.
A sample variance--covariance matrix, or a list/array of such matrices for multiple groups. Make sure the covtype argument is set correctly to the type of covariances used.
Not supported in panelvar(); supplying correlation input produces an error.
A vector of sample means, or a list/matrix containing such vectors for multiple groups. This argument always refers to the sample means; the mean structure of the model is set with stationary_means.
The number of observations used in covs and means, or a vector of such numbers of observations for multiple groups.
Logical. Not supported in panelvar() and will produce an error if set to TRUE.
Start value specification. Can be either a string or a psychonetrics model. If it is a string, "version2" indicates the latest version of start value computation, "version1" indicates start values as they were computed up to version 0.11, and "simple" indicate simple starting values. If this is a psychonetrics model the starting values will be based on the output. This can be useful, for example, if you first estimate a model with matrices set to a Cholesky decomposition, then use those values as start values for estimating Gaussian graphical models.
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.
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.
A character vector indicating which matrices should be constrained equal across groups.
A logical indicating if the baseline and saturated model should be included. Mostly used internally and NOT Recommended to be used manually.
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.
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".
Logical, should the raw data be stored? Needed for bootstrapping (see bootstrap).
Logical, should progress be printed to the console?
An optional sample statistics object. Mostly used internally.
What baseline model should be used? "stationary_random_intercept" includes both within- and between person variances constrained equal across time (default), "stationary" only includes within-person variances constrained equal across time, "independence" (default up to version 0.11) includes a variance for every variable at every time point (not constrained equal across time), and "none" includes no baseline model. The baseline model always uses the same mean structure as the model itself (see the stationary_means argument and the details below).
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).
Proportion of cases to be subsampled (round(boot_sub * N)).
Logical, should the bootstrap be with replacement (TRUE) or without replacement (FALSE)
Preferred name for within_latent (there are no latent variables in panelvar). If both within and within_latent are explicitly specified, within_latent takes precedence and a warning is issued.
Preferred name for between_latent (there are no latent variables in panelvar). If both between and between_latent are explicitly specified, between_latent takes precedence and a warning is issued.
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.
Elastic net mixing parameter: 1 = LASSO (default), 0 = ridge.
Character vector of matrix names to penalize. If missing, defaults are selected based on the model type.
Only used when within_latent = "cor". Correlation matrix specification ("full", "zero", or a pattern matrix; see omega_zeta_within).
Only used when within_latent = "cor". Diagonal standard-deviation matrix specification ("diag", "full", or a pattern matrix; see delta_zeta_within).
Only used when between_latent = "cor". Correlation matrix specification ("full", "zero", or a pattern matrix).
Only used when between_latent = "cor". Diagonal standard-deviation matrix specification ("diag", "full", or a pattern matrix).
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.
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.
Logical. TRUE (default, and the behavior of all earlier versions) models a single stationary mean per variable, constrained equal across waves. FALSE models a separate mean for every variable at every wave (the mu matrix becomes nVar by nTime), which saturates the mean structure and adds nVar * (nTime - 1) parameters. See the details below.
Arguments sent to panelvar.
Sacha Epskamp
Stationary versus wave-specific means. By default the panel VAR model is fully stationary: every variable has a single mean, and the implied means are therefore equal across waves. With stationary_means = FALSE the means become wave-specific, which is the classical trend-stationary case. The level shifts are then period effects that are exogenous to the dynamics: they are deterministic offsets added to all cases at a given wave, the within-person deviation process
$$y_{i,t} = \mu_t + \eta_{i,t} + \zeta_i, \quad \eta_{i,t} = B \eta_{i,t-1} + \varepsilon_{i,t}$$
remains covariance-stationary, and wave 1 is still drawn from that stationary distribution. Consequently the implied covariance matrix is identical to that of the stationary-means model, and only the mean part of the model (and of its Jacobian) changes. Fitting stationary_means = FALSE to raw data is equivalent to fitting the default model to per-wave mean-centered data, except that the free means are estimated (and tested) rather than removed, and the degrees of freedom are adjusted accordingly.
Identification. The mean structure under stationary_means = FALSE is exactly saturated: one free mean per variable--wave combination present in the vars design matrix (variable--wave combinations that are absent have no mean parameter). The estimated means therefore equal the observed per-wave means, and the likelihood-ratio test of stationary_means = TRUE against FALSE is a test of equal means over time with nVar * (nTime - 1) degrees of freedom. Cross-group equality constraints (equal = "mu") equate the means per wave across groups.
No free innovation variances. Only the means can be freed over time in this framework, not the innovation (co)variances: a wave-varying innovation covariance has no stationary distribution for wave 1 to be drawn from, so the Kronecker/stationarity machinery on which panelvar is built does not apply. Models with wave-specific innovation structures are covered by the RI-CLPM / RI-CLPN framework instead, see ri_clpm (which also leaves the exogenous wave-1 covariance block free).
Baseline model. The baseline model always uses the same mean structure as the fitted model: with stationary_means = FALSE the baseline gets wave-specific means as well. A baseline model must be nested in the fitted model, otherwise it can fit better than the fitted model on the mean side (a free-means baseline against an equal-means target would beat it), which makes incremental fit indices meaningless or negative. A corollary is that incremental fit indices (CFI, TLI, NFI, ...) are not directly comparable between a stationary_means = TRUE and a stationary_means = FALSE fit of the same data, because the two use different baselines by design; compare those two models with a likelihood-ratio test (compare) instead. The saturated model already has a saturated mean structure and is unchanged.
dlvm1 for the latent variable generalization of this model family.
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
)
# \donttest{
# 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):
# \donttest{
# 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")
# The same model with wave-specific (trend-stationary) means, which can be
# compared to the stationary-means model with a chi-square difference test:
mod_free <- panelgvar(data = longdata,
vars = c("x", "y"),
idvar = "id",
beepvar = "time",
stationary_means = FALSE)
# }
Run the code above in your browser using DataLab