Type1ErrorSimulation

gistfile1.txt

sims <- replicate(1000,{ 
  X <- as.data.frame(matrix(rnorm(200*25), ncol = 25, nrow = 200))
  colnames(X)[1] <- "Y"
  
  has.cor <- sapply(2:ncol(X), function(i){cor.test(X$Y, X[,i])$p.value < 0.2})
  lm.out.1 <- lm(Y~., data = X)
  lm.out.2 <- lm(Y~., data = X[,c(1, which(has.cor))])
  
  ftest <- function(lmrslt){
    fs <- summary(lmrslt)$fstatistic
    pf(fs[1], fs[2], fs[3], lower.tail = FALSE)
  }
  
  pv.1 <- ftest(lm.out.1)
  pv.2 <- ftest(lm.out.2)
  c(pv.1, pv.2)
})

# Type 1 error rate for full model
mean(sims[1,] < 0.05)

# Type 1 error rate for model with selection
mean(sims[2,] < 0.05)