Differential Item Functioning & Explanatory IRT

Python · PyMC · R  ·  Download DIF module

Impact Is Not Bias

The previous projects fitted items. This one interrogates them, and it is where IRT stops being descriptive and becomes the statistical machinery behind test fairness. Two questions have to be kept apart. Impact asks whether the groups differ in the trait — a legitimate difference, modelled by letting ability depend on a covariate. Differential item functioning asks whether an item behaves differently across groups at the same trait level. An item with DIF is biased: two equally aggressive people of different groups endorse it at different rates. A group can genuinely differ on the trait and some items can still be unfair, which is why the two parameters have to be estimated together.

P(xij=1θi)=Φ ⁣(ajθi+djδjgi),θiN(βgi,1)P(x_{ij}=1\mid\theta_i)=\Phi\!\big(a_j\theta_i+d_j-\delta_j g_i\big),\qquad \theta_i\sim N(\beta g_i,\,1)

Both from one augmentation

Both fall out of the same augmentation. Adding a per-item group shift δj\delta_j and a latent-regression coefficient β\beta leaves every conditional conjugate Gaussian, so ability, item parameters, DIF shifts and the group trait difference are all ordinary Gibbs draws. The DIF parameters are identified by centring them to sum to zero — DIF is inherently relative to the average item, since a shift common to every item is indistinguishable from a difference in the trait.

Verbal aggression

On the verbal-aggression data — 316 respondents answering 24 items about wanting or doing aggressive acts — the split does real work. The overall trait difference is β=+0.22\beta=+0.22, 95% CrI [−0.05, +0.50]: not distinguishable from zero, so men and women do not differ much in verbal aggressiveness overall. But five items function differentially, and their pattern is legible. The one item women over-endorse is a want-to-shout item; every item men over-endorse is a do-curse or do-scold item. The genders differ less in how aggressive they are than in how that aggression is expressed — exactly the distinction a total score cannot make.

Which Items Count as Flagged

Which items get named matters more here than anywhere else in the collection. Ranking by the largest and smallest δ\delta is not the same as selecting the set whose credible interval excludes zero, and on these data the two differ: S4WantShout and S2WantCurse have among the largest shifts but are not flagged, while S3DoScold is flagged and does not make the top three either way. In an analysis of test fairness, naming an unbiased item as biased is the error that matters most, so the flagged set is reported directly, split by sign, with the counts reconciling: 1 + 4 = 5.

itemdirectionflagged?among the largest δ|\delta|?
S2WantShoutwomen over-endorseyesyes
S4WantShoutnoyes
S2WantCursenoyes
S2DoCursemen over-endorseyesyes
S2DoScoldmen over-endorseyesyes
S3DoCursemen over-endorseyesyes
S3DoScoldmen over-endorseyesno

What the result rests on

Two things are worth adding alongside it. R's difR reaches the question by an entirely different route — a logistic DIF test with a Benjamini–Hochberg correction rather than a Bayesian credible interval — and flags four items, all four of which appear among the five. Two methods sharing no machinery agreeing on which items misbehave is worth more than either flag count alone. And the sample deserves stating plainly: 73 men against 243 women. That imbalance is why most intervals in the forest plot span zero. The study has limited power, so an absent flag is weak evidence of fairness, and the flagged set is the conservative end of what may be there.

What Makes an Item Hard

The closing section turns the question around: not is this item biased but what makes an item hard. Regressing item easiness on the items' own design features — the linear logistic test model — explains R2=0.88R^2=0.88 of the difficulty variation from three properties. Merely wanting to act is far easier to admit than doing it; cursing is easier than scolding and shouting hardest of all; blaming someone else is easier than blaming yourself. Item difficulties are not arbitrary facts to be estimated one by one but a function of how the items were built, and lme4 recovers the same ordering in R.

Notebooks

Downloads

DIF Module — Source Code

"""
dif.py -- DIFFERENTIAL ITEM FUNCTIONING and explanatory IRT (from scratch).

Backs the notebooks in  "Differential Item Functioning & Explanatory IRT".

Two questions about groups in a test:

  * IMPACT (latent regression): do the groups differ in the TRAIT itself? (Are men more
    verbally aggressive than women?) This is a legitimate group difference, modelled by letting
    the ability distribution depend on a covariate:  theta_i ~ N(beta * g_i, 1).

  * DIFFERENTIAL ITEM FUNCTIONING (DIF): does an ITEM behave differently across groups AT THE
    SAME trait level? An item with DIF is biased -- two equally-aggressive people of different
    groups endorse it at different rates. Uniform DIF is a group-specific shift in item difficulty:

        P(x_ij = 1 | theta_i) = Phi( a_j theta_i + d_j - delta_j * g_i ),

    where delta_j is item j's DIF (how much harder it is for group g=1 at the same theta).

Separating impact (beta) from bias (delta) is the whole point: a group can be genuinely more
aggressive AND some items can still be biased. From scratch it is the ALBERT-CHIB augmentation
again: draw z_ij truncated by the sign of x_ij, after which theta, the item parameters (a_j,d_j),
the DIF shifts delta_j and the regression coefficient beta all have conjugate Gaussian conditionals.
Ability scale/location are fixed by theta ~ N(beta g, 1); the DIF parameters are identified by
centring them to sum to zero (DIF is relative to the average item). Items whose delta_j credible
interval excludes zero are flagged as functioning differentially -- the psychometric fairness check.
Explanatory IRT (the LLTM idea) then regresses the item difficulties on ITEM FEATURES to explain
what makes items hard.
"""

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


def _rtrunc_sign(mean, pos, rng):
    lo = np.where(pos, Phi(-mean), 0.0); hi = np.where(pos, 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_dif(N, a, d, delta, beta, rng, g=None):
    """a,d,delta:(J,); beta scalar (group trait difference); g:(N,) group indicator 0/1."""
    a = np.asarray(a, float); d = np.asarray(d, float); delta = np.asarray(delta, float); J = len(a)
    if g is None:
        g = (rng.random(N) < 0.5).astype(float)
    theta = beta * g + rng.standard_normal(N)
    eta = np.outer(theta, a) + d[None, :] - delta[None, :] * g[:, None]
    X = (Phi(eta) > rng.random((N, J))).astype(float)
    return X, theta, g


# --------------------------------------------------------------------------- #
#  Gibbs sampler                                                                #
# --------------------------------------------------------------------------- #

def dif_gibbs(X, g, rng, draws=4000, burn=2000, a_sd=3.0, dif_sd=1.0, beta_sd=3.0):
    """2PL with uniform DIF (per-item shift by group) and a latent-regression group effect.
    g:(N,) group indicator 0/1 (the focal group = 1). Returns posterior draws of discriminations
    a, easinesses d, DIF shifts delta:(draws,J), the group trait difference beta, and abilities."""
    X = np.asarray(X, float); g = np.asarray(g, float); N, J = X.shape
    n1 = g.sum()
    theta = rng.standard_normal(N); a = np.ones(J); d = np.zeros(J); delta = np.zeros(J); beta = 0.0
    A = np.empty((draws, J)); D = np.empty((draws, J)); DE = np.empty((draws, J)); B = np.empty(draws); TH = np.empty((draws, N))
    for it in range(draws + burn):
        mu = np.outer(theta, a) + d[None, :] - delta[None, :] * g[:, None]
        z = _rtrunc_sign(mu, X > 0.5, rng)
        # theta_i | z, a, d, delta, beta   (prior N(beta g_i, 1))
        prec = 1.0 + np.sum(a ** 2)
        adjr = (z - d[None, :] + delta[None, :] * g[:, None]) @ a       # sum_j a_j (z - d + delta g)
        theta = (adjr + beta * g) / prec + rng.standard_normal(N) / np.sqrt(prec)
        # (a_j, d_j) | z, theta   (regress z + delta g on [theta, 1])
        zc = z + delta[None, :] * g[:, None]
        Xd = np.column_stack([theta, np.ones(N)]); Pp = np.diag([1.0 / a_sd ** 2, 1e-6])
        V = np.linalg.inv(Xd.T @ Xd + Pp); M = V @ (Xd.T @ zc); L = np.linalg.cholesky(V)
        dr = M + L @ rng.standard_normal((2, J)); a, d = dr[0], dr[1]
        if a.sum() < 0:
            a = -a; theta = -theta; beta = -beta
        # delta_j | z, theta, a, d   (only focal-group residuals inform it)
        r = z - np.outer(theta, a) - d[None, :]                        # residual = -delta_j g_i + eps
        r1sum = (r * g[:, None]).sum(0)                                # sum over g=1
        precd = n1 + 1.0 / dif_sd ** 2
        delta = -r1sum / precd + rng.standard_normal(J) / np.sqrt(precd)
        delta -= delta.mean()                                          # identify: DIF relative to average item
        # beta | theta, g   (regress theta on g)
        precb = g @ g + 1.0 / beta_sd ** 2
        beta = (theta @ g) / precb + rng.standard_normal() / np.sqrt(precb)
        if it >= burn:
            i = it - burn; A[i] = a; D[i] = d; DE[i] = delta; B[i] = beta; TH[i] = theta
    return dict(a=A, d=D, delta=DE, beta=B, theta=TH)

References