3-level binomial model
Hi Iasonas, a few observations, questions, suggestions: 1) it's not clear to me what you are trying to do, so I can't comment on the technique. What is your goal for estimating VPCschool and VPCpupil? Those percentages don't make sense to me. I would ordinarily assume that if such a decomposition is possible, then VPCschool = VARschool/(VARschool+VARpupil+3.29) and VPCpupil = VARpupil/(VARschool+VARpupil+3.29) 2) does the MlWin manual not provide a citation for their suggestion? If so you could use that. If not, you could contact the authors and explain the problem. Perhaps they can help. 3) do the reviewers not provide more guidance? "not sure if this is the best way" is not the same as "sure that this is not the best way", and in any case, sometimes the not-the-best way is just fine. Can you contact the Associate Editor via the journal and ask for more guidance? Eg emphasize that you are willing to use a different technique but you are unaware of any, and if the reviewers could provide a pointer to a suitable technique it would be very helpful. 4) 3 months is a doddle. I just got a paper accepted that I started 7 years ago. I hope that these thoughts help, Andrew
On Mon, Apr 14, 2008 at 11:34:34PM -0700, Iasonas Lamprianou wrote:
Hi all of you,
I have a problem. I run a model having pupils nested within schools. Each pupil took a number of different tests e.g. maths, science etc. All in all I had 7 different tests, most of the pupils took all of them. So, my model is school:pupil:tests (7 tests, the same for all pupils). The dependent variable is binomial (0/1). I got the results in lme4. This was the model:
ABERRANT/NON-ABERRANT ~ CONSTANT + PAPERFixed + SCHOOLRandom + PUPILRandom
I computed the school-level and the pupil-level variance like that (as described for 2-level models in MlWin manual): I assumed that my dependent variable is based on a continuous unobserved variable (perfectly valid according to my theoretical model). Therefore, eijk follows a logistic distribution with variance pi2/3=3.29. So,
VPCschool=VARschool/(VARschool+3.29)= 0.17577/(0.17577+3.29)=6.4% and VPCpupil=VPCpupil /(VPCpupil+3.29)=0.19977/(0.19977+3.29)=7.3%. The reviewers of my paper are not sure if this is the best way to do it. They may reject my paper and I worry because I have spent 3months!!!! writing it. Any ideas to support my method or to use a better one? These are the results of lme4:
AIC BIC logLik deviance
6372 6440 -3177 6354
Random effects:
Groups Name Variance Std.Dev.
Pupils (Intercept) 0.19977 0.44696
Schools (Intercept) 0.17577 0.41925
number of observations: 14829, groups: Pupils= 2230; Schools= 151
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.04953 0.10565 -28.866 <0.001 ***
PAPER[Reading] 0.37060 0.13023 2.846 0.00443 **
PAPER[Maths A] 0.17508 0.13536 1.293 0.19585
PAPER[Maths B] -0.07881 0.14298 -0.551 0.58148
PAPER[Mental Maths] 0.01880 0.14130 0.133 0.89418
PAPER[Science A] 0.10525 0.13842 0.760 0.44703
PAPER[Science B] -0.18359 0.14958 -1.227 0.21970
---
Signif. codes: 0 ???***??? 0.001 ???**??? 0.01 ???*??? 0.05 ???.??? 0.1 ??? ??? 1
Correlation of Fixed Effects:
(Intr) Read MathA MathB MathMental SciA
PAPER[Read] -0.715
PAPER[MatA] -0.688 0.558
PAPER[MatB] -0.652 0.528 0.510
PAPER[MatM] -0.658 0.534 0.514 0.487
PAPER[SciA] -0.672 0.545 0.525 0.497 0.504
PAPER[SciB] -0.622 0.505 0.486 0.460 0.467 0.476
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
----- Original Message ----
From: "r-sig-mixed-models-request at r-project.org" <r-sig-mixed-models-request at r-project.org>
To: r-sig-mixed-models at r-project.org
Sent: Friday, 4 April, 2008 1:00:02 PM
Subject: R-sig-mixed-models Digest, Vol 16, Issue 12
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: random effects specification (Douglas Bates) 2. Re: random effects specification (Sebastian P. Luque) ---------------------------------------------------------------------- Message: 1 Date: Thu, 3 Apr 2008 17:32:40 -0500 From: "Douglas Bates" <bates at stat.wisc.edu> Subject: Re: [R-sig-ME] random effects specification To: "Sebastian P. Luque" <spluque at gmail.com> Cc: r-sig-mixed-models at r-project.org Message-ID: <40e66e0b0804031532q12880b01mcba6173fc3e32c10 at mail.gmail.com> Content-Type: text/plain; charset=ISO-8859-1 On Thu, Apr 3, 2008 at 3:45 PM, Sebastian P. Luque <spluque at gmail.com> wrote: Hi, In the past I've used lme to fit simple mixed models to longitudinal data (growth), but now I'm trying to learn lmer and its different syntax to fit a more complex model. I have a structure with subjects (id, random factor) exposed to 4 different treatments and a continuous response variable is measured (n). The subjects come from 2 different communities, so it's a nested design very much like the Oats data in the nlme package. The interest is in the fixed effects of community and treatment, and their interaction, so I thought this could be modelled in lmer with this call: lmer(n ~ treatment + community + (1 | id/treatment), mydata) but got this error: Error in lmerFactorList(formula, mf, fltype) : number of levels in grouping factor(s) 'treatment:id' is too large Am I using the right formula here? Thanks. It seems that the observations are indexed by subject and treatment so the number of levels in the factor treatment:id equals the number of observations. You can't estimate a variance for such a term and also estimate a residual variance. I would start with n ~ treatment * community +(1|id) ------------------------------ Message: 2 Date: Thu, 03 Apr 2008 18:32:23 -0500 From: "Sebastian P. Luque" <spluque at gmail.com> Subject: Re: [R-sig-ME] random effects specification To: r-sig-mixed-models at r-project.org Message-ID: <878wzuzf94.fsf at patagonia.sebmags.homelinux.org> Content-Type: text/plain; charset=us-ascii On Thu, 3 Apr 2008 17:32:40 -0500, "Douglas Bates" <bates at stat.wisc.edu> wrote: [...] It seems that the observations are indexed by subject and treatment so the number of levels in the factor treatment:id equals the number of observations. You can't estimate a variance for such a term and also estimate a residual variance. I would start with n ~ treatment * community +(1|id) Yes, the observations are indexed by subject and treatment in the sense that id levels are the same within treatments of the same community, but are different among communities. This is a subset of the data: ---<---------------cut here---------------start-------------->--- id community treatment n 1 A 1 13.93 2 A 1 14.42 3 A 1 13.56 1 A 2 14.61 2 A 2 14.74 3 A 2 15.59 1 A 3 13.95 2 A 3 15.21 3 A 3 14.51 1 A 4 13.61 2 A 4 14.99 3 A 4 15.13 4 B 1 14.79 5 B 1 13.41 6 B 1 14.71 4 B 2 14.69 5 B 2 13.46 6 B 2 14.28 4 B 3 14.30 5 B 3 13.18 6 B 3 13.58 4 B 4 14.54 5 B 4 13.25 6 B 4 14.09 ---<---------------cut here---------------end---------------->--- Of course, there are many more individuals, but the levels of id differ among communities, and are the same among treatments. lmer did converge rapidly with your suggested formula though: ---<---------------cut here---------------start-------------->--- Linear mixed-effects model fit by REML Formula: n ~ treatment * community + (1 | id) Data: isotope.m.ph AIC BIC logLik MLdeviance REMLdeviance 450 481 -216 410 432 Random effects: Groups Name Variance Std.Dev. id (Intercept) 0.193 0.439 Residual 0.232 0.481 number of obs: 240, groups: id, 61 Fixed effects: Estimate Std. Error t value (Intercept) 14.9748 0.1170 128.0 treatment2 0.0884 0.1222 0.7 treatment3 -0.2829 0.1222 -2.3 treatment4 0.3568 0.1222 2.9 communitysanikiluaq -0.5749 0.1678 -3.4 treatment2:communitysanikiluaq -0.2471 0.1763 -1.4 treatment3:communitysanikiluaq -0.7479 0.1763 -4.2 treatment4:communitysanikiluaq -0.6169 0.1763 -3.5 Correlation of Fixed Effects: (Intr) trtmn2 trtmn3 trtmn4 cmmnty trtm2: trtm3: treatment2 -0.522 treatment3 -0.522 0.500 treatment4 -0.522 0.500 0.500 cmmntysnklq -0.697 0.364 0.364 0.364 trtmnt2:cmm 0.362 -0.693 -0.347 -0.347 -0.524 trtmnt3:cmm 0.362 -0.347 -0.693 -0.347 -0.524 0.503 trtmnt4:cmm 0.362 -0.347 -0.347 -0.693 -0.524 0.503 0.503 ---<---------------cut here---------------end---------------->--- However, I don't understand how (1 | id) accounts for treatment being nested within community. Maybe it's time for me to re-read some more examples from "Mixed-effects models in S and S-plus". Thanks. Cheers, -- Seb ------------------------------ _______________________________________________ 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 16, Issue 12 ************************************************** ___________________________________________________________ Yahoo! For Good helps you make a difference _______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
Andrew Robinson Department of Mathematics and Statistics Tel: +61-3-8344-6410 University of Melbourne, VIC 3010 Australia Fax: +61-3-8344-4599 http://www.ms.unimelb.edu.au/~andrewpr http://blogs.mbs.edu/fishing-in-the-bay/