---
title: 'Rasch IRT Model fit as Multilevel Logistic Regression: LSAT6 Data'
author: "Carolyn J. Anderson"
date: "4/19/2022"
output: 
  html_document:
    toc: yes
     
---

# Introduction

The LSAT6 data are responses by 1000 examinees to five items on the LSAST6 exam.  This data set can be found in many packages and places.  

The following code will fit the Rasch model to data, which is 
$$P(Y_{ij} =1) =  \frac{\exp(\theta_i + b_j)}{1 + \exp(\theta_i + b_j)}),$$
where $\theta_i$ is the person parameter for item j that measures how much knowledge, ability, etc a person has-- it is a latent variable.  The parameter $b_j$ is the item parameter that indicates the "difficulty" parameter (when parameterized as "- b_j" is it an "easiness" parameter).

This model is a multilevel logistic regression model.  To see this, set indicator variable for each of the five items,
$x_{ij} = 1$ for person $i$ responds to item $j$, 0 otherwise. 

Level 1 Model:

$$ P(Y_{ij}=1) = \frac{\exp(\theta_i + b_{1i}x_{i1} +  \ldots + b_{5i}x_{i5})}{1 + \exp(\theta_i + b_{1i}x_{i1} +  \ldots + b_{5i}x_{i5})}$$

Level 2 Models

$$ \theta_i = U_{0i} $$
$$ b_{ij} = \gamma_{i} $$ for j=1,\ldots, 5$.
The person parameter is $\theta_i=U_{0i} is random intercept with no fixed value for intercept (actually it's 0) and the b's are fixed parameters.  

# Set Up

The packages we will use
```{r package, message=FALSE, warning=FALSE,include=TRUE}
library(lme4)
```

Read in data: 
```{r indata}
lsat <- read.table(file="C:/Users/cja/Dropbox/edps587/lectures/10 logreg/lsat6_long_data.txt", header=TRUE)
head(lsat, n=25)
tail(lsat, n=10)
```


# Guass quadrature 
This is the MLE gold standard. Note that nAGQ equals the number of quadrature points, and $i_{ij}= x_{ij}$ define above.

```{r gauss}
rasch.quad <- glmer(y ~ -1 + i1 + i2 + i3 +i4 + i5 + (1 | id),data=lsat,family=binomial, nAGQ=10 )
summary(rasch.quad)
```
Unfortunately quadrature it tends to have convergence problems for more than one random effect. LaPlace gives nearly the same results and is a good alternative. 

# Laplace (default)

```{r laplace}
rasch.laplace <- glmer(y ~ -1 + i1 + i2 + i3 +i4 + i5 + (1 | id),data=lsat, family=binomial)
summary(rasch.laplace)
```
# Alernative Parameterizations

In the model above P(Y=0) is modeled.  To get P(Y=1) you xN either

a    flip the signs of the estimated parameters, the "estimates or )
b    recode the $x_{ij}$s

For example, recoding so that you modle $P(Y_{ij}=1$)
```{r}
lsat$ycorrect <- ifelse (lsat$y==0, 1, 0)

rasch2.quad <- glmer(ycorrect ~ -1 + i1 + i2 + i3 +i4 + i5 + (1 | id),data=lsat,family=binomial, nAGQ=10 )
summary(rasch2.quad)
```

Or is you want typical paramterization of Rasch model, re-code so that $x_{ij} = -1 or 0$.


# The Item Characteristic Curves:

```{r icc}
x <- seq(from=0, to=8, by=0.1)
theta<- x-4
item1 <-  exp(theta +2.73001)/(1 +exp(theta +2.73001))
item2 <-  exp(theta +0.99862)/(1 +exp(theta +0.99862))
item3 <-  exp(theta +0.23986)/(1 +exp(theta +0.23986))
item4 <-  exp(theta +1.30646)/(1 +exp(theta +1.30646))
item5 <-  exp(theta +2.09940)/(1 +exp(theta +2.09940))


plot(theta,item1, type='l', 
     col='blue', ylim=c(0,1),
     ylab='Probability Correct', 
     main='Item Characteristic Curves: Rasch (lsat6)')
lines(theta,item2,col='red')
lines(theta,item3,col='green')
lines(theta,item4,col='black')
lines(theta,item5,color='cyan')
abline(h=.5)
```

