R/scr_gamma_frailty.R

In lcomm/rsemicompstan: Semicompeting Risks in Stan

#' Fit semicompeting risks with gamma frailty
#'
#' @param x N x P design matrix, no intercept
#' @param z Length-N vector of binary treatment indicators
#' @param yr Length-N vector of non-terminal event times
#' @param yt Length-N vector of terminal event times
#' @param dyr Length-N vector binary indicators for having observed the 
#' non-terminal event
#' @param dyt Length-N vector binary indicators for having observed the 
#' terminal event
#' @param shared_beta Whether to use transition-specific hazards which force
#' adjustment covariate coefficients to be the same across treatment arms 
#' (\code{shared_beta = TRUE}) or allow all models to have different coefficients
#' (\code{shared_beta = FALSE})
#' @param use_priors Whether to use weakly informative/data-driven priors
#' @param sigma_pa Hyperparameter alpha for inverse gamma prior on sigma
#' @param sigma_pb Hyperparameter beta for inverse gamma prior on sigma. Prior 
#' mean for sigma is beta/(alpha - 1) for alpha > 1 and prior mode is 
#' beta/(alpha + 1).
#' @param ... Additional parameters to pass to `rstan::sampling`
#' @return an object of class `stanfit` returned by `rstan::sampling`
#' @examples 
#' \dontrun{
#' rstan_options(auto_write = TRUE)
#' options(mc.cores = 4)
#' library("rsemicompstan")
#' set.seed(123)
#' N <- 5000
#' x1 <-matrix(rnorm(N), ncol = 1)
#' dat <- SemiCompRisks::simID(x1 = x1, x2 = x1, x3 = x1, 
#'                             beta1.true = 0.1, beta2.true = 0.2, beta3.true = 0.3, 
#'                             alpha1.true = 1, alpha2.true = 0.95, alpha3.true = 1,
#'                             kappa1.true = 0.2, kappa2.true = 0.3, kappa3.true = 0.4,
#'                             theta.true = 0.5, SigmaV.true = NULL,
#'                             cens = c(0.5, 10))
#' z <- rbinom(N, size = 1, prob = 0.5)
#' resg <- scr_gamma_frailty_stan(x = x1, z = z, yr = dat$y1, yt = dat$y2,
#'                                dyr = dat$delta1, dyt = dat$delta2,
#'                                use_priors = TRUE,
#'                                shared_beta = TRUE,
#'                                sigma_pa = 0.6, sigma_pb = 0.6,
#'                                iter = 2000, chains = 4)
#' }
#' @export
scr_gamma_frailty_stan <- function(x, z, yr, yt, dyr, dyt, 
                                   shared_beta = FALSE,
                                   use_priors = TRUE, 
                                   sigma_pa = 0.7, sigma_pb = 0.7, 
                                   mc.cores = 1, ...) {
  mc.cores <- min(mc.cores, parallel::detectCores())
  options(mc.cores = mc.cores)
  if (use_priors) {
    pm <- make_prior_means(yr = yr, yt = yt, dyr = dyr, dyt = dyt)  
  } else {
    pm <- list(log_alpha_pmean = rep(0, 6), log_kappa_pmean = rep(0, 6))
  }
  out <- rstan::sampling(stanmodels$scr_gamma_frailty,
                         data = list(N = NROW(x), 
                                     z = z,
                                     P = NCOL(x),
                                     yr = yr,
                                     yt = yt,
                                     dyr = dyr,
                                     dyt = dyt,
                                     shared_beta = shared_beta * 1,
                                     use_priors = use_priors * 1,
                                     log_alpha_pmean = pm$log_alpha_pmean,
                                     log_kappa_pmean = pm$log_kappa_pmean,
                                     sigma_pa = sigma_pa,
                                     sigma_pb = sigma_pb),
                         ...)
  return(out)
  
}