PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of amjepidLink to Publisher's site
 
Am J Epidemiol. 2016 August 15; 184(4): 325–335.
Published online 2016 July 25. doi:  10.1093/aje/kwv445
PMCID: PMC4983651

Comparison of Statistical Approaches for Dealing With Immortal Time Bias in Drug Effectiveness Studies

Mohammad Ehsanul Karim,* Paul Gustafson, John Petkau, and Helen Tremlett, the Long-Term Benefits and Adverse Effects of Beta-Interferon for Multiple Sclerosis (BeAMS) Study Group
the Long-Term Benefits and Adverse Effects of Beta-Interferon for Multiple Sclerosis (BeAMS) Study Group, Mohammad Ehsanul Karim, Paul Gustafson, John Petkau, Helen Tremlett, Afsaneh Shirani, Yinshan Zhao, Charity Evans, Elaine Kingwell, Mia L. van der Kop, Joel Oger, A. Traboulsee, A.-L. Sayao, V. Devonshire, S. Hashimoto, J. Hooge, L. Kastrukoff, J. Oger, D. Adams, D. Craig, S. Meckling, L. Daly, O. Hrebicek, D. Parton, and K. Pope

Abstract

In time-to-event analyses of observational studies of drug effectiveness, incorrect handling of the period between cohort entry and first treatment exposure during follow-up may result in immortal time bias. This bias can be eliminated by acknowledging a change in treatment exposure status with time-dependent analyses, such as fitting a time-dependent Cox model. The prescription time-distribution matching (PTDM) method has been proposed as a simpler approach for controlling immortal time bias. Using simulation studies and theoretical quantification of bias, we compared the performance of the PTDM approach with that of the time-dependent Cox model in the presence of immortal time. Both assessments revealed that the PTDM approach did not adequately address immortal time bias. Based on our simulation results, another recently proposed observational data analysis technique, the sequential Cox approach, was found to be more useful than the PTDM approach (Cox: bias = −0.002, mean squared error = 0.025; PTDM: bias = −1.411, mean squared error = 2.011). We applied these approaches to investigate the association of β-interferon treatment with delaying disability progression in a multiple sclerosis cohort in British Columbia, Canada (Long-Term Benefits and Adverse Effects of Beta-Interferon for Multiple Sclerosis (BeAMS) Study, 1995–2008).

Keywords: bias (epidemiology), confounding factors (epidemiology), epidemiologic methods, immortal time bias, longitudinal studies, models, survival analysis

In many observational studies of drug effectiveness, there may be a delay or waiting period before a subject begins to receive a treatment. Let us use multiple sclerosis (MS) as an example, with the disease-modifying drug β-interferon being the study drug of interest. If treated subjects entered a cohort at the time of their first β-interferon administration but untreated subjects entered the cohort at an earlier time, such as their first assessment at an MS clinic or study center conducting the observational study (as in the paper by Trojano et al. (1)), then the time period between the first clinic visit (baseline) and the first administration of β-interferon corresponds to “immortal time” (2, 3). This delay period is labeled “immortal” because, by design, subjects in the treated group cannot develop an outcome before receiving the treatment. Therefore, if some subjects develop the outcome before getting the opportunity to initiate treatment, they will be assigned to the untreated group.

If treatment status is misclassified during this delay period or if this delay period is excluded from the analysis, “immortal time bias” is introduced. Classifying subjects as “treated” or “untreated” with a time-invariant variable results in a spurious survival advantage (protective association) in favor of the treated group. This can considerably distort the estimated hazard ratio if a large number of failures (events or reaching outcome) occur before the initiation of treatment or if the duration of the immortal time is large (4). Although this bias was first identified in the 1970s (5), many studies still fail to account for this source of bias (2, 6, 7). Recently, the issue of immortal time bias has resurfaced in pharmacoepidemiology studies in the context of propensity score adjustment (3, 811), where treatment status must be defined at baseline.

The most common approach to accounting for immortal time bias is the adoption of a proper treatment exposure definition via time-dependent analyses (12), such as time-dependent Cox proportional hazards models (1316). Prescription time-distribution matching (PTDM) (13) is a simple approach to adjusting for immortal time bias that has been frequently cited (1721) and used (2228) in recent years. The performance and statistical efficiency of this approach have been claimed to be similar to those of the Cox model with a binary time-dependent exposure (hereafter called the “time-dependent Cox model”) (13).

Several studies have quantified the amount and direction of bias due to “misclassifying” or “ignoring” immortal time by means of simulation (14, 2931) or theoretically (15, 3234); however, to the best of our knowledge, no attempt has been made to explore the appropriateness of the PTDM approach in minimizing immortal time bias. Our primary focus in this study was to assess the performance of the PTDM for dealing with this bias. We quantified the bias due to PTDM via simulation studies as well as theoretically (see Web Appendix 1 and Web Figure 1, available at http://aje.oxfordjournals.org/) For comparison, we considered a simplified version of the sequential Cox approach (35), which was originally proposed for dealing with time-dependent confounding (36). Additionally, misclassified and excluded immortal time approaches were also compared. We applied these methods to investigate the association of time-varying β-interferon treatment with delaying disability progression in subjects from the database of a British Columbia, Canada, MS study, the Long-Term Benefits and Adverse Effects of Beta-Interferon for Multiple Sclerosis (BeAMS) Study (1995–2008) (37, 38).

METHODS

Notation

Consider a hypothetical longitudinal study consisting of n subjects (i = 1, 2, …, n). Let t0 = 0 be the start of follow-up (i.e., the time of the baseline visit). Data on baseline covariates L0 (binary) are recorded at this time. Follow-up continues until the time of failure (reaching outcome) T or censoring TC. Regular measurements of the binary treatment status Atm (Atm = 1 for treated and 0 otherwise) are recorded at intervals m = 0, 1, 2, …, K. Let [tm, tm+1) denote the mth interval. Because this study is focusing on the implications of immortal time, we assume that the subjects may initiate treatment no more than once and continue taking the treatment thereafter until the study ends (i.e., the value for treatment exposure changes only once, from 0 to 1). Let treatment initiation occur at time TA.

For group-based comparison, suppose k = 1 indicates the ever-treated group, whereas k = 0 indicates the never-treated group. Let nk indicate the size of these groups. Further, let Nk and Tk (k = 0, 1) indicate the observed number of failures and follow-up person-time in these groups, respectively. Let r = T0:T1, the ratio of the observed person-times in the never-treated and ever-treated subjects. Denote TIT as the observed immortal time—that is, the aggregated follow-up time not under treatment in the ever-treated group—and set f = TIT/T1. We describe the analytical approaches in detail in the following sections.

Cox models with a time-invariant treatment variable

Two imprecise treatment exposure definitions are often used to estimate the log hazard ratio (log HR), where a Cox model with time-invariant treatment variable A (which does not depend on the time interval) is fitted while adjusting for the potential confounders measured at the original baseline as follows:

λT(tm|L0)=λ0(tm)exp(ψ1A+ψ2L0),

where m is the visit index, λ0(tm) is the unspecified baseline hazard function, ψ1 is the log HR of the treatment status (A), and ψ2 is the vector of log HRs for the baseline covariates L0.

Included immortal time

In this approach, subjects who were ever exposed to treatment (“ever-treated”) are classified as treated for their whole duration of follow-up. Under the assumption of constant hazard, comparing the failure rate ratio from this approach (RR′) with the correct rate ratio (RR) obtained from a time-dependent analysis (see Suissa (15) and Web Appendix 2) yields

RRRR=N1/T1N0/T0N1/(T1TIT)N0/(T0+TIT)=(1f)×r(r+f).

From this equation, we can see that this approach always underestimates the correct failure rate ratio (also see Web Figure 2).

Excluded immortal time

In this approach, the immortal time—that is, time from cohort entry to the initiation of treatment—is excluded from the ever-treated subjects' follow-up, and time zero for these subjects is taken to be the time of treatment initiation, TA. However, the follow-up period for the subjects who were never exposed to treatment (“never-treated”) remains the same; that is, time zero is the time of cohort entry t0 = 0. Comparing the failure rate ratio from this approach (RR″) with the correct rate ratio obtained from a time-dependent analysis (see Suissa (15) and Web Appendix 2) yields

RRRR=N1/(T1TIT)N0/T0N1/(T1TIT)N0/(T0+TIT)=r(r+f).

From this equation, we can see that if the person-times in the ever-treated cohort are much smaller than the person-times in the never-treated cohort (i.e., if r is large), the bias from this approach may be negligible, even for large fractions of immortal time f (also see Web Figure 3).

Cox model with a time-dependent treatment variable

In the presence of the baseline covariates L0, the hazard function can be expressed as the following time-dependent Cox model:

λT(tm|L0)=λ0(tm)exp(ψ1Atm+ψ2L0),
(1)

where ψ1 is the log HR of the time-dependent treatment status (Atm).

The PTDM approach

The essence of the PTDM approach is to redefine time zero in both the ever-treated and the never-treated groups (13) by shifting the start of follow-up to the time of treatment initiation (the end of the immortal time period) TA for the ever-treated subjects. First the immortal time (wait) periods for all of the ever-treated subjects are listed. For each never-treated subject, a time zero is selected at random (with replacement) from this list and assigned as the new time zero for the corresponding never-treated subjects. If a never-treated subject has an event before the end of the assigned immortal time (new time zero), he or she is excluded from further analysis. The analysis is performed on the new time zero points—that is, the newly defined baselines after exclusion of the observed or assigned immortal time from the follow-up for the ever-treated and never-treated groups, respectively (see Web Figure 1). This is expected to eliminate the imbalance in the excluded time distribution between the two treatment groups.

Let us denote the immortal time for the jth treated subject as Tj,IT. Similarly, the matched-immortal time for the j′th untreated subject is denoted by Tj′,IT. The aggregated immortal time is denoted by TITj=1n1Tj,IT, and we denote the aggregated matched-immortal time by TITj=1n0Tj,IT. The number of failures observed during the aggregated immortal time TIT is zero. Let NIT be the number of failures observed during the aggregated matched-immortal time TIT. For NIT0, the total untreated person-time Tx of the NIT patients who had failures within their assigned immortal times is excluded. Set x = Tx/T1 ≥ 0. In addition, set TIT=q×TIT; that is, q is the ratio of assigned and observed immortal times. Here, q > 1 for the setting in which the aggregated matched-immortal time TIT is larger than the aggregated observed immortal time TIT; otherwise, 0 < q < 1. Comparing the failure rate ratio from PTDM (RR[triple prime]) with the correct rate ratio obtained from a time-dependent analysis (detailed derivation shown in Web Appendix 1) yields

RRRR=N1/(T1TIT)(N0NIT)/(T0TITTx)N1/(T1TIT)N0/(T0+TIT)=N0N0NITrq×fxr+f.

Varying r, f, x, and q yields Figure Figure1.1. This approach shows a downward bias for increasing values of f, the fraction of the immortal person-time in the ever-treated subjects. The level of bias increases with both q and x. As with the excluded immortal time approach (described in Web Appendix 2), the bias from this approach is also small if the ever-treated cohort (or its person-time contribution) is much smaller than the never-treated cohort (i.e., r is large), even for large fractions of immortal time f (see Figure Figure11).

Figure 1.
Rate ratios for the prescription time-distribution matching method (RR[triple prime]) of dealing with immortal time bias in drug effectiveness studies as compared with rate ratios for a time-dependent analysis (RR) as a function of the fraction of immortal ...

Further theoretical assessment of this approach is provided in Web Appendix 3.

The sequential Cox approach

The sequential Cox approach (35) mimics a randomized clinical trial for each of the intervals m (m = 0, 1, 2, …, K) during which subjects initiate treatment. Based on treatment initiation at the mth interval, the mth minitrial is created as follows: Only subjects who have not received any treatment before the mth interval are considered. Among these subjects, those initiating treatment during the mth interval (tm < TAtm+1) are considered the treated group, while the remaining subjects are considered the control group (for an illustration, see Web Figure 4 in Web Appendix 4). This mode of data restructuring closely resembles the “sequence of trials” approach (39) but differs with regard to some details; for example, in the current approach, treatment weights are not used to avoid potential instability in the estimated weights (35, 36) and the control subjects are artificially censored at the time of later treatment initiation (TA > tm+1) to avoid confounding due to misclassification of treatment status. To address such artificial censoring, the analysis results are adjusted using inverse-probability-of-censoring weights (36, 40).

We create the pseudodata by aggregating the data from all minitrials (see Web Appendix 4 for implementation details). Under the assumption that the different minitrials may have different baseline hazard functions but all subjects in the same minitrial will have the same baseline hazard function, we can fit a weighted stratified (time-invariant) Cox model to the pseudodata, stratified by the treatment initiation interval and weighted by the inverse-probability-of-censoring weights. Because the same subject may be included in more than 1 minitrial, we use the robust (sandwich) estimator to obtain the standard error (41, p. 170; 42, 43).

Design of simulation

To simulate survival times, we adapt the permutation algorithm (44), which has been validated for generating survival times conditional on time-dependent variables (45, 46) (see Web Appendix 5).

Simulation specifications

In all of our simulations, we generated N = 1,000 data sets, each with n = 2,000 subjects followed for up to m = 30 subsequent visits. For a given percentage of subjects, treatment initiation time TA is generated from a uniform distribution U(0, 30) (in months). In our simulations 1–3, the percentages of subjects who remain unexposed to treatment for the entire duration of follow-up (infinite TA) are set to 0%, 25%, and 75%, respectively, for the case in which there is no failure or censoring during follow-up. When subjects fail (reach the outcome) or become censored before they get a chance to initiate treatment, the percentages of never-exposed subjects increase. In all of these simulation settings, we assume an exponential distribution for generating failure times T with the constant λ0 = 0.01 rate of monthly events throughout the follow-up period. A uniform distribution U(1, 60) months is assumed to generate censoring times TC. Administrative censoring is set at 30 months of follow-up for all subjects who remain at risk at that time. For each simulated data set, we can calculate the exact percentages of ever-treated and never-treated subjects. For any subject, treatment status is binary for any given month, either exposed or unexposed, where a person is considered exposed throughout the month of the treatment initiation. To focus on the immortal time issue, we assume there are no discontinuations or interruptions for persons who initiate treatment. Additionally, we consider sex as a baseline confounder. A subject's sex is generated based on a Bernoulli distribution where the probability of being male is 0.3 (as might occur in an MS population).

After generating values for the survival time Ti, the censoring time TiC, and the treatment and covariate matrix Xim = (Aim, Li0) for each subject i = 1, 2, …, n for up to m = 30 time periods, the permutation algorithm (44) is used to generate survival data, where treatment Atm is time-dependent but the confounder L0 is fixed at the baseline value. Arbitrarily, the parameters for the effects of treatment and sex on the survival outcome are set such that the treatment is associated with the outcome in a harmful direction (a log HR of ψ1 = 0.5) and males are at a lower risk than females (a log HR of ψ2 = −0.7).

In order to consider more frequent events, we repeat the Monte Carlo study (similar to simulation 1) with λ0 = 0.10 (on a monthly scale) (simulation 4). We further assume a gamma distribution (scale = 1/0.01, shape = 0.4) and a Weibull distribution (scale = 1/0.01, shape = 2) for generating failure times T (simulations 5 and 6, respectively; all other settings are similar to simulation 1).

Performance metrics

We assess the performance of the various approaches by means of the following measures:

  • Bias(ψ1ˆ)=i=1N(ψ1iˆψ1)/N: the average difference between the true parameters and N = 1,000 estimated parameters (log HR).
  • Standarddeviation(ψ1ˆ)=i=1N(ψ1iˆψ1¯)2/(N1), where ψ1¯=i=1Nψˆ1iN.
  • Model-based (average) standard error: The average of N = 1,000 estimated standard errors of the estimated causal association parameter.
  • Coverage probabilities of model-based nominal 95% confidence intervals: proportion of N = 1,000 data sets in which the true parameter is contained in the nominal 95% confidence interval.
  • Meansquarederror(ψ1ˆ)=i=1N(ψ1iˆψ1)2/N.

SIMULATION RESULTS

Description of the simulated data

To describe the data sets generated from each of the simulations, we generated 1 data set for each simulation with a larger number of subjects (n = 25,000). Figure Figure22 shows the cumulative percentages of the cohort having initiated treatment over the 30 follow-up periods.

Figure 2.
Cumulative percentages of a simulated cohort of 25,000 multiple sclerosis patients having initiated treatment over 30 follow-up periods. A) Simulation 1, with percentages of “failure,” “never-treated,” and “ever-treated” ...

Rare event condition

The results from the simulations of the rare event condition (λ0 = 0.01 on a monthly time scale), simulations 1–3, are given in Tables Tables113, respectively.

Table 1.
Results Derived From Analytical Approaches to Adjusting for Immortal Time Bias in a Simulation of 1,000 Data Sets, Each Containing 2,000 Multiple Sclerosis Patients Followed for Up to 30 Time Intervals (Rare Event Case λ0 = 0.01) (Simulation 1) ...
Table 3.
Results Derived From Analytical Approaches to Adjusting for Immortal Time Bias in a Simulation of 1,000 Data Sets, Each Containing 2,000 Multiple Sclerosis Patients Followed for Up to 30 Time Intervals (Rare Event Case λ0 = 0.01) (Simulation 3) ...

When the immortal time is misclassified as exposed time, we see a substantial downward bias. The bias improves slightly when excluded immortal time and PTDM approaches are used, but these estimate are still off the target. As the ratio for the subjects in the never-treated group versus the ever-treated group increases, the bias monotonically goes down (see Tables Tables22 and and3).3). The variability of the estimator is larger for PTDM than for the time-dependent Cox results. Interestingly, the standard errors of the estimates from both the misclassified and excluded immortal-time approaches are smaller than those for the time-dependent Cox results, but the mean squared errors are generally higher. When the sequential Cox approach is applied, in all 3 of these simulations, the biases are negligible and the average coverage probabilities of the model-based nominal 95% confidence intervals are close to 0.95; both are comparable to the time-dependent Cox results. This estimator is generally more variable than that from the time-dependent Cox model fit (from simulations 1 and 2). However, from simulation 3 (Table (Table3),3), we see that exceptions may be possible.

Table 2.
Results Derived From Analytical Approaches to Adjusting for Immortal Time Bias in a Simulation of 1,000 Data Sets, Each Containing 2,000 Multiple Sclerosis Patients Followed for Up to 30 Time Intervals (Rare Event Case λ0 = 0.01) (Simulation 2) ...

When more events are available

Results from the more frequent event condition (λ0 = 0.10 on a monthly time scale) are presented in Table Table4.4. The bias for misclassified, excluded immortal-time, and PTDM approaches exhibiting large bias in Table Table11 is considerably reduced but is still large, whereas the time-dependent Cox and sequential Cox approaches again exhibit negligible bias. For the PTDM approach, less bias is expected when the event rate is higher, as is discussed via an example of the analytical results in Web Appendix 3. In addition, the standard errors are much larger when failure rates are low.

Table 4.
Results Derived From Analytical Approaches to Adjusting for Immortal Time Bias in a Simulation of 1,000 Data Sets, Each Containing 2,000 Multiple Sclerosis Patients Followed for Up to 30 Time Intervals (More Frequent Event λ0 = 0.10) (Simulation ...

When event times are generated from gamma and Weibull distributions

When event times are generated from a gamma distribution and a Weibull distribution, the results (Web Tables 1 and 2) closely resemble those obtained when event times are generated from an exponential distribution (Table (Table11).

Application to MS

We applied the methods described in this paper to data on the British Columbia MS cohort used earlier in the BeAMS Study (1995–2008) (38, 4751). Web Appendix 6 shows the baseline characteristics of this MS cohort (see Web Table 3) and gives the study eligibility and exclusion criteria. In this retrospective observational study, the outcome measure was time from β-interferon eligibility to a “confirmed and sustained” Expanded Disability Status Scale score of 6, considered irreversible disability when all subsequent such scores were at least 6, with 1 or more measurements being taken at least 150 days later (52, 53). The goal was to assess the association of β-interferon with irreversible disease progression in MS patients.

Results are reported in Table Table5.5. None of the approaches showed a significant hazard ratio. The hazard ratio estimated from the included immortal time approach was the lowest, and this was expected from the simulation. The hazard ratio estimated from the sequential Cox approach was slightly lower than that estimated from the time-dependent Cox analysis, while that estimated from the PTDM was even lower. The distributions of excluded times from both groups when using the PTDM method looked reasonably similar (Figure (Figure3).3). Because the PTDM approach produces different estimates from the same data due to the random sampling of the control subjects, we estimated the hazard ratio from the British Columbia MS data 1,000 times (average hazard ratio = 1.26; listed in Table Table5).5). The distribution of the estimated hazard ratios was moderately symmetrical (see Web Figure 5). Density plots of the estimated inverse-probability-of-censoring weights are shown in Web Figure 6.

Table 5.
Estimated Parameters (Using the Hazard Scale) for Multiple Sclerosis Patients From British Columbia, Canada, 1995–2008
Figure 3.
Matched wait periods (or matched immortal times), in years, from the prescription time-distribution matching approach in a simulated cohort of 25,000 multiple sclerosis patients initiating β-interferon treatment, British Columbia, Canada, 1995–2008. ...

To better understand the possible level of bias from each approach in this context, we conducted another Monte Carlo study, choosing data simulation settings similar to those of the MS data set. Simulation specifications are shown in Web Appendix 7 (simulated cohort characteristics are shown in Web Table 4 and Web Figure 7). The results are reported in Table Table6.6. Here, the level of bias from each approach is much less than that reported for the previous simulations (Table (Table66 vs. Tables Tables11 and and2)2) but similar to that in simulation 3 (Table (Table33).

Table 6.
Comparison of Analytical Approaches to Adjusting for Immortal Time Bias in a Simulation Inspired by Multiple Sclerosis Data From 1,000 Data Sets, Each Containing 1,700 Subjects Followed for Up to 150 Monthsa

DISCUSSION

In observational studies of drug effectiveness, subjects may not be exposed to the treatment at the beginning of follow-up. Statistical procedures, such as the time-dependent Cox model, are known to deal with time-dependent treatment. However, in the medical literature, it is common to see the Cox model applied on the basis of baseline treatment information, likely for the convenience of model-fitting and group-based interpretation (54). Therefore, there is a need for methods that are capable of handling the time-dependent nature of longitudinal data and will also help us to better understand the treatment-outcome mechanism. To this end, the PTDM approach has been proposed for the situation where treatment initiation occurs later than cohort entry (13). The sequential Cox approach has been proposed for more complex longitudinal settings (35), but we considered a simplified version of this approach to deal with the immortal time bias.

In the PTDM approach, treatment exposure is converted into a time-independent variable so that a simple Cox model for treatment-group comparison can be applied. New time zeros are defined after excluding the observed and assigned immortal times in the ever- and never-treated groups. The excluded times (wait periods) for the ever-treated and never-treated groups follow the same distribution. However, the baseline is not the same for all subjects. The question of whether this approach (assigning the baseline borrowed from the ever-treated cohort) adequately addresses the immortal time bias has not previously been investigated in a systematic fashion, to our knowledge. From our bias-quantifying equations (Web Appendices 1 and 3) and the results of our simulations, we can see that the bias is still substantial in the PTDM analysis and the direction of the bias is negative, highlighting the value of setting a well-defined time zero or baseline. The results of simulations 1–3 (also in Web Appendix 1) identify an increase in the ratio of subjects in the never-treated group to subjects in the ever-treated group as an important design parameter for reducing bias associated with the PTDM method (biases were −1.411, −0.558, and −0.101 for the ratios 0.49, 0.98, and 4.95, respectively). From simulation 4, we learned that less bias is expected for the PTDM approach when the event rate is higher. We also considered 2 other popular event-time-generating distributions, namely the gamma and Weibull distributions (simulations 5 and 6, respectively), and the results were consistent with those from simulation 1.

In contrast, the sequential Cox approach allows a consistent definition of baseline and utilizes all subjects. The baseline is set on the basis of treatment initiation, and the subjects are stratified accordingly. This approach recreates the covariate process at each treatment start using the minitrial approach (36) and works very well in our simulation studies in terms of bias. In practice, as the time-dependent Cox model approach can easily deal with the immortal time bias and is simpler to implement, one may not be inclined to use the sequential Cox approach, because additional steps are involved in the data-restructuring for the latter approach (i.e., to create minitrials and then pseudodata). The data associated with a given minitrial in a sequential Cox approach, however, can be extracted and separated quite easily from the pseudodata, and it is straightforward to compare the treatment associations (of early versus late treatment initiation) with the outcome. It is also possible to estimate the treatment association for patients with a specific level of a covariate at treatment initiation (35). Such detailed scrutiny using this minitrial approach could yield insights about the data which may be hard to gain using the time-dependent Cox approach.

We applied the methods under consideration to estimate the association of β-interferon with disease progression in the British Columbia MS cohort. In comparison with results from the time-dependent Cox model (hazard ratio (HR) = 1.29, 95% confidence interval (CI): 0.91, 1.81), the sequential Cox approach provided a slightly lower association estimate with a wider confidence interval (HR = 1.22, 95% CI: 0.72, 2.08). Repeating the PTDM analysis 1,000 times, the average results (HR = 1.26, average 95% CI: 0.86, 1.85) were quite similar to those from the time-dependent Cox approach, although the average confidence interval for the estimated treatment association was slightly wider. We designed another Monte Carlo study using simulation settings reflecting our MS data set to explore further the level of bias from each approach. The smaller amount of bias in this simulation for each approach may have been a consequence of the fact that the ever-treated group was much smaller than the never-treated group in this simulation (19% vs. 81%; ratio: 4.26). This phenomenon also explains why the PTDM approach provided estimates close to those from the time-dependent Cox model fit (29% user vs. 71% nonuser) (13). Even in this simulation setting, the time-dependent Cox model and sequential Cox approach were still associated with less bias than the PTDM approach.

Similar to other simulation studies, we investigated only a few scenarios. However, the assumptions underlying the data simulation were consistent with patterns typical of epidemiologic observational survival studies, where treatment initiation may happen later for some subjects than others (7). Furthermore, our assumption of no discontinuations or interruptions in the treatment was restrictive and may not be suitable for more complex disease scenarios in which an investigator might wish to assess different treatment strategies (i.e., switching between therapies) over the course of time. Substantial immortal time bias was introduced by group-based Cox model analyses in the scenarios investigated (analytical expressions for PTDM approach in Web Appendix 3 also support this finding). However, even in these scenarios, the sequential Cox method appeared to estimate the target parameter adequately. Variance estimation is still a challenge with the sequential Cox model and similar methods (55, 56), and the resampling methods suggested for estimating the variance are computationally demanding (35, 57). Future research could focus on assessing the ability of the sequential Cox approach to adjust for time-dependent confounders (37).

Supplementary Material

Web Material:

ACKNOWLEDGMENTS

Author affiliations: Department of Epidemiology, Biostatistics and Occupational Health, Faculty of Medicine, McGill University (Mohammad Ehsanul Karim); Department of Statistics, Faculty of Science, University of British Columbia, Vancouver, British Columbia, Canada (Mohammad Ehsanul Karim, Paul Gustafson, John Petkau); Division of Neurology, Department of Medicine, Faculty of Medicine, University of British Columbia, Vancouver, British Columbia, Canada (Helen Tremlett); and Brain Research Centre, University of British Columbia, Vancouver, British Columbia, Canada (Helen Tremlett).

Members of the BeAMS Study Group: Afsaneh Shirani, Yinshan Zhao, Charity Evans, Elaine Kingwell, Mia L. van der Kop, and Joel Oger.

This work was supported by a studentship from the MS [Multiple Sclerosis] Society of Canada (MSSOC) (M.E.K.) and by grants from the National Multiple Sclerosis Society (NMSS) (grant RG 4202-A-2; Principal Investigator (PI): H.T.) and the Canadian Institutes of Health Research (CIHR) (grant MOP-93646; PI: H.T.). J.P. holds research grants from the NMSS and the Natural Sciences and Engineering Research Council of Canada. P.G. is supported by the Natural Sciences and Engineering Research Council of Canada. H.T. is supported by the Canada Research Chair Program and a MSSOC Don Paty Career Development Award and is a Michael Smith Foundation for Health Research scholar. She has also received research support from the NMSS, the CIHR, and the Multiple Sclerosis Trust (United Kingdom). Y.Z. receives research funding from the CIHR, the MSSOC, and the NMSS. The work of A.S. was funded through a postdoctoral fellowship from the MSSOC and grants from the CIHR (grant MOP-93646; PI: H.T.) and the NMSS (grant RG 4202-A-2; PI: H.T.). E.K. was supported by postdoctoral fellowships from the Michael Smith Foundation for Health Research and the MSSOC. The work of C.E. was funded through grants from the CIHR (grant MOP-93646; PI: H.T.) and the NMSS (grant RG 4202-A-2; PI: H.T.) and the Michael Smith Foundation for Health Research. M.L.v.d.K. is supported by a CIHR Doctoral Award-Doctoral Foreign Study Award (October 2012), offered in partnership with the CIHR Strategy for Patient-Oriented Research and the CIHR HIV/AIDS Research Initiative. J.O. receives support from the Christopher Foundation and the University of British Columbia.

We gratefully acknowledge the neurologists of the University of British Columbia (UBC) Hospital MS Clinic and Research Group (current members listed here by primary clinic), who contributed to the study through patient examination and data collection: UBC MS Clinic—Drs. A. Traboulsee (clinic director and head of UBC MS programs), A.-L. Sayao, V. Devonshire, S. Hashimoto (UBC and Victoria MS Clinics), J. Hooge, L. Kastrukoff (UBC and Prince George MS Clinics), and J. Oger; Kelowna MS Clinic—Drs. D. Adams, D. Craig, and S. Meckling; Prince George MS Clinic—Dr. L. Daly; Victoria MS Clinic—Drs. O. Hrebicek, D. Parton, and K. Pope.

The views expressed in this paper do not necessarily reflect the views of each individual acknowledged. None of the authors received compensation for their role in the study.

M.E.K. has had travel and accommodation costs covered by the endMS Research and Training Network (ERTN) to make presentations at conferences (2011, 2012) and has received financial support from the Pacific Institute for the Mathematical Sciences to attend a workshop (2013). P.G. has received consulting fees from Biogen to serve on an advisory board. Over the past 3 years, J.P. has received consulting fees and/or fees for service on Data Safety Monitoring Boards from EMD Serono, the Myelin Research Foundation, and Novartis. H.T. has received speaker honoraria and/or money for travel expenses to attend conferences of the Consortium of MS Centres (2013), the MS Society of Canada (2013), the National Multiple Sclerosis Society (2012, 2014), Bayer Pharmaceutical (speaker, 2010; honoraria declined), Teva Pharmaceutical (speaker, 2011), the European Committee for Treatment and Research in Multiple Sclerosis (ECTRIMS) (2011, 2012, 2013, 2014), the Multiple Sclerosis Trust (2011), the Chesapeake Health Education Program, the US Department of Veterans Affairs (2012, honorarium declined), Novartis Canada (2012), Biogen Idec (2014; honorarium declined), and the American Academy of Neurology (2013, 2014, 2015). Unless otherwise stated, H.T. donates all speaker honoraria to either an MS charity or an unrestricted grant fund for use by her research group. A.S. has received travel grants to attend and make presentations at conferences of the ERTN (2010, 2011), ECTRIMS (2010, 2011), and the Consortium of MS Centres (2012). E.K. has had travel and accommodation costs covered to attend and make presentations at conferences of the ERTN (2008, 2011), the International Society for Pharmacoepidemiology (2010, 2013), and Bayer Schering Pharma (2010). C.E. has received travel grants to attend and make presentations at conferences of the ERTN (2011) and ECTRIMS (2011). Over the past 5 years, J.O. has received speaker honoraria, consulting fees, travel grants, research grants, and/or educational grants from Aventis, Bayer, Biogen Idec, BioMS, Corixa, Genentech, GlaxoSmithKlein, Merck-Serono, Novartis, Serono, Shering, Talecris, and Teva Neuroscience. J.O. has received fees for services from Bayer, Novartis, and Biogen Idec to serve on advisory committees.

REFERENCES

1. Trojano M, Pellegrini F, Fuiani A et al. New natural history of interferon-beta-treated relapsing multiple sclerosis. Ann Neurol. 2007;614:300–306. [PubMed]
2. Lévesque LE, Hanley JA, Kezouh A et al. Problem of immortal time bias in cohort studies: example using statins for preventing progression of diabetes. BMJ. 2010;340:b5087. [PubMed]
3. Renoux C, Suissa S Immortal time bias in the study of effectiveness of interferon-beta in multiple sclerosis [letter] Ann Neurol. 2008;641:109–110. [PubMed]
4. Wolkewitz M, Allignol A, Harbarth S et al. Time-dependent study entries and exposures in cohort studies can easily be sources of different and avoidable types of bias. J Clin Epidemiol. 2012;6511:1171–1180. [PubMed]
5. Gail MH. Does cardiac transplantation prolong life? A reassessment. Ann Intern Med. 1972;765:815–817. [PubMed]
6. van Walraven C, Davis D, Forster AJ et al. Time-dependent bias was common in survival analyses published in leading clinical journals. J Clin Epidemiol. 2004;577:672–682. [PubMed]
7. Suissa S. Immortal time bias in observational studies of drug effects. Pharmacoepidemiol Drug Saf. 2007;163:241–249. [PubMed]
8. Kiri VA, Pride NB, Soriano JB et al. Inhaled corticosteroids in chronic obstructive pulmonary disease: results from two observational designs free of immortal time bias. Am J Respir Crit Care Med. 2005;1724:460–464. [PubMed]
9. Austin PC, Platt RW Survivor treatment bias, treatment selection bias, and propensity scores in observational research. J Clin Epidemiol. 2010;632:136–138. [PubMed]
10. Ravi B, Croxford R, Austin PC et al. The relation between total joint arthroplasty and risk for serious cardiovascular events in patients with moderate-severe osteoarthritis: propensity score matched landmark analysis. BMJ. 2013;347:f6187. [PubMed]
11. Karim ME. Can joint replacement reduce cardiovascular risk? [editorial] BMJ. 2013;347:f6651. [PubMed]
12. Suissa S. Effectiveness of inhaled corticosteroids in chronic obstructive pulmonary disease: immortal time bias in observational studies. Am J Respir Crit Care Med. 2003;1681:49–53. [PubMed]
13. Zhou Z, Rahme E, Abrahamowicz M et al. Survival bias associated with time-to-treatment initiation in drug effectiveness evaluation: a comparison of methods. Am J Epidemiol. 2005;16210:1016–1023. [PubMed]
14. Austin PC, Mamdani MM, Van Walraven C et al. Quantifying the impact of survivor treatment bias in observational studies. J Eval Clin Pract. 2006;126:601–612. [PubMed]
15. Suissa S. Immortal time bias in pharmacoepidemiology. Am J Epidemiol. 2008;1674:492–499. [PubMed]
16. Wolkewitz M, Allignol A, Schumacher M et al. Two pitfalls in survival analyses of time-dependent exposure: a case study in a cohort of Oscar nominees. Am Stat. 2010;643:205–211.
17. Sylvestre MP, Huszti E, Hanley JA Do OSCAR winners live longer than less successful peers? A reanalysis of the evidence. Ann Intern Med. 2006;1455:361–363. [PubMed]
18. Ho PM, Fihn SD, Wang L et al. Clopidogrel and long-term outcomes after stent implantation for acute coronary syndrome. Am Heart J. 2007;1545:846–851. [PubMed]
19. Karp I, Behlouli H, Lelorier J et al. Statins and cancer risk. Am J Med. 2008;1214:302–309. [PubMed]
20. Ho PM, Maddox TM, Wang L et al. Risk of adverse outcomes associated with concomitant use of clopidogrel and proton pump inhibitors following acute coronary syndrome. JAMA. 2009;3019:937–944. [PubMed]
21. Snyder CW, Weinberg JA, McGwin G Jr et al. The relationship of blood product ratio to mortality: survival benefit or survival bias? J Trauma. 2009;662:358–362. [PubMed]
22. Keyhan G, Chen SF, Pilote L The effectiveness of beta-blockers in women with congestive heart failure. J Gen Intern Med. 2007;227:955–961. [PMC free article] [PubMed]
23. Keyhan G, Chen SF, Pilote L Angiotensin-converting enzyme inhibitors and survival in women and men with heart failure . Eur J Heart Fail. 2007;9(6-7):594–601. [PubMed]
24. Morin S, Rahme E, Behlouli H et al. Effectiveness of antiresorptive agents in the prevention of recurrent hip fractures. Osteoporos Int. 2007;1812:1625–1632. [PubMed]
25. Walkey AJ, Wiener RS Utilization patterns and patient outcomes associated with use of rescue therapies in acute lung injury. Crit Care Med. 2011;396:1322–1328. [PMC free article] [PubMed]
26. Manari A, Varani E, Guastaroba P et al. Long-term outcome in patients with ST segment elevation myocardial infarction and multivessel disease treated with culprit-only, immediate, or staged multivessel percutaneous revascularization strategies: insights from the REAL registry. Catheter Cardiovasc Interv. 2014;846:912–922. [PubMed]
27. Prieto-Alhambra D, Javaid MK, Judge A et al. Hormone replacement therapy and mid-term implant survival following knee or hip arthroplasty for osteoarthritis: a population-based cohort study. Ann Rheum Dis. 2015;743:557–563. [PubMed]
28. Li G, Holbrook A, Delate T et al. Prediction of individual combined benefit and harm for patients with atrial fibrillation considering warfarin therapy: a study protocol. BMJ Open. 2015;511:e009518. [PMC free article] [PubMed]
29. Shintani AK, Girard TD, Eden SK et al. Immortal time bias in critical care research: application of time-varying cox regression for observational cohort studies. Crit Care Med. 2009;3711:2939–2945. [PMC free article] [PubMed]
30. Liu J, Weinhandl ED, Gilbertson DT et al. Issues regarding ‘immortal time’ in the analysis of the treatment effects in observational studies. Kidney Int. 2012;814:341–350. [PubMed]
31. Ho AM, Dion PW, Yeung JH et al. Simulation of survivorship bias in observational studies on plasma to red blood cell ratios in massive transfusion for trauma. Br J Surg. 2012;99(suppl 1):132–139. [PubMed]
32. Buyse M, Piedbois P On the relationship between response to treatment and survival time. Stat Med. 1996;1524:2797–2812. [PubMed]
33. Beyersmann J, Gastmeier P, Wolkewitz M et al. An easy mathematical proof showed that time-dependent bias inevitably leads to biased effect estimation. J Clin Epidemiol. 2008;6112:1216–1221. [PubMed]
34. Beyersmann J, Wolkewitz M, Schumacher M The impact of time-dependent bias in proportional hazards modelling. Stat Med. 2008;2730:6439–6454. [PubMed]
35. Gran JM, Røysland K, Wolbers M et al. A sequential Cox approach for estimating the causal effect of treatment in the presence of time-dependent confounding applied to data from the Swiss HIV Cohort Study. Stat Med. 2010;2926:2757–2768. [PubMed]
36. Aalen O. Armitage Lecture 2010. Understanding treatment effects: the value of integrating longitudinal data and survival analysis. Stat Med. 2012;3118:1903–1917. [PubMed]
37. Hernán MA, Brumback B, Robins JM Marginal structural models to estimate the causal effect of zidovudine on the survival of HIV-positive men. Epidemiology. 2000;115:561–570. [PubMed]
38. Shirani A, Zhao Y, Karim ME et al. Association between use of interferon beta and progression of disability in patients with relapsing-remitting multiple sclerosis. JAMA. 2012;3083:247–256. [PubMed]
39. Hernán MA, Alonso A, Logan R et al. Observational studies analyzed like randomized experiments: an application to postmenopausal hormone therapy and coronary heart disease. Epidemiology. 2008;196:766–779. [PMC free article] [PubMed]
40. Robins JM, Hernán MA, Brumback B Marginal structural models and causal inference in epidemiology. Epidemiology. 2000;115:550–560. [PubMed]
41. Therneau T. Modeling Survival Data: Extending the Cox Model. New York, NY: Springer-Verlag New York; 2000.
42. Xiao Y, Abrahamowicz M, Moodie EE Accuracy of conventional and marginal structural Cox model estimators: a simulation study. Int J Biostat. 2010;62:Article 13. [PubMed]
43. Taylor JM, Shen J, Kennedy EH et al. Comparison of methods for estimating the effect of salvage therapy in prostate cancer when treatment is given by indication. Stat Med. 2014;332:257–274. [PMC free article] [PubMed]
44. Abrahamowicz M, Mackenzie T, Esdaile J Time-dependent hazard ratio: modeling and hypothesis testing with application in lupus nephritis. J Am Stat Assoc. 1996;91436: 1432–1439.
45. Mackenzie T, Abrahamowicz M Marginal and hazard ratio specific random data generation: applications to semi-parametric bootstrapping. Stat Comput. 2002;123:245–252.
46. Sylvestre MP, Abrahamowicz M Comparison of algorithms to generate event times conditional on time-dependent covariates. Stat Med. 2008;2714:2618–2634. [PubMed]
47. Karim ME, Gustafson P, Petkau J et al. Marginal structural Cox models for estimating the association between β-interferon exposure and disease progression in a multiple sclerosis cohort. Am J Epidemiol. 2014;1802:160–171. [PMC free article] [PubMed]
48. Karim M. Causal Inference Approaches for Dealing With Time-Dependent Confounding in Longitudinal Studies, With Applications to Multiple Sclerosis Research [doctoral thesis] Vancouver, BC, Canada: University of British Columbia; 2015.
49. Shirani A, Zhao Y, Karim ME et al. Investigation of heterogeneity in the association between interferon beta and disability progression in multiple sclerosis: an observational study. Eur J Neurol. 2014;216:835–844. [PubMed]
50. Zhang T, Shirani A, Zhao Y et al. Beta-interferon exposure and onset of secondary progressive multiple sclerosis. Eur J Neurol. 2015;226:990–1000. [PubMed]
51. Shirani A, Zhao Y, Petkau J et al. Multiple sclerosis in older adults: the clinical profile and impact of interferon beta treatment. Biomed Res Int. 2015;2015:451912. [PMC free article] [PubMed]
52. Tremlett H, Paty D, Devonshire V Disability progression in multiple sclerosis is slower than previously reported. Neurology. 2006;662:172–177. [PubMed]
53. Tremlett H, Yousefi M, Devonshire V et al. Impact of multiple sclerosis relapses on progression diminishes with time. Neurology. 2009;7320:1616–1623. [PMC free article] [PubMed]
54. Leffondré K, Abrahamowicz M, Siemiatycki J Evaluation of Cox's model and logistic regression for matched case-control data with time-dependent covariates: a simulation study. Stat Med. 2003;2224:3781–3794. [PubMed]
55. Schaubel DE, Wolfe RA, Port FK A sequential stratification method for estimating the effect of a time-dependent experimental treatment in observational studies. Biometrics. 2006;623:910–917. [PubMed]
56. Schaubel D, Wolfe R, Sima C et al. Estimating the effect of a time-dependent treatment by levels of an internal time-dependent covariate: application to the contrast between liver wait-list and posttransplant mortality. J Am Stat Assoc. 2009;104485:49–59.
57. Røysland K, Gran JM, Ledergerber B et al. Analyzing direct and indirect effects of treatment using dynamic path analysis applied to data from the Swiss HIV Cohort Study. Stat Med. 2011;3024:2947–2958. [PubMed]
58. Aalen OO. A linear regression model for the analysis of life times. Stat Med. 1989;88:907–925. [PubMed]
59. Aalen OO. Further results on the non-parametric linear model in survival analysis. Stat Med. 1993;1217:1569–1588. [PubMed]

Articles from American Journal of Epidemiology are provided here courtesy of Oxford University Press