source("http://www.uvm.edu/~rsingle/stat5210/data/scripts-5210.R")

#High school and beyond (hsb2) dataset: 200 obs
#scores on standardized tests of reading (read), writing (write), mathematics (math) and social studies (socst).
#demographic info: gender (female), socio-economic status (ses) ethnic background (race)
 dat <- otherdata("hsb2.dat")

#re-code SES as a character variable
 dat$ses2 <- NULL
 dat$ses2[dat$ses==1] <- "low"
 dat$ses2[dat$ses==2] <- "med"
 dat$ses2[dat$ses==3] <- "high"
  
#re-code SES as a factor
 dat$ses2 <- factor(dat$ses2, levels=c("low","med","high"))
 class(dat$ses2)
 
#---Multiple Comparisons
 pairwise.t.test(dat$write, dat$ses2, p.adj = "none")
 pairwise.t.test(dat$write, dat$ses2, p.adj = "bonf")
 
 mc.none <- pairwise.t.test(dat$write, dat$ses2, p.adj = "none")
 mc.none$p.value * 3 #same results 
 
 TukeyHSD(mod2) #NOTE: Tukey's HSD does not work on lm() objects
 
 #alternative ANOVA syntax using aov() which is a wrapper to lm() 
  mod2.aov <- aov(write ~ ses2, data=dat)
  mod2.aov
  summary(mod2.aov)
 
  TukeyHSD(mod2.aov)
 
#---MC procedures - detail  
 m <- tapply(dat$write, dat$ses, mean)  
 n <- tapply(dat$write, dat$ses, length) 
 df<- sum(n)-3
 mse <- 86.4
 
 #None (no correction)
  pairwise.t.test(dat$write, dat$ses2, 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(dat$write, dat$ses2, p.adj = "bonf")
  t.bon <- t.none
  ((1-pt(t.bon, df))*2)*3 #correcting for all 3 pairwise comparisons
 
 #Tukey-Kramer
  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(mod2.aov)
 
  plot(TukeyHSD(mod2.aov))
 
 
 
 
 
 
 #Scheffe
 N    <- tapply(dat$write, dat$ses2, length)
 MEAN <- tapply(dat$write, dat$ses2, mean) 
 SD   <- tapply(dat$write, dat$ses2, sd) 
 round(rbind(N, MEAN, SD), 2)
 
 #compute MSE based on summary data
 #Use it to compute a 95% CI for the avg of (low & med) vs high using Scheffe's method
 
 ((N[1]-1)*SD[1]^2 + (N[2]-1)*SD[2]^2 + (N[3]-1)*SD[3]^2) / (sum(N)-3)
 L.hat <- (MEAN[1]+MEAN[2])/2 - MEAN[3]
 F.crit <- 3.09 #qf(.95, 3-1, df)
 sqrt((3-1)*F.crit)
 SE <- sqrt(mse*(.5^2/N[1] + .5^2/N[2] + 1^2/N[3]))
 
 L.hat - sqrt((3-1)*F.crit) * SE
 L.hat + sqrt((3-1)*F.crit) * SE