[R-sig-ME] 3-level binomial model
Andrew Robinson
A.Robinson at ms.unimelb.edu.au
Tue Apr 15 21:22:42 CEST 2008
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/
More information about the R-sig-mixed-models
mailing list