This study utilized data from the CHARLS, which implements multistage stratified sampling with probability-proportional-to-size weighting based on demographic stratification. The survey encompasses 150 county-level units across 28 provincial administrative regions in China. The baseline survey was conducted in 2011, with follow-up waves completed in 2013, 2015, 2018, and 2020, collecting comprehensive data on demographic characteristics, socioeconomic status, health behaviors, and medical history. The study protocol obtained ethical certification from Peking University’s Biomedical Ethics Committee (Approval ID: IRB00001052-11015). Sample selection followed three inclusion criteria: (1) age ≥60 years at baseline; (2) exclusion of individuals with pre-existing depression diagnosis or without basic or instrumental activities of daily living (BADL/IADL) disability at baseline; (3) completion of at least two consecutive follow-up assessments. Through integration of baseline (2011–2013, N = 2,440) and follow-up (2013–2015, N = 2,943) data, we constructed a longitudinal panel dataset containing 5,383 observations (2011–2015). The dataset was partitioned using stratified random sampling, allocating samples in a 7:3 ratio to training (N = 3,768) and testing (N = 1,615) sets. An independent 2018–2020 follow-up cohort (N = 3,254) served as the external validation set. The study flowchart is presented in Figure 1.

Depression was assessed using the 10-item Center for Epidemiological Studies Depression Scale (CESD-10), a widely validated screening tool for depressive symptoms in older adults (Zhang et al., 2024; Wang et al., 2023). The scale demonstrates good reliability and validity (Chen and Mui, 2014; Cheng and Chan, 2005). In CHARLS, the CESD-10 evaluates the frequency of 10 symptoms experienced in the past week: feeling bothered, trouble concentrating, feeling depressed, difficulty in doing things, feeling hopeful about the future (reverse-coded), feeling fearful, restless sleep, feeling happy (reverse-coded), feeling lonely, feeling unable to carry on. Each item is scored from 0 to 3, yielding a total score ranging from 0 to 30. Following established international criteria, a score ≥10 was used to define clinically significant depressive symptoms.

Disability in this study specifically refers to functional limitations in BADL or IADL, consistent with geriatric assessment standards (Katz et al., 1963; Lawton and Brody, 1969). BADL evaluates six fundamental self-care functions (dressing, bathing, eating, bed transferring, toileting, and continence control). IADL assesses five complex daily living skills (housekeeping, cooking, shopping, medication management, and financial management), where each item was scored 1 for inability to perform independently or 0 otherwise, resulting in total score ranges of 0–6 for BADL and 0–5 for IADL. The BADL/IADL-based criteria were applied during data screening to select eligible participants with BADL≥2 or IADL≥2 scores. This operational definition excludes sensory or cognitive disabilities alone, ensuring a homogeneous cohort with physical functional impairments.

The definitions and measurement approaches of covariates encompassed four domains: (1) demographic characteristics (gender, age, registered residence, educational level, marital status, number of children, and region); (2) health behaviors, including chronic disease history (14 conditions such as hypertension and diabetes), sensory functions (visual, auditory, and oral health assessed through assistive device use, functional scores, and tooth loss status), bodily pain, sleep time, physical activity intensity, social engagement, and lifestyle factors; (3) subjective perceptions comprising episodic memory, cognitive ability, life satisfaction, and self-rated health; (4) health care and insurance, incorporating health insurance type, healthcare utilization (frequency, duration, and costs of inpatient and outpatient services), as well as pension status, with detailed variable specifications and coding schemes provided in Supplementary Table 1.

Data preprocessing included four key steps: (1) Outlier handling: We applied the interquartile range (IQR) method to detect and truncate outliers for all continuous variables. Values beyond ±1.5IQR of the 25th-75th percentile range were clipped to the lower/upper bounds. This mitigated the impact of extreme values on tree-based models while preserving data distribution integrity. (2) One-hot encoding: Categorical variables (e.g., gender, region) were converted into binary dummy variables to avoid misinterpreting ordinal relationships. (3) Normalization: Continuous features were standardized using z-score normalization (mean = 0, variance = 1) to enhance convergence speed for linear models. (4) Missing value imputation: Among 78 candidate variables, 4 variables (5.13%) with missing rates >30% were excluded; (5) Missing value imputation: The remaining 74 variables had an average missing rate of 8.27% (range: 0.09–25.18%). A total of 15,918 missing records (8.82% of total training observations) were iteratively imputed using the MissForest algorithm. (6) Class imbalance adjustment: We implemented the SMOTE-Tomek hybrid sampling technique, combining synthetic minority oversampling (SMOTE) with Tomek links under sampling. This approach effectively enhanced the model’s sensitivity in detecting depression risk and improved clinical utility by generating synthetic samples.

We refer to existing studies for feature selection (Huang et al., 2025; Zheng et al., 2024; Zhu et al., 2025). A three-stage serial consensus approach was utilized to identify robust predictors in this study. This serial consensus approach integrates complementary strengths of distinct selection paradigms. (1) LASSO (L1 regularization): Efficiently screens out zero-importance features (73 variables) by imposing sparsity constraints. Although its linearity assumption may oversimplify relationships, it serves as a high-recall initial filter. (2) Elastic Net (L1 + L2 regularization): Reduces multicollinearity-induced instability by retaining correlated but biologically plausible features. The α = 0.5 setting balances sparsity and grouping effects, mitigating LASSO’s limitation in correlated feature selection, retaining 42 stable features. (3) The Elastic Net output variables were fed into the Boruta algorithm, which identified 28 significant predictors by comparing random forest importance scores with shadow variables (p < 0.01). To ensure reproducibility, a random seed (random_state = 42) was set for both the MissForest imputation and Boruta’s shadow variable generation.

Integration of feature selection results via a strict intersection strategy. The three feature selection outcomes were consolidated through a stringent intersection strategy. Specifically, we quantified the selection frequency of each variable across LASSO, elastic net, and Boruta algorithms, retaining only variables unanimously selected by all three methods (i.e., frequency ≥3). This approach yielded 21 high-confidence predictors, including age, self-rated health, arthritis, renal disease, stomach, asthma, memory-related disorders, observe the situation up close, hearing ability, self-reported pain in head, wrist, leg, toes, neck, sleep time, social activities, episodic memory, life satisfaction, medical insurance types, hospitalization expenses (total expenses), outpatient expenses (out of pocket expenses) (as shown in Figure 2). Compared to individual methods, this strategy significantly enhanced feature stability.

Ten ML algorithms were implemented, including logistic regression (LR), support vector machine (SVM), extreme gradient boosting (XGBoost), light gradient boosting machine (LightGBM), categorical boosting (CatBoost), random forest (RF), bootstrap aggregating (Bagging), histogram-based gradient boosting machine (HistGBM), multilayer perceptron (MLP), and decision tree (DT). To optimize model generalizability, hyperparameter tuning was performed using grid search with 3-fold stratified cross-validation (specific hyperparameter configurations are provided in Supplementary Table 2). Model performance was comprehensively evaluated through five metrics: (1) the area under the receiver operating characteristic curve (AUC), measuring the model’s ability to discriminate between positive and negative cases (Formulas 1–3); (2) accuracy, representing the proportion of correctly classified samples (Formula 4); (3) precision, indicating the ratio of true positives among all predicted positives (Formula 5); (4) recall, reflecting the model’s capacity to identify actual positive cases (Formula 6); and (5) the F1-score, the harmonic mean of precision and recall, which provides a balanced assessment of the model’s performance on the positive class (Formula 7). The mathematical formulations were derived from established methodologies (Theerthagiri, 2022; Yang and Ying, 2022).

In these formulations, TPR denotes the true positive rate, FPR represents the false positive rate, TP indicates true positives, FP signifies false positives, TN refers to true negatives, and FN stands for false negatives.

Statistical analyses were performed using Stata 18.0 for data description and Python 3.13 for subsequent modeling. Continuous variables were characterized differentially based on their distribution: normally distributed variables were presented as mean ± standard deviation, while non-normally distributed variables were summarized using median and interquartile range, with normality assessed via the Shapiro–Wilk test. Categorical data were expressed as cardinality measures (absolute frequencies) with proportional composition. The statistical significance threshold was set at p < 0.05 for all analyses.

The analysis included samples from the training sets (N = 3,768), testing sets (N = 1,615), and external validation sets (N = 3,254). Table 2 summarizes the baseline demographic characteristics, disease status, and depression status across the three cohorts. The median ages were 71, 71, and 72 years in the training, testing, and validation sets, respectively, with statistically significant inter-group differences (all p < 0.05). In terms of gender, the proportion of females is similar in the training sets (57.94%), testing sets (57.28%), and external validation sets (60.97%). There was no statistically significant difference in gender distribution between groups (all p > 0.05), indicating that the gender ratio remained balanced in the data partitioning. In terms of marital status, the married group accounts for 70.25, 70.59, and 69.61% of the three groups, respectively, which is much higher than the unmarried group. In terms of registered residence, the proportion of rural registered residence registration slightly decreased in training sets (79.91%), testing sets (79.57%), and external verification sets (78.83%), but all exceeded 78%. There was a significant difference in the distribution of registered residence among groups (p < 0.05). The proportion of “3 or more children” in the three groups was 72.24, 73.32, and 66.04%, respectively, with a significant decrease in the proportion of external validation sets. There were significant differences in distribution between groups (all p < 0.05). There was no significant difference in regional distribution between groups (all p > 0.05). In terms of education level, the proportion of people who have not received formal education gradually decreased in the training sets (68.21%), testing sets (69.05%), and external validation sets (65.43%), while the proportion of high school and above education increased from 4.03 to 5.53%. There was a significant difference between the groups (all p < 0.05). In terms of disease characteristics, there was no significant difference (p > 0.05) in the prevalence of memory-related diseases and stroke diseases among the training sets, testing sets, and external validation sets. The incidence of heart disease was significant in the training and testing sets (p < 0.05), but not significant in the validation sets (p > 0.05). The incidence of arthritis disease remained stable among the three groups (56.22–56.61%), with no significant difference between the groups (p > 0.05). The proportion of depression showed a significant increasing trend among the training sets (56.32%), testing sets (56.35%), and validation sets (64.20%), with no significant difference (p > 0.05).

Figure 3 reveals the demographic differences in the prevalence of depression among disabled individuals. The gender distribution shows that the prevalence of depression in the female population (80.80%) is significantly higher than that in the male population (37.15%). Analysis of marital status shows that unmarried individuals have a higher risk of depression (61.18%) compared to married individuals (53.35%). In age stratification, the prevalence of depression in the elderly group aged 80 and above reached 60.06%, which was higher than that in the 70–80 age group (52.49%) and the 60–70 age group (56.41%). The regional distribution shows that the incidence rate in the western region (59.24%) and rural areas (57.29%) is significantly higher than that in the eastern region (52.50%) and urban areas (49.07%). Education level analysis shows that the illiterate population has the highest incidence of disease (57.86%), and there is a non-linear relationship between educational attainment and morbidity probability. The dimension of family support shows that the risk of depression in the childless group (75.18%) is significantly higher than that in the childbearing group (53.07–56.50%).

This study systematically evaluated the performance of 10 ML algorithms in predicting depression risk among older adults with disability across training, testing, and external validation sets (as shown in Table 3). In terms of accuracy, RF (0.741), LightGBM (0.728), and HistGBM (0.713) demonstrated the highest performance in the testing sets. While LR and DT exhibited relatively stable performance between training and testing sets, their overall accuracy was the lowest (LR: 0.667; DT: 0.633). For AUC metrics, RF (0.797) achieved the strongest discriminative capacity, followed closely by LightGBM (0.785), XGBoost (0.781), and HistGBM (0.779). XGBoost, LightGBM, and HistGBM showed superior generalizability, whereas DT (0.636) performed the poorest. Regarding the F1-score, RF (0.762) exhibited the optimal balance between precision and recall, with LightGBM (0.749), CatBoost (0.741), and HistGBM (0.735) maintaining stable performance in the testing sets. For precision, RF (0.791), LightGBM (0.778), XGBoost (0.767), and HistGBM (0.723) achieved the highest positive predictive values and lowest false positive rates, significantly outperforming DT (0.691). In recall analysis, RF (0.735), HistGBM (0.723), and HistGBM (0.707) demonstrated the strongest ability to identify true positive cases, while DT (0.633) exhibited markedly higher missed-detection risks compared to ensemble methods.

Through comprehensive evaluation of 10 ML models, HistGBM was selected as the optimal model based on three key criteria: (1) superior performance on testing sets metrics (AUC = 0.779, F1-score = 0.735, accuracy = 0.713), (2) excellent generalizability demonstrated by a minimal training–testing AUC gap (8.5%), and (3) consistent performance across validation sets. HistGBM exhibited well-balanced predictive capabilities, showing above-average performance across all evaluation metrics without significant weaknesses in either precision (0.766) or recall (0.707), indicating robust discriminative power between positive and negative cases along with stable predictive performance. The model’s exceptional generalizability was particularly noteworthy, with only a 10% difference in AUC between the testing sets and validation sets, as shown in Figure 4, significantly outperforming other models and demonstrating strong robustness across different data distributions. Although RF achieved the highest individual metrics on the specific testing set used in Table 3, the substantial performance degradation observed on the external validation set raised concerns about its real-world applicability and stability. HistGBM, while having marginally lower peak testing set scores than RF, offered the best overall package of strong predictive performance, minimal overfitting (small train-test gap), and exceptional stability across the independent validation set. This superior generalizability was the primary reason for selecting HistGBM as the optimal model for potential clinical application, where reliability across diverse data sources is paramount.

XGBoost demonstrated strong performance on testing sets metrics (AUC = 0.781, F1-score = 0.735, accuracy = 0.713), though it exhibited a relatively large training–testing AUC gap (10.1%) compared to LR, HistGBM, MLP, and CatBoost. However, its validation set’s AUC showed a 10.7% difference from testing sets’ performance (Figures 5–7), suggesting reasonable stability across data partitions and potential suitability for resource-constrained scenarios. RF achieved excellent testing set results (AUC = 0.797, F1-score = 0.762, accuracy = 0.741), but displayed concerning generalization issues with substantial training–testing and testing-validation AUC differences, indicating potential overfitting to training data noise or specific patterns. LightGBM ranked second in testing sets AUC (0.785) with stable validation performance (0.654). However, the difference between the testing set AUC and the validation set is too large (13.1%). CatBoost performed comparably to top models in testing set metrics (AUC = 0.774, F1-score = 0.741) with excellent categorical feature handling, though its relatively lower F1-score (0.741) and recall (0.720) scores suggested weaker minority class identification, limiting its utility for imbalanced datasets. Among the remaining models, SVM exhibited severe overfitting, while MLP and DT significantly underperformed ensemble methods in both AUC and F1-score metrics.

The SHAP values quantify the absolute average impact of each feature on model predictions across all possible feature combinations, revealing their global importance. As shown in Figure 8, the SHAP analysis of the HistGBM model demonstrated significant variability in feature contributions. Sleep time (mean SHAP = 0.344), life satisfaction (0.339), episodic memory (0.220), and self-rated health (0.197) emerged as the top four predictive features, indicating that health behaviors, subjective perceptions, and cognitive function were the core drivers of model predictions. The high contribution of sleep time likely reflects its well-established associations with chronic diseases, metabolic disorders, and cognitive decline. Life satisfaction and self-rated health, as subjective health indicators, capture the interplay between psychosocial factors and physiological health. Episodic memory directly influences prediction through cognitive and sensory pathways. Moderate contributions were observed for features such as stomach diseases, observing the situation up close, and memory-related disorders. While self-reported pain in the head, wrist, leg, toes, neck, and mental health conditions showed limited predictive importance, suggesting either weak signals or sparse data distributions. These results validate the model’s multidimensional feature selection approach and provide actionable insights for intervention prioritization. Health management strategies targeting high-contribution features could enhance the model’s real-world utility. Additionally, domain knowledge should guide the evaluation of low-contribution features to optimize the balance between model complexity and interpretability.

Figure 9 presents the SHAP value distributions, revealing the heterogeneous directional effects and magnitudes of various features on depression probability predictions among older adults with disability. The x-axis (SHAP value) indicates each feature’s influence on model output, where positive values increase and negative values decrease predicted risk. The color gradient (red means high feature value, blue means low feature value) demonstrates that: (1) higher values of sleep time, life satisfaction, and self-rated health (red clusters with negative SHAP) were strongly protective against depression, consistent with established epidemiological mechanisms; (2) better episodic memory performance (blue with positive SHAP) correlated with reduced depression risk, potentially through preserved cognitive resilience; (3) stomach diseases (red with positive SHAP) elevated risk through chronic somatic burden and psychological stress pathways; and (4) bodily pain (head, wrist, leg, toes, neck; red with positive SHAP) increased depression vulnerability in this population. These findings highlight the central role of health behaviors and psychosocial factors in depression comorbidity risk while identifying specific physiological pain features as contributory predictors.