Multivariate Meta-Analysis of Discrete-Time VAR Estimates (Random-Effects Model)

Ivan Jacob Agaloos Pesigan

2025-08-18

Note: When the discrete-time vector autoregressive model is estimated separately for each ID, any of the parameters can vary across individuals. In this example, however, we are only varying $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$.

Model

The measurement model is given by $$\begin{equation} \mathbf{y}_{i, t} = \boldsymbol{\Lambda} \boldsymbol{\eta}_{i, t} + \boldsymbol{\varepsilon}_{i, t}, \quad \mathrm{with} \quad \boldsymbol{\varepsilon}_{i, t} \sim \mathcal{N} \left( \mathbf{0}, \boldsymbol{\Theta} \right) \end{equation}$$ where $\mathbf{y}_{i, t}$, $\boldsymbol{\eta}_{i, t}$, and $\boldsymbol{\varepsilon}_{i, t}$ are random variables and $\boldsymbol{\Lambda}$, and $\boldsymbol{\Theta}$ are model parameters. $\mathbf{y}_{i, t}$ represents a vector of observed random variables, $\boldsymbol{\eta}_{i, t}$ a vector of latent random variables, and $\boldsymbol{\varepsilon}_{i, t}$ a vector of random measurement errors, at time $t$ and individual $i$. $\boldsymbol{\Lambda}$ denotes a matrix of factor loadings, and $\boldsymbol{\Theta}$ the covariance matrix of $\boldsymbol{\varepsilon}$. In this model, $\boldsymbol{\Lambda}$ is an identity matrix and $\boldsymbol{\Theta}$ is a diagonal matrix.

The dynamic structure is given by $$\begin{equation} \boldsymbol{\eta}_{i, t} = \boldsymbol{\alpha}_{i} + \boldsymbol{\beta}_{i} \boldsymbol{\eta}_{i, t - 1} + \boldsymbol{\zeta}_{i, t}, \quad \mathrm{with} \quad \boldsymbol{\zeta}_{i, t} \sim \mathcal{N} \left( \mathbf{0}, \boldsymbol{\Psi} \right) \end{equation}$$ where $\boldsymbol{\eta}_{i, t}$, $\boldsymbol{\eta}_{i, t - 1}$, and $\boldsymbol{\zeta}_{i, t}$ are random variables, and $\boldsymbol{\alpha}_{i}$, $\boldsymbol{\beta}_{i}$, and $\boldsymbol{\Psi}$ are model parameters. Here, $\boldsymbol{\eta}_{i, t}$ is a vector of latent variables at time $t$ and individual $i$, $\boldsymbol{\eta}_{i, t - 1}$ represents a vector of latent variables at time $t - 1$ and individual $i$, and $\boldsymbol{\zeta}_{i, t}$ represents a vector of dynamic noise at time $t$ and individual $i$. $\boldsymbol{\alpha}_{i}$ denotes a vector of intercepts for individual $i$, $\boldsymbol{\beta}_{i}$ a matrix of autoregression and cross regression coefficients for individual $i$, and $\boldsymbol{\Psi}$ the covariance matrix of $\boldsymbol{\zeta}_{i, t}$. In this model, $\boldsymbol{\Psi}$ is a symmetric matrix.

Data Generation

The parameters used in this example were based on Bringmann et al. (2013).

Notation

Let $t = 1000$ be the number of time points and $n = 1000$ be the number of individuals.

Let the measurement model intecept vector $\boldsymbol{\nu}$ be given by

$$\begin{equation} \boldsymbol{\nu} = \left( \begin{array}{c} 0 \\ 0 \\ \end{array} \right) . \end{equation}$$

Let the factor loadings matrix $\boldsymbol{\Lambda}$ be given by

$$\begin{equation} \boldsymbol{\Lambda} = \left( \begin{array}{cc} 1 & 0 \\ 0 & 1 \\ \end{array} \right) . \end{equation}$$

Let the measurement error covariance matrix $\boldsymbol{\Theta}$ be given by

$$\begin{equation} \boldsymbol{\Theta} = \left( \begin{array}{cc} 0.2 & 0 \\ 0 & 0.2 \\ \end{array} \right) . \end{equation}$$

Let the initial condition $\boldsymbol{\eta}_{0}$ be given by

$$\begin{equation} \boldsymbol{\eta}_{0} \sim \mathcal{N} \left( \boldsymbol{\mu}_{\boldsymbol{\eta} \mid 0}, \boldsymbol{\Sigma}_{\boldsymbol{\eta} \mid 0} \right) \end{equation}$$

$$\begin{equation} \boldsymbol{\mu}_{\boldsymbol{\eta} \mid 0} = \left( \begin{array}{c} 3.8143546 \\ 2.5763481 \\ \end{array} \right) \end{equation}$$

$$\begin{equation} \boldsymbol{\Sigma}_{\boldsymbol{\eta} \mid 0} = \left( \begin{array}{cc} 1.3978842 & 0.5782369 \\ 0.5782369 & 1.6636513 \\ \end{array} \right) . \end{equation}$$

Let the intercept vector $\boldsymbol{\alpha}$ be normally distributed with the following means

$$\begin{equation} \left( \begin{array}{c} 2.87 \\ 2.04 \\ \end{array} \right) \end{equation}$$

and covariance matrix

$$\begin{equation} \left( \begin{array}{cc} 1.2 & 0.459565 \\ 0.459565 & 1.1 \\ \end{array} \right) . \end{equation}$$

Let the transition matrix $\boldsymbol{\beta}$ be normally distributed with the following means

$$\begin{equation} \left( \begin{array}{cc} 0.28 & -0.048 \\ -0.035 & 0.26 \\ \end{array} \right) \end{equation}$$

and covariance matrix

$$\begin{equation} \left( \begin{array}{cccc} 0.0169 & 0.00468 & 0.001456 & 0.00832 \\ 0.00468 & 0.0081 & 0.001008 & 0.00576 \\ 0.001456 & 0.001008 & 0.000784 & 0.001792 \\ 0.00832 & 0.00576 & 0.001792 & 0.0256 \\ \end{array} \right) . \end{equation}$$

The SimAlphaN and SimBetaN functions from the simStateSpace package generates random intercept vectors and transition matrices from the multivariate normal distribution. Note that the SimBetaN function generates transition matrices that are weakly stationary.

Let the dynamic process noise $\boldsymbol{\Psi}$ be given by

$$\begin{equation} \boldsymbol{\Psi} = \left( \begin{array}{cc} 1.3 & 0.5696315 \\ 0.5696315 & 1.56 \\ \end{array} \right) . \end{equation}$$

R Function Arguments

n
#> [1] 1000
time
#> [1] 1000
mu0
#> [[1]]
#> [1] 3.814355 2.576348
sigma0
#> [[1]]
#>           [,1]      [,2]
#> [1,] 1.3978842 0.5782369
#> [2,] 0.5782369 1.6636513
sigma0_l # sigma0_l <- t(chol(sigma0))
#> [[1]]
#>           [,1]     [,2]
#> [1,] 1.1823215 0.000000
#> [2,] 0.4890691 1.193509
# first alpha in the list of length n
alpha[[1]]
#> [1] 1.382816 2.056527
# first beta in the list of length n
beta[[1]]
#>            [,1]       [,2]
#> [1,] 0.21090478 0.00660174
#> [2,] 0.01360809 0.37057228
psi
#>           [,1]      [,2]
#> [1,] 1.3000000 0.5696315
#> [2,] 0.5696315 1.5600000
psi_l # psi_l <- t(chol(psi))
#> [[1]]
#>           [,1]     [,2]
#> [1,] 1.1401754 0.000000
#> [2,] 0.4995998 1.144727
nu
#> [[1]]
#> [1] 0 0
lambda
#> [[1]]
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
theta
#>      [,1] [,2]
#> [1,]  0.2  0.0
#> [2,]  0.0  0.2

Visualizing the Dynamics Without Process Noise (n = 5 with Different Initial Condition)

Using the SimSSMIVary Function from the simStateSpace Package to Simulate Data

library(simStateSpace)
sim <- SimSSMIVary(
  n = n,
  time = time,
  mu0 = mu0,
  sigma0_l = sigma0_l,
  alpha = alpha,
  beta = beta,
  psi_l = psi_l,
  nu = nu,
  lambda = lambda,
  theta_l = theta_l
)
data <- as.data.frame(sim)
head(data)
#>   id time       y1       y2
#> 1  1    0 3.040533 3.202203
#> 2  1    1 3.630653 4.691667
#> 3  1    2 2.621350 5.240797
#> 4  1    3 2.073843 5.555230
#> 5  1    4 2.013430 5.269451
#> 6  1    5 1.677296 2.836631
plot(sim)

Model Fitting

The FitDTVARIDMx function fits a DT-VAR model on each individual $i$.

library(fitDTVARMx)
fit <- FitDTVARIDMx(
  data = data,
  observed = paste0("y", seq_len(k)),
  id = "id",
  mu0_values = mu0[[1]],
  sigma0_values = sigma0[[1]],
  ncores = parallel::detectCores()
)

Multivariate Meta-Analysis

The MetaVARMx function performs multivariate meta-analysis using the estimated parameters and the corresponding sampling variance-covariance matrix for each individual $i$. Estimates with the prefix b0 correspond to the estimates of beta_mu and alpha_mu. Estimates with the prefix t2 correspond to the estimates of beta_sigma and alpha_sigma. Estimates with the prefix i2 correspond to the estimates of heterogeniety.

library(metaVAR)
meta <- MetaVARMx(
  object = fit,
  intercept = TRUE,
  ncores = parallel::detectCores()
)
#> Running Model with 27 parameters
#> 
#> Beginning initial fit attempt
#> Running Model with 27 parameters
#> 
#>  Lowest minimum so far:  10219.2546092858
#> 
#> Solution found

#> 
#>  Solution found!  Final fit=10219.255 (started at 36755.244)  (1 attempt(s): 1 valid, 0 errors)
#>  Start values from best fit:
#> 0.172176323087149,0.0872331149200129,0.0506409957838052,0.195435618555425,2.98880430042743,1.73526790798659,0.556381434337175,-0.353109510117613,0.100488432796781,-0.416966768192526,-2.62888604739992,2.37372360926509,0.356751790512116,-0.0590609265736989,0.0165839459274479,0.279295888611737,-1.36502655242848,0.201233008455563,-0.232579211238073,-0.52462463618858,0.859705901533594,0.248594903277555,0.247482332249011,-0.232085883462846,1.33824250830376,0.0819217334694941,1.16687501844385
summary(meta)
#>            est     se         z p    2.5%   97.5%
#> b0_1    0.1722 0.0180    9.5893 0  0.1370  0.2074
#> b0_2    0.0872 0.0164    5.3256 0  0.0551  0.1193
#> b0_3    0.0506 0.0079    6.3895 0  0.0351  0.0662
#> b0_4    0.1954 0.0174   11.2023 0  0.1612  0.2296
#> b0_5    2.9888 0.0965   30.9798 0  2.7997  3.1779
#> b0_6    1.7353 0.0996   17.4213 0  1.5400  1.9305
#> t2_1_1  0.3096 0.0162   19.1483 0  0.2779  0.3412
#> t2_2_1 -0.1965 0.0123  -15.9209 0 -0.2206 -0.1723
#> t2_3_1  0.0559 0.0072    7.7291 0  0.0417  0.0701
#> t2_4_1 -0.2320 0.0141  -16.4597 0 -0.2596 -0.2044
#> t2_5_1 -1.4627 0.0790  -18.5092 0 -1.6175 -1.3078
#> t2_6_1  1.3207 0.0765   17.2695 0  1.1708  1.4706
#> t2_2_2  0.2520 0.0156   16.1817 0  0.2214  0.2825
#> t2_3_2 -0.0566 0.0051  -11.0389 0 -0.0666 -0.0465
#> t2_4_2  0.1532 0.0112   13.6651 0  0.1312  0.1751
#> t2_5_2  1.0279 0.0645   15.9463 0  0.9016  1.1543
#> t2_6_2 -1.3252 0.0807  -16.4169 0 -1.4834 -1.1670
#> t2_3_3  0.0541 0.0044   12.4057 0  0.0455  0.0626
#> t2_4_3 -0.0897 0.0076  -11.7977 0 -0.1046 -0.0748
#> t2_5_3 -0.3862 0.0417   -9.2723 0 -0.4679 -0.3046
#> t2_6_3  0.4922 0.0378   13.0211 0  0.4181  0.5662
#> t2_4_4  0.2900 0.0158   18.4113 0  0.2592  0.3209
#> t2_5_4  1.2843 0.0752   17.0737 0  1.1369  1.4318
#> t2_6_4 -1.2700 0.0752  -16.8949 0 -1.4174 -1.1227
#> t2_5_5  9.1164 0.4470   20.3967 0  8.2404  9.9924
#> t2_6_5 -7.0203 0.4048  -17.3444 0 -7.8136 -6.2270
#> t2_6_6  9.6591 0.4851   19.9116 0  8.7083 10.6099
#> i2_1    0.9751 0.0013  770.5674 0  0.9727  0.9776
#> i2_2    0.9802 0.0012  817.2680 0  0.9779  0.9826
#> i2_3    0.9392 0.0046  204.0035 0  0.9302  0.9482
#> i2_4    0.9746 0.0013  724.4164 0  0.9719  0.9772
#> i2_5    0.9956 0.0002 4634.8032 0  0.9952  0.9960
#> i2_6    0.9966 0.0002 5872.6729 0  0.9963  0.9969

References

Bringmann, L. F., Vissers, N., Wichers, M., Geschwind, N., Kuppens, P., Peeters, F., Borsboom, D., & Tuerlinckx, F. (2013). A network approach to psychopathology: New insights into clinical longitudinal data. PLoS ONE, 8(4). https://doi.org/10.1371/journal.pone.0060188

Cheung, M. W.-L. (2015). Meta‐analysis: A structural equation modeling approach. Wiley. https://doi.org/10.1002/9781118957813

Hunter, M. D. (2017). State space modeling in an open source, modular, structural equation modeling environment. Structural Equation Modeling: A Multidisciplinary Journal, 25(2), 307–324. https://doi.org/10.1080/10705511.2017.1369354

Neale, M. C., Hunter, M. D., Pritikin, J. N., Zahery, M., Brick, T. R., Kirkpatrick, R. M., Estabrook, R., Bates, T. C., Maes, H. H., & Boker, S. M. (2015). OpenMx 2.0: Extended structural equation and statistical modeling. Psychometrika, 81(2), 535–549. https://doi.org/10.1007/s11336-014-9435-8

R Core Team. (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/