##############################################################################
# dae.R
# scripts used in dae.tex to illustrate Montgomery's Design and Analysis
# of Experiments.
# Time-stamp: <2006-12-30 14:01:17 chl>
##############################################################################

##############################################################################
# Chapter 2 (pp. 23--59)
#   This chapter covers:
#   - probability distribution
#   - random sampling
#   - test of hypotheses
##############################################################################
x <- rnorm(10)
set.seed(891)

# Tension Bond Strength data (Tab. 2-1, p. 24)
y1 <- c(16.85,16.40,17.21,16.35,16.52,17.04,16.96,17.15,16.59,16.57)
y2 <- c(16.62,16.75,17.37,17.12,16.98,16.87,17.34,17.02,17.08,17.27)
y <- data.frame(Modified=y1,Unmodified=y2) 
y.means <- as.numeric(apply(y,2,mean))
opar <- par(mfrow=c(2,1),mar=c(5,7,4,2),las=1)
stripchart(y,xlab=expression("Strength (kgf/cm^2)"),pch=19)
arrows(y.means,rep(1.5,2),y.means,c(1.1,1.9),length=.1)
text(y.means,c(1.2,1.8),round(y.means,2),pos=4,cex=.8)
# Random deviates (instead of data from metal recovery used in the book)
rd <- rnorm(200,mean=70,sd=5)
hist(rd,xlab="quantile",ylab="Relative frequency",
     main="Random Normal Deviates\n N(70,5)")
par(opar)

boxplot(y,ylab="Strength (kgf/cm^2)",las=1)

x <- seq(-3.5,3.5,by=0.01)
y <- dnorm(x)
plot(x,y,xlab="",ylab="",type="l",axes=F,lwd=2)
axis(side=1,cex.axis=.8); axis(2,pos=0,las=1,cex.axis=.8)
mtext(expression(mu),side=1,line=2,at=0)
mtext(expression(paste(frac(1, sigma*sqrt(2*pi)), " ",
                       plain(e)^{frac(-(x-mu)^2,
                       2*sigma^2)})),side=3,line=0)
# highlight a specific area (drawing is from left to right, 
# then from right to left)
polygon(c(x[471:591],rev(x[471:591])),c(rep(0,121),rev(y[471:591])),
        col="lightgray",border=NA)

##############################################################################
# Chapter 3 (pp. 60--118)
#   This chapter reviews:
#   - the basis of ANOVA
#   - checking of model assumptions
#   - equivalence with the regression approach
#   - non-parametric alternative

etch.rate <- read.table("etchrate.txt",header=T)
grp.means <- with(etch.rate, tapply(rate,RF,mean))
with(etch.rate, stripchart(rate~RF,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("Observed Etch Rate (",
    ring(A),"/min)")),xlab="RF Power (W)")
legend("bottomright","Group Means",pch="x",bty="n")

# first, we convert each variable to factor
etch.rate$RF <- as.factor(etch.rate$RF)
etch.rate$run <- as.factor(etch.rate$run)
# next, we run the model
etch.rate.aov <- aov(rate~RF,etch.rate)
summary(etch.rate.aov)

# overall mean
(erate.mean <- mean(etch.rate$rate))
# treatment effects
with(etch.rate, tapply(rate,RF,function(x) mean(x)-erate.mean))

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)

with(etch.rate, t.test(rate[RF==160],rate[RF==180],var.equal=TRUE))

mean(tapply(etch.rate$rate,etch.rate$RF,var))/5

mean(c(var(etch.rate$rate[etch.rate$RF==160]),
       var(etch.rate$rate[etch.rate$RF==180])))

confint(etch.rate.aov)

as.numeric(grp.means[4]-grp.means[1])+c(-1,1)*qt(.975,16)*sqrt(2*MSe/5)

opar <- par(mfrow=c(2,2),cex=.8)
plot(etch.rate.aov)
par(opar)

require(car)
durbin.watson(etch.rate.aov)

bartlett.test(rate~RF,data=etch.rate)

levene.test(etch.rate.aov)

shapiro.test(etch.rate$rate[etch.rate$RF==160])

pairwise.t.test(etch.rate$rate,etch.rate$RF,p.adjust.method="bonferroni") 
pairwise.t.test(etch.rate$rate,etch.rate$RF,p.adjust.method="hochberg")

TukeyHSD(etch.rate.aov)
plot(TukeyHSD(etch.rate.aov),las=1)

grp.means <- c(575,600,650,675)
power.anova.test(groups=4,between.var=var(grp.means),within.var=25^2,
                 sig.level=.01,power=.90)

sd <- seq(20,80,by=2)
nn <- seq(4,20,by=2)
beta <- matrix(NA,nr=length(sd),nc=length(nn))
for (i in 1:length(sd)) 
  beta[i,] <- power.anova.test(groups=4,n=nn,between.var=var(grp.means),
                               within.var=sd[i]^2,sig.level=.01)$power
colnames(beta) <- nn; rownames(beta) <- sd
opar <- par(las=1,cex=.8)
matplot(sd,beta,type="l",xlab=expression(sigma),ylab=expression(1-beta),col=1,
        lty=1)
grid()
text(rep(80,10),beta[length(sd),],as.character(nn),pos=3)
title("Operating Characteristic Curve\n for a=4 treatment means")
par(opar)

kruskal.test(rate~RF,data=etch.rate)

library(npmc)
# we need to reformat the data.frame with var/class names
etch.rate2 <- etch.rate
names(etch.rate2) <- c("class","run","var")
summary(npmc(etch.rate2),type="BF")

##############################################################################
# Chapter 4 (pp. 119--159)
##############################################################################

x <- scan("vascgraft.txt")
PSI.labels <- c(8500,8700,8900,9100)
vasc.graft <- data.frame(PSI=gl(4,6,24),block=gl(6,1,24),x)
vasc.graft.aov <- aov(x~block+PSI,vasc.graft)

boxplot(x~PSI,data=vasc.graft,xlab="PSI",
        main="Responses averaged over blocks")
boxplot(x~block,data=vasc.graft,xlab="Blocks",
        main="Responses averaged over treatments")

with(vasc.graft, interaction.plot(PSI,block,x,col=1:6))

summary(vasc.graft.aov)

opar <- par(mfrow=c(2,2),cex=.8)
plot(vasc.graft.aov)
par(opar)

# we delete the 10th observation
x2 <- x
x2[10] <- NA
vasc.graft2 <- data.frame(PSI=gl(4,6,24),block=gl(6,1,24),x2)

rocket <- read.table("rocket.txt",header=T)

matrix(rocket$treat,nr=5,byrow=T)

plot(y~op+batch+treat,rocket)

rocket.lm <- lm(y~factor(op)+factor(batch)+treat,rocket)
anova(rocket.lm)

plot.design(y~factor(op)+factor(batch)+treat,data=rocket)

tab.4.21 <- matrix(c(73,NA,73,75,74,75,75,NA,NA,67,68,72,71,72,NA,75),nc=4)
tab.4.21.df <- data.frame(rep=as.vector(tab.4.21),
                          treat=factor(rep(1:4,4)),
                          block=factor(rep(1:4,each=4)))
summary(aov(rep~treat+block+Error(block),tab.4.21.df))

anova(lm(rep~block+treat,tab.4.21.df))

tab.4.21.lm <- lm(rep~block+treat,tab.4.21.df)
treat.coef <- tab.4.21.lm$coef[5:7]
# effect for catalyst 4 (baseline) is missing, so we add it
treat.coef <- c(0,treat.coef)
pairwise.diff <- outer(treat.coef,treat.coef,"-")

summary(tab.4.21.lm)

crit.val <- qtukey(0.95,4,5)

ic.width <- crit.val*0.6982/sqrt(2)
xx <- pairwise.diff[lower.tri(pairwise.diff)]
plot(xx,1:6,xlab="Pairwise Difference
      (95% CI)",ylab="",xlim=c(-5,5),pch=19,cex=1.2,axes=F)
axis(1,seq(-5,5))
mtext(c("4-1","4-2","4-3","1-2","1-3","2-3"),side=2,at=1:6,line=2,las=2)
segments(xx-ic.width,1:6,xx+ic.width,1:6,lwd=2)
abline(v=0,lty=2,col="lightgray")

tab.4.21.lm.crd <- lm(rep~treat,tab.4.21.df)
(summary(tab.4.21.lm.crd)$sig/summary(tab.4.21.lm)$sig)^2

require(lattice)
xyplot(rep~treat|block,tab.4.21.df,
       aspect="xy",xlab="Catalyst",ylab="Response",
       panel=function(x,y) {
         panel.xyplot(x,y)
         panel.lmline(x,y)
       })

require(lme4)
print(tab.4.21.lm <- lmer(rep~treat+(1|block),tab.4.21.df),corr=F)

print(tab.4.21.lm0 <- lmer(rep~-1+treat+(1|block),tab.4.21.df))

coef(tab.4.21.lm)[[1]]$`(Intercept)`
mean(coef(tab.4.21.lm)[[1]][,1])

col <- c(1,1,0,1,1,1,0,0,0,1,0)
perm <- function(x) { 
  s <- length(x)
  m <- matrix(nc=s,nr=s)
  y <- rep(x,2)
  m[,1] <- x
  for (i in 2:s) { m[,i] <- y[i:(s+i-1)] }
  m
}
col.perm <- perm(col)
bib11 <- rbind(rep(0,11),col.perm)
# check that the design is well balanced
apply(bib11[-1,],1,sum)
apply(bib11,2,sum)