Structured Estimation of Heterogeneous Time Series

A number of approaches are capable of accounting for heterogeneity in multivariate time series arising from multiple individuals. Here we briefly highlight a few of these approaches. In doing so we make a distinction between two types of heterogeneity at the time series model level: quantitative heterogeneity and qualitative heterogeneity. We say models allow for quantitative heterogeneity when parameters are allowed to vary in magnitude across individuals in the sample (typically according to a distribution). On the other hand, we say models allow for qualitative heterogeneity when the model structure itself is allowed to vary across individuals. For example, individuals might differ in the lag order of their model, the direction and sign of directed relations, or the overall pattern of zero and nonzero parameters. Finally, models may also accommodate both quantitative and qualitative heterogeneity simultaneously.

The most popular approaches for modeling multiple-subject time series in psychological research are designed to accommodate quantitative heterogeneity (e.g. multilevel time series; bringmann2013a; epskamp2018a; lafit2021; li2022). In these approaches parameters in the model are typically treated as being sampled from a distribution, where fixed and random effects are recovered, and subsequently used to fashion individual-level models from between-person variation in parameter estimates. When all individuals in the sample possess homogenous dynamics, empirical Bayes estimates of the dynamic parameters are generally accurate and reliable (liu2017; liu2018). Less work has focused on how these estimates perform in situations where individuals meaningfully vary in their dynamics.

Another group of approaches account for heterogeneity at the subgroup or cluster level (bulteel2016; takano2021; park2022b; park2022a). These approaches typically allow for groups of individuals to differ qualitatively in terms of their underlying time series model structure. Approaches vary on whether parameters governing individual-level models within the same cluster can vary qualitatively. Furthermore, any model that implies a Gaussian distribution can be used to define the components of a mixture model. Recently, hunter2022 extended this idea to state space models where it becomes possible to have state space mixtures with qualitatively different parameterizations across groups, for example.

Another approach, designed to allow for both quantitative and qualitative heterogeneity across individual time series models, is GIMME (gates2012). The GIMME approach is built on the Structural-Vector Autoregressive (S-VAR) model and uses a stepwise model search algorithm based on modification indices (MI; sorbom1989) in conjunction with test statistics developed in the structural equation modeling (SEM) framework. GIMME is designed to recover both group (generalizable) and individual-level (idiosyncratic) models from multivariate time series data and is available in a well-developed R package (gimme). Recently, the GIMME algorithm has been extended to accommodate the standard and hybrid VAR (luo2022) models, in addition to the original S-VAR. Although not relevant to the current work, there are additional capabilities of GIMME worth mentioning, including the ability to model direct and modulatory effects of tasks (duffy2021), exogenous variables (arizmendi2021), data-driven (gates2017unsupervised) and confirmatory (henry2019) subgroups, multiple solutions (beltz2016) ¹¹1In the original GIMME algorithm each individual begins with a null model (all coefficient values are set to zero). Modification indices are then used to determine which parameter, if freed, would most improve model fit. Multiple solutions refers to the instance where multiple parameters would lead to numerically equivalent improvements in model fit. To address this situation (beltz2016) developed an algorithm, termed multiple solutions GIMME, that accounts for these equivalent solutions in terms of final model selection., and latent variables (gates2020latent).

The multi-VAR approach originally described in fisher2022, and extended in the current work, is similar in spirit to GIMME in that both produce group and individual-level models from qualitatively and quantitatively heterogeneous time series data. However, beyond this common objective, the approaches share little in common in their operatlization. For example, GIMME induces sparsity in individual-level models via forward-selection, a stepwise procedure where paths are freed according to modification indices until fit indices suggest a well-fitting model has been achieved.hastie2015 describe a sparse model as having only a small number of nonzero parameters or weights. In the case of GIMME and multi-VAR we are typically concerned with inducing sparsity in the dynamic part of the model. Furthermore, in GIMME individual datasets are analyzed separately and sequentially, and common paths are determined by user-specified thresholds. The multi-VAR approach induces sparsity in individual-level models via structured penalization of common and unique dynamics, where all individual-level models are estimated simultaneously. What constitutes a common and unique path in the multi-VAR framework is determined using competing penalties in the objective function and cross-validation. Although additional differences exist, a detailed comparison is beyond the scope of this work. Readers interested in a more technical comparison should see lane2017 for a step-by-step description of the GIMME algorithm and fisher2022 for the multi-VAR algorithm, including pseudocode.

Finally, it is also possible to accommodate both quantitative and qualitative heterogeneity across individual time series model by simply fitting individual-level models separately. In this way each individual-level model is fit without pooling or sharing any information across individuals. In fact, individual-level models represent an important benchmark for any proposed joint-modeling framework. If joint modeling approaches, such as multi-VAR or GIMME, do not show an advantage over individual-level approaches, then sharing of information across individuals provides no additional benefit and one should prefer the most parsimonious modeling approach. For these reasons we also include popular individual-level modeling approaches as comparisons when evaluating the performance of the proposed estimators.

fisher2022 introduced the multi-VAR framework, a method for estimating multiple-subject multivariate time series models using structured regularization. The original proposal contained descriptions of a standard multi-VAR implementation based on the Least Absolute Shrinkage and Selection Operator (LASSO; tibshirani1996) and an adaptive multi-VAR based on the adaptive LASSO (zou2006). In the remainder of the paper we use the term multi-VAR to refer to the general modeling framework proposed in fisher2022, and when relevant for a given context we specifically refer to the standard or adaptive variants.

In this paper we extend the adaptive multi-VAR framework described by fisher2022 in a number of meaningful ways. First, we propose new estimation methods for the adaptive multi-VAR. Specifically, we introduce new methods for computing adaptive weights and a novel cross-validation procedure for selecting the model hyperparameters. Second, we evaluate the utility of these newly proposed weighting schemes and cross-validation procedures in a small simulation study. In this simulation study we evaluate the performance of the proposed approach across various levels of heterogeneity, comparing it to alternative methods, and illustrating its application using an fMRI study example assessing habitual control. Embedded in the paper are illustrations demonstrating how the multi-VAR approach handles commonly encountered modeling situations, with companion code from the multivar package (multivar).

fisher2022 described one approach for solving (3) using a modified Lasso penalty,

where $\|\mathbf{A}\|_{1}$ denotes the $\ell_{1}$ norm of $\mathrm{vec}(\mathbf{A})$. Here we induce sparsity in the individual transition matrices, $\boldsymbol{\Phi}^{k}$, through the decomposition of common and unique effects. Sparsity in the multi-VAR solution is governed by the penalty parameters $\lambda_{1}$ and $\lambda_{2,k}$ chosen using cross-validation. Importantly, heterogeneity of the solution is also determined by the competition of the two penalty parameters.

Here it is helpful to consider three common situations. First, suppose individuals share little in common. In this case a sensible approach to estimating $K$ related VAR models would return essentially independent solutions. Specifically, for large enough values of $\lambda_{1}$, $\hat{\boldsymbol{\Gamma}}^{0}=\mathbf{0}$, and we obtain $\hat{\boldsymbol{\Phi}}^{k}=\hat{\boldsymbol{\Gamma}}^{k}$. Second, suppose there was very little heterogeneity in a given process across individuals. Here, a sensible approach would pool the time series and estimate a single transition matrix. In this case, large enough values of $\lambda_{2,k}$ will lead to $\hat{\boldsymbol{\Gamma}}^{k}=\mathbf{0}$, and we will have $\hat{\boldsymbol{\Phi}}^{k}=\hat{\boldsymbol{\Gamma}}^{0}$. Third, if the heterogeneity in dynamics among the $K$ individuals is unknown in advance, and both common and unique features are plausible, a sensible approach would return a model in which the common $\hat{\boldsymbol{\Gamma}}^{0}$ and unique effects $\hat{\boldsymbol{\Gamma}}^{k}$ are developed in accordance with how similar the $K$ individuals are.

It is also worth noting that just because a path is shared by a subset of individuals does not mean the estimated value of the corresponding parameter will be equal across individuals. In this way multi-VAR accommodates both qualitative and quantitative heterogeneity. See Figure 1 for a simulated toy example of these three settings (in the order of the figure, (1) heterogeneous features, (2) homogeneous features, and (3) heterogeneous and homogenous features) with $T=100$, $K=3$ and $d=10$. All models were estimated using the multivar R package (multivar). Code producing these results is available in the Supplementary Materials.

It is worth noting the standard Lasso procedure is characterized by some important limitations, nicely discussed and summarized in a recent review paper by freijeiro2022. First, for consistent path selection by Lasso, the design matrix needs to satisfy additional assumptions such as the so-called irrepresentable condition (zhao2006). These assumptions require the covariates not to be too strongly dependent, and are difficult to verify in practice. Second, even when the covariates are independent, Lasso is well-known to suffer from an excessive number of false positives. This is a delicate point having to do with the shrinkage procedure itself, thought to add pseudo-noise to the model (see section 2.1.3 in freijeiro2022, and references therein). Moreover, bias can be a problem with $\ell_{1}$ penalization, primarily resulting from the shrinkage applied to coefficient estimates corresponding to the true signal (buhlmann2011). In certain cases, the bias associated with large estimates is equivalent to the penalty parameter used in soft-thresholding. We will revisit this final point in the following section.

zou2006 proposed the adaptive Lasso to overcome some of the limitations associated with the original Lasso procedure. The adaptive Lasso essentially replaces the $\ell_{1}$ penalty with a re-weighted version, where the weights are determined by some initial consistent estimate of the model parameters. This weighting allows for differential penalization across the model coefficients, such that if an initial coefficient estimate is large, a smaller penalty is applied. The opposite holds for initial coefficient estimates that are small in magnitude, where a larger penalty is applied, and a more sparse solution is obtained. Remembering the bias in the $\ell_{1}$ penalty for large coefficients is often proportional to the penalty parameter, it is straightforward to see applying a smaller penalty to large coefficients will reduce the bias of the estimator. Likewise, it is also intuitive that applying larger penalties to coefficients with small initial estimates will reduce the number of false positives.

In addition to the standard approach, fisher2022 also proposed an adaptive version of the multi-VAR problem based on a between-person weighting scheme described by ollier2017. The adaptive weights are built using a preliminary estimate $\tilde{\boldsymbol{\Phi}}^{k}=(\tilde{\phi}_{i,j}^{k})$ according to the following penalty function

where $\Gamma^{0}_{i,j}$ corresponds to $\{i,j\}^{th}$ element of the common effects matrix, $\boldsymbol{\Gamma}^{0}$, and $\Gamma^{k}_{i,j}$ corresponds to the $\{i,j\}^{th}$ element of the unique effects matrix, $\boldsymbol{\Gamma}^{k}$. The divisors $|\tilde{\boldsymbol{\Phi}}_{median}|=(|\tilde{\phi}_{i,j,\textrm{median}}|)$ and $|\tilde{\boldsymbol{\Phi}}^{k}-\tilde{\boldsymbol{\Phi}}_{median}|=(|\tilde{\phi}_{i,j}^{k}-\tilde{\phi}_{i,j,\textrm{median}}|)$ are defined using the entrywise weighted medians of $(\tilde{\boldsymbol{\Phi}}^{1},\ldots,\tilde{\boldsymbol{\Phi}}^{K})$ and $\alpha\geq 1$.

In the regression setting, zhou2009 suggest any initial estimate $\tilde{\theta}$ of $\theta$ with a small bound on $||\tilde{\theta}-\theta||_{\infty}$. In practice, the choice of initial estimates can have a large impact on the bias and overall sparsity of the adaptive Lasso solution. fisher2022 used the unpenalized MLE of $\tilde{\boldsymbol{\Phi}}^{k}$ to construct the adaptive weights. However, other initial estimators have been explored in the literature. For example, zou2006 suggested ridge regression and zhou2009 suggested the Lasso as promising estimators for constructing the initial weights needed for the adaptive approach. Indeed the multi-VAR approach is even more sensitive to the choice of adaptive weights as adjustments are needed for coefficients not immediately estimable from the data (common and unique effects). For this reason it is important to look at how these alternative methods for weight construction impact the adaptive multi-VAR results.

The multi-VAR solution quality depends on the selection of the unknown penalty parameters $\lambda_{1}$ and $\lambda_{2,k}$, $k=1,\dots,K$. In the multi-VAR setting we construct a grid of plausible values for these parameters using the approach described in fisher2022. To identify the optimal parameter values from a grid of plausible values we use cross-validation. fisher2022 proposed adapting a rolling-window cross-validation (RWCV) procedure for high-dimensional VAR models (banbura2010; song2011) to the multiple-subject setting. However, this approach relies on 1-step ahead cross-validation and is computationally expensive. In the current work we propose a computationally efficient adaptation of traditional fold-based cross-validation using blocked sampling (BCV) to the multiple-subject time series context. Although the BCV approach does not preserve the temporal ordering of the component time series (a model can be tested on data that precedes chronologically the training data), it makes more efficient use of the available data compared to RWCV, where many timepoints are essentially removed from testing and training. Numerous authors have found BCV methods work well for single-subject autoregressive models and stationary time series data (bulteel2018a; bergmeir2012; bergmeir2014; bergmeir2018).

Figure 2 provides a visual depiction of the BCV approach. For each value of the $\lambda_{1}$ and $\lambda_{2,k}$ grid we perform the following sequence. First, we split each individual’s times series into $F$ equivalently sized folds. Next, we remove one of the $1,\dots,F$ folds from each individual’s multivariate time series (marked as the Test block in Figure 2). Using all remaining folds (marked as the Train blocks in Figure 2) we solve the problem in (5) using the training blocks from each individual to obtain $\hat{\boldsymbol{\Gamma}}^{0}$ and $\hat{\boldsymbol{\Gamma}}^{k}$. Separately for each individual, these estimates are then used to forecast the testing block and obtain the MSFE. We continue in this fashion, removing each of the $F$ folds, forecasting the omitted test data, and calculating the error, at which time the forecast performance is aggregated across the $K$ individuals and $F$ folds for each combination of $\lambda_{1}$ and $\lambda_{2,k}$ as in

and the values of $\lambda_{1}$ and $\lambda_{2,k}$ which correspond to the smallest MSFE are chosen for the final model.

Although a thorough discussion of solving the multi-VAR optimization problem is beyond the scope of the current paper interested readers should see fisher2022, describing a proximal gradient algorithm based on Fast Iterative Shrinkage-Thresholding Algorithm (FISTA; beck2009). This approach is implemented in the multivar R package (multivar).

Data was generated according to time series length of $30$, $50$, and $100$ for $15$ individuals and $10$ variables, in line with our empirical example. Although our empirical dataset had $100$ timepoints per individual, we also considered shorter time series lengths, as these are also quite common in the literature. Each dataset was then analyzed using the multi-VAR approach, GIMME-VAR and single-subject ($K=1$) methods. The multi-VAR approaches consisted of adaptive multi-VAR with either maximum likelihood, ridge, or lasso adaptive weights. The single-subject approaches consisted of unpenalized maximum likelihood and adaptive Lasso with either maximum likelihood, ridge, or lasso adaptive weights (labeled as $K=1$ approaches in the results). For all penalized approaches we used blocked cross-validation with 10 equal folds and mean-square forecast error (MSFE) calculated on each test fold to choose the hyper parameters. We also compared the blocked cross-validation (BCV) approach described above to the rolling-window cross validation (RWCV) approach described by fisher2022 for the adaptive multi-VAR approaches. This was done to determine whether the proposed BCV method performs at least as well as the established method of RWCV. Additional details on these and other cross-validation procedures for time series can be found in cerqueira2020. Finally, for each cell of the simulation design $750$ datasets were generated and analyzed.

To compare the performance of the selected approaches we considered a number of metrics. To characterize variable selection performance, how well the different approaches recovered the data generating model, we used Matthew’s correlation coefficient (MCC). MCC is a robust single-measure summary of path recovery that gives equal weight to positive and negative cases and incorporates both sensitivity and specificity in its definition (chicco2021matthews; chicco2020advantages). MCC ranges from perfect disagreement ($MCC=-1$) to perfect agreement ($MCC=1$), with a value of $0$ indicating chance performance. For each multiple-subject dataset, the mean MCC is calculated as

where the measure in (7) is averaged over $K$ individuals and TP is the number of nonzero-valued parameters in the data generating model correctly recovered as nonzero-valued in the fitted model, TN is the number of zero-valued parameters in the data generating model correctly recovered as zero-valued in the fitted model, FP is the number of zero-valued parameters in the data generating model incorrectly recovered as nonzero-valued in the fitted model, and FN is the number of nonzero-valued parameters in the data generating model incorrectly recovered as zero-valued in the fitted model.
We also looked at the quality and variability of the estimated coefficients using absolute bias and root means square error (RMSE). For each multiple-subject dataset the mean absolute bias and RMSE were calculated as

where ${\phi}^{k}_{i,j}$ and $\hat{\phi}^{k}_{i,j}$ are the $(k,i,j)^{th}$ elements of the data generating, and estimated transition matrices, respectively, for individual $k$ in a given design condition.

Below we show how one might use adaptive multi-VAR to analyze multivariate time series data from multiple subjects. The goal of this analysis was to better understand the time-dependent relations underlying brain networks thought to be involved in habitual control. Empirical data from $12$ participants ($6$ adolescents, ages $15-17$ years, M_age = 16, SD_age = $0.89$; $6$ adults, ages $25-32$ years, M_age = $28.5$, SD_age = $2.43$) were included in the analyses. $10$ regions of interest (ROIs) were of interest as they represent 3 functional networks derived from the Power $264$ atlas (power2011). These networks include the cingulo-opercular task control network, the fronto-parietal task control network, and the default mode network (DMN) as they are broadly implicated in habitual control (fox2005; uddin2009). Time series data from each ROI were concatenated across participants to produce a single matrix with $10$ columns representing each ROI and $100$ timepoints.

Below we include the code used to generate Figure 1.

To better illustrate the data generation algorithm presented earlier it may be helpful to consider a single condition in more detail. Let us supposed a 10-dimensional VAR(1) model with 15 individuals. Let us also consider the high heterogeneity condition where we have path proportions vector, $\pi_{P}=(1/5,1/10,1/20)$, and individual proportions vector, $\pi_{I}=(1/3,2/3,1)$. We will simultaneously iterate through these two vectors in the following manner to construct the individual transition matrices a given instance of the high heterogeneity condition.
For iteration 1 we begin by taking the first element from the path proportions vector, $\pi_{P,1}=1/5$. This number indicates the proportion of nonzero paths that will be shared by the proportion of individuals in $\pi_{I,1}$. As we have $100$ possible nonzero paths in an arbitrary $10\times 10$ transition matrix in this example, we randomly select the location of $20$ nonzero elements ($\pi_{P,1}=1/5=20/100$). We now assign those nonzero paths to $\pi_{I,1}=1/3$ of our sample. As there are $15$ individuals in our current example, we randomly select $5$ individuals ($\pi_{I,1}=1/3=5/15$) and assign those $5$ individuals to the same $20$ nonzero paths. Note, although the location of these $20$ non-zero paths are uniform across the $5$ individuals selected, the numeric values of the non-zero paths are distinct.
In iteration 2 we take the second element from the path proportions vector, $\pi_{P,2}=1/10$ and the second element from the individual proportions vector $\pi_{I,2}=2/3$. Using these proportions we randomly select $10$ new nonzero paths that were not previously selected in iteration 1 ($\pi_{P,2}=1/10=10/100$). We then assign these nonzero paths to $10$ randomly selected individuals ($\pi_{I,2}=2/3=10/15$). Note, the randomly selected individuals from iteration 2 may overlap with the individuals selected in iteration 1.
In iteration 3 we take the third and final element from the path proportions vector, $\pi_{P,3}=1/20$ and the second element from the individual proportions vector $\pi_{I,3}=1$. Using these proportions we randomly select $5$ new nonzero paths ($\pi_{P,3}=1/20=5/100$) that were not previously selected in iteration 1 or 2, and assign them to all $15$ individuals in the sample ($\pi_{I,3}=1=15/15$). Again we note that although the location of these non-zero paths are uniform across the individuals selected, the numeric values of those non-zero paths are distinct across individuals.
Finally, for all individuals separately, choose parameter values at random from the uniform distribution $U(0.1,0.9)$ for all nonzero paths identified in the previous iterations and, if needed, re-scale the resulting transition matrices, $\boldsymbol{\Phi}^{k}$ to ensure a stable solution. Using their respective path and individual proportion vectors this exact procedure is then followed to construct the data generating matrices for each replication of no and low heterogeneity conditions.

In the paper we occasionally reference sensitivity and specificity when it provides additional insight into path recovery performance. Sensitivity and specificity are measures of the probability a path is recovered, given it was in the data generating model, and the probability a path is not recovered, given it was not in the data generating model, respectively. Here we show how these measures were calculated,

where aggregate measures were obtained by averaging across Monte Carlo iterations. Here, as in the calculation of MCC above, TP is the number of nonzero-valued parameters in the data generating model correctly recovered as nonzero-valued in the fitted model, TN is the number of zero-valued parameters in the data generating model correctly recovered as zero-valued in the fitted model, FP is the number of zero-valued parameters in the data generating model incorrectly recovered as nonzero-valued in the fitted model, and FN is the number of nonzero-valued parameters in the data generating model incorrectly recovered as zero-valued in the fitted model.

Simulation results for the adaptive multi-VAR models using both the block and rolling-window cross-validation approaches are provided in Figure 5. No discernible pattern of differences emerged across the two cross-validation approaches for the set of simulation conditions considered.