tseLCA Workflow
Overview
tseLCA implements the BCH and ML bias-adjusted
three-step estimators for latent class analysis (LCA) with covariates
and distal outcomes, following the methodological framework for both BCH
and Vermunt’s ML approaches from Bakk, Tekle & Vermunt (2013).
tseLCA also builds on top of the two-step LCA estimation
procedure outlined by Bakk & Kuha (2018), and using the R package
multilevLCA for efficient measurement model estimation from
Lyrvall et al. (2025). tseLCA provides analytic sandwich
variance estimation that propagates measurement uncertainty through the
classification-error correction in the final step.
The three-step approach separates the model into:
- Step 1 — Estimate the LCA measurement model
(class-conditional item probabilities and class prevalences).
- Step 2 — Assign posterior class probabilities and
compute the misclassification matrix.
- Step 3 — Estimate the structural model (covariate
effects or distal outcome means) using the bias-adjusted weights.
Synthetic data
The built-in data-generating process replicates the design of Bakk
& Kuha (2018). Each dataset has six binary indicators (\(Y_1, \ldots, Y_6\)) drawn from a
three-class LCA, plus either a covariate \(Z_p
\sim \text{Uniform}\{1,\ldots,5\}\) predicting class membership,
or a continuous distal outcome \(Z_o\)
predicted by class membership.
# High separation: P(Y_h = 1 | class) = 0.9 / 0.1
d <- generate_data(
n = 500,
separation = "high",
scenario = "covariate",
seed = 1
)
head(d)
#> Y1 Y2 Y3 Y4 Y5 Y6 X Zp
#> 1 1 1 1 0 0 0 2 1
#> 2 0 0 0 0 0 0 3 4
#> 3 1 0 1 0 0 0 2 1
#> 4 1 1 0 1 1 1 1 2
#> 5 0 0 0 0 0 0 3 5
#> 6 1 1 1 1 1 1 1 3
# Low separation: P(Y_h = 1 | class) = 0.7 / 0.3
# Zp and X are identical to 'd' because seed = 1
d.low <- generate_data(
n = 500,
separation = "low",
scenario = "covariate",
seed = 1
)
head(d.low)
#> Y1 Y2 Y3 Y4 Y5 Y6 X Zp
#> 1 1 1 1 0 0 0 2 1
#> 2 1 0 0 1 0 1 3 4
#> 3 0 0 1 0 0 0 2 1
#> 4 1 1 0 0 1 1 1 2
#> 5 0 0 1 1 1 0 3 5
#> 6 1 1 1 1 1 1 1 3
Step 1: Measurement model
three_step() with no Zp.names or
Zo.name fits the measurement model only, returning a
tseLCA_measurement object. Internally this calls
multilevLCA::multiLCA() with random restarts when entropy
\(R^2\) is low.
d.measurement <- three_step(
data = d,
Y.names = paste0("Y", 1:6),
n_classes = 3,
measurement.tol = 1e-8
)
summary(d.measurement)
#> -- tseLCA Measurement Model --------------------------------
#> Latent classes : 3
#> Log-likelihood : -1455.5052
#> AIC : 2951.0104
#> BIC : 3035.3025
#> Entropy R² : 0.8780
#>
#> Class prevalences:
#>
#> P(C1) 0.3570
#> P(C2) 0.3308
#> P(C3) 0.3122
#> attr(,"names")
#> [1] "C1" "C2" "C3"
#>
#> Item-response probabilities (P(Y=1|class)):
#> C1 C2 C3
#> P(Y1|C) 0.9237 0.8621 0.1187
#> P(Y2|C) 0.9083 0.9219 0.1178
#> P(Y3|C) 0.9148 0.9571 0.0731
#> P(Y4|C) 0.8843 0.1481 0.0875
#> P(Y5|C) 0.8817 0.1340 0.1118
#> P(Y6|C) 0.9174 0.0889 0.1252
With low separation the measurement model can struggle to find the
global maximum. Use iter.measurement to trigger the number
of random restarts whenever entropy \(R^2\) falls below
R2.threshold.
d.low.measurement <- three_step(
data = d.low,
Y.names = paste0("Y", 1:6),
n_classes = 3,
iter.measurement = 10,
R2.threshold = 0.9
)
summary(d.low.measurement)
#> -- tseLCA Measurement Model --------------------------------
#> Latent classes : 3
#> Log-likelihood : -2019.2458
#> AIC : 4078.4916
#> BIC : 4162.7837
#> Entropy R² : 0.3327
#>
#> Class prevalences:
#>
#> P(C1) 0.2753
#> P(C2) 0.4551
#> P(C3) 0.2696
#> attr(,"names")
#> [1] "C1" "C2" "C3"
#>
#> Item-response probabilities (P(Y=1|class)):
#> C1 C2 C3
#> P(Y1|C) 0.7328 0.6418 0.3345
#> P(Y2|C) 0.5723 0.7549 0.3223
#> P(Y3|C) 0.6937 0.7101 0.3036
#> P(Y4|C) 0.6846 0.4588 0.2105
#> P(Y5|C) 0.6947 0.4058 0.3787
#> P(Y6|C) 0.8651 0.3551 0.2456
The plot() S3 method delegates to
multilevLCA’s item-profile plot.
Two-step estimates
fitZ_from_fit0() fixes the measurement parameters at
their Step-1 values and estimates multinomial logit coefficients \(\gamma\) via EM. These two-step estimates
serve as starting values for Step 3 and are generally close to the final
three-step estimates.
d.fitZ <- fitZ_from_fit0(
fit0 = d.measurement$measurement_model$fit0,
data = d,
Y.names = paste0("Y", 1:6),
Zp.names = "Zp"
)
# True slopes: -1 (C2) and +1 (C3) relative to C1
d.fitZ$mGamma
#> C2 C3
#> Intercept 2.1934130 -3.4524271
#> Zp -0.9411383 0.8971774
Starting values from the high-separation fit can be passed to the
low-separation fit to help it converge.
d.low.fitZ <- fitZ_from_fit0(
fit0 = d.low.measurement$measurement_model$fit0,
data = d.low,
Y.names = paste0("Y", 1:6),
Zp.names = "Zp",
starting_val = d.fitZ$mGamma
)
d.low.fitZ$mGamma
#> C2 C3
#> Intercept 3.0446359 -3.6948016
#> Zp -0.9832601 0.9487392
Three-Step estimation
ML estimator (default)
A single three_step() call handles all three steps. By
default it uses the ML correction of Vermunt (2010) and modal class
assignment.
d.three_step <- three_step(
data = d,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp"
)
summary(d.three_step)
#> -- tseLCA Three-step Covariate Model -----------------------
#> Latent classes : 3
#> Estimator : ML
#> Log-likelihood : -1339.0650
#> AIC : 2758.1299
#> BIC : 2926.7143
#> Entropy R² : 0.8693 (covariate-adjusted)
#>
#> Two-step (starting) estimates:
#> C2 C3
#> Intercept 2.1934 -3.4524
#> Zp -0.9411 0.8972
#>
#> Three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C2 2.0411 0.3237 6.3050 < 0.001 ***
#> Zp:C2 -0.8821 0.1406 -6.2730 < 0.001 ***
#> Intercept:C3 -3.4836 0.5913 -5.8913 < 0.001 ***
#> Zp:C3 0.8985 0.1435 6.2606 < 0.001 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The standard coef() and vcov() S3 methods
work on any tseLCA object.
coef(d.three_step)
#> C2 C3
#> Intercept 2.0410764 -3.4835616
#> Zp -0.8820801 0.8984978
vcov(d.three_step)
#> Intercept:C2 Zp:C2 Intercept:C3 Zp:C3
#> Intercept:C2 0.104798063 -0.041343485 -0.01048872 0.001510799
#> Zp:C2 -0.041343485 0.019772531 0.01121397 -0.001932282
#> Intercept:C3 -0.010488717 0.011213968 0.34964328 -0.082842095
#> Zp:C3 0.001510799 -0.001932282 -0.08284210 0.020597096
Proportional assignment
With modal assignment (use.modal.assignment = TRUE, the
default), the Jacobian in the measurement-uncertainty correction is not
mathematically defined. Setting
use.modal.assignment = FALSE uses soft posterior weights
throughout, giving an analytic Jacobian and is recommended when
separation is moderate or low. When
use.modal.assignment = TRUE, the Jacobian \(\frac{\partial\theta_2}{\partial\theta_1}\)
computed using the full posterior weights (e.g., behaving as if
use.modal.assignment = FALSE) to maintain well-defined
derivatives, though three-step estimates would still be computed with
modal assignment as specified. The different is negligible when
separation is high.
d.three_step.prop <- three_step(
data = d,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp",
use.modal.assignment = FALSE
)
summary(d.three_step.prop)
#> -- tseLCA Three-step Covariate Model -----------------------
#> Latent classes : 3
#> Estimator : ML
#> Log-likelihood : -1339.0617
#> AIC : 2758.1234
#> BIC : 2926.7078
#> Entropy R² : 0.8680 (covariate-adjusted)
#>
#> Two-step (starting) estimates:
#> C2 C3
#> Intercept 2.1934 -3.4524
#> Zp -0.9411 0.8972
#>
#> Three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C2 2.1096 0.3287 6.4183 < 0.001 ***
#> Zp:C2 -0.9121 0.1460 -6.2456 < 0.001 ***
#> Intercept:C3 -3.6508 0.6430 -5.6777 < 0.001 ***
#> Zp:C3 0.9367 0.1547 6.0563 < 0.001 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Simple (robust) standard errors
Setting use.simple.cov = TRUE skips the
measurement-uncertainty correction and returns the robust sandwich SEs
from Step 3 only. When separation is high the correction is negligible,
so this is a useful computational shortcut for large samples.
d.three_step.simple <- three_step(
data = d,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp",
use.simple.cov = TRUE
)
summary(d.three_step.simple)
#> -- tseLCA Three-step Covariate Model -----------------------
#> Latent classes : 3
#> Estimator : ML
#> Log-likelihood : -1339.0650
#> AIC : 2758.1299
#> BIC : 2926.7143
#> Entropy R² : 0.8693 (covariate-adjusted)
#>
#> Two-step (starting) estimates:
#> C2 C3
#> Intercept 2.1934 -3.4524
#> Zp -0.9411 0.8972
#>
#> Three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C2 2.0411 0.3214 6.3504 < 0.001 ***
#> Zp:C2 -0.8821 0.1394 -6.3257 < 0.001 ***
#> Intercept:C3 -3.4836 0.5876 -5.9281 < 0.001 ***
#> Zp:C3 0.8985 0.1426 6.3008 < 0.001 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BCH estimator
The BCH correction of Bolck, Croon & Hagenaars (2004) is
available via use.bch = TRUE. It works well with high
separation but can produce an ill-conditioned Hessian when separation is
low (resulting in a covariance matrix that is not positive
semi-definite), in which case the ML estimator is preferred.
d.three_step.bch <- three_step(
data = d,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp",
use.bch = TRUE
)
summary(d.three_step.bch)
#> -- tseLCA Three-step Covariate Model -----------------------
#> Latent classes : 3
#> Estimator : BCH
#> Log-likelihood : -1339.2863
#> AIC : 2758.5726
#> BIC : 2927.1569
#> Entropy R² : 0.8700 (covariate-adjusted)
#>
#> Two-step (starting) estimates:
#> C2 C3
#> Intercept 2.1934 -3.4524
#> Zp -0.9411 0.8972
#>
#> Three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C2 1.9554 0.3111 6.2844 < 0.001 ***
#> Zp:C2 -0.8424 0.1304 -6.4613 < 0.001 ***
#> Intercept:C3 -3.4634 0.5697 -6.0790 < 0.001 ***
#> Zp:C3 0.8923 0.1385 6.4412 < 0.001 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BCH with low-separation data can fail to produce a positive
semi-definite Hessian. The ML estimator with proportional assignment is
more reliable in this setting.
# Not run in vignette build (slow and and produces warnings)
bch.fail <- three_step(
data = d.low,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp",
use.bch = TRUE,
maxIter.measurement = 2000,
iter.measurement = 10
)
# Preferred approach for low separation
d.low.three_step.prop <- three_step(
data = d.low,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp",
use.modal.assignment = FALSE
)
summary(d.low.three_step.prop)
#> -- tseLCA Three-step Covariate Model -----------------------
#> Latent classes : 3
#> Estimator : ML
#> Log-likelihood : -1979.3372
#> AIC : 4038.6744
#> BIC : 4207.2588
#> Entropy R² : 0.3518 (covariate-adjusted)
#>
#> Two-step (starting) estimates:
#> C2 C3
#> Intercept 3.0267 -3.6919
#> Zp -0.9767 0.9482
#>
#> Three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C2 3.2034 2.2929 1.3971 0.1624
#> Zp:C2 -1.0761 1.9088 -0.5638 0.5729
#> Intercept:C3 -3.8431 2.9955 -1.2830 0.1995
#> Zp:C3 0.9554 0.6034 1.5832 0.1134
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Choosing the reference class
By default, class 1 ("C1") is the reference category for
the multinomial logit parameterization. The rebase argument
changes this. Estimates are reparameterized consistently:
log-likelihoods are invariant, and the coefficients satisfy the
transitivity relation \(\log(\pi_t / \pi_j) =
\log(\pi_t / \pi_1) - \log(\pi_j / \pi_1)\).
# Default: C1 as reference
summary(d.three_step.simple)
#> -- tseLCA Three-step Covariate Model -----------------------
#> Latent classes : 3
#> Estimator : ML
#> Log-likelihood : -1339.0650
#> AIC : 2758.1299
#> BIC : 2926.7143
#> Entropy R² : 0.8693 (covariate-adjusted)
#>
#> Two-step (starting) estimates:
#> C2 C3
#> Intercept 2.1934 -3.4524
#> Zp -0.9411 0.8972
#>
#> Three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C2 2.0411 0.3214 6.3504 < 0.001 ***
#> Zp:C2 -0.8821 0.1394 -6.3257 < 0.001 ***
#> Intercept:C3 -3.4836 0.5876 -5.9281 < 0.001 ***
#> Zp:C3 0.8985 0.1426 6.3008 < 0.001 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
d.three_step.simpleC2 <- three_step(
data = d,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp",
use.simple.cov = TRUE,
rebase = "C2"
)
summary(d.three_step.simpleC2)
#> -- tseLCA Three-step Covariate Model -----------------------
#> Latent classes : 3
#> Estimator : ML
#> Log-likelihood : -1339.0650
#> AIC : 2758.1299
#> BIC : 2926.7143
#> Entropy R² : 0.8693 (covariate-adjusted)
#>
#> Two-step (starting) estimates:
#> C1 C3
#> Intercept -2.1941 -5.6433
#> Zp 0.9413 1.8377
#>
#> Three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C1 -2.0411 0.3214 -6.3504 < 0.001 ***
#> Zp:C1 0.8821 0.1394 6.3257 < 0.001 ***
#> Intercept:C3 -5.5246 0.6823 -8.0976 < 0.001 ***
#> Zp:C3 1.7806 0.2078 8.5696 < 0.001 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
d.three_step.simpleC3 <- three_step(
data = d,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp",
use.simple.cov = TRUE,
rebase = "C3"
)
summary(d.three_step.simpleC3)
#> -- tseLCA Three-step Covariate Model -----------------------
#> Latent classes : 3
#> Estimator : ML
#> Log-likelihood : -1339.0650
#> AIC : 2758.1299
#> BIC : 2926.7143
#> Entropy R² : 0.8693 (covariate-adjusted)
#>
#> Two-step (starting) estimates:
#> C1 C2
#> Intercept 3.4492 5.6433
#> Zp -0.8964 -1.8377
#>
#> Three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C1 3.4836 0.5876 5.9281 < 0.001 ***
#> Zp:C1 -0.8985 0.1426 -6.3008 < 0.001 ***
#> Intercept:C2 5.5246 0.6823 8.0976 < 0.001 ***
#> Zp:C2 -1.7806 0.2078 -8.5696 < 0.001 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Passing a pre-fitted measurement model
The step1 argument accepts any previously fitted
tseLCA object or the raw output of
lca_step1(). This is useful when you want to:
- Reuse an expensive measurement model across multiple structural
models.
- Estimate the measurement model on a large reference sample and apply
it to a smaller analysis sample.
- Inject custom two-step starting values computed via
fitZ_from_fit0().
# Reuse the measurement model estimated above
d.three_step.prop2 <- three_step(
data = d,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp",
use.modal.assignment = FALSE,
step1 = d.measurement$measurement_model
)
summary(d.three_step.prop2)
#> -- tseLCA Three-step Covariate Model -----------------------
#> Latent classes : 3
#> Estimator : ML
#> Log-likelihood : -1339.0617
#> AIC : 2758.1234
#> BIC : 2926.7078
#> Entropy R² : 0.8680 (covariate-adjusted)
#>
#> Two-step (starting) estimates:
#> C2 C3
#> Intercept 2.1934 -3.4524
#> Zp -0.9411 0.8972
#>
#> Three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C2 2.1096 0.3287 6.4183 < 0.001 ***
#> Zp:C2 -0.9121 0.1460 -6.2456 < 0.001 ***
#> Intercept:C3 -3.6508 0.6430 -5.6777 < 0.001 ***
#> Zp:C3 0.9367 0.1547 6.0563 < 0.001 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Measurement model from a larger low-separation sample
d.low2000 <- generate_data(
n = 2000,
separation = "low",
scenario = "covariate",
seed = 2
)
d.low.measurement2000 <- three_step(
data = d.low2000,
Y.names = paste0("Y", 1:6),
n_classes = 3
)
# Apply to the smaller sample; get.twostep.vcov returns multilevLCA's
# bias-corrected vcov for the two-step estimates
d.low.three_step.prop2 <- three_step(
data = d.low,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp",
use.modal.assignment = FALSE,
step1 = d.low.measurement2000$measurement_model,
get.twostep.vcov = TRUE
)
summary(d.low.three_step.prop2)
#> -- tseLCA Three-step Covariate Model -----------------------
#> Latent classes : 3
#> Estimator : ML
#> Log-likelihood : -1983.8234
#> AIC : 4047.6468
#> BIC : 4216.2311
#> Entropy R² : 0.3770 (covariate-adjusted)
#>
#> Two-step (starting) estimates:
#> C2 C3
#> Intercept 2.5851 -4.2908
#> Zp -1.3556 1.0806
#>
#> Three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C2 2.7051 1.1848 2.2831 0.0224 *
#> Zp:C2 -1.3769 0.9785 -1.4072 0.1594
#> Intercept:C3 -3.9488 2.0927 -1.8870 0.0592 .
#> Zp:C3 1.0303 0.4742 2.1730 0.0298 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
You can also compute two-step starting values separately and inject
them before calling three_step().
d.low.fitZ2 <- fitZ_from_fit0(
fit0 = d.low.measurement2000$measurement_model$fit0,
data = d.low,
Y.names = paste0("Y", 1:6),
Zp.names = "Zp"
)
d.low.measurement2000$measurement_model$fitZ <- d.low.fitZ2
d.low.three_step.prop3 <- three_step(
data = d.low,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp",
use.modal.assignment = FALSE,
step1 = d.low.measurement2000$measurement_model
)
summary(d.low.three_step.prop3)
#> -- tseLCA Three-step Covariate Model -----------------------
#> Latent classes : 3
#> Estimator : ML
#> Log-likelihood : -1983.8234
#> AIC : 4047.6468
#> BIC : 4216.2311
#> Entropy R² : 0.3770 (covariate-adjusted)
#>
#> Two-step (starting) estimates:
#> C2 C3
#> Intercept 2.5851 -4.2908
#> Zp -1.3556 1.0806
#>
#> Three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C2 2.7051 1.1848 2.2831 0.0224 *
#> Zp:C2 -1.3769 0.9785 -1.4072 0.1594
#> Intercept:C3 -3.9488 2.0927 -1.8870 0.0592 .
#> Zp:C3 1.0303 0.4742 2.1730 0.0298 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Missing data
tseLCA uses a two-pass row-filtering strategy that
matches multilevLCA’s approach for the measurement model
while allowing more observations into Steps 1 and 2 than Step 3.
set.seed(42)
d.new <- generate_data(500, separation = "high", seed = 3)
sparsity <- 0.1
missing <- 1 -
matrix(
rbinom(prod(dim(d.new)), size = 1, prob = sparsity),
nrow = nrow(d.new),
ncol = ncol(d.new)
)
missing[missing == 0] <- NA_real_
d.sparse <- d.new * missing
head(d.sparse)
#> Y1 Y2 Y3 Y4 Y5 Y6 X Zp
#> 1 0 0 NA 0 1 0 3 5
#> 2 1 1 NA 0 0 0 2 2
#> 3 1 1 1 1 1 0 1 4
#> 4 0 0 0 0 0 0 3 4
#> 5 1 1 1 0 0 0 NA 2
#> 6 1 1 NA 0 0 0 2 3
With incomplete = FALSE (the default), any row with a
missing indicator is dropped before the measurement model is
estimated.
d.sparse.measurement <- three_step(
data = d.sparse,
Y.names = paste0("Y", 1:6),
n_classes = 3,
incomplete = FALSE,
verbose = TRUE
)
#> 242 row(s) dropped from measurement/classification steps (missing Y).
# Rows dropped = number of rows with at least one missing Y
sum(apply(d.sparse[, paste0("Y", 1:6)], 1, \(x) any(is.na(x))))
#> [1] 242
summary(d.sparse.measurement)
#> -- tseLCA Measurement Model --------------------------------
#> Latent classes : 3
#> Log-likelihood : -742.0656
#> AIC : 1524.1311
#> BIC : 1595.1903
#> Entropy R² : 0.9027
#>
#> Class prevalences:
#>
#> P(C1) 0.2995
#> P(C2) 0.3967
#> P(C3) 0.3037
#> attr(,"names")
#> [1] "C1" "C2" "C3"
#>
#> Item-response probabilities (P(Y=1|class)):
#> C1 C2 C3
#> P(Y1|C) 0.8241 0.8869 0.0834
#> P(Y2|C) 0.8600 0.9104 0.0811
#> P(Y3|C) 0.9110 0.9344 0.0633
#> P(Y4|C) 0.9035 0.0486 0.1302
#> P(Y5|C) 0.9349 0.1517 0.0795
#> P(Y6|C) 0.9331 0.1555 0.0762
With incomplete = TRUE, only fully-missing rows are
dropped; partially observed rows contribute to the measurement model via
FIML.
d.sparse.measurement2 <- three_step(
data = d.sparse,
Y.names = paste0("Y", 1:6),
n_classes = 3,
incomplete = TRUE,
verbose = TRUE
)
summary(d.sparse.measurement2)
#> -- tseLCA Measurement Model --------------------------------
#> Latent classes : 3
#> Log-likelihood : -1342.8026
#> AIC : 2725.6052
#> BIC : 2809.8974
#> Entropy R² : 0.8425
#>
#> Class prevalences:
#>
#> P(C1) 0.3049
#> P(C2) 0.3652
#> P(C3) 0.3299
#> attr(,"names")
#> [1] "C1" "C2" "C3"
#>
#> Item-response probabilities (P(Y=1|class)):
#> C1 C2 C3
#> P(Y1|C) 0.8797 0.8916 0.0925
#> P(Y2|C) 0.8888 0.8858 0.0677
#> P(Y3|C) 0.9337 0.8859 0.1453
#> P(Y4|C) 0.9079 0.0819 0.1339
#> P(Y5|C) 0.9536 0.1359 0.1056
#> P(Y6|C) 0.9583 0.1590 0.1176
Regardless of incomplete, Step 3 drops any row with a
missing covariate. The rows used in Step 3 are a subset of those used in
Steps 1 and 2.
d.sparse.three_step <- three_step(
data = d.sparse,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp",
incomplete = TRUE,
verbose = TRUE
)
#> 43 row(s) excluded from covariate step (missing Z).
#> fitZ EM converged in 9 iterations.
#> 43 row(s) excluded from covariate step (missing Z).
#> EM converged in 8 iterations.
# Additional rows dropped from Step 3 due to missing Zp
sum(is.na(d.sparse$Zp))
#> [1] 43
summary(d.sparse.three_step)
#> -- tseLCA Three-step Covariate Model -----------------------
#> Latent classes : 3
#> Estimator : ML
#> Log-likelihood : -1088.3344
#> AIC : 2256.6688
#> BIC : 2421.6562
#> Entropy R² : 0.8672 (covariate-adjusted)
#>
#> Two-step (starting) estimates:
#> C2 C3
#> Intercept 2.5023 -4.6555
#> Zp -1.0494 1.2929
#>
#> Three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C2 2.4803 0.3807 6.5149 < 0.001 ***
#> Zp:C2 -1.0136 0.1667 -6.0789 < 0.001 ***
#> Intercept:C3 -4.7754 0.7267 -6.5716 < 0.001 ***
#> Zp:C3 1.3164 0.1796 7.3290 < 0.001 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A FIML measurement model can be passed in and then reused for the
covariate step on the same sparse data.
d.sparse.three_step2 <- three_step(
data = d.sparse,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp",
incomplete = TRUE,
step1 = d.sparse.measurement2$measurement_model,
verbose = TRUE
)
#> 43 row(s) excluded from covariate step (missing Z).
#> fitZ EM converged in 9 iterations.
#> 43 row(s) excluded from covariate step (missing Z).
#> EM converged in 8 iterations.
summary(d.sparse.three_step2)
#> -- tseLCA Three-step Covariate Model -----------------------
#> Latent classes : 3
#> Estimator : ML
#> Log-likelihood : -1088.3344
#> AIC : 2256.6688
#> BIC : 2421.6562
#> Entropy R² : 0.8672 (covariate-adjusted)
#>
#> Two-step (starting) estimates:
#> C2 C3
#> Intercept 2.5023 -4.6555
#> Zp -1.0494 1.2929
#>
#> Three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C2 2.4803 0.3807 6.5149 < 0.001 ***
#> Zp:C2 -1.0136 0.1667 -6.0789 < 0.001 ***
#> Intercept:C3 -4.7754 0.7267 -6.5716 < 0.001 ***
#> Zp:C3 1.3164 0.1796 7.3290 < 0.001 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Polytomous items
tseLCA supports polytomous indicators, following
multilevLCA’s convention that item categories are coded as
consecutive integers starting at 0.
Here we reproduce the example from the poLCA
package.
data(election, package = "poLCA")
elec <- election
elec.items <- colnames(election)[1:12]
# Recode to 0-based integers as required by multilevLCA
elec[, elec.items] <- lapply(elec[, elec.items], \(x) as.integer(x) - 1L)
elec.measurement <- three_step(
data = elec,
Y.names = elec.items,
n_classes = 3,
#The poLCA example drops any row with a missing cell
incomplete = FALSE
)
elec.three_step <- three_step(
data = elec,
Y.names = elec.items,
n_classes = 3,
Zp.names = c("PARTY"),
step1 = elec.measurement$measurement_model,
incomplete = FALSE,
#With the neutral group as the base-category
rebase = "C3"
)
#> Warning: lca_indiv_varmat: Infomat is singular even after removing boundary
#> parameters; returning NA matrix. Check for near-empty classes.
summary(elec.three_step)
#> -- tseLCA Three-step Covariate Model -----------------------
#> Latent classes : 3
#> Estimator : ML
#> Log-likelihood : -16278.0242
#> AIC : 32852.0485
#> BIC : 33617.2262
#> Entropy R² : 0.7956 (covariate-adjusted)
#>
#> Two-step (starting) estimates:
#> C1 C2
#> Intercept -2.5781 1.8687
#> PARTY 0.4289 -0.6983
#>
#> Three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C1 -2.4701 NA NA NA
#> PARTY:C1 0.4077 NA NA NA
#> Intercept:C2 1.7324 NA NA NA
#> PARTY:C2 -0.6727 NA NA NA
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
party.x <- seq(from = 1, to = 7, length.out = 101)
pidmat <- cbind(1, party.x)
exb <- exp(pidmat %*% coef(elec.three_step))
matplot(
party.x,
(cbind(1, exb)) / (1 + rowSums(exb)),
ylim = c(0, 1),
type = "l",
lwd = 3,
col = 1,
xlab = "Party ID: strong Democratic (1) to strong Republican (7)",
ylab = "Probability of latent class membership",
main = "Party ID as a predictor of candidate affinity class",
)
text(3.9, 0.60, "Other")
text(6.2, 0.6, "Bush affinity")
text(2.0, 0.65, "Gore affinity")
Distal outcomes
For distal outcomes (\(Z_o \leftarrow X
\rightarrow Y\)), supply Zo.name and a
family argument. The available families are
"gaussian" (default), "poisson", and
"binomial". Both ML and BCH estimators are available.
d.distal <- generate_data(
n = 500,
separation = "high",
scenario = "distal",
seed = 4
)
# True class means: mu = (0, 1, -1) for C1, C2, C3
d.distal.measurement <- three_step(
data = d.distal,
Y.names = paste0("Y", 1:6),
n_classes = 3
)
# ML estimator
d.distal.three_step.ml <- three_step(
data = d.distal,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zo.name = "Zo",
step1 = d.distal.measurement$measurement_model,
use.modal.assignment = FALSE,
family = "gaussian"
)
# BCH estimator: closed-form M-step for distal outcomes
d.distal.three_step.bch <- three_step(
data = d.distal,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zo.name = "Zo",
step1 = d.distal.measurement$measurement_model,
use.modal.assignment = FALSE,
use.bch = TRUE,
family = "gaussian"
)
summary(d.distal.three_step.ml)
#> -- tseLCA Three-step Distal Outcome Model -------------------
#> Latent classes : 3
#> Estimator : ML
#> Family : gaussian
#> Log-likelihood : -2169.0110
#> AIC : 4384.0220
#> BIC : 4480.9580
#>
#> Distal outcome estimates by class:
#> Estimate Std.Error z.value p.value
#> mu_C1 (mean) -1.0821 0.0817 -13.2495 < 0.001 ***
#> mu_C2 (mean) 1.0172 0.0819 12.4151 < 0.001 ***
#> mu_C3 (mean) 0.0254 0.0878 0.2887 0.7728
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(d.distal.three_step.bch)
#> -- tseLCA Three-step Distal Outcome Model -------------------
#> Latent classes : 3
#> Estimator : BCH
#> Family : gaussian
#> Log-likelihood : -2168.8685
#> AIC : 4383.7370
#> BIC : 4480.6730
#>
#> Distal outcome estimates by class:
#> Estimate Std.Error z.value p.value
#> mu_C1 (mean) -1.0941 0.0867 -12.6134 < 0.001 ***
#> mu_C2 (mean) 0.9751 0.0823 11.8534 < 0.001 ***
#> mu_C3 (mean) 0.0578 0.0859 0.6733 0.5008
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Three-step estimation with both covariates (Zp) and distal outcomes
(Zo)
Consistent with how most research in the social sciences construct
the relationships between \(Z_p\) and
\(X\), and \(X\) and \(Z_o\), the relationship between \(Z_p\) and \(X\) is estimated first,
followed by estimation between \(X\)
and \(Z_o\), adjusting for the
covariate-adjusted posteriors in the estimation procedures for the
distal outcome model in step 3.
d.covariate <- generate_data(
n = 500,
separation = "high",
scenario = "covariate",
seed = 4
)
d.covariate$Zo <- draw_Zo(d.covariate$X, bk2018_params$distal_params)
head(d.covariate)
#> Y1 Y2 Y3 Y4 Y5 Y6 X Zp Zo
#> 1 1 1 1 1 0 0 2 3 -0.1624650
#> 2 1 1 1 1 1 1 1 3 -1.1591833
#> 3 1 1 1 1 1 1 1 3 -1.2055132
#> 4 0 0 0 0 0 0 3 4 1.8752276
#> 5 1 1 1 1 1 1 1 3 -2.5582369
#> 6 1 1 1 1 1 1 1 5 -0.4723262
d.covariate.three_step <- three_step(
data = d.covariate,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp",
Zo.name = "Zo",
use.modal.assignment = FALSE
)
summary(d.covariate.three_step)
#> -- tseLCA Three-step Model: Covariate + Distal Outcome -----
#> Latent classes : 3
#> Estimator : ML
#> Family : gaussian
#> Log-likelihood : -1315.6596
#> AIC : 2711.3193
#> BIC : 2879.9036
#>
#> Covariate -- two-step (starting) estimates:
#> C2 C3
#> Intercept 2.4973 -4.2196
#> Zp -1.0177 1.1159
#>
#> Covariate -- three-step estimates:
#> Estimate Std.Error z.value p.value
#> Intercept:C2 2.6602 0.3998 6.6538 < 0.001 ***
#> Zp:C2 -1.0790 0.1632 -6.6134 < 0.001 ***
#> Intercept:C3 -4.6917 0.7070 -6.6365 < 0.001 ***
#> Zp:C3 1.2278 0.1735 7.0773 < 0.001 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Distal outcome -- three-step estimates:
#> Estimate Std.Error z.value p.value
#> mu_C1 (mean) -0.9851 0.0828 -11.8995 < 0.001 ***
#> mu_C2 (mean) 0.9298 0.0851 10.9217 < 0.001 ***
#> mu_C3 (mean) 0.1188 0.0722 1.6458 0.0998 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note that with covariates in a model with high separation, the
standard errors above should, on average, by systematically smaller for
distal outcome estimation than if there were no covariates in the model
(see below).
three_step(
data = d.covariate,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zo.name = "Zo",
use.modal.assignment = FALSE
) |>
vcov() |>
diag() |>
sqrt()
#> mu_C1 mu_C2 mu_C3
#> 0.08302639 0.08922331 0.07595896
References
Bakk, Z., Tekle, F. B., & Vermunt, J. K. (2013). Estimating the
association between latent class membership and external variables using
bias-adjusted three-step approaches. Sociological Methodology,
43(1), 272–311. https://doi.org/10.1177/0081175012470644
Bakk, Z., & Kuha, J. (2018). Two-step estimation of models
between latent classes and external variables. Psychometrika,
83(4), 871–892. https://doi.org/10.1007/s11336-017-9592-7
Bolck, A., Croon, M., & Hagenaars, J. (2004). Estimating latent
structure models with categorical variables: One-step versus three-step
estimators. Political Analysis, 12(1), 3–27. https://doi.org/10.1093/pan/mph001
Lyrvall, J., Di Mari, R., Bakk, Z., Oser, J., & Kuha, J. (2025).
Multilevel latent class analysis: State-of-the-art methodologies and
their implementation in the R package multilevLCA. Multivariate
Behavioral Research, 60(4), 731–747. https://doi.org/10.1080/00273171.2025.2473935
Vermunt, J. K. (2010). Latent class modeling with covariates: Two
improved three-step approaches. Political Analysis, 18(4),
450–469. https://doi.org/10.1093/pan/mpq025
sessionInfo()
#> R version 4.5.1 (2025-06-13 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26200)
#>
#> Matrix products: default
#> LAPACK version 3.12.1
#>
#> locale:
#> [1] LC_COLLATE=C
#> [2] LC_CTYPE=English_United States.utf8
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=English_United States.utf8
#>
#> time zone: America/Denver
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] tseLCA_1.0.2
#>
#> loaded via a namespace (and not attached):
#> [1] sass_0.4.10 generics_0.1.4 tidyr_1.3.1 pracma_2.4.6
#> [5] hms_1.1.4 digest_0.6.39 magrittr_2.0.3 evaluate_1.0.5
#> [9] RColorBrewer_1.1-3 iterators_1.0.14 fastmap_1.2.0 foreach_1.5.2
#> [13] jsonlite_2.0.0 combinat_0.0-8 promises_1.5.0 purrr_1.0.4
#> [17] codetools_0.2-20 jquerylib_0.1.4 cli_3.6.6 shiny_1.14.0
#> [21] labelled_2.16.0 rlang_1.2.0 cachem_1.1.0 yaml_2.3.12
#> [25] otel_0.2.0 klaR_1.7-4 parallel_4.5.1 tools_4.5.1
#> [29] dplyr_1.1.4 httpuv_1.6.17 forcats_1.0.1 vctrs_0.6.5
#> [33] R6_2.6.1 mime_0.13 lifecycle_1.0.5 multilevLCA_2.1.4
#> [37] tictoc_1.2.1 MASS_7.3-65 miniUI_0.1.2 cluster_2.1.8.1
#> [41] pkgconfig_2.0.3 pillar_1.11.1 bslib_0.11.0 later_1.4.8
#> [45] glue_1.8.0 Rcpp_1.1.1-1.1 haven_2.5.5 xfun_0.52
#> [49] tibble_3.2.1 tidyselect_1.2.1 highr_0.12 rstudioapi_0.19.0
#> [53] knitr_1.51 xtable_1.8-8 htmltools_0.5.9 rmarkdown_2.31
#> [57] clustMixType_0.5-1 compiler_4.5.1 questionr_0.8.2