[R-sig-ME] 3-level binomial model
Iasonas Lamprianou
lamprianou at yahoo.com
Tue Apr 15 08:34:34 CEST 2008
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
More information about the R-sig-mixed-models
mailing list