#
# Exercise 11.2
#
# 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 P Lutman, Rothamsted Research.
# Version 1, 3/05/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
density <- read.table("density.dat",sep="",header=T)
density$Block <- as.factor(density$Block)
density$Plot <- as.factor(density$Plot)
density$Brate <- as.factor(density$B)
density$Crate <- as.factor(density$C)
summary(density)

# Plot data - against barley and chickweed rates separately
qplot(y=Grain, x=B, colour=Block, facets=Crate~., data=density)
qplot(y=Grain, x=C, colour=Block, facets=Brate~., data=density)

# Structural component: Block/Plot
# Explanatory component: Brate*Crate - note: not all combinations present
# Check allocation of treaments
table(density$Brate, density$Crate)

# Fit crossed model using intra-block analysis
# Barley before chickweed
density.lm.1 <- lm(Grain~Block+Brate*Crate, data=density)
anova(density.lm.1)

# Check residual plots
plot(density.lm.1, ask=FALSE)

# Chickweed before barley
density.lm.2 <- lm(Grain~Block+Crate*Brate, data=density)
anova(density.lm.2)

# Predictive model includes all explanatory terms
density.means <- lsmeans(density.lm.1, c("Brate","Crate"))
density.means.df <- summary(density.means)
# Remove missing predictions (combinations not available)
density.means.df <- subset(density.means.df,!is.na(lsmean))
density.means.df
str(density.means.df)


# Plot predictions with 95% CIs
# Set up basic plot against chickweed rate
plot.means <- ggplot(data=density.means.df, aes(y=lsmean, x=Crate))
# Plot means with error bars 
plot.means + geom_point(shape=1, size=4) + 
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL, width=0.1)) +
    facet_wrap(~Brate) +
  ggtitle("Predictions vs chickweed density, grouped by barley rate") 

# Set up basic plot against barley rate
plot.means <- ggplot(data=density.means.df, aes(y=lsmean, x=Brate))
# Plot means with error bars 
plot.means + geom_point(shape=1, size=4) + 
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL, width=0.1)) +
  facet_wrap(~Crate) +
  ggtitle("Predictions vs barley density, grouped by chickweed rate") 


# End of file