casebase Brier score Development

Question

casebase Brier score Development

Closed this issue 4 years ago · 6 comments

The first issue I wanted to get feedback on, is what method we should use to calculate the censoring probabilities. You will see the difference between using KM and using Weibull below. It is minimal between them, but there is a significant difference when comparing to a score from Cox. I'm not sure what would be ideal in this scenario. My gut says Weibull since its comparing apples to apples, but then we may lose the comparison between parametric and semi-parametric. Beyond that though, the brier score from a cox model may always perform better than a parametric one, from the looks of these results.

NOTE: this is before the choice in null model for IPA, which I will match to the censoring distribution's model.

library(casebase)
#> See example usage at http://sahirbhatnagar.com/casebase/
library(survival)
library(riskRegression)
#> riskRegression version 2019.11.03
library(reprex)
library(reshape2)
library(ggplot2)
library(prodlim)
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following objects are masked from 'package:reshape2':
#> 
#>     dcast, melt


brierScoreKMCens <- function(predSurvs, times, newData) {
  # inverse probability censoring weights
  # Note that we take probabilities right before the drops.
  
  #DIFFERENCES START HERE
  fitCens <- prodlim::prodlim(Hist(time, event != 0) ~ 1, newData, reverse = TRUE)
  IPCW.subject.times <- prodlim::predictSurvIndividual(fitCens, lag = 1) # G(t-|X)
  #DIFFERENCES END HERE
  
  # Empty matrix that will be filled in with the following loop
  Score <- matrix(NA, nrow(predSurvs), ncol(predSurvs))
  IPCWMatrix <- matrix(NA, nrow(predSurvs), ncol(predSurvs))
  # for each point in time we have predicted
  for (i in 1:length(times)) {
    # get number of censored individuals so long as their survival time is less than times[i]
    # these individuals do not have an effect on right-censored brier score.
    CensBefore <- newData$event == 0 & newData$time < times[i]
    # y encompasses all survival times larger than t[i] with a 1, 0 otherwise
    y <- drop(t(newData$time > times[i]))
    # above permits the two parts of right-censored brier score to be calculated, without IPCW, in one line
    Score[i, ] <- (y - predSurvs[i, ])^2

    # Generate IPCWMatrix
    IPCWMatrix[i, y == 0] <- IPCW.subject.times[y == 0] # G(t-|X) filled in corresponding positions
    #DIFFERENCES START HERE
    IPCW.time <- predict(fitCens, newdata = newData, times = times[i], level.chaos = 1, mode = "matrix", type = "surv", lag = 1)
    #DIFFERENCES END HERE
    IPCWMatrix[i, y == 1] <- IPCW.time # G(t) filled, same value, for remaining positions.
    # above calculated individuals who don't have an effect on brier score at specfied time. they're scores are set to 0.
    Score[i, CensBefore] <- 0
  }

  # apply IPCW to all scores
  Err <- Score / IPCWMatrix

  # Average curve demonstrating right-censored brier averaged over test-set, for each time of interest.
  Err <- apply(Err, 1, mean)
  return(Err)
}

brierScoreweibCens <- function(predSurvs= absRiskcb, times=times, newDataX=newDataX,newDataY=newDataY) {
  # inverse probability censoring weights
  
  #DIFFERENCES START HERE
  yCens <- newDataY
  yCens[,1] <- yCens[,1]==0
  fitCens <- fitSmoothHazard(event~log(time) ,data= as.data.frame(yCens), time="time",ratio=100)
  IPCW.subject.times <- absoluteRisk(object = fitCens, time = times, s = "lambda.1se")
  IPCW.subject.times[, -c(1)] <- 1 - IPCW.subject.times[, -c(1)] # make it survival probabilities
  IPCW.subject.times <- rowMeans(IPCW.subject.times[-c(1), -c(1)]) # remove extra time at 0. Also remove the first column (times) as it matches perfectly now
  #DIFFERENCES END HERE
  
  # Empty matrix that will be filled in with the following loop
  Score <- matrix(NA, nrow(predSurvs), ncol(predSurvs))
  IPCWMatrix <- matrix(NA, nrow(predSurvs), ncol(predSurvs))
  # for each point in time we have predicted
  for (i in 1:length(times)) {
    # get number of censored individuals so long as their survival time is less than times[i]
    # these individuals do not have an effect on right-censored brier score.
    CensBefore <- newData$event == 0 & newData$time < times[i]
    # y encompasses all survival times larger than t[i] with a 1, 0 otherwise
    y <- drop(t(newData$time > times[i]))
    # above permits the two parts of right-censored brier score to be calculated, without IPCW, in one line
    Score[i, ] <- (y - predSurvs[i, ])^2
    
    # Generate IPCWMatrix
    IPCWMatrix[i, y == 0] <- IPCW.subject.times[y == 0] # G(t-|X) filled in corresponding positions
    #DIFFERENCES START HERE
    IPCW.time <- IPCW.subject.times[i]
    #DIFFERENCES END HERE
    IPCWMatrix[i, y == 1] <- IPCW.time # G(t) filled, same value, for remaining positions.
    # above calculated individuals who don't have an effect on brier score at specfied time. they're scores are set to 0.
    Score[i, CensBefore] <- 0
  }
  
  # apply IPCW to all scores
  Err <- Score / IPCWMatrix
  # Average curve demonstrating right-censored brier averaged over test-set, for each time of interest.
  Err <- apply(Err, 1, mean)
  return(Err)
}


#Simulate/prepare data
set.seed(18)

astrain <- simActiveSurveillance(278)
astrain$event <- astrain$event != 0 #take care of competing risks
data.table::setorder(astrain, time, -event)
astrainX <- model.matrix(event ~ . - time, data = astrain)[, -c(1)] #covariates
astrainY <- data.matrix(astrain[, c(9, 8)]) #response

newData <- simActiveSurveillance(208)
newData$event <- newData$event != 0 #take care of competing risks
data.table::setorder(newData, time, -event)
newDataX <- model.matrix(event ~ . - time, data = newData)[, -c(1)] #covariates
newDataY <- data.matrix(newData[, c(9, 8)]) #response

times <- sort(unique(newData$time))


# Fit a cox model
coxfit <- coxph(Surv(time, event != 0) ~ ., data = astrain, x = TRUE)
# Get Cox survivals
coxPred <- summary(survfit(coxfit, newdata = newData), times = times)$surv
# calculate the IPA/brier scores with riskRegression
X2 <- Score(list("PredictionModel" = coxfit), data = newData, formula = Surv(time, event != 0) ~ 1, summary = "ipa", se.fit = 0L, metrics = "brier", contrasts = FALSE, times = times)


#Fit Weibull casebase model
cbModel <- fitSmoothHazard.fit(astrainX, astrainY, family = "glmnet", time = "time", event = "event", alpha = 0,formula_time = ~log(time), ratio = 100, standardize = TRUE)
#we require a specific format for the probabilities
absRiskcb <- absoluteRisk(object = cbModel, newdata = newDataX, time = times, s = "lambda.1se")
absRiskcb[, -c(1)] <- 1 - absRiskcb[, -c(1)] # make it survival probabilities
absRiskcb <- absRiskcb[-c(1), -c(1)] # remove extra time at 0. Also remove the first column (times) as it matches perfectly now


#data is prepared to calculate brier scores
kmCensoredBrierScores <- brierScoreKMCens(predSurvs = absRiskcb, times = times, newData = newData)
weibullCensoredBrierScores<-brierScoreweibCens(predSurvs=absRiskcb, times, newDataX,newDataY)
#> Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

# restructuring to make plotting easier
resultsCB<-rbind(rbind(data.frame(brierScore=kmCensoredBrierScores,times=times,model="KMCensoring")),rbind(data.frame(brierScore=weibullCensoredBrierScores,times=times,model="WeibullCensoring")),data.frame(brierScore = X2$Brier$score$Brier[X2$Brier$score$model == "PredictionModel"],times=times,model="riskRegressionCOX"))




#brierScorePlot
ggplot(data = resultsCB, mapping = aes(x = times, y = brierScore, col = model)) +
  geom_line() +
  geom_rug(sides = "b") +
  xlab("Times") +
  ylab("Brier-Score")+
  ggtitle("KMcensoring vs. WeibullCensoring")

#to visualize, here are the survival probabilities, note red is cox, something weird is happening with the legend.
plot(times, rowMeans(coxPred), type = "l", col = "red", ylab = "survival")
lines(times, rowMeans(absRiskcb))
legend(x = "topright", legend = c("cox", "cbWeibull"), col = c("red", "black"))
title("RED IS COX")

^{Created on 2020-03-12 by the reprex package (v0.3.0)}

Answer 1 · 2020-03-12T18:25:04.000Z

How does IPA compare to likelihood ratio from logistic regression?

Answer 2 · 2020-03-17T14:28:16.000Z

fixing the above issue for brier score only. I no longer use regularization. that is the only change.

If the following brier scores are deemed correct, please close this issue. I'll then open a new one concerning IPA. I plan on using K-M for the censoring distribution

library(casebase)
#> See example usage at http://sahirbhatnagar.com/casebase/
library(survival)
#> Warning: package 'survival' was built under R version 3.6.2
library(riskRegression)
#> riskRegression version 2020.02.05
library(reprex)
library(reshape2)
library(ggplot2)
library(prodlim)
#> Warning: package 'prodlim' was built under R version 3.6.3
library(data.table)
#> Warning: package 'data.table' was built under R version 3.6.2
#> 
#> Attaching package: 'data.table'
#> The following objects are masked from 'package:reshape2':
#> 
#>     dcast, melt



brierScoreKMCens <- function(predSurvs, times, newData) {
  # inverse probability censoring weights
  # Note that we take probabilities right before the drops.
  
  #DIFFERENCES START HERE
  fitCens <- prodlim::prodlim(Hist(time, event != 0) ~ 1, newData, reverse = TRUE)
  IPCW.subject.times <- prodlim::predictSurvIndividual(fitCens, lag = 1) # G(t-|X)
  #DIFFERENCES END HERE
  
  # Empty matrix that will be filled in with the following loop
  Score <- matrix(NA, nrow(predSurvs), ncol(predSurvs))
  IPCWMatrix <- matrix(NA, nrow(predSurvs), ncol(predSurvs))
  # for each point in time we have predicted
  for (i in 1:length(times)) {
    # get number of censored individuals so long as their survival time is less than times[i]
    # these individuals do not have an effect on right-censored brier score.
    CensBefore <- newData$event == 0 & newData$time < times[i]
    # y encompasses all survival times larger than t[i] with a 1, 0 otherwise
    y <- drop(t(newData$time > times[i]))
    # above permits the two parts of right-censored brier score to be calculated, without IPCW, in one line
    Score[i, ] <- (y - predSurvs[i, ])^2
    
    # Generate IPCWMatrix
    IPCWMatrix[i, y == 0] <- IPCW.subject.times[y == 0] # G(t-|X) filled in corresponding positions
    #DIFFERENCES START HERE
    IPCW.time <- predict(fitCens, newdata = newData, times = times[i], level.chaos = 1, mode = "matrix", type = "surv", lag = 1)
    #DIFFERENCES END HERE
    IPCWMatrix[i, y == 1] <- IPCW.time # G(t) filled, same value, for remaining positions.
    # above calculated individuals who don't have an effect on brier score at specfied time. they're scores are set to 0.
    Score[i, CensBefore] <- 0
  }
  
  # apply IPCW to all scores
  Err <- Score / IPCWMatrix
  
  # Average curve demonstrating right-censored brier averaged over test-set, for each time of interest.
  Err <- apply(Err, 1, mean)
  return(Err)
}




brierScoreweibCens <- function(predSurvs= absRiskcb, times=times, newDataX=newDataX,newDataY=newDataY) {
  # inverse probability censoring weights
  
  #DIFFERENCES START HERE
  yCens <- newDataY
  yCens[,1] <- yCens[,1]==0
  fitCens <- fitSmoothHazard(event~log(time) ,data= as.data.frame(yCens), time="time",ratio=100)
  IPCW.subject.times <- absoluteRisk(object = fitCens, time = times, s = "lambda.1se")
  IPCW.subject.times[, -c(1)] <- 1 - IPCW.subject.times[, -c(1)] # make it survival probabilities
  IPCW.subject.times <- rowMeans(IPCW.subject.times[-c(1), -c(1)]) # remove extra time at 0. Also remove the first column (times) as it matches perfectly now
  #DIFFERENCES END HERE
  
  # Empty matrix that will be filled in with the following loop
  Score <- matrix(NA, nrow(predSurvs), ncol(predSurvs))
  IPCWMatrix <- matrix(NA, nrow(predSurvs), ncol(predSurvs))
  # for each point in time we have predicted
  for (i in 1:length(times)) {
    # get number of censored individuals so long as their survival time is less than times[i]
    # these individuals do not have an effect on right-censored brier score.
    CensBefore <- newData$event == 0 & newData$time < times[i]
    # y encompasses all survival times larger than t[i] with a 1, 0 otherwise
    y <- drop(t(newData$time > times[i]))
    # above permits the two parts of right-censored brier score to be calculated, without IPCW, in one line
    Score[i, ] <- (y - predSurvs[i, ])^2
    
    # Generate IPCWMatrix
    IPCWMatrix[i, y == 0] <- IPCW.subject.times[y == 0] # G(t-|X) filled in corresponding positions
    #DIFFERENCES START HERE
    IPCW.time <- IPCW.subject.times[i]
    #DIFFERENCES END HERE
    IPCWMatrix[i, y == 1] <- IPCW.time # G(t) filled, same value, for remaining positions.
    # above calculated individuals who don't have an effect on brier score at specfied time. they're scores are set to 0.
    Score[i, CensBefore] <- 0
  }
  
  # apply IPCW to all scores
  Err <- Score / IPCWMatrix
  # Average curve demonstrating right-censored brier averaged over test-set, for each time of interest.
  Err <- apply(Err, 1, mean)
  return(Err)
}




#Simulate/prepare data
set.seed(18)

astrain <- simActiveSurveillance(278)
astrain$event <- astrain$event != 0 #take care of competing risks
data.table::setorder(astrain, time, -event)
astrainX <- model.matrix(event ~ . - time, data = astrain)[, -c(1)] #covariates
astrainY <- data.matrix(astrain[, c(9, 8)]) #response

newData <- simActiveSurveillance(208)
newData$event <- newData$event != 0 #take care of competing risks
data.table::setorder(newData, time, -event)
newDataX <- model.matrix(event ~ . - time, data = newData)[, -c(1)] #covariates
newDataY <- data.matrix(newData[, c(9, 8)]) #response

times <- sort(unique(newData$time))


# Fit a cox model
coxfit <- coxph(Surv(time, event != 0) ~ ., data = astrain, x = TRUE)
# Get Cox survivals
coxPred <- summary(survfit(coxfit, newdata = newData), times = times)$surv
# calculate the IPA/brier scores with riskRegression
X2 <- Score(list("PredictionModel" = coxfit), data = newData, formula = Surv(time, event != 0) ~ 1, summary = "ipa", se.fit = 0L, metrics = "brier", contrasts = FALSE, times = times)


astrain$event=as.numeric(astrain$event)
cbModel <- fitSmoothHazard(event==1 ~log(time) + ., data = astrain, ratio = 100,time="time",event="event")
#Fit Weibull casebase model (ridge)
#cbModel <- fitSmoothHazard.fit(astrainX, astrainY, family = "glmnet", time = "time", event = "event", alpha = 0,formula_time = ~log(time), ratio = 100, standardize = TRUE)
#we require a specific format for the probabilities
absRiskcb <- absoluteRisk(object = cbModel, newdata = newData, time = times)
absRiskcb[, -c(1)] <- 1 - absRiskcb[, -c(1)] # make it survival probabilities
absRiskcb <- absRiskcb[-c(1), -c(1)] # remove extra time at 0. Also remove the first column (times) as it matches perfectly now


#data is prepared to calculate brier scores
kmCensoredBrierScores <- brierScoreKMCens(predSurvs = absRiskcb, times = times, newData = newData)
weibullCensoredBrierScores<-brierScoreweibCens(predSurvs=absRiskcb, times, newDataX,newDataY)

# restructuring to make plotting easier
resultsCB<-rbind(rbind(data.frame(brierScore=kmCensoredBrierScores,times=times,model="KMCensoring")),rbind(data.frame(brierScore=weibullCensoredBrierScores,times=times,model="WeibullCensoring")),data.frame(brierScore = X2$Brier$score$Brier[X2$Brier$score$model == "PredictionModel"],times=times,model="riskRegressionCOX"))




#brierScorePlot
ggplot(data = resultsCB, mapping = aes(x = times, y = brierScore, col = model)) +
  geom_line() +
  geom_rug(sides = "b") +
  xlab("Times") +
  ylab("Brier-Score")+
  ggtitle("KMcensoring vs. WeibullCensoring")

^{Created on 2020-03-17 by the reprex package (v0.3.0)}

Answer 3 · 2020-03-17T18:29:39.000Z

@Jesse-Islam: This is very hard to follow. I'm glad that you're making reproducible examples, but I think you now have to focus on the minimal part too. I'm not even sure what I'm supposed to look at: the implementation? The theory underlying IPA?

There are too many moving pieces right now. And a lot of redundancies in your code: you have three functions for computing Brier scores, and quickly looking at the code, it seems the code is very similar. Maybe you can try putting it all into one function, focusing on the points you're actually asking feedback on.

And then, I also saw this line:

expModelNull <- fitSmoothHazard(event==1 ~log(time), data = astrain, 
                                ratio = 100,time="time",event="event")

Either you forgot to change the name of your object, or you misunderstood what the exponential hazard is like: the hazard should be constant over time, so there shouldn't be a function of time at all. What you fitted is a Weibull model.

Answer 4 · 2020-03-17T18:36:02.000Z

@Jesse-Islam: This is very hard to follow. I'm glad that you're making reproducible examples, but I think you now have to focus on the minimal part too. I'm not even sure what I'm supposed to look at: the implementation? The theory underlying IPA?

There are too many moving pieces right now. And a lot of redundancies in your code: you have three functions for computing Brier scores, and quickly looking at the code, it seems the code is very similar. Maybe you can try putting it all into one function, focusing on the points you're actually asking feedback on.

And then, I also saw this line:
expModelNull <- fitSmoothHazard(event==1 ~log(time), data = astrain, 
                                ratio = 100,time="time",event="event")
Either you forgot to change the name of your object, or you misunderstood what the exponential hazard is like: the hazard should be constant over time, so there shouldn't be a function of time at all. What you fitted is a Weibull model.

Yes now that I look at it its definitely not minimal at all, part of it was a response to the previous post i made, another part was continuing. I'll restructure this!

Answer 5 · 2020-03-17T19:39:36.000Z

I've updated the above post I made to stay within brierscore only. If deemed correct, I'll move onto IPA in another enhancement issue. I think its close enough now! The fix was removing regularization.

Answer 6 · 2020-03-17T20:05:22.000Z

They do look close, and I'm not sure there is any reason to expect anything closer, unless the data generating mechanism is truly Weibull--did you check?