---
title: "tseLCA Workflow"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{tseLCA Workflow}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5
)
```

```{r setup, message = FALSE}
library(tseLCA)
```

## Overview

Latent class analysis (LCA) groups observations into unobserved classes from
a set of categorical indicators. Researchers usually also want to know how the
classes relate to other variables: *covariates* that predict class membership,
and *distal outcomes* that the classes predict.

`tseLCA` does this with **three-step estimation**:

1. **Measurement model.** Estimate the latent classes from the indicators
   alone (`tse_lca()`).
2. **Classification.** Assign observations to classes and quantify the
   classification error (`tse_classify()`).
3. **Structural model.** Relate the classes to covariates (`tse_covariate()`)
   or distal outcomes (`tse_distal()`), correcting for the classification
   error of Step 2.

Because the measurement model is fixed before any structural variable enters,
the covariates and outcomes cannot change what the classes mean. That is the
main reason to prefer three-step over one-step estimation, in which
indicators, covariates, and outcomes are modeled jointly and the class
solution can shift with every change in the structural specification.
Assigning observations to classes and then analyzing the assignments, as if
they were the true classes, biases the structural estimates toward zero;
the bias-adjusted estimators of Bolck, Croon, and Hagenaars (2004; "BCH") and
Vermunt (2010; "ML") remove that bias, and `tseLCA` adds standard errors that
account for the uncertainty of the Step-1 measurement model (Bakk, Oberski,
and Vermunt 2014).

Each step returns an object that can be inspected with the usual R tools
(`print()`, `summary()`, `coef()`, `vcov()`, `confint()`, `logLik()`,
`AIC()`, `BIC()`, `predict()`, `plot()`) before moving to the next. The
one-call interface `tseLCA()` runs all steps at once.

## Example data

`generate_data()` simulates data from the design of Bakk and Kuha (2018):
three classes, six binary indicators `Y1`--`Y6`, and either a covariate `Zp`
(`scenario = "covariate"`) or a continuous distal outcome `Zo`
(`scenario = "distal"`). `separation` controls how well the indicators
separate the classes; the true class `X` is included for reference.

```{r data}
d <- generate_data(n = 1000, separation = "high", scenario = "covariate", seed = 1)
d$Zo <- draw_Zo(d$X, bk2018_params$distal_params) # add a distal outcome
head(d)
```

## Step 1: the measurement model

### Choosing the number of classes

The number of classes is chosen from the measurement model alone, before any
covariates or outcomes are considered. With several values of `nclass`,
`tse_lca()` returns a class-enumeration table of fit statistics (see Nylund,
Asparouhov, and Muthén 2007, and Masyn 2013, for guidance on using them).

```{r enumeration}
f_items <- cbind(Y1, Y2, Y3, Y4, Y5, Y6) ~ 1
sel <- tse_lca(f_items, data = d, nclass = 1:4)
sel
```

```{r enumeration-plot, fig.height = 4}
plot(sel)
```

The BIC and the sample-size adjusted BIC (SABIC) favor three classes. The
AIC, which penalizes additional parameters less, marginally prefers four, but
the fourth class holds less than 2% of the sample, a common sign of
over-extraction; the three-class model is also the more interpretable. In
practice the choice combines these criteria with class sizes, separation
(entropy), and substantive interpretation. `best_model()` extracts the model
selected by a criterion; `sel[[k]]` extracts the `k`-class model.

```{r best}
m <- best_model(sel, criterion = "BIC")
```

### Inspecting the measurement model

```{r measurement}
summary(m)
```

```{r measurement-plot}
plot(m)
```

`class_sizes()` and `item_probs()` give the parameters on the probability
scale; `coef()` and `vcov()` give them on the unconstrained log-ratio scale in
which they are estimated.

```{r measurement-accessors}
class_sizes(m)
item_probs(m)
```

Latent class models can have several local optima. By default `tse_lca()`
uses the k-means initialization of **multilevLCA**, with additional random
starts when the entropy is low; `control = tse_control(n_init = 20)` fits the
model from 20 random starts, and `start` fixes a starting
classification.

## Step 2: classification

`tse_classify()` assigns each observation to its most likely class (modal
assignment) or, with `assignment = "proportional"`, to every class with its
posterior probability as weight. It reports the classification-error
probabilities $P(W = s \mid X = t)$ between the true class $X$ and the
assigned class $W$, which the Step-3 estimators correct for.

```{r classify}
cl <- tse_classify(m)
cl
```

With these well-separated classes, 92--97% of the members of each class are
assigned to it. `posterior()` and `classes()` return the posterior
probabilities and the modal classes.

## Step 3: covariates

`tse_covariate()` estimates a multinomial logistic regression of class
membership on the covariates. The default estimator is ML with standard
errors corrected for the Step-1 uncertainty.

```{r covariate}
fc <- tse_covariate(cl, ~ Zp)
summary(fc)
```

The coefficients are named `covariate:class` and the first class is the
reference. The generic tools work as usual:

```{r covariate-tools}
confint(fc)
AIC(fc)
anova(fc) # Wald test of each covariate term, across all classes
```

`predict()` gives the class-membership probabilities at given covariate
values:

```{r covariate-predict}
predict(fc, newdata = data.frame(Zp = 1:5))
```

### Estimators and standard errors

The BCH estimator and the uncorrected three-step estimator (which analyzes
the assigned classes as if they were the true classes) are available for
comparison; so is the two-step estimator of Bakk and Kuha (2018).

```{r estimators}
fc_bch <- tse_covariate(cl, ~ Zp, method = "BCH")
fc_raw <- tse_covariate(cl, ~ Zp, method = "none")
ft <- tse_twostep(m, ~ Zp)
round(cbind(
  ML = coef(fc), BCH = coef(fc_bch), uncorrected = coef(fc_raw), two.step = coef(ft)
), 3)
```

The uncorrected estimates are attenuated toward zero (the true slopes are
$-1$ and $1$); the bias-adjusted estimates are not.

`se = "robust"` omits the Step-1 correction. With well-separated classes the
two are close; the correction matters more when classes are less distinct.

```{r se}
round(cbind(
  corrected = sqrt(diag(vcov(fc))),
  robust = sqrt(diag(vcov(tse_covariate(cl, ~ Zp, se = "robust"))))
), 4)
```

Under low separation, proportional assignment (`tse_classify(m, assignment =
"proportional")`) with the ML estimator is generally the most reliable choice.

### Reference class and covariate formulas

`ref` (or `relevel()` on a fitted model) changes the reference class.
Covariates follow the usual formula syntax, including factors, interactions,
and transformations. A variable created after classification is supplied with
`data`, which must hold the classified rows (it may add columns).

```{r ref}
coef(relevel(fc, ref = "C3"), matrix = TRUE)

d$group <- factor(ifelse(d$Zp > 3, "high", "low"))
anova(tse_covariate(cl, ~ Zp + group, data = d))
```

## Step 3: distal outcomes

`tse_distal()` estimates the distribution of a distal outcome in each class.
Outcomes can be `"gaussian"` (class means with a common variance),
`"poisson"`, `"binomial"`, or `"multinomial"` (nominal).

```{r distal}
fd <- tse_distal(cl, Zo ~ 1)
summary(fd)
```

`omnibus_test()` tests whether the outcome's distribution differs across
classes:

```{r omnibus}
omnibus_test(fd)
```

For a nominal outcome, `coef(fit, matrix = TRUE)` gives the class-by-category
probability matrix:

```{r multinomial}
d$Zcat <- cut(d$Zo, c(-Inf, -0.5, 0.5, Inf), labels = c("low", "mid", "high"))
fm <- tse_distal(cl, Zcat ~ 1, family = "multinomial", data = d)
round(coef(fm, matrix = TRUE), 3)
omnibus_test(fm)
```

### Covariates and a distal outcome

Passing a covariate model to `tse_distal()` fits both parts. The class prior
then depends on the covariates, and the uncertainty of the covariate model is
propagated to the distal estimates.

```{r combined}
fb <- tse_distal(fc, Zo ~ 1)
fb
```

## All steps in one call

`tseLCA()` runs the three steps from one formula,
`indicators ~ covariates | distal outcome`. The components remain available
through `measurement()`, `classification()`, `covariate()`, and `distal()`.

```{r one-call}
fit <- tseLCA(cbind(Y1, Y2, Y3, Y4, Y5, Y6) ~ Zp | Zo, data = d, nclass = 3)
all.equal(coef(fit), coef(fb))
round(classification(fit)$D, 3)
```

## A measurement model from another sample

Because the measurement model is fixed in Step 1, it can be estimated on one
sample and applied to another, for example when covariates are observed only
in a subsample. `tse_classify(m, newdata = ...)` classifies the new sample
with the existing measurement model; the Step-1 uncertainty in later steps is
that of the sample the model was estimated on.

```{r multisample}
sub <- d[1:300, ]
fc_sub <- tse_covariate(tse_classify(m, newdata = sub), ~ Zp)
coef(fc_sub)
nobs(fc_sub)
```

## Missing data and indicator coding

Indicators can be factors, logicals, character variables, or numeric codes in
any coding; their categories are stored with the model and reused on new
data. With `missing = "fiml"`, observations with some missing indicators are
kept (full-information maximum likelihood); the default drops them. Rows with
a missing covariate or distal outcome are dropped from that Step-3 model only.

```{r missing}
d_miss <- d
set.seed(2)
d_miss$Y1[sample(nrow(d), 100)] <- NA
d_miss$Y2 <- factor(d_miss$Y2, labels = c("no", "yes"))
m_fiml <- tse_lca(f_items, data = d_miss, nclass = 3, missing = "fiml")
nobs(m_fiml)
nobs(tse_lca(f_items, data = d_miss, nclass = 3))
```

## Estimation settings

`tse_control()` collects the numerical settings: iteration limits and
tolerances for Step 1 and Step 3, random starts, the boundary tolerance for
the Step-1 variance, and the information matrix used for Step-3 standard
errors.

```{r control}
tse_control(step1.maxit = 10000, n_init = 20)
```

## Migrating from tseLCA 1.x

`three_step()` still works (with the same estimates) but is deprecated. The
table maps its arguments to the new interface.

| `three_step()` | tseLCA 2.0 |
|---|---|
| `Y.names`, `n_classes` | `tse_lca(cbind(...) ~ 1, nclass = )` |
| `Zp.names` | `tse_covariate(, ~ ...)` |
| `Zo.name`, `family` | `tse_distal(, outcome ~ 1, family = )` |
| `step1` (measurement model from another sample) | `tse_classify(, newdata = )` |
| `startval`, `n_init` | `tse_lca(start = )`, `tse_control(n_init = )` |
| `use.modal.assignment` | `tse_classify(assignment = )` |
| `use.bch`, `use.simple.cov` | `method = "BCH"`, `se = "robust"` |
| `rebase` | `ref`, or `relevel()` |
| `incomplete` | `tse_lca(missing = "fiml")` |
| other tuning arguments | `tse_control()` |
| `get.twostep.vcov` | `tse_twostep(se = TRUE)` |

`coef()` now returns a named vector matching `vcov()` (use
`coef(fit, matrix = TRUE)` for the coefficient matrix), so `confint()` works.
See `NEWS.md` for all changes, including bug fixes affecting results for
polytomous indicators, Gaussian distal outcomes, and proportional-assignment
ML distal models.

## References

Bakk, Z., & Kuha, J. (2018). Two-step estimation of models between latent
classes and external variables. *Psychometrika*, 83(4), 871--892.

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.

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.

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.

Masyn, K. E. (2013). Latent class analysis and finite mixture modeling. In
T. D. Little (Ed.), *The Oxford Handbook of Quantitative Methods*, Vol. 2,
551--611. Oxford University Press.

Nylund, K. L., Asparouhov, T., & Muthén, B. O. (2007). Deciding on the number
of classes in latent class analysis and growth mixture modeling: A Monte Carlo
simulation study. *Structural Equation Modeling*, 14(4), 535--569.

Vermunt, J. K. (2010). Latent class modeling with covariates: Two improved
three-step approaches. *Political Analysis*, 18(4), 450--469.
