Блог

Я хотел бы (и не смог) пересчитать взвешенную оценку шиповника (уравнение (3) и приведенное ниже в этой статье), рассчитанную с помощью пакета pec.

Вот мои попытки до сих пор:

 set.seed(130971)
library(dplyr)
library(survival)
library(prodlim) ## for SimSurv()
library(pec)


## generate simulated data
dat <- SimSurv(100)

## compute the apparent prediction error
Models <- list("Cox.X1"=coxph(Surv(time,status)~X1,data=dat,y=TRUE))
PredError <- pec(object=Models,
                 formula=Surv(time,status)~X1 X2,
                 data=dat,
                 verbose=TRUE)

## I'd like to re-calc the brier score at time 5
mytime <- PredError$time[58]
pecbs <- PredError$AppErr$Cox.X1[58]
mytime
## [1] 5.072562
pecbs
## [1] 0.2159007

## so, I need the weights for the weighted brier score
dat <- dat %>%
  mutate(ipcw5 = ipcw(
        ##formula=Surv(time,status)~X1 X2,
        formula=Surv(time,status)~1,
        data = dat,
        ##method = "cox",
        method = "marginal",
        what = "IPCW.times",
        times = mytime)$IPCW.times)

## I also need the predicted survival probabilities at time = 5
coxmod <- coxph(Surv(time,status)~X1,data=dat,y=TRUE)
dat <- dat %>%
  mutate(risk5 = as.numeric(1-summary(survfit(
                                coxmod,
                                newdata = dat,
                                se.fit = F,
                                conf.int = F
                              ),
                              times = mytime)$surv))

## and finally, I need the status at time = 5
dat <- dat %>%
  mutate(status5 = ifelse(status == 1 amp; time <= mytime, 1,
                   ifelse(status == 0 amp; time  > mytime, 0,
                   ifelse(status == 1 amp; time  > mytime, 0,
                   ifelse(status == 0 amp; time <= mytime, 0,  ## this is censored
                          NA)))))

## brier score
mybs <- dat %>%
  summarize(bs = sum((status5-risk5)^2*ipcw5)/n())

## compare
pecbs
## [1] 0.2159007
mybs$bs
## [1] 0.1612572
pecbs - mybs$bs
## [1] 0.05464348
  

Итак, я недооцениваю показатель шиповника, когда вычисляю его вручную (в mybs ).

Мой вопрос: как я могу правильно рассчитать взвешенный показатель шиповника в некоторый момент времени t (правильно = как в pec )?