# Exercise 8.4

# 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, 26/06/2015

# Set working directory - use setwd() function or from Session menu in RStudio
# e.g. setwd("d:/stats4biol/data)

# load external packWeeds - availabel from CRAN
library(ggplot2)
library(lsmeans)  # tables of predicted means & contrasts

# Read data & assign factors
weed <- read.table('weedcompetition.dat', sep="", header=TRUE)
weed$Block <- as.factor(weed$Block)
weed$Plot <- as.factor(weed$Plot)
weed$Weed <- as.factor(weed$Weed)
summary(weed)

# Plot data
qplot(data=weed, y=Yield, x=Density, colour=Weed)

# Create factor for level of seed density within each weed species
# Find density levels used for each species
densBG <- sort(unique(weed$Density[weed$Weed=='BG']))
densCL <- sort(unique(weed$Density[weed$Weed=='CL']))
densCW <- sort(unique(weed$Density[weed$Weed=='CW']))
# Show densities
cbind(densBG,densCL,densCW)
# Create factor for each weed
weed$BG <- match(weed$Density,densBG)  # find position in list of densities
weed$BG[weed$Weed!='BG'] <- 7          # insert 7 where 'not BG'
weed$BG <- as.factor(weed$BG)          # convert to factor
weed$CL <- match(weed$Density,densCL)
weed$CL[weed$Weed!='CL'] <- 7
weed$CL <- as.factor(weed$CL)
weed$CW <- match(weed$Density,densCW)
weed$CW[weed$Weed!='CW'] <- 7
weed$CW <- as.factor(weed$CW)

# Create combined factor with levels 1-6 within each weed species
weed$fDensity <- rep(0,length(weed$Density)) # initialize to zero
weed$fDensity[weed$Weed=='BG'] <- weed$BG[weed$Weed=='BG'] # add values for BG
weed$fDensity[weed$Weed=='CL'] <- weed$CL[weed$Weed=='CL'] # add values for CL
weed$fDensity[weed$Weed=='CW'] <- weed$CW[weed$Weed=='CW'] # add values for CW
weed$fDensity <- as.factor(weed$fDensity)
cbind(weed$Weed,weed$fDensity,weed$BG,weed$CL,weed$CW)
summary(weed)

# Multi-stratum ANOVA
weed.msaov <- aov(Yield ~ Weed/fDensity + Error(Block/Plot), data=weed)
summary(weed.msaov)

# Equivalent ANOVA to get residuals
weed.aov <- aov(Yield ~ Block + Weed/fDensity, data=weed)
summary(weed.aov)

# Residual plots
plot(weed.aov, ask=FALSE)
hist(residuals(weed.aov))

# Treatment means and comparisons
weed.lsm <- lsmeans(weed.aov, ~Weed:fDensity)
# Treatment means
weed.lsm
# Comparisons vs control within each weed type
contrast(weed.lsm, method='trt.vs.ctrl1', by="Weed")

# Split nested factors into three individual weed types
weed2.msaov <- aov(Yield ~ Weed/(BG+CL+CW) + Error(Block/Plot), data=weed)
summary(weed2.msaov)

# End of file