Radiomic analysis as a field is predicated on the idea that radiological imaging contains meaningful biological information contained in the patterns of intensity values within regions of interest, which could contribute to a better understanding of patient health (Gillies et al., 2016). The typical approach in radiomic studies is to extract a large number of pre-defined quantitative imaging features from a region of interest, and then use machine learning methods to reduce the dimensionality of the feature set and build models of a biological correlate or a medical outcome of interest (Horvat et al., 2022). The promise of radiomics to develop quantitative imaging biomarkers is of broad interest because it poses a non-invasive means to characterize patient disease using routine clinical imaging. However, to be clinically deployed, radiomic models must be widely validated, and their robustness to variable imaging settings well-established.

In the present study, we are focused on contrast-enhanced abdominal computed tomography (CT) of the liver for patients with colorectal liver metastases (CRLM). These patients have an overall poor prognosis, which could potentially be improved by prognostic radiomic signatures that could better target patients for surgery or chemotherapy. In patients with liver metastases, radiomic models derived from contrast-enhanced CT have shown substantial prognostic capability in both survival modeling and prediction of chemotherapy response (Fiz et al., 2020). Furthermore, although most studies focus on radiomic features extracted from the metastases themselves, studies of CRLM have shown that features from the liver parenchyma also contain important information when predicting hepatic disease-free survival or overall survival after hepatic resection (Simpson et al., 2017), or progression-free survival after radiotherapy (Hu et al., 2022b). With these applications in mind, in this study we are seeking to gain a better understanding of the robustness of radiomic models in contrast-enhanced CT of CRLM.

One way to study the robustness of radiomic models is to study the reproducibility of the model inputs—that is, to understand how consistent the radiomic features are under real-world variations that occur in image acquisition and reconstruction. Treating the radiomic features as measurements drawn from radiological images, the reproducibility of the features can be studied from a metrological perspective (Raunig et al., 2014), using statistical measures such as the concordance correlation coefficient (CCC). Understanding the reproducibility of radiomic features, however, is difficult given the many factors affecting the results. Restricting the discussion to CT, scans are acquired from different medical centers, using different scanners, with different imaging acquisition parameters and protocols, after which the images are reconstructed with different algorithms, different resolution parameters, and different kernels, all of which have been shown to affect feature reproducibility (Zhao, 2021). Once the images are in hand, feature extraction itself is not without concerns regarding reproducibility; while standardization efforts such as the Image Biomarker Standardisation Initiative (IBSI) (Zwanenburg et al., 2020) have been instrumental in creating a common nomenclature and a standard set of well-defined features, they do not provide a final answer in how to set the many configurable parameters which can affect the final feature values computed by common software such as pyradiomics (van Griethuysen et al., 2017).

In this study, we are considering feature reproducibility particularly when the slice thickness used to reconstruct the CT scans is varied. When comparing features extracted from CT images with thinner or thicker slices, studies have found that features from thinner images are more reproducible across variations in segmentation (Hu et al., 2022a), or across repeat imaging (Zhao et al., 2016), and may also produce more accurate models (He et al., 2016; Li et al., 2018; Xu et al., 2022). In direct comparisons, many features had poor reproducibility when comparing features from images with different slice thicknesses in a variety of phantom studies (Zhao et al., 2014; Ger et al., 2018; Berenguer et al., 2018; Kim et al., 2019; Varghese et al., 2019; Ligero et al., 2021; Ibrahim et al., 2022). Prospective studies producing multiple reconstructions for each patient have reproduced this result on patient images for lung cancer (Lu et al., 2016; Park et al., 2019; Erdal et al., 2020; Yang et al., 2021; Emaminejad et al., 2021), and liver metastases (Meyer et al., 2019). Poor feature reproducibility with respect to slice thickness is concerning because many real-world retrospective or multi-site data sets include images with a range of different slice thicknesses due to variations in local protocols (Ger et al., 2018). Although image interpolation to a common, isotropic voxel size is considered a best practice for preprocessing during feature extraction (Zwanenburg et al., 2020) in order to ensure the image features are comparable between images with different voxel sizes, the optimal choice of resolution and resampling algorithm is undecided. Furthermore, resampling to a common voxel size appears, on its own, to be insufficient to overcome the inconsistency of feature values due to slice thickness variation, except in a small subset of features (Shafiq-Ul-Hassan et al., 2017; Shafiq-ul Hassan et al., 2018).

To further complicate matters, it has been shown that reproducibility is not necessarily consistent across cancer types, even for a single modality such as CT (van Timmeren et al., 2016). However, in a systematic review of radiomics reproducibility studies, Traverso et al. (2018) found that the literature was limited to a small number of cancer types, with the greatest number of studies addressing lung cancers, amongst which CT was the most common imaging modality. Across these studies there was not a clear consensus on the most reproducible features, although in CT they found agreement that first-order features tended to be more reproducible than higher-order texture features. Ultimately, it seems that the reproducibility of features is not easily generalized across different anatomies or cancers. Therefore, reproducibility studies for the region of interest (ROI) and disease under consideration is an important part of the validation any radiomics-based imaging biomarker.

In this study we present an analysis of the relationship between the reproducibility and prognostic value of radiomic features drawn from contrast enhanced CT of CRLM. We present a reproducibility analysis of radiomic features on a cohort of 81 prospectively enrolled patients from two major US cancer centers, who underwent contrast enhanced abdominal imaging with a controlled and systematically varied protocol. Our analysis is primarily focused on the effects of slice thickness chosen at reconstruction time. Features were extracted from the largest liver metastasis and the liver parenchyma from each subject using a variety of different configurations, varying the level of resampling, and the method of aggregation used in computing the higher-order texture features. To investigate the relationship between reproducibility and prognostic value, we conducted an in depth univariate and multivariable survival modeling analysis on an independent, publicly available data set of 197 preoperative contrast enhanced CT scans of patients who underwent hepatic resection to treat CRLM (Simpson et al., 2023, 2024).

This paper is a significant expansion of our previous work (Peoples et al., 2023), with a greater focus on the relationship between the reproducibility and prognostic value of the radiomic features under consideration. In this paper, we are less concerned with finding the feature extraction settings that produce the most reproducible features. Instead, we focus on the integration of reproducibility information into the development of prognostic radiomic signatures. We present a joint, univariate analysis of both the reproducibility and prognostic discriminative ability of the features, taking a multi-objective optimization point of view. The multivariable analysis was revised to use a standard feature selection algorithm, and to conduct many iterations of the cross-validation to ensure our results were stable. Additional details are provided throughout the paper, giving greater context on the methodology, more visualisation and discussion of the results, and more interpretation of how our results fit into the greater context of the literature on the reproducibility of radiomics. Code for all analysis is available at github.com/jpeoples/melba2024.

The results of the first phase of the reproducibility analysis are summarized in Figure 3. Here, the distribution of CCCs across all feature values for each extractor is compared across each pair of slice thicknesses. Intuitively, the largest difference in slice thickness (2.5 mm vs. 5 mm) shows a significantly lower distribution of CCCs than the other two pairs, which compare slice thicknesses which are closer together (2.5 mm vs. 3.75 mm, and 3.75 mm vs 5 mm). These differences were highly statistically significant based on a Wilcoxon sign-rank test between all CCC pairs ($p<3.8\times 10^{-30}$ in the bottom two rows of the table in Figure 3). Interestingly, the 3.75 mm vs. 5 mm pairs also seem to have slightly higher CCC than the 2.5 mm vs. 3.75 mm pairs, at a lower level of significance ($p<6.4\times 10^{-4}$ in the top row of the table in Figure 3).

We now turn to the second phase of the reproducibility analysis, based on the generalized CCCs computed across all slice thicknesses using the LMM model. Figure 4 (A) summarizes the results across all features, extraction settings, and ROI. The results are visualized as a heat map, where each row is a feature, and each column corresponds to a combination of ROI and extractor setting, and the values are CCCs. Both the rows, and columns underwent hierarchical clustering based on Euclidean distance and Ward linkage. We can see, near the top of the map, a cluster of features that are highly reproducible across all extraction settings. As can be seen in the left-most column of the map, these highly reproducible features all belong to the first-order feature class. With some exceptions, within both ROI, the CCCs across extractors tend to follow similar patterns; however, overall, the liver features visually appear to tend toward lower CCCs. This tendency is confirmed by a Wilcoxon sign-rank test, which shows that feature CCCs between the two ROI are statistically significantly lower in the liver parenchyma, with $p\in[8.8\times 10^{-11},3.6\times 10^{-5}]$ across all extractors. Interestingly, the CCCs between ROI and extractor settings are sufficiently different that they cluster logically: the two ROI form two higher level clusters; and within each ROI, L2i and S2i form their own cluster, and within the other cluster, the 2.5D and 3D settings cluster together.

Figure 4 (B) and (C) show similar cluster maps of the univariate C-indexes, computed against overall survival in the publicly available cohort, for each feature, from the tumor and liver ROI, respectively. The ROI were clustered separately because the ROI are biologically different, and therefore the discriminative ability of a given feature may be completely different in a different ROI. Both maps cluster features into three broad groups—a low C-index group at the top, a high C-index group in the middle, and a middle group at the bottom. Both the high C-index and low C-index groups are larger in the tumor ROI than in the liver parenchyma. Broadly speaking, the features tend to show similar trends in terms of predictive value across the extractor settings. On the other hand, there are some exceptions: at the bottom of the tumor cluster map (B), there is a cluster of features where the L2i and S2i features tend to be less predictive than the others. Similarly, in the high C-index cluster of the liver ROI, L2i and S2i again tend to be less predictive. The column clusters in both maps reproduce the finding for CCC in Figure 4 (A): three high-level clusters corresponding to (L2i and S2i), (A2, L2, and S2), and (A3, L3, and S3).

To examine the relationship between feature reproducibility and discriminative ability, we considered the relationships between the CCC and C-index for each feature. For every feature, we checked which extractors produced the highest CCC, the highest C-index, and which extractors were Pareto efficient for both CCC and C-index in that feature. The results are summarized in Figure 5 (A)–(C). Comparing Figure 5 (A) and (B), we see that the results are inconsistent—there is no set of extractors that produce both the highest CCC and the highest C-index. With the exceptions of A2 and A3, Figure 5 (C) shows that all feature extraction settings are well-represented on the Pareto fronts across features. Finally, grouping all features from all extractors together, and computing the set of Pareto efficient features for CCC and C-index shows that all extractors except A2 contribute (see Figure 5 (D)), with S3 contributing the greatest number of features (6). A scatter plot of all features is shown, color coded by extractor setting, with the Pareto front highlighted, in Figure 5 (E). Finally, the set of all 23 Pareto efficient features is listed in Table 5 in the appendix. It is worth noting that several features that are listed are actually equal, and therefore over-counted in Figure 5 (D). Furthermore, note that first-order features across S2 and S3, or L2 and L3 are also equal, and therefore some features that appear in both sets are still identical. These groups of equal features are marked in Table 5. Figure 5 (A)–(D) also indicates that the optimal feature extraction method for a given criterion (C-index, CCC, or Pareto front) appears to vary not only across specific features, but across ROI. Indeed, most of the extractors, and both ROI are represented, at least to some extent, on the overall Pareto front across all features and extractors (see Figure 5 (D)).

The multivariable CPH cross-validation experiments produced results for 315 combinations of three parameters (feature extraction setting (8 extractors plus 1 for all), feature count (1, 2, 4, 8, 16, 32, 64), CCC threshold (0, 0.8, 0.85, 0.9, 0.95)). The test-set C-index and 95% confidence intervals across all 100 repetitions are listed for the ten highest performing parameter combinations in Table 4. Models based on the features from extractor L2i attained the highest C-index of 0.630 (0.603–0.649), with 4 features being selected, and no CCC thresholding. The second highest performance of 0.629 (0.605–0.645) was attained using features from all extractors, with 4 features being selected, and a threshold of $\mathrm{CCC}\geq 0.85$. Based on these top two performing cases, the performance of the models using 4 features, across all extractors, and $\mathrm{CCC_{t}}=0$ and $\mathrm{CCC}_{t}=0.85$ are visualized in Figure 6. For most extractors, the thresholding of features to those with CCC over 0.85 appears to have little effect on the performance of the resulting models. In a few cases, L2i, A2, and A3, the performance is reduced after thresholding. When including features from all extractors (left-most in the figure), the performance increases after CCC thresholding. The performance on the training set tends to be lower after thresholding, suggesting that the models using reproducible features may be less prone to overfitting.

To visualize the sensitivity of these results with respect to the CCC threshold and the number of features selected, we have visualized the results from all parameter settings in Figure 7. Each plot shows the C-index for both the training and test sets, plotted against the number of features selected in the model. In some cases, after thresholding, and applying the univariate feature filter (removing features with C-index below 0.55 or with $p$-value over 0.1), there may have been fewer than the desired number of features remaining for feature selection. In these cases, all features passing the threshold and univariate filter were used. The 95% confidence intervals for the number of available features at the time of applying mRMR for the different threshold levels is indicated by the gray regions, where either limit was less than 64. Typically, beyond this point, requesting more features does not change the performance of either the train or test set, because the number of features is saturated. This effect is clearly visible in the plots for $\mathrm{CCC}\geq 0.9$ and $\mathrm{CCC}\geq 0.95$, where fewer features tend to remain after thresholding, due to the stringent reproducibility requirement. Broadly, we observe that the best performance tends to be obtained using 4 or 8 features, insofar as there are sufficiently many features available after thresholding. Again we can observe that, although the effect appears small, the training set performance is slightly reduced when the features are thresholded on $\mathrm{CCC}$, indicating a potential reduction in overfitting. For higher thresholds, this effect is primarily due to the number of features being saturated, but for $\mathrm{CCC}\geq 0.8$ and $\mathrm{CCC}\geq 0.85$ the number of available features is high enough to observe the difference even before feature count saturation.

A visualization of tumors with different values for features with high/low C-index (prognostic value), and high/low CCC (reproducibility) is given in Figure 8.

Although in the present study we have restricted our analysis to the reproducibility of radiomic features when slice thickness is varied at reconstruction time, the ultimate purpose of the data collection effort from which this study draws is to look at the effects of contrast timing and reconstruction parameters using a test-retest paradigm with two portal venous phase images collected within 15s of each other. As part of a multicenter prospective study systematically varying contrast timing and reconstruction parameters, this paper presents early results from a research effort that will be a substantial step forward in the understanding of the reproducibility of radiomics for contrast-enhanced CT of CRLM.

Our results show that when comparing the 20% ASiR images pairwise at different slice thicknesses, the CCCs for features from the 5 mm and 2.5 mm images are significantly lower than when we compare 5 mm and 3.75 mm, or 2.5 and 3.75 mm pairs. The degradation in consistency of features as the change in slice thickness increases is consistent with findings in the literature from phantom studies (Zhao et al., 2014; Kim et al., 2019; Ligero et al., 2021), as well as lung cancer CT (Lu et al., 2016; Park et al., 2019; Erdal et al., 2020). The consistency of this finding across cancers is noteworthy given that feature reproducibility is anatomy/disease-specific, even within CT (van Timmeren et al., 2016). We also found, at a lower level of statistical significance, that the features had higher CCC in the 5 mm and 3.75 mm pairs, compared to the 2.5 mm and 3.75 mm pairs. We hypothesize that this result could be due to the lower noise level in thicker sliced images, compared to thinner slices, or due to resampling issues when upsampling the ROI segmentations to 2.5 mm.

Moving to liver cancer more specifically, Perrin et al. (2018) examined a database of consecutive patients from MSK diagnosed with liver malignancy, with the additional inclusion requirement that they had two contrast-enhanced abdominal CT scans taken within no more than 14 days of each other. By including two scans per patient with a small separation in time, this data set served as an approximation of a test-retest study, allowing the study of the reproducibility of radiomic features across two consecutive scans. The database included patients with multiple liver cancers including liver metastases (n=22), intrahepatic cholangiocarcinoma (n=10), and hepatocellular carcinoma (n=6). Because image acquisition and reconstruction parameters naturally varied between the scans for a given patient, this data set allowed the exploration of the effects of these different parameters on the overall concordance of the radiomic features between the consecutive scans, as measured by CCC. One parameter that Perrin et al. (2018) reported on was pixel spacing, which is closely related to slice thickness, given that both variables directly affect the resulting voxel volume of the 3D image. They found that scan pairs with a greater difference in pixel spacing had lower feature agreement for features extracted from both tumor and liver parenchyma ROIs, highlighting along with our own results that the feature reproducibility across variations of any parameters affecting voxel volume, including slice thickness and field of view, is an important consideration when drawing inference from radiomic features. For all the variables they considered, including pixel spacing, Perrin et al. (2018) also found that the tumor tended to have a larger number of reproducible features than the liver parenchyma. Similarly, our reproducibility analysis showed that features extracted from the liver parenchyma were overall less reproducible across slice thickness variation than those drawn from the largest metastasis, regardless of the extractor setting under consideration.

Despite the lower reproducibility of the features drawn from the liver parenchyma, our univariate survival analysis showed that the liver features still contained a relevant signal with regard to overall patient survival in the independent retrospective data set. Indeed, the highest overall univariate C-index was attained by a feature from the liver parenchyma—GLDM SmallDependenceLowGrayLevelEmphasis (C-index=0.6176)—and features from the liver parenchyma were well-represented on the Pareto front for C-index and CCC (see Table 5). Similarly, in the top-performing multivariable models, across all cross-validation runs, the average proportion of selected features coming from the liver was usually over 50% (see Table 4). The presence of important information relevant to patient outcomes in the radiological texture of the liver parenchyma is biologically feasible, and consistent with previous analyses in the literature (Simpson et al., 2017; Hu et al., 2022b). Given the apparent importance of the liver parenchyma features, future studies should investigate methods of tailoring the feature extraction to improve the reproducibility of the features extracted from the liver parenchyma, without degrading their prognostic value.

The IBSI reference manual (Zwanenburg et al., 2016) recommends resampling to a common voxel size to ensure the features are comparable between images with different resolutions. However, several empirical studies have shown that resampling has a limited ability to mitigate the effects of slice thickness, or resolution generally, on radiomic features. Ligero et al. (2021) showed in a phantom study that resampling could improve feature concordance across varying slice thickness to a limited extent, but found that harmonization of features treating slice thickness as a batch effect had a stronger effect. In lung CT, studies have found limited improvement when resampling images to the same voxel volume with linear interpolation (Yang et al., 2021), although deep-learning-based super-resolution algorithms have shown promise (Park et al., 2019). In a study varying pixel spacing by varying field of view during reconstruction, Mackin et al. (2017) found that resampling on its own actually worsened reproducibility of features across pixel sizes; however, combining resampling with a low-pass Butterworth filter improved feature concordance. Shafiq-Ul-Hassan et al. (2017) found in a phantom study that resampling improved feature concordance across images acquired with different voxel sizes for only a small subset of features. The majority of these features could also be corrected without resampling, by including a correction factor based on the voxel volume and/or ROI volume—a result which they subsequently validated in a lung cancer cohort (Shafiq-ul Hassan et al., 2018). Using a similar approach to feature correction, Escudero Sanchez et al. (2021), in a study of contrast enhanced CT imaging of liver cancer, found that while resampling to the mean spacing in all directions improved the number of features that were reproducible across changes in slice thickness, after feature correction, the result was modest in terms of absolute improvement over extracting features from the original images. In aggregate, these studies show that isotropic resampling is, on its own, not sufficient to mitigate the effects of slice thickness variation on the resulting radiomic feature values.

Our results also indicate limitations of resampling to mitigate voxel resolution effects. In the present study, we were further motivated to question the use of isotropic resampling—as recommended by IBSI (Zwanenburg et al., 2016) to preserve rotational invariance of feature definitions—due to the inherently anisotropic nature of our abdominal images. In our previous study on this topic (Peoples et al., 2023), we showed that aggregating features in 2.5D with no z-axis resampling (i.e. L2i and S2i) tended to produce more reproducible texture features, and still produced good prognostic models. In the present work, we conducted a more in-depth analysis of the relationship between feature extraction method, reproducibility, and prognostic value. From Figure 4, we can see that no one extraction method produces more reproducible or more prognostic features in every case. Indeed, both Figure 4 and Figure 5 illustrate that there is no consistently “best” feature extraction approach across all features and ROI, when we consider either C-index or CCC in isolation, or when we consider both by looking at the Pareto front. Furthermore, while the best overall performance in the multivariable survival models was obtained by using features from L2i, the combination of all extractors was able to rival this performance when we removed features with $\mathrm{CCC}<0.85$ across slice thicknesses. These results suggest that when feature reproducibility across important variations in imaging protocol is known a priori (such as from a test-retest study), including reproducibility into the feature selection process can improve model performance. In this case, highly reproducible and prognostic models can be achieved without optimizing the feature extraction process, by including features extracted using a variety of settings, and allowing the best features to be selected in a data-driven manner accounting for both reproducibility and discriminative ability. This data-driven approach is reminiscent of a method used by Vallières et al. (2017), wherein features extracted with multiple settings, such as resampling resolutions, were pooled into a larger table, and allowed to be considered as separate features during model building. In the absence of reproducibility scores for each feature, choosing settings that can produce features that are more robust to protocol variations present in the data set may be more important, but it is unclear how to choose these settings without a reproducibility study, given the study-specific nature of feature CCCs (van Timmeren et al., 2016).

One limitation of this study is the use of only one fixed bin count of 24 when computing the texture features. The bin count is a key parameter for computing texture features, and has a large effect on the results; though, for many features this effect is at least somewhat predictable (Shafiq-Ul-Hassan et al., 2017; Shafiq-ul Hassan et al., 2018). Furthermore, phantom studies have suggested that variation of the discretization, although it affects the feature values, does not have a large effect on feature reproducibility (Larue et al., 2017). Due to the large number of modeling experiments, and already large number of dimensions under consideration in the present study, we chose to avoid adding bin count as an additional variable at this time. In a study of contrast enhanced CT of hepatocellular carcinoma, Escudero Sanchez et al. (2021) found that the optimal bin count for feature reproducibility across slice thicknesses was in the range of (32–64), which given the step-size in their analysis, is nearly encompassing our chosen value of 24. Depite this, future work should consider the effect of bin count on feature reproducibility and prognostic value.

In conclusion, our results demonstrate the strong effect of slice thickness on feature reproducibility for contrast enhanced CT imaging of CRLM. Although some methods of feature extraction may mitigate the effects of slice thickness on some features, overall, we found that the extractor producing the most reproducible value for a given feature can vary across features and ROI. Similarly, the greatest discriminative ability for a given feature may be attained by different extractors, dependent on feature, and ROI. Given this, our results support a data-driven approach, where features from a variety of extractor settings are all considered, and selected in a manner accounting for reproducibility across relevant variations in protocol. Where disease-specific reproducibility metrics are not available, some methods of feature extraction may perform better due to improvements in reproducibility across certain prognostic features (such as the L2i features in the present study), however it is unclear how to determine this without a reproducibility study. Overall, our results demonstrate that we can find radiomic features that are both reproducible across slice thickness variation, and prognostic in patients undergoing hepatic resection, in the context of contrast enhanced CT of colorectal liver metastases.

Table 5 lists all features on the Pareto front for CCC and C-index when combining all feature extractor settings and ROI.

Table LABEL:tab:feature_list lists all features used in this study.