Multivariate Meta-Analysis of Discrete-Time VAR Estimates (Fixed-Effect Model)

Ivan Jacob Agaloos Pesigan

2025-05-08

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} + \boldsymbol{\beta} \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}$, $\boldsymbol{\beta}$, 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}$ denotes a vector of intercepts, $\boldsymbol{\beta}$ a matrix of autoregression and cross regression coefficients, 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 transition matrix $\boldsymbol{\beta}$ be given by

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

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] 3.814355 2.576348
sigma0
#>           [,1]      [,2]
#> [1,] 1.3978842 0.5782369
#> [2,] 0.5782369 1.6636513
sigma0_l # sigma0_l <- t(chol(sigma0))
#>           [,1]     [,2]
#> [1,] 1.1823215 0.000000
#> [2,] 0.4890691 1.193509
alpha
#> [1] 2.87 2.04
beta
#>        [,1]   [,2]
#> [1,]  0.280 -0.048
#> [2,] -0.035  0.260
psi
#>           [,1]      [,2]
#> [1,] 1.3000000 0.5696315
#> [2,] 0.5696315 1.5600000
psi_l # psi_l <- t(chol(psi))
#>           [,1]     [,2]
#> [1,] 1.1401754 0.000000
#> [2,] 0.4995998 1.144727
nu
#> [1] 0 0
lambda
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
theta
#>      [,1] [,2]
#> [1,]  0.2  0.0
#> [2,]  0.0  0.2
theta_l # theta_l <- t(chol(theta))
#>           [,1]      [,2]
#> [1,] 0.4472136 0.0000000
#> [2,] 0.0000000 0.4472136

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

Using the SimSSMFixed Function from the simStateSpace Package to Simulate Data

library(simStateSpace)
sim <- SimSSMFixed(
  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.455114 1.9223092
#> 2  1    1 5.696312 0.6452660
#> 3  1    2 4.364664 1.4954592
#> 4  1    3 4.900281 2.6730811
#> 5  1    4 2.325515 0.2139305
#> 6  1    5 3.340405 1.9930331
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,
  sigma0_values = sigma0,
  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 and alpha.

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

#> 
#>  Solution found!  Final fit=14630.543 (started at 25387129)  (1 attempt(s): 1 valid, 0 errors)
#>  Start values from best fit:
#> 0.185482523593782,0.0121788228619116,0.0115184897567764,0.184817876404306,3.10216459456662,2.02115148863649
summary(meta)
#>         est     se        z p   2.5%  97.5%
#> b0_1 0.1855 0.0017 109.4794 0 0.1822 0.1888
#> b0_2 0.0122 0.0014   8.9432 0 0.0095 0.0148
#> b0_3 0.0115 0.0013   8.9464 0 0.0090 0.0140
#> b0_4 0.1848 0.0017 108.4805 0 0.1815 0.1882
#> b0_5 3.1022 0.0064 483.2179 0 3.0896 3.1147
#> b0_6 2.0212 0.0060 334.1024 0 2.0093 2.0330

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/