import numpy as np, pandas as pd, matplotlib.pyplot as plt
import lcasel as L, warnings; warnings.filterwarnings("ignore")
from scipy.special import logsumexp

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)

def loo_ic(X, g):
    # LOO information criterion (-2*elpd) by importance-sampling LOO on the Gibbs pointwise log-lik:
    # elpd_loo_i = -log mean_s exp(-ll[s,i]) ;  LOOIC = -2 sum_i elpd_loo_i.
    ll = L.pointwise_loglik(X, g["lam"], g["delta"])           # (S, N)
    S = ll.shape[0]
    elpd_i = np.log(S) - logsumexp(-ll, axis=0)
    return -2 * float(elpd_i.sum())

def scan_C(X, Cs, rng, n_starts=30, draws=1500, burn=500):
    rows = []
    for C in Cs:
        f = L.em_lca(X, C, rng, n_starts=n_starts); ic = L.info_criteria(f["loglik"], f["npar"], len(X))
        g = L.gibbs_lca(X, C, rng, draws=draws, burn=burn); w = L.waic(L.pointwise_loglik(X, g["lam"], g["delta"]))
        rows.append(dict(C=C, loglik=f["loglik"], npar=f["npar"], **ic, WAIC=w["WAIC"], LOOIC=loo_ic(X, g),
                         entropy=L.entropy_R2(f["resp"]) if C > 1 else np.nan))
    return pd.DataFrame(rows).set_index("C")
print("ready")

ready

lam_t = np.array([0.5, 0.3, 0.2])
delta_t = np.array([[0.9,0.9,0.8,0.8,0.2,0.2,0.1,0.1],
                    [0.1,0.2,0.1,0.2,0.9,0.8,0.9,0.8],
                    [0.8,0.2,0.8,0.2,0.8,0.2,0.8,0.2]])
Xsim, _ = L.simulate_lca(800, lam_t, delta_t, rng)
sim = scan_C(Xsim, range(1, 6), rng, n_starts=25)
print(sim[["loglik","AIC","BIC","SABIC","WAIC","LOOIC","entropy"]].round(1).to_string())
fig, ax = plt.subplots(1, 2, figsize=(13, 4.3))
for name,c in [("AIC",ORANGE),("BIC",BLUE),("SABIC",GREEN),("WAIC",RED),("LOOIC",PURP)]:
    y = sim[name] - sim[name].min()
    ax[0].plot(sim.index, y, "o-", color=c, lw=1.8, label=name)
ax[0].axvline(3, color="k", ls=":", lw=1.2); ax[0].text(3.05, ax[0].get_ylim()[1]*.8, "true C=3", fontsize=9)
ax[0].set_title("Criteria (offset to min) vs #classes — simulation"); ax[0].set_xlabel("number of classes C"); ax[0].set_ylabel("criterion − min"); ax[0].legend(fontsize=8)
mins = {n: int(sim[n].idxmin()) for n in ["AIC","BIC","SABIC","WAIC","LOOIC"]}
ax[1].bar(range(len(mins)), list(mins.values()), color=[ORANGE,BLUE,GREEN,RED,PURP])
ax[1].axhline(3, color="k", ls=":", lw=1.2); ax[1].set_xticks(range(len(mins))); ax[1].set_xticklabels(mins.keys())
ax[1].set_title("Selected C by each criterion (truth = 3)"); ax[1].set_ylabel("selected C"); ax[1].set_ylim(0,5.5)
plt.tight_layout(); plt.show()
print("BIC, SABIC, WAIC and LOO recover the true 3 classes; AIC's lighter penalty makes it prone to over-extraction.")

   loglik     AIC     BIC   SABIC    WAIC   LOOIC  entropy
C                                                         
1 -4324.7  8665.4  8702.9  8677.5  8665.6  8665.6      NaN
2 -3555.4  7144.7  7224.3  7170.4  7146.2  7146.2      0.9
3 -3409.0  6870.0  6991.8  6909.3  6870.5  6870.5      0.8
4 -3399.8  6869.6  7033.5  6922.4  6872.5  6872.6      0.8
5 -3392.5  6872.9  7079.0  6939.3  6872.7  6872.9      0.7

BIC, SABIC, WAIC and LOO recover the true 3 classes; AIC's lighter penalty makes it prone to over-extraction.

# the bootstrap LRT on the simulated data: C=2 vs 3 should reject; C=3 vs 4 should not
b23 = L.blrt(Xsim, 2, rng, B=100); b34 = L.blrt(Xsim, 3, rng, B=100)
fig, ax = plt.subplots(1, 2, figsize=(13, 3.8))
for k,(b,c) in enumerate([(b23,"2 vs 3"), (b34,"3 vs 4")]):
    ax[k].hist(b[0]["LR_null"] if isinstance(b,tuple) else b["LR_null"], bins=25, color=GREY, edgecolor="white", lw=.3, label="null LR (bootstrap)")
    obs = b["LR_obs"]; ax[k].axvline(obs, color=RED, lw=2, label=f"observed LR = {obs:.1f}")
    ax[k].set_title(f"BLRT: {c} classes   (p = {b['pval']:.3f})"); ax[k].set_xlabel("likelihood-ratio statistic"); ax[k].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"C=2 vs 3: p={b23['pval']:.3f} (reject — a third class is needed).  C=3 vs 4: p={b34['pval']:.3f} (do not reject — stop at 3).")
print("The bootstrap LRT pinpoints the true number of classes; note its observed LR falls far in the tail only for the real extra class.")

C=2 vs 3: p=0.010 (reject — a third class is needed).  C=3 vs 4: p=0.366 (do not reject — stop at 3).
The bootstrap LRT pinpoints the true number of classes; note its observed LR falls far in the tail only for the real extra class.

df = pd.read_csv("carcinoma.csv")               # raters A..G coded 1/2
X = (df.values - 1).astype(int)                 # 0/1 ; 1 = carcinoma call
raters = list(df.columns); N, J = X.shape
print(f"{N} slides rated by {J} pathologists {raters}")
print("carcinoma-call rate per pathologist:", dict(zip(raters, X.mean(0).round(2))))
print("pairwise agreement (fraction of slides two raters agree on):")
agr = np.array([[np.mean(X[:,i]==X[:,j]) for j in range(J)] for i in range(J)])
print(pd.DataFrame(agr.round(2), index=raters, columns=raters).to_string())

118 slides rated by 7 pathologists ['A', 'B', 'C', 'D', 'E', 'F', 'G']
carcinoma-call rate per pathologist: {'A': np.float64(0.56), 'B': np.float64(0.67), 'C': np.float64(0.38), 'D': np.float64(0.27), 'E': np.float64(0.6), 'F': np.float64(0.21), 'G': np.float64(0.56)}
pairwise agreement (fraction of slides two raters agree on):
      A     B     C     D     E     F     G
A  1.00  0.84  0.82  0.71  0.86  0.65  0.90
B  0.84  1.00  0.69  0.60  0.88  0.54  0.87
C  0.82  0.69  1.00  0.79  0.78  0.76  0.82
D  0.71  0.60  0.79  1.00  0.65  0.84  0.71
E  0.86  0.88  0.78  0.65  1.00  0.61  0.91
F  0.65  0.54  0.76  0.84  0.61  1.00  0.65
G  0.90  0.87  0.82  0.71  0.91  0.65  1.00

sel = scan_C(X, range(1, 6), rng, n_starts=40, draws=2500, burn=700)
show = sel[["loglik","AIC","BIC","SABIC","CAIC","WAIC","LOOIC","entropy"]].round(1)
print(show.to_string())
gof3 = L.gof(X, *[L.em_lca(X,3,rng,40)[k] for k in ("lam","delta")])
print(f"\n(full-table fit at C=3: G²={gof3['G2']:.1f}, df={gof3['df']}, but {gof3['sparsity']*100:.0f}% of cells have expected<1 — G² unreliable.)")
# two df conventions are in circulation; state ours so the R notebook does not look contradictory
_npar3 = (3 - 1) + 3 * X.shape[1]
print(f" df counts the 2^J−1 = {2**X.shape[1]-1} free cells minus {_npar3} parameters. poLCA instead reports")
print(f" residual df as N−npar = {X.shape[0]}−{_npar3} = {X.shape[0]-_npar3}, so the same G² is quoted against a")
print(" different df in the R notebook — a difference of convention, not of result.")
fig, ax = plt.subplots(1, 2, figsize=(13, 4.3))
for name,c in [("AIC",ORANGE),("BIC",BLUE),("SABIC",GREEN),("CAIC",GREY),("WAIC",RED),("LOOIC",PURP)]:
    ax[0].plot(sel.index, sel[name]-sel[name].min(), "o-", color=c, lw=1.7, label=f"{name} (min at C={int(sel[name].idxmin())})")
ax[0].set_title("Information criteria vs #classes — carcinoma"); ax[0].set_xlabel("number of classes C"); ax[0].set_ylabel("criterion − min"); ax[0].legend(fontsize=8)
ax[1].plot(sel.index[1:], sel["entropy"][1:], "s-", color=PURP, lw=2)
ax[1].set_title("Entropy R² (class separation)"); ax[1].set_xlabel("number of classes C"); ax[1].set_ylabel("entropy R²"); ax[1].set_ylim(0.85,1.0)
plt.tight_layout(); plt.show()
print("Every penalised criterion (BIC, SABIC, CAIC, WAIC, LOO) bottoms out at C=3; AIC agrees here. The 3-class model is favoured.")

   loglik     AIC     BIC   SABIC    CAIC    WAIC   LOOIC  entropy
C                                                                 
1  -524.5  1062.9  1082.3  1060.2  1089.3  1063.2  1063.2      NaN
2  -317.3   664.5   706.1   658.7   721.1   665.5   665.8      1.0
3  -293.7   633.4   697.1   624.4   720.1   633.3   633.5      0.9
4  -289.3   640.6   726.5   628.5   757.5   635.1   635.4      0.9
5  -286.9   651.8   759.8   636.5   798.8   636.8   637.0      0.9

(full-table fit at C=3: G²=15.3, df=104, but 86% of cells have expected<1 — G² unreliable.)
 df counts the 2^J−1 = 127 free cells minus 23 parameters. poLCA instead reports
 residual df as N−npar = 118−23 = 95, so the same G² is quoted against a
 different df in the R notebook — a difference of convention, not of result.

Every penalised criterion (BIC, SABIC, CAIC, WAIC, LOO) bottoms out at C=3; AIC agrees here. The 3-class model is favoured.

# BLRT confirms C=3: reject 2, do not reject 3->4
c23 = L.blrt(X, 2, rng, B=150); c34 = L.blrt(X, 3, rng, B=150)
fig, ax = plt.subplots(1, 2, figsize=(13, 3.8))
for k,(b,lab) in enumerate([(c23,"2 vs 3"),(c34,"3 vs 4")]):
    ax[k].hist(b["LR_null"], bins=28, color=GREY, edgecolor="white", lw=.3)
    ax[k].axvline(b["LR_obs"], color=RED, lw=2, label=f"observed = {b['LR_obs']:.1f}")
    ax[k].set_title(f"BLRT {lab} classes (p={b['pval']:.3f})"); ax[k].set_xlabel("LR statistic"); ax[k].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"C=2 vs 3: p={c23['pval']:.3f} — a third class is needed.  C=3 vs 4: p={c34['pval']:.3f} — a fourth is not. The BLRT selects C=3.")

C=2 vs 3: p=0.007 — a third class is needed.  C=3 vs 4: p=0.159 — a fourth is not. The BLRT selects C=3.

# validate Chib against the exact 1-class marginal likelihood, then scan C
exact1 = L.exact_logml_c1(X); chib1 = L.chib_logml(X, 1, rng, draws=4000, burn=1000)["logml"]
print(f"C=1 marginal likelihood — exact {exact1:.3f}  vs  Chib {chib1:.3f}  (difference {abs(exact1-chib1):.4f}): the estimator is correct.\n")
ml = [];
for C_ in range(1, 6):
    m = L.chib_logml(X, C_, rng, draws=6000, burn=1500)["logml"]; ml.append(m)
ml = np.array(ml); Cs = np.arange(1, 6)
mBIC = -0.5 * sel["BIC"].values                                   # BIC approximates -2 log ML
tab = pd.DataFrame({"log p(x|C)": ml.round(2), "−BIC/2": mBIC.round(2),
                    "2·log BF vs C−1": np.concatenate([[np.nan], 2*np.diff(ml)]).round(1)}, index=Cs)
tab.index.name = "C"; print(tab.to_string())
fig, ax = plt.subplots(1, 2, figsize=(13, 4.2))
ax[0].plot(Cs, ml, "o-", color=BLUE, lw=2, ms=8, label="Chib log marginal likelihood")
ax[0].plot(Cs, mBIC, "s--", color=ORANGE, lw=1.6, label="−BIC/2 (the approximation)")
ax[0].axvline(Cs[ml.argmax()], color="k", ls=":", lw=1.2)
ax[0].set_title("Marginal likelihood peaks at the chosen C"); ax[0].set_xlabel("number of classes C"); ax[0].set_ylabel("log p(x | C)"); ax[0].legend(fontsize=8)
bf = 2*np.diff(ml)
ax[1].bar(Cs[1:], bf, color=np.where(bf>0, GREEN, RED))
ax[1].axhline(10, color="k", ls=":", lw=1); ax[1].axhline(-10, color="k", ls=":", lw=1)
ax[1].set_title("2·log Bayes factor: C vs C−1  (|·|>10 = decisive)"); ax[1].set_xlabel("C"); ax[1].set_ylabel("2·log BF")
plt.tight_layout(); plt.show()
# label the Bayes factors on the Kass & Raftery (1995) scale rather than by hand,
# so the wording can never drift from the number (and matches the |.|>10 line above)
def _kr(v):
    a = abs(v)
    return "barely worth a mention" if a < 2 else "positive" if a < 6 else "strong" if a < 10 else "very strong"
b32, b43 = 2*(ml[2]-ml[1]), 2*(ml[3]-ml[2])
print(f"The evidence is maximised at C={Cs[ml.argmax()]}: 2·log BF(3:2)={b32:+.1f} — {_kr(b32)} evidence FOR a 3rd class;")
print(f"2·log BF(4:3)={b43:+.1f} — {_kr(b43)} against a 4th. −BIC/2 tracks the exact log-evidence — BIC is its large-N approximation.")

C=1 marginal likelihood — exact -540.068  vs  Chib -540.068  (difference 0.0000): the estimator is correct.

   log p(x|C)  −BIC/2  2·log BF vs C−1
C                                     
1     -540.07 -541.16              NaN
2     -353.65 -353.04            372.8
3     -349.77 -348.57              7.8
4     -350.74 -363.23             -1.9
5     -350.97 -379.91             -0.5

The evidence is maximised at C=3: 2·log BF(3:2)=+7.8 — strong evidence FOR a 3rd class;
2·log BF(4:3)=-1.9 — barely worth a mention against a 4th. −BIC/2 tracks the exact log-evidence — BIC is its large-N approximation.

f3 = L.em_lca(X, 3, rng, n_starts=60)
prof = pd.DataFrame(f3["delta"].round(2), columns=raters,
                    index=[f"class {c+1} (λ={f3['lam'][c]:.2f})" for c in range(3)])
print("P(carcinoma call | class), by pathologist:\n", prof.to_string())
names = ["clear carcinoma", "clear benign", "ambiguous / disputed"]
fig, ax = plt.subplots(1, 2, figsize=(13.5, 4.4))
xp = np.arange(J)
for c,(col,nm) in enumerate(zip([RED,GREEN,ORANGE], names)):
    ax[0].plot(xp, f3["delta"][c], "o-", color=col, lw=2, ms=7, label=f"{nm} (λ={f3['lam'][c]:.2f})")
ax[0].set_xticks(xp); ax[0].set_xticklabels(raters); ax[0].set_ylim(0,1); ax[0].set_ylabel("P(carcinoma call)")
ax[0].set_title("Three latent slide types"); ax[0].legend(fontsize=8.5)
# bivariate-residual fit: are pairwise agreements reproduced?
bvr = L.bivariate_residuals(X, f3["lam"], f3["delta"])
im = ax[1].imshow(bvr, cmap="OrRd", vmin=0, vmax=6)
ax[1].set_xticks(range(J)); ax[1].set_xticklabels(raters); ax[1].set_yticks(range(J)); ax[1].set_yticklabels(raters)
for i in range(J):
    for j in range(J):
        if i!=j: ax[1].text(j,i,f"{bvr[i,j]:.1f}",ha="center",va="center",fontsize=7)
ax[1].set_title(f"Bivariate residuals (χ²₁; max {bvr.max():.1f})"); plt.colorbar(im, ax=ax[1], shrink=.8)
plt.tight_layout(); plt.show()
print("Class interpretation: ~44% clear carcinoma (all raters call it), ~37% clear benign (none do), ~18% ambiguous")
print("(raters split — high on B,E,G, low on C,D,F). Bivariate residuals are near the χ²₁ scale: the 3-class model")
print("reproduces the pairwise rater agreements well (only a pair or two mildly exceed 3.8).")

P(carcinoma call | class), by pathologist:
                      A     B     C     D     E     F     G
class 1 (λ=0.44)  1.00  0.98  0.86  0.59  1.00  0.48  1.00
class 2 (λ=0.37)  0.06  0.14  0.00  0.00  0.06  0.00  0.00
class 3 (λ=0.18)  0.51  1.00  0.00  0.06  0.75  0.00  0.63

Class interpretation: ~44% clear carcinoma (all raters call it), ~37% clear benign (none do), ~18% ambiguous
(raters split — high on B,E,G, low on C,D,F). Bivariate residuals are near the χ²₁ scale: the 3-class model
reproduces the pairwise rater agreements well (only a pair or two mildly exceed 3.8).

Latent Class Analysis — Choosing the Number of Classes¶

Model selection for LCA, validated and applied¶

The toolkit¶

Why validation first¶

1. Validation — can the criteria recover a known number of classes?¶

2. The carcinoma data — how many kinds of slide?¶

3. The full toolkit on carcinoma¶

4. The marginal likelihood & Bayes factors — Chib's method¶

5. The chosen model — three kinds of slide¶

6. Summary¶