查看原文
其他

NRI和 IDI的计算方法,再也不用为它们的计算发愁了

2017-10-06 张华 赵一鸣 Freescience联盟

本文转自临床流行病学和循证医学公众号

我们前面几期介绍过NRINet Reclassification Improvement,重分类改善指标,)和IDIIntegrated Discrimination Improvement,综合判别改善指数),两个指标均用于比较两个指标的诊断能力,一个指标比另外一个指标诊断准确率是否提高。

在诊断试验中,NRIIDI越来越受到专家的重视,很多诊断指标的比较时要求报告NRIIDI,但现行软件多不支持NRIIDI的计算,虽然根据公式可以手工计算,但计算起来还是有些麻烦,特别是对NRIIDI的假设检验,计算公式较为复杂。

我们应大家的需要,编写了两个指标计算的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);只将修改相应的变量变量名即可,发v1v2为两个检验变量,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代表阴性,否则程序会报错。

以上分享希望对大家有帮助。

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存