Estimating abilities: MLE and EAP

Module: Item response theory
Tranche: core (a first course)
Before this: Information, precision, and short forms
Builds toward: Estimating item parameters: JML and EM · Guessing and priors · Item banks and adaptive testing · Person fit: does a respondent’s pattern fit the model? · Missing responses: omitted, not reached, and ignorable
Starts threads: EAP estimates shrink toward the mean (→ Item banks and adaptive testing · Guessing and priors · Missing responses: omitted, not reached, and ignorable · Person fit: does a respondent’s pattern fit the model?); No maximum-likelihood estimate for all-right or all-wrong patterns (→ Item banks and adaptive testing · Estimating item parameters: JML and EM · Missing responses: omitted, not reached, and ignorable · Person fit: does a respondent’s pattern fit the model?)
Picks up threads: The sum score is sufficient under the Rasch model (from The Rasch model); Under the 2PL the weighted score is sufficient (from From the 1PL to the 4PL)
IRW tables: c19prc_uk_mcbride_2021_wordsum
From EDUC 252: slides c4, c5; problem sets PS3#2, PS5#2; code ps3/1_abilities.R, ps5/eap_versus_mle.R
Status: draft
Authorship: Written largely by Claude (Anthropic), from Ben Domingue’s EDUC 252 materials and under Ben’s direction. Ben reviews each lesson before it is marked done.
Comments: open an issue on GitHub about this lesson.

What this is for

A respondent answers all ten items of a vocabulary test correctly. What is their ability? The likelihood, which has given us every estimate so far, says it is larger than any number you care to name. Abilities come to us only through an estimator, and the estimators disagree most for the respondents a test measures worst, at the top and the bottom. With the item parameters known, this lesson starts from maximum likelihood, shows where it has no answer, brings in a prior to get the expected a posteriori (EAP) estimate, and fits four estimators to the Wordsum vocabulary test, where 12% of respondents have no maximum-likelihood estimate at all.

Goals

By the end of this lesson you should be able to:

  • write down the likelihood for one respondent’s \(\theta\) with the item parameters known, find its maximum, and give its standard error;
  • explain why a perfect pattern has no maximum-likelihood estimate, and how the alternatives stay finite;
  • compute an EAP estimate on a grid of nodes, and explain why it is pulled toward the mean and by how much;
  • say what the MLE, WLE, MAP and EAP each assume, and what each assumption means for reporting scores to individuals and for comparing groups;
  • show that under the Rasch model every one of these estimates is a function of the sum score, and that under the 2PL it isn’t.

Core ideas

One respondent, one unknown

Recall from The Rasch model: mirt estimates the item parameters first, then abilities treating them as known.

This lesson is the second step. With the item parameters known, a respondent’s answers are independent Bernoulli trials whose probabilities depend on one unknown, \(\theta\): a logistic regression with the respondent’s items as the observations (Lord, 1980). For the 2PL, \[\ell(\theta) = \sum_{i} \Big[x_i \log P_i(\theta) + (1 - x_i)\log\big(1 - P_i(\theta)\big)\Big], \qquad \frac{\partial \ell}{\partial \theta} = \sum_i a_i\big(x_i - P_i(\theta)\big),\] and the Rasch model sets every \(a_i = 1\).

Recall from The likelihood and logistic regression: at the maximum the derivative of the log likelihood (the score equation) is zero, and the curvature gives the standard error. With an intercept, the score equation says the predicted number of successes equals the observed number.

Here \(\theta\) is that intercept. Under the Rasch model the maximum-likelihood estimate (MLE) solves \(\sum_i P_i(\hat\theta) = r\): the expected number correct equals the observed number. Software solves it by Newton–Raphson, stepping from a guess by the ratio of the derivative to the curvature (Baker & Kim, 2004).

Recall from Information, precision, and short forms: the negative curvature of the log likelihood is the test information, \(I(\theta) = \sum_i a_i^2 P_i(1 - P_i)\), and the conditional standard error of measurement is \(SE(\theta) = 1/\sqrt{I(\theta)}\). For one respondent’s estimate, evaluate it at their \(\hat\theta\).

So the standard error of the MLE is \(1/\sqrt{I(\hat\theta)}\), a large-sample approximation that holds up better in the middle of a test than at its ends.

Try it: the likelihood of a response pattern

Three Rasch items (\(b\) = −1, 0, 1.5). Choose which the respondent gets right; the bold curve is that pattern’s likelihood, the pale ones the other seven.

The default pattern, 1-1-0, peaks at 1.02, with a standard error of 1.37: three items locate a respondent only roughly. The other two-right patterns peak at the same place, though lower. Tick all three: the likelihood climbs toward 1 and never turns over.

Perfect patterns have no maximum

Why does the all-correct pattern never peak? Each item’s probability of a correct answer rises with \(\theta\), so their product does too. The score equation asks the expected number correct to equal the number of items, which it approaches but never reaches. The MLE is \(+\infty\), and for an all-wrong pattern \(-\infty\).

Recall from The likelihood and logistic regression, its Go deeper callout: if a predictor perfectly separates successes from failures, the logistic regression MLE doesn’t exist (Albert & Anderson, 1984). It named a respondent who answers every item correctly as the item-response version.

It is the coin with ten heads in ten flips: \(\hat p = 1\), whose log odds is infinite. mirt reports \(\pm\)Inf; optimize() on \([-10, 10]\) returns the bound, a number that looks like an estimate (the simulation below shows it). One fix that needs no prior is the weighted likelihood estimate (WLE), which is finite for every pattern under the Rasch model, the 2PL and even the 3PL and 4PL (Magis & Verhelst, 2017); it has its own subsection below.

Assumptions: one respondent answers \(I\) dichotomous items under the 2PL (the Rasch model is \(a_i = 1\)), with the item parameters known and responses locally independent given \(\theta\). Write \(P_i = \text{logit}^{-1}\big(a_i(\theta - b_i)\big)\).

No maximum. If every \(x_i = 1\), the log likelihood is \(\ell(\theta) = \sum_i \log P_i(\theta)\), and \[\frac{d\ell}{d\theta} = \sum_i a_i\,\big(1 - P_i(\theta)\big) > 0 \quad\text{for every } \theta.\] The log likelihood is strictly increasing, bounded above by 0, and approaches 0 as \(\theta \to \infty\). It has a least upper bound but no maximum. The all-wrong pattern is the mirror image. This is separation in the sense of Albert and Anderson (1984): with \(\theta\) as the only parameter, a pattern with no variation is perfectly separated.

The WLE stays finite. Warm (1989) adds a term to the score equation that removes the first-order bias of the MLE. For the Rasch model it is \[\sum_i (x_i - P_i) + \frac{J(\theta)}{2I(\theta)} = 0, \qquad J(\theta) = \sum_i P_i(1 - P_i)(1 - 2P_i),\quad I(\theta) = \sum_i P_i(1 - P_i).\] For a perfect score the first term is \(\sum_i (1 - P_i)\), which falls to 0 as \(\theta\) grows, while the second tends to \(-1/2\): as every \(P_i \to 1\), \(1 - 2P_i \to -1\) and \(J \to -I\). As \(\theta\) falls, the first term tends to \(I\), the number of items, and the second to \(+1/2\). The left-hand side is positive at the bottom of the scale and negative at the top, so it has a root: the WLE is finite. Magis and Verhelst (2017) prove this for every response pattern under the one- to four-parameter logistic models.

The EAP and MAP stay finite. With a normal prior, the log posterior is \(\ell(\theta) - (\theta - \mu)^2/(2\sigma^2)\) plus a constant. The first term is bounded above by 0 and the second falls without limit, so the log posterior has a maximum (the MAP). The posterior is a bounded function times a normal density, so it has a finite mean (the EAP).

Why it matters later. Joint maximum likelihood, which estimates items and respondents together, has to set perfect patterns aside before it starts (Estimating item parameters). An adaptive test estimates \(\theta\) after one or two items, when all-correct patterns are the rule (Item banks and adaptive testing). And a reported score has to say something about everyone at the ceiling (What does a score mean?).

Bring in a prior: the EAP estimate

The likelihood lesson set aside the Bayesian reading, in which \(\theta\) gets a prior distribution, \(g(\theta)\), and the data turn it into a posterior, \[p(\theta \mid \mathbf{x}) = \frac{L(\theta \mid \mathbf{x})\,g(\theta)}{\int L(t \mid \mathbf{x})\,g(t)\,dt}.\] For ability, the natural prior is the distribution of \(\theta\) in the respondent’s population, which mirt estimates with the items (a normal with mean 0 and an SD fixed at 1 or estimated). The expected a posteriori (EAP) estimate is the posterior mean, and its standard error is the posterior SD (Bock & Mislevy, 1982). The integral needs no closed form: evaluate likelihood × prior at nodes \(t_1, \dots, t_Q\) and average, \[\hat\theta_{\text{EAP}} = \frac{\sum_q t_q\,L(t_q)\,g(t_q)}{\sum_q L(t_q)\,g(t_q)}.\] Nothing iterates and nothing fails: a normal prior times a bounded likelihood always has a mean, even when the likelihood never peaks. The posterior’s peak is the maximum a posteriori (MAP) estimate, also called Bayes modal (Mislevy, 1986).

Try it: prior × likelihood = posterior

The ten Wordsum items analysed below, with their Rasch difficulties. Choose a sum score and the prior’s SD (the fitted value is 1.88); each curve is scaled to a peak of 1.

At a score of 9 the MLE is 2.55, the MAP 1.95 and the EAP 2.08: the prior pulls about half a logit toward 0. At 10 the likelihood never falls, yet the posterior has a peak (2.83) and a mean (3.05). With a prior SD of 1 the EAP for a score of 9 is 1.34; with 4 it is 2.66, just beyond the MLE, as the posterior is skewed toward the open end. At a score of 6 the EAP stays within about 0.1 of the MLE at any prior SD. That pull toward the prior mean is shrinkage, largest where the likelihood is flattest: at the ends of the scale, and on short tests.

Try it: shrinkage by score, and by test length

Estimates for every score on the Wordsum items. “Test length” repeats the ten items (40 items is four copies of each); the horizontal axis is the proportion correct.

The MLE at 90% correct is 2.55 whatever the length, because repeating every item multiplies both sides of the score equation by the same number. The EAP moves toward it as the test grows, from 2.08 with 10 items to 2.42 with 40, as the likelihood sharpens. The WLE, taken up below, sits nearer the MLE and gives a perfect score a finite value (3.59 with 10 items).

Shrinkage carries into anything computed from the estimates. If you’ve done Model fit and out-of-sample prediction, it is why outfit computed from mirt’s default EAP abilities averaged about 0.9 on simulated 20-item data where the model was true.

MAP, EAP and plausible values

The MAP sits close to the EAP; it is found by Newton–Raphson with the prior adding a term to the score equation, and suits software that finds modes more easily than means. If the prior describes the population, the EAP has the smallest mean squared error averaged over that population, since squared error is what a posterior mean minimizes; shrinkage is the other side of that coin.

No individual estimate carries the spread of \(\theta\) itself (EAPs vary too little, MLEs too much), so variances, correlations and regressions built on them are off. Large-scale assessments such as NAEP use plausible values (Mislevy, 1991; Mislevy, Beaton, Kaplan & Sheehan, 1992): several random draws from each respondent’s posterior, under a prior that is a regression of \(\theta\) on background variables. Each set has the spread and group differences of \(\theta\); the analysis is run on each and combined as in multiple imputation. They are a tool for population quantities, not scores to report for individuals (Wu, 2005). What does a score mean? takes them up in full.

The weighted likelihood estimate

What if we want an estimate that borrows nothing from the population but exists for every pattern? Warm (1989) weights the likelihood by a function of \(\theta\) chosen to cancel the leading term of the MLE’s bias, the bias that Lord (1983) derives. This adds \(J(\theta)/2I(\theta)\) to the score equation (the Go deeper callout writes it out). The WLE has the MLE’s large-sample variance, less bias, and a finite value for every pattern (Magis & Verhelst, 2017). It is Firth’s (1993) bias reduction, which penalizes the likelihood by the Jeffreys prior and also gives finite logistic regressions under separation (Heinze & Schemper, 2002); for the Rasch model and 2PL the WLE is the peak of the likelihood times that prior, while under the 3PL the two part (Magis & Raîche, 2012).

How does it compare? Warm’s Monte Carlo studies favoured it over the MLE and Bayes modal estimates on several criteria across a wide range of abilities, and Kim and Nicewander (1993) found the WLE, MAP and EAP all close to unbiased where the test’s items are. That is the evidence: less bias than the MLE at the same variance, and finiteness, not optimality in general. The cost is the EAP’s lower error. The WLE doesn’t shrink, so in the simulation below its RMSE is 0.77 against the EAP’s 0.61, and its SD 1.28 against a true 1. What it offers is independence from a population prior, which matters when respondents come from different populations.

Sufficiency, revisited

Recall from The Rasch model: under the Rasch model the sum score \(r\) is sufficient for \(\theta\). The likelihood factorizes into \(e^{r\theta}\), times a term free of \(\theta\), over a term free of the pattern.

Every estimator here uses the responses only through the likelihood, so under the Rasch model every one depends on them only through \(r\): the score equations contain the data only as \(r\), and in the posterior the pattern’s own factor cancels. Ten Wordsum items give eleven possible EAPs, which is how the widgets above could list every estimate.

Recall from From the 1PL to the 4PL: under the 2PL the likelihood depends on the pattern only through the weighted score \(\sum_i a_i x_i\), which is therefore sufficient; the score equation is \(\sum_i a_i (x_i - P_i) = 0\).

So respondents with the same sum score but different items right get different estimates, ordered by \(\sum_i a_i x_i\).

Recall from From the 1PL to the 4PL, its Go deeper callout: under the 3PL the weight of each correct answer depends on \(\theta\), so no simple statistic is sufficient. There an estimate depends on the whole pattern, not on a weighted sum.

Simulate

Which estimator lands closest to the truth? The code below simulates Rasch responses to ten items (\(b\) from −2 to 2) from 5,000 respondents in two groups whose means differ by 1 (overall mean 0, SD 1), estimates every \(\theta\) four ways with the difficulties known and a normal(0, 1) prior, and reports bias and root mean squared error (RMSE) by true \(\theta\), the SD of the estimates and the group difference. It is base R and runs in a second.

With the seed as written, 112 of the 5,000 respondents (2.2%) have no finite MLE, and optimize() on a perfect score returns 10, its upper bound. The MLE and WLE are close to unbiased across the middle; the EAP and MAP are pulled toward 0, by about 0.16 within one SD of the mean, 0.5 one to two SDs out, and 1 beyond. In return the EAP’s RMSE is 0.52 in the middle, against 0.78 to 0.80 for the MLE, and 0.61 overall against 0.79 (without the 112). At the ends the order reverses. The SDs are 1.26 (MLE) and 0.80 (EAP) against a true 1, and the group difference is 0.98 (MLE), 1.03 (WLE) and 0.65 (EAP). With length.out at 20 or 40, the EAP’s group difference rises to 0.80 and 0.89.

Download this code to run it in your own R session.

With real data

Our main example is c19prc_uk_mcbride_2021_wordsum: the Wordsum vocabulary test as answered by 2,058 UK adults in the sixth wave of the COVID-19 Psychological Research Consortium study (Bentall et al., 2022; the study’s design is described by McBride et al., 2021). Wordsum is the ten-item vocabulary test of the US General Social Survey, drawn from Thorndike (1942) and analysed with item response theory by Cor, Haertel, Krosnick and Malhotra (2012): for each target word, choose the closest in meaning from five words, or “don’t know”. Its job here is to be a short test on which many respondents reach the floor or the ceiling. A response is 1 for the right word and 0 otherwise, “don’t know” included. Download the code to run every chunk yourself; it needs mirt (Chalmers, 2012).

Code 
library(mirt)
# Each IRW table has a landing page (itemresponsewarehouse.org/tables/<name>/) with a
# CSV download. This link is pinned to one version of the data.
ws_url <- "https://redivis.com/api/v1/tables/datapages.item_response_warehouse_4:v7_0.c19prc_uk_mcbride_2021_wordsum/rows?format=csv"
df <- read.csv(ws_url)
table(df$wave)   # one wave (6), so one response per respondent and item

    6 
20580
Code 
# resp is 1 when the respondent chose the listed word closest in meaning to the target
# word, and 0 for any other choice, including "don't know".
wide <- tapply(df$resp, list(df$id, df$item), function(x) x[1])
wide <- as.data.frame(wide[, order(as.numeric(gsub("\\D", "", colnames(wide))))])
c(respondents = nrow(wide), items = ncol(wide), complete = sum(complete.cases(wide)))
respondents       items    complete 
       2058          10        2058
Code 
round(colMeans(wide), 2)   # proportion correct, item by item
 Wordsum_1  Wordsum_2  Wordsum_3  Wordsum_4  Wordsum_5  Wordsum_6  Wordsum_7 
      0.69       0.78       0.27       0.84       0.67       0.82       0.37 
 Wordsum_8  Wordsum_9 Wordsum_10 
      0.30       0.67       0.29
Code 
r <- rowSums(wide)          # sum score, 0 to 10
table(r)
r
  0   1   2   3   4   5   6   7   8   9  10 
 95  83 114 152 176 282 333 253 211 203 156

There is one wave, every respondent answered all ten items, and the items run from 27% to 84% correct.

Code 
# The Rasch model. mirt fixes every slope at 1 and the mean ability at 0, and
# estimates the SD of ability; its d is an easiness, so our difficulty is b = -d.
m <- mirt(wide, 1, itemtype = "Rasch", verbose = FALSE)
b <- -coef(m, simplify = TRUE)$items[, "d"]
round(sort(b), 2)
 Wordsum_4  Wordsum_6  Wordsum_2  Wordsum_1  Wordsum_9  Wordsum_5  Wordsum_7 
     -2.51      -2.32      -1.98      -1.24      -1.13      -1.10       0.88 
 Wordsum_8 Wordsum_10  Wordsum_3 
      1.37       1.39       1.55
Code 
c(sd_theta = round(sqrt(coef(m)$GroupPars[, "COV_11"]), 2))
sd_theta 
    1.88

The difficulties fall in two groups, six easy items from −2.51 to −1.10 and four hard ones from 0.88 to 1.55, with nothing between. The SD of ability, 1.88, is the prior’s SD for the EAP and MAP estimates below.

Code 
# How many respondents have no finite maximum-likelihood estimate?
ml <- fscores(m, method = "ML", full.scores.SE = TRUE)
data.frame(zero = sum(r == 0), perfect = sum(r == 10), no_mle = sum(r %in% c(0, 10)),
           share = round(mean(r %in% c(0, 10)), 3))
zero perfect no_mle share
95 156 251 0.122
Code 
table(ml[, "F1"][r %in% c(0, 10)])   # mirt reports them as -Inf and Inf

-Inf  Inf 
  95  156

The 251 respondents with a score of 0 or 10 are 12.2% of the sample, and mirt gives them \(-\infty\) and \(+\infty\). The simulation above, with an ability SD of 1, had 2.2%.

Code 
# Four estimators. Under the Rasch model each is a function of the sum score alone,
# so one row per score says everything (the last line checks it).
meth <- c("ML", "WLE", "MAP", "EAP")
est <- lapply(meth, function(k) fscores(m, method = k, full.scores.SE = TRUE))
names(est) <- meth
by_r <- function(col) sapply(est, function(e) tapply(e[, col], r, function(z) z[1]))
tab <- data.frame(score = 0:10, n = as.vector(table(r)), round(by_r("F1"), 2),
                  SE = round(by_r("SE_F1"), 2), check.names = FALSE)
print(tab, row.names = FALSE)
 score   n    ML   WLE   MAP   EAP SE.ML SE.WLE SE.MAP SE.EAP
     0  95  -Inf -4.48 -3.49 -3.65    NA   1.48   0.96   0.94
     1  83 -3.48 -3.17 -2.72 -2.83  1.11   0.95   0.82   0.86
     2 114 -2.53 -2.41 -2.10 -2.16  0.88   0.81   0.76   0.79
     3 152 -1.83 -1.79 -1.55 -1.57  0.81   0.77   0.73   0.75
     4 176 -1.20 -1.21 -1.02 -1.02  0.79   0.77   0.73   0.74
     5 282 -0.57 -0.60 -0.48 -0.47  0.80   0.80   0.74   0.75
     6 333  0.09  0.07  0.08  0.10  0.83   0.83   0.76   0.76
     7 253  0.79  0.76  0.66  0.69  0.85   0.83   0.77   0.78
     8 211  1.56  1.46  1.27  1.33  0.91   0.85   0.80   0.82
     9 203  2.55  2.26  1.95  2.07  1.13   0.96   0.86   0.91
    10 156   Inf  3.59  2.83  3.03    NA   1.49   1.03   1.05
Code 
c(largest_spread_within_a_score = round(max(sapply(est, function(e)
  max(tapply(e[, "F1"], r, function(z) diff(range(z))), na.rm = TRUE))), 4))
largest_spread_within_a_score 
                            0

Each estimator gives one value per sum score (the largest spread within a score is 0). In the middle they nearly agree: a score of 5 gives −0.57 (MLE), −0.60 (WLE), −0.48 (MAP) and −0.47 (EAP). They part at the ends. A score of 9 gives 2.55, 2.26, 1.95 and 2.07, and a score of 10 gives \(+\infty\), 3.59, 2.83 and 3.03. The MLE’s standard error is 0.79 at a score of 4 and 1.13 at 9; the posterior SDs are smaller everywhere.

Code 
# The ML standard error is 1 / sqrt(test information) at the estimate.
fin <- is.finite(ml[, "F1"])
all.equal(unname(ml[fin, "SE_F1"]), 1 / sqrt(testinfo(m, ml[fin, "F1", drop = FALSE])),
          tolerance = 1e-4)
[1] TRUE

mirt’s MLE standard errors are \(1/\sqrt{I(\hat\theta)}\), the test information evaluated at each estimate (all.equal returns TRUE).

Code 
# EAPs are pulled toward the mean: their SD is below the model's SD of ability.
# The posterior variances make up the difference.
eap <- est$EAP
sd_model <- sqrt(coef(m)$GroupPars[, "COV_11"])
round(c(sd_model = sd_model, sd_eap = sd(eap[, "F1"]), sd_wle = sd(est$WLE[, "F1"])), 2)
sd_model   sd_eap   sd_wle 
    1.88     1.69     1.96
Code 
round(c(var_model = sd_model^2, var_eap = var(eap[, "F1"]),
        mean_posterior_var = mean(eap[, "SE_F1"]^2),
        sum = var(eap[, "F1"]) + mean(eap[, "SE_F1"]^2)), 2)
         var_model            var_eap mean_posterior_var                sum 
              3.54               2.86               0.67               3.53

The EAPs have an SD of 1.69, below the model’s 1.88; the WLEs, at 1.96, are above it. The shortfall isn’t lost: the variance of the EAPs (2.86) plus the average posterior variance (0.67) is 3.53, the model’s variance (3.54) to rounding.

Code 
# The 2PL. Now the weighted score sum(a * x) carries the information about theta,
# so respondents with the same sum score can get different estimates.
m2 <- mirt(wide, 1, itemtype = "2PL", verbose = FALSE)
c(BIC_Rasch = round(extract.mirt(m, "BIC")), BIC_2PL = round(extract.mirt(m2, "BIC")))
BIC_Rasch   BIC_2PL 
    19887     19398
Code 
a <- coef(m2, simplify = TRUE)$items[, "a1"]
round(a, 2)
 Wordsum_1  Wordsum_2  Wordsum_3  Wordsum_4  Wordsum_5  Wordsum_6  Wordsum_7 
      1.62       3.75       1.21       5.20       3.64       3.94       1.09 
 Wordsum_8  Wordsum_9 Wordsum_10 
      1.21       1.59       2.20
Code 
eap2 <- fscores(m2, method = "EAP")[, "F1"]
five <- r == 5
w <- as.matrix(wide) %*% a   # weighted score
data.frame(respondents = sum(five),
  patterns = nrow(unique(wide[five, ])),
  weighted_scores = length(unique(round(w[five], 6))),
  EAPs = length(unique(round(eap2[five], 6))),
  lowest = round(min(eap2[five]), 2), highest = round(max(eap2[five]), 2),
  rank_cor = round(cor(w[five], eap2[five], method = "spearman"), 3))
respondents patterns weighted_scores EAPs lowest highest rank_cor
282 40 40 40 -1.05 -0.04 1

BIC (Schwarz, 1978) prefers the 2PL to the Rasch model (19,398 against 19,887), and the slopes run from 1.09 to 5.20. The 282 respondents with a score of 5 gave 40 different response patterns, which have 40 different weighted scores and get 40 different EAPs, from −1.05 to −0.04, in exactly the order of their weighted scores (rank correlation 1). Under the Rasch model all 282 shared one EAP, −0.47: a full logit of difference among respondents with the same number right comes from the choice of model.

Which estimate to report

The choice rests on what we’ll assume and on what the numbers are for.

  • The MLE and WLE assume nothing about the population: the estimate comes from the respondent’s own responses, which suits comparisons across populations. The price is noise: on Wordsum the MLE’s standard errors run from 0.79 to 1.13, and 12% of respondents have no MLE.
  • The EAP and MAP assume the respondent comes from the population the prior describes. If so, they are closer to the truth on average (in the simulation, an RMSE of 0.61 against 0.79). The price is shrinkage toward that population’s mean, largest at the ends and on short tests, which carries into spreads and group differences.

For a score reported to an individual from a test like this one, the EAP with its posterior SD may be sufficient for some purposes, read as pulled toward the mean of the population the prior describes. It always exists, it is the most accurate on average for the people the test was given to, and its standard error grows where the items run out. Averaged to compare groups, EAPs understate the difference: the simulation’s group difference of 1 came back as 0.65. The MLEs and WLEs came closer, but they vary too much for spreads and correlations. Where a prior is unwelcome, say because respondents come from different populations, the WLE keeps the respondents at the floor and the ceiling whom the MLE loses. Group questions can also go into the model itself, or be answered with plausible values. And if the Rasch model fits and everyone took the same items, The Rasch model has already argued for the sum score: every estimator in the table above is a relabelling of it.

Problems

  1. Three items by hand. For Rasch items with \(b = -1, 0, 1.5\), plot the likelihood of 1-1-0. At what \(\theta\) does 1-1-0 become more likely than 1-0-0, and why is that a difficulty? (PS4#1)
  2. The Rasch score equation. Show that the Rasch MLE solves \(\sum_i P_i(\theta) = r\). Then use the factorization from the Rasch lesson’s Go deeper callout to show that the EAP depends on the pattern only through \(r\) too.
  3. Change the prior. Rescore Wordsum by EAP with a prior SD of 1 instead of 1.88 (fscores(m, method = "EAP", mean = 0, cov = 1)). Which respondents move most? What happens to the SD of the EAPs, and why?
  4. Classroom averages. A school compares classroom means of students’ EAPs from a 15-item test, with one normal prior for everyone. What goes wrong, in which direction, and what would you report instead?
  5. Does EAP reduce error? Extend the simulation: 1,000 respondents at each true \(\theta\) from −3 to 3 in steps of 0.5. Where does the EAP’s RMSE beat the MLE’s, and where not? Then score normal(0, 1) respondents with a prior of mean 1. What changes? (PS5#2)
  6. Challenge (open): the floor and the ceiling. A survey analyst wants Wordsum as a predictor of, say, trust in science. What should they do with the 251 respondents at 0 or 10: drop them, use WLEs, EAPs or plausible values, or fit the regression inside the item response model? Does it depend on the question? I don’t think there is a single answer; make the case for yours.

Stuck, or curious about something this lesson doesn’t cover? Ask Claude about this lesson opens a new conversation that already knows what the lesson is about; add your question at the end. Claude shows a caution notice on any prompt that arrives through a link; that’s expected.

Going further

  • Bock, R. D., & Mislevy, R. J. (1982). Adaptive EAP estimation of ability in a microcomputer environment. Applied Psychological Measurement, 6(4), 431–444. https://doi.org/10.1177/014662168200600405. The EAP estimate on a grid of nodes, and its posterior SD as a standard error.
  • Warm, T. A. (1989). Weighted likelihood estimation of ability in item response theory. Psychometrika, 54(3), 427–450. https://doi.org/10.1007/bf02294627. The WLE, compared with the MLE and Bayes modal estimates.
  • Magis, D., & Verhelst, N. (2017). On the finiteness of the weighted likelihood estimator of ability. Psychometrika, 82(3), 637–647. https://doi.org/10.1007/s11336-016-9518-9. The WLE is finite for every pattern under the 1PL to 4PL and many polytomous models.
  • Magis, D., & Raîche, G. (2012). On the relationships between Jeffreys modal and weighted likelihood estimation of ability under logistic IRT models. Psychometrika, 77(1), 163–169. https://doi.org/10.1007/s11336-011-9233-5. Where the WLE and the Jeffreys-prior mode coincide, and where they part.
  • Firth, D. (1993). Bias reduction of maximum likelihood estimates. Biometrika, 80(1), 27–38. https://doi.org/10.1093/biomet/80.1.27. The general form of Warm’s correction.
  • Heinze, G., & Schemper, M. (2002). A solution to the problem of separation in logistic regression. Statistics in Medicine, 21(16), 2409–2419. https://doi.org/10.1002/sim.1047. Firth’s correction as a fix for separation, the logistic-regression version of perfect patterns.
  • Kim, J. K., & Nicewander, W. A. (1993). Ability estimation for conventional tests. Psychometrika, 58(4), 587–599. https://doi.org/10.1007/bf02294829. Five estimators compared for bias, standard error and reliability.
  • Lord, F. M. (1983). Unbiased estimators of ability parameters, of their variance, and of their parallel-forms reliability. Psychometrika, 48(2), 233–245. https://doi.org/10.1007/BF02294018. The bias of the MLE, and estimators without it.
  • Mislevy, R. J. (1986). Bayes modal estimation in item response models. Psychometrika, 51(2), 177–195. https://doi.org/10.1007/bf02293979. Priors in item response models.
  • Mislevy, R. J. (1991). Randomization-based inference about latent variables from complex samples. Psychometrika, 56(2), 177–196. https://doi.org/10.1007/bf02294457. Plausible values as multiple imputation, with NAEP data.
  • Mislevy, R. J., Beaton, A. E., Kaplan, B., & Sheehan, K. M. (1992). Estimating population characteristics from sparse matrix samples of item responses. Journal of Educational Measurement, 29(2), 133–161. https://doi.org/10.1111/j.1745-3984.1992.tb00371.x. Plausible values, from the concepts to NAEP.
  • Wu, M. (2005). The role of plausible values in large-scale surveys. Studies in Educational Evaluation, 31(2–3), 114–128. https://doi.org/10.1016/j.stueduc.2005.05.005. Plausible values and individual estimates, side by side.
  • Baker, F. B., & Kim, S.-H. (2004). Item response theory: Parameter estimation techniques (2nd ed.). Marcel Dekker. https://doi.org/10.1201/9781482276725. Newton–Raphson for abilities, and the rest of the machinery, step by step.
  • Lord, F. M. (1980). Applications of item response theory to practical testing problems. Erlbaum. https://doi.org/10.4324/9780203056615. The likelihood for ability with the items known.
  • Albert, A., & Anderson, J. A. (1984). On the existence of maximum likelihood estimates in logistic regression models. Biometrika, 71(1), 1–10. https://doi.org/10.1093/biomet/71.1.1. Separation, of which perfect patterns are an instance.
  • Andersen, E. B. (1977). Sufficient statistics and latent trait models. Psychometrika, 42(1), 69–81. https://doi.org/10.1007/BF02293746. Why only the Rasch model makes the sum score sufficient.
  • Cor, M. K., Haertel, E., Krosnick, J. A., & Malhotra, N. (2012). Improving ability measurement in surveys by following the principles of IRT: The Wordsum vocabulary test in the General Social Survey. Social Science Research, 41(5), 1003–1016. https://doi.org/10.1016/j.ssresearch.2012.05.007. Wordsum analysed with item response theory.
  • Thorndike, R. L. (1942). Two screening tests of verbal intelligence. Journal of Applied Psychology, 26(2), 128–135. https://doi.org/10.1037/h0060053. Where the Wordsum items come from.
  • McBride, O., Murphy, J., Shevlin, M., Gibson-Miller, J., Hartman, T. K., Hyland, P., et al. (2021). Monitoring the psychological, social, and economic impact of the COVID-19 pandemic in the population: Context, design and conduct of the longitudinal COVID-19 Psychological Research Consortium (C19PRC) Study. International Journal of Methods in Psychiatric Research, 30(1), e1861. https://doi.org/10.1002/mpr.1861. The study that gave Wordsum to these respondents.
  • Schwarz, G. (1978). Estimating the dimension of a model. The Annals of Statistics, 6(2), 461–464. https://doi.org/10.1214/aos/1176344136. The BIC.
  • Chalmers, R. P. (2012). mirt: A multidimensional item response theory package for the R environment. Journal of Statistical Software, 48(6). https://doi.org/10.18637/jss.v048.i06. fscores() computes every estimator in this lesson.
  • Johnson, M. S., & Sinharay, S. (2017). Bayesian estimation. In W. J. van der Linden (Ed.), Handbook of item response theory: Vol. 2. Statistical tools (pp. 237–258). Chapman and Hall/CRC. https://doi.org/10.1201/b19166-13. A handbook chapter on Bayesian estimation for item response models.

Data sources

Data from the Item Response Warehouse (each table name links to its IRW page); references from IRW biblio.

  • c19prc_uk_mcbride_2021_wordsum: Bentall, Richard and Shevlin, Mark and McBride, Orla and Murphy, Jamie and Hartman, Todd K and Levita, Liat and Gibson-Miller, Jilly and Mason, Liam and Bennett, Kate (2022). The COVID-19 Psychological Research Consortium Study, 2020-2021. [Data Collection]. Colchester, Essex: UK Data Service. 10.5255/UKDA-SN-855552. https://osf.io/v2zur/. Licence: CC BY 4.0.

For instructors

  • Source: EDUC 252 class 4 (slides c4, 44–48) and class 5 (slides c5, 64–66), problem sets 3 (#2), 4 (#1) and 5 (#2), and the code in ben-domingue/252: ps3/1_abilities.R, ps5/eap_versus_mle.R.
  • Timing: the core ideas and widgets take about 60 minutes; the Wordsum analysis works as a 25-minute lab. Every fit takes a few seconds in local R.
  • Changed from 252: the 252 code finds the MLE with optim() on a bounded interval, which returns the bound for a perfect pattern; the simulation here solves the score equation, counts the perfect patterns, and adds the WLE, the MAP and a group difference to PS5#2’s comparison of MLE and EAP. The real data (Wordsum) are new.
  • Held back: where the item parameters come from, and the EAP’s nodes and weights as the E step of EM (Estimating item parameters); EAP in adaptive testing (Item banks and adaptive testing); priors on item parameters (Guessing and priors); plausible values in full, and which estimate to put on a reported score scale (What does a score mean?). Out-of-sample prediction with MLE and EAP abilities (PS5#2b) is left as an extension of problem 5.
  • Solutions: held back in solutions/ability-estimation.qmd (not published).