Latent Class Analysis — Foundations¶
The mixture model for categorical data, from scratch¶
Project 1 of the Latent Class Analysis arc. Latent class analysis (LCA) is the discrete-data mixture model: a population is a blend of a few unobserved classes, and each subject's categorical responses are generated by their class. It is the tool behind market segmentation, psychological typologies, and — as we will see in later projects — diagnostic testing without a gold standard. This notebook builds the core model from scratch two ways (maximum likelihood by EM, and Bayesian by a data-augmentation Gibbs sampler), fits it with PyMC, and validates all three against each other; the companion R notebook uses poLCA and BayesLCA.
The model¶
For $N$ subjects, $J$ binary items and $C$ latent classes: $$T_i\sim\text{Categorical}(\boldsymbol\lambda),\qquad x_{ij}\mid T_i=c\ \sim\ \text{Bernoulli}(\delta_{cj}).$$ The single defining assumption is local independence: within a class the items are independent, so all the association between items is explained by the class. Marginalising the latent class gives a mixture of product-Bernoulli components: $$\Pr(\mathbf x_i)=\sum_{c=1}^{C}\lambda_c\prod_{j=1}^{J}\delta_{cj}^{x_{ij}}(1-\delta_{cj})^{1-x_{ij}}.$$ The parameters are the class prevalences $\boldsymbol\lambda$ (how big each class is) and the item-response probabilities $\delta_{cj}$ (the "profile" of each class — the chance a member endorses item $j$).
Three ways to fit it¶
- Maximum likelihood by EM. The E-step computes each subject's posterior class membership (the responsibilities); the M-step updates $\boldsymbol\lambda$ and $\delta$. The likelihood is multimodal, so we use random restarts.
- Bayesian by data augmentation (Gibbs). Treat the latent class $T_i$ as a variable to impute. Then $\boldsymbol\lambda\mid T\sim\text{Dirichlet}$ and $\delta_{cj}\mid T\sim\text{Beta}$ are conjugate — three clean draws per sweep. This is exactly Congdon's Bayesian Models for Categorical Data Program 6.8.
- PyMC. NUTS cannot sample the discrete $T_i$, so we marginalise the class analytically (a log-sum-exp over the $C$ components) and let NUTS sample $\boldsymbol\lambda,\delta$ directly.
Both Bayesian routes hit label switching — independent chains can permute the class labels, so a pooled posterior is bimodal — which we demonstrate and fix by relabelling.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import lca, warnings; warnings.filterwarnings("ignore")
plt.rcParams.update({"figure.figsize": (9,4), "axes.grid": True, "grid.alpha": .22,
"axes.spines.top": False, "axes.spines.right": False, "font.size": 10.5})
BLUE, RED, GREEN, ORANGE, GREY, PURP = "#2b6cb0","#c53030","#2f855a","#dd6b20","#718096","#6b46c1"
rng = np.random.default_rng(0)
print("ready")
ready
1. Validation on simulated data¶
Before touching real data we check the estimators on data with a known truth: a 2-class, 4-item population where class 1 tends to answer "yes" on the first three items and "no" on the fourth, and class 2 the reverse. Both EM and the Gibbs sampler should recover the true $\boldsymbol\lambda$ and $\delta$.
lam_true = np.array([0.6, 0.4])
delta_true = np.array([[0.90, 0.85, 0.80, 0.10],
[0.15, 0.20, 0.10, 0.90]])
Xs, Ts = lca.simulate_lca(4000, lam_true, delta_true, rng)
em = lca.em_lca(Xs, 2, rng, n_starts=15)
gb = lca.gibbs_lca(Xs, 2, rng, draws=3000, burn=800)
tab = pd.DataFrame({
"λ true": lam_true, "λ EM": em["lam"].round(3), "λ Gibbs": gb["lam"].mean(0).round(3)})
print(tab.to_string(index=False))
print("\nδ (item-response probabilities): rows = class, cols = item")
print("true:\n", delta_true, "\nEM:\n", em["delta"].round(3), "\nGibbs (posterior mean):\n", gb["delta"].mean(0).round(3))
print("\nBoth estimators recover the truth to within Monte-Carlo error — the model and code are correct.")
λ true λ EM λ Gibbs
0.6 0.614 0.614
0.4 0.386 0.386
δ (item-response probabilities): rows = class, cols = item
true:
[[0.9 0.85 0.8 0.1 ]
[0.15 0.2 0.1 0.9 ]]
EM:
[[0.89 0.846 0.798 0.105]
[0.145 0.187 0.092 0.918]]
Gibbs (posterior mean):
[[0.89 0.846 0.798 0.105]
[0.146 0.187 0.093 0.917]]
Both estimators recover the truth to within Monte-Carlo error — the model and code are correct.
2. The Stouffer–Toby data — role conflict¶
The classic LCA teaching dataset (Stouffer & Toby, 1951; the values data in poLCA). 216 respondents each faced 4 scenarios (items A–D) posing a conflict between a universalistic norm (obligation to society/rules) and a particularistic one (obligation to a friend). We code the particularistic choice as $1$. The question LCA answers: are respondents a smooth continuum, or a mixture of a few types — say a "universalistic" type and a "particularistic" type?
df = pd.read_csv("stouffer_toby.csv") # items A,B,C,D coded 1/2
X = (df.values - 1).astype(int) # recode to 0/1 ; 1 = particularistic choice
N, J = X.shape; items = list(df.columns)
print(f"{N} respondents, {J} items {items}")
print("particularistic-response rate per item:", dict(zip(items, X.mean(0).round(3))))
# the 16 response patterns and their counts
pats, counts = np.unique(X, axis=0, return_counts=True)
order = np.argsort(counts)[::-1]
print("\nmost common response patterns (1=particularistic):")
for p, c in zip(pats[order][:6], counts[order][:6]):
print(" ", "".join(map(str, p)), "×", c)
216 respondents, 4 items ['A', 'B', 'C', 'D']
particularistic-response rate per item: {'A': np.float64(0.792), 'B': np.float64(0.5), 'C': np.float64(0.514), 'D': np.float64(0.31)}
most common response patterns (1=particularistic):
1111 × 42
1000 × 38
1100 × 25
1010 × 24
1110 × 23
0000 × 20
3. Maximum likelihood by EM¶
We fit the 2-class model by EM with 40 random restarts (the likelihood surface has multiple modes). The output is two class profiles — the probability that a member of each class makes the particularistic choice on each item — and the class prevalences. The profiles are what make the classes interpretable.
fit2 = lca.em_lca(X, 2, rng, n_starts=40)
print("class prevalences λ:", fit2["lam"].round(3))
prof = pd.DataFrame(fit2["delta"].round(3), columns=items,
index=[f"class {c+1} (λ={fit2['lam'][c]:.2f})" for c in range(2)])
print("\nitem-response probabilities δ (P[particularistic | class]):"); print(prof.to_string())
print(f"\nlog-likelihood {fit2['loglik']:.2f}, entropy R² {lca.entropy_R2(fit2['resp']):.3f} (class separation)")
# the signature LCA figure: class response profiles
fig, ax = plt.subplots(1, 2, figsize=(13, 4.2))
xpos = np.arange(J)
for c,col in zip(range(2), [BLUE, RED]):
ax[0].plot(xpos, fit2["delta"][c], "o-", color=col, lw=2, ms=8, label=f"class {c+1} (λ={fit2['lam'][c]:.2f})")
ax[0].set_xticks(xpos); ax[0].set_xticklabels(items); ax[0].set_ylim(0,1); ax[0].set_ylabel("P(particularistic)")
ax[0].set_title("Latent class profiles: the two 'types'"); ax[0].legend(fontsize=9)
# posterior classification: histogram of P(class 2 | responses)
p2 = fit2["resp"][:, 1]
ax[1].hist(p2, bins=25, color=GREY, edgecolor="white", lw=.3)
ax[1].set_title("Posterior class membership P(class 2 | responses)"); ax[1].set_xlabel("P(particularistic class)"); ax[1].set_ylabel("# respondents")
plt.tight_layout(); plt.show()
print("Class 1 (~72%) leans universalistic; class 2 (~28%) is strongly particularistic on every item — the classic 2-type result.")
class prevalences λ: [0.721 0.279]
item-response probabilities δ (P[particularistic | class]):
A B C D
class 1 (λ=0.72) 0.714 0.33 0.354 0.132
class 2 (λ=0.28) 0.993 0.94 0.927 0.769
log-likelihood -504.47, entropy R² 0.719 (class separation)
Class 1 (~72%) leans universalistic; class 2 (~28%) is strongly particularistic on every item — the classic 2-type result.
4. Bayesian LCA by data-augmentation Gibbs¶
The Bayesian fit treats the latent class $T_i$ as missing data to impute. With priors $\boldsymbol\lambda\sim\text{Dirichlet}(1,\dots,1)$ and $\delta_{cj}\sim\text{Beta}(1,1)$, each Gibbs sweep is three conjugate draws: $$T_i\sim\text{Categorical}(r_i),\qquad \boldsymbol\lambda\sim\text{Dirichlet}(1+n_c),\qquad \delta_{cj}\sim\text{Beta}(1+s_{cj},\,1+f_{cj}),$$ where $r_i$ are the responsibilities, $n_c$ the imputed class sizes, and $s_{cj},f_{cj}$ the within-class successes/failures. Unlike EM it returns full posterior distributions (with credible intervals) rather than point estimates.
# Label switching surfaces ACROSS chains: a single well-separated Gibbs chain stays put in
# one labelling, but independent random starts adopt different (mirror-image) labellings — so a
# POOLED posterior is bimodal until every draw is relabelled.
fix = lca.gibbs_lca(X, 2, rng, draws=4000, burn=1000, relabel=True) # working (relabelled) posterior
chains = [lca.gibbs_lca(X, 2, np.random.default_rng(100+s), draws=1500, burn=400, relabel=False) for s in range(6)]
fig, ax = plt.subplots(1, 2, figsize=(13, 4))
for ch in chains:
ax[0].plot(ch["lam"][:, 0], lw=.6, alpha=.85)
ax[0].axhline(0.5, color="k", ls=":", lw=1)
ax[0].set_title("λ₁ from 6 independent chains — no relabelling"); ax[0].set_xlabel("iteration"); ax[0].set_ylabel("λ₁"); ax[0].set_ylim(0,1)
pooled = np.concatenate([ch["lam"][:, 0] for ch in chains])
ax[1].hist(pooled, bins=45, color=GREY, alpha=.7, label="pooled, no relabelling")
ax[1].hist(np.maximum(pooled, 1-pooled), bins=45, color=GREEN, alpha=.7, label="after relabelling (λ decreasing)")
ax[1].set_title("Pooled λ₁ posterior: bimodal → unimodal"); ax[1].set_xlabel("λ₁"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print("Left: each chain stays in ONE labelling, but different chains disagree on which class is 'class 1' — some sit")
print("near 0.72, others near 0.28 (the mirror image). Right: pooling them gives a bimodal λ₁; relabelling every draw")
print("so λ is decreasing collapses it to one interpretable mode. (Congdon 6.8 instead builds the ordering into the")
print("sampler via constraints on δ, so its chains cannot switch in the first place.)")
Left: each chain stays in ONE labelling, but different chains disagree on which class is 'class 1' — some sit near 0.72, others near 0.28 (the mirror image). Right: pooling them gives a bimodal λ₁; relabelling every draw so λ is decreasing collapses it to one interpretable mode. (Congdon 6.8 instead builds the ordering into the sampler via constraints on δ, so its chains cannot switch in the first place.)
# posterior summaries (relabelled) with 95% credible intervals, vs EM
lam_pm = fix["lam"].mean(0); lam_ci = np.percentile(fix["lam"], [2.5, 97.5], axis=0)
print("class prevalences λ:")
for c in range(2):
print(f" class {c+1}: posterior mean {lam_pm[c]:.3f} 95% CI [{lam_ci[0,c]:.3f}, {lam_ci[1,c]:.3f}] (EM {fit2['lam'][c]:.3f})")
print("\nitem-response profiles δ with 95% CIs:")
dm = fix["delta"].mean(0); dlo, dhi = np.percentile(fix["delta"], [2.5, 97.5], axis=0)
fig, ax = plt.subplots(figsize=(9.5, 4.3)); xpos = np.arange(J)
for c,col,off in zip(range(2), [BLUE, RED], [-0.06, 0.06]):
ax.errorbar(xpos+off, dm[c], yerr=[dm[c]-dlo[c], dhi[c]-dm[c]], fmt="o", color=col, capsize=3, ms=7, lw=1.5, label=f"class {c+1}")
ax.plot(xpos+off, fit2["delta"][c], "_", color="k", ms=14, mew=2)
ax.set_xticks(xpos); ax.set_xticklabels(items); ax.set_ylim(0,1); ax.set_ylabel("P(particularistic)")
ax.set_title("Posterior profiles ±95% CI (points), with EM estimates (black dashes)"); ax.legend(fontsize=9); plt.show()
print("The Bayesian posterior means sit right on the EM estimates; the credible intervals quantify what EM leaves as points.")
class prevalences λ: class 1: posterior mean 0.684 95% CI [0.549, 0.794] (EM 0.721) class 2: posterior mean 0.316 95% CI [0.206, 0.451] (EM 0.279) item-response profiles δ with 95% CIs:
The Bayesian posterior means sit right on the EM estimates; the credible intervals quantify what EM leaves as points.
5. PyMC — marginalising the latent class¶
Sampling the discrete indicators $T_i$ with a gradient sampler (NUTS) is impossible, so the standard PyMC approach is to marginalise the class analytically. For each subject we compute the log-probability under every class and combine them with logsumexp:
$$\log\Pr(\mathbf x_i)=\operatorname{logsumexp}_c\Big[\log\lambda_c+\textstyle\sum_j\big(x_{ij}\log\delta_{cj}+(1-x_{ij})\log(1-\delta_{cj})\big)\Big],$$
added to the model as a pm.Potential. NUTS then samples the continuous $\boldsymbol\lambda,\delta$ efficiently. (This is the same marginalised likelihood the EM E-step and the Gibbs responsibilities are built from — three views of one object.)
import pymc as pm, pytensor.tensor as pt, arviz as az
C = 2; Xf = X.astype(float)
with pm.Model() as lca_model:
lam = pm.Dirichlet("lam", np.ones(C))
delta = pm.Beta("delta", 1.0, 1.0, shape=(C, J))
ld = pt.log(delta); l1 = pt.log1p(-delta)
comp = pt.log(lam)[None, :] + (Xf[:, None, :]*ld[None] + (1-Xf[:, None, :])*l1[None]).sum(-1) # (N,C)
pm.Potential("lik", pt.logsumexp(comp, axis=1).sum())
idata = pm.sample(1000, tune=1000, chains=4, cores=1, target_accept=0.92, random_seed=1, progressbar=False)
# relabel each posterior draw by prevalence (marginalised model is permutation-invariant)
la = idata.posterior["lam"].values.reshape(-1, C); de = idata.posterior["delta"].values.reshape(-1, C, J)
for i in range(len(la)):
o = np.argsort(la[i])[::-1]; la[i] = la[i][o]; de[i] = de[i][o]
print("PyMC posterior λ:", la.mean(0).round(3), " (Gibbs", fix["lam"].mean(0).round(3), ", EM", fit2["lam"].round(3), ")")
print("PyMC posterior δ:\n", de.mean(0).round(3))
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (4 chains in 1 job)
NUTS: [lam, delta]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 4 seconds.
PyMC posterior λ: [0.69 0.31] (Gibbs [0.684 0.316] , EM [0.721 0.279] ) PyMC posterior δ: [[0.707 0.322 0.348 0.128] [0.96 0.894 0.884 0.732]]
# the three engines agree — overlay the class-1 profile from EM, Gibbs and PyMC
fig, ax = plt.subplots(figsize=(9.5, 4.2)); xpos = np.arange(J)
ax.plot(xpos, fit2["delta"][0], "o-", color=BLUE, lw=2, ms=8, label="EM (MLE)")
ax.plot(xpos, fix["delta"].mean(0)[0], "s--", color=GREEN, lw=1.6, ms=7, label="Gibbs (from scratch)")
ax.plot(xpos, de.mean(0)[0], "^:", color=RED, lw=1.6, ms=8, label="PyMC (marginalised)")
ax.set_xticks(xpos); ax.set_xticklabels(items); ax.set_ylim(0,1); ax.set_ylabel("P(particularistic)")
ax.set_title("Class 1 profile — three independent engines agree"); ax.legend(fontsize=9); plt.show()
print("EM, the from-scratch Gibbs, and PyMC's marginalised NUTS give the same class profiles — mutual validation.")
EM, the from-scratch Gibbs, and PyMC's marginalised NUTS give the same class profiles — mutual validation.
6. How many classes, and does the model fit?¶
Two questions close the analysis. How many classes? We compare $C=1,2,3$ by the BIC ($-2\hat\ell+p\log N$), which penalises the extra parameters of larger models. Does it fit? Because $J=4$ gives only $2^4=16$ response patterns, we can compare the observed pattern counts to those expected under the fitted model, and summarise the discrepancy with the likelihood-ratio $G^2$ and Pearson $X^2$ statistics. (Project 2 develops model selection in depth on a larger dataset.)
rows = []
for C_ in [1, 2, 3]:
f = lca.em_lca(X, C_, rng, n_starts=40)
bic = -2*f["loglik"] + f["npar"]*np.log(N); aic = -2*f["loglik"] + 2*f["npar"]
g = lca.gof(X, f["lam"], f["delta"])
rows.append((C_, round(f["loglik"],2), f["npar"], round(aic,1), round(bic,1),
round(g["G2"],2), g["df"], round(lca.entropy_R2(f["resp"]),3) if C_>1 else np.nan))
sel = pd.DataFrame(rows, columns=["C","loglik","npar","AIC","BIC","G²","df","entropy R²"]).set_index("C")
print(sel.to_string()); best = sel["BIC"].idxmin()
print(f"\nBIC is minimised at C={best} classes — the classic 2-class Stouffer–Toby solution.")
# observed vs expected pattern frequencies at C=2
f2 = lca.em_lca(X, 2, rng, n_starts=40); g2 = lca.gof(X, f2["lam"], f2["delta"])
labels = ["".join(map(str, p.astype(int))) for p in g2["patterns"]]
o = np.argsort(g2["observed"])[::-1]
fig, ax = plt.subplots(figsize=(11, 4.2)); xp = np.arange(len(labels))
ax.bar(xp-0.2, g2["observed"][o], 0.4, color=BLUE, label="observed")
ax.bar(xp+0.2, g2["expected"][o], 0.4, color=ORANGE, label="expected (2-class model)")
ax.set_xticks(xp); ax.set_xticklabels([labels[i] for i in o], rotation=90, fontsize=7)
ax.set_title(f"Observed vs expected pattern counts (G²={g2['G2']:.2f}, df={g2['df']}): the model fits"); ax.set_xlabel("response pattern (A B C D)"); ax.legend(); plt.tight_layout(); plt.show()
print(f"G²={g2['G2']:.2f} on df={g2['df']} is small — the 2-class model reproduces the 16 response-pattern frequencies well.")
loglik npar AIC BIC G² df entropy R² C 1 -543.65 4 1095.3 1108.8 81.08 11 NaN 2 -504.47 9 1026.9 1057.3 2.72 6 0.719 3 -503.30 14 1034.6 1081.9 0.39 1 0.589 BIC is minimised at C=2 classes — the classic 2-class Stouffer–Toby solution.
G²=2.72 on df=6 is small — the 2-class model reproduces the 16 response-pattern frequencies well.
7. Summary¶
Latent class analysis fits a mixture of product-Bernoulli distributions: a few unobserved classes, each with its own item-endorsement profile, under the local-independence assumption that the classes explain all inter-item association. On the classic Stouffer–Toby role-conflict data, all three engines — EM (maximum likelihood), a from-scratch data-augmentation Gibbs sampler (Congdon 6.8), and PyMC with the class marginalised — agree on a two-class solution: a large (~72%) universalistic type and a smaller (~28%) strongly particularistic type, with BIC selecting $C=2$ and $G^2$ confirming the fit.
The transferable lessons. The latent class is a data-augmentation variable, which makes the Gibbs sampler fully conjugate (Dirichlet + Beta), tying LCA to the augmentation toolkit. Its discreteness forces the marginalise-the-class trick in gradient-based samplers like PyMC. And mixture models carry label switching, handled here by relabelling.
Next in the arc: choosing the number of classes rigorously (WAIC/LOO, bootstrap LRT, posterior-predictive checks) on Congdon's larger 17-item dataset; then latent class regression (class membership driven by covariates), and the Hui–Walter model for diagnostic testing without a gold standard.