Rasch model

Dear friends,
I am running a simple model where 424 pupils took 13 test items. I model the item easiness as fixed effects and the pupil abilities as random effects (this is basically the 'so-called' marginal maximum likelihood estimation used in some Rasch software packages. I also run the analysis using two 'traditional' Rasch packages. I have found that there is a near-perfect correlation between item estimates from the 'traditional' packages and lme4. See the results below. However, when I use ranef(model$id) to  get the 'ability' estimates of the pupils, the correlation is just 0.9! Shouldnt the correlation be much bigger? I mean, how would you estimate the ability estimates of the pupils in this context?

Thanks for any help

Generalized linear mixed model fit by the Laplace approximation
Formula: score ~ 0 + item + (1 | id)
  Data: rasch_data
 AIC  BIC logLik deviance
 6502 6595  -3237     6474
Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 1.9821   1.4079
Number of obs: 5512, groups: id, 424

Fixed effects:
        Estimate Std. Error z value Pr(>|z|)
item   1  0.05819    0.13004   0.447 0.654524
item   2  1.26791    0.14126   8.976  < 2e-16 ***
item   3  0.93972    0.13615   6.902 5.12e-12 ***
item   4  0.70219    0.13345   5.262 1.43e-07 ***
item   5  0.36877    0.13099   2.815 0.004874 **
item   6  0.52699    0.13197   3.993 6.52e-05 ***
item   7 -0.79568    0.13404  -5.936 2.92e-09 ***
item   8  0.38186    0.13106   2.914 0.003572 **
item   9 -0.61867    0.13239  -4.673 2.97e-06 ***
item  10  0.30358    0.13069   2.323 0.020184 *
item  11  0.23870    0.13044   1.830 0.067262 .
item  12 -0.79568    0.13404  -5.936 2.92e-09 ***
item  13 -0.47278    0.13137  -3.599 0.000320 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
        item1 item2 item3 item4 item5 item6 item7 item8 item9 item10 item11 item12
item   2 0.258
item   3 0.269 0.252
item   4 0.275 0.256 0.265
item   5 0.281 0.258 0.268 0.274
item   6 0.278 0.257 0.267 0.273 0.277
item   7 0.273 0.243 0.255 0.262 0.269 0.266
item   8 0.280 0.258 0.268 0.274 0.279 0.277 0.269
item   9 0.277 0.248 0.259 0.266 0.273 0.270 0.271 0.273
item  10 0.281 0.258 0.269 0.275 0.280 0.278 0.270 0.280 0.274
item  11 0.282 0.258 0.269 0.275 0.280 0.278 0.271 0.280 0.275 0.281
item  12 0.273 0.243 0.255 0.262 0.269 0.266 0.268 0.269 0.271 0.270  0.271
item  13 0.279 0.251 0.263 0.269 0.276 0.273 0.272 0.276 0.275 0.277  0.278  0.272

MML_Rasch <- lmer(score ~ 0+item+(1|id), rasch_data, family=binomial(link="logit"))
summary(MML_Rasch)

Generalized linear mixed model fit by the Laplace approximation
Formula: score ~ 0 + item + (1 | id)
  Data: rasch_data
 AIC  BIC logLik deviance
 6502 6595  -3237     6474
Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 1.9821   1.4079
Number of obs: 5512, groups: id, 424

Fixed effects:
        Estimate Std. Error z value Pr(>|z|)
item   1  0.05819    0.13004   0.447 0.654524
item   2  1.26791    0.14126   8.976  < 2e-16 ***
item   3  0.93972    0.13615   6.902 5.12e-12 ***
item   4  0.70219    0.13345   5.262 1.43e-07 ***
item   5  0.36877    0.13099   2.815 0.004874 **
item   6  0.52699    0.13197   3.993 6.52e-05 ***
item   7 -0.79568    0.13404  -5.936 2.92e-09 ***
item   8  0.38186    0.13106   2.914 0.003572 **
item   9 -0.61867    0.13239  -4.673 2.97e-06 ***
item  10  0.30358    0.13069   2.323 0.020184 *
item  11  0.23870    0.13044   1.830 0.067262 .
item  12 -0.79568    0.13404  -5.936 2.92e-09 ***
item  13 -0.47278    0.13137  -3.599 0.000320 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
        item1 item2 item3 item4 item5 item6 item7 item8 item9 item10 item11 item12
item   2 0.258
item   3 0.269 0.252
item   4 0.275 0.256 0.265
item   5 0.281 0.258 0.268 0.274
item   6 0.278 0.257 0.267 0.273 0.277
item   7 0.273 0.243 0.255 0.262 0.269 0.266
item   8 0.280 0.258 0.268 0.274 0.279 0.277 0.269
item   9 0.277 0.248 0.259 0.266 0.273 0.270 0.271 0.273
item  10 0.281 0.258 0.269 0.275 0.280 0.278 0.270 0.280 0.274
item  11 0.282 0.258 0.269 0.275 0.280 0.278 0.271 0.280 0.275 0.281
item  12 0.273 0.243 0.255 0.262 0.269 0.266 0.268 0.269 0.271 0.270  0.271
item  13 0.279 0.251 0.263 0.269 0.276 0.273 0.272 0.276 0.275 0.277  0.278  0.272


Dr. Iasonas Lamprianou
Department of Education
The University of Manchester
Oxford Road, Manchester M13 9PL, UK
Tel. 0044  161 275 3485
iasonas.lamprianou at manchester.ac.uk


--- On Sat, 21/2/09, r-sig-mixed-models-request at r-project.org <r-sig-mixed-models-request at r-project.org> wrote:

From: r-sig-mixed-models-request at r-project.org <r-sig-mixed-models-request at r-project.org>
Subject: R-sig-mixed-models Digest, Vol 26, Issue 33
To: r-sig-mixed-models at r-project.org
Date: Saturday, 21 February, 2009, 11:00 AM
Send R-sig-mixed-models mailing list submissions to
      r-sig-mixed-models at r-project.org

To subscribe or unsubscribe via the World Wide Web, visit

From: Doran, Harold <HDoran at air.org>
Subject: RE: [R-sig-ME] Rasch model
To: "Douglas Bates" <bates at stat.wisc.edu>, lamprianou at yahoo.com
Cc: r-sig-mixed-models at r-project.org
Date: Sunday, 22 February, 2009, 2:33 PM
What "traditional" method (and what packages) did
you use for ability estimation in the other packages? When
you treat students as random as you did and then get the
ability estimates, this is equivalent to MAP estimation.
Most traditional packages use ML estimation, which limits
one's ability to get estimates for those individuals
with perfect or 0 scores since the likelihood function is
unbounded. Arbitrary scoring rules are typically applied in
these cases.

You might check out the irt.ability() function in the
MiscPsycho package for generating abilities as it offer you
MAP, EAP and ML estimates in R.


-----Original Message-----
From: r-sig-mixed-models-bounces at r-project.org on behalf of
Douglas Bates
Sent: Sat 2/21/2009 9:19 AM
To: lamprianou at yahoo.com
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] Rasch model
 
On Sat, Feb 21, 2009 at 5:31 AM, Iasonas Lamprianou
<lamprianou at yahoo.com> wrote:

Dear friends,
I am running a simple model where 424 pupils took 13

test items. I model the item easiness as fixed effects and
the pupil abilities as random effects (this is basically the
'so-called' marginal maximum likelihood estimation
used in some Rasch software packages. I also run the
analysis using two 'traditional' Rasch packages. I
have found that there is a near-perfect correlation between
item estimates from the 'traditional' packages and
lme4. See the results below. However, when I use
ranef(model$id) to  get the 'ability' estimates of
the pupils, the correlation is just 0.9! Shouldnt the
correlation be much bigger? I mean, how would you estimate
the ability estimates of the pupils in this context?

Thanks for any help

As far as I know the conditional modes of the random
effects for
students should correspond to the student ability estimates
from
marginal maximum likelihood.  However, I don't know
much about
commercial Rasch packages calculate these parameters.

With more items is becomes reasonable to model both the
item easiness
and the pupil abilities as random effects which is not
easily done in
commercial packages.  There is a chapter in the book
"Exploratory Item
Response Models" describing the model and why it is
desirable but
without indication of how it could be fit.  In the paper

@article{Dowling:Bliese:Bates:Doran:2007:JSSOBK:v20i02,
  author =	"Harold Doran and Douglas Bates and Paul
Bliese and Maritza
 Dowling",
  title =	"Estimating the Multilevel Rasch Model: With
the lme4 Package",
  journal =	"Journal of Statistical Software",
  volume =	"20",
  number =	"2",
  pages =	"1--18",
  day =  	"22",
  month =	"2",
  year = 	"2007",
  CODEN =	"JSSOBK",
  ISSN = 	"1548-7660",
  bibdate =	"2007-02-22",
  URL =  	"http://www.jstatsoft.org/v20/i02",
  accepted =	"2007-02-22",
  acknowledgement = "",
  keywords =	"",
  submitted =	"2006-10-01",
}

we show how to fit that model and generalizations of it
where items
and/or subjects fall into groups.  Also, there are several
slides in
the section "Item Response Models as GLMMs" of
the workshop
presentation at
http://www.stat.wisc.edu/~bates/UseR2008/WorkshopD.pdf
related to fitting Rasch models with crossed random
effects.

Generalized linear mixed model fit by the Laplace

approximation

Formula: score ~ 0 + item + (1 | id)
  Data: rasch_data
 AIC  BIC logLik deviance
 6502 6595  -3237     6474
Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 1.9821   1.4079
Number of obs: 5512, groups: id, 424

Fixed effects:
        Estimate Std. Error z value Pr(>|z|)
item   1  0.05819    0.13004   0.447 0.654524
item   2  1.26791    0.14126   8.976  < 2e-16 ***
item   3  0.93972    0.13615   6.902 5.12e-12 ***
item   4  0.70219    0.13345   5.262 1.43e-07 ***
item   5  0.36877    0.13099   2.815 0.004874 **
item   6  0.52699    0.13197   3.993 6.52e-05 ***
item   7 -0.79568    0.13404  -5.936 2.92e-09 ***
item   8  0.38186    0.13106   2.914 0.003572 **
item   9 -0.61867    0.13239  -4.673 2.97e-06 ***
item  10  0.30358    0.13069   2.323 0.020184 *
item  11  0.23870    0.13044   1.830 0.067262 .
item  12 -0.79568    0.13404  -5.936 2.92e-09 ***
item  13 -0.47278    0.13137  -3.599 0.000320 ***
---
Signif. codes:  0 '***' 0.001 '**'

0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
        item1 item2 item3 item4 item5 item6 item7

item8 item9 item10 item11 item12

item   2 0.258
item   3 0.269 0.252
item   4 0.275 0.256 0.265
item   5 0.281 0.258 0.268 0.274
item   6 0.278 0.257 0.267 0.273 0.277
item   7 0.273 0.243 0.255 0.262 0.269 0.266
item   8 0.280 0.258 0.268 0.274 0.279 0.277 0.269
item   9 0.277 0.248 0.259 0.266 0.273 0.270 0.271

0.273

item  10 0.281 0.258 0.269 0.275 0.280 0.278 0.270

0.280 0.274

item  11 0.282 0.258 0.269 0.275 0.280 0.278 0.271

0.280 0.275 0.281

item  12 0.273 0.243 0.255 0.262 0.269 0.266 0.268

0.269 0.271 0.270  0.271

item  13 0.279 0.251 0.263 0.269 0.276 0.273 0.272

0.276 0.275 0.277  0.278  0.272

MML_Rasch <- lmer(score ~ 0+item+(1|id),

rasch_data, family=binomial(link="logit"))

summary(MML_Rasch)

Generalized linear mixed model fit by the Laplace

approximation

Formula: score ~ 0 + item + (1 | id)
  Data: rasch_data
 AIC  BIC logLik deviance
 6502 6595  -3237     6474
Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 1.9821   1.4079
Number of obs: 5512, groups: id, 424

Fixed effects:
        Estimate Std. Error z value Pr(>|z|)
item   1  0.05819    0.13004   0.447 0.654524
item   2  1.26791    0.14126   8.976  < 2e-16 ***
item   3  0.93972    0.13615   6.902 5.12e-12 ***
item   4  0.70219    0.13345   5.262 1.43e-07 ***
item   5  0.36877    0.13099   2.815 0.004874 **
item   6  0.52699    0.13197   3.993 6.52e-05 ***
item   7 -0.79568    0.13404  -5.936 2.92e-09 ***
item   8  0.38186    0.13106   2.914 0.003572 **
item   9 -0.61867    0.13239  -4.673 2.97e-06 ***
item  10  0.30358    0.13069   2.323 0.020184 *
item  11  0.23870    0.13044   1.830 0.067262 .
item  12 -0.79568    0.13404  -5.936 2.92e-09 ***
item  13 -0.47278    0.13137  -3.599 0.000320 ***
---
Signif. codes:  0 '***' 0.001 '**'

0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
        item1 item2 item3 item4 item5 item6 item7

item8 item9 item10 item11 item12

item   2 0.258
item   3 0.269 0.252
item   4 0.275 0.256 0.265
item   5 0.281 0.258 0.268 0.274
item   6 0.278 0.257 0.267 0.273 0.277
item   7 0.273 0.243 0.255 0.262 0.269 0.266
item   8 0.280 0.258 0.268 0.274 0.279 0.277 0.269
item   9 0.277 0.248 0.259 0.266 0.273 0.270 0.271

0.273

item  10 0.281 0.258 0.269 0.275 0.280 0.278 0.270

0.280 0.274

item  11 0.282 0.258 0.269 0.275 0.280 0.278 0.271

0.280 0.275 0.281

item  12 0.273 0.243 0.255 0.262 0.269 0.266 0.268

0.269 0.271 0.270  0.271

item  13 0.279 0.251 0.263 0.269 0.276 0.273 0.272

0.276 0.275 0.277  0.278  0.272

Dr. Iasonas Lamprianou
Department of Education
The University of Manchester
Oxford Road, Manchester M13 9PL, UK
Tel. 0044  161 275 3485
iasonas.lamprianou at manchester.ac.uk


--- On Sat, 21/2/09,

r-sig-mixed-models-request at r-project.org
<r-sig-mixed-models-request at r-project.org> wrote:

From: r-sig-mixed-models-request at r-project.org

<r-sig-mixed-models-request at r-project.org>

Subject: R-sig-mixed-models Digest, Vol 26, Issue

33

To: r-sig-mixed-models at r-project.org
Date: Saturday, 21 February, 2009, 11:00 AM
Send R-sig-mixed-models mailing list submissions

to

      r-sig-mixed-models at r-project.org

To subscribe or unsubscribe via the World Wide

Web, visit

sessionInfo()

pupils<-read.table("http://www.jkyleroberts.com/rfiles/mlm/HoxData/pupils.txt", header=T)

m0 <- lmer(ACHIEV ~ 1 + (1 | PSCHOOL) + (1 | SSCHOOL), pupils)
summary(m0)

m1 <- lmer(ACHIEV ~ PUPSEX + PUPSES + (1 | PSCHOOL) + (1 | SSCHOOL), pupils)
summary(m1)

m2 <- lmer(ACHIEV~PUPSEX+PUPSES+PDENOM+SDENOM+(1 | PSCHOOL) + (1 | SSCHOOL), pupils)
summary(m2)

anova(m0, m1, m2)

-----Original Message-----
From: r-sig-mixed-models-bounces at r-project.org 
[mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf 
Of Roberts, Kyle
Sent: Tuesday, February 24, 2009 9:29 AM
To: 'Douglas Bates'; r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] AIC in Cross-Classified Model

Dear All,

I was working an example for my class and came upon a 
problem. This example is the cross-classified model presented 
by Hox on p 127 of "Multilevel Analysis." The problem is that 
the AIC and BIC in the summary of the model are not the same 
AIC and BIC when the anova function is used. Any idea as to 
why? I think that the anova function gives the correct AIC, 
but I don't know why. I threw the dataset up on my website if 
anyone wants to replicate.

sessionInfo()

R version 2.7.2 (2008-08-25)
i386-pc-mingw32 

locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United 
States.1252;LC_MONETARY=English_United 
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods 
  base     

other attached packages:
[1] foreign_0.8-29     mlmRev_0.99875-1   lme4_0.999375-27   
Matrix_0.999375-16
[5] lattice_0.17-15   

loaded via a namespace (and not attached):
[1] grid_2.7.2  tools_2.7.2

pupils<-read.table("http://www.jkyleroberts.com/rfiles/mlm/HoxData/pup

ils.txt", header=T)

m0 <- lmer(ACHIEV ~ 1 + (1 | PSCHOOL) + (1 | SSCHOOL), pupils)
summary(m0)

Linear mixed model fit by REML
Formula: ACHIEV ~ 1 + (1 | PSCHOOL) + (1 | SSCHOOL) 
   Data: pupils
  AIC  BIC logLik deviance REMLdev
 2329 2349  -1161     2318    2321
Random effects:
 Groups   Name        Variance Std.Dev.
 PSCHOOL  (Intercept) 0.171900 0.41461
 SSCHOOL  (Intercept) 0.066652 0.25817 
 Residual             0.513128 0.71633 
Number of obs: 1000, groups: PSCHOOL, 50; SSCHOOL, 30

Fixed effects:
            Estimate Std. Error t value
(Intercept)  6.34861    0.07891   80.45

m1 <- lmer(ACHIEV ~ PUPSEX + PUPSES + (1 | PSCHOOL) + (1 | 

SSCHOOL), 

pupils)
summary(m1)

Linear mixed model fit by REML
Formula: ACHIEV ~ PUPSEX + PUPSES + (1 | PSCHOOL) + (1 | SSCHOOL) 
   Data: pupils
  AIC  BIC logLik deviance REMLdev
 2270 2299  -1129     2243    2258
Random effects:
 Groups   Name        Variance Std.Dev.
 PSCHOOL  (Intercept) 0.171530 0.41416
 SSCHOOL  (Intercept) 0.064809 0.25458 
 Residual             0.475236 0.68937 
Number of obs: 1000, groups: PSCHOOL, 50; SSCHOOL, 30

Fixed effects:
            Estimate Std. Error t value
(Intercept)  5.75548    0.10576   54.42
PUPSEXgirl   0.26132    0.04569    5.72
PUPSES       0.11407    0.01612    7.08

Correlation of Fixed Effects:
           (Intr) PUPSEX
PUPSEXgirl -0.254       
PUPSES     -0.641  0.075

m2 <- lmer(ACHIEV~PUPSEX+PUPSES+PDENOM+SDENOM+(1 | PSCHOOL) + (1 | 
SSCHOOL), pupils)
summary(m2)

Linear mixed model fit by REML 
Formula: ACHIEV ~ PUPSEX + PUPSES + PDENOM + SDENOM + (1 | 
PSCHOOL) +      (1 | SSCHOOL) 
   Data: pupils
  AIC  BIC logLik deviance REMLdev
 2273 2312  -1128     2238    2257
Random effects:
 Groups   Name        Variance Std.Dev.
 PSCHOOL  (Intercept) 0.165721 0.40709
 SSCHOOL  (Intercept) 0.058446 0.24176 
 Residual             0.475207 0.68935 
Number of obs: 1000, groups: PSCHOOL, 50; SSCHOOL, 30

Fixed effects:
            Estimate Std. Error t value
(Intercept)  5.51888    0.14275   38.66
PUPSEXgirl   0.26309    0.04568    5.76
PUPSES       0.11352    0.01612    7.04
PDENOMyes    0.20400    0.12623    1.62
SDENOMyes    0.17572    0.09625    1.83

Correlation of Fixed Effects:
           (Intr) PUPSEX PUPSES PDENOM
PUPSEXgirl -0.200                     
PUPSES     -0.466  0.075              
PDENOMyes  -0.526  0.021  0.004       
SDENOMyes  -0.429  0.004 -0.025 -0.014

anova(m0, m1, m2)

Data: pupils
Models:
m0: ACHIEV ~ 1 + (1 | PSCHOOL) + (1 | SSCHOOL)
m1: ACHIEV ~ PUPSEX + PUPSES + (1 | PSCHOOL) + (1 | SSCHOOL)
m2: ACHIEV ~ PUPSEX + PUPSES + PDENOM + SDENOM + (1 | PSCHOOL) + 
m2:     (1 | SSCHOOL)
   Df     AIC     BIC  logLik   Chisq Chi Df Pr(>Chisq)    
m0  4  2325.9  2345.5 -1158.9                              
m1  6  2255.5  2284.9 -1121.7 74.3678      2    < 2e-16 ***
m2  8  2253.5  2292.8 -1118.8  5.9684      2    0.05058 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

*********************************************************
Dr. J. Kyle Roberts
Department of Teaching and Learning
Annette Caldwell Simmons School of Education 
   and Human Development
Southern Methodist University
P.O. Box 750381
Dallas, TX? 75275
214-768-4494
http://www.hlm-online.com/