R/fit_single_exp_covar.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 one predictor of the intercurrent event
#'
#' @param data Data frame of a structure as generated by `sim_dat_one_trial_exp_covar()`.
#' @param params List, containing model parameters:
#'   * `tg` Positive integer value, number of intervals to calculate restricted mean survival time using the trapezoidal rule.
#'   * `p`  Positive integer value, number of predictors of the intercurrent event of interest (i.e. the event that determines. the principal stratum membership).
#'   * `prior_delta`  `p`x2 matrix of positive numerical values, containing normal priors (mean and standard deviation) of the model parameter delta.
#'   * `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 (S4 class), 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_nocovar()] and [rstan::sampling()]
#'
#' @examples
#' d_params_covar <- list(
#'   n = 1000,        
#'   nt = 500,       
#'   prob_X1 = 0.4, 
#'   prob_ice_X1 = 0.5, 
#'   prob_ice_X0 = 0.2,
#'   fu_max = 48*7,       
#'   T0T_rate = 0.2,     
#'   T0N_rate = 0.2,     
#'   T1T_rate = 0.15,     
#'   T1N_rate = 0.1
#'  )
#' dat_single_trial <- sim_dat_one_trial_exp_covar(
#'   n = d_params_covar[["n"]], 
#'   nt = d_params_covar[["nt"]],
#'   prob_X1 = d_params_covar[["prob_X1"]],
#'   prob_ice_X1 = d_params_covar[["prob_ice_X1"]],
#'   prob_ice_X0 = d_params_covar[["prob_ice_X0"]],
#'   fu_max = d_params_covar[["fu_max"]],  
#'   T0T_rate = d_params_covar[["T0T_rate"]],
#'   T0N_rate = d_params_covar[["T0N_rate"]],
#'   T1T_rate = d_params_covar[["T1T_rate"]],
#'   T1N_rate = d_params_covar[["T1N_rate"]] 
#' )
#' m_params_covar <- list(
#'   tg = 48,
#'   p = 2, 
#'   prior_delta = matrix(
#'     c(0, 5, 0, 5),
#'     nrow = 2, byrow = TRUE),
#'   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 = 2,
#'   n_iter = 3000,
#'   warmup = 1500,
#'   cores = 2,
#'   open_progress = FALSE,
#'   show_messages = FALSE   
#' )
#' \donttest{
#' fit_single <- fit_single_exp_covar(
#'   data = dat_single_trial,
#'   params = m_params_covar,
#'   summarize_fit = FALSE
#' )
#' print(fit_single) 
#' }
fit_single_exp_covar <- function(data, params, summarize_fit = TRUE) {
  # input data for model
  data_stan <- list(
    n = nrow(data),
    nt = sum(data$Z==1),
    p = params[["p"]], 
    tg = params[["tg"]],
    Z = data$Z,
    S = data$S,
    S_trt = data$S[data$Z==1],
    TIME = data$TIME/30,
    EVENT = data$EVENT,
    X = matrix(c(rep(1,nrow(data)), data$X), nrow = nrow(data)),
    X_trt = matrix(c(rep(1,sum(data$Z==1)), data$X[data$Z==1]), nrow = sum(data$Z==1)),
    prior_delta = params[["prior_delta"]],
    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_covar,
    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)
}