Association between benzodiazepine receptor agonist use and increased mortality among patients hospitalized for COVID-19: results from an observational study

2.1. Setting and Cohort Assembly

We conducted a multicenter observational retrospective cohort study at 36 AP-HP hospitals. We included all adults aged 18 years or over who have been admitted to these medical centers from the beginning of the epidemic in France, i.e. January 24^th, until May 1^st. COVID-19 was ascertained by a positive reverse-transcriptase–polymerase-chain-reaction (RT-PCR) test from analysis of nasopharyngeal or oropharyngeal swab specimens. This observational study using routinely collected data received approval from the Institutional Review Board of the AP-HP clinical data warehouse (decision CSE-20-20_COVID19, IRB00011591). AP-HP clinical Data Warehouse initiatives ensure patient information and informed consent regarding the different approved studies through a transparency portal in accordance with European Regulation on data protection and authorization n°1980120 from National Commission for Information Technology and Civil Liberties (CNIL). All procedures related to this work adhered to the ethical standards of the relevant national and institutional committees on human experimentation and with the Helsinki Declaration of 1975, as revised in 2008.

2.2. Data sources

We used data from the AP-HP Health Data Warehouse (‘Entrepôt de Données de Santé (EDS)’). This warehouse contains all available clinical data on all inpatient visits for COVID-19 to any Greater Paris University hospital. The data included patient demographic characteristics, vital signs, laboratory test and RT-PCR test results, medication administration data, medication lists during current and past hospitalizations in AP-HP hospitals, current diagnoses, discharge disposition, and death certificates.

2.3. Variables assessed

We obtained the following data for each patient at the time of the hospitalization through electronic health records (17,18): sex, age, hospital, obesity, self-reported current smoking status, any medical condition associated with increased clinical severity related to COVID-19 or BZRA use (21–24), any current psychiatric disorder, any medication prescribed according to compassionate use or as part of a clinical trial (25), any psychotropic medication, including any antidepressant (26), mood stabilizer (i.e. lithium or antiepileptic medications with mood stabilizing effects), and antipsychotic medication, respiratory depression, defined by a respiratory rate < 12 breaths/min or a resting peripheral capillary oxygen saturation in ambient air < 90% (25), and other clinical markers of disease severity, defined as having temperature > 40°C or systolic blood pressure < 100 mmHg or respiratory rate > 24 breaths/min or plasma lactate levels higher than 2 mmol/L (25). For these two later variables, a third category for missing data was added. These variables are detailed in eMethods.

2.4. Benzodiazepine receptor agonist (BZRA) use

Study baseline was defined as the date of hospital admission. BZRA use was defined as receiving these medications at baseline, i.e. within the first 48 hours of hospital admission, and before the end of the index hospitalization or death, to minimize potential confounding effects of late prescription of BZRAs, which can be used in palliative care. We used this delay because we considered that, in a context of overwhelmed hospital units during the COVID-19 peak incidence, all patients may not have received or been prescribed their usual medication regimens the first day of their hospital admission, or this treatment may not have been recorded at this time.

2.5. Primary endpoint

The primary endpoint was the time from study baseline to death. Patients without an end-point event had their data censored on May 1^st, 2020.

2.6. Statistical analysis

We calculated frequencies of all baseline characteristics described above in patients receiving or not receiving BZRAs and compared them using standardized mean differences (SMD).

To examine the association between BZRA use and the endpoint of death, we performed Cox proportional-hazards regression models (27). To help account for the nonrandomized prescription of BZRAs and reduce the effects of confounding, the primary analysis used propensity score analysis with inverse probability weighting (IPW) (19,28). The individual propensities for receiving any BZRA were estimated using a multivariable logistic regression model that included patient characteristics (i.e., sex, age, hospital, obesity, current smoking status, any significant medical condition, and any current psychiatric disorder) and other other medications, including any medication prescribed according to compassionate use or as part of a clinical trial and other psychotropic medications, comprising any antidepressant, any mood stabilizer, and any antipsychotic medication (see eMethods). In the inverse-probability-weighted analyses, the predicted probabilities from the propensity-score models were used to calculate the stabilized inverse-probability-weighting weights (28). The association between BZRA use and the endpoint was then estimated using an IPW Cox regression model. In case of non-balanced covariates, an IPW multivariable Cox regression model adjusting for the non-balanced covariates was also performed. Kaplan-Meier curves were performed using the inverse-probability-weighting weights (29,30), and their pointwise 95% confidence intervals were estimated using the nonparametric bootstrap method (30).

We conducted five sensitivity analyses. First, we performed a multivariable Cox regression model including as covariates the same variables as in the IPW analysis. Second, we used a univariate Cox regression model in a matched analytic sample using a 1:1 ratio, based on the same variables used for the IPW analysis and the multivariable Cox regression analysis. To reduce the effects of confounding, optimal matching was used in order to obtain the smallest average absolute distance across all clinical characteristics between exposed patients and non-exposed matched controls. Third, to examine a potential indication bias of prescription of BZRAs in ICUs as a possible treatment for palliative care or as an aid to oral intubation, we reproduced the main analyses after excluding all patients who had been hospitalized in ICUs. Fourth, we examined whether our findings were similar in models imputing missing data using multiple imputation (31) instead of excluding patients with any missing data as done in the main analyses. Finally, to examine the potential influence on the results of excluding patients who received BZRA after 48 hours from hospital admission, we reproduced the main analysis while including all participants who received any BZRA at any time from hospital admission until the end of the index hospitalization or death, and considering BZRA use as a time-dependent variable (27). In that analysis, patients who received any BZRA after hospital admission enter the analysis risk-sets at the time of actual first initiation of the treatment. In all these analyses, weighted Cox regression model was used when proportional hazards assumption was not met (32).

We performed several additional exploratory analyses. First, we searched for a potential dose-dependent relationship by testing the association between daily dose of BZRA received at baseline (converted into diazepam-equivalent dose (33) and dichotomized at the mean value) and the endpoint, adjusted for the same covariates as used in the main analysis, among patients who received a BZRA. Second, we examined whether respiratory depression or other clinical markers of disease severity may at least partly explain this association by adjusting successively for these variables. Finally, we examined the relationships between each BZRA and mortality using the same statistical approach as described for the main analysis.

For all associations, we performed residual analyses to assess the fit of the data, check assumptions, including proportional hazards assumption using proportional hazards tests and diagnostics based on weighted residuals (27,34), and examined the potential influence of outliers. To improve the quality of result reporting, we followed the recommendations of The Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) Initiative. Because our main hypothesis focused on the association between BZRA use and mortality, statistical significance was fixed a priori at two-sided p-value <0.05. Only if a significant association were found in the main analysis, we planned to perform additional exploratory analyses as described above. All analyses were conducted in R software version 2.4.3 (R Project for Statistical Computing).

3.1. Characteristics of the cohort

Of the 17,131 patients with a positive COVID-19 RT-PCR test who had been hospitalized for COVID-19, 1,963 patients (11.5%) were excluded because of missing data or their young age (i.e. less than 18 years old of age). In addition, of the 1,473 patients who received a BZRA at any time during the visit, 787 were excluded because they were prescribed it more than 48 hours after hospital admission. Of the remaining 14,381 adult patients, 686 (4.8%) received a BZRA in the first 48 hours of hospitalization at a mean diazepam-equivalent dose of 19.7 mg per day (SD=25.4, median=10.0, lower quartile=5.0, higher quartile=22.9, range=1.85-240.0 mg) (eFigure 1). Of these 686 patients, 41.1% (N=282) had a prescription of BZRA during a prior hospitalization at AP-HP in the past 6 months, and 36.4% (N=183) received at least two BZRA medications. The median delay from hospital admission to first prescription of BZRA was 0.94 days (SD=0.55).

PT-PCR test results were obtained after a median delay of 1.2 days (SD=12.7) from hospital admission date. This median delay was of 1.0 day (SD=11.6) in the exposed group and of 1.2 day (SD=12.8) in the non-exposed group.

Over a mean follow-up of 14.5 days (SD=18.1; median=7 days), 1,320 patients (9.2%) had an end-point event at the time of data cutoff on May 1^st. Among patients who received a BZRA, the mean follow-up was 12.8 days (SD=13.1, median=8 days), while it was of 14.5 days (SD=18.3, median=7 days) in those who did not.

All baseline characteristics were independently and significantly associated with the endpoint. A multivariable Cox regression model showed that age, sex, hospital, any medical condition, any current psychiatric disorder, respiratory depression and any other clinical marker of severity were significantly and independently associated with the endpoint (eTable 1).

The distributions of patient characteristics according to BZRA use are shown in Table

1. In the full sample, BZRA use significantly differed according to all characteristics except for sex ratio, and the direction of the associations indicated older age and greater medical severity of patients receiving any BZRA. After applying the propensity score weights, these differences were substantially reduced. In the matched analytic sample, there were no significant differences in any characteristic, except the greater prevalence of men in the non-exposed group (Table 1).

3.2. Study endpoint

The endpoint of death occurred in 186 patients (27.1%) who received a BZRA at baseline and 1,134 patients (8.3%) who did not. The crude, unadjusted analysis (hazard ratio (HR=3.20; 95% CI=2.74-3.74; p<0.001), the primary analysis with inverse probability weighting (HR=1.61; 95% CI=1.31-1.98; p<0.001) and the multivariable IPW Cox regression adjusting for unbalanced covariates (HR=1.56; 95% CI=1.29-1.89; p<0.001) showed a significant association between BZRA use and increased mortality (Figure 1; Table 2). A post-hoc analysis indicated that we had 80% power in the crude analysis to detect a hazard ratio of at least 1.27.

• Download figure
• Open in new tab

Figure 1. Kaplan-Meier curves for mortality in the full sample crude analysis (N=14,381) (A), in the full sample analysis with IPW (N=14,381) (B), and in the matched analytic sample (N=1,372) (C), according to BZRA use at baseline.

The shaded areas represent pointwise 95% confidence intervals.

In sensitivity analyses, the multivariable Cox regression model in the full sample also a significant association (HR=1.94; 95% CI=1.45-2.59; p<0.001), as did the univariate Cox regression model in a matched analytic sample using a 1:1 ratio (HR=1.34; 95% CI=1.08-1.67; p=0.009) (Table 2). Similarly, the primary analysis using imputed data yielded significant results (eTable 2), as did that considering BZRA use as a time-dependent variable and including all patients who received a BZRA at any time during the visit (eTable 3). The exclusion from the analyses of the patients who had been admitted to ICUs did not alter the significance of the association (eTable 4).

Additional analyses showed a significant dose-dependent relationship between baseline daily BZRA dose and the endpoint (HR=1.55; 95% CI=1.08-2.22; p=0.017), based on the primary IPW analysis. Additional adjustments for respiratory depression, any other clinical markers of disease severity, or both, resulted in still significant associations, which were of similar magnitude to that observed in the primary analysis (eTable 5). We found that all individual BZRAs were significantly associated with increased mortality, except for diazepam (eTable 6). Following adjustments, there were significant associations of any BZRA other than diazepam, midazolam, and any BZRA other than diazepam and midazolam with increased mortality, as compared to not receiving BZRAs (eTable 6). Finally, compared with any other BZRA treatment, diazepam use was significantly associated with reduced mortality in the crude analysis (HR, 0.44; 95% CI, 0.23 to 0.86; p=0.017), in the primary analysis (HR, 0.31; 95% CI, 0.13 to 0.74; p=0.008), and in the sensitivity analyses (eFigure 2; eTables 6 to8).