使用 pROC包进行ROC分析

PROC

可以进行ROC分析和ROC 曲线展示的R包。


#1. 安装

> install.packages("pROC")

#2. 数据导入

> library(pROC)
> data(aSAH)
> head(aSAH)
   gos6 outcome gender age wfns s100b  ndka
29    5    Good Female  42    1  0.13  3.01
30    5    Good Female  37    1  0.14  8.54
31    5    Good Female  42    1  0.10  8.09
32    5    Good Female  27    1  0.04 10.42
33    1    Poor Female  42    3  0.13 17.40
34    1    Poor   Male  48    2  0.10 12.75

#3. ROC分析

##3.1 使用roc()进行ROC分析

> roc(aSAH$outcome, aSAH$s100b,plot = T)
> roc(outcome ~ s100b, aSAH,plot=T,levels=c("Good", "Poor"))
> roc(controls=aSAH$s100b[aSAH$outcome=="Good"], cases=aSAH$s100b[aSAH$outcome=="Poor"])
Call:
roc.formula(formula = outcome ~ s100b, data = aSAH, plot = T)

Data: s100b in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.7314
Rplot roc()

##3.2 roc()参数详细解释

roc(...)
# S3 method for formula
roc(formula, data, ...)
# S3 method for default
roc(response, predictor, controls, cases,
density.controls, density.cases,
levels=base::levels(as.factor(response)), percent=FALSE, na.rm=TRUE,
direction=c("auto", "<", "="">"), algorithm = 5, quiet = TRUE, 
smooth=FALSE, auc=TRUE, ci=FALSE, plot=FALSE, smooth.method="binormal",
ci.method=NULL, density=NULL, ...)
roc1 <- roc(aSAH$outcome,
            aSAH$s100b, percent=TRUE,
            # 设置auc参数
            partial.auc=c(100, 90), partial.auc.correct=TRUE,
            partial.auc.focus="sens",
            # 设置ci参数
            ci=TRUE, boot.n=100, ci.alpha=0.9, stratified=FALSE,
            # 设置画图参数
            plot=TRUE, auc.polygon=TRUE, max.auc.polygon=TRUE, grid=TRUE,
            show.thres=TRUE
            )

#在原有的图上加ROC曲线
roc2 <- roc(aSAH$outcome, aSAH$wfns,
            plot=TRUE, add=TRUE, percent=roc1$percent)   
  • response:样本类型,一般只有两类( 0 (controls) 和 1 (cases)),可做ROC曲线。当类型过多时,需要使用levels参数指定那些值作为control ,那些作为 case。
  • predictor:样本的预测值,可以是概率、排名之类。
  • controls, cases: 直接提供controls和cases,可以是数值向量,也可以是排好序的向量。
  • formula, data: 通过表达式传入数据框中的值 。
  • levels: controls 和 cases 对应的值,默认为levels(as.factor(response)),剩下的忽略;当输入数据为0、1,默认0为controls。
  • percent:sensitivities, specificities 和 AUC返回形式为百分数。
  • direction:根据两组数据中位数大小确定;“>”: control组中位数值大于cases组;“<”:control组中位数值小于或等于cases组。
  • algorithm:1,也是默认,数量较少;2,数量大于1000时,速度更快;3,C++实现算法,快3-5x; 4 (debug only, slow): 三个算法都运行一遍,确认返回值知否一样。5,迅速选择2或3。
  • smooth:ROC 曲线修饰为平滑曲线。
  • auc:计算AUC,默认TRUE。
  • ci:是否计算置信区间。
  • plot:是否作图。
  • smooth.method, ci.method:smooth 算法。

##3.3 返回ROC计算对象

coords(roc1, "best", ret=c("threshold", "specificity", "1-npv"))
coords(roc2, "local maximas", ret=c("threshold", "sens", "spec", "ppv", "npv"))

##3.4 置信区间

# Of the AUC
ci(roc2)

# Of the curve
sens.ci <- ci.se(roc1, specificities=seq(0, 100, 5))
plot(sens.ci, type="shape", col="lightblue")
plot(sens.ci, type="bars")

# need to re-add roc2 over the shape
plot(roc2, add=TRUE)

# CI of thresholds
plot(ci.thresholds(roc2))

##3.5 roc.test()对ROC进行统计检验

# Test on the whole AUC
> roc.test(roc1, roc2, reuse.auc=FALSE)
DeLong's test for two correlated ROC curves

data:  roc1 and roc2
Z = -2.209, p-value = 0.02718
alternative hypothesis: true difference in AUC is not equal to 0
sample estimates:
AUC of roc1 AUC of roc2 
   73.13686    82.36789 

# Test on a portion of the whole AUC
> roc.test(roc1, roc2, reuse.auc=FALSE, partial.auc=c(100, 90),
         partial.auc.focus="se", partial.auc.correct=TRUE)

    # With modified bootstrap parameters
> roc.test(roc1, roc2, reuse.auc=FALSE, partial.auc=c(100, 90),
         partial.auc.correct=TRUE, boot.n=1000, boot.stratified=FALSE)

##3.6 roc.test()参数

alternative:“two.sided”, “less” ,“greater”。对于method="venkatraman",只能使用 “two.sided”
paired: 是否时配对,会自动检测。
boot.n:指定method="bootstrap"中自检举重复次数和 method="venkatraman"中置换次数;默认,2000。
boot.stratified:每次自检举过程中,cases/controls 比例与原始样本中比例一致。
method

  • “delong”,默认,不适用于partial AUC, smoothed curves,curves with different direction。
  • “bootstrap” :
  • “venkatraman”:不适用于smoothed ROC curves, or curves with partial AUC specifications

##3.7 roc.test()中统计学方法

  • method="bootstrap"
  1. 自检举抽样,通过boot.n指定自检举次数。

  2. 计算AUC。

  3. 计算标准偏差。

    D=\frac{AUC1-AUC2}{s}

  4. D与正态分布进行比较

  • method="delong"
    method="delong" 在文章DeLong *et al.* (1988) 中提到用于配对ROC 曲线的检验,实际利用是在文章 Sun and Xu (2014);现在已经延伸用于非配对的 ROC 曲线研究中,

    D=\frac{V^r(\theta^r) - V^s(\theta^s) }{ \sqrt{S^r + S^s}}

  • method="venkatrama"
    method="venkatraman"在文章Venkatraman and Begg (1996) (for paired ROC curves)Venkatraman (2000) (for unpaired ROC curves)对样本排名进行置换检验。通过boot.n指定置换次数。

##3.7 样本量

# Two ROC curves
power.roc.test(roc1, roc2, reuse.auc=FALSE)
power.roc.test(roc1, roc2, power=0.9, reuse.auc=FALSE)

# One ROC curve
power.roc.test(auc=0.8, ncases=41, ncontrols=72)
power.roc.test(auc=0.8, power=0.9)
power.roc.test(auc=0.8, ncases=41, ncontrols=72, sig.level=0.01)
power.roc.test(ncases=41, ncontrols=72, power=0.9)

#4 pROC使用更多实例

EXPASY基于R代码上给出了pROC的6个示例,见pROC: Screenshots,下面看一个例子:

library(pROC)
data(aSAH)
rocobj <- plot.roc(aSAH$outcome, aSAH$s100b, percent = TRUE, main="Smoothing")
lines(smooth(rocobj), # smoothing (default: binormal)
col = "#1c61b6")
lines(smooth(rocobj, method = "density"), # density smoothing
col = "#008600")
lines(smooth(rocobj, method = "fitdistr", # fit a distribution
density = "lognormal"), # let the distribution be log-normal
col = "#840000")
legend("bottomright", legend = c("Empirical", "Binormal", "Density", "Fitdistr\n(Log-normal)"), col = c("black", "#1c61b6", "#008600", "#840000"),lwd = 2)
Smoothing

#5 参考:

pROC: Display and Analyze ROC Curves
pROC: an open-source package for R and S+ to analyze and compare ROC curves
例子参考

系列文章:
使用yardstick包进行ROC分析
使用 pROC包进行ROC分析

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。