Smooth Calibration Curves With Uncertainty

Why Smooth Calibration Curves?

A binary prediction \(p\) is calibrated when, among cases assigned that prediction, the event occurs with probability \(p\). Binned calibration tables summarize this relationship in a few groups, but their appearance can depend strongly on where the bin boundaries fall. A smooth calibration curve instead estimates the event probability across the prediction range without imposing a small set of bins. The curve is useful only when its uncertainty and the amount of data supporting each part of it are visible too.

Let \(Y_i\) indicate whether the event occurred for observation \(i\), and let \(p_i\) be its predicted risk, with \(Y_i\in\{0,1\}\) and \(0<p_i<1\). The predictions should be evaluated on observations not used to fit the prediction model, or otherwise obtained in a way that avoids evaluating the model on its own training outcomes. Define the calibration function

m(p)=P(Y=1\mid \text{predicted risk}=p).

Perfect calibration corresponds to \(m(p)=p\). A curve above the diagonal has observed event probabilities greater than the predictions at those risks; a curve below it has lower probabilities. In practice, \(m\) is unknown and must be estimated from the observed pairs \((p_i,Y_i)\).

Logistic Recalibration as a Smooth Curve

One parametric way to estimate the curve is logistic recalibration. Write \(\operatorname{logit}(p)=\log\{p/(1-p)\}\) and \(\operatorname{expit}(z)=1/(1+e^{-z})\). The recalibration model is

\operatorname{logit}\{m(p)\}=\alpha+\beta\operatorname{logit}(p), \qquad \widehat m(p)=\operatorname{expit}\{\hat\alpha+\hat\beta\operatorname{logit}(p)\}.

Here \(\alpha\) shifts the curve on the log-odds scale, while \(\beta\) changes its slope against the predicted log odds. Perfect calibration is the special case \(\alpha=0,\ \beta=1\). Unlike a flexible smoother, this model restricts the entire curve to a two-parameter family. Its smooth appearance does not establish that this shape is adequate.

Suppose the fitted coefficient vector \((\hat\alpha,\hat\beta)^\mathsf{T}\) has estimated covariance matrix \(V\). At a chosen risk \(p_0\), put \(x_0=\operatorname{logit}(p_0)\) and \(g=(1,x_0)^\mathsf{T}\). The fitted log odds are \(g^\mathsf{T}(\hat\alpha,\hat\beta)^\mathsf{T}\), and their estimated variance is \(g^\mathsf{T}Vg\). A pointwise Wald interval on the log-odds scale can be transformed by the expit function to obtain an interval for \(m(p_0)\).

Worked Example: A Pointwise Interval From Logistic Recalibration

Suppose the fitted coefficients are \(\hat\alpha=-0.20\) and \(\hat\beta=0.80\), and their estimated covariance matrix is

V= \begin{pmatrix} 0.04 & -0.01\\ -0.01 & 0.09 \end{pmatrix}.

Find a pointwise 95% Wald interval for the calibrated event probability at predicted risk \(p_0=1/4\).

First, \(x_0=\log\{(1/4)/(3/4)\}=\log(1/3)\approx-1.0986\). The fitted log odds at this risk are

\hat\eta=-0.20+0.80(-1.0986)\approx-1.0789.

The estimated variance includes the covariance contribution:

\widehat{\operatorname{Var}}(\hat\eta) =0.04+2(-1.0986)(-0.01)+(-1.0986)^2(0.09) \approx0.17060.

Thus the standard error is \(\sqrt{0.17060}\approx0.4130\). The 95% log-odds limits are \(-1.0789\pm1.96(0.4130)\), or approximately \((-1.8884,-0.2693)\). Applying expit to each endpoint gives

\bigl(\operatorname{expit}(-1.8884),\operatorname{expit}(-0.2693)\bigr) \approx(0.1314,0.4331).

The interval is pointwise: it concerns the curve at this one specified risk, under the fitted model and its usual large-sample Wald approximation. It does not account for possible misspecification of the logistic recalibration shape.

Nonparametric Smoothing

A nonparametric smoother estimates the curve without requiring it to be linear on the logit scale. For a local-constant kernel estimate at \(p_0\), choose a bandwidth \(h>0\) and a nonnegative kernel function \(K\). Define weights \(w_i=K\{(p_i-p_0)/h\}\). When their sum is positive, the estimate is

\hat m(p_0)=\frac{\sum_{i=1}^n w_iY_i}{\sum_{i=1}^n w_i}.

Observations with predictions near \(p_0\) receive more weight, according to the chosen kernel and bandwidth. A small bandwidth gives a more local but more variable estimate; a large bandwidth gives a smoother, less variable curve that can conceal local departures from calibration. The bandwidth is therefore part of the analysis, not a cosmetic plotting choice.

Conditional on the prediction values, if outcomes are independent Bernoulli variables with probabilities \(m(p_i)\), the variance of this weighted average is

\operatorname{Var}\{\hat m(p_0)\mid p_1,\ldots,p_n\} =\frac{\sum_i w_i^2m(p_i)\{1-m(p_i)\}}{(\sum_i w_i)^2}.

A common approximate standard error substitutes a fitted local probability, such as \(\hat m(p_0)\), for each \(m(p_i)\) in this expression. This is a local approximation: if the true probabilities vary substantially across the window, treating them as equal can misstate uncertainty. Also, this standard error describes sampling variability around the smoother's local target; it does not include smoothing bias.

Worked Example: A Kernel Estimate and Its Approximate Standard Error

Consider a triangular kernel \(K(u)=1-|u|\) for \(|u|<1\), and \(K(u)=0\) otherwise. At \(p_0=0.45\), use bandwidth \(h=0.20\). Suppose the observations with nonzero weights are \((p_i,Y_i)=(0.35,1),(0.45,0),(0.55,1)\).

The standardized distances are \(-0.5,0,0.5\), so the weights are \(0.5,1,0.5\). Their sum is \(2\). Therefore

\hat m(0.45)=\frac{0.5(1)+1(0)+0.5(1)}{2}=0.50.

Using the local approximation \(m(p_i)\approx\hat m(0.45)=0.50\), each Bernoulli variance is approximated by \(0.50(1-0.50)=0.25\). The estimated variance is

\widehat{\operatorname{Var}}\{\hat m(0.45)\} \approx \frac{(0.5^2+1^2+0.5^2)(0.25)}{2^2} =\frac{1.5(0.25)}{4} =0.09375.

The approximate standard error is \(\sqrt{0.09375}\approx0.3062\). This is large relative to the estimate, reflecting the limited effective information in these three weighted observations. A normal interval based on this standard error may be poor with such a small local sample; the calculation is a useful uncertainty diagnostic, not a guarantee of accurate coverage.

Why the Boundary Needs Special Care

Near the minimum or maximum predicted risk, a symmetric smoothing window is truncated: there are observations on one side of the target but few or none on the other. A local-constant estimate then averages outcomes over an asymmetric window. Its target is generally not \(m(p_0)\), even with a very large sample.

The population target can be derived directly. Suppose predictions have density \(f\), and use a uniform kernel with bandwidth \(h\) at the lower boundary \(p_0=0\). The window is \([0,h]\). The expected local-constant target is the prediction-density-weighted average

\frac{\int_0^h m(p)f(p)\,dp}{\int_0^h f(p)\,dp},

provided the denominator is positive. It is not generally the unweighted average of \(m(p)\) across the window. For example, if \(m(p)=p\) and predictions are uniform on \([0,1]\), the target is \(h/2\). With a nonuniform prediction density, that value changes.

Worked Example: A Boundary Target Depends on Prediction Density

Let \(m(p)=p\), and consider the local-constant uniform-kernel target at \(p_0=0\) with bandwidth \(h\), where \(0<h\leq1\). First suppose predictions are uniform on \([0,1]\), so \(f(p)=1\). Then

\frac{\int_0^h p\,dp}{\int_0^h 1\,dp} =\frac{h^2/2}{h} =\frac h2.

Now suppose instead that prediction density is \(f(p)=2p\) on \([0,1]\). This is a valid density because \(\int_0^1 2p\,dp=1\). The target becomes

\frac{\int_0^h p(2p)\,dp}{\int_0^h 2p\,dp} =\frac{2h^3/3}{h^2} =\frac{2h}{3}.

Thus the value \(h/2\) requires the uniform-density assumption in this example. Under the second density the target is \(2h/3\), despite the same calibration function and the same window. Both targets differ from \(m(0)=0\) for positive \(h\): this is boundary smoothing bias.

One alternative is local-linear smoothing. At \(p_0\), fit a weighted straight line \(a+b(p_i-p_0)\) to the outcomes, minimizing the weighted sum of squared residuals, and report the fitted intercept \(a\). Its boundary advantage can be seen from an exact reproduction property.

Result: If the local-linear weighted least-squares problem has a unique solution, its fitted value at \(p_0\) equals \(m(p_0)\) whenever the observed conditional means satisfy \(m(p_i)=c+d(p_i-p_0)\) for constants \(c,d\). This holds whether \(p_0\) is in the interior or at a boundary.

Write the weighted objective.
Under the stated linear-mean condition, the objective is \(\sum_i w_i\{Y_i-a-b(p_i-p_0)\}^2\). For the reproduction claim, consider its conditional mean version, replacing \(Y_i\) by \(m(p_i)=c+d(p_i-p_0)\).

Identify a zero-residual fit.
Setting \(a=c\) and \(b=d\) makes every residual \(m(p_i)-a-b(p_i-p_0)\) equal zero, so the objective is zero.

Use uniqueness.
The objective is a sum of nonnegative weighted squares. A zero value is its minimum, and uniqueness of the weighted least-squares solution means the fitted coefficients must be \(a=c,b=d\). The fitted value at \(p_0\) is its intercept \(a=c=m(p_0)\). \(\square\)

This result concerns the conditional mean fit under an exactly linear local relationship; it does not say that local-linear estimates are unbiased for every curved calibration function. It explains why local-linear methods can reduce the leading boundary bias when the function is smooth, but they may have higher variance and can sometimes produce estimates below zero or above one. The local-constant estimate, by contrast, is a weighted average of binary outcomes and therefore lies in \([0,1]\) whenever the weights are nonnegative.

Uncertainty, Interpretation, and Reporting

A pointwise uncertainty interval describes uncertainty at a single chosen prediction value. For a logistic recalibration curve, the coefficient covariance matrix and a delta-method or transformed Wald interval provide a standard approach. For a nonparametric curve, a pointwise standard error can be calculated from the local weighted observations, but its reliability depends on local sample size, bandwidth, and the adequacy of the variance approximation. Neither pointwise method automatically gives a simultaneous band that covers the entire curve with a specified probability.

The uncertainty band should be interpreted alongside the prediction distribution. Where few observations have predictions, the curve and its interval are unstable; a narrow-looking line from an overly smooth fit does not create information. It is often sensible to restrict the displayed curve to the range with adequate data support rather than imply certainty in sparsely populated tails. A useful plot shows the identity line, the estimated smooth curve, an uncertainty display, and some indication of where predictions are concentrated.

For kernel smoothers, distinguish three issues: variability from observing a finite sample, bias from smoothing across different risks, and boundary bias from a one-sided neighborhood. Increasing the bandwidth generally lowers variability but can increase smoothing bias. A confidence ribbon based only on a standard error accounts for variability, not necessarily either kind of bias. State the smoothing method and bandwidth, explain how uncertainty was obtained, and avoid interpreting a boundary estimate as though it were based on a balanced neighborhood.

Common Mistakes

Assuming a smooth line is automatically reliable. Smoothness may be imposed by the method; sparse data can still make the estimate uncertain.
Calling a local average an unweighted average over risks. Its population target weights by the prediction density as well as the kernel.
Using a symmetric-window intuition at an endpoint. Boundary neighborhoods are one-sided, and local-constant estimates can be biased even in large samples.
Treating a pointwise interval as a whole-curve guarantee. Coverage at each fixed risk does not imply simultaneous coverage across all risks.
Ignoring smoothing bias when displaying a standard-error ribbon. A narrow variability interval does not rule out bias from the chosen bandwidth or boundary behavior.

Key takeaway: A smooth calibration curve estimates the conditional event probability across predicted risks, but its shape and uncertainty depend on the method, bandwidth, data density, and boundary location. Show uncertainty and support, and be especially cautious at the edges.

Check Your Understanding

Try these questions without looking back.

What does the logistic recalibration model reduce to under perfect calibration?
Why does a local-constant smoother's population target depend on the density of predicted risks?
For the triangular-kernel example, explain why the variance calculation uses squared weights.
When does the local-linear reproduction result imply that the fitted value at \(p_0\) equals the conditional mean there?
What is the difference between a pointwise uncertainty interval and a simultaneous band for a curve?