Fitting Multilevel Vector Autoregressive Models in Stan, JAGS, and Mplus

ABSTRACT

The influx of intensive longitudinal data creates a pressing need for complex modeling tools that help enrich our understanding of how individuals change over time. Multilevel vector autoregressive (mlVAR) models allow for simultaneous evaluations of reciprocal linkages between dynamic processes and individual differences, and have gained increased recognition in recent years. High-dimensional and other complex variations of mlVAR models, though often computationally intractable in the frequentist framework, can be readily handled using Markov chain Monte Carlo techniques in a Bayesian framework. However, researchers in social science fields may be unfamiliar with ways to capitalize on recent developments in Bayesian software programs. In this paper, we provide step-by-step illustrations and comparisons of options to fit Bayesian mlVAR models using Stan, JAGS and Mplus, supplemented with a Monte Carlo simulation study. An empirical example is used to demonstrate the utility of mlVAR models in studying intra- and inter-individual variations in affective dynamics.

KEYWORDS:

• Multilevel vector autoregressive models
• Bayesian modeling
• missing data
• affective dynamics

Acknowledgments

We would like to thank Dr. Linhai Song for permission to use his lab server as additional computational resources.

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website.

Correction Statement

This article has been republished with minor changes. These changes do not impact the academic content of the article.

Notes

¹ Since the data set was used only for illustration purposes, we pre-processed the raw data following the same procedure as implemented in previous studies. Other data pre-processing and exclusion criteria may be used in other studies as appropriate.

² Extraversion was measured using three selected questions from the revised NEO personality inventory (NEO-PI-R) on which the participants were asked to rate on a 4-point scale the extent to which they felt Passive vs. Active; Unenergetic vs. Energetic; and Dominant vs. Submissive since the last time they responded to the survey.

³ Although the presented covariance matrix decomposition and the subsequent Fisher-z transformation of the correlation do not generalize to more than 2 dimensions, it was adopted in this specific context to aid computational efficiency, and to provide simple, intuitive interpretations of the effects of predictors on the person-specific correlations (see also Oravecz et al., 2011).

⁴ Note that this treatment of missing data is limited to two dimensions and do not readily generalize to higher dimensions. Such programming was implemented only to work around the restriction of handling partially observed multivariate nodes in JAGS, and is not recommended in principle.

⁵ With the R package, MplusAutomation, users can extract all posterior samples and include part of the first half for estimation, but caution should be taken when deciding the number of burn-in iterations to make sure that the chains are not influenced by the starting values anymore.

⁶ The estimated population mean of person-specific baselines of PA was calculated based on ${\hat{\mu }}_1+{\hat{\alpha }}_{{\mu }_1}{\bar{x}}_1+{\hat{\beta }}_{{\mu }_1}{\bar{x}}_2$ where ${\bar{x}}_1$ and ${\bar{x}}_2$ were empirical means of $x_{1,i}$ and $x_{2,i}$ across persons. Since $x_{1,i}$ and $x_{2,i}$ have been standardized to zero mean and unit variance, the estimated population mean was equal to ${\hat{\mu }}_1$ (i.e., the estimate of ${\mu }_1$ in Table 1). The estimated population means of other person-specific parameters were calculated in a similar way.

⁷ Ideally, more Monte Carlo replications need to be run to generate more reliable results. In this study, we chose 100 replications due to the high computational time. To verify that 100 replications were sufficient for the simulation purpose, we ran 500 replications with JAGS under the N = 100, T = 60, and high-stability condition, and the Monte Carlo standard error (MCSE; see Equation (15)(15) $\mathrm{M}\mathrm{C}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{M-1}\sum _{m=1}^M{({\widehat{\text{thetaa}}}_m-\overline{\text{thetaa}})}^2},$(15) ) did not change substantially compared with that with 100 replications. We note that some previous studies on Bayesian VAR models also used 100 replications (see, e.g., Huber & Feldkircher, 2019; Ji et al., 2020).

⁸ By checking the samplers used in JAGS, we found that when all person-specific parameters were specified to be correlated, JAGS would use an adaptive random walk Metropolis algorithm called bugs::MNormal for these correlated parameters. According to the JAGS manual, this algorithm can be very inefficient and may require an extremely long adaptation period, which has been verified by our simulation, where about 20% of the replications with JAGS did not finish the adaptation phase within 4000 iterations. In addition, we found that when only person-specific intercepts were specified to be correlated to yield a 2-dimensional covariance matrix, JAGS would use a Gibbs sampling algorithm called bugs::ConjugateMNormal, which was more efficient, as indicated by much higher ESSs.

Additional information

Funding

Research reported in this publication was supported by the Intensive Longitudinal Health Behavior Cooperative Agreement Program funded by the National Institutes of Health [U24AA027684]; National Science Foundation [IGE-1806874]. Part of the computations for this research were performed on the Pennsylvania State University’s Institute for Computational and Data Sciences’ Roar supercomputer.