H^-1 M H^-1. The right-hand factor should be the
transpose of H^-1, and for the multinomial ML estimator H (the Jacobian of
the estimating equation) is not symmetric. The variance matrix was not
symmetric, and the standard errors were wrong, both robust and corrected, and
under both modal and proportional assignment. For a two-category outcome, the
standard errors of the two categories of a class now agree, and they match
those of the equivalent Binomial model. The point estimates, the BCH
estimator (whose H is diagonal), and the other distal families (whose H
is symmetric) are unchanged.step3.maxit iterations, whatever the
class separation. It now checks the negative Hessian. Fits that do not
converge now stop with "BCH NR reached maximum iterations" instead.tse_distal() of a
tse_covariate() fit), the respondents were reassigned to classes with the
covariate-adjusted posteriors, and the classification-error matrix was
computed from those assignments. The assignments then depended on the
covariates given the class, which biased the class means of the distal
outcome (about -0.08 for a mean of 1 at n = 1000, mid separation). The model
now keeps the Step-2 assignments and classification errors, which use the
indicators only, and uses the covariate-adjusted class probabilities as the
prior, as the corrected standard errors assume.tse_classify() now reports the entropy R^2 as the measurement model does
(relative to the entropy of the class proportions, as in Latent GOLD),
instead of relative to log(T). as_tse_lca() models use the same definition.as_tse_lca() has a new argument vcov: the variance matrix of the Step-1
parameters, in the layout of vcov() of a tse_lca() model. Corrected
standard errors in Step 3 now use it. Before, the Step-1 variance of such a
model was recomputed at the given parameters on data, which is only valid
when they are the maximum likelihood estimates for data. For parameters
estimated on another sample, the correction was wrong: too large when data
is smaller than that sample, too small when it is larger.vcov, Step 3 (tse_covariate(), tse_distal()) of an
as_tse_lca() model gives a warning and reports robust standard errors
(se = "robust"), which do not include the Step-1 uncertainty. vcov() of
such a model is an error, since it has no variance matrix.newdata argument of tse_classify() on the original tse_lca() model
keeps the original Step-1 variance.tse_lca() fits the Step-1 measurement model from a formula,
cbind(Y1, Y2, ...) ~ 1, with a control = tse_control() argument.
With several numbers of classes (e.g. nclass = 1:6) it returns a
tseLCA_select class-enumeration table: log-likelihood, parameters, AIC,
BIC, SABIC, entropy R^2, and smallest class. The table has print(),
plot(), as.data.frame(), best_model(), and [[ methods. The
one-class (independence) model is fitted in closed form.tse_classify() (Step 2) assigns observations to classes with a
fixed measurement model (assignment = "modal" or "proportional"). It
reports the classification-error probabilities P(W = s | X = t) and
entropy R^2. With newdata it classifies another sample, replacing the
step1 = route for using a measurement model from one sample on another.
tse_lca() models now keep their data for this purpose.tse_covariate() fits the multinomial logit of class membership on a
covariate formula (factors, interactions, transformations). Arguments:
method = "ML", "BCH", or "none" (the uncorrected three-step
estimator, for comparison); se = "corrected" or "robust"; ref
(reference class); start. Methods: predict() (class probabilities
given covariates), relevel(), and anova() (Wald tests per term).tse_distal() fits a distal outcome (Zo ~ 1, any supported
family, as a string or family object). Given a tse_covariate()
model it fits the combined model.tse_twostep() gives two-step estimates (Bakk & Kuha 2018). With
se = TRUE they come with multilevLCA's corrected standard errors.tseLCA() fits a whole model in one call:
cbind(indicators) ~ covariates | distal outcome (Formula package). It
chains the step-wise functions. Accessors measurement(),
classification(), covariate(), and distal() extract the components
of any fitted model.tse_covariate() and tse_distal() accept data: the classified data
(same rows) with additional columns, e.g. variables created after
classification.omnibus_test() returns a standard htest object.three_step() arguments to the new functions.predict() for measurement models gives posterior class probabilities
(or modal classes) for new data, coded with the model's stored
categories. fitted(), formula(), and update() also work.item_probs(se = TRUE) and class_sizes(se = TRUE) also return standard
errors (delta method from the measurement model's variance).start in tse_lca() and tseLCA() accepts item_probs() of a fitted
model directly, e.g. to refit a chosen model with covariates.as_tse_lca() builds a measurement model from given class sizes and
item-response probabilities (e.g. estimates reported elsewhere, or saved
from an earlier fit) without re-estimating it. The result can be passed to
tse_classify() like any tse_lca() model.tseLCA is the parent class of
every fit; tseLCA_structural is the parent of tseLCA_covariate,
tseLCA_distal, and tseLCA_both. Methods are defined once on the parent
classes, not separately for each subclass.logLik() (with df and nobs, so AIC()
and BIC() work), nobs(), posterior(), classes(),
class_sizes(), and item_probs().summary() now returns a summary.tseLCA_structural or
summary.tseLCA_measurement object. Its tables print with
printCoefmat() and are extracted with coef(summary(fit)).plot() now works on distal-outcome fits. These objects previously did
not store the measurement model, so the plot failed.$measurement_model$Y.levels). Previously they had to be integers
starting at 0. Codes such as 1..K were accepted but silently misread: the
highest category was never matched.step1 applies its own categories to
new data. The new data may omit categories, and values outside them are
an error.model.frame()/model.matrix(), so
factor covariates are dummy coded and unused factor levels are dropped.
The intercept is named "(Intercept)" (was "Intercept").tse_control() collects the numerical estimation settings.tse_lca(missing = "fiml") with random starts (n_init) failed with a
cryptic error from multilevLCA ("sort_index(): detected NaN") when
indicator values were missing: multilevLCA cannot start from a given
classification on incomplete data. It now warns and uses multilevLCA's
default initialization; a user-supplied start gives a clear error.
Corrected (ML) standard errors: the Jacobian of the classification-error matrix with respect to the Step-1 parameters, which carries the Step-1 uncertainty into Step 3, had its item-parameter columns ordered item by item, while the Step-1 variance is ordered class by class; with modal assignment it also differentiated the assignment as if it were the posterior probabilities. Both are fixed (checked against numerical derivatives). Point estimates, robust, and BCH standard errors are unchanged; corrected covariate standard errors typically increase by a few percent, and corrected standard errors no longer depend on the choice of reference class.
Structural models with a reference class other than the first
(ref =): posterior(), classes(), item_probs(), class_sizes(), and
plot() reported the classes in the rebased order under the original
labels, and in combined covariate + distal models the distal parameters
of class t belonged to another class. All now use the measurement
model's class order.
Step-1 variance of polytomous indicators when the first (reference)
category of an item has a boundary probability in some class: the free
parameters of that item are log-ratios against this category, and their
common shift, log P(Y = first | class), is not informed by the data. It is
now treated as fixed, like other boundary parameters. Previously the
variance of these log-ratios was huge (the information matrix was nearly
singular) or all NA (singular), which made the corrected standard errors
of the Step-3 models NA.
If the Step-1 information matrix is still singular, Step-3 models warn
and report robust standard errors (was NA).
tse_twostep(se = TRUE) with a reference class other than the first
returned the class-1-reference estimates and variance under the new class
labels. They are now transformed to the requested reference class. It also
no longer warns, wrongly, that multilevLCA's measurement model differs.
For measurement-only fits, posterior() / $posteriors and classes() /
$classifications were not in data-row order. They were read from
multilevLCA's fit0$mU, which is sorted by response pattern. They are now
computed from the Step-1 data in data-row order. Code that used
$classifications from a measurement model, e.g. as startval, received
mismatched rows in 1.1.x.
Polytomous indicators were decoded from fit0$mU as category 0 whenever
the measurement model was fitted with listwise deletion
(incomplete = FALSE). The Step-1 information matrix was then singular, so
vcov() of a polytomous measurement model and the Step-1-corrected Step-3
variance (use.simple.cov = FALSE) were all NA. The measurement
model now stores its Step-1 sample in data-row order and uses it for all
of these computations. The fallback decoder for objects without stored data
handles both of multilevLCA's codings.
Parameter counts behind the reported AIC/BIC were wrong. Covariate models
counted one-hot indicator columns, not free item parameters (e.g.
40, not 22, for six binary items, three classes, one covariate).
Combined covariate + distal models counted class sizes, not the
covariate coefficients. Reported AIC/BIC for these models change;
estimates and standard errors do not.
Gaussian distal outcomes: the within-class variance was fixed at
sigma2 = 1 in the Step-3 likelihood. Three-step ML now estimates a common
sigma2 jointly with the class means, as in Bakk, Tekle & Vermunt (2013).
Unlike ordinary regression, sigma2 does not factor out of the mean
estimates here: it enters the posterior weights P(X = t | W, Zo). With the
variance fixed at 1, ML class means were biased whenever the outcome's
variance differed from 1 (e.g. roughly 3.5x too far apart for a residual SD
of 3) and were not equivariant to rescaling the outcome. The Hessian,
score, and Step-1/Step-2 uncertainty propagation now include sigma2. It is
reported as $sigma2 (estimate and, under ML, standard error), and counted
in the log-likelihood's degrees of freedom. BCH class means were
unaffected; their log-likelihood now uses the estimated variance.
Three-step ML for distal outcomes with proportional assignment
(use.modal.assignment = FALSE) maximized the wrong likelihood. It used
log sum_t P(t) f(z|t) sum_s w_s P(W=s|t), with the assignment weights inside
the log. The likelihood of Bakk, Tekle & Vermunt (2013), for the expanded
data file weighted by the posterior assignment probabilities, is
sum_s w_s log sum_t P(t) f(z|t) P(W=s|t). The covariate model already used
this form. The E-step, score, Hessian/Jacobian, and Step-1/Step-2
uncertainty propagation now use the correct likelihood for all distal
families. Proportional-assignment ML distal estimates were biased away from
zero (e.g. class means of -1.07 and 1.06, not -1 and 1, at mid
separation). Modal-assignment fits, whose single assignment makes the two
forms identical, BCH, and covariate models are unchanged.
Covariate models with include.intercept = FALSE failed with a
"length of 'dimnames'" error. The coefficient rows were always labeled
with an Intercept row.
In combined covariate + distal models, rows with a missing covariate but
an observed distal outcome were kept in the distal model with NA class
priors. The distal model now uses rows with complete covariates.
When the Step-3 Hessian of an ML covariate model could not be inverted
without an error being raised (non-finite entries), standard errors fell
back to the outer product of the scores without a warning. Every fallback
now warns. The documentation of hessian = "opg" (correct.spec in
three_step()) described the outer-product estimator as model-robust; it
is valid only under a correctly specified Step-3 model.
three_step() is deprecated in favor of tseLCA() and the step-wise
functions. It keeps working, with the same estimates, and warns once per
session. Its help page maps each argument to the 2.0 interface.lca_step1(), lca_step1_startval(), fitZ_from_fit0(), and
fitZ_from_multiLCA() are deprecated in favor of tse_lca() and
tse_twostep(). They warn once per session when called directly.options(tseLCA.warn.deprecated = FALSE) silences these warnings.coef() returns a named vector whose names match vcov(), so
confint() works. coef(fit, matrix = TRUE) gives the previous
Q x (T-1) (covariate) or T x C (multinomial distal) layout.coef() on a measurement model returns the log-ratio parameters that
vcov() describes. The previous list is now available as
class_sizes() and item_probs()."(Intercept)" (was "Intercept"),
e.g. "(Intercept):C2".which argument of coef() and vcov() is replaced by
component ("covariate" or "distal") and step ("two_step").
Passing which gives an error, so old code cannot silently return a
different quantity.vcov() on a tseLCA_both object returns one matrix. The covariate
and distal blocks are on the diagonal; the cross-covariances are not
computed and are NA.Authors@R from contributor ("ctb") to
author and copyright holder (c("aut", "cph")), reflecting his contribution
to the package. No code changes.lca_step1_startval(), a wrapper around multilevLCA::multiLCA()
that injects a user-supplied Step-1 starting value and sets kmea = FALSE,
bypassing multilevLCA's default deterministic k-means-on-principal-
components initialization. This addresses feedback that the default
initialization can consistently converge to a local optimum of the Step-1
log-likelihood on some datasets. startval accepts either:
1..n_classes, one entry per row of
data), e.g. the modal class from an external solver run with many
random starts ('StepMix', 'poLCA'); orP(Y_h = k | X = t), from which a classification is derived internally
(naive-Bayes argmax under a flat class prior). This is the natural
format for an externally estimated Step-1 solution that isn't tied to
the current sample, e.g. poLCA's probs output or a published
item-response table.
Either way, users can pass their external Step-1 solution straight through
instead of hand-assembling a tseLCA measurement object.startval argument to lca_step1() and
fitZ_from_multiLCA(), so every place a measurement model is fit from raw
data can use an external starting value instead of multilevLCA's k-means
initialization. Supplying startval skips the
iter.measurement/R2.threshold random-restart logic, since restarting
from fresh k-means seeds would defeat the purpose of a user-vetted start.startval argument to three_step(), forwarded to lca_step1()
(and to fitZ_from_multiLCA() when get.twostep.vcov = TRUE) so the
measurement model can be fit from an external starting value in a single
call.n_init argument to lca_step1(), fitZ_from_multiLCA(), and
three_step(): the unconditional multi-random-start analog of n_init in
StepMix or nrep in poLCA. When supplied, the measurement model is
fit n_init times from independent uniform-random classifications
(kmea = FALSE, not multilevLCA's deterministic k-means-on-PCA path),
and the highest-log-likelihood fit is kept. This is a separate mechanism
from the existing iter.measurement/R2.threshold restart logic, which
reruns multilevLCA's own k-means initialization and only when entropy
R^2 is low; n_init restarts always run and never use k-means.step1, startval, and n_init are mutually exclusive ways of
controlling Step 1 in three_step() (and startval/n_init in
lca_step1() and fitZ_from_multiLCA()); supplying more than one errors.use.bch = TRUE estimation (both covariate and distal-outcome models).
pwx[s, t] = P(W = s | X = t) is column-stochastic; the correct BCH
weight matrix is w.is %*% t(pwx)^-1 (Mplus Web Note 21: the row of
H^-1 for each case's most likely class, where H = t(pwx) is
row-stochastic), not w.is %*% pwx^-1, which the code had been computing.
The two orientations only agree when pwx is symmetric, so the bug was
largely invisible in well-separated, balanced test cases; with genuine
classification-error asymmetry it biased BCH point estimates and could
produce the "negative column sums" error the package warns about (users
were advised to fall back to use.bch = FALSE). With the corrected
orientation, every row of the weight matrix sums to 1 and the class
totals it implies exactly match the posterior class sizes. The
weight-matrix computation is now consolidated into a single internal
helper (bch_weight_matrix()) used by all four call sites that
previously duplicated it, with a regression test asserting
rowSums(w.it) == 1 under a deliberately asymmetric classification-error
matrix.family = "multinomial" to three_step()'s distal-outcome
estimation, for a nominal categorical outcome with 2 or more categories
(Zo.name may be a factor, character, or integer column). Estimates a
saturated model -- the T x C matrix of class-conditional category
probabilities pi_hat[t, c] = P(Zo = c | X = t) -- with the closed-form
weighted-proportion estimator pi_hat[t, c] = sum_i w_it * 1(y_i = c) / sum_i w_it, for both use.bch = TRUE (BCH weights) and use.bch = FALSE (ML, through EM with the same closed-form M-step and
responsibility-weighted E-step). This replaces the previous workaround of
fitting one family = "binomial" model per category and renormalizing
the resulting probabilities by hand, which wasn't constrained to the
simplex before renormalizing and had no joint covariance across
categories.
coef() returns the T x C probability matrix (rows sum to 1) instead
of a length-T vector; vcov() returns its (T*C) x (T*C) sandwich
covariance, which is necessarily rank-deficient (each class's row sums
to 1).use.bch = FALSE, use.simple.cov = FALSE (the default) fully
propagates both Step-1 measurement uncertainty and, when Zp.names is
also supplied, Step-3 covariate uncertainty into the SEs, matching the
existing gaussian/poisson/binomial ML paths exactly (the T x C
generalization of each chain-rule term -- C1_mat for Step 1, C_mat
for Step 2 -- only required expanding the "unit score" g_it across
categories, since neither term's derivation otherwise depends on the
outcome's dimensionality). The ML bread (multinomial_ml_jacobian())
is a closed-form Jacobian of the estimating equation, generalizing
ml_hessian_distal()'s "observed = complete - missing information"
correction to the T x C case. Both the bread and the full propagation
(including the Step-2 term) were cross-validated against independent,
from-scratch numerical differentiation (explicit per-case loops sharing
no code with the package, optim() from multiple starting points for
the point estimate) and matched to machine precision.omnibus_test(), a generalized Wald test of
H0: theta_1 = theta_2 = ... = theta_T (the distal outcome's class-t
parameter vector is the same for every class) for any tseLCA_distal or
tseLCA_both object, regardless of family. Uses a Moore-Penrose
pseudo-inverse of the contrast covariance so it remains valid when that
covariance is singular, as it always is for family = "multinomial"; the
resulting degrees of freedom recover the textbook (T-1)*(C-1) for a
T x C chi-squared test of homogeneity in that case, and T-1 for the
scalar-parameter families.family = "binomial", whose coef()/vcov() are on the logit
scale, family = "multinomial" reports coef()/vcov() directly on the
probability scale -- Std.Error is directly interpretable without a
delta-method back-transform, but a symmetric interval
Estimate +/- 1.96*Std.Error can fall outside the unit interval for a
probability near a boundary, the same known limitation as a naive Wald
interval for a sample proportion, and the per-cell z.value/p.value
(testing each probability against 0) are rarely the question of interest.
print()/summary() now print a one-line reminder of this after a
multinomial distal-outcome table, pointing to omnibus_test() for the
intended, boundary-safe test of whether the distribution differs across
classes.T/F with TRUE/FALSE throughout internal codebase\value tags to all exported functions missing them, including bk2018_params.inst/examples: examples now write to tempdir() instead of the home
filespace.inst/examples: commented out rm(list = ls()) calls.inst/examples: commented out install.packages() calls.use.simple.cov = TRUE) that bypasses the measurement-uncertainty correction for faster computation in large, well-separated samples.incomplete = TRUE argument (using a two-pass row-filtering strategy).step1 argument) to reuse across multiple structural models or apply to different sample subsets.use.modal.assignment).rebase argument to allow users to easily change the reference latent class for the multinomial logit parameterization while maintaining invariant log-likelihoods.fitZ_from_fit0()) to generate stable starting values for the three-step structural model.tseLCA objects: summary(), coef(), vcov(), and plot() (which delegates to 'multilevLCA' for item-profile visualization).generate_data()) that replicates the Bakk & Kuha (2018) simulation study design for both covariates and distal outcomes under varying separation conditions.