Repository: sebp/scikit-survival (language: python)

Issue:
RandomSurvivalForest is not consistent making predictions of survival functions at certain times
<!--
Before submitting a bug, please make sure the issue hasn't been already
addressed by searching through the past issues.
-->

**Describe the bug**

RandomSurvivalForest is not consistent making predictions of survival functions (predict_survival_function) at certain times. I would expect that RandomSurvivalForest can predict survival functions when the times are within the intervals of the training and validation times. However, while the intervals of both the training and validation times are [0,182], RandomSurvivalForest cannot predict the survival function on times between [0,175]. Note that this doesn't happen with CoxPHSurvivalAnalysis
 or GradientBoostingSurvivalAnalysis

**Code Sample to Reproduce the Bug**

```python
import numpy as np
import pandas as pd
import statsmodels.api as sm
pharmacoSmoking = sm.datasets.get_rdataset("pharmacoSmoking", "asaur")
data = pharmacoSmoking.data
data = data.drop(columns=["id","ageGroup2","ageGroup4"]) # Drop redundant information and ids
from sksurv.datasets import get_x_y
X, y = get_x_y(data, attr_labels=["relapse", "ttr"], pos_label=True)
for c in X.columns:
    if X[c].dtype.kind not in ['i', 'f']:
        X[c] = X[c].astype("category")
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sksurv.preprocessing import OneHotEncoder
from sksurv.ensemble import RandomSurvivalForest
from sksurv.metrics import integrated_brier_score
def get_ibs(seed, times):
  X_trn, X_val, y_trn, y_val = train_test_split(X, y, random_state=seed)
  print(f"Minimum training time: {y_trn['ttr'].min()}")
  print(f"Maximum training time: {y_trn['ttr'].max()}")
  print(f"Minimum validation time: {y_val['ttr'].min()}")
  print(f"Maximum validation time: {y_val['ttr'].max()}")
  enc = OneHotEncoder()
  scaler = StandardScaler()
  X_trn = enc.fit_transform(X_trn)
  X_trn = pd.DataFrame(scaler.fit_transform(X_trn), columns=X_trn.columns)
  X_val = enc.transform(X_val)
  X_val = pd.DataFrame(scaler.transform(X_val), columns=X_val.columns)
  rsf = RandomSurvivalForest(random_state=42)
  rsf.fit(X_trn, y_trn)
  survs = rsf.predict_survival_function(X_val)
  preds = np.asarray([[fn(t) for t in times] for fn in survs])
  return integrated_brier_score(y_trn, y_val, preds, times)
print(f"Integrated Brier Score: {get_ibs(0, np.arange(0,170))}")
print(f"Integrated Brier Score: {get_ibs(0, np.arange(0,175))}")
```

**Expected Results**

Minimum training time: 0.0
Maximum training time: 182.0
Minimum validation time: 0.0
Maximum validation time: 182.0
Integrated Brier Score: 0.21497683508486282
Minimum training time: 0.0
Maximum training time: 182.0
Minimum validation time: 0.0
Maximum validation time: 182.0
Integrated Brier Score: 0.21497683508486282

**Actual Results**

Minimum training time: 0.0
Maximum training time: 182.0
Minimum validation time: 0.0
Maximum validation time: 182.0
Integrated Brier Score: 0.21497683508486282
Minimum training time: 0.0
Maximum training time: 182.0
Minimum validation time: 0.0
Maximum validation time: 182.0

```
ValueError                                Traceback (most recent call last)
[<ipython-input-3-fbfef721f54b>](https://localhost:8080/#) in <module>()
     33   return integrated_brier_score(y_trn, y_val, preds, times)
     34 print(f"Integrated Brier Score: {get_ibs(0, np.arange(0,170))}")
---> 35 print(f"Integrated Brier Score: {get_ibs(0, np.arange(0,175))}")

3 frames
[/usr/local/lib/python3.7/dist-packages/sksurv/functions.py](https://localhost:8080/#) in __call__(self, x)
     65         if numpy.min(x) < self.x[0] or numpy.max(x) > self.x[-1]:
     66             raise ValueError(
---> 67                 "x must be within [%f; %f]" % (self.x[0], self.x[-1]))
     68         i = numpy.searchsorted(self.x, x, side='left')
     69         not_exact = self.x[i] != x

ValueError: x must be within [0.000000; 170.000000]
```

**Versions**
Please execute the following snippet and paste the output below.

```python
import sklearn; sklearn.show_versions()
import sksurv; print("sksurv:", sksurv.__version__)
import cvxopt; print("cvxopt:", cvxopt.__version__)
import cvxpy; print("cvxpy:", cvxpy.__version__)
import numexpr; print("numexpr:", numexpr.__version__)
import osqp; print("osqp:", osqp.OSQP().version())
```
System:
    python: 3.7.12 (default, Jan 15 2022, 18:48:18)  [GCC 7.5.0]
executable: /usr/bin/python3
   machine: Linux-5.4.144+-x86_64-with-Ubuntu-18.04-bionic

Python dependencies:
          pip: 21.1.3
   setuptools: 57.4.0
      sklearn: 1.0.2
        numpy: 1.21.5
        scipy: 1.4.1
       Cython: 0.29.27
       pandas: 1.3.5
   matplotlib: 3.2.2
       joblib: 1.1.0
threadpoolctl: 3.1.0

Built with OpenMP: True
sksurv: 0.17.0
cvxopt: 1.2.7
cvxpy: 1.0.31
numexpr: 2.8.1
osqp: 0.6.2
RandomSurvivalForest is not consistent making predictions of survival functions at certain times
<!--
Before submitting a bug, please make sure the issue hasn't been already
addressed by searching through the past issues.
-->

**Describe the bug**

RandomSurvivalForest is not consistent making predictions of survival functions (predict_survival_function) at certain times. I would expect that RandomSurvivalForest can predict survival functions when the times are within the intervals of the training and validation times. However, while the intervals of both the training and validation times are [0,182], RandomSurvivalForest cannot predict the survival function on times between [0,175]. Note that this doesn't happen with CoxPHSurvivalAnalysis
 or GradientBoostingSurvivalAnalysis

**Code Sample to Reproduce the Bug**

```python
import numpy as np
import pandas as pd
import statsmodels.api as sm
pharmacoSmoking = sm.datasets.get_rdataset("pharmacoSmoking", "asaur")
data = pharmacoSmoking.data
data = data.drop(columns=["id","ageGroup2","ageGroup4"]) # Drop redundant information and ids
from sksurv.datasets import get_x_y
X, y = get_x_y(data, attr_labels=["relapse", "ttr"], pos_label=True)
for c in X.columns:
    if X[c].dtype.kind not in ['i', 'f']:
        X[c] = X[c].astype("category")
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sksurv.preprocessing import OneHotEncoder
from sksurv.ensemble import RandomSurvivalForest
from sksurv.metrics import integrated_brier_score
def get_ibs(seed, times):
  X_trn, X_val, y_trn, y_val = train_test_split(X, y, random_state=seed)
  print(f"Minimum training time: {y_trn['ttr'].min()}")
  print(f"Maximum training time: {y_trn['ttr'].max()}")
  print(f"Minimum validation time: {y_val['ttr'].min()}")
  print(f"Maximum validation time: {y_val['ttr'].max()}")
  enc = OneHotEncoder()
  scaler = StandardScaler()
  X_trn = enc.fit_transform(X_trn)
  X_trn = pd.DataFrame(scaler.fit_transform(X_trn), columns=X_trn.columns)
  X_val = enc.transform(X_val)
  X_val = pd.DataFrame(scaler.transform(X_val), columns=X_val.columns)
  rsf = RandomSurvivalForest(random_state=42)
  rsf.fit(X_trn, y_trn)
  survs = rsf.predict_survival_function(X_val)
  preds = np.asarray([[fn(t) for t in times] for fn in survs])
  return integrated_brier_score(y_trn, y_val, preds, times)
print(f"Integrated Brier Score: {get_ibs(0, np.arange(0,170))}")
print(f"Integrated Brier Score: {get_ibs(0, np.arange(0,175))}")
```

**Expected Results**

Minimum training time: 0.0
Maximum training time: 182.0
Minimum validation time: 0.0
Maximum validation time: 182.0
Integrated Brier Score: 0.21497683508486282
Minimum training time: 0.0
Maximum training time: 182.0
Minimum validation time: 0.0
Maximum validation time: 182.0
Integrated Brier Score: 0.21497683508486282

**Actual Results**

Minimum training time: 0.0
Maximum training time: 182.0
Minimum validation time: 0.0
Maximum validation time: 182.0
Integrated Brier Score: 0.21497683508486282
Minimum training time: 0.0
Maximum training time: 182.0
Minimum validation time: 0.0
Maximum validation time: 182.0

```
ValueError                                Traceback (most recent call last)
[<ipython-input-3-fbfef721f54b>](https://localhost:8080/#) in <module>()
     33   return integrated_brier_score(y_trn, y_val, preds, times)
     34 print(f"Integrated Brier Score: {get_ibs(0, np.arange(0,170))}")
---> 35 print(f"Integrated Brier Score: {get_ibs(0, np.arange(0,175))}")

3 frames
[/usr/local/lib/python3.7/dist-packages/sksurv/functions.py](https://localhost:8080/#) in __call__(self, x)
     65         if numpy.min(x) < self.x[0] or numpy.max(x) > self.x[-1]:
     66             raise ValueError(
---> 67                 "x must be within [%f; %f]" % (self.x[0], self.x[-1]))
     68         i = numpy.searchsorted(self.x, x, side='left')
     69         not_exact = self.x[i] != x

ValueError: x must be within [0.000000; 170.000000]
```

**Versions**
Please execute the following snippet and paste the output below.

```python
import sklearn; sklearn.show_versions()
import sksurv; print("sksurv:", sksurv.__version__)
import cvxopt; print("cvxopt:", cvxopt.__version__)
import cvxpy; print("cvxpy:", cvxpy.__version__)
import numexpr; print("numexpr:", numexpr.__version__)
import osqp; print("osqp:", osqp.OSQP().version())
```
System:
    python: 3.7.12 (default, Jan 15 2022, 18:48:18)  [GCC 7.5.0]
executable: /usr/bin/python3
   machine: Linux-5.4.144+-x86_64-with-Ubuntu-18.04-bionic

Python dependencies:
          pip: 21.1.3
   setuptools: 57.4.0
      sklearn: 1.0.2
        numpy: 1.21.5
        scipy: 1.4.1
       Cython: 0.29.27
       pandas: 1.3.5
   matplotlib: 3.2.2
       joblib: 1.1.0
threadpoolctl: 3.1.0

Built with OpenMP: True
sksurv: 0.17.0
cvxopt: 1.2.7
cvxpy: 1.0.31
numexpr: 2.8.1
osqp: 0.6.2


Discussion / hints:
There's indeed an inconsistency. For `RandomSurvivalForest` and `SurvivalTree`, only the domain of the survival function is limited to the range of *event* times, whereas for `CoxPHSurvivalAnalysis` it is limited to *any* time point, whether an event has been observed or not.
#363 adds a test describing (my understanding of) the current behaviour described in this issue.

Wondering what a fix might look like?
* All models to support all-times range?
* All models to only support event-times range?
* Some models supporting all-times range and other models only supporting event-times range with some way to tell which is supported before making the `predict_survival_function` or `fn(t)` call?
* Something else?

@cpoerschke I guess the question goes to @sebp however for me all models should support all-times range. If they give wrong predictions outside the trained ranges so be it.
If with "all-times ranges" you mean the range of the training data (disregarding censoring status), I would agree.

It would also be worthwhile to check how random survival forest or gradient boosting implementation in R handle this.
Thank you very much @sebp I don't know how those implementations do it in R, however, scikit-learn in random forest model doesn't put constraints on users if they want to go outside the training data (even though the results are obviously not good outside the training region). Here is a gist of that behavior that you can run in colab as well: https://gist.github.com/alonsosilvaallende/ef813de35a8b0f4328b451aea46b9c48
https://github.com/sebp/scikit-survival/pull/363/commits/414354ac10c69609b05fc8f725aa3255b86bcf56 explores clipping inputs on an opt-in basis.
@cpoerschke Thanks for adding the clipping option. This might be a good option for advanced users that are aware of the risks. For time points that are smaller than the smallest time point in the training data, it would actually be okay to predict a survival probability of 1. For time points beyond the far end of the time axis, predictions would be very speculative.

Regarding that `RandomSurvivalForest` currently limits the interval to *event* times, I'm probably going to fix that together with #343.
There's indeed an inconsistency. For `RandomSurvivalForest` and `SurvivalTree`, only the domain of the survival function is limited to the range of *event* times, whereas for `CoxPHSurvivalAnalysis` it is limited to *any* time point, whether an event has been observed or not.
#363 adds a test describing (my understanding of) the current behaviour described in this issue.

Wondering what a fix might look like?
* All models to support all-times range?
* All models to only support event-times range?
* Some models supporting all-times range and other models only supporting event-times range with some way to tell which is supported before making the `predict_survival_function` or `fn(t)` call?
* Something else?

@cpoerschke I guess the question goes to @sebp however for me all models should support all-times range. If they give wrong predictions outside the trained ranges so be it.
If with "all-times ranges" you mean the range of the training data (disregarding censoring status), I would agree.

It would also be worthwhile to check how random survival forest or gradient boosting implementation in R handle this.
Thank you very much @sebp I don't know how those implementations do it in R, however, scikit-learn in random forest model doesn't put constraints on users if they want to go outside the training data (even though the results are obviously not good outside the training region). Here is a gist of that behavior that you can run in colab as well: https://gist.github.com/alonsosilvaallende/ef813de35a8b0f4328b451aea46b9c48
https://github.com/sebp/scikit-survival/pull/363/commits/414354ac10c69609b05fc8f725aa3255b86bcf56 explores clipping inputs on an opt-in basis.
@cpoerschke Thanks for adding the clipping option. This might be a good option for advanced users that are aware of the risks. For time points that are smaller than the smallest time point in the training data, it would actually be okay to predict a survival probability of 1. For time points beyond the far end of the time axis, predictions would be very speculative.

Regarding that `RandomSurvivalForest` currently limits the interval to *event* times, I'm probably going to fix that together with #343.

Fix patch:
diff --git a/sksurv/functions.py b/sksurv/functions.py
index 5df79df1..ee978cc8 100644
--- a/sksurv/functions.py
+++ b/sksurv/functions.py
@@ -38,18 +38,46 @@ class StepFunction:
 
     b : float, optional, default: 0.0
         Constant offset term.
+
+    domain : tuple, optional
+        A tuple with two entries that sets the limits of the
+        domain of the step function.
+        If entry is `None`, use the first/last value of `x` as limit.
     """
 
-    def __init__(self, x, y, *, a=1.0, b=0.0):
+    def __init__(self, x, y, *, a=1.0, b=0.0, domain=(0, None)):
         check_consistent_length(x, y)
         self.x = x
         self.y = y
         self.a = a
         self.b = b
+        domain_lower = self.x[0] if domain[0] is None else domain[0]
+        domain_upper = self.x[-1] if domain[1] is None else domain[1]
+        self._domain = (float(domain_lower), float(domain_upper))
+
+    @property
+    def domain(self):
+        """Returns the domain of the function, that means
+        the range of values that the function accepts.
+
+        Returns
+        -------
+        lower_limit : float
+            Lower limit of domain.
+
+        upper_limit : float
+            Upper limit of domain.
+        """
+        return self._domain
 
     def __call__(self, x):
         """Evaluate step function.
 
+        Values outside the interval specified by `self.domain`
+        will raise an exception.
+        Values in `x` that are in the interval `[self.domain[0]; self.x[0]]`
+        get mapped to `self.y[0]`.
+
         Parameters
         ----------
         x : float|array-like, shape=(n_values,)
@@ -63,8 +91,12 @@ def __call__(self, x):
         x = np.atleast_1d(x)
         if not np.isfinite(x).all():
             raise ValueError("x must be finite")
-        if np.min(x) < self.x[0] or np.max(x) > self.x[-1]:
-            raise ValueError(f"x must be within [{self.x[0]:f}; {self.x[-1]:f}]")
+        if np.min(x) < self._domain[0] or np.max(x) > self.domain[1]:
+            raise ValueError(f"x must be within [{self.domain[0]:f}; {self.domain[1]:f}]")
+
+        # x is within the domain, but we need to account for self.domain[0] <= x < self.x[0]
+        x = np.clip(x, a_min=self.x[0], a_max=None)
+
         i = np.searchsorted(self.x, x, side="left")
         not_exact = self.x[i] != x
         i[not_exact] -= 1


Test patch:
diff --git a/tests/test_functions.py b/tests/test_functions.py
index 1467d094..7414ae14 100644
--- a/tests/test_functions.py
+++ b/tests/test_functions.py
@@ -3,6 +3,7 @@
 import pytest
 
 from sksurv.functions import StepFunction
+from sksurv.testing import all_survival_estimators
 
 
 @pytest.fixture()
@@ -13,6 +14,25 @@ def a_step_function():
     return f
 
 
+@pytest.fixture()
+def toy_data_exponential():
+    rnd = np.random.RandomState(2)
+    n_samples = 100
+    x = rnd.randn(n_samples, 2)
+    y = np.empty(n_samples, dtype=[("event", bool), ("time", float)])
+    y["time"] = rnd.exponential(scale=np.exp(x[:, 0]), size=n_samples)
+    y["event"] = rnd.binomial(1, 0.5, size=n_samples) == 1
+
+    # ensure at least 2 uncensored events exist
+    y["event"][:2] = True
+
+    # mark entry with largest time as censored
+    # see https://github.com/sebp/scikit-survival/issues/249
+    idxmax = np.argmax(y["time"])
+    y["event"][idxmax] = False
+    return x, y
+
+
 class TestStepFunction:
     @staticmethod
     def test_exact(a_step_function):
@@ -26,18 +46,14 @@ def test_not_exact(a_step_function):
         assert_array_equal(actual, a_step_function.y[:-1])
 
     @staticmethod
-    def test_out_of_bounds(a_step_function):
-        eps = np.finfo(float).eps * 8
-        values = [
-            a_step_function.x[0] - 100,
-            a_step_function.x[-1] + 100,
-            a_step_function.x[0] - eps,
-            a_step_function.x[-1] + eps,
-        ]
-
-        for v in values:
-            with pytest.raises(ValueError, match=r"x must be within \[0.0+; 9.0+\]"):
-                a_step_function(v)
+    @pytest.mark.parametrize("value", [-100, 100, -np.finfo(float).eps * 8, np.finfo(float).eps * 8])
+    def test_out_of_bounds(a_step_function, value):
+        v = value
+        if v > 0:
+            v += a_step_function.domain[1]
+
+        with pytest.raises(ValueError, match=r"x must be within \[0\.0+; 9\.0+\]"):
+            a_step_function(v)
 
     @staticmethod
     def test_not_finite(a_step_function, non_finite_value):
@@ -56,3 +72,87 @@ def test_equal(a_step_function):
         assert a_step_function != different_step_function
 
         assert a_step_function != x
+
+
+@pytest.mark.parametrize(
+    "estimator_cls", [est for est in all_survival_estimators() if hasattr(est, "predict_cumulative_hazard_function")]
+)
+def test_predict_cumulative_hazard_function_range(estimator_cls, toy_data_exponential):
+    x, y = toy_data_exponential
+
+    estimator = estimator_cls()
+    if "fit_baseline_model" in estimator.get_params():
+        estimator.set_params(fit_baseline_model=True)
+    estimator.fit(x, y)
+
+    t_min = y["time"].min()
+    t_max = y["time"].max()
+    t_mid = (t_max - t_min) / 2.0
+
+    for fn in estimator.predict_cumulative_hazard_function(x):
+        v = fn(t_min)
+        assert np.isfinite(v)
+        assert v == 0
+
+    for fn in estimator.predict_cumulative_hazard_function(x):
+        v = fn(t_mid)
+        assert np.isfinite(v)
+        assert v >= 0
+
+    t_smaller_min = t_min / 2
+    for fn in estimator.predict_cumulative_hazard_function(x):
+        v = fn(t_smaller_min)
+        assert np.isfinite(v)
+        assert v == 0
+
+    for fn in estimator.predict_cumulative_hazard_function(x):
+        v = fn(t_max)
+        assert np.isfinite(v)
+        assert v >= 0
+
+    t_bigger_max = t_max + 1
+    for fn in estimator.predict_cumulative_hazard_function(x):
+        with pytest.raises(ValueError, match=r"x must be within \[[0-9.]+; [0-9.]+\]"):
+            fn(t_bigger_max)
+
+
+@pytest.mark.parametrize(
+    "estimator_cls", [est for est in all_survival_estimators() if hasattr(est, "predict_survival_function")]
+)
+def test_predict_survival_function_range(estimator_cls, toy_data_exponential):
+    x, y = toy_data_exponential
+
+    estimator = estimator_cls()
+    if "fit_baseline_model" in estimator.get_params():
+        estimator.set_params(fit_baseline_model=True)
+    estimator.fit(x, y)
+
+    t_min = y["time"].min()
+    t_max = y["time"].max()
+    t_mid = (t_max - t_min) / 2.0
+
+    for fn in estimator.predict_survival_function(x):
+        v = fn(t_min)
+        assert np.isfinite(v)
+        assert v == 1
+
+    for fn in estimator.predict_survival_function(x):
+        v = fn(t_mid)
+        assert np.isfinite(v)
+        assert 0 <= v <= 1
+
+    t_smaller_min = t_min / 2
+    for fn in estimator.predict_survival_function(x):
+        v = fn(t_smaller_min)
+        assert np.isfinite(v)
+        assert v == 1
+
+    for fn in estimator.predict_survival_function(x):
+        v = fn(t_max)
+        assert np.isfinite(v)
+        assert 0 <= v <= 1
+
+    t_bigger_max = t_max + 1
+    for fn in estimator.predict_survival_function(x):
+        with pytest.raises(ValueError, match=r"x must be within \[[0-9.]+; [0-9.]+\]"):
+            fn(t_bigger_max)