Diagnostic Testing Without a Gold Standard¶

The Hui–Walter latent class model¶

The diagnostic-testing project of the Latent Class Analysis arc. How accurate is a medical test when there is no gold standard — no way to observe who is truly diseased? You cannot compute sensitivity and specificity the usual way, because the "truth" column is missing. Hui & Walter (1980) solved this by treating the true disease status as a latent class: apply several imperfect tests, and estimate the disease prevalence and every test's accuracy jointly, with the unobserved status integrated out. It is a two-class LCA read in the language of epidemiology, and it turns the carcinoma rater-agreement data of Projects 2–3 into a working example: the seven pathologists become seven imperfect tests of true carcinoma.

The model¶

For subject $i$ in population $g$, with latent disease status $T_i\in\{0,1\}$ and $J$ binary tests: $$T_i=1\ \text{with prob. }\pi_g\ \text{(prevalence)},\qquad x_{ij}\mid T_i=1\sim\text{Bernoulli}(\text{Se}_j),\qquad x_{ij}\mid T_i=0\sim\text{Bernoulli}(1-\text{Sp}_j),$$ where Se$_j$ = sensitivity = $\Pr(\text{test }j\ +\mid\text{diseased})$ and Sp$_j$ = specificity = $\Pr(\text{test }j\ -\mid\text{healthy})$, under conditional independence of the tests given true status. This is exactly a two-class LCA: the classes are diseased and healthy, the item-response probabilities are the Se and $1-$Sp, and the class prevalence is the disease prevalence.

Sensitivity and specificity, precisely. Both are properties of a test, defined by conditioning on the (unobserved) true status — which is exactly why a latent class model can estimate them. Writing the $2\times2$ confusion counts as true-positive TP, false-negative FN, false-positive FP, true-negative TN:

  • Sensitivity $\text{Se}_j$ — the true-positive rate: among subjects who truly have the disease, the fraction test $j$ correctly calls positive, $\text{Se}_j=\Pr(\text{test }+\mid\text{diseased})=\dfrac{\text{TP}}{\text{TP}+\text{FN}}$. A highly sensitive test rarely misses a case, so a negative result helps rule the disease out (mnemonic SnNout).
  • Specificity $\text{Sp}_j$ — the true-negative rate: among subjects who are truly healthy, the fraction test $j$ correctly calls negative, $\text{Sp}_j=\Pr(\text{test }-\mid\text{healthy})=\dfrac{\text{TN}}{\text{TN}+\text{FP}}$. A highly specific test rarely raises a false alarm, so a positive result helps rule the disease in (mnemonic SpPin).

Because both condition on the truth, they are intrinsic to the test and do not depend on prevalence — unlike the predictive values (PPV $=\Pr(\text{diseased}\mid+)$, NPV $=\Pr(\text{healthy}\mid-)$), which do. The two quantities trade off as a reader's calling threshold shifts — the ROC curve of §4: an aggressive reader buys sensitivity at the cost of specificity, a conservative one the reverse. (A perfect test has $\text{Se}=\text{Sp}=1$; a coin flip has $\text{Se}+\text{Sp}=1$, the ROC diagonal.)

The catch: identifiability¶

The unobserved truth is not free. Counting degrees of freedom:

  • 2 tests, 1 population — the $2\times2$ table has only 3 free cells, but the model has 5 parameters ($\pi,\text{Se}_1,\text{Se}_2,\text{Sp}_1,\text{Sp}_2$). Under-identified.
  • Hui & Walter's fix: use two populations of different prevalence sharing the same test characteristics — $2\times3=6$ df for $6$ parameters. Identified.
  • Or use $\ge3$ tests in one population — $2^3-1=7$ df for $7$ parameters. Identified.
  • Or supply informative priors on Se/Sp — the Bayesian route (Joseph, Gyorkos & Coupal 1995), which regularises even the under-identified design.

We fit by a conjugate data-augmentation Gibbs sampler (impute the disease status; then $\pi_g$, Se, Sp are Beta-conjugate), and demonstrate each of these identifiability facts directly.

Sections¶

  1. Identifiability — when can the truth be recovered?
  2. Recovery on identified designs
  3. Priors to the rescue — the under-identified design
  4. Carcinoma — seven pathologists as seven tests
  5. Summary
In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import huiwalter as H, warnings; warnings.filterwarnings("ignore")

plt.rcParams.update({"figure.figsize": (9,4), "axes.grid": True, "grid.alpha": .22,
                     "axes.spines.top": False, "axes.spines.right": False, "font.size": 10.5})
BLUE, RED, GREEN, ORANGE, GREY, PURP = "#2b6cb0","#c53030","#2f855a","#dd6b20","#718096","#6b46c1"
rng = np.random.default_rng(0)
Se_t = np.array([0.90, 0.80, 0.75]); Sp_t = np.array([0.95, 0.85, 0.90])   # true test accuracies
print("ready")

ready

1. Identifiability — when can the truth be recovered?¶

The single most important fact about diagnostic-test LCA is that more parameters than data cannot be estimated, however clever the sampler. We show it three ways from the same true tests. With two tests in one population the posterior fails to concentrate at the truth — it drifts, because the design is under-identified. Adding a third test, or a second population, supplies the missing degrees of freedom and the posterior locks onto the true sensitivities.

In [2]:
# same true tests, three designs
X1,p1,_ = H.simulate_hw([1500], [0.30], Se_t[:2], Sp_t[:2], rng)          # 2 tests, 1 pop (under-identified)
X3,p3,_ = H.simulate_hw([1500], [0.30], Se_t,     Sp_t,     rng)          # 3 tests, 1 pop (identified)
X2,p2,_ = H.simulate_hw([1200,1200], [0.20,0.55], Se_t[:2], Sp_t[:2], rng) # 2 tests, 2 pops (identified)
post1 = H.hw_gibbs(X1, p1, rng, draws=3000, burn=1000)
post3 = H.hw_gibbs(X3, p3, rng, draws=3000, burn=1000)
post2 = H.hw_gibbs(X2, p2, rng, draws=3000, burn=1000)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.3))
# posterior of Se_1 under each design vs the truth
for lab,pst,c in [("2 tests, 1 pop (under-identified)",post1,RED),("3 tests, 1 pop",post3,GREEN),("2 tests, 2 pops",post2,BLUE)]:
    ax[0].hist(pst["Se"][:,0], bins=40, density=True, histtype="step", lw=2, color=c, label=lab)
ax[0].axvline(Se_t[0], color="k", ls="--", lw=1.5, label=f"true Se₁={Se_t[0]}")
ax[0].set_title("Posterior of test-1 sensitivity by design"); ax[0].set_xlabel("Se₁"); ax[0].legend(fontsize=8); ax[0].set_xlim(0.4,1)
# the under-identified ridge: Se1 vs Sp1 draws wander; the identified design is a tight blob at truth
ax[1].scatter(post1["Se"][:,0], post1["Sp"][:,0], s=4, alpha=.15, color=RED, label="2 tests, 1 pop")
ax[1].scatter(post3["Se"][:,0], post3["Sp"][:,0], s=4, alpha=.20, color=GREEN, label="3 tests, 1 pop")
ax[1].plot(Se_t[0], Sp_t[0], "k*", ms=16, label="truth")
ax[1].set_title("Joint posterior (Se₁, Sp₁): ridge vs tight blob"); ax[1].set_xlabel("Se₁"); ax[1].set_ylabel("Sp₁"); ax[1].legend(fontsize=8); ax[1].set_xlim(0.4,1); ax[1].set_ylim(0.7,1)
plt.tight_layout(); plt.show()
print(f"Under-identified (2 tests,1 pop): Se₁ posterior mean {post1['Se'][:,0].mean():.2f} ± {post1['Se'][:,0].std():.2f} — wide and biased from {Se_t[0]}.")
print(f"Identified (3 tests): {post3['Se'][:,0].mean():.2f} ± {post3['Se'][:,0].std():.2f}. Identified (2 pops): {post2['Se'][:,0].mean():.2f} ± {post2['Se'][:,0].std():.2f}. Both lock onto the truth.")
No description has been provided for this image 
Under-identified (2 tests,1 pop): Se₁ posterior mean 0.84 ± 0.09 — wide and biased from 0.9.
Identified (3 tests): 0.91 ± 0.03. Identified (2 pops): 0.88 ± 0.03. Both lock onto the truth.

2. Recovery on an identified design¶

With an identifiable design, the Hui–Walter model recovers the whole panel of test accuracies and the prevalence — all without ever observing who is diseased. Below, the three-test/one-population design: posterior means and 95% credible intervals against the known truth.

In [3]:
post = post3
tab = H.summary(post)
print(tab.round(3).to_string())
truth = {"prevalence[pop1]": 0.30, **{f"Se[test{j+1}]": Se_t[j] for j in range(3)}, **{f"Sp[test{j+1}]": Sp_t[j] for j in range(3)}}
fig, ax = plt.subplots(figsize=(9.5, 4.6)); y = np.arange(len(tab))
ax.errorbar(tab["mean"], y, xerr=[tab["mean"]-tab["lo95"], tab["hi95"]-tab["mean"]], fmt="o", color=BLUE, capsize=3, ms=7, label="posterior mean ±95% CI")
ax.plot([truth[p] for p in tab.index], y, "D", color=RED, ms=8, label="true value")
ax.set_yticks(y); ax.set_yticklabels(tab.index); ax.invert_yaxis(); ax.set_xlim(0,1)
ax.set_title("Hui–Walter recovers prevalence, sensitivity and specificity"); ax.set_xlabel("probability"); ax.legend(fontsize=8.5)
plt.tight_layout(); plt.show()
print("Every true value (red) falls inside its 95% credible interval — the latent-class model reconstructs the missing 'truth'.")

                   mean   lo95   hi95
parameter                            
prevalence[pop1]  0.290  0.260  0.322
Se[test1]         0.905  0.854  0.954
Se[test2]         0.810  0.759  0.856
Se[test3]         0.744  0.692  0.794
Sp[test1]         0.958  0.936  0.977
Sp[test2]         0.863  0.838  0.887
Sp[test3]         0.899  0.877  0.921
No description has been provided for this image 
Every true value (red) falls inside its 95% credible interval — the latent-class model reconstructs the missing 'truth'.

3. Priors to the rescue — the under-identified design¶

Often only two tests in one population are available, which we saw is not identified. The Bayesian remedy (Joseph, Gyorkos & Coupal, 1995) is to bring external knowledge — an informative $\text{Beta}$ prior on each Se and Sp, from the literature or a manufacturer's insert. The prior supplies the degrees of freedom the data lack, turning a hopeless likelihood into a usable posterior. It is not magic: the result leans on the prior, so honest, defensible priors matter.

In [4]:
Xu, pu, _ = H.simulate_hw([1500], [0.30], Se_t[:2], Sp_t[:2], rng)     # under-identified design
vague = H.hw_gibbs(Xu, pu, rng, se_ab=(1,1), sp_ab=(1,1), draws=3000, burn=1000)
info  = H.hw_gibbs(Xu, pu, rng, se_ab=(8,2), sp_ab=(8,2), draws=3000, burn=1000)   # Beta(8,2): mean 0.8
fig, ax = plt.subplots(1, 2, figsize=(13, 4))
for k,(lab,c) in enumerate([("Se₁",BLUE),("Se₂",GREEN)]):
    ax[0].hist(vague["Se"][:,k], bins=40, density=True, histtype="step", lw=2, color=c, label=f"{lab} vague prior")
    ax[0].hist(info["Se"][:,k], bins=40, density=True, histtype="stepfilled", lw=0, alpha=.25, color=c, label=f"{lab} informative prior")
    ax[0].axvline(Se_t[k], color=c, ls="--", lw=1)
ax[0].set_title("Sensitivities: vague vs informative priors (2 tests, 1 pop)"); ax[0].set_xlabel("sensitivity"); ax[0].legend(fontsize=7.5); ax[0].set_xlim(0.4,1)
xx = np.linspace(0,1,200)
from scipy.stats import beta as Bdist
ax[1].plot(xx, Bdist.pdf(xx,1,1), color=GREY, lw=2, label="vague  Beta(1,1)")
ax[1].plot(xx, Bdist.pdf(xx,8,2), color=PURP, lw=2, label="informative  Beta(8,2)")
ax[1].set_title("The two prior choices for Se/Sp"); ax[1].set_xlabel("probability"); ax[1].legend(fontsize=8.5)
plt.tight_layout(); plt.show()
sd_v, sd_i = vague["Se"].std(0).mean(), info["Se"].std(0).mean()
print(f"Vague Beta(1,1):       Se posterior means {np.round(vague['Se'].mean(0),3)}, mean sd {sd_v:.3f}")
_chg = 100 * (1 - sd_i / sd_v)
print(f"Informative Beta(8,2): Se posterior means {np.round(info['Se'].mean(0),3)}, mean sd {sd_i:.3f}"
      f"  ({abs(_chg):.0f}% {'tighter' if _chg > 0 else 'wider'})")
print(f"The true sensitivities are {Se_t[:2]} and the prior mean is 0.8. With only two tests in one")
print("population the likelihood cannot separate sensitivity from specificity, so the posterior location")
print("is set by the prior rather than by the data — which is what 'under-identified' means in practice.")
No description has been provided for this image 
Vague Beta(1,1):       Se posterior means [0.822 0.862], mean sd 0.089
Informative Beta(8,2): Se posterior means [0.81 0.86], mean sd 0.076  (15% tighter)
The true sensitivities are [0.9 0.8] and the prior mean is 0.8. With only two tests in one
population the likelihood cannot separate sensitivity from specificity, so the posterior location
is set by the prior rather than by the data — which is what 'under-identified' means in practice.

4. Carcinoma — seven pathologists as seven tests¶

Return to the carcinoma data: 118 slides, each read by seven pathologists for carcinoma. With no gold standard, treat each pathologist as an imperfect binary test and let the Hui–Walter model estimate the disease prevalence among the slides and each reader's sensitivity and specificity. Seven tests in one population is comfortably identified. The result is a data-driven accuracy audit of the raters — some are aggressive callers (high Se, low Sp), some conservative (low Se, high Sp).

In [5]:
df = pd.read_csv("carcinoma.csv"); Xc = (df.values - 1).astype(int); raters = list(df.columns)
pc = H.hw_gibbs(Xc, np.zeros(len(Xc), int), rng, draws=5000, burn=1500)
Se = pc["Se"].mean(0); Sp = pc["Sp"].mean(0); Jc = H.youden(pc).mean(0)
res = pd.DataFrame({"Se": Se.round(2), "Sp": Sp.round(2), "Youden J": Jc.round(2)}, index=raters)
print(f"estimated carcinoma prevalence among slides: {pc['prev'].mean():.2f}  [{np.percentile(pc['prev'],2.5):.2f}, {np.percentile(pc['prev'],97.5):.2f}]\n")
print(res.to_string())
fig, ax = plt.subplots(1, 2, figsize=(13.5, 4.4))
# forest plot of Se and Sp with CIs
y = np.arange(len(raters))
ax[0].errorbar(Se, y+0.12, xerr=[Se-np.percentile(pc["Se"],2.5,0), np.percentile(pc["Se"],97.5,0)-Se], fmt="o", color=BLUE, capsize=2, label="sensitivity")
ax[0].errorbar(Sp, y-0.12, xerr=[Sp-np.percentile(pc["Sp"],2.5,0), np.percentile(pc["Sp"],97.5,0)-Sp], fmt="s", color=RED, capsize=2, label="specificity")
ax[0].set_yticks(y); ax[0].set_yticklabels(raters); ax[0].invert_yaxis(); ax[0].set_xlim(0,1)
ax[0].set_title("Per-pathologist sensitivity & specificity (±95% CI)"); ax[0].set_xlabel("probability"); ax[0].legend(fontsize=8.5)
# ROC space: each rater at (1-Sp, Se); distance above the diagonal = Youden's J
ax[1].plot([0,1],[0,1],"--",color=GREY,lw=1)
ax[1].scatter(1-Sp, Se, s=90, color=PURP, zorder=3)
for j,r in enumerate(raters): ax[1].annotate(r, (1-Sp[j], Se[j]), fontsize=9, ha="left", va="bottom")
ax[1].set_xlim(0,1); ax[1].set_ylim(0,1); ax[1].set_xlabel("1 − specificity (false-positive rate)"); ax[1].set_ylabel("sensitivity (true-positive rate)")
ax[1].set_title("Pathologists in ROC space"); ax[1].set_aspect("equal")
plt.tight_layout(); plt.show()
print("About half the slides carry 'true' carcinoma. Readers differ sharply: some (top-left of ROC) are accurate; aggressive")
print("callers sit toward the right (high Se, more false positives), conservative ones toward the bottom (high Sp, misses).")

estimated carcinoma prevalence among slides: 0.53  [0.43, 0.62]

     Se    Sp  Youden J
A  0.95  0.88      0.82
B  0.97  0.68      0.65
C  0.72  0.98      0.70
D  0.51  0.98      0.50
E  0.96  0.81      0.77
F  0.40  0.98      0.39
G  0.98  0.91      0.89
No description has been provided for this image 
About half the slides carry 'true' carcinoma. Readers differ sharply: some (top-left of ROC) are accurate; aggressive
callers sit toward the right (high Se, more false positives), conservative ones toward the bottom (high Sp, misses).

A caveat carried from Project 2. The Hui–Walter model assumes the tests are conditionally independent given true status. But Project 2 found that these ratings are best described by three latent classes — a sign that, beyond "carcinoma vs benign", the pathologists share a residual ambiguous structure they agree on, i.e. conditional dependence. When tests are correlated given the truth, the two-class Se/Sp estimates are biased (typically Se and Sp are over-stated). The principled fixes are to add a dependence term or a shared random effect between correlated readers (the Dendukuri–Joseph and fixed-/random-effects extensions) — the diagnostic-testing analogue of relaxing local independence in LCA.

5. Summary¶

The Hui–Walter model estimates test accuracy when the truth is never observed, by treating true status as a latent class and fitting a two-class LCA in the language of prevalence, sensitivity and specificity. Its defining lesson is identifiability: two tests in one population cannot be estimated, but a second population (Hui & Walter's original device) or a third test supplies the missing degrees of freedom — both demonstrated here, and both recovering the known truth on simulated data. Informative priors make the under-identified design estimable but not identified: the posterior tightens and stops drifting, yet its location is set by the prior rather than the data, and it does not return to the truth.

Applied to the carcinoma ratings, it audited seven pathologists without a gold standard, estimating the slide prevalence and each reader's sensitivity, specificity and Youden index — placing them in ROC space and separating the accurate from the aggressive and the conservative. The conditional-independence assumption is the model's soft spot, and the three-class structure of Project 2 is exactly the warning sign that a dependence-aware extension may be needed.

This completes a four-project tour of latent class analysis — foundations, choosing the number of classes, the nonparametric (Dirichlet-process) view, and diagnostic testing — all built on the same conjugate data-augmentation core. A natural next step is latent class regression, where covariates predict class membership, and its longitudinal cousin, latent transition analysis.