# R-program for Exercise 26.33 of PSLS 3e

gfish <- read.csv("r:/Chapter 26/ex26_033.csv")
library(doBy)
summaryBy( Ventilation ~ Acclimation+Temperature, data=gfish, FUN=c(length,mean,sd))

gfish.full <- lm(Ventilation ~ as.factor(Acclimation)*as.factor(Temperature), data=gfish)
summary(gfish.full)
anova(gfish.full)
# residuals analysis
gfish.diag <- gfish # all diagnostics to be stored in adv.diag (for convenience, not at all necessary)
gfish.diag$fit <- fitted(gfish.full) 
gfish.diag$stdres <- rstandard(gfish.full)
plot(stdres~fit, data=gfish.diag)
qqnorm(gfish.diag$stdres)
qqline(gfish.diag$stdres)
library(lattice)
xyplot(fit ~ Temperature, groups=Acclimation, type="b", data=gfish.diag)

# square-root transformation
gfish$rootvent <- sqrt(gfish$Ventilation)
gfish.rfull <- lm(rootvent ~ as.factor(Acclimation)*as.factor(Temperature), data=gfish)
summary(gfish.rfull)
anova(gfish.rfull)
# residuals analysis
gfish.rdiag <- gfish # all diagnostics to be stored in adv.diag (for convenience, not at all necessary)
gfish.rdiag$fit <- fitted(gfish.rfull) 
gfish.rdiag$stdres <- rstandard(gfish.rfull)
plot(stdres~fit, data=gfish.rdiag)
qqnorm(gfish.rdiag$stdres)
qqline(gfish.rdiag$stdres)
xyplot(fit ~ Temperature, groups=Acclimation, type="b", data=gfish.rdiag)
# for the post-ANOVA analysis, we just compute LSD and Bonferroni-adjusted LSD values
LSD <- qt(1-.025,170)*1.44441*sqrt(2/18); LSD
BSD <- qt(1-.025/25,170)*1.44441*sqrt(2/18); BSD
# the easiest way to get pairwise comparisons within the interaction
# is to reparametrize the full model as a one-way ANOVA (see templ13R.txt for details)