#
# Exercise 17.6
#
# 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
# Version 1, 11/10/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

# Read data & create factor
fertilizer <- read.table("fertilizer.dat",sep="",header=T)
fertilizer$Block <- as.factor(fertilizer$Block)
fertilizer$Plot <- as.factor(fertilizer$Plot)
fertilizer$fN <- as.factor(fertilizer$N)
summary(fertilizer)

# Plot data
qplot(y=Yield, x=N, colour=Block, data=fertilizer)

# Create dummy variates for blocks 2-4 (first-level-zero parameterization)
fertilizer$d2 <- 1*(fertilizer$Block==2)
fertilizer$d3 <- 1*(fertilizer$Block==3)
fertilizer$d4 <- 1*(fertilizer$Block==4)
fertilizer

# Shape looks exponential but with possible block differences

# First fit exponential curve without block differences, using starting values
# Asymptote = 22 => alpha = 22
# Intercept = 10 => alpha-beta = 10 => beta = 12
# Curve should go through point (Yield=18,N=100) => gamma = 0.011
fert.exp  <- nls(Yield~alpha-beta*exp(-gamma*N), 
                       start=list(alpha=22, beta=10, gamma=0.011), 
                       data=fertilizer)
summary(fert.exp)

# No pre-initialised version for factor + exponential available
# Use dummy variates and assume block differences are zero
fert.exp.2  <- nls(Yield~alpha-beta*exp(-gamma*N)+delta2*d2+delta3*d3+delta4*d4, 
                       start=list(alpha=23.92, beta=13.20, gamma=0.0097, delta2=0, delta3=0, delta4=0), 
                       data=fertilizer)
summary(fert.exp.2)

# Check whether adding blocks has made an improvement to the model
anova(fert.exp, fert.exp.2)

# Check for lack of fit - ResMS equivalent to Blocks + exponential + fac(N)
fert.lm <- lm(Yield~Block+fN, data=fertilizer)
summary(fert.lm)

# Check whether adding lack of fit term makes an improvement to the model (nested models)
anova(fert.exp.2, fert.lm)