Nested mixed models
Markus Löw
28 July 2016
This is my attempt at providing an overview on statistical methods that can be used for our research.
Know your data
Before any analysis can take place, it is important to know the data. Let’s explore the data. This is an example how to analyse yield data [g m^2]. Make sure that all data that are needed for the analysis are available!
# check the data
head(df)## Plot.Order CO2.TOS.Irr.Cultivar.DATABASE TrialID Year RingID PlotID
## 3 1 aCO2 TOS1 Sup YitpiN0 Horsham 2007 1 A
## 27 2 aCO2 TOS2 Sup YitpiN0 Horsham 2008 1 D
## 33 3 aCO2 TOS1 Sup Janz Horsham 2007 1 F
## 36 3 aCO2 TOS2 Sup Janz Horsham 2008 1 F
## 54 4 aCO2 TOS2 Sup YitpiN0 Horsham 2009 1 D
## 63 5 aCO2 TOS2 Rain YitpiN0 Horsham 2008 1 I
## HalfRingID RingTrt CO2 Irrigation Nitrogen TOS PlotTrt Sampling
## 3 W aCO2 Sup aCO2 Sup N0 TOS1 YitpiN0 growth
## 27 W TOS 2 aCO2 aCO2 Sup N0 TOS2 YitpiN0 Growth
## 33 W aCO2 Sup aCO2 Sup N0 TOS1 Janz growth
## 36 W TOS 2 aCO2 aCO2 Sup N0 TOS2 JanzN0 Growth
## 54 W TOS 2 aCO2 aCO2 Sup N0 TOS2 YitpiN0 Growth
## 63 E TOS 2 aCO2 aCO2 Rain N0 TOS2 YitpiN0 Growth
## Cultivar Stage Yield YearRing Environment
## 3 Yitpi DC90 327.88162 2007.1 Horsham.2007.Sup.TOS1
## 27 Yitpi DC90 163.65654 2008.1 Horsham.2008.Sup.TOS2
## 33 Janz DC90 299.06542 2007.1 Horsham.2007.Sup.TOS1
## 36 Janz DC90 63.30607 2008.1 Horsham.2008.Sup.TOS2
## 54 Yitpi DC90 140.89744 2009.1 Horsham.2009.Sup.TOS2
## 63 Yitpi DC90 189.84813 2008.1 Horsham.2008.Rain.TOS2
## Ord.Environment
## 3 Horsham.2007.Sup.TOS1
## 27 Horsham.2008.Sup.TOS2
## 33 Horsham.2007.Sup.TOS1
## 36 Horsham.2008.Sup.TOS2
## 54 Horsham.2009.Sup.TOS2
## 63 Horsham.2008.Rain.TOS2tail(df)## Plot.Order CO2.TOS.Irr.Cultivar.DATABASE TrialID Year RingID PlotID
## 1461 350 eCO2 TOS1 Rain YitpiN0 Horsham 2009 15 N
## 1470 355 eCO2 TOS1 Rain Janz Horsham 2009 15 S
## 1497 367 aCO2 TOS2 Sup Janz Horsham 2009 16 G
## 1506 371 aCO2 TOS2 Sup YitpiN0 Horsham 2009 16 K
## 1521 378 aCO2 TOS2 Rain YitpiN0 Horsham 2009 16 R
## 1530 382 aCO2 TOS2 Rain Janz Horsham 2009 16 V
## HalfRingID RingTrt CO2 Irrigation Nitrogen TOS PlotTrt Sampling
## 1461 E TOS 1 eCO2 eCO2 Rain N0 TOS1 YitpiN0 Growth
## 1470 E TOS 1 eCO2 eCO2 Rain N0 TOS1 Janz Growth
## 1497 W TOS 2 aCO2 aCO2 Sup N0 TOS2 Janz Growth
## 1506 W TOS 2 aCO2 aCO2 Sup N0 TOS2 YitpiN0 Growth
## 1521 E TOS 2 aCO2 aCO2 Rain N0 TOS2 YitpiN0 Growth
## 1530 E TOS 2 aCO2 aCO2 Rain N0 TOS2 Janz Growth
## Cultivar Stage Yield YearRing Environment
## 1461 Yitpi DC90 381.6026 2009.15 Horsham.2009.Rain.TOS1
## 1470 Janz DC90 283.5385 2009.15 Horsham.2009.Rain.TOS1
## 1497 Janz DC90 137.5897 2009.16 Horsham.2009.Sup.TOS2
## 1506 Yitpi DC90 101.9103 2009.16 Horsham.2009.Sup.TOS2
## 1521 Yitpi DC90 145.9615 2009.16 Horsham.2009.Rain.TOS2
## 1530 Janz DC90 183.0385 2009.16 Horsham.2009.Rain.TOS2
## Ord.Environment
## 1461 Horsham.2009.Rain.TOS1
## 1470 Horsham.2009.Rain.TOS1
## 1497 Horsham.2009.Sup.TOS2
## 1506 Horsham.2009.Sup.TOS2
## 1521 Horsham.2009.Rain.TOS2
## 1530 Horsham.2009.Rain.TOS2summary(df)## Plot.Order CO2.TOS.Irr.Cultivar.DATABASE TrialID
## Min. : 1.00 aCO2 TOS1 Rain Janz : 12 Horsham:192
## 1st Qu.: 32.75 aCO2 TOS1 Rain YitpiN0: 12 Walpeup: 0
## Median : 64.50 aCO2 TOS1 Sup Janz : 12
## Mean : 96.36 aCO2 TOS1 Sup YitpiN0 : 12
## 3rd Qu.: 97.50 aCO2 TOS2 Rain Janz : 12
## Max. :382.00 aCO2 TOS2 Rain YitpiN0: 12
## (Other) :120
## Year RingID PlotID HalfRingID RingTrt
## 2007:64 Min. : 1.00 D :32 E:96 TOS 1 aCO2:32
## 2008:64 1st Qu.: 4.75 I :25 W:96 TOS 1 eCO2:32
## 2009:64 Median : 8.50 A :23 TOS 2 aCO2:32
## Mean : 8.50 G :23 TOS 2 eCO2:32
## 3rd Qu.:12.25 F :22 aCO2 Rain :16
## Max. :16.00 L :22 aCO2 Sup :16
## (Other):45 (Other) :32
## CO2 Irrigation Nitrogen TOS PlotTrt Sampling
## aCO2:96 Rain:96 N0:192 TOS1:96 YitpiN0 :96 : 0
## eCO2:96 Sup :96 N+: 0 TOS2:96 Janz :64 growth: 64
## JanzN0 :32 Growth:128
## : 0
## Drysdale: 0
## Gladius : 0
## (Other) : 0
## Cultivar Stage Yield YearRing
## Janz :96 DC90:192 Min. : 19.52 2007.1 : 4
## Yitpi:96 1st Qu.:176.02 2008.1 : 4
## Median :256.90 2009.1 : 4
## Mean :257.12 2007.2 : 4
## 3rd Qu.:328.90 2008.2 : 4
## Max. :556.17 2009.2 : 4
## NA's :1 (Other):168
## Environment Ord.Environment
## Horsham.2007.Rain.TOS1:16 Horsham.2009.Rain.TOS2:16
## Horsham.2008.Rain.TOS1:16 Horsham.2008.Rain.TOS2:16
## Horsham.2009.Rain.TOS1:16 Horsham.2009.Sup.TOS2 :16
## Horsham.2007.Sup.TOS1 :16 Horsham.2008.Sup.TOS2 :16
## Horsham.2008.Sup.TOS1 :16 Horsham.2007.Rain.TOS2:16
## Horsham.2009.Sup.TOS1 :16 Horsham.2007.Sup.TOS2 :16
## (Other) :96 (Other) :96This simple overview already reveals that there are several missing values in the yield data. This might cause problems with the statistical analysis later.
# Frequency tables of samples
# check the number of samples
with(df, table(Year, TOS, Irrigation, CO2, Cultivar))## , , Irrigation = Rain, CO2 = aCO2, Cultivar = Janz
##
## TOS
## Year TOS1 TOS2
## 2007 4 4
## 2008 4 4
## 2009 4 4
##
## , , Irrigation = Sup, CO2 = aCO2, Cultivar = Janz
##
## TOS
## Year TOS1 TOS2
## 2007 4 4
## 2008 4 4
## 2009 4 4
##
## , , Irrigation = Rain, CO2 = eCO2, Cultivar = Janz
##
## TOS
## Year TOS1 TOS2
## 2007 4 4
## 2008 4 4
## 2009 4 4
##
## , , Irrigation = Sup, CO2 = eCO2, Cultivar = Janz
##
## TOS
## Year TOS1 TOS2
## 2007 4 4
## 2008 4 4
## 2009 4 4
##
## , , Irrigation = Rain, CO2 = aCO2, Cultivar = Yitpi
##
## TOS
## Year TOS1 TOS2
## 2007 4 4
## 2008 4 4
## 2009 4 4
##
## , , Irrigation = Sup, CO2 = aCO2, Cultivar = Yitpi
##
## TOS
## Year TOS1 TOS2
## 2007 4 4
## 2008 4 4
## 2009 4 4
##
## , , Irrigation = Rain, CO2 = eCO2, Cultivar = Yitpi
##
## TOS
## Year TOS1 TOS2
## 2007 4 4
## 2008 4 4
## 2009 4 4
##
## , , Irrigation = Sup, CO2 = eCO2, Cultivar = Yitpi
##
## TOS
## Year TOS1 TOS2
## 2007 4 4
## 2008 4 4
## 2009 4 4In principle, the number of samples is balanced across years, and the number of samples reflects the number of repeats in the experiment. But what about the actual values - there could be missing values for some samples!
library(plyr) # loading the plyr library for split, apply, combine processing
myNA <- ddply(df,
.(Year, TOS, Irrigation, CO2, Cultivar),
summarise,
No_of_NA = sum(is.na(Yield)))
myNA## Year TOS Irrigation CO2 Cultivar No_of_NA
## 1 2007 TOS1 Rain aCO2 Janz 0
## 2 2007 TOS1 Rain aCO2 Yitpi 1
## 3 2007 TOS1 Rain eCO2 Janz 0
## 4 2007 TOS1 Rain eCO2 Yitpi 0
## 5 2007 TOS1 Sup aCO2 Janz 0
## 6 2007 TOS1 Sup aCO2 Yitpi 0
## 7 2007 TOS1 Sup eCO2 Janz 0
## 8 2007 TOS1 Sup eCO2 Yitpi 0
## 9 2007 TOS2 Rain aCO2 Janz 0
## 10 2007 TOS2 Rain aCO2 Yitpi 0
## 11 2007 TOS2 Rain eCO2 Janz 0
## 12 2007 TOS2 Rain eCO2 Yitpi 0
## 13 2007 TOS2 Sup aCO2 Janz 0
## 14 2007 TOS2 Sup aCO2 Yitpi 0
## 15 2007 TOS2 Sup eCO2 Janz 0
## 16 2007 TOS2 Sup eCO2 Yitpi 0
## 17 2008 TOS1 Rain aCO2 Janz 0
## 18 2008 TOS1 Rain aCO2 Yitpi 0
## 19 2008 TOS1 Rain eCO2 Janz 0
## 20 2008 TOS1 Rain eCO2 Yitpi 0
## 21 2008 TOS1 Sup aCO2 Janz 0
## 22 2008 TOS1 Sup aCO2 Yitpi 0
## 23 2008 TOS1 Sup eCO2 Janz 0
## 24 2008 TOS1 Sup eCO2 Yitpi 0
## 25 2008 TOS2 Rain aCO2 Janz 0
## 26 2008 TOS2 Rain aCO2 Yitpi 0
## 27 2008 TOS2 Rain eCO2 Janz 0
## 28 2008 TOS2 Rain eCO2 Yitpi 0
## 29 2008 TOS2 Sup aCO2 Janz 0
## 30 2008 TOS2 Sup aCO2 Yitpi 0
## 31 2008 TOS2 Sup eCO2 Janz 0
## 32 2008 TOS2 Sup eCO2 Yitpi 0
## 33 2009 TOS1 Rain aCO2 Janz 0
## 34 2009 TOS1 Rain aCO2 Yitpi 0
## 35 2009 TOS1 Rain eCO2 Janz 0
## 36 2009 TOS1 Rain eCO2 Yitpi 0
## 37 2009 TOS1 Sup aCO2 Janz 0
## 38 2009 TOS1 Sup aCO2 Yitpi 0
## 39 2009 TOS1 Sup eCO2 Janz 0
## 40 2009 TOS1 Sup eCO2 Yitpi 0
## 41 2009 TOS2 Rain aCO2 Janz 0
## 42 2009 TOS2 Rain aCO2 Yitpi 0
## 43 2009 TOS2 Rain eCO2 Janz 0
## 44 2009 TOS2 Rain eCO2 Yitpi 0
## 45 2009 TOS2 Sup aCO2 Janz 0
## 46 2009 TOS2 Sup aCO2 Yitpi 0
## 47 2009 TOS2 Sup eCO2 Janz 0
## 48 2009 TOS2 Sup eCO2 Yitpi 0This reveals that there is one missing value in the year 2007. This also means that our experiment is unbalanced. There is one treatment combination that has less valid samples! A full data set would be preferred. There are ways to deal with missing data, and in the worst case, gap-filling options exist.
Figures allow a better overview on the data
# boxplot provide an excellent overview of data as they show distribution and potential outliers.
library(ggplot2) # loading the Grammar of graphics library for easy creation of figures
p <- ggplot(df, aes(x = Cultivar, y = Yield))
p <- p + geom_boxplot(aes(fill = CO2, linetype = TOS))
p <- p + facet_grid(Year ~ Irrigation)
p## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
This gives a first idea of what information could be hiding in the data. We are mostly interested on the effect of CO2 on yield in different cutlivars, at different times of sowing and different irrigation levels.
Analysis using aov
Analysis of Variance is a well-established procedure to check for effects of explanatory factors on a response variable.
It is important to ensure that the explanatory variables are factors! If the explanatory variables are integer or numeric, the following tests will be a regression analysis, not a strict ANOVA!
# to see the class of each vector in the data frame have a look at the structure of df
# make sure the explanatory variables to be used in the follwoing tests are indeed factors!!
str(df)## 'data.frame': 192 obs. of 20 variables:
## $ Plot.Order : int 1 2 3 3 4 5 5 6 6 7 ...
## $ CO2.TOS.Irr.Cultivar.DATABASE: Factor w/ 76 levels "aCO2 TOS1 Rain Drysdale",..: 17 36 14 33 36 27 36 24 33 33 ...
## $ TrialID : Factor w/ 2 levels "Horsham","Walpeup": 1 1 1 1 1 1 1 1 1 1 ...
## $ Year : Factor w/ 3 levels "2007","2008",..: 1 2 1 2 3 2 1 2 1 3 ...
## $ RingID : int 1 1 1 1 1 1 1 1 1 1 ...
## $ PlotID : Factor w/ 24 levels "A","B","C","D",..: 1 4 6 6 4 9 9 12 12 7 ...
## $ HalfRingID : Factor w/ 2 levels "E","W": 2 2 2 2 2 1 1 1 1 2 ...
## $ RingTrt : Factor w/ 13 levels "","aCO2 Rain",..: 3 12 3 12 12 12 3 12 3 12 ...
## $ CO2 : Factor w/ 2 levels "aCO2","eCO2": 1 1 1 1 1 1 1 1 1 1 ...
## $ Irrigation : Factor w/ 2 levels "Rain","Sup": 2 2 2 2 2 1 2 1 2 2 ...
## $ Nitrogen : Factor w/ 2 levels "N0","N+": 1 1 1 1 1 1 1 1 1 1 ...
## $ TOS : Factor w/ 2 levels "TOS1","TOS2": 1 2 1 2 2 2 2 2 2 2 ...
## $ PlotTrt : Factor w/ 11 levels "","Drysdale",..: 10 10 6 7 10 10 10 7 6 6 ...
## $ Sampling : Factor w/ 3 levels "","growth","Growth": 2 3 2 3 3 3 2 3 2 3 ...
## $ Cultivar : Factor w/ 2 levels "Janz","Yitpi": 2 2 1 1 2 2 2 1 1 1 ...
## $ Stage : Factor w/ 1 level "DC90": 1 1 1 1 1 1 1 1 1 1 ...
## $ Yield : num 327.9 163.7 299.1 63.3 140.9 ...
## $ YearRing : Factor w/ 48 levels "2007.1","2008.1",..: 1 2 1 2 3 2 1 2 1 3 ...
## $ Environment : Factor w/ 12 levels "Horsham.2007.Rain.TOS1",..: 4 11 4 11 12 8 10 8 10 12 ...
## $ Ord.Environment : Ord.factor w/ 12 levels "Horsham.2009.Rain.TOS2"<..: 12 4 12 4 3 2 6 2 6 3 ...# list if the parametes in the data frame are factors
sapply(df, is.factor)## Plot.Order CO2.TOS.Irr.Cultivar.DATABASE
## FALSE TRUE
## TrialID Year
## TRUE TRUE
## RingID PlotID
## FALSE TRUE
## HalfRingID RingTrt
## TRUE TRUE
## CO2 Irrigation
## TRUE TRUE
## Nitrogen TOS
## TRUE TRUE
## PlotTrt Sampling
## TRUE TRUE
## Cultivar Stage
## TRUE TRUE
## Yield YearRing
## FALSE TRUE
## Environment Ord.Environment
## TRUE TRUEMore than three factors are difficult to interpret simultaneously. Therefore, it’s adviced to limit the research question, data, and explanatory factor accordingly.
Subset the data
This subset keeps only data for the first time of sowing TOS, and we limit the data to the rainfed irrigation treatment. To create a smaller, more specific data set:
small <- df[df$TOS == "TOS1" &
df$Irrigation == "Rain", ]Simple AOV
Now, to test the effects of the three factors Year, CO2, Cultivar
# a simple AOV model
my.aov <- aov(Yield ~ Year * CO2 * Cultivar,
data = small)
summary(my.aov)## Df Sum Sq Mean Sq F value Pr(>F)
## Year 2 4058 2029 0.313 0.7334
## CO2 1 14848 14848 2.289 0.1392
## Cultivar 1 20364 20364 3.140 0.0851 .
## Year:CO2 2 3681 1840 0.284 0.7547
## Year:Cultivar 2 41793 20896 3.222 0.0519 .
## CO2:Cultivar 1 21 21 0.003 0.9546
## Year:CO2:Cultivar 2 2339 1169 0.180 0.8358
## Residuals 35 226988 6485
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1 observation deleted due to missingnessBefore we can interpret the model result, it is important to check if the variances are homogeneous across the tested groups.
library(car) # load package "Companion to Applied Regression"
leveneTest(Yield ~ Year * CO2 * Cultivar,
data = small)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 11 1.6483 0.1279
## 35In this case, the Levene test does not detect a statistically significant deviation of the homogeneity, so far, so good.
However, this anova model is in principle not statistically correct, unfortunately: It does not represent the experimental split-plot design!! However, the levene test for homogeneity of variances does not work with nested designs! See the test via the nlme package below for an other option.
In each year the experiment was set up as rings that contain plots with the cultivars (and irrigated / non-irrigated plots). Each ring had a fixed CO2 treatment. The model needs to be nested, the unnested model is wrong!
Care must be taken regarding the nesting levels. In this experiment, the plots are within a ring. Rings are identified by their RingID. However, due to the randomisation of the rings each year, rings can have a different treatment every year:
# create combinations of all possible cases for RingId, CO2 and Year
# only show unique entries, only show the first few entries
head(unique(interaction(small$RingID, small$CO2, small$Year)))## [1] 2.eCO2.2008 2.aCO2.2007 4.aCO2.2008 4.eCO2.2007 5.eCO2.2007 6.eCO2.2008
## 90 Levels: 2.aCO2.2007 3.aCO2.2007 4.aCO2.2007 5.aCO2.2007 ... 16.eCO2.2009The list above shows that Ring “2” had an “eCO2” treatment in 2008, but had “aCO2” in 2007! That means that a RingID might change it’s CO2 treatment beween the years. Ths will cause problems in the statistical test. Therefore the RingID must include the information on the year! An example is shown below.
Then the actual nested test can be done.
nested AOV
# a nested AOV model
# the nesting follows the actual experiment structure. What element is inside of another one? The largest 'container' is on the left of the error term.
my.aov <- aov(Yield ~ Year * CO2 * Cultivar + Error(Year/RingID),
data = small)
summary(my.aov)##
## Error: Year
## Df Sum Sq Mean Sq
## Year 2 4058 2029
##
## Error: Year:RingID
## Df Sum Sq Mean Sq
## CO2 1 66718 66718
## Cultivar 1 13219 13219
## Year:CO2 1 7068 7068
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## CO2 1 6998 6998 1.482 0.2324
## Cultivar 1 20116 20116 4.260 0.0472 *
## Year:CO2 2 546 273 0.058 0.9439
## Year:Cultivar 2 41879 20939 4.434 0.0200 *
## CO2:Cultivar 1 43 43 0.009 0.9246
## Year:CO2:Cultivar 2 2339 1170 0.248 0.7821
## Residuals 32 151108 4722
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1See the problem in the above model summary? CO2 is shown twice! Once for Year:RingID and again for Within! The CO2 effect shown for within is unexpected, as there is only one CO2 treatment ‘within’ one ring. how can there be a CO2 effect. This points to a flawed design!
As mentioned above, the rings can be named the same in different years!
# creating a new variable to identifies each Ring in each year
small$YearRing <- interaction(small$Year, small$RingID, drop = TRUE)
# now ring and year are coded as one nesting factor
my.aov <- aov(Yield ~ Year * CO2 * Cultivar + Error(YearRing),
data = small)
summary(my.aov)##
## Error: YearRing
## Df Sum Sq Mean Sq F value Pr(>F)
## Year 2 4058 2029 0.197 0.823
## CO2 1 14848 14848 1.440 0.247
## Cultivar 1 2011 2011 0.195 0.664
## Year:CO2 2 2409 1205 0.117 0.890
## Residuals 17 175230 10308
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Cultivar 1 22623 22623 7.881 0.01212 *
## Year:Cultivar 2 41584 20792 7.243 0.00531 **
## CO2:Cultivar 1 11 11 0.004 0.95082
## Year:CO2:Cultivar 2 2520 1260 0.439 0.65180
## Residuals 17 48798 2870
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# This is much better.
# In the case of the previous model with many more factors like time of sowing and irrigation, the design would be Cultivars, Irrigation, and Time of sowing are nested in rings and half rings each year. Half of each ring was irrigated, as indicated by the HalfringID.
my.aov2 <- aov(Yield ~ CO2 * Cultivar * Irrigation * TOS + Error(YearRing / HalfRingID),
data = df)
summary(my.aov2)##
## Error: YearRing
## Df Sum Sq Mean Sq F value Pr(>F)
## CO2 1 153449 153449 11.317 0.00168 **
## Cultivar 1 7770 7770 0.573 0.45340
## Irrigation 1 101686 101686 7.499 0.00909 **
## TOS 1 826553 826553 60.959 1.22e-09 ***
## CO2:Irrigation 1 13761 13761 1.015 0.31965
## CO2:TOS 1 24148 24148 1.781 0.18940
## Residuals 41 555928 13559
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: YearRing:HalfRingID
## Df Sum Sq Mean Sq F value Pr(>F)
## Cultivar 1 26216 26216 5.468 0.02432 *
## Irrigation 1 51747 51747 10.794 0.00209 **
## TOS 1 131221 131221 27.371 5.31e-06 ***
## CO2:Irrigation 1 29671 29671 6.189 0.01701 *
## CO2:TOS 1 221 221 0.046 0.83112
## Irrigation:TOS 1 9333 9333 1.947 0.17046
## CO2:Irrigation:TOS 1 1784 1784 0.372 0.54523
## Residuals 41 196561 4794
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## Cultivar 1 15480 15480 3.908 0.0512 .
## CO2:Cultivar 1 1258 1258 0.317 0.5746
## Cultivar:Irrigation 1 2827 2827 0.714 0.4006
## Cultivar:TOS 1 8415 8415 2.125 0.1485
## CO2:Cultivar:Irrigation 1 761 761 0.192 0.6622
## CO2:Cultivar:TOS 1 204 204 0.051 0.8211
## Cultivar:Irrigation:TOS 1 1143 1143 0.288 0.5926
## CO2:Cultivar:Irrigation:TOS 1 286 286 0.072 0.7889
## Residuals 87 344591 3961
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Be careful when interpreting nested models, keep the physical experimental plot layout in mind:
# CO2 effect are determined between rings (YearRing table). The irrigation and TOS difference between years and rings is not very relevant to our research question.
# Irrigation effect is estimated between Halfrings (YearRing:HalfRingID table). A cultivar effect is reported here as well.
# At the last level of nesting, "Within"" the halfring, there is the comparison between cultivars growing in the same halfring (= irrigation). No difference between cultivars within the halfring. But between the halfrings the cultivars respond differently. Therefore an interaction is reported in the YearRing:HalfRingID table.Linear mixed effect models
Linear mixed effect model make it easier to work with nested statistical models. This is an example using the package “nlme”.
library(nlme)
# complex model
my.lme <- lme(Yield ~ Year * CO2 * Irrigation * Cultivar,
random = ~ 1 | Year/YearRing/Irrigation,
data = df[df$TOS == "TOS1", ],
na.action = na.exclude)
# This model complains about missing values!
# na.action allows to specify what to do with missing data for the fitted object
# This still does not help much, as there is no automatic gap-filling, just the length of residuals and predictions get padded to the correct length of the data. See help page for ?na.omit
# However, excluding missing values causes problems with the diagnostic figures later. na.omit is a solution for that. Again, "good" balanced datas is important for any statistical analysis!
# Detaild result of the test:
summary(my.lme)## Warning in pt(-abs(tVal), fDF): NaNs produced## Linear mixed-effects model fit by REML
## Data: df[df$TOS == "TOS1", ]
## AIC BIC logLik
## 906.7741 970.1291 -425.387
##
## Random effects:
## Formula: ~1 | Year
## (Intercept)
## StdDev: 30.41083
##
## Formula: ~1 | YearRing %in% Year
## (Intercept)
## StdDev: 28.81868
##
## Formula: ~1 | Irrigation %in% YearRing %in% Year
## (Intercept) Residual
## StdDev: 36.93418 64.17431
##
## Fixed effects: Yield ~ Year * CO2 * Irrigation * Cultivar
## Value Std.Error DF
## (Intercept) 263.62928 50.03066 35
## Year2008 -33.26713 70.75403 0
## Year2009 14.94443 70.75403 0
## CO2eCO2 12.46106 56.18270 26
## IrrigationSup 85.08567 56.18270 10
## CultivarYitpi 40.37701 50.21628 35
## Year2008:CO2eCO2 25.12948 79.45434 26
## Year2009:CO2eCO2 36.38509 79.45434 26
## Year2008:IrrigationSup -82.31114 76.79670 10
## Year2009:IrrigationSup -83.36772 76.79670 10
## CO2eCO2:IrrigationSup 73.59813 79.45434 10
## Year2008:CultivarYitpi 58.14887 67.68194 35
## Year2009:CultivarYitpi -51.83214 67.68194 35
## CO2eCO2:CultivarYitpi -3.96735 67.68194 35
## IrrigationSup:CultivarYitpi -31.61533 67.68194 35
## Year2008:CO2eCO2:IrrigationSup 71.03680 108.60694 10
## Year2009:CO2eCO2:IrrigationSup -11.48916 108.60694 10
## Year2008:CO2eCO2:CultivarYitpi 41.40781 93.26943 35
## Year2009:CO2eCO2:CultivarYitpi -26.95893 93.26943 35
## Year2008:IrrigationSup:CultivarYitpi 57.71166 93.26943 35
## Year2009:IrrigationSup:CultivarYitpi 8.07793 93.26943 35
## CO2eCO2:IrrigationSup:CultivarYitpi -34.19464 93.26943 35
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -49.24293 130.13789 35
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi 65.98203 130.13789 35
## t-value p-value
## (Intercept) 5.269355 0.0000
## Year2008 -0.470180 NaN
## Year2009 0.211217 NaN
## CO2eCO2 0.221795 0.8262
## IrrigationSup 1.514446 0.1609
## CultivarYitpi 0.804062 0.4268
## Year2008:CO2eCO2 0.316276 0.7543
## Year2009:CO2eCO2 0.457937 0.6508
## Year2008:IrrigationSup -1.071806 0.3090
## Year2009:IrrigationSup -1.085564 0.3031
## CO2eCO2:IrrigationSup 0.926295 0.3761
## Year2008:CultivarYitpi 0.859149 0.3961
## Year2009:CultivarYitpi -0.765819 0.4489
## CO2eCO2:CultivarYitpi -0.058618 0.9536
## IrrigationSup:CultivarYitpi -0.467116 0.6433
## Year2008:CO2eCO2:IrrigationSup 0.654072 0.5278
## Year2009:CO2eCO2:IrrigationSup -0.105787 0.9178
## Year2008:CO2eCO2:CultivarYitpi 0.443959 0.6598
## Year2009:CO2eCO2:CultivarYitpi -0.289044 0.7743
## Year2008:IrrigationSup:CultivarYitpi 0.618763 0.5401
## Year2009:IrrigationSup:CultivarYitpi 0.086609 0.9315
## CO2eCO2:IrrigationSup:CultivarYitpi -0.366622 0.7161
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -0.378390 0.7074
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi 0.507016 0.6153
## Correlation:
## (Intr) Yr2008 Yr2009 CO2eCO2
## Year2008 -0.707
## Year2009 -0.707 0.500
## CO2eCO2 -0.561 0.397 0.397
## IrrigationSup -0.561 0.397 0.397 0.500
## CultivarYitpi -0.410 0.290 0.290 0.365
## Year2008:CO2eCO2 0.397 -0.561 -0.281 -0.707
## Year2009:CO2eCO2 0.397 -0.281 -0.561 -0.707
## Year2008:IrrigationSup 0.411 -0.543 -0.290 -0.366
## Year2009:IrrigationSup 0.411 -0.290 -0.543 -0.366
## CO2eCO2:IrrigationSup 0.397 -0.281 -0.281 -0.707
## Year2008:CultivarYitpi 0.304 -0.430 -0.215 -0.271
## Year2009:CultivarYitpi 0.304 -0.215 -0.430 -0.271
## CO2eCO2:CultivarYitpi 0.304 -0.215 -0.215 -0.542
## IrrigationSup:CultivarYitpi 0.304 -0.215 -0.215 -0.271
## Year2008:CO2eCO2:IrrigationSup -0.290 0.384 0.205 0.517
## Year2009:CO2eCO2:IrrigationSup -0.290 0.205 0.384 0.517
## Year2008:CO2eCO2:CultivarYitpi -0.221 0.312 0.156 0.393
## Year2009:CO2eCO2:CultivarYitpi -0.221 0.156 0.312 0.393
## Year2008:IrrigationSup:CultivarYitpi -0.221 0.312 0.156 0.196
## Year2009:IrrigationSup:CultivarYitpi -0.221 0.156 0.312 0.196
## CO2eCO2:IrrigationSup:CultivarYitpi -0.221 0.156 0.156 0.393
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi 0.158 -0.224 -0.112 -0.282
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi 0.158 -0.112 -0.224 -0.282
## IrrgtS CltvrY Yr2008:CO2CO2
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi 0.365
## Year2008:CO2eCO2 -0.354 -0.258
## Year2009:CO2eCO2 -0.354 -0.258 0.500
## Year2008:IrrigationSup -0.732 -0.267 0.483
## Year2009:IrrigationSup -0.732 -0.267 0.259
## CO2eCO2:IrrigationSup -0.707 -0.258 0.500
## Year2008:CultivarYitpi -0.271 -0.742 0.383
## Year2009:CultivarYitpi -0.271 -0.742 0.191
## CO2eCO2:CultivarYitpi -0.271 -0.742 0.383
## IrrigationSup:CultivarYitpi -0.542 -0.742 0.191
## Year2008:CO2eCO2:IrrigationSup 0.517 0.189 -0.683
## Year2009:CO2eCO2:IrrigationSup 0.517 0.189 -0.366
## Year2008:CO2eCO2:CultivarYitpi 0.196 0.538 -0.556
## Year2009:CO2eCO2:CultivarYitpi 0.196 0.538 -0.278
## Year2008:IrrigationSup:CultivarYitpi 0.393 0.538 -0.278
## Year2009:IrrigationSup:CultivarYitpi 0.393 0.538 -0.139
## CO2eCO2:IrrigationSup:CultivarYitpi 0.393 0.538 -0.278
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -0.282 -0.386 0.398
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi -0.282 -0.386 0.199
## Yr2009:CO2CO2 Yr2008:IS
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup 0.259
## Year2009:IrrigationSup 0.483 0.535
## CO2eCO2:IrrigationSup 0.500 0.517
## Year2008:CultivarYitpi 0.191 0.396
## Year2009:CultivarYitpi 0.383 0.198
## CO2eCO2:CultivarYitpi 0.383 0.198
## IrrigationSup:CultivarYitpi 0.191 0.396
## Year2008:CO2eCO2:IrrigationSup -0.366 -0.707
## Year2009:CO2eCO2:IrrigationSup -0.683 -0.378
## Year2008:CO2eCO2:CultivarYitpi -0.278 -0.287
## Year2009:CO2eCO2:CultivarYitpi -0.556 -0.144
## Year2008:IrrigationSup:CultivarYitpi -0.139 -0.575
## Year2009:IrrigationSup:CultivarYitpi -0.278 -0.287
## CO2eCO2:IrrigationSup:CultivarYitpi -0.278 -0.287
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi 0.199 0.412
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi 0.398 0.206
## Yr2009:IS CO2CO2:IrS Y2008:CY
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup
## Year2009:IrrigationSup
## CO2eCO2:IrrigationSup 0.517
## Year2008:CultivarYitpi 0.198 0.191
## Year2009:CultivarYitpi 0.396 0.191 0.550
## CO2eCO2:CultivarYitpi 0.198 0.383 0.550
## IrrigationSup:CultivarYitpi 0.396 0.383 0.550
## Year2008:CO2eCO2:IrrigationSup -0.378 -0.732 -0.280
## Year2009:CO2eCO2:IrrigationSup -0.707 -0.732 -0.140
## Year2008:CO2eCO2:CultivarYitpi -0.144 -0.278 -0.726
## Year2009:CO2eCO2:CultivarYitpi -0.287 -0.278 -0.399
## Year2008:IrrigationSup:CultivarYitpi -0.287 -0.278 -0.726
## Year2009:IrrigationSup:CultivarYitpi -0.575 -0.278 -0.399
## CO2eCO2:IrrigationSup:CultivarYitpi -0.287 -0.556 -0.399
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi 0.206 0.398 0.520
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi 0.412 0.398 0.286
## Y2009:CY CO2CO2:C IrS:CY
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup
## Year2009:IrrigationSup
## CO2eCO2:IrrigationSup
## Year2008:CultivarYitpi
## Year2009:CultivarYitpi
## CO2eCO2:CultivarYitpi 0.550
## IrrigationSup:CultivarYitpi 0.550 0.550
## Year2008:CO2eCO2:IrrigationSup -0.140 -0.280 -0.280
## Year2009:CO2eCO2:IrrigationSup -0.280 -0.280 -0.280
## Year2008:CO2eCO2:CultivarYitpi -0.399 -0.726 -0.399
## Year2009:CO2eCO2:CultivarYitpi -0.726 -0.726 -0.399
## Year2008:IrrigationSup:CultivarYitpi -0.399 -0.399 -0.726
## Year2009:IrrigationSup:CultivarYitpi -0.726 -0.399 -0.726
## CO2eCO2:IrrigationSup:CultivarYitpi -0.399 -0.726 -0.726
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi 0.286 0.520 0.520
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi 0.520 0.520 0.520
## Yr2008:CO2CO2:IS
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup
## Year2009:IrrigationSup
## CO2eCO2:IrrigationSup
## Year2008:CultivarYitpi
## Year2009:CultivarYitpi
## CO2eCO2:CultivarYitpi
## IrrigationSup:CultivarYitpi
## Year2008:CO2eCO2:IrrigationSup
## Year2009:CO2eCO2:IrrigationSup 0.535
## Year2008:CO2eCO2:CultivarYitpi 0.407
## Year2009:CO2eCO2:CultivarYitpi 0.203
## Year2008:IrrigationSup:CultivarYitpi 0.407
## Year2009:IrrigationSup:CultivarYitpi 0.203
## CO2eCO2:IrrigationSup:CultivarYitpi 0.407
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -0.583
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi -0.291
## Yr2009:CO2CO2:IS
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup
## Year2009:IrrigationSup
## CO2eCO2:IrrigationSup
## Year2008:CultivarYitpi
## Year2009:CultivarYitpi
## CO2eCO2:CultivarYitpi
## IrrigationSup:CultivarYitpi
## Year2008:CO2eCO2:IrrigationSup
## Year2009:CO2eCO2:IrrigationSup
## Year2008:CO2eCO2:CultivarYitpi 0.203
## Year2009:CO2eCO2:CultivarYitpi 0.407
## Year2008:IrrigationSup:CultivarYitpi 0.203
## Year2009:IrrigationSup:CultivarYitpi 0.407
## CO2eCO2:IrrigationSup:CultivarYitpi 0.407
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -0.291
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi -0.583
## Y2008:CO2CO2:C Y2009:CO2CO2:C
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup
## Year2009:IrrigationSup
## CO2eCO2:IrrigationSup
## Year2008:CultivarYitpi
## Year2009:CultivarYitpi
## CO2eCO2:CultivarYitpi
## IrrigationSup:CultivarYitpi
## Year2008:CO2eCO2:IrrigationSup
## Year2009:CO2eCO2:IrrigationSup
## Year2008:CO2eCO2:CultivarYitpi
## Year2009:CO2eCO2:CultivarYitpi 0.527
## Year2008:IrrigationSup:CultivarYitpi 0.527 0.290
## Year2009:IrrigationSup:CultivarYitpi 0.290 0.527
## CO2eCO2:IrrigationSup:CultivarYitpi 0.527 0.527
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -0.717 -0.377
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi -0.377 -0.717
## Y2008:IS: Y2009:IS:
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup
## Year2009:IrrigationSup
## CO2eCO2:IrrigationSup
## Year2008:CultivarYitpi
## Year2009:CultivarYitpi
## CO2eCO2:CultivarYitpi
## IrrigationSup:CultivarYitpi
## Year2008:CO2eCO2:IrrigationSup
## Year2009:CO2eCO2:IrrigationSup
## Year2008:CO2eCO2:CultivarYitpi
## Year2009:CO2eCO2:CultivarYitpi
## Year2008:IrrigationSup:CultivarYitpi
## Year2009:IrrigationSup:CultivarYitpi 0.527
## CO2eCO2:IrrigationSup:CultivarYitpi 0.527 0.527
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -0.717 -0.377
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi -0.377 -0.717
## CO2CO2:IS: Y2008:CO2CO2:IS:
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup
## Year2009:IrrigationSup
## CO2eCO2:IrrigationSup
## Year2008:CultivarYitpi
## Year2009:CultivarYitpi
## CO2eCO2:CultivarYitpi
## IrrigationSup:CultivarYitpi
## Year2008:CO2eCO2:IrrigationSup
## Year2009:CO2eCO2:IrrigationSup
## Year2008:CO2eCO2:CultivarYitpi
## Year2009:CO2eCO2:CultivarYitpi
## Year2008:IrrigationSup:CultivarYitpi
## Year2009:IrrigationSup:CultivarYitpi
## CO2eCO2:IrrigationSup:CultivarYitpi
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -0.717
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi -0.717 0.514
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.80337239 -0.40297439 -0.04572708 0.45929044 1.85779558
##
## Number of Observations: 95
## Number of Groups:
## Year YearRing %in% Year
## 3 32
## Irrigation %in% YearRing %in% Year
## 48# anova table of the test
anova(my.lme)## Warning in pf(Fval[i], nDF[i], dDF[i]): NaNs produced## numDF denDF F-value p-value
## (Intercept) 1 35 263.04571 <.0001
## Year 2 0 0.34644 NaN
## CO2 1 26 12.65555 0.0015
## Irrigation 1 10 12.90896 0.0049
## Cultivar 1 35 5.37444 0.0264
## Year:CO2 2 26 1.00370 0.3803
## Year:Irrigation 2 10 1.26983 0.3225
## CO2:Irrigation 1 10 5.29164 0.0442
## Year:Cultivar 2 35 9.97569 0.0004
## CO2:Cultivar 1 35 0.24302 0.6251
## Irrigation:Cultivar 1 35 0.82133 0.3710
## Year:CO2:Irrigation 2 10 0.14463 0.8671
## Year:CO2:Cultivar 2 35 0.03304 0.9675
## Year:Irrigation:Cultivar 2 35 0.23334 0.7931
## CO2:Irrigation:Cultivar 1 35 0.29089 0.5931
## Year:CO2:Irrigation:Cultivar 2 35 0.40571 0.6696# test normality for each nesting group
# level indicates nesting level, there are four nesting levels in the model
# check for normally distributed data via a qq plot:
# e.g. for level YearRing
qqnorm(my.lme, ~ranef(., level = 2))
# This looks rather straight, normally distributed data can be assumed here.
# Another benefit of lme is, that it provides methods to check the homogeneity of residuals for sub-groups visually by using the model itself.
# lme, similar to aov expects that the variance of the residuals are homogeneous across groups. This can be checked comparing the figures below.
plot(my.lme, resid(.) ~ fitted(.) | YearRing, id = 0.05)
# the error message in the figure is due to the missing data.
# this does not look very homogeneous, especially as there only few data per subgroup
plot(my.lme, resid(.) ~ fitted(.) | CO2, id = 0.05)
# Again, the error message in the figure is due to the missing data. As the fitted model data is padded with missing values, the error message is created. This can be omitted by changing na.action = na.omit during model creation.
plot(my.lme, resid(.) ~ fitted(.) | Irrigation, id = 0.05)
plot(my.lme, resid(.) ~ fitted(.) | Cultivar, id = 0.05)
# The figure for cultivar looks rather sketchy, not very homogeneous.
# The model can take the variable variance across subgroups into account. In this case via the option VarIdent. The existing model is updated with aditional information on the underlying variance structure.
my.lme2 <- update(my.lme, weights = varIdent(form = ~ 1 | Cultivar) )
summary(my.lme2)## Warning in pt(-abs(tVal), fDF): NaNs produced## Linear mixed-effects model fit by REML
## Data: df[df$TOS == "TOS1", ]
## AIC BIC logLik
## 908.6012 974.2189 -425.3006
##
## Random effects:
## Formula: ~1 | Year
## (Intercept)
## StdDev: 28.81416
##
## Formula: ~1 | YearRing %in% Year
## (Intercept)
## StdDev: 29.75889
##
## Formula: ~1 | Irrigation %in% YearRing %in% Year
## (Intercept) Residual
## StdDev: 36.27708 60.80472
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Cultivar
## Parameter estimates:
## Yitpi Janz
## 1.00000 1.10708
## Fixed effects: Yield ~ Year * CO2 * Irrigation * Cultivar
## Value Std.Error DF
## (Intercept) 263.62928 50.13493 35
## Year2008 -33.26713 70.90150 0
## Year2009 14.94443 70.90150 0
## CO2eCO2 12.46106 58.02164 26
## IrrigationSup 85.08567 58.02164 10
## CultivarYitpi 40.05264 49.88719 35
## Year2008:CO2eCO2 25.12948 82.05500 26
## Year2009:CO2eCO2 36.38509 82.05500 26
## Year2008:IrrigationSup -82.31114 79.31095 10
## Year2009:IrrigationSup -83.36772 79.31095 10
## CO2eCO2:IrrigationSup 73.59813 82.05500 10
## Year2008:CultivarYitpi 58.47324 67.42319 35
## Year2009:CultivarYitpi -51.50777 67.42319 35
## CO2eCO2:CultivarYitpi -3.64299 67.42319 35
## IrrigationSup:CultivarYitpi -31.29096 67.42319 35
## Year2008:CO2eCO2:IrrigationSup 71.03680 112.16262 10
## Year2009:CO2eCO2:IrrigationSup -11.48916 112.16262 10
## Year2008:CO2eCO2:CultivarYitpi 41.08345 93.06017 35
## Year2009:CO2eCO2:CultivarYitpi -27.28330 93.06017 35
## Year2008:IrrigationSup:CultivarYitpi 57.38729 93.06017 35
## Year2009:IrrigationSup:CultivarYitpi 7.75357 93.06017 35
## CO2eCO2:IrrigationSup:CultivarYitpi -34.51901 93.06017 35
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -48.91857 129.95697 35
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi 66.30640 129.95697 35
## t-value p-value
## (Intercept) 5.258395 0.0000
## Year2008 -0.469202 NaN
## Year2009 0.210777 NaN
## CO2eCO2 0.214766 0.8316
## IrrigationSup 1.466447 0.1733
## CultivarYitpi 0.802864 0.4275
## Year2008:CO2eCO2 0.306252 0.7619
## Year2009:CO2eCO2 0.443423 0.6611
## Year2008:IrrigationSup -1.037828 0.3238
## Year2009:IrrigationSup -1.051150 0.3179
## CO2eCO2:IrrigationSup 0.896937 0.3908
## Year2008:CultivarYitpi 0.867257 0.3917
## Year2009:CultivarYitpi -0.763947 0.4500
## CO2eCO2:CultivarYitpi -0.054032 0.9572
## IrrigationSup:CultivarYitpi -0.464098 0.6455
## Year2008:CO2eCO2:IrrigationSup 0.633338 0.5407
## Year2009:CO2eCO2:IrrigationSup -0.102433 0.9204
## Year2008:CO2eCO2:CultivarYitpi 0.441472 0.6616
## Year2009:CO2eCO2:CultivarYitpi -0.293179 0.7711
## Year2008:IrrigationSup:CultivarYitpi 0.616669 0.5414
## Year2009:IrrigationSup:CultivarYitpi 0.083318 0.9341
## CO2eCO2:IrrigationSup:CultivarYitpi -0.370932 0.7129
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -0.376421 0.7089
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi 0.510218 0.6131
## Correlation:
## (Intr) Yr2008 Yr2009 CO2eCO2
## Year2008 -0.707
## Year2009 -0.707 0.500
## CO2eCO2 -0.579 0.409 0.409
## IrrigationSup -0.579 0.409 0.409 0.500
## CultivarYitpi -0.453 0.320 0.320 0.391
## Year2008:CO2eCO2 0.409 -0.579 -0.289 -0.707
## Year2009:CO2eCO2 0.409 -0.289 -0.579 -0.707
## Year2008:IrrigationSup 0.423 -0.559 -0.299 -0.366
## Year2009:IrrigationSup 0.423 -0.299 -0.559 -0.366
## CO2eCO2:IrrigationSup 0.409 -0.289 -0.289 -0.707
## Year2008:CultivarYitpi 0.335 -0.474 -0.237 -0.290
## Year2009:CultivarYitpi 0.335 -0.237 -0.474 -0.290
## CO2eCO2:CultivarYitpi 0.335 -0.237 -0.237 -0.579
## IrrigationSup:CultivarYitpi 0.335 -0.237 -0.237 -0.290
## Year2008:CO2eCO2:IrrigationSup -0.299 0.395 0.212 0.517
## Year2009:CO2eCO2:IrrigationSup -0.299 0.212 0.395 0.517
## Year2008:CO2eCO2:CultivarYitpi -0.243 0.343 0.172 0.420
## Year2009:CO2eCO2:CultivarYitpi -0.243 0.172 0.343 0.420
## Year2008:IrrigationSup:CultivarYitpi -0.243 0.343 0.172 0.210
## Year2009:IrrigationSup:CultivarYitpi -0.243 0.172 0.343 0.210
## CO2eCO2:IrrigationSup:CultivarYitpi -0.243 0.172 0.172 0.420
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi 0.174 -0.246 -0.123 -0.300
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi 0.174 -0.123 -0.246 -0.300
## IrrgtS CltvrY Yr2008:CO2CO2
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi 0.391
## Year2008:CO2eCO2 -0.354 -0.277
## Year2009:CO2eCO2 -0.354 -0.277 0.500
## Year2008:IrrigationSup -0.732 -0.286 0.483
## Year2009:IrrigationSup -0.732 -0.286 0.259
## CO2eCO2:IrrigationSup -0.707 -0.277 0.500
## Year2008:CultivarYitpi -0.290 -0.740 0.410
## Year2009:CultivarYitpi -0.290 -0.740 0.205
## CO2eCO2:CultivarYitpi -0.290 -0.740 0.410
## IrrigationSup:CultivarYitpi -0.579 -0.740 0.205
## Year2008:CO2eCO2:IrrigationSup 0.517 0.202 -0.683
## Year2009:CO2eCO2:IrrigationSup 0.517 0.202 -0.366
## Year2008:CO2eCO2:CultivarYitpi 0.210 0.536 -0.593
## Year2009:CO2eCO2:CultivarYitpi 0.210 0.536 -0.297
## Year2008:IrrigationSup:CultivarYitpi 0.420 0.536 -0.297
## Year2009:IrrigationSup:CultivarYitpi 0.420 0.536 -0.148
## CO2eCO2:IrrigationSup:CultivarYitpi 0.420 0.536 -0.297
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -0.300 -0.384 0.425
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi -0.300 -0.384 0.212
## Yr2009:CO2CO2 Yr2008:IS
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup 0.259
## Year2009:IrrigationSup 0.483 0.535
## CO2eCO2:IrrigationSup 0.500 0.517
## Year2008:CultivarYitpi 0.205 0.424
## Year2009:CultivarYitpi 0.410 0.212
## CO2eCO2:CultivarYitpi 0.410 0.212
## IrrigationSup:CultivarYitpi 0.205 0.424
## Year2008:CO2eCO2:IrrigationSup -0.366 -0.707
## Year2009:CO2eCO2:IrrigationSup -0.683 -0.378
## Year2008:CO2eCO2:CultivarYitpi -0.297 -0.307
## Year2009:CO2eCO2:CultivarYitpi -0.593 -0.153
## Year2008:IrrigationSup:CultivarYitpi -0.148 -0.614
## Year2009:IrrigationSup:CultivarYitpi -0.297 -0.307
## CO2eCO2:IrrigationSup:CultivarYitpi -0.297 -0.307
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi 0.212 0.440
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi 0.425 0.220
## Yr2009:IS CO2CO2:IrS Y2008:CY
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup
## Year2009:IrrigationSup
## CO2eCO2:IrrigationSup 0.517
## Year2008:CultivarYitpi 0.212 0.205
## Year2009:CultivarYitpi 0.424 0.205 0.547
## CO2eCO2:CultivarYitpi 0.212 0.410 0.547
## IrrigationSup:CultivarYitpi 0.424 0.410 0.547
## Year2008:CO2eCO2:IrrigationSup -0.378 -0.732 -0.300
## Year2009:CO2eCO2:IrrigationSup -0.707 -0.732 -0.150
## Year2008:CO2eCO2:CultivarYitpi -0.153 -0.297 -0.725
## Year2009:CO2eCO2:CultivarYitpi -0.307 -0.297 -0.397
## Year2008:IrrigationSup:CultivarYitpi -0.307 -0.297 -0.725
## Year2009:IrrigationSup:CultivarYitpi -0.614 -0.297 -0.397
## CO2eCO2:IrrigationSup:CultivarYitpi -0.307 -0.593 -0.397
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi 0.220 0.425 0.519
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi 0.440 0.425 0.284
## Y2009:CY CO2CO2:C IrS:CY
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup
## Year2009:IrrigationSup
## CO2eCO2:IrrigationSup
## Year2008:CultivarYitpi
## Year2009:CultivarYitpi
## CO2eCO2:CultivarYitpi 0.547
## IrrigationSup:CultivarYitpi 0.547 0.547
## Year2008:CO2eCO2:IrrigationSup -0.150 -0.300 -0.300
## Year2009:CO2eCO2:IrrigationSup -0.300 -0.300 -0.300
## Year2008:CO2eCO2:CultivarYitpi -0.397 -0.725 -0.397
## Year2009:CO2eCO2:CultivarYitpi -0.725 -0.725 -0.397
## Year2008:IrrigationSup:CultivarYitpi -0.397 -0.397 -0.725
## Year2009:IrrigationSup:CultivarYitpi -0.725 -0.397 -0.725
## CO2eCO2:IrrigationSup:CultivarYitpi -0.397 -0.725 -0.725
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi 0.284 0.519 0.519
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi 0.519 0.519 0.519
## Yr2008:CO2CO2:IS
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup
## Year2009:IrrigationSup
## CO2eCO2:IrrigationSup
## Year2008:CultivarYitpi
## Year2009:CultivarYitpi
## CO2eCO2:CultivarYitpi
## IrrigationSup:CultivarYitpi
## Year2008:CO2eCO2:IrrigationSup
## Year2009:CO2eCO2:IrrigationSup 0.535
## Year2008:CO2eCO2:CultivarYitpi 0.434
## Year2009:CO2eCO2:CultivarYitpi 0.217
## Year2008:IrrigationSup:CultivarYitpi 0.434
## Year2009:IrrigationSup:CultivarYitpi 0.217
## CO2eCO2:IrrigationSup:CultivarYitpi 0.434
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -0.622
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi -0.311
## Yr2009:CO2CO2:IS
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup
## Year2009:IrrigationSup
## CO2eCO2:IrrigationSup
## Year2008:CultivarYitpi
## Year2009:CultivarYitpi
## CO2eCO2:CultivarYitpi
## IrrigationSup:CultivarYitpi
## Year2008:CO2eCO2:IrrigationSup
## Year2009:CO2eCO2:IrrigationSup
## Year2008:CO2eCO2:CultivarYitpi 0.217
## Year2009:CO2eCO2:CultivarYitpi 0.434
## Year2008:IrrigationSup:CultivarYitpi 0.217
## Year2009:IrrigationSup:CultivarYitpi 0.434
## CO2eCO2:IrrigationSup:CultivarYitpi 0.434
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -0.311
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi -0.622
## Y2008:CO2CO2:C Y2009:CO2CO2:C
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup
## Year2009:IrrigationSup
## CO2eCO2:IrrigationSup
## Year2008:CultivarYitpi
## Year2009:CultivarYitpi
## CO2eCO2:CultivarYitpi
## IrrigationSup:CultivarYitpi
## Year2008:CO2eCO2:IrrigationSup
## Year2009:CO2eCO2:IrrigationSup
## Year2008:CO2eCO2:CultivarYitpi
## Year2009:CO2eCO2:CultivarYitpi 0.525
## Year2008:IrrigationSup:CultivarYitpi 0.525 0.287
## Year2009:IrrigationSup:CultivarYitpi 0.287 0.525
## CO2eCO2:IrrigationSup:CultivarYitpi 0.525 0.525
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -0.716 -0.376
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi -0.376 -0.716
## Y2008:IS: Y2009:IS:
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup
## Year2009:IrrigationSup
## CO2eCO2:IrrigationSup
## Year2008:CultivarYitpi
## Year2009:CultivarYitpi
## CO2eCO2:CultivarYitpi
## IrrigationSup:CultivarYitpi
## Year2008:CO2eCO2:IrrigationSup
## Year2009:CO2eCO2:IrrigationSup
## Year2008:CO2eCO2:CultivarYitpi
## Year2009:CO2eCO2:CultivarYitpi
## Year2008:IrrigationSup:CultivarYitpi
## Year2009:IrrigationSup:CultivarYitpi 0.525
## CO2eCO2:IrrigationSup:CultivarYitpi 0.525 0.525
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -0.716 -0.376
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi -0.376 -0.716
## CO2CO2:IS: Y2008:CO2CO2:IS:
## Year2008
## Year2009
## CO2eCO2
## IrrigationSup
## CultivarYitpi
## Year2008:CO2eCO2
## Year2009:CO2eCO2
## Year2008:IrrigationSup
## Year2009:IrrigationSup
## CO2eCO2:IrrigationSup
## Year2008:CultivarYitpi
## Year2009:CultivarYitpi
## CO2eCO2:CultivarYitpi
## IrrigationSup:CultivarYitpi
## Year2008:CO2eCO2:IrrigationSup
## Year2009:CO2eCO2:IrrigationSup
## Year2008:CO2eCO2:CultivarYitpi
## Year2009:CO2eCO2:CultivarYitpi
## Year2008:IrrigationSup:CultivarYitpi
## Year2009:IrrigationSup:CultivarYitpi
## CO2eCO2:IrrigationSup:CultivarYitpi
## Year2008:CO2eCO2:IrrigationSup:CultivarYitpi -0.716
## Year2009:CO2eCO2:IrrigationSup:CultivarYitpi -0.716 0.513
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.78849627 -0.41675242 -0.05124116 0.44837204 1.85933371
##
## Number of Observations: 95
## Number of Groups:
## Year YearRing %in% Year
## 3 32
## Irrigation %in% YearRing %in% Year
## 48# to account for potential variance across the other groups, additional models can be designed.
my.lme3 <- update(my.lme, weights = varIdent(form = ~ 1 | Year) )
my.lme4 <- update(my.lme, weights = varIdent(form = ~ 1 | Irrigation) )
# compare the models
anova(my.lme, my.lme2, my.lme3, my.lme4)## Model df AIC BIC logLik Test L.Ratio p-value
## my.lme 1 28 906.7741 970.1291 -425.3870
## my.lme2 2 29 908.6012 974.2189 -425.3006 1 vs 2 0.1728385 0.6776
## my.lme3 3 30 908.9004 976.7808 -424.4502 2 vs 3 1.7008405 0.1922
## my.lme4 4 29 908.0127 973.6304 -425.0063 3 vs 4 1.1122820 0.2916Comparing the four models reveals no clear difference between them. In this case, the model with the lowest Akaike criterion (AIC) is prefferred. I.e. the original model:
# therefore, the final result (using the original model) is
anova(my.lme)## Warning in pf(Fval[i], nDF[i], dDF[i]): NaNs produced## numDF denDF F-value p-value
## (Intercept) 1 35 263.04571 <.0001
## Year 2 0 0.34644 NaN
## CO2 1 26 12.65555 0.0015
## Irrigation 1 10 12.90896 0.0049
## Cultivar 1 35 5.37444 0.0264
## Year:CO2 2 26 1.00370 0.3803
## Year:Irrigation 2 10 1.26983 0.3225
## CO2:Irrigation 1 10 5.29164 0.0442
## Year:Cultivar 2 35 9.97569 0.0004
## CO2:Cultivar 1 35 0.24302 0.6251
## Irrigation:Cultivar 1 35 0.82133 0.3710
## Year:CO2:Irrigation 2 10 0.14463 0.8671
## Year:CO2:Cultivar 2 35 0.03304 0.9675
## Year:Irrigation:Cultivar 2 35 0.23334 0.7931
## CO2:Irrigation:Cultivar 1 35 0.29089 0.5931
## Year:CO2:Irrigation:Cultivar 2 35 0.40571 0.6696Second example
Repeated measure lme
Soil moisure data is collected weekly in eight rings for four soil cores within the rings. Each profile probe measurement consists of six readings, representing six depths within the plot. Half of the soil cores have wicks installed. There are two soiltypes in this experiment, but this example only uses one of them for simplicity. The random effects in this case is Date / Ring / Core / Depth. At each date samples are taken in each ring. In each ring, samples are taken in each core. In each core six measurements of soil moisture take place. In this experiment we want to know if the wicks change the soil moisture in the cores. This is a one-way test with four split levels. Date is used as the first split level, as all measurement take place within each date during each repeat. The example below uses Date, Ring, CoreID and Depth as grouping elements. However, we could assume that Date is a random effect (the weather on each date might affect the soil moisture). In this case, Date is put in front of the “|” in the random term. See ?lme for details on the structure of the random term. In this example though, date is kept as a grouping factor, as it is clear that the soil moisture increases currently with steady rainfall and therefore each date will have a different absolute soil moisture. The Date is therefore not a random element, but has, in context of soil moisture, a distinct order that is used in the model below.
my.lme <- lme(Percent_Vol ~ has.wick,
random = ~ 1 | Date / Ring / CoreID / Depth,
data = sm[sm$Soil == "HOR", ], # only testing one soiltype
na.action = na.omit) # there are missing values in the data.
# in contrast to the previous example, na.omit is used. This means that residuals and predictions will not have missing data and there will be no further error messages about it. The problem of a missing data still exists, though!.
# summary(my.lme) # gives additional diagnostic information
# the anova table
summary(my.lme)## Linear mixed-effects model fit by REML
## Data: sm[sm$Soil == "HOR", ]
## AIC BIC logLik
## 6042.316 6074.143 -3014.158
##
## Random effects:
## Formula: ~1 | Date
## (Intercept)
## StdDev: 6.243158
##
## Formula: ~1 | Ring %in% Date
## (Intercept)
## StdDev: 0.0005711042
##
## Formula: ~1 | CoreID %in% Ring %in% Date
## (Intercept)
## StdDev: 0.0006928839
##
## Formula: ~1 | Depth %in% CoreID %in% Ring %in% Date
## (Intercept) Residual
## StdDev: 17.87937 0.07122964
##
## Fixed effects: Percent_Vol ~ has.wick
## Value Std.Error DF t-value p-value
## (Intercept) 32.97151 2.288324 555 14.408585 0.0000
## has.wickWick 2.91357 1.352708 62 2.153876 0.0351
## Correlation:
## (Intr)
## has.wickWick -0.293
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -0.0079255957 -0.0031808479 -0.0004136288 0.0030574552 0.0140279655
##
## Number of Observations: 699
## Number of Groups:
## Date
## 9
## Ring %in% Date
## 72
## CoreID %in% Ring %in% Date
## 144
## Depth %in% CoreID %in% Ring %in% Date
## 699# investigating the homogeneity of variances in the subgroups visually
plot(my.lme, resid(., type = "p") ~ fitted(.) | Date, id = 0.05)
plot(my.lme, resid(., type = "p") ~ fitted(.) | Ring, id = 0.05)
plot(my.lme, resid(., type = "p") ~ fitted(.) | CoreID, id = 0.05)
plot(my.lme, resid(., type = "p") ~ fitted(.) | Depth, id = 0.05)
# the variance across depths does not seem to be homogeneous!!
# updating the model to incorporate changing (in-homogeneous) variance per depth
# the parameter VarIdent allows a different variance per factor level
# VarFixed would allow to use a covariate to describe the variance, i.e. if the soil mositure is a function of air temperature, etc.
# other Var* functions allow to transform the data on the fly.
# see ?VarClasses for details
my.lme2 <- update(my.lme, weights = varIdent(form = ~ 1 | Depth))
summary(my.lme2) # detailed diagnostic output## Linear mixed-effects model fit by REML
## Data: sm[sm$Soil == "HOR", ]
## AIC BIC logLik
## 5909.34 5959.354 -2943.67
##
## Random effects:
## Formula: ~1 | Date
## (Intercept)
## StdDev: 5.708419
##
## Formula: ~1 | Ring %in% Date
## (Intercept)
## StdDev: 0.0001490954
##
## Formula: ~1 | CoreID %in% Ring %in% Date
## (Intercept)
## StdDev: 0.001266887
##
## Formula: ~1 | Depth %in% CoreID %in% Ring %in% Date
## (Intercept) Residual
## StdDev: 7.699998 15.2863
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | Depth
## Parameter estimates:
## 100 200 300 500 900
## 1.0000000 0.1594373 0.6709731 1.6589926 1.5283727
## Fixed effects: Percent_Vol ~ has.wick
## Value Std.Error DF t-value p-value
## (Intercept) 39.18517 2.030118 555 19.301919 0.0000
## has.wickWick 2.78810 1.002167 62 2.782071 0.0071
## Correlation:
## (Intr)
## has.wickWick -0.246
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.5308463 -0.9370398 -0.3408183 0.3044531 2.2635804
##
## Number of Observations: 699
## Number of Groups:
## Date
## 9
## Ring %in% Date
## 72
## CoreID %in% Ring %in% Date
## 144
## Depth %in% CoreID %in% Ring %in% Date
## 699anova(my.lme2) # show anova result only## numDF denDF F-value p-value
## (Intercept) 1 555 425.2186 <.0001
## has.wick 1 62 7.7399 0.0071# Please note how the model result changed compared to the previous simple model when the variance structure is taken into account! See output of the anova(my.lme) above!
# comparing the two models
anova(my.lme, my.lme2)## Model df AIC BIC logLik Test L.Ratio p-value
## my.lme 1 7 6042.316 6074.143 -3014.158
## my.lme2 2 11 5909.340 5959.354 -2943.670 1 vs 2 140.9762 <.0001There is a differene between the two models. The one with the lower AIC value is preferred, in this case the model that includes the term for variable variance. But this does not indicates if the model fits the data well regarding the null hypothesis or if it makes sense in the first place. It only tells us about one “qualitatively better” model, it does not tell us anything quantitatively. AIC is only meaningful comparing two models based on the same data.