R/survTreeLaplaceHazards.R

In discSurv: Discrete Time Survival Analysis

#' Laplace Hazards Of Survival Tree For Competing Risks
#'
#' Predicts the Laplace-smoothed hazards of discrete survival tree based on
#' fitted objects of class "rpart" or "ranger". Can be used for single-risk or 
#' competing risk discrete survival data.
#' 
#' @param treeModel Fitted tree object as generated by function 
#' \code{\link[rpart]{rpart}} (class "rpart").
#' @param newdata Data in long format for which hazards are to be computed. Must 
#' contain the same columns that were used for tree fitting (class "data.frame").
#' @param alpha Smoothing parameter for laplace-smoothing. Must be a non-negative 
#' number. A value of 0 corresponds to no smoothing (class "numeric").
#' @param rangerData Original training data (class "data.frame") for fitting 
#' \emph{treeModel} of class "ranger". Must be in long format. 
#' @return A m by k matrix with m being the length of newdata and k being the 
#' number of classes in treeModel. Each row corresponds to the smoothed hazard 
#' of the respective observation.
#' @seealso \code{\link[ranger]{ranger}}
#' @keywords survival
#' @examples
#' ############################################
#' # Example with rpart discrete survival tree
#' library(pec)
#' library(caret)
#'
#' # Example data
#' data(cost)
#'
#' # Convert time to years and select training and testing subsample
#' cost$time <- ceiling(cost$time/365)
#' costTrain <- cost[1:100, ]
#' costTest  <- cost[101:120, ]
#'
#' # Convert to long format
#' timeColumn <- "time"
#' eventColumn <- "status"
#' costTrainLong <- dataLong(dataShort=costTrain, timeColumn = "time", 
#'                           eventColumn = "status")
#' costTestLong  <- dataLong(dataShort=costTest, timeColumn = "time", 
#'                           eventColumn = "status")
#' head(costTrainLong)
#'
#' # Fit a survival tree
#' costTree <- rpart(formula = y ~ timeInt + prevStroke + age + sex, data = costTrainLong, 
#'                   method = "class")
#'
#' # Compute smoothed hazards for test data
#' predictedhazards <- survTreeLaplaceHazard(costTree, costTestLong, 1)
#' predictedhazards
#'
#' ############################################
#' # Example with ranger discrete survival tree
#' library(pec)
#' library(caret)
#' library(ranger)
#' data(cost)
#'
#' # Take subsample and convert time to years
#' cost$time <- ceiling(cost$time/365)
#' costSubTrain <- cost[1:50,]
#' costSubTest <- cost[51:70,]
#'
#' # Specify column names for data augmentation
#' timeColumn<-"time"
#' eventColumn<-"status"
#' costSubTrainLong <- dataLong(costSubTrain, timeColumn, eventColumn)
#' costSubTestLong <- dataLong(costSubTest, timeColumn, eventColumn)
#'
#' # Estimate discrete survival tree
#' formula <- y ~ timeInt + diabetes + prevStroke + age + sex
#' rangerTree <- ranger(formula, costSubTrainLong, num.trees = 1, mtry = 5, 
#' classification = TRUE, splitrule = "hellinger", replace = FALSE, 
#' sample.fraction = 1, max.depth = 5)
#'
#' # Compute laplace-smoothed hazards
#' laplHaz <- survTreeLaplaceHazard(rangerTree, 
#' costSubTestLong, alpha = 1, costSubTrainLong)
#' laplHaz
#' @export survTreeLaplaceHazard
survTreeLaplaceHazard <- function(treeModel, newdata, alpha, rangerData=NULL){
  
  # Input Checks
  if(alpha < 0 | !is.numeric(alpha) | length(alpha) != 1)
  {
    stop("Alpha must be a non-negative number.")
  }

  if( !is.data.frame(newdata) ) {stop("Argument *newdata* is not in the correct format! Please specify as class data.frame.")}
  if( !is.data.frame(rangerData) & !is.null(rangerData) ) {stop("Argument *rangerData* is not in the correct format! Please specify as class data.frame.")}
  
  if( inherits(treeModel, "rpart") ){
    
    if(is.null(treeModel$frame)){
      stop("Incorrect model. Please provide an object of class rpart or ranger.")
    }
    
    if(!all(unique(treeModel$frame$var)[-which(unique(treeModel$frame$var) == "<leaf>")]
            %in% colnames(newdata))){
      stop("Newdata does not contain the same covariates as the tree model.")
    }
    
    #derive number of risks
    n_events <- length(unique(treeModel$y))
    #derive index of terminal nodes
    leaf_index <- factor(rownames(treeModel$frame[which(treeModel$frame$var == "<leaf>"), ]))
    #predict node for new data
    predicted_values <- factor(predict_leaves(treeModel, newdata), levels = leaf_index)
    #compute laplace-smoothed hazards for new data
    y_table <- treeModel$frame$yval2
    y_table <- y_table[leaf_index,2:(1 + n_events)]
    if(is.null(nrow(y_table))) y_table = matrix(y_table, nrow = 1)
    hazards <- t(apply(y_table, 1, function(x) (x + alpha)/(sum(x) + alpha * n_events)))
    hazards_fitted <- hazards[predicted_values, ]
    return(hazards_fitted)
  } else{
    if( inherits(treeModel, "ranger") ){
      
      if(alpha < 0|!is.numeric(alpha)|length(alpha)!=1)
      {
        stop("Alpha must be a non-negative number.")
      }
      if(!all(unique(treeModel$frame$var)[-which(unique(treeModel$frame$var) == "<leaf>")]
              %in% colnames(newdata)))
      {
        stop("Newdata does not contain the same covariates as the tree model.")
      }
      if(is.null(rangerData$y))
      {
        rangerData$y <- as.numeric(factor(rangerData$responses)) - 1
      }
      if( is.factor(rangerData$y) ){
        rangerData$y <- as.numeric(as.character(rangerData$y))
      }
      #get predicted nodes of original data frame
      pred_nodes_rangerData <- cbind(predict(treeModel, rangerData, type = "terminalNodes")$predictions, rangerData$y)
      #get event and node count
      n_ev <- sort(unique(rangerData$y)) - 1
      uq_row <- sort(unique(pred_nodes_rangerData[, 1]))
      #get absolute and laplace-smoothed relative frequency of observations per event per node
      ev_node <- expand.grid(uq_row, n_ev)
      freq_vals <- apply(ev_node, 1, function(x) length(which(rowSums(sweep(pred_nodes_rangerData, 2, x, "!="))==0)))
      freq_table <- matrix(freq_vals + alpha, nrow = length(uq_row), byrow = TRUE)
      freq_table <- t(apply(freq_table, 1, function(x) x / sum(x)))
      #get predicted nodes for new data
      pred_nodes <- factor(predict(treeModel,data = newdata, type = "terminalNodes")$predictions,levels = uq_row)
      hazards_fitted <- freq_table[pred_nodes, ]
      return(hazards_fitted)
    } else{
      stop("Incorrect model. Please provide an object of class 'rpart' or 'ranger'.")
    }
  }
}

#predict nodes from ranger tree object
predict_leaves <-
  function (object, newdata, na.action = na.pass) {
    where <-
      if (missing(newdata)) 
        object$where
    else {
      if (is.null(attr(newdata, "terms"))) {
        Terms <- delete.response(object$terms)
        newdata <- model.frame(Terms, newdata, na.action = na.action, 
                               xlev = attr(object, "xlevels"))
        if (!is.null(cl <- attr(Terms, "dataClasses"))) 
          .checkMFClasses(cl, newdata, TRUE)
      }
      pred.rpart(object, rpart.matrix(newdata))
    }
    as.integer(row.names(object$frame))[where]
  }

pred.rpart <- getFromNamespace("pred.rpart", "rpart")
rpart.matrix <- getFromNamespace("rpart.matrix", "rpart")