Polytomous IRT — the Graded Response Model¶
Ordered Likert items: a curve for every category¶
The 2PL handles right/wrong items. Survey and attitude items are usually ordered categories — strongly disagree … strongly agree — and each category needs its own curve. Samejima's graded response model (1969) is the ordered-probit item response model: a person with ability $\theta_i$ produces a latent propensity $z_{ij}=a_j\theta_i+\text{noise}$, and the observed category is the interval $z$ falls in, $$x_{ij}=k \iff \gamma_{j,k-1}<z_{ij}\le\gamma_{j,k},\qquad z_{ij}\sim N(a_j\theta_i,1),$$ with item discrimination $a_j$ and ordered thresholds $\gamma_{j,1}<\dots<\gamma_{j,K-1}$. The cumulative curve is $P(x_{ij}\ge k)=\Phi(a_j\theta_i-\gamma_{j,k-1})$, and differencing adjacent cumulatives gives the category response curves. With $K=2$ this is the 2PL.
From scratch it is the ordered-probit data augmentation (Albert & Chib 1993) — the same truncated-normal trick as the dichotomous IRT and the ordinal/sequential-probit models: draw $z_{ij}$ truncated to the observed category's interval, then $\theta$ and $a$ are conjugate Gaussian and the thresholds have order-constrained conditionals. We validate recovery, fit the Neuroticism scale (five six-point Big-Five personality items), read the category curves and information, and cross-check in PyMC with an ordered-logistic model. Higher scores mean more neurotic.
import numpy as np, pandas as pd, matplotlib.pyplot as plt
import grm as G
rng = np.random.default_rng(4)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; PURP="#6b46c1"; GREY="#718096"
CAT=[BLUE,GREEN,ORANGE,RED]
print("Graded response model: an ordered-probit item response model with a curve per category.")
Graded response model: an ordered-probit item response model with a curve per category.
1. Category response curves¶
For a four-category item the ability axis is cut into four ordered regions by three thresholds. Low ability makes the lowest category most likely, high ability the highest, and the middle categories peak in between. The cumulative curves $P(x\ge k)$ are parallel ogives (shifted by the thresholds, steepness set by discrimination); their differences are the category curves.
th=np.linspace(-4,4,300); a=1.4; gam=np.array([-1.2,0.0,1.5])
P=G.category_probs(th,a,gam)
fig,ax=plt.subplots(1,2,figsize=(12,4.2))
for k in range(4): ax[0].plot(th, P[:,k], color=CAT[k], lw=2.2, label=f"category {k+1}")
ax[0].set_title("Category response curves (4-point item)"); ax[0].set_xlabel(r"ability $\theta$"); ax[0].set_ylabel("probability")
ax[0].legend(frameon=False, fontsize=8)
from scipy.special import ndtr as Phi
for k in range(1,4): ax[1].plot(th, Phi(a*th-gam[k-1]), color=CAT[k], lw=2.2, label=f"P(x>={k+1})")
ax[1].set_title("Cumulative curves P(x >= k) -- parallel ogives"); ax[1].set_xlabel(r"ability $\theta$"); ax[1].set_ylabel("probability"); ax[1].legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print("Each category owns a stretch of the ability scale. Discrimination sets how sharply the item separates adjacent")
print("levels; the thresholds place the boundaries. This is the 2PL generalised from two categories to K.")
Each category owns a stretch of the ability scale. Discrimination sets how sharply the item separates adjacent levels; the thresholds place the boundaries. This is the 2PL generalised from two categories to K.
2. Does it work? — recovering discriminations and thresholds¶
Simulated responses of 1500 people to 7 four-category items with known discriminations and thresholds. The ordered-probit augmentation should recover both, plus the abilities.
J,K=7,4; a_t=rng.uniform(0.7,2.0,J); gam_t=np.sort(rng.uniform(-2,2,(J,K-1)),axis=1)
X,th_t=G.simulate_grm(1500,a_t,gam_t,rng)
r=G.grm_gibbs(X,rng,draws=1800,burn=900); a_h=r["a"].mean(0); gam_h=r["gamma"].mean(0)
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(gam_t.ravel(),gam_h.ravel(),color=GREEN); ax[1].plot([-3,3],[-3,3],"k--",lw=1); ax[1].set_xlabel("true threshold"); ax[1].set_ylabel("estimated"); ax[1].set_title(f"thresholds (RMSE={np.sqrt(np.mean((gam_t-gam_h)**2)):.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("Discriminations, thresholds and abilities all recovered -- the ordered-probit augmentation reconstructs the")
print("full polytomous item bank.")
Discriminations, thresholds and abilities all recovered -- the ordered-probit augmentation reconstructs the full polytomous item bank.
3. A Neuroticism scale — five items, one trait¶
2694 respondents rate five Big-Five neuroticism statements — angers easily, irritated easily, has mood swings, feels blue, panics easily — on a six-point scale (higher = more neurotic). All five discriminate; the graded model shows how sharply each separates neuroticism levels.
d=pd.read_csv("neuro.csv"); X=d.to_numpy(); Kc=int(X.max()); _lab={"N1":"angers easily","N2":"irritated easily","N3":"mood swings","N4":"feels blue","N5":"panics easily"}; items=[_lab.get(c,c) for c in d.columns]
# Chain length matters a great deal here. The order-constrained threshold update draws each
# cutpoint between max(z | x=c) and min(z | x=c+1); at N=2694 those order statistics sit O(1/N)
# apart, so the cutpoints crawl and the sharpest items' discriminations are still climbing at
# 3000 draws. Section 6 below fits this same data at both lengths and diagnoses the difference
# against simulated data with known parameters. So we run the longer chain.
r=G.grm_gibbs(X,rng,draws=12000,burn=6000); a=r["a"].mean(0); gam=r["gamma"].mean(0)
tab=pd.DataFrame({"discrimination a":a.round(2)}, index=items).sort_values("discrimination a",ascending=False)
print(tab.to_string()); print(f"\nall five items discriminate on neuroticism; anger and irritation are the sharpest indicators, panic the bluntest.")
strong=int(np.argmax(a)); weak=int(np.argmin(a))
th=np.linspace(-4,4,200)
fig,ax=plt.subplots(1,2,figsize=(12,4.3))
for k in range(Kc): ax[0].plot(th, G.category_probs(th,a[strong],gam[strong])[:,k], lw=2.0, label=f"cat {k+1}")
ax[0].set_title(f"Sharpest item: {items[strong]} (a={a[strong]:.2f})"); ax[0].set_xlabel(r"neuroticism $\theta$"); ax[0].set_ylabel("P(category)"); ax[0].legend(frameon=False,fontsize=7,ncol=2)
for k in range(Kc): ax[1].plot(th, G.category_probs(th,a[weak],gam[weak])[:,k], lw=2.0)
ax[1].set_title(f"Bluntest item: {items[weak]} (a={a[weak]:.2f})"); ax[1].set_xlabel(r"neuroticism $\theta$"); ax[1].set_ylabel("P(category)")
plt.tight_layout(); plt.show()
print("The sharp item's six category curves are well separated -- each response marks a distinct neuroticism level;")
print("the blunt item overlaps more.")
D = 1.702
a_mirt = np.array([3.13, 2.89, 2.03, 1.28, 1.12]) # mirt/ltm, logistic metric, same item order
ratio = a_mirt / (a * D)
print("\nAgainst the R packages (mirt and ltm agree with each other to three decimals), converted to the")
print("same logistic metric:")
for nm, ours, theirs, rt in zip(items, a*D, a_mirt, ratio):
print(" %-18s ours %.2f mirt %.2f ratio %.2f" % (nm, ours, theirs, rt))
print("\nThe ORDERING is identical and the three blunt items match closely, but the two sharpest sit")
print("below mirt by around %.0f%%. That gap is not a modelling difference, it is residual Monte Carlo" % (100*(ratio[:2].mean()-1)))
print("bias: the order-constrained threshold update is slow at N=%d, and the items whose thresholds are" % len(X))
print("most spread out are the last to settle. Section 5 puts numbers on both halves of that claim:")
print("the same data fitted at 3,000 draws against 12,000, and a simulation at these dimensions with")
print("known parameters. Worth knowing before trusting a discrimination from a short run: this sampler")
print("is exact in the limit and visibly biased before it.")
discrimination a angers easily 1.74 irritated easily 1.52 mood swings 1.14 feels blue 0.72 panics easily 0.63 all five items discriminate on neuroticism; anger and irritation are the sharpest indicators, panic the bluntest.
The sharp item's six category curves are well separated -- each response marks a distinct neuroticism level; the blunt item overlaps more. Against the R packages (mirt and ltm agree with each other to three decimals), converted to the same logistic metric: angers easily ours 2.96 mirt 3.13 ratio 1.06 irritated easily ours 2.58 mirt 2.89 ratio 1.12 mood swings ours 1.95 mirt 2.03 ratio 1.04 feels blue ours 1.23 mirt 1.28 ratio 1.04 panics easily ours 1.07 mirt 1.12 ratio 1.04 The ORDERING is identical and the three blunt items match closely, but the two sharpest sit below mirt by around 9%. That gap is not a modelling difference, it is residual Monte Carlo bias: the order-constrained threshold update is slow at N=2694, and the items whose thresholds are most spread out are the last to settle. Section 5 puts numbers on both halves of that claim: the same data fitted at 3,000 draws against 12,000, and a simulation at these dimensions with known parameters. Worth knowing before trusting a discrimination from a short run: this sampler is exact in the limit and visibly biased before it.
Reading a discrimination and a threshold¶
Both numbers are on the latent-response scale, and neither is a probability, so it is worth saying what they amount to.
A discrimination $a_j$ is the slope of the item's latent response on the trait. Because the residual has unit variance, $a_j$ converts directly into the correlation between the item and the trait, $\lambda_j = a_j/\sqrt{1+a_j^2}$ — the same standardisation a factor analysis would print as a loading. So angers easily at $a=1.74$ correlates about 0.87 with neuroticism, while panics easily at $a=0.63$ correlates about 0.53. Both are real indicators; one is roughly twice as informative per response as the other, which is what "sharp" and "blunt" mean here. A useful reference point: $a=0$ is an item that tells you nothing, and $a\to\infty$ is a perfect deterministic indicator that nobody ever writes.
A threshold $\gamma_{jc}$ is a cutpoint on that same latent scale. Dividing by the discrimination puts it back on the trait scale, $b_{jc} = \gamma_{jc}/a_j$, where it reads directly: $b_{jc}$ is the neuroticism level at which a respondent becomes more likely than not to answer above category $c$. Five thresholds carve the trait axis into six regions, one per response option, and the spacing between them is what makes the category curves in the figures wide or narrow. An item whose thresholds are bunched together distinguishes finely in a narrow band of the trait and poorly everywhere else; an item whose thresholds are spread out grades the whole range coarsely.
That distinction matters for the section below: the spread-out items are the ones whose thresholds the sampler is slowest to place.
All five items contribute to the test information, but because information grows with $a_j^2$ the sharpest items (anger, irritation) dominate. The information function shows the range of neuroticism the scale pins down most precisely.
th=np.linspace(-4,4,200); TIF=G.test_information(th,a,gam)
order=np.argsort(-a)
fig,ax=plt.subplots(figsize=(8.5,4.2))
for j in order: ax.plot(th, G.item_information(th,a[j],gam[j]), lw=1.6, label=f"{items[j]} (a={a[j]:.2f})")
ax.plot(th, TIF, color="k", lw=2.6, label="test information")
ax.set_xlabel(r"neuroticism $\theta$"); ax.set_ylabel("information"); ax.set_title("Item and test information (Neuroticism scale)"); ax.legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print(f"the two sharpest items supply {G.test_information(th,a[order[:2]],gam[order[:2]]).max()/TIF.max():.0%} of the peak information; the scale is most precise near theta={th[TIF.argmax()]:+.1f}.")
the two sharpest items supply 70% of the peak information; the scale is most precise near theta=+0.2.
4. Cross-check in PyMC — an ordered-logistic graded model¶
The graded model in the logistic metric: each item is an ordered logistic with a discrimination-scaled linear predictor $a_j\theta_i$ and its own ordered cutpoints. NUTS samples the continuous abilities and item parameters directly. We compare discriminations (normal-ogive $a\times1.7$ maps to the logistic metric) with the from-scratch sampler.
import pymc as pm, pytensor.tensor as pt
from pymc.distributions.transforms import ordered
N,J=X.shape; K=int(X.max()); Xi=X-1
with pm.Model() as mod:
theta=pm.Normal("theta",0,1,shape=N)
a_=pm.HalfNormal("a",2.0,shape=J)
for j in range(J):
cut=pm.Normal(f"c{j}", 0, 3, shape=K-1, transform=ordered, initval=np.linspace(-2,2,K-1))
pm.OrderedLogistic(f"x{j}", eta=a_[j]*theta, cutpoints=cut, observed=Xi[:,j])
idata=pm.sample(700, tune=1200, chains=4, target_accept=0.9, random_seed=5, progressbar=False)
a_pm=idata.posterior["a"].mean(("chain","draw")).values
cmp=pd.DataFrame({"from-scratch a x1.7":(a*1.7).round(2),"PyMC a (logistic)":a_pm.round(2)}, index=items)
print(cmp.to_string()); print(f"\ndiscrimination-ordering agreement: correlation {np.corrcoef(a*1.7,a_pm)[0,1]:.3f}")
fig,ax=plt.subplots(figsize=(7.5,4)); ax.scatter(a*1.7,a_pm,color=BLUE)
for i,it in enumerate(items): ax.annotate(it,(a[i]*1.7,a_pm[i]),fontsize=7)
ax.plot([0,3],[0,3],"k--",lw=1); ax.set_xlabel("from-scratch a x1.7"); ax.set_ylabel("PyMC ordered-logistic a"); ax.set_title("Discrimination: from-scratch vs PyMC")
plt.tight_layout(); plt.show()
print("Same three strong items, same near-zero weak ones: the ordered-probit Gibbs and the PyMC ordered-logistic")
print("model agree on which statements measure the neuroticism trait.")
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, c0, c1, c2, c3, c4]
Sampling 4 chains for 1_200 tune and 700 draw iterations (4_800 + 2_800 draws total) took 87 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
from-scratch a x1.7 PyMC a (logistic) angers easily 2.95 3.12 irritated easily 2.58 2.89 mood swings 1.94 2.04 feels blue 1.23 1.28 panics easily 1.07 1.12 discrimination-ordering agreement: correlation 0.997
Same three strong items, same near-zero weak ones: the ordered-probit Gibbs and the PyMC ordered-logistic model agree on which statements measure the neuroticism trait.
5. Chain length, diagnosed¶
The discriminations above came from a 12,000-draw chain, and section 3 claimed that length mattered. This section stops asserting that and measures it: the same neuroticism data fitted at 3,000 draws as well, and then a simulation at the data's exact dimensions — 2,694 respondents, 5 items, 6 categories — where the true discriminations are known and recovery can be checked directly.
For the simulation the "truth" is mirt's own fit converted to the normal-ogive metric ($a/D$), with mirt's thresholds placed on the latent scale as $\gamma_{jc} = a_j b_{jc}$. Five independently simulated datasets are used, so the conclusion does not rest on one lucky seed.
# The same data at the shorter chain -- the length that looked adequate.
r_short = G.grm_gibbs(X, np.random.default_rng(21), draws=3000, burn=1500)
a_short = r_short["a"].mean(0)
cmp3 = pd.DataFrame({"3,000 draws": (a_short*D).round(2),
"12,000 draws": (a*D).round(2),
"mirt": a_mirt.round(2)}, index=items)
cmp3["gap before"] = [f"{v:.0f}%" for v in (a_mirt - a_short*D)/a_mirt*100]
cmp3["gap after"] = [f"{v:.0f}%" for v in (a_mirt - a*D)/a_mirt*100]
print("discrimination, logistic metric -- short chain against long:\n")
print(cmp3.to_string())
sharp = np.argsort(-a_mirt)[:2]; blunt = np.argsort(-a_mirt)[2:]
print("\nThe gap is patterned, not random: at 3,000 draws the two sharpest items sit %.0f%% and %.0f%%"
% tuple((a_mirt[sharp] - a_short[sharp]*D)/a_mirt[sharp]*100))
print("below mirt while the three blunt ones are within %.0f%%. Lengthening the chain moves the sharp"
% max(abs((a_mirt[blunt] - a_short[blunt]*D)/a_mirt[blunt]*100)))
print("pair to %.0f%% and %.0f%%. The blunt three barely move -- they sit within a few percent of mirt"
% tuple((a_mirt[sharp] - a[sharp]*D)/a_mirt[sharp]*100))
print("on either side at both chain lengths, which is what having already converged looks like.")
discrimination, logistic metric -- short chain against long:
3,000 draws 12,000 draws mirt gap before gap after
angers easily 2.63 2.96 3.13 16% 5%
irritated easily 2.33 2.58 2.89 19% 11%
mood swings 2.08 1.95 2.03 -3% 4%
feels blue 1.32 1.23 1.28 -3% 4%
panics easily 1.12 1.07 1.12 -0% 4%
The gap is patterned, not random: at 3,000 draws the two sharpest items sit 16% and 19%
below mirt while the three blunt ones are within 3%. Lengthening the chain moves the sharp
pair to 5% and 11%. The blunt three barely move -- they sit within a few percent of mirt
on either side at both chain lengths, which is what having already converged looks like.
# Simulate at the real data's exact dimensions with KNOWN parameters, and fit at both lengths.
b_mirt = np.array([[-0.82,-0.10, 0.33, 0.97, 1.70], # mirt thresholds, same item order
[-1.37,-0.56,-0.12, 0.64, 1.47], # (from the R notebook's graded fit)
[-1.19,-0.30, 0.11, 0.87, 1.76],
[-1.57,-0.37, 0.23, 1.21, 2.25],
[-1.30,-0.13, 0.48, 1.45, 2.51]])
a_true = a_mirt / D # mirt's discriminations in the normal-ogive metric
gam_true = a_true[:,None] * b_mirt # thresholds on the latent-response scale
SEEDS = [1,2,3,4,5]
acc = {3000: [], 12000: []}
for s in SEEDS:
Xs,_ = G.simulate_grm(len(X), a_true, gam_true, np.random.default_rng(100+s))
for dr in (3000, 12000):
rr = G.grm_gibbs(Xs, np.random.default_rng(200+s), draws=dr, burn=dr//2)
acc[dr].append((rr["a"].mean(0), np.sqrt(np.mean((rr["gamma"].mean(0)-gam_true)**2))))
print("simulation at N=%d, J=%d, K=%d over %d datasets\n" % (len(X), X.shape[1], Kc, len(SEEDS)))
rows = {}
for dr in (3000, 12000):
A = np.array([x[0] for x in acc[dr]]); RM = np.array([x[1] for x in acc[dr]])
rows[dr] = (A.mean(0), (A.mean(0)-a_true)/a_true*100, RM.mean(), RM.min(), RM.max())
sim = pd.DataFrame({"true a": a_true.round(2),
"3,000 draws": rows[3000][0].round(2),
"bias": [f"{v:+.0f}%" for v in rows[3000][1]],
"12,000 draws": rows[12000][0].round(2),
"bias ": [f"{v:+.0f}%" for v in rows[12000][1]]}, index=items)
print(sim.to_string())
print("\nthreshold RMSE: 3,000 draws %.2f (range %.2f-%.2f) 12,000 draws %.2f (range %.2f-%.2f)"
% (rows[3000][2], rows[3000][3], rows[3000][4], rows[12000][2], rows[12000][3], rows[12000][4]))
print("\nThe short chain is biased low on exactly the two sharpest items (%.0f%% and %.0f%%) while the"
% tuple(rows[3000][1][np.argsort(-a_true)[:2]]))
print("three blunt ones are already within %.0f%%. Four times the chain removes it: every item lands"
% max(abs(rows[3000][1][np.argsort(-a_true)[2:]])))
print("within %.0f%% of truth and the threshold RMSE falls by a factor of %.1f. So the real-data gap"
% (max(abs(rows[12000][1])), rows[3000][2]/rows[12000][2]))
print("above was chain length, not the model -- and the items it hits are the ones whose thresholds")
print("are most spread out, which is exactly what the O(1/N) order-statistic argument predicts.")
simulation at N=2694, J=5, K=6 over 5 datasets
true a 3,000 draws bias 12,000 draws bias
angers easily 1.84 1.76 -4% 1.87 +2%
irritated easily 1.70 1.51 -11% 1.71 +1%
mood swings 1.19 1.20 +1% 1.19 -0%
feels blue 0.75 0.77 +3% 0.76 +1%
panics easily 0.66 0.68 +3% 0.66 -0%
threshold RMSE: 3,000 draws 0.20 (range 0.18-0.21) 12,000 draws 0.08 (range 0.05-0.10)
The short chain is biased low on exactly the two sharpest items (-4% and -11%) while the
three blunt ones are already within 3%. Four times the chain removes it: every item lands
within 2% of truth and the threshold RMSE falls by a factor of 2.5. So the real-data gap
above was chain length, not the model -- and the items it hits are the ones whose thresholds
are most spread out, which is exactly what the O(1/N) order-statistic argument predicts.
# --- Fig: the chain-length diagnosis, from the runs above ---
A3 = np.array([x[0] for x in acc[3000]]); A12 = np.array([x[0] for x in acc[12000]])
R3 = np.array([x[1] for x in acc[3000]]); R12 = np.array([x[1] for x in acc[12000]])
fig, (axL, axR) = plt.subplots(1, 2, figsize=(12.4, 4.6))
# left: recovery against known truth, both chain lengths, error bars over seeds
lim = [0.5, 2.1]
axL.plot(lim, lim, "k--", lw=1, zorder=1)
axL.errorbar(a_true, A3.mean(0), yerr=A3.std(0), fmt="o", ms=7, color=RED, capsize=3,
label="3,000 draws", zorder=3)
axL.errorbar(a_true, A12.mean(0), yerr=A12.std(0), fmt="s", ms=6, color=BLUE, capsize=3,
label="12,000 draws", zorder=3)
for t, v in zip(a_true, A3.mean(0)):
if (t - v) / t > 0.03:
axL.annotate("", xy=(t, v), xytext=(t, t), zorder=2,
arrowprops=dict(arrowstyle="-", color=RED, lw=.9, alpha=.6))
axL.set_xlim(lim); axL.set_ylim(lim)
axL.set_xlabel("true discrimination"); axL.set_ylabel("posterior mean")
axL.set_title("Recovery at known truth (%d simulated datasets)" % len(SEEDS))
axL.legend(frameon=False, fontsize=8, loc="upper left"); axL.grid(alpha=.25)
axL.text(0.97, 0.05, "threshold RMSE\n3,000: %.2f\n12,000: %.2f" % (R3.mean(), R12.mean()),
transform=axL.transAxes, ha="right", va="bottom", fontsize=8,
bbox=dict(boxstyle="round,pad=0.4", fc="white", ec="#cbd5e0"))
# right: the real-data gap against mirt, per item, at both lengths
o = np.argsort(-a_mirt)
gap3 = (a_mirt - a_short*D)/a_mirt*100
gap12 = (a_mirt - a*D)/a_mirt*100
y = np.arange(len(o))
axR.axvline(0, color="0.5", lw=1)
for i, j in enumerate(o):
axR.plot([gap3[j], gap12[j]], [i, i], color="0.75", lw=1.4, zorder=1)
axR.scatter(gap3[o], y, s=55, color=RED, zorder=3, label="3,000 draws")
axR.scatter(gap12[o], y, s=45, color=BLUE, marker="s", zorder=3, label="12,000 draws")
axR.set_yticks(y); axR.set_yticklabels([items[j] for j in o], fontsize=8.5)
axR.invert_yaxis()
axR.set_xlabel("gap below mirt (%)")
axR.set_title("Neuroticism scale: the gap closes only where it was open")
axR.legend(frameon=False, fontsize=8, loc="lower right"); axR.grid(alpha=.25, axis="x")
plt.tight_layout(); plt.show()
print("The left panel is the diagnosis: at 3,000 draws the two sharpest items sit visibly below the")
print("identity line and the three blunt ones are already on it, so the bias is a property of the")
print("items' threshold spread rather than of the model. At 12,000 every point is on the line and the")
print("threshold RMSE has fallen by a factor of %.1f. The right panel shows the same thing on the real" % (R3.mean()/R12.mean()))
print("data, where there is no truth to check against. Read it by distance from zero rather than by")
print("how far each point travels: the two sharp items start %.0f%% and %.0f%% outside and come in, while"
% tuple(np.sort(gap3)[::-1][:2]))
print("the blunt three sit within %.0f%% of mirt at BOTH lengths, on one side or the other. The chain"
% max(abs(np.r_[gap3[o[2:]], gap12[o[2:]]])))
print("length changes the answer only where the answer had not yet settled.")
The left panel is the diagnosis: at 3,000 draws the two sharpest items sit visibly below the identity line and the three blunt ones are already on it, so the bias is a property of the items' threshold spread rather than of the model. At 12,000 every point is on the line and the threshold RMSE has fallen by a factor of 2.5. The right panel shows the same thing on the real data, where there is no truth to check against. Read it by distance from zero rather than by how far each point travels: the two sharp items start 19% and 16% outside and come in, while the blunt three sit within 4% of mirt at BOTH lengths, on one side or the other. The chain length changes the answer only where the answer had not yet settled.
6. Summary¶
The graded response model gives every ordered category its own curve: a discrimination and a set of thresholds turn a person's latent attitude into a distribution over the response options. From scratch it is ordered-probit data augmentation — the same Albert–Chib truncated-normal machinery as the dichotomous IRT and the ordinal/sequential-probit models — recovering discriminations, thresholds and abilities from simulated data. On a five-item Neuroticism scale all items discriminated on the trait — anger and irritation the sharpest, panic the bluntest — carrying the bulk of the test information, exactly the pattern ltm and mirt report (they agree to three decimals), confirmed by a PyMC ordered-logistic fit.
The connections: with two categories the graded model is the 2PL; its ordered-probit engine is the same one behind Bayesian Sequential Probit (ordinal outcomes), now with a latent ability; and the Rasch-family alternative — the Partial Credit and Generalised Partial Credit models (adjacent-category rather than cumulative logits) — is fitted with mirt in the companion R notebook. Next in the arc: multidimensional IRT and the factor-analysis bridge.