Using precomputed regression function estimates in vimp

Brian D. Williamson

2021-03-01

library("vimp")
library("SuperLearner")

Introduction

In the main vignette, we analyzed the South African heart disease study data (Hastie, Tibshirani, and Friedman 2009) (freely available from the Elements of Statistical Learning website; more information about these data is available here). In each of the analyses, I used run_regression = TRUE. In this vignette, I discuss how to use precomputed regression function estimates with vimp. The results of this analysis replicate the analysis in the main vignette.

# read in the data from the Elements website
library("RCurl")
heart_data <- read.csv(text = getURL("http://web.stanford.edu/~hastie/ElemStatLearn/datasets/SAheart.data"), header = TRUE, stringsAsFactors = FALSE)
# minor data cleaning
heart <- heart_data[, 2:dim(heart_data)[2]]
heart$famhist <- ifelse(heart$famhist == "Present", 1, 0)
# sample-splitting folds for hypothesis testing
heart_folds <- sample(rep(seq_len(2), length = dim(heart)[1]))

Using precomputed regression function estimates without cross-fitting

A first approach: linear regression

As in the main vignette, we first start by fitting only linear regression models. In this section, we use the function vim(); this function does not use cross-fitting to estimate variable importance, and greatly simplifies the code for precomputed regression models.

full_mod <- lm(chd ~ ., data = subset(heart, heart_folds == 1))
full_fit <- predict(full_mod)

# estimate the reduced conditional means for each of the individual variables
# remove the outcome for the predictor matrix
X <- as.matrix(heart[, -dim(heart)[2]])[heart_folds == 2, ] 
# get the full regression fit and use as outcome for reduced regressions; this may provide some stability in this analysis
full_fit_2 <- predict(lm(chd ~ ., data = subset(heart, heart_folds == 2)))
red_mod_sbp <- lm(full_fit ~ X[,-1])
red_fit_sbp <- predict(red_mod_sbp)
red_mod_tob <- lm(full_fit ~ X[,-2])
red_fit_tob <- predict(red_mod_tob)
red_mod_ldl <- lm(full_fit ~ X[,-3])
red_fit_ldl <- predict(red_mod_ldl)
red_mod_adi <- lm(full_fit ~ X[,-4])
red_fit_adi <- predict(red_mod_adi)
red_mod_fam <- lm(full_fit ~ X[,-5])
red_fit_fam <- predict(red_mod_fam)
red_mod_tpa <- lm(full_fit ~ X[, -6])
red_fit_tpa <- predict(red_mod_tpa)
red_mod_obe <- lm(full_fit ~ X[,-7])
red_fit_obe <- predict(red_mod_obe)
red_mod_alc <- lm(full_fit ~ X[, -8])
red_fit_alc <- predict(red_mod_alc)
red_mod_age <- lm(full_fit ~ X[,-9])
red_fit_age <- predict(red_mod_age)

# plug the regression function estimates into vim
lm_vim_sbp <- vim(Y = heart$chd, f1 = full_fit, f2 = red_fit_sbp, indx = 1, run_regression = FALSE, type = "r_squared", folds = heart_folds)
lm_vim_tob <- vim(Y = heart$chd, f1 = full_fit, f2 = red_fit_tob, indx = 2, run_regression = FALSE, type = "r_squared", folds = heart_folds)
lm_vim_ldl <- vim(Y = heart$chd, f1 = full_fit, f2 = red_fit_ldl, indx = 3, run_regression = FALSE, type = "r_squared", folds = heart_folds)
lm_vim_adi <- vim(Y = heart$chd, f1 = full_fit, f2 = red_fit_adi, indx = 4, run_regression = FALSE, type = "r_squared", folds = heart_folds)
lm_vim_fam <- vim(Y = heart$chd, f1 = full_fit, f2 = red_fit_fam, indx = 5, run_regression = FALSE, type = "r_squared", folds = heart_folds)
lm_vim_tpa <- vim(Y = heart$chd, f1 = full_fit, f2 = red_fit_tpa, indx = 6, run_regression = FALSE, type = "r_squared", folds = heart_folds)
lm_vim_obe <- vim(Y = heart$chd, f1 = full_fit, f2 = red_fit_obe, indx = 7, run_regression = FALSE, type = "r_squared", folds = heart_folds)
lm_vim_alc <- vim(Y = heart$chd, f1 = full_fit, f2 = red_fit_alc, indx = 8, run_regression = FALSE, type = "r_squared", folds = heart_folds)
lm_vim_age <- vim(Y = heart$chd, f1 = full_fit, f2 = red_fit_age, indx = 9, run_regression = FALSE, type = "r_squared", folds = heart_folds)

# make a table with the estimates using the merge_vim() function
library("dplyr")
lm_mat <- merge_vim(lm_vim_sbp, lm_vim_tob, lm_vim_ldl, lm_vim_adi,
                lm_vim_fam, lm_vim_tpa, lm_vim_obe, lm_vim_alc, lm_vim_age)
# print out the matrix
lm_mat
#> Variable importance estimates:
#>       Estimate  SE         95% CI                 VIMP > 0 p-value     
#> s = 9 0.2728031 0.04881682 [0.1771239, 0.3684823] TRUE     1.731885e-08
#> s = 6 0.2657658 0.04838169 [0.1709394, 0.3605921] TRUE     3.221815e-08
#> s = 3 0.2584517 0.04851628 [0.1633615, 0.3535418] TRUE     1.008172e-07
#> s = 8 0.2573102 0.04829716 [0.1626495, 0.3519708] TRUE     1.157730e-07
#> s = 7 0.2569762 0.04849923 [0.1619194, 0.3520329] TRUE     1.193987e-07
#> s = 4 0.2569604 0.04849745 [0.1619071, 0.3520137] TRUE     1.185522e-07
#> s = 1 0.2561405 0.04791674 [0.1622254, 0.3500556] TRUE     9.795631e-08
#> s = 5 0.2515921 0.04845900 [0.1566142, 0.3465700] TRUE     2.063392e-07
#> s = 2 0.2416433 0.04815970 [0.1472520, 0.3360346] TRUE     4.168884e-07

Estimating the importance of a single variable with a large library of learners

Say I have already computed one regression function estimate for the full set of features; call the fitted values from this procedure f1. Further, suppose that I have done estimated a second regression function by removing the indx of interest; call the fitted values from this procedure f2. Then I can use vimp to compute variable importance based on these fitted values. There five main arguments to vim() when I have already computed:

I can estimate the variable importance for family history only using one base learner (for illustration; also, using a small number of cross-fitting folds and a small number of Super Learner cross-validation folds) as follows:

# small learners library
learners.2 <- c("SL.ranger")

# small number of folds
V <- 2
sl_cvcontrol <- list(V = 2)

# now estimate variable importance
fam_vim <- vim(Y = heart$chd, X = heart[, -dim(heart)[2]], indx = 5, SL.library = learners.2, na.rm = TRUE, env = environment(), cvControl = sl_cvcontrol)

I said earlier that I want to obtain estimates of all individual features in these data, so let’s choose type A behavior (tpa) next. Now that I have estimated variable importance for family history, the full_fit object contains our estimate of \(\mu_{P_0}\). Since I have spent the time to estimate this using SuperLearner(), there is no reason to estimate this function again. This leads me to choose (2) above for estimating the importance of type A behavior, since I have already estimated variable importance on one feature in this dataset. Using the small learners library (again only for illustration) yields

# specify that full_fit doesn't change; but note that this was fit on the first half of the data
full_fit <- fam_vim$full_fit
# get the full fit on the second half of the data; for stability later
full_fit_2 <- predict(SuperLearner::SuperLearner(Y = heart$chd[fam_vim$folds == 2], X = heart[fam_vim$folds == 2, -dim(heart)[2]], SL.library = learners.2, cvControl = sl_cvcontrol))$pred

# estimate variable importance for the average number of rooms
reduced_fit <- SuperLearner::SuperLearner(Y = full_fit_2, X = heart[fam_vim$folds == 2, -c(6, dim(heart)[2]), drop = FALSE], SL.library = learners.2, cvControl = sl_cvcontrol)
red_fit <- predict(reduced_fit)$pred
tpa_vim <- vim(Y = heart$chd, f1 = full_fit_2, f2 = red_fit, indx = 6, run_regression = FALSE, type = "r_squared", folds = heart_folds)
#> Warning in vim(Y = heart$chd, f1 = full_fit_2, f2 = red_fit, indx = 6,
#> run_regression = FALSE, : Original estimate < 0; returning zero.

tpa_vim
#> Variable importance estimates:
#>       Estimate SE        95% CI         VIMP > 0 p-value 
#> s = 6 0        0.1177334 [0, 0.2307532] FALSE    0.974724

This takes approximately 5 seconds — now rather than estimating both conditional means, I am only estimating one.

I can achieve the same results by also estimating the regression function based on all variables myself, and plugging these estimates into vim(). Note that I am using the folds from the first fit above so that the estimates are consistent with each other; if I was doing this analysis from scratch, I could use heart_folds as I did in the linear regression example. Then vim() returns variable importance estimates based on the inputted fitted values. For example, let’s estimate variable importance for current alcohol consumption using this approach.

# set up the data, removing the columns for alcohol use and chd
x <- heart[, -c(8, dim(heart)[2])]

# fit an RF model using SuperLearner
reduced_mod <- SuperLearner(Y = full_fit_2, X = x[fam_vim$folds == 2, ], SL.library = learners.2, cvControl = sl_cvcontrol)
reduced_fit <- predict(reduced_mod)$pred
# this takes 2 seconds

# estimate variable importance
alc_vim <- vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_fit, indx = 8, run_regression = FALSE, type = "r_squared", folds = fam_vim$folds)
#> Warning in vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_fit, indx = 8, :
#> Original estimate < 0; returning zero.

I can obtain estimates for the remaining individual features in the same way (again using a small library for illustration):

reduced_sbp <- predict(SuperLearner(Y = full_fit_2, X = heart[fam_vim$folds == 2, -c(1, dim(heart)[2])], SL.library = learners.2, cvControl = sl_cvcontrol))$pred
sbp_vim <- vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_sbp,
                indx = 1, run_regression = FALSE, type = "r_squared", folds = fam_vim$folds)
#> Warning in vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_sbp, indx = 1, :
#> Original estimate < 0; returning zero.

reduced_tob <- predict(SuperLearner(Y = full_fit_2, X = heart[fam_vim$folds == 2, -c(2, dim(heart)[2])], SL.library = learners.2, cvControl = sl_cvcontrol))$pred
tob_vim <- vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_tob,
                indx = 2, run_regression = FALSE, type = "r_squared", folds = fam_vim$folds)
#> Warning in vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_tob, indx = 2, :
#> Original estimate < 0; returning zero.

reduced_ldl <- predict(SuperLearner(Y = full_fit_2, X = heart[fam_vim$folds == 2, -c(3, dim(heart)[2])], SL.library = learners.2, cvControl = sl_cvcontrol))$pred
ldl_vim <- vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_ldl,
                indx = 3, run_regression = FALSE, type = "r_squared", folds = fam_vim$folds)
#> Warning in vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_ldl, indx = 3, :
#> Original estimate < 0; returning zero.

reduced_adi <- predict(SuperLearner(Y = full_fit_2, X = heart[fam_vim$folds == 2, -c(4, dim(heart)[2])], SL.library = learners.2, cvControl = sl_cvcontrol))$pred
adi_vim <- vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_adi,
                indx = 4, run_regression = FALSE, type = "r_squared", folds = fam_vim$folds)
#> Warning in vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_adi, indx = 4, :
#> Original estimate < 0; returning zero.

reduced_obe <- predict(SuperLearner(Y = full_fit_2, X = heart[fam_vim$folds == 2, -c(7, dim(heart)[2])], SL.library = learners.2, cvControl = sl_cvcontrol))$pred
obe_vim <- vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_obe,
                indx = 7, run_regression = FALSE, type = "r_squared", folds = fam_vim$folds)
#> Warning in vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_obe, indx = 7, :
#> Original estimate < 0; returning zero.

reduced_age <- predict(SuperLearner(Y = full_fit_2, X = heart[fam_vim$folds == 2, -c(9, dim(heart)[2])], SL.library = learners.2, cvControl = sl_cvcontrol))$pred
age_vim <- vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_age,
                indx = 9, run_regression = FALSE, type = "r_squared", folds = fam_vim$folds)
#> Warning in vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_age, indx = 9, :
#> Original estimate < 0; returning zero.

Now that I have estimates of each of individual feature’s variable importance, I can view them all simultaneously by plotting:

library("dplyr")
library("tibble")
library("ggplot2")
library("cowplot")
theme_set(theme_cowplot())
# combine the objects together
ests <- merge_vim(sbp_vim, tob_vim, ldl_vim, adi_vim,
                fam_vim, tpa_vim, obe_vim, alc_vim, age_vim)
all_vars <- c("Sys. blood press.", "Tobacco consump.", "LDL cholesterol",
              "Adiposity", "Family history", "Type A behavior", "Obesity",
              "Alcohol consump.", "Age")

est_plot_tib <- ests$mat %>% 
  mutate(
    var_fct = rev(factor(s, levels = ests$mat$s,
                     labels = all_vars[as.numeric(ests$mat$s)], 
                     ordered = TRUE))
  )

# plot
est_plot_tib %>%
  ggplot(aes(x = est, y = var_fct)) +
  geom_point() +
  geom_errorbarh(aes(xmin = cil, xmax = ciu)) +
  xlab(expression(paste("Variable importance estimates: ", R^2, sep = ""))) +
  ylab("") +
  ggtitle("Estimated individual feature importance") +
  labs(subtitle = "in the South African heart disease study data")

Estimating variable importance for a group of variables

Now that I have estimated variable importance for each of the individual features, I can estimate variable importance for each of the groups that I mentioned above: biological and behavioral features.

The only difference between estimating variable importance for a group of features rather than an individual feature is that now I specify a vector for s; I can use any of the options listed in the previous section to compute these estimates.

# get the estimates
reduced_behav <- predict(SuperLearner(Y = heart$chd[fam_vim$folds == 2], X = heart[fam_vim$folds == 2, -c(2, 6, 8, dim(heart)[2])], SL.library = learners.2, cvControl = sl_cvcontrol))$pred
behav_vim <- vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_behav, indx = c(2, 6, 8), run_regression = FALSE, type = "r_squared", folds = fam_vim$folds)
#> Warning in vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_behav, indx = c(2, :
#> Original estimate < 0; returning zero.

reduced_bios <- predict(SuperLearner(Y = heart$chd[fam_vim$folds == 2], X = heart[fam_vim$folds == 2, -c(1, 3, 4, 5, 7, 9, dim(heart)[2])], SL.library = learners.2, cvControl = sl_cvcontrol))$pred
bios_vim <- vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_bios, indx = c(1, 3, 4, 5, 7, 9), run_regression = FALSE, type = "r_squared", folds = fam_vim$folds)
#> Warning in vim(Y = heart$chd, f1 = full_fit_2, f2 = reduced_bios, indx = c(1, :
#> Original estimate < 0; returning zero.

# combine and plot
groups <- merge_vim(behav_vim, bios_vim)
all_grp_nms <- c("Behavioral features", "Biological features")

grp_plot_tib <- groups$mat %>% 
  mutate(
    grp_fct = factor(case_when(
      s == "2,6,8" ~ "1",
      s == "1,3,4,5,7,9" ~ "2"
    ), levels = c("1", "2"),  labels = all_grp_nms, ordered = TRUE)
  )
grp_plot_tib %>%
  ggplot(aes(x = est, y = grp_fct)) +
  geom_point() +
  geom_errorbarh(aes(xmin = cil, xmax = ciu)) +
  xlab(expression(paste("Variable importance estimates: ", R^2, sep = ""))) +
  ylab("") +
  ggtitle("Estimated feature group importance") +
  labs(subtitle = "in the South African heart disease study data")

Using precomputed regression function estimates with cross-fitting

In this section, we will use cross-fitting and pre-computed estimates of the regression functions. This can be especially useful if you have already run a call to CV.SuperLearner – that function returns estimates based on each observation being part of the hold-out set. However, while this approach can save you some computation time, it requires a hefty amount of mental overhead.

To replicate the Super Learning analysis from above using cross-fitting, we have to first create a set of folds for hypothesis testing (this involves simply splitting the dataset in two halves) and then create two sets of folds (one for each half) specifying the indices for cross-fitting. For this analysis, we will use V = 5 folds for cross-fitting (as we did in the main vignette).

# the outer folds for hypothesis testing
outer_folds <- sample(rep(seq_len(2), length = length(heart$chd)))
# the first set of inner folds for cross-fitting
V <- 2
inner_folds_1 <- sample(rep(seq_len(V), length = sum(outer_folds == 1)))
inner_folds_2 <- sample(rep(seq_len(V), length = sum(outer_folds == 2)))

First, we estimate the regression function based on all variables:

# set up some objects to simplify the code
y_1 <- heart$chd[outer_folds == 1]
x_1 <- heart[outer_folds == 1, -ncol(heart), drop = FALSE]
y_2 <- heart$chd[outer_folds == 2]
x_2 <- heart[outer_folds == 2, -ncol(heart), drop = FALSE]
# set up a list to hold the fitted values
fhat_ful <- list()
for (v in 1:V) {
  # fit the Super Learner
  fit <- SuperLearner::SuperLearner(Y = y_1[inner_folds_1 != v],
                                    X = x_1[inner_folds_1 != v, , drop = FALSE], 
                                    SL.library = learners.2, cvControl = list(V = 2))
  fhat_ful[[v]] <- SuperLearner::predict.SuperLearner(fit, newdata = x_1[inner_folds_1 == v, , drop = FALSE])$pred
}

Next, to estimate the importance of each variable, we need to estimate the reduced regression function. To reduce the number of loops, we will estimate all reduced regression functions at once as follows:

# set up lists to hold the fitted values
fhat_red_sbp <- fhat_red_tob <- fhat_red_ldl <- fhat_red_adi <- fhat_red_fam <- fhat_red_tpa <- fhat_red_obe <- fhat_red_alc <- fhat_red_age <- list()
vars <- c("sbp", "tob", "ldl", "adi", "fam", "tpa", "obe", "alc", "age")
for (i in 1:length(vars)) {
  for (v in 1:V) {
    # fit the Super Learner
    fit <- SuperLearner::SuperLearner(Y = y_2[inner_folds_2 != v],
                                      X = x_2[inner_folds_2 != v, -i, drop = FALSE], SL.library = learners.2, cvControl = list(V = 2))
    # note the use of "eval" and "parse" to assign to the correct list
    eval(parse(text = paste0("fhat_red_", vars[i], "[[v]] <- SuperLearner::predict.SuperLearner(fit, newdata = x_2[inner_folds_2 == v, -i, drop = FALSE])$pred")))
  }  
}

Then we can plug these values into vimp_rsquared() (or equivalently, cv_vim() with type = "r_squared") as follows:

# note that folds is a list! also, V must equal the length of folds
cf_sbp_vim <- vimp_rsquared(Y = heart$chd, f1 = fhat_ful, f2 = fhat_red_sbp, indx = 1, run_regression = FALSE, V = V, folds = list(outer_folds, list(inner_folds_1, inner_folds_2)))
cf_tob_vim <- vimp_rsquared(Y = heart$chd, f1 = fhat_ful, f2 = fhat_red_tob, indx = 2, run_regression = FALSE, V = V, folds = list(outer_folds, list(inner_folds_1, inner_folds_2)))
cf_ldl_vim <- vimp_rsquared(Y = heart$chd, f1 = fhat_ful, f2 = fhat_red_ldl, indx = 3, run_regression = FALSE, V = V, folds = list(outer_folds, list(inner_folds_1, inner_folds_2)))
cf_adi_vim <- vimp_rsquared(Y = heart$chd, f1 = fhat_ful, f2 = fhat_red_adi, indx = 4, run_regression = FALSE, V = V, folds = list(outer_folds, list(inner_folds_1, inner_folds_2)))
cf_fam_vim <- vimp_rsquared(Y = heart$chd, f1 = fhat_ful, f2 = fhat_red_fam, indx = 5, run_regression = FALSE, V = V, folds = list(outer_folds, list(inner_folds_1, inner_folds_2)))
cf_tpa_vim <- vimp_rsquared(Y = heart$chd, f1 = fhat_ful, f2 = fhat_red_tpa, indx = 6, run_regression = FALSE, V = V, folds = list(outer_folds, list(inner_folds_1, inner_folds_2)))
cf_obe_vim <- vimp_rsquared(Y = heart$chd, f1 = fhat_ful, f2 = fhat_red_obe, indx = 7, run_regression = FALSE, V = V, folds = list(outer_folds, list(inner_folds_1, inner_folds_2)))
cf_alc_vim <- vimp_rsquared(Y = heart$chd, f1 = fhat_ful, f2 = fhat_red_alc, indx = 8, run_regression = FALSE, V = V, folds = list(outer_folds, list(inner_folds_1, inner_folds_2)))
cf_age_vim <- vimp_rsquared(Y = heart$chd, f1 = fhat_ful, f2 = fhat_red_age, indx = 9, run_regression = FALSE, V = V, folds = list(outer_folds, list(inner_folds_1, inner_folds_2)))

cf_ests <- merge_vim(cf_sbp_vim, cf_tob_vim, cf_ldl_vim, cf_adi_vim, cf_fam_vim, cf_tpa_vim, cf_obe_vim, cf_alc_vim, cf_age_vim)

And we can view them all simultaneously by plotting:

cf_est_plot_tib <- cf_ests$mat %>% 
  mutate(
    var_fct = rev(factor(s, levels = cf_ests$mat$s,
                     labels = all_vars[as.numeric(cf_ests$mat$s)], 
                     ordered = TRUE))
  )

# plot
cf_est_plot_tib %>%
  ggplot(aes(x = est, y = var_fct)) +
  geom_point() +
  geom_errorbarh(aes(xmin = cil, xmax = ciu)) +
  xlab(expression(paste("Variable importance estimates: ", R^2, sep = ""))) +
  ylab("") +
  ggtitle("Cross-fitted individual feature importance estimates") +
  labs(subtitle = "in the South African heart disease study data")

Finally, we can estimate and plot group importance:

fhat_red_bios <- fhat_red_behav <- list()
for (v in 1:V) {
    # fit the Super Learner
    fit_behav <- SuperLearner::SuperLearner(Y = y_2[inner_folds_2 != v],
                                      X = x_2[inner_folds_2 != v, -c(2, 6, 8), drop = FALSE], SL.library = learners.2, cvControl = list(V = 2))
    fit_bios <- SuperLearner::SuperLearner(Y = y_2[inner_folds_2 != v],
                                      X = x_2[inner_folds_2 != v, -c(1, 3, 4, 5, 7, 9), drop = FALSE], SL.library = learners.2, cvControl = list(V = 2))
    fhat_red_behav[[v]] <- SuperLearner::predict.SuperLearner(fit_behav, newdata = x_2[inner_folds_2 == v, -c(2, 6, 8), drop = FALSE])$pred
    fhat_red_bios[[v]] <- SuperLearner::predict.SuperLearner(fit_bios, newdata = x_2[inner_folds_2 == v, -c(1, 3, 4, 5, 7, 9), drop = FALSE])$pred
}
cf_behav_vim <- vimp_rsquared(Y = heart$chd, f1 = fhat_ful, f2 = fhat_red_behav, indx = c(2, 6, 8), run_regression = FALSE, V = V, folds = list(outer_folds, list(inner_folds_1, inner_folds_2)))
cf_bios_vim <- vimp_rsquared(Y = heart$chd, f1 = fhat_ful, f2 = fhat_red_bios, indx = c(1, 3, 4, 5, 7, 9), run_regression = FALSE, V = V, folds = list(outer_folds, list(inner_folds_1, inner_folds_2)))
cf_groups <- merge_vim(cf_behav_vim, cf_bios_vim)

cf_grp_plot_tib <- cf_groups$mat %>% 
  mutate(
    grp_fct = factor(case_when(
      s == "2,6,8" ~ "1",
      s == "1,3,4,5,7,9" ~ "2"
    ), levels = c("1", "2"),  labels = all_grp_nms, ordered = TRUE)
  )
cf_grp_plot_tib %>%
  ggplot(aes(x = est, y = grp_fct)) +
  geom_point() +
  geom_errorbarh(aes(xmin = cil, xmax = ciu)) +
  xlab(expression(paste("Variable importance estimates: ", R^2, sep = ""))) +
  ylab("") +
  ggtitle("Cross-fitted feature group importance estimates") +
  labs(subtitle = "in the South African heart disease study data")

Conclusion

In this document, we learned a second method for computing variable importance estimates: rather than having vimp run all regression functions for you, you can compute your own regressions and pass these to vimp. The results are equivalent, but there is a tradeoff: what you save in computation time by only computing the full regression once must be balanced with the mental overhead of correctly computing the regressions. Additionally, this task is more difficult when using cross-fitted variable importance, which I recommend in nearly all cases when using flexible machine learning tools.

References

Hastie, T, R Tibshirani, and J Friedman. 2009. The Elements of Statistical Learning.