# # Exercise 15.1 # # 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 R Harrington & C Shortall, Rothamsted Research # Version 1, 01/01/2016 # Set working directory - use setwd() function or from Session menu in RStudio # e.g. setwd("d:/stats4biol/data") # Set up packages to be used later - available from CRAN library(ggplot2) # for plotting # Read data allsites <- read.table("allsites.dat",sep="",header=T,na="*") summary(allsites) # Get log transformation of wet weights allsites$logWt <- log10(allsites$WetWeight+0.5) # Plot data qplot(y=logWt, x=Year, facets=Site~., data=allsites) # Fit separate lines model on log scale allsites.sl <- lm(logWt~Site*Year, data=allsites) anova(allsites.sl) summary(allsites.sl) # Check residuals plot(allsites.sl, ask=FALSE) # Suggestion of variance heterogeneity - check further qplot(y=residuals(allsites.sl), x=fitted(allsites.sl), colour=Site, data=allsites.sl) # Get fitted model plot # Set up basic plot using model object plot.obs <- ggplot(data=allsites.sl, aes(y=logWt, x=Year) ) # Plot observations and fitted line plot.obs + geom_point() + geom_line(aes(y=fitted(allsites.sl))) + facet_wrap(~Site) ggtitle("Observed data with separate lines model on log scale") # Compare with residuals from untransformed data allsites.sl2 <- lm(WetWeight~Site*Year, data=allsites) plot(allsites.sl2, ask=FALSE) # Extra: taking things further # Combine across non-Hereford sites? allsites$Hereford <- as.numeric(allsites$Site=="Hereford")+1 allsites$Hereford <- as.factor(allsites$Hereford) allsites.sl3 <- lm(logWt~(Hereford/Site)*Year, data=allsites) anova(allsites.sl3) # Final model: separate intercepts, common slope for non-Hereford sites allsites.sl4 <- lm(logWt~Year*Hereford+Hereford:Site, data=allsites) anova(allsites.sl4) summary(allsites.sl4) # End of file