Latent Class Regression¶
Letting covariates predict class membership — the concomitant-variable LCA¶
Ordinary latent class analysis gives every subject the same prior class weights $\lambda_c$. But often we have covariates — age, education, party identification — that ought to shift who lands in which class. Latent class regression (the concomitant-variable model of Dayton & Macready 1988; exactly what poLCA fits when you hand it a formula) replaces the fixed weights with a multinomial logit of the covariates:
$$\Pr(T_i = c \mid w_i) \;=\; \frac{\exp(w_i^\top \gamma_c)}{\sum_{k=1}^{C}\exp(w_i^\top \gamma_k)}\qquad\text{(membership model)}$$ $$x_{ij}\mid T_i=c \;\sim\; \text{Categorical}(\delta_{c,j,\cdot})\qquad\text{(measurement model)}$$
The $\gamma_c$ are membership log-odds coefficients; $\delta_{c,j,\cdot}$ is class $c$'s response-probability vector for item $j$. Setting every $\gamma_c=0$ recovers ordinary LCA. The measurement part is unchanged — the classes are still defined by the response patterns — but now the mix of classes bends with the covariate.
From scratch we sample it with a data-augmentation Gibbs: $\delta\mid T$ is conjugate Dirichlet, $T\mid\delta,\gamma$ is conjugate Categorical, and $\gamma\mid T$ is a multinomial-logit regression of the imputed labels on $W$, updated by a small random-walk Metropolis step (the one non-conjugate block). We validate on simulated data, then fit the classic poLCA election data — 12 trait ratings of Gore and Bush — with party identification as the covariate, and watch the latent vote-classes reorganise across the political spectrum. A PyMC marginalised-class fit and (in the companion notebook) poLCA itself confirm the result.
Data: 2000 American National Election Study, 1294 complete respondents, 12 items (6 traits × 2 candidates) each rated 1–4 (1 = describes extremely well), plus 7-point PARTY (1 = strong Democrat … 7 = strong Republican).
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import lcareg as L
rng = np.random.default_rng(7)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; PURP="#6b46c1"; GREY="#718096"
d = pd.read_csv("election.csv")
items = ["MORALG","CARESG","KNOWG","LEADG","DISHONG","INTELG",
"MORALB","CARESB","KNOWB","LEADB","DISHONB","INTELB"]
X = d[items].to_numpy() - 1 # recode 1..4 -> 0..3
party = d["PARTY"].to_numpy().astype(float) # 1=strong Dem .. 7=strong Rep
W = np.column_stack([np.ones(len(d)), party - 4.0]) # intercept + centred party
print(f"{X.shape[0]} respondents, {X.shape[1]} items (4 levels each), covariate = 7-point PARTY")
print("party distribution 1..7:", np.bincount(d['PARTY'])[1:])
1294 respondents, 12 items (4 levels each), covariate = 7-point PARTY party distribution 1..7: [256 191 196 111 180 162 198]
1. Does it work? — recovering a known membership regression¶
Simulated data with $C=3$ classes, 8 four-level items, and a covariate that pushes class 0 toward low values and class 2 toward high values. We fit from scratch and check that the estimated membership curves $\Pr(T=c\mid w)$ and the item profiles land on the truth (raw $\gamma$ is only identified up to a constant shift across classes, so we compare the identified probabilities, and align labels by item profile).
from itertools import permutations
C,J,Lv = 3,8,4
gamma_true = np.array([[1.5,-0.6],[0.0,0.0],[-1.5,0.6]])
delta_true = np.zeros((C,J,Lv))
for c in range(C):
for j in range(J):
b=np.ones(Lv); b[(c+j)%Lv]=7.0; delta_true[c,j]=b/b.sum()
pcov = rng.integers(1,8,1600).astype(float)
Xs,Ws,Ts = L.simulate_lcareg(1600, gamma_true, delta_true, rng, covariate=pcov-4.0)
sim = L.lcareg_gibbs(Xs, Ws, C, rng, draws=1500, burn=1500)
G = sim['gamma'].mean(0); D = sim['delta'].mean(0)
# align estimated labels to truth by item profile
pm_best=min(permutations(range(C)), key=lambda p:((D[list(p)]-delta_true)**2).sum()); pm=list(pm_best)
grid = np.column_stack([np.ones(7), np.arange(1,8)-4.0])
tc = L.membership_probs(gamma_true, grid); ec = L.membership_probs(G[pm], grid)
lp = np.log(L.membership_probs(G,Ws))+L._item_loglik(Xs,D); ari=L.adjusted_rand(Ts,lp.argmax(1))
print(f"membership-curve max error {np.abs(tc-ec).max():.3f} | profile max error {np.abs(D[pm]-delta_true).max():.3f} | class-recovery ARI {ari:.2f}")
fig,ax=plt.subplots(1,2,figsize=(12,4.2))
cols=[BLUE,GREEN,RED]; xp=np.arange(1,8)
for c in range(C):
ax[0].plot(xp,tc[:,c],'-',color=cols[c],lw=2)
ax[0].plot(xp,ec[:,c],'o--',color=cols[c],ms=6,mfc='white')
ax[0].set_title("Membership $P(T=c\\,|\\,$covariate$)$: truth (line) vs estimated (dots)")
ax[0].set_xlabel("covariate value"); ax[0].set_ylabel("class probability"); ax[0].set_ylim(0,1)
ax[0].legend([plt.Line2D([],[],color=cols[c],lw=2) for c in range(C)],[f"class {c}" for c in range(C)],frameon=False,fontsize=8)
ax[1].scatter(delta_true.ravel(), D[pm].ravel(), s=14, color=PURP, alpha=.6)
ax[1].plot([0,1],[0,1],'k--',lw=1); ax[1].set_title("Item response probabilities: true vs estimated")
ax[1].set_xlabel("true $\\delta$"); ax[1].set_ylabel("estimated $\\delta$")
plt.tight_layout(); plt.show()
print("The estimated membership curves and item profiles sit on the truth: the sampler recovers both the")
print("measurement model AND how the covariate re-weights the classes.")
membership-curve max error 0.022 | profile max error 0.058 | class-recovery ARI 0.94
The estimated membership curves and item profiles sit on the truth: the sampler recovers both the measurement model AND how the covariate re-weights the classes.
2. The election data — do party lines reshape the latent vote-classes?¶
We fit $C=3$ classes to the 12 candidate-trait ratings with party as the membership covariate. The three classes turn out to be a pro-Gore, a pro-Bush, and an ambivalent/mixed typology (read off each class's average rating of the two candidates — lower level = describes better). The regression then tells us how membership in each shifts as we move from strong Democrat to strong Republican.
fit = L.lcareg_gibbs(X, W, 3, rng, draws=3000, burn=3000)
G = fit['gamma'].mean(0); D = fit['delta'].mean(0) # (3,2) and (3,12,4)
levels = np.arange(1,5)
gore = np.array([ (D[c,:6]*levels).sum(1).mean() for c in range(3) ]) # mean rating of Gore items (low=favourable)
bush = np.array([ (D[c,6:]*levels).sum(1).mean() for c in range(3) ])
lab = {}
lab[int(np.argmin(gore-bush))]="pro-Gore"; lab[int(np.argmax(gore-bush))]="pro-Bush"
lab[[c for c in range(3) if c not in lab][0]]="ambivalent"
names=[lab[c] for c in range(3)]
print("class mean Gore rating mean Bush rating (1=best .. 4=worst)")
for c in range(3): print(f" {names[c]:12s} {gore[c]:.2f} {bush[c]:.2f}")
grid = np.column_stack([np.ones(7), np.arange(1,8)-4.0])
curve = L.membership_probs(G, grid) # P(class | party) for party 1..7
cols={'pro-Gore':BLUE,'pro-Bush':RED,'ambivalent':GREY}
fig,ax=plt.subplots(1,2,figsize=(12.5,4.4))
for c in range(3):
ax[0].plot(range(1,8), curve[:,c], 'o-', color=cols[names[c]], lw=2.2, ms=6, label=names[c])
ax[0].set_xlabel("PARTY (1 = strong Democrat → 7 = strong Republican)"); ax[0].set_ylabel("P(class | party)")
ax[0].set_title("Latent vote-class membership across the party spectrum"); ax[0].set_ylim(0,1); ax[0].legend(frameon=False)
# candidate favourability by class
xb=np.arange(3); ax[1].bar(xb-0.18,4-gore,0.36,color=BLUE,label="Gore favourability"); ax[1].bar(xb+0.18,4-bush,0.36,color=RED,label="Bush favourability")
ax[1].set_xticks(xb); ax[1].set_xticklabels(names); ax[1].set_ylabel("favourability (4 − mean rating)")
ax[1].set_title("How each latent class rates the two candidates"); ax[1].legend(frameon=False)
plt.tight_layout(); plt.show()
gorec = [c for c in range(3) if names[c]=="pro-Gore"][0]
print(f"\nParty log-odds slope, pro-Bush vs pro-Gore: {G[[c for c in range(3) if names[c]=='pro-Bush'][0],1]-G[gorec,1]:+.2f} per party point")
print("(the identified contrast). The pro-Bush membership curve rises with PARTY while pro-Gore falls —")
print("exactly the concomitant effect a fixed-weight LCA cannot express.")
class mean Gore rating mean Bush rating (1=best .. 4=worst) pro-Bush 2.58 2.02 pro-Gore 1.81 2.69 ambivalent 2.33 2.55
Party log-odds slope, pro-Bush vs pro-Gore: +1.38 per party point (the identified contrast). The pro-Bush membership curve rises with PARTY while pro-Gore falls — exactly the concomitant effect a fixed-weight LCA cannot express.
The three classes are substantively clean — one rates Gore highly and Bush poorly, one the reverse, one is lukewarm on both — and party identification sharply re-weights them: strong Democrats are overwhelmingly in the pro-Gore class, strong Republicans in the pro-Bush class, with the ambivalent class most common in the middle. A fixed-weight LCA would report a single average mix and miss this entirely; the regression turns the class weights into a function of who the respondent is.
3. What the covariate buys — regression vs fixed-weight LCA¶
To make the gain explicit we refit with the covariate switched off (intercept only, $W=\mathbf 1$) and compare the implied membership. The fixed-weight model collapses every respondent onto the same class mix; the regression lets it swing across the spectrum. We also report the fitted party log-odds slopes $\gamma_{c,\text{party}}$ — the heart of the model.
fit0 = L.lcareg_gibbs(X, np.ones((len(d),1)), 3, rng, draws=1500, burn=1500)
G0 = fit0['gamma'].mean(0); D0 = fit0['delta'].mean(0)
# align the intercept-only classes to the regression classes by profile
pm0 = list(min(permutations(range(3)), key=lambda p: ((D0[list(p)]-D)**2).sum()))
base = L.membership_probs(G0[pm0], np.ones((1,1)))[0] # single fixed mix
gorec = [c for c in range(3) if names[c]=="pro-Gore"][0]
print("party log-odds slope RELATIVE to the pro-Gore class (identified contrast, per party point):")
for c in range(3):
if c!=gorec: print(f" {names[c]:12s} vs pro-Gore {G[c,1]-G[gorec,1]:+.2f}")
print(f"\nfixed-weight LCA mix (no covariate): " + ", ".join(f"{names[c]} {base[c]:.2f}" for c in range(3)))
fig,ax=plt.subplots(figsize=(8.5,4.3))
for c in range(3):
ax.plot(range(1,8), curve[:,c], 'o-', color=cols[names[c]], lw=2.2, ms=5, label=f"{names[c]} (regression)")
ax.hlines(base[c], 1, 7, color=cols[names[c]], ls=':', lw=1.6)
ax.set_xlabel("PARTY"); ax.set_ylabel("P(class)"); ax.set_ylim(0,1)
ax.set_title("Regression membership (solid) vs one fixed mix (dotted)"); ax.legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print("The dotted lines are what fixed-weight LCA assumes for everyone; the solid curves are the regression's")
print("respondent-specific membership. The party slopes are large and opposite for the pro-Gore/pro-Bush classes.")
party log-odds slope RELATIVE to the pro-Gore class (identified contrast, per party point): pro-Bush vs pro-Gore +1.38 ambivalent vs pro-Gore +0.57 fixed-weight LCA mix (no covariate): pro-Bush 0.30, pro-Gore 0.28, ambivalent 0.42
The dotted lines are what fixed-weight LCA assumes for everyone; the solid curves are the regression's respondent-specific membership. The party slopes are large and opposite for the pro-Gore/pro-Bush classes.
4. The same model in PyMC — marginalising the class¶
NUTS cannot sample the discrete label $T_i$, so we marginalise it: the observed-data log-likelihood sums over classes,
$$\log p(x_i\mid w_i,\gamma,\delta)=\log\!\sum_{c=1}^{C}\underbrace{\text{softmax}_c(w_i^\top\gamma)}_{\text{membership}}\;\prod_{j}\delta_{c,j,x_{ij}},$$
added with pm.Potential via logsumexp. We pin the reference class ($\gamma_0=0$) for identification, put a Dirichlet on each $\delta_{c,j,\cdot}$, and recover the same party gradient. (One-hot encoding the responses makes the item likelihood a single tensordot.)
import pymc as pm, pytensor.tensor as pt
Xoh = np.eye(4)[X].astype(float) # (N,12,4) one-hot
Wc = W.copy() # intercept + centred party
with pm.Model() as mod:
g_free = pm.Normal("g_free", 0, 3, shape=(2,2)) # C-1 classes x p ; reference class = 0
gamma = pt.concatenate([pt.zeros((1,2)), g_free], axis=0) # (3,2)
delta = pm.Dirichlet("delta", a=np.ones(4), shape=(3,12,4))
eta = pt.dot(Wc, gamma.T) # (N,3)
log_mem = eta - pt.logsumexp(eta, axis=1, keepdims=True)
item_ll = pt.tensordot(Xoh, pt.log(delta), axes=[[1,2],[1,2]]) # (N,3)
tot = log_mem + item_ll
pm.Potential("like", pt.sum(pt.logsumexp(tot, axis=1)))
idata = pm.sample(1000, tune=1000, chains=4, target_accept=0.9, random_seed=11, progressbar=False)
# PyMC mixtures label-switch ACROSS chains, so relabel EVERY draw to the from-scratch
# profiles D before averaging (aligning the pooled mean would blend switched chains).
gdraw = idata.posterior["g_free"].stack(s=("chain","draw")).transpose("s",...).values # (S,2,2)
ddraw = idata.posterior["delta"].stack(s=("chain","draw")).transpose("s",...).values # (S,3,12,4)
perms = [list(p) for p in permutations(range(3))]
Gacc = np.zeros((3,2))
for s in range(ddraw.shape[0]):
gs = np.vstack([np.zeros((1,2)), gdraw[s]]) # (3,2) full gamma this draw
pm = min(perms, key=lambda p: ((ddraw[s][p]-D)**2).sum())
Gacc += gs[pm]
Gp = Gacc/ddraw.shape[0]
curve_pm = L.membership_probs(Gp, grid)
fig,ax=plt.subplots(figsize=(8.5,4.3))
for c in range(3):
ax.plot(range(1,8), curve[:,c], 'o-', color=cols[names[c]], lw=2.2, ms=6, label=f"{names[c]} (Gibbs)")
ax.plot(range(1,8), curve_pm[:,c],'s--',color=cols[names[c]], lw=1.6, ms=5, mfc='white')
ax.set_xlabel("PARTY"); ax.set_ylabel("P(class | party)"); ax.set_ylim(0,1)
ax.set_title("Membership curves: from-scratch Gibbs (solid) vs PyMC marginalised (dashed)"); ax.legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
gorec=[c for c in range(3) if names[c]=="pro-Gore"][0]; bushc=[c for c in range(3) if names[c]=="pro-Bush"][0]
print(f"PyMC party log-odds slope, pro-Bush vs pro-Gore: {Gp[bushc,1]-Gp[gorec,1]:+.2f} (matches the from-scratch Gibbs contrast)")
print("Two samplers, one story: membership swings from the pro-Gore to the pro-Bush class across the party scale.")
print("\nOn any sampler warnings above: the class labels in a mixture are not identified, so per-component")
print("R-hat and ESS look poor by construction. The quantity reported here — the pro-Bush vs pro-Gore party")
print("slope — is a CONTRAST between classes, which is identified, and it matches the from-scratch Gibbs.")
g++ not available, if using conda: `conda install gxx`
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [g_free, delta]
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 18 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters. A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
PyMC party log-odds slope, pro-Bush vs pro-Gore: +1.38 (matches the from-scratch Gibbs contrast) Two samplers, one story: membership swings from the pro-Gore to the pro-Bush class across the party scale. On any sampler warnings above: the class labels in a mixture are not identified, so per-component R-hat and ESS look poor by construction. The quantity reported here — the pro-Bush vs pro-Gore party slope — is a CONTRAST between classes, which is identified, and it matches the from-scratch Gibbs.
5. Summary¶
Latent class regression keeps LCA's measurement model but lets the class weights depend on covariates through a multinomial logit — the concomitant-variable model. On the 2000 election data the three latent vote-classes (pro-Gore, pro-Bush, ambivalent) are exactly what candidate-trait ratings should produce, and party identification re-weights them sharply: the membership curves cross as we move from strong Democrat to strong Republican, a respondent-specific mix that fixed-weight LCA cannot represent.
From scratch the fit needed only one new ingredient over ordinary LCA — a Metropolis step for the membership coefficients $\gamma$ inside the same data-augmentation Gibbs (Dirichlet item updates, Categorical label updates). The identified quantities are the membership probabilities $\Pr(T=c\mid w)$, which we validated against a known simulation (membership-curve error under $0.05$ and class recovery ARI $\approx 0.95$) and reproduced with a PyMC marginalised-class model and, in the companion notebook, poLCA's native formula interface.
This is the bridge from LCA to the regression world: the covariate can equally be education, age, or income; the machinery is identical. It also sets up latent transition analysis (the covariate becomes time, the classes evolve) and connects back to the multinomial-logit arc, whose likelihood is exactly the membership model here.