R/post_surv_bgam.R

In csetraynor/rstanhaz: Bayesian Survival Models (rstanhaz)

#' Posterior Hazard
#'
#'
#' Bayesian semiparametric survival model with GAMM
#'
#' @param x data
#' @param surv_mod vector with covariates
#' @return mod
#' @seealso [coxph]
#' @keywords coxph
#'
#' @author Carlos S Traynor
#' @references
#'
#'  Terry M. Therneau and Patricia M. Grambsch (2000).
#'   Modeling Survival Data: Extending the Cox Model.
#'   Springer, New York. ISBN 0-387-98784-3.
#'@export post_surv

post_surv <- function(x, status = "status", time = "time", surv_form = NULL, prior = NULL, form = NULL,  ...) {
  rstan::rstan_options(auto_write = TRUE)
  
  if(!is.null(surv_form)) {
    input.formula <- as.formula( paste0(c( "status~1+offset(log_dtime)+s(tend)+", paste(surv_form, collapse = "+")), collapse = ""))
  } else{
    if(!is.null(form) ) {
      input.formula <- form
    } else {
      input.formula <- as.formula( "status~1+offset(log_dtime)+s(tend)")
    }
  }
  
  #prepare long format dataset
  x$status <- unlist(x[status])
  x$time <- unlist(x[time])
  long_x <- gen_long_dat(x, status = status, time = time)
  
  if(!is.null(prior)){
    m1.stan_gam <- rstanarm::stan_gamm4(input.formula, data = long_x ,
                                        #set likelihood (poisson)
                                        family= poisson(),
                                        #set priors (default)
                                         prior = prior,
                                        # prior_intercept = normal(),
                                        # prior_smooth = exponential(autoscale = FALSE),
                                        # prior_aux = exponential(),
                                        # prior_covariance = decov(),
                                        # prior_PD = FALSE,
                                        # #arguments passed to sampling algorithm
                                        # seed = 7,
                                        # # iter = 10000,
                                        # # thin = 5,
                                        refresh = 0,
                                        # # control = list(adapt_delta = 0.9),
                                        chains = 4 )
  } else {
    m1.stan_gam <- rstanarm::stan_gamm4(input.formula, data = long_x ,
                                        #set likelihood (poisson)
                                        family= poisson(),
                                        #set priors (default)
                                        # prior = normal(),
                                        # prior_intercept = normal(),
                                        # prior_smooth = exponential(autoscale = FALSE),
                                        # prior_aux = exponential(),
                                        # prior_covariance = decov(),
                                        # prior_PD = FALSE,
                                        # #arguments passed to sampling algorithm
                                        # seed = 7,
                                        # # iter = 10000,
                                        # # thin = 5,
                                         refresh = 0,
                                        # # control = list(adapt_delta = 0.9),
                                         chains = 4 )
  }
  return(m1.stan_gam)
}

#' Prediction of Survival probability from a Bayesian model
#'
#'
#' Predicts the survival probability for each drawn of the posterior distribution.
#'
#' @param x data
#' @param bgam Bayesian Generalized Additive Models
#' @param preProcess Pre-processing information from training set to Test set.
#' @param ... Ignored
#' @return mod
#' @seealso [gamair]
#' @keywords gamair
#'
#' @author Carlos S Traynor
#' @references
#'
#'  Wood, S.N. (2006) Generalized Additive Models: An Introduction with R 
#'  Chapman & Hall/CRC, Boca Raton, Florida. ISBN 1-58488-474-6
#'@export surv_pred_bgam

surv_pred_bgam <- function(x, bgam, preProcValues = NULL, ...) {
  train_data <- rsample::analysis(x)
  
  #get timepoints
  timepoints <-  seq(0, max(train_data$time),
                     length.out = 100L)
  
  test_data <- rsample::assessment(x)
  
  #Standardise test data 
  if(is.null(preProcValues)){
    preProcValues <- caret::preProcess(train_data[c("age_at_diagnosis", "npi")],
                                       method = c("center", "scale") )
  }

  testTransformed <- predict(preProcValues, test_data[c("age_at_diagnosis", "npi")])
  colnames(testTransformed) = c("age_std", "npi_std")
  
  test_data <- cbind(test_data, testTransformed)
  
  #get new dataset
  newdat <-  gen_new.frame(dat = test_data, timepoints = timepoints)
  
  # get surv formula
  form <-  names(bgam$coefficients)
  toMatch <- c(":", "Intercept", "s\\(time")
  toMatch <- form[ !grepl(paste(toMatch,collapse="|"), 
                               form) ]
  surv_form <- gsub("TRUE", "", toMatch)
  
  newdat <- newdat[ ,match(c(surv_form, "log_dtime", "time", "sample_id"), colnames(newdat))]
  #predict linear predictor
  post <- suppressWarnings(posterior_linpred(bgam, newdata = newdat) )
  pred_frame <- pred_surv(post = post, longdata = newdat)
  
  pred_frame
}
#' Prediction of survival from model
#'
#'
#' Fits bayesian semi-parametric survival models
#'
#' @param x data
#' @param mod Coxph model object fitted with coxph (survival).
#' @return mod
#' @seealso [coxph]
#' @keywords coxph
#'
#' @author Carlos S Traynor
#' @references
#'
#'  Terry M. Therneau and Patricia M. Grambsch (2000).
#'   Modeling Survival Data: Extending the Cox Model.
#'   Springer, New York. ISBN 0-387-98784-3.
#'@export pred_sbm

pred_sbm <- function(x, surv_form, prior = NULL, ...) {
  train_data <- rsample::analysis(x)
  
  model.map2stan <- post_surv(x = train_data , surv_form = surv_form, prior = prior, iter = 12000, thin = 10, adapt_delta = 0.99)

  timepoints <-  seq(0, max(train_data$time),
                     length.out = 100L)
  
  test_data <- rsample::assessment(x)
  
  newdat <-  gen_new.frame(dat = test_data, timepoints = timepoints)
  newdat <- newdat[ ,match(c(surv_form[!grepl(":", surv_form)], "log_dtime", "time", "sample_id"), colnames(newdat))]
  post <- suppressWarnings(posterior_linpred(model.map2stan, newdata = newdat) )
  pred_frame <- pred_surv(post = post, longdata = newdat)

  pred_frame
}


 

#' Get concordance index
#'
#'
#' Get concordance index from a survbayes model
#'
#' @param x data
#' @param mod Coxph model object fitted with coxph (survival).
#' @return mod
#' @seealso [coxph]
#' @keywords coxph
#'
#' @author Carlos S Traynor
#' @references
#'
#'  Terry M. Therneau and Patricia M. Grambsch (2000).
#'   Modeling Survival Data: Extending the Cox Model.
#'   Springer, New York. ISBN 0-387-98784-3.
#'@export  get_ci2

get_ci2 <- function(pred_frame, x ) {
  
  train_data <- rsample::analysis(x)
  
  timepoints <-  seq(0, max(train_data$time),
                     length.out = 100L)
  
  test_data <- rsample::assessment(x)
  newdat <-  gen_new.frame(dat = test_data, timepoints = timepoints)
  newdat <- newdat[ ,match(c(surv_form[!grepl(":", surv_form)], "log_dtime", "time", "sample_id"), colnames(newdat))]
  
  concordance <- sapply(seq_along(colnames(pred_frame)), function(i){
    newdat$surv <- as.vector( pred_frame[ ,i] ) 
    
    probs <- get_survProb(newdat = newdat)
    
    
    tail(suppressWarnings(suppressMessages( (pec::cindex(probs, 
                                           Surv(time, status) ~ 1,
                                        data = test_data,
                                  pred.times = timepoints,
                                  eval.times = timepoints))))$AppCindex$matrix , n=1)
    
  })
}