##################################################################################################
### File runs a 2 arm fixed trial from example_2arm_Fixed_Ordinal.R ###
### Uses Bayesian ordinal regression as primary analysis, assuming proportional odds effects ###
### Prepared by Liz Lorenzi and Lindsay Berry ####
##################################################################################################

source("example_2arm_Fixed_Ordinal_stanarm.R")

# Alternative options -- these run much slower than the stan_polr function but are more customizable
#source("example_2arm_Fixed_Ordinal_stan.R")     # option to run in Stan
#source("example_2arm_Fixed_Ordinal_JAGS.R")     # option to run with JAGS


### Set up assumptions ###
## Control rate assumptions for 8-point WHO endpoint ##
cont.rates = c(0.02, 0.01, 0.02, 0.07, 0.08, 0.10, 0.30, 0.40) # 1 is worst outcome (death), 8 is best
N = 300       # total sample size
nsims = 100   # number of simulations to run

## Simulate trials under different treatment effects ##
ORs = c(1, 1.5, 2)
sims1 = getPower(maxN = N, n_it = nsims, OR = ORs[1])
sims1.5 = getPower(maxN = N, n_it = nsims, OR = ORs[2])
sims2 = getPower(maxN = N, n_it = nsims, OR = ORs[3])

#Compare power under a few different treatment effects
allPowers = c(sims1$power, sims1.5$power, sims2$power)

plot(ORs, allPowers, type = 'b', pch = 19,
     ylim = c(0, 1),
     xlab = "Odds Ratio", 
     ylab = "Power")
grid()

#Summarize number of deaths on each arm
allSimOutput = rbind(data.frame(sims1$interims, OR=1),
                     data.frame(sims1.5$interims, OR=1.5),
                     data.frame(sims2$interims, OR=2))

allSimOutput %>% group_by(OR) %>%
  summarize(meanDeaths1 = mean(numDeaths1),
            meanDeaths2 = mean(numDeaths2))