LRS_question_null_simulations

n <- 1e3
t <- 1e3

covar_effect <- 1
trt_effect <- 0 #0.1

null_no_covar <- alt_no_covar <- null_yes_covar <- alt_yes_covar <- rep(NA, t)
for (i in 1:t) {
  
  covar <- rnorm(n = n)
  snp <- rnorm(n = n)
  y <- rnorm(n = n, mean = covar*covar_effect + snp*trt_effect)
  
  null_no_covar[i] <- logLik(object = lm(formula = y ~ 1))
  alt_no_covar[i] <- logLik(lm(formula = y ~ snp))
  null_yes_covar[i] <- logLik(lm(formula = y ~ covar))
  alt_yes_covar[i] <- logLik(lm(formula = y ~ covar + snp))
}

my_line <- function(x, y, ...){
  points(x, y, ...)
  abline(a = 0, b = 1,...)
}

no_covar_LRS <- 2*(alt_no_covar - null_no_covar)
yes_covar_LRS <- 2*(alt_yes_covar - null_yes_covar)

pairs(x = data.frame(no_covar_LRS = sort(pchisq(q = no_covar_LRS, df = 1)),
                     yes_covar_LRS = sort(pchisq(q = yes_covar_LRS, df = 1)),
                     unif = seq(from = 0.5/t, to = 1 - 0.5/t, length.out = t)),
      panel = my_line,
      upper.panel = NULL)