This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0/).
The R markdown source for this document is https://www.stat.umn.edu/geyer/5421/notes/rasch.Rmd.
Rasch models are a special case of item response models (IRT).
CRAN has a task view https://cran.r-project.org/view=Psychometrics that has several packages listed in its IRT section.
We use one of them.
The version of R used to make this document is 4.6.1.
The version of the bookdown package used to make this document is
0.47.
The version of the rmarkdown package used to make this document is
2.31.
The version of the eRm package used to make this document is
1.0.10.
library(eRm)Item response theory (IRT) begins with Rasch (1961). It is about test data. We have a test with various test items of varying difficulty and various test takers with various abilities.
There is an enormous literature on IRT and many statistical methods of analysis have been proposed for such data. We will just look at the simplest model for such data where every test item is marked right or wrong, so we have a matrix of zero-or-one-valued (Bernoulli) random variables. The rows are test takers (subjects) and the columns are test items. The Rasch model is just like the Bradley-Terry model used for the frequentist and Bayesian analysis of volleyball data except instead of the rows and columns denoting the same things (teams) for Bradley-Terry, they denote different things (subjects and items) for Rasch. Thus Rasch used the log-linear model with canonical parameter of the form \[ \theta_{a b} = \alpha_a + \beta_b \] where \(\alpha\) is the subject ability parameter vector and \(\beta\) is the item difficulty parameter vector.
There is a problem with such data. If we hold one dimension (subjects) fixed and let the other dimension (items) go to infinity, which is what we do to derive large-sample-size, asymptotic, approximations to sampling distributions, we do not even get consistent estimation. This phenomenon was first pointed out by Neyman and Scott (1948).
Within the theory of exponential families, there is a complete solution to the problem.
Theorem 5.1 In an exponential family, conditioning on some components of the canonical statistic vector results in another exponential family model whose canonical statistic is the other components of the original canonical statistic vector (those not conditioned on) and whose canonical parameters are the corresponding components of the original canonical parameter vector.
Proof. The conditional log likelihood has the form \[ l(\theta) = \left( \sum_{i \in A} y_i \theta_i \right) + \left( \sum_{i \in I \setminus A} y_i \theta_i \right) - c(\theta) \] where \(I\) is the index set of the canonical statistic and parameter vectors and \(A \subset I\). When we condition on \(y_i\), \(i \in I \setminus A\), those components of the canonical statistic vector are treated as known constants (fixed at their observed values), so this gives us a new exponential family with
canonical statistic vector \(y_A\),
canonical parameter vector \(\theta_A\)
and cumulant function whatever the conditioning process results in. But we also have the formula for cumulant functions in the exponential family notes that shows that the cumulant function is a function of the canonical parameter vector only. Hence, whatever the functional form of the cumulant function may be, we can write the log likelihood as \[ l(\theta_A) = \langle y_A, \theta_A \rangle - c_A(\theta_A) \] and we see the assertion of the theorem is correct: \(\theta_{I \setminus A}\) is gone.
The other shoe dropping is that Andersen (1970) proved consistency and asymptotic normality of these conditional maximum likelihood estimates. It needed a separate proof because the distribution is changing as the number of items goes to infinity (which is the data the conditional distribution is conditioning on). We could also use the theory presented in Geyer (2013) (based on the theory of Le Cam cited therein) that says as long as the conditional log-likelihood is approximately quadratic, then the conditional MLE are approximately normal.
So everything is great! Except we may or may not know how to calculate the conditional exponential family.
The conditional probability distribution of the data conditioning on column sums (item canonical statistic) has independent columns (because of the conditioning). So we can just derive the distribution for one column.
So, for this section only, let \(y\) denote a single column of the matrix, and let \(s\) denote the sum. We have the unconditional distribution proportional to \(\pi_i^{y_i} (1 - \pi_i)^{1 - y_i}\), where the \(\pi_i\) are between zero and one (probabilities) and the \(y_i\) are Bernoulli (zero or one). And the conditional distribution just normalizes \[ f(y | s) = \frac{\prod_i \pi_i^{y_i} (1 - \pi_i)^{1 - y_i}} {\sum_{\substack{y \in \{0, 1\}^k \\ \sum_{i = 1}^k y_i = s}} \prod_i \pi_i^{y_i} (1 - \pi_i)^{1 - y_i}} \] The tricky part is that the sum in the denominator is not over all possible values of the vector \(y\) but rather only over those values whose components sum to \(s\).
There is no easy way to do this. We are going to use R package eRm
and the vignette for this package says it uses the algorithms of
Andersen (1972).
A glance at that paper says it does the denominator by recursion. Denote the denominator by \(D(\pi, k, s)\). Then \[\begin{alignat*}{2} D(\pi, k, s) & = \pi_i D(\pi, k - 1, s - 1) + (1 - \pi_i) D(\pi, k - 1, s), \qquad 0 < s < k \\ D(\pi, k, 0) & = \prod_{i = 1}^k (1 - \pi_i) \\ D(\pi, k, k) & = \prod_{i = 1}^k \pi_i \end{alignat*}\]
If this is calculated in the most obvious way, it still requires an exponential in \(k\) number of evaluations of the function \(D\). But if we cache evaluations, there are only \(s + 1\) values \(D(\pi, k, t)\) because \(0 \le t \le s\), hence only \(k (s + 1)\) values in all. So this computation can be made very efficient if done in language that supports this algorithm.
Unfortunately, this sort of recursive computation is not very efficiently
done in R. So we won’t bother to implement it ourselves. We will use the
implementation in CRAN package eRm.
First we need some data. We will get that from CRAN package eRm.
data(raschdat1)
dim(raschdat1)## [1] 100 30And now we need a conditional maximum likelihood fitter.
rout <- RM(raschdat1)
rout$beta## beta I1 beta I2 beta I3 beta I4 beta I5 beta I6
## 1.565269744 0.051171710 0.782190083 -0.650231923 -1.300578949 0.099296295
## beta I7 beta I8 beta I9 beta I10 beta I11 beta I12
## 0.681696811 0.731734141 0.533662281 -1.107727187 -0.650231922 0.387903887
## beta I13 beta I14 beta I15 beta I16 beta I17 beta I18
## -1.511191792 -2.116116955 0.339649387 -0.597111088 0.339649389 -0.093927362
## beta I19 beta I20 beta I21 beta I22 beta I23 beta I24
## -0.758721103 0.681696809 0.936549352 0.989173483 0.681696813 0.002949604
## beta I25 beta I26 beta I27 beta I28 beta I29 beta I30
## -0.814227470 1.207133494 -0.093927363 -0.290443211 -0.758721107 0.731734147We see that contrary to the way we have been talking so far, this package is estimating the item parameters while conditioning away the subject parameters. Clearly, this does not matter, because just transposing the data matrix makes items into subjects and vice versa.
trout <- RM(t(raschdat1))
trout$beta## beta I1 beta I2 beta I3 beta I4 beta I5 beta I6
## -3.44818876 -2.70255239 -2.24317188 -1.90007430 -1.90007430 -1.90007430
## beta I7 beta I8 beta I9 beta I10 beta I11 beta I12
## -1.61968371 -1.37838597 -1.37838597 -1.37838597 -1.37838597 -1.37838595
## beta I13 beta I14 beta I15 beta I16 beta I17 beta I18
## -1.16314205 -1.16314205 -1.16314205 -1.16314203 -1.16314203 -0.96625830
## beta I19 beta I20 beta I21 beta I22 beta I23 beta I24
## -0.96625830 -0.96625830 -0.78270521 -0.78270521 -0.78270521 -0.78270521
## beta I25 beta I26 beta I27 beta I28 beta I29 beta I30
## -0.78270521 -0.60893685 -0.60893685 -0.60893682 -0.60893685 -0.60893691
## beta I31 beta I32 beta I33 beta I34 beta I35 beta I36
## -0.60893687 -0.44231944 -0.44231952 -0.44231952 -0.44231952 -0.28080450
## beta I37 beta I38 beta I39 beta I40 beta I41 beta I42
## -0.28080459 -0.28080459 -0.28080450 -0.28080450 -0.28080447 -0.12272214
## beta I43 beta I44 beta I45 beta I46 beta I47 beta I48
## -0.12272211 -0.12272212 -0.12272212 -0.12272212 -0.12272212 0.03335708
## beta I49 beta I50 beta I51 beta I52 beta I53 beta I54
## 0.03335707 0.03335708 0.03335710 0.03335710 0.03335706 0.03335705
## beta I55 beta I56 beta I57 beta I58 beta I59 beta I60
## 0.18871973 0.18871971 0.18871971 0.18871972 0.18871969 0.18871970
## beta I61 beta I62 beta I63 beta I64 beta I65 beta I66
## 0.34458700 0.34458696 0.34458699 0.34458698 0.34458699 0.50218149
## beta I67 beta I68 beta I69 beta I70 beta I71 beta I72
## 0.50218149 0.50218149 0.50218149 0.50218149 0.50218149 0.50218149
## beta I73 beta I74 beta I75 beta I76 beta I77 beta I78
## 0.50218149 0.66279491 0.66279491 0.66279489 0.82786641 0.82786641
## beta I79 beta I80 beta I81 beta I82 beta I83 beta I84
## 0.99908025 0.99908026 0.99908026 0.99908026 0.99908026 1.17850494
## beta I85 beta I86 beta I87 beta I88 beta I89 beta I90
## 1.17850496 1.17850495 1.17850491 1.17850497 1.17850497 1.36879500
## beta I91 beta I92 beta I93 beta I94 beta I95 beta I96
## 1.36879503 1.36879500 1.36879500 1.57351632 1.79773908 2.04928502
## beta I97 beta I98 beta I99 beta I100
## 2.34133649 2.34133651 2.34133651 2.34133653The estimates come with standard errors. There is not what IMHO is a complete suite of functions. There is no way, for example, to get conditional mean value parameters, even though they must have been produced as part of the likelihood equations that are solved to get CMLE.
It does do many extensions of the Rasch model. That is what the “e” in the package name is for. And the “rm” in the package name is for Rasch model. And capitalizing the R is for R.