##Code used for the review for Exp Design final.


## data-set

y <- c(0.5, 0.634, 0.487, 0.329, 0.512, 0.535, 0.675, 0.520, 0.435, 0.540, 0.513, 0.595, 0.488, 0.400, 0.510)
truck <- rep(1:5, 3)
oil <- rep(1:3, each=5) 

data <- data.frame(y, truck, oil)


## two sample t.test  vs one factor with two levels
y.oil1 <- y[oil==1]
y.oil2 <- y[oil==2]
n.1 <- length(y.oil1)
n.2 <- length(y.oil2)


y.oil1.bar <- mean(y.oil1)
y.oil2.bar <- mean(y.oil2)

var.oil1 <- var(y.oil1)
var.oil2 <- var(y.oil2)

s.sq.pooled <- (n.1/(n.1+n.2))*var.oil1 + (n.2/(n.1+n.2))*var.oil2

t.obs <- (y.oil1.bar-y.oil2.bar)/sqrt(s.sq.pooled*(1/n.1 + 1/n.2))

p.value <- 2*(1-pt(abs(t.obs), n.1+n.2-2))

t.test(y.oil1, y.oil2, var.equal=T)

## one factor
y.oil <- c(y.oil1, y.oil2)
x <- rep(1:2, each=c(n.1, n.2))
anova(lm(y.oil~as.factor(x)))

abs(t.obs)^2

## one factor with block
anova(lm(y~as.factor(oil) + as.factor(truck)))
anova(lm(y~as.factor(truck) + as.factor(oil)))

## two factor model
## data from Kirk chapter 9 question 5.

y <- c(14, 0, -6, 9, 10, -8, -2, -7, 2, -4, -10, 5, -2, 8, -18, 1, 6, 4, -14, -3)
shock <- rep(c("m", "m", "s", "s"), 5)
race <- rep(c("b", "w", "b", "w"),5)

options(contrasts=c("contr.sum", "contr.poly"))
lm.full <- lm(y~ as.factor(shock) + as.factor(race) + as.factor(shock):as.factor(race))
anova(lm.full)

par(mfrow=c(1,2))
plot(jitter(lm.full$fit), jitter(lm.full$res), pch=16, col="blue")
qqnorm(lm.full$res,pch=16, col="blue"); qqline(lm.full$res)

## 
summary(lm.full)


## lets get a confidence interval for the interaction: (mu11-mu12) - (mu21-mu22)

cell.means <- tapply(y, as.factor(shock):as.factor(race), mean)
c.hat <- (cell.means[1]-cell.means[2])-  (cell.means[3]-cell.means[4])
se.c.hat <- sqrt(4*40/5)

t.obs <- c.hat/se.c.hat
p.value <- 1-pt(t.obs, 16)
p.value