Rasch IRT Model fit as Multilevel Logistic Regression: LSAT6 Data
Carolyn J. Anderson
4/19/2022
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,, 5$. The person parameter is $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
library(lme4)Read in data:
lsat <- read.table(file="C:/Users/cja/Dropbox/edps587/lectures/10 logreg/lsat6_long_data.txt", header=TRUE)
head(lsat, n=25)## count j id i1 i2 i3 i4 i5 y
## 1 3 1 1 1 0 0 0 0 0
## 2 3 1 1 0 1 0 0 0 0
## 3 3 1 1 0 0 1 0 0 0
## 4 3 1 1 0 0 0 1 0 0
## 5 3 1 1 0 0 0 0 1 0
## 6 3 2 2 1 0 0 0 0 0
## 7 3 2 2 0 1 0 0 0 0
## 8 3 2 2 0 0 1 0 0 0
## 9 3 2 2 0 0 0 1 0 0
## 10 3 2 2 0 0 0 0 1 0
## 11 3 3 3 1 0 0 0 0 0
## 12 3 3 3 0 1 0 0 0 0
## 13 3 3 3 0 0 1 0 0 0
## 14 3 3 3 0 0 0 1 0 0
## 15 3 3 3 0 0 0 0 1 0
## 16 6 1 4 1 0 0 0 0 0
## 17 6 1 4 0 1 0 0 0 0
## 18 6 1 4 0 0 1 0 0 0
## 19 6 1 4 0 0 0 1 0 0
## 20 6 1 4 0 0 0 0 1 1
## 21 6 2 5 1 0 0 0 0 0
## 22 6 2 5 0 1 0 0 0 0
## 23 6 2 5 0 0 1 0 0 0
## 24 6 2 5 0 0 0 1 0 0
## 25 6 2 5 0 0 0 0 1 1tail(lsat, n=10)## count j id i1 i2 i3 i4 i5 y
## 4991 298 297 999 1 0 0 0 0 1
## 4992 298 297 999 0 1 0 0 0 1
## 4993 298 297 999 0 0 1 0 0 1
## 4994 298 297 999 0 0 0 1 0 1
## 4995 298 297 999 0 0 0 0 1 1
## 4996 298 298 1000 1 0 0 0 0 1
## 4997 298 298 1000 0 1 0 0 0 1
## 4998 298 298 1000 0 0 1 0 0 1
## 4999 298 298 1000 0 0 0 1 0 1
## 5000 298 298 1000 0 0 0 0 1 1Guass quadrature
This is the MLE gold standard. Note that nAGQ equals the number of quadrature points, and \(i_{ij}= x_{ij}\) define above.
rasch.quad <- glmer(y ~ -1 + i1 + i2 + i3 +i4 + i5 + (1 | id),data=lsat,family=binomial, nAGQ=10 )
summary(rasch.quad)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: y ~ -1 + i1 + i2 + i3 + i4 + i5 + (1 | id)
## Data: lsat
##
## AIC BIC logLik deviance df.resid
## 4945.9 4985.0 -2466.9 4933.9 4994
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9889 0.2042 0.3436 0.5108 1.5268
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.5702 0.7551
## Number of obs: 5000, groups: id, 1000
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## i1 2.73001 0.13044 20.929 < 2e-16 ***
## i2 0.99861 0.07918 12.612 < 2e-16 ***
## i3 0.23985 0.07178 3.342 0.000833 ***
## i4 1.30645 0.08464 15.435 < 2e-16 ***
## i5 2.09940 0.10545 19.909 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## i1 i2 i3 i4
## i2 0.123
## i3 0.080 0.115
## i4 0.132 0.146 0.111
## i5 0.137 0.139 0.095 0.148Unfortunately 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)
rasch.laplace <- glmer(y ~ -1 + i1 + i2 + i3 +i4 + i5 + (1 | id),data=lsat, family=binomial)
summary(rasch.laplace)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: y ~ -1 + i1 + i2 + i3 + i4 + i5 + (1 | id)
## Data: lsat
##
## AIC BIC logLik deviance df.resid
## 4950.8 4989.9 -2469.4 4938.8 4994
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9366 0.2108 0.3468 0.5130 1.4622
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.502 0.7085
## Number of obs: 5000, groups: id, 1000
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## i1 2.70468 0.12872 21.012 < 2e-16 ***
## i2 0.99366 0.07855 12.650 < 2e-16 ***
## i3 0.23736 0.07123 3.333 0.000861 ***
## i4 1.29893 0.08390 15.482 < 2e-16 ***
## i5 2.08222 0.10417 19.989 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## i1 i2 i3 i4
## i2 0.115
## i3 0.071 0.103
## i4 0.124 0.138 0.100
## i5 0.126 0.132 0.086 0.141Alernative 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\))
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)## Generalized linear mixed model fit by maximum likelihood (Adaptive
## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
## Family: binomial ( logit )
## Formula: ycorrect ~ -1 + i1 + i2 + i3 + i4 + i5 + (1 | id)
## Data: lsat
##
## AIC BIC logLik deviance df.resid
## 4945.9 4985.0 -2466.9 4933.9 4994
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5268 -0.5108 -0.3436 -0.2042 3.9889
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.5702 0.7551
## Number of obs: 5000, groups: id, 1000
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## i1 -2.73002 0.13044 -20.929 < 2e-16 ***
## i2 -0.99861 0.07918 -12.612 < 2e-16 ***
## i3 -0.23986 0.07177 -3.342 0.000832 ***
## i4 -1.30646 0.08464 -15.436 < 2e-16 ***
## i5 -2.09940 0.10543 -19.912 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## i1 i2 i3 i4
## i2 0.123
## i3 0.080 0.115
## i4 0.132 0.146 0.111
## i5 0.137 0.139 0.095 0.148Or is you want typical paramterization of Rasch model, re-code so that \(x_{ij} = -1 or 0\).
The Item Characteristic Curves:
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')## Warning in plot.xy(xy.coords(x, y), type = type, ...): "color" is not a
## graphical parameterabline(h=.5)