# Exercise 8.8 # 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, 10/07/2015 # Set working directory - use setwd() function or from Session menu in RStudio # e.g. setwd("d:/stats4biol/data) # load external packages - available from CRAN library(ggplot2) # Read data & assign factors calcium <- read.table('calcium.dat', sep="", header=TRUE) summary(calcium) # Plot data qplot(data=calcium, y=Length, x=Calcium) # define vCalcium as a variate of concentrations calcium$vCalcium <- rep(0,length(calcium$Calcium)) calcium$vCalcium[calcium$Calcium=="A"] <- 1 calcium$vCalcium[calcium$Calcium=="B"] <- 5 calcium$vCalcium[calcium$Calcium=="C"] <- 10 calcium$vCalcium[calcium$Calcium=="D"] <- 20 # Baseline ANOVA calcium.aov <- aov(Length ~ Calcium, data=calcium) summary(calcium.aov) # Residual plots plot(calcium.aov, ask=FALSE) hist(residuals(calcium.aov)) # To partition Calcium SS in ANOVA table, use orthogonal polynomials z <- poly(calcium$vCalcium, degree=3) calcium$lin <- z[,1] calcium$quad <- z[,2] calcium$cub <- z[,3] # ANOVA with Calcium factor replaced by orthogonal polynomials calcium2.aov <- aov(Length ~ 1 + lin + quad + cub, data=calcium) summary(calcium2.aov) # Get coefficients coef <- coef(calcium2.aov) coef # Calculate predictions from quadratic polynomial c <- c(1:20) fit <- coef[1] + coef[2]*c + coef[3]*c*c # Plot predictions with data and fitted means # Set up basic plot using model object plot.obs <- ggplot(data=calcium, aes(y=Length, x=vCalcium) ) # Plot observations and fitted line plot.obs + geom_point() + geom_point(aes(y=fitted(calcium.aov)),colour="blue",size=2,shape=2) + geom_line(aes(y=fit,x=c)) + ggtitle("Observed data with predicted means and quadratic polynomial") # To get coefficients with respect to untransformed variable, refit quadratic calcium3.aov <- aov(Length ~ 1 + vCalcium + I(vCalcium^2), data=calcium) coef(calcium3.aov) # End of file