# Inconsistent results calculating the integrated brier score in R

I would like to calculate the integrated brier score as a measure of model performance for a cox survival model I am fitting.

There are multiple packages and functions to do this:

pec (pec and crps)

Initially I have simulated some fake data to test the functions out. However I get inconsistent results when using the different packages. Code to produce the estimates (from a reproducible example) inside the following table can be found at the bottom of this question. In my simulated data I try to calculate

(1) The brier score (BS) at time = 30

(2) The integrated brier score (IBS) from time = 1 to time = 30

| Package | BS | IBS |

| survcomp | 0.242 | 0.0575 |

| survAUC | 0.242 | 0.0917 |

| pec | NA | 1.53 |

Note the pec function crps only does integrated brier scores.

The functions in survAUC and survcomp both calculate the same brier score at the time point = 30, however the integrated brier scores (which I want) do not match up. In survcomp, there is no option to choose a time at which to integrate up to, so I hypothesise maybe it is integrated over all survival times? But then it makes no sense that the integral should be smaller than the other two integrals, which both integrate up to time point 30. The two functions which can integrate up to a specified time point (from survAUV and pec), give very different answers.

One hypothesis is that I should be using predicted probabilities instead of the values of the linear predictor, or I am implementing the software incorrectly.

### Load relevant packages:

source("https://bioconductor.org/biocLite.R")
biocLite("survcomp")
library(survcomp)

install.packages("survival")
library(survival)

install.packages("pec")
library(pec)

install.packages("survAUC")
library(survAUC)

## This function creates mixed effects survival data, with fixed effects beta, and a random effect from sigma
simulWeib <- function(N, lambda, rho, beta1, beta2, beta3, beta4, rateC, sigma)
{
# covariate --> N Bernoulli trials
x1 <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5))
x2 <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5))
x3 <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5))
x4 <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5))

# Now create random effect stuff
# Create one vector of length N, all drawn from same normal distribution
rand.effect <- rnorm(N,0,sigma)

# Weibull latent event times
v <- runif(n=N)
Tlat <- ceiling((- log(v) / (lambda * exp(x1 * beta1 + x2 * beta2 + x3 * beta3 + x4 * beta4 + rand.effect)))^(1 / rho))

# censoring times
#C <-rep(100000,N)
C <- rexp(n=N, rate=rateC)

# follow-up times and event indicators
time <- pmin(Tlat, C)
#status <- as.numeric(rep(1,N))
status <- as.numeric(Tlat <= C)

# data set
data.frame(id=1:N,
time=time,
status=status,
x1 = x1,
x2 = x2,
x3 = x3,
x4 = x4)
}

set.seed(10001)

## Generate the data, it has no random effects, and three fixed effects, and the censoring distribution has such a low rate that no observations are censored
data.sim<-simulWeib(10000,lambda=0.005,rho=1.3,beta1=0.25,beta2=0.33,beta3=0.125,beta4=0,rateC=0.0000000001, sigma = 0)

## Fit the correct model
model<-coxph(Surv(time,status) ~ x1 + x2 +x3, data=data.sim)
summary(model)

## Create risk scores (these are not probabilities, but are the linear predictor centered around the mean of the covarites)
pred.model<-predict(model, newdata=data.sim)

### 1 Get BS and IBS from survcomp function sbrier.score2proba
dd<-data.frame("time"=data.sim$time,"event"=data.sim$status,"score"=(pred.model))
BS<-sbrier.score2proba(data.tr=dd, data.ts = dd, method = "cox")
BS$bsc[BS$time==30]
BS$bsc.integrated ### 2 Get BS and IBS from survAUC function predErr BS2<-predErr(Surv.rsp=Surv(data.sim$time,data.sim$status),Surv.rsp.new=Surv(data.sim$time,data.sim$status),lp=pred.model, lpnew=pred.model, times=1:30,type="brier",int.type='unweighted') BS2 ### 3 Get BS and IBS from pec functions pec and crps pec.obj<-pec(model,formula=Surv(data.sim$time,data.sim$status)~x1+x2+x3,data=data.sim) BS3<-crps(pec.obj,times=c(30),start=c(1)) BS3 ### Results: BS$bsc[BS$time==30] BS$bsc.integrated
BS2
BS3