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.

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
  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)), 

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!


1 Answer 1


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.

  • $\begingroup$ Thank you EdM. I will try the boot package, since even only 100 predictors with 50 bootstrap repetitions takes a long time to run. Your second point seems important - I take it as the sample function should remain within the inner loop? I get different results when I put it in the outside loop (with seed set), not sure where this difference comes from? Thanks for the course notes. He there puts the predictors in one regression fit - that's not meaningful for >1000 predictors? In his example here (hbiostat.org/blog/post/badb) I read it as running 213 models with 600 bootstraps? $\endgroup$ Commented Oct 14, 2022 at 9:49
  • $\begingroup$ @rstats_enthusiast now that I look more closely at the code details, the resampling should be in the outer loop. Presumably your variable selection for the full data set is based on that single data set without resampling. So all of the inner-loop variable selection should be done one bootstrap sample at a time--unless your variable selection process itself involves a bootstrap. That should also speed up the code considerably. You could evaluate 1000 predictors at once with LASSO; see Section 6.2 of SLS. $\endgroup$
    – EdM
    Commented Oct 14, 2022 at 13:43
  • $\begingroup$ @rstats_enthusiast revised answer to incorporate much of what I said in the comment. $\endgroup$
    – EdM
    Commented Oct 14, 2022 at 13:50
  • $\begingroup$ Thanks again. I ran it in the outer loop. A subset of 100 predictors and 200 bootstrap repetitions took around 1h to run. The CIs for ranks were surprisingly narrow - I guess for illustrative purposes I will need to downsize the data to get more unreliable results. $\endgroup$ Commented Oct 15, 2022 at 17:57
  • $\begingroup$ @rstats_enthusiast you don't want to do confidence intervals for this if they are based on means and standard errors of 200 bootstrap samples. You typically only get one data set to play with at a time, not 200 of them, so what's important is the spread in relative ranks among the resamples. The example in the Harrell course notes uses the 0.025 and 0.975 quantiles of the ranks for each predictor, which better describes the risk of relying on ranks from a single data sample. $\endgroup$
    – EdM
    Commented Oct 15, 2022 at 19:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.