# Exercise 6.5

# 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 Rothamsted Research (A. Karp)
# Version 1, 13/08/2014

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

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

# Read data & assign factors
parasitoids <- read.table('parasitoids.dat', sep="", header=TRUE)
summary(parasitoids)

# Plot data
qplot(data=parasitoids, y=OpticalDensity, x=Treatment)

# Get mean and variance for each treatment group
t.mean <- aggregate(OpticalDensity~Treatment, data=parasitoids, FUN=mean)
t.var <- aggregate(OpticalDensity~Treatment, data=parasitoids, FUN=var)
cbind(t.mean,t.var$OpticalDensity)

# Bartlett test for homogeneity of variances
b.test <- bartlett.test(OpticalDensity~Treatment, data=parasitoids); b.test

# One-way ANOVA of raw data
parasitoids.aov <- aov(OpticalDensity ~ Treatment, data=parasitoids)
summary(parasitoids.aov)

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

# Take log10-transformation
parasitoids$logOD <- log10(parasitoids$OpticalDensity)
summary(parasitoids)

# One-way ANOVA of log10-transformed data excluding Day0
parasitoids.aov.2 <- aov(logOD ~ Treatment, data=parasitoids, subset=(Treatment!="Day0"))
summary(parasitoids.aov.2)

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

# One-way ANOVA of log10-transformed data excluding Day0 and Control
parasitoids.aov.3 <- aov(logOD ~ Treatment, data=parasitoids, 
                         subset=(Treatment!="Day0")&(Treatment!="Control"))
summary(parasitoids.aov.3)

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

# Get table of predcited means with back-transform
parasitoids.lsm <- lsmeans(parasitoids.aov.3, ~Treatment)
parasitoids.lsm.df <- summary(parasitoids.lsm)
parasitoids.lsm.df$bt <- 10^(parasitoids.lsm.df$lsmean)
parasitoids.lsm.df

# End of file