Stepwise estimation of latent variable models: An overview of approaches

Abstract

Stepwise approaches for the estimation of latent variable models are becoming increasingly popular, both in the context of models for continuous (factor analysis and latent trait models) and discrete (latent class and latent profile models) latent variables. Examples include two-stage path analysis, structural-after-measurement and Croon’s bias-corrected estimation of structural equation models, and two- and three-step latent class and latent Markov modelling. These methods have in common that the measurement/clustering part of the model is estimated first, followed by the estimation of a—possibly complex—structural model. In this article, we review the existing approaches, which differ in how the information on the latent variable(s) is used when estimating the structural model. We show that based on these differences, stepwise latent variable modelling approaches can be classified into three main types: the fixed parameters, the single indicator and the bias adjustment approach. We discuss similarities and differences between these approaches, as well as between approaches proposed specifically for either continuous or discrete latent variables. Special attention is paid to heterogeneous measurement error resulting from missing data or measurement non-invariance, standard error estimation and software implementations.

Keywords

Classification error Croon’s bias correction latent class analysis measurement error structural equation modelling

1 Introduction

Recently, we have seen a renewed interest in stepwise latent variable modelling approaches, which involve separating the estimation of the measurement and the structural parts of the model of interest. Examples for continuous latent variables include Structural After Measurement (SAM) estimation (Rosseel and Loh, 2024), two-stage path analysis (Lai and Hsiao, 2022), two-step latent trait modelling (Kuha and Bakk, 2023) and Measurement and Uncertainty Preserving ParamETric (MUPPET) modelling (Levy, 2023; Levy and McNeish, 2024). Similar approaches have been proposed for discrete latent variables, such as three-step latent class analysis (Bolck et al., 2004; Vermunt, 2010) and two-step latent class analysis (Bakk and Kuha, 2018).

Typically, two types of arguments are mentioned by authors contributing to the field of stepwise latent variable modelling; that is, arguments related to model estimation and to model building. Regarding model estimation, separate estimation of the measurement and the structural parameters yields fewer convergence problems is less affected by model misspecifications, and is computationally less demanding when dealing with large models (Bartolucci et al., 2015; Perez Alonso et al., 2024; Rosseel and Loh, 2024). In addition, stepwise model building, where one constructs measurements for the variables of interest and subsequently uses these to investigate their relationships, is also what most applied researchers prefer doing, among others to comply with the standard practice in their field, to prevent interpretational confounding, or to circumvent the need to remove cases with missing values on covariates (Burt, 1976; Vermunt, 2010; Clouth et al., 2022; Levy, 2023). Given that most researchers use such a strategy, stepwise latent variable modelling approaches have the potential to yield great improvement over naive use of estimated latent variable scores in subsequent analyses without accounting for their uncertainty (Rein et al., 2025).

As shown in Table 1, stepwise approaches have been developed for each of the four types of latent variables models described by Lazarsfeld and Henry (1968) and Bartholomew and Knott (1999). These are (a) factor analytic models for continuous latent variables and continuous observed indicators, (b) latent trait models for continuous latent variables and discrete observed indicators, (c) latent profile models for discrete latent variables and continuous observed indicators and (d) latent class models for discrete latent variables and discrete observed indicators. Latent trait models are also referred to as item response theory (IRT) models or categorical factor analytic models, and latent profile and latent class models are also referred to as mixture models.

Table 1

Types of stepwise approaches for different measurement models.

	Measurement model
Stepwise approach	Factor analysis	Latent trait and categorical factor analysis	Latent class, latent profile and mixture analysis
Fixed parameters	global SAM (Rosseel and Loh, 2024); MUPPET (Levy, 2023)	global SAM (Rosseel and Loh, 2024); two-step latent trait (Kuha and Bakk, 2024); MUPPET (Levy, 2023)	two-step latent class (Bakk and Kuha, 2018), three-step latent Markov (Bartolucci et al., 2015), two-step multilevel latent class (Di Mari et al., 2023)
Single indicator	two-step path analysis (Lai and Hsiao, 2022), improved regression calibration (Skrondal and Kuha, 2012), stepwise latent vector autoregressive modelling (Rein et al., 2025)	two-step path analysis (Lai and Hsiao, 2022), latent trait measurement error correction (Wang et al., 2019)	three-step ML latent class (Vermunt, 2010; Bakk et al. 2013), three-step ML mixture model clustering (Gudicha and Vermunt, 2013), three-step latent Markov (Di Mari et al., 2018), three-step multilevel latent class (Lyrvall et al., 2024), three-step mixture modelling (Asparouhov and Muthén 2014)
Bias adjustment	local SAM (Rosseel and Loh, 2024), Croon’s correction (Croon, 2002), factor score regression (Hoshino and Bentler, 2013; Hayes and Usami, 2020), factor score path analysis (Devlieger and Rosseel, 2017), multilevel factor score regression (Croon and van Veldhoven, 2007; Devlieger and Rosseel, 2017); C roon’s bias corrected latent interactions (Cox and Kelcey, 2021), mixture multigroup SEM (Perez Alonso, 2024)	local SAM (Rosseel and Loh, 2024), factor score regression (Hoshino and Bentler, 2013), Croon’s correction (this article)	three-step BCH latent class (Bolck et al., 2004; Vermunt, 2010; Bakk and Vermunt, 2016), three-step BCH mixture model clustering (Gudicha and Vermunt, 2013)

The structural models for which stepwise approaches have been developed are sometimes simple linear or logistic regression models in which a single latent variable serves as dependent or independent variable, but can also concern more complex models, such as path models containing multiple latent variables, latent Markov models or dynamic factor models for longitudinal data, or latent variable models for multilevel data.

Stepwise latent variable approaches appear in the literature under different names such as Structural After Measurement (SAM) estimation, two-stage analysis, two-step analysis, Croon’s bias-corrected estimation, Croon’s method and three-step analysis. However, the first step of all these approaches deals with the estimation of the measurement model for the latent variables without accounting for the structural model. The approaches differ from one another with respect to how the measurement model parameters are used when estimating the parameters of the structural model. Based on these differences, stepwise latent variable modelling approaches can be classified into three main types, which we will refer to as the fixed parameters, the single indicator and the bias adjustment approach. Table 1 lists the most important contributions to each of these.

In the fixed parameters approach, when estimating the parameters of the structural model, the measurement model parameters are fixed to the values obtained in the first step. Examples of this approach include Burt’s approach (Burt, 1976), global SAM estimation (Rosseel and Loh, 2024), two-step latent trait analysis (Kuha and Bakk, 2023), MUPPET modelling (Levy, 2023; Levy and McNeish, 2024), two-step latent class analysis (Bakk and Kuha, 2018) and three-step latent Markov modelling (Bartolucci et al., 2015).

The single-indicator approach involves obtaining predictions for the latent variables using the parameters of the measurement model. The predicted scores are used as single indicators when estimating the structural model, while accounting for their unreliability or misclassifications. This is a rather old idea within the SEM framework, and referred to as two-step SEM (Bollen, 1989). Lai and Hsiao (2022) showed how to use this approach when the measurement model is a latent trait model instead of factor analytic model, and refer to this approach as two-stage path analysis. Moreover, Lai et al. (2023) showed how to apply two-stage path analysis when there is measurement non-invariance. In the context of latent class and latent profile analysis, this approach was proposed by Vermunt (2010) and Gudicha and Vermunt (2012), who refer to it as three-step analysis with maximum likelihood adjustment for classification errors (see also Asparouhov and Muthén, 2014). Somewhat related single-indicator approaches have been proposed by Skrondal and Kuha (2012) and Savalei (2019).

The bias adjustment approach involves obtaining predictions for the latent variables using the parameters of the measurement model, and subsequently adjusting their covariances (or associations) so that these represent the true latent variable covariances. The adjusted covariances are used as if they were observed covariances when estimating the structural model parameters. This method was originally developed by Croon (2002) who derived the required adjustment for both factor analytic and latent class models (see also the related work by Wall and Amemiya, 2000). Devlieger and Rosseel (2017) used this approach in what they refer to as factor score path analysis. Although derived and implemented in a slightly different manner, the recently proposed local SAM approach by Rosseel and Loh (2024) can be seen as a special case of Croon’s method for continuous latent variables. Croon’s method for discrete latent variables was extended in various ways by Bolck et al. (2004), Vermunt (2010), Gudicha and Vermunt (2011), and Bakk et al. (2013), and is typically referred to as the Bolck-Croon-Hagenaars (BCH) approach.

The next three sections provide more details on the implementation of the three types of stepwise estimation methods with factor analytic models (Section 2), latent trait and categorical factor analytic models (Section 3), and latent class, latent profile and mixture models (Section 4). Subsequently, in Section 5, we discuss topics such as heterogeneous measurement error resulting from missing data or measurement non-invariance, standard error estimation and software implementations. We conclude with a short discussion in Section 6.

2 Continuous latent and observed variables

2.1 Background

Assuming the data are centred, and therefore ignoring the mean structure and representing the vector of latent variables by the vector η, the regression equations defining the measurement and structural parts of a structural equation model (SEM; Bollen, 1989) for the response vector y are:

\begin{matrix} y = Λ η + ϵ_{i} \\ η = B η + ζ_{i} . \end{matrix}

(2.1)

The free parameters are the factor loadings, residual covariances, regression coefficients and residual factor covariances, which are collected in the matrices Λ, Θ, B and Ψ. These yield the expected covariance matrix Σ as follows:

Σ = Λ V a r (η) Λ^{'} + Θ

with

V a r (η) = {(I - B)}^{- 1} Ψ {(I - B^{'})}^{- 1} .

The various stepwise approaches have in common that they first estimate the measurement parameters Λ and Θ, while the structural parameters B and Ψ are estimated in a next step of the analysis. The approaches differ in how the information on the measurement parameters is used when estimating the structural parameters. Two of the approaches make use of factor scores, either as single indicators or as an intermediate step to obtain an adjusted factor-score covariance matrix.

Let f be the factor scores which are obtained as f = Ay, with A being the factor-score matrix. The most common type of factor score is the regression factor score, with factor-score matrix

A_{R} = Φ Λ' Σ^{- 1},

(2.2)

or, equivalently,

A_{R} = {(Φ^{- 1} + Λ' Θ^{- 1} Λ)}^{- 1} Λ' Θ^{- 1} .

Here Φ denotes the factor covariance matrix obtained from the measurement model. This yields modal a posteriori (MAP) or expected a posteriori (EAP) estimates of the factor scores. The posterior variance of the regression factor scores equals:

V a r_{f} = Φ - Φ' Λ' Σ^{- 1} Λ Φ .

(2.3)

Maximum likelihood or Barlett factor scores (Barlett, 1938) are obtained with factor-score matrix

A_{B} = {(Λ' Θ^{- 1} Λ)}^{- 1} Λ' Θ^{- 1} .

(2.4)

Instead of using a factor-score matrix based on the measurement model parameters, one may also use a simple sum score of the items as an estimate of the latent variable, in which case the factor-score matrix equals A _S = 1.

2.2 Fixed parameters

To prevent what he refers to as interpretational confounding, Burt (1976) proposed estimating the measurement model separately for each latent factor, and subsequently estimating the structural model of interest with the measurement parameters—loadings and residual covariances—fixed. A similar but slightly more general fixed parameters approach was recently proposed by Rosseel and Loh (2024), who refer to it as global SAM estimation. An important difference is that global SAM was presented as an estimation method, whereas Burt presented his approach as a modelling method. In addition, the global SAM approach is more flexible in the sense that the measurement model may be estimated jointly for all latent variables, separately for each latent variable, or separately for blocks of latent variables. Moreover, the estimator used for the measurement model can differ from the one used for the structural model.

Another fixed measurement parameters approach called MUPPET modelling was proposed by Levy (2023), who also stressed the issue of interpretational confounding. MUPPET is a Bayesian two-step MCMC estimation approach for factor analytic models with covariates and outcome variables. The Bayesian estimation framework allows accounting for the uncertainty in the measurement parameters when estimating the structural parameters; that is, the saved draws of measurement model parameters from their posterior distribution are reused when drawing the structural parameters.

2.3 Single indicator

The single-indicator approach has been advocated among others by Lai and Hsiao (2022) who refer to it as two-stage path modelling. Rein et al. (2025) proposed using this approach in the context of latent vector autoregressive models for intensive longitudinal data subject to measurement error.

The single-indicator approach assumes that factor scores f are computed after estimating the measurement model, where A can be any type of factor-score matrix. When estimating the structural parameters, the factor scores are used as indicators for the latent variables η. Using f = Ay and substituting y by the measurement part of the SEM defined in equation 2.1 yields:

f = A Λ η + A ϵ = Λ^{*} η + ϵ^{*},

(2.5)

which shows how the factor scores are related to the true latent variable scores. Moreover, this equation demonstrates factor scores f can serve as single indicators of the latent variables η by using Λ^∗ = AΛ as their ‘loading’ matrix’ and AΘA^′ as the covariance matrix Θ^∗ of their residuals ϵ^∗ = Aϵ. In other words, the structural parameters can be estimated using a SEM with a single indicator per factor, where loadings and residual variances depend on the parameters of the measurement model and the type of factor scores used.

The measurement model for the factor scores can also be written in terms of factor covariances, which yields

V a r (f) = A Λ V a r (η) Λ' A' + A Θ A' .

(2.6)

The reliability of the factor scores equals:

ρ = diag (A Λ V a r (η) Λ' A' V a r {(f)}^{- 1}),

(2.7)

that is, the ratio between the (appropriately weighed) true and the estimated factor variance.

Let a _R be the regression factor score vector from a single factor. A useful feature of regression factor scores obtained from the measurement model for a single factor (or by ignoring the estimated true factor correlations when creating the factor scores) is that their fixed loading and residual variance $a_{R}^{'} λ$ and $a_{R}^{'} Θ a_{R}$ can be computed from the individual factor scores and their variances. Moreover, the fixed loading $a_{R}^{'} λ$ equals the factor score reliability ρ, which in turn equals the ratio between the variance of the factor scores and the true factor variance, var(f)/ϕ. Using the definition of a _R from equation 2.2, the expected variance of regression factor scores can be expressed as follows:

v a r (f) = a_{R}^{'} Σ a_{R} = ϕ λ^{'} Σ^{- 1} {ΣΣ}^{- 1} λ ϕ = ϕ λ^{'} Σ^{- 1} λ ϕ = a_{R}^{'} λ ϕ .

(2.8)

Using equation 2.7 for a single factor and noting that ϕ = var(η), we obtain

ρ = \frac{{(a_{R}^{'} λ)}^{2} ϕ}{v a r (f)} = \frac{{(a_{R}^{'} λ)}^{2} ϕ}{a_{R}^{'} λ ϕ} = a_{R}^{'} λ .

Multiply by ϕ/ϕ yields

ρ = \frac{a_{R}^{'} λ ϕ}{ϕ} = \frac{v a r (f)}{ϕ} .

This shows the reliability of regression factor scores equals the ratio between the estimated and the true factor variance. Using this equation, the factor score variance var(f) can be decomposed as follows:

v a r (f) = ϕ ρ = ϕ ρ ρ + ϕ ρ (1 - ρ),

which is a special case of the more general equation 2.6, with $ρ = a_{R}^{'} λ$ and the $ϕ ρ (1 - ρ) = a_{R}^{'} Θ a_{R},$ where, as above, Θ represents the residual covariance matrix of the indicators. This shows that using the single-indicator approach with regression factor scores involves fixing their loading to ρ and their residual variance to ϕρ (1 − ρ).

Another important characteristic of regression factor scores is that their true factor variances ϕ can be derived from the factor scores and their variances as follows:

ϕ = v a r (f) + v a r_{f} .

(2.9)

where var_f represents the posterior variance of the individual’s factor scores (see equation 2.3). Equation 2.9 follows from the definition of var(f) in equation 2.8 and the application of the definition in equation 2.3 to a single factor, which yields var_f. Note that this equation is equivalent to the Expectation step for obtaining the expected factor covariance matrix when using the Expectation-Maximization algorithm for maximum likelihood factor analysis (Dempster et al., 1977).

In summary, stepwise SEM using regression factor scores as single indicators involves fixing their loadings and residual variances to ρ and ϕρ(1 − ρ) respectively. The quantities ρ and ϕ can be computed from the factor scores and their variances, thus without knowledge of original measurement model parameters. This is a very convenient since factor scores and their variances are standard output provided by software for factor analysis.

With Barlett factor scores, the fixed loading and residual variance equal 1 and $a_{B}^{'} Θ a_{B}$ respectively. When applied to a single factor at a time, Barlett and regression scores are proportional to one another; that is, Barlett scores are 1/ρ times larger than regression scores. This implies their loading is 1/ρ times larger and their residual variance 1/ρ² times larger. This yields a loading of 1 and a residual variance of ϕ(1 − ρ)/ρ.

Instead of using factor scores, the single-indicator approach can also be applied using a simple sum score of the items. In that case, a _S = 1, implying $λ^{*} = a_{S}^{'} λ$ and $θ^{*} = a_{S}^{'} Θ a_{S}$ are simple sums of loadings and residual covariances from the measurement model for the factor concerned. The sum score reliability computed by applying equation 2.6 for a single factor yields the familiar omega reliability for composite scores (Raykov, 1997; McDonald, 1999). Note that this correction for the unreliability of sum scores requires estimating the measurement model. Somewhat related is the work by Oberski and Satorra (2013) and Savalei (2019), who proposed using the single-indicator approach with composite scores with a known reliability ρ. This involves fixing their loading and residual variance to 1 and var(f)(1 − ρ) respectively.

Other related work concerns the improved regression calibration by Skrondal and Kuha (2012), which can be used to correct for covariate measurement error in regression models. It involves using EAP factor scores as predictors while accounting for their uncertainty using their posterior variances. In fact, these authors make use of the fact that the total variance of the true latent factor equals the sum of the variance of the estimated factor scores and their posterior variance (see equation 2.9).

2.4 Bias adjustment

The approach proposed by Croon (2002) involves obtaining factor scores and their covariance matrix Var(f), and subsequently computing the adjusted factor covariance matrix Var(η). This matrix is used as input when estimating the parameters of the structural model. Croon’s method can be best understood by writing the relationship between the factor scores and the true factors defined in equation 2.5 as follows:

A Λ η = f - A ϵ,

or, in terms of covariances, as

A Λ V a r (η) Λ' A' = V a r (f) - A Θ A' .

This shows Var(η) can be obtained from Var(f) as follows:

V a r (η) = {(A Λ)}^{- 1} (V a r (f) - A Θ A') {(Λ' A')}^{- 1} .

Note that Croon (2002) derived this formula in scalar form for each of the elements of the factor covariance matrix Var(η), and, moreover, assumed the measurement models were estimated separately for each factor involved. The more general form given here allows for the simultaneous estimation of the measurement models of all factors, which is required when there are cross-loadings or correlated residuals among items loading on different factors (Hayes and Usami, 2020). As can be seen, the terms AΛ and AΘA^′ are the same as those appearing in the derivation of the single-indicator approach, which shows among others that for regression factor scores, Croon’s correction can be applied without knowledge of the original measurement model parameters when there are no cross-loadings or correlated errors. Note also that with Barlett factor scores, Croon’s correction formula simplifies to

V a r (η) = V a r (f) - A_{B} Θ A_{B}^{'} .

Based on this equation, Hoshino and Bentler (2013) proposed equating the diagonal elements of Var(η) to the estimated factor variances from the first step and the off-diagonal elements to the corresponding entries in Var(f). This works because the off-diagonal elements of $A_{B} Θ A_{B}^{'}$ equal 0 as long as there are no cross-loadings or residual correlations between items loading of different factors. Somewhat related is the work by Skrondal and Laake (2001) who showed that unbiased linear regression coefficients can be obtained without any bias adjustment of Var(f) by using regression factor scores for predictors and Barlett factor scores for outcome variables.

Croon’s bias adjustment has been extended for dealing with various types of more complex structural models, such as for micro-macro multilevel analysis (Croon and van Veldhoven, 2007), multilevel path analysis (Devlieger and Rosseel, 2019), and path models with latent interactions (Cox and Kelcey, 2021; Rosseel et al., 2022).

Although the local SAM approach by Rosseel and Loh (2024) is similar to Croon’s approach, it uses a slightly different starting point. In local SAM, we rewrite equation 2.5 as follows

A Λ η = A (y - ϵ),

where A is referred to as the mapping matrix. It should be chosen such that AΛ = I, which yields

η = A (y - ϵ),

or, in terms of covariances,

V a r (η) = A (V a r (y) - Θ) A^{'} .

(2.10)

The form of A depends on the estimation method used to find Var(η); for example, with maximum likelihood (ML) estimation, it coincides with the Barlett factor-score matrix A _B (see equation 2.6). It can easily be verified that $A_{B} Λ = I$ . Note that local SAM differs from Croon’s approach in that it does not obtain factor scores as an intermediate step, but that instead Var(η) is computed directly from y. As in Croon’s approach, the structural model is estimated using Var(η) as input data.

Recently, Perez Alonso et al. (2024) proposed a rather advanced application of the SAM approach; that is, they used it for estimation of a mixture multigroup SEM. Their measurement model consists of a multiple group factor analysis (for many groups) and their structural model uses the group-specific factor covariances in a mixture model aimed at clustering the groups based on their structural parameters.

Similar to Croon’s method and local SAM, Gerbing and Anderson (1987) and Lance et al. (1988) proposed estimating the structural parameters using Var(η) as data matrix. However, rather than using bias-adjusted factor scores or a mapping matrix, they used the fact Var(η) equals the step-one factor covariance matrix Φ when the measurement models for all latent variables appearing in the structural model are estimated as a single block with an unrestricted factor covariance matrix.

3 Continuous latent and categorical observed variables

3.1 Background

When the latent variables are continuous and the observed variables are ordinal categorical, the measurement model may take on the form of a categorical factor analysis (CatFA). As factor analysis with continuous observed variables, CatFA involves estimating a linear factor model, but with the observed variables’ polychoric correlation matrix instead of their covariance matrix as input data (see, e.g., Forero et al., 2009). Parameters are typically estimated by diagonally weighted least squares (DWLS; Christoffersson, 1977; Muthén et al., 1997).

With continuous latent variables and categorical observed variables, the measurement model may also be a latent trait or IRT model (Bartholomew and Knott, 1999; van der Linden, 2017), typically estimated using marginal maximum likelihood. A unidimensional latent trait model for response pattern y _i can be defined as follows:

P (y_{i}) = \int P (y_{i} | η) f (η) d η = \int \prod_{j = 1}^{J} P (y_{i j} ∣ η) f (η) d η,

where the probability of the response of person i on item j conditional on the latent trait value η, $P (y_{i j} | η)$ , is modeled using a logit or a probit model in which the latent trait enters as predictor. EAP estimates of the latent trait scores and their posterior variances are defined as follows:

\begin{matrix} f_{i} = E (η ∣ y_{i}) = \frac{\int η P (y_{i} ∣ η) f (η) d η}{P (y_{i})} \\ v a r_{f_{i}} = v a r (η ∣ y_{i}) = \frac{\int η^{2} P (y_{i} ∣ η) f (η) d η}{P (y_{i})} - E {(η ∣ y_{i})}^{2} . \end{matrix}

Since there is no closed form expression for the integral, it has to be solved either numerically or by Monte Carlo simulation.

3.2 Fixed parameters

One particular fixed parameters approach for dealing with continuous latent variables and categorical observed variables is the global SAM approach. Although not mentioned explicitly by Rosseel and Loh (2024), their global SAM approach can also be applied in combination with a CatFA for dichotomous or ordinal items. This method is implemented in the lavaan package in R.

Kuha and Bakk (2023) proposed a two-step approach which involves fixing the latent trait measurement parameters defining the conditional response probabilities $P (y_{i j} | η)$ when estimating the structural part of the model. In its most general form, their structural model consisted of a path model for multiple latent variables (with fixed measurement models), with in addition observed predictors affecting the latent variables and observed outcomes affected by the latent variables.

The MUPPET approach for factor analytic models has recently been extended to allow for latent trait type measurement models (Levy and Neish, 2024). MUPPET uses a two-stage Bayesian estimation method for dealing with covariates and outcome variables.

3.3 Single indicator

Lai and Hsiao (2022) showed how the single-indicator approach they refer to as two-stage path analysis can be used with latent variable scores obtained from a latent trait model. They proposed using EAP trait scores, which are IRT based equivalents of regression factor scores. Related work was done by Wang et al. (2019), who applied the single-indicator approach to deal with measurement error in MAP latent trait scores used as the outcome variable in a linear mixed effects model.

As explained above for the continuous indicators case, with regression (and thus EAP) factor scores we can fix the loadings of the factor scores to their reliability ρ and their residual variance to φρ(1 − ρ). However, different from the factor analytic case, the reliability of IRT-based trait scores varies across individuals because the EAP scores’ posterior variance, $v a r_{f_{i}}$ , varies across response patterns. More specifically, Lai and Hsiao (2022) proposed using

λ_{i}^{*} = ρ_{i} = 1 - \frac{v a r_{f_{i}}}{φ}

(3.1)

and

θ_{i}^{*} = φ ρ_{i} (1 - ρ_{i})

as fixed loading and residual variance for person i. Note that, similar to the factor analysis case, φ can be obtained from the individual factor scores and their variances; that is, by using $v a r_{f} = E (v a r_{f_{i}})$ in equation 2.9.

Note that the single-indicator approach using latent trait scores would be much simpler and, moreover, computationally less demanding if we replace $v a r_{f_{i}}$ in equation 3.1 by $v a r_{f} = E (v a r_{f_{i}})$ ; that is, if we ignore the heterogeneity of the measurement error when estimating the structural parameters. Lai and Hsiao (2022) indicated that unmodeled heterogeneity may affect estimation in a path model, but at this time it is unclear how this might occur. Most likely, it has little impact as long as the measurement error heterogeneity is uncorrelated with the other variables appearing in the structural model.

3.4 Bias adjustment

Although not described in the Rosseel and Loh (2024) paper, their local SAM approach is available in the lavaan package for categorical variables. First CatFA models are estimated using DWLS and subsequently equation 2.10 is applied with a Barlett mapping matrix and Var(y) replaced by the polychoric correlation matrix. The structural model is estimated with Var(η) serving as the input data matrix.

Although not proposed in the literature yet, it also seems possible to implement Croon’s approach with EAP trait scores from IRT measurement models. We saw that with regression factor scores $a^{'} λ = λ^{*} = ρ$ and $a^{'} Θ a = θ^{*} = ϕ ρ (1 - ρ)$ . This suggests that the trait score covariance matrix can be corrected using λ^∗ and θ^∗, defined as

\begin{matrix} λ^{*} = ρ = 1 - \frac{E (v a r_{f_{i}})}{φ} \\ θ^{*} = φ ρ (1 - ρ), \end{matrix}

with

φ = v a r (f) + E (v a r_{f_{i}}) .

While such an adjustment does not account for the heterogeneous measurement error, it is also possible to apply the correction with individual-specific $λ_{i}^{*}$ and $θ_{i}^{*}$ ; that is, by applying Croon’s correction formula separately to the individual contributions to the factor-score covariance matrix and taking their average. For a variance, this would be

v a r (η) = E \{{[f_{i} - E (f)]}^{2} / ρ_{i}^{2} - φ_{i} (1 - ρ_{i}) / ρ_{i}\},

and for a covariance between two factors η₁ and η₂

v a r (η_{1}, η_{2}) = E \{[f_{1 i} - E (f_{1})] [f_{2 i} - E (f_{2})] / [ρ_{1 i} ρ_{2 i}]\} .

Hoshino and Bentler (2013) indicated their stepwise procedure using maximum likelihood factor scores can also be applied with IRT type measurement models. As was explained above for factor analytic models, their approach involves equating the diagonal elements of Var(η) to the estimated true factor variances from the first step and the off-diagonal elements to the factor score covariances.

4 Categorical latent variables

4.1 Background

Denoting the single discrete latent variable by η, a particular class by c, and the number of classes by C, a latent class measurement model for P(y _i ), the probability of the response pattern of individual i, can be defined as follows (Goodman, 1974; Hagenaars and McCutcheon, 2002):

P (y_{i}) = \sum_{c = 1}^{C} P (η = c) P (y_{i} | η = c) = \sum_{c = 1}^{C} P (η = c) \prod_{j = 1}^{J} P (y_{i j} | η = c) .

Here $P (η = c), P (y_{i} | η = c)$ , and $P (y_{i j} | η = c)$ represent the class proportion, the probability of the response pattern, and the probability of a particular response on item j for class c. A latent profile or a mixture model for continuous responses has the same structure, but with $P (y_{i}), P (y_{i} | η = c)$ , and $P (y_{i j} | η = c)$ being probability densities instead of probabilities. See Bartolucci et al. (2022) for a recent review of discrete latent variable models.

The above formula is the one for a simple latent class measurement model. The most common extension containing a structural part involves including covariates x _i affecting class membership:

P (y_{i} | x_{i}) = \sum_{c = 1}^{C} P (η = c | x_{i}) \prod_{j = 1}^{J} P (y_{i j} | η = c) .

Other common extensions containing a structural part are models with a distal outcome z_i affected by class membership,

P (y_{i}) = \sum_{c = 1}^{C} P (η_{1} = c_{1}) P (η_{2} = c_{2} | η_{1} = c_{1}) \prod_{j = 1}^{J_{1}} P (y_{i j} | η_{1} = c_{1}) \prod_{k = J_{1} + 1}^{J_{1} + J_{2}} P (y_{i j} ∣ η_{2} = c_{2}),

as well as models with multiple latent variables affecting each other, for example, η₁ affecting η₂,

P (y_{i}) = \sum_{c = 1}^{C} P (η_{1} = c_{1}) P (η_{2} = c_{2} | η_{1} = c_{1}) \prod_{j = 1}^{J_{1}} P (y_{i j} | η_{1} = c_{1}) \prod_{k = J_{1} + 1}^{J_{1} + J_{2}} P (y_{i j} | η_{2} = c_{2}),

The structural parameters of interest are the regression coefficients defining $P (η = c | x_{i}), P (z_{i j} | η = c)$ and $P (η_{2} = c_{2} | η_{1} = c_{1})$ , respectively, which as explained next can be estimated in a stepwise manner.

Latent class predictions f_i based on the latent class measurement model can be obtained based on the posterior membership probabilities $P (η = c | y_{i})$ , which are obtained as follows:

P (η = c | y_{i}) = \frac{P (η = c) P (y_{i} | η = c)}{P (y_{i})} .

The standard practice is to assign individuals to the class with the largest posterior probability, which yields MAP estimates. Other options are random assignment or proportional assignment (Bolck et al., 2004). As shown by Bolck et al. (2004), the assignment method used defines the values of $P (f = d | y_{i})$ , the probability of being assigned class d conditional on the response pattern. Modal assignment gives a value of 1 for the class with the largest posterior probability and 0 for the other classes, whereas random and proportional assignment give values equal to the posterior probabilities.

4.2 Fixed parameters

Bakk and Kuha (2018) proposed a two-step latent class analysis approach, which involves fixing latent class measurement parameters defining the conditional response probabilities $P (y_{i j} | η = c)$ when estimating the structural model. In fact, one fixes $P (y_{i} | η = c)$ , the probability of a full response pattern given class membership to its estimate from the first step. The structural models considered by Bakk and Kuha were regression models in which a single latent variable serves as either predictor or outcome variable. Xue and Bandeen-Roche (2002) proposed the same approach for a latent variable serving as the outcome variable, but motivated by applications where the first step was based on a much larger sample than the second step.

Applications of the fixed parameters approach with more complex structural models concern the three-step latent Markov model by Bartolucci et al. (2015) and the two-step multilevel latent class model by Di Mari et al. (2023). In both applications, the first step involves estimating the measurement model parameters, while ignoring the (longitudinal or multilevel) dependence structure in the data.

Although not yet mentioned in the literature, the fixed measurement parameters approach can be used without any modification with continuous response variables, that is, with latent profile and other types of mixture models. In that case, the fixed parameters in the estimation of the structural model of interest are the class-specific means and (co)variances defining the multivariate densities $f (y_{i} | η = c)$ .

4.3 Single indicator

For latent class models with covariates, Bolck et al. (2004) derived the relationship between $P (η = c | x_{i})$ and $P (f_{i} = d | x_{i})$ ; that is, between the true class distribution given covariates x _i and the corresponding assigned class distribution given x _i . It has the following form:

P (f_{i} = d | x_{i}) = \sum_{c = 1}^{C} P (η = c | x_{i}) P (f_{i} = d | η = c) .

(4.1)

As noted by Vermunt (2010), this is the equation of a latent class model with covariates x _i , where the class assignments f_i serve as the single indicator with $P (f_{i} = d | η = c)$ defining its quality (the amount of classification errors). Vermunt (2010) proposed estimating this single-indicator model by maximum likelihood treating $P (f_{i} = d | η = c)$ as fixed parameters, and called it the three-step ML method. It uses another important result from Bolck et al. (2004); that is, $P (f_{i} = d | η = c)$ can be obtained from the measurement model parameters as follows:

P (f = d | η = c) = \sum_{y} P (y | η = c) P (f = d | y) = \frac{\sum_{y} P (y) P (η = c | y) P (f = d | y)}{P (η = c)} .

Vermunt (2010) proposed a slight modification of this formula, which involves replacing P(y) by its empirical distribution in the dataset at hand. This yields:

P (f = d | η = c) = \frac{\sum_{i = 1}^{N} \frac{1}{N} P (η = c | y_{i}) P (f = d | y_{i})}{P (η = c)} .

(4.2)

Important advantages of this modification are: (a) it prevents the possibility of summing over a very large numbers of data patterns; (b) as shown by Gudicha and Vermunt (2013), it also works with latent profile and other mixture models for continuous observed variables, where the summation over all possible y takes the form of a multivariate integral; (c) irrespective of the type of classification rule, it allows computing the classification error probabilities $P (f = d | η = c)$ using the posteriors $P (η = c | y_{i})$ of the N persons in the data file, which is standard output provided by latent class analysis software.

As shown by Bakk et al. (2013), the three-step ML approach can not only be used for latent class models with covariates, but also for latent class models with a distal outcome and SEM-like models with multiple discrete latent variables. Moreover, this stepwise modelling approach has been used with latent Markov models (Di Mari et al., 2016; Vogelsmeier et al., 2023a), multilevel latent class models (Lyrvall et al., 2024), latent class tree models (van den Bergh and Vermunt, 2019), and inverse propensity weighting for the estimation of the causal effect of a treatment on latent class membership (Clouth et al., 2022, 2023).

A similar stepwise approach for investigating the relationship between class membership and a distal outcome z_i was proposed by Lanza et al. (2013). They proposed including z_i as covariate in the step-one latent class model, which gives posterior class membership probabilities containing information on its association with the latent classes. For the structural model estimation, they use in fact a single-indicator approach with proportional class assignment, but without the need to correct classification errors.

4.4 Bias adjustment

The BCH method (Bolck et al., 2004) is a bias adjustment approach for latent class and latent profile models. Similar to the single-indicator approach, it is based on equation 4.1 which shows the relationship between $P (η = c | x_{i})$ and $P (f_{i} = d | x_{i})$ . What the BCH approach does is to adjust $P (f_{i} = d | x_{i})$ so that it becomes $P (η_{i} = c | x_{i})$ , which is the same as what Croon’s method does for continuous latent variables. The derivation of the adjustment by Bolck et al. (2004) starts by writing equation 4.1 in matrix form: F = HD, with F, H and D containing the $P (f_{i} = d | x_{i}), P (η = c | x_{i})$ , and $P (f = d | η = c)$ respectively. It can now easily be seen that H = FD⁻¹. This shows the adjustment involves multiplying the $P (f_{i} = d | x_{i})$ with the inverse of the matrix with elements $P (f = d | η = c)$ .

As shown by Vermunt (2010), the BCH approach can be conceptualized as an analysis of an expanded dataset with C records per person with weight w_ic. Denoting the elements of D⁻¹ by g_dc, the weight w_ic is obtained as follows:

w_{i c} = \sum_{d = 1}^{C} P (f = d | y_{i}) g_{d c} .

Note that the class assignment probabilities $P (f = d | y_{i})$ enter not only in this formula, but also in the formula for $P (f = d | η = c)$ from equation 4.2. Vermunt (2010) proposed estimating the structural parameters – the logistic regression coefficients defining $P (η_{i} = c | x_{i})$ – by pseudo maximum likelihood using weights w_ic, and using cluster robust standard errors to account for the multiple records per person and the weighting.

The three-step BCH approach can also be used with distal outcome models (Bakk et al., 2013). In fact, it has been shown to be superior to the three-step ML approach when the outcome variable is a continuous variable (Bakk and Vermunt, 2016). Bolck et al. (2004) and Bakk et al. (2013) showed how to apply this approach with multiple latent variables. Lê et al. (2025) proposed using the BCH approach together with inverse propensity weighting for the estimation of the causal effect of latent class membership on an outcome variable.

5 Remaining topics

In this section, we discuss various remaining topics which are relevant for the practical application of stepwise latent variable modelling approaches. These include heterogeneous measurement or classification error resulting from measurement non-invariance and missing data on the latent variable indicators, correcting standard error estimates for the stepwise modelling, and software implementing stepwise modelling approaches.

5.1 Heterogeneous measurement or classification error

When discussing the implementation of the single-indicator approach with IRT type measurement models (Lai and Hsiao, 2022), we already touched upon the topic of heterogeneous measurement error. Our tentative conclusion was that ignoring the heterogeneous nature of the measurement errors probably has little impact on the structural parameter estimates as long as it is uncorrelated with the other variables appearing in the structural model.

Another situation in which measurement or classification errors are heterogeneous occurs when there is missing data on the response variables y. More specifically, each possible missing data pattern has its own λ^∗ and Θ^∗ values or its own $P (f_{i} = d | η = c)$ values, where the quality of the estimated latent variable scores decreases when the proportion of missing values increases. Alagöz and Vermunt (2022) investigated whether one can ignore the classification error heterogeneity caused by missing values in stepwise latent class modelling. Based on a simulation study, their conclusion was that classification error heterogeneity can indeed be ignored as long as a missing at random (MAR) mechanism holds for the y variables concerned. An example of a situation in which the MAR assumption is violated occurs when, conditional of the observed y variables, covariates (which are not used in the measurement model) affect the probability of having missing values. For the continuous latent variable case, the effect of ignoring measurement error heterogeneity caused by missing data has not been investigated. But most likely the same condition holds here; that is, it can be ignored as long as we can assume the missing data is MAR when estimating the step-one model.

Heterogeneous measurement and classification errors also occurs when the step-one model accounts for measurement non-invariance or differential item functioning (van der Linden, 2017). In that case, the terms λ^∗ and Θ^∗ or the $P (f_{i} = d | η = c)$ used in the single indicator and bias adjustment approaches vary across groups with differences in their measurement model parameters. Lai et al. (2023) proposed a two-stage multiple-group path analysis in which the measurement model is a multiple-group factor analysis model and in which group-specific values are used for λ^∗ and Θ^∗ when estimating the multiple-group path model. Vermunt and Magidson (2021) looked into the latent class analysis case with measurement non-invariance (see also Clouth et al., 2023). Their conclusion was that group-specific classification errors $P (f_{i} = d | η = c)$ should be accounted for when the grouping variable is itself part of the structural model. In other words, if the grouping variable is not used in the structural model, the heterogeneous classification errors can be ignored. Most likely, this same conclusion applies in the continuous latent variable case.

It should be noted that whether to account for or to ignore heterogeneous measurement or classification error is not an issue in the fixed parameters approach, since it does not involve estimation of the latent variables scores as an intermediate step. Instead, one fixes the measurement parameters, which may differ per group if one accounted for measurement non-invariance in the first step.

5.2 Standard error computation

In the fixed parameters and the single-indicator approach, we treat either the estimated parameters from the measurement model or functions of these as known when estimating the structural model parameters. As shown by Oberski and Satorra (2013) and Rosseel and Loh (2024) for the continuous latent variable case and by Bakk et al. (2014) and Bakk and Kuha (2019) for the discrete latent variable case, standard errors can be corrected for the fact that these are in fact estimates with their own sampling variation when the structural parameters are estimated by maximizing a log likelihood, which known as pseudo-maximum likelihood estimation (Gong and Samaniego, 1981).

Let Σ _S contain the uncorrected covariances of the parameters of the structural model, Σ _M the covariances of the measurement parameters or, in the single-indicator case, of functions of these, H the second derivatives towards the free parameters, and C towards the free parameters with respect to the fixed parameters. The adjusted covariance matrix of the structural parameters, $Σ_{S}^{*}$ , is obtained as follows:

Σ_{S}^{*} = Σ_{S} + H^{- 1} C Σ_{M} C^{'} H^{- 1} .

In the single-indicator approach Σ_M can be obtained from the covariances of the measurement parameters using the delta method. Bakk et al. (2014) illustrated this for the latent class analysis case; that is, for the three-step ML procedure in which Σ_M contains the covariances of the log $P (f_{i} = d | η = c) / P (f_{i} = c | η = c)$ terms. It should be noted that when dealing with dependent observations such as in multilevel latent class models (Lyrvall et al., 2024), latent Markov models (Di Mari et al., 2016; Vogelsmeier et al., 2023a), and latent vector autoregressive models (Rein et al., 2025), Σ_M will be a cluster-robust parameter covariance matrix.

Similar types of standard error corrections are not available for the bias adjustment approach. This is why, for example, the lavaan package reports global SAM based standard errors when using the local SAM procedure. However, when applying the BCH approach for discrete latent variables, one should account for the fact that the expanded dataset contains multiple records per person and, moreover, BCH weights. As proposed by Vermunt (2010), this can be dealt with using complex sampling (or cluster-robust) standard errors. Robust standard errors were also proposed by Bartolucci et al. (2015) for their three-step latent Markov approach.

As an alternative to the frequentist approach of fixing measurement model parameters to their first-stage point estimates and subsequently adjusting the standard error estimates, one may use Bayesian multiple-stage MCMC estimation, as done in MUPPET modelling. This involves dealing with parameter uncertainty by accounting for the full first-stage posterior parameter distribution when estimating the structural parameters (Levy, 2023; Levy and McNeish, 2024). Similar Bayesian estimation methods could be developed for single indicator and bias adjustment approaches. For example, in the single-indicator approach, instead of using a single set of latent variable scores obtained with point estimates of the measurement parameters, one could use multiple sets of latent variable scores drawn from their posterior distribution while accounting for the uncertainty of the measurement parameters.

5.3 Software implementations

Most of the stepwise approaches we discussed are available in latent variable modelling software. However, with some additional some effort, each of the approaches can be implemented without specific stepwise modelling routines. The fixed measurement and single-indicator approaches can be implemented in SEM or latent class analysis software that allows imposing fixed value constraints on the model parameters. Slightly more tedious to implement is the heterogeneous measurement case resulting from an IRT type measurement model, for which Lai and Hsiao (2022) provided example code for OpenMx (Neale et al., 2016). The bias adjustment approach for continuous latent variables involves creating an adjusted factor covariance matrix, which can subsequently be used as input data for any structural modelling software. For categorical latent variables, one should create an expanded dataset containing the BCH weights, after which any type of structural model can be estimated using a weighted analysis. The only requirement is that the routine used to estimate the structural model accepts negative weights.

LatentGOLD (Vermunt and Magidson, 2016, 2025) is one of the programs that contains options for stepwise latent variable modelling. For latent class models, it implements each of the three approaches, which it refers to as the Bakk-Kuha, ML, and BCH adjustment methods. In version 6.1, it also implements the single-indicator approach for continuous latent variables, where the step-one model can be either a factor analysis or latent trait analysis. After estimation of the measurement model of interest, one saves the posterior probabilities or the logdensities of the latent class model, or the EAP factor scores and their standard errors from the factor analytic or IRT model to an output data file. This information suffices to set up the structural model, where it is possible to indicate that the measurement or classification errors are heterogeneous across levels of a grouping variable. For latent class models, the most common types of structural models (including latent Markov models) are available via the Step3 point-and-click module. Other stepwise latent class models, such as path models for multiple discrete latent variables and multilevel latent class models, can be specified via the syntax system. The structural model for continuous latent variables should always be specified with LatentGOLD’s syntax system, where not only a simple regression or path model can be specified, but also a more complex model, such as a mixture multiple group SEM (Perez Alonso et al., 2024) or a dynamic factor model (Rein et al., 2025).

Mplus (Asparouhov and Muthén, 2014) implements the ML (single indicator) and BCH (bias adjustment) approaches for latent class models in an automated form for simple structural models (covariates and distal outcomes). However, it also allows saving BCH weights to an output file and using these weights in subsequent analysis for the estimation of more complex structural models with latent classes. The BCH approach is also available as SAS (SAS Institute Inc., 2010) and Stata (Statacorp, 2025) functions for dealing with covariates and distal outcomes (Dziak et al., 2017, 2022).

R packages (R Core Team, 2024) implementing stepwise latent variable models include lavaan (Rosseel, 2012), multilevLCA (Lyrvall et al., 2023), lmfa (Vogelsmeier et al., 2023b) and tidySEM (van Lissa, 2019). Lavaan implements the fixed parameters and bias adjustment approaches for factor analysis and CatFA measurement models, and refers to these as global and local SAM. The multilevLCA package implements standard and multilevel latent class analysis with covariates for binary response variables estimated using a two-step (or fixed parameters) approach. The lmfa package implements the single-indicator approach for a continuous-time latent Markov model with covariates for intensive longitudinal data, where the first step is a mixture factor analysis. Finally, tidySEM contains a BCH function for estimating structural models in which the latent classes are treated as the grouping variable.

In addition to full R packages, various recent publications include R code which allows researchers to apply the newly proposed method concerned with their own data. Examples include two-step path modelling (Lai and Hsiao, 2022), MUPPET (Levy, 2023; Levy and McNeish, 2024), mixture multiple-group SEM (Perez Alonso et al., 2024) and stepwise latent vector autoregressive modelling (Rein et al., 2025).

6 Discussion

An overview of stepwise latent variable approaches was provided. It was shown that similar approaches have been proposed for factor analytic, latent trait, and latent class type measurement models. These involve using fixed measurement parameters, latent variable predictions as single indicators, or Croon’s or similar bias adjustments of the latent variable predictions. We explained the logic underlying these approaches with the appropriate references, as well as references to applications in combination with more complex structural models. Moreover, we explained how Croon’s method can also be implemented with IRT-based measurement models.

Special attention was paid to the issue of heterogeneous measurement and classification errors, which occurs with IRT models, with missing values on the latent variable indicators, and with measurement non-invariance. When a grouping variable causing the measurement invariance is also used in the structural model, the heterogeneity should clearly be taken into account. In most other situations, it appears that the heterogeneity can be ignored, but this is something that needs further investigation.

We showed how to obtain standard errors corrected for the stepwise estimation. However, it has also been reported that this correction is not needed or that it may even yield too conservative tests if the measurement or classification errors are not large, and the sample is not small. This issue is also a topic requiring further investigation. Moreover, work remains to be done for the bias adjustment approaches for which standard error correction is not yet available. Interesting is the Bayesian MUPPET modelling approach, which accounts for the stepwise estimation in a rather straightforward manner. Similar Bayesian estimation methods may be developed for application in conjunction with single indicator or bias adjustment approaches. Another popular approach for improved standard error estimation in stepwise modelling is bootstrapping, which surprisingly has not yet been investigated in the context of stepwise latent variable modelling.

We discussed software implementations of the stepwise approach, which include LatentGOLD, Mplus, packages and code written for R, as well as functions written for Stata and SAS. Given that this is a lively field of research, surely, more software for stepwise latent variable modelling will become available in the near future. An important application of the stepwise estimation methods discussed in this article is in big data analytics. Significant computational efficiency gains can be achieved using these stepwise approaches, provided the available software is capable of handling very large datasets.

Footnotes

Declaration of Conflicting Interests

Jeroen Vermunt is a co-developer (along with Jay Magidson) of the LatentGOLD software. While LatentGOLD is a commercial program, as of January 2025, it is free for academic use.

Funding

The authors received no financial support for the research, authorship and/or publication of this article.

References

Alagöz

ÖEC

and Vermunt

(2022) Stepwise latent class analysis in the presence of missing values on the class indicators. Structural Equation Modeling , 29, 784–90.

Asparouhov

and Muthén

(2014) Auxiliary variables in mixture modelling: Three-step approaches using Mplus. Structural Equation Modelling , 21, 329341. doi: 10.1080/10705511.2014.915181.

Bakk

and Kuha

(2018) Two-step estimation of models between latent classes and external variables. Psychometrika , 83, 871–92. doi: 10.1007/s11336-017-9592-7.

Bakk

and Kuha

(2021) Relating latent class membership to external variables: An overview. British Journal of Mathematical and Statistical Psychology , 74, 340–62. doi: 10.1111/bmsp.12227.

Bakk

, Oberski

and Vermunt

(2014) Relating latent class assignments to external variables: Standard errors for correct inference. Political Analysis , 22, 520–40. doi: 10.1093/pan/mpu003.

Bakk

, Tekle

and Vermunt

(2013) Estimating the association between latent class membership and external variables using bias-adjusted three-step approaches. Sociological Methodology , 43, 272–311. doi: 10.1177/0081175012470644.

Bakk

and Vermunt

(2016) Robustness of stepwise latent class modelling with continuous distal outcomes. Structural Equation Modelling , 23, 20–31.

Bartlett

(1938) Methods of estimating mental factors. Natue , 141, 609–10.

Bartholomew

and Knott

(1999) Latent Variable Models and Factor Analysis . London: Arnold.

10.

Bartolucci

, Montanari

and Pandolfi

(2015) Three-step estimation of latent Markov models with covariates. Computational Statistics and Data Analysis , 83, 287–301.

11.

Bartolucci

, Pandolfi

and Pennoni

(2022) Discrete latent variable models. Annual Review of Statistics and Its Application , 9, 425–52.

12.

Bolck

, Croon

and Hagenaars

(2004) Estimating latent structure models with categorical variables: One-step versus three-step estimators. Political Analysis , 12, 3–27.

13.

Bollen

(1989) Structural Equations with Latent Variables . New York: Wiley.

14.

Burt

(1976) Interpretational confounding of unobserved variables in structural equation models. Sociological Methods and Research, 5, 3–52. doi: 10.1177/004912417600500101.

15.

Christoffersson

(1977) Two-step weighted least squares factor analysis of dichotomized variables. Psychometrika , 40, 433–38.

16.

Clouth

, Pauws

, Mols

and Vermunt

(2022) A new three-step method for using inverse propensity weighting with latent class analysis. Advances in Classification and Data Analysis , 16, 351–71.

17.

Clouth

, Pauws

and Vermunt

(2023) Three-step latent class analysis with inverse propensity weighting in the presence of differential item functioning. Structural Equation Modelling , 30, 737–48.

18.

Cox

and Kelcey

(2021) Croon’s bias corrected estimation of latent interactions. Structural Equation Modelling , 28, 863–74.

19.

Croon

(2002) Using predicted latent scores in general latent structure models. In Latent Variable and Latent Structure Models , edited by Marcoulides

and Moustaki

Pages 195–224. Lawrence Erlbaum.

20.

Croon

and van Veldhoven

(2007) Predicting group-level outcome variables from variables measured at the individual level: A latent variable multi-level model. Psychological Methods , 12, 45–57. doi: 10.1037/1082-989X.12.1.45.

21.

Dempster

, Laird

and Rubin

(1977) Maximum likelihood from incomplete data via the EM algorithm (with discussion). Journal of the Royal Statistical Society, Series B , 39, 1–22.

22.

Devlieger

and Rosseel

(2017) Factor score path analysis: An alternative for SEM? Methodology , 13, 31–8. doi: 10.1027/1614-2241/a000130.

23.

Devlieger

and Rosseel

(2019) Multilevel factor score regression. Multivariate Behavioral Research , 55, 600–24. doi: 10.1080/00273171.2019.1661817.

24.

Di Mari

, Bakk

, Oser

and Kuha

(2023) A two-step estimator for multilevel latent class analysis with covariates. Psychometrika , 88, 1144–70. doi: 10.1007/s11336-023-09929-2.

25.

Di Mari

, Oberski

and Vermunt

(2016) Bias-adjusted three-step latent Markov modelling with covariates. Structural Equation Modelling , 23, 649–60. doi: 10.1080/10705511.2016.1191015.

26.

Dziak

, Bray

and Wagner

(2017) LCA_Distal_BCH SAS Macro Users’ Guide (Version 1.1) . University Park, PA: The Methodology Center, Penn State. Retrieved from http://methodology.psu.edu

27.

Dziak

, Bray

and Wagner

(2020) LCA_Covariates_3Step SAS Macro Users’ Guide (Version 1.0) . University Park, PA: The Methodology Center, Penn State. Retrieved from http://methodology.psu.edu

28.

Forero

, Maydeu-Olivares

and Gallardo-Pujol

(2009) Factor analysis with ordinal indicator: A Monte Carlo study comparing DWLS and ULS estimation. Structural Equation Modelling , 16, 625–41.

29.

Gerbing

and Anderson

(1987) Improper solutions in the analysis of covariance structures: Their interpretability and a comparison of alternate respecifications. Psychometrika , 52, 99–111. doi: 10.1007/BF02293958.

30.

Gong

and Samaniego

(1981). Pseudo maximum likelihood estimation: Theory and applications. Annals of Statistics , 9, 861–69.

31.

Goodman

(1974) Exploratory latent structure analysis using both identifiable and unidentifiable models. Biometrika , 61, 215–31.

32.

Gudicha

and Vermunt

(2013) Mixture model clustering with covariates using adjusted three-step approaches. In Algorithms from and for Nature and Life; Studies in Classification, Data Analysis, and Knowledge Organization , edited by Lausen

, van den Poel

and Ultsch

Pages 87–93. Heidelberg: Springer-Verlag GmbH.

33.

Hagenaars

and McCutcheon

(2002) Applied Latent Class Analysis . Cambridge, UK: Cambridge University Press.

34.

Hayes

and Usami

(2020) Factor score regression in the presence of correlated unique factors. Educational and Psychological Measurement, 80, 5–40. doi: 10.1177/0013164419854492.

35.

Hoshino

and Bentler

(2013) Bias in factor score regression and a simple solution. In Analysis of Mixed Data: Methods and Applications , edited by de Leon

and Chough

KC.

Pages 43–61. Boca Raton, FL: Chapman and Hall. doi: 10.1201/b14571-5.

36.

Kuha

and Bakk

(2023) Two-step estimation of latent trait models. Retrieved at https://arxiv.org/pdf/2303.16101.pdf

37.

Lai

MHC

and Hsiao

Y-Y

(2022) Two-stage path analysis with definition variables: An alternative framework to account for measurement error. Psychological Methods , 27, 568–88. doi: 10.1037/met0000410.

38.

Lai

MHC

, WW-Y

Tse

, Zhang

, Li

and Y-Y

Hsiao

(2023) Correcting for unreliability and partial invariance: A two-stage path analysis approach. Structural Equation Modelling , 30, 258–71. doi: 10.1080/10705511.2022.2125397.

39.

Lance

, Cornwell

and Mulaik

(1988) Limited information parameter estimates for latent or mixed manifest and latent variable models. Multivariate Behavioral Research , 23, 171–87. doi: 10.1207/s15327906mbr2302_3.

40.

Lanza

, Tan

and Bray

(2013) Latent class analysis with distal outcomes: A flexible model-based approach. Structural Equation Modelling , 20, 1–26.

41.

Lê

, Clouth

and Vermunt

(2025) Causal latent class analysis with distal outcomes: A modified three-step method using inverse propensity weighting. Multivariate Behavioral Research , 60, 30–60.

42.

Levy

(2023) Precluding interpretational confounding in factor analysis with a covariate or outcome via measurement and uncertainty preserving parametric modelling. Structural Equation Modelling , 30, 719–36. doi: 10.1080/10705511.2022.2154214.

43.

Levy

and McNeish

(2024) Measurement and uncertainty preserving parametric modelling for continuous latent variables with discrete indicators and external variables. Journal of Educational and Behavioral Statistics . Advance online publication.

44.

Lyrvall

, Di Mari

, Bakk

, Oser

and Kuha

(2023) multilevLCA: An R package for single-level and multilevel latent class analysis with covariates. https://arxiv.org/abs/2305.07276

45.

Lyrvall

, Bakk

, Oser

and Di Mari

(2024) Bias-adjusted three-step multilevel latent class modelling with covariates, Structural Equation Modelling , 31, 592–603. doi: 10.1080/10705511.2023.2300087.

46.

McDonald

(1999) Test Theory: A Unified Treatment . Hillsdale, NJ: Erlbaum.

47.

Muthén

, du Toit

SHC

and Spisic

(1997) Robust inference using weighted least squares and quadratic estimating equations in latent variable modelling with categorical and continuous outcomes. Retrieved at https://www.statmodel.com/download/Article_075.pdf

48.

Neale

, Hunter

, Pritikin

, Zahery

, Brick

, Kirkpatrick

, Estabrook

, Bates

, Maes

and Boker

(2016) OpenMx 2.0: Extended structural equation and statistical modelling. Psychometrika , 81, 535–49. doi: 10.1007/s11336-014-9435-8.

49.

Oberski

and Satorra

(2013) Measurement error models with uncertainty about the Error variance. Structural Equation Modelling , 20, 409–28. doi: 10.1080/10705511.2013.797820.

50.

Perez Alonso

, Rosseel

, Vermunt

and De Roover

(2024) Mixture multigroup structural equation modelling: A novel method for comparing structural relations across many groups. Psychological Methods . Advance online publication.

51.

Raykov

(1997) Estimation of composite reliability for congeneric measures. Applied Psychological Measurement , 21, 173–84. doi: 10.1177/01466216970212006.

52.

R Core Team (2024) R: A Language and Environment for Statistical Computing . R Foundation for Statistical Computing, Vienna, Austria. URL https://www.Rproject.org

53.

Rein

, Vermunt

, De Roover

and Vogelsmeier

LVDE

(2025) Latent vector autoregressive modelling: A stepwise estimation approach. Structural Equations Modelling , 32, 105–16.

54.

Rosseel

(2012) lavaan: An R package for structural equation modelling. Journal of Statistical Software , 48, 1–6. doi: 10.18637/ jss.v048.i02.

55.

Rosseel

, Burghgraeve

, Loh

and Schermelleh-Engel

KWW

(2022) Structural after Measurement (SAM) approaches for accommodating latent quadratic and interaction effects . Unpublished manuscript.

56.

Rosseel

and Loh

(2024) A structural after measurement approach to structural equation modelling. Psychological Methods , 29, 561–88. doi: 10.1037/met0000503.

57.

Institute Inc

SAS

(2010) STAT-SAS, Version 9.4 . SAS Institute, Cary, North Carolina.

58.

Savalei

(2019) A comparison of several approaches for controlling measurement error in small samples. Psychological Methods , 24, 352–70.

59.

Skrondal

and Laake

(2001) Regression among factor scores. Psychometrika , 66, 563–75. doi: 10.1007/BF02296196.

60.

Skrondal

and Kuha

(2012). Improved regression calibration. Psychometrika , 77, 649–69.

61.

StataCorp (2025) Stata Statistical Software: Release 19 . College Station, TX: StataCorp LLC.

62.

Van den Bergh

and Vermunt

(2019) Latent class trees with the three-step approach. Structural Equation Modelling , 26, 481–92.

63.

van der Linden

(2017) Handbook of Item Response Theory, Volume One: Models . Boca Raton: Chapman and Hall/CRC.

64.

Van Lissa

(2019) tidySEM: Tidy Structural Equation Modelling. R Package Version 0.2.1 . https://github.com/cjvanlissa/tidySEM/

65.

Vermunt

(2010) Latent class modelling with covariates: Two improved three-step approaches. Political Analysis , 18, 450–69.

66.

Vermunt

and Magidson

(2016) Technical Guide for Latent GOLD 5.1: Basic, Advanced, and Syntax . Belmont MA: Statistical Innovations Inc. Retrieved from www.statisticalinnovations.com

67.

Vermunt

and Magidson

(2021) How to perform three-step latent class analysis in the presence of measurement non-invariance or differential item functioning. Structural Equation Modelling , 28, 356–64.

68.

Vermunt

and Magidson

(2025) Upgrade Manual for Latent GOLD Basic, Advanced, Syntax, and Choice Versions 6.0 and 6.1 . Tilburg, NL: Statistical Innovations Europe. Retrieved from www.statisticalinnovations.com

69.

Vogelsmeier

LVDE

, Vermunt

, Bülow

and De Roover

(2023a) Evaluating covariate effects on ESM measurement model changes with latent Markov factor analysis: A three-step approach. Multivariate Behavioral Research , 58, 262–91.

70.

Vogelsmeier

LVDE

, Vermunt

and De Roover

(2023b) How to explore within-person and between-person measurement model differences in intensive longitudinal data with the R package lmfa. Behavior Research Methods , 55, 2387–2422.

71.

Wall

and Amemiya

(2000) Estimation for polynomial structural equation models. Journal of the American Statistical Association , 95, 929–40.

72.

Wang

, Xu

and Zhang

(2019) Correction for item response theory latent trait measurement error in linear mixed effects models. Psychometrika , 84, 673–700.

73.

Xue

and Bandeen-Roche

(2002) Combining complete multivariate outcomes with incomplete covariate information: A latent class approach. Biometrics , 58, 110–20.