The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
42,899 characters
Honest calibration assessment for binary outcome predictions
\title{Honest calibration assessment for binary outcome predictions}
\author[1]{Timo Dimitriadis}
\author[2]{Lutz D\"umbgen}
\author[3]{Alexander Henzi}
\author[4]{Marius Puke}
\author[2]{Johanna Ziegel}
\affil[1]{Heidelberg University, Germany}
\affil[2]{University of Bern, Switzerland}
\affil[3]{ETH Zürich, Switzerland}
\affil[4]{University of Hohenheim, Germany}
\affil[ ]{\linebreak
[email removed],
[email removed],
[email removed],
[email removed],
[email removed]
}
\maketitle
\begin{abstract}
Probability predictions from binary regressions or machine learning methods ought to be calibrated: If an event is predicted to occur with probability $x$, it should materialize with approximately that frequency, which means that the so-called calibration curve $p(\cdot)$ should equal the identity, $p(x) = x$ for all $x$ in the unit interval.
We propose honest calibration assessment based on novel confidence bands for the calibration curve, which are valid only subject to the natural assumption of isotonicity.
Besides testing the classical goodness-of-fit null hypothesis of perfect calibration, our bands facilitate inverted goodness-of-fit tests whose rejection allows for the sought-after conclusion of a sufficiently well specified model.
We show that our bands have a finite sample coverage guarantee, are narrower than existing approaches, and adapt to the local smoothness of the calibration curve $p$ and the local variance of the binary observations.
In an application to model predictions of an infant having a low birth weight, the bounds give informative insights on model calibration. \\[0.8em]
\noindent \textit{Keywords:}
Binary regression, calibration validation, isotonic regression, confidence band, goodness-of-fit, universally valid inference
\end{abstract}
\section{Introduction}
\label{sec: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}$,
\begin{equation}
\label{eq:isotonicity}
p(x) \ \le \ p(x'), \quad x \leq x'.
\end{equation}
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., \citet[Section 4.2]{Tutz2011} and \citet[Chapter 5]{Hosmer2013Book}.
It is still regularly carried out by the classical test of \cite{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 \citep{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 \citep{Nattino2020, Paul2013}, resulting in calls for a goodness-of-fit tests with inverted hypotheses \citep{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
\begin{equation}
\label{eq:coverage}
\mathrm{pr}\{L^\alpha \le p \le U^\alpha \ \text{on} \ \mathbb{R}\}
\geq 1 - \alpha.
\end{equation}
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{fig:ApplicationIntro} 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{sec:ApplicationLBW} 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 \texttt{R} \citep{R2022} is available under
\href{https://cran.r-project.org/package=calibrationband}{https://cran.r-project.org/package=calibrationband}.
\begin{figure}[tb]
\centering
\includegraphics[width=\textwidth]{Fig_1.pdf}
\caption{Left: Confidence band for the calibration curve (denoted calibration band in the application section) for the first model specification of the low birth weight application in Section \ref{sec:ApplicationLBW}. The blue band shows the confidence band based on the non-crossing method in \eqref{eq:NonCrossingBands} together with the rounding in \eqref{eqn:Rounding} with $K=10^3$, and the grey step function shows the isotonic regression estimate.
Right: Magnified version focusing on predicted probabilities below $10\%$.}
\label{fig:ApplicationIntro}
\end{figure}
Our confidence bands are valid in finite samples subject only to the mild monotonicity assumption at \eqref{eq:isotonicity}, implying that higher predictions entail a higher probability for $\{Y=1\}$, which is natural in the context of assessing calibration as already argued in \citet{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{app:test_iso}. 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 \citet{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 \citet{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 \citet{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, \cite{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 \citet{Yang2019} are valid, in a modified sense even in settings where the isotonicity assumption \eqref{eq:isotonicity} 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 \citep{Wright1981}.
Resampling schemes are theoretically found to be inconsistent for the isotonic regression \citep{Sen2010, GuntuboyinaSen2018}.
Other non-parametric approaches in the literature for constructing confidence bands for functions, many of them presented in the review by \citet{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 \citet{Stodden2016, Yu2020}.
\section{Construction of the confidence bands}
\label{sec:construction}
Within the regression setting, we construct confidence bands for the isotonic regression function $p$ by means of the classical confidence bounds of \cite{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
\begin{align*}
u^\delta(Z,m)
&= \max\{\xi \in [0,1] \colon \mathrm{pbin}(Z,m,\xi) \ge \delta\} \\
&= \begin{cases}
\mathrm{qbeta}(1-\delta,Z+1,m-Z),
& \ Z < m , \\
1,
& \ Z = m ,
\end{cases} \\
\ell^\delta(Z,m)
&= \min\{\xi \in [0,1] \colon \mathrm{pbin}(Z-1,m,\xi) \le 1 - \delta\} \\
&= \begin{cases}
\mathrm{qbeta}(\delta,Z,m+1-Z),
& \ Z > 0 , \\
0,
& \ Z = 0 .
\end{cases}
\end{align*}
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 \cite{Johnson2005}.
Assumption \eqref{eq:isotonicity} 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{lem:hoeffding},
\begin{equation}
\label{eq:single.bounds}
\mathrm{pr}\{p_j \le u^\delta(Z_{jk},n_{jk})\} \ge 1 - \delta,
\qquad
\mathrm{pr}\{p_k \ge \ell^\delta(Z_{jk},n_{jk})\} \ge 1 - \delta .
\end{equation}
If we combine these bounds for all pairs $(j,k)$ in a given set $\mathcal{J}$ and use the assumption at \eqref{eq:isotonicity}, 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.
\begin{theorem}
\label{thm:Coverage}
For $x \in \mathbb{R}$, let
\begin{align}
\label{eq:upper_bound}
U^{\alpha, \mathrm{raw}}(x)
&= \inf_{(j,k) \in \mathcal{J}\colon x_j \ge x}
u^{\alpha/(N^2+N)}(Z_{jk},n_{jk}) , \\
\label{eq:lower_bound}
L^{\alpha, \mathrm{raw}}(x)
&= \sup_{(j,k) \in \mathcal{J}\colon x_k \le x}
\ell^{\alpha/(N^2+N)}(Z_{jk},n_{jk}) ,
\end{align}
where $\inf_\emptyset := 1$ and $\sup_\emptyset := 0$. If $p$ satisfies the isotonicity assumption at \eqref{eq:isotonicity}, then the resulting confidence band $(L^{\alpha,\mathrm{raw}},U^{\alpha,\mathrm{raw}})$ satisfies requirement \eqref{eq:coverage}.
\end{theorem}
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
\begin{align*}
U^{\alpha,\mathrm{raw}}(x) & = U^{\alpha,\mathrm{raw}}(x_i), \ \ x \in (x_{i-1},x_i] , \\
L^{\alpha,\mathrm{raw}}(x) & = L^{\alpha,\mathrm{raw}}(x_i), \ \ x \in [x_i,x_{i+1}) ,
\end{align*}
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{thm:Coverage} has two potential drawbacks. First, a natural nonparametric estimator for the function $p$ under the assumption \eqref{eq:isotonicity} 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]$ \citep{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 \eqref{eq:upper_bound} and \eqref{eq:lower_bound} 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:
\begin{align}
\label{eq:NonCrossingBands}
L^{\alpha,\mathrm{nc}} = \min(L^{\alpha,\mathrm{raw}}, \hat{p}),
\qquad
U^{\alpha,\mathrm{nc}} = \max(L^{\alpha,\mathrm{raw}}, \hat{p}) .
\end{align}
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 \eqref{eq:isotonicity}; see Section \ref{sec:Simulations} for details.
The potential crossing of the two bounds in Theorem \ref{thm:Coverage} also has an advantage. It allows for inference about the non-isotonicity of $p$, see Appendix \ref{app:test_iso}.
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
\begin{align}
\label{eqn:Rounding}
\widetilde{\mathcal{J}} = \{(j,k)\colon \{x_j\dots, x_k\} = \{x_1, \dots, x_n\} \cap [r/K, s/K] \text{ for some } r, s \in \mathbb{Z} \},
\end{align}
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 \eqref{eq:upper_bound} and \eqref{eq:lower_bound} 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 \eqref{eqn:Rounding} 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$.
\section{Relation to \cite{Yang2019}}
\label{sec:Relation_to_YB}
The methods of \cite{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
\begin{align}
\label{eq:upper_bound_YB}
U^{\alpha, \mathrm{YB}}(x)
& = \inf_{(j,k)\in \mathcal{J}\colon x_j \ge x}
\Bigl[ \frac{Z_{jk}^{\mathrm{iso}}}{n_{jk}}
+ \sqrt{\frac{\log\{(N^2 + N)/\alpha\}}{2 n_{jk}}} \Bigr] , \\
\label{eq:lower_bound_YB}
L^{\alpha, \mathrm{YB}}(x)
& = \sup_{(j,k)\in \mathcal{J}\colon x_k \le x}
\Bigl[ \frac{Z_{jk}^{\mathrm{iso}}}{n_{jk}}
- \sqrt{\frac{\log\{(N^2 + N)/\alpha\}}{2 n_{jk}}} \Bigr] .
\end{align}
This defines a confidence band $(L^{\alpha,\mathrm{YB}},U^{\alpha,\mathrm{YB}})$ with the following property:
\begin{equation}
\label{eq:coverage_YB}
\mathrm{pr} \{
L^{\alpha,\mathrm{YB}} \le \tilde{p} \le U^{\alpha,\mathrm{YB}}
\ \ \text{on} \ \ \mathbb{R} \}
\ge 1 - \alpha ,
\end{equation}
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 \eqref{eq:isotonicity} is violated. The proof of \eqref{eq:coverage_YB} follows from the arguments of \cite{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 \cite{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}})$.
\begin{theorem}
\label{thm:iso}
For all $\alpha \in (0, 1)$ and any data vector $\mathcal{Y} \in \{0,1\}^n$,
\[
L^{\alpha,\mathrm{YB}} \le L^{\alpha,\mathrm{nc}} \le L^{\alpha,\mathrm{raw}},
\quad
U^{\alpha,\mathrm{raw}} \le U^{\alpha,\mathrm{nc}} \le U^{\alpha,\mathrm{YB}}
\ \ \text{on} \ \ \mathbb{R}.
\]
\end{theorem}
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{app:test_iso}. In contrast, the bands by \citet{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 \eqref{eq:upper_bound_YB}, 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{prop:isobounds}. Likewise, for the computation of $L_i^{\alpha,\mathrm{YB}}$ in \eqref{eq:lower_bound_YB}, 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.
\begin{lemma}
\label{lem:isopart}
The cardinality of $\{\hat{p}(x_i)\colon i = 1, \dots, n\}$ is smaller than $3 n^{2/3}$.
\end{lemma}
\section{Theoretical properties of the confidence bands}
\label{sec:iso}
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 \cite{Yang2019} and consider the following assumption.
\begin{assumption}
\label{ass:A}
Let $\mathrm{Leb}(\cdot)$ denote Lebesgue measure, and let $W_n(B) = \#\{i\colon x_{ni} \in B\}$ for $B \subset \mathbb{R}$. There exist a non-degenerate interval $[a_o,b_o]$ and constants $C_1, C_2 > 0$ such that for sufficiently large $n$,
\[W_n(B) \ge C_1 n \mathrm{Leb}(B)\]
for arbitrary intervals $B \subset [a_o,b_o]$ such that $\mathrm{Leb}(B) \ge C_2 \log(n)/n$.
\end{assumption}
This assumption comprises the setting of \cite{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 \citet[Section~4.3]{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{ass:A} is satisfied almost surely, provided that $C_1, C_2 > 0$ are chosen appropriately.
\begin{theorem}
\label{thm:asymptotics}
Suppose that Assumption \ref{ass:A} is satisfied. Let $\rho_n = \log(n)/n$. There exist constants $C > 0$ depending only on $C_1,C_2$ with the following properties:
\noindent
(i) Suppose that $p$ is constant on $[a_o,b_o]$. With asymptotic probability one,
\begin{align*}
U_n^{\alpha,\mathrm{raw}}(x) &\le p(x) + C \sqrt{\rho_n/(b_o - x)},
\quad x \in [a_o,b_o) , \\
L_n^{\alpha,\mathrm{raw}}(x) &\ge p(x) - C \sqrt{\rho_n/(x - a_o)},
\quad x \in (a_o,b_o] .
\end{align*}
\noindent
(ii) Suppose that $p$ is Lipschitz-continuous on $[a_o,b_o]$ with Lipschitz constant $L > 0$. With asymptotic probability one,
\begin{align*}
U_n^{\alpha,\mathrm{raw}}(x) &\le p(x) + C (L\rho_n)^{1/3} ,
\quad x \in [a_o, b_o - \rho_n^{1/3} L^{-2/3}] , \\
L_n^{\alpha,\mathrm{raw}}(x) &\ge p(x) - C (L \rho_n)^{1/3} ,
\quad x \in [a_o + \rho_n^{1/3} L^{-2/3}, b_o] .
\end{align*}
\noindent
(iii) Suppose that $p$ is discontinuous at some point $x_o \in (a_o,b_o)$. With asymptotic probability one,
\begin{align*}
U_n^{\alpha,\mathrm{raw}}(x) &\le p(x_o-) + C \sqrt{\rho_n /(x_o - x)} ,
\quad x \in [a_o,x_o) , \\
L_n^{\alpha,\mathrm{raw}}(x) &\ge p(x_o+) - C \sqrt{\rho_n /(x - x_o)} ,
\quad x \in (x_o,b_o] .
\end{align*}
\noindent
(iv) Suppose that $\lim_{x \to a_o} p(x) = 0$. For sufficiently large $n$,
\[
\mathbb{E} \{U_n^{\alpha,\mathrm{raw}}(x)\}
\le C \inf_{y \in (x,b_o]} \{ p(y) + \rho_n/(y - x) \} ,
\quad x \in [a_o, b_o) .
\]
Analogously, if $\lim_{x \to b_o} p(x) = 1$, then for sufficiently large $n$,
\[
\mathbb{E} \{1 - L_n^{\alpha,\mathrm{raw}}(x)\}
\le C \inf_{y \in [a_o,x)} \{ 1 - p(y) + \rho_n/(x - y) \} ,
\quad x \in (a_o,b_o] .
\]
\end{theorem}
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 \citet[Sections 4.4 and 4.6]{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}})$.
\section{Simulations}
\label{sec:Simulations}
Here, we illustrate that our confidence bands have correct coverage in the sense of \eqref{eq:coverage} and are narrower than existing techniques.
We consider both, the raw method in \eqref{eq:upper_bound} and \eqref{eq:lower_bound} and the non-crossing variant in \eqref{eq:NonCrossingBands}.
Both methods are combined with the rounding technique in \eqref{eqn:Rounding} with $K=10^3$ in order to facilitate faster computation at a minimal cost in accuracy.
For comparison, we use the bands of \cite{Yang2019} given in \eqref{eq:upper_bound_YB} and \eqref{eq:lower_bound_YB} with a minimal variance factor of $\sigma^2 = 1/4$ and the parametric bands of \cite{Nattino2014}, implemented in the \texttt{GivitiR} package in the statistical software \texttt{R} \citep{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}.
\begin{figure}[tb]
\centering
\includegraphics[width=\textwidth]{Fig_2.pdf}
\caption{Illustration of the five simulated regression functions $p_s(\cdot)$, where the solid red line corresponds to the shape parameter value $s=0.3$ and the dashed blue line to $s=0.7$.}
\label{fig:p_alternatives}
\end{figure}
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 \eqref{eq:isotonicity} 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{fig:p_alternatives} for two exemplary shape values $s \in \{0.3, 0.7\}$.
\begin{enumerate}
\item
Monomial: \
First, we use the regression function $p_s(x) = x^{1-s}$, where $s \in \mathcal{S} \setminus \{1\}$.
This function is already used in the simulations in \citet[Appendix A]{DGJ_2021}.
\item
S-shaped: \
Second, the regression function follows an S-shaped form $p_s(x) = \left( 1+((1-x)/x)^{1+s} \right)^{-1}$, where $s \in \mathcal{S}$ pronounces the curves for larger values of $s$.
\item
Kink: \
Third, $p_s(x)$ linearly interpolates the points $(0,0), (0.2 + 0.8s, 0.2)$ and $(1,1)$ for $s \in \mathcal{S}$, resulting in a kink at the point $(0.2 + 0.8s, 0.2)$ for all $s > 0$.
\item
Step: \
Fourth, we use a step function with $s^\star \in \{5,6,\dots,14\}$ equidistant steps in the unit interval.
It is given by $p_s(x) = \big\{ \lfloor s^\star x \rfloor + \mathds{1}(x \not= 1) \big\} / s^\star$, where $s^\star = 15 - 10s$ and $s \in \mathcal{S} \setminus \{0\}$.
It doesn't nest the diagonal, but the deviation from it increases with $s$.
\item
Wave: \
Fifth, we use the cubic function $p_s(x) = 0.5 - (2s-1)(x-0.5) + 8s (x-0.5)^3$ that violates the isotonicity assumption in \eqref{eq:isotonicity} for any $s > 0.5$.
\end{enumerate}
Figure \ref{fig:coverage_unif} 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 \eqref{eq:upper_bound} and \eqref{eq:lower_bound} 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 \cite{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 \cite{Yang2019} bands are still larger as these are shown to be wider by Theorem \ref{thm:iso}.
\begin{figure}[tb]
\centering
\includegraphics[width=\textwidth]{Fig_3.pdf}
\caption{Empirical coverage rates of our confidence bands, the bands of \citet{Yang2019}, and the GiVitI bands for $1-\alpha = 0.95$, averaged over all covariate values for the five specifications of the regression function $p_s(\cdot)$, different shape values $s$ and a range of sample sizes $n$.
For our bands, we use the raw method in \eqref{eq:upper_bound} and \eqref{eq:lower_bound}, together with rounding in \eqref{eqn:Rounding} with $K=10^3$.
The choices $s=1$ in the Monomial, and $s=0$ in Step specification are not defined.}
\label{fig:coverage_unif}
\end{figure}
The parametric bands of \cite{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.
\begin{figure}[tb]
\centering
\includegraphics[width=\textwidth]{Fig_4_top.pdf}
\includegraphics[width=\textwidth]{Fig_4_bottom.pdf}
\caption{Top: Average widths of the $95\%$ confidence bands by sample size for each of the five specifications of $p_s(x)$ given in the main text for a fixed value $s=0.5$.
Bottom: Average widths by covariate value $x$ for two sample sizes.
In both panels, the solid red line corresponds to our bands based on the non-crossing method in \eqref{eq:NonCrossingBands} together with rounding in \eqref{eqn:Rounding} with $K=10^3$, and the dashed blue line corresponds to the \cite{Yang2019} bands.}
\label{fig:width_unif}
\end{figure}
Figure \ref{fig:width_unif} displays the average widths of our and the \cite{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{thm:iso}.
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 \cite{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 \cite{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.
\section{Application: Predicting low birth weight probabilities}
\label{sec:ApplicationLBW}
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 \citep{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.
\begin{figure}[tb]
\centering
\includegraphics[width=\textwidth]{Fig_5.pdf}
\caption{Calibration bands for the second model specification on the left and for the third specification on the right for the low birth weight application.
The blue band denotes the calibration band based on the non-crossing method in \eqref{eq:NonCrossingBands} together with the rounding in \eqref{eqn:Rounding} with $K=10^3$, and the grey step function shows the isotonic regression estimate. The diagonal line is given in red color whenever it is not contained in the calibration band.}
\label{fig:appl_lbw}
\end{figure}
We use U.S.~Natality Data from the \citet{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 \citep{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 \citep{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 \eqref{eqn:Rounding}, with a confidence level of $1-\alpha = 95\%$ for the first model in Figure \ref{fig:ApplicationIntro} and for the other two models in Figure \ref{fig:appl_lbw}.
We constantly extrapolate the bands on the unit interval which preserves their theoretical coverage guarantees as discussed after Theorem \ref{thm:Coverage}.
Figure S3 in the Supplementary Material illustrates that the bands of \citet{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{app:test_iso}.
The test only rejects isotonicity at the $5\%$ level for the second model specification displayed on the left side of Figure \ref{fig:appl_lbw} 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{app:test_iso} 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{fig:ApplicationIntro}, 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{fig:appl_lbw}.
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.
\section*{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.
\section*{Supplementary material}
The Supplementary Material further illustrates the rounding method in simulations, gives details on the low birth weight application and contains additional proofs.