isbrier: Calculate Brier Score Based on code in pecs

Description Usage Arguments Details Author(s) References Examples

View source: R/isbrier.R

Description

Calculate Brier Score Based on code in pecs

Usage

1
isbrier(obj, pred, btime = range(obj[, 1]))

Arguments

obj
pred
btime

Details

Calculate Brier Score Based on code in pecs

Author(s)

pec

References

pec

Examples

 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
  }

whcsu/rsfse documentation built on Dec. 4, 2019, 2:10 p.m.