1- brier <- function (fit , times , newdata , ties = TRUE ) {
2- # Baseline predicted probabilities, use the same approximations as the
3- # Cox model
4- if (fit $ method == " efron"
5- s0 <- survfit(fit $ y ~ 1 , se.fit = FALSE , ctype = 2 , stype = 2 )
6- else s0 <- survfit(fit $ y ~ 1 , se.fit = FALSE , stype = 1 )
1+ brier <- function (fit , times , newdata , ties = TRUE , detail = FALSE , timefix = TRUE ,
2+ efron = FALSE ) {
3+ Call <- match.call()
4+ if (! inherits(fit , " coxph" " fit must be a coxph object"
5+
6+ if (missing(newdata )) mf <- stats :: model.frame(fit )
7+ else mf <- stats :: model.frame(fit , data = newdata )
8+ Y <- mf [[1 ]] # the survival object
9+ if (! is.Surv(Y )) stop(" response must be a Surv object"
10+ type <- attr(Y , " type"
11+ if (! (type %in% c(" right" " mright" " counting" " mcounting"
12+ stop(" response must be right censored"
13+ n <- nrow(Y )
14+ ny <- ncol(Y ) # 3 = time1, time2 data
15+
16+ if (! is.logical(timefix ) || length(timefix ) > 1 )
17+ stop(" invalid value for timefix option"
18+ if (timefix ) Y <- aeqSurv(Y )
19+
20+ casewt <- model.weights(mf )
21+ if (is.null(casewt )) casewt <- rep(1 , n )
22+ else {
23+ if (! is.numeric(casewt )) stop(" weights must be numeric"
24+ if (any(! is.finite(casewt ))) stop(" weights must be finite"
25+ if (any(casewt < 0 )) stop(" weights must be non-negative"
26+ casewt <- as.numeric(casewt ) # transform integer to numeric
27+ }
28+
29+ id <- model.extract(mf , " (id)"
30+ if (ny == 3 && (is.null(id ))) stop(" id is required for start-stop data"
31+ if (! is.null(id )) {
32+ if (is.null(attr(Y , ' states'
33+ ytemp <- Y
34+ attr(ytemp , ' states' <- ' event' # survcheck2 wants a states attr
35+ check <- survcheck2(ytemp , id )
36+ }
37+ else check <- survcheck2(Y , id )
38+
39+ if (any(check $ flag > 0 ))
40+ stop(" one or more flags are >0 in survcheck"
41+ n.startstate <- sum(check $ transitions [,1 ] > 1 )
42+ if (ny == 2 ) samestart = TRUE # everyone starts at the same time
43+ else {
44+ etemp <- tapply(Y [,1 ], id , min )
45+ samestart <- all(etemp == etemp [1 ])
46+ }
47+ } else check <- NULL
48+ if (length(id )== 0 ) id <- seq.int(n )
49+
50+ # Finally, it's time to do the work.
51+ if (is.null(check ) || (n.startstate == 1 & samestart )) simple <- TRUE
52+ else simple <- FALSE
53+ if (! simple ) stop(" delayed entry is not yet implemented"
54+
55+ # For baseline predicted probabilities, use the same approximations as the
56+ # Cox model, if allowed by the 'efron' option.
57+ if (efron && fit $ method == " efron"
58+ s0 <- survfit(Y ~ 1 , weights = casewt , se.fit = FALSE , ctype = 2 , stype = 2 )
59+ else s0 <- survfit(Y ~ 1 , weights = casewt ,se.fit = FALSE , stype = 1 )
760 if (missing(times )) times <- s0 $ time [s0 $ n.event > 0 ]
861 p0 <- 1 - summary(s0 , times , extend = TRUE )$ surv
962
1063 # model predictions
11- s1 <- survfit(fit , newdata = fit $ call $ data , se.fit = FALSE )# FALSE is faster
64+ if (missing(newdata ))
65+ s1 <- survfit(fit , newdata = fit $ call $ data , se.fit = FALSE )# FALSE is faster
66+ else s1 <- survfit(fit , newdata = newdata , se.fit = FALSE )
1267 p1 <- 1 - summary(s1 , times , extend = TRUE )$ surv
1368
14- ny <- ncol(fit $ y )
15- dtime <- fit $ y [,ny - 1 ] # time and status in the data
16- dstat <- fit $ y [,ny ]
69+ dtime <- Y [,ny - 1 ] # time and status in the data
70+ dstat <- Y [,ny ]
1771 n <- length(dtime )
1872 ntime <- length(times )
1973
@@ -24,19 +78,25 @@ brier <- function(fit, times, newdata, ties=TRUE) {
2478 dtime <- dtime + ifelse(dstat == 0 , mindiff / 2 , 0 )
2579 }
2680
27- c0 <- survfit0(survfit(Surv(dtime , 1 - dstat ) ~ 1 ))
81+ c0 <- survfit0(survfit(Surv(dtime , 1 - dstat ) ~ 1 , weights = casewt ))
2882 b0 <- b1 <- matrix (0 , nrow = n , ncol = ntime )
29- wt <- rep( 1 / n , n ) # everyone starts out with a weight of 1/n
83+ casewt <- casewt / sum( casewt )
3084 brier <- matrix (0 , ntime , 2 )
85+ eff.n <- double(ntime ) # the effective n
3186 for (i in 1 : ntime ) {
3287 indx <- findInterval(pmin(dtime , times [i ]), c0 $ time )
33- wt <- ifelse(dtime < times [i ] & dstat == 0 , 0 , 1 / c0 $ surv [indx ])
88+ wt <- ifelse(dtime < times [i ] & dstat == 0 , 0 , casewt / c0 $ surv [indx ])
89+ eff.n [i ] <- 1 / sum(wt ^ 2 ) # the sum of the weights is always 1
3490
3591 b0 [,i ] <- ifelse(dtime > times [i ], p0 [i ]^ 2 , (dstat - p0 [i ])^ 2 )
3692 b1 [,i ] <- ifelse(dtime > times [i ], p1 [i ,]^ 2 , (dstat - p1 [i ,])^ 2 )
3793 brier [i ,] <- c(sum(wt * b0 [,i ]), sum(wt * b1 [,i ]))/ sum(wt )
3894 }
39-
40- dimnames(brier ) <- list (times , c(" NULL" " Model"
41- list (brier = brier , b0 = b0 , b1 = b1 , time = times )
95+ ret <- list (rsquared = 1 - brier [,2 ]/ brier [,1 ], brier = brier [,2 ], times = times )
96+ if (detail ) {
97+ ret $ p0 <- p0
98+ ret $ phat <- p1
99+ ret $ eff.n <- eff.n
100+ }
101+ ret
42102}
0 commit comments