Add the c-index with IPCW

[Vincent-Maladiere]

What does this PR propose?

This PR proposes to add the c-index as defined in [1]. I think this is ready to be reviewed for merging, with some questions/suggestions in the TODO section below.

show maths

where:

and

and

where $M$ is the probability of incidence of the event of interest.

• The concordance_index_incidence function is inspired by the concordance_index function in lifelines, with some significant differences:
  • It computes the above formula, designed for competing risks and IPCW
  • It accepts truncated times taus
  • It can use virtually any estimator to compute the IPCW (only KM for now)
  • It uses increasing probabilities of incidence instead of decreasing survival probabilities
• The main advantage of the lifelines design is the use of a balanced tree, leading to a time complexity in $O ( n \times log (n))$ instead of $O(n^2)$ like in scikit-survival concordance_index_ipcw.
• To support IPCW, I extended the _BTree class from lifelines.utils.btree.py by adding a weighting count mechanism. I referenced lifelines in hazardous.metrics._btree.py, but I can reference it also in the hazardous.metrics._concordance_index.py file if necessary.
• I added an extensive test suite on a few data points, allowing me to check the results manually.
• A good effort on documentation has been started by @judithabk6

TODO

• Keep improving the documentation
• Should we add tests for bigger datasets and compare our results with scikit-survival for the survival case? And with results from a specific r-package dealing with competing risk?
• In the formula of $\hat{W}_{ij, 1}$ (see maths above), I didn't differentiate between $\tilde{T}_i-$ and $\tilde{T}_i$.
• Use the Cox IPCW estimator. This depends on first adding Cox to the IPCW estimator (to be done in another PR).
• Should we use the tied_tol parameter for ties in predictions?

cc @ogrisel @GaelVaroquaux @juAlberge @glemaitre

[1] Wolbers, M., Blanche, P., Koller, M. T., Witteman, J. C., & Gerds, T. A. (2014). Concordance for prognostic models with competing risks.

[Vincent-Maladiere]

The CI for the doc fails because the previous boosting tree model is missing. This should be fixed when #53 is merged.

[Vincent-Maladiere]

Update on performance

Our implementation is 100x slower than scikit-survival concordance_index_ipcw. This is due to the weight computing (the IPCWs) inside the BalancedTree, which lifelines doesn't perform.

code benchmark

import numpy as np
import pandas as pd
from time import time
from lifelines import CoxPHFitter
from lifelines.datasets import load_kidney_transplant
from sklearn.model_selection import train_test_split

from hazardous.metrics._concordance_index import _concordance_index_incidence_report

df = load_kidney_transplant()

# make the dataset 100x times longer for benchmarking purposes
df = pd.concat([df] * 100, axis=0)

df_train, df_test = train_test_split(df, stratify=df["death"])
cox = CoxPHFitter().fit(df_train, duration_col="time", event_col="death")

t_min, t_max = df["time"].min(), df["time"].max()
time_grid = np.linspace(t_min, t_max, 20)
y_pred = 1 - cox.predict_survival_function(df_test, times=time_grid).T.to_numpy()

y_train = df_train[["death", "time"]].rename(columns=dict(
    death="event", time="duration"
))
y_test = df_test[["death", "time"]].rename(columns=dict(
    death="event", time="duration"
))

tic = time()
result = _concordance_index_incidence_report(
    y_test=y_test,
    y_pred=y_pred,
    time_grid=time_grid,
    taus=None,
    y_train=y_train,
) 
print(f"our implementation: {time() - tic:.2f}s")

# scikit-survival
from sksurv.metrics import concordance_index_ipcw

def make_recarray(y):
    event, duration = y["event"].values, y["duration"].values
    return np.array(
        [(event[i], duration[i]) for i in range(len(event))],
        dtype=[("e", bool), ("t", float)],
    )

tic = time()
concordance_index_ipcw(
    make_recarray(y_train),
    make_recarray(y_test),
    y_pred[:, -1],
    tau=None,
)
print(f"scikit-survival: {time() - tic:.2f}s")

# lifelines
from lifelines.utils import concordance_index

concordance_index(
    event_times=y_test["duration"],
    predicted_scores=1 - y_pred[:, -1],
    event_observed=y_test["event"],
)
print(f"lifelines: {time() - tic:.2f}s")

On a dataset with 20k rows:

our implementation: 18.10s
scikit-survival: 0.24s
lifelines: 0.27s

The flamegraph is quite clear about the culprit, being the list comprehension that computes the IPCW weight for each pair. When I remove the IPCWs, the performance becomes similar to lifelines.

Speedscope views of our implementation

I tried to fix this performance issue using numba @jitclass on the BTree, but it is still very slow. I put the numba BTree on a separate draft branch for reference.

Conclusion

I only see two ways forward:

1. Either my computation of the IPCW in the BTree is flawed, and we can fix the performance issue
2. or the BTree is not adapted for our metric and we have to look at a non-optimized pairwise implementation like scikit-survival with a $O(n^2)$ instead of $n \log (n)$ time complexity. This would simplify the code base though.

[jjerphan]

Pinged by @Vincent-Maladiere, but have no time for it.

Random pile of pieces of advice:

• find if a better algorithm exist first
• profile to see what's the bottleneck
• see if tree-based structures can be used from another library (e.g. pydatastructures
• use another language (like Cython or C++) to implement the costly algorithmic part

[GaelVaroquaux]

No, don't use compiled languages, please. It will make release and distribution much harder.

…

On Jul 26, 2024, 13:46, at 13:46, Julien Jerphanion ***@***.***> wrote: Pinged by @Vincent-Maladiere, but have no time for it. Random pile of pieces of advice: - find if a better algorithm exist first - profile to see what's the bottleneck - see if tree-based structures can be used from another library (e.g. [`pydatastructures`](https://github.com/codezonediitj/pydatastructs/tree/main/pydatastructs/trees) - use another language (like Cython or C++) to implement the costly algorithmic part -- Reply to this email directly or view it on GitHub: #71 (comment) You are receiving this because you were mentioned. Message ID: ***@***.***>

[Vincent-Maladiere]

After giving it some more thought, there is room for improvement with the current balanced tree design :

1. When we don't use an IPCW estimator (like lifelines):
   $$W_{ij,1} = W_{ij,2} = 1$$
2. When we use a non-conditional IPCW estimator (Kaplan-Meier, like scikit-survival):
   $$W_{ij,1} = W_{i,1} = \hat{G}(T_i) ^ 2 \space \mathrm{and} \space W_{ij,2} = \hat{G}(T_i) \hat{G}(T_j) $$

However, when we use a conditional IPCW estimator (like Cox or SurvivalBoost), we have:
$$W_{ij,1} = \hat{G}(T_i | X_i) \hat{G}(T_i | X_j) \space \mathrm{and} \space W_{ij,2} = \hat{G}(T_i | X_i) \hat{G}(T_j | X_j)$$

In this case, the balanced tree is not adapted anymore, and we should use the naive implementation.

So, to make things simpler, I suggest we only implement the naive version for now, and eventually return to the balanced tree later, for the non-conditional and unweighted cases.

WDYT?

[GaelVaroquaux]

Sounds good to me. We can always iterate if needed

…

On Jul 26, 2024, 18:38, at 18:38, Vincent M ***@***.***> wrote: After giving it some more thought, there is room for improvement with the current balanced tree design : 1. When we don't use an IPCW estimator (like lifelines): $$W_{ij,1} = W_{ij,2} = 1$$ 2. When we use a **non-conditional** IPCW estimator (Kaplan-Meier, like scikit-survival): $$W_{ij,1} = W_{i,1} = \hat{G}(T_i) ^ 2 \space \mathrm{and} \space W_{ij,2} = \hat{G}(T_i) \hat{G}(T_j) $$ However, when we use a **conditional** IPCW estimator (like Cox or SurvivalBoost), we have: $$W_{ij,1} = \hat{G}(T_i | X_i) \hat{G}(T_i | X_j) \space \mathrm{and} \space W_{ij,2} = \hat{G}(T_i | X_i) \hat{G}(T_j | X_j)$$ In this case, the balanced tree is not adapted anymore, and we should use the naive implementation. So, to make things simpler, **I suggest we only implement the naive version for now**, and eventually return to the balanced tree later, for the non-conditional and unweighted cases. WDYT? -- Reply to this email directly or view it on GitHub: #71 (comment) You are receiving this because you were mentioned. Message ID: ***@***.***>

[Vincent-Maladiere]

Here is the revised version. When used on a survival dataset, it gives identical results to scikit-survival, with a slightly better time complexity.

[Vincent-Maladiere]

This PR is now ready to be reviewed :)

[antoinebaker]

Does type A / type B refer to the A_ij and B_ij matrices in Eq (3.4) of ref [1] ?

• type=A: j experienced any event (competing, censoring, event of interest) strictly after i $(T_j > T_i)$
• type=B: j experienced a competing event before i $(T_j \leq T_i, \Delta_j \neq 0, k)$

Are competing events after i missing in the list for pair of type A ? Should censoring events only be in type A if T_j > T_i ?

[Vincent-Maladiere]

Yes, A and B refer to Eq (3.4). You're right, this docstring is slightly inaccurate; let's rewrite it as you suggest.

[Vincent-Maladiere]

Looking deeper, that part was actually correct. Both scikit survival and reference R implementation consider ties in times with censoring (D_i = 1, D_j = 0 and T_i = T_j) as comparable pairs, although this is not mentioned in the paper [1].

Also looking a bit more at our implementation, we don't consider these.
Therefore:

1. I'll revert the documentation
2. I'll fix the implementation to consider ties in time for censored j

[antoinebaker]

Suggested change

	n_ties_times_a: list of int
	Number of tied pairs of type A with D_i = D_j = 1.
	n_ties_times_a: list of int
	Number of tied pairs of type A with T_i = T_j = 1.

[Vincent-Maladiere]

So here I believe the docstring is accurate, we consider ties of type A only individuals having experienced the event of interest simultaneously.
If you think it would bring clarity, we could write:

D_i = D_j = 1 and T_j = T_i

[antoinebaker]

This tie in time pair is not considered in ref [1] ?

[antoinebaker]

I think the type A with T_i=T_j is not consistent with ref [1] (where type A is T_i < T_j , and any event for j). Maybe you have a mix of ref [1] and sksurv.metrics.concordance_index_ipcw conventions ?

Are you maybe adding events (T_i=T_j and D_i=D_j=event_of_interest), discarded in ref [1], but here counted as "tied in time" with 1/2 weight ?

[Vincent-Maladiere]

Yes, you are right, we chose to include them to be closer to the results of sksurv in the survival analysis setting (i.e. no competing events). Maybe @judithabk6 can add a little bit more context for our understanding?

[Vincent-Maladiere]

Actually, let's remove this section from the docstring because ties in time when D_j = 1 are not acceptable pairs, and we don't count them. We only use the ties in predictions.

[antoinebaker]

It could be nice to add the explicit math formula in the docstring for the API reference as was done for the Brier score.

[Vincent-Maladiere]

Thank you very much for this thorough review @antoinebaker! I updated the docstring and returned a numpy array instead of a list.

The two most important topics left to address are:

• The accounting of ties in the C-index, which deviates from [1]
• Passing an IPCW estimator or a string in the case of the KM

[juAlberge]

@Vincent-Maladiere For the difference is G(t-) = P(C >= t) and G(t) = P(C > t).
It only matters when T is discrete

[Vincent-Maladiere]

@Vincent-Maladiere For the difference is G(t-) = P(C >= t) and G(t) = P(C > t). It only matters when T is discrete

Yes, but the question is should we stay consistent with [1], where we have theoretical evidence that we're doing the right thing, or derive our own method?

[Vincent-Maladiere]

So, after discussing this with @judithabk6, the paper doesn't mention the ties because they are implementation details. The ties are present in the R version of this metric: https://github.com/cran/pec/blob/fb32746f6119450e8435deb52aa8583a91b05ba5/src/ccr.c#L105

So, I'll make sure to update the documentation accordingly.

[Vincent-Maladiere]

Let's also add the C-index in an example :)

[antoinebaker]

I'm still a bit confused about the tied in time events, I think there are some inconsistencies in the docstrings.

Pairs of type A are according to

1. the math formula for \tilde{A}_{ij} : only $T_i &lt; T_j$
2. the docstring of concordance_index_incidence, as well as the docstring of n_pairs_a in _concordance_index_incidence_report : $T_i &lt; T_j$ or $T_i = T_j, D_j = 0$
3. the docstring of n_ties_times_a in _concordance_index_incidence_report : possibly $T_i = T_j, D_j = 1$

I guess 2. is the correct one ? If so, what is the use of n_ties_times_a (it doesn't seem to appear in the c-index computation) ?

[antoinebaker]

The example is nice! The C-index (mathematical formulation) should be replaced by the C-index in time. Should we add also a version with IPWC ? [but maybe the example will be a bit confusing, with two unrelated KM estimators, one as a baseline to compare with SurvivalBoost, the other as IPWC to debias the C-index in time]

[Vincent-Maladiere]

Hey @antoinebaker, thanks for this additional feedback.

I'm still a bit confused about the tied in time events, I think there are some inconsistencies in the docstrings [...]

• Right, I updated 1. with the definition of 2.
• The definition of n_ties_times_a looks correct to me. Ties in times are pairs having both experienced the event of interest simultaneously. They are not comparable, so they don't contribute to the C-index. We only put them in the report for debugging purposes, akin to what scikit-survival does. This debugging mode, using _concordance_index_incidence_report, is not documented.

The example is nice! The C-index (mathematical formulation) should be replaced by the C-index in time. Should we add also a version with IPWC ? [but maybe the example will be a bit confusing, with two unrelated KM estimators, one as a baseline to compare with SurvivalBoost, the other as IPWC to debias the C-index in time]

I wonder if we should remove these math formulation entirely from the example, to only keep the ones in the function docstring. WDYT?

[Vincent-Maladiere]

Ping @antoinebaker :)

[antoinebaker]

I add a few comments to make the code easier to follow, but otherwise LGTM :)

[judithabk6]

it's probably defined elsewhere, but on standalone reading, might be useful to redefine $\tilde{T}$ and $\tilde{D}$.

[judithabk6]

the docstring should specify that this is the event-specific concordance index, to explicit why we distinguish the event and the competing events (maybe at the beginning)

[judithabk6]

those links are clickable, but refer to the same page, maybe they should link to the article?

[judithabk6]

maybe explicit why cox would be relevant here (just say that it would allow the censoring model to be conditional)?

[Vincent-Maladiere]

Thank you very much @judithabk6 and @antoinebaker for this thorough round of reviews! I added your suggestions and both the doc and the codebase are much clearer now :)