Getting confidence interval for test data from model fitted on bootstrapped train data

Question

I am using Boostrappping to fit a cox regression model in R. I am using boot function from R. I have divided my data to test and train. I generate bootstraps of train data and compute concordance index(statistics) on those bootstraps.My question is how I can apply those 100 bootstrapped models to test data and get Confidence intervals for test data as well. Following is my example

library(splitTools)
library(survival)
library(dplyr)
library(intsurv)
library(boot)

bootstrap_train <- function(data, indices)
{
set.seed(3451)
train<- data[indices,]
res.cox<- coxph(Surv(time,status)~., data=train) # fitting cox model
train_predict <- predict(res.cox,newdata = train,type = 'risk') #getting train prediction
cindex_train<- cIndex(train$time, event = train$status, train_predict, weight = NULL) #getting concordance index
return(cindex_train[1])}

cox_train<- X_train   #train data frame with variables and endpoint
cox_test<- X_test # test data frame with variables and endpoint 
train_boots <- boot(data=cox_train, statistic=bootstrap_train, R=100)

This returns me computed statistics and I can get confidence intervals with boot.ci function. However I am not sure how to fit all those 100 train models to test data and get 100 prediction for test data as well.I also tried passing my test data to boot function assuming that each time fitted model in statistic function will be applied to test data and I can get predictions on new data and store those predictions (in my case concordance index) in a list and then I use that list to calculate 95% confidence interval manually using the normal distribution. However, this gives me an error. Here is what I tried

  bootstrap_train <- function(data, indices, test)
    {
    set.seed(3451)
    list_test_cIndex<- NULL # list to store test results 
    train<- data[indices,]
    res.cox<- coxph(Surv(time,status)~., data=train) # fitting model on train data 
    train_predict <- predict(res.cox,newdata = train,type = 'risk') 
    cindex_train<- cIndex(train$time, event = train$status, train_predict, weight = NULL) 
    test_predict <- predict(res.cox,newdata = test,type = 'risk') # fitting trained model on test 
    cindex_test <- cIndex(test$time, event = test$status, test_predict, weight = NULL)
    list_test_cIndex<- cbind(list_test_CI, cindex_test[1]) # storing test cIndex in list 
    return(list(cindex_train[1]),list_test_cIndex )}

cox_train<- X_train   #train data frame with variables and endpoint
cox_test<- X_test # test data frame with variables and endpoint 
c(train_boots, cIndex_test) %<-% boot(data=cox_train, statistic=bootstrap_train, R=100, test= 
                                           cox_test)

 #calculating Confidence intervals of concordance index for test data 
 error_test <- qnorm(0.975)*sd(cIndex_test) 
 c(CI_low_test, CI_hi_test) %<-% list(median(cIndex_test)-error_test,median(cIndex_test)+error_test)

Error in t.star[r, ] <- res[[r]] : 
     incorrect number of subscripts on matrix

Note that I do not want to bootstrap my test data.

Answer 1

I'll provide suggestions to help solve your error later, but first please consider whether you should alter your strategy a bit.

I have divided my data to test and train.

Frank Harrell (who popularized the concordance index that you use to evaluate model performance) recommends against separate train/test data splits unless you have more than 20,000 or so cases. Otherwise you run the risk of imprecision both in the model and in the testing.

Unless you have that many cases, much of what you want to accomplish would better be done by using all your data to build the model. You then evaluate the likely performance on new data similarly to how you are suggesting: modeling on each of multiple bootstrapped data sets but evaluating performance in each case on the entire original data set. The validate() function in Harrell's rms package does that evaluation on a large number of measures of model performance, including the $D_{xy}$ that is easily converted to concordance.

how I can apply those 100 bootstrapped models to test data and get Confidence intervals

If you only run 100 bootstraps, your estimates of 95% confidence intervals (CI) might be imprecise: the 2.5% of cases at each end of the distribution will only be 2 or 3 cases each. You are depending on an assumption of normality in the distribution of estimated concordance values to get the CI. For empirical 95% confidence intervals, a standard approach is to run 999 bootstraps, sort the values, and take the 25th and 975th as the 95% CI. That makes the estimated CI less susceptible to the vagaries of just a couple of bootstrap samples and doesn't require a normality assumption. Even if you just submit your set of values to boot.ci(), providing more values should give you better estimates.

prediction (in my case concordance index)

Concordance is a fine way to evaluate the ability of a model to discriminate among cases, but it has its limits. It doesn't provide an estimate of calibration--how well predicted and observed survival probabilities agree. The calibrate() function in rms evaluates that aspect of the model via resampling. Also, concordance is not very sensitive for comparing different models.

In terms of the error you are getting, that seems to be a coding error and thus is technically off topic on this site. Your code doesn't explicitly include either a t.star or a res object, so the error is coming from some function that your code is calling that isn't getting properly formatted data.* That said, I think you could simplify your model a lot (and thus minimize the chance of such errors) by relying on standard survival package functions. Its concordance() function can return the concordance from a model directly or, with a newdata argument, evaluate the concordance of a given model on a new data set.

---

*It's easy to get confused between lists and vectors, data frames and matrices, etc. I notice that one line of your code

list_test_cIndex<- cbind(list_test_CI, cindex_test[1]) # storing test cIndex in list 

says that you are trying to store the cIndex in a list. I believe that cbind() can work with lists, but I don't see where list_test_CI is ever defined (particularly as a list), and as cindex_test[1] is a single numeric value what you are generating might be a vector rather than a list. That's an important distinction. Along that line of thought, the structure of what you intend to return in your second bootstrap_train() function isn't clear; do you want that to be a list or a vector? Be careful where you place your parentheses.

Answer 2

Bootstrapping in a loop is a much simpler way of applying bootstrapped train models on original test data. Nevertheless, the code could be simplified further.

library(splitTools)
library(survival)
library(dplyr)
library(intsurv)

model_fit <- function(cox_train, cox_test){
      
      
      res.cox<- coxph(Surv(time,status)~., data=cox_train)
      train_predict <- predict(res.cox,newdata = cox_train,type = 'risk') # fit on train
      cindex_train<- cIndex(cox_train$time, event = cox_train$status, train_predict, weight = NULL)
      cIndex_train <- cindex_train[1]

      test_predict <- predict(res.cox,newdata = cox_test,type = 'risk') # fit on test
      cindex_test<- cIndex(cox_test$time, event = cox_test$status,test_predict, weight = NULL)
      cIndex_test <- cindex_test[1]
      
      return(list(cIndex_train, cIndex_test))
      
    }
    
    
    
    set.seed(3451)
    cIndex_train<- NULL # empty list to store k train statistics, in this case concordance index 
    cIndex_train<- NULL  # empty list to store k test statistics
    
    for (i in 1:k) # k should be number of desired bootstraps 
    {
      
      boot_strap <- sample(rownames(cox_train), replace = TRUE) # apply resampling to train data 
      cox_train %<-% cox_train[boot_strap,]  #get resampled observation from train data 
    
      # fit model on train and evaluate on original test data
      c(cindex_train, cindex_test)%<-% model_fit(cox_train, cox_test) 
      cIndex_train<- rbind(cIndex_train, cindex_train)  # append your k train statistics to list 
      
      cIndex_test<- rbind(cIndex_test, cindex_test) # append your k test statistics to list 
      
      
    }
    
    # Apply your desired method of computing CI's e.g. Normal bootstrap confidence limits
    error_train <- qnorm(0.975)*sd(cIndex_train) 
    c(Conf_low_train, Conf_hi_train) %<-% list(median(cIndex_train)-error_train, median(cIndex_train)+error_train)
    
    
    error_valid <- qnorm(0.975)*sd(cIndex_test)
    c(Conf_low_valid, Conf_hi_valid) %<-% list(median(cIndex_test)-error_valid, median(cIndex_test)+error_valid)