From d0ad848f2725dcb6ee8b0d174116005cf62d2ec7 Mon Sep 17 00:00:00 2001
From: Lukas Burk <burk@leibniz-bips.de>
Date: Sat, 28 Feb 2026 02:10:39 +0100
Subject: [PATCH 1/3] Add competing risks support via cause-specific hazards

Extend survival forests to handle competing risks data where the status
variable encodes censoring (0) or one of K event types (1, 2, ..., K).
Uses composite sum of squared log-rank statistics across event types for
splitting, following Ishwaran et al. (2014). Returns cause-specific CHFs
as a list of K matrices and cumulative incidence functions (CIF) via the
Aalen-Johansen estimator. Fully backward compatible when K=1.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
---
 R/predict.R                           |  32 ++-
 R/print.R                             |   3 +
 R/ranger.R                            |  92 ++++++-
 R/utility.R                           |  40 ++++
 src/ForestSurvival.cpp                |  67 ++++--
 src/ForestSurvival.h                  |  15 +-
 src/TreeSurvival.cpp                  | 332 ++++++++++++--------------
 src/TreeSurvival.h                    |  25 +-
 src/rangerCpp.cpp                     |   6 +
 tests/testthat/test_competing_risks.R | 166 +++++++++++++
 tests/testthat/test_survival.R        |  13 +-
 11 files changed, 579 insertions(+), 212 deletions(-)
 create mode 100644 tests/testthat/test_competing_risks.R

diff --git a/R/predict.R b/R/predict.R
index e5230db02..9f1aad22d 100644
--- a/R/predict.R
+++ b/R/predict.R
@@ -323,9 +323,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
diff --git a/R/print.R b/R/print.R
index 8f02bf096..14821f00c 100644
--- a/R/print.R
+++ b/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")
diff --git a/R/ranger.R b/R/ranger.R
index 3acf9167d..dfebc6cc4 100644
--- a/R/ranger.R
+++ b/R/ranger.R
@@ -297,7 +297,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 +319,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 +1040,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 +1105,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 +1184,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 +1195,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
diff --git a/R/utility.R b/R/utility.R
index 3d83208ee..d4f3e9d91 100644
--- a/R/utility.R
+++ b/R/utility.R
@@ -31,6 +31,46 @@ integer.to.factor <- function(x, labels) {
   factor(x, levels = seq_along(labels), labels = labels)
 }
 
+# Compute cumulative incidence functions from cause-specific CHFs
+# using the Aalen-Johansen estimator
+# chf_list: list of K matrices (n x num_timepoints), each containing cause-specific CHF
+# Returns: list of K matrices (n x num_timepoints), each containing CIF
+compute_cif <- function(chf_list, num_event_types, num_timepoints) {
+  n <- nrow(chf_list[[1]])
+
+  # Extract cause-specific hazard increments from cumulative hazard
+  # h_e(t) = CHF_e(t) - CHF_e(t-1)
+  hazard_list <- lapply(chf_list, function(chf_mat) {
+    h <- chf_mat
+    if (num_timepoints > 1) {
+      h[, 2:num_timepoints] <- chf_mat[, 2:num_timepoints] - chf_mat[, 1:(num_timepoints - 1)]
+    }
+    h
+  })
+
+  # Overall hazard at each time: sum across event types
+  overall_hazard <- Reduce("+", hazard_list)
+
+  # Compute event-free survival: S(t) = prod(1 - overall_hazard(s), s <= t)
+  # S(t-) = S(t-1) for t > 1, S(0-) = 1
+  one_minus_h <- 1 - overall_hazard
+  surv <- t(apply(one_minus_h, 1, cumprod))
+  if (!is.matrix(surv)) surv <- matrix(surv, nrow = 1)
+
+  # S(t-1): shift right, with S(0-) = 1
+  surv_prev <- cbind(1, surv[, -num_timepoints, drop = FALSE])
+
+  # CIF_e(t) = sum_{s<=t} S(s-) * h_e(s)
+  cif_list <- lapply(hazard_list, function(h) {
+    cif_increments <- surv_prev * h
+    cif <- t(apply(cif_increments, 1, cumsum))
+    if (!is.matrix(cif)) cif <- matrix(cif, nrow = 1)
+    cif
+  })
+
+  cif_list
+}
+
 # Save version of sample() for length(x) == 1
 # See help(sample)
 save.sample <- function(x, ...) {
diff --git a/src/ForestSurvival.cpp b/src/ForestSurvival.cpp
index 9a31f3da3..9f5b64afa 100644
--- a/src/ForestSurvival.cpp
+++ b/src/ForestSurvival.cpp
@@ -30,12 +30,28 @@ void ForestSurvival::loadForest(size_t num_trees, std::vector<std::vector<std::v
   this->unique_timepoints = unique_timepoints;
   data->setIsOrderedVariable(is_ordered_variable);
 
+  // Detect num_event_types from CHF size vs timepoints
+  if (!forest_chf.empty()) {
+    for (const auto& tree_chf : forest_chf) {
+      for (const auto& node_chf : tree_chf) {
+        if (!node_chf.empty()) {
+          num_event_types = node_chf.size() / unique_timepoints.size();
+          break;
+        }
+      }
+      if (num_event_types > 1) break;
+    }
+  }
+  if (num_event_types < 1) {
+    num_event_types = 1;
+  }
+
   // Create trees
   trees.reserve(num_trees);
   for (size_t i = 0; i < num_trees; ++i) {
     trees.push_back(
         std::make_unique<TreeSurvival>(forest_child_nodeIDs[i], forest_split_varIDs[i], forest_split_values[i],
-            forest_chf[i], &this->unique_timepoints, &response_timepointIDs));
+            forest_chf[i], &this->unique_timepoints, &response_timepointIDs, num_event_types));
   }
 
   // Create thread ranges
@@ -107,6 +123,23 @@ void ForestSurvival::initInternal() {
     min_bucket[0] = DEFAULT_MIN_BUCKET_SURVIVAL;
   }
 
+  // Detect number of event types from data if not already set (only in grow mode)
+  if (!prediction_mode && num_event_types <= 1) {
+    double max_status = 0;
+    for (size_t i = 0; i < num_samples; ++i) {
+      double status = data->get_y(i, 1);
+      if (status > max_status) {
+        max_status = status;
+      }
+    }
+    num_event_types = std::max((size_t) 1, (size_t) max_status);
+  }
+
+  // Initialize cause weights to equal if not set
+  if (cause_weights.empty()) {
+    cause_weights.resize(num_event_types, 1.0);
+  }
+
   // Sort data if extratrees and not memory saving mode
   if (splitrule == EXTRATREES && !memory_saving_splitting) {
     data->sort();
@@ -114,16 +147,16 @@ void ForestSurvival::initInternal() {
 }
 
 void ForestSurvival::growInternal() {
-  
+
   // If unique time points not set, use observed times
   if (unique_timepoints.empty()) {
     setUniqueTimepoints(std::vector<double>());
   }
-  
-  
+
   trees.reserve(num_trees);
   for (size_t i = 0; i < num_trees; ++i) {
-    trees.push_back(std::make_unique<TreeSurvival>(&unique_timepoints, &response_timepointIDs));
+    trees.push_back(std::make_unique<TreeSurvival>(&unique_timepoints, &response_timepointIDs,
+        num_event_types, &cause_weights));
   }
 }
 
@@ -132,20 +165,22 @@ void ForestSurvival::allocatePredictMemory() {
   size_t num_timepoints = unique_timepoints.size();
   if (predict_all) {
     predictions = std::vector<std::vector<std::vector<double>>>(num_prediction_samples,
-        std::vector<std::vector<double>>(num_timepoints, std::vector<double>(num_trees, 0)));
+        std::vector<std::vector<double>>(num_event_types * num_timepoints, std::vector<double>(num_trees, 0)));
   } else if (prediction_type == TERMINALNODES) {
     predictions = std::vector<std::vector<std::vector<double>>>(1,
         std::vector<std::vector<double>>(num_prediction_samples, std::vector<double>(num_trees, 0)));
   } else {
     predictions = std::vector<std::vector<std::vector<double>>>(1,
-        std::vector<std::vector<double>>(num_prediction_samples, std::vector<double>(num_timepoints, 0)));
+        std::vector<std::vector<double>>(num_prediction_samples, std::vector<double>(num_event_types * num_timepoints, 0)));
   }
 }
 
 void ForestSurvival::predictInternal(size_t sample_idx) {
-  // For each timepoint sum over trees
+  size_t num_timepoints = unique_timepoints.size();
+
+  // For each timepoint and event type sum over trees
   if (predict_all) {
-    for (size_t j = 0; j < unique_timepoints.size(); ++j) {
+    for (size_t j = 0; j < num_event_types * num_timepoints; ++j) {
       for (size_t k = 0; k < num_trees; ++k) {
         predictions[sample_idx][j][k] = getTreePrediction(k, sample_idx)[j];
       }
@@ -155,7 +190,7 @@ void ForestSurvival::predictInternal(size_t sample_idx) {
       predictions[0][sample_idx][k] = getTreePredictionTerminalNodeID(k, sample_idx);
     }
   } else {
-    for (size_t j = 0; j < unique_timepoints.size(); ++j) {
+    for (size_t j = 0; j < num_event_types * num_timepoints; ++j) {
       double sample_time_prediction = 0;
       for (size_t k = 0; k < num_trees; ++k) {
         sample_time_prediction += getTreePrediction(k, sample_idx)[j];
@@ -168,12 +203,13 @@ void ForestSurvival::predictInternal(size_t sample_idx) {
 void ForestSurvival::computePredictionErrorInternal() {
 
   size_t num_timepoints = unique_timepoints.size();
+  size_t total_chf_size = num_event_types * num_timepoints;
 
   // For each sample sum over trees where sample is OOB
   std::vector<size_t> samples_oob_count;
   samples_oob_count.resize(num_samples, 0);
   predictions = std::vector<std::vector<std::vector<double>>>(1,
-      std::vector<std::vector<double>>(num_samples, std::vector<double>(num_timepoints, 0)));
+      std::vector<std::vector<double>>(num_samples, std::vector<double>(total_chf_size, 0)));
 
   for (size_t tree_idx = 0; tree_idx < num_trees; ++tree_idx) {
     for (size_t sample_idx = 0; sample_idx < trees[tree_idx]->getNumSamplesOob(); ++sample_idx) {
@@ -187,17 +223,20 @@ void ForestSurvival::computePredictionErrorInternal() {
     }
   }
 
-  // Divide sample predictions by number of trees where sample is oob and compute summed chf for samples
+  // Divide sample predictions by number of trees where sample is oob and compute summed chf for first event type
   std::vector<double> sum_chf;
   sum_chf.reserve(predictions[0].size());
   std::vector<size_t> oob_sampleIDs;
   oob_sampleIDs.reserve(predictions[0].size());
   for (size_t i = 0; i < predictions[0].size(); ++i) {
     if (samples_oob_count[i] > 0) {
-      double sum = 0;
       for (size_t j = 0; j < predictions[0][i].size(); ++j) {
         predictions[0][i][j] /= samples_oob_count[i];
-        sum += predictions[0][i][j];
+      }
+      // Sum CHF for first event type only (indices 0..num_timepoints-1) for concordance
+      double sum = 0;
+      for (size_t t = 0; t < num_timepoints; ++t) {
+        sum += predictions[0][i][t];
       }
       sum_chf.push_back(sum);
       oob_sampleIDs.push_back(i);
diff --git a/src/ForestSurvival.h b/src/ForestSurvival.h
index d2efe4a82..76ea15913 100644
--- a/src/ForestSurvival.h
+++ b/src/ForestSurvival.h
@@ -34,15 +34,26 @@ class ForestSurvival: public Forest {
       std::vector<std::vector<size_t>>& forest_split_varIDs, std::vector<std::vector<double>>& forest_split_values,
       std::vector<std::vector<std::vector<double>> >& forest_chf, std::vector<double>& unique_timepoints,
       std::vector<bool>& is_ordered_variable);
-  
+
   void setUniqueTimepoints(const std::vector<double>& time_interest);
 
+  void setNumEventTypes(size_t num_event_types) {
+    this->num_event_types = num_event_types;
+  }
+  void setCauseWeights(const std::vector<double>& cause_weights) {
+    this->cause_weights = cause_weights;
+  }
+
   std::vector<std::vector<std::vector<double>>> getChf() const;
 
   const std::vector<double>& getUniqueTimepoints() const {
     return unique_timepoints;
   }
 
+  size_t getNumEventTypes() const {
+    return num_event_types;
+  }
+
 private:
   void initInternal() override;
   void growInternal() override;
@@ -57,6 +68,8 @@ class ForestSurvival: public Forest {
 
   std::vector<double> unique_timepoints;
   std::vector<size_t> response_timepointIDs;
+  size_t num_event_types = 1;
+  std::vector<double> cause_weights;
 
 private:
   const std::vector<double>& getTreePrediction(size_t tree_idx, size_t sample_idx) const;
diff --git a/src/TreeSurvival.cpp b/src/TreeSurvival.cpp
index 7ae947681..f25d2eab7 100644
--- a/src/TreeSurvival.cpp
+++ b/src/TreeSurvival.cpp
@@ -21,23 +21,26 @@
 
 namespace ranger {
 
-TreeSurvival::TreeSurvival(std::vector<double>* unique_timepoints, std::vector<size_t>* response_timepointIDs) :
-    unique_timepoints(unique_timepoints), response_timepointIDs(response_timepointIDs), num_deaths(0), num_samples_at_risk(
-        0) {
+TreeSurvival::TreeSurvival(std::vector<double>* unique_timepoints, std::vector<size_t>* response_timepointIDs,
+    size_t num_event_types, const std::vector<double>* cause_weights) :
+    unique_timepoints(unique_timepoints), response_timepointIDs(response_timepointIDs),
+    num_event_types(num_event_types), cause_weights(cause_weights),
+    num_samples_at_risk(0) {
   this->num_timepoints = unique_timepoints->size();
 }
 
 TreeSurvival::TreeSurvival(std::vector<std::vector<size_t>>& child_nodeIDs, std::vector<size_t>& split_varIDs,
     std::vector<double>& split_values, std::vector<std::vector<double>> chf, std::vector<double>* unique_timepoints,
-    std::vector<size_t>* response_timepointIDs) :
+    std::vector<size_t>* response_timepointIDs, size_t num_event_types) :
     Tree(child_nodeIDs, split_varIDs, split_values), unique_timepoints(unique_timepoints), response_timepointIDs(
-        response_timepointIDs), chf(chf), num_deaths(0), num_samples_at_risk(0) {
+        response_timepointIDs), chf(chf), num_event_types(num_event_types), cause_weights(nullptr),
+    num_samples_at_risk(0) {
   this->num_timepoints = unique_timepoints->size();
 }
 
 void TreeSurvival::allocateMemory() {
-  // Number of deaths and samples at risk for each timepoint
-  num_deaths.resize(num_timepoints);
+  // Number of deaths per event type and samples at risk for each timepoint
+  num_deaths.resize(num_event_types, std::vector<size_t>(num_timepoints, 0));
   num_samples_at_risk.resize(num_timepoints);
 }
 
@@ -62,24 +65,32 @@ void TreeSurvival::createEmptyNodeInternal() {
 
 void TreeSurvival::computeSurvival(size_t nodeID) {
   std::vector<double> chf_temp;
-  chf_temp.reserve(num_timepoints);
-  double chf_value = 0;
-  for (size_t i = 0; i < num_timepoints; ++i) {
-    if (num_samples_at_risk[i] != 0) {
-      chf_value += (double) num_deaths[i] / (double) num_samples_at_risk[i];
+  chf_temp.reserve(num_event_types * num_timepoints);
+  for (size_t e = 0; e < num_event_types; ++e) {
+    double chf_value = 0;
+    for (size_t t = 0; t < num_timepoints; ++t) {
+      if (num_samples_at_risk[t] != 0) {
+        chf_value += (double) num_deaths[e][t] / (double) num_samples_at_risk[t];
+      }
+      chf_temp.push_back(chf_value);
     }
-    chf_temp.push_back(chf_value);
   }
   chf[nodeID] = chf_temp;
 }
 
 double TreeSurvival::computePredictionAccuracyInternal(std::vector<double>* prediction_error_casewise) {
 
-  // Compute summed chf for samples
+  // Compute summed chf for samples (use first event type only for concordance)
   std::vector<double> sum_chf;
   for (size_t i = 0; i < prediction_terminal_nodeIDs.size(); ++i) {
     size_t terminal_nodeID = prediction_terminal_nodeIDs[i];
-    sum_chf.push_back(std::accumulate(chf[terminal_nodeID].begin(), chf[terminal_nodeID].end(), 0.0));
+    const auto& node_chf = chf[terminal_nodeID];
+    // Sum only the first event type's CHF (indices 0..num_timepoints-1)
+    double sum = 0;
+    for (size_t t = 0; t < num_timepoints; ++t) {
+      sum += node_chf[t];
+    }
+    sum_chf.push_back(sum);
   }
 
   // Return concordance index
@@ -331,9 +342,13 @@ bool TreeSurvival::findBestSplitMaxstat(size_t nodeID, std::vector<size_t>& poss
 void TreeSurvival::computeDeathCounts(size_t nodeID) {
 
   // Initialize
-  for (size_t i = 0; i < num_timepoints; ++i) {
-    num_deaths[i] = 0;
-    num_samples_at_risk[i] = 0;
+  for (size_t e = 0; e < num_event_types; ++e) {
+    for (size_t t = 0; t < num_timepoints; ++t) {
+      num_deaths[e][t] = 0;
+    }
+  }
+  for (size_t t = 0; t < num_timepoints; ++t) {
+    num_samples_at_risk[t] = 0;
   }
 
   for (size_t pos = start_pos[nodeID]; pos < end_pos[nodeID]; ++pos) {
@@ -346,11 +361,13 @@ void TreeSurvival::computeDeathCounts(size_t nodeID) {
       ++t;
     }
 
-    // Now t is the survival time, add to at risk and to death if death
+    // Now t is the survival time, add to at risk and to death if event
     if (t < num_timepoints) {
       ++num_samples_at_risk[t];
-      if (data->get_y(sampleID, 1) == 1) {
-        ++num_deaths[t];
+      double status = data->get_y(sampleID, 1);
+      if (status > 0) {
+        size_t event_type = (size_t) status - 1;
+        ++num_deaths[event_type][t];
       }
     }
   }
@@ -358,13 +375,14 @@ void TreeSurvival::computeDeathCounts(size_t nodeID) {
 
 void TreeSurvival::computeChildDeathCounts(size_t nodeID, size_t varID, std::vector<double>& possible_split_values,
     std::vector<size_t>& num_samples_right_child, std::vector<size_t>& delta_samples_at_risk_right_child,
-    std::vector<size_t>& num_deaths_right_child, size_t num_splits) {
+    std::vector<std::vector<size_t>>& num_deaths_right_child, size_t num_splits) {
 
-  // Count deaths in right child per timepoint and possbile split
+  // Count deaths in right child per timepoint and possible split
   for (size_t pos = start_pos[nodeID]; pos < end_pos[nodeID]; ++pos) {
     size_t sampleID = sampleIDs[pos];
     double value = data->get_x(sampleID, varID);
     size_t survival_timeID = (*response_timepointIDs)[sampleID];
+    double status = data->get_y(sampleID, 1);
 
     // Count deaths until split_value reached
     for (size_t i = 0; i < num_splits; ++i) {
@@ -372,8 +390,9 @@ void TreeSurvival::computeChildDeathCounts(size_t nodeID, size_t varID, std::vec
       if (value > possible_split_values[i]) {
         ++num_samples_right_child[i];
         ++delta_samples_at_risk_right_child[i * num_timepoints + survival_timeID];
-        if (data->get_y(sampleID, 1) == 1) {
-          ++num_deaths_right_child[i * num_timepoints + survival_timeID];
+        if (status > 0) {
+          size_t event_type = (size_t) status - 1;
+          ++num_deaths_right_child[event_type][i * num_timepoints + survival_timeID];
         }
       } else {
         break;
@@ -382,6 +401,97 @@ void TreeSurvival::computeChildDeathCounts(size_t nodeID, size_t varID, std::vec
   }
 }
 
+double TreeSurvival::computeLogRankSplit(size_t i, size_t num_samples_node,
+    std::vector<size_t>& num_samples_right_child,
+    std::vector<size_t>& delta_samples_at_risk_right_child,
+    std::vector<std::vector<size_t>>& num_deaths_right_child) {
+
+  // Stop if minimal bucket size reached
+  size_t num_samples_left_child = num_samples_node - num_samples_right_child[i];
+  if (num_samples_right_child[i] < (*min_bucket)[0] || num_samples_left_child < (*min_bucket)[0]) {
+    return -1;
+  }
+
+  // Composite log-rank across event types
+  double composite = 0;
+  for (size_t e = 0; e < num_event_types; ++e) {
+    double numerator = 0;
+    double denominator_squared = 0;
+
+    size_t num_samples_at_risk_right_child_e = num_samples_right_child[i];
+    for (size_t t = 0; t < num_timepoints; ++t) {
+      if (num_samples_at_risk[t] < 2 || num_samples_at_risk_right_child_e < 1) {
+        break;
+      }
+
+      if (num_deaths[e][t] > 0) {
+        double di = (double) num_deaths[e][t];
+        double di1 = (double) num_deaths_right_child[e][i * num_timepoints + t];
+        double Yi = (double) num_samples_at_risk[t];
+        double Yi1 = (double) num_samples_at_risk_right_child_e;
+        numerator += di1 - Yi1 * (di / Yi);
+        denominator_squared += (Yi1 / Yi) * (1.0 - Yi1 / Yi) * ((Yi - di) / (Yi - 1)) * di;
+      }
+
+      // Reduce number of samples at risk for next timepoint (shared across event types)
+      num_samples_at_risk_right_child_e -= delta_samples_at_risk_right_child[i * num_timepoints + t];
+    }
+
+    if (denominator_squared > 0) {
+      double z_e = numerator / sqrt(denominator_squared);
+      double w = (cause_weights != nullptr) ? (*cause_weights)[e] : 1.0;
+      composite += w * z_e * z_e;
+    }
+  }
+
+  return (composite > 0) ? sqrt(composite) : -1;
+}
+
+double TreeSurvival::computeLogRankUnorderedSplit(size_t num_samples_node,
+    size_t num_samples_right_child_count,
+    std::vector<size_t>& delta_samples_at_risk_right_child,
+    std::vector<size_t>& num_deaths_right_child_flat) {
+
+  // Stop if minimal bucket size reached
+  size_t num_samples_left_child = num_samples_node - num_samples_right_child_count;
+  if (num_samples_right_child_count < (*min_bucket)[0] || num_samples_left_child < (*min_bucket)[0]) {
+    return -1;
+  }
+
+  // For unordered splits, num_deaths_right_child_flat is flat: [e * num_timepoints + t]
+  double composite = 0;
+  for (size_t e = 0; e < num_event_types; ++e) {
+    double numerator = 0;
+    double denominator_squared = 0;
+
+    size_t num_samples_at_risk_right_child_e = num_samples_right_child_count;
+    for (size_t t = 0; t < num_timepoints; ++t) {
+      if (num_samples_at_risk[t] < 2 || num_samples_at_risk_right_child_e < 1) {
+        break;
+      }
+
+      if (num_deaths[e][t] > 0) {
+        double di = (double) num_deaths[e][t];
+        double di1 = (double) num_deaths_right_child_flat[e * num_timepoints + t];
+        double Yi = (double) num_samples_at_risk[t];
+        double Yi1 = (double) num_samples_at_risk_right_child_e;
+        numerator += di1 - Yi1 * (di / Yi);
+        denominator_squared += (Yi1 / Yi) * (1.0 - Yi1 / Yi) * ((Yi - di) / (Yi - 1)) * di;
+      }
+
+      num_samples_at_risk_right_child_e -= delta_samples_at_risk_right_child[t];
+    }
+
+    if (denominator_squared > 0) {
+      double z_e = numerator / sqrt(denominator_squared);
+      double w = (cause_weights != nullptr) ? (*cause_weights)[e] : 1.0;
+      composite += w * z_e * z_e;
+    }
+  }
+
+  return (composite > 0) ? sqrt(composite) : -1;
+}
+
 void TreeSurvival::findBestSplitValueLogRank(size_t nodeID, size_t varID, double& best_value, size_t& best_varID,
     double& best_logrank) {
 
@@ -400,7 +510,7 @@ void TreeSurvival::findBestSplitValueLogRank(size_t nodeID, size_t varID, double
   size_t num_splits = possible_split_values.size() - 1;
 
   // Initialize
-  std::vector<size_t> num_deaths_right_child(num_splits * num_timepoints);
+  std::vector<std::vector<size_t>> num_deaths_right_child(num_event_types, std::vector<size_t>(num_splits * num_timepoints, 0));
   std::vector<size_t> delta_samples_at_risk_right_child(num_splits * num_timepoints);
   std::vector<size_t> num_samples_right_child(num_splits);
 
@@ -409,40 +519,8 @@ void TreeSurvival::findBestSplitValueLogRank(size_t nodeID, size_t varID, double
 
   // Compute logrank test for all splits and use best
   for (size_t i = 0; i < num_splits; ++i) {
-    double numerator = 0;
-    double denominator_squared = 0;
-
-    // Stop if minimal bucket size reached
-    size_t num_samples_left_child = num_samples_node - num_samples_right_child[i];
-    if (num_samples_right_child[i] < (*min_bucket)[0] || num_samples_left_child < (*min_bucket)[0]) {
-      continue;
-    }
-
-    // Compute logrank test statistic for this split
-    size_t num_samples_at_risk_right_child = num_samples_right_child[i];
-    for (size_t t = 0; t < num_timepoints; ++t) {
-      if (num_samples_at_risk[t] < 2 || num_samples_at_risk_right_child < 1) {
-        break;
-      }
-
-      if (num_deaths[t] > 0) {
-        // Numerator and demoninator for log-rank test, notation from Ishwaran et al.
-        double di = (double) num_deaths[t];
-        double di1 = (double) num_deaths_right_child[i * num_timepoints + t];
-        double Yi = (double) num_samples_at_risk[t];
-        double Yi1 = (double) num_samples_at_risk_right_child;
-        numerator += di1 - Yi1 * (di / Yi);
-        denominator_squared += (Yi1 / Yi) * (1.0 - Yi1 / Yi) * ((Yi - di) / (Yi - 1)) * di;
-      }
-
-      // Reduce number of samples at risk for next timepoint
-      num_samples_at_risk_right_child -= delta_samples_at_risk_right_child[i * num_timepoints + t];
-
-    }
-    double logrank = -1;
-    if (denominator_squared != 0) {
-      logrank = fabs(numerator / sqrt(denominator_squared));
-    }
+    double logrank = computeLogRankSplit(i, num_samples_node, num_samples_right_child,
+        delta_samples_at_risk_right_child, num_deaths_right_child);
 
     // Regularization
     regularize(logrank, varID);
@@ -492,12 +570,10 @@ void TreeSurvival::findBestSplitValueLogRankUnordered(size_t nodeID, size_t varI
       }
     }
 
-    // Initialize
-    std::vector<size_t> num_deaths_right_child(num_timepoints);
-    std::vector<size_t> delta_samples_at_risk_right_child(num_timepoints);
+    // Initialize: flat layout [event_type * num_timepoints + t]
+    std::vector<size_t> num_deaths_right_child(num_event_types * num_timepoints, 0);
+    std::vector<size_t> delta_samples_at_risk_right_child(num_timepoints, 0);
     size_t num_samples_right_child = 0;
-    double numerator = 0;
-    double denominator_squared = 0;
 
     // Count deaths in right child per timepoint
     for (size_t pos = start_pos[nodeID]; pos < end_pos[nodeID]; ++pos) {
@@ -507,47 +583,19 @@ void TreeSurvival::findBestSplitValueLogRankUnordered(size_t nodeID, size_t varI
       size_t factorID = floor(value) - 1;
 
       // If in right child, count
-      // In right child, if bitwise splitID at position factorID is 1
       if ((splitID & (1ULL << factorID))) {
         ++num_samples_right_child;
         ++delta_samples_at_risk_right_child[survival_timeID];
-        if (data->get_y(sampleID, 1) == 1) {
-          ++num_deaths_right_child[survival_timeID];
+        double status = data->get_y(sampleID, 1);
+        if (status > 0) {
+          size_t event_type = (size_t) status - 1;
+          ++num_deaths_right_child[event_type * num_timepoints + survival_timeID];
         }
       }
-
     }
 
-    // Stop if minimal bucket size reached
-    size_t num_samples_left_child = num_samples_node - num_samples_right_child;
-    if (num_samples_right_child < (*min_bucket)[0] || num_samples_left_child < (*min_bucket)[0]) {
-      continue;
-    }
-
-    // Compute logrank test statistic for this split
-    size_t num_samples_at_risk_right_child = num_samples_right_child;
-    for (size_t t = 0; t < num_timepoints; ++t) {
-      if (num_samples_at_risk[t] < 2 || num_samples_at_risk_right_child < 1) {
-        break;
-      }
-
-      if (num_deaths[t] > 0) {
-        // Numerator and demoninator for log-rank test, notation from Ishwaran et al.
-        double di = (double) num_deaths[t];
-        double di1 = (double) num_deaths_right_child[t];
-        double Yi = (double) num_samples_at_risk[t];
-        double Yi1 = (double) num_samples_at_risk_right_child;
-        numerator += di1 - Yi1 * (di / Yi);
-        denominator_squared += (Yi1 / Yi) * (1.0 - Yi1 / Yi) * ((Yi - di) / (Yi - 1)) * di;
-      }
-
-      // Reduce number of samples at risk for next timepoint
-      num_samples_at_risk_right_child -= delta_samples_at_risk_right_child[t];
-    }
-    double logrank = -1;
-    if (denominator_squared != 0) {
-      logrank = fabs(numerator / sqrt(denominator_squared));
-    }
+    double logrank = computeLogRankUnorderedSplit(num_samples_node, num_samples_right_child,
+        delta_samples_at_risk_right_child, num_deaths_right_child);
 
     // Regularization
     regularize(logrank, varID);
@@ -791,7 +839,7 @@ void TreeSurvival::findBestSplitValueExtraTrees(size_t nodeID, size_t varID, dou
   size_t num_splits = possible_split_values.size();
 
   // Initialize
-  std::vector<size_t> num_deaths_right_child(num_splits * num_timepoints);
+  std::vector<std::vector<size_t>> num_deaths_right_child(num_event_types, std::vector<size_t>(num_splits * num_timepoints, 0));
   std::vector<size_t> delta_samples_at_risk_right_child(num_splits * num_timepoints);
   std::vector<size_t> num_samples_right_child(num_splits);
 
@@ -800,40 +848,8 @@ void TreeSurvival::findBestSplitValueExtraTrees(size_t nodeID, size_t varID, dou
 
   // Compute logrank test for all splits and use best
   for (size_t i = 0; i < num_splits; ++i) {
-    double numerator = 0;
-    double denominator_squared = 0;
-
-    // Stop if minimal node size reached
-    size_t num_samples_left_child = num_samples_node - num_samples_right_child[i];
-    if (num_samples_right_child[i] < (*min_bucket)[0] || num_samples_left_child < (*min_bucket)[0]) {
-      continue;
-    }
-
-    // Compute logrank test statistic for this split
-    size_t num_samples_at_risk_right_child = num_samples_right_child[i];
-    for (size_t t = 0; t < num_timepoints; ++t) {
-      if (num_samples_at_risk[t] < 2 || num_samples_at_risk_right_child < 1) {
-        break;
-      }
-
-      if (num_deaths[t] > 0) {
-        // Numerator and demoninator for log-rank test, notation from Ishwaran et al.
-        double di = (double) num_deaths[t];
-        double di1 = (double) num_deaths_right_child[i * num_timepoints + t];
-        double Yi = (double) num_samples_at_risk[t];
-        double Yi1 = (double) num_samples_at_risk_right_child;
-        numerator += di1 - Yi1 * (di / Yi);
-        denominator_squared += (Yi1 / Yi) * (1.0 - Yi1 / Yi) * ((Yi - di) / (Yi - 1)) * di;
-      }
-
-      // Reduce number of samples at risk for next timepoint
-      num_samples_at_risk_right_child -= delta_samples_at_risk_right_child[i * num_timepoints + t];
-
-    }
-    double logrank = -1;
-    if (denominator_squared != 0) {
-      logrank = fabs(numerator / sqrt(denominator_squared));
-    }
+    double logrank = computeLogRankSplit(i, num_samples_node, num_samples_right_child,
+        delta_samples_at_risk_right_child, num_deaths_right_child);
 
     // Regularization
     regularize(logrank, varID);
@@ -906,12 +922,10 @@ void TreeSurvival::findBestSplitValueExtraTreesUnordered(size_t nodeID, size_t v
       splitID |= 1ULL << idx;
     }
 
-    // Initialize
-    std::vector<size_t> num_deaths_right_child(num_timepoints);
-    std::vector<size_t> delta_samples_at_risk_right_child(num_timepoints);
+    // Initialize: flat layout [event_type * num_timepoints + t]
+    std::vector<size_t> num_deaths_right_child(num_event_types * num_timepoints, 0);
+    std::vector<size_t> delta_samples_at_risk_right_child(num_timepoints, 0);
     size_t num_samples_right_child = 0;
-    double numerator = 0;
-    double denominator_squared = 0;
 
     // Count deaths in right child per timepoint
     for (size_t pos = start_pos[nodeID]; pos < end_pos[nodeID]; ++pos) {
@@ -921,47 +935,19 @@ void TreeSurvival::findBestSplitValueExtraTreesUnordered(size_t nodeID, size_t v
       size_t factorID = floor(value) - 1;
 
       // If in right child, count
-      // In right child, if bitwise splitID at position factorID is 1
       if ((splitID & (1ULL << factorID))) {
         ++num_samples_right_child;
         ++delta_samples_at_risk_right_child[survival_timeID];
-        if (data->get_y(sampleID, 1) == 1) {
-          ++num_deaths_right_child[survival_timeID];
+        double status = data->get_y(sampleID, 1);
+        if (status > 0) {
+          size_t event_type = (size_t) status - 1;
+          ++num_deaths_right_child[event_type * num_timepoints + survival_timeID];
         }
       }
-
-    }
-
-    // Stop if minimal node size reached
-    size_t num_samples_left_child = num_samples_node - num_samples_right_child;
-    if (num_samples_right_child < (*min_bucket)[0] || num_samples_left_child < (*min_bucket)[0]) {
-      continue;
     }
 
-    // Compute logrank test statistic for this split
-    size_t num_samples_at_risk_right_child = num_samples_right_child;
-    for (size_t t = 0; t < num_timepoints; ++t) {
-      if (num_samples_at_risk[t] < 2 || num_samples_at_risk_right_child < 1) {
-        break;
-      }
-
-      if (num_deaths[t] > 0) {
-        // Numerator and demoninator for log-rank test, notation from Ishwaran et al.
-        double di = (double) num_deaths[t];
-        double di1 = (double) num_deaths_right_child[t];
-        double Yi = (double) num_samples_at_risk[t];
-        double Yi1 = (double) num_samples_at_risk_right_child;
-        numerator += di1 - Yi1 * (di / Yi);
-        denominator_squared += (Yi1 / Yi) * (1.0 - Yi1 / Yi) * ((Yi - di) / (Yi - 1)) * di;
-      }
-
-      // Reduce number of samples at risk for next timepoint
-      num_samples_at_risk_right_child -= delta_samples_at_risk_right_child[t];
-    }
-    double logrank = -1;
-    if (denominator_squared != 0) {
-      logrank = fabs(numerator / sqrt(denominator_squared));
-    }
+    double logrank = computeLogRankUnorderedSplit(num_samples_node, num_samples_right_child,
+        delta_samples_at_risk_right_child, num_deaths_right_child);
 
     // Regularization
     regularize(logrank, varID);
diff --git a/src/TreeSurvival.h b/src/TreeSurvival.h
index f27cef3c8..8dce54ddf 100644
--- a/src/TreeSurvival.h
+++ b/src/TreeSurvival.h
@@ -21,12 +21,13 @@ namespace ranger {
 
 class TreeSurvival: public Tree {
 public:
-  TreeSurvival(std::vector<double>* unique_timepoints, std::vector<size_t>* response_timepointIDs);
+  TreeSurvival(std::vector<double>* unique_timepoints, std::vector<size_t>* response_timepointIDs,
+      size_t num_event_types = 1, const std::vector<double>* cause_weights = nullptr);
 
   // Create from loaded forest
   TreeSurvival(std::vector<std::vector<size_t>>& child_nodeIDs, std::vector<size_t>& split_varIDs,
       std::vector<double>& split_values, std::vector<std::vector<double>> chf, std::vector<double>* unique_timepoints,
-      std::vector<size_t>* response_timepointIDs);
+      std::vector<size_t>* response_timepointIDs, size_t num_event_types = 1);
 
   TreeSurvival(const TreeSurvival&) = delete;
   TreeSurvival& operator=(const TreeSurvival&) = delete;
@@ -70,8 +71,17 @@ class TreeSurvival: public Tree {
 
   void computeDeathCounts(size_t nodeID);
   void computeChildDeathCounts(size_t nodeID, size_t varID, std::vector<double>& possible_split_values,
-      std::vector<size_t>& num_samples_right_child, std::vector<size_t>& num_samples_at_risk_right_child,
-      std::vector<size_t>& num_deaths_right_child, size_t num_splits);
+      std::vector<size_t>& num_samples_right_child, std::vector<size_t>& delta_samples_at_risk_right_child,
+      std::vector<std::vector<size_t>>& num_deaths_right_child, size_t num_splits);
+
+  double computeLogRankSplit(size_t i, size_t num_samples_node,
+      std::vector<size_t>& num_samples_right_child,
+      std::vector<size_t>& delta_samples_at_risk_right_child,
+      std::vector<std::vector<size_t>>& num_deaths_right_child);
+  double computeLogRankUnorderedSplit(size_t num_samples_node,
+      size_t num_samples_right_child,
+      std::vector<size_t>& delta_samples_at_risk_right_child,
+      std::vector<size_t>& num_deaths_right_child);
 
   void computeAucSplit(double time_k, double time_l, double status_k, double status_l, double value_k, double value_l,
       size_t num_splits, std::vector<double>& possible_split_values, std::vector<double>& num_count,
@@ -105,8 +115,13 @@ class TreeSurvival: public Tree {
   // For all terminal nodes CHF for all unique timepoints. For other nodes empty vector (except if save_node_stats).
   std::vector<std::vector<double>> chf;
 
+  // Number of competing event types (1 = standard survival)
+  size_t num_event_types;
+  const std::vector<double>* cause_weights;
+
   // Fields to save to while tree growing
-  std::vector<size_t> num_deaths;
+  // num_deaths[event_type][timepoint]
+  std::vector<std::vector<size_t>> num_deaths;
   std::vector<size_t> num_samples_at_risk;
 };
 
diff --git a/src/rangerCpp.cpp b/src/rangerCpp.cpp
index e990fd41b..37721f28d 100644
--- a/src/rangerCpp.cpp
+++ b/src/rangerCpp.cpp
@@ -177,6 +177,10 @@ Rcpp::List rangerCpp(uint treetype, Rcpp::NumericMatrix& input_x, Rcpp::NumericM
         std::vector<std::vector<std::vector<double>> > chf = loaded_forest["chf"];
         std::vector<double> unique_timepoints = loaded_forest["unique.death.times"];
         auto& temp = dynamic_cast<ForestSurvival&>(*forest);
+        if (loaded_forest.containsElementNamed("num.event.types")) {
+          size_t loaded_num_event_types = loaded_forest["num.event.types"];
+          temp.setNumEventTypes(loaded_num_event_types);
+        }
         temp.loadForest(num_trees, child_nodeIDs, split_varIDs, split_values, chf,
             unique_timepoints, is_ordered);
       } else if (treetype == TREE_PROBABILITY) {
@@ -232,6 +236,7 @@ Rcpp::List rangerCpp(uint treetype, Rcpp::NumericMatrix& input_x, Rcpp::NumericM
     if (treetype == TREE_SURVIVAL) {
       auto& temp = dynamic_cast<ForestSurvival&>(*forest);
       result.push_back(temp.getUniqueTimepoints(), "unique.death.times");
+      result.push_back(temp.getNumEventTypes(), "num.event.types");
     }
     if (!prediction_mode) {
       result.push_back(forest->getMtry(), "mtry");
@@ -287,6 +292,7 @@ Rcpp::List rangerCpp(uint treetype, Rcpp::NumericMatrix& input_x, Rcpp::NumericM
         auto& temp = dynamic_cast<ForestSurvival&>(*forest);
         forest_object.push_back(temp.getChf(), "chf");
         forest_object.push_back(temp.getUniqueTimepoints(), "unique.death.times");
+        forest_object.push_back(temp.getNumEventTypes(), "num.event.types");
       }
       result.push_back(forest_object, "forest");
     }
diff --git a/tests/testthat/test_competing_risks.R b/tests/testthat/test_competing_risks.R
new file mode 100644
index 000000000..43370f71a
--- /dev/null
+++ b/tests/testthat/test_competing_risks.R
@@ -0,0 +1,166 @@
+## Tests for competing risks survival forests
+
+library(ranger)
+library(survival)
+context("ranger_competing_risks")
+
+## Create test data with competing risks
+set.seed(42)
+n <- 200
+x1 <- rnorm(n)
+x2 <- rnorm(n)
+time <- rexp(n, rate = exp(0.3 * x1))
+event_type <- sample(0:2, n, replace = TRUE, prob = c(0.3, 0.4, 0.3))
+dat_cr <- data.frame(time = time, status = event_type, x1 = x1, x2 = x2)
+
+## Standard single-event survival data for backward compatibility
+dat_surv <- dat_cr
+dat_surv$status <- ifelse(dat_cr$status > 0, 1, 0)
+
+test_that("competing risks forest grows without error", {
+  rf <- ranger(Surv(time, status) ~ ., data = dat_cr, num.trees = 20)
+  expect_is(rf, "ranger")
+  expect_equal(rf$treetype, "Survival")
+})
+
+test_that("correct number of event types detected", {
+  rf <- ranger(Surv(time, status) ~ ., data = dat_cr, num.trees = 20)
+  expect_equal(rf$num.event.types, 2)
+})
+
+test_that("CHF output is list of matrices for competing risks", {
+  rf <- ranger(Surv(time, status) ~ ., data = dat_cr, num.trees = 20)
+  expect_is(rf$chf, "list")
+  expect_equal(length(rf$chf), 2)
+  expect_is(rf$chf[[1]], "matrix")
+  expect_is(rf$chf[[2]], "matrix")
+  expect_equal(nrow(rf$chf[[1]]), n)
+  expect_equal(ncol(rf$chf[[1]]), length(rf$unique.death.times))
+})
+
+test_that("CIF output is computed and valid", {
+  rf <- ranger(Surv(time, status) ~ ., data = dat_cr, num.trees = 20)
+  expect_is(rf$cif, "list")
+  expect_equal(length(rf$cif), 2)
+  expect_equal(nrow(rf$cif[[1]]), n)
+
+  # CIF values should be between 0 and 1
+  for (e in 1:2) {
+    expect_true(all(rf$cif[[e]] >= -1e-10))
+    expect_true(all(rf$cif[[e]] <= 1 + 1e-10))
+  }
+
+  # Sum of CIFs should be <= 1 at all timepoints
+  cif_sum <- rf$cif[[1]] + rf$cif[[2]]
+  expect_true(all(cif_sum <= 1 + 1e-10))
+})
+
+test_that("predictions work for new data with competing risks", {
+  rf <- ranger(Surv(time, status) ~ ., data = dat_cr, num.trees = 20)
+  pred <- predict(rf, data = dat_cr[1:10, ])
+  expect_is(pred$chf, "list")
+  expect_equal(length(pred$chf), 2)
+  expect_equal(nrow(pred$chf[[1]]), 10)
+  expect_is(pred$cif, "list")
+  expect_equal(length(pred$cif), 2)
+})
+
+test_that("backward compatibility: single-event survival identical behavior", {
+  rf <- ranger(Surv(time, status) ~ ., data = dat_surv, num.trees = 20, seed = 123)
+  expect_is(rf$chf, "matrix")
+  expect_is(rf$survival, "matrix")
+  expect_null(rf$cif)
+  expect_equal(rf$num.event.types, 1)
+})
+
+test_that("formula interface with Surv() works for competing risks", {
+  rf <- ranger(Surv(time, status) ~ x1 + x2, data = dat_cr, num.trees = 20)
+  expect_equal(rf$num.event.types, 2)
+})
+
+test_that("alternative interface works for competing risks", {
+  rf <- ranger(dependent.variable.name = "time", status.variable.name = "status",
+               data = dat_cr, num.trees = 20)
+  expect_equal(rf$num.event.types, 2)
+  expect_is(rf$chf, "list")
+})
+
+test_that("OOB error is computed for competing risks", {
+  rf <- ranger(Surv(time, status) ~ ., data = dat_cr, num.trees = 50)
+  expect_true(!is.null(rf$prediction.error))
+  expect_true(rf$prediction.error >= 0 && rf$prediction.error <= 1)
+})
+
+test_that("extratrees splitrule works with competing risks", {
+  rf <- ranger(Surv(time, status) ~ ., data = dat_cr, num.trees = 20,
+               splitrule = "extratrees")
+  expect_is(rf, "ranger")
+  expect_equal(rf$num.event.types, 2)
+})
+
+test_that("maxstat splitrule works with competing risks", {
+  rf <- ranger(Surv(time, status) ~ ., data = dat_cr, num.trees = 20,
+               splitrule = "maxstat")
+  expect_is(rf, "ranger")
+  expect_equal(rf$num.event.types, 2)
+})
+
+test_that("C-index splitrule errors for competing risks", {
+  expect_error(
+    ranger(Surv(time, status) ~ ., data = dat_cr, num.trees = 20, splitrule = "C"),
+    "C-index.*splitting not supported for competing risks"
+  )
+})
+
+test_that("variable importance works with competing risks", {
+  rf <- ranger(Surv(time, status) ~ ., data = dat_cr, num.trees = 50,
+               importance = "permutation")
+  expect_true(!is.null(rf$variable.importance))
+  expect_equal(length(rf$variable.importance), 2)
+})
+
+test_that("impurity importance works with competing risks", {
+  rf <- ranger(Surv(time, status) ~ ., data = dat_cr, num.trees = 50,
+               importance = "impurity")
+  expect_true(!is.null(rf$variable.importance))
+  expect_equal(length(rf$variable.importance), 2)
+})
+
+test_that("invalid status values give error", {
+  dat_bad <- dat_cr
+  dat_bad$status[1] <- -1
+  expect_error(
+    ranger(dependent.variable.name = "time", status.variable.name = "status",
+           data = dat_bad, num.trees = 20),
+    "consecutive integers"
+  )
+})
+
+test_that("three event types work", {
+  dat3 <- dat_cr
+  dat3$status <- sample(0:3, n, replace = TRUE, prob = c(0.25, 0.25, 0.25, 0.25))
+  rf <- ranger(Surv(time, status) ~ ., data = dat3, num.trees = 20)
+  expect_equal(rf$num.event.types, 3)
+  expect_equal(length(rf$chf), 3)
+  expect_equal(length(rf$cif), 3)
+})
+
+test_that("forest object stores num.event.types", {
+  rf <- ranger(Surv(time, status) ~ ., data = dat_cr, num.trees = 20)
+  expect_equal(rf$forest$num.event.types, 2)
+})
+
+test_that("print works for competing risks", {
+  rf <- ranger(Surv(time, status) ~ ., data = dat_cr, num.trees = 20)
+  expect_output(print(rf), "Number of event types")
+})
+
+test_that("CHF values are non-decreasing for competing risks", {
+  rf <- ranger(Surv(time, status) ~ ., data = dat_cr, num.trees = 50)
+  for (e in 1:2) {
+    for (i in 1:min(10, nrow(rf$chf[[e]]))) {
+      diffs <- diff(rf$chf[[e]][i, ])
+      expect_true(all(diffs >= -1e-10))
+    }
+  }
+})
diff --git a/tests/testthat/test_survival.R b/tests/testthat/test_survival.R
index 358a40960..8ecea7b86 100644
--- a/tests/testthat/test_survival.R
+++ b/tests/testthat/test_survival.R
@@ -8,9 +8,9 @@ context("ranger_surv")
 rg.surv <- ranger(Surv(time, status) ~ ., data = veteran, num.trees = 10)
 
 ## Basic tests (for all random forests equal)
-test_that("survival result is of class ranger with 18 elements", {
+test_that("survival result is of class ranger with 19 elements", {
   expect_is(rg.surv, "ranger")
-  expect_equal(length(rg.surv), 18)
+  expect_equal(length(rg.surv), 19)
 })
 
 test_that("results have right number of trees", {
@@ -119,10 +119,11 @@ test_that("Survival error without covariates", {
                "Error: No covariates found.")
 })
 
-test_that("Survival error for competing risk data", {
-  sobj <- Surv(veteran$time, factor(sample(1:3, nrow(veteran), replace = TRUE)))
-  expect_error(ranger(y = sobj, x = veteran[, 1:2], num.trees = 5), 
-               "Error: Competing risks not supported yet\\. Use status=1 for events and status=0 for censoring\\.")
+test_that("Competing risk data works with y/x interface", {
+  sobj <- Surv(veteran$time, factor(sample(0:2, nrow(veteran), replace = TRUE)))
+  rf <- ranger(y = sobj, x = veteran[, 1:2], num.trees = 5)
+  expect_is(rf, "ranger")
+  expect_equal(rf$treetype, "Survival")
 })
 
 test_that("Right unique time points without time.interest", {

From 4432bb37657849ff90a630aee5b5dc653d23b49a Mon Sep 17 00:00:00 2001
From: Lukas Burk <burk@leibniz-bips.de>
Date: Sat, 28 Feb 2026 02:14:26 +0100
Subject: [PATCH 2/3] Fix partial argument match warnings in tests and source

- R/utility.R: prob -> probs in quantile() call
- test_interface.R: dependent.variable -> dependent.variable.name
- test_ranger.R: $prediction -> $predictions

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
---
 R/utility.R                     |  2 +-
 tests/testthat/test_interface.R |  4 ++--
 tests/testthat/test_ranger.R    | 22 +++++++++++-----------
 3 files changed, 14 insertions(+), 14 deletions(-)

diff --git a/R/utility.R b/R/utility.R
index d4f3e9d91..7eb511e4b 100644
--- a/R/utility.R
+++ b/R/utility.R
@@ -106,7 +106,7 @@ largest.quantile <- function(formula) {
   
   ## Use median survival if available or largest quantile available in all strata if median not available
   max_quant <- max(aggregate(smry$surv ~ smry$strata, FUN = min)[, "smry$surv"])
-  quantiles <- quantile(fit, conf.int = FALSE, prob = min(0.5, 1 - max_quant))[, 1]
+  quantiles <- quantile(fit, conf.int = FALSE, probs = min(0.5, 1 - max_quant))[, 1]
   names(quantiles) <- gsub(".+=", "", names(quantiles))
   
   ## Return ordered levels
diff --git a/tests/testthat/test_interface.R b/tests/testthat/test_interface.R
index 940dc4b83..38da7edf8 100644
--- a/tests/testthat/test_interface.R
+++ b/tests/testthat/test_interface.R
@@ -60,13 +60,13 @@ test_that("Error if interaction of factor variable included", {
 test_that("Working if dependent variable has attributes other than names", {
   iris2 <- iris
   attr(iris2$Sepal.Width, "aaa") <- "bbb"
-  expect_silent(ranger(data = iris2, dependent.variable = "Sepal.Width"))
+  expect_silent(ranger(data = iris2, dependent.variable.name = "Sepal.Width"))
 })
 
 test_that("Working if dependent variable is matrix with one column", {
   iris2 <- iris
   iris2$Sepal.Width = scale(iris$Sepal.Width)
-  expect_silent(ranger(data = iris2, dependent.variable = "Sepal.Width"))
+  expect_silent(ranger(data = iris2, dependent.variable.name = "Sepal.Width"))
 })
 
 test_that("Same result with x/y interface, classification", {
diff --git a/tests/testthat/test_ranger.R b/tests/testthat/test_ranger.R
index 73f7a63ae..834533642 100644
--- a/tests/testthat/test_ranger.R
+++ b/tests/testthat/test_ranger.R
@@ -347,28 +347,28 @@ test_that("min.bucket creates nodes of correct size", {
   # Size 2
   rf <- ranger(Species ~ ., iris, num.trees = 5, replace = FALSE, 
                min.bucket = 2, keep.inbag = TRUE)
-  pred <- predict(rf, iris, type = "terminalNodes")$prediction
+  pred <- predict(rf, iris, type = "terminalNodes")$predictions
   inbag <- sapply(rf$inbag.counts, function(x) x == 1)
   smallest_node <- min(sapply(1:ncol(pred), function(i) {
     min(table(pred[inbag[, i], i]))
   }))
   expect_gte(smallest_node, 2)
-  
+
   # Size 10
-  rf <- ranger(Species ~ ., iris, num.trees = 5, replace = FALSE, 
+  rf <- ranger(Species ~ ., iris, num.trees = 5, replace = FALSE,
                min.bucket = 10, keep.inbag = TRUE)
-  pred <- predict(rf, iris, type = "terminalNodes")$prediction
+  pred <- predict(rf, iris, type = "terminalNodes")$predictions
   inbag <- sapply(rf$inbag.counts, function(x) x == 1)
   smallest_node <- min(sapply(1:ncol(pred), function(i) {
     min(table(pred[inbag[, i], i]))
   }))
   expect_gte(smallest_node, 10)
-  
+
   # Random size
   min.bucket <- round(runif(1, 1, 40))
-  rf <- ranger(Species ~ ., iris, num.trees = 5, replace = FALSE, 
+  rf <- ranger(Species ~ ., iris, num.trees = 5, replace = FALSE,
                min.bucket = min.bucket, keep.inbag = TRUE)
-  pred <- predict(rf, iris, type = "terminalNodes")$prediction
+  pred <- predict(rf, iris, type = "terminalNodes")$predictions
   inbag <- sapply(rf$inbag.counts, function(x) x == 1)
   smallest_node <- min(sapply(1:ncol(pred), function(i) {
     min(table(pred[inbag[, i], i]))
@@ -381,7 +381,7 @@ test_that("Vector min.bucket creates nodes of correct size", {
   # Size 2,3,4
   rf <- ranger(Species ~ ., iris, num.trees = 5, replace = FALSE, 
                min.bucket = c(2, 3, 4), keep.inbag = TRUE)
-  pred <- predict(rf, iris, type = "terminalNodes")$prediction
+  pred <- predict(rf, iris, type = "terminalNodes")$predictions
   inbag <- sapply(rf$inbag.counts, function(x) x == 1)
   
   smallest_nodes <- sapply(1:ncol(pred), function(i) {
@@ -400,7 +400,7 @@ test_that("Vector min.bucket creates nodes of correct size", {
   # Size 4,3,2
   rf <- ranger(Species ~ ., iris, num.trees = 5, replace = FALSE, 
                min.bucket = c(4, 3, 2), keep.inbag = TRUE)
-  pred <- predict(rf, iris, type = "terminalNodes")$prediction
+  pred <- predict(rf, iris, type = "terminalNodes")$predictions
   inbag <- sapply(rf$inbag.counts, function(x) x == 1)
   
   smallest_nodes <- sapply(1:ncol(pred), function(i) {
@@ -420,7 +420,7 @@ test_that("Vector min.bucket creates nodes of correct size", {
   min.bucket <- round(runif(3, 1, 10))
   rf <- ranger(Species ~ ., iris, num.trees = 5, replace = FALSE, 
                min.bucket = min.bucket, keep.inbag = TRUE)
-  pred <- predict(rf, iris, type = "terminalNodes")$prediction
+  pred <- predict(rf, iris, type = "terminalNodes")$predictions
   inbag <- sapply(rf$inbag.counts, function(x) x == 1)
   
   smallest_nodes <- sapply(1:ncol(pred), function(i) {
@@ -439,7 +439,7 @@ test_that("Vector min.bucket creates nodes of correct size", {
   # No factor outcome
   rf <- ranger(Species ~ ., data.matrix(iris), num.trees = 5, replace = FALSE, 
                min.bucket = c(2, 3, 4), keep.inbag = TRUE, classification = TRUE)
-  pred <- predict(rf, iris, type = "terminalNodes")$prediction
+  pred <- predict(rf, iris, type = "terminalNodes")$predictions
   inbag <- sapply(rf$inbag.counts, function(x) x == 1)
   
   smallest_nodes <- sapply(1:ncol(pred), function(i) {

From eb014c1e857b9a031c9e6378f9418373edb12182 Mon Sep 17 00:00:00 2001
From: Lukas Burk <burk@leibniz-bips.de>
Date: Sat, 28 Feb 2026 14:57:19 +0100
Subject: [PATCH 3/3] Add competing risks documentation and vignette

- Update roxygen docs for ranger() and predict.ranger() with competing
  risks parameters, return values, and examples
- Add Ishwaran et al. (2014) reference for cause-specific hazard approach
- Add vignettes/competing_risks.Rmd with method description and usage
- Add inst/examples/competing_risks_comparison.R (ranger vs rfsrc)
- Add knitr/rmarkdown to Suggests and VignetteBuilder to DESCRIPTION

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
---
 DESCRIPTION                                |   7 +-
 R/predict.R                                |  10 +-
 R/ranger.R                                 |  32 +++-
 inst/examples/competing_risks_comparison.R | 208 +++++++++++++++++++++
 man/predict.ranger.Rd                      |   5 +-
 man/predict.ranger.forest.Rd               |   5 +-
 man/ranger.Rd                              |  30 ++-
 vignettes/competing_risks.Rmd              | 208 +++++++++++++++++++++
 8 files changed, 482 insertions(+), 23 deletions(-)
 create mode 100644 inst/examples/competing_risks_comparison.R
 create mode 100644 vignettes/competing_risks.Rmd

diff --git a/DESCRIPTION b/DESCRIPTION
index 452b453bd..1d9be39b6 100644
--- a/DESCRIPTION
+++ b/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/,
diff --git a/R/predict.R b/R/predict.R
index 9f1aad22d..d03222250 100644
--- a/R/predict.R
+++ b/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
@@ -484,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
diff --git a/R/ranger.R b/R/ranger.R
index dfebc6cc4..7b1d34419 100644
--- a/R/ranger.R
+++ b/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}.
diff --git a/inst/examples/competing_risks_comparison.R b/inst/examples/competing_risks_comparison.R
new file mode 100644
index 000000000..c97b3fd3b
--- /dev/null
+++ b/inst/examples/competing_risks_comparison.R
@@ -0,0 +1,208 @@
+## Competing risks: ranger vs randomForestSRC comparison
+##
+## Both packages are fit with cause-specific log-rank splitting
+## (splitrule = "logrank") and otherwise similar settings.
+## The resulting CIF curves should look similar for a large number of trees.
+
+library(ranger)
+library(randomForestSRC)
+library(survival)
+
+## ---------------------------------------------------------------------------
+## 1. Simulate competing risks data
+## ---------------------------------------------------------------------------
+set.seed(42)
+n <- 800
+
+x1 <- rnorm(n)
+x2 <- rnorm(n)
+x3 <- rbinom(n, 1, 0.5)
+
+## Latent event times for two competing causes
+t1 <- rexp(n, rate = exp(0.5 * x1 + 0.2 * x3))   # cause 1
+t2 <- rexp(n, rate = exp(0.3 * x2 - 0.1 * x3))   # cause 2
+tc <- rexp(n, rate = 0.3)                          # censoring
+
+## Observed time and status
+time   <- pmin(t1, t2, tc)
+status <- ifelse(time == tc, 0L,
+          ifelse(time == t1, 1L, 2L))
+
+dat <- data.frame(time = time, status = status, x1 = x1, x2 = x2, x3 = x3)
+
+cat("Event distribution:\n")
+print(table(status = dat$status))
+
+## ---------------------------------------------------------------------------
+## 2. Fit forests with matched settings
+## ---------------------------------------------------------------------------
+
+## Common settings
+ntree    <- 1000
+nodesize <- 15
+mtry     <- 2
+
+## ranger: uses composite cause-specific log-rank by default
+## replace = TRUE is the default (bootstrap with replacement)
+t_ranger <- system.time(
+  rf_ranger <- ranger(
+    Surv(time, status) ~ .,
+    data      = dat,
+    num.trees = ntree,
+    min.node.size = nodesize,
+    mtry      = mtry,
+    replace   = TRUE,
+    seed      = 1
+  )
+)
+
+## randomForestSRC: use splitrule = "logrank" for cause-specific log-rank
+## (the default "logrankCR" uses Gray's test, which differs)
+## Match ranger defaults: bootstrap w/ replacement, all split points, all timepoints
+t_rfsrc <- system.time(
+  rf_rfsrc <- rfsrc(
+    Surv(time, status) ~ .,
+    data      = dat,
+    ntree     = ntree,
+    nodesize  = nodesize,
+    mtry      = mtry,
+    splitrule = "logrank",
+    nsplit    = 0,
+    samptype  = "swr",
+    ntime     = 0,
+    seed      = -1
+  )
+)
+
+cat("\nTiming:\n")
+cat("  ranger:          ", round(t_ranger["elapsed"], 2), "s\n")
+cat("  randomForestSRC: ", round(t_rfsrc["elapsed"], 2), "s\n")
+
+## ---------------------------------------------------------------------------
+## 3. Extract CIF predictions (OOB)
+## ---------------------------------------------------------------------------
+
+## ranger: cif is a list of K matrices (n x num_timepoints)
+cif_ranger_1 <- rf_ranger$cif[[1]]
+cif_ranger_2 <- rf_ranger$cif[[2]]
+times_ranger <- rf_ranger$unique.death.times
+
+## randomForestSRC: cif.oob is a 3d array (n x num_timepoints x K)
+cif_rfsrc_1 <- rf_rfsrc$cif.oob[, , 1]
+cif_rfsrc_2 <- rf_rfsrc$cif.oob[, , 2]
+times_rfsrc <- rf_rfsrc$time.interest
+
+## ---------------------------------------------------------------------------
+## 4. Compare mean CIF curves across all subjects
+## ---------------------------------------------------------------------------
+
+mean_cif_ranger_1 <- colMeans(cif_ranger_1)
+mean_cif_ranger_2 <- colMeans(cif_ranger_2)
+
+mean_cif_rfsrc_1 <- colMeans(cif_rfsrc_1)
+mean_cif_rfsrc_2 <- colMeans(cif_rfsrc_2)
+
+## ---------------------------------------------------------------------------
+## 5. Plot mean CIF curves
+## ---------------------------------------------------------------------------
+
+pdf("competing_risks_comparison.pdf", width = 12, height = 10)
+
+par(mfrow = c(2, 2), mar = c(4.5, 4.5, 3, 1))
+
+## Cause 1
+plot(times_ranger, mean_cif_ranger_1, type = "s", col = "steelblue", lwd = 2,
+     xlab = "Time", ylab = "Mean CIF", main = "Cause 1 -Mean CIF",
+     ylim = c(0, max(mean_cif_ranger_1, mean_cif_rfsrc_1) * 1.1))
+lines(times_rfsrc, mean_cif_rfsrc_1, type = "s", col = "firebrick", lwd = 2, lty = 2)
+legend("bottomright", legend = c("ranger", "randomForestSRC"),
+       col = c("steelblue", "firebrick"), lwd = 2, lty = c(1, 2), bty = "n")
+
+## Cause 2
+plot(times_ranger, mean_cif_ranger_2, type = "s", col = "steelblue", lwd = 2,
+     xlab = "Time", ylab = "Mean CIF", main = "Cause 2 -Mean CIF",
+     ylim = c(0, max(mean_cif_ranger_2, mean_cif_rfsrc_2) * 1.1))
+lines(times_rfsrc, mean_cif_rfsrc_2, type = "s", col = "firebrick", lwd = 2, lty = 2)
+legend("bottomright", legend = c("ranger", "randomForestSRC"),
+       col = c("steelblue", "firebrick"), lwd = 2, lty = c(1, 2), bty = "n")
+
+## ---------------------------------------------------------------------------
+## 6. Compare CIF for individual subjects (high/low risk for cause 1)
+## ---------------------------------------------------------------------------
+
+## Pick a high-risk (large x1) and low-risk (small x1) subject
+i_high <- which.max(dat$x1)
+i_low  <- which.min(dat$x1)
+
+## High-risk subject
+plot(times_ranger, cif_ranger_1[i_high, ], type = "s", col = "steelblue", lwd = 2,
+     xlab = "Time", ylab = "CIF (Cause 1)",
+     main = sprintf("High-risk subject (x1 = %.1f)", dat$x1[i_high]),
+     ylim = c(0, 1))
+lines(times_rfsrc, cif_rfsrc_1[i_high, ], type = "s", col = "firebrick", lwd = 2, lty = 2)
+legend("bottomright", legend = c("ranger", "randomForestSRC"),
+       col = c("steelblue", "firebrick"), lwd = 2, lty = c(1, 2), bty = "n")
+
+## Low-risk subject
+plot(times_ranger, cif_ranger_1[i_low, ], type = "s", col = "steelblue", lwd = 2,
+     xlab = "Time", ylab = "CIF (Cause 1)",
+     main = sprintf("Low-risk subject (x1 = %.1f)", dat$x1[i_low]),
+     ylim = c(0, 1))
+lines(times_rfsrc, cif_rfsrc_1[i_low, ], type = "s", col = "firebrick", lwd = 2, lty = 2)
+legend("bottomright", legend = c("ranger", "randomForestSRC"),
+       col = c("steelblue", "firebrick"), lwd = 2, lty = c(1, 2), bty = "n")
+
+## ---------------------------------------------------------------------------
+## 7. Prediction on new data
+## ---------------------------------------------------------------------------
+
+newdat <- data.frame(
+  x1 = c(-1, 0, 1),
+  x2 = c( 0, 0, 0),
+  x3 = c( 0, 1, 1)
+)
+
+pred_ranger <- predict(rf_ranger, data = newdat)
+pred_rfsrc  <- predict(rf_rfsrc, newdata = newdat)
+
+par(mfrow = c(1, 3), mar = c(4.5, 4.5, 3, 1))
+for (i in 1:3) {
+  plot(pred_ranger$unique.death.times, pred_ranger$cif[[1]][i, ],
+       type = "s", col = "steelblue", lwd = 2,
+       xlab = "Time", ylab = "CIF (Cause 1)",
+       main = sprintf("x1=%.0f, x2=%.0f, x3=%.0f", newdat$x1[i], newdat$x2[i], newdat$x3[i]),
+       ylim = c(0, 1))
+  lines(pred_rfsrc$time.interest, pred_rfsrc$cif[i, , 1],
+        type = "s", col = "firebrick", lwd = 2, lty = 2)
+  legend("bottomright", legend = c("ranger", "rfsrc"),
+         col = c("steelblue", "firebrick"), lwd = 2, lty = c(1, 2), bty = "n")
+}
+
+dev.off()
+
+## ---------------------------------------------------------------------------
+## 8. Numeric comparison
+## ---------------------------------------------------------------------------
+
+## Interpolate rfsrc CIF onto ranger's time grid for direct comparison
+interp <- function(x_from, y_from, x_to) {
+  stepfun(x_from, c(0, y_from))(x_to)
+}
+
+cif1_rfsrc_interp <- t(sapply(1:n, function(i) {
+  interp(times_rfsrc, cif_rfsrc_1[i, ], times_ranger)
+}))
+
+## Per-subject mean absolute difference
+mad_per_subject <- rowMeans(abs(cif_ranger_1 - cif1_rfsrc_interp))
+
+cat("\nCause-1 CIF comparison (ranger vs rfsrc):\n")
+cat("  Per-subject mean absolute difference:\n")
+cat("    median: ", round(median(mad_per_subject), 4), "\n")
+cat("    95th:   ", round(quantile(mad_per_subject, 0.95), 4), "\n")
+cat("    max:    ", round(max(mad_per_subject), 4), "\n")
+cat("  Correlation of subject-level CIF at final timepoint:\n")
+cat("   ", round(cor(cif_ranger_1[, ncol(cif_ranger_1)],
+                     cif1_rfsrc_interp[, ncol(cif1_rfsrc_interp)]), 4), "\n")
+
+cat("\nPlots saved to competing_risks_comparison.pdf\n")
diff --git a/man/predict.ranger.Rd b/man/predict.ranger.Rd
index 2f9c63ac3..76c58f526 100644
--- a/man/predict.ranger.Rd
+++ b/man/predict.ranger.Rd
@@ -49,8 +49,9 @@ Object of class \code{ranger.prediction} with elements
   \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
diff --git a/man/predict.ranger.forest.Rd b/man/predict.ranger.forest.Rd
index 805effda4..a0bc04b98 100644
--- a/man/predict.ranger.forest.Rd
+++ b/man/predict.ranger.forest.Rd
@@ -46,8 +46,9 @@ Object of class \code{ranger.prediction} with elements
   \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
diff --git a/man/ranger.Rd b/man/ranger.Rd
index 66b7abd29..cf1a647c8 100644
--- a/man/ranger.Rd
+++ b/man/ranger.Rd
@@ -82,7 +82,7 @@ ranger(
 
 \item{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.}
 
 \item{num.random.splits}{For "extratrees" splitrule.: Number of random splits to consider for each candidate splitting variable.}
 
@@ -133,13 +133,13 @@ If a terminal node has only 0 responses, the estimate is set to \eqn{\alpha 0 +
 
 \item{dependent.variable.name}{Name of dependent variable, needed if no formula given. For survival forests this is the time variable.}
 
-\item{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.}
+\item{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.}
 
 \item{classification}{Set to \code{TRUE} to grow a classification forest. Only needed if the data is a matrix or the response numeric.}
 
 \item{x}{Predictor data (independent variables), alternative interface to data with formula or dependent.variable.name.}
 
-\item{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.}
+\item{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.}
 
 \item{...}{Further arguments passed to or from other methods (currently ignored).}
 }
@@ -153,8 +153,10 @@ Object of class \code{ranger} with elements
   \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.}
@@ -181,6 +183,12 @@ For regression, the estimated response variances or maximally selected rank stat
 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.
@@ -274,6 +282,15 @@ require(survival)
 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])
@@ -296,7 +313,8 @@ ranger(trait ~ ., data = dat.gwaa)
   \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}.
diff --git a/vignettes/competing_risks.Rmd b/vignettes/competing_risks.Rmd
new file mode 100644
index 000000000..61c27d9e1
--- /dev/null
+++ b/vignettes/competing_risks.Rmd
@@ -0,0 +1,208 @@
+---
+title: "Competing Risks Survival Forests with ranger"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Competing Risks Survival Forests with ranger}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r setup, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  fig.width = 7,
+  fig.height = 5
+)
+```
+
+## Overview
+
+Ranger supports competing risks survival forests, where subjects can experience one of several mutually exclusive event types.
+The implementation follows the cause-specific hazard approach described in Ishwaran et al. (2014), similar to the `randomForestSRC` package.
+
+In a competing risks setting, the status variable takes values:
+
+- `0` = censored (no event observed)
+- `1, 2, ..., K` = event of type 1, 2, ..., K
+
+For example, in a study of cancer patients, death from cancer (status = 1) and death from other causes (status = 2) are competing events.
+
+## Method
+
+**Splitting:** At each node, ranger computes a cause-specific log-rank statistic for each event type, treating all other events as censored.
+These are combined into a composite statistic: $\sqrt{\sum_{j=1}^K w_j Z_j^2}$, where $Z_j$ is the log-rank Z-statistic for event type $j$ and $w_j$ are equal cause weights.
+This composite test finds splits that are informative for discriminating risk across all event types simultaneously.
+
+**Predictions:** Each terminal node stores a cause-specific cumulative hazard function (CHF) per event type via the Nelson-Aalen estimator.
+The cumulative incidence functions (CIF) are then computed from the cause-specific CHFs using the Aalen-Johansen estimator:
+$$\text{CIF}_j(t) = \sum_{s \le t} \hat{S}(s^{-}) \cdot \hat{h}_j(s)$$
+where $\hat{S}(t^{-})$ is the event-free survival probability just before time $t$, and $\hat{h}_j(s)$ is the cause-specific hazard increment for event $j$.
+
+## Quick start
+
+```{r basic}
+library(ranger)
+library(survival)
+
+## Simulate competing risks data
+set.seed(42)
+n <- 500
+x1 <- rnorm(n)
+x2 <- rnorm(n)
+
+## Latent event times for two competing causes
+t1 <- rexp(n, rate = exp(0.5 * x1))   # cause 1 depends on x1
+t2 <- rexp(n, rate = exp(0.3 * x2))   # cause 2 depends on x2
+tc <- rexp(n, rate = 0.3)             # censoring
+
+time   <- pmin(t1, t2, tc)
+status <- ifelse(time == tc, 0L, ifelse(time == t1, 1L, 2L))
+
+dat <- data.frame(time = time, status = status, x1 = x1, x2 = x2)
+table(dat$status)
+```
+
+Fit a competing risks forest using the standard `Surv()` formula interface:
+
+```{r fit}
+rf <- ranger(Surv(time, status) ~ ., data = dat, num.trees = 500)
+rf
+```
+
+The `num.event.types` field confirms competing risks were detected:
+
+```{r check}
+rf$num.event.types
+```
+
+## Output structure
+
+For competing risks, `chf` and `cif` are lists of matrices (one per event type), each of dimension (n_samples x n_timepoints):
+
+```{r output}
+## Cause-specific CHFs
+length(rf$chf)
+dim(rf$chf[[1]])
+
+## Cumulative incidence functions
+length(rf$cif)
+dim(rf$cif[[1]])
+```
+
+Access the CIF for event type 1 for the first sample:
+
+```{r access}
+head(rf$cif[[1]][1, ])
+```
+
+## Plotting CIF curves
+
+```{r plot_mean, fig.cap = "Mean CIF across all subjects for both event types."}
+times <- rf$unique.death.times
+
+par(mar = c(4.5, 4.5, 2, 1))
+plot(times, colMeans(rf$cif[[1]]), type = "s", col = "steelblue", lwd = 2,
+     xlab = "Time", ylab = "Cumulative Incidence", ylim = c(0, 0.7))
+lines(times, colMeans(rf$cif[[2]]), type = "s", col = "firebrick", lwd = 2)
+legend("bottomright", legend = c("Event 1", "Event 2"),
+       col = c("steelblue", "firebrick"), lwd = 2, bty = "n")
+```
+
+Individual subject CIF curves for a high-risk and low-risk subject:
+
+```{r plot_individual, fig.cap = "CIF curves for individual subjects."}
+i_high <- which.max(dat$x1)
+i_low  <- which.min(dat$x1)
+
+par(mfrow = c(1, 2), mar = c(4.5, 4.5, 3, 1))
+
+plot(times, rf$cif[[1]][i_high, ], type = "s", col = "steelblue", lwd = 2,
+     xlab = "Time", ylab = "CIF (Event 1)",
+     main = sprintf("High risk (x1 = %.1f)", dat$x1[i_high]), ylim = c(0, 1))
+
+plot(times, rf$cif[[1]][i_low, ], type = "s", col = "steelblue", lwd = 2,
+     xlab = "Time", ylab = "CIF (Event 1)",
+     main = sprintf("Low risk (x1 = %.1f)", dat$x1[i_low]), ylim = c(0, 1))
+```
+
+## Prediction on new data
+
+```{r predict}
+newdat <- data.frame(x1 = c(-1, 0, 1), x2 = c(0, 0, 0))
+pred <- predict(rf, data = newdat)
+
+## Predictions have the same structure
+length(pred$cif)
+dim(pred$cif[[1]])
+```
+
+```{r plot_pred, fig.cap = "Predicted CIF for event 1 for three covariate profiles."}
+par(mar = c(4.5, 4.5, 2, 1))
+cols <- c("forestgreen", "orange", "firebrick")
+plot(NULL, xlim = range(pred$unique.death.times), ylim = c(0, 1),
+     xlab = "Time", ylab = "CIF (Event 1)")
+for (i in 1:3) {
+  lines(pred$unique.death.times, pred$cif[[1]][i, ],
+        type = "s", col = cols[i], lwd = 2)
+}
+legend("bottomright", legend = paste0("x1 = ", c(-1, 0, 1)),
+       col = cols, lwd = 2, bty = "n")
+```
+
+## Alternative interface
+
+The `dependent.variable.name` / `status.variable.name` interface also works:
+
+```{r alternative}
+rf2 <- ranger(dependent.variable.name = "time",
+              status.variable.name = "status",
+              data = dat, num.trees = 100)
+rf2$num.event.types
+```
+
+## Variable importance
+
+Variable importance works as usual:
+
+```{r importance}
+rf_imp <- ranger(Surv(time, status) ~ ., data = dat,
+                 num.trees = 500, importance = "permutation")
+rf_imp$variable.importance
+```
+
+## Split rules
+
+The following split rules are supported for competing risks:
+
+- `"logrank"` (default) -- composite cause-specific log-rank
+- `"extratrees"` -- extremely randomized trees with composite log-rank
+- `"maxstat"` -- maximally selected rank statistics
+
+The C-index (`"C"`) split rule is not supported for competing risks.
+
+```{r splitrules}
+rf_et <- ranger(Surv(time, status) ~ ., data = dat,
+                num.trees = 100, splitrule = "extratrees")
+rf_et$num.event.types
+```
+
+## Backward compatibility
+
+For standard survival data (binary status 0/1), ranger behaves exactly as before -- `chf` and `survival` are matrices and `cif` is `NULL`:
+
+```{r compat}
+dat_std <- dat
+dat_std$status <- ifelse(dat$status > 0, 1, 0)
+
+rf_std <- ranger(Surv(time, status) ~ ., data = dat_std, num.trees = 100)
+is.matrix(rf_std$chf)
+is.matrix(rf_std$survival)
+is.null(rf_std$cif)
+```
+
+## References
+
+- 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(4):757-773. doi:10.1093/biostatistics/kxu010.
+- Ishwaran, H., Kogalur, U. B., Blackstone, E. H. & Lauer, M. S. (2008). Random survival forests. *Annals of Applied Statistics*, 2(3):841-860.
+- Wright, M. N. & Ziegler, A. (2017). ranger: A fast implementation of random forests for high dimensional data in C++ and R. *Journal of Statistical Software*, 77(1):1-17. doi:10.18637/jss.v077.i01.