R/post_summary_bin_2arm.R

In SAMprior: Self-Adapting Mixture (SAM) Priors

## internal cache for pairwise beta-beta posterior probabilities
.bin_postprob_cache <- new.env(parent = emptyenv())

.beta_mean_vec <- function(a, b) a / (a + b)

.beta_var_vec <- function(a, b) a * b / ((a + b)^2 * (a + b + 1))

## Computes P(theta_t - theta > margin | data)
.post_prob_pair_bin <- function(a.t, b.t, a.c, b.c, margin, rel.tol) {
  key <- paste(a.t, b.t, a.c, b.c, format(margin, digits = 16),
               format(rel.tol, digits = 16), sep = "|")

  if (exists(key, envir = .bin_postprob_cache, inherits = FALSE)) {
    return(get(key, envir = .bin_postprob_cache, inherits = FALSE))
  }

  lower <- max(0, margin)
  upper <- 1

  val <- if (lower >= upper) {
    0
  } else {
    integrand <- function(pt) {
      pbeta(pt - margin, shape1 = a.c, shape2 = b.c) *
        dbeta(pt, shape1 = a.t, shape2 = b.t)
    }

    integrate(
      f = integrand,
      lower = lower,
      upper = upper,
      rel.tol = rel.tol,
      subdivisions = 200L
    )$value
  }

  assign(key, val, envir = .bin_postprob_cache)
  val
}

## Computes P(theta_t - theta > margin | data) from posterior mixtures
.post_prob_bin_from_posts <- function(post_t, post_c, margin, rel.tol) {
  out <- 0

  for (j in seq_along(post_t$w)) {
    wt_j <- post_t$w[j]
    a_tj <- post_t$a[j]
    b_tj <- post_t$b[j]

    for (k in seq_along(post_c$w)) {
      out <- out + wt_j * post_c$w[k] *
        .post_prob_pair_bin(
          a.t = a_tj,
          b.t = b_tj,
          a.c = post_c$a[k],
          b.c = post_c$b[k],
          margin = margin,
          rel.tol = rel.tol
        )
    }
  }

  out
}

#' Posterior Summary for Two-Arm Comparative Trial with Binary Endpoint
#'
#' The \code{post_summary_bin_2arm} function is designed to compute the
#' posterior summary for the treatment effect in a two-arm comparative trial
#' with a binary endpoint under one of three borrowing strategies:
#' self-adapting mixture prior (SAM), robust MAP prior with fixed weight
#' (rMAP), or non-informative prior (NP).
#'
#' The treatment effect is defined as \eqn{\tau = \theta_t - \theta}, where
#' \eqn{\theta_t} and \eqn{\theta} denote the response rates in the treatment
#' and control arms, respectively. Inference is based on the posterior
#' distribution of \eqn{\tau} given the observed response counts from the two
#' arms.
#'
#' @param x.t Observed number of responses in the treatment arm.
#' @param x Observed number of responses in the control arm.
#' @param if.prior Informative prior constructed based on historical data for
#' the control arm, represented (approximately) as a beta mixture prior.
#' @param nf.prior Non-informative prior used as the robustifying component
#' for the control arm prior.
#' @param prior.t Prior used for the treatment arm. If missing, the default
#' value is set to be \code{nf.prior}.
#' @param n.t Sample size for the treatment arm.
#' @param n Sample size for the control arm.
#' @param delta Clinically significant difference used for the SAM prior.
#' This argument is only used when \code{method = "SAM"}.
#' @param cutoff Posterior probability cutoff used for decision making.
#' The null hypothesis is rejected if the posterior tail probability
#' exceeds \code{cutoff}.
#' @param method Borrowing strategy for the control arm. Must be one of
#' \code{"SAM"}, \code{"rMAP"}, or \code{"NP"}.
#' @param alternative Direction of the posterior decision. Must be one of
#' \code{"greater"} (for superiority) or \code{"less"} (for inferiority).
#' @param margin Clinical margin. Must be a non-negative scalar. The default
#' value is \code{0}.
#' @param weight_rMAP Weight assigned to the informative prior component
#' (\eqn{0 \leq} \code{weight_rMAP} \eqn{\leq 1}) for the robust MAP prior.
#' This argument is only used when \code{method = "rMAP"}. The default value is
#' 0.5.
#' @param method.w Methods used to determine the mixture weight for SAM priors.
#' The default method is "LRT" (Likelihood Ratio Test), the alternative option
#' is "PPR" (Posterior Probability Ratio). See \code{\link{SAM_weight}} for
#' more details.
#' @param prior.odds The prior probability of \eqn{H_0} being true compared to
#' the prior probability of \eqn{H_1} being true using PPR method. The default
#' value is 1. See \code{\link{SAM_weight}} for more details.
#' @param rel.tol Relative tolerance for numerical integration used to evaluate
#' the posterior probability.
#'
#' @details
#' The posterior for the treatment arm is obtained by updating
#' \code{prior.t} using the observed response count \code{x.t}. The posterior
#' for the control arm depends on the selected borrowing strategy:
#'
#' \itemize{
#' \item \code{SAM}: the prior for the control arm is a mixture of
#' \code{if.prior} and \code{nf.prior}, with adaptive mixture weight
#' determined by the SAM borrowing rule.
#' \item \code{rMAP}: the prior for the control arm is a fixed mixture of
#' \code{if.prior} and \code{nf.prior}, with fixed weight
#' \code{weight_rMAP}.
#' \item \code{NP}: the prior for control arm is \code{nf.prior}.
#' }
#'
#' When \code{alternative = "greater"}, inference is based on
#' \eqn{P(\theta_t - \theta > margin \mid X_t, X)}. When
#' \code{alternative = "less"}, inference is based on
#' \eqn{P(\theta_t - \theta < -margin \mid X_t, X)}.
#'
#' The posterior mean and posterior variance of \eqn{\tau} are defined as
#' \deqn{E(\tau \mid X_t, X) = E(\theta_t \mid X_t) - E(\theta \mid X),}
#' and
#' \deqn{\mathrm{Var}(\tau \mid X_t, X) = \mathrm{Var}(\theta_t \mid X_t) + \mathrm{Var}(\theta \mid X),}
#' where independence between the treatment and control arm posteriors is
#' assumed conditional on the current and historical data.
#'
#' @return A list containing the posterior probability in the requested
#' direction, posterior mean and variance of \eqn{\tau},
#' decision indicator, borrowing weight used for the control arm prior,
#' and the corresponding trial data and method information.
#'
#' @seealso \code{\link{SAM_weight}}, \code{\link{SAM_prior}}
#'
#' @export
post_summary_bin_2arm <- function(x.t, x,
                                  if.prior, nf.prior, prior.t = nf.prior,
                                  n.t, n,
                                  delta,
                                  cutoff,
                                  method = c("SAM", "rMAP", "NP"),
                                  alternative = c("greater", "less"),
                                  margin = 0,
                                  weight_rMAP = 0.5,
                                  method.w = "LRT",
                                  prior.odds = 1,
                                  rel.tol = 1e-8) {
  method <- match.arg(method)
  alternative <- match.arg(alternative)

  if (!is.numeric(x.t) || length(x.t) != 1 || !is.finite(x.t) || x.t < 0 || x.t > n.t) {
    stop("`x.t` must be a scalar between 0 and `n.t`.")
  }
  if (!is.numeric(x) || length(x) != 1 || !is.finite(x) || x < 0 || x > n) {
    stop("`x` must be a scalar between 0 and `n`.")
  }
  if (!is.numeric(cutoff) || length(cutoff) != 1 || cutoff <= 0 || cutoff >= 1) {
    stop("`cutoff` must be a scalar in (0, 1).")
  }
  if (!is.numeric(margin) || length(margin) != 1 || !is.finite(margin) || margin < 0) {
    stop("`margin` must be a non-negative scalar.")
  }
  if (!is.numeric(n.t) || length(n.t) != 1 || !is.finite(n.t) || n.t <= 0) {
    stop("`n.t` must be a positive scalar.")
  }
  if (!is.numeric(n) || length(n) != 1 || !is.finite(n) || n <= 0) {
    stop("`n` must be a positive scalar.")
  }
  if (!method.w %in% c("LRT", "PPR")) {
    stop("`method.w` must be either \"LRT\" or \"PPR\".")
  }
  if (!is.numeric(prior.odds) || length(prior.odds) != 1 ||
      !is.finite(prior.odds) || prior.odds <= 0) {
    stop("`prior.odds` must be a positive scalar.")
  }
  if (!is.numeric(rel.tol) || length(rel.tol) != 1 ||
      !is.finite(rel.tol) || rel.tol <= 0) {
    stop("`rel.tol` must be a positive scalar.")
  }

  post_t <- .posterior_mixbeta_update(
    prior = prior.t,
    x = x.t,
    n = n.t
  )

  if (method == "SAM") {
    post_c <- .control_sam_update_bin(
      if.prior = if.prior,
      nf.prior = nf.prior,
      x = x,
      n = n,
      delta = delta,
      weight_fun = weight_fun_betamix,
      method.w = method.w,
      prior.odds = prior.odds
    )
    weight_used <- post_c$wSAM

  } else if (method == "rMAP") {
    if (is.null(weight_rMAP) || length(weight_rMAP) != 1 ||
        !is.finite(weight_rMAP) || weight_rMAP < 0 || weight_rMAP > 1) {
      stop("`weight_rMAP` must be a scalar in [0, 1] when `method = \"rMAP\"`.")
    }

    post_c <- .control_fixed_update_bin(
      if.prior = if.prior,
      nf.prior = nf.prior,
      x = x,
      n = n,
      weight = weight_rMAP
    )
    weight_used <- weight_rMAP

  } else {
    post_c <- .posterior_mixbeta_update(
      prior = nf.prior,
      x = x,
      n = n
    )
    weight_used <- 0
  }

  post_prob_greater <- .post_prob_bin_from_posts(
    post_t = post_t,
    post_c = post_c,
    margin = margin,
    rel.tol = rel.tol
  )

  post_prob <- if (alternative == "greater") {
    post_prob_greater
  } else {
    1 - .post_prob_bin_from_posts(
      post_t = post_t,
      post_c = post_c,
      margin = -margin,
      rel.tol = rel.tol
    )
  }

  mean_t <- sum(post_t$w * .beta_mean_vec(post_t$a, post_t$b))
  var_t  <- sum(post_t$w * (.beta_var_vec(post_t$a, post_t$b) +
                              .beta_mean_vec(post_t$a, post_t$b)^2)) - mean_t^2

  mean_c <- sum(post_c$w * .beta_mean_vec(post_c$a, post_c$b))
  var_c  <- sum(post_c$w * (.beta_var_vec(post_c$a, post_c$b) +
                              .beta_mean_vec(post_c$a, post_c$b)^2)) - mean_c^2

  post_mean <- mean_t - mean_c
  post_var  <- var_t + var_c

  list(
    post_prob = post_prob,
    post_prob_greater = post_prob_greater,
    post_mean = post_mean,
    post_var = post_var,
    decision = as.numeric(post_prob > cutoff),
    weight = weight_used,
    method = method,
    alternative = alternative,
    margin = margin,
    x.t = x.t,
    x = x
  )
}