Item Response Theory — 2PL & 3PL

Python · PyMC · R  ·  Download IRT module

A Curve for Every Item

Item response theory models a test rather than a score. Each person carries a latent ability θ\theta, and each item has a curve giving the probability of a correct answer as a smooth function of that ability. The 2PL gives every item a difficulty and a discrimination — how sharply it separates ability levels — and the 3PL adds a guessing floor for multiple-choice items. Constraining all discriminations equal recovers the Rasch/1PL model, which is why this connects directly to the random-effects logit elsewhere in the collection; here it is read psychometrically.

P(xij=1∣θi)=cj+(1−cj) Φ ⁣(ajθi+dj),bj=−dj/ajP(x_{ij}=1\mid\theta_i)=c_j+(1-c_j)\,\Phi\!\big(a_j\theta_i+d_j\big),\qquad b_j=-d_j/a_j

The sampler, and where it came from

The from-scratch sampler is Albert & Chib data augmentation — and worth noting, the 1992 paper that introduced the trick did so for IRT, before it spread to probit, Tobit and multivariate probit. Introducing a truncated normal latent variable makes both ability and item parameters conjugate Gaussian draws. The 3PL extends it with a latent "knows-it" indicator (Béguin & Glas): a correct answer from a non-knower is a guess, so the guessing rate gets a Beta posterior and the augmentation is built on knowing rather than on the response.

The LSAT, at the Rasch boundary

On the LSAT — 1,000 examinees, 5 items — the model reproduces classical test theory and then goes past it. IRT difficulty ranks the items exactly as proportion-correct does and discrimination tracks the point-biserial, but now items and people sit on one common scale. The discriminations span a narrow 0.37–0.51, which is the LSAT's near-Rasch reputation showing up as a number; R's likelihood-ratio test confirms it, at LR = 0.57 on 4 df, p=0.97p=0.97 — the 2PL buys nothing over Rasch here.

Two Engines, Two Rulers

The two engines appear to disagree on discrimination by a factor of two — 0.37–0.51 in Python against 0.66–0.89 in R — and the reason is not that either is wrong. They are in different metrics. The from-scratch sampler is normal-ogive (probit); ltm and mirt report the logistic parameterisation, and the two differ by the scaling constant D=1.702D=1.702. The item-wise ratios come out at 1.93, 1.57, 1.75, 1.64, 1.78 — mean 1.73. It is the kind of constant that is easy to carry in a docstring and forget to apply.

LSAT item12345agreement
discrimination, normal-ogive (Python)0.430.460.510.420.37mean gap 0.05
after rescaling
  × D=1.702D=1.702 → logistic0.730.780.870.710.63
discrimination, logistic (R)0.830.720.890.690.66
difficulty bb (Python) — metric-free−3.67−1.32−0.29−1.86−3.27correlation 0.998
no conversion needed
difficulty bb (R)−3.36−1.37−0.28−1.87−3.12

Rescaled, the mean absolute gap is 0.05 on values around 0.7 — the same ballpark, not the same number, because one side is a posterior mean under a HalfNormal prior and the other is marginal maximum likelihood on five items. Both notebooks name their metric and report the measured agreement rather than claiming identity. And the asymmetry that falls out is itself the lesson: difficulty b=−d/ab=-d/a is a ratio, so it is metric-free, needs no conversion, and agrees across engines at correlation 0.998. Difficulty is what a test is for and it is robustly estimated; discrimination is metric-dependent and prior-sensitive, and deserves more care when quoted.

Guessing, and a Parameter the Data Cannot Pin Down

The 3PL section produces the sharper finding, and it turns on checking a convergence flag rather than reading the estimates. The SAT12 guessing rates — ltm 0.15, mirt 0.15, from-scratch 0.19 — look close enough to call agreement. Neither maximum-likelihood fit converged: ltm's EM loop hits its iteration cap while its convergence code reports 0 (the two contradict each other), and mirt returns converged = FALSE outright.

That is not a defect of the packages so much as a fact about the 3PL. With 32 items, 600 examinees, and a guessing parameter identified almost entirely by the handful of lowest-ability examinees, the likelihood is nearly flat in cc and there is little for an optimiser to climb. The Bayesian augmentation sampler fits the same data without complaint precisely because its Beta(1,4) prior supplies the curvature the likelihood lacks — and that prior has mean exactly 0.20, which is why its posterior mean sits above the MLEs. So the 0.15-against-0.19 gap is not two estimators disagreeing on equal terms: one pair never converged, and the other converges because a prior is holding it up. These data do not pin guessing down at all, and any 3PL guessing rate quoted from either side should carry that caveat.

How well do the samplers really agree?

The PyMC cross-check deserves the same scrutiny, because the agreement statistic is itself computed from a chain that has not cleanly converged. Pooling all four chains gives a discrimination correlation of 0.957, with the worst item differing 0.51 against 0.55 — an 8% gap — at max R-hat 1.03 and a smallest bulk ESS of 349, both past the thresholds PyMC warns at. Recomputing that same correlation from subsets of the identical four chains shows how much of it is noise: any two of them would have reported somewhere between 0.941 and 0.970. That spread is the same order as the distance from the pooled figure to a perfect 1.0, so the last two digits are not a finding, and a cross-check summarised by one correlation is worth recomputing before it is quoted. Why NUTS struggles here is worth naming: 1,000 ability parameters against 5 items, a scale fixed only by the θ∼N(0,1)\theta\sim N(0,1) prior, and a sign that is not identified at all — nothing stops a chain exploring (−a,−θ)(-a,-\theta). The Albert–Chib sampler sidesteps every bit of that by drawing from exact conjugate conditionals and pinning the sign each sweep, and it is the reason difficulty — a ratio, and far better behaved — agrees at 0.999 where discrimination does not.

Notebooks

Python ICCs and the 1PL/2PL contrast, parameter recovery against simulated truth, the LSAT beside classical test theory, a 3PL with guessing on SAT12, test information and standard error, and a PyMC cross-check with its convergence reported (Python) Open → R Independent R implementation — ltm rasch/ltm/tpm with a 1PL-vs-2PL likelihood-ratio test, mirt item parameters, test information and EAP abilities, and both 3PL convergence flags checked rather than assumed (R) Open →

Downloads

irt_2pl3pl.py Albert–Chib augmentation Gibbs samplers for the normal-ogive 2PL and 3PL (the latter with a latent knows-it indicator), item characteristic curves, and Fisher item and test information with standard errors (NumPy / SciPy) lsat.csv LSAT section 6 — 1,000 examinees, 5 dichotomous items, the standard Rasch benchmark sat12.csv SAT12 — 600 examinees, 32 five-option multiple-choice items scored correct/incorrect

IRT Module — Source Code

"""
irt.py -- ITEM RESPONSE THEORY: 2PL and 3PL models (from scratch).

Backs the notebooks in  "Item Response Theory -- 2PL & 3PL".

Item response theory models a test: each person i has a continuous latent ABILITY theta_i, and
each binary item j has an item characteristic curve giving the probability of a correct answer as
a smooth function of ability. In the normal-ogive parameterisation

    2PL:  P(x_ij = 1 | theta_i) = Phi( a_j theta_i + d_j )
    3PL:  P(x_ij = 1 | theta_i) = 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 (the ability at which P = (1+c)/2), and c_j is the GUESSING floor (3PL only).
Setting all a_j equal is the RASCH / 1PL model -- a random-effects logistic regression, which is
why this connects to the panel-logit notebooks; here we read it psychometrically. (Normal-ogive
a's map to the logistic metric by a x 1.7.)

The from-scratch sampler is exactly ALBERT & CHIB's data augmentation (Albert 1992 -- the paper
that INTRODUCED this trick, for IRT, before it spread to probit/Tobit/multivariate-probit). Introduce
z_ij ~ N(a_j theta_i + d_j, 1) truncated by the sign of x_ij; then theta and (a_j, d_j) have
conjugate Gaussian full conditionals (a factor-analysis Gibbs). Scale/location are fixed by the
theta ~ N(0,1) prior, sign by keeping sum_j a_j > 0. The 3PL adds a latent "knows-it" indicator
K_ij ~ Bernoulli(Phi(a theta + d)) (Beguin & Glas 2001): correct answers by non-knowers are guesses,
so c_j gets a Beta posterior and the augmentation z is built on K rather than x.
"""

import numpy as np
from scipy.special import ndtr as Phi, ndtri as Phinv


def _rtrunc_sign(mean, positive, rng):
    """draw N(mean,1) truncated to (0,inf) where positive else (-inf,0), by inverse-CDF."""
    lo = np.where(positive, Phi(-mean), 0.0); hi = np.where(positive, 1.0, Phi(-mean))
    u = lo + rng.random(mean.shape) * (hi - lo)
    return mean + Phinv(np.clip(u, 1e-12, 1 - 1e-12))


# --------------------------------------------------------------------------- #
#  simulation                                                                   #
# --------------------------------------------------------------------------- #

def simulate_irt(n, a, d, rng, c=None, theta=None):
    """simulate a response matrix from a 2PL (c=None) or 3PL item bank."""
    a = np.asarray(a, float); d = np.asarray(d, float); J = len(a)
    if theta is None:
        theta = rng.standard_normal(n)
    P = Phi(np.outer(theta, a) + d[None, :])
    if c is not None:
        c = np.asarray(c, float); P = c[None, :] + (1 - c[None, :]) * P
    X = (rng.random((n, J)) < P).astype(float)
    return X, theta


# --------------------------------------------------------------------------- #
#  2PL Gibbs (Albert-Chib augmentation)                                         #
# --------------------------------------------------------------------------- #

def irt2pl_gibbs(X, rng, draws=3000, burn=1500, a_sd=3.0):
    """normal-ogive 2PL via Albert-Chib augmentation. Returns posterior draws of discriminations
    a:(draws,J), intercepts d:(draws,J) and abilities theta:(draws,N). Difficulty = -d/a."""
    X = np.asarray(X, float); N, J = X.shape
    theta = rng.standard_normal(N); a = np.ones(J); d = np.zeros(J)
    Pp = np.diag([1.0 / a_sd ** 2, 1e-6])              # prior precision for (a_j, d_j)
    A = np.empty((draws, J)); D = np.empty((draws, J)); TH = np.empty((draws, N))
    for it in range(draws + burn):
        mu = np.outer(theta, a) + d[None, :]
        z = _rtrunc_sign(mu, X > 0.5, rng)
        # theta_i | z, a, d  (N(0,1) prior fixes the scale)
        prec = 1.0 + np.sum(a ** 2)
        theta = ((z - d[None, :]) @ a) / prec + rng.standard_normal(N) / np.sqrt(prec)
        # (a_j, d_j) | z, theta  -- Bayesian regression of z[:,j] on [theta, 1]
        Xd = np.column_stack([theta, np.ones(N)])
        V = np.linalg.inv(Xd.T @ Xd + Pp); L = np.linalg.cholesky(V)
        M = V @ (Xd.T @ z)                             # (2, J)
        draw = M + L @ rng.standard_normal((2, J))
        a, d = draw[0], draw[1]
        if a.sum() < 0:                                # identify sign of the trait
            a = -a; theta = -theta
        if it >= burn:
            A[it - burn] = a; D[it - burn] = d; TH[it - burn] = theta
    return dict(a=A, d=D, theta=TH)


# --------------------------------------------------------------------------- #
#  3PL Gibbs (add a latent "knows-it" indicator for guessing)                   #
# --------------------------------------------------------------------------- #

def irt3pl_gibbs(X, rng, draws=4000, burn=2000, a_sd=3.0, cg=(1.0, 4.0)):
    """normal-ogive 3PL. cg = Beta prior on guessing c_j (default Beta(1,4), mean 0.2).
    Returns a, d, theta and guessing c:(draws,J)."""
    X = np.asarray(X, float); N, J = X.shape
    theta = rng.standard_normal(N); a = np.ones(J); d = np.zeros(J); c = np.full(J, 0.2)
    Pp = np.diag([1.0 / a_sd ** 2, 1e-6])
    A = np.empty((draws, J)); D = np.empty((draws, J)); TH = np.empty((draws, N)); C = np.empty((draws, J))
    for it in range(draws + burn):
        Pknow = Phi(np.outer(theta, a) + d[None, :])   # prob of KNOWING each item
        # latent knows-it indicator K: x=0 -> K=0 (would have answered right); x=1 -> Bernoulli
        pK1 = Pknow / (Pknow + (1 - Pknow) * c[None, :] + 1e-12)
        K = np.where(X > 0.5, (rng.random((N, J)) < pK1).astype(float), 0.0)
        # augment z on K (the probit "knows" process)
        mu = np.outer(theta, a) + d[None, :]
        z = _rtrunc_sign(mu, K > 0.5, rng)
        prec = 1.0 + np.sum(a ** 2)
        theta = ((z - d[None, :]) @ a) / prec + rng.standard_normal(N) / np.sqrt(prec)
        Xd = np.column_stack([theta, np.ones(N)])
        V = np.linalg.inv(Xd.T @ Xd + Pp); L = np.linalg.cholesky(V)
        M = V @ (Xd.T @ z); draw = M + L @ rng.standard_normal((2, J)); a, d = draw[0], draw[1]
        if a.sum() < 0: a = -a; theta = -theta
        # guessing c_j | data: among non-knowers (K=0), fraction answering correctly
        non = K < 0.5; corr = ((X > 0.5) & non).sum(0); tot = non.sum(0)
        c = rng.beta(cg[0] + corr, cg[1] + (tot - corr))
        if it >= burn:
            A[it - burn] = a; D[it - burn] = d; TH[it - burn] = theta; C[it - burn] = c
    return dict(a=A, d=D, theta=TH, c=C)


# --------------------------------------------------------------------------- #
#  item characteristic curve + information functions                            #
# --------------------------------------------------------------------------- #

def icc(theta, a, d, c=0.0):
    """item characteristic curve P(correct | theta)."""
    P = Phi(a * theta + d)
    return c + (1 - c) * P

def item_information(theta, a, d, c=0.0):
    """Fisher item information I_j(theta) = (P')^2 / (P(1-P)); normal-ogive metric."""
    phi = np.exp(-0.5 * (a * theta + d) ** 2) / np.sqrt(2 * np.pi)
    P = icc(theta, a, d, c); Pp = (1 - c) * a * phi
    return Pp ** 2 / np.clip(P * (1 - P), 1e-9, None)

def test_information(theta, a, d, c=None):
    """total test information (sum over items) and the standard error 1/sqrt(TIF)."""
    a = np.asarray(a); d = np.asarray(d); c = np.zeros_like(a) if c is None else np.asarray(c)
    TIF = np.zeros_like(theta, float)
    for j in range(len(a)):
        TIF += item_information(theta, a[j], d[j], c[j])
    return TIF, 1.0 / np.sqrt(np.clip(TIF, 1e-9, None))

References