R bootstrap boot sampling leads to some predictors having only one value The context is not very important but still, here it is: i am using the boot package to calculate bootstrap confidence interval for different statistics (coefficients of doubly robust method estimates, C-index, etc…) and I happen to face the same problem each time.
If I understood it correctly, the issue is that, no matter which statistic, I need to fit some glm model on each sample and at some point it throws a contrast error because some of the factors predictors in my model have only one level on this sample.
I can think of 3 solutions but either I don’t know how to code it easily or I don’t like it :

*

*I can include some test in my statistic function to exclude a predictor with less than one value but this means that i’m not fitting the same model on each sample, which I don’t like…


*I could just drop the samples where it happens (and I can increase the number of replicates - R parameter of boot function - to compensate the loss of replicates) but I can’t do it easily with boot function since the error breaks the run. I could use some try or tryCatch code inside my statistic function but I feel like there must be an easier way


*I could use the strata parameter of boot function: at the moment, I only stratified the sampling on my dependant variable y="event occured" but I could use a new strata variable that would consist of crossing all categorical predictors so as to have one strata for each combination of factors’ levels. However, this may result in a very large number of stratas...
I hope someone might have some insights on this issue because I can’t find anything online.
Any idea, advice or comment will be much appreciated.
It is not trivial but if needed I'll try to build some reproducible example.

Following @Roland's advice and after some research about how to generalize residuals resampling to GLM models, here is what i did in the end.
# removing na
data_tmp <- df %>% 
  filter(if_all(everything(), ~ !is.na(.)))

# fit the GLM model
mod1 <- glm(y ~ .,
            data = data_tmp %>% select(-ind1),
            family = "binomial")

mod2 <- glm(y ~ .,
            data = data_tmp,
            family = "binomial")


# Bootstrapping C-index
cind_boot_stat <- function(data, indices){
  # generate new y sample with Bernouilli each time
  d <- data %>% mutate(
    y1 = rbinom(nrow(data_tmp), 1, mod1$fitted.values),
y2 = rbinom(nrow(data_tmp), 1, mod2$fitted.values))
  
  # fitting the logistic and calculating C-index
  c1 <- DescTools::Cstat(glm(y1 ~ .,
                             data = d %>% select(-y, -ind1),
                             family = "binomial"))
  
  c2 <- DescTools::Cstat(glm(y2 ~ .,
                             data = d %>% select(-y),
                             family = "binomial"))
  
  return(c(c1, c2))
}

cind_boot <- boot::boot(data_tmp, cind_boot_stat, R = 1000)

Hopefully this is statistically right but at least it seems to work !
