# Evaluating the Impact of Sequence Combinations on Breast Tumor Segmentation
in Multiparametric MRI

###### Abstract

Multiparametric magnetic resonance imaging (mpMRI) is a key tool for assessing
breast cancer progression. Although deep learning has been applied to automate
tumor segmentation in breast MRI, the effect of sequence combinations in mpMRI
remains under-investigated. This study explores the impact of different
combinations of T2-weighted (T2w), dynamic contrast-enhanced MRI (DCE-MRI) and
diffusion-weighted imaging (DWI) with apparent diffusion coefficient (ADC) map
on breast tumor segmentation using nnU-Net. Evaluated on a multicenter mpMRI
dataset, the nnU-Net model using DCE sequences achieved a Dice similarity
coefficient (DSC) of 0.69 $\pm$ 0.18 for functional tumor volume (FTV)
segmentation. For whole tumor mask (WTM) segmentation, adding the predicted
FTV to DWI and ADC map improved the DSC from 0.57 $\pm$ 0.24 to 0.60 $\pm$
0.21. Adding T2w did not yield significant improvement, which still requires
further investigation under a more standardized imaging protocol. This study
serves as a foundation for future work on predicting breast cancer treatment
response using mpMRI.

Index Terms—  Breast cancer, multiparametric MRI, tumor segmentation, nnU-Net

##  1 Introduction

Breast cancer is one of the most prevalent cancers among women worldwide [1].
Multiparametric magnetic resonance imaging (mpMRI) is a valuable tool in
breast imaging that facilitates not only the diagnosis of breast cancer but
also the assessment of disease progression and treatment response [2].
Commonly incorporating T2-weighted (T2w), dynamic contrast-enhanced MRI (DCE-
MRI) and diffusion-weighted imaging (DWI), mpMRI offers both anatomical and
functional information of tumors. Specifically, T2w imaging delivers
anatomical and fluid content information, while DCE and DWI quantify tumor
perfusion and diffusion patterns respectively [3, 4]. Tumor segmentation is an
important step in the quantitative analysis of breast MRI data and is commonly
carried out manually by radiologists [5]. However, manual segmentation of
mpMRI data is both challenging and labor-intensive, given that mpMRI generates
a large volume of data across multiple image sequences and temporal
resolutions.

In recent years, deep learning (DL) techniques have gained significant
attention in the field of medical image analysis. They are increasingly
applied to automate breast cancer segmentation on breast MRI. Chen et al. [6]
introduced a U-Net [7] with ConvLSTM structure [8] to analyze both spatial and
temporal features in DCE-MRI for breast tumor segmentation. However, this
approach only operates at the slice level rather than the entire 4D data
(which involve 3D volumetric DCE images with a temporal dimension). Both Xu et
al. [5] and Janse et al. [9] utilized nnU-Net [10], a self-adaptive network,
to segment breast tumors on DCE-MRI. Despite the prevalence of DL applications
in DCE-MRI, there is a noticeable gap in research focusing on tumor
segmentation on mpMRI. Zhu et al. [11] proposed applying a 2D VNet [12] on
DCE-MRI and attention-UNet [13] on DWI for breast tumor segmentation. However,
this work also only operates on 2D slice-level instead of the multimodal 3D
image data with the temporal dimension generated by mpMRI. Currently, there
remains a lack of studies which investigate how to optimally combine mpMRI
sequences to fully harness the rich information they offer for enhancing
breast tumor segmentation.

This study aims to investigate the impact of different combinations of
anatomical and functional sequences in mpMRI on the performance of breast
tumor segmentation. The evaluation was conducted on a publicly available,
multicenter mpMRI dataset, using different combinations of T2w, DCE, and DWI
sequences with apparent diffusion coefficient (ADC) map as inputs for nnU-Net,
U-Net with a ConvLSTM structure, and 3D U-Net. The segmentation outcomes were
quantitatively analyzed to identify the most effective sequence combinations.

##  2 Methodology

###  2.1 Datasets

The Breast Multiparametric MRI for prediction of neoadjuvant chemotherapy
(NAC) Response Challenge (BMMR2 Challenge) dataset [14, 15], which is a
multicenter breast mpMRI dataset, was used in this study. It contains 191
subjects divided into approximately 60% for training and 40% for testing. The
dataset contains both pre- and post-treatment mpMRI data and only pre-
treatment data were used in this study to avoid possible alterations to the
lesions induced by NAC. Each subject contains T1-weighted (T1w) DCE sequences
($80s<\text{phase duration}<100s$ with at least 8 minutes continuous post-
injection acquisition), 4 b-value (b=0, 100, 600, 800$s/mm^{2}$) DWI with ADC
map and T2w imaging that varies in acquisition sequences and fat suppression
methods. One case from the training set and one from the testing set were
excluded since not all sequences described in the imaging protocol can be
identified. The dataset also includes single-breast ipsilateral cropped DCE
images that were used to generate the functional tumor volume (FTV) mask.
After the FTV was identified on post-contrast DCE subtraction images, a whole
tumor mask (WTM) was subsequently segmented on the ADC map manually.

###  2.2 Data preparation

To prepare the mpMRI for the DL networks, the T2w, DWI and ADC images were
resampled to align with the single breast cropped DCE images. The DCE images
were temporally resampled into seven DCE frames. The interval between the pre-
contrast and the first post-contrast image was set as 150s, while a 90s
interval was maintained among the subsequent post-contrast images. Six DCE
subtraction images DCE${}_{sub}=\\{DCE_{i}-DCE_{0}|i=1,2,...6\\}$ were then
generated for the DL segmentation networks, where $DCE_{0}$ stands for the
pre-contrast image and $DCE_{i}$ for the post-contrast ones.

###  2.3 Segmentation model training and evaluation

The nnU-Netv2, a recent release of the nnU-Net [10], was trained to segment
FTV and WTM separately. The nnU-Net is a self-configuring approach for
biomedical image segmentation that includes preprocessing, network
architecture, training, and post-processing. It has been applied to a wide
range of biomedical image segmentation tasks and demonstrated outstanding
performance [16]. In this study, the 3D full resolution nnU-Net with
multichannel input was used.

To segment the FTV and WTM, various mpMRI sequence combinations were explored
as listed in Table 1. For FTV segmentation, the initial combination started
with using all DCEsub images, and then adding T2w. For WTM segmentation, the
ADC map, which was used as the primary reference during WTM manual
segmentation, was the initial input explored for nnU-Net. Other inputs
including DWI, predicted FTV (FTVpred1 and FTVpred2 as shown in Table 1 row 1
and 2) and T2w were subsequently introduced. The nnU-Net models were trained
via 5-fold cross validation within the training set for 100 epochs with a
batch number of 2 on two NVidia Tesla P100 with 16GB memory.

A 3D U-Net with ConvLSTM structure (ConvLSTM-UNet) and a 3D U-Net with
multichannel inputs were used for FTV segmentation on DCEsub images and WTM
segmentation on the DWI sequence with ADC map respectively for comparison. The
ConvLSTM-UNet, which is an extension of Chen et al. [6], attaches a ConvLSTM
block with 16 output features to the first decoding layer of a 3D U-Net to
enable the network to process 3D DCEsub images across 6 temporal frames. The
3D U-Net consists of an encoding and a decoding path, each containing three
layers, and a bridge layer in between. The numbers of features for these
layers are 16, 32, 64 and 128 respectively. Both networks were trained via
5-fold cross validation within the training set for 100 epochs with a batch
number of 2.

Table 1: MpMRI sequence combinations. FTV stands for functional tumor volume and WTM for whole tumor mask. | Sequence combinations | Prediction output
---|---|---
1 | DCEsub | FTVpred1
2 | DCEsub, T2w | FTVpred2
3 | ADC | WTMpred1
4 | DWI, ADC | WTMpred2
5 | DWI, ADC, T2w | WTMpred3
6 | DWI, ADC, FTVpred1 | WTMpred4
7 | DWI, ADC, FTVpred2 | WTMpred5
8 | DWI, ADC, FTVpred1, T2w | WTMpred6
9 | DWI, ADC, FTVpred2, T2w | WTMpred7

After model training, the best configuration plan was automatically determined
by nnU-Net on the validation sets, which was used to generate FTV and WTM
predictions on the test set. Several metrics on the voxel level including
precision, recall, Dice similarity coefficient (DSC) and 95% Hausdorff
distance (HD95) were used to evaluate the FTV and WTM segmentation
performance. For ConvLSTM-UNet and U-Net, the segmentation results were
obtained by averaging the prediction outcomes from all models.

##  3 Results

Table 2: FTV and WTM segmentation performance using nnU-Net, ConvLSTM-UNet, U-Net with difference mpMRI sequence combinations. | Targets |  Networks | Sequence combinations | Metrics |  |  |
---|---|---|---|---|---|---|---
|  |  |  | Precision | Recall | DSC | $HD_{95}^{\mathrm{a}}$
1 | FTV |  ConvLSTM-UNet | DCEsub | 0.66 $\pm$ 0.19 | 0.67 $\pm$ 0.21 | 0.62 $\pm$ 0.16 | 15.49 $\pm$ 14.08
2 | FTV |  nnU-Net | DCEsub | 0.72 $\pm$ 0.21 | 0.75 $\pm$ 0.21 | 0.69 $\pm$ 0.18 | 14.83 $\pm$ 15.34
3 | FTV |  nnU-Net | DCEsub, T2w | 0.72 $\pm$ 0.20 | 0.75 $\pm$ 0.21 | 0.69 $\pm$ 0.18 | 13.67 $\pm$ 13.16
4 | WTM |  U-Net | DWI, ADC | 0.57 $\pm$ 0.26 | 0.61 $\pm$ 0.30 | 0.53 $\pm$ 0.24 | 9.93 $\pm$ 11.90
5 | WTM |  nnU-Net | ADC | 0.60 $\pm$ 0.29 | 0.59 $\pm$ 0.30 | 0.54 $\pm$ 0.26 | 10.57 $\pm$ 12.40
6 | WTM |  nnU-Net | DWI, ADC | 0.61 $\pm$ 0.26 | 0.63 $\pm$ 0.27 | 0.57 $\pm$ 0.24 | 7.70 $\pm$ 7.22
7 | WTM |  nnU-Net | DWI, ADC, T2w | 0.61 $\pm$ 0.27 | 0.60 $\pm$ 0.31 | 0.56 $\pm$ 0.25 | 8.44 $\pm$ 8.03
8 | WTM |  nnU-Net | DWI, ADC, FTVpred1 | 0.63 $\pm$ 0.25 | 0.68 $\pm$ 0.22 | 0.60 $\pm$ 0.21 | 6.75 $\pm$ 7.33
9 | WTM |  nnU-Net | DWI, ADC, FTVpred2 | 0.64 $\pm$ 0.24 | 0.67 $\pm$ 0.23 | 0.60 $\pm$ 0.21 | 7.06 $\pm$ 7.64
10 | WTM |  nnU-Net | DWI, ADC, FTVpred1, T2w | 0.64 $\pm$ 0.23 | 0.67 $\pm$ 0.23 | 0.61 $\pm$ 0.20 | 6.84 $\pm$ 6.41
11 | WTM |  nnU-Net | DWI, ADC, FTVpred2, T2w | 0.65 $\pm$ 0.22 | 0.66 $\pm$ 0.23 | 0.61 $\pm$ 0.20 | 7.42 $\pm$ 7.96
a$HD_{95}$ was calculated between the true positive regions (i.e. regions that overlap with the groundtruth) and the groundtruth annotation. |  |  |  |  |  |  |
Fig. 1: FTV and WTM segmentation examples using ConvLSTM-UNet, U-Net and nnU-
Net with different sequence combinations in axial, coronal and sagittal views
(from top to bottom). The segmentation examples are outlined on the first
post-contrast DCE sequence for FTV segmentation and DWI
(b-value=800$s/mm^{2}$) for WTM segmentation. The network segmentation is
represented in green and the ground truth in red.

Table 2 lists the performance of nnU-Net, ConvLSTM-UNet and 3D U-Net using
different mpMRI sequence combinations for FTV and WTM segmentation. For FTV
segmentation, nnU-Net significantly outperformed ConvLSTM-Unet using DCEsub
images (p-value $<$ 0.05 for precision, recall and DSC using Wilcoxon test) as
shown in Table 2 row 1 and 2. Compared with using only the DCEsub images with
nnU-Net, the addition of T2w sequence resulted in similar outcomes (p-value
$>$ 0.05 for all metrics) as reflected in Table 2 row 2 and 3.

For WTM segmentation, nnU-Net outperformed U-Net using DWI sequences and ADC
map (p-value $<$ 0.05 for precision and DSC) as shown in Table 2 row 4 and 6.
It can be observed in Table 2 row 5, 6 and 7 that incorporating DWI sequences
with ADC map yielded better performance than using ADC alone with nnU-Net
(p-value $<$ 0.05 for recall and HD95), while the addition of T2w sequence did
not result in improvement (p-value $>$ 0.05 for all metrics). As indicated in
Table 2 row 6, 8 and 9, integrating the predicted FTV mask with DWI and ADC
map improved the segmentation performance (p-value $<$ 0.05 for DSC and HD95),
while no significant difference was observed between incorporating FTVpred1
and FTVpred2 (p-value $>$ 0.05 for all metrics). The addition of T2w, as shown
in Table 2 row 10 and 11, did not yield a significant difference in
performance (p-value $>$ 0.05 for all metrics). The nnU-Net models took
approximately 2-3 minutes to segment each image, whereas the ConvLSTM-UNet and
3D U-Net took less than a minute.

Considering both the segmentation accuracy and model complexity, the optimal
approach for segmenting mpMRI appears to be as follows: The nnU-Net models
with DCEsub as inputs are firstly utilized for FTV segmentation (generating
FTVpred1), and the nnU-Net models with DWI+ADC+ FTVpred1 as inputs are then
employed for WTM segmentation. Fig. 1 illustrates the FTV and WTM segmentation
examples using ConvLSTM-UNet, U-Net, nnU-Net with several selected sequence
combinations.

##  4 Discussion

In this study, we assessed the effect of mpMRI input combinations on breast
tumor segmentation performance using nnU-Net, which demonstrated superior
results compared with ConvLSTM-UNet and 3D U-Net. Our experiments suggest that
the most effective practice is to employ DCEsub images for FTV segmentation
and a combination of DWI, ADC and the predicted FTV mask for WTM segmentation.

In clinical practice, radiologists commonly employ a combination of T2w, DCE,
and DWI along with ADC maps for precise tumor contouring. However, in many
studies exploring DL methods for breast tumor segmentation, either the DCE
sequence was utilized [5, 6, 9], or segmentation was conducted independently
on DCE and DWI sequences [11]. In this study, the FTV predicted from the DCE
sequence was introduced as an additional input channel, along with DWI and ADC
map, to guide the segmentation process. This approach significantly improved
the WTM segmentation outcome. The addition of T2w sequence did not result in a
significant improvement for FTV and WTM segmentation. This could be attributed
to the fact that the imaging protocol for the T2w sequence in the dataset was
inconsistent, leading to data variability. Consequently, the potential
benefits of incorporating a T2w sequence into the segmentation process cannot
be entirely dismissed.

The current study still has several limitations and presents opportunities for
further extension. This study only focuses on the pre-treatment mpMRI data and
only involved the cropped breast images with tumor provided in the dataset.
Segmenting tumors in post-treatment mpMRI poses additional challenges due to
potential changes in tumor size, shape, and texture. For future work, we aim
to extend our methodology to segment tumors in post-treatment mpMRI,
leveraging pre-treatment segmentation as attention guidance to improve the
accuracy of post-treatment segmentation. To determine the potential benefits
of including T2w sequences in tumor segmentation, further comparative studies
should be carried out using an mpMRI cohort with a more standardized imaging
protocol. This study also serves as a preliminary investigation into
predicting NAC treatment responses. Going forward, we aim to conduct
quantitative analysis of breast tumors using mpMRI to identify imaging
biomarkers that are predictive of NAC treatment outcomes.

##  5 Conclusion

In this study, we conducted a comprehensive analysis comparing the
segmentation outcomes of FTV and WTM on a multicenter breast mpMRI dataset
using different sequence combinations. For FTV segmentation, the nnU-Net using
DCE subtraction volumes achieved satisfactory results. While the inclusion of
the T2w sequence led to slight improvement in some evaluation metrics, these
are not statistically significant. For WTM segmentation, the addition of
predicted FTV to DWI sequences and ADC map as network inputs significantly
improved the segmentation performance. Although the inclusion of T2w did not
yield a statistically significant improvement, further experiments are still
needed to evaluate its potential for enhancing tumor segmentation performance,
particularly under a more consistent imaging protocol. This work serves as a
foundational study for future research aimed at predicting breast NAC
treatment responses.

##  6 COMPLIANCE WITH ETHICAL STANDARDS

This study used the publicly available BMMR2 dataset [14, 15] (with CC BY 4.0
license), which is a subset of the ACRIN trial 6698 (NCT01564368). The ethical
approval was acquired by the original data providers.

##  7 Acknowledgments

This study is supported by Metro South Health Study, Education and Research
Trust Account grant (grant number: RSS_2020_004), and Australian Government
Medical Research Future Fund (grant number: GHFM76722).

## References

  * [1] A. G. Waks and E. P. Winer, “Breast cancer treatment: a review,” _Jama_ , vol. 321, no. 3, pp. 288–300, 2019.
  * [2] M. A. Marino, T. Helbich, P. Baltzer, and K. Pinker-Domenig, “Multiparametric mri of the breast: A review,” _Journal of Magnetic Resonance Imaging_ , vol. 47, no. 2, pp. 301–315, 2018.
  * [3] K. A. Wahid, S. Ahmed, R. He, L. V. van Dijk, J. Teuwen, B. A. McDonald, V. Salama, A. S. Mohamed, T. Salzillo, C. Dede _et al._ , “Evaluation of deep learning-based multiparametric mri oropharyngeal primary tumor auto-segmentation and investigation of input channel effects: Results from a prospective imaging registry,” _Clinical and translational radiation oncology_ , vol. 32, pp. 6–14, 2022.
  * [4] R. E. Hendrick, _Breast MRI: fundamentals and technical aspects_. Springer, 2008, vol. 2.
  * [5] Z. Xu, D. E. Rauch, R. M. Mohamed, S. Pashapoor, Z. Zhou, B. Panthi, J. B. Son, K.-P. Hwang, B. C. Musall, B. E. Adrada _et al._ , “Deep learning for fully automatic tumor segmentation on serially acquired dynamic contrast-enhanced mri images of triple-negative breast cancer,” _Cancers_ , vol. 15, no. 19, p. 4829, 2023.
  * [6] M. Chen, H. Zheng, C. Lu, E. Tu, J. Yang, and N. Kasabov, “A spatio-temporal fully convolutional network for breast lesion segmentation in dce-mri,” in _Neural Information Processing: 25th International Conference, ICONIP 2018, Siem Reap, Cambodia, December 13–16, 2018, Proceedings, Part VII 25_. Springer, 2018, pp. 358–368.
  * [7] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18_. Springer, 2015, pp. 234–241.
  * [8] X. Shi, Z. Chen, H. Wang, D.-Y. Yeung, W.-K. Wong, and W.-c. Woo, “Convolutional lstm network: A machine learning approach for precipitation nowcasting,” _Advances in neural information processing systems_ , vol. 28, 2015.
  * [9] M. H. Janse, L. M. Janssen, B. H. van Der Velden, M. R. Moman, E. J. Wolters-van der Ben, M. C. Kock, M. A. Viergever, P. J. van Diest, and K. G. Gilhuijs, “Deep learning-based segmentation of locally advanced breast cancer on mri in relation to residual cancer burden: A multi-institutional cohort study,” _Journal of Magnetic Resonance Imaging_ , 2023.
  * [10] F. Isensee, P. F. Jaeger, S. A. Kohl, J. Petersen, and K. H. Maier-Hein, “nnu-net: a self-configuring method for deep learning-based biomedical image segmentation,” _Nature methods_ , vol. 18, no. 2, pp. 203–211, 2021.
  * [11] J. Zhu, J. Geng, W. Shan, B. Zhang, H. Shen, X. Dong, M. Liu, X. Li, and L. Cheng, “Development and validation of a deep learning model for breast lesion segmentation and characterization in multiparametric mri,” _Frontiers in oncology_ , vol. 12, p. 946580, 2022.
  * [12] F. Milletari, N. Navab, and S.-A. Ahmadi, “V-net: Fully convolutional neural networks for volumetric medical image segmentation,” in _2016 fourth international conference on 3D vision (3DV)_. Ieee, 2016, pp. 565–571.
  * [13] O. Oktay, J. Schlemper, L. L. Folgoc, M. Lee, M. Heinrich, K. Misawa, K. Mori, S. McDonagh, N. Y. Hammerla, B. Kainz _et al._ , “Attention u-net: Learning where to look for the pancreas,” _arXiv preprint arXiv:1804.03999_ , 2018.
  * [14] D. C. Newitt, S. C. Partridge, Z. Zhang, J. Gibbs, T. Chenevert, M. Rosen, P. Bolan, H. Marques, J. Romanoff, L. Cimino, B. N. Joe, H. Umphrey, H. Ojeda-Fournier, B. Dogan, K. Y. Oh, H. Abe, J. Drukteinis, L. J. Esserman, and N. M. Hylton, “ACRIN 6698/I-SPY2 Breast DWI [Data set],” The Cancer Imaging Archive https://doi.org/10.7937/TCIA.KK02-6D95, 2021.
  * [15] D. C. Newitt, Z. Zhang, J. E. Gibbs, S. C. Partridge, T. L. Chenevert, M. A. Rosen, P. J. Bolan, H. S. Marques, S. Aliu, W. Li _et al._ , “Test–retest repeatability and reproducibility of adc measures by breast dwi: Results from the acrin 6698 trial,” _Journal of Magnetic Resonance Imaging_ , vol. 49, no. 6, pp. 1617–1628, 2019.
  * [16] M. Alqaoud, J. Plemmons, E. Feliberti, S. Dong, K. Kaipa, G. Fichtinger, Y. Xiao, and M. A. Audette, “nnunet-based multi-modality breast mri segmentation and tissue-delineating phantom for robotic tumor surgery planning,” in _2022 44th Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC)_. IEEE, 2022, pp. 3495–3501.


# Recalibration of Aleatoric and Epistemic
Regression Uncertainty in Medical Imaging

Max-Heinrich Laves<EMAIL_ADDRESS>
Institute of Medical Technology and Intelligent Systems, Hamburg University of
Technology, Hamburg, Germany
Institute of Mechatronic Systems, Leibniz Universität Hannover, Hannover,
Germany Sontje Ihler<EMAIL_ADDRESS>
Institute of Mechatronic Systems, Leibniz Universität Hannover, Hannover,
Germany Jacob F. Fast<EMAIL_ADDRESS>
Institute of Mechatronic Systems, Leibniz Universität Hannover, Hannover,
Germany
Department of Phoniatrics and Pediatric Audiology, Hannover Medical School,
Hannover, Germany Lüder A. Kahrs<EMAIL_ADDRESS>
Centre for Image Guided Innovation and Therapeutic Intervention, The Hospital
for Sick Children, Toronto, Canada
Department of Mathematical and Computational Sciences, University of Toronto
Mississauga, Mississauga, Canada Tobias Ortmaier<EMAIL_ADDRESS>
Institute of Mechatronic Systems, Leibniz Universität Hannover, Hannover,
Germany

###### Abstract

The consideration of predictive uncertainty in medical imaging with deep
learning is of utmost importance. We apply estimation of both aleatoric and
epistemic uncertainty by variational Bayesian inference with Monte Carlo
dropout to regression tasks and show that predictive uncertainty is
systematically underestimated. We apply $\sigma$ scaling with a single scalar
value; a simple, yet effective calibration method for both types of
uncertainty. The performance of our approach is evaluated on a variety of
common medical regression data sets using different state-of-the-art
convolutional network architectures. In our experiments, $\sigma$ scaling is
able to reliably recalibrate predictive uncertainty. It is easy to implement
and maintains the accuracy. Well-calibrated uncertainty in regression allows
robust rejection of unreliable predictions or detection of out-of-distribution
samples. Our source code is available at: github.com/mlaves/well-calibrated-
regression-uncertainty

Keywords: Bayesian approximation, variational inference

##  1 Introduction

Predictive uncertainty should be considered in any medical imaging task that
is approached with deep learning. Well-calibrated uncertainty is of great
importance for decision-making and is anticipated to increase patient safety.
It allows to robustly reject unreliable predictions or out-of-distribution
samples. In this paper, we address the problem of miscalibration of regression
uncertainty with application to medical image analysis.

For the task of regression, we aim to estimate a continuous target value
$\bm{y}\in\mathbb{R}^{d}$ given an input image $\bm{x}$. Regression in medical
imaging with deep learning has been applied to forensic age estimation from
hand CT/MRI (Halabi et al., 2019; Štern et al., 2016), natural landmark
localization (Payer et al., 2019), cell detection in histology (Xie et al.,
2018), or instrument pose estimation (Gessert et al., 2018). By predicting the
coordinates of object boundaries, segmentation can also be performed as a
regression task. This has been done for segmentation of pulmonary nodules in
CT (Messay et al., 2015), kidneys in ultrasound (Yin et al., 2020), or left
ventricles in MRI (Tan et al., 2017). In registration of medical images, a
continuous displacement field is predicted for each coordinate of $\bm{x}$,
which has also recently been addressed by CNNs for regression (Dalca et al.,
2019).

In medical imaging, it is crucial to consider the predictive uncertainty of
deep learning models. Bayesian neural networks (BNN) and their approximation
provide mathematical tools for reasoning the uncertainty (Bishop, 2006; Kingma
and Welling, 2014). In general, predictive uncertainty can be split into two
types: aleatoric and epistemic uncertainty (Tanno et al., 2017; Kendall and
Gal, 2017). This distinction was first made in the context of risk management
(Hora, 1996). Aleatoric uncertainty arises from the data directly; e. g.
sensor noise or motion artifacts. Epistemic uncertainty is caused by
uncertainty in the model parameters due to a limited amount of training data
(Bishop, 2006). A well-accepted approach to quantify epistemic uncertainty is
variational inference with Monte Carlo (MC) dropout, where dropout is used at
test time to sample from the approximate posterior (Gal and Ghahramani, 2016).

Figure 1: Calibration plots and uncertainty calibration error (UCE) for
EfficientNet-B4 on BreastPathQ test set. Uncalibrated uncertainty is
underestimated and does not correspond well with the model error (left).
Uncertainty can be calibrated most effectively with $\sigma$ scaling (right).
Solid lines show the mean and shaded areas show standard deviation from 5
repeated runs. Dashed lines denote perfect calibration.

Uncertainty quantification in regression problems in medical imaging has been
addressed by prior work. Medical image enhancement with image quality transfer
(IQT) has been extended to a Bayesian approach to obtain pixel-wise
uncertainty (Tanno et al., 2016). Additionally, CNN-based IQT was used to
estimate both aleatoric and epistemic uncertainty in MRI super-resolution
(Tanno et al., 2017). Dalca et al. (2019) estimated uncertainty for a
deformation field in medical image registration using a probabilistic CNN.
Registration uncertainty has also been addressed outside the deep learning
community (Luo et al., 2019). Schlemper et al. (2018) used sub-network
ensembles to obtain uncertainty estimates in cardiac MRI reconstruction.
Aleatoric and epistemic uncertainty was also used in multitask learning for
MRI-based radiotherapy planning (Bragman et al., 2018).

Uncertainty obtained by deep BNNs tends to be miscalibrated, i. e. it does not
correlate well with the model error (Laves et al., 2019). Fig. 1 shows
calibration plots (observed uncertainty vs. expected uncertainty) for
uncalibrated and calibrated uncertainty estimates. The predicted uncertainty
(taking into account both epistemic and aleatoric uncertainty) is
underestimated and does not allow robust detection of uncertain predictions at
test time.

Calibration of uncertainty in regression has been addressed in prior work
outside medical imaging. In (Kuleshov et al., 2018), inaccurate uncertainties
from Bayesian models for regression are recalibrated using a technique
inspired by Platt scaling. Given a pre-trained, miscalibrated model $\bm{H}$,
an auxiliary model $\bm{R}:[0,1]^{d}\rightarrow[0,1]^{d}$ is trained, that
yields a calibrated regressor $\bm{R}\circ\bm{H}$. In (Phan et al., 2018),
this method was applied to bounding box regression. However, an auxiliary
model with enough capacity will always be able to recalibrate, even if the
predicted uncertainty is completely uncorrelated with the real uncertainty.
Furthermore, Kuleshov et al. (2018) state that calibration via $\bm{R}$ is
possible if enough independent and identically distributed (i.i.d.) data is
available. In medical imaging, large data sets are usually hard to obtain,
which can cause $\bm{R}$ to overfit the calibration set. This downside was
addressed in (Levi et al., 2019), which is most related to our work. They
proposed to scale the standard deviation of a Gaussian model to recalibrate
aleatoric uncertainty. In contrast to our work, they do not take into account
epistemic uncertainty, which is an important source of uncertainty, especially
when dealing with small data sets in medical imaging.

This paper extends a preliminary version of this work presented at the Medical
Imaging with Deep Learning (MIDL) 2020 conference (Laves et al., 2020). We
continue this work by providing a new derivation of our definition of perfect
calibrtaion, new experimental results, analysis and discussion. Additionally,
prediction intervals are computed to further assess the quality of the
estimated uncertainty. We find that prediction intervals are estimated too
narrow and that recalibration can mitigate this problem.

To the best of our knowledge, calibration of predictive uncertainty for
regression tasks in medical imaging has not been addressed. Our main
contributions are: (1) We suggest to use $\sigma$ scaling in a separate
calibration phase to tackle underestimation of aleatoric and epistemic
uncertainty, (2) we propose to use the uncertainty calibration error and
prediction intervals to assess the quality of the estimated uncertainty, and
(3) we perform extensive experiments on four different data sets to show the
effectiveness of the proposed method.

##  2 Methods

In this section, we discuss estimation of aleatoric and epistemic uncertainty
for regression and show why uncertainty is systematically miscalibrated. We
propose to use $\sigma$ scaling to jointly calibrate aleatoric and epistemic
uncertainty.

###  2.1 Conditional Log-Likelihood for Regression

We revisit regression under the maximum posterior (MAP) framework to derive
direct estimation of heteroscedastic aleatoric uncertainty. That is, the
aleatoric uncertainty varies with the input and is not assumed to be constant.
The goal of our regression model is to predict a target value $\bm{y}$ given
some new input $\bm{x}$ and a training set $\mathcal{D}$ of $m$ inputs
$\bm{X}=\\{\bm{x}_{1},\ldots,\bm{x}_{m}\\}$ and their corresponding (observed)
target values $\bm{Y}=\\{\bm{y}_{1},\ldots,\bm{y}_{m}\\}$. We assume that
$\bm{y}$ has a Gaussian distribution
$\mathcal{N}\left(\bm{y};\hat{\bm{y}}(\bm{x}),\hat{\sigma}^{2}(\bm{x})\right)$
with mean equal to $\hat{\bm{y}}(\bm{x})$ and variance
$\hat{\sigma}^{2}(\bm{x})$. A neural network with parameters $\bm{\theta}$

$\bm{f}_{\bm{\theta}}\left(\bm{x}\right)=\left[\hat{\bm{y}}(\bm{x}),\hat{\sigma}^{2}(\bm{x})\right],~{}\hat{\bm{y}}\in\mathbb{R}^{d},~{}\hat{\sigma}^{2}\in\mathbb{R},\hat{\sigma}^{2}\geq
0$ (1)

outputs these values for a given input (Nix and Weigend, 1994). By assuming a
Gaussian prior over the parameters
$\bm{\theta}\sim\mathcal{N}(\bm{\theta};\bm{0},\lambda^{-1}\bm{I})$, MAP
estimation becomes maximum-likelihood estimation with added weight decay
(Bishop, 2006). With $m$ i.i.d. random samples, the conditional log-likelihood
$\log p(\bm{Y}\,|\,\bm{X},\bm{\theta})$ is given by

$\displaystyle\log p(\bm{Y}\,|\,\bm{X},\bm{\theta})=$
$\displaystyle\sum_{i=1}^{m}\log\left(\frac{1}{\sqrt{2\pi}\hat{\sigma}^{(i)}_{\bm{\theta}}}\exp\left\\{-\frac{\big{\|}\bm{y}^{(i)}-\hat{\bm{y}}_{\bm{\theta}}^{(i)}\big{\|}^{2}}{2\big{(}\hat{\sigma}^{(i)}_{\bm{\theta}}\big{)}^{2}}\right\\}\right)$
(2) $\displaystyle=$
$\displaystyle-\dfrac{m}{2}\log\left(2\pi\right)-\sum_{i=1}^{m}\log\big{(}\hat{\sigma}^{(i)}_{\bm{\theta}}\big{)}+\frac{1}{2\big{(}\hat{\sigma}_{\bm{\theta}}^{(i)}\big{)}^{2}}\big{\|}\bm{y}^{(i)}-\hat{\bm{y}}_{\bm{\theta}}^{(i)}\big{\|}^{2}~{}.$
(3)

The dependence on $\bm{x}$ has been omitted to simplify the notation.
Maximizing the log-likelihood in Eq. (3) w.r.t. $\bm{\theta}$ is equivalent to
minimizing the negative log-likelihood (NLL), which leads to the following
optimization criterion (with weight decay)

$\mathcal{L}_{\mathrm{G}}(\bm{\theta})=\sum_{i=1}^{m}\big{(}\hat{\sigma}^{(i)}_{\bm{\theta}}\big{)}^{-2}\big{\|}\bm{y}^{(i)}-\hat{\bm{y}}_{\bm{\theta}}^{(i)}\big{\|}^{2}+\log\big{(}(\hat{\sigma}_{\bm{\theta}}^{(i)})^{2}\big{)}~{}.$
(4)

Here, $\hat{\bm{y}}_{\bm{\theta}}$ and $\hat{\sigma}_{\bm{\theta}}$ are
estimated jointly by finding $\bm{\theta}$ that minimizes Eq. (4). This can be
achieved using gradient descent in a standard training procedure. In this
case, $\hat{\sigma}_{\bm{\theta}}$ captures the uncertainty that is inherent
in the data (aleatoric uncertainty). To avoid numerical instability due to
potential division by zero, we directly estimate
$\log\hat{\sigma}^{2}(\bm{x})$ and implement Eq. (4) in similar practice to
Kendall and Gal (2017).

###  2.2 Biased estimation of $\sigma$

Figure 2: Biased estimation of aleatoric uncertainty $\sigma^{2}$. The deep
model overfits estimation of $\bm{y}$ on the training set. On unseen test
data, the MSE of predictive mean is higher and $\sigma^{2}$ is underestimated.
Early stopping (e.g. at epoch 50) would result in an unbiased estimator, but
this would not be optimal in terms of test MSE.

Ignoring their dependence through $\bm{\theta}$, the solution to Eq. (4)
decouples estimation of $\hat{\bm{y}}$ and $\hat{\sigma}$. In case of a
Gaussian likelihood, minimizing Eq. (4) w.r.t. $\hat{\bm{y}}^{(i)}$ yields

$\hat{\bm{y}}^{(i)}=\operatorname*{arg\,min}_{\hat{\bm{y}}^{(i)}}\mathcal{L}_{\mathrm{G}}=\bm{y}^{(i)}~{}\operatorname{\forall}i~{}.$
(5)

Minimizing (4) w.r.t. $(\hat{\sigma}^{(i)})^{2}$ yields

$\big{(}\hat{\sigma}^{(i)}\big{)}^{2}=\operatorname*{arg\,min}_{(\hat{\sigma}^{(i)})^{2}}\mathcal{L}_{\mathrm{G}}=\|\bm{y}^{(i)}-\hat{\bm{y}}^{(i)}\|^{2}~{}\operatorname{\forall}i~{}.$
(6)

That is, estimation of $\sigma^{2}$ should perfectly reflect the squared
error. However, in Eq. (6) $\sigma^{2}$ is estimated relative to the estimated
mean $\hat{\bm{y}}$ and therefore biased. In fact, the maximum likelihood
solution systematically underestimates $\sigma^{2}$, which is a phenomenon of
overfitting the training set (Bishop, 2006). The squared error
$\|\bm{y}-\hat{\bm{y}}\|^{2}$ will be lower on the training set and
$\hat{\sigma}^{2}$ on new samples will be systematically too low (see Fig. 2).
This is a problem especially in deep learning, where large models have
millions of parameters and tend to overfit. To solve this issue, we introduce
a simple learnable scalar parameter $s$ to rescale the biased estimation of
$\sigma^{2}$.

###  2.3 $\sigma$ Scaling for Aleatoric Uncertainty

We first derive $\sigma$ scaling for aleatoric uncertainty. Using a Gaussian
model, we scale the standard deviation $\sigma$ with a scalar value $s$ to
recalibrate the probability density function

$p\left(\bm{y}|\bm{x};\hat{\bm{y}}(\bm{x}),\hat{\sigma}^{2}(\bm{x})\right)=\mathcal{N}\left(\bm{y};\hat{\bm{y}}(\bm{x}),(s\cdot\hat{\sigma}(\bm{x}))^{2}\right)~{}.$
(7)

This results in the following minimization objective:

$\mathcal{L}_{\mathrm{G}}(s)=m\log(s)+\tfrac{1}{2}s^{-2}\sum_{i=1}^{m}\big{(}\hat{\sigma}_{\bm{\theta}}^{(i)}\big{)}^{-2}\big{\|}\bm{y}^{(i)}-\hat{\bm{y}}_{\bm{\theta}}^{(i)}\big{\|}^{2}~{}.$
(8)

Eq. (8) is optimized w.r.t. $s$ with fixed $\bm{\theta}$ using gradient
descent in a separate calibration phase after training to calibrate aleatoric
uncertainty measured by $\hat{\sigma}_{\bm{\theta}}^{2}$. In case of a single
scalar, the solution to Eq. (8) can also be written in closed form as

$s=\pm\sqrt{\frac{1}{m}\sum_{i=1}^{m}\big{(}\hat{\sigma}_{\bm{\theta}}^{(i)}\big{)}^{-2}\big{\|}\bm{y}^{(i)}-\hat{\bm{y}}_{\bm{\theta}}^{(i)}\big{\|}^{2}}~{}.$
(9)

We apply $\sigma$ scaling to jointly calibrate aleatoric and epistemic
uncertainty in the next section.

###  2.4 Well-Calibrated Estimation of Predictive Uncertainty

So far we have assumed a MAP point estimate for $\bm{\theta}$ which does not
consider uncertainty in the parameters. To quantify both aleatoric and
epistemic uncertainty, we extend $\bm{f}_{\bm{\theta}}$ into a fully Bayesian
model under the variational inference framework with Monte Carlo dropout (Gal
and Ghahramani, 2016). In MC dropout, the model $\bm{f}_{\tilde{\bm{\theta}}}$
is trained with dropout (Srivastava et al., 2014) and dropout is applied at
test time by performing $N$ stochastic forward passes to sample from the
approximate Bayesian posterior $\tilde{\bm{\theta}}\sim q(\bm{\theta})$.
Following (Kendall and Gal, 2017), we use MC integration to approximate the
predictive variance

$\hat{\Sigma}^{2}=\underbrace{\frac{1}{N}\sum_{n=1}^{N}\left(\hat{\bm{y}}_{n}-\frac{1}{N}\sum_{n=1}^{N}\hat{\bm{y}}_{n}\right)^{2}}_{\mathrm{epistemic}}+\underbrace{\frac{1}{N}\sum_{n=1}^{N}\hat{\sigma}^{2}_{n}}_{\mathrm{aleatoric}}$
(10)

and use $\hat{\Sigma}^{2}$ as a measure of predictive uncertainty. If the
neural network has multiple outputs ($d>1$), the predictive variance is
calculated per output and the mean across $d$ forms the final uncertainty
value. Eq. (10) is an unbiased estimator of the approximate predictive
variance (see proof in Appendix C). From Eq. (33) of our proof follows, that
$\hat{\Sigma}^{2}$ is expected to equal the true variance
$\Sigma=\mathbb{E}[(\hat{\bm{y}}-\bm{y})^{2}]$. Thus, we define perfect
calibration of regression uncertainty as

$\mathbb{E}_{\bm{x},\bm{y}}\left[\mathbb{E}[(\hat{\bm{y}}-\bm{y})^{2}]\,\big{|}\,\hat{\Sigma}^{2}=\alpha^{2}\right]=\alpha^{2}\quad\operatorname{\forall}\left\\{\alpha^{2}\in\mathbb{R}\,|\,\alpha^{2}\geq
0\right\\}~{},$ (11)

which extends the definition of (Levi et al., 2019) to both aleatoric and
epistemic uncertainty. We expect that additionally accounting for epistemic
uncertainty is particularly beneficial for smaller data sets. However, even in
deep learning with Bayesian principles, the approximate posterior predictive
distribution can overfit on small data sets. In practice, this leads to
underestimation of the predictive uncertainty.

One could regularize overfitting by early stopping that prevents large
differences between training and test loss, which would circumvent
underestimation of $\sigma^{2}$. However, our experiments show that early
stopping is not optimal with regard to accuracy, i.e. the squared error of
$\hat{\bm{y}}$ on both training and testing data (see Fig. 2). In contrast,
the model with lowest mean error on the validation set underestimates
predictive uncertainty considerably. Therefore, we apply $\sigma$ scaling to
recalibrate the predictive uncertainty $\hat{\Sigma}^{2}$. This allows a lower
squared error while reducing underestimation of uncertainty as shown
experimentally in the following section.

###  2.5 Expected Uncertainty Calibration Error for Regression

We extend the definition of miscalibrated uncertainty for classification
(Laves et al., 2019) to quantify miscalibration of uncertainty in regression

$\mathbb{E}_{\hat{\Sigma}^{2}}\left[\big{|}\big{(}\mathbb{E}[(\hat{\bm{y}}-\bm{y})^{2}]\,\big{|}\,\hat{\Sigma}^{2}=\alpha^{2}\big{)}-\alpha^{2}\big{|}\right]\quad\operatorname{\forall}\left\\{\alpha^{2}\in\mathbb{R}\,|\,\alpha^{2}\geq
0\right\\}~{},$ (12)

using the second moment of the error. On finite data sets, this can be
approximated with the expected uncertainty calibration error (UCE) for
regression. Following (Guo et al., 2017), the uncertainty output
$\hat{\Sigma}^{2}$ of a deep model is partitioned into $K$ bins with equal
width. A weighted average of the difference between the variance and
predictive uncertainty is used:

$\mathrm{UCE}:={\sum_{k=1}^{K}}\frac{|B_{k}|}{m}\big{|}{\mathrm{var}}(B_{k})-\mathrm{uncert}(B_{k})\big{|}~{},$
(13)

with number of inputs $m$ and set of indices $B_{k}$ of inputs, for which the
uncertainty falls into the bin $k$. The variance per bin is defined as

$\mathrm{var}(B_{k}):=\frac{1}{|B_{k}|}\sum_{i\in
B_{m}}{\frac{1}{N}\sum_{n=1}^{N}\left(\hat{\bm{y}}_{i,n}-\bm{y}_{i}\right)^{2}}~{},$
(14)

with $N$ stochastic forward passes, and the uncertainty per bin is defined as

$\mathrm{uncert}(B_{k}):=\frac{1}{|B_{k}|}\sum_{i\in
B_{k}}\hat{\Sigma}_{i}^{2}~{}.$ (15)

Note that computing the second moment from Eq. (12) also incorporates MC
samples, which can introduce some bias in the evaluation. The UCE considers
both aleatoric and epistemic uncertainty and is given in % throughout this
work. Additionally, we plot $\mathrm{var}(B_{k})$ vs. $\mathrm{uncert}(B_{k})$
to create calibration diagrams.

##  3 Experiments

We use four data sets and three common deep network architectures to evaluate
recalibration with $\sigma$ scaling. The data sets were selected to represent
various regression tasks in medical imaging with different dimension $d$ of
target value $\bm{y}\in\mathbb{R}^{d}$:

(1) Estimation of tumor cellularity in histology whole slides of cancerous
breast tissue from the BreastPathQ SPIE challenge data set ($d=1$) (Martel et
al., 2019). The public data set consists of 2579 images, from which
1379/600/600 are used for training/validation/testing. The ground truth label
is a single scalar $y\in[0,1]$ denoting the ratio of tumor cells to non-tumor
cells.

(2) Hand CT age regression from the RSNA pediatric bone age data set ($d=1$)
(Halabi et al., 2019). The task is to infer a person’s age in months from CT
scans of the hand. This data set is the largest used in this paper and has
12,811 images, from which we use 6811/2000/4000 images for
training/validation/testing.

(3) Surgical instrument tracking on endoscopic images from the EndoVis
endoscopic vision challenge 2015111endovissub-instrument.grand-challenge.org
data set ($d=2$). This data set contains 8,984 video frames from 6 different
robot-assisted laparoscopic interventions showing surgical instruments with
ground truth pixel coordinates of the instrument’s center point
$\bm{y}\in\mathbb{R}^{2}$. We use 4483/2248/2253 frames for
training/validation/testing. As the public data set is only sparsely
annotated, we created our own ground truth labels, which can be found in our
code repository.

(4) 6DoF needle pose estimation on optical coherence tomography (OCT) scans
from our own data set222Our OCT pose estimation data set is publicly available
at github.com/mlaves/3doct-pose-dataset. This data set contains 5,000 3D-OCT
scans with the accompanying needle pose $\bm{y}\in\mathbb{R}^{6}$, from which
we use 3300/850/850 for training/validation/testing. Additional details on
this data set can be found in Appendix E.

All outputs are normalized such that $\bm{y}\in[0,1]^{d}$. The employed
network architectures are ResNet-101, DenseNet-201 and EfficientNet-B4 (He et
al., 2016; Huang et al., 2017; Tan and Le, 2019). The last linear layer of all
networks is replaced by two linear layers predicting $\hat{\bm{y}}$ and
$\hat{\sigma}^{2}$ as described in § 2.1. For MC dropout, we use dropout
before the last linear layers. Dropout is further added after each of the four
layers of stacked residual blocks in ResNet. In DenseNet and EfficientNet, we
use the default configuration of dropout during training and testing. The
networks are trained until no further decrease in mean squared error (MSE) on
the validation set can be observed. More details on the training procedure can
be found in Appendix D.

Calibration is performed after training in a separate calibration phase using
the validation data set. We plug the predictive uncertainty $\hat{\Sigma}^{2}$
into Eq. (8) and minimize w.r.t. $s$. Additionally, we compare $\sigma$
scaling to a more powerful auxiliary recalibration model $\bm{R}$ consisting
of a two-layer fully-connected network with 16 hidden units and ReLU
activations (inspired by (Kuleshov et al., 2018), see § 1).

Figure 3: Calibration plots for ResNet-101 on BreastPathQ (top row) and
EfficientNet-B4 on EndoVis (bottom row). Aux scaling tends to overfit the
calibration set, which results in higher UCE compared to simple $\sigma$
scaling. Dashed lines denote perfect calibration.

##  4 Results

To quantify miscalibration, we use the proposed expected uncertainty
calibration error for regression. We visualize (mis-)calibration in Fig. 1 and
Fig. 3 using calibration diagrams, which show expected uncertainty vs.
observed uncertainty. The discrepancy to the identity function reveals
miscalibration. The calibration diagrams clearly show the underestimation of
uncertainty for the uncalibrated models. After calibration with both aux and
$\sigma$ scaling, the estimated uncertainty better reflects the actual
uncertainty. Figures for all configurations are listed in Appendix G.

Tab. 1 reports UCE values of all data set/model combinations on the respective
test sets. The negative log-likelihood also measures miscalibration; the
values on the test set can be found in Tab. 2 in the appendix. In general,
recalibration considerably reduces miscalibration. On the data sets BoneAge,
EndoVis and OCT, both scaling methods perform similarly well. However, on the
BreastPathQ data set, $\sigma$ scaling clearly outperforms aux scaling in
terms of UCE. BreastPathQ is the smallest data set and thus has the smallest
calibration set size. We hypothesize that the more powerful auxiliary model
$\bm{R}$ overfits the calibration set (see BreastPathQ/DenseNet-201 in Tab.
1), which leads to an increase of UCE on the test set. An ablation study on
BreastPathQ for the auxiliary model can be found in Appendix F.

We also compare our approach to Levi et al. (2019) in Tab. 1, which only
considers aleatoric uncertainty. The aleatoric uncertainty is well-calibrated
if it reflects the bias
$\left(\mathbb{E}\left[\hat{\bm{y}}_{n}\right]-\bm{y}\right)^{2}$, which is
given by the squared error between the expectation of the stochastic
predictions $\hat{\bm{y}}_{n}$ and the ground truth. Therefore, the UCE for
aleatoric-only is computed by
$\mathrm{UCE}={\sum_{k=1}^{K}}\frac{|B_{k}|}{m}\big{|}{\mathrm{err}}(B_{k})-\mathrm{uncert}(B_{k})\big{|}~{}$,
where $\mathrm{err}(\cdot)$ is the mean squared error and
$\mathrm{uncert}(\cdot)$ is the mean aleatoric uncertainty per bin.
Consideration of epistemic uncertainty is especially beneficial on smaller
data sets (BreastPathQ), where our approach outperforms Levi et al. (2019). On
larger data sets, the benefit diminishes and both approaches are equally
calibrated.

Additionally, we report UCE values from a DenseNet ensemble for comparison. In
contrast to what is reported by Lakshminarayanan et al. (2017), the deep
ensemble tends to be calibrated worse. Only on BoneAge, the ensemble is better
calibrated prior to recalibration of the other methods. After recalibration,
both approaches outperform the deep ensemble.

Fig. 4 shows the result of intra-training calibration of aleatoric
uncertainty. It indicates that the gap between training and test loss is
successfully closed.

Figure 4: (Left) Intra-training calibration of aleatoric uncertainty with
$\sigma$ scaling. The deep model no longer underestimates $\hat{\sigma}^{2}$
on unseen test data. (Right) The MSE of predictive mean is higher and
$\sigma^{2}$ is underestimated. Note: Calibration is only applied at test
time.

pixel coordinates

Figure 5: Example result from EndoVis test set. The task is to predict pixel coordinates of the forceps shaft center. Before calibration, the uncertainty is underestimated and the true instrument position $\bm{y}$ does not fall into the prediction region $\hat{\bm{y}}\pm\hat{\Sigma}$. After calibration with $\sigma$ scaling, the uncertainty better reflects the predictive error. Table 1: Uncertainty calibration error test set results for different datasets and model architectures (averaged over 5 runs). High UCE values indicate miscalibration. In addition, the resulting $s$ for $\sigma$ scaling is given. We also report UCE values for an ensemble of DenseNets. Bold font indicates lowest values in each experiment. |  | Levi et al. | ours |
---|---|---|---|---
Data Set | Model | MSE | none | aux | $\sigma$ | $s$ | none | aux | $\sigma$ | $s$ | ensemble
| ResNet-101 | 6.4e-3 | 0.51 | 0.35 | 0.28 | 2.91 | 0.49 | 0.31 | 0.20 | 2.37 |
BreastPathQ | DenseNet-201 | 7.0e-3 | 0.21 | 0.38 | 0.15 | 1.62 | 0.11 | 0.36 | 0.15 | 1.33 | 0.51
| EfficientNet-B4 | 6.4e-3 | 0.49 | 0.65 | 0.10 | 2.30 | 0.46 | 0.47 | 0.17 | 1.77 |
| ResNet-101 | 5.3e-3 | 0.28 | 0.07 | 0.06 | 1.46 | 0.28 | 0.02 | 0.06 | 1.40 |
BoneAge | DenseNet-201 | 3.5e-3 | 0.31 | 0.05 | 0.05 | 2.98 | 0.31 | 0.05 | 0.05 | 2.54 | 0.09
| EfficientNet-B4 | 3.5e-3 | 0.30 | 0.05 | 0.10 | 4.83 | 0.30 | 0.03 | 0.12 | 3.98 |
| ResNet-101 | 4.0e-4 | 0.04 | 0.10 | 0.09 | 6.07 | 0.04 | 0.04 | 0.04 | 3.50 |
EndoVis | DenseNet-201 | 1.1e-3 | 0.09 | 0.05 | 0.05 | 3.24 | 0.04 | 0.04 | 0.04 | 2.57 | 0.08
| EfficientNet-B4 | 8.9e-4 | 0.06 | 0.05 | 0.06 | 2.25 | 0.06 | 0.04 | 0.04 | 1.79 |
| ResNet-101 | 2.0e-3 | 0.17 | 0.02 | 0.02 | 2.74 | 0.17 | 0.01 | 0.02 | 2.14 |
OCT | DenseNet-201 | 1.3e-3 | 0.08 | 0.01 | 0.02 | 1.60 | 0.04 | 0.03 | 0.02 | 1.26 | 0.67
| EfficientNet-B4 | 1.4e-3 | 0.12 | 0.01 | 0.01 | 2.65 | 0.12 | 0.01 | 0.01 | 1.94 |

###  4.1 Posterior Prediction Intervals

In addition to the calibration diagrams, we compute prediction intervals from
the uncalibrated and calibrated posterior predictive distribution. Well-
calibrated prediction intervals provide a reliable measure of precision of the
estimated target value. In Bayesian inference, prediction intervals define an
interval within which the true target value $\bm{y}^{\ast}$ of a new,
unobserved input $\bm{x}^{\ast}$ is expected to fall with a specific
probability (Heskes, 1997; Held and Sabanés Bové, 2014). This is also referred
to as the credible interval of the posterior predictive distribution. For
$\gamma\in(0,1)$, a $\gamma\cdot 100\,\%$ prediction interval is defined
through $z_{l}$ and $z_{u}$ such that

$\int_{z_{l}}^{z_{u}}p(\bm{y}^{\ast}\,|\,\bm{x}^{\ast},\mathcal{D})\,\mathrm{d}\bm{y}^{\ast}=\gamma~{},$
(16)

with posterior predictive distribution
$p(\bm{y}^{\ast}\,|\,\bm{x}^{\ast},\mathcal{D})$. We compute the 50 %, 90 %,
95 %, and 99 % prediction interval using the root of the predictive variance
from Eq. (10); that is, the $\hat{\bm{y}}\pm z\hat{\Sigma}$ intervals with
$z\in\\{\Phi(0.5),\Phi(0.9),\Phi(0.95),\Phi(0.99)\\}$ (estimated interval),
with probit function $\Phi(p)=\sqrt{2}\mathrm{erf}^{-1}(p)$ and
$\mathrm{erf}(p)$ is the Gaussian error function. This assumes that the
posterior predictive distribution is Gaussian, which is not generally the
case. To assess the calibration of the posterior prediction interval, we
compute the percentage of how many of the ground truth values of the test set
actually fall within the respective intervals (observed interval). In Fig. 6,
selected plots of observed vs. estimated prediction intervals are shown. A
complete list of prediction intervals can be found in Appendix G.1.

In general, the uncalibrated prediction intervals are estimated to be too
narrow, which is a direct consequence of the underestimated predictive
variance. For example, the uncalibrated 90 % interval on DenseNet-201/BoneAge
actually only contains approx. 50 % of the ground truth values. On this data
set, the prediction intervals are considerably improved after recalibration
(Fig. 6 left). If a network is already well-calibrated, recalibration can lead
to overestimation of the lower prediction intervals (Fig. 6 right). However,
in all cases, the 99 % prediction interval contains approx. 99 % of the ground
truth test set values after recalibration. This is not the case without the
proposed calibration methods. Fig. 5 shows a practical example of the
$\hat{\bm{y}}\pm\hat{\Sigma}$ prediction region from the EndoVis test set.
Even though the posterior predictive distribution is not necessarily Gaussian,
the calibrated results fit the prediction intervals well. This is especially
the case for BoneAge, which is the largest data set used in this paper.

Figure 6: Observed vs. estimated posterior prediction intervals on the test
sets. (Left & center) The uncalibrated prediction interval is too narrow due
to underestimation of uncertainty. (Right) Calibration can lead to
overestimation of predictive intervals, if the network is already well-
calibrated. Dashed lines denote 1:1 mapping.

###  4.2 Detection of Out-of-Distribution Data and Unreliable Predictions

Deep neural networks only yield reliable predictions for data which follow the
same distribution as the training data. A shift in distribution could occur
when a model trained on CT data from a specific CT device is applied to data
from another manufacturer’s CT device, for example. This could potentially
lead to wrong predictions with low uncertainty, which we tackle with
recalibration. To create a moderate distribution shift, we preprocess images
from the BoneAge data set using Contrast Limited Adapative Histogram
Equalization (CLAHE) (Pizer et al., 1987) with a clip-limit of 0.03 and report
histograms of the uncertainties (see Fig. 7). Additionally, a severe
distribution shift is created by presenting images from the BreastPathQ data
set to the models trained on BoneAge. Lakshminarayanan et al. (2017) state
that deep ensembles provide better-calibrated uncertainty than Bayesian neural
networks with MC dropout variational inference. We therefore train an ensemble
of 5 randomly initialized DenseNet-201 and compare Bayesian uncertainty with
$\sigma$ scaling to ensemble uncertainty under distribution shift. The results
with $\sigma$ scaling are comparable to those from a deep ensemble for a
moderate shift, but without the need to train multiple models on the same data
set. A severe shift leads to only slightly increased uncertainties from the
calibrated MC dropout model, while the deep ensemble is more sensitive.

Additionally, we apply the well-calibrated models to detect and reject
uncertain predictions, as crucial decisions in medical practice should only be
made on the basis of reliable predictions. An uncertainty threshold
$\Sigma_{\mathrm{max}}^{2}$ is defined and all predictions from the test set
are rejected where $\hat{\Sigma}^{2}>\Sigma_{\mathrm{max}}^{2}$ (see Fig. 8).
From this, a decrease in overall MSE is expected. We additionally compare
rejection on the basis of $\sigma$ scaled uncertainty to uncertainty from the
aforementioned ensemble. In case of $\sigma$ scaling, the test set MSE
decreases monotonically as a function of the uncertainty threshold, whereas
the ensemble initially shows an increasing MSE (see Fig. 8).

Figure 7: Histograms of the uncertainties for out-of-distribution detection
with DenseNet-201 on BoneAge test set. (Left) Uncertainties from a non-
Bayesian ensemble of five DenseNets and (right) Bayesian uncertainties
calibrated with $\sigma$ scaling. The distribution shifts have been created
with pre-processing by CLAHE (moderate) and images from a different domain
(severe).

Figure 8: Rejection of uncertain predictions with DenseNet-201 on BoneAge test
set with $\hat{\Sigma}^{2}>\Sigma^{2}_{\mathrm{max}}$. The shaded area width
visualizes the percentage of rejected samples. The dashed line visualizes
linear relationship.

##  5 Discussion & Conclusion

In this paper, well-calibrated predictive uncertainty in medical imaging
obtained by variational inference with deep Bayesian models is discussed. Both
aux and $\sigma$ scaling calibration methods considerably reduce
miscalibration of predictive uncertainty in terms of UCE. If the deep model is
already well-calibrated, $\sigma$ scaling does not negatively affect the
calibration, which results in $s\rightarrow 1$. More complex calibration
methods such as aux scaling have to be used with caution, as they can overfit
the data set used for calibration. If the calibration set is sufficiently
large, they can outperform simple scaling. However, models trained on large
data sets are generally better calibrated and the benefit diminishes. Compared
to the work of Levi et al. (2019), accounting for epistemic uncertainty is
particularly beneficial for smaller data sets, which is helpful in medical
practice where access to large labeled data sets is less common and is
associated with great costs.

Posterior prediction intervals provide another insight into the calibration of
deep models. After recalibration, the 99 % posterior prediction intervals
correctly contain approx. 99 % of the ground truth test set values. In some
cases, lower prediction intervals are estimated to be too wide after
calibration. This is especially the case for smaller data sets and we
conjecture that small calibration sets may not contain enough i.i.d. data for
calibrating lower prediction intervals and that the assumption of a Gaussian
predictive distribution is too strong in this case. On the smallest data set
BreastPathQ, aux scaling seems to perform better in terms of prediction
intervals, but not in terms of UCE.

Well-calibrated uncertainties from MC dropout are able to detect a moderate
shift in the data distribution. However, deep ensembles perform better under a
severe distribution shift. BNNs with calibrated uncertainty by $\sigma$
scaling outperform ensemble uncertainty in the rejection task, which we
attribute to the generally poorer calibration of ensembles on in-distribution
data.

$\sigma$ scaling is simple to implement, does not change the predictive mean
$\hat{\bm{y}}$, and therefore guarantees to conserve the model’s accuracy. It
is preferable to regularization (e. g. early stopping) or more complex
recalibration methods in calibrated uncertainty estimation with Bayesian deep
learning. The disconnection between training and test NLL can successfully be
closed, which creates highly accurate models with reliable uncertainty
estimates. However, there are many factors (e. g. network capacity, weight
decay, dropout configuration) influencing the uncertainty that have not been
discussed here and will be addressed in future work.

Acknowledgments

We thank Vincent Modes and Mark Wielitzka for their insightful comments. This
research has received funding from the European Union as being part of the
ERDF OPhonLas project.

## References

  * Bishop (2006) Christopher M. Bishop.  _Pattern Recognition and Machine Learning_.  Springer, 2006.  ISBN 978-0-387-31073-2.
  * Bragman et al. (2018) Felix JS Bragman, Ryutaro Tanno, Zach Eaton-Rosen, Wenqi Li, David J Hawkes, Sebastien Ourselin, Daniel C Alexander, Jamie R McClelland, and M Jorge Cardoso.  Uncertainty in multitask learning: joint representations for probabilistic mr-only radiotherapy planning.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 3–11, 2018.
  * Dalca et al. (2019) Adrian V. Dalca, Guha Balakrishnan, John Guttag, and Mert R. Sabuncu.  Unsupervised learning of probabilistic diffeomorphic registration for images and surfaces.  _Med Image Anal_ , 57:226–236, 2019.  doi: https://doi.org/10.1016/j.media.2019.07.006.
  * Gal and Ghahramani (2016) Yarin Gal and Zoubin Ghahramani.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In _ICML_ , pages 1050–1059, 2016.
  * Gessert et al. (2018) Nils Gessert, Matthias Schlüter, and Alexander Schlaefer.  A deep learning approach for pose estimation from volumetric OCT data.  _Med Image Anal_ , 46:162–179, 2018.  doi: 10.1016/j.media.2018.03.002.
  * Guo et al. (2017) Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q. Weinberger.  On calibration of modern neural networks.  In _ICML_ , pages 1321–1330, 2017.
  * Halabi et al. (2019) Safwan S. Halabi, Luciano M. Prevedello, Jayashree Kalpathy-Cramer, Artem B. Mamonov, Alexander Bilbily, Mark Cicero, Ian Pan, Lucas Araújo Pereira, Rafael Teixeira Sousa, Nitamar Abdala, Felipe Campos Kitamura, Hans H. Thodberg, Leon Chen, George Shih, Katherine Andriole, Marc D. Kohli, Bradley J. Erickson, and Adam E. Flanders.  The RSNA pediatric bone age machine learning challenge.  _Radiol_ , 290(2):498–503, 2019.  doi: 10.1148/radiol.2018180736.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _CVPR_ , pages 770–778, 2016.
  * Held and Sabanés Bové (2014) Leonhard Held and Daniel Sabanés Bové.  _Applied Statistical Inference_.  Springer, 1 edition, 2014.  ISBN 978-3-642-37887-4.  doi: 10.1007/978-3-642-37887-4.
  * Heskes (1997) Tom Heskes.  Practical confidence and prediction intervals.  In _NeurIPS_ , pages 176–182, 1997.
  * Hora (1996) Stephen C Hora.  Aleatory and epistemic uncertainty in probability elicitation with an example from hazardous waste management.  _Reliability Engineering & System Safety_, 54(2-3):217–223, 1996.
  * Huang et al. (2017) Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q. Weinberger.  Densely connected convolutional networks.  In _CVPR_ , pages 4700–4708, 2017.
  * Kendall and Gal (2017) Alex Kendall and Yarin Gal.  What uncertainties do we need in bayesian deep learning for computer vision?  In _NeurIPS_ , pages 5574–5584, 2017.
  * Kingma and Welling (2014) Diederik P. Kingma and Max Welling.  Auto-encoding variational bayes.  In _ICLR_ , 2014.
  * Kuleshov et al. (2018) Volodymyr Kuleshov, Nathan Fenner, and Stefano Ermon.  Accurate uncertainties for deep learning using calibrated regression.  In _ICML_ , volume 80, pages 2796–2804, 2018.
  * Lakshminarayanan et al. (2017) Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell.  Simple and scalable predictive uncertainty estimation using deep ensembles.  In _NeurIPS_ , pages 6402–6413, 2017.
  * Laves et al. (2017) Max-Heinrich Laves, Andreas Schoob, Lüder A. Kahrs, Tom Pfeiffer, Robert Huber, and Tobias Ortmaier.  Feature tracking for automated volume of interest stabilization on 4D-OCT images.  In _SPIE Medical Imaging_ , volume 10135, pages 256–262, 2017.  doi: 10.1117/12.2255090.
  * Laves et al. (2019) Max-Heinrich Laves, Sontje Ihler, Karl-Philipp Kortmann, and Tobias Ortmaier.  Well-calibrated model uncertainty with temperature scaling for dropout variational inference.  In _Bayesian Deep Learning Workshop (NeurIPS)_ , 2019.  arXiv:1909.13550.
  * Laves et al. (2020) Max-Heinrich Laves, Sontje Ihler, Jacob F Fast, Lüder A Kahrs, and Tobias Ortmaier.  Well-calibrated regression uncertainty in medical imaging with deep learning.  In _Medical Imaging with Deep Learning_ , 2020.
  * Levi et al. (2019) Dan Levi, Liran Gispan, Niv Giladi, and Ethan Fetaya.  Evaluating and calibrating uncertainty prediction in regression tasks.  In _arXiv_ , 2019.  arXiv:1905.11659.
  * Luo et al. (2019) Jie Luo, Alireza Sedghi, Karteek Popuri, Dana Cobzas, Miaomiao Zhang, Frank Preiswerk, Matthew Toews, Alexandra Golby, Masashi Sugiyama, William M Wells, et al.  On the applicability of registration uncertainty.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 410–419, 2019.
  * Martel et al. (2019) A. L. Martel, S. Nofech-Mozes, S. Salama, S. Akbar, and M. Peikari.  Assessment of residual breast cancer cellularity after neoadjuvant chemotherapy using digital pathology [data set].  _The Cancer Imaging Archive_ , 2019.  doi: 10.7937/TCIA.2019.4YIBTJNO.
  * Messay et al. (2015) Temesguen Messay, Russell C. Hardie, and Timothy R. Tuinstra.  Segmentation of pulmonary nodules in computed tomography using a regression neural network approach and its application to the lung image database consortium and image database resource initiative dataset.  _Med Image Anal_ , 22(1):48–62, 2015.  doi: 10.1016/j.media.2015.02.002.
  * Nix and Weigend (1994) David A Nix and Andreas S Weigend.  Estimating the mean and variance of the target probability distribution.  In _Proceedings of IEEE International Conference on Neural Networks_ , volume 1, pages 55–60, 1994.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In _NeurIPS_ , pages 8024–8035, 2019.
  * Payer et al. (2019) Christian Payer, Darko Štern, Horst Bischof, and Martin Urschler.  Integrating spatial configuration into heatmap regression based CNNs for landmark localization.  _Med Image Anal_ , 54:207–219, 2019.  doi: 10.1016/j.media.2019.03.007.
  * Phan et al. (2018) Buu Phan, Rick Salay, Krzysztof Czarnecki, Vahdat Abdelzad, Taylor Denouden, and Sachin Vernekar.  Calibrating uncertainties in object localization task.  In _Bayesian Deep Learning Workshop (NeurIPS)_ , 2018.  arXiv:1811.11210.
  * Pizer et al. (1987) Stephen M Pizer, E Philip Amburn, John D Austin, Robert Cromartie, Ari Geselowitz, Trey Greer, Bart ter Haar Romeny, John B Zimmerman, and Karel Zuiderveld.  Adaptive histogram equalization and its variations.  _Computer vision, graphics, and image processing_ , 39(3):355–368, 1987.
  * Schlemper et al. (2018) Jo Schlemper, Guang Yang, Pedro Ferreira, Andrew Scott, Laura-Ann McGill, Zohya Khalique, Margarita Gorodezky, Malte Roehl, Jennifer Keegan, Dudley Pennell, et al.  Stochastic deep compressive sensing for the reconstruction of diffusion tensor cardiac mri.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 295–303, 2018.
  * Srivastava et al. (2014) Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  _JMLR_ , 15:1929–1958, 2014.
  * Štern et al. (2016) Darko Štern, Christian Payer, Vincent Lepetit, and Martin Urschler.  Automated age estimation from hand mri volumes using deep learning.  In _MICCAI_ , pages 194–202, 2016.
  * Tan et al. (2017) Li Kuo Tan, Yih Miin Liew, Einly Lim, and Robert A. McLaughlin.  Convolutional neural network regression for short-axis left ventricle segmentation in cardiac cine MR sequences.  _Med Image Anal_ , 39:78–86, 2017.  doi: 10.1016/j.media.2017.04.002.
  * Tan and Le (2019) Mingxing Tan and Quoc V Le.  EfficientNet: Rethinking model scaling for convolutional neural networks.  In _ICML_ , 2019.
  * Tanno et al. (2016) Ryutaro Tanno, Aurobrata Ghosh, Francesco Grussu, Enrico Kaden, Antonio Criminisi, and Daniel C Alexander.  Bayesian image quality transfer.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 265–273, 2016.
  * Tanno et al. (2017) Ryutaro Tanno, Daniel E Worrall, Aurobrata Ghosh, Enrico Kaden, Stamatios N Sotiropoulos, Antonio Criminisi, and Daniel C Alexander.  Bayesian image quality transfer with cnns: exploring uncertainty in dmri super-resolution.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 611–619, 2017.
  * Xie et al. (2018) Yuanpu Xie, Fuyong Xing, Xiaoshuang Shi, Xiangfei Kong, Hai Su, and Lin Yang.  Efficient and robust cell detection: A structured regression approach.  _Med Image Anal_ , 44:245–254, 2018.  doi: 10.1016/j.media.2017.07.003.
  * Yin et al. (2020) Shi Yin, Qinmu Peng, Hongming Li, Zhengqiang Zhang, Xinge You, Katherine Fischer, Susan L. Furth, Gregory E. Tasian, and Yong Fan.  Automatic kidney segmentation in ultrasound images using subsequent boundary distance regression and pixelwise classification networks.  _Med Image Anal_ , 60:101602, 2020.  doi: 10.1016/j.media.2019.101602.

##  A Laplacian Model

Using $\mathsf{Laplace}(\hat{\bm{y}}(\bm{x}),\hat{\sigma}(\bm{x}))$ as model,
the conditional log-likelihood is given by

$\displaystyle\log p(\bm{Y}\,|\,\bm{X},\bm{\theta})=$
$\displaystyle\sum_{i=1}^{m}\log\left(\frac{1}{2\hat{\sigma}^{(i)}_{\bm{\theta}}}\exp\left\\{-\frac{\big{|}\bm{y}^{(i)}-\hat{\bm{y}}_{\bm{\theta}}^{(i)}\big{|}}{\hat{\sigma}^{(i)}_{\bm{\theta}}}\right\\}\right)~{},$
(17)

which results in the following minimization criterion:

$\mathcal{L}_{\mathrm{L}}(\bm{\theta})=\sum_{i=1}^{m}\frac{1}{\hat{\sigma}^{(i)}_{\bm{\theta}}}\big{|}\bm{y}^{(i)}-\hat{\bm{y}}_{\bm{\theta}}^{(i)}\big{|}+\log\big{(}\hat{\sigma}_{\bm{\theta}}^{(i)}\big{)}~{}.$
(18)

Using $\mathcal{L}_{\mathrm{L}}(\bm{\theta})$ instead of
$\mathcal{L}_{\mathrm{G}}(\bm{\theta})$ results in applying an L1 metric on
the predictive mean. In some cases, this led to better results. However, we
have not conducted extensive experiments with it and leave it to future work.

##  B Derivation of $\sigma$ Scaling

See § 2.3. Using a Gaussian model, we scale the standard deviation $\sigma$
with a scalar value $s$ to calibrate the probability density function

$p\left(\bm{y}|\bm{x};\hat{\bm{y}}(x),\hat{\sigma}^{2}(x)\right)=\mathcal{N}\left(\bm{y};\hat{\bm{y}}(x),(s\cdot\hat{\sigma}(x))^{2}\right)~{}.$
(19)

The conditional log-likelihood is given by

$\displaystyle\log p(\bm{Y}\,|\,\bm{X},\bm{\theta})$
$\displaystyle=\sum_{i=1}^{m}\log\left(\frac{1}{\sqrt{2\pi}s\hat{\sigma}_{\theta}^{(i)}}\exp\left(\frac{\big{\|}\bm{y}^{(i)}-\hat{\bm{y}}^{(i)}_{\bm{\theta}}\big{\|}^{2}}{2\left(s\hat{\sigma}_{\bm{\theta}}^{(i)}\right)^{2}}\right)\right)$
(20)
$\displaystyle=-\dfrac{m}{2}\log\left(2\pi\right)-\sum_{i=1}^{m}\log\left(s\hat{\sigma}^{(i)}_{\bm{\theta}}\right)+\frac{1}{2}\left(s\hat{\sigma}_{\bm{\theta}}^{(i)}\right)^{-2}\cdot\big{\|}\bm{y}^{(i)}-\hat{\bm{y}}_{\bm{\theta}}^{(i)}\big{\|}^{2}$
(21)

This results in the following optimization objective (ignoring constants):

$\mathcal{L}_{\mathrm{G}}(s)=m\log(s)+\tfrac{1}{2}s^{-2}\sum_{i=1}^{m}(\hat{\sigma}_{\bm{\theta}}^{(i)})^{-2}\big{\|}\bm{y}^{(i)}-\hat{\bm{y}}_{\bm{\theta}}^{(i)}\big{\|}^{2}~{}.$
(22)

Using a Laplacian model, the optimization criterion follows as

$\mathcal{L}_{\mathrm{L}}(s)=m\log(s)+s^{-1}\sum_{i=1}^{m}\frac{1}{\hat{\sigma}_{\bm{\theta}}^{(i)}}\left|\bm{y}^{(i)}-\hat{\bm{y}}_{\bm{\theta}}^{(i)}\right|~{}.$
(23)

Eq. (22) and (23) are optimized w.r.t. $s$ with fixed $\bm{\theta}$ using
gradient descent in a separate calibration phase after training. The solution
to Eq. (22) can also be written in closed form as

$s_{\mathrm{G}}=\pm\sqrt{\frac{1}{m}\sum_{i=1}^{m}\big{(}\hat{\sigma}_{\bm{\theta}}^{(i)}\big{)}^{-2}\big{\|}\bm{y}^{(i)}-\hat{\bm{y}}_{\bm{\theta}}^{(i)}\big{\|}^{2}}$
(24)

and the solution to Eq. (23) follows as

$s_{\mathrm{L}}=\frac{1}{m}\sum_{i=1}^{m}\frac{1}{\hat{\sigma}_{\bm{\theta}}^{(i)}}\left|\bm{y}^{(i)}-\hat{\bm{y}}_{\bm{\theta}}^{(i)}\right|~{},$
(25)

respectively. We apply $\sigma$ scaling to jointly calibrate aleatoric and
epistemic uncertainty as described in § 2.4.

##  C Unbiased Estimator of the Approximate Predictive Variance

We show that the expectation of the predictive sample variance from MC
dropout, as given in (Kendall and Gal, 2017), equals the true variance of the
approximate posterior predictive distribution.

###### Proposition 1

Given $N$ MC dropout samples
$\bm{f}_{\bm{\theta}_{n}}=[\hat{\bm{y}}_{n},\hat{\sigma}^{2}_{n}]$ from our
approximate predictive distribution
$p(\bm{y}^{\ast}|\bm{x}^{\ast},\mathcal{D})=\mathcal{N}(\bm{y}^{\ast};\bm{y},\Sigma^{2})$,
the predictive sample variance

$\hat{\Sigma}^{2}=\frac{1}{N}\sum_{n=1}^{N}\left(\hat{\bm{y}}_{n}-\frac{1}{N}\sum_{n=1}^{N}\hat{\bm{y}}_{n}\right)^{2}+\frac{1}{N}\sum_{n=1}^{N}\hat{\sigma}^{2}_{n}$
(26)

is an unbiased estimator of the approximate predictive variance.

Proof

$\displaystyle\mathbb{E}\left[\hat{\Sigma}^{2}\right]$
$\displaystyle=\mathbb{E}\left[\frac{1}{N}\sum_{n=1}^{N}\left(\hat{\bm{y}}_{n}-\frac{1}{N}\sum_{n=1}^{N}\hat{\bm{y}}_{n}\right)^{2}+\frac{1}{N}\sum_{n=1}^{N}\hat{\sigma}^{2}_{n}\right]$
(27)
$\displaystyle=\mathbb{E}\left[\frac{1}{N}\sum_{n=1}^{N}\left(\hat{\bm{y}}_{n}-\frac{1}{N}\sum_{n=1}^{N}\hat{\bm{y}}_{n}\right)^{2}\right]+\mathbb{E}\left[\frac{1}{N}\sum_{n=1}^{N}\hat{\sigma}^{2}_{n}\right]$
(28)
$\displaystyle\mathrm{with}\quad\frac{1}{N}\sum_{n=1}^{N}\hat{\bm{y}}_{n}=\bar{\bm{y}}\quad\mathrm{follows}$
(29)
$\displaystyle=\mathbb{E}\left[\frac{1}{N}\sum_{n=1}^{N}\left(\hat{\bm{y}}_{n}-\bar{\bm{y}}\right)^{2}\right]+\hat{\sigma}^{2}$
(30)
$\displaystyle=\mathbb{E}\left[\frac{1}{N}\sum_{n=1}^{N}\left(\hat{\bm{y}}_{n}-\bar{\bm{y}}\right)^{2}+\bar{\bm{y}}^{2}-\bar{\bm{y}}^{2}+\bm{y}^{2}-\bm{y}^{2}+2\bar{\bm{y}}\bm{y}-2\bar{\bm{y}}\bm{y}\right]+\hat{\sigma}^{2}$
(31)
$\displaystyle=\mathbb{E}\left[\frac{1}{N}\sum_{n=1}^{N}\left(\hat{\bm{y}}_{n}-\bm{y}\right)^{2}-\left(\bar{\bm{y}}-\bm{y}\right)^{2}\right]+\hat{\sigma}^{2}$
(32)
$\displaystyle=\mathbb{E}\left[\left(\hat{\bm{y}}-\bm{y}\right)^{2}\right]-\mathbb{E}\left[\left(\bar{\bm{y}}-\bm{y}\right)^{2}\right]+\hat{\sigma}^{2}$
$\displaystyle=\Sigma^{2}-\hat{\sigma}^{2}+\hat{\sigma}^{2}$ (33)
$\displaystyle\mathbb{E}\left[\hat{\Sigma}^{2}\right]$
$\displaystyle=\Sigma^{2}$ (34)

Note that the predicted heteroscedastic aleatoric uncertainty
$\hat{\sigma}^{2}$ equals the bias $\mathbb{E}[(\bar{\bm{y}}-\bm{y})^{2}]$ in
Eq. (33) when the aleatoric uncertainty is perfectly calibrated, thus
$\mathbb{E}[(\bar{\bm{y}}-\bm{y})^{2}]=\hat{\sigma}^{2}$.

##  D Training Procedure

The model implementations from PyTorch 1.3 (Paszke et al., 2019) are used and
trained with the following settings:

  * •

training for 500 epochs with batch size of 16

  * •

Adam optimizer with initial learn rate of $3\cdot 10^{-4}$ and weight decay
with $\lambda=10^{-7}$

  * •

reduce-on-plateau learn rate scheduler (patience of 20 epochs) with factor of
0.1

  * •

in MC dropout, $N=25$ forward passes were performed with dropout with $p=0.5$
used for ResNet (as described in (Gal and Ghahramani, 2016)). In DenseNet
($p=0.2$) and EfficientNet ($p=0.4$) standard dropout $p$ of the architecture
is used.

  * •

Additional validation and test sets are used if provided by the data sets;
otherwise, a train/validation/test split of approx. 50% / 25% / 25% is used

  * •

Source code for all experiments is available at github.com/mlaves/well-
calibrated-regression-uncertainty

##  E 3D OCT Needle Pose Data Set

pixel coordinates

Figure 9: Example image from OCT data set showing $\operatorname*{arg\,max}$
projections of a surgical needle tip acquired by optical coherence tomography.

Our data set was created by attaching a surgical needle to a high-precision
six-axis hexapod robot (H-826, Physik Instrumente GmbH & Co. KG, Germany) and
observing the needle tip with 3D optical coherence tomography (OCS1300SS,
Thorlabs Inc., USA). The data set consists of 5,000 OCT acquisitions with
$(64\times 64\times 512)$ voxels, covering a volume of approx. $(3\times
3\times 3)$ $\mathrm{mm}^{3}$. Each acquisition is taken at a different robot
configuration and labeled with the corresponding 6DoF pose
$\bm{y}\in\mathbb{R}^{6}$. To process the volumetric data with CNNs for planar
images, we calculate 3 planar projections along the spatial dimensions using
the $\operatorname*{arg\,max}$ operator, scale them to equal size and stack
them together as three-channel image (see Fig. 9). A similar approach was
presented in (Laves et al., 2017) and (Gessert et al., 2018). The data are
characterized by a high amount of speckle noise, which is a typical phenomenon
in optical coherence tomography. The data set is publicly available at
github.com/mlaves/3doct-pose-dataset.

##  F Ablation Study on Auxiliary Model Scaling

Figure 10: Calibration diagrams for aux scaling with different number of
hidden layer units $h$ on BreastPathQ/DenseNet-201.

Here, we investigate the overfitting behavior of aux scaling by reducing the
number of hidden layer units $h$ of the two-layer auxiliary model with ReLU
activations. Aux scaling is more powerful than $\sigma$ scaling, which can
lead to overfitting the calibration set. Fig. 10 shows calibration diagrams
for the auxiliary model ablations. Reducing $h$ leads to a minor calibration
improvement, but at $h=2$, the model outputs a constant uncertainty, which is
close to the overall mean of the observed uncertainty. A single-layer single-
unit model without bias would be equivalent to $\sigma$ scaling.

##  G Additional Results and Calibration Diagrams

All test set runs have been repeated 5 times. Solid lines denote mean and
shaded areas denote standard deviation calculated from the repeated runs.

Table 2: Negative log-likelihood test set results for different datasets and model architectures (averaged over 5 runs). High NLL values indicate miscalibration. We also report NLL values for an ensemble of DenseNets. Bold font indicates lowest values in each experiment. |  | Levi et al. | ours |
---|---|---|---|---
Data Set | Model | MSE | none | aux | $\sigma$ | none | aux | $\sigma$ | ensemble
| ResNet-101 | 6.4e-3 | -0.78 | -5.06 | -5.06 | -2.89 | -5.17 | -5.16 |
BreastPathQ | DenseNet-201 | 7.0e-3 | -5.16 | -5.84 | -5.70 | -5.67 | -6.03 | -5.78 | 0.11
| EfficientNet-B4 | 6.4e-3 | -3.11 | -5.99 | -5.53 | -4.73 | -6.16 | -5.62 |
| ResNet-101 | 5.3e-3 | -3.90 | -4.34 | -4.34 | -3.99 | -4.34 | -4.34 |
BoneAge | DenseNet-201 | 3.5e-3 | 1.74 | -4.70 | -4.69 | -0.75 | -4.70 | -4.69 | 0.07
| EfficientNet-B4 | 3.5e-3 | 13.61 | -4.74 | -4.67 | 6.40 | -4.75 | -4.64 |
| ResNet-101 | 4.0e-4 | -0.53 | -6.32 | -6.33 | -3.85 | -6.76 | -6.72 |
EndoVis | DenseNet-201 | 1.1e-3 | -0.72 | -6.10 | -5.99 | -4.94 | -6.05 | -6.04 | 0.04
| EfficientNet-B4 | 8.9e-4 | -5.10 | -6.06 | -6.07 | -5.94 | -6.17 | -6.17 |
| ResNet-101 | 2.0e-3 | -1.08 | -5.24 | -5.24 | -3.38 | -5.24 | -5.24 |
OCT | DenseNet-201 | 1.3e-3 | -5.05 | -5.61 | -5.61 | -5.51 | -5.62 | -5.61 | 0.10
| EfficientNet-B4 | 1.4e-3 | -1.72 | -5.58 | -5.57 | -4.25 | -5.58 | -5.57 |

###  G.1 Additional Prediction Intervals

Figure 11: Observed vs. estimated posterior prediction intervals for all
networks/data sets.