最终，我找到了用于创建校准图的算法的全面描述

J. Esarey, A. Pierce：“评估拟合质量和测试二元因变量模型中的错误指定。” 政治分析 20.4，第 480-500 页，2012

本文将其与基于分类的评估进行了比较。以下是这些想法的摘要以及我的评论和用于创建校准图的 R 代码。

当将模型预测的概率与“观察到的”概率进行比较时，存在没有观察到概率而只有零或一的问题，即（不）出现响应。这些值可以通过每个值的“邻域”中的距离加权平均来平滑到概率，例如使用 LOESS 局部回归。建立“邻域”的距离和权重可以在不同的空间进行测量。两个明显的可能选择是

1. 链接尺度上的距离，即在上，其中是第个观察的预测变量值，而是模型参数。$\eta_i=\beta_0 + \langle\vec{\beta},\vec{x}_i\rangle$$x_i$$i$$\beta$
2. 概率尺度上的距离，即$p_i=P(Y=1|\vec{x}_i) = 1 / (1+e^{-\eta_i})$

通过点或拟合的 LOESS将为每个产生一个估计量，可以将其与模型预测 进行比较： 有两个注意事项， 然而：$(y_i,\eta_i)$$(y_i,p_i)$$\hat{p}_i$$y_i$$p_i$

• 对于大于零的度数，LOESS 拟合可以产生 [0,1] 之外的值。由于这个原因，上述两个图中都缺少第一个值：它的估计概率为负。这可以通过切断零和一的概率来轻松纠正。$\hat{p}_i$
• LOESS 只考虑一定百分比（参数span）的邻居。

上面的图是使用默认span=0.75创建的。Esarey & Pierce 提出了两种不同的优化方法，并在脚注中链接到参考实现，但该链接同时停滞不前。因此，我实现了一个非常简单的优化标准：和之间的最小 MSE ，即。挑战者号航天飞机 O 形环数据集的结果可以在这里看到： 我还包括了$\hat{p}_i$$p_i$$\sum_i(\hat{p}_i - p_i)^2$$p_i$正如模型所预测的那样。Esarey & Pierce 还通过参数引导程序计算位于 80% 置信区间之外的值的百分比，但这可能更容易直接从的置信区间计算。这是在链接级别（右侧）生成校准图的代码：$p_i$

# Challenger Space Shuttle O-ring data:  ok vs temp
data <- data.frame(y=c(0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1,
                       0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1),
                   x=c(53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 
                       70, 70, 70, 70, 72, 73, 75, 75, 76, 76, 
                       78, 79, 81))
fit <- glm(y ~ x, data=data, family=binomial)

#
# calibration plot on link level
#
link.model <- predict(fit, data, type="link", se.fit=TRUE)
sort.key <- order(link.model$fit)
x <- link.model$fit[sort.key]

# prediction interval for probability
plot(link.model$fit, data$y, main="link level")
p.lower <- plogis(link.model$fit - qnorm(1-0.05/2) * 
            link.model$se.fit)[sort.key]
p.upper <- plogis(link.model$fit + qnorm(1-0.05/2) * 
            link.model$se.fit)[sort.key]
polygon(c(x,rev(x)), c(p.lower, rev(p.upper)), col="#dddddd", 
          border=NA)
points(link.model$fit, data$y) # replot overplotted points
lines(x, plogis(x), col="red")

# LOESS fit
optim.span <- optimize(resub.mse, c(0.1,1.0),
                       y=data$y, x=link.model$fit, 
                       p.model=p.model)
span <- optim.span$minimum
p.fit <- loess(y ~ x, data=data.frame(y=data$y, 
          x=link.model$fit), family="gaussian", degree=1, 
                span=span)
p.cutfit <- predict(p.fit, data.frame(x=x))
p.cutfit[p.cutfit < 0] <- 0
p.cutfit[p.cutfit > 1] <- 1
lines(x, p.cutfit)

legend("topleft", c("model", sprintf("LOESS (span=%4.2f)", 
       span)),
       col=c("red","black"), lty=1)

# the optimization function for estimating span
resub.mse <- function(span, y, x, p.model) {
  fit <- loess(y ~ x, family="gaussian", degree=1, span=span)
  return(sum((fit$fitted - p.model)^2))
}

另一种方法，显然没有在文献中讨论，是由 R 函数开箱即用提供的条件密度图cdplot。

条件密度图直接估计的任意数量的水平非参数，而不假设统计模型。在逻辑回归的情况下，只有两个级别（0 和 1）并且回归拟合的参数模型。因此可以直接比较这两个估计量，以查看逻辑模型是否与数据匹配。$P(Y=\omega_i|x)$$\omega_i$$P(Y=1|x)$

cdplot通过贝叶斯定理 估计 其中表示概率密度，由数据中的核密度估计器估计。这个估计中唯一棘手的部分是和的估计器必须使用相同的内核带宽。与 LOESS 方法相比，这有两个概念优势：$P(Y=1|x)$

• 结果保证产生概率，并且不会发生（如对于 LOESS）该值超出范围$[0,1]$
• 它不需要将级别的数字解释为 0 和 1，以便在数字上有意义或完全适用。

与 LOESS 方法一样，预测变量必须是一个标量值，在完全类比中，可以使用内核密度估计器需要选择一个带宽，在文献中通常推荐使用“插件方法”（在 R 函数中）（以及在 的文档中，尽管它使用不同的默认值）。$\eta$bw="SJ"cdplotdensity

为了比较，我实现了一个额外的带宽选择方法，该方法选择使 cdplot 最接近逻辑预测的带宽。这可以作为关于逻辑模型最好的说法的基线；-)

为了更好的易读性，这里省略了逻辑模型中 95% 置信带图的代码：

# Space Shuttle Challenger temp vs oring-ok
data <- data.frame(y=c(0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1,
                       0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1),
                   x=c(53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70, 70,
                       70, 70, 72, 73, 75, 75, 76, 76, 78, 79, 81))
fit <- glm(y ~ x, data=data, family=binomial)

# helper function for finding the bandwidth
# that is closest to the logistic model
resub.mse <- function(bw, y, x, p.model) {
  cdfit <- cdplot(x, y, bw=bw, plot=FALSE)
  return(sum((cdfit[[levels(y.factor)[1]]](x) - p.model)^2))
}

#
# logistic prediction vs. link
#
link.model <- predict(fit, data, type="link", se.fit=TRUE)
p.model <- plogis(link.model$fit)
sort.key <- order(link.model$fit)
x <- link.model$fit[sort.key]
plot(link.model$fit, data$y, main="link level")
lines(x, plogis(x), col="red")

# cdplot vs. link
# note that we must code the level of interest
# as FIRST level (for cdplot)
y.factor <- factor(data$y, levels=c(1,0))
optim.bw <- optimize(resub.mse, c(bw.nrd0(x)/10, (max(x)-min(x))/2),
                     y=y.factor, x=link.model$fit, p.model=p.model)
bw <- optim.bw$minimum
p.kernel <- cdplot(link.model$fit, y.factor, bw="SJ", plot=FALSE)
lines(x, p.kernel$'1'(x))
p.kernel <- cdplot(link.model$fit, y.factor, bw=bw, plot=FALSE)
lines(x, p.kernel$'1'(x), col="blue", lty=2)

legend("topleft",
       c("model", "cdplot (bw='SJ')", sprintf("closest cdplot (bw=%4.2f)", bw)),
       col=c("red", "black", "blue"), lty=c(1,1,2))

校准

为了完整起见，这里有另外两种生成校准图的方法：第一种是使用 Nattino 等人介绍的校准带。(2014)，纳蒂诺等人。(2016)和 Nattino 等人。(2017)。简而言之，他们使用要评估的模型的预测概率将阶多项式逻辑函数（使用参数是使用由似然比统计控制的标准前向选择程序来选择的，该统计量说明了用于选择的前向过程。校准带可用于内部和外部校准。该过程在Stata中实现（$^{1}$$^{2}$$^{3}$$m$$m\geq 2$$m$$m$calibrationbelt) 和 R (包givitiR)。以下是使用挑战者数据的示例：

# Challenger Shuttle Challenger temp vs oring-ok
dat <- data.frame(y=c(0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1,
                       0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1),
                   x=c(53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70, 70,
                       70, 70, 72, 73, 75, 75, 76, 76, 78, 79, 81))

fit <- glm(y ~ x, data=dat, family=binomial)

preds_p <- predict(fit, type = "response")

cb <- givitiCalibrationBelt(o = dat$y, e = preds_p , devel = "internal", confLevels = c(0.95, 0.8))

plot(cb, main = "", las = 1, ylab = "Observed probabilities", xlab = "Model probabilities", table = FALSE)

身份线显示为红色。浅灰色区域是 80% 置信校准区间，而深灰色区域是 95% 置信区间。理想情况下，红线在整个概率范围内都在腰带内。在示例中，置信区间很大：校准带在校准方面显示出很大的不确定性。由于红线位于区间内，我们不能拒绝良好校准模型的假设。

为了更好地说明校准带，让我们看一下校准良好模型的校准带：

在这里，置信区间要窄得多。因为红色标识线在整个范围内都位于带内，所以它几乎没有提供错误校准的证据。

作为其他一些答案，在 R 包中实现的第二种方法rms依赖于拟合预测和观察概率的非参数平滑器。它还绘制了基于 bootstrap 的偏差校正估计。详细信息可以在 Harrel (2015)中找到。$^{4}$

library(rms)

mod <- lrm(y~x, dat = dat, x = TRUE, y = TRUE)

res <- calibrate(mod, B = 10000)

plot(res)

该模型似乎低估了低于的概率，而高估了的概率。但同样，由于样本量小，不确定性很大。$0.75$$0.75$

残差

广义线性模型有许多类型的残差，但它们的解释通常很困难。一种可能性是查看 R 包中实现的基于模拟的分位数残差DHARMa。这是使用与上述相同数据的示例：

fit <- glm(y ~ x, data=dat, family=binomial)

simres <- simulateResiduals(fit, n = 1e4, seed = 142857)

plot(simres)

这些残差的好处是它们可以被解释为线性回归模型的“通常”残差。左侧显示了残差的 QQ 图。在右侧，残差是针对预测值绘制的。在这两种情况下，似乎都没有证据表明存在问题。

$[1]:$ Nattino, G.、Finazzi, S. 和 Bertolini, G. (2014)。一种新的校准测试和校准带的重新评估，用于评估基于二分结果的预测模型。医学统计，33（14），2390-2407。

$[2]:$ Nattino, G.、Finazzi, S. 和 Bertolini, G. (2016)。一种新的测试和图形工具，用于评估逻辑回归模型的拟合优度。医学统计，35（5），709-720。

$[3]:$ Nattino, G.、Lemeshow, S.、Phillips, G.、Finazzi, S. 和 Bertolini, G.（2017 年）。使用校准带评估二分结果模型的校准。统计杂志，17（4），1003-1014。

$[4]:$哈雷尔，FE（2015 年）。回归建模策略：应用于线性模型、逻辑和序数回归以及生存分析（第 3 卷）。纽约：斯普林格。