Choosing the Number of Classes (R)¶
Model selection for LCA with poLCA¶
The R counterpart to modelsel_python.ipynb. We select the number of latent classes for the classic carcinoma data (7 pathologists × 118 slides) using poLCA: the information criteria across models, a bootstrap likelihood-ratio test (implemented with poLCA.simdata), and the fitted 3-class profiles. This cross-checks the Python result — that three latent slide types (clear carcinoma, clear benign, ambiguous) are supported — with the standard R package.
Recall the two obstacles: adding a class always raises the likelihood (so raw likelihood cannot choose $C$), and with $J=7$ items the $2^7=128$-cell table is far too sparse for the classical $\chi^2$ reference distribution of the likelihood-ratio statistic — hence the bootstrap LRT.
options(repr.plot.width=12, repr.plot.height=4.2)
.libPaths(Sys.getenv("R_LIBS_USER")); suppressMessages(library(poLCA))
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; ORANGE<-"#dd6b20"; GREY<-"#718096"; PURP<-"#6b46c1"
d <- read.csv("carcinoma.csv") # raters A..G coded 1/2 (2 = carcinoma call)
raters <- names(d); N <- nrow(d)
cat(N, "slides,", ncol(d), "pathologists:", raters, "\n")
cat("carcinoma-call rate per rater:\n"); print(round(colMeans(d==2),2))
Warning message: "package 'poLCA' was built under R version 4.6.1"
118 slides, 7 pathologists: A B C D E F G
carcinoma-call rate per rater:
A B C D E F G 0.56 0.67 0.38 0.27 0.60 0.21 0.56
1. Information criteria across models¶
We fit $C=1,\dots,5$ with 30 random restarts each and tabulate the log-likelihood, AIC, BIC and the likelihood-ratio $G^2$. BIC (and the small-penalty check via AIC) point to the number of classes; the sparsity caveat means we do not trust the $G^2$ $p$-value, only its trend.
set.seed(1)
f <- cbind(A,B,C,D,E,F,G) ~ 1
res <- data.frame()
fits <- list()
for(C in 1:5){ m <- poLCA(f, data=d, nclass=C, nrep=30, verbose=FALSE); fits[[C]] <- m
res <- rbind(res, data.frame(C=C, loglik=round(m$llik,1), npar=m$npar, AIC=round(m$aic,1), BIC=round(m$bic,1), Gsq=round(m$Gsq,1), df=m$resid.df)) }
print(res)
cat(sprintf("\nAIC min at C=%d, BIC min at C=%d.\n", res$C[which.min(res$AIC)], res$C[which.min(res$BIC)]))
matplot(res$C, cbind(res$AIC-min(res$AIC), res$BIC-min(res$BIC)), type="b", pch=19, lwd=2, col=c(ORANGE,BLUE), lty=1,
xlab="number of classes C", ylab="criterion - min", main="Information criteria vs #classes (carcinoma)")
legend("top", c(paste0("AIC (min C=",which.min(res$AIC),")"), paste0("BIC (min C=",which.min(res$BIC),")")), col=c(ORANGE,BLUE), lwd=2, pch=19, bty="n")
C loglik npar AIC BIC Gsq df 1 1 -524.5 7 1062.9 1082.3 476.8 111 2 2 -317.3 15 664.5 706.1 62.4 103 3 3 -293.7 23 633.4 697.1 15.3 95 4 4 -289.3 31 640.6 726.5 6.4 87 5 5 -286.9 39 651.8 759.8 1.6 79
AIC min at C=3, BIC min at C=3.
2. Bootstrap likelihood-ratio test¶
The rigorous test of "$C$ vs $C+1$ classes". The likelihood-ratio statistic $\text{LR}=2(\hat\ell_{C+1}-\hat\ell_{C})$ does not follow the usual $\chi^2$ here, so we simulate its null distribution: draw many datasets from the fitted $C$-class model (poLCA.simdata), refit both models to each, and locate the observed LR in that bootstrap distribution.
blrt <- function(C, B=40){
m0 <- poLCA(f, data=d, nclass=C, nrep=20, verbose=FALSE)
m1 <- poLCA(f, data=d, nclass=C+1, nrep=20, verbose=FALSE)
LRobs <- 2*(m1$llik - m0$llik)
LRnull <- numeric(B)
for(b in 1:B){
sim <- poLCA.simdata(N=N, probs=m0$probs, P=m0$P)
sd <- as.data.frame(sim$dat); names(sd) <- raters
g0 <- poLCA(f, sd, nclass=C, nrep=5, verbose=FALSE)
g1 <- poLCA(f, sd, nclass=C+1, nrep=5, verbose=FALSE)
LRnull[b] <- 2*(g1$llik - g0$llik)
}
list(LRobs=LRobs, LRnull=LRnull, pval=(1+sum(LRnull>=LRobs))/(B+1))
}
set.seed(2); b23 <- blrt(2, B=40); b34 <- blrt(3, B=40)
par(mfrow=c(1,2), mar=c(4,4,3,1))
for(bb in list(list(b23,"2 vs 3"), list(b34,"3 vs 4"))){ r<-bb[[1]]
hist(r$LRnull, breaks=20, col=GREY, border="white", xlim=range(c(r$LRnull,r$LRobs)), main=sprintf("BLRT %s (p=%.3f)", bb[[2]], r$pval), xlab="LR statistic")
abline(v=r$LRobs, col=RED, lwd=2); legend("topright", sprintf("observed=%.1f", r$LRobs), col=RED, lwd=2, bty="n", cex=.8) }
par(mfrow=c(1,1))
cat(sprintf("C=2 vs 3: p=%.3f (reject -> need a 3rd class). C=3 vs 4: p=%.3f (do not reject -> stop at 3).\n", b23$pval, b34$pval))
C=2 vs 3: p=0.024 (reject -> need a 3rd class). C=3 vs 4: p=0.146 (do not reject -> stop at 3).
3. The chosen 3-class model¶
Every criterion favours three classes. The poLCA profiles give the same substantive reading as the Python notebook: two clear classes and one ambiguous class on which the pathologists disagree.
m3 <- fits[[3]]; ord <- order(m3$P, decreasing=TRUE)
Pc <- m3$P[ord]; prof <- sapply(m3$probs, function(m) m[ord, 2]) # P(carcinoma call) by class x rater
rownames(prof) <- paste0("class", 1:3, " (lambda=", round(Pc,2), ")"); colnames(prof) <- raters
cat("class prevalences:", round(Pc,3), "\n\nP(carcinoma call | class):\n"); print(round(prof,2))
nm <- c("clear carcinoma","clear benign","ambiguous / disputed")
matplot(t(prof), type="b", pch=19, lwd=2, col=c(RED,GREEN,ORANGE), lty=1, ylim=c(0,1), xaxt="n", xlab="pathologist", ylab="P(carcinoma call)", main="Three latent slide types (poLCA)")
axis(1, at=1:7, labels=raters); legend("right", paste0(nm," (",round(Pc,2),")"), col=c(RED,GREEN,ORANGE), lwd=2, pch=19, bty="n", cex=.8)
class prevalences: 0.445 0.374 0.182 P(carcinoma call | class):
A B C D E F G class1 (lambda=0.44) 1.00 0.98 0.86 0.59 1.00 0.48 1.00 class2 (lambda=0.37) 0.06 0.14 0.00 0.00 0.06 0.00 0.00 class3 (lambda=0.18) 0.51 1.00 0.00 0.06 0.75 0.00 0.63
4. Summary¶
poLCA reaches the same conclusion as the from-scratch and PyMC engines: the carcinoma ratings support three latent slide types — clear carcinoma, clear benign, and an ambiguous class — selected by BIC and confirmed by the bootstrap likelihood-ratio test (reject two classes, do not reject four). The full-table $G^2$ was read only for its trend, since the 128-cell table is too sparse for its asymptotic $p$-value.
The Python notebook adds the fully-Bayesian view — WAIC/LOO and the marginal likelihood via Chib's method (with Bayes factors decisively selecting three classes and BIC seen as its large-$N$ approximation). Next in the arc: latent class regression with covariates, and the Hui–Walter diagnostic-testing model.