Latent Class Analysis — Foundations (R)¶
The mixture model for categorical data, via poLCA & BayesLCA¶
The R counterpart to lca_python.ipynb. We fit the same latent class model to the Stouffer–Toby role-conflict data two ways with established R packages: poLCA (maximum likelihood by EM — the standard LCA package) and BayesLCA (the Bayesian fit by Gibbs, the analogue of our from-scratch sampler). Together they cross-validate the Python results in a second language and against reference implementations.
Recall the model: $T_i\sim\text{Categorical}(\boldsymbol\lambda)$, $x_{ij}\mid T_i=c\sim\text{Bernoulli}(\delta_{cj})$, under local independence. poLCA expects items coded as consecutive integers starting at 1 (our data is 1/2, with 2 = particularistic); BayesLCA expects a 0/1 matrix.
options(repr.plot.width=12, repr.plot.height=4.2)
.libPaths(Sys.getenv("R_LIBS_USER"))
suppressMessages({library(poLCA); library(BayesLCA)})
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; ORANGE<-"#dd6b20"; GREY<-"#718096"
d <- read.csv("stouffer_toby.csv") # items A,B,C,D coded 1/2 (2 = particularistic)
items <- names(d); N <- nrow(d)
cat(N, "respondents,", ncol(d), "items:", items, "\n")
cat("particularistic-response rate per item:\n"); print(round(colMeans(d==2),3))
Warning message: "package 'poLCA' was built under R version 4.6.1"
Warning message: "package 'BayesLCA' was built under R version 4.6.1"
216 respondents, 4 items: A B C D
particularistic-response rate per item:
A B C D 0.792 0.500 0.514 0.310
1. Maximum likelihood by EM — poLCA¶
poLCA fits the latent class model by EM. We request 2 classes with 30 random restarts (nrep) to guard against local maxima, and read off the class prevalences ($P$), the item-response probabilities (probr), and the fit statistics. We reorder the classes by prevalence so the labelling matches the Python notebook.
set.seed(1)
f <- cbind(A,B,C,D) ~ 1
m2 <- poLCA(f, data=d, nclass=2, nrep=30, verbose=FALSE)
# reorder classes by prevalence (descending) for a canonical labelling
ord <- order(m2$P, decreasing=TRUE)
Pc <- m2$P[ord]
# probs[[j]] is class x outcome; column 2 = P(response 2 = particularistic)
prof <- sapply(m2$probs, function(m) m[ord, 2]) # classes x items
rownames(prof) <- paste0("class", 1:2, " (lambda=", round(Pc,2), ")"); colnames(prof) <- items
cat("class prevalences lambda:", round(Pc,3), "\n\nitem-response probabilities P(particularistic | class):\n"); print(round(prof,3))
cat(sprintf("\nlog-lik %.2f | df %d | G^2 %.2f | X^2 %.2f | BIC %.1f | AIC %.1f\n", m2$llik, m2$resid.df, m2$Gsq, m2$Chisq, m2$bic, m2$aic))
# class-profile figure
matplot(t(prof), type="b", pch=19, lwd=2, col=c(BLUE,RED), lty=1, ylim=c(0,1), xaxt="n", xlab="item", ylab="P(particularistic)", main="poLCA class profiles: the two 'types'")
axis(1, at=1:4, labels=items); legend("topright", rownames(prof), col=c(BLUE,RED), lwd=2, pch=19, bty="n", cex=.85)
class prevalences lambda: 0.721 0.279 item-response probabilities P(particularistic | class):
A B C D class1 (lambda=0.72) 0.714 0.33 0.354 0.132 class2 (lambda=0.28) 0.993 0.94 0.927 0.769
log-lik -504.47 | df 6 | G^2 2.72 | X^2 2.72 | BIC 1057.3 | AIC 1026.9
2. How many classes? — BIC across models¶
Fitting $C=1,2,3$ and comparing the BIC selects the number of classes. As in the Python notebook, the BIC is minimised at two classes.
set.seed(2)
res <- data.frame()
for(C in 1:3){ m <- poLCA(f, data=d, nclass=C, nrep=30, verbose=FALSE)
res <- rbind(res, data.frame(C=C, loglik=round(m$llik,2), npar=m$npar, AIC=round(m$aic,1), BIC=round(m$bic,1), Gsq=round(m$Gsq,2), df=m$resid.df)) }
print(res)
cat(sprintf("\nBIC minimised at C = %d classes -- the classic 2-class Stouffer-Toby solution.\n", res$C[which.min(res$BIC)]))
C loglik npar AIC BIC Gsq df 1 1 -543.65 4 1095.3 1108.8 81.08 11 2 2 -504.47 9 1026.9 1057.3 2.72 6 3 3 -503.30 14 1034.6 1081.9 0.39 1
BIC minimised at C = 2 classes -- the classic 2-class Stouffer-Toby solution.
3. Bayesian LCA by Gibbs — BayesLCA¶
BayesLCA::blca fits the same model in the Bayesian framework. With method="gibbs" it runs the data-augmentation Gibbs sampler — the package analogue of our from-scratch sampler — returning posterior means (and, with the MAP/EM and variational options, alternative fits) for the class probabilities and item-response probabilities. blca wants a 0/1 matrix, so we recode 2→1.
set.seed(3)
X01 <- as.matrix(d) - 1 # 0/1, 1 = particularistic
bg <- blca(X01, 2, method="gibbs", iter=6000, burn=1000, verbose=FALSE)
ordb <- order(bg$classprob, decreasing=TRUE)
Pb <- bg$classprob[ordb]; ipb <- bg$itemprob[ordb, , drop=FALSE]
rownames(ipb) <- paste0("class", 1:2); colnames(ipb) <- items
cat("BayesLCA(Gibbs) class prevalences:", round(Pb,3), "\n\nposterior item-response probabilities P(particularistic | class):\n"); print(round(ipb,3))
BayesLCA(Gibbs) class prevalences: 0.694 0.306 posterior item-response probabilities P(particularistic | class):
A B C D class1 0.708 0.323 0.348 0.126 class2 0.961 0.898 0.889 0.745
# poLCA (EM) vs BayesLCA (Gibbs): the two engines agree on the class profiles
matplot(t(prof), type="b", pch=19, lwd=2, col=c(BLUE,RED), lty=1, ylim=c(0,1), xaxt="n", xlab="item", ylab="P(particularistic)", main="poLCA (solid) vs BayesLCA Gibbs (dashed) -- same profiles")
matlines(t(ipb), type="b", pch=1, lwd=2, col=c(BLUE,RED), lty=2)
axis(1, at=1:4, labels=items)
legend("topright", c("class1 poLCA","class2 poLCA","class1 BayesLCA","class2 BayesLCA"), col=c(BLUE,RED,BLUE,RED), lwd=2, lty=c(1,1,2,2), pch=c(19,19,1,1), bty="n", cex=.8)
cat(sprintf("Max |poLCA - BayesLCA| item-prob difference: %.3f -- EM (poLCA) and Gibbs (BayesLCA) agree.\n", max(abs(prof-ipb))))
Max |poLCA - BayesLCA| item-prob difference: 0.042 -- EM (poLCA) and Gibbs (BayesLCA) agree.
4. Summary¶
poLCA (EM / maximum likelihood) and BayesLCA (Gibbs / Bayesian) fit the same latent class model to the Stouffer–Toby data and reach the same conclusion as the from-scratch and PyMC engines in the Python notebook: a two-class solution (BIC-selected), a large universalistic type and a smaller strongly particularistic type, with matching item-response profiles. A cross-language, cross-implementation corroboration of the foundations.
Next in the arc: choosing the number of classes rigorously on Congdon's larger 17-item dataset, then latent class regression with covariates, and the Hui–Walter diagnostic-testing model.