Radiomics Based Bayesian Inversion Method for Prediction of Cancer and Pathological Stage

Objective: To develop a Bayesian inversion framework on longitudinal chest CT scans which can perform efficient multi-class classification of lung cancer. Methods: While the unavailability of large number of training medical images impedes the performance of lung cancer classifiers, the purpose built deep networks have not performed well in multi-class classification. The presented framework employs particle filtering approach to address the non-linear behaviour of radiomic features towards benign and cancerous (stages I, II, III, IV) nodules and performs efficient multi-class classification (benign, early stage cancer, advanced stage cancer) in terms of posterior probability function. A joint likelihood function incorporating diagnostic radiomic features is formulated which can compute likelihood of cancer and its pathological stage. The proposed research study also investigates and validates diagnostic features to discriminate accurately between early stage (I, II) and advanced stage (III, IV) cancer. Results: The proposed stochastic framework achieved 86% accuracy on the benchmark database which is better than the other prominent cancer detection methods. Conclusion: The presented classification framework can aid radiologists in accurate interpretation of lung CT images at an early stage and can lead to timely medical treatment of cancer patients.


I. INTRODUCTION
Image analysis techniques have successfully provided personalized prognosis and treatment plans for cancer patients with a greater accuracy in the recent years. There are several studies available that show the importance of imaging features for treatment, monitoring and outcome prediction in lung cancer as well as other cancer types [1], [2]. Besides prognosis, the quantitative imaging features also known as radiomic features have been extracted from the medical images to be further used in clinical research for lung cancer diagnosis [3], [4]. Most of the prominent radiomics based classification methods for lung cancer have incorporated neural networks and deep learning approaches [5]- [13]. A multiview convolutional network is proposed in [9] for binary (benign and malignant) and ternary classification (binary, primary malignant and metastatic malignant). The proposed model achieved error rate of 5.41% for binary classification and 13.91% for ternary classification respectively. Another multi-view knowledge-based collaborative (MV-KBC) deep model classified nodules on limited chest CT data in [11]. The prediction error computed in terms of Root Mean Square Error (RMSE) of the model is 0.62. Authors in [12] investigated the usefulness of a deep convolutional neural network (DCNN) for benign, primary lung cancer and metastatic cancer classification. The best averaged accuracy of 68% was achieved with transfer learning for image size of 224.
Another common approach adopted to address classification problem is the application of machine learning classifiers. A structural co-occurrence matrix (SCM) based method was proposed by [13]. The proposed technique did not only classify the nodules as benign or malignant but also classified the malignant nodules into 5 prescribed levels with an accuracy of 74.5%. Yu et al. in [5] presented a study to detect pathological stages I, II, III and IV respectively in non-small cell lung cancer (NSCLC) using Random Forest algorithm. Out of the bag (OOB) error was computed to measure classification accuracy of the Random Forest algorithm whereas the predicted cancer stages of test data sets achieved ROCs between 0.6 to 1.
However, there is a dearth of scientific approaches to diagnose the benign or cancerous nature of a lung nodule along with its pathological stage in case of malignancy. Furthermore, the insufficient size of training data sets required for neural networks training is a challenge for accurate cancer stage detection of the tumor. Keeping in view the aforementioned challenges, there is an emerging need to investigate new sophisticated solutions for the detection of malignancy status in lung nodules. A useful approach towards joint estimation of tumor malignancy and its stage can involve tracking of tumor characteristics using temporal analysis of CT images [14], [15]. Such an approach also offers good potential to address the problem of smaller data sets in medical imaging.
In this research article, a Bayesian inversion method is proposed which employs longitudinal radiomics data extracted from CT scans of lung nodules to detect cancer and its pathological stage. The presented system adopts particle filtering approach to compute the posterior Probability Density Function (PDF) of cancer occurrence and cancer pathological stage respectively. For PDF estimation, a joint likelihood function incorporating diagnostic radiomic features is developed which can diagnose malignant nodules and the pathological cancer stage. The cancer stage is determined as early stage cancer (stage I & stage II) or advanced stage cancer (stage III & stage IV) respectively. Furthermore, radiomic features which can discriminate between cancer stages are also investigated and evaluated.
The presented Bayesian inversion method was validated on a total of 200 nodules with 130 malignant and 70 benign nodules. There were 109 malignant and 63 benign nodules which were successfully diagnosed showing better performance of the proposed framework than the other prominent methods in literature.
In the following sections, research problem is formulated and proposed stochastic frame work is described. Next, implementation of the presented model is discussed. Then, results and discussion are presented followed by a conclusion. The notations used in this research work are listed in Table 1 for reader's understanding.

II. PROBLEM FORMULATION
A number of radiomic features depicting shape and gray level intensities of lung nodule have shown diagnostic power towards lung, colon and head and neck cancer [16]. The aim of the presented research is to utilize such radiomics features measurements from periodic chest screening images in a stochastic framework to detect lung cancer. The sequential evaluation of nodules' malignancy status given the radiomics data from periodic screening images is an inverse problem and can be solved using Bayesian inversion approach [17]. There are several methods available to solve the aforesaid problem but a recursive Bayesian filter based on sequential Monte Carlo method is proposed since it uses nonlinear models [18]. The need of nonlinear modeling arises due to the nonlinear behavior of measured radiomics features and the associated likelihood function for lung cancer diagnosis [16].
The proposed Bayesian inversion framework requires a likelihood function for PDF estimation. Since multi-class classification using a single mathematical model has not been a preferred approach due to the complex relationship between benign and cancerous (multiple stages) nodules [19], [20], an intuitive approach has been adopted to combine two mathematical likelihood equations into one for PDF estimation. The proposed system for cancer prediction is shown in Fig. 1. In practice, an inverse problem can be represented by two models namely state transition model and measurement model [21] and is shown in Fig. 2 with states represented by u and measurements by z. A description of these two models and their application in Bayesian inversion method for malignancy estimation is given in the following section.

III. DIAGNOSTIC STOCHASTIC FRAMEWORK A. STOCHASTIC STATE-SPACE MODELS
To start with, a state vector U = {u 1 ,u 2 ,u 3 ,u 4 , . . . . ., u k−1 , u k } is defined where each element u k is a diagnostic score between 0 and 1 that represents the benign/malignancy state of a lung tumor at a screening interval k. A state transition model represents the evolution of state u k with respect to time instants and can be expressed as [22]: (1) VOLUME 9, 2021 Here x k (.) is a non-linear and time varying function that models the state transition process and v k is process noise respectively.
The q number of diagnostic radiomic feature measurements of a lung nodule in a CT image at each screening interval k are assumed to be z 1k , z 2k , . . . , z qk and are represented by a measurement vector Z such that Z = {z 11 , z 21 , . . . ., z q1 , z 21 , z 22 , . . . z q2 ,. . . , z 1k , z 2k , . . . , z qk }. The measurement function relates the features measurements to the true malignancy states of a lung tumor and can be expressed as: Here y k (.) is a time varying and independent function that represents measurement model and n k is the measurement noise. Since, there has been a lack of discussion on the linearity of radiomics based state transition and measurement functions in the literature, a more generalized approach is adopted to perform the Bayesian inversion estimation which is offered by particle filter to estimate the posterior probability of malignancy and pathological cancer stage respectively in lung nodules.

B. PARTICLE FILTERING
A particle filter is an algorithm which represents the desired posterior PDF of an unknown state as weighted samples, also known as particles. Using a large number of samples (particles), the estimated PDF is approximated to the true PDF. In particle filtering approach to dynamic state estimation of a lung nodule, the posterior PDF of the unknown malignancy state u k is constructed based on the set of received feature measurements Z k .
Given the likelihood p(Z k |u k ) of measured diagnostic features Z k at time k, the posterior density p(u k |Z k ) of unknown state u k for the set of received features Z k can be computed using Bayes' Theorem [23] as follows: Here the normalizing constant p(z k | Z k−1 ) is given by: The normalization constant depends upon the likelihood function p(Z k |u k ), defined by the measurement model in Eq. (2). Recursively, the filter predicts and updates to estimate the sequential unknown malignancy states of nodules. In summary, the proposed particle filter expresses the PDF p(u k |Z k ) of a nodule's unknown state at screening interval k as a set of weighted samples. The total number of samples drawn from an importance density are assumed N s . We represent a sample point as u i k and its associated weight as w i k where i = {1, 2, . . . ., N s }. A randomly measured weighted sample at any screening interval can be represented as a set {u i k , w i k }. Using the principle of importance sampling, the weights are normalized such that w i k = 1. Then, the posterior density at a screening time k can be approximated as: The samples can be drawn from the importance density q(Z i k |U k ) and the weights can be assigned using the following association: If the samples at screening time k − 1 are available constituting to probability p(z 1:k−1 | u 1:k−1 ), then p(z 1:k | u 1:k ) is approximated with a new set of samples once the observation z k is known at the nodule screening time k.
The prior or state transition probability p(u k |u k−1 ) of malignancy in a nodule is modeled using exponential distribution because the progression rate of lung cancer follows exponential distribution [24]. Since the transition density depends upon the malignancy states of cancer u k−1 and u k at any consecutive times k − 1 and k respectively, it can be represented as follows: where α is a scalar value that controls the variability in the predicted state value.

IV. IMPLEMENTATION OF PARTICLE FILTER
The test datasets which are given as input to the proposed model comprises of longitudinal radiomics data of patients' CT images acquired from National Lung Screening Trial (NLST) [25].

A. NLST DATASETS
NLST is a randomized multi-site trial which enrolled 53,454 people with high cancer risk to assess the impact of using low dose CT screening in mortality reduction from lung cancer. In NLST project, low dose CT screenings of the enrolled people were performed for three consecutive years which are termed as T0, T1 and T2 respectively. For the implementation of proposed filter, longitudinal data of participants who had cancer detected after the third year (T2) of screening were included. This ensured availability of CT data for all the three years in order to apply the prediction model. A total of 130 cases with positive screening at either 4300208 VOLUME 9, 2021 year T0 or T1 followed by screen detected cancer at year T2 were acquired. Another 70 CT data sets with benign nodules detected at year T2 and with positive screening at either year T0 or T1 were accessed from the database.
The detected nodules from all the 3 consecutive CT screenings of each patient were segmented and 105 3-D features were computed. This resulted in a set of 3 longitudinal radiomic readings for each patient. In order to tackle the constraint of fewer data points, single data point was interpolated using cubic-spline method [17] between data obtained from year T0 and year T1, and between year T1 and year T2 respectively. The total number of screening intervals after interpolation were 5(K = 5) within a duration of 3 years. Starting with an initial CT scan, there are 4 sample points following and each point represents a half-yearly CT screening data of the patient which was fed into the proposed system for PDF estimation.

B. LIKELIHOOD FUNCTION FOR CANCER AND CANCER STAGE DETECTION
For PDF estimation using particle filter, a likelihood function has to be developed. In the literature reviewed, authors in [16] formulated a likelihood function for cancer detection (benign or malign) incorporating two highly diagnostic radiomic features namely Surface volume ratio and Sum entropy values of a nodule and is given as: Here z 1k represents the Surface volume ratio and z 2k denotes the Sum entropy value of nodule at k th quarter. The value of p 1 (Z k |u k ) is below 0.51 for benign state and equal or above 0.51 for malignancy in the nodule. The coefficients of Eq. (8)  As reported in a research study [26], a single multi-class likelihood function is prone to large errors hence cancer detection and cancer stage detection have to be performed using two distinct mathematical functions.
Malignancy with cancer pathological stage I or II was collectively termed as early stage cancer whereas stage III and IV cancer was termed as advanced stage cancer [27]. The training cohort comprised of 73 nodules with early stage cancer and 127 nodules with advanced stage cancer from Lung1 [28] database. A pre-processing step of features extraction and features reliability test [16] was applied to the training data sets which resulted in 51 stable and reliable features reported in Table SI (Supplementary material). Furthermore, the redundant features were removed and highly diagnostic features were selected using Least absolute shrinkage and selection operator (LASSO) [10], a method of regression analysis. LASSO technique performs L1 regularization of the features via penalized estimation function to achieve a reduced set of diagnostic features which were validated using 5-fold cross validation. Fig. 3 shows the LASSO coefficients of predictors on y-axis and the L1 norm (the penalization parameter also known as lambda) of coefficients on x-axis for a regression fit. The top x-axis of the plot represents degree of freedom denoted by df and it shows the non-zero coefficients of predictors. Following 7 diagnostic features were obtained after the application of LASSO technique (lambda = 0.0282 with minimum MSE) on pre-processed features: 1

) Sphericity 2) Dependence Variance (DV) 3) Large Dependence High Gray Level Emphasis (LDHGLE) 4) Cluster Prominence (CP) 5) Large Area Low Gray Level Emphasis (LALGLE) 6) Small Area Emphasis (SAE) 7) Strength
The LASSO coefficients of the aforementioned diagnostic features are reported in Table SII. Their diagnostic power ranking was computed using mean scores of 3 supervised feature ranking approaches namely feature based Neighborhood Component Analysis (fNCA), ReliefF network and Infinite Latent Feature Selection (ILFS) respectively [29]. Out of these 7 coefficients, five most discriminating features were used to formulate the likelihood function for early stage or advanced stage cancer detection as follows: +3.0780e −08 (z 5k ) + 0.3917(z 6k ) − 0.00361(z 7k ) (9) Here p 2 (Z k |u k ) indicates likelihood of cancer stage as early stage or advanced stage cancer. It produces a value between 0.51 and below 0.70 for early stage cancer and assumes a value between 0.70 and 0.90 for the detection of advanced stage cancer. The features z 3k , z 4k , z 5k , z 6k and z 7k represent Sphericity, LDHGLE, Cluster Prominence, Small Area Emphasis and Strength respectively for k th CT screening. The coefficients z 4k , z 5k of features LDHGLE and Cluster Prominence respectively are very small and can be excluded. The proposed likelihood function was validated on 47 nodules with early stage cancer and 143 nodules with advanced stage cancer. The function detected 40 nodules with early stage cancer and 125 nodules with advance stage cancer following the designated score ranges. Overall, the proposed function achieved an accuracy of 86.84% (CI:81.19% to 91.30%), specificity of 85.11% (CI:71.69% to 93.80%) and sensitivity of 87.41% (CI:80.84% to 92.37%) respectively with an AUCROC of 84.61%. The ROC plot of the nodule pathological stage classification using the presented likelihood function is shown in Fig. S1. The correlation of these top 5 diagnostic features (LASSO method) with two cancer groups (early stage and advanced stage) is further demonstrated using hierarchical clustering in a heatmap in Fig. 4. The numerical values of plotted radiomic features are standardized with mean 0 and standard deviation 1. The features are represented on x-axis and their hierarchical clustering is represented on y-axis. Two distinct clusters were formed showing strong association of CP and LDHGLE with the defined cancer groups. However, the heatmap could not fully capture the association of other features with two cancer stages besides CP and LDHGLE.
Finally, in order to integrate cancer and its stage detection equations into one likelihood function, highly diagnostic features namely Surface volume ratio (SVR) and Sum entropy (SE) were chosen as the selection criterion [16]. It was observed through an empirical study carried out on the training data that values of SVR (z 1k ) and SE (z 2k ) are always above threshold values T a and T b respectively for malignant nodules and therefore can be used to estimate likelihood of early and advanced stage cancer from Eq. (9). On the contrary, the values of z 1k and z 2k below the threshold T a and T b respectively indicate that nodule is benign/non-malignant and its likelihood can be computed using Eq. (8) Using the above mentioned criterion and Eq. (8) and Eq. (9), the proposed likelihood function for multi-class classification is given as follows: The sequence of steps carried out for the proposed diagnostic framework implementation is demonstrated in Fig. 5 and is summarized as follows: 1) To start, N s samples are generated from the prior/state transition PDF in Eq. (7) for the first time screening of a lung nodule. The value of N s is equal to 4000. 2) Next, the features are computed using the measurement function for the known malignancy state from Step 1. For this, an M th order polynomial is proposed as a measurement function as follows: Here z qk represents q th diagnostic radiomic feature at k th chest screening. The coefficients p m of polynomial were determined from a training database of known states and the corresponding diagnostic feature measurements. As discussed, two features (SVR(z 1k ), SE(z 2k )) were required to estimate benign /malignant state of a nodule [16] and their measurement functions were formulated as: 4300208 VOLUME 9, 2021 The coefficients p 1m and p 2m were computed for z 1k and z 2k respectively. Similarly, measurement functions of cancer stage estimation were also developed.

3) Using the features computed in
Step 2 and the test radiomics features, two likelihood PDFs were estimated from Eq. (10). The closeness of the two likelihood values obtained determines the sample weight.
Large difference refers to small weight assignment and vice versa. Subsequently, the weights were updated using association in Eq.(6). 4) The problem of degeneracy with particles is addressed by performing resampling [30]. After resampling, the posterior PDF of malignancy occurrence was estimated using samples and their associated weights. 5) The PDF estimation process outlined above was repeated for all the cancer screening times from 1 up till K . 6) The expectation of state computed from the PDF estimates at k was represented as posterior mean estimate (PME) as follows: 7) If the error in single point PDF estimates in two consecutive iterations was within a preset threshold ξ , each iteration was denoted by N and the estimated malignancy state of a nodule (K screening times) evaluated at N was u N 1:K , then the convergence criterion was given by:

V. DATA SETS SUMMARY
The proposed particle filtering approach was validated using 200 low dose CT data sets collected from NLST database. From these, there are 130 longitudinal datasets suffering from cancer stage I, II, III or IV respectively at the third year (T3) of screening. Other 70 datasets include patients' longitudinal CT data diagnosed with benign nodules at the third year (T3) of screening. The dimensions of the acquired images were 512 x 512 pixels. A summary of the test nodules is reported in Table SIII.

VI. RESULTS
The 3-D segmentation of nodules in lung CT volumes was carried out using Grow Cut algorithm after the verification from a senior radiologist. Pyradomics [4] was used for the extraction of 3-D radiomic feature extraction. The diagnosis results of the filter were evaluated by setting the value of M as 2 and 3 respectively in measurement functions, described in subsection IV-C, Step 2.
With the order of measurement equations set to 2 (M = 2), the particle filter diagnosed 65 out of 70 benign tumors successfully. In malignant cases, the proposed framework diagnosed 51 out of 70 nodules having early stage cancer and 37 out of 60 nodules with advanced stage cancer correctly. For a comparison, the performance of the proposed framework was assessed using measurement equations of order 3 (M = 3). A total of 63 out of 70 benign tumors were successfully detected whereas in malignant cases, there were 60 out of 70 nodules with early stage cancer and 49 out of 60 nodules with advanced cancer stage which were correctly diagnosed. Increasing the order M to 4 did not show any significant improvement in the detection results and are excluded in the discussion.
The described diagnostic performance of the filter is expressed in terms of precision and recall for 95% confidence interval along with confusion matrix in Fig. 6. The recall values determine the sensitivity of classification model and precision values compute the positive predictive value of model. The average precision of the filter decreased from 92.86% to 90.0% for benign nodule classification but increased from 72.857% to 85.714% and 61.667% to 86.667% respectively for early and advanced stage cancer estimation when M value was incremented from 2 to 3. The recall metric improved significantly for both benign as well as early stage classification with M = 3. The overall accuracy of the filter is 86%.
The estimated diagnostic scores of benign as well as malignant test nodules at final CT screening (k = 5) which is the 4 th half yearly reading are plotted against each nodule count for M = 2 & 3 in Fig. 7. Through visual inspection, an overall improvement in the state estimation of malignant lung nodules can be seen using measurement functions of order 3 in PF. To summarize, the diagnostic scores using M = 2 showed a tendency to shift towards the lower diagnostic ranges causing a spill over in the adjacent classification group. For example, 5 more nodules with early stage cancer obtained a score below 0.51 with M = 2, hence were diagnosed as benign nodules. Similar trend was observed in nodules with advanced stage cancer where 9 more nodules got diagnosed having early stage cancer stage due to low scores for M = 2. The above discussed trend in state estimation is illustrated by plotting the state evolution (diagnostic scores) of 1 benign, 1 malignant (cancer stage IB) and 1 malignant (cancer stage IIIA) against the half-yearly CT screenings in Fig. S2. Furthermore, the ROC curves demonstrating the performance of proposed filter for classification of benign, early stage and advanced stage cancer are shown in Fig. S3. The AUCs achieved for benign nodule, early stage cancer and advanced stage cancer detection are 0.8962, 0.8849 and 0.8691 respectively. Computed values of diagnostic features and likelihood scores for stage I B(early stage) and stage III A(advanced stage) cancer are shown in Fig. 8.

VII. DISCUSSION
The binary classification performance of proposed model into benign and malignant nodules is compared with two other significant models tested on NLST database in Table SIV. The research work in [31] investigated the performance of different classifiers to predict malignancy in screen detected cancer at first year (T0) and second year (T1) of Lung screening trial. A total of 208 nodules were detected among which 104 were diagnosed with cancer. An accuracy (ACC) of 80%, specificity (SPS) of 91% and AUC = 0.83 were obtained. Authors in [6] analyzed deep learning algorithm generated heat-maps to predict the malignancy in nodules. The presented explainable AI approach to detect malignancy achieved a weighted accuracy of 85%. In comparison, the proposed PF achieved an accuracy of 92.0%, sensitivity of 87.50% and specificity of 94.5% respectively showing excellent performance of the filter in binary classification.
Since no multi-class classification techniques have been reported on NLST datasets in the literature, multi-class classification was performed using SVM with Radial Basis Function (RBF) kernel, Random Forest and K nearest neighbors KNN) classifiers for a comparative performance analysis. The data sets were chosen from third year (T2) screening with 180 datasets used for training and the remaining 20 datasets were chosen for the testing of classifiers. The number of benign and malignant nodules in the experiments were same as used in the proposed filter and reported in Table SIII. The classifier was trained and tested 10 times with 10-fold cross validation. The evaluation metric included averages (for 10 iterations) of precision, recall, accuracy and Mathew's Correlation Coefficient (MCC) and is reported in Table 2. Classification of nodules as benign, with early stage cancer and advanced stage cancer using SVM, KNN and Random Forest classifiers achieved average precision of 79.08%, 77.20% and 82.79%, recall of 79.09%, 78.25% and 83.23% and an accuracy of 79%, 78% and 83% respectively. In comparison, the proposed particle filter resulted in an average precision of 85.79%, recall of 86.11% and an accuracy of 86% respectively. Furthermore, MCC [32] was computed which achieved values of 71.12%, 70.51%, 76.44% and 80.08% for SVM, KNN, Random Forest and PF respectively. Evidently, the achieved multi-class diagnosis results using diagnostic radiomic features and particle filter are better when compared with the trained SVM, KNN and Random Forest classifier on low dose datasets. The proposed filter has classified 86% of the nodules correctly with fewer input longitudinal data points proving particle filter a good choice for multi-class classification. However, the classification results could be improved with the availability of more chest CT data but is non-attainable due to health hazards attached with multiple CT screenings of a patient. Moreover, a limitation faced to implement the proposed approach is poor management of medical images database currently at many hospitals.
With proper management of patients' longitudinal CT data, this research work can be extended to cancer prognostics. Besides current status, the future malignancy status of a nodule can also be estimated using PF approach which re-iterates its efficacy.

VIII. CONCLUSION
An accurate and efficient radiomics based Bayesian inversion framework is proposed to detect cancer and its pathological stage as early or advanced stage. This contribution is superior to other prominent cancer detection algorithms whose performances are dependent on large data sets training, a limitation in medical imaging. Moreover, the published contemporary methods do not provide a complete picture of malignancy and its stage using a single research approach as presented in this research work. Furthermore, highly discriminative features towards cancer stages were also identified. The high accuracy of the presented method shows that the proposed model can be used by clinicians in cancer diagnostics.