ANOVA with crossed Error Structure splityield
Note: This article shows an R example on how to conduct an analysis of variance (ANOVA) with a split plot experimental design on the dataset splityield. For general information about ANOVA please refer to the following article: ANOVA.
In short: In this article a split plot ANOVA is performed on the dataset splityield to test the effects of irrigation, density and fertilizer and their interactions on the yield. The results are visualized in interaction plots and boxplots. The model is reduced to the minimum adequate model and is evaluated by plotting the residuals.
The Dataset[edit]
The dataset splityield is used as an example in the RBook by Crawley (Crawley, 2013, Analyis of Variance) and can be downloaded on Github under the link provided below. It contains data of a split plot field experiment which contains levels of irrigation, density and fertilizer as well as the yield of the field. Different treatments are nested within each other [Fig. 1]. The variables contained in the dataset are:
- block: A factor containing the block design of the study, which provide the sufficient number of replicates. There are 4 blocks in total.
- irrigation: A factor containing the irrigation which the subplot received (control or irrigated). Irrigation is nested within blocks. Each block received both treatments which leads to one half of the block being irrigated, the other half acting as a control.
- density: A factor containing the density of the soil .Density is nested within irrigation. The 3 densities low, medium, high are randomly distributed across each level of irrigation in each plot.
- fertilizer: A factor with the levels N, P and NP, containing the fertilizer treatment. Fertilizer is nested within density. Each level of density in each level of irrigation in each block is treated with the three fertilizers. These are randomly distributed within each density level.
- yield: A numeric variable containing the yield of each sub sub plot (The blocks of fertilizer treatment.)
Examining the Dataset[edit]
#load necessary libraries
library(ggplot2)
#load data from csv file
splityield <- read.csv("splityield.csv")
str(splityield)
summary(splityield
#visual data inspection
par(mfrow = c(2,2))
hist(splityield$yield)
boxplot(yield ~ fertilizer, data = splityield)
boxplot(yield ~ density, data = splityield)
boxplot(yield ~ irrigation, data = splityield)
par(mfrow = c(1,1))
Output[edit]
> str(splityield)
'data.frame': 72 obs. of 5 variables:
$ yield : int 90 95 107 92 89 92 81 92 93 80 ...
$ block : chr "A" "A" "A" "A" ...
$ irrigation: chr "control" "control" "control" "control" ...
$ density : chr "low" "low" "low" "medium" ...
$ fertilizer: chr "N" "P" "NP" "N" ...
> summary(splityield)
yield block irrigation density fertilizer
Min. : 60.00 Length :72 Length :72 Length :72 Length :72
1st Qu.: 86.00 N.unique : 4 N.unique : 2 N.unique : 3 N.unique : 3
Median : 95.00 N.blank : 0 N.blank : 0 N.blank : 0 N.blank : 0
Mean : 99.72 Min.nchar: 1 Min.nchar: 7 Min.nchar: 3 Min.nchar: 1
3rd Qu.:114.00 Max.nchar: 1 Max.nchar: 9 Max.nchar: 6 Max.nchar: 2
Max. :136.00
>
The variable and different levels can be examined with the str() and summary() commands. What does not get apparent through these commands is the three stage split plot design of the experiment. Knowledge of how the data was obtained is crucial to apply an appropriate statistical design, obtain valid results, and avoid alpha errors resulting from pseudoreplication. See below for more details on pseudoreplication.
To get a quick overview of the data, a histogram of yield (the response variable) can be helpful. To get a more adequate overview one can look at the distribution of yield over the different factor levels to check if the data is strongly skewed at one factor level. Yield is more or less normally distributed across the factor levels [Fig. 2].
ANOVA[edit]
Fitting Model with crossed Error Structure[edit]
An ANOVA is fitted where yield is the dependent variable and irrigation, density and fertilizer as well as their interaction terms are used as predictors. To take the experimental design into account, an error term with a crossed error structure is implemented. The error structure is crossed (not nested, like similar designs) because every factor level of each variable is observed with every factor level of the other variables. The slashes in the error term of the aov() command indicate the hierarchy in which the factors are ordered (density within irrigation and irrigation within block). The factor with the lowest order (fertilizer) does not need to be specified in the error structure. Only the command aov() includes the option for this type of error structure while lm() does not.
#fit ANOVA model with crossed error structure ANOVA_splityield <- aov(yield~irrigation*density*fertilizer+Error(block/irrigation/density), data = splityield) summary(ANOVA_splityield)
Output[edit]
> summary(ANOVA_splityield)
Error: block
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 3 194.4 64.81
Error: block:irrigation
Df Sum Sq Mean Sq F value Pr(>F)
irrigation 1 8278 8278 17.59 0.0247 *
Residuals 3 1412 471
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error: block:irrigation:density
Df Sum Sq Mean Sq F value Pr(>F)
density 2 1758 879.2 3.784 0.0532 .
irrigation:density 2 2747 1373.5 5.912 0.0163 *
Residuals 12 2788 232.3
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
fertilizer 2 1977.4 988.7 11.449 0.000142 ***
irrigation:fertilizer 2 953.4 476.7 5.520 0.008108 **
density:fertilizer 4 304.9 76.2 0.883 0.484053
irrigation:density:fertilizer 4 234.7 58.7 0.680 0.610667
Residuals 36 3108.8 86.4
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
For each level of the hierachichal design, a separate ANOVA table is printed.
There are 4 different blocks, thus the error stratum of block has 3 degrees of freedom. There is no F-value or corresponding p-value as block is not implemented as a fixed effect. This level represents the variation among blocks.
In the second level, the error term is the interaction of irrigation and block. The model has 3 residual degrees of freedom as the degrees of freedom of irrigation (1) and degrees of freedom of block multiplied result in 3. The effect of irrigation can be considered significant.
In the third level, the error is the interaction of block, crossed with irrigation, crossed with density. The residual degrees of freedom are 12 as the degrees of freedom of density and degrees of freedom of blocks multiplied with the two levels of irrigation result in 12. The effect of density can be considered marginally significant. It is however of importance for the model as the interaction if irrigation and density is significant.
Error: Within refers to the lowest error stratum in the model. The effect of fertilizer as well as the interaction between fertilizer and irrigation is significant. The other interaction terms are not.
Model Reduction[edit]
In the lowest hierarchical level, the plots at which the fertilizer treatments were applied, non-significant interactions can be removed. In this example stepwise backwards model reduction based on p-values was applied in which the most complex interactions that were least significant were removed first. After each removal, the model was re-evaluated. Note: Model reduction should typically be done via the AIC. For more information regarding model reduction see: Model reduction. For an example of model reduction via AIC see: Linear Model Reduction LifeCycleSavings.
#remove most most complex, least significant interaction at the lowest level
reduction1 <- aov(yield ~ irrigation*density*fertilizer - irrigation:density:fertilizer
+ Error(block/irrigation/density),
data = splityield)
summary(reduction1)
#remove least significant most complex interaction at lowest level
ANOVA_splityield2 <- aov(yield ~ irrigation*density*fertilizer
- irrigation:density:fertilizer
- density:fertilizer
+ Error(block/irrigation/density),
data = splityield)
summary(ANOVA_splityield2)
Output[edit]
The minimum adequate model obtained by this action is the following:
> summary(ANOVA_splityield2)
Error: block
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 3 194.4 64.81
Error: block:irrigation
Df Sum Sq Mean Sq F value Pr(>F)
irrigation 1 8278 8278 17.59 0.0247 *
Residuals 3 1412 471
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error: block:irrigation:density
Df Sum Sq Mean Sq F value Pr(>F)
density 2 1758 879.2 3.784 0.0532 .
irrigation:density 2 2747 1373.5 5.912 0.0163 *
Residuals 12 2788 232.3
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
fertilizer 2 1977 988.7 11.924 7.28e-05 ***
irrigation:fertilizer 2 953 476.7 5.749 0.00605 **
Residuals 44 3648 82.9
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
Plotting Results[edit]
To visualise the data, interaction plots and a boxplot are printed.
interaction.plot(splityield$fertilizer, splityield$irrigation, splityield$yield)
interaction.plot(splityield$density, splityield$irrigation, splityield$yield)
#boxplot
ggplot(splityield,
aes(x = fertilizer,
y = yield,
fill = density)) +
geom_boxplot() +
facet_wrap(~ irrigation) +
theme_bw()
The interaction of fertilizer and yield is visualised in an interaction plot [Fig.3]. If the lines in an interaction plot are non-parallel, an interaction can be assumed. the plot shows that in the irrigated plots the effects of NP fertilizer increased much more drastically than in the control plots. The model also includes a significant interaction between density and irrigation. Yield is highest on the medium density on the control plots but on the irrigated plots it is highest on the high density plots.
The boxplots visualise the yields for the irrigated group and the control group, grouped by fertilizer and density.
The ANOVA implies that fertilizer density and irrigation all have an influence on the yield. The highest yields can be achieved on the irrigated plots with high density and NP or P fertilizer.
Exemplifying the Need for a crossed Error Structure[edit]
To exemplify the need for a crossed error structure, two ANOVA models are built that do not use an appropriate error structure.
Irrigation <- aov (yield~irrigation, data = splityield) summary(Irrigation) fertilizer <- aov (yield~fertilizer, data = splityield) summary(fertilizer)
Output[edit]
> summary(Irrigation)
Df Sum Sq Mean Sq F value Pr(>F)
irrigation 1 8278 8278 37.43 4.84e-08 ***
Residuals 70 15479 221
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
When the model evaluates the significance of irrigation without the implimentation of an appropriate error structure, the p-value changes and irrigation is regarded as more significant. While it is significant in both cases, there could be models where the predictor becomes significant when leaving out the crossed error structure. This is due a phenomenon called “pseudoreplication”. The table that the ANOVA uses contains 72 observations and thus a high amount of replications. The model assumes 71 overall degrees of freedom based on which the F-value and the p-value are calculated. Replication however specifies how many times the treatment was independently randomised. In this case irrigation was replicated independently only across blocks which equals to four replicates. The assumption of a higher number of replicates than were actually the case in the experiment is called pseudoreplication and can lead to inflated p-values.
> summary(fertilizer)
Df Sum Sq Mean Sq F value Pr(>F)
fertilizer 2 1977 988.7 3.132 0.0499 *
Residuals 69 21779 315.6
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
When testing for the significance of fertilizer without implementing an appropriate error structure, pseudoreplication does not become a problem as fertilizer is applied at the smallest sub plot level. The p-value changes and fertilizer is regarded as less significant. The variation that is explained by plots and other treatments is not included in the model. This leads to inflated residual variation and less significance.
Evaluating the model[edit]
General Residual Inspection[edit]
To evaluate the model, the distribution of residuals can be inspected additionally to the values provided in the ANOVA table. because of the implemented error structure, the plot() command cannot be used with this model. With square brackets the different hierarchy levels of the model can however be extracted and the residuals can be plotted in histograms and normal Q-Q plots.
#extract the different hierachichal levels to plot residuals
model1 <- ANOVA_splityield2[[1]]
summary(model1)
hist(model1$residuals)
model2 <- ANOVA_splityield2[[2]]
summary(model2)
hist(model2$residuals)
model3 <- ANOVA_splityield2[[3]]
summary(model3)
hist(model3$residuals)
par(mfrow = c(2,2))
model4 <- ANOVA_splityield2[[4]]
summary(model4)
hist(model4$residuals)
qqnorm(model4$residuals,
main = "Normal Q-Q Plot model 4")
shapiro.test(model4$residuals)
model5 <- ANOVA_splityield2[[5]]
summary(model5)
hist(model5$residuals)
qqnorm(model5$residuals,
main = "Normal Q-Q Plot model 5")
shapiro.test(model5$residuals)
plot(model5$residuals ~ model5$fitted.values,
main = "residuals vs. fitted model 5")
Output[edit]
- Note: not every output is printed here*
> model4 <- ANOVA_splityield2[[4]]
> summary(model4)
Df Sum Sq Mean Sq F value Pr(>F)
density 2 1758 879.2 3.784 0.0532 .
irrigation:density 2 2747 1373.5 5.912 0.0163 *
Residuals 12 2788 232.3
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> hist(model4$residuals)
> qqnorm(model4$residuals)
> shapiro.test(model4$residuals)
Shapiro-Wilk normality test
data: model4$residuals
W = 0.98846, p-value = 0.9981
> model5 <- ANOVA_splityield2[[5]]
> summary(model5)
Df Sum Sq Mean Sq F value Pr(>F)
fertilizer 2 1977.4 988.7 11.449 0.000142 ***
irrigation:fertilizer 2 953.4 476.7 5.520 0.008108 **
density:fertilizer 4 304.9 76.2 0.883 0.484053
irrigation:density:fertilizer 4 234.7 58.7 0.680 0.610667
Residuals 36 3108.8 86.4
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> hist(model5$residuals)
> qqnorm(model5$residuals)
>
The histograms and normal Q-Q plots show quite normally distributed residuals [Fig. 5]. This implies that there is no big factor influencing data distribution that was missed in the analysis. In the residuals vs. fitted values, no pattern can be observed.
Residual inspection across factor levels[edit]
The residuals should also be inspected across the different factor levels of the model [Fig. 6].
index_5 <- as.numeric(names(model5$residuals)) par(mfrow = c(2,2)) boxplot( model5$residuals ~ splityield$fertilizer[index_5], xlab = "Fertilizer", ylab = "Residuals" ) boxplot( model5$residuals ~ splityield$density[index_5], xlab = "Density", ylab = "Residuals" ) boxplot( model5$residuals ~ splityield$irrigation[index_5], xlab = "Irrigation", ylab = "Residuals" ) boxplot( model5$residuals ~ splityield$block[index_5], xlab = "Block", ylab = "Residuals" )
The model with the within error stratum does not include all the residuals as they are used up for the degrees of freedom. To be able to plot the residuals against the different factors, an index is created that takes the assigned number of all the observations used in model5. This can then be used to only plot the values which are used in this model.
The residuals of the model are approximately normally distributed and have similar variance. It is noteworthy that there seems to be quite the difference in variance among the different blocks. As block A is used for the degrees of freedom, it does not appear in the boxplot.
Further Resources[edit]
Information on ANOVA and split plot designs as well as the splityield dataset can be found in the R book: Analysis of Variance. (2013). In M. J. Crawley, *The R book* (2. ed). Wiley.
The dataset splityield can be downloaded on Github: https://github.com/lmmx/crawley-r-statistics-book-notes/blob/1c291c98d182bb7f8aee3af957d1ccd6cd2efa7a/code/splityield.csv
The author of this entry is Jana Simon