Simulation for Valeria ANOVA

valisim.r

#Vali-simulation
rm(list = ls())
library(car)
# 2x2 simulation with 1 group all == 0

# n per group
n <- 15

meanWTct <- 1 
varWTct <- 0 # the all - equal group

meanWTELS <- 1 # no effect
varWTELS <- 1

meanKOct <- 1
varKOct <- 1

meanKOELS <- 1
varKOELS <- 1

set.seed(3665364)

k <- 10000
results <- matrix(numeric(k*3),ncol=3)
colnames(results) <- c("Main effect gen", "Main effect con", "Interaction")

for (i in 1:k){
  df <- data.frame(val = c(rnorm(n,meanWTct,sqrt(varWTct)),
                           rnorm(n,meanWTELS,sqrt(varWTELS)),
                           rnorm(n,meanKOct,sqrt(varKOct)),
                           rnorm(n,meanKOELS,sqrt(varKOELS))), 
                   gen = factor(c(rep("WT", 2*n), 
                                  rep("KO", 2*n))),
                   con = factor(c(rep("ctrl",n),
                                  rep("ELS",n),
                                  rep("ctrl",n),
                                  rep("ELS",n))))
  
  
  
  obj <- lm(val~gen*con, data = df)
  res <- Anova(obj, type = "III")
  results[i,] <- res$`Pr(>F)`[-c(1,5)]
}

sum(results<0.05)/length(results)