[c++] Add survival_cox objective for Cox proportional hazards modelling
#7212
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,51 @@ | ||
| # coding: utf-8 | ||
| import numpy as np | ||
| from sklearn.datasets import fetch_openml | ||
| from sklearn.model_selection import train_test_split | ||
|
|
||
| import lightgbm as lgb | ||
|
|
||
| # Load FLCHAIN dataset (serum free light chain and mortality) | ||
| data = fetch_openml("flchain", version=1, as_frame=True, parser="auto") | ||
|
Member
There was a problem hiding this comment. The AppVeyor builds are failing like this:
https://ci.appveyor.com/project/guolinke/lightgbm/builds/53791302/job/oj3cfvbuifsjc7au?fullLog=true Those jobs use a very old
Author
There was a problem hiding this comment. Thanks for pointing this out. I wasn't sure how to debug the appveyor failing tests.
Author
There was a problem hiding this comment. AppVeyor fails again with
I will look into it. It seems that
Member
There was a problem hiding this comment. In general, we try to support a wide range of We'd prefer to have a compelling reason to bump a runtime floor, and "makes this example in documentation easier to test" isn't that compelling, in my opinion. That said we do already have a Linux job testing an even older So I wouldn't be opposed to updating the pin for Python 3.9 environments like the one on Appveyor. That could be done here: LightGBM/.ci/conda-envs/ci-core-py39.txt Line 31 in a7d00a9 I'd support trying to bump that up to a newer
Author
There was a problem hiding this comment. I tested this on my local machine with These are the api calls that are made (generated by adding a print statement here) In any case I can just replace the example to use synthetic data.
Member
There was a problem hiding this comment.
That job uses Python 3.9, not 3.10. It's the standard Appveyor runner for open source projects and should have full access to the internet. I suspect that maybe that "not found" error is actually from a broad try-catch and that something else in the environment (like some other dependency version) is causing it to fail. The approach with synthetic data looks good to me! Thanks for working through that. |
||
| df = data.data.copy() | ||
| time = data.target.values | ||
| event = df.pop("status").values | ||
|
|
||
| # Encode strings as integers | ||
| df["sex"] = (df["sex"] == "F").astype(int) | ||
| df["flc.grp"] = df["flc.grp"].astype(int) | ||
| df["mgus"] = df["mgus"].astype(int) | ||
|
|
||
| # Encode labels: positive = event (death), negative = censored | ||
| y = np.where(event == 1, time, -time) | ||
| X = df.values | ||
|
|
||
| X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.2, random_state=42) | ||
|
|
||
| lgb_train = lgb.Dataset(X_train, label=y_train) | ||
| lgb_val = lgb.Dataset(X_val, label=y_val, reference=lgb_train) | ||
|
|
||
| params = { | ||
| "objective": "survival_cox", | ||
| "metric": ["survival_cox_nll", "concordance_index"], | ||
| "num_leaves": 31, | ||
| "learning_rate": 0.05, | ||
| "verbose": 0, | ||
| } | ||
|
|
||
| evals_result = {} | ||
| gbm = lgb.train( | ||
| params, | ||
| lgb_train, | ||
| num_boost_round=200, | ||
| valid_sets=[lgb_val], | ||
| valid_names=["val"], | ||
| callbacks=[ | ||
| lgb.early_stopping(stopping_rounds=20, first_metric_only=True), | ||
| lgb.record_evaluation(evals_result), | ||
| ], | ||
| ) | ||
|
|
||
| # Predictions are log-hazard ratios (higher = more risk) | ||
| preds = gbm.predict(X_val, num_iteration=gbm.best_iteration) | ||
| print(f"\nPrediction range: [{preds.min():.3f}, {preds.max():.3f}]") | ||
| Original file line number | Diff line number | Diff line change | ||||||
|---|---|---|---|---|---|---|---|---|
| @@ -0,0 +1,188 @@ | ||||||||
| /*! | ||||||||
| * Copyright (c) 2016-2026 Microsoft Corporation. All rights reserved. | ||||||||
| * Copyright (c) 2016-2026 The LightGBM developers. All rights reserved. | ||||||||
|
Comment on lines
+2
to
+3
Member
There was a problem hiding this comment.
Suggested change
(and XGBoost, if this has anything copied verbatim from there) |
||||||||
| * Licensed under the MIT License. See LICENSE file in the project root for license information. | ||||||||
| */ | ||||||||
| #ifndef LIGHTGBM_SRC_METRIC_COX_SURVIVAL_METRIC_HPP_ | ||||||||
| #define LIGHTGBM_SRC_METRIC_COX_SURVIVAL_METRIC_HPP_ | ||||||||
|
|
||||||||
| #include <LightGBM/metric.h> | ||||||||
| #include <LightGBM/utils/log.h> | ||||||||
|
|
||||||||
| #include <string> | ||||||||
| #include <algorithm> | ||||||||
| #include <cmath> | ||||||||
| #include <vector> | ||||||||
|
|
||||||||
| namespace LightGBM { | ||||||||
|
|
||||||||
| /*! | ||||||||
| * \brief Negative partial log-likelihood metric for Cox PH models (Breslow's method). | ||||||||
| * | ||||||||
| * Labels encode censoring via sign: +t = event at time t, -t = censored at t. | ||||||||
| * Lower is better. | ||||||||
| */ | ||||||||
| class CoxNLLMetric : public Metric { | ||||||||
| public: | ||||||||
| explicit CoxNLLMetric(const Config&) {} | ||||||||
|
|
||||||||
| ~CoxNLLMetric() {} | ||||||||
|
|
||||||||
| void Init(const Metadata& metadata, data_size_t num_data) override { | ||||||||
| name_.emplace_back("survival_cox_nll"); | ||||||||
| num_data_ = num_data; | ||||||||
| label_ = metadata.label(); | ||||||||
|
|
||||||||
| // Build sorted indices by ascending |label| (survival time) | ||||||||
| sorted_indices_.resize(num_data_); | ||||||||
| for (data_size_t i = 0; i < num_data_; ++i) { | ||||||||
| sorted_indices_[i] = i; | ||||||||
| } | ||||||||
| std::stable_sort(sorted_indices_.begin(), sorted_indices_.end(), | ||||||||
| [this](data_size_t a, data_size_t b) { | ||||||||
| return std::fabs(label_[a]) < std::fabs(label_[b]); | ||||||||
| }); | ||||||||
| } | ||||||||
|
|
||||||||
| const std::vector<std::string>& GetName() const override { | ||||||||
| return name_; | ||||||||
| } | ||||||||
|
|
||||||||
| double factor_to_bigger_better() const override { | ||||||||
| return -1.0; // lower is better | ||||||||
| } | ||||||||
|
|
||||||||
| std::vector<double> Eval(const double* score, const ObjectiveFunction*) const override { | ||||||||
| // Breslow forward-pass to compute negative partial log-likelihood | ||||||||
| double max_p = score[sorted_indices_[0]]; | ||||||||
| for (data_size_t k = 1; k < num_data_; ++k) { | ||||||||
| const data_size_t idx = sorted_indices_[k]; | ||||||||
| if (score[idx] > max_p) { | ||||||||
| max_p = score[idx]; | ||||||||
| } | ||||||||
| } | ||||||||
|
|
||||||||
| double exp_p_sum = 0.0; | ||||||||
| for (data_size_t k = 0; k < num_data_; ++k) { | ||||||||
| exp_p_sum += std::exp(score[sorted_indices_[k]] - max_p); | ||||||||
| } | ||||||||
|
|
||||||||
| double last_exp_p = 0.0; | ||||||||
| double last_abs_y = 0.0; | ||||||||
| double accumulated_sum = 0.0; | ||||||||
| double pll = 0.0; | ||||||||
| int n_events = 0; | ||||||||
|
|
||||||||
| for (data_size_t k = 0; k < num_data_; ++k) { | ||||||||
| const data_size_t idx = sorted_indices_[k]; | ||||||||
| const double p = score[idx]; | ||||||||
| const double exp_p = std::exp(p - max_p); | ||||||||
| const double y = static_cast<double>(label_[idx]); | ||||||||
| const double abs_y = std::fabs(y); | ||||||||
|
|
||||||||
| accumulated_sum += last_exp_p; | ||||||||
| if (last_abs_y < abs_y) { | ||||||||
| exp_p_sum -= accumulated_sum; | ||||||||
| accumulated_sum = 0.0; | ||||||||
| } | ||||||||
|
|
||||||||
| if (y > 0) { | ||||||||
| // p - log(sum exp(p_k)) = (p - max_p) - log(sum exp(p_k - max_p)) | ||||||||
| pll += (p - max_p) - std::log(exp_p_sum); | ||||||||
| n_events += 1; | ||||||||
| } | ||||||||
|
|
||||||||
| last_abs_y = abs_y; | ||||||||
| last_exp_p = exp_p; | ||||||||
| } | ||||||||
|
|
||||||||
| double loss = -pll / std::max(n_events, 1); | ||||||||
| return std::vector<double>(1, loss); | ||||||||
| } | ||||||||
|
|
||||||||
| private: | ||||||||
| data_size_t num_data_; | ||||||||
| const label_t* label_; | ||||||||
| std::vector<data_size_t> sorted_indices_; | ||||||||
| std::vector<std::string> name_; | ||||||||
| }; | ||||||||
|
|
||||||||
| /*! | ||||||||
| * \brief Harrell's concordance index metric for Cox PH models. | ||||||||
| * | ||||||||
| * Higher predictions = higher risk. A pair (i, j) is comparable if subject i | ||||||||
| * had an event and T_i < |T_j|. | ||||||||
| * Returns value in [0, 1] where 0.5 = random, 1.0 = perfect. | ||||||||
| * Higher is better. | ||||||||
| */ | ||||||||
| class ConcordanceIndexMetric : public Metric { | ||||||||
| public: | ||||||||
| explicit ConcordanceIndexMetric(const Config&) {} | ||||||||
|
|
||||||||
| ~ConcordanceIndexMetric() {} | ||||||||
|
|
||||||||
| void Init(const Metadata& metadata, data_size_t num_data) override { | ||||||||
| name_.emplace_back("concordance_index"); | ||||||||
| num_data_ = num_data; | ||||||||
| label_ = metadata.label(); | ||||||||
| } | ||||||||
|
|
||||||||
| const std::vector<std::string>& GetName() const override { | ||||||||
| return name_; | ||||||||
| } | ||||||||
|
|
||||||||
| double factor_to_bigger_better() const override { | ||||||||
| return 1.0; // higher is better | ||||||||
| } | ||||||||
|
|
||||||||
| std::vector<double> Eval(const double* score, const ObjectiveFunction*) const override { | ||||||||
| int64_t concordant = 0; | ||||||||
| int64_t discordant = 0; | ||||||||
| int64_t tied_risk = 0; | ||||||||
|
|
||||||||
| for (data_size_t i = 0; i < num_data_; ++i) { | ||||||||
| const double y_i = static_cast<double>(label_[i]); | ||||||||
| if (y_i <= 0) { | ||||||||
| continue; // i must have an event | ||||||||
| } | ||||||||
| const double t_i = y_i; | ||||||||
| for (data_size_t j = 0; j < num_data_; ++j) { | ||||||||
| if (i == j) { | ||||||||
| continue; | ||||||||
| } | ||||||||
| const double y_j = static_cast<double>(label_[j]); | ||||||||
| const double t_j = std::fabs(y_j); | ||||||||
| if (t_j < t_i) { | ||||||||
| continue; // j's time is strictly earlier | ||||||||
| } | ||||||||
| if (t_j == t_i && y_j > 0) { | ||||||||
| continue; // j also had event at the same time — not comparable | ||||||||
| } | ||||||||
| // comparable: j survived past t_i, or was censored at t_i | ||||||||
| if (score[i] > score[j]) { | ||||||||
| concordant += 1; | ||||||||
| } else if (score[i] < score[j]) { | ||||||||
| discordant += 1; | ||||||||
| } else { | ||||||||
| tied_risk += 1; | ||||||||
| } | ||||||||
| } | ||||||||
| } | ||||||||
|
|
||||||||
| const int64_t total = concordant + discordant + tied_risk; | ||||||||
| double c_index = 0.5; | ||||||||
| if (total > 0) { | ||||||||
| c_index = (static_cast<double>(concordant) + 0.5 * static_cast<double>(tied_risk)) | ||||||||
| / static_cast<double>(total); | ||||||||
| } | ||||||||
| return std::vector<double>(1, c_index); | ||||||||
| } | ||||||||
|
|
||||||||
| private: | ||||||||
| data_size_t num_data_; | ||||||||
| const label_t* label_; | ||||||||
| std::vector<std::string> name_; | ||||||||
| }; | ||||||||
|
|
||||||||
| } // namespace LightGBM | ||||||||
| #endif // LIGHTGBM_SRC_METRIC_COX_SURVIVAL_METRIC_HPP_ | ||||||||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -8,6 +8,7 @@ | |
| #include <string> | ||
|
|
||
| #include "binary_metric.hpp" | ||
| #include "cox_survival_metric.hpp" | ||
| #include "map_metric.hpp" | ||
| #include "multiclass_metric.hpp" | ||
| #include "rank_metric.hpp" | ||
|
|
@@ -81,6 +82,12 @@ Metric* Metric::CreateMetric(const std::string& type, const Config& config) { | |
| } else if (type == std::string("r2")) { | ||
| Log::Warning("Metric r2 is not implemented in cuda version. Fall back to evaluation on CPU."); | ||
| return new R2Metric(config); | ||
| } else if (type == std::string("survival_cox_nll")) { | ||
| Log::Warning("Metric survival_cox_nll is not implemented in cuda version. Fall back to evaluation on CPU."); | ||
| return new CoxNLLMetric(config); | ||
| } else if (type == std::string("concordance_index") || type == std::string("c_index")) { | ||
|
Comment on lines
+85
to
+88
Member
There was a problem hiding this comment. Please update this mapping in the R package as well: LightGBM/R-package/R/metrics.R Line 9 in a7d00a9 If you're comfortable writing R code we'd welcome new tests in the R package too, but at a minimum that mapping should be updated so the R package's early stopping behavior will be correct. |
||
| Log::Warning("Metric concordance_index is not implemented in cuda version. Fall back to evaluation on CPU."); | ||
| return new ConcordanceIndexMetric(config); | ||
| } | ||
| } else { | ||
| #endif // USE_CUDA | ||
|
|
@@ -132,6 +139,10 @@ Metric* Metric::CreateMetric(const std::string& type, const Config& config) { | |
| return new TweedieMetric(config); | ||
| } else if (type == std::string("r2")) { | ||
| return new R2Metric(config); | ||
| } else if (type == std::string("survival_cox_nll")) { | ||
| return new CoxNLLMetric(config); | ||
| } else if (type == std::string("concordance_index") || type == std::string("c_index")) { | ||
| return new ConcordanceIndexMetric(config); | ||
| } | ||
| #ifdef USE_CUDA | ||
| } | ||
|
|
||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Are these aliases used in other projects or research?
If not, let's please not use any aliases for this objective. Aliases add complexity and maintenance burden, and I'd especially like to avoid committing
survivallike this in case other survival objectives are added in the future.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I agree about removing the name
survival.survival:coxandcox-nloglik.CoxPHSurvivalAnalysiscoxphwithties=“breslow”CoxPHFitterThere was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Awesome, thanks for those links! That's exactly the type of thing I was looking for. Based on that, I'm happy with dropping
survivalbut keepingcoxandcox_ph.