Logistic regression feature ranking reliability using bootstrapping [R] I wish to display the (un)reliability in ranking of variables in an omics dataset in the classification of disease vs healthy. In this case it’s mostly of academic interest, as many use rank-based univariable screening to identify the (apparent) “best” markers in e.g. prognostic biomarker studies. I wish to show, in my data (containing more than 1000 predictors measured in 200 samples), that exploratory rank-based feature selection (by e.g. AIC or the Chi square statistic) is unstable in my small sample and that the “best” apparent predictors are unlikely to generalize well.
Say for simplicity that I have 3 variables and a binary outcome variable, and I wish to use 100 bootstrap repetitions to visualize the distributions of ranks (e.g. their confidence intervals) evaluated by the AIC metric (lowest AIC = best rank etc). Keep in mind that I wish to evaluate each predictor by itself, i.e. not multiple regression.
library(dplyr)
library(tidyr)
set.seed(12345)
d <- data.frame(
  y = rbinom(50,1, .5),
  x1 = rnorm(50, 10, 5),
  x2 = rnorm(50, 15, 4),
  x3 = rnorm(50, 20, 1)
)

Now in my code below I am really unsure how to do this correctly, especially whether the sample function should occur within the inner or outer loop, and whether this code in fact represents bootstrapping.
aic_vec <- vector()
aic_df <- list()
B <- 100
variables <- c("x1", "x2", "x3")

for(j in 1:B){
for(i in 1:length(variables)) {
    sample_d <- d[sample(1:nrow(d), nrow(d), replace=TRUE),]
    glm_loop <- purrr::map(variables, ~ glm(as.formula(paste("y ~ ", .x)), data=sample_d, family=binomial))
    aic_vec[i] <- glm_loop[[i]]$aic
}
aic_df[[j]] <- aic_vec   
}

# Store and tidy data
aic_df <- do.call(rbind.data.frame, aic_df)
names(aic_df) <- variables
aic_df <- aic_df %>% mutate(bootstrap_rep = 1:nrow(.)) %>% relocate(bootstrap_rep, .before=variables[1])

# Calculate ranks
aic_df_ranks <- 
aic_df %>% pivot_longer(cols=-bootstrap_rep) %>%
    group_by(bootstrap_rep) %>%
    mutate(rank=dense_rank(desc(-value))) %>%
    select(-value) %>%
    pivot_wider(names_from=name, values_from=rank)

# Display rank frequency per variable
data.frame(
  x1 = table(aic_df_ranks$x1),
  x2 = table(aic_df_ranks$x2),
  x3 = table(aic_df_ranks$x3)) %>%
  select(rank=x1.Var1, x1_freq=x1.Freq, x2_freq=x2.Freq, x3_freq=x3.Freq)

In the next step I wish to repeat this pipeline but also add two fixed covariates on the form
glm_loop <- purrr::map(variables, ~ glm(as.formula(paste("y ~ age + sex + ", .x)), 
    data=sample_d, 
    family=binomial))

and within each bootstrap calculate an LR chi square statistic and use this as a ranking metric. The apparent rank and bootstrapped confidence intervals can then be visually displayed. But first I need to get the code right!

The coding details are off-topic on this site (and I'm not fluent in tidyverse), but you seem to have the critical points in the analysis correct (except perhaps for the placement of the resampling).
First, take a sample of size nrow(data) with replacement from the rows of your data for a bootstrap sample. That starts in the outer loop, once for each bootstrap sample.
Second, repeat all of the steps of your analysis on the bootstrap sample, with whatever feature-selection algorithm you would use on the full data set. This is whatever code you would use to rank features on the full data set, applied to the bootstrap sample. If the feature-selection methods applied to the full original data set don't involve resampling, then you shouldn't resample from your bootstrap sample here. If the feature-selection method you use on the full data sample itself involves resampling, you need to do that on each bootstrap sample as well. Similarly for any other methods like cross-validation to choose a penalty for LASSO. Whatever you would do to the full data set, repeat exactly on each bootstrap sample.
Third, store whatever statistic(s) you want to evaluate. Again, that's whatever you would report for a model on the full data set.
Fourth, repeat the above steps for each of a large number of bootstrap samples.
Fifth, summarize the results over all the bootstrap samples to show the variability in ranking.
The sample code shows only 100 bootstrap samples. You should use more. Frank Harrell uses 1000 in a very similar application in Section 5.4 of his course notes.
With this type of data set (1000 predictors) it would be more efficient to use the boot() function in the R boot package, which allows for parallel processing. You write a function that produces the statistic(s) you want to keep from each bootstrap sample, with the full data set and the row indices selected for that sample accepted as arguments, and submit that function as the statistic argument to boot(). See the help page for the choice of random number generator and timing of seed setting that ensures reproducibility with parallel processing.
