CODEWITHSUNDEEP
rpredictionsurvival-analysiscox-regression
tomka
tomka

Reputation: 2638

Estimating prediction accuracy of a Cox survival model using sbrier (R)

The integrated Brier score (IBS) has been suggested in a paper by Graf et al (1999) as a good measure for prediction accuracy in survival models (see e.g. overview paper by Wiering et al., page 23).

It was implemented in the package ipred as function sbrier. However, whereas the brier score definition obviously applies to Cox proportional hazard models, I cannot get sbrier to return the Brier score for a coxph model.

Here is the problem set up.

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

Now I want to estimate the IBS. Following ?sbrier 

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

So we have a list of survfit objects

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

or survival probabilties

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

The first returning

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

The second

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

The examples do not list a coxph model. Perhaps it's not supported. Otherwise, where do I go wrong?

Upvotes: 1

Views: 2601

Answers (2)

Yimin
Yimin

Reputation: 11

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)

Only 38 observations were used . The others needed to be removed before predicting.

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

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

result

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

The btime must be specified, which means ibs can't be calculated. I have not figured it out yet.

Upvotes: 1

PhilippPro
PhilippPro

Reputation: 698

You can use the pec package, instead.

Example: 

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)

Upvotes: 2

Related Questions