# R program for additional exercise 6.1 (VHM 802), slightly updated for contrasts

rat.diet <- read.csv("r:/hs06_1.csv")
library(MASS)
library(doBy) # package to facilitate processing by groups
summaryBy( iron~ diet, data=rat.diet, FUN=c(length,mean,var,min,max))
boxplot( iron~ diet, data=rat.diet)
summary( lm(iron~ as.factor(diet), data=rat.diet))
bc <- boxcox(iron~ as.factor(diet), data=rat.diet, lambda=seq(-1, 1, 0.01))
 
# same thing for a couple of ln and inverse square-root transformations
rat.diet$lniron <- log(rat.diet$iron)
rat.diet$invrootiron <- -1/sqrt(rat.diet$iron)
summaryBy( lniron~ diet, data=rat.diet, FUN=c(length,mean,var,min,max))
boxplot( lniron~ diet, data=rat.diet)
summary( lm(lniron~ as.factor(diet), data=rat.diet))
summaryBy( invrootiron~ diet, data=rat.diet, FUN=c(length,mean,var,min,max))
boxplot( invrootiron~ diet, data=rat.diet)
summary( lm(invrootiron~ as.factor(diet), data=rat.diet))

# work with lniron
rat.aov<- lm(lniron~ as.factor(diet), data=rat.diet)
summary(rat.aov)
anova(rat.aov)

# contrasts; print the different objects to follow the calculations
g <- 5 # number of groups
coef <- coef(rat.aov)
vcov <- vcov(rat.aov)
w<- matrix(c(0,-1,0,0,0, 0,1,-2,0,0, 0,2,2,-3,-3, 0,0,0,1,-1), ncol=4); w
est<- t(w)%*%coef
var <- diag(t(w)%*%vcov%*%w)
se <- matrix(sqrt(var), ncol=1)
t<- est/se
p<- 2*(1-pt(abs(t),df.residual(rat.aov)))
ss<- (t^2)*(mse<- deviance(rat.aov)/df.residual(rat.aov))
F <- t^2/(g-1) # Scheffe test
pF <- 1-pf(F, g-1, df.residual(rat.aov))
data.frame(est,se,t,p,ss,F,pF)

# multiple testing
t(p) # unadjusted P-values
p.adjust(p, method="bonferroni")
p.adjust(p, method="holm") # note: different formula for NS values

# added code with more intuitive contrast specification
rat.aov0<- lm(lniron~ as.factor(diet)-1 , data=rat.diet)
summary(rat.aov0)
anova(rat.aov0) # note that ANOVA table is no longer meaningful!

# contrasts; print the different objects to follow the calculations
g <- 5 # number of groups
coef <- coef(rat.aov0)
vcov <- vcov(rat.aov0)
w<- matrix(c(1,-1,0,0,0, 1,1,-2,0,0, 2,2,2,-3,-3, 0,0,0,1,-1), ncol=4); w
est<- t(w)%*%coef
var <- diag(t(w)%*%vcov%*%w)
se <- matrix(sqrt(var), ncol=1)
t<- est/se
p<- 2*(1-pt(abs(t),df.residual(rat.aov)))
ss<- (t^2)*(mse<- deviance(rat.aov)/df.residual(rat.aov))
F <- t^2/(g-1) # Scheffe test
pF <- 1-pf(F, g-1, df.residual(rat.aov))
data.frame(est,se,t,p,ss,F,pF)