Penalized Subgrouping of Heterogeneous Time Series
Christopher M. Crawford Jonathan J. Park Sy-Miin Chow Anja F. Ernst Vladas Pipiras Zachary F. Fisher

Emergence of two inertial sub-ranges in solar wind turbulence: 
 dependence on heliospheric distance and solar activity
    [
    
School of Physics, Nankai University, Tianjin, 300071, China 

=======================================================================================================================



Recent technological advances have decreased the burden associated with collecting intensive longitudinal data (ILD) in the social, behavioral, and health sciences. The availability of such data has catalyzed interest in the study of constructs defined by the complex interplay between dynamic biopsychosocial processes. Despite these increases in both access to ILD and interest in the study of dynamic processes Hamaker2017, how best to model multivariate time series data arising from multiple individuals is still an open question. Central to this question is how researchers should best accommodate the persistent heterogeneity observed in many aspects of human behavior.  

Current methods for analyzing dynamic processes vary in the degree to which this heterogeneity is addressed [e.g.,][]Liu2023. Multilevel modeling approaches, for example, allow individuals to differ quantitatively on a limited set of dynamic features through the inclusion of random effects [e.g.,][]Bringmann2013. However, because standard multilevel models assume no qualitative differences in the pattern of relations among the dynamic processes, these approaches may be overly restrictive. For example, data generating processes characterized by individual differences in the patterning of zero and nonzero dynamics are poorly represented by such approaches. Moreover, violation of assumptions regarding the structure and distribution of random effects can adversely impact estimation and inference McNeish2017. Idiographic methods, conversely, allow for a great deal of flexibility through the specification of person-specific models (<cit.>; <cit.>). Indeed, recent work has shown that the dynamics of commonly studied constructs—such as depression and anxiety Fisher2017, externalizing behavior Wright2015, and personality <cit.>—are characterized by marked heterogeneity in both the magnitude of the features and the patterning of dynamics, suggesting that person-specific methods may be needed to fully capture the complexity of the data generating processes. The increased flexibility afforded by idiographic approaches, however, may limit generalizability and inhibit the identification of shared dynamics critical for areas such as intervention development.

Recently, a number of approaches have emerged that offer alternatives for characterizing heterogeneity in dynamic processes by seeking to bridge the divide between nomothetic and idiographic approaches to modeling multiple-subject, multivariate time series. One of these approaches, the multi-VAR framework (<cit.>; <cit.>), is built upon the vector autoregressive (VAR) model and simultaneously estimates group- and individual-level models. Importantly, the multi-VAR approach accommodates both quantitative and qualitative heterogeneity in dynamics across individuals—that is, heterogeneity in both the magnitude and pattern of zero and nonzero dynamics—and is compatible with a number of penalization methods for structuring how information is shared across individuals. Currently, multi-VAR is only capable of estimating a single group-level model, presumably by all individuals. It may be the case, however, that for many processes shared patterns of dynamics also exist among subgroups, or clusters, of individuals. 

To address this limitation, the current project extends the multi-VAR framework to allow for data-driven identification of subgroups and penalized estimation of subgroup-level dynamics. The approach detailed herein is characterized by a number of advantageous features. First, both quantitative and qualitative heterogeneity are accommodated in the estimation of group-, subgroup-, and individual-level models. In contrast to most existing subgrouping frameworks, this allows for individual differences in the magnitude of dynamics both within and between subgroups. Moreover, these effects are not assumed to follow any specific distributional form, thereby capitalizing on the flexibility inherent in fully idiographic approaches while also prioritizing generalizability through the group- and subgroup-level models. Second, sparsity in the group-, subgroup-, and individual-level dynamics is induced through a penalized estimation procedure, addressing challenges associated with overparameterization in both single-subject [e.g.,][]Sims1980 and multilevel approaches. Indeed, recent work on the development of a subgrouping approach incorporating finite mixture modeling into the multilevel VAR framework notes that estimation can be burdensome when specifying a large number of random effects Ernst2024. Finally, the simultaneous estimation procedure for group- and individual-level dynamics Fisher2022 extends to subgroup-level effects. Conversely, iterative approaches to the identification and estimation of subgroup dynamics may encounter issues associated with similarly iterative variable selection methods, including overfitting and suboptimal solutions McNeish2015. These features are described in detail in the following sections.

The remainder of the article is organized as follows. First, we review existing methods for subgrouping multiple-subject, multivariate time series. We then introduce the multi-VAR framework and the penalized subgrouping extension. We conclude with a presentation of results from both a simulation study and empirical example conducted to examine the efficacy and utility of the multi-VAR subgrouping extension. 



§ SUBGROUPING METHODS FOR TIME SERIES


The identification of subgroups of individuals for whom shared patterns of dynamics exist, and the estimation of parameters that provide meaningful information about said subgroups, are of substantive interest in a variety of domains. Recent work, for example, has noted that current taxometric procedures for the diagnosis of psychopathological syndromes fail to account for the observed heterogeneity within putatively homogeneous diagnostic categories Kotov2017. Appropriately parsing such heterogeneity can accelerate the development of more effective treatment paradigms [e.g.,][]Fisher2019, echoing calls for more personalized approaches to diagnosis and treatment in both clinical science Wright2020 and medicine Hamburg2010. Importantly, however, there is not a clear consensus regarding best practices for the identification of subgroups in multiple-subject, multivariate time series. Liao (), for example, notes that time series clustering methods can be organized according to which aspects of the data are used in the subgrouping procedure. For the purposes of the current work, we focus our attention on methods that derive clusters based on parameter estimates in the VAR framework.

One popular VAR-based clustering method is the alternating least squares (ALS) approach Bulteel2016, an iterative algorithm consisting of three steps. First, individuals are assigned to initial clusters, which can be done randomly or through the use of alternative clustering approaches, such as hierarchical clustering using Ward's criterion Ward1963. Next, a VAR model is fit to each cluster to obtain subgroup-specific parameter estimates. Finally, subgroup membership is updated by reassigning individuals to the cluster for which the sum of their squared prediction errors are minimized. These steps are repeated until there are no changes in subgroup membership. A notable limitation of this framework is that only between-cluster quantitative heterogeneity is modeled—that is, the ALS approach assumes a priori that individuals within a cluster are governed by the same data generating process, and that there are no qualitative differences between clusters. Thus, different clusters are assumed to share the same pattern of zero and nonzero dynamics.

Work by Ernst and colleagues () relaxes the within-group homogeneity assumption with a VAR-based clustering approach that incorporates the Gaussian finite mixture model (GMM; <cit.>; <cit.>). Proceeding in a two-step fashion, person-specific VAR models are first fit to data for each individual, with the restriction that there are no between-person differences in lag order. A GMM is then applied to the person-specific parameter estimates, thereby allowing for quantitative differences in within-cluster estimates. With this added flexibility, however, comes increased complexity in both the estimation of the parameters that define the multivariate Gaussian distribution and selection of the model that best represents the structure of the data. Indeed, the use of relative fit statistics to compare models in the GMM framework often requires bootstrap-based methods Mclachlan2000. Simulation studies comparing these two VAR-based clustering methods have found that the GMM approach outperforms the ALS framework when within-group quantitative differences are present Ernst2021, whereas the opposite is true when the homogeneity assumption is met Takano2021.

Another approach is the subgrouping extension within the group iterative multiple model estimation (S-GIMME) algorithm Gates2017. The S-GIMME algorithm can be described in three sequential stages. In the group stage, the S-GIMME framework estimates a structural VAR model Lutkepohl2005 for each individual. Modification indices Sorbom1989 are used to identify paths to be added to the model if they significantly improve model fit for a majority of individuals, thereby returning estimates that are shared across individuals. Next, the subgroup stage uses these group-level results to construct a similarity matrix such that each element of said matrix represents the number of parameter estimates shared between two individuals in terms of sign and significance. The community detection algorithm Walktrap Pons2006 is then applied to the similarity matrix to identify clusters of individuals with shared dynamics. Once subgroups have been identified, estimates for individuals within each subgroup are obtained in the same manner as in the group stage. Finally, the third stage estimates additional dynamics for each individual by iteratively adding paths until the model is deemed acceptable via commonly used fit indices and associated cutoff values Lane2019. The S-GIMME approach therefore returns estimates for each individual and whether each estimate is unique, subgroup-, or group-specific, thereby accommodating both quantitative and qualitative heterogeneity. Moreover, challenges associated with the estimation of VAR processes comprised of many variables are addressed through the use of stopping criteria at each stage, which promotes sparsity and parsimony. However, forward-selection procedures, such as S-GIMME, can be limited by their sequential nature, and thresholds for the fit indices used in the individual-stage search procedure have been shown to be inadequate in many settings Mcneish2023.

Finally, a recently developed subgrouping approach, the subgrouped chain graphical VAR (scGVAR), extends the graphical VAR framework epskamp2018 by estimating group-, subgroup-, and individual-level dynamics in a three-stage process Park2024. First, individual-level results are obtained by estimating a graphical VAR model for each individual. Next, in a procedure similar to that employed by the S-GIMME framework, results from the first stage are used to create an adjacency matrix quantifying structural similarities in estimated dynamics across individuals. However, whereas S-GIMME defines each element of said matrix as the number of shared nonzero paths between each pair of individuals, scGVAR includes shared zero paths in this definition. Following the construction of this similarity matrix, the Walktrap Pons2006 algorithm is used to identify subgroups of individuals. These subgroup assignments are then used in the estimation of subgroup-level dynamics, wherein a graphical VAR model is fit to the aggregated time series of individuals in each subgroup. Group-level dynamics are then obtained by fitting a single graphical VAR model to the aggregated time series of all individuals. Similar to S-GIMME, scGVAR is therefore able to accommodate both quantitative and qualitative heterogeneity through the estimation of group-, subgroup-, and individual-level dynamics. However, whereas the S-GIMME procedure allows group- and subgroup-level dynamics to vary in magnitude between individuals, these effects are fixed in the scGVAR approach—that is, all individuals in the same subgroup are characterized by identical parameter estimates. Moreover, the two frameworks differ with respect to the challenges associated with estimation in high-dimensional settings. Indeed, in contrast to the use of fit indices in S-GIMME, scGVAR relies on a penalized estimation procedure to promote sparsity and parsimony.



§ MULTI-VAR


The motivation for the development of the multi-VAR framework is twofold. First, ILD collected in the social and behavioral sciences are generally composed of data from multiple individuals. From the single-subject VAR perspective, this now involves the estimation of transition matrices for each individual. It is unlikely, however, that the dynamic processes of interest are entirely heterogeneous across all individuals in a given sample. Indeed, a key feature of common ILD modeling approaches, such as multilevel modeling, is that information is shared across individuals in the parameter estimation procedure. A purely idiographic approach, wherein each of the individual transition matrices are estimated independently, disregards this shared information. Thus, the single-subject VAR, though capable of modeling both qualitative and quantitative heterogeneity, is unable to account for shared dynamics and may be suboptimal in many commonly encountered data settings where sharing information across individuals is beneficial (e.g., when time series length is small).

Second, as noted previously, the single-subject VAR is plagued by profligate parameterization Sims1980, such that the number of parameters grows quadratically with each additional component series. Given typical sample sizes (i.e., number of time points) in the social, behavioral, and health sciences [e.g.,][]Rot2012, this could result in the estimation of a large number of parameters given the available data, with such estimates lacking precision Bruggemann2012. Solutions to this dimensionality issue in the single-subject VAR have focused on reducing the parameter space using methods ranging from sequential search procedures to regularization via the least absolute shrinkage and selection operator (Lasso; <cit.>). Central to these remedies is the assumption that the true model is sparse in nature. The feasibility of this assumption is open to interpretation; however, the bet on sparsity principle Hastie2001 argues that methods assuming sparsity are preferable because if this assumption is false—that is, the true model is dense—then existing approaches are unable to recover the true model without a large amount of data. The multi-VAR framework aims to address the challenges associated with common modeling approaches, such as the single-subject VAR and multilevel modeling, while retaining their desirable features.

To motivate the construction of the multi-VAR framework, we first consider a multivariate (d-variate) time series for a single individual, {𝐗_t}_t ∈ℤ = {(X_j,t)_j=1,…,d}_t ∈ℤ, where 𝐗_t follows the canonical VAR model of order p, VAR(p), if


    𝐗_t = 
      Φ_1𝐗_t-1 + 
      … + 
      Φ_p𝐗_t-p + 
      𝐄_t,    t ∈ℤ,


for d × d transition matrices Φ_1, …,Φ_p containing autoregressive and cross-regressive parameters and a white noise series {𝐄_t}_t ∈ℤ∼WN(0, Σ_𝐄) characterized by 𝔼(𝐄_t)=0 and 𝔼(𝐄_t𝐄^'_s)=0 for s ≠ t. We assume here that 𝐗_t is of zero mean for simplicity, though developments that follow can easily accommodate time series with nonzero means. Generally, a unique causal stationary solution to (<ref>) can be ensured when the roots of det( Φ(z)), where Φ(z)=𝐈_d -Φ_1z - … - Φ_p z^p, all have moduli greater than unity. With observations 𝐗_1, …, 𝐗_T, we can more concisely express (<ref>) in the familiar regression format:  


    𝐘=Φ𝐙 + 𝐔,


where 𝐘 = (𝐗_p+1, …, 𝐗_T) is a d × (T-p) outcome matrix, Φ = (Φ_1, …, Φ_p) is a d × (dp) transition matrix, 𝐙 is a (dp) × (T-p) design matrix, and 𝐔 is a d × (T-p) matrix of process noise. In the remainder of this work we only consider first-order VAR models—that is, VAR(1). However, all arguments can be extended to accommodate models with arbitrary lag orders without any loss of generality.

The multi-VAR framework extends the above single-subject VAR representation to readily accommodate multiple-subject, multivariate time series through the estimation of Φ^1, …,Φ^K sparse transition matrices for K individuals, each composed of common and unique effects. To do so, the multi-VAR approach relies on the following decomposition of Φ^k,


    Φ^k = Γ + Υ^k, 
         k = 1, … , K,


where Γ∈ℝ^d × d corresponds to the common effects across K individuals and Υ^k∈ℝ^d × d represents the unique effects for individual k. This decomposition allows for heterogeneity in the structure of the dynamics through the inclusion of unique, person-specific effects while allowing for some degree of homogeneity in the common effects. Further, as no distributional assumptions are imposed on (<ref>), individual transition matrices are free to vary both quantitatively and qualitatively across individuals. Notably, shared paths in the common effects matrix are allowed to vary in magnitude—that is, the effects are freely estimated while preserving the structure of Γ for all K individuals.

One approach for sparse estimation of (<ref>) was proposed by Fisher and colleagues () using the Lasso penalization paradigm,


    Γ, Υ^1,...,Υ^Kargmin1/N∑_k=1^K 𝐘^k - (Γ + Υ^k) 𝐙^k _2^2 + P_standard
    
        P_standard = λ_1 Γ_1 + ∑_k = 1^K λ_2,kΥ^k _1,


where P_standard is the standard Lasso penalty, 𝐀_1 denotes the ℓ_1 norm of vec(𝐀), and N = (T - p). Sparsity and heterogeneity in the multi-VAR solution is therefore determined and governed by the two penalty parameters, λ_1 and λ_2,k, which are chosen using cross-validation. The inclusion of penalty parameters on both the common and unique effects allows for flexibility in the approximation of potentially heterogeneous data generating processes. If, for example, individuals share little in common, it would be expected that the solution would be equivalent to estimating independent, single-subject VAR models. This would correspond to large values of λ_1, such that Γ̂ = 0 and the solution returns Φ̂^k = Υ̂^k. Conversely, if individuals are highly homogeneous, large values of λ_2,k would result in Υ̂^k = 0 and Φ̂^k = Γ̂, paralleling approaches that pool the time series and estimate a single transition matrix. It is likely that the reality falls between these two extremes, thereby resulting in a solution in which the common and unique effects reflect the degree of heterogeneity observed. 

Despite widespread use of the standard Lasso penalization framework, it is characterized by several important limitations [for a review, see][]Freijeiro‐González2022. Indeed, it is well known that the Lasso exhibits drawbacks with respect to consistent path selection, excessive false positives, and bias in many situations, such as when the variables are strongly dependent. The adaptive Lasso—developed, in part, to address these limitations—replaces the standard ℓ_1 penalty term with a re-weighted version, where the weights are determined by consistent initial estimates of the model parameters  Zou2006. This weighting enables differential penalization across parameters of interest, such that larger initial estimates correspond to smaller weights (and vice versa for small initial estimates). Using this idea, the adaptive multi-VAR was also proposed by Fisher and colleagues (), wherein the weights are constructed using initial estimates Φ^k as in


    P_adaptive = λ_1 ∑_i,j=1^d 1/|ϕ_i,j,median|^α
        |Γ_i,j| + ∑_k=1^K λ_2,k∑_i,j=1^d
        1/|ϕ_i,j^k - ϕ_i,j,median|^α|Υ_i,j^k|,


where Γ_i,j corresponds to the {i,j}^th element of Γ, Υ_i,j^k corresponds to the {i,j}^th element of Υ^k, and similarly with ϕ_i,j^k for Φ^k and ϕ_i,j,median for Φ_median, with Φ_median representing the matrix of median coefficient estimates for all K individuals, and α≥ 1. Substituting the standard Lasso penalty, P_standard, in (<ref>) for the adaptive Lasso, P_adaptive, in (<ref>) results in the adaptive multi-VAR objective function. Prior work has shown that multi-VAR with adaptive Lasso performs well across a range of factors Fisher2024. 



§ SUBGROUPING MULTI-VAR


The subgrouping multi-VAR procedure consists of two steps. First, subgroup enumeration and classification occur. That is, prior to estimation of subgroup-specific effects, the number of relevant subgroups is determined, and each individual is assigned to a subgroup. Next, this information is used in both the construction of the design matrix, Z, and the decomposition and penalized estimation of the individual transition matrices, Φ^k. These two steps are described in detail in the following subsections.



 §.§ Identifying Subgroups


To identify the number and membership of clusters, subgrouping multi-VAR first estimates Φ^k for K individuals using the standard multi-VAR framework. The individual-level effects from these matrices, Υ^k, are then used to construct a K × K similarity matrix, where the off-diagonal elements represent the number of shared dynamics between each pair of individuals in terms of presence and sign. If, for example, a specific path is both nonzero and of the same sign (i.e., positive or negative) for individuals i and j, then the {i,j}^th element of the similarity matrix increments by one.  The use of individual-level effects in the construction of the similarity matrix ensures that deviations from common effects contribute to subgroup identification. The construction of the similarity matrix in subgrouping multi-VAR is similar to the procedures implemented in S-GIMME and scGVAR (<cit.>; <cit.>). 

The community detection algorithm Walktrap Pons2006 is then applied to this similarity matrix. The Walktrap algorithm uses the information in the similarity matrix to compute a transition matrix, such that each element corresponds to the probability of transitioning from one individual to another for a random walk of a given length. Intuitively, this allows the Walktrap algorithm to identify densely connected areas (i.e., communities). Ward's () hierarchical clustering procedure is then used to determine the optimal number of clusters by iteratively merging communities until all individuals are in a single cluster. Subgroup enumeration then proceeds by identifying the configuration with the maximum modularity Newman2004, which indicates the degree to which individuals within a cluster are similar relative to those from different clusters. The Walktrap algorithm has been found to perform well when applied to count matrices Gates2016, and has been successfully implemented in similar VAR-based methodological frameworks (e.g., <cit.>; <cit.>).



 §.§ Estimating Subgrouping Effects


To extend the standard multi-VAR framework to allow for estimation of subgroup-specific effects, we consider the further decomposition of Φ^k,


    Φ^k = Γ + Π^s +
        Υ^k,
         s = 1, … , S,
         k = 1, … , K,


where Γ and Υ^k continue to correspond to the common and unique effects, respectively, and Π^s∈ℝ^d × d represents the subgroup effects for subgroup s. Paralleling the initial decomposition in (<ref>), individual transition matrices in (<ref>) are free to vary both in structure and magnitude across individuals. This decomposition can be further incorporated into the objective function detailed in (<ref>) within the standard Lasso penalization framework


    Γ, Π^1, …,Π^S, …, Υ^Kargmin1/N∑_k=1^K 𝐘^k - (Γ + Π^s + Υ^k) 
        𝐙^k _2^2 + P_standard
    
        P_standard = λ_1 Γ_1 + 
        ∑_s = 1^S α_sΠ^s _1 +
        ∑_k = 1^K λ_2,kΥ^k _1,


where the addition of α_s in P_standard indicates that the sparsity and heterogeneity of the solution is now determined by three penalty parameters, λ_1, α_s, and λ_2,k. Importantly, this suggests that the competition of the three penalty terms affords even greater flexibility in approximating the underlying data generating processes. In addition to the scenarios noted previously, for example, the subgrouping multi-VAR can accommodate data characterized by within-subgroup homogeneity and a high degree of between-subgroup heterogeneity, resulting in large values of λ_1 and λ_2,k, such that Φ̂^k = Π̂^s.

This decomposition can also be incorporated into the penalty function in (<ref>) for the adaptive multi-VAR


    P_adaptive =    λ_1  ∑_i,j=1^d 1/|ϕ_i,j,median|^α
        |Γ_i,j|  + 
       ∑_s=1^S α_s∑_i,j=1^d
        1/|ϕ_i,j,median^s|^α|Π_i,j^s|  + 
       ∑_k=1^K λ_2,k∑_i,j=1^d
        1/|ϕ_i,j^k - ϕ_i,j,median|^α|Υ_i,j^k|,


where Π_i,j^s corresponds to the {i,j}^th element of Π^s, and similarly for ϕ^s_i,j,median and Φ^s_median. In the current project—as with the standard multi-VAR—the weights in the adaptive Lasso are constructed using initial estimates Φ^k, Φ_median, and Φ^s_median, with Φ^s_median representing the matrix of median coefficient estimates across all individuals in a given subgroup. Note that for the oracle properties of the adaptive Lasso to hold—that is, identification of the correct subset model and optimal estimation rate—consistent estimates for the weights must be chosen Zou2006. In practice, however, selecting an appropriate estimator can be nontrivial in many situations, such as high dimensional settings or when covariates are strongly dependent. In such cases, estimation of adaptive weights via ℓ_1 or ℓ_2 penalization, respectively, may convey benefits over standard procedures, such as ordinary least squares or maximum likelihood. For more information regarding the performance and availability of various adaptive weights estimators in the multi-VAR framework, see Fisher et al. ().

Selection of the unknown penalty parameters—λ_1, α_s, and λ_2,k—in the subgrouping multi-VAR framework is done using a blocked cross-validation approach [BCV;][]Bulteel2018, which proceeds as follows. First, each of the K multivariate time series is divided into F equally sized folds. Next, one of the F folds is removed from each time series to serve as the testing block. Using the remaining folds comprising the training block, estimates of Γ, Π^s, and Υ^k are obtained, which are subsequently used to predict the testing block and obtain the mean squared error [MSE; though see][for an alternative to the MSE]Revol2024. This procedure continues such that each of the F folds serves as a testing block on which the MSE is computed. The performance is then aggregated across the K individuals and F folds for each combination of λ_1, α_s, and λ_2,k, wherein total error is calculated as 


    MSE_λ_1,α_s,λ_2,k =
        1/K∑_k=1^K 1/F∑_f=1^F
        Ŷ_f^k - Y_f^k _2^2,


and the values of the penalty parameters corresponding to the smallest total MSE are selected for the final model. Note that though the subscripts associated with penalty parameters α_s and λ_2,k suggest the possibility of subgroup- and person-specific penalization, the current project only examines a single parameter for each—that is, α_s = α and λ_2,k = λ_2. This more general notation is nonetheless included to emphasize the flexibility of the subgrouping multi-VAR framework.



§ SIMULATION STUDY


To evaluate the performance of the subgrouping multi-VAR framework, we conducted a Monte Carlo simulation study, wherein factors of interest were selected to represent scenarios typically encountered in ILD applications. Factors examined in the current study included number of individuals, K = (50, 100); length of individual time series, T = (50, 100); number of subgroups, S = (2, 3); and composition of the subgroups, balanced or unbalanced. For each condition, 50 replications were generated, yielding a fully factorial design with 16 conditions (2 × 2 × 2 × 2) and 800 (16 × 50) unique data sets. In the K = 50 conditions, for example, Monte Carlo estimates were computed using 2500 (50 × 50) multivariate time series. Number of variables was held constant at d = 10 for each condition. The performance of the subgrouping extension was assessed in comparison to alternative methods for modeling multi-subject, multivariate time series, namely, S-GIMME and scGVAR (<cit.>; <cit.>). As a benchmark, performance was also compared to the standard multi-VAR framework and multi-VAR with confirmatory subgrouping—that is, subgrouping multi-VAR when subgroup membership is known a priori.

The motivation for the selection of the comparison methods, S-GIMME and scGVAR, assessed in the current simulation study was twofold. First, like subgrouping multi-VAR, both frameworks seek to accommodate quantitative and qualitative heterogeneity through the estimation of group-, subgroup-, and individual-level dynamics. Second, both methods address the challenges associated with modeling VAR processes comprised of many variables by inducing sparsity and parsimony in the estimation procedure. As noted previously, however, the manner in which these features are incorporated varies considerably between approaches. For example, whereas subgrouping multi-VAR estimates group-, subgroup-, and individual-level dynamics simultaneously, S-GIMME and scGVAR are characterized by iterative and stage-based estimation procedures. In addition to assessing the performance of subgrouping multi-VAR across a range of factors, the current simulation study therefore represents a rigorous evaluation of competing approaches for the analysis of multiple-subject, multivariate time series. 



 §.§ Data Generation and Model Estimation


Across all conditions, 10 × 10 transition matrices were generated for each individual in the following manner. First, subgroup membership was specified. For the balanced conditions, subgroups were constructed such that each contained the same number of individuals. For the unbalanced conditions with two subgroups, 30% of individuals were placed into the first subgroup, and the remaining 70% were specified as members of the second subgroup. For the unbalanced conditions with three subgroups, the first and second subgroups consisted of 20% of individuals, and the third subgroup contained the remaining 60%. Next, common effects, Γ, were specified as the 10 diagonal elements of each transition matrix—that is, the autoregressive effects. The location of subgroup, Π^s, and unique, Υ^k, effects were then chosen at random, wherein the number of effects for each represented 5% of possible paths. Thus, each transition matrix consisted of 10 common effects, five subgroup effects, and five unique effects, yielding a sparse matrix with 20% nonzero entries. All effects were drawn from a 𝒰(0,1) distribution until a stationary solution was obtained, as determined by the stationarity condition noted previously. The data generation procedure therefore incorporated both qualitative and quantitative heterogeneity through variation in the location of the subgroup and unique effects and variation in the magnitude of group, subgroup, and unique effects, respectively. 

All models analyzed within the multi-VAR framework—that is, subgrouping multi-VAR, multi-VAR with confirmatory subgrouping, and standard multi-VAR—were estimated using the adaptive Lasso procedure detailed above, with adaptive weights obtained via ℓ_1 penalization. Paralleling the group-level stage in the S-GIMME framework, however, all autoregressive paths were assumed to be nonzero, and were therefore not subjected to penalization in either the subgroup enumeration or estimation stages. All multi-VAR and S-GIMME analyses were conducted using the  multivar and  gimme R packages, respectively.



 §.§ Outcome Measures


The performance of the subgrouping multi-VAR framework and comparison methods was evaluated in three ways: model recovery, quality of estimated effects, and accuracy of subgroup identification. A number of metrics were used to assess model recovery, including sensitivity, specificity, and Matthews correlation coefficient [MCC;][]Matthews1975. Sensitivity and specificity can be interpreted as the true positive and true negative rates, respectively, and were computed as follows: 


    Mean sensitivity = 1/K∑_k = 1^K( ∑_i,j (ϕ̂_i,j^k ≠ 0  and ϕ_i,j^k ≠ 0 ) /∑_i,j (ϕ_i,j^k ≠ 0) ),



    Mean specificity = 1/K∑_k = 1^K( ∑_i,j (ϕ̂_i,j^k = 0  and ϕ_i,j^k = 0 ) /∑_i,j (ϕ_i,j^k = 0) ),


where ϕ̂_i,j^k and ϕ_i,j^k correspond to the {i,j}^th element of the estimated and true transition matrices for the k^th individual, respectively. Mean sensitivity and specificity were then averaged across all replications, resulting in sensitivity and specificity values for each condition. MCC similarly provides a metric for assessing the recovery of the data generating model by incorporating both sensitivity and specificity in its calculation:


    Mean MCC = 1/K∑_k = 1^K
        ( TP_k × TN_k - FP_k × FN_k/√((TP_k + FP_k)(TP_k + FN_k)(TN_k + FP_k)(TN_k + FN_k))),


where TP is the number of parameters correctly estimated as nonzero, TN is the number of parameters correctly estimated as zero, FP is the number of parameters incorrectly estimated as nonzero, and FN is the number of parameters incorrectly estimated as zero. Notably, MCC ranges from perfect disagreement between estimated and true models (MCC = -1) to perfect agreement between the two (MCC = 1), and can therefore be interpreted as a discretization of the standard correlation coefficient in the context of binary classification Boughorbel2017. As such, the MCC can be interpreted with respect to commonly employed benchmarks [e.g.,][]Cohen1988.

To evaluate the quality and variability of the estimated parameters, we computed the mean absolute bias and root mean squared error (RMSE):


    Mean Absolute Bias = 1/K∑_k=1^K 1/d^2∑_i,j=1^d
        |ϕ̂_i,j^k - ϕ_i,j^k|,



    RMSE = 1/K∑_k=1^K √(1/d^2∑_i,j=1^d
        (ϕ̂_i,j^k - ϕ_i,j^k)^2),


where ϕ̂_i,j^k and ϕ_i,j^k again correspond to the {i,j}^th element of the estimated and true transition matrices for the k^th individual, respectively, and d represents the number of variables (i.e., 10). These values were then averaged over all replications to provide mean absolute bias and RMSE metrics for each condition.

To assess the accuracy of subgroup identification, we computed the Hubert-Arabie adjusted Rand index [ARI;][]Hubert1985:


    ARI = K  2 (a+d) - [(a+b)(a+c) + (c+d)(b+d)]/K  2^2 - [(a+b)(a+c) + (c+d)(b+d)],


where a represents the number of pairs of individuals correctly placed in the same cluster, b is the number of pairs incorrectly placed in different in different cluster, c indicates the number of pairs incorrectly placed in the same cluster, and d is the number of pairs correctly placed in different clusters. The ARI therefore incorporates information about the number of true positive, false negative, false positive, and true negative classifications. The ARI has an upper bound of 1, which corresponds to perfect subgroup identification. ARI values of 0, conversely, indicate that subgroup assignments were equal to chance. ARI values greater than or equal to 0.90 were considered excellent, values between 0.80 and 0.89 were considered good, values between 0.65 and 0.79 indicated moderate recovery, and values below 0.64 suggested poor subgroup identification Steinley2004. Monte Carlo errors (MCE), defined as the standard deviation of the Monte Carlo estimates across all replications, were computed to quantify the uncertainty in all outcome measure estimates Koehler2009.



 §.§ Simulation Results




  
Model Recovery

Sensitivity, specificity, and MCC values across all conditions for each estimator can be seen in Tables <ref> and <ref>, and are visualized in Figure <ref>. Given that MCC incorporates both sensitivity and specificity, assessment of model recovery focused on this metric. Across all conditions examined in the current study, subgrouping multi-VAR demonstrated more accurate model recovery than S-GIMME, scGVAR, or standard multi-VAR, and displayed only a slight decrement in performance when compared to the confirmatory subgrouping framework, where the true subgroup membership is assumed to be known. When K = 50 and S = 2, mean MCC for subgrouping multi-VAR increased as a function of time series length, with values of 0.81 and 0.88 when T = 50 and T = 100, respectively, for balanced subgroups, and values of 0.79 and 0.84 when subgroup membership was unbalanced. Subgrouping multi-VAR therefore demonstrated good model recovery for these conditions, with some MCC values reflecting near perfect agreement between true and estimated models. Differences in model recovery between multi-VAR with data-driven subgrouping and confirmatory subgrouping were minimal across these conditions (less than 0.02). Model recovery for S-GIMME, scGVAR, and standard multi-VAR similarly increased as a function of time series length. Mean MCC values for standard multi-VAR were 0.68 and 0.77 when T = 50 and T = 100, respectively, when subgroups were balanced, and 0.69 and 0.78 when subgroups were unbalanced. Mean MCC values for S-GIMME were 0.66 and 0.80 for both balanced and unbalanced subgroups, whereas values for scGVAR were 0.56 and 0.64 for balanced subgroups and 0.55 and 0.63 for unbalanced subgroups. 

Similar results were observed when K = 50 and S = 3. Indeed, mean MCC values for subgrouping multi-VAR increased from 0.80 to 0.88 as T increased from 50 to 100 for subgroups with balanced membership, and from 0.78 to 0.84 for unbalanced subgroups. Observed results therefore suggest that model recovery for subgrouping multi-VAR was largely unaffected by the addition of a third subgroup when K = 50. Discrepancies between data-driven and confirmatory subgrouping approaches were again minimal across these conditions. Mean MCC values for standard multi-VAR were 0.68 and 0.78 when subgroups were balanced, and 0.68 and 0.76 when membership was unbalanced. Comparable metrics were observed for both S-GIMME and scGVAR. Thus, model recovery for subgrouping multi-VAR and comparison approaches was largely unaffected by the inclusion of an additional subgroup when K = 50.

When K = 100 and S = 2, mean MCC values for subgrouping multi-VAR similarly increased as a function of time series length, with values of 0.82 when T = 50 and 0.88 for T = 100 when subgroups were balanced, and values of 0.79 and 0.82 when unbalanced. The addition of a third balanced subgroup resulted in similar model recovery metrics for subgrouping multi-VAR, with mean MCC values 0.81 and 0.87. Notably, model recovery was better for three unbalanced subgroups than two unbalanced subgroups, with mean MCC values of 0.81 when T = 50 and 0.86 when T = 100. In general, observed results suggest that increasing K from 50 to 100 resulted in a slight improvement in model recovery. Discrepancies between data-driven and confirmatory approaches remained minimal. Model recovery for S-GIMME mirrored results observed in previously reported conditions, with mean MCC values of 0.66 and 0.80 when T = 50 and T = 100, respectively, for each combination of number of subgroups and membership composition. Conversely, model recovery for scGVAR improved when K = 100, with mean MCC values ranging from 0.58 to 0.59 when T = 50 and from 0.66 to 0.67 when T = 100. Standard multi-VAR displayed a slight decrease in model recovery when K = 100, with mean MCC values of approximately 0.67 for T = 50 and 0. 75 for T = 100 across all subgroup number and membership composition conditions.





  
Parameter Estimates

Mean absolute bias and RMSE values can be found in Table <ref> and are visualized in Figure <ref>. Mean absolute bias for subgrouping multi-VAR was largely consistent across all conditions examined in the current study, with values ranging from 0.018 to 0.023. In general, bias was at a minimum when K = 50 and T = 100, mirroring results observed for model recovery. Conversely, bias was largest when K = 100 and T = 50, suggesting that an increase in sample size without a corresponding increase in time series length results in reduced performance. Thus, large sample sizes may result in a high degree of heterogeneity, thereby diminishing model performance when the number of time points is small [e.g.,][]Fisher2024. It is worth noting, however, that these differences were small in magnitude, and may be a function of sampling variability, consistent with both observed MCE estimates and visualization of estimates in Figure <ref>. There were no meaningful differences observed between data-driven and confirmatory approaches with respect to mean absolute bias. Differences were observed, however, between subgrouping and standard multi-VAR, such that failing to account for the presence of subgroups increased mean absolute bias across all conditions, with values ranging from 0.021 to 0.032. Notably, compared to subgrouping multi-VAR, mean absolute bias for S-GIMME estimates was both smaller and less variable across all conditions, with values ranging from 0.014 to 0.018. Mean absolute bias for S-GIMME was smallest when K = 100 and T = 100, and largest when K = 50 and T = 100. Conversely, mean absolute bias for scGVAR estimates was larger than that observed for subgrouping multi-VAR across all conditions, with values ranging from 0.030 to 0.037.

A similar pattern of results was observed for the RMSE of subgrouping multi-VAR estimates. Indeed, RMSE was at a minimum when K = 50 and T = 100, and at a maximum when K = 100 and T = 50, with values ranging from 0.087 to 0.109 across all conditions. Differences in RMSE values between data-driven and confirmatory subgrouping approaches were less than 0.01 across all conditions. RMSE values for standard multi-VAR were larger across all conditions, with values ranging from 0.091 to 0.117. RMSE values for S-GIMME estimates were larger than those observed for subgrouping multi-VAR when both K and T were at a minimum, and smaller for all other conditions. Thus, the subgrouping multi-VAR estimates were less variable than S-GIMME estimates when both sample size and time series length were small. In general, subgrouping multi-VAR estimates were less variable than scGVAR estimates. However, RMSE values for scGVAR estimates were smaller than values for subgrouping multi-VAR estimates when both K and T were at a maximum.





  
Subgroup Recovery

ARI values for subgrouping multi-VAR, S-GIMME, and scGVAR are shown in Table <ref>. In general, subgroup recovery was considered good or excellent for subgrouping multi-VAR across all conditions Steinley2004. Subgroup recovery was at a maximum (ARI = 0.98) when K = 50, T = 100, and subgroup composition was balanced, indicating near perfect subgroup recovery. Subgroup recovery was at a minimum when K = 100, T = 50, with three unbalanced subgroups (ARI = 0.88). Conversely, subgroup recovery for S-GIMME was considered poor across all conditions, with ARI values ranging from 0.01 to 0.64. The minimum ARI value for S-GIMME was observed when K = 50, T = 50, with two unbalanced subgroups, and the maximum value occurred when K = 100, T = 100, with two balanced subgroups. Subgroup recovery for scGVAR ranged from moderate to excellent across all conditions, such that ARI values were at a minimum of 0.73 when K = 50, T = 50, with three unbalanced subgroups, and a maximum of 0.98 when K = 100, T = 100, with two subgroups. 

Mirroring observed ARI values, subgrouping multi-VAR generally recovered the correct number of subgroups. Observed differences in ARI values between subgrouping multi-VAR and S-GIMME were similarly consistent with discrepancies in the number of subgroups identified by both frameworks. When K = 100 and T = 50, for example, the number of subgroups recovered ranged from 12 to 17 for S-GIMME and three to six for subgrouping multi-VAR. Thus, though both approaches overestimated the number of subgroups in these conditions, S-GIMME was adversely impacted to a greater degree. This is consistent with prior work showing that recovery of nonzero dynamics in GIMME is reduced when the number of time points is small Nestler2021. Notably, scGVAR subgroup enumeration was not affected by time series length, suggesting that it may be particularly well-suited to settings in which number of time points is small. 



§ EMPIRICAL EXAMPLE


To demonstrate the utility of the subgrouping multi-VAR framework, we present an empirical example using data from Fisher et al. (). Data consisted of 40 individuals with a primary diagnosis of either major depressive disorder (MDD) or generalized anxiety disorder (GAD) who were assessed four times per day for 30 days. For the purposes of the current application, we restricted our analyses to the 10 variables related to MDD (e.g., down and depressed) and GAD (e.g., worried) symptomatology, thereby ensuring that the number of variables examined mirrored those assessed in the simulation study. All variables were measured using a visual analogue scale ranging from 0-100, such that 0 indicated that the participant did not experience the symptom at all during the preceding hours and 100 indicated that the symptom was experienced as much as possible. As the multi-VAR framework does not currently accommodate missing data, linear imputation via the  package imputeTS was employed. For more information regarding data characteristics and study procedures, see Fisher et al. ().

The motivation for the selection of the current empirical dataset was influenced by several factors. First, both the number of individuals and length of time series were consistent with data analyzed in the simulation study. The quality of the group-, subgroup-, and individual-level estimates were therefore unlikely to be impacted by characteristics of the dataset, such as sample size. Second, the symptomatology of many psychopathological syndromes, including MDD and GAD, is characterized by persistent heterogeneity [e.g.,][]Kotov2017. Indeed, prior work identified 1030 unique symptom profiles in a sample of individuals diagnosed with MDD Fried2015. Thus, the current example represents an evaluation of the degree to which subgrouping multi-VAR accommodates quantitative and qualitative heterogeneity in an empirical setting. Finally, all participants had a primary diagnosis of either MDD or GAD. Comparisons between diagnostic status obtained via structured clinical interview and subgroup membership derived by subgrouping multi-VAR were therefore possible. 

Estimation of group-, subgroup-, and individual-level effects proceeded in the manner detailed above. First, standard multi-VAR without subgrouping was fit to the data. Estimates of the individual-level effects, Υ̂^k, were then used to derive subgroup membership. Using these subgroup labels, subgrouping multi-VAR was then fit to the data. Paralleling the modeling procedure for the simulated data, the autoregressive effects were not subjected to penalization via the adaptive Lasso in either the subgroup enumeration or estimation stages. In addition to analysis of the estimated group-, subgroup-, and individual-level effect, results were examined with respect to several between-person variables. That is, associations between subgroup membership and constructs of interest, such as diagnostic status, were assessed. All analyses were conducted using the  R package multivar.



 §.§ Empirical Results


Estimated effects are depicted in Figure <ref>. Estimates of common effects consisted of all autoregressive paths, with magnitudes ranging from 0.23 to 0.31. This autoregressive behavior is consistent with findings from the original study Fisher2017. With respect to subgroup enumeration, subgrouping multi-VAR identified four subgroups. The first three subgroups consisted of 25, 10, and four individuals, respectively, and the fourth subgroup was defined by a single individual. Given its singleton status, the fourth subgroup was excluded from subsequent analyses. Estimated effects for the first subgroup consisted of all autoregressive paths, with magnitudes ranging from 0.17 to 0.26, and a bidirectional relationship between anhedonia and down and depressed. The second subgroup was also characterized by estimated autoregressive effects, though most of these were negative, indicating that participants in this subgroup had autoregressive dynamics that were smaller in magnitude than those estimated at the group level. A relation between down and depressed and anhedonia was also observed for participants in this subgroup, though it was unidirectional in nature. Thus, the first two subgroups were distinguished by differences in both the sign of the estimated autoregressive effects and the nature of the dynamic relation between down and depressed and anhedonia.

The third subgroup also featured estimates of all autoregressive effects, with magnitudes ranging from -0.01 to 0.08, suggesting that, in general, participants in this subgroup were characterized by slightly stronger autoregressive dynamics than those derived at the group level. Notably, differences between this subgroup and the first two were observed with respect to the number of estimated cross-lagged effects. Indeed, out of 90 possible cross-lagged paths, 19 were estimated as nonzero, suggesting that the network of dynamic relations between MDD and GAD symptoms was much denser for individuals in this subgroup. Despite the observed qualitative and quantitative heterogeneity in estimated subgroup effects, subgroup membership was not associated with diagnostic status or between-person measures of depression and anxiety. However, given the degree of comorbidity in the present sample—75% of individuals had at least one comorbid diagnosis—identifying clear distinctions between participants with respect to diagnostic status and symptom level is a challenging endeavour.

Estimates of individual-level effects varied both qualitatively and quantitatively between individuals. For some individuals, these effects corresponded to the autoregressive paths, thereby serving to increase or decrease the magnitude of these dynamics, depending on the sign, compared to group- and subgroup-level estimates. Estimated effects for other individuals epitomized the persistent qualitative heterogeneity inherent in many dynamic processes of interest, such that nonzero cross-lagged paths represented an increase in both the complexity and density of symptom relations. Notably, estimates of individual-level effects for one individual were all zero, suggesting that their pattern of dynamics was fully captured by the group- and subgroup-level estimates. 





§ DISCUSSION


Many dynamic processes are characterized by persistent heterogeneity in both the magnitude of the effects of interest and the functional form of the process itself. Failing to appropriately account for such heterogeneity when present limits the degree to which observed results can be generalized across various levels of analysis Molenaar2004. The current paper sought to address this challenge through the introduction of a subgrouping extension to the multi-VAR framework. Built on the VAR model, subgrouping multi-VAR is a method for analyzing multiple-subject, multivariate time series that incorporates both data-driven subgroup identification and penalized estimation of group-, subgroup-, and individual-level effects. In contrast to similar approaches, the incorporation of penalty parameters at each level of analysis allows for both sparsity in the estimated transition matrices and flexible approximation of potentially heterogeneous data generating processes. In doing so, it helps to bridge the divide between nomothetic and idiographic approaches to the analysis of ILD. The efficacy and utility of the subgrouping multi-VAR approach was demonstrated in a simulation study and empirical example.

The simulation study assessed the performance of the subgrouping multi-VAR framework across a range of design factors. Moreover, the current approach was compared to two alternative multiple-subject, multivariate time series methods, as well as benchmark comparisons with the standard multi-VAR framework and multi-VAR with confirmatory subgrouping (i.e., when true subgroup membership was known). With respect to model recovery, subgrouping multi-VAR outperformed S-GIMME, scGVAR, and the standard multi-VAR approach across all conditions, and demonstrated only a slight decrease in performance compared to multi-VAR with confirmatory subgrouping. Model recovery was strongest for the condition corresponding to 50 individuals, 100 time points, and three balanced subgroups, and weakest for 50 individuals, 50 time points, and three unbalanced subgroups. Even in this weakest condition, however, subgrouping multi-VAR demonstrated good recovery of the true data generating process. Observed differences in model recovery between subgrouping multi-VAR and comparison methods are likely a function of both the unique manner in which the current approach models heterogeneity in dynamic processes and the quality of the estimates used in the subgroup enumeration stage. Indeed, the structured nature of the penalization procedure and simultaneous estimation of effects of interest ensures that the heterogeneity of the solution is determined by the competition of the three penalty parameters corresponding to the group-, subgroup-, and individual-level effects. Conversely, estimation of heterogeneous dynamics in the S-GIMME framework proceeds iteratively and is a function of user-specified thresholds, modification indices, and fit statistics. The scGVAR approach is similarly iterative in nature, such that effects at different levels of analysis are estimated separately. Though prior work has shown that both S-GIMME and scGVAR demonstrate good model recovery (e.g., <cit.>; <cit.>), the subgrouping multi-VAR estimation procedure described herein may better approximate data generating processes characterized by a high degree of quantitative and qualitative heterogeneity, consistent with observed results. 

Notably, this was also observed in the subgroup enumeration stage, such that subgrouping multi-VAR exhibited excellent identification of subgroup membership across nearly all conditions. Thus, standard multi-VAR estimates used in the construction of the Walktrap adjacency matrix Pons2006 can effectively inform the derivation of shared patterns of dynamics. These results contrast those observed for S-GIMME, wherein subgroup recovery was considered poor across all conditions. Indeed, S-GIMME often overestimated the number of subgroups, whereas subgrouping multi-VAR and scGVAR more consistently recovered the correct subgroup structure. These results therefore suggest that subgroup enumeration in the S-GIMME framework, wherein group-stage estimates are used to derive subgroup membership, is less effective than approaches that employ individual-level estimates, such as subgrouping multi-VAR and scGVAR. It is worth noting, however, that observed differences in subgroup identification between subgrouping multi-VAR and comparison methods may not generalize to situations in which subgroup-level effects are driven by contemporaneous dynamics, as the multi-VAR framework does not accommodate such dynamics. Prior S-GIMME work, for example, observed subgroup-level effects that consisted entirely of directed contemporaneous relations Lane2019. Future work should examine the behavior of subgrouping multi-VAR under such conditions.

Mirroring model recovery results, absolute bias and RMSE for subgrouping multi-VAR estimates were lowest for conditions corresponding to 50 individuals, 100 time points, and three subgroups. Largest values of absolute bias and RMSE were observed for conditions with 100 individuals and 50 time points, suggesting that analyses characterized by a large number of individuals may require a correspondingly large number of time points. In general, however, the quality and variability of estimated parameters—as indicated by absolute bias and RMSE—remained relatively stable across conditions. More work is needed to determine the behavior of these metrics as the discrepancy between number of individuals and number of time points increases. In contrast to model recovery results, absolute bias and RMSE values for subgrouping multi-VAR estimates were not consistently preferable when compared those obtained for S-GIMME estimates. Across all conditions, for example, subgrouping multi-VAR parameter estimates exhibited a greater degree of bias than S-GIMME estimates, though these differences were small in magnitude. Similar results were observed for RMSE values when data consisted of either 100 individuals or 100 time points. However, when both number of time points and number of individuals were small, multi-VAR parameter estimates were less variable than those obtained from S-GIMME. Observed differences in the quality and variability of parameter estimates between subgrouping multi-VAR and S-GIMME are generally consistent with prior work showing that GIMME estimates exhibit lower absolute bias and RMSE values than standard multi-VAR when the data are highly heterogeneous Fisher2024. Differences in absolute bias and RMSE values between subgrouping multi-VAR and scGVAR were more consistent with model recovery results. That is, subgrouping multi-VAR estimates were less biased and variable than scGVAR estimates across most conditions. 

Results from an empirical example using data from Fisher et al. () demonstrated the utility of the subgrouping multi-VAR framework. At the group level, symptom relations were characterized entirely by autoregressive dynamics, and were consistent with results observed in other clinical samples (e.g., <cit.>; <cit.>). Symptoms of MDD and GAD therefore exhibited some degree of inertia [e.g.,][]Kuppens2010 for all individuals. In addition to group-level effects, three primary subgroups were identified, each distinguished by differences in both the magnitude and structure of symptom dynamics. Though subgroup membership was not associated with between-person variables of interest, such as diagnostic status, identification of shared patterns of symptom relations is not without utility. Indeed, observed findings are consistent with prior work showing that meaningful patterns of subgroup-specific dynamics are not necessarily associated with mean levels of symptoms Lane2019. Such dynamics could, for example, inform future research and treatment efforts, consistent with the recent shift in focus toward the development of personalized models of psychopathology Wright2020. The degree of quantitative and qualitative heterogeneity observed with respect to estimated individual-level effects similarly indicates a need for person-specific approaches to prevention and intervention [e.g.,][]Fisher2019.

Simulation and empirical results observed in the current study should be interpreted in the context of several limiting factors, each of which represents an area for future research. First, as noted above, the multi-VAR framework cannot estimate contemporaneous dynamics. As such, subgrouping multi-VAR may not be an appropriate methodological approach for processes characterized primarily by said dynamics. Next, despite work showing that Walktrap is a reliable community detection algorithm Gates2016, the efficacy of alternative clustering approaches was not assessed. Recent work, for example, suggests that Walktrap could be improved by replacing Ward's hierarchical clustering algorithm Ward1963 with methods for K-means clustering Brusco2024. Moreover, approaches relying on modularity optimization may fail to identify subgroups consisting of a small number of individuals Fortunato2007. More work is needed to determine if alternative approaches improve the performance of subgrouping multi-VAR. Finally, observed results may not generalize to conditions not assessed in the current simulation and empirical studies. Indeed, prior work has demonstrated that model recovery for standard multi-VAR is impacted by the number of variables comprising the multivariate time series, such that performance improves as number of variables increases Fisher2022. 

The current study introduced a novel methodological framework for analyzing multiple-subject, multivariate time series characterized by persistent quantitative and qualitative heterogeneity. Results from both an extensive simulation study and empirical example suggest that subgrouping multi-VAR is an effective approach for estimating group-, subgroup-, and individual-level dynamics. Notably, the current approach demonstrated good model recovery under commonly encountered conditions, such as when both sample size and time series length are small. Moreover, the advantage of subgrouping multi-VAR over popular alternatives with respect to model recovery was greatest when time series length was small, suggesting that the current approach may be particularly well suited for the types of data frequently collected in the social, behavioral, and health sciences.









