#
# Example 18.9
#
# 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 S Foster, Rothamsted Research
# 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
clone <- read.table("clone.dat",sep="",header=T, na.strings="*")
summary(clone)

# Plot data
qplot(data=clone, y=Moving/Total, x=LogDose, facets=Clone~Marker)

# Fit preliminary model
clone.glm <- glm(cbind(Moving,Total-Moving)~LogDose*Clone*Marker, family=binomial(), data=clone)
# Check for over-dispersion
ResDev <- summary(clone.glm)$deviance; ResDev
ResDF <- summary(clone.glm)$df.residual; ResDF
p <- 1-pchisq(ResDev, df=ResDF); p

# Re-fit with estimated dispersion parameter
clone.glm <- glm(cbind(Moving,Total-Moving)~LogDose*Clone*Marker, family=quasibinomial(), data=clone)

# Check residual plots
plot(clone.glm, ask=FALSE)

# Look at sequential ANODEV table
anova(clone.glm, test="F")

# Refine model using marginal F-tests

# Drop 3-way interaction
clone.glm.2 <- glm(cbind(Moving,Total-Moving)~LogDose*Clone*Marker-LogDose:Clone:Marker, 
                   family=quasibinomial(), data=clone)
drop1(clone.glm.2, test="F")

# Drop LogDose:Clone
clone.glm.3 <- glm(cbind(Moving,Total-Moving)~LogDose*Marker+Clone+Clone:Marker, 
                   family=quasibinomial(), data=clone)
drop1(clone.glm.3, test="F")

# Drop LogDose:Marker
clone.glm.4 <- glm(cbind(Moving,Total-Moving)~LogDose+Marker+Clone+Clone:Marker, 
                   family=quasibinomial(), data=clone)
drop1(clone.glm.4, test="F")

# Drop Clone:Marker
clone.glm.5 <- glm(cbind(Moving,Total-Moving)~LogDose+Marker+Clone, 
                   family=quasibinomial(), data=clone)
drop1(clone.glm.5, test="F")

# Drop Marker
clone.glm.6 <- glm(cbind(Moving,Total-Moving)~LogDose+Clone, 
                   family=quasibinomial(), data=clone)
drop1(clone.glm.6, test="F")

# Get predictions from fitted model
clone.glm.6.lsmeans <- lsmeans(clone.glm.6, ~LogDose+Clone, at=list(LogDose=seq(-1.2,0.5,0.1)))
clone.glm.6.lsm.df <- summary(clone.glm.6.lsmeans); clone.glm.6.lsm.df

# Plot data with fitted values on linear predictor scale 
plot.obs <- ggplot(data=clone, aes(y=log(Moving/(Total-Moving)), x=LogDose) ) 
plot.obs + geom_point() + facet_grid(Clone~.) +
  geom_line(data=clone.glm.6.lsm.df, aes(y=lsmean, x=LogDose)) + 
  ggtitle("Predictions with data on linear predictor scale")

# and back-transformed (proportion scale)
plot.obs <- ggplot(data=clone, aes(y=Moving/Total, x=LogDose) ) 
plot.obs + geom_point() + facet_grid(Clone~.) +
  geom_line(data=clone.glm.6.lsm.df, aes(y=exp(lsmean)/(1+exp(lsmean)), x=LogDose)) + 
  ggtitle("Predictions with data on linear predictor scale")

# End of file