2.1.

Patient Selection, Annotation, and Chest Computed Tomography Acquisition

In this retrospective study, we included patients from two sites, who had a suspicious lung nodule on CT scan and underwent surgical resection. All scans acquired for this study were collected as part of an Institutional Review Board-approved, HIPAA-compliant protocol. These scans were obtained as part of standard of care clinical management for these patients, all of whom subsequently went on to have a surgical wedge resection for excision of a suspicious nodule. All the information provided was deidentified and the need for an informed consent was waived. Histology was confirmed by an anatomical pathologist based off visual interrogation of the surgical specimen.

One hunderd and thirty-nine solitary nodules were obtained from site 1 (70 carcinomas and 69 granulomas) and were used as the training set. Fifty six solitary nodules (34 carcinomas and 22 granulomas) were acquired from site 2 and composed of the independent test set. Patients with multiple solitary nodules were excluded. Figure 1 illustrates the inclusion and exclusion criteria.

Fig. 1

Inclusion and exclusion criteria for patient selection in two cohorts.

2.2.

Lesion Segmentation and Computer-Extracted Texture and Shape Features

The region of interest (RoI) containing the lesion was manually segmented across contiguous slices by an expert cardiothoracic radiologist with 20 years of experience in interpreting chest CT scans via a hand-annotation tool in 3D-Slicer^® software.¹⁹ To evaluate the interobserver variability in the lesion segmentation and also to show the effect of it on the stability of the top radiomics, a subset of 10 adenocarcinoma and 10 granuloma lesions is randomly picked and resegmented by an expert radiologist with 14 years of experience in cardiothoracic radiology. Note that the second reader was completely blinded to the segmentation of the first reader. A total of 645 two-dimensional (2-D) texture and 24 three-dimensional (3-D) shape features were extracted from the lesion area. Texture features were extracted in 2-D instead of 3-D, since the available retrospective CT volumes were all anisotropic. After extracting per-voxel based features within the nodule of interest, five statistics (mean, variance, minimum, maximum, and entropy) were calculated for each nodule. All feature calculations were implemented using MATLAB^® 2014b platform (Mathworks, Natick, Massachusetts). The description of the extracted features is given in the following subsections.

2.3.

Stability Evaluation of Radiomic Features

While feature discriminability is critical to be able to distinguish the two pathologies of interest, an important consideration in constructing generalizable and stable classifiers is to ensure that the features invoked in constructing the machine learning classifier are relatively robust across sites and slice thickness for a specific pathology (i.e., granulomas or adenocarcinomas). The criticality of invoking stable features has been well documented in the radiology domain,^27–29 however, relatively little attention has to be paid to this issue in the context of radiomics.^30,31 In this work, we also investigated the sensitivity of radiomic feature expression to variations in slice thickness and across the different acquisition sites and scanners. In this paper, we use the preparation-induced instability (PI) score defined by Leo et al.³⁰ to measure robustness of the feature value. The PI value is a number between 0 and 1 and calculated per radiomics. The PI closer to 0 implies that the corresponding radiomics is more stable between two cohorts of interests.

In this paper, the PI score was calculated between the testing cohort and the training cohort to quantitatively demonstrate the robustness of the top radiomic features. Additionally, the PI score was calculated for a subset of the training data to investigate the variance in radiomic features between scans with a slice thickness of less than 3 mm versus those scans with a slice thickness more than 3 mm.

To provide a better evaluation and more valid conclusion, a subset of 20 lesions (10 from each pathology) is randomly selected and segmented by an independent cardiothoracic radiologist. The Dice coefficient,³² undersegmentation error, and oversegmentation error between the segmentation of the readers, and also the area under the curve (AUC) of diagnostics of each one, was reported as a measure of inter-reader variability. In addition, the PI value of top radiomics, from the segmentation of two experts, is expressed as a measure of radiomics stability versus interreader variability.

2.4.

Statistical Analysis

We employed a consensus clustering approach to graphically compare the within-class and between-classes correlation for the features identified as discriminatory by feature selection within the training set.³³ To perform the consensus clustering, the similarity between different nodules, assessed according to a distance in the space of the top ranked features, was first determined. The closer two nodules were to each other, the higher the likelihood that they both belong to the same cluster. Clusters were determined by invoking the idea that nodules within a cluster should have a high intraclass correlation, whereas nodules belonging to different clusters should have minimal correlation. Figure 2(a) shows the consensus clustering result generated for the nodules in the training set within the space of the top ranked radiomic features. Figure 2(b) shows the corresponding consensus clustering obtained for the nodules in the training set when only invoking the Haralick texture features.

Fig. 2

Consensus clustering of the cases in the training set when the data were clustered into six partitions. The clustering was done (a) on the most predictive shape and texture features identified during feature selection and (b) for Haralick texture features alone. While the result of clustering in panel (a) reveals two primarily distinct clusters of adenocarcinomas and granulomas and (b) reveals multiple disjointed clusters of varying sizes. The comparison of (a) and (b) shows good inherent discriminability between adenocarcinoma and granuloma when a combination of the most predictive shape and texture features is considered, as opposed to Haralick features alone.

2.5.

Features-Based Classification and Validation on Independent Test Set

Following feature extraction, sequential forward feature selection was employed to identify the most discriminating subset of features.³⁴ To mitigate bias in feature selection and classifier training, a threefold (onefold held-out for testing), patient-stratified, cross-validation scheme was used for classifier constructing using the instances within the training set, and cross validation was repeated a 1000 times.

The top six features identified by the feature selection approach were then used to train and lock down linear discriminant analysis (LDA), quadratic discriminant analysis (QDA), and support vector machine (SVM) classifiers. These classifiers were then applied to predict the class labels of the nodules in the validation set. Each nodule was assigned a probability of being an adenocarcinoma. Area under the receiver operating characteristic (AUC) curve was used to evaluate the performance for each of the SVM, LDA, and QDA classifiers on the validation set. A Bayesian cross-validation extensive search on the normalized data was employed to find the best box constraint and the sigma value of radial basis function (RBF) kernel of SVM classifier. The employed box constraint and sigma are 0.431 and 14.07, respectively.

2.6.

Human–Machine Comparison

The classification performance of our SVM classifier was compared against the diagnosis of two experts, one a board certified attending radiologist with 12 years of experience in thoracic radiology and a pulmonologist with 2 years of experience in reading chest CT scans. Both experts were blinded to the true histopathologic diagnosis of the 56 cases, which comprised the test set. Each reader was asked to assign a score between 1 and 5 to each nodule, with 1 referring to a high confidence that the nodule is “benign,” 2 referring to a diagnosis of “mostly benign,” 3 being “not sure,” 4 being “mostly malignant,” and 5 being “malignant.” To evaluate the performance of the experts, the classifier probability output was compared to diagnostic ground truth determined from the pathology reports. From these comparisons, a receiver operating characteristic (ROC) curve was obtained and the AUC was calculated. The AUC for the human readers on the test set was then compared against the corresponding values for the SVM classifiers.

3.1.

Patient Characteristics

Table 3 illustrates the patient characteristics in the two patient cohorts considered in this study.

Table 3

The patient characteristics in the two cohorts.

A single- and double-tailed paired student’s $t$-test was applied between the lesion size of adenocarcinoma and granuloma in the testing set and the training set. The results of this comparison are reported in Table 4. With the double-tailed test, no significant difference between the lesion size of two classes across the two cohorts was observed, though a weak correlation ($p=0.03$) was observed for the single-tailed paired student’s $t$-test.

Table 4

Results of performing a one-tail and two-tail paired student’s t-test between the lesion diameters of two classes in testing and training cohort.

3.2.

Top Features Identified by Feature Selection

The top features selected during the cross-validation process in the training phase of cohort-1 and three different classifiers (LDA, QDA, and SVM) are shown in Table 5. Based on empirical evaluation, the RBF kernel was considered the optimal kernel for use in conjunction with the SVM, and polynomial kernels of degrees 1, 2, and 3 were also evaluated. The Table suggests that among texture features, “skewness of Laws features (L5 × E5) and (L5 × R5),” “skewness of gradient features (Sobel in both direction and diagonal),” and “Gabor texture features” were the most consistent and top performing texture features across the three different classifiers. The top performing shape features were identified as “mean of extend,” “mean of convexity,” and “variance of eccentricity,” though the texture features substantially outperformed the shape features. The most informative statistics were skewness and variance, and none of the top features include kurtosis or entropy of texture features. As a new experiment, without causing any contamination to the validation set for subsequent experiments, we swapped the training and the testing cohorts. In this new experiment, three classifiers—SVM, LDA, and QDA—were trained on the new training (previously validation) and testing (previously training) sets. The top features for this experiment are listed in Table 6. While the top features for this experiment are not identical to the top features identified using the original training set, the same feature classes (Gabor and Haralick) and statistics (skewness and kurtosis) were again represented.

Table 5

The most discriminative texture and shape features (based off AUC) when training the LDA, QDA, and SVM classifiers with cohort-1 (N=139) and testing on cohort-2.

Table 6

The most discriminative texture and shape features (based off AUC) when training the LDA, QDA, and SVM classifiers with cohort-2 (N=56) and testing on cohort-1.

Figure 3 visually illustrates the discriminability of the top two features, i.e., “energy” and “Gabor $S=1/4$, $\mathrm{ϴ}=3\pi /8$” for a granuloma and an adenocarcinoma. The texture heat maps of the nodule in Fig. 3 appear to suggest more heterogeneity in adenocarcinomas compared to granulomas, patterns not immediately obvious on the original CT scans in Fig. 3. This trend also appears reflected in the box and whisker plots shown in Fig. 4, where the average Hounsfield units of the nodules within the CT scans do not appear to show any clear separation between granulomas and adenocarcinomas. Note that the corresponding texture features shown in Fig. 3 do appear to show statistically significant separation on the validation set.

Fig. 3

An illustration of the discriminability of the texture features. Despite the homogeneity of HU in CT of the adenocarcinoma and granuloma, the texture features demonstrate informative heterogeneity in adenocarcinoma in comparison to the granuloma.

Fig. 4

Box and whisker plots corresponding to the mean of (a) CT Hounsfield units, (b) variance of energy, and (c) variance of Gabor texture features extracted from within the nodules in the training set.

Figure 5 shows the AUC values for the top features identified by the sequential feed-forward feature selection algorithm for the cases in the training set. As the number of employed features increases, the AUC value increases. However, to avoid overfitting,³⁵ the cardinality of the feature set was restricted to 10 features. Figure 6 illustrates the 2-D scatter plot of top features, sum of standard deviation and skewness of Laws L5 × E5 texture features for the adenocarcinomas and granulomas on the training set. The green line suggests the false-negative free domain for adenocarcinoma (equivalent to false omission rate = 0) in the training set. In this scenario, the classifier can achieve a positive predictive value as high as 72%.

Fig. 5

The AUC values (with the standard deviation bars) for three different machine learning classifiers (LDA, QDA, and SVM) versus the cardinality of the top features subset achieved via sequential feature selection method on the training set ($N=139$).

Fig. 6

A 2-D scatter plot of adenocarcinomas (red dots) and granulomas (black dots) in the training set plotted in the space of the top two texture features (Laws L5 × E5 and sum of standard deviation) identified by the feature selection algorithm. The green line was identified as the optimal linear boundary separating the granulomas from the adenocarcinomas in the training set.

3.3.

PI Stability Evaluation of Top Radiomic Features

As a measure of segmentation agreements, the Dice coefficient between the segmentation of two radiologists, $R1$ and $R2$, is defined as follows:

The results of the stability experiments for the top ranked radiomic features are illustrated in Table 7. Table 7 suggests that the top discriminating features are not necessarily the most stable and reproducible features across the different sites and between scans with greater than and less than 3-mm slice thickness. Interestingly, the shape feature “mean of extend” appeared to be the most stable feature while the “variance of sum variance” feature was the most stable texture feature.

Table 7

The values of PI instability measure for top ranked radiomics.

The stability of the top radiomics features across the segmentation of two radiologists is also measured by PI values. The PI values of the most stable shape feature and the most stable texture feature are 0.891 and 0.924, respectively. The comparison between the PI values of the most stable radiomics in the different sites and thickness with the PI values for different readers reveals that the stable features across sites and thickness were remained stable across different readers.

3.4.

Statistical Analysis

Figure 2 shows the consensus clustering for radiomic features extracted from the cases in the training set. This figure suggests that the combination of top texture and shape features yields two distinct clusters, the individual clusters corresponding almost exclusively to either granulomas or adenocarcinomas. However, the corresponding consensus clustering results obtained using the Haralick features employed by Dennie et al.¹⁸ yield multiple disjointed clusters of varying sizes.

Figure 7(a) illustrates the discriminability of the combination of texture and shape features for the LDA, QDA, and SVM classifiers in the training set. Results for these three different classifiers on the testing set are shown in Fig. 8. The best AUC, corresponding to the SVM classifier, on the training set for a combination of four texture and two shape features was $92.9\%\pm 1.14\%$. On the independent test set, the resulting AUC of the locked down SVM classifier was 77.8%. When the training and testing sets were swapped, the QDA classifier was identified as the top ranked classifier with an $\mathrm{AUC}=82.5\%$ on the test set. To extend the generalizability of the most informative radiomics, their stability between the training and testing cohorts and between different slice thicknesses was evaluated.

Fig. 7

(a) Receiver operating characteristic curves and corresponding AUC values for the LDA, QDA, and SVM classifiers for discriminating adenocarcinoma from granulomas on the training set on cohort-1 ($n=139$). The left panel shows a zoomed in version of the ROC curves for the different classifiers in the specificity and sensitivity range of 60% to 100%. (b) Receiver operating characteristic curves and corresponding AUC values for the training set on cohort-2.

Fig. 8

(a) Receiver operating characteristic curves and corresponding AUC values using the LDA, QDA, and SVM classifiers for discriminating adenocarcinomas from granulomas on the test set ($n=56$). (b) Receiver operating characteristic curves and corresponding AUC values for training on cohort-2 and testing on cohort-1.

For the human–machine comparison in the same holdout set, the AUCs for an attending radiologist with over 12 years of experience and a pulmonology fellow were found to be 69.72% and 72.39%, respectively. Pearson’s correlation and single- and double-tailed paired student’s $t$-test were performed between the prediction results on the test set of the SVM classifier and reader 1 (radiologist). However, no statistically significant differences were found. The results of this comparison are reported in Table 8.

Table 8

Results of performing a Pearson correlation and one-tail and two-tail paired student’s t-test between the prediction results of the machine classifier (support vector machine) with the two human readers (radiologist and pulmonology fellow) on the test set (n=56).