This retrospective study used a delinked clinical research database from 3 hospitals in the Taipei Medical University health care system, including Taipei Medical University Hospital, Wan Fang Hospital, and Shuang Ho Hospital. The study was approved by the Taipei Medical University Joint Institutional Review Board (N202107054). As the data were deidentified, the requirement for informed consent was waived. This study adhered to the TRIPOD (Transparent Reporting of a Multivariate Prediction Model for Individual Prognosis or Diagnosis) checklist [40].

Figure 1 shows the study design, including data collection, feature selection, model construction, and 5-fold internal validation using a training set, external validation by a test set, and model interpretation. The training set of the study comprised patients from Wan Fang Hospital and Taipei Medical University Hospital, whereas the external test set comprised patients from Shuang Ho Hospital. There were 16 machine learning models built by the training set. Model performance was compared with the external testing set.

Patients older than 18 years who had first been prescribed oral amiodarone between January 2005 and December 2017 were included. Patients were excluded if they were pregnant or had a history of a thyroid disorder diagnosis, thyroid surgery, or subclinical thyroid laboratory data (Multimedia Appendix 1) within 1 year before the index day, which was the day of the first amiodarone prescription. The predefined inclusion and exclusion criteria were modified from previous medication-induced thyroid dysfunction research [6,27]. The study subjects were followed up for 2 years after the index date. This study collected features until the end of the study period and recorded patient loss to follow-up and the occurrence of AITD events.

The study collected stationary and dynamic features from the clinical research database. Stationary features, including sex, age, and BMI, were collected once at the index date. Dynamic features, including laboratory tests, comorbidities, and comedications, were continuously collected to reflect the patient’s clinical condition. The collection points for dynamic features were at baseline and the first, second, third, sixth, ninth, 12th, 15th, 18th, and 21st months after amiodarone initiation. Details of the dynamic data collection are shown as a diagram in Multimedia Appendix 2. The features reengineered by the research team included the accumulated or average dose of amiodarone and annual time-series variations in laboratory tests.

Robust scaler algorithms were used to normalize the data and reduce the effect of extreme numeric variables and large differences in range between laboratory values [41]. Variables with more than 90% missing data were deleted from the machine learning programs. Missing values were first imputed with the last observation carried forward, and then the remainders were imputed with Multivariate Imputation by Chained Equations (Scikit-learn), an open-source imputation software package [42-44]. Zero was imputed for laboratory data with changing rates, such as the trend or slope of lab values. The codebook and missing rates in the training and test sets are provided in Multimedia Appendix 3.

The outcome for prediction in this study was the occurrence of AITD, which was defined by a thyroid function test and with diagnostic criteria from previous studies [6,22,26,45]. Cases of AITD were identified if the thyroid-stimulating hormone (TSH) titer was <0.1 mU/l and the free thyroxin (fT4) level was higher than the normal range, while cases of amiodarone-induced hypothyroidism were ascertained by a serum TSH titer of >10 mU/l, regardless of the fT4 level; a TSH titer of 4 to 10 mU/l with a lower than normal fT4 level; an International Classification of Diseases (ICD)-9 code 242, 243, or 244 or an ICD-10 code E02, E03, E05, or E06; having received pharmacotherapy (eg, levothyroxine; Propylthiouracil, carbimazole, or methimazole) for thyroid disease; or an ICD-9 procedure code for a thyroidectomy. As TSH and fT4 are definite diagnostic criteria for AITD, the last data values before the prediction point were masked to avoid data leakage in the test set. Once a patient developed AITD, the data at that time point were coded 1, and it was designated the earliest AITD onset date.

Hyperparameter tuning of the XGBoost, AdaBoost, KNN, and LR algorithms was performed with an exhausted-grid search toward maximizing F₁-score metrics. Five-fold cross-validation was performed inside each grid option, and the optimal hyperparameter set was chosen based on the model in the grid search with the highest F₁-score. XGBoost was tuned on 7 hyperparameters, including max_depth, min_child_weight, gamma, subsample, colsample_bytree, n_estimators, and learning_rate for 94,080 grid options. Three hyperparameters of AdaBoost that were tuned included n_estimators, learning rate, and algorithm, with 160 combinations. As the KNN is based on the KNN of the prediction point, n_neighbors, weights, and metric, combining 120 sets of parameters, were tested. Finally, the study used an LR established with the scikit-learn module for binary outcome classification. P Penalty, solvers, and C, hyperparameters of LR, were calculated in 140 selections. Details of the hyperparameters and the final best combination of the above 16 models are shown in Multimedia Appendix 4.

Recursive feature elimination (RFE) with cross-validation was used for the training set. Pseudocodes of the grid search and cross-validation are presented in Multimedia Appendix 5. The minimal number of feature sets was generated by XGBoost and AdaBoost. As the incidence of AITD in the training and test sets was 14.3% and 11.4%, respectively, the imbalance issue was managed by (1) oversampling with the borderline synthetic minority oversampling technique (B-SMT) [46]; (2) undersampling the majority class with ENN [47]; and (3) a combination of oversampling and undersampling with B-SMT-ENN. The raw strategy and 3 resampling strategies were applied to XGBoost, AdaBoost, KNN, and LR, as shown in Figure 1. This study finally constructed 16 models through the 5-fold cross-validation process of the training sets.

The performance of the 16 models generated with the training data set were validated and evaluated on the test data set. Model performance was compared using accuracy; Precision, recall, F₁-score, geometric mean, AUROC, and AUPRC [48]. The AUPRC, G-mean, and F₁-score were major metrics over AUROC due to the imbalanced data in this study [49,50]. The study also prioritized recall over precision as the major performance index, to minimize the cost of failing to detect AITD [51], while the accuracy; Precision, and AUROC were minor indices. The formulas of each evaluation metric are provided in Multimedia Appendix 6.

This study further analyzed individualized feature importance and survival curves of different thresholds to assess risk factors and differentiate high-risk patients. The Shapley additive explanations (SHAP) python package was used to understand the importance and influence of each risk factor that caused AITD [52]. The contribution of each feature was computed and plotted to interpret the model. The precision-recall (PR) curve of the best model was plotted to determine the optimal cutoff based on the maximum F₁-score. A threshold-moving system was further used by placing different cutoff points on the PR curve for binary classification. Five cutoff points were selected for analysis, including the points to forecast the top 1%, 5%, 15%, and 25% of patients with AITD risk, as well as the one determined by the threshold with the maximized F₁-score [53]. The recall and sensitivity for the above 5 cutoff points were then calculated and compared. A Kaplan-Meier (KM) survival curve was plotted using different cutoff thresholds to compare the actual survival of high- and low-risk groups for statistical comparison.

Baseline characteristics were evaluated with the chi-square test or Fisher exact test for categorical variables, and independent 2-tailed t tests were used for continuous variables. The Wilcoxon rank-sum test was used when the data were not normally distributed. The cumulative thyroid dysfunction incidence was compared with the log-rank test. Data were analyzed using SAS (version 9.4; SAS Institute); Python (version 3.9.5; Python Software Foundation), and R studio (version 1.3.1093; R Studio). The statistical significance of the AUPRC was calculated using MedCalc (MedCalc Software).

The study included 6497 amiodarone users. The results of a univariate analysis of their demographics and other features are shown in Table 1. The Strengthening The Reporting of Observational Studies in Epidemiology (STROBE) flowchart for patient selection is presented in Multimedia Appendix 7. The training set contained 4075 subjects, among whom 583 (14.3%) developed AITD, while the test set had 2422 patients, among whom 275 (11.35%) had AITD. The distribution of gender and mean age did not significantly differ between the training and test sets. The AITD group had a higher proportion of female patients, and its median age was older than that of the non-AITD group.

Feature selection by RFE with 5-fold cross-validation generated 19 features with an accuracy of 0.895 with AdaBoost, while 46 features with an accuracy of 0.914 were generated by XGBoost. Considering the simplicity and accuracy of the model, the 19 features selected by AdaBoost were used for further model development. A figure showing RFE and the features selected is provided in Multimedia Appendix 8.

Figure 2 compares the major model performance indices for the test set: AUPRC, F₁-score, G-mean, and recall. The internal validation performance of the training set is provided in Multimedia Appendix 9. The 4 major performance metrics for the XGBoost and AdaBoost models were consistently higher than for the KNN and LR models, with higher AUPRC and F₁-scores for the XGBoost-based model. Among different resampling methods, the best AUPRC was for XGBoost-B-SMT (0.751, 95% CI 0.697-0.799), which was significantly higher than for XGBoost-Raw (0.742, 95% CI 0.687-0.790; P<.05), XGBoost-ENN (0.741, 95% CI 0.686-0.790; P<.05), and XGBoost-B-SMT-ENN (0.730, 95% CI 0.675-0.779; P<.05). Multimedia Appendix 10 summarizes the results of the statistical comparisons for AUPRCs.

XGBoost-based models also produced higher G-means, with values of 0.688, 0.661, and 0.641 for XGBoost-B-SMT, XGBoost-ENN, and XGBoost-B-SMT-ENN, respectively. The 3 resampling methods all increased G-mean performance for XGBoost and AdaBoost but did not consistently increase the G-mean for KNN or LR. Finally, the recall levels of XGBoost-B-SMT and AdaBoost-B-SMT reached 0.756 and 0.858, respectively, while KNN-ENN only rounded to 0.255.

Table 2 further lists the performance on all metrics, including the minor indices of accuracy; Precision, and AUROC, for the test sets. All models had an AUROC >0.8, except the LR model without resampling. Like the major metrics, the AdaBoost- and XGBoost-based models had higher accuracy and AUROC values compared to the KNN- and LR-based models. The AUROC values for the XGBoost- and AdaBoost-based models were >0.9, while only the accuracy of the XGBoost-based models exceeded 0.9.

Figure 3 shows the SHAP summary plot. As shown in this graph, the TSH level had the highest contribution to AITD risk, with a SHAP value of 1.68. The fT4 level, amiodarone treatment duration, alkaline phosphatase level, high-density lipoprotein (HDL) level, and low-density lipoprotein (LDL) level were associated with a higher predicted probability of AITD, with respective SHAP values of 0.95, 0.76, 0.73, 0.52 and 0.37. Furthermore, therapeutic days, cumulative dose, age, and BMI, with SHAP values of 0.45, 0.41, 0.33, and 0.25, respectively, were important global predictors. The local explanation summary plot demonstrated the direction of relationships between clinical variables and AITD, with positive SHAP values indicating higher AITD risk. A higher TSH level was the most informative feature in determining AITD, with a lower fT4 level, shorter treatment duration, lower alkaline phosphatase level, higher HDL level, and lower LDL level increasing the AITD risk. In addition, longer therapeutic days, higher cumulative dose, older age, and lower BMI also raised the AITD risk.

Figure 4 shows KM curves based on using different cutoff points on the PR curve for the XGBoost–borderline SMOTE model. The thresholds of the 5 cutoff points were 0.995 for the top 1% (point A), 0.953 for the top 5% high-risk patients (point B), 0.627 for the optimal point determined by the maximized F₁-score (point C), 0.5 for the top 15% as the default value (point D), and 0.142 for the top 25% (point E) of patients predicted to be at risk of AITD. Moving from default point D, with a threshold of 0.5, to point A, with a threshold of 0.995, the recall significantly decreased from 0.756 to 0.116. When changing the threshold from 0.5 (default; Point D) to 0.142 (point E), recall increased from 0.756 to 0.88. The optimal cutoff (point C), with a threshold of 0.627, yielded better performance, with an accuracy of 0.93 and a precision of 0.685, and the best F₁-score was achieved at 0.699. The corresponding KM curves for the 5 cutoff points were able to differentiate high- and low-risk patients in the log-rank test (P<.001). Point A had only 33 within 2422 patients (1.4%) in the high-risk group, and point E predicted 617 within 2422 patients (25.4%) at high risk of AITD. Point C, with 286 high-risk patients within 2422 patients (11.8%), demonstrated an optimal prediction of 275 true positives for patients with AITD.

Figure 5 visualizes the prediction distribution of AITD for the 2422 subjects in the external test set. A color change from blue to red indicates increased predicted risk, with a dramatic change occurring around the cutoff point of 0.627. Among 275 patients who developed AITD, 196 (71.3%) had a predicted risk above 0.627 and were determined to be at high risk of AITD according to the XGBoost-borderline SMOTE model. Among the remaining 2147 non-AITD patients, 2057 (95.8%) had a predicted risk lower than the optimal threshold and were thus true negatives estimated by the model.

This study constructed an explainable and threshold-modifiable machine learning model with the resampling method for AITD risk stratification using dynamic clinical features from the Taipei Medical University clinical research database. The model with the best prediction performance, XGBoost-borderline SMOTE, was validated using external data from another hospital to ensure the credibility and generalizability of the results. It remained robust under conditions of different physicians; Prescription patterns, and hospitals. Resampling methods effectively tackled the imbalanced data and enhanced the model performance. There were 19 clinical features selected by the RFE. Time-series input for dynamic clinical features allowed for real-time assessment and prediction according to a patient’s changing disease state. The SHAP plot provided a better visualization tool to understand the contributions of features to AITD. Modifying the threshold on the PR curve by comparison to the KM curve could improve the help provided for clinical decision-making by determining the percentage of the AITD risk population in different practice settings.

The outstanding performance of XGBoost-borderline SMOTE in this study resulted from ensemble learning boosting algorithms and a resampling-oversampling technique. As a tree-based ensemble learning algorithm, XGBoost has been shown to be able to detect the complex and potentially nonlinear relationships in imbalance classification, such as in predicting adverse outcomes of chronic heart failure or the adverse effects of analgesics for osteoarthritis [37,54]. In contrast, KNN was sensitive to the majority of instances and thus performed poorly for imbalance classification [55], while the traditional LR model led to biased parameter estimates and classification performance and was less suitable for handling imbalanced data [55-57]. Interestingly, the study also found that oversampling B-SMT consistently outperformed undersampling ENN and hybrid sampling B-SMT-ENN on AUPRC, recall, and G-mean. The undersampling ENN and hybrid sampling methods used a process that deletes samples from the majority class; therefore, valuable information determining the decision boundary between the classes might have been lost. Borderline SMOTE created synthetic data only along the decision boundary, which was the best-performing strategy in this study; it has previously been used to successfully predict chronic kidney disease and hepatitis B virus infection with class imbalance [58,59].

This study used SHAP for model interpretation, which allowed precise ranking of variables with clinical reasoning and justification. Top features, such as a high TSH level, short treatment duration, low FT4 level, advanced age, and a higher cumulative amiodarone dose, have been well explained by previous epidemiological studies [6,7,16]. Low alkaline phosphatase and high alanine aminotransferase levels were also found to be associated with hypothyroidism [60]. However, HDL and LDL, which were selected by RFE in this research, were not statistically significantly correlated with thyroid function in previous pharmacoepidemiological reports [61]. Their contributions might be masked by factorial interactions in statistical approaches. Similar phenomena were reported for thyroid disease, chronic kidney disease, and analgesic adverse-effect prediction [54,62,63]. The machine learning models identified relevant features with nonlinear relationships and complex interactions between factors and outcomes, such as HDL and LDL, in the present study; Promising to help doctors and pharmacists to pay special attention when checking amiodarone users’ lipid panels to prevent adverse thyroid events.

This study used a moving threshold system on the PR curve to select optimal cutoff points for assessments with the KM curve. This innovative approach not only allowed comparisons of model performance at different decision thresholds, but also further ensured the capability of the model to differentiate clinically significant high- and low-risk groups. As decision threshold adjustment is a known strategy to deal with imbalanced classification, the best cutoff was based on the maximum F₁-score, as in a previous study of diabetes risk prediction [64]. Selecting an extremely low threshold allows capture of all potential AITD events, but a high false alarm rate can overwhelm clinicians. Conversely, an extremely high decision threshold can greatly reduce the false alarm rate with the cost of failing to detect AITD cases. With an unequal class distribution and high misclassification costs in adverse effect predictions, this study will increase attention paid to future studies to determine the optimal threshold based on the maximum F₁-score with a threshold-moving approach and a KM curve to ensure differentiation ability and clinical justification, rather than using the default cutoff (0.5) directly provided by the machine learning software.

This machine learning model could be used as a clinical decision-making aid for the early and real-time prediction of AITD by incorporating it into computerized physician order entry systems to optimize amiodarone use. This study collected patients’ time-series data to build the model, giving it the capability to provide assessment not only of new amiodarone users but also patients who have used amiodarone for a long time. Patients’ disease status changes over time, so only dynamically analyzing the cumulative dose, duration of therapy, and changes in multiple recent laboratory data will enable simultaneous surveillance. The present model does not provide a one-time glimpse of patient status, but a long-term, real-time prediction. Previous studies that used time-series concepts for intradialytic adverse-event risk prediction also provide evidence for this methodology, as their models had better prediction performance than models lacking features extracted from time-series data [65-67]. Using a threshold-moving system on the PR curve allowed us to further visualize threshold adjustments, balancing high noise and the cost of missing cases for clinical consideration. The model in this study should increase the safety of amiodarone use by enabling individualized risk prediction of AITD.

The study used the Taipei Medical University clinical research database, incorporating data from 3 hospitals in Taiwan. This database provided detailed clinical information, such as laboratory results, cause of death, and time-points for each medical treatment; it has previously been used to successfully predict mortality or classify patients with end-stage liver disease [68]. However, the study was potentially affected by the loss to follow-up of patients, unrecorded disease status, or unrecorded medications due to its nature as a retrospective data analysis. Family history, genetic data, and dietary intake were not documented in the database, but these factors might be integral to the occurrence of AITD. Extrapolation of the study model is thus restricted. Future multicenter, multicountry data are needed to further train and test the model before applying it in a broader clinical setting.

This study found that XGBoost with the borderline SMOTE resampling technique achieved the best model performance to predict AITD among amiodarone users. Feature selection by RFE and interpretation by SHAP demonstrated good predictive abilities and an explainable model. The optimal point of the threshold was determined to be the one with the maximal F₁-score, found by moving the threshold on the PR curve and differentiating risk groups assessed by the Kaplan-Meier curve. This time-series predictive model can serve as a preliminary tool to support clinicians with individualized AITD risk stratification among amiodarone users.