Introduction
Survival analysis allows us to examine the time to an event such as death. However, there are potential problems when the time interval between the exposure and outcome are improperly classified. One serious issue is immortal time bias where the exposure for one group requires them to be “alive” long enough to receive the treatment thereby giving them an advantage in terms of longer time at risk. For a succinct summary of immortal time bias, I recommend reading Samy Suissa’s excellent paper. Figure 1 illustrates the immortal time bias where the index date is set at the time of enrollment, but the group classificaiton occurs afterwards when the subject receive the treatment.
Figure 1. Immortal time bias.
Notice the time interval between enrollment and treatment; the subject must remain alive long enough to receive treatment. In survival analysis, this time interval may be improperly included as part of the “exposure” time or “at-risk” time, which can provide the treatment group with an advantage over the control group.
To address immortal time bias, it will be necessary to change the index date from the time of enrollment to the time when the subject receives the treatment. However, this would also mean that the time between enrollment to treatment receipt will somehow need to be accounted for.
Time-varying covariate
To address potential immortal time bias, it is recommended that the analysis includes a time-varying predictor. This is where the main predictor of interest (e.g., treatment receipt) is a function of time. Hence, when the subject is enrolled and does not receive treatment, they are considered part of the no-treatment (control) group. Once they receive the treatment, then the subject will be coded as the treatment group at that specific point in time. Figure 2 illustrates this concept.
Figure 2. Time-varying predictor; period before the treatment is considered part of the unexposed group.
To create a time-varying predictor, you will need to get the data
into the long format and create an variable where the treatment exposure
is a function of time. Figure 3 illustrates the time-varying variable
where the subject received the treatment at time = 4.
Figure 3. Time-varying predictor; subject received treatment at
time = 4.
Once you have the time-varying covariate, you can use this grouping variable (that is a function of time) in your survival model.
Motivating Example: Academy Award Winners versus Control - A study on survival
In 2021, a retrospective cohort study by Redelmeier and Singh that was publised in the Annals of Internal Medicine, investigated the differences in survival between actors/actresses who won and did not win an Academy Award. They concluded that actors/actress who won an Academy Award had a significantly greater survival compared to those who did not. This sparked some interesting discussion about the link between status and winning the Academy Award. So much so that Sylvestre and colleagues requested data from the original authors to redo the study. (Note: Sylvestre and colleagues report was published in 2022 in the Annals of Internal Medicine link.)
Sylvestre and colleagues were concerned with immortal time bias and sought to redo the analysis with this in mind. The only difference was the dataset. Redelmeier and Singh performed their analysis will data up until 2020. However, the data they provided to Sylvertre, et al, contained data up to 2021, which may result in values that are not exact with that of the orignal publications. However, Sylvestre and colleagues were able to replicate many of the original findings and adding some adjustments to address the immortal time bias in the original study.
Figure 4 is an excerpt from Sylvestre and colleagues findings where they compared their results with those of the original paper. Notice that the new data yieled a non-significant reduction in mortarlity when the time-varying predictor is applied.
Figure 4. Comparison of Proportional Hazards models between the old and new data.
We will recreate the results by Sylvestre and colleague using data they received from the original authors in Stata.
Step 1: Import the data
I saved the data on my GitHub
site, which you can importt using the Stata
import command.
import delimited "https://raw.githubusercontent.com/mbounthavong/Survival-analysis-and-immortal-time-bias/main/Data/data1.csv"Step 2: Prepare the data
Then, we will need to apply our inclusion criteria. We only want to compare the Winners to those who did not win the Academy Award. Therefore, we will drop the actors/actress who received a nomination but did not win an award.
// Drop nominees that didn't win (Choose this option to make the comparison between winners and controls)
drop if noms >= 1 & wins == 0Step 3: Transform data to long format
In order to generate the time-varying exposure term, we will need to change the data from wide to long format.
I used the following code, which I modified from the Stata FAQ blog.
// Generate a repeated measures file (This will allow for time-varying predictor [e.g., group])
expand time /* expands the time to include all missing value between start and end */
sort pid /* sort by patient identifier */
qui by pid: gen t = _n /* populate the values for the time "t" */
gen outcome = 0 /* create an outcome variable */
qui by pid: replace outcome = death if _n == _N /* indicator for the end point */. // Generate a repeated measures file (This will allow for time-varying predic
> tor [e.g., group])
. expand time /* expands the time to include all missing value between start an
> d end */
(74,465 observations created)
. sort pid /* sort by patient identifier */
. qui by pid: gen t = _n /* populate the values for the time "t" */
. gen outcome = 0 /* create an outcome variable */
. qui by pid: replace outcome = death if _n == _N /* indicator for the end poin
> t */
. Step 4: Perform the survival analysis without a time-varying predictor
Once we have the data set up in the long format, we can perform the survival analysis.
****************************************************
// Analysis for the Winners v. Controls:
****************************************************
*** STSET
stset t, failure(outcome) id(pid). // Analysis for the Winners v. Controls:
. ****************************************************
. *** STSET
. stset t, failure(outcome) id(pid)
Survival-time data settings
ID variable: pid
Failure event: outcome!=0 & outcome<.
Observed time interval: (t[_n-1], t]
Exit on or before: failure
--------------------------------------------------------------------------
75,606 total observations
0 exclusions
--------------------------------------------------------------------------
75,606 observations remaining, representing
1,141 subjects
565 failures in single-failure-per-subject data
75,606 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 104
. After setting data data for survival analysis, we can generate the Kaplan-Meier curves for the Winners and Controls.
// Survival analysis (without the time-varying group predictor)
*** KM curve and life expectency (without the time-varying group predictor)
sts graph, by(group) legend(label(1 "Controls") label(2 "Winners")) risktable(, order(1 "Controls" 2 "Winners")) /* kaplan meier curve */
quietly graph export kaplan1.svg, replace /* save graphic */. // Survival analysis (without the time-varying group predictor)
. *** KM curve and life expectency (without the time-varying group predictor)
. sts graph, by(group) legend(label(1 "Controls") label(2 "Winners")) risktable
> (, order(1 "Controls" 2 "Winners")) /* kaplan meier curve */
Failure _d: outcome
Analysis time _t: t
ID variable: pid
. quietly graph export kaplan1.svg, replace /* save graphic */
. 
We can also get the restricted mean survival time (e.g., life
expectancy) for each group. There are two ways to do this. The first
involves using the stci command with the
rmeans option. This is not very accurate from my
experience, and it doesn’t give you the statistical comparisons between
the two life expectancies Hence, I prefer to use the
strmst2 command. The strmst2 command will
provide a more accurate life expectancy and the statistical comparisons
between the groups. The life expectancy for the Winners is 79.66 years
and for the Controls is 75.96 years. The difference is 3.69 years, which
is statistically significant (P = 0.03). This is close to the value that
Sylvestre and colleagues reported in their re-analysis of the original
paper (see Figure 4).
stci, rmean by(group) /* Estimate the restricted mean survival time using rmean */
strmst2 group /* Estimate the restricted mean survival time using strmst2 package */. stci, rmean by(group) /* Estimate the restricted mean survival time using rm
> ean */
Failure _d: outcome
Analysis time _t: t
ID variable: pid
| Number of Restricted
group | subjects mean Std. err. [95% conf. interval]
-------------+-------------------------------------------------------------
0 | 902 76.05677(*) .5375206 75.0033 77.1103
1 | 239 79.66331 .9903624 77.7222 81.6044
-------------+-------------------------------------------------------------
Total | 1141 76.8056(*) .4763389 75.872 77.7392
(*) largest observed analysis time is censored, mean is underestimated
. strmst2 group /* Estimate the restricted mean survival time using strmst2 pac
> kage */
Number of observations for analysis = 75606
The truncation time, tau, was not specified. Thus, the default tau (the minimum
> of the
largest observed event time within each group), 100.000, is used.
Restricted Mean Survival Time (RMST) by arm
-----------------------------------------------------------
Group | Estimate Std. Err. [95% Conf. Interval]
---------+-------------------------------------------------
arm 1 | 79.663 1.089 77.529 81.798
arm 0 | 75.973 0.567 74.861 77.085
-----------------------------------------------------------
Between-group contrast (arm 1 versus arm 0)
------------------------------------------------------------------------
Contrast | Estimate [95% Conf. Interval] P>|z|
---------------------+--------------------------------------------------
RMST (arm 1 - arm 0) | 3.690 1.284 6.097 0.003
RMST (arm 1 / arm 0) | 1.049 1.017 1.081 0.002
------------------------------------------------------------------------
. We can perform the log-rank test, which is statistically significant (P = 0.004).
*** Log-rank test
sts test group, logrank detailNumber of observations for analysis = 75606
The truncation time, tau, was not specified. Thus, the default tau (the minimum
> of the
largest observed event time within each group), 100.000, is used.
. *** Log-rank test
. sts test group, logrank detail
note: option detail relevant only in combination with option strata().
Failure _d: outcome
Analysis time _t: t
ID variable: pid
Equality of survivor functions
Log-rank test
| Observed Expected
group | events events
------+-------------------------
0 | 461 433.06
1 | 104 131.94
------+-------------------------
Total | 565 565.00
chi2(1) = 8.17
Pr>chi2 = 0.0043
. Lastly, we construct the Cox Proportional Hazard model and compare our findings with those of Sylvestre and colleagues. In our results, the hazard ratio is 0.74 with a 95% confidence interval of 0.60 and 0.92. This translates to a mortality reduction of 26% (95% CI: 8% and 40%). This is exactly the same values that Sylvestre and colleagues reported in their paper (see Figure 4).
*** Cox PH model
stcox groupNumber of observations for analysis = 75606
The truncation time, tau, was not specified. Thus, the default tau (the minimum
> of the
largest observed event time within each group), 100.000, is used.
. *** Cox PH model
. stcox group
Failure _d: outcome
Analysis time _t: t
ID variable: pid
Iteration 0: Log likelihood = -3297.9291
Iteration 1: Log likelihood = -3293.8686
Iteration 2: Log likelihood = -3293.8436
Iteration 3: Log likelihood = -3293.8436
Refining estimates:
Iteration 0: Log likelihood = -3293.8436
Cox regression with Breslow method for ties
No. of subjects = 1,141 Number of obs = 75,606
No. of failures = 565
Time at risk = 75,606
LR chi2(1) = 8.17
Log likelihood = -3293.8436 Prob > chi2 = 0.0043
------------------------------------------------------------------------------
_t | Haz. ratio Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
group | .739926 .0804105 -2.77 0.006 .5979779 .9155699
------------------------------------------------------------------------------
. Step 5: Perform the survival analysis with a time-varying predictor
Now, we can perform the survival analysis using the time-varying
predictor. In our data, we created a time-varying group variable called
group_win.
// Generate the time-varying predictor
gen group_win = .
replace group_win = 0 if t < firstwin
replace group_win = 1 if t >= firstwin
tab group_win, m Number of observations for analysis = 75606
The truncation time, tau, was not specified. Thus, the default tau (the minimum
> of the
largest observed event time within each group), 100.000, is used.
. // Generate the time-varying predictor
. gen group_win = .
(75,606 missing values generated)
. replace group_win = 0 if t < firstwin
(69,095 real changes made)
. replace group_win = 1 if t >= firstwin
(6,511 real changes made)
. tab group_win, m
group_win | Freq. Percent Cum.
------------+-----------------------------------
0 | 69,095 91.39 91.39
1 | 6,511 8.61 100.00
------------+-----------------------------------
Total | 75,606 100.00
. Once we have this time-varying grouping variable, we can repeat the
exercises above but with the group_win instead of the
group variable. Essentially, we just replace the
group variable with the group_win
variable.
We can compare the Kaplan-Meier curves for the Winners and Controls using the time-varying predictor.
// Survival analysis (with time-varying group predictor)
*** KM curve and life expectency (with the time-varying group predictor)
sts graph, by(group_win) legend(label(1 "Controls") label(2 "Winners")) risktable(, order(1 "Controls" 2 "Winners")) /* kaplan meier curve */
quietly graph export kaplan2.svg, replace /* save graphic */Number of observations for analysis = 75606
The truncation time, tau, was not specified. Thus, the default tau (the minimum
> of the
largest observed event time within each group), 100.000, is used.
. // Survival analysis (with time-varying group predictor)
. *** KM curve and life expectency (with the time-varying group predictor)
. sts graph, by(group_win) legend(label(1 "Controls") label(2 "Winners")) riskt
> able(, order(1 "Controls" 2 "Winners")) /* kaplan meier curve */
Failure _d: outcome
Analysis time _t: t
ID variable: pid
. quietly graph export kaplan2.svg, replace /* save graphic */
. 
The life expectancy of the Winners and Controls are 78.77 and 76.31 years, respectively. This is not statistically significant (P = 0.056).
stci, rmean by(group_win) /* Estimate the restricted mean survival time using rmean */
strmst2 group_win /* Estimate the restricted mean survival time using strmst2 package */Number of observations for analysis = 75606
The truncation time, tau, was not specified. Thus, the default tau (the minimum
> of the
largest observed event time within each group), 100.000, is used.
. stci, rmean by(group_win) /* Estimate the restricted mean survival time usin
> g rmean */
Failure _d: outcome
Analysis time _t: t
ID variable: pid
| Number of Restricted
group_win | subjects mean Std. err. [95% conf. interval]
-------------+-------------------------------------------------------------
0 | 1141 76.39926(*) .5261978 75.3679 77.4306
1 | 239 69.76821 1.063185 67.6844 71.852
-------------+-------------------------------------------------------------
Total | 1141 76.8056(*) .4763389 75.872 77.7392
(*) largest observed analysis time is censored, mean is underestimated
. strmst2 group_win /* Estimate the restricted mean survival time using strmst2
> package */
Number of observations for analysis = 75606
The truncation time, tau, was not specified. Thus, the default tau (the minimum
> of the
largest observed event time within each group), 100.000, is used.
Restricted Mean Survival Time (RMST) by arm
-----------------------------------------------------------
Group | Estimate Std. Err. [95% Conf. Interval]
---------+-------------------------------------------------
arm 1 | 78.768 1.157 76.500 81.037
arm 0 | 76.314 0.556 75.224 77.403
-----------------------------------------------------------
Between-group contrast (arm 1 versus arm 0)
------------------------------------------------------------------------
Contrast | Estimate [95% Conf. Interval] P>|z|
---------------------+--------------------------------------------------
RMST (arm 1 - arm 0) | 2.455 -0.062 4.971 0.056
RMST (arm 1 / arm 0) | 1.032 1.000 1.066 0.054
------------------------------------------------------------------------
. The log-rank test is not statistically significant (P = 0.07).
*** Log-rank test
sts test group_win, logrank detailNumber of observations for analysis = 75606
The truncation time, tau, was not specified. Thus, the default tau (the minimum
> of the
largest observed event time within each group), 100.000, is used.
Number of observations for analysis = 75606
The truncation time, tau, was not specified. Thus, the default tau (the minimum
> of the
largest observed event time within each group), 100.000, is used.
. *** Log-rank test
. sts test group_win, logrank detail
note: option detail relevant only in combination with option strata().
Failure _d: outcome
Analysis time _t: t
ID variable: pid
Equality of survivor functions
Log-rank test
| Observed Expected
group_win | events events
----------+-------------------------
0 | 461 443.67
1 | 104 121.33
----------+-------------------------
Total | 565 565.00
chi2(1) = 3.37
Pr>chi2 = 0.0666
. Lastly, we construct the Cox Proportional Hazard model and compare our findings with those of Sylvestre and colleagues. In our results, the hazard ratio is 0.82 with a 95% confidence interval of 0.67 and 1.02. This translates to a mortality reduction of 18% (95% CI: -2% and 3%). This is almost the same values that Sylvestre and colleagues reported in their paper (see Figure 4). The point estimate is off by 1%-point, but I believe this could be due to a rounding error.
Unlike the original paper, the reduction in mortality rate is not statistically significant when the grouping variable is time-varying. This indicates that immortal time bias may have resulted in a false positive or type I error.
*** Cox PH model
stcox group_winNumber of observations for analysis = 75606
The truncation time, tau, was not specified. Thus, the default tau (the minimum
> of the
largest observed event time within each group), 100.000, is used.
Number of observations for analysis = 75606
The truncation time, tau, was not specified. Thus, the default tau (the minimum
> of the
largest observed event time within each group), 100.000, is used.
. *** Cox PH model
. stcox group_win
Failure _d: outcome
Analysis time _t: t
ID variable: pid
Iteration 0: Log likelihood = -3297.9291
Iteration 1: Log likelihood = -3296.285
Iteration 2: Log likelihood = -3296.2803
Iteration 3: Log likelihood = -3296.2803
Refining estimates:
Iteration 0: Log likelihood = -3296.2803
Cox regression with Breslow method for ties
No. of subjects = 1,141 Number of obs = 75,606
No. of failures = 565
Time at risk = 75,606
LR chi2(1) = 3.30
Log likelihood = -3296.2803 Prob > chi2 = 0.0694
------------------------------------------------------------------------------
_t | Haz. ratio Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
group_win | .8236429 .0897609 -1.78 0.075 .6652344 1.019772
------------------------------------------------------------------------------
. Here is a comparison of our tutorial findings with those of Redelmeier & Singh and Sylvertre and colleagues.
Figure 7. Comparison of mortality reduction between the analyses.
Conclusions
Address immortal time bias in surival analysis is possible with a time-varying predictor. In this tutorial, we review how the immortal time bias can impact the association between Academy Award Winners and Controls. With the time-varying predictor, the association is non-significant due to the amount of “immortal time” the subject in the Winners group is allowed. However, when we assign this time interval to the Control group, the differences in mortality rates is attentuated.
Acknowledgement
Redelmeier and Singh’s study is avaialble at the Annals of Internal Medicine link.
In 2020, they redid this analysis with updated data and published their findings in Plos One link.
Sylvestre and colleagues re-did the 2001 analysis by requesting the data from the original authors in 2022. link. The data from Sylvestre and colleagues were used for this tutorial.
William Gould from StataCorp wrote a Stata Blog on creating a time-varying covariate for use in a survival analysis, which was incredibly helpful link.
Disclosures and Disclaimers
This is a work in progress, and I will likely update this as I learn more efficient methods to produce 95% CI around average values.
This is for educational purposes only.
The Github code and data for this R Markdown tutorial is located here.