| Type: | Package |
| Title: | Three-Step Estimation for Latent Class Analysis |
| Version: | 1.0.2 |
| Description: | Implements BCH (Bolck-Croon-Hagenaars) <doi:10.1093/pan/mph001> and ML (Vermunt's maximum likelihood) <doi:10.1093/pan/mpq025> approaches for three-step estimation of latent class models with covariates and distal outcomes, following Bakk, Tekle & Vermunt (2013) <doi:10.1177/0081175012470644>, Bakk, Oberski & Vermunt (2014) https://www.jstor.org/stable/24573086, and Bakk & Kuha (2018) <doi:10.1007/s11336-017-9592-7>. Built on 'multilevLCA' (Lyrvall et al., 2025) <doi:10.1080/00273171.2025.2473935> for Step-1 measurement model estimation, this package extends it with support for Gaussian, Poisson, and binomial distal outcome families. Unlike 'poLCA', which relies on one-step estimation and cannot accommodate a measurement model from a different sample, this package uses a stepwise approach to prevent the structural model from influencing latent class formation. Implements correct sandwich variance estimation that propagates measurement uncertainty from the first-step through classification-error correction in the final step (Bakk, Oberski & Vermunt, 2014). Supports polytomous items and missing data in the measurement model with full information maximum likelihood. A data-generating process replicating the Bakk & Kuha (2018) simulation study is included. |
| License: | GPL (≥ 3) |
| Encoding: | UTF-8 |
| Language: | en-US |
| Depends: | R (≥ 4.1.0) |
| Imports: | multilevLCA, cli |
| Suggests: | poLCA, testthat (≥ 3.0.0), parallel, knitr, rmarkdown, spelling |
| Config/testthat/edition: | 3 |
| Config/roxygen2/version: | 8.0.0 |
| URL: | https://samleebyu.github.io/tseLCA/, https://github.com/SamLeeBYU/tseLCA |
| BugReports: | https://github.com/SamLeeBYU/tseLCA/issues |
| VignetteBuilder: | knitr, rmarkdown |
| NeedsCompilation: | no |
| Packaged: | 2026-07-05 03:07:06 UTC; samle |
| Author: | Sam Lee |
| Maintainer: | Sam Lee <samlee@arizona.edu> |
| Repository: | CRAN |
| Date/Publication: | 2026-07-11 08:40:02 UTC |
tseLCA: Three-Step Estimation for Latent Class Analysis
Description
tseLCA implements bias-adjusted three-step estimators for structural latent class models with covariates and distal outcomes. Building on the efficient measurement-model estimation in multilevLCA (Lyrvall et al., 2025), tseLCA adds modern three-step estimators, classification-error corrections, and analytic sandwich variance estimation that propagates measurement uncertainty from the latent class measurement stage through to the final structural parameter estimates.
In contrast to one-step approaches such as poLCA, where including covariates may alter the underlying latent class definitions, three-step methods fix the measurement model before estimating structural relationships and adjust for classification error at the final stage. tseLCA also allows the measurement and structural models to be estimated on different datasets, enabling researchers to calibrate a measurement model on a large reference sample and apply it to a separate analysis sample.
The three-step approach
-
Measurement model: estimate class-conditional item-response probabilities
\phiand class prevalences\piusing multilevLCA (Lyrvall et al., 2025 ). -
Classification-error matrix: assign posterior class probabilities and compute the T x T misclassification matrix
P(W = s \mid X = t), with standard errors corrected for classification-error propagation (Bakk, Oberski & Vermunt, 2014). -
Structural model: estimate covariate effects using two-step starting values (Bakk & Kuha, 2018) and/or distal outcome (Bakk, Tekle & Vermunt) means with either the ML correction (Vermunt, 2010) or the BCH correction (Bolck, Croon & Hagenaars, 2004).
Main functions
three_stepFull three-step estimation pipeline. Accepts covariates (
Zp.names), distal outcomes (Zo.name), or both. Handles Steps 1–3 in a single call, with optional pre-fitted Step-1 input throughstep1.lca_step1Standalone Step-1 measurement model estimation via multilevLCA. Returns a reusable fit object that can be passed to
three_step()to avoid re-estimating the measurement model across multiple structural specifications.fitZ_from_fit0Two-step covariate estimation by fixing measurement parameters at their Step-1 values and estimating multinomial logit coefficients
\gammawith an EM algorithm. Returns starting values for Step 3. Custom starting values can be supplied viastarting_val.fitZ_from_multiLCATwo-step covariate estimation via
multiLCA(fixedpars = 1), returning multilevLCA's bias-corrected standard errors. Called automatically whenget.twostep.vcov = TRUEinthree_step.generate_dataSimulate data replicating the Bakk & Kuha (2018) three-class design: six binary indicators across three separation levels (
"low","medium","high") and two scenarios ("covariate"or"distal").
Estimators
- ML (default,
use.bch = FALSE) The Vermunt (2010) ML correction uses a weighted pseudo-likelihood with the misclassification matrix as a bias adjustment. Preferred when class separation is low or moderate.
- BCH (
use.bch = TRUE) The Bolck, Croon & Hagenaars (2004) correction inverts the misclassification matrix to obtain direct class weights. Works well under high separation but may produce an ill-conditioned Hessian (non-positive semi-definite covariance matrix) when separation is low; use the ML estimator in that case.
Variance estimation
- Full correction (
use.simple.cov = FALSE, default) Analytic propagation of Step-1 measurement uncertainty through the classification-error correction, following Bakk, Oberski & Vermunt (2014). Uses soft (proportional) posteriors for the Jacobian
\partial\theta_2/\partial\theta_1regardless ofuse.modal.assignment. Recommended when separation is moderate or low.- Simple/robust (
use.simple.cov = TRUE) Sandwich SEs from Step 3 only, ignoring measurement uncertainty. A useful computational shortcut when separation is high and the correction is negligible.
Class assignment
- Modal (
use.modal.assignment = TRUE, default) Each observation is assigned to its most probable class (hard assignment). The Jacobian for variance correction is still computed from soft posteriors.
- Proportional (
use.modal.assignment = FALSE) Soft posterior weights are used throughout Steps 2 and 3. Recommended when separation is moderate or low, and required for a mathematically well-defined analytic Jacobian.
Supported features
Binary and polytomous indicators, following multilevLCA coding conventions.
Gaussian, Poisson, and binomial distal outcome families.
Full-information maximum likelihood (FIML) for partially observed indicator patterns (
incomplete = TRUE). Step 3 always performs listwise deletion on missing covariates or distal outcomes.Flexible measurement and structural samples: fit the measurement model on a reference sample and apply it to a different analysis sample via the
step1argument.Arbitrary reference class for the multinomial logit parameterization via the
rebaseargument. Log-likelihoods are invariant to this choice.Joint covariate and distal outcome estimation (
Zp.namesandZo.namesupplied together). The covariate model is estimated first; covariate-adjusted posteriors are then used as priors in the distal outcome step.S3 methods (
print,summary,coef,vcov,plot) for all four return subclasses:tseLCA_measurement,tseLCA_covariate,tseLCA_distal,tseLCA_both.
Getting started
# Introductory vignette
vignette("tseLCA-workflow", package = "tseLCA")
Author(s)
Sam Lee samlee@arizona.edu
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. doi:10.1177/0081175012470644
Bakk, Z., Oberski, D. L., & Vermunt, J. K. (2014). Relating latent class assignments to external variables: Standard errors for correct inference. Political Analysis, 22(4), 520–540. https://www.jstor.org/stable/24573086
Bakk, Z., & Kuha, J. (2018). Two-step estimation of models between latent classes and external variables. Psychometrika, 83(4), 871–892. doi: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. doi: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. doi: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. doi:10.1093/pan/mpq025
See Also
Useful links:
Report bugs at https://github.com/SamLeeBYU/tseLCA/issues
Default population parameters for the Bakk & Kuha (2018) simulation
Description
A list of pre-specified parameters used by generate_data() and related
functions. All elements correspond to the values stated in the paper
(Section 3, p. 879).
Usage
bk2018_params
Details
class_propsLength-3 vector of equal class proportions (1/3 each).
separation_levelsNamed vector mapping
"low","mid","high"to the probability of a "likely" response (0.70, 0.80, 0.90).covariate_paramsList with
$b0(intercepts) and$b(slopes) for the multinomial logit P(X=t | Zp). Interceptsb02andb03are set so that marginal class sizes average to 1/3 when Zp ~ Uniform{1..5}.distal_paramsList with
$mu(class means, c(-1, 1, 0)) and$sigma(residual SD = 1) for the distal outcome model.
Value
A named list with four elements:
class_propsLength-3 numeric vector of equal class proportions (1/3 each).
separation_levelsNamed numeric vector mapping
"low","mid","high"to 0.70, 0.80, 0.90.covariate_paramsList with
$b0(intercepts) and$b(slopes) for the multinomial logit P(X=t|Zp).distal_paramsList with
$mu(class means) and$sigma(residual SD).
Examples
# True item-response probabilities for high separation
bk2018_params$rho_high
# Covariate model parameters (intercepts and slopes)
bk2018_params$covariate_params
# Distal outcome parameters
bk2018_params$distal_params
Prepare and validate data for tseLCA estimation
Description
Prepare and validate data for tseLCA estimation
Usage
clean_data(
data,
Y.names,
Zp.names = NULL,
Zo.name = NULL,
incomplete = FALSE,
include.intercept = TRUE,
verbose = FALSE
)
Arguments
data |
A data.frame. |
Y.names |
Character vector of item column names. |
Zp.names |
Character vector of covariate column names, or |
Zo.name |
Single distal outcome column name, or |
incomplete |
Logical. If |
include.intercept |
Logical. Prepend intercept column to Z. |
verbose |
Logical. Print row-drop messages. |
Value
A named list with:
- Y.obs
N_Y x K expanded one-hot indicator matrix for Steps 1 & 2.
- mDesign
N_Y x K design/mask matrix (all 1s when incomplete = FALSE).
- ivItemcat
Integer vector of category counts per item.
- keep_Y
Integer indices of rows kept for Steps 1 & 2 (into original N).
- Z_mat
N_Z x Q covariate design matrix, or NULL.
- keep_step3_Z_in_Y
Positions of Z-complete rows within keep_Y.
- Zo_mat
N_Zo x 1 distal outcome matrix, or NULL.
- keep_step3_Zo_in_Y
Positions of Zo-complete rows within keep_Y.
Extract coefficients from a tseLCA model object
Description
Extract coefficients from a tseLCA model object
Usage
## S3 method for class 'tseLCA_measurement'
coef(object, ...)
## S3 method for class 'tseLCA_covariate'
coef(object, which = c("three_step", "two_step"), ...)
## S3 method for class 'tseLCA_distal'
coef(object, ...)
## S3 method for class 'tseLCA_both'
coef(
object,
which = c("both", "covariate", "distal"),
step = c("three_step", "two_step"),
...
)
Arguments
object |
A |
... |
Further arguments (currently unused). |
which |
Character. For covariate and both models: |
step |
Character. For |
Value
The coefficient matrix (covariate models), named numeric vector (distal models), or a named list of both (measurement or both models).
Examples
d <- generate_data(100, "high", "covariate", seed = 1)
fit_m <- three_step(d, paste0("Y", 1:6), n_classes = 3)
coef(fit_m) # returns list with $prevalences and $item_probs
d <- generate_data(200, "high", "covariate", seed = 1)
fit <- three_step(d, paste0("Y", 1:6), n_classes = 3,
Zp.names = "Zp", use.simple.cov = TRUE)
coef(fit) # three-step estimates
coef(fit, which = "two_step") # two-step starting values
d <- generate_data(200, "high", "distal", seed = 2)
fit <- three_step(d, paste0("Y", 1:6), n_classes = 3,
Zo.name = "Zo", use.simple.cov = TRUE)
coef(fit) # named vector of class means
d <- generate_data(200, "high", "covariate", seed = 1)
d$Zo <- rnorm(200, mean = c(-1, 0, 1)[d$X], sd = 1)
fit <- three_step(d, paste0("Y", 1:6), n_classes = 3,
Zp.names = "Zp", Zo.name = "Zo",
use.simple.cov = TRUE)
coef(fit, which = "covariate")
coef(fit, which = "distal")
coef(fit, which = "both")
Compress a one-hot expanded Y matrix back to integer codes
Description
Inverse of expand_Y. Takes a one-hot expanded matrix where each
item occupies K_h consecutive columns (one per category, 0-based)
and returns an N x H integer matrix of category codes (0, 1, ..., K_h-1).
Usage
compress_Y(mY_exp, ivItemcat)
Arguments
mY_exp |
N x sum(K_h) one-hot matrix (as stored in |
ivItemcat |
Integer vector of category counts per item (length H). |
Details
Rows where all K_h columns for an item are NA are returned as
NA for that item.
Value
N x H integer matrix of category codes.
Draw a continuous distal outcome given true class memberships (scenario "distal")
Description
Draw a continuous distal outcome given true class memberships (scenario "distal")
Usage
draw_Zo(X, params)
Arguments
X |
Integer vector of length n. Latent class (1-indexed). |
params |
List with |
Value
Numeric vector of length n.
Examples
X <- draw_classes(100, c(1/3, 1/3, 1/3))
Zo <- draw_Zo(X, bk2018_params$distal_params)
tapply(Zo, X, mean) # should be close to true mu
Draw the covariate Zp ~ Uniform{1, 2, 3, 4, 5}
Description
Draw the covariate Zp ~ Uniform{1, 2, 3, 4, 5}
Usage
draw_Zp(n)
Arguments
n |
Integer. Sample size. |
Value
Integer vector of length n.
Examples
Zp <- draw_Zp(100)
table(Zp)
Draw latent class memberships from their marginal distribution
Description
Draw latent class memberships from their marginal distribution
Usage
draw_classes(n, pi_)
Arguments
n |
Integer. Sample size. |
pi_ |
Numeric vector of length T. Class proportions (must sum to 1). |
Value
Integer vector of length n with values in 1:T.
Examples
# Draw 100 class labels from equal prevalences
draw_classes(100, c(1/3, 1/3, 1/3))
Draw latent classes conditional on the covariate (scenario "covariate")
Description
Draw latent classes conditional on the covariate (scenario "covariate")
Usage
draw_classes_given_Zp(Zp, params)
Arguments
Zp |
Numeric vector of length n. Covariate values. |
params |
Multinomial logistic parameter list (see |
Value
Integer vector of length n with class labels in 1:T.
Examples
Zp <- draw_Zp(1000)
X <- draw_classes_given_Zp(Zp, bk2018_params$covariate_params)
table(X) # Should be roughly uniform
Draw binary indicators given true class memberships
Description
Draw binary indicators given true class memberships
Usage
draw_indicators(X, rho)
Arguments
X |
Integer vector of length n. True latent class (1-indexed). |
rho |
T x K matrix. |
Value
An n x K integer matrix of 0/1 values.
Examples
rho <- make_rho(0.9)
X <- draw_classes(50, c(1/3, 1/3, 1/3))
draw_indicators(X, rho)
Extract Y.exp, mDesign, posteriors from a multilevLCA mU matrix
Description
fit0$mU from multilevLCA stores data already in one-hot expanded
form: each item h occupies K_h consecutive columns (one per category),
followed by T columns of posterior class probabilities.
Usage
extract_Y_from_mU(fit0, ivItemcat = NULL)
Arguments
fit0 |
Raw multilevLCA fit object with |
ivItemcat |
Integer vector of category counts per item (length H).
If |
Details
For dichotomous items (K_h=2) the two columns are stored. For polytomous
items (K_h>2) all K_h columns are stored. This function first compresses
the expanded Y back to integer codes through compress_Y, then re-expands
consistently via expand_Y so downstream functions receive the correct
N x K_total matrix.
Value
A list with:
- Y.exp
N x K_total expanded one-hot matrix (NAs replaced with 0).
- mDesign
N x K_total design/mask matrix.
NULLif no missing.- ivItemcat
Integer vector of category counts per item.
- u_post
N x T posterior class probability matrix from
mU.
Estimate covariate effects with measurement parameters fixed (two-step EM)
Description
Fixes mPhi at fit0$mPhi and estimates multinomial logit coefficients
mGamma (Q x (T-1)) via an EM algorithm with a BFGS M-step.
Usage
fitZ_from_fit0(
fit0,
data,
Y.names,
Zp.names,
tol = 1e-06,
maxIter = 200L,
incomplete = FALSE,
include.intercept = TRUE,
rebase = "C1",
starting_val = NULL,
verbose = FALSE
)
Arguments
fit0 |
Output of |
data |
A data.frame. |
Y.names |
Character vector of item column names. |
Zp.names |
Character vector of covariate column names. |
tol |
Convergence tolerance. Default |
maxIter |
Maximum EM iterations. Default |
incomplete |
Logical. FIML for partially missing indicators. See the
|
include.intercept |
Logical. Prepend intercept to covariate design matrix. Default |
rebase |
Character or integer. Reference class for the multinomial logit
parameterization (e.g. |
starting_val |
Optional Q x (T-1) starting value matrix for |
verbose |
Logical. Print convergence messages. Default |
Value
A list with the following elements:
mGammaQ x (T-1) numeric matrix of multinomial logit coefficients, where Q is the number of columns in the covariate design matrix (including intercept if
include.intercept = TRUE). Rows are named by covariate, columns by non-reference class (e.g."C2","C3").mPhiExpanded item parameter matrix (items x classes), fixed at
fit0$mPhithroughout estimation.vOmegaLength-T vector of marginal class proportions implied by the final
mGamma, computed as column means of the fitted class probability matrix.LLKSeriesSingle-column matrix of observed-data log-likelihoods, one row per EM iteration. Useful for diagnosing convergence.
convergedLogical.
TRUEif the EM loop exited beforemaxIteriterations or if the final log-likelihood change was belowtol.n_obsInteger. Number of observations used in estimation after listwise deletion on covariates.
Examples
d <- generate_data(200, "high", "covariate", seed = 1)
s1 <- lca_step1(d, Y.names = paste0("Y", 1:6), n_classes = 3)
# Estimate two-step gamma with mPhi fixed at Step-1 values
fZ <- fitZ_from_fit0(
fit0 = s1$fit0,
data = d,
Y.names = paste0("Y", 1:6),
Zp.names = "Zp",
verbose = TRUE
)
fZ$mGamma # Q x (T-1) coefficient matrix
fZ$converged
Estimate two-step covariate model via multilevLCA (optional reference path)
Description
Calls multilevLCA::multiLCA with fixedpars = 1 and Z = Zp.names to
fit the two-step covariate model. This is the original multilevLCA approach
and is used when get.twostep.vcov = TRUE in three_step() to obtain
multilevLCA's corrected standard errors for the two-step gamma estimates.
Usage
fitZ_from_multiLCA(
data,
Y.names,
n_classes,
Zp.names,
maxIter.measurement,
measurement.tol,
covariate.tol,
iter.measurement,
R2.threshold,
incomplete = FALSE,
rebase = "C1",
verbose = FALSE
)
Arguments
data |
A data.frame. |
Y.names |
Character vector of item column names. |
n_classes |
Integer. Number of latent classes. |
Zp.names |
Character vector of covariate column names. |
maxIter.measurement |
Maximum EM iterations. |
measurement.tol |
Convergence tolerance. |
covariate.tol |
NR tolerance for the covariate model. |
iter.measurement |
Number of random restarts. |
R2.threshold |
Entropy R |
incomplete |
Logical. FIML for partially missing indicators. See the
|
rebase |
Character or integer. Reference class for column naming of
|
verbose |
Logical. |
Value
A list with the following elements:
mGammaQ x (T-1) numeric matrix of multinomial logit coefficients. Rows are named by covariate (including
"Intercept"), columns by non-reference class (e.g."C2","C3").mPhiItem parameter matrix (items x classes) from the fixed-parameter multilevLCA fit.
vOmegaLength-T vector of marginal class proportions, computed as the average of the fitted class probability matrix (
vPi_avgin multilevLCA output).LLKSeriesMatrix of observed-data log-likelihoods across EM iterations, passed through directly from the multilevLCA fit.
raw_fitThe full
multilevLCA::multiLCA()output object, including$Varmat_cor(corrected variance matrix) and$SEs_cor_gamma(corrected standard errors formGamma) if available.
Examples
d <- generate_data(200, "high", "covariate", seed = 1)
# Two-step estimation via multiLCA (fixedpars = 1)
fZ_ml <- fitZ_from_multiLCA(
data = d,
Y.names = paste0("Y", 1:6),
n_classes = 3,
Zp.names = "Zp",
maxIter.measurement = 5000L,
measurement.tol = 1e-8,
covariate.tol = 1e-6,
iter.measurement = 10L,
R2.threshold = 0.70
)
fZ_ml$mGamma # two-step estimates
fZ_ml$raw_fit$Varmat_cor # multilevLCA corrected vcov
Generate datasets for all 18 conditions in the simulation design
Description
Iterates over the 2 scenarios x 3 separation levels x 3 sample sizes,
generating n_rep independent replications per condition. Seeds are
derived deterministically from base_seed so the entire experiment is
reproducible from a single integer.
Usage
generate_all_conditions(
n_rep = 500L,
base_seed = 5262026L,
params = bk2018_params,
scenarios = c("covariate", "distal"),
sep_levels = c("low", "mid", "high"),
sample_sizes = c(500L, 1000L, 2000L),
verbose = TRUE
)
Arguments
n_rep |
Integer. Replications per condition (paper uses 500). |
base_seed |
Integer. Base seed for reproducibility. |
params |
Population parameters list. Defaults to bk2018_params. |
scenarios |
Character. Lists the scenario(s) ("covariate" and/or "distal")
wanting to be simulated. Passed into |
sep_levels |
Character. Lists the separation level(s) ("low", "mid", "high")
wanting to be simulated. Passed into |
sample_sizes |
Integer. Lists the sample size(s) wanting to be generated
for each replication condition. Passed into |
verbose |
Logical. If |
Value
Nested list indexed as
datasets[[scenario]][[separation]][[as.character(n)]],
each element a list of n_rep data frames.
Examples
# Generate 5 replicates for mid and high separation only
datasets <- generate_all_conditions(n_rep = 5L, base_seed = 1L,
sep_levels = c("mid", "high"))
# Access a single replicate
head(datasets[["covariate"]][["high"]][["500"]][[1]])
Generate one dataset following the Bakk & Kuha (2018) simulation design
Description
Generate one dataset following the Bakk & Kuha (2018) simulation design
Usage
generate_data(
n,
separation = c("low", "mid", "high"),
scenario = c("covariate", "distal"),
params = bk2018_params,
seed = NULL
)
Arguments
n |
Integer. Sample size (paper uses 500, 1000, or 2000). |
separation |
Character. One of |
scenario |
Character. One of:
|
params |
List of population parameters. Defaults to bk2018_params. |
seed |
Integer or |
Value
A data.frame with columns:
Y1..Y6Binary indicators (always present).
XTrue latent class, integer 1-3 (not observed in practice).
ZpInteger covariate 1-5 (scenario
"covariate"only).ZoContinuous distal outcome (scenario
"distal"only).
Examples
# Covariate scenario with high separation
d <- generate_data(n = 200, separation = "high", scenario = "covariate",
seed = 1)
head(d)
colMeans(d)
# Distal outcome scenario
d2 <- generate_data(n = 200, separation = "high", scenario = "distal",
seed = 2)
head(d2)
Individual-level BHHH variance matrix for binary and polytomous LCA
Description
Computes the outer-product (BHHH) information matrix and variance-covariance
matrix for LCA measurement model parameters in the unconstrained
(logit/log-ratio) space, matching multilevLCA's $Varmat.
Usage
lca_indiv_varmat(
Y.exp,
mDesign.exp,
fit0,
ivItemcat,
boundary.tol = 0.01,
use.freq = TRUE,
u_post = NULL
)
Arguments
Y.exp |
N x sum(K_h) expanded one-hot indicator matrix. |
mDesign.exp |
Expanded design matrix (same dimensions as |
fit0 |
Step-1 fit object with |
ivItemcat |
Integer vector of category counts per item. |
boundary.tol |
Scalar tolerance for boundary detection. Default
|
use.freq |
Logical. Collapse duplicate score rows before computing
the cross-product, weighting by frequency. Default |
u_post |
Optional N x T matrix of posterior class probabilities.
When supplied (e.g. extracted from |
Details
The score in unconstrained space is
s_{it} = u_{it}(y_i - d_i \circ p_{it}),
where d_i is the missing-data design indicator matrix.
Assumes fit0$mPhi follows the multilevLCA storage convention:
Dichotomous item h (
ivItemcat[h] == 2): 1 row =P(Y=1|C); the base levelP(Y=0|C)is excluded.Polytomous item h (
ivItemcat[h] > 2):K_hrows =P(Y=0|C), \ldots, P(Y=K_h-1|C); the base level is included.
expand_Y produces one-hot columns in the same order so that
expand_Phi(fit0$mPhi, ivItemcat) aligns column-wise with
expand_Y(mY, ivItemcat). Free (estimable) parameters per item are
the single P(Y=1|C) row for dichotomous items, and rows 2 through
K_h for polytomous items (row 1, P(Y=0|C), is the reference).
Boundary parameters (within boundary.tol of 0 or 1) are treated as
fixed: their score columns are zeroed and they do not contribute to the
information matrix.
Value
A list with the following elements:
InfomatSquare BHHH information matrix of dimension p x p, where p = (iT-1) + sum(ivItemcat - 1) * iT is the total number of free parameters. Boundary parameters have zero rows and columns.
VarmatInverse of
Infomatdivided by N, giving the asymptotic variance-covariance matrix on the same scale as multilevLCA's$Varmat. Boundary parameters have zero rows and columns.SEsNumeric vector of length p. Square root of the diagonal of
Varmat; zero for boundary parameters.mScoreN x p matrix of individual score contributions in the unconstrained parameterization, used for sandwich variance propagation in
lca_vcovandlca_vcov_distal.
Fit the LCA measurement model (Step 1)
Description
Estimates the latent class measurement model with multilevLCA and
optionally, fixes mPhi and estimates covariate effects (two-step
initialization) with fitZ_from_fit0().
Usage
lca_step1(
data,
Y.names,
n_classes,
Zp.names = NULL,
maxIter.measurement = 5000L,
measurement.tol = 1e-08,
covariate.tol = 1e-06,
iter.measurement = 10L,
R2.threshold = 0.7,
use.two.step = TRUE,
estimate.one.step = TRUE,
incomplete = FALSE,
maxIter.fitZ = 200L,
include.intercept = TRUE,
rebase = "C1",
verbose = FALSE
)
Arguments
data |
A data.frame containing at minimum the indicator columns. |
Y.names |
Character vector of item column names. |
n_classes |
Integer. Number of latent classes. |
Zp.names |
Character vector of covariate column names, or |
maxIter.measurement |
Maximum EM iterations before giving up on convergence. Default |
measurement.tol |
Convergence tolerance. Default |
covariate.tol |
Convergence tolerance for the |
iter.measurement |
Number of random restarts when entropy R |
R2.threshold |
Entropy R |
use.two.step |
Logical. If |
estimate.one.step |
Logical. If |
incomplete |
Logical. FIML for partially missing indicators. See the
|
maxIter.fitZ |
Maximum EM iterations for |
include.intercept |
Logical. Prepend intercept to covariate design matrix. Default |
rebase |
Character or integer specifying the reference latent class.
Use |
verbose |
Logical. Print progress messages. Default |
Value
A list with $fit0 (multilevLCA::multiLCA() measurement model) and $fitZ
(two-step covariate model from fitZ_from_fit0(), or NULL).
Examples
d <- generate_data(200, "high", "covariate", seed = 1)
# Measurement model only
s1 <- lca_step1(d, Y.names = paste0("Y", 1:6), n_classes = 3)
s1$fit0$vPi # estimated class prevalences
s1$fit0$mPhi # item-response probabilities
# With two-step covariate initialization
s1z <- lca_step1(d, Y.names = paste0("Y", 1:6), n_classes = 3,
Zp.names = "Zp", use.two.step = TRUE, verbose = TRUE)
s1z$fitZ$mGamma # two-step gamma estimates
Build the item-response probability matrix for the simulation
Description
Returns a T x K matrix rho where rho[t, k] = P(Y_k = 1 | X = t).
Usage
make_rho(pi_)
Arguments
pi_ |
Numeric scalar in (0.5, 1). Probability of the "likely" response.
Use |
Details
The three-class structure is:
Class 1:
pi_on all 6 items (high responders).Class 2:
pi_on items 1-3,1 - pi_on items 4-6 (mixed).Class 3:
1 - pi_on all 6 items (low responders).
Value
A 3 x 6 numeric matrix.
Examples
# High separation: P(Y=1|class) = 0.9 for the "high" class
make_rho(0.9)
# Low separation
make_rho(0.7)
Compute multinomial logistic class probabilities given covariates
Description
Evaluates P(X = t | Zp) for each observation using a multinomial logit with one or more covariates and class-specific intercepts and slopes.
Usage
mnl_probs(Zp, params)
Arguments
Zp |
Numeric vector of length n, or numeric matrix of dimension n x P, where P is the number of covariates. A vector is treated as a single covariate (P = 1). |
params |
List with elements |
Value
An n x T matrix of class probabilities (rows sum to 1).
Examples
# Single covariate: class membership probabilities for Zp = 1..5
mnl_probs(1:5, bk2018_params$covariate_params)
# Multiple covariates (n = 5, P = 2)
Zp_mat <- matrix(rnorm(10), nrow = 5, ncol = 2)
params2 <- list(b0 = c(0, 0.5, -0.5), b = matrix(rnorm(6), nrow = 2, ncol = 3))
mnl_probs(Zp_mat, params2)
Normalize row/column names of a fitZ$mGamma matrix
Description
A plain multiLCA object uses rownames like "gamma(Intercept|C)" and
"gamma(Zp|C)". This function strips the gamma(...) wrapper so names
match the clean format used throughout tseLCA ("Intercept", "Zp", etc.)
and ensures column names are "C2", "C3", etc.
Usage
normalize_fitZ_names(fitZ, Zp.names = NULL, n_classes = NULL)
Arguments
fitZ |
A fitZ-like list with at least |
Zp.names |
Character vector of covariate column names (used to set
clean rownames when the raw names can't be parsed). If |
n_classes |
Integer. Total number of classes (used to derive clean column names if they are non-standard). |
Value
fitZ with normalized $mGamma row/col names.
Parse and validate the rebase argument
Description
Parse and validate the rebase argument
Usage
parse_rebase(rebase, iT)
Arguments
rebase |
Character like "C2" or integer class index. |
iT |
Total number of classes. |
Value
Integer class index (1-based) to use as reference.
Permute class columns of a fit0 object so that class ref_idx is first
Description
Reorders columns of mPhi and vPi so that the desired reference class becomes column 1 before estimation. This ensures the multinomial logit is parameterized with the correct baseline from the start.
Usage
permute_fit0_classes(fit0, ref_idx)
Arguments
fit0 |
Raw multilevLCA fit object (has $mPhi and $vPi). |
ref_idx |
Integer. Class index to move to position 1. |
Value
fit0 with columns permuted.
Permute class columns of a fitZ object to match a new reference class
Description
Rebases a fitZ object (output of fitZ_from_fit0 or
fitZ_from_multiLCA) so that ref_idx becomes the reference class.
This involves:
Rebasing
$mGamma: reconstructing the full T-column log-ratio matrix, subtracting the new reference column, and dropping it.Propagating through
$Varmat_corvia the delta method: the rebasing transformation is linear (gamma_new = A * gamma_old) so the vcov transforms exactly asA %*% V %*% t(A).Updating all column names.
Usage
permute_fitZ_classes(fitZ, ref_idx)
Arguments
fitZ |
Output of |
ref_idx |
Integer. New reference class (1-based index into the T
classes as currently ordered in |
Value
fitZ with $mGamma, $Varmat_cor, and names updated.
Plot item-response probability profiles for a tseLCA model
Description
Delegates to plot.multiLCA from multilevLCA, which draws the
class-specific item-response probability profiles from the Step-1
measurement model.
Usage
## S3 method for class 'tseLCA_measurement'
plot(x, horiz = FALSE, clab = NULL, ...)
## S3 method for class 'tseLCA_covariate'
plot(x, horiz = FALSE, clab = NULL, ...)
## S3 method for class 'tseLCA_distal'
plot(x, horiz = FALSE, clab = NULL, ...)
## S3 method for class 'tseLCA_both'
plot(x, horiz = FALSE, clab = NULL, ...)
Arguments
x |
A |
horiz |
Logical. If |
clab |
Optional character vector of length T giving class labels. |
... |
Further arguments passed to |
Value
Called for its side effect (a base-graphics plot). Invisibly
returns NULL.
Examples
d <- generate_data(100, "high", "covariate", seed = 1)
fit_m <- three_step(d, paste0("Y", 1:6), n_classes = 3)
plot(fit_m)
# Custom class labels
plot(fit_m, clab = c("Low risk", "Mixed", "High risk"))
Print a tseLCA model object
Description
Compact one-line or table summary printed to the console.
Usage
## S3 method for class 'tseLCA_measurement'
print(x, ...)
## S3 method for class 'tseLCA_covariate'
print(x, digits = 4, ...)
## S3 method for class 'tseLCA_distal'
print(x, digits = 4, ...)
## S3 method for class 'tseLCA_both'
print(x, digits = 4, ...)
Arguments
x |
A |
... |
Further arguments. |
digits |
Integer. Number of decimal places for coefficient tables. |
Value
Invisibly returns x.
Examples
d <- generate_data(100, "high", "covariate", seed = 1)
fit_m <- three_step(d, paste0("Y", 1:6), n_classes = 3)
print(fit_m)
d <- generate_data(200, "high", "covariate", seed = 1)
fit <- three_step(d, paste0("Y", 1:6), n_classes = 3,
Zp.names = "Zp", use.simple.cov = TRUE)
print(fit)
d <- generate_data(200, "high", "distal", seed = 2)
fit <- three_step(d, paste0("Y", 1:6), n_classes = 3,
Zo.name = "Zo", use.simple.cov = TRUE)
print(fit)
Summarize a tseLCA model object
Description
Verbose summary including model fit, class prevalences, item-response probabilities, and coefficient tables with standard errors and p-values.
Usage
## S3 method for class 'tseLCA_measurement'
summary(object, ...)
## S3 method for class 'tseLCA_covariate'
summary(object, digits = 4, ...)
## S3 method for class 'tseLCA_distal'
summary(object, digits = 4, ...)
## S3 method for class 'tseLCA_both'
summary(object, digits = 4, ...)
Arguments
object |
A |
... |
Further arguments (currently unused). |
digits |
Integer. Number of decimal places for coefficient tables. |
Value
Invisibly returns object.
Examples
d <- generate_data(100, "high", "covariate", seed = 1)
fit_m <- three_step(d, paste0("Y", 1:6), n_classes = 3)
summary(fit_m)
d <- generate_data(200, "high", "covariate", seed = 1)
fit <- three_step(d, paste0("Y", 1:6), n_classes = 3,
Zp.names = "Zp", use.simple.cov = TRUE)
summary(fit)
d <- generate_data(200, "high", "distal", seed = 2)
fit <- three_step(d, paste0("Y", 1:6), n_classes = 3,
Zo.name = "Zo", use.simple.cov = TRUE)
summary(fit)
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,
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 |
Character vector of indicator column names. Need to be coded as consecutive integers with base level starting at |
n_classes |
Integer. Number of latent classes. |
Zp.names |
Character vector of covariate column names, or |
Zo.name |
Single character name of the distal outcome column, or
|
step1 |
Pre-fitted Step-1 object (output of |
use.two.step |
Logical. Initialize Step-3 from two-step estimates.
Default |
use.modal.assignment |
Logical. Use modal (hard) class assignments in
Step 2 and 3. |
include.intercept |
Logical. Prepend an intercept column to the
covariate design matrix. Default |
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
|
incomplete |
Logical. FIML for partially missing indicators. See the
|
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
|
maxIter.measurement |
Integer. Maximum EM iterations for Step 1.
Default |
measurement.tol |
Scalar. Convergence tolerance for the Step-1 EM
algorithm. Default |
covariate.tol |
Scalar. Convergence tolerance for the Step-3
Newton-Raphson or EM algorithm. Default |
iter.measurement |
Integer. Number of random restarts triggered when
the Step-1 entropy R |
R2.threshold |
Scalar. Entropy R |
use.bch |
Logical. Use BCH-corrected weights instead of the ML
estimator in Step 3. May error if BCH weights induce a non-positive semi-definite Hessian in the third step (common in cases of low separation). Default |
em.maxIter |
Integer. Maximum EM iterations for the Step-3 covariate
or distal outcome model. Default |
get.twostep.vcov |
Logical. If |
rebase |
Character (e.g. |
family |
Character. Distal outcome family: one of |
correct.spec |
Logical. Use the model-robust outer-product Hessian for
Step-3 standard errors rather than the observed-data Hessian. Not appropriate
when the Step-3 model may be misspecified. Default |
verbose |
Logical. Print convergence messages. Default |
Value
An S3 object of class tseLCA. The subclass depends on which
models were estimated:
tseLCA_measurementReturned when neither
Zp.namesnorZo.nameis supplied. Contains the following elements:measurement_modelStep-1 output list from
lca_step1().llikFinal Step-1 log-likelihood.
AIC,BICInformation criteria from the measurement model.
R2entrEntropy R
^2of the measurement model.n_classesNumber of latent classes.
posteriorsN x T matrix of soft posterior class probabilities.
classificationsLength-N integer vector of modal class assignments.
tseLCA_covariateReturned when
Zp.namesis supplied andZo.nameisNULL. Contains all elements oftseLCA_measurementplus:three_stepQ x (T-1) matrix of Step-3 gamma coefficients.
three_step_vcovQ(T-1) x Q(T-1) variance-covariance matrix for
three_step, with measurement-uncertainty correction unlessuse.simple.cov = TRUE.two_stepQ x (T-1) matrix of two-step starting values, or
NULLifuse.two.step = FALSE.two_step_vcovmultilevLCA bias-corrected vcov for the two-step estimates, or
NULL.estimatorCharacter:
"ML"or"BCH".entropy.R2Covariate-adjusted entropy R
^2.llikProfile 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_distalReturned when
Zo.nameis supplied andZp.namesisNULL. Contains:three_stepNamed length-T vector of Step-3 distal outcome parameters (means, log-rates, or logits depending on
family).three_step_vcovT x T variance-covariance matrix for
three_step, namedmu_C1throughmu_CT.three_step.llikStep-3 distal log-likelihood
\log P(Z_o|X=t)at converged estimates.llikProfile 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.AICAkaike information criterion based on
llik.BICBayesian information criterion based on
llik, using the number of distal-complete observations.familyCharacter. The distal outcome family used.
estimatorCharacter:
"ML"or"BCH".posteriorsN x T soft posterior matrix.
classificationsLength-N modal class assignment vector.
tseLCA_bothReturned when both
Zp.namesandZo.nameare supplied. Contains:covariateA
tseLCA_covariate-structured sub-list (see above), includingllik,AIC,BIC,entropy.R2.distalA
tseLCA_distal-structured sub-list (see above), includingllik,AIC,BIC,three_step.llik.family,n_classes,estimatorShared top-level fields.
posteriors,classificationsShared N x T posterior matrix and length-N modal class vector.
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. 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. 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. doi:10.1080/00273171.2025.2473935
See Also
vignette("tseLCA", package = "tseLCA") for a full worked
example; lca_step1() for standalone Step-1 estimation;
fitZ_from_fit0() and fitZ_from_multiLCA() for two-step covariate
estimation.
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)
# 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)
# Plot item-response profiles from the measurement model
plot(fit)
Extract the variance-covariance matrix from a tseLCA model object
Description
For measurement models, returns the BHHH variance-covariance matrix in
the unconstrained log-ratio parameterization (NOT the probability scale).
Row and column names identify each parameter as
log(pi_t/pi_1) (class prevalences) or
log(P(Y=k|C_t)/P(Y=0|C_t)) (item-response probabilities).
An attribute "parameterization" is attached to remind the user of the
scale.
Usage
## S3 method for class 'tseLCA_measurement'
vcov(object, boundary.tol = 0.01, ...)
## S3 method for class 'tseLCA_covariate'
vcov(object, which = c("three_step", "two_step"), ...)
## S3 method for class 'tseLCA_distal'
vcov(object, ...)
## S3 method for class 'tseLCA_both'
vcov(
object,
which = c("both", "covariate", "distal"),
step = c("three_step", "two_step"),
...
)
Arguments
object |
A |
boundary.tol |
Scalar. Parameters within this tolerance of 0 or 1
are treated as fixed. Default |
... |
Further arguments (currently unused). |
which |
Character. |
step |
Character. For |
Value
A named square matrix in the unconstrained log-ratio
parameterization. Row/column names identify each parameter as
log(pi_t/pi_1) or log(P(Y=k|C_t)/P(Y=0|C_t)). An attribute
"parameterization" is attached as a reminder. Returns NULL
invisibly if fit0$mU is not available. For structural models,
returns the Step-3 vcov matrix; the two-step vcov is only available
when get.twostep.vcov = TRUE.
Examples
d <- generate_data(100, "high", "covariate", seed = 1)
fit_m <- three_step(d, paste0("Y", 1:6), n_classes = 3)
V <- vcov(fit_m)
# Names show log-ratio parameterization:
rownames(V)
attr(V, "parameterization")
d <- generate_data(200, "high", "covariate", seed = 1)
fit <- three_step(d, paste0("Y", 1:6), n_classes = 3,
Zp.names = "Zp", use.simple.cov = TRUE)
vcov(fit) # Q*(T-1) x Q*(T-1) vcov matrix with named rows/cols
d <- generate_data(200, "high", "distal", seed = 2)
fit <- three_step(d, paste0("Y", 1:6), n_classes = 3,
Zo.name = "Zo", use.simple.cov = TRUE)
vcov(fit) # T x T vcov matrix with mu_C1..mu_CT row/col names
d <- generate_data(200, "high", "covariate", seed = 1)
d$Zo <- rnorm(200, mean = c(-1, 0, 1)[d$X], sd = 0.5)
fit <- three_step(d, paste0("Y", 1:6), n_classes = 3,
Zp.names = "Zp", Zo.name = "Zo",
use.simple.cov = TRUE)
vcov(fit, which = "covariate")
vcov(fit, which = "distal")