# Relative variable importance/explained variation from a single model fit

I am seeking a measure of relative variable importance or relative explained variation that will apply to all types of linear and nonlinear regression models and that requires only fitting one model. As an example suppose we have the model y ~ x1 + x2 + x1:x2 with the last term representing an interaction term. The model may be of any type including logistic, ordinal, and Cox proportional hazards regression models. Various importance measures have been developed based on Wald $$\chi^2$$ statistics, and one can compute a measure of importance of x1 by comparing the 2 degree of freedom Wald statistic for x1 and x1:x2 with the 3 d.f. Wald $$\chi^2$$ for the complete model. This is all done using the full model fit. This will not transport to Bayesian models though.

One simple thought is to compute the linear correlation coefficient between lp1 and lp where lp1 is the portion of the linear predictor that involves x1 (which will have two terms combined into one sum, using two estimated regression coefficients) and lp is the overall model fit's linear predictor $$X\hat{\beta}$$. I know that traditional $$R^2$$-based measures will use two $$R^2$$: from the large and reduced models and this will handle collinearities appropriately, but I'm seeking measures that work reasonable well from a single model fit.

Does anyone have an idea of how well the corr(lp1, lp) approach will work or have a better idea for estimating relative variable importance from a single model fit? Measures that operate on the estimated full-model regression coefficients will work in both frequentist and Bayesian settings.

Update: Another thing to try. For a relative explained variation in the response that is due to x1, let $$\hat{y}$$ be the full model linear predictor $$X\hat{\beta}$$. Fit a linear model to the linear predictor: $$\hat{y} \sim x1 + x2 + x1:x2$$. The $$R^2$$ from this fit will be 1.0. Quantify the relative importance of x1 by the partial $$R^2$$ that is due to the 2 degrees of freedom involving x1, i.e., omitting the x2 main effect from the contrast. In a linear model one can compute any needed partial $$R^2$$ from contrasts without refitting anything. Are there downsides or strange collinearity issues?

Update 2: This second method yields an exactly correct result when the initial model is a linear model, as the relative $$R^2$$ values from predicting the linear predictor, when multiplied by the $$R^2$$ in the original fit, give exactly the partial $$R^2$$ in the original fit. So I'm pretty happy with method 2. For nonlinear models it is using the linear model as a bridge to get relative explained variation.

Update 3: @Michael M discussed global surrogate models below. This would extend the update in Update 2 along the following lines:

1. Fit a predictive model using any desired method, including random forests (if the sample size is incredibly large)
2. Predict the predictions from that fit using a fairly saturated linear model, i.e., one that uses regression splines for all continuous predictors and includes all two-way interactions (tensor splines). If this fit has $$R^{2} > 0.9$$ we would call it an adequate approximation to the random forest.
3. Use that fitted linear model to decompose explained outcome variation (relative $$R^2$$) in the random forest in an interpretable way. In the R rms package this decomposition combines any effect with higher-order effects involving that variable, so the partial $$R^2$$ for a predictor would count its nonlinear and interactive effects.

Update 4: I've implemented the new general purpose relative explained variation method. Explanation and examples are here. It would be nice to compare its performance with random permutation-based methods.

• Two model-agnostic measures: Permutation importance; SHAP importance (via permutation SHAP or some other agnostic-approximation). They both have their ups and downs, but they require only one model fit. Commented Sep 13, 2023 at 12:58
• And here a quite cool paper of Mr. Kruskal: jstor.org/stable/2684310, where he rediscovered Shapley values as fairest variable importance measure in models. But his approach requires refitting the model a zillion times. Commented Sep 13, 2023 at 13:30
• Why aren't Shapley/SHAP values what you want? Commented Sep 13, 2023 at 15:43
• Not at all @FrankHarrell. I might not have explained the approaches very well. But there is no refitting involved, even if you might have other such methods in mind. For SHAP, one uses permutation SHAP or KernelSHAP. Again, no refit is required. (I have implemented all mentioned methods in R.) Commented Sep 13, 2023 at 15:46
• To me, the estimators discussed in her papers are fungible. In other words, the metrics can be adapted to other methods. Commented Sep 13, 2023 at 16:45

I don't think there is the perfect measure to do so. Still, here are two options from explainable ML that do not need to refit:

## Permutation feature importance

The average loss (or some other meaningful performance measure) is calculated on the reference data, typically the validation data to get rid of model optimism. Then, for each feature, its data column is randomly shuffled and the increase in average loss is calculated. If the increase is small, it means the feature is unimportant. The measure was introduced by Leo Breiman for random forests and then later generalized.

### The algorithm

Screenshot from our ETH Zurich lecture

($$\hat m$$ is the fitted model)

• The process can be repeated a couple of times to get more robust results (and standard errors).
• It also works for feature groups (shuffling multiple feature columns together).
• One can also study relative increases in average loss.
• It's easy to implement. The only annoying thing is to specify the loss.

Major drawback: Especially when features are collinear/dependent, shuffling creates uncommon or even impossible values. So the model is forced to extrapolate, which may lead to strange results.

It can be calculated by R packages {iml}, {DALEX}, {hstats} (my new one dedicated to Friedman's model-agnostic interaction statistics).

With {hstats}:

library(shapviz)  # to get a cool dataset
library(ggplot2)
library(hstats)

fit <- glm(
SALE_PRC ~ log(LND_SQFOOT) + TOT_LVG_AREA +
log(OCEAN_DIST) + log(CNTR_DIST) + factor(structure_quality),
data = miami,
)

summary(fit)

v <- setdiff(all.vars(formula(fit)), "SALE_PRC")

system.time(  # 0.3 seconds - note the type = "response" of predict()
imp <- perm_importance(
fit, X = miami[v], y = miami$SALE_PRC, loss = "gamma", type = "response" ) ) plot(imp) # Relative increases instead of absolute ones rel_imp <- perm_importance( fit, X = miami[v], y = miami$SALE_PRC, loss = "gamma",
normalize = TRUE, type = "response"
)
plot(rel_imp) +
xlab("Rel. increase in average loss")


Shuffling living area increases the average unit Gamma deviance by almost 300%. (The error bars are created by shuffling four times by default.)

References:

1. Breiman, Leo. 2001. “Random Forests.” Machine Learning 45 (1): 5–32. https://doi.org/10.1023/A:1010933404324.
2. https://christophm.github.io/interpretable-ml-book/feature-importance.html

## SHAP importance

The other frequently used model-agnostic method is based on SHAP values. The SHAP value of feature $$j$$ and observation $$i$$ equals the fair additive contribution of the feature to the prediction. Taking the average absolute SHAP value over a set of observations is a very natural measure of importance. A value of 0.2, e.g., means that the features increases or decreases the prediction on average by 0.2.

The problem with SHAP values: While their game theoretic pendants (Shapley value) are perfect in all desireable aspects, their translation to statistics is not. Model-agnostic algorithms are either slow or imprecise, and like with partial dependence plots (and permutation importance), the model needs to be applied on sometimes unnatural/impossible feature combinations.

• A comment above clarifies some of what I was confused about. If you permute a variable's rows and recompute predictions that is a very fast step when you save the design matrix at the front end. I was thinking you had to get new regression coefficients on permuted data. Thanks. Commented Sep 13, 2023 at 17:08
• This goes into the direction of "global surrogate models". Fit simple model (linear regression or a single tree) to the model predictions of a black box model and then study the "white box" model instead. Its R-squared is used as an indicator of approximation quality. Regarding our discussion above: I am not so experienced with Bayesian models, so calculating predictions as modes from some posterior distribution might indeed be expensive. Commented Sep 13, 2023 at 17:15
• Thanks Michael. The Bayesian procedure would be to repeat the variable importance measure calculation 1000 times over 1000 posterior draws of the regression coefficients, get 1000 importance measures for each predictor, and feed those 1000 through a fast algorithm for getting an uncertainty interval for importances. Commented Sep 13, 2023 at 17:34
• I've learned a lot from this discussion. My current thinking is that random permutation methods are quite general but probably lead to more volatile (higher standard errors) importance estimates. The global surrogate model @MichaelM spoke about is more appealing to me. I'll expand upon that thought with an update to the OP. Commented Sep 14, 2023 at 11:21
• I didn’t see it mentioned in the responses (apologies if it was), but the permutation method does not make much sense when you have multicollinearity issues since it requires individually permuting the individual predictor columns (I.e., you’ll likely end up with invalid or unlikely data instances like age=3 and height=6 ft, which is essentially extrapolating). Commented Sep 19, 2023 at 18:18

Update 4: I've implemented the new general purpose relative explained variation method. Explanation and examples are here. It would be nice to compare its performance with random permutation-based methods.

To perform the comparison, I first wrote an R function which computes the random permutation-based relative explained variation (REV) metric. This metric quantifies the impact of shuffling a given predictor on model performance, indexed by the model linear predictor ($$X\hat{\beta}$$).

rexVar_perm <- function (object, nperm = 20) {
#' compute permutation-based REV from ols and orm objects

# extract predictor indices from anova.rms(fit) (great suggestion by Prof Harrell)
# codes taken from rex() inner function of rexVar()  https://github.com/harrelfe/rms/blob/master/R/rexVar.r
a <- anova(object)
pind <- attributes(a)\$which
rm <- c("TOTAL NONLINEAR","TOTAL NONLINEAR + INTERACTION",
"TOTAL INTERACTION","TOTAL",
" Nonlinear"," All Interactions", "ERROR",
" f(A,B) vs. Af(B) + Bg(A)")
rn <- rownames(a)
rm <- c(rm, rn[substring(rn, 2, 10) == "Nonlinear"])
pind <- pind[rn %nin% rm]
names(pind) <- sub('  (Factor+Higher Order Factors)', '', names(pind), fixed=TRUE)

#' write a nested loop to compute rev of individual predictors and
#  repeat process nperm times

lp_full <- predict(object)
rx <- matrix(NA, nrow=nperm, ncol=length(pind),
dimnames=list(NULL, names(pind)))

for (i in names(pind)) {
for(j in 1:nperm){
X     <- object$$x draws <- object$$coefficients
set.seed(j)  # ensure reproducibility
srows <- sample(seq_len(nrow(X)))
X[, pind[[i]]]   <- X[srows, pind[[i]], drop = FALSE]
lp_shuffle <- matxv(X, draws)
rev <- 1- cor(lp_shuffle, lp_full)
rm(X)
rx[j, i] <- round(rev,3)

}
}

return(apply(rx, 2, mean))

}


Next, I used the MASS package to simulate datasets with different collinearities and compute the corresponding general purpose-based and random permutation-based REVs. Focussing on the simplest case (and to make things easy for myself), each dataset has 2 predictors (v1 and v2) and a continuous outcome (y).

library(rms)
library(here)
library(MASS)
library(ggplot2)
library(data.table)
source("rexVar.R")  ## https://github.com/harrelfe/rms/blob/master/R/rexVar.r
source("rexVar_perm.R")

# write function to generate dataframe with correlated predictors
sim_cdata <- function(cor_val, n = 100){
Sigma <- matrix(c(1, cor_val, cor_val, 1), ncol = 2)
predictors <- mvrnorm(n, mu = c(0, 0), Sigma = Sigma)

# assume predictors share equal predictive power
y <- 1 * predictors[, 1] + 1 * predictors[, 2] +
rnorm(n, mean = 0, sd = 1)
data <- data.frame(v1 = predictors[, 1], v2 = predictors[, 2], y)
return(data)
}


As it turns out and as cautioned by @Michael M and @bgreenwell, when the predictors are strongly collinear (correlation > 0.80 in our example), the random permutation-based REV results become somewhat erratic. In contrast, the general purpose-based REV values decline steadily with increasing collinearity, reflecting the diminishing incremental predictive value that a given predictor has over its (correlated) counterpart.

@Michael M and @bgreenwell, given that datasets with collinearities are not a rarity in biomedical research, could you advise on alternative variable importance measures that are used in the machine learning world?

gen_rev <- function (data){
mod_ols <- ols(y ~ v1 + v2, x=T, data = data)
perm_rev <-   rexVar_perm(mod_ols, nperm = 500)
rexVar_rev <- rexVar(mod_ols, data = data)
r <- c(perm_rev, rexVar_rev)
return(r)
}

cor_seq <- seq(0, 0.98, by = 0.02)

rx <- matrix(NA, nrow= length(cor_seq),
ncol=5,
dimnames=list(NULL,
c("cor_val", "v1_perm","v2_perm",
"v1_rexvar", "v2_rexvar")))
s <- 0
for(i in cor_seq){
s <- s + 1
mydata <- sim_cdata(cor_val = i, n = 500)
rx[s,] <- c(i, gen_rev(mydata))
}
dt <- data.table (rx)
long_dt <- melt(dt, id.vars = "cor_val", variable.name = "variable")
long_dt[, c("variable", "method") := tstrsplit(variable, "_",
fixed = TRUE)]

ggplot(long_dt, aes(x = cor_val, y = value, color = method)) +
geom_line() +
geom_point() +
facet_grid(. ~ variable) +
labs(
x = "Correlation between predictors",
y = "REV value",
title = "Comparison between rexVar and permutation algorithm"
) +
theme_bw(base_size = 13)