Outcomes in Patients with Acute Hypoxemic Respiratory Failure Secondary to COVID-19 Treated with Noninvasive Respiratory Support versus Invasive Mechanical Ventilation

Data Preprocessing

Data were preprocessed to exclude participants under age 18, visits for which a phenotype was not available, and participants who exited the hospital (either alive or dead) before treatment. Some patients appeared to have multiple applicable hospital visits. Hospital transfers were identified by looking for multiple entries with the same ID where the patients had the same gender, race, ethnicity, age (max one year difference), the patient was alive at transfer, and consecutive hospital discharge and admission times were fewer than 24 hours apart. Identified transfers were combined into a single visit entry, and for six visits, phenotypes were adjusted as appropriate. Overall, 48 transfers were identified, and 122 later visits were discarded. Later visits were discarded in order to sample from the population of first COVID-19 visits. Height and weight were defined preferentially as the first recorded height and weight and secondarily as the admission height and weight (available for ICU data only); among ICU visits, there were 72 instances of using an admission instead of a time-stamped height and 11 instances of using an admission instead of a time-stamped weight. Heights below 48 inches and weights below 50 pounds were classified as potential data errors, and the next valid height and/or weight was sought. For nine people, this resulted in no height data being available; for two people these rules resulted in no weight data being available. In the ICU data, 72 visits had an admission height but not a time-stamped height, and 11 visits had an admission weight but not a time-stamped weight. In order to work with a deidentified data set, for the 42 participants with age over 89 years (18 participants with “>89”, and 24 with specific ages >89), their age was marked as 90. Race and ethnicity were standardized. For the respiratory rate, FiO₂, oxygen flow rate, and SpO₂ variables, the value closest to and before first treatment were identified. For people with phenotype 8 (evidence of all three treatments but unclear order), the first treatment start time was taken to be the earliest treatment start time identified by the phenotyping algorithm. Oxygen flow rates were used to calculate FiO₂ (FiO₂ = 100(0.21 + oxygen flow [L/min^-1] x 0.03), and the FiO₂ value (either raw or calculated) selected for later use was the latter one.

Propensity Score Estimation via Generalized Boosted Models

We follow McCaffrey et al.’s (2013) approach to estimating propensity scores by using generalized boosted models (GBMs).¹ GBMs estimate propensity scores for the three treatments by fitting the models to the dichotomized outcomes of each treatment versus the other two. This results in three GBMs being fit each time propensity scores are estimated for models comparing NIPPV, HFNO, and IMV. Because each GBM consists of an iterative process of combining simple regression trees into a piecewise constant function, for each GBM, we must choose a point at which to stop adding onto the regression tree. Four criteria can readily be used to choose a stopping point; all four attempt to describe the difference between the distributions of covariates for the comparison groups for each dichotomized outcome (e.g., between patients who received IMV and all patients); the goal is to minimize the differences between the groups. Note that we wish to estimate the average treatment effect between pairs of treatments within the context of our entire population. This is considered in each stopping criterion:¹

1. Mean Absolute Standardized Bias: mean across covariates of the absolute standardized difference between the weighted mean of the covariate for the treatment group of interest and the unweighted mean of the covariate for the entire sample. Weights are the inverse of the estimated propensity scores. Because we are interested in estimating the treatment effect in the general population of people receiving any of the three treatments, we standardize (i.e., divide) here by the estimated standard deviation from the unweighted, pooled, full sample.

2. Maximum Absolute Standardized Bias: maximum across covariates of the absolute standardized difference between the weighted mean of the covariate for the treatment group of interest and the unweighted mean of the covariate for the entire sample. Weights and standardization are as described above.

3. Mean Kolmogorov-Smirnov mean across covariates of the supremum of pointwise absolute differences between the weighted empirical distribution function for the covariate among those treated with the treatment of interest and the unweighted empirical distribution function for the entire sample.

4. Maximum Kolmogorov-Smirnov: maximum across covariates of the supremum of pointwise absolute differences between the weighted empirical distribution function for the covariate among those treated with the treatment of interest and the unweighted empirical distribution function for the entire sample.

For each GBM of a dichotomized outcome, we allow a maximum of 5000 simple regression trees. Empirically, this number appeared to allow all four criteria above to achieve their minimum and thus was large enough without being excessively large. For each outcome and each stopping criterion, the optimal number of trees was selected. Thus, for each treatment, we estimated four propensity scores for having been given the treatment in question. Inverse probability of treatment weighting was used in subsequent analyses by taking the inverse of the propensity scores.

For all non-interaction and non-sensitivity models, a separate analysis was done for each stopping criterion to ensure that results were not dependent on the particular criterion used. No result was dependent on stopping criterion.

Sensitivity Analyses for Time to Intubation

The competing risks analyses of time-to-intubation made a series of analytical decisions that could impact the subsequent inference. We thus designed a series of sensitivity analyses that systematically varied our assumptions and analytical approach to ensure that results were not dependent on specific decisions.

The primary reason for sensitivity analyses was that 326 visits could not be reliably classified using the phenotyping algorithm. These unclassifiable visits always contained evidence of all three treatments (IMV, HFNO, and NIPPV). However, in these cases, conflicting information in the medical record often made it unclear to the algorithm which treatment came first. For example, a patient might appear to be intubated, switched to NIPPV, then reintubated in a 30-minute period. Such information would cause the algorithm to fail to assign a treatment ordering to that visit. While these same cases were present in the other analyses, in those analyses, we were able to impute treatment, estimate propensity scores for each imputed data set, model each data set, then combine model estimates across data sets using Rubin’s Rules.

The weighted Gray’s test of competing risks, however, used a bootstrapped null distribution. There was thus no obvious way to use Rubin’s Rules. A common alternative approach to combining propensity scores and multiple imputation --- estimating propensity scores for each imputed data set, averaging them, and choosing a single imputed data set at random to use in the analysis along with the averaged propensity scores --- wasn’t applicable because of the need to impute treatment. When treatment is imputed, it is unclear how well a particular data set approximates the true but unknown treatment identities, and thus using the average propensity score approach is inappropriate.² Moreover, different observations could be included in each data set since people imputed to have an invasive first treatment would be excluded. Determining how to compute average propensity scores for people who were not consistently included in the analyses and how to appropriately handle the data from those individuals left out of the randomly sampled data set is beyond the scope of this paper. A final challenge of using the average propensity score approach came from the need to impute the time to event for 52 people; if the imputed outcomes were particularly outlying and influential for a randomly selected data set, then the inference may not be correct. An alternative to Rubin’s Rules — judging statistical significance via the median p-value from the pooled p-values across all the imputed data sets, while appealing, is less common and would need to be utilized in an analysis scenario not covered in Eekhout et al.³

The sensitivity analyses thus encompass a range of ways to handle the above challenges. Each combines a data decision, an outcome decision, an imputation decision, and an inference approach decision to arrive at an analysis pipeline. Note that some combinations of the four decisions were not possible. Moreover, where it was possible to take either the average propensity score or median p-value approach, we did not always conduct both corresponding analyses. Rather, we defaulted to the average propensity score approach except when we are imputing time-to-event information; in those cases, we also conduct the corresponding analysis with the median p-value inference approach as an extra confirmation that the specific data set chosen in the average propensity score approach didn’t have particularly influential imputed times-to-event.

Data Decision: How and how much of the unclassifiable data to include.

- Option A: We assume that we know everything about the unclassifiable data and base our knowledge on treatment start information. We exclude patients with IMV as the first marked treatment and assign patients with HFNO or NIPPV marked as the first treatment to the corresponding treatment group. If the two noninvasive therapies tie for first treatment (7 people), we randomly select one.

- Option B: We assume that we know whether the first treatment was IMV or NIRS based on the treatment time stamps. Patients with IMV as the first treatment are excluded from analysis. However, because of potential nomenclature inconsistencies between NIPPV and HFNO, we assume that the first noninvasive treatment could have been either NIPPV or HFNO. We impute the initial noninvasive treatment type.

- Option C: We assume that we do not know what treatment came first and thus whether or not anyone with the unclassifiable phenotype should be included in this competing risks analysis. Thus, we exclude them all from the analysis.

Unknown Outcome Decision:

• - Yes: Include patients with unknown time-to-event.

• - No: Do not include patients with unknown time-to-event.

Imputation Decision: The decision whether to do imputation impacts which predictors are used in the propensity score calculations. If imputation is conducted, the propensity score algorithm uses the previously described predictors; otherwise, it uses a restricted subset for which all information is available (age, gender, hospital size).

• - Yes: Conduct multiple imputation by chained equations.

• - No: Do not impute missing data.

Inference Approach Decision:

• - Average Propensity Score Approach: In this approach, the propensity scores for each visit are averaged across imputed data sets. These averaged scores are used with a randomly-selected imputation dataset. This approach is not appropriate when treatment is imputed.

• - Median P-Value Approach: In this approach, the weighted Gray’s test is conducted on each imputed data set (with its own set of propensity scores). Then, the median p-value across data sets is taken for inference.

• - Single Propensity Score Approach: This case arises when we are not imputing data. Then, we estimate propensity scores once (based on the limited predictor set) and use them in the weighted Gray’s test.

The ten proposed competing risks analyses correspond to the following combinations of the four decision points (all p-values were <0.001):

We considered the analysis Data Decision B + Unknown Outcome Decision: No + Imputation Decision: Yes + Inference Approach Decision: Median P-Value to be primary and the others to be sensitivity analyses. Additionally, given the lack of inferential difference associated with the different propensity score estimation stopping rules, for all sensitivity analyses, we only used the mean absolute standardized bias stopping rule. We did, however, use all four stopping rules for the primary competing risks analysis.

• Download figure
• Open in new tab

Figure 1: Time-to-intubation sensitivity analyses.

In the estimation of these weighted survival plots, people who die remain in the risk set after their death because they have neither been censored nor experienced intubation. All time-to-intubation analyses showed statistically significant differences between initial treatment with HFNO and NIPPV (p < 0.001).

• Download figure
• Open in new tab

Figure 2: Time-to-In-Hospital Death (NIRS vs IMV).

Time-to-in-hospital death weighted Kaplan-Meier curves for noninvasive respiratory support versus invasive mechanical ventilation. Curves use propensity scores estimated via the mean absolute standardized bias stopping rule in the generalized boosted models.

• Download figure
• Open in new tab

Figure 3: Time-to-Hospital Discharge Alive (NIRS vs IMV).

Time-to-hospital discharge alive weighted Kaplan-Meier curves for noninvasive respiratory support versus invasive mechanical ventilation. Time zero is time of first treatment assignment. In these Kaplan-Meier estimates, people experiencing death are treated as censored. Curves use propensity scores estimated via the mean absolute standardized bias stopping rule in the generalized boosted models.

• Download figure
• Open in new tab

Figure 4: Time-to-Hospital Discharge Alive (NIRS vs IMV) Sensitivity Analysis.

Time-to-hospital discharge alive weighted Kaplan-Meier curves for noninvasive respiratory support versus invasive mechanical ventilation. Time zero is time of hospital admission. In these Kaplan-Meier estimates, people experiencing death are treated as censored. Curves use propensity scores estimated via the mean absolute standardized bias stopping rule in the generalized boosted models.

• Download figure
• Open in new tab

Figure 5: Time-to-ICU Discharge Alive Secondary Outcome (NIRS vs IMV).

Weighted Kaplan-Meier curves for noninvasive respiratory support versus invasive mechanical ventilation. Time zero is time of ICU admission. In these Kaplan-Meier estimates, people experiencing death are treated as censored. Curves use propensity scores estimated via the mean absolute standardized bias stopping rule in the generalized boosted models.

• Download figure
• Open in new tab

Figure 6: Time-to-In-Hospital Death (HFNO vs NIPPV vs IMV).

Time-to-in-hospital death weighted Kaplan-Meier curves for HFNO versus NIPPV versus invasive mechanical ventilation. Curves use propensity scores estimated via the mean absolute standardized bias stopping rule in the generalized boosted models.

• Download figure
• Open in new tab

Figure 7: Time-to-Hospital Discharge Alive (HFNO vs NIPPV vs IMV).

Time-to-hospital discharge alive weighted Kaplan-Meier curves for HFNO versus NIPPV versus invasive mechanical ventilation. Time zero is time of first treatment assignment. In these Kaplan-Meier estimates, people experiencing death are treated as censored. Curves use propensity scores estimated via the mean absolute standardized bias stopping rule in the generalized boosted models.

• Download figure
• Open in new tab

Figure 8: Time-to-Hospital Discharge Alive (HFNO vs NIPPV vs IMV) Sensitivity Analysis.

Time-to-hospital discharge alive weighted Kaplan-Meier curves for HFNO versus NIPPV versus invasive mechanical ventilation. Time zero is time of hospital admission. In these Kaplan-Meier estimates, people experiencing death are treated as censored. Curves use propensity scores estimated via the mean absolute standardized bias stopping rule in the generalized boosted models.

• Download figure
• Open in new tab

Figure 9: Time-to-ICU Discharge Alive (HFNO vs NIPPV vs IMV).

Weighted Kaplan-Meier curves for HFNO versus NIPPV versus invasive mechanical ventilation. Time zero is time of ICU admission. In these Kaplan-Meier estimates, people experiencing death are treated as censored. Curves use propensity scores estimated via the mean absolute standardized bias stopping rule in the generalized boosted models.

STROBE Checklist