# EX3.1, p. 75
etch.rate <- read.table("etchrate.txt",header=T)
grp.means <- with(etch.rate, tapply(EtchRate,Power,mean))
with(etch.rate, stripchart(EtchRate~Power,vert=T,method="overplot",pch=1))
stripchart(as.numeric(grp.means)~as.numeric(names(grp.means)),pch="x",
cex=1.5,vert=T,add=T)
title(main="Etch Rate data",ylab=expression(paste("Etch Rate (",
ring(A),"/min)")),xlab="RF Power (W)")
legend("bottomright","Group Means",pch="x",bty="n")

etch.rate$RF <- as.factor(etch.rate$Power)
etch.rate.aov <- aov(EtchRate~RF,etch.rate)
summary(etch.rate.aov)  # Table3.4 p. 76

#Ex. 3.3  p. 79
# overall mean
(erate.mean <- mean(etch.rate$EtchRate))
# treatment effects
with(etch.rate, tapply(EtchRate,RF,function(x) mean(x)-erate.mean))
# Or by this
model.tables(etch.rate.aov)

(MSe <- summary(etch.rate.aov)[[1]][2,3])
(SD.pool <- sqrt(MSe/5))
(t.crit <- c(-1,1)*qt(.975,16))
# .95 conf. int for \mu4 , p.79    
grp.means[4]+t.crit*SD.pool  

# Model Adequacy Checking: residual plots
# indep? constant variance? normal? influential points?
# Checking independence : residual v.s. index (time, run order.....)
# Fig. 3.5  p. 83 
windows()
plot(1:20,resid(etch.rate.aov))

# constant variance? (standardized, sqrt) of resid.  v.s. fitted  Fig. 3.6, p. 83 
# normality: QQ plot  Fig. 3.4, p. 81
# influential? resid Leverage plot
# NOTE: 4 standard residual plots given in R for model checking: plot(fittedmodel)
windows()
opar <- par(mfrow=c(2,2),cex=.8)
plot(etch.rate.aov)    
par(opar)

#Looking at the plot in Figure 3.5, no such pattern are visible thus we have no reason
#to reject the independence hypothesis. A more formal test, and an historical ones, is
#called the Durbin-Watson. This procedures aims at testing the serial autocorrelation of
#errors and by default makes use of constant lag of 1. It is readily available in the car
#and lmtest packages.
require(car)
durbinWatsonTest(etch.rate.aov)  # Not given in Textbook, just FYI.
#the most widely recommended test of homoscedasticity is
# Bartlett¡¦s test.  rmk: sensitive to normality
# H0: constant variance...  so LARGE P value => such assumption is not violated
# p. 85 Ex. 3.4
 bartlett.test(EtchRate~RF,data=etch.rate)  

#In case one suspect strong departures from normality, we may use Levene¡¦s test as an
# alternative test for homogeneity of variances. This test is available in the car package.
# LARGE P value => "good"  (not given in Textbook for this example)
leveneTest(etch.rate.aov)  #FYI#

#Finally, the normality of the residuals can be assessed directly using a Q-Q plot as
#in Figure 3.4 (the so-called droite de Henry, in French) where we expect the values to
#lie approximately on the first bisecting line, or using the Shapiro-Wilk¡¦s test.
shapiro.test(resid(etch.rate.aov))

#checking if data from each RF setting normally distributed.
shapiro.test(etch.rate$EtchRate[etch.rate$RF==160])
shapiro.test(etch.rate$EtchRate[etch.rate$RF==180])
shapiro.test(etch.rate$EtchRate[etch.rate$RF==200])
shapiro.test(etch.rate$EtchRate[etch.rate$RF==220])

# Multiple comparisons 
#Given our 4 treatments, we have a set of 4(4-1)/2 comparisons, the null hypothesis
#being H0 : £gi = £gj for a given (i, j) pair of treatment means. There are several
# ways to carry out parametric multiple comparisons within R.
# Tukey's test p.98~99

TukeyHSD(etch.rate.aov)

# Family-wise .95 conf. intervals in graph
windows()
plot(TukeyHSD(etch.rate.aov),las=1)