一、ROC分析的基本概念

上面介绍了诸如Youden指数、阳性似然比、阴性似然比、阳性预测值、阴性预测值等综合考虑灵敏度和特异度的诊断性评价指标，但这些指标只能表达指定某一特定诊断界点时所对应的指标，当改变诊断界点时，就会得到不同的指标值，不便于评价整个诊断系统的准确性。因而有人引入ROC分析方法。

ROC是受试者工作特征（Receiver Operating Characteristic）的缩写或相对工作特征（Relative Operating Characteristic）的缩写。它起源于上世纪50年代统计决策理论，最早用于描述雷达信号和噪声之间的关系，并用来比较不同雷达之间的性能差异，后在气象、材料检验、心理物理学及医学检验诊断等领域都有广泛的应用。

ROC分析的定义是：对于可能或将会存在混淆的两种条件或自然状态，需要试验者、专业诊断学工作者以及预测工作者作出精确判别或者准确决策的一种定量方法。它结合不同诊断界点下所得到的灵敏度和特异度对整个诊断系统进行综合评价，根据所绘制出的曲线的形状和面积对诊断试验作定量分析，并进一步评价检测方法的诊断价值大小。

二、ROC曲线的计算与构建

下面结合例子说明如何绘制和计算ROC曲线。

【例1】有人欲探讨果糖胺（FTA）对糖尿病（DM）的诊断价值，收集无糖尿病的健康中老年人血液标本55例和有糖尿病的中老年人血液标本74例，测定FTA的具体含量。研究者选取四个诊断界点（1.30、1.50、1.70、1.90mmol/L），依据这四个界点，将FTA检测值从小到大分成五部分，分别按正常（<1.30）、大致正常（1.30-）、可疑（1.50-）、大致异常（1.70-）、异常（31.90）5个等级分类评估患糖尿病的可能性。试作ROC曲线图，计算曲线下面积（A）与标准误（SE）。

通过诊断试验所获得的资料可分为连续性资料和有序分类资料，本例所测得的FTA检测值即为连续性资料；如在影像诊断中，比如通过看X光片从而判断结果为“正常、大致正常、可疑、大致异常、异常”，则这种资料为有序分类资料。在实际工作中，如果连续性资料样本量较大，且相同值较少时，可将资料整理成频数表的形式，进行简化计算。但由于计算机的普及和统计软件的广泛使用，对于连续性资料没有必要再转化为频数表的形式进行计算，直接输入原始检测值即可，这样就避免了由于转化而损失大量有用的信息，导致检验效能的降低。本例出于讲解及简化计算的需要，故按转化为有序分类资料的形式进行计算。

分别考虑将诊断界值取从高到低的前4个分类为界值计算ROC工作点。该类及以上例数的和为阳性，该类以下例数的和为阴性。诊断界值取1.90、1.70、1.50及1.30时，可以整理出以下4个四格表（以简化的形式表达）。

在不同的诊断界值作为判断标准的情况下，都能计算出相应的真阳性率（即灵敏度）和假阳性率（即1-特异度）。整理后见表5-4。

以假阳性率为横轴，以真阳性率为纵轴，横轴与纵轴的长度相等且均为1，形成正方形，在坐标系上分别标上4个诊断界值所对应的工作点，即（0.0000，0.5946）、（0.0182，0.8108）、（0.1818，0.9596）、（0.5091，0.9865），同时标上（0,0）和（1,1）两点，这两点分别对应于两种极端情况，即灵敏度为0而特异度为1和灵敏度为1而特异度为0，将以上六点各相邻两点用线段连接，即构建出ROC曲线。见图5-2。

图5-2       ROC曲线图式样

ROC曲线对诊断系统的准确性以直观的印象，曲线下面积反映诊断系统的准确性。理论上，此面积的取值范围为0.5至1，完全无价值的诊断系统，其真阳性率与假阳性率相等且始终为0.5，相当于从原点到（1,1）点的对角线，这条线又称为机会线，其下面积为0.5；完善的诊断系统相当于金标准，其真阳性率始终为1，假阳性率始终为0，相当于从原点垂直上升到（0,1）点，然后水平到达（1,1）点，其下面积为1。一般认为曲线下面积在0.50～0.70之间，表示诊断价值较低；在0.70～0.90之间表示诊断价值为中等；0.90以上表示诊断价值较高。

曲线下面积（A_Z）可以采用Hanley和McNeil非参数法进行计算，其计算公式为：

假设观察值越大，异常者越多，异常组有n_a个观察值，记为；正常组有n_n个观测值，记为。异常组的某个与正常组的某个相比，如果前者大于后者，则得分为1，如果相等，则得分为0.5，如果前者小于后者，则得分为0。将n_a×n_n次比较得分加和后取平均，即得A_Z。如果观察值越小，异常者越多，此时，改变计算公式中的大于和小于号即可。

A_Z的标准误SE(A_Z)计算公式为：

（3）

Q₁是两个随机选择的病例组观察值比一个随机选择的对照组观察值都将有更大可能判为异常的概率。Q₂是一个随机选择的病例组观察值比两个随机选择的对照组观察值将有更大可能划归为病例的概率。Q1和Q2的计算详见下表及公式。

已知ROC曲线图中从原点到右上角机会线下的面积为0.5，已计算出的ROC曲线下面积是否与0.5有统计学差异还需进一步作假设检验，以评价其诊断价值。

表5-5   进行ROC曲线分析的实验数据及中间计算结果

将A_z、Q₁、Q₂、n_n、n_a代入SE(A_Z)计算公式：

对ROC曲线下面积A_z与0.5的差别进一步作假设检验，H₀：A_z=0.5；H₁：A_z≠0.5，a=0.5。

，P<0.0001。A_z的95%置信区间为0.9636±1.96×0.0147，即（0.9348，0.9924），大于0.90，表明果糖胺对糖尿病的诊断价值较高。

利用SAS软件logistic过程中的模型选项“outroc=”可以建立包含各临界点所对应的灵敏度和1-特异度等信息的数据集，并进而调用gplot过程绘制ROC曲线。所需的SAS程序如下。

data rocplot;

do y=1,0;

do a=1 to 5;

input f@@;

do i=1 to f;

output;

end;end;end;cards;

1 2 11 16 44

27 18 9 1 0

;

run;

proc logistic descending;

model y=a/outroc=roc1;

run;

goptions dev=win;

symbol1 i=join v=none;

proc gplot data=roc1;

title ‘roc curve’;

axis1 length=4 in order=0 to 1 by 0.2;

plot _sensit_ * _1mspec_=1/vaxis=axis1 haxis=axis1;

run;

【程序修改指导】本程序借用logistic过程产生绘制ROC曲线所需要的数据，在proc logistic语句后面加上descending语句，是为了使SAS按照y=1的概率拟合模型，否则SAS要按照y=0的概率拟合模型。model语句中涉及到两个变量，即结果变量y（1代表金标准异常，0代表金标准正常）和原因变量a,若为有序分类资料，可将资料整理成程序中数据的形式，注意第一行代表金标准为异常的而判断为各个等级的频数，第二行代表金标准为正常的而判断为各个等级的频数；若为连续性资料，则不需要在数据步中用do循环语句，只需用“input y a  @@；”语句，然后在cards语句后输入两列数据，第一列表示金标准结果，第二列表示实际测量结果，这样程序就会自动绘制出相应资料的ROC曲线来。

【主要输出结果及其解释】

图5-3       与本例资料对应的ROC曲线

SAS绘制的ROC曲线不太美观，实际绘制时也可将SAS程序建立的ROC1数据集导出的结果形成EXCEL文件，然后利用EXCEL的绘图功能，将（0,0）、成对的灵敏度和1-特异度、（1,1）在坐标轴上所对应的点依次连接起来，绘制成线图，即可得到较为美观的ROC曲线图。

计算有序分类资料的ROC曲线下面积及其标准误的SAS程序如下：

%LET K=5;

DATA E;

INPUT XA1-XA&K XN1-XN&K @@;

CARDS;

1       2       11      16      44

27      18      9       1       0

;

RUN;

DATA F;

SET E;

ARRAY XN(*) XN1-XN&K;   ARRAY XA(*) XA1-XA&K;

ARRAY NN1(&K);          ARRAY NA1(&K);

ARRAY YN(&K);           ARRAY YA(&K);

ARRAY AREAS(&K);

ARRAY Q1S(&K);            ARRAY Q2S(&K);

NN=SUM(OF XN1-XN&K); NA=SUM(OF XA1-XA&K);

DO I=1 TO &K;           DO J=1 TO I;

NN1(I)+XN(J);           NA1(I)+XA(J);

YN(I)+LAG(XN(J));           /* YA(I)+LAG(XA(J));*/

END;

YA(I)=NA-NA1(I);            /* YN(I)=NN-NN1(I);*/

IF YN(I)=. THEN YN(I)=0;    /* IF YA(I)=. THEN YA(I)=0;*/

AREAS(I)=XN(I)*YA(I)+1/2*XN(I)*XA(I);

Q1S(I)=XN(I)*(YA(I)**2+XA(I)*YA(I)+1/3*XA(I)**2);

Q2S(I)=XA(I)*(YN(I)**2+XN(I)*YN(I)+1/3*XN(I)**2);

AREA=SUM(OF AREAS1-AREAS&K)/(NN*NA);

Q1=SUM(OF Q1S1-Q1S&K)/(NN*NA**2);

Q2=SUM(OF Q2S1-Q2S&K)/(NA*NN**2);

END;

SE_AREA=SQRT((AREA*(1-AREA)+(NA-1)*(Q1-AREA**2)+

(NN-1)*(Q2-AREA**2))/(NA*NN));

Z=(AREA-0.5)/SE_AREA;  P=(1-PROBNORM(Z))*2;

LCL95_A=AREA-PROBIT(0.975)*SE_AREA;

UCL95_A=AREA+PROBIT(0.975)*SE_AREA;

AREA=ROUND(AREA,0.0001);

SE_AREA=ROUND(SE_AREA,0.0001);

Z=ROUND(Z,0.0001);

P=ROUND(P,0.0001);

LCL95_A=ROUND(LCL95_A,0.0001);

UCL95_A=ROUND(UCL95_A,0.0001);

FILE PRINT;

PUT #2  @5 'ROC曲线下面积AZ=' AREA

/@5 '面积的标准误SE(AZ)=' SE_AREA

/@5 '面积与0.5差异的假设检验统计量Z值=' Z

/@5 '假设检验对应的P值=' P

/@5 '面积的95%置信区间('LCL95_A','UCL95_A')' ;

RUN;

【程序修改指导】程序第一行K=5，表明有5个分类标准，实际修改时应改为相应的分类数；在CARDS语句后，输入两行数据，第一行按类别从小到大输入异常组频数，第二行输入正常组频数；程序中部有三行包含”/*”和”*/”符号的语句，位于“/*”和“*/”中的语句被屏蔽，并不被执行。如果实际测量值越大，结果异常的概率越大，则此部分不必修改，如果实际测量值越大，结果异常的概率越小，则必须要修改此三行语句，具体把三行中左侧的语句用“/*”和“*/”屏蔽掉，而把此三行中右侧的“/*”和“*/”删掉。该程序输出ROC曲线下面积及其标准误、Z统计量及对应的P值、面积的95%置信区间等。

【主要输出结果及其解释】

ROC曲线下面积AZ=0.9636

面积的标准误SE(AZ)=0.0149

面积与0.5差异的假设检验统计量Z值=31.0567

假设检验对应的P值=0

面积的95%置信区间(0.9344 ,0.9929 )

专业结论：因ROC曲线下的面积为0.9636，此值与无效假设下的面积0.5之间的差别有统计学意义（Z=31.0567、P=0），说明此诊断方法的诊断效果是令人满意的。

上例中是人为地把连续性资料转为有序分类资料，这样做简化了计算过程，但也降低了计算的精度，因而最好还是采用原始数据进行计算，具体的SAS程序如下。

DATA ROC;

INPUT G$ NUMS;

DO I=1 TO NUMS;

INPUT VALUE@@;

OUTPUT; END;

CARDS;

异常组 74

1.21 1.30 1.41 1.63 1.54 1.69 1.65 1.62 1.68 1.61 1.52

1.67 1.59 1.63 1.84 1.78 1.89 1.76 1.83 1.78 1.72 1.87

1.75 1.78 1.72 1.83 1.86 1.78 1.83 1.76 1.87 1.88 1.88

2.04 1.92 1.93 1.99 2.00 1.90 1.97 1.94 1.89 1.96 1.95

1.87 1.91 2.03 1.93 1.97 2.02 1.93 2.01 1.86 1.91 2.02

2.00 2.03 1.87 1.95 2.02 1.99 2.04 1.89 1.96 2.00 1.88

1.99 1.96 1.89 1.91 1.97 2.02 1.98 1.87

正常组 55

1.24 1.23 1.30 1.27 1.27 1.29 1.18 1.21 1.19 1.24 1.31

1.16 1.23 1.27 1.15 1.27 1.31 1.18 1.15 1.32 1.30 1.28

1.29 1.30 1.31 1.21 1.18 1.40 1.46 1.39 1.35 1.36 1.41

1.46 1.45 1.42 1.44 1.47 1.46 1.48 1.33 1.42 1.41 1.48

1.43 1.60 1.68 1.53 1.68 1.59 1.68 1.57 1.60 1.59 1.73

;

PROC FREQ ORDER=FORMATTED PAGE;

TABLES G*VALUE /NOPRINT SPARSE OUT=A ;

DATA B;

SET A NOBS=KK; K=KK/2;  PUT K=;

IF G='正常组' THEN XNN=COUNT;  ELSE XAA=COUNT;

PROC MEANS NWAY NOPRINT; VAR XNN XAA;CLASS  VALUE;

OUTPUT OUT=C SUM=S1 S2 ;

PROC TRANSPOSE DATA=C OUT=D PREFIX=XN;   VAR S1 ;

PROC TRANSPOSE DATA=C OUT=E PREFIX=XA;   VAR S2 ;

DATA F;

SET D;  SET E;

RUN;

%LET K=65;

DATA F;

SET F;

ARRAY XN(*) XN1-XN&K;   ARRAY XA(*) XA1-XA&K;

ARRAY NN1(&K);          ARRAY NA1(&K);

ARRAY YN(&K);           ARRAY YA(&K);

ARRAY AREAS(&K);

ARRAY Q1S(&K);            ARRAY Q2S(&K);

NN=SUM(OF XN1-XN&K); NA=SUM(OF XA1-XA&K);

DO I=1 TO &K;           DO J=1 TO I;

NN1(I)+XN(J);           NA1(I)+XA(J);

YN(I)+LAG(XN(J));       /* YA(I)+LAG(XA(J));*/

END;

YA(I)=NA-NA1(I);         /* YN(I)=NN-NN1(I);*/

IF YN(I)=. THEN YN(I)=0; /*  IF YA(I)=. THEN YA(I)=0;*/

AREAS(I)=XN(I)*YA(I)+1/2*XN(I)*XA(I);

Q1S(I)=XN(I)*(YA(I)**2+XA(I)*YA(I)+1/3*XA(I)**2);

Q2S(I)=XA(I)*(YN(I)**2+XN(I)*YN(I)+1/3*XN(I)**2);

AREA=SUM(OF AREAS1-AREAS&K)/(NN*NA);

Q1=SUM(OF Q1S1-Q1S&K)/(NN*NA**2);

Q2=SUM(OF Q2S1-Q2S&K)/(NA*NN**2);

END;

SE_AREA=SQRT((AREA*(1-AREA)+(NA-1)*(Q1-AREA**2)+

(NN-1)*(Q2-AREA**2))/(NA*NN));

Z=(AREA-0.5)/SE_AREA;  P=(1-PROBNORM(Z))*2;

LCL95_A=AREA-PROBIT(0.975)*SE_AREA;

UCL95_A=AREA+PROBIT(0.975)*SE_AREA;

AREA=ROUND(AREA,0.0001);

SE_AREA=ROUND(SE_AREA,0.0001);

Z=ROUND(Z,0.0001);

P=ROUND(P,0.0001);

LCL95_A=ROUND(LCL95_A,0.0001);

UCL95_A=ROUND(UCL95_A,0.0001);

FILE PRINT;

PUT #2  @5 'ROC曲线下面积AZ=' AREA

/@5 '面积的标准误SE(AZ)=' SE_AREA

/@5 '面积与0.5差异的假设检验统计量Z值=' Z

/@5 '假设检验对应的P值=' P

/@5 '面积的95%置信区间('LCL95_A','UCL95_A')' ;

RUN;

【程序修改指导】对于连续性资料ROC曲线的计算，也可理解为对诊断界值很多的有序分类资料进行计算，此时每一个不重复的测量值都要作为诊断界值考虑，程序思路与有序分类资料相似，但需要先算出不重复的测量值个数或诊断分类个数。需将数据部分中两组的74及55个数据删掉，换为自已的数据，并将表示各组数量的“74”和“55”换为自已各组数据的个数。首先运行一遍程序，此时不必看结果，而要先从程序记录窗口中读出k值，即不重复的测量值个数，然后重新调回刚发送的程序，程序中的“%LET K=65;”中的“65”表示不重复的测量值个数，将此数值替换掉，重新运行程序，即可得到正确的结果。如果实际测量值越大，结果异常的概率越小，则还要修改程序中含有“/*”和“*/”的三行语句，具体把三行中左侧的语句用“/*”和“*/”屏蔽掉，而把此三行中右侧的“/*”和“*/”删掉。

【主要输出结果及其解释】

ROC曲线下面积AZ=0.9612

面积的标准误SE(AZ)=0.0165

面积与0.5差异的假设检验统计量Z值=27.9776

假设检验对应的P值=0

面积的95%置信区间(0.9289 ，0.9935 )

专业结论：因ROC曲线下的面积为0.9612，此值与无效假设下的面积0.5之间的差别有统计学意义（Z=27.9776、P=0），说明此诊断方法的诊断效果是令人满意的。

三、ROC曲线下面积的比较

不同的诊断系统都能获得相应的ROC曲线，ROC曲线下面积的差别是否具有统计学意义是评价各诊断系统优劣的重要标志。可以计算两面积之差及其标准误，构造出Z统计量，进行假设检验，具体计算见式（4）。

（4）

r是两个ROC曲线下面积间的相关系数，r的计算需要首先求得两个相关系数r_n和r_a，r_n和r_a分别为金标准正常组的两诊断试验间的相关系数和金标准异常组的两诊断试验的相关系数。其计算对于连续性资料可采用Pearson积差法，对于有序分类资料可采用Kendal￠s tau等级相关法。以两诊断试验平均相关(r_n+r_a)/2和平均面积(Az₁+Az₂)/2查“两个ROC曲线下面积估计值间的相关系数表”，即可得到r值。若两种方法分别测自独立样本，可令r=0，从而简化为公式（5）进行计算。

（5）

【例2】两诊断方法的Az₁、Az₂、SE₁、SE₂分别为0.9045、0.7351、0.0165、0.0423，且两诊断方法分别测自独立样本，问两个ROC曲线下面积间差异有无统计学意义？

查表得P<0.01，可以认为两个ROC曲线下面积间差异具有统计学意义。

专业结论：因第一种诊断方法所对应的曲线下的面积（0.9045）大于第二种诊断方法所对应的曲线下的面积（0.7351），且P<0.01，说明第一种诊断方法优于第二种诊断方法。

实现上述计算所需的SAS程序如下。

DATA ROCCOMP;

INPUT AZ1 AZ2 SE1 SE2;

Z=(AZ1-AZ2)/SQRT(SE1**2+SE2**2);

P=(1-PROBNORM(Z))*2;

Z=ROUND(Z,0.0001);

P=ROUND(P,0.0001);

CARDS;

0.9045 0.7351 0.0165 0.0423

;

DATA COMP;

SET ROCCOMP;

FILE PRINT;

PUT #2  @5 '统计量Z=' Z

/@5 'P值=' P;

RUN;

【程序修改指导】只需将CARDS语句后的4个数值换为自己的数据即可。

【主要输出结果及其解释】

统计量Z=3.7309

P值=0.0002

专业结论：因第一种诊断方法所对应的曲线下的面积（0.9045）大于第二种诊断方法所对应的曲线下的面积（0.7351），且P=0.0002，说明第一种诊断方法优于第二种诊断方法。