pdf("kuehl_02_03.pdf")

t <- 3
tau1 <- 0.5
tau2 <- -0.2
tau3 <- -0.3
sigma2 <- 0.22

sumtau2 <- (tau1^2) + (tau2^2) + (tau3^2)

r_power <- matrix(rep(0,6*29),ncol=6)

for (r in 2:30) {
trtdf <- t-1            # Treatment df
errdf <- t*(r-1)        # Error df
f_alpha <- qf(.95,trtdf,errdf)     # Critical Value for F-test
lambda <- r*sumtau2/sigma2           # Non-centrality parameter
power <- 1-pf(f_alpha,trtdf,errdf,lambda)
r_power[r-1,] <- cbind(r,trtdf,errdf,f_alpha,lambda,power)
}

power.names <- c("r","trtdf","errdf","f_0.05","lambda","power")
colnames(r_power) <- power.names

r_power

dev.off()