7.评价预测模型——ROC曲线_roc曲线下面积多少才有预测价值-CSDN博客

library(readxl)
data <- read_excel("data.xlsx")
data<-na.omit(data)
data<-as.data.frame(data)

form.bestglm<-as.formula(group~age+BMI+ToS+CA153+CDU+transfusion+stage)
form.all<-as.formula(group~.)

使用最优子集发筛选的自变量，构成模型form.bestglm，然后构建一个所有自变量的模型公式form.all

对数据集data进行打包：

library(rms)
dd=datadist(data)
options(datadist="dd")

模型进行拟合：

fit.glm<- lrm(formula=form.bestglm,data=data,x=TRUE,y=TRUE)

（1）roc函数法

ROC拟合：

predvalue<-predict(fit.glm)
#install.packages("pROC")
library(pROC)
ROC <- roc(data$group,predvalue)

predict（）函数进行预测，将预测结果储存在predvalue中

roc（）函数，第一个为结局变量，第二个为预测结果。

查看ROC结果：

ci(auc(ROC))
coords(ROC, x="best", ret="all", transpose = FALSE)

auc（）查看曲线下面积；ci（）查看曲线下面积的置信区间95%CI；

此处的AUC及其置信区间即C指数及其可信区间。

查看更多的ROC结果：

coords（）函数中输入模型ROC，选项x设置为“best”，选项ret设置为“all”，选项transpose设置为F。

输出的结果可以看到截断值（threshold）为-1.30197，特异度（specificity）为0.808，灵敏度（sensitivity）为0.655，准确度（accuracy）为0.783，真阳性（True Positive, TP）57例，真阴性（True Negative, TN）346例，假阴性（False Negative, FN）30例，假阳性（False Positive, FP） 82例，约登指数（r, Youden’s index）1.46，阳性预测值（positive predictive value, ppv）0.41，阴性预测值（npv) 0.92,

诊断的结果可整理成四格表的形式。表中有4个可能结果，其中2个结果表明被评价诊断方法的诊断结果是正确的，即病例被诊断为阳性（真阳性；True Positive, TP）和对照被诊断为阴性（真阴性；True Negative, TN）；两个是错误的，即病例被诊断为阴性（假阴性；False Negative, FN）和对照被诊断为阳性（假阳性；False Positive, FP）。

1. 患病率

患病率即研究对象中患有该疾病的人所占的比例。不同于发病率（incidence），患病率的计算既包括新病例也包括旧病例。

2. 准确度

准确度（符合率），即与金标准相比，待评价的诊断方法获得正确结果的比例，即真阳性和真阴性的比例。

准确度受患病率的影响，且没有揭示假阳性和假阴性各自的概率。因此，准确度是比较粗略的评价指标。

3. 灵敏度和特异度

灵敏度，即实际患病而按该诊断方法被正确地判定为患病的百分率（真阳性率）。灵敏度反映该诊断方法发现病人的能力。

        特异度，即实际无病按该诊断方法被正确地判定为无病的百分率（真阴性率）。特异度反映该诊断方法确定不患病的能力。

灵敏度与特异度具有不受患病率影响的优点，其取值范围均在（0, 1）之间，其值越接近于1，说明该项诊断试验的价值越好，对于有病、无病的识别能力越强。但由于灵敏度和特异度随界值的影响此消彼长，所以无法比较两种方法的诊断价值。

4. 误诊率和漏诊率

误诊是指将不患病错误地判定为患病，即假阳性。

误诊率（%）=假阳性例数/实际阴性例数×100%

漏诊是指将患病错误地判定为不患病，即假阴性。

漏诊率（%）=假阴性例数/实际阳性例数×100%

5. 约登指数（r, Youden’s index）

约登指数，又称为正确诊断指数，表示诊断方法的真实度。其值越接近1，诊断试验的真实性越好。当漏诊和误诊对疾病诊断的危害性具有同等意义的时候，可以用约登指数来比较两种诊断方法的优劣。

约登指数=1-（假阳性率+假阴性率）=（灵敏度+特异度）-1

6. 似然比

阳性似然比（positive likelihood ratios, LR+）反映被评价的诊断方法正确判断阳性的可能性是错误判断阳性可能性的倍数，即真阳性率与假阳性率之比。阳性似然比越大，说明诊断结果为阳性时，患者真的患病的概率越大，诊断价值越大。

        阴性似然比（negative likehood ratios, LR-）反映被评价的诊断方法错误判断阴性的可能性是正确判断阴性的可能性的倍数，即假阴性率与真阴性率之比。阴性似然比越小，说明能够否定患有该病的可能性越大，诊断价值越大。

7. 预测值

阳性预测值（positive predictive value, PV+）：真阳性在待评价方法判定的阳性患者中所占的比例；反映待评价的诊断方法判定为阳性时，该受试对象确实患某病的概率。

         阴性预测值（negative predictive value, PV-）：真阴性在待评价方法判定的阴性非患者中所占的比例；反映待评价的诊断方法判定为阴性时，该受试者不会患某病的概率。

预测值受患病率的影响很大。

8. 优势比(odds ratio, OR)

病例组中诊断阳性数与阴性数比值与对照组中诊断阳性数与阴性数比值的比值。

OR>1，表示试验组比对照组更容易获得诊断阳性结果；反之亦然。

ROC绘图

par(mar=c(5,5,2,1))
plot(1-ROC$specificities,
     ROC$sensitivities,
     type="l",col="black",
     lty=1,lwd=2,
     xlab="1-Specificities",ylab="Sensitivities",
     xlim = c(0,1),
     ylim = c(0,1),
     xaxs = "i",
     yaxs = "i",
     cex=1.5,
     cex.lab=1.5,
     cex.axis=1.5)
abline(0,1,lty=2,lwd=2)

（2）AUC函数法

AUC曲线

#install.packages("modEvA")
library(modEvA)
pred<-exp(predvalue)/(1+exp(predvalue))

将处理的数据进行logit变换，处理成为概率。

AUC(obs=data$group,
    pred=pred,
    plot.values = T,
    curve = "ROC", 
    main = "ROC",
    diag.col="black",
    diag.lty=1,
    curve.col="black",
    curve.lty=1,
    curve.lwd=2,     
    xaxs = "i",
    yaxs = "i")

obs（）指定分组变量；pred指定预测概率；plot.values=F表示不显示AUC数值结果；curve指定绘制ROC曲线，main标题；diag.col参考线颜色，diag.lty参考线类型，curve.col, curve.lty, curve.lwd 设置ROC曲线的颜色，类型，宽度，xaxs，yaxs设置原点相交。

PR曲线

绘制PR曲线：

AUC(obs=data$group,
    pred=pred,
    plot.values = FALSE,
    curve = "PR", 
    main = "PR curve",
    diag.col="black",
    diag.lty=1,
    curve.col="black",
    curve.lty=1,
    curve.lwd=2,     
    xaxs = "i",
    yaxs = "i")

PR曲线是以precision（准确率，作为纵坐标）和recall（灵敏度（召回率），作为横坐标）这两个变量作出的曲线。

$R=\frac{TP}{TP+FN}$ $P=\frac{TP}{TP+FP}$

P和R是一对矛盾的度量，P高时往往R偏低，R高时往往P偏低。

1.若曲线a完全“包住”曲线b，则a性能优于b。

2.若曲线a和曲线b发生了相交，则无法比较。

3.比较曲线a和曲线b曲线下面积的大小。面积大则性能优。

4.比较“平衡点”。用的较少。

“平衡点”（Break-Event Point，简称BEP，F1）时P值和R值相同的时的取值。F1值越大，我们可以认为PR曲线较好。

F1=2*P*R/（P+R）

（3）roc.curve函数法

#install.packages("PRROC")
library(PRROC)
x=pred[data$group==1]
y=pred[data$group==0]

将数据中group==1的样本预测概率传递给x，group==0的样本预测概率传递给y

PR曲线

pr<-pr.curve(scores.class0=x,
             scores.class1=y,
             curve = TRUE )
print(pr)
plot(pr,
     xaxs = "i",
     yaxs = "i")

ROC曲线

roc <-roc.curve(scores.class0=x,
                scores.class1=y,
                curve = TRUE )
print(roc)
plot(roc,
     xlab="1-Specificities",
     ylab="Sensitivities",
     xaxs = "i",
     yaxs = "i")

（4）Score函数

注意：score（）函数只能使用glm（）拟合的模型

fit2.glm = glm(formula=form.bestglm,data=data,family=binomial())

创建Score函数的结果xb

library(riskRegression)
xb <-Score(list("Model Bestglm"=fit2.glm),
           formula=group~1,
           plots=c("calibration","ROC"),
           metrics = c("auc", "brier"),
           B=1000,M=50,
           data=data)

object指定fit2.glm，模型通过list包装，命名为Model Bestglm。

formula设置公式左侧的内容，需要与前文保持一致。

plots设置可绘制哪些图形；

metric设置可得到哪些统计量；

选项B设置bootstrap次数；

选项M设置每次bootstrap的样本大小。

绘制ROC曲线：

plotROC(x=xb,
        xlab="1-Specificities",
        ylab="Sensitivities",
        col=c("#E64B35B2"),
        brier.in.legend=FALSE,
        auc.in.legend=FALSE,
        cex=1.2)

x指定Score（）函数拟合xb；xlab，ylab设置x，y轴的名称。选项brier.in .legend=F，表示不显示brier得分。选项auc.in.legend=F，表示不显示auc。

2. 多模型ROC曲线

将多个ROC绘制在一个图中，创建不同的两个模型

form.bestglm<-as.formula(group~age+BMI+ToS+CA153+CDU+transfusion+stage)
form.all<-as.formula(group~.)
fit.glm<- lrm(formula=form.bestglm,data=data,x=TRUE,y=TRUE)
fit3.glm<-lrm(formula=form.all,data=data,x=TRUE,y=TRUE)
predvalue <- predict(fit.glm)
predvalue3<-predict(fit3.glm)

（1）roc（）函数

创建两个ROC的结果

ROC <- roc(data$group,predvalue)
ROC2 <- roc(data$group,predvalue2)

绘制ROC曲线

par(mar=c(5,5,2,1))
plot(1-ROC$specificities,ROC$sensitivities,
     type="l",col="#E64B35B2",lty=1,lwd=2,
     xlab="1-Specificities",ylab="Sensitivities",
     xlim = c(0,1),
     ylim = c(0,1),
     xaxs = "i",
     yaxs = "i",
     cex=1.5,
     cex.lab=1.5,
     cex.axis=1.5)
lines(1-ROC2$specificities,ROC2$sensitivities,col="#4DBBD5B2",lty=2,lwd=2)
abline(0,1)
legend(0.58,0.25,
       legend=c("Model Bestglm","Model All"),
       lty = c(1,2),
       lwd = c(2,2),
       col = c("#E64B35B2","#4DBBD5B2"),
       bty="n") #"o"为加边框

比较两个模型：

roc.test(ROC,ROC2)

（2）Score（）函数

glm（）函数拟合模型：

form.bestglm<-as.formula(group~age+BMI+ToS+CA153+CDU+transfusion+stage)
form.all<-as.formula(group~.)
fit.glm<- glm(formula=form.bestglm,data=data,family = binomial())
fit2.glm = glm(formula=form.all,data=data,family=binomial())

绘制ROC曲线

library(riskRegression)
xb <-Score(list("Model Bestglm"=fit.glm,
                "Model All"=fit2.glm),
           formula=group~1,
           plots=c("calibration","ROC"),
           metrics = c("auc", "brier"),
           B=1000,M=50,
           data=data)
plotROC(x=xb,
        xlab="1-Specificities",
        ylab="Sensitivities",
        col=c("#E64B35B2","#4DBBD5B2"),
        brier.in.legend=FALSE,
        auc.in.legend=FALSE,
        cex=1.2)

xb2$AUC

生存资料

多模型多时间点ROC曲线通常不绘制在一幅图中，整体下来较推荐使用timeROC函数，可以进行time-AUC模型，可以compare比较两组ROC曲线

1. timeROC包中的timeROC（）函数

2. survivalROC包中的survivalROC（）函数

案例：原发性胆汁性肝硬化的生存分析研究

1.单模型单时间点ROC曲线

（1）TimeROC函数

创建一个生存模型fit.cox

load("pbc.Rdata")
pbc<-na.omit(pbc)
library(rms)
dd <- datadist(pbc)
options(datadist = "dd")
fit.cox <- cph(formula=Surv(days,status) ~ascites+edema+bili+albumin+copper+prothrombin+chol,
               data=pbc,
               x=TRUE,y=TRUE,surv=TRUE)

产生一个预测值

predvalue<-predict(fit.cox)

创建一个ROC

#install.packages("timeROC")
library(timeROC)
ROC.marginal<-timeROC(T=pbc$days, 
                           delta=pbc$status,
                           marker=predvalue, 
                           cause=1, 
                           weighting="marginal", 
                           times=c(365*1,365*3,365*5), 
                           iid=TRUE)

T指定随访时间，delta指定随访结局，marker指定预测的变量predvalue，cause指定1为结局时间，weighting指定计算方法，marginal表示KM模型。time指定计算的ROC曲线的时间节点，iid=T表示保存置信区间。

显示ROC曲线下面积及置信区间

ROC.marginal$AUC
confint(ROC.marginal)$CI_AUC

confint（）函数计算AUC的可信区间

从输出结果可以看到，1年的AUC为0.9016998，95%CI可信区间为80.58%-99.76%，2年的AUC为0.8937，3年AUC为0.8864

绘制ROC曲线

plot(ROC.marginal,
     time=365*3,
     lty=1,lwd=2,
     xlab="1-Specificities",ylab="Sensitivities",
     col = "#E64B35B2",
     title="")

表示为模型3年的ROC。

（2）survivalROC函数

library(survivalROC)
survROC= survivalROC(Stime=pbc$days,  
                     status=pbc$status,      
                     marker = predvalue,     
                     predict.time =365, method="KM")

marker指定预测变量，stime指定生存时间，status指定生存结局，predict.time指定计算时点

查看AUC值，看不了95%CI

survROC$AUC

绘制ROC曲线

plot(survROC$FP, survROC$TP,
     type="l",col="#E64B35B2",lty=1,lwd=2,
     xlab="1-Specificities",ylab="Sensitivities",
     xlim = c(0,1),
     ylim = c(0,1),
     xaxs = "i",
     yaxs = "i",
     cex=1.5,
     cex.lab=1.5,
     cex.axis=1.5)
abline(0,1,lty=2,lwd=2)

FP假阳性率，特异度为真阴性率，1-特异度为假阳性率，灵敏度为真阳性率。

2.单模型多时间点ROC曲线

（1）TimeROC函数

#载入数据
load("pbc.Rdata")
pbc<-na.omit(pbc)#删除缺失值
#打包
library(rms)
dd=datadist(pbc)
options(datadist="dd")
#创建生存模型
fit.cox <- cph(formula=Surv(days,status) ~ascites+edema+bili+albumin+copper+prothrombin+chol,
               data=pbc,
               x=TRUE,y=TRUE,surv=TRUE)
#创建预测值
predvalue<-predict(fit.cox)

创建ROC

library(timeROC)
ROC.marginal<-timeROC(T=pbc$days, 
                           delta=pbc$status,
                           marker=predvalue, 
                           cause=1, 
                           weighting="marginal", 
                           times=c(365*1,365*3,365*5), 
                           iid=TRUE)

多个时间点

plot(ROC.marginal,time=365*1,col = "#E64B35B2",lwd=2,title="")        
plot(ROC.marginal,time=365*3,add=TRUE,col="#4DBBD5B2",lwd=2)
plot(ROC.marginal,time=365*5,add=TRUE,col="#00A087B2",lwd=2)
legend("bottomright",
       legend=c("t=365","t=1095","t=1825"),
       col=c("#E64B35B2","#4DBBD5B2","#00A087B2"),
       lty=1,lwd=2,
       bty = "n")

time设置需要绘制的时间点。add=T表示在原图的基础上增加曲线

绘制时间-AUC曲线

plotAUCcurve(ROC.marginal,conf.int=TRUE,col="#00A087B2")

提示随着时间的改变ROC曲线下面积的变化。

（2）survival函数

创建ROC：

library(survivalROC)
#创建第一个时间点ROC
survROC= survivalROC(Stime=pbc$days,  
                     status=pbc$status,      
                     marker = predvalue,     
                     predict.time =365, method="KM")
#创建第二个时间点ROC 365*3天
survROC2= survivalROC(Stime=pbc$days,  
                      status=pbc$status,      
                      marker = predvalue,     
                      predict.time =365*3, method="KM")
#创建第三个时间点ROC365*5
survROC3= survivalROC(Stime=pbc$days,  
                      status=pbc$status,      
                      marker = predvalue,     
                      predict.time =365*5, method="KM")

绘制曲线

plot(survROC$FP, survROC$TP,
     type="l",col="#E64B35B2",lty=1,lwd=2,
     xlab="1-Specificities",ylab="Sensitivities",
     xlim = c(0,1),
     ylim = c(0,1),
     xaxs = "i",
     yaxs = "i",
     cex=1.5,
     cex.lab=1.5,
     cex.axis=1.5)
lines(survROC2$FP, survROC2$TP,col="#4DBBD5B2",lty=1,lwd=2)
lines(survROC3$FP, survROC3$TP,col="#00A087B2",lty=1,lwd=2)
abline(0,1)
legend(0.58,0.25,
       legend=c("time 365","time 365*3","time 365*5"),
       lty = c(1,1,1),
       lwd = c(2,2,2),
       col = c("#E64B35B2","#4DBBD5B2","#00A087B2"),
       bty="n")

3.多模型同时点ROC曲线

（1）TimeROC函数

创建两个不同的生存模型

#载入数据
load("pbc.Rdata")
pbc<-na.omit(pbc)#删除缺失值
#打包
library(rms)
dd=datadist(pbc)
options(datadist="dd")
#创建生存模型
fit1.cox <- cph(formula=Surv(days,status) ~ascites+edema+bili+albumin+copper+prothrombin+chol,
               data=pbc,
               x=TRUE,y=TRUE,surv=TRUE)
fit2.cox <- cph(formula=Surv(days,status) ~ascites+edema+bili+albumin,
                data=pbc,
                x=TRUE,y=TRUE,surv=TRUE)
#创建预测值
predvalue<-predict(fit1.cox)
predvalue2<-predict(fit2.cox)

根据生存模型，创建ROC

ROC1.marginal <- timeROC(T=pbc$days,
                         delta = pbc$status,
                         marker = predvalue1,
                         cause = 1,
                         weighting = "marginal",
                         times = c(365*1,365*3,365*5),
                         iid = T)
ROC2.marginal<-timeROC(T=pbc$days,   
                              delta=pbc$status,  
                              marker=predvalue2, 
                              cause=1,
                              weighting="marginal", 
                              times=c(365*1,365*3,365*5),
                              iid=TRUE)

两个模型的ROC曲线下面积比较

ROC1.marginal$AUC
ROC2.marginal$AUC
compare(ROC1.marginal, 
        ROC2.marginal, 
        adjusted = TRUE)

Non-adjusted显示的是未校正的p值，Adjusted显示的是经校正的p值。

对第一年的两个模型相比，校正p为0.7137743，不存在统计学差异。

绘制ROC曲线

plot(ROC1.marginal,time=365*3,col = "#E64B35B2",lwd=2,title = "")        
plot(ROC2.marginal,time=365*3,add=TRUE,col="#4DBBD5B2",lwd=2)
legend("bottomright",
       legend=c("Model 1","Model 2"),
       col=c("#E64B35B2","#4DBBD5B2"),lty=1,lwd=2,bty = "n")

（2）survival函数

#载入数据
load("pbc.Rdata")
pbc<-na.omit(pbc)#删除缺失值
#打包
library(rms)
dd=datadist(pbc)
options(datadist="dd")
#创建生存模型
fit1.cox <- cph(formula=Surv(days,status) ~ascites+edema+bili+albumin+copper+prothrombin+chol,
               data=pbc,
               x=TRUE,y=TRUE,surv=TRUE)
fit2.cox <- cph(formula=Surv(days,status) ~ascites+edema+bili+albumin,
                data=pbc,
                x=TRUE,y=TRUE,surv=TRUE)
#创建预测值
predvalue1<-predict(fit1.cox)
predvalue2<-predict(fit2.cox)

创建ROC模型

survROC.1=survivalROC(Stime = pbc$days,
                      status = pbc$status,
                      marker = predvalue1,
                      predict.time = 365*3,
                      method = "KM")
survROC.2= survivalROC(Stime=pbc$days,  
                      status=pbc$status,      
                      marker = predvalue2,     
                      predict.time =365*3, method="KM")

绘图曲线

plot(survROC2$FP, survROC2$TP,
     type="l",col="#E64B35B2",lty=1,lwd=2,
     xlab="1-Specificities",ylab="Sensitivities",
     xlim = c(0,1),
     ylim = c(0,1),
     xaxs = "i",
     yaxs = "i",
     cex=1.5,
     cex.lab=1.5,
     cex.axis=1.5)
lines(survROC.2$FP, survROC.2$TP,col="#4DBBD5B2",lty=1,lwd=2)
abline(0,1)
legend(0.58,0.25,
       legend=c("Model1:time 365*3","Model2:time 365*3"),
       lty = c(1,1),
       lwd = c(2,2),
       col = c("#E64B35B2","#4DBBD5B2"),
       bty="n")

模型1的面积较模型3年的预测准确性面积更大，但统计学没有差异。

竞争风险模型

案例：探究骨髓移植相较血液移植治疗白血病的疗效，结局事件为复发，某些患者因移植后不良反应死亡，这些死亡阻碍了复发的观察，移植相关性死亡与复发存在竞争风险。

单模型单时点ROC曲线

拟合竞争风险模型

library(riskRegression)
library(prodlim)
fgr<-FGR(formula=Hist(ftime,Status)~Sex+D+Phase+Age,data = bmt,cause = 1)
summary(fgr)

利用riskResgression包中的FGR（）函数和prodlim包中的函数Hist创建一个竞争风险模型。

拟合ROC结果

xc <- Score(list('Model1'=fgr),
            formula=Hist(ftime,Status)~1,
            se.fit=TRUE,
            times=c(12,36),
            plots=c("calibration","ROC"),
            metrics ="auc",
            B=1000,M=50,
            data = bmt)

xc

object 指定fgr，模型通过list包装，并命名为model1；formula设置公式，只包含截距；se.fit=T表示估计标准误；times设定计算的时点；plot表示需要绘制的图形；metric指定可得到哪些统计量；B设置bootstrap次数；M设置每次bootstrap的样本大小

12个月的曲线下面积为0.703，36个月的曲线下面积为68.4

绘制ROC曲线

plotROC(xc, times = 12,col="#E64B35B2",lwd=2,
        xlab="1-Specificity",ylab="Sensitivity",
        legend="",
        auc.in.legend = FALSE)

auc.in.legend = T表示显示AUC的值

单模型多时点ROC曲线

plotROC(x=xc, times = 12,col="#E64B35B2",lwd=2,
        xlab="1-Specificity",ylab="Sensitivity",
        legend="",
        auc.in.legend = FALSE)
plotROC(xc,times = 36,col="#4DBBD5B2",lwd=2,
        add=TRUE, 
        legend = '',
        auc.in.legend = F)
legend(0.8,0.2,
       legend=c("1 Year","3-Year"),
       col=c("#E64B35B2","#4DBBD5B2"),
       lwd=2,
       bty='n',
       title='Model1')

绘制Time-AUC曲线

xc2 <- Score(list('Model1'=fgr),
             formula=Hist(ftime,Status)~1,
             se.fit=TRUE,
             times=seq(0,36,1),
             plots=c("calibration","ROC"),
             metrics ="auc",
             B=1000,M=50,
             data = bmt)
plotAUC(xc2,
        conf.int =TRUE,
        col="#E64B35B2")

注意Times的设定

多模型同时点ROC曲线

fgr2<-FGR(Hist(ftime,Status)~Sex+D+Phase,data = bmt,cause = 1)

xc3 <- Score(list('Model1'=fgr,"Model2"=fgr2),
            formula=Hist(ftime,Status)~1,
            se.fit=TRUE,
            times=c(12,36),
            plots=c("calibration","ROC"),
            metrics ="auc",
            B=1000,M=50,
            data = bmt)

xc3

plotROC(x=xc3,times = 12,
        xlab="1-Specificity",
        ylab="Sensitivity",
        col=c("#E64B35B2","#4DBBD5B2"),
        legend=c('Model1','Model2'), 
        auc.in.legend = FALSE)