# ---------------------------------------------------------------
 # Problem 18.6
 # ---------------------------------------------------------------

#  first read the data directly from Dr. Oehlerts website

Prb18.6 <- read.table("http://www.stat.umn.edu/~gary/book/fcdae.data/pr18.6",header=TRUE)

Prb18.6$A  <- as.factor(Prb18.6$A)
Prb18.6$B  <- as.factor(Prb18.6$B)
Prb18.6$C  <- as.factor(Prb18.6$C)
Prb18.6$D  <- as.factor(Prb18.6$D)
Prb18.6$E  <- as.factor(Prb18.6$E)
Prb18.6$F  <- as.factor(Prb18.6$F)
Prb18.6$G  <- as.factor(Prb18.6$G)
Prb18.6$H  <- as.factor(Prb18.6$H)



Prb18.6.aov  <-  aov(lm(thickness ~ A*B*C*E, data = Prb18.6))

library(gplots)
qqnorm(Prb18.6.aov)

 # with this command, you can click on and label points
 # to get out of this mode, you can select the "R Console" window and hit the Escape button
if (dev.interactive()) qqnorm(Prb18.6.aov, label=TRUE)


 # now fitting the model with only "active-appearing" effects

Prb18.6.lm2 <- lm(thickness ~ A*B*C, data = Prb18.6)

anova(Prb18.6.lm2)

par(mfrow=c(2,2))  
plot(Prb18.6.lm2)
par(mfrow=c(1,1))  





 # ---------------------------------------------------------------
 # ANCOVA problem 
 # ---------------------------------------------------------------

#  first read the data directly from my website

Forbes <- read.table("http://webpages.uidaho.edu/~chrisw/stat507live/companies.txt",header=TRUE)

head(Forbes)


 # look at the data

plot(Forbes$lassets,Forbes$lsales,pch=as.numeric(Forbes$sector),
      col=as.numeric(Forbes$sector),ylab="Log of Sales",xlab="Log of Assets")
legend("bottomright",legend=c("sector 1","sector 2","sector 3","sector 4","sector 5","sector 6"),
        levels(Forbes$sector),pch=1:6,col=1:6)
title("Log Sales by Sector")


 # standard ANCOVA assuming equal slopes

Forbes.lm1 <- lm(lsales ~ lassets +sector , data=Forbes)
anova(Forbes.lm1)
summary(Forbes.lm1)



 # calculate adjustedMeans and standard errors from the effects package

library(effects)
adjustedMeans <- effect("sector",Forbes.lm1) 
adjustedMeans
adjustedMeans$se


 # multiple comparisons using the multcomp package

library(multcomp)
summary(glht(Forbes.lm1,linfct=mcp(sector="Tukey")))


 # check assumptions

windows()
par(mfrow=c(2,2))  
plot(Forbes.lm1)
par(mfrow=c(1,1))  


 # check equal slopes assumption

Forbes.lm2 <- lm(lsales ~ lassets +sector +lassets:sector , data=Forbes)
anova(Forbes.lm2)

 # information for calculating adjusted means

summary(Forbes)

tapply(Forbes$lsales,Forbes$sector,summary)