Abstract
With the introduction of immunotherapy with non-small lung cancer, prognosis of these patients has improved. However, socio-economic differences in access to various immunotherapy treatments have been reported. In the Netherlands, such differences are not expected due to universal insurance coverage. We investigated the existence of differential susceptibility by socioeconomic status (SES) of the effect of distance to treatment hospital on access to Durvalumab in patients with stage III non-small lung cancer who received chemoradiation, and the influence of differential mortality. We used data from the Netherlands Cancer Registry (n = 3774) from the period 2017-2021. First, we fitted Bayesian discrete failure time models and compared SES-by-distance-to-hospital interaction to a baseline model including age, distance, SES and performance score. We then fitted a time to mortality model and used both models in a g-formula to simulate a scenario where mortality levels were equalized. Our results showed that the high SES group received Durvalumab more often than the low SES group (HR = 1.23, CI95% = [1.07, 1.53]), and showed a strong negative effect of distance (OR = 0.63, CI95% = [0.54, 0.72]). Bayes factor < 3 indicated inconclusive evidence for a SES by distance interaction effect, while g-formula results showed virtually no influence from differential mortality on SES differences. Secondary analyses showed strong evidence that SES differences in using Durvalumab were constant over the years (Bayes factor > 17). Overall, these results are significant for understanding how socio-economic inequality affects proper care and can be vital for public policy.
Introduction
Worldwide, lung cancer is the most commonly diagnosed cancer, and is the leading cause of cancer- related mortality.1 Incidence rates of lung cancer are higher among individuals with low socio-economic status (SES), and patients with low SES have poorer overall survival rates. Even after adjusting for smoking and comorbidities, this SES differential persists.2,3 Research shows that patients belonging to lower SES groups have less access to therapy overall, not only traditional but also next-generation therapies.4 In the Netherlands, health insurance coverage is virtually universal and does not include (economically differentiated) tiers of service for cancer therapy. This makes the Netherlands a useful case-study for causes of SES differences in access to various cancer therapies, as it eases investigation of non-institutional drivers of such differences.
Durvalumab is a recently developed immunotherapy that substantially increases progression-free survival in patients with stage III non-small cell lung cancer (NSCLC) after chemoradiotherapy.5 The therapy was introduced in 2017 in the Netherlands for clinical studies, where it became part of the insurance package in 2019. The number of hospitals prescribing Durvalumab increased over time. Durvalumab is an expensive therapy, with its cost estimated at about 61,000 euros per patient in 2019. 6
Two potential mechanisms of interest are travel distance to treatment facility and differential mortality. Firstly, while individuals with a lower SES might be eligible for treatment with Durvalumab, hospitals that provide Durvalumab treatment might be located far enough away that overcoming this distance is a greater burden to individuals with less economic and social capital to bridge the distance. Hence, these individuals might opt for more traditional treatment in a less distant hospital. This would be a case of differential susceptibility by SES7. Secondly, as socio-economic differentials in lung cancer survival also exist, individuals with a higher SES might survive long enough to receive Durvalumab treatment, whereas death would be a relatively stronger competing risk for individuals with a lower SES.
The primary aim of this study is to determine whether a possible effect of distance to Durvalumab treating hospital on use of Durvalumab is moderated by SES for stage III non-small cell lung cancer (NSCLC) patients diagnosed in the period 2017-2021 in the Netherlands. To determine this, we must also establish whether SES differentials in Durvalumab usage exist, and whether a distance effect exists. The secondary aim of this study is to determine to what extent possible SES differences in use of Durvalumab is influenced by differential survival between SES groups. In addition, we evaluate whether SES modifies how use of Durvalumab unfolds across time (based on year of diagnosis) per SES group.
Data and Methods
Data source
We extracted information for lung cancer patients from the Netherlands Cancer Registry (NCR), which is a population-based nationwide registry keeping records of each individual cancer patient since 1989 in the Netherlands. The registry contains demographic, tumor and primary treatment information. Information on vital status and date of death is updated annually using a computerized link with the national civil registry, and complete until Feb 1st, 2023.
Net median household income at 6-digit postal code level in 2016 was used as proxy for SES, obtained from Statistics Netherlands. SES was clustered into 9 groups, covering 99% of the postal codes of the Dutch population. Each postal codes covers on average 17 households.
Study population and sample
In the NCR, we identified 49504 patients with non-small cell lung cancer (NSCLC) diagnosed in the years 2017 –2021 . Following clinical trials on Durvalumab effectiveness, we restricted the sample to patients with clinical stage III cancer8 who received at least one cycle of chemotherapy and radiation therapy before potentially being administered Durvalumab (n=11114).5 Furthermore, patients with surgery were excluded from the sample (n = 3345). We further excluded patients with WHO performance status score higher than 1 (n = 729) or with missing performance status score (n = 1208). This resulted in a sample with patients able to carry out normal unrestricted activities, or restricted in strenuous activity but fully ambulant patients, because Durvalumab is given to relatively healthy patients. Patients who ended up in a different (immuno)therapy regime than Durvalumab after chemoradiotherapy were also excluded (n=2040). Finally, we further excluded patients who are registered as receiving Durvalumab more than once (n = 18) as this likely indicates a registration error. The final sample size was 3774 patients.
Outcome, censoring and competing risk
Our outcome variable is time from diagnosis until first treatment with Durvalumab, binned in 25 day units. Patients who do not receive Durvalumab by the end of follow-up are censored at the time of loss to follow- up. Patients who died before receiving Durvalumab instead are categorized as having died and in the model for time-to-durvalumab are, like censored patients, treated as has having left the risk-set at the time of loss-to-follow up.
Primary and secondary exposures and effect modifier
Our primary exposures is SES. This variable was originally coded in deciles from 10% to 90% by Statistics Netherlands, and in this study was further categorized into 10% to 30%, 40% to 60% and 70% to 90%, representing low, medium and high SES respectively.
Our secondary exposure is minimum driving distance to the closest hospital providing Durvalumab (henceforth ‘driving distance’). Driving distance was calculated based on the postal code (PC4) of the patients’ residence at diagnosis and the postal code of hospitals. Availability of Durvalumab at each hospital was based on the earliest year that the therapy was recorded in the NCR. For cases of different hospital locations of the same hospital the average driving distance was taken.
We treat SES also as a potential effect modifier of the relationship between minimum-driving distance and time-to-Durvalumab treatment.
Potential confounders and clustering
We consider age at diagnosis, WHO performance status score and year of diagnosis as potential confounding variables. Age was measured in years. Performance status followed the WHO scale, where ‘0’ indicates asymptomatic patients and ‘1’ indicates symptomatic patients but fully ambulant. Year of diagnosis was based on semesters (defined from January to June, and from July to December) and treated as a continuous variable. Furthermore, we consider observations to be clustered within the hospital which the patients first visited with symptoms related to lung cancer. This is the hospital where the patients is usually diagnosed, and in this hospital the patient might receive a referral for treatment in another hospital. Hospital variation can potentially impact the type of treatment.
Statistical analysis
Our estimands of interest are the restricted mean survival times (RMSTs) if all individuals had low, medium, and high SES, and the differences between these quantities, as well as the marginal hazard ratios between these groups.
Furthermore, to gain insight into the influence of death as an outcome that would prevent receiving Durvalumab treatment, also known as a competing risk, we also approximated these estimands under a scenario where death distributions were equalized between SES groups. We set all groups to have the mortality levels, conditional on their other covariates, of the low SES group. These scenarios are estimated using a Bayesian g-formula with Monte Carlo integration.9
In order to approximate these quantities, we first assume the extended DAG shown in Figure 1. Secondly, following this DAG, we fit a Bayesian hierarchical discrete failure time model (DFT) with time to Durvalumab as the outcome variable, and the minimum driving distance, SES, and potential confounders as covariates. The STAN probabilistic programming language was used to approximate posterior distributions of quantities through Hamiltonian Monte Carlo sampling.10 Inclusion of an interaction term between SES and driving distance, and between diagnostic year and SES, were determined using Bayes factor with a factor of 3 as the cutoff for inclusion.11,12 Patients are assumed nested within hospitals of first contact. This procedure was also followed for a second DFT with time to death as the outcome variable. The modelling procedure is described in more detail in appendix A.
Subsequently, by taking draws from the posterior distributions of these two fitted models, we first simulated survival times and choose the earliest survival time (Durvalumab or death) as the one observed. This was initially done to produce the Natural Course scenario, an approximation of the empirical data using the observed covariate information as input. The estimated Kaplan-Meier curve of this scenario is compared to the empirical Kaplan-Meier curve to guard against gross model misspecification. Afterwards, this procedure was performed with input data where all individuals were counterfactually set to have low SES, then medium SES, and finally high SES. Based on this, restricted mean survival times and marginal hazard ratios were calculated (scenario A). Our inferences are based on 3000 draws from the posterior, and we performed 25 Monte Carlo iterations per posterior draw to reduce Monte Carlo error. Finally, this procedure was repeated, but now in a simulation where all individuals belonging to medium and high SES were given the death distribution corresponding to the low SES simulation in order to equalize the potential influence of this competing risk (scenario B). Unless otherwise indicated, reported point estimates refer to means taken over their respective posterior distribution.
Results
Table 1 shows the descriptive statistics of the variables included in the statistical models. Of interest is that 37% of included patients received Durvalumab and that 10% of patients died/were censored before they could be prescribed Durvalumab (i.e., within time window of 250 days since diagnosis after which patients are administratively censored). Most included patients, 54% lived in a neighborhood with a medium SES, as compared to low (27%) and high (17%). The average time to receive Durvalumab was 160 days. Within the time window of 250 days, the average time-to-death for our cohort was 168 days with about 10% of patients having died.
Table 2 presents the results of the model comparison to determine if interactions between SES and driving distance, and between SES and year of diagnosis, should be included in our final model. We found neither evidence for nor against the inclusion of the interaction term between SES and driving distance; the baseline model was 2.38 times more likely than the model with interaction (3 was the decision threshold). Furthermore, we found strong evidence against the inclusion of the interaction term between SES and year of diagnosis; the baseline model was 17.06 times more likely. As a result, we used the baseline model for the Monte Carlo integration to approximate our estimands of interests. Findings for the model for time-to-death were similar and are shown in appendix B. The comparison between the natural course simulation from the baseline model and the empirical Kaplan Meier curve showed some underestimation of survival in days 50 to 100, but no major divergences (appendix C).
Regarding the main effects of the baseline time-to-durvalumab discrete failure time model, the coefficient of driving distance indicated a negative association between SES and receiving Durvalumab (OR = 0.63, CI95% = [0.54, 0.72]), while that of year of diagnosis indicated a positive association (OR = 2.22, CI95% = [2.08, 2.38]). Additionally, our covariates of performance status and age both showed a negative association with respect to receiving Durvalumab (OR = 0.71, CI95% = [0.62, 0.80], and OR = 0.80, CI95% = [0.76, 0.85] respectively).
In line with the Bayes factor for inclusion, In the expanded model that included an interaction term between SES and driving distance, we found no evidence of an interaction, neither between low and high SES (OR difference = 0.10, CI95% = [-0.29, 0.52]), nor between low and medium SES (OR difference = -0.13, CI95% = [-0.43, 0.15]), nor between medium and high SES (OR difference = 0.23, CI95% = [-0.12, 0.60]).
Figure 2 shows the estimated Kaplan-Meier curves from the Monte Carlo integration based on the baseline models. Specifically, average proportions of population not yet treated with Durvalumab are shown under the two simulated scenarios; on the left the scenario A where everyone is set to low, medium and high SES; on the right the scenario B where all are set to low, medium and high SES for their estimated time-to-durvalumab, but all are set to low SES for their estimated time-to-death. After having equalized the death distribution of all SES groups, the effect of SES group remains very similar across the two scenarios.
To formally test the effect of SES group (corresponding to the average treatment effect of SES), we calculated Hazard Ratios (HR) and Restricted Mean Survival Time (RMST) based on the time-to- Durvalumab probability curves of the counterfactual scenario B simulation (see Table 3). Overall, both HRs and RMSTs suggested that the high SES group receives Durvalumab earlier than the low SES group. HRs and RMSTs based on scenario A derived very similar values to those shown in Table 3 from scenario B.
Discussion
Our primary research question was “To what extent does a potential effect of driving distance to the closest Durvalumab-offering hospital differ between SES groups in stage III NSCLC patients?”. Our findings neither support nor negate effect modification of SES on driving distance. Before establishing this finding, we first had to establish main effects of SES and driving distance. Firstly, we found a strong negative main effect of driving distance (OR = 0.63, CI95% = [0.54, 0.72]). Secondly, in the scenarios where all patients were set to high, medium and low SES we found that in the high SES scenario had a 26% higher hazard to receive Durvalumab than in the low SES scenario (CI95% 1.07 to 1.53). ). As a secondary analysis, we determined whether the effect of SES on Durvalumab can be (partly) explained by the differential early mortality of lower income earners, which could, in turn, prevent use of Durvalumab. Our findings suggest that the effect of SES on using Durvalumab after equalizing the death distribution of the different SES groups (scenario B) was very similar to the effect given the observed death distributions (scenario A). Finally, although we found an overall increase of the usage of Durvalumab from 2017 to 2021 by about 60% (OR = 2.22, CI95% = [2.08, 2.38]), we found strong evidence against effect modification of SES on year of diagnosis, indicating no evidence for closing or increasing the gap between high SES and low SES over the years. Overall, these findings suggest that the direct effect of SES on using Durvalumab cannot be explained away by either the differential effect of driving distance to the closest center offering Durvalumab, or by differential survival.
Evaluation of data and methods
Our study has a number of strengths. Firstly, internal validity is strong as the number of relevant confounding variables is limited; conditional on SES, mobility decisions tend to be taken with distance to work or family in mind and not with distance to a hospital (specifically hospital treating with Durvalumab) in mind.13 Hence, few variables exist that influence distance to such a hospital while also affecting treatment decisions regarding Durvalumab. The number of hospitals that provide Durvalumab treatment was greatly expanded in our data set, starting at 10 hospitals in 2017 to 72 in 2021, strengthening the exogeneity of the distance variable. In addition, we covered all confounding variables following the medical guidelines in the Netherlands on which patients may receive Durvalumab and who may not. In the same vein, the assumed DAG shows no collider bias as the outcome cannot cause other variables for the simple reason that receiving Durvalumab temporally comes later than any other variable included in the DAG.
Secondly, we consider the functional form of our DFT model as an additional strength. Random effects, through partial pooling, allow the exchange of information between the different hierarchical levels of clustered data.15 In addition, we provided a flexible functional (natural spline) form of the time predictor to account for the non-linear effect of time.14 This spline effectively models the equivalent of what would be the baseline hazard curve in other common time-to-event models. We opted for linear relationships for the rest of the covariates and the secondary exposure for simplicity as well as after having visually examined the resulted (quasi-linear) form of such relationships when applying a spline transformation to them.
Finally, death is a competing risk for time-to-Durvalumab. By censoring patients at the time of death, it is implicitly assumed that they would have the same hazards of receiving Durvalumab as another individual with the same covariate values still alive at that time. However, this represents an untestable assumption. For this reason, in our equalized death scenario (scenario B), we gave all individuals the survival probabilities of individuals with low SES, conditional on their other covariate values; this represents a scenario in which we made medium and high SES individuals die sooner, rather than later, thereby avoiding the aforementioned untestable assumption.
An additional possible limitation of our work consists of the exclusion of a sizeable (n=2040) number of patients who were given a different (immuno)therapy regime than Durvalumab after receiving chemoradiotherapy . Systematic variation of ending up in either of these therapy regimes could be a potential bias for our results: for example, if patients belonging to the low SES group were systematically given other (immuno)therapy regimes than Durvalumab (and vice-versa). Given that the choice between other (immuno)therapies and Durvalumab is mutually exclusive, future research could go a step further by equalizing the reception of immunotherapy distribution (thereby, treating the reception of other (immuno)therapies as a competing risk). However, at this time there is no evidence that the low SES group has higher usage (than other SES groups) of other types of immunotherapy.
Finally, a limitation concerns the follow-up period studied; patients are selected if they completed at least one round of chemotherapy and radiation therapy. However, our follow-up time starts at diagnosis. This effectively means that patients are functionally immortal until they finish their first round of chemotherapy and radiation therapy15. Differential delays in reaching these inclusion conditions would likely affect timing of first Durvalumab treatment as well. In this paper, we effectively study the accumulation of SES differences in both such a potential delay and the one to Durvalumab. Future research should determine the degree to which SES differences in time-to-Durvalumab are explained by initial differences in time-to-first chemotherapy and radiotherapy between the SES groups.
Connection to previous studies
Our study showed patients of higher SES tend to have more use of immunotherapy for NSCLC. This is in line with previous studies showing similar results for treatments in other cancer types.16,17 Research has shown the general trend of elevated end-of-life costs for higher SES groups.18
In addition, our finding of a negative association between administration of Durvalumab and driving distance to the closest Durvalumab-offering hospital can be seen in the broader context of regional variations of cancer treatment in the Netherlands. For example, regional variations have been observed in the use of radiotherapy for prostate cancer19 and rectal cancer20. Specifically for NSCLC, recent evidence suggests that regional variations in the Netherlands are associated with differential use of radiotherapy and type of chemoradiotherapy (i.e., sequential vs. concurrent).21 Finally, our data remain inconclusive regarding modification of the effect of driving distance by SES, as shown by Bayes Factors and parameter estimates. A larger sample size is required to establish a conclusion with confidence.
Further possible pathways
There are many possible pathways through which SES affects use of Durvalumab. For example, individuals with a higher SES tend to have more developed social capital, which can make them more informed or aware of new cancer therapies. The relationship between SES and health literacy is well documented in the literature.22–24 This, in turn, may cause individuals with a higher SES to be more demanding of the latest and advanced therapies as compared to individuals with a low SES. One way to indirectly test this pathway would be to assess if individuals with a high SES switch hospitals of therapy between diagnosis and administration of Durvalumab, and therefore, bypassing the initial plan of treatment provided by the medical group at the time of diagnosis.
Moreover, other studies have stressed the link between race, ethnicity and SES as a factor for receiving appropriate care and treatment, not only in lung cancer25, but also more generally.26 For example, in the Netherlands non-Western ethnic groups have been shown to receive adjuvant chemotherapy for colon cancer less, which cannot be explained away by SES.27 Therefore, a potentially mutual influence between SES and ethnicity and their combined effect on Durvalumab (or other immunotherapies) could be studied in future research.
Finally, in this study, we focused on a patient cohort with good performance status, and thereby excluding patients with likely multimorbidity. Nevertheless, future research could pay more attention to the interplay between multimorbidity and receiving immunotherapy for lung cancer, also given that evidence suggests that multimorbidity and ethic differences are related in the Netherlands and not fully explained away by SES.28 Similar findings have been reported worldwide for a variety of clinical factors.29,30
Conclusion
Our work provides evidence that patients with high SES are more likely to receive the immunotherapy Durvalumab in the Netherlands for the treatment of stage III NSCLC, after chemoradiation. We showed that this effect cannot be explained by differential early mortality of patients by socio-economic status (SES). Our data remain inconclusive regarding the differential effect of driving distance by SES.
Data Availability
Data are not directly available, but need to be requested from Integraal Kankercentrum Nederland (IKNL) through a formal application procedure (www.iknl.nl)
Appendix
A. Bayesian g-formula and model building algorithm
Prior to the Monte Carlo integration, we performed a model comparison procedure to test whether (i) SES modifies the effect of driving distance similarly across the three SES groups, and (ii) whether year of diagnosis modifies the effect of SES, therefore testing if the main effect of SES changed over the years. The comparison was done between the baseline model (Eq. 1 below) and two additional nested models (Eq. 2-3):
where , indicate the coefficients of SES, hospital of first contact, performance status, distance, age, year of diagnosis and natural spline of time-to- Durvalumab with 5 degrees of freedom respectively; Dist, Age, DiagYear, signifies the respective numeric variables. Index-based parametrization of the categorical variables31 and z-score scaling of the numerical variables was followed. In addition, bhospital was modeled as a random intercept, drawn from a normal distribution with estimated hyperparameters μ and σ, both drawn from a standard normal distribution: μ ∼ N(0,1) and σ ∼ N(0,1). The hyperparameters were estimated through a non-centered parametrization. The rest of the parameters had a standard normal prior as well.