library(data.table)
library(ggplot2)
library(gridExtra)
library(boot)
library(knitr)
library(printr)
library(doParallel)
library(parallel)

R = 200

Load TIDE data

sim_dt = fread('~/Google Drive/My Drive/Stanford/CGM/sim_dt.csv')

eval_dt = sim_dt[in_first_year_since_onset==0 & !is.na(dTIR_2w) & !is.na(TIR_lag)]

eval_dt[,.(uniqueN(pat_id), .N, mean(eligible_for_review), sum(tide_message_sent))]

V1	N	V3	V4
78	2740	0.2481752	236

Power Simulations with Constant treatment effects

# Define fn to return coefficient and pval given a constant treatment effect and sample size

boostrap_constant_eval = function(seed, d, t, n) {
  
  # bootstrap sample (n) observations
  set.seed(seed)
  boot_dt = d[sample.int(nrow(d), n, replace = T)]
  
  # randomize 50% of sample as treated
  boot_dt[, treated := rbinom(.N, 1, 0.5)]
  
  # apply constant treatment effect (t) if treated
  boot_dt[, outcome := dTIR_2w + treated*t]
  
  # fit OLS to estimate effect and p-value
  m = lm(outcome ~ treated, data=boot_dt)
  
  # return coefficient and p-value
  return(data.table(
    n = n,
    t = t,
    coef = summary(m)$coefficients["treated", "Estimate"],
    pval = summary(m)$coefficients["treated", "Pr(>|t|)"]))
}


# Define grid of effect sizes and sample sizes to bootstrap over
sim_grid_dt = as.data.table(expand.grid(
  tau = seq(0,0.02,0.0025),
  size = c(1200, 2400, 52*130, 52*200)
))

cl <- parallel::makeCluster(detectCores())
clusterExport(cl, varlist=c('boostrap_constant_eval','eval_dt', 'sim_grid_dt'))

doParallel::registerDoParallel(cl)

out_dt = rbindlist(lapply(1:nrow(sim_grid_dt), function(i) {
  
  foreach::foreach(s = 1:R, .combine = rbind) %dopar% {
      library(data.table)
      boostrap_constant_eval(s, eval_dt, sim_grid_dt$tau[i], sim_grid_dt$size[i])
  }
  
}))

parallel::stopCluster(cl)

# Aggregate
agg_out_dt = out_dt[, .(
  mean_coef = mean(coef),
  median_coef = quantile(coef, 0.5),
  power = mean(pval<0.05),
  power_me = 1.96*sqrt(mean(pval<0.05)*mean(pval>=0.05)/.N)
  ), by=c('n','t')]

ggplot(agg_out_dt, aes(x=t, y=power, color=factor(n), ymin=power-power_me, ymax=power+power_me)) + geom_line() + geom_point() + geom_errorbar(width=0.0005, position=position_dodge(0.0005)) +
  theme_bw() + 
  xlab("Effect size") + 
  # random targeting estimates from paper
  annotate("point", x=0.00912, color="red", y=.282, size=5, shape=10) +
  annotate("point", x=.00889, color="red", y=0.264, size=5, shape=10) +
  annotate("point", x=0.00906, color="green", y=0.468, size=5, shape=10) +
  annotate("point", x=0.00919, color="green", y=0.482, size=5, shape=10) +
  annotate("point", x=0.00919, color="blue", y=0.902, size=5, shape=10) +
  annotate("point", x=0.00923, color="purple", y=0.976, size=5, shape=10) +
  # TIR_lag targeting estimates from paper
  annotate("point", x=.01018, color="red", y=0.276, size=5, shape=11) +
  annotate("point", x=0.01517, color="green", y=0.780, size=5, shape=11) +
  annotate("point", x=0.01011, color="green", y=0.466, size=5, shape=11) +
  # smarter targeting estimates from paper
  annotate("point", x=0.01935, color="red", y=.754, size=5, shape=12) +
  annotate("point", x=0.01734, color="green", y=0.928, size=5, shape=12)