使用nricens包分析Cox模型NRI时如何计算p值？

你使用nricens包时已经设置了niter=1000做bootstrap抽样，包已经完成了重复抽样计算，只是默认没输出置信区间和p值，你可以通过以下两种方法获取：

方法1：利用bootstrap抽样分布直接计算 #

先把nricens的运行结果存入对象：

nri_result <- nricens(time = dat$DaysDSCtoFOreadmORlastpresc, event = dat$ReAdmFOHosp90D, 
                      mdl.std = modelSurvClin90d, mdl.new = modelSurvClinPsySoc90d, 
                      cut = c(0.15,0.25), niter = 1000, updown = 'category', t0 = 90)

结果对象里的$boot元素存储了所有bootstrap迭代的NRI相关估计值，基于此计算双侧p值：

# 提取bootstrap得到的NRI值
boot_nri <- nri_result$boot[, "NRI"]
# 计算双侧p值：统计抽样结果中偏离0的极端比例，乘以2
p_value <- 2 * min(mean(boot_nri <= 0), mean(boot_nri >= 0))

方法2：通过Z统计量计算 #

用点估计和bootstrap的标准差构建Z值，再推导p值：

# 计算bootstrap样本的标准差
nri_se <- sd(boot_nri)
# 计算Z统计量
z_score <- nri_result$est["NRI"] / nri_se
# 计算双侧p值
p_value <- 2 * pnorm(-abs(z_score))

补充说明 #

• 若要计算NRI+、NRI-或其他指标的p值，只需替换提取的bootstrap列名（比如boot_nri_plus <- nri_result$boot[, "NRI+"]）即可。
• niter=1000的设置已经足够保证bootstrap结果的稳定性，无需调整。

内容的提问来源于stack exchange，提问作者ultimac