Using flashlight

Model performance

How precise are models if applied to unseen data? This aspect is of key interest of basically any supervised machine learning model and helps to identify the best models or, if applied to subgroups, identify problematic segments with low performance.

Variable importance

Which variables are particularly relevant for the model? This aspect is helpful in different ways. Firstly, it might help to simplify the full modelling process by eliminating difficult to assess input variables with low explanatory power. Secondly, its pure information. Thirdly, it might help to identify problems in data structure: if one variable is extremely relevant and all others are not, then there might be some sort of information leakage from the response. Finally, variable importance considerations are very relevant for quality assurance as well.

Different modelling techniques offer different ways of variable importance. In linear models, we consider for example F-test statistics or p values, in tree-based methods the number of splits or average split gains etc. A model agnostic way to assess this is called permutation importance: For each input variable \(X\), its values are randomly shuffled and the drop in performance with respect to a scoring function is calculated. The more important a variable, the larger the drop. If a variable can be shuffled without any impact on model precision, it is completely irrelevant. The method is described in (Fisher, Rudin, and Dominici 2018). Multiple permutations lead to more stable estimates of variable importance and allow the calculation of standard errors.

Effects of input variables

In linear regression, the fitted model consists of an affine linear function in the inputs. Its coefficients immediately tell us how the response is expected to change if the value of one single input variable \(X\) is systematically being adapted. How to describe such effects of a variable \(X\) for more complex models that include non-linearities and high-order interactions?

One approach is to study Individual Conditional Expectation (ICE) profiles of selected observations: They show how predictions of observation \(i\) react when the input variable \(X\) is systematically being changed, see (Goldstein et al. 2015). The more different the profiles are in shape or slope, the stronger are the interaction effects. For a linear regression without interactions, all such profiles would be parallel. Centered ICE profiles (“c-ICE”) help to detect interactions even better.

If many ICE profiles are averaged, we get partial dependence profiles which can be viewed as the average effect of variable \(X\), pooled over all interactions. Partial dependence plots where introduced in Friedman’s seminal 2001 article on gradient boosting (Friedman 2001). Partial dependence profiles (as well as ICE profiles) are not limited to one single variable \(X\). Also multiple variables can be systematically varied on a grid in order to study the multivariate impact on the typical response.

Studying ICE and partial dependence profiles means investigating the effect of \(X\) while holding all other predictors fixed. The more unnatural this Ceteris Paribus assumption is, the less reasonable are the results of ICE and partial dependence. Accumulated local effects or ALE profiles (Apley and Zhu 2016) try to overcome this weakness of ICE and partial dependence considerations.

ALE profiles at positions \(x_i\) are calculated as follows:

1. First, left derivatives \(\Delta_i\) are estimated for all \(i\): Calculate slope of partial dependence line from \(x_i-\varepsilon\) to \(x\) based on data points in \([x_i-\varepsilon, x_i]\).

2. The uncalibrated ALE value at \(x_i\) is the (integrated/accumulated) partial sum \(\sum_{j \le i} \Delta_j\).

3. In the calibration step, the effects are shifted to the level of the response variable or its predictions.

An alternative to partial dependence and ALE profiles is to look at combined effects of \(X\) including the effects from all other predictors. Such effects are estimated by averaging the predictions within values of the predictor \(X\) of interest. If such prediction profile differs considerably from the observed average response, this might be a sign of model underfit. This visualization is sometimes called “marginal plot” or “M plot,” see (Apley and Zhu 2016). Residual profiles immediately show such misfits as well. In classical statistical modelling, this sort of plots are called “fitted versus covariable” plots or “residual versus covariable” plots.

If SHAP values have been calculated, plotting their average contributions of \(X\) against \(X\) is an elegant way to visualize the effect of \(X\), see Subsection “Variable contribution breakdown for single observations.”

Besides looking at average profiles, it is often also revealing to consider quartiles or to visualize partial dependence, response and prediction profiles in the same plot.

Instead of studying profiles, it might be revealing to look at scatter plots of responses, predictions, residuals, or SHAP values versus \(X\).

Interaction strength

Besides measuring overall variable importance, an interesting aspect is to measure the strength of non-additivity associated with each covariable (i.e. the overall interaction strength) and/or between pairs of covariables. A model-agnostic way to assess such measures is based on partial dependence curves, see (Friedman and Popescu 2008). A computationally less demanding alternative to their H-statistic of measuring total interaction strength per covariable can be extracted from centered ICE: How much of their variability is unexplained by the main effect?

Variable contribution breakdown for single observations

Instead of studying global variable importance and effects, there are different techniques to entangle how a single prediction can be explained: LIME, LIVE, SHAP and breakdown, see (Molnar 2019) and (Gosiewska and Biecek 2019) for an overview. The flashlight package currently supports the breakdown method as well as approximate SHAP. These two methods are based on additive decomposition of a prediction.

The breakdown algorithm works as follows: First, the visit order \((x_1, ..., x_m)\) of variables is specified. Then, in the query data, the column \(x_1\) is set to the constant value of \(x_1\) of the observation to be explained. The change in the (weighted) average predicted value on the query data measures the contribution of \(x_1\) on the prediction. This procedure is iterated over all \(x_i\) until eventually, all rows in the query data are identical to the observation to be explained.

A complication with this approach is that the visit order is relevant, at least for non-additive models. Ideally, the algorithm could be repeated for all \(m!\) possible visit orders and its results averaged per variable. This is what SHAP does, see e.g. (Gosiewska and Biecek 2019) for an explanation. Unfortunately, there is no efficient way to do this in a model agnostic way. Thus we offer two approximations. The first one is the short-cut described in (Gosiewska and Biecek 2019) and implemented in the ibreakdown package (https://CRAN.R-project.org/package=iBreakDown): There, the variables \(x_i\) are sorted by the size of their contribution in the same way as the breakdown algorithm but without iteration, i.e. starting from the original query data for each variable \(x_i\). We call this approach “breakdown.” The second, computationally more intensive way to approximate SHAP values is based on running the breakdown algorithm on a small number (e.g. 12) random permutations and then averaging the results. We call this approach (approximate) “SHAP.”

Global surrogate models

As suggested in (Molnar 2019), one way to explain a complex model is to fit an easy to interprete model like a decision tree to the predictions of the complex model and then study this simple surrogate model. The quality of its approximation can be measured by the R-squared of the surrogate model.

flashlight

The flashlight package offers these tools in a very simple way.

Data preparation

Let’s do some data preparation common for all models under consideration.

prep <- mutate(
  house_prices,
  log_price = log(price),
  grade = as.integer(as.character(grade)),
  year = as.numeric(format(date, '%Y')),
  age = year - yr_built,
  year = factor(year),
  zipcode = as.factor(as.character(zipcode)),
  waterfront = factor(waterfront, levels = c(FALSE, TRUE), 
                      labels = c("no", "yes")))

x <- c("grade", "year", "age", "sqft_living", "sqft_lot", "zipcode",
       "condition", "waterfront")

Modelling

The random forest can directly work with this data structure. However, for the linear model, we need a small function with additional feature engineering. We log transform some input variables and categorize the zipcode in large groups. Similarly, for XGBoost, such a wrapper function turns non-numeric input variables to numeric. We will make use of these functions for both data preparation and prediction.

# Data wrapper for the linear model
prep_lm <- function(data) {
  data %>%
    mutate(sqrt_living = log(sqft_living),
           sqrt_lot = log(sqft_lot))
}

# Data wrapper for xgboost
prep_xgb <- function(data, x) {
  data %>%
    select_at(x) %>%
    mutate_if(Negate(is.numeric), as.integer) %>%
    data.matrix()
}

Then, we split the data and train our models.

# Train / valid / test split (70% / 20% / 10%)
set.seed(56745)
ind <- sample(10, nrow(prep), replace = TRUE)

train <- prep[ind >= 4, ]
valid <- prep[ind %in% 2:3, ]
test <- prep[ind == 1, ]

(form <- reformulate(x, "log_price"))
#> log_price ~ grade + year + age + sqft_living + sqft_lot + zipcode + 
#>     condition + waterfront
fit_lm <- lm(update.formula(form, . ~ . + I(sqft_living^2)), 
             data = prep_lm(train))

# Random forest
fit_rf <- ranger(form, data = train, respect.unordered.factors = TRUE, 
                 num.trees = 100, seed = 8373)
cat("R-squared OOB:", fit_rf$r.squared)
#> R-squared OOB: 0.8654302

# Gradient boosting
if (has_xgb) {
  dtrain <- xgboost::xgb.DMatrix(prep_xgb(train, x), 
                                 label = train[["log_price"]])
  dvalid <- xgboost::xgb.DMatrix(prep_xgb(valid, x), 
                                 label = valid[["log_price"]])
  
  params <- list(learning_rate = 0.2,
                 max_depth = 5,
                 alpha = 1,
                 lambda = 1,
                 colsample_bytree = 0.8)
  
  fit_xgb <- xgboost::xgb.train(
    params,
    data = dtrain,
    watchlist = list(train = dtrain, valid = dvalid),
    nrounds = 250,
    print_every_n = 25,
    objective = "reg:squarederror"
  )
} else {
  fit_xgb <- rpart(form, data = train, 
                   control = list(xval = 0, cp = -1, minsplit = 100))
  prep_xgb <- function(data, x) {
    data
  }
}
#> [1]  train-rmse:10.048225    valid-rmse:10.056726 
#> [26] train-rmse:0.240650 valid-rmse:0.249891 
#> [51] train-rmse:0.191478 valid-rmse:0.205166 
#> [76] train-rmse:0.177957 valid-rmse:0.195167 
#> [101]    train-rmse:0.169047 valid-rmse:0.189242 
#> [126]    train-rmse:0.162973 valid-rmse:0.186441 
#> [151]    train-rmse:0.158342 valid-rmse:0.184533 
#> [176]    train-rmse:0.154759 valid-rmse:0.183916 
#> [201]    train-rmse:0.151803 valid-rmse:0.183548 
#> [226]    train-rmse:0.149121 valid-rmse:0.182889 
#> [250]    train-rmse:0.146648 valid-rmse:0.182994

Creating the flashlights

Let’s initialize one flashlight per model. Thanks to individual prediction functions, any model can be used in flashlight, even keras models.

fl_mean <- flashlight(
  model = mean(train$log_price), 
  label = "mean",
  predict_function = function(mod, X) rep(mod, nrow(X))
)
fl_lm <- flashlight(
  model = fit_lm, 
  label = "lm",
  predict_function = function(mod, X) predict(mod, prep_lm(X))
)
fl_rf <- flashlight(
  model = fit_rf, 
  label = "rf",
  predict_function = function(mod, X) predict(mod, X)$predictions
)
fl_xgb <- flashlight(
  model = fit_xgb, 
  label = "xgb",
  predict_function = function(mod, X) predict(mod, prep_xgb(X, x))
)
print(fl_xgb)
#> 
#> Flashlight xgb 
#> 
#> Model:            Yes
#> y:            No
#> w:            No
#> by:           No
#> data dim:         No
#> predict_fct default:  FALSE
#> linkinv default:  TRUE
#> metrics:      rmse
#> SHAP:             No

What about all other relevant elements of a flashlight like the underlying data, the response name, metrics, retransformation functions etc? We could pass them to each of our flashlights. Or better, we can combine the flashlights to a multiflashlight and pass additional common arguments there.

fls <- multiflashlight(
  list(fl_mean, fl_lm, fl_rf, fl_xgb), 
  y = "log_price", 
  linkinv = exp,
  data = valid, 
  metrics = list(rmse = rmse, `R-squared` = r_squared)
)

We could even extract these completed flashlights from the multiflashlight as if the latter is a list (actually it is a list with additional class multiflashlight).

fl_lm <- fls$lm

Assess performance

Let’s compare the models regarding their validation performance.

perf <- light_performance(fls)
perf
#> 
#> I am an object of class light_performance_multi 
#> 
#> data.frames:
#> 
#>  data 
#> # A tibble: 8 x 3
#>   label metric        value
#>   <fct> <fct>         <dbl>
#> 1 mean  rmse       0.527   
#> 2 mean  R-squared -0.000267
#> 3 lm    rmse       0.188   
#> 4 lm    R-squared  0.872   
#> 5 rf    rmse       0.192   
#> 6 rf    R-squared  0.868   
#> 7 xgb   rmse       0.183   
#> 8 xgb   R-squared  0.879
plot(perf)

plot(perf, fill = "darkred") +
  xlab(element_blank())

The plot “politics” of flashlight is to provide simple graphics with minimal ggplot-tuning, so you are able to add your own modifications. If you are completely unhappy about the proposed plot (e.g. rather favour a scatterplot over a barplot), extract the data slot of perf and create the figure from scratch:

head(perf$data)
#> # A tibble: 6 x 3
#>   label metric        value
#>   <fct> <fct>         <dbl>
#> 1 mean  rmse       0.527   
#> 2 mean  R-squared -0.000267
#> 3 lm    rmse       0.188   
#> 4 lm    R-squared  0.872   
#> 5 rf    rmse       0.192   
#> 6 rf    R-squared  0.868

perf$data %>%
  ggplot(aes(x = label, y = value, group = metric, color = metric)) +
  geom_point()

The same logic holds for all other main functions in the flashlight package.

For performance considerations, the minimum required info in the (multi-)flashlight are: “y,” “predict_function,” “model,” “data” and “metrics.” The latter two can also be passed on the fly.

Variable importance

Now let’s study variable importance of the explainers.

(imp <- light_importance(fls, v = x))
#> 
#> I am an object of class light_importance_multi 
#> 
#> data.frames:
#> 
#>  data 
#> # A tibble: 32 x 5
#>    label metric variable      value error
#>    <fct> <fct>  <chr>         <dbl> <lgl>
#>  1 mean  rmse   grade       0       NA   
#>  2 mean  rmse   year        0       NA   
#>  3 mean  rmse   age         0       NA   
#>  4 mean  rmse   sqft_living 0       NA   
#>  5 mean  rmse   sqft_lot    0       NA   
#>  6 mean  rmse   zipcode     0       NA   
#>  7 mean  rmse   condition   0       NA   
#>  8 mean  rmse   waterfront  0       NA   
#>  9 lm    rmse   grade       0.0828  NA   
#> 10 lm    rmse   year        0.00187 NA   
#> # ... with 22 more rows
plot(imp, fill = "darkred")

Note for permutation importance: If the calculation takes too long (e.g. large query data), set n_max to some reasonable value. light_importance will then randomly pick n_max rows and use only these for assessment of importance. On the other hand, if the calculations is very fast, consider setting m_repetitions to a value > 1 in order to repeat the permutations multiple times. The resulting variable importance corresponds to the average drop in performance and standard errors are added.

If SHAP values have been precomputed by add_shap, use type = "shap" to show mean absolute shap values per covariable.

Individual conditional expectation

How do predictions change when sqft_living changes alone? We can investigate this question by looking at “Individual Conditional Expectation” (ICE) profiles of a couple of observations.

cp <- light_ice(fls, v = "sqft_living", n_max = 30, seed = 35)
plot(cp, alpha = 0.2)

Note: Setting seed to a fixed value will ensure that the flashlights will consider the same rows. An alternative would be to pass a small subset of the data to light_ice and calculate all profiles or by passing row indices through indices for fixed selection.

Centered ICE profiles (“c-ICE”) can help to increase visibility of interactions.

cp <- light_ice(fls, v = "sqft_living", n_max = 30, seed = 35, center = "first")
plot(cp, alpha = 0.2)

Partial dependence profiles

If many ICE profiles (in our case 1000) are averaged, we get an impression on the average effect of the considered variable. Such curves are called partial dependence profiles (PD) resp. partial dependence plots.

pd <- light_profile(fls, v = "sqft_living")
pd
#> 
#> I am an object of class light_profile_multi 
#> 
#> data.frames:
#> 
#>  data 
#> # A tibble: 40 x 5
#>    sqft_living counts   value label type              
#>          <dbl>  <int>   <dbl> <fct> <fct>             
#>  1         500   1000 462845. mean  partial dependence
#>  2        1500   1000 462845. mean  partial dependence
#>  3        2500   1000 462845. mean  partial dependence
#>  4        3500   1000 462845. mean  partial dependence
#>  5        4500   1000 462845. mean  partial dependence
#>  6        5500   1000 462845. mean  partial dependence
#>  7        6500   1000 462845. mean  partial dependence
#>  8        7500   1000 462845. mean  partial dependence
#>  9        8500   1000 462845. mean  partial dependence
#> 10        9500   1000 462845. mean  partial dependence
#> # ... with 30 more rows
plot(pd)

The light_profile function offers different ways to specify the evaluation points of the profiles. One option is to explicitly pass such points.

pd <- light_profile(fls, v = "sqft_living", 
                    pd_evaluate_at = seq(1000, 4000, by = 100))
plot(pd)

Two-dimensional:

pd <- light_profile2d(fls, v = c("condition", "grade"))
plot(pd)

Accumulated local effects (ALE)

An approximation of main effects without Ceteris Paribus clause are accumulated local effects profiles (Fisher, Rudin, and Dominici 2018). They are based on accumulating local partial dependence slopes.

ale <- light_profile(fls, v = "sqft_living", type = "ale")
ale
#> 
#> I am an object of class light_profile_multi 
#> 
#> data.frames:
#> 
#>  data 
#> # A tibble: 32 x 5
#>    sqft_living counts   value label type 
#>          <dbl>  <int>   <dbl> <fct> <fct>
#>  1        1500   1000 462845. mean  ale  
#>  2        2500   1000 462845. mean  ale  
#>  3        3500    859 462845. mean  ale  
#>  4        4500    231 462845. mean  ale  
#>  5        5500     52 462845. mean  ale  
#>  6        6500     16 462845. mean  ale  
#>  7        7500      9 462845. mean  ale  
#>  8        8500      5 462845. mean  ale  
#>  9        1500   1000 360940. lm    ale  
#> 10        2500   1000 483272. lm    ale  
#> # ... with 22 more rows
plot(ale)

Note: While equally sized x-breaks are easy to read, quantile binning usually leads to more stable results.

plot(light_profile(fls, v = "sqft_living", type = "ale", cut_type = "quantile"))

In order to calculate ICEs, PDs and ALEs, the following elements need to be available in the (multi-)flashlight: “predict_function,” “model,” “linkinv” and “data.” “data” can also be passed on the fly.

Profiles of predicted values, residuals, and response

We can use the function light_profile not only to create partial dependence profiles but also to get profiles of predicted values (“M plots”), responses or residuals. Additionally, we can either use averages or quartiles as summary statistics.

Average predicted values versus the living area are as follows:

format_y <- function(x) format(x, big.mark = "'", scientific = FALSE)

pvp <- light_profile(fls, v = "sqft_living", type = "predicted", 
                     format = "fg", big.mark = "'")
plot(pvp) +
  scale_y_continuous(labels = format_y)

Note the formatting of y values as well as the formatC option format = "fg" and big.mark passed to the constructor of the x labels in order to improve basic appearance. We will recycle some of these settings for the next plots.

Similarly, the average response profiles (identical for all flashlights, to we only show one of them):

rvp <- light_profile(fl_lm, v = "sqft_living", type = "response", format = "fg")
plot(rvp) +
  scale_y_continuous(labels = format_y)

Same, but quartiles:

rvp <- light_profile(fl_lm, v = "sqft_living", type = "response",
                     stats = "quartiles", format = "fg")
plot(rvp) +
  scale_y_continuous(labels = format_y)

What about residuals? First, we remove the “mean” flashlight by setting it NULL.

fls$mean <- NULL
rvp <- light_profile(fls, v = "sqft_living", type = "residual",
                     stats = "quartiles", format = "fg")
plot(rvp) +
  scale_y_continuous(labels = format_y)

While the tree-based models have smaller residuals and medians close to 0, the linear model shows residual curvature that could be captured by representing sqft_living by more parameters.

If your are unhappy about the “group by” strategy, set swap_dim to TRUE.

plot(rvp, swap_dim = TRUE) +
  scale_y_continuous(labels = format_y)

For less bars, set n_bins in light_profile:

rvp <- light_profile(fls, v = "sqft_living", type = "residual",
                     stats = "quartiles", format = "fg", n_bins = 5)
plot(rvp, swap_dim = TRUE) +
  scale_y_continuous(labels = format_y)

Just as diverging ICE profiles give you a clou about the presence of interactions, we can use the option stats = "quartiles" (with pd_center) to show divergence of the mean-centered ICE profiles as boxes (predictions at log-scale to suppress interaction-like effects of the retransformation function):

rvp <- light_profile(fls, v = "sqft_living", use_linkinv = FALSE,
                     stats = "quartiles", pd_center = "mean")
plot(rvp)

For prediction profiles, the same elements as for ICE/PDs are required, while for response profiles we need “y,” “linkinv” and “data.” “data” can also be passed on the fly.

Visualizing different types of profiles as “effects” plot

In assessing the model quality, it is often useful to visualize

• response profile (quartiles or means),

• average predictions, and

• model effects (partial dependence profiles, or accumulated local effects)

in the same plot and for each input variable. The flashlight package offers the function light_effects to combine such profile plots:

eff <- light_effects(fl_lm, v = "condition")
p <- plot(eff) +
  scale_y_continuous(labels = format_y)
p

Let’s add counts to see if the gaps between response and predicted profiles are problematic or just due to small samples.

plot_counts(p, eff, alpha = 0.2)

The biggest gaps occur at very rare conditions, so the model looks quite fine.

Note: Due to retransformation from log scale, the response profile is slightly higher than the profile of predicted values. If we would evaluate on the modelled log scale, that gap would vanish.

eff <- light_effects(flashlight(fl_lm, linkinv = I), v = "condition")
p <- plot(eff, use = "all") +
  scale_y_continuous(labels = format_y) +
  ggtitle("Effects plot on modelled log scale")
p

Besides adding counts to the figure, representing observed responses as boxplots (no whiskers and outliers in order to avoid too large y scale) might help to judge if there is a problematic misfit.

eff <- light_effects(fl_lm, v = "condition", stats = "quartiles")
p <- plot(eff, rotate_x = FALSE) +
   scale_y_continuous(labels = format_y)
plot_counts(p, eff, fill = "blue", alpha = 0.2, width = 0.3)

The plot method of light_effects allows to hide certain plot elements if it looks too dense. If we would e.g. want to compare partial dependence and accumulated local effects plots, we can use the following approach.

eff <- light_effects(fls, v = "sqft_living", v_labels = FALSE, 
                     cut_type = "quantile", n_bins = 25)
plot(eff, use = c("pd", "ale"), show_points = FALSE) + 
  scale_y_continuous(labels = format_y) +
  coord_cartesian(xlim = c(1000, 4000), ylim = c(3e5, 9e5))

As an alternative to profile plots, we could also consider scatter plots.

pr <- light_scatter(fls, v = "sqft_living", n_max = 300)
plot(pr, alpha = 0.2)

Interaction strength

Based on Friedman’s H statistic (Friedman and Popescu 2008), we can investigate which variable has strongest interaction effect relative to its full effect.

st <- light_interaction(fls, v = x, grid_size = 30, n_max = 50, seed = 42)
plot(st)

Since we have not used interaction terms in the linear model, all values are zero there.

Pairwise interactions are expensive to compute as each variable pair has to be assessed. We limit the search for interactions to the four variables with highest relative interaction strength.

st_pair <- light_interaction(fls, v = most_important(st, 4), 
                             pairwise = TRUE, n_max = 50, seed = 42)
plot(st_pair)

Variable contribution breakdown for single observation

Besides investigating global effects, we can use light_breakdown to calculate variable impact on one single (log) prediction.

bd <- light_breakdown(fl_lm, new_obs = valid[1, ], 
                      v = x, n_max = 1000, seed = 74)
plot(bd, size = 3)

We have set the size of the reference data set to n_max = 1000 in order to save time. The only variable with positive impact on the prediction is “age.” The other variables have a negative impact compared to the 1000 reference observations.

Global surrogate trees

A very different, but intuitive approach is to fit the original predictions of the models by simple decision trees and then visualize these trees.

surr <- light_global_surrogate(fls$xgb, v = x)
print(surr$data$r_squared)
#> [1] 0.7001643
plot(surr)

Even if it is extremely simple to explain, the surrogate model is not a bad approximation (see R-squared).