Item Response Theory — 2PL & 3PL¶
Ability, discrimination, guessing, and where a test measures¶
Item response theory models a test one item at a time. Each person $i$ has a continuous latent ability $\theta_i$, and each item $j$ has an item characteristic curve — the probability of a correct answer as a smooth function of ability. In the normal-ogive form, $$\text{2PL: } P(x_{ij}=1\mid\theta_i)=\Phi(a_j\theta_i+d_j),\qquad \text{3PL: } P=c_j+(1-c_j)\Phi(a_j\theta_i+d_j).$$ $a_j$ is the item's discrimination (how sharply it separates ability levels), the difficulty is $b_j=-d_j/a_j$, and $c_j$ is the guessing floor. Making all $a_j$ equal is the Rasch / 1PL model — which is exactly a random-effects logistic regression (person random intercept $=$ ability, item effect $=$ difficulty), the model in the panel-logit notebooks. Here we read it psychometrically and let discrimination and guessing vary.
The from-scratch sampler is Albert & Chib's data augmentation — introduced in Albert (1992) for this very problem before it spread to probit, Tobit and multivariate probit across this portfolio: augment $z_{ij}\sim N(a_j\theta_i+d_j,1)$ truncated by the sign of $x_{ij}$, after which $\theta$ and $(a_j,d_j)$ are conjugate Gaussian. We validate recovery, fit the LSAT (a 2PL near the Rasch boundary), add guessing on the SAT12 multiple-choice test, and read the information functions that say where the test is precise. A PyMC fit and (companion notebook) ltm/mirt confirm it.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy.stats import norm, pointbiserialr
import irt as I
rng = np.random.default_rng(3)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; PURP="#6b46c1"; GREY="#718096"
print("IRT: a continuous ability theta, an item characteristic curve for each item.")
IRT: a continuous ability theta, an item characteristic curve for each item.
1. The item characteristic curve — 1PL vs 2PL¶
The ICC maps ability to the chance of a correct answer. Difficulty $b_j$ slides the curve left/right (harder items need more ability); discrimination $a_j$ sets its steepness (a high-$a$ item sharply distinguishes people near its difficulty). The Rasch/1PL model forces every item to the same steepness — a single random-effects logistic regression — while the 2PL lets discrimination vary.
th=np.linspace(-4,4,200)
fig,ax=plt.subplots(1,2,figsize=(12,4.2))
for b,c_ in [(-1.5,BLUE),(0,GREEN),(1.5,RED)]:
ax[0].plot(th, I.icc(th,1.0,-b*1.0), color=c_, lw=2, label=f"difficulty b={b}")
ax[0].set_title("1PL / Rasch: same discrimination, different difficulty"); ax[0].set_xlabel(r"ability $\theta$"); ax[0].set_ylabel("P(correct)"); ax[0].legend(frameon=False,fontsize=8)
for a_,c_ in [(0.4,GREY),(1.0,GREEN),(2.2,PURP)]:
ax[1].plot(th, I.icc(th,a_,0.0), color=c_, lw=2, label=f"discrimination a={a_}")
ax[1].set_title("2PL: same difficulty, different discrimination"); ax[1].set_xlabel(r"ability $\theta$"); ax[1].set_ylabel("P(correct)"); ax[1].legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print("A steep (high-a) item is a sharp ruler near its difficulty but says little far from it. The Rasch model")
print("assumes every item is equally sharp; the 2PL measures that sharpness -- the extra parameter that separates")
print("IRT from a plain random-effects logit (Bayesian Random-Effects Panel Logit / Bayesian Hierarchical Binary Logit).")
A steep (high-a) item is a sharp ruler near its difficulty but says little far from it. The Rasch model assumes every item is equally sharp; the 2PL measures that sharpness -- the extra parameter that separates IRT from a plain random-effects logit (Bayesian Random-Effects Panel Logit / Bayesian Hierarchical Binary Logit).
2. Does it work? — recovering item parameters¶
Simulated responses of 1500 people to 12 items with known discriminations and difficulties. The Albert–Chib Gibbs sampler should recover both, plus the abilities (up to the sign/scale the $\theta\sim N(0,1)$ prior fixes).
J=12; a_t=rng.uniform(0.5,2.2,J); d_t=rng.uniform(-1.5,1.5,J); b_t=-d_t/a_t
X,th_t=I.simulate_irt(1500,a_t,d_t,rng)
r=I.irt2pl_gibbs(X,rng,draws=1800,burn=900)
a_h=r["a"].mean(0); b_h=-r["d"].mean(0)/a_h
fig,ax=plt.subplots(1,3,figsize=(13,3.9))
ax[0].scatter(a_t,a_h,color=BLUE); ax[0].plot([0,2.5],[0,2.5],"k--",lw=1); ax[0].set_xlabel("true a"); ax[0].set_ylabel("estimated a"); ax[0].set_title(f"discrimination (r={np.corrcoef(a_t,a_h)[0,1]:.2f})")
ax[1].scatter(b_t,b_h,color=GREEN); ax[1].plot([-4,4],[-4,4],"k--",lw=1); ax[1].set_xlabel("true b"); ax[1].set_ylabel("estimated b"); ax[1].set_title(f"difficulty (r={np.corrcoef(b_t,b_h)[0,1]:.2f})")
ax[2].scatter(th_t,r["theta"].mean(0),s=6,color=GREY,alpha=.4); ax[2].plot([-3,3],[-3,3],"k--",lw=1); ax[2].set_xlabel("true θ"); ax[2].set_ylabel("estimated θ"); ax[2].set_title(f"ability (r={np.corrcoef(th_t,r['theta'].mean(0))[0,1]:.2f})")
plt.tight_layout(); plt.show()
print("Discrimination, difficulty and ability all recovered -- the augmentation Gibbs reconstructs the whole item bank")
print("and the people's abilities on one scale.")
Discrimination, difficulty and ability all recovered -- the augmentation Gibbs reconstructs the whole item bank and the people's abilities on one scale.
3. The LSAT — a 2PL at the Rasch boundary, and the classical-test-theory view¶
The five-item Law School Admission Test (1000 examinees). We fit the 2PL and read discrimination and difficulty, then compare with classical test theory (the frequentist psychometric standard that predates IRT): item difficulty as the proportion correct, item discrimination as the point-biserial correlation with total score. IRT reproduces the CTT ordering but places items and people on a common latent scale.
L=pd.read_csv("lsat.csv"); X=L.to_numpy().astype(float); N,J=X.shape
r=I.irt2pl_gibbs(X,rng,draws=3000,burn=1500); a=r["a"].mean(0); d=r["d"].mean(0); b=-d/a
tot=X.sum(1); ctt_disc=np.array([pointbiserialr(X[:,j],tot)[0] for j in range(J)]); ctt_diff=X.mean(0)
tab=pd.DataFrame({"p-correct":ctt_diff.round(2),"IRT difficulty b":b.round(2),"point-biserial":ctt_disc.round(2),"IRT discrimination a":a.round(2)}, index=[f"item {j+1}" for j in range(J)])
print(tab.to_string()); print(f"\ndiscriminations range {a.min():.2f}-{a.max():.2f} -- narrow, so the LSAT is close to Rasch/1PL (equal a).")
th=np.linspace(-4,4,200)
fig,ax=plt.subplots(1,2,figsize=(12,4.2))
for j in range(J): ax[0].plot(th, I.icc(th,a[j],d[j]), lw=2, label=f"item {j+1}")
ax[0].set_xlabel(r"ability $\theta$"); ax[0].set_ylabel("P(correct)"); ax[0].set_title("LSAT item characteristic curves"); ax[0].legend(frameon=False,fontsize=8)
ax[1].scatter(ctt_disc,a,color=PURP); ax[1].set_xlabel("classical point-biserial"); ax[1].set_ylabel("IRT discrimination a"); ax[1].set_title("IRT discrimination vs classical test theory")
plt.tight_layout(); plt.show()
print("The IRT difficulty ranks items exactly as proportion-correct does, and IRT discrimination tracks the point-")
print("biserial -- IRT is the model-based generalisation of classical test theory, with items and abilities on one")
print("scale. The narrow discriminations confirm the LSAT's near-Rasch reputation.")
D = 1.702
print("\nWHICH METRIC. This sampler is NORMAL-OGIVE: P = Phi(a*theta + d). R's ltm and mirt report the")
print("LOGISTIC parameterisation, and the two differ by the scaling constant D = 1.702 -- so the same items")
print("legitimately carry two different numbers, and comparing them without converting invites a phantom")
print("disagreement of about a factor of two. In the logistic metric these discriminations are:")
print(" " + " ".join("item %d: %.2f" % (j+1, a[j]*D) for j in range(J)))
a_R = np.array([0.83, 0.72, 0.89, 0.69, 0.66]) # ltm/mirt, logistic metric
b_R = np.array([-3.36, -1.37, -0.28, -1.87, -3.12])
print("R reports " + " ".join("item %d: %.2f" % (j+1, a_R[j]) for j in range(J)))
print("\nSame ballpark, not the same number: mean absolute gap %.2f on a scale where the values are ~0.7."
% np.abs(a*D - a_R).mean())
print("Rescaling removes the factor-of-two artefact; what is left is a genuine estimator difference, since")
print("these are posterior means under a HalfNormal(2) prior and those are marginal maximum likelihood, on")
print("five items. Difficulty b = -d/a is a RATIO and therefore metric-free, which is why it needs no")
print("conversion at all -- and it agrees far more tightly: correlation %.3f, mean absolute gap %.2f."
% (np.corrcoef(b, b_R)[0,1], np.abs(b - b_R).mean()))
print("That asymmetry is itself the lesson. Difficulty is what a test is FOR and it is robustly estimated;")
print("discrimination is more prior-sensitive and metric-dependent, and deserves more care when quoted.")
p-correct IRT difficulty b point-biserial IRT discrimination a item 1 0.92 -3.67 0.36 0.43 item 2 0.71 -1.32 0.57 0.46 item 3 0.55 -0.29 0.62 0.51 item 4 0.76 -1.86 0.53 0.42 item 5 0.87 -3.27 0.44 0.37 discriminations range 0.37-0.51 -- narrow, so the LSAT is close to Rasch/1PL (equal a).
The IRT difficulty ranks items exactly as proportion-correct does, and IRT discrimination tracks the point- biserial -- IRT is the model-based generalisation of classical test theory, with items and abilities on one scale. The narrow discriminations confirm the LSAT's near-Rasch reputation. WHICH METRIC. This sampler is NORMAL-OGIVE: P = Phi(a*theta + d). R's ltm and mirt report the LOGISTIC parameterisation, and the two differ by the scaling constant D = 1.702 -- so the same items legitimately carry two different numbers, and comparing them without converting invites a phantom disagreement of about a factor of two. In the logistic metric these discriminations are: item 1: 0.73 item 2: 0.78 item 3: 0.87 item 4: 0.71 item 5: 0.63 R reports item 1: 0.83 item 2: 0.72 item 3: 0.89 item 4: 0.69 item 5: 0.66 Same ballpark, not the same number: mean absolute gap 0.05 on a scale where the values are ~0.7. Rescaling removes the factor-of-two artefact; what is left is a genuine estimator difference, since these are posterior means under a HalfNormal(2) prior and those are marginal maximum likelihood, on five items. Difficulty b = -d/a is a RATIO and therefore metric-free, which is why it needs no conversion at all -- and it agrees far more tightly: correlation 0.998, mean absolute gap 0.11. That asymmetry is itself the lesson. Difficulty is what a test is FOR and it is robustly estimated; discrimination is more prior-sensitive and metric-dependent, and deserves more care when quoted.
What these numbers actually say about the test¶
The three item parameters are each on a different scale, and none of them is a probability.
Difficulty $b_j$ lives on the same scale as ability, which is the whole point of IRT: it is the ability at which an examinee has a 50% chance of answering the item correctly. So item 1 at $b=-3.67$ is not merely "easy" — it is answered correctly by anyone within three and a half standard deviations of average ability, which on a normal population is essentially everyone. Item 3 at $b=-0.29$ is the only item on this test pitched near the middle of the group. A test whose difficulties all sit far below the examinees measures almost nobody well, and the information curve later in this notebook is where that shows up.
Discrimination $a_j$ is the slope of the curve at that midpoint: how fast the probability of a correct answer rises as ability increases. Standardised as $a_j/\sqrt{1+a_j^2}$ it is the correlation between the item and the trait, so the LSAT's range of $0.37$–$0.51$ corresponds to item–trait correlations of about $0.35$–$0.45$ — modest, similar across items, and the reason the Rasch model (which assumes they are all equal) is not rejected here. Discrimination is what separates IRT from a plain random-effects logit: the 1PL assumes every item is an equally sharp ruler, and the 2PL measures how sharp each one really is.
Ability $\theta_i$ is fixed to mean 0 and variance 1 by the prior, so it is a z-score within this sample and carries no absolute meaning — an ability of $+1$ means one standard deviation above these 1,000 examinees, not a score out of anything. The scale's zero point and unit are conventions; what is real is the ordering and the distances between people and items, which is exactly what being on one common scale buys.
One consequence worth stating plainly: because people and items share a scale, two examinees who answered different items can still be compared. That is what a total score cannot do, and it is why item response theory underpins adaptive testing and test equating.
4. Guessing — the 3PL on a multiple-choice test¶
On a five-option multiple-choice test even a low-ability examinee gets $\approx$ 1 in 5 by chance, so the ICC should not fall to zero. The 3PL adds a lower asymptote $c_j$, estimated (Béguin–Glas augmentation) from the correct answers of people who did not know the item. We fit the 32-item SAT12 and compare the 3PL and 2PL curves.
S=pd.read_csv("sat12.csv"); XS=S.to_numpy().astype(float)
r3=I.irt3pl_gibbs(XS,rng,draws=2500,burn=1500); a3=r3["a"].mean(0); d3=r3["d"].mean(0); c3=r3["c"].mean(0)
r2=I.irt2pl_gibbs(XS,rng,draws=2500,burn=1500); a2=r2["a"].mean(0); d2=r2["d"].mean(0)
print(f"SAT12 3PL guessing: mean c = {c3.mean():.2f} (a 5-option item guesses at ~0.20), range {c3.min():.2f}-{c3.max():.2f}")
jg=int(np.argmax(c3)) # the item with the most guessing
th=np.linspace(-4,4,200)
fig,ax=plt.subplots(1,2,figsize=(12,4.2))
ax[0].plot(th, I.icc(th,a3[jg],d3[jg],c3[jg]), color=BLUE, lw=2.4, label=f"3PL (c={c3[jg]:.2f})")
ax[0].plot(th, I.icc(th,a2[jg],d2[jg]), color=RED, lw=1.8, ls="--", label="2PL (c=0)")
ax[0].axhline(c3[jg], color=GREY, ls=":", lw=1); ax[0].set_ylim(0,1)
ax[0].set_xlabel(r"ability $\theta$"); ax[0].set_ylabel("P(correct)"); ax[0].set_title(f"SAT12 item {jg+1}: guessing lifts the floor"); ax[0].legend(frameon=False)
ax[1].hist(c3,bins=15,color=PURP,alpha=.7); ax[1].axvline(0.2,color="k",ls="--",lw=1,label="1/5 chance"); ax[1].set_xlabel("estimated guessing c"); ax[1].set_ylabel("items"); ax[1].set_title("Guessing across the 32 items"); ax[1].legend(frameon=False)
plt.tight_layout(); plt.show()
print("The 3PL curve flattens to the guessing floor instead of zero; ignoring it (2PL) forces the low-ability tail")
print("down and distorts difficulty. Most items sit near the 1/5 chance rate a five-option question implies.")
print("\nOn the R side ltm and mirt both put mean guessing near 0.15 against the %.2f here. That gap is not" % c3.mean())
print("noise and it has a specific cause: this is a POSTERIOR MEAN under a Beta(1,4) prior on c, whose own")
print("mean is exactly 0.20. Guessing is the least identified parameter in the 3PL -- it is estimated almost")
print("entirely from low-ability examinees, of whom any test has few -- so the posterior leans on that prior")
print("and sits above the maximum-likelihood estimate. Neither is wrong; a Bayesian 3PL reports where the")
print("prior and the data meet, and on this parameter the prior is doing visible work. That is worth knowing")
print("before quoting a guessing floor to three decimal places.")
SAT12 3PL guessing: mean c = 0.19 (a 5-option item guesses at ~0.20), range 0.08-0.38
The 3PL curve flattens to the guessing floor instead of zero; ignoring it (2PL) forces the low-ability tail down and distorts difficulty. Most items sit near the 1/5 chance rate a five-option question implies. On the R side ltm and mirt both put mean guessing near 0.15 against the 0.19 here. That gap is not noise and it has a specific cause: this is a POSTERIOR MEAN under a Beta(1,4) prior on c, whose own mean is exactly 0.20. Guessing is the least identified parameter in the 3PL -- it is estimated almost entirely from low-ability examinees, of whom any test has few -- so the posterior leans on that prior and sits above the maximum-likelihood estimate. Neither is wrong; a Bayesian 3PL reports where the prior and the data meet, and on this parameter the prior is doing visible work. That is worth knowing before quoting a guessing floor to three decimal places.
Reading the guessing parameter¶
$c_j$ is a lower asymptote, not a probability of guessing: it is where the curve flattens as ability falls, i.e. the chance an examinee with effectively no knowledge still marks the right option. On a five-option item chance alone implies $c \approx 0.20$, and the estimates below cluster there, which is what a credible guessing rate looks like.
Three things follow, and they matter more than the point estimates. A non-zero floor means a correct answer from a low-ability examinee carries less evidence about ability than the 2PL assumes, so modelling guessing honestly lowers the precision the test appears to have. Difficulty also changes meaning: with a floor at $c_j$, the 50% point is no longer $b_j$, since the curve now runs from $c_j$ to 1 rather than 0 to 1. And $c_j$ is the least identified parameter in the model — it is informed almost entirely by the small number of very low-ability examinees any test has — which is why the posterior below leans visibly on its prior, and why classical fits of the same model struggle to converge at all.
5. Information — where the test measures precisely¶
The reason discrimination matters: an item's information $I_j(\theta)=P'(\theta)^2/[P(\theta)(1-P(\theta))]$ is largest near its difficulty and grows with $a_j^2$. Summed over items it is the test information function, and $1/\sqrt{\text{TIF}}$ is the standard error of the ability estimate — so the TIF shows at which abilities the test is sharp. Guessing drains information from the low end.
th=np.linspace(-4,4,200)
fig,ax=plt.subplots(1,2,figsize=(12.5,4.3))
for j in range(min(J,5)): ax[0].plot(th, I.item_information(th,a[j],d[j]), lw=2, label=f"item {j+1}")
TIF,se=I.test_information(th,a,d)
ax[0].plot(th, TIF, color="k", lw=2.5, label="test information"); ax[0].set_xlabel(r"ability $\theta$"); ax[0].set_ylabel("information"); ax[0].set_title("LSAT: item and test information"); ax[0].legend(frameon=False,fontsize=8)
ax2=ax[0].twinx(); ax2.plot(th, se, color=RED, ls=":", lw=1.6); ax2.set_ylabel("SE(θ) = 1/√TIF", color=RED)
# SAT12 3PL vs 2PL total information
T3,_=I.test_information(th,a3,d3,c3); T2,_=I.test_information(th,a2,d2)
ax[1].plot(th,T2,color=RED,lw=2,ls="--",label="ignoring guessing (2PL)"); ax[1].plot(th,T3,color=BLUE,lw=2.3,label="with guessing (3PL)")
ax[1].set_xlabel(r"ability $\theta$"); ax[1].set_ylabel("test information"); ax[1].set_title("SAT12: guessing drains low-end information"); ax[1].legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print("The LSAT's items are easy (all difficulties below 0), so information peaks at somewhat-below-average ability --")
print("that is where SE is smallest -- and the test is blunt at the extremes, especially for high-ability examinees.")
print("On the SAT12 the 3PL shows less information at low ability than the 2PL claims -- a guessed-correct answer")
print("carries little about ability, so honest guessing modelling lowers the measurement precision the test really has.")
The LSAT's items are easy (all difficulties below 0), so information peaks at somewhat-below-average ability -- that is where SE is smallest -- and the test is blunt at the extremes, especially for high-ability examinees. On the SAT12 the 3PL shows less information at low ability than the 2PL claims -- a guessed-correct answer carries little about ability, so honest guessing modelling lowers the measurement precision the test really has.
6. Cross-check in PyMC¶
The same 2PL, sampled directly by NUTS: a standard-normal ability per person, a positive discrimination and a free intercept per item, and a Bernoulli likelihood through the normal-ogive link invprobit. The continuous latent abilities are sampled (nothing to marginalise), so it runs cleanly. We compare discriminations and difficulties with the from-scratch sampler on the LSAT.
import pymc as pm, pytensor.tensor as pt
X=L.to_numpy().astype(float); N,J=X.shape
with pm.Model() as mod:
theta=pm.Normal("theta",0,1,shape=N)
a_=pm.HalfNormal("a",2.0,shape=J); d_=pm.Normal("d",0,3,shape=J)
eta=theta[:,None]*a_[None,:]+d_[None,:]
pm.Bernoulli("x", p=pm.math.invprobit(eta), observed=X)
idata=pm.sample(800,tune=1200,chains=4,target_accept=0.95,random_seed=4,progressbar=False)
import arviz as az
a_pm=idata.posterior["a"].mean(("chain","draw")).values; d_pm=idata.posterior["d"].mean(("chain","draw")).values
rh = az.summary(idata, var_names=["a","d"])["r_hat"]; max_rhat = float(rh.max())
ess = az.summary(idata, var_names=["a","d"])["ess_bulk"]; min_ess = float(ess.min())
print("item a (from-scratch / PyMC) difficulty b (from-scratch / PyMC)")
for j in range(J): print(f" {j+1} {a[j]:.2f} / {a_pm[j]:.2f} {b[j]:.2f} / {(-d_pm[j]/a_pm[j]):.2f}")
print(f"\ndiscrimination agreement: correlation {np.corrcoef(a,a_pm)[0,1]:.3f}")
fig,ax=plt.subplots(figsize=(8,4)); th=np.linspace(-4,4,200)
for j in range(J):
ax.plot(th, I.icc(th,a[j],d[j]), color=BLUE, lw=2); ax.plot(th, I.icc(th,a_pm[j],d_pm[j]), color=RED, lw=1, ls="--")
ax.plot([],[],color=BLUE,lw=2,label="from-scratch"); ax.plot([],[],color=RED,lw=1,ls="--",label="PyMC")
ax.set_xlabel(r"ability $\theta$"); ax.set_ylabel("P(correct)"); ax.set_title("LSAT ICCs: from-scratch vs PyMC"); ax.legend(frameon=False)
plt.tight_layout(); plt.show()
worst = int(np.argmax(np.abs(a - a_pm)))
print("\nHow good is that agreement really? Max R-hat on the item parameters is %.2f and the smallest bulk" % max_rhat)
print("ESS is %.0f. Both are past the thresholds PyMC warns at, so this chain is not cleanly converged and" % min_ess)
print("some of the gap below is sampling noise rather than a difference of model. The worst item is %d:"
% (worst+1))
print("a = %.2f from the Gibbs sampler against %.2f from PyMC, a %.0f%% gap."
% (a[worst], a_pm[worst], 100*abs(a[worst]-a_pm[worst])/a[worst]))
print("Difficulty, again, holds up much better (correlation %.3f) than discrimination (%.3f)."
% (np.corrcoef(b, -d_pm/a_pm)[0,1], np.corrcoef(a, a_pm)[0,1]))
print("\nWhy NUTS struggles is worth naming: this model has 1000 ability parameters and only 5 items, the")
print("scale is fixed solely by the theta ~ N(0,1) prior, and the sign is not identified at all -- nothing")
print("stops a chain exploring (-a, -theta). The Albert-Chib Gibbs sampler sidesteps all of it by drawing")
print("from exact conjugate full conditionals and pinning the sign each sweep. That is the practical case")
print("for the augmentation scheme: on this geometry the purpose-built sampler beats general-purpose HMC.")
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 1_200 tune and 800 draw iterations (4_800 + 3_200 draws total) took 13 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
item a (from-scratch / PyMC) difficulty b (from-scratch / PyMC) 1 0.43 / 0.43 -3.67 / -3.64 2 0.46 / 0.43 -1.32 / -1.39 3 0.51 / 0.55 -0.29 / -0.28 4 0.42 / 0.40 -1.86 / -1.93 5 0.37 / 0.36 -3.27 / -3.32 discrimination agreement: correlation 0.957
How good is that agreement really? Max R-hat on the item parameters is 1.03 and the smallest bulk ESS is 349. Both are past the thresholds PyMC warns at, so this chain is not cleanly converged and some of the gap below is sampling noise rather than a difference of model. The worst item is 3: a = 0.51 from the Gibbs sampler against 0.55 from PyMC, a 8% gap. Difficulty, again, holds up much better (correlation 0.999) than discrimination (0.957). Why NUTS struggles is worth naming: this model has 1000 ability parameters and only 5 items, the scale is fixed solely by the theta ~ N(0,1) prior, and the sign is not identified at all -- nothing stops a chain exploring (-a, -theta). The Albert-Chib Gibbs sampler sidesteps all of it by drawing from exact conjugate full conditionals and pinning the sign each sweep. That is the practical case for the augmentation scheme: on this geometry the purpose-built sampler beats general-purpose HMC.
How stable is that agreement number?¶
The correlation just printed is itself a statistic computed from a chain that has not converged, so it inherits the same problem. Recomputing it from subsets of the very same four chains shows how much of it is sampling noise: if two chains would have given a visibly different answer, the number is not reporting agreement between two models, it is reporting which chains happened to run.
import itertools
a_by_chain = idata.posterior["a"].mean("draw").values # (chain, J)
print("discrimination correlation with the Gibbs fit, by which of the 4 chains are pooled:\n")
spread = {}
for k in (1, 2, 3, 4):
vals = [np.corrcoef(a, a_by_chain[list(s)].mean(0))[0,1] for s in itertools.combinations(range(4), k)]
spread[k] = (min(vals), max(vals), float(np.mean(vals)))
print(" %d chain(s), %2d subset(s): %.3f to %.3f (mean %.3f)" % (k, len(vals), *spread[k]))
lo, hi = spread[2][0], spread[2][1]
print("\nA two-chain run of this model could have reported anywhere from %.3f to %.3f for the same" % (lo, hi))
print("data and the same two samplers -- a spread of %.3f, which is the same order as the distance" % (hi-lo))
print("between the pooled figure of %.3f and a perfect 1.0. So the agreement statistic carries real" % spread[4][2])
print("sampling error of its own, and quoting it to three decimals overstates what it establishes.")
print("This is worth doing whenever a cross-check is summarised by a single correlation: recompute it")
print("from subsets before treating the last two digits as a finding. Difficulty, being a ratio and far")
print("better behaved, is the parameter on which the two engines can actually be compared.")
discrimination correlation with the Gibbs fit, by which of the 4 chains are pooled: 1 chain(s), 4 subset(s): 0.929 to 0.968 (mean 0.950) 2 chain(s), 6 subset(s): 0.941 to 0.970 (mean 0.956) 3 chain(s), 4 subset(s): 0.951 to 0.963 (mean 0.957) 4 chain(s), 1 subset(s): 0.957 to 0.957 (mean 0.957) A two-chain run of this model could have reported anywhere from 0.941 to 0.970 for the same data and the same two samplers -- a spread of 0.029, which is the same order as the distance between the pooled figure of 0.957 and a perfect 1.0. So the agreement statistic carries real sampling error of its own, and quoting it to three decimals overstates what it establishes. This is worth doing whenever a cross-check is summarised by a single correlation: recompute it from subsets before treating the last two digits as a finding. Difficulty, being a ratio and far better behaved, is the parameter on which the two engines can actually be compared.
7. Summary¶
Item response theory puts people and test items on one latent scale. The 2PL gives each item a difficulty and a discrimination; the 3PL adds a guessing floor for multiple-choice questions. From scratch it is Albert & Chib augmentation — the trick that Albert (1992) invented for this model and that recurs across the portfolio (Tobit, multivariate probit, probit LCA) — with $\theta$ and item parameters conjugate given the augmented data. We recovered a known item bank, found the LSAT sitting near the Rasch boundary (narrow discriminations), estimated SAT12 guessing near the $\approx 1/5$ rate a five-option item implies, and read the information functions that say where a test measures precisely and how guessing erodes that precision.
The connections run throughout: the Rasch/1PL special case is the random-effects binary logit of Bayesian Random-Effects Panel Logit and Bayesian Hierarchical Binary Logit, read psychometrically; IRT is the continuous-trait sibling of latent class analysis (discrete class) in the Bartholomew latent-variable taxonomy, and the Conditional-Dependence LCA notebook already put a continuous trait inside the classes; and IRT generalises classical test theory's proportion-correct and point-biserial into a probability model. The frequentist marginal-ML fits (ltm, mirt) are in the companion R notebook. Next in the arc: polytomous IRT (graded-response and partial-credit models) for ordinal items.