Comparative performance of a modified landmark approach when no time of treatment data are available within oncological databases: exemplary cohort study among resected pancreatic cancer patients.

PURPOSE: The Mantel-Byar method is the gold standard analytical approach to avoid immortal time bias, but requires information on the time between start of follow-up and exposure initiation. Alternatively, a modified landmark approach might be used to mitigate the amount of immortal time bias, which assumes exposure initiation at a predefined landmark time. In the context of an expected positive association between adjuvant chemotherapy (ACT) and overall survival among resected pancreatic cancer (PCa) patients, this study aims to empirically assess the performance of this approach relative to the Mantel-Byar method. PATIENTS AND METHODS: Data from resected PCa patients diagnosed between 2003 and 2014 and registered in the national cancer registries of Belgium, the Netherlands, and Slovenia were used to estimate the association between ACT and overall survival using a Cox proportional hazards model by country and overall. Results derived from the immortal time bias (misclassifying the time to ACT initiation), Mantel-Byar method, and conventional and modified landmark analyses with assumed cutoff times of ACT initiation at 9, 12 and 15 weeks post-diagnosis were compared.
RESULTS: In total, 5,668 patients with a total of 10,921 person-years of follow-up were eligible. All analytical approaches showed a significant survival benefit for resected PCa patients who received ACT, but immortal time bias analyses led to strong overestimation of ACT benefits compared to the Mantel-Byar method (immortal time bias: overall HR [95% CI] 0.68 [0.62-0.75] vs Mantel-Byar method: 0.82 [0.71-0.93]), whereas the conventional landmark approach generally provided rather conservative estimates (0.86 [0.75-1.00], 15 weeks landmark). HRs derived from modified landmark analyses depended on the cutoff time, but were similar compared to the Mantel-Byar method at 15 weeks (0.81 [0.70-0.94]).
CONCLUSION: A modified landmark approach might be a valid alternative to the Mantel-Byar method if no time of treatment information is available. The performance depends on the chosen cutoff time.

Introduction

Associations between treatment interventions and outcomes such as mortality or recurrence are of major interest in observational cancer research. Therefore, studies are often conceptualized and analyzed as emulated target trials, which holds the potential of serious biases when the start of follow-up (time zero) and the time of treatment initiation are not aligned.1 This often leads to immortal time bias – a time-dependent bias still frequently seen in medical literature.2 In brief, immortal time bias refers to a cohort study design scenario where the exposure of a cohort starts after the start of follow-up, so that patients who received treatment are “immortal” up to that time, which results in a biased mortality rate ratio between the exposed and unexposed in favor of the exposure. The bias is introduced when the unexposed person-time between start of follow-up and exposure initiation is not modeled as such.3,4

However, immortal time bias can be avoided if appropriate statistical models are used. Recent studies confirmed that the Mantel–Byar method – a time-dependent approach to model each patient’s person-time under exposure – is the gold standard when analyzing time-varying interventions.5–7 However, depending on the available data source, the time of treatment information is often not registered, which precludes a time-dependent approach. In this case, a landmark analysis may be considered as an alternative analytical strategy.

The conventional landmark approach, introduced by Anderson et al,8 stratifies patients into two cohorts according to whether they have initiated the treatment before a predefined landmark time and follows them from this cutoff time whereas patients who died before are excluded. However, this still requires the temporal information about the exposure status of each patient at the landmark time. An alternative approach –referred to, in the following, as a modified landmark approach – might still give the opportunity to assess treatment associations without requiring any time of treatment data. Using this approach, clinical criteria (e.g., guidelines) can be exploited to make assumptions on possible landmark cutoff times and the patient’s probability to have initiated the treatment at this time.

As only little is known about the validity of the derived results relative to the Mantel–Byar method as the gold standard, this study aims to assess and illustrate the comparative performance of a modified landmark approach using an exemplary population-based study with the Mantel–Byar approach as gold standard to which all other analytical approaches will be compared to. Therefore, the different analytical approaches are applied within a population-based setting to assess the well-known positive clinical benefit of adjuvant chemotherapy (ACT) for resected pancreatic cancer (PCa) patients, for which randomized controlled trials (RCTs) have found a significant reduction in mortality.9

Patients and methods

Data sources

For this methodological study, data from three different European population-based cancer registries (Belgium, the Netherlands, and Slovenia) were used.

The registries provided data for the years of diagnosis between 2003 and 2014 with end of follow-up between February 2015 and May 2016. According to a predefined protocol, all registries collected information on treatment, date of diagnosis, start date of ACT, and further possible confounder variables and prognostic factors including age at diagnosis, sex, calendar year, tumor location, morphology, staging, previous other tumors, and harvested lymph nodes. All of the important variables were categorized and harmonized for statistical analysis, and a detailed overview of these variables is provided in Table S1.

Ethical approval for the analysis of the study data was granted by the Medical Faculty of the University of Heidelberg (code S-319/2015).

Study population

Eligible patients were included in the study if they had a diagnosis of PCa (ICD-10 code C25.x) throughout any of the available calendar years and underwent surgery. The date of diagnosis was defined as either the date of pathological or clinical diagnosis (whichever was available). Patients with diagnosis of a neuroendocrine tumor (C25.4) or who survived less than a week after the diagnosis were excluded. Patients with missing or implausible values on important variables were excluded (8.9%; details are shown in Table S1).

Clinical research question and study design

Major RCTs have shown that the only potentially curative treatment approach in PCa is the employment of surgical tumor resection followed by ACT with either 5-fluorouracil (5-FU) or gemcitabine.9–15 To emulate these trials, a retrospective cohort design was used within all analytical approaches, with ACT as the exposure variable of interest. Because the ESPAC-3 v2 trial showed that with regard to the efficacy of ACT there is no sufficient difference between 5-FU and gemcitabine, exposure was not restricted to any specific antineoplastic compound.16 The outcome was defined as time to death due to any reason.

Analytical approaches

The different analytical approaches to model exposure are described in detail in the following sections and are graphically summarized in Figure 1.

Figure 1

Illustration of study design choices including the immortal time-biased model (A), the Mantel–Byar model (B), the conventional landmark model (C), and the modified landmark model (D).

Analysis model 1: immortal time-biased approach

In order to assess the magnitude of how much an immortal time bias might have affected the results of the above-described study by neither considering a time-dependent nor a landmark analysis, naive estimates were calculated using a time-fixed approach by defining the chemotherapy exposure as a dichotomous variable (chemotherapy vs no chemotherapy; Figure 1A). This artificially leads to a misclassification of the exposed person-time of those patients who received chemotherapy as they are considered “ever exposed” (Figure 1A continuous green line patients 1–3), whereas those who did not receive chemotherapy are considered “never exposed” (Figure 1A, continuous red line Patient 4). This leads to a biased rate ratio in favor of the chemotherapy, as those patients who received chemotherapy must have survived the time between surgery and chemotherapy initiation whereas those who did not might have died immediately after surgery.

Analysis model 2: Mantel–Byar approach

Using a time-dependent approach according to Mantel–Byar,5 every patient contributed unexposed and exposed person-time according to the patient’s exact exposure status at every time point during follow-up. Therefore, the whole follow-up time for every patient was divided into 1-week intervals, indicating the patient’s exposure status (chemotherapy vs no chemotherapy) in each of these intervals. This means that a patient who received chemotherapy, after a certain time after surgery, first contributed unexposed person-time from the start of follow-up until the interval when ACT was initiated, and was then modeled as exposed thereafter until the end of follow-up (Figure 1B).

Analysis models 3 and 4: conventional and modified landmark approaches with cutoff at 9, 12, and 15 weeks

If time of treatment information is available, conventional landmark analyses might be conducted alternatively to the Mantel–Byar approach, although precision might be lost and the derived results might be more prone to be biased toward the null due to exposure misclassification of patients who started the exposure after the landmark time but are falsely modeled as unexposed throughout the follow-up. In this study, analyses using a conventional landmark approach were undertaken by assigning each patient’s exposure status conditional on whether the patient had received the exposure before the respective landmark time (Figure 1C).8 The landmark cutoff time was then set as the new start of follow-up and patients who died before were excluded.

However, as time of treatment initiation is often not recorded within cancer registries, a modified landmark analysis was carried out defining landmark cutoff times at which it can be assumed that the majority of patients have initiated ACT. Thus, patients were modeled as exposed or unexposed starting from this landmark time based on the information as to whether they actually ever received ACT (yes/no), whereas patients who died before were excluded from the analysis. Contrary to the conventional landmark approach, patients who had started ACT after the landmark time were still modeled as exposed (Figure 1D).

The landmark cutoff times for both conventional and modified landmark analysis were predefined at 9, 12, and 15 weeks. The 9-week cutoff was chosen based on published literature,17 and the 12 weeks based on the ESPAC-3 trial which assessed the optimal timing of ACT initiation.18,19 Finally, the landmark cutoff time at 15 weeks was also chosen based on the ESPAC-3 trial but with an additional grace period of 25% of the time, resulting in 3 additional weeks.

Statistical analysis

For descriptive analyses, the median time and interquartile range (IQR) and mean time and SD between diagnosis and ACT initiation were calculated. Additionally, density functions of the distribution of this time interval were computed and plotted for every registry.

The association between ACT exposure and overall survival was described for all analytical approaches using hazard ratios (HR) and corresponding 95% confidence intervals (CI) derived from multivariable Cox proportional hazards regression. All resected patients were followed from diagnosis (time zero) until death due to any reason, censoring, or end of follow-up.

All analyses were conducted for each country individually and adjusted for the confounders and prognostic factors listed in Table S1. Additionally, analyses were carried out using the pooled data of all countries using all overlapping adjustment variables. To account for potential correlation due to clustering effects of selected patients from the same country, a robust sandwich estimation was used.20 This was repeated in a sensitivity analysis using a shared frailty model, introducing random effects.21 An additional sensitivity analysis was carried out for patients with known pathological TNM staging.

Moreover, the overall amount of possible misclassified and unobserved person-time and events coming from each analytical approach was analyzed and quantified. Therefore, “misclassified person-time” was computed as the sum of the wrongly modeled (un)exposed person-time, and “misclassified events” as the sum of deaths occurring in these misclassified person-time intervals. The “unobserved person-time” was calculated as the sum of (un)exposed person-times before each respective landmark cutoff, and “unobserved events” as the sum of deaths that occurred before the respective landmark cutoff times.

The proportional hazards assumption was assessed by plotting weighted Schoenfeld residuals for each variable included in the regression models according to Grambsch and Therneau.22

Statistical analyses were conducted with R (R Foundation for Statistical Computing, Vienna, Austria)23–26 and SAS software, version 9.4 (SAS Institute Inc., Cary, NC, USA). Statistical significance was defined by a two-sided P<0.05. Simplified exemplary SAS codes for all analytical approaches are attached as an additional appendix in the Supplementary materials.

Results

In total, 2,449 (Belgium), 2,563 (the Netherlands), and 656 (Slovenia) resected PCa patients with median survival times (IQR) of 1.5 (0.8–3.2), 1.4 (0.8–3.0), and 1.1 (0.6–2.5) years, respectively, were eligible to be included in the study (Table 1).

Table 1

Basic characteristics of resected pancreatic cancer patients who were eligible to be included in the analysis

	Belgium		The Netherlands		Slovenia
	ACT		ACT		ACT
	Yes	No	Yes	No	Yes	No
Available years of diagnosis	2004–2013		2003–2014		2003–2013
Follow-up until	July 1, 2015		February 1, 2015		May 25, 2016
N (%)	1,394 (56.9%)	1,055 (43.1%)	909 (35.5%)	1,654 (64.5%)	174 (26.5%)	482 (73.5%)
Median time of ACT (IQR)	9 (7–12) weeks	NA	8 (6–11) weeks	NA	9 (7–12) weeks	NA
Mean time of ACT (SD)	10.2 (4.8) weeks	NA	9.2 (4.6) weeks	NA	10.6 (6.2) weeks	NA
Age (median [IQR])	64 (57–71)	71 (63–76)	64 (57–69)	69 (61–74)	62 (54–68)	67 (58–74)
Age groups (%)
<60	462 (33.1)	176 (16.7)	291 (32.0)	347 (21.0)	72 (41.4)	138 (28.6)
60–69	510 (36.6)	291 (27.6)	398 (43.8)	549 (33.2)	69 (39.7)	143 (29.7)
70–79	393 (28.2)	455 (43.1)	215 (23.7)	652 (39.4)	33 (19.0)	171 (35.5)
≥80	29 (2.1)	133 (12.6)	5 (0.6)	106 (6.4)	0 (0.0)	30 (6.2)
Sex (%)
Male	736 (52.8)	567 (53.7)	454 (49.9)	887 (53.6)	97 (55.7)	242 (50.2)
Female	658 (47.2)	488 (46.3)	455 (50.1)	767 (46.4)	77 (44.3)	240 (49.8)
Tumor location (%)
Head of pancreas	785 (56.3)	547 (51.8)	746 (82.1)	1,321 (79.9)	126 (72.4)	325 (67.4)
Body of pancreas	69 (4.9)	50 (4.7)	28 (3.1)	63 (3.8)	13 (7.5)	33 (6.8)
Tail of pancreas	129 (9.3)	76 (7.2)	72 (7.9)	132 (8.0)	11 (6.3)	33 (6.8)
Other/unspecified	411 (29.5)	382 (36.2)	63 (6.9)	138 (8.3)	24 (13.8)	91 (18.9)
Histology (%)
Acinar	13 (0.9)	14 (1.3)	2 (0.2)	13 (0.8)	0 (0.0)	3 (0.6)
Adenocarcinoma, NOS	661 (47.4)	494 (46.8)	501 (55.1)	991 (59.9)	35 (20.1)	104 (21.6)
Ductal	647 (46.4)	453 (42.9)	358 (39.4)	487 (29.4)	125 (71.8)	258 (53.5)
Intraductal	23 (1.6)	28 (2.7)	10 (1.1)	34 (2.1)	0 (0.0)	1 (0.2)
Mucinous	5 (0.4)	11 (1.0)	6 (0.7)	15 (0.9)	1 (0.6)	4 (0.8)
Other	45 (3.2)	55 (5.2)	32 (3.5)	114 (6.9)	13 (7.5)	112 (23.2)
Grade (%)
Well	189 (13.6)	154 (14.6)	91 (10.0)	153 (9.3)	15 (8.6)	61 (12.7)
Intermediate	613 (44.0)	401 (38.0)	370 (40.7)	718 (43.4)	57 (32.8)	131 (27.2)
Poor	388 (27.8)	271 (25.7)	293 (32.2)	467 (28.2)	78 (44.8)	174 (36.1)
Undifferentiated	18 (1.3)	19 (1.8)	2 (0.2)	4 (0.2)	3 (1.7)	8 (1.7)
Unknown	186 (13.3)	210 (19.9)	153 (16.8)	312 (18.9)	21 (12.1)	108 (22.4)
Spread (%)
Regional	1,186 (85.1)	760 (72.0)	816 (89.8)	1,284 (77.6)	140 (80.5)	351 (72.8)
Localized	94 (6.7)	230 (21.8)	83 (9.1)	309 (18.7)	2 (1.1)	60 (12.4)
Distant	114 (8.2)	65 (6.2)	10 (1.1)	61 (3.7)	32 (18.4)	71 (14.7)
Radiotherapy (%)	185 (13.3)	40 (3.8)	5 (0.6)	7 (0.4)	14 (8.0)	0 (0.0)
Previous other tumors (%)	183 (13.1)	159 (15.1)	106 (11.7)	253 (15.3)	NA	NA

Abbreviations: ACT, adjuvant chemotherapy; IQR, interquartile range; NOS, not otherwise specified; SD, standard deviation; NA, not available.

The proportion of patients who received ACT following resection was highest in Belgium (56.9%), followed by the Netherlands (35.5%) and Slovenia (26.5%). The median age at diagnosis ranged from 62 (ACT patients Slovenia) to 71 (non-ACT patients Belgium) years. Patients who received ACT were, on average, approximately 5 years younger than patients not receiving ACT. The majority of tumors were located at the head of the pancreas (51.8%–82.1%), were intermediately differentiated (27.2%–44.0%), and exhibited regional spread (72.0%–85.1%). Radiotherapy administration was least often observed in the Netherlands (0.4%–0.6%), whereas 13.3% of the patients in Belgium who received chemotherapy were also radiated.

On average, patients initiated ACT between 9.2 (SD=4.6; the Netherlands) and 10.6 (6.2; Slovenia) weeks. The median (IQR) time from diagnosis to ACT initiation was found to be lower and quite consistent throughout all registries at 8 (6–11; the Netherlands) and 9 (7–12; Belgium and Slovenia) weeks, respectively (Figure 2). The probability density functions showed that, in all registries, the distribution of the time from diagnosis to ACT initiation was right skewed, indicating that the majority (75%) of PCa patients initiated chemotherapy within 11–12 weeks or less, but the upper quartile was distributed more broadly toward longer time intervals.

Figure 2

Distribution of time between diagnosis and ACT initiation among resected pancreatic cancer patients (by country).

Abbreviation: ACT, adjuvant chemotherapy.

The results of the different analytical approaches for each country and overall for the study population are summarized in the forest plots displayed in Figure 3.

Figure 3

HR and corresponding 95% CI to assess the association between ACT and overall survival among resected pancreatic cancer patients calculated with different study design choices by country (yellow squares indicate the inverse of the variance of the estimate).

Abbreviations: ACT, adjuvant chemotherapy; HR, hazard ratio.

When the analyses were carried out by country, the HRs (95% CI) were generally lower in the Netherlands compared to Belgium and Slovenia, but the same patterns were observed consistently throughout each country (Figure 3). As expected, the immortal time-biased analyses resulted in much lower HRs than the ones derived using the Mantel–Byar method, with HRs of 0.69 vs 0.86 (Belgium), 0.59 vs 0.69 (the Netherlands), and 0.73 vs 0.86 (Slovenia), thereby overestimating protective effects of ACT approximately two-fold in Belgium and Slovenia. Results derived from the conventional landmark analysis were generally more conservative than the Mantel–Byar method, but performed especially well in the Dutch database and yielded the same result as the Mantel–Byar method when using the landmark time at 12 weeks. In comparison, results derived from the modified landmark approach resulted in HRs (95% CI) which were, in general, marginally smaller but still close to the results from the Mantel–Byar method and increased consistently with increasing landmark cutoff time. Whereas the estimates at 9 weeks were approximately half-way between the immortal time bias and the Mantel–Byar method results, the landmark cutoff time at 15 weeks provided nearly the same results as the Mantel–Byar method (HRs [95% CI] Mantel–Byar method vs modified landmark at 15 weeks of 0.86 [0.78–0.96] vs 0.87 [0.78–0.97] for Belgium, 0.69 [0.61–0.78] vs 0.67 [0.59–0.76] for the Netherlands, and 0.86 [0.70–1.06] vs 0.85 [0.68–1.06] for Slovenia, respectively).

When all countries were analyzed as a pooled dataset, overall the same results and patterns were found (Figure 4). Again, the immortal time-biased analysis provided a much lower HR than the Mantel–Byar method (0.68 vs 0.82) and resulted in a total of 469 misclassified exposed person-years (4.3% of total follow-up time, Table S2). The modified landmark analyses yielded HRs that were much less biased (0.76, 0.78, and 0.81, respectively), especially when setting the landmark cutoff at 15 weeks (HR [95% CI] 0.82 [0.71–0.93] Mantel–Byar vs 0.81 [0.70–0.94]). Here, the amount of misclassified exposed person-years (% of total follow-up) decreased with later landmark cutoff times (97 [1.0%], 48 [0.5%], and 26 [0.3%] with a 9-, 12-, and 15-week cutoff time, respectively). Results derived from the conventional landmark analysis were again found more conservative but with less variance between the chosen landmark times compared to the modified landmark approach. However, with 2,396 (58.8%), 1,161 (12.0%), and 489 (5.2%) person-years (% of total follow-up), the amount of misclassified unexposed person-years exceeded the amount of misclassified exposed person-years in the modified landmark analysis. Repeating the analyses using a shared frailty model or only including patients with a recorded pathological tumor-node-metastasis (pTNM) stage did not meaningfully alter the results (Figures S1 and 2).

Figure 4

HR and corresponding 95% CI to assess the association between ACT and overall survival among resected pancreatic cancer patients calculated with different study design choices overall (yellow squares indicate the inverse of the variance of the estimate).

Abbreviations: ACT, adjuvant chemotherapy; HR, hazard ratio.

Discussion

In this methodologically focused population-based cohort study, we aimed to systematically analyze the comparative performance of a modified landmark approach relative to the Mantel–Byar method as the gold standard by applying these analytical approaches to a well-known positive association between ACT among resected PCa patients and overall survival. As expected, all analyses showed a significantly increased survival among resected PCa patients who received ACT relative to surgery only. Compared to the Mantel–Byar method, immortal time-biased analyses overestimated the reduction in mortality for patients treated with ACT almost two-fold, whereas the conventional landmark approach yielded consistent and precise, but more conservative, estimates compared to the Mantel–Byar method. In comparison, results derived from the modified landmark approach resulted in HRs similar to the Mantel–Byar method, and were especially accurate when the landmark cutoff time was set at 15 weeks.

To our knowledge, this is the first methodological study to investigate the comparative performance of a modified landmark approach in a population-based setting. So far, published studies on that topic have rather focused on the conventional landmark approach which, however, requires time of treatment information.6,27–30 Given that this information is often not available and the modified version of the landmark approach is usually intuitively used,31 exploring the performance of the derived results relative to the Mantel–Byar method is of high relevance.

Recently, studies on PCa, which investigated the same or similar research questions as in this study, have been found to be affected by immortal time bias.32–34 Especially, cohort studies in the adjuvant setting are prone to be immortal time biased,35 but studies in the neoadjuvant setting can be biased as well. For example, a recently published study assessed the overall survival among pancreatic head cancer patients who received neoadjuvant treatment followed by surgery as compared to patients with surgery only.34 The authors found a 28% lower mortality among those patients with neoadjuvant chemotherapy followed by resection, which might be a result likely to be influenced by immortal time bias.36 As patients who received neoadjuvant treatment had to remain alive at least until the time of neoadjuvant treatment, they were artificially assigned a survival advantage compared to those who received surgery only. As the authors modeled the exposure as a time-fixed variable, this led to a biased rate ratio in favor of the neoadjuvant treatment. Again, approaches to avoid immortal time bias would have been to model the time between diagnosis and surgery as unexposed and exposed thereafter (Mantel–Byar), or to define a landmark at the time of surgery and follow the patients, who were still alive at that time, from this landmark cutoff.

In this study, the Mantel–Byar method provided results similar to results from RCTs, also indicating a significant survival benefit among resected PCa patients receiving ACT.9 The observed patterns were consistent in all registries, but with lower HRs in the Dutch cohort, which might be likely due to strong changes in treatment rates in the last decade. Expectedly, immortal time-biased analyses led to results which were in comparison to the Mantel–Byar method almost two-fold lower. Given that all eligible patients were chosen conditional on having received surgery, the only possible “immortal” time window was the time between surgery and potential chemotherapy initiation which, from a clinical perspective, is a rather short time window. However, given that the magnitude of immortal time bias is higher the worse a prognosis of the underlying disease is,37 the results presented herein are sensitive even for short time intervals, because survival among PCa patients is low (as displayed in the crude 1- and 5-year survival estimates in Table S3).

However, one has to be aware of limitations coming with a modified landmark approach. A critical decision in conducting a modified landmark analysis is the choice of the landmark cutoff time. In this study, landmark cutoff times of 9, 12, and 15 weeks were chosen based on published literature and information derived from RCTs. It was found that adding a 25% grace period to the cutoff time suggested by RCTs led to the best performance relative to the Mantel– Byar method. Therefore, the choice of the best cutoff time has to be pre-evaluated for every specific research question and cancer entity. Good sources for potential landmark times might be expert opinions, published literature, and guidelines, if available. Methodological studies investigating the conventional landmark approach suggested to present results using a sequential landmark analysis, which means that results from multiple cutoff times are reported.29 Accordingly, the sequential analysis can also be recommended to present results derived from a modified landmark analysis, which should always go along with a critical interpretation of the results.

Moreover, taking the average time of treatment initiation as a landmark cutoff – which could be available on an aggregated data level from cancer registries – might be applicable in the conventional landmark approach due to the generally higher precision, but not recommendable in the modified landmark approach as the results with the average time to ACT initiation (consistently between 9 and 12 weeks) overestimated the potential true association. This might be explained by the right-skewed distribution of the time to ACT initiation which was consistently observed throughout all registries. However, if the time of treatment information is available anyway on an individual patient level, the Mantel–Byar method should always be preferred over the conventional landmark analysis. The only potential advantage of the conventional landmark analysis might be the easier computation compared to the more complex programming of the Mantel–Byar method, but this approach was found to be biased if treatment effects are large.6

When, in the conventional landmark approach, a patient starts the exposure after the landmark time, the patient will be modeled as unexposed throughout the whole follow-up, which explains the conservative estimates in our analyses. Using the modified landmark approach, patients might be modeled as exposed at the landmark time, although the patient’s treatment has, in reality, started sometime later, due to the absence of treatment timing information. This is why it is important to emphasize that the modified landmark approach only mitigates – but does not completely remove – immortal time bias, because the gap of the expected and the actual treatment initiation among some patients still leads to a small immortal time gap. However, the quantification of this possible immortal time interval showed that the amount of misclassified exposed person-time decreased with later landmark cutoff times to an almost negligible proportion of 0.3% of the total follow-up time at the 15-week landmark cutoff, which explains why the HRs derived from the modified landmark approach were generally a bit lower and steadily increased with increasing the landmark time (Table S2). Therefore, a good balance in the choice of the landmark time has to be found, because the later a landmark cutoff time is defined, the less possible it is that “immortal” person-time is included but also that more patients will be excluded (unobserved) from the analysis. This not only leads to decreased statistical power, but may also result in selection bias as many patients with an early event are being excluded, which is problematic (e.g., in cancer entities with a bad prognosis). A recently published simulation study, which investigated the performance of the conventional landmark approach, also found that the bias and the precision depend on the choice of the landmark time as well as the event and treatment rates, given that there is a large treatment effect.6 The authors concluded that with varying landmark times, the differences in results were small; however, in general, longer landmark times performed better but with decreasing precision over time when the treatment rate was high in early follow-up. This once more emphasizes the need of a careful consideration of an appropriate landmark cutoff time and underlines the recommended sequential approach using multiple landmark times.

In other scenarios than the one illustrated in this study, the modified landmark approach might also perform differently, depending on the prognosis of each cancer entity, the effectiveness of the treatment intervention, and the time when the treatment is started during the course of the disease. In a previous study37 it was found that the magnitude of immortal time bias conclusively correlates with the prognosis of each cancer entity; with a larger impact, the worse the prognosis was (with PCa as the cancer entity with the worst prognosis and, hence, the largest bias). This implies that in cancer entities with better prognosis, such as breast cancer or prostate cancer, the residual misclassified exposed person-time in the modified landmark analysis will most likely have less impact on the derived estimates, and the differences between landmark cutoff times will most likely decrease. In another scenario where the treatment effect is rather small, little variations in the landmark time might not reveal the potential treatment association. Additionally, if the landmark cutoff is set very late during the course of the disease, results might not be generalizable anymore as the landmark approach implies a conditional survival to the time of the landmark cutoff.38

This study has strengths and limitations. Firstly, residual confounding by indication cannot be ruled out, because guidelines have increasingly recommended chemotherapy whenever possible among resected PCa patients, which means that resected patients who have not received post-surgical chemotherapy might have been too ill to receive chemotherapy or refused chemotherapy due to other reasons. However, this potential bias would be expected to apply to all of the analytical approaches in a similar manner, and should not have had a major impact on our results. Secondly, there was not sufficient information available to control for comorbidities or concomitantly used drugs, and the results may not be generalizable as no simulations were undertaken. However, due to the bad prognosis of PCa, deaths due to competing risks and further analyses on cancer-specific survival might be negligible.

Nevertheless, it was possible to use data from multiple cancer registries from different countries which was done in order to emulate a multicentric RCT and to decrease the chance of random false-positive findings compared to only using one database. Additionally, with the provided data, the most important confounders and prognostic factors could be controlled for in multivariable regression analysis.

Another strength is that the study question investigated here was undertaken in a very sensitive setting and was, therefore, found ideal due to the following reasons. Firstly, the evidence of major RCTs on this question is quite strong, serving as a good estimate of the “true” treatment association by ACT. Secondly, resected PCa patients represent a small but homogeneous group of patients, which provides a better basis to study treatment associations as there are less influential factors that might distort a potential association.39 Thirdly, the ESPAC-3 trial revealed that survival of PCa patients does not depend on the time of initiation of ACT after surgery, which would have introduced bias had there been a causal relationship.18,19 Finally, recommendations were recently made to combine adjuvant gemcitabine and capecitabine (ESPAC-4 trial); however, in general, PCa unfortunately still has a bad prognosis and no major therapeutic breakthroughs were approved in the last decade, which means that only small time windows can yet cause severe immortal time bias in this setting.37,40 As the modified landmark approach with a landmark time at 15 weeks led to nearly the same results as the Mantel–Byar model, this strengthens confidence in the modified landmark approach as a useful analytical tool which might also be used in studies with cancer entities other than PCa.

In conclusion, this PCa-focused study provides some evidence that a modified landmark approach leads to results that are close to the ones derived by the Mantel–Byar method, which is considered the gold standard analytical approach to avoid immortal time bias. The information on time of exposure initiation is often lacking in many oncological databases, especially in larger cancer registry databases such as those in the US (Surveillance, Epidemiology, and End Results Program [SEER]) or Germany.41 This gap in the linkage to detailed and high-quality treatment data for research purposes might be due to restricted access to such data, logistical difficulties, or data protection policies. The conduct of this study was motivated by the large amount of studies that -although not having detailed information on the time of exposure initiation-tried to investigate treatment associations, which has frequently led to tremendously (immortal time) biased results.42 The modified landmark approach provides a fast and easy to use way to meaningfully mitigate immortal time bias when the time of exposure initiation is not recorded in pertinent databases. If researchers are aware of the potential factors that might influence a modified landmark analysis and critically interpret the derived results, this analytical approach might be a valuable tool to assess treatment associations in databases that do not provide any time of treatment data. Due to the very easy computational realization and application of this method, validation of these findings and more experience with the modified landmark analysis in other cancer entities and indications is of utmost interest for future research.

Supplementary materials

Simplified exemplary SAS code for different analytical approaches including conventional landmark approach (A1), modified landmark approach (A2) and time-dependent Mantel–Byar analysis (A3)

/*Documentation exemplary dataset (“cohort”): */

/*The exemplary dataset of resected pancreatic cancer patients with all variables is called “cohort”*/

/*fu_time=follow-up time (in weeks)*/

/*death=event indicator (1=death, 0=censoring)*/

/*t_cx=time from diagnosis to adjuvant chemotherapy start*/

/*start_cx=week of follow-up, in which adjuvant chemotherapy was initiated*/

/*stop_cx=week of follow-up, in which adjuvant chemotherapy was terminated (in this study end of follow-up since the exposure was modeled according to a first-treatment-carried-forward (intention-to-treat) analysis*/

/*cht=chemotherapy (1=yes, 0=no)*/

/*age=age (continuous)*/

/*sex=sex (male, female)*/

/*diag_year=year of diagnosis*/

/*loc=location of tumor (head, body, tail, other)*/

/*hist=histology (adenocarcinoma, ductal, etc.)*/

/*staging=stage at diagnosis (I, II/III, IV)*/

A1. Conventional landmark approach

/*Conventional landmark approach - cut-off time 1 (9 weeks)*/

data conv_landmark1_cohort;

set cohort;

if start_cx ne . and start_cx < 9 then cht=1; else cht=0;

if fu_time < 9 then delete;

fu_time1=fu_time-9;

run;

proc phreg data= conv_landmark1_cohort;

class sex diag_year loc hist staging cht(ref=”0”);

model fu_time1*death(0) = age sex diag_year loc hist

staging cht/ties=efron rl;

run;

/* Conventional landmark approach - cut-off time 2 (12 weeks)*/

data conv_landmark2_cohort;

set cohort;

if start_cx ne . and start_cx < 12 then cht=1; else cht=0;

if fu_time < 12 then delete;

fu_time2=fu_time-12;

run;

proc phreg data= conv_landmark2_cohort;

class sex diag_year loc hist staging cht(ref=”0”);

model fu_time2*death(0) = age sex diag_year loc hist

staging cht/ties=efron rl;

run;

/* Conventional landmark approach - cut-off time 3 (15 weeks)*/

data conv_landmark3_cohort;

set cohort;

if start_cx ne . and start_cx < 15 then cht=1; else cht=0;

if fu_time < 15 then delete;

fu_time3=fu_time-15;

run;

proc phreg data= conv_landmark3_cohort;

class sex diag_year loc hist staging cht(ref=”0”);

model fu_time3*death(0) = age sex diag_year loc hist

staging cht/ties=efron rl;

run;

A2. Modified landmark approach

/*Modified landmark approach - cut-off time 1 (9 weeks)*/

data mod_landmark1_cohort;

set cohort;

if fu_time < 9 then delete;

fu_time1=fu_time-9;

run;

proc phreg data= mod_landmark1_cohort;

class sex diag_year loc hist staging cht(ref=”0”);

model fu_time1*death(0) = age sex diag_year loc hist

staging cht/ties=efron rl;

run;

/*Modified landmark approach - cut-off time 2 (12 weeks)*/

data mod_landmark2_cohort;

set cohort;

if fu_time < 12 then delete;

fu_time2=fu_time-12;

run;

proc phreg data= mod_landmark2_cohort;

class sex diag_year loc hist staging cht(ref=”0”);

model fu_time2*death(0) = age sex diag_year loc hist

staging cht/ties=efron rl;

run;

/*Modified landmark approach - cut-off time 3 (15 weeks)*/

data mod_landmark3_cohort;

set cohort;

if fu_time < 15 then delete;

fu_time3=fu_time-15;

run;

proc phreg data= mod_landmark3_cohort;

class sex diag_year loc hist staging cht(ref=”0”);

model fu_time3*death(0) = age sex diag_year loc hist

staging cht/ties=efron rl;

run;

For the Mantel–Byar analysis there are several approaches in SAS which lead to the same results as long as the variable for the follow-up time is expressed in the same unit and rounded the same way (e.g., in weeks). If one wants to learn more, the following resources are recommended for further reading: Allison, P. D. (2010). Survival Analysis Using SAS: A Practical Guide: Sas Inst.

Powell T, Bagnell M. Your “Survival” Guide to Using Time-Dependent Covariates. SAS Global Forum 2012 2012.

A3.1. Time-dependent Mantel–Byar analysis (programming statement approach)

proc phreg data= cohort;

class sex diag_year loc hist staging;

model fu_time*death(0) = age sex diag_year loc hist

staging ACT /ties=efron rl;

if cht = 1 and fu_time >= start_cx then ACT=1; else

ACT=0 ;

run;

A3.2. Time-dependent Mantel–Byar analysis (array approach)

/*Using the Mantel–Byar array approach some preliminary work has to be performed first*/

/*Let’s assume the longest follow-up for a patient in our dataset “cohort” is 686 weeks*/

*Compute maximum follow-up time to create time intervals; title “Average/Median time from diagnosis to death/end of fu (in weeks)”;

proc means data=cohort N MEAN MIN MAX MEDIAN

Q1 Q3 QRANGE NMISS;

var fu_time;

run;

title;

/*Let’s assume the maximum follow-up was 686 weeks => create 686 (yes/no) variables indicating exposure for each week*/

*Compute statistics for time to adjuvant chemotherapy initiation (optional);

title “Average/Median time from diagnosis to Cx initiation of fu (in weeks)”;

proc means data=cohort N MEAN MIN MAX MEDIAN

Q1 Q3 QRANGE NMISS;

var t_cx;

run;

title;

/*Prepare time-dependent variables and initially set all to 0*/

%macro dx_vars(intv);

data cohort;

set cohort;

Cx_t&intv.=0; /*Cx_t&intv. Indicates exposure in each week for every patient for the whole follow-up*/

run;

%mend;

%macro loop_intake();

%do i = 0 %to 686;/*MAX FU was 686 weeks*/

%(&i.);

%end;

%mend;

%();

/*Check plausibility!*/

data cohort;

set cohort;

if start_Cx ne. and stop_Cx ne. and start_Cx>stop_Cx

then start_stop_error=”Yes”;

else start_stop_error=”No”;

label start_stop_error=”Start date after stop date (error)”;

run;

proc freq data= cohort;

table start_stop_error;

run;

/*Set all time intervals to 1 between start and stop date*/

%macro dx_vars(intv);

data cohort;

set cohort;

if start_Cx <= &intv. <= stop_Cx then Cx_t&intv.=1;

run;

%mend;

%macro loop_intake();

%do i = 0 %to 686;

%(&i.);

%end;

%mend;

%();

/*Now the actual analysis (Cox proportional hazards regression) can be performed*/

proc phreg data= cohort;

class sex diag_year loc hist staging;

model fu_time*death(0) = age sex diag_year loc hist

staging cht Cx_time/ties=efron rl;

array Cx_array(*) Cx_t0-Cx_t686;

Cx_time = Cx_array[fu_time];

run;

Hazard ratio (HR) and corresponding 95% confidence interval (95% CI) to assess the association between adjuvant chemotherapy (ACT) and overall survival among resected pancreatic cancer patients calculated with different study design choices overall using a shared frailty model (sensitivity analysis) (yellow squares indicate the inverse of the variance of the estimate).

Hazard ratio (HR) and corresponding 95% confidence interval (95% CI) to assess the association between adjuvant chemotherapy (ACT) and overall survival among resected pancreatic cancer patients calculated with different study design choices when only including patients with pathological TNM staging (sensitivity analysis; yellow squares indicate the inverse of the variance of the estimate).

Kaplan-Meier plots (life-table [acturial] method) for (crude) 1- and 5-year survival among all resected PaC patients.

Table S1

Availability of confounders and important prognostic factors to adjust for in multivariable Cox proportional hazards regression (by country)

Variables adjusted for*	BCR (Belgium)	NCR (Netherlands)	CRS (Slovenia)
Age (at diagnosis)	Y	Y	Y
Sex (male/female)	Y	Y	Y
Calendar year of diagnosis (2003-2014)	Y	Y	Y
Topography/tumor location (body, head, tail, other)	Y	Y	Y
Morphology (adenocarcinoma NOS, acinar, ductal, intraductal, mucinous, other)	Y	Y	Y
Staging/grouping (I, II/III, IV)	Y (7.0% missing)	Y (0.5% missing)	Y (1.9% missing)
Grade (well, intermediate, poor, undifferentiated)	N (16.1% missing)	N (18.2% missing)	N (19.7% missing)
Comorbidities (score)	N (24.2% missing)	N (85.2% missing)	NA
Other/former tumors (yes, no)	Y	Y	NA
Harvested lymph nodes (<8, 9-12, 13-20, >20	NA	Y (4.9% missing)	NA

Notes:

Defined a-priori. Y: adjusted for; patients with missing values were excluded from analysis. N: not available or not adjusted for due to too many missing values.

Abbreviations: BCR, Belgium Cancer Registry; CRS, Cancer Registry of Slovenia; NCR, Netherlands Cancer Registry; NOS, not otherwise specified; NA, not available.

Table S2

Quantification of potential unobserved and misclassified person-times and events through different analytical approaches.

Method	Exposed			Unexposed
	PY/events	Misclassified PYa/eventsb	Unobserved PY/eventsc	PY/events	Misclassified PYa/eventsb	Unobserved PY/eventsc
Immortal time bias	5073/1785	469/−	−/−	5848/2588	−	−/−
Mantel–Byar	4602/1785	−/−	−/−	6319/2588	−	−/−
Conventional landmark
9 weeks	2154/813	−/−	60/6	7814/3265	2396/782	893/289
12 weeks	3295/1303	−/−	154/11	6367/2692	1161/332	1105/367
15 weeks	3847/1565	−/−	274/23	5516/2340	489/145	1284/445
Modified landmark
9 weeks	2154/813	97/−	60/6	7814/3265	−	893/289
12 weeks	3295/1303	48/−	154/11	6367/2692	−	1105/367
15 weeks	3847/1565	26/−	274/23	5516/2340	−	1284/445

Notes: PY=person-years,

Misclassified (un)exposed person-years were calculated for the exposed person-years that were misclassified as exposed (immortal time bias analysis), or as the unexposed person-years that were misclassified as unexposed between start of chemotherapy and end of follow-up (conventional landmark analysis), or as the person-years that were misclassified as exposed between the landmark time and the actual start of chemotherapy (modified landmark analysis),

Misclassified events were defined as events that occurred during a period of follow-up time where patients were misclassified as unexposed (conventional landmark analysis),

Unobserved person-years/events were calculated as the (un)exposed person-times/events before the respective landmark time (including those patients who deceased before the landmark time).

Table S3

Crude 1- and 5-year survival rates for all resected PaC patients by country, and overall

Cancer registry	1-year survival (%) (95% CI)	5-year survival (%) (95% CI)
BCR (Belgium)	66.5 (64.6–68.4)	16.4 (14.8–18.0)
NCR (Netherlands)	64.9 (62.9–66.7)	15.2 (13.5–17.0)
CRS (Slovenia)	55.6 (51.7–59.3)	16.5 (13.6–19.6)
Overall	64.5 (63.2–65.8)	16.0 (14.9–17.1)

Abbreviations: CI, Confidence interval; PaC, Pancreatic cancer; BCR, Belgium Cancer Registry; CRS, Cancer Registry of Slovenia; NCR, Netherlands Cancer Registry.