R/GNGFunctions.R

In lylyf1987/GNGpkg:

# Package Documentataion --------------------------------------------------
#' @title DSS Bayesian Quantitative Decision Making Focus Group's Go-NoGo
#' @description This packages houses functions used by the BQDM's GNG Shiny Applications
#' @details These function will support Shiny apps with the same goal
#' @author Greg Cicconetti
#' @docType package
#' @name GNGSep2019
NULL

# Helpers ------
#' @title gcurve
#' @param expr expression to be plotted
#' @param from lower bound of x-values
#' @param to upper bound of x-values
#' @param n number of points to plot
#' @param category optional text column appends to data.frame returned
#' @return A data.frame is returned with x and y-values and an optional column called category
#' @examples
#' \dontrun{
#' curve(dnorm(x, mean=0, sd=1), from=-4, to = 4, n= 1001)
#' ggplot(gcurve(expr = dnorm(x, mean=0, sd=1),from=-4, to = 4, n= 1001,
#' category= "Standard Normal"), aes(x=x, y=y)) + geom_line()
#' }
#' @description Returns a data.frame associated with a call to base::curve.
gcurve <- function (expr, from = NULL, to = NULL, n = 101, add = FALSE,
                    type = "l", xname = "x", xlab = xname, ylab = NULL,
                    log = NULL, xlim = NULL, category = NULL, ...)
{
  sexpr <- substitute(expr)
  if (is.name(sexpr)) {
    expr <- call(as.character(sexpr), as.name(xname))
  }
  else {
    if (!((is.call(sexpr) || is.expression(sexpr)) && xname %in%
          all.vars(sexpr)))
      stop(gettextf("'expr' must be a function, or a call or an expression containing '%s'",
                    xname), domain = NA)
    expr <- sexpr
  }
  if (dev.cur() == 1L && !identical(add, FALSE)) {
    warning("'add' will be ignored as there is no existing plot")
    add <- FALSE
  }
  addF <- identical(add, FALSE)
  if (is.null(ylab))
    ylab <- deparse(expr)
  if (is.null(from) || is.null(to)) {
    xl <- if (!is.null(xlim))
      xlim
    else if (!addF) {
      pu <- par("usr")[1L:2L]
      if (par("xaxs") == "r")
        pu <- extendrange(pu, f = -1 / 27)
      if (par("xlog"))
        10  ^  pu
      else pu
    }
    else c(0, 1)
    if (is.null(from))
      from <- xl[1L]
    if (is.null(to))
      to <- xl[2L]
  }
  lg <- if (length(log))
    log
  else if (!addF && par("xlog"))
    "x"
  else ""
  if (length(lg) == 0)
    lg <- ""
  if (grepl("x", lg, fixed = TRUE)) {
    if (from <= 0 || to <= 0)
      stop("'from' and 'to' must be > 0 with log=\"x\"")
    x <- exp(seq.int(log(from), log(to), length.out = n))
  }
  else x <- seq.int(from, to, length.out = n)
  ll <- list(x = x)
  names(ll) <- xname
  y <- eval(expr, envir = ll, enclos = parent.frame())
  for.return <- data.frame(x = x, y = y)
  if (is.null(category) == FALSE)
    for.return$category <- factor(category)
  return(for.return)
}

#' @title pooled.sd
#' @param s expression to be plotted
#' @param n number of points to plot
#' @param category optional text column appends to data.frame returned
#' @return A vector holding the pooled standard deviation is returned.
#' @description Computes to the pooled standard deviation.
#' @examples
#' \dontrun{
#' pooled.sd()
#' }
pooled.sd <- function(s = c(4, 5), n = c(14, 20)){
  return(sqrt(sum(s  ^  2 * (n - 1))/(sum(n) - (length(n) - 1))))}

#' @title Density function for the Normal-Gamma distribution
#' @param mu likelihood evaluated when mean takes on value mu
#' @param tau likelihood evaluated when precision takes on value tau
#' @param mu0 hyperparameter describing mu
#' @param n0 hyperparameter describing effective sample size associated with mu0
#' @param a0 hyperparameter describing shape parameter of precision parameter
#' @param b0 hyperparameter describing rate parameter of precision parameter
#' @return Returns the value of the Normal-gamma density function at the point passed.
#' @examples
#' \dontrun{
#' dnorgam()
#' }
#' @description Density function for the Normal-Gamma distribution
dnorgam <- function(mu, tau, mu0, n0, a0, b0) {
  Zng <- gamma(a0)/b0 ^ (a0) * sqrt(2 * pi / n0)
  den <- 1 / Zng * tau ^ (a0 - 1 / 2) *exp(-tau / 2 * (n0*(mu - mu0) ^ 2 + 2 * b0))
  return(den)
}

#' @title Normal-Gamma Posterior Updating
#' @param mu.0 prior mean
#' @param n.0 prior effective sample size
#' @param alpha.0 prior alpha parameter
#' @param beta.0 prior beta parameter
#' @param xbar observed sampled mean
#' @param s observed sample standard deviation
#' @param n sample size
#' @param group text string for group label
#' @return Returns a data.frame with prior, data, and posterior parameters.
#' @examples
#' \dontrun{
#' get.NG.post()
#' }
#' @description Normal-Gamma Posterior Updating
get.NG.post <- function(mu.0 = 0, n.0 = 10, alpha.0 = .25, beta.0 = 1,
                        xbar = .25, s = 3, n = 15, group = "Control"){
  data.frame(group = group) %>%
    mutate(
      # Prior NG parameters
      mu.0 = mu.0,
      n.0 = n.0,
      alpha.0 = alpha.0,
      beta.0 = beta.0,
      # Marginal parameters of t-distribution describing mu.0
      tdf.0 = 2 * alpha.0,
      sigma.0 = beta.0/(alpha.0 * n.0),
      # data
      xbar = xbar,
      s = s,
      n = n,
      # Posterior parameters
      mu.n = (n.0 * mu.0 + n * xbar)/(n.0 + n),
      n.n = n.0 + n,
      alpha.n = alpha.0 + n / 2,
      beta.n = beta.0 + (n - 1)/2 * s ^ 2 + n.0 * n * (xbar - mu.0) ^ 2/(2 * (n.0 + n)),
      # Marginal parameters of t-distribution describing mu.n
      tdf.n = 2 * alpha.n,
      sigma.n = beta.n/(alpha.n * n.n)
    )
}


# Random normal-gamma generation

rnormgam <- function(n=100000, mu.0 = 0, n.0 = 1, alpha.0 = .01, beta.0 = .01){
  tau <- rgamma(n=n, shape = alpha.0, rate = beta.0)
  x <- rnorm(n=n, mean=mu.0, sd=(1/n.0*tau)^.5)
  return(x)
}


# Reparameterized.normal ------
#' @name Reparameterized.normal
#' @title Reparameterized normal functions
#' @param mean mean
#' @param tau precision parameter
#' @param x likelihood function evaluates at point x
#' @param q quantile
#' @param p cumulative probability
#' @param n sample size
#' @return Returns density values, cumulative probabilities, quantiles and random samples from a normal distribution.
#' @description Reparameterized normal functions
#' @examples
#' \dontrun{
#' dnorm.rp()
#' pnorm.rp()
#' qnorm.rp()
#' rnorm.rp()
#' }
dnorm.rp <- function(x, mean = 0, tau = 1, log = FALSE){
  dnorm(x, mean = mean, sd = 1 / sqrt(tau), log = log)
}

#' @rdname Reparameterized.normal
pnorm.rp <- function(q, mean = 0, tau = 1, lower.tail = TRUE, log.p = FALSE){
  pnorm(q, mean = mean, sd = 1 / sqrt(tau), lower.tail = lower.tail, log.p = log.p)}

#' @rdname Reparameterized.normal
qnorm.rp <- function(p, mean = 0, tau = 1, lower.tail = TRUE, log.p = FALSE){
  qnorm(p, mean = mean, sd = 1 / sqrt(tau), lower.tail = lower.tail, log.p = log.p)}

#' @rdname Reparameterized.normal
rnorm.rp <- function(n, mean = 0, tau = 1){
  rnorm(n, mean = mean, sd = 1 / sqrt(tau))}

# Location.scale.t ------
#' @name Location.scale.t
#' @title Location-scale t distribution functions
#' @param mu mean
#' @param df degrees of freedom
#' @param sigma scale parameter
#' @param x likelihood function evaluates at point x
#' @param q quantile
#' @param p cumulative probability
#' @param n sample size
#' @return Returns density values, cumulative probabilities, quantiles and random samples from a Location-scale t distribution.
#' @examples
#' \dontrun{
#' dt_ls()
#' pt_ls()
#' qt_ls()
#' rt_ls()
#' }
#' @description Location-scale t distribution functions
dt_ls <- function(x, df, mu, sigma) 1 / sigma * dt((x - mu)/sigma, df)

#' @rdname Location.scale.t
pt_ls <- function(x, df, mu, sigma) pt(q = (x - mu)/sigma, df)

#' @rdname Location.scale.t
qt_ls <- function(prob, df, mu, sigma) qt(p = prob, df)*sigma + mu

#' @rdname Location.scale.t
rt_ls <- function(n, df, mu, sigma) rt(n, df)*sigma + mu

# Reparameterized.beta ------
#' @name Reparameterized.beta
#' @title Reparameterized Beta distribution functions
#' @param mean mean
#' @param effective.ss effective sample size
#' @param ncp non-centrality parameter
#' @param log as in stats::*beta
#' @param x likelihood function evaluates at point x
#' @param q quantile
#' @param p cumulative probability
#' @param n sample size
#' @return Returns density values, cumulative probabilities, quantiles and random samples from a Beta distribution.
#' @examples
#' \dontrun{
#' dbeta.rp()
#' pbeta.rp()
#' qbeta.rp()
#' }
#' @description Reparameterized Beta distribution functions
dbeta.rp <- function(x, mean = .5, effective.ss = 1, ncp = 0, log = FALSE){
  dbeta(x = x,
        shape1 = mean * effective.ss,
        shape2 = effective.ss - mean * effective.ss,
        ncp = ncp, log = log)
}

#' @rdname Reparameterized.beta
pbeta.rp <- function(q, mean = .5, effective.ss = 1, ncp = 0, lower.tail = TRUE,
                     log.p = FALSE){
  pbeta(q = q,
        shape1 = mean * effective.ss,
        shape2 = effective.ss - mean * effective.ss,
        ncp = ncp, lower.tail = lower.tail, log.p = log.p)
}

#' @rdname Reparameterized.beta
qbeta.rp <- function(p, mean = .5, effective.ss = 1, ncp = 0, lower.tail = TRUE,
                     log.p = FALSE){
  qbeta(p = p,
        shape1 = mean * effective.ss,
        shape2 = effective.ss - mean * effective.ss,
        ncp = ncp, lower.tail = lower.tail, log.p = log.p)
}

# Reparameterized.gamma ------
#' @name Reparameterized.gamma
#' @title Reparameterized Gamma distribution functions
#' @param mean mean
#' @param var variance
#' @param log as in stats::*gamma
#' @param x likelihood function evaluates at point x
#' @param q quantile
#' @param p cumulative probability
#' @param n sample size
#' @return Returns density values, cumulative probabilities, quantiles and random samples from a gamma distribution.
#' @examples
#' \dontrun{
#' dgamma.rp()
#' pgamma.rp()
#' qgamma.rp()
#' rgamma.rp()
#' }
#' @description Reparameterized Gamma distribution functions
dgamma.rp <- function(x, mean, var, log = FALSE){
  dgamma(x, shape = mean * mean / var, rate = mean / var)}

#' @rdname Reparameterized.gamma
pgamma.rp <- function(q, mean, var, lower.tail = TRUE, log.p = FALSE){
  pgamma(q, shape = mean * mean / var, rate = mean / var, lower.tail = lower.tail,
         log.p = log.p)}

#' @rdname Reparameterized.gamma
qgamma.rp <- function(p, mean, var, lower.tail = TRUE, log.p = FALSE){
  qgamma(p, shape = mean * mean / var, rate = mean / var, lower.tail = lower.tail,
         log.p = log.p)}

#' @rdname Reparameterized.gamma
rgamma.rp <- function(n, mean, var){
  rgamma(n, shape = mean * mean / var, rate = mean / var)}



# SingleSampleBinary ------

#' @name SingleSampleBinary
#' @title Functions to Support the One Sample Binary Scenario
#' @param a.trt prior alpha parameter
#' @param b.trt prior beta parameter
#' @param n.trt observed sample size
#' @param x.trt observed number of responders
#' @param Delta.lrv TPP Lower Reference Value aka Min TPP
#' @param Delta.tv TPP Target Value aka Base TPP
#' @param tau.tv threshold associated with Base TPP
#' @param tau.lrv threshold associated with Min TPP
#' @param tau.ng threshold associated with No-Go
#' @param seed random seed
#' @param nlines Control for text spacing
#' @param tsize Control for text size
#' @param nlines.ria Control for text spacing
#' @param Delta_OC_LB treatment effect lower bound
#' @param Delta_OC_UB treatment effect upper bound
#' @param my.df data.frame returned by get.ss.bin.trt.oc.df
#' @param beta.mean range of values for prior mean
#' @param eff.ss range of values for effective sample size
#' @param ulow lower sample size value
#' @param uupp upper sample size value
#' @param for.plot data.frame from get.ss.bin.ssize.oc.df
#' @examples
#' \dontrun{
#' make.ss.bin.ppp()
#' make.ss.bin.spp()
#' get.ss.bin.trt.oc.df()
#' make.ss.bin.trt.oc1()
#' make.ss.bin.trt.oc2()
#' get.ss.bin.ssize.oc.df()
#' make.ss.bin.ssize.oc()
#' }
#' @description
#'
#'   make.ss.bin.ppp: Make Single Sample Binary Prior/Posterior Plot. Returns a ggplot object.
#'
#'   make.ss.bin.spp: Make Single Sample Binary Shaded Posterior Plot.  Returns a graphic built using grid.arrange.
#'
#'   get.ss.bin.trt.oc.df: Get Single Sample Binary Treatment Effect OC. Returns a data.frame.
#'
#'   make.ss.bin.trt.oc1: Make Single Sample Binary Treatment Effect. Returns a graphic built using grid.arrange.
#'
#'   make.ss.bin.trt.oc2: Make Single Sample Binary Treatment Effect.  Returns a graphic built using grid.arrange.
#'
#'   get.ss.bin.ssize.oc.df:  Get Single Sample Binary sample size OC data.frame. Returns a data.frame.
#'
#'   make.ss.bin.ssize.oc: Make Single Sample Binary Sample size OC plot

make.ss.bin.ppp <- function(a.trt = 1, b.trt = 1, n.trt = 40, x.trt = 20){

  prior.df <- rbind(
    gcurve(expr = dbeta(x, shape1 = a.trt,shape2 = b.trt), from = 0, to = 1, n = 1001,
           category = "Treatment") %>%
      mutate(alpha = a.trt, beta = b.trt, group="Prior", density="Treatment Prior"),
    gcurve(expr = dbeta(x, shape1 = a.trt + x.trt, shape2  = b.trt + n.trt - x.trt),
           from = 0, to = 1, n = 1001, category = "Treatment") %>%
      mutate(alpha = a.trt + x.trt, beta = b.trt + n.trt - x.trt, group="Posterior",
             density="Treatment Posterior")) %>%
    mutate(density = factor(density, c("Treatment Prior", "Treatment Posterior")))
  levels(prior.df$density) <- c(paste0("Treatment Prior: Beta(", a.trt, ", ", b.trt,
                                       ")"),
                                paste0("Treatment Posterior: Beta(", a.trt + x.trt,
                                       ", ",  b.trt + (n.trt - x.trt), ")"))

  # This ggplot call is same for both plots http://127.0.0.1:7505/#tab-1955-1 and downloads too
  ggplot(data = prior.df, aes(x = x,y = y, color = category))+geom_line() +
    labs(x = "Response Rate (%)",
         y = NULL,
         title="Prior and posterior distributions of treatment response rate",
         subtitle = paste0("Treatment data: ",round(x.trt/n.trt,2)*100,
                           "% (",  x.trt, "/", n.trt, ")"),
         color="Group", linetype="Density") +
    facet_wrap(~density)+
    theme(legend.position = "bottom")+guides(color=F)+
    scale_y_continuous(breaks=NULL, label=NULL)+
    scale_x_continuous(labels = scales::percent)
}

#' @rdname SingleSampleBinary
make.ss.bin.spp <- function(a.trt = 1, b.trt = 1, n.trt = 40, x.trt = 9,
                            Delta.lrv = .2, Delta.tv = .35,
                            tau.tv = 0.10, tau.lrv = .80, tau.ng = .65,
                            seed = 1234, nlines = 25, tsize = 4, nlines.ria=20){

  results <- get.ss.bin.df(a.trt = a.trt, b.trt = b.trt, n.trt = n.trt, x.trt = 0:n.trt,
                           Delta.lrv = Delta.lrv, Delta.tv=Delta.tv, tau.lrv=tau.lrv,
                           tau.tv=tau.tv, tau.ng=tau.ng)

  result.go <- results %>% dplyr::filter(result== "Go") %>% slice(1)

  # Determine max number of TRT responders for No-Go
  result.ng <- results %>% dplyr::filter(result== "No-Go") %>% slice(n())
  my.df <- gcurve(expr = dbeta(x,shape1 = a.trt + x.trt, shape2 = b.trt + n.trt - x.trt),
                  from = 0, to =1, category="Posterior") %>%
    mutate(alpha = a.trt + x.trt, beta = b.trt + n.trt - x.trt, group="Posterior")
  my.df$facet <- paste0("Prior: Beta(", a.trt, ", ", b.trt, ")")

  P.R1 = 1 - pbeta(Delta.lrv, a.trt + x.trt, b.trt + (n.trt - x.trt))
  P.R3 = 1 - pbeta(Delta.tv, a.trt + x.trt, b.trt + (n.trt - x.trt))
  result = ifelse(P.R1 >= tau.lrv & P.R3 >= tau.tv, "Go",
                  ifelse(P.R1 < tau.ng  & P.R3 < tau.tv, "No-Go", "Consider"))

  result.color <- ifelse(result=="Go", "darkgreen",
                         ifelse(result=="No-Go", "red", "black"))

  for.subtitle <- paste0("Treatment data: ", round(x.trt/n.trt,2)*100,
                         "% (",  x.trt, "/", n.trt, "). ",
                         "Given control data, responders needed for Go: ",
                         result.go$x[1],
                         " (", round(result.go$x[1]/n.trt,2)*100,
                         "%). Needed for No-Go: ",  result.ng$x[1], " (",
                         round(result.ng$x[1]/n.trt, 2)*100,"%)")

  # Annotation line 1: Decision ----
  for.decision <- paste0("Decision: ", result)

  # Annotation line 2: Decision interval ----
  if(result!="No-Go"){
    for.decision.interval <- paste0("Decision interval: (",
                                    round(qbeta(p = 1 - tau.lrv, shape1 = a.trt + x.trt,
                                                shape2 = b.trt + n.trt - x.trt),3)*100,
                                    "%, ",
                                    round(qbeta(p = 1 - tau.tv, shape1 = a.trt + x.trt,
                                                shape2 = b.trt + n.trt - x.trt),3)*100,
                                    "%)")
  } else {
    for.decision.interval <- paste0("Decision interval: (",
                                    round(qbeta(p = 1 - tau.ng, shape1 = a.trt + x.trt,
                                                shape2 = b.trt + n.trt - x.trt),3)*100, "%, ",
                                    round(qbeta(p = 1 - tau.tv, shape1 = a.trt + x.trt,
                                                shape2 = b.trt + n.trt - x.trt),3)*100, "%)")
  }

  # Annotation line 3: P(Delta >= Min TPP)

  annotate.P1 <- ifelse(
    result=="Go",
    TeX(paste0("$P(\\Delta\\, \\geq$ Min TPP) = $",
               round(1 - pbeta(q = Delta.lrv, shape1 = a.trt + x.trt,
                               shape2 = b.trt + n.trt - x.trt),3)*100,
               "$% > ", tau.lrv*100, "%")),
    ifelse(result=="No-Go",
           TeX(paste0("$P(\\Delta\\, \\geq$ Min TPP) = $",
                      round(1 - pbeta(q = Delta.lrv, shape1 = a.trt + x.trt,
                                      shape2 = b.trt + n.trt - x.trt),3)*100,
                      "$% $\\leq$ ", tau.ng*100, "%")),
           TeX(paste0("$P(\\Delta\\, \\geq$ Min TPP) = $",
                      round(1 - pbeta(q = Delta.lrv, shape1 = a.trt + x.trt,
                                      shape2 = b.trt + n.trt - x.trt),3)*100, "$%"))))

  # Annotation line 4: P(Delta >= Base TPP)
  annotate.P2 <- ifelse(
    result=="Go",
    TeX(paste0("$P(\\Delta\\, \\geq$ Base TPP) = $",
               round(1 - pbeta(q = Delta.tv, shape1 = a.trt + x.trt,
                               shape2 = b.trt + n.trt - x.trt),3)*100, "$% > ",
               tau.tv*100, "%")),
    ifelse(result=="No-Go", TeX(paste0("$P(\\Delta\\, \\geq$ Base TPP) = $",
                                       round(1 - pbeta(q = Delta.tv, shape1 = a.trt +
                                                         x.trt, shape2 = b.trt + n.trt
                                                       - x.trt),3)*100, "$% $\\leq$",
                                       tau.tv*100, "%")),
           TeX(paste0("$P(\\Delta\\, \\geq$ Base TPP) = $",
                      round(1 - pbeta(q = Delta.tv, shape1 = a.trt + x.trt,
                                      shape2 = b.trt + n.trt - x.trt),3)*100, "$%"))))

  # Initialize a ggplot
  dplot <- ggplot() + geom_line(data = my.df, aes(x = x,y = y)) + facet_wrap(~facet)
  # Access the ggplot to get goodies to help accomplish shading
  dpb <- ggplot_build(dplot)
  x1.1 <- min(which(dpb$data[[1]]$x >=Delta.lrv))
  x2.1 <- max(which(dpb$data[[1]]$x <=Delta.tv))+1
  x1.2 <- min(which(dpb$data[[1]]$x >=Delta.tv))
  x2.2 <- max(which(dpb$data[[1]]$x <=1))

  if(which(dpb$data[[1]]$y== max(dpb$data[[1]]$y)) > length(dpb$data[[1]]$x)/2){
    annotate.x = min(min(dpb$data[[1]]$x), Delta.lrv, 0)
    annotate.j = 0
  } else {annotate.x = max(max(dpb$data[[1]]$x), Delta.tv, 1)
  annotate.j = 1}

  go.segment <- data.frame(x = qbeta(p = 1 - tau.lrv, shape1 = a.trt + x.trt,
                                     shape2 = b.trt + n.trt - x.trt),
                           xend = qbeta(p = 1 - tau.tv, shape1 = a.trt + x.trt,
                                        shape2 = b.trt + n.trt - x.trt),
                           y= min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines * 3,
                           yend = min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines * 3,
                           group=my.df$group[1])
  nogo.segment <- data.frame(x = qbeta(p = 1 - tau.ng, shape1 = a.trt + x.trt,
                                       shape2 = b.trt + n.trt - x.trt),
                             xend = qbeta(p = 1 - tau.tv, shape1 = a.trt + x.trt,
                                          shape2 = b.trt + n.trt - x.trt),
                             y= min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines * 3,
                             yend = min(dpb$data[[1]]$y) +
                               max(dpb$data[[1]]$y)/nlines * 3,
                             group=my.df$group[1])

  # Introduce shading
  main.plot <- dplot +
    geom_area(data = data.frame(x = dpb$data[[1]]$x[1:x1.1],
                                y = dpb$data[[1]]$y[1:x1.1]),
              aes(x = x, y = y), fill=alpha("red", 0))+
    geom_area(data = data.frame(x = dpb$data[[1]]$x[x1.1:x2.1],
                                y = dpb$data[[1]]$y[x1.1:x2.1]),
              aes(x = x, y = y), fill=alpha("grey80", 0))+
    geom_area(data = data.frame(x = dpb$data[[1]]$x[x1.2:x2.2],
                                y = dpb$data[[1]]$y[x1.2:x2.2]),
              aes(x = x, y = y), fill=alpha("green", 0))+
    geom_line(data = my.df, aes(x = x,y = y))+
    scale_x_continuous(limits = c(0,1),
                       breaks = unique(c(
                         pretty(x = c(0, Delta.lrv-.06),
                                n = diff(c(0, Delta.lrv))/2 * 10)[
                                  pretty(x = c(0, Delta.lrv-.06),
                                         n = diff(c(0, Delta.lrv))/2 * 10) <
                                    Delta.lrv-.06],
                         c(Delta.lrv, Delta.tv),
                         pretty(x = c(Delta.tv, 1), n = diff(c(Delta.tv,1))/2 * 10)[
                           pretty(x = c(Delta.tv, 1), n = diff(c(Delta.tv,1))/2 * 10) >
                             Delta.tv+.06]
                       )),
                       labels = scales::percent)+
    scale_y_continuous(breaks=NULL, labels = NULL)

  # Add Annotations -----
  main.plot <- main.plot+
    labs(title = TeX("Posterior distribution for the proportion of treatment responders"),
         subtitle = for.subtitle,
         x = TeX("$\\Delta\\,$ = Proportion of treatment responders"), y = NULL,
         caption = NULL)+
    annotate("text", label = for.decision,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 0,
             size = tsize+1, colour = result.color, hjust = annotate.j)+
    annotate("text", label = for.decision.interval,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 1,
             size = tsize, colour = result.color, hjust = annotate.j)+
    annotate("text", label = annotate.P1, color=result.color,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 2,
             size = tsize,  hjust = annotate.j)+
    annotate("text", label = annotate.P2, color=result.color,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 3,
             size = tsize,  hjust = annotate.j)

  # Add reference lines and Credible interval

  if(result.color == "red"){
    main.plot <- main.plot +
      geom_segment(data=nogo.segment, aes(x=x, xend=xend, y=y, yend=yend, group=group),
                   arrow = arrow(ends="both", angle = 90), color=result.color, size=.75)
  } else {
    main.plot <- main.plot +
      geom_segment(data=go.segment, aes(x=x, xend=xend, y=y, yend=yend, group=group),
                   arrow = arrow(ends="both", angle = 90), color=result.color, size=.75)}
  main.plot <- main.plot +
    geom_vline(xintercept = c(Delta.lrv, Delta.tv), linetype = 2,
               color = c("blue", "blue")) +
    annotate("label", label = TeX(paste0("Min TPP = $", Delta.lrv*100,"$%")),
             x = Delta.lrv, y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize,
             colour = "black", hjust = 1, label.size=NA, fill="grey92", alpha=.8)+
    annotate("text", label = TeX(paste0("Base TPP = $",  Delta.tv*100,"$%")),
             x = Delta.tv, y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize,
             colour = "black", hjust = 0)

  table.plot2 <- ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) > ", tau.lrv*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) > ", tau.tv*100,"%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) $\\leq$ $", tau.ng*100, "$$%\\;\\;\\;\\;$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) $\\leq$ $", tau.tv*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")

  # plot_grid(main.plot,
  #            table.plot2, ncol=1, align = "v", rel_heights=c(.78,.22), axis="l")
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
}

#' @rdname SingleSampleBinary
get.ss.bin.trt.oc.df <- function(a.trt = 1, b.trt = 1, n.trt = 40,
                                 Delta.tv = 0.35, Delta.lrv = 0.2,
                                 Delta.OC.LB= 0, Delta.OC.UB=.5,
                                 tau.tv = 0.1, tau.lrv = 0.8, tau.ng = 0.65) {
  my.df <- expand.grid(successes=0:n.trt, Delta=seq(0,1,.01)) %>%
    mutate(prob=dbinom(successes, n.trt, Delta)) %>%
    mutate(a.post = a.trt + successes, b.post = b.trt + n.trt - successes) %>%
    mutate(P.R1 = 1 - pbeta(Delta.lrv, a.post, b.post),
           P.R3 = 1 - pbeta(Delta.tv, a.post, b.post),
           result = ifelse(P.R1 >= tau.lrv & P.R3 >= tau.tv, "Go",
                           ifelse(P.R1 < tau.ng & P.R3 < tau.tv, "No-Go",
                                  "Consider"))) %>%
    mutate(result = factor(result, c("Go", "Consider", "No-Go"))) %>%
    group_by(Delta, result, .drop=FALSE) %>%
    summarize(p=sum(prob)) %>%
    mutate(a.trt = a.trt, b.trt = b.trt, n.trt = n.trt,
           Delta.tv = Delta.tv, Delta.lrv = Delta.lrv,
           tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng)

  return(my.df)
}

#' @rdname SingleSampleBinary
make.ss.bin.trt.oc1 <- function(my.df = get.ss.bin.trt.oc.df(),
                                tsize=4, nlines=25, nlines.ria=20,
                                Delta_OC_LB = 0, Delta_OC_UB = .41){
  Delta.tv <- my.df$Delta.tv[1]
  Delta.lrv <- my.df$Delta.lrv[1]
  tau.tv <- my.df$tau.tv[1]
  tau.lrv <- my.df$tau.lrv[1]
  tau.ng <- my.df$tau.ng[1]
  my.df2 <- subset(my.df, result=="Go") %>%
    dplyr::filter(Delta <= Delta_OC_UB & Delta >= Delta_OC_LB)
  my.df3 <- subset(my.df, result=="No-Go") %>%
    dplyr::filter(Delta <= Delta_OC_UB & Delta >= Delta_OC_LB)
  # In the case where No-Go is never attained...
  if(nrow(my.df3) ==0){
    my.df3 <- my.df2
    my.df3$result <- "No-Go"
    my.df3$p <- 0
  }
  my.df3$p <- 1 - my.df3$p

  total.n <- my.df$n.trt[1]

  my.plot <- ggplot()+geom_line(data=my.df2, aes(x=Delta, y=p))+
    geom_line(data=my.df3, aes(x=Delta, y=p))
  dpb <- ggplot_build(my.plot)
  main.plot <- my.plot +
    geom_area(data = dpb$data[[1]], aes(x = x, y = y), fill=alpha("lightgreen",0.5))+
    geom_ribbon(data = dpb$data[[1]], aes(x = x, ymin = y, ymax=1),
                fill=alpha("grey", .5))+
    geom_ribbon(data = dpb$data[[2]], aes(x = x, ymin = y, ymax = 1),
                fill=alpha("red", 0.5))+
    geom_line(data=my.df2, aes(x=Delta, y=p), size=.75, color="green")+
    geom_line(data=my.df3, aes(x=Delta, y=p), size=.75, color="red")+
    geom_vline(xintercept = c(Delta.tv, Delta.lrv), linetype=2, color="blue")+
    annotate("text", label = TeX(paste0("Min TPP = $", Delta.lrv*100,"%$")),
             x = Delta.lrv, y = 0 + max(dpb$data[[1]]$y)/nlines.ria,
             size = tsize, colour = "black", hjust = 1)+
    annotate("text", label = TeX(paste0("Base TPP = $",  Delta.tv*100,"%$")),
             x = Delta.tv, y = 0 + max(dpb$data[[1]]$y)/nlines.ria,
             size = tsize, colour = "black", hjust = 0)+
    labs(x = TeX("$\\Delta\\,$ = Underlying proportion of treatment responders"),
         y = "Probability",
         title = "Operating characteristics as a function of treatment effect",
         subtitle = paste0("Total randomized to treatment: ", total.n))+
    scale_x_continuous(expand = c(0, 0), limits=c(Delta_OC_LB,
                                                  min(Delta_OC_UB, max(dpb$data[[1]]$x))),
                       label=scales::percent)+
    scale_y_continuous(expand = c(0,0), breaks=seq(0,1,.2), label=scales::percent)+
    theme(panel.spacing.x = unit(6, "mm"),
          axis.text.x = element_text(angle=45, hjust=1,vjust=1))


  table.plot2 <- ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) > ", tau.lrv*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) > ", tau.tv*100,"%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) $\\leq$ $", tau.ng*100, "$$%\\;\\;\\;\\;$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) $\\leq$ $", tau.tv*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")

  # plot_grid(main.plot,
  #            table.plot2, ncol=1, align = "v", rel_heights=c(.78,.22), axis="l")
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
}


#' @rdname SingleSampleBinary
make.ss.bin.trt.oc2 <- function(my.df = get.ss.bin.trt.oc.df(), tsize=4, nlines=25,
                                nlines.ria=20, ud.low = 0, ud.upp = 1){

  Delta.tv <- my.df$Delta.tv[1]
  Delta.lrv <- my.df$Delta.lrv[1]
  tau.tv <- my.df$tau.tv[1]
  tau.lrv <- my.df$tau.lrv[1]
  tau.ng <- my.df$tau.ng[1]
  my.df$result <- factor(my.df$result, c("Go", "Consider", "No-Go"))
  total.n <- my.df$n.trt[1]

  main.plot <- ggplot(data=my.df, aes(x=Delta, y=p, color=result))+
    geom_line()
  dpb <- ggplot_build(main.plot)
  main.plot <- main.plot +
    scale_color_manual(values=c("green", "grey50", "red"))+
    scale_x_continuous(expand = c(0,0), breaks=seq(0,1,.1), limits=c(ud.low, ud.upp)) +
    scale_y_continuous(expand = c(0,0), breaks=seq(0,1,.2),
                       limits=c(0,1), label=scales::percent)+
    theme(panel.spacing.x = unit(6, "mm"),
          axis.text.x = element_text(angle=45, hjust=1,vjust=1))+
    geom_vline(xintercept = c(Delta.tv, Delta.lrv), color="blue", linetype=2)+
    annotate("text", label = TeX(paste0("Min TPP = $", Delta.lrv*100,"%$")),
             x = Delta.lrv, y = 0 + max(dpb$data[[1]]$y)/nlines.ria, size = tsize,
             colour = "black", hjust = 1)+
    annotate("text", label = TeX(paste0("Base TPP = $",  Delta.tv*100,"%$")),
             x = Delta.tv, y = 0 + max(dpb$data[[1]]$y)/nlines.ria, size = tsize,
             colour = "black", hjust = 0)+
    labs(x = TeX("$\\Delta\\,$ = Underlying proportion of treatment responders"),
         y = "Probability",
         title = "Operating characteristics as a function of treatment effect",
         color = "Decision",
         subtitle = paste0("Total randomized to treatment: ", total.n)) +
    theme(legend.position = "bottom")

  table.plot2 <- ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) > ", tau.lrv*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) > ", tau.tv*100,"%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) $\\leq$ $", tau.ng*100, "$$%\\;\\;\\;\\;$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) $\\leq$ $", tau.tv*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
}

#' @rdname SingleSampleBinary
get.ss.bin.df = function(
  a.trt = seq(0.5,1,0.5), b.trt = seq(0.5,1,0.5),
  beta.mean = seq(0.3,0.7,.01), eff.ss = 1:40,
  x.trt = 0:80, n.trt = c(40:80),
  Delta.tv = .4, Delta.lrv = .3,
  tau.tv = 0.10, tau.lrv = .80, tau.ng = .65, rp = FALSE)
{
  if(rp==FALSE) {
    my.grid <- expand.grid(n = n.trt,
                           x = x.trt,
                           a.prior = a.trt,
                           b.prior = b.trt,
                           Delta.tv = Delta.tv,
                           Delta.lrv = Delta.lrv,
                           tau.tv = tau.tv,
                           tau.lrv = tau.lrv,
                           tau.ng = tau.ng) %>%
      mutate(beta.mean = a.prior / (a.prior + b.prior),
             eff.ss = (a.prior + b.prior))
  } else {
    my.grid <- expand.grid(n = n.trt,
                           x = x.trt,
                           beta.mean = beta.mean,
                           eff.ss = eff.ss,
                           Delta.tv = Delta.tv,
                           Delta.lrv = Delta.lrv,
                           tau.tv = tau.tv,
                           tau.lrv = tau.lrv,
                           tau.ng = tau.ng) %>%
      mutate(a.prior= beta.mean * eff.ss,
             b.prior= eff.ss - beta.mean * eff.ss)
  }

  my.grid <- my.grid %>%
    mutate(a.post = a.prior + x, b.post = b.prior + (n - x)) %>%  # Get posterior parameters
    mutate(
      P.R1 = 1 - pbeta(Delta.lrv, a.post, b.post),
      P.R3 = 1 - pbeta(Delta.tv, a.post, b.post),

      result = ifelse(P.R1 >= tau.lrv & P.R3 >= tau.tv, "Go",
                      ifelse(P.R1 < tau.ng  & P.R3 < tau.tv, "No-Go", "Consider")))

  my.grid$result <- factor(my.grid$result,
                           c("Go", "Consider", "No-Go",
                             "NA: Successes > Subjects"))
  my.grid$result[is.na(my.grid$result)==T] <- "NA: Successes > Subjects"
  return(my.grid)
}

#' @rdname SingleSampleBinary
get.ss.bin.ssize.oc.df <- function(a.trt = 1, b.trt = 1, n.trt = 40,
                                   Delta.lrv = .35, Delta.tv=.35, Delta.user = .40,
                                   tau.tv = 0.10, tau.lrv = 0.8, tau.ng = 0.65,
                                   ulow = 20, uupp = 60){

  # Treatment effect as entered, under null, Base, and Min
  # For a spectrum of sample sizes combination grab the Go/NoGo probabilities
  # For each row we'll run

  specs <- expand.grid(n.trt=ulow:uupp) %>%
    mutate(a.trt = a.trt, b.trt = b.trt,
           Delta.lrv = Delta.lrv, Delta.tv=Delta.tv, Delta.user=Delta.user,
           tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng)

  # These are custom functions that randomly sample the treatment arm's data
  # Placebo prior and data will be AS ENTERED BY USER
  # Treatment data is sampled from binomial distributons with priors provided by
  # user and prob = DATA ENTERED BY USER

  check <- bind_rows(apply(X = matrix(1:nrow(specs)), MARGIN = 1, FUN = function(x) {
    get.ss.bin.df(a.trt = specs$a.trt[x], b.trt = specs$b.trt[x], n.trt = specs$n.trt[x],
                  x.trt = 0:specs$n.trt[x], Delta.tv = specs$Delta.tv[x],
                  Delta.lrv = specs$Delta.lrv[x],
                  tau.tv = specs$tau.tv[x], tau.lrv = specs$tau.lrv[x],
                  tau.ng = specs$tau.ng[x])})) %>% mutate(
                    p.x.trt.null = dbinom(x = x, size = n, prob = 0),
                    p.x.trt.lrv = dbinom(x = x, size = n, prob = Delta.lrv),
                    p.x.trt.tv = dbinom(x = x, size = n, prob = Delta.tv),
                    p.x.trt.user = dbinom(x = x, size = n, prob = Delta.user))

  check <- check %>% group_by(n, result) %>% summarize(
    s.null = sum(p.x.trt.null),
    s.lrv = sum(p.x.trt.lrv),
    s.tv = sum(p.x.trt.tv),
    s.user = sum(p.x.trt.user)) %>%
    gather(value = value, key=key, -c(n,result), factor_key = T) %>%
    dplyr::filter(result!="Consider") %>%
    mutate(value = ifelse(result=="No-Go", 1-value, value))

  levels(check$key) <- c(
    TeX(paste0("$\\Delta\\,$ = Null = ", round(0,2)*100,'%')),
    TeX(paste0("$\\Delta\\,$ = Min TPP = ", round(Delta.lrv,2)*100,'%')),
    TeX(paste0("$\\Delta\\,$ = Base TPP = ", round(Delta.tv,2)*100,'%')),
    TeX(paste0("$\\Delta\\,$ = User defined = ", round(Delta.user,2)*100,'%')))
  check$key <- factor(check$key, levels(check$key)[order(c(0, Delta.lrv, Delta.tv,
                                                           Delta.user))])

  check <- check %>% mutate(a.trt = a.trt, b.trt = b.trt, n.trt = n.trt,
                            Delta.lrv = Delta.lrv, Delta.tv=Delta.tv, Delta.user = Delta.user,
                            tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng)
  return(check)
}

#' @rdname SingleSampleBinary
make.ss.bin.ssize.oc <- function(for.plot=get.ss.bin.ssize.oc.df(), nlines=25, tsize=4){
  Delta.lrv <- for.plot$Delta.lrv[1]
  Delta.tv <- for.plot$Delta.tv[1]
  Delta.user <- for.plot$Delta.user[1]
  tau.tv <- for.plot$tau.tv[1]
  tau.lrv <- for.plot$tau.lrv[1]
  tau.ng <- for.plot$tau.ng[1]
  # Simply takes a dataframe from get.ss.bin.ssize.oc.df and plots
  main.plot <-  ggplot() +
    geom_line(data=for.plot %>% dplyr::filter(result=="Go"),
              aes(x=n, y=value, group=result), color="green") +
    geom_line(data=for.plot %>% dplyr::filter(result=="No-Go"),
              aes(x=n, y=value, group=result), color="red")+
    facet_wrap(~key, labeller="label_parsed")+
    geom_ribbon(data=for.plot %>% dplyr::filter(result=="Go"),
                aes(x=n, ymin=0, ymax=value), fill="green", alpha=.5)+
    geom_ribbon(data=for.plot %>% dplyr::filter(result=="Go"),
                aes(x=n, ymin=value, ymax=1), fill="grey", alpha=.5)+
    geom_ribbon(data=for.plot %>% dplyr::filter(result=="No-Go"),
                aes(x=n, ymin=value, ymax=1), fill="red", alpha=.5)+
    labs(x="Sample size per arm", y="Decision Probabilities",
         title="Operating characteristics as a function of sample size")+
    scale_x_continuous(expand=c(0,0))+
    theme(panel.spacing = unit(2, "lines")) +
    scale_y_continuous(expand = c(0,0), breaks=seq(0,1,.2), limits=c(0,1.01),
                       label=scales::percent)+
    theme(panel.spacing.x = unit(6, "mm"),
          axis.text.x = element_text(angle=45, hjust=1,vjust=1))

  table.plot2 <- ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) > ", tau.lrv*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) > ", tau.tv*100,"%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) $\\leq$ $", tau.ng*100, "$$%\\;\\;\\;\\;$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) $\\leq$ $", tau.tv*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")

  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
}

# SverdlovFunctions ------
## SUPPLEMENT for the paper "Exact Bayesian Inference Comparing Binomial Proportions,
##                           with Application to Proof of Concept Clinical Trials"
## R CODE
## Probability Density Function (PDF)
## Cumulative Distribution Function (CDF)
## Quantiles of the Distribution
## Random Generator
#  of the RELATION of two independent betas
#  RELATION means:
# -- ('DIFF') risk difference (X2 - X1)
# -- ('RR') relative risk (X2/X1)
# -- ('OR') odds ratio ((X2/(1 - X2))/(X1/(1 - X1)))

## Auxiliary functions used to calculate distribution functions and options

# 1) Appell's hypergeometric function:
#     appel.hypgeom(a, b1, b2, c, x, y) = integral from 0 to 1 of f(a, b1, b2, c, x, y) by dt
# where f(a, b1, b2, c, x, y) = 1/beta(a,c-a)*t^(a-1)*(1-t)^(c-a-1)*(1-t*x)^(-b1)*(1-t*y)^(-b2)

#' @name SverdlovFunction
#' @title Sverdlov's functions
#' @param a a
#' @param b b
#' @param a1 a1
#' @param b1 b1
#' @param a2 a2
#' @param b2 b2
#' @param c c
#' @param x x
#' @param y y
#' @param t t
#' @param fun fun
#' @param left left
#' @param right right
#' @param tol tol
#' @param relation relation
#' @param approach approach
#' @param method method
#' @param n n
#' @param left0 left0
#' @param right0 right0
#' @param alpha alpha
#' @description Sverdlov's f function
f <- function(a, b1, b2, c, x, y, t){
  return(gamma(c)/gamma(a)/gamma(c-a)*t^(a-1)*(1-t)^(c-a-1)*(1-t*x)^(-b1)*(1-t*y)^(-b2))
}

#' @rdname SverdlovFunction
appel.hypgeom <- function(a, b1, b2, c, x, y){
  res <- rep(0, length(x))
  for (r in 1:length(res)){
    res[r] <- integrate(function(t) f(a, b1, b2, c, x[r], y[r], t), 0, 1)$value
  }
  return(res)
}

#' @rdname SverdlovFunction
g <- function(a, b, c, x, t){
  return(1/beta(a, c-a)*t^(a-1)*(1-t)^(c-a-1)*(1-t*x)^(-b))
}

#' @rdname SverdlovFunction
gauss.hypgeom <- function(a, b, c, x){
  res <- rep(0, length(x))
  for (r in 1:length(res)){
    res[r] <- integrate(function(t) g(a, b, c, x[r], t), 0, 1)$value
  }
  return(res)
}

#' @rdname SverdlovFunction
root <- function(fun, left, right, tol){
  rt <- 0.5*(left + right)
  while (abs(fun(rt)) > tol){
    left <- rt*(fun(left)*fun(rt) > 0) + left*(fun(left)*fun(rt) < 0)
    right <- rt*(fun(right)*fun(rt) > 0) + right*(fun(right)*fun(rt) < 0)
    rt <- 0.5*(left + right)
  }
  return(rt)
}

#' @rdname SverdlovFunction
r2beta <- function(relation, n, a1, b1, a2, b2){

  set.seed(1)
  X1 <- rbeta(n, shape1 = a1, shape2 = b1)
  X2 <- rbeta(n, shape1 = a2, shape2 = b2)
  # difference
  if (relation == 'DIFF'){
    rnd <- X2 - X1
  }
  # relative risk
  else if (relation == 'RR'){
    rnd <- X2/X1
  }
  # odds ratio
  else if (relation == 'OR'){
    rnd <- (X2/(1-X2))/(X1/(1-X1))
  }
  return(rnd)
}

#' @rdname SverdlovFunction
d2beta <- function(relation, x, a1, b1, a2, b2){
  pdf <- rep(0, length(x))
  # difference
  if (relation == 'DIFF'){
    x.neg <- x[x <= 0]
    if (length(x.neg) != 0){
      pdf[x <= 0] <- beta(a2,b1)/beta(a1,b1)/beta(a2,b2)*(-x.neg)^(b1+b2-1)*
        (1+x.neg)^(a2+b1-1)*
        appel.hypgeom(b1, a1+a2+b1+b2-2, 1-a1, a2+b1, 1+x.neg, 1-x.neg^2)
    }
    x.pos <- x[x > 0]
    if (length(x.pos) != 0){
      pdf[x > 0] <- beta(a1,b2)/beta(a1,b1)/beta(a2,b2)*(x.pos)^(b1+b2-1)*
        (1-x.pos)^(a1+b2-1)*
        appel.hypgeom(b2, a1+a2+b1+b2-2, 1-a2, a1+b2, 1-x.pos, 1-x.pos^2)
    }
  }
  # relative risk
  else if (relation == 'RR'){
    x.01 <- x[x <= 1]
    if (length(x.01) != 0){
      pdf[x <= 1] <- beta(a1+a2,b1)/beta(a1,b1)/beta(a2,b2)*(x.01)^(a2-1)*
        gauss.hypgeom(a1+a2, 1-b2, a1+a2+b1, x.01)
    }
    x.1 <- x[x > 1]
    if (length(x.1) != 0){
      pdf[x > 1] <- beta(a1+a2,b2)/beta(a1,b1)/beta(a2,b2)*(x.1)^(-a1-1)*
        gauss.hypgeom(a1+a2, 1-b1, a1+a2+b2, 1/x.1)
    }
  }
  # odds ration
  else if (relation == 'OR'){
    x.01 <- x[x <= 1]
    if (length(x.01) != 0){
      pdf[x <= 1] <- beta(a1+a2,b1+b2)/beta(a1,b1)/beta(a2,b2)*(x.01)^(a2-1)*
        gauss.hypgeom(a2+b2, a1+a2, a1+b1+a2+b2, 1-x.01)
    }
    x.1 <- x[x > 1]
    if (length(x.1) != 0){
      pdf[x > 1] <- beta(a1+a2,b1+b2)/beta(a1,b1)/beta(a2,b2)*(x.1)^(-b2-1)*
        gauss.hypgeom(a2+b2, b1+b2, a1+b1+a2+b2, 1-1/x.1)
    }
  }
  return(pdf)
}


#' @rdname SverdlovFunction
p2beta <- function(relation, approach, x, a1, b1, a2, b2, n = 1000000){
  cdf <- rep(0, length(x))
  for (r in 1:length(cdf)){
    if (approach == 'SIMULATION'){
      rnd <- r2beta(relation, n, a1, b1, a2, b2)
      cdf[r] <- sum(rnd < x[r])/length(rnd)
    }
    else if (approach == 'DIRECT'){
      # difference    
      if (relation == 'DIFF'){
        if (x[r] < -1) {
          cdf[r] <- 0
        }
        else if ((x[r] > -1) & (x[r] <= 0)) {
          cdf[r] <- integrate(function(t) pbeta(x[r]+t, a2, b2)*dbeta(t, a1, b1), -x[r], 1)$value
        }
        else if ((x[r] > 0) & (x[r] <= 1)) {
          cdf[r] <- integrate(function(t) pbeta(x[r]+t, a2, b2)*dbeta(t, a1, b1), 0, 1 - x[r])$value + 
            integrate(function(t) dbeta(t, a1, b1), 1 - x[r], 1)$value
        }
        else if (x[r] > 1) {
          cdf[r] <- 1
        }
      }
      # relative risk  
      else if (relation == 'RR'){
        if (x[r] < 0) {
          cdf[r] <- 0
        }
        else if ((x[r] >= 0) & (x[r] <= 1)) {
          cdf[r] <- integrate(function(t) pbeta(x[r]*t, a2, b2)*dbeta(t, a1, b1), 0, 1)$value
        }
        else if (x[r] > 1) {
          cdf[r] <- integrate(function(t) pbeta(x[r]*t, a2, b2)*dbeta(t, a1, b1), 0, 1/x[r])$value + 
            integrate(function(t) dbeta(t, a1, b1), 1/x[r], 1)$value
        }
      }
      # odds ratio  
      else if (relation == 'OR'){
        if (x[r] < 0) {
          cdf[r] <- 0
        }
        else {
          cdf[r] <- integrate(function(t) pf(a1*b2/a2/b1*x[r]*t, 2*a2, 2*b2)*df(t, 2*a1, 2*b1), 0, Inf)$value
        }
      }
    }
  }
  return(cdf)
}

#' @rdname SverdlovFunction
q2beta <- function(relation, a1, b1, a2, b2, alpha, tol = 10^(-5)){
  quantile <- rep(0, length(alpha))
  for (r in 1:length(quantile)){
    if (relation == 'DIFF'){
      fun <- function(t){p2beta(relation, approach = 'DIRECT', t, a1, b1, a2, b2) - alpha[r]}
      quantile[r] <- root(fun, -0.9999, 0.9999, tol)
    }
    else {
      fun <- function(x){p2beta(relation, approach = 'DIRECT', tan(pi*x/2), a1, b1, a2, b2) - alpha[r]}
      x <- root(fun, -0.9999, 0.9999, tol)
      quantile[r] <- tan(pi*x/2)
    }
    
  }
  return(quantile)
}

## Credible interval based on Nelder-Mead algorithm or inverse CDF approach
#' @rdname SverdlovFunction
ci2beta <- function(relation, method, a1, b1, a2, b2, alpha, left0, right0){
  if (method == 'neldermead'){
    optim.fun <- function(x) {
      abs(p2beta(relation, approach = 'DIRECT', x[2], a1, b1, a2, b2) -
            p2beta(relation, approach = 'DIRECT', x[1], a1, b1, a2, b2) -
            1 + alpha) +
        abs(d2beta(relation, x[2], a1, b1, a2, b2) -
              d2beta(relation, x[1], a1, b1, a2, b2))
    }
    ci <- optim(c(left0, right0), optim.fun)$par
  }
  else if (method == 'inv.cdf'){
    ci <- q2beta(relation, a1, b1, a2, b2, c(alpha/2, 1 - alpha/2))
  }
  return(ci)
}

# TwoSampleBinary ------
#' @name TwoSampleBinary
#' @title Functions to Support the Two Sample Binary Scenario
#' @param a.con prior alpha parameter for control group
#' @param b.con prior beta parameter for control group
#' @param n.con sample size for control group
#' @param x.con number of responders for control group
#' @param a.trt prior alpha parameter for treatment group
#' @param b.trt prior beta parameter for treatment group
#' @param n.trt sample size for control treatment group
#' @param x.trt number of responders for treatment group
#' @param Delta.lrv TPP Lower Reference Value aka Min TPP
#' @param Delta.tv TPP Target Value aka Base TPP
#' @param tau.tv threshold associated with Base TPP
#' @param tau.lrv threshold associated with Min TPP
#' @param tau.ng threshold associated with No-Go
#' @param seed random seed
#' @param nlines Control for text spacing
#' @param tsize Control for text size
#' @param nlines.ria Control for text spacing
#' @param dcurve.con response rate assumed for control group
#' @param b.trt prior beta parameter for treatment group
#' @param TE.OC.N total sample size for treatment effect OC
#' @param Aratio Allocation ratio
#' @param SS.OC.N.LB sample size lower bound
#' @param SS.OC.N.UB sample size upper bound
#' @param SS.OC.Delta user's TPP
#' @param npoints number of points for sample size OC curve
#' @examples
#' \dontrun{
#' make.ts.bin.ppp()
#' make.ts.bin.spp()
#' get.ts.bin.trt.oc.df()
#' make.ts.bin.trt.oc1()
#' make.ts.bin.trt.oc2()
#' get.ts.bin.ssize.oc.df()
#' make.ts.bin.ssize.oc()
#' }
#' @description
#'
#'   make.ts.bin.ppp: Make Two Sample Binary Prior/Posterior Plot. Returns a ggplot object.
#'
#'   make.ts.bin.spp: Make Two Sample Binary Shaded Posterior Plot.  Returns a graphic built using grid.arrange.
#'
#'   get.ts.bin.trt.oc.df: Get Two Sample Binary Treatment Effect OC. Returns a data.frame.
#'
#'   make.ts.bin.trt.oc1: Make Two Sample Binary Treatment Effect. Returns a graphic built using grid.arrange.
#'
#'   make.ts.bin.trt.oc2: Make Two Sample Binary Treatment Effect.  Returns a graphic built using grid.arrange.
#'
#'   get.ts.bin.ssize.oc.df:  Get Two Sample Binary sample size OC data.frame. Returns a data.frame.
#'
#'   make.ts.bin.ssize.oc: Make Two Sample Binary Sample size OC plot
make.ts.bin.ppp <- function (a.con = 1, b.con = 1, n.con = 40, x.con = 5,
                             a.trt = 1, b.trt = 1, n.trt = 40, x.trt = 20)
{
  prior.df <- rbind(gcurve(expr = dbeta(x, shape1 = a.con, shape2 = b.con),
                           from = 0, to = 1, n = 1001, category = "Control") %>%
                      mutate(alpha = a.con, beta = b.con, group = "Prior",
                             density = "Control Prior"),
                    gcurve(expr = dbeta(x, shape1 = a.trt, shape2 = b.trt),
                           from = 0, to = 1, n = 1001, category = "Treatment") %>%
                      mutate(alpha = a.trt, beta = b.trt, group = "Prior",
                             density = "Treatment Prior"),
                    gcurve(expr = dbeta(x, shape1 = a.con + x.con,
                                        shape2 = b.con + n.con - x.con),
                           from = 0, to = 1, n = 1001, category = "Control") %>%
                      mutate(alpha = a.con + x.con, beta = b.con + n.con - x.con,
                             group = "Posterior",
                             density = "Control Posterior"),
                    gcurve(expr = dbeta(x, shape1 = a.trt + x.trt,
                                        shape2 = b.trt + n.trt - x.trt),
                           from = 0, to = 1, n = 1001, category = "Treatment") %>%
                      mutate(alpha = a.trt + x.trt, beta = b.trt + n.trt - x.trt,
                             group = "Posterior",
                             density = "Treatment Posterior")) %>%
    mutate(group = factor(group, c("Prior", "Posterior"))) %>%
    mutate(category = factor(category, c("Treatment", "Control"))) %>%
    mutate(density = factor(density, c("Control Prior", "Control Posterior",
                                       "Treatment Prior", "Treatment Posterior")))
  levels(prior.df$density) <- c(paste0("Control Prior: Beta(",
                                       a.con, ", ", b.con, ")"),
                                paste0("Control Posterior: Beta(",
                                       a.con + x.con, ", ", b.con + (n.con - x.con), ")"),
                                paste0("Treatment Prior: Beta(",
                                       a.trt, ", ", b.trt, ")"),
                                paste0("Treatment Posterior: Beta(",a.trt + x.trt,
                                       ", ", b.trt + (n.trt - x.trt), ")"))
  levels(prior.df$group) <- c(paste0("Control Prior: Beta(", a.con, ", ",
                                     b.con, ")\nTreatment Prior: Beta(", a.trt,
                                     ", ", b.trt, ")"),
                              paste0("Control Posterior: Beta(", a.con + x.con, ", ",
                                     b.con + (n.con - x.con),
                                     ")\nTreatment Posterior: Beta(", a.trt + x.trt,
                                     ", ", b.trt + (n.trt - x.trt), ")"
                              )
  )

  ggplot(data = prior.df, aes(x = x, y = y, color = category, linetype=category)) +
    geom_line(size = 0.75) + facet_wrap(~group) +
    labs(x = "Response Rate (%)",
         y = NULL,
         color="Group",
         linetype="Group",
         title = "Prior and posterior distributions for response rates",
         subtitle = paste0("Control data: ", round(x.con/n.con,2)*100,
                           "% (", x.con, "/", n.con, "). ",
                           "Treatment data: ", round(x.trt/n.trt,2)*100,
                           "% (", x.trt, "/", n.trt, "). ",
                           "Observed difference: ", (round(x.trt/n.trt,2)
                                                     - round(x.con/n.con,2))*100,"%"))+
    scale_linetype_manual(values=c("solid", "dashed"))+
    scale_y_continuous(breaks=NULL)+
    scale_x_continuous(labels = scales::percent)
}

#' @rdname TwoSampleBinary
get.ts.bin.decision = function(a.con = 1, b.con = 1, n.con = 40, x.con = 5,
                               a.trt = 1, b.trt = 1, n.trt = 40, x.trt = 20,
                               Delta.tv = 0.25, Delta.lrv = 0.2,
                               tau.tv = 0.10, tau.lrv = 0.8, tau.ng = 0.65)
{
  # Reporting P(Delta > Min.tpp), P(Delta > Base.tpp)
  probs <- 1 - p2beta(relation="DIFF", approach="DIRECT", x=c(Delta.lrv, Delta.tv), 
                      a1=a.con+x.con, b1=b.con+n.con-x.con, 
                      a2=a.trt+x.trt, b2=b.trt+n.trt-x.trt)
  # NOW COMPARE HERE
  my.df <- data.frame(P.R1 = probs[1], P.R3 = probs[2]) %>%
    mutate(result = ifelse(P.R1 >= tau.lrv & P.R3 >= tau.tv, "Go",
                           ifelse(P.R1 < tau.ng & P.R3 < tau.tv, "No-Go", "Consider")))
  return(my.df)
}


#' @rdname TwoSampleBinary
get.ts.bin.decision.df = function(a.con = 1, b.con = 1, n.con = 40, x.con = 0:40,
                                  a.trt = 1, b.trt = 1, n.trt = 40, x.trt = 0:40,
                                  Delta.tv = .25, Delta.lrv = .2,
                                  tau.tv = 0.10, tau.lrv = 0.80, tau.ng = 0.65)
{
  
  # create simulation grid
  
  my.grid <- expand.grid(n.con = n.con, n.trt = n.trt,
                         x.con = x.con,
                         x.trt = x.trt,
                         a.con= a.con, a.trt = a.trt,
                         b.con = b.con, b.trt = b.trt,
                         Delta.tv = Delta.tv, Delta.lrv = Delta.lrv,
                         tau.tv = tau.tv,tau.lrv = tau.lrv,tau.ng = tau.ng)
  
  # Define function to be run on each row of grid
  my.function1 <- function(n.con = my.grid$n.con[1], n.trt = my.grid$n.trt[1], x.con = my.grid$x.con[1], x.trt = my.grid$x1[1],
                           a.con = my.grid$a.con[1], a.trt = my.grid$a.trt[1], b.con = my.grid$b.con[1], b.trt = my.grid$b1[1],
                           Delta.tv = my.grid$Delta.tv[1], Delta.lrv = my.grid$Delta.lrv[1],
                           tau.tv = my.grid$tau.tv[1], tau.lrv = my.grid$tau.lrv[1], tau.ng = my.grid$tau.ng[1]){
    
    return(get.ts.bin.decision(x.con = x.con, x.trt = x.trt, n.con = n.con, n.trt = n.trt, a.con = a.con,
                               a.trt = a.trt, b.con = b.con, b.trt = b.trt, 
                               Delta.tv = Delta.tv, Delta.lrv = Delta.lrv,
                               tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng))
  }
  
  # Make call for each row of grid
  listOfDataFrames <- do.call(Map, c(f = my.function1, my.grid))
  # Collapse list of data.frames into a single data.frame
  df <- do.call("rbind", listOfDataFrames)
  my.grid <- cbind(my.grid, df)
  my.grid$result <- factor(my.grid$result, c("Go", "Consider", "No-Go"))
  return(my.grid)
}


#' @rdname TwoSampleBinary
make.ts.bin.spp <- function(a.con = 4, b.con = 36, n.con = 40, x.con = 4,
                            a.trt = 1, b.trt = 1, n.trt = 40, x.trt = 14,
                            Delta.lrv = 0.3, Delta.tv = .4,
                            tau.tv = 0.10, tau.lrv = 0.80, tau.ng = 0.8,
                            nlines.ria = 20, tsize = 4, nlines = 25)
{
  # Run this for fixed value of CON
  results <- get.ts.bin.decision.df(a.con = a.con, b.con = b.con, n.con = n.con, x.con = x.con,
                                    a.trt = a.trt, b.trt = b.trt, n.trt = n.trt, x.trt = 0:n.trt,
                                    Delta.tv = Delta.tv, Delta.lrv = Delta.lrv,
                                    tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng)
  
  # Determine mi/max number of TRT responders for Go/No G0 ----
  result.go <- results %>% dplyr::filter(result== "Go") %>% slice(1)
  result.ng <- results %>% dplyr::filter(result== "No-Go") %>% slice(n())
  
  probs <- 1 - p2beta(relation="DIFF", approach="DIRECT", x=c(Delta.lrv, Delta.tv),
                      a1=a.con+x.con, b1=b.con+n.con-x.con,
                      a2=a.trt+x.trt, b2=b.trt+n.trt-x.trt, n = 1000000)
  
  # comput result -----
  result <- ifelse(probs[1] >= tau.lrv & probs[2] >= tau.tv, "Go",
                   ifelse(probs[1] < tau.ng & probs[2] < tau.tv, "No-Go", "Consider"))
  result.color <- ifelse(result=="Go", "darkgreen", ifelse(result=="No-Go", "red", "black"))
  
  # Subtitle ----
  for.subtitle <- paste0("Control data: ", round(x.con/n.con*100,2), "% (",  x.con, "/", n.con, "). ",
                         "Treatment data: ",round(x.trt/n.trt*100,2), "% (",  x.trt, "/", n.trt, "). ",
                         "Observed difference: ", round(x.trt/n.trt*100,2) - round(x.con/n.con*100,2),
                         "%\nGiven control data, treatment responders needed for Go: ", result.go$x.trt[1], " (",round(result.go$x.trt[1]/n.trt*100,2), "%).",
                         "Needed for No-Go ", result.ng$x.trt[1]," (", round(result.ng$x.trt/n.trt*100,2), "%)")
  
  # Annotation line 1: Decision ----
  for.decision <- paste0("Decision: ", result)
  
  
  # Annotation line 3: P(Delta >= Min TPP)
  annotate.P1 <- ifelse(result=="Go",
                        TeX(paste0("$P(\\Delta\\, \\geq$ Min TPP) = $",  round((probs[1])*100,1), "$% > ", (tau.lrv)*100,"%" )),
                        ifelse(result=="No-Go",
                               TeX(paste0("$P(\\Delta\\, \\geq$ Min TPP) = $",  round((probs[1])*100,1), "% $\\leq$ ", (tau.ng)*100,"%")),
                               TeX(paste0("$P(\\Delta\\, \\geq$ Min TPP) = $",  round((probs[1])*100,1), "%"))))
  # Annotation line 4: P(Delta >= Base TPP)
  
  annotate.P2 <- ifelse(result=="Go",
                        TeX(paste0("$P(\\Delta\\, \\geq$ Base TPP) = $",  round(probs[2]*100,1), "$% >", tau.tv*100,"%")),
                        ifelse(result=="No-Go", TeX(paste0("$P(\\Delta\\, \\geq$ Base TPP$) = $",  round(probs[2]*100, 1), "% $\\leq$", tau.tv*100,"$%")),
                               TeX(paste0("$P(\\Delta\\, \\geq$ Base TPP) = $",  round(probs[2]*100, 1),"$%"))))
  
  
  my.df <- try(data.frame(x = c(seq(-1, -0.001, by = 0.001), seq(0.001, 1, by = 0.001))) %>%
                 mutate(y = (d2beta(relation='DIFF', x=x, a1=a.con+x.con, b1=b.con+n.con-x.con, a2=a.trt+x.trt, b2=b.trt+n.trt-x.trt)))
  )
  
  if(class(my.df) != "try-error"){
    my.df$group <- c(paste0("Control Prior: Beta(", a.con, ", ", b.con, "). Treatment Prior: Beta(", a.trt, ", ", b.trt, ")"))
    
    # Compute probabilities
    # Probs are reported as exceeding
    probs <- 1 - p2beta(relation="DIFF", approach="DIRECT", x=c(Delta.lrv, Delta.tv),
                        a1=a.con+x.con, b1=b.con+n.con-x.con,
                        a2=a.trt+x.trt, b2=b.trt+n.trt-x.trt, n = 100000)
    
    # comput result -----
    result <- ifelse(probs[1] >= tau.lrv & probs[2] >= tau.tv, "Go",
                     ifelse(probs[1] < tau.ng & probs[2] < tau.tv, "No-Go", "Consider"))
    result.color <- ifelse(result=="Go", "darkgreen", ifelse(result=="No-Go", "red", "black"))
    
    # Subtitle ----
    for.subtitle <- paste0("Control data: ", round(x.con/n.con*100,2), "% (",  x.con, "/", n.con, "). ",
                           "Treatment data: ",round(x.trt/n.trt*100,2), "% (",  x.trt, "/", n.trt, "). ",
                           "Observed difference: ", round(x.trt/n.trt*100,2) - round(x.con/n.con*100,2),
                           "%\nGiven control data, treatment responders needed for Go: ", result.go$x.trt[1], " (",round(result.go$x.trt[1]/n.trt*100,2), "%). ",
                           "Needed for No-Go ", result.ng$x.trt[1]," (", round(result.ng$x.trt/n.trt*100,2), "%)")
    
    # Annotation line 1: Decision ----
    for.decision <- paste0("Decision: ", result)
    
    
    # Annotation line 3: P(Delta >= Min TPP)
    annotate.P1 <- ifelse(result=="Go",
                          TeX(paste0("$P(\\Delta\\, \\geq$ Min TPP) = $",  round((probs[1])*100,1), "$% > ", (tau.lrv)*100,"%" )),
                          ifelse(result=="No-Go",
                                 TeX(paste0("$P(\\Delta\\, \\geq$ Min TPP) = $",  round((probs[1])*100,1), "% $\\leq$ ", (tau.ng)*100,"%")),
                                 TeX(paste0("$P(\\Delta\\, \\geq$ Min TPP) = $",  round((probs[1])*100,1), "%"))))
    # Annotation line 4: P(Delta >= Base TPP)
    
    annotate.P2 <- ifelse(result=="Go",
                          TeX(paste0("$P(\\Delta\\, \\geq$ Base TPP) = $",  round(probs[2]*100,1), "$% >", tau.tv*100,"%")),
                          ifelse(result=="No-Go", TeX(paste0("$P(\\Delta \\geq$ Base TPP$) = $",  round(probs[2]*100, 1), "% $\\leq$", tau.tv*100,"$%")),
                                 TeX(paste0("$P(\\Delta\\, \\geq$ Base TPP) = $",  round(probs[2]*100, 1),"$%"))))
    
    # Initialize a ggplot
    dplot <- ggplot() + geom_line(data = my.df, aes(x = x, y=y)) + facet_wrap(~group)
    # Access the ggplot to get goodies to help accomplish shading
    dpb <- ggplot_build(dplot)
    x1.1 <- min(which(dpb$data[[1]]$x >=Delta.lrv))
    x2.1 <- max(which(dpb$data[[1]]$x <=Delta.tv))+1
    x1.2 <- min(which(dpb$data[[1]]$x >=Delta.tv))
    x2.2 <- max(which(dpb$data[[1]]$x <=1))
    
    
    beta.diff.quants <- q2beta(relation="DIFF", a1=a.con+x.con, b1=b.con+n.con-x.con,
                               a2=a.trt+x.trt, b2=b.trt+n.trt-x.trt, alpha=c(1-tau.ng, 1-tau.lrv, 1 - tau.tv), tol = 10^(-5))
    
    if(beta.diff.quants[3]>=Delta.tv){
      # Annotation line 2: Decision interval ----
      for.decision.interval <- paste0("Decision interval: (", round(beta.diff.quants[2]*100,1), "%, ", round(beta.diff.quants[3]*100,1), "%)")
      report.segment <- data.frame(x = beta.diff.quants[2],
                                   xend = beta.diff.quants[3],
                                   y= min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                                   yend = min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                                   group=my.df$group[1])
    } else {  for.decision.interval <- paste0("Decision interval: (", round(beta.diff.quants[1]*100,1), "%, ", round(beta.diff.quants[3]*100,1), "%)")
    report.segment <- data.frame(x = beta.diff.quants[1],
                                 xend = beta.diff.quants[3],
                                 y= min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                                 yend = min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                                 group=my.df$group[1])
    }
    
    if(which(dpb$data[[1]]$y== max(dpb$data[[1]]$y)) > length(dpb$data[[1]]$x)/2){
      annotate.x = min(min(dpb$data[[1]]$x), Delta.lrv,-1)
      annotate.j = 0
    } else {annotate.x = max(max(dpb$data[[1]]$x), Delta.tv, 1)
    annotate.j = 1}
    
    # Introduce shading ----
    main.plot <- dplot +
      geom_area(data = data.frame(x = dpb$data[[1]]$x[1:x1.1],
                                  y = dpb$data[[1]]$y[1:x1.1]),
                aes(x = x, y = y), fill=alpha("red", 0))+
      geom_area(data = data.frame(x = dpb$data[[1]]$x[x1.1:x2.1],
                                  y = dpb$data[[1]]$y[x1.1:x2.1]),
                aes(x = x, y = y), fill=alpha("grey80", 0))+
      geom_area(data = data.frame(x = dpb$data[[1]]$x[x1.2:x2.2],
                                  y = dpb$data[[1]]$y[x1.2:x2.2]),
                aes(x = x, y = y), fill=alpha("green", 0))+
      geom_line(data = my.df, aes(x = x, y=y))+
      scale_x_continuous(limits = c(-1,1),
                         breaks = pretty(x=c(-1,1), n=10),
                         labels = scales::percent)+
      scale_y_continuous(breaks=NULL, labels=NULL)
    
    # Add Annotations -----
    main.plot <- main.plot +
      labs(title = TeX("Posterior distribution for the treatment effect"),
           subtitle = for.subtitle,
           x = TeX("$\\Delta\\,$ = Treatment difference (Treatment - Control)"), y = NULL,
           caption = NULL)+
      annotate("text", label = for.decision,
               x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 0, size = tsize+1, colour = result.color, hjust = annotate.j)+
      annotate("text", label = for.decision.interval,
               x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 1, size = tsize, colour = result.color, hjust = annotate.j)+
      annotate("text", label = annotate.P1, color=result.color,
               x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 2, size = tsize,  hjust = annotate.j)+
      annotate("text", label = annotate.P2, color=result.color,
               x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 3, size = tsize,  hjust = annotate.j)
    
    # Add reference lines and Credible interval
    main.plot <- main.plot +
      geom_segment(data=report.segment, aes(x=x, xend=xend, y=y, yend=yend, group=group), arrow = arrow(ends="both", angle = 90), color=result.color, size=.75)
    
    main.plot <- main.plot +
      geom_vline(xintercept = c(Delta.lrv, Delta.tv), linetype = 2, color = c("blue", "blue")) +
      annotate("text", label = TeX(paste0("Min TPP = $", Delta.lrv*100,"%$")), x = Delta.lrv, y = 0 + max(dpb$data[[1]]$y)/nlines.ria, size = tsize, colour = "black", hjust = 1)+
      annotate("text", label = TeX(paste0("Base TPP = $",  Delta.tv*100,"%$")), x = Delta.tv, y = 0 + max(dpb$data[[1]]$y)/nlines.ria, size = tsize, colour = "black", hjust = 0)
  } else {
    X1 <- rbeta(n=10^6, shape1 = a.con+x.con, shape2 = b.con+n.con-x.con)
    X2 <- rbeta(n=10^6, shape1 = a.trt+x.trt, shape2 = b.trt+n.trt-x.trt)
    my.df <- data.frame(x=X2 - X1, group=c(paste0("Control Prior: Beta(", a.con, ", ", b.con, "). Treatment Prior: Beta(", a.trt, ", ", b.trt, ")")))
    
    dplot <- ggplot() + geom_density(data = my.df, aes(x = x)) + facet_wrap(~group)
    # Access the ggplot to get goodies to help accomplish shading
    dpb <- ggplot_build(dplot)
    x1.1 <- ifelse(length(which(dpb$data[[1]]$x >=Delta.lrv)) == 0, length(dpb$data[[1]]$x >=Delta.lrv),  min(which(dpb$data[[1]]$x >=Delta.lrv)))
    x2.1 <- max(which(dpb$data[[1]]$x <=Delta.tv))
    x1.2 <- ifelse(length(which(dpb$data[[1]]$x >=Delta.tv)) == 0, length(dpb$data[[1]]$x >=Delta.tv),  min(which(dpb$data[[1]]$x >=Delta.tv)))
    x2.2 <- max(which(dpb$data[[1]]$x <=1))
    
    beta.diff.quants <- q2beta(relation="DIFF", a1=a.con+x.con, b1=b.con+n.con-x.con,
                               a2=a.trt+x.trt, b2=b.trt+n.trt-x.trt, alpha=c(1-tau.ng, 1-tau.lrv, 1 - tau.tv), tol = 10^(-5))
    
    if(beta.diff.quants[3]>=Delta.tv){
      # Annotation line 2: Decision interval ----
      for.decision.interval <- paste0("Decision interval: (", round(beta.diff.quants[2]*100,1), "%, ", round(beta.diff.quants[3]*100,1), "%)")
      report.segment <- data.frame(x = beta.diff.quants[2],
                                   xend = beta.diff.quants[3],
                                   y= min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                                   yend = min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                                   group=my.df$group[1])
    } else {  for.decision.interval <- paste0("Decision interval: (", round(beta.diff.quants[1]*100,1), "%, ", round(beta.diff.quants[3]*100,1), "%)")
    report.segment <- data.frame(x = beta.diff.quants[1],
                                 xend = beta.diff.quants[3],
                                 y= min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                                 yend = min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                                 group=my.df$group[1])
    }
    
    if(which(dpb$data[[1]]$y== max(dpb$data[[1]]$y)) > length(dpb$data[[1]]$x)/2){
      annotate.x = min(min(dpb$data[[1]]$x), Delta.lrv,-1)
      annotate.j = 0
    } else {annotate.x = max(max(dpb$data[[1]]$x), Delta.tv, 1)
    annotate.j = 1}
    
    # Introduce shading ----
    main.plot <- dplot +
      geom_area(data = data.frame(x = dpb$data[[1]]$x[1:x1.1],
                                  y = dpb$data[[1]]$y[1:x1.1]),
                aes(x = x, y = y), fill=alpha("red", 0.0))+
      geom_area(data = data.frame(x = dpb$data[[1]]$x[x1.1:x2.1],
                                  y = dpb$data[[1]]$y[x1.1:x2.1]),
                aes(x = x, y = y), fill=alpha("grey80", 0.0))+
      geom_area(data = data.frame(x = dpb$data[[1]]$x[x1.2:x2.2],
                                  y = dpb$data[[1]]$y[x1.2:x2.2]),
                aes(x = x, y = y), fill=alpha("green", 0.0))+
      geom_density(data = my.df, aes(x = x))+
      scale_x_continuous(limits = c(-1,1),
                         breaks = pretty(x=c(-1,1), n=10),
                         labels = scales::percent)+
      scale_y_continuous(breaks=NULL, labels=NULL)
    
    # Add Annotations -----
    main.plot <- main.plot +
      labs(title = TeX("Posterior distribution for the treatment effect"),
           subtitle = for.subtitle,
           x = TeX("$\\Delta\\,$ = Treatment difference (Treatment - Control)"), y = NULL,
           caption = NULL)+
      annotate("text", label = for.decision,
               x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 0, size = tsize+1, colour = result.color, hjust = annotate.j)+
      annotate("text", label = for.decision.interval,
               x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 1, size = tsize, colour = result.color, hjust = annotate.j)+
      annotate("text", label = annotate.P1, color=result.color,
               x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 2, size = tsize,  hjust = annotate.j)+
      annotate("text", label = annotate.P2, color=result.color,
               x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 3, size = tsize,  hjust = annotate.j)
    
    # Add reference lines and Credible interval
    main.plot <- main.plot +
      geom_segment(data=report.segment, aes(x=x, xend=xend, y=y, yend=yend, group=group), arrow = arrow(ends="both", angle = 90), color=result.color, size=.75)
    
    main.plot <- main.plot +
      geom_vline(xintercept = c(Delta.lrv, Delta.tv), linetype = 2, color = c("blue", "blue")) +
      annotate("text", label = TeX(paste0("Min TPP = $", Delta.lrv*100,"%$")), x = Delta.lrv, y = 0 + max(dpb$data[[1]]$y)/nlines.ria, size = tsize, colour = "black", hjust = 1)+
      annotate("text", label = TeX(paste0("Base TPP = $",  Delta.tv*100,"%$")), x = Delta.tv, y = 0 + max(dpb$data[[1]]$y)/nlines.ria, size = tsize, colour = "black", hjust = 0)
  } 
  
  table.plot2 <- ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) > ", tau.lrv*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) > ", tau.tv*100,"%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) $\\leq$ $", tau.ng*100, "$$%\\;\\;\\;\\;$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) $\\leq$ $", tau.tv*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")
  
  # plot_grid(main.plot,
  #            table.plot2, ncol=1, align = "v", rel_heights=c(.78,.22), axis="l")
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
}

#' @rdname TwoSampleBinary
get.ts.bin.trt.oc.df  <- function(a.con = 1, b.con = 1,  dcurve.con = 0.12,
                                  a.trt = 1, b.trt = 1,
                                  TE.OC.N = 80, Aratio=1,
                                  TE.OC.Delta.LB = 0, TE.OC.Delta.UB = 1 - 0.12,
                                  Delta.tv =  .35, Delta.lrv = 0.2,
                                  tau.tv = .01, tau.lrv = 0.8, tau.ng = 0.65){

  # Determine number of control subjects
  n.con = round(TE.OC.N / (1 + Aratio))
  n.trt = TE.OC.N - n.con

  results <- get.ts.bin.decision.df(a.con = a.con, b.con = b.con, n.con = n.con,
                                    a.trt = a.trt, b.trt = b.trt, n.trt = n.trt,
                                    x.trt = 0:n.trt, x.con = 0:n.con,
                                    Delta.tv = Delta.tv, Delta.lrv = Delta.lrv,
                                    tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng)


  my.df <-  bind_rows(apply(X=matrix(seq(TE.OC.Delta.LB, TE.OC.Delta.UB, 0.01)),
                            MARGIN = 1, FUN = function(x) {
                              results %>% mutate(
                                prob = dbinom(x = x.con, size = n.con, prob = dcurve.con) *
                                  dbinom(x = x.trt, size = n.trt, prob = dcurve.con + x)) %>%
                                group_by(result) %>%
                                summarize(prob=sum(prob)) %>% mutate(dcurve.con=dcurve.con, x = x)
                            })) %>%       mutate(a.con = a.con, b.con = b.con,
                                                 a.trt = a.trt, b.trt = b.trt,
                                                 n.con = n.con,  n.trt = n.trt, Aratio = Aratio,
                                                 Delta.tv = Delta.tv, Delta.lrv = Delta.lrv,
                                                 tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng,
                                                 TE.OC.Delta.LB=TE.OC.Delta.LB, TE.OC.Delta.UB=TE.OC.Delta.UB)

  my.df
}


#' @rdname TwoSampleBinary
make.ts.bin.trt.oc1 <- function(my.df=get.ts.bin.trt.oc.df(), nlines=25, tsize=4){

  # Take complement of NoGo for purpose of graphic
  my.df <- my.df %>% mutate(prob = ifelse(result=="No-Go", 1 - prob, prob ) )

  # Pull items for annotation from data.frame
  Delta.tv <- my.df$Delta.tv[1]
  Delta.lrv <- my.df$Delta.lrv[1]
  tau.tv <- my.df$tau.tv[1]
  tau.lrv <- my.df$tau.lrv[1]
  tau.ng <- my.df$tau.ng[1]
  dcurve.con <- my.df$dcurve.con[1]
  Aratio <- my.df$Aratio[1]
  total.n <- my.df$n.con[1] + my.df$n.trt[1]

  # Build plot
  my.plot <- ggplot() +
    geom_line(data = my.df %>%
                dplyr::filter(result=="Go"), aes(x = x, y = prob), color="green")+
    geom_line(data = my.df %>%
                dplyr::filter(result=="No-Go"), aes(x = x, y = prob), color="red")
  dpb <- ggplot_build(my.plot)
  main.plot <- my.plot +
    geom_area(data = dpb$data[[1]], aes(x = x, y = y), fill=alpha("green", .5))+
    geom_ribbon(data = dpb$data[[1]], aes(x = x, ymin = y, ymax=1),
                fill=alpha("grey", .5))+
    geom_ribbon(data = dpb$data[[2]], aes(x = x, ymin = y, ymax = 1),
                fill=alpha("red", .5))+
    geom_line(data = my.df %>% dplyr::filter(result=="Go"), aes(x = x, y = prob),
              color="green", size=.75)+
    geom_line(data = my.df %>% dplyr::filter(result=="No-Go"), aes(x = x, y = prob),
              color="red", size=.75)+
    geom_vline(xintercept = c(Delta.tv, Delta.lrv), color="blue", linetype=2)+
    annotate("text", label = TeX(paste0("Min TPP = $", Delta.lrv*100, "%$")),
             x = Delta.lrv, y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize,
             colour = "black", hjust = 1)+
    annotate("text", label = TeX(paste("Base TPP = $", Delta.tv*100, "%$")),
             x = Delta.tv, y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize,
             colour = "black", hjust = 0)+
    labs(x = "Underlying treatment effect",
         y="Probability",
         title = "Operating characteristics as a function of treatment effect",
         subtitle = paste0("Total randomized: ", total.n, " with ", my.df$n.con[1],
                           " to Control and ", my.df$n.trt[1], " to Treatment"))+
    scale_x_continuous(expand = c(0,0),
                       breaks=seq(my.df$TE.OC.Delta.LB[1],
                                  my.df$TE.OC.Delta.UB[1], length.out=5),
                       limits=c(head(dpb$data[[1]]$x,1),tail(dpb$data[[1]]$x,1)),
                       label=scales::percent) +
    scale_y_continuous(expand = c(0,0), breaks=seq(0,1,.2), minor_breaks=seq(0,1,.1),
                       label=scales::percent)+
    theme(panel.spacing.x = unit(6, "mm"), axis.text.x = element_text(angle=45,
                                                                      hjust=1,vjust=1))

  table.plot2 <- ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) > ", tau.lrv*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) > ", tau.tv*100,"%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) $\\leq$ $", tau.ng*100, "$$%\\;\\;\\;\\;$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) $\\leq$ $", tau.tv*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))

}

#' @rdname TwoSampleBinary
make.ts.bin.trt.oc2 <- function(my.df=get.ts.bin.trt.oc.df(), tsize=4, nlines=25){
  # Pull items for annotation from data.frame
  Delta.tv <- my.df$Delta.tv[1]
  Delta.lrv <- my.df$Delta.lrv[1]
  tau.tv <- my.df$tau.tv[1]
  tau.lrv <- my.df$tau.lrv[1]
  tau.ng <- my.df$tau.ng[1]
  dcurve.con <- my.df$dcurve.con[1]
  Aratio <- my.df$Aratio[1]
  total.n <- my.df$n.con[1] + my.df$n.trt[1]

  my.plot <- ggplot(data = my.df, aes(x = x, y = prob, color = result)) +
    geom_line(size=.75)
  dpb <- ggplot_build(my.plot)
  main.plot <-  my.plot +
    labs(x = "Underlying treatment effect",
         y="Probability", color="Decision",
         title = "Operating characteristics as a function of treatment effect",
         subtitle = paste0("Total randomized: ", total.n, " with ", my.df$n.con[1],
                           " to Control and ", my.df$n.trt[1], " to Treatment"))+
    scale_color_manual(values = c("Go" = "green", "No-Go" = "red",
                                  "Consider" = "darkgrey"))+
    geom_vline(xintercept = c(Delta.tv, Delta.lrv), linetype=2, color="blue")+
    annotate("text", label = TeX(paste0("Min TPP = $", Delta.lrv*100, "%$")),
             x = Delta.lrv, y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize,
             colour = "black", hjust = 1)+
    annotate("text", label = TeX(paste("Base TPP = $", Delta.tv*100, "%$")),
             x = Delta.tv, y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize,
             colour = "black", hjust = 0)+
    theme(legend.position = "bottom")+
    scale_x_continuous(expand = c(0,0),
                       breaks=seq(my.df$TE.OC.Delta.LB[1],
                                  my.df$TE.OC.Delta.UB[1], length.out=5),
                       limits=c(head(dpb$data[[1]]$x,1),tail(dpb$data[[1]]$x,1)),
                       label=scales::percent) +
    scale_y_continuous(expand = c(0,0), breaks=seq(0,1,.2), minor_breaks=seq(0,1,.1),
                       label=scales::percent)+
    theme(panel.spacing.x = unit(6, "mm"),
          axis.text.x = element_text(angle=45, hjust=1,vjust=1))

  table.plot2 <- ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) > ", tau.lrv*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) > ", tau.tv*100,"%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) $\\leq$ $", tau.ng*100, "$$%\\;\\;\\;\\;$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) $\\leq$ $", tau.tv*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")

  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
}


#' @rdname TwoSampleBinary
get.ts.bin.ssize.oc.df <- function(a.con = 1, b.con = 1,
                                   a.trt = 1, b.trt = 1,
                                   dcurve.con = .12,
                                   TE.OC.N = 80, Aratio=2,
                                   SS.OC.N.LB = 40, SS.OC.N.UB = 160,
                                   Delta.lrv = .2, Delta.tv=.25, SS.OC.Delta = .25,
                                   tau.tv = 0.10, tau.lrv = .8, tau.ng = .65,
                                   Delta.user = .30, nlines = 15, tsize = 4, npoints=3){

  # Run total sample size over all integers between bounds
  specs <- data.frame(TE.OC.N = seq(SS.OC.N.LB, SS.OC.N.UB, 1)) %>%
    mutate(n.con = round(TE.OC.N/(1+Aratio)),
           n.trt = TE.OC.N - n.con,
           dcurve.con = dcurve.con,
           p.con1 = dcurve.con + Delta.lrv,
           p.con2 = dcurve.con + Delta.tv,
           p.con3 = dcurve.con + SS.OC.Delta,
           a.con = a.con, b.con = b.con,
           a.trt = a.trt, b.trt = b.trt,
           Delta.lrv = Delta.lrv, Delta.tv=Delta.tv,
           tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng, Aratio=Aratio) %>%
    dplyr::filter(n.trt == Aratio*n.con)
  # This function sets the observed placebo rate based on rounding the user's dcurve.con*n.con
  runit <- function(x){
    get.ts.bin.decision.df(a.con = specs$a.con[x], b.con = specs$b.con[x],
                           a.trt = specs$a.trt[x], b.trt = specs$b.trt[x],
                           n.con = specs$n.con[x], x.con = 0:specs$n.con[x],
                           n.trt = specs$n.trt[x], x.trt = 0:specs$n.trt[x],
                           Delta.tv = specs$Delta.tv[x], Delta.lrv = specs$Delta.lrv[x],
                           tau.tv = specs$tau.tv[x], tau.lrv = specs$tau.lrv[x],
                           tau.ng = specs$tau.ng[x]) %>%
      mutate(result = factor(result, c("Go", "Consider", "No-Go")),
             total.n=n.con+n.trt,
             prob.null = dbinom(x=x.trt, size = n.trt, prob = specs$dcurve.con[x])*
               dbinom(x=x.con, size = n.con, prob = specs$dcurve.con[x]),
             prob.lrv = dbinom(x=x.trt, size = n.trt, prob = specs$p.con1[x])*
               dbinom(x=x.con, size = n.con, prob = specs$dcurve.con[x]),
             prob.tv = dbinom(x=x.trt, size = n.trt, prob = specs$p.con2[x])*
               dbinom(x=x.con, size = n.con, prob = specs$dcurve.con[x]),
             prob.user = dbinom(x=x.trt, size = n.trt, prob = specs$p.con3[x])*
               dbinom(x=x.con, size = n.con, prob = specs$dcurve.con[x]),
             Aratio = specs$Aratio[x]) %>%
      group_by(result) %>%
      summarize(prob.null = sum(prob.null),
                prob.lrv = sum(prob.lrv),
                prob.tv = sum(prob.tv),
                prob.user = sum(prob.user)) %>%
      mutate(a.con = specs$a.con[x], b.con = specs$b.con[x],
             a.trt = specs$a.trt[x], b.trt = specs$b.trt[x],
             n.con = specs$n.con[x], n.trt = specs$n.trt[x],
             n.total=n.trt + n.con, Delta.tv = specs$Delta.tv[x],
             Delta.lrv = specs$Delta.lrv[x], tau.tv = specs$tau.tv[x],
             tau.lrv = specs$tau.lrv[x], tau.ng = specs$tau.ng[x],
             Aratio=specs$Aratio[x], dcurve.con=specs$dcurve.con[x]) }

  big.results <- bind_rows(apply(X = matrix(seq(1,nrow(specs), length.out = npoints)),
                                 MARGIN=1, FUN= function(x) runit(x) )) %>%
    gather(key=key, value=value, prob.null, prob.lrv, prob.tv, prob.user, factor_key = T)

  levels(big.results$key) <-c(TeX("$\\Delta\\,$ = NULL = 0%"),
                              TeX(paste("$\\Delta\\,$ = Min TPP = $", Delta.lrv*100, "%$ ")),
                              TeX(paste("$\\Delta\\,$ = Base TPP = $", Delta.tv*100, "%$ ")),
                              TeX(paste("$\\Delta\\,$ = User defined = $", SS.OC.Delta*100,
                                        "%$")))
  big.results$key <- factor(big.results$key,
                            levels(big.results$key)[order(c(0, Delta.lrv, Delta.tv,
                                                            SS.OC.Delta))])

  return(big.results)
}


#' @rdname TwoSampleBinary
make.ts.bin.ssize.oc <- function(for.plot = get.ts.bin.ssize.oc.df(), tsize=4,
                                 nlines=25, npoints=5){

  for.plot <- for.plot %>% mutate(value = ifelse(result=="No-Go", 1 - value, value))
  main.plot <-  for.plot %>%
    dplyr::filter(result!="Consider") %>%
    ggplot(aes(x=n.total, y=value, color=result)) +
    geom_line(alpha=.25, size=.75) +
    facet_wrap(~key, labeller="label_parsed")+
    scale_color_manual(values=c("green", "red"))+
    geom_ribbon(data=for.plot %>% dplyr::filter(result=="Go"),
                aes(x=n.total, ymin=0, ymax=value),
                fill=alpha("green",.5), color=alpha("green",.25))+
    geom_ribbon(data=for.plot %>% dplyr::filter(result=="Go"),
                aes(x=n.total, ymin=value, ymax=1),
                fill=alpha("grey",.5), color=alpha("grey",0))+
    geom_ribbon(data=for.plot %>% dplyr::filter(result=="No-Go"),
                aes(x=n.total, ymin=value, ymax=1), fill=alpha("red",.5),
                color=alpha("red", .25))+
    guides(color=F, fill=F)+
    labs(x="Total sample size", y="Decision Probabilities",
         title=paste0("Operating Characteristics as a function of sample size when control response rate is ",
                      round(for.plot$dcurve.con[1]*100,2), "%"))+
    scale_x_continuous(expand=c(0,0))+
    scale_y_continuous(expand=c(0,0), breaks=seq(0,1,.2), minor_breaks=seq(0,1,.1),
                       label=scales::percent)+
    theme(panel.spacing.x = unit(6, "mm"), axis.text.x =
            element_text(angle=45, hjust=1,vjust=1))

  head(for.plot)
  tau.lrv <- for.plot$tau.lrv[1]
  tau.tv <- for.plot$tau.tv[1]
  tau.ng <- for.plot$tau.lrv[1]
  table.plot2 <- ggplot()+
    
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) > ", tau.lrv*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) > ", tau.tv*100,"%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) $\\leq$ $", tau.ng*100, "$$%\\;\\;\\;\\;$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) $\\leq$ $", tau.tv*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
}

# OneSampleNormalGamma ------
#' @name OneSampleNormalGamma
#' @title Functions to Support the One Sample Continuous Scenario
#' @param mu.0.t prior mean
#' @param n.0.t prior effective sample size
#' @param alpha.0.t prior alpha parameter
#' @param beta.0.t prior beta parameter
#' @param xbar.t observed sample mean
#' @param s.t observed sample standard deviation
#' @param n.t sample size
#' @param Delta.lrv TPP Lower Reference Value aka Min TPP
#' @param Delta.tv TPP Target Value aka Base TPP
#' @param tau.tv threshold associated with Base TPP
#' @param tau.lrv threshold associated with Min TPP
#' @param tau.ng threshold associated with No-Go
#' @param seed random seed
#' @param nlines Control for text spacing
#' @param tsize Control for text size
#' @param nlines.ria Control for text spacing
#' @param from.here Lower bound for treatment effect
#' @param to.here Upper bound for treatment effect
#' @param my.df data.frame returned by get.ss.ng.trt.oc.df
#' @param n_LB_OC Lower bound for sample size
#' @param n_UB_OC Upper bound for sample size
#' @param npoints Number of points to use in plot
#' @examples
#' \dontrun{
#' make.ss.ng.ppp()
#' make.ss.ng.spp()
#' get.ss.ng.trt.oc.df()
#' make.ss.ng.trt.oc1()
#' make.ss.ng.trt.oc2()
#' get.ss.ng.ssize.oc.df()
#' make.ss.ng.ssize.oc()
#' }
#' @description
#'
#'   make.ss.ng.ppp: Make One Sample Normal-Gamma Prior/Posterior Plot. Returns a ggplot object.
#'
#'   make.ss.ng.spp: Make One Sample Normal-Gamma Shaded Posterior Plot.  Returns a graphic built using grid.arrange.
#'
#'   get.ss.ng.trt.oc.df: Get One Sample Normal-Gamma Treatment Effect OC. Returns a data.frame.
#'
#'   make.ss.ng.trt.oc1: Make One Sample Normal-Gamma Treatment Effect. Returns a graphic built using grid.arrange.
#'
#'   make.ss.ng.trt.oc2: Make One Sample Normal-Gamma Treatment Effect.  Returns a graphic built using grid.arrange.
#'
#'   get.ss.ng.ssize.oc.df:  Get One Sample Normal-Gamma sample size OC data.frame. Returns a data.frame.
#'
#'   make.ss.ng.ssize.oc: Make One Sample Normal-Gamma Sample size OC plot
make.ss.ng.ppp <- function(mu.0.t = 0, n.0.t = 10, alpha.0.t = .25, beta.0.t = 1,
                           xbar.t = 1.75, s.t = 2, n.t = 50) {
  pp <-   get.NG.post(mu.0 = mu.0.t, n.0 = n.0.t, alpha.0 = alpha.0.t, beta.0 = beta.0.t,
                      xbar = xbar.t, s = s.t, n = n.t, group = "Treatment")
  tau.limits <- c(min(pmax(.1, qgamma(p = .005,shape = alpha.0.t, rate = beta.0.t)),
                      pmax(.1, qgamma(p = .005,shape = pp$alpha.n, rate = pp$beta.n))),
                  max(qgamma(p = .995,shape = alpha.0.t, rate = beta.0.t),
                      qgamma(p = .995,shape = pp$alpha.n, rate = pp$beta.n)))

  mu.limits <- c(min(qnorm(p = .005, mean = mu.0.t,
                           sd = (alpha.0.t / beta.0.t) ^ (-1)/sqrt(n.0.t)),
                     qnorm(p = .005, mean = pp$mu.n,
                           sd = (pp$alpha.n / pp$beta.n) ^ (-1)/sqrt(pp$n.n))),
                 max(qnorm(p = .995, mean = mu.0.t,
                           sd = (alpha.0.t / beta.0.t) ^ (-1)/sqrt(n.0.t)),
                     qnorm(p = .995, mean = pp$mu.n,
                           sd = (pp$alpha.n / pp$beta.n) ^ (-1)/sqrt(pp$n.n))))

  prior <- expand.grid(
    tau = seq(tau.limits[1],tau.limits[2], length.out = 100),
    mu = seq(mu.limits[1],mu.limits[2], length.out = 100)) %>%
    mutate(n0 = n.0.t, a0 = alpha.0.t, b0 = beta.0.t, category ="Treatment",
           dens = dnorgam(mu = mu, tau = tau, mu0 = mu.0.t, n0 = n.0.t,
                          a0 = alpha.0.t, b0 = beta.0.t),
           color = as.numeric(cut((dens),100)),
           group = "Prior",
           density = "Treatment Prior")

  post <- expand.grid(
    tau = seq(tau.limits[1],tau.limits[2], length.out = 100),
    mu = seq(mu.limits[1],mu.limits[2], length.out = 100)) %>%
    mutate(n0 = pp$n.n, a0 = pp$alpha.n, b0 = pp$beta.n, category ="Treatment",
           dens = dnorgam(mu = mu, tau = tau, mu0 = pp$mu.n, n0 = pp$n.n,
                          a0 = pp$alpha.n, b0 = pp$beta.n),
           color = as.numeric(cut((dens),100)),
           group = "Posterior",
           density="Treatment Posterior")

  my.df <- rbind(prior,post) %>% mutate(density=factor(density, c("Treatment Prior",
                                                                  "Treatment Posterior")))
  levels(my.df$density) <- c(
    paste0("Treatment Prior: NG(", round(mu.0.t, 2), ", ", round(n.0.t, 2), ", ",
           round(alpha.0.t, 2), ", ", round(beta.0.t, 2), ")"),
    paste0("Treatment Posterior NG(", round(pp$mu.n, 2),", ", round(pp$n.n, 2), ", ",
           round(pp$alpha.n, 2), ", ", round(pp$beta.n, 2), ")" )  )

  my.df$color2 <- as.numeric(cut((my.df$dens),150))

  my.colors <- colorRampPalette(c("grey80", "red",  "yellow"))(150)
  joint.pdf <- ggplot(data= my.df, aes(x = mu, y = tau, fill = factor(color2)))+
    geom_tile() +facet_wrap(~density)+guides(fill=F)+
    scale_x_continuous(expand = c(0,0))+
    scale_y_continuous(expand = c(0,0),
                       breaks = seq(0,max(my.df$tau),.25),
                       labels= round((seq(0,max(my.df$tau),.25)) ^ -.5,2),
                       sec.axis = sec_axis(~., name = TeX("precision parameter, $\\tau$")))+
    labs("Mean treatment response",
         y = TeX("standard deviation, $\\sigma$"),
         title="Prior and posterior joint distributions for mean and precision parameters",
         subtitle=TeX(paste0("Treatment data (n, $\\bar{x}$, s): (",
                             round(n.t,2), ", ", round(xbar.t,2), ", ", round(s.t,2),")")))+
    scale_fill_manual(values= my.colors)

  marginal.pdf <- rbind(
    gcurve(dt_ls(x,df = 2 * alpha.0.t, mu = mu.0.t,
                 sigma = beta.0.t/(alpha.0.t * n.0.t)) ,
           from = mu.limits[1],
           to = mu.limits[2], n = 1001, category = "Treatment") %>%
      mutate(group="Prior", density="Treatment Prior"),
    gcurve(dt_ls(x,df = 2 * pp$alpha.n, mu = pp$mu.n,
                 sigma = pp$beta.n/(pp$alpha.n * pp$n.n )),
           from = mu.limits[1],
           to = mu.limits[2], n = 1001, category = "Treatment") %>%
      mutate(group="Posterior", density="Treatment Posterior")) %>%
    mutate(density = factor(density, c("Treatment Prior", "Treatment Posterior")))

  levels(marginal.pdf$density) <- c(
    paste0("Treatment Prior: t(", 2 * round(alpha.0.t, 2), ", ", round(mu.0.t, 2),
           ", ", round(beta.0.t/(alpha.0.t * n.0.t), 2), ")"),
    paste0("Treatment Posterior: t(", 2 * round(pp$alpha.n), ", ", round(pp$mu.n, 2),
           ", ", round(pp$beta.n/(pp$alpha.n * pp$n.n ), 2), ")")  )

  marginal.plot <-  ggplot(data = marginal.pdf,
                           aes(x = x,y = y, color = category))+
    geom_line(size=.75)+ facet_wrap(~density)+guides(color=F)+
    theme(legend.position = "bottom")+
    labs(title="Marginal distribution for the mean",color = NULL,
         x = TeX("Mean treatment response"), y="")+
    scale_y_continuous(breaks=NULL)

  grid.arrange(joint.pdf, marginal.plot, ncol = 1)

}
#' @rdname OneSampleNormalGamma
make.ss.ng.spp <- function(mu.0.t = 0, alpha.0.t=.25, beta.0.t = 1, n.0.t = 10,
                           xbar.t = 1.97, s.t = 2, n.t = 20,
                           Delta.lrv = 1.25, Delta.tv = 1.75,
                           tau.tv=.1, tau.lrv=.8, tau.ng=.65,
                           tsize=4, nlines=25, nlines.ria=20){

  # Run this for fixed values of
  results <- get.ss.ng.df(mu.0.t = mu.0.t, n.0.t = n.0.t, alpha.0.t = alpha.0.t,
                          beta.0.t = beta.0.t,
                          xbar.t = seq(Delta.lrv*.75, Delta.tv*1.5, length.out=1000),
                          s.t = s.t, n.t = n.t, Delta.tv = Delta.tv, Delta.lrv = Delta.lrv,
                          tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng)

  # Determine min number of TRT responders for Go
  result.go <- results %>% dplyr::filter(result== "Go") %>% slice(1)
  # Determine max number of TRT responders for No-Go
  result.ng <- results %>% dplyr::filter(result== "No-Go") %>% slice(n())

  # Get post parameters
  pp.t <- get.NG.post(mu.0 = mu.0.t, n.0 = n.0.t, alpha.0 = alpha.0.t,
                      beta.0 = beta.0.t,
                      xbar = xbar.t, s = s.t, n = n.t, group = "Treatment")
  my.df <- gcurve(expr = dt_ls(x = x,df = 2 * pp.t$alpha.n, mu = pp.t$mu.n,
                               sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n)),
                  from = pp.t$mu.n - 5 * s.t / sqrt(pp.t$n.n), to = pp.t$mu.n +
                    5 * s.t / sqrt(pp.t$n.n),n = 1001)

  # Compute result
  P.R1 = 1 - pt_ls(x = Delta.lrv, df = pp.t$tdf.n, mu = pp.t$mu.n, sigma = pp.t$sigma.n)
  P.R3 = 1 - pt_ls(Delta.tv, df = pp.t$tdf.n, mu = pp.t$mu.n, sigma = pp.t$sigma.n)
  result = ifelse(P.R1 > tau.lrv & P.R3 > tau.tv, "Go",
                  ifelse(P.R1 < tau.ng  & P.R3 < tau.tv, "No-Go", "Consider"))
  result.color <- ifelse(result=="Go", "darkgreen", ifelse(result=="No-Go", "red", "black"))



  for.subtitle <- TeX(paste0("Treatment data ", "(n, $\\bar{x}$, s): (",
                             round(n.t,2), ", $", round(xbar.t,2), "$, ",
                             round(s.t,2),"). Sample mean needed for Go: $",
                             round(result.go$xbar,2), "$. Needed for No-Go: $",
                             round(result.ng$xbar,2), "$."))

  my.df$group = paste0("Prior: NG(", round(mu.0.t,2), ", ", round(n.0.t,2),
                       ", ", round(alpha.0.t,2), ", ", round(beta.0.t,2), "); ")

  # Annotation line 1: Decision ----
  for.decision <- paste0("Decision: ", result)

  # Annotation line 2: Decision interval ----
  if(result!="No-Go"){
    for.decision.interval <- paste0(
      "Decision interval: (",
      round(qt_ls(p = 1 - tau.lrv,
                  df = 2 * pp.t$alpha.n, mu = pp.t$mu.n,
                  sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n)),3), ", ",
      round(qt_ls(p = 1 - tau.tv,
                  df = 2 * pp.t$alpha.n, mu = pp.t$mu.n,
                  sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n)),3), ")")
  } else {
    for.decision.interval <- paste0(
      "Decision interval: (",
      round(qt_ls(p = 1 - tau.ng,
                  df = 2 * pp.t$alpha.n, mu = pp.t$mu.n,
                  sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n)),3)*100, "%, ",
      round(qt_ls(p = 1 - tau.tv, df = 2 * pp.t$alpha.n, mu = pp.t$mu.n,
                  sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n)),3)*100, "%)")
  }

  # Annotation line 3: P(Delta >= Min TPP)
  annotate.P1 <- ifelse(
    result=="Go",
    TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) = $",
               round(1 - pt_ls(x = Delta.lrv,
                               df = 2 * pp.t$alpha.n,
                               mu = pp.t$mu.n,
                               sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n)),
                     4)*100, "%$ >", tau.lrv*100, "%")),
    ifelse(result=="No-Go",
           TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) = $",
                      round(1 - pt_ls(x = Delta.lrv, df = 2 * pp.t$alpha.n,
                                      mu = pp.t$mu.n,
                                      sigma = pp.t$beta.n/(pp.t$alpha.n* pp.t$n.n)),4)*100,
                      "%$ < ", tau.ng*100, "%")),
           TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) = $",
                      round(1 - pt_ls(x = Delta.lrv, df = 2 * pp.t$alpha.n,
                                      mu = pp.t$mu.n,
                                      sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n)),
                            4)*100, "%$"))))

  annotate.P2 <- ifelse(
    result=="Go",
    TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) = $",
               round(1 - pt_ls(x = Delta.tv,
                               df = 2 * pp.t$alpha.n,
                               mu = pp.t$mu.n,
                               sigma = pp.t$beta.n/(pp.t$alpha.n *pp.t$n.n)),4)*100,
               "$ > ", tau.tv*100, "%")),
    ifelse(
      result=="No-Go",
      TeX(paste0("P($\\Delta\\,$ > Base TPP$) = $",
                 round(1 - pt_ls(x = Delta.tv,
                                 df = 2 * pp.t$alpha.n,
                                 mu = pp.t$mu.n,
                                 sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n)),4)*100,
                 "$ < ", tau.tv*100,"%")),
      TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) = $",
                 round(1 - pt_ls(x = Delta.tv,
                                 df = 2 * pp.t$alpha.n,
                                 mu = pp.t$mu.n,
                                 sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n)),
                       4)*100, "%"))))

  dplot <- ggplot() + geom_line(data = my.df, aes(x = x, y = y)) +
    facet_wrap(~group)
  # Access the ggplot to get goodies to help accomplish shading
  dpb <- ggplot_build(dplot)
  x1.1 <- min(which(dpb$data[[1]]$x >=Delta.lrv))
  x2.1 <- max(which(dpb$data[[1]]$x <=Delta.tv))+1
  x1.2 <- pmin(min(which(dpb$data[[1]]$x >=Delta.tv)),
               max(which(dpb$data[[1]]$x <=Inf)))
  x2.2 <- max(which(dpb$data[[1]]$x <=Inf))

  # Custom graphic parameters
  x.limits <- c(min(min(dpb$data[[1]]$x), Delta.lrv),
                max(max(dpb$data[[1]]$x), Delta.tv))
  ticks <- pretty(x = c(x.limits[1], x.limits[2]), n = 15)
  # When the mound is on the right... we want to display annotation on the left
  if(which(dpb$data[[1]]$y== max(dpb$data[[1]]$y)) > length(dpb$data[[1]]$x)/2){
    annotate.x = min(min(dpb$data[[1]]$x), Delta.lrv)
    annotate.j = 0
  } else {annotate.x = max(max(dpb$data[[1]]$x), Delta.tv)
  annotate.j = 1}


  go.segment <- data.frame(x = qt_ls(p = 1 - tau.lrv,
                                     df = 2 * pp.t$alpha.n, mu = pp.t$mu.n,
                                     sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n)),
                           xend = qt_ls(p = 1 - tau.tv,
                                        df = 2 * pp.t$alpha.n, mu = pp.t$mu.n,
                                        sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n)),
                           y= min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                           yend = min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                           group=my.df$group[1])
  nogo.segment <- data.frame(x = qt_ls(p = 1 - tau.ng,
                                       df = 2 * pp.t$alpha.n,
                                       mu = pp.t$mu.n,
                                       sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n)),
                             xend = qt_ls(p = 1 - tau.lrv,
                                          df = 2 * pp.t$alpha.n, mu = pp.t$mu.n,
                                          sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n)),
                             y= min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                             yend = min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                             group=my.df$group[1])
  # Introduce shading
  main.plot <- dplot +
    geom_area(data = data.frame(x = dpb$data[[1]]$x[1:x1.1],
                                y = dpb$data[[1]]$y[1:x1.1]),
              aes(x = x, y = y), fill=alpha("red",0))+
    geom_area(data = data.frame(x = dpb$data[[1]]$x[x1.1:x2.1],
                                y = dpb$data[[1]]$y[x1.1:x2.1]),
              aes(x = x, y = y), fill=alpha("grey80",0))+
    geom_area(data = data.frame(x = dpb$data[[1]]$x[x1.2:x2.2],
                                y = dpb$data[[1]]$y[x1.2:x2.2]),
              aes(x = x, y = y), fill=alpha("green",0))+
    geom_vline(xintercept = c(Delta.lrv, Delta.tv), linetype = 2, color = c("blue", "blue"))+
    scale_y_continuous(breaks=NULL) +
    scale_x_continuous(limits = x.limits, breaks=pretty(x = x.limits, n=10))

  # Add Annotations -----

  main.plot <- main.plot +
    labs(title = "Posterior distribution for the treatment effect",
         x = "Mean treatment response", y = NULL,
         subtitle= for.subtitle,
         caption = NULL)+
    annotate("text", label = for.decision,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 0,
             size = tsize+1, colour = result.color, hjust = annotate.j)+
    annotate("text", label = for.decision.interval,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 1,
             size = tsize, colour = result.color, hjust = annotate.j)+
    annotate("text", label = annotate.P1, color=result.color,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 2,
             size = tsize,  hjust = annotate.j)+
    annotate("text", label = annotate.P2, color=result.color,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 3,
             size = tsize,  hjust = annotate.j)

  # Add reference lines and Credible interval
  if(result.color == "red"){
    main.plot <- main.plot +
      geom_segment(data=nogo.segment, aes(x=x, xend=xend, y=y, yend=yend, group=group),
                   arrow = arrow(ends="both", angle = 90), color=result.color, size=.75)
  } else {
    main.plot <- main.plot +
      geom_segment(data=go.segment, aes(x=x, xend=xend, y=y, yend=yend, group=group),
                   arrow = arrow(ends="both", angle = 90), color=result.color, size=.75)
  }
  main.plot <- main.plot +
    geom_vline(xintercept = c(Delta.lrv, Delta.tv), linetype = 2, color = c("blue", "blue")) +
    annotate("text", label = paste0("Min TPP = ", Delta.lrv), x = Delta.lrv,
             y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize, colour = "black", hjust = 1)+
    annotate("text", label = paste0("Base TPP = ",  Delta.tv), x = Delta.tv,
             y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize, colour = "black", hjust = 0)

  table.plot2 <- ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) > ", tau.lrv*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) > ", tau.tv*100,"%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) $\\leq$ $", tau.ng*100, "$$%\\;\\;\\;\\;$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) $\\leq$ $", tau.tv*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")

  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
}
#' @rdname OneSampleNormalGamma
get.ss.ng.df <- function(mu.0.t = 0, n.0.t = 10, alpha.0.t = .25, beta.0.t = 1,
                         xbar.t = seq(-1,5,.1), s.t = seq(1,6,.1), n.t = 50,
                         Delta.tv = 1.75, Delta.lrv = 1.5, tau.tv = 0.10,
                         tau.lrv = 0.80, tau.ng = 0.65)
{
  # Create a simulation grid
  my.grid <- expand.grid(
    mu.0 = mu.0.t, n.0 = n.0.t, alpha.0 = alpha.0.t, beta.0 = beta.0.t,
    xbar = xbar.t,  s = s.t, n = n.t,
    Delta.tv = Delta.tv, Delta.lrv = Delta.lrv,
    tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng)

  my.results <- bind_rows(apply(X = matrix(1:nrow(my.grid)), MARGIN = 1,
                                FUN = function(x){
                                  get.NG.post(mu.0 = my.grid$mu.0[x], n.0 = my.grid$n.0[x],
                                              alpha.0 = my.grid$alpha.0[x],
                                              beta.0 = my.grid$beta.0[x],
                                              xbar = my.grid$xbar[x], s = my.grid$s[x], n = my.grid$n[x],
                                              group = "Treatment")})) %>%
    mutate(Delta.tv =Delta.tv, Delta.lrv=Delta.lrv, tau.tv=tau.tv,
           tau.lrv=tau.lrv, tau.ng=tau.ng) %>%
    mutate(
      P.R1 = 1 - pt_ls(x = Delta.lrv, df = tdf.n, mu = mu.n, sigma = sigma.n),
      P.R3 = 1 - pt_ls(x = Delta.tv, df = tdf.n, mu = mu.n, sigma = sigma.n),
      result = ifelse(P.R1 > tau.lrv & P.R3 > tau.tv, "Go",
                      ifelse(P.R1 <= tau.ng  & P.R3 <= tau.tv, "No-Go", "Consider"))) %>%
    mutate(result = factor(result, c("Go", "Consider", "No-Go")))

  return(my.results)
}

#' @rdname OneSampleNormalGamma

get.ss.ng.trt.oc.df <- function(mu.0.t = 0, n.0.t = 10, alpha.0.t = 0.25, beta.0.t = 1,
                                s.t = 2, n.t = 40, from.here = 0, to.here =4,
                                Delta.tv = 1.75, Delta.lrv = 1,
                                tau.tv = 0.1, tau.lrv = 0.8, tau.ng = 0.65)  {

  # Create a grid to pass to decision.ss.continuous
  # This grid really just runs over a sequence of values for underlying treatment effect

  results <- get.ss.ng.df(mu.0.t = mu.0.t, n.0.t = n.0.t, alpha.0.t = alpha.0.t,
                          beta.0.t = beta.0.t,
                          xbar.t = seq(from.here, to.here, length.out = 1000),
                          s.t = s.t, n.t = n.t,
                          Delta.tv = Delta.tv, Delta.lrv = Delta.lrv, tau.tv = tau.tv,
                          tau.lrv = tau.lrv, tau.ng = tau.ng)

  result.go <- results %>% dplyr::filter(result== "Go") %>% slice(1)
  # Determine max number of TRT responders for No-Go
  result.ng <- results %>% dplyr::filter(result== "No-Go") %>% slice(n())


  my.df <- data.frame(xbar=seq(from.here, to.here, length.out = 1000)) %>%
    mutate(Go= 1 - pnorm(q=result.go$xbar, mean=xbar, sd=s.t/sqrt(n.t)),
           NoGo= pnorm(q = result.ng$xbar, mean=xbar, sd=s.t/sqrt(n.t))) %>%
    mutate(mu.0.t = mu.0.t, n.0.t = n.0.t, alpha.0.t = alpha.0.t, beta.0.t = beta.0.t,
           s.t = s.t, n.t = n.t,
           Delta.tv = Delta.tv, Delta.lrv = Delta.lrv,
           tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng)

  my.df
}

#' @rdname OneSampleNormalGamma
make.ss.ng.trt.oc1 <- function(my.df = get.ss.ng.trt.oc.df(), nlines=25, tsize=4){
  
  
  my.plot <- my.df  %>%
    ggplot() + geom_line(aes(x = xbar, y = Go), color="green")+
    geom_line(aes(x = xbar, y = 1-NoGo), color="red")
  dpb <- ggplot_build(my.plot)
  
  main.plot <- my.plot+
    geom_ribbon(data = my.df, aes(x = xbar, ymin=0, ymax=Go), fill=alpha("lightgreen", .5))+
    geom_ribbon(data = my.df, aes(x = xbar, ymin=Go, ymax=1-NoGo), fill=alpha("grey", .5))+
    geom_ribbon(data = my.df, aes(x = xbar, ymin=1-NoGo, ymax=1), fill=alpha("red", .5))+
    geom_line(data = my.df, aes(x = xbar, y = Go), color="green", size=.75)+
    geom_line(data = my.df, aes(x = xbar, y = 1-NoGo), color="red", size=.75)+
    geom_vline(xintercept = c(my.df$Delta.tv[1], my.df$Delta.lrv[1]), color="blue", linetype=2)+
    annotate("text", label = paste0("Min TPP = ", my.df$Delta.lrv[1]), x = my.df$Delta.lrv[1], 
             y = 0 + max(dpb$data[[2]]$y,dpb$data[[2]]$y)/nlines, size = tsize, colour = "black", hjust = 1)+
    annotate("text", label = paste("Base TPP = ", my.df$Delta.tv[1]), x = my.df$Delta.tv[1], 
             y = 0 + max(dpb$data[[2]]$y,dpb$data[[2]]$y)/nlines, size = tsize,
             colour = "black", hjust = 0)+
    labs(x = "Underlying treatment effect",
         y="Probability")+
    scale_x_continuous(expand = c(0,0), breaks=pretty(my.df$xbar, 10))+
    scale_y_continuous(expand = c(0,0), breaks=seq(0,1,.2), limits=c(0,1), label=scales::percent)+
    theme(panel.spacing.x = unit(6, "mm"), axis.text.x = element_text(angle=45, hjust=1,vjust=1))+
    labs(title="Operating characteristics as a function of treatment effect")
  
  table.plot2 <- ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ > Min TPP) > ", my.df$tau.lrv[1]*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ > Base TPP) > ", my.df$tau.tv[1]*100,"$%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ > Min TPP) $\\leq$ ", my.df$tau.ng[1]*100, "$%\\;\\;\\;\\,\\,$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ > Base TPP) $\\leq$ ", my.df$tau.tv[1]*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+ 
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")
  
  # plot_grid(main.plot, 
  #            table.plot2, ncol=1, align = "v", rel_heights=c(.78,.22), axis="l")
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
  
}

#' @rdname OneSampleNormalGamma
make.ss.ng.trt.oc2 <- function(my.df = get.ss.ng.trt.oc.df(), nlines=25, tsize=4){
  
  main.plot <-  my.df %>% mutate(Consider = 1 - Go - NoGo) %>% 
    gather(value= value, key=key, factor_key = T, c(Go, NoGo, Consider)) %>% 
    mutate(key=factor(key, c("Go", "Consider", "NoGo"))) %>%
    ggplot(aes(x=xbar, y=value, color=key))+
    geom_line(size=.75)
  dpb <- ggplot_build(main.plot)
  
  main.plot <- main.plot+
    labs(x = "Underlying treatment effect",
         y="Probability",
         color="Decision",
         title="Operating characteristics as a function of treatment effect")+
    theme(legend.position = "bottom")+
    scale_color_manual(values=c( "green", "grey20","red"))+
    scale_x_continuous(expand = c(0,0), breaks=pretty(my.df$xbar, 10))+
    scale_y_continuous(expand = c(0,0), breaks=seq(0,1,.2), limits=c(-0.0005,1), label=scales::percent)+
    theme(panel.spacing.x = unit(6, "mm"), axis.text.x = element_text(angle=45, hjust=1,vjust=1))+
    geom_vline(xintercept = c(my.df$Delta.tv[1], my.df$Delta.lrv[1]), color="blue", linetype=2)+
    annotate("text", label = paste0("Min TPP = ", my.df$Delta.lrv[1]), x = my.df$Delta.lrv[1], 
             y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize, colour = "black", hjust = 1)+
    annotate("text", label = paste("Base TPP = ", my.df$Delta.tv[1]), x = my.df$Delta.tv[1], 
             y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize,
             colour = "black", hjust = 0)+
    labs(x = "Underlying treatment effect",
         y="Probability")
  
  
  
  table.plot2 <- ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ > Min TPP) > ", my.df$tau.lrv[1]*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ > Base TPP) > ", my.df$tau.tv[1]*100,"$%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ > Min TPP) $\\leq$ ", my.df$tau.ng[1]*100, "$%\\;\\;\\;\\,\\,$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ > Base TPP) $\\leq$ ", my.df$tau.tv[1]*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+ 
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")
  
  # plot_grid(main.plot, 
  #            table.plot2, ncol=1, align = "v", rel_heights=c(.78,.22), axis="l")
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
  
}

#' @rdname OneSampleNormalGamma
get.ss.ng.ssize.oc.df <- function(mu.0.t = 3, n.0.t = 10, alpha.0.t = 0.25,
                                  beta.0.t = 1,
                                  s.t = 5, n.t = 50, n_LB_OC = floor(50*0.75),
                                  n_UB_OC=floor(50*2), npoints=15,
                                  Delta.lrv=2.5, Delta.tv=4, Delta.user = 3,
                                  tau.tv=.1, tau.lrv=.8, tau.ng=.65){
  # Create a grid
  specs <- expand.grid(n.t=floor(seq(n_LB_OC, n_UB_OC, length.out=npoints)),
                       Deltas.from=c(0, Delta.lrv, Delta.tv, Delta.user)) %>%
    mutate(mu.0.t = mu.0.t, n.0.t = n.0.t, alpha.0.t = alpha.0.t, beta.0.t = beta.0.t,
           s.t = s.t, n.t = n.t,
           Delta.lrv=Delta.lrv, Delta.tv=Delta.tv, Delta.user = Delta.user,
           tau.tv=tau.tv, tau.lrv=tau.lrv, tau.ng=tau.ng)


  for.plot <- bind_rows(apply(X = matrix(1:nrow(specs)), MARGIN = 1,
                              FUN = function(x){
                                # For each row of specs, get all possible probs along a fine grid
                                results <- get.ss.ng.df(mu.0.t = specs$mu.0.t[x], n.0.t = specs$n.0.t[x],
                                                        alpha.0.t = specs$alpha.0.t[x], beta.0.t = specs$beta.0.t[x],
                                                        s.t=specs$s.t[x], n.t=specs$n.t[x],
                                                        xbar.t=seq(min(min(specs$Delta.lrv, specs$Delta.user, 0)
                                                                       -specs$s.t/sqrt(specs$n.t)*5),
                                                                   max(max(specs$Delta.tv, specs$Delta.user, 0)
                                                                       +specs$s.t/sqrt(specs$n.t)*5), length.out=75),
                                                        Delta.lrv=specs$Delta.lrv[x], Delta.tv = specs$Delta.tv[x],
                                                        tau.tv=specs$tau.tv[x], tau.lrv=specs$tau.lrv[x],
                                                        tau.ng=specs$tau.ng[x])
                                # Identify max criteria for Go and min criteria for no-go
                                result.go <- results %>% dplyr::filter(result== "Go") %>% slice(1)
                                # Determine max number of TRT responders for No-Go
                                result.ng <- results %>% dplyr::filter(result== "No-Go") %>% slice(n())

                                refine.go <-  get.ss.ng.df(mu.0.t = specs$mu.0.t[x], n.0.t = specs$n.0.t[x],
                                                           alpha.0.t = specs$alpha.0.t[x], beta.0.t = specs$beta.0.t[x],
                                                           s.t=specs$s.t[x], n.t=specs$n.t[x],
                                                           xbar.t=seq(result.go$xbar - 0.25*result.go$s/sqrt(result.go$n),
                                                                      result.go$xbar + 0.25*result.go$s/sqrt(result.go$n),
                                                                      length.out=50),
                                                           Delta.lrv=specs$Delta.lrv[x], Delta.tv = specs$Delta.tv[x],
                                                           tau.tv=specs$tau.tv[x], tau.lrv=specs$tau.lrv[x],
                                                           tau.ng=specs$tau.ng[x])

                                refine.ng <-  get.ss.ng.df(mu.0.t = specs$mu.0.t[x], n.0.t = specs$n.0.t[x],
                                                           alpha.0.t = specs$alpha.0.t[x], beta.0.t = specs$beta.0.t[x],
                                                           s.t=specs$s.t[x], n.t=specs$n.t[x],
                                                           xbar.t=seq(result.ng$xbar - 0.25*result.ng$s/sqrt(result.ng$n),
                                                                      result.ng$xbar + 0.25*result.ng$s/sqrt(result.ng$n),
                                                                      length.out=50),
                                                           Delta.lrv=specs$Delta.lrv[x], Delta.tv = specs$Delta.tv[x],
                                                           tau.tv=specs$tau.tv[x], tau.lrv=specs$tau.lrv[x],
                                                           tau.ng=specs$tau.ng[x])
                                # Identify max criteria for Go and min criteria for no-go
                                result.go <- refine.go %>% dplyr::filter(result== "Go") %>% slice(1)
                                # Determine max number of TRT responders for No-Go
                                result.ng <- refine.ng %>% dplyr::filter(result== "No-Go") %>% slice(n())

                                my.df <- data.frame(Go = 1 - pnorm(q = result.go$xbar,
                                                                   mean = c(0, Delta.lrv, Delta.tv, Delta.user),
                                                                   sd=(result.go$s/sqrt(result.ng$n))),
                                                    NoGo = pnorm(q = result.ng$xbar,
                                                                 mean = c(0, Delta.lrv, Delta.tv, Delta.user),
                                                                 sd=result.ng$s/sqrt(result.ng$n)),
                                                    Delta=c(0, Delta.lrv, Delta.tv, Delta.user))  %>%
                                  mutate(n.t = specs$n.t[x], mu.0.t = specs$mu.0.t[x], n.0.t = specs$n.0.t[x],
                                         beta.0.t = specs$beta.0.t[x], s.t = specs$s.t[x],
                                         Delta.lrv=specs$Delta.lrv[x], Delta.tv=specs$Delta.tv[x],
                                         Delta.user = specs$Delta.user[x], tau.tv=specs$tau.tv[x],
                                         tau.lrv=specs$tau.lrv[x], tau.ng=specs$tau.ng[x])
                                my.df
                              }))

  for.plot$key <- factor(as.character(for.plot$Delta), as.character(for.plot$Delta)[1:4])
  levels(for.plot$key) <-  c(
    TeX("$\\Delta\\,$ = NULL = 0"),
    TeX(paste("$\\Delta\\,$ = Min TPP = $", Delta.lrv, "$ ")),
    TeX(paste("$\\Delta\\,$ = Base TPP = $", Delta.tv, "$ ")),
    TeX(paste("$\\Delta\\,$ = User defined = $", Delta.user, "$")))
  for.plot$key <- factor(for.plot$key, levels(for.plot$key)[order(c(0, Delta.lrv, Delta.tv,
                                                                    Delta.user))])

  for.plot
}


#' @rdname OneSampleNormalGamma
make.ss.ng.ssize.oc <- function(for.plot=get.ss.ng.ssize.oc.df(), nlines=25, tsize=4){
  
  main.plot <- ggplot() + 
    geom_line(data=for.plot, aes(x=n.t, y=Go), color="green") + 
    geom_line(data=for.plot, aes(x=n.t, y=1 - NoGo), color="red") +
    facet_wrap(~key, labeller="label_parsed")+
    geom_ribbon(data=for.plot, aes(x=n.t, ymin=0, ymax=Go), fill="green", alpha=.5)+
    geom_ribbon(data=for.plot, aes(x=n.t, ymin=Go, ymax=1-NoGo), fill="grey", alpha=.5)+
    geom_ribbon(data=for.plot, aes(x=n.t, ymin=1-NoGo, ymax=1), fill="red", alpha=.5)+
    labs(x="Total number of observed events", 
         y="Decision Probabilities", 
         title="Operating characteristics as a function of sample size",
         subtitle=TeX(paste0("Prior and nuisance parameter and standard deviation held fixed." )))+
    scale_x_continuous(expand=c(0,0))+
    scale_y_continuous(expand=c(0,0), breaks=seq(0,1,.2), minor_breaks=seq(0,1,.1), label=scales::percent)+
    theme(panel.spacing.x = unit(6, "mm"), axis.text.x = element_text(angle=45, hjust=1,vjust=1))
  
  
  table.plot2 <- ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ > Min TPP) > ", for.plot$tau.lrv[1]*100, "%$\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ > Base TPP) > $", for.plot$tau.tv[1]*100,"$%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ > Min TPP) $\\leq$ ", for.plot$tau.ng[1]*100, "%$\\;\\;\\;\\,\\,$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ > Base TPP) $\\leq$ ", for.plot$tau.tv[1]*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+ 
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits = c(0.75, .975), expand=c(0,0), labels = NULL, breaks = NULL)+
    scale_x_continuous(limits = c(-1, 3.5), expand=c(0,0), labels = NULL, breaks = NULL)+
    labs(x = "\n", y = "")
  
  # plot_grid(main.plot, 
  #            table.plot2, ncol=1, align = "v", rel_heights=c(.78,.22), axis="l")
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
}

# TwoSampleNormalGamma ------
#' @name TwoSampleNormalGamma
#' @title Functions to Support the Two Sample Continuous Scenario
#' @param mu.0.t prior mean for treatment group
#' @param n.0.t prior effective sample size for treatment group
#' @param alpha.0.t prior alpha parameter for treatment group
#' @param beta.0.t prior beta parameter for treatment group
#' @param xbar.t observed sample mean for treatment group
#' @param s.t observed sample standard deviation for treatment group
#' @param n.t sample size for treatment group
#' @param mu.0.c prior mean for control group
#' @param n.0.c prior effective sample size for control group
#' @param alpha.0.c prior alpha parameter for control group
#' @param beta.0.c prior beta parameter for control group
#' @param xbar.c observed sample mean for control group
#' @param s.c observed sample standard deviation for control group
#' @param n.c sample size for control group
#' @param Delta.lrv TPP Lower Reference Value aka Min TPP
#' @param Delta.tv TPP Target Value aka Base TPP
#' @param tau.tv threshold associated with Base TPP
#' @param tau.lrv threshold associated with Min TPP
#' @param tau.ng threshold associated with No-Go
#' @param seed random seed
#' @param nlines Control for text spacing
#' @param tsize Control for text size
#' @param nlines.ria Control for text spacing
#' @param seed random seed
#' @param n.MC Number of MC samples
#' @param Delta.LB Lower bound for treatment effect
#' @param Delta.UB Upper bound for treatment effect
#' @param ARatio Allocation ratio
#' @param npoints number of points
#' @param n_LB_OC Lower bound for sample size
#' @param n_UB_OC Upper bound for sample size
#' @examples
#' \dontrun{
#' make.ts.ng.ppp()
#' make.ts.ng.spp()
#' get.ts.ng.trt.oc.df()
#' make.ts.ng.trt.oc1()
#' make.ts.ng.trt.oc2()
#' get.ts.ng.ssize.oc.df()
#' make.ts.ng.ssize.oc()
#' }
#' @description
#'
#'   make.ts.ng.ppp: Make Two Sample Normal-Gamma Prior/Posterior Plot. Returns a ggplot object.
#'
#'   make.ts.ng.spp: Make Two Sample Normal-Gamma Shaded Posterior Plot.  Returns a graphic built using grid.arrange.
#'
#'   get.ts.ng.trt.oc.df: Get Two Sample Normal-Gamma Treatment Effect OC. Returns a data.frame.
#'
#'   make.ts.ng.trt.oc1: Make Two Sample Normal-Gamma Treatment Effect. Returns a graphic built using grid.arrange.
#'
#'   make.ts.ng.trt.oc2: Make Two Sample Normal-Gamma Treatment Effect.  Returns a graphic built using grid.arrange.
#'
#'   get.ts.ng.ssize.oc.df:  Get Two Sample Normal-Gamma sample size OC data.frame. Returns a data.frame.
#'
#'   make.ts.ng.ssize.oc: Make Two Sample Normal-Gamma Sample size OC plot
make.ts.ng.ppp <- function(mu.0.c = 0, alpha.c = .25, beta.c = 1, n.0.c = 1,
                           mu.0.t = 0, alpha.t = .25, beta.t = 1, n.0.t = 1,
                           xbar.c = 1.5, s.c = 4, n.c = 40,
                           xbar.t = 3, s.t = 4, n.t = 40) {
  
  pp.t <-   get.NG.post(mu.0 = mu.0.t, n.0 = n.0.t, alpha.0 = alpha.t, beta.0 = beta.t,
                        xbar = xbar.t, s = s.t, n = n.t, group = "Treatment")
  
  pp.c <-   get.NG.post(mu.0 = mu.0.c, n.0 = n.0.c, alpha.0 = alpha.c, beta.0 = beta.c,
                        xbar = xbar.c, s = s.c, n = n.c, group = "Control")
  
  tau.limits <- c(min(pmax(.1, qgamma(p = .005, shape = alpha.t, rate = beta.t)),
                      pmax(.1, qgamma(p = .005, shape = pp.t$alpha.n, rate = pp.t$beta.n)),
                      pmax(.1, qgamma(p = .005, shape = alpha.c, rate = beta.c)),
                      pmax(.1, qgamma(p = .005, shape = pp.c$alpha.n, rate = pp.c$beta.n))),
                  max(qgamma(p = .995,shape = alpha.t, rate = beta.t),
                      qgamma(p = .995,shape = pp.t$alpha.n, rate = pp.t$beta.n),
                      qgamma(p = .995,shape = alpha.c, rate = beta.c),
                      qgamma(p = .995,shape = pp.c$alpha.n, rate = pp.c$beta.n)
                  ))
  
  mu.limits <- c(min(qnorm(p = .005, mean = mu.0.t, sd = (alpha.t / beta.t) ^ (-1)/sqrt(n.0.t)),
                     qnorm(p = .005, mean = pp.t$mu.n, sd = (pp.t$alpha.n / pp.t$beta.n) ^ (-1)/sqrt(pp.t$n.n)),
                     qnorm(p = .005, mean = mu.0.c, sd = (alpha.c / beta.c) ^ (-1)/sqrt(n.0.c)),
                     qnorm(p = .005, mean = pp.c$mu.n, sd = (pp.c$alpha.n / pp.c$beta.n) ^ (-1)/sqrt(pp.c$n.n))
  ),
  max(qnorm(p = .995, mean = mu.0.t, sd = (alpha.t / beta.t) ^ (-1)/sqrt(n.0.t)),
      qnorm(p = .995, mean = pp.t$mu.n, sd = (pp.t$alpha.n / pp.t$beta.n) ^ (-1)/sqrt(pp.t$n.n)),
      qnorm(p = .995, mean = mu.0.c, sd = (alpha.c / beta.c) ^ (-1)/sqrt(n.0.c)),
      qnorm(p = .995, mean = pp.c$mu.n, sd = (pp.c$alpha.n / pp.c$beta.n) ^ (-1)/sqrt(pp.c$n.n))
  ))
  
  CON.prior <- expand.grid(
    tau = seq(tau.limits[1],tau.limits[2], length.out = 100),
    mu = seq(mu.limits[1],mu.limits[2], length.out = 100)) %>%
    mutate(n0 = pp.c$n.0, a0 = pp.c$alpha.0, b0 = pp.c$beta.0, category="Control", group="Prior", density="Control Prior",
           dens = dnorgam(mu = mu, tau = tau, mu0 = pp.c$mu.0, n0 = pp.c$n.0, a0 = pp.c$alpha.0, b0 = pp.c$beta.0))
  
  CON.post <- expand.grid(
    tau = seq(tau.limits[1],tau.limits[2], length.out = 100),
    mu = seq(mu.limits[1],mu.limits[2], length.out = 100)) %>%
    mutate(n0 = pp.c$n.n, a0 = pp.c$alpha.n, b0 = pp.c$beta.n, category="Control", group="Posterior", density="Control Posterior",
           dens = dnorgam(mu = mu, tau = tau, mu0 = pp.c$mu.n, n0 = pp.c$n.n, a0 = pp.c$alpha.n, b0 = pp.c$beta.n))
  
  TRT.prior <- expand.grid(
    tau = seq(tau.limits[1],tau.limits[2], length.out = 100),
    mu = seq(mu.limits[1],mu.limits[2], length.out = 100)) %>%
    mutate(n0 = pp.t$n.0, a0 = pp.t$alpha.0, b0 = pp.t$beta.0, category="Treatment", group="Prior", density="Treatment Prior",
           dens = dnorgam(mu = mu, tau = tau, mu0 = pp.t$mu.0, n0 = pp.t$n.0, a0 = pp.t$alpha.0, b0 = pp.t$beta.0))
  
  TRT.post <- expand.grid(
    tau = seq(tau.limits[1],tau.limits[2], length.out = 100),
    mu = seq(mu.limits[1],mu.limits[2], length.out = 100)) %>%
    mutate(n0 = pp.t$n.n, a0 = pp.t$alpha.n, b0 = pp.t$beta.n, category="Treatment", group="Posterior", density="Treatment Posterior",
           dens = dnorgam(mu = mu, tau = tau, mu0 = pp.t$mu.n, n0 = pp.t$n.n, a0 = pp.t$alpha.n, b0 = pp.t$beta.n))
  
  my.df <- rbind(CON.prior,CON.post,TRT.prior, TRT.post)
  my.df$color <-   as.numeric(cut((my.df$dens),150))
  my.colors <- colorRampPalette(c("grey80", "red",  "yellow"))(150)
  
  my.df$density <- factor(my.df$density, c("Control Prior", "Control Posterior", "Treatment Prior", "Treatment Posterior"))
  levels(my.df$density) <- c(
    paste0("Control Prior: NG(", round(mu.0.c,2), ", ", round(n.0.c,2), ", ", round(alpha.c,2), ", ", round(beta.c,2), ")"),
    paste0("Control Posterior: NG(", round(pp.c$mu.n), ", ", round(pp.c$n.n), ", ", round(pp.c$alpha.n), ", ", round(pp.c$beta.n), ")"),
    paste0("Treatment Prior: NG(", round(mu.0.t,2), ", ", round(n.0.t,2), ", ", round(alpha.t,2), ", ", round(beta.t,2), ")"),
    paste0("Treatment Posterior: NG(", round(pp.t$mu.n), ", ", round(pp.t$n.n), ", ", round(pp.t$alpha.n), ", ", round(pp.t$beta.n), ")")
  )
  
  joint.pdf <- ggplot(data= my.df, aes(x = mu, y = tau, fill = factor(color)))+
    geom_tile() + facet_wrap(~density)+
    scale_x_continuous(expand = c(0,0))+
    scale_y_continuous(expand = c(0,0),
                       breaks = pretty(c(min(my.df$tau),max(my.df$tau)), n=10), labels= round(pretty(c(min(my.df$tau),max(my.df$tau)), n=10)^ -.5,2),
                       sec.axis = sec_axis(~., name = TeX("Response precision, $\\tau$")))+
    scale_fill_manual(values= my.colors)+
    guides(fill=F)+
    labs(x = "Mean response",
         y = TeX("Response standard deviation, $\\sigma$"),
         title="Prior and posterior joint distributions for mean and precision parameters",
         subtitle=TeX(paste0("Control Data (n, $\\bar{x}$, s): (", round(n.c,2), ", ", round(xbar.c,2), ", ", round(s.c,2),"), ",
                             "Treatment Data: (", round(n.t,2), ", ", round(xbar.t,2), ", ", round(s.t,2),"). Observed treatment effect: $", 
                             round(xbar.t,2) - round(xbar.c,2), "$."))
    )
  
  
  marginal.densities <- rbind(
    gcurve(dt_ls(x,df = 2 * alpha.c, mu = mu.0.c, sigma = beta.c/(alpha.c * n.0.c)),
           from = mu.limits[1],
           to = mu.limits[2], n = 1001, category = "Control") %>% mutate(group="Prior", density="Control Prior"),
    gcurve(dt_ls(x,df = 2 * alpha.t, mu = mu.0.t, sigma = beta.t/(alpha.t * n.0.t)) ,
           from = mu.limits[1],
           to = mu.limits[2], n = 1001, category = "Treatment") %>% mutate(group="Prior", density="Treatment Prior"),
    gcurve(dt_ls(x,df = 2 * pp.c$alpha.n, mu = pp.c$mu.n, sigma = pp.c$beta.n/(pp.c$alpha.n * pp.c$ n.n )),
           from = mu.limits[1],
           to = mu.limits[2], n = 1001, category = "Control") %>% mutate(group="Posterior", density="Control Posterior"),
    gcurve(dt_ls(x,df = 2 * pp.t$alpha.n, mu = pp.t$mu.n, sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n )),
           from = mu.limits[1],
           to = mu.limits[2], n = 1001, category = "Treatment") %>% mutate(group="Posterior", density="Treatment Posterior")) %>%
    mutate(density=factor(density, c("Control Prior", "Control Posterior", "Treatment Prior", "Treatment Posterior")))
  
  
  
  levels(marginal.densities$density) <- c(
    paste0("Control Prior: t(", round(2 * alpha.c, 2), ", ", round(mu.0.c, 2), ", ", round(beta.c/(alpha.c * n.0.c),2), ")"),
    paste0("Control Posterior: t(",  round(2 * pp.c$alpha.n, 2), ", ", round(pp.c$mu.n, 2), ", ", round(pp.c$beta.n/(pp.c$alpha.n * pp.c$ n.n ),2), ")"),
    paste0("Treatment Prior t(",  round(2 * alpha.t, 2), ", ", round(mu.0.t, 2), ", ", round(beta.t/(alpha.t * n.0.t),2), ")"),
    paste0("Treatment Posterior: t(",  round(2 * pp.t$alpha.n, 2), ", ", round(pp.t$mu.n, 2), ", ", round(pp.t$beta.n/(pp.t$alpha.n * pp.t$ n.n ),2), ")")
  )
  
  marginal.pdf <- ggplot(data=marginal.densities, aes(x = x,y = y, color = category))+
    geom_line(size=.75)+
    theme(legend.position = "bottom")+
    labs(title="Marginal distributions for the means",color = NULL,
         x = "Mean response", y="")+facet_wrap(~density)+
    scale_y_continuous(breaks=NULL, labels = NULL) + guides(color=F)
  
  return(grid.arrange(joint.pdf, marginal.pdf, ncol=1))
  
}

#' @rdname TwoSampleNormalGamma
make.ts.ng.spp <- function(mu.0.c = 0, alpha.c = .25, beta.c = 1, n.0.c = 1,
                           mu.0.t = 0, alpha.t = .25, beta.t = 1, n.0.t = 1,
                           xbar.c = 1.5, s.c = 4, n.c = 40,
                           xbar.t = 26, s.t = 4, n.t = 40,
                           Delta.lrv = 1, Delta.tv = 1.5,
                           tau.ng = .65, tau.lrv = .8, tau.tv = .1,
                           seed = 1234, n.MC = 1000,
                           nlines=25, nlines.ria = 20, tsize=4){
  set.seed(seed=seed)
  # Smart two-stage search-----
  my.means <- seq(-Delta.tv*4, Delta.tv*4, length.out=50)
  stage1 <- get.ts.ng.mc.df(mu.0.c = mu.0.c, n.0.c = n.0.c, alpha.0.c = alpha.c,
                            beta.0.c = beta.c,
                            xbar.c = xbar.c, s.c = s.c, n.c = n.c, group.c = "Control",
                            mu.0.t = mu.0.t, n.0.t = n.0.t, alpha.0.t = alpha.t,
                            beta.0.t = beta.t,
                            xbar.t = xbar.c + my.means, s.t = s.t, n.t = n.t,
                            group.t="treat",
                            Delta.tv = Delta.tv,Delta.lrv = Delta.lrv,tau.tv = tau.tv,
                            tau.lrv = tau.lrv,tau.ng = tau.ng, n.MC = n.MC)
  stage1.go <- stage1 %>% dplyr::filter(result== "Go") %>% slice(1)
  # Determine max number of TRT responders for No-Go
  stage1.ng <- stage1 %>% dplyr::filter(result== "No-Go") %>% slice(n())

  my.means <- c(seq(stage1.go$xbar.t - stage1.go$s.t/4, stage1.go$xbar.t, length.out=100),
                seq(stage1.ng$xbar.t , stage1.ng$xbar.t + stage1.ng$s.t/4, length.out=100))
  stage2 <- get.ts.ng.mc.df(mu.0.c = mu.0.c, n.0.c = n.0.c, alpha.0.c = alpha.c,
                            beta.0.c = beta.c,
                            xbar.c = xbar.c, s.c = s.c, n.c = n.c, group.c = "Control",
                            mu.0.t = mu.0.t, n.0.t = n.0.t, alpha.0.t = alpha.t,
                            beta.0.t = beta.t,
                            xbar.t = my.means, s.t = s.t, n.t = n.t, group.t="treat",
                            Delta.tv = Delta.tv,Delta.lrv = Delta.lrv,tau.tv = tau.tv,
                            tau.lrv = tau.lrv,tau.ng = tau.ng,
                            n.MC = n.MC)
  result.go <- stage2 %>% dplyr::filter(result== "Go") %>% slice(1)
  result.ng <- stage2 %>% dplyr::filter(result== "No-Go") %>% slice(n())


  # Get post parameters
  pp.c <- get.NG.post(mu.0 = mu.0.c, n.0 = n.0.c, alpha.0 = alpha.c, beta.0 = beta.c,
                      xbar = xbar.c, s = s.c, n = n.c, group="Control")
  pp.t <- get.NG.post(mu.0 = mu.0.t, n.0 = n.0.t, alpha.0 = alpha.t, beta.0 = beta.t,
                      xbar = xbar.t, s = s.t, n = n.t, group="Treatment")

  # Get samples
  samples <- data.frame(samples =
                          rt_ls(n = 100000, df = 2 * pp.t$alpha.n, mu = pp.t$mu.n,
                                sigma = pp.t$beta.n/(pp.t$alpha.n * pp.t$n.n))-
                          rt_ls(n = 100000, df = 2 * pp.c$alpha.n, mu = pp.c$mu.n,
                                sigma = pp.c$beta.n/(pp.c$alpha.n * pp.c$n.n))
  )
  samples$group <- paste0(
    paste0("Control Prior: NG(",  round(pp.c$mu.0, 2), ", ", round(pp.c$n.0, 2), ", ",
           round(pp.c$alpha.0,2), ", ", round(pp.c$beta.0,2), ")"),
    paste0("; Treatment Prior: NG(",  round(pp.t$mu.0, 2), ", ", round(pp.t$n.0, 2), ", ",
           round(pp.t$alpha.0,2), ", ", round(pp.t$beta.0,2), ")")
  )

  P.R1 = mean(samples$samples > Delta.lrv)
  P.R3 = mean(samples$samples > Delta.tv)

  # comput result -----

  result = ifelse(P.R1 > tau.lrv & P.R3 > tau.tv, "Go",
                  ifelse(P.R1 <= tau.ng & P.R3 <= tau.tv, "No-Go",
                         "Consider"))
  result.color <- ifelse(result=="Go", "darkgreen", ifelse(result=="No-Go",
                                                           "red", "black"))

  # Subtitle ----

  for.subtitle <- TeX(paste0("Control, Treatment Data (n, $\\bar{x}_{T}$, s): (",
                             round(n.c,2), ", ", round(xbar.c,2), ", ", round(s.c,2),"), ",
                             "(", round(n.t,2), ", ", round(xbar.t,2), ", ", round(s.t,2),
                             "); $\\bar{x}_{T}$ needed for Go: $",round(result.go$xbar.t[1],2),
                             "$; for No-Go: $", round(result.ng$xbar.t[1],2) , "$"))

  # Annotation line 1: Decision ----
  for.decision <- paste0("Decision: ", result)

  # Annotation line 2: Decision interval----

  for.decision.interval <- ifelse(mean(samples$samples > Delta.tv) > tau.tv,
                                  paste0("Decision interval: (",
                                         round(quantile(samples$samples,
                                                        1 - tau.lrv), 2), ", ",
                                         round(quantile(samples$samples,
                                                        1 - tau.tv), 2), ")"),
                                  paste0("Decision interval: (",
                                         round(quantile(samples$samples,
                                                        1 - tau.ng), 2), ", ",
                                         round(quantile(samples$samples,
                                                        1 - tau.tv), 2), ")"))

  # Annotation line 3: P(Delta >= Min TPP)

  annotate.P1 <- ifelse(result=="Go",
                        TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) = $",
                                   round(mean(samples$samples > Delta.lrv),4)*100,
                                   "%$ > $", tau.lrv*100,"$%")),
                        ifelse(result=="No-Go",
                               TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) = $",
                                          round(mean(samples$samples > Delta.lrv),4)*100 ,
                                          "%$ $\\leq$ $", tau.ng*100,"$%")),
                               TeX(paste0("$P(\\Delta\\, \\geq$ Min TPP$) = ",
                                          round(mean(samples$samples > Delta.lrv),4)*100,
                                          "%$" ))))

  # Annotation line 3: P(Delta >= Base TPP)

  annotate.P2 <- ifelse(result=="Go",
                        TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) = $",
                                   round(mean(samples$samples > Delta.tv),4)*100 ,
                                   "%$ $\\geq$", tau.tv*100, "%")),
                        ifelse(result=="No-Go",
                               TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) = $",
                                          round(mean(samples$samples > Delta.tv),
                                                4)*100
                                          , "%$ <", tau.tv*100, "%")),
                               TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) = $",
                                          round(mean(samples$samples > Delta.tv),
                                                4)*100, "%$"  ))))

  dplot <- ggplot() + geom_density(data = samples, aes(x = samples))+
    facet_wrap(~group)
  # Access the ggplot to get goodies to help accomplish shading
  dpb <- ggplot_build(dplot)

  x1.1 <- min(which(dpb$data[[1]]$x >=Delta.lrv))
  x2.1 <- max(which(dpb$data[[1]]$x <=Delta.tv))+1
  x1.2 <- pmin(min(which(dpb$data[[1]]$x >=Delta.tv)),
               max(which(dpb$data[[1]]$x <=Inf)))
  x2.2 <- max(which(dpb$data[[1]]$x <=Inf))

  go.segment <- data.frame(x = quantile(samples$samples, 1 - tau.lrv),
                           xend =quantile(samples$samples, 1 - tau.tv),
                           y= min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                           yend = min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3) %>%
    mutate(group= paste0(
      paste0("Control Prior: NG(",  round(pp.c$mu.0, 2), ", ", round(pp.c$n.0, 2), ", ",
             round(pp.c$alpha.0,2), ", ", round(pp.c$beta.0,2), ")"),
      paste0("; Treatment Prior: NG(",  round(pp.t$mu.0, 2), ", ", round(pp.t$n.0, 2),
             ", ", round(pp.t$alpha.0,2), ", ", round(pp.t$beta.0,2), ")")
    ))

  nogo.segment <- data.frame(x = quantile(samples$samples, 1 - tau.ng),
                             xend =quantile(samples$samples, 1 - tau.tv),
                             y= min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                             yend = min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3) %>%
    mutate(group= paste0(
      paste0("Control Prior: NG(",  round(pp.c$mu.0, 2), ", ", round(pp.c$n.0, 2), ", ",
             round(pp.c$alpha.0,2), ", ", round(pp.c$beta.0,2), ")"),
      paste0("; Treatment Prior: NG(",  round(pp.t$mu.0, 2), ", ", round(pp.t$n.0, 2),
             ", ", round(pp.t$alpha.0,2), ", ", round(pp.t$beta.0,2), ")")
    ))

  # Custom graphic parameters
  x.limits <- c(min(min(dpb$data[[1]]$x), Delta.lrv), max(max(dpb$data[[1]]$x),
                                                          Delta.tv))
  ticks <- pretty(x = c(x.limits[1], x.limits[2]), n = 15)

  # When the mound is on the right... we want to display annotation on the left
  if(which(dpb$data[[1]]$y== max(dpb$data[[1]]$y)) > length(dpb$data[[1]]$x)/2){
    annotate.x = min(min(dpb$data[[1]]$x), Delta.lrv)
    annotate.j = 0
  } else {annotate.x = max(max(dpb$data[[1]]$x), Delta.tv)
  annotate.j = 1}

  # Introduce shading ----
  main.plot <- dplot +
    # NOTE: GREG SHUT THIS OPTION OFF 3/27/20 BECAUSE LARGE XBAR_TRT VALUES WOULD CAUSE GRAPHIC TO CRASH
    # ALSO ISSUE IS RELATED TO SHADING WHICH HAD BEEN SHUT OFF WITH CODE IN PLACE BY SETTING FILL ALPHA TO 0.
    # geom_area(data = data.frame(x = dpb$data[[1]]$x[1:x1.1],
    #                             y = dpb$data[[1]]$y[1:x1.1]),
    #           aes(x = x, y = y), fill=alpha("red", 0))+
    # geom_area(data = data.frame(x = dpb$data[[1]]$x[x1.1:x2.1],
    #                             y = dpb$data[[1]]$y[x1.1:x2.1]),
    #           aes(x = x, y = y), fill=alpha("grey80", 0))+
    # geom_area(data = data.frame(x = dpb$data[[1]]$x[x1.2:x2.2],
    #                             y = dpb$data[[1]]$y[x1.2:x2.2]),
    #           aes(x = x, y = y), fill=alpha("green", 0))+
    scale_x_continuous(limits = x.limits, breaks=pretty(x=x.limits, n=10))+
    scale_y_continuous(breaks=NULL)

  # Add Annotations -----
  main.plot <- main.plot +
    labs(title = TeX("Posterior Density of Treatment Effect"),
         x = TeX("$\\Delta\\,$ = Mean treatment difference (Treatment - Control)"),
         y = NULL,
         subtitle=for.subtitle)+
    annotate("text", label = for.decision,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 0,
             size = tsize+1, colour = result.color, hjust = annotate.j)+
    annotate("text", label = for.decision.interval,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 1,
             size = tsize, colour = result.color, hjust = annotate.j)+
    annotate("text", label = annotate.P1, color=result.color,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 2,
             size = tsize,  hjust = annotate.j)+
    annotate("text", label = annotate.P2, color=result.color,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 3,
             size = tsize,  hjust = annotate.j)

  # Add reference lines and Credible interval
  if(mean(samples$samples > Delta.tv) > tau.tv) {
    main.plot <- main.plot +
      geom_segment(data=go.segment, aes(x=x, xend=xend, y=y, yend=yend),
                   arrow = arrow(ends="both", angle = 90), color=result.color, size=.75)
  } else {
    main.plot <- main.plot +
      geom_segment(data=nogo.segment, aes(x=x, xend=xend, y=y, yend=yend, group=group),
                   arrow = arrow(ends="both", angle = 90), color=result.color, size=.75)

  }
  main.plot <- main.plot +
    geom_vline(xintercept = c(Delta.lrv, Delta.tv), linetype = 2, color = c("blue", "blue")) +
    annotate("text", label = TeX(paste0("$", Delta.lrv,"$ = Min TPP")), x = Delta.lrv, y = 0 + max(dpb$data[[1]]$y)/nlines.ria, size = tsize, colour = "black", hjust = 1)+
    annotate("text", label = TeX(paste0("Base TPP = $",  Delta.tv,"$")), x = Delta.tv, y = 0 + max(dpb$data[[1]]$y)/nlines.ria, size = tsize, colour = "black", hjust = 0)


  table.plot2 <- ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) > ", tau.lrv*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) > ", tau.tv*100,"%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Min TPP) $\\leq$ $", tau.ng*100, "$$%\\;\\;\\;\\;$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ $\\geq$ Base TPP) $\\leq$ $", tau.tv*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")

  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))

}

#' @rdname TwoSampleNormalGamma
get.ts.ng.mc.df <- function(
  mu.0.c = 0, n.0.c = 10, alpha.0.c=.25 * 4, beta.0.c = 1 * 4,
  xbar.c = seq(-3,3,length.out=20), s.c = 3, n.c = 25, group.c="Control",
  mu.0.t = 0, n.0.t = 10, alpha.0.t=.25 * 4, beta.0.t = 1 * 4,
  xbar.t = seq(0, 6, length.out=20), s.t = 2, n.t = 25, group.t="Treatment",
  Delta.tv = 1.75, Delta.lrv = 1.5,  tau.tv = 1, tau.lrv = 1, tau.ng = .65,
  n.MC = 1000, seed=1234){

  # Create a simulation grid
  my.grid <- expand.grid(
    mu.0.c = mu.0.c, n.0.c = n.0.c, alpha.0.c = alpha.0.c, beta.0.c = beta.0.c,
    xbar.c = xbar.c, s.c = s.c, n.c = n.c,
    mu.0.t = mu.0.t, n.0.t = n.0.t, alpha.0.t = alpha.0.t, beta.0.t = beta.0.t,
    xbar.t = xbar.t, s.t = s.t, n.t = n.t,
    Delta.tv = Delta.tv, Delta.lrv = Delta.lrv,
    tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng,
    n.MC = n.MC, seed = seed)

  # Create a function that makes a single evaluation
  my.function1 <- function(
    mu.0.c = my.grid$mu.0.c, n.0.c = my.grid$n.0.c, alpha.0.c = my.grid$alpha.0.c,
    beta.0.c = my.grid$beta.0.c,
    xbar.c = my.grid$xbar.c, s.c = my.grid$s.c, n.c = my.grid$n.c,
    mu.0.t = my.grid$mu.0.t, n.0.t = my.grid$n.0.t, alpha.0.t = my.grid$alpha.0.t,
    beta.0.t = my.grid$beta.0.t,
    xbar.t = my.grid$xbar.t, s.t = my.grid$s.t, n.t = my.grid$n.t,
    Delta.tv = my.grid$Delta.tv, Delta.lrv = my.grid$Delta.lrv,
    tau.tv = my.grid$tau.tv, tau.lrv = my.grid$tau.lrv, tau.ng = my.grid$tau.ng,
    seed = my.grid$seed, n.MC = my.grid$n.MC){

    return(get.ts.ng.mc(mu.0.c = mu.0.c, n.0.c = n.0.c, alpha.0.c = alpha.0.c,
                        beta.0.c = beta.0.c,
                        xbar.c = xbar.c, s.c = s.c, n.c = n.c, group.c="Control",
                        mu.0.t = mu.0.t, n.0.t = n.0.t, alpha.0.t = alpha.0.t,
                        beta.0.t = beta.0.t,
                        xbar.t = xbar.t, s.t = s.t, n.t = n.t, group.t="Treatment",
                        Delta.tv = Delta.tv, Delta.lrv = Delta.lrv, seed = seed,
                        # This next line was missing, so defaults were being used.
                        tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng,
                        n.MC = n.MC))
  }

  # Apply the function to each row of the grid
  my.grid <- cbind(my.grid, plyr:::list_to_dataframe(do.call(Map, c(f= my.function1,
                                                                    my.grid)))) %>%
    mutate(result = factor(result, c("Go", "Consider", "No-Go")))
  return(my.grid)
}

#' @rdname TwoSampleNormalGamma
get.ts.ng.mc <- function(mu.0.c = 1.5, n.0.c = 1, alpha.0.c=.25 ,
                         beta.0.c = 1 ,
                         xbar.c = 1.5, s.c = 1, n.c = 25, group.c="Control",
                         mu.0.t = 1.75, n.0.t = 1, alpha.0.t=.25,
                         beta.0.t = 1 ,
                         xbar.t = 0, s.t = 1, n.t = 25, group.t="Treatment",
                         Delta.tv = 1.5, Delta.lrv = 1, tau.tv = 1,
                         tau.lrv = 1, tau.ng = 0.35,
                         seed = 1234, n.MC = 5000)
{
  # Compute individual posterior distribution parameters
  CON.results <- get.NG.post(mu.0 = mu.0.c, n.0 = n.0.c, alpha.0 = alpha.0.c,
                             beta.0 = beta.0.c, xbar = xbar.c, s = s.c,
                             n = n.c, group = group.c)
  TRT.results <- get.NG.post(mu.0 = mu.0.t, n.0 = n.0.t, alpha.0 = alpha.0.t,
                             beta.0 = beta.0.t, xbar = xbar.t, s = s.t,
                             n = n.t, group = group.t)

  if(is.null(seed) ==F) set.seed(seed = seed)
  # Estimating prior probabilities and posterior probabilities
  Z.0 <- rt_ls(n = n.MC, df = TRT.results$tdf.0, mu = TRT.results$mu.0,
               sigma = TRT.results$sigma.0) - rt_ls(n = n.MC, df = CON.results$tdf.0,
                                                    mu = CON.results$mu.0,
                                                    sigma = CON.results$sigma.0)
  Z.n <- rt_ls(n = n.MC, df = TRT.results$tdf.n, mu = TRT.results$mu.n,
               sigma = TRT.results$sigma.n) -
    rt_ls(n = n.MC, df = CON.results$tdf.n, mu = CON.results$mu.n,
          sigma = CON.results$sigma.n)

  #### COMPARE HERE: I synced this with TS binomial
  # TIANYU: I made the following changes:
  # 1. Reporting P(Delta > Min.tpp), P(Delta > Base.tpp)
  # 2. I had to change P.R1 >= tau.lrv & P.R3 >= tau.tv because MC step encounters cases where 1 is returned. 
  my.df <- data.frame(P.R1 = mean(Z.n > Delta.lrv),
                      P.R3 = mean(Z.n > Delta.tv)) %>%
    mutate(result = ifelse(P.R1 > tau.lrv & P.R3 > tau.tv, "Go",
                           ifelse(P.R3 <= tau.tv & P.R1 <= tau.ng, "No-Go", "Consider")))
  
  
  return(my.df)
}

#' @rdname TwoSampleNormalGamma

# This must be deprecated!!!  Its not used elsewhere
get.ts.ng.decision <-   function(mu.0.c = 1.5, n.0.c = 10, alpha.0.c=.25 * 4, beta.0.c = 1 * 4,
           xbar.c = 1.5, s.c = 1, n.c = 25, group.c="Control",
           mu.0.t = 1.75, n.0.t = 10, alpha.0.t=.25 * 4, beta.0.t = 1 * 4,
           xbar.t = 1.85, s.t = 1, n.t = 25, group.t="Treatment",
           Delta.tv = 0.5, Delta.lrv = 0.25, tau.tv = 1, tau.lrv = 1,
           tau.ng = 0.35,
           seed = 1234, n.MC = 1000)  {

    from.here=Delta.lrv-s.t/4
    to.here=Delta.tv+s.t/2
    npoints=20

    # This grid really just runs over a sequence of values for underlying treatment effect
    #
    my.grid <- expand.grid(mu.0.c = mu.0.c, n.0.c = n.0.c, alpha.0.c=alpha.0.c,
                           beta.0.c = beta.0.c,
                           xbar.c = xbar.c, s.c = s.c, n.c = n.c, group.c="Control",
                           mu.0.t = mu.0.t, n.0.t = n.0.t, alpha.0.t=alpha.0.t,
                           beta.0.t = beta.0.t,
                           xbar.t = xbar.c + seq(from.here, to.here, npoints),
                           s.t = s.t, n.t = n.t, group.t="Treatment",
                           Delta.tv = Delta.tv, Delta.lrv = Delta.lrv, tau.tv = tau.tv,
                           tau.lrv = tau.lrv, tau.ng = tau.ng,
                           seed = seed, n.MC = n.MC)

    # After running over the grid, posterior probabiltities are calculated.
    # For each value of underlying treatment effect (group_by(x)), the relative frequency of Go, Consider and No-Go is computed
    # The complement of nogo is also computed

    # time:
    start.time <- Sys.time()
    my.list <- apply(X=matrix(1:nrow(my.grid)), MARGIN = 1,
                     FUN = function(x) decision.ts.continuous(as.list(my.grid[x,])))
    elapsed <- Sys.time() - start.time; print(elapsed)
    my.list <- parApply(cl = cl, X=matrix(1:nrow(my.grid)), MARGIN = 1,
                        FUN = function(x) decision.ts.continuous(as.list(my.grid[x,])))


    plot.df <- plyr:::list_to_dataframe(my.list)

    plot.df <- plot.df %>%
      mutate(
        P.R1 = pt_ls(x = Delta.lrv, df = tdf.n, mu = mu.n, sigma = sigma.n),
        P.R3 = 1 - pt_ls(Delta.tv, df = tdf.n, mu = mu.n, sigma = sigma.n),
        result = ifelse(P.R1 <= tau.lrv & P.R3 > tau.tv, "Go",
                        ifelse(P.R1 > tau.ng  & P.R3 < tau.tv, "No-Go", "Consider"))) %>%
      group_by(x) %>%
      summarize(Go=mean(result=="Go"), Consider=mean(result=="Consider"), NoGo=mean(result=="No-Go")) %>%
      mutate(NoGo2=1-NoGo)


    my.plot <- plot.df  %>%
      ggplot() + geom_line(aes(x = x, y = Go), color="green")+
      geom_line(aes(x = x, y = NoGo2), color="red")
    dpb <- ggplot_build(my.plot)

    my.2plot <- grid.arrange(
      my.plot+
        geom_area(data = dpb$data[[1]], aes(x = x, y = y), fill=alpha("lightgreen", .5))+
        geom_ribbon(data = dpb$data[[2]], aes(x = x, ymin = y, ymax = 1), fill=alpha("red", .5))+
        geom_line(data = plot.df, aes(x = x, y = Go), color="green", size=.75)+
        geom_line(data = plot.df, aes(x = x, y = NoGo2), color="red", size=.75)+
        geom_vline(xintercept = c(Delta.tv, Delta.lrv), color="blue", linetype=2)+
        annotate("text", label = TeX(paste0(Delta.lrv,"$=\\Delta_{Min}$")), x = Delta.lrv, y = 0 + max(dpb$data[[1]]$y)/nlines, size =
                   tsize, colour = "black", hjust = 1)+
        annotate("text", label = TeX(paste("$\\Delta_{Base}=$", Delta.tv)), x = Delta.tv, y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize,
                 colour = "black", hjust = 0)+
        labs(x = TeX("$\\Delta\\,$ = Underlying treatment effect"),
             y="Probability")+
        scale_x_continuous(expand = c(0,0), breaks=pretty(plot.df$x, 10))+
        scale_y_continuous(expand = c(0,0), breaks=seq(0,1,.2), limits=c(0,1))+
        labs(title="Decision Curves View 1",
             subtitle="y-axis plots Go probabilities and the complement of No-Go probabilities\n"),
      plot.df %>% gather(key = key, value=value, -x,factor_key = T) %>%
        dplyr::filter(key !="NoGo2") %>%
        ggplot(aes(x=x, y=value, color=key))+
        geom_line(size=.75)+
        geom_point()+
        labs(x = TeX("$\\Delta\\,$ = Underlying treatment effect"),
             y="Probability",
             color="Decision",
             title="Decision Curves View 2",
             subtitle="y-axis plots Go, Consider, and No-Go probabilities")+
        theme(legend.position = "bottom")+
        scale_color_manual(values=c( "green", "grey20","red"))+
        scale_x_continuous(expand = c(0,0), breaks=pretty(plot.df$x, 10))+
        scale_y_continuous(expand = c(0,0), breaks=seq(0,1,.2), limits=c(0,1))+
        geom_vline(xintercept = c(Delta.tv, Delta.lrv), color="blue", linetype=2)+
        annotate("text", label = TeX(paste0(Delta.lrv,"$=\\Delta_{Min}$")), x = Delta.lrv, y = 0 + max(dpb$data[[1]]$y)/nlines, size =
                   tsize, colour = "black", hjust = 1)+
        annotate("text", label = TeX(paste("$\\Delta_{Base}=$", Delta.tv)), x = Delta.tv, y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize,
                 colour = "black", hjust = 0)

    )

    return(list(my.2plot, plot.df))
  }

#' @rdname TwoSampleNormalGamma
get.ts.ng.trt.oc.df <- function(mu.0.c = 0, n.0.c = 1, alpha.0.c=.25, beta.0.c = 1,
                                xbar.c = 1.5, s.c = 2,  group.c="Control",
                                mu.0.t = 3.75, n.0.t = 1, alpha.0.t=.25, beta.0.t = 1,
                                xbar.t = 1.85, s.t = 2,  group.t="Treatment",
                                Delta.LB=0, Delta.UB=5, ARatio=1, N=50,
                                Delta.tv = 2.5, Delta.lrv = 1.5, Delta.user = 4,
                                tau.tv = 1, tau.lrv = 1, tau.ng = .65,
                                npoints=5, n.MC = 500,
                                seed = 1234){
  cat("\nstarting decision.ts.continuous\n")
  my.specs = data.frame(treatment.effect=seq(Delta.LB, Delta.UB, length.out=npoints))
  my.specs$mu.0.c <- mu.0.c
  my.specs$n.0.c = n.0.c
  my.specs$alpha.0.c=alpha.0.c
  my.specs$beta.0.c = beta.0.c
  my.specs$xbar.c = xbar.c
  my.specs$s.c = s.c
  my.specs$n.c = N - floor(N*(ARatio/(1+ARatio)))
  my.specs$mu.0.t = mu.0.t
  my.specs$n.0.t = n.0.t
  my.specs$alpha.0.t=alpha.0.t
  my.specs$beta.0.t = beta.0.t
  my.specs$s.t = s.t
  my.specs$n.t = floor(N*(ARatio/(1+ARatio)))
  my.specs$Delta.tv = Delta.tv
  my.specs$Delta.lrv = Delta.lrv
  my.specs$tau.tv = tau.tv
  my.specs$tau.lrv = tau.lrv
  my.specs$tau.ng = tau.ng
  my.specs$n.MC = n.MC
  my.specs$npoints = npoints
  my.specs$seed = seed

  # Fixed placebo data and all priors + observed standard deviation and sample size.

  my.func <- function(x){
    #cat("my.treatmenteffect start")
    # TIANYU: This previously used: sd=my.specs$s.t[x]
    # I now divide by sqrt of treatment sample size
    my.treatmeant.effect <- rnorm(n = my.specs$n.MC[x],
                                  mean = my.specs$treatment.effect[x],
                                  sd=my.specs$s.t[x]/sqrt(my.specs$n.t[x]))

    for.return <- get.ts.ng.mc.df(mu.0.c = my.specs$mu.0.c[x], n.0.c = my.specs$n.0.c[x],
                                  alpha.0.c=my.specs$alpha.0.c[x], beta.0.c = my.specs$beta.0.c[x],
                                  xbar.c = my.specs$xbar.c[x], s.c = my.specs$s.c[x],
                                  n.c = my.specs$n.c[x], group.c="Control",
                                  mu.0.t = my.specs$mu.0.t[x], n.0.t = my.specs$n.0.t[x],
                                  alpha.0.t=my.specs$alpha.0.t[x], beta.0.t = my.specs$beta.0.t[x],
                                  xbar.t = my.specs$xbar.c[x] + my.treatmeant.effect ,
                                  s.t = my.specs$s.t[x], n.t = my.specs$n.t[x], group.t="Treatment",
                                  Delta.tv = my.specs$Delta.tv[x], Delta.lrv = my.specs$Delta.lrv[x],
                                  tau.tv = my.specs$tau.tv[x], tau.lrv = my.specs$tau.lrv[x],
                                  tau.ng = my.specs$tau.ng[x],
                                  n.MC = my.specs$n.MC[x]) %>%
      mutate(treatment.effect = my.specs$treatment.effect[x])
    #cat("\nget.NG.ts.df.mc finished\n")
    return(for.return)
  }

  # Add parallel parApply here
   for.plot <- bind_rows(apply(X = matrix(1:nrow(my.specs)), MARGIN = 1, FUN = my.func)) %>%
  #for.plot <- bind_rows(apply(X = matrix(1:1), MARGIN = 1, FUN = my.func)) %>%
    
   mutate(result = factor(result, c("Go", "Consider", "No-Go"))) %>%
    group_by(treatment.effect, result) %>%
    summarize(N = n()) %>%
    mutate(freq = N / sum(N)) %>% dplyr::select(result, treatment.effect, freq) %>%
    ungroup() %>%
    complete(result, fill = list(N = 0, freq = 0))%>%
    mutate(mu.0.c = mu.0.c, n.0.c = n.0.c, alpha.0.c=alpha.0.c, beta.0.c = beta.0.c,
           xbar.c = xbar.c, s.c = s.c,  group.c="Control",
           mu.0.t = mu.0.t, n.0.t = n.0.t, alpha.0.t=alpha.0.t, beta.0.t = beta.0.t,
           xbar.t = xbar.t, s.t = s.t,  group.t="Treatment",
           Delta.LB=Delta.LB, Delta.UB=Delta.UB, ARatio=ARatio, N=N,
           Delta.lrv = Delta.lrv, Delta.tv = Delta.tv, Delta.user = Delta.user,
           tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng,
           npoints=npoints, n.MC = n.MC,
           seed = seed)
  return(for.plot)
}

# check1 <- get.ts.ng.trt.oc.df()
# make.ts.ng.trt.oc1(for.plot=check1)


#' @rdname TwoSampleNormalGamma
make.ts.ng.trt.oc1 <- function(for.plot=get.ts.ng.trt.oc.df(), nlines=20, tsize=4){
  
  for.plot1 <- for.plot %>%
    dplyr::filter(result!="Consider") %>% 
    mutate(freq=ifelse(result=="Go", freq, 1 - freq)) 
  
  my.plot <- ggplot() + 
    geom_line(data=for.plot1 %>% dplyr::filter(result=="Go"), aes(x=treatment.effect, y=freq), color="darkgreen") 
  dpb <- ggplot_build(my.plot)
  
  # Simply takes a dataframe from return.ssize.df and plots
  main.plot <- my.plot +
    geom_ribbon(data=for.plot1 %>% dplyr::filter(result=="Go"), aes(x=treatment.effect, ymin=0, ymax=freq), fill="darkgreen")+
    geom_ribbon(data=for.plot1 %>% dplyr::filter(result=="Go"), aes(x=treatment.effect, ymin=freq, ymax=1), fill="grey")+
    geom_line(data=for.plot1 %>% dplyr::filter(result=="No-Go"), aes(x=treatment.effect, y=freq), color="red")+
    geom_ribbon(data=for.plot1 %>% dplyr::filter(result=="No-Go"), aes(x=treatment.effect, ymin=freq, ymax=1), fill="red")+
    labs(x="Underlying Treatment effect", 
         y="Probability", 
         title="Operating characteristics as a function of treatment effect",
         subtitle = paste0("Total randomized with ", 
                           for.plot$N[1], 
                           " with ", 
                           for.plot$N[1] - floor(for.plot$N[1]*(for.plot$ARatio[1]/(1+for.plot$ARatio[1]))), 
                           " to Control and ", 
                           floor(for.plot$N[1]*(for.plot$ARatio[1]/(1+for.plot$ARatio[1]))), 
                           " to Treatment")) +
    scale_x_continuous(expand=c(0,0))+
    scale_y_continuous(expand = c(0,0), breaks=seq(0,1,.2), limits=c(0,1), label=scales::percent)+
    theme(panel.spacing.x = unit(6, "mm"), axis.text.x = element_text(angle=45, hjust=1,vjust=1))+
    geom_vline(xintercept = c(for.plot$Delta.tv[1], for.plot$Delta.lrv[1]), linetype=2, color="blue")+
    annotate("text", label = paste0("Min TPP = ", for.plot$Delta.lrv[1]),
             x = for.plot$Delta.lrv[1], y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize, colour = "black", hjust = 1, vjust=0)+
    annotate("text", label = paste("Base TPP = ", for.plot$Delta.tv[1]),
             x = for.plot$Delta.tv[1], y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize, colour = "black", hjust = 0, vjust=0)
  
  table.plot2 <- ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ > Min TPP) > ", for.plot$tau.lrv[1]*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ > Base TPP) >  ", for.plot$tau.tv[1]*100,"$%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ > Min TPP) $\\leq$ ", for.plot$tau.ng[1]*100, "$%\\;\\;\\;\\,\\,$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ > Base TPP) $\\leq$ $", for.plot$tau.tv[1]*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+ 
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")
  
  # plot_grid(main.plot, 
  #            table.plot2, ncol=1, align = "v", rel_heights=c(.78,.22), axis="l")
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
}

#' @rdname TwoSampleNormalGamma
make.ts.ng.trt.oc2 <- function(for.plot=get.ts.ng.trt.oc.df(), nlines=25, tsize=4){
  my.plot2 <- ggplot(data=for.plot, aes(x=treatment.effect, y=freq, color=result))+
    geom_line(size=.75)+  scale_x_continuous(expand = c(0,0))+
    scale_y_continuous(expand = c(0,0), breaks=seq(0,1,.2), limits=c(0,1), label=scales::percent)+
    theme(panel.spacing.x = unit(6, "mm"), axis.text.x = element_text(angle=45, hjust=1,vjust=1))
  dpb <- ggplot_build(my.plot2)
  
  main.plot <- my.plot2+
    scale_color_manual(values=c("green", "grey50", "red"))+
    geom_vline(xintercept = c(for.plot$Delta.tv[1], for.plot$Delta.lrv[1]), color="blue", linetype=2)+
    annotate("text", label = paste0("Min TPP = ", for.plot$Delta.lrv[1]),
             x = for.plot$Delta.lrv[1], y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize, colour = "black", hjust = 1, vjust=0)+
    annotate("text", label = paste("Base TPP = ", for.plot$Delta.tv[1]),
             x = for.plot$Delta.tv[1], y = 0 + max(dpb$data[[1]]$y)/nlines, size = tsize, colour = "black", hjust = 0, vjust=0)+
    labs(x = "Underlying treatment effect",
         y = "Probability",
         title = "Operating characteristics as a function of treatment effect",
         subtitle = paste0("Total randomized with ", for.plot$N[1], " with ", for.plot$N[1] - floor(for.plot$N[1]*(for.plot$ARatio[1]/(1+for.plot$ARatio[1]))), " to Control and ", floor(for.plot$N[1]*(for.plot$ARatio[1]/(1+for.plot$ARatio[1]))), " to Treatment"),
         color = "Decision")+
    theme(legend.position = "bottom")
  
  table.plot2 <- ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ > Min TPP) > ", for.plot$tau.lrv[1]*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ > Base TPP) >  ", for.plot$tau.tv[1]*100,"$%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\,$ > Min TPP) $\\leq$ ", for.plot$tau.ng[1]*100, "$%\\;\\;\\;\\,\\,$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\,$ > Base TPP) $\\leq$ $", for.plot$tau.tv[1]*100,"$%")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+ 
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")
  
  # plot_grid(main.plot, 
  #            table.plot2, ncol=1, align = "v", rel_heights=c(.78,.22), axis="l")
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
  
}

#' @rdname TwoSampleNormalGamma
get.ts.ng.ssize.oc.df <- function(
  mu.0.c = 0, n.0.c = 1, alpha.0.c=.25, beta.0.c = 1, s.c = 2,  group.c="Control",
  mu.0.t = 3.75, n.0.t = 1, alpha.0.t=.25, beta.0.t = 1, s.t = 2,  group.t="Treatment",
  ARatio=2, N=50,
  n_LB_OC = floor(50*0.75), n_UB_OC=floor(50*2),
  Delta.tv = 2.5, Delta.lrv = 1.5, Delta.user = 4,
  tau.tv = 0.10, tau.lrv = 0.20, tau.ng = 0.35,
  npoints=10, n.MC = 500,
  seed = 1234){

  #cat("\nstarting decision.ts.continuous\n")
  ## We need to do this for a variety of sample sizes
  my.specs <- expand.grid(N=floor(seq(n_LB_OC, n_UB_OC, length.out=npoints))) %>%
    mutate(mu.0.c = mu.0.c,
           n.0.c = n.0.c,
           alpha.0.c=alpha.0.c,
           beta.0.c = beta.0.c,
           s.c = s.c,
           n.c = N - floor(N*(ARatio/(1+ARatio))),
           mu.0.t = mu.0.t,
           n.0.t = n.0.t,
           alpha.0.t=alpha.0.t,
           beta.0.t = beta.0.t,
           s.t = s.t,
           n.t = floor(N*(ARatio/(1+ARatio))),
           Delta.tv = Delta.tv,
           Delta.lrv = Delta.lrv,
           Delta.user=Delta.user,
           tau.tv = tau.tv,
           tau.lrv = tau.lrv,
           tau.ng = tau.ng,
           n.MC = n.MC,
           npoints = npoints,
           seed = seed + 0:(n()-1))

  # Now nrow(my.specs should be small, say npoints=10)
  my.FUN <- function(x){
    # Get x-bar samples pulled from the x'th row of my.specs
    # Reminder: Across rows of specs it is n.c and n.t that are changing!
    my.control.null <- rnorm(n = my.specs$n.MC[x], mean = my.specs$mu.0.c[x],
                             sd=my.specs$s.c[x]/sqrt(my.specs$n.c[x]))
    my.treatment.null <- rnorm(n = my.specs$n.MC[x], mean = my.specs$mu.0.c[x],
                               sd=my.specs$s.t[x]/sqrt(my.specs$n.t[x]))
    my.treatment.tv <- rnorm(n = my.specs$n.MC[x], mean = my.specs$mu.0.c[x] +
                               my.specs$Delta.tv[x], sd=my.specs$s.t[x]/sqrt(my.specs$n.t[x]))
    my.treatment.lrv <- rnorm(n = my.specs$n.MC[x], mean = my.specs$mu.0.c[x] +
                                my.specs$Delta.lrv[x], sd=my.specs$s.t[x]/sqrt(my.specs$n.t[x]))
    my.treatment.user <- rnorm(n = my.specs$n.MC[x], mean = my.specs$mu.0.c[x] +
                                 my.specs$Delta.user[x], sd=my.specs$s.t[x]/sqrt(my.specs$n.t[x]))

    # Merge with specs
    for.apply.null <- cbind(my.control.null, my.treatment.null, my.specs[x,])
    for.apply.tv <- cbind(my.control.null, my.treatment.tv, my.specs[x,])
    for.apply.lrv <- cbind(my.control.null, my.treatment.lrv, my.specs[x,])
    for.apply.user <- cbind(my.control.null, my.treatment.user, my.specs[x,])

    # Now evaluate decision outcomes based on the random sampling
    # These will be tallied later
    bind_rows(apply(X = matrix(1:nrow(for.apply.lrv)), MARGIN = 1, FUN = function(x) {
      for.return.null <- get.ts.ng.mc(mu.0.c = for.apply.null$mu.0.c[x],
                                      n.0.c = for.apply.null$n.0.c[x],
                                      alpha.0.c=for.apply.null$alpha.0.c[x],
                                      beta.0.c = for.apply.null$beta.0.c[x],
                                      xbar.c = for.apply.null$my.control.null[x],
                                      s.c = for.apply.null$s.c[x],
                                      n.c = for.apply.null$n.c[x], group.c="Control",
                                      mu.0.t = for.apply.null$mu.0.t[x],
                                      n.0.t = for.apply.null$n.0.t[x],
                                      alpha.0.t=for.apply.null$alpha.0.t[x],
                                      beta.0.t = for.apply.null$beta.0.t[x],
                                      xbar.t =  for.apply.null$my.treatment.null[x],
                                      s.t = for.apply.null$s.t[x],
                                      n.t = for.apply.null$n.t[x], group.t="Treatment",
                                      Delta.tv = for.apply.null$Delta.tv[x],
                                      Delta.lrv = for.apply.null$Delta.lrv[x],
                                      tau.tv = for.apply.null$tau.tv[x],
                                      tau.lrv = for.apply.null$tau.lrv[x],
                                      tau.ng = for.apply.null$tau.ng[x],
                                      n.MC = for.apply.null$n.MC[x]) %>%
        mutate(treatment.effect = "null",
               mu.0.c = for.apply.null$mu.0.c[x], n.0.c = for.apply.null$n.0.c[x],
               alpha.0.c=for.apply.null$alpha.0.c[x], beta.0.c = for.apply.null$beta.0.c[x],
               xbar.c = for.apply.null$my.control.null[x], s.c = for.apply.null$s.c[x],
               n.c = for.apply.null$n.c[x],
               group.c="Control",
               mu.0.t = for.apply.null$mu.0.t[x], n.0.t = for.apply.null$n.0.t[x],
               alpha.0.t=for.apply.null$alpha.0.t[x],
               beta.0.t = for.apply.null$beta.0.t[x],
               xbar.t =  for.apply.null$my.treatment.null[x], s.t = for.apply.null$s.t[x],
               n.t = for.apply.null$n.t[x],
               group.t="Treatment",
               Delta.tv = for.apply.null$Delta.tv[x], Delta.lrv = for.apply.null$Delta.lrv[x],
               tau.tv = for.apply.null$tau.tv[x], tau.lrv = for.apply.null$tau.lrv[x],
               tau.ng = for.apply.null$tau.ng[x],
               n.MC = for.apply.null$n.MC[x])

      for.return.tv <- get.ts.ng.mc(mu.0.c = for.apply.tv$mu.0.c[x],
                                    n.0.c = for.apply.tv$n.0.c[x],
                                    alpha.0.c=for.apply.tv$alpha.0.c[x],
                                    beta.0.c = for.apply.tv$beta.0.c[x],
                                    xbar.c = for.apply.tv$my.control.null[x],
                                    s.c = for.apply.tv$s.c[x], n.c = for.apply.tv$n.c[x],
                                    group.c="Control",
                                    mu.0.t = for.apply.tv$mu.0.t[x],
                                    n.0.t = for.apply.tv$n.0.t[x],
                                    alpha.0.t=for.apply.tv$alpha.0.t[x],
                                    beta.0.t = for.apply.tv$beta.0.t[x],
                                    xbar.t =  for.apply.tv$my.treatment.tv[x],
                                    s.t = for.apply.tv$s.t[x],
                                    n.t = for.apply.tv$n.t[x],
                                    group.t="Treatment",
                                    Delta.tv = for.apply.tv$Delta.tv[x],
                                    Delta.lrv = for.apply.tv$Delta.lrv[x],
                                    tau.tv = for.apply.tv$tau.tv[x],
                                    tau.lrv = for.apply.tv$tau.lrv[x],
                                    tau.ng = for.apply.tv$tau.ng[x],
                                    n.MC = for.apply.tv$n.MC[x]) %>%
        mutate(treatment.effect = "tv",
               mu.0.c = for.apply.tv$mu.0.c[x],
               n.0.c = for.apply.tv$n.0.c[x],
               alpha.0.c=for.apply.tv$alpha.0.c[x],
               beta.0.c = for.apply.tv$beta.0.c[x],
               xbar.c = for.apply.tv$my.control.null[x],
               s.c = for.apply.tv$s.c[x],
               n.c = for.apply.tv$n.c[x],
               group.c="Control",
               mu.0.t = for.apply.tv$mu.0.t[x],
               n.0.t = for.apply.tv$n.0.t[x],
               alpha.0.t=for.apply.tv$alpha.0.t[x],
               beta.0.t = for.apply.tv$beta.0.t[x],
               xbar.t =  for.apply.tv$my.treatment.tv[x],
               s.t = for.apply.tv$s.t[x],
               n.t = for.apply.tv$n.t[x],
               group.t="Treatment",
               Delta.tv = for.apply.tv$Delta.tv[x],
               Delta.lrv = for.apply.tv$Delta.lrv[x],
               tau.tv = for.apply.tv$tau.tv[x],
               tau.lrv = for.apply.tv$tau.lrv[x],
               tau.ng = for.apply.tv$tau.ng[x],
               n.MC = for.apply.tv$n.MC[x])
      for.return.lrv <- get.ts.ng.mc(mu.0.c = for.apply.lrv$mu.0.c[x],
                                     n.0.c = for.apply.lrv$n.0.c[x],
                                     alpha.0.c=for.apply.lrv$alpha.0.c[x],
                                     beta.0.c = for.apply.lrv$beta.0.c[x],
                                     xbar.c = for.apply.lrv$my.control.null[x],
                                     s.c = for.apply.lrv$s.c[x], n.c = for.apply.lrv$n.c[x],
                                     group.c="Control",
                                     mu.0.t = for.apply.lrv$mu.0.t[x],
                                     n.0.t = for.apply.lrv$n.0.t[x],
                                     alpha.0.t=for.apply.lrv$alpha.0.t[x],
                                     beta.0.t = for.apply.lrv$beta.0.t[x],
                                     xbar.t =  for.apply.lrv$my.treatment.lrv[x],
                                     s.t = for.apply.lrv$s.t[x],
                                     n.t = for.apply.lrv$n.t[x], group.t="Treatment",
                                     Delta.tv = for.apply.lrv$Delta.tv[x],
                                     Delta.lrv = for.apply.lrv$Delta.lrv[x],
                                     tau.tv = for.apply.lrv$tau.tv[x],
                                     tau.lrv = for.apply.lrv$tau.lrv[x],
                                     tau.ng = for.apply.lrv$tau.ng[x],
                                     n.MC = for.apply.lrv$n.MC[x]) %>%
        mutate(treatment.effect = "lrv",
               mu.0.c = for.apply.lrv$mu.0.c[x],
               n.0.c = for.apply.lrv$n.0.c[x],
               alpha.0.c=for.apply.lrv$alpha.0.c[x],
               beta.0.c = for.apply.lrv$beta.0.c[x],
               xbar.c = for.apply.lrv$my.control.null[x],
               s.c = for.apply.lrv$s.c[x],
               n.c = for.apply.lrv$n.c[x],
               group.c="Control",
               mu.0.t = for.apply.lrv$mu.0.t[x],
               n.0.t = for.apply.lrv$n.0.t[x],
               alpha.0.t=for.apply.lrv$alpha.0.t[x],
               beta.0.t = for.apply.lrv$beta.0.t[x],
               xbar.t =  for.apply.lrv$my.treatment.lrv[x],
               s.t = for.apply.lrv$s.t[x],
               n.t = for.apply.lrv$n.t[x],
               group.t="Treatment",
               Delta.tv = for.apply.lrv$Delta.tv[x],
               Delta.lrv = for.apply.lrv$Delta.lrv[x],
               tau.tv = for.apply.lrv$tau.tv[x],
               tau.lrv = for.apply.lrv$tau.lrv[x],
               tau.ng = for.apply.lrv$tau.ng[x],
               n.MC = for.apply.lrv$n.MC[x])

      for.return.user <- get.ts.ng.mc(mu.0.c = for.apply.user$mu.0.c[x],
                                      n.0.c = for.apply.user$n.0.c[x],
                                      alpha.0.c=for.apply.user$alpha.0.c[x],
                                      beta.0.c = for.apply.user$beta.0.c[x],
                                      xbar.c = for.apply.user$my.control.null[x],
                                      s.c = for.apply.user$s.c[x],
                                      n.c = for.apply.user$n.c[x],
                                      group.c="Control",
                                      mu.0.t = for.apply.user$mu.0.t[x],
                                      n.0.t = for.apply.user$n.0.t[x],
                                      alpha.0.t=for.apply.user$alpha.0.t[x],
                                      beta.0.t = for.apply.user$beta.0.t[x],
                                      xbar.t =  for.apply.user$my.treatment.user[x],
                                      s.t = for.apply.user$s.t[x],
                                      n.t = for.apply.user$n.t[x],
                                      group.t="Treatment",
                                      Delta.tv = for.apply.user$Delta.tv[x],
                                      Delta.lrv = for.apply.user$Delta.lrv[x],
                                      tau.tv = for.apply.user$tau.tv[x],
                                      tau.lrv = for.apply.user$tau.lrv[x],
                                      tau.ng = for.apply.user$tau.ng[x],
                                      n.MC = for.apply.user$n.MC[x]) %>%
        mutate(treatment.effect = "user",
               mu.0.c = for.apply.user$mu.0.c[x],
               n.0.c = for.apply.user$n.0.c[x],
               alpha.0.c=for.apply.user$alpha.0.c[x],
               beta.0.c = for.apply.user$beta.0.c[x],
               xbar.c = for.apply.user$my.control.null[x],
               s.c = for.apply.user$s.c[x],
               n.c = for.apply.user$n.c[x],
               group.c="Control",
               mu.0.t = for.apply.user$mu.0.t[x],
               n.0.t = for.apply.user$n.0.t[x],
               alpha.0.t=for.apply.user$alpha.0.t[x],
               beta.0.t = for.apply.user$beta.0.t[x],
               xbar.t =  for.apply.user$my.treatment.user[x],
               s.t = for.apply.user$s.t[x],
               n.t = for.apply.user$n.t[x],
               group.t="Treatment",
               Delta.tv = for.apply.user$Delta.tv[x],
               Delta.lrv = for.apply.user$Delta.lrv[x],
               tau.tv = for.apply.user$tau.tv[x],
               tau.lrv = for.apply.user$tau.lrv[x],
               tau.ng = for.apply.user$tau.ng[x],
               n.MC = for.apply.user$n.MC[x])

      for.return <- rbind(for.return.null, for.return.lrv, for.return.tv, for.return.user)
      #cat("\nget.NG.ts.df.mc finished\n")
      return(for.return)
    }))
  }

  # Now for each row of my.specs
  # Call my.FUN:
  get.results <- bind_rows(apply(X = matrix(1:nrow(my.specs)), MARGIN = 1, FUN = my.FUN))

  for.plot <- get.results %>% group_by(treatment.effect, mu.0.c, n.0.c, alpha.0.c,
                                       beta.0.c, s.c, n.c, group.c, mu.0.t, n.0.t,
                                       alpha.0.t, beta.0.t, s.t, n.t, group.t, Delta.tv,
                                       Delta.lrv, tau.tv, tau.lrv, tau.ng, n.MC) %>%
    summarize(NoGo=1 - sum(result=="No-Go")/my.specs$n.MC[1],
              Go = sum(result=="Go")/my.specs$n.MC[1]) %>% gather(value=value, key=result,
                                                                  Go, NoGo) %>%
    mutate(n.total=n.c+n.t)
  for.plot$treatment.effect <- factor(for.plot$treatment.effect,
                                      c("null", "lrv", "tv", "user"))
  levels(for.plot$treatment.effect) <-  c(TeX("$\\Delta\\,$ = NULL = 0%"),
                                          TeX(paste("$\\Delta\\,$ = Min TPP = $", Delta.lrv, "$ ")),
                                          TeX(paste("$\\Delta\\,$ = Base TPP = $", Delta.tv, "$ ")),
                                          TeX(paste("$\\Delta\\,$ = User defined = $", Delta.user, "$")))
  for.plot$treatment.effect <- factor(for.plot$treatment.effect,
                                      levels(for.plot$treatment.effect)[order(c(0, Delta.lrv, Delta.tv,
                                                                                Delta.user))])


  return(for.plot)
}

#' @rdname TwoSampleNormalGamma
make.ts.ng.ssize.oc <- function(for.plot=get.ts.ng.ssize.oc.df(), tsize=4, nlines=25){
  
  # Simply takes a dataframe from return.ssize.df and plots
  main.plot <-  for.plot %>% 
    ggplot(aes(x=n.total, y=value, color=result)) + 
    geom_line(alpha=.25, size=.75) + 
    facet_wrap(~treatment.effect, labeller="label_parsed")+
    scale_color_manual(values=c("green", "red"))+
    geom_ribbon(data=for.plot %>% dplyr::filter(result=="Go"), aes(x=n.total, ymin=0, ymax=value), fill=alpha("green",.5), color=alpha("green",.25))+
    geom_ribbon(data=for.plot %>% dplyr::filter(result=="Go"), aes(x=n.total, ymin=value, ymax=1), fill=alpha("grey",.5), color=alpha("grey",0))+
    geom_ribbon(data=for.plot %>% dplyr::filter(result=="NoGo"), aes(x=n.total, ymin=value, ymax=1), fill=alpha("red",.5), color=alpha("red", .25))+
    guides(color=F, fill=F)+
    labs(x="Total sample size", y="Decision Probabilities", 
         title=paste0("Operating Characteristics as a function of sample size"))+
    scale_x_continuous(expand=c(0,0))+
    scale_y_continuous(expand=c(0,0), breaks=seq(0,1,.2), minor_breaks=seq(0,1,.1), label=scales::percent)+
    theme(panel.spacing.x = unit(6, "mm"), axis.text.x = element_text(angle=45, hjust=1,vjust=1))
  
  table.plot2 <-  ggplot()+
    annotate("text", label = paste0("Decision Criteria"),
             x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
    annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
             x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\, \\geq$ Min TPP) > ", for.plot$tau.lrv[1]*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\, \\geq$ Base TPP) > ", for.plot$tau.tv[1]*100,"$%")),
             x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
    annotate("text", label = TeX(paste0("No-Go if:")), color="red",
             x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = TeX(paste0("P($\\Delta\\, \\geq$ Min TPP) $\\leq$", for.plot$tau.ng[1]*100, "$%\\;\\;\\;\\;\\,$ &")),
             x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
    annotate("text", label =  TeX(paste0("P($\\Delta\\, \\geq$ Base TPP) $\\leq$ ", for.plot$tau.tv[1]*100,"$%$")),
             x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+ 
    annotate("text", label = TeX(paste0("Consider:")), color="black",
             x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
    annotate("text", label = "Otherwise",
             x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
    scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
    scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
    labs(x="\n", y="")
  
  # plot_grid(main.plot, 
  #            table.plot2, ncol=1, align = "v", rel_heights=c(.78,.22), axis="l")
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
}

# TwoSampleTTE ------
#' @name TwoSampleTTE
#' @title Functions to Support the Two Sample Time to Event Scenario
#' @param m.con.prior prior number of control events
#' @param m.trt.prior prior number of treatment events
#' @param HR.prior prior estimate for HR
#' @param ARatio Allocation ratio
#' @param HR.obs observed hazard ratio
#' @param m.obs observed number of events
#' @param HR.lrv TPP Lower Reference Value aka Max TPP (large HRs lead to No-Go)
#' @param HR.tv TPP Target Value aka Base TPP
#' @param tau.tv threshold associated with Base TPP
#' @param tau.lrv threshold associated with Min TPP
#' @param tau.ng threshold associated with No-Go
#' @param seed random seed
#' @param nlines Control for text spacing
#' @param tsize Control for text size
#' @param nlines.ria Control for text spacing
#' @param for.plot data.frame returned by get.ts.ng.trt.oc.df
#' @examples
#' \dontrun{
#' make.ts.ng.ppp()
#' get.tte.post.param()
#' make.ts.ng.spp()
#' get.tte.trt.oc.df()
#' make.tte.trt.oc1()
#' make.tte.trt.oc2()
#' get.tte.ssize.oc.df()
#' make.tte.ssize.oc()
#' }
#' @description
#'
#'   make.tte.ppp: Make Two Sample Time to Event Prior/Posterior Plot. Returns a ggplot object.
#'
#'   make.tte.spp: Make Two Sample Time to Event Shaded Posterior Plot.  Returns a graphic built using grid.arrange.
#'
#'   get.tte.trt.oc.df: Get Two Sample Time to Event Treatment Effect OC. Returns a data.frame.
#'
#'   make.tte.trt.oc1: Make Two Sample Time to Event Treatment Effect. Returns a graphic built using grid.arrange.
#'
#'   make.tte.trt.oc2: Make Two Sample Time to Event Treatment Effect.  Returns a graphic built using grid.arrange.
#'
#'   get.tte.ssize.oc.df:  Get Two Sample Time to Event sample size OC data.frame. Returns a data.frame.
#'
#'   make.tte.ssize.oc: Make Two Sample Time to Event Sample size OC plot
make.tte.ppp <- function(m.con.prior = 10, m.trt.prior = 10, HR.prior = .8,
                         ARatio = 1, HR.obs = .75, m.obs = 50){
  ARatio <- 1/(1+ARatio)
  post.params = get.tte.post.param(m.con.prior = m.con.prior,
                                   m.trt.prior = m.trt.prior, HR.prior = HR.prior,
                                   ARatio = ARatio, HR.obs = HR.obs, m.obs = m.obs)
  pdfs <- rbind(
    gcurve(expr = dnorm(x, mean = log(HR.prior), sd = sqrt(4/(m.con.prior + m.trt.prior))),
           from = log(0.01), to = log(3), n = 1001, category = "Hazard ratio Prior") %>%
      mutate(HR = HR.prior, events = m.con.prior + m.trt.prior, group="Prior"),
    gcurve(expr = dnorm(x, mean = post.params[1,1], sd = post.params[1,2]),
           from = log(0.01), to = log(3), n = 1001, category="Hazard ratio posterior") %>%
      mutate(HR = post.params[1,1], events = m.con.prior + m.trt.prior + m.obs[2],
             group="Posterior")
  ) %>% mutate(group=factor(group, c("Prior", "Posterior")))
  levels(pdfs$group) <- c(paste0("Prior for HR: ", HR.prior, " worth ",
                                 m.trt.prior+ m.con.prior, " events"),
                          paste0("Posterior for HR: ", round(exp(post.params[1]),4),
                                 " worth ", m.trt.prior+ m.con.prior + m.obs, " events"))

  ggplot(data = pdfs, aes(x = exp(x),y = y, color = category))+geom_line(size=.75) +
    facet_wrap(~group)+guides(color = F)+
    labs(x = "Hazard Ratio", y=NULL,
         title="Prior and posterior distribtuions for the hazard ratio",
         subtitle = paste0("Data: Observed hazard ratio: ", round(HR.obs,2)," based on ",
                           m.obs, " events."),
         color="Density") +
    theme(legend.position = "bottom")+
    scale_x_continuous(breaks = seq(0,3,.2))+
    scale_y_continuous(breaks=NULL, labels = NULL)
}

#' @rdname TwoSampleTTE
make.tte.spp <- function(m.con.prior = 10,m.trt.prior = 10, HR.prior=1.1, ARatio=1, HR.obs=.89, m.obs = 50,
                         HR.tv= 1.3, HR.lrv = 1.1, tau.tv=.1, tau.lrv=.8, tau.ng=.65, tsize = 4, nlines = 25, nlines.ria=20){
  # Run this for fixed values of 
  results <- get.tte.df(m.con.prior = m.con.prior, m.trt.prior = m.trt.prior, 
                        HR.prior = HR.prior, ARatio = ARatio, m.obs=m.obs, 
                        HR.obs=seq(0.01, 1.9,.01), HR.tv=HR.tv, HR.lrv = HR.lrv, 
                        tau.tv=tau.tv, tau.lrv=tau.lrv, tau.ng=tau.ng)
  
  # Determine min number of TRT responders for Go
  if(HR.tv < HR.lrv){
    result.go <- results %>% dplyr::filter(result== "Go") %>% slice(n())
    # Determine max number of TRT responders for No-Go
    result.ng <- results %>% dplyr::filter(result== "No-Go") %>% slice(1)
  } else {
    result.go <- results %>% dplyr::filter(result== "Go") %>% slice(1)
    # Determine max number of TRT responders for No-Go
    result.ng <- results %>% dplyr::filter(result== "No-Go") %>% slice(n())
  }
  
  post.params <- get.tte.post.param(m.con.prior = m.con.prior,m.trt.prior = m.trt.prior, HR.prior = HR.prior, ARatio = ARatio, HR.obs = HR.obs, m.obs = m.obs)
  
  my.df <- rbind(
    gcurve(expr = dnorm(x, mean = post.params[1,1],sd = post.params[1,2]), from = log(0.01), to = log(3),
           n = 1001, category=paste0("Posterior hazard ratio: ", round(exp(post.params[1]),2), " worth ", m.trt.prior+ m.con.prior + m.obs, " events")) %>%
      mutate(HR = post.params[1,1], events = m.con.prior + m.trt.prior + m.obs, group="Posterior"))
  
  if(HR.tv < HR.lrv){
    P.R1 = pnorm(log(HR.tv), post.params[1,1], post.params[1,2])
    P.R3 = pnorm(log(HR.lrv),post.params[1,1], post.params[1,2])
    
  } else {
    P.R1 = 1 - pnorm(log(HR.tv), post.params[1,1], post.params[1,2])
    P.R3 = 1 - pnorm(log(HR.lrv),post.params[1,1], post.params[1,2])
  }
  
  # Annotation line 1: Decision ----
  
  result = ifelse(P.R1 > tau.tv  & P.R3 > tau.lrv, "Go",
                  ifelse(P.R1 < tau.tv  & P.R3 < tau.ng, "No-Go", "Consider"))
  for.decision <- paste0("Decision: ", result) 
  result.color <- ifelse(result=="Go", "darkgreen", ifelse(result=="No-Go", "red", "black"))
  
  
  # Annotation line 2: Decision interval ----
  if(HR.tv < HR.lrv){
    if(result!="No-Go"){
      a <- exp(qnorm(p = tau.tv,mean = post.params[1,1], sd = post.params[1,2]))
      b <- exp(qnorm(p = tau.ng,mean = post.params[1,1], sd = post.params[1,2]))
      for.decision.interval <- paste0("Decision interval: (",  
                                      round(a,3), ", ", 
                                      round(b,3),")")
    } else {
      a <- exp(qnorm(p = tau.tv,mean = post.params[1,1], sd = post.params[1,2]))
      b <- exp(qnorm(p = tau.lrv,mean = post.params[1,1], sd = post.params[1,2]))
      for.decision.interval <- paste0("Decision interval: (",  
                                      round(a,3), ", ", 
                                      round(b,3),")")
    }
  } else {
    if(result!="No-Go"){
      a <- exp(qnorm(p = 1-tau.ng,mean = post.params[1,1], sd = post.params[1,2]))      
      b <- exp(qnorm(p = 1-tau.tv,mean = post.params[1,1], sd = post.params[1,2]))
      for.decision.interval <- paste0("Decision interval: (", 
                                      round(a,3), ", ", 
                                      round(b,3),")")
    } else {
      a <- exp(qnorm(p = 1- tau.lrv,mean = post.params[1,1], sd = post.params[1,2]))
      b <- exp(qnorm(p = 1 - tau.tv,mean = post.params[1,1], sd = post.params[1,2]))
      for.decision.interval <- paste0("Decision interval: (",  
                                      round(a,3), ", ", 
                                      round(b,3),")")
    }
    
  }
  for.subtitle <- paste0("Data ", "(Obs HR, Total events): (", round(HR.obs, 2), ", ", m.obs, "). Observed HR needed for Go: ",
                         round(result.go$HR.obs,2), ". Needed for No-Go: ", round(result.ng$HR.obs,2))
  
  # Annotation line 3: P(Delta < HR Max)
  
  if(HR.tv < HR.lrv){
    annotate.P1 <- ifelse(result=="Go", 
                          TeX(paste0("P($\\Delta$ $\\leq$ $\\HR_{Max}$) = $", round(P.R3,2), "$ >", tau.lrv)), 
                          ifelse(result=="No-Go", 
                                 TeX(paste0("P($\\Delta$ $\\leq$ $\\HR_{Max}$) = $", round(P.R3,2), "$ < ", tau.ng)),
                                 TeX(paste0("P($\\Delta$ $\\leq$ $\\HR_{Max}$) = $", round(P.R3,2),"$"  ))))
    # Annotation line 4: P(Delta < HR Base)
    
    annotate.P2 <- ifelse(result=="Go", 
                          TeX(paste0("P($\\Delta$ $\\leq$ $\\HR_{Base}$) = $",  round(P.R1 ,2), "$ > ", tau.tv)),
                          ifelse(result=="No-Go", TeX(paste0("P($\\Delta$ $\\leq$ $\\HR_{Base}$) = $",  round(P.R1,2), "$ < ", tau.tv)),
                                 TeX(paste0("P($\\Delta$ $\\leq$ $\\HR_{Base}$) = $",  round(P.R1,2),"$"))))
  } else {
    annotate.P1 <- ifelse(result=="Go", 
                          TeX(paste0("P($\\Delta$ $\\geq$ $\\HR_{Max}$) = $", round(P.R3,2), "$ >", tau.lrv)), 
                          ifelse(result=="No-Go", 
                                 TeX(paste0("P($\\Delta$ $\\geq$ $\\HR_{Max}$) = $", round(P.R3,2), "$ < ", tau.ng)),
                                 TeX(paste0("P($\\Delta$ $\\geq$ $\\HR_{Max}$) = $", round(P.R3,2),"$"  ))))
    # Annotation line 4: P(Delta < HR Base)
    
    annotate.P2 <- ifelse(result=="Go", 
                          TeX(paste0("P($\\Delta$ $\\geq$ $\\HR_{Base}$) = $",  round(P.R1 ,2), "$ > ", tau.tv)),
                          ifelse(result=="No-Go", TeX(paste0("P($\\Delta$ $\\geq$ $\\HR_{Base}$) = $",  round(P.R1,2), "$ < ", tau.tv)),
                                 TeX(paste0("P($\\Delta$ $\\geq$ $\\HR_{Base}$) = $",  round(P.R1,2),"$")))) 
    
  }
  # Initialize a ggplot
  
  dplot <- ggplot(data = my.df, aes(x = exp(x),y = y)) + geom_line() + facet_wrap(~category)
  
  # Access the ggplot to get goodies to help accomplish shading
  dpb <- ggplot_build(dplot)
  x1.1 <- max(which(dpb$data[[1]]$x <= HR.lrv))+1
  x2.1 <- max(which((dpb$data[[1]]$x) <= (HR.tv)))
  x1.2 <- max(which((dpb$data[[1]]$x) <= (HR.tv)))
  x2.2 <- min(which((dpb$data[[1]]$x) >= 0))
  
  if(which(dpb$data[[1]]$y== max(dpb$data[[1]]$y)) < length(dpb$data[[1]]$x)/2){
    annotate.x = min(min(dpb$data[[1]]$x), HR.tv, 0)
    annotate.j = 0
  } else {annotate.x = max(max(dpb$data[[1]]$x), HR.lrv, 3)
  annotate.j = 1}
  
  
  for.decision.interval.df <- data.frame(x = round(a,3), 
                                         xend = round(b,3),
                                         y= min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3, 
                                         yend = min(dpb$data[[1]]$y) + max(dpb$data[[1]]$y)/nlines.ria * 3,
                                         group=my.df$group[1])
  
  
  # Introduce shading
  main.plot <- dplot +
    geom_area(data = data.frame(x = dpb$data[[1]]$x[x1.1:length(dpb$data[[1]]$x)],
                                y = dpb$data[[1]]$y[x1.1:length(dpb$data[[1]]$y)]),
              aes(x=x, y = y), fill=alpha("red", 0))+
    geom_area(data = data.frame(x = dpb$data[[1]]$x[x1.1:x2.1],
                                y = dpb$data[[1]]$y[x1.1:x2.1]),
              aes(x=(x), y = y), fill=alpha("grey80", 0))+
    geom_area(data = data.frame(x = dpb$data[[1]]$x[x1.2:x2.2],
                                y = dpb$data[[1]]$y[x1.2:x2.2]),
              aes(x=(x), y = y), fill=alpha("green", 0))+
    geom_line(data = my.df, aes(x = exp(x),y = y), size=.75)
  
  
  main.plot <- main.plot + 
    labs(title="Posterior distribution for the hazard ratio",
         subtitle = for.subtitle,
         x="Hazard Ratio", y = NULL,
         caption = NULL)+
    annotate("text", label = for.decision,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 0, size = tsize+1, colour = result.color, hjust = annotate.j)+
    annotate("text", label = for.decision.interval,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 1, size = tsize, colour = result.color, hjust = annotate.j)+
    annotate("text", label = annotate.P1, color=result.color,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 2, size = tsize,  hjust = annotate.j)+
    annotate("text", label = annotate.P2, color=result.color,
             x = annotate.x, y = max(dpb$data[[1]]$y)-max(dpb$data[[1]]$y)/nlines.ria * 3, size = tsize,  hjust = annotate.j) +
    scale_y_continuous(breaks=NULL, labels=NULL)+
    scale_x_continuous(limits = c(0,3.01),
                       breaks = pretty(x=c(0,3), n=10))
  # Add reference lines and Credible interval
  
  if(result.color == "red"){
    main.plot <- main.plot + 
      geom_segment(data=for.decision.interval.df, aes(x=x, xend=xend, y=y, yend=yend, group=group), arrow = arrow(ends="both", angle = 90), color=result.color, size=.75)
  } else {
    main.plot <- main.plot + 
      geom_segment(data=for.decision.interval.df, aes(x=x, xend=xend, y=y, yend=yend, group=group),  arrow = arrow(ends="both", angle = 90), color=result.color, size=.75)
  }
  
  if(HR.tv < HR.lrv){
  main.plot <- main.plot +
    geom_vline(xintercept = c(HR.lrv, HR.tv), linetype = 2, color = c("blue", "blue"))+
    annotate("text", label =paste0("Max TPP = ", HR.lrv), x = HR.lrv, y = 0 + max(dpb$data[[1]]$y)/nlines.ria, size = tsize, colour = "black", hjust = 0)+
    annotate("text", label = paste0("Base TPP = ",  HR.tv), x = HR.tv, y = 0 + max(dpb$data[[1]]$y)/nlines.ria, size = tsize, colour = "black", hjust = 1)
  } else {
    main.plot <- main.plot +
      geom_vline(xintercept = c(HR.lrv, HR.tv), linetype = 2, color = c("blue", "blue"))+
      annotate("text", label =paste0("Max TPP = ", HR.lrv), x = HR.lrv, y = 0 + max(dpb$data[[1]]$y)/nlines.ria, size = tsize, colour = "black", hjust = 1)+
      annotate("text", label = paste0("Base TPP = ",  HR.tv), x = HR.tv, y = 0 + max(dpb$data[[1]]$y)/nlines.ria, size = tsize, colour = "black", hjust = 0)
  }
  if(HR.tv < HR.lrv){
    table.plot2 <- ggplot()+
      annotate("text", label = paste0("Decision Criteria"),
               x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
      annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
               x = -1, y = 1-2.5/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = TeX(paste0("P($\\Delta$ < $HR_{Max}$) > ", tau.lrv*100, "%$\\;\\;\\;\\;\\,$ &")),
               x = 0, y = 1-2.5/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
      annotate("text", label =  TeX(paste0("P($\\Delta$ < $HR_{Base}$) > $", tau.tv*100,"$%")),
               x = 1.5, y = 1-2.5/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
      annotate("text", label = TeX(paste0("No-Go if:")), color="red",
               x = -1, y = 1-3.25/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = TeX(paste0("P($\\Delta$ < $HR_{Max}$) $\\leq$ ", tau.ng*100, "%$\\;\\;\\;\\,\\,$ &")),
               x = 0, y = 1-3.25/nlines, size = tsize+1, colour = "red", hjust = 0) +
      annotate("text", label =  TeX(paste0("P($\\Delta$ < $HR_{Base}$) $\\leq$ ", tau.tv*100,"$%")),
               x = 1.5, y = 1-3.25/nlines, size = tsize+1, colour = "red", hjust = 0)+ 
      annotate("text", label = TeX(paste0("Consider:")), color="black",
               x = -1, y = 1-4./nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = "Otherwise",
               x = 0, y = 1-4./nlines, size = tsize+1, colour = "black", hjust = 0)+
      scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
      scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
      labs(x="", y="")} else {
        table.plot2 <- ggplot()+
          annotate("text", label = paste0("Decision Criteria"),
                   x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
          annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
                   x = -1, y = 1-2.5/nlines, size = tsize+1, hjust = 0)+
          annotate("text", label = TeX(paste0("P($\\Delta$ > $HR_{Max}$) > ", tau.lrv*100, "%$\\;\\;\\;\\;\\,$ &")),
                   x = 0, y = 1-2.5/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
          annotate("text", label =  TeX(paste0("P($\\Delta$ > $HR_{Base}$) > $", tau.tv*100,"$%")),
                   x = 1.5, y = 1-2.5/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
          annotate("text", label = TeX(paste0("No-Go if:")), color="red",
                   x = -1, y = 1-3.25/nlines, size = tsize+1, hjust = 0)+
          annotate("text", label = TeX(paste0("P($\\Delta$ > $HR_{Max}$) $\\leq$ ", tau.ng*100, "%$\\;\\;\\;\\,\\,$ &")),
                   x = 0, y = 1-3.25/nlines, size = tsize+1, colour = "red", hjust = 0) +
          annotate("text", label =  TeX(paste0("P($\\Delta$ > $HR_{Base}$) $\\leq$ ", tau.tv*100,"$%")),
                   x = 1.5, y = 1-3.25/nlines, size = tsize+1, colour = "red", hjust = 0)+ 
          annotate("text", label = TeX(paste0("Consider:")), color="black",
                   x = -1, y = 1-4./nlines, size = tsize+1, hjust = 0)+
          annotate("text", label = "Otherwise",
                   x = 0, y = 1-4./nlines, size = tsize+1, colour = "black", hjust = 0)+
          scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
          scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
          labs(x="", y="")
      }
  
  
  
  # plot_grid(main.plot, 
  #            table.plot2, ncol=1, align = "v", rel_heights=c(.78,.22), axis="l")
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
  
}

#' @rdname TwoSampleTTE
get.tte.post.param <- function(m.con.prior=10, m.trt.prior=10, HR.prior=.7,
                               ARatio=1, HR.obs=.8, m.obs=50){
  ARatio <- 1/(1+ARatio)
  data.frame(
    post.mean = (1 / m.con.prior + 1 / m.trt.prior) ^ (-1)/
      ((1 / m.con.prior + 1 / m.trt.prior) ^ (-1) +
         (1/(m.obs * ARatio*(1 - ARatio))) ^ (-1))*log(HR.prior) +
      (1/(m.obs * ARatio*(1 - ARatio))) ^ (-1)/((1 / m.con.prior + 1 / m.trt.prior) ^ (-1) +
                                                  (1/(m.obs * ARatio*(1 - ARatio))) ^ (-1))*log(HR.obs),
    post.sd = sqrt(((1 / m.con.prior + 1 / m.trt.prior)*
                      (1/(m.obs * ARatio*(1 - ARatio))))/((1 / m.con.prior + 1 / m.trt.prior)+
                                                            (1/(m.obs * ARatio*(1 - ARatio))))))
}
#' @rdname TwoSampleTTE
get.tte.df <- function(m.con.prior = 50, m.trt.prior = 50, HR.prior = .845,
                       HR.obs = seq(0.3, 1, .01), m.obs = seq(10, 200, 5),
                       ARatio = 0.5,
                       HR.tv = .80, HR.lrv = 0.9,
                       tau.tv = .10, tau.lrv = .20, tau.ng = 0.35) {
  
  #if this holds we assume that smaller HR's are preferred
  P = ARatio / (ARatio + 1)
  my.grid <- expand.grid(m.con.prior = m.con.prior, m.trt.prior = m.trt.prior,
                         HR.prior = HR.prior, ARatio = ARatio, HR.obs = HR.obs, m.obs = m.obs)  
  my.results <- plyr:::list_to_dataframe(do.call(Map, c(f = get.tte.post.param, my.grid)))
  if(HR.tv < HR.lrv){
    my.grid <- cbind(my.grid, my.results) %>%
      mutate(m.prior = m.con.prior + m.trt.prior,
             ARatio=ARatio,
             HR.tv = HR.tv,
             HR.lrv = HR.lrv,
             tau.tv = tau.tv,
             tau.lrv = tau.lrv,
             tau.ng = tau.ng) %>% 
      mutate(
        P.R1 = pnorm(log(HR.tv), post.mean, post.sd),  # P(HR < HR.tv)  larger values here preferred
        P.R3 = pnorm(log(HR.lrv),post.mean, post.sd),  # P(HR < HR.tv)
        result = factor(
          ifelse(P.R1 > tau.tv & P.R3 > tau.lrv, "Go", 
                 ifelse(P.R1 < tau.tv & P.R3 < tau.ng, "No-Go", "Consider")),
          c("Go", "Consider", "No-Go")))
  } else 
  {
    my.grid <- cbind(my.grid, my.results) %>%
      mutate(m.prior = m.con.prior + m.trt.prior,
             ARatio=ARatio,
             HR.tv = HR.tv,
             HR.lrv = HR.lrv,
             tau.tv = tau.tv,
             tau.lrv = tau.lrv,
             tau.ng = tau.ng) %>%
      mutate(
        P.R1 = 1 - pnorm(log(HR.tv), post.mean, post.sd), # P(HR < HR.tv)   Should be big for a go
        P.R3 = 1 - pnorm(log(HR.lrv),post.mean, post.sd), # P(HR < HR.lrv)  Should be big for a go
        result = factor(
          ifelse(P.R1 > tau.tv & P.R3 > tau.lrv, "Go",
                 ifelse(P.R1 < tau.tv & P.R3 < tau.ng, "No-Go", "Consider")),
          c("Go", "Consider", "No-Go")))
  }
  
  return(my.grid)
}

#' @rdname TwoSampleTTE
get.tte.trt.oc.df <- function(m.con.prior = 10, m.trt.prior = 10, HR.prior = .8,
                              ARatio = .5, m.obs = 50,
                              HR.tv = 1.3, HR.lrv = 1.1,
                              HR.lower=0.3, HR.upper=2,
                              tau.tv = 0.1, tau.lrv = 0.8, tau.ng = 0.65){
  
  if(HR.tv < HR.lrv) {
    # Get results for fine grid of HR.obs values
    results <- get.tte.df(m.con.prior = m.con.prior, m.trt.prior = m.trt.prior,
                          HR.prior = HR.prior,
                          HR.obs = seq(HR.lower, HR.upper, .005), m.obs = m.obs,
                          ARatio = ARatio,
                          HR.tv = HR.tv, HR.lrv = HR.lrv,
                          tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng)
    
    # Identify min/max needed for no-go, go
    result.go <- results %>% dplyr::filter(result== "Go") %>% slice(n())
    # Determine max number of TRT responders for No-Go
    result.ng <- results %>% dplyr::filter(result== "No-Go") %>% slice(1)
    
    my.df <- results %>%
      mutate(Go=pnorm(q =log(result.go$HR.obs), mean=log(HR.obs), sd=sqrt(4/m.obs)),
             NoGo=1 - pnorm(q = log(result.ng$HR.obs), mean = log(HR.obs), sqrt(4/m.obs))) %>%
      mutate(m.con.prior = m.con.prior, m.trt.prior = m.trt.prior, HR.prior = HR.prior,
             m.obs = m.obs, ARatio = ARatio,
             HR.tv = HR.tv, HR.lrv = HR.lrv,
             tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng,
             HR.upper=HR.upper, HR.lower=HR.lower)} else {
               results <- get.tte.df(m.con.prior = m.con.prior, m.trt.prior = m.trt.prior,
                                     HR.prior = HR.prior,
                                     HR.obs = seq(HR.lower, HR.upper, .005), m.obs = m.obs,
                                     ARatio = ARatio,
                                     HR.tv = HR.tv, HR.lrv = HR.lrv,
                                     tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng)
               
               # Identify min/max needed for no-go, go
               result.go <- results %>% dplyr::filter(result== "Go") %>% slice(1)
               # Determine max number of TRT responders for No-Go
               result.ng <- results %>% dplyr::filter(result== "No-Go") %>% slice(n())
               
               my.df <- results %>%
                 mutate(Go=1 - pnorm(q =log(result.go$HR.obs), mean=log(HR.obs), sd=sqrt(4/m.obs)),
                        NoGo= pnorm(q = log(result.ng$HR.obs), mean = log(HR.obs), sqrt(4/m.obs))) %>%
                 mutate(m.con.prior = m.con.prior, m.trt.prior = m.trt.prior, HR.prior = HR.prior,
                        m.obs = m.obs, ARatio = ARatio,
                        HR.tv = HR.tv, HR.lrv = HR.lrv,
                        tau.tv = tau.tv, tau.lrv = tau.lrv, tau.ng = tau.ng,
                        HR.upper=HR.upper, HR.lower=HR.lower)        
               
             }
  
  
  my.df
}

#' @rdname TwoSampleTTE
make.tte.trt.oc1 <- function(plot.df = get.tte.trt.oc.df(), nlines=25, tsize=4){
  HR.tv = plot.df$HR.tv[1]
  HR.lrv = plot.df$HR.lrv[1]
  tau.tv = plot.df$tau.tv[1]
  tau.lrv = plot.df$tau.lrv[1]
  tau.ng = plot.df$tau.ng[1]
  m.obs = plot.df$m.obs[1]
  ARatio = plot.df$ARatio[1]
  HR.upper = plot.df$HR.upper[1]
  HR.lower = plot.df$HR.lower[1]
  
  if(HR.tv < HR.lrv) {
    main.plot <- ggplot() + geom_line(data=plot.df, aes(x=HR.obs, y=Go),
                                      color="lightgreen") +
      geom_line(data=plot.df, aes(x=HR.obs, y=1 - NoGo), color="red")
    dpb <- ggplot_build(main.plot)
    
    main.plot <- main.plot + geom_line(data=plot.df, aes(x=HR.obs, y=Go),
                                       color="lightgreen") +
      geom_line(data=plot.df, aes(x=HR.obs, y=1 - NoGo), color="red") +
      geom_ribbon(data=plot.df, aes(x=HR.obs, ymin=0, ymax=Go),
                  fill="lightgreen", alpha=.5)+
      geom_ribbon(data=plot.df, aes(x=HR.obs, ymin=Go, ymax=1-NoGo),
                  fill="grey", alpha=.5)+
      geom_ribbon(data=plot.df, aes(x=HR.obs, ymin=1-NoGo, ymax=1),
                  fill="red", alpha=.5)+
      labs(title="Operating characteristics as a function of treatment effect",
           subtitle=paste0("Total events: ", m.obs,
                           ". Randomization ratio (Control:Treatment): (1:", ARatio, ")"))+
      annotate("text", label = TeX(paste0(HR.lrv,"$=\\HR_{Max}$")), x = HR.lrv,
               y = 0 + 2*max(dpb$data[[1]]$y)/nlines, size =
                 tsize, colour = "black", hjust = 0)+
      annotate("text", label = TeX(paste( HR.tv, "$\\HR_{Base}=$")), x = HR.tv,
               y = 0 + 2*max(dpb$data[[1]]$y)/nlines, size = tsize,
               colour = "black", hjust = 1)+
      labs(x = "Underlying Hazard Ratio",
           y="Probability")+
      geom_vline(xintercept=(c(HR.tv, HR.lrv)), color="blue",linetype=2)+
      scale_x_continuous(expand = c(0,0), breaks=seq(0,2,.2),
                         limits=c(HR.lower, HR.upper))+
      scale_y_continuous(expand=c(0,0), breaks=seq(0,1,.2), minor_breaks=seq(0,1,.1),
                         label=scales::percent)+
      theme(panel.spacing.x = unit(6, "mm"), axis.text.x =
              element_text(angle=45, hjust=1,vjust=1))
    
    table.plot2 <- ggplot()+
      annotate("text", label = paste0("Decision Criteria"),
               x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
      annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
               x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = TeX(paste0("P($\\Delta\\,$", " ", "$\\leq$", " ", "$HR_{Max}$) > ", plot.df$tau.lrv[1]*100, "%$\\;\\;\\;\\;\\,$ &")),
               x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
      annotate("text", label =  TeX(paste0("P($\\Delta\\,$", " ", "$\\leq$", " ", "$HR_{Base}$) > $", plot.df$tau.tv[1]*100,"$%")),
               x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
      annotate("text", label = TeX(paste0("No-Go if:")), color="red",
               x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = TeX(paste0("P($\\Delta\\,$", " ", "$\\leq$", " ", "$HR_{Max}$) $\\leq$ ", plot.df$tau.ng[1]*100, "%$\\;\\;\\;\\,$ &")),
               x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
      annotate("text", label =  TeX(paste0("P($\\Delta\\,$", " ", "$\\leq$", " ", "$HR_{Base}$) $\\leq$ ", plot.df$tau.tv[1]*100,"$%")),
               x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+ 
      annotate("text", label = TeX(paste0("Consider:")), color="black",
               x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = "Otherwise",
               x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
      scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
      scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
      labs(x="\n", y="")
    grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
  } else {
    main.plot <- ggplot() + geom_line(data=plot.df, aes(x=HR.obs, y=Go),
                                      color="lightgreen") +
      geom_line(data=plot.df, aes(x=HR.obs, y=1 - NoGo), color="red")
    dpb <- ggplot_build(main.plot)
    
    main.plot <- main.plot + geom_line(data=plot.df, aes(x=HR.obs, y=Go),
                                       color="lightgreen") +
      geom_line(data=plot.df, aes(x=HR.obs, y=1 - NoGo), color="red") +
      geom_ribbon(data=plot.df, aes(x=HR.obs, ymin=0, ymax=Go),
                  fill="lightgreen", alpha=.5)+
      geom_ribbon(data=plot.df, aes(x=HR.obs, ymin=Go, ymax=1-NoGo),
                  fill="grey", alpha=.5)+
      geom_ribbon(data=plot.df, aes(x=HR.obs, ymin=1-NoGo, ymax=1),
                  fill="red", alpha=.5)+
      labs(title="Operating characteristics as a function of treatment effect",
           subtitle=paste0("Total events: ", m.obs,
                           ". Randomization ratio (Control:Treatment): (1:", ARatio, ")"))+
      annotate("text", label = TeX(paste0("$\\HR_{Max}$=",HR.lrv)), x = HR.lrv,
               y = 0 + 2*max(dpb$data[[1]]$y)/nlines, size =
                 tsize, colour = "black", hjust = 1)+
      annotate("text", label = TeX(paste( "$\\HR_{Base}=$",HR.tv)), x = HR.tv,
               y = 0 + 2*max(dpb$data[[1]]$y)/nlines, size = tsize,
               colour = "black", hjust = 0)+
      labs(x = "Underlying Hazard Ratio",
           y="Probability")+
      geom_vline(xintercept=(c(HR.tv, HR.lrv)), color="blue",linetype=2)+
      scale_x_continuous(expand = c(0,0), breaks=seq(0,2,.2),
                         limits=c(HR.lower, HR.upper))+
      scale_y_continuous(expand=c(0,0), breaks=seq(0,1,.2), minor_breaks=seq(0,1,.1),
                         label=scales::percent)+
      theme(panel.spacing.x = unit(6, "mm"), axis.text.x =
              element_text(angle=45, hjust=1,vjust=1))
    
    table.plot2 <- ggplot()+
      annotate("text", label = paste0("Decision Criteria"),
               x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
      annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
               x = -1, y = 1-2.5/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = TeX(paste0("P($\\Delta$ > $HR_{Max}$) > ", tau.lrv*100, "%$\\;\\;\\;\\;\\,$ &")),
               x = 0, y = 1-2.5/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
      annotate("text", label =  TeX(paste0("P($\\Delta$ > $HR_{Base}$) > $", tau.tv*100,"$%")),
               x = 1.5, y = 1-2.5/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
      annotate("text", label = TeX(paste0("No-Go if:")), color="red",
               x = -1, y = 1-3.25/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = TeX(paste0("P($\\Delta$ > $HR_{Max}$) $\\leq$ ", tau.ng*100, "%$\\;\\;\\;\\,\\,$ &")),
               x = 0, y = 1-3.25/nlines, size = tsize+1, colour = "red", hjust = 0) +
      annotate("text", label =  TeX(paste0("P($\\Delta$ > $HR_{Base}$) $\\leq$ ", tau.tv*100,"$%")),
               x = 1.5, y = 1-3.25/nlines, size = tsize+1, colour = "red", hjust = 0)+ 
      annotate("text", label = TeX(paste0("Consider:")), color="black",
               x = -1, y = 1-4./nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = "Otherwise",
               x = 0, y = 1-4./nlines, size = tsize+1, colour = "black", hjust = 0)+
      scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
      scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
      labs(x="", y="")
  }
  
  
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
  
}
#' @rdname TwoSampleTTE
make.tte.trt.oc2 <- function(plot.df = get.tte.trt.oc.df(), nlines=25, tsize=4){
  HR.tv = plot.df$HR.tv[1]
  HR.lrv = plot.df$HR.lrv[1]
  tau.tv = plot.df$tau.tv[1]
  tau.lrv = plot.df$tau.lrv[1]
  tau.ng = plot.df$tau.ng[1]
  m.obs = plot.df$m.obs[1]
  ARatio = plot.df$ARatio[1]
  HR.upper = plot.df$HR.upper[1]
  HR.lower = plot.df$HR.lower[1]
  
  
  
  if(HR.tv < HR.lrv) {
    main.plot <- plot.df %>% mutate(Consider = 1 - Go - NoGo) %>%
      gather(key = key, value=value, c(Go, NoGo, Consider), factor_key = T) %>%
      dplyr::filter(key !="NoGo2") %>%
      ggplot(aes(x= HR.obs, y=value, color=key))+
      geom_line(size=.75)+
      scale_color_manual(values = c("Go" = "green", "NoGo" = "red",
                                    "Consider" = "darkgrey"))
    dpb <- ggplot_build(main.plot)
    main.plot <- main.plot +
      geom_vline(xintercept = c(HR.tv, HR.lrv), color="blue", linetype=2)+
      labs(x="Underlying Hazard Ratio", y="Probability",
           color="Decision",
           title="Operating characteristics as a function of treatment effect",
           subtitle=paste0("Total events: ", m.obs,
                           ". Randomization ratio (Control:Treatment): (1:",
                           ARatio, ")"))+
      theme(legend.position = "bottom")+
      annotate("text", label = TeX(paste0("$HR_{Max}$ = ",HR.lrv)), x = HR.lrv,
               y = 0 + 2*max(dpb$data[[1]]$y)/nlines, size =
                 tsize, colour = "black", hjust = 0)+
      annotate("text", label = TeX(paste0("$\\HR_{Base}$ = ", HR.tv)), x = HR.tv,
               y = 0 + 2*max(dpb$data[[1]]$y)/nlines, size = tsize,
               colour = "black", hjust = 1)+
      scale_x_continuous(expand = c(0,0), breaks=seq(0,2,.2),
                         limits=c(HR.lower, HR.upper))+
      scale_y_continuous(expand = c(0,0), breaks=seq(0,1,.2),
                         minor_breaks=seq(0,1,.1), limits=c(0,1), label=scales::percent)+
      theme(panel.spacing.x = unit(6, "mm"), axis.text.x = element_text(angle=45, hjust=1,vjust=1))
    
    table.plot2 <- ggplot()+
      annotate("text", label = paste0("Decision Criteria"),
               x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
      annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
               x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = TeX(paste0("P($\\Delta\\,$", " ", "$\\leq$", " ", "$HR_{Max}$) > ", plot.df$tau.lrv[1]*100, "%$\\;\\;\\;\\;\\,$ &")),
               x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
      annotate("text", label =  TeX(paste0("P($\\Delta\\,$", " ", "$\\leq$", " ", "$HR_{Base}$) > $", plot.df$tau.tv[1]*100,"$%")),
               x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
      annotate("text", label = TeX(paste0("No-Go if:")), color="red",
               x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = TeX(paste0("P($\\Delta\\,$", " ", "$\\leq$", " ", "$HR_{Max}$) $\\leq$ ", plot.df$tau.ng[1]*100, "%$\\;\\;\\;\\,$ &")),
               x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
      annotate("text", label =  TeX(paste0("P($\\Delta\\,$", " ", "$\\leq$", " ", "$HR_{Base}$) $\\leq$ ", plot.df$tau.tv[1]*100,"$%")),
               x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+ 
      annotate("text", label = TeX(paste0("Consider:")), color="black",
               x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = "Otherwise",
               x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
      scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
      scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
      labs(x="\n", y="")
  } else {   
    main.plot <- plot.df %>% mutate(Consider = 1 - Go - NoGo) %>%
      gather(key = key, value=value, c(Go, NoGo, Consider), factor_key = T) %>%
      dplyr::filter(key !="NoGo2") %>%
      ggplot(aes(x= HR.obs, y=value, color=key))+
      geom_line(size=.75)+
      scale_color_manual(values = c("Go" = "green", "NoGo" = "red",
                                    "Consider" = "darkgrey"))
    dpb <- ggplot_build(main.plot)
    main.plot <- main.plot +
      geom_vline(xintercept = c(HR.tv, HR.lrv), color="blue", linetype=2)+
      labs(x="Underlying Hazard Ratio", y="Probability",
           color="Decision",
           title="Operating characteristics as a function of treatment effect",
           subtitle=paste0("Total events: ", m.obs,
                           ". Randomization ratio (Control:Treatment): (1:",
                           ARatio, ")"))+
      theme(legend.position = "bottom")+
      annotate("text", label = TeX(paste0("$HR_{Max}$ = ",HR.lrv)), x = HR.lrv,
               y = 0 + 2*max(dpb$data[[1]]$y)/nlines, size =
                 tsize, colour = "black", hjust = 1)+
      annotate("text", label = TeX(paste0("$\\HR_{Base}$ = ", HR.tv)), x = HR.tv,
               y = 0 + 2*max(dpb$data[[1]]$y)/nlines, size = tsize,
               colour = "black", hjust = 0)+
      scale_x_continuous(expand = c(0,0), breaks=seq(0,2,.2),
                         limits=c(HR.lower, HR.upper))+
      scale_y_continuous(expand = c(0,0), breaks=seq(0,1,.2),
                         minor_breaks=seq(0,1,.1), limits=c(0,1), label=scales::percent)+
      theme(panel.spacing.x = unit(6, "mm"), axis.text.x = element_text(angle=45, hjust=1,vjust=1))
    
    table.plot2 <- ggplot()+
      annotate("text", label = paste0("Decision Criteria"),
               x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
      annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
               x = -1, y = 1-2.5/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = TeX(paste0("P($\\Delta$ > $HR_{Max}$) > ", tau.lrv*100, "%$\\;\\;\\;\\;\\,$ &")),
               x = 0, y = 1-2.5/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
      annotate("text", label =  TeX(paste0("P($\\Delta$ > $HR_{Base}$) > $", tau.tv*100,"$%")),
               x = 1.5, y = 1-2.5/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
      annotate("text", label = TeX(paste0("No-Go if:")), color="red",
               x = -1, y = 1-3.25/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = TeX(paste0("P($\\Delta$ > $HR_{Max}$) $\\leq$ ", tau.ng*100, "%$\\;\\;\\;\\,\\,$ &")),
               x = 0, y = 1-3.25/nlines, size = tsize+1, colour = "red", hjust = 0) +
      annotate("text", label =  TeX(paste0("P($\\Delta$ > $HR_{Base}$) $\\leq$ ", tau.tv*100,"$%")),
               x = 1.5, y = 1-3.25/nlines, size = tsize+1, colour = "red", hjust = 0)+ 
      annotate("text", label = TeX(paste0("Consider:")), color="black",
               x = -1, y = 1-4./nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = "Otherwise",
               x = 0, y = 1-4./nlines, size = tsize+1, colour = "black", hjust = 0)+
      scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
      scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
      labs(x="", y="")
    
  }
  
  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
}

#' @rdname TwoSampleTTE
get.tte.ssize.oc.df <- function(m.con.prior=10, m.trt.prior=10, HR.prior=.75,
                                ARatio=1, m.obs=50,
                                m.lower=40, m.upper=120,
                                HR.lrv=1.1, HR.tv=1.4, HR.user = 0.845, tau.tv=.1,
                                tau.lrv=.8, tau.ng=.65){
  
  specs <- expand.grid(m.obs=floor(seq(m.lower, m.upper, length.out=20)),
                       HRs.from=c(HR.lrv, HR.tv, HR.user, 1)) %>%
    mutate(m.con.prior=m.con.prior, m.trt.prior=m.trt.prior, HR.prior=HR.prior,
           ARatio=ARatio, HR.lrv=HR.lrv, HR.tv=HR.tv, HR.user = HR.user, tau.tv=tau.tv,
           tau.lrv=tau.lrv, tau.ng=tau.ng)
  
  
  if(HR.tv < HR.lrv){
    for.plot <- bind_rows(apply(X = matrix(1:nrow(specs)), MARGIN = 1,
                                FUN = function(x){
                                  # For each row of specs, get all possible probs along a fine grid
                                  results <- get.tte.df(m.con.prior = specs$m.con.prior[x], m.trt.prior = specs$m.trt.prior[x],
                                                        HR.prior = specs$HR.prior[x],
                                                        ARatio = specs$ARatio[x], m.obs=specs$m.obs[x],
                                                        HR.obs=seq(0.1, 2,.005), HR.tv=specs$HR.tv[x],
                                                        HR.lrv = specs$HR.lrv[x], tau.tv=specs$tau.tv[x],
                                                        tau.lrv=specs$tau.lrv[x], tau.ng=specs$tau.ng[x])
                                  # Identify max criteria for Go and min criteria for no-go
                                  result.go <- results %>% dplyr::filter(result== "Go") %>% slice(n())
                                  # Determine max number of TRT responders for No-Go
                                  result.ng <- results %>% dplyr::filter(result== "No-Go") %>% slice(1)
                                  
                                  my.df <- data.frame(Go = pnorm(q = log(result.go$HR.obs),
                                                                 mean = log(c(HR.lrv, HR.tv, HR.user, 1)),
                                                                 sd=sqrt(4/result.go$m.obs)),
                                                      NoGo = 1 - pnorm(q = log(result.ng$HR.obs),
                                                                       mean = log(c(HR.lrv, HR.tv, HR.user, 1)),
                                                                       sd=sqrt(4/result.go$m.obs)),
                                                      HR=c(HR.lrv, HR.tv, HR.user, 1))  %>%
                                    mutate(m.con.prior = specs$m.con.prior[x],
                                           m.trt.prior = specs$m.trt.prior[x], HR.prior = specs$HR.prior[x],
                                           ARatio = specs$ARatio[x], m.obs=specs$m.obs[x],
                                           HR.tv=specs$HR.tv[x], HR.lrv = specs$HR.lrv[x], HR.user=specs$HR.user[x],
                                           tau.tv=specs$tau.tv[x], tau.lrv=specs$tau.lrv[x], tau.ng=specs$tau.ng[x])
                                  my.df
                                }))
    for.plot$key.label = c("HR = Max TPP = ", "HR = Base TPP = ", "HR = User's TPP = ", "HR = Null = ")
    for.plot$key = factor(paste0(for.plot$key.label, for.plot$HR),
                          paste0(for.plot$key.label, for.plot$HR)[1:4][order(c(HR.lrv, HR.tv, HR.user, 1))])
    
  } else{
    for.plot <- bind_rows(apply(X = matrix(1:nrow(specs)), MARGIN = 1,
                                FUN = function(x){
                                  # For each row of specs, get all possible probs along a fine grid
                                  results <- get.tte.df(m.con.prior = specs$m.con.prior[x], m.trt.prior = specs$m.trt.prior[x],
                                                        HR.prior = specs$HR.prior[x],
                                                        ARatio = specs$ARatio[x], m.obs=specs$m.obs[x],
                                                        HR.obs=seq(0.1, 2,.005), HR.tv=specs$HR.tv[x],
                                                        HR.lrv = specs$HR.lrv[x], tau.tv=specs$tau.tv[x],
                                                        tau.lrv=specs$tau.lrv[x], tau.ng=specs$tau.ng[x])
                                  # Identify max criteria for Go and min criteria for no-go
                                  result.go <- results %>% dplyr::filter(result== "Go") %>% slice(1)
                                  # Determine max number of TRT responders for No-Go
                                  result.ng <- results %>% dplyr::filter(result== "No-Go") %>% slice(n())
                                  
                                  my.df <- data.frame(Go = 1 - pnorm(q = log(result.go$HR.obs),
                                                                     mean = log(c(HR.lrv, HR.tv, HR.user, 1)),
                                                                     sd=sqrt(4/result.go$m.obs)),
                                                      NoGo = pnorm(q = log(result.ng$HR.obs),
                                                                   mean = log(c(HR.lrv, HR.tv, HR.user, 1)),
                                                                   sd=sqrt(4/result.go$m.obs)),
                                                      HR=c(HR.lrv, HR.tv, HR.user, 1))  %>%
                                    mutate(m.con.prior = specs$m.con.prior[x],
                                           m.trt.prior = specs$m.trt.prior[x], HR.prior = specs$HR.prior[x],
                                           ARatio = specs$ARatio[x], m.obs=specs$m.obs[x],
                                           HR.tv=specs$HR.tv[x], HR.lrv = specs$HR.lrv[x], HR.user=specs$HR.user[x],
                                           tau.tv=specs$tau.tv[x], tau.lrv=specs$tau.lrv[x], tau.ng=specs$tau.ng[x])
                                  my.df
                                }))
    for.plot$key.label = c("HR = Max TPP = ", "HR = Base TPP = ", "HR = User's TPP = ", "HR = Null = ")
    for.plot$key = factor(paste0(for.plot$key.label, for.plot$HR),
                          paste0(for.plot$key.label, for.plot$HR)[1:4][order(c(HR.lrv, HR.tv, HR.user, 1))])
    
  }
  
  
  for.plot
}

#' @rdname TwoSampleTTE
make.tte.ssize.oc <- function(for.plot=get.tte.ssize.oc.df(HR.lrv=1.1, HR.tv=1.4), tsize=4, nlines=25){
  
  # Simply takes a dataframe from get.tte.ssize.oc.df and plots
  main.plot <- ggplot() + geom_line(data=for.plot, aes(x=m.obs, y=Go), color="green") + 
    geom_line(data=for.plot, aes(x=m.obs, y=1-NoGo), color="red") +
    facet_wrap(~key)+
    geom_ribbon(data=for.plot, aes(x=m.obs, ymin=0, ymax=Go), fill="green", alpha=.5)+
    geom_ribbon(data=for.plot, aes(x=m.obs, ymin=1-NoGo, ymax=1), fill="red", alpha=.5)+
    geom_ribbon(data=for.plot, aes(x=m.obs, ymin=Go, ymax=1-NoGo), fill="grey", alpha=.5)+
    
    labs(x="Total number of observed events", y="Decision Probabilities", 
         title="Operating Characteristics as a function of number of events")+
    scale_x_continuous(expand=c(0,0))+
    scale_y_continuous(expand=c(0,0), breaks=seq(0,1,.2), minor_breaks=seq(0,1,.1), label=scales::percent)+
    theme(panel.spacing.x = unit(6, "mm"), axis.text.x = element_text(angle=45, hjust=1,vjust=1))
  
  if(for.plot$HR.tv[1] < for.plot$HR.lrv[1]){
    table.plot2 <- ggplot()+
      annotate("text", label = paste0("Decision Criteria"),
               x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
      annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
               x = -1, y = 1-3/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = TeX(paste0("P($\\Delta\\,$", " ", "$\\leq$", " ", "$HR_{Max}$) > ", for.plot$tau.lrv[1]*100, "%$\\;\\;\\;\\;\\,$ &")),
               x = 0, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
      annotate("text", label =  TeX(paste0("P($\\Delta\\,$", " ", "$\\leq$", " ", "$HR_{Base}$) > $", for.plot$tau.tv[1]*100,"$%")),
               x = 1.5, y = 1-3/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
      annotate("text", label = TeX(paste0("No-Go if:")), color="red",
               x = -1, y = 1-4/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = TeX(paste0("P($\\Delta\\,$", " ", "$\\leq$", " ", "$HR_{Max}$) $\\leq$ ", for.plot$tau.ng[1]*100, "%$\\;\\;\\;\\,$ &")),
               x = 0, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0) +
      annotate("text", label =  TeX(paste0("P($\\Delta\\,$", " ", "$\\leq$", " ", "$HR_{Base}$) $\\leq$ ", for.plot$tau.tv[1]*100,"$%")),
               x = 1.5, y = 1-4/nlines, size = tsize+1, colour = "red", hjust = 0)+ 
      annotate("text", label = TeX(paste0("Consider:")), color="black",
               x = -1, y = 1-5/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = "Otherwise",
               x = 0, y = 1-5/nlines, size = tsize+1, colour = "black", hjust = 0)+
      scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
      scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
      labs(x="\n", y="")
  } else {
    table.plot2 <- ggplot()+
      annotate("text", label = paste0("Decision Criteria"),
               x = -1, y = .95, size = tsize+2, colour = "black", hjust=0)+
      annotate("text", label = TeX(paste0("Go if:")), color="darkgreen",
               x = -1, y = 1-2.5/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = TeX(paste0("P($\\Delta$ > $HR_{Max}$) > ", for.plot$tau.lrv[1]*100, "%$\\;\\;\\;\\;\\,$ &")),
               x = 0, y = 1-2.5/nlines, size = tsize+1, colour = "darkgreen", hjust = 0) +
      annotate("text", label =  TeX(paste0("P($\\Delta$ > $HR_{Base}$) > $", for.plot$tau.tv[1]*100,"$%")),
               x = 1.5, y = 1-2.5/nlines, size = tsize+1, colour = "darkgreen", hjust = 0)+
      annotate("text", label = TeX(paste0("No-Go if:")), color="red",
               x = -1, y = 1-3.25/nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = TeX(paste0("P($\\Delta$ > $HR_{Max}$) $\\leq$ ", for.plot$tau.ng[1]*100, "%$\\;\\;\\;\\,\\,$ &")),
               x = 0, y = 1-3.25/nlines, size = tsize+1, colour = "red", hjust = 0) +
      annotate("text", label =  TeX(paste0("P($\\Delta$ > $HR_{Base}$) $\\leq$ ", for.plot$tau.tv[1]*100,"$%")),
               x = 1.5, y = 1-3.25/nlines, size = tsize+1, colour = "red", hjust = 0)+ 
      annotate("text", label = TeX(paste0("Consider:")), color="black",
               x = -1, y = 1-4./nlines, size = tsize+1, hjust = 0)+
      annotate("text", label = "Otherwise",
               x = 0, y = 1-4./nlines, size = tsize+1, colour = "black", hjust = 0)+
      scale_y_continuous(limits=c(0.75,.975), expand=c(0,0), labels=NULL, breaks=NULL)+
      scale_x_continuous(limits=c(-1,3.5), expand=c(0,0), labels=NULL, breaks=NULL)+
      labs(x="", y="")
  }

  grid.arrange(main.plot, table.plot2, nrow=2, heights=c(.78,.22))
  
}