Study design

This secondary dataset diagnostic-classification modeling study with internal validation utilized cross-sectional data from the NHANES spanning three consecutive survey cycles from 2007 to 2012. This specific timeframe was selected because only these three cycles contained complete pulmonary function test data, including measurements of forced expiratory volume in 1 second (FEV₁) and forced vital capacity (FVC). NHANES employs a complex multistage probability sampling design to ensure national representativeness of the non-institutionalized civilian population of the United States. The NHANES study protocol received approval from the Research Ethics Review Board of the National Center for Health Statistics (NCHS), and all participants provided written informed consent.

Inclusion criteria required participants to have completed pulmonary function testing that met ATS/ ERS standards, with test quality grades of A or B. Exclusion criteria encompassed missing pulmonary function data or quality grades of C or D, aged <40 years, and pregnant women. The latter group was excluded because physiological pulmonary function changes during pregnancy may compromise the accuracy of FEV₁/FVC measurements.

From the initial sample of 30442 individuals across the three NHANES cycles from 2007 to 2012, we first excluded 11366 participants with missing pulmonary function data, followed by the exclusion of 11015 participants who were either younger than 40 years or pregnant women. This yielded a final analytical sample of 8061 eligible participants.

Analytical strategy

This study employed three methodological approaches for COPD classification. First, LASSO regression with repeated 10-fold cross-validation was used to systematically select optimal features from 33 candidate variables spanning demographics, anthropometrics, lifestyle factors, comorbidities, and laboratory biomarkers. LASSO’s L1 regularization prevents overfitting while enhancing model parsimony. Second, 12 machine learning algorithms were compared, including traditional methods (logistic regression, linear discriminant analysis) and ensemble approaches (Random Forest, XGBoost, CatBoost, LightGBM), to identify the optimal classification framework. Third, SHAP analysis was applied to quantify individual feature contributions, revealing non-linear relationships that traditional methods cannot detect. Clinical utility was evaluated through Decision Curve Analysis.

Outcome variable definition

This study employed Global Initiative for Chronic Obstructive Lung Disease (GOLD) guideline-recommended spirometric criteria to define COPD^5,6, whereby a post-bronchodilator FEV₁/FVC ratio of <0.70 was classified as COPD and a ratio of ≥0.70 as non-COPD. Spirometry was performed using calibrated Ohio 822/827 dry-rolling seal spirometers following ATS standardized protocols. Eligible participants received albuterol and underwent post-bronchodilator spirometry after a minimum 10-minute waiting period. Participants performed repeated maximal forced expiratory maneuvers until achieving reproducible results, with a maximum of eight attempts permitted. The highest FEV₁ and FVC values from acceptable maneuvers were selected to calculate the FEV₁/FVC ratio. To ensure data quality, this study exclusively included spirometry measurements with quality grades of A or B, where Grade A exceeded ATS standards and Grade B met ATS standards for reproducibility.

Feature variable selection

Based on existing literature and clinical expertise, this study systematically extracted 36 candidate feature variables associated with COPD from the NHANES database^7,8. These variables encompassed multiple dimensions, including demographic characteristics (age, sex, race/ethnicity, education level, and poverty-income ratio), anthropometric measurements (body mass index and waist circumference), lifestyle factors (smoking status, smoking years, and physical activity level), comorbidities (asthma, cardiovascular disease, hypertension, diabetes, and chronic kidney disease [CKD]), and laboratory indicators.

Smoking status was categorized as never smoker (<100 cigarettes in lifetime), current smoker (≥100 cigarettes in lifetime and currently smoking), or former smoker (≥100 cigarettes in lifetime but not currently smoking). Laboratory indicators were further subdivided into hematological parameters (hemoglobin, hematocrit, white blood cell count, neutrophil percentage, lymphocyte percentage, and eosinophil percentage), inflammatory markers (neutrophil-to-lymphocyte ratio [NLR], platelet-to-lymphocyte ratio [PLR], and systemic immune-inflammation index [SII]), nutritional indicators (serum albumin, vitamin D, advanced lung cancer inflammation index [ALI], geriatric nutritional risk index [GNRI], and prognostic nutritional index [PNI]), lipid profile (total cholesterol, high-density lipoprotein cholesterol [HDL-C], and triglycerides), hepatorenal function indicators (alanine aminotransferase [ALT] and estimated glomerular filtration rate [eGFR]), and glucose metabolism indicators (fasting glucose and glycated hemoglobin). All laboratory examinations were performed at the NHANES Mobile Examination Center by certified technicians following standardized operating procedures, with rigorous quality control measures implemented throughout.

Data preprocessing

This study employed a systematic approach to handle missing data. After assessing missingness, three variables (triglycerides, fasting glucose, glycated hemoglobin) exceeded 10% missing rate and were excluded.

The dataset was randomly partitioned into training and validation sets at a 7:3 ratio. The training set comprised 70% of the total sample and was utilized for model training and hyperparameter optimization. The validation set comprised the remaining 30% of the sample and served as an independent validation set for final model performance evaluation. The validation set data remained completely unseen during the model training and hyperparameter optimization processes and was employed solely for independent performance verification after all model training was completed, thereby ensuring the objectivity and unbiasedness of the evaluation results.

For the remaining variables with minimal missing values, different imputation methods were applied according to variable type, with all imputation and standardization parameters derived exclusively from the training set and subsequently applied to the validation set. For continuous variables, distribution characteristics were first assessed using the Shapiro-Wilk test. Variables following a normal distribution were imputed using mean imputation, while those exhibiting skewed distributions were imputed using median imputation to enhance robustness. For categorical variables, the K-Nearest Neighbors (KNN) imputation method was employed with K set to 5, selecting the most similar neighboring samples based on Euclidean distance to infer missing values. These conservative single imputation methods were chosen because: 1) missing rates were minimal across all retained variables; and 2) neither mean/median nor KNN imputation incorporates outcome information, thereby minimizing the risk of information leakage. Additionally, to eliminate the influence of different measurement scales and numerical ranges across variables, all continuous variables underwent Z-score standardization, while categorical variables retained their original coding format.

LASSO-based variable selection

To avoid overfitting and multicollinearity issues arising from high-dimensional data, this study employed the Least Absolute Shrinkage and Selection Operator (LASSO) logistic regression to perform feature selection on the 33 candidate variables. LASSO introduces an L1 regularization penalty term into the loss function, which automatically shrinks the regression coefficients of unimportant variables to zero, thereby achieving feature dimensionality reduction and model simplification. This study implemented LASSO regression using the glmnet package in R software, with the optimal regularization parameter λ determined through repeated 10-fold cross-validation (five repetitions). According to the one-standard-error rule, lambda.1se was selected as the final regularization parameter, which corresponds to the largest λ value within one standard error of lambda.min (the λ value that minimizes cross-validation error). This strategy yielded a more parsimonious and robust model while sacrificing minimal classification accuracy. To visualize the variable selection process, LASSO coefficient path plots and cross-validation error curves were generated. Based on the lambda.1se criterion, 16 variables with non-zero regression coefficients were successfully selected from the 33 candidate variables as the final feature set for subsequent machine learning model construction.

Machine learning model construction

Based on the 16 features variables selected by LASSO from the training set, this study constructed 12 different types of machine learning models to comprehensively compare their performance in COPD classification tasks. These models included logistic regression (LR) as the traditional baseline model, linear discriminant analysis (LDA) for linear classification, support vector machine with radial basis function kernel (SVM-RBF) to handle nonlinear relationships, K-Nearest Neighbors (KNN) as an instance-based nonparametric method, and neural network (NN) to learn complex feature representations through multilayer perceptrons. Additionally, this study emphasized the evaluation of multiple ensemble learning algorithms, including Random Forest (RF), Gradient Boosting Machine (GBM), Extreme Gradient Boosting (XGBoost), Light Gradient Boosting Machine (LightGBM), Categorical Boosting (CatBoost), Adaptive Boosting (AdaBoost), and Partial Least Squares Discriminant Analysis (PLS-DA).

All hyperparameter tuning was performed exclusively on the training set using grid search with 5-fold cross-validation; the validation set remained completely unseen during model development. For models implemented within the caret framework, hyperparameters were optimized using classification accuracy as the selection criterion, with algorithm-specific tuning grids defined for each model (e.g. mtry for Random Forest, interaction depth and number of trees for GBM, K for KNN). CatBoost and LightGBM were tuned using their native implementations with custom grid search over key hyperparameters (tree depth, learning rate, and regularization coefficients), using 5-fold cross-validation with AUC as the optimization metric and early stopping (patience of 50 rounds, maximum 500 iterations) to prevent overfitting. Final models were retrained on the full training set with optimal hyperparameters and then evaluated on the held-out validation set.

NHANES sampling weights were incorporated in all descriptive analyses (Tables 1 and 2) to ensure nationally representative population estimates, but were not applied during machine learning model training. The primary modeling objective was to optimize individual-level classification accuracy rather than to estimate population-level parameters, and incorporating complex survey weights into iterative machine learning algorithms lacks established theoretical frameworks and standardized implementations.

Model evaluation and validation

Performance evaluation of all models was conducted on the independent validation set. This study employed the area under the receiver operating characteristic curve (AUC) as the core metric for model selection, as AUC comprehensively reflects a model’s ability to discriminate between COPD patients and non-COPD individuals across all possible classification thresholds. To comprehensively evaluate and visualize model performance, this study generated multiple graphical representations.

Receiver operating characteristic (ROC) curves were plotted for all 12 models with their respective AUC values annotated, and all curves were consolidated into a single figure to facilitate direct comparison. Forest plots summarized the AUC values and their 95% confidence intervals for all 12 models, enabling straightforward horizontal comparison of model performance differences. Uncertainty around classification metrics was quantified using 1000-iteration bootstrap resampling of the validation and training sets, with the optimal threshold re-derived via Youden’s index within each iteration; 95% confidence intervals were calculated using the percentile method. Residual box plots were generated to assess the distribution characteristics and stability of classification errors across models. Residuals were defined as the actual value (0 or 1) minus the predicted probability, with root mean square error (RMSE) calculated for each model. Calibration was assessed by dividing samples into 10 groups according to predicted probabilities, with scatter plots generated using mean predicted probability as the x-axis and observed COPD incidence rate as the y-axis. The diagonal line represents ideal calibration.

Decision curve analysis (DCA) evaluated the net benefit of models at different risk thresholds, with ‘treat all’ and ‘treat none’ strategies serving as reference baselines. Clinical impact curves were constructed to display the number of cases predicted as high-risk and the actual number of high-risk cases at different thresholds for the optimal model, providing an intuitive assessment of the model’s potential impact in clinical decision-making contexts.

Model interpretation

To elucidate the underlying mechanisms of the optimal classification model and identify key risk factors, this study performed SHapley Additive exPlanations (SHAP) analysis on the model with the highest AUC. The SHAP method, grounded in Shapley value theory from game theory, quantifies the marginal contribution of each feature to model classifications, thereby providing feature importance assessments that are both globally meaningful and locally interpretable. Multiple visualization plots were generated to interpret model behavior from different perspectives.

The SHAP summary plot (bee swarm plot) displays the distribution of SHAP values for all features, where each point represents a sample, color indicates the magnitude of the feature value, and the x-axis represents the SHAP value (positive values indicate increased COPD risk, while negative values indicate decreased risk). This plot intuitively presents feature importance and the directionality of their relationship with the outcome. The SHAP feature importance plot (bar chart) ranks features according to their mean absolute SHAP values, displayed as horizontal bars with specific numerical annotations, clearly identifying the most important risk factors.

SHAP dependence plots were generated for the top three important features, illustrating the dependence relationship between feature values and SHAP values, which facilitates the identification of nonlinear effects and threshold effects. The SHAP waterfall plot was constructed for representative cases to decompose individual predictions, demonstrating how each feature cumulatively contributes from the baseline value to the final predicted probability, thereby providing a transparent interpretation of prediction results for individual samples.

Statistical analysis

For descriptive statistics, continuous variables were expressed as median and interquartile range (IQR), while categorical variables were presented as frequencies and percentages. Group comparisons employed appropriate statistical tests, with weighted Mann-Whitney U tests applied to continuous variables and weighted chi-squared tests applied to categorical variables. To ensure that analytical results were representative of the United States population aged ≥40 years, statistical analyses incorporated NHANES-provided sampling weights, with weighted analyses under complex sampling design conducted using the survey package in R⁹.

All data processing, statistical analysis, and machine learning modeling in this study were performed using R software (version 4.4.1). Primary R packages utilized included tidyverse for data cleaning and processing, haven for reading XPT format data files, survey for weighted statistical analysis under complex sampling design, VIM for missing data visualization, glmnet for LASSO regression and variable selection, and caret for providing a unified framework and hyperparameter optimization for machine learning models^10,11. Additional packages including randomForest, e1071, nnet, gbm, adabag, and pls were employed to implement various machine learning algorithms. Furthermore, the xgboost, lightgbm, and catboost packages were used to implement extreme gradient boosting, light gradient boosting machine, and categorical boosting algorithms, respectively. For model evaluation and visualization, the pROC package was utilized for ROC curve analysis and AUC calculation, the rms package for generating calibration curves, the dcurves package for decision curve analysis, the ggplot2 package for high-quality graphical visualization, and the gridExtra and grid packages for figure layout and arrangement. For model interpretation, the kernelshap and shapviz packages were employed for SHAP value calculation and visualization analysis. All statistical tests were two-sided, with p<0.05 considered to indicate statistical significance.

Baseline characteristics

Table 1 presents a comparison of baseline characteristics across the overall population, training set, and validation set. Overall, the study participants had a median age of 54 years (IQR: 47–62), with males comprising 48.44% and females 51.56% of the sample. The training and validation sets demonstrated no significant differences in any baseline characteristics (all p>0.05), indicating excellent balance in dataset partitioning and ensuring the reliability of subsequent model evaluation.

Table 2 presents a comparison of baseline characteristics between the COPD and non-COPD groups. Regarding demographic characteristics, COPD patients were significantly older than non-COPD individuals (p<0.001). The proportion of males was significantly higher in the COPD group compared to the non-COPD group (p<0.001). Body mass index (BMI) was significantly lower in the COPD group than in the non-COPD group (p<0.001). Smoking status differed significantly between the two groups (p<0.001). In the COPD group, current smokers comprised 34.01% and former smokers comprised 38.28%, whereas never smokers comprised 57.60% in the non-COPD group. The number of smoking years was significantly higher in the COPD group compared to the non-COPD group, with median values of 22 years and 0 years, respectively (p<0.001).

LASSO variable selection results

This section presents the variable selection process and final results of LASSO regression. The optimal regularization parameter was determined through repeated 10-fold cross-validation (five repetitions), with lambda.1se selected as the final parameter according to the one-standard-error rule. Supplementary file Figure 1 illustrates the complete LASSO variable selection process. Supplementary file Figure 1a depicts the LASSO coefficient path plot, clearly demonstrating the dynamic trajectory of how regression coefficients for each variable gradually shrink toward zero as the penalty strength (log lambda) increases. Supplementary file Figure 1b shows the cross-validation error curve, displaying the binomial deviance corresponding to different lambda values along with their 95% confidence intervals, with two vertical dashed lines marking the positions of lambda. min and lambda.1se. Based on the lambda.1se criterion, LASSO regression successfully identified 16 variables with non-zero regression coefficients from the 33 candidate feature variables as the final feature set. Table 3 provides a detailed list of these 16 variables and their corresponding LASSO regression coefficients.

Performance comparison of machine learning models

This section systematically compared the classification performance of 12 machine learning models on both training and validation sets. Supplementary file Figure 2 presents the AUC value comparisons and ROC curves for all models. Supplementary file Figures 2a and 2b display forest plots of AUC values for the training and validation sets, respectively, with diamond-shaped points representing the AUC values of each model and horizontal error bars indicating 95% CI, facilitating intuitive comparison of performance differences across models. On the training set, LightGBM demonstrated optimal performance (AUC=0.840), followed by CatBoost (AUC=0.811) and Gradient Boosting Machine (AUC=0.788). On the validation set, CatBoost achieved the highest AUC value (AUC=0.788), outperforming LightGBM (AUC=0.753). Although LightGBM exhibited the best performance on the training set, its inferior performance on the validation set compared to CatBoost suggested potential overfitting. Supplementary file Figures 2c and 2d present the ROC curves for the training and validation sets, respectively. From the morphology of the ROC curves, it is evident that the curves for CatBoost and LightGBM are closest to the upper left corner, indicating the best classification performance.

Supplementary file Figures 3a and 3b present residual box plots for the training and validation sets, respectively. The validation set residual distribution revealed that CatBoost and Logistic Regression models exhibited the most compact residual boxes with smaller RMSE values, suggesting superior classification stability, consistent with their excellent AUC performance.

On the validation set, CatBoost achieved the highest AUC of 0.788 (95% CI: 0.766–0.810). DeLong’s test for pairwise AUC comparison demonstrated that CatBoost significantly outperformed seven models: SVM-RBF (0.727, p<0.001), KNN (0.711, p<0.001), PLS-DA (0.750, p<0.001), XGBoost (0.741, p<0.001), LightGBM (0.753, p<0.001), Random Forest (0.772, p=0.004), and AdaBoost (0.779, p=0.011). No statistically significant differences were observed between CatBoost and Logistic Regression (0.785, p=0.375), Neural Network (0.785, p=0.526), LDA (0.785, p=0.484), or GBM (0.782, p=0.070). Complete pairwise DeLong test results are provided in Supplementary file Table S1.

Table 4 presents performance metrics calculated at model-specific optimal thresholds for both training and validation sets. On the validation set, CATBoost achieved the highest Youden’s index of 0.434 with sensitivity of 0.706 and specificity of 0.728. Validation set sensitivity ranged from 0.607 to 0.763 and specificity from 0.609 to 0.786 across all models.

Clinical evaluation of the optimal model

Based on validation set AUC values, CatBoost was identified as the optimal classification model (AUC=0.788), and this section provides an in-depth clinical performance evaluation of this model. Supplementary file Figures 3c and 3d present the calibration curves for the CatBoost model on the training and validation sets, respectively. As shown in the Supplementary file Figures, the CatBoost model demonstrated excellent calibration performance on both datasets, with the vast majority of scatter points falling near the diagonal line, indicating high concordance between the model-predicted COPD risk probabilities and actual incidence rates.

Supplementary file Figures 4a and 4b present the decision curve analysis (DCA) results for all 12 models on the training and validation sets, respectively. From the validation set DCA curve (Supplementary file Figure 4b), it is evident that across a wide range of risk thresholds from 0.1 to 0.6, the decision curves for ensemble learning models such as CatBoost, LightGBM, and Random Forest are positioned above both baseline strategies, indicating that using these models to guide clinical decisions can achieve higher clinical net benefit compared to extreme strategies. The decision curve for the CatBoost model is the smoothest and positioned at the top, suggesting its strongest clinical utility. To facilitate clinical interpretation, we examined the CatBoost model’s performance at specific probability thresholds on the validation set (Supplementary file Table 2). At a 20% threshold, the model identified 823 individuals as high risk, capturing 314 of 456 COPD cases (sensitivity 68.9%) while generating 509 false positives, corresponding to 1.6 false positives per true positive and a positive predictive value of 38.2%. At a more conservative 30% threshold, 485 individuals were flagged, with 225 true positives (sensitivity 49.3%), 260 false positives, and 1.2 false positives per true positive. Compared with a ‘screen all’ strategy, the CatBoost model at the 20% threshold reduced unnecessary referrals for spirometry by approximately 37 per 100 individuals screened while maintaining a positive net benefit.

Supplementary file Figures 4c and 4d present the clinical impact curves for the CatBoost model on the training and validation sets, respectively. These Supplementary file Figures display the number of cases predicted as high-risk by the model (red solid line) and the actual number of high-risk cases (blue dashed line) at different risk thresholds. From the validation set clinical impact curve (Supplementary file Figure 4d), it can be observed that within the risk threshold range of 0.0 to 0.4, the red curve (predicted high-risk number) consistently remains slightly above the blue curve (actual high-risk number), indicating a mild tendency toward overestimation but with overall consistent trends. As the threshold increases above 0.6, both curves gradually approach zero, indicating that the number of cases identified as high-risk decreases sharply at high thresholds, which is consistent with the epidemiological characteristics of relatively low COPD prevalence in the general population.

Model interpretability analysis

To elucidate the classification mechanisms of the CatBoost model and identify key risk factors, this section employed the SHAP method for model interpretation. Supplementary file Figure 5a presents the SHAP feature importance plot, ranking features according to their mean absolute impact on model predictions The results revealed that age and smoking years were jointly the most important classification features, far exceeding other variables, indicating that these two factors exert the most significant influence on COPD risk. These were followed sequentially by BMI, race/ethnicity, sex, asthma, and smoking status.

Supplementary file Figure 5b displays the SHAP bee swarm plot, illustrating the distribution of SHAP values for all features and their relationship with COPD risk. Key patterns emerged from the analysis. Age and smoking years both demonstrated strong positive correlations with COPD risk, with higher values concentrated in positive SHAP regions. Conversely, high BMI (yellow points) appeared predominantly in negative SHAP regions, suggesting a protective effect of obesity against COPD. Female sex (yellow points) clustered in negative SHAP regions while male sex (purple points) was distributed in positive regions, validating the protective effect of female sex. Both asthma and current smoking status increased COPD risk.

Supplementary file Figure 5c presents SHAP dependence plots for the three most important features (age, smoking years, and BMI). For age, SHAP values increased approximately linearly from negative values at age 40 years to high positive values at age 80 years, with samples having longer smoking duration (yellow points) showing higher SHAP values at the same age, indicating a synergistic effect between age and smoking exposure. Smoking years similarly demonstrated a linear increasing trend, with elderly smokers exhibiting higher SHAP values than younger smokers. The BMI variable revealed a significant U-shaped nonlinear relationship with COPD risk. Low BMI represented a risk factor (positive SHAP values), moderate BMI (30–40 kg/m²) conferred protection (lowest SHAP values), while morbid obesity (>40 kg/m²) showed diminished protective effects. Color coding by asthma status revealed interaction effects: asthma significantly increased COPD risk at low BMI but appeared to reduce risk in the moderate BMI range.

Supplementary file Figure 5d presents the SHAP waterfall plot for a representative COPD case (Sample 6), demonstrating individual classification decomposition. Smoking years (29 years) contributed most substantially to COPD risk, followed by age (64 years). Female sex reduced the predicted risk (-0.177). After cumulative contribution of all features, the predicted COPD risk for this sample reached approximately 28% [f(x)= -0.917].

Limitations

This study has several important limitations. Most critically, the lack of external validation in independent cohorts, limits conclusions about generalizability. All performance metrics derive from internal validation using train-test splitting from the same NHANES cycles, which cannot assess performance when applied to different populations, geographical regions, or time periods. External validation in diverse settings is essential before clinical implementation. The cross-sectional design precludes establishing temporal relationships or inferring causality. The model classifies concurrent rather than incident COPD, limiting its utility for primary prevention. Outcome misclassification is possible despite Grade A/B spirometry quality standards, as the GOLD fixed-ratio criterion (FEV₁/FVC <0.70) may overdiagnose COPD in older adults. Additionally, machine learning models were trained without survey sampling weights, which may limit the representativeness of model parameters for the broader US population, although the impact on individual-level discriminative performance is expected to be modest. Additional limitations include unmeasured confounders (detailed occupational exposures, residential air pollution, biomass fuel exposure, childhood respiratory infections, and genetic variants), which likely contribute to the moderate AUC of 0.788.