# ANOVA 2010 EX4-Solution


#1)
## ## Read in the data.
feed <- read.table(file="http://stat.ethz.ch/Teaching/Datasets/feed.txt",header=TRUE)
                         
feed$Feeding <- as.factor(feed$Feeding)

## 1a)

modF <- aov(Final~Feeding,data=feed)
summary(modF)

TukeyHSD(modF,"Feeding", conf.level=0.95)

## 1b)

modF2 <- aov((Final-Initial)~Feeding,data=feed)
summary(modF2)

## 1c)

modF3 <- aov(Final~Feeding+Initial,data=feed)
summary(modF3)


par(mfrow=c(2,2))
boxplot(feed[,2]~feed$Feeding,ylab="Final",xlab="TR")
boxplot(feed[,2]-feed[,1]~feed$Feeding,ylab="Final-Initial",xlab="TR")
boxplot(feed$Initial~feed$Feeding,ylab="Initial",xlab="TR")
plot(feed[,1],feed[,2],xlab="Initial",ylab="Final")
lines(190:240,190:240)



par(mfrow=c(2,2))
plot(feed[,1],feed[,2]/feed[,1],xlab="initial",ylab="Final/Initial")
lines(190:240,190:240*0+1)
plot(feed[1:9,1],feed[1:9,2]/feed[1:9,1],xlab="initial",ylab="Final/Initial",main="Food 1")
lines(sort(feed[1:9,1]),feed[1:9,2]*0+1)
plot(feed[10:21,1],feed[10:21,2]/feed[10:21,1],xlab="initial",ylab="Final/Initial",main="Food 2")
lines(sort(feed[10:21,1]),feed[10:21,2]*0+1)
plot(feed[22:32,1],feed[22:32,2]/feed[22:32,1],xlab="initial",ylab="Final/Initial",main="Food 3")
lines(sort(feed[22:32,1]),feed[22:32,1]*0+1)


# 2)
soft<-read.table("http://stat.ethz.ch/Teaching/Datasets/softdrinkANOVA.txt",header=TRUE)
soft <- soft[1:16*2-1,]

soft$sugar <- as.factor(soft$sugar)
soft$soda <- as.factor(soft$soda)
soft$water <- as.factor(soft$water)
soft$temp <- as.factor(soft$temp)
str(soft)

sg <- as.numeric(soft[,2])*2-3
sd <- as.numeric(soft[,3])*2-3
wt <- as.numeric(soft[,4])*2-3
tp <- as.numeric(soft[,5])*2-3
BLOCK <- (sg*sd*wt*tp+3)/2

softBL <- data.frame(soft,BLOCK=BLOCK)
softBL$BLOCK <- as.factor(softBL$BLOCK)
str(softBL)

sB.fit <- aov(score~sugar+soda+water+temp+BLOCK,data=softBL)
summary(sB.fit)

sB.2k <- aov(score~sugar*soda*water*temp*BLOCK,data=softBL)
summary(sB.2k)

sB.2k$coef


# 3)
a <- c(-1,1,-1,1,-1,1,-1,1)
b <- c(-1,-1,1,1,-1,-1,1,1)
c <- c(-1,-1,-1,-1,1,1,1,1)
d <- a*b
e <- a*d


# 4)
v <- c(24.50,22.05, 24.52, 25.00, 25.68, 24.51, 24.68, 24.23)
a <- c(-1,1,-1,1,-1,1,-1,1)
b <- c(-1,-1,1,1,-1,-1,1,1)
c <- c(-1,-1,-1,-1,1,1,1,1)
d <- a*b*c
dat <- data.frame(Y=v,A=a,B=b,C=c,D=d)
dat$A <- as.factor(dat$A)
dat$B <- as.factor(dat$B)
dat$C <- as.factor(dat$C)
dat$D <- as.factor(dat$D)

str(dat)

# 4c)
EA <- sum(a*v/4)
EB <- sum(b*v/4)
EC <- sum(c*v/4)
ED <- sum(d*v/4)
EAB <- sum(a*b*v/4)
EAC <- sum(a*c*v/4)
EAD <- sum(a*d*v/4)

z <- data.frame(Effect=c("A","B","C","D","AB;CD","AC;BD","AD;BC"),Value=c(EA,EB,EC,ED,EAB,EAC,EAD))

# 4d)
abs(z$Value)>0.35

# 4e)
z$Value[1:4]>0