R/conduct.R

In DTEBOP2: Bayesian Optimal Phase II Randomized Clinical Trial Design with Delayed Outcomes

#' Computes the hazard ratio between post- and pre-separation periods, provides the 95% Bayesian credible interval,
#' and makes a Go/No-Go decision based on the trial data. This function is intended for interim and final analyses
#' in two-arm survival trials with delayed treatment effects.
#'
#' Computes the hazard ratio between post- and pre-separation hazard rates, along with the 95% Bayesian credible interval. Based on the user-supplied survival data, the function also provides a Go/No-Go decision for interim or final analysis in two-arm survival trials with delayed treatment effects
#' @importFrom stats pbeta rexp runif
#' @importFrom truncdist rtrunc
#' @param data_E A data frame containing the survival data for the experimental (treatment) arm. The first column should be the observed survival times, and the second column should be the event indicator (1 = event, 0 = censored).
#' @param data_C A data frame containing the survival data for the control arm. The format is the same as data_E: first column for observed time, second column for status (1 = event, 0 = censored).
#' @param median.1 Numeric. The overall median survival time for the standard-of-care (SOC) arm.
#' @param median.2 Numeric. The overall median survival time for the experimental arm.
#' @param gprior.E_1 Optional. A numeric vector of length two specifying the shape and scale parameters of the inverse-gamma prior for the pre-separation mean survival time (i.e., 1/hazard rate). If NULL, the default is  \code{c(4,3/log(2)*median.1)}
#' @param gprior.E_2 Optional. A numeric vector of length two specifying the shape and scale parameters of the inverse-gamma prior for the post-separation mean survival time (i.e., 1/hazard rate). If NULL, the default is  \code{c(4,6/log(2)*median.1)}
#' @param L Numeric. The lower bound of the delayed treatment effect (DTE) separation timepoint.
#' @param U Numeric. The upper bound of the delayed treatment effect (DTE) separation timepoint.
#' @param trunc.para Numeric vector of length two. Specifies the shape and scale parameters of the truncated Gamma prior for the DTE timepoint.
#' @param lambda Numeric. Weighting parameter used in the posterior decision function:\eqn{1-\lambda(n/N)^{\gamma}}
#' @param gamma   Numeric. Exponent parameter used in the posterior decision function:\eqn{1-\lambda(n/N)^{\gamma}}
#' @param nmax Integer. Maximum sample size per arm.
#' @param S_likely Numeric. The most likely value of the DTE separation timepoint. Defaults to the midpoint of \code{L} and \code{U}
#' @param value If TRUE, returns posterior means, 95% credible intervals, and event rates. Default is FALSE.
#' @return A list with the following components:
#' \item{posterior_mean}{Posterior mean of the hazard ratio (post- / pre-separation).}
#' \item{credible_interval}{95% Bayesian credible interval for the hazard ratio.}
#' \item{event_rates}{Event rates for the experimental arm, control arm, and overall.}
#' \item{decision}{
#'   One of the following trial decisions based on observed data:
#'   \describe{
#'     \item{Go}{Trial proceeds to next stage.}
#'     \item{No-Go}{Trial stops early for futility.}
#'     \item{Reject Null Hypothesis}{Final analysis leads to rejecting H0.}
#'     \item{Not Reject Null Hypothesis}{Final analysis fails to reject H0.}
#'   }
#' }
#' @examples
#' \donttest{
#' conduct(data_E, data_C, median.1 = 2.8, median.2 = 3.5, L = 2, U = 2.5,
#'         trunc.para = c(1, 1), lambda = 0.95, gamma = 1, nmax = 63, S_likely = 2.3)
#'         }
#' @export
conduct <- function(data_E, data_C, median.1, median.2, gprior.E_1 = NULL, gprior.E_2 = NULL,
                    L, U, trunc.para = c(1, 1), lambda, gamma, nmax, S_likely=(L+U)/2,value=FALSE) {

  median_inuse <- function(median_0, median_1, S) {
    if (median_0 < S) median_1 <- median_0 else median_1 <- (median_1 - S) / (1 - S / median_0)
    return(c(median_0, median_1))
  }

  # Median calculation
  median_vals <- median_inuse(median.1, median.2, S = S_likely)
  median_1 <- median_vals[1]
  median_2 <- median_vals[2]

  # Prior setting
  gprior.E_1 <- if (is.null(gprior.E_1)) c(4, 3 * median_1 / log(2)) else gprior.E_1
  gprior.E_2 <- if (is.null(gprior.E_2)) c(4, 6 * median_1 / log(2)) else gprior.E_2

  current.n <- dim(data_E)[1]
  time_data <- data.frame(Time = data_E[, 1], Status = data_E[, 2])
  TTOT <- TTOT_fun(time_data, t = S_likely)

  para1 <- gprior.E_1[1] + TTOT$event_count_1 + sum(data_C[, 2])
  para2 <- gprior.E_2[1] + TTOT$event_count_2

  cutoff <- (gprior.E_1[2] + TTOT$TTOT_1 + sum(data_C[, 1])) /
    (gprior.E_1[2] + gprior.E_2[2] + TTOT$TTOT_1 + TTOT$TTOT_2 + sum(data_C[, 1]))

  temp <- (pbeta(cutoff, para1, para2) > 1 - lambda * (current.n / nmax)^gamma)

  a_0 <- gprior.E_1[1] + TTOT$event_count_1 + sum(data_C[, 2])
  a_1 <- gprior.E_2[1] + TTOT$event_count_2
  b_0 <- log(2) * (gprior.E_1[2] + TTOT$TTOT_1 + sum(data_C[, 1]))
  b_1 <- log(2) * (gprior.E_2[2] + TTOT$TTOT_2)

  posterior <- rinvgamma(10000, a_0, b_0) / rinvgamma(10000, a_1, b_1)
  ratio_value <- mean(posterior)
  ci_value <- quantile(posterior, probs = c(0.025, 0.975))

  death_all <- sum(data_C[, 2]) + sum(data_E[, 2])
  r_all <- death_all / (nrow(data_C) + nrow(data_E))
  r_1 <- sum(data_E[, 2]) / nrow(data_E)
  r_0 <- sum(data_C[, 2]) / nrow(data_C)

  # Decision logic
  decision <- if (temp == 1 && current.n == nmax) {
    "Not Reject Null hypothesis"
  } else if (temp == 0 && current.n == nmax) {
    "Reject Null hypothesis"
  } else if (temp == 1) {
    "No-Go"
  } else {
    "Go"
  }
  if (!value) {
    result <- list(
      sprintf("Posterior mean of ratio \u0303\u03bc0 to \u0303\u03bc1: %.6f", ratio_value),
      sprintf("95%% credible interval of the ratio:     %.4f    %.4f", ci_value[1], ci_value[2]),
      sprintf("Experimental/SOC/Overall event rate: %.4f/%.4f/%.4f",
              r_1, r_0, r_all),
      sprintf("Decision: %s", decision)
    )
    return(result)
  }
  if(value){
    result_1 <- list(
      ratio=ratio_value,
      CI= c(ci_value[1], ci_value[2]),
      event_rate=c(r_1, r_0, r_all)
    )
    return(result_1)
  }

}