This family models a multi-level (random intercept) lag-1 vector auto-regression
(VAR) / graphical VAR for intensive longitudinal data (many repeated measures nested
in subjects), estimated by a two-level summary-statistics pseudo maximum likelihood.
Each subject follows a within-person stationary VAR(1) process (temporal matrix
beta and contemporaneous structure sigma_zeta_within) around a
person-specific random intercept, and the random intercepts have a between-person
covariance structure (sigma_zeta_between). The ml_var1 function models
the contemporaneous and between-person structures as variance--covariance matrices by
default, and the ml_gvar1 wrapper models them as Gaussian graphical models
(GGMs), yielding the three multi-level GVAR networks: the temporal network
(beta / partial directed correlations), the contemporaneous network
(omega_zeta_within) and the between-person network (omega_zeta_between).
This model targets the same population model as panelvar /
panelgvar (a multi-level VAR with random intercepts), but is designed for
intensive longitudinal / experience-sampling (ESM) data, where the number of repeated
measures per subject is large and the full-likelihood panelvar route is
infeasible. See Details for the estimator and its statistical caveats.
ml_var and ml_gvar are the primary names of this family
(ml_var1 / ml_gvar1 are identical and remain fully supported; the
historical "1" referred to the lag order). In psychonetrics 0.15 ml_var /
ml_gvar were wrappers around dlvm1 with
estimator = "FIML"; the default estimator = "auto" reproduces that
behavior for short (panel-like) measurement-occasion grids (at most 10 occasions,
now through the faster panelvar framework with identical estimates)
and uses the two-level summary-statistics estimator for longer (ESM-like) series,
where full-information ML is not feasible. Their argument names
contemporaneous and between are mapped to within_latent and
between_latent.
This estimator is experimental (a message is shown once per session).
ml_var1(data, idvar, vars, dayvar, beepvar, groups, groupvar,
within_latent = c("cov", "chol", "prec", "ggm", "cor"),
between_latent = c("cov", "chol", "prec", "ggm", "cor"),
beta = "full",
omega_zeta_within = "full", delta_zeta_within = "diag",
kappa_zeta_within = "full", sigma_zeta_within = "full",
lowertri_zeta_within = "full", rho_zeta_within = "full",
SD_zeta_within = "full",
omega_zeta_between = "full", delta_zeta_between = "diag",
kappa_zeta_between = "full", sigma_zeta_between = "full",
lowertri_zeta_between = "full", rho_zeta_between = "full",
SD_zeta_between = "full",
mu, equal = "none", baseline_saturated = TRUE,
estimator = c("auto", "ML", "FIML"),
optimizer, storedata = FALSE,
standardize = c("none", "z", "quantile"),
sampleStats, verbose = FALSE, toeplitz = TRUE,
within, between)ml_gvar1(..., within_latent = "ggm", between_latent = "ggm")
ml_var(..., contemporaneous = c("cov", "chol", "prec", "ggm", "cor"),
between = c("cov", "chol", "prec", "ggm", "cor"))
ml_gvar(..., contemporaneous = c("ggm", "cov", "chol", "prec", "cor"),
between = c("ggm", "cov", "chol", "prec", "cor"))
An object of the class psychonetrics (psychonetrics-class).
A data frame encoding the (long-format) data used in the analysis. Must be a raw dataset with one row per measurement occasion, containing the subject id variable (idvar) and the variables in vars.
Required string indicating the subject/cluster ID variable in data. Each subject's rows form one within-person time series; the random intercepts vary over subjects. (Use var1 for single-subject time-series models.)
An optional character vector encoding the variables used in the analysis. Must equal names in data. If missing, all columns except idvar, dayvar, beepvar and the group variable are used.
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.
Optional string indicating assessment beep per day. Adding this argument will cause non-consecutive beeps to be treated as missing (the corresponding lag pairs are dropped).
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 for multi-group analysis. Replaces the deprecated groups argument.
The type of within-person contemporaneous model to be used: "cov" (covariance, default), "chol" (Cholesky), "prec" (precision), "ggm" (Gaussian graphical model, giving the contemporaneous network omega_zeta_within), or "cor" (correlation).
The type of between-person (random-intercept) model to be used: "cov" (default), "chol", "prec", "ggm" (giving the between-person network omega_zeta_between), or "cor".
A model matrix encoding the temporal relationships (transpose of the temporal network; from = row, to = column, as in var1). A 0 encodes a fixed to zero element, a 1 a free to estimate element, and higher integers equality constraints. For multiple groups, this can be a list or array. Can also be "full" (default) for a full temporal network or "zero" for an empty temporal network.
Only used when within_latent = "ggm". Can be "full", "zero", or a model matrix with 0s indicating parameters constrained to zero, 1s indicating free parameters, and higher integers indicating equality constraints. For multiple groups, this can be a list or array.
Only used when within_latent = "ggm". Scaling (delta) matrix. Can be "diag", "zero" (not recommended), or a model matrix.
Only used when within_latent = "prec". Precision matrix. Can be "full", "diag", or a model matrix.
Only used when within_latent = "cov". Covariance matrix. Can be "full", "diag", or a model matrix.
Only used when within_latent = "chol". Cholesky (lower triangular) matrix. Can be "full", "diag", or a model matrix.
Only used when within_latent = "cor". Correlation matrix specification ("full", "zero", or a pattern matrix).
Only used when within_latent = "cor". Diagonal standard-deviation matrix specification ("full", "diag", or a pattern matrix).
Only used when between_latent = "ggm". Between-person network. Can be "full", "zero", or a model matrix.
Only used when between_latent = "ggm". Scaling (delta) matrix. Can be "diag", "zero" (not recommended), or a model matrix.
Only used when between_latent = "prec". Precision matrix. Can be "full", "diag", or a model matrix.
Only used when between_latent = "cov". Covariance matrix. Can be "full", "diag", or a model matrix.
Only used when between_latent = "chol". Cholesky (lower triangular) matrix. Can be "full", "diag", or a model matrix.
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 ("full", "diag", or a pattern matrix).
Optional vector encoding the stationary means. Set elements to 0 to indicate fixed to zero constraints, 1 to indicate free means, and higher integers to indicate equality constraints. For multiple groups, this can be a list or matrix.
A character vector indicating which matrices should be constrained equal across groups (e.g. "beta").
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. The default "auto" selects "FIML" when the measurement-occasion grid is short (at most 10 occasions per subject, the panel-data regime, where full-information ML is feasible and matches the behavior of ml_var in psychonetrics 0.15) and "ML" otherwise, with a message stating the choice. "ML" uses the two-level summary-statistics pseudo maximum likelihood described in Details. "FIML" dispatches to the panelvar framework: the measurement occasions become the waves of a panel model that is estimated with full-information maximum likelihood (missing responses and missed beeps are handled by FIML; a night break, when dayvar is used, is bridged by one structurally missing occasion, so the model links the last beep of a day to the first beep of the next day through two VAR steps). The FIML route is exact but its cost grows with the number of measurement occasions, so it is only feasible for a modest number of occasions per subject.
The optimizer to be used. Can be one of "nlminb" (the default R nlminb function), "ucminf", "nloptr_TNEWTON", and "LBFGS++".
Logical, should the raw data be stored?
Which standardization method should be used on the raw variables before the lag embedding? "none" (default), "z" for z-scores, and "quantile" for a non-parametric quantile transformation. Standardization is grand (over all observations), not within-person, so that the between-person variance survives.
An optional sample statistics object. Mostly used internally (for the recursive baseline/saturated calls).
Logical, should progress be printed to the console?
Internal. If TRUE (default) the stationary Toeplitz VAR(1) within structure is imposed; FALSE selects the free-within (saturated) mode used internally to build the reference model. Not intended for direct use.
Arguments sent to ml_var1 (used by ml_gvar1, ml_var and ml_gvar).
Alias for within_latent (there are no latent variables in this model, so the shorter name is preferred).
Alias for between_latent; also the name used by ml_var and ml_gvar since psychonetrics 0.15.
Alias for within_latent, used by ml_var and ml_gvar (matching the mlVAR package).
Sacha Epskamp <mail@sachaepskamp.com>
For each subject the observed \(p\)-variate series is augmented into complete lag
pairs \(z_{it} = (y_{i,t-1}', y_{it}')'\) (the \(2p\)-variate "Toeplitz" embedding
also used by var1; pairs are broken across days/beeps via dayvar /
beepvar). These lag pairs are then treated, within subject, as
conditionally independent given the subject's random intercept \(b_i\):
$$z_{it} \mid b_i \sim N\!\left(\mu_z + (b_i', b_i')',\ \Sigma_W\right), \qquad b_i \sim N(0, \Sigma_B),$$
with \(\mu_z = (\mu', \mu')'\) the stationary mean repeated across
the two halves, the within-person block a stationary Toeplitz structure
\(\Sigma_W = [[\Sigma_0, \Sigma_0 \beta'], [\beta \Sigma_0, \Sigma_0]]\) with
\(\mathrm{vec}(\Sigma_0) = (I - \beta \otimes \beta)^{-1} \mathrm{vec}(\Sigma_\zeta)\),
and the between-person block \(\Sigma_{Bf} = 1_{2 \times 2} \otimes \Sigma_B\) (a rank-\(p\), \(2p \times 2p\) matrix). Minus twice the
log-likelihood is the two-level Gaussian objective of McDonald & Goldstein (1989) /
Muthen (1990) evaluated at \((\mu_z, \Sigma_W, \Sigma_{Bf})\) with the subject as cluster.
Because consecutive lag pairs share \(y_{it}\), this is a working-independence
pseudo (composite) likelihood, exactly the move var1 makes at a single
subject. Its key advantage is that the data are compressed once into the two-level
sufficient statistics of the lag pairs, so every fit / gradient / Fisher-information
evaluation works only with \(2p \times 2p\) matrices, and the per-evaluation cost is
independent of the number of time points. ml_var1 therefore uses the same
data augmentation as var1 and the same two-level estimator as
ml_lvm (estimator = "ML").
ml_var1 shares the estimand with panelvar / panelgvar (a
multi-level VAR(1) with random intercepts) but not the estimator: panelvar forms
the full \((T p) \times (T p)\) block-Toeplitz likelihood, which is statistically ideal
for panel data (few waves, many subjects) but grows with the number of waves \(T\) and
becomes infeasible for intensive longitudinal / ESM data with dozens to hundreds of
measurements per subject. Use panelvar / panelgvar for few-wave panel
data, and ml_var1 / ml_gvar1 for many repeated measures where panelvar
is infeasible. In this first version, incomplete lag pairs (rows with any missing value,
including the structural lag rows at day/beep breaks) are dropped (complete-pairs, listwise
analysis); a full-information treatment of the augmented rows is planned for a later version.
Statistical caveats. Because the estimator treats consecutive lag pairs as
independent (working-independence pseudo-likelihood, as in var1), two caveats
apply. First, under positive temporal autocorrelation the between-person covariance
\(\Sigma_B\) carries a small upward finite-length bias of order \(O(1/T)\): in
simulation (\(p = 4\), \(J = 60\), \(\beta\) spectral radius 0.5) the
between-person variances were biased about +9% at \(T = 25\) and +3% at \(T = 50\),
and negligibly by \(T = 100\); a full-likelihood panelvar / panelgvar
fit does not show this bias. Second, the Fisher-based standard errors are pseudo-likelihood
standard errors and can be anti-conservative -- most noticeably for the contemporaneous
(within) network, where reported standard errors were about 75--80% of the empirical
standard deviation, an under-statement that grows with temporal dependence and does not
vanish as \(T\) grows. Point estimates of the temporal network beta and the
contemporaneous network omega_zeta_within are consistent and closely track the
full-likelihood panelvar estimates.
The saturated (unrestricted) and baseline (independence) reference models used by
runmodel for the fit measures are estimated within the ml_var1
framework itself (a free two-level model over the augmented \(2p\) variables and a
two-level independence model, respectively). The number of independent observations used in
sample-size-dependent fit measures (BIC, RMSEA) is the number of subjects. The model
layer is currently implemented in R only (no C++ acceleration), but estimation is still fast
because every evaluation touches only \(2p \times 2p\) matrices.
McDonald, R. P., & Goldstein, H. (1989). Balanced versus unbalanced designs for linear structural relations in two-level data. British Journal of Mathematical and Statistical Psychology, 42(2), 215-232.
Muthen, B. O. (1990). Mean and covariance structure analysis of hierarchical data. UCLA Statistics Series #62.
Epskamp, S. (2020). Psychometric network models from time-series and panel data. Psychometrika, 85(1), 206-231. tools:::Rd_expr_doi("10.1007/s11336-020-09697-3")
var1 for the single-subject (or pooled) VAR(1)/GVAR, panelvar / panelgvar for the full-likelihood random-intercept VAR on panel data, ml_lvm for the two-level latent variable model that supplies the estimator, and dlvm1 for the dynamic latent variable model.
if (FALSE) {
library("dplyr")
# Simulate a tiny two-level VAR(1) data set (random intercepts):
set.seed(1)
p <- 3 # nodes
J <- 40 # subjects
TT <- 25 # measurements per subject
# Temporal (beta) and contemporaneous (within) structure:
beta <- diag(0.3, p); beta[1, 2] <- 0.2; beta[3, 1] <- -0.15
SigZ <- diag(p) * 0.8; SigZ[1, 2] <- SigZ[2, 1] <- 0.2
Sig0 <- matrix(solve(diag(p^2) - kronecker(beta, beta), as.vector(SigZ)), p, p)
SigB <- diag(p) * 0.4; SigB[2, 3] <- SigB[3, 2] <- 0.1
# Generate:
long <- do.call(rbind, lapply(seq_len(J), function(i){
bi <- MASS::mvrnorm(1, rep(0, p), SigB) # random intercept
eta <- MASS::mvrnorm(1, rep(0, p), Sig0)
Y <- matrix(NA_real_, TT, p)
Y[1, ] <- eta + bi
for (t in 2:TT){
eta <- as.vector(beta %*% eta) + MASS::mvrnorm(1, rep(0, p), SigZ)
Y[t, ] <- eta + bi
}
data.frame(id = i, beep = seq_len(TT), Y)
}))
names(long)[3:5] <- paste0("V", 1:p)
# Fit the multi-level graphical VAR (GGM within and between):
mod <- ml_gvar1(long, vars = paste0("V", 1:p), idvar = "id", beepvar = "beep")
mod <- runmodel(mod)
# The three multi-level GVAR networks:
getmatrix(mod, "beta") # temporal network
getmatrix(mod, "omega_zeta_within") # contemporaneous network
getmatrix(mod, "omega_zeta_between") # between-person network
}
Run the code above in your browser using DataLab