如何手动计算曲线下面积 (AUC) 或 c 统计量

机器算法验证 回归 物流 分类 奥克
2022-01-28 21:54:38

我对手动计算二元逻辑回归模型的曲线下面积 (AUC) 或 c 统计量感兴趣。

例如,在验证数据集中,我有因变量的真实值,保留(1 = 保留;0 = 不保留),以及我的回归分析生成的每个观察值的预测保留状态,该模型是使用训练集构建(范围从 0 到 1)。

我最初的想法是确定模型分类的“正确”数量,然后简单地将“正确”观察的数量除以总观察的数量来计算 c 统计量。通过“正确”,如果观察的真实保留状态 = 1 并且预测的保留状态 > 0.5,那么这是一个“正确”分类。此外,如果观察的真实保留状态 = 0 并且预测的保留状态 < 0.5,那么这也是“正确”分类。我假设当预测值 = 0.5 时会出现“平局”,但这种现象不会出现在我的验证数据集中。另一方面,如果观察的真实保留状态 = 1 并且预测的保留状态 < 0,则“不正确”分类将是。5 或者如果结果的真实保留状态 = 0 并且预测的保留状态 > 0.5。我知道 TP、FP、FN、TN,但不知道如何根据这些信息计算 c 统计量。

4个回答

我会推荐 Hanley 和 McNeil 1982 年的论文“接收器操作特征 (ROC) 曲线下面积的意义和用途”。

例子

它们具有下表的疾病状态和测试结果(例如,对应于逻辑模型的估计风险)。右边的第一个数字是真实疾病状态“正常”的患者数量,第二个数字是真实疾病状态“异常”的患者数量:

(1) 绝对正常:33/3
(2) 可能正常:6/2
(3) 可疑:6/2
(4) 可能异常:11/11
(5) 绝对异常:2/33

因此,总共有 58 名“正常”患者和“51”名异常患者。我们看到,当预测变量为 1 时,“绝对正常”,患者通常正常(36 名患者中有 33 人为真),而当预测变量为 5,“绝对异常”时,患者通常异常(33 名患者为真) 35 名患者),因此预测器是有意义的。但是我们应该如何判断一个得分为 2、3 或 4 的患者呢?我们设置了判断患者异常或正常的临界值,以确定结果测试的敏感性和特异性。

敏感性和特异性

我们可以计算不同临界值的估计敏感性和特异性。(从现在开始,我将只写“敏感性”和“特异性”,让值的估计性质是隐含的。)

如果我们选择我们的截止值,以便我们将所有患者分类为异常,无论他们的测试结果如何(即我们选择截止值 1+),我们将获得 51/51 = 1 的敏感性。特异性将为 0 /58 = 0。听起来不太好。

好的,所以让我们选择一个不太严格的截止值。我们仅将测试结果为 2 或更高的患者归类为异常。然后我们漏掉了 3 名异常患者,灵敏度为 48/51 = 0.94。但我们的特异性大大提高,为 33/58 = 0.57。

我们现在可以继续这个,选择各种截止值(3、4、5、>5)。(在最后一种情况下,我们不会将任何患者归类为异常,即使他们的最高测试分数为 5。)

ROC曲线

如果我们对所有可能的截止值都这样做,并将灵敏度与 1 减去特异性作图,我们得到 ROC 曲线。我们可以使用以下 R 代码:

# Data
norm     = rep(1:5, times=c(33,6,6,11,2))
abnorm   = rep(1:5, times=c(3,2,2,11,33))
testres  = c(abnorm,norm)
truestat = c(rep(1,length(abnorm)), rep(0,length(norm)))

# Summary table (Table I in the paper)
( tab=as.matrix(table(truestat, testres)) )

输出是:

        testres
truestat  1  2  3  4  5
       0 33  6  6 11  2
       1  3  2  2 11 33

我们可以计算各种统计数据:

( tot=colSums(tab) )                            # Number of patients w/ each test result
( truepos=unname(rev(cumsum(rev(tab[2,])))) )   # Number of true positives
( falsepos=unname(rev(cumsum(rev(tab[1,])))) )  # Number of false positives
( totpos=sum(tab[2,]) )                         # The total number of positives (one number)
( totneg=sum(tab[1,]) )                         # The total number of negatives (one number)
(sens=truepos/totpos)                           # Sensitivity (fraction true positives)
(omspec=falsepos/totneg)                        # 1 − specificity (false positives)
sens=c(sens,0); omspec=c(omspec,0)              # Numbers when we classify all as normal

使用它,我们可以绘制(估计的)ROC 曲线:

plot(omspec, sens, type="b", xlim=c(0,1), ylim=c(0,1), lwd=2,
     xlab="1 − specificity", ylab="Sensitivity") # perhaps with xaxs="i"
grid()
abline(0,1, col="red", lty=2)

AUC曲线

手动计算 AUC

我们可以很容易地计算 ROC 曲线下的面积,使用梯形面积的公式:

height = (sens[-1]+sens[-length(sens)])/2
width = -diff(omspec) # = diff(rev(omspec))
sum(height*width)

结果是 0.8931711。

一致性测量

AUC 也可以看作是一种一致性度量。如果我们取所有可能的一患者,其中一个是正常的,另一个是异常的,我们可以计算出具有最高(最“异常”)测试结果的异常患者的频率(如果它们具有相同的值,我们把这算作“半场胜利”):

o = outer(abnorm, norm, "-")
mean((o>0) + .5*(o==0))

答案还是 0.8931711,即 ROC 曲线下的面积。情况将永远如此。

一致性的图形视图

正如 Harrell 在他的回答中指出的那样,这也有一个图形解释。让我们在y轴上绘制测试分数(风险估计),在x轴上绘制真实疾病状态(这里有一些抖动,以显示重叠点):

plot(jitter(truestat,.2), jitter(testres,.8), las=1,
     xlab="True disease status", ylab="Test score")

风险评分与真实疾病状态的散点图。

现在让我们在左侧的每个点(“正常”患者)和右侧的每个点(“异常”患者)之间画一条线。具有正斜率的线的比例(即一致对的比例)是一致指数(平线算作“50% 一致”)。

由于关系的数量(风险评分相等),这个例子的实际线条有点难以可视化,但是通过一些抖动和透明度,我们可以得到一个合理的图:

d = cbind(x_norm=0, x_abnorm=1, expand.grid(y_norm=norm, y_abnorm=abnorm))
library(ggplot2)
ggplot(d, aes(x=x_norm, xend=x_abnorm, y=y_norm, yend=y_abnorm)) +
  geom_segment(colour="#ff000006",
               position=position_jitter(width=0, height=.1)) +
  xlab("True disease status") + ylab("Test\nscore") +
  theme_light()  + theme(axis.title.y=element_text(angle=0))

风险评分与真实疾病状态的散点图,所有可能的观察对之间都有线。

我们看到大部分线条向上倾斜,因此一致性指数会很高。我们还看到了每种类型的观察对对指数的贡献。大部分来自风险评分为 1 的正常患者和风险评分为 5 的异常患者(1-5 对),但也有不少来自 1-4 对和 4-5 对。根据斜率定义计算实际的一致性指数非常容易:

d = transform(d, slope=(y_norm-y_abnorm)/(x_norm-x_abnorm))
mean((d$slope > 0) + .5*(d$slope==0))

答案又是 0.8931711,即 AUC。

Wilcoxon-Mann-Whitney 检验

一致性测量与 Wilcoxon-Mann-Whitney 检验之间存在密切联系。实际上,后者测试一致性的概率(即,随机的正常-异常对中的异常患者将获得最“异常”的测试结果)是否恰好为 0.5。它的检验统计量只是估计的一致性概率的简单变换:

> ( wi = wilcox.test(abnorm,norm) )
    Wilcoxon rank sum test with continuity correction

data:  abnorm and norm
W = 2642, p-value = 1.944e-13
alternative hypothesis: true location shift is not equal to 0

检验统计量 ( W = 2642) 计算一致对的数量。如果我们将它除以可能对的数量,我们会得到一个熟悉的数字:

w = wi$statistic
w/(length(abnorm)*length(norm))

是的,它是 0.8931711,即 ROC 曲线下的面积。

计算 AUC 的更简单方法(在 R 中)

但是,让我们让自己的生活更轻松。有多种软件包可以自动为我们计算 AUC。

Epi 包

Epi软件包创建了一个漂亮的 ROC 曲线,其中嵌入了各种统计数据(包括 AUC):

library(Epi)
ROC(testres, truestat) # also try adding plot="sp"

Epi 包中的 ROC 曲线

pROC 软件包

我也喜欢这个pROC包,因为它可以平滑 ROC 估计(并根据平滑的 ROC 计算 AUC 估计):

来自 pROC 包的 ROC 曲线(未平滑和平滑)

(红线是原始 ROC,黑线是平滑后的 ROC。还要注意默认的 1:1 纵横比。使用它是有意义的,因为灵敏度和特异性都有 0-1 的范围。)

平滑ROC的估计 AUC为 0.9107,与未平滑 ROC 的 AUC 相似但略大(如果您查看该图,您可以很容易地看出它为什么更大)。(尽管我们确实没有太多可能的不同测试结果值来计算平滑的 AUC)。

有效值包

Harrell 的包可以使用该函数rms计算各种相关的一致性统计。rcorr.cens()它的C Index输出是 AUC:

> library(rms)
> rcorr.cens(testres,truestat)[1]
  C Index 
0.8931711

caTools 包

最后,我们有了caTools包和它的colAUC()功能。与其他软件包相比,它有一些优势(主要是速度和处理多维数据的能力——请参阅?colAUC),这些优势有时会有所帮助。但当然,它给出了与我们一遍又一遍计算的相同答案:

library(caTools)
colAUC(testres, truestat, plotROC=TRUE)
             [,1]
0 vs. 1 0.8931711

caTools 包中的 ROC 曲线

最后的话

许多人似乎认为 AUC 告诉我们测试有多“好”。有些人认为 AUC 是测试正确分类患者的概率。不是从上面的示例和计算中可以看出,AUC 告诉我们有关一系列测试的信息,每个可能的截止值都有一个测试。

并且 AUC 是根据在实践中永远不会使用的截止值计算的。为什么我们应该关心“无意义”截止值的敏感性和特异性?尽管如此,这就是 AUC 的(部分)基础。(当然,如果 AUC非常接近 1,几乎所有可能的测试都会有很大的辨别力,我们都会很高兴。)

AUC 的“随机正常 - 异常”对解释很好(并且可以扩展到生存模型,我们可以看到它是否是具有最高(相对)危险的人最早死亡)。但人们永远不会在实践中使用它。这是一种罕见的情况,人们知道一个人有一个健康的人和一个生病的人,不知道哪个人是生病的人,并且必须决定治疗哪个人。(无论如何,决定很容易;对待估计风险最高的人。)

因此,我认为研究实际的ROC 曲线将比仅查看 AUC 汇总度量更有用。如果您将 ROC 与(估计的)误报和误报成本以及您正在研究的内容的基本比率一起使用,您就可以有所收获。

另请注意,AUC 仅测量歧视,而不是校准。也就是说,它根据风险评分衡量您是否可以区分两个人(一个生病和一个健康)。为此,它只查看相对风险值(或排名,如果您愿意,参见 Wilcoxon-Mann-Whitney 检验解释),而不是您应该感兴趣的绝对值。例如,如果您划分每个风险从您的逻辑模型估计 2,您将获得完全相同的 AUC(和 ROC)。

在评估风险模型时,校准也非常重要。要检查这一点,您将查看风险评分约为 0.7 的所有患者,并查看其中大约 70% 是否确实患病。对每个可能的风险评分执行此操作(可能使用某种平滑/局部回归)。绘制结果,您将获得校准的图形测量。

如果有一个模型具有良好的校准和良好的辨别能力,那么你就开始有一个好的模型。:)

看看这个问题:了解ROC曲线

以下是如何构建 ROC 曲线(来自那个问题):

绘制 ROC 曲线

给定由您的排名分类器处理的数据集

  • 对分数递减的测试示例进行排名
  • (0,0)
  • 对于每个示例(按降序排列) x
    • 如果为正,则向上移动1x1/pos
    • 如果是负数,向右x1/neg

其中分别是正例和负例的分数。posneg

您可以使用此想法使用以下算法手动计算 AUC ROC:

auc = 0.0
height = 0.0

for each training example x_i, y_i
  if y_i = 1.0:
    height = height + tpr
  else 
    auc = auc + height * fpr

return auc

这张漂亮的 gif 动画图片应该更清楚地说明这个过程

构建曲线

卡尔的帖子有很多很好的信息。但在过去的 20 年里,我还没有看到一个 ROC 曲线的例子,它改变了任何人的想法,朝着好的方向发展。在我看来,ROC 曲线的唯一价值是它的面积恰好等于一个非常有用的一致性概率。ROC 曲线本身会诱使读者使用截止值,这是不好的统计实践。

就手动计算指数而言,在轴上绘制并在的连续预测变量或预测概率的图如果将的每个点与的每个点连接起来,具有正斜率的线的比例就是一致性概率。cY=0,1xY=1yY=0Y=1

的任何度量都是不正确的准确度评分规则,应避免使用。这包括正确分类的比例、敏感性和特异性。n

对于 RHmiscrcorr.cens函数,打印整个结果以查看更多信息,尤其是标准错误。

这是计算 AUC 的另一种自然方法,只需使用梯形规则来获得 ROC 曲线下的面积。

AUC 等于随机采样的阳性观察值的预测概率(为阳性)大于随机采样的阴性观察值的概率。您可以使用它在任何编程语言中通过检查正面和负面观察的所有成对组合来非常容易地计算 AUC。如果样本量太大,您也可以随机抽样观察。如果您想使用笔和纸计算 AUC,这可能不是最好的方法,除非您的样本非常少/时间很多。例如在 R 中:

n <- 100L

x1 <- rnorm(n, 2.0, 0.5)
x2 <- rnorm(n, -1.0, 2)
y <- rbinom(n, 1L, plogis(-0.4 + 0.5 * x1 + 0.1 * x2))

mod <- glm(y ~ x1 + x2, "binomial")

probs <- predict(mod, type = "response")

combinations <- expand.grid(positiveProbs = probs[y == 1L], 
        negativeProbs = probs[y == 0L])

mean(combinations$positiveProbs > combinations$negativeProbs)
[1] 0.628723

我们可以使用pROC包进行验证:

library(pROC)
auc(y, probs)
Area under the curve: 0.6287

使用随机抽样:

mean(sample(probs[y == 1L], 100000L, TRUE) > sample(probs[y == 0L], 100000L, TRUE))
[1] 0.62896