Implications of Alternative Parameterizations in Structural Equation Models for Longitudinal Categorical Variables

Abstract

When analyzing scaling conditions in latent variable structural equation models (SEMs) with continuous observed variables, analysts scaling a latent variable typically set the factor loading of one indicator to one and either set its intercept to zero or the mean of its latent variable to zero. When binary and ordinal observed variables are part of SEMs, the identification and scaling choices are more varied and multifaceted. Longitudinal data further complicate this. In SEM software, such as lavaan and Mplus, fixing the underlying variables’ variances or the error variances to one are two primary scaling conventions. As demonstrated in this paper, choosing between these constraints can significantly impact longitudinal analysis, affecting model fit, degrees of freedom, and assumptions about the dynamic process and error structure. We explore alternative parameterizations and conditions of model equivalence with categorical repeated measures.

Using data from the National Longitudinal Survey of Youth 1997, we empirically explore how different parameterizations lead to varying conclusions in longitudinal categorical analysis. More specifically, we provide insights into the specifications of the autoregressive latent trajectory model and its special cases—the linear growth curve and first-order autoregressive models—for categorical repeated measures. These findings have broader implications for a wide range of longitudinal models.

1 Introduction

Researchers increasingly use categorical endogenous observed variables in their structural equation models (SEMs). Binary and ordinal variables are among the most common. Examples include binary variables registering life events such as divorce, job loss, or pregnancy or ordinal variables such as general health status (e.g., poor, fair, good, very good, excellent). Although SEMs with and without such categorical endogenous variables share many similarities, there are differences. Identification issues, which are absent with continuous endogenous variables, arise with their categorical counterparts when assuming a continuous latent variable underlies the observed categorical response. This issue is particularly pronounced in the context of longitudinal data analysis, where the choice of constraints becomes pivotal. The common diagonally Weighted Least Squares with Mean and Variance adjustments (WLSMV) estimator, such as that in Mplus 8.6 (Muthén & Muthén, 1998–2017) and lavaan 0.6-16 (Rosseel, 2012), ¹ offers various techniques for estimating models with categorical data. Many textbooks and papers recommend imposing constraints on thresholds, errors and/or underlying variable variances to establish model identification when using WLSMV. This advice implies that researchers can apply these constraints almost automatically without a clear preference for one over the other. Our paper demonstrates that such constraints are not universally applicable, especially with longitudinal data. Different constraints can lead to alternative model specifications that are not always equivalent (Lee & Hershberger, 1990; Levy & Hancock, 2007; Raykov & Penev, 1999; Stelzl, 1986).

Furthermore, relying on default parameterizations in SEM software that set the variance of the normally distributed variable underlying the observed categorical response to one is not always appropriate. For instance, Grimm and Liu (2016) argue that this is inappropriate with linear growth curve models, and they suggest setting the error variances to one. While researchers agree that the choice of setting the variances of the underlying variable or its error variance to one does not impact the fit of factor models in cross-sectional applications (Kamata & Bauer, 2008; Muthén & Asparouhov, 2002; Paek et al., 2018; Wang et al., 2023), this interchangeability need not carry over to repeated measures. Several papers empirically compare the two specifications for linear latent growth models but arrive at contrasting conclusions. For Grimm and Liu (2016) and Lee et al. (2018), parameterization affects both model fit and parameter estimation because different assumptions about the dynamic process and its errors are embedded in the chosen specification. On the other hand, for Newsom and Smith (2020), the choice between these two parameterizations is otherwise arbitrary as they are equivalent, and one solution can be transformed into another.

Various authors have proposed alternative parameterizations for latent variable SEMs for longitudinal data, mainly focusing on linear growth models. Muthén and Asparouhov (2002) and Grimm and Liu (2016) propose to allow the means and variances of the underlying variables to be freely estimated except on the first occasion. Thresholds are kept invariant over time, but for a multi-group confirmatory factor model, Millsap and Yun-Tein (2004) relaxed this assumption for all but two thresholds. Mehta et al. (2004) propose a different parameterization for the linear growth model applied to categorical data. It consists of freely estimating all the latent response means and variances and fixing the first two thresholds to zero and one, respectively, keeping the remaining thresholds invariant over time. Prior to their work, Jöreskog (2001) relaxed the threshold invariance condition in this parameterization.

Taken together, we find conflicting advice on the identifying constraints to impose when analyzing longitudinal data with endogenous categorical variables and incompatible claims on the equivalency of the constraints. The result is confusion and possible misapplication of categorical SEM techniques.

This paper examines the technical aspects of parameterization and modeling in analyzing binary and ordered categorical repeated measures. We provide a clear and easily accessible summary of the closed-form relations between longitudinal model specifications under alternative parameterizations. The role of these parameterizations is examined for the Autoregressive Latent Trajectory (ALT) model introduced by Bollen and Curran (2004) for continuous outcomes and its special cases, the linear latent growth curve and the autoregressive model. Unlike existing research, our study provides a comprehensive overview of possible parameterizations, evaluating necessary and sufficient conditions for model equivalence. Drawing on established comparisons, we determine when different models are empirically indistinguishable, shedding light on the practical implications of choosing specific constraints in the analysis of longitudinal categorical data.

2 Motivating data and example

The autoregressive latent trajectory model has proven invaluable for exploring longitudinal dynamics, especially with regard to changes in multiple variables over time. This is evidenced by its application to a diverse range of phenomena, including developmental trajectories of anxiety and depression (Connell et al., 2021; Lee & Vaillancourt, 2020; McLaughlin & King, 2015), trends in disposition and happiness (Caprara et al., 2017), and between aggressive behaviors and peer victimization in preadolescence (Yao & Enright, 2022). In the literature, it has been compared with other models that incorporate lagged effects alongside random components, such as the Random-Intercept Cross-Lagged Panel Model (RI-CLPM) (Hamaker et al., 2015) and the general cross-lagged panel model (Zyphur et al., 2020). A comprehensive discussion can be found in Usami (2021). A recent study by Andersen (2022) has also clarified and proven that the RI-CLPM is a constrained version of the ALT model.

While the model has been extensively applied to continuous variables, its application to categorical variables remains underexplored in the literature. This paper addresses this gap by presenting alternative specifications of the ALT model and two special cases—the linear growth curve and autoregressive of order one model—tailored for scenarios involving categorical variables.

To illustrate the practical implications of our study, we turn to a concrete example examining the co-occurrence of illegal drug use, depressive symptoms, and general health over time. Previous research by Silver et al. (2023) touched upon this nexus but applied default lavaan parameterization, treating illegal drug use as binary and depressive symptoms and general health as continuous variables. They separately analyzed autoregressive and cross-lagged effects and multivariate latent growth patterns for the three variables. Utilizing data from the National Longitudinal Survey of Youth 1997 (NLSY97), our study focuses on repeated measures of these variables as respondents transition from adolescence to adulthood, covering waves from 2000 to 2010. NLSY97 is a national study of 8984 respondents born in the US between January 1, 1980, and December 31, 1984. The respondents participated in in-person or phone interviews annually from 1997 (wave 1) to 2013 (wave 16), with two additional interviews in 2015 (wave 17) and 2017 (wave 18). For information on NLSY97 sampling and interviewing methods, refer to Moore et al. (2000).

Due to the lack of consensus on how to handle missing data with categorical variables, we followed the default option in lavaan 0.6-16 and applied listwise deletion, reducing the sample from 8,984 to 5,309 cases. ² Participants’ age ranged from 13 to 17 (average age: 14.92) during wave 4 and from 23 to 27 during wave 14. Data were analyzed for 2000 (wave 4), 2002 (wave 6), 2004 (wave 8), 2006 (wave 10), 2008 (wave 12), and 2010 (wave 14).

Illegal drug use was assessed by a single item asking respondents if they had used any illegal drugs (excluding marijuana and alcohol) since the previous interview. Responses were coded as one if the respondent had used illegal drugs and zero if he/she had not.

Depressive symptoms were evaluated based on respondents’ answer to the question: “How much of the time during the last month have you felt so down in the dumps that nothing could cheer you up?”. Ordinal responses were “all of the time,” “most of the time,” “some of the time,” and “none of the time,” with higher scores indicating more frequent depressive symptoms.

General health status was self-reported by respondents in answer to the question: “In general, how is your health?”. Choices were “excellent,” “very good,” “fair,” and “poor.” The measures of general health were coded with higher scores, indicating better perceived general health.

Table 1 presents the marginal proportions for each category of the observed variables at each time point. The frequency of adolescents experiencing depressive symptoms indicates a decrease in the likelihood of being more depressed as they age into adulthood. Similarly, the propensity to use illegal drugs and the perception of excellent health decreases over time, while the perception of good and very good health remains relatively constant, especially in the adulthood period (waves 10–14).

Table 1 Proportions for each year of the categories related to depressive symptoms, general health status, and illegal drug use

In contrast to conventional approaches, we explore various specifications of the ALT model, considering different constraints on thresholds, means of the underlying variables, and error variances. Our findings reveal that the choice of these identification constraints significantly influences the results, challenging the notion of interchangeable specifications.

Figure 1 illustrates different estimated ALT specifications, demonstrating the impact of different parameterizations on cross-lagged and autoregressive relationships (on the right) as well as growth components (on the left). Only statistically significant $(p-value<0.05)$ path coefficients are presented for simplicity. The two ALT components—the multivariate growth and the cross-lagged and autoregressive part—are shown separately for illustrative purposes, although they are estimated simultaneously to describe the temporal dynamic of the underlying latent variables.

Figure 1 Path diagram depicting the cross-lagged and autoregressive component (right) and multivariate growth part (left) of alternative autoregressive latent trajectory model specifications for illegal drug use (drug), depressive symptoms (depr), and general health status (health). The top panel is based on the standard (theta) parameterization of the auxiliary model linking observed variables to the underlying continuous ones, as adopted by Mplus 8.6 and lavaan 0.6-16. Dashed lines indicate statistically significant effects ( $p-value< 0.05$ ) only under this specification. The middle panel showcases an alternative parameterization proposed by Muthén and Asparouhov (2002). Paths in dark gray signify significance under the alternative parameterizations but not under the standard one. The bottom panel employs the parameterization introduced by Jöreskog (2001). Light gray paths indicate a significant pattern under this specification but not in the others.

The upper panel in Figure 1 adopts the standard (theta) parameterization of Mplus 8.6 and lavaan 0.6-16. In this configuration, all the error variances are set to one, all the thresholds in the model linking observed variables to the underlying continuous ones are freely estimated, and the mean of the underlying variables is fixed to zero.

Regarding autoregressive and cross-lagged relationships, prior illegal drug use, depressive symptoms, and previous health status perception consistently impact general health status on each occasion. However, only depressive symptoms and health status on the preceding occasion influence the propensity to experience depression on each occasion.

Concerning illegal drug use, no significant cross-lagged effects are observed; instead, there is a direct effect of the illegal drug use propensity on the preceding occasion on subsequent illegal drug use propensity in the initial three waves (4, 6, and 8 - corresponding to the years 2000, 2002, and 2004), associated with late adolescence/early adulthood. This contrasts with the middle panel in Figure 1, where a significant influence of general health status at waves 8 and 10 on illegal drug use propensities in waves 10 and 12 (associated with the years 2006 and 2008), respectively, is observed. This ALT specification uses an alternative parameterization proposed by Muthén and Asparouhov (2002) for categorical longitudinal data. In this approach, instead of fixing the error variances to one on all occasions, they are freely estimated on all occasions but one, with thresholds assumed to be time-invariant.

Conversely, in the lower panel of Figure 1, the propensities for illegal drug use during waves associated with late adolescence (6 and 8) are significantly influenced by depressive symptoms in the preceding occasion. This influence is specific to the late adolescence/early adulthood phase and not applicable during adulthood. The model parameterization in this lower panel aligns with the alternative approach proposed by Jöreskog (2001). It involves freely estimating the error variances and the means of the underlying variables on each occasion while fixing the first two thresholds to zero and one on all occasions.

The impact of alternative parameterizations of the ALT model is also noticeable in the estimated multivariate growth component. In the standard parameterization (top left panel), a correlation is observed between the intercept for the illegal drug use variable and the slope of general health status. Conversely, under the alternative parameterizations (middle and bottom panels), a correlation is estimated between the slope of illegal drug use propensity and the intercept of the growth component for health status perception. Notably, in the bottom panel, there is no significant covariance between the intercept and slope growth factors for general health status, distinguishing it from other specifications.

These findings highlight the critical importance of carefully selecting the most suitable parameterization for models involving categorical longitudinal data. The following sections offer a theoretical interpretation and practical application of these results.

3 The auxiliary measurement model linking $Y_{it}$ and $Y_{it}^{*}$

Let $Y_{it}$ be the ordered categorical measure for the ith person on the tth occasion, with $i = 1, 2, \ldots , n; t=1,2, \ldots , T$ . The values of $Y_{it}$ range from 0 to $C-1$ , where C is the number of response categories across all occasions.

In SEMs, a common approach is to consider the observed values of $Y_{it}$ as discretized manifestations of an underlying continuous variable, denoted as $Y_{it}^{*}$ . Depending on the measure, this underlying variable represents the level of understanding, attitude, or propensity to respond in a particular category. In our application on the NLSY data, it reflects the propensity to respond for all three variables under investigation. The observed category is then determined through

(1) $$ \begin{align}Y_{it}=c \quad \text{if} \quad \tau_{c,t} \leq Y_{it}^{*} < \tau_{c+1,t}, \quad c=0,1, \ldots, C-1.\end{align} $$

Here, $\{\tau _{0,t}, \tau _{1,t}, \ldots , \tau _{C,t}\}$ are threshold parameters for $Y_{it}^{*}, t=1, \ldots , T$ , where two of them are predefined, $\tau _{0,t}=-\infty $ and $\tau _{C,t}=\infty $ , whereas the remaining $(C-1)$ may vary across occasions. Let $\mathbf {Y}_{i}'=(Y_{i1}, Y_{i2}, \ldots , Y_{iT})'$ be a T-dimensional vector representing the observed scores for the ith person on the T occasions, and let $\mathbf {Y}_{i}^{*'}=(Y_{i1}^{*}, Y_{i2}^{*}, \ldots , Y_{iT}^{*})'$ be the corresponding vector for the underlying variables. The latter is typically assumed to follow a multivariate normal distribution with mean vector $\boldsymbol \mu _{\mathbf {Y}^{*}}$ and covariance matrix $\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ of size $T \times T$ , such that the probabilities associated with the observed values of $\mathbf {Y}_{i}$ can be determined by the probability distribution of $\mathbf {Y}_{i}^{*}$ .

The auxiliary model establishes the connection between the observed ordinal response $Y_{it}$ and the corresponding unobserved or latent continuous variable $Y_{it}^{*}$ at each time point. This linkage is achieved through unknown cut-points or thresholds. Although $Y_{it}^{*}$ is assumed to follow a normal distribution, its mean and variance remain unidentified due to the limited availability of ordinal information.

For an ordinal variable $Y_{it}$ with C categories measured at T fixed occasions, there are $C^{T}$ possible observed response patterns, not all of which may be observed in a given dataset. In this case, the sample data consists of the number of individuals with each of these response patterns. Any hypothesized model must explain (1) the univariate or marginal proportions, that is, the proportion of individuals in each of the C response categories for each of the T variables, and (b) the bivariate or joint proportions, that is, the proportions of individuals with each of the $C^{T}$ possible response pattern. The $T(C-1)$ observed marginal proportions are insufficient for estimating the thresholds $\tau _{ct}, c=1, \ldots , C-1, t=1,\ldots , T$ , or the parameters of the underlying variables, namely $\boldsymbol \mu _{\mathbf {Y}^{*}}$ and diag( $\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ ), without imposing certain restrictions.

3.1 Identification issues

Constraints must be applied to achieve model identification, and this is where different parameterizations of the auxiliary model come into play. Various sets of conditions exist that are sufficient to ensure identification, and they vary based on assumptions made about the thresholds (whether they are time-varying, time-invariant, or fixed) and the means and variances of the underlying variables (whether they are fixed on all occasions, all free but one, or freely estimated on all time points). Different software may adopt distinct identification conditions, leading to variations in parameter estimates and model fit across programs. Table 2 details the main parameterizations in the literature that we discuss.

Table 2 Alternative sets of identification constraints for the auxiliary model in presence of ordinal data

• Standard parameterization. The first row of Table 2 refers to the set of constraints commonly adopted in single-population or cross-sectional applications. It consists of fixing the expected value and standard deviation of $Y_{it}^{*}$ equal to zero and one, respectively, on each occasion. When this standard parameterization is adopted, the thresholds are freely estimated and assumed to be time-varying. These $T(C-1)$ parameters are determined as percentiles of the standard normal distribution, whereas the $T(T-1)/2$ off-diagonal elements of $\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ are estimated as polychoric correlations. This is the default parametrization applied in Mplus and lavaan, known as standard delta parameterization. Jöreskog (2001) also refers to it as standard parameterization, whereas Kamata and Bauer (2008) term it marginal parameterization since the marginal distribution of the continuous underlying variables is standardized.

One weakness of standard parameterization is that all underlying variables are standardized to have zero means and unit standard deviations. However, in the context of longitudinal data, where the response alternatives are the same across multiple time points, differences in the distribution of these variables can reflect differences in the means and/or variances of the underlying latent variables.

Researchers have proposed alternative identification constraints for categorical repeated measures, recognizing the limitations of the standard parameterization. These alternative parameterizations, which allow the means and variances of the underlying variables to vary freely, are a significant step toward a more accurate representation of the underlying data structure. Additional constraints are imposed, mainly on the thresholds of these variables, to ensure that the latent propensities are identified and comparable across time points. By imposing these constraints, analysts can establish a fixed reference point for the scale of the underlying response variables. This ensures that observed changes in the distributions reflect differences in the underlying variables, not inconsistencies in measurement.

• Alternative parameterization 1. The first alternative freely estimates the underlying variable’s means and variances on all occasions except on the first one (Millsap & Yun-Tein, 2004; Muthén & Asparouhov, 2002). That is, the underlying variable mean and variance at the first occasion, $\mu _{Y_{1}^{*}}^{alt1}$ and ${\sigma _{Y_{1}^{*}}^{2alt1}}$ , are fixed to zero and one, respectively, while $\mu _{Y_{t}^{*}}^{alt1}$ and $\sigma _{Y_{t}^{*}}^{2alt1}$ are estimated for $t> 1$ . ³ These assumptions are sufficient to identify all the thresholds on the first occasion. On subsequent occasions, two thresholds are assumed to be time-invariant, e.g., $\tau _{1t}^{alt1}=\tau _{1}^{alt1}$ and $\tau _{2t}^{alt1}=\tau _{2}^{alt1}$ for all t. A notable feature is that complete invariance of all threshold parameters is not required. As discussed before, to estimate differences in the means and variances of latent variables over time—as done with latent growth models—one must ascertain that the underlying variables are on the same scale on different occasions. This is achieved by defining a common across-time metric in terms of the standard deviation on the first occasion.

  As shown in the second row of Table 2, the number of parameters to be estimated is the same as in the standard parameterization, with $T(C-3)$ thresholds, $2(T-1)$ means and variances of the underlying variables, and $T(T-1)/2$ polychoric covariances, for a total of $T(C-1) + T(T-1)/2$ parameters. A one-to-one correspondence exists between the parameters in this alternative parameterization and those from the standard one. That is, on the first occasion,

  (2) $$ \begin{align} \begin{array}{ccc} \tau_{1}^{alt1}=\tau_{11}^{std},&\tau_{2}^{alt1}=\tau_{21}^{std},&\tau_{c1}^{alt1}=\tau_{c1}^{std}, \quad c=3, \ldots C,\\ \end{array}\end{align} $$

  whereas, on subsequent occasions,

  (3) $$ \begin{align} \begin{array}{ccc} \sigma_{Y_{t}^{*}}^{alt1}=\frac{\tau_{21}^{std}-\tau_{11}^{std}}{\tau_{2t}^{std}-\tau_{1t}^{std}},&\mu_{Y_{t}^{*}}^{alt1}=\frac{\tau_{11}^{std}\tau_{2t}^{std}-\tau_{1t}^{std}\tau_{21}^{std}}{\tau_{2t}^{std}-\tau_{1t}^{std}},&\tau_{ct}^{alt1}=\frac{\tau_{11}^{std}(\tau_{2t}^{std}-\tau_{ct}^{std})-\tau_{21}^{std}(\tau_{1t}^{std}-\tau_{ct}^{std})}{\tau_{2t}^{std}-\tau_{1t}^{std}}, \\&&\qquad \qquad\qquad\qquad \qquad\qquad c=3, \ldots C.\\ \end{array}\end{align} $$

  Some authors claim that analysts should also estimate the thresholds but constrain the same threshold to be equal over time (Muthén & Muthén, 1998–2017, Example 6.5). That is, $\tau _{ct}^{alt1}=\tau _{c}^{alt1}$ for $c=1, \ldots , C-1$ , and $t=1, \ldots , T$ , and this parameterization is detailed in the third row of Table 2. ⁴

  If the assumption of threshold invariance holds, then the alternative thresholds $\tau _{c}^{alt1}$ for each category, $c=1, \ldots , C-1$ , should be equal to the corresponding threshold estimated on the first occasion under the standard parameterization ( $\tau _{c1}^{std}$ ), since the mean $\mu _{Y_{1}^{*}}^{alt1}$ and standard deviation $\sigma _{Y_{1}^{*}}^{alt1}$ on the first occasion are set equal to zero and one, respectively. Based on this relationship,

  $$ \begin{align*}\tau_{ct}^{std}=\frac{\tau_{c1}^{std}-\mu_{Y_{t}^{*}}^{alt1}}{\sigma_{Y_{t}^{*}}^{alt1}},\end{align*} $$

  or, equivalently,

  (4) $$ \begin{align}\sigma_{Y_{t}^{*}}^{alt1}\tau_{ct}^{std}=\tau_{c1}^{std}-\mu_{Y_{t}^{*}}^{alt1}.\end{align} $$

  Based on eq. (4), if thresholds are assumed to be equal over time, the ratio of the standard deviation of the underlying variable on two different occasions must be equal to the ratio of the differences between any two thresholds in the standard parameterization on these two occasions:

  $$ \begin{align*}\frac{\sigma_{Y_{t}^{*}}^{alt1}}{\sigma_{Y_{t'}^{*}}^{alt1}}=\frac{\tau_{c't'}^{std}-\tau_{ct'}^{std}}{\tau_{c't}^{std}-\tau_{ct}^{std}},\quad c < c', t \neq t'.\end{align*} $$

  In other words, if thresholds under the alternative 1 parameterization are truly invariant at times t and $t'$ , the ratio of differences between any two standardized thresholds at the two time points should be equal. If any one of the thresholds for a time point is not invariant, then the equality mentioned above should not hold for the thresholds in the standard parameterization for that point. For the alternative 1 parameterization, for any pairs ( $c',c$ ) and ( $c",c'$ ), with $c < c' < c"$ , we also obtain

  (5) $$ \begin{align}\frac{\sigma_{Y_{t}^{*}}^{alt1}}{\sigma_{Y_{1}^{*}}^{alt1}}=\sigma_{Y_{t}^{*}}^{alt1}=\frac{\tau_{c'1}^{std}-\tau_{c1}^{std}}{\tau_{c't}^{std}-\tau_{ct}^{std}}=\frac{\tau_{c"1}^{std}-\tau_{c'1}^{std}}{\tau_{c"t}^{std}-\tau_{c't}^{std}}, \quad \quad \quad t=1, \ldots, T.\end{align} $$

  Therefore, this equation can be considered as a test of the assumption regarding the lack of invariance of at least one of the thresholds for a given time point.

• Alternative parameterization 2. Jöreskog (2001) proposed an alternative specification for longitudinal data. It is based on defining the origin and unit of measurement of $\mathbf {Y}_{i}^{*}$ in terms of thresholds (Bollen & Curran, 2006; Fisher & Bollen, 2020; Mehta et al., 2004). The common practice is to fix the distance between the first two thresholds as one unit on the new scale, that is $\tau _{2t}^{alt2} - \tau _{1t}^{alt2} = 1$ , for all $t=1, \dots , T$ . Equivalently, the first threshold could be fixed to zero and the second at one on each occasion, indicating that $Y_{it}^{*}$ on this new scale must be greater than zero for an individual to be in an ordinal category greater than one. This allows us to recover, on each occasion, both the mean and variance of $Y_{it}^{*}$ as well as the other ( $C-3$ ) thresholds in this new measurement scale. Given a number of ordinal categories equal to or greater than three, these conditions are necessary and sufficient to ensure the identification of the model linking $\mathbf {Y}_{i}$ to $\mathbf {Y}_{i}^{*}$ . This set of constraints also represents an alternative but equivalent parameterization of the standard one. Indeed, it can be easily shown that

  (6) $$ \begin{align} \begin{array}{ccc} \sigma_{Y_{t}^{*}}^{alt2}=\frac{1}{\tau_{2t}^{std}-\tau_{1t}^{std}},&\mu_{Y_{t}^{*}}^{alt2}=-\frac{\tau_{1t}^{std}}{\tau_{2t}^{std}-\tau_{1t}^{std}},&\tau_{ct}^{alt2}=\frac{\tau_{ct}^{std}-\tau_{1t}^{std}}{\tau_{2t}^{std}-\tau_{1t}^{std}},\quad c=3, \ldots, C.\\ \end{array} \end{align} $$

  A one-to-one relationship also stands between the two alternative parameterizations, as illustrated in Appendix A. In applying this parameterization in conjunction with a linear growth model for categorical data, Mehta et al. (2004) additionally imposed the assumption of threshold invariance over time. That is, $\tau _{ct}=\tau _{c}$ , for $c=3, \ldots , C$ , and all t. As discussed by Jöreskog (2001) and Mehta et al. (2004), this assumption implies that

  (7) $$ \begin{align} \tau_{c}^{alt2}= \mu_{Y_{t}^{*}}^{alt2}+\sigma_{Y_{t}^{*}}^{alt2}\tau_{ct}^{std}, \quad c =3, \ldots, C-1; t=1, \ldots T, \end{align} $$

  where $\tau _{ct}^{std}$ is the unconstrained threshold for the cth category at the t time point estimated under the standard parameterization. Eq. (7) sets constraints on $\mu _{Y_{t}^{*}}^{alt2}$ and $\sigma _{Y_{t}^{*}}^{alt2}$ because the right-hand side varies with t, whereas the left-hand side does not. If $C \geq 3$ , the common thresholds, $\mu _{Y_{t}^{*}}^{alt2}$ and $\sigma _{Y_{t}^{*}}^{alt2}$ can be estimated from the univariate marginal data of those variables whose thresholds are supposed to be equal. Eq. (7) can be rewritten for any two thresholds c and $c'$ for any $Y_{t}^{*}$ as $\tau _{ct}^{std}=(\tau _{c}^{alt2}-\mu _{Y_{t}^{*}}^{alt2})/\sigma _{Y_{t}^{*}}^{alt2}$ , and $\tau _{c't}^{std}=(\tau _{c'}^{alt2}-\mu _{Y_{t}^{*}}^{alt2})/\sigma _{Y_{t}^{*}}^{alt2},$ such that

  $$ \begin{align*}\sigma_{Y_{t}^{*}}^{alt2}=\frac{\tau_{c'}^{alt2}-\tau_{c}^{alt2}}{\tau_{c't}^{std}-\tau_{ct}^{std}}, \quad c < c'.\end{align*} $$

  It follows that the ratio of standard deviations on two different occasions is equal to the reciprocal of the ratio of the differences between any two standardized thresholds at those time points. For example, there are two independent equations with three thresholds at each time point. That is, for the pairs (2,1) and (3,2), we get

  (8) $$ \begin{align}\frac{\sigma_{Y_{t}^{*}}^{alt2}}{\sigma_{Y_{t'}^{*}}^{alt2}}=\frac{\tau_{2t'}^{std}-\tau_{1t'}^{std}}{\tau_{2t}^{std}-\tau_{1t}^{std}}=\frac{\tau_{3t'}^{std}-\tau_{2t'}^{std}}{\tau_{3t}^{std}-\tau_{2t}^{std}}.\end{align} $$

  If any one threshold for a time point under the alternative 2 parameterization is not invariant, eq. (8) fails to be true for that time point. Hence, this equation evaluates the lack of invariance of at least one of the thresholds at a given time point.

3.2 Binary case

Binary variables ($C=2)$ are special cases that require further comments. Indeed, with dichotomous variables, there are T available proportions, $\pi _{1t}=P(Y_{it}=1), t = 1, \ldots , T$ , to estimate the only threshold $\tau _{1t}$ , the mean and variance of the underlying variables at each occasion.

When the standard parameterization is adopted, the only available threshold is freely estimated on each occasion, whereas the means and variances of the underlying variables are always set to zero and one, respectively. Hence, no different constraints are placed with respect to the categorical case (see the first row of Table 3).

Table 3 Alternative sets of identification constraints for the auxiliary model in presence of binary data

On the other hand, additional restrictions have to be placed for the alternative parameterizations. Millsap and Yun-Tein (2004) suggest freely estimating the means of the underlying variables on each occasion except on the first one, setting the variances of the underlying variables to one on each occasion and estimating the only available threshold under the assumption of time-invariance. These constraints are detailed in the second row of Table 3. Differently, Muthén and Muthén (1998–2017) suggest freely estimating the underlying variable variances on all the occasions except on the first one, where $\sigma _{Y_{i1}^{*}}^{2}$ is fixed to one. They also suggest keeping the only available threshold invariant on each occasion. In this regard, the means of the underlying variables $Y_{it}^{*}$ have to be fixed to zero for the auxiliary model to be identified (see the third row in Table 3).

In the presence of binary data, it is impossible to implement the threshold constraints proposed by Jöreskog (2001), being only one threshold available per occasion. Even when holding thresholds invariant over time, this does not identify both means and variances on all occasions. Jöreskog (2001) proposes to fix the thresholds equal to zero and the variances $\sigma _{Y_{t}^{*}}^{2}$ to one on all occasions, only allowing the means $\mu _{Y_{t}^{*}}$ to vary over time. This is detailed in the last row in Table 3.

3.3 Estimation

SEMs for continuous endogenous variables analyze the mean vector $\boldsymbol \mu _{\mathbf {Y}}$ and covariance matrix $\boldsymbol \Sigma _{\mathbf {Y}\mathbf {Y}}$ of the observed indicators. However, when one or more of the endogenous observed variables are categorical, the analysis shifts to the mean vector $\boldsymbol \mu _{\mathbf {Y}^{*}}$ and covariance matrix $\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ corresponding to $\mathbf {Y}_{i}^{*}$ . If consistent estimators for $\boldsymbol \mu _{\mathbf {Y}^{*}}$ and $\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ are available, researchers can analyze them similarly to continuous indicators. Consequently, the estimation procedure comprises two distinct steps. The first step obtains consistent estimates of the means $\boldsymbol \mu _{\mathbf {Y}^{*}}$ and covariance matrix $\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ for $\mathbf {Y}_{i}^{*}$ . To perform significance testing, analysts also need the asymptotic covariance matrix of the elements in $\hat {\boldsymbol \mu }_{\mathbf {Y}^{*}}$ and $\hat {\boldsymbol \Sigma }_{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ . Once the means, variances, and covariances of $\mathbf {Y}_{i}^{*}$ are in hand, researchers can estimate the parameters of any longitudinal model they apply to $\mathbf {Y}_{i}^{*}$ .

Here, we focus on the first step needed to estimate $\boldsymbol \mu _{\mathbf {Y}^{*}}$ and the unconstrained covariance matrix $\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ . As we discussed in Section 3.1, distributional assumptions and identification constraints are necessary for estimation. Assuming a bivariate standard normal distribution for each pair of variables in $\mathbf {Y}_{i}^{*}$ facilitates the estimation of thresholds and polychoric correlations/covariances for noncontinuous variables. Univariate standard normality for each underlying variable $Y_{it}^{*}$ enables threshold estimation as percentiles of the standard normal distribution:

$$ \begin{align*}\tau_{ct}^{std}=\Phi^{-1}\left(\sum_{r=1}^{c}\frac{n_{rt}}{n}\right),\end{align*} $$

where $\Phi ^{-1}$ is the quantile function of the standard normal distribution, $n_{rt}$ is the number of cases in the rth category at time t, and n is the sample size. While univariate margins are valuable for estimating the thresholds in the standard parameterization, bivariate tables are essential for estimating the correlation between $Y_{it}^{*}$ and $Y_{it'}^{*}$ . Following Olsson (1979), the log-likelihood for the polychoric correlation is given by

$$ \begin{align*}\ln L=A+\sum_{l=1}^{C}\sum_{r=1}^{C}n_{lr}\ln(\pi_{lr}),\end{align*} $$

where C is the number of categories for both $Y_{it}^{*}$ and $Y_{it'}^{*}$ , $n_{lr}$ is the number of cases in the $lr$ th cell of the bivariate table, and A is an irrelevant constant that does not influence the values that maximize the likelihood. Hold the thresholds fixed at the estimates from the univariate margins,

$$ \begin{align*}\pi_{lr}=\Phi_{2}(\hat{\tau}_{lt},\hat{\tau}_{rt'})-\Phi_{2}(\hat{\tau}_{l-1,t}-\hat{\tau}_{rt'})-\Phi_{2}(\hat{\tau}_{lt},\hat{\tau}_{r-1,t'})+\Phi_{2}(\hat{\tau}_{l-1,t}-\hat{\tau}_{r-1,t'}),\end{align*} $$

where $\Phi _{2}(\cdot )$ is the bivariate normal distribution function with correlation $\rho _{Y_{it}^{*}Y_{it'}^{*}}$ . As a result of this conditional estimation procedure, $\hat {\rho }_{Y_{it}^{*}Y_{it'}^{*}}$ , for $t \neq t'$ , is a pseudo-maximum likelihood estimate.

When the means and variances of each $Y_{it}^{*}$ are of interest, it is possible to estimate these quantities rather than constraining them to 0 and 1, respectively. The alternative parameterizations of the auxiliary model allow this option. In the standard parameterization, assuming $Y_{it}^{*} \sim N(0,1)$ at each occasion allows for estimation for the mean (keeping the variance $\sigma _{Y_{t}^{*}}^{2}$ fixed to one) by shifting the distribution of $Y_{it}^{*}$ by a constant ${\mu _{Y_{it}^{*}}}^{alt}$ . This results in the transformed distribution as ${Y_{it}^{*}} \sim N({\mu _{Y_{it}^{*}}}^{alt},1)$ . If, in addition to the mean, the variance of $Y_{it}^{*}$ is estimated, the distribution of $Y_{it}^{*}$ is shifted and scaled, such that $Y_{it}^{*} \sim N({\mu _{Y_{it}^{*}}}^{alt}, {\sigma _{Y_{it}}^{*}}^{alt})$ . Hence, in the alternative parameterization of the auxiliary model, the corresponding thresholds, means, and variances of the underlying variables are based on their relationships with the thresholds of the standard parameterization. These relationships are given in Section 3.1 and detailed in eqs. (2) and (3) for the alternative parameterization proposed by Millsap and Yun-Tein (2004), and in eq. (6) for the parameterization proposed by Jöreskog (2001). If thresholds are assumed to be invariant over time, additional constraints are placed on the means and variances of the underlying variables under the alternative parameterizations, as discussed in Section 3.1.

The polychoric covariance matrix of the underlying variables under the alternative parameterization is then determined by scaling the polychoric correlation matrix $\hat {\mathbf {R}}_{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ estimated for the standard parameterization as follows:

$$ \begin{align*}\hat{\boldsymbol \Sigma}_{\mathbf{Y}^{*}\mathbf{Y}^{*}}=\mathbf{D}^{-1}\hat{\mathbf{R}}_{\mathbf{Y}^{*}\mathbf{Y}^{*}}\mathbf{D}^{-1},\end{align*} $$

where $\mathbf {D}$ is a $T \times T$ diagonal matrix with generic element $1/{\sigma _{Y_{t}^{*}}}^{alt}, t=1, \ldots , T$ .

3.4 Illustrative example: auxiliary model for NLSY97 data

Returning to our motivating NLSY97 example, we estimate the auxiliary measurement model connecting the binary and ordinal variables to their underlying continuous variables. The underlying variables $\mathbf {Y}_{i}^{*}$ represent propensities such as the inclination to use illegal drugs, experience depression, or perceive good/excellent health status. In the case of the illegal drug use variable, the threshold is the point that separates use from nonuse of illegal drugs. An individual whose propensity exceeds the threshold uses illegal drugs, whereas those who fall below it do not. Similarly, for the other two categorical variables, when the latent propensity or perception falls between thresholds $\tau _{ct}$ and $\tau _{c+1,t}$ on a given occasion, the observed ordinal response corresponds to category c.

As discussed in Section 3.1, various parameterizations of the auxiliary model are available. Both the standard auxiliary model and alternative parameterizations, where no threshold invariance constraints are imposed, result in just identified models (with zero degrees of freedom) that perfectly fit the data. When researchers fit these to the data, they obtain estimates of unknown thresholds, means of underlying variables, their variances, and polychoric correlations/covariances. For illustrative purposes, we report—in the Supplementary Material, due to space constraints—the estimated means and polychoric correlations/covariances for the auxiliary model that jointly considers all three observed variables based on the standard and alternative parameterizations. For the latter, we consider the models with and without threshold invariance.

The assumption of threshold invariance is tested by estimating the auxiliary model in which the same threshold at all time points is restricted to be equal. It is worth noting that both alternative (1 or 2) parameterizations of the auxiliary model are equivalent in their chi-square and degrees of freedom. Table 4 provides fit statistics for each categorical variable’s auxiliary model based on the alternative (1 or 2) parameterization with threshold invariance. For the binary variable, threshold invariance is mandatory for the auxiliary model to be just identified and is not testable. Under the assumption of threshold invariance, we have also estimated the multivariate auxiliary model for the three variables, which includes the estimation of thresholds, means, and variances for each variable, along with the polychoric covariances among all the underlying variables. Note that the threshold invariance hypothesis is not rejected for depressive symptoms at any conventional level and that for general health status is marginally significant. With over 5000 cases, statistical power should be high. We also checked the fit indexes and found that the CFI, TLI, and RMSEA all suggest excellent fit. The negative value of the BIC ( $= \chi ^{2} - \ln (n)df$ ) also support the invariant threshold models (Raftery, 1995).

Table 4 Fit statistics for the alternative(1 and 2) auxiliary models based on threshold invariance estimated for depressive symptoms, general health status, and for all the three variables jointly considered

These results suggest that the auxiliary model under alternative 1 or 2 with invariant thresholds is the most promising structure to use for these data when moving to the longitudinal model. However, these results are only with regard to the auxiliary model. In the following Sections, we consider the possible identification constraints of the different longitudinal models, and we turn to this topic next.

4 A general longitudinal model for categorical repeated measures

Given the latent response variates $\mathbf {Y}_{i}^{*}$ , we can apply the ALT model (Bollen & Curran, 2004). It encompasses the classical autoregressive and latent growth models as special cases. Following the general representation provided by the authors, the unconditional model is specified through two equations

(9) $$ \begin{align} \boldsymbol \eta_{i}&= \boldsymbol \mu_{\boldsymbol \eta} + \mathbf{B} \boldsymbol \eta_{i} + \boldsymbol \varsigma_{i}, \end{align} $$

(10) $$ \begin{align} \mathbf{Y}_{i}^{*}&=\mathbf{P}\boldsymbol \eta_{i}.~~~~~~~~~~~~~~~~~ \end{align} $$

In eq. (9), $\boldsymbol \eta _{i}$ is a vector that includes both the underlying variables $\mathbf {Y}_{i}^{*}$ and the random growth components $\boldsymbol \alpha _{i}' = (\alpha _{i0}, \alpha _{i1})$ . $\boldsymbol \mu _{\boldsymbol \eta }'=(\boldsymbol \nu _{\mathbf {Y}^{*}}', \boldsymbol \mu _{\boldsymbol \alpha }')$ represents a vector of intercepts and means, and the $(T+2) \times (T+2)$ matrix $\mathbf {B}$ specifies the coefficients for the relationships between the elements of $\boldsymbol \eta _{i}$ . It is divided into sub-matrices as follows

$$ \begin{align*}\mathbf{B}=\left[\begin{array}{ccc}\mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}}& \mathbf{B}_{\mathbf{Y}^{*}\boldsymbol \alpha}\\ \mathbf{O}&\mathbf{O}\end{array}\right],\end{align*} $$

where $\mathbf {B}_{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ contains the autoregressive effects among the underlying variables $\mathbf {Y}^{*}$ , that is

$$ \begin{align*}\mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}}=\left[\begin{array}{ccccc}0&0&\ldots&0&0\\\phi_{21}&0&\ldots&0&0\\0&\phi_{32}&\ldots&0&0\\\vdots&\vdots&\ldots&\vdots&\vdots\\0&0&\ldots&\phi_{T(T-1)}&0\end{array}\right],\end{align*} $$

and $\mathbf {B}_{\mathbf {Y}^{*}\boldsymbol \alpha }$ relates the underlying variables $\mathbf {Y}_{i}^{*}$ to the random growth components.

The autoregressive relations in the ALT model give rise to an initial condition problem since the variable at the start of the observation period, $Y_{i1}^{*}$ , should be affected by the random intercept and slope as well as unavailable pre-sample latent responses, say $Y_{i0}^{*}$ . To handle this problem, we can assume $Y_{i1}^{*}$ to be predetermined, correlated with the growth components, or treated as endogenous. We consider $Y_{i1}^{*}$ predetermined, knowing that, using rules from Lee and Hershberger (1990) and Hershberger (2006), the unconditional autoregressive latent trajectory model with $Y_{i1}^{*}$ predetermined or $Y_{i1}^{*}$ endogenous are (globally and covariance) equivalent (Bianconcini & Bollen, 2018). Hence, considering a linear growth curve, $\mathbf {B}_{\mathbf {Y}^{*}\boldsymbol \alpha }'=\left [\begin {array}{cccc}0&1&\ldots &1\\0&1&\ldots &T-1\end {array}\right ]$ .

The disturbance vector for $\boldsymbol \eta _{i}$ is $\boldsymbol \varsigma _{i}$ , that is $\boldsymbol \varsigma _{i}'=[\boldsymbol \varepsilon _{i}' \quad \boldsymbol \varsigma _{\boldsymbol \alpha }']$ . It is assumed to have a zero mean vector, and its covariances depend on the model. Under the assumption that $Y_{i1}^{*}$ is predetermined, $\boldsymbol \Sigma _{\boldsymbol \varsigma \boldsymbol \varsigma }$ is a block-diagonal matrix of the form $\left [\begin {array}{cc}\boldsymbol \Theta _{\boldsymbol \varepsilon }&\boldsymbol \Sigma _{\mathbf {Y^{*}}\boldsymbol \alpha }\\ \boldsymbol \Sigma _{\mathbf {Y^{*}}\boldsymbol \alpha }'&\boldsymbol \Sigma _{\boldsymbol \alpha \boldsymbol \alpha }\end {array}\right ]$ , where $\boldsymbol \Theta _{\boldsymbol \varepsilon }$ is a diagonal matrix with error variances that can differ over time, such that $V(\varepsilon _{it})=\sigma _{\varepsilon _{t}}^{2}$ , for $t=1, \ldots , T$ , $\boldsymbol \Sigma _{\mathbf {Y^{*}}\boldsymbol \alpha }$ contains the covariances between the predetermined $Y_{i1}^{*}$ and the random components $\boldsymbol \alpha _{i}$ , such that $\boldsymbol \Sigma _{\mathbf {Y^{*}}\boldsymbol \alpha }'=\left [\begin {array}{cccc}\sigma _{Y_{1}^{*}\alpha _{0}}&0&\ldots &0\\\sigma _{Y_{1}^{*}\alpha _{1}}&0&\ldots &0\end {array}\right ]$ , and $\boldsymbol \Sigma _{\boldsymbol \alpha \boldsymbol \alpha }$ is the full and symmetric covariance matrix of the growth components.

Equation (10) links the underlying response variates $\mathbf {Y}_{i}^{*}$ to the latent variables in $\boldsymbol \eta _{i}$ through the matrix $\mathbf {P}=[\mathbf {I}_{T} \quad \mathbf {O}]$ , where $\mathbf {I}_{T}$ is an identity matrix with dimensions that depend on the number of repeated measures and $\mathbf {O}$ is a zero matrix of dimensions $T \times 2$ . These assumptions lead to the following implied moments for the underlying variables in the autoregressive latent trajectory model of $\mathbf {Y}_{i}^{*}$

(11) $$ \begin{align} &\boldsymbol \mu_{\mathbf{Y}^{*}}=(\mathbf{I}_{T}- \mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}})^{-1}\nu_{\mathbf{Y}^{*}}+(\mathbf{I}_{T}- \mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}})^{-1}\mathbf{B}_{\mathbf{Y}^{*}\boldsymbol \alpha}\boldsymbol \mu_{\boldsymbol \alpha}, \end{align} $$

(12) $$ \begin{align}\boldsymbol \Sigma_{\mathbf{Y}^{*}\mathbf{Y}^{*}}&=(\mathbf{I}_{T}- \mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}})^{-1}\boldsymbol \Theta_{\boldsymbol \varepsilon}(\mathbf{I}_{T}- \mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}})^{-1'}+(\mathbf{I}_{T}- \mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}})^{-1}\mathbf{B}_{\mathbf{Y}^{*}\boldsymbol \alpha}\boldsymbol \Sigma_{\mathbf{Y^{*}}\boldsymbol \alpha}'(\mathbf{I}_{T}- \mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}})^{-1'}\\\nonumber& \quad +(\mathbf{I}_{T}- \mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}})^{-1}\boldsymbol \Sigma_{\mathbf{Y^{*}}\boldsymbol \alpha}\mathbf{B}_{\mathbf{Y}^{*}\boldsymbol \alpha}'(\mathbf{I}_{T}- \mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}})^{-1'}+(\mathbf{I}_{T}- \mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}})^{-1}\mathbf{B}_{\mathbf{Y}^{*}\boldsymbol \alpha}\boldsymbol \Sigma_{\boldsymbol \alpha \boldsymbol \alpha}\mathbf{B}_{\mathbf{Y}^{*}\boldsymbol \alpha}'(\mathbf{I}_{T}- \mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}})^{-1'}. \end{align} $$

The first term on the right-hand side of eqs. (11) and (12) are the moments implied by the autoregressive of order one component of the model, whereas the other terms account for the interaction between the autoregressive and growth components. The implied moments help determine the identification of the model parameters, test the model fit, and prove if the choice of different parameterizations is arbitrary or not.

Multivariate autoregressive latent trajectory model. If two or more binary or ordinal repeated measures are observed, the ALT model has to be generalized to deal with multiple series. To illustrate and clarify this model specification, consider that we have two series of repeated ordinal variables and that their underlying variables have autoregressive and cross-lagged relations with each other. Say that $\mathbf {Y}_{a,i}^{*}$ contains the series of longitudinal underlying variables for the first series and $\mathbf {Y}_{b,i}^*$ contains the longitudinal underlying variables for the second one. Based on our empirical application, we can consider $\mathbf {Y}_{a,i}^{*}$ to be the underlying variables associated with the observed depressive symptoms, while $\mathbf {Y}_{b,i}^{*}$ relates to the observed general health status variables. Let $\boldsymbol \alpha _{i0}$ and $\boldsymbol \alpha _{i1}$ be vectors of the growth components for the two series. In this situation, we write $\boldsymbol \eta _{i}^{T} = \left [\mathbf {Y}_{a,i}^* \hspace {0.1cm} \mathbf {Y}_{b,i}^*\hspace {0.1cm} \boldsymbol \alpha _{i0}\hspace {0.1cm} \boldsymbol \alpha _{i1} \right ]$ , and the vector $\boldsymbol \mu ^{T} = \left [\boldsymbol \mu _{\mathbf {Y}_a^*} \hspace {0.1cm } \boldsymbol \mu _{\mathbf {Y}_b^*} \hspace {0.1cm} \boldsymbol \mu _{\boldsymbol \alpha _{0}} \hspace {0.1cm} \boldsymbol \mu _{\boldsymbol \alpha _{1}} \right ]$ , where the $\boldsymbol \mu $ ’s are intercepts of endogenous variables or means of latent variables, and the subscript signifies which variable. Next is the $\mathbf {B}$ matrix,

$$ \begin{align*}\mathbf{B}=\left[\begin{array}{cccc} \mathbf{B}_{\mathbf{Y}_{a}^{*}\mathbf{Y}_{a}^{*}}&\mathbf{B}_{\mathbf{Y}_{a}^{*}\mathbf{Y}_{b}^{*}}&\mathbf{B}_{\mathbf{Y}_{a}^{*}\boldsymbol \alpha_{0}}&\mathbf{B}_{\mathbf{Y}_{a}^{*}\boldsymbol \alpha_{1}}\\ \mathbf{B}_{\mathbf{Y}_{b}^{*}\mathbf{Y}_{a}^{*}}&\mathbf{B}_{\mathbf{Y}_{b}^{*}\mathbf{Y}_{b}^{*}}&\mathbf{B}_{\mathbf{Y}_{b}^{*}\boldsymbol \alpha_{0}}&\mathbf{B}_{\mathbf{Y}_{b}^{*}\boldsymbol \alpha_{1}}\\ \mathbf{B}_{\boldsymbol \alpha_{0}\mathbf{Y}_{a}^{*}}&\mathbf{B}_{\boldsymbol \alpha_{0}\mathbf{Y}_{b}^{*}}&\mathbf{B}_{\boldsymbol \alpha_{0}\boldsymbol \alpha_{0}}&\mathbf{B}_{\boldsymbol \alpha_{0}\boldsymbol \alpha_{1}}\\ \mathbf{B}_{\boldsymbol \alpha_{1}\mathbf{Y}_{a}^{*}}&\mathbf{B}_{\boldsymbol \alpha_{1}\mathbf{Y}_{b}^{*}}&\mathbf{B}_{\boldsymbol \alpha_{1}\boldsymbol \alpha_{0}}&\mathbf{B}_{\boldsymbol \alpha_{1}\boldsymbol \alpha_{1}}\\ \end{array}\right],\end{align*} $$

which is partitioned. As before, the first subscript of the coefficients shows the variables receiving effects, while the second subscript is the variables emitting effects. For instance, $\mathbf {B}_{\mathbf {Y}_a^*\mathbf {Y}_b^*}$ gives the coefficients of the direct effects of the second series of repeated variables on the first series, this would include any cross-lagged effects from Y _{b, i} ^* to Y _{a, i} ^*. Similarly, $\mathbf {B}_{\mathbf {Y}_b^*\mathbf {Y}_a^*}$ gives the direct effects of the first series on the second and $\mathbf {B}_{\mathbf {Y}_b^*\boldsymbol \alpha _{i1}}$ contains the coefficients of the direct effects of the random slopes of both series on the repeated measures of the second series. Knowing that none of the repeated measures has direct effects on the random intercepts or random slopes, and no direct effects are estimated for the random components, all coefficients that correspond to such effects are set to zero; that is,

$$ \begin{align*}\mathbf{B}=\left[\begin{array}{cccc} \mathbf{B}_{\mathbf{Y}_{a}^{*}\mathbf{Y}_{a}^{*}}&\mathbf{B}_{\mathbf{Y}_{a}^{*}\mathbf{Y}_{b}^{*}}&\mathbf{B}_{\mathbf{Y}_{a}^{*}\boldsymbol \alpha_{0}}&\mathbf{B}_{\mathbf{Y}_{a}^{*}\boldsymbol \alpha_{1}}\\ \mathbf{B}_{\mathbf{Y}_{b}^{*}\mathbf{Y}_{a}^{*}}&\mathbf{B}_{\mathbf{Y}_{b}^{*}\mathbf{Y}_{b}^{*}}&\mathbf{B}_{\mathbf{Y}_{b}^{*}\boldsymbol \alpha_{0}}&\mathbf{B}_{\mathbf{Y}_{b}^{*}\boldsymbol \alpha_{1}}\\ \mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\ \end{array}\right].\end{align*} $$

The vector $\mathbf {\zeta }_{i}^{T}$ is now given by $\left [ \boldsymbol \varepsilon _{\mathbf {Y}_{a,i}} \hspace {0.1cm}\boldsymbol \varepsilon _{\mathbf {Y}_{b,i}} \hspace {0.1cm}\boldsymbol \varsigma _{\boldsymbol \alpha _{i0}} \hspace {0.1cm}\boldsymbol \varsigma _{\boldsymbol \alpha _{i1}}\right ]$ . The first two elements, $\boldsymbol \varepsilon _{\mathbf {Y}_{a,i}}$ and $\boldsymbol \varepsilon _{\mathbf {Y}_{b,i}}$ , are disturbances or errors of the two repeated measures series. These errors are correlated when referring to the same occasion (concomitant effects). The last two vectors, $\boldsymbol \varsigma _{\boldsymbol \alpha _{i0}}$ and $\boldsymbol \varsigma _{\boldsymbol \alpha _{i1}}$ , are the disturbances or errors from the equations for $\boldsymbol \alpha _{i0}$ and $\boldsymbol \alpha _{i1}$ , assumed to be all correlated. In brief, the preceding equations enable the analysis of two repeated measures with binary or ordinal variables. Generalizing the model to analyze three or more series of repeated measures is straightforward.

4.1 Identification issues

The ALT model for continuous measures requires constraints to ensure its identification, particularly with fewer than five waves of data (Bollen & Curran, 2004). Extending this model—and its special cases—to ordered categorical measures presents additional challenges. We can approach the identification problem by dividing the constraints needed into two parts. The first part contains constraints to identify $(\boldsymbol \mu _{\mathbf {Y}^{*}}, \boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}})$ . The second part contains conditions to identify the autoregressive latent trajectory model parameters. For the latter, we will also discuss the identification issues for the linear latent growth and autoregressive of order one models.

Once the auxiliary model is accurately identified and estimated, we treat the means and covariances between the underlying variables $\mathbf {Y}_{i}^{*}$ as known and use them to identify the longitudinal model for $\mathbf {Y}_{i}^{*}$ .

This section outlines identification conditions for the linear latent growth, the first-order autoregressive, and ALT models. We explore alternative parameterizations for these models based on the previously discussed auxiliary model specifications. Our assessment focuses on whether distinct specifications lead to equivalent models. In SEMs, equivalent models generate identical model-implied moment matrices and equally fit the data, with equal test statistics, fit indexes, and degrees of freedom (Lee & Hershberger, 1990; Stelzl, 1986). Equivalent models impose the same constraints on the population covariance matrix, known as $\boldsymbol \Sigma $ constraints (Steiger, 2002). Models with the same $\boldsymbol \Sigma $ constraints are $\boldsymbol \Sigma $ -equivalent since they cannot be empirically distinguished.

A formal definition of model equivalence has been provided by Raykov and Penev (1999), who established a necessary and sufficient condition based on the existence of a parameter transformation that preserves the implied covariance matrix and covers the entire parameter space of the other model. This implies that equivalent models remain invariant under this transformation. Validating this condition involves deriving implied covariance matrices, equating corresponding elements, and solving for parameters to confirm the $\boldsymbol \Sigma $ -equivalence. If there is no mapping, the two models are not equivalent. An extension of the rule by Raykov and Penev (1999) for identifying equivalent models with a mean structure has been discussed by Levy and Hancock (2007) and Losardo (2009) with a specific focus on the equivalence of latent curve model specifications. We follow the approach of Losardo (2009), accommodating the presence of categorical repeated measures.

4.1.1 Linear growth model

We begin by considering the linear latent growth model for $\mathbf {Y}_{i}^{*}$ . Three different specifications correspond to the standard and alternative parameterizations of the auxiliary model. Despite the time-specific scale of thresholds in the standard parameterization, we consider it being the default option in Mplus 8.6 and lavaan 0.6-16.

To derive identification constraints for each alternative specification, as detailed in Table 5, we need the mean and covariance matrix implied by the linear growth model. These are given by

(13) $$ \begin{align} \boldsymbol \mu_{\mathbf{Y}^{*}}&=\mathbf{B}_{\mathbf{Y}^{*}\boldsymbol \alpha}\boldsymbol \mu_{\boldsymbol \alpha},~~~~~~~~~~~~~~~ \end{align} $$

(14) $$ \begin{align} \boldsymbol \Sigma_{\mathbf{Y}^{*}\mathbf{Y}^{*}}&=\mathbf{B}_{\mathbf{Y}^{*}\boldsymbol \alpha}\boldsymbol \Sigma_{\boldsymbol \alpha\boldsymbol \alpha}\mathbf{B}_{\mathbf{Y}^{*}\boldsymbol \alpha}' + \boldsymbol \Theta_{\boldsymbol \varepsilon}. \end{align} $$

Table 5 Identification constraints and degrees of freedom for the linear growth model for categorical repeated measures

For a deeper understanding of the relationship between these specifications, Table 5 also reports the corresponding degrees of freedom.

The WLSMV estimator only uses information coming from the first- and second-order sample statistics corresponding to the $T(C-1)$ univariate proportions and the $T(T-1)/2$ polychoric correlations. The degrees of freedom of a given specification are the difference in the number of available information, $T(C-1)+T(T-1)/2$ , and the number of parameters in the considered model. It is important to note that not all three specifications have the same degrees of freedom. Hence, the standard parameterization would not provide an equivalent specification to the alternative ones. The assumptions made for the identification of the auxiliary model impact the specification and interpretation of the linear latent growth model defined for $\mathbf {Y}_{i}^{*}$ . The two alternative parameterizations provide equivalent specifications to each other since a one-to-one relationship between their parameters exists, as detailed in Table 6.

Table 6 Parameter transformations for the linear growth models based on the parameterization proposed by Muthén and Muthén (1998–2017) (Model A) and the one illustrated by Mehta et al. (2004) (Model B)

When adopting the standard parametrization (second column of Table 5), Stage 1 constraints are sufficient to identify the auxiliary model. Setting ( $\boldsymbol \mu _{\mathbf {Y}^{*}},\texttt {diag}(\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ )) to ( $\mathbf {0},\mathbf {I}$ ) identifies the thresholds $\tau _{ct}, c=1, \ldots , C-1$ , on all occasions. Stage 2 constraints involve fixing $\boldsymbol \mu _{\boldsymbol \alpha }$ to $\mathbf {0}$ to satisfy the implied mean condition in eq. (13). Appendix B shows that these conditions are sufficient to identify the linear latent growth model., followed by the fact that since the implied mean is a function of the constant matrix $\mathbf {B}_{\mathbf {Y}^{*}\boldsymbol \alpha }$ and of the vector of growth means $\boldsymbol \mu _{\boldsymbol \alpha }$ , the latter has to be fixed equal to zero for the implied mean condition in eq. (13) to be satisfied.

A key distinction between SEMs for continuous outcomes and models for categorical data lies in the identification of the error $\boldsymbol \varepsilon _{i}$ variances, $\boldsymbol \Theta _{\boldsymbol \varepsilon }$ . The variances of the errors cannot be independently identified from the variances of the underlying variables, diag( $\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ ), as $\mathbf {Y}_{i}^{*}$ is a latent vector without an inherent scale. If diag( $\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}})=\mathbf {I}$ , the variance–covariance matrix of the growth components $\boldsymbol \Sigma _{\boldsymbol \alpha \boldsymbol \alpha }$ is correctly identified, while the error variances $\boldsymbol \Theta _{\boldsymbol \varepsilon }$ are determined as a remainder based on eq. (14).

An alternative set of identification constraints replaces diag( $\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}})=\mathbf {I}$ with $\boldsymbol \Theta _{\boldsymbol \varepsilon }=\mathbf {I}$ . It is denoted as (standard) theta parameterization in Mplus 8.6 and lavaan 0.6-16 and termed conditional parameterization by Kamata and Bauer (2008). This approach assumes standardized conditional distributions of $\mathbf {Y}_{i}^{*}$ given $\boldsymbol \eta _{i}$ , with marginal variances of $\mathbf {Y}_{i}^{*}$ obtained as the remainder based on eq. (14). In longitudinal data, this assumption could be more suitable in that the variances of $\mathbf {Y}_{i}^{*}$ are permitted to vary over time while the error variances are not, and this latter assumption might be more plausible. Of course, if there are reasons to think that these error variances differ over time, this assumption is also not desirable.

To address the conflicting conclusions drawn by Grimm and Liu (2016), Lee et al. (2018) and Newsom and Smith (2020) regarding the relationship between these standard parameterizations of the linear growth model, we provide a theoretical illustration. We limit our analyses to models with four time points to avoid complex and unproductive mathematical details. We compare the two linear growth model parameterizations (based on diag $(\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}})=\mathbf {I}$ or $\boldsymbol \Theta _{\boldsymbol \varepsilon }=\mathbf {I}$ ) in terms of $\boldsymbol \Sigma $ and $\boldsymbol \mu $ constraints. Following Steiger (2002), the covariance matrix and the mean vector implied by each standard (delta and theta) model specification are first derived and reported in the Supplementary Material. Both model specifications have no $\boldsymbol \mu $ constraints and $T(T+1)/2-3$ (seven) $\boldsymbol \Sigma $ constraints corresponding to the degrees of freedom on the covariance structure of the model, as illustrated in Table 7.

Table 7 $\boldsymbol \Sigma $ constraints for the linear growth model based on the (standard) delta parameterization (diag( $\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}})=\mathbf {I}$ ) and the (standard) theta parameterization ( $\boldsymbol \Theta _{\boldsymbol \varepsilon }=\mathbf {I}$ )

It is evident that these two standard specifications of the linear growth model are not empirically equivalent, as the constraints on the underlying variable variances are different, and the $\boldsymbol \Sigma $ constraints derived when all the underlying variable variances are fixed to one are defined in terms of polychoric correlations $\rho _{t,t'}$ .

When adopting the alternative parameterization proposed by Muthén and Asparouhov (2002) and by Millsap and Yun-Tein (2004) for the auxiliary model, which fixes the mean of the underlying variable to zero only on the first occasion, the implied mean condition in eq. (13) necessitates constraining only the mean of the intercept component, $\mu _{\alpha _{0}}$ , to zero (refer to the third column of Table 5). The variance of the underlying variable on the first occasion is set to one, while the other ( $T-1$ ) variances are freely estimated. Due to the interdependence of $\sigma _{Y_{t}^{*}}^{2}$ and $\sigma _{\varepsilon _{t}}^{2}$ , an alternative specification is obtained by placing these variance restrictions on the errors rather than on the underlying variable. In Appendix B, we prove that these conditions are sufficient for the model identification.

In contrast to the standard parameterization, it can be demonstrated that placing restrictions on diag( $\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ ) or $\boldsymbol \Theta _{\boldsymbol \varepsilon }$ defines two specifications of the linear growth model that are equivalent. Table 8 presents the parameter transformations for the model where all underlying variable variances but one are freely estimated (Model C) and the model where all error variances are free parameters except on the first occasion (Model D).

Table 8 Parameter transformations for the linear growth models where all underlying variable variances are free except on the first occasion (Model C) and when all error variances are free except at the first time point (Model D)

While these transformations are not identity functions, it is evident that parameters in Model C are functions of the corresponding parameters in Model D divided by the underlying variable variance on the first occasion under that specification, that is, ${\sigma _{{Y}^{*}_{1}}^{2D}}=1+ \sigma _{\alpha _{0}}^{2D}$ . Conversely, parameters in the theta specification are functions of the corresponding parameters in the Model C parameterization divided by $1- \sigma _{\alpha _{0}}^{2C}$ , where $\sigma _{\alpha _{0}}^{2C}$ represents the intercept variance in the model where all underlying variable variances are freely estimated except on the first occasion.

Finally, following the alternative 2 parameterization suggested by Jöreskog (2001), where all means and variances of the underlying variables are freely estimated, no restrictions are needed for the identification of the growth parameters, as detailed in the fourth column on Table 5 and proven in Appendix B. For this specification, focusing on either the $\mathbf {Y}_{i}^{*}$ or error variances yields the same $\boldsymbol \Sigma $ and $\boldsymbol \mu $ constraints, as illustrated in Table 9. Hence, these two alternative linear growth model specifications are empirically equivalent, with parameters of the two specifications related via identity transformations.

Table 9 $\boldsymbol {\Sigma }$ and $\boldsymbol \mu $ constraints for the linear growth model based on the alternative 2 delta parameterization (diag $(\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}})$ freely estimated) and the alternative 2 theta parameterization $(\boldsymbol \Theta _{\boldsymbol \varepsilon }$ freely estimated)

Binary data. A special case that requires further comments arises when the measured variables are dichotomous $(C=2)$ . Four different specifications of the linear latent growth model can be derived based on the various parameterizations of the auxiliary model detailed in Table 3. Identification constraints for these specifications are provided in Table 10, along with the corresponding degrees of freedom. The degrees of freedom are computed as the difference between the number of available information, equal to $T+T(T-1)/2$ , and the number of parameters in each model.

Table 10 Identification constraints and degrees of freedom for the linear growth model for binary data

All the constraints in Table 10 are sufficient for identifying the corresponding linear latent growth model. The proof follows the same line as detailed in Appendix B in the presence of categorical data. Attention should be paid to the alternative specification suggested by Muthén and Muthén (1998–2017), whose constraints are reported in the fourth column in Table 10. Due to the dependence of the underlying variable moments $(\boldsymbol \mu _{\mathbf {Y}^{*}}, \boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ ) on the parameters of the growth model $(\boldsymbol \mu _{\boldsymbol \alpha }, \boldsymbol \Sigma _{\boldsymbol \alpha \boldsymbol \alpha })$ , based on eqs. (13) and (14), we can employ fewer Stage 1 constraints than those outlined in Table 10 by taking advantage of the model structure. In this case, the identification problem is treated simultaneously for the thresholds, $(\boldsymbol \mu _{\mathbf {Y}^{*}}, \boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ ), and for the latent growth parameters $(\boldsymbol \mu _{\boldsymbol \alpha }, \boldsymbol \Sigma _{\boldsymbol \alpha \boldsymbol \alpha }$ ). A formal proof is provided in Appendix B.

It is important to note that only the alternative parameterization proposed by Millsap and Yun-Tein (2004) and the one suggested by Jöreskog (2001) share the same number of degrees of freedom. Furthermore, they are equivalent specifications of the linear growth model, since a surjective transformation exists that expresses the parameters of one specification as a function of those of the other and vice versa. All the parameters ( $\sigma _{\alpha _{0}}^{2},\sigma _{\alpha _{0}\alpha _{1}},\sigma _{\alpha _{1}}^{2}, \mu _{\alpha _{1}}$ ) are related by identity functions except for the threshold $\tau _{1}$ in the alternative parameterization by Millsap and Yun-Tein (2004) that is equal to minus the intercept mean $\mu _{\alpha _{0}}$ in the Jöreskog (2001) parameterization and vice versa.

For each parameterization, an alternative set of sufficient conditions can be derived by replacing the constraints on diag( $\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ ) with corresponding constraints on the error variances $\boldsymbol \Theta _{\boldsymbol \varepsilon }$ . The results derived for categorical data directly apply when the models are fitted to binary observations: anytime the underlying variable or error variances are all fixed to one, the two alternative parameterizations of the linear growth model are not empirically equivalent. This occurs when the standard parameterization, the alternative one proposed by Millsap and Yun-Tein (2004), and that suggested by Jöreskog (2001) are considered. On the other hand, for the parameterization suggested by Muthén and Muthén (1998–2017), the two parameterizations based on freely estimating all but the first underlying variable or error variances are empirically equivalent. The parameters of the two specifications can be expressed as one function of the other, as detailed in Table 8.

4.1.2 First-order autoregressive model

We now consider alternative parameterizations for the first-order autoregressive model. In situations where the dependent variable is influenced by or influences other dependent variables, assumptions need to be placed on the error variances due to the improper parameter constraints that come into play when focusing on the variances of underlying variables (Muthén and Muthén, 1998–2017, pp. 485–486). Specifically, when autoregressive components are present, Mplus exclusively supports theta parameterizations.

The mean and covariance matrix implied by the first-order autoregressive model can be expressed as follows

(15) $$ \begin{align} \boldsymbol \mu_{\mathbf{Y}^{*}}&= (\mathbf{I}_{T}- \mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}})^{-1}\nu_{\mathbf{Y}^{*}}, ~~~~~~~~~~~~~~~~~~~~~\end{align} $$

(16) $$ \begin{align} \boldsymbol \Sigma_{\mathbf{Y}^{*}\mathbf{Y}^{*}}&=(\mathbf{I}_{T}- \mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}})^{-1}\boldsymbol \Theta_{\boldsymbol \varepsilon}(\mathbf{I}_{T}- \mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}})^{-1'}.\end{align} $$

Here, $ \mathbf {B}_{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ contains all autoregressive effects among underlying variables, such that $(\mathbf {I}-\mathbf {B}_{\mathbf {Y}^{*}\mathbf {Y}^{*}})^{-1}$ has the following specific structure

$$ \begin{align*}(\mathbf{I}-\mathbf{B}_{\mathbf{Y}^{*}\mathbf{Y}^{*}})^{-1}=\left[\begin{array}{ccccc}1&0&\ldots&0&0\\\phi_{21}&1&\ldots&0&0\\\phi_{32}\phi_{21}&\phi_{32}&\ldots&0&0\\\vdots&\vdots&\cdots&\vdots&\vdots\\\phi_{T(T-1)} \ldots \phi_{21}&\phi_{T(T-1)} \ldots \phi_{32}&\ldots&\phi_{T(T-1)}&1\end{array}\right].\end{align*} $$

Identification conditions for different autoregressive model specifications, along with corresponding degrees of freedom, are presented in Table 11. Appendix B proves the sufficiency of these constraints for model identification. Degrees of freedom are computed based on the difference between information from the first- and second-order statistics, equal to $T(C-1)+T(T-1)/2$ , and the number of parameters in each specification.

Table 11 Identification constraints and degrees of freedom for the autoregressive model for categorical repeated measures

Coherently with the analysis performed for the linear latent growth model, we still consider alternative parameterizations based on the assumption of threshold invariance. Understanding how constraints on the auxiliary model imply different or equivalent autoregressive model specifications is crucial. Alternative parameterizations endow the model with $(C-4)(T-1)+T(T-1)/2$ degrees of freedom, while the standard (theta) parameterization, which involves fixing the error means and variances on all occasions, has $(T-2)(T-1)/2$ degrees of freedom.

When the assumption of threshold invariance is relaxed, a notable distinction is observed between the linear latent growth and autoregressive models. The standard and alternative autoregressive specifications are equivalent, characterized by $(T-2)(T-1)/2$ degrees of freedom. This is unique to the autoregressive model, as the standard specification of the linear growth model for categorical data still differs from the alternative ones.

In the Supplementary Material, we illustrate the model-implied mean and covariance structures for each autoregressive model, focusing on a simplified scenario with four observed time points. Based on these implied moments, we show that the alternative 1 (Model E) and 2 (Model F) parameterizations, both based on threshold invariance, share the same $\boldsymbol \Sigma $ constraints, represented by $\sigma _{3,2}\sigma _{4,1}=\sigma _{4,2}\sigma _{3,1}$ , $\sigma _{3,2}\sigma _{2,1}=\sigma _{3,1}\sigma _{Y_{2}^{*}}^{2}$ , and $\sigma _{4,3}\sigma _{3,1}=\sigma _{4,1}\sigma _{Y_{3}^{*}}^{2}$ . Table 12 outlines the relationships between parameters of these equivalent autoregressive model specifications.

Table 12 Parameter transformations for the autoregressive models based on the parameterization proposed by Muthén and Muthén (1998–2017) (Model E) and the one illustrated by Mehta et al. (2004) (Model F)

It is evident that Model F parameters are functions of the estimated error variance on the first occasion in Model E, while Model E parameters depend on the distance between the two first thresholds in Model F.

Binary data. The distinct characteristics of various parameterizations of the autoregressive model of order one become more apparent when applied to binary data. In this context, all considered specifications yield the same degrees of freedom, as presented in Table 13. It is straightforward to demonstrate that the parameters of each model can be expressed as functions of the others. To simplify the presentation, Table 14 just illustrates how parameters in the standard parameterization (Model G) relate to those in each of the alternative parameterizations.

Table 13 Identification constraints and degrees of freedom for the autoregressive model for binary data

Table 14 Parameter transformations for the autoregressive model of order one for binary data based on the standard parameterization of the auxiliary model (Model G), following the alternative parameterization proposed by Millsap and Yun-Tein (2004) (Model H), that proposed by Muthén and Muthén (1998–2017) (Model I), and the one proposed by Jöreskog (2001) (Model L)

Autoregressive latent trajectory model

The ALT model integrates components from the previously discussed models, the linear latent growth and the first-order autoregressive model. Including autoregressive effects requires constraints to be specifically placed on the error variances rather than on the underlying variable variances, as discussed earlier.

In eqs. (11) and (12), the mean vector and covariance matrix implied by the ALT model are directly influenced by the autoregressive component. This influence is evident through terms such as $(\mathbf {I}_{T}- \mathbf {B}_{\mathbf {Y}^{*}\mathbf {Y}^{*}})^{-1}\nu _{\mathbf {Y}^{*}}$ and $(\mathbf {I}_{T}- \mathbf {B}_{\mathbf {Y}^{*}\mathbf {Y}^{*}})^{-1}\boldsymbol \Theta _{\boldsymbol \varepsilon }(\mathbf {I}_{T}- \mathbf {B}_{\mathbf {Y}^{*}\mathbf {Y}^{*}})^{-1'},$ respectively.

Additionally, both moments depend on the interaction between the autoregressive and growth components. The Supplementary Material provides a detailed presentation of these moments for various ALT model specifications under the simplified scenario of a stationary autoregressive process with four observed occasions.

Table 15 presents identification conditions for different ALT parameterizations. These conditions are sufficient for identification. They can be easily proven in a manner similar to that for the linear latent growth and autoregressive model in Appendix B. The degrees of freedom for each specification are also reported. Consistent with the previous analysis, we consider alternative parameterizations based on the assumption of threshold invariance.

Table 15 Identification constraints and degrees of freedom for the autoregressive latent trajectory model for categorical repeated measures

As observed for the linear growth model, both alternative parameterizations share the same degrees of freedom, equal to $C(T-1)+\frac {T(T-7)}{2}-4$ , while the standard parameterization has $T(T-3)/2-4$ degrees of freedom. The model requires five occasions for identification without imposing additional parameter constraints beyond those outlined in Table 15.

A notable difference from both the linear latent growth and autoregressive model is that the two alternative parameterizations of the ALT model are not empirically equivalent. They exhibit distinct $\boldsymbol \mu $ and $\boldsymbol \Sigma $ constraints. To illustrate this point clearly, we consider the simplified case of a stationary first-order autoregressive model with five observed occasions. The $\boldsymbol \mu $ constraints for the different autoregressive latent trajectory specifications are illustrated in the following Table 16.

Table 16 $\boldsymbol \mu $ constraints for the different autoregressive latent trajectory model specifications

The three ALT specifications also exhibit different $\boldsymbol \Sigma $ constraints, which are detailed in the Supplementary Material due to space constraints. Attempting to establish a correspondence between the parameters of the two alternative parameterizations, as requested by Raykov and Penev (1999), proves unsuccessful when equating corresponding elements in the implied covariances matrices. Consequently, the assumptions made for the identification of the auxiliary model strongly influence the specification and interpretation of the ALT model for $\mathbf {Y}_{i}^{*}$ .

Binary data. Moving on to the analysis of binary data, the ALT model inherits the discrepancies observed in the alternative specifications for categorical data. In this context, none of the proposed parameterizations of the auxiliary model leads to equivalent specifications of the ALT model. Despite the fact that the alternative specifications proposed by Millsap and Yun-Tein (2004) and Jöreskog (2001) yield the same degrees of freedom, as presented in Table 17, the parameters of each model cannot be expressed as functions of the others, as required by Raykov and Penev (1999).

Table 17 Identification constraints and degrees of freedom for the autoregressive latent trajectory model for binary data

Table 17 provides sufficient identification conditions for each model, along with corresponding degrees of freedom. Five waves of data are necessary—under each parameterization—to identify the model without placing additional constraints on the model parameters. This is coherent with what was observed in the presence of continuous observations (Bollen & Curran, 2004). Based on these findings, researchers should carefully consider the implications of selecting specific parameterizations of the autoregressive latent trajectory model, as evidenced by the NLSY97 data analysis, where different choices significantly influenced the interpretation of the results.

4.2 Estimation

The correlation matrix (standard parameterization) or the mean vector and unconstrained covariance matrix (alternative parameterizations) estimated based on the auxiliary model specification—as described in Section 3.3—are utilized to derive point estimates for the parameters of the dynamic model specified for $\mathbf {Y}_{i}^{*}$ . Various estimators apply to estimate the structural parameters in a SEM. We specifically focus on the Diagonally Weighted Least Squares (DWLS) estimator, which is the default choice in Mplus 8.6 and lavaan 0.6-16. Once $\hat {R}_{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ or ( $\hat {\boldsymbol \mu }_{\mathbf {Y}^{*}},\hat {\boldsymbol \Sigma }_{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ ) are available, all parameters of SEM are estimated simultaneously.

To begin, we organize the estimated means $\hat {\boldsymbol \mu }_{\mathbf {Y}^{*}}$ and all the diagonal and below diagonal elements in $\hat {\boldsymbol \Sigma }_{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ (or $\hat {R}_{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ ) into a vector $\hat {\boldsymbol \rho }$ . Similarly, we place the implied moments $\boldsymbol \mu _{Y^{*}}$ and $\boldsymbol \Sigma _{\mathbf {Y}^{*}\mathbf {Y}^{*}}$ in a vector $\boldsymbol \rho (\boldsymbol \theta )$ , where $\boldsymbol \theta $ contains the SEM parameters. The DWLS estimator is then determined by minimizing

$$ \begin{align*}F_{DWLS}=[\hat{\boldsymbol \rho}-\boldsymbol \rho(\boldsymbol \theta)]'\texttt{diag}(\boldsymbol \Sigma_{\hat{\boldsymbol \rho}\hat{\boldsymbol \rho}})^{-1}[\hat{\boldsymbol \rho}-\boldsymbol \rho(\boldsymbol \theta)],\end{align*} $$

where $\boldsymbol \Sigma _{\hat {\boldsymbol \rho }\hat {\boldsymbol \rho }}$ is the asymptotic covariance matrix of $\hat {\boldsymbol \rho }$ whose derivation has been widely detailed in Muthén (1984), Jöreskog (1994), and Muthén and Satorra (1995).

The parameters $\boldsymbol \theta $ are chosen to minimize the weighted sum of squared deviations of $[\hat {\boldsymbol \rho }-\boldsymbol \rho (\boldsymbol \theta )]$ . The $F_{DWLS}$ is consistent, asymptotically unbiased, and normally distributed (Browne, 1984). However, it lacks asymptotic efficiency. Default standard errors are no longer accurate, and goodness of fit tests no longer follow a $\chi ^{2}$ distribution. To address this, robust standard errors are obtained by considering the following sandwich-type asymptotic covariance matrix of the parameter estimates $\hat {\boldsymbol \theta }$ (Muthén et al., 1997)

$$ \begin{align*}\nonumber \boldsymbol \Sigma_{\hat{\boldsymbol \theta}\hat{\boldsymbol \theta}}&=n^{-1}\left[\left(\frac{\partial \boldsymbol \rho(\boldsymbol \theta)}{\partial \boldsymbol \theta}\right)'\texttt{diag}(\boldsymbol \Sigma_{\hat{\boldsymbol \rho}\hat{\boldsymbol \rho}})^{-1}\left(\frac{\partial \boldsymbol \rho(\boldsymbol \theta)}{\partial \boldsymbol \theta}\right)\right]^{-1}\left(\frac{\partial \boldsymbol \rho(\boldsymbol \theta)}{\partial \boldsymbol \theta}\right)'\texttt{diag}(\boldsymbol \Sigma_{\hat{\boldsymbol \rho}\hat{\boldsymbol \rho}})^{-1}\boldsymbol \Sigma_{\hat{\boldsymbol \rho}\hat{\boldsymbol \rho}}^{-1}\\\nonumber &\quad\times\texttt{diag}(\boldsymbol \Sigma_{\hat{\boldsymbol \rho}\hat{\boldsymbol \rho}})^{-1}\left(\frac{\partial \boldsymbol \rho(\boldsymbol \theta)}{\partial \boldsymbol \theta}\right)\left[\left(\frac{\partial \boldsymbol \rho(\boldsymbol \theta)}{\partial \boldsymbol \theta}\right)'\texttt{diag}(\boldsymbol \Sigma_{\hat{\boldsymbol \rho}\hat{\boldsymbol \rho}})^{-1}\left(\frac{\partial \boldsymbol \rho(\boldsymbol \theta)}{\partial \boldsymbol \theta}\right)\right]^{-1}. \end{align*} $$

The square root of the main diagonal at the estimated parameters represents the robust standard errors of the parameter estimates. Mean and variance-adjusted chi-square statistics have been proposed to approximate the shape of the test statistics to the reference chi-square distribution with the associated degrees of freedom. The WLSMV estimator, which is the default estimator in Mplus 8.6 and lavaan 0.6-16 for models with endogenous categorical variables, relies on the Satterthwaite (1941) type correction. We refer the reader to Satorra and Bentler (1994) and Muthén et al. (1997) for its detailed description.

5 Illustrative example: longitudinal model parameterizations

To substantiate the theoretical insights presented in the paper, we applied the multivariate ALT models using both the standard and alternative (1 and 2) parameterizations, where threshold invariance is assumed. The results for illegal drug use (drug), depressive symptoms (depr), and general health status (health) are presented in Table 18. The different fit statistics for these models highlight their distinct specifications.

Table 18 (Significant) parameter estimates for different parameterizations of the ALT model.[ $(^{*})$ : significant at 5% level. $(^{**})$ : significant at 1% levelNote: $(^{***})$ : significant at 0.1% level.]

Significant estimates are selectively displayed in the table for clarity, and parameters are grouped based on their relevance to either the cross-lagged/autoregressive component of ALT or the multivariate latent growth part. Thresholds and error variances are not reported here due to space constraints but are all found to be significantly different from zero. All models demonstrate a comparable fit to the data, with a superior performance observed for the ALT based on the alternative 2 parameterization.

Examining the cross-lagged relationships among the underlying variables, we consistently find that the perception of health status at each wave is influenced by both the prior perception and the propensities of using illegal drugs and being depressed in the preceding period. The propensity for being depressed on a given occasion is influenced by the propensity at the previous occasion and the prior perception of health status. Notably, the standard parameterization indicates a non-significant autoregressive effect in wave 6.

Divergent conclusions emerge regarding the cross-lagged and autoregressive relationship for the drug variable in different parameterizations. While all models detect a significant influence of prior drug use propensity on the current level during the late adolescence period (waves 4, 6, and 8), additional influences of the prior depressive symptoms propensity are found only in the model based on the alternative 2 parameterization proposed by Mehta et al. (2004).

In the adulthood period (waves 10-12-14), the two alternative parameterizations estimate a significant impact of the previous general health status on the propensity of using drugs, aligning with the findings by Silver et al. (2023).

Consistent conclusions are drawn across all parameterizations in terms of concomitant effects. Regarding the associative multivariate growth component of ALT models, both the growth components (intercept and slope) associated with the underlying variables of depressive symptoms and general health status are correlated.

On the other hand, the intercept and slope associated with the propensity to use drugs correlate with the intercept and slope of the health status variable, respectively. None of the growth components specific to each variable correlate with each other, except for the perception of the health status, but not in the alternative 2 parameterization.

In the standard parameterization, the intercept associated with the drug’s underlying variable correlates with the slope of the general health variable; conversely, in the alternative parameterizations, the slope of the drug variable correlates with the intercept of the health status variable.

Hence, these empirical results underscore the critical impact of alternative parameterizations in autoregressive latent trajectory models on the conclusions that researchers can draw when applying these models to real data.

6 Conclusions

In this paper, we have undertaken a thorough examination of how various scaling constraints, implemented to ensure model identification, influence the specification of SEMs used in the analysis of longitudinal categorical data. Our study combines theoretical considerations with empirical validations, providing essential insights into the potential consequences of different parameterization choices.

We focused our attention on ALT models and their special cases: the linear latent growth and first-order autoregressive model. Theoretical investigations have revealed that different parameterizations of the auxiliary model can yield different specifications of the linear latent growth model and of the autoregressive model that, in some cases, are equivalent. Equivalence between alternative specifications has been proven following the approach proposed by Raykov and Penev (1999), such that one-to-one relationships between the parameters of these equivalent specifications have been derived. However, we have shown that when the latent growth model and first-order autoregressive component are jointly considered in the ALT model, different specifications of the auxiliary models imply nonequivalent ALT specifications that are characterized by different constraints on $\boldsymbol \mu $ and $\boldsymbol \Sigma $ .

The implications of the different nonequivalent parameterizations of the ALT model have been evaluated empirically using data from the NLSY97 cohort in examining the relationship between illegal drug use, depressive symptoms, and general health status. By fitting different specifications of the autoregressive latent trajectory model to these data, temporal influences on health perception emerge as a consistent pattern, revealing significant associations between prior perceptions, propensities for drug use, depressive symptoms, and perceived health status.

The examination of the causes affecting the propensity to use illegal drugs reveals sensitivity to parameterization choices since different specifications can provide us with different answers. Using alternative parametrizations showed us that there are extra factors influencing the propensity of using drugs beyond what the standard model suggests. This emphasizes how researchers need to be careful when picking and imposing different identification constraints. We suggest that the alternative parameterization of the auxiliary model proposed by Jöreskog (2001) be more widely adopted, as it closely aligns with the case of observed continuous variables and provides a more coherent approach than what is often used in practice.

Drawing strength from the comprehensive exploration of both theoretical considerations and empirical validations, our study leverages data from the NLSY97 cohort. This dataset provides a robust empirical foundation, aligning seamlessly with the theoretical framework and enhancing the credibility of our findings. Future research could extend these insights by exploring diverse dataset to validate observed patterns across different populations. Additionally, while our focus on specific models provides depth in understanding these structures, further research is needed to explore the applicability of our findings to other SEM configurations.

In conclusion, this paper contributes valuable insights into the complexities of parameterization choices in SEMs for longitudinal categorical data.