Implement (cause-specific) competing risks

DESCRIPTION

@@ -17,9 +17,12 @@ License: GPL-3
 Imports: Rcpp (>= 0.11.2), Matrix
 LinkingTo: Rcpp, RcppEigen
 Depends: R (>= 3.1)
-Suggests: 
+Suggests:
     survival,
-    testthat
+    testthat,
+    knitr,
+    rmarkdown
+VignetteBuilder: knitr
 Encoding: UTF-8
 RoxygenNote: 7.3.3
 URL: https://imbs-hl.github.io/ranger/,

R/predict.R

@@ -56,8 +56,9 @@
 ##'   \tabular{ll}{
 ##'       \code{predictions}    \tab Predicted classes/values (only for classification and regression)  \cr
 ##'       \code{unique.death.times} \tab Unique death times (only for survival). \cr
-##'       \code{chf} \tab Estimated cumulative hazard function for each sample (only for survival). \cr
-##'       \code{survival} \tab Estimated survival function for each sample (only for survival). \cr
+##'       \code{chf} \tab Estimated cumulative hazard function for each sample (only for survival). For competing risks, a list of matrices, one per event type. \cr
+##'       \code{survival} \tab Estimated survival function for each sample (only for survival, not for competing risks). \cr
+##'       \code{cif} \tab Estimated cumulative incidence function for each sample (only for competing risks). A list of matrices, one per event type. \cr
 ##'       \code{num.trees}   \tab Number of trees. \cr
 ##'       \code{num.independent.variables} \tab Number of independent variables. \cr
 ##'       \code{treetype}    \tab Type of forest/tree. Classification, regression or survival. \cr
@@ -323,9 +324,35 @@ predict.ranger.forest <- function(object, data, predict.all = FALSE,
       ## Empty
     } else if (forest$treetype == "Survival") {
       result$unique.death.times <- forest$unique.death.times
-      result$chf <- result$predictions
-      result$predictions <- NULL
-      result$survival <- exp(-result$chf)
+      num.event.types <- forest$num.event.types
+      if (is.null(num.event.types) || num.event.types <= 1) {
+        result$chf <- result$predictions
+        result$predictions <- NULL
+        result$survival <- exp(-result$chf)
+      } else {
+        num.timepoints <- length(forest$unique.death.times)
+        if (predict.all) {
+          # predict.all: predictions is 3d array (sample x (event*time) x tree)
+          result$chf <- lapply(seq_len(num.event.types), function(e) {
+            cols <- ((e - 1) * num.timepoints + 1):(e * num.timepoints)
+            result$predictions[, cols, , drop = FALSE]
+          })
+        } else {
+          # Standard prediction: matrix (sample x (event*time))
+          result$chf <- lapply(seq_len(num.event.types), function(e) {
+            cols <- ((e - 1) * num.timepoints + 1):(e * num.timepoints)
+            result$predictions[, cols, drop = FALSE]
+          })
+        }
+        result$cif <- compute_cif(
+          lapply(result$chf, function(m) {
+            if (length(dim(m)) == 3) apply(m, c(1, 2), mean) else m
+          }),
+          num.event.types, num.timepoints
+        )
+        result$predictions <- NULL
+        result$num.event.types <- num.event.types
+      }
     } else if (forest$treetype == "Probability estimation") {
       if (predict.all) {
         ## Set colnames and sort by levels
@@ -458,8 +485,9 @@ predict.ranger.forest <- function(object, data, predict.all = FALSE,
 ##'   \tabular{ll}{
 ##'       \code{predictions}    \tab Predicted classes/values (only for classification and regression)  \cr
 ##'       \code{unique.death.times} \tab Unique death times (only for survival). \cr
-##'       \code{chf} \tab Estimated cumulative hazard function for each sample (only for survival). \cr
-##'       \code{survival} \tab Estimated survival function for each sample (only for survival). \cr
+##'       \code{chf} \tab Estimated cumulative hazard function for each sample (only for survival). For competing risks, a list of matrices, one per event type. \cr
+##'       \code{survival} \tab Estimated survival function for each sample (only for survival, not for competing risks). \cr
+##'       \code{cif} \tab Estimated cumulative incidence function for each sample (only for competing risks). A list of matrices, one per event type. \cr
 ##'       \code{num.trees}   \tab Number of trees. \cr
 ##'       \code{num.independent.variables} \tab Number of independent variables. \cr
 ##'       \code{treetype}    \tab Type of forest/tree. Classification, regression or survival. \cr

R/print.R

@@ -48,6 +48,9 @@ print.ranger <- function(x, ...) {
   cat("Splitrule:                       ", x$splitrule, "\n")
   if (x$treetype == "Survival") {
     cat("Number of unique death times:    ", length(x$unique.death.times), "\n")
+    if (!is.null(x$num.event.types) && x$num.event.types > 1) {
+      cat("Number of event types:           ", x$num.event.types, "\n")
+    }
   }
   if (!is.null(x$splitrule) && x$splitrule == "extratrees" && !is.null(x$num.random.splits)) {
     cat("Number of random splits:         ", x$num.random.splits, "\n")

R/ranger.R

@@ -29,7 +29,7 @@
 ##' Ranger is a fast implementation of random forests (Breiman 2001) or recursive partitioning, particularly suited for high dimensional data.
 ##' Classification, regression, and survival forests are supported.
 ##' Classification and regression forests are implemented as in the original Random Forest (Breiman 2001), survival forests as in Random Survival Forests (Ishwaran et al. 2008).
-##' Includes implementations of extremely randomized trees (Geurts et al. 2006) and quantile regression forests (Meinshausen 2006). 
+##' Includes implementations of extremely randomized trees (Geurts et al. 2006) and quantile regression forests (Meinshausen 2006).
 ##'
 ##' The tree type is determined by the type of the dependent variable.
 ##' For factors classification trees are grown, for numeric values regression trees and for survival objects survival trees.
@@ -38,6 +38,12 @@
 ##' For Survival the log-rank test, a C-index based splitting rule (Schmid et al. 2015) and maximally selected rank statistics (Wright et al. 2016) are available.
 ##' For all tree types, forests of extremely randomized trees (Geurts et al. 2006) can be grown.
 ##'
+##' Competing risks survival forests are supported for status variables with more than one event type (status = 0 for censoring, 1, 2, ..., K for K event types).
+##' The implementation follows the cause-specific hazard approach of Ishwaran et al. (2014), similar to the \code{randomForestSRC} package.
+##' Splitting uses a weighted composite of cause-specific log-rank statistics across all event types.
+##' Predictions include cause-specific cumulative hazard functions (CHF) and cumulative incidence functions (CIF) estimated via the Aalen-Johansen estimator.
+##' See the \code{competing_risks} vignette for details and examples.
+##'
 ##' With the \code{probability} option and factor dependent variable a probability forest is grown.
 ##' Here, the node impurity is used for splitting, as in classification forests.
 ##' Predictions are class probabilities for each sample.
@@ -122,7 +128,7 @@
 ##' @param class.weights Weights for the outcome classes (in order of the factor levels) in the splitting rule (cost sensitive learning). Classification and probability prediction only. For classification the weights are also applied in the majority vote in terminal nodes.
 ##' @param splitrule Splitting rule. For classification and probability estimation "gini", "extratrees" or "hellinger" with default "gini".
 ##'   For regression "variance", "extratrees", "maxstat", "beta" or "poisson" with default "variance".
-##'   For survival "logrank", "extratrees", "C" or "maxstat" with default "logrank". 
+##'   For survival "logrank", "extratrees", "C" or "maxstat" with default "logrank". For competing risks, "C" is not supported.
 ##' @param num.random.splits For "extratrees" splitrule.: Number of random splits to consider for each candidate splitting variable.
 ##' @param alpha For "maxstat" splitrule: Significance threshold to allow splitting.
 ##' @param minprop For "maxstat" splitrule: Lower quantile of covariate distribution to be considered for splitting.
@@ -148,10 +154,10 @@
 ##' @param seed Random seed. Default is \code{NULL}, which generates the seed from \code{R}. Set to \code{0} to ignore the \code{R} seed. 
 ##' @param na.action Handling of missing values. Set to "na.learn" to internally handle missing values (default, see below), to "na.omit" to omit observations with missing values and to "na.fail" to stop if missing values are found.
 ##' @param dependent.variable.name Name of dependent variable, needed if no formula given. For survival forests this is the time variable.
-##' @param status.variable.name Name of status variable, only applicable to survival data and needed if no formula given. Use 1 for event and 0 for censoring.
+##' @param status.variable.name Name of status variable, only applicable to survival data and needed if no formula given. Use 1 for event and 0 for censoring. For competing risks, use 0 for censoring and 1, 2, ..., K for K event types.
 ##' @param classification Set to \code{TRUE} to grow a classification forest. Only needed if the data is a matrix or the response numeric. 
 ##' @param x Predictor data (independent variables), alternative interface to data with formula or dependent.variable.name.
-##' @param y Response vector (dependent variable), alternative interface to data with formula or dependent.variable.name. For survival use a \code{Surv()} object or a matrix with time and status.
+##' @param y Response vector (dependent variable), alternative interface to data with formula or dependent.variable.name. For survival use a \code{Surv()} object or a matrix with time and status. For competing risks, status values should be 0 (censored) or 1, 2, ..., K for K event types.
 ##' @param ... Further arguments passed to or from other methods (currently ignored).
 ##' @return Object of class \code{ranger} with elements
 ##'   \item{\code{forest}}{Saved forest (If write.forest set to TRUE). Note that the variable IDs in the \code{split.varIDs} object do not necessarily represent the column number in R.}
@@ -162,8 +168,10 @@
 ##'   \item{\code{r.squared}}{R squared. Also called explained variance or coefficient of determination (regression only). Computed on out-of-bag data.}
 ##'   \item{\code{confusion.matrix}}{Contingency table for classes and predictions based on out-of-bag samples (classification only).}
 ##'   \item{\code{unique.death.times}}{Unique death times (survival only).}
-##'   \item{\code{chf}}{Estimated cumulative hazard function for each sample (survival only).}
-##'   \item{\code{survival}}{Estimated survival function for each sample (survival only).}
+##'   \item{\code{chf}}{Estimated cumulative hazard function for each sample (survival only). For competing risks, a list of matrices, one per event type.}
+##'   \item{\code{survival}}{Estimated survival function for each sample (survival only). Not returned for competing risks.}
+##'   \item{\code{cif}}{Estimated cumulative incidence function for each sample (competing risks only). A list of matrices, one per event type, computed via the Aalen-Johansen estimator.}
+##'   \item{\code{num.event.types}}{Number of competing event types (survival only, 1 for standard survival).}
 ##'   \item{\code{call}}{Function call.}
 ##'   \item{\code{num.trees}}{Number of trees.}
 ##'   \item{\code{num.independent.variables}}{Number of independent variables.}
@@ -201,6 +209,15 @@
 ##' rg.veteran <- ranger(Surv(time, status) ~ ., data = veteran)
 ##' plot(rg.veteran$unique.death.times, rg.veteran$survival[1,])
 ##'
+##' ## Competing risks survival forest
+##' ## Status: 0 = censored, 1 = event type 1, 2 = event type 2
+##' require(survival)
+##' dat_cr <- data.frame(time = rexp(100), status = sample(0:2, 100, replace = TRUE),
+##'                      x1 = rnorm(100), x2 = rnorm(100))
+##' rg.cr <- ranger(Surv(time, status) ~ ., data = dat_cr)
+##' rg.cr$num.event.types
+##' rg.cr$cif  ## list of CIF matrices, one per event type
+##'
 ##' ## Alternative interfaces (same results)
 ##' ranger(dependent.variable.name = "Species", data = iris)
 ##' ranger(y = iris[, 5], x = iris[, -5])
@@ -223,7 +240,8 @@
 ##'   \item Wright, M. N., Dankowski, T. & Ziegler, A. (2017). Unbiased split variable selection for random survival forests using maximally selected rank statistics. Stat Med 36:1272-1284. \doi{10.1002/sim.7212}.
 ##'   \item Nembrini, S., Koenig, I. R. & Wright, M. N. (2018). The revival of the Gini Importance? Bioinformatics. \doi{10.1093/bioinformatics/bty373}.
 ##'   \item Breiman, L. (2001). Random forests. Mach Learn, 45:5-32. \doi{10.1023/A:1010933404324}. 
-##'   \item Ishwaran, H., Kogalur, U. B., Blackstone, E. H., & Lauer, M. S. (2008). Random survival forests. Ann Appl Stat 2:841-860. \doi{10.1097/JTO.0b013e318233d835}. 
+##'   \item Ishwaran, H., Kogalur, U. B., Blackstone, E. H., & Lauer, M. S. (2008). Random survival forests. Ann Appl Stat 2:841-860. \doi{10.1097/JTO.0b013e318233d835}.
+##'   \item Ishwaran, H., Gerds, T. A., Kogalur, U. B., Moore, R. D., Gange, S. J. & Lau, B. M. (2014). Random survival forests for competing risks. Biostatistics 15:757-773. \doi{10.1093/biostatistics/kxu010}.
 ##'   \item Malley, J. D., Kruppa, J., Dasgupta, A., Malley, K. G., & Ziegler, A. (2012). Probability machines: consistent probability estimation using nonparametric learning machines. Methods Inf Med 51:74-81. \doi{10.3414/ME00-01-0052}.
 ##'   \item Hastie, T., Tibshirani, R., Friedman, J. (2009). The Elements of Statistical Learning. Springer, New York. 2nd edition.
 ##'   \item Geurts, P., Ernst, D., Wehenkel, L. (2006). Extremely randomized trees. Mach Learn 63:3-42. \doi{10.1007/s10994-006-6226-1}.
@@ -297,7 +315,17 @@ ranger <- function(formula = NULL, data = NULL, num.trees = 500, mtry = NULL,
           y <- data[, dependent.variable.name, drop = TRUE]
           x <- data[, !(colnames(data) %in% dependent.variable.name), drop = FALSE]
         } else {
-          y <- survival::Surv(data[, dependent.variable.name], data[, status.variable.name]) 
+          status_vals <- data[, status.variable.name]
+          ## Validate status values before creating Surv object
+          max_status <- max(status_vals)
+          if (max_status >= 1 && !all(status_vals %in% 0:max_status)) {
+            stop("Error: Status values must be 0 (censored) or consecutive integers 1, 2, ..., K for K event types.")
+          }
+          if (max_status > 1) {
+            y <- survival::Surv(data[, dependent.variable.name], status_vals, type = "mstate")
+          } else {
+            y <- survival::Surv(data[, dependent.variable.name], status_vals)
+          }
           x <- data[, !(colnames(data) %in% c(dependent.variable.name, status.variable.name)), drop = FALSE]
         }
       }
@@ -309,10 +337,33 @@ ranger <- function(formula = NULL, data = NULL, num.trees = 500, mtry = NULL,
       if (ncol(data) > 10000) {
         warning("Avoid the formula interface for high-dimensional data. If ranger is slow or you get a 'protection stack overflow' error, consider the x/y or dependent.variable.name interface (see examples).")
       }
+      ## Check if formula involves Surv() with potential competing risks status
+      ## If so, pre-detect and modify the call to use type = "mstate"
+      all_formula_vars <- all.vars(formula)
+      dependent.variable.name <- all_formula_vars[1]
+      lhs <- formula[[2]]
+      is_surv_call <- is.call(lhs) && as.character(lhs[[1]]) == "Surv"
+
+      if (is_surv_call && length(all_formula_vars) >= 2) {
+        status_var <- all_formula_vars[2]
+        if (status_var %in% colnames(data)) {
+          raw_status <- data[, status_var]
+          max_status <- max(raw_status)
+          ## Validate status values before creating Surv object
+          if (max_status >= 1 && !all(raw_status %in% 0:max_status)) {
+            stop("Error: Status values must be 0 (censored) or consecutive integers 1, 2, ..., K for K event types.")
+          }
+          if (max_status > 1) {
+            ## Modify formula LHS to use Surv(..., type = "mstate")
+            lhs[["type"]] <- "mstate"
+            formula[[2]] <- lhs
+          }
+        }
+      }
+
       data.selected <- parse.formula(formula, data, env = parent.frame())
-      dependent.variable.name <- all.vars(formula)[1]
       if (inherits(data.selected[, 1], "Surv")) {
-        status.variable.name <- all.vars(formula)[2]
+        status.variable.name <- all_formula_vars[2]
       }
       y <- data.selected[, 1]
       x <- data.selected[, -1, drop = FALSE]
@@ -1007,10 +1058,19 @@ ranger <- function(formula = NULL, data = NULL, num.trees = 500, mtry = NULL,
     order.snps <- FALSE
   }
 
-  ## No competing risks check
+  ## Competing risks: detect number of event types
+  num.event.types <- 1L
   if (treetype == 5) {
-    if (!all(y.mat[, 2] %in% 0:1)) {
-      stop("Error: Competing risks not supported yet. Use status=1 for events and status=0 for censoring.")
+    status_vals <- y.mat[, 2]
+    num.event.types <- as.integer(max(status_vals))
+    if (num.event.types < 1) {
+      stop("Error: No events in survival data. All status values are 0.")
+    }
+    if (!all(status_vals %in% 0:num.event.types)) {
+      stop("Error: Status values must be 0 (censored) or consecutive integers 1, 2, ..., K for K event types.")
+    }
+    if (num.event.types > 1 && splitrule.num %in% c(2, 3)) {
+      stop("Error: C-index (AUC) splitting not supported for competing risks. Use 'logrank', 'extratrees', or 'maxstat'.")
     }
   }
 
@@ -1063,13 +1123,35 @@ ranger <- function(formula = NULL, data = NULL, num.trees = 500, mtry = NULL,
   } else if (treetype == 5 && oob.error) {
     if (is.list(result$predictions)) {
       result$predictions <- do.call(rbind, result$predictions)
-    } 
+    }
     if (is.vector(result$predictions)) {
       result$predictions <- matrix(result$predictions, nrow = 1)
     }
-    result$chf <- result$predictions
-    result$predictions <- NULL
-    result$survival <- exp(-result$chf)
+
+    num.event.types.result <- result$num.event.types
+    if (is.null(num.event.types.result) || num.event.types.result <= 1) {
+      # Standard single-event survival
+      result$chf <- result$predictions
+      result$predictions <- NULL
+      result$survival <- exp(-result$chf)
+    } else {
+      # Competing risks: predictions matrix has num_event_types * num_timepoints columns
+      num.timepoints <- length(result$unique.death.times)
+      n <- nrow(result$predictions)
+
+      # Split into per-event-type CHF matrices
+      result$chf <- lapply(seq_len(num.event.types.result), function(e) {
+        cols <- ((e - 1) * num.timepoints + 1):(e * num.timepoints)
+        result$predictions[, cols, drop = FALSE]
+      })
+
+      # Compute CIF using Aalen-Johansen estimator
+      result$cif <- compute_cif(result$chf, num.event.types.result, num.timepoints)
+
+      result$predictions <- NULL
+      # Overall event-free survival: product of (1 - sum of cause-specific hazards) over time
+      # Not stored as $survival for competing risks since it's ambiguous
+    }
   } else if (treetype == 9 && oob.error) {
     if (is.list(result$predictions)) {
       result$predictions <- do.call(rbind, result$predictions)
@@ -1120,6 +1202,9 @@ ranger <- function(formula = NULL, data = NULL, num.trees = 500, mtry = NULL,
     }
     result$forest$independent.variable.names <- independent.variable.names
     result$forest$treetype <- result$treetype
+    if (treetype == 5) {
+      result$forest$num.event.types <- num.event.types
+    }
     class(result$forest) <- "ranger.forest"
 
     ## Save covariate levels
@@ -1128,6 +1213,11 @@ ranger <- function(formula = NULL, data = NULL, num.trees = 500, mtry = NULL,
     }
   }
 
+  ## Number of event types (for competing risks)
+  if (treetype == 5) {
+    result$num.event.types <- num.event.types
+  }
+
   ## Dependent (and status) variable name
   ## will be NULL only when x/y interface is used
   result$dependent.variable.name <- dependent.variable.name