#
# Exercise 14.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 Exercise 18.8 of Bliss, C.I. (1970) Statistics in Biology: Statistical 
# Methods for Research in the Natural Sciences, Volume 2. 
# Version 1, 01/01/2016

# 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 
foliar <- read.table("foliar.dat",sep="",header=T)
summary(foliar)

# Calculate relative concentration of Ca and K in each sample
foliar$Caconc <- foliar$Ca/(foliar$SampleWt/20)
foliar$Kconc <- foliar$K/(foliar$SampleWt/20)

# Scatterplot matrix & correlations
pairs(IncHt~IncBA+Caconc+Kconc, data=foliar)
round(cor(foliar[5:8]),3)

#
# Fit MLR for IncHt
#
IncHt.mlr <- lm(IncHt~Caconc+Kconc, data=foliar)
summary(IncHt.mlr)
# Check residuals
plot(IncHt.mlr, ask=FALSE)
# Marginal t-tests in table of coefficients => drop Kconc
IncHt.mlr.2 <- lm(IncHt~Caconc, data=foliar)
summary(IncHt.mlr.2)
# Check residuals
plot(IncHt.mlr.2, ask=FALSE)
# Check fitted model plot
plot.obs <- ggplot(data=IncHt.mlr.2, aes(y=IncHt, x=Caconc)) 
plot.obs + geom_point() + geom_line(aes(y=fitted(IncHt.mlr.2))) +
  ggtitle("Observed data with fitted SLR for IncHt")
# Predict for 20mg Ca in 20g sample: Caconc = 20
IncHt.pred <- predict(IncHt.mlr.2, new=data.frame(Caconc=20), 
                      interval="confidence", se.fit=TRUE)
# Predicted value with 95% CI
IncHt.pred$fit

#
# Fit MLR for IncBA
#
IncBA.mlr <- lm(IncBA~Caconc+Kconc, data=foliar)
summary(IncBA.mlr)
# Check residuals
plot(IncBA.mlr, ask=FALSE)
# Marginal t-tests in table of coefficients => drop Caconc
IncBA.mlr.2 <- lm(IncBA~Kconc, data=foliar)
summary(IncBA.mlr.2)
# Check residuals
plot(IncBA.mlr.2, ask=FALSE)
# Check fitted model plot
plot.obs <- ggplot(data=IncBA.mlr.2, aes(y=IncBA, x=Kconc)) 
plot.obs + geom_point() + geom_line(aes(y=fitted(IncBA.mlr.2))) +
  ggtitle("Observed data with fitted SLR for IncBA")
# Predict for 100mg K in 20g sample: Kconc = 100
IncBA.pred <- predict(IncBA.mlr.2, new=data.frame(Kconc=100), 
                      interval="confidence", se.fit=TRUE)
# Predicted value with 95% CI
IncBA.pred$fit

# End of file