Does lifelines have a competitive risk model for survival analysis

[Abraxas2071]

As title show

[pzivich]

lifelines has two different options for competing risk.

1) AalenJohansenFitter

The Aalen-Johansen (AJ) estimator is a non-parametric alternative to the Kaplan-Meier in scenarios with competing risks. AJ estimates the cause-specific cumulative incidence function. The online docs currently don't have an example (my fault for not creating one), but below is an example (with v0.17.5). First, I will create a simple data set, where 'events' indicates the type of event. In the example, 1 indicates the event of interest and 2 indicates the competing risk

import pandas as pd
from lifelines import AalenJohansenFitter

df = pd.DataFrame()
df['events'] = [0, 1, 1, 1, 2, 1, 0, 1, 0, 2, 2]
df['times'] = [1, 5, 3, 6, 1, 3, 2, 5, 2, 3, 6]

We can then fit the AJ estimator using

aj = AalenJohansenFitter()
aj.fit(durations=df['times'], event_observed=df['events'], event_of_interest=1)
print(aj.cumulative_density_)

where the cumulative incidence function is estimated for event type . Note that in my example a warning will be thrown. This is because the AJ estimator cannot deal with tied event times. In the background, AalenJohansenFitter randomly pushes the event times around by a small amount to break ties (by 0.0001). We can compare the above to a Kaplan-Meier with competing events as censored

df['events'] = np.where(df['events']==1, 1, 0)
km = KaplanMeierFitter()
km.fit(durations=df['times'], event_observed=df['events'])

Below are the different results

          KM_estimate     CIF_1
0.000000     0.000000  0.000000
1.000000     0.000000  0.000000
1.000021     0.000000  0.000000
1.999944     0.000000  0.100000
2.000000     0.111111  0.100000
3.000000     0.259259  0.100000
3.000001     0.259259  0.233333
3.000036     0.259259  0.233333
4.999919     0.259259  0.366667
5.000000     0.629630  0.366667
5.000075     0.629630  0.500000
5.999954     0.629630  0.500000
6.000000     0.814815  0.500000
6.000074     0.814815  0.633333

So the Kaplan-Meier overestimates the cumulative indcidence versus the Aalen-Johansen

Disadvantage: no covariate adjustment or estimation is possible. If you are interested in a single treatment, you can use inverse probability weights to adjust for confounders (Cole Hernan 2004)

2) Cox Proportional Hazard

You can also use what is called the subdistribution hazard model in the Cox Proportional Hazard Model. lifelines can do this because no special functions are needed. For the subdistribution hazards CoxPH model, all you need to do is recode the competing event as censored then estimate the CoxPH model.

Why this works: this approach works since the CoxPH likelihood function is based on the proportional hazards, i.e. it does not directly estimate the shape of the hazard function. So, by treating the competing events as censored, you can estimate the proportionality of the sub-distribution (event without the competing event) hazards.

Problems: This approach only allows you to make claims about the hazard ratios from the model. You can no longer estimate cumulative incidence, survival, etc. from the CoxPH through this approach. This is because you are having CoxPH model estimate the actual hazard function, which no longer is correctly estimated.

Quoting from Zhang 2017

The effect of covariates on hazard rate can be estimated using COX proportional hazard regression model. The exponentiation of the coefficient gives the hazard ratio (HR) which is the ratio of rate in epidemiology. Because hazard function has a one-to-one correspondence to the cumulative incidence, the HR also reflects the risk (cumulative incidence) of the study population. However, this relationship does not exist in the presence of competing risks. Although the cause-specific hazard regression model represents the impact of covariates on the cause-specific hazard, it does not necessarily reflect the impact on the cumulative incidence.

3) ???

I don' believe that lifelines has any other approaches to deal with competing events. Other options that might be worth considering in the future are Fine and Gray models. There is also a variation of the CoxPH model that allows for competing risks, but I can't remember the terminology currently

[CamDavidsonPilon]

Thanks @pzivich for that write-up, TIL about using cox models for competing risks.

Like Paul said, there's sparse support of competing risks in lifelines aside from the two solutions he mentioned.

Searching Fine and Gray in google, the first article displayed is a cautionary tale of using it and a promotion of point 2. above: https://statisticalhorizons.com/for-causal-analysis-of-competing-risks (though they use opposite nomenclature... maybe I confused myself though)

[pzivich]

Maybe outside this discussion, but the following quote from the blog post is misleading

Unfortunately, there’s no way to test the assumption that censoring is non-informative (Tsiatis 1975). And even if there were, there’s no good method available for estimating causal effects when censoring is informative

There are methods to account for informative censoring. Inverse probability of censoring weights (IPCW) were made to account for informative censoring. You do so under that assumption that censoring is noninformative conditional on the covariates included in your IPCW model. While referred to as monotone missing data in Sun & Tchetgen Tchetgen 2018, section 3 details these weights briefly (don't have a better source on hand).

When Allison says,

However, as I now show, the subdistribution method does no better (and actually somewhat worse) than the cause-specific method when competing events are informative.

I agree, but I would word it slightly differently. For causal analysis in competing risk scenarios, you need to adjust for both the confounders of the treatment-event of interest relationship but also the treatment-competing risk relationship. At least this is how it has been explained to me

EDIT: this paper explicitly demonstrates the equivalence of Kaplan-Meier and IPCW estimators under non-informative censoring Satten & Datta 2001. They also refer to the 1992 paper where Robins and Rotnitzky propose IPCW

[Abraxas2071]

@pzivich Our Dear Doc Pzivich~Thank you very much for your reply and I learned a lot from your answer But I still have some questions： 1】“AalenJohansenFitter()” can be used to predict？ 2】“all you need to do is recode the competing event as censored then estimate the CoxPH model”。This seems to be similar to the way in which competing events are treated as censored in KM estimation。According to your first point, KM estimates overestimate the incidence because they treat competing events as censored .But is it really feasible for COX to do the same?I don't quite understand 3】I don't know what “predict_partial_hazard” means. What does that mean “partial hazard”，link：https://lifelines.readthedocs.io/en/latest/Survival%20Regression.html#cox-s-proportional-hazard-model ——Because I am a Chinese, English is not my strong point

PS：What is "tied event times"?

[Abraxas2071]

@CamDavidsonPilon Our Dear Cameron~Nice to meet you！Has anyone ever said you look like Jake Gyllenhaal？ I have the Chinese version of your book and I noticed your relationship to bayesian analysis。 Then，How is Lifelines the package combined with pymc3 for bayesian survival analysis？

[pzivich]

@Abraxas2071 Happy to help!

1】“AalenJohansenFitter()” can be used to predict？

You can predict the cumulative incidence function [F(t)] at t much like how you would predict survival times from a Kaplan Meier curve. As for AalenJohansenFitter() I did not create a predict() functionality. The way around this is to look at the cumulative_incidence_ output for the time you are interested in.

An aside, AalenJohansenFitter should not be used to predict survival [S(t)]. S(t) is not defined well in competing risk scenarios. For example, think about influenza infection with death as a competing risk. The incidence of influenza (accounting for death) makes sense. However, the survival of avoiding influenza but not dying is not clear. This problem for S(t) extends to all competing risk estimation methods.

If you want to predict survival, the other option is to create a combined outcome. You could create the event influenza infection or death. S(t) is makes more sense in this scenario. However, you interpretation changes to influenza infection or death. It is no longer for each event specifically. Combined outcomes would be compatible with any survival analysis method

2】“all you need to do is recode the competing event as censored then estimate the CoxPH model”。This seems to be similar to the way in which competing events are treated as censored in KM estimation。

Correct. It is, but unlike the Kaplan-Meier, this trick in the CoxPH model is valid. It is valid because of some underlying mathematics regarding the partial likelihood function

According to your first point, KM estimates overestimate the incidence because they treat competing events as censored .But is it really feasible for COX to do the same?I don't quite understand

Kaplan-Meier will overestimate the incidence. Similarly CoxPH will as well, if you use it to predict incidence. This trick for the CoxPH will not allow you to predict specific values at times (like hazards or incidence). The only valid measure will be the hazard ratios from the CoxPH model. The hazard ratios will describe the relation between the partial hazards, but not the particular value of the hazards. Going back to the influenza example, imagine you put age into the CoxPH model. The estimated hazard ratio is 1.3 (HR = 1.3). This means that the hazard of influenza is 1.3 times higher for every year increase in age. I cannot tell you the hazard at a specific time with this method.

3】I don't know what “predict_partial_hazard” means. What does that mean “partial hazard”，link：https://lifelines.readthedocs.io/en/latest/Survival%20Regression.html#cox-s-proportional-hazard-model

See previous. lifelines does not have a predict_partial_hazard() because the trick I mentioned above does not give you particular hazard values. Rather it describes the relative difference between hazards by some variable. It is a confusing trick and took some time to make sense to me

Because I am a Chinese, English is not my strong point

You're doing great! I hope my reply is helpful and clear

PS：What is "tied event times"?

Tied event times are when at least two individuals have events occur at the same time in your data. If person 1 has one of the outcomes at t=5 and person 2 has one of the outcomes at t=5, then your data has a tied event times at t=5

To summarize you have several options to deal with competing risks depending on your specific scenario:

1) Estimate cumulative incidence with Aalen-Johansen. This approach will allow you to estimate the cumulative incidence. However, it is cumbersome if you want to predict values for a bunch of variable combinations. You would need to stratify your data by the variables and estimate Aalen-Johansen separately. Worse, it cannot handle continuous variables you might be interested in as predicting incidence

2) Estimate hazard ratios with CoxPH trick. You can estimate hazard ratios to describe what is correlated to increased/decreased hazards. You cannot validly estimate the specific hazards/incidence with this approach

3) Create a combined outcome. This will allow you to use standard estimation methods and predict the hazard/survival/incidence. The issue is that your interpretation will need to change to reflect the combined events.

4) Fine-Gray model. This model is not in lifelines currently but can be used to predict incidence. I am unaware of any Python implementations of Fine-Gray

5) Treat all competing events as censoring for any estimation method. I would not recommend this generally. However, if there are few competing events occurring in your data, then the bias in your estimates will be small. The proportion of people experience the competing event should be less than 5% for this option

@CamDavidsonPilon I will work on writing up a description of different competing risk approaches for ReadTheDocs sometime in the next week. I will go over the five options above and the various advantages/disadvantages

[Abraxas2071]

@pzivich Our Dear Doc Pzivich~I have some minor problems： 1】What does the "create a combined outcome"？Is it the combined probability of influenza infection and accidental death? 2】How to operate “create a combined outcome” in Lifelines？ 3】“Estimate hazard ratios with CoxPH trick”。Can the independent variable be continuous? 4】“This means that the hazard of influenza is 1.3 times higher for every year increase in age”。I learned from wikipedia that HR is a comparison between two groups, but in your case it is a similar interpretation of the coefficient in regression analysis. What should I make of it? ——I only have one group at the moment 5】https://en.wikipedia.org/wiki/Hazard_ratio。This link said " When a study reports one hazard ratio per time period, it is assumed that difference between groups was proportional. Hazard ratios become meaningless when this assumption of proportionality is not met” ——Does this mean that accelerated failure models are needed?

Thank you very much~ヾ(ﾟ∀ﾟゞ)

[pzivich]

1】What does the "create a combined outcome"？Is it the combined probability of influenza infection and accidental death? 2】How to operate “create a combined outcome” in Lifelines？

Combined outcome would be recoding your data like the following example. You have a column in your data set labelled as Y, which has the outcomes. Using the previous example, 0 indicates censoring, 1 indicates influenza infection, and 2 indicates death. To create a combined column, you would create a new Y (Y_new), where 0 is censoring, and 1 is either influenza infetion or death. You can use numpy.where() to do this

3】“Estimate hazard ratios with CoxPH trick”。Can the independent variable be continuous?

For CoxPH the independent variable can be continuous

4】“This means that the hazard of influenza is 1.3 times higher for every year increase in age”。I learned from wikipedia that HR is a comparison between two groups, but in your case it is a similar interpretation of the coefficient in regression analysis. What should I make of it?

The HR compares two groups. It does have a similar interpretation to regression. CoxPH is a type of regression model for survival data

5】This link said " When a study reports one hazard ratio per time period, it is assumed that difference between groups was proportional. Hazard ratios become meaningless when this assumption of proportionality is not met” ——Does this mean that accelerated failure models are needed?

lifelines has the ability to determine whether the proportional hazards assumption is reasonable. Here is a great detailed example. When the assumption is violated, accelerated failure time (AFT) models are one option. As far as I know, lifelines does not have AFT models yet

[Abraxas2071]

@pzivich I am very glad to get your serious and professional reply

This link：https://lifelines.readthedocs.io/en/latest/Survival%20Regression.html#cox-s-proportional-hazard-model In this section【Prediction】,said：“After fitting, you can use use the suite of prediction methods: .predict_partial_hazard, .predict_survival_function, etc” The data of this document is about crime, there is only one group, there is no control group, but it can still be used to get “partial_hazard” Our Dear Doc Pzivich~That's where my question comes from. How do I understand it?

[pzivich]

Yes, so the example you linked is for a setting without competing risks. While you can use CoxPH model to calculate the partial hazard, it will be incorrect when competing risks occur. To get the partial hazard with competing risks, you will need to use a model like Fine-Gray, which is not available in lifelines

[Abraxas2071]

@pzivich Our Dear Doc Pzivich~It's so late and you are still replying. You have worked so hard

3】I don't know what “predict_partial_hazard” means. What does that mean “partial hazard”，link：https://lifelines.readthedocs.io/en/latest/Survival%20Regression.html#cox-s-proportional-hazard-model

You see, I don't quite understand what the partial hazard is。The formula in the documentation seems to say relative to the average of the covariates, right?

[CamDavidsonPilon]

The formula in the documentation seems to say relative to the average of the covariates, right?

Yea that's correct. The partial hazard excludes the baseline hazard (hence "partial"), and centers all the covariates by the average.

If using predict_partial_hazard, the interpretation is different from standard prediction models. For example, a linear regression prediction will give you the model's best estimate to E[Y], and one can use this directly in prediction applications. On the other hand, predict_partial_hazard gives you values that are really only interpretable against each other prediction, that is, the ratio of two predictions is important, not their absolute value. For example:

Subject 1 has predicted partial hazard 2.5, and Subject 1 has predicted partial hazard 0.5. We can say: 1) Subject 1 has larger hazard => more likely to die before Subject 2 2) Subject 1 has 5 times larger hazard than Subject 2.

[Abraxas2071]

The formula in the documentation seems to say relative to the average of the covariates, right?

Yea that's correct. The partial hazard excludes the baseline hazard (hence "partial"), and centers all the covariates by the average.

If using predict_partial_hazard, the interpretation is different from standard prediction models. For example, a linear regression prediction will give you the model's best estimate to E[Y], and one can use this directly in prediction applications. On the other hand, predict_partial_hazard gives you values that are really only interpretable against each other prediction, that is, the ratio of two predictions is important, not their absolute value. For example:

Subject 1 has predicted partial hazard 2.5, and Subject 1 has predicted partial hazard 0.5. We can say:

1. Subject 1 has larger hazard => more likely to die before Subject 2
2. Subject 1 has 5 times larger hazard than Subject 2.

@CamDavidsonPilon That's great！Our Dear Cameron！Thank you for your reply, which is the same as what I guessed

[Abraxas2071]

@CamDavidsonPilon @pzivich Our Dear Cameron and Our Dear Doc Pzivich I again come！ “To get the partial hazard with competing risks, you will need to use a model like Fine-Gray” It seems that in order to calculate the partial hazard, I need to use fine-gray. I plan to implement it in Python by myself. Is there any relevant information that can be provided to me?

[cvesteghem]

@pzivich Dear Paul,

Thanks a lot for your addition to support competing risks. I was trying to use the Aalen Johansen estimator with lifelines. You said that this estimator does not support tied events, and then your implementation randomly shifts the time of the events around the actual values to avoid the problem of concomitant events. I have nevertheless a rather large data set, notably with a lot of events at t=0, and I can see that some of the shifted times are actually identical. I was wondering if that could be a problem for the interpretation of the results.

Regards,

Charles

[pzivich]

So Aalen-Johansen does not support tied event times for events of different types. My apologies for not clearly explaining that in this thread. As long as the events are of the same type (for example j=1 for all events at that time=t) then the estimator is fine with tied events. The problem occurs when events of different types are recorded as happening at the same time. For the estimator, one of the event types must occur before the other. The random jitter process should resolve any tied event times of different types

My question would be is how are events occurring at t=0? Shouldn't the events occur some time after the origin?

[cvesteghem]

Thanks for the clarification (and the very fast answer). Then it's most probably less a problem in my case. The time unit I'm using is the day and I have some patients who dies on the very day they are included, I therefore count them as having an event on day 0.

[pzivich]

Ahhh I understand. I don't think Aalen-Johansen should have an issue with that in lifelines. However, all the t=0 events, I would update them to something like t=0.001 just to check.

You can also manually break ties in your events, if you would rather have some decision process rather than leaving it to chance (since it sounds like your time units are days). However, it may change the inference. The random jitter process is best for small time units. For example of a decision process, it may be reasonable to have all deaths occur before competing events each day. This will be different from jitter, which will have approx half of all deaths occur before the competing event

[cvesteghem]

Ok, Indeed, I could add "manual" shifting according to the type of event. Then I can make sure they can not be concomitant. Thanks for the advice.

[pzivich]

No problem! One more thought I had, I would try both versions of extreme manual shifts; (1) all deaths occur right before competing events, and (2) all deaths occur right after competing events. Since the persons-at-risk will differ under the two options, it may be important to see how it influences your results. It will give you an idea how sensitive your results are to the time-ordering of events. With sufficiently fine-grained time data (like maybe to the minute), it should make no difference, but it might for larger time units.

The random jitter process would be somewhere in-between the two above extremes

[cvesteghem]

That's a good idea.

By the way, I was thinking it could be nice to have an option for the fitter, for example:

ajf = AalenJohansenFitter()
ajf.fit(
    durations=df['durations'],
    event_observed=df['events'],
    event_of_interest=1,
    shift=[1,2]
)

where shift can be a list to know in which order to shift the corresponding events or a string for a specific method, for example "jitter", "none" ... At least it would make the test you are mentioning easier to implement ;)

[arinbasu]

Excellent discussion!

There seems to be a cmprsk fork of Fine and Gray's approach to competing risk model, here:

https://pypi.org/project/cmprsk/

From the webpage,

Regression modeling of sub-distribution functions in competing risks. A python wrapper around the cmprsk R package. Description: Estimation, testing and regression modeling of subdistribution functions in competing risks, as described in Gray (1988), A class of K-sample tests for comparing the cumulative incidence of a competing risk, Ann. Stat. 16:1141-1154, and Fine JP and Gray RJ (1999), A proportional hazards model for the subdistribution of a competing risk, JASA, 94:496-509.

[calvin-zcx]

lifelines has two different options for competing risk.

1. AalenJohansenFitter

2. Cox Proportional Hazard

You can also use what is called the subdistribution hazard model in the Cox Proportional Hazard Model. lifelines can do this because no special functions are needed. For the subdistribution hazards CoxPH model, all you need to do is recode the competing event as censored then estimate the CoxPH model.

Why this works: this approach works since the CoxPH likelihood function is based on the proportional hazards, i.e. it does not directly estimate the shape of the hazard function. So, by treating the competing events as censored, you can estimate the proportionality of the sub-distribution (event without the competing event) hazards.

Problems: This approach only allows you to make claims about the hazard ratios from the model. You can no longer estimate cumulative incidence, survival, etc. from the CoxPH through this approach. This is because you are having CoxPH model estimate the actual hazard function, which no longer is correctly estimated.

Quoting from Zhang 2017

The effect of covariates on hazard rate can be estimated using COX proportional hazard regression model. The exponentiation of the coefficient gives the hazard ratio (HR) which is the ratio of rate in epidemiology. Because hazard function has a one-to-one correspondence to the cumulative incidence, the HR also reflects the risk (cumulative incidence) of the study population. However, this relationship does not exist in the presence of competing risks. Although the cause-specific hazard regression model represents the impact of covariates on the cause-specific hazard, it does not necessarily reflect the impact on the cumulative incidence.

Regarding Point 2:
" For the subdistribution hazards CoxPH model, all you need to do is recode the competing event as censored then estimate the CoxPH model." I thinks this recoding approach leads to a CAUSAE-specific hazard model instead of a sub-distribution hazard model. See Ref. Austin, Peter C., Douglas S. Lee, and Jason P. Fine. "Introduction to the analysis of survival data in the presence of competing risks." Circulation 133, no. 6 (2016): 601-609. Section-statistical software for competing risks analyses. "Cause-specific hazard models can be fit in any statistical software package that permits estimation of the conventional Cox proportional hazards model. One simply treats those subjects who experience a competing event as being censored at the time of the occurrence of the competing event. In R, one can use the coxph function in the survival package, in SAS, one can use PROC PHREG, and in Stata, one can use the stcox function."