NRI和 IDI的计算方法,再也不用为它们的计算发愁了
本文转自临床流行病学和循证医学公众号
我们前面几期介绍过NRI(Net Reclassification Improvement,重分类改善指标,)和IDI(Integrated Discrimination Improvement,综合判别改善指数),两个指标均用于比较两个指标的诊断能力,一个指标比另外一个指标诊断准确率是否提高。
在诊断试验中,NRI和IDI越来越受到专家的重视,很多诊断指标的比较时要求报告NRI和IDI,但现行软件多不支持NRI和IDI的计算,虽然根据公式可以手工计算,但计算起来还是有些麻烦,特别是对NRI和IDI的假设检验,计算公式较为复杂。
我们应大家的需要,编写了两个指标计算的R程序。在此介绍给大家,供有需要者使用。
我们知道NRI及其检验z值的计算公式是:
对应的R程序为:
NRIcalculate=function(m1="dia1",m2="dia2",gold="gold"){
datanri=datanri[complete.cases(datanri),];
for (i in 1:length(names(datanri))){
if (names(datanri)[i]==m1)nm1=as.numeric(i);
if (names(datanri)[i]==m2)nm2=as.numeric(i);
if(names(datanri)[i]==gold)ngold=as.numeric(i);
};
if(names(table(datanri[,nm1]))[1]!="0" ||
names(table(datanri[,nm1]))[2]!="1")stop("指标1诊断值不是0和1");
if(names(table(datanri[,nm2]))[1]!="0" ||
names(table(datanri[,nm2]))[2]!="1")stop("指标2诊断值不是0和1");
if(names(table(datanri[,ngold]))[1]!="0" ||
names(table(datanri[,ngold]))[2]!="1")stop("金标准诊断值不是0和1");
datanri1=datanri[datanri[,ngold]==1,]
table1=table(datanri1[,nm1],datanri1[,nm2]);
datanri2=datanri[datanri[,ngold]==0,]
table2=table(datanri2[,nm1],datanri2[,nm2]);
p1=as.numeric(table1[2,1]/table(datanri[,ngold])[2]);
p2=as.numeric(table1[1,2]/table(datanri[,ngold])[2]);
p3=as.numeric(table2[2,1]/table(datanri[,ngold])[1]);
p4=as.numeric(table2[1,2]/table(datanri[,ngold])[1]);
NRI=round(p1-p2-p3+p4,3);
z=NRI/sqrt((p1+p2)/table(datanri[,ngold])[2]+(p3+p4)/table(datanri[,ngold])[1]);
z=round(as.numeric(z),3);
pvalue=round((1-pnorm(abs(z)))*2,3);
if(pvalue<0.001)pvalue="<0.001";
result=paste("NRI=",NRI,",z=",z,",p=",pvalue,sep= "");
return(result)
}
#以上程序不需要更改
library(foreign);
datanri=as.data.frame(read.spss("D:\\datanri.sav")); #请更改数据存在路径
NRIcalculate(m1="v1",m2="v2",gold="gold");#m1为第一诊断变量名,m2第二诊断变量名,gold为金标准
看不懂上面的程序怎么办,不用着急,您复制上面的程序,做一些简单的修改即可。
程序先写了一个函数,将数据类型判断和计算全部写进函数,不需要使用者做任何修改。
倒数第3行,加载foreign程序包,为了导入其它类型的数据到R里;
倒数第2行,导入数据,如果是SPSS数据,直接修改SPSS数据的存放路径即可,其它类型数据请选择适当的方法导入;
倒数第1行,使用函数,计算NRI及检验检验(检验NRI是否等于0);只将修改相应的变量变量名即可,发v1和v2为两个检验变量,gold为金标准的变量名;
注:变量名必须为英文或英文+数字;三个变量必须为二分类变量,均用1代表阳性,0代表阴性,否则程序会报错。
例如下面这个数据,将数据保存在D盘根目录下,命名文件名为datanri.sav,变量名为为v1,v2,gold;
使用上方程序运行后,结果如下:NRI=0.42,z=1.787,p=0.074。
同样的方法,我们使用IDI其及z值的计算公式编写程序。
程序如下:
IDIcalculate=function(m1="v1",m2="v2",gold="gold"){
dataidi= dataidi [complete.cases(dataidi),];
for (i in 1:length(names(dataidi))){
if(names(dataidi)[i]==m1)nm1=as.numeric(i);
if(names(dataidi)[i]==m2)nm2=as.numeric(i);
if(names(dataidi)[i]==gold)ngold=as.numeric(i);
};
if(names(table(dataidi[,ngold]))[1]!="0" ||
names(table(dataidi[,ngold]))[2]!="1")stop("金标准诊断值不是0和1");
logit1=glm(dataidi[,ngold]~dataidi[,nm1],family=binomial(link='logit'),data=dataidi)
dataidi$pre1=logit1$fitted.values;
logit2=glm(dataidi[,ngold]~dataidi[,nm2],family=binomial(link='logit'),data=dataidi)
dataidi$pre2=logit2$fitted.values;
dataidi$predif=dataidi$pre1-dataidi$pre2;
dataidi1=dataidi[dataidi[,ngold]==1,];
dataidi2=dataidi[dataidi[,ngold]==0,];
p1=mean(dataidi1$pre1);
p2=mean(dataidi1$pre2);
p3=mean(dataidi2$pre1);
p4=mean(dataidi2$pre2);
IDI=round(p1-p2-p3+p4,3);
z=IDI/sqrt(sd(dataidi1$predif)/length(dataidi1$predif)+sd(dataidi2$predif)/length(dataidi2$predif));
z=round(as.numeric(z),3);
pvalue=round((1-pnorm(abs(z)))*2,3);
if(pvalue<0.001)pvalue="<0.001";
result=paste("IDI=",IDI,",z=",z,",p=",pvalue,sep= "");
return(result)
}
#以上程序不需要更改
library(foreign);
dataidi=as.data.frame(read.spss("D:\\dataidi.sav"));#请更改数据存在路径
IDIcalculate(m1="v1",m2="v2",gold="gold");#m1为第一诊断变量名,m2第二诊断变量名,gold为金标准
程序与NRI计算程序相似,不作解释。需要注意的是两个检验变量可以为连续变量或二分类变量,金标准变量必须为二分类变量,用1代表阳性,0代表阴性,否则程序会报错。
以上分享希望对大家有帮助。