Power simulations with randomized enrollment month and TIR as outcome.

full_dt = fread("~/Downloads/CGM/deid_daily_timeseries.csv")

# Aggregate to month
month_dt = full_dt[months_since_onset>0, .(
  p_worn = mean(fcoalesce(p_worn,0)),
  TIR = mean(TIR, na.rm=T),
  epic_read = max(fcoalesce(epic_read,0L)),
  public_insurance = max(public_insurance)
  ), by=c('pop','pat_id','months_since_onset')]


# Define when patients are starting program
month_dt[,start_month := 1 + 3*(pat_id %% 4)]
month_dt[pat_id>=150, start_month:= 1E6]

month_dt[,started := months_since_onset>=start_month]
month_dt[,time_since_started := months_since_onset - start_month]

month_dt[, uniqueN(pat_id)]

## [1] 198

Define treatment effect by weeks_in_program.

Will assume 5pp treatment effect until 6 months and then 10pp from 6 months on.

treatment_effects = data.table(time_since_started = seq(0,18))
treatment_effects[, treatment_effect := 0.00]
treatment_effects[time_since_started>=0 & time_since_started<6, treatment_effect := 0.05]
treatment_effects[time_since_started>=6, treatment_effect := 0.10]

# Apply treatment effects
month_dt_w_te = merge(month_dt, treatment_effects, by='time_since_started', all.x=T)
month_dt_w_te[is.na(treatment_effect), treatment_effect:=0]
month_dt_w_te[, TIR_adj := pmin(1, TIR + treatment_effect)]

Try estimating treatment effects

# Find patients with full time series

month_dt_w_te[months_since_onset<=18 & months_since_onset>=1 & !is.na(TIR), unique_obs := uniqueN(months_since_onset), by='pat_id']
filtered_dt = copy(month_dt_w_te[unique_obs==18])
filtered_dt[, uniqueN(pat_id)]

## [1] 94

eventPlotResults <- staggered(df = filtered_dt,
                              i = "pat_id",
                              t = "months_since_onset",
                              g = "start_month",
                              y = "TIR_adj",
                              estimand = "eventstudy",
                              eventTime = -1:12)

## Warning in compute_se_Thetahat_beta(beta = beta, Ybar_g_list, A_theta_list, :
## var_conservative is less than adjustmentFactor

ggplot(eventPlotResults, aes(x=eventTime, y=estimate, ymin=estimate-2*se, ymax=estimate+2*se)) + geom_pointrange() + theme_bw() + geom_vline(xintercept = -0.5, linetype="dashed", color="red")  + geom_vline(xintercept = 5.5, linetype="dashed", color="blue") + annotate("segment", x = 0, xend = 5, y = 0.05, yend = 0.05, colour = "red")  + annotate("segment", x = 6, xend = 12, y = 0.1, yend = 0.1, colour = "blue")

Try estiamting a single effect across months

staggered(df = filtered_dt,
          i = "pat_id",
          t = "months_since_onset",
          g = "start_month",
          y = "TIR_adj",
          estimand = "simple")

##     estimate         se  se_neyman
## 1 0.09584492 0.02080486 0.02080486

Pretty low SE so can maybe be used on our full population to get overall effect (not estimated for each month in program).