roberto toro
2008-Sep-14 12:20 UTC
[R] Help please! How to code a mixed-model with 2 within-subject factors using lme or lmer?
Hello,
I'm using aov() to analyse changes in brain volume between males and
females. For every subject (there are 331 in total) I have 8 volume
measurements (4 different brain lobes and 2 different tissues
(grey/white matter)). The data looks like this:
Subject Sex Lobe Tissue Volume
subect1 1 F g 262374
subect1 1 F w 173758
subect1 1 O g 67155
subect1 1 O w 30067
subect1 1 P g 117981
subect1 1 P w 85441
subect1 1 T g 185241
subect1 1 T w 83183
subect2 1 F g 255309
subect2 1 F w 164335
subect2 1 O g 71769
subect2 1 O w 31879
subect2 1 P g 120518
subect2 1 P w 90334
subect2 1 T g 168413
subect2 1 T w 75790
subect3 0 F g 243621
subect3 0 F w 167025
subect3 0 O g 65998
subect3 0 O w 29758
subect3 0 P g 118026
subect3 0 P w 91903
subect3 0 T g 156279
subect3 0 T w 82349
....
I'm trying to see if there is an interaction Sex*Lobe*Tissue. This is
the command I use with aov():
mod1<-aov(Volume~Sex*Lobe*Tissue+Error(Subject/(Lobe*Tissue)),data.vslt)
Subject is a random effect, Sex, Lobe and Tissue are fixed effects;
Sex is an outer factor (between subjects), and Lobe and Tissue are
inner factors (within-subjects); and there is indeed a significant
3-way interaction.
I was told, however, that the results reported by aov() may depend on
the order of the factors
(type I anova), and that is better to use lme() or lmer() with type
II, but I'm struggling to find the right syntaxis...
To begin, how should I write the model using lme() or lmer()??
I tried this with lme():
gvslt<-groupedData(Volume~1|Subject,outer=~Val,inner=list(~Lobe,~Tissue),data=vslt)
mod2<-lme(Volume~Val*Lobe*Tissue,random=~1|Subject,data=gvslt)
but I have interaction terms for every level of Lobe and Tissue, and 8
times the number of DF I should have... (around 331*8 instead of
~331).
Using lmer(), the specification of Subject as a random effect is
straightforward:
mod2<-lmer(Volume~Sex*Lobe*Tissue+(1|Subject),data.vslt)
but I can't figure out the /(Lobe*Tissue) part...
Thank you very much in advance!
roberto
Mark Difford
2008-Sep-14 13:36 UTC
[R] Help please! How to code a mixed-model with 2 within-subject factors using lme or lmer?
Hi Roberto,>> but I can't figure out the /(Lobe*Tissue) part...This type of nesting is easier to do using lmer(). To do it using lme() you have to generate the crossed factor yourself. Do something like this: ## tfac <- with(vslt, interaction(Lobe, Tissue, drop=T)) str(tfac); head(tfac) mod2<-lme(Volume ~ Val*Lobe*Tissue, random = ~1|Subject/tfac, data = vslt) Pre-Scriptum: You can also use ?":" but ?interaction is more flexible and powerful. Regards, Mark. roberto toro wrote:> > Hello, > > I'm using aov() to analyse changes in brain volume between males and > females. For every subject (there are 331 in total) I have 8 volume > measurements (4 different brain lobes and 2 different tissues > (grey/white matter)). The data looks like this: > > Subject Sex Lobe Tissue Volume > subect1 1 F g 262374 > subect1 1 F w 173758 > subect1 1 O g 67155 > subect1 1 O w 30067 > subect1 1 P g 117981 > subect1 1 P w 85441 > subect1 1 T g 185241 > subect1 1 T w 83183 > subect2 1 F g 255309 > subect2 1 F w 164335 > subect2 1 O g 71769 > subect2 1 O w 31879 > subect2 1 P g 120518 > subect2 1 P w 90334 > subect2 1 T g 168413 > subect2 1 T w 75790 > subect3 0 F g 243621 > subect3 0 F w 167025 > subect3 0 O g 65998 > subect3 0 O w 29758 > subect3 0 P g 118026 > subect3 0 P w 91903 > subect3 0 T g 156279 > subect3 0 T w 82349 > .... > > I'm trying to see if there is an interaction Sex*Lobe*Tissue. This is > the command I use with aov(): > > mod1<-aov(Volume~Sex*Lobe*Tissue+Error(Subject/(Lobe*Tissue)),data.vslt) > > Subject is a random effect, Sex, Lobe and Tissue are fixed effects; > Sex is an outer factor (between subjects), and Lobe and Tissue are > inner factors (within-subjects); and there is indeed a significant > 3-way interaction. > > I was told, however, that the results reported by aov() may depend on > the order of the factors > (type I anova), and that is better to use lme() or lmer() with type > II, but I'm struggling to find the right syntaxis... > > To begin, how should I write the model using lme() or lmer()?? > > I tried this with lme(): > > gvslt<-groupedData(Volume~1|Subject,outer=~Val,inner=list(~Lobe,~Tissue),data=vslt) > mod2<-lme(Volume~Val*Lobe*Tissue,random=~1|Subject,data=gvslt) > > but I have interaction terms for every level of Lobe and Tissue, and 8 > times the number of DF I should have... (around 331*8 instead of > ~331). > > Using lmer(), the specification of Subject as a random effect is > straightforward: > > mod2<-lmer(Volume~Sex*Lobe*Tissue+(1|Subject),data.vslt) > > but I can't figure out the /(Lobe*Tissue) part... > > Thank you very much in advance! > roberto > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > >-- View this message in context: http://www.nabble.com/Help-please%21-How-to-code-a-mixed-model-with-2-within-subject-factors-using-lme-or-lmer--tp19479860p19480387.html Sent from the R help mailing list archive at Nabble.com.
roberto toro
2008-Sep-14 14:25 UTC
[R] Help please! How to code a mixed-model with 2 within-subject factors using lme or lmer?
Thanks for answering Mark! I tried with the coding of the interaction you suggested:> tfac<-with(vlt,interaction(Lobe,Tissue,drop=T)) > mod<-lme(Volume~Sex*Lobe*Tissue,random=~1|Subject/tfac,data=vlt)But is it normal that the DF are 2303? DF is 2303 even for the estimate of LobeO that has only 662 values (331 for Tissue=white and 331 for Tissue=grey). I'm not sure either that Sex, Lobe and Tissue are correctly handled.... why are there different estimates called Sex:LobeO, Sex:LobeP, etc, and not just Sex:Lobe as with aov()?. Why there's Tissuew, but not Sex1, for example? Thanks again! roberto ps1. How would you code this with lmer()? ps2. this is part of the output of mod<-lme:> summary(mod)Linear mixed-effects model fit by REML Data: vlt AIC BIC logLik 57528.35 57639.98 -28745.17 Random effects: Formula: ~1 | Subject (Intercept) StdDev: 11294.65 Formula: ~1 | tfac %in% Subject (Intercept) Residual StdDev: 10569.03 4587.472 Fixed effects: Volume ~ Sex * Lobe * Tissue Value Std.Error DF t-value p-value (Intercept) 245224.61 1511.124 2303 162.27963 0.0000 Sex 2800.01 1866.312 329 1.50029 0.1345 LobeO -180794.83 1526.084 2303 -118.46975 0.0000 LobeP -131609.27 1526.084 2303 -86.23984 0.0000 LobeT -73189.97 1526.084 2303 -47.95932 0.0000 Tissuew -72461.05 1526.084 2303 -47.48168 0.0000 Sex:LobeO -663.27 1884.789 2303 -0.35191 0.7249 Sex:LobeP -2146.08 1884.789 2303 -1.13863 0.2550 Sex:LobeT 1379.49 1884.789 2303 0.73191 0.4643 Sex:Tissuew 5387.65 1884.789 2303 2.85849 0.0043 LobeO:Tissuew 43296.99 2158.209 2303 20.06154 0.0000 LobeP:Tissuew 50952.21 2158.209 2303 23.60856 0.0000 LobeT:Tissuew -15959.31 2158.209 2303 -7.39470 0.0000 Sex:LobeO:Tissuew -5228.66 2665.494 2303 -1.96161 0.0499 Sex:LobeP:Tissuew -1482.83 2665.494 2303 -0.55631 0.5781 Sex:LobeT:Tissuew -6037.49 2665.494 2303 -2.26506 0.0236