Rasch model
Thank you DH, I used bigsteps and Analysis which both give exactly the same results. The ltm package also gave similar results to the BigSteps and Analysis. However, why does lme4 give different ability estimates, while it gives practically the same difficulty estimates? In case you are interested, I attach the raw data as well as a comparison with the results of various packages in an excel file thanks 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 Sun, 22/2/09, Doran, Harold <HDoran at air.org> wrote:
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
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 ************************************************** _______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models _______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models