Var1 Family

Lag-1 vector autoregressive models for time-series data

Contents

Overview

The var1 family models time-series data using a lag-1 vector autoregressive (VAR) model. This framework decomposes multivariate time-series into two network structures:

The model is fitted to the Toeplitz (block) covariance matrix of the augmented data, where both the lagged ($t-1$) and current ($t$) measurements are included. Importantly, the block of covariances of the lagged variables is modeled completely separately (using a Cholesky decomposition) and does not enter the temporal or contemporaneous model estimation.

The gvar() wrapper function sets contemporaneous = "ggm", resulting in a graphical VAR model where the contemporaneous structure is represented as a Gaussian Graphical Model (network of partial correlations among innovations).

The VAR(1) Model

The VAR(1) model specifies that each variable at time $t$ is a linear function of all variables at time $t-1$, plus an innovation term:

$$\boldsymbol{y}_t = \boldsymbol{\mu} + \boldsymbol{B}(\boldsymbol{y}_{t-1} - \boldsymbol{\mu}) + \boldsymbol{\zeta}_t$$

Where:

Under stationarity, define $\boldsymbol{\Sigma}_0 = \text{var}(\boldsymbol{y}_t)$ as the $p \times p$ stationary (lag-0) covariance matrix and $\boldsymbol{\Sigma}_1 = \text{cov}(\boldsymbol{y}_t,\, \boldsymbol{y}_{t-1})$ as the $p \times p$ lag-1 cross-covariance matrix. The model then implies:

The model is fitted by maximum likelihood estimation on the combined Toeplitz block covariance matrix formed from $\boldsymbol{\Sigma}_0$ and $\boldsymbol{\Sigma}_1$.

Contemporaneous Parameterizations

The contemporaneous covariance matrix $\boldsymbol{\Sigma}_\zeta$ (the covariance among innovations) can be parameterized in four ways using the contemporaneous argument:

Type Key Matrix Equation
"cov" (default) sigma_zeta $\boldsymbol{\Sigma}_\zeta$ (directly estimated)
"chol" lowertri_zeta $\boldsymbol{\Sigma}_\zeta = \boldsymbol{L}_\zeta\boldsymbol{L}_\zeta'$
"prec" kappa_zeta $\boldsymbol{\Sigma}_\zeta = \boldsymbol{K}_\zeta^{-1}$
"ggm" omega_zeta, delta_zeta $\boldsymbol{\Sigma}_\zeta = \boldsymbol{\Delta}_\zeta(\boldsymbol{I} - \boldsymbol{\Omega}_\zeta)^{-1}\boldsymbol{\Delta}_\zeta$
Note: The gvar() wrapper is shorthand for var1(..., contemporaneous = "ggm"). This is the recommended parameterization for network analysis, as it represents the contemporaneous structure as a network of partial correlations among innovations.

Data Preparation

Time-series data typically comes from experience sampling methodology (ESM) or ecological momentary assessment (EMA) studies, where participants are measured repeatedly over time. The var1() function handles the data augmentation (creating lagged variables) internally.

Key Arguments for Time-Series Data

Example Data Structure

# Typical ESM data format:
#   id  day  beep  relaxed  sad  nervous  concentration
#   1   1    1     3        2    1        4
#   1   1    2     4        1    2        3
#   1   1    3     2        3    1        4
#   1   2    1     3        2    2        3   # New day: not lagged from previous row
#   ...

# Fit model with day structure:
model <- gvar(data, vars = c("relaxed", "sad", "nervous", "concentration"),
              dayvar = "day", beepvar = "beep")

Example: Graphical VAR on Simulated Data

This example demonstrates the basic workflow using simulated data from the graphicalVAR package. We simulate a 2-variable VAR(1) process with known parameters and recover them using gvar(). 

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

# Simulate a 2-variable VAR(1) process:
beta <- matrix(c(
  0,   0.5,
  0.5, 0
), 2, 2, byrow = TRUE)
kappa <- diag(2)
simData <- graphicalVARsim(50, beta, kappa)

# Fit graphical VAR model:
model <- gvar(simData)
model <- model %>% runmodel

# Parameter estimates:
model %>% parameters

# Confidence interval plot for temporal effects:
CIplot(model, "beta")

# Extract the temporal network:
getmatrix(model, "beta")

# Extract the contemporaneous correlations/network:
getmatrix(model, "omega_zeta")

Example: Clinical Time-Series (ESM Data)

This example uses ESM data from a single clinical patient (Epskamp, van Borkulo, van der Veen, Servaas, Isvoranu, Riese, & Cramer, 2018), available at https://osf.io/c8wjz/. The dataset contains 7 variables measured multiple times per day. 

# Download the data from https://osf.io/c8wjz/
tsdata <- read.csv("Supplementary2_data.csv")

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

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

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

# Fit graphical VAR with FIML for missing data:
model <- gvar(tsdata, vars = vars, dayvar = "Day",
              estimator = "FIML")
model <- model %>% runmodel

# Model search:
model_pruned <- model %>% prune(alpha = 0.01)
model_final <- model_pruned %>% stepup(criterion = "bic")

# Compare models:
compare(saturated = model, pruned = model_pruned, final = model_final)

# Extract networks:
temporal <- getmatrix(model_final, "PDC")
contemporaneous <- getmatrix(model_final, "omega_zeta")

# Visualize both networks side by side:
library("qgraph")
labs <- gsub("\\.", "\n", vars)
L <- averageLayout(temporal, contemporaneous)

layout(t(1:2))
qgraph(temporal, layout = L, theme = "colorblind", cut = 0,
       directed = TRUE, diag = TRUE,
       title = "Temporal", vsize = 12, mar = rep(6, 4),
       asize = 5, labels = labs)
qgraph(contemporaneous, layout = L, theme = "colorblind", cut = 0,
       title = "Contemporaneous", vsize = 12, mar = rep(6, 4),
       asize = 5, labels = labs)
Note: The dayvar = "Day" argument ensures the first measurement of each day is not regressed on the last measurement of the previous day, properly handling overnight gaps in ESM data.
Warning: Do not use dayvar when there is only one measurement per day. Because the first measurement of each day is excluded from the lagged regression, using dayvar with daily data will result in all observations being removed.

Derived Network Matrices (PDC & PCC)

In addition to the raw model parameters (beta, omega_zeta), psychonetrics can compute standardized network matrices that are easier to interpret and compare:

PDC: Partial Directed Correlations

The Partial Directed Correlations (PDC) are the standardized version of the temporal regression coefficients. Specifically, the PDC matrix is the transpose of the beta matrix after standardization. This transposition is needed because in the beta matrix, element $b_{ij}$ represents the effect of variable $j$ on variable $i$, whereas qgraph() expects element $(i,j)$ to represent an edge from node $i$ to node $j$. Extract with: 

temporal_network <- getmatrix(model, "PDC")
Warning: Do not use getmatrix(model, "beta") directly in qgraph(). Because of the row/column convention difference described above, the beta matrix will display edges in the wrong direction. Always use getmatrix(model, "PDC") for visualization, which transposes and standardizes the temporal effects correctly.

When visualizing the PDC matrix with qgraph(), use directed = TRUE and diag = TRUE to display the directed edges and auto-regressive effects.

PCC: Partial Contemporaneous Correlations

The Partial Contemporaneous Correlations (PCC) are the standardized precision matrix of the innovations. When contemporaneous = "ggm", the PCC is equivalent to the omega_zeta matrix. Extract with: 

contemporaneous_network <- getmatrix(model, "PCC")

# Equivalent to omega_zeta when using gvar():
contemporaneous_network2 <- getmatrix(model, "omega_zeta")

Summary

The var1 family provides the foundation for time-series network modeling in psychonetrics. Key takeaways:

Next Steps

Now that you understand var1 models, you can explore:

Further Reading