| Title: | Three-Step Estimation for Latent Class Analysis |
|---|---|
| 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. |
| Authors: | Sam Lee [aut, cre, cph] (ORCID: <https://orcid.org/0009-0007-2977-6661>), Jay Goodliffe [ctb] |
| Maintainer: | Sam Lee <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 1.0.2 |
| Built: | 2026-07-11 10:58:31 UTC |
| Source: | https://github.com/cran/tseLCA |
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).
bk2018_paramsbk2018_params
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). Intercepts b02 and b03 are
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.
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).
# 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# 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
Extract coefficients from a tseLCA model object
## 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"), ... )## 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"), ... )
object |
A |
... |
Further arguments (currently unused). |
which |
Character. For covariate and both models: |
step |
Character. For |
The coefficient matrix (covariate models), named numeric vector (distal models), or a named list of both (measurement or both models).
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")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")
Draw latent class memberships from their marginal distribution
draw_classes(n, pi_)draw_classes(n, pi_)
n |
Integer. Sample size. |
pi_ |
Numeric vector of length T. Class proportions (must sum to 1). |
Integer vector of length n with values in 1:T.
# Draw 100 class labels from equal prevalences draw_classes(100, c(1/3, 1/3, 1/3))# 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")
draw_classes_given_Zp(Zp, params)draw_classes_given_Zp(Zp, params)
Zp |
Numeric vector of length n. Covariate values. |
params |
Multinomial logistic parameter list (see |
Integer vector of length n with class labels in 1:T.
Zp <- draw_Zp(1000) X <- draw_classes_given_Zp(Zp, bk2018_params$covariate_params) table(X) # Should be roughly uniformZp <- 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
draw_indicators(X, rho)draw_indicators(X, rho)
X |
Integer vector of length n. True latent class (1-indexed). |
rho |
T x K matrix. |
An n x K integer matrix of 0/1 values.
rho <- make_rho(0.9) X <- draw_classes(50, c(1/3, 1/3, 1/3)) draw_indicators(X, rho)rho <- make_rho(0.9) X <- draw_classes(50, c(1/3, 1/3, 1/3)) draw_indicators(X, rho)
Draw a continuous distal outcome given true class memberships (scenario "distal")
draw_Zo(X, params)draw_Zo(X, params)
X |
Integer vector of length n. Latent class (1-indexed). |
params |
List with |
Numeric vector of length n.
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 muX <- 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}
draw_Zp(n)draw_Zp(n)
n |
Integer. Sample size. |
Integer vector of length n.
Zp <- draw_Zp(100) table(Zp)Zp <- draw_Zp(100) table(Zp)
Fixes mPhi at fit0$mPhi and estimates multinomial logit coefficients
mGamma (Q x (T-1)) via an EM algorithm with a BFGS M-step.
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 )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 )
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 |
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$mPhi throughout 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. TRUE if the EM loop exited before
maxIter iterations or if the final log-likelihood change was below
tol.
n_obsInteger. Number of observations used in estimation after listwise deletion on covariates.
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$convergedd <- 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
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.
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 )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 )
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. |
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_avg in
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 for mGamma) if
available.
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 vcovd <- 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
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.
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 )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 )
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 |
Nested list indexed as
datasets[[scenario]][[separation]][[as.character(n)]],
each element a list of n_rep data frames.
# 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 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
generate_data( n, separation = c("low", "mid", "high"), scenario = c("covariate", "distal"), params = bk2018_params, seed = NULL )generate_data( n, separation = c("low", "mid", "high"), scenario = c("covariate", "distal"), params = bk2018_params, seed = NULL )
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 |
A data.frame with columns:
Y1 .. Y6
Binary 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).
# 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)# 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)
Estimates the latent class measurement model with multilevLCA and
optionally, fixes mPhi and estimates covariate effects (two-step
initialization) with fitZ_from_fit0().
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 )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 )
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 |
A list with $fit0 (multilevLCA::multiLCA() measurement model) and $fitZ
(two-step covariate model from fitZ_from_fit0(), or NULL).
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 estimatesd <- 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
Returns a T x K matrix rho where rho[t, k] = P(Y_k = 1 | X = t).
make_rho(pi_)make_rho(pi_)
pi_ |
Numeric scalar in (0.5, 1). Probability of the "likely" response.
Use |
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).
A 3 x 6 numeric matrix.
# High separation: P(Y=1|class) = 0.9 for the "high" class make_rho(0.9) # Low separation make_rho(0.7)# High separation: P(Y=1|class) = 0.9 for the "high" class make_rho(0.9) # Low separation make_rho(0.7)
Evaluates P(X = t | Zp) for each observation using a multinomial logit with one or more covariates and class-specific intercepts and slopes.
mnl_probs(Zp, params)mnl_probs(Zp, params)
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 |
An n x T matrix of class probabilities (rows sum to 1).
# 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)# 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)
Delegates to plot.multiLCA from multilevLCA, which draws the
class-specific item-response probability profiles from the Step-1
measurement model.
## 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, ...)## 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, ...)
x |
A |
horiz |
Logical. If |
clab |
Optional character vector of length T giving class labels. |
... |
Further arguments passed to |
Called for its side effect (a base-graphics plot). Invisibly
returns NULL.
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"))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"))
Compact one-line or table summary printed to the console.
## 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, ...)## 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, ...)
x |
A |
... |
Further arguments. |
digits |
Integer. Number of decimal places for coefficient tables. |
Invisibly returns x.
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)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)
Verbose summary including model fit, class prevalences, item-response probabilities, and coefficient tables with standard errors and p-values.
## 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, ...)## 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, ...)
object |
A |
... |
Further arguments (currently unused). |
digits |
Integer. Number of decimal places for coefficient tables. |
Invisibly returns object.
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)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)
Fits a three-step latent class model through the following steps:
Measurement model: estimates latent class parameters
(, ) using multilevLCA
(Lyrvall et al., 2025).
Classification-error matrix: computes posterior class
probabilities and the T x T misclassification probability matrix
, 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.
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 )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 )
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 |
An S3 object of class tseLCA. The subclass depends on which
models were estimated:
tseLCA_measurementReturned when neither Zp.names nor
Zo.name is supplied. Contains the following elements:
measurement_modelStep-1 output list from lca_step1().
llikFinal Step-1 log-likelihood.
AIC, BIC
Information criteria from the measurement model.
R2entrEntropy R of 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.names is supplied and
Zo.name is NULL. Contains all elements of
tseLCA_measurement plus:
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
unless use.simple.cov = TRUE.
two_stepQ x (T-1) matrix of two-step starting values, or
NULL if use.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.
llikProfile log-likelihood
,
with Step-1 parameters held fixed. By construction
smaller than the equivalent one-step MLE likelihood.
tseLCA_distalReturned when Zo.name is supplied and
Zp.names is NULL. 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, named mu_C1 through mu_CT.
three_step.llikStep-3 distal log-likelihood
at converged estimates.
llikProfile log-likelihood
,
with Step-1 parameters 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.names and
Zo.name are supplied. Contains:
covariateA tseLCA_covariate-structured sub-list
(see above), including llik, AIC, BIC,
entropy.R2.
distalA 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.
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
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.
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)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)
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.
## 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"), ... )## 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"), ... )
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 |
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.
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")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")