IBS greater than 1 issue

[thebasheersj]

I have tried SurvMetrics on many datasets in my work, but in two datasets, It gives me IBS value greater than 1 which can't be so.
The code where I get wrong IBS is below for verification purpose.

`# packages
library(survival)
library(randomForestSRC)
library(glmnet)
library(caTools)
library(tidyverse)
library(SurvMetrics)

data(peakVO2, package = "randomForestSRC")

Data splitting

set.seed(2345)
peak.sample <- sample.split(peakVO2$ttodead, 0.80)
tr.peak<- subset(peakVO2, peak.sample == T)
ts.peak<- subset(peakVO2, peak.sample == F)
dim(ts.peak)

Fit RSF model

set.seed(23456)
peak.rsf<- rfsrc(Surv(ttodead, died)~., ntree = 1000, importance = T, statistics = T, tr.peak)
peak.rsf

peak.rsf.pred<- predict(peak.rsf, ts.peak)$survival

peakful_time<- peak.rsf$time.interest
peakful_time

surv_peakful = Surv(ts.peak$ttodead, ts.peak$died)
surv_peakful

set.seed(75943038)
peakful_IBS<- IBS(surv_peakful, sp_matrix = peak.rsf.pred, peakful_time)
peakful_IBS

Result: 1.203971`

[LauricF]

I had a similar problem and the Brier function would return a value of 134 in my dataset.

I think the mistake comes from the Brier function,

particularly between the line 76 to 92;

The current code is:
Gtstar <- Gt(object, t_star)
for (i in c(1:length(time))) {
if (time[i] < t_star & (status[i] == 1)) {
Gti <- Gt(Surv(time, status), time[i])
if (is.na(Gti)) {
next
}
sum_before_t <- sum_before_t + 1/Gti * (pre_sp[i])^2
next
}
if (time[i] >= t_star) {
if (is.na(Gtstar)) {
next
}
sum_after_t <- sum_after_t + 1/Gt(Surv(time, status),
t_star) * (1 - pre_sp[i])^2
}
}

For a reason that I cannot yet understand Gt(Surv(time, status), t_star) return a different result that Gt(object, t_star). Gt(Surv(time, status), t_star) return a very small value which inflate the Brier score and so the IBS too.
When I corrected the code with the use of Gt(object, t_star) instead of Gt(Surv(time, status), t_star) I had a meaningful result below 1.

corrected code
Gtstar <- Gt(object, t_star)
for (i in c(1:length(time))) {
if (time[i] < t_star & (status[i] == 1)) {
Gti <- Gt(object, time[i])
if (is.na(Gti)) {
next
}
sum_before_t <- sum_before_t + 1/Gti * (pre_sp[i])^2
next
}
if (time[i] >= t_star) {
if (is.na(Gtstar)) {
next
}
sum_after_t <- sum_after_t + 1/Gtstar * (1 - pre_sp[i])^2
}
}

[anlo328]

Hello, any update on this? I'm having the same issue; I gave the exact same data with the same integration interval to two different implementations and got:

SurvMetrics::IBS -> 4.121665
survex::integrated_brier_score -> 0.1830576

Seems the second makes more sense

[comisora]

Gt() uses row-indices from the deduplicated KM table to subscript the raw (non-deduplicated) time vector, producing invalid Brier/IBS values (>1) when event/censoring times contain ties

Summary

Brier() and IBS() can return values greater than 1, which is mathematically
impossible for a proper Brier score (it is bounded in [0, 1] by construction).
This happens whenever the survival data contains tied times, which is the
common case for any data recorded on a discrete scale (days, weeks, visit
numbers, etc.) rather than as continuous, distinct timestamps.

Root cause

In Gt(), when the requested timepoint falls between two points of the
(deduplicated) KM summary table, the function does this:

res.sum <- survminer::surv_summary(fit)
res.sum <- na.omit(res.sum)
...
index1 <- ifelse(length(which(res.sum$time < timepoint)),
                  max(which(res.sum$time < timepoint)), 1)
index2 <- ifelse(length(which(res.sum$time > timepoint)),
                  min(which(res.sum$time > timepoint)), length(res.sum$time))
Gtemp <- res.sum$surv
if (index1 == index2) {
  Gvalue <- Gtemp[index1]
} else {
  Gvalue <- ((time[index2] - timepoint) * Gtemp[index1] +
             (timepoint - time[index1]) * Gtemp[index2]) /
            (time[index2] - time[index1])
}

index1/index2 are row positions inside res.sum (the deduplicated KM
table, one row per unique event/censoring time). But the interpolation step
then uses time[index1] and time[index2] — where time is the raw,
per-subject vector (object[, 1]), which has a different length and a
different value at each index whenever there are ties (so
length(unique(time)) < length(time)).

Concretely: index1/index2 are valid positions in res.sum (length =
number of unique times), but are used to index time (length = number of
subjects). As soon as any tie has occurred earlier in the timeline, position
k in res.sum no longer corresponds to position k in the sorted
per-subject time vector — the two vectors drift apart. Whenever this drift
makes time[index1] and time[index2] coincide (or otherwise land on the
wrong values), the interpolation denominator time[index2] - time[index1]
becomes 0 (or the numerator/denominator combination becomes nonsensical),
producing NaN/Inf. The subsequent safety check

if (is.na(Gvalue) | Gvalue <= 0) {
  Gvalue <- min(res.sum$surv[res.sum$surv != 0])
}

then silently replaces the broken value with the smallest non-zero survival
probability found anywhere in the whole censoring KM table — which for an
early query time can be a value belonging to the very end of follow-up,
i.e. a tiny, completely unrelated number. This is why Gt() doesn't error
out or return NA, but instead returns a small-but-plausible-looking wrong
number, which then propagates into Brier()'s 1/Gti weight, inflating the
score, sometimes past the theoretical upper bound of 1.

Reproducible example

library(survival)
library(SurvMetrics)

set.seed(123)
n <- 208

## Right-skewed event times with heavy early ties (e.g. recorded in whole
## days) + censoring starting only after day 90 -- a very ordinary shape
## for clinical time-to-event data.
event_time <- pmax(1, round(rexp(n, 1 / 150)))
is_event   <- rbinom(n, 1, 0.55)
cens_time  <- pmax(90, round(rexp(n, 1 / 500)))

time   <- ifelse(is_event == 1, pmin(event_time, 1800), cens_time)
status <- ifelse(is_event == 1 & event_time <= 1800, 1, 0)
status[is_event == 1 & event_time > cens_time] <- 0
time[is_event == 1 & event_time > cens_time]   <- cens_time[is_event == 1 & event_time > cens_time]

obj <- Surv(time, status)

km <- survfit(Surv(time, 1 - status) ~ 1)
rs <- survminer::surv_summary(km)
cat("n subjects:", n, " | unique times in KM table:", nrow(rs), "\n")
#> n subjects: 208  | unique times in KM table: 152      <-- ties present

## Any predicted survival matrix will do -- use the same monotone decreasing
## curve for every subject, purely to isolate the Gt()/Brier() bug.
grid <- seq(1, 365, 7)
sp   <- matrix(rep(seq(0.95, 0.5, length.out = length(grid)), each = n), nrow = n)

briers <- sapply(seq_along(grid), function(j) Brier(obj, sp[, j], t_star = grid[j]))

data.frame(t = grid, brier = round(briers, 3))[briers > 1, ]
#>      t  brier
#> ...  85  1.900
#> ...  99  1.818
#> ... 113  2.039
#> ... 169  3.050
#> ... 344  5.131
#> ... 351  5.182
#> ... 358  5.338
## 19 of 53 evaluated time points give a Brier score above 1 -- impossible
## by definition. max(briers) = 5.34 here.

You can inspect the broken value directly:

SurvMetrics:::Gt(obj, 344)   # returns a tiny, wrong value
summary(km, times = 344, extend = TRUE)$surv  # the actual, correct G(344)

On this seed, Gt(obj, 344) and the two other affected points all collapse
to min(res.sum$surv[res.sum$surv != 0]) — the fallback branch described
above — while summary(km, times = 344)$surv gives the value that should
have been returned.

Suggested fix

Use res.sum$time (matching the vector index1/index2 were computed
against) for the interpolation boundaries instead of the raw time vector:

} else {
  t1 <- res.sum$time[index1]
  t2 <- res.sum$time[index2]
  Gvalue <- ((t2 - timepoint) * Gtemp[index1] +
             (timepoint - t1) * Gtemp[index2]) / (t2 - t1)
}

More robust alternative: since Gt() is meant to evaluate a step function
(the KM estimate of the censoring distribution), linear interpolation
between KM points isn't standard practice anyway (Graf et al. 1999, and
essentially every other implementation, e.g. pec::predictSurvProb,
riskRegression, use the right-continuous step-function value, not a
linear interpolation). Replacing the whole else branch with a simple
step-function lookup would both fix this bug and match the conventional
definition:

Gvalue <- summary(fit, times = timepoint, extend = TRUE)$surv

which is a one-line replacement for the entire index1/index2/interpolation
block, and cannot suffer from an index/length mismatch since it queries the
survfit object directly rather than manually walking surv_summary()'s
output.