Data source
We used the National Health Insurance Service-National Health Insurance Database (NHIS-NHID) of South Korea [15], which contains health insurance claims data for the entire Korean population of 50 million inhabitants collected between 1 January 2002 and 31 December 2017. The NHIS-NHID is comprised of anonymized patient identifier and corresponding information on associated sociodemographic characteristics (age, sex, health insurance type, and income level), healthcare utilization history, diagnostic codes based on the International Classification of Diseases 10th Revision (ICD-10), and drug prescription information (national drug codes [NDC], days’ supply, dosage, and administration route). The NDCs are based on the drug’s active ingredient and mapped to the World Health Organization’s Anatomical Therapeutic Chemical classification codes. Moreover, the date of death was linked to the national vital statistics maintained by Statistics Korea.
Study population
We identified patients with either a primary or secondary recorded diagnosis of incident depression (ICD-10: F32, F33, F34.1) from both inpatient and outpatient settings between 1 January 2004 and 31 December 2017 (Fig. 1). Cohort entry was defined as the date of incident diagnosis with depression. The following patients were excluded from our study: (1) those diagnosed with depression (ICD-10: F32, F33), bipolar disorders (F31), manic episodes (F30), or persistent mood disorders (F34) within 2 years prior to cohort entry, to restrict the analysis to incident patients with depression; (2) those prescribed any ADs or BZDs within 2 years prior to cohort entry to restrict to new users of ADs or BZDs; (3) those aged < 18 years at cohort entry as depression is unlikely to be prevalent and because AD or BZD use is not advised in this age group; and (4) missing demographic information.
Exposure definition
We used prescription records of ADs and BZDs from both inpatient and outpatient settings to ascertain exposure, where an intention-to-treat approach was used to define follow-up (Additional file 1: Table S1). Prescriptions of ADs or BZDs within 6 months after cohort entry were eligible. Exposure was classified into two groups and their corresponding index dates were defined as follows: (1) AD monotherapy, date when one class of AD was prescribed in a single prescription; (2) AD+BZD therapy, date when ADs were prescribed together with BZDs in a single prescription on the same day (Additional file 2: Fig. S1 and Additional file 3: Fig. S2).
Outcome definition
Time to all-cause mortality was our primary outcome of interest. Our secondary outcomes were the time to incident diagnosis with suicide attempt/self-harm (from both inpatient and outpatient settings) and time to all-cause hospitalization, defined as a visit to a medical institution resulting in admission. Follow-up began on the index date and ended on the earliest of the outcome occurrence, or end of the study period (31 December 2017).
Potential confounders
Age, sex, health insurance type, residential district, and income level were assessed at cohort entry. Furthermore, comorbidities (anxiety, cancer, cerebrovascular disease, chronic kidney disease, chronic obstructive pulmonary disease, dementia, diabetes mellitus, epilepsy, fractures, hypertension, hyperlipidemia, insomnia, ischemic heart disease, osteoarthritis, Parkinson’s disease, rheumatoid arthritis, and substance abuse) and history of medication use (angiotensin-converting-enzyme inhibitors, angiotensin II receptor blockers, anticholinergics, antiplatelets and anticoagulants, antipsychotics, anticonvulsants, digoxin, non-insulin glucose-lowering agents, anti-inflammatory analgesics, beta-blockers, calcium channel blockers, insulin, lipid-lowering agents, narcotic analgesics, nonsteroidal anti-inflammatory drugs, other anxiolytics, and thiazide diuretics) were assessed within the year before cohort entry. Comorbidities were defined using ICD-10 codes and use of medications was defined using NDCs (Additional file 1: Table S1). The Charlson comorbidity index (CCI) score was also estimated to determine the overall comorbidity burden [16].
To obtain comparability between treatment groups, we conducted propensity score (PS) matching, where the PS of receiving AD therapy was estimated using multivariable logistic regression. Upon conducting chi-square tests to determine statistical significance between each confounder and all-cause mortality, only confounders that had p value < 0.2 were included as independent variables into the multivariable logistic regression model (history of other anxiolytic use was not included in the model) [17]; age and sex were always included regardless to their p value. Matching based on PS was done in a 1:1 ratio for the two groups using the Greedy 5→1 digit matching macro [18, 19], where the c statistic (0.6–0.8) was used to assess model discrimination [20].
Statistical analyses
The baseline characteristics of our study subjects are represented as counts (proportions) and means (standard deviations) for categorical and continuous variables, respectively. We estimated the absolute standardized difference (aSD) to determine imbalances between groups (aSD > 0.1 indicated an important imbalance) [21].
We first aimed to examine the risk of study outcomes by calculating the incidence of study outcomes per 1000 person-years and estimating the absolute risk and risk differences with 95% confidence intervals (CIs). Adjusted hazard ratios (HRs) with 95% CIs for study outcomes associated with AD+BZD therapy versus AD monotherapy were estimated using multivariable Cox proportional hazards regression, adjusted for all potential confounders. We tested for the proportional hazards assumption by using the Schoenfeld residuals and found a p value of 0.0003. As the proportional hazards assumption was violated and, thus, indicated non-proportional hazards, we stratified into different, non-overlapping time periods in 1-year intervals to utilize all available follow-up. Moreover, we plotted the cumulative incidence for the risk of study outcomes associated with the type of AD therapy received.
Subgroup analyses
Stratification was based on age (< 65 years versus ≥ 65 years), and sex, where a single model with interaction terms, was used to observe whether the association between exposure and outcome differed significantly among subgroups. Moreover, among AD+BZD recipients, we stratified on the type of BZD received (short-acting versus long-acting; patients who received both BZD types with antidepressants were excluded in this subgroup analysis) and patients who discontinued BZDs, defined as those who did not use benzodiazepines continuously for 6 months (180 days) after the index date; this definition was previously used by several studies and is considered a common definition [9, 22]. Benzodiazepine treatment length was defined using a 30-day grace period that was added to the days’ supply to allow for gaps in between prescriptions.
Sensitivity analyses
First, we varied the definition of concomitant therapy to at least one prescription of BZD within 2 weeks of the first AD prescription (Additional file 4: Fig. S3). Second, exposure was ascertained using an as-treated approach, where patients were censored when they have either discontinued treatment (defined as when no new AD or BZD prescription was given within 30 days of the end of the previous prescription [23]) or added BZDs throughout follow-up among AD monotherapy recipients. Third, we applied two other PS methods of inverse probability of treatment weighting and model adjustment [19, 24]. Lastly, we estimated the E-value to assess the potential impact of unmeasured confounders on our study findings [25], where in brief, a large E-value would imply that large unmeasured confounding would be needed to explain the observed association. All analyses were performed using the SAS Enterprise Guide program provided by the NHIS (Release 9.71, SAS Institute Inc., Cary, NC, USA).