#
# Exercise 4.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 K. Hammond-Kosack, Rothamsted Research
# Version 1, 22/09/2014

# 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)  # tables of predicted means & contrasts

# Read data 
grains <- read.table("grains.dat",sep="",header=T)
summary(grains)

# One-way ANOVA
grains.aov <- aov(Grains~Treatment, data=grains)
summary(grains.aov)

# Tables of predicted means
grains.lsm <- lsmeans(grains.aov, ~Treatment)
grains.lsm
# Get contrasts between means with SE = SED
grains.con <- contrast(grains.lsm, method="pairwise")
grains.con.df <- summary(grains.con)
# Calculate LSD from SED
grains.con.df$LSD <- grains.con.df$SE * qt(0.975,grains.con.df$df)
grains.con.df

# Rest of file verifies results by direct calculations

# Get overall mean
grand.mean <- mean(grains$Grains)

# Get treatment group means
means.df <- aggregate(Grains~Treatment, FUN=mean, data=grains); means.df
# Get number of replicates for each treatment
count.df <- aggregate(Grains~Treatment, FUN=length, data=grains); count.df

# Total SS
TotSS <- sum((grains$Grains-grand.mean)^2); TotSS 
# Treatment SS
grains$means <- means.df$Grains[grains$Treatment]
TrtSS <- sum((grains$means-grand.mean)^2); TrtSS
# Residual SS by subtraction 
ResSS <- TotSS - TrtSS; ResSS
# or by direct calculation
ResSS <- sum((grains$Grains-grains$means)^2); ResSS

# Set up parameters 
N <- length(grains$Grains); N
t <- length(levels(grains$Treatment)); t

# Total df
TotDF <- N-1; TotDF
# Treatment df
TrtDF <- t-1; TrtDF
# Residual df
ResDF <- TotDF - TrtDF; ResDF

# Mean squares
# Treatment mean square 
TrtMS <- TrtSS/TrtDF; TrtMS
# Residual mean square = estimate of s^2
ResMS <- ResSS/ResDF; ResMS

# Variance ratio
VR <- TrtMS/ResMS; VR

# Observed significance level
P <- 1-pf(VR, TrtDF, ResDF); P

# SED
n <- mean(count.df$Grains) # only valid for equal replication
SED <- sqrt(2*ResMS/n); SED

# LSD 
LSD <- SED * qt(0.975, ResDF); LSD

# End of file