有点坑的survIDINRI包超级bug

今天在用survIDINRI包做cox模型的预测性能指标， IDI 95%CI（综合区分改善度（Integrated Discrimination Improvement）分析的时候，一直报错！折腾了一天，搜遍google、stackoverflow、baidu、微信、知乎、简书、CSDN都没有答案，这个bug广大水友们都有留言却石沉大海。这个bug就是“

Error in unoecdf(cc, pdiff[case], Wi[case] * PTB.Vi[case]) :   NA/NaN/Inf in foreign function call (arg 5)”。

小编倒腾了一天终于下班前灵光乍现，解决了这个天坑。期间几经波澜，有的学生能跑也是运气，关了rstudio换个数据又不能跑。究竟这个超级bug怎么发生的，又怎么解决的呢？

1. survIDINRI 自带pbc数据重现

自带pbc数据为原发性硬化性胆管炎的竞争风险数据，status=2 (原发性硬化性胆管炎), status=1 (死亡)，status=0 (生存)，下列代码可顺利计算生存数据IDI NRI的95%CI。

library('survIDINRI')
library('survival')
#--- sample data (pbc in survival package) ---
D=subset(pbc,select=c("time","status","age","albumin","edema","protime","bili"))

D$status=as.numeric(D$status==2) # 先根据status的数值是否等于2，将数据转换成布尔值TRUE或者FALSE，再利用as.numeric将TRUE转换成1，FALSE转换成0。
D=D[!is.na(apply(D,1,mean)),] ; dim(D)  #！is.na去空值，apply(D,1,mean)对行取平均值用来找任意变量NA，
mydata=D[1:100,]  #选第1~100行
t0=365*5 
indata1=mydata;
indata0=mydata[,-7] ; n=nrow(D) ; #去掉倒数第7列
covs1<-as.matrix(indata1[,c(-1,-2)]) #去掉倒第1，2列
covs0<-as.matrix(indata0[,c(-1,-2)])
#让我们来看一下最后建立的模型,
covs1
covs0
#可以看到，模型covs1有5个指标，模型covs0有4个指标
#--- inference ---
x<-IDI.INF(mydata[,1:2], covs0, covs1, t0, npert=200) ;
#详解：第一个参数输入生存时间，和生存结局组成的data.frame，
#第二，第三个模型指标两个模型指标，t0时间，第三个参数，时间节点，第五个参数，置换检验的次数
#--- results ---
IDI.INF.OUT(x) 
#可以看计算结果
#M1:IDI
#M2：NRI
#M3：中值改善度（median improvement)
#第一行，估计值，第二行，第三行2.5% 97.5%  置信区间，最后一行，p值
#--- Graphical presentaion of the estimates ---
IDI.INF.GRAPH(x) 
#根据计算结果输出图片

2.Bug解决

在应用自己的数据时候，发现考虑了status01编码、数据tibble、缺失值NA问题后，依旧出现恶心的error提示。那么这个bug究竟怎么回事呢？答案是

rm(list = ls())
setwd("D:\\聂个人文件\\我的公众号\\IDI")
# setwd("C:/Users/83718/Desktop/肺功能与心衰/数据整理与统计分析/第二版统计分析")
pacman::p_load(haven,survIDINRI,survival)
#tb格式 
library(tibble)  
dat<-read_sav("6.sav")
dat=as_tibble(dat)
#status，2=event事件，1=HF，0=live，
dat$HF_final=ifelse(dat$HF_final==1,1, 0)
str(dat)
#精简数据集，93161
#-----方法1 无缺失空
dat <- subset(dat, select = c(HF_scxf_timecox,HF_final,Age,Female,smoke_status,IS_HR))
dat <- na.omit(dat)
#-----方法2 无缺失空
# dat <- dat[!is.na(apply(dat,1,mean)),]
#全变量indata1，精简-1 var indata0


#---计算方法1  模型变量递增， IDI
#ts也必须是time，status 顺序，ts必须是matrix！不然容易出错
ts = as.matrix(subset(dat, select = c(HF_scxf_timecox,HF_final)))
#参考原模型，协变量矩阵
old = as.matrix(subset(dat, select = c(Age,Female,smoke_status)))
#新模型新增HR，协变量矩阵
new = as.matrix(subset(dat, select = c(Age,Female,smoke_status,IS_HR)))
res<-IDI.INF(ts,old,new,t0=4,npert = 20,seed1 =1234) 
IDI.INF.OUT(res)
##t0为设置的时间。npert设置迭代次数。
#m1：IDI的值，可信区间，P值；m2：NRI的值，可信区间，P值；m3：Median improvement in risk score，可信区间，p值。

#-------ts 不是matrix报错！
#Error in unoecdf(cc, pdiff[case], Wi[case] * PTB.Vi[case]) : 
#  NA/NaN/Inf in foreign function call (arg 5)

#---计算方法2  模型变量递减，IDI
indata1=dat
indata0=dat[,-6] 
##旧模型中预测变量矩阵,ts也必须是time，status 顺序
new<-as.matrix(indata1[,c(-1,-2)])
old<-as.matrix(indata0[,c(-1,-2)])
ts <- as.matrix(dat[,1:2])
res<-IDI.INF(ts,old, new, t0=4, npert=20,seed1 =1234)
IDI.INF.OUT(res)

#---计算方法2 发散模型不加matrix
new<- indata1[,c(-1,-2)] 
old<- dat[,c(-1,-2)] 
ts <- as.matrix(dat[,1:2])
res<-IDI.INF(ts,old, new, t0=4, npert=20,seed1 =1234)
IDI.INF.OUT(res)

注意几处细节，小编标记如下。

小结

原创不易，转载请说明来自本公众号。