[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