multiple_comparisons.R

#
# Start with Basic ANOVA
#
p16.7 <- read.table("CH16PR07.txt",col.names=c("improvement","rd","observation"))
attach(p16.7)
rdf <- factor(rd,levels=c(1:3),labels=c("Low","Moderate","High"))
result <- lm(improvement~rdf)
anova(result)
boxplot(improvement~rdf)

TukeyHSD(result)
TukeyHSD(aov(result))
plot(TukeyHSD(aov(result)))

#
# By Hand?
#

mu.i <- tapply(improvement,rdf,mean)
(n.i <- tapply(improvement,rdf,length))
MSE <- .6401
(se.i <- sqrt(MSE/n.i))

(mult <- 1/sqrt(2)*qtukey(.95,3,24))

mu.i[2]-mu.i[1] + c(-1,1)*mult*sqrt(MSE*(1/n.i[1]+1/n.i[2]))

#
# Other Pairwise Comparisons
# 

pairwise.t.test(improvement,rdf)
pairwise.t.test(improvement,rdf,p.adjust="none")
pairwise.t.test(improvement,rdf,p.adjust="bonferroni")
?p.adjust

#
# Scheffe Contrasts
# 17.4
install.packages("gmodels")
library(gmodels)

cmat <- rbind( "D1" = c(1,-1,0),
		 "D2" = c(0,1,-1),
		 "D3" = c(1,0,-1),
		 "L1" = c(.5,.5,-1))
fit.contrast(result,rdf,cmat)

cmat <- rbind(  "D3" = c(1,0,-1),
		 "L1" = c(.5,.5,-1))
fit.contrast(result,rdf,cmat)
(S <- sqrt((3-1)*qf(.95,3-1,27-3)))
#
# By hand
#
-2.322 + c(-1,1)*S*.4216
-1.694 + c(-1,1)*S*.3712
#
# Let R Help
#
(cres <- as.data.frame(fit.contrast(result,rdf,cmat)))
names(cres)
cres$Est
cres$Est - cres$Std*S  #Lower Bounds
cres$Est + cres$Std*S  #Upper Bounds
#
# Bonferroni
# Family alpha = .05
# Individual?
#
.05/4

cmat1 <- rbind(  "D3" = c(1,0,-1),
		 "L1" = c(.5,.5,-1))
fit.contrast(result,rdf,cmat1)
round(fit.contrast(result,rdf,cmat),3)
cmat2 <- rbind(  "D1" = c(1,-1,0),
		 "D2" = c(0,1,-1))
round(fit.contrast(result,rdf,cmat2),3)

#
# Cofindence Itervals are even eaiser
#
# NOTE THE CHANGE TO 4, NOT 8, IN THE DENOMINATOR
round(fit.contrast(result,rdf,cmat1,conf.int=1-.05/4),3)
round(fit.contrast(result,rdf,cmat2,conf.int=1-.05/4),3)

# Compare
S
qt(1-.05/8,27-3)
#
# Linear Combinations
# Say another post-hoc test is needed, the weighted mean
# g = 5
#
# NOTE THE CHANGE TO 5, from 10, IN THE DENOMINATOR
#
round(fit.contrast(result,rdf,cmat1,conf.int=1-.05/5),3)
round(fit.contrast(result,rdf,cmat2,conf.int=1-.05/5),3)

c.i  <- c(.5,.25,.25)
mu.i*c.i
sum(mu.i*c.i)
(se <- sqrt(MSE*sum(c.i^2/n.i)))
sum(mu.i*c.i)+c(-1,1)*qt(1-.05/10,24)*se