2.1.

Databases

DCE-MR images of breast lesions were retrospectively collected under Internal Review Board and HIPAA compliance from an institution in the United States (during the period of 2005 to 2017) and in China (during the period of 2015 to 2017) (Table 1). These images constituted the two “populations.” Images collected from the United States were acquired in the axial plane, whereas images collected from China were acquired in the sagittal plane. The time interval between postcontrast images was typically 60 s in the United States database and typically 90 s in the China database. Most images collected in the United States were acquired using Philips scanners (Best, The Netherlands) (1145 out of 1163 cases), whereas all images collected in the China database were acquired using GE Discovery 750 scanners (Waukesha, Wisconsin).

Table 1

Composition of the databases: Number of lesions by status as benign or cancer, as well as maximum linear size and age.

2.2.

Computerized Breast Lesion Segmentation and Radiomic Feature Extraction

Lesions were segmented using a fuzzy c-means method.⁶ Thirty-two computer-extracted radiomic features were automatically calculated, covering categories of size, shape, morphology, texture, and kinetics of contrast dynamics.^7–9 Features that were candidates for harmonization were identified for their dependence upon gray level and potential for intrinsic variation related to gain, image resolution, image noise, and imaging protocol. Features describing lesion geometry (i.e., categories of size and shape) and the kinetic curve assessment feature of volume of most enhancing voxels were deemed to not be candidates for harmonization due to these being characteristics directly of the lesions and not dependent on the image acquisition. Additionally, two other kinetic curve features, washout rate and curve shape index, were excluded from harmonization because they are semi-categorical variables, in which the algorithm sets the feature value to zero if criteria for indication of washout of contrast agent are not met, as is true for many benign lesions.⁹ Because features were harmonized by feature category, these two features were not included in the harmonization of the other kinetic curve features even when values were not set to zero, as doing so would have resulted in unmatched sets of features of cases after removal of features that were equal to zero. Based upon the consideration of features and their description, features that were candidates for harmonization were determined, whereas the other features did not undergo harmonization (Table 2).

Table 2

Description of radiomic features.

2.3.

Radiomic Feature Harmonization

Harmonization was conducted using the parametric version of the ComBat method, which is based on additive and multiplicative batch effects using empirical Bayes estimates, with which to transform the features.⁵ In our study, the two populations (i.e., patients imaged in the United States and patients imaged in China) served as the two “batches,” yielding the source of differences the ComBat algorithm seeks to reduce. The harmonization was applied separately within categories features, i.e., morphology, texture, and kinetics, as described above. The ComBat harmonization method can accommodate covariates, i.e., confounding variables, which for our study was malignancy status of the lesions. The cancer prevalence was similar between the two populations (see Table 1), an important consideration for the implementation of ComBat harmonization for classification of lesions as benign or malignant.¹⁰

2.4.

Effect of Harmonization on Radiomic Features

The effect of harmonization on feature value distributions in each population in benign lesions and in cancers was evaluated for each eligible feature using the Kolmogorov–Smirnov (K-S) test to compare the distributions.^11,12 The K-S test statistic was used to characterize the change in feature value distribution due to harmonization.

To visualize the impact of harmonization methods on clustering of features by population, t-distributed stochastic neighbor embedding (t-SNE) methods were used to reduce the dimensionality of the feature sets from 32 to two.¹³ The Davies–Bouldin index¹⁴ was used after k-means clustering of the t-SNE values to assess inter- and intra-cluster agreement across the two populations, separately for benign lesions and for cancers.

2.5.

Effect of Harmonization on CAD Performance

Classification of lesions as benign or malignant was performed using two sets of 32 radiomic features: (a) a feature set with all features in their original form (called as “preharmonization” in this work) and (b) a feature set comprised of features not eligible for harmonization combined with those that had undergone harmonization in their respective feature groups (called as “postharmonization” in this work).

The classification performance evaluations were performed separately for lesions imaged in the United States, for lesions imaged in China, and the combined dataset. For each classification evaluation on each set, 10-fold cross validation was performed using a random forest classifier with 100 decision trees, with the posterior probability of malignancy used as the classifier output for receiver operating characteristic (ROC) curve analysis. No feature selection was conducted prior to the use of the random forest classifier. The area under the ROC curve (AUC), determined using the proper binormal model,¹⁵ served as figure of merit, with its value and 95% CI determined using ROCkit software.¹⁶ The difference in AUC for each set of cases was deemed to be statistically significantly different if $p<0.05$.

3.1.

Effect of Harmonization on Radiomic Features

The harmonization method changed the distribution of values of the features that had been determined to be eligible for harmonization (Fig. 1). In general, for each radiomic feature eligible for harmonization, the difference in median feature value for benign lesions and for cancers was reduced after harmonization, and the distributions became more similar. The K-S test statistic for comparing feature value distributions pre- and postharmonization showed that in benign lesions, larger changes in most feature value distributions were seen in the lesions imaged in the United States. For cancers, some feature value distributions demonstrated more change in feature values in lesions imaged in China (features T2, T4, T9, T10, T12, and T13), and others demonstrated more change in feature values in lesions imaged in the United States (mostly, morphology and kinetic curve features, as well as the remaining texture features, Fig. 2). T-SNE space visualized pre- and postharmonization in benign lesions and in cancers (Fig. 3) show that the degree of clustering between populations was reduced after harmonization was applied to eligible features. The Davies–Bouldin index for degree of interclustering and intraclustering of t-SNE values increased after harmonization was applied to eligible features, demonstrating that the features were more similar between populations postharmonization (Fig. 4). The Davies–Bouldin index for benign lesions increased by 27%, whereas it increased for cancers by 43%.

Fig. 1

Example raincloud plots for the distribution of the values of (a) the texture feature of difference entropy and (b) the kinetic curve feature of enhancement at first postcontrast time point, before and after harmonization for benign lesions and cancers within each population. The black circles at the base of each distribution indicate the median of the distribution, and the connecting lines between the medians of the distributions pre- and postharmonization demonstrate the change of the median with harmonization of the feature. The dots below each distribution show individual data points.

Fig. 2

The K-S test statistic for feature value distributions compared pre- and postharmonization, for ($x$ axis) lesions imaged in the United States and ($y$ axis) lesions imaged in China [(a) benign lesions and (b) cancers]. Data points below the diagonal indicate that there is a larger difference between the feature value distributions compared pre- and postharmonization in lesions imaged in the United States, and data points above the diagonal indicate that there is a larger difference between the feature value distributions compared pre- and postharmonization in lesions imaged in China. Feature abbreviations are the same as in Table 2.

Fig. 3

t-SNE space for the 32 features reduced to 2 features for (a), (c) preharmonization state and (b), (d) postharmonization. (a), (b) t-SNE space for benign features and (c), (d) t-SNE space for cancers.

Fig. 4

The Davies–Bouldin index, a measure of inter- and intra-clustering, for the t-SNE values (32 features reduced to 2) for benign lesions and cancers. A higher Davies–Bouldin index indicates improved harmonization across the features between populations, as is also seen in Fig. 3.

3.2.

Effect of Harmonization on CAD Performance

In the task of classification of lesions as benign or malignant, classification performances, as measured by AUC, demonstrated statistically significant increase when postharmonization features were used compared with preharmonization features (Fig. 5, Table 3)

Fig. 5

ROC curves in the task of classification of lesions as benign or malignant, for the three datasets pre- and postharmonization of selected radiomic features. Within each dataset, ROC analysis using the posterior probability of malignancy was performed after 10-fold cross validation with random forest classifier.

Table 3

The AUC in the task of classification of lesions as benign or malignant, for classification using preharmonization features and for classification using postharmonization features. Within each dataset, the AUC was determined using classification with random forest classifier via 10-fold cross validation, with the posterior probability of malignancy used as the classifier output for ROC curve analysis.