IRT Model Comparison — 1PL vs 2PL vs 3PL¶
Do the extra item parameters earn their keep?¶
The three workhorse dichotomous IRT models are nested: $$\text{1PL: }\Phi(a\theta_i+d_j)\ \subset\ \text{2PL: }\Phi(a_j\theta_i+d_j)\ \subset\ \text{3PL: }c_j+(1-c_j)\Phi(a_j\theta_i+d_j).$$ The 2PL lets each item have its own discrimination $a_j$; the 3PL adds a lower asymptote $c_j$ for guessing. Because each model contains the last, the in-sample fit can only improve as parameters are added — a richer model never fits the training data worse. So a raw likelihood cannot choose between them. WAIC and PSIS-LOO can: they estimate the expected log pointwise predictive density (elpd) on new data, subtracting a complexity penalty, so a bigger model is preferred only when it genuinely predicts better.
One subtlety decides whether the comparison is trustworthy. Scoring the models by the conditional likelihood — treating each of the $N\times J$ responses as the unit and carrying the $N$ ability parameters — is badly behaved: the penalty is swamped by the abilities. The fix is to integrate $\theta$ out so the person is the exchangeable unit, by Gauss–Hermite quadrature — the same marginalisation used to score the latent-class models in Latent Class Analysis — Choosing the Number of Classes. We fit all three models from scratch (Albert–Chib augmentation), compare them on two datasets that pull in opposite directions — the near-Rasch LSAT where parsimony should win and the multiple-choice SAT12 where guessing should pay — check item fit with posterior predictive ICCs, and confirm the ranking in PyMC.
import numpy as np, pandas as pd, matplotlib.pyplot as plt, arviz as az
import irtcompare as M
rng = np.random.default_rng(7)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; PURP="#6b46c1"; GREY="#718096"
def fit_all(X, dr, bn):
return {"1PL":(M.irt1pl_gibbs(X,rng,draws=dr,burn=bn),"1pl"),
"2PL":(M.irt2pl_gibbs(X,rng,draws=dr,burn=bn),"2pl"),
"3PL":(M.irt3pl_gibbs(X,rng,draws=dr,burn=bn),"3pl")}
def marg_waic(X, fits):
return {k:(M.waic(M.marginal_loglik(X,r,m)), r, m) for k,(r,m) in fits.items()}
def loo_of(LLm): # arviz PSIS-LOO from a (draws,N) marginal log-lik
idata = az.from_dict({"log_likelihood":{"obs":LLm[None]}}) # (chain=1, draw, N)
return az.loo(idata, reff=1.0)
print("1PL < 2PL < 3PL are nested; more parameters always fit training better, so WAIC/LOO -- which estimate")
print("out-of-sample prediction with a complexity penalty on the theta-integrated (marginal) likelihood -- decide.")
1PL < 2PL < 3PL are nested; more parameters always fit training better, so WAIC/LOO -- which estimate out-of-sample prediction with a complexity penalty on the theta-integrated (marginal) likelihood -- decide.
1. Why the likelihood must be marginal¶
Take the LSAT data and score the 1PL by WAIC two ways. The conditional version counts each response as an observation and carries the 1000 ability parameters; the marginal version integrates $\theta$ out and counts each person once. The effective number of parameters $p_{\text{WAIC}}$ tells the story: conditional WAIC reports hundreds (it is penalising the abilities), while marginal WAIC reports a handful — the actual item-parameter count. Only the marginal version is a fair basis for comparison.
LSAT = pd.read_csv("lsat.csv").to_numpy().astype(float)
r1 = M.irt1pl_gibbs(LSAT, rng, draws=1500, burn=800)
wc = M.waic(r1["loglik"]) # conditional: unit = response, carries theta
wm = M.waic(M.marginal_loglik(LSAT, r1, "1pl")) # marginal: unit = person, theta integrated out
print(f"1PL on LSAT ({LSAT.shape[0]} persons x {LSAT.shape[1]} items):")
print(f" CONDITIONAL WAIC: p_eff = {wc['p_eff']:6.1f} (swamped by the {LSAT.shape[0]} ability parameters)")
print(f" MARGINAL WAIC: p_eff = {wm['p_eff']:6.1f} (~ the {LSAT.shape[1]} easinesses + 1 slope = the real item parameters)")
print("The conditional penalty is meaningless for comparing item structure; the marginal one recovers the true")
print("model complexity. Every comparison below uses the marginal (theta-integrated) likelihood by Gauss-Hermite.")
1PL on LSAT (1000 persons x 5 items): CONDITIONAL WAIC: p_eff = 278.3 (swamped by the 1000 ability parameters) MARGINAL WAIC: p_eff = 6.5 (~ the 5 easinesses + 1 slope = the real item parameters) The conditional penalty is meaningless for comparing item structure; the marginal one recovers the true model complexity. Every comparison below uses the marginal (theta-integrated) likelihood by Gauss-Hermite.
Why the likelihood has to be integrated, and how¶
The choice being made here is about the unit of prediction. The conditional likelihood asks how well the model predicts one more response from an examinee whose ability has already been estimated — but that examinee's ability was estimated from those very responses, so the question flatters the model and the penalty ends up counting one parameter per person. The marginal likelihood asks the question a test developer actually cares about: how well does this item set predict the responses of a new person, whose ability is unknown and must be treated as drawn from the population?
Making that switch means removing $\theta_i$ by integration:
$$p(\mathbf{x}_i \mid \text{items}) = \int \prod_j p(x_{ij} \mid \theta,\ \text{item } j)\ \phi(\theta)\, d\theta$$
The integral has no closed form, but it is one-dimensional and its weight function is exactly the standard normal, which is the textbook case for Gauss–Hermite quadrature: approximate the integral by evaluating the integrand at a fixed set of nodes and taking a weighted sum, $\int f(\theta)\phi(\theta)\,d\theta \approx \sum_k w_k f(x_k)$. The nodes and weights are chosen so the approximation is exact for polynomials up to a given degree, so a few dozen nodes give many digits of accuracy here — far cheaper and far more stable than sampling the integral would be.
The result is a likelihood in which the only free parameters are the item parameters, which is what makes the penalty count them and the comparison mean what it says.
2. LSAT — where parsimony wins¶
The five LSAT items are famously near-Rasch: the discriminations barely differ, so the 2PL's per-item slopes should be wasted, and the items are easy multiple-guess-proof logic items, so guessing should not help either. Fit all three, rank by marginal WAIC and PSIS-LOO, and read the pairwise elpd differences.
fL = fit_all(LSAT, 1500, 800); WL = marg_waic(LSAT, fL)
rows=[]
for k,(w,r,m) in WL.items():
lo = loo_of(M.marginal_loglik(LSAT, r, m))
rows.append((k, w["elpd"], w["p_eff"], w["waic"], float(lo.elpd), float(lo.p), float(lo.pareto_k.max())))
tab=pd.DataFrame(rows, columns=["model","elpd_waic","p_waic","WAIC","elpd_loo","p_loo","max_k"]).set_index("model")
print(tab.round(1).to_string())
print("(max_k = worst Pareto-k; all < 0.7 means PSIS-LOO is reliable and agrees with WAIC.)")
d21=M.compare_elpd(WL["2PL"][0],WL["1PL"][0]); d32=M.compare_elpd(WL["3PL"][0],WL["2PL"][0])
print(f"\n2PL - 1PL elpd: {d21['delta']:+.1f} +/- {d21['se']:.1f} (adding per-item slopes: {'no gain' if d21['delta']<2*d21['se'] else 'gain'})")
print(f"3PL - 2PL elpd: {d32['delta']:+.1f} +/- {d32['se']:.1f} (adding guessing: {'no gain' if abs(d32['delta'])<2*d32['se'] else 'gain'})")
elpd_waic p_waic WAIC elpd_loo p_loo max_k model 1PL -2473.0 5.9 4946.1 -2473.1 5.9 0.2 2PL -2476.8 10.0 4953.6 -2476.8 10.1 0.2 3PL -2477.6 10.6 4955.2 -2477.6 10.6 0.1 (max_k = worst Pareto-k; all < 0.7 means PSIS-LOO is reliable and agrees with WAIC.) 2PL - 1PL elpd: -3.8 +/- 1.3 (adding per-item slopes: no gain) 3PL - 2PL elpd: -0.8 +/- 1.3 (adding guessing: no gain)
Reading elpd, the penalty, and the $\pm$¶
Every comparison in this notebook is reported as elpd — expected log pointwise predictive density. It is the average log-probability the model would assign to data it has not seen, summed over units, so it is on a log-probability scale: higher is better, it has no absolute meaning, and only differences between models on the same data are interpretable. WAIC is the same quantity on the deviance scale, $\text{WAIC} = -2\,\text{elpd}$, which is why the WAIC column runs the opposite way.
The penalty $p_{\text{eff}}$ is an effective number of parameters, estimated from how much the log-likelihood varies across the posterior rather than counted from the model definition. That is what makes it meaningful for hierarchical models, where parameters are partially pooled and the count of free parameters overstates the flexibility actually used — and it is why the conditional version above returns 278.3 for a model with six item parameters, having counted 1,000 abilities that the marginal likelihood integrates away.
The $\pm$ is a standard error of the difference, computed from the pointwise variation across units, not a posterior interval. It is the right yardstick for these comparisons: a difference of $-3.8 \pm 1.3$ is under three standard errors and points the wrong way for complexity, while $+42.0 \pm 11.5$ is more than three and points toward it. Differences smaller than one or two standard errors are not evidence for either model.
It also helps to convert to a per-unit scale, because elpd totals look larger than they are. The SAT12 gain of 42 elpd is spread over 600 examinees — about 0.07 in log-probability per person, a real but small improvement in predicting one person's whole response pattern. Model comparison is telling you the 3PL predicts better, not that it predicts well; whether it fits at all is the separate posterior-predictive question the next section asks.
Pareto-$k$ is the diagnostic for the second criterion. PSIS-LOO estimates leave-one-out prediction by importance-reweighting the fit rather than refitting 600 times, and $k$ measures how heavy-tailed those weights are for each unit. Below 0.7 the approximation is trustworthy; above it, that unit's contribution is unreliable and should be refitted exactly. All the values here are below 0.5, which is why LOO and WAIC agree throughout and neither needs a caveat.
best_L = tab["elpd_waic"].idxmax()
fig,ax=plt.subplots(1,2,figsize=(12,4.2))
ks=list(WL); el=[WL[k][0]["elpd"] for k in ks]; se=[WL[k][0]["se"] for k in ks]
ax[0].errorbar(el, range(len(ks)), xerr=se, fmt="o", color=BLUE, capsize=4)
ax[0].scatter([WL[best_L][0]["elpd"]],[ks.index(best_L)], color=GREEN, zorder=3, s=90, label="best")
ax[0].set_yticks(range(len(ks))); ax[0].set_yticklabels(ks); ax[0].set_xlabel("elpd (higher = better prediction)")
ax[0].set_title(f"LSAT marginal WAIC — winner: {best_L}"); ax[0].legend(frameon=False,fontsize=8)
a2 = fL["2PL"][0]["a"].mean(0)
ax[1].bar(range(len(a2)), a2, color=GREY); ax[1].axhline(a2.mean(),color=RED,ls="--",label=f"mean {a2.mean():.2f}")
ax[1].set_xticks(range(len(a2))); ax[1].set_xticklabels([f"i{j+1}" for j in range(len(a2))])
ax[1].set_ylabel("2PL discrimination $a_j$"); ax[1].set_title("LSAT discriminations barely differ → near-Rasch"); ax[1].legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print(f"WAIC/LOO prefer the {best_L}: the discriminations are nearly equal (near-Rasch), so the 2PL's extra slopes do")
print("not improve prediction, and the 3PL's guessing parameter is idle. On simple data, parsimony wins -- exactly")
print("what the complexity penalty is for.")
WAIC/LOO prefer the 1PL: the discriminations are nearly equal (near-Rasch), so the 2PL's extra slopes do not improve prediction, and the 3PL's guessing parameter is idle. On simple data, parsimony wins -- exactly what the complexity penalty is for.
3. Does the chosen model actually fit? — posterior-predictive ICCs¶
A model can win a comparison and still fit poorly in absolute terms. The item-fit check bins people by ability and overlays the observed proportion-correct on the model's expected curve. If the winning LSAT model tracks the empirical points, its fit is adequate — not merely less bad than the alternatives.
rB, mB = fL[best_L][0], fL[best_L][1]
if mB=="1pl": Pm = M.Phi(rB["a"].mean()*rB["theta"].mean(0)[:,None] + rB["d"].mean(0)[None,:])
elif mB=="2pl": Pm = M.prob2(rB["a"].mean(0), rB["d"].mean(0), rB["theta"].mean(0))
else: Pm = M.prob3(rB["a"].mean(0), rB["d"].mean(0), rB["c"].mean(0), rB["theta"].mean(0))
th = rB["theta"].mean(0)
fig,ax=plt.subplots(1,5,figsize=(14,3),sharey=True)
for j in range(5):
xb,ob,mb = M.item_fit_curve(th, Pm, LSAT, j, nbin=6)
ax[j].plot(xb, mb, color=BLUE, lw=2, label="model")
ax[j].scatter(xb, ob, color=RED, zorder=3, label="observed")
ax[j].set_title(f"item {j+1}", fontsize=9); ax[j].set_xlabel(r"$\theta$")
if j==0: ax[j].set_ylabel("P(correct)"); ax[j].legend(frameon=False,fontsize=7)
plt.suptitle(f"LSAT item fit — observed vs {best_L} expected proportion-correct by ability group", y=1.04)
plt.tight_layout(); plt.show()
print(f"The observed points sit on the {best_L} curves across the ability range: the parsimonious model is not just")
print("the comparative winner, it fits the LSAT items well in absolute terms.")
The observed points sit on the 1PL curves across the ability range: the parsimonious model is not just the comparative winner, it fits the LSAT items well in absolute terms.
4. SAT12 — where guessing pays off¶
The SAT12 items are multiple-choice, so a low-ability examinee can get an item right by guessing — precisely the behaviour the 3PL's lower asymptote $c_j$ captures. Here the ordering should reverse: the extra parameters should improve prediction. Fit all three and read the marginal WAIC.
SAT = pd.read_csv("sat12.csv").to_numpy().astype(float)
fS = fit_all(SAT, 1200, 600); WS = marg_waic(SAT, fS)
rows=[]
for k,(w,r,m) in WS.items():
lo=loo_of(M.marginal_loglik(SAT,r,m))
rows.append((k,w["elpd"],w["p_eff"],w["waic"],float(lo.elpd),float(lo.p),float(lo.pareto_k.max())))
tabS=pd.DataFrame(rows,columns=["model","elpd_waic","p_waic","WAIC","elpd_loo","p_loo","max_k"]).set_index("model")
print(f"SAT12 ({SAT.shape[0]} persons x {SAT.shape[1]} multiple-choice items):"); print(tabS.round(1).to_string())
best_S=tabS["elpd_waic"].idxmax()
d21=M.compare_elpd(WS["2PL"][0],WS["1PL"][0]); d32=M.compare_elpd(WS["3PL"][0],WS["2PL"][0])
print(f"\n2PL - 1PL elpd: {d21['delta']:+.1f} +/- {d21['se']:.1f} (per-item slopes matter here)")
print(f"3PL - 2PL elpd: {d32['delta']:+.1f} +/- {d32['se']:.1f} ({'>2 SE: guessing earns its keep' if d32['delta']>2*d32['se'] else 'marginal'})")
SAT12 (600 persons x 32 multiple-choice items):
elpd_waic p_waic WAIC elpd_loo p_loo max_k
model
1PL -9681.8 33.7 19363.6 -9681.9 33.8 0.3
2PL -9561.1 66.0 19122.2 -9561.3 66.2 0.5
3PL -9519.1 77.2 19038.1 -9519.3 77.4 0.4
2PL - 1PL elpd: +120.7 +/- 19.0 (per-item slopes matter here)
3PL - 2PL elpd: +42.0 +/- 11.5 (>2 SE: guessing earns its keep)
cj = fS["3PL"][0]["c"].mean(0)
fig,ax=plt.subplots(1,2,figsize=(12,4.2))
ks=list(WS); el=[WS[k][0]["elpd"] for k in ks]; se=[WS[k][0]["se"] for k in ks]
ax[0].errorbar(el, range(len(ks)), xerr=se, fmt="o", color=BLUE, capsize=4)
ax[0].scatter([WS[best_S][0]["elpd"]],[ks.index(best_S)], color=GREEN, zorder=3, s=90, label="best")
ax[0].set_yticks(range(len(ks))); ax[0].set_yticklabels(ks); ax[0].set_xlabel("elpd (higher = better)")
ax[0].set_title(f"SAT12 marginal WAIC — winner: {best_S}"); ax[0].legend(frameon=False,fontsize=8)
ax[1].hist(cj, bins=12, color=PURP, alpha=.8); ax[1].axvline(cj.mean(),color=RED,ls="--",label=f"mean {cj.mean():.2f}")
ax[1].axvline(0.2,color=GREY,ls=":",label="1/5 chance")
ax[1].set_xlabel("estimated guessing $c_j$"); ax[1].set_ylabel("items"); ax[1].set_title("SAT12 guessing parameters cluster near chance"); ax[1].legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print(f"Now the {best_S} wins: the multiple-choice items have real guessing (c_j clustered near the ~1/5 chance rate),")
print("so the lower asymptote genuinely improves prediction and the extra parameters are justified. Same criterion,")
print("opposite verdict from LSAT — model comparison adapts the complexity to the data, not to a fixed preference.")
Now the 3PL wins: the multiple-choice items have real guessing (c_j clustered near the ~1/5 chance rate), so the lower asymptote genuinely improves prediction and the extra parameters are justified. Same criterion, opposite verdict from LSAT — model comparison adapts the complexity to the data, not to a fixed preference.
5. Cross-check in PyMC¶
We refit the 2PL and 3PL in PyMC, then score them with the same marginal-WAIC machinery applied to PyMC's item-parameter posterior — an apples-to-apples check that the ranking does not depend on the sampler. We use SAT12, where the models genuinely differ.
import pymc as pm
def pymc_fit(X, three):
N,J=X.shape
with pm.Model():
th=pm.Normal("theta",0,1,shape=N); a=pm.HalfNormal("a",2,shape=J); d=pm.Normal("d",0,3,shape=J)
p=pm.math.invprobit(th[:,None]*a[None,:]+d[None,:])
if three:
c=pm.Beta("c",1,4,shape=J); p=c[None,:]+(1-c[None,:])*p
pm.Bernoulli("x", p=p, observed=X)
idata=pm.sample(500, tune=800, chains=4, target_accept=0.9, random_seed=1, progressbar=False)
po=idata.posterior; S=po.dims["chain"]*po.dims["draw"]
res={"a":po["a"].stack(s=("chain","draw")).values.T, "d":po["d"].stack(s=("chain","draw")).values.T,
"theta":po["theta"].stack(s=("chain","draw")).values.T}
if three: res["c"]=po["c"].stack(s=("chain","draw")).values.T
# keep the sampler diagnostics so they can be REPORTED rather than left in the warning stream
res["_ndiv"] = int(idata.sample_stats["diverging"].values.sum())
res["_ndraw"] = int(idata.sample_stats["diverging"].size)
return res
res2=pymc_fit(SAT, False); res3=pymc_fit(SAT, True)
w2=M.waic(M.marginal_loglik(SAT,res2,"2pl")); w3=M.waic(M.marginal_loglik(SAT,res3,"3pl"))
print("PyMC on SAT12 (marginal WAIC):")
print(f" 2PL elpd {w2['elpd']:8.1f} 3PL elpd {w3['elpd']:8.1f}")
print(f" PyMC prefers: {'3PL' if w3['elpd']>w2['elpd'] else '2PL'} | from-scratch preferred: {best_S}")
print(f" 3PL - 2PL elpd (PyMC) {w3['elpd']-w2['elpd']:+.1f} (from-scratch {WS['3PL'][0]['elpd']-WS['2PL'][0]['elpd']:+.1f})")
print("Both samplers reach the same verdict on SAT12: the guessing parameter improves out-of-sample prediction.")
print("\nThe agreement is on the RANKING, and that is the part worth trusting, because NUTS had a hard")
print("time here: %d of %d draws diverged on the 3PL against %d on the 2PL." % (res3["_ndiv"], res3["_ndraw"], res2["_ndiv"]))
print("The 3PL geometry is the reason -- guessing is weakly identified, so c trades off against d and the")
print("sampler is dragged through a funnel that the conjugate Gibbs sampler never enters, since it draws")
print("each block from its exact full conditional. A gap of %.0f elpd survives that comfortably; a" % abs(w3["elpd"]-w2["elpd"]))
print("comparison decided by a handful of elpd points would not, and should not be attempted from this fit.")
g++ not available, if using conda: `conda install gxx`
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [theta, a, d]
Sampling 4 chains for 800 tune and 500 draw iterations (3_200 + 2_000 draws total) took 23 seconds.
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
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [theta, a, d, c]
Sampling 4 chains for 800 tune and 500 draw iterations (3_200 + 2_000 draws total) took 55 seconds.
There were 499 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
PyMC on SAT12 (marginal WAIC): 2PL elpd -9561.1 3PL elpd -9522.5 PyMC prefers: 3PL | from-scratch preferred: 3PL 3PL - 2PL elpd (PyMC) +38.6 (from-scratch +42.0) Both samplers reach the same verdict on SAT12: the guessing parameter improves out-of-sample prediction. The agreement is on the RANKING, and that is the part worth trusting, because NUTS had a hard time here: 499 of 2000 draws diverged on the 3PL against 0 on the 2PL. The 3PL geometry is the reason -- guessing is weakly identified, so c trades off against d and the sampler is dragged through a funnel that the conjugate Gibbs sampler never enters, since it draws each block from its exact full conditional. A gap of 39 elpd survives that comfortably; a comparison decided by a handful of elpd points would not, and should not be attempted from this fit.
Reading the divergences¶
A divergence is a specific HMC failure, not a generic warning. The sampler simulates a trajectory through parameter space using a discrete step size; where the posterior has a region of high curvature relative to that step, the simulation becomes unstable and its energy is not conserved. PyMC detects that and marks the draw. Divergences matter more than their count suggests because they are not random: they concentrate in exactly the region the sampler is failing to explore, so the affected part of the posterior is under-represented rather than noisily represented.
Here they are informative rather than merely alarming. The 3PL divergences arise because $c_j$ and $d_j$ trade off — a higher guessing floor and a lower intercept produce nearly the same curve — which creates a narrow ridge that steepens where $c_j$ approaches zero. The conjugate Gibbs sampler never meets it, because it draws each block from an exact full conditional and does not simulate trajectories at all. That the 2PL, which has no $c_j$, produces none is the confirmation that the geometry rather than the data is responsible.
What it means for the verdict: an elpd gap of 39 is far larger than the distortion a partially explored ridge could plausibly introduce, so the conclusion stands. A comparison decided by a few elpd points would not survive the same diagnostic, and should not be attempted from a fit in this state.
6. Summary¶
Nested IRT models cannot be chosen by fit alone — the 3PL always fits the training data at least as well as the 1PL. WAIC and PSIS-LOO, computed on the $\theta$-integrated (marginal) likelihood so the person is the unit of prediction, resolve the choice by estimating out-of-sample performance with a complexity penalty. The two datasets gave opposite verdicts from the same criterion: on the near-Rasch LSAT the 1PL won — the discriminations barely differ and there is no guessing, so the extra parameters only add variance; on the multiple-choice SAT12 the 3PL won by more than two standard errors — the lower asymptote captures real guessing near the $1/5$ chance rate. Posterior-predictive ICCs confirmed the winning LSAT model fits in absolute terms, and PyMC reproduced the SAT12 ranking.
The connections: the marginal-likelihood WAIC/LOO is the same Gauss–Hermite scoring used for the latent-class models in Latent Class Analysis — Choosing the Number of Classes — model comparison is model-agnostic. The conditional-vs-marginal distinction is the IRT face of the general cluster/leave-one-group-out issue in hierarchical models. And the three engines here are the 1PL/2PL/3PL from Item Response Theory — 2PL & 3PL, now judged rather than just fit. This capstone closes the IRT arc: 2PL/3PL → polytomous GRM → multidimensional factor model → DIF/explanatory → model comparison — building the models, extending them to ordered and multi-trait data, testing them for bias, and finally choosing among them.