# # Example 18.11 # # 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 Clarkson, University of Warwick # Version 1, 31/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(ggplot2) # for plotting library(lsmeans) # for getting predictions with averaging # Read data & assign factors - note repeat is function so don't use as data frame name! carrot <- read.table("carrot.dat",sep="",header=T, na.strings="*") summary(carrot) # Plot counts over time qplot(data=carrot, y=Count, x=Days, facets=Condition~.) # Plot counts over time on logit scale qplot(data=carrot, y=log((Count+1)/(101-Count)), x=Days, facets=Condition~.) # Fit preliminary model carrot.glm <- glm(cbind(Count,100-Count)~Condition*Days, family=binomial(), data=carrot) # Check for over-dispersion ResDev <- summary(carrot.glm)$deviance; ResDev ResDF <- summary(carrot.glm)$df.residual; ResDF p <- 1-pchisq(ResDev, df=ResDF); p # Re-fit with estimated dispersion parameter carrot.glm <- glm(cbind(Count,100-Count)~Condition*Days, family=quasibinomial(), data=carrot) # Check residual plots plot(carrot.glm, ask=FALSE) # Plot fitted model on natural scale plot.obs <- ggplot(data=carrot, aes(y=Count, x=Days)) plot.obs + geom_point() + facet_grid(Condition~.) + geom_line(aes(y=100*fitted(carrot.glm)), data=carrot.glm) + ggtitle("Data with fitted model from glm") # Plot residuals against time qplot(data=carrot.glm, y=residuals(carrot.glm), x=Days, facets=Condition~. ) # Try quadratic function of time carrot.glm.2 <- glm(cbind(Count,100-Count)~Condition*(Days+I(Days^2)), family=quasibinomial(), data=carrot) plot(carrot.glm.2, ask=FALSE) # Still over-dispersion present? ResDev <- summary(carrot.glm.2)$deviance; ResDev ResDF <- summary(carrot.glm.2)$df.residual; ResDF p <- 1-pchisq(ResDev, df=ResDF); p # Plot fitted model on natural scale plot.obs <- ggplot(data=carrot, aes(y=Count, x=Days)) plot.obs + geom_point() + facet_grid(Condition~.) + geom_line(aes(y=100*fitted(carrot.glm.2)), data=carrot.glm.2) + ggtitle("Data with fitted model from glm") # Plot residuals qplot(data=carrot.glm.2, y=residuals(carrot.glm.2), x=Days, facets=Condition~. ) # Look at model more closely anova(carrot.glm.2, test="F") summary(carrot.glm.2) # End of file