PDF
Citation
Print
INTRODUCTION
Thyroid cancer is one of the most prevalent endocrine malignancies worldwide and is the most frequently diagnosed cancer in South Korea [1]. It primarily affects individuals in their 40s, with a 5-year survival rate exceeding 98%, enabling most patients to achieve long-term survival [2]. As a result, the longterm management of post-treatment complications has become a major concern, affecting not only individual patients but also the healthcare system burden.
Osteoporosis is a frequent complication among thyroid cancer survivors who undergo thyroid-stimulating hormone (TSH) suppression therapy, which is used to inhibit cancer cell growth and prevent recurrence [3-5]. Although osteoporosis is commonly associated with aging, thyroid cancer survivors display distinct incidence patterns compared with the general population. A nationwide cohort study reported a significantly increased risk of osteoporosis in female thyroid cancer patients, particularly in premenopausal women [6]. In addition, a meta-analysis revealed reductions in bone mineral density (BMD) at the lumbar spine and hip in both premenopausal and menopausal women receiving long-term TSH suppression therapy [7]. TSH exerts substantial effects on bone metabolism beyond its role in thyroid hormone regulation. It decreases bone resorption by inhibiting osteoclast development and activity, while simultaneously promoting bone formation by stimulating osteoblast growth and differentiation [8]. Consequently, levothyroxine administration for TSH suppression accelerates bone remodeling and resorption, ultimately leading to reduced BMD [9].
Early diagnosis and management of osteoporosis in patients with thyroid cancer are crucial. Osteoporosis-induced fractures cause pain, dysfunction, decreased quality of life, delayed cancer rehabilitation, and increased medical expenditures. A recent cohort study further demonstrated that fracture incidence varies depending on whether patients receive active osteoporosis treatment [10]. This highlights the urgent need for precise management through a prediction–selection–intervention approach, which requires early identification of bone health risks in thyroid cancer survivors and the application of targeted pharmacological and lifestyle interventions for high-risk groups.
Previous research on osteoporosis risk in thyroid cancer patients has largely focused on prevalence and associated risk factors, typically using regression-based statistical analyses, but has not extended to predictive modeling [6,11]. Although recent studies have developed prediction models for osteoporosis in other cancer populations, these have mainly used cross-sectional designs with baseline classification machine learning approaches [12-14].
In this nationwide cohort study, we aimed to develop and validate three time-to-event prediction models for osteoporosis in thyroid cancer survivors: random survival forest (RSF), Cox proportional hazards with Boruta feature selection (Boruta-Cox), and least absolute shrinkage and selection operator (LASSO)-penalized Cox proportional hazards (LASSO-Cox). Risk stratification was performed to evaluate the clinical applicability of these models. This work is expected to inform strategies for managing long-term complications in thyroid cancer patients and to guide the design of personalized treatment approaches.
METHODS
Data source
This study used data from the Korea National Health Insurance Service–Sample Cohort, which includes 1,134,108 individual records containing demographic qualifications, health examination results, medical diagnoses (based on the International Classification of Diseases, 10th Revision), and prescription information spanning January 1, 2002, to December 31, 2019.
Study design
We identified patients newly diagnosed with thyroid cancer (C73) between January 1, 2004, and December 31, 2014. Patients with a prior diagnosis of thyroid cancer before 2004 were excluded to ensure inclusion of only newly diagnosed cases. Additional exclusion criteria were applied (Fig. 1): (1) history of osteoporosis; (2) history of any other cancer; (3) death between the cohort entry date and the index date; (4) age younger than 19 or older than 80 years; (5) no health checkup before the first thyroid cancer diagnosis; (6) development of osteoporosis between the cohort entry date and the index date; and (7) missing or outlier data. Patients were followed from diagnosis until the onset of osteoporosis, death, or December 31, 2019.
The cohort was divided into training and test datasets by stratified splitting at a 7:3 ratio while maintaining the distribution of event occurrences. A recent machine learning study reported that allocating 70% of the data for training provided the best results [15]. Definitions of variables and outcomes are summarized in Supplemental Table S1, and details of the study design are presented in Supplemental Fig. S1.
Model design
We developed and compared three different time-to-event prediction models: RSF, Boruta-Cox, and LASSO-Cox. Model development and feature selection were performed exclusively on the training set, while performance was evaluated using both training and test sets. To prevent scale-dependent penalty bias, all continuous variables were standardized to z-scores using the distributions of the training and test datasets.
Before modeling, hyperparameter tuning was conducted for the RSF and LASSO-Cox models. For the RSF model, a grid search strategy was implemented on the training dataset, exploring mtry (4, 5, 6, 7, 8), ntree (50, 100, 300, 500), and nodesize (5, 10, 15). Five-fold cross-validation was performed for each parameter set, and the configuration yielding the highest mean C-index (mtry=7, ntree=50, nodesize=15) was selected. For the LASSO-Cox model, the regularization parameter lambda (λ) was tuned using 10-fold cross-validation on the training set, with the optimal λ=0.009902 chosen to minimize the cross-validated partial likelihood error.
Three distinct feature-selection strategies were applied, each tailored to the corresponding model. First, an RSF model was trained on the training set, and the top 10 most important variables were selected for prediction. Second, Boruta was applied to construct an optimized Cox proportional hazards model, using only variables classified as ‘Confirmed.’ From these, the 10 features with the highest importance scores were retained. For both RSF and Boruta, limiting the final models to the top 10 variables followed the events-per-variable (EPV) ≥10 rule and ensured parsimony [16,17]. Finally, the LASSO penalty was applied to the training set for variable selection, which was then used to construct a Cox proportional hazards model with fivefold cross-validation.
Model performance
Model performance was evaluated on both training and test datasets. We calculated the time-dependent area under the curve (AUC) at 1-year intervals from 2 to 5 years after the cohort entry date. Discriminative ability was also assessed using the concordance index (C-index). A C-index of 0.5 indicated random prediction, whereas a value of 1.0 indicated perfect prediction [18].
Risk stratification
Risk stratification was conducted to ensure that the findings were not model-specific. Risk scores for each patient were calculated from the three models. The median predicted risk score of the training set was used as the cutoff, which was also applied to the test set. Patients were stratified into low-risk and high-risk groups based on this cutoff, facilitating clinically actionable decision-making and clearer differentiation between risk groups [18].
Statistical analysis
Differences in demographic and clinical characteristics were compared using the Wilcoxon test for continuous variables and the chi-square test or Fisher exact test for categorical variables in both the training and test sets. A two-tailed P>0.05 indicated no significant difference between groups, supporting the appropriateness of the stratified split. Baseline characteristics were described using frequency and percentage for categorical variables, and median and interquartile range for continuous variables. For stratification into high- and low-risk groups, survival distributions were estimated using the Kaplan–Meier method, and differences in survival probability were tested with the log-rank test. Preprocessing was performed with SAS software version 9.4 (SAS Institute Inc., Cary, NC, USA). Statistical analyses and predictive modeling, including model development and evaluation, were conducted in R version 4.3.2 (R Foundation for Statistical Computing, Vienna, Austria).
Approval
This study was approved by the Institutional Review Board of Konkuk University Medical Center (KUMC 2024-08-018). The study was conducted in strict accordance with the Declaration of Helsinki. Because the National Health Insurance Service data were anonymized and researcher access was strictly controlled, the requirement for informed consent was waived by the National Health Insurance Service.
RESULTS
Participants’ characteristics
A total of 3,089 patients were included in the study. No statistically significant differences were observed between the training and test sets in terms of baseline characteristics, confirming their suitability for comparative analysis (Table 1). The median age of participants was in the mid-40s, and males accounted for approximately 25% of the cohort. Osteoporosis occurred in 21.4% of cases overall (21.4% in the training set; 21.6% in the test set), and the median follow-up duration was 4.2 years.
Model performance
All variables were incorporated into feature selection to determine the relative importance of predictors for osteoporosis risk (Fig. 2). In the test set, the Boruta-Cox model consistently demonstrated superior discrimination accuracy across all time points (Fig. 3). Time-dependent AUC values for the three models in the training set are presented in Supplemental Fig. S2. Moreover, the Boruta-Cox model achieved the highest C-index in the test set, while the RSF model exhibited the best performance in the training set (Table 2).
Risk stratification
Risk scores for individual patients were calculated with each model, and participants were stratified into high-risk and low-risk groups using the designated cutoff point. Kaplan–Meier survival analysis showed a significant difference between the two groups, as confirmed by the log-rank test (P<0.001) (Fig. 4). Hazard ratio analysis further demonstrated that the high-risk group experienced a significantly higher incidence of events (Table 3). Results of risk stratification in the training set are provided in Supplemental Table S2, Supplemental Fig. S3. Overall, stratification effectively distinguished the potential risk of osteoporosis between low- and high-risk groups across all three models.
DISCUSSION
This study developed and evaluated three time-to-event prediction models for assessing long-term osteoporosis risk in patients with thyroid cancer. Our findings demonstrated that the Boruta-Cox model provided superior performance in predicting the 2-, 3-, 4-, and 5-year risk of osteoporosis in this population. Although the RSF model showed the highest C-index and AUC across all time points in the training set, it appeared to be overfitted. We excluded cases of 1-year events to minimize the potential influence of baseline fragility, as therapy-related bone loss can begin within months of TSH suppression. Several studies have reported no significant risk of osteoporosis or clinically meaningful BMD loss within the first year of TSH-suppressive therapy, provided that baseline bone status is normal [19-21]. Therefore, the predictive scope of our models applies to patients followed for longer than 1 year, making them unsuitable for evaluating early-onset risk within the initial year.
Across the three feature-selection models, variables selected two or more times included age, sex, number of outpatient visits, height, diastolic blood pressure, and diabetes (Supplemental Fig. S4). Prior research has shown that osteoporosis prevalence increases significantly with age and remains consistently higher in women than in men across all age groups [22,23]. Hypertension has also been reported to be significantly associated with osteoporosis risk [24-26]. Furthermore, diabetes contributes to osteoporosis through mechanisms involving hyperglycemia, inflammation, and hormonal dysregulation [27,28]. The number of outpatient visits may be a relevant factor, as frequent visits increase opportunities for diagnostic services such as dual-energy X-ray absorptiometry (DXA), which measures BMD. Height also emerged as an independent predictor in our models. Nevertheless, these variables primarily function as proxy markers of risk rather than direct biological mechanisms of osteoporosis.
The models demonstrated strong risk stratification, supporting early identification of high-risk patients and prioritization of further diagnostic assessment. Looking forward, our prediction models could serve as initial screening tools for estimating osteoporosis risk in thyroid cancer survivors. A two-step approach could be implemented, in which DXA testing is reserved for those identified as high risk by the prediction model [28,29]. Such a strategy would help translate epidemiologic evidence into clinical practice. To enable clinical application of the Boruta-Cox model, we plan to develop an interactive web-based calculator using R/Shiny [30]. Clinicians would enter the 10 predictors to receive both a continuous 5-year risk estimate and a traffic-light risk classification, alongside prompts to order DXA scans, check vitamin D levels, and initiate first-line antiresorptive therapy in high-risk patients. In parallel, we will containerize the model as a Fast Healthcare Interoperability Resources (FHIR)-compliant Representational State Transfer (RESTful) Application Programming Interface (API) to facilitate seamless Electronic Medical Record (EMR) integration [31]. With this integration, osteoporosis risk could be calculated automatically at the time of levothyroxine prescription or DXA order and displayed as an inline decision-support gauge without requiring additional clicks. Finally, we plan to conduct a stepped-wedge trial involving approximately 1,000 thyroid cancer survivors to prospectively evaluate the real-world impact of this model on guideline-concordant bone-density screening and initiation of preventive therapies [32]. Ultimately, this work will validate the predictive capacity of our model while also demonstrating its ability to improve patient outcomes and enhance clinician workflow efficiency in routine practice.
This study has several limitations. First, it utilized National Health Insurance Service claims data, which lack certain variables such as menopausal status, TSH level, and BMD. Furthermore, we characterized TSH suppression treatment based on prescriptions of levothyroxine rather than confirmed TSH suppression. The claims data do not specify whether levothyroxine prescriptions were intended for replacement or suppressive therapy. However, because levothyroxine is most commonly prescribed to thyroid cancer patients for TSH suppression rather than simple replacement, we used its prescription as a proxy indicator of TSH suppression in this study. These missing variables may represent critical predictors of osteoporosis risk. Consequently, our findings should be interpreted as applicable to the general thyroid cancer population rather than as precise estimates for subgroups such as pre- or postmenopausal women. Future research should investigate differences in osteoporosis risk factors across subgroups of thyroid cancer patients. Such subgroup analyses would facilitate the selection of relevant predictors and the development of models tailored to premenopausal women, postmenopausal women, and men. This approach would help identify the most influential risk factors in each subgroup and enable construction of subgroup-specific prediction models. Because osteoporosis risk varies substantially across these groups, subgroup-specific models could provide more clinically meaningful insights. Second, the constructed prediction models were not designed to determine causality. Therefore, the study does not provide causal inferences regarding the impact of specific predictors on osteoporosis risk in thyroid cancer patients. Third, the study may be subject to selection bias due to the exclusion of 1-year events. We excluded osteoporosis events occurring within the first 12 months after thyroid cancer diagnosis to minimize reverse causation. However, this criterion could have introduced bias by omitting patients with subclinical bone loss or rapid treatment-related bone changes. Thus, our risk estimates are valid only for osteoporosis developing after the first year following diagnosis. Finally, this study lacked an external validation cohort, such as hospital-based datasets. Although the models underwent 5-fold cross-validation, all folds were drawn from the same source population. Future studies should perform external validation using independent cohorts and update the models before clinical implementation. We plan to deploy the prediction model as a RESTful API service that can be integrated into hospital EMR systems. The model will be prospectively applied to newly diagnosed thyroid cancer patients, and predicted outcomes will be compared with observed 5-year osteoporosis incidence. Conducting external validation is crucial for ascertaining model applicability across various healthcare environments. Until such evidence is available, the tool should be regarded as exploratory and applied cautiously. Only after consistent performance is demonstrated in independent cohorts should the model be integrated into standard clinical practice.
This study also has notable strengths. First, we developed prediction models for long-term osteoporosis risk in thyroid cancer survivors, addressing an important clinical gap. Second, by employing longitudinal data and time-to-event survival models, we enabled dynamic prediction of osteoporosis onset across multiple clinically meaningful time points. Third, we constructed and compared three distinct survival-based prediction models, demonstrating robustness of findings and offering clinicians algorithmic options that balance interpretability with predictive accuracy.
In this study, we employed the RSF, Boruta-Cox, and LASSO-Cox algorithms to develop predictive models for time-to-osteoporosis in thyroid cancer survivors. Among these, the Boruta-Cox model demonstrated superior performance in both C-index and time-dependent AUC. Furthermore, we conducted risk stratification with all three models, showing their potential as primary screening tools to identify thyroid cancer patients at high risk of developing osteoporosis. However, before clinical integration, external validation and subgroup analyses stratified by sex and menopausal status in independent cohorts are essential to ensure generalizability and safe application in routine practice.
XML Download