Data source

We used Taiwan’s National Health Insurance Research Database (NHIRD) that covers 23.6 million individuals of the Taiwanese population. The NHIRD was derived from a mandatory single-payer health care insurance program, which includes complete health claim data of all residents with limited loss to follow-up (e.g., emigration). Denosumab is reimbursed by the national insurance program with very low copayment and can be identified by a unique national drug code in the database. It is an ideal setting for evaluation of the real-world effect of denosumab in reducing risk of fractures.

Another data source was the Chang Gung Research Database (CGRD) based upon electronic health record (EHR) data within the medical system of Chang Gung Medical Foundation (CGMF) [7]. CGMF is the largest medical system in Taiwan providing about 10% of all health care. In contrast to claims data, EHR data contain more clinical information, including weight, height, and results of bone mineral density (BMD) testing. CGRD and the main study data source of NHIRD are not directly linkable by researchers. However, since the data source of NHIRD covers 100% of Taiwan, patients in CGRD are by definition a subset of patients in NHIRD.

Study population

All men initiating denosumab in NHIRD were available for study inclusion [8]. Of note, the national health insurance program in Taiwan provides reimbursement for the use of denosumab among patients at high risk for fracture, defined by (1) a T-score ≤  − 2.5 SD and a vertebral or hip fracture, or (2) T-score − 1.0 SD to <  − 2.5 SD plus two vertebral or hip fractures [9], which means the study cohorts are all with fracture history (i.e., secondary fracture prevention population).

This study included men aged 50 years or older with osteoporosis who received at least one administration of denosumab, inclusive of years 2014 (initial year of denosumab use in men in Taiwan) through 2018, and with follow-up through 2019. To ensure that included men were receiving denosumab only for the indication of male osteoporosis, patients with a history of Paget’s disease or malignancy were excluded.

Study design

This was a retrospective cohort study. Similar to published literature [9], we measured the effectiveness of denosumab by comparing the incidence rates of clinical fractures between patients on denosumab 60 mg subcutaneously every 6 months (on-treatment) and patients ending therapy after one administration (off-treatment) (Fig. 1). One administration of denosumab, which has no known effect on fracture reduction [10], serves as a proxy for placebo exposure and also removes the possibility of bias resulting from differences in the initial treatment decision (i.e., limit confounding by indication).

Fig. 1

We defined day 1 as the day of patient’s first dose of denosumab. For the fracture analysis, index date was defined as day 225 (180 days for expected date of second administration + 45-day grace period). We followed patients from the index date to the change of assigned treatment, fracture endpoints, end date of available data, or death. The change of assigned treatment was discontinuation of denosumab for the on-treatment cohort and re-initiation of osteoporosis medication for the off-treatment cohort (Fig. 1). The 45-day grace period allowed for potential administrative challenges associated with return clinic visits and was based on the duration of denosumab activity.

Endpoints

The primary endpoint was hip fracture. Secondary endpoints included clinical vertebral fracture and nonvertebral fracture (hip, humerus, wrist, and distal forearm). All fractures were identified from inpatient claims to increase the likelihood of high specificity, albeit with decreased sensitivity. In the NHIRD, diagnosis codes were based on the International Classification of Diseases, Ninth Revision, Clinical Modification (ICD-9-CM) through December 2015 and on the International Classification of Diseases, Tenth Revision, Clinical Modification (ICD-10-CM) through the end of the study period. To increase the likelihood that fractures were due to osteoporosis, fractures concurrent with a motor vehicle accident on the same date were not included as endpoints.

The diagnosis of hip fracture was previously validated by comparing the operational definition of hip fracture versus medical chart review in NHIRD. Medical charts and radiographs of 300 randomly selected patients having either an ICD-9-CM or ICD-10-CM code for hip fracture at a large academic institution were previously reviewed by two physicians independently. Of these 300 patients, 298 were confirmed to have a hip fracture by the physicians (i.e., a positive predictive value of 99.3%) [9].

In addition, unadjudicated events of interest, including osteonecrosis of the jaw (ONJ), atypical femur fracture (AFF), and hypocalcemia, were identified using previously published algorithms [11,12,13]. A descriptive analysis of the incidence rate of the above events of interest was conducted that included all patients who had received at least one dose of denosumab and were followed up from day 1 through the earliest of denosumab discontinuation (6 months [recommended dosage] + 45-day grace period]), event of interest, death, or end of study period (December 2019).

To supplement the primary results of fracture, change in BMD was described in an on-treatment subgroup with both claims and EHRs from CGRD.

Covariates

Covariates were included to characterize the study population and for the application of statistical approaches accounting for any prognostic differences between study cohorts for the study endpoints. Fifty-nine fracture-related risk factors, including demographic characteristics and histories of comorbidities, medication use, and health-seeking behavior, in patients were assessed based on published literature and clinical knowledge. The 1-year period before patients’ initial use of denosumab was used to exclude patients with malignancy and Paget’s disease. The 1-year period before the index date (day 225) was used to describe all the covariates in the analysis except for previous fractures and osteoporosis drug use, which were described in the 3-year period before the index date.

Statistical analysis

Incidence rates were calculated as the number of patients with clinical osteoporotic fractures or specific unadjudicated events of interest per 100 person-years for fracture endpoints or per 10,000 person-years for unadjudicated events of interest, in the on-treatment and off-treatment cohorts. Patient characteristics were described, and standard mean differences between study cohorts were calculated.

The comparison of fracture incidence between the on-treatment and off-treatment cohorts used inverse probability of treatment weighting (IPTW) propensity score (PS) models [14] to adjust for the measured differences in baseline prognostic values. Two extreme weights were capped as 4. All patients were included to estimate the adjusted hazard ratio (HR) in the IPTW cohort using a Cox proportional hazard model. Results for treatment effectiveness (i.e., fracture risk reduction) were presented as HRs with 95% confidence intervals (CIs).

Assessment of unmeasured confounding

Five sensitivity analyses were conducted to address potential bias from unmeasured confounding.

For the first analysis, as an extension of PS in the primary analysis, the high-dimensional PS (HDPS) was calculated. The HDPS algorithm is an automated technique that empirically identifies potential confounders or proxies for confounders in longitudinal data sets; the algorithm assesses thousands of diagnosis, procedures, and drug-dispensing codes recorded in administrative databases and then selects the several hundred of those codes, as transformed into binary covariates, that appear to be confounders [15].

Two quantitative bias analysis were conducted to evaluate the extent of residual confounding that would be required to refute an observed difference in fracture incidence between cohorts. One was the rule-out method, described previously [16] and publicly available (www.drugepi.org), and has been applied extensively in the literature [17]. Another is the E-value method. The E-value method, unlike the rule-out method, does not assume the unmeasured confounder is binary, nor does it assume only a single unmeasured confounder, or any specified confounder. The E-value is defined as the minimum strength of association that an unmeasured confounder would need to have with both the treatment and outcome to refute a specific treatment-outcome association [18].

For the fourth analysis, a negative control outcome, which has no causal relationship to osteoporosis treatments of interest, was used to detect residual confounding [19]. The negative control outcome variable was early hip fracture during the initial 3 months of follow-up (a period with no expectation of fracture reduction attributable to denosumab).

Finally, an additional fifth analysis was conducted to evaluate the baseline balance of BMD and body mass index (BMI), as well as smoking and alcohol use between denosumab on-treatment and off-treatment cohorts in a subset of the study population within NHIRD claims who have EHRs in the CGRD [7].

Subgroup analysis and other sensitivity analysis

The results of fracture analysis were examined by subgroups, including (1) age (50–64 years, 65–74 years, ≥ 75 years), (2) history of bisphosphonate use before the first dose of denosumab (yes or no), and (3) fracture between the first and second doses of denosumab (yes or no).

To examine the robustness of controlling for measured variables in the primary endpoint, an alternative method of statistical analysis included the use of PS-matched cohorts with a 1:1 ratio.

The algorithm to identify the clinical vertebral fracture endpoint intentionally favored specificity over sensitivity. As a sensitivity analysis, any position of diagnosis in the inpatient claim was used to identify a clinical vertebral fracture endpoint to favor sensitivity.