Description Usage Arguments Details Author(s) References Examples
Calculate Brier Score Based on code in pecs
1 |
obj |
|
pred |
|
btime |
Calculate Brier Score Based on code in pecs
pec
pec
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 | ##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (obj, pred, btime = range(obj[, 1]))
{
if (!inherits(obj, "Surv"))
stop("obj is not of class Surv")
class(obj) <- NULL
if (attr(obj, "type") != "right")
stop("only right-censoring allowed")
N <- nrow(obj)
time <- obj[, 1]
ot <- order(time)
cens <- obj[ot, 2]
time <- time[ot]
if (is.null(btime))
stop("btime not given")
if (length(btime) < 1)
stop("btime not given")
if (length(btime) == 2) {
if (btime[1] < min(time))
warning("btime[1] is smaller than min(time)")
if (btime[2] > max(time))
warning("btime[2] is larger than max(time)")
btime <- time[time >= btime[1] & time <= btime[2]]
}
ptype <- class(pred)
if (is.null(ptype)) {
if (is.vector(pred))
ptype <- "vector"
if (is.list(pred))
ptype <- "list"
}
if (ptype == "numeric" && is.vector(pred))
ptype <- "vector"
survs <- NULL
switch(ptype, survfit = {
survs <- getsurvprob(pred, btime)
survs <- matrix(rep(survs, N), nrow = length(btime))
}, list = {
if (!inherits(pred[[1]], "survfit")) stop("pred is not a list of survfit objects")
if (length(pred) != N) stop("pred must be of length(time)")
pred <- pred[ot]
survs <- matrix(unlist(lapply(pred, getsurvprob, times = btime)),
nrow = length(btime), ncol = N)
}, vector = {
if (length(pred) != N) stop("pred must be of length(time)")
if (length(btime) != 1) stop("cannot compute integrated Brier score with pred")
survs <- pred[ot]
}, matrix = {
if (all(dim(pred) == c(length(btime), N))) survs <- pred[,
ot] else stop("wrong dimensions of pred")
})
if (is.null(survs))
stop("unknown type of pred")
hatcdist <- prodlim(Surv(time, cens) ~ 1, reverse = TRUE)
csurv <- predict(hatcdist, times = time, type = "surv")
if (length(csurv) < length(btime)) {
addsurv = rep(0, length(btime) - length(csurv))
csurv = c(csurv, addsurv)
}
csurv[csurv == 0] <- Inf
bsc <- rep(0, length(btime))
if (length(btime) > 1) {
for (j in 1:length(btime)) {
help1 <- as.integer(time <= btime[j] & cens == 1)
help2 <- as.integer(time > btime[j])
bsc[j] <- mean((0 - survs[j, ])^2 * help1 * (1/csurv) +
(1 - survs[j, ])^2 * help2 * (1/csurv[j]))
}
idx <- 2:length(btime)
RET <- diff(btime) %*% ((bsc[idx - 1] + bsc[idx])/2)
RET <- RET/diff(range(btime))
names(RET) <- "integrated Brier score"
attr(RET, "time") <- range(btime)
}
else {
help1 <- as.integer(time <= btime & cens == 1)
help2 <- as.integer(time > btime)
cs <- predict(hatcdist, times = btime, type = "surv")
if (cs == 0)
cs <- Inf
RET <- mean((0 - survs)^2 * help1 * (1/csurv) + (1 -
survs)^2 * help2 * (1/cs))
names(RET) <- "Brier score"
attr(RET, "time") <- btime
}
RET
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.