j'ai aps réussi dans un seul commit

hazardous/_fine_and_gray.py

@@ -8,6 +8,10 @@
 from sklearn.base import BaseEstimator, check_is_fitted
 from sklearn.utils import check_random_state
 
+from hazardous.metrics._brier_score import (
+    integrated_brier_score_incidence,
+    integrated_brier_score_incidence_oracle,
+)
 from hazardous.utils import check_y_survival
 
 r_cmprsk = importr("cmprsk")
@@ -106,11 +110,9 @@ class FineGrayEstimator(BaseEstimator):
 
     def __init__(
         self,
-        event_of_interest=1,
         max_fit_samples=10_000,
         random_state=None,
     ):
-        self.event_of_interest = event_of_interest
         self.max_fit_samples = max_fit_samples
         self.random_state = random_state
 
@@ -144,20 +146,31 @@ def fit(self, X, y):
 
         y = y.loc[X.index]
         event, duration = check_y_survival(y)
-        self.times_ = np.unique(duration[event == self.event_of_interest])
+        self.times_ = np.unique(duration[event > 0])
+        self.event_ids_ = np.array(sorted(list(set([0]) | set(event))))
+
         event, duration = r_vector(event), r_vector(duration)
+
         X = r_dataframe(X)
 
-        self.r_crr_result_ = r_cmprsk.crr(
-            duration,
-            event,
-            X,
-            failcode=self.event_of_interest,
-            cencode=0,
-        )
-        parsed = parse_r_list(self.r_crr_result_)
-        self.coef_ = np.array(parsed["coef"])
+        self.r_crr_results_ = []
+        self.coefs_ = []
+        for event_id in self.event_ids_[1:]:
+            r_crr_result_ = r_cmprsk.crr(
+                duration,
+                event,
+                X,
+                failcode=int(event_id),
+                cencode=0,
+            )
+
+            parsed = parse_r_list(r_crr_result_)
+            coef_ = np.array(parsed["coef"])
 
+            self.r_crr_results_.append(r_crr_result_)
+            self.coefs_.append(coef_)
+
+        self.y_train = y
         return self
 
     def predict_cumulative_incidence(self, X, times=None):
@@ -178,35 +191,51 @@ def predict_cumulative_incidence(self, X, times=None):
         y_pred : ndarray of shape (n_samples, n_times)
             The conditional cumulative cumulative incidence at times.
         """
-        check_is_fitted(self, "r_crr_result_")
+        check_is_fitted(self, "r_crr_results_")
 
         self._check_feature_names(X, reset=False)
         self._check_n_features(X, reset=False)
 
         X = r_matrix(X)
 
-        y_pred = r_cmprsk.predict_crr(
-            self.r_crr_result_,
-            X,
-        )
-        y_pred = np_matrix(y_pred)
+        all_event_y_pred = []
+        for event_id in self.event_ids_[1:]:
+            # Interpolate each sample
+            y_pred = r_cmprsk.predict_crr(
+                self.r_crr_results_[event_id - 1],
+                X,
+            )
+            y_pred = np_matrix(y_pred)
+            times_event = np.unique(
+                self.y_train["duration"][self.y_train["event"] == event_id]
+            )
+            y_pred = y_pred[:, 1:].T  # shape (n_samples, n_unique_times)
 
-        unique_times = y_pred[:, 0]  # durations seen during fit
-        y_pred = y_pred[:, 1:].T  # shape (n_samples, n_unique_times)
+            y_pred_at_0 = np.zeros((y_pred.shape[0], 1))
+            y_pred_t_max = y_pred[:, [-1]]
+            y_pred = np.hstack([y_pred_at_0, y_pred, y_pred_t_max])
+
+            times_event = np.hstack([[0], times_event, [np.inf]])
+
+            if times is None:
+                times = self.times_
 
-        # Interpolate each sample
-        if times is not None:
             all_y_pred = []
-            for idx in y_pred.shape[0]:
+            for idx in range(y_pred.shape[0]):
                 y_pred_ = interp1d(
-                    x=unique_times,
+                    x=times_event,
                     y=y_pred[idx, :],
                     kind="linear",
                 )(times)
                 all_y_pred.append(y_pred_)
+
             y_pred = np.vstack(all_y_pred)
+            all_event_y_pred.append(y_pred)
 
-        return y_pred
+        surv_curve = 1 - np.sum(all_event_y_pred, axis=0)
+        all_event_y_pred = [surv_curve] + all_event_y_pred
+        all_event_y_pred = np.asarray(all_event_y_pred)
+        return all_event_y_pred
 
     def _check_input(self, X, y):
         if not hasattr(X, "__dataframe__"):
@@ -230,3 +259,36 @@ def _check_input(self, X, y):
             raise ValueError(f"Constant columns {constant_columns} need jittering.")
 
         return X
+
+    def score(self, X, y, shape_censoring=None, scale_censoring=None):
+        """Return
+
+        #TODO: implement time integrated NLL.
+        """
+        predicted_curves = self.predict_cumulative_incidence(X)
+        ibs_events = []
+
+        for idx, event_id in enumerate(self.event_ids_[1:]):
+            predicted_curves_for_event = predicted_curves[idx]
+            if scale_censoring is not None and shape_censoring is not None:
+                ibs_event = integrated_brier_score_incidence_oracle(
+                    y_train=self.y_train,
+                    y_test=y,
+                    y_pred=predicted_curves_for_event,
+                    times=self.times_,
+                    shape_censoring=shape_censoring,
+                    scale_censoring=scale_censoring,
+                    event_of_interest=event_id,
+                )
+
+            else:
+                ibs_event = integrated_brier_score_incidence(
+                    y_train=self.y_train,
+                    y_test=y,
+                    y_pred=predicted_curves_for_event,
+                    times=self.times_,
+                    event_of_interest=event_id,
+                )
+
+            ibs_events.append(ibs_event)
+        return -np.mean(ibs_events)