R/apply_mcmc.R

In psborrow: Bayesian Dynamic Borrowing with Propensity Score

#' Fit Dynamic Borrowing MCMC Model
#'
#' Fit a dynamic borrowing Weibull survival model to the given dataset and extract the posterior
#' samples using MCMC.
#' See the user guide for more information on the model formulation.
#' See [run_mcmc()] for more information on the available parameters for tuning the MCMC sampling
#' process
#'
#' @param dt A data.frame containing data required for modelling. See details
#'
#' @param formula_cov A one sided formula specifying which non-treatment covariates should
#' be included into the model. See details
#'
#' @param object A `apply_mcmc` object created by [apply_mcmc()]
#'
#' @param ... Additional arguments passed onto [run_mcmc()]. Only exception being the
#' `path` argument which is not supported by this function
#'
#' @details
#'
#' ## `apply_mcmc()`
#'
#' The `dt` data.frame must contain 1 row per subject with the following variables:
#'   - **time** - A continuous non-zero number specifying the time that the subject had an event at
#'   - **cnsr** - A column of 0/1's where 1 indicates that the event was censored/right truncated
#'   - **ext** - A column of 0/1's where 1 indicates that the subject was part of the external
#' control
#'   - **trt** - A column of 0/1's where 1 indicates that the subject was receiving the experimental
#' treatment
#'
#' The `dt` data.frame may also contain any additional covariates to be used in the Weibull model
#' as specified by `formula_cov`. In order to fit a valid model `formula_cov` must contain
#' the intercept term. The formula will be automatically adjusted to include the treatment term
#' and as such should not be included here, if you want to include a treatment interaction term
#' this should be done by using `~ trt:covariate` and NOT via `~ trt*covariate`.
#'
#' ## `extract_samples()`
#'
#' This function can be used to extract the samples generated by `apply_mcmc()`
#'
#' ## `summary()`
#'
#' This function provides summary statistics about the samples generated by `apply_mcmc()`
#'
#' ## Extracted Samples
#'
#' The extracted samples can be roughly defined as follows (see the user guide for full details):
#'
#' - **`HR_cc_hc`** - The hazard ratio between the concurrent control arm and the historical
#' control arm. This can be
#' be thought of as the ratio of the scale parameter between the baseline trial distribution and
#' the baseline
#' external control distribution. This is equivalent to `exp(alpha[2] - alpha[1])`
#'
#' - **`HR_trt_cc`** - The hazard ratio between the treatment arm and the concurrent control arm.
#' This is equivalent
#' to `exp(beta_trt)`
#'
#' - **`alpha[1]`** - The shape parameter for the trial's baseline distribution
#'
#' - **`alpha[2]`** - The shape parameter for the historical control's baseline distribution
#'
#' - **`beta_trt`** - The log-hazard ratio for the treatment effect. This is equivalent to
#' `log(HR_trt_cc)`
#'
#' - **`beta_<var>`** - The log-hazard ratio for any other covariate provided to the model via
#' `formula_cov`
#'
#' - **`r0`** - The scale parameter for the baseline distribution of both the trial and the
#' historical control
#'
#' - **`tau/sigma`** - The precision/variance for `alpha[1]` i.e. controls how much
#' information is borrowed from the
#' historical control arm
#'
#' @name apply_mcmc
#' @importFrom stats model.matrix sd quantile
#' @export
apply_mcmc <- function(dt, formula_cov, ...) {
    design_mat <- model.matrix(formula_cov, dt)

    if (!"(Intercept)" %in% colnames(design_mat)) {
        stop("Covariate formula must contain the intercept term")
    }

    if (ncol(design_mat) == 1) {
      stop(paste(
          "The function cannot be run without additional covariates",
          "i.e. `formula_cov = ~ 1 ` is not supported"
      ))
    }

    keep_columns <- colnames(design_mat)[!grepl("\\(Intercept\\)", colnames(design_mat))]
    design_mat <- design_mat[, keep_columns]

    colnames(design_mat) <- paste0("cov", seq_len(ncol(design_mat)))

    dt2 <- dt[, c("time", "ext", "trt", "cnsr")]

    design_df <- as.data.frame(design_mat) %>%
        as_tibble() %>%
        bind_cols(dt2)

    mcmc_results <- add_mcmc(
        dt = as.data.frame(design_df),
        ...
    )[[1]]

    mcmc_results$mcmc.list <- lapply(
        mcmc_results$mcmc.list,
        fix_col_names,
        column_names = keep_columns
    )

    class(mcmc_results) <- "apply_mcmc"

    return(mcmc_results)
}


#' Fix Column Names
#'
#' Utility function to make the mcmc column names more human friendly
#'
#' @param x a mcmc results object created by [add_mcmc()]
#' @param column_names The names to change the beta columns to
#'
fix_col_names <- function(x, column_names) {
    index_trt <- which(colnames(x) == "beta[1]")
    colnames(x)[index_trt] <- "beta_trt"

    index_beta <- grep("beta\\[[0-9]+\\]", colnames(x))
    colnames(x)[index_beta] <- paste0("beta_", column_names)
    return(x)
}



#' @rdname apply_mcmc
#' @export
extract_samples <- function(object) {
    df <- do.call(rbind, object$mcmc.list)
    df
}



#' @rdname apply_mcmc
#' @export
summary.apply_mcmc <- function(object, ...) {
    samples <- extract_samples(object)
    data.frame(
        parameter = colnames(samples),
        mean = apply(samples, 2, mean),
        sd = apply(samples, 2, sd),
        lci = apply(samples, 2, quantile, 0.025),
        uci = apply(samples, 2, quantile, 0.975),
        row.names = NULL
    )
}