使用 pROC包进行ROC分析

2019-06-06 本文已影响26人 JeremyL

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"

自检举抽样，通过boot.n指定自检举次数。
计算AUC。
计算标准偏差。
$D=\frac{AUC1-AUC2}{s}$
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