Prediction model for radiation pneumonitis based on computed tomography

Abstract

The invention relates to the technical field of medical treatment, and particularly discloses a radiation pneumonitis prediction model based on computed tomography, which comprises the following steps of: defining an intra-pulmonary region receiving radiation dose of more than 20Gy and 30Gy on a CT image (CT1) as V20 and V30 as a region of interest, extracting radiomics characteristics and dosiomics characteristics in the region of interest, and calculating the radiation pneumonitis in the region of interest according to the radiomics characteristics and the dosiomics characteristics; the method comprises the following steps: selecting patients subjected to 40Gy-50Gy radiotherapy and subjected to a field contraction plan CT image (CT2) from the patients, extracting radiomics characteristics of regions of interest V20 and V30 in the field contraction plan CT image (CT2), obtaining delta-radiomics characteristics by using delta-RF = RFCT2-RFCT1, evaluating the performance of a model by using an area under an ROC curve, calculating a confidence interval by using statistics and a bootstrap method, and calculating the reliability of the model according to the confidence interval. The significant difference of ROC curve areas of different models is detected by using Delon detection, and the pneumonia is predicted by constructing the models in a manner of combining three omics, namely radiomics, dosiomics and delta-radiomics.

Description

Technical Field

[0001] The present invention relates to the field of medical technology, and specifically to a prediction model for radiation pneumonitis based on computed tomography. Background Art

[0002] Radiation pneumonitis is a lung inflammation caused by radiation exposure to the lungs. Patients usually present with symptoms such as cough, sputum production, chest tightness, and shortness of breath. Severe patients may also experience symptoms such as hemoptysis and dyspnea. If a patient exhibits the above symptoms, they can be preliminarily diagnosed with radiation pneumonitis.

[0003] However, currently, there are models for predicting radiation pneumonitis constructed by separately applying radiomics, separately applying dosomics, separately applying delta-radiomics, or by combining radiomics and dosomics. There has been no application of a model constructed by combining the three omics methods. Summary of the Invention

[0004] The purpose of the present invention is to provide a prediction model for radiation pneumonitis based on computed tomography to solve the problems raised in the above background art.

[0005] To achieve the above purpose, the present invention provides the following technical solution: A prediction model for radiation pneumonitis based on computed tomography, including the following steps:

[0006] S1: Image segmentation module, respectively collect the initial plan CT image (CT1), reduced field plan CT image (CT2) of the case, and the initial plan dose distribution map, and obtain the corresponding regions of interest;

[0007] S2: Perform feature extraction;

[0008] S3: Perform feature selection;

[0009] S4: Perform model construction, constructing a single model, a combined model, and a nomogram;

[0010] S5: Perform model evaluation;

[0011] S6: Perform model interpretation.

[0012] Preferably, in S2, the lung regions on the CT image (CT1) receiving irradiation doses above 20 Gy and 30 Gy are defined as V20 and V30 as regions of interest, and the radiomics features and dosomics features within the regions of interest are extracted. Among the patients, those who have undergone reduced field plan CT image (CT2) after radiotherapy of 40 Gy - 50 Gy are selected, and the radiomics features of the regions of interest V20 and V30 in the reduced field plan CT image (CT2) are extracted. Using delta-RF = RFCT2 - RFCT1, delta-radiomics features are obtained.

[0013] Preferably, in step S3, the variance selection method, Pearson correlation coefficient, least absolute shrinkage and selection operator are used to screen the features with non-zero feature screening coefficients.

[0014] Preferably, in step S5, the area under the ROC curve is used to evaluate the performance of the model, the confidence interval is calculated using statistics and the bootstrap method, and the DeLong test is used to test the significant difference in the ROC curve area of different models.

[0015] The present invention has at least the following beneficial effects:

[0016] 1. By applying a combination of three omics, namely radiomics, dosomics and delta-radiomics, this solution constructs a model to predict pneumonia. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 It is a schematic diagram of the overall process of the present invention;

[0018] Figure 2 In it, ① is a schematic diagram of image segmentation of the present invention;

[0019] Figure 2 In it, ② is a schematic diagram of feature extraction of the present invention;

[0020] Figure 2 In it, ③ is a schematic diagram of feature selection of the present invention;

[0021] Figure 2 In it, ④ is a schematic diagram of model construction of the present invention;

[0022] Figure 2 In it, ⑤ is a schematic diagram of model evaluation of the present invention;

[0023] Figure 2 In it, ⑥ is a schematic diagram of explaining the model intention of the present invention;

[0024] Figure 3 It is a schematic diagram of ROC training of the present invention, where a and c are the ROCs of different models in the training set, and b and d are the ROCs of different models in the test set. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0025] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0026] Embodiment 1

[0027] Please refer toFigure 1 - Figure 3 , a prediction model for radiation pneumonitis based on computed tomography, as shown in Table 1 (Figure 1), which is a technical flow chart summarizing the key method steps of this study:

[0028] ①. Image segmentation: (1), (2), and (3) are respectively the initial planning CT image (CT1), reduced field planning CT image (CT2), and initial planning dose distribution map of the collected cases, and the corresponding regions of interest (ROI) are obtained;

[0029] ②. Feature extraction: Outline the lung regions (defined as V20 and V30) in the initial planning CT image (CT1) that receive irradiation doses of 20 Gy, 30 Gy, and above as regions of interest, extract the radiomics features and dosomics features within the regions of interest. Among 147 patients, 76 patients who underwent reduced field planning CT image (CT2) after radiotherapy with 40 Gy - 50 Gy were selected, and the radiomics features (Radiomics feature, RF) of regions of interest V20 and V30 in the reduced field planning CT image (CT2) were extracted. Use delta - RF = RFCT2 - RFCT1 to obtain delta - radiomics features.

[0030] ③. Feature selection: Use the variance selection method, Pearson correlation coefficient, and least absolute shrinkage and selection operator (LASSO) to screen features with non - zero feature screening coefficients;

[0031] ④. Model construction: Construct a single model, a combined model, and a nomogram;

[0032] ⑤. Model evaluation: Use indicators such as Area Under the ROC Curve / Area Under the ROC Curve / AUC, accuracy, C - index, F1 - score, sensitivity, specificity, and precision to evaluate the performance of the model. Use statistics and the bootstrap method

[0033] / Bootstrap to calculate the confidence interval. Use the Delong Test to test the significant differences in AUC of different models.

[0034] ⑥. Explain the model: Use Shapley Additive Explanations / SHAP to explain the model.

[0035] Specific detailed methods:

[0036] I. Image segmentation: A retrospective analysis was conducted on 147 lung cancer patients in our hospital (117 in the training set and 30 in the test set). The pulmonary regions (defined as V20 and V30) receiving irradiation doses of 20 Gy, 30 Gy and above on the initial planning CT images (CT1) were semi - automatically outlined as regions of interest. Radiomics features and dosomics features within the regions of interest were extracted. Among the 147 patients, 76 patients (61 in the training set and 15 in the test set) who underwent reduced - field planning CT images (CT2) after radiotherapy with 40 Gy - 50 Gy were selected. The radiomics features (Radiomics feature, RF) of the regions of interest V20 and V30 in the reduced - field planning CT images (CT2) were extracted. The delta - radiomics features were obtained using delta - RF = RFCT2 - RFCT1 (delta represents Δ, meaning the difference before and after; delta - RF represents delta - radiomics feature; RFCT2 represents the radiomics feature from CT2; RFCT1 represents the radiomics feature from CT1).

[0037] II. Feature extraction: A total of 2790 radiomics features, including 186 original image features and 2604 image - filtering - based features, were extracted from CT1 and CT2 using AccuContour3.0 in the Manteia scientific data management system (version 3.1, http: / / www.manteiatech.com). A total of six types of features were extracted, including FirstOrderStatistics (540), Glcm (GrayLevelCooccurenceMatrix, 720), Glszm (GrayLevelSizeZoneMatrix, 480), Glrlm (GrayLevelRunLengthMatrix, 480), Ngtdm (NeigbouringGrayToneDifferenceMatrix, 150), Gldm (GrayLevelDependenceMatrix, 420). Shape, which is a descriptor of the three - dimensional size and shape of the ROI, was not extracted here. A total of 214 dosomics features, including seven feature classes of FirstOrderStatistics (36), Shape (28), Glcm (48), Glszm (32), Glrlm (32), Ngtdm (10), and Gldm (28), were extracted from the 3D dose distribution maps using PyRadiomics (version 3.0.1) in 3Dslicer (version 4.11, an open - source python3.7 package, website: https: / / www.slicer.org).

[0038] III. Feature Screening: Variance selection method, Pearson correlation coefficient (PCC), and least absolute shrinkage and selection operator (LASSO) are used for feature screening.

[0039] IV. Model Building: Eight models are respectively built for each type of omics feature after screening. The models are Logistic regression (LR), least absolute shrinkage and selection operator (LASSO), random forest (RF), k-nearest neighbors (KNN), support vector machine (SVM), extreme gradient boosting (Xgboost), light gradient boosting machine (LightGBM), and gradient boosting decision tree (GBDT). The model with the highest receiver operating characteristic curve is selected. Five-fold cross-validation is used to evaluate the training set. At the same time, the radiomics risk score / Dosiomics risk score, / D-score, the dosomics risk score / Radiomics risk score / R-score, and the delta-radiomics risk score / Delta-radiomics risk score / Delta-Rscore are calculated. A combined model with clinical factors and dosimetric parameters is established (two combined models are established in total, namely the combination of radiomics score, dosomics score, clinical factors, and dosimetric factors (Model A+B+D), and the combination of radiomics score / R-score, dosomics score / D-score, delta-radiomics score / Delta-R-score, clinical factors, and dosimetric factors (Model A+B+C+D)) and a nomogram.

[0040] V. Model Evaluation: Finally, the area under the receiver operating characteristic curve, decision curve analysis, and calibration curve were used to evaluate the performance of the model. The Shapley Additive Explanations (SHAP) was used to interpret the model. The method adopted in a of the above figure is different from previous studies. A prediction model was constructed using omics features extracted from three types of omics. Radiomics: Radiomics features (RFCT1) within the region of interest of the initial planning CT1 and radiomics features (RFCT2) within the region of interest of the reduced field planning CT2 were extracted. Dosomics: Dosomics features within the region of interest of the initial planning dose distribution map were extracted. Delta-radiomics: Delta-radiomics feature = RFCT2 - RFCT1.

[0041] Specifically, the detailed calculation steps involved are as follows:

[0042] I. Variance Selection Method: Two different statistical test methods were used for feature selection: T-test and U-test. The following are the calculation formulas for these two tests:

[0043] T-test (Independent Samples T-test):

[0044] The T-test is used to test whether there is a significant difference in the means of two independent samples. Its calculation formula is as follows:

[0045] [t = \frac{\bar{X}_1 - \bar{X}_2}{\sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}}]

[0046] Where:

[0047] (t) is the T statistic.

[0048] (\bar{X}_1) and (\bar{X}_2) are the means of the two samples respectively.

[0049] (s_1^2) and (s_2^2) are the variances of the two samples respectively.

[0050] (n_1) and (n_2) are the sizes of the two samples respectively.

[0051] U-test (Mann-Whitney U-test):

[0052] The U-test, also known as the Mann-Whitney U-test, is used to test whether two independent samples have the same distribution. Its calculation formula is as follows:

[0053] [U = \sum_{i = 1}^{n_1}\sum_{j = 1}^{n_2}I(X_{ij} < Y_{ik})]

[0054] Where: (U) is the U - statistic.

[0055] (I) is the indicator function, if (X_{ij} < Y_{ik}) is true, then (I = 1), otherwise (I = 0).

[0056] (X_{ij}) and (Y_{ik}) are the elements in two samples respectively.

[0057] II. The Pearson correlation coefficient formula, which is used to calculate the strength of the linear relationship between two variables. The formula is as follows:

[0058] [r=\frac{n(\sum xy)-(\sum x)(\sum y)}{\sqrt{[n\sum x^{2}-(\sum x)^{2}][n\sum y^{2}-(\sum y)^{2}]}}]

[0059] Where: (r) is the correlation coefficient. (n) is the sample size. (x) and (y) are the values of two variables respectively. (\sum xy) is the sum of the products of two variables. (\sum x) and (\sum y) are the sums of each variable. (\sum x^{2}) and (\sum y^{2}) are the sums of the squares of each variable.

[0060] III. The calculation formula of LASSO regression is:

[0061] [\hat{y}=X\beta+\epsilon]

[0062] [\beta=\arg\min_{\beta}\frac{1}{n}\sum_{i = 1}^{n}(y_{i}-X_{i}\beta)^{2}+\alpha\sum_{j = 1}^{p}|\beta_{j}|]

[0063] Where: (\hat{y}) is the predicted value. (X) is the feature matrix. (\beta) is the regression coefficient. (\epsilon) is the error term. (\alpha) is the regularization parameter. (p) is the number of features.

[0064] IV. Select features to build models. Each omics builds a model, including the following: (I). Logistic Regression (LR) is a generalized linear model used to estimate the probability of an event occurring. Its formula is:

[0065] [P(y = 1)=\frac{1}{1 + e^{-\beta_0 - \beta_1X_1 - \beta_2X_2 - \ldots - \beta_nX_n}}]

[0066] The meanings of the symbols are as follows:

[0067] (P(y = 1)) is the probability that the target variable (y) is 1. (X_1, X_2, …, X_n) are the feature variables.

[0068] (\beta_0, \beta_1, …, \beta_n) are the model parameters, also known as weights. (e) is the base of the natural logarithm.

[0069] The meaning of this formula is that given a set of features (X_1, X_2, …, X_n), the logistic regression model will calculate a probability value (P(y = 1)), which is the probability of a certain event occurring. This probability is transformed from the linear combination (\beta_0 + \beta_1X_1 + \beta_2X_2 + … + \beta_nX_n) through a logistic function (also known as the sigmoid function).

[0070] The linear combination part is a linear equation in the form of (\beta_0 + \beta_1X_1 + \beta_2X_2 + …

[0071] + \beta_nX_n). The value of this linear combination determines the input of the sigmoid function and thus determines the value of (P(y = 1)).

[0072] The formula for the sigmoid function is:

[0073] [\sigma(z)=\frac{1}{1 + e^{-z}}]

[0074] where (z) is the value of the linear combination. This function maps any real number to a probability value between 0 and 1.

[0075] (2). The calculation formula of the Random Forest (RF) is not the formula of a single decision tree, but a combination of multiple decision trees.

[0076] For each decision tree in the random forest, its calculation formula is as follows: (1) Decision tree construction: Select a subset of the dataset and a subset of features. For each node in the dataset, find the optimal feature for splitting to maximize the information gain or Gini impurity. Recursively divide the dataset into child nodes until a stopping condition is reached (such as the number of samples in the leaf node is less than a certain threshold, or all features have been used). (2) Prediction: For a new data instance, input its feature values into each decision tree. For each decision tree, calculate the path by which the instance reaches the leaf node according to the path found in its internal nodes. Each decision tree will output a prediction result (classification or regression), and combine the prediction results of all decision trees. For a classification task, use majority voting to determine the final prediction label. For a regression task, take the average of all decision tree prediction results as the final predicted value. The calculation formula of the random forest can be expressed as:

[0077] [\text{RandomForestPrediction}=\text{Mode}(\text{TreePrediction}_1,\text{Tree Prediction}_2,\ldots,\text{TreePrediction}_T)];

[0078] Among them, (\text{TreePrediction}_i) is the prediction result of the i-th decision tree, and (T) is the number of decision trees. The random forest constructs multiple decision trees, makes predictions for each tree, and then uses majority voting to obtain the final prediction result. This method can improve the accuracy and robustness of the prediction because if a certain decision tree has a bias in a certain feature, other decision trees may have different prediction results, thus correcting this bias through majority voting.

[0079] K-Nearest Neighbors (KNN) is an instance-based learning method that predicts the label of a new instance by finding the K nearest instances to the new instance in the training dataset. The calculation formula of KNN is not complicated and mainly consists of two parts: feature distance calculation and label prediction.

[0080] (3). Feature distance calculation: For each feature of the new instance, calculate the square of the difference between it and the corresponding features of all instances in the training set. Sum up these squared differences to obtain the total distance between each training instance and the new instance. Use Euclidean distance or other distance metrics (such as Manhattan distance) to calculate the distance. (2) Label prediction: Sort all instances in the training set according to the distance, and select the nearest K instances. Calculate the mode of the labels of these K instances as the predicted label for the new instance. The calculation formula of KNN can be expressed as: [\hat{y}=\text{mode}(y_1,y_2,\ldots,y_K)] where (\hat{y}) is the predicted label, and (y_1,y_2,\ldots,y_K) are the labels of the K nearest neighbor instances. KNN finds the K nearest neighbors by calculating the distance between the new instance and each instance in the training set, and then makes a prediction based on the labels of these neighbor instances. Here, the "mode" refers to the label that appears most frequently. The performance of KNN depends to a large extent on the distance metric and the value of K. The distance metric determines how to measure the similarity between two instances, and the value of K determines the number of neighbors considered. Selecting appropriate distance metrics and K values is crucial for the performance of the KNN model.

[0081] (4). Support Vector Machine (SVM) is a binary classification model, and the mathematical expression of its optimization problem is:

[0082] [\min_{\beta,\gamma}\frac{1}{2}|\beta|^2\text{subjectto}y_i(\beta^Tx_i+\gamma)\geq1,\foralli];

[0083] Among them, the meanings of each symbol are as follows:

[0084] · (\beta) is the normal vector of the hyperplane, that is, (\beta) in the hyperplane equation (\beta^Tx+\gamma = 0).

[0085] · (\gamma) is the intercept of the hyperplane, that is, (\gamma) in the hyperplane equation.

[0086] · (y_i) is the label of the i-th sample in the training dataset. For binary classification problems, it is a binary variable, usually represented as (+1) or (-1).

[0087] · (x_i) is the feature vector of the i-th sample in the training dataset.

[0088] The meaning of this formula is that SVM attempts to find a hyperplane (\(\beta^Tx+\gamma=0\)) such that the distance from all sample points in the training dataset to this hyperplane is at least 1. To achieve this, SVM controls the complexity of the hyperplane by minimizing the L2 norm of the normal vector (\(\beta\)) (i.e., \((\frac{1}{2}|\beta|^2)\)), while ensuring that the distance from all training sample points to the hyperplane is at least 1. To solve this optimization problem, SVM usually adopts the Lagrange multiplier method, introducing Lagrange multipliers (\(\alpha_i\)) to handle the inequality constraints and transforming the original problem into a Lagrangian problem. Then, by solving the Lagrangian problem, the optimal values of \(\beta\) and \(\gamma\) are obtained.

[0089] (V). Gradient Boosting Decision Tree (GBDT) is a powerful machine learning algorithm that improves prediction performance by constructing multiple decision trees and combining them. The core idea of GBDT is to iteratively train decision trees, and each training is performed based on the prediction results of the previous model to minimize the loss function. The calculation formula of GBDT can be expressed as: \([f_t(x)=f_{t - 1}(x)+h_t(x)]\);

[0090] Among them, \((f_t(x))\) is the prediction result of the \(t\) - th model (i.e., the combination of the first \(t\) decision trees), \((f_{t - 1}(x))\) is the prediction result of the \((t - 1)\) - th model, and \((h_t(x))\) is the \(t\) - th decision tree. The training process of GBDT can be divided into the following steps:

[0091] Initializing the model:

[0092] Initialize the model as a constant prediction (\(f_0(x)=0\)).

[0093] Iteratively training decision trees:

[0094] For the \(t\) - th iteration, calculate the gradient \(\nabla L(\hat{y},y)\) of the loss function \(L(\hat{y},y)\) with respect to \(f_{t - 1}(x)\).

[0095] Select the negative gradient of a loss function as the new feature (\(g_t(x)\)), i.e., \(g_t(x)=-\nabla L(\hat{y},y)\).

[0096] Use \(g_t(x)\) as the training data to train a new decision tree (\(h_t(x)\)).

[0097] Add \((h_t(x))\) to the current model \((f_{t - 1}(x))\) to obtain a new model \((f_t(x))\).

[0098] Final model:

[0099] Repeat the above steps until the stopping condition is reached (such as the number of iterations reaches the preset maximum number, or the improvement of the loss function is less than a certain threshold).

[0100] (VI). Light Gradient Boosting Machine (LightGBM) is an effective gradient boosting framework that adopts some optimization techniques to improve the training speed and prediction efficiency of gradient boosting decision trees.

[0101] The loss function used by LightGBM is a tree-structured logarithmic loss function, and it accelerates the training process by minimizing the second derivative (i.e., the trace of the Hessian matrix). The calculation formula of LightGBM can be expressed as:

[0102] \[\min_{\beta}\sum_{i=1}^{n}-\log(P(y_i|X_i))+\frac{\lambda}{2}|\beta|^2\];

[0103] Among them, the meanings of each symbol are as follows:

[0104] -\(\beta\) is the model parameter, corresponding to the weight of each leaf node in the decision tree.

[0105] -\(X_i\) is the feature vector of the \(i\)-th sample.

[0106] -\(y_i\) is the target value of the \(i\)-th sample.

[0107] -\(P(y_i|X_i)\) is the conditional probability of the target value \(y_i\) given the feature vector \(X_i\).

[0108] -\(\lambda\) is the regularization parameter, used to control the complexity of the model.

[0109] (VII). Extreme Gradient Boosting / XGBoost (eXtreme Gradient Boosting) is a popular machine learning algorithm that improves the prediction performance by constructing multiple decision trees and combining them. The core idea of XGBoost is to iteratively train decision trees, and each training is carried out on the prediction results of the previous model to minimize the loss function. To control the complexity of the model, XGBoost introduces a regularization term. The calculation formula of XGBoost can be expressed as:

[0110] \[\min_{\beta}\sum_{i = 1}^{n}L(y_i,\hat{y}_i^{(t - 1)}+f_t(x_i))+\Omega(\beta)\];

[0111] Among them, the meanings of each symbol are as follows:

[0112] - \(\beta\) is the model parameter, corresponding to the weight of each leaf node in the decision tree.

[0113] - \(\hat{y}_i^{(t - 1)}\) is the prediction result of the first \(t - 1\) decision trees for the \(i\)-th sample.

[0114] - \(f_t(x_i)\) is the prediction result of the \(t\)-th decision tree for the \(i\)-th sample.

[0115] - \(L(y_i,\hat{y}_i)\) is the loss function, used to measure the difference between the predicted value and the true value.

[0116] - \(\Omega(\beta)\) is the regularization term, used to control the complexity of the model and prevent overfitting.

[0117] - \(n\) is the number of training samples.

[0118] - \(t\) is the number of iterations, indicating the \(t\)-th decision tree trained.

[0119] XGBoost uses the squared error loss function and the regularization term, so the formula can be further expanded as:

[0120] \[\min_{\beta}\sum_{i = 1}^{n}(y_i-\hat{y}_i^{(t - 1)}-f_t(x_i))^2+\lambda\sum_{j = 1}^{m}\beta_j^2\];

[0121] Among them, \(\lambda\) is the regularization parameter, used to control the complexity of the model.

[0122] By minimizing this loss function, XGBoost can learn the complex features of the data and gradually optimize the prediction performance of the model. This feature of XGBoost makes it an effective tool for dealing with regression and classification problems.

[0123] (8). LASSO (Least Absolute Shrinkage and Selection Operator) is a regularization technique for linear regression models. It realizes feature selection and parameter shrinkage by adding an absolute value regularization term to the process of minimizing the loss function. LASSO can make some coefficients become zero, thus automatically performing feature selection and only retaining important features. The goal of LASSO is to minimize the sum of the loss function and the regularization term, and its calculation formula can be expressed as:

[0124] \[\min_{\beta}\frac{1}{n}\sum_{i=1}^{n}\ell(y_i,X_i^T\beta)+\lambda\sum_{j=1}^{p}|\beta_j|\]; where the meanings of each symbol are as follows:

[0125] - \(\beta\) is the model parameter, representing the coefficient vector in the linear regression model.

[0126] - \(\ell(y_i,X_i^T\beta)\) is the loss function, usually the squared loss function, which is used to measure the difference between the predicted value and the true value.

[0127] - \(X_i\) is the feature matrix of the \(i\)th sample.

[0128] - \(y_i\) is the target value of the \(i\)th sample.

[0129] - \(\lambda\) is the regularization parameter, which is used to control the size of the absolute value regularization term.

[0130] - \(p\) is the number of features.

[0131] In LASSO, the regularization term is the absolute value function, which makes the coefficient \(\beta_j\) either remain unchanged or become zero. This characteristic enables LASSO to impose sparsity on the coefficients and only retain the features that make a significant contribution to the prediction.

[0132] The essence of the combined model is a multivariate Logistics Regression model. The basic form of the logistic regression model is:

[0133] \[P(Y=1|X)=\frac{1}{1+e^{-\beta_0-\beta_1X_1-\beta_2X_2-...-\beta_pX_p}}\];

[0134] Where:

[0135] - \(P(Y = 1|X)\) is the probability that event Y occurs given the independent variable X.

[0136] - \(X_1, X_2, \ldots, X_p\) are independent variables.

[0137] - \(\beta_0, \beta_1, \ldots, \beta_p\) are model parameters, also known as coefficients or slopes.

[0138] - \(e\) is the base of the natural logarithm, approximately equal to 2.71828.

[0139] The log-likelihood function of the logistic regression model is:

[0140] \[\log L(\beta_0, \beta_1, \ldots, \beta_p)=\sum_{i = 1}^{n}[y_i\cdot(\beta_0+\beta_1X_{1i}+\beta_2X_{2i}+\ldots+\beta_pX_{pi})-\log(1 + e^{\beta_0+\beta_1X_{1i}+\beta_2X_{2i}+\ldots+\beta_pX_{pi}})]\];

[0141] Where:

[0142] - \(y_i\) is the value of the target variable (0 or 1).

[0143] - \(n\) is the number of samples.

[0144] By maximizing the log-likelihood function, we can estimate the parameters of the logistic regression model. This is usually achieved by using methods such as Maximum Likelihood Estimation (MLE) or, in some cases, regularization methods (such as Lasso for L1 regularization or Ridge regression for L2 regularization).

[0145] V. Specifically mentioned, the following are the calculation formulas for these metrics:

[0146] (I). AUC (Area Under the ROC Curve): AUC is the area under the Receiver Operating Characteristic (ROC) curve. It is a metric for measuring the ability of a model to distinguish between positive and negative samples, with a value range from 0 to 1. The closer the value is to 1, the better the performance of the model. The calculation formula for AUC is:

[0147] \[AUC = \frac{1}{n}\sum_{i = 1}^{n}\max(0, 1 - |y_i - y_j|)\]; where \(y_i\) and \(y_j\) are the true values of two randomly selected samples.

[0148] (2) Accuracy: Accuracy is the proportion of correctly predicted samples in the total samples. Its calculation formula is: \[Accuracy = \frac{TP + TN}{TP + TN + FP + FN}\] where TP is true positive, TN is true negative, FP is false positive, and FN is false negative.

[0149] (3) C-index: The C-index is the square root of the area under the receiver operating characteristic curve (ROC curve). It is an index to evaluate the prediction stability of the model, and the value range is from 0 to 1. The closer it is to 1, the better the performance of the model. The calculation formula of the C-index is: \[C = \sqrt{AUC}\] Therefore, the calculation of the C-index is actually to calculate the square root of AUC.

[0150] (4) F1-score: The F1-score is the harmonic mean of precision and recall. Its calculation formula is: \[F1 = 2\times\frac{Precision\times Recall}{Precision + Recall}\] where Precision is precision and Recall is recall.

[0151] (5) Sensitivity (Recall): Sensitivity or recall is the proportion of true positives in all positive samples. Its calculation formula is:

[0152] \[Sensitivity = Recall = \frac{TP}{TP + FN}\]; where TP is true positive and FN is false negative.

[0153] (6) Specificity: Specificity is the proportion of true negatives in all negative samples. Its calculation formula is:

[0154] \[Specificity = \frac{TN}{TN + FP}\]; where TN is true negative and FP is false positive.

[0155] (7) Precision: Precision is the proportion of true positives in all samples predicted as positive. Its calculation formula is:

[0156] [Precision=\frac{TP}{TP + FP}]; where TP is true positive and FP is false positive.

[0157] VI. Calculate the confidence interval using statistics and resampling / Bootstrap method:

[0158] Specifically, to calculate the confidence interval using statistics and the Bootstrap method, the method is as follows: (1) Z - score method: If the sample size is large enough (usually greater than 30), the distribution of the sample mean is approximately normal. In this case, the Z - score method can be used to calculate the confidence interval. The formula is: [\bar{x}\pm Z\times\frac{\sigma}{\sqrt{n}}];

[0159] where (\bar{x}) is the sample mean, (Z) is the Z - score (a constant related to the confidence level), (\sigma) is the population standard deviation, and (n) is the sample size. (2) The calculation steps of the Bootstrap confidence interval usually include: (1) Randomly sample with replacement from the original sample multiple times, each time drawing a sample of the same size as the original sample. (2) Calculate the statistic (such as sample mean, sample proportion, etc.) for each sampling. (3) Record the values of the statistic obtained from each sampling. (4) Calculate the confidence interval based on these values. Usually, the upper and lower limits of the confidence interval are the 2.5% and 97.5% quantiles of these values.

[0160] VII. Delong Test / DeLong test. The Delong test is used to compare whether there is a significant difference in the AUC values of two models. Its calculation formula involves the calculation of the Mann - Whitney statistic, covariance, and correlation coefficient. The specific formula is relatively complex and involves multiple steps. The following is an overview of its main calculation steps:

[0161] (I). Mann - Whitney statistic:

[0162] [\text{Mann - Whitney statistic}=\frac{1}{n}\sum_{i = 1}^{n}\text{Kernel}(X_i,Y_i)];

[0163] where (\text{Kernel}(X_i,Y_i)) is an indicator function that equals 0.5 if (Y_i\leq X_i) and 1 otherwise.

[0164] Calculation of the covariance matrix S:

[0165] [S_{11}=\frac{1}{n}\sum_{i = 1}^{n}(X_i - \text{AUC}A)(Y_i - \text{AUC}A)]

[0166] [S{12}=\frac{1}{n}\sum{i = 1}^{n}(X_i - \text{AUC}A)(Y_i - \text{AUC}B)]

[0167] [S{21}=S{12}]

[0168] [S_{22}=\frac{1}{n}\sum_{i = 1}^{n}(X_i - \text{AUC}_B)(Y_i - \text{AUC}_B)]

[0169] Where (\text{AUC}_A) and (\text{AUC}_B) are the AUC values of the two models respectively.

[0170] Calculate the z - score and p - value:

[0171] [z=\frac{\text{AUC}_A - \text{AUC}_B}{\sqrt{S}}]

[0172] [p = 2\times\text{CDF}(z)]

[0173] Where (S) is the covariance matrix and (\text{CDF}) is the cumulative distribution function;

[0174] Find_Optimal_Cutoff;

[0175] This function is used to find the optimal decision threshold to maximize the Youden index. The Youden index is the harmonic mean of sensitivity (true positive rate) and specificity (true negative rate). Its calculation formula is as follows:

[0176] [\text{Youdenindex}=\text{sensitivity}-\text{specificity}];

[0177] In the function, the optimal threshold is determined by finding the threshold that maximizes the Youden index.

[0178] VIII. 5-fold cross-validation: It is a technique for evaluating the performance of machine learning models. The dataset is divided into 5 equal parts. Each time, 4 of these parts are used as the training set and 1 part as the test set, and this is repeated 5 times, with a different test set each time. This allows the performance of the model to be evaluated on different subsets of training and test data, resulting in a more robust model evaluation. In 5-fold cross-validation, there is no single calculation formula, but rather it involves multiple steps of model training and evaluation. Specifically, this process can be divided into the following steps:

[0179] Divide the dataset into 5 equal parts;

[0180] For each part, repeat the following steps: Select 4 parts as the training set. Select the remaining 1 part as the test set. Train the model on the training set. Evaluate the performance of the model on the test set;

[0181] Record the results of each evaluation; Calculate the average of all the results as an estimate of the model's performance;

[0182] IX. Shapley Additive Explanations / SHAP explains predictions by calculating the contribution of each feature in the cooperation. The specific calculation process is as follows: Calculate the Shapley value of the feature: For each feature, calculate its Shapley value in all feature combinations. The calculation formula for the Shapley value is as follows:

[0183] [\text{SHAP}(f_i)=\frac{1}{n!}\sum_{\substack{S\subseteq\text{dom}(f)\S\neq\emptyset}}\phi(S)\cdot\text{Imp}(f_i|S)];

[0184] Where:

[0185] (f_i) is the (i)-th feature.

[0186] (\text{dom}(f)) is the set of all features.

[0187] (n) is the total number of features.

[0188] (\phi(S)) is the Shapley value of the feature combination (S).

[0189] (\text{Imp}(f_i|S)) is the measure of the importance of feature (f_i) in set (S).

[0190] Among the eight single models constructed, the radiomics model constructed using LASSO has the best performance.

[0191] Calculate the radiomics score for each patient based on the coefficients of the LASSO model. The score is usually a linear combination of the model coefficients and feature values. Therefore, the calculated radiomics score can be expressed as:

[0192] y = 0.453809521 - 0.0000000328 * original_glszm_LargeAreaHighGrayLevelEmphasis_1 - 0.00216 * original_glszm_SmallAreaHighGrayLevelEmphasis_1 - 0.0000222 * original_gldm_LargeDependenceHighGrayLevelEmphasis_1 + 0.00186 * wavelet-HHL_firstorder_Maximum_1 - 0.00384 * wavelet-HHH_firstorder_Maximum_1 - 0.000000000197 * wavelet-HHH_glszm_LargeAreaHighGrayLevelEmphasis_1 + 0.000000000152 * wavelet-HHL_firstorder_TotalEnergy_2 - 0.00127 * squareroot_firstorder_Kurtosis_2 + 0.000784 * squareroot_glszm_GrayLevelNonUniformity_2 - 0.00000122 * logarithm_glrlm_RunLengthNonUniformity_2 - 0.000000147 * logarithm_glszm_LargeAreaLowGrayLevelEmphasis_2. The purpose of this formula is to calculate the radiomics score for each patient, and this score can be used to evaluate the risk of a patient having cancer. The higher the score, the higher the risk of the patient having cancer. This score is based on the prediction of the model, which takes into account the coefficients and feature values of all features.

[0193] Among the eight models constructed for single model building, the dosomics model constructed by Random Forest (RF) has the best performance. However, RF is a non-linear model that improves prediction accuracy by constructing multiple decision trees and combining their predictions. Due to the complexity and non-linear characteristics of the RF model, there is usually no direct formula to represent the dosomics risk score. Although the random forest model itself does not have a direct formula to represent the dosomics risk score, a score representing the patient's risk level can be obtained by calculating the prediction probability of the model on the test set. This score can be used as the dosomics risk score to predict the class probability: Use the trained model to predict new patient data to obtain the probability that each sample belongs to the positive class. This probability can be calculated using the following formula: P(positive class) = (number of samples predicted as positive class by the model / total number of samples) Prediction probability: In some cases, you may want to obtain a more refined probability estimate rather than a simple binary classification probability. This can be done by calculating the probability that each sample is predicted as the positive class in each decision tree and then averaging or aggregating these probabilities. This usually involves calculating the path probability of each sample in each decision tree and then taking a weighted average of these probabilities. P(positive class) = (weighted average of the probability that the sample is predicted as the positive class in the decision tree) The weighted average here can be a simple arithmetic mean or a weighted average considering the importance of each decision tree in the model. These calculations are usually implemented by writing code (such as using Python, R, or MATLAB, etc.) rather than through a simple formula. In practical applications, specialized machine learning libraries (such as scikit-learn, the randomForest package in R, etc.) are needed to calculate these probabilities.

[0194] In the constructed delta-radiomics model, XGBoost performs the best. When dealing with classification problems, the XGBoost model usually outputs the predicted probabilities for each class. These probability values are obtained through the softmax function or other forms of probability transformation, which represent the confidence of the model in predicting each class. In regression problems, the XGBoost model outputs continuous values, that is, the values of the predicted target variables. In XGBoost, you can specify the objective of the model by setting the objective parameter. For example, for binary classification, you can set objective='binary:logistic'. This will cause the model to output the predicted probabilities for each class. This predicted value is the obtained delta-radiomics risk score, and what the predict function returns is the probability that each sample belongs to the positive class. The predicted class is determined by comparing the predicted probability with a certain threshold (such as 0.5). For example, if predictions[i] is greater than 0.5, the predicted class is the positive class, otherwise it is the negative class. When dealing with regression problems, the XGBoost model directly outputs the predicted continuous values. You can specify the objective function for the regression problem by modifying the objective parameter. For example, objective='reg:squarederror' is used to minimize the squared error between the predicted value and the true value. What the predict function returns is the predicted target value for each sample. These predicted values can be directly used as radiomics risk scores or for further analysis.)

[0195] The essence of the combined models is all logistic regression / LogisticsRegression. The calculation expression of the first combined model (ModelA+B+D) is

[0196] log(p / (1-p)) = -2.9011 + 4.7678*D-score + 6.0579*R-score + 0.9507*SI ≤

[0197] 400 + 1.4669*SI > 400 - 2.3986*NoInducchemo - 1.8985*NoCCRT + 0.936*NoConsochemo where p represents the probability that the sample belongs to the positive class, and the clinical factors and dosimetric parameters involved include SI / smoking index

[0198] / smokingindex; V40 represents the lung volume when receiving 40 Gy;

[0199] NoInducchemo / concurrentchemoradiotherapy represents no induction chemotherapy;

[0200] NoCCRT / concurrent chemoradiotherapy represents no concurrent chemoradiotherapy;

[0201] The expression of the second combined model (Model A+B+C+D) is:

[0202] log(p / (1-p)) = -5.1077 + 4.6731*R-score + 5.5775*D-score + 2.6583*Delta-R-score - 2.2927*AJCCStageII - 2.5096*AJCCStageIII - 2.9755*AJCCStageIV + 0.1847*V20 - 1.6841*NoIndu cchemo, where p represents the probability that the sample belongs to the positive class, and the clinical factors and dosimetric parameters involved include AJCCStageII representing the American Joint Committee on Cancer (AJCC) stage II; AJCCStageIII representing the AJCC stage III; V20 representing the lung volume at 20 Gy;

[0203] NoInducchemo / concurrent chemoradiotherapy represents no induction chemotherapy).

[0204] a and c are the ROCs of different models in the training set, and b and d are the ROCs of different models in the test set. After constructing 8 models respectively, the best models are selected from each omics. They are the radiomics model (Model A) LASSO ROC = 0.718 (95% CI 0.598, 0.839), the dosomics model (Model B): RF ROC = 0.702 (95% CI 0.580, 0.824), and the delta-radiomics model (Model C): Xgboost ROC = 0.709 (95% CI 0.558, 0.860). The ROC of the combined model A+B+D is 0.840 (95% CI 0.740, 0.940), and the ROC of the three-omics combined model A+B+C+D is 0.839 (95% CI 0.721, 0.926). The effect of the combined model is better than that of the individual omics models.

[0205] It should be noted that, in this text, relational terms such as first and second are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the term "comprising", "including" or any other variant thereof is intended to cover non-exclusive inclusion, so that a process, method, article or device comprising a series of elements not only includes those elements, but also includes other elements not expressly listed, or further includes elements inherent to such process, method, article or device.

[0206] Although the embodiments of the present invention have been shown and described, those of ordinary skill in the art can understand that various changes, modifications, substitutions and variations can be made to these embodiments without departing from the principles and spirit of the present invention.

Claims

1. A prediction model for radiation pneumonitis based on computed tomography, characterized in that: It includes the following steps: S1: Image segmentation module, which collects the initial planned CT image (CT1), the reduced field planned CT image (CT2) and the initial planned dose distribution map of the case, and obtains the corresponding region of interest; S2: perform feature extraction; S3: perform feature selection; S4: Conduct model building, construct single models, combined models, and nomograms; S5: Conduct model evaluation; S6: Perform model interpretation.

2. The prediction model for radiation pneumonitis based on computed tomography according to claim 1, characterized in that: In S2, the lung areas that received irradiation doses of 20 Gy and 30 Gy or more on the CT image (CT1) were defined as V20 and V30 as regions of interest, and the radiomic features and dosimetric features within the regions of interest were extracted. Patients who underwent 40 Gy-50 Gy radiotherapy and underwent reduced-field planning CT images (CT2) were selected from the patients, and the radiomic features of the regions of interest V20 and V30 in the reduced-field planning CT images (CT2) were extracted. Delta-RF=RFCT2-RFCT1 was used to obtain delta-radiomics features.

3. The prediction model for radiation pneumonitis based on computed tomography according to claim 1, characterized in that: In S3, features with non-zero coefficients are screened using variance selection method, Pearson correlation coefficient, least absolute shrinkage and selection operator.

4. The prediction model for radiation pneumonitis based on computed tomography according to claim 1, characterized in that: In S5, the area under the ROC curve is used to evaluate the performance of the model, the confidence interval is calculated using statistics and the bootstrap method, and the DeLong test is used to test the significant differences in the ROC curve areas of different models.

Images