three_step: Three-step LCA estimation with covariates and/or distal...
In tseLCA: Three-Step Estimation for Latent Class Analysis

three_step

R Documentation

Three-step LCA estimation with covariates and/or distal outcomes

Description

Fits a three-step latent class model through the following steps:

Measurement model: estimates latent class parameters (\pi, \phi) using multilevLCA (Lyrvall et al., 2025).
Classification-error matrix: computes posterior class probabilities and the T x T misclassification probability matrix P(W = s \mid X = t), with standard errors corrected for classification-error propagation (Bakk, Oberski & Vermunt, 2014).
Structural model: estimates covariate effects using two-step starting values (Bakk & Kuha, 2018) and/or distal outcome means following Bakk, Tekle & Vermunt (2013), with the ML correction (Vermunt, 2010) or BCH correction (Bolck, Croon & Hagenaars, 2004). See vignette("tseLCA", package = "tseLCA") for a worked example.

Usage

three_step(
  data,
  Y.names,
  n_classes,
  Zp.names = NULL,
  Zo.name = NULL,
  step1 = NULL,
  startval = NULL,
  n_init = NULL,
  use.two.step = TRUE,
  use.modal.assignment = TRUE,
  include.intercept = TRUE,
  use.simple.cov = FALSE,
  incomplete = FALSE,
  boundary.tol = 0.01,
  maxIter.measurement = 5000,
  measurement.tol = 1e-08,
  covariate.tol = 1e-06,
  iter.measurement = 10L,
  R2.threshold = 0.7,
  use.bch = FALSE,
  em.maxIter = 200L,
  get.twostep.vcov = FALSE,
  rebase = "C1",
  family = "gaussian",
  correct.spec = FALSE,
  verbose = FALSE
)

Arguments

`data`	A data.frame containing all columns referenced by `Y.names`, `Zp.names`, and `Zo.name`.
`Y.names`	Character vector of indicator column names. Indicators may be factors, logicals, character, or numeric codes (see `tse_lca()`).
`n_classes`	Integer. Number of latent classes.
`Zp.names`	Character vector of covariate column names, or `NULL` for a measurement-only fit. Default `NULL`.
`Zo.name`	Single character name of the distal outcome column, or `NULL`. Default `NULL`.
`step1`	Pre-fitted Step-1 object (output of `lca_step1()` or a prior `three_step()` call), or `NULL` to run Step 1 internally. Default `NULL`.
`startval`	Optional starting classification for the Step-1 measurement model, either an integer vector of length `nrow(data)` (a class assignment `1..n_classes` for every row) or a numeric matrix of conditional item-response probabilities `P(Y_h = k \mid X = t)` (one row per item-category pair in `Y.names` order, one column per class) from which a classification is derived internally. See `lca_step1_startval()` for the full description of both forms and typical sources (an external solver run with many random starts, or a published item-response table). multilevLCA's default initialization (k-means on principal components) is deterministic given the data and can converge to a local optimum of the Step-1 log-likelihood; supplying `startval` bypasses it (`kmea = FALSE` with the classification injected as multilevLCA's `startval`). Mutually exclusive with `step1` and `n_init`. Default `NULL`.
`n_init`	Optional positive integer. If supplied, fits the Step-1 measurement model `n_init` times from independent uniform-random classifications (`kmea = FALSE`, not multilevLCA's k-means-on-PCA path) and keeps the fit with the highest log-likelihood – the unconditional multi-start analog of `n_init` in StepMix or `nrep` in poLCA. This is distinct from `iter.measurement`, which reruns multilevLCA's own k-means initialization, and only when the entropy R`^2` of the default fit is below `R2.threshold`; `n_init` restarts always run. Mutually exclusive with `step1` and `startval`. Default `NULL`.
`use.two.step`	Logical. Initialize Step-3 from two-step estimates. Default `TRUE`.
`use.modal.assignment`	Logical. Use modal (hard) class assignments in Step 2 and 3. `FALSE` uses soft posterior weights. Default `TRUE`.
`include.intercept`	Logical. Prepend an intercept column to the covariate design matrix. Default `TRUE`.
`use.simple.cov`	Logical. Skip the Step-1 measurement-uncertainty correction and return only the robust sandwich variance. Faster but underestimates standard errors when class separation is low. Default `FALSE`.
`incomplete`	Logical. FIML for partially missing indicators. See the `Missing Data` section of `vignette("tseLCA", package = "tseLCA")`. Default `FALSE`.
`boundary.tol`	Scalar. Parameters within this tolerance of 0 or 1 are treated as fixed when computing the Step-1 variance matrix for numerical stability. Default `1e-2`.
`maxIter.measurement`	Integer. Maximum EM iterations for Step 1. Default `5000L`.
`measurement.tol`	Scalar. Convergence tolerance for the Step-1 EM algorithm. Default `1e-8`.
`covariate.tol`	Scalar. Convergence tolerance for the Step-3 Newton-Raphson or EM algorithm. Default `1e-6`.
`iter.measurement`	Integer. Number of random restarts triggered when the Step-1 entropy R`^2` falls below `R2.threshold`. Default `10L`.
`R2.threshold`	Scalar. Entropy R`^2` threshold below which Step-1 random restarts are triggered. Default `0.70`.
`use.bch`	Logical. Use the BCH estimator in Step 3 (default: the ML estimator). May error if BCH weights induce a non-positive semi-definite Hessian in the third step (common in cases of low separation). Default `FALSE`.
`em.maxIter`	Integer. Maximum EM iterations for the Step-3 covariate or distal outcome model. Default `200L`.
`get.twostep.vcov`	Logical. If `TRUE`, obtain multilevLCA's bias-corrected variance-covariance matrix for the two-step gamma estimates and store it in `$two_step_vcov`. If the `fitZ` object passed through `step1` already contains a `Varmat_cor` (from a prior `fitZ_from_multiLCA()` or plain `multiLCA` call), it is attached automatically even when `get.twostep.vcov = FALSE`. Default `FALSE`.
`rebase`	Character (e.g. `"C1"`, `"C2"`) or integer specifying which latent class to use as the reference category in the multinomial logit. The measurement model is permuted so this class becomes column 1 before any structural estimation. Default `"C1"`.
`family`	Character. Distal outcome family: one of `"gaussian"` (class means), `"poisson"` (log-rates), `"binomial"` (logits), or `"multinomial"` (a saturated model for a nominal categorical outcome with 2 or more categories – `Zo.name` may be a factor, character, or integer column; categories are taken from `sort(unique(data[[Zo.name]]))` with `factor()`). For `"multinomial"`, `coef()` returns a `T x C` matrix of class-conditional category probabilities `\hat\pi_{tc} = P(Zo = c \mid X = t)` (rows sum to 1), not a length-`T` vector, and `vcov()` returns its `(TC) x (TC)` sandwich covariance (necessarily singular, since each class's row sums to 1 – see `omnibus_test()` for a Wald test that accounts for this). Unlike `"binomial"`, whose `coef()`/`vcov()` are on the logit scale, `"multinomial"` reports `coef()`/`vcov()` directly on the probability scale, so `Std.Error` is directly interpretable without a delta-method back-transform – but a symmetric interval `Estimate +/- 1.96*Std.Error` can fall outside `[0, 1]` for a probability near a boundary, the same well-known limitation as a naive Wald interval for any sample proportion. The `z.value`/`p.value` columns `summary()`/`print()` show for this family test each probability against 0, which is rarely the question of interest; `omnibus_test()` is the intended, boundary-safe test of whether the outcome's distribution differs across classes. Combining `family = "multinomial"` with both `Zp.names` and `Zo.name` fully propagates both Step-1 measurement and Step-3 covariate uncertainty under `use.simple.cov = FALSE`, the same as the other families. Default `"gaussian"`.
`correct.spec`	Logical. Estimate the Step-3 information matrix by the outer product of the case-wise scores, not the observed-data Hessian. Valid only when the Step-3 model is correctly specified; the default observed-Hessian sandwich is robust to misspecification. Default `FALSE`.
`verbose`	Logical. Print convergence messages. Default `FALSE`.

Value

An S3 object of class tseLCA. The subclass depends on which models were estimated:

tseLCA_measurement

Returned when neither Zp.names nor Zo.name is supplied. Contains the following elements:

measurement_model: Step-1 output list from lca_step1().
llik: Final Step-1 log-likelihood.
AIC, BIC: Information criteria from the measurement model.
R2entr: Entropy R^2 of the measurement model.
n_classes: Number of latent classes.
posteriors: n x T matrix of soft posterior class probabilities.
classifications: Length-n integer vector of modal class assignments.

tseLCA_covariate

Returned when Zp.names is supplied and Zo.name is NULL. Contains all elements of tseLCA_measurement plus:

three_step: (Q+1) x (T-1) matrix of Step-3 gamma coefficients.
three_step_vcov: (Q+1)(T-1) x (Q+1)(T-1) variance-covariance matrix for three_step, with measurement-uncertainty correction unless use.simple.cov = TRUE.
two_step: (Q+1) x (T-1) matrix of two-step starting values, or NULL if use.two.step = FALSE.
two_step_vcov: multilevLCA bias-corrected vcov for the two-step estimates, or NULL.
estimator: Character: "ML" or "BCH".
entropy.R2: Covariate-adjusted entropy R^2.
llik: Profile log-likelihood \sum_i \log \sum_t P(X=t|Z_{p,i};\hat{\gamma}) P(Y_i|X=t;\hat{\phi}), with Step-1 parameters \hat{\phi} held fixed. By construction smaller than the equivalent one-step MLE likelihood.

tseLCA_distal

Returned when Zo.name is supplied and Zp.names is NULL. Contains:

three_step: Named length-T vector of Step-3 distal outcome parameters (means, log-rates, or logits depending on family) – or, for family = "multinomial", a T x C matrix of class-conditional category probabilities (rows sum to 1).
three_step_vcov: T x T variance-covariance matrix for three_step, named mu_C1 through mu_CT – or, for family = "multinomial", a (T*C) x (T*C) (necessarily rank-deficient) matrix named "C{t}:{category}".
three_step.llik: Step-3 distal log-likelihood \log P(Z_o|X=t) at converged estimates.
llik: Profile log-likelihood \sum_i \log \sum_t P(X=t|\hat{\pi}) P(Z_{o,i}|X=t;\hat{\mu}) P(Y_i|X=t;\hat{\phi}), with Step-1 parameters \hat{\pi}, \hat{\phi} held fixed. By construction smaller than the equivalent one-step MLE likelihood.
AIC: Akaike information criterion based on llik.
BIC: Bayesian information criterion based on llik, using the number of distal-complete observations.
family: Character. The distal outcome family used.
estimator: Character: "ML" or "BCH".
posteriors: n x T soft posterior matrix.
classifications: Length-n modal class assignment vector.

tseLCA_both

Returned when both Zp.names and Zo.name are supplied. Contains:

covariate: A tseLCA_covariate-structured sub-list (see above), including llik, AIC, BIC, entropy.R2.
distal: A tseLCA_distal-structured sub-list (see above), including llik, AIC, BIC, three_step.llik.
family, n_classes, estimator: Shared top-level fields.
posteriors, classifications: Shared n x T posterior matrix and length-n modal class vector.

Deprecated

Deprecated as of tseLCA 2.0.0. It keeps working (and gives the same estimates) but warns once per session; set options(tseLCA.warn.deprecated = FALSE) to silence the warning. Use tseLCA() or the step-wise functions:

`three_step()`	tseLCA 2.0
`Y.names`, `n_classes`	`tse_lca(cbind(...) ~ 1, nclass = )`
`Zp.names`	`tse_covariate(, ~ ...)` or `tseLCA(... ~ covariates)`
`Zo.name`, `family`	`tse_distal(, outcome ~ 1, family = )`, or `tseLCA()` with the outcome after the bar
`step1` (measurement model from another sample)	`tse_classify(, newdata = )`
`startval`	`tse_lca(start = )`
`use.modal.assignment`	`tse_classify(assignment = )`
`use.bch`	`method = "BCH"`
`use.simple.cov`	`se = "robust"`
`rebase`	`ref` argument, or `relevel()`
`incomplete`	`tse_lca(missing = "fiml")`
`n_init`, `maxIter.measurement`, `measurement.tol`, `iter.measurement`, `R2.threshold`, `em.maxIter`, `covariate.tol`, `boundary.tol`, `correct.spec`	`tse_control()`
`get.twostep.vcov`	`tse_twostep(se = TRUE)`

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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1177/0081175012470644")}

Bakk, Z., & Kuha, J. (2018). Two-step estimation of models between latent classes and external variables. Psychometrika, 83(4), 871–892. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/s11336-017-9592-7")}

Bakk, Z., Pohle, M. J., & Kuha, J. (2025). Bias-adjusted three-step estimation of structural models for latent classes. Multivariate Behavioral Research. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/00273171.2025.2473935")}

Examples

d <- generate_data(n = 200, separation = "high",
                   scenario = "covariate", seed = 1)

# Measurement model only
fit_m <- three_step(d, Y.names = paste0("Y", 1:6), n_classes = 3)
summary(fit_m)

# ML three-step with simple SEs (fast)
fit <- three_step(d, Y.names = paste0("Y", 1:6), n_classes = 3,
                  Zp.names = "Zp", use.simple.cov = TRUE)
summary(fit)
coef(fit)
vcov(fit)

# Full measurement-uncertainty correction (see vignette for interpretation)
fit_cor <- three_step(d, Y.names = paste0("Y", 1:6), n_classes = 3,
                      Zp.names = "Zp", use.simple.cov = FALSE,
                      use.modal.assignment = FALSE)
summary(fit_cor)

# BCH estimator
fit_bch <- three_step(d, Y.names = paste0("Y", 1:6), n_classes = 3,
                      Zp.names = "Zp", use.bch = TRUE,
                      use.simple.cov = TRUE)
summary(fit_bch)

# Change reference class
fit_c2 <- three_step(d, Y.names = paste0("Y", 1:6), n_classes = 3,
                     Zp.names = "Zp", use.simple.cov = TRUE,
                     rebase = "C2")
summary(fit_c2)

# Gaussian distal outcome
d2 <- generate_data(200, "high", "distal", seed = 2)
fit_dis <- three_step(d2, Y.names = paste0("Y", 1:6), n_classes = 3,
                      Zo.name = "Zo", family = "gaussian",
                      use.simple.cov = TRUE)
summary(fit_dis)

# Nominal categorical distal outcome (3+ categories): coef() returns a
# T x C matrix of class-conditional category probabilities; omnibus_test()
# gives a single Wald test of whether the category distribution differs
# across classes at all.
d2$Zcat <- factor(sample(c("low", "mid", "high"), nrow(d2), replace = TRUE))
fit_cat <- three_step(d2, Y.names = paste0("Y", 1:6), n_classes = 3,
                      Zo.name = "Zcat", family = "multinomial",
                      use.simple.cov = TRUE)
coef(fit_cat)
omnibus_test(fit_cat)

# Pass a pre-fitted measurement model to skip Step 1
fit_step1 <- three_step(d, Y.names = paste0("Y", 1:6), n_classes = 3)
fit2 <- three_step(d, Y.names = paste0("Y", 1:6), n_classes = 3,
                   Zp.names = "Zp", step1 = fit_step1,
                   use.simple.cov = TRUE)
summary(fit2)

# Supply an external starting classification for Step 1 (bypasses
# multilevLCA's k-means-on-PCA initialization; here we use the DGP's own
# true classes as a stand-in for e.g. a StepMix solution with many
# random starts)
fit_ext <- three_step(d, Y.names = paste0("Y", 1:6), n_classes = 3,
                      startval = d$X, use.simple.cov = TRUE)
summary(fit_ext)

# Many random-classification restarts for Step 1, keeping the best
# (analogous to n_init in StepMix or nrep in poLCA)
fit_ninit <- three_step(d, Y.names = paste0("Y", 1:6), n_classes = 3,
                        n_init = 20L, use.simple.cov = TRUE)
summary(fit_ninit)

# Plot item-response profiles from the measurement model
plot(fit)

tseLCA documentation built on Nov. 11, 2026, 9:07 a.m.