Latent Class Regression (R)¶
The concomitant-variable model — poLCA's native formula, plus a from-scratch Gibbs¶
The R counterpart to lcareg_python.ipynb. Latent class regression lets the class weights depend on covariates through a multinomial logit, $\Pr(T_i=c\mid w_i)=\text{softmax}_c(w_i^\top\gamma_c)$, while the measurement model $x_{ij}\mid T_i=c\sim\text{Categorical}(\delta_{c,j,\cdot})$ is unchanged. poLCA is the reference implementation — hand it a formula with the covariate on the right and it fits exactly this concomitant model. We also port the data-augmentation Gibbs from scratch in base R (Dirichlet item updates + a Metropolis step for the membership coefficients) so the mechanics reproduce, then cross-check the two membership curves. Data: the poLCA election study (12 Gore/Bush trait ratings, 4 levels each; covariate = 7-point PARTY).
options(repr.plot.width=12, repr.plot.height=4.5)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths())); suppressMessages(library(poLCA))
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; GREY<-"#718096"; PURP<-"#6b46c1"
d <- read.csv("election.csv") # items coded 1..4, PARTY 1..7
items <- c("MORALG","CARESG","KNOWG","LEADG","DISHONG","INTELG",
"MORALB","CARESB","KNOWB","LEADB","DISHONB","INTELB")
cat(nrow(d), "respondents,", length(items), "items (4 levels), covariate PARTY 1..7\n")
print(table(d$PARTY))
Warning message: "package 'poLCA' was built under R version 4.6.1"
1294 respondents, 12 items (4 levels), covariate PARTY 1..7
1 2 3 4 5 6 7 256 191 196 111 180 162 198
1. poLCA with a formula — the concomitant model¶
cbind(item1, …, item12) ~ PARTY tells poLCA to make class membership a multinomial logit of party. nrep restarts guard against local optima. We read off the three classes' candidate ratings to label them (pro-Gore / pro-Bush / ambivalent; lower level = describes better) and plot the fitted membership curve $\Pr(\text{class}\mid\text{party})$.
set.seed(1)
f <- cbind(MORALG,CARESG,KNOWG,LEADG,DISHONG,INTELG,
MORALB,CARESB,KNOWB,LEADB,DISHONB,INTELB) ~ PARTY
m <- poLCA(f, d, nclass=3, nrep=10, maxiter=4000, verbose=FALSE)
# label classes by mean Gore vs Bush rating (probs is a list: item -> class x level)
lev <- 1:4
rate <- sapply(1:12, function(j) as.vector(m$probs[[j]] %*% lev)) # (class x item)
gore <- rowMeans(rate[,1:6]); bush <- rowMeans(rate[,7:12])
lab <- rep("ambivalent",3); lab[which.min(gore-bush)]<-"pro-Gore"; lab[which.max(gore-bush)]<-"pro-Bush"
cat("class labels:", paste(lab,collapse=", "), "\n")
cat("mean Gore rating:", round(gore,2), " | mean Bush rating:", round(bush,2), "(1=best..4=worst)\n\n")
# predicted membership across party from the fitted logit coefficients (class 1 = reference)
grid <- cbind(1, 1:7)
eta <- cbind(0, grid %*% m$coeff) # (7 x 3)
curveP <- exp(eta)/rowSums(exp(eta))
colP <- c("pro-Gore"=BLUE,"pro-Bush"=RED,"ambivalent"=GREY)
par(mfrow=c(1,2), mar=c(4,4,3,1))
matplot(1:7, curveP, type="n", ylim=c(0,1), xlab="PARTY (1=strong Dem .. 7=strong Rep)", ylab="P(class | party)",
main="poLCA membership regression")
for(c in 1:3) lines(1:7, curveP[,c], type="b", pch=19, lwd=2.2, col=colP[lab[c]])
legend("topright", lab, col=colP[lab], lwd=2, pch=19, bty="n")
barplot(rbind(4-gore,4-bush), beside=TRUE, names.arg=lab, col=c(BLUE,RED),
ylab="favourability (4 - mean rating)", main="Candidate favourability by class")
legend("topright", c("Gore","Bush"), fill=c(BLUE,RED), bty="n")
par(mfrow=c(1,1))
cat("The pro-Bush class rises steeply with PARTY and the pro-Gore class falls: party identification\n")
cat("re-weights the latent vote-classes -- the concomitant effect a fixed-weight LCA cannot express.\n")
class labels: pro-Gore, pro-Bush, ambivalent
mean Gore rating: 1.81 2.58 2.33 | mean Bush rating: 2.7 2.01 2.54 (1=best..4=worst)
The pro-Bush class rises steeply with PARTY and the pro-Gore class falls: party identification
re-weights the latent vote-classes -- the concomitant effect a fixed-weight LCA cannot express.
2. From-scratch Gibbs in base R¶
The same data-augmentation sampler as the Python engine: $\delta\mid T$ is conjugate Dirichlet, $T\mid\delta,\gamma$ is Categorical, and the membership coefficients $\gamma$ get a per-class random-walk Metropolis update (multinomial-logit regression of the imputed labels on the covariate). Membership is identified by a ridge prior $\gamma_c\sim N(0,\sigma^2)$; we centre party so the intercept is at a moderate voter.
rdir <- function(a){ g<-rgamma(length(a),a,1); g/sum(g) }
softmax_rows <- function(A){ A<-A-apply(A,1,max); E<-exp(A); E/rowSums(E) }
item_loglik <- function(X, delta){ # X:(N,J) in 1..L ; delta:(C,J,L) -> (N,C)
N<-nrow(X); J<-ncol(X); C<-dim(delta)[1]; out<-matrix(0,N,C)
for(c in 1:C){ ld<-log(delta[c,,]); out[,c]<-rowSums(sapply(1:J, function(j) ld[j, X[,j]])) }
out }
mnl_ll <- function(gamma, Tt, W){ eta<-W %*% t(gamma); sum(eta[cbind(1:length(Tt),Tt)] - log(rowSums(exp(eta)))) }
lcareg_gibbs <- function(X, W, C, draws=2000, burn=2000, sd_gamma=5, step=0.08){
N<-nrow(X); J<-ncol(X); p<-ncol(W); L<-max(X)
Tt<-sample(1:C, N, replace=TRUE); gamma<-matrix(0,C,p); Bp<-1/sd_gamma^2
steps<-rep(step,C); acc<-rep(0,C); att<-rep(0,C)
GS<-array(0,c(draws,C,p)); DS<-array(0,c(draws,C,J,L)); llc<-mnl_ll(gamma,Tt,W)
for(it in 1:(draws+burn)){
delta<-array(0,c(C,J,L))
for(c in 1:C){ Xc<-X[Tt==c,,drop=FALSE]
for(j in 1:J){ cnt<-tabulate(Xc[,j],nbins=L); delta[c,j,]<-rdir(1+cnt) } }
lp<-W %*% t(gamma); lp<-lp-log(rowSums(exp(lp))); lp<-lp+item_loglik(X,delta)
P<-softmax_rows(lp); u<-runif(N); Tt<-max.col(-(matrix(u,N,C) > t(apply(P,1,cumsum)))*0 + (t(apply(P,1,cumsum))>=u), "first")
llc<-mnl_ll(gamma,Tt,W)
for(c in 1:C){ prop<-gamma; prop[c,]<-gamma[c,]+steps[c]*rnorm(p); llp<-mnl_ll(prop,Tt,W)
logr<-(llp-llc)-0.5*Bp*(sum(prop[c,]^2)-sum(gamma[c,]^2)); att[c]<-att[c]+1
if(log(runif(1))<logr){ gamma<-prop; llc<-llp; acc[c]<-acc[c]+1 }
if(it<burn && att[c]>=50){ r<-acc[c]/att[c]; steps[c]<-min(max(steps[c]*exp((r-0.3)*0.5),1e-3),2); acc[c]<-0; att[c]<-0 } }
if(it>burn){ GS[it-burn,,]<-gamma; DS[it-burn,,,]<-delta } }
# relabel to first draw by item profile
ref<-DS[1,,,]; perms<-list(c(1,2,3),c(1,3,2),c(2,1,3),c(2,3,1),c(3,1,2),c(3,2,1))
for(dd in 1:draws){ best<-Inf; bp<-1:C
for(pm in perms){ dist<-sum((DS[dd,pm,,]-ref)^2); if(dist<best){best<-dist; bp<-pm} }
DS[dd,,,]<-DS[dd,bp,,]; GS[dd,,]<-GS[dd,bp,] }
list(gamma=apply(GS,c(2,3),mean), delta=apply(DS,c(2,3,4),mean)) }
X <- as.matrix(d[,items]); Wc <- cbind(1, d$PARTY-4)
set.seed(2); fit <- lcareg_gibbs(X, Wc, 3, draws=2500, burn=2500)
Gm<-fit$gamma; Dm<-fit$delta
rate2<-sapply(1:12,function(j) as.vector(Dm[,j,] %*% lev)); g2<-rowMeans(rate2[,1:6]); b2<-rowMeans(rate2[,7:12])
lab2<-rep("ambivalent",3); lab2[which.min(g2-b2)]<-"pro-Gore"; lab2[which.max(g2-b2)]<-"pro-Bush"
gridc<-cbind(1, (1:7)-4); curveG<-softmax_rows(gridc %*% t(Gm))
cat("from-scratch class labels:", paste(lab2,collapse=", "), "\n")
cat("party log-odds slopes (raw — NOT individually identified; only differences between classes are):",
paste(sprintf("%s %.2f",lab2,Gm[,2]),collapse=" | "), "\n")
gref <- which(lab2=="pro-Gore")
cat("identified contrasts vs the pro-Gore class:",
paste(sprintf("%s %+.2f",lab2[-gref],Gm[-gref,2]-Gm[gref,2]),collapse=" | "),
" <- compare with the Python notebook\n")
from-scratch class labels: ambivalent, pro-Gore, pro-Bush
party log-odds slopes (raw — NOT individually identified; only differences between classes are): ambivalent -0.87 | pro-Gore -1.44 | pro-Bush -0.06
identified contrasts vs the pro-Gore class: ambivalent +0.57 | pro-Bush +1.38 <- compare with the Python notebook
3. Do the two engines agree?¶
Overlay poLCA's membership curve (its native concomitant fit) with the from-scratch Gibbs curve. Same shape, same crossing point — the vote-classes swing from pro-Gore to pro-Bush across the party scale in both.
ord<-match(lab, lab2) # align from-scratch classes to poLCA labels for overlay
par(mfrow=c(1,1), mar=c(4,4,3,1))
matplot(1:7, curveP, type="n", ylim=c(0,1), xlab="PARTY (1=strong Dem .. 7=strong Rep)",
ylab="P(class | party)", main="Membership curves: poLCA (solid) vs from-scratch Gibbs (dashed)")
for(c in 1:3){ lines(1:7, curveP[,c], type="b", pch=19, lwd=2.4, col=colP[lab[c]])
lines(1:7, curveG[,ord[c]], type="b", pch=1, lty=2, lwd=1.7, col=colP[lab[c]]) }
legend("topright", c(lab,"poLCA","Gibbs"), col=c(colP[lab],"black","black"),
lwd=c(2,2,2,2,1.7), pch=c(19,19,19,19,1), lty=c(1,1,1,1,2), bty="n", cex=.85)
cat("poLCA and the from-scratch Gibbs give the same party gradient: strong Democrats sit in the pro-Gore class,\n")
cat("strong Republicans in the pro-Bush class, the ambivalent class peaking in the middle.\n")
poLCA and the from-scratch Gibbs give the same party gradient: strong Democrats sit in the pro-Gore class,
strong Republicans in the pro-Bush class, the ambivalent class peaking in the middle.
4. Summary¶
poLCA's formula interface fits latent class regression directly — cbind(items) ~ PARTY makes membership a multinomial logit of the covariate — and the from-scratch base-R Gibbs (Dirichlet items + a Metropolis step for the membership coefficients) reproduces the same three-class solution and the same party gradient. Both agree with the Python engine and PyMC: the 2000 electorate splits into pro-Gore, pro-Bush, and ambivalent latent classes, and party identification sharply re-weights them — a respondent-specific membership that ordinary fixed-weight LCA cannot express.
The covariate could as easily be education or age; poLCA also allows several covariates and interactions in the formula. This closes the latent class regression project and links the arc to the multinomial-logit models, whose likelihood is exactly the membership model used here.