#
# Example 18.10
#
# 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 Horticulture Research International
# Version 1, 30/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!
rept <- read.table("repeat.dat",sep="",header=T, na.strings="*")
rept$Block <- as.factor(rept$Block)
rept$Cage <- as.factor(rept$Cage)
summary(rept)

# Plot pre-treatment numbers
qplot(data=rept, y=Pre, x=Insecticide, colour=Clone, facets=Block~.)

# Fit preliminary model to pre-treatment numbers
pre.glm <- glm(Pre~Block+Insecticide*Clone, family=poisson(), data=rept)
# Check for over-dispersion
ResDev <- summary(pre.glm)$deviance; ResDev
ResDF <- summary(pre.glm)$df.residual; ResDF
p <- 1-pchisq(ResDev, df=ResDF); p

# Re-fit with estimated dispersion parameter
pre.glm <- glm(Pre~Block+Insecticide*Clone, family=quasipoisson(), data=rept)
# Check residual plots
plot(pre.glm, ask=FALSE)
# Check for treatment differences - none, but clone numbers differ 
anova(pre.glm, test="F")
# Drop Insecticide:Clone
pre.glm.2 <- glm(Pre~Block+Insecticide+Clone, family=quasipoisson(), data=rept)
# Check for treatment differences - differences between clone types
drop1(pre.glm.2, test="F")
# Drop Insecticide
pre.glm.3 <- glm(Pre~Block+Clone, family=quasipoisson(), data=rept)
# Check for treatment differences - differences between clone types
drop1(pre.glm.3, test="F")
summary(pre.glm.3)

# Plot day 2 numbers
qplot(data=rept, y=Day2, x=Insecticide, colour=Clone)

# We will use Pre as an offset in analysis of Day 2 numbers
day2.glm <- glm(Day2~Block+Insecticide*Clone, offset=log(Pre), family=quasipoisson(), 
                data=rept)
# Check residual plots
plot(day2.glm, ask=FALSE)
# Check for treatment differences 
anova(day2.glm, test="F")
# Drop interaction
day2.glm.2 <- glm(Day2~Block+Insecticide+Clone, offset=log(Pre), 
                  family=quasipoisson(), data=rept)
anova(day2.glm.2, test="F")
drop1(day2.glm.2, test="F")
summary(day2.glm.2)

# Get control factor to compare Control vs treated
rept$Control <- as.factor(as.numeric(rept$Insecticide!='Control')+1)
table(rept$Control, rept$Insecticide)

# Fit model with insecticide partitioned into control vs treated
day2.glm.3 <- glm(Day2~Block+Clone+Control/Insecticide, offset=log(Pre), 
                  family=quasipoisson(), data=rept)
anova(day2.glm.3, test="F")

# Simpler model with common insecticde effect
day2.glm.4 <- glm(Day2~Block+Clone+Control, offset=log(Pre), 
                  family=quasipoisson(), data=rept)
summary(day2.glm.4, test="F")

# Get predictions from fitted model on log scale
day2.glm.4.lsmeans <- lsmeans(day2.glm.4, ~Control+Clone )
day2.glm.4.lsm.df <- summary(day2.glm.4.lsmeans); day2.glm.4.lsm.df

# We might use Pre or Day2 as an offset for Day6 and follow a similar procedure

# End of file