# R-program for Supplementary Exercises 12.54 and 12.55 of IPS7e

vit <- read.csv("r:/Chapter 12/ex12_054.csv")
# vitamin E
vite.fit <- lm(vite ~ as.factor(condition), data=vit)
summary(vite.fit)
anova(vite.fit)
# vitamin A
vita.fit <- lm(vita ~ as.factor(condition), data=vit)
summary(vita.fit)
anova(vita.fit)
# same model but different parametrization
vita.fit2 <- lm(vita ~ as.factor(condition)-1, data=vit)
summary(vita.fit2)
# pairwise comparisons: simple method for 1-way ANOVA (does not work for higher order ANOVA)
pairwise.t.test(vit$vita, vit$condition, p.adj="none")
pairwise.t.test(vit$vita, vit$condition, p.adj="bonf")

# general method: manual calculation of pairwise comparisons
g <- 5 # number of groups
coef <- coef(vita.fit2)
vcov <- vcov(vita.fit2)
w<- matrix(c(1,-1,0,0,0, 1,0,-1,0,0, 1,0,0,-1,0, 1,0,0,0,-1, 0,1,-1,0,0, 
             0,1,0,-1,0, 0,1,0,0,-1, 0,0,1,-1,0, 0,0,1,0,-1, 0,0,0,1,-1), ncol=10); 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(vita.fit2)))
# multiple testing
t(p) # unadjusted P-values
p.adjust(p, method="bonferroni")