Power calculations

power_calcs.R

library(tidyverse)
library(lme4)

possible.ns <- 400
courseworks <- 8
prob_inclusion <- 0.1
powers <- rep(NA, length(possible.ns))           # Empty object to collect simulation estimates
alpha <- 0.05                                    # Standard significance level
sims <- 500

for (j in 1:length(possible.ns)){
  N <- possible.ns[j]                              # Pick the jth value for N
  
  significant_experiments <- rep(NA, sims) # Empty object to count significant experiments
  estimate <- rep(NA, sims)
  #### Inner loop to conduct experiments "sims" times over for each N ####
  for (i in 1:sims){
      random_effects <- rnorm(n = N, mean = 0, sd = 5)
      Y0 <- rnorm(n = N * courseworks, mean = 60, sd = 11.5)
      Y0 <- Y0 + rep(random_effects, times = courseworks)
      tau <- 3 ## treatment effect
      Y1 <- Y0 + tau
      is_chosen <- rbinom(n = N, size = 1, prob = prob_inclusion)
      ## Rep this for multiple courseworks
      is_chosen <- rep(is_chosen, times = courseworks)
      Y.sim <- Y1*is_chosen + Y0*(1-is_chosen)
      ## Subject ID
      id <- factor(rep(1:N, times = courseworks))
      fit <- lmer(Y.sim ~ is_chosen + (1|id))
      t_value <- summary(fit)$coefficients["is_chosen", "t value"]
      p_value <- dt(abs(t_value), df = Inf)
      significant_experiments[i] <- (p_value <= alpha) # Determine significance according to p <= 0.05
      estimate[i] <- summary(fit)$coefficients["is_chosen", 1]
  } 
}

mean(significant_experiments)
mean(estimate)