#Splus handout: 2-factor ANOVA 

Data:  Tensile strength (psi) of asphaltic concrete specimens for two
aggregate types with each of four compaction methods.  

Note: I have shown the command line code for a two-factor ANOVA
below.  Another (easier) approach is to use the menus, but this does
not offer the studentized residual plots.   

________________________________________________________________________

>concrete<-read.table("DATA/Kuehl/EXPL6-2.DAT",header=T)
>concrete$compact<-factor(concrete$compact)
>concrete$aggregate<-factor(concrete$aggregate)

> concrete
   compact aggregate psi      
 1       2         1 126      Compaction method: 1=static     
 2       2         1 128      			 2=regular    
 3       2         1 133      			 3=low	      
 4       2         2 107      			 4=very low   
 5       2         2 110      				      
 6       2         2 116      Aggregate type:    1=basalt     
 7       3         1  93      			 2=silicious  
 8       3         1 101
 9       3         1  98
10       3         2  63
11       3         2  60
12       3         2  59
13       4         1  56
14       4         1  59
15       4         1  57
16       4         2  40
17       4         2  41
18       4         2  44
19       1         1  68
20       1         1  63
21       1         1  65
22       1         2  71
23       1         2  66
24       1         2  66


#____________________________________________________________________________
#Some EDA

> summary(concrete)
 compact   aggregate         psi        
 1:6       1:12        Min.   : 40.00  
 2:6       2:12        1st Qu.: 59.00  
 3:6                   Median : 66.00  
 4:6                   Mean   : 78.75  
                       3rd Qu.:102.50  
                       Max.   :133.00  
motif()
par(mfrow=c(3,3))

#See also the command:  plot.factor
par(mfrow=c(3,3))
boxplot(split(concrete$psi,concrete$aggregate),xlab="aggregate",
	 ylab="psi",names=c("Basalt","Silicious"))
title("Figure 1")

boxplot(split(concrete$psi,concrete$compact),xlab="compaction method",
	 ylab="psi",names=c("static","regular","low","very low"))
title("Figure 2")

a<-tapply(concrete$psi,list(concrete$aggregate,concrete$compact),mean)    

> a
         1   2        3        4 
1 65.33333 129 97.33333 57.33333
2 67.66667 111 60.66667 41.66667


plot(1:4,a[1,],ylim=range(a),xlab="compaction method",ylab="psi",xaxt="n")
axes(xaxp=c(1,4,3))      
lines(1:4,a[1,],lty=1)
points(1:4,a[2,],pch=2)
lines(1:4,a[2,],lty=2)
title("Figure 3")

#I would add a legend to this plot if it were included in a paper.
#See also the commands: plot.factor and plot.design.

#____________________________________________________________________________
#ANOVA

>con.aov<-aov(psi~aggregate*compact,concrete)
#This is equivalent to:  aov(psi~aggregate+compact+aggregate*compact,concrete)

>summary(con.aov)

                  Df Sum of Sq  Mean Sq  F Value        Pr(F) 
aggregate          1    1734.0 1734.000 182.5263 3.628000e-10
compact            3   16243.5 5414.500 569.9474 0.000000e+00
aggregate:compact  3    1145.0  381.667  40.1754 1.124293e-07
Residuals         16     152.0    9.500                      

#____________________________________________________________________________
#Some diagnostics

e.star<-con.aov$residual/sqrt(sum(con.aov$residual^2)/con.aov$df.residual)
boxplot(split(e.star,concrete$aggregate),xlab="aggregate")
title("Figure 4")

boxplot(split(e.star,concrete$compact),xlab="compaction method")
title("Figure 5")

plot(unclass(concrete$aggregate),e.star)
abline(h=0)
title("Figure 6")

plot(unclass(concrete$compact),e.star)
abline(h=0)
title("Figure 7")

plot(fitted(con.aov),e.star)
abline(h=0)
title("Figure 8")

qqnorm(e.star)
title("Figure 9")