sequential_trials_novavax.Rmd

---
title: "Sequential trials"
output: html_document
date: '2022-12-08'
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(rjags)
library(HDInterval)
library(dplyr, warn.conflicts = FALSE)
options(dplyr.summarise.inform = FALSE)

library(patchwork)
library(stringr)
library(pbapply)
library(parallel)
library(multidplyr)

library(plotly)

library(ggplot2)
library(extraDistr)
source('./R/prior.data.func.R')
source('./R/call_jags.R')
source('./R/sim.data.R')
source('./R/jags_basic_pois.R')
source('./R/jags_modified_pois_commensurate_gamma.R') #NOTE, this is experimental version
source('./R/jags_modified_pois_commensurate_uniform.R') #NOTE, this is experimental version
source('./R/jags_modified_pois_commensurate_half_cauchy.R') #NOTE, this is experimental version

source('./R/caterplot.func.R')
source('./R/prior_post_compare.R')
source('./R/extract_parm.R')

source('./R/run_all_models.R')

N.sim=500


```

Prior for VE proposed <https://hbiostat.org/proj/covid19/bayesplan.html#operating-characteristics> (low probability of extreme effects, symmetric around 0) A skeptical prior as defined by Frank harrell has low probability of extremes (e.g., 5% prior probability that the effect is \>95% or 5% that it is the inverse )

```{r}

# p(RR<0.05) = P(RR>1/0.05) 

#IRR=0.05 = 95%VE 
#log(0.05) log(IRR)

#IRR=1.95   #-95%VE




hist(rnorm(10000, mean=0, sd=1.821))

mean(rnorm(10000, mean=0, sd=1.821) <log(0.05) ) #~5% of people have a VE of >95%
mean(rnorm(10000, mean=0, sd=1.821) >log(1/0.05) ) # ~5% of people have a VE of >-95%


prior.sd.logrr <- (abs(log(1-44.4/100)  -  log(1-19.6/100)))/1.96
prior.mean.logrr <- log(1-44.4/100)

hist(rnorm(10000, mean=prior.mean.logrr, sd=prior.sd.logrr))  #informative prior

```

For commensurate prior, need to set hyperior for variance as well. If we were to use the historical data fully, we would want the hyperprior for posterior variance for commensurate to approach 0. If it doesn't use it at all, we would want variance to be sufficiently large to cover negative and neutral effects

```{r}



 hist(rgamma(1000,0.01, 0.01)) # gamma hyperprior from psborrow
 mean(rgamma(1000,0.01, 0.01)) # gamma hyperprior from psborrow

  hist(rgamma(1000,0.1, 0.1)) # gamma hyperprior from psborrow
 mean(rgamma(1000,0.1, 0.1)) # gamma hyperprior from psborrow

 
   hist(rgamma(1000,1, 1)) # gamma hyperprior from psborrow
   
 mean(rgamma(1000,1, 1)) # gamma hyperprior from psborrow

    hist(rgamma(1000,10, 10)) # gamma hyperprior from psborrow
 mean(rgamma(1000,10, 10)) # gamma hyperprior from psborrow



 mean(rgamma(1000,10, 10)) # gamma hyperprior from psborrow
 var(rgamma(1000,10, 10)) # gamma hyperprior from psborrow

hist(rnorm(10000,0, rgamma(1000,1.821, 1)) )



hist(rnorm(10000,-1.2, rgamma(10000,1, 3)) )
 


#If there is no information in the new trial, then we would want it to revert to something like the skeptical prior from the original trial N(0,sd=1.821) in terms of coverage of extremes P(VE>95%)

#if prior mean= -0.5, want to have a variance that is large enough to cover null with decent probability, maybe something like a SD=2

hist(rnorm(10000, -1.2 ,2))
mean(rnorm(10000, -1.2, 2.0)< log(0.05) )
mean(rnorm(10000, -1.2, 2.0)> log(1/0.05) )

#how does this combine with a prior mean?
hist(rnorm(10000,-1.2, rgamma(10000,1, 1)) )


hist(rnorm(10000,-1.2, rgamma(10000,5, 2)) )

hist(rnorm(10000, mean=prior.mean.logrr, sd=1.8)) #infrmative mean

```

VE from Novavax. This function takes the mean and 95% CI of vaccine effectiveness, converts to a log(RR), and estimates the precision fro the CIs, assuming it is symmetric around the mean (on scale of log(RR))

```{r}
#This is main novavax result
#prior.data <- prior.data.func(ve.obs.mean=39.4, ve.obs.lcl=5.3, ve.obs.ucl=61.2, N_cases_orig =c(35,41), pop_orig= c(1430, 2765))

#Table 2: per protocol 90 day VE against RSV LRTI with hospitalization
prior.data <- prior.data.func(ve.obs.mean=44.4, ve.obs.lcl=19.6, ve.obs.ucl=61.5, N_cases_orig =c(53,57), pop_orig= c(1430, 2765))
```

Set parameters for new trial population

```{r}


#proportion of kids in placebo group who get the outcome set to 0.037 for Noavavx secondary endpoint
set.rate.control <- 0.037 

#VE in the new trial, if it agrees with original
set.ve.new.trial.agree = prior.data$ve.obs.mean

#discrepant VE in new trial (higher)
set.ve.new.trial.higher = 90

```

## Scenario 1: Same VE as original trial

Novavax had 3045 vaccinees and 1581 controls; 2.4% in placebo group had primary endpoint in 90 days

```{r}
set.seed(123)

sim1 <- sim.data.func(N.vax = 3000,
                     N.control = 3000,
                     rate.control = set.rate.control,
                     ve.new.trial=set.ve.new.trial.agree,
                     n.sim=N.sim)
  
```

## Scenario 2: no actual effect

```{r}
set.seed(456)

sim2 <- sim.data.func(N.vax = 3000,
                     N.control = 3000,
                     rate.control = set.rate.control,
                     ve.new.trial=0,
                     n.sim=N.sim)
```

## Scenario 3: discrepant actual effect

```{r}
set.seed(456)

sim3 <- sim.data.func(N.vax = 3000,
                     N.control = 3000,
                     rate.control = set.rate.control,
                     ve.new.trial=set.ve.new.trial.higher,
                     n.sim=N.sim)
```

## At which points do we want to evaluate VE?

How many vaccinees are enrolled?

```{r}

eval_points = seq(from=200, to=3000, by=100)

```

test run. 500 vaccinees+500 controls. It takes 0.27 seconds for 1 run. Need to test 20 cut points, 3 true VE levels, and 500 sims (30,000 combos). This would take 135 minutes if run sequentially.

generate combined dataset for modeling--need to fill 0s if stratum not represented

```{r}
mod.ds1 <- lapply(eval_points , function(set.pop){
    ds1 <- bind_rows(sim1, sim2, sim3) %>%
        group_by(repN,ve.new.trial, treatment) %>%
      mutate(pop=set.pop) %>%
      filter( ((ID<pop &treatment=='placebo')|(ID<pop &treatment=='vax') )) %>%
          summarize(N_cases= n(), pop=mean(pop)) %>%
      ungroup() 
    
return(ds1)
    }) %>%
  bind_rows() %>%
  tidyr::complete(repN,ve.new.trial, treatment , pop,   fill=list(N_cases=0)) %>%
    mutate( vax=if_else(treatment=='vax',1,0)) 

```

Accrual of cases in the 500 trials

```{r}
mod.ds1 %>%
  filter(ve.new.trial==44.4) %>%
  ggplot(aes(x=pop, y=N_cases, group=interaction(repN, treatment), color=treatment))+
  geom_line(alpha=0.1) +
  theme_classic()
```

VE of 90

```{r}
mod.ds1 %>%
  filter(ve.new.trial==90) %>%
  ggplot(aes(x=pop, y=N_cases, group=interaction(repN, treatment), color=treatment))+
  geom_line(alpha=0.1) +
  theme_classic()
```

total accrual across both arms

```{r}
mod.ds1 %>%
  filter(ve.new.trial==44.4) %>%
  group_by(repN, pop) %>%
  summarize( tot_cases=sum(N_cases)) %>%
  ggplot(aes(x=pop, y=tot_cases, group=repN)) +
  geom_line(alpha=0.1) +
  theme_classic()+
  ylim(0,NA)
```

total accrual across both arms if VE=0

```{r}
mod.ds1 %>%
  filter(ve.new.trial==0) %>%
  group_by(repN, pop) %>%
  summarize( tot_cases=sum(N_cases)) %>%
  ggplot(aes(x=pop, y=tot_cases, group=repN)) +
  geom_line(alpha=0.1) +
  theme_classic()+
  ylim(0,NA) 
```

Prior selection based on: <https://hbiostat.org/proj/covid19/bayesplan.html>

```{r}
ptm <- proc.time()

mod1 <- mod.ds1 %>%
  filter( repN==3  & pop==1000 &  ve.new.trial==0 ) %>%
  call_jags( prior.mean=0, prior.prec=1/1.821^2) %>%
  filter(parm=='beta1')

proc.time() - ptm


mod2 <- mod.ds1 %>%
  filter( repN==3  & pop==1000 &  ve.new.trial==0 ) %>%
  call_jags( prior.mean=0, prior.prec=1/1.821^2, model.select=model_string_commensurate_gamma) %>%
  filter(parm=='beta1')

mod3 <- mod.ds1 %>%
  filter( repN==3  & pop==1000 &  ve.new.trial==0 ) %>%
  call_jags( prior.sd.upper=100, model.select=model_string_commensurate_uniform) %>%
  filter(parm=='beta1')


mod4 <- mod.ds1 %>%
  filter( repN==3  & pop==1000 &  ve.new.trial==0 ) %>%
  call_jags( prior.mean=prior.data$log_irr.obs[1], prior.prec=1/1.821^2, model.select=model_string_commensurate_half_cauchy) %>%
  filter(parm=='beta1')



bind_rows(mod1, mod2, mod3, mod4)
```

Check dataset--shows all strata have 2 obs

```{r}

mod.ds1 %>%
  group_by(repN, ve.new.trial, pop) %>%
  summarize(n_obs=n()) %>%
  ungroup() %>%
  summarize(n_obs_range =range(n_obs))

mod.ds1.c <- mod.ds1 %>%
  reshape2::dcast(repN+ve.new.trial + pop ~ treatment, value.var='N_cases') %>%
  mutate(N_events =placebo+vax, raw_RR= vax/placebo)

mod.ds1.c %>%
  filter(repN<500 & ve.new.trial==44.4 &N_events>20) %>%
ggplot(aes(x=N_events, y=raw_RR, group=repN))+
  geom_line(alpha=0.1)+
  theme_classic()

mod.ds1.c %>%
  filter(repN<500 & ve.new.trial==44.4&N_events>20) %>%
ggplot(aes(x=pop, y=raw_RR, group=repN))+
  geom_line(alpha=0.1)+
  theme_classic()
```

```{r}
 obs_data <- mod.ds1.c %>%
  mutate( obs_RR = vax/placebo )
```

Call all models for both simulated datasets. Use parallel processing; takes \~60 min to run on 7 cores

```{r, eval=F}
ptm <- proc.time()


ptm <- proc.time()
numCores = detectCores() -1
cluster1 <- new_cluster(numCores)
cluster_library(cluster1, "dplyr")
cluster_library(cluster1, "rjags")
cluster_library(cluster1, "HDInterval")
cluster_copy(cluster1, "call_jags")
cluster_copy(cluster1, "prior.data")
cluster_copy(cluster1, "model_string_basic_pois")
cluster_copy(cluster1, "model_string_commensurate_gamma")
cluster_copy(cluster1, "model_string_commensurate_uniform")
cluster_copy(cluster1, "model_string_commensurate_half_cauchy")


#baseline model without pooling, uninformative prior

mod0 <- mod.ds1 %>%
  group_by(repN,pop,ve.new.trial) %>%
  multidplyr::partition(cluster=cluster1)  %>%
 do(call_jags(., prior.mean=0,  prior.prec=1e-4 )) %>%
  collect() %>%
  ungroup() %>%
  mutate(Pooling_type='none', prior.info= 'N(0,1e-4)' ) 

#baseline model without pooling
mod1 <- mod.ds1 %>%
  group_by(repN,pop,ve.new.trial) %>%
  multidplyr::partition(cluster=cluster1)  %>%
 do(call_jags(. , model.select=model_string_basic_pois, prior.mean=0,  prior.prec=1/1.821^2 ) ) %>%
  collect() %>%
  ungroup() %>%
  mutate(Pooling_type='none', prior.info= 'N(0,1.821^2)')

proc.time() - ptm

ptm <- proc.time()
#Commensurate prior pooling
mod3 <- mod.ds1 %>%
  group_by(repN,pop,ve.new.trial) %>%
  multidplyr::partition(cluster=cluster1)  %>%
 do(call_jags(., prior.mean=prior.data$log_irr.obs[1], prior.prec=1/1.821^2, model.select=model_string_commensurate_gamma, set_tau_shp=3, set_tau_rate=2)) %>%
  collect() %>%
  ungroup() %>%
  mutate(Pooling_type='Commensurate', prior.info= 'N(u1,tau); IG(3,2)')

proc.time() - ptm

#Commensurate prior pooling
mod4 <- mod.ds1 %>%
  group_by(repN,pop,ve.new.trial) %>%
  multidplyr::partition(cluster=cluster1)  %>%
 do(call_jags(., prior.mean=prior.data$log_irr.obs[1], prior.prec=1/1.821^2, model.select=model_string_commensurate_gamma, set_tau_shp=0.01, set_tau_rate=0.01)) %>%
  collect() %>%
  ungroup() %>%
  mutate(Pooling_type='Commensurate', prior.info= 'N(u1,tau); IG(0.01,0.01)')

mod5 <- mod.ds1 %>%
  group_by(repN,pop,ve.new.trial) %>%
  multidplyr::partition(cluster=cluster1)  %>%
 do(call_jags(., prior.mean=prior.data$log_irr.obs[1], prior.prec=1/1.821^2, model.select=model_string_commensurate_uniform, prior.sd.upper=2)) %>%
  collect() %>%
  ungroup() %>%
  mutate(Pooling_type='Commensurate', prior.info= 'N(\u03B4\u2081,\u03C3\u00B2); \u03C3~U(0,2)')

mod6 <- mod.ds1 %>%
  group_by(repN,pop,ve.new.trial) %>%
  multidplyr::partition(cluster=cluster1)  %>%
 do(call_jags(., prior.mean=prior.data$log_irr.obs[1], prior.prec=1/1.821^2, model.select=model_string_commensurate_uniform, prior.sd.upper=100)) %>%
  collect() %>%
  ungroup() %>%
  mutate(Pooling_type='Commensurate', prior.info= 'N(u1,tau); U(0,100)')

mod7 <- mod.ds1 %>%
  group_by(repN,pop,ve.new.trial) %>%
  multidplyr::partition(cluster=cluster1)  %>%
 do(call_jags(. , prior.mean=prior.data$log_irr.obs[1], prior.prec=1/1.821^2, model.select=model_string_commensurate_half_cauchy)) %>%
  collect() %>%
  ungroup() %>%
  mutate(Pooling_type='Commensurate', prior.info= 'N(u1,tau);t(0,1,1)T(0,)')


proc.time() - ptm

all.mods <- bind_rows(mod0,mod1, mod3, mod4, mod5, mod6, mod7)

saveRDS(all.mods, './Results/all.mods_novavax.rds')

ptm <- proc.time()

```

```{r}
all.mods <- readRDS( './Results/all.mods_novavax.rds') %>%
  rename(parm_mean=mean) %>%
  mutate(  parm_mean = if_else( parm_mean< (-4), (-4), 
                  if_else(parm_mean> 4,4, parm_mean)),
          lcl = if_else(lcl< (-4), (-4), 
                  if_else(lcl>4,4, lcl)),
           ucl = if_else(ucl< (-4), (-4), 
                  if_else(ucl>4,4, ucl)),
          total_pop=pop*2
         
  ) 


```

## How much agreement is there between prior and observed data

```{r}
N=100000
prec.sim <- rgamma(n=N, shape=0.01, rate=0.01)
var.sim <- 1/prec.sim
alpha.sim <- exp(-1*var.sim)
hist(alpha.sim)

full.prior <- rnorm(N, mean=rep(0,N) , sd=sqrt(var.sim) )
hist(full.prior)
```

In these models tau is precision = 1/var
If precision is high, this indicates a lot of pooling on commensurate prior

```{r}

check.tau <- readRDS( './Results/all.mods_novavax.rds') %>%
  filter(parm=='beta1' & Pooling_type !='none') 

ggplot(data=check.tau, aes(x=pop, y=alpha1_mean) )+
         geom_point(alpha=0.01) +
  theme_classic() +
  facet_grid(ve.new.trial~prior.info)


# check.tau %>% 
#   filter(prior.info== "N(u1,tau); IG(0.01,0.01)") %>%
#   ggplot( aes(x=pop, y=alpha1_mean) )+
#          geom_point(alpha=0.01) +
#   theme_classic() +
#   facet_grid(~ve.new.trial)

check.tau %>% 
  filter(prior.info== "N(u1,tau); IG(0.01,0.01)") %>%
  ggplot( aes(x=pop, y=alpha1_median) )+
    geom_hline(yintercept=c(0.2, 0.4, 0.6, 0.8,1.0), lty=2, col='black', alpha=0.1)+
         geom_point(alpha=0.01) +
  theme_classic() +
  facet_grid(~ve.new.trial)+
  ylim(0,1) +
  ylab('Agreement between prior and data')
```

Figure S2 
```{r}
check.tau %>% 
  filter(prior.info== "N(u1,tau); IG(0.01,0.01)") %>%
  mutate(     new.ve.label = if_else(ve.new.trial==0, 'A) VE=0%',
                            if_else(ve.new.trial==44.4, 'B) VE=44.4%',
                            if_else(ve.new.trial==90, 'C) VE=90%', NA_character_)))) %>%
  ggplot( aes(x=100*(1-exp(mean)), y=alpha1_median) )+
    geom_hline(yintercept=c(0.2, 0.4, 0.6, 0.8,1.0), lty=2, col='black', alpha=0.1)+
      geom_vline(xintercept=c(44.4), lty=2, col='gray')+

         geom_point(alpha=0.01) +
  theme_classic() +
  facet_grid(~new.ve.label)+
  ylim(0,1) +
  xlim(c(-100,100))+
  ylab('Agreement between prior and data') +
  xlab ('estimate of VE (%)')

check.tau %>% 
  filter(prior.info== "N(u1,tau); IG(0.01,0.01)") %>%
  ggplot( aes(x=100*(1-exp(mean)), y=alpha1_median, color=as.factor(ve.new.trial)) )+
    geom_hline(yintercept=c(0.2, 0.4, 0.6, 0.8,1.0), lty=2, col='black', alpha=0.1)+
         geom_point(alpha=0.01) +
  theme_classic() +
  #facet_grid(~ve.new.trial)+
  ylim(0,1) +
  xlim(c(-100,100))+
  ylab('Agreement between prior and data') +
  xlab ('estimate of VE')



```

## Compare delta and beta estimates
```{r}
check.beta <- readRDS( './Results/all.mods_novavax.rds') %>%
  filter(parm=='beta1' & Pooling_type !='none')   %>%
  select(prior.info, mean, repN, pop, ve.new.trial) %>%
  rename(beta_est=mean)

check.delta <- readRDS( './Results/all.mods_novavax.rds') %>%
  filter(parm=='delta' & Pooling_type !='none') %>%
  select(prior.info, mean, repN, pop, ve.new.trial) %>%
  rename(delta_est=mean)

comp1 <- check.delta %>%
  left_join(check.beta, by=c('prior.info','repN','pop', 've.new.trial'))

comp1 %>% 
  filter(prior.info== "N(u1,tau); IG(0.01,0.01)") %>%
  ggplot( aes(x=100*(1-exp(delta_est)), y=100*(1-exp(beta_est)), color=as.factor(ve.new.trial)) )+
    geom_hline(yintercept=c(0.2, 0.4, 0.6, 0.8,1.0), lty=2, col='black', alpha=0.1)+
      geom_vline(xintercept=c(30,40,50), lty=2, col='black', alpha=0.1)+

         geom_point(alpha=0.01) +
  theme_classic() +
  #facet_grid(~ve.new.trial)+
  ylim(0,100) +
  xlim(c(0,100))+
  ylab('estimate of VE original') +
  xlab ('estimate of VE original')
  
```



## Bayesian estimates of performance

```{r}

#the p_xx are replicated for each parameter so filter to just 1 parameter
probs <- all.mods %>% 
  group_by(repN, pop, ve.new.trial,Pooling_type,prior.info) %>%
  filter(parm=='beta1')  %>%
  ungroup() %>%
  mutate(prior.info = if_else(prior.info=="N(0,1.821^2)", "N(0,3.32)", prior.info ),
         prior.info = if_else(prior.info=="N(u1,tau); IG(3,2)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(3,2)", prior.info ),
         prior.info = if_else(prior.info=="N(u1,tau); IG(0.01,0.01)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(0.01,0.01)", prior.info ),
         prior.info = if_else(prior.info=="N(u1,tau); tau~U(0,2)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3~U(0,2)", prior.info ),
         prior.info = if_else(prior.info=="N(u1,tau); U(0,100)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3~U(0,100)", prior.info ),
         prior.info = if_else(prior.info=="N(u1,tau);t(0,1,1)T(0,)", "N(\u03B4\u2081,\u03C3\u00B2);\u03C3~t(0,1,1)T(0,)", prior.info ),
         
         prior.label = if_else(prior.info=="N(0,3.32)", 'No pooling', 
                  if_else(prior.info== "N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(3,2)", 'With pooling (\u03C3\u00B2~IG(3,2))',
                  if_else(prior.info== "N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(0.01,0.01)", 'With pooling (\u03C3\u00B2~IG(0.01,0.01))',
                  if_else(prior.info== "N(\u03B4\u2081,\u03C3\u00B2); \u03C3~U(0,2)", 'With pooling (\u03C3~U(0,2))',
                 if_else(prior.info== "N(\u03B4\u2081,\u03C3\u00B2); \u03C3~U(0,100)", 'With pooling  (\u03C3~U(0,100))',
      if_else(prior.info== "N(u1,tau);t(0,1,1)T(0,)", "N(\u03B4\u2081,\u03C3\u00B2);\u03C3~t(0,1,1)T(0,)", 'With pooling  (\u03C3~t(0,1,1)T(0,))'
      ))))))                                                                         
      )

         
```

## Figure S1: Number of participants vs Number of events

```{r}
obs_data %>%
  filter(ve.new.trial==44.4)  %>%
  ggplot( aes(x=N_events, y=pop*2, group=repN)) +
  geom_point(alpha=0.1) +
  ylab('Total Participants')+
  xlab('Number of events') +
  theme_classic() +
  geom_hline(yintercept=4195, lty=2, col='red') +
  geom_vline(xintercept=110, lty=2, col='red')
```

## How do posteriors from a single trial evolve as data accrue?

### Figure 1


```{r, fig.width=4, fig.height=7}
fig2a <- probs %>%
    left_join(obs_data, by=c('repN','pop','ve.new.trial')) %>%
  filter( ve.new.trial==44.4 & Pooling_type != 'Full' & repN==1 & prior.info %in% c("N(0,3.32)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(0.01,0.01)"  )) %>%
  mutate(Pooling_type = as.factor(Pooling_type),
         VE = 100*(1-exp(parm_mean)),
         VE_ucl=100*(1-exp(ucl)),
         VE_lcl=100*(1-exp(lcl)),
         obs_VE = 100*(1-obs_RR),
         #VE_ucl=if_else(rr_ucl>=4,3.99,rr_ucl),
        # VE_lcl=if_else(rr_lcl<0.01,0.01,rr_lcl),
         Prior.label = if_else(prior.info == "N(0,3.32)", 'No pooling', 
                               'With pooling')) %>%
ggplot( aes(x=pop*2, y=(VE), group=interaction(repN,prior.info))) +
  geom_line( alpha=1) +
  theme_classic() +
   geom_ribbon( aes(x=pop*2, ymin=VE_lcl, 
                   ymax=VE_ucl), alpha=0.3) +
  facet_wrap(~Prior.label) +
  ylab('Vaccine effectiveness (%)') +
  xlab('Total Participants')+
  geom_line(aes(x=pop*2, y=(obs_VE)), color='black', lty=3)+
  coord_cartesian( ylim=c(-10,100)) +
  geom_hline(yintercept=44.4, lty=2)+
    theme(text=element_text( family="serif"))
  

fig2b <- probs %>%
    left_join(obs_data, by=c('repN','pop','ve.new.trial')) %>%
  filter( ve.new.trial==90 & Pooling_type != 'Full' & repN==1 & prior.info %in% c("N(0,3.32)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(0.01,0.01)"  )) %>%
  mutate(Pooling_type = as.factor(Pooling_type),
         VE = 100*(1-exp(parm_mean)),
         VE_ucl=100*(1-exp(ucl)),
         VE_lcl=100*(1-exp(lcl)),
         obs_VE = 100*(1-obs_RR),
         #VE_ucl=if_else(rr_ucl>=4,3.99,rr_ucl),
        # VE_lcl=if_else(rr_lcl<0.01,0.01,rr_lcl),
         Prior.label = if_else(prior.info == "N(0,3.32)", 'No pooling', 
                               'With pooling')) %>%
ggplot( aes(x=pop*2, y=(VE), group=interaction(repN,prior.info))) +
  geom_line( alpha=1) +
  theme_classic() +
   geom_ribbon( aes(x=pop*2, ymin=VE_lcl, 
                   ymax=VE_ucl), alpha=0.3) +
  facet_wrap(~Prior.label) +
  ylab('Vaccine effectiveness (%)') +
  xlab('Total Participants')+
  geom_line(aes(x=pop*2, y=(obs_VE)), color='black', lty=3)+
  coord_cartesian( ylim=c(-10,100)) +
  geom_hline(yintercept=90, lty=2)+
    theme(text=element_text( family="serif"))

fig2c <- probs %>%
    left_join(obs_data, by=c('repN','pop','ve.new.trial')) %>%
  filter( ve.new.trial==0 & Pooling_type != 'Full' & repN==1 & prior.info %in% c("N(0,3.32)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(0.01,0.01)"  )) %>%
  mutate(Pooling_type = as.factor(Pooling_type),
         VE = 100*(1-exp(parm_mean)),
         VE_ucl=100*(1-exp(ucl)),
         VE_lcl=100*(1-exp(lcl)),
         obs_VE = 100*(1-obs_RR),
         #VE_ucl=if_else(rr_ucl>=4,3.99,rr_ucl),
        # VE_lcl=if_else(rr_lcl<0.01,0.01,rr_lcl),
         Prior.label = if_else(prior.info == "N(0,3.32)", 'No pooling', 
                               'With pooling')) %>%
ggplot( aes(x=pop*2, y=(VE), group=interaction(repN,prior.info))) +
  geom_line( alpha=1) +
  theme_classic() +
   geom_ribbon( aes(x=pop*2, ymin=VE_lcl, 
                   ymax=VE_ucl), alpha=0.3) +
  facet_wrap(~Prior.label) +
  ylab('Vaccine effectiveness (%)') +
  xlab('Total Participants')+
  geom_line(aes(x=pop*2, y=(obs_VE)), color='black', lty=3)+
  coord_cartesian( ylim=c(-50,50)) +
  geom_hline(yintercept=0, lty=2)+
    theme(text=element_text( family="serif"))

fig1 <- (fig2a+ ggtitle('True effectiveness in new trial=44.4%') +     theme( strip.background = element_blank() ))/
  (fig2b+ ggtitle('True effectiveness in new trial=90%')  +  theme( strip.background = element_blank() ))/
  (fig2c+ ggtitle('True effectiveness in new trial=0%')+     theme( strip.background = element_blank() ) ) + 
  plot_annotation(tag_levels = 'A')

fig1

tiff("./Results/Figs/fig1.tiff",  res=600,width=4, height=7, units='in')
fig1
dev.off()
```

### Trajectories of P(benefit) for a subset of trials

This shows how the probabilities evolve as more data accrue. This shows probability of any benefit: P(VE\>0\|data)

```{r}
probs %>%
  filter(ve.new.trial==44.4 & Pooling_type != 'Full' & repN %in% c(1:10)) %>%
  mutate(prior.info = as.factor(prior.info)) %>%
  
ggplot( aes(x=total_pop, y=p_0, group=interaction(repN,prior.info), color=Pooling_type)) +
  geom_line( alpha=0.1) +
  theme_classic() +
  facet_wrap(~prior.info) +
    theme(text=element_text( family="serif"))

```

## Power calculations

What proportion of trials would have a P(benefit\>0.95\|data) at different sample sizes?

for each trial, call 0/1 based on if we have hit 0.95 at this sample size or earlier

```{r, fig.width=5, fig.height=3}
power1 =probs %>%
  arrange(repN, ve.new.trial,prior.info ,total_pop) %>%
 group_by(repN,  ve.new.trial,prior.info) %>%
  mutate( eff_0 = if_else(p_0>0.95,1,0),
          eff_0_p99 = if_else(p_0>0.99,1,0),
          eff_0_p97 = if_else(p_0>0.98,1,0),
          eff_0_p98 = if_else(p_0>0.97,1,0),

                    eff_0_p96 = if_else(p_0>0.96,1,0),
          eff_5 = if_else(p_0_95>0.95,1,0),
          eff_10 = if_else(p_0_90>0.95,1,0),
          eff_15 = if_else(p_0_85>0.95,1,0),
          eff_20 = if_else(p_0_80>0.95,1,0),
          eff_25 = if_else(p_0_75>0.95,1,0),
          eff_30 = if_else(p_0_70>0.95,1,0),
          futile_30 = if_else(p_futile>0.95,1,0),
            ) %>%
  ungroup() %>%
  group_by(total_pop, ve.new.trial,prior.info) %>%
    summarize(prop_eff_0 = mean(eff_0),
              prop_eff_0_p99=mean(eff_0_p99),
              prop_eff_0_p98=mean(eff_0_p98),
              prop_eff_0_p97=mean(eff_0_p97),
              prop_eff_0_p96=mean(eff_0_p96),

                            prop_eff_5 = mean(eff_5),
            prop_eff_10 = mean(eff_10),
            prop_eff_15 = mean(eff_15),
            prop_eff_20 = mean(eff_20),
              prop_eff_25 = mean(eff_25),
              prop_eff_30 = mean(eff_30),
              prop_futile_30 = mean(futile_30),
              )        
```

### Power for sequential trial with control of type 1 error

We want to control Type I error for overall trial. i.e., by the final expected data look. this will lead to lower type I error and power at earlier looks. This uses a simple rule with same cutoff at all time points. There might be more efficient spending functions that have very low prob at early points and higher at later points

First find the threshold, then apply it

This first function iterates through different alpha to look at Type I error

```{r}
threshold_finder_sequential <- function(test.alpha=0.05, eval.points ){

find_alpha =probs %>%
    filter(total_pop %in% eval.points & ve.new.trial==0) %>%
  arrange(repN, prior.info ,total_pop) %>%
 group_by(repN,  prior.info, prior.label) %>%
  mutate( eff0 = 1*(p_0>(1-test.alpha)),
          eff30 = 1*(p_0_70>(1-test.alpha)) ,
            stop0 = cumsum(eff0) ,
          stop0 = if_else(stop0>=1,1,0),
          stop30 = cumsum(eff30) ,
          stop30 = if_else(stop30>=1,1,0)) %>%
  ungroup() %>%
  group_by(total_pop, prior.info,prior.label) %>%
    filter(total_pop==eval.points[length(eval.points)] ) %>%
       summarize(prop_stopped0 = mean(stop0),prop_stopped30 = mean(stop30)) %>%
   mutate(test.threshold= (1-test.alpha))
}

```

Function applies threshold identified by threshold_finder_sequential

```{r}
threshold_application <- function( set.eval.points){
    seq.tresholds <- pbapply::pblapply(seq(0.001,0.05, by=0.001), threshold_finder_sequential, eval.points=set.eval.points) %>%   bind_rows()
   
    #find (1-alpha) that returns acceptable type 1 error  
   min.tresholds1 <- seq.tresholds %>%
    arrange(prior.info, test.threshold) %>%
    group_by(prior.info,prior.label) %>%
    mutate(  acceptable30 = (prop_stopped30 <=0.05), rank_acceptable30=cumsum(acceptable30) ,
             acceptable0 = (prop_stopped0 <=0.05), rank_acceptable0=cumsum(acceptable0) ) 
   
   min.tresholdS3a <- min.tresholds1 %>%
    filter(rank_acceptable0==1) %>%
     select(prior.info,test.threshold,prior.label) %>%
     rename(test.threshold0=test.threshold)
   
      min.tresholdS3b <- min.tresholds1  %>%
    filter(rank_acceptable30==1) %>%
     select(prior.info,prior.label,test.threshold) %>%
     rename(test.threshold30=test.threshold)
   
      min.tresholds <- min.tresholdS3a %>%
        left_join(min.tresholdS3b, by='prior.info')
   
   stopped_efficacy1 =probs %>%
    filter(total_pop %in% set.eval.points ) %>%
     left_join(min.tresholds, by=c('prior.info')) %>%
    arrange(repN, ve.new.trial,prior.info ,total_pop) %>%
   group_by(repN,  ve.new.trial,prior.info,prior.label) %>%
    mutate( eff30 = 1*(p_0_70>test.threshold30),
              stop30 = cumsum(eff30) ,
            stop30 = if_else(stop30>=1,1,0),
            eff0 = 1*(p_0>test.threshold0),
              stop0 = cumsum(eff0) ,
            stop0 = if_else(stop0>=1,1,0),
            
            new.ve.label = if_else(ve.new.trial==0, 'A) VE=0%',
                            if_else(ve.new.trial==44.4, 'B) VE=44.4%',
                            if_else(ve.new.trial==90, 'C) VE=90%', NA_character_)))
            ) %>%
    ungroup() %>%
    group_by(total_pop, ve.new.trial,new.ve.label,prior.info,prior.label) %>%
      summarize(prop_stopped0 = mean(stop0),
                prop_stopped30 = mean(stop30))
   
   stopped_efficacy_Pop1 =probs %>%
    filter(total_pop %in% set.eval.points ) %>%
     left_join(min.tresholds, by=c('prior.info')) %>%
    arrange(repN, ve.new.trial,prior.info ,total_pop) %>%
   group_by(repN,  ve.new.trial,prior.info,prior.label) %>%
    mutate( eff30 = 1*(p_0_70>test.threshold30),
              stopN30 = cumsum(eff30), 
                stopN30=if_else(stopN30==0 & total_pop==6000,1,stopN30), 
            eff0 = 1*(p_0>test.threshold0),
              stopN0 = cumsum(eff0), 
                stopN0=if_else(stopN0==0 & total_pop==6000,1,stopN0),
            new.ve.label = if_else(ve.new.trial==0, 'A) VE=0%',
                            if_else(ve.new.trial==44.4, 'B) VE=44.4%',
                            if_else(ve.new.trial==90, 'C) VE=90%', NA_character_)))
            )  %>%
    ungroup() %>%
         group_by(total_pop, ve.new.trial,new.ve.label,prior.info) 

   
  stopped_efficacy_Pop2a <- stopped_efficacy_Pop1 %>%
       filter(stopN0==1)%>%
       ungroup() 
  
    stopped_efficacy_Pop2b <- stopped_efficacy_Pop1 %>%
       filter(stopN30==1)%>%
       ungroup() 
   
   ave_pop_stop0 <- stopped_efficacy_Pop2a %>%
       ungroup() %>%
       group_by(ve.new.trial,prior.info,prior.label) %>%
       summarize(ave_pop_stop0=mean(total_pop))
   
  ave_pop_stop30 <- stopped_efficacy_Pop2b %>%
       ungroup() %>%
       group_by(ve.new.trial,prior.info,prior.label) %>%
       summarize(ave_pop_stop30=mean(total_pop))
   
  ave_pop_stop <- ave_pop_stop0 %>%
    full_join(ave_pop_stop30, by=c('ve.new.trial','prior.info'))
     
   p1a <- ggplot(seq.tresholds, aes( x=test.threshold, y=prop_stopped0, group=prior.label, color=prior.label, shape=prior.label)) + 
  geom_line() +
  theme_classic() +
  geom_hline(yintercept=0.05)+
       guides(color = guide_legend(title = "Prior type") ,shape = guide_legend(title = "Prior type"))+
    theme(text=element_text( family="serif"))
p1a


p2a <- stopped_efficacy1 %>%
  filter( prior.info %in% c("N(0,3.32)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(0.01,0.01)"  )) %>%
ggplot(aes( x=total_pop, y=prop_stopped0, group=prior.label, color=prior.label, shape=prior.label, linetype=prior.label )) +
  geom_line() +
    geom_point() +
  facet_wrap(~ new.ve.label) +
  theme_classic() +
  guides(color = guide_legend(title = "Prior type") ,shape = guide_legend(title = "Prior type"),linetype = guide_legend(title = "Prior type"))+
  geom_hline(yintercept = c(0.05, 0.8), lty=2, col='gray')+
  ggtitle('Cumulative proportion with >X% prob that effect is >0%VE')+
    theme(text=element_text( family="serif")) +
     theme( strip.background = element_blank() )

p2a1 <- stopped_efficacy1 %>%
ggplot(aes( x=total_pop, y=prop_stopped0, group=prior.label, color=prior.label, shape=prior.label)) +
  geom_line() +
    geom_point() +
  facet_wrap(~ new.ve.label) +
     theme( strip.background = element_blank() )+
  theme_classic() +
    guides(color = guide_legend(title = "Prior type") ,shape = guide_legend(title = "Prior type"))+
  geom_hline(yintercept = c(0.05, 0.8), lty=2, col='gray')+
  ggtitle('Cumulative proportion with >X% prob that effect is >0%VE')+
    theme(text=element_text( family="serif"))


p2b <- stopped_efficacy1 %>%
ggplot(aes( x=total_pop, y=prop_stopped30, group=prior.label, color=prior.label, shape=prior.label)) +
  geom_line() +
    geom_point() +
  facet_wrap(~ new.ve.label) +
  theme_classic() +
     theme( strip.background = element_blank() )+
    guides(color = guide_legend(title = "Prior type") ,shape = guide_legend(title = "Prior type"))+
  geom_hline(yintercept = c(0.05, 0.8), lty=2, col='gray')+
  ggtitle('Cumulative proportion with >X% prob that effect is >30%VE')+
    theme(text=element_text( family="serif"))


   out.obj = list('min.tresholds'=min.tresholds, 'stopped_efficacy1'=stopped_efficacy1,'ave_pop_stop'=ave_pop_stop,'seq.tresholds'=seq.tresholds,'p1a'=p1a,'p2a'=p2a,'p2a1'=p2a1,'p2b'=p2b)
   return(out.obj)
}
```

call the functions and look at output

Every 1000, starting at 1000 Participants

```{r}

seq.results1 <-threshold_application( set.eval.points=c(1000, 2000,3000, 4000, 5000, 6000))

#seq.results1$p1


tab1 <- seq.results1$ave_pop_stop %>%
  full_join(seq.results1$min.tresholds, by=c('prior.info')) %>%
  mutate(ave_pop_stop30=round(ave_pop_stop30,-1),ave_pop_stop0=round(ave_pop_stop0,-1)) %>%
  filter(!is.na(ve.new.trial))

ggplotly(seq.results1$p2a)

ggplotly(seq.results1$p2b)


write.csv(tab1, './Results/Table_sequential_probs.csv')


```

### Figure 2 for manuscript

```{r, fig.width=6, fig.height=3.5}
fig2 <- seq.results1$p2a+
  ggtitle('') +
  ylab('Proportion of trials stopped for efficacy') + 
  xlab('Total Number of Participants')+
  theme(axis.text.x=element_text(angle=45, hjust=1))


fig2

tiff("./Results/Figs/fig2.tiff",  res=600,width=6, height=3.5, units='in')
fig2
dev.off()
```

### Figure S3

```{r, fig.width=6, fig.height=3.5}
seq.results1$p2a1+
  ggtitle('') +
  ylab('Proportion of trials stopped for efficacy') + 
  xlab('Total Number of Participants')+
  theme(axis.text.x=element_text(angle=45, hjust=1))+
     theme( strip.background = element_blank() )

```

### Figure S4

```{r, fig.width=6, fig.height=3.5}
seq.results1$p2b+
  ggtitle('') +
  ylab('Proportion of trials stopped for efficacy') + 
  xlab('Total Number of Participant')+
  theme(axis.text.x=element_text(angle=45, hjust=1))+
     theme( strip.background = element_blank() )

```

### Power for fixed sample size, no type I error control

```{r, fig.width=5, fig.height=3}

p1 <- power1 %>%
  filter(prior.info %in% c("N(0,3.32)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(0.01,0.01)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3~U(0,2)" )) %>%
ggplot(aes( x=total_pop, y=prop_eff_0, group=prior.info, color=prior.info, lty=prior.info)) +
  geom_line() +
  facet_wrap(~ ve.new.trial) +
  theme_classic() +
  geom_hline(yintercept = c(0.05, 0.8), lty=2, col='gray')+
  ylab('Proportion') +
  ggtitle('Proportion of trials with >95% probability that effect is >0% VE')+
    theme(text=element_text( family="serif"))
p1
ggplotly(p1)
```

```{r, fig.width=5, fig.height=3}

power1 %>%
ggplot(aes( x=total_pop, y=prop_eff_0, group=prior.info, color=prior.info)) +
  geom_line() +
  facet_wrap(~ ve.new.trial) +
  theme_classic() +
  geom_hline(yintercept = c(0.05, 0.8), lty=2, col='gray')+
  ggtitle('Proportion of trials with >95% probability that effect is >0% VE')+
    theme(text=element_text( family="serif"))
```

### Posterior probability that the trial is Futile (VE \<30%)

```{r, fig.width=5, fig.height=3}
power1 %>%
        filter(prior.info %in% c("N(0,3.32)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(0.01,0.01)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3~U(0,2)"  )) %>%

ggplot(aes( x=total_pop, y=prop_futile_30, group=prior.info, color=prior.info, lty=prior.info)) +
  geom_line() +
  facet_wrap(~ ve.new.trial) +
  theme_classic() +
  geom_hline(yintercept = c(0.05, 0.8), lty=2, col='gray')+
  ggtitle('Proportion of trials with >95% probability that effect is <30% VE')+
    theme(text=element_text( family="serif"))
```

```{r, fig.width=5, fig.height=3}
power1 %>%
ggplot(aes( x=total_pop, y=prop_futile_30, group=prior.info, color=prior.info)) +
  geom_line() +
  facet_wrap(~ ve.new.trial) +
  theme_classic() +
  geom_hline(yintercept = c(0.05, 0.8), lty=2, col='gray')+
  ggtitle('Proportion of trials with >95% probability that effect is <30% VE')+
    theme(text=element_text( family="serif"))
```

### Power for sequential trial with looks every 1000 people, no control of type I error Figure S5

```{r, fig.width=6, fig.height=3.5}
stopped_efficacy1 =probs %>%
  filter(total_pop %in% c(1000,2000, 3000,4000, 5000, 6000) & prior.info %in% c("N(0,3.32)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(0.01,0.01)" )) %>%
  arrange(repN, ve.new.trial,prior.info ,total_pop) %>%
 group_by(repN,  ve.new.trial,prior.info,prior.label) %>%
  mutate( eff = 1*(p_0>0.95),
            stop = cumsum(eff) ,
          stop = if_else(stop>=1,1,0)) %>%
  ungroup() %>%
  group_by(total_pop, ve.new.trial,prior.info,prior.label) %>%
    summarize(prop_stopped = mean(stop))    %>%
  mutate( new.ve.label = if_else(ve.new.trial==0, 'A) VE=0%',
                            if_else(ve.new.trial==44.4, 'B) VE=44.4%',
                            if_else(ve.new.trial==90, 'C) VE=90%', NA_character_)))
          )


p1a <- stopped_efficacy1 %>%
ggplot(aes( x=total_pop, y=prop_stopped, group=prior.label, color=prior.label)) +
  geom_line() +
    geom_point() +
  facet_wrap(~ new.ve.label) +
  theme_classic() +
  geom_hline(yintercept = c(0.05, 0.8), lty=2, col='gray')+
  ggtitle('Cumulative proportion with >95% prob that effect is >0%VE')+
    theme(text=element_text( family="serif"))+
     theme( strip.background = element_blank() )+
  ylab("Proportion of trials stopped for efficacy")+
  xlab('Total number of participants')
ggplotly(p1a)
p1a
```

Futility (P VE\<30%)

```{r,fig.width=6, fig.height=3}
stopped_futility =probs %>%
    filter(total_pop %in% c(1000,2000, 3000,4000, 5000, 6000)) %>%
  arrange(repN, ve.new.trial,prior.info ,pop) %>%
 group_by(repN, total_pop, ve.new.trial,prior.info) %>%
  mutate( fut = 1*((1-p_0_70)>0.8),
            stop = cumsum(fut) ,
          stop = if_else(stop>=1,1,0)) %>%
  ungroup() %>%
  group_by(total_pop, ve.new.trial,prior.info) %>%
    summarize(prop_stopped = mean(stop))        


p2 <- ggplot(stopped_futility, aes( x=total_pop, y=prop_stopped, group=prior.info, color=prior.info)) +
  geom_line() +
  geom_point() +
  facet_wrap(~ ve.new.trial) +
  theme_classic() +
  geom_hline(yintercept = c(0.05, 0.8), lty=2, col='gray')+
    theme(text=element_text( family="serif"))

p2 
ggplotly(p2)
```

## Adjust threshold for Type I error control, static trial

First consider a static trial (one look at the data at specified point)

```{r, fig.width=6, fig.height=3}

#First figure out what cutoff should be when true VE=0 (ve.new.trial) to get 5% Type 1 error rate
#looks at posterior probs for VE>0|data, VE>30|data when true VE=0, look at what threshold captures 95% of the simulated trials
cutoff =probs %>%
    ungroup() %>%
#    filter(total_pop %in% seq(200, 6000, by=100)) %>%
  arrange( ve.new.trial,prior.info ,total_pop,repN) %>%
  group_by(total_pop, ve.new.trial,prior.info) %>%
  filter(ve.new.trial==0) %>%
  summarize(p_cutoff0= quantile(p_0,probs=0.95), p_cutoff30= quantile(p_0_70,probs=0.95)) %>%  #what quantile of P(VE>0%|Data) captures 95% of the trials? 
  ungroup() %>%
  select(p_cutoff0,p_cutoff30,total_pop, prior.info)
  
cutoff

    #Then merge in and apply this cutoff for different VE
stopped_efficacy1 =probs %>%
   right_join(cutoff, by=c( 'prior.info', 'total_pop')) %>%
  arrange(repN, ve.new.trial,prior.info ,total_pop) %>%
 group_by(repN,  ve.new.trial,prior.info,total_pop) %>%

  mutate( eff0 = 1*(p_0>p_cutoff0),
          eff30 = 1*(p_0_70>p_cutoff30))%>%
  ungroup() %>%
  group_by(total_pop, ve.new.trial,prior.info) %>%
    summarize(prop_stopped_eff0 = mean(eff0),
              prop_stopped_eff30 = mean(eff30))


p1a <- stopped_efficacy1 %>%
  filter(  prior.info %in% c("N(0,3.32)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(0.01,0.01)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3~U(0,2)"  )) %>%
ggplot(aes( x=total_pop, y=prop_stopped_eff0, group=prior.info, color=prior.info)) +
  geom_line() +
 #   geom_point() +
  facet_wrap(~ ve.new.trial) +
  theme_classic() +
  geom_hline(yintercept = c(0.05, 0.8), lty=2, col='gray')+
  ggtitle('Power and type I error with control for type I error, single evaluation')+
    theme(text=element_text( family="serif"))
ggplotly(p1a)
p1a


p1b <- stopped_efficacy1 %>%
  #filter(  prior.info %in% c("N(0,3.32)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(0.01,0.01)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3~U(0,2)"  )) %>%
ggplot(aes( x=total_pop, y=prop_stopped_eff0, group=prior.info, color=prior.info)) +
  geom_line() +
 #   geom_point() +
  facet_wrap(~ ve.new.trial) +
  theme_classic() +
  geom_hline(yintercept = c(0.05, 0.8), lty=2, col='gray')+
  ggtitle('Power and type I error with control for type I error, single evaluation')+
    theme(text=element_text( family="serif"))
ggplotly(p1b)
p1b

p1c <- cutoff %>%
ggplot(aes( x=total_pop, y=p_cutoff0, group=prior.info, color=prior.info)) +
  geom_line() +
   # geom_point() +
  theme_classic() +
  geom_hline(yintercept = c(0.95), lty=2, col='gray')+
  ggtitle('P cutoff with control for type I error; single evaluation')+
    theme(text=element_text( family="serif"))
ggplotly(p1c)
p1c

p1d <- stopped_efficacy1 %>%
  #filter(  prior.info %in% c("N(0,3.32)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(0.01,0.01)", "N(\u03B4\u2081,\u03C3\u00B2); \u03C3~U(0,2)"  )) %>%
ggplot(aes( x=total_pop, y=prop_stopped_eff30, group=prior.info, color=prior.info)) +
  geom_line() +
 #   geom_point() +
  facet_wrap(~ ve.new.trial) +
  theme_classic() +
  geom_hline(yintercept = c(0.05, 0.8), lty=2, col='gray')+
  ggtitle('Power and type I error with control for type I error, VE>30% single evaluation')+
    theme(text=element_text( family="serif"))
ggplotly(p1d)
p1d



```

### Observed data vs estimate

```{r}
obs_data %>%
  filter(ve.new.trial==44.4)  %>%
  ggplot( aes(x=N_events, y=obs_RR, group=repN)) +
  geom_line(alpha=0.1) +
  theme_classic() +
  ylim(0.2, 2)
```

```{r}
obs_data %>%
  filter(ve.new.trial==44.4)  %>%
  ggplot( aes(x=N_events, y=pop*2, group=repN)) +
  geom_point(alpha=0.1) +
  ylab('Total Participants')+
  xlab('Number of events') +
  theme_classic() +
  geom_hline(yintercept=4195, lty=2, col='red') +
  geom_vline(xintercept=110, lty=2, col='red')
```

```{r}
obs_data %>%
  filter(ve.new.trial==44.4)  %>%
  ggplot( aes(x=pop, y=obs_RR, group=repN)) +
  geom_line(alpha=0.1) +
  theme_classic() +
  ylim(0.2, 2)
```

```{r}
obs_data %>%
  filter(ve.new.trial==0)  %>%
  ggplot( aes(x=pop, y=obs_RR, group=repN)) +
  geom_line(alpha=0.1) +
  theme_classic() +
  ylim(0.2, 2)

obs_data %>%
  filter(ve.new.trial==0 & pop==3000 ) %>%
   ggplot( aes(x=N_events, y=obs_RR)) +
  geom_point()+
  theme_classic() 
```

```{r}
obs_exp_data <-   obs_data %>%
  left_join(probs, by=c('repN','pop','ve.new.trial')) 

obs_exp_data %>%
  filter(ve.new.trial==44.4 & prior.info=="N(0,1e-4)") %>%
  ggplot(aes(x=log(obs_RR), y=parm_mean)) +
  geom_point() +
  theme_classic()

obs_exp_data %>%
  filter(ve.new.trial==44.4 & prior.info=="N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(0.01,0.01)") %>%
  ggplot(aes(x=log(obs_RR), y=parm_mean)) +
  geom_point() +
  theme_classic()

```

Probability of success as a function of number of events

```{r}
obs_exp_data %>%
  filter(ve.new.trial==44.4 & prior.info=="N(\u03B4\u2081,\u03C3\u00B2); \u03C3\u00B2~IG(0.01,0.01)") %>%
  mutate(tot_events=placebo+vax) %>%
   ggplot(aes(x=tot_events, y=p_0)) +
  geom_point() +
  theme_classic()

```

```{r}
probs %>%
  filter(ve.new.trial==44.4 ) %>%
  mutate(Pooling_type = as.factor(Pooling_type)) %>%
ggplot( aes(x=pop, y=parm_mean, group=repN, color=Pooling_type)) +
  geom_line( alpha=0.1) +
  theme_classic() +
  ylim(-4,4)+
  facet_wrap(prior.info~Pooling_type)
```