医学相关

ROCR:R包学习(一)

2020-04-25  本文已影响0人  庄小尽

一、预测对象

一组预测和标签

#一个来自prediction帮助文件的简单示例,该示例使用预测和标签矢量(即不是矩阵)。我们看到数据是一些连续的预测和二进制标签:
rm(list = ls())
library(ROCR)
data(ROCR.simple)
head(cbind(ROCR.simple$predictions, ROCR.simple$labels), 5)
#predictions 包含预测的向量、矩阵、列表或数据帧。
#labels 包含真实类别标签的向量、矩阵、列表或数据框。必须与“预测”具有相同的维度
#predictions是“真相”的一些预测度量(通常是连续的),即labels。在许多应用中,predictions估计的概率(或对数赔率)和labels是二进制值。两个参数都可以采用向量,矩阵或data.frame进行预测,但dim(predictions)必须相等dim(labels)

#制作预测对象并显示它的内容:
pred <- prediction(ROCR.simple$predictions,ROCR.simple$labels)
class(pred)
slotNames(pred)#这些函数返回或设置关于对象中各个槽的信息。
sn = slotNames(pred)
sapply(sn, function(x) length(slot(pred, x)))
sapply(sn, function(x) class(slot(pred, x)))
#我们看到每个插槽的长度为1,并且是一个列表。

多组预测和标签

#使用ROCR.hiv数据集来说明如果提供了一组以上的预测和标签,这将如何工作。在这里,我们向prediction函数传递了(10)个预测列表和标签列表:
data(ROCR.hiv)
manypred = prediction(ROCR.hiv$hiv.nn$predictions, ROCR.hiv$hiv.nn$labels)
class(manypred)
sapply(sn, function(x) length(slot(manypred, x)))
sapply(sn, function(x) class(slot(manypred, x)))
#我们看到所有插槽仍然是列表,但是现在它们的长度为(10),对应于(10)个预测/标签。如果2个参数是矩阵,我们将得到相同的结果,但是这将要求所有预测和标签具有相同的长度。使用预测/标签列表更加灵活。

二、Performance objects

performance(prediction.obj, measure, x.measure="cutoff", ...)我们看到第一个参数是一个prediction对象,第二个参数是一个measure。如果运行?performance,则可以看到已实施的所有性能指标。

ROC曲线


#一组预测和标签
#制作一条ROC曲线,该曲线在x轴上绘制误报率(FPR),在y轴上绘制真实误报率(TPR)
roc.perf = performance(pred, measure = "tpr", x.measure = "fpr")
plot(roc.perf)
abline(a=0, b= 1)

#在每个截止点,TPR和FPR都被计算并绘制出来。图越平滑,预测的截断点就越多。我们还绘制了一条45度的直线,它平均表示一个均匀(0,1)随机变量的性能。离对角线越远越好。总的来说,我们看到我们在敏感性(真实阳性率,(>80%))上得到了提高,抵消了假阳性率(1-特异性),直到约15%的FPR。在经历了15%的FPR之后,我们并没有看到TPR在FPR增加的情况下有明显的改善。

#多组预测和标签
many.roc.perf = performance(manypred, measure = "tpr", x.measure = "fpr")
plot(many.roc.perf, col=1:10)
abline(a=0, b= 1)
#具有多个预测和标签plot的performance对象上的函数将遍历列表并为每个列表绘制ROC

限制在一个FPR:部分ROC曲线

#您可能只想接受一定水平的误报率,例如10%。pROC下面的函数将仅使值小于或等于您设置的FPR。
pROC = function(pred, fpr.stop){
    perf <- performance(pred,"tpr","fpr")
    for (iperf in seq_along(perf@x.values)){
        ind = which(perf@x.values[[iperf]] <= fpr.stop)
        perf@y.values[[iperf]] = perf@y.values[[iperf]][ind]
        perf@x.values[[iperf]] = perf@x.values[[iperf]][ind]
    }
    return(perf)
}

# 绘制部分ROC曲线:
proc.perf = pROC(pred, fpr.stop=0.1)
plot(proc.perf)
abline(a=0, b= 1)
#如果我们只能接受10%的FPR,则该模型在10%FPR(1-特异性)下仅给出50%的敏感性(TPR)

获得“最佳”切割点

# 在某些ROC曲线应用中,您希望该点最接近(1)的TPR和(0)的FPR。该切点在“灵敏度”和“特异性”均等权衡的意义上是“最佳”的。要确定此临界值,可以使用以下代码。该代码同时吸收了performance对象和prediction对象,并给出了预测的最佳截止值:
opt.cut = function(perf, pred){
    cut.ind = mapply(FUN=function(x, y, p){
        d = (x - 0)^2 + (y-1)^2
        ind = which(d == min(d))
        c(sensitivity = y[[ind]], specificity = 1-x[[ind]], 
            cutoff = p[[ind]])
    }, perf@x.values, perf@y.values, pred@cutoffs)
}
print(opt.cut(roc.perf, pred))

#现在,cost可以在ROCR软件包中使用一种措施来创建performance对象。如果您使用它来寻找最低成本,那么它将给您与相同的临界值opt.cut,但不会给您带来敏感性和特异性。
cost.perf = performance(pred, "cost")
pred@cutoffs[[1]][which.min(cost.perf@y.values[[1]])]


三、准确性

#一组预测和标签
#另一个流行的度量是总体准确性。此度量可优化正确结果,但如果负数比正数多,则可能会偏斜,反之亦然。让我们获得简单预测的整体准确性并将其绘制出来
acc.perf = performance(pred, measure = "acc")
plot(acc.perf)
#如果我们实际上想提取最大精度和与之相对应的临界值怎么办?在performance对象中,我们有slot x.values,cutoff在这种情况下与y.values相对应,在处,与每个切除的精度相对应。我们将获取索引以获得最大的准确性,然后获取相应的截止值:
ind = which.max( slot(acc.perf, "y.values")[[1]] )
acc = slot(acc.perf, "y.values")[[1]][ind]
cutoff = slot(acc.perf, "x.values")[[1]][ind]
print(c(accuracy= acc, cutoff = cutoff))#您可以使用cutoff(希望)测试数据中的最大准确性来对模型进行阈值处理。

#多组预测和标签

#我们将对许多预测和标签执行相同的操作,但是必须遍历结果(使用mapply语句)
many.acc.perf = performance(manypred, measure = "acc")
sapply(manypred@labels, function(x) mean(x == 1))
mapply(function(x, y){
  ind = which.max( y )
  acc = y[ind]
  cutoff = x[ind]
  return(c(accuracy= acc, cutoff = cutoff))
}, slot(many.acc.perf, "x.values"), slot(many.acc.perf, "y.values"))

曲线下面积(AUC)和部分AUC(pAUC)

#一组预测和标签
#曲线下的面积仅通过获取曲线和x轴之间的面积即可概括ROC曲线。让我们为简单的预测获取曲线下的面积
auc.perf = performance(pred, measure = "auc")
auc.perf@y.values
#您所见的结果是标量数,即曲线下的面积(AUC)。该数字的范围是(0)至(1),其中(1)表示100%的特异性和100%的灵敏度。

#如果您只想接受固定的FPR,则可以使用以下fpr.stop参数来计算部分AUC :
pauc.perf = performance(pred, measure = "auc", fpr.stop=0.1)
pauc.perf@y.values

#现在,我们看到pAUC要低得多。值得注意的是,这个值的范围从(0)到fpr.stop。为了使它标准化为(1),可以将它除以fpr.stop来测量([0,1])
pauc.perf@y.values = lapply(pauc.perf@y.values, function(x) x / 0.1)
pauc.perf@y.values

#尽管这一指标与完整的AUC指标相比更具可比性,但仍处于较低水平。注意,AUC或pAUC没有“一个”截止时间,因为它衡量的是所有截止时间的性能。同样,为标量结果度量(如AUC)绘制函数也不适用于性能对象。多个预测的代码是相同的。
manypauc.perf = performance(manypred, measure = "auc", fpr.stop=0.1)
manypauc.perf@y.values = lapply(manypauc.perf@y.values, function(x) x / 0.1)
manypauc.perf@y.values

上一篇下一篇

猜你喜欢

热点阅读