Honest calibration assessment for binary outcome predictions

Honest calibration assessment for binary outcome predictions

\affil[1]{Heidelberg University, Germany} \affil[2]{University of Bern, Switzerland} \affil[3]{ETH Zürich, Switzerland} \affil[4]{University of Hohenheim, Germany}

\affil[ ]{ [email removed], [email removed], [email removed], [email removed], [email removed] }

Introduction

Consider first a univariate regression setting with fixed real covariates $x_1 \le \cdots \le x_n$ and independent binary observations $Y_1, \dots, Y_n \in \{0, 1\}$, where $\mathrm{pr}(Y_i = 1) = p(x_i)$ for some unknown regression function $p : \mathbb{R} \to [0,1]$. Standard parametric models for this setting, e.g.\ logistic or probit regression, involve monotone regression functions $p$. Thus, an interesting nonparametric alternative would be to draw inference on $p$ under the sole assumption that it is isotonic on $\mathbb{R}$,

In the specific applications we have in mind, the $x_i$ are themselves probability predictions for the binary outcomes, i.e. $x_i \in [0,1]$ is a prediction for the probability of the event $\{Y_i = 1\}$. In practice, the predictions can be obtained from a test sample of binary regressions, machine learning methods, or any other statistical model for binary data. A reliable interpretation of these predictions relies on the property of calibration, meaning that if the value $x_i$ is predicted, the corresponding event should indeed occur with probability $x_i$. In this setting, the regression function $p$ is called calibration curve, and it maps the predicted probabilities $x_i$ to the actual, or recalibrated, event probabilities $p(x_i) = \mathrm{pr}(Y_i = 1)$. For calibrated predictions, the calibration curve equals the diagonal, $p(x) = x$ for all $x \in [0,1]$. Drawing inference about $p$ thus allows to assess the calibration of the predictions.

Testing the null hypothesis of calibration, $\mathbb{H}_0 \colon p(x) = x$ for all $x \in [0,1]$, is closely related to goodness-of-fit testing, which is crucial in applications, see e.g., Tutz2011 and Hosmer2013Book. It is still regularly carried out by the classical test of HosmerLemeshow1980, which groups the predictions $x_i$ into bins and applies a $\chi^2$-test. It is however subject to multiple criticisms: First, its ad hoc choice of bins can result in untenable instabilities Bertolini2000, Allison2014. Second, placing the hypothesis of calibration in the null only allows for rejecting calibration rather than showing that a model is sufficiently well calibrated, where the latter would be highly desirable for applied researchers. Third, the test rejects essentially all, even acceptably well-specified models in large samples Nattino2020, Paul2013, resulting in calls for a goodness-of-fit tests with inverted hypotheses Nattino2020Rejoinder, that is, tests where the hypothesis $p(x) = x$ is contained in the alternative.

We propose a statistically sound solution to these criticisms by constructing honest, simultaneous confidence bands $(L^\alpha,U^\alpha)$ for the function $p$. That is, for a given small number $\alpha \in (0,1)$ and $\mathcal{Y} := (Y_i)_{i=1}^n$, we compute data-dependent functions $L^\alpha = L^\alpha(\cdot,\mathcal{Y})$ and $U^\alpha = U^\alpha(\cdot,\mathcal{Y})$ on $\mathbb{R}$ such that

In the context of calibration assessment, the functions $p, L^\alpha, U^\alpha$ are defined on $[0,1]$, and we call $(L^\alpha,U^\alpha)$ a calibration band, which is hence a confidence band for the calibration curve. It allows for the desirable conclusion that with confidence $1-\alpha$, the true calibration curve $p$ lies inside the band, simultaneously for all values of the predicted probabilities. This nests a classical goodness-of-fit test with $\mathbb{H}_0 \colon p(x) = x$ by checking whether the band contains the diagonal $d(x) = x$ for all relevant values $x \in [0,1]$, but also any other hypothesis on the calibration curve such as e.g., an inverted goodness-of-fit test with $\mathbb{H}_0 \colon |p(x) - x| > \epsilon$ for some small $\epsilon>0$. Hence, this band resolves the above mentioned criticisms of classical goodness-of-fit tests.

Figure (ref) shows the bands in a large data example for probit model predictions for the binary outcome of a fetus having a low birth weight. See Section (ref) for additional details. The test of Hosmer and Lemeshow clearly rejects calibration even though our bands indicate a well-calibrated model by including the diagonal line for all values in the unit interval. The magnified right panel of the figure shows that with confidence $1-\alpha$, the model is remarkably well calibrated for the most important region of small probability predictions in this application. It is important to notice that even though we build our bands on the model predictions, the methodology applies equally to both, causal and predictive regressions. An open-source implementation in the statistical software R R2022 is available under \href{https://cran.r-project.org/package=calibrationband}{https://cran.r-project.org/package=calibrationband}.

Our confidence bands are valid in finite samples subject only to the mild monotonicity assumption at (ref), implying that higher predictions entail a higher probability for $\{Y=1\}$, which is natural in the context of assessing calibration as already argued in DGJ_2021, Roelofs2020. For classical goodness-of-fit tests, the null hypothesis $p(x) = x$ already nests the monotonicity assumption and if a researcher aims to demonstrate calibration, i.e., $p(x) = x$ holds at least approximately, it is unlikely that there are strong deviations from isotonicity. Moreover, our confidence bands allow to detect and quantify violations of monotonicity as described in Appendix (ref). A non-monotonic $p$ may lead to a crossing of the lower and upper bound, i.e., $U^{\alpha}(x) < L^{\alpha}(x)$ for some $x$, which allows to reject monotonicity at level $\alpha$. This is supported by the graphical display that reacts to non-isotonicity by generating elongated horizontal segments in both, the isotonic regression estimate and the confidence bands. Finally, deriving confidence bands without any assumption on $p$ seems unrealistic and the assumption of monotonicity is relatively weak, e.g., in comparison to the parametric one used in Nattino2014.

As expected for a non-parametric, pathwise and almost universally valid confidence band, we require large data sets of at least $5\,000$ observations to obtain sensibly narrow bands. These are exactly the sample sizes where the classical goodness-of-fit tests become uninformative by rejecting all models in applications, see the simulation study of Kramer2007.

A theoretical analysis shows that the proposed confidence band adapts locally to the smoothness of the function $p$ and to the variance of the observations. Adaptivity to the smoothness means that the width of the bands decreases faster with the sample size $n$ in regions where $p$ is constant, and at a slower rate where $p$ is steeper. This property is known for more general confidence bands for a monotone mean function developed by Yang2019. Adaptivity to the variance means that the band is substantially narrower at $x$ if $p(x)$ is close to zero or one, compared to $p(x)$ near $0.5$. In many practical applications, including the low birth weight predictions analyzed in this article, predicted probabilities close to zero or one are of most relevance and a sharp assessment of calibration in these regions is particularly important.

Existing methods for the construction of confidence bands in this setting are rare with the following two exceptions: First, Nattino2014 propose the use of confidence bands based on a parametric assumption on the function $p$, which we show to have incorrect coverage in almost all of our simulation settings. Second, the nonparametric bands of Yang2019 are valid, in a modified sense even in settings where the isotonicity assumption (ref) is violated. But they are shown to be wider than our bands in theory and simulations.

We explain the absence of competing methods by their theoretical difficulties. Using asymptotic theory of the isotonic regression estimator is complicated as it requires the estimation of nuisance quantities such as the derivative of the unknown function $p$, the convergence rate depends on the functional form of $p$, it is subject to more restrictive assumptions and only results in bands with a pointwise interpretation Wright1981. Resampling schemes are theoretically found to be inconsistent for the isotonic regression Sen2010, GuntuboyinaSen2018. Other non-parametric approaches in the literature for constructing confidence bands for functions, many of them presented in the review by Hall2013, are often pointwise, not simultaneous, and require the selection of tuning parameters that may lead to instabilities, similar to the choice of the bins in the Hosmer-Lemeshow test. In contrast, the confidence bands proposed here are simple to compute and do not involve any implementation decisions resulting in a stable and reproducible method as called for by Stodden2016, Yu2020.

Construction of the confidence bands

Within the regression setting, we construct confidence bands for the isotonic regression function $p$ by means of the classical confidence bounds of Clopper1934 for a binomial parameter. Suppose that $Z$ is a binomial random variable with parameters $m$ and $q \in [0,1]$. For $\delta \in (0,1)$ let

Here $\mathrm{pbin}(\cdot,m,\xi)$ denotes the distribution function of the binomial distribution with parameters $m$ and $\xi$, while $\mathrm{qbeta}(\cdot,a,b)$ stands for the quantile function of the beta distribution with parameters $a,b > 0$. Then \[ \mathrm{pr}\{q \le u^\delta(Z,m)\} \ge 1 - \delta \quad\text{and}\quad \mathrm{pr}\{q \ge \ell^\delta(Z,m)\} \ge 1 - \delta . \] For the representation of $\ell^\delta(Z,m)$ and $u^\delta(Z,m)$ in terms of beta quantiles, we refer to Johnson2005.

Assumption (ref) allows to construct confidence bands for $p$ as follows. With $p_i := p(x_i)$, for arbitrary indices $1 \le j \le k \le n$, the random sum \[ Z_{jk} = \sum_{i=j}^k Y_i \] is stochastically larger than a binomial random variable with parameters $n_{jk} = k-j+1$ and $p_j$, and it is stochastically smaller than a binomial variable with parameters $n_{jk}$ and $p_k$. Thus, as explained in Lemma (ref),

If we combine these bounds for all pairs $(j,k)$ in a given set $\mathcal{J}$ and use the assumption at (ref), then we may claim with confidence $1-2|\mathcal{J}| \delta$ that simultaneously for all $(j,k) \in \mathcal{J}$, \[ p(x) \le u^\delta(Z_{jk},n_{jk}) \ \ \forall \ x \le x_j, \qquad p(x) \ge \ell^\delta(Z_{jk},n_{jk}) \ \ \forall \ x \ge x_k . \] Specifically, let $\mathcal{J}$ be the set of all index pairs $(j,k)$ such that $j \le k$ and $x_{j-1} < x_j$ and $x_k < x_{k+1}$. If there are tied values in $(x_i)_{i=1}^n$, $\mathcal{J}$ selects the outermost indices of the tied values. Hence, if $\{x_1,\ldots,x_n\}$ contains $N \le n$ different points, then $|\mathcal{J}| = (N^2 + N)/2$. Consequently, for a given confidence level $1 - \alpha \in (0,1)$, we may combine the bounds $u^\delta(Z_{jk},n_{jk})$ and $\ell^\delta(Z_{jk},n_{jk})$ with $\delta = \alpha/(N^2+N)$ to obtain a first confidence band.

The functions $U^{\alpha,\mathrm{raw}}, L^{\alpha,\mathrm{raw}}$ are isotonic and piecewise constant. Precisely, with $x_0 := -\infty$ and $x_{n+1} := \infty$, we know that $U^{\alpha,\mathrm{raw}} = 1$ on $(x_n,\infty)$, $L^{\alpha,\mathrm{raw}} = 0$ on $(-\infty,x_1)$, and

for $i = 1, \dots, n$. Consequently, computing the band $(L^{\alpha,\mathrm{raw}}, U^{\alpha,\mathrm{raw}})$ amounts to determining the $2n$ numbers $L^{\alpha,\mathrm{raw}}(x_i)$ and $U^{\alpha,\mathrm{raw}}(x_i)$, $i = 1, \dots, n$.

The confidence band proposed in Theorem (ref) has two potential drawbacks. First, a natural nonparametric estimator for the function $p$ under the assumption (ref) is given by a minimizer $\hat{p}$ of $\sum_{i=1}^n \{h(x_i) - Y_i\}^2$ over all isotonic functions $h\colon [0,1] \to [0,1]$ DGJ_2021. This minimizer is unique on the set $\{x_1,\ldots,x_n\}$. But there is no guarantee that $L^{\alpha,\mathrm{raw}} \le \hat{p} \le U^{\alpha,\mathrm{raw}}$. Second, the upper and lower bounds in (ref) and (ref) may even cross, resulting in an empty, and hence, nonsensical confidence band. These problems can be dealt with by using the non-crossing confidence band $(L^{\alpha,\mathrm{nc}}, U^{\alpha,\mathrm{nc}})$ given by pointwise minima and maxima:

Obviously, $L^{\alpha,\mathrm{nc}} \le \hat{p} \le U^{\alpha,\mathrm{nc}}$ on $\mathbb{R}$. Our simulation experiments indicate that $(L^{\alpha,\mathrm{raw}}, U^{\alpha,\mathrm{raw}}) = (L^{\alpha,\mathrm{nc}},U^{\alpha,\mathrm{nc}})$ holds in almost all cases whenever $p$ satisfies (ref); see Section (ref) for details. The potential crossing of the two bounds in Theorem (ref) also has an advantage. It allows for inference about the non-isotonicity of $p$, see Appendix (ref).

A potential obstacle in the practical application of the confidence bands proposed in this section is that their computation requires $\mathcal{O}(N^2)$ steps. This can be relieved by using a smaller family of index pairs $(j, k)$ in the definition of the confidence band. Specifically, if for some fixed integer $K \geq 1$ differences in the covariate smaller than $K^{-1}$ are regarded as negligible, then one could define

such that only blocks of covariate values between $r/K$ and $s/K$, $r,s \in \mathbb{Z}$, are considered. The resulting band is still honest, can be computed in $\mathcal{O}(|\widetilde{\mathcal{J}}|)$ steps, and one can reduce the correction factor of the significance level in (ref) and (ref) from $N^2+N$ to $|\widetilde{\mathcal{J}}|$. The drawback is that the constant regions in $L^{\alpha}$ and $U^{\alpha}$ become larger, thereby limiting the adaptivity of the band, so the number $K$ should not be too small. We henceforth refer to the restricted choice of $\widetilde{\mathcal{J}}$ in (ref) as the rounding method. Section 1 in the Supplementary Material illustrates in simulations that the rounding method drastically decreases the computation time and even results in narrower bands for all but very steep regions of $p$.

Relation to Yang2019

The methods of Yang2019 may be adapted to the present regression setting with covariates $x_1 \le \cdots \le x_n$ as follows: With the isotonic estimator $\hat{p}$ introduced before, let \[ Z_{jk}^{\mathrm{iso}} = \sum_{i=j}^k \hat{p}(x_i) . \] Set

This defines a confidence band $(L^{\alpha,\mathrm{YB}},U^{\alpha,\mathrm{YB}})$ with the following property:

where $\tilde{p}\colon \mathbb{R} \to [0,1]$ is any fixed isotonic function minimizing $\sum_{i=1}^n \{\tilde{p}(x_i) - p_i\}^2$. Thus one obtains a confidence band with guaranteed coverage probability $1 - \alpha$ for an isotonic approximation of $p$, even if (ref) is violated. The proof of (ref) follows from the arguments of Yang2019, noting that the random variables $Y_i$ are sub-Gaussian with scale parameter $\sigma = 1/2$. Thus, $\mathbb{E} \exp(t(Y_i - p_i)) \le \exp(\sigma^2t^2/2)$ for all $t \in \mathbb{R}$, implying that for arbitrary $\eta \ge 0$, \[ \mathrm{pr}\{ \pm (Z_{jk} - \mathbb{E} Z_{jk}) \ge \eta\} \le \exp(-2 n_{jk} \eta^2) , \] see Hoeffding1963. The following result shows that the confidence bands $(L^{\alpha,\mathrm{raw}},U^{\alpha,\mathrm{raw}})$ and $(L^{\alpha,\mathrm{nc}},U^{\alpha,\mathrm{nc}})$ are always contained in the band $(L^{\alpha,\mathrm{YB}},U^{\alpha,\mathrm{YB}})$.

Recall that the inequalities $L^{\alpha,\mathrm{raw}} \leq U^{\alpha,\mathrm{raw}}$ do not hold in general, and a crossing of the bounds allows to reject isotonicity at level $\alpha$, see Appendix (ref). In contrast, the bands by Yang2019 always contain the isotonic estimator $\hat{p}$, and are guaranteed to cover an optimal isotonic approximation of $p$ with probability at least $1-\alpha$. For calibration testing, the possibility of rejecting isotonicity seems more desirable than information about an isotonic approximation of the calibration curve, whose interpretation may be unclear in practice. It should be mentioned, however, that the band $(L^{\alpha,\mathrm{YB}},U^{\alpha,\mathrm{YB}})$ has a computational advantage. For the computation of $U_i^{\alpha,\mathrm{YB}}$ in (ref), it suffices to take the minimum over endpoints of constancy regions of $\hat{p}$, that is, all $(j,k) \in \mathcal{J}$ such that $j = \min(s\colon x_s\ge x_i)$ and $\hat{p}(x_k) < \hat{p}(x_{k+1})$ or $k = n$, see Proposition (ref). Likewise, for the computation of $L_i^{\alpha,\mathrm{YB}}$ in (ref), it suffices to take the maximum over all $(j,k) \in \mathcal{J}$ such that $\hat{p}(x_{j-1}) < \hat{p}(x_j)$ or $j = 1$ and $k = \max(s\colon x_s \le x_i)$. While the computation of $(L^{\alpha,\mathrm{raw}},U^{\alpha,\mathrm{raw}})$ or $(L^{\alpha,\mathrm{nc}},U^{\alpha,\mathrm{nc}})$ requires $\mathcal{O}(N^2)$ steps, the following lemma, whose proof is in the Supplementary Material, implies that the computation of $(L^{\alpha,\mathrm{YB}},U^{\alpha,\mathrm{YB}})$ requires only $\mathcal{O}(N \min\{n^{2/3},N\})$ steps.

Theoretical properties of the confidence bands

This section illustrates consistency and adaptivity properties of the confidence band $(L_n^{\alpha,\mathrm{raw}},U_n^{\alpha,\mathrm{raw}})$, where the subscript $n$ indicates the sample size, and we consider a triangular scheme of observations $(x_i,Y_i) = (x_{ni},Y_{ni})$, $i = 1, \dots, n$. We are interested in situations in which the observed covariates $x_{ni}$ could be the realizations of the order statistics of a random sample. Thus we extend the framework of Yang2019 and consider the following assumption.

This assumption comprises the setting of Yang2019. Let $G$ be a differentiable distribution function on $[0,1]$ such that $G'$ is bounded away from $0$ on $[a_o,b_o]$. If $x_{ni} = G^{-1}(i/n)$ for $i = 1, \dots, n$, then it is satisfied for any $C_1 < \inf_{[a_o,b_o]} G'$ and arbitrary $C_2 > 0$. The arguments in Moesching2020 can be modified to show that if $x_{n1},\ldots,x_{nn}$ are the order statistics of $n$ independent random variables with distribution function $G$, then Assumption (ref) is satisfied almost surely, provided that $C_1, C_2 > 0$ are chosen appropriately.

Part (i) implies that if $p$ is constant on $[a_o,b_o]$, then for arbitrary fixed $a_o < a < b < b_o$, \[ \sup_{x \in [a_o,b]} \{U_n^{\alpha,\mathrm{raw}}(x) - p(x)\}^+ + \sup_{x \in [a_o,b]} \{p(x) - L_n^{\alpha,\mathrm{raw}}(x)\}^+ = \mathcal{O}_p(\rho_n^{1/2}). \] Thus, parts (i-ii) of this theorem are analogous to results of Yang2019. Part (iii) implies that with asymptotic probability one, \[ U_n^{\alpha,\mathrm{raw}}(x) < \frac{p(x_o-) + p(x_o+)}{2} < L_n^{\alpha,\mathrm{raw}}(y) \] for $x < x_o - D \rho_n$, $y > x_o + D \rho_n$ and $D = 4 C^2 \{p(x_o+) - p(x_o-)\}^{-2}$. Thus, at points of discontinuity of $p$, the confidence band crosses a horizontal line on an interval of length $\mathcal{O}_p(\rho_n)$. Part (iv) demonstrates that our bounds are particularly accurate in regions where $p(x)$ is close to $0$ or $1$. Specifically, suppose that for some $\gamma > 0$, $p(x) = \mathcal{O}\{(x - a_o)^\gamma\}$ for $x \in [a_o,b_o]$. Then plugging in $y(x) = x + \rho_n^{1/(\gamma + 1)}$ reveals that \[ \mathbb{E}\{U_n^{\alpha,\mathrm{raw}}(x)\} \le D \{(x - a_o)^\gamma + \rho_n^{\gamma/(\gamma + 1)}\} , \quad x \in [a_o,b_o] , \] where $D = D(C_1,C_2,p)$. Analogously, if $1 - p(x) = \mathcal{O}\{(b_o - x)^\gamma\}$ for $x \in [a_o,b_o]$, then \[ \mathbb{E}\{1 - L_n^{\alpha,\mathrm{raw}}(x)\} \le D \{(b_o - x)^\gamma + \rho_n^{\gamma/(\gamma + 1)}\} , \quad x \in [a_o,b_o] . \] Presumably, the conclusions in part (iv) are not satisfied for the confidence band $(L_n^{\alpha,\mathrm{YB}},U_n^{\alpha,\mathrm{YB}})$.

Simulations

Here, we illustrate that our confidence bands have correct coverage in the sense of (ref) and are narrower than existing techniques. We consider both, the raw method in (ref) and (ref) and the non-crossing variant in (ref). Both methods are combined with the rounding technique in (ref) with $K=10^3$ in order to facilitate faster computation at a minimal cost in accuracy. For comparison, we use the bands of Yang2019 given in (ref) and (ref) with a minimal variance factor of $\sigma^2 = 1/4$ and the parametric bands of Nattino2014, implemented in the GivitiR package in the statistical software R R2022. Replication material for the simulations and applications is available under \href{https://github.com/marius-cp/replication_DDHPZ22}{https://github.com/marius-cp/replication_DDHPZ22}.

We use 1000 replications, a significance level of $\alpha = 0.05$ and simulate the covariates $X \sim \operatorname{U}[0,1]$. The binary outcomes are generated by $Y \sim \text{Bern}\{p_s(X)\}$ based on five distinct functional forms of the regression function $p_s(x)$ for $x \in [0,1]$ depending on a shape parameter $s \in \mathcal{S} := \{0, 0.1, \dots, 1\}$. The first four specifications of $p_s(x)$ satisfy the isotonicity assumption at (ref) and cover smooth, non-smooth as well as discontinuous setups. The last one contains non-isotonic functions $p_s(x)$ for $s > 0.5$. The choice $s=0$ results in the diagonal line $p_0(x) = x$ whereas the deviation from the diagonal increases with $s$. In particular, we consider the following specifications, which are illustrated in Figure (ref) for two exemplary shape values $s \in \{0.3, 0.7\}$.

Figure (ref) presents the average coverage rates for a range of sample sizes between 512 and 32\,768. We use the raw method for our bands in (ref) and (ref) as the raw bands are contained in the non-crossing ones. We find that, as predicted by the theory, our confidence bands have conservative coverage throughout all isotonic simulation setups and sample sizes. We observe coverage rates above 0.998 with the majority of 162 out of the 192 displayed coverage values being exactly one. The unreported non-crossing bands differ from the raw ones in less than one out of a hundred thousand instances. These deviations occur exclusively for large values of $s$ in the Step specification within constancy regions of the function $p$. As expected, our method as well as the bands of Yang2019 have incorrect coverage rates for the values $s > 0.5$ that violate isotonicity in the Wave specification when the sample size increases. The coverage rates of the Yang2019 bands are still larger as these are shown to be wider by Theorem (ref).

The parametric bands of Nattino2014 rarely achieve correct coverage rates unless in the cases $s=0$ and for the S-shaped regression functions. This can be explained as these bands are based on the assumption of a certain parametric form of $p_s(x)$, which is rarely satisfied. The results get worse for the non-smooth, the discontinuous and the non-isotonic specifications.

Figure (ref) displays the average widths of our and the Yang2019 bands. We present the theoretically wider non-crossing bands instead of the raw versions thereof. Their average widths is however non-distinguishable in these displays. We fix a medium degree of miscalibration $s=0.5$. The upper plot panel displays the widths averaged over all simulation runs and values $x \in [0,1]$ depending on the sample size $n$. We find that the size of both bands shrinks with $n$ and that we can reconfirm the ordering established in Theorem (ref). We further see that our bands are only narrow enough for practical use in large samples. The relative gain in width of our bands is the highest for large sample sizes, exactly for which we propose the application of our method for calibration validation. It is worth noting that the bands of Yang2019 are more generally valid than for the special case of binary observations.

The lower plot panel shows the widths averaged over the simulation replications, but depending on the values $x \in [0,1]$ for two selected sample sizes. It shows that the relative gains in width upon the bands of Yang2019 are particularly pronounced close to the edges of the unit interval. In applications to calibration assessment, these regions of predicted probabilities close to zero or one are often of the highest interest as for example in the subsequent section assessing the goodness-of-fit of low birth weight probability predictions.

Application: Predicting low birth weight probabilities

We apply our confidence bands to assess calibration of three binary regression specifications predicting the probability of a fetus having a low birth weight, defined as weighting less than 2500 grams at birth def_who. Recall that in the setting of calibration assessment, we call the function $p$ the calibration curve and our confidence bands are denoted as calibration bands. This follows the interpretation that for an event predicted with probability $x$, $p(x)$ denotes its true but unknown event probability. Perfectly calibrated predictions entail a calibration curve matching the diagonal line, $d(x) = x$. As the calibration band is a simultaneously valid confidence band for $p$, deviations of the calibration band from the diagonal line imply significantly miscalibrated predictions in this region.

We use U.S. Natality Data from the Natality2017, which provides demographic and health data for 3\,864\,754 births in the year 2017. For the data set at hand, a low birth weight is observed in 8.1% of the cases. We estimate three binary regression models by maximum likelihood on the same randomly drawn subset that contains all but 1\,000\,000 observations that we leave for external model validation. All three models contain standard risk factors such as the mother's age, body mass index and smoking behavior but they differ as follows. The first model uses a probit link function, and the explanatory variable week of gestation is categorized into four left-closed and right-open intervals with lower interval limits of 0, 28, 32 and 37 weeks, pertaining to the standard definitions of the World Health Organization of extremely, very, moderate and non preterm Quinn2016PretermDef. Through this categorization, the model specification can capture the week of gestation in a non-linear fashion. In contrast, the second model uses the week of gestation as a continuous explanatory variable and the third specification employs the cauchit instead of the probit link function, which is known to produce less confident predictions close to zero and one Koenker2009. Additional details of the model specifications are given in the Supplementary Material.

The classical Hosmer-Lemeshow test rejects perfect calibration of all three models with p-values of essentially zero for both, internal and external model validation, which leaves an applied researcher without any useful conclusions on model calibration. We show our calibration bands based on the non-crossing method with rounding to three digits, i.e., $K=10^3$ in (ref), with a confidence level of $1-\alpha = 95\%$ for the first model in Figure (ref) and for the other two models in Figure (ref). We constantly extrapolate the bands on the unit interval which preserves their theoretical coverage guarantees as discussed after Theorem (ref). Figure S3 in the Supplementary Material illustrates that the bands of Yang2019 are considerably wider in this application.

Recall that the validity of our bands relies on the isotonicity assumption of $p$, which we test for as detailed on in Appendix (ref). The test only rejects isotonicity at the $5\%$ level for the second model specification displayed on the left side of Figure (ref) with a crossing of the lower and upper bounds for probability predictions between $0.1\%$ and $2.7\%$. Hence, we can directly reject calibration for this model in the critical area of small predictions and furthermore, the remaining calibration band has to be interpreted carefully for this model. As the simulations in Appendix (ref) show that the isotonicity test can even detect slight violations of isotonicity with high power for much smaller sizes as considered in this application, type II test errors are barely a problem here and we can be confident about the isotonicity assumption for the other two model specifications.

For the first model, the calibration band encompasses the diagonal line for all forecast values, meaning that we cannot reject the null hypothesis of perfect calibration $p(x) = x$ at the $5\%$ level. More importantly, we are $95\%$ certain that the true calibration curve lies within the the band at any point $x \in [0,1]$, implying that we are confident that the model is at least as well calibrated as specified by the band. This is especially notable in the important region of predictions below $10\%$ in the magnified right panel of Figure (ref), where the calibration bands are remarkably close to the diagonal implying a particularly well calibrated model. E.g., we can conclude that for a prediction of $x=5\%$, a low birth weight occurs with a probability between $4.6\%$ and $6.7\%$.

In contrast, we reject calibration for both, the second and third model specifications as shown in Figure (ref). However, these bands are much more informative than a simple test rejection as they directly show the exact form of model miscalibration. For the second model specification, we can conclude that the predicted probabilities are particularly miscalibrated for the non-isotonic region discussed above and for values larger than $20\%$. The third specification entails miscalibrated probabilities for predictions below $10\%$ that are presumably of the highest importance for medical decision making. Finally notice that the wide bands for the third model specification between predicted probabilities of $5\%$ and $20\%$ are caused by little predictions in this interval.

Acknowledgement

T. Dimitriadis gratefully acknowledges financial support from the German Research Foundation (DFG) through grant number 502572912. A. Henzi and J. Ziegel gratefully acknowledge financial support from the Swiss National Science Foundation.

Supplementary material

The Supplementary Material further illustrates the rounding method in simulations, gives details on the low birth weight application and contains additional proofs.