11Measuring Like Ben Wright: Rasch in R, Python, and Julia
In the reliability chapter we flagged a road past classical test theory. Here we walk it, in the company of the person who paved it. Benjamin D. Wright (1926–2015) of the University of Chicago spent a career insisting that psychology could have real measurement - the kind physicists take for granted, where the ruler does not depend on what you happen to measure with it, and the measurement does not depend on which ruler you happened to grab. His software BIGSTEPS (and its descendant WINSTEPS) put the Rasch model in the hands of a generation.
Wright liked to say that if you understood the arithmetic, you could compute a Rasch analysis on the back of an envelope. Let’s honor that: we will build a working Rasch calibration - the same computation BIGSTEPS performs - from scratch, in about thirty lines, and then show it is the same thirty lines in R, Python, and Julia. Ben would have wanted his legacy to run in whatever language you love.
11.1 Learning Objectives
State the dichotomous Rasch model and why its estimates are “objective.”
Understand joint maximum likelihood estimation (JMLE) - the BIGSTEPS engine.
Calibrate items and measure persons from scratch in R.
Recognize that the method is the mathematics, not the software - the same code runs in three languages.
11.2 The Rasch Model
The dichotomous Rasch model says the probability that person \(v\) (with ability \(\theta_v\)) answers item \(i\) (with difficulty \(\delta_i\)) correctly depends on one thing: the difference between ability and difficulty.
That is the logistic curve from the z-distribution and GLM chapters, now measuring people against items on a single shared scale (the logit). Two consequences make Wright’s eyes light up:
Sufficiency. A person’s raw score (number correct) carries all the information about their ability; an item’s raw score (number of people who passed) carries all the information about its difficulty. Nothing else matters - not which items, not which people.
Objectivity. Because ability and difficulty live on the same subtractive scale, item calibrations come out (in principle) the same whatever sample you use, and person measures come out the same whatever items you use. Wright called this sample-free item calibration and test-free person measurement (Wright and Stone 1979). It is the closest psychology gets to a physicist’s ruler.
11.3 Wright’s Data: The Knox Cube Test
Wright’s running example in Best Test Design(Wright and Stone 1979) is the Knox Cube Test (KCT): a child watches the examiner tap a sequence on four cubes, then tries to repeat it. Each of the 18 items is a tapping pattern, and the items are named by the taps - 1-4, 2-3, 1-3-2-4, and so on. Longer sequences are harder, so the KCT is very nearly a perfect ladder of difficulty - which is exactly why Wright loved it for teaching.
We use Wright’s actual data: 35 children, 18 items, the very responses from Best Test Design, stored as a plain text table (data/knox.dat). The children even keep their names - Richard, Tracie, Walter, and the rest.
## Wright's own Knox Cube Test responses (Best Test Design, 1979):## 35 children (rows) by 18 tapping items (columns), 1 = repeated correctly.kct <-read.table("data/knox.dat", header =TRUE, stringsAsFactors =FALSE)taps <-c("1-4","2-3","1-2-4","1-3-4","2-1-4","3-4-1","1-4-3-2","1-4-2-3","1-3-2-4","2-4-3-1","1-3-1-2-4","1-3-2-4-3","1-4-3-2-4","1-4-2-3-4-1","1-3-2-4-1-3","1-4-2-3-1-4","1-4-3-1-2-4","4-1-3-4-2-1-4")X <-as.matrix(kct[, paste0("i", 1:18)])rownames(X) <- kct$Namecolnames(X) <- tapsX[1:6, 1:9] # a corner: the first children, the easiest items
Everyone passes the first three items and no one passes the last. Those items carry no information about differences - and you cannot compute a finite logit for a perfect or zero score (it is \(\log\) of \(\infty\)). So the calibration engine will set them aside automatically. Watch for that in a moment: 18 items go in, but only the estimable ones come out.
11.4 The Engine: Joint Maximum Likelihood (What BIGSTEPS Does)
BIGSTEPS estimates every ability and every difficulty at once by joint maximum likelihood (Wright’s “UCON”). The recipe is a beautifully simple loop of Newton steps:
Start each person’s ability and each item’s difficulty from its raw score (a log-odds).
Nudge each ability so the person’s expected score matches their observed score.
Nudge each difficulty the same way, using the item’s score.
Re-center the item difficulties (the scale needs an origin; convention puts it at the mean item).
Repeat until nothing moves.
Every “nudge” is a one-step Newton update: (observed − expected) divided by the information (the sum of \(P(1-P)\)). Here it is, whole:
rasch_jmle <-function(X, tol =1e-6, maxit =300) {# BIGSTEPS drops persons/items with perfect or zero scores (no finite estimate)repeat { ci <-colSums(X); ri <-rowSums(X) keepI <- ci >0& ci <nrow(X); keepP <- ri >0& ri <ncol(X)if (all(keepI) &&all(keepP)) break X <- X[keepP, keepI, drop =FALSE] } N <-nrow(X); L <-ncol(X) r <-rowSums(X); s <-colSums(X) b <-log(r / (L - r)) # person abilities, from raw scores d <-log((N - s) / s); d <- d -mean(d) # item difficulties, centeredfor (it in1:maxit) { d0 <- dfor (v in1:N) { Pv <-plogis(b[v] - d); b[v] <- b[v] + (r[v] -sum(Pv)) /sum(Pv * (1- Pv)) }for (i in1:L) { Pi <-plogis(b - d[i]); d[i] <- d[i] - (s[i] -sum(Pi)) /sum(Pi * (1- Pi)) } d <- d -mean(d) # fix the origin at the mean itemif (max(abs(d - d0)) < tol) break }list(difficulty = d, ability = b, items =colnames(X), iterations = it)}fit <-rasch_jmle(X)data.frame(item = fit$items, difficulty =round(fit$difficulty, 2))
item
difficulty
1-3-4
1-3-4
-4.55
2-1-4
2-1-4
-3.97
3-4-1
3-4-1
-3.51
1-4-3-2
1-4-3-2
-3.97
1-4-2-3
1-4-2-3
-2.44
1-3-2-4
1-3-2-4
-3.51
2-4-3-1
2-4-3-1
-1.63
1-3-1-2-4
1-3-1-2-4
0.83
1-3-2-4-3
1-3-2-4-3
2.33
1-4-3-2-4
1-4-3-2-4
2.03
1-4-2-3-4-1
1-4-2-3-4-1
3.50
1-3-2-4-1-3
1-3-2-4-1-3
4.96
1-4-2-3-1-4
1-4-2-3-1-4
4.96
1-4-3-1-2-4
1-4-3-1-2-4
4.96
Read that difficulty column top to bottom: the short tapping sequences sit low (easy), the long ones sit high (hard), climbing steadily just as Wright’s KCT does. We just calibrated a test the way BIGSTEPS calibrates a test - in thirty lines, and every number rebuildable by hand. The persons are measured on the very same logit scale, so a child’s ability and an item’s difficulty are directly comparable: if \(\theta_v > \delta_i\), the odds favor a correct tap.
summary(fit$ability) # the 35 children, measured in logits on the item scale
Min. 1st Qu. Median Mean 3rd Qu. Max.
-4.4750 -1.4619 -0.2759 -0.1965 0.9886 3.8891
11.5 The Same Thirty Lines in Python
Nothing above was special to R. Here is the identical algorithm in Python with NumPy - same updates, same convergence, same estimates:
import numpy as np# Wright's Knox Cube Test data: skip the header, drop the Name & sex columns.rows = [ln.split() for ln inopen("data/knox.dat")][1:]X = np.array([[int(v) for v in r[2:]] for r in rows], float)def rasch_jmle(X, tol=1e-6, maxit=300): X = np.asarray(X, float)whileTrue: # drop perfect/zero rows & cols ci, ri = X.sum(0), X.sum(1) keepI = (ci >0) & (ci < X.shape[0]) keepP = (ri >0) & (ri < X.shape[1])if keepI.all() and keepP.all():break X = X[keepP][:, keepI] N, L = X.shape r, s = X.sum(1), X.sum(0) b = np.log(r / (L - r)) # abilities from raw scores d = np.log((N - s) / s); d -= d.mean() # difficulties, centeredfor _ inrange(maxit): d0 = d.copy()for v inrange(N): Pv =1/ (1+ np.exp(-(b[v] - d))) b[v] += (r[v] - Pv.sum()) / (Pv * (1- Pv)).sum()for i inrange(L): Pi =1/ (1+ np.exp(-(b - d[i]))) d[i] -= (s[i] - Pi.sum()) / (Pi * (1- Pi)).sum() d -= d.mean() # fix the originif np.max(np.abs(d - d0)) < tol:breakreturn d, bdifficulty, ability = rasch_jmle(X)
11.6 The Same Thirty Lines in Julia
And in Julia, where the math reads almost like the equations themselves:
usingStatisticsσ(z) =1/ (1+exp(-z))# Wright's Knox Cube Test data: skip the header, drop the Name & sex columns.rows =split.(readlines("data/knox.dat")[2:end])X = [parse(Int, rows[i][j]) for i in1:length(rows), j in3:20]functionrasch_jmle(X; tol=1e-6, maxit=300) X =float.(X)whiletrue# drop perfect/zero rows & cols ci =vec(sum(X, dims=1)); ri =vec(sum(X, dims=2)) keepI = (ci .>0) .& (ci .<size(X, 1)) keepP = (ri .>0) .& (ri .<size(X, 2)) (all(keepI) &&all(keepP)) &&break X = X[keepP, keepI]end N, L =size(X) r =vec(sum(X, dims=2)); s =vec(sum(X, dims=1)) b =log.(r ./ (L .- r)) # abilities from raw scores d =log.((N .- s) ./ s); d .-=mean(d) # difficulties, centeredfor _ in1:maxit d0 =copy(d)for v in1:N Pv =σ.(b[v] .- d) b[v] += (r[v] -sum(Pv)) /sum(Pv .* (1.- Pv))endfor i in1:L Pi =σ.(b .- d[i]) d[i] -= (s[i] -sum(Pi)) /sum(Pi .* (1.- Pi))end d .-=mean(d) # fix the originmaximum(abs.(d .- d0)) < tol &&breakendreturn d, benddifficulty, ability =rasch_jmle(X)
Three languages, one algorithm, identical numbers. That is the whole point Wright kept making: the measurement is in the mathematics, not the software. BIGSTEPS was never magic - it was this loop, run carefully.
ImportantWhat the Real Tools Add
Our thirty lines calibrate items and measure persons - the heart of BIGSTEPS. Production tools add three things worth knowing about: (1) fit statistics (Wright’s infit and outfit mean-squares) that flag items and people who do not behave as the model predicts; (2) better small-sample bias corrections than raw JMLE; and (3) standard errors for every estimate. When you do this for real, reach for a maintained package - eRm or TAM in R, py-irt in Python, or a Julia IRT package - all of which are, underneath, doing what you just watched.
NoteFrom the Field
This is not a museum piece. The authors have run this exact machinery on real clinical measures - including, fittingly, a paper co-written with Ben Wright himself: Conrad, Wright, McKnight, and colleagues (Conrad et al. 2004) used Rasch analysis on the Mississippi PTSD Scale and found that five of its seven worst-fitting items were reverse-scored - the model exposing a flaw a sum score would have hidden. In another study, McKnight and Babcock-Parziale (McKnight and Babcock-Parziale 2007) paired Rasch with generalizability theory to evaluate a vision-rehabilitation outcome. Rasch is not exotic; it is a working tool for anyone who needs a measure to behave.
11.7 Challenge
TipDo One Yourself
Wright’s data already hands you the lesson: items 1-4, 2-3, 1-2-4 were passed by all 35 children and 4-1-3-4-2-1-4 by none. Confirm from colSums(X) that these are exactly the items the estimator dropped, and explain in one sentence why no finite difficulty exists for them.
Plot person ability against raw score. What shape is it, and why does every child with the same score (there are several with 10) get the same measure?
Translate one of the two for loops (ability or difficulty updates) into the language you know best, and check it against the R output. You will have written a piece of BIGSTEPS.
11.8 Where We Go Next
Wright’s dream - objective, sample-free measurement - is the aspiration behind every scale you will ever build or trust. With reliable, valid, well-measured variables in hand, we are almost ready to model relationships among them. One preparation step remains: combining our many measures into fewer, more dependable ones. That is data reduction, the subject of the next part, and it is where we turn before the modeling begins. Ben Wright would tell you: get the measurement right first, and the models take care of themselves.
Conrad, Kendon J., Benjamin D. Wright, Patrick McKnight, Miles McFall, Alan Fontana, and Robert Rosenheck. 2004. “Comparing Traditional and Rasch Analyses of the Mississippi PTSD Scale: Revealing Limitations of Reverse-Scored Items.”Journal of Applied Measurement 5 (1): 15–30.
McKnight, Patrick E., and Judith Babcock-Parziale. 2007. “Respondent Impact on Functional Ability Outcome Measures in Vision Rehabilitation.”Optometry and Vision Science 84 (8): 721–28.
Wright, Benjamin D., and Mark H. Stone. 1979. Best Test Design. MESA Press.