# 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)