`vimp`

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]))
```

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.2728248 0.04897922 [0.1768273, 0.3688223] TRUE 1.936813e-08
#> s = 6 0.2657869 0.04855174 [0.1706273, 0.3609466] TRUE 3.590847e-08
#> s = 3 0.2584722 0.04868162 [0.1630580, 0.3538865] TRUE 1.111272e-07
#> s = 8 0.2573306 0.04845916 [0.1623524, 0.3523088] TRUE 1.275047e-07
#> s = 7 0.2569966 0.04866336 [0.1616182, 0.3523750] TRUE 1.314361e-07
#> s = 4 0.2569808 0.04866073 [0.1616076, 0.3523541] TRUE 1.305091e-07
#> s = 1 0.2561609 0.04807789 [0.1619299, 0.3503918] TRUE 1.081972e-07
#> s = 5 0.2516121 0.04862427 [0.1563103, 0.3469140] TRUE 2.263829e-07
#> s = 2 0.2416625 0.04832979 [0.1469378, 0.3363872] TRUE 4.563608e-07
```

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:

`Y`

, the outcome`f1`

and`f2`

, the fitted values from a regression procedure (`f1`

used all covariates,`f2`

drops the covariates in`indx`

)`indx`

, which determines the feature I want to estimate variable importance for`run_regression = FALSE`

`folds`

, a list of folds for hypothesis testing and cross-fitting

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.1179355 [0, 0.2311494] FALSE 0.9745693
```

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")
```

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")
```

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")
```

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.

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