# # Exercise 18.3 # # 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 # Version 1, 21/08/2015 # 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(lsmeans) # for calculating predicted means library(ggplot2) # for plotting # Read data & assign factors wetness <- read.table("wetness.dat",sep="",header=T) wetness$Shelf <- as.factor(wetness$Shelf) wetness$Tray <- as.factor(wetness$Tray) wetness$fWetness <- as.factor(wetness$Wetness) summary(wetness) # Plot data on natural scale qplot(y=NInf, x=Wetness, colour=Shelf, shape=Shelf, data=wetness) # and logit scale wetness$logitNInf <- log( (wetness$NInf+1)/(17-wetness$NInf) ) qplot(y=logitNInf, x=Wetness, colour=Shelf, shape=Shelf, data=wetness) # GLM with binomial distribution and logit link - fitting blocks and treatment means wetness.glm <- glm(cbind(NInf,16-NInf) ~ Shelf + fWetness, family=binomial(), data=wetness) summary(wetness.glm) # Check residual plots plot(wetness.glm, ask=FALSE) # Check for over-dispersion - no evidence ResDev <- summary(wetness.glm)$deviance; ResDev ResDF <- summary(wetness.glm)$df.residual; ResDF p <- 1-pchisq(ResDev, df=ResDF); p # Check whether infection related to wetness anova(wetness.glm, test="Chisq") # Fit linear function of wetness (on logit scale) with test for lack of fit wetness.glm.2 <- glm(cbind(NInf,16-NInf) ~ Shelf + Wetness + fWetness, family=binomial(), data=wetness) anova(wetness.glm.2, test="Chisq") # Drop lack of fit term to investigate linear relationship wetness.glm.3 <- glm(cbind(NInf,16-NInf) ~ Shelf + Wetness, family=binomial(), data=wetness) anova(wetness.glm.3, test="Chisq") summary(wetness.glm.3) # Check residuals from linear response model plot(wetness.glm.3, ask=FALSE) # Use lsmeans package to predict for Wetness, averaging across shelves wetness.glm.3.lsmeans <- lsmeans(wetness.glm.3, ~Shelf+Wetness, at=list(Wetness=c(16:72))) wetness.glm.3.lsm.df <- summary(wetness.glm.3.lsmeans); wetness.glm.3.lsm.df # Plot data with fitted values on linear predictor scale plot.obs <- ggplot(data=wetness, aes(y=log((NInf+1)/(17-NInf)), x=Wetness, colour=Shelf) ) plot.obs + geom_point() + facet_wrap(~Shelf) + geom_line(data=wetness.glm.3.lsm.df, aes(y=lsmean, x=Wetness)) + ggtitle("Predictions with data on linear predictor scale") # Plot data with fitted values on proportion scale plot.obs <- ggplot(data=wetness, aes(y=NInf, x=Wetness, colour=Shelf) ) plot.obs + geom_point() + facet_wrap(~Shelf) + geom_line(data=wetness.glm.3.lsm.df, aes(y=16*exp(lsmean)/(1+exp(lsmean)), x=Wetness)) + ggtitle("Predictions with data on natural scale") # Predict response after 36 hours wetness in an average block on linear predictor scale wetness.glm.3.pv <- lsmeans(wetness.glm.3, ~Wetness, at=list(Wetness=36)) wetness.glm.3.pv.df <- summary(wetness.glm.3.pv) wetness.glm.3.pv.df # Back transform prediction and confidence interval wetness.glm.3.pv.df$back <- exp(wetness.glm.3.pv.df$lsmean)/(1+exp(wetness.glm.3.pv.df$lsmean)) wetness.glm.3.pv.df$blwr <- exp(wetness.glm.3.pv.df$asymp.LCL)/(1+exp(wetness.glm.3.pv.df$asymp.LCL)) wetness.glm.3.pv.df$bupr <- exp(wetness.glm.3.pv.df$asymp.UCL)/(1+exp(wetness.glm.3.pv.df$asymp.UCL)) wetness.glm.3.pv.df # End of file