# R program for additional exercise 7.3 (VHM 802)
# demonstration of: pairwise comparisons

fish <- read.csv("r:/hs07_3.csv")
fish.full <- lm(invsurv ~ as.factor(rep)+as.factor(conc)*as.factor(prep), data=fish)
summary(fish.full)
fish.red <- lm(invsurv ~ as.factor(conc)+as.factor(prep)+as.factor(rep), data=fish)
summary(fish.red)

# conc comparisons, using contrasts
coef <- coef(fish.red)
vcov <- vcov(fish.red)
w<- matrix(c(0,1,0,0,0,0,
             0,0,1,0,0,0,
             0,1,-1,0,0,0), ncol=3); 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(fish.red)))
data.frame(est,se,t,p)
p.adjust(p, method="bonferroni")
p.adjust(p, method="holm")

# prep comparisons, using contrasts
coef <- coef(fish.red)
vcov <- vcov(fish.red)
w<- matrix(c(0,0,0,1,0,0,
             0,0,0,0,1,0,
             0,0,0,1,-1,0), ncol=3); 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(fish.red)))
data.frame(est,se,t,p)
p.adjust(p, method="bonferroni")
p.adjust(p, method="holm")