simulate_power.R

library(tidyverse)

logit <- function(p) log(p/(1-p))

simulate_trial<-function(N, effect_size){  
  catheter_diameter = sample(c(-1,0,1), replace = T, size = N)
  vaso_band = rbinom(N, 1, 0.5)
  
  X = model.matrix(~catheter_diameter*vaso_band)
  
  # baseline effect
  # corresponding to 5/6 catheter size
  # from sam's email
  base_line_effect = logit(0.011)
  
  # effect of increasing diameter
  beta_diameter = (logit(0.048) - base_line_effect)/2
  
  # effect vasoband
  # from sam's email
  effect_vaso = logit(0.01) - base_line_effect
  betas = c(base_line_effect, beta_diameter, effect_vaso, effect_size )
  
  # linear predictor
  eta = X%*%betas
  p = 1/(1+exp(-eta))
  y = rbinom(length(eta), 1, p=p)
  
  # Now model
  model = glm(y ~ catheter_diameter*vaso_band, family = binomial(link = 'logit'))
  
  # p value
  p = as.numeric(broom::tidy(model)[4, 'p.value'])
  
  return(p<0.05)
}