Lung Radiomics Features Selection for COPD Stage Classification Based on Auto-Metric Graph Neural Network

Chronic obstructive pulmonary disease (COPD) is a preventable, treatable, progressive chronic disease characterized by persistent airflow limitation. Patients with COPD deserve special consideration regarding treatment in this fragile population for preclinical health management. Therefore, this paper proposes a novel lung radiomics combination vector generated by a generalized linear model (GLM) and Lasso algorithm for COPD stage classification based on an auto-metric graph neural network (AMGNN) with a meta-learning strategy. Firstly, the parenchyma images were segmented from chest high-resolution computed tomography (HRCT) images by ResU-Net. Second, lung radiomics features are extracted from the parenchyma images by PyRadiomics. Third, a novel lung radiomics combination vector (3 + 106) is constructed by the GLM and Lasso algorithm for determining the radiomics risk factors (K = 3) and radiomics node features (d = 106). Last, the COPD stage is classified based on the AMGNN. The results show that compared with the convolutional neural networks and machine learning models, the AMGNN based on constructed novel lung radiomics combination vector performs best, achieving an accuracy of 0.943, precision of 0.946, recall of 0.943, F1-score of 0.943, and ACU of 0.984. Furthermore, it is found that our method is effective for COPD stage classification.


Introduction
As a common and non-infectious lung disease, chronic obstructive pulmonary disease (COPD) presents a preventable, treatable, and progressive chronic disease with debilitating lung conditions characterized by persistent airflow limitation [1].Severe COPD can cause chronic morbidity and eventually lead to death.In 2030, it will become the third-largest death factor worldwide [2].
Because of this characterization, the COPD stage is diagnosed from stage 0 to IV according to Global Initiative for Chronic Obstructive Lung Disease (GOLD) criteria accepted by the American Thoracic Society and the European Respiratory Society [1,3].The assessment parameters in the GOLD criteria are the forced expiratory volume in 1 s/forced vital capacity (FEV 1 /FVC) and FEV 1 % predicted [1,4].The FEV 1 /FVC and FEV 1 % predicted can explain the impact on symptoms and life quality of COPD patients [5,6], but they cannot reflect the change of the lung tissue in COPD patients with COPD stage evolution.The predicted FEV 1 /FVC and FEV 1 % changes only occur when lung tissue is destroyed to a certain extent.However, the measurement accuracy of PFT is limited by the compliance degree of patients.Specifically, the measurement process of PFT is very complex, and it is difficult for patients to understand and comply with the requirements put forward by doctors [7].In addition, PFT cannot intuitively provide detailed anatomical information and morphological changes, such as subtypes of emphysema and bronchial wall thickening [8,9].
Compared with the GOLD criteria and other imaging equipment, computed tomography (CT) has been regarded as the most effective modality for characterizing and quantifying COPD [10].For example, compared with PFT, chest CT images can indicate the patients have suffered from mild lobular central emphysema and decreased exercise tolerance in smokers without airflow limitation [11].In addition, with the significant progress of CT imaging, especially high-resolution CT (HRCT), it has become an effective method for quantitative analysis of COPD, such as measuring the severity of air trapping, emphysema, airflow obstruction, and airway diseases [12,13].However, the quantitative analysis of the bronchial and vascular flow are limited by the resolution of the HRCT.Therefore, it is challenging to automatically, semi-automatically, or manually segment the tiny trachea (small airway) and vascular from chest CT images.In particular, small airways (diameter <2 mm) and associated vessels can hardly be observed from chest CT images.Based on the above, parametric response mapping (PRM) [14] was proposed to locate the small airway lesion region, the emphysema region, and the healthy lung region.First, PRM needs to register the expiratory and inspiratory chest CT images.Then, the mall airway lesion, the emphysema, and the healthy lung region are located by the set thresholds.Therefore, the registration method directly affects the location of the above regions.
Radiomics features in lung disease imaging have been regarded as the state of the art for clinicians [15].However, radiomics features in COPD develop slower than in other lung diseases, such as lung cancer and pulmonary nodules.The diffuse distribution of COPD in the lung limits the application of radiomics features in COPD.Radiomics features should be extracted from the region of interest (ROI) of the chest CT images.However, the diffuse distribution of COPD makes it difficult to determine ROI.Until 2020, Refaee T. et al. point out that radiomics features have not been extensively investigated yet in COPD [16].There are potential applications of radiomics features in COPD for the diagnosis, treatment, and follow-up of COPD and future directions [16].Meanwhile, the value of lung radiomics features in COPD assessment has also been confirmed [17].
Lung radiomics features extracted from the peripheral airway, pulmonary parenchyma, and pulmonary vessels in the chest CT images are suitable for reflecting the change of the lung tissue in COPD patients with COPD stage evolution.Specifically, the characteristic pathological changes of COPD exist in the central airway (trachea, bronchus, bronchioles with an inner diameter greater than 2-4 mm), peripheral airway (bronchioles and bronchioles with an inner diameter less than 2 mm), and pulmonary parenchyma and pulmonary vessels.Chronic inflammation causes the airway wall to repeatedly be damaged and repaired as the COPD stage evolves, resulting in airway blockage and a narrow air cavity in the peripheral airway [18].In addition, the destruction of pulmonary parenchyma in patients with COPD involves the expansion and collapse of respiratory bronchioles.With COPD stage evolution, this expansion and collapse of respiratory bronchioles spread from the upper region to the whole lung, and the pulmonary capillary beds destroy [19].In the earlier COPD stages, the changes in pulmonary vessels are characterized by the thickening of the vascular wall.With the continuous development of COPD, the increase of smooth muscle, proteoglycan, and collagen further thickens pulmonary vessels' walls, which may lead to cor pulmonale [20].Therefore, COPD results from the joint action of the peripheral airway, pulmonary parenchyma, and pulmonary vessels.Thus, the peripheral airway, pulmonary parenchyma, and pulmonary vessels as ROI [17] to extract lung radiomics features are reasonable for COPD stage classification.
Currently, radiomics features have also been used in COPD for survival prediction [21,22], spirometric assessment of emphysema presence and severity [23], COPD exacerbations [24], COPD early decision [3], COPD stage classification [25,26], COPD prediction [27,28], and analysis of COPD and resting heart rate [29].The convolutional neural networks (CNN) and machine learning (ML) models can implement the COPD stage classification task.Compared with the CNN based on chest HRCT images, the multi-layer perceptron (MLP, a kind of ML model) classifier performs better, achieving an accuracy of 0.83, precision of 0.83, recall of 0.83, F1-score of 0.82, and AUC of 0.95 [25].However, the classification performance needs to be further improved.The graph neural network (GNN) was first proposed in 2005 [30].The GNN also developed into the graph convolutional network (GCN) inspired by CNN [31].However, the fixed-graph structure of the GCN by using the entire dataset limits its application and development.To overcome the limitation of the GCN and maintain the advantages of the GNN, an auto-metric graph neural network (AMGNN) based on a meta-learning strategy is proposed [32].The proposed AMGNN has been applied to the Alzheimer's disease classification, achieving a good performance [32,33].However, the risk factors and node features are critical for the classification performance of the AMGNN.Compared with node features, the risk factors with a lower dimension.The selection of the risk factors may have an important impact on the classification effect.Therefore, we focus on COPD and construct a novel lung radiomics combination vector based on the AMGNN for COPD stage classification.
Lung radiomics features have been applied to COPD stage classification, and compared to the CNN models, the ML models with lung radiomics features selected by least absolute shrinkage and selection operator (Lasso) algorithm perform better [25].However, the risk factors and node features of the AMGNN limit the application in COPD stage classification.Therefore, this paper applies the AMGNN with lung radiomics features to improve COPD stage classification performance.Our contributions in this paper are briefly described as follows: (1) This paper proposes a lung radiomics features selection method for risk factors and node features of the AMGNN with the advantages of modeling the correlation between samples and building a small graph structure to classify the COPD stage.The lung radiomics features are only extracted from routine chest HRCT images, eliminating the limitation that risk factors and node features based on prior knowledge are difficult to obtain, such as gene information.The Lasso algorithm with 10-fold cross-validation selects the independent variables related to the dependent variables, determining the critical, independent variables to simplify the classification model.The Lasso algorithm for improving COPD classification has been confirmed [25].The Cox model [3] and the generalized linear model (GLM) are often used for determining the risk factors.Compared with the Cox model, the GLM eliminates the limitation of follow-up features with multiple time series.Last, a novel lung radiomics combination vector (3 + 106) is constructed by GLM and Lasso algorithm for determining the radiomics risk factors (K = 3) and radiomics node features (d = 106) of the AMGNN; (2) Compared with previous work on COPD identification (binary classification: COPD and without COPD) [34,35], our work (COPD stage classification) has more clinical significance.Therefore, our proposed model eliminates the limitations of PFT and may become an effective tool for COPD management.

Materials and Methods
Materials and methods are described in detail in Sections 2.1 and 2.2, respectively.

Materials
Figure 1 shows the participants'selection flow diagram and GOLD distribution of the participants in this study.Specifically, Figure 1a shows the participants' selection flow.Our study cohort was enrolled in the national clinical research center of respiratory diseases, China, from 25 May 2009, to 11 January 2011.Four hundred and sixty-five Chinese participants aged 40-49 were included in this study after being strictly selected by the inclusion and exclusion criteria [36].Then, these 465 participants underwent chest HRCT scans (TOSHIBA, kVp: 120 kV, X-ray tube current: 40 mA, slice thickness: 1.0 mm) at the deep inspiration state and PFT on the same day.The COPD stage 0-III-IV (GOLD 0-III-IV) is diagnosed by GOLD 2008 using FEV 1 /FVC and FEV 1 % predicted PFT [1,4].Figure 1b shows that our study cohort has 129, 108, 121, and 107 participants in the GOLD 0, I, II, and III-IV, respectively.The 465 participants are divided into the train set (70%) and the test set (30%).Furthermore, Figure 1c,d shows the detailed training and test sets in each COPD stage.

Materials
Figure 1 shows the participants'selection flow diagram and GOLD distribution of the participants in this study.Specifically, Figure 1a shows the participants' selection flow.Our study cohort was enrolled in the national clinical research center of respiratory diseases, China, from 25 May 2009, to 11 January 2011.Four hundred and sixty-five Chinese participants aged 40-49 were included in this study after being strictly selected by the inclusion and exclusion criteria [36].Then, these 465 participants underwent chest HRCT scans (TOSHIBA, kVp: 120 kV, X-ray tube current: 40 mA, slice thickness: 1.0 mm) at the deep inspiration state and PFT on the same day.The COPD stage 0-III-IV (GOLD 0-III-IV) is diagnosed by GOLD 2008 using FEV1/FVC and FEV1% predicted PFT [1,4].Figure 1b shows that our study cohort has 129, 108, 121, and 107 participants in the GOLD 0, I, II, and III-IV, respectively.The 465 participants are divided into the train set (70%) and the test set (30%).Furthermore, Figure 1c,d shows the detailed training and test sets in each COPD stage.
This study was approved by the ethics committee in the national clinical research center of China's respiratory diseases.In addition, all participants have provided written informed consent to the first affiliated hospital of Guangzhou Medical University before chest HRCT scans and PFT.

Methods
Figure 2 shows the detailed flow chart of our proposed method in this study.

Lung Parenchyma Segmentation and Radiomics Feature Extraction
Figure 2A(a) shows that a state-of-the-art ResU-Net (U-net (R231)) [37] trained by data diversity (diverse lung disease images), which has been a robust and standard segmentation model of pathological lungs, is transferred to segment lung parenchyma images from 465 sets of chest HRCT images.The network architecture of the ResU-Net has been described in detail in our previous study [38].Then, Figure 2A(b) shows that 1316 lung radiomics features of each participant are extracted from lung parenchyma images with the Hounsfield unit [39] by PyRadiomics [40].
Figure 2B shows the lung radiomics feature extraction model: PyRadiomics.Specifically, two steps are performed to extract the lung radiomics features from lung parenchyma images.The two steps include (1) original lung parenchyma images are filtered by wavelet and Laplacian of Gaussian (LoG) filters, generating derived lung parenchyma images; (2) the original and derived lung parenchyma images are further used to calculate This study was approved by the ethics committee in the national clinical research center of China's respiratory diseases.In addition, all participants have provided written informed consent to the first affiliated hospital of Guangzhou Medical University before chest HRCT scans and PFT.

Methods
Figure 2 shows the detailed flow chart of our proposed method in this study.

Lung Parenchyma Segmentation and Radiomics Feature Extraction
Figure 2A(a) shows that a state-of-the-art ResU-Net (U-net (R231)) [37] trained by data diversity (diverse lung disease images), which has been a robust and standard segmentation model of pathological lungs, is transferred to segment lung parenchyma images from 465 sets of chest HRCT images.The network architecture of the ResU-Net has been described in detail in our previous study [38].Then, Figure 2A(b) shows that 1316 lung radiomics features of each participant are extracted from lung parenchyma images with the Hounsfield unit [39] by PyRadiomics [40].
Figure 2B shows the lung radiomics feature extraction model: PyRadiomics.Specifically, two steps are performed to extract the lung radiomics features from lung parenchyma images.The two steps include (1) original lung parenchyma images are filtered by wavelet and Laplacian of Gaussian (LoG) filters, generating derived lung parenchyma images; (2) the original and derived lung parenchyma images are further used to calculate the lung radiomics features based on the preset classes.More detailed descriptions can be found in our previous studies [25,29].

Radiomics Feature Combination
The risk factors and node features are critical for the classification performance of the AMGNN.Therefore, a novel lung radiomics combination vector (3 + 106) is constructed by the GLM [41] and the Lasso algorithm [42] for determining the radiomics risk factors (K = 3) and radiomics node features (d = 106) of the AMGNN.

Radiomics Feature Combination
The risk factors and node features are critical for the classification performance of the AMGNN.Therefore, a novel lung radiomics combination vector (3 + 106) is constructed by the GLM [41] and the Lasso algorithm [42] for determining the radiomics risk factors (K = 3) and radiomics node features (d = 106) of the AMGNN.
The results of a previous study [25] have shown that the Lasso algorithm helps improve the classification effect of COPD.Specifically, Figure 2A(c) shows that the Lasso algorithm with 10-fold cross-validation selects the radiomics node features from 1316 lung radiomics features.Meanwhile, the radiomics risk factors are selected from 1316 lung radiomics features by calculating the top three maximal R 2 values generated from the GLM.Finally, the radiomics risk factors (K = 3) and radiomics node features (d = 106) are concatenated, obtaining the proposed lung radiomics combination vector (3 + 106) for COPD classification.
where the link function g(µ y i ) = η relates the mean µ y i = E(y i ) to the linear predictor . y i denotes the dependent variable (the 465 COPD stages: GOLD 0, I, II, and III-IV).x j denotes the independent variable (the 465 × 1316 normalized lung radiomics features).β j denotes the regression coefficients i∈ [1, n] and j∈[0, p].

COPD Stage Classification Based on AMGNN
Figure 2A(d) shows that the COPD stage is classified based on AMGNN with the proposed lung radiomics combination vector.Specifically, Figure 2C shows the pipeline of AMGNN based on meta-learning.Specifically, T graphs are separately established by selecting the 40 known nodes in the training set and 1 unknown node in the training set.The 40 (10 × 4) known nodes of each graph include 10 known nodes in each COPD stage, and four COPD stages are included in our study.Subsequently, T-1 meta tasks train the AMGNN by randomly to the constantly updated network parameter P and loss value , establishing the AMGNN (P T ).Lastly, the established AMGNN (P T ) realizes the unknown nodes in the test set for COPD stage classification.
In addition, Figure 2D shows the detailed network structure of AMGNN (P i ) for calculating the auto-metics connection with probability constraints and updating the node.Specifically, two stages are parallelly performed in each AMGNN (P i ) for generating the adjacency matrix A. First, the adjacency matrix A multiplies the proposed lung radiomics combination vector of the nodes in the l th layer V (l) , and then updating the node Gn(V (l) ) is obtained by a nonlinear activation function FC (Leaky ReLU).Finally, the output V (l+1) of the AMGNN (P i ) is obtained by concatenating V (l) and Gn(V (l) ).However, the output of the last AMGNN (P T-1 ) is input into the softmax instead of the concatenating operation.
Furthermore, one stage is constructing the edge weight matrix W using the 106 radiomics nodes features.Meanwhile, the other stage is constructing the edge constraint matrix E using the 3 radiomics risk factors.Finally, the edge weight matrix W element wisely multiplies the edge constraint matrix E, obtaining the adjacency matrix A. In constructing the W stage, the N × N × 106 feature block C is generated by copying N times of the N × 106 feature map.Then, each N × N × 1 in the N × N × 106 feature block C is transposed, generating a new N × N × 106 feature block C .The difference feature block C is further obtained by calculating the absolute difference between each feature of the two nodes in C and C .Subsequently, a CNN with a 1 × 1 kernel generates the edge weight matrix W. The detailed description of calculating the edge constraint matrix E refers to the study [32].

Experiments
Figure 3 shows four experiments to verify the effectiveness of our proposed method.The radiomics/CNN combination vectors (K risk factors + d node features) based on the AMGNN for COPD classification in Experiments 2-4 can be obtained in Supplementary Materials.The open-source code of AMGNN with 5-fold cross-validation is directly applied to COPD stage classification.Similarly, other ML models are also trained with the 5-fold cross-validation.Finally, the trained models with the best AUC are loaded and used on the test set.
The CNN models based on the chest HRCT images perform unsatisfactory classification [25].However, the model fusion strategy based on the CNN and ML models should be further studied.Specifically, the features are extracted by the CNN model based on transfer learning, and then the ML classifiers are used for classification.Therefore, 3D CNN features are extracted based on the encoder backbone (ResNet 10) of a pre-trained Med3d [41] using truncated transfer learning for comparing the lung radiomics features.Med3d, a heterogeneous 3D network, is used to extract general medical 3D features by building a 3DSeg-8 dataset with diverse modalities, target organs, and pathologies.Five hundred and twelve 3D CNN feature maps with the size of 3 × 3 × 3 are generated by ResNet 10.Therefore, each participant has 13,824 3D CNN features (512 × 3 × 3 × 3 = 13,824) by flattening the 512 3D CNN feature maps.
Experiments 1 is designed based on ML classifiers with lung radiomics or 3D CNN features.The support vector machine (SVM) [42], multi-layer perceptron (MLP) [25,43], random forest (RF) [44], logistic regression (LR) [45], and the proposed AMGNN are applied to diagnose Alzheimer's disease [32].In addition, SVM, LR, MLP, RF, gradient boosting (GB) [46], and linear discriminant analysis (LDA) [47] are applied to classify the COPD stages [25].Therefore, the above six classic ML classifiers and the AMGNN are used for COPD stage classification.Specifically, the original lung radiomics (1316) and 3D CNN features (13,824) are directly and separately used to classify the COPD stage based on ML classifiers.In addition, the radiomics features (106) and 3D CNN features (60) of each participant are automatically selected by the Lasso algorithm with the 10-fold cross-validation.For comparison of the same number of features, the radiomics features (106) and 3D CNN features (60) of each participant are separately selected or fused by the GLM/principal component analysis (PCA) algorithm.These selected or fused features are separately used to classify the COPD stage based on ML classifiers.However, the combination feature vector of the lung radiomics or 3D CNN features selected by GLM/Lasso and these features separately fused by PCA are also considered in our experimental design in Experiment 1. Specifically, because the GLM and Lasso algorithm perform the feature selection task, the features selected by the GLM and Lasso algorithm are not further combined.
Experiments 2-4 are designed based on the AMGNN with lung radiomics or 3D CNN features for determining the risk factors and node features.Specifically, compared with node features of AMGNN, the risk factors have a lower dimension.Therefore, the PCA algorithm and GLM are separately used to determine the risk factors (K = 2-6).In addition, the PCA algorithm and GLM are separately used to determine the risk factors (K = 2-6) based on the original lung radiomics (1316) or 3D CNN features (13,824), and the corresponding node features are original lung radiomics (1316) or 3D CNN features (13,824) in Experiment 2. The PCA and GLM are separately used to determine the risk factors (K = 2-6) based on the selected lung radiomics (106) or 3D CNN features (60) by the Lasso algorithm, and the corresponding node features also are the selected lung radiomics (106) or 3D CNN features (60) in Experiment 3. The PCA and GLM are separately used to determine the risk factors (K = 2-6) based on 1316 original lung radiomics or 13,824 3D CNN features, and the corresponding node features are the lung radiomics (106) or 3D CNN features (60) selected by the Lasso algorithm in Experiment 4.  Table 1 reports the different ML classifiers with their definitions.The 4-way-10-shot of the few-shot learning [48] is set in the AMGNN (iterations = 600, batch_size_test = 28, batch_size_train = 28, and random_seed = 42).All the source codes run on PyCharm Table 1 reports the different ML classifiers with their definitions.The 4-way-10-shot of the few-shot learning [48] is set in the AMGNN (iterations = 600, batch_size_test = 28, batch_size_train = 28, and random_seed = 42).All the source codes run on PyCharm 2020.3.5 (professional edition) on Windows 10 Pro 64-bit with two 2080 Ti GPUs, 32 GB RAM, 1 TB mechanical storage, and a 256 G SSD.

Results
This section reports the results of Experiments 1-4, including the five standard evaluation metrics (accuracy, precision, recall, F1-score, and area under the curve (AUC)).The evaluation metric AUC for multi-classification is calculated by the receiver operating characteristic curve (ROC) [25].Compared with the CNN and ML models, the AMGNN based on the novel lung radiomics combination vector (3 + 106) performs best, achieving an accuracy of 0.943, precision of 0.946, recall of 0.943, F1-score of 0.943, and ACU of 0.984.

COPD Stage Classification Based on Different ML Classifiers
Tables 2-7 report the ML classifiers' performance in Experiment 1.Meanwhile, Figures 4 and 5 visually show the evaluation metrics and ROC curves of the ML classifier in Experiment 1.The MLP classifier performs better than the other ML classifiers for COPD stage classification.Specifically, the MLP classifier with 13,824 original 3D CNN/1316 original lung radiomics features performs best, achieving an accuracy of 0.793/0.786,precision of 0.798/0.784,recall of 0.793/0.784,F1-score of 0.790/0.784,and ACU of 0.938/0.919.In addition, 60 3D CNN/106 lung radiomics features separately selected or fused by the Lasso algorithm, GLM, and PCA algorithm are further used to classify the COPD stage.Again, the MLP classifier with 60 3D CNN/106 lung radiomics features selected by the Lasso algorithm performs best, achieving an accuracy of 0.821/0.829,precision of 0.826/0.828,recall of 0.821/0.829,F1-score of 0.821/0.824,and ACU of 0.946/0.950.Even for the CNN/radiomics combination vector, the MLP classifier with 60 3D CNN/106 lung radiomics features only selected by the Lasso algorithm also performs best.The GLM fails to improve the classification performance of all ML classifiers with 3D CNN features.Compared with the GLM, the PCA algorithm only improves the classification performance of the RF classifier with 3D CNN features.The Lasso algorithm only improves the classification performance of the MLP classifier with the 3D CNN features.In addition, the GLM fails to improve the classification performance of the SVM classifier, but it improves the classification performance of other ML classifiers with lung radiomics features.Compared with the GLM, the PCA algorithm only improves the classification performance of the RF and LDA classifiers with lung radiomics features.However, the Lasso algorithm improves the classification performance of all ML classifiers with lung radiomics features.
Compared to the evaluation metrics of the MLP classifier with 13,824 original 3D CNN features, that of the MLP classifier with 60 selected 3D CNN features has been improved by 2.8% in accuracy, 2.8% in precision, 2.8% in recall, 3.1% in F1-score, and 0.8% in ACU.Similar to 3D CNN features, compared to the evaluation metrics of the MLP classifier with 1316 original lung radiomics features, that of the MLP classifier with 106 selected lung radiomics features has been improved by 4.3% in accuracy, 4.4% in precision, 4.5% in recall, 4.0% in F1-score, and 3.1% in ACU.

COPD Stage Classification Based on the AMGNN Classifier
Tables 8-10 report the AMGNN classifier's performance in Experiments 2-4.Meanwhile, Figures 6 and 7 visually show the evaluation metrics and ROC curves of the AMGNN classifier in Experiments 2-4.The AMGNN classifier with the proposed lung radiomics combination vector constructed by 3 radiomics risk factors (selected from 1316 original lung radiomics features by GLM) and 106 radiomics node features (selected from 1316 original lung radiomics features by the Lasso algorithm) performs best for COPD stage classification, achieving an accuracy of 0.943, precision of 0.946, recall of 0.943, F1-score of 0.943, and ACU of 0.984.
Tables 8-10 and Figures 6 and 7 also show that the evaluation metrics of the AMGNN classifier with lung radiomics combination vector are better than the AMGNN classifier with the 3D CNN combination vector.
Table 11 compares the results of our proposed method with other previous methods.Furthermore, it fully demonstrates the superiority of our method over other previous methods.

Discussion
Based on the experimental results, we give the following discussion and point out the limitations in this study and the future direction.Compared with multimodal sleep data [49] and electromyography [50], lung radiomics features extracted from chest HRCT images are more suitable for COPD stage classification.
First, the MLP classifier with 3D CNN or lung radiomics features performs best for the COPD stage classification in the ML classifiers.The MLP classifier's structure and ability to handle complex nonlinear features determine its excellent ability.The MLP classifier is composed of three full connection layers, which is more efficient and more suitable for modeling long-range dependencies [51].Meanwhile, COPD patients have high heterogeneity and different phenotypes [1], resulting in complex nonlinear 3D CNN or lung radiomics features extracted from their chest HRCT images.However, the MLP classifier can handle these complex nonlinear features and discover dependencies between different input features by approximating the nonlinear map globally to realize the COPD stage classification [25,43].
Second, the AMGNN classifier performs better than the best MLP classifier in the COPD stage classification.This is because the AMGNN classifier, which overcomes the lack of flexibility of existing GNN models, introduces a meta-learning strategy, and is insensitive to graphics size [32].In addition, the AMGNN classifier needs fewer training samples than traditional ML classifiers while maintaining stability in reducing training samples [32].Therefore, even if a small number of graphs are used, it also can show good performance for the COPD stage classification.The AMGNN classifier also solves the problems that the standard medical image cohort is usually challenging to obtain, which results in a small data volume.On the other hand, the MPL classifier needs sufficient data to train its parameters.Therefore, the MPL classifier's performance with small data may deteriorate very seriously.Meanwhile, the AMGNN classifier performs well on the test set, showing no underfitting and overfitting problems.Specifically, too deep/complex network structure often causes underfitting.However, the AMGNN classifier only with two layers effectively avoids underfitting.In addition, the AMGNN classifier introduces the meta-learning strategy.Meta-learning-based approaches can increasingly train powerful CNNs on small datasets in many vision problems [52].Therefore, the meta-learning strategy effectively avoids overfitting.
Third, the Lasso algorithm improves the MLP and AMGNN classifiers' performance.Compared with the Lasso algorithm, the PCA algorithm reduces the original features' dimension by identifying the orthogonal linear combinations with the largest covariance.Therefore, the PCA algorithm can compress the original features into a small number of features and obtain new fused features without losing too much information.However, because the PCA algorithm needs to retain most of the features, the fused features still affect the COPD classification effect.The Lasso algorithm and GLM perform the feature selection task.Compared with the GLM, the Lasso algorithm is a penalized likelihood approach.Therefore, the Lasso algorithm automatically selects the original features, obtaining the fixed numbers of the selected features.However, GLM generates R 2 values corresponding to each feature instead of the fixed numbers of the selected features.In addition, the Lasso algorithm is often used with survival analysis models to determine variables and eliminate the collinearity problem between variables [3].The Lasso algorithm is also applied to select the features to improve the MLP classifier's performance by establishing the relationship between the independent variables (3D CNN or lung radiomics features) and dependent variables (the COPD stages) [25].Therefore, the complexity of the 3D CNN or lung radiomics features is reduced.As a result, the MLP classifier can focus on the 60 selected 3D CNN or 106 selected lung radiomics features to improve the classifier's performance.The Lasso algorithm also applies to the AMGNN classifier to construct the excellent edge weight matrix by reducing redundant collinearity 3D CNN or lung radiomics features.Meanwhile, the Lasso algorithm determines node features of the AMGNN classifier, eliminates collinearity features, and further avoids overfitting.
Fourth, lung radiomics features remove the limitation of risk factors and node features that are difficult to obtain.Meanwhile, the AMGNN classifier with the lung radiomics combination vector performs better than the 3D CNN combination vector.The risk factors of diseases, such as the risk factors (age, gender, year of education, and APOe4 gene information) of Alzheimer's disease [32], are hardly obtained or determined, which limits the application of the AMGNN classifier.However, because a large number of lung radiomics features can be extracted from the ROI of chest HRCT images, we may overcome the limitation of obtaining and determining risk factors.Compared with the 3D CNN features, the lung radiomics features are calculated by preset formulas, which are easier to explain in subsequent studies.
Lastly, this study also has some limitations, and we point out the future direction.Although the AMGNN classifier with the proposed novel lung radiomics combination vector achieves promising results, the model must be retrained with a newcomer to be classified.Furthermore, we only explained the method from the engineering and algorithm, but the clinical significance of three radiomics risk factors needs further analysis by COPD experts and doctors.However, the combination vector of 3D CNN and lung radiomics based on ML or AMGNN classifiers for the COPD stage classification should be studied in the future.

Conclusions
This paper constructs a novel lung radiomics combination vector generated by GLM and Lasso algorithm for the COPD stage classification based on the auto-metric graph neural network.Compared with the CNN and ML models, the AMGNN based on the novel lung radiomics combination vector (K = 3 (GLM) + d = 106, 3 radiomics risk factors are selected by GLM and 106 radiomics node features are selected by the Lasso algorithm) performs best, achieving an accuracy of 0.943, precision of 0.946, recall of 0.943, F1-score of 0.943, and ACU of 0.984.Therefore, our proposed model eliminates the limitations of PFT and may become an effective tool for COPD management.

Patents
The method and device, electronic device and storage medium for stage classification of chronic obstructive pulmonary disease, CN2022104685981, Shenzhen Technology University and Northeastern University, China.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/diagnostics12102274/s1.Informed Consent Statement: Written informed consent has been obtained from the patients to publish this paper.

Figure 1 .
Figure 1.The participants' selection flow diagram and GOLD distribution of the participants in this study.(a) The participants' selection flow diagram; (b) GOLD distribution in our study cohort; (c) training set distribution; (d) test set distribution.

Figure 1 .
Figure 1.The participants' selection flow diagram and GOLD distribution of the participants in this study.(a) The participants' selection flow diagram; (b) GOLD distribution in our study cohort; (c) training set distribution; (d) test set distribution.

Figure 2 .
Figure 2. The proposed method in this study.(A) Our proposed method: COPD stage classification using auto-metric graph neural network; specifically, our proposed method includes (a) lung parenchyma segmentation, (b) radiomics feature extraction, (c) radiomics feature combination, and

Figure 2 .
Figure 2. The proposed method in this study.(A) Our proposed method: COPD stage classification using auto-metric graph neural network; specifically, our proposed method includes (a) lung parenchyma segmentation, (b) radiomics feature extraction, (c) radiomics feature combination, and (d) COPD stage classification.(B) Lung radiomics feature extraction model: PyRadiomics [25,29,33,40].(C) The pipeline of the AMGNN based on meta-learning.(D) Detailed network structure of the AMGNN (P i ).

Figure 3 .
Figure 3. Experimental design in this paper.

Figure 3 .
Figure 3. Experimental design in this paper.

Figure 4 .
Figure 4. Visual evaluation metrics of different ML classifiers with different features in Experiment 1. (a,b) The evaluation metrics of different ML classifiers with 13,824 original 3D CNN/1316 original lung radiomics features.(c-h) The evaluation metrics of different ML classifiers with 60 3D CNN/106 lung radiomics features separately selected or fused by the Lasso algorithm, GLM, and PCA algorithm.(i-l)The evaluation metrics of different ML classifiers with different combinations of feature vectors of the selected and fused lung radiomics features or the selected and fused 3D CNN features.

Diagnostics 2022 ,Figure 5 .
Figure 5. Visual ROC curves of different ML classifiers with different features in Experiment 1. (a,b) The ROC curves of different ML classifiers with 13,824 original 3D CNN/1316 original lung radiomics features.(c-h) The ROC curves of different ML classifiers with 60 3D CNN/106 lung radiomics features separately selected or fused by the Lasso algorithm, GLM, and PCA algorithm.(i-l) The ROC curves of different ML classifiers with different combinations of feature vectors of the selected and fused lung radiomics features or the selected and fused 3D CNN features.

Figure 5 .
Figure 5. Visual ROC curves of different ML classifiers with different features in Experiment 1. (a,b) The ROC curves of different ML classifiers with 13,824 original 3D CNN/1316 original lung radiomics features.(c-h) The ROC curves of different ML classifiers with 60 3D CNN/106 lung radiomics features separately selected or fused by the Lasso algorithm, GLM, and PCA algorithm.(i-l) The ROC curves of different ML classifiers with different combinations of feature vectors of the selected and fused lung radiomics features or the selected and fused 3D CNN features.

K = 3 (Figure 6 .
Figure 6.Visual evaluation metrics and ROC curves of the AMGNN classifier in Experiments 2-4.(a-c) The mean evaluation metrics, and (d-f) the best evaluation metrics.Figure 6. Visual evaluation metrics and ROC curves of the AMGNN classifier in Experiments 2-4.(a-c) The mean evaluation metrics, and (d-f) the best evaluation metrics.

Figure 6 .Figure 7 .
Figure 6.Visual evaluation metrics and ROC curves of the AMGNN classifier in Experiments 2-4.(a-c) The mean evaluation metrics, and (d-f) the best evaluation metrics.Figure 6. Visual evaluation metrics and ROC curves of the AMGNN classifier in Experiments 2-4.(a-c) The mean evaluation metrics, and (d-f) the best evaluation metrics.Diagnostics 2022, 12, x FOR PEER REVIEW 17 of 23

Figure 7 .
Figure 7. ROC curves of the AMGNN classifier in Experiments 2-4.(a-f) ROC curves of the AMGNN classifiers with different combination vectors.

Author Contributions:
Conceptualization, R.C. and Y.K.; methodology, Y.Y., S.W. and N.Z.; software, Y.Y. and S.W.; validation, N.Z., W.D. and Z.C.; formal analysis, Y.Y., S.W., W.L. and N.Z.; investigation, Y.Y., S.W., Y.G. and H.C.; resources, H.C., X.L. and R.C.; data curation, H.C., X.L. and R.C.; writingoriginal draft preparation, Y.Y.; writing-review and editing, Y.L., W.D. and Y.G.; visualization, Y.Y., S.W., Y.G. and N.Z.; supervision, R.C. and Y.K.; project administration, R.C. and Y.K.; funding acquisition, Y.K., W.L. and H.C. All authors have read and agreed to the published version of the manuscript.Funding: This research was funded by the National Natural Science Foundation of China, grant number 62071311; the Stable Support Plan for Colleges and Universities in Shenzhen of China, grant number SZWD2021010; the Scientific Research Fund of Liaoning Province of China, grant number JL201919; the Natural Science Foundation of Guangdong Province of China, grant number 2019A1515011382; the special program for key fields of colleges and universities in Guangdong Province (biomedicine and health) of China, grant number 2021ZDZX2008.Institutional Review Board Statement: All patients provided written informed consent, and this study was approved by the Guangzhou medical university Ethics Committee (grant number: 2009-09, Approval Date: 5 August 2009) and registered at http://www.chictr.org.cn(registration number: ChiCTR2000034586, accessed on 11 July 2020).

Table 1 .
The different ML classifiers with their definitions.

Table 2 .
The different ML classifiers' performance on the test set of original 3D CNN features (13,824)/lung radiomics features (1316) in Experiment 1.

Table 3 .
The different ML classifiers' performance on the test set of selected 3D CNN features (60) or lung radiomics features (106) generated by the Lasso algorithm in Experiment 1.

Table 4 .
The different ML classifiers' performance on the test set of selected 3D CNN features (60) or lung radiomics features (106) generated by the GLM in Experiment 1.

Table 5 .
The different ML classifiers' performance on the test set of fused 3D CNN features (60) or lung radiomics features (106) generated by the PCA algorithm in Experiment 1.

Table 6 .
The different ML classifiers' performance on the test set of the CNN combination vector (Lasso + PCA/GLM + PCA) in Experiment 1.

Table 7 .
The different ML classifiers' performance on the test set of lung radiomics combination vector (Lasso + PCA/GLM + PCA) in Experiment 1.

Table 8 .
The AMGNN classifier's performance on the test set in Experiment 2.

Table 9 .
The AMGNN classifier's performance on the test set in Experiment 3.

Table 10 .
The AMGNN classifier's performance on the test set in Experiment 4.

Table 11 .
Comparison of the results of our proposed method with other previous methods.