Although LTA accounts for the collection of data from the same individuals across time (i.e., time nested within person), the model can also be extended to account for individuals being nested within higher level units, such as schools, hospitals, organizations, etc. In education research, statistical methods are commonly used that model students nested within schools, which is the context for the illustration in this paper. In such cases, the hierarchical structure of organizational data must be taken into account because independence between observations is not tenable; that is, students (i.e., level-1 units) are nested within and share influence of schools (i.e., level-2 units). Thus, multilevel LTA can be used to estimate the number of latent classes in each structured unit and the transitions participants make between classes across time. Finally, the transitions/stability between latent classes across time can be treated as the outcome in and of itself, or the transitions/stability can be used as a correlate or predictor of some other, distal outcome.

The multilevel LTA can be expressed as a series of multinomial logistic regressions at level-1 and as a linear regression at level-2 (Asparouhov and Muthén, 2008). To illustrate, consider the model below that has two latent classes across two time points. At level-1 the multinomial logistic regression for the latent classification variable at time 1, C 1 = 1,2 , can be expressed as P ( C 1 i g = 1 ) = exp ( α 1 g ) exp ( α 1 g ) + 1 , (1) where C 1 i g represents the latent class at time 1 for individual i in group g, and α 1 g is the intercept of latent class 1 for group g and is assumed to be normally distributed. The intercept for latent class 2 at time 1 is set to zero for identification because only one intercept is needed to distinguish two latent classes. The multinomial logistic regression of latent class at time 2 ( C 2 ) on latent class at time 1 ( C 1 ) can be expressed as P ( C 2 i g = 1 | C 1 i g = c ) = exp [ α 2 g + γ I ( C 1 i g = 1 ) ] exp [ α 2 g + γ I ( C 1 i g = 1 ) ] + 1 , (2) where α 2 g represents the latent class intercept at time 2 for individual i in group g, and γ represents the expected change in logits from the multinomial logistic regression predicting latent class at time 2 from latent class at time 1. The indicator function ( I ( C 1 i g = 1 ) ) in Eq. 2 demonstrates how the latent regression parameter γ is specific to class 1 which is how class specific transition probabilities are captured in the model (Asparouhov and Muthén, 2008; Kaplan et al., 2011). That is, the indicator function takes on values { 0,1 } depending on whether or not latent class membership at time 1 was equal to 1. Like the latent class intercept for class 2 at time 1, the latent class intercept for class 2 at time 2 is set to zero for identification. It should be noted the number of regression weights, γ, increases as the number of latent classes increases. For example, two latent classes implies one γ whereas three classes implies up to four γs. Additionally, the model above can be extended to incorporate student level covariates into the transitional structure (Vermunt et al., 1999).

The multinomial autoregression model above is akin to what is used in single-level LTA; however, a unique contribution of multilevel LTA is the incorporation of a latent regression model of latent class intercepts over time. At level-2, the random effects of latent class size across schools can be explained as part of a series of latent linear regression models such as α 1 g = μ α 1 + ε α 1 (3) α 2 g = μ α 2 + β α 1 g + ε α 2 . (4) The regression in Eq. 3 models the difference in latent class size across schools at time 1 where α 1 g is the latent class size in logits for school g at time 1, μ α 1 is the average latent class size at time 1, and ε α 1 is random effect of latent class size across schools. Similarly, the regression in Eq. 4 models the differences in latent class size across schools at time 2, where α 2 g is the latent class size in logits for school g at time 2, μ α 2 is the average latent class size at time 2 unconditional on latent class size at time 1, β is the fixed effect of latent class size at time 1 on time 2 latent class size, and ε α 2 is the random effect of latent class size across schools unique to time 2. The random effects are commonly assumed to be normally distributed with unique variance estimates at each timepoint [e.g., var ( α 1 ) and var ( α 2 ) ]. The level-2 latent regression of the multilevel LTA model expresses how latent class sizes change over time among schools. Factors that influence differences in latent class size over time among schools can be studied in more detail if level-2 covariates are included in the model. The incorporation of covariates can be guided by substantive interest and by information about how much information can be accounted for by these covariates. The amount of information that is contained in the level-2 portion of the model can be expressed by different R 2 -like measure that can be computed.

The multilevel LTA model has many parameters to describe the process which generated the differences in observed characteristics across time, and some of the model features are directly interpretable whereas other features are less easily interpreted. In order to help explain the complex features of the model, various R 2 -like measures can be computed to provide information about how the variability in latent class membership is influenced by 1) time, 2) nested data structure, and/or 3) individual latent class membership. For example, Asparouhov and Muthén (2008) and Kaplan et al. (2011) explicitly described an R 2 measure for the proportion of variance in latent class membership at time 2 that is accounted for by latent class membership at time 1. And is readily available for use in single-level LTA as well. This R 2 measure is R 2 = γ 2 P ( C 1 = 1 ) ( 1 − P ( C 1 = 1 ) ) γ 2 P ( C 1 = 1 ) ( 1 − P ( C 1 = 1 ) ) + π 2 / 3 , (5) where P ( C 1 = 1 ) is the probability of being latent class 1 at time 1 and π 2 / 3 is the residual variance associated with the logistic regression performed at level-1. P ( C 1 = 1 ) may also be viewed as the relative size of latent class 1 at time 1. Although the result in Eq. 5 is useful, there is more information in a multilevel LTA model that can be used to gain additional insights into the process under investigation.

Asparouhov and Muthén (2008) used these other R 2 -like measures, such as the proportion of variance in C 1 explained by α 1 , the proportion of variance in C 2 explained by the group effect at time 1, among others, but they did not describe the steps necessary to calculate these values. A detailed explanation on how to obtain such R 2 -like measures is therefore one of the major contributions of this work.

One potential limitation of the R 2 measure in Eq. 5 is that the hierarchical structure of the data is ignored, which means that it may overestimate the effect latent class membership at time 1 has on latent class membership at time 2. The methods we demonstrate for decomposing the variance in multilevel LTA explicitly account for this feature of the data. That is, the variance decomposition we describe accounts for the nested data structure by incorporating all model components into the variability in latent class membership at time 2.

There are several considerations specific to multilevel LTA that extend beyond those associated with LCA and LTA. First, one must consider whether the research question posed requires the multilevel aspect of the data explicitly incorporated into a mutlilevel LTA model. Not all questions require that the nested nature of the data be explicitly modeled (McNeish et al., 2017). For example, a researcher primarily interested in transitions of students among latent classes over time may not need to explicitly account for a school effect if differences among schools does not influence the students’ transitions. Instead, the multilevel aspect of the data can be incorporated implicitly through the use of sampling weights (Stapleton, 2013) or alternatives such as cluster-robust standard errors (McNeish et al., 2017). However, the use of multilevel LTA is likely warranted when researchers believe that characteristics of the group or school are related to differences in latent class membership. This is commonly encountered in education and healthcare applications where between-school and between-hospital differences, respectively, influence large groups of participants simultaneously.

In additional to the nested feature of one´s data, another important consideration is the time scale in which data were collected. The time scale of data collection may, or may not, adhere time scale of the transitions that individual may experience. Collins and Lanza (2009, p. 209–211) expressed how the transition structure estimated may reveal only chance transitions due to a underlying structure that transitions very rapidly (e.g., the example of indicators of depression in the last week but data were collected one year apart). Therefore, researchers must think carefully about how observed transitions among latent classes are related to transitions in the underlying construct of interest. In multilevel LTA, in particular, an additional consideration is whether the time scale of the transition is equal across level-2 units, such as schools. In healthcare settings, for example, the time scale of transitioning among depression latent classes may depend in part on the care received across different clinics if clinics were to have a general approach to helping patients with, say, depressive symptoms. As noted above, these considerations should be applied in addition to those important considerations that have been identified for LCA and LTA, such as model selection (Nylund et al., 2007; Tofighi and Enders, 2008; Morgan, 2015), label switching (Chung et al., 2004; Tueller et al., 2011), nature of the latent variables (Lubke and Neale, 2008), and incorporation of distal outcome (Lanza et al., 2013; Bakk and Vermunt, 2016; Nylund-Gibson et al., 2019). An excellent collection of applied and methodological papers using these procedures can be found on the Mplus website (www.statmodel.com/paper.shtml).

Next, we illustrate the use of multilevel LTA and explicitly model the multilevel nature of the data.

The data used are a subset of the ECSL-K national dataset (Tourangeau et al., 2009). The analytic sample for this demonstration was approximately 7,080 students nested within approximately 1,100 schools (sample sizes have been rounded to the nearest 10 in compliance with federal restricted-use data reporting guidelines). Prior to estimating the multilevel LTA model, we subset the ECLS-K data file on the students who remained in the same school from at least Grade 3 to Grade 5. The average number of students per school was 6.4 ( S D = 5.3 ) and ranged from 1 to about 30 students.

In order to demonstrate the model output and subsequent decomposition of the model variance, we used the five Social Rating Scale subscales from the Early Childhood Longitudinal Survey–Kindergarten (ECLS-K) data. The five major constructs of interests are: Approaches to Learning (AtL), Self-Control (SC), Interpersonal Skills (IPS), Externalizing Problem Behaviors (EPB), and Internalizing Problem Behaviors (IPB) (Tourangeau et al., 2009). These five constructs of child behaviors/characteristics are modeled as being reflective of a child’s need for possible additional behavioral intervention. The reliability estimates (coefficient α) for these constructs in the full ECLS-K in spring of fifth grade ranged from 0.77 (Internalizing Problem Behaviors) to 0.91 (Approaches to Learning). Reading teachers were asked to report how frequently students exhibited the social skill or behavior identified by each item. The response scale used a four-point frequency scale ranging from 1 (Never) to 4 (Very Often). The same 26 SRS items administered in Grade 3 and 5. A summary of these raw subscale scores is shown in Table 1.

The raw subscale scores were computed as the average of the responses to the items on each subscale.

The model was estimated using maximum likelihood estimator with robust standard errors (MLR) in Mplus v8.4 (L. Muthén and Muthén, 2017) using 2,000 random starting values and 50 final stage optimizations. For illustrative purposes, we estimated only a two-class solution. In practice, additional class enumeration models would be estimated and compared. For this demonstration, we elected to not use sampling weights to reduce the complexity of the example analysis. All inferences from the following model are restricted to this sample of students and is not necessarily a representation of the characteristics of students more broadly.

The path diagram for the multilevel LTA model is presented in Figure 1. We should note that the path diagram includes variance components to aid in interpretation of variance decomposition discussion below.

The major inferential goals are the evaluation of the transition parameters ( γ ,   β ) and the variability in latent class size across schools ( var ( α 1 ) ,   var ( α 2 ) ).

The resulting latent class patterns are shown in Table 2. In the estimation, the latent class structure was fixed to be invariant across time. Latent class 1 is characterized by students who had lower ratings on the three positive constructs (i.e., AtL, SC, and IPS) and higher scores on the constructs reflecting problem behaviors (i.e., EPB and IPB). Latent class 2 was characterized as having higher scores on the three positive constructs (i.e., AtL, SC, and IPS) and lower ratings on the problem behavior constructs (i.e., EPB and IPB).

The structural model parameters are described in Table 3. At Time 1, Class 1 was the smaller of the two latent classes, making up about 32% of the sample, whereas Class 2 made up about 68% of the sample. Due to the multilevel nature of the data, the parameter estimate, var ( α 1 ) = 0.64 , offers additional insights into the latent class structure at time 1. That is, the estimate of 0.64 suggests the proportion of students in Class 1 and Class 2 at time 1 varies depending on the school. In other words, Class 1 contains about 32% of the students at time 1, on average, but this percentage differs across schools with a 95% probable range of 9–69%. The larger the variance estimate, the greater the school effect and greater range of relative class sizes across schools.

The transition component of the multilevel LTA model is characterized by the parameters γ ( γ = 2.93 , S E = 0.11 , p < 0.001 ) and β ( β = − 0.19 , S E = 0.11 , p = 0.077 ). From these two parameters, a transition matrix ( τ ) is constructed to help explain the overall effect of time. The details of computing these values are given in the Multilevel LTA Variance Decomposition section; but for now, these results are reported in Table 4 along with the interpretation. We found that the, on average, about 13% of students who were classified in Class 2 at Time 1 (i.e., third grade) transitioned into the Class 1 at Time 2 (i.e., fifth grade). Of those students classified in the Class 1 at Time 1, approximately 26% transitioned into Class 2.

As alluded to above, there are numerous calculations necessary to extract important modeling results that guide interpretation. The intraclass correlation (ICC) estimate for this model indicates that about 16% proportion of variability in class assignment at Time 1 can be explained by the school effect. Using the variance decomposition, 10.4% of the variability in latent class membership at Time 2 can be accounted for by the school effect at Time 1. However,15.8% of the variability in latent class membership at Time 2 can be accounted for by the school effect at Time 2. The incorporation of school level covariates into Eq. 3 and Eq. 4 could give insight into what school characteristics are associated with latent class membership at Time 2 by investigating the change in the previous two R 2 measures in models estimated with and without the covariate included. At the student level, the class assignment at Time 1 accounted for about 30.3% of the variability in class assignment at Time 2. The unexplained variability in class assignment at Time 1 also accounted for about 53% of the variability in class assignment at Time 2. Overall, the multilevel LTA model without any covariates explained approximately 46.5% of the variance in class assignment at Time 2. Finally, the inclusion of the multilevel structure explained about 16.2% of the variability in latent class membership at Time 2.

To derive the variance components, the structural equations associated with the path diagram are needed. The structural equation are: α 1 g = μ α 1 + ε α 1 g (6) α 2 g = μ α 2 + β α 1 g + ε α 2 g (7) C 1 i g * = α 1 g + ε C 1 * i g (8) C 2 i g ∗ = α 2 g + γ I ( C 1 i g = 1 ) + ε C 2 ∗ i g . (9) The variances associated with these structural component are defined as follows. The variance of α 1 g , the random effect at time 1, reduced the variance of the error term only, as μ α 1 is a constant. V ( α 1 g ) = V ( μ α 1 + ε α 1 g ) = σ α 1 2 . (10) The remaining pieces are slightly more complex. For the variance of the random effect at time 2, a long form derivation is V ( α 2 g ) = V ( μ α 2 + β α 1 g + ε α 2 g ) = V ( β α 1 g ) + V ( ε α 2 g ) + 2 ℂ O V ( β α 1 g , ε α 2 g ) = β 2 V ( α 1 g ) + V ( ε α 2 g ) = β 2 σ α 1 2 + σ α 2 2 . (11) It should be noted that we assumed that the covariance between the time 1 random effect and the time 2 random effect is 0.

The variance of the latent response tendency variable relative to the reference class 2 is defined as follows. V ( C 1 i g * ) = V ( α 1 g + ε C 1 i g * ) = V ( α 1 g ) + V ( C 1 i g * ) = σ α 1 2 + π 2 3 . (12) Again, we assumed that the error terms between the random effect at level 2 and the latent response tendency variable residual variance for the logistic regression have a covariance of 0. The residual variance of the latent response residual variance is a known constant of π 2 3 ≈ 3.29 .

Lastly, the variance of the latent response tendency variable for time 2 is defined as V ( C 2 i g ∗ ) = V ( α 2 g + γ I ( C 1 i g = 1 ) + ε C 2 i g ∗ ) = V ( α 2 g ) + V ( γ I ( C 1 i g = 1 ) ) + V ( ε C 2 i g ∗ ) = β 2 σ α 1 2 + σ α 2 2 + γ 2 Pr ( C 1 i g = 1 ) ( 1 − Pr ( C 1 i g = 1 ) ) + π 2 3 . (13) Again, the assumption of a covariance of 0 among the terms is imposed. The unique part of obtaining the variance of the latent response variable at time 2 is that an indicator function is a part of the structural equation. An indicator function, I ( . ) , is Bernoulli random variable with variance of Pr ( c o n d i t i o n   t r u e ) × ( 1 − Pr ( c o n d i t i o n   t r u e ) ) . Therefore, the variance of the indicator function in this case is a function of the size of class 1 (i.e., V ( I ( C 1 i g = 1 ) ) = Pr ( C 1 i g = 1 ) × ( 1 − Pr ( C 1 i g = 1 ) ) ).

To summarize, the variance components are V ( α 1 g ) = σ α 1 2 V ( α 2 g ) = β 2 σ α 1 2 + σ α 2 2 V ( C 1 i g ∗ ) = σ α 1 2 + π 2 3 V ( C 2 i g ∗ ) = β 2 σ α 1 2 + σ α 2 2 + γ 2 Pr ( C 1 i g = 1 ) ( 1 − Pr ( C 1 i g = 1 ) ) + π 2 3 .

The R 2 -like measures that we can compute to help interpret the results from ML-LTA can therefore be defined as follows. First, a useful initial measure is the intraclass correlation, defined at time 1 as I C C = σ α 1 2 σ α 1 2 + π 2 3 . (14) The ICC above will be a useful component to disentangle the variance of the latent response tendency variable at time 1. The R 2 -like measures are as follows.

The estimate of the proportion of variance in latent class membership at time 2 ( C 2 * , the latent response tendency on logit scale) explained by the random effect at time 1 is R 1 2 = V ( α 1 g ) V ( C 2 i g * ) . (15) The estimate of the proportion of variance in latent class membership at time 2 ( C 2 * ) explained by residual variance of α 2 is R 2 2 = σ α 2 2 V ( C 2 i g * ) . (16) The estimate of the proportion of variance in ( C 2 * ) explained by residual of ( C 1 * ) is R 3 2 = ( 1 − I C C ) × V ( C 1 i g * ) V ( C 2 i g * ) . (17) The estimate of the proportion of variance in C 2 * explained by C 1 * is R C 1 2 = γ 2 Pr ( C 1 i g = 1 ) ( 1 − Pr ( C 1 i g = 1 ) ) V ( C 2 i g ∗ ) . (18) The proportion of variance in C 2 * explained by the model is the combination of all the variance components in the denominator minus the residual variance, that is R model 2 = β 2 σ α 1 2 + σ α 2 2 + γ 2 Pr ( C 1 i g = 1 ) ( 1 − Pr ( C 1 i g = 1 ) ) β 2 σ α 1 2 + σ α 2 2 + γ 2 Pr ( C 1 i g = 1 ) ( 1 − Pr ( C 1 i g = 1 ) ) + π 2 3 . (19) Lastly, the proportion of variance in C 2 * explained by adding the level-2 structure can be estimated as R a d d 2 = β 2 σ α 1 2 + σ α 2 2 β 2 σ α 1 2 + σ α 2 2 + γ 2 Pr ( C 1 i g = 1 ) ( 1 − Pr ( C 1 i g = 1 ) ) + π 2 3 . (20) It should be noted that similar decomposition is possible for higher number of latent classes at each time point. However, the decomposition is more involved given the complexity of more transitions and random effects at level-2. Methods for expanding the results described above to k-class solutions are built on ideas similar to the random effects models for multinomial outcomes (Hedeker, 2003). We are currently developing the extension to three latent classes and intend to identify some concise patterns that will allow for relatively straightforward variance decomposition with more latent classes.