Non-parametric inference on calibration of predicted risks

Calibration for a risk prediction model refers to the ability of the model to generate predicted probabilities that are close to their true counterparts. Calibration is a critical element of a risk model’s performance. It has been shown that the clinical utility of a risk prediction model is more sensitive to its calibration than to its discrimination [1]. Poorly calibrated predictions also provide incorrect information to patients and hinder informed shared decision-making. Despite this, the critical role of risk model calibration remains underappreciated, so much so that it is called the Achilles’ heel of predictive analytics [2].

For the prediction of binary outcomes, Van Calster et al. [3] proposed a hierarchy of calibration definitions. According to this hierarchy, mean calibration (aka calibration-in-the-large) refers to the closeness of the average of predicted and observed risks, and is often the very first step of evaluating a model [4]. Weak calibration refers to the intercept and slope of the ‘calibration line’, obtained by fitting a logit model associating the response with logit-transformed predicted risks, being equal to, respectively, zero and one. The corresponding likelihood ratio (LR) test enables inference on weak calibration. Moderate calibration refers to the average risk among all individuals with a given predicted risk being equal to that predicted risk. Finally, strong calibration requires that the predicted and observed risks are equal within each subgroup defined by all distinct covariate patterns. The authors argue that strong calibration is neither achievable nor strictly desirable, and that moderate calibration should be the ultimate metric of merit for risk prediction models.

Moderate calibration is often visually assessed via a flexible calibration plot (henceforth referred to as calibration plot for brevity) [3], which is the plot of the conditional mean of observed risk (y-axis) as a function of predicted risk (x-axis). When the response is binary, in the common situation that predicted risks are continuous with no or few ties, calculating this conditional mean requires grouping the data into bins or applying local smoothing methods. Moderate calibration is then assessed by the closeness of this conditional mean to the line of perfect calibration (identity line). Scalar indices have been proposed for summarizing this discrepancy. Examples include Emax, E50, Estimated Calibration Index [5] and Integrated Calibration Index (ICI) [6, 7]. As such metrics require estimating the conditional mean function, their calculation similarly requires a tuning parameter (e.g., loess bandwidth).

In addition to visual assessment, answering the question “is this model moderately calibrated in this population?” can be approached as a formal inference problem. The historical approach towards inference on moderate calibration for binary outcomes has been the Hosmer-Lemeshow (HL) test. This test is based on grouping predicted risks into quantiles (conventionally deciles) and comparing predicted and observed event frequencies, with the resulting statistic approximated as chi-square [8]. The HL test has been criticized due to its sensitivity to the arbitrary specification of the number of groups and lack of information it provides on the direction of miscalibration [9]. An alternative test for moderate calibration has been proposed that does not require tuning parameters, but it is simulation-based and thus is subject to Monte Carlo error [10].

Hosmer et al.[9] reviewed several alternative goodness-of-fit (GOF) tests for binary responses. However, none of the alternative methods examined seems to be a test for moderate calibration. They defined the purpose of GOF testing for logistic regression models to be three-fold: to examine whether the logit transformation is the correct link function, whether the linear predictor is of the correct form (i.e., there is no need for the inclusion of additional covariates, transformations, or interactions), and whether the variance of the response variable conditional on covariates is that of Bernoulli distribution [9]. Similarly, Allison [11] considered the objective of GOF testing to be to examine whether a more complicated model is required. These are fundamentally different pursuits than examining moderate calibration.

Given this, classical methods do not provide an opportunity for an objective and unequivocal assessment of moderate calibration. In this work, we discuss recently proposed, and propose novel, methods for inference on moderate calibration and associated metrics for quantifying the degree of miscalibration, that do not involve regularization. Our novel method is motivated by standard practice in predictive analytics which involves evaluating mean calibration before, or in tandem with, moderate calibration [4]. As mean calibration is a necessary condition for moderate calibration, proper inference on moderate calibration should incorporate the results of mean calibration evaluation. As we demonstrate, such a two-pronged approach also increases statistical power.

The rest of this manuscript is structured as follows: after describing the context, we briefly review recently proposed graphical and statistical methods for moderate calibration based on partial sums of prediction errors. Motivated by these developments, we construct a stochastic process that converges asymptotically to standard Brownian motion, enabling us to take advantage of known properties of Brownian motion to propose joint inference on mean and moderate calibration. We conduct simulation studies comparing the performance of these tests. We showcase these methods in a case study in prediction of mortality after a heart attack, where we propose graphical illustration of results and interpretation of findings. We conclude by discussing how the proposed test can be added to the toolbox for examining model calibration, and suggesting future research in this area.

We conducted a series of simulation studies to evaluate the finite-sample behavior of the Brownian motion and bridge tests under the null hypothesis, and their comparative power in detecting various forms of model miscalibration. All numerical results in this work are produced in R [23]. Our implementation of $F()$ is provided in the cumulcalib package (function pMAD_BM)[24]. For the Kolmogorov distribution ($G()$), we used the implementation in the CPAT package [25], but one can trick any implementation of the asymptotic Kolmogorov-Smirnov test to obtain the desired results (an example is provided in Appendix B).

The first set of simulations evaluated the finite-sample properties of both tests under $H_{0}$. This was performed to assess at what effective sample size the null behavior of these tests is substantially different from their asymptotic behavior. These simulations were performed by generating predicted risks as $\mbox{logit}(\pi)=\beta_{0}+X$, with $X\sim Normal(0,1)$, and responses $Y\sim Bernoulli(\pi)$. We varied $\beta_{0}$ in $\{-2,-1,0\}$ (resulting in average event probabilities of, respectively, 0.155, 0.303, and 0.500) and sample size in $\{50,100,250,1000\}$ in a fully factorial design. Each scenario was run 10,000 times. Empirical CDF plots of p-values are provided in Figure 1. As p-values should follow the standard uniform distribution under $H_{0}$, empirical CDFs should be close to the identity line. The proportion of simulated p-values<0.05 is also reported on each panel. Generally, both tests behaved similarly. In small samples, the p-values seem to be slightly biased upwards, so the corresponding tests are conservative. The proportion of p-values<0.05 was within 10% of the desired proportion (i.e., 0.045-0.055) only in the bottom half of the panels ($n\geq 250$). Taking this as a subjective criterion, the total variance $\sum_{i=1}^{n}\pi_{i}(1-\pi_{i})$ (as a measure proportional to the effective sample size) should be $\geq 30$ for the asymptotic behavior to be acceptable. This corresponds to 40, 51, and 72 events for $\beta_{0}$ of -2, -1, and 0, respectively. However, we leave deeper investigation of the small-sample properties of the tests to future studies and also note that for the Kolmogorov distribution small-sample corrections are available [21, 26], although the appropriate form of such a correction is not obvious for the current context. No such correction was used here to provide more comparable results between the two tests. In general, if any doubt exists about the adequacy of the sample size, one can carry out a simulation-based version of these tests where the null distribution of $S^{*}$ and $B^{*}$ (and $S_{n}$) is obtained by simulating many vectors of responses in each $Y_{i}\sim Bernoulli(\pi_{i})$.

Next, we evaluated the power of the tests in detecting miscalibration. Following recommendations on objectively deciding on the number of simulations [27], we obtained the results through 2,500 Monte Carlo samples such that the maximum S.E. of the estimated probability of rejecting $H_{0}$ would be 0.01. In these simulations, we modeled a single predictor $X\sim Normal(0,1)$ and the true risk as $p=1/(1+\mbox{exp}(-X))$. We evaluated the performance of the tests in simulated independent samples of $n$ observations when the predicted risks suffered from various degrees of miscalibration.

Two sets of simulations were performed. In the first set, we assumed the risk model generated potentially miscalibrated predictions in the form of $\mbox{logit}(\pi)=a+b\,\mbox{logit}(p)=a+bX$. Given the linear association on the logit scale between predicted and actual risks, weak and moderate calibration are equivalent in these scenarios and, therefore, the LR test (simultaneously testing whether $a=0$ and $b=1$), which is otherwise a test for weak calibration, has the maximum theoretical power in detecting miscalibration. As such, this simple setup provides an opportunity to judge the performance of the tests against a gold standard. We simulated response values and predicted risks in a fully factorial design with values $a=\{-1/4,-1/8,0,1/8,1/4\}$ and $b=\{1/2,3/4,1,4/3,2\}$, and for three sample sizes: $n=\{100,250,1000\}$. Results are provided in Figure 2. All tests rejected the null hypothesis around 5% when the model was calibrated (the central panel). As expected, the LR test had the highest power in this set of simulations. In almost all scenarios, the bridge test had the highest power among the other tests. One exception was when $b=1$ (middle column). In such instances, miscalibration is entirely in the mean calibration component ($a$), and thus the bridge component of the bridge test cannot be expected to contribute to the sensitivity of the assessment.

In the second set, the true risk model remained the same as above, and we modeled non-linear miscalibrations as $\mbox{logit}(\pi)=a+b\,\mbox{sign}(X)|X|^{1/b}$. Here, $a$ affects mean calibration, while the term involving $b$ is an odd function that flexibly changes the calibration slope. Similar to the previous setting, we simulated response values and predicted risks with values $a=\{-1/4,-1/8,0,1/8,1/4\}$, $b=\{1/2,3/4,1,4/3,2\}$, and $n=\{100,250,1000\}$, in a fully factorial design. Figure 3 presents the results. Here, the LR test could have low power due to non-linear miscalibrations (e.g., panels on the left side of the figure). The bridge test was more powerful than the Brownian motion test in almost all scenarios, again except those where the shape of the calibration plot was unchanged ($b=1$) and only the calibration intercept was affected ($a\neq 1$).

We used data from GUSTO-I, a clinical trial of multiple thrombolytic strategies for acute myocardial infarction (AMI), with the primary endpoint being 30-day mortality [28]. This dataset has frequently been used to study methodological aspects of developing or validating risk prediction models [29, 30, 31]. In line with previous studies, we used the non-US sample ($n$=17,796) of GUSTO-I to fit a prediction model for 30-day post-AMI mortality, and the US sample ($n$=23,034) for validating the prediction model [32]. We developed two models: one based on the entire development sample, and one based on a subset of 500 observations from that sample (to create a potentially miscalibrated model due to small development sample size). 30-day mortality was 7.2% in the full development sample, 8.2% in the small development sample, and 6.8% in the validation sample. All the analyses were conducted in R [23]. Ethics approval was not required because the anonymized data are publicly available for research. Our candidate risk prediction model is similar to those previously developed using this dataset [32]. The models were of the following form (coefficients are provided in the footnote of Figures 4 and 5)

$\mbox{logit}(\pi)=\beta_{0}+\beta_{1}[\text{age}]+\beta_{2}[\text{AMI location other}]+\beta_{3}[\text{AMI locaton anterior}]+\beta_{4}[\text{previous AMI history}]+\beta_{5}[\text{Killip score}]+\beta_{6}min[\text{blood pressure}],100)+\beta_{7}[\text{heart rate}]$

Figure 4 – panels (a) and (b) are calibration plots of the full model, which seem to be compatible with good calibration. Panels (c) and (d) present the graph of the $\{S_{i}\}$ process over time $t$, with additional information pertaining to the Brownian motion and the bridge tests, respectively. The corresponding $\pi$s are provided as the secondary x-axis on top. To indicate the size of the fluctuations, Tygert [15] suggested adding a triangle, centered at origin, with height equal to the standard deviation of the terminal value of the random walk being plotted as an aide for viewing the fluctuations of the process. We also find this useful here, even though as our random-walk is standardized, the height of the triangle will always be fixed at 1. In panel (c), the vertical (red) line indicates the location and value of the test statistic $S^{*}$. In addition, if hypothesis testing is intended, the critical value can be highlighted on the graph (the dashed red horizontal line, for a 5% significance level in this case). In our example, the line indicates that the null hypothesis (that the model is moderately calibrated) is not rejected at this level. $C^{*}$ was $0.0020$; the test statistic ($S^{*}$) was $1.2973$, corresponding to a p-value of $0.3889$. Panel (d) is our suggested presentation of the bridge test. Here, the gray ‘bridge’ line connecting the start and end of the random walk is provided, and $S_{n}$ and $B^{*}$ are drawn as vertical distances from this line (in, respectively, blue and red), along with the above-mentioned triangle and dashed lines of significance (at 5% for each component - note the line for $B^{*}$ is parallel to the bridge line). The mean calibration, $C_{n}$, was $-0.0016$; the corresponding Z-score ($S_{n}$), was $-1.0091$ (p=$0.3129$), while the bridge test statistic, $B^{*}$, had a value of $1.0284$ (p=$0.2407$). The bridge test’s unified p-value was $0.2701$, indicating no substantial miscalibration of this prediction model.

Figure 5 – panells(a) and (b) are the calibration plots of the small-sample model. The model seems to be making optimistic predictions: underestimating the risk among low/medium risk groups and overestimating it among high-risk groups. This is clearer in the partial sum plot (panels (c) and (d)), which shows an inverse U-shape that is typical of optimistic predictions. Again, the Brownian motion and bridge tests are shown on the graph in panels (c) and (d), respectively. $C^{*}$ was $0.0042$, corresponding to the test statistic ($S^{*}$) of $2.8151$ and a p-value of $0.0098$. Maximum cumulative prediction error occurred at $\pi^{*}=0.232$ ($t=0.752$). This is an estimation of the location ($\pi$) where the calibration function $\mathbb{E}(Y|\pi)$ crosses the identity line. Around this point, predictions are calibrated. For individuals with predicted risk above $\pi^{*}$, the model overestimates the risk by an average of $0.086$ (note the downward slope of the partial sum on the right side of the vertical red line in the middle panel). For others, it underestimated the risk by an average of $0.004$. Bridge test results are shown in panel (d). Mean calibration was $-0.0010$; the Z-score for mean calibration was $-0.6950$ (p=$0.4871$), while $B^{*}$ was $3.3381$ (p$<0.001$). The unified p-value was $<0.001$, clearly indicating miscalibration of this prediction model.

In this work, we built on recently proposed non-parametric methods for the assessment of moderate calibration for models for binary outcomes. We provided suggestions on visual and statistical assessment of model calibration on the partial sum domain, proposed a novel non-parametric test, and evaluated its performance against another recently proposed test. A major advantage of these methods is no requirements for binning or smoothing of observations. Our main innovation was a stochastic process representation of the partial sum of prediction errors yielding weak convergence to standard Brownian motion. This in turn enabled us to propose a unified ‘bridge’ test of moderate calibration based on two asymptotically independent test statistics. These developments are implemented as the R package cumulcalib[24].

Our simulations point towards the higher power of the (Brownian) bridge test over the recently proposed Brownian motion test. This is not surprising as the bridge test evaluates the behavior of the underlying stochastic process in two independent ways, thus having more opportunities for identifying miscalibration. However, the bridge test also has other interesting properties. Asymptotically, $B^{*}$ follows the Kolmogorov distribution, which is ubiquitously implemented as part of the Kolmogorov-Smirnov test. The joint inference on mean and moderate calibration is also of particular relevance in clinical prediction modeling, as oftentimes mean calibration is of independent relevance and is the very first step in examining model calibration [4]. As mean calibration is a necessary condition for moderate calibration, the joint inference helps control the overall error rate of inference on calibration. We used Fisher’s method for generating a unified p-value due to its optimality properties [22]. One can also use alternative methods for retaining an overall familywise error, such as Bonferroni’s correction for each component of the test.

Our focus in this work was to propose inferential methods and evaluate their performance. However, we do not suggest inattentive use of p-values for making true / false conclusions at an arbitrary significance level on whether a model has good calibration in a population. Inferential methods quantify the compatibility of data with the null hypothesis (moderate calibration in this case), and as such complement the existing visual and numerical assessments of miscalibration. As mentioned earlier, evaluation of calibration on the partial sum domain also results in ‘distance’ metrics (e.g., $C^{*}$ or ECCE-MAD as referred to by Arrieta-Ibarra et al. [16]), similar to metrics such as Emax and ICI for calibration plots, with convergence to a finite value that is affected by the degree of miscalibration. $C^{*}$ can thus be used in a similar fashion to ICI and Emax, with the added benefit that its estimation does not require any tuning parameters. However, there are situations where one can reasonably expect that the model is calibrated, and thus hypothesis-testing becomes particularity relevant. An example is when the target population is known to be identical or very close to the development population (e.g., another random sample from the same population, or different study sites participating in a clinical trial with a strict protocol). Another relevant context is after a model is recalibrated in a new sample. Here, it is justified to expect that the revised model can provide calibrated predictions, and $H_{0}$ has a plausible chance of being true. Strictly speaking, using the same data to revise the model and evaluate its performance will result in violation of the assumption underlying the proposed tests (that the predicted risks are fixed quantities). Nevertheless, as long as the degrees of freedom used in the model updating process is small (e.g., changing the intercept, rather than re-estimating coefficients [33]), the suggested inference procedures would still provide useful insight. Indeed, this is a typical approximation underlying many GOF tests that follow model fitting. Interestingly, this is one instance where the $B^{*}$ test can be used in isolation: maximum likelihood estimation of generalized linear models (e.g., logistic regression) guarantees that the average of predicted and observed risks will be equal. In such instances, $S_{n}=0$ by the estimation process, and thus $B^{*}$ can be used in isolation, to avoid loss of statistical power for testing another component of the null hypothesis. Rejection of the null hypothesis in this case would indicate the need to modify the linear component (e.g., adding interaction terms) or the link function. In this sense, the bridge test acts as a GOF test for logistic regression.

There are several areas for future research. The weak convergence of the random walk $\{S_{i}\}$ to Brownian motion opens the door for possibly more powerful tests. We identified two aspects of the stochastic process that provided independent opportunities for examining moderate calibration, but one can expand this logic towards comparing multiple aspects. Examples of possible attributes include minimum, maximum, range, or integrated square value of the scaled partial sum process. An important extension of this method would be for other types of outcome, including survival outcomes (as proposed, for example, for ICI [34]). To what extent small sample correction will improve the approximation of the Kolmogorov distribution to that of the bridge test statistic $B^{*}$ can be investigated (and similar small-sample corrections can be devised for the Brownian motion test). The finite sample performance of the distance metrics based on partial sum of prediction errors versus metrics based on the calibration plots should be investigated.

Non-parametric approaches based on partial sums of prediction errors lead to visual, numerical, and inferential tools that do not require any grouping or smoothing of data, and thus enable objective, unequivocal, and reproducible assessments. Simulations show that the proposed bridge test for moderate calibration has comparatively high power in detecting miscalibration. We recommend such assessments to accompany conventional tools (such as calibration plots and distance metrics) to provide a more complete picture of the calibration of a risk prediction model in a target population.

Let $e(z)=\mathbb{E}(Y|\pi=z)-z$ be the calibration error at predicted risk $z$. The model is moderately calibrated if $\forall z:e(z)=0$. Let $P(.)$ be the cumulative distribution function of the infinite population of predicted risks. The argument is for when $P(.)$ is strictly monotone in [0,1].

In the general case where the null hypothesis is not known to hold, the random walk constructed in the main text might have a non-zero expected value (drift). Indeed, as $n\rightarrow\infty$, at time $t$, $\frac{1}{\sqrt{n}}\mathbb{E}(S(t))\rightarrow\frac{1}{\sqrt{\tau}}\int_{0}^{\pi_{t}}e(x)dP(x)$, where $\tau=\int_{0}^{1}x(1-x)dP(x)$, and $\pi_{t}$ is the unique solution for $\pi$ in $\frac{1}{\tau}\int_{0}^{\pi}x(1-x)dP(x)=t$. Correspondingly, the bridged random walk will have an expected value of:

The second component of the bridge test ($H_{0B}$) demands that $S(t)-tS(1)$ converges to a Brownian bridge, thus requiring it to have zero expected value everywhere:

Given the 1:1 relation between $\pi_{t}$ and $t$ and the monotonicity of $P(.)$, we can map this equality to the $\pi$ domain:

We use birth weight data in the MASS package and an exemplary model for predicting low birth weight. Below, as a generic implementation we use the ks.test implementation of Kolmogorov CDF by creating a vector of n (=100) observation whose empirical CDF has distance B*/sqrt(n) from the CDF of the standard uniform distribution, such that when multiplied by sqrt(n) as part of the KS test it generates the desired metric. In R, one can request asymptotic test via ks.test() with exact=F argument. Of note, the /pKolmogorov() function in the cumulcalib package provides a direct implementation of this test.

library(MASS)
data("birthwt")
pi <- 1/(1+exp(-(2.15-0.050*birthwt$age-0.015*birthwt$lwt))) #The model

#Ordering predictions (and responses) from smallest to largest pi
o <- order(pi)
pi <- pi[o]
Y <- birthwt$low[o] #Outcome is low birth weight (1 v 0)

#Constructing the {S} process
T <- sum(pi*(1-pi))
t <- cumsum(pi*(1-pi))/T
S <- cumsum(Y-pi)/sqrt(T)

plot(t,S, type=’l’)

Sn <- S[length(S)] #This is Sn, the first component of the bridge test
pA <- 2*pnorm(-abs(Sn)) #p-value for first component of H0

Bstar <- max(abs(S-t*Sn)) #This is B*, the second component of the bridge test
d <- Bstar/10
X <- seq(from=0, to=1-d, length.out=100)
pB <- ks.test(X, punif, exact=F)$p.value #p-value for second component of H0

#Fisher’s method for generating a unified p-value
p <- 1-pchisq(-2*(log(pA)+log(pB)),4) #0.8381805