R/SurvEval.R

In BayesSurvival: Bayesian Survival Analysis for Right Censored Data

#' Evaluate whether a true survival function is contained
#' in the credible set.
#'
#' This function is intended to evaluate the Bayesian procedure in a
#' simulation study. To that end, this function can be used to check
#' whether the true (user-defined) survival function is
#' contained in the credible set generated by the function
#' \link{BayesSurv}.
#'
#' @seealso \link{BayesSurv}, which computes the posterior mean of the
#' survival function as well as the radius for its credible set.
#'
#' @references Castillo and Van der Pas (2020). Multiscale Bayesian survival
#'   analysis. <arXiv:2005.02889>.
#'
#' @param time.grid The time grid on which to evaluate the survival
#' function.
#' @param true.surv The true survival function.
#' @param post.mean The posterior mean of the survival function, given
#' as a function.
#' @param radius The radius of the credible set for the survival function
#'
#' @return \item{covered}{Indicator whether the true survival function is
#' completely covered by the credible set on the times contained in
#' \code{time.grid}. 0 = not completely covered, 1 = completely covered.}
#'
#' @examples #Demonstration on a simulated data set
#' library(simsurv)
#' library(ggplot2)
#' hazard.true <- function(t,x, betas, ...){1.2*(5*(t+0.05)^3 - 10*(t+0.05)^2 + 5*(t+0.05) ) + 0.7}
#' cumhaz.true <- Vectorize( function(t){integrate(hazard.true, 0, t)$value} )
#' surv.true <- function(t){exp(-cumhaz.true(t))}
#'
#' sim.df <- data.frame(id = 1:1000)
#' df <- simsurv(x = sim.df, maxt = 1, hazard = hazard.true)
#'
#' bs <- BayesSurv(df, "eventtime", "status")
#' surv.pm <- approxfun(bs$surv.eval.grid, bs$surv.post.mean)
#' SurvEval(bs$surv.eval.grid, surv.true, surv.pm, bs$surv.radius)
#'
#'
#' @export


SurvEval <- function(time.grid, true.surv, post.mean, radius){
  true <- true.surv(time.grid)
  pm <- post.mean(time.grid)
  low <- pmax(0, pm - radius)
  up <- pmin(1, pm + radius)

  lower.check <- as.numeric(low <= true)
  upper.check <- as.numeric(up >= true)

  covered.vec <- as.numeric(lower.check == 1 & upper.check == 1)
  covered <- as.numeric(sum(covered.vec) == length(time.grid))

  return(covered)
}