# Exercise 9.6 # Statistical Methods in Biology: Design and Analysis of Experiments and Regression # by S.J. Welham, S.A. Gezan, S.J. Clark & A. Mead (2014) # Chapman & Hall/CRC Press, Boca Raton, Florida. ISBN: 978-1-4398-0878-8 # Data from J. Jenkyn, Rothamsted Research # Version 1, 30/06/2015 # Set working directory - use setwd() function or from Session menu in RStudio # e.g. setwd("d:/stats4biol/data) # load external packages - available from CRAN library(ggplot2) # Read data & assign factors springbarley <- read.table('springbarley.dat', sep="", header=TRUE) springbarley$Block <- as.factor(springbarley$Block) springbarley$MainPlot <- as.factor(springbarley$MainPlot) springbarley$Subplot <- as.factor(springbarley$Subplot) springbarley$NRate <- as.factor(springbarley$NRate) summary(springbarley) # Plot data qplot(data=springbarley, y=Yield, x=MildewF, facets=NForm~NRate) # Multi-stratum ANOVA springbarley.msaov <- aov(Yield ~ NForm*NRate*MildewF+Error(Block/MainPlot/Subplot), data=springbarley) summary(springbarley.msaov) # Tables of effects with SE model.tables(springbarley.msaov, se=TRUE) # Tables of means - not possible to get SEDs?? model.tables(springbarley.msaov, type="means") # Equivalent single-stratum ANOVA (wrong variance ratios) to get residuals springbarley.aov <- aov(terms(Yield ~ Block + NForm*NRate + Block:MainPlot + MildewF/(NForm*NRate), keep.order=TRUE), data=springbarley) summary(springbarley.aov) # Residual plots plot(springbarley.aov, ask=FALSE) hist(residuals(springbarley.aov)) # To partition NForm SS in ANOVA table, need to define new factor to partion # liquid vs solid forms springbarley$LvS <- rep(1,length(springbarley$NForm)) springbarley$LvS[springbarley$NForm=='SS'] <- 2 springbarley$LvS[springbarley$NForm=='SST'] <- 2 springbarley$LvS[springbarley$NForm=='ST'] <- 2 springbarley$LvS <- as.factor(springbarley$LvS) # Multi-stratum ANOVA with partitioned NForm factor - NForm part is the remainder in each case springbarley.msaov2 <- aov(Yield ~ (LvS+NForm)*NRate*MildewF+Error(Block/MainPlot/Subplot), data=springbarley) summary(springbarley.msaov2) # End of file