使用 sbrier (R) 估计 Cox 生存模型的预测准确性答案

【问题标题】：Estimating prediction accuracy of a Cox survival model using sbrier (R)使用 sbrier (R) 估计 Cox 生存模型的预测准确性
【发布时间】：2017-03-03 15:46:33
【问题描述】：

Graf 等人 (1999) 在一篇论文中建议使用综合 Brier 评分 (IBS) 作为生存模型中预测准确性的良好衡量标准（例如，参见 overview paper by Wiering et al.，第 23 页）。

它在包ipred 中实现为函数sbrier。然而，虽然 Brier 分数定义显然适用于 Cox 比例风险模型，但我无法让 sbrier 返回 coxph 模型的 Brier 分数。

这是设置的问题。

library(survival)
library(ipred)
data("DLBCL", package = "ipred")

#Fit coxph model    
smod   <- Surv(DLBCL$time, DLBCL$cens)
coxmod <- coxph(smod ~ IPI, data = DLBCL) # I just chose a significant covariate from DLBCL

现在我想估算 IBS。关注?sbrier

obj  : an object of class Surv.
pred : predicted values. Either a probability or a list of survfit objects.

所以我们有一个 survfit 对象列表

sbrier(smod, list(survfit(coxmod) ))

或生存概率

sbrier(smod, survfit(coxmod,newdata=DLBCL)$surv )

第一个返回

Error in sbrier(smod, list(survfit(coxmod))) : 
  pred must be of length(time)

第二个

Error in sbrier(smod, survfit(coxmod, newdata = DLBCL)$surv) : 
  wrong dimensions of pred

示例未列出coxph 模型。也许它不受支持。不然我哪里错了？

【问题讨论】：

标签： r prediction survival-analysis cox-regression

【解决方案1】：

您可以改用 pec 包。

例子：

library(pec)
set.seed(18713)
library(prodlim)
library(survival)
dat=SimSurv(100)
pmodel=coxph(Surv(time,status)~X1+X2,data=dat)
perror=pec(list(Cox=pmodel),Hist(time,status)~1,data=dat)
## cumulative prediction error
crps(perror) # between min time and 1
## same thing:
ibs(perror)

【讨论】：

【解决方案2】：

library(survival)
library(ipred)
data("DLBCL", package = "ipred")



smod   <- Surv(DLBCL$time, DLBCL$cens)
coxmod <- coxph(smod ~ IPI, data = DLBCL)
coxmod
Call:
coxph(formula = smod ~ IPI, data = DLBCL)
     coef exp(coef) se(coef)    z      p
IPI 0.505     1.657    0.181 2.79 0.0053
Likelihood ratio test=8.15  on 1 df, p=0.0043
n= 38, number of events= 20 
   (2 observations deleted due to missingness)

只使用了 38 个观察值。在预测之前需要删除其他的。

DLBCL <- DLBCL[!(is.na(DLBCL$IPI)|is.na(DLBCL$time)|is.na(DLBCL$cens)),]

pred  <- predict(coxmod)
 sbrier(smod, pred,btime=2)

结果

Brier score 
  1.457241 
attr(,"time")
[1] 2

必须指定btime，即无法计算ibs。我还没有弄清楚。

【讨论】：