R/CumhazEval.R

In BayesSurvival: Bayesian Survival Analysis for Right Censored Data

#' Evaluate whether a true cumulative hazard 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) cumulative hazard function is
#' contained in the credible set generated by the function
#' \link{BayesSurv}.
#'
#' @seealso \link{BayesSurv}, which computes the posterior mean of the
#' cumulative hazard 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 cumulative
#' hazard.
#' @param true.cumhaz The true cumulative hazard function.
#' @param post.mean The posterior mean of the cumulative hazard, given
#' as a function.
#' @param radius The radius of the credible set for the cumulative hazard.
#'
#' @return
#' \item{covered}{Indicator whether the true cumulative hazard 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} )
#' sim.df <- data.frame(id = 1:1000)
#' df <- simsurv(x = sim.df, maxt = 1, hazard = hazard.true)
#'
#' bs <- BayesSurv(df, "eventtime", "status")
#' K <- length(bs$haz.post.mean)
#' cumhaz.pm <- approxfun(c(0, (bs$time.max/K)*(1:K) ), c(0, cumsum(bs$haz.post.mean*bs$time.max/K)))
#' CumhazEval(bs$surv.eval.grid, cumhaz.true, cumhaz.pm, bs$cumhaz.radius)
#'
#'
#' @export

CumhazEval <- function(time.grid, true.cumhaz, post.mean, radius){
  true <- true.cumhaz(time.grid)
  pm <- post.mean(time.grid)
  low <- pmax(0, pm - radius)
  up <- 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)
}