file <- "http://www.uvm.edu/~rsingle/stat211/data/extra/ch6_BatteryLife.txt"
 dat<-read.table(file, header=T, na.strings=c("NA",".")) 

#ANOVA model
 mod.aov <- aov(lifetime ~ as.factor(battery), data=dat)
 summary(mod.aov)

#---Multiple Comparisons
 pairwise.t.test(dat$lifetime, dat$battery, p.adj = "none") #this function does not have a 'data=' option
 pairwise.t.test(dat$lifetime, dat$battery, p.adj = "bonf")
 
 mc.none <- pairwise.t.test(dat$lifetime, dat$battery, p.adj = "none")
 mc.none$p.value * 3 #same results 
 
 TukeyHSD(mod.aov)
   mod.aov <- aov(lifetime ~ as.factor(battery), data=dat)
   TukeyHSD(mod.aov)
   plot(TukeyHSD(mod.aov))


#---Non-parametric test
#Kruskal-Wallis Rank Sum test
 kruskal.test(lifetime ~ as.factor(battery), data=dat)
 
 
#---Unequal Variance ANOVA (Welch's ANOVA) 
 oneway.test(lifetime ~ as.factor(battery), data=dat, var.equal=FALSE)
 
 
 
 
 
 
 
 
 
 
 
 
 
 
#---MC procedures - detail (NOTE: not necessary for class)
#assign local variables accessible without dat$ prefix
 battery <- dat$battery
 lifetime <- dat$lifetime
 
 m <- tapply(lifetime, battery, mean)   #vector of means 
 n <- tapply(lifetime, battery, length) #vector of sample sizes
 df<- sum(n)-3
 mse <- 86.4
 
 #None (no correction)
  pairwise.t.test(lifetime, battery, p.adj = "none")
  t.none <- (m[3]-m[2]) / sqrt(mse*( (1/n[3])+(1/n[2]) ))
  (1-pt(t.none, df))*2 #2-sided p.value

 #Bonferroni
  pairwise.t.test(lifetime, battery, p.adj = "bonf")
  ((1-pt(t.none, df))*2)*3 #correcting for all 3 pairwise comparisons
 
 #Tukey's HSD
  crit.hsd <- qtukey(.95,3,df)/sqrt(2)
  m[3]-m[2]
  m[3]-m[2] - crit.hsd * sqrt(mse*( (1/n[3])+(1/n[2]) ))
  m[3]-m[2] + crit.hsd * sqrt(mse*( (1/n[3])+(1/n[2]) ))
  w.hsd <- (m[3]-m[2])*sqrt(2) / sqrt(mse*( (1/n[3])+(1/n[2]) ))
  1-ptukey(w.hsd, nmeans=3, df) 

  TukeyHSD(mod.aov)