################################################################################
## Chapter 11
## January 21, 2009
## Alan T. Arnholt
options(width=65)
library(PASWR)
################################################################################
# Motivational Example: Tires
population <- rep(LETTERS[1:4],6)
Treatment <- sample(population)
datf <- data.frame(Run=1:24,Treatment)
datf[1:5,]
#
attach(Tire)
oneway.plots(StopDist, tire)
detach(Tire)
################################################################################
# Example 11.1
attach(Tire)
TreatmentMean <- tapply(StopDist,tire,mean)
TreatmentMean
a <- length(TreatmentMean)
N <- length(StopDist)
dft <- a-1
dfe <- N-a
GrandMean <- mean(StopDist)
GrandMean
SStreat <- 6*sum((TreatmentMean-GrandMean)^2)
SStreat
SStotal <- sum((StopDist - GrandMean)^2)
SStotal
SSerror <- SStotal - SStreat
SSerror
MStreat <- SStreat/dft
MStreat
MSerror <- SSerror/dfe
MSerror
Fobs <- MStreat/MSerror
Fobs
pvalue <- 1-pf(Fobs,3,20)
pvalue
summary(aov(StopDist~tire))
#
model.tables(aov(StopDist~tire),type="means")
detach(Tire)
################################################################################
# Example 11.2
# part a.
MEANS <- c(405,390)
a <- length(MEANS)
n <- 6
N <- a*n
dfe <- N-a
SD <- 10
alpha <- .05
Y <- rep(MEANS,rep(n,a))
treat <- factor(rep(1:a,rep(n,a)))
SStreat <- summary(aov(Y~treat))[[1]][1,2]
lambda <- SStreat/SD^2
Gamma <- sqrt(lambda)
Gamma
cv <- qt(1-alpha,dfe)
Power <- 1-pt(cv,dfe,ncp=Gamma)
Power
#
power.t.test(n=6,delta=15,sd=10,alternative="one.sided")
################################################################################
# part b.
### Sampling Distribution of MST/MSE
set.seed(10)
sims <- 10000       # number of simulations
n1 <- 6; n2 <- 6; n3 <- 6; n4 <- 6
a <- 4              # number of treatments
N <- n1+n2+n3+n4
df.treat <- a - 1   # dof treatment
df.error <- N - a   # dof error
alpha <- 0.05       # alpha level
### Normal Distribution
mu1 <- 390; sig1 <- 20 # pop1 mean and sigma
mu2 <- 405; sig2 <- 20 # pop2 mean and sigma
mu3 <- 415; sig3 <- 20 # pop3 mean and sigma
mu4 <- 410; sig4 <- 20 # pop4 mean and sigma
t1 <- matrix(rnorm(sims*n1,mu1,sig1),nrow=sims,byrow=TRUE)
t2 <- matrix(rnorm(sims*n2,mu2,sig2),nrow=sims,byrow=TRUE)
t3 <- matrix(rnorm(sims*n3,mu3,sig3),nrow=sims,byrow=TRUE)
t4 <- matrix(rnorm(sims*n4,mu4,sig4),nrow=sims,byrow=TRUE)
MUS <- c(mu1,mu2,mu3,mu4)
MUB <- mean(MUS)
lambda <-(n1*(mu1-MUB)^2 + n2*(mu2-MUB)^2 + n3*(mu3-MUB)^2 +
n4*(mu4-MUB)^2)/(sig1^2)
mt1 <- apply(t1,1,mean)
mt2 <- apply(t2,1,mean)
mt3 <- apply(t3,1,mean)
mt4 <- apply(t4,1,mean)
##############################################################################
mmean <- cbind(mt1,mt2,mt3,mt4)
TT <- cbind(t1,t2,t3,t4)
gm <- apply(TT,1,mean)
SStreat <- n1*((mt1-gm)^2) + n2*((mt2-gm)^2) + n3*((mt3-gm)^2) +
n4*((mt4-gm)^2)
JU2 <- (TT - gm)^2
SStotal <- apply(JU2,1,sum)
SSerror <- SStotal - SStreat
Fobs <- (SStreat/df.treat)/(SSerror/df.error)
q995 <- quantile(Fobs,.995)
##############################################################################
# Figure 11.5
hist(Fobs,col="pink",prob=TRUE,breaks="Scott",main="",xlim=c(0,q995))
title(main="Simulated Sampling Distribution")
curve(df(x,df.treat,df.error,lambda),0,q995,col="red",
add=TRUE,lwd=3)   # only R
val <- c(.80,.85,.90,.95,.99)
Theoretical <- qf(val,df.treat,df.error,lambda)
Simulated <- quantile(Fobs,val)
SimSigLev <- c( sum(Fobs>Theoretical[1])/sims,
sum(Fobs>Theoretical[2])/sims,
sum(Fobs>Theoretical[3])/sims,
sum(Fobs>Theoretical[4])/sims,
sum(Fobs>Theoretical[5])/sims )
TheSigLev <- 1-val
ANS <- rbind(Theoretical,Simulated,TheSigLev,SimSigLev)
round(ANS,4)
################################################################################
Simulated.Power <- sum(Fobs > qf(1-alpha,df.treat,df.error))/sims
Simulated.Power
####
MEANS <- c(390,405,415,410)
a <- length(MEANS)
n <- 6
N <- a*n
SD <- 20
alpha <- .05
Y <- rep(MEANS,rep(n,a))
treat <- factor(rep(1:a,rep(n,a)))
SStreat <- summary(aov(Y~treat))[[1]][1,2] # For R
# SStreat <- summary(aov(Y~treat))[1,2] # For S-PLUS
lambda <- SStreat/SD^2
lambda
cv <- qf(1-alpha,a-1,N-a)
Power <- 1-pf(cv,a-1,N-a,ncp=lambda)
Power
#
power.anova.test(groups=a,n=n, between.var=var(MEANS),within.var=SD^2)
# part c.
MEANS <- c(390,405,415,410)
a <- length(MEANS)
n <- 6
N <- a*n
SD <- 10
alpha <- .05
Y <- rep(MEANS,rep(n,a))
treat <- factor(rep(1:a,rep(n,a)))
SStreat <- summary(aov(Y~treat))[[1]][1,2] # For R
# SStreat <- summary(aov(Y~treat))[1,2] # For S-Plus
lambda <- SStreat/SD^2
lambda
cv <- qf(1-alpha,a-1,N-a)
Power <- 1-pf(cv,a-1,N-a,ncp=lambda)
Power
# Or
SD <- 10
power.anova.test(groups=a,n=n, between.var=var(MEANS),within.var=SD^2)
# part d.
SD <- 20
Powerr <- 0
nr <- 1
MEANS <- c(390,405,415,410)
a <- length(MEANS)
while(Powerr < .80)
{
nr <- nr+1
Nr <- a*nr
alpha <- .05
Yr <- rep(MEANS,rep(nr,a))
treatr <- factor(rep(1:a,rep(nr,a)))
SStreatr <- summary(aov(Yr~treatr))[[1]][1,2] # R
# SStreatr <- summary(aov(Yr~treatr))[1,2] # S-PLUS
lambdar <- SStreatr/SD^2
cvr <- qf(1-alpha,a-1,Nr-a)
Powerr <- 1-pf(cvr,a-1,Nr-a,ncp=lambdar)
}
c(nr,lambdar,Powerr)
# Or
power.anova.test(groups=4, power=.80,
between.var=var(MEANS), within.var=20^2)
#
nrf <- ceiling(power.anova.test(groups=4,power=.80,
between.var=var(MEANS),within.var=20^2)$n)
nrf
# part e.
MEANS <- c(390,405,415,410)
SD <- 14
a <- length(MEANS)
n1=6; n2=6; n3=12; n4=12;
N <- n1+n2+n3+n4
alpha <- .05
Y <- rep(MEANS,c(n1,n2,n3,n4))
treat <- factor(rep(1:a,c(n1,n2,n3,n4)))
SStreat <- summary(aov(Y~treat))[[1]][1,2] # R
# SStreat <- summary(aov(Y~treat))[1,2] # S-PLUS
lambda <- SStreat/SD^2
lambda
cv <- qf(1-alpha,a-1,N-a)
Power <- 1-pf(cv,a-1,N-a,ncp=lambda)
Power
################################################################################
# Checking for Independence of Errors
# Figure 11.7
attach(Tire)
par(pty="s")
mod.aov <- aov(StopDist~tire)
library(MASS)
r <- stdres(mod.aov)
n <- length(StopDist)
plot(1:n,r,ylab="Standardized Residual",xlab="Ordered Value")
detach(Tire)

################################################################################
# Checking for Normality of Errors
# Figure 11.8
attach(Tire)
mod.aov <- aov(StopDist~tire)
library(MASS)
r <- stdres(mod.aov)
par(pty="s")
qqnorm(r)
abline(a=0,b=1)
shapiro.test(r)
detach(Tire)
################################################################################
# Checking for Constant Variance
attach(Tire)
mod.aov <- aov(StopDist~tire)
library(MASS)
r <- stdres(mod.aov)
tm <- fitted(mod.aov)
plot(tm,r,xlab="Fitted Value",ylab="Standardized Residual")
med <- tapply(StopDist,tire,median)
ZIJ <- abs(StopDist-med[tire])
summary(aov(ZIJ~tire))
checking.plots(mod.aov)
detach(Tire)
################################################################################
# Example 11.3
attach(FCD)
FCD.aov <- aov(Weight~Diet)
# Figure 11.12
checking.plots(FCD.aov)
# Figure 11.13
library(MASS)
boxcox(FCD.aov)
detach(FCD)
# Welch's Test
attach(FCD)
ni <- tapply(Weight,Diet,length)
a <- length(ni)
si2 <- tapply(Weight,Diet,var)
wi <- ni/si2
yb <- tapply(Weight,Diet,mean)
ytild <- sum(wi*yb)/sum(wi)
wlamb <- 3*sum((1 - (wi/sum(wi)))^2 / (ni-1) )/(a^2-1)
dfn <- (a-1)
dfd <- 1/wlamb
W <- sum(wi*(yb-ytild)^2/(3-1)) / (1 + 2/3*(3-2)*wlamb)
W
pvalue <- 1-pf(W,dfn,dfd)
pvalue
oneway.test(Weight~Diet)
detach(FCD)
################################################################################
# Example 11.4
attach(Drosophila)
summary(aov(Fecundity~Line))
#
checking.plots(aov(Fecundity~Line))
#
MSE <- summary(aov(Fecundity~Line))[[1]][2,3] #Remove [[1]] for S-PLUS
SSTreat <- summary(aov(Fecundity~Line))[[1]][1,2] # [[1]] for R
ybari <- tapply(Fecundity,Line,mean)
dfe <- 75-3
ni <- c(25,25,25)
ci <- c(1,-.5,-.5)
di <- c(0,1,-1)
ORTHO <- sum(ci*di/ni) # verify orthogonality
ORTHO
SSC1 <- (sum(ci*ybari))^2/sum((ci^2/ni))
SSC2 <- (sum(di*ybari))^2/sum((di^2/ni))
OSSC <- c(SSC1,SSC2)
c(SSC1,SSC2,SSC1+SSC2,SSTreat)
FC1 <- SSC1/MSE
FC2 <- SSC2/MSE
Fs <-c(FC1,FC2)
pval <- 1-pf(Fs,1,dfe)
cbind(OSSC,Fs,pval)
contrasts(Line)[,1]<- ci
contrasts(Line)[,2]<- di
CO <- contrasts(Line)
colnames(CO) <- c("C1","C2")
CO
summary(aov(Fecundity~C(Line,CO,1)+C(Line,CO,2)))
#
contrasts(Line) <- contr.treatment(levels(Line)) # R defaults
contrasts(Line)
# Helmert contrasts
contrasts(Line) <- contr.helmert(levels(Line))
contrasts(Line)
# 
contrasts(Line) <- contr.helmert(levels(Line))[3:1,2:1]
contrasts(Line)
#
CO <- contrasts(Line)
colnames(CO) <- c("Contrast 1","Contrast 2")
CO
summary(aov(Fecundity~C(Line,CO,1)+C(Line,CO,2)))
# Simultaneous p-values
library(multcomp)
summary(glht(aov(Fecundity~Line), linfct = mcp(Line = t(CO))))
CI <- confint(glht(aov(Fecundity~Line), linfct = mcp(Line = t(CO))))
CI
#
##############
# Some Code for Figure 11.15
library(gregmisc)
par(mfrow=c(2,2),mgp=c(3,2,0),pty="m")  # Making room for labels
# All means with 95% CIs
MEANS <- tapply(Fecundity,Line,mean)
mod.dro <- aov(Fecundity~Line)
MSE <- summary(mod.dro)[[1]][2,3]
NS <- tapply(Fecundity,Line,length)
SE <- sqrt(MSE)/sqrt(NS)
t.v <- qt(.975,sum(NS)-length(NS))
ci.l <- MEANS - t.v*SE
ci.u <- MEANS + t.v*SE
barplot2(MEANS,plot.ci=T,ci.l=ci.l,ci.u=ci.u,col="skyblue",
ylim=c(0,40),ci.lwd=2)
title(main="Means with 95% CIs")
# Contrast 1
MEANS <- c(mean(Fecundity[Line=="Nonselected"]),
mean(Fecundity[Line=="Resistant"|Line=="Susceptible"]))
NS <- c(length(Fecundity[Line=="Nonselected"]),
length(Fecundity[Line=="Resistant"|Line=="Susceptible"]))
SE <- sqrt(MSE)/sqrt(NS)
t.v <- qt(.975,sum(NS)-length(NS))
ci.l <- MEANS - t.v*SE
ci.u <- MEANS + t.v*SE
NAMES <- c(expression(hat(mu)[{Nonselected}]),expression(frac(hat(mu)[{Resistant}]+hat(mu)[{Selected}],2)))
barplot2(MEANS,plot.ci=T,ci.l=ci.l,ci.u=ci.u,col="skyblue",names.arg=NAMES,
ylim=c(0,40),ci.lwd=2)
title(main="Individual 95% CIs \n (Contrast 1)")
# Contrast 2
MEANS <- tapply(Fecundity,Line,mean)
mod.dro <- aov(Fecundity~Line)
MSE <- summary(mod.dro)[[1]][2,3]
NS <- tapply(Fecundity,Line,length)
SE <- sqrt(MSE)/sqrt(NS[2:3])
t.v <- qt(.975,sum(NS)-length(NS))
ci.l <- MEANS[2:3] - t.v*SE[1:2]
ci.u <- MEANS[2:3] + t.v*SE[1:2]
NAMES <- c(expression(hat(mu)[{Resistant}]),expression(hat(mu)[{Susceptible}]))
barplot2(MEANS[2:3],plot.ci=T,ci.l=ci.l,ci.u=ci.u,col="skyblue",
names.arg=NAMES,ylim=c(0,40),ci.lwd=2)
title(main="Individual 95% CIs \n (Contrast 2)")
par(mgp=c(3,1,0))
# Simultaneous CI
library(multcomp)
CI <- confint(glht(aov(Fecundity~Line), linfct = mcp(Line = t(CO))))
plot(CI)
par(mfrow=c(1,1))
# name: BARPLOTci.ps
detach(Drosophila)
################################################################################
# Example 11.5
attach(Tire)
CI <- TukeyHSD(aov(StopDist~tire),which="tire")
CI
# Figure 11.16
plot(CI,las=1)
#
library(multcompView)
# Figure 11.17
multcompBoxplot(StopDist~tire,data=Tire)
# Code to compute LSD, BSD, HSD, and Scheffe statistics
tire.aov <- aov(StopDist~tire)
alpha.c <- 0.05
summary(tire.aov)
MSE <- summary(aov(tire.aov))[[1]][2,3] # For R
ybari <- tapply(StopDist,tire,mean)
a <- length(ybari)
N <- length(StopDist)
dfe <- N-a
sort(ybari)
TcritLSD <- qt(1-alpha.c/2,dfe)
LSD <- TcritLSD*sqrt(MSE)*sqrt(2/6)
TcritBON <- qt(1-alpha.c/(choose(a,2)*2),dfe)
BON <- TcritBON*sqrt(MSE)*sqrt(2/6)
TcritTUK <- qtukey(1 - alpha.c,a,dfe)/sqrt(2)
HSD <- TcritTUK*sqrt(MSE)*sqrt(2/6)
CSF <- sqrt((a-1)*qf(1-alpha.c,a-1,dfe))
SCH <- CSF*sqrt(MSE*2/6)
c(LSD,BON,HSD,SCH)
outer(sort(ybari),sort(ybari),"-")
#
# Figure 11.18
library(gregmisc)
NS <- tapply(StopDist,tire,length)
SE <- sqrt(MSE)/sqrt(NS)
t.v <- qt(.975,dfe)
ci.l <- ybari - t.v*SE
ci.u <- ybari + t.v*SE
barplot2(ybari,plot.ci=TRUE,ci.l=ci.l,ci.u=ci.u,col="skyblue",
ylim=c(0,450),ci.lwd=2)
title(main="Mean Stoppping Distance by Tire \n with Individual 95% CIs")
detach(Tire)
################################################################################
# Example 11.6
attach(food)
summary(aov(shear~freezer))
MSC <- summary(aov(shear~freezer))[[1]][1,3] # omit [[1]] for S-PLUS
MSE <- summary(aov(shear~freezer))[[1]][2,3] # omit [[1]] for S-PLUS
sig2tau <- (MSC-MSE)/4
sig2tau
detach(food)
################################################################################
# Randomized Complete Block Design - Sampling Plan
car <- rep(c("Car1","Car2","Car3","Car4"),c(4,4,4,4))
tire <- rep(LETTERS[1:4],c(4,4,4,4))
tireCRD <- sample(tire)
tireCRBD <- c(sample(LETTERS[1:4]),sample(LETTERS[1:4]),
sample(LETTERS[1:4]),sample(LETTERS[1:4]))
Designs <- cbind(car,tire,tireCRD,tireCRBD)
Designs
################################################################################
# Example 11.7
# Figure 11.19
attach(TireWear)
par(mfrow=c(1,2), cex=.8)
interaction.plot(Treat,Block,Wear, type="b", legend=FALSE)
interaction.plot(Block,Treat,Wear, type="b", legend=FALSE)
par(mfrow=c(1,1), cex=1)
# Figure 11.20
library(lattice) # R
A <- stripplot(Treat~Wear|Block,layout=c(4,1))
B <- stripplot(Block~Wear|Treat,layout=c(4,1))
print(A,split=c(1,1,1,2),more=TRUE)
print(B,split=c(1,2,1,2),more=FALSE)
# Figure 11.21
plot.design(Wear~Treat+Block)
# part b.
mod.aov <- aov(Wear~Treat+Block)
summary(mod.aov)
# part c.
yidotbar <- tapply(Wear,Treat,mean)
ydotjbar <- tapply(Wear,Block,mean)
gm <- mean(Wear)
taui <- yidotbar - gm
blockj <- ydotjbar - gm
GM <- matrix(rep(gm,16),nrow=4)
GM
treatm <- matrix(rep(taui,4),nrow=4,byrow=FALSE)
treatm
blockm <- matrix(rep(blockj,4),nrow=4,byrow=TRUE)
blockm
residm <- matrix(resid(aov(Wear~Treat+Block)),
nrow=4,byrow=FALSE)
residm
GM+treatm+blockm+residm
# Or one might use proj()
proj(mod.aov)
# part d.
checking.plots(mod.aov)
# part e.
CI <- TukeyHSD(mod.aov,which="Treat")
CI
# Figure 11.23
plot(CI, las=1)
# Figure 11.24
library(gregmisc)
ybari <- tapply(Wear,Treat,mean)
NS <- tapply(Wear,Treat,length)
MSE <- summary(mod.aov)[[1]][3,3]
SE <- sqrt(MSE)/sqrt(NS)
t.v <- qt(.975,9)
ci.l <- ybari - t.v*SE
ci.u <- ybari + t.v*SE
barplot2(ybari,plot.ci=T,ci.l=ci.l,ci.u=ci.u,col="skyblue",ci.lwd=2)
title(main="Mean Wear by Tire \n with Individual 95% CIs")
# name wearBARPLOT.ps
detach(TireWear)
################################################################################
# Example 11.8
# part a.
Microamps <- c(280,290,285,300,310,295,270,285,290,230,
235,240,260,240,235,220,225,230)
Glass <- factor(c(rep("Glass I",9),rep("Glass II",9)))
Phosphor <- factor(rep(rep(c(rep("Phosphor A",3),
rep("Phosphor B",3),rep("Phosphor C",3)),2)))
# part b.  The data is examined with the function twoway.plots()
# Figure 11.25
twoway.plots(Microamps,Glass,Phosphor)
# part c.
mod.TVB <- aov(Microamps ~ Glass + Phosphor + Glass:Phosphor)
model.tables(mod.TVB,type="means")
#
model.tables(mod.TVB,type="effects")
# part d.
anova(mod.TVB)
# part e.
checking.plots(mod.TVB)
# part f.
# Figure 11.27
par(mfrow=c(1,2),pty="s")
interaction.plot(Glass,Phosphor,Microamps,type="b",legend=FALSE)
interaction.plot(Phosphor,Glass,Microamps,type="b",legend=FALSE)
par(mfrow=c(1,1),pty="m")
# part g.
# Figure 11.28
par(mar=c(5.1,10.1,4.1,2.1),cex.axis=.5)
CI <- TukeyHSD(mod.TVB)
plot(CI,las=1)
par(mar=c(5.1,4.1,4.1,2.1),cex.axis=1)
# Figure 11.29
library(gregmisc)
meanM <- tapply(Microamps,list(Glass,Phosphor),mean)
nsM <- tapply(Microamps,list(Glass,Phosphor),length)
MSE <- anova(mod.TVB)[4,3]
t.c <- qt(.975,12)
lci <- meanM - t.c*sqrt(MSE/nsM)
uci <- meanM + t.c*sqrt(MSE/nsM)
barplot2(meanM,beside=TRUE,legend=TRUE,ylim=c(0,400),
plot.ci=TRUE,ci.l=lci,ci.u=uci,ci.lwd=2,col=c("#A9E2FF","#0080FF"))
title(main="Treatment Means with Individual 95% CIs")
################################################################################