code801.Rmd

---
title: "Code801"
author: "Marianne Huebner"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: word_document
---

```{r, echo=FALSE}
setwd("~/Desktop/Teaching/801/Lectures/Rcode")

library(ggplot2)
library(ggthemes)
library(extrafont)
library(graphics)
library(lattice)
library(scatterplot3d)
library(reshape2)
library(lme4)
library(lsmeans)

library(gplots)
library(car)
library(alr3)
library(faraway)

library(lattice)

library(xtable)
```


*Much of this code is based on Christophe Lalanne's R Companion to Montgomery's Design and Analysis of Experiments (2005)*

## Statistical concepts: Probability distributions

For example, for the F distribution the R code format is \tt{df(x, df1, df2, log = FALSE)} where
\tt{df} gives the density, \tt{pf}  the distribution function, \tt{qf} is 
the quantile function, and \tt{rf} generates random deviates.
```{r, prob}

curve(dnorm(x,0,1), -3,3, ylab="normal density", col="blue")
curve(dnorm(x,1,1), -3,3, col="green", add=TRUE)
curve(dnorm(x,0,2), -3,3, col="red", add=TRUE)

#dt(x, df, ncp=0, log = FALSE)
#pt(q, df, ncp=0, lower.tail = TRUE, log.p = FALSE)
#qt(p, df,        lower.tail = TRUE, log.p = FALSE)
#rt(n, df)

curve(dt(x,10), -3,3, ylab="t density", col="green")
curve(dt(x,1), -3,3, col="blue", add=TRUE)
 legend(1,.35, c("df=10", "df=1"),col=c("green","blue"), lty=1)

#dchisq(x, df, ncp=0, log = FALSE)
curve(dchisq(x,2), 0,30, lty=1, ylab="chisquare density", col="green")
curve(dchisq(x,5),  lty=2, col="blue", add=TRUE)
curve(dchisq(x,10),  lty=3, col="red", add=TRUE)
legend(20,.5, c("df=2","df=5","df=10"),cex=1, col=c("green","blue", "red"), lty=1:3)
```



## Chapter 2: Comparison of two groups


```{r, chap2tension}
set.seed(801)
# 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",nclass=200, ylab="Relative frequency",
main="Random Normal Deviates\n N(70,5)") 
#par(opar)


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

```{r chap2ttest}
t.test(y1,y2,var.equal=TRUE)
as.numeric(diff(apply(y,2,mean)))
t.test(y1,y2)

qqnorm(y1)
qqline(y1)

```


```{r chap2hardness}
tmp<-c(7,3,3,4,8,3,2,9,5,4,6,3,5,3,8,2,4,9,4,5)
hardness <- data.frame(y = tmp, tip = gl(2,10))
t.test(y ~ tip, data = hardness, paired = TRUE)

with(hardness, plot(y[tip==1],y[tip==2],xlab="Tip 1",ylab="Tip 2"))
  abline(0,1)
  with(hardness, plot(y[tip==1]+y[tip==2],y[tip==1]-y[tip==2],
                      xlab="Tip 1 + Tip 2",ylab="Tip 1 - Tip 2",ylim=c(-3,3)))
abline(h=0)

#ignoring pairing
t.test(y ~ tip, data = hardness, var.equal = TRUE)

#nonparametric
wilcox.test(y1,y2)
wilcox.test(y~tip,data=hardness,paired=TRUE)
```


## Chapter 3: ANOVA - Comparison of multiple groups


**Example**: Unwanted material on silion coated wafers is removed by an etching process. The figure below visualizes the association of the radio-frequency (RF) power setting with the etch rate.

```{r etchrate}
etch.rate <- read.table("~/Desktop/Teaching/801/Lectures/Rcode/data/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, ylab=""));
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")
```


Radio-frequency and run are factors in the ANOVA analysis.
```{r aov}
# 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)
```

```{r lambs}
diet1<-c(8,16,9)
diet2<-c(9,16,21,11,18)
diet3<-c(15,10,17,6)
lambs<-cbind(c(rep(1,3),rep(2,5),rep(3,4)),
+ c(diet1,diet2,diet3))
colnames(lambs)<-c("diet","wtgain")
lambs<-data.frame(lambs)

 lambs$diet<-factor(lambs$diet, labels=c(1,2,3))
 anova(lm(lambs$wtgain ~ factor(lambs$diet)))
```

Plots
```{r lampbsplots}
 stripchart(lambs$wtgain~lambs$diet, vert=T, pch=16)
#pdf(file="~/Desktop/Teaching/801/Lectures/Anova_Chap3/lambs_diagnostics.pdf")
par(mfrow=c(2,2))
plot(lm(lambs$wtgain ~ factor(lambs$diet)))
par(mfrow=c(1,1))
#dev.off()
```

```{r modelcheck}
par(mfrow=c(2,2), cex=0.8)
plot(etch.rate.aov)
par(mfrow=c(1,1))

# for a subset of predictors use individual plots
plot(fitted(etch.rate.aov), residuals(etch.rate.aov))  #residuals
qqnorm(residuals(etch.rate.aov)); qqline(residuals(etch.rate.aov))  #normal qq-plot   

durbinWatsonTest(etch.rate.aov)  #test for independence of residuals
bartlett.test(rate~RF,data=etch.rate)  #test for homoscedasticity (constant variance)
leveneTest(etch.rate.aov)   #test for homogeneity of variances

shapiro.test(etch.rate$rate[etch.rate$RF==160])  #normality in subgroups

shapiro.test(etch.rate$rate)  
```



**Power and sample size for ANOVA models**
```{r power}

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)


#operating characteristic curve (OCC) = plot power against a parameter
# here a=4, 
# how does sd and sample size influence power?

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

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")
```

**Multiple Comparisons**

```{r comparison}
#comparison of treatment means
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")

#taking into account inflation of type I error
TukeyHSD(etch.rate.aov) 
plot(TukeyHSD(etch.rate.aov),las=1)
```

*Unbalanced design*
There are four different package designs for a new breakfast cereal. Each was tested in 5 stores (20 stores total), but one store had a fire and was dropped from the study. Sales [number of cases] were recorded.

```{r sales}
sales<-c(11, 17, 16, 14, 15,12, 10, 15, 19,11, 23, 20, 18, 17, NA, 27,33,22,26,28)
trt<-rep(1:4, each=5); trt<-factor(trt)
sales.df<-data.frame(sales=sales, trt=trt)
yi=with(sales.df, tapply(sales,list(trt),sum, na.rm=TRUE))
yibar=with(sales.df, tapply(sales,list(trt),mean, na.rm=TRUE))
ni=c(5,5,4,5)

sales.aov<-aov(sales~trt); summary(sales.aov)

#manually
SStrt <- sum(yi^2/ni) - sum(yi)^2/sum(ni)
SST <- sum(sales^2, na.rm=TRUE) - sum(yi)^2/sum(ni)
SSE<-SST-SSTrt

# contrasts design3 vs design4
yibar %*% c(0,0,1,-1)


```

Tensile strength and cotton weight percent (problem 3.10 and 3.11 Montogomery)
```{r cotton}
cotton<-c(15, 20,25,30,35,15,20,25,30,35,15,20,25,30,35,15,20,25,30,35,15,20,25,30,35)
tensile<-c(7,12,14,19,7,7,17,19,25,10,15,12,19,22,11,11,18,18,19,15,9,18,18,23,11)
cloth<-data.frame(cotton=factor(cotton), tensile=tensile)

# 1. calculate group means
cloth.means<-tapply(cloth$tensile, cloth$cotton, mean)
cloth.means

nn<-tapply(cloth$tensile, cloth$cotton, length)

# 2. ANOVA model
cloth.aov<-aov(tensile ~ cotton, cloth)
summary(cloth.aov, intercept=T)

# 3. Pairwise test of mean differences
pairwise.t.test(cloth$tensile, cloth$cotton, p.adjust="none", pool.sd=T)

# 4. Tukey test on all possible pairs
TukeyHSD(cloth.aov, conf.level=0.95)
plot(TukeyHSD(cloth.aov, conf.level=0.95))

# 5. Tukey method step-by-step
       # from ANOVA table
ntrt<-5        # number of treatments
dferror<-20    # N-a
mserror<-8.06  # MSE residuals
contrastcoef<-c(-1,0,0,1,0)  #contrast (30% cotton as control)
meandiff<-sum(contrastcoef*cloth.means); meandiff  #difference in means
sediff<-sqrt( (mserror/2)* sum( (1/nn)*abs(contrastcoef)) ); sediff #standard error of the difference

critq<-meandiff/sediff
qtukey(0.95, ntrt, dferror)  #studentized range distribution quantile
# 95% confidence interval for the contrast
lowci<-meandiff - qtukey(0.95, ntrt, dferror)*sediff
upci<-meandiff + qtukey(0.95, ntrt, dferror)*sediff
lowci;upci

# pvalue
ptukey(meandiff/sediff, ntrt, dferror, lower.tail=F)


# 6. Scheffe method step-by-step
       # from ANOVA table
ntrt<-5        # number of treatments
dfnum<- 4      # a-1
dferror<-20    # N-a
mserror<-8.06  # MSE residuals
contrastcoef<-c(-1,0,0,1,0)  #contrast (30% cotton as control)
meandiff<-sum(contrastcoef*cloth.means); meandiff  #estimated contrast
sscoef<-sum(contrastcoef*contrastcoef/nn) # sum of squared contrast coefficients

msdiff<-meandiff*meandiff/sscoef    # mean squared for contrast
Fcontrast<-msdiff/mserror; Fcontrast

critval <- dfnum* qf(0.05, dfnum, dferror, lower.tail=F)  #Scheffe F

lowci<-meandiff - sqrt(critval)*sqrt(mserror*sscoef)
upci<-meandiff + sqrt(critval)*sqrt(mserror*sscoef)
lowci; upci

# 7. Dunnett method step-by-step

ntrt<-5        # number of treatments
dfnum<- 4      # a-1
dferror<-20    # N-a
mserror<-8.06  # MSE residuals
contrastcoef<-c(-1,0,0,1,0)  #contrast (30% cotton as control)
meandiff<-sum(contrastcoef*cloth.means);   #difference in means
sediff<-sqrt( (mserror/ntrt)); 

# critical value
# library(nCDunnett)
# qNCDun(0.95, nu=dferror, rho=(rep(0.5,times=dfnum)), delta=rep(0,times=dfnum), two.sided=F)

critval<-2.65

# 95% confidence interval for the contrast
lowci<-meandiff - critval*sediff
upci<-meandiff + critval*sediff
lowci;upci

```



## Chapter 4: RCBD - Randomized Complete Block Design

A product developer decides to investigate the effect of four different levels of extrusion pressure on flicks using a RCBD considering batches of resin as blocks.

```{r rcbd}
y<-c(90.3, 89.2, 98.2, 93.9, 87.4, 97.9,
92.5, 89.5, 90.6, 94.7, 87.0, 95.8,
85.5, 90.8, 89.6, 86.2, 88.0, 93.4,
82.5, 89.5, 85.6, 87.4, 78.9, 90.7)
psi.labels <- c(8500,8700,8900,9100)
vasc <- data.frame(psi=gl(4,6,24),batch=gl(6,1,24),y)

boxplot(y~psi, vasc)
boxplot(y~batch, vasc)

interaction.plot(vasc$psi, vasc$batch, vasc$y,  fun=mean, col=1:6, ylab="mean of y", xlab="PSI")

vasc.aov <- aov(y~batch+psi,vasc)
summary(vasc.aov)
```


**Latin square design: Rocket propellant problem**
```{r latin}
rocket<-data.frame(y=c(24, 20, 19, 24, 24, 17, 24, 30, 27, 36, 18, 38, 26, 27, 21, 26, 31, 26, 23, 22, 22, 30, 20, 29, 31), batch=rep(1:5, each=5), op=rep(1:5, 5),        
                   treat=c("A","B","C","D","E",  "B","C","D","E","A", "C", "D","E", "A", "B", "D","E", "A","B","C", "E", "A","B","C","D"))
 
plot(y~op+batch+treat,rocket)
rocket.lm <- lm(y~factor(op)+factor(batch)+treat,rocket) 
anova(rocket.lm)

```


**Paired t-test is  a special case of an RCBD**

```{r paired}
## Paired t-test vs RCBD

#The following data are from a matched pairs design
df<-data.frame(a = c(12.9, 13.5, 12.8, 15.6, 17.2, 19.2, 12.6, 15.3, 14.4, 11.3),
b = c(12.0, 12.2, 11.2, 13.0, 15.0, 15.8, 12.2, 13.4, 12.9, 11.0))

# testing this as a two-sample t-test (incorrect) or a paired t-test (correct) gives different results
t.test(df$a,df$b)
t.test(df$a,df$b, paired=T)

# Recreate this analysis as CRD and RCBD
y =  c(t(as.matrix(df)))  #response
f = c("pre", "post")  # treatment levels
k=length(f)  # number of treatment levels
n=length(df$a)  # number of blocks
tm = gl(k, 1, n*k, factor(f))    # paired treatments
blk = gl(n, k, k*n)    # blocking factor 

# CRD compare to two sample t-test results
summary(aov(y ~ tm) )

#RCBD compare to paired t-test results
summary(aov(y ~ tm + blk) )


```

**BIB: Balanced incomplete block design**

In a catalyst experiment the time of reaction for a chemical process is studied as a function of catalyst type administered to four different batches of raw material. These batches are considered as the blocking elements.

```{r bib}
y <- matrix(c(73,NA,73,75,74,75,75,NA,NA,67,68,72,71,72,NA,75),nc=4) 
chemproc <- data.frame(rep=as.vector(y),
treat=factor(rep(1:4,4)),
batch=factor(rep(1:4,each=4))) 

summary(aov(rep~treat+batch+Error(batch),chemproc))
anova(lm(rep~batch+treat,chemproc))  #treatment effect adjusted for the blocking factor

# batch effect adjusted for treatment
summary(aov(rep~treat+batch+Error(treat),chemproc))

#Tukey pairwise differences
chemproc.lm <- lm(rep~batch+treat,chemproc)
treat.coef <- chemproc.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(chemproc.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")

# Does this BIB perform better than a complete randomized design (without blocking)?
# relative efficiency sigma^2(CRD)/sigma^2 (RCBD)

chemproc.lm.crd <- lm(rep~treat,chemproc) 
(summary(chemproc.lm.crd)$sig/summary(chemproc.lm)$sig)^2
#Thus CRD would require 13% more bservations to obtain the same level of precision as the BIB

#interbloc variation
require(lattice) 
xyplot(rep~treat|batch,chemproc,
aspect="xy",xlab="Catalyst",ylab="Response", panel=function(x,y) {
panel.xyplot(x,y)
panel.lmline(x,y) })

```

```{r rabbit}
library(faraway)
data(rabbit)

summary(aov(gain~treat+block+Error(block),rabbit))
anova(lm(gain~block+treat,rabbit))

g<- lm(gain~block+treat,rabbit)
g1<-lm(gain~treat, rabbit)

releff<- (summary(g1)$sig/summary(g)$sig)^2
```


## Chapter 5: Factorial Design

```{r battery}
#battery <- read.table("~/Desktop/Teaching/801/Lectures/Rcode/Data/battery.txt",header=TRUE) 
y<-c(130,74,155,180, 150, 159, 188, 126, 138, 168, 110, 160,  34,  80,  40,  75, 136, 106, 122, 115, 174, 150, 120, 139,  20,  82,  70,  58,  25,  58,  70,  45,
                96,  82, 104,  60)
x1 <- as.factor(c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3)) 
x2 <- as.factor(c(rep(15,12), rep(70,12), rep(125,12))) 
battery<-data.frame(life=y, material=x1, temperature=x2)
summary(battery)

#two-way ANOVA including interaction
battery.aov <- aov(life~material*temperature,data=battery)
summary(battery.aov)

#interaction plot
with(battery, interaction.plot(temperature,material,life,type="b",pch=19, fixed=T,xlab="Temperature (◦F)",ylab="Average life"))

# effects plot
plot.design(life~material*temperature,data=battery)
```
From the interaction plot we see that average life decreases as temperature increases, with Material type 3 leading to extended battery life compared to the other, especially at higher temperature.

From the effects plot we see that the large Temperature effect is reflected in the range of battery life variation induced by its manipulation.

```{r batterycomparisons}
TukeyHSD(battery.aov,which="material")

# BUT interaction is significant! So we need to compute means contrast in presence of a significant interaction.
# compute the three means at Temperature=70◦F 
mm <- with(subset(battery, temperature==70),aggregate(life,list(M=material),mean))
# next the studentized t quantile times the error type (based on pooled SD
# from ANOVA)
val.crit <- qtukey(.95,3,27)*sqrt(unlist(summary(battery.aov))[["Mean Sq4"]]/4) # finally compare the observed difference of means with the critical value
diff.mm <- c(d.3.1=mm$x[3]-mm$x[1],d.3.2=mm$x[3]-mm$x[2],d.2.1=mm$x[2]-mm$x[1])
names(which(diff.mm > val.crit))

```
In conclusion, only Material type 3 vs. type 1 and Material type 2 vs. type 1 appear to be significantly different when Temperature is fixed at 70 degrees F.

```{r batteryresiduals}

plot(battery.aov)  # residual plots


with(battery, tapply(life,list(material,temperature),var))

```


Examining the plot of residuals vs. fitted values, we can see that a larger variance is associated to larger fitted value, and two observations (2 and 4) are highlighted. At 15 degrees F material type 1 contains outliers that account for the the inequality in variance.


```{r batteryoutlier}
summary(battery.aov1 <- aov(life~material*temperature,data=battery, subset=-c(2,4)))

plot(battery.aov1)

with(battery[-c(2,4),], tapply(life,list(material,temperature),var))

with(battery[-c(2,4),], interaction.plot(temperature,material,life,type="b",pch=19, fixed=T,xlab="Temperature (◦F)",ylab="Average life"))
```



```{r batterynointeraction}
# no interactions
summary(battery.aov2 <- aov(life~material+temperature,data=battery))

# Plot difference between the observed cell means and the estimated cell means assuming no interaction: any pattern is suggestive of the presence of an interaction.
mm2 <- with(battery, tapply(life,list(material,temperature),mean))
mm2 <- as.vector(mm2) 
plot(fitted(battery.aov2)[seq(1,36,by=4)],mm2-fitted(battery.aov2)[seq(1,36,by=4)], xlab="",ylab=expression(bar(y)[ij.]-hat(y)[ijk]), pch=19,axes=FALSE,ylim=c(-30,30)) 
axis(1,at=seq(0,200,by=50),pos=0) 
text(155,4, expression(hat(y)[ijk]),pos=4) 
axis(2,at=seq(-30,30,by=10),las=1)
# superimpose loess line
yy <- order(fitted(battery.aov2)[seq(1,36,by=4)])
xx.fit <- fitted(battery.aov2)[seq(1,36,by=4)]
yy.fit <- mm2-fitted(battery.aov2)[seq(1,36,by=4)] 
lines(xx.fit[yy],predict(loess(yy.fit[yy]~xx.fit[yy])),col="lightgray",lwd=2)
```

There is a tendency toward alternated fitted values differences with increasing predicted values. That means a no interaction model is not appropriate.

**Adding a quadratic term**

Model: $$y \sim A + B + B^2 + AB + AB^2 $$

```{r batteryquadratic}
# adding a quadratic term for temperature
battery$temp2 <- as.numeric(as.character(battery$temperature))


battery.aov3 <- aov(life~material+temp2+I(temp2^2) +material:temp2+material:I(temp2^2),data=battery) 
summary(battery.aov3)

```

Is there an optimal temperature that depends on material type where battery life reaches its maximum?

```{r batteryint}
#interpolate predictions over the range of temperatures
new <- data.frame(temp2=rep(seq(15,125,by=5),3),material=gl(3,23)) 
new$fit <- predict(battery.aov3,new)

# we first plot the fitted values
with(new, interaction.plot(temp2,material,fit,legend=FALSE, xlab="Temperature",ylab="Life",ylim=c(20,190)))
txt.leg <- paste("Material type",1:3) 
text(5:7,new$fit[new$temp2==c(45,55,65)]-c(3,3,-20),txt.leg,pos=1) 
# next the observed values 
points(rep(c(1,15,23),each=12),battery$life,pch=19)

```

*Contour plot*
This is plotting fitted values against continous predictor values after interpolation. Since material is not continous, this plot is not usual, but it gives some idea about the influence of material type and temperature on battery life.

```{r batterycontour}
# require(lattice)
fit<-lm(life~material+temp2+I(temp2^2) +material:temp2+material:I(temp2^2),data=battery) 

battery.tmp <- data.frame(temp2=rep(seq(15,125,by=5),3),material=gl(3,23)) 
#tmp <- list(material=tmp$material,temp2=tmp$temp2) 
#new <- expand.grid(battery.tmp)
#new$fit <- c(predict(fit,new))

#contourplot(fit~material*temp2,data=new,cuts=8,region=T,col.regions=gray(7:16/16))
```

*Example: Hormone study*
An experiment is conducted to study the effect of hormones injected into test rats. There are two distinct hormones, each with two levels: A, a, B, b. Each treatment is applied to six rats with the response being the amount of glycogen (in mg) in the liver.


```{r hormone}
hlevel<-rep(gl(2,6, labels=c("high", "low")),2)
horm<-gl(2,12, labels=c("A", "B"))
glyc<-c(106, 101,120, 86,132, 97, 51, 98, 85, 50, 111, 72, 103, 84, 100, 83, 110, 91, 50, 66, 61, 72, 85, 60 )
rats<-data.frame(hlevel=hlevel, hormone=horm, glyc=glyc)

rats.lm<-lm(glyc~hormone*hlevel, rats)
anova(rats.lm)

with(rats, tapply(glyc,list(hlevel, hormone),mean))

# interaction plot
with(rats, interaction.plot(hlevel,hormone,glyc,type="b",pch=19, fixed=T,xlab="Level",ylab="Glycogen"))

residualPlots(rats.lm)

rats.lm2<-lm(glyc~hormone+hlevel, rats)
anova(rats.lm2)

residualPlots(rats.lm2)

TukeyHSD(aov(glyc~hormone*hlevel, rats))


# sample size for power 80%, alpha=0.05
p1<-(summary(rats.lm)$sigma)^2
a=2
b=2
D=25
n=4  #n=5
phiA=sqrt(n*b*D^2/(2*a*p1))
phiB=sqrt(n*a*D^2/(2*b*p1))
phiAB=sqrt(n*D^2/(2*((a-1)*(b-1)+1)*p1))
dfnum=1
dfE=a*b*(n-1)
# continue with Table V Appendix Montgomery 

```



*Example: Winery*
A grape grower is interested in maximizing the number of bushels per acre on her winery. She limits her study to the combinations of 3 varieties and 4 pesticides (12 combinations). For each combination, two replicates will be obtained. If there is an interaction the grower wants to compare the pesticides for each fixed variety of grape.

```{r winery}
x1<-rep(1:3, 8)
x2<-rep(1:4, each=6)
y<-c(49, 39, 50 , 55 , 43 ,38 , 85, 73, 55, 41, 67, 58 ,53 , 42 , 53, 48, 66, 68, 85,92 ,69 , 62 , 85 , 99)
wine<-data.frame(bushel=y, variety=x1, pesticide=x2)

with(wine, tapply(bushel,list(variety,pesticide),mean))

wine.lm<-lm(bushel~variety*pesticide, wine)
summary(wine.lm)
anova(wine.lm)

#test interaction   library(alr3)
residualPlots(wine.lm)
#This includes Tukey's test for nonadditivity when plotting against fitted values.

```


**General factorial design**


*Example: Hormone study*
A soft drink bottler is interested in obtaining more uniform fill heights in the bottles produced by his manufacturing process. The process engineer can control three variables during the filling process: the percent carbonation (A), the operating pressure in the filler (B), and the bottles produced per minute on the line speed (C).

$$y \sim A + B + C + AB + AC + BC + ABC$$

```{r bottling}
bottling <- read.table("~/Desktop/Teaching/801/Lectures/Rcode/data/bottling.txt",header=TRUE, colClasses=c("numeric",rep("factor",3)))
summary(bottling)

par(mfrow=c(2,2),cex=.8)
boxplot(Deviation~.,data=bottling,las=2,cex.axis=.8,ylab="Deviation") 
abline(h=0,lty=2)
par(las=1)
mm <- with(bottling, tapply(Deviation,Carbonation,mean))
ss <- with(bottling, tapply(Deviation,Carbonation,sd))
bp <- barplot(mm,xlab="Carbonation",ylab="Deviation",ylim=c(-2,9)) 
arrows(bp,mm-ss/sqrt(4),bp,mm+ss/sqrt(4),code=3,angle=90,length=.1) 
with(bottling, interaction.plot(Carbonation,Pressure,Deviation,type="b")) 
with(bottling, interaction.plot(Carbonation,Speed,Deviation,type="b"))
par(mfrow=c(1,1),cex=1)
```

Top right panel presents aggregated data on B and C: The means observed, together with their standard errors, are plotted against the three level of Carbonation. The bottom interaction plots displays only two factors at the same time.

```{r bottlingaov}
summary(bottling.aov <- aov(Deviation~.^3,bottling))
```

The main effects are significant, while none of the four interaction effects are.
Remove interaction effects:

```{r bottlingaov2}
bottling.aov2 <- aov(Deviation~.,bottling) 
anova(bottling.aov2,bottling.aov)
```

Since the F-test is not significant, we can assume the interaction terms need not be part of the model.

*Example: Tool life*

The effect of cutting speed (A) and tool angle (B) is studied on the effective life of a cutting tool. This is an example with several continuous factors.

The hierarchical principle means that if a higher order term i sincluded, then teh lower order terms should also be included.

$$y \sim A + B + A^2 + B^2 + AB^2 + A^2B + AB$$

```{r toollife}
tool <- read.table("~/Desktop/Teaching/801/Lectures/Rcode/Data/toollife.txt",header=TRUE)
tool.lm <- lm(Life~Angle*Speed+I(Angle^2)*I(Speed^2)+Angle:I(Speed^2)+I(Angle^2):Speed,tool) 
summary(tool.lm)

#response surface
tmp.angle <- seq(15,25,by=.1)
tmp.speed <- seq(125,175,by=.5)
tmp <- list(Angle=tmp.angle,Speed=tmp.speed) 
new <- expand.grid(tmp)
new$fit <- c(predict(tool.lm,new))

contourplot(fit~Angle*Speed,data=new,cuts=8,region=T,col.regions=gray(7:16/16))

```

**Blocking in factorial designs**
Run a factorial model within each block. Do not consider factor by block interactions.


*Example: Intensity*
Experiment to improve detecting targets on a radar scope. The two factors are the amount of background noise (“ground clutter”) on the scope (3 levels) and the type of filter (2 types) placed over the screen. A blocking factor is operator (availability and knowledge). Hence there are 3x2 treatment combinations.


```{r intensity}
intensity <- read.table("~/Desktop/Teaching/801/Lectures/Rcode/Data/intensity.txt",header=TRUE, colClasses=c("numeric",rep("factor",3)))

xyplot(Intensity~Ground|Operator,data=intensity,groups=Filter,
panel=function(x,y,...){
subs <- list(...)$subscripts
panel.xyplot(x,y,pch=c(1,19),...) 
panel.superpose(x,y,panel.groups="panel.lmline",lty=c(1,2),...)
},key=list(text=list(lab=as.character(1:2)),lines=list(lty=1:2,col=1), corner=c(1,.95),title="Filter Type",cex.title=.8),col=1)
#Note large inter-individual variation with respect to the response variable


# intensity.aov <- aov(Intensity~Ground*Filter+Error(Operator),intensity)
intensity.aov <- aov(Intensity~Operator+Ground*Filter,intensity) 
summary(intensity.aov)

```

The blocking factor SS is rather large compared to other main effects SS.


## Chapter 6: $2^K$ Factorial designs

*Example: Reactant yield*
The objective is to study how reactant concentration (15 or 25%) and the catalyst (1 or 2 pounds) impact the conversion (yield) in a chemical process, with three replicates.

```{r yield}
chem <- read.table("~/Desktop/Teaching/801/Lectures/Rcode/Data/yield.txt",header=T)

yield.sums <- aggregate(chem$yield,list(reactant=chem$reactant,catalyst=chem$catalyst),sum)
yield.sums

summary(aov(yield~reactant*catalyst, chem))

attach(chem)
reactant.num <- chem$reactant
levels(reactant.num) <- c(25,15)
reactant.num <- as.numeric(as.character(reactant.num)) 
catalyst.num <- chem$catalyst
levels(catalyst.num) <- c(2,1)
catalyst.num <- as.numeric(as.character(catalyst.num)) 
yield.lm <- lm(yield~reactant.num+catalyst.num, chem) 
yield.lm ## gives the coefficients of the LM

s3d <- scatterplot3d(reactant.num,catalyst.num,chem$yield,type="n", angle=135,scale.y=1,xlab="Reactant",ylab="Catalyst")
s3d$plane3d(yield.lm,lty.box="solid",col="darkgray")

tmp <- list(reactant.num=seq(15,25,by=.5),catalyst.num=seq(1,2,by=.1)) 
new.data <- expand.grid(tmp)
new.data$fit <- predict(yield.lm,new.data) 
contourplot(fit~reactant.num+catalyst.num,new.data,xlab="Reactant",ylab="Catalyst")

rm(chem)
```


For a first-order model, which includes only the main effects, the response curve is a plane.


*Example:  Filtration rate*

A chemical product is produced in a  pressure vessel. Four factors are thought to influence the filtration rate, namely temperature (A), pressure (B), concentration of formaldehyde (C), and stirring rate (D).

```{r filter}
# response
y<- c(45, 71, 48, 65, 68, 60, 80, 65, 43, 100, 45, 104, 75, 86, 70, 96)


#factor
A<-gl(n=2, k=1, length=16, labels=c(-1,1), ordered=F)
B<-gl(n=2, k=2, length=16, labels=c(-1,1), ordered=F)
C<-gl(n=2, k=4, length=16, labels=c(-1,1), ordered=F)
D<-gl(n=2, k=8, labels=c(-1,1), ordered=T)

filter.df<-data.frame(A=A, B=B, C=C, D=D, y=y)

filter.aov<-aov(y~A*B*C*D, filter.df)
filter.aov



#numeric
A<-rep(c(-1,1),  len = 16)
B<-rep(c(-1,1), each = 2, len = 16)
C<-rep(c(-1,1), each = 4, len = 16)
D<-rep(c(-1,1), each = 8, len = 16)

# show interactions
filterdata<-data.frame(y=y, A=A, B=B, C=C, D=D, AB=A*B, AC=A*C, AD=A*D, BC=B*C, BD=B*D, CD=C*D, ABC=A*B*C,ABD=A*B*D, ACD=A*C*D, BCD=B*C*D, ABCD=A*B*C*D)
filterdata

# regression model to get effects (all other is  NA!)
filter.lm<-lm(y~A*B*C*D, filterdata); summary(fit)


filter.effect<-as.matrix(summary(filter.lm)$coefficients[-1, 1:2])
factornames<-dimnames(summary(filter.lm)$coefficients[-1, 1:2])[[1]]
tmp<-data.frame(factors=factornames, estimate=filter.effect[,1], effect=2*filter.effect[,1])
tmp[order(tmp$effect),]
```


Potentially significant effects: A, AD, C, D, AC

```{r filterinteraction}
# check interaction
interaction.plot(filterdata$A, filterdata$C, filterdata$y,  fun=mean, col=1:6, ylab="mean of y", xlab="A")

interaction.plot(filter.df$A, filter.df$C, filter.df$y,  fun=mean, col=1:6, ylab="mean of y", xlab="A")

# ANOVA model involving only A, C, D and their interactions
fit.aov2<-aov(y~A*C*D, filter.df)
fit.aov2

# Regression model with A, C, D, AC, and AD
fit.lm2<-lm(y~A+C+D+AC+AD, filterdata)
summary(fit.lm2)

contrasts(filter.df$A)
contrasts(filter.df$C)

X<-model.matrix(fit.lm2); X
coef <-  solve(t(X)%*%X)%*%t(X)%*%filterdata$y
coef


# model diagnostics
plot(fit)

# Contourplot 
# Set D=1 (high level) to maximize response, see >fit 

xx<-as.numeric(fit$coef)

threevar <- function(A,C,D){ 
 #  xx[1] + xx[2]*A + xx[3]*C +  xx[4]*D +xx[5]*A*C + xx[6]*A*D}
  70.0625 + 10.812*A + 4.938*C +  7.313*D -9.063*A*C + 8.312*A*D}

mat <- outer( seq(-1, 1, by=0.1),  
              seq(-1, 1, by=0.1), 
              Vectorize( function(x,y) threevar(A=x, C=y, D=1) ) )

 contourplot(mat, xlab="A", ylab="C", row.values=seq(-1, 1, by=0.1), column.values=seq(-1, 1, by=0.1))

```


*Example:  Plasma*

This is a $2^3$ design used to develop a nitride etch process on a single-wafer plasma etching tool. The design factor are the gap [in cm]  between the electrodes, the gas flow ($C_2 F_6$ is used a the reactant gas), and the RF power [in W] applied to the cathode. Each factor is run at two levels, and the design is replicated twice.

```{r plasma}
plasma <- read.table("~/Desktop/Teaching/801/Lectures/Rcode/Data/plasma.txt",header=T)
plasma

# convert to long format with library(reshape2)
plasma.df<-reshape(plasma, varying=c("R1", "R2"),  idvar = "Run", v.names="etch", direction = "long" )

plasma.df.aov <- aov(etch~A*B*C, data=plasma.df)
summary(plasma.df.aov)



```

Factors A, C and their two-way interaction are significant.

The levels of the factors are as follows:

Factor  | Low (-1) | High (+1) |
---| ---| ---|
A (gap) | 0.80 | 1.20 |
B (flow) | 125 | 200 |
C (power) | 275 | 325 |


```{r plasmasurface}
plasma.num <- plasma.df
levels(plasma.num$A) <- c(0.8,1.2) 
levels(plasma.num$C) <- c(275,325)
plasma.num$A <- as.numeric(as.character(plasma.num$A)) 
plasma.num$C <- as.numeric(as.character(plasma.num$C)) 

plasma.num.lm <- lm(etch~A*C, plasma.num)
#require(scatterplot3d)
# s3d <- scatterplot3d(plasma.num$A,plasma.num$C, plasma.num$etch,type="n", angle=135,scale.y=1,xlab="Gap",ylab="Power") 
# doesn't work?
# s3d$plane3d(plasma.num.lm,lty.box="solid",col="darkgray")

tmp <- list(C=seq(275,325,by=1),A=seq(0.8,1.2,by=.1)) 
new.data <- expand.grid(tmp)
new.data$fit <- predict(plasma.num.lm, new.data) 

contourplot(fit~A+C,new.data,xlab="gap",ylab="power")

```



## Chapter 13: Random effects
 
Montogomery Example 13.1: A gauge is used to measure a critical dimension on a part. In order to investigate the sources of variability 20 parts have been selected from the prodcution process, and three operators measure each part twice.


```{r gauge}
#library(lme4)
part<-rep(c(1:20),each=6)
operator<-rep(1:3,each=2)
resp<-c(
21,20,20,20,19,21, 24,23,24,24,23,24, 20,21,19,21,20,22,
27,27,28,26,27,28, 19,18,19,18,18,21, 23,21,24,21,23,22,
22,21,22,24,22,20, 19,17,18,20,19,18, 24,23,25,23,24,24,
25,23,26,25,24,25, 21,20,20,20,21,20, 18,19,17,19,18,19,
23,25,25,25,25,25, 24,24,23,25,24,25, 29,30,30,28,31,30,
26,26,25,26,25,27, 20,20,19,20,20,20, 19,21,19,19,21,23,
25,26,25,24,25,25, 19,19,18,17,19,17)
gauge<-cbind(part,operator,resp);gauge<-data.frame(gauge)


#library(car)
scatterplot(resp~part|operator, data=gauge, boxplots=FALSE, smooth=TRUE, reg.line=FALSE)
abline(lm(gauge$resp~gauge$part), lwd=3, col="grey")
fixed.dummy<-lm(resp~part +factor(operator) -1, data=gauge)
summary(fixed.dummy)

# library(gplots)
plotmeans(resp~operator, gauge, main="Heterogeneity across operators")
plotmeans(resp~part, gauge, main="Heterogeneity across parts")

xyplot(gauge$resp~gauge$part|factor(gauge$operator), pch=16, col=4, xlab="part", ylab="response")

# both part and operator are random effects
gauge.lmer <- lmer(resp ~ (1|part) + (1|operator) , gauge)
summary(gauge.lmer)

# operator is a  fixed effect, part a random effect
gauge.lmer2 <- lmer(resp ~ operator+ (1|part)  , gauge)
summary(gauge.lmer2)

sig2part<-as.vector(VarCorr(gauge.lmer)$part)
sig2op<-as.vector(VarCorr(gauge.lmer)$operator)
sig.gauge<-summary(gauge.lmer)$sigma

plot(fitted(gauge.lmer), residuals(gauge.lmer))


# fitted values are not sample means! These are "shrinkage estimates" closer to the grand mean
means=sapply(split(resp, operator), mean)
operatorfit<-sapply(split(fitted(gauge.lmer), operator), mean)
data.frame(mean=means, fitted=operatorfit)

```

## Chapter 14: Nested designs
A company is interested in testing the uniformity of their film-coated pain tablets. A random sample of three batches were collected from each of their two blending sites. Five tablets were assayed from each batch.

```{r tablets}
site<-gl(n=2, k=3, length=30, labels=c(1,2))
batch<-gl(n=6, k=1, length=30, labels=1:6)
y<-c(5.03, 4.64, 5.1, 5.05, 5.46, 4.9,
     5.1, 4.73, 5.15, 4.96, 5.15, 4.95, 
     5.25, 4.82, 5.20, 5.12, 5.18, 4.86,
     4.98, 4.95, 5.08, 5.12, 5.18, 4.86,
     5.05, 5.06, 5.14, 5.11, 5.11, 5.07)
tablets<-data.frame(site=site, batch=batch, y=y)

boxplot(y~site)
boxplot(y~site:batch)

#library(lattice)
dotplot(y ~ batch | site)


#numeric factors
site<-as.numeric(rep(c(1,2), each=3, len=30))
batch<-as.numeric(rep(1:6, len=30))

#factorial without interaction
tablets.lm1 <- lm(y ~ site+batch)
summary(tablets.lm1)
anova(tablets.lm1)

# full factorial
tablets.lm2 <- lm(y ~ site*batch)
summary(tablets.lm2)
anova(tablets.lm2)


#random batch within site
tablets.lm3 <- lm(y ~ site + (site|batch))
summary(tablets.lm3)
anova(tablets.lm3)

#nested
tablets.lm4 <- lm(y ~ site/batch )
summary(tablets.lm4)
anova(tablets.lm4)

```

Hower the F values should be  MSA/MSB(A) for site (site is fixed)  and  MSB(A)/MSE for batch nested within site.

```{r tabletsmodel}
cat("F test statistic for site=", 0.02133/0.11338, "p value =", 1-pf(0.02133/0.11338, 1,4))

cat("F test statistic for batch nested within site=", 0.11338/0.01216, "p value =", 1-pf(0.11338/0.01216, 4, 24))

tablets.lmer<-lmer(y ~ site+ (1|batch:site) , tablets)
summary(tablets.lmer)
```

$\hat{sigma}^2_\beta= (MS_{B(A)}-MSE)/n = 0.02024$

$\hat{sigma}^2= MSE=0.01216$

So we can conclude that variation due to site can be ignored, but batch varies a lot and sites should be encouraged to reduce the batch-to-batch variation.



```{r tablets2}
TukeyHSD(aov(tablets.lm), "site:batch")

# library(lsmeans)
# same as TukeyHSD
lsmeans(tablets.lm, pairwise ~ site:batch)
# this gives unadjusted pvalues
lsmeans(tablets.lm, pairwise ~ site:batch, adjust="none")

```

## Chapter 15: Repeated measures 

In this exercise study  (UCLA wesbsite on epeated measures) people were randomly assigned to two types of diet (low-fat, regular), 
three types of exercises (at rest, walking, running). 
Pulse rate was measured at three different time points (1 min, 15 min, 30 min).

```{r exercise}
exer<-read.csv("~/Desktop/Teaching/801/Lectures/Rcode/Data/exer.csv",header=T)

#FYI transforming data from long to wide format
#library(reshape2)
exer.w <- reshape(exer,
   timevar=c( "time"),
  idvar = c("id","diet", "exertype"),
  direction = "wide")

head(exer.w)


#FYI transforming data from wide to long format
exer.l <- reshape(exer.w,
   timevar=c( "time"),
   varying=c( "pulse.1", "pulse.2", "pulse.3"),
   v.names = "pulse",
   times=1:3,  #this orders the data by ID
   new.row.names = 1:1000,
  direction = "long")

head(exer.l)

interaction.plot(exer$time,exer$diet,exer$pulse, 
ylim=c(90,110),lty=c(1,12),ylab="mean of pulse",
xlab="time")
```

Pulse rate of both diet types increase over time; for the regular diet the pulse rate increases more than for the low-fat diet.

```{r exercisemodel1}
exer$subject<-factor(exer$id, c(1:30))
diet.aov<-aov(pulse ~diet*time + Error(subject), exer)
summary(diet.aov)

diet.aov<-aov(pulse ~diet*time + Error(subject), exer)
summary(diet.aov)
```

Model 2: group variable=diet
```{r exercisemodel2}
interaction.plot(exer$time,exer$exertype,exer$pulse, 
ylim=c(80,130),lty=c(1,12),ylab="mean of pulse",xlab="time")

exertype.aov<-aov(pulse ~exertype*time + Error(subject), exer)
summary(exertype.aov)

```

Compare covariance structures
```{r covariance}
#compound symmetry
fit.cs<-gls(pulse ~ diet*time, data=exer, 
corr=corCompSymm(,form=~1|id))
summary(fit.cs)
intervals(fit.cs)
plot(fit.cs)

plot(fit.cs, id~resid(.), abline=0)

#autoregressive
fit.ar1<-gls(pulse ~ diet*time, data=longg, 
corr=corAR1(,form=~1|id))
summary(fit.ar1)

#Autoregressive with heterogeneous variances
fit.arh1<-gls(pulse ~ diet*time, data=longg, 
corr=corAR1(,form=~1|id), weight=varFixed(~time))
summary(fit.arh1)

#model comparisons
anova(fit.cs, fit.ar1, fit.arh1)
```


Model 3: time, diet, exertype
```{r exercisemodel3}
interaction.plot(exer$time[exer$diet==1], 
exer$exertype[exer$diet==1], exer$pulse[exer$diet==1],
ylim=c(80, 150), lty=c(1, 12, 8), ylab="mean of pulse", 
xlab="time")
title("Diet=1")

interaction.plot(exer$time[exer$diet==2], 
exer$exertype[exer$diet==2], exer$pulse[exer$diet==2],
ylim=c(80, 150), lty=c(1, 12, 8), ylab="mean of pulse", 
xlab="time")
title("Diet=2")

longa <- groupedData(pulse~exertype*diet*time|id, data=exer)
both.arh1 <- gls(pulse ~ exertype*diet*time, 
data=longg, corr=corAR1(, form= ~ 1 | id), 
weight=varFixed(~time),  )
summary(both.arh1)

#contrasts and interactions
longa <- groupedData(pulse~exertype*diet*time | id, data=exer)
longa$exertype.f <- factor(longa$exertype, c(1,2,3))
longa$diet.f <- factor(longa$diet, c(1,2))
longa$time.f <- factor(longa$time, c(1,2,3))

m <- matrix( c( c(-1/2, 1/2, 0), c(-1/3, -1/3, 2/3) ), ncol=2)
contrasts(longa$exertype.f) <- m
contrasts(longa$diet.f) <- c(-1/2, 1/2)

model.cs <- gls(pulse ~ exertype.f*diet.f , 
data=longa, corr=corCompSymm(, form= ~ 1 | id) )
summary(model.cs)


# other contrasts e.g. difference in mean pulse rate for runners in the two diet groups
longa$e1d12 <- -1/2*(longa$exertype==1 & longa$diet==1)
longa$e1d12[longa$exertype==1 & longa$diet==2] <- 1/2
longa$e2d12 <- -1/2*(longa$exertype==2 & longa$diet==1)
longa$e2d12[longa$exertype==2 & longa$diet==2] <- 1/2
longa$e3d12 <- -1/2*(longa$exertype==3 & longa$diet==1)
longa$e3d12[longa$exertype==3 & longa$diet==2] <- 1/2
modela.cs <- gls(pulse ~ exertype.f + e1d12+ e2d12+e3d12 , 
data=longa, corr=corCompSymm(, form= ~ 1 | id) )
```