R/CGaPloth.R

In BGPhazard: Markov Beta and Gamma Processes for Modeling Hazard Rates

#' Plots for the Hazard and Survival Funcion Estimates for the Bayesian
#' non-parametric Markov gamma model with covariates in survival analysis.
#' 
#' Plots the resulting hazard function along with the survival function
#' estimate defined by the Markov gamma process with covariates (Nieto-Barajas,
#' 2003).
#' 
#' This function return plots for the resulting hazard rate as it is computed
#' by \code{\link{CGaMRes} and the quantile of Tao specified by the user aswell
#' as an annotation}. In the same plot the credible intervals for both
#' variables are plotted; The mean of Pi is also annotated. Additionally, it
#' plots the survival function with their corresponding credible intervals.
#' 
#' @param M tibble. Contains the output generated by \code{CuMRres}.
#' @param new_obs tibble. The function calculates the hazard rates and survival
#' function estimates for specific individuals expressed in a tibble, the names of the
#' columns have to be the same as the data input.
#' @param type.h character. "segment"= use segments to plot hazard rates,
#' "line" = link hazard rates by a line
#' @param coxSurv logical. Add estimated Survival function with the Cox-Model
#' @param intervals logical. If TRUE, plots confidence bands for the selected functions including Cox-Model.
#' @param confidence Numeric. Confidence level.
#' @param summary logical. If \code{TRUE}, a summary for hazard and survival
#' functions is returned as a tibble.
#' @return \item{SUM.h}{Numeric tibble. Summary for the mean, median, and a
#' \code{confint / 100} confidence interval for each segment of the hazard
#' function. If \code{summary = TRUE}} \item{SUM.S}{Numeric tibble. Summary for
#' the mean, median, and a \code{confint / 100} confidence interval for each
#' segment of the survival function. If \code{summary = TRUE}}
#' @seealso \link{CGaMRes},
#' @references - Nieto-Barajas, L. E. (2003). Discrete time Markov gamma
#' processes and time dependent covariates in survival analysis. \emph{Bulletin
#' of the International Statistical Institute 54th Session}. Berlin. (CD-ROM).
#' 
#' - Nieto-Barajas, L. E. & Walker, S. G. (2002). Markov beta and gamma
#' processes for modelling hazard rates. \emph{Scandinavian Journal of
#' Statistics} \strong{29}: 413-424.
#' @examples
#' 
#' 
#' 
#' ## Simulations may be time intensive. Be patient.
#' 
#'    # ## Example 1
#'    #  data(leukemiaFZ)
#'    #  leukemia1 <- leukemiaFZ
#'    #  leukemia1$wbc <- log(leukemiaFZ$wbc)
#'    #  CGEX1 <- CGaMRes(data = leukemia1, K = 10, iterations = 100, thinning = 1)
#'    #  CGaPloth(CGEX1)
#' 
#' 
#' 
#' @export CGaPloth
CGaPloth <-
  function(M, new_obs = NULL, type.h= "segment", coxSurv = T, intervals = T,
           confidence = 0.95, summary = FALSE) {
    SUM <- CGaLambdaSumm(M, new = new_obs, confidence)
    h <- extract(SUM, "SUM.h")
    S <- extract(SUM, "SUM.S")
    
    v <- list("tao",
              "K",
              "times",
              "delta",
              "data"
    ) %>% purrr::map(~extract(M,.x)) %>% rlang::set_names(c("tao","K","times","delta","data"))
    tao <- v$tao
    K <- v$K
    delta <- v$delta
    times <- v$times
    data <- v$data
    
    if(type.h == "segment") h.graf <- purrr::imap(h,~{
      if(.y == 1) tit <- "Estimate of Baseline Hazard Rates"
      if(.y == 2) tit <- "Estimate of Hazard Rates for median observation"
      if(!.y %in% c(1,2)) tit <- sprintf("Estimate of Hazard Rates for observation %s",.y-2)
      out <- .x %>% ggplot2::ggplot() + 
        ggplot2::geom_segment(ggplot2::aes(x = tao[-(K+1)], xend = tao[-1], 
                                           y = mean, yend = mean, color = "Hazard Function")) +
        ggplot2::scale_color_manual(values = c("black"), limits = "Hazard Function") +
        ggplot2::guides(color = ggplot2::guide_legend(title = "")) +
        ggplot2::xlab("Time") + ggplot2::ylab("Hazard rate") + ggplot2::scale_alpha_continuous(guide = F) + 
        ggplot2::ggtitle(paste0(tit," with intervals at ",confidence * 100,"% of credibility")) +
        ggthemes::theme_tufte() +
        ggplot2::theme(axis.line = ggplot2::element_line(colour = "black"),
                       legend.position="bottom")
      if(intervals){
        out <- out + ggplot2::geom_errorbar(ggplot2::aes(ymin = lower, ymax = upper, x = (tao[-(K+1)] + tao[-1])/2, width = tao[-1]-tao[-(K+1)]), 
                                            alpha = 0.5, color = "gray50")
      }
      return(out)
    }
    )
    
    if(type.h == "line") h.graf <- purrr::imap(h,~{
      if(.y == 1) tit <- "Estimate of Baseline Hazard Rates"
      if(.y == 2) tit <- "Estimate of Hazard Rates for median observation"
      if(!.y %in% c(1,2)) tit <- sprintf("Estimate of hazard rates for observation %s",.y-2)
      out <- .x %>% ggplot2::ggplot() + 
        ggplot2::geom_line(ggplot2::aes(x = (tao[-(K+1)] + tao[-1])/2, y = mean, color = "Hazard Function")) +
        ggplot2::scale_color_manual(values = c("black"), limits = "Hazard Function") +
        ggplot2::guides(color = ggplot2::guide_legend(title = "")) +
        ggplot2::xlab("Time") + ggplot2::ylab("Hazard rate") + ggplot2::scale_alpha_continuous(guide = F) + 
        ggplot2::ggtitle(paste0(tit," with intervals at ",confidence * 100,"%  of credibility")) +
        ggthemes::theme_tufte() +
        ggplot2::theme(axis.line = ggplot2::element_line(colour = "black"),
                       legend.position="bottom")
      if(intervals){
        out <- out + ggplot2::geom_ribbon(ggplot2::aes(x = (tao[-(K+1)] + tao[-1])/2, ymin = lower, ymax = upper), alpha = .5, fill = "gray70")
      }
      return(out)
    }
    
    )
    
    
    S.graf <- purrr::imap(S,~{
      if(.y == 1) tit <- "Estimate of Baseline Survival Function"
      if(.y == 2) tit <- "Estimate of Survival Function for median observation"
      if(!.y %in% c(1,2)) tit <- sprintf("Estimate of Survival Function for observation %s",.y-2)
      out <- .x %>% ggplot2::ggplot() + ggplot2::geom_line(ggplot2::aes(x = t, y = `S^(t)`,color = "Model estimate")) + 
        ggplot2::scale_color_manual(limits = c("Model estimate"),values = c("black")) +
        ggplot2::guides(color = ggplot2::guide_legend(title = "")) +
        ggplot2::scale_y_continuous(limits = c(0,1)) + 
        ggplot2::ggtitle(paste0(tit," with intervals at ", confidence * 100,"%  of credibility")) +
        ggplot2::labs(x = "t",
                      y = expression(S^{(t)})) +
        ggthemes::theme_tufte() +
        ggplot2::theme(axis.line = ggplot2::element_line(colour = "black"),
                       legend.position = "bottom")
      if(intervals){
        out <- out + ggplot2::geom_ribbon(ggplot2::aes(x = t, ymin = lower, ymax = upper), fill = "gray50", alpha = 0.3)
      }
      return(out)
    }
    )
    
    if(coxSurv){
      t<-survival::Surv(times,delta)
      xfitc<-survival::coxph(as.formula(sprintf("t~%s",paste(names(data),collapse = "+"))),data=data)
      data.b <- rep(0,ncol(data)) %>% matrix(nrow = 1,byrow = T) %>% 
        tibble::as_tibble() %>% rlang::set_names(names(data))
      data.m <- data %>% dplyr::summarise_all(median)
      new <- rbind(data.b,data.m,new_obs)
      
      S.graf<- purrr::map(seq_len(nrow(new)),~{
        fit <- survival::survfit(xfitc,newdata=new[.x,],conf.type="none")
        
        
        km.data <- tibble::tibble(time = fit$time,surv = fit$surv)
        if(km.data$time[1]!= 0){
          km.data <- dplyr::bind_rows(tibble::tibble(time = 0, surv = 1),km.data)
        }
        S.graf[[.x]] + ggplot2::geom_step(data = km.data,na.rm = T, ggplot2::aes(x = time,y = surv), color = "#b22222") + 
          ggplot2::scale_color_manual(limits = c("Model estimate","Cox proportional hazards"),
                                      values = c("black","#b22222")) 
      })
    }
    
    if (summary == TRUE) {
      return(list(h.graf, S.graf, SUM))
    } else{
      return(list(h.graf, S.graf))
    }
  }