## R-code for HW #4
## AHW3
## Feb 27, 2008
rm(list = ls())

y <- matrix(c(100, 76, 108, 
	96, 80, 100, 
	92, 75, 96, 
	96, 84, 98, 
	92, 82, 100), 5,3, byrow=T)

t <- 3
r <- 5
N <- r*t		

## ANOVA table
y.batt <- c(y[,1], y[,2], y[,3])
x.batt <- rep(1:3, each=5)
anova(lm(y.batt~as.factor(x.batt)))
		
## randominization test
F.obs <- anova(lm(y.batt~as.factor(x.batt)))[1,4]
F.null <- real()
for(nsim in 1:5000){
x.sim <- sample(x.batt)
F.null[nsim] <- anova(lm(y.batt~as.factor(x.sim)))[1,4]
}
p.value <- mean(F.null>=F.obs)
p.value

## population test
1-pf(38.338, 3-1,15-3)

## 95% CIs
mu.i <- tapply(y.batt, x.batt, mean)
se.mu.i <- rep(sqrt(15.60/5),3)
CI <-  round(cbind(mu.i - se.mu.i*qt(1-0.05/2, N-3), mu.i + se.mu.i*qt(1-0.05/2, N-3)),3)
CI

## comparisons
k1 <- c(1,-1,0)
k2 <- c(1,0,-1)
k3 <- c(0,1, -1)

C.hat.1 <- k1%*%mu.i
C.hat.2 <- k2%*%mu.i
C.hat.3 <- k3%*%mu.i

SE.C.hat.1 <- sqrt(mse*sum(k1^2)/5)
SE.C.hat.2 <- sqrt(mse*sum(k2^2)/5)
SE.C.hat.3 <- sqrt(mse*sum(k3^2)/5)

C.hat <- c(C.hat.1, C.hat.2, C.hat.3)
SE.C.hat <- c(SE.C.hat.1,SE.C.hat.2,SE.C.hat.3)
CI.Bon <-  round(cbind(C.hat - SE.C.hat*qt(1-(0.05/3)/2, N-3), C.hat + SE.C.hat*qt(1-(0.05/3)/2, N-3)),3)
CI.Bon

## residual plots and other diagnostics
mod <- lm(y.batt~factor(x.batt))
pdf("diag.pdf",family="Times",height=3,width=6)
#par(mar=c(3,3,1,1))
par(mfrow=c(1,2))
qqnorm(mod$res); qqline(mod$res)
plot(mod$fit, mod$res)
dev.off()

max(tapply(y.batt, x.batt, var))/min(tapply(y.batt, x.batt, var))

## power 
pow <- function(r){
alpha <- 0.05
sigma.sq <- 20
t <- 3
N <- r*t
lambda <- (r*t*9)/sigma.sq
pow <-1-pf(qf(1-alpha, t-1, N-t), t-1, N-t, ncp=lambda)
return(pow)
}

pow(5)
pow(7)
pow(8)
pow(9)