# # Exercise 13.4 # # 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, 1/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 hereford <- read.table("hereford.dat",sep="",header=T) summary(hereford) # Get log transformation of wet weights - as baseline for comparison hereford$logWt <- log10(hereford$WetWeight+0.5) # Fit SLR model on log scale hereford.slr <- lm(logWt~Year, data=hereford) # Check residuals plot(hereford.slr, ask=FALSE) # Get fitted model plot # Set up basic plot using model object plot.obs <- ggplot(data=hereford.slr, aes(y=logWt, x=Year) ) # Plot observations and fitted line plot.obs + geom_point() + geom_line(aes(y=fitted(hereford.slr))) + ggtitle("Observed data with fitted SLR on log scale") # Now compare with original scale # Fit SLR model on original scale hereford.slr.2 <- lm(WetWeight~Year, data=hereford) # Check residuals plot(hereford.slr.2, ask=FALSE) # Get fitted model plot # Set up basic plot using model object plot.obs <- ggplot(data=hereford.slr.2, aes(y=WetWeight, x=Year) ) # Plot observations and fitted line plot.obs + geom_point() + geom_line(aes(y=fitted(hereford.slr.2))) + ggtitle("Observed data with fitted SLR - untransformed") # Log-transform seems slightly preferable as outlier is less pronounced # and less suggestion of heterogeneity # Check for temporal correlation # Index plot qplot(data=hereford.slr, y=rstandard(hereford.slr), x=Year, geom=c("point","line")) # Plot of residuals against previous residual N <- length(hereford$Year) Nm1 <- N-1 plot(y=rstandard(hereford.slr)[1:Nm1], x=rstandard(hereford.slr)[2:N], ,xlab="Previous residual",ylab="Standardized residual") lines(c(-3,3),c(0,0)) lines(c(0,0),c(-3,3)) # Check serial correlation cor(y=rstandard(hereford.slr)[1:Nm1], x=rstandard(hereford.slr)[2:N]) # End of file - but see also Exercise 15.1