library(tidyverse) # dplyr, ggplot2, purrr, tibble, readr, stringr, forcats
source("_common.R") # book-wide helpers: round2(), fmt_p(), tidy2()
set.seed(2026)
N <- 800; J <- 8
prev <- c(0.35, 0.65)
rho <- rbind(rep(0.80, J), rep(0.20, J))
rho[1, c(7, 8)] <- 0.30; rho[2, c(7, 8)] <- 0.70 # reverse two items: pattern, not severity
cls <- sample(1:2, N, replace = TRUE, prob = prev)
# one response vector per person; the estimator below wants them as a matrix
X <- do.call(rbind, map(cls, \(k) rbinom(J, 1, rho[k, ])))22 Latent Classes: Which Kind of Person Is This?
Factor analysis, in the last chapter, asked a dimensional question: where does each person fall on a continuum of depression, or curiosity? Latent class analysis (LCA) asks a categorical one: which kind of person is this? Instead of a score along an axis, it looks for hidden types - subgroups that respond in qualitatively different ways. A dimension sorts people high to low; a class says they come in distinct kinds.
This matters because some of the most important findings in psychology are subgroups that a single score would blur together. A total symptom count cannot tell you that some anxious people are also risk-prone approachers rather than shy avoiders - but a latent-class model can, and the authors’ own work found exactly that (more below).
22.1 Learning Objectives
- Contrast a dimensional latent variable (a factor) with a categorical one (a class).
- State the latent-class model as a mixture of response patterns.
- Fit one from scratch with the EM algorithm, and read the classes off the item probabilities.
- Recognize that classes can differ in pattern, not just severity - which is the whole point.
22.2 The Model: A Mixture of Kinds
Suppose there are \(K\) hidden classes. Each class \(k\) has (a) a size \(\pi_k\) (its share of the population) and (b) its own set of item-endorsement probabilities \(\rho_{jk}\) - the chance a member of class \(k\) says “yes” to item \(j\). A person belongs to one class, but we do not see which; we see only their answers. The probability of a whole response pattern is a mixture - a weighted blend across the classes:
\[P(\mathbf{x}_i) = \sum_{k=1}^{K} \pi_k \prod_{j=1}^{J} \rho_{jk}^{x_{ij}} (1-\rho_{jk})^{1-x_{ij}}\]
The classes are the profiles \(\rho_{\cdot k}\); membership is hidden. Our job is to recover both from the responses alone.
22.3 Fitting It From Scratch with EM
The difficulty is a chicken-and-egg: if we knew who was in each class we could estimate the class profiles, and if we knew the profiles we could assign people to classes - but we know neither. The EM algorithm breaks the loop by alternating, exactly as Rasch estimation did:
- E-step (Expectation): given current profiles, compute each person’s posterior probability of belonging to each class.
- M-step (Maximization): given those soft memberships, re-estimate each class’s size and item probabilities as membership-weighted averages.
Repeat until nothing moves. First, simulate two classes that differ in pattern - items 7–8 are reversed, so the classes are not simply “more” vs. “less” symptomatic:
NoteWorking in SPSS, Julia, or Python?
The code tabs below assume this chapter’s data is already loaded. Grab the one-file setup for your language from Getting the Book’s Data, run it once, then load what you need by name - this chapter uses ch22-lca. For example, book_data("ch22-lca") in R, Julia, or Python, or !bookdata name = "ch22-lca". in SPSS. Every language reads the same shipped files, so your numbers will match the ones printed here exactly.
* The book's two-class data: 800 people, 8 binary items, with items 7-8
* reversed so the classes differ in PATTERN rather than severity.
INSERT FILE='data/sim/ch22-lca.sps'.
DESCRIPTIVES VARIABLES=x1 TO x8 /STATISTICS=MEAN. Descriptive Statistics
+--------------------+---+----+
| | N |Mean|
+--------------------+---+----+
|x1 |800| .41|
|x2 |800| .39|
|x3 |800| .42|
|x4 |800| .41|
|x5 |800| .41|
|x6 |800| .43|
|x7 |800| .55|
|x8 |800| .57|
|Valid N (listwise) |800| |
|Missing N (listwise)| 0| |
+--------------------+---+----+using CSV, DataFrames
X = Matrix(CSV.read("data/sim/ch22-lca.csv", DataFrame))
size(X) # 800 people by 8 binary items(800, 8)import pandas as pd
X = pd.read_csv("data/sim/ch22-lca.csv").values
X.shape # 800 people by 8 binary items(800, 8)Now the entire estimator - about a dozen lines:
lca_em <- function(X, K = 2, tol = 1e-6, maxit = 500) {
N <- nrow(X); J <- ncol(X)
set.seed(1); pi <- rep(1/K, K); p <- matrix(runif(K * J, .2, .8), K, J)
ll_old <- -Inf
# EM is inherently sequential - each sweep starts from the previous estimates -
# so this stays a loop. The per-class E-step inside it does not, and maps.
for (it in 1:maxit) {
# E-step: log P(x_i, class k) for every person and class, then normalize to a posterior
loglik <- map(1:K, \(k)
log(pi[k]) + X %*% log(p[k, ]) + (1 - X) %*% log(1 - p[k, ])) |>
reduce(cbind)
m <- apply(loglik, 1, max)
post <- exp(loglik - m); post <- post / rowSums(post)
# M-step: re-estimate sizes and item probabilities as posterior-weighted averages
pi <- colMeans(post)
p <- t(post) %*% X / colSums(post)
ll <- sum(m + log(rowSums(exp(loglik - m))))
if (abs(ll - ll_old) < tol) break; ll_old <- ll
}
list(prevalence = pi, item_prob = p, iterations = it)
}
fit <- lca_em(X, K = 2)
cat("converged in", fit$iterations, "iterations\n")converged in 22 iterations# recovered class sizes (true: .35 / .65)
tibble(class = seq_along(fit$prevalence), prevalence = fit$prevalence) |> round2()| class | prevalence |
|---|---|
| 1 | 0.64 |
| 2 | 0.36 |
* Base SPSS has no latent class analysis, and neither does PSPP. A genuine
* LCA needs Latent GOLD, Mplus, or R's poLCA; SPSS's own nearest relative is
* TwoStep Cluster, which groups cases but estimates no class-conditional item
* probabilities and no membership posteriors.
*
* What SPSS CAN show is the raw endorsement rate of each item, which is the
* marginal the mixture model decomposes:
INSERT FILE='data/sim/ch22-lca.sps'.
DESCRIPTIVES VARIABLES=x1 TO x8 /STATISTICS=MEAN. Descriptive Statistics
+--------------------+---+----+
| | N |Mean|
+--------------------+---+----+
|x1 |800| .41|
|x2 |800| .39|
|x3 |800| .42|
|x4 |800| .41|
|x5 |800| .41|
|x6 |800| .43|
|x7 |800| .55|
|x8 |800| .57|
|Valid N (listwise) |800| |
|Missing N (listwise)| 0| |
+--------------------+---+----+using CSV, DataFrames, Statistics
X = Matrix(CSV.read("data/sim/ch22-lca.csv", DataFrame))
function lca_em(X, K = 2; tol = 1e-6, maxit = 500)
N, J = size(X)
pi_ = fill(1/K, K)
p = reshape(range(0.3, 0.7, length = K*J), K, J)
ll_old = -Inf; res = nothing
for it in 1:maxit
loglik = [log(pi_[k]) + sum(X[i,j]*log(p[k,j]) + (1-X[i,j])*log(1-p[k,j])
for j in 1:J) for i in 1:N, k in 1:K]
m = maximum(loglik, dims = 2)
post = exp.(loglik .- m); post ./= sum(post, dims = 2)
pi_ = vec(mean(post, dims = 1))
p = (post' * X) ./ sum(post, dims = 1)'
ll = sum(m .+ log.(sum(exp.(loglik .- m), dims = 2)))
res = (prevalence = pi_, item_prob = p, iterations = it)
abs(ll - ll_old) < tol && break
ll_old = ll
end
res
end
fit = lca_em(X, 2)
fit.prevalence # recovered class sizes (true: .35 / .65)2-element Vector{Float64}:
0.6434586522896812
0.35654134771031876import numpy as np, pandas as pd
X = pd.read_csv("data/sim/ch22-lca.csv").values
def lca_em(X, K=2, tol=1e-6, maxit=500):
N, J = X.shape
pi_ = np.full(K, 1 / K)
p = np.linspace(0.3, 0.7, K * J).reshape(K, J)
ll_old = -np.inf
for it in range(1, maxit + 1):
loglik = np.log(pi_) + X @ np.log(p).T + (1 - X) @ np.log(1 - p).T
m = loglik.max(axis=1, keepdims=True)
post = np.exp(loglik - m); post /= post.sum(axis=1, keepdims=True)
pi_ = post.mean(axis=0)
p = (post.T @ X) / post.sum(axis=0)[:, None]
ll = float((m + np.log(np.exp(loglik - m).sum(axis=1, keepdims=True))).sum())
if abs(ll - ll_old) < tol:
break
ll_old = ll
return dict(prevalence=pi_, item_prob=p, iterations=it)
fit = lca_em(X, 2)
fit["prevalence"] # recovered class sizes (true: .35 / .65)array([0.64345765, 0.35654235])The class sizes come back essentially exact. And the item-probability profiles - the definitions of the two kinds - recover the reversal we planted:
# rows = classes, columns = items 1..8
fit$item_prob |>
as_tibble(.name_repair = \(x) paste0("item", seq_along(x))) |>
mutate(class = row_number(), .before = 1) |>
round2()| class | item1 | item2 | item3 | item4 | item5 | item6 | item7 | item8 |
|---|---|---|---|---|---|---|---|---|
| 1 | 0.19 | 0.18 | 0.20 | 0.20 | 0.21 | 0.21 | 0.72 | 0.68 |
| 2 | 0.81 | 0.78 | 0.81 | 0.78 | 0.78 | 0.83 | 0.24 | 0.36 |
* The item-probability profiles - the DEFINITIONS of the two kinds - are
* exactly what an LCA produces and SPSS cannot. See the note above; the R,
* Julia and Python tabs estimate them from scratch with EM.using CSV, DataFrames, Statistics
X = Matrix(CSV.read("data/sim/ch22-lca.csv", DataFrame))
function lca_em(X, K = 2; tol = 1e-6, maxit = 500)
N, J = size(X)
pi_ = fill(1/K, K)
p = reshape(range(0.3, 0.7, length = K*J), K, J)
ll_old = -Inf; res = nothing
for it in 1:maxit
loglik = [log(pi_[k]) + sum(X[i,j]*log(p[k,j]) + (1-X[i,j])*log(1-p[k,j])
for j in 1:J) for i in 1:N, k in 1:K]
m = maximum(loglik, dims = 2)
post = exp.(loglik .- m); post ./= sum(post, dims = 2)
pi_ = vec(mean(post, dims = 1))
p = (post' * X) ./ sum(post, dims = 1)'
ll = sum(m .+ log.(sum(exp.(loglik .- m), dims = 2)))
res = (prevalence = pi_, item_prob = p, iterations = it)
abs(ll - ll_old) < tol && break
ll_old = ll
end
res
end
fit = lca_em(X, 2)
fit.item_prob # rows = classes, columns = items 1..82×8 Matrix{Float64}:
0.192319 0.176371 0.203465 0.2006 0.205411 0.211559 0.722236 0.684457
0.809866 0.775541 0.810785 0.777392 0.779226 0.834744 0.239161 0.356424import numpy as np, pandas as pd
X = pd.read_csv("data/sim/ch22-lca.csv").values
def lca_em(X, K=2, tol=1e-6, maxit=500):
N, J = X.shape
pi_ = np.full(K, 1 / K)
p = np.linspace(0.3, 0.7, K * J).reshape(K, J)
ll_old = -np.inf
for it in range(1, maxit + 1):
loglik = np.log(pi_) + X @ np.log(p).T + (1 - X) @ np.log(1 - p).T
m = loglik.max(axis=1, keepdims=True)
post = np.exp(loglik - m); post /= post.sum(axis=1, keepdims=True)
pi_ = post.mean(axis=0)
p = (post.T @ X) / post.sum(axis=0)[:, None]
ll = float((m + np.log(np.exp(loglik - m).sum(axis=1, keepdims=True))).sum())
if abs(ll - ll_old) < tol:
break
ll_old = ll
return dict(prevalence=pi_, item_prob=p, iterations=it)
fit = lca_em(X, 2)
fit["item_prob"].round(2) # rows = classes, columns = items 1..8array([[0.19, 0.18, 0.2 , 0.2 , 0.21, 0.21, 0.72, 0.68],
[0.81, 0.78, 0.81, 0.78, 0.78, 0.83, 0.24, 0.36]])Read the two rows: one class endorses items 1–6 heavily but items 7–8 rarely; the other is its mirror image. These classes are not “high” and “low” depression - they are two different shapes of response. A sum score, which only counts yeses, would blend them together; the mixture model pulls them apart. That is what LCA is for.
22.4 How Many Classes?
We told the estimator \(K = 2\), but in real life the number of classes is itself unknown - and it is the crux of the analysis. You fit \(K = 1, 2, 3, \dots\) and compare, usually with an information criterion (BIC or AIC) that rewards fit and penalizes complexity, plus interpretability (a class you cannot describe is probably an artifact). In practice you would use a maintained tool that also handles this:
# The real tool (needs the poLCA package — shown, not run on the book's server)
library(poLCA)
f <- cbind(i1, i2, i3, i4, i5, i6, i7, i8) ~ 1
m2 <- poLCA(f, data, nclass = 2) # fit 2, 3, 4 classes ...
m3 <- poLCA(f, data, nclass = 3) # ... and compare
c(m2$bic, m3$bic) # lower BIC = better; also judge interpretability
NoteFrom the Field
Kashdan, McKnight, and colleagues (Kashdan et al. 2009) applied latent-class analysis to the National Comorbidity Survey and found a neglected subtype of social-anxiety disorder: people who are anxious and risk-prone/approach-oriented, not the shy-avoidant stereotype. That subgroup is invisible to a symptom count - you can only see it once you stop asking “how much?” and start asking “what kind?” It is the same lesson as the antitypes case study: structure that a total score throws away can be the most clinically important thing in the data.
22.5 What This Chapter Teaches
- Dimensions vs. types. Factor analysis places people on continua; latent-class analysis sorts them into kinds. Ask which question your science is really asking.
- A latent class model is a mixture of response-pattern probabilities - classes you cannot see, recovered from responses you can.
- EM solves the chicken-and-egg by alternating soft assignment (E) and re-estimation (M), the same iterate-until-stable trick behind Rasch.
- Classes can differ in pattern, not just level - and those pattern-defined subgroups are exactly what a sum score is built to hide.
We have now met latent variables as dimensions (factors) and as categories (classes). The next chapter returns to the dimensional view but drills deeper into it - extending the Rasch model into the wider world of item response theory and generalizability, where the authors have published extensively.
Kashdan, Todd B., Patrick E. McKnight, J. Anthony Richey, and Stefan G. Hofmann. 2009. “When Social Anxiety Disorder Co-Exists with Risk-Prone, Approach Behavior: Investigating a Neglected, Meaningful Subset of People in the National Comorbidity Survey–Replication.” Behaviour Research and Therapy 47 (7): 559–68.