---
title: "Stat 5421 Lecture Notes: Rasch Models"
author: "Charles J. Geyer"
date: "`r format(Sys.time(), '%B %d, %Y')`"
output:
  bookdown::html_document2:
    number_sections: true
    md_extensions: -tex_math_single_backslash
    mathjax: https://cdn.jsdelivr.net/npm/mathjax@3.2.2/es5/tex-mml-chtml.min.js
    css: bar.css
  bookdown::pdf_document2:
    extra_dependencies: "amscd"
    number_sections: true
    md_extensions: -tex_math_single_backslash
linkcolor: blue
urlcolor: blue
bibliography: ["rasch.bib", "expfam.bib"]
csl: journal-of-the-royal-statistical-society.csl
link-citations: true
---

# License

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.

# CRAN

Rasch models are a special case of item response models (IRT).

[CRAN](https://cran.r-project.org) has a task view
[https://cran.r-project.org/view=Psychometrics](https://cran.r-project.org/view=Psychometrics) that has several packages listed in its IRT section.

We use one of them.

# R

 * The version of R used to make this document is `r getRversion()`.

 * The version of the `bookdown` package used to make this document is
   `r packageVersion("bookdown")`.

 * The version of the `rmarkdown` package used to make this document is
   `r packageVersion("rmarkdown")`.

 * The version of the `eRm` package used to make this document is
   `r packageVersion("eRm")`.

```{r libraries}
library(eRm)
```

# Item Response Theory

Item response theory (IRT) begins with @rasch.  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](https://www.stat.umn.edu/geyer/3701/notes/mcmc-bayes.html#volleyball) 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.

# An Infinite Number of Nuisance Parameters

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-scott.

Within the theory of exponential families, there is a complete solution to the
problem.

::: {.theorem #conditional}
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](expfam.html#eq:cumfun) 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 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-no-n
(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.

# Computation

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-compute.

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*}

::: {#rasch-recursion .paragraph}
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`.

# Estimation

First we need some data.  We will get that from CRAN package `eRm`.
```{r data}
data(raschdat1)
dim(raschdat1)
```

And now we need a conditional maximum likelihood fitter.
```{r fit}
rout <- RM(raschdat1)
rout$beta
```
We 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.
```{r fit.tr}
trout <- RM(t(raschdat1))
trout$beta
```

The estimates come with standard errors.  There is not what
[IMHO](https://en.wiktionary.org/wiki/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.

# Bibliography