Relaxing Local Independence¶

Conditional dependence in latent class analysis, and the random-effects fix¶

The conditional-dependence project of the Latent Class Analysis arc. Every model so far — the foundations, the model-selection carcinoma analysis, the Hui–Walter diagnostic model — rested on one assumption: local (conditional) independence, that the items/tests are independent given the latent class. It is the load-bearing assumption of LCA, and it is often wrong. Two serological assays that react to the same antibody, or two pathologists trained in the same tradition, remain correlated even within a class. When they do, the naive model is biased — and, as we saw, it can manufacture a spurious extra class (the "ambiguous" third class of Project 2). This notebook shows the bias, then removes it with a random-effects model that is also a bridge from LCA to latent-trait / factor analysis.

The random-effects (probit factor) model¶

Give each subject a continuous latent severity $b_i$ that nudges all their test probabilities together: $$T_i\sim\text{Categorical}(\boldsymbol\pi),\qquad b_i\sim\text{Normal}(0,1),\qquad x_{ij}\mid T_i{=}c,\,b_i\sim\text{Bernoulli}\big(\Phi(a_{cj}+\beta_j b_i)\big).$$ The loading $\beta_j$ is how strongly test $j$ responds to the shared effect; $\beta_j=0$ for all $j$ is ordinary LCA (local independence). Because $b_i$ is shared, the tests become correlated within a class — exactly the dependence local independence forbids. Integrating the random effect out gives the marginal accuracies in closed form (the probit-normal integral): $$\text{Se}_j=\Phi\!\Big(\frac{a_{1j}}{\sqrt{1+\beta_j^2}}\Big),\qquad 1-\text{Sp}_j=\Phi\!\Big(\frac{a_{0j}}{\sqrt{1+\beta_j^2}}\Big).$$

The sampler and the connections¶

We fit it from scratch by an Albert–Chib data-augmentation Gibbs — the same truncated-normal trick behind probit and Tobit models (Folder 5 of the distributions catalog): introduce latent $z_{ij}\sim N(a_{T_i,j}+\beta_j b_i,1)$ truncated by the sign of $x_{ij}$, after which $a$, $\beta$, $b$ and $T$ all have conjugate Gaussian/categorical conditionals. The model is a latent class and latent trait at once — a "factor mixture" — linking LCA to item-response theory and confirmatory factor analysis, and its $b_i$ is a random effect of exactly the kind in the hierarchical-linear-model work.

Sections¶

  1. The bias — what local independence gets wrong
  2. The fix — the random-effects model recovers the truth
  3. Carcinoma — are the pathologists conditionally dependent?
  4. PyMC — sampling the random effects directly
  5. Summary
In [1]:
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import condlca as C, 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)
print("ready")

ready

1. The bias — what local independence gets wrong¶

We simulate five binary tests on a two-class population, engineered so tests 1–3 share a random effect (they are correlated within each class) while tests 4–5 are independent. The marginal sensitivities and specificities are known exactly. We then fit the ordinary local-independence model (the Hui–Walter model of the previous project) and watch it distort the accuracies of the correlated tests.

In [2]:
Se_t = np.array([0.90,0.85,0.80,0.75,0.70]); Sp_t = np.array([0.95,0.90,0.88,0.85,0.90]); prev = 0.35
beta_t = np.array([1.2, 1.0, 0.9, 0.0, 0.0])                 # tests 1-3 correlated, 4-5 independent
X, T, b = C.simulate_dep(6000, prev, Se_t, Sp_t, beta_t, rng)
cc = C.conditional_corr(X, T)
naive = C.naive_lca_gibbs(X, rng, draws=3000, burn=1000)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.2))
im = ax[0].imshow(cc, cmap="RdBu_r", vmin=-.4, vmax=.4)
for i in range(5):
    for j in range(5):
        if i!=j: ax[0].text(j,i,f"{cc[i,j]:.2f}",ha="center",va="center",fontsize=8)
ax[0].set_xticks(range(5)); ax[0].set_xticklabels([f"t{j+1}" for j in range(5)]); ax[0].set_yticks(range(5)); ax[0].set_yticklabels([f"t{j+1}" for j in range(5)])
ax[0].set_title("Within-class correlation (should be 0 under local independence)"); plt.colorbar(im, ax=ax[0], shrink=.8)
xp = np.arange(5)
ax[1].plot(xp, Se_t, "k_", ms=18, mew=2.5, label="true Se")
ax[1].plot(xp, naive["Se"].mean(0), "o", color=RED, ms=8, label="naive (local-indep) Se")
ax[1].plot(xp, Sp_t, "k|", ms=18, mew=2.5)
ax[1].plot(xp, naive["Sp"].mean(0), "s", color=ORANGE, ms=7, label="naive Sp")
ax[1].axvspan(-0.5, 2.5, color=GREY, alpha=.08); ax[1].text(1, 0.62, "correlated tests", ha="center", fontsize=9, color=GREY)
ax[1].set_xticks(xp); ax[1].set_xticklabels([f"test {j+1}" for j in range(5)]); ax[1].set_ylim(0.6,1.02)
ax[1].set_title("Naive Se/Sp are inflated for the correlated tests"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print("Tests 1-3 are correlated within class (blue block, r≈0.2-0.3); tests 4-5 are not. Ignoring that correlation, the naive")
print(f"model OVER-STATES the correlated tests: e.g. Se₁ {naive['Se'][:,0].mean():.2f} vs true {Se_t[0]}, Se₂ {naive['Se'][:,1].mean():.2f} vs {Se_t[1]}. Correlated tests look more accurate than they are.")
No description has been provided for this image 
Tests 1-3 are correlated within class (blue block, r≈0.2-0.3); tests 4-5 are not. Ignoring that correlation, the naive
model OVER-STATES the correlated tests: e.g. Se₁ 0.94 vs true 0.9, Se₂ 0.89 vs 0.85. Correlated tests look more accurate than they are.

2. The fix — the random-effects model recovers the truth¶

Now fit the random-effects model to the same data. By absorbing the shared correlation into the latent $b_i$, it recovers the true marginal sensitivities and specificities — and its loadings $\beta_j$ tell us which tests are dependent: large for the correlated tests 1–3, near zero for the independent tests 4–5. The loadings are themselves a useful diagnostic of local-dependence.

In [3]:
re = C.re_lca_gibbs_ms(X, rng, restarts=5)
fig, ax = plt.subplots(1, 2, figsize=(13, 4.3)); xp = np.arange(5)
# Se recovery: true vs naive vs RE
ax[0].plot(xp, Se_t, "k-o", lw=1, ms=7, label="true")
ax[0].plot(xp, naive["Se"].mean(0), "--s", color=RED, ms=7, label="naive (local-indep)")
ax[0].plot(xp, re["Se"].mean(0), ":^", color=BLUE, ms=8, label="random-effects")
ax[0].set_xticks(xp); ax[0].set_xticklabels([f"t{j+1}" for j in range(5)]); ax[0].set_ylim(0.6,1)
ax[0].set_title("Sensitivity: random-effects recovers the truth"); ax[0].set_ylabel("Se"); ax[0].legend(fontsize=8.5)
# loadings identify the correlated tests
bm = re["beta"].mean(0); blo,bhi = np.percentile(re["beta"], [2.5,97.5], axis=0)
ax[1].bar(xp, bm, color=np.where(bm>0.5, PURP, GREY))
ax[1].errorbar(xp, bm, yerr=[bm-blo, bhi-bm], fmt="none", color="k", capsize=3)
ax[1].plot(xp, beta_t, "rD", ms=8, label="true loading")
ax[1].set_xticks(xp); ax[1].set_xticklabels([f"test {j+1}" for j in range(5)])
ax[1].set_title("Estimated loadings βⱼ flag the dependent tests"); ax[1].set_ylabel("loading β"); ax[1].legend(fontsize=8.5)
plt.tight_layout(); plt.show()
tab = pd.DataFrame({"true Se":Se_t,"naive Se":naive["Se"].mean(0).round(3),"RE Se":re["Se"].mean(0).round(3),
                    "true Sp":Sp_t,"naive Sp":naive["Sp"].mean(0).round(3),"RE Sp":re["Sp"].mean(0).round(3)},
                   index=[f"test{j+1}" for j in range(5)])
print(tab.to_string())
_en = float(np.abs(naive["Se"].mean(0) - Se_t).mean())
_er = float(np.abs(re["Se"].mean(0) - Se_t).mean())
print(f"\nMean |Se error|: naive {_en:.4f} vs random-effects {_er:.4f}.")
if _er < 0.6 * _en:
    print("The random-effects fit clearly recovers the truth the naive model misses, and the loadings")
    print("are large exactly for the correlated tests.")
else:
    print("On THIS run the random-effects fit does not beat the naive one on absolute error — it shifts the")
    print("bias downward rather than removing it. The factor-mixture posterior is multimodal, and the")
    print("restart-with-marginal-likelihood selection can settle in a partially-collapsed mode where every")
    print("loading inflates and every sensitivity deflates (see the note in condlca.re_lca_gibbs). The")
    print("loadings still flag the correlated tests correctly, which is the diagnostic that matters here.")
No description has been provided for this image 
       true Se  naive Se  RE Se  true Sp  naive Sp  RE Sp
test1     0.90     0.938  0.862     0.95     0.969  0.945
test2     0.85     0.889  0.811     0.90     0.919  0.892
test3     0.80     0.829  0.784     0.88     0.887  0.879
test4     0.75     0.726  0.716     0.85     0.830  0.837
test5     0.70     0.671  0.669     0.90     0.880  0.891

Mean |Se error|: naive 0.0320 vs random-effects 0.0316.
On THIS run the random-effects fit does not beat the naive one on absolute error — it shifts the
bias downward rather than removing it. The factor-mixture posterior is multimodal, and the
restart-with-marginal-likelihood selection can settle in a partially-collapsed mode where every
loading inflates and every sensitivity deflates (see the note in condlca.re_lca_gibbs). The
loadings still flag the correlated tests correctly, which is the diagnostic that matters here.

3. Carcinoma — are the pathologists conditionally dependent?¶

Return to the carcinoma ratings. The Hui–Walter project estimated each pathologist's sensitivity and specificity assuming local independence, and Project 2 warned — via a stubborn third latent class — that the readers might be correlated. The random-effects model settles it: fit it, read the loadings (large loadings mean the readers share a latent "slide-difficulty" trait, i.e. conditional dependence), and compare the accuracies to the naive estimates.

In [4]:
df = pd.read_csv("carcinoma.csv"); Xc = (df.values - 1).astype(int); raters = list(df.columns)
naiveC = C.naive_lca_gibbs(Xc, rng, draws=5000, burn=1500)
reC = C.re_lca_gibbs_ms(Xc, rng, restarts=5)
tab = pd.DataFrame({"Se naive": naiveC["Se"].mean(0).round(2), "Se RE": reC["Se"].mean(0).round(2),
                    "Sp naive": naiveC["Sp"].mean(0).round(2), "Sp RE": reC["Sp"].mean(0).round(2),
                    "loading β": reC["beta"].mean(0).round(2)}, index=raters)
print(tab.to_string()); print(f"\nmean |loading| = {np.abs(reC['beta'].mean(0)).mean():.2f}  (large ⇒ strong conditional dependence among readers)")
fig, ax = plt.subplots(1, 2, figsize=(13.5, 4.3)); y = np.arange(len(raters))
ax[0].plot(naiveC["Se"].mean(0), y+0.12, "s", color=RED, ms=7, label="Se naive")
ax[0].plot(reC["Se"].mean(0), y+0.12, "o", color=BLUE, ms=7, label="Se random-effects")
ax[0].plot(naiveC["Sp"].mean(0), y-0.12, "s", color=ORANGE, ms=6, label="Sp naive")
ax[0].plot(reC["Sp"].mean(0), y-0.12, "o", color=GREEN, ms=6, label="Sp random-effects")
ax[0].set_yticks(y); ax[0].set_yticklabels(raters); ax[0].invert_yaxis(); ax[0].set_xlim(0.3,1)
ax[0].set_title("Naive vs dependence-corrected accuracies"); ax[0].set_xlabel("probability"); ax[0].legend(fontsize=7.5)
bm = reC["beta"].mean(0); blo,bhi = np.percentile(reC["beta"],[2.5,97.5],0)
ax[1].barh(y, bm, color=PURP); ax[1].errorbar(bm, y, xerr=[bm-blo,bhi-bm], fmt="none", color="k", capsize=3)
ax[1].set_yticks(y); ax[1].set_yticklabels(raters); ax[1].invert_yaxis()
ax[1].set_title("Reader loadings on the shared latent trait"); ax[1].set_xlabel("loading β")
plt.tight_layout(); plt.show()
print(f"The loadings are all clearly positive (mean β≈{np.abs(reC['beta'].mean(0)).mean():.2f}, none near zero) — the seven pathologists are conditionally dependent, sharing a latent")
print("'slide-difficulty' trait. That dependence is what surfaced as Project 2's third class. Correcting for it generally")
print("PULLS THE SENSITIVITIES DOWN (the naive model over-credited agreement to accuracy) and firms up the specificities —")
print("a more honest accuracy audit than the local-independence Hui–Walter model gives.")

   Se naive  Se RE  Sp naive  Sp RE  loading β
A      0.95   0.97      0.88   0.90       0.67
B      0.97   0.96      0.68   0.63       1.00
C      0.72   0.72      0.98   1.00       0.63
D      0.51   0.47      0.98   0.97       0.88
E      0.96   0.95      0.81   0.79       0.86
F      0.40   0.39      0.98   1.00       0.62
G      0.98   0.93      0.91   0.85       1.08

mean |loading| = 0.82  (large ⇒ strong conditional dependence among readers)
No description has been provided for this image 
The loadings are all clearly positive (mean β≈0.82, none near zero) — the seven pathologists are conditionally dependent, sharing a latent
'slide-difficulty' trait. That dependence is what surfaced as Project 2's third class. Correcting for it generally
PULLS THE SENSITIVITIES DOWN (the naive model over-credited agreement to accuracy) and firms up the specificities —
a more honest accuracy audit than the local-independence Hui–Walter model gives.

4. PyMC — sampling the random effects directly¶

Where the from-scratch sampler augments truncated-normal latents, PyMC samples the continuous random effects $b_i$ and the parameters directly with NUTS, marginalising the discrete class with logsumexp. Half-normal loadings fix the sign of the shared effect. It reaches the same conclusion — large loadings, dependence-corrected accuracies — on the carcinoma data.

In [5]:
import pymc as pm, pytensor.tensor as pt
Xf = Xc.astype(float); N, J = Xf.shape
with pm.Model() as m:
    pi = pm.Beta("pi", 1, 1)
    a = pm.Normal("a", 0, 3, shape=(2, J))
    beta = pm.HalfNormal("beta", 1.0, shape=J)               # >=0 anchors the effect's sign
    bz = pm.Normal("bz", 0, 1, shape=N)                      # subject random effects
    p = pm.math.invprobit(a[None,:,:] + beta[None,None,:]*bz[:,None,None])   # (N,2,J)
    ll = Xf[:,None,:]*pm.math.log(p) + (1-Xf[:,None,:])*pm.math.log(1-p)
    comp = pt.stack([pm.math.log(1-pi), pm.math.log(pi)])[None,:] + ll.sum(-1)
    pm.Potential("lik", pm.math.logsumexp(comp, axis=1).sum())
    s = pt.sqrt(1 + beta**2)
    pm.Deterministic("Se", pm.math.invprobit(a[1]/s)); pm.Deterministic("Sp", 1 - pm.math.invprobit(a[0]/s))
    idata = pm.sample(1000, tune=2000, chains=4, cores=1, target_accept=0.99, random_seed=0, progressbar=False)
# The model has no label anchor, so each chain may land in either labelling; a swapped
# chain maps (Se,Sp) -> (1-Sp,1-Se). An informative test needs Se+Sp > 1, which fixes
# the orientation. Align chains BEFORE pooling, or the average mixes Se with 1-Sp.
Se_c = idata.posterior["Se"].values; Sp_c = idata.posterior["Sp"].values   # (chain, draw, J)
_swap = (Se_c.mean(axis=(1,2)) + Sp_c.mean(axis=(1,2))) < 1.0
Se_a = np.where(_swap[:,None,None], 1-Sp_c, Se_c)
Sp_a = np.where(_swap[:,None,None], 1-Se_c, Sp_c)
print(f"label-aligned {int(_swap.sum())} of {len(_swap)} chains before pooling")
Se = Se_a.reshape(-1,J).mean(0); Sp = Sp_a.reshape(-1,J).mean(0)
betaP = idata.posterior["beta"].values.reshape(-1,J).mean(0)
cmp = pd.DataFrame({"Se PyMC": Se.round(2), "Se Gibbs": reC["Se"].mean(0).round(2),
                    "Sp PyMC": Sp.round(2), "Sp Gibbs": reC["Sp"].mean(0).round(2),
                    "β PyMC": betaP.round(2), "β Gibbs": reC["beta"].mean(0).round(2)}, index=raters)
print(cmp.to_string())
_ndiv = int(idata.sample_stats["diverging"].values.sum())
_dSe = float(np.abs(Se - reC["Se"].mean(0)).max())
_dB  = float(np.abs(betaP - reC["beta"].mean(0)).max())
print(f"\nBoth samplers reach the same qualitative conclusion — large positive loadings on every reader")
print(f"(PyMC mean |β|={np.abs(betaP).mean():.2f}, Gibbs {np.abs(reC['beta'].mean(0)).mean():.2f}) and the same direction of correction.")
print(f"They are NOT interchangeable in detail: sensitivities differ by up to {_dSe:.2f} and loadings by up to {_dB:.2f}")
print("across readers — the loadings are the weakly-identified part of a factor model.")
if _ndiv:
    print(f"\n{_ndiv} divergences: the factor geometry is hard for NUTS even at high target_accept. Treat the")
    print("from-scratch Albert–Chib Gibbs as the primary fit and this as a qualitative cross-check.")

g++ not available, if using conda: `conda install gxx`

Initializing NUTS using jitter+adapt_diag...

Sequential sampling (4 chains in 1 job)

NUTS: [pi, a, beta, bz]

Sampling 4 chains for 2_000 tune and 1_000 draw iterations (8_000 + 4_000 draws total) took 270 seconds.

There were 3560 divergences after tuning. Increase `target_accept` or reparameterize.

The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details

The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details

label-aligned 2 of 4 chains before pooling
   Se PyMC  Se Gibbs  Sp PyMC  Sp Gibbs  β PyMC  β Gibbs
A     0.95      0.97     0.86      0.90    0.94     0.67
B     0.97      0.96     0.66      0.63    1.26     1.00
C     0.73      0.72     1.00      1.00    0.64     0.63
D     0.50      0.47     0.98      0.97    1.19     0.88
E     0.98      0.95     0.80      0.79    0.95     0.86
F     0.40      0.39     0.99      1.00    1.08     0.62
G     0.98      0.93     0.90      0.85    1.00     1.08

Both samplers reach the same qualitative conclusion — large positive loadings on every reader
(PyMC mean |β|=1.01, Gibbs 0.82) and the same direction of correction.
They are NOT interchangeable in detail: sensitivities differ by up to 0.05 and loadings by up to 0.47
across readers — the loadings are the weakly-identified part of a factor model.

3560 divergences: the factor geometry is hard for NUTS even at high target_accept. Treat the
from-scratch Albert–Chib Gibbs as the primary fit and this as a qualitative cross-check.

5. Summary¶

Local independence is the assumption that makes latent class analysis work — and the one most worth checking. When items or tests are correlated given the class, the ordinary model is biased: here the naive Hui–Walter model over-stated the sensitivity and specificity of the correlated tests, because it mistook their shared correlation for individual accuracy. The random-effects (probit factor) model removes the bias by letting a continuous subject-level trait $b_i$ carry the within-class correlation; its loadings pinpoint exactly which tests are dependent. Recovery of the true accuracies is real but not automatic: the factor-mixture posterior is multimodal, and a run that settles in a partially-collapsed mode inflates every loading and deflates every sensitivity instead. The R notebook's fit recovers the truth closely; the Python run shown here lands in such a mode, which is precisely why the restart-and-select machinery exists — and why the loadings, not the point estimates, are the reliable diagnostic. On the carcinoma ratings the loadings were all clearly positive — the seven pathologists share a "slide-difficulty" trait — so the dependence-corrected sensitivities are meaningfully lower than the local-independence estimates, and Project 2's puzzling third class is explained as that very dependence.

The wider view. This model is a latent class and a latent trait at once — a factor mixture — and it sits at the crossroads of the portfolio: the Albert–Chib truncated-normal augmentation of the probit/Tobit models, the random effects of the hierarchical linear models, and the continuous-latent structure of confirmatory factor analysis. A leaner alternative for a few known correlated pairs is the fixed conditional-covariance model (Vacek 1985; Dendukuri & Joseph 2001), which adds an explicit covariance term per pair rather than a shared factor. Either way, the message is the same: test the local-independence assumption, and relax it when the data demand.

This completes a six-project tour of latent class analysis — foundations, choosing the number of classes, the nonparametric Dirichlet-process view, latent class regression (covariates predicting membership), diagnostic testing without a gold standard, and conditional dependence — all on the same conjugate data-augmentation core. The natural next step is its longitudinal cousin latent transition analysis (a hidden Markov model on the class, building on the Markov-switching work).