The Northeast China Rural Cardiovascular Health Study (NCRCHS) is a prospective cohort study conducted in the rural areas of Northeast China, whose inclusion criteria and design are described in the previous article (Guo et al., 2021a). In all, 11,956 participants aged ≥35 years were recruited from Dawa, Zhangwu, and Liaoyang counties of Liaoning province between 2012–2013 using a multi-stage, randomly stratified cluster-sampling scheme. The participants were followed up in 2015 and 2017. The NCRCHS was approved by the Ethics Committee of China Medical University (Shenyang, China). All participants provided written informed consent.

In this study, we included data from 5,260 participants of the NCRCHS study. Of these, 63 patients with missing or abnormal values were excluded. As a result, 5,197 hypertensive patients were included in the final analyses (Supplementary Figure S1).

The data collection procedure is described in detail in our previous paper (Li Z. et al., 2021). The patients' demographic characteristics, medical history and medication use, lifestyle factors, and other information were obtained at baseline through interviews by trained research staff using a standardized validated questionnaire. Indices such as weight, height, and waist circumference (WC). Body mass index (BMI) was computed as body weight in kilograms divided by the square of the height in meters. Blood pressure (BP) was assessed thrice after a 5-min rest using an automatic electronic sphygmomanometer (HEM-907; Omron, Kyoto, Japan) and averaged. HTN was defined according to the JNC 7 report as systolic blood pressure (SBP) ≥140 mmHg, diastolic blood pressure (DBP) ≥90 mmHg, and/or the use of antihypertensive medications (Chobanian et al., 2003). We collected blood samples after at least 12 H of fasting to determine the plasma levels of fasting glucose (FPG), triglycerides (TG), high-density lipoprotein cholesterol (HDL-C), uric acid, estimated glomerular filtration rate (eGFR), and blood routine biochemical indicators. Standard 12-lead electrocardiograms (ECGs) were used with a MAC 5500 (GE Healthcare, Little Chalfont, UK) as previously described (Li Z. et al., 2021). The ECGs were analyzed automatically, including QRS duration, PR duration, P axis, R axis, T axis, left ventricular hypertrophy (LVH) ECG (define per Sokolow–Lyon criteria), and QT interval (Framingham). Atrial fibrillation (AF) was defined as having a previous history of AF or an ECG suggestive of AF. Echocardiography was performed for all participants based on the American Society of Echocardiography guidelines, which were consistent with our previous study (Li T. et al., 2021). A Doppler echocardiography (Vivid; GE Healthcare, Connecticut, USA) with a 3.0-MHz transducer (Vivid, GE Healthcare, USA), including M-mode, two-dimensional, spectral, and color Doppler was used. Aortic dimension (AD), left atrial diameter (LAD), left ventricular end-diastolic internal dimension (LVIDd), left ventricular end—systolic internal dimension (LVIDs), interventricular septal thickness (IVSd), posterior wall thickness (PWTd), left ventricular ejection fraction (LVEF), E wave, and A wave were measured. In all, we selected 68 variables for model construction (Table 1), including subsets of characteristics related to the course, prognosis, hypertension-related target organ damage, and complications.

During follow-up, we collected the end point events of new fatal or non-fatal strokes. According to the World Health Organization (WHO) Multinational Monitoring of Trends and Determinants in CVD criteria, stroke was defined as rapidly developing signs of focal or global disturbance of cerebral functions lasting for >24 H (unless interrupted by surgery or death) with no apparent non-vascular cause (Asplund et al., 1988). Chronic cerebral vascular disease and transient ischemic attack were excluded. We collected medical records and death certificates for all participants who were possibly diagnosed or died. All information was independently reviewed and judged by the end point assessment committee.

We established four Cox regression models and compared them with the ML model. For one Cox proportional hazards model, Cox proportional hazards analysis was performed on all variables, and redundant variables were eliminated via the forward conditional stepwise selection method in Cox regression. Two regression models selected variables established in the China-PAR Score (Yang et al., 2016) and the Framingham Stroke Risk Score (D'Agostino et al., 1994). The last Cox model used LASSO regression to filter the variables.

To avoid data leakage, all data preprocessing and model construction procedures were conducted strictly within each training subset of a 10-fold cross-validation framework, while the corresponding validation subset was used only for model evaluation. Missing data were handled using multivariate imputation by chained equations (MICE), and outlier detection procedures were performed on the training data in each fold, with the derived parameters applied to the validation data.

The machine learning workflow, as shown in Supplementary Figure S1, involved feature selection using RFE, model training, and testing with XGBoost within a 10-fold cross-validation cycle. RFE was performed independently in each training subset, and validation errors across all folds were calculated to identify the optimal feature combination with the lowest average error. After feature selection, Synthetic Minority Over-sampling Technique (SMOTE) was applied only to the training data to balance class distribution, whereas validation data retained their original distribution, and SMOTE was implemented using the k-nearest neighbors approach (k = 5). Subsequently, XGBoost was used to develop the prediction model, a gradient tree boosting-based classifier that aggregates multiple weak learners into a strong learner, and was trained using the gbtree booster (max_depth = 4, learning rate = 0.05, n_estimators = 300, subsample = 0.8, and colsample_bytree = 0.8). All machine learning analyses were implemented in the open-source R software (version 4.1.1).

Missing laboratory values were imputed using the mice package in R (m = 5 imputations, max it = 10 iterations), with the imputation model fitted on the training data in each fold and then applied to the corresponding validation data. For outliers, the local outlier factor (LOF) was used for numerical variables and the attribute value frequency (AVF) algorithm for categorical variables. Continuous variables were represented as mean ± standard deviation (SD) and compared using the t-test or Mann–Whitney U test. Categorical variables were represented as frequency (n) and proportion (%) and compared using the chi-square test. C-statistics was used to evaluate the performance of the models (DeLong et al., 1988). Calibration of the models was evaluated by the Brier score method (range, 0– 1) (Brier, 1950) and the numbers of observed and predicted events proportion were grouped according to the decile of predicted risk (Liu et al., 2004). A P-value < 0.05 was considered statistically significant. Decision curve analysis (DCA) was performed to assess the clinical net benefit of each model at different threshold probabilities, providing a visual comparison of their potential clinical utility.

SHapley Additive exPlanations (SHAP) is a framework based on the additive feature attribution method that explains the output of the XGBoost model. A positive SHAP value indicates that the feature has a positive effect, while a negative SHAP value indicates that the feature reduces the outcome value and has a negative effect. This method can output the importance ranking of the features as well as the relationship between the features and the outcomes. SHAP-force plot was used to visualize the impact of individual feature values on the model's prediction for each observation. Descriptive analyses and comparisons between clinically defined groups were performed using R 4.1.1.

In all, 5,197 individuals with hypertension were included in the study. Of these, 49.4% were males, 50.6% were females, and the mean age was 57.16 years. End point events occurred in 294 (5.7%) patients during a mean follow-up period of 4.26 ± 1.03 years, among which 185 were males and 109 were females. Table 2 shows the distribution of the risk factors. Individuals with end points were older and had higher WC, SBP, DBP, FBG, QTc Framingham, AD, LAD, IVSd, LVIDs, PWTd, and A wave compared with those without end points. Furthermore, they take antihypertensive medications more frequently and were more likely to have a stroke history and left ventricular hypertrophy (LVH) ECG. In contrast, eGFR, R axis, LVEF, and E wave were lower in individuals with end points than those without end points.

Finally, we selected and analyzed 10 variables for the ML model by the RFE and XGBoost combination. The 10-variable combination included PWTd, IVSd, age, SBP, eGFR, QTc Framingham, platelet, calcium, uric acid, and FBG. We constructed a SHAP summary plot of the XGBoost model (Figure 1) to identify the importance of each feature in the prediction model. We identified that PWTd, IVSd, age, and SBP were the most important risk factors for stroke (Figure 1). In contrast, eGFR was associated with a decreased risk of stroke. We further generated a SHAP force-style plot for a representative patient to illustrate the individualized prediction of the XGBoost model, showing how each feature contributed to the predicted stroke recurrence risk (Supplementary Figure S2).

The results of the performance assessment were compared using C-statistics for discrimination and Brier score for calibration (Table 3) with 10-fold cycle of cross-validation. For predicting the end point events, the C-statistics was highest for the ML model [0.967 (95% CI, 0.956, 0.978)] among the five models. The four Cox regression models were similar to each other and their C-statistics were: Cox Regression (Framingham Stroke Risk Score) [0.747 (95% CI, 0.734, 0.762)], Cox Regression (China-PAR score) [0.725 (95% CI, 0.707, 0.731)], Cox Regression (Stepwise) [0.781 (95% CI, 0.772, 0.785)], and Cox Regression (LASSO) [0.764 (95% CI, 0.757, 0.772)]. Therefore, our ML model based on XGBoost had a better classification effect and better performance compared with the traditional risk scales.

The Brier score for the ML model was 0.053, indicating a good calibration between the estimated predicted risk and the observed 4.26 years risk. Calibration was also assessed by comparing the predicted and observed risks in each decile (Table 3). The largest difference for ML was small (2.9% in the 6th decile) compared with that of Cox (Framingham Stroke Risk Score) (42.7% in the 6th decile), Cox (China-PAR score) (33.8% in the 6th decile), Cox (stepwise) (27.9% in the 8th decile), and Cox (LASSO) (32.9% in the 8th decile).

Decision curve analysis (DCA) showed that all five models provided net clinical benefit mainly at low threshold probabilities (approximately 0–0.15), whereas net benefit rapidly approached zero as the threshold increased (0.15–0.20), indicating limited utility at moderate-to-high thresholds (Figure 3). The DCA curves largely overlapped across the evaluable range, suggesting that the differences in net benefit among models were small and not clinically apparent; importantly, the ML model demonstrated net benefit comparable to the Cox regression models with only marginal variations across thresholds.