Multi-level VAR Family

Overview Experimental (0.16.6+)

The ml_var family models a multi-level (random intercept) lag-1 vector autoregression for intensive longitudinal data: many repeated measures nested in subjects, as collected in experience sampling (ESM) or ecological momentary assessment (EMA) studies. The ml_gvar() wrapper parameterizes the covariance structures as Gaussian graphical models, yielding the three familiar multi-level GVAR networks (as in the mlVAR package):

• Temporal network: directed lag-1 relationships within persons (the beta matrix; extract the standardized network with getmatrix(model, "PDC")).
• Contemporaneous network: partial correlations among the within-person innovations (omega_zeta_within).
• Between-persons network: partial correlations among the person-specific means (omega_zeta_between).

The estimand is the same multi-level VAR with random intercepts that panelvar() / panelgvar() target, but the family is designed for long time series, where the full-information panel likelihood is computationally infeasible. Its default summary-statistics estimator compresses each subject's time series once into two-level sufficient statistics, so estimation cost is independent of the number of time points.

The Multi-level VAR(1) Model

Each subject $i$ follows a within-person stationary VAR(1) process around a person-specific mean:

$$\boldsymbol{y}_{it} = \boldsymbol{\mu} + \boldsymbol{b}_i + \boldsymbol{B}\,(\boldsymbol{y}_{i,t-1} - \boldsymbol{\mu} - \boldsymbol{b}_i) + \boldsymbol{\zeta}_{it}, \qquad \boldsymbol{b}_i \sim N(\boldsymbol{0}, \boldsymbol{\Sigma}_B)$$

• $\boldsymbol{B}$ (beta): the fixed temporal effects (the temporal network after standardization/transposition).
• $\boldsymbol{\Sigma}_\zeta$ (sigma_zeta_within / omega_zeta_within): the within-person innovation structure (contemporaneous network under the GGM parameterization).
• $\boldsymbol{\Sigma}_B$ (sigma_zeta_between / omega_zeta_between): the covariance of the random intercepts (between-persons network under the GGM parameterization).

Both covariance structures accept the usual five parameterizations through the within and between arguments ("cov", "chol", "prec", "ggm", "cor"); ml_gvar() sets both to "ggm". Note that this model estimates fixed (population) networks with random intercepts; unlike mlVAR(), it does not model random effects in the temporal or contemporaneous structures themselves.

Estimators: "auto", "ML" and "FIML"

Two estimators are available, selected with the estimator argument:

The default estimator = "auto" chooses "FIML" when the measurement-occasion grid has at most 10 occasions per subject and "ML" otherwise, and always prints a message stating the choice. Set estimator explicitly to override.

Data Preparation

ml_var() takes long-format raw data (one row per subject per measurement occasion), exactly like mlVAR():

• idvar (required): the subject ID column. Each subject's rows form one within-person time series.
• vars: character vector of variable names.
• dayvar: assessment-day column; ensures the first measurement of a day is not regressed on the last measurement of the previous day.
• beepvar: beep-number column; non-consecutive beeps are not linked by a lag-1 regression.
• groupvar: optional grouping variable for multi-group analysis (with equal constraints across groups).

Under estimator = "ML", incomplete lag pairs (rows with any missing value among the current and lagged variables) are currently dropped with a warning (complete-pairs analysis); under estimator = "FIML", missing data are handled by full-information ML.

Example: Multi-level Graphical VAR on ESM Data

library("psychonetrics")
library("dplyr")

# Long-format ESM data:
#   id  day  beep  Stressed  Relaxed  Sad  ...
#   1   1    1     2         4        1
#   1   1    2     3         3        1
#   ...

vars <- c("Stressed", "Relaxed", "Sad", "Nervous", "Concentration")

# Multi-level graphical VAR (the message will report the selected estimator):
model <- ml_gvar(data, vars = vars, idvar = "id",
                 dayvar = "day", beepvar = "beep")
model <- model %>% runmodel

# Fit measures (AIC / BIC etc.):
model %>% summary
model %>% fit

# The three networks:
temporal        <- getmatrix(model, "PDC")                 # standardized temporal
contemporaneous <- getmatrix(model, "omega_zeta_within")   # contemporaneous
between         <- getmatrix(model, "omega_zeta_between")  # between-persons

# Visualize:
library("qgraph")
L <- averageLayout(temporal, contemporaneous, between)
layout(t(1:3))
qgraph(temporal, layout = L, theme = "colorblind", directed = TRUE,
       title = "Temporal", vsize = 12)
qgraph(contemporaneous, layout = L, theme = "colorblind",
       title = "Contemporaneous", vsize = 12)
qgraph(between, layout = L, theme = "colorblind",
       title = "Between-persons", vsize = 12)

# The usual workflow applies:
model_pruned <- model %>% prune(alpha = 0.01)
compare(full = model, pruned = model_pruned)

Statistical Caveats

The summary-statistics estimator ("ML") is a working-independence pseudo (composite) likelihood. Two caveats follow (see ?ml_var1 for simulation details):

• Between-person variances carry a small upward finite-length bias of order $O(1/T)$: around +9% at $T = 25$ and +3% at $T = 50$ measurements per subject, negligible by $T = 100$. The full-likelihood panelvar() route does not show this bias.
• Standard errors are pseudo-likelihood standard errors and can be anti-conservative, most noticeably for the contemporaneous network (reported SEs around 75–80% of the empirical sampling variation).
• Information criteria (AIC/BIC) are pseudo-likelihood criteria: fine for comparing ml_var models with each other, but not against models fit with other estimators.

Point estimates of the temporal and contemporaneous networks are consistent and closely track both the full-likelihood panelvar() estimates and mlVAR() estimates on the same data.

Naming & History

• ml_var() / ml_gvar() are the primary names; 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 (with identical estimates, through the faster panelvar framework) for short occasion grids, and switches to the summary-statistics estimator for longer series where full-information ML was never feasible.
• The old argument names contemporaneous and between are still accepted and map to within / between (i.e., within_latent / between_latent).
• For latent-variable multi-level time-series models, see dlvm1() (the deprecated wrappers ml_tsdlvm1() / ml_ts_lvgvar() route there).

Summary

• ml_gvar() estimates the three multi-level GVAR networks — temporal, contemporaneous, and between-persons — from long-format ESM data with confirmatory model tests, model search, and multi-group support
• The default estimator = "auto" picks full-information ML (via panelvar()) for short occasion grids and the fast summary-statistics estimator for long time series
• Fixed networks with random intercepts: for per-subject random network weights, use mlVAR
• For few-wave panel data use panelvar() / panelgvar(); for meta-analytic pooling of VAR networks across studies use meta_var() / meta_gvar() (see the meta-analysis page)

Further Reading

• Epskamp, S. (2020). Psychometric network models from time-series and panel data. Psychometrika, 85(1), 206–231. DOI: 10.1007/s11336-020-09697-3
• 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.
• Bringmann, L. F., et al. (2013). A network approach to psychopathology: New insights into clinical longitudinal data. PLoS ONE, 8(4), e60188.