anova_power.R

#
# Built in ANOVA
#
p16.7 <- read.table("CH16PR07.txt")
names(p16.7) <- c("improvement","rd","observation")
result <- lm(improvement~rd,data=p16.7)
anova(result)
#
result <- lm(improvement~factor(rd),data=p16.7)
anova(result)
p16.7$rd[p16.7$rd==1] <- "low"
p16.7$rd[p16.7$rd==2] <- "moderate"
p16.7$rd[p16.7$rd==3] <- "high"
result2 <- lm(improvement~rd,data=p16.7)
anova(result2)
#
# Using the "pwr" Package
#
install.packages("pwr")
library(pwr)
mu <- c(550,598,598,646)
sigma <- 80
(ef.size <- sqrt(sum(1/length(mu)*(mu-mean(mu))^2)/sigma^2))
pwr.anova.test(k=4,f=ef.size,sig.level=.05,power=.8)
pwr.anova.test(k=4,f=ef.size,sig.level=.05,n=15)
#
#  Simulation
#
mu <- c(550,598,598,646)
sigma <- 80
n <- c(15,15,15,15)
nt <- sum(n)
(y <- rnorm(nt,mean=rep(mu,n),sigma))
(x <- rep(1:length(mu),n))
boxplot(y~x)
anova(lm(y~factor(x)))
names(anova(lm(y~factor(x))))
anova(lm(y~factor(x)))$Pr
anova(lm(y~factor(x)))$Pr[1]

pvalues <- c()
for (i in 1:1000) {
  y <- rnorm(nt,mean=rep(mu,n),sigma)
  x <- rep(1:4,n)
  pvalues <- c(pvalues,anova(lm(y~factor(x)))$Pr[1])
}
pvalues
plot(pvalues)
abline(.05,0,col="red")
sum(pvalues<=.05)
sum(pvalues<=.05)/length(pvalues)

#
# Effect of different sigmas
# 

mu <- c(550,598,598,646)
sigma <- c(80,60,60,40)
n <- c(15,15,15,15)
nt <- sum(n)
pvalues <- c()
for (i in 1:1000) {
  y <- rnorm(nt,mean=rep(mu,n),sd=rep(sigma,n))
  x <- rep(1:4,n)
  pvalues <- c(pvalues,anova(lm(y~factor(x)))$Pr[1])
}
plot(pvalues)
abline(.05,0,col="red")
sum(pvalues<=.05)
sum(pvalues<=.05)/length(pvalues)