##############################################
###  Lec2.r   STA 4211  Script file with code for lecture notes pp. 106 - 162
###                                Needs dataset corn-yield.txt
##############################################
##
##
## Corn yield data
yield <- read.table("corn-yield.txt",header=TRUE)
##  This data file only contains the response variables.
##   Next create the group variable from scratch.
yield$fertilizer <- factor(rep(c("high","low"),c(10,10)), levels=c("low","high"))
yield$manure <- factor(rep(rep(c("high","low"),c(5,5)),2), levels=c("low","high"))
##  Check that the data setup is correct
yield[1:3,]
##
## One way to get summary statistics is with tapply()
##   Use round() to print summary stats to just one or two signif. places
tapply(yield$yield,   list(yield$manure,yield$fertilizer),   mean)
tapply(yield$yield,list(yield$manure,yield$fertilizer),sd)
## Function tapply() is similar to aggregate().
##  You decide which you prefer:  tapply(), or aggregate()
aggregate(yield$yield,list(yield$manure,yield$fertilizer),mean)
aggregate(yield$yield,list(yield$manure,yield$fertilizer),sd)
## First do dotplots of yield for each of the four trt groups
stripchart(yield ~ fertilizer*manure, data=yield)
## Comparative boxplots for each factor separately
##   Note: Boxplots for manure include both levels of fertilizer
##      and boxplots for fertilizer include both levels of manure.
plot(yield ~ manure*fertilizer, data=yield)
## Factor level means
plot.design(yield ~ manure*fertilizer, data=yield)
## Interaction plot, both ways.  X- and Y- labels need shortening.
interaction.plot(yield$manure,yield$fertilizer,response=yield$yield)
interaction.plot(yield$fertilizer,yield$manure,response=yield$yield)
####
##
## Fit two-way ANOVA model with interaction
##   Get ANOVA table
##   and look at residual plots for a check of assumptions
m1 <- aov(yield ~ manure*fertilizer,data=yield)
anova(m1)
plot(fitted(m1), rstandard(m1))
qqnorm(rstandard(m1))
##
## Fit the additive model, do the residual plots.
m2 <- aov(yield ~ manure + fertilizer, data=yield)
anova(m2)
plot(fitted(m2), rstandard(m2))
qqnorm(rstandard(m2))
##  Estimate main effects of fertilizer (done here) and manure (not
##  included).
fert.h <- mean(c(15.1,13.92)); fert.l <- mean(c(14.0,11.26))
fert.main.eff <- fert.h - fert.l
s <- sqrt(2.7907)
se <- sqrt(2.7907*(2/10))
cc <- qt(.975,df=17)
allow <- cc*se
ci.l <- fert.main.eff - allow; ci.u <- fert.main.eff + allow
ci.l; ci.u
## Clean up from corn yield data.
rm(allow,cc,ci.l,ci.u,fert.h,fert.l,fert.main.eff,s,se,yield,m1,m2)
##
#### Poisons data
library("boot")
data(poisons)
attach(poisons)
rate <- 1/poisons$time
poisons.rate <-data.frame(rate=rate,poisons)
##
detach(poisons) # we are finished with the original data frame
                                        #  so detach it
attach(poisons.rate) # for ease of doing the interaction plot
##
cell.means <-
round(tapply(poisons.rate$rate,list(poisons.rate$poison,poisons.rate$treat),
mean),digits=2)
cell.means
cell.sds <-
round(tapply(poisons.rate$rate,list(poisons.rate$poison,poisons.rate$treat),FUN=sd),digits=2)
cell.sds
plot.design(poisons.rate,y=rate)
interaction.plot(treat,poison,response=rate)
interaction.plot(poison,treat,response=rate)
rate.aov <- aov(rate ~ poison*treat,data=poisons.rate)
anova(rate.aov)
##  Chart on p. 160
##  Get row means (Poison)
P.means <-round(aggregate(rate,list(poison),FUN=mean)$x,digits=2); P.means#row
                                                                          #means,
                                                                          #bn=16
##  Get column means (Treatment)
T.means <- round(aggregate(rate,list(treat),FUN=mean)$x,digits=2); T.means
##  Get grand mean
grandmu <- mean(rate); round(grandmu,digits=2)
####
####
## Clean up
detach(poisons)
rm(P.means,T.means,cell.means,cell.sds,grandmu,poisons,poisons.rate,rate,rate.aov)
####