Effects of bisphosphonates after denosumab discontinuation and treatment effect heterogeneity using causal machine learning.

Текст сохранён. Часть иллюстраций пока доступна только в первоисточнике; загрузка будет проверена повторно.

Introduction

Osteoporotic fractures (OFs) significantly impair quality of life, restrict physical activity, and impose a substantial economic burden 1 , 2 . This burden is particularly severe in older adults, where OFs are associated with increased mortality and a higher risk of subsequent fractures 3 – 6 . Preventing OFs is therefore a key clinical and public health priority. Clinical guidelines recommend assessing fracture risk in older adults and initiating appropriate pharmacologic interventions 7 , 8 .

Denosumab, a widely used antiresorptive agent, enhances bone mineralization by inhibiting bone resorption 9 . Unlike bisphosphonates (BPs), denosumab does not accumulate in bone tissue and can lead to rapid bone loss and increased fracture risk upon discontinuation due to a rebound effect 10 , 11 . Sequential treatment with zoledronate or alendronate has been proposed to counteract this rebound by preserving bone mineral density (BMD) 12 , 13 .

Despite this clinical rationale, real-world evidence on the effectiveness of BPs following denosumab discontinuation remains limited. Grassi et al. reported a significantly lower incidence of vertebral fractures in patients receiving BPs compared to untreated patients (odds ratio: 13.9; 95% CI 1.7–111.1; p = 0.014) 14 . However, the study’s retrospective chart review design introduced potential selection bias. Similarly, Burckhardt et al. demonstrated a protective effect of BPs in a Swiss cohort (Hazard ratio [HR]: 0.14; p < 0.001), but the study was monocentric and lacked adequate adjustment for confounding 15 .

While traditional statistical methods are well-established and widely used, they may be inadequate for adjusting complex confounding structures—particularly when accounting for interactions between covariates—and are especially limited in high-dimensional settings with many covariates 16 , 17 . In such contexts, machine learning techniques offer greater flexibility and improved performance in estimating causal effects. Specifically, in high-dimensional data, machine learning can systematically identify relevant covariates and capture complex interactions among them, thereby reducing model uncertainty and improving adjustment for confounding. However, conventional machine learning techniques are optimized for prediction rather than causal inference, limiting their ability to estimate treatment effects or account for counterfactual outcomes 16 – 18 . To address these limitations, recent advances in causal machine learning—such as causal forests and causal survival forests (CSFs)—enable direct estimation of individual-level treatment effects by combining flexible modeling with causal identification strategies 18 , 19 . These causal machine learning models enable estimation of the average treatment effect (ATE), providing an overall measure of treatment benefit across the population. More importantly, they can estimate conditional average treatment effects (CATEs), which quantify how treatment effects vary across subgroups defined by patient characteristics. This ability to uncover treatment effect heterogeneity is critical for advancing personalized medicine, as it allows the identification of patient subpopulations that may benefit more—or less—from a given therapy. Such insights are particularly valuable for clinicians and policymakers seeking to optimize treatment decisions and allocate healthcare resources more efficiently. Despite these promising capabilities, the application of causal machine learning to health-related outcomes studies remains relatively limited.

Therefore, we aimed to estimate the average treatment effects and heterogeneous effects of BPs on the risk of recurrent OFs after the discontinuation of denosumab using causal machine learning (Fig. 1 ).

Results

Baseline characteristics

Of a total of 778,768 patients aged ≥ 50 years diagnosed with OFs in 2015, 57,979 individuals received denosumab treatment between January 1, 2015 and May 31, 2021. Among them, 53.9% (N = 31,263) discontinued denosumab. The process of selecting the study population is shown in Fig. 2 . After matching for the index date, the BP and non-BP cohorts comprised 7138 (22.9%) and 24,079 (77.1%) patients, respectively. Mean (SD) age of the study population was 74.8 (8.7) years, with no statistically significant difference between the cohorts (74.7 [7.9] years for the BP cohort vs. 74.9 [8.9] years for the non-BP cohort; p = 0.199). Both the cohorts were predominantly composed of women, accounting for 94.4% and 90.7% of the BP and non-BP groups ( p < 0.001). Mean number of denosumab administrations before the cohort entry date differed between the cohorts, with 1.7 (1.1) and 1.8 (1.1) in the BP and non-BP cohorts, respectively, demonstrating a statistically significant difference ( p < 0.001). A significant difference was observed in the presence of a history of BPs, with 62.8% and 27.9% in the BP and non-BP groups ( p < 0.001), respectively. Detailed material on baseline characteristics is presented in Table 1 . Notably, several covariates—including duration of prior bisphosphonate use, history of prior steroid use, and history of vertebral fractures—showed substantial differences between the groups; however, these imbalances were corrected after weighting (Supplementary Figure S2 ).

Incidence of OFs recurrence

Recurrent OFs occurred in 472 (23.0%) and 1383 (26.2%) patients in the BP and non-BP cohorts ( p < 0.001), respectively. Mean fracture-free times from the index date were 1.7 (0.6) and 1.6 (0.7) years in the BP and non-BP cohorts ( p < 0.001), respectively. Incidence rate per person-year of recurrent OFs in the BP cohort was significantly lower than that in the non-BP cohort (0.14 [95% CI 0.13–0.15] vs. 0.17 [95% CI 0.16–0.18], p < 0.001). Based on the Kaplan–Meier curve analysis, the fracture-free probability of the BP cohort differed significantly from that of the non-BP cohort (log-rank test, p < 0.001) (Figure S1 in the Supplementary Material).

Average treatment effects by CSF

The initial model revealed that 43 of 60 covariates influenced both tree splitting and forest growth. The variance importance of the covariates is shown in Fig. 3 . Out of these, the top 50%, which amounted to 21 covariates, were selected for the final CSF. In the initial and final CSF samples, the characteristics of the BP and non-BP cohorts were compared because the standardized mean differences (SMD) was < 0.1. The SMDs of the final CSF are shown in Figure S2 in the Supplementary Material. Fracture-free probabilities and incidence rates for the BP cohort were higher than those of the non-BP cohort (difference in fracture-free probabilities, 0.032 [95% CI 0.022–0.042; p < 0.001]; incidence rate ratio (IRR), 0.795 [95% CI 0.734–0.859; p < 0.001]).

Inverse propensity-weighted multivariable Cox proportional hazard model

The proportional hazard assumption remained unviolated, as confirmed by examining the log–log plot. The BP use demonstrated a significant association with a reduced risk of recurrent OFs (HR, 0.703 [95% CI 0.641–0.771; p < 0.001]) (Table 2 ). Factors, such as a history of diabetes without complications, vertebral fracture, and age were significantly associated with an increased risk of recurrent OFs. Conversely, a higher number of prior denosumab administrations and history of dementia were significantly associated with a decreased risk of recurrent OFs.

Heterogeneous treatment effects

Individual CATE values ranged between − 0.096 and 0.137, indicating potential heterogeneity (Figure S3) . However, there was no significant difference between negative and positive treatment effects ( p = 0.393). Area under the targeting operator characteristic (TOC) curve was 0.010 (95% CI − 0.005 to 0.031), closely approaching zero. This implies minimal heterogeneity in the rank-weighted average treatment effect. The TOC curve is illustrated in Figure S4 in the Supplementary Material. Meanwhile, the graphical assessment stratified by covariates showed the potential presence of heterogeneity in several covariates: age, duration of prior BPs, history of diabetes without complications, peptic ulcer disease, and steroid use, and type of institution for the last denosumab administration (Figure S5 ). The CATEs for patients aged ≥ 80 years exceeded those for individuals < 80 years. Furthermore, the CATEs scores were lower in patients with a history of peptic ulcer disease, diabetes without complications, or prior steroid use than in those without such medical histories. Additionally, patients without prior BP treatment exhibited higher CATE levels than those with prior BPs treatment. Moreover, CATEs were higher among patients receiving their last denosumab dose at tertiary and general hospitals than among those administered at hospitals, convalescent hospitals, and clinics.

The results of the best linear projection BLP test for the CATE by covariates were consistent with the graphical assessment results (Table S3 ).

Sensitivity analysis

We conducted an as-treated sensitivity analysis in which follow-up was discontinued at the time of any treatment change after the index date. The results were consistent with those of the main analysis, supporting the robustness of our findings: IRR was 0.740 (95% CI 0.663–0.825; p < 0.001), and the HR was 0.705 (95% CI 0.528–0.941; p < 0.001). In the negative control outcome analysis, no significant differences were observed between the BP and non-BP cohorts. The IRR was 1.020 (95% CI 0.835–1.240; p = 0.834), and the HR was 0.986 (95% CI 0.789–1.232; p = 0.901), supporting the validity of our study design and analytic approach.

Discussion

In this nationwide retrospective cohort study, we applied a CSF model to evaluate the effect of BP use on the risk of recurrent OFs following denosumab discontinuation. BP use was causally associated with a significantly reduced risk of recurrent OFs, with a fracture-free probability difference of 0.032 between the BP and non-BP groups. This corresponded to an IRR of 0.795, indicating a 20.5% relative risk reduction. These findings support the use of BPs as an effective post-denosumab intervention to reduce the burden of recurrent fractures. Furthermore, we observed heterogeneity in treatment effects based on age, prior BP use, diabetes status, peptic ulcer disease, steroid use, and the type of treatment institution, underscoring the need for individualized therapy.

According to a Swiss retrospective study by Burckhardt et al., the cumulative incidence of vertebral fractures after denosumab discontinuation in women was reported as 10.3% 15 . The fracture-free probability for patients without BP after denosumab discontinuation was approximately 0.8, whereas that for patients with BP after denosumab discontinuation was approximately 0.9. In the present study, evaluating all types of recurrent OFs in both men and women, the cumulative incidence was 25.3%, exceeding the figures reported by Burckhardt et al. Although direct comparisons with previous studies are challenging due to differences in study design and population characteristics, the incidence reported in our study appears reasonable, considering the high incidence of recurrent OFs 6 , 20 , 21 . Furthermore, the observed trend toward higher odds of fracture-free status in the BP cohort than that in the non-BP cohort aligns with findings from previous investigations. Burckhardt et al. (2021) reported that BP lowered the risk of vertebral fractures seven-fold (HR, 0.14; p < 0.001), which was higher than the 1.4-fold (HR, 0.70; p < 0.001) reduction estimated in our study. Given that their study was a retrospective chart review, there is a possibility of bias in patient selection.

To validate the robustness of our findings and enable comparison with previous studies, we translated the absolute risk difference derived from the CSF model into an IRR, yielding a value of 0.795. This indicates a more conservative estimate than that of the inverse propensity-weighted Cox model. This effect estimate was more conservative than that obtained using an inverse probability–weighted Cox model, which may be attributed to the doubly robust nature of the CSF model. Specifically, the CSF simultaneously adjusts for both the probability of treatment assignment and the outcome (fracture-free time) based on covariates 18 , 19 . In addition, the CSF provides a key methodological advantage over traditional Cox models and standard machine learning algorithms: it corrects for censoring when estimating survival time 19 . This feature helps mitigate bias that may arise when censoring is associated with patient characteristics that differ between treatment groups 22 . Unlike traditional Cox models or typical machine learning approaches, the CSF model incorporates the probability of censoring, thereby mitigating this form of bias and improving the credibility of treatment effect estimates.

Beyond improving internal validity, the CSF model enabled us to estimate individual-level CATEs and explore treatment effect heterogeneity across subgroups. While there was no evidence of a consistent linear gradient in treatment effects across the entire population, subgroup analyses revealed greater benefit from BP therapy in patients over 80 years of age and in those without prior BP use within one year before their index fracture. In addition, the association between the number of denosumab doses and fracture risk may be partly explained by the long-term therapeutic effects of denosumab reported in previous studies, as well as greater treatment persistence and continuity of care among patients who received more injections 23 , 24 . The type of institution where the last denosumab dose was administered was also associated with fracture outcomes. Patients who discontinued denosumab in clinics may have been relatively healthier, which could partly explain the lower fracture risk observed in Cox regression models. However, when we applied the CSF model to estimate counterfactual average treatment effects, bisphosphonate initiation did not show additional benefit in patients who discontinued denosumab in clinics. This suggests that the association in the Cox model may reflect underlying patient differences in risk rather than a true treatment effect. Whereas the Cox model can only estimate the relative hazard of fracture at the population level, the CSF model allows estimation of counterfactual treatment effects at the individual level, thereby providing a more nuanced understanding of heterogeneity in treatment response. These findings demonstrate the ability of causal machine learning to uncover clinically meaningful variation in treatment response, which may help guide targeted intervention strategies.

Although this study offers valuable clinical and methodological insights, it has several limitations. First, it was based on insurance claims data, which limited the available covariates. Specifically, BMD data, a critical factor in previous studies on refracture risk, could not be collected. However, considering that the risk of OFs can be assessed without BMD measurements in clinical practice, the absence of BMD is expected to have a minimal effect on the study findings. Second, the definition of recurrent OFs, which is a primary outcome, relied solely on diagnostic codes, potentially leading to an overestimation of the incidence. To mitigate this, algorithms were applied that considered fractures in the same area after six months and in different areas after three months 25 – 27 . This approach was validated by previous study and expert opinions. Additionally, because this study focused on the difference in treatment effects between the two groups, the effect of refracture overestimation was expected to be balanced between them. Third, misclassification of the BP and non-BP cohorts was possible. Insurance coverage for BP medications requires a T-score of − 2.5 or lower. If the T-score improves beyond − 2.5 after discontinuing denosumab, but BP medication is deemed necessary based on clinical judgment, patients might continue treatment without insurance coverage. Future studies should use databases to verify information on medications not covered by insurance. Fourth, this study focused on the Asian population, which may limit the generalizability of the results to other populations. Future studies should use data from diverse regions to enhance the generalizability of our findings. Fifth, the definition of denosumab discontinuation was operationalized as the absence of prescriptions within one year of the last administration. While this approach aligns with prior claims-based studies and clinical input, rebound-associated fractures can occur as early as 60 days after discontinuation. Therefore, this definition may underestimate fracture risk and preferentially select healthier patients who survive longer after discontinuation. Sensitivity analyses using alternative definitions were not feasible in this dataset, and future studies should address this issue. Sixth, this study uses an early version of a novel methodology. Sixth, this study uses an early version of this new methodology. Although the causal survival forest offers advantages over conventional machine learning by providing effect sizes with confidence intervals, it remains subject to biases inherent in the dataset, such as distributional imbalance or modeling errors. The relative contribution of these sources of bias to the observed results cannot be fully determined. Therefore, while the CSF results are informative, they should be interpreted with caution. To support robustness, we validated the CSF findings using Cox regression, which showed consistent results. Further research will be needed to address these methodological challenges. As a result, it may have undiscovered limitations, and external validation is required to confirm the robustness and reliability of the findings. Further studies will be needed to validate these findings using recent data beyond the study period or external datasets from other countries.

These findings have important implications for clinical decision-making in the management of patients discontinuing denosumab. While current guidelines highlight the risk of rebound fractures, they offer limited guidance on tailoring sequential therapy. Our results suggest that BPs are not only effective but may offer greater benefits in specific subgroups—such as older adults or those without recent BP exposure. These insights support a more risk-stratified approach to post-denosumab therapy selection.

By identifying patient-level factors associated with differential BP response, our study contributes to the advancement of personalized osteoporosis care. Causal machine learning offers an innovative tool for subgroup-specific treatment effect estimation, providing actionable evidence beyond population averages. Such an approach may enable clinicians to better match patients with the most effective interventions, reducing fracture risk while avoiding unnecessary treatment.

Methods

An overall schematic of the study design is provided in Fig. 4 28 .

Ethics statement

The Institutional Review Board of Kyung Hee University (KHSIRB-23–568[EA]) approved this study an exemption from review. All methods were carried out in accordance with relevant guidelines and regulations. Informed consent was waived by the Institutional Review Board of Kyung Hee University, as this study used only de-identified secondary data obtained from the Health Insurance Review and Assessment Service (HIRA). The study strictly adhered to the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) reporting guidelines 29 .

Data source

This study used the HIRA database between January 1, 2015 and May 31, 2022 (M20230221004). The HIRA database in South Korea serves as a comprehensive repository for claims data generated in a healthcare reimbursement process 30 . Covering the entirety of the South Korean population through a universal coverage system, this database provides a wealth of healthcare service information, encompassing treatments, medications, procedures, and diagnoses for approximately 50 million beneficiaries.

Study population

The study cohort included patients aged ≥ 50 years who had an OF and subsequently initiated denosumab between January 1, 2015, and December 31, 2015, and discontinued denosumab treatment between January 1, 2015 and May 31, 2021. During this period, in Korea, the denosumab reimbursement criteria primarily included patients with a history of OFs. Therefore, we focused specifically on individuals who received denosumab treatment after OF. In Korea, denosumab is marketed under two formulations with distinct indications and drug codes. The formulation coded as 629001BIJ is indicated for the treatment of bone loss and OFs, whereas the formulation coded as 629002BIJ is indicated for the prevention of skeletal-related events in patients with multiple myeloma or bone metastases from solid tumors. To ensure a homogeneous study population, our analysis included only patients who received the formulation indicated for osteoporosis. The OFs encompassed fractures of the hip/femur, vertebrae, humerus, radius, tibia/fibula, and ankle 21 , 26 , 31 . Based on the opinion of orthopedic clinicians, we considered a low-energy fracture in individuals aged ≥ 50 years to be an OF, even without specific osteoporosis testing or diagnosis. Detailed diagnostic codes for OFs using the International Classification of Diseases 10 (ICD-10) are provided in Table S1 in the Supplementary Material. Discontinuation of denosumab was defined as the absence of any denosumab prescriptions within one year of the last denosumab administration 14 , 32 . To minimize confounding, we excluded patients with a diagnosis of cancer (excluding non-melanoma skin cancer), Paget’s disease (ICD-10 code: M88), or other metabolic bone diseases (ICD-10 codes: M89, M90, and Q78), including disorders of bone development and growth, hypertrophy of bone, osteolysis, osteonecrosis, osteitis deformans, or osteopathy.

Exposure

Exposure was defined as the use of BP at least once during the period from the last denosumab administration to nine months. This exposure assessment timeframe was established in accordance with the guidelines recommending the immediate administration of antiresorptive agents following denosumab discontinuation 7 . The practical guidelines recommend the administration of antiresorptive agents on the date of the next visit, occurring six months after the last denosumab dose. We established nine months as the exposure assessment window to account for potential delays in the administration schedule. Patients exposed to BPs, including alendronate, risedronate, ibandronate, pamidronate, and zoledronate, were categorized into a BP cohort. In contrast, those who were not exposed were categorized into a non-BP cohort. For the BP cohort, the index date was defined as the date of the first BP prescription. To ensure comparability, the index date for the non-BP cohort was aligned with that of the BP cohort by matching according to age, sex, and year of last denosumab administration 33 , 34 .

Outcomes of interest

Primary outcome of interest was OF recurrence. The definition of recurrent OFs was aligned with that of OFs using diagnostic codes (Table S1 in the Supplementary Material). To avoid potential duplication of recurrent OFs, we considered cases occurring at the same anatomical site as recurrent OFs only if they occurred at least six months after the previous OF. Additionally, recurrent OFs at different anatomical sites were counted if they occurred at least three months after the prior OF. These definitions were based on prior studies and expert opinions from five orthopedic specialists 25 – 27 . The outcome assessment window was from the index date to the end of the follow-up, defined as the earliest occurrence date among recurrent OFs, death, resumption of denosumab, 1 year after the index date, or the end of the study period. Given our emphasis on capturing fractures resulting from rebound effects following denosumab discontinuation, we observed the outcomes up to one year after the index date. Death was defined using ICD-10 codes I461, R96, R98, or R99, or if the result of treatment was coded as death 21 . The total number of recurrent OFs was divided by the sum of person-times in each cohort to calculate the incidence rates per person-year in the BP and non-BP cohorts.

Average treatment effects

The ATE of BPs for recurrent OFs was assessed using two models. We used CSF to estimate the differences in fracture-free probabilities between the BP and non-BP cohorts. The CSF is a doubly robust model that combines outcome and treatment models to enhance causal effect estimation in the presence of censoring and covariates 19 , 35 . Figure 1 illustrates the overall methodology used to estimate the treatment effect using the CSF. We then calculated the relative incidence rate, i.e., the IRR from the fracture-free probabilities in each cohort 36 . Moreover, we utilized conventional survival analysis, the multivariable Cox proportional hazard model (Cox model) to estimate the HR for the BP. The Cox model was adjusted using the inverse propensity score derived from the CSF to consider the probability of treatment assignment to either the BP or non-BP cohorts.

Heterogeneous treatment effects

We evaluated the heterogeneity in ATE based on the CATE in various ways 18 , 19 . First, the distribution of individual CATE was depicted using a histogram, and a t-test was conducted to verify the heterogeneity between the negative and positive treatment effects. Second, a TOC curve was used to confirm the presence of heterogeneity in ranked weighted average treatment effects. The TOC curve was interpreted as follows: the x-axis represents the upper fraction of descending-ranked individual CATEs, and the y-axis illustrates the differences in ATE between the subgroup corresponding to that fraction and entire patient population. For instance, if the Y-axis is one and X-axis 0.2, it signifies that the difference in the ATE of the subgroup comprising the top 20% of CATE and that of the entire population would be one. Consequently, the area under the TOC curve approaching zero suggests minimal heterogeneity in the rank-weighted average treatment effects. Third, box plots were used to visually assess the heterogeneity. Finally, the BLP test was applied to assess the heterogeneous treatment effect stratified by covariates 18 , 37 , 38 .

Covariates

Covariates were collected based on previous studies on fractures 14 , 15 . Most covariates were assessed over the year preceding the index date. However, the number of prior denosumab administrations was tracked from January 1, 2015, until the date of the last denosumab administration. The duration of BP administration was evaluated for the year before the last denosumab administration. In total, 45 covariates, including sociodemographic factors, medication, and disease factors, were used for the analysis. Disease-related covariates were defined using the ICD-10 codes 39 , 40 . Specific details regarding the types and codes of these covariates can be found in Table S2 in the Supplementary Material.

Statistical analysis

Continuous variables were summarized using means and standard deviations (SD), while categorical variables were expressed as frequencies and proportions. Differences between groups were assessed using independent t-tests for continuous variables and chi-square tests for categorical variables. We calculated SMD, deeming an SMD below 0.1 as indicative of comparability. Furthermore, visual balance assessment was conducted using an overlapping histogram.

A Kaplan–Meier curve was plotted to visually assess the fracture-free probability in the BP and non-BP cohorts. The log-rank test was used to compare the fracture-free curves between the two cohorts. The CSF, akin to a surviving forest, iteratively divides the sample into progressively smaller subgroups, generating a sequence of decision trees that result in subgroups. While survival forests split data to maximize the survival time between subgroups, CSF makes splitting decisions to maximize the treatment effect, given its capability to directly estimate treatment effects. To address potential overfitting and biased treatment effect estimation, an honest CSF was employed, with the forest comprising 2000 trees 19 . We identified the top 50% of covariates based on their importance variance scores obtained from the initial CSF that employed all covariates. This step aimed to exclude unnecessary covariates from the estimation of treatment effects. Subsequently, we regenerated the CSF using this subset of selected covariates. A log–log plot was used to assess the proportional hazards assumption for the Cox model. The Cox model employs covariates selected by the CSF to enhance model performance.

We performed an as-treated sensitivity analysis in which follow-up was discontinued at the time of any treatment change after the index date to account for potential bias introduced by switching therapies. In addition, model validation was undertaken through an analysis using a negative control outcome, specifically “tinea pedis,” defined by the ICD-10 code B353. All statistical analyses were conducted using two-sided tests, with the significance level set at p < 0.05. The analysis was performed using SAS Enterprise Guide version 7.1 (SAS Institute Inc., Cary, NC) and R, version 4.3.1, (R Foundation for Statistical Computing). The CSF was run using the grf package version 2.3.0.

Supplementary Information

Below is the link to the electronic supplementary material.

Supplementary Material 1

Acknowledgements

This study was supported by a grant (21153MFDS601)from the Ministry of Food and Drug Safety in 2025. This work was supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government (RS-2024-00331719 and RS-2024-00345981). This study acknowledges Editage ( www.editage.com ) for providing professional language editing services to enhance the clarity and readability of the manuscript.

Author contributions

*Concept and design*: M Kim, N Jeon, and HS Suh. *Acquisition, analysis, or interpretation of data*: M Kim. *Drafting of the manuscript*: M Kim. *Critical review of the manuscript for important intellectual content*: M Kim, N Jeon, and HS Suh. *Statistical analysis*: M Kim. *Administrative, technical, or material support*: M Kim. *Supervision*: HS Suh.

Data availability

The datasets generated and/or analyzed during the current study are not publicly available due to privacy and ethical restrictions, but are available from the corresponding author, Hae Sun Suh (haesun.suh@khu.ac.kr), on reasonable request.

Declarations

Competing interests

The authors declare no competing interests.

Footnotes

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Burge, R. et al. Incidence and economic burden of osteoporosis-related fractures in the United States, 2005–2025. J. Bone Miner. Res. 22 , 465–475 (2007).

2. Mohd-Tahir, N. & Li, S. Economic burden of osteoporosis-related hip fracture in Asia: A systematic review. Osteoporos. Int. 28 , 2035–2044 (2017).

3. Berry, S. D., Kiel, D. P. & Colón-Emeric, C. Hip fractures in older adults in 2019. JAMA 321 , 2231–2232 (2019).

4. Bliuc, D. et al. Mortality risk associated with low-trauma osteoporotic fracture and subsequent fracture in men and women. JAMA 301 , 513–521 (2009).

5. Crandall, C. J. et al. Risk of subsequent fractures in postmenopausal women after nontraumatic vs traumatic fractures. JAMA Intern. Med. 181 , 1055–1063 (2021).

6. Söreskog, E. et al. Risk of major osteoporotic fracture after first, second and third fracture in Swedish women aged 50 years and older. Bone 134 , 115286 (2020).

7. Compston, J. et al. UK clinical guideline for the prevention and treatment of osteoporosis. Arch. Osteoporos. 12 , 1–24 (2017).

8. Qaseem, A. et al. Pharmacologic treatment of primary osteoporosis or low bone mass to prevent fractures in adults: A living clinical guideline from the American college of physicians. Ann. Intern. Med. 176 , 224–238 (2023).

9. Baron, R., Ferrari, S. & Russell, R. G. G. Denosumab and bisphosphonates: Different mechanisms of action and effects. Bone 48 , 677–692 (2011).

10. Bone, H. G. et al. Effects of denosumab treatment and discontinuation on bone mineral density and bone turnover markers in postmenopausal women with low bone mass. J. Clin. Endocrinol. Metab. 96 , 972–980 (2011).

11. Cummings, S. R. et al. Vertebral fractures after discontinuation of denosumab: A post hoc analysis of the randomized placebo-controlled FREEDOM trial and its extension. J. Bone Miner. Res. 33 , 190–198 (2018).

12. Anastasilakis, A. D., Papapoulos, S. E., Polyzos, S. A., Appelman-Dijkstra, N. M. & Makras, P. Zoledronate for the prevention of bone loss in women discontinuing denosumab treatment. A prospective 2-year clinical trial. J. Bone Miner. Res. 34 , 2220–2228 (2019).

13. Tutaworn, T. et al. Bone loss after denosumab discontinuation is prevented by alendronate and zoledronic acid but not risedronate: A retrospective study. Osteoporos. Int. 34 , 573–584 (2023).

14. Grassi, G. et al. Bisphosphonates after denosumab withdrawal reduce the vertebral fractures incidence. Eur. J. Endocrinol. 185 , 387–396 (2021).

15. Burckhardt, P., Faouzi, M., Buclin, T., Lamy, O. & Group, S. D. S. Fractures after denosumab discontinuation: A retrospective study of 797 cases. J. Bone Miner. Res. 36 , 1717–1728 (2021).

16. Lang, M. et al. Automatic model selection for high-dimensional survival analysis. J. Stat. Comput. Simul. 85 , 62–76 (2015).

17. Wang, P., Li, Y. & Reddy, C. K. Machine learning for survival analysis: A survey. ACM Comput. Surv. (CSUR) 51 , 1–36 (2019).

18. Athey, S., Tibshirani, J. & Wager, S. Generalized random forests. (2019).

19. Cui, Y., Kosorok, M. R., Sverdrup, E., Wager, S. & Zhu, R. Estimating heterogeneous treatment effects with right-censored data via causal survival forests. J. R. Stat. Soc. Ser. B Stat Methodol. 85 , 179–211 (2023).

20. Hansen, L. et al. Subsequent fracture rates in a nationwide population-based cohort study with a 10-year perspective. Osteoporos. Int. 26 , 513–519 (2015).

21. Shim, Y.-B. et al. Incidence and risk factors of subsequent osteoporotic fracture: A nationwide cohort study in South Korea. Arch. Osteoporos. 15 , 1–7 (2020).

22. Howe, C. J., Cole, S. R., Lau, B., Napravnik, S. & Eron, J. J. Jr. Selection bias due to loss to follow up in cohort studies. Epidemiology 27 , 91–97 (2016).

23. Bone, H. G. et al. 10 years of denosumab treatment in postmenopausal women with osteoporosis: Results from the phase 3 randomised freedom trial and open-label extension. Lancet Diabetes Endocrinol. 5 , 513–523 (2017).

24. Ferrari, S. & Langdahl, B. Mechanisms underlying the long-term and withdrawal effects of denosumab therapy on bone. Nat. Rev. Rheumatol. 19 , 307–317 (2023).

25. Banefelt, J. et al. Risk of imminent fracture following a previous fracture in a Swedish database study. Osteoporos. Int. 30 , 601–609 (2019).

26. Lee, S.-B. et al. Association between mortality risk and the number, location, and sequence of subsequent fractures in the elderly. Osteoporos. Int. 32 , 233–241 (2021).

27. Weaver, J., Sajjan, S., Lewiecki, E. M., Harris, S. T. & Marvos, P. Prevalence and cost of subsequent fractures among US patients with an incident fracture. J. Manag. Care Spec. Pharm. 23 , 461–471 (2017).

28. Schneeweiss, S. et al. Graphical depiction of longitudinal study designs in health care databases. Ann. Intern. Med. 170 , 398–406 (2019).

29. Vandenbroucke, J. P. et al. Strengthening the reporting of observational studies in epidemiology (STROBE): Explanation and elaboration. Ann. Intern. Med. 147 , W-163-W−194 (2007).

30. Kim, J.-A., Yoon, S., Kim, L.-Y. & Kim, D.-S. Towards actualizing the value potential of Korea health insurance review and assessment (HIRA) data as a resource for health research: Strengths, limitations, applications, and strategies for optimal use of HIRA data. J. Korean Med. Sci. 32 , 718–728 (2017).

31. Wright, N. C. et al. The design and validation of a new algorithm to identify incident fractures in administrative claims data. J. Bone Miner. Res. 34 , 1798–1807 (2019).

32. Ha, J. et al. Effect of follow-up raloxifene therapy after denosumab discontinuation in postmenopausal women. Osteoporos. Int. 33 , 1591–1599 (2022).

33. Iwagami, M. & Shinozaki, T. Introduction to matching in case-control and cohort studies. Ann. Clin. Epidemiol. 4 , 33–40 (2022).

34. Song, J. W. & Chung, K. C. Observational studies: Cohort and case-control studies. Plast. Reconstr. Surg. 126 , 2234 (2010).

35. Funk, M. J. et al. Doubly robust estimation of causal effects. Am. J. Epidemiol. 173 , 761–767 (2011).

36. Briggs, A. & Sculpher, M. An introduction to Markov modelling for economic evaluation. Pharmacoeconomics 13 , 397–409 (1998).

37. Semenova, V. & Chernozhukov, V. Debiased machine learning of conditional average treatment effects and other causal functions. Economet. J. 24 , 264–289 (2021).

38. Zhang, Y., Li, H. & Ren, G. Estimating heterogeneous treatment effects in road safety analysis using generalized random forests. Accid. Anal. Prev. 165 , 106507 (2022).

39. Li, B., Evans, D., Faris, P., Dean, S. & Quan, H. Risk adjustment performance of Charlson and Elixhauser comorbidities in ICD-9 and ICD-10 administrative databases. BMC Health Serv. Res. 8 , 1–7 (2008).

40. Thygesen, S. K., Christiansen, C. F., Christensen, S., Lash, T. L. & Sørensen, H. T. The predictive value of ICD-10 diagnostic coding used to assess Charlson comorbidity index conditions in the population-based Danish National registry of patients. BMC Med. Res. Methodol. 11 , 1–6 (2011).

Associated Data

Supplementary Materials

Supplementary Material 1

Data Availability Statement

The datasets generated and/or analyzed during the current study are not publicly available due to privacy and ethical restrictions, but are available from the corresponding author, Hae Sun Suh (haesun.suh@khu.ac.kr), on reasonable request.