32 lines
1.1 KiB
R
Raw Permalink Normal View History

2025-01-12 00:52:51 +08:00
# Tests of the Brier score.
# Start with the example in the vignette
library(survival)
rott2 <- rotterdam
ignore <- with(rott2, recur ==0 & death==1 & rtime < dtime)
rott2$rfs <- with(rott2, ifelse(recur==1 | ignore, recur, death))
rott2$rfstime <- with(rott2, ifelse(recur==1 | ignore, rtime, dtime))/365.25
rsurv <- survfit(Surv(rfstime, rfs) ~1, rott2) #KM
rfit <- coxph(Surv(rfstime, rfs) ~ pspline(age) + meno + size + pmin(nodes,12),
rott2)
tau <- c(2,4,6, 8) # four tau values
bfit <- brier(rfit, times=tau)
# Now by hand
wtmat <- rttright(Surv(rfstime, rfs) ~ 1, rott2, times=tau)
psurv <- survfit(rfit, newdata= rott2) # one curve per subject
yhat <- 1- summary(psurv, times=tau)$surv
ybar <- 1- summary(rsurv, times=tau)$surv
y <- with(rott2, cbind(rfstime <=tau[1] & rfs==1,
rfstime <=tau[2] & rfs==1,
rfstime <=tau[3] & rfs==1,
rfstime <=tau[4] & rfs==1)) * 1L
ss1 <- colSums(wtmat * (y - t(yhat))^2)
ss2 <- colSums(wtmat * (y - rep(ybar, each=nrow(y)))^2)
all.equal(unname(1- ss1/ss2), bfit$rsquared)
all.equal(unname(ss1), bfit$brier)