set.seed(139)
k <- 1000
m1 <- 10
m2 <- 10
m3 <- 10
power <- c(rep(0,5))
pvalue <- c(rep(0,k))
for(i in 1:k){
res1 <- -3 + rnorm(m1,0,1)
res2 <- -2 + rnorm(m2,0,1)
res3 <- -1.5 + rnorm(m3,0,1)
res <- c(res1,res2,res3)
fac <- c(rep("A",m1),rep("B",m2),rep("C",m3))
mod <- aov(res~fac)
mod.summary <- summary(mod)
pvalue[i] <- mod.summary[[1]][,5][[1]]
}
power[1] <- length(pvalue[which(pvalue < 0.05)])/k
set.seed(139)
pvalue <- c(rep(0,k))
for(i in 1:k){
res1 <- -rexp(m1,1/3)
res2 <- -2 + rnorm(m2,0,1)
res3 <- -1.5 + rnorm(m3,0,1)
res <- c(res1,res2,res3)
fac <- c(rep("A",m1),rep("B",m2),rep("C",m3))
mod <- aov(res~fac)
mod.summary <- summary(mod)
pvalue[i] <- mod.summary[[1]][,5][[1]]
}
power[2] <- length(pvalue[which(pvalue < 0.05)])/k
set.seed(139)
pvalue <- c(rep(0,k))
for(i in 1:k){
res1 <- -3 + rnorm(m1,0,1)
res2 <- -2 + rnorm(m2,0,1)
res3 <- -rexp(m3,2/3)
res <- c(res1,res2,res3)
fac <- c(rep("A",m1),rep("B",m2),rep("C",m3))
mod <- aov(res~fac)
mod.summary <- summary(mod)
pvalue[i] <- mod.summary[[1]][,5][[1]]
}
power[3] <- length(pvalue[which(pvalue < 0.05)])/k
set.seed(139)
pvalue <- c(rep(0,k))
for(i in 1:k){
res1 <- -rexp(m1,1/3)
res2 <- -2 + rnorm(m2,0,1)
res3 <- -rexp(m3,2/3)
res <- c(res1,res2,res3)
fac <- c(rep("A",m1),rep("B",m2),rep("C",m3))
mod <- aov(res~fac)
mod.summary <- summary(mod)
pvalue[i] <- mod.summary[[1]][,5][[1]]
}
power[4] <- length(pvalue[which(pvalue < 0.05)])/k
set.seed(139)
pvalue <- c(rep(0,k))
for(i in 1:k){
res1 <- -rexp(m1,1/3)
res3 <- -rexp(m3,1/2)
res3 <- -rexp(m3,2/3)
res <- c(res1,res2,res3)
fac <- c(rep("A",m1),rep("B",m2),rep("C",m3))
mod <- aov(res~fac)
mod.summary <- summary(mod)
pvalue[i] <- mod.summary[[1]][,5][[1]]
}
power[5] <- length(pvalue[which(pvalue < 0.05)])/k
data.frame(Power=power, Skew.no.temp=c("No","No","No","No","Yes"), Skew.pos.treat=c("No","Yes","No","Yes","Yes"), Skew.neg.treat=c("No","No","Yes","Yes","Yes"))
set.seed(139)
res1 <- -rexp(m1,1/3)
res3 <- -rexp(m3,1/2)
res3 <- -rexp(m3,2/3)
res <- c(res1,res2,res3)
fac <- c(rep("Heat",m1),rep("No.temp",m2),rep("Cool",m3))
boxplot(res~fac)