# load library
library(tidyverse)

# Function to create probability distribution of survival data (PFS)
prob_dist <- function(surv, lcl, ucl, ci = 0.95, arm) {

  lci <- log(-log(lcl))
  uci <- log(-log(ucl))
  lsurv <- log(-log(surv))

  z <- qnorm((1 + ci)/2)
  lsd <- abs(lsurv - lci)/z
  usd <- abs(uci - lsurv)/z

  sd <- (lsd + usd)/2
  xs <- log(-log(seq(from = 0, to = 1, length.out = 1000)))

  samp <- dnorm(x = xs, mean = lsurv, sd = sd)

  tibble::tibble(arm = arm, xs = exp(-exp(xs)), samp = samp)
}

# Function to calculate vector of events (no or yes) for a particular survival value
sim_event <- function(surv, lcl, ucl, ci = 0.95, N_sim = 1000) {

  lci <- log(-log(lcl))
  uci <- log(-log(ucl))
  lsurv <- log(-log(surv))

  z <- qnorm((1 + ci)/2)
  lsd <- abs(lsurv - lci)/z
  usd <- abs(uci - lsurv)/z

  sd <- (lsd + usd)/2

  pop_surv <- exp(-exp(rnorm(N_sim, mean = lsurv, sd = sd)))

  events <- vector("numeric", N_sim)

  for (i in seq_len(N_sim)) {
    events[i] <- sample(c(0, 1), size = 1, replace = TRUE,
                        prob = c(pop_surv[i], 1 - pop_surv[i]))
  }

  events
}

# Function which is the main simulation and calculates the proportion of times
# Arm 1 and Arm 2 achieve any event
simulation <- function(s1, lcl1, ucl1, s2, lcl2, ucl2, n_sim) {

  foo <- function() {
    ev1 <- sim_event(s1, lcl1, ucl1, 0.95, 1000)
    ev2 <- sim_event(s2, lcl2, ucl2, 0.95, 1000)

    c(prop_01 = mean(1*((ev1 == 0) & (ev2 == 1))),
      prop_10 = mean(1*((ev1 == 1) & (ev2 == 0))),
      prop_11 = mean(1*((ev1 == 1) & (ev2 == 1))),
      prop_00 = mean(1*((ev1 == 0) & (ev2 == 0))))
  }

  rerun(n_sim, foo())
}

# CODES TO CREATE THE PROBABILITY DISTRIBUTION GRAPH
df1 <- prob_dist(surv = 0.867, lcl = 0.826, ucl = 0.899, arm = "Arm 1")
df2 <- prob_dist(surv = 0.76, lcl = 0.71, ucl = 0.802, arm = "Arm 2")

df <- bind_rows(df1, df2)

ggplot(df, aes(x = xs, y = samp, colour = arm)) +
  geom_line(size = 1.5) + xlim(0.65, 1) + labs(x = "PFS", y = "density", colour = "Arm", title = "There is some overlap between the PFSs") + theme_bw()

# CODES TO CREATE THE TABLE OF PROPORTIONS AND GRAPH OF PROPORTIONS

ss <- simulation(0.867, 0.826, 0.899, 0.76, 0.71, 0.802, n_sim = 500)

df_prop <- tibble::as_tibble(do.call(rbind, ss))
df_prop <- df_prop %>% gather("scenario", "proportion") %>% mutate(scenario =
                                case_when(scenario == "prop_01" ~ "Arm1 event- Arm2 event+",
                                        scenario == "prop_00" ~ "Arm1 event- Arm2 event-",
                                        scenario == "prop_10" ~ "Arm1 event+ Arm2 event-",
                                        scenario == "prop_11" ~ "Arm1 event+ Arm2 event+"))

knitr::kable(df_prop %>%  group_by(scenario) %>%
               summarise(mean = mean(proportion),
                         sd = sd(proportion),
                         minimum = min(proportion),
                         maximum = max(proportion)), digits = 4)

ggplot(df_prop, aes(x = proportion, y = ..density.., fill = scenario)) +
geom_density() + labs(x = "Proportion", y = "Density", colour = "Scenario",
title = "Arm 1 is better than Arm 2 in only around 20% of times.\n
In 80% of times, Arm 1 is either same as Arm 2 or\n worse than Arm 2") +
theme_bw() + xlim(0, 0.9)