# simple ANOVA example
chipsahoy=read.table("chips.ahoy",header=T)
safeway=scan("chips.safeway")
keebler=scan("chips.keebler")
cookies.c=cbind(chipsahoy$chips,rep(1,length(chipsahoy$chips)))
cookies.s=cbind(safeway,rep(2,length(safeway)))
cookies.k=cbind(keebler,rep(3,length(keebler)))
cookies=rbind(cookies.c,cookies.s,cookies.k)
cookies=as.data.frame(cookies)
names(cookies)=c("chips","brand")
summary(cookies)

# need to code categorical variables as "factor" in R
cookies$brand=as.factor(cookies$brand)
attach(cookies)
boxplot(split(chips,brand))
anova(lm(chips~brand))

     Analysis of Variance Table
     
     Response: chips
               Df  Sum Sq Mean Sq F value    Pr(>F)    
     brand      2  767.26  383.63  15.559 1.913e-06
     Residuals 81 1997.16   24.66                      



# multilevel factorial desgin example
eggs=read.table("eggs.txt",header=T,skip=12)
summary(eggs)
data.class(eggs$Fat_Content)
data.class(eggs$Lab)

# need to code categorical variables as "factor" in R
eggs$Lab=as.factor(eggs$Lab)
eggs$Technician=as.factor(eggs$Technician)
eggs$Sample=as.factor(eggs$Sample)
data.class(eggs$Lab)
summary(eggs)

# visualize main effects
attach(eggs)
boxplot(split(Fat_Content,Lab))
boxplot(split(Fat_Content,Technician))
boxplot(split(Fat_Content,Sample))

# fit the model
eggs.lm=lm(Fat_Content~Lab*Technician*Sample)
coef(eggs.lm)

# find fitted value for Lab 3, Tech 1, Sample 1
0.585 - 0.165
# find fitted value for Lab 3, Tech 2, Sample 2
0.585 - 0.165 + 0.155 - 0.295 - 0.180 + 0.195 + 0.260 - 0.060
fitted(eggs.lm)  # rows 17, 18, and 23, 24
mean(eggs$Fat_Content[17:18])

# find significant effects from ANOVA table
summary(eggs.lm)
anova(eggs.lm)

      Response: Fat_Content
                            Df  Sum Sq Mean Sq F value    Pr(>F)    
      Lab                    5 0.44302 0.08860 12.3134 5.486e-06 ***
      Technician             1 0.00441 0.00441  0.6126 0.4414570    
      Sample                 1 0.02901 0.02901  4.0313 0.0560491 .  
      Lab:Technician         5 0.24307 0.04861  6.7558 0.0004636 ***
      Lab:Sample             5 0.04907 0.00981  1.3638 0.2727139    
      Technician:Sample      1 0.00163 0.00163  0.2270 0.6380787    
      Lab:Technician:Sample  5 0.08019 0.01604  2.2288 0.0843659 .  
      Residuals             24 0.17270 0.00720                      

# note: residual standard error equals the square root of the mean square error
sqrt(.00720)

# model building
eggs.lm2=lm(Fat_Content~Lab*Technician+Sample)
anova(eggs.lm2)
eggs.lm3=lm(Fat_Content~Lab*Technician)
anova(eggs.lm3)