[R-sig-ME] Split plot with repeated measures- lme & aov
irantzu.primicia at unavarra.es
irantzu.primicia at unavarra.es
Tue Sep 13 10:28:48 CEST 2011
Dear Mailing List,
I know there have been multiple discussions about this topic before, but
I’m just starting using R and the more I read the more confused I become,
so I would really appreciate any help you could give me.
I have a balanced Split plot design with repeated measures. Thus, I have
DOY (Day of the year) nested in Canopy, nested in Thinning, nested in
Block (random variable). Example of the data (mean MEW by subplot):
Block Thinning Canopy DOY MEW
1 0 M 157 1.25
1 0 M 184 8.4
1 0 M 246 12.8
1 0 P 157 2.2
1 0 P 184 12.6
1 0 P 246 21.2
1 20 M 157 1
1 20 M 184 11.8
1 20 M 246 19
1 20 P 157 6.333333
1 20 P 184 19.4
1 20 P 246 27.4
1 30 M 157 1.2
1 30 M 184 10.4
1 30 M 246 17.8
1 30 P 157 1.6
1 30 P 184 14
1 30 P 246 18.8
2 0 M 157 0.6
2 0 M 184 13
2 0 M 246 18.2
2 0 P 157 0.8
2 0 P 184 13.2
2 0 P 246 18.2
2 20 M 157 1.4
2 20 M 184 13.25
2 20 M 246 23.2
2 20 P 157 1.4
2 20 P 184 12.6
2 20 P 246 22.2
2 30 M 157 1
2 30 M 184 10.4
2 30 M 246 16.4
2 30 P 157 1.2
2 30 P 184 20
2 30 P 246 24.6
3 0 M 157 1.75
3 0 M 184 13
3 0 M 246 18
3 0 P 157 1.2
3 0 P 184 12
3 0 P 246 18
3 20 M 157 2
3 20 M 184 17.6
3 20 M 246 26.4
3 20 P 157 2.25
3 20 P 184 17.6
3 20 P 246 17.75
3 30 M 157 1.4
3 30 M 184 14.6
3 30 M 246 17.75
3 30 P 157 1.6
3 30 P 184 14.4
3 30 P 246 21
My doubts are:
1. If I fit the model (without taking into account the temporal
correlation) with lme and aov, they produce different results for some
variables, but only for the whole-plot (Thinning) and subplot (Canopy)
levels. The codes are:
MEW.aov<-aov(log10(MEW)~Thinning*Canopy*DOY+Error(Block/Thinning/Canopy),data=xylo)
summary(MEW.aov)
Error: Block
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 2 0.061083 0.030542
Error: Block:Thinning
Df Sum Sq Mean Sq F value Pr(>F)
Thinning 2 0.188641 0.094320 17.583 0.01043 *
Residuals 4 0.021457 0.005364
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error: Block:Thinning:Canopy
Df Sum Sq Mean Sq F value Pr(>F)
Canopy 1 0.118541 0.118541 2.8392 0.1430
Thinning:Canopy 2 0.006205 0.003103 0.0743 0.9292
Residuals 6 0.250511 0.041752
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
DOY 2 13.4620 6.7310 517.0589 <2e-16 ***
Thinning:DOY 4 0.0360 0.0090 0.6911 0.6053
Canopy:DOY 2 0.0214 0.0107 0.8213 0.4519
Thinning:Canopy:DOY 4 0.0518 0.0129 0.9944 0.4297
Residuals 24 0.3124 0.0130
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
MEW.lme<-lme(log10(MEW)~Thinning*Canopy*DOY,random=~1|Block/Thinning/Canopy,data=xylo,method="ML")
anova(MEW.lme)
numDF denDF F-value p-value
(Intercept) 1 24 1310.8941 <.0001
Thinning 2 4 3.4681 0.1338
Canopy 1 6 4.3586 0.0819
DOY 2 24 517.0588 <.0001
Thinning:Canopy 2 6 0.1141 0.8941
Thinning:DOY 4 24 0.6911 0.6053
Canopy:DOY 2 24 0.8213 0.4519
Thinning:Canopy:DOY 4 24 0.9944 0.4297
2. Once I have the correct lme code, I think it’s better to not aggregate
the data by subplots (like when using aov), isn’t it? I mean, if for
example I measure the diameter of the same trees repeatedly and I have the
design: Block>Thinning>Canopy>Tree code> DOY then, I think that it
would be better to avoid the loos of the inter-tree variability and
include the corrector for temporal correlation. The tree CODEs are unique
and in this data frame there are some missing data:
Block Thinning Canopy CODE DOY MEW
1 20 M A1074 157 1
1 20 M A1074 184 7
1 20 M A1074 246 29
1 20 P A1090 157 14
1 20 P A1090 184 18
1 20 P A1090 246 32
1 20 P A1120 157 3
1 20 P A1120 184 9
1 20 P A1120 246 16
1 20 M A1164 157 2
1 20 M A1164 184 20
1 20 M A1164 246 20
1 20 M A1202 157 1
1 20 M A1202 184 14
1 20 M A1202 246 30
1 20 M A1208 157 0
1 20 M A1208 184 2
1 20 M A1208 246 4
1 20 P A1237 157 NA
1 20 P A1237 184 29
1 20 P A1237 246 43
1 20 P A1247 157 NA
1 20 P A1247 184 22
1 20 P A1247 246 15
1 20 P A1445 157 2
1 20 P A1445 184 19
1 20 P A1445 246 31
1 20 M A1478 157 1
1 20 M A1478 184 16
1 20 M A1478 246 12
1 30 P A2024 157 1
1 30 P A2024 184 23
1 30 P A2024 246 19
1 30 P A2075 157 3
1 30 P A2075 184 19
1 30 P A2075 246 31
1 30 P A2154 157 4
MEW.lme<-lme(log10(MEW)~Thinning*Canopy*DOY, random=~1|CODE,
correlation=corAR1(form=~as.numeric(DOY)| CODE), data=xylo,
na.action=na.omit)
anova(MEW.lme)
numDF denDF F-value p-value
(Intercept) 1 237 451.6690 <.0001
Thinning 2 84 2.4694 0.0908
Canopy 1 84 2.7178 0.1030
DOY 3 237 252.3011 <.0001
Thinning:Canopy 2 84 0.1822 0.8338
Thinning:DOY 6 237 1.4853 0.1839
Canopy:DOY 3 237 1.1750 0.3199
Thinning:Canopy:DOY 6 237 0.7248 0.6300
But, should I include the “Block” variable as well in the random effects?
Thank you very much for your help,
Irantzu Primicia
PhD Student
Area de Ecología
Departamento de Ciencias del Medio Natural
Universidad Pública de Navarra
More information about the R-sig-mixed-models
mailing list