Emulated trial investigating effects of multiple treatments: estimating combined effects of mucoactive nebulisers in cystic fibrosis using registry data

Introduction People with cystic fibrosis (CF) are often on multiple long-term treatments, including mucoactive nebulisers. In the UK, the most common mucoactive nebuliser is dornase alfa (DNase). A common therapeutic approach for people already on DNase is to add hypertonic saline (HS). The effects of DNase and HS used alone have been studied in randomised trials, but their effects in combination have not. This study investigates whether, for people already prescribed DNase, adding HS has additional benefit for lung function or use of intravenous antibiotics. Methods Using UK CF Registry data from 2007 to 2018, we emulated a target trial. We included people aged 6 years and over who were prescribed DNase without HS for 2 years. We investigated the effects of combinations of DNase and HS over 5 years of follow-up. Inverse-probability-of-treatment weighting was used to control confounding. The period predated triple combination CF transmembrane conductance regulator modulators in routine care. Results 4498 individuals were included. At baseline, average age and forced expiratory volume in 1 s (FEV1%) predicted were 21.1 years and 69.7 respectively. During first year of follow-up, 3799 individuals were prescribed DNase alone; 426 added HS; 57 switched to HS alone and 216 were prescribed neither. We found no evidence that adding HS improved FEV1% at 1–5 years, or use of intravenous antibiotics at 1–4 years, compared with DNase alone. Conclusion For individuals with CF prescribed DNase, we found no evidence that adding HS had an effect on FEV1% or prescription of intravenous antibiotics. Our study illustrates the emulated target trial approach using CF Registry data.


Additional notes on methodology 1.1 Causal estimands
The treatments of interest are dornase alfa and hypertonic saline.In the UK CF Registry, treatment use at each annual review is recorded as a yes/no indicating whether treatment was prescribed over the past year.Time is measured in years since baseline.Let  , and  , denote whether dornase alfa and hypertonic saline, respectively, was recorded for the ith person at time k (k=1,2,3,4,5).Let  , denote which treatment combination the ith person was on at time k (i.e. between times k-1 and k).Then  , is defined as:  , = { 0   , = 0   , = 0 1   , = 0   , = 1 2   , = 1   , = 0 3   , = 1   , = 1 Henceforth we suppress the subscript .The outcome at time  is denoted   .The outcomes were FEV1% (continuous, measured on the day of the annual review visit) and IV antibiotic use (binary, denoting whether or not any IV antibiotics were prescribed since the last annual review visit).Recall that at baseline (time 1) all individuals had been using DNase for 2 years, according to our inclusion criteria.Let  ̅  = { 1 , … ,   } denote the treatment history from time 1 to time point  and let    ̅  denote the potential outcome that would be observed for an individual with a particular treatment history  ̅  .Our primary aim was to compare the strategies of adding hypertonic saline to DNase (denoted DN&HS) up to follow-up time of interest and continuing DNase alone up to the follow-up time of FEV1% interest (denoted DN).Using our above notation, for the continuous outcome of the main estimands of interest are defined as: 3 year: 4 year: 5 year: ( 5  ̅ 5 =(3,3,3,3,3) ) − ( 5  5 =(2,2,2,2,2) ) Comparisons between other treatment strategies were of secondary interest.We also compared the treatment strategies of switching from DNase to hypertonic saline and continuing HS alone to the follow-up time of interest (denoted HS) versus continuing DNase alone (denoted DN): We also compared the treatment strategies of dropping DNase and not adding hypertonic saline (denoted Nil) versus continuing DNase alone (DN): For the binary outcome of whether a person was prescribed any days of IV antibiotic treatment the estimands are odds ratios instead of mean differences.These are discussed in more detail in section 2.1.

Confounding variables
To obtain unbiased estimates of the treatment effects, we needed control for both time-invariant and time-varying confounders.Figures S.1 and S.2 are the directed acyclic graphs (DAGs) which show the assumed relationships between variables in our data for the analyses with FEV1% and IV days as the outcome, respectively.Both DAGs are simplified versions of reality.We have not included long-BMJ Publishing Group Limited (BMJ) disclaims all liability and responsibility arising from any reliance Supplemental material placed on this supplemental material which has been supplied by the author(s)    As can be seen from Figures S.1 and S.2, the variables included as time-invariant confounders were: age at baseline, genotype, sex, ethnicity, rate of decline in FEV1%, BMI z-score at baseline and FEV1% at baseline.The variables included as time-varying confounders were: pancreatic insufficiency, ivacaftor use, p. Aeruginosa infection, staphylococcus aureus infection, nontuberculous mycobacteria infection, hospital admissions for intravenous antibiotics, days on intravenous antibiotics, past BMI z-score and past FEV1%.
Genotype was classed as either high risk, low risk or not assigned as previously defined 1 .Ethnicity was classed as white or non-white due to small numbers in non-white ethnic groups in this population.Rate of decline in FEV1% represented the change in FEV1% observed prior to baseline.We defined the following linear mixed model with random slope and intercept: BMJ Publishing Group Limited (BMJ) disclaims all liability and responsibility arising from any reliance Supplemental material placed on this supplemental material which has been supplied by the author(s)

Granger E
Where  ∈ {0,1,2,3,4} is the number of years before baseline ( = 0 is the baseline year).The estimate of the slope parameter ( 1 +  1 ) for each individual is used as a time-invariant variable representing rate of change in FEV1%.
Pancreatic insufficiency was a yes/no indicator where individuals were assigned "yes" if they were prescribed pancreatic enzyme supplements.IV hospital admissions was a yes/no indicator were yes indicated individuals had at least one hospital admission for IV antibiotics over the past year.IV days included home and hospital admissions and was categorised as: 0, 1-4, 15-28 and 29+.BMI z-scores were calculated using the WHO reference distribution 2 and FEV1% was calculated using the Global Lung Initiative equations 3 .

Inverse-probability-of-treatment weighted estimation of marginal structural models.
Let   denote the set of time-invariant confounders and   denote the set of time-varying confounders recorded at time k.In both the FEV1% and IV days analyses,   = { 0 , , , ℎ,     1%, 1% 0  0 } For the FEV1% analysis,   is defined as: For the IV days analysis,   is defined as: Then, the stabilised inverse-probability-of-treatment weights for individual  at time  (.  ) were defined as: Weights were also used to account for missing data in FEV1% or BMI (.  ), loss-to-follow-up (.  ), censoring due to death or transplant (.  ) and time-varying eligibility due to use of lumacaftor/ivacaftor or texacaftor/ivacaftor (.  ) or mannitol (.  ).
The probabilities required for each set of weights were obtained using logistic regression.For .  , and .  , the outcomes were indicators for use of the relevant treatments.
For each individual, we excluded time-points with missing data for FEV1% or BMI.To account for the missing data, the remaining individuals were re-weighted by the inverse of their probability of remaining in the study at a given time.The weights, .  , were defined using a similar equation as the one for .  , but the outcome was an indicator for missingness in FEV1% or BMI for the i th individual at time .
Individuals who were lost to follow-up, died or had an organ transplant were censored at the time of the event.For the loss to follow-up weights, the outcome at time k was an indicator for whether the individual was lost to follow-up at time  + 1.For the censoring weights due to death or organ transplant (whichever occurred first), the outcome at time k was an indicator for whether the individual died or had a transplant between times  and  + 1.
All weights were stabilised and probabilities were conditioned on the same variables as the probabilities defined in the inverse-probability-of-treatment weights.
The combined weight for individual  at time point k (.  ) was defined as a product of all of the above weights:

Granger E
For our main analysis, we specified the following linear marginal structural model (MSM) for the continuous outcome of FEV1%: The parameters of the MSM are estimated by fitting the model using the observed data weighted using the combined weight.This enables estimation of the estimands specified in supplementary section 1.1.We note that these are marginal mean differences, as the conditional and marginal mean differences coincide for the linear MSM above. .
For the binary outcome of whether the individual was prescribed any IV antibiotics over the past year, the marginal structural model (MSM) used for the main analysis was: This can be fitted using the observed data weighted using the combined weight.This results in estimates of conditional odds ratios.For example, our primary odds ratios of interest are For the analyses investigating whether the treatment effects differed by FEV1% measured at baseline, the above MSMs were extended to include an interaction between FEV1% (a component of   ) and (  = ).

Missing data
We found 5360 individuals with CF who were documented as having been prescribed dornase alfa and not hypertonic saline for at least two consecutive years between 2007 and 2017, and who had at least one baseline visit and one follow-up year.After excluding individuals who were under the age of 6 years, had received a solid organ transplant by baseline, or were prescribed mannitol, tezacaftor/ivacaftor or lumacaftor/ivacaftor at baseline, we were left with 4810 individuals who were eligible for inclusion.Table S.1 shows the amount of missing data by year for those individuals.Note that this includes people who were transplanted, or prescribed mannitol, tezacaftor/ivacaftor or lumacaftor/ivacaftor post-baseline.We excluded individuals with missing data on time-invariant variables (genotype and ethnicity) and individuals with missing FEV1% data at baseline ( = 0).
The last observation carried forward was used to impute the missing infection data (P.aeruginosa infection, Staphylococcus aureus infection, Bukholderia cepa infection and NTM).This was considered a valid approach as these infections are usually long-term and there was no missing data in these variables at time 0.
Missing data weights were used to account for missing BMI and FEV1% (except for individuals who had missing FEV1% at baseline, who were excluded).

Summary of exclusions due to loss to follow-up, death, transplant, ineligibility and missing data.
Table S.2 gives the numbers of individuals who were excluded or censored each year for different reasons.Individuals who were censored due to loss-to-follow-up, death or transplant, and these individuals account for the decreasing number observed by follow-up year.For example, in between visits 1 and 2, 250 individuals were lost-to-follow-up, 58 died and 29 received an organ transplant.By visit 2, 4498-(250+58+29)=4161 individuals remained in the study.
Individuals who were temporarily excluded due to missing data or temporary ineligibility (due to initiating treatment with CFTR modulators or mannitol) were allowed to re-enter the study, and these numbers account for the differences between the number of people observed in each follow up year and the number of people included in the final analysis (final N).
Column percentages are given with respect to the sample sizes in row 1.

Summary of the numbers of people prescribed each treatment combination and flow of participants between treatment combinations by year
Figure S.3 shows the number of people prescribed each treatment combination by year.Across all follow-up years, the percent of individuals using neither DNase nor hypertonic saline ranged between 4.2% and 5.2%.The percentages prescribed DNase only and hypertonic saline only ranged from 51.7%-81.2%and 1.2%-2.8%respectively.The percentage prescribed both DNase and hypertonic saline ranged from 9.1%-29.6%.BMJ Publishing Group Limited (BMJ) disclaims all liability and responsibility arising from any reliance Supplemental material placed on this supplemental material which has been supplied by the author(s)     BMJ Publishing Group Limited (BMJ) disclaims all liability and responsibility arising from any reliance Supplemental material placed on this supplemental material which has been supplied by the author(s)

Distribution of weights
term arrows (e.g. from a variable recorded at time  − 2 to one recorded at time ), or relationships between the time-dependent variables measured at a given visit.FEV1% and BMI are recorded at the annual review visit.The following covariates (included together in a box in the DAGs) at a given annual visit refer to whether infection, pancreatic insufficiency, IV hospitalisation or ivacaftor prescription were recorded since the previous annual visit: pancreatic insufficiency, Pseudomonas aeruginosa infection, Staphylococcus aureus infection, Nontuberculous Mycobacteria infection, IV hospitalisation and Ivacaftor use.

Figure S. 1 :
Figure S. 1: Directed Acyclic Graph depicting the assumed short-term confounding paths of the treatment-outcome association when FEV1% was the outcome (  denotes FEV1% at time ).

Figure S. 2 :
Figure S. 2: Directed Acyclic Graph depicting the assumed short-term confounding paths of the treatment-outcome association when IV days was the outcome (  denotes binary indicator for IV prescription recorded at time ).

Figure S. 4
Figure S.4 describes the flow of participants between treatment combinations by year.Of the 143 individuals who were using neither DNase nor hypertonic salin in the first year and had 5 years of followup, 31 (21.7%)continued to use neither treatment for 5 years.Of the 51 individuals using hypertonic saline only in the first year (i.e, who switched from DNase to hypertonic saline) and had 5 years of follow-up, 16 (31.4%)remained on hypertonic saline only for 5 years.Of the 2521 individuals who continued to be prescribed DNase only in the first year and had 5 years of follow-up, 1615 (64.1%) remained on DNase only for 5 years.Of the 260 individuals who added hypertonic saline to DNase in the first year and had 5 years of follow-up, 185 (71.2%) remained on this combination for 5 years.

Figure S. 3 :
Figure S.3: Flowchart showing the number of participants in the study and number of participants prescribed each treatment combination by year.

Figure
Figure S. 5: Average FEV1% and proportion of people on IV antibiotics in the whole cohort, by followup visit.Note that the vertical axes are truncated and the changes over time are small.

Figures
Figures S.6 and S.7 show the distribution of inverse-probability-of-treatment (IPT) weights and combined weights by year, respectively (weights are defined in Section 1.3).Boxplots show that the median weights are approximately 1 for each year, as expected.The variance of weights tends to increase by year but there are no extreme values.
.  = .  × .  × .  × .  × .  × .  BMJ Publishing Group Limited (BMJ) disclaims all liability and responsibility arising from any reliance Supplemental material placed on this supplemental material which has been supplied by the author(s)