Introduction

Starting with the simplest non-trivial case, we assume that the test consists of \(n\) equally difficult Rasch items with \[ \pi(\theta) = \frac{\exp(\theta)}{1+\exp(\theta)}, \] the probability to find the correct answer on any of the items. The ML estimate (MLE) of \(\pi(\theta)\) is \(\hat{\pi} = \frac{x_+}{n}\) where \(x_+\) is the number of items answered correctly. From the invariance property of the MLE, it follows that \[ \hat{\theta} = \ln \frac{\frac{x_+}{n}}{1-\frac{x_+}{n}} \] is the MLE of ability.

Unfortunately, \(\hat{\theta}\) is infinite for extreme scores; that is, \(-\infty\) when \(x_+=0\) or \(+\infty\) when \(x_+=n\). One could say that \(\hat{\theta}\) is extremely biased and (under) over-estimates the ability corresponding to (\(x_+=0\)) \(x_+=n\). To reduce bias, Haldane (1956) suggested that we replace \(\hat{\pi}\) by \[ \hat{\pi}^*=\frac{x_++\frac{1}{2}}{n+1}. \] If we do so, we find \[ \hat{\theta}^* = \ln \frac{x_++\frac{1}{2}}{n+1-\left(x_+ + \frac{1}{2}\right)} \] which we will call Haldane’s estimate. As Haldene’s estimate is finite even for the extreme scores, it is obviously less biased than the MLE.

The Warm estimate was similarly developed to prevent small-sample bias but, apparently, using a very different approach. Specifically, the Warm estimator is a weighted maximum likelihood estimate (WMLE) that maximizes a weighted likelihood function \(L(\theta)W(\theta)\). We obtain the Warm estimate when \(W(\theta) = \sqrt{I(\theta)}\), where \(I(\theta)\) is the test information function. Here,
\[ L(\theta)\sqrt{I(\theta)} \propto \frac{e^{x_+\theta}}{(1+e^\theta)^{n}}\sqrt{ \frac{ne^{\theta}}{(1+e^{\theta})^2}} \propto\frac{e^{(x_++\frac{1}{2})\theta}}{(1+e^\theta)^{n+1}} \] As the right-hand side is the original likelihood function with \(x_+\) replaced by \(x_++1/2\) and \(n\) replaced by \(n+1\), it follows that the Warm estimate of ability equals Haldane’s estimate; see also Verhelst, Verstralen, and Jansen (1997).

This blog intends to explain the Warm estimate. We assume that the IRT model is a nominal response model. That is, a person with ability \(\theta\) answers independently to different items and has a probability
\[ P_{ij}(\theta) = \frac{b_{ij} e^{a_{ij}\theta}}{\sum_h b_{ih} e^{a_{ih}\theta}} \] to earn an item score of \(a_{ij}\) on item \(i\) where \(b_{ij}\) is an item category parameter with \(j=0, \dots, m_i\), \(a_{i0}=0\) and \(b_{i0}=1\). We focus on the estimation of \(\theta\) and we assume both the \(a_{ij}\) and \(b_{ij}\) to be fixed. The IRT model is then a member of the exponential family with canonical parameter \(\theta\) where the test score is a sufficient statistic for \(\theta\). This includes the Rasch model and the OPLM as special case. If the \(a_{ij}\) are fixed, it also includes the 2PL.

Preliminaries

We have a scalar parameter \(\theta\) and a statistic \(\hat{\theta}\). When estimating \(\theta\), \(\hat{\theta}\) is called an estimator and its observed value an estimate. For simplicity, we use the single term estimate and a single symbol; i.e., \(\hat{\theta}\) for the MLE and \(\hat{\theta}^*\) for the Warm estimate.

First-order bias of the MLE

Using a Taylor-series expansion to approximate \(S(\hat{\theta})\) around the true value \(\theta\), Lord (1983) found: \[ E[\hat{\theta}-\theta] = \frac{-J(\theta)}{2I^2(\theta)} +o(n^{-1}) \] where \(J(\theta)=I'(\theta)=-S''(\theta)\). Lord’s result is an instance of a general result derived, for example, in Cox and Hinkley (1974) or, slightly easier, in Mardia, Southworth, and Taylor (1999).

Note that \(I(\theta)=\sum_i I_i(\theta)\) with \(I_i(\theta)\) the item information function. If we write \(\bar{I}(\theta) = n^{-1}\sum_i I(\theta)\) for the average information per item, we find that \[ b(\theta) = -\frac{\bar{J}(\theta)}{2\bar{I}^2(\theta)} =-\frac{d}{d\theta} \ln \sqrt{\bar{I}(\theta)} \] Thus, the first-order bias is zero at the value of \(\theta\) where the (average) information function reaches its maximum value. It follows that first-order bias is smallest for those persons whose ability matches the difficulty of the test. For example, in the case of equally difficult Rasch items, \[ b(\theta) = -\frac{1-e^{\theta-\delta}}{2(1+e^{\theta-\delta})} \] which is zero when \(\theta=\delta\) and the ability of the test taker matches exactly the difficulty of the items.

Correction or prevention

Lord suggests to correct the MLE for first-order bias. That is, to use \(\tilde{\theta} = \hat{\theta} - \frac{b(\theta)}{n}\) whose bias \(B(\tilde{\theta})=o(n^{-1})\). Efron and others (1975) shows that this still holds if we substitute \(\theta=\hat{\theta}\) for the true value. The problem is that this ‘bias-corrected’ estimate would be undefined when the ML estimates are infinite.

Like Haldane, Warm suggests instead a method to prevent bias. Namely, to weight (or bias) the score function to reduce the bias in the MLE. He found that: \[ B(\theta^*) = \frac{1}{I(\theta)}\left[w'(\theta) - \frac{J(\theta)}{2I(\theta)}\right] + o(n^{-1}) \] such that \(B(\theta^*) = o(n^{-1})\) if we define \(w(\theta) =\frac{1}{2}\ln I(\theta)\).

Warm also succeeded in showing that the asymptotic distribution of the WMLE equals that of the MLE. As \(n\) grows to infinity, both estimates are normally distributed around the true value with variance \(I(\theta)^{-1}\).

The connection to Jeffrey’s prior

If \(W(\theta)\) is a prior density function of \(\theta\), the weighted likelihood is (proportional to) a posterior density function. It follows that \(\hat{\theta}^*\) is a Bayesian modal estimate of \(\theta\) when the prior \(W(\theta) \propto \sqrt{\ln I(\theta)}\). This prior is is known as Jeffreys prior.

Jeffreys prior was named after the statistician who introduced it as a prior distribution that results in a posterior that is invariant under reparametrization Jeffreys (1946). Warm notes the connection but doesn’t seem to believe that Jeffreys prior makes sense. He writes:

Suppose the same test was given to two groups of examinees, whose distributions of \(\theta\) are different and known. Bayesians will be obliged to use different priors for the two groups. For WLE the same \(w(\theta)\) would be used for both groups, because \(w(\theta)\) is a function of the test only. Similarly, if two different tests of the same ability are given to a single group, Bayesians would use the same prior for both tests, whereas for WLE a different \(w(\theta)\) would be used for each test.

Personally, we see no objection to being ignorant of group differences and base our ability estimate on the test only.

Calculating the Warm estimate in R

The Warm estimate is not easy to calculate.

Warm_Rasch = function(test_score, delta)
{
  Sw = function(theta) 
  {
    Pi = 1/(exp(delta - theta)+1)
    dPi = Pi*(1 - Pi)
    I = sum(dPi)
    J = sum(dPi*(1 - 2*Pi))
    return(test_score - sum(Pi) + J/(2*I))
  }
  return(uniroot(Sw, interval=c(-10,10))$root)
}

Warm = function(test_score, parms)
{
  a = parms$inputs$ssIS$item_score
  b = parms$est$b
  first = parms$inputs$ssI$first
  last = parms$inputs$ssI$last
  if (!(test_score %in% (0:sum(a[last])))) error("wrong test score")
  
  weighted_score = function(theta)
  {  
    nI=length(first)
    I = J = matrix(0,nI, length(theta))
    E = rep(0,length(theta))
    for (i in 1:nI)
    {
      Fij = b[first[i]:last[i]]*exp(a[first[i]:last[i]]*theta)
      Pi = Fij/sum(Fij)
      M1 = Pi*a[first[i]:last[i]]
      M2 = M1*a[first[i]:last[i]]
      M3 = M2*a[first[i]:last[i]]
      M1 = sum(M1)
      M2 = sum(M2)
      M3 = sum(M3)
      E = E + M1
      I[i,] = M2 - M1^2
      J[i,] = M3 - 3*M1*M2 + 2*M1^3
    }
    test_score - E + sum(J)/(2*sum(I))
  }
  uniroot(weighted_score, interval=c(-10,10))$root
}

Conclusion

In closing, we mention that Firth (1993), independently of Warm, developed the idea beyond the exponential family. In these very capable hands, Warm’s legacy was taken up and generalized by Kosmidis and Firth (2009) and implemented in an R package called brgl.

Appendix: First-order bias in the MLE

We derive an expression for the first-order bias term \(b(\theta)/n\). To this effect, we use a Taylor-series expansion to approximate \(S(\hat{\theta})\) around the true ability \(\theta\). \[ S(\hat{\theta}) = S(\theta)+(\hat{\theta}-\theta)S'(\theta)+\frac{1}{2}(\hat{\theta}-\theta)^2S''(\theta) + o(n^{-1/2})= 0 \] By definition, \(S(\hat{\theta})=0\). Thus, if we take expectations on both sides we find: \[ E[\hat{\theta}-\theta]S'(\theta) + \frac{1}{2}E[(\hat{\theta}-\theta)^2]S''(\theta) = 0 \] Using:

• \(E[S(\theta)] = 0\)
• We consider the situation where all derivatives of \(S\) are independent of data.
• \(E[(\hat{\theta}-\theta)^2]=Var(\hat{\theta}) + o(n^{-1}) = \frac{1}{I(\theta)} + o(n^{-1})\),

we find that \[ E[\hat{\theta}-\theta]S'(\theta) + \frac{S''(\theta)}{2I(\theta)} +o(n^{-1})= 0 \] Rearranging shows that: \[ E[\hat{\theta}-\theta] = \frac{-S''(\theta)}{2I^2(\theta)} +o(n^{-1}) \] Since \(I(\theta)=-S'(\theta)\), it follows that \(S''(\theta)=-I'(\theta)\).