[R-sig-ME] Rasch model

Iasonas Lamprianou lamprianou at yahoo.com
Sat Feb 21 12:31:44 CET 2009


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
> 	https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> or, via email, send a message with subject or body
> 'help' to
> 	r-sig-mixed-models-request at r-project.org
> 
> You can reach the person managing the list at
> 	r-sig-mixed-models-owner at r-project.org
> 
> When replying, please edit your Subject line so it is more
> specific
> than "Re: Contents of R-sig-mixed-models
> digest..."
> 
> 
> Today's Topics:
> 
>    1. Re: Power calculations for random effect models (Greg
> Snow)
> 
> 
> ----------------------------------------------------------------------
> 
> Message: 1
> Date: Fri, 20 Feb 2009 18:12:52 -0700
> From: Greg Snow <Greg.Snow at imail.org>
> Subject: Re: [R-sig-ME] Power calculations for random
> effect models
> To: "christos at nuverabio.com"
> <christos at nuverabio.com>,
> 	"r-sig-mixed-models at r-project.org"
> <r-sig-mixed-models at r-project.org>
> Message-ID:
> 	<B37C0A15B8FB3C468B5BC7EBC7DA14CC61CA3E2EBE at LP-EXMBVS10.CO.IHC.COM>
> Content-Type: text/plain; charset="us-ascii"
> 
> Since you already have a function to simulate the data, it
> should be fairly simple to use simulations to compare the
> methods you suggest below, or to do a simple test based on
> the likelihood ratio, but simulating the cutoff rather than
> depending on the chi-squared approximation.
> 
> See:
> http://finzi.psych.upenn.edu/R/Rhelp08/archive/156522.html
> 
> for some examples of assessing these types of questions
> using simulation (I have an updated version of the last
> example if you want it as well).
> 
> Hope this helps,
> 
> -- 
> Gregory (Greg) L. Snow Ph.D.
> Statistical Data Center
> Intermountain Healthcare
> greg.snow at imail.org
> 801.408.8111
> 
> 
> > -----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 Christos
> Hatzis
> > Sent: Friday, February 20, 2009 11:01 AM
> > To: r-sig-mixed-models at r-project.org
> > Subject: [R-sig-ME] Power calculations for random
> effect models
> > 
> > Hello,
> > 
> > I have a nested random effects model of the form
> > 
> > Yijk = M + Ai + Bj(i) + Ek(ij)
> > 
> > where biopsies (B) are nested within persons (A) and
> arrays (E) are
> > nested
> > within biopsies.
> > I am interested in estimating the power of a study of
> a given size to
> > determine whether the variance associated with B is
> trivial, i.e. H0:
> > var_b
> > = 0 vs Ha: var_b > 0 at a fixed Type-I error rate.
> > 
> > I have written a function to simulate data from this
> model and used
> > lmer to
> > estimate the random effects as shown below.  What is
> the recommended
> > way of
> > going about testing the null hypothesis?  One approach
> would be to use
> > the
> > estimated standard deviation of the random effect and
> assuming
> > normality to
> > test whether the appropriate CI contains zero. Another
> would be to use
> > the
> > theoretical chi-square distribution for var_b, but
> that would require
> > the
> > appropriate degrees of freedom.  Or to use mcmc to
> estimate the
> > distribution
> > of var_b and use this distribution for inference.
> > 
> > I would think that the mcmc approach is the
> recommended one, but I
> > would
> > appreciate any advise on this. In this case, I have
> tried to run the
> > MCMC
> > simulation (which runs fine), but have not been able
> yet to figure out
> > how
> > to use the results of MCMC to test the above
> hypothesis.  Any hints or
> > pointing out materials that explain how to use MCMC
> for inference on
> > random
> > effects will be very much appreciated.
> > 
> > Thank you.
> > -Christos
> > 
> > Christos Hatzis, Ph.D.
> > Nuvera Biosciences, Inc.
> > 400 West Cummings Park
> > Suite 5350
> > Woburn, MA 01801
> > Tel: 781-938-3830
> > www.nuverabio.com <http://www.nuverabio.com/>
> > 
> > > fake.dt <- tumor.heter.dt(10, 3, 2, k=0, sa=4,
> sb=.6, se=.2)
> > > fake.dt
> >    person biopsy array      resp bioWper
> > 1       1      1     1  4.308434     1:1
> > 2       1      1     2  4.293186     1:1
> > 3       1      2     1  5.503841     1:2
> > 4       1      2     2  5.640362     1:2
> > 5       1      3     1  5.579461     1:3
> > 6       1      3     2  5.416201     1:3
> > 7       2      1     1 12.479513     2:1
> > 8       2      1     2 12.311426     2:1
> > 9       2      2     1 13.283566     2:2
> > 10      2      2     2 13.138276     2:2
> > 11      2      3     1 12.954277     2:3
> > 12      2      3     2 13.059925     2:3
> > 13      3      1     1 11.649726     3:1
> > 14      3      1     2 11.694472     3:1
> > 15      3      2     1 12.342050     3:2
> > 16      3      2     2 12.316214     3:2
> > 17      3      3     1 12.762695     3:3
> > 18      3      3     2 12.451338     3:3
> > 19      4      1     1 17.168248     4:1
> > 20      4      1     2 17.573315     4:1
> > 21      4      2     1 16.885176     4:2
> > 22      4      2     2 16.819536     4:2
> > 23      4      3     1 15.924120     4:3
> > 24      4      3     2 16.038491     4:3
> > 25      5      1     1  7.981279     5:1
> > 26      5      1     2  7.777929     5:1
> > 27      5      2     1  6.701801     5:2
> > 28      5      2     2  6.978284     5:2
> > 29      5      3     1  7.136417     5:3
> > 30      5      3     2  7.378498     5:3
> > 31      6      1     1 21.499796     6:1
> > 32      6      1     2 21.397173     6:1
> > 33      6      2     1 22.551116     6:2
> > 34      6      2     2 22.521006     6:2
> > 35      6      3     1 21.027998     6:3
> > 36      6      3     2 21.416872     6:3
> > 37      7      1     1 13.312857     7:1
> > 38      7      1     2 13.559062     7:1
> > 39      7      2     1 13.753016     7:2
> > 40      7      2     2 13.608642     7:2
> > 41      7      3     1 13.556446     7:3
> > 42      7      3     2 13.400672     7:3
> > 43      8      1     1 12.758548     8:1
> > 44      8      1     2 12.486574     8:1
> > 45      8      2     1 13.388409     8:2
> > 46      8      2     2 13.263029     8:2
> > 47      8      3     1 12.991308     8:3
> > 48      8      3     2 12.962116     8:3
> > 49      9      1     1 11.420214     9:1
> > 50      9      1     2 11.308010     9:1
> > 51      9      2     1 13.186774     9:2
> > 52      9      2     2 12.778966     9:2
> > 53      9      3     1 11.625079     9:3
> > 54      9      3     2 11.637015     9:3
> > 55     10      1     1  9.108635    10:1
> > 56     10      1     2  8.895658    10:1
> > 57     10      2     1  9.718046    10:2
> > 58     10      2     2  9.552510    10:2
> > 59     10      3     1  8.794893    10:3
> > 60     10      3     2  9.047379    10:3
> > > heter.lmer <- lmer(resp ~ (1 | person) + (1 |
> bioWper), fake.dt)
> > > heter.lmer
> > Linear mixed model fit by REML
> > Formula: resp ~ (1 | person) + (1 | bioWper)
> >    Data: fake.dt
> >    AIC   BIC logLik deviance REMLdev
> >  98.24 106.6 -45.12    92.85   90.24
> > Random effects:
> >  Groups   Name        Variance  Std.Dev.
> >  bioWper  (Intercept)  0.312327 0.55886
> >  person   (Intercept) 21.780993 4.66701
> >  Residual              0.020633 0.14364
> > Number of obs: 60, groups: bioWper, 30; person, 10
> > 
> > Fixed effects:
> >             Estimate Std. Error t value
> > (Intercept)   12.368      1.479    8.36
> > 
> > 
> > > VarCorr(heter.lmer)[["bioWper"]]
> >             (Intercept)
> > (Intercept)   0.3123273
> > attr(,"stddev")
> > (Intercept)
> >   0.5588625
> > attr(,"correlation")
> >             (Intercept)
> > (Intercept)           1
> > 
> > > heter.mcmc <- mcmcsamp(heter.lmer, n=1000)
> > > str(heter.mcmc)
> > Formal class 'merMCMC' [package
> "lme4"] with 9 slots
> >   ..@ Gp      : int [1:3] 0 30 40
> >   ..@ ST      : num [1:2, 1:1000] 3.89 32.49 1.44
> 27.06 1.24 ...
> >   ..@ call    : language lmer(formula = resp ~ (1 |
> person) + (1 |
> > bioWper),
> > data = fake.dt)
> >   ..@ deviance: num [1:1000] 92.9 92.9 114.5 119.1 134
> ...
> >   ..@ dims    : Named int [1:18] 2 60 1 40 1 2 0 1 2 5
> ...
> >   .. ..- attr(*, "names")= chr [1:18]
> "nt" "n" "p" "q" ...
> >   ..@ fixef   : num [1, 1:1000] 12.4 14.1 11.9 12.6
> 12.6 ...
> >   .. ..- attr(*, "dimnames")=List of 2
> >   .. .. ..$ : chr "(Intercept)"
> >   .. .. ..$ : NULL
> >   ..@ nc      : int [1:2] 1 1
> >   ..@ ranef   : num[1:40, 0 ]
> >   ..@ sigma   : num [1, 1:1000] 0.144 0.128 0.126
> 0.212 0.228 ...
> > 
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> >
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
> 
> 
> ------------------------------
> 
> _______________________________________________
> R-sig-mixed-models mailing list
> R-sig-mixed-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
> 
> End of R-sig-mixed-models Digest, Vol 26, Issue 33
> **************************************************







More information about the R-sig-mixed-models mailing list