Polytomous IRT — the Graded Response Model

A Curve for Every Category

Real questionnaires rarely ask yes-or-no. They ask strongly disagree through strongly agree, and a 2PL has nowhere to put the middle. Samejima's graded response model is the ordered-probit answer: a person's latent propensity on an item is $z_{ij}=a_j\theta_i+\varepsilon$, and the observed category is whichever interval $z$ falls into, cut by ordered thresholds $\gamma_{j,1}<\dots<\gamma_{j,K-1}$. Differencing adjacent cumulative curves gives a category response curve for every option, and two categories recovers the 2PL exactly.

The same augmentation, one level up

The sampler is the same Albert–Chib trick one level up: draw the latent $z$ truncated to the observed category's interval, after which ability and discrimination are conjugate Gaussian draws and the thresholds have order-constrained conditionals. It is the identical augmentation used by the sequential and ordinal-probit models elsewhere in the collection, read psychometrically.

A neuroticism scale

Applied to a five-item neuroticism scale — 2,694 respondents, six-point items — the model ranks the indicators clearly. Angers easily and irritated easily are the sharpest, panics easily the bluntest, and the sharp item's six category curves are visibly better separated: each response marks a distinct level of the trait, while the blunt item's overlap. R's mirt and ltm reproduce the item bank and agree with each other to three decimals, and a graded-versus-partial-credit comparison picks the graded family on BIC.

Biased Low, and Only Where the Items Are Sharp

The from-scratch discriminations come out biased low on exactly the sharpest items, which is the kind of pattern that is easy to read as agreement. Fitting the same data at 3,000 draws and again at 12,000 separates the two effects. At the shorter length the three blunt items sit within 3% of mirt in the common metric while the two sharp ones sit 16% and 19% below it; lengthening the chain moves that pair to 5% and 11% and leaves the blunt three where they already were. A gap that shrinks with chain length on some items and not others is not a disagreement about the model.

Diagnosed against known truth

Diagnosing it meant simulating at the real data's exact dimensions — 2,694 respondents, 5 items, 6 categories — with known discriminations, taking mirt's own fit as the truth. Five independent simulated datasets are used, so the answer does not rest on one lucky seed. At 3,000 draws the two sharpest items come back 4% and 11% low while the three blunt ones are already within 3%; at four times the length every item lands within 2% of truth and the threshold RMSE falls by a factor of 2.5. The bias is real, it is confined to the items the real-data comparison flagged, and it disappears with chain length. So it was chain length, not the model.

The mechanism is specific and worth naming. The order-constrained update draws each threshold uniformly between $\max(z\mid x=c)$ and $\min(z\mid x=c+1)$. At $N=2694$ those order statistics sit $O(1/N)$ apart, so the cutpoints crawl, and the items whose thresholds are most spread out — the sharp ones — are the last to settle. This is the well-known slow mixing of Albert–Chib ordinal augmentation at large $N$. Running the longer chain closed the gap on the two sharp items from 16% and 19% to 5% and 11%, while the rest stayed within a few percent throughout, and the residual gap is reported with its cause rather than smoothed over. The sampler is exact in the limit and visibly biased before it — worth knowing before quoting a discrimination from a short run.

Notebooks

Downloads

GRM Module — Source Code

"""
grm.py -- POLYTOMOUS IRT: the Graded Response Model (from scratch).

Backs the notebooks in  "Polytomous IRT -- the Graded Response Model".

When test items are ORDERED categories (Likert: strongly disagree ... strongly agree), the
2PL is not enough -- we need a curve for each category. 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 + noise, and the observed category is the interval that z falls in,

    x_ij = k   iff   gamma_{j,k-1} < z_ij <= gamma_{j,k},     z_ij ~ N(a_j theta_i, 1),

with item DISCRIMINATION a_j and ordered category thresholds gamma_{j,1} < ... < gamma_{j,K-1}
(and gamma_{j,0} = -inf, gamma_{j,K} = +inf). Equivalently the cumulative curve is
P(x_ij >= k) = Phi(a_j theta_i - gamma_{j,k-1}); differencing adjacent cumulatives gives the
CATEGORY RESPONSE CURVES. This is the polytomous generalisation of the 2PL (K=2 recovers it).

The from-scratch sampler is the ORDERED-PROBIT data augmentation (Albert & Chib 1993) -- the same
truncated-normal trick as the dichotomous IRT and the sequential/ordinal-probit models: draw the
latent z_ij truncated to the observed category's interval, after which theta and a are conjugate
Gaussian and the thresholds gamma have simple order-constrained conditionals. Scale and location
are fixed by theta ~ N(0,1), the discrimination sign by keeping sum_j a_j > 0. The Rasch-family
alternative (Partial Credit / Generalised Partial Credit) uses adjacent-category rather than
cumulative logits -- fitted with the packages in the R notebook.
"""

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


def _rtruncnorm(mean, lo, hi, rng):
    """draw N(mean, 1) truncated to (lo, hi], elementwise."""
    a = Phi(lo - mean); b = Phi(hi - mean)
    u = a + rng.random(mean.shape) * (b - a)
    return mean + Phinv(np.clip(u, 1e-12, 1 - 1e-12))


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

def simulate_grm(n, a, gamma, rng, theta=None):
    """simulate ordered responses (categories 1..K). a:(J,), gamma:(J,K-1) interior thresholds."""
    a = np.asarray(a, float); gamma = np.asarray(gamma, float); J = len(a)
    if theta is None:
        theta = rng.standard_normal(n)
    z = np.outer(theta, a) + rng.standard_normal((n, J))
    X = np.ones((n, J), int)
    for j in range(J):
        for c in gamma[j]:
            X[:, j] += (z[:, j] > c).astype(int)
    return X, theta


# --------------------------------------------------------------------------- #
#  Gibbs sampler (ordered-probit augmentation)                                  #
# --------------------------------------------------------------------------- #

def grm_gibbs(X, rng, draws=3000, burn=1500, a_sd=3.0):
    """graded response model via ordered-probit augmentation. X:(N,J) integer categories 1..K.
    Returns posterior draws of discriminations a:(draws,J), interior thresholds
    gamma:(draws,J,K-1) and abilities theta:(draws,N)."""
    X = np.asarray(X, int); N, J = X.shape; K = int(X.max())
    theta = rng.standard_normal(N); a = np.ones(J)
    gam = np.tile(np.linspace(-1.0, 1.0, K - 1), (J, 1)).astype(float)   # (J, K-1) interior
    A = np.empty((draws, J)); G = np.empty((draws, J, K - 1)); TH = np.empty((draws, N))
    for it in range(draws + burn):
        # full cutpoint arrays with +-inf sentinels: (J, K+1)
        full = np.concatenate([np.full((J, 1), -1e6), gam, np.full((J, 1), 1e6)], axis=1)
        lo = full[np.arange(J)[None, :], X - 1]
        hi = full[np.arange(J)[None, :], X]
        mu = np.outer(theta, a)
        z = _rtruncnorm(mu, lo, hi, rng)
        # theta | z, a  (N(0,1) prior)
        prec = 1.0 + np.sum(a ** 2)
        theta = (z @ a) / prec + rng.standard_normal(N) / np.sqrt(prec)
        # a_j | z, theta  (regression through the origin, positive)
        sa = 1.0 / a_sd ** 2
        prec_a = sa + np.sum(theta ** 2)
        a = (theta @ z) / prec_a + rng.standard_normal(J) / np.sqrt(prec_a)
        if a.sum() < 0:
            a = -a; theta = -theta
        # thresholds gamma_{j,c} | z, X  (order-constrained uniform)
        for j in range(J):
            zj = z[:, j]; xj = X[:, j]
            for c in range(1, K):
                L = gam[j, c - 2] if c >= 2 else -1e6
                U = gam[j, c] if c <= K - 2 else 1e6
                below = zj[xj == c]; above = zj[xj == c + 1]
                if below.size: L = max(L, below.max())
                if above.size: U = min(U, above.min())
                if U > L:
                    gam[j, c - 1] = rng.uniform(L, U)
        if it >= burn:
            A[it - burn] = a; G[it - burn] = gam; TH[it - burn] = theta
    return dict(a=A, gamma=G, theta=TH, K=K)


# --------------------------------------------------------------------------- #
#  category response curves + information                                       #
# --------------------------------------------------------------------------- #

def category_probs(theta, a, gamma):
    """P(x = k | theta) for k = 1..K.  gamma: (K-1,) interior thresholds. Returns (len(theta), K)."""
    gamma = np.asarray(gamma, float); K = len(gamma) + 1
    cuts = np.concatenate([[-1e6], gamma, [1e6]])
    cum = Phi(a * theta[:, None] - cuts[None, :])        # P(z > cut) = P(x >= k+1)
    return cum[:, :-1] - cum[:, 1:]                       # differences give P(x=k)

def item_information(theta, a, gamma):
    """graded-model item information I_j(theta) = sum_k (dP_k/dtheta)^2 / P_k."""
    gamma = np.asarray(gamma, float); cuts = np.concatenate([[-1e6], gamma, [1e6]])
    z = a * theta[:, None] - cuts[None, :]
    phi = np.exp(-0.5 * z ** 2) / np.sqrt(2 * np.pi)
    dcum = a * phi                                        # d/dtheta P(z>cut)
    P = category_probs(theta, a, gamma)
    dP = dcum[:, :-1] - dcum[:, 1:]
    return np.sum(dP ** 2 / np.clip(P, 1e-9, None), axis=1)

def test_information(theta, a, gamma):
    a = np.asarray(a); TIF = np.zeros_like(theta, float)
    for j in range(len(a)):
        TIF += item_information(theta, a[j], gamma[j])
    return TIF

References

• Samejima, F. (1969). Estimation of latent ability using a response pattern of graded scores. Psychometrika Monograph Supplement 34(4, Pt. 2). — the graded response model
• Albert, J. H. & Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. JASA 88(422), 669–679. — the ordered-probit augmentation implemented here
• Cowles, M. K. (1996). Accelerating Monte Carlo Markov chain convergence for cumulative-link generalized linear models. Statistics and Computing 6(2), 101–111. — the slow threshold mixing this project's discrimination bias traces to, and the Metropolis fix for it
• Muraki, E. (1992). A generalized partial credit model. Applied Psychological Measurement 16(2), 159–176. — the adjacent-category alternative compared in the R notebook
• Chalmers, R. P. (2012). mirt: A multidimensional item response theory package for the R environment. Journal of Statistical Software 48(6). — the reference fit the from-scratch sampler is measured against