This retrospective study used discharge data from a tertiary hospital in Northeast China. The data collection period was from January 1, 2021, to December 31, 2023. The dataset contains comprehensive information on inpatients, including demographic characteristics (such as age and gender), diagnostic information at admission, including primary and secondary diagnoses, admission and discharge status, and total hospitalization costs. On the basis of the ICD-10 coding system, this study screened for patients with stroke and ultimately included inpatient cases with a primary diagnosis code of I60, I61, or I63. To minimize information bias and ensure data quality, this study established a rigorous data screening and preprocessing workflow. During the data quality control phase, extreme cases with a LOS of less than 24 hours or more than 90 days were first excluded. Subsequently, a systematic evaluation of the missingness of key features was conducted. To balance data integrity and sample representativeness, as well as to prevent the introduction of potential systematic bias from overreliance on data imputation, variables demonstrating a missingness rate greater than 20% were eliminated [32]. Among the initially extracted raw features, a total of 2 variables were excluded for exceeding this threshold; the remaining variables exhibited good data completeness, with missing proportions ranging from 0.13% to 1.95%. Statistical evaluation confirmed that the missing data in this study satisfied the assumption of missing at random, and no statistically significant differences in missingness rates were detected across the various cost groups. To impute the remaining missing values, this research used the K-Nearest Neighbors (KNN) imputation technique, a method that estimates missing data by leveraging local structural similarities across samples, thereby effectively retaining the original data’s distributional properties and the underlying associations among variables. To avoid data leakage from distance-based imputation, we fitted the KNN imputer using only the 2021 to 2022 network construction cohort and the 2023 training set and then applied the trained fixed model to the 2023 test set for missing value imputation, which ensured complete isolation of data between the training and test sets.

As readmission within 30 days of discharge typically exhibits a strong clinical correlation with the initial hospitalization, it is considered to reflect continuous disease progression following the primary admission. Therefore, for patients with multiple hospitalization records, this study referred to the processing methods of previous health services research and merged the data using a 30-day time window rule based on unique inpatient numbers and admission and discharge dates [33,34]. The specific record integration rule was as follows: if the time interval between a patient’s previous discharge and subsequent admission was ≤30 days, they were merged into a single hospitalization event; if the time interval was greater than 30 days, they were treated as independent inpatient cases. For cases with multiple consecutive hospitalizations where the intervals between all adjacent admissions met the ≤30-day condition, this study adopted a continuous sequential merging strategy, unifying them into one complete medical event. When merging specific variables, the following principles were adhered to: for features such as gender, age, insurance type, admission year, and primary diagnosis, data from the first record within the merged period were extracted; the admission route was determined using a priority rule (ie, prioritizing emergency admissions, followed by outpatient, and finally other routes); regarding the planned readmission within 31 days indicator, if any record within the merged group was marked as “yes,” the entire merged event was classified as “yes”; and the discharge disposition was based on the status at the time of the final discharge. For secondary diagnostic codes, the system aggregated all records within the group and performed deduplication to construct a complete, comprehensive diagnostic set for the patient. Furthermore, the total hospitalization costs and LOS were calculated by summing all records within the merged group. In this merging process, a total of 396 patients accounted for 1096 repeat hospitalization records; through the aforementioned rules, 518 records were retained individually, and the remaining records were merged into 132, ultimately forming 650 hospitalization event records.

Following the aforementioned screening and processing, this study ultimately included 10,556 patients with stroke, of which 6618 were from 2021 to 2022 and 3938 were from 2023. To prevent data leakage, this study adopted a time span–based data partitioning strategy: data from 2021 to 2022 were used to construct the comorbidity network and extract topological features (the definition of high-cost nodes within the network was also determined entirely based on data from this period); meanwhile, data from 2023 were independently divided into training and testing sets. It should be specifically noted that the network features extracted on this basis are essentially supervised risk representations based on historical data, but no information related to the 2023 data was used during their derivation process. In addition, to eliminate measurement bias caused by macroeconomic fluctuations and inflation, this study uniformly adjusted the total hospitalization costs for 2021 and 2022 to the 2023 value level based on the Consumer Price Index. The complete research flowchart is shown in Figure 1.

This study was approved by the Scientific Research Ethics Committee of Hangzhou Normal University (approval 2025‐1026). All procedures strictly complied with the Declaration of Helsinki and relevant ethical guidelines. Data were collected retrospectively through medical record reviews. To ensure patient privacy, all data were deidentified prior to analysis.

Prior to network construction, this study established a minimum prevalence threshold, based exclusively on the 2021 to 2022 network construction cohort. Specifically, only diseases with ≥5 affected individuals within this cohort were included as network nodes, thereby ensuring the stability of the constructed network and the statistical representativeness of the extracted features. Additionally, this study incorporated the patients’ primary and secondary diagnosis codes in the early stage after admission into the network construction process. Considering the extremely high prevalence of the primary diagnosis within the study cohort, to avoid its potential computational bias on the assessment of disease associations, this study performed corresponding statistical corrections in the subsequent calculation formulas for feature indicators.

In this network, nodes represent the included disease categories, and edges represent the pairwise correlations between them. Previous studies have shown that the choice of association measure can significantly affect network topology and research findings [26]. Therefore, this study used a co-occurrence correlation as the indicator to quantify the strength of association between disease pairs [35]. Compared with traditional measurement methods, the core advantage of this calculation formula lies in its denominator, which uses the arithmetic square root of the sum of the squared occurrence frequencies of each disease. This characteristic makes the indicator more robust when dealing with extremely unbalanced occurrence frequencies of 2 diseases; even if a certain disease is relatively rare, as long as it has a very high co-occurrence ratio with another high-frequency disease, this formula can still objectively assign it a reasonable association weight. To filter out noisy edges caused by coincidental co-occurrences, this study set a minimum edge weight screening mechanism, retaining only associations with an edge weight >0.01. The final constructed PCN contained 139 valid nodes and 2404 valid edges, with a network density of 0.251. The specific calculation formula is defined as follows:

Here, C_xy refers to the co-occurrence frequency of diseases x and y within each patient, while P_x and P_y represent the prevalence rates of diseases x and y, respectively. Figure 2 illustrates the generated PCN; to optimize visualization, the network displays only edges with weights ≥0.1.

Before constructing the distance-based disease-cost network (DDCN), this study first built an initial network, specifically using the OR to quantify the strength of association between nodes. Unlike the aforementioned network, this network introduced a specific “high-cost” node in addition to the conventional disease nodes, aiming to further explore the association patterns between specific diseases and high medical expenditures. The calculation formula is defined as follows:

In this framework, E represents the specific high-cost outcome. D_E and D_N represent the number of high-cost and non–high-cost patients, respectively, in the patient subgroup with the specific comorbidity d. Accordingly, H_E and H_N, respectively, represent the number of patients corresponding to high-cost and non–high-cost states in the control subgroup (ie, patients without the comorbidity d). As disease co-occurrence relationships are symmetric, the constructed DDCN is an undirected graph. To address the issue of undefined OR values caused by zero cells in the 2×2 contingency table, this study introduced the Haldane-Anscombe correction, which adds 0.5 to all 4 cells of the contingency table to avoid division-by-zero errors and smooth extreme values generated by small samples. To ensure the high specificity of the network edges, the network only retained significant association edges with an OR >1 and a lower bound of its 95% CI >1. Subsequently, this study used the min-max normalization formula to transform the ORs into distance metrics, thereby completing the construction of the DDCN. Through this formula, all association strengths were mapped into the standardized distance space of [0, 1], effectively smoothing the extreme numerical fluctuations associated with the extremely high prevalence of the primary diagnosis. The finally constructed DDCN contained 140 valid nodes and 1713 valid edges, with a network density of 0.176.

In addition to the topological features extracted from the comorbidity network, this study also extracted comorbidity-related features similarly derived from patients’ diagnostic coding data to compare the impact of different feature subsets on the model’s identification efficacy in subsequent analyses. Specifically, these features include the number of comorbidities, the CCI, and the ECI.

This study extracted a total of 8 baseline features from the medical record front pages, specifically including gender, age, insurance type, stroke type, admission route, discharge disposition, LOS, and planned readmission within 31 days. Among these, the insurance type was divided into 4 categories based on the local medical insurance pooling management structure: municipal insurance (applicable to the insured population in the municipal pooling area where the medical institution is located), provincial insurance (corresponding to the insured population in the provincial pooling area), out-of-town insurance (referring to the situation where the insured location and the treatment location belong to different administrative regions), and others (covering various special forms of medical security). The stroke type was defined based on the primary diagnosis code: cases with codes I60 or I61 were classified as hemorrhagic stroke, and cases with code I63 were classified as ischemic stroke. The admission route was divided into emergency, outpatient, and others; the discharge disposition covered routine discharge upon medical advice, discharge against medical advice, and death. The planned readmission within 31 days was dichotomized as “yes” or “no” based on the discharge record in the hospital discharge data.

It is worth noting that among the 8 baseline features mentioned earlier, 5 features—gender, age, insurance type, stroke type, and admission way—can be obtained at the initial stage of patient admission, whereas discharge disposition, LOS, and planned readmission within 31 days can only be finally determined at the hospital discharge stage, exhibiting an obvious information lag. Accordingly, in the subsequent analysis, we excluded these 3 lagging features and additionally developed a prediction model specifically for patients in the early stage of hospitalization to investigate its performance in the early identification of high-cost patients.

Although previous literature has not reached a complete consensus on the definition of high-cost status, using the top 10th percentile as a classification benchmark has been widely recognized and adopted [36,37]. Consequently, this study defined cases ranking in the top 10% of total hospitalization costs within the study cohort as high-cost patients. Furthermore, the detailed components of inpatient costs were solely used to define the outcome variable and were strictly excluded from subsequent analyses. The selection of this threshold was mainly based on the following considerations: on the one hand, this criterion can effectively alleviate the extreme data imbalance problem associated with setting a threshold too strictly, thereby ensuring the efficacy of subsequent statistical tests and classification evaluations; on the other hand, this criterion can also prevent the dilution of the typical distribution characteristics of the high-cost group caused by a definition scope that is too broad. Furthermore, from the macro perspective of medical policy, the top 10% classification can accurately cover the high-multiplier case groups closely monitored under China’s current medical insurance payment system, highly aligning with the practical governance needs of current hospital cost control and refined medical insurance management. To further verify the robustness of the research findings, we conducted expanded testing in the sensitivity analysis using other alternative thresholds.

This study used the 2023 data subset for the construction and validation of models. To ensure the consistency of the target variable distribution between the training and testing sets, this study adopted the hold-out method, performing a stratified random splitting based on the outcome variable at a ratio of 8:2. After splitting, a training set of 3151 cases (including 315 high-cost cases) and an independent testing set of 787 cases (including 78 high-cost cases) were ultimately generated. Aiming at the class imbalance problem existing between high-cost and non–high-cost cases, this study strictly applied the Synthetic Minority Over-sampling Technique (SMOTE) only within the training set before formally training the models. To rigorously prevent overestimation of model performance due to data leakage, no global oversampling was performed prior to cross-validation. Instead, the SMOTE algorithm was specifically embedded within the resampling pipeline. Specifically, during each iteration of the 5-fold cross-validation, SMOTE was dynamically applied only to the training folds, whereas the validation fold always retained the original real-world data distribution. This strategy effectively balanced the class distribution, mitigated the classification bias of machine learning algorithms toward the majority class, and ensured the objectivity of model evaluation. On the basis of the processed data, this study constructed 5 machine learning models to identify high-cost stroke inpatient cases, specifically covering Decision Tree (DT), Support Vector Machine (SVM), Neural Network (NN), Random Forest (RF), and Extreme Gradient Boosting (XGBoost). To guarantee the reproducibility of the experiment and effectively control the risk of overfitting, the model’s hyperparameter optimization process used a grid search strategy combined with 5-fold cross-validation strictly executed within the training set. This was done to objectively screen and determine the optimal parameter combination for each model. During the hyperparameter tuning phase, specific search spaces were tailored to the architecture of each model. For DT, the cost-complexity parameter was primarily adjusted. For SVM, both the Radial Basis Function kernel parameter and the regularization penalty were optimized. For NN, the optimization focused on the number of hidden nodes and the weight decay. For RF, the number of candidate variables per split and the minimum node size were fine-tuned. Finally, for XGBoost, a simultaneous optimization was performed on core parameters, including the number of boosting iterations, maximum tree depth, learning rate, and minimum child weight. To balance computational efficiency and mitigate the risk of overfitting, the grid search was restricted to predefined boundaries. Consequently, the identified parameters represent local optima within these ranges rather than exhaustive global optima. The identification performance of all final models was objectively evaluated on the single-split independent hold-out testing set, combined with 1000 Bootstrap resampling to obtain the 95% CIs for various evaluation metrics. To comprehensively measure the overall identification ability of the models, this study adopted the following multidimensional evaluation metric matrix: the area under the curve (AUC), sensitivity, specificity, accuracy, G-mean index, and F₁-score [38]. Addressing the class imbalance characteristic commonly found in medical data, this study specifically introduced the G-mean and F₁-score to more accurately and robustly reflect the overall balanced performance of the model in simultaneously handling the identification of both minority and majority classes. Furthermore, to objectively evaluate the potential gain of network topological features on the model’s identification performance, this study used the DeLong test to conduct a statistical examination of the AUC differences between models containing only baseline features and composite models integrating both baseline and network features.

In addition to the aforementioned evaluation metrics, this study further introduced calibration curves and the Brier score, aiming to evaluate the consistency between the high-cost identification probabilities output by the model and the actual observed true proportions. Finally, this study used decision curve analysis (DCA) to evaluate the net clinical benefit that the model can generate across a wide range of threshold probabilities, thereby verifying its practical application value in medical management and decision-making. Furthermore, to deeply explore the incremental value of the extracted network topological features, this study also implemented an ablation comparative analysis to systematically compare the objective differences in identification performance between “models containing only baseline clinical features” and “models further integrating traditional comorbidity features or network features.”

To deeply analyze the key associated features behind the high medical expenditures of inpatients with stroke, this study introduced the SHAP technique to provide post-hoc interpretability analysis for the optimal-performing machine learning model. Given that the XGBoost model ultimately used for feature attribution is a tree-based ensemble algorithm, this study specifically selected the TreeExplainer optimized for such algorithms in the analysis to accurately calculate the local and global explanations of features. On the basis of co-operative game theory, the SHAP method can theoretically and fairly distribute the model’s identification results among all input variables. Under this framework, the global importance of a feature is quantified by the average level of its absolute SHAP values across the entire study cohort. Specifically, a positive SHAP value indicates that the feature is associated with an increased probability of a case being classified as “high cost,” while a negative SHAP value denotes a decrease in this classification probability.

This study used R (version 4.4.2; R Core Team) for related statistical analyses and plotting. In the comparative analysis of features between groups, the distribution characteristics of categorical variables were expressed as frequencies and percentages, and the differences between the high-cost and non–high-cost groups were compared using the chi-square test; continuous variables underwent a normality test prior to analysis, where data conforming to a normal distribution were presented as mean (SD), and intergroup comparisons were conducted using the independent samples t test; continuous variables not conforming to a normal distribution were represented by the median and IQR, and intergroup comparisons adopted the Mann-Whitney U test. The visualization of the PCN was completed in Cytoscape (Cytoscape Consortium).

According to the previously defined high-cost threshold criteria (the top 10%, equating to ≥94,696.79 RMB), of the 3938 inpatient with stroke incorporated into the 2023 analysis, 393 cases were assigned to the high-cost cohort, and the remaining 3545 cases comprised the non–high-cost cohort. Figure 3 displays the distribution of inpatient costs among patients with stroke, with the total health care expenditures of the high-cost group making up 55.1% of the overarching cohort’s total costs.

The results of the intergroup analysis presented in Table 1 reveal that the high-cost and non–high-cost groups exhibit statistically significant differences across various features. Specifically, regarding demographics, the proportion of female patients in the high-cost group was significantly higher than that in the non–high-cost group; in terms of age distribution, the proportions of patients aged 51 to 60 years and 61 to 70 years in the high-cost group were slightly higher than those in the non–high-cost group, whereas the proportion of the population aged over 71 years in the high-cost group was distinctly lower than in the non–high-cost group. In terms of clinical features, stroke types exhibited substantial intergroup differences. Ischemic stroke predominated in the overall population, while hemorrhagic stroke accounted for 57% of the high-cost group, a proportion significantly higher than the 14.9% observed in the non–high-cost group. Furthermore, the LOS showed a strong positive correlation with high-cost status. The proportion of patients with provincial insurance among high-cost cases was markedly higher than in the non–high-cost group. Regarding network features, the intergroup comparison revealed that the NHCP value was significantly elevated in the high-cost group; conversely, the short distance value in this group was significantly reduced. As for comorbidity features, the results indicated that the high-cost patient group had a higher proportion of lower CCI scores and fewer comorbidities compared to the non–high-cost group.

Table 2 presents the identification performance of the 5 machine learning models after incorporating the 2 network features. Overall, after incorporating the network-derived features, the identification performance of all models achieved a certain degree of improvement. Among them, XGBoost exhibited the largest AUC improvement, increasing from 0.814 to 0.899; moreover, the remaining evaluation metrics of this model also demonstrated the optimal performance among the 5 models.

Figure 4 displays the receiver operating characteristic curves of the 5 machine learning models combined with network features; the DeLong test results in Figure 4B further confirmed that the difference in AUC gain between the “composite model fusing network features” and the “initial model containing only baseline features” is statistically significant. Specifically, the difference in AUC between the 2 models was 0.085 (95% CI 0.052‐0.123), with an exact P value of 1.42×10⁻⁶ according to the DeLong test.

The calibration curve results in Figure 5A show that the curve trajectory of the XGBoost model is closest to the ideal reference line, indicating a high degree of consistency between the high-cost identification probabilities output by the model and the actual observed true proportions. Simultaneously, XGBoost achieved a Brier score of 0.087, belonging to the low-error tier alongside RF. Its overall calibration performance was significantly better than that of the SVM, DT, and NN, further confirming the robust reliability of its output identification probability estimates. In the DCA shown in Figure 5B, XGBoost similarly highlighted significant practical application advantages. Specifically, within the broad threshold probability interval of clinical and management significance, the net benefit curve of XGBoost consistently remained at the highest level; this means that compared to other machine learning models participating in the evaluation, using this model to guide medical resource allocation can yield a more substantial comprehensive net benefit. Given the outstanding comprehensive identification performance described earlier, this study ultimately selected XGBoost as the core model for subsequent feature attribution and mechanism analysis. The optimal hyperparameter combinations and their specific configurations for each model are detailed in Table A1 in the Multimedia Appendix 1.

To objectively evaluate the specific contributions of traditional comorbidity features and single network topological features to the model’s identification performance, this study further implemented an ablation analysis relying on the XGBoost model. The outcomes are presented in Table 3: in terms of traditional comorbidity features, regardless of the introduction of CCI, ECI, or a basic comorbidity count, the model’s identification performance demonstrated a certain degree of enhancement; however, its overall performance still fell markedly short of the model integrating network features. In contrast, the incorporation of CCI led to a certain degree of performance degradation. Furthermore, the 2 network topological features also had different focuses regarding their gains in identification performance. Specifically, the “baseline+SD” model demonstrated better discrimination in the AUC metric; in contrast, although the “baseline+NHCP” model had a slightly lower AUC value, it performed better on the F₁-score. Finally, by rigorously excluding lagging variables that could only be determined after discharge and retaining only features available in the early stage of hospitalization, we independently developed a corresponding early identification model. The results showed that, although the model achieved an acceptable AUC, its F₁-score was markedly low, suggesting a potential imbalance in model performance when based solely on the currently available features.

Previous intergroup analysis of baseline features showed a statistically significant difference in stroke type between the 2 groups, indicating a high degree of clinical relevance between this feature and high-cost status. Given that all diagnostic codes, including the principal diagnosis, were incorporated in the construction of the comorbidity network, the resulting network topological features to some extent captured information inherent to stroke subtypes. To deeply explore the potential multicollinearity and information overlap effects among features, this study additionally constructed a model that excluded the stroke subtype and retained only the network features. The results showed that after excluding this subtype feature, the model’s identification performance was basically unaffected (relevant evaluation data are detailed in Table A2 in the Multimedia Appendix 1). In other words, the value of network features lies more in complementing existing clinical information than in providing discriminative ability that is entirely independent of stroke subtypes. Therefore, the incremental value of network features should be interpreted as a marginal improvement beyond established clinical variables, rather than as a complete replacement for conventional features. Table A2 in Multimedia Appendix 1 also summarizes the performance comparison of models fusing all feature combinations. The comprehensive analytical results corroborated that the feature combination of “baseline characteristics+dual network features” delivered the most optimal overall identification capabilities across all evaluation metrics. Therefore, this study ultimately selected this feature subset as the standard input for the subsequent SHAP feature attribution analysis.

Meanwhile, to verify the applicability and robustness of the aforementioned identification framework in specific clinical subgroups, this study further evaluated the model’s performance within the hemorrhagic and ischemic stroke subgroups in the independent testing set, based on this optimal feature combination. The stratified analysis results (Table 4) showed that although the AUC value of each model within a single subgroup slightly decreased by approximately 0.03 to 0.04 compared to the overall cohort, they still demonstrated stable identification performance across different clinical subgroups.

To enhance the transparency and interpretability of the model in the context of clinical decision-making, this study introduced the SHAP method to objectively quantify the marginal contribution of each input variable to the final identification result. This method covers 2 dimensions: global interpretability at the feature level and local interpretability at the individual level. In terms of global interpretability, the bar chart in Figure 6A is sorted in descending order based on the mean absolute SHAP values of the features, intuitively illustrating the relative importance of different variables in the identification of high-cost patients. The analysis results showed that the top 5 core associated features contributing to the model’s performance were, in order, short distance, LOS, NHCP, age, and insurance type. As shown in the pie chart in Figure 6A, the cumulative contribution of network features accounts for 50.4% of the model’s total output. Among the baseline features, baseline features available at admission account for 19.7%, while baseline features available at hospital discharge reach 29.9%. The SHAP beeswarm plot in Figure 6B intuitively presents the distribution trend of SHAP values across various features for the entire sample. Figure 6C displays the SHAP dependence plots for the top 3 features, showing how these variables regulate the model’s output within different value ranges. Specifically, a lower value of shortest distance, along with higher values of LOS and NHCP, is significantly associated with a high probability of cases being identified as a high-cost status. The SHAP interaction analysis in Figure 6D reveals that the interaction effects between network features and LOS exhibit complex nonlinear synergistic associations. Specifically, patients with a LOS ≥14 days show a significant inverted U-shaped trend, with the interaction value reaching a positive peak when shortest distance is at a moderate level, significantly elevating the high-cost risk; the 11‐ to 13-day group displays a U-shape with mostly negative interactions; whereas the ≤10-day group has a flat interaction effect, fluctuating around the baseline. On the other hand, when examining the interaction between NHCP and LOS: for patients with an LOS of 8 to 10 days, the interaction effect shows an upward trend, turning from a negative value to a positive value as NHCP increases. Conversely, for the subgroup with an LOS ≥14 days, the interaction effect presents a downward trend, turning negative with the increase of NHCP. In contrast, for the subgroups with an LOS ≤7 days and 11 to 13 days, the interaction values remain relatively stable across different levels of NHCP.

At the local explanation level, we generated SHAP force plots to elucidate the contribution of each feature to the identification result for a specific individual. Figure 7A displays a patient who was identified by the model as highly likely to be in a high-cost status. The core features making the primary positive contribution to this identification probability were shortest distance and NHCP, while the LOS played a certain mitigating role. Conversely, the patient in Figure 7B was identified as having an extremely low probability of being in a high-cost status. The dominant associated features supporting the identification of this case into the non–high-cost group were, in order, a higher shortest distance, a shorter LOS, a lower NHCP, and a specific insurance type.

As presented in Table 5, at the extreme threshold of 5%, the model’s AUC value reached its highest level; however, affected by the extremely imbalanced data, the model’s F₁-score was below 0.4, meaning that the model performance was poor and the identification ability was relatively weak. As the threshold was relaxed, although the model’s AUC value experienced a slight decline, it showed a relatively stable trend within the 10% to 20% range. In addition, due to the steady increase in the number of positive samples, the class imbalance phenomenon was effectively alleviated, prompting the F₁-scores of all models to show a gradually increasing trend. Among them, the threshold adjustment from 5% to 10% optimized the model performance most significantly, and the F₁-scores of all models achieved substantial improvements. This indicates that the 10% threshold effectively improved the problem of impaired model performance caused by extreme class imbalance under the low threshold.

Unlike most previous models, which incorporated clinical laboratory, imaging, or genomic data [39-41], this study used hospital discharge data to develop a high-performance model for identifying patients with stroke who have high hospitalization costs by combining diagnostic network analysis with machine learning algorithms. Compared with these data sources, hospital discharge data offer several important advantages, including greater standardization, lower cost, and broad availability across hospitals at all levels, making it a practical basis for implementation in diverse health care settings. Although this data source has inherent limitations in dimensionality, incorporating topological features derived from the comorbidity network substantially improved the performance of all models, consistent with previous studies [30,31]. Moreover, compared with other models for identifying high-need, high-cost patients that did not include network features, our model showed superior performance [17,42-44].

Our results indicate that, although both are derived from diagnostic data, network features from the comorbidity network contributed more to model performance than conventional comorbidity measures such as CCI, ECI, and comorbidity count. By quantifying the pathway from a specific disease node to high-cost outcomes and integrating the disease’s intrinsic risk with the combined effects of its local neighborhood, these features captured latent information that is difficult to detect using traditional statistical methods, thereby substantially improving the model’s ability to identify high-cost patients. Additionally, subgroup analyses across different stroke subtypes further verify the stability and reliability of the comorbidity network features extracted in this study. The model maintained satisfactory predictive performance under the macroscopic classification of hemorrhagic and ischemic stroke, demonstrating good applicability across major clinical subgroups. Nevertheless, this classification is relatively coarse and ignores substantial clinical heterogeneity among different subtypes within hemorrhagic stroke. We therefore conducted a further refined subgroup analysis focusing on subarachnoid hemorrhage (I60) and intracerebral hemorrhage (I61), as presented in Table A3 in Multimedia Appendix 1. The results showed that the model performed stably in the I61 subgroup but yielded unsatisfactory outcomes in the I60 subgroup. This discrepancy can be explained by the sample characteristics of the I60 population. This subgroup had a small sample size and a high degree of internal homogeneity. Among the total of 247 patients with I60, 170 were categorized as high-cost cases, indicating highly consistent hospitalization costs in this subtype. Limited information from routine hospital discharge data failed to capture subtle individual differences, which impaired the model’s discriminative ability. In future research, we will expand the sample size and incorporate more clinical indicators, treatment-related information, and other variables to develop refined models for individual stroke subtypes.

In the global feature importance ranking, LOS was the second most important feature, after shortest distance. Previous studies have likewise highlighted the critical role of LOS in identifying high-cost pediatric inpatients [42]. SHAP interaction analysis further showed that the association between LOS and comorbidity network features was strongly time dependent, suggesting that the trajectory of medical resource use in stroke may evolve dynamically over the course of the disease. In the subgroup with an LOS of 8 to 10 days, higher NHCP showed a clear positive synergistic interaction, which may reflect the clinical profile of complex cases requiring highly resource-intensive interventions during the acute phase. However, when LOS was ≥14 days, high NHCP instead showed a strong negative interaction; at the same time, this long-stay subgroup exhibited a positive interaction peak at moderate levels of shortest distance. This time-varying pattern may suggest a law of diminishing marginal costs: as hospitalization becomes substantially prolonged, the main drivers of high expenditure may shift from fluctuations in acute disease severity to the cumulative baseline costs of ongoing care. Meanwhile, the 11 to 13 days subgroup showed a U-shaped pattern, with negative interactions across most of the value range, which may indicate a transitional stage in cost generation as patients move from the acute phase to the longer-term recovery phase.

Age and insurance type were also important determinants of the model’s identification performance. In particular, the effect of age on high hospitalization costs showed a nonlinear pattern, with older age exerting a negative effect on the identification of high-cost patients. One possible explanation is that because comorbid degenerative conditions increase surgical risk [45], older patients are more likely to receive conservative treatment and, compared with younger patients, are less likely to undergo high-cost interventional procedures, such as mechanical thrombectomy [46]. This pattern of clinical management is consistent with the baseline finding that patients in the high-cost group were, on average, younger. In addition, this study found a relatively high proportion of patients with a low comorbidity burden in the high-cost group. As patients in this group were generally younger, the prevalence of chronic underlying diseases would be expected to be lower, which may partly explain the clustering of relatively mild comorbidity in this group. However, given the inherent limitations of retrospective observational studies, future research should further examine the causal relationships among age, comorbidity, and high hospitalization costs in patients with stroke through prospective study designs. Notably, the performance of CCI, comorbidity count, and ECI in between-group comparisons was not entirely consistent, which may be attributable to differences in disease composition and weighting structure among these comorbidity measures. Comorbidity count reflects only the number of coexisting conditions and does not incorporate disease-specific weights. In contrast, the CCI applies weights based on a relatively limited set of comorbidity categories, whereas the ECI encompasses a broader range of comorbidity conditions and is based on a different scoring framework. Regarding insurance type, patients covered by provincial medical insurance generally benefit from higher reimbursement rates and broader formulary coverage. As a result, they may have greater financial access to high-cost diagnostic and therapeutic services and, when treated at tertiary grade A hospitals, may be more likely to accept more comprehensive or resource-intensive interventions. This may explain the significant positive association between provincial medical insurance and high-cost status. In contrast, patients with cross-provincial medical insurance or other payment methods may be constrained by reimbursement limits or out-of-pocket affordability. In clinical decision-making, these patients may therefore be more likely to adopt more economical or conservative treatment strategies, which is reflected in the model as a lower probability of high medical expenditure. Clinical management in these groups may prioritize cost-effective care, resulting in a negative association with high-cost status. This observed pattern is consistent with the findings of Yang et al in the Chinese ischemic stroke population [47].

Sensitivity analysis showed that when the outcome threshold was set at 5%, the constructed model performed poorly, with limited practical value for management purposes. As the threshold was relaxed from 10% to 20%, the model’s AUC tended to stabilize, while the F₁-score increased. However, this improvement in F₁-score may mainly reflect the passive effect of a larger number of positive samples and a more balanced class distribution, rather than a true improvement in the model’s discriminative ability. In addition, a broader threshold of 15% or 20% may dilute the core characteristics of the truly high-expenditure group, thereby creating potential challenges for the precise allocation of clinical resources and the refined management of medical insurance cost control. In contrast, a 10% threshold allows greater concentration on the core high-cost group and more accurately targets the disproportionately high-expenditure cases prioritized under the current medical insurance payment system. It should be noted, however, that at this 10% prevalence level, the optimal model yields a positive predictive value of approximately 0.41. While this implies a certain false-positive fraction, these misclassified cases typically represent moderately complex patients who may still benefit from subsequent care management. Ultimately, the practical implementation of this model relies on available administrative capacity. If hospital resources are relatively sufficient, health care managers can effectively absorb the interventions for these false positives and may even appropriately increase the threshold to expand coverage, thereby minimizing the risk of overlooking other potential high-cost cases.

First, the training and test sets in this study were generated using a single stratified random split, and the model evaluation results may therefore be influenced by variability arising from random sample partitioning. Although bootstrap resampling was used to estimate CIs, it could not fully eliminate the bias inherent to a single split. Due to the constraints of the study design for temporal validation and to preserve the original temporal characteristics of the cohort, repeated data splitting was not performed. Future studies could further enhance the reliability of the results through repeated split-sample validation. Second, the data were primarily sourced from a single center, which may introduce biases related to region-specific medical practices or insurance policies. Future research should prioritize multicenter external validation to confirm the model’s robustness; furthermore, the generalizability of this network feature-based framework to other disease spectra remains to be validated. Third, given the retrospective nature of this study and its reliance on hospital discharge data, some features exhibit information lag, which limits the utility of the model for preadmission screening. Although we used features available at admission to construct an early identification model, achieving more precise screening will require incorporating other features available at admission based on this study. However, SHAP-based interpretability analysis still holds significant value in providing decision support for clinical process management and guiding personalized interventions. Additionally, a specific limitation exists regarding the “planned readmission within 31 days” indicator. Owing to our 30-day aggregation strategy, merged readmission episodes accrue higher total costs, introducing a potential circularity issue that could theoretically inflate this feature’s apparent contribution. Although our SHAP analysis revealed it had the lowest importance contribution—indicating a minimal practical impact on the model’s identification performance in this cohort—future studies should carefully account for this potential confounding effect when constructing cost-related variables. Fourth, although the combination of comorbidity network features enriches the feature space, the current absence of biomarkers may limit the granularity of disease severity representation. Future research should prioritize the development of multimodal fusion models to evaluate the incremental value of biomarkers relative to network topological features, thereby establishing a more comprehensive and multidimensional risk profile for patient expenditures. Fifth, the static network constructed in this study may not fully capture high-order dependencies or the dynamic trajectories of disease evolution. Therefore, future research should integrate advanced algorithms, such as Graph Neural Networks and temporal network analysis [48,49], to unlock new opportunities for mining deep comorbidity associations and optimizing identification performance. Sixth, certain limitations remain in our data processing. Continuous variables were neither standardized nor normalized prior to KNN imputation and model training. Although our optimal model, XGBoost, is invariant to feature scaling, the lack of normalization may have impaired the performance of distance-based algorithms, such as SVM and NN. Additionally, despite a low missing data rate (0.13%‐1.95%), applying KNN imputation without prior scaling is a methodological limitation that may introduce minor computational biases. Future studies using distance-sensitive algorithms should incorporate a comprehensive data scaling pipeline. Finally, total hospitalization costs were adjusted for economic fluctuations using the Consumer Price Index. As the general Consumer Price Index primarily tracks a standard basket of consumer goods and medical inflation often outpaces general inflation, this approach may slightly misestimate the true medical cost inflation.

By integrating comorbidity network analysis with machine learning, this research used hospital discharge data to develop a high-performance framework for the identification of high-cost stroke patients. Our results confirm that, compared to stroke type within baseline features and conventional comorbidity features, network topological features are able to capture deeper information regarding disease complexity. The inclusion of network features improved model performance, with the XGBoost model demonstrating the best identification performance. Furthermore, SHAP analysis identified strong associations between high costs and relevant features while also revealing complex nonlinear interactions among these variables. These findings indicate that this framework demonstrates promising application prospects for cost-risk stratification in patients with stroke. If validated externally and prospectively in the future, this model could provide valuable decision support for exploring early identification of medical costs and stratified intervention strategies, thereby assisting in the optimization of medical resource allocation.