R/MLB_SLB_anlaysis.R

In BDIB: Bayesian Dynamic Information Borrowing for Longitudinal Data with Informative Dropout

#' The mixture prior for longitudinal data borrowing (MLB) and self-adapting mixture prior for longitudinal data information borrowing (SLB) Analysis under SPM with MNAR dropout
#'
#' \strong{Arms and notation (H/C/T):}
#' \itemize{
#'   \item \code{H} — Historical control arm. External/previous-study control,
#'         used to build an \emph{informative prior} for the concurrent control.
#'   \item \code{C} — Concurrent control arm. In-trial control; first fitted
#'         under a non-informative prior to assess prior–data conflict and
#'         compute an adaptive weight; then re-fitted under an MLB/SLB prior.
#'   \item \code{T} — Treatment arm. In-trial treatment; typically fitted under
#'         a non-informative prior, then compared against the (re)fitted control.
#' }
#'
#' \strong{Target parameters \eqn{\theta_H,\theta_C,\theta_T}:}
#' \itemize{
#'   \item \code{family = "gaussian"} (continuous outcome):
#'         \eqn{\theta} is the total change at the final visit \eqn{t_K},
#'         i.e., \eqn{\theta=\beta_1 t_K} (thus \eqn{\beta_1=\theta/t_K}).
#'   \item \code{family = "bernoulli"} (binary, probit link):
#'         \eqn{\theta} is the marginal response probability at the final visit,
#'         \eqn{\Pr(Y=1 \text{ at } t_K)\in (0,1)}.
#' }
#' We write \eqn{\theta_H}, \eqn{\theta_C}, \eqn{\theta_T} for the historical,
#' concurrent-control, and treatment arms, respectively.
#'
#' \strong{Workflow (continuous or binary endpoints):}
#' \enumerate{
#'   \item Fit \strong{H} under SPM with a vague prior on \eqn{\theta} to obtain posterior draws.
#'   \item Build an \emph{informative prior} \eqn{g_0} for \eqn{\theta_C} via
#'         \code{RBesT::mixfit} on the \eqn{\theta_H} draws.
#'   \item Build a \emph{non-informative prior} \eqn{g_1}:
#'         Gaussian: unit-information prior from \eqn{g_0} (inflate SD by \eqn{\sqrt{n_\mathrm{hist}}});
#'         Bernoulli: \eqn{\mathrm{Beta}(1,1)}.
#'   \item Fit \strong{C} under \eqn{g_1} to get posterior samples.
#'   \item Compute SLB weight via posterior-probability-ratio KDE (see \code{\link{SLB_weight}}).
#'   \item Construct MLB prior (fixed weight) and SLB prior (data-adaptive weight) via \code{\link{Mixprior}}.
#'   \item Fit \strong{C} under MLB/SLB priors; fit \strong{T} under \eqn{g_1}.
#'   \item Posterior comparisons \eqn{\theta_T} vs \eqn{\theta_C} and decisions per \code{direction}.
#' }
#'
#' @param historical_df Long-format data.frame for H with columns \code{id,time,R,Y}.
#' @param concurrent_df Long-format data.frame for C with columns \code{id,time,R,Y}.
#' @param treatment_df  Long-format data.frame for T with columns \code{id,time,R,Y}.
#' @param family Character, \code{"gaussian"} or \code{"bernoulli"} (probit link), passed to SPM.
#' @param delta Clinically significant difference (CSD) used in SLB weight.
#' @param n_hist Positive integer: historical sample size for unit-information prior (continuous only).
#' @param fixed_weight MLB mixing weight \eqn{\omega} (default 0.5).
#' @param mcmc_settings JAGS settings list with fields \code{n.chains,n.iter,n.burnin,n.thin,parallel}.
#' @param decision_MLB A REQUIRED decision function for MLB prior with threshold p that takes two numeric vectors of
#'   posterior samples \eqn{(\theta_T,\theta_C)} and returns logical/0-1 indicating treatment superiority.
#'
#' \strong{Direction of test (how probabilities are computed/decisions typically defined):}
#' \itemize{
#'   \item Larger-is-better:  \code{function(theta_t, theta_c, p=0.95) mean(theta_t > theta_c) > p}
#'   \item Smaller-is-better: \code{function(theta_t, theta_c, p=0.95) mean(theta_t < theta_c) > p}
#' }
#'
#' @param direction_MLB Character scalar specifying the clinical direction of superiority for MLB prior:
#'   \itemize{
#'     \item \code{"lower"}: smaller is better; report \eqn{\Pr(\theta_T < \theta_C)}.
#'     \item \code{"higher"}: larger is better; report \eqn{\Pr(\theta_T > \theta_C)}.
#'   }
#' @param decision_SLB A REQUIRED decision function for SLB prior with threshold p that takes two numeric vectors of
#'   posterior samples \eqn{(\theta_T,\theta_C)} and returns logical/0-1 indicating treatment superiority.
#'
#' \strong{Direction of test (how probabilities are computed/decisions typically defined):}
#' \itemize{
#'   \item Larger-is-better:  \code{function(theta_t, theta_c, p=0.95) mean(theta_t > theta_c) > p}
#'   \item Smaller-is-better: \code{function(theta_t, theta_c, p=0.95) mean(theta_t < theta_c) > p}
#' }
#'
#' @param direction_SLB Character scalar specifying the clinical direction of superiority for SLB prior:
#'   \itemize{
#'     \item \code{"lower"}: smaller is better; report \eqn{\Pr(\theta_T < \theta_C)}.
#'     \item \code{"higher"}: larger is better; report \eqn{\Pr(\theta_T > \theta_C)}.
#'   }
#' @return An object of class \code{"mlbslb_result"}:
#' \itemize{
#'   \item \code{estimates}: posterior means of \eqn{\theta} for each arm/prior (MLB/SLB).
#'   \item \code{weights}: list with \code{MLB_weight} and \code{SLB_weight}.
#'   \item \code{probabilities}: posterior superiority probabilities (per \code{direction}) under MLB/SLB.
#'   \item \code{decisions}: 0/1 decisions under MLB/SLB according to your \code{decision()}. (1 indicates test superiority is established).
#'   \item \code{posterior_samples}: posterior draws used in the analysis.
#' }
#'
#' @export
#'
#' @examples
#' set.seed(1000)
#' t_vec <- seq(4, 48, 4)
#' # --- Continuous outcome ---
#' data_H <- generate_spm(n=60, t_vec, family="gaussian",
#'                        alpha0=-4, alpha1=0.05, rho=-1, sigma_b=1,
#'                        theta= -3.2, sigma_e=1)$dat
#' data_C <- generate_spm(n=30,  t_vec, family="gaussian",
#'                        alpha0=-4, alpha1=0.05, rho=-1, sigma_b=1,
#'                        theta= -3.2, sigma_e=1)$dat
#' data_T <- generate_spm(n=60,  t_vec, family="gaussian",
#'                        alpha0=-5.3, alpha1=0.05, rho=-1.5, sigma_b=1,
#'                        theta= -3.7, sigma_e=1)$dat
#'
#' \donttest{
#' # Full MCMC fitting takes several minutes.
#' ans_1 <- MLB_SLB_anlaysis(
#'   historical_df = data_H, concurrent_df = data_C, treatment_df = data_T,
#'   family="gaussian", delta=0.5, n_hist=length(unique(data_H$id)), fixed_weight=0.5,
#'   mcmc_settings = list(n.chains=3, n.iter=10000, n.burnin=5000, n.thin=2, parallel=FALSE),
#'   decision_MLB = function(theta_t, theta_c, p = 0.95) mean(theta_t < theta_c) > 0.95,
#'   # Smaller-is-better (treatment < control), p = 0.95
#'   direction_MLB =  "lower",
#'   decision_SLB = function(theta_t, theta_c, p = 0.95) mean(theta_t < theta_c) > 0.95,
#'   # Smaller-is-better (treatment < control), p = 0.95
#'   direction_SLB =  "lower"
#' )
#' ans_1$decisions; ans_1$probabilities
#' }
#'
#' # --- Binary outcome ---
#' set.seed(1000)
#' data_H <- generate_spm(n=200, t_vec, family="bernoulli",
#'                        alpha0=-4, alpha1=0.05, rho=-1, sigma_b=1,
#'                        theta=0.40, beta0=-2)$dat
#' data_C <- generate_spm(n=100, t_vec, family="bernoulli",
#'                        alpha0=-4, alpha1=0.05, rho=-1, sigma_b=1,
#'                        theta=0.40, beta0=-2)$dat
#' data_T <- generate_spm(n=200, t_vec, family="bernoulli",
#'                        alpha0=-5.3, alpha1=0.05, rho=-1.5, sigma_b=1,
#'                        theta=0.50, beta0=-1.5)$dat
#'
#' \donttest{
#' # Full MCMC fitting takes several minutes.
#' ans_2 <- MLB_SLB_anlaysis(historical_df = data_H, concurrent_df = data_C, treatment_df = data_T,
#'    family="bernoulli", delta=0.1, n_hist=length(unique(data_H$id)), fixed_weight=0.5,
#'    mcmc_settings = list(n.chains=3, n.iter=10000, n.burnin=5000, n.thin=2, parallel=FALSE),
#'    decision_MLB = function(theta_t, theta_c, p = 0.95) mean(theta_t > theta_c) > 0.95,
#'    # Larger-is-better (treatment > control), p = 0.95
#'    direction_MLB =  "higher",
#'    decision_SLB = function(theta_t, theta_c, p = 0.95) mean(theta_t > theta_c) > 0.95,
#'    # Larger-is-better (treatment > control), p = 0.95
#'    direction_SLB =  "higher"
#' )
#' ans_2$decisions; ans_2$probabilities
#' }
MLB_SLB_anlaysis <- function(historical_df,
                             concurrent_df,
                             treatment_df,
                             family = c("gaussian","bernoulli"),
                             delta,
                             n_hist = NULL,
                             fixed_weight = 0.5,
                             mcmc_settings = list(n.chains=3, n.iter=10000, n.burnin=5000, n.thin=2, parallel=FALSE),
                             decision_MLB,
                             direction_MLB = c("lower","higher"),
                             decision_SLB,
                             direction_SLB = c("lower","higher")) {

  family    <- match.arg(family)
  direction_MLB <- match.arg(direction_MLB)
  direction_SLB <- match.arg(direction_SLB)

  if (!is.numeric(delta) || delta <= 0) stop("`delta` must be a positive numeric CSD.", call. = FALSE)

  # --- Step 1: Fit historical H under a vague prior (via fit_spm_jags default weak prior) ---
  message("Step 1: Fit historical (SPM, vague prior) ...")
  theta_H_draws <- fit_spm_jags(
    df = historical_df,
    family = family,
    prior_theta = NULL,
    params_to_save = NULL,
    n.chains = mcmc_settings$n.chains,
    n.iter   = mcmc_settings$n.iter,
    n.burnin = mcmc_settings$n.burnin,
    n.thin   = mcmc_settings$n.thin,
    parallel = isTRUE(mcmc_settings$parallel)
  )

  # --- Step 2: Build informative prior g0 from historical posterior via mixfit ---
  message("Step 2: Build informative prior from historical posterior draws ...")
  if (family == "gaussian") {
    inf_prior <- RBesT::mixfit(theta_H_draws, type = "norm", Nc = 3)
  } else { # bernoulli
    inf_prior <- RBesT::mixfit(theta_H_draws, type = "beta", Nc = 3)
  }

  # --- Step 3: Build non-informative prior g1 ---
  message("Step 3: Construct non-informative (unit-information / uniform) prior ...")
  if (family == "gaussian") {
    if (is.null(n_hist) || !is.numeric(n_hist) || n_hist <= 0)
      stop("For continuous endpoints, please provide positive `n_hist` for unit-information prior.")
    nf_prior <- unit_information_prior(informative_prior = inf_prior, n_hist = n_hist)
  } else {
    nf_prior <- RBesT::mixbeta(c(1, 1, 1))  # Uniform(0,1)
  }

  # --- Step 4: Fit concurrent control under NP prior (g1) ---
  message("Step 4: Fit concurrent control under NP prior ...")
  cc_NP_draws <- fit_spm_jags(
    df = concurrent_df,
    family = family,
    prior_theta = nf_prior,
    params_to_save = NULL,
    n.chains = mcmc_settings$n.chains,
    n.iter   = mcmc_settings$n.iter,
    n.burnin = mcmc_settings$n.burnin,
    n.thin   = mcmc_settings$n.thin,
    parallel = isTRUE(mcmc_settings$parallel)
  )

  # Historical mean estimate for SLB
  theta_H_mean <- as.numeric(summary(inf_prior)["mean"])

  # --- Step 5: Compute SLB weight via KDE posterior-probability ratio ---
  message("Step 5: Compute SLB weight ...")
  w_SLB <- SLB_weight(theta_cc_np_samples = cc_NP_draws, theta_h = theta_H_mean, delta = delta)

  # --- Step 6: Construct MLB & SLB priors and fit concurrent control under them ---
  message("Step 6: Fit concurrent control under MLB/SLB priors ...")
  MLB_prior <- Mixprior(informative_prior = inf_prior, noninformative_prior = nf_prior, weight = fixed_weight)
  SLB_prior <- Mixprior(informative_prior = inf_prior, noninformative_prior = nf_prior, weight = w_SLB)

  cc_MLB_draws <- fit_spm_jags(
    df = concurrent_df, family = family, prior_theta = MLB_prior,
    n.chains = mcmc_settings$n.chains, n.iter = mcmc_settings$n.iter,
    n.burnin = mcmc_settings$n.burnin, n.thin = mcmc_settings$n.thin,
    parallel = isTRUE(mcmc_settings$parallel)
  )
  cc_SLB_draws <- fit_spm_jags(
    df = concurrent_df, family = family, prior_theta = SLB_prior,
    n.chains = mcmc_settings$n.chains, n.iter = mcmc_settings$n.iter,
    n.burnin = mcmc_settings$n.burnin, n.thin = mcmc_settings$n.thin,
    parallel = isTRUE(mcmc_settings$parallel)
  )

  # --- Step 7: Fit treatment arm under NP prior (g1) ---
  message("Step 7: Fit treatment arm under NP prior ...")
  trt_draws <- fit_spm_jags(
    df = treatment_df,
    family = family,
    prior_theta = nf_prior,
    n.chains = mcmc_settings$n.chains, n.iter = mcmc_settings$n.iter,
    n.burnin = mcmc_settings$n.burnin, n.thin = mcmc_settings$n.thin,
    parallel = isTRUE(mcmc_settings$parallel)
  )

  # --- Step 8: Posterior comparisons ---
  message("Step 8: Posterior comparisons & decisions ...")
  if (direction_MLB == "lower") {
    pr_MLB <- mean(trt_draws < cc_MLB_draws)
  } else {
    pr_MLB <- mean(trt_draws > cc_MLB_draws)
  }
  if (direction_SLB == "lower") {
    pr_SLB <- mean(trt_draws < cc_SLB_draws)
  } else {
    pr_SLB <- mean(trt_draws > cc_SLB_draws)
  }

  dec_MLB <- as.integer(decision_MLB(trt_draws, cc_MLB_draws))
  dec_SLB <- as.integer(decision_SLB(trt_draws, cc_SLB_draws))

  # Collect results
  out <- list(
    estimates = list(
      concurrent_MLB_theta = mean(cc_MLB_draws),
      concurrent_SLB_theta = mean(cc_SLB_draws),
      treatment_theta      = mean(trt_draws)
    ),
    weights = list(
      MLB_weight = fixed_weight,
      SLB_weight   = w_SLB
    ),
    probabilities = list(
      posterior_superiority_prob_MLB = pr_MLB,
      posterior_superiority_prob_SLB = pr_SLB
    ),
    decisions = list(
      decision_MLB = dec_MLB,
      decision_SLB = dec_SLB
    ),
    posterior_samples = list(
      concurrent_MLB = cc_MLB_draws,
      concurrent_SLB = cc_SLB_draws,
      treatment      = trt_draws
    )
  )
  class(out) <- "mlbslb_result"
  return(out)
}