Clinical data were retrospectively obtained from the Taipei Medical University Clinical Research Database (TMUCRD) in a deidentified format to ensure patient privacy. This database integrates electronic medical records from 3 affiliated hospitals within the Taipei Medical University system: Taipei Medical University Hospital, Taipei Municipal Wanfang Hospital, and Shuang Ho Hospital (managed by the Ministry of Health and Welfare, Taiwan). The TMUCRD provides comprehensive longitudinal data, including structured clinical information (demographics, diagnostic records, pharmacological treatments, and laboratory results) and unstructured narratives (physician progress notes, radiology and pathology reports, and discharge summaries), which encompass clinical features that support risk prediction modeling [24]. Data spanning 2011 to 2022 were analyzed. The study framework consisted of the following phases: cohort identification, feature engineering, model development (including training, evaluation, and interpretation), and risk prediction.

Patients aged 18 years or older with newly diagnosed RA were identified from the TMUCRD using International Classification of Diseases, Ninth Revision, Clinical Modification (ICD-9-CM), codes 714.0 to 714.2 and International Classification of Diseases, Tenth Revision, Clinical Modification (ICD-10-CM), codes M05 to M06. To ensure the inclusion of incident RA cases, we applied a minimum 12-month washout period prior to the first RA diagnosis, requiring at least 12 months of RA-free observation after database entry. Individuals with fewer than 3 rheumatology outpatient visits were excluded to reduce diagnostic misclassification. The occurrence of subsequent incident MACEs during follow-up was identified using ICD-9-CM and ICD-10-CM codes for the following events: stroke (ICD-9-CM: 430-437.xx; ICD-10-CM: G45.x, H34.1, and I60.x-I69.x), acute myocardial infarction (ICD-9-CM: 410.xx; ICD-10-CM: I21.xx and I22.xx), heart failure (ICD-9-CM: 428.xx; ICD-10-CM: I11.0, I50.xx, and I97.1), and acute coronary syndrome and ischemic heart disease (ICD-9-CM: 410-414.xx; ICD-10-CM: I21-I24.xx). If patients had any MACEs documented before the RA diagnosis, they were excluded.

To provide a benchmark for comparison with ML-based survival models, a conventional Cox PH regression approach was used. The patient cohort was stratified according to MACE occurrence. Univariable and multivariable Cox PH regression analyses were performed to evaluate the associations between clinical features and MACE risk. Kaplan-Meier (KM) analysis was further conducted to illustrate the cumulative incidence of MACEs over time, and log-rank tests were applied to compare risk differences between groups [25].

Continuous variables were summarized as means and SDs, whereas categorical variables were presented as counts and percentages. Differences in clinical characteristics were assessed using chi-square tests for categorical variables and 2-tailed independent-sample t tests for continuous variables, with statistical significance defined as a P value of less than .05. All statistical analyses and computations were performed using Python (version 3.11.4) and SAS (version 9.4; SAS Institute Inc).

This study was approved by the Taipei Medical University Joint Institutional Review Board (approval number N202312006) and authorized by the Data Management Center and the Institutional Data Governance Committee of Taipei Medical University.

Figure 1 shows the cohort identification and study workflow. From approximately 4.2 million individuals in the TMUCRD, 9771 patients with RA were initially identified. After applying the exclusion criteria, 2461 patients were included in the final analysis, of whom 253 (10.3%) experienced MACEs. The overall mean age was 51.9 (SD 13.9) years, and most patients were female (1917/2461, 77.9%). The dataset was stratified according to MACE occurrence and randomly divided into an 80% training set and a 20% test set. No significant differences in variable distributions were observed between the 2 subsets, ensuring their comparability for model development and evaluation (Table 1).

For each model, candidate hyperparameters were tuned using 5-fold cross-validation with GridSearchCV on the training set. The configuration yielding the highest average C-index across validation folds was retained as the optimal setting, as summarized in Table 2. Each model was subsequently retrained on the full training set using these optimized hyperparameters and evaluated on the independent test set, with detailed results shown in Table 3. Among all models, RSF achieved the best performance during 5-fold cross-validation, with an average C-index of 0.8443 (SD 0.0373; 95% CI 0.8116-0.8771). Consistent with the cross-validation results, the RSF model maintained the highest performance on the independent test set (C-index=0.8771; IBS=0.0775).

The RSF model was applied, and permutation importance was evaluated to rank feature contributions, as shown in Multimedia Appendix 1. Creatinine, csDMARDs, CRP, ALT, age, and RF were identified as the most critical features for predicting MACE risk, demonstrating markedly higher importance than other variables. Figure 2 shows the SHAP summary plots illustrating the contribution of individual features to the prediction of MACEs. Creatinine, which had the largest mean absolute SHAP values, exhibited the greatest impact on the model output, followed by ALT, csDMARDs, age, CRP, RF, AST, and ESR (left panel). In contrast, the use of csDMARDs, bDMARDs, glucocorticoids, and lipid-lowering agents showed negative SHAP values, suggesting protective factors against MACEs (right panel).

For comparison with ML- and DL-based survival models, a traditional Cox PH regression was used. Multimedia Appendix 2 summarizes the clinical characteristics stratified by MACE occurrence. Patients with MACEs were older and showed a higher prevalence of hypertension, DM, hyperlipidemia, COPD, and ILD. Laboratory markers were also elevated in the MACE group, including RF, ESR, and creatinine. In contrast, patients without MACEs had longer treatment durations with bDMARDs and csDMARDs.

In the Cox PH model (Multimedia Appendix 3), univariable regression showed that an older age at RA diagnosis; comorbidities such as DM, hypertension, hyperlipidemia, COPD, and ILD; elevated ESR; and higher creatinine levels were associated with increased MACE risk . In contrast, the defined medication combination was associated with lower MACE risk. Multivariate regression further identified age at RA diagnosis, hypertension, and creatinine levels as independent predictors of MACEs. Among these, hypertension was the most prominent risk factor for MACEs.

For further model application and personalized interpretation, the KM survival curve for the entire cohort was overlaid with the RSF model’s average survival function based on the testing set, as shown in Figure 3, demonstrating highly consistent 5-year survival rates of 0.8725 for the KM model and 0.8549 for the RSF model.

To the best of our knowledge, this study is the first to use multi-hospital electronic health record (EHR) data to develop and validate interpretable time-to-event ML and DL models for predicting MACE risk in Asian patients with RA. Using real-world data from a large health care system, the RSF model achieved the best performance in the independent test set (C-index=0.8771; IBS=0.0775), outperforming the Cox PH, DeepSurv, and Cox-Time models. The strong concordance between the RSF-predicted and KM-estimated survival curves confirmed the robustness of the model. Moreover, by integrating the SHAP framework, the model provided transparent interpretability and identified key features, including creatinine, csDMARD use, CRP, ALT, and RA diagnosis age, as the most influential determinants of MACE risk in RA.

This study offers several important advances over previous prognostic research on RA. First, compared with conventional Cox models that typically rely on linear effects and PH assumptions, ML- and DL-based survival algorithms such as RSF, DeepSurv, and Cox-Time extend traditional Cox modeling by capturing nonlinear effects and, in some cases, relaxing the PH assumption. RSF and Cox-Time in particular provide more flexible modeling frameworks that may better accommodate the complexity of real-world clinical data. Second, through 5-fold cross-validation and independent test set evaluation, the RSF model achieved the best predictive performance among the models examined, supporting the feasibility of this modeling framework for MACE risk prediction in patients with RA. Third, we combined permutation importance, which quantified the overall contribution of each feature to model performance, with SHAP analysis, which enhanced individual-level interpretability, thereby linking model predictions to personalized risk assessment.

Analysis of the RSF model using SHAP interpretability techniques identified creatinine, ALT, CRP, csDMARD use, and RA diagnosis age as key features associated with MACE risk. The observed associations were largely consistent with previous evidence: elevated creatinine reflects the adverse cardiovascular impact of renal dysfunction [39,40]. Higher ALT levels may reflect fatty liver disease [41], including nonalcoholic fatty liver disease and metabolic dysfunction–associated fatty liver disease, both of which have been associated with an increased risk of MACEs [42-44]. Similarly, higher CRP highlights the critical role of systemic inflammation in RA-related cardiovascular complications [45,46], and an older age at RA diagnosis has been linked to an increased risk of adverse cardiovascular outcomes [7,47]. Moreover, SHAP analysis supported prior evidence by showing that longer exposure to csDMARDs, bDMARDs [48,49], and lipid-lowering agents [50,51] aligns with reports highlighting the prognostic importance for cardiovascular outcomes of initiating effective treatment and inflammation control, as illustrated by the red-colored dots clustering toward the left (negative SHAP values) in Figure 2 (right panel).

Interestingly, previous studies have reported conflicting results regarding the association between glucocorticoids and MACEs. Cardiovascular protection appears limited to low-dose use (≤5 mg of prednisone daily) [52], whereas higher doses, prolonged use, older age, and comorbidities are associated with increased risk [53]. In SHAP analysis, glucocorticoids were linked to reduced MACE risk. Additional analysis of follow-up duration and exposure days in this study revealed that MACE cases had brief daily exposure, suggesting short-term, low-dose use. Nonetheless, MACE risk should be interpreted in relation to all clinical features, as reflected by the overall distribution of SHAP values.

In summary, our study, by integrating the RSF model with permutation importance and SHAP analyses, offers valuable insights into high-risk factors and provides interpretable predictions for clinical application. This demonstrates the feasibility of ML models for predicting MACE risk in patients with RA, potentially aiding case managers in improving patient prognosis and developing personalized treatment plans.

Our study offers several significant advances over existing prognostic models. By leveraging large-scale, multi-hospital EHR data to develop survival ML models for MACE risk in patients with RA, this study demonstrates superior predictive accuracy through the RSF model’s ability to capture complex nonlinear relationships and complex feature interactions without relying on the PH assumption. By integrating SHAP analysis, it provides transparent, individual-level risk interpretations that bridge the gap between complex algorithms and personalized clinical decision-making.

Despite this study’s strengths, several limitations should be considered. First, this retrospective study involved substantial missingness in the TMUCRD dataset. Variables with more than 50% missing data were excluded as such missingness often suggests limited clinical relevance in routine care. For the remaining important variables, missing values were addressed using an iterative chained equations imputation approach informed by the multiple imputation by chained equations to preserve data structure and intervariable relationships. Although this may reduce bias compared with single imputation, excluding highly incomplete variables may still have limited our ability to capture patient heterogeneity. Missing laboratory covariates were further imputed under a missing at random assumption; however, the imputation model did not incorporate survival outcome information, which may have affected covariate-outcome associations. Second, the model was internally evaluated using a stratified train-test split and cross-validation, but neither external nor temporal validation was performed in this study. Accordingly, the generalizability of the findings to other hospitals, health care settings, ethnic populations, and calendar periods remains uncertain. In addition, calibration plots, decision curve analysis, and threshold-based performance evaluation were not included in this study. Therefore, the findings should be interpreted as preliminary rather than definitive evidence of clinical usefulness. Third, although ML-based survival models were applied, the predictors were incorporated as patient-level covariates rather than as truly time-varying covariates updated during follow-up. Therefore, the models may not fully capture dynamic changes in clinical status, laboratory results, and treatment exposure over time. Fourth, despite the use of cross-validation, the relatively modest number of events (n=253) and the complexity of the ML models may still introduce a risk of overfitting. Fifth, while disease and medication may contribute to MACE risk prediction, medication exposure in observational EHR data may be influenced by disease severity, prescribing patterns, and other potential confounders. Finally, the EHR database lacked key lifestyle covariates such as smoking status, BMI, and physical activity, which may have resulted in incomplete adjustment for potential confounding.

In conclusion, this study demonstrates that the RSF model significantly outperformed traditional Cox PH and DL-based survival models in predicting MACE risk among patients with RA. By effectively capturing nonlinear relationships and feature interactions, our model provides superior predictive accuracy and robustness. The integration of permutation importance and SHAP analyses bridges the gap between complex “black-box” algorithms and clinical practice, offering a transparent framework for personalized risk assessment. These findings empower clinicians and case managers to identify high-risk individuals early and develop tailored preventive strategies. Future work should include calibration assessment, decision curve analysis, and threshold-based performance evaluation, as well as temporal validation using data from later calendar periods; external validation using data from other health care systems across diverse ethnic populations; and the incorporation of time-varying variables for clinical status, laboratory measurements, and treatment exposure to improve predictive accuracy, clinical utility, and generalizability.