Diagnostic Testing Without a Gold Standard (R)¶
The Hui–Walter latent class model, from scratch¶
The R counterpart to huiwalter_python.ipynb. We implement the Hui–Walter data-augmentation Gibbs sampler from scratch in base R — the same algorithm as the Python engine — and confirm it recovers test accuracies on simulated data and audits the seven carcinoma pathologists as seven imperfect tests.
Recall the model: latent disease status $T_i\sim\text{Bernoulli}(\pi_g)$; $x_{ij}\mid T_i=1\sim\text{Bernoulli}(\text{Se}_j)$; $x_{ij}\mid T_i=0\sim\text{Bernoulli}(1-\text{Sp}_j)$, under conditional independence. Here sensitivity $\text{Se}_j=\Pr(\text{test }+\mid\text{diseased})$ is the true-positive rate (a sensitive test rarely misses a case — a negative result rules disease out) and specificity $\text{Sp}_j=\Pr(\text{test }-\mid\text{healthy})$ is the true-negative rate (a specific test rarely gives a false alarm — a positive result rules disease in). Both condition on the true status, so they are intrinsic to the test and do not depend on prevalence (unlike predictive values). The identifiability rule: two tests in one population is under-identified; a second population (Hui & Walter's device) or a third test identifies it. Self-contained; base R only. The standard applied tool in R is runjags::template.huiwalter(), which auto-generates a JAGS Hui–Walter model (any number of tests/populations, with optional conditional-dependence terms); we note it at the end.
options(repr.plot.width=12, repr.plot.height=4.4)
BLUE<-"#2b6cb0"; RED<-"#c53030"; GREEN<-"#2f855a"; ORANGE<-"#dd6b20"; GREY<-"#718096"; PURP<-"#6b46c1"
hw_gibbs <- function(X, pop, se_ab=c(1,1), sp_ab=c(1,1), prev_ab=c(1,1), draws=4000, burn=1500, anchor=TRUE){
N<-nrow(X); J<-ncol(X); G<-max(pop)
Se<-runif(J,0.5,0.9); Sp<-runif(J,0.5,0.9); prev<-runif(G,0.2,0.6)
PREV<-matrix(0,draws,G); SE<-matrix(0,draws,J); SP<-matrix(0,draws,J)
for(t in 1:(draws+burn)){
l1<-log(prev[pop]) + X%*%log(Se) + (1-X)%*%log(1-Se)
l0<-log(1-prev[pop]) + X%*%log(1-Sp) + (1-X)%*%log(Sp)
p1<-1/(1+exp(pmin(pmax(l0-l1,-700),700))); Tt<-as.integer(runif(N)<p1)
dis<-Tt==1; hea<-!dis
for(g in 1:G){ m<-pop==g; nd<-sum(Tt[m]); prev[g]<-rbeta(1, prev_ab[1]+nd, prev_ab[2]+sum(m)-nd) }
sd<-colSums(X[dis,,drop=FALSE]); nd<-sum(dis); Se<-rbeta(J, se_ab[1]+sd, se_ab[2]+nd-sd)
sh<-colSums(X[hea,,drop=FALSE]); nh<-sum(hea); Sp<-rbeta(J, sp_ab[1]+(nh-sh), sp_ab[2]+sh)
if(anchor && (mean(Se)+mean(Sp)<1)){ s2<-1-Sp; Sp<-1-Se; Se<-s2; prev<-1-prev }
if(t>burn){ i<-t-burn; PREV[i,]<-prev; SE[i,]<-Se; SP[i,]<-Sp }
}
list(prev=PREV, Se=SE, Sp=SP)
}
sim_hw <- function(n_per_pop, prev, Se, Sp){ G<-length(prev); J<-length(Se); if(length(n_per_pop)==1) n_per_pop<-rep(n_per_pop,G)
X<-NULL; pop<-NULL; for(g in 1:G){ n<-n_per_pop[g]; Tt<-as.integer(runif(n)<prev[g])
ppos<-ifelse(matrix(Tt,n,J)==1, matrix(Se,n,J,byrow=TRUE), matrix(1-Sp,n,J,byrow=TRUE))
X<-rbind(X, matrix(as.integer(runif(n*J)<ppos),n,J)); pop<-c(pop,rep(g,n)) }
list(X=X, pop=pop) }
ci <- function(m) c(mean(m), quantile(m, c(.025,.975)))
cat("from-scratch R Hui-Walter Gibbs ready\n")
from-scratch R Hui-Walter Gibbs ready
1. Recovery on an identified design¶
Three tests in one population is identified. With known true accuracies ($Se=0.90,0.80,0.75$; $Sp=0.95,0.85,0.90$; prevalence $0.30$), the sampler should recover them all — without ever seeing the true status.
set.seed(1)
Se_t<-c(0.90,0.80,0.75); Sp_t<-c(0.95,0.85,0.90)
s<-sim_hw(1500, 0.30, Se_t, Sp_t); post<-hw_gibbs(s$X, s$pop, draws=3000, burn=1000)
est<-rbind(prevalence=ci(post$prev[,1]),
setNames(as.data.frame(t(sapply(1:3, function(j) ci(post$Se[,j])))), c("mean","lo","hi")),
setNames(as.data.frame(t(sapply(1:3, function(j) ci(post$Sp[,j])))), c("mean","lo","hi")))
lab<-c("prevalence", paste0("Se_test",1:3), paste0("Sp_test",1:3)); tru<-c(0.30, Se_t, Sp_t)
E<-t(sapply(list(post$prev[,1], post$Se[,1],post$Se[,2],post$Se[,3], post$Sp[,1],post$Sp[,2],post$Sp[,3]), ci))
rownames(E)<-lab; colnames(E)<-c("mean","lo95","hi95"); print(round(E,3))
y<-length(lab):1
plot(E[,1], y, xlim=c(0,1), pch=19, col=BLUE, yaxt="n", ylab="", xlab="probability", main="Hui-Walter recovers prevalence, sensitivity, specificity")
arrows(E[,2], y, E[,3], y, angle=90, code=3, length=.03, col=BLUE); points(tru, y, pch=18, col=RED, cex=1.4)
axis(2, at=y, labels=lab, las=1, cex.axis=.8); legend("bottomright", c("posterior mean +/-95% CI","true value"), col=c(BLUE,RED), pch=c(19,18), bty="n", cex=.8)
cat("Every true value (red) lies within its 95% credible interval -- the latent-class model reconstructs the missing truth.\n")
mean lo95 hi95 prevalence 0.294 0.263 0.331 Se_test1 0.927 0.871 0.983 Se_test2 0.753 0.701 0.804 Se_test3 0.730 0.672 0.784 Sp_test1 0.950 0.924 0.977 Sp_test2 0.822 0.795 0.848 Sp_test3 0.906 0.883 0.929
Every true value (red) lies within its 95% credible interval -- the latent-class model reconstructs the missing truth.
2. Carcinoma — seven pathologists as seven tests¶
The 118 carcinoma slides read by seven pathologists: treat each reader as an imperfect test and estimate the slide prevalence and every reader's sensitivity and specificity (seven tests in one population is identified).
d<-read.csv("carcinoma.csv"); Xc<-as.matrix(d)-1; raters<-names(d)
set.seed(2); pc<-hw_gibbs(Xc, rep(1,nrow(Xc)), draws=5000, burn=1500)
Se<-colMeans(pc$Se); Sp<-colMeans(pc$Sp); J<-Se+Sp-1
res<-data.frame(Se=round(Se,2), Sp=round(Sp,2), Youden=round(J,2), row.names=raters)
cat(sprintf("estimated carcinoma prevalence: %.2f [%.2f, %.2f]\n\n", mean(pc$prev), quantile(pc$prev,.025), quantile(pc$prev,.975)))
print(res)
par(mfrow=c(1,2), mar=c(4,4,3,1)); y<-7:1
plot(Se, y+0.12, xlim=c(0,1), pch=19, col=BLUE, yaxt="n", ylab="", xlab="probability", main="Per-pathologist Se & Sp (+/-95% CI)")
arrows(apply(pc$Se,2,quantile,.025), y+0.12, apply(pc$Se,2,quantile,.975), y+0.12, angle=90, code=3, length=.02, col=BLUE)
points(Sp, y-0.12, pch=15, col=RED); arrows(apply(pc$Sp,2,quantile,.025), y-0.12, apply(pc$Sp,2,quantile,.975), y-0.12, angle=90, code=3, length=.02, col=RED)
axis(2, at=y, labels=raters, las=1); legend("bottomleft", c("sensitivity","specificity"), col=c(BLUE,RED), pch=c(19,15), bty="n", cex=.8)
plot(1-Sp, Se, xlim=c(0,1), ylim=c(0,1), pch=19, col=PURP, cex=1.6, xlab="1 - specificity", ylab="sensitivity", main="Pathologists in ROC space", asp=1)
abline(0,1,lty=2,col=GREY); text(1-Sp, Se, raters, pos=4, cex=.9)
par(mfrow=c(1,1))
cat("About half the slides carry 'true' carcinoma; readers range from accurate (top-left) to aggressive (high Se, more false positives)\n")
cat("or conservative (high Sp, misses). NOTE: Hui-Walter assumes conditional independence; Project 2's third latent class signals\n")
cat("residual dependence among readers, which can bias these Se/Sp -- a dependence-aware extension (e.g. runjags) would address it.\n")
estimated carcinoma prevalence: 0.53 [0.43, 0.62]
Se Sp Youden A 0.95 0.88 0.83 B 0.97 0.68 0.65 C 0.71 0.98 0.70 D 0.51 0.98 0.49 E 0.96 0.81 0.77 F 0.40 0.98 0.39 G 0.98 0.91 0.89
About half the slides carry 'true' carcinoma; readers range from accurate (top-left) to aggressive (high Se, more false positives)
or conservative (high Sp, misses). NOTE: Hui-Walter assumes conditional independence; Project 2's third latent class signals
residual dependence among readers, which can bias these Se/Sp -- a dependence-aware extension (e.g. runjags) would address it.
3. Summary¶
The from-scratch R Hui–Walter Gibbs sampler reproduces the Python results: on an identified design it recovers prevalence, sensitivity and specificity without a gold standard, and on the carcinoma data it audits the seven pathologists, placing them in ROC space. A cross-language corroboration of the diagnostic-testing latent class model.
For applied work, runjags::template.huiwalter() (Matt Denwood) is the standard R tool: given a data frame of test results and a population column, it auto-writes a JAGS Hui–Walter model — handling any number of tests and populations, informative Se/Sp priors, and conditional-dependence terms between correlated tests, the very extension flagged by the three-class structure of Project 2. This completes the latent class arc's diagnostic-testing project; a natural next step is latent class regression with covariates.