R/fit_single_exp_nocovar.R

In BPrinStratTTE: Causal Effects in Principal Strata Defined by Antidrug Antibodies

#' Fit single model to data from a two-arm trial with an exponentially distributed time-to-event endpoint and no predictor of the intercurrent event
#'
#' @param data Data frame of a structure as generated by `sim_dat_one_trial_exp_nocovar()`.
#' @param params List, containing model parameters:
#'   * `tg` Positive integer value, number of intervals to calculate restricted mean survival time using the trapezoidal rule.
#'   * `prior_piT` Numeric vector of length 2, containing parameters (alpha, beta) of the beta prior on pi, indicating the probability of belonging to the stratum of subjects developing the intercurrent event if given treatment.  
#'   * `prior_0N` Numeric vector of length 2, containing parameters (alpha, beta) of the gamma prior on lambda_0N.
#'   * `prior_1N` Numeric vector of length 2, containing parameters (alpha, beta) of the gamma prior on lambda_1N.
#'   * `prior_0T` Numeric vector of length 2, containing parameters (alpha, beta) of the gamma prior on lambda_0T.
#'   * `prior_1T` Numeric vector of length 2, containing parameters (alpha, beta) of the gamma prior on lambda_1T.
#'   * `t_grid` Numeric vector of length `tg`, containing time points defining the time grid (in months) to calculate restricted mean survival time using the trapezoidal rule.
#'   * `chains` Positive integer value, specifying the number of Markov chains.
#'   * `n_iter` Positive integer value, specifying the number of iterations for each chain (including warmup).
#'   * `warmup` Positive integer value, specifying the number of warmup (aka burnin) iterations per chain.
#'   * `cores` Positive integer value, specifying the number of cores to use when executing the chains in parallel.
#'   * `open_progress` Logical value, indicating whether the progress of the chains will be redirected to a file that is automatically opened for inspection.
#'   * `show_messages` Logical value, indicating whether to print the summary of informational messages.
#' @param summarize_fit Logical, if `TRUE` (default), the output is restricted to a summary of results on key parameters over all chains, if `FALSE`, the complete `stanfit` object is returned.
#'
#' @return `tibble()` containing a summary of results on key parameters, or a `stanfit` object, depending on `summarize_fit`.
#' @export
#'
#' @details
#' The data supplied as `params` are used either as priors (`prior_delta`, `prior_0N`, `prior_1N`, `prior_1T`), to inform the model setup (`tg`, `p`, `t_grid`), or as parameters to `rstan::sampling()` which is invoked internally (`chains`, `n_iter`, `warmup`, `cores`, `open_progress`, `show_messages`).
#' 
#' @seealso [fit_single_exp_covar()] and [rstan::sampling()]
#' 
#' @examples
#' d_params_nocovar <- list(
#'   n = 500L,
#'   nt = 250L,
#'   prob_ice = 0.5,
#'   fu_max = 336L,
#'   T0T_rate = 0.2,
#'   T0N_rate = 0.2,
#'   T1T_rate = 0.15,
#'   T1N_rate = 0.1
#' )
#' dat_single_trial <- sim_dat_one_trial_exp_nocovar(
#'   n = d_params_nocovar[["n"]], 
#'   nt = d_params_nocovar[["nt"]],
#'   prob_ice = d_params_nocovar[["prob_ice"]],
#'   fu_max = d_params_nocovar[["fu_max"]],  
#'   T0T_rate = d_params_nocovar[["T0T_rate"]],
#'   T0N_rate = d_params_nocovar[["T0N_rate"]],
#'   T1T_rate = d_params_nocovar[["T1T_rate"]],
#'   T1N_rate = d_params_nocovar[["T1N_rate"]] 
#' )
#' m_params_nocovar <- list(
#'   tg = 48L,
#'   prior_piT = c(0.5, 0.5),
#'   prior_0N = c(1.5, 5),
#'   prior_1N = c(1.5, 5),
#'   prior_0T = c(1.5, 5),
#'   prior_1T = c(1.5, 5),
#'   t_grid =  seq(7, 7 * 48, 7) / 30,
#'   chains = 2L,
#'   n_iter = 3000L,
#'   warmup = 1500L,
#'   cores = 2L,
#'   open_progress = FALSE,
#'   show_messages = TRUE
#' )
#' \donttest{
#' fit_single <- fit_single_exp_nocovar(
#'   data = dat_single_trial,
#'   params = m_params_nocovar,
#'   summarize_fit = TRUE
#' )
#' print(fit_single)
#' }
fit_single_exp_nocovar <- function(data, params, summarize_fit = TRUE) {
  # input data for model
  data_stan <- list(
    n = nrow(data),
    tg = params[["tg"]],
    Z = data$Z,
    S = data$S,
    TIME = data$TIME/30,
    EVENT = data$EVENT,
    prior_piT = params[["prior_piT"]],
    prior_0N = params[["prior_0N"]],
    prior_1N = params[["prior_1N"]],
    prior_0T = params[["prior_0T"]],
    prior_1T = params[["prior_1T"]],
    t_grid = params[["t_grid"]]
  )
  # fit model
  fit_stan <- rstan::sampling(
    object = stanmodels$m_exp_nocovar,
    data   = data_stan,
    iter   = params[["n_iter"]],
    warmup = params[["warmup"]],
    chains = params[["chains"]],
    cores  = params[["cores"]],
    open_progress = params[["open_progress"]],
    show_messages = params[["show_messages"]]
  )
  # for use with .stan files:
  # fit_stan <- rstan::stan(
  #  file   = model,
  #  data   = data_stan,
  #  iter   = params[["n_iter"]],
  #  warmup = params[["warmup"]],
  #  chains = params[["chains"]],
  #  cores  = params[["cores"]]
  #)
  if(isTRUE(summarize_fit)) {
    fit_stan <- fit_stan %>% rstan::summary() %>% magrittr::extract2("summary")
    patterns <- c("S_", "lp", "n_eff")
    fit_stan <- tibble::as_tibble(fit_stan, rownames="var") %>%
      dplyr::filter(!grepl(paste(patterns, collapse="|"), var)) %>%
      dplyr::select(!c("se_mean", "sd", "25%", "75%"))   
    }
  return(fit_stan)
}