How to address 'immortal time bias' using R - equivalent to Stata stset?

Question

We have a dataset on cancer patients who have consented to join a study after their diagnosis, which could be months or even years later. After some follow-up, an event occurs. We can fit this data using the following code in R:

fit <- survfit(Surv(time, status) ~ group, data = data)

Here, the variable time represents the time between the event and the diagnosis, while status indicates whether the event occurred or not.

However, it is important to note that any patient who enters the study must have survived from the time of diagnosis until the time of entry. If they had died during this period, they would not have been included in the study, and we would not have any information about them. This phenomenon is known as "immortal time bias."

To avoid this bias, we need to define the time period during which we are "watching" each subject. In other words, we need to set the origin as the date of diagnosis and the entry time as the date of consent. In Stata, this can be done using the stset command with the enter() and origin() options. What is the equivalent in R?

In short, in the stset command Stata has the options, 'start', 'stop' and 'origin'. We have a cohort where the start dates come after the origin date. What is the equivalent in R?

Answer 1

You are correct that this needs to be taken into account. Anyone entering the study at some time after diagnosis provides no information about what might have happened during the intervening time. If diagnosis date is the time=0 reference for the survival model, then time to study entry needs to be treated as a left truncation of survival time.

This situation can be handled with the counting-process data format for survival. The simplest way to proceed is to set time=0 for each individual to the date of diagnosis and express other times relative to that time. That sets the time origin for the observations. If calendar date of diagnosis is a potential predictor, that calendar date can be coded separately.

Then you can specify the outcome variable for the regression model as Surv(startTime, stopTime, event). That is interpreted as left truncation at startTime and either right censoring or an event at stopTime, depending on the value of event. This counting-process data format also simplifies many other types of survival analysis, as explained in Section 3.7 of Therneau and Grambsch. Section 3.7.3 in particular discusses a situation similar to what you describe.

Answer 2

This is an interesting problem. As an R question it would, of course, be off topic for the site. You already have one solution from EdM regarding setting the starttime and stoptime variables in the Surv command: this effectively defines the time-dependent risk set treating diagnosis date as Day 0 or Time 0. What you are modeling with this approach is a flexible non-parametric baseline hazard function as a function of time from diagnosis and not time from entry into study.

The problem is that there are many "times" which can be predictive of outcomes: there is for instance the patient's age, or the calendar year, or time from entry into study (either initiation of treatment or signing informed consent). These are broadly referred to as the age, period, and cohort effects.

You might be surprised to learn left truncation as a solution to immortal time bias is less common than you expect. For instance, in every cancer study I have supported, not one has modeled event time from diagnosis or time from last treatment. The truth is, if the "gap" exists for every patient in the sample, the associated Nelson-Aelan curve is very misleading because the incidence is, by definition, 0 per person-time from time 0 until the first observed failure. Note: you can't even estimate a Kaplan Meier because not all subjects are at risk (for a measurable event) at time 0. Adjusting the start and stop times affects whom you compare to whom, but this approach still cannot account for subjects who died prior to inclusion in the sample. If you are lucky enough to sample a few subjects in the early time period, it's important to handle these events appropriately, and in this case, EdM's approach is a valid --but not the only-- approach to get a stable estimate of incidence over time. The approach should also be justified by the proportionality of the hazard ratios as assessed by plots and tests if need be.

I actually disagree with the way most time-to-event analyses are conducted for cancer studies, but not for that reason! As you know, a Cox model is a very flexible model in that it is semi-parametric in how it handles an arbitrary baseline hazard function (as a default or in multiple groups via stratification), and it handles censoring. A cox model also allows for adjustments like in a regression model. I actually believe time from entry into study should be the time 0 because this time is the most "artificially" constructed time; that is to say, when generalizing the results there is no need to "predict" time from entry into study for a patient (not subject) who's not actually in a study. For generalizable times like time from diagnosis, calendar year, etc. you can simply adjust for these as variables in the model. Depending on the sample size, you can use more and more sophisticated approaches like splines to estimate other time-dependent incidence parameters.

Stephen Senn has written extensively on the benefits of adjusting for prognostic variables, even in randomized ITT analyses. I agree with him and would suggest, depending on the analysis, and the nature of the functional shape of the hazard function, entry-into-study is a reasonable choice of time 0, and other predictive "times" can be handled as variables in the model. Note: these would not be time varying functions since, from study entry, age, cohort, and period all move in unison, so the typical adjustments would be "age at study entry" or "years since initial diagnosis at study entry" or "drug-free interval" all measured in clinically relevant units.