- That they provide information about the model.
- That they are highly comprehensible.

In attempts to white-box complex machine learning models, one dichotomy stands out:

**Global models**seek to determine the relative role of features in the construction of the predictions once the model has been trained. This is done at the global level, so that the patterns that are shown in the interpretation hold on average over the whole training set.**Local models**aim to characterize how the model behaves around one particular instance by considering small variations around this instance. The way these variations are processed by the original model allows to simplify it by approximating it, e.g., in a linear fashion. This approximation can for example determine the sign and magnitude of the impact of each relevant feature in the vicinity of the original instance.

Let us start with the simplest example of all. In a linear model,

\begin{align}
y_i=\alpha+\sum_{k=1}^K\beta_kx_i^k+\epsilon_i,
\end{align}

the following elements are usually extracted from the estimation of the $\beta_k$:

- the $R^2$, which appreciates the
**global fit**of the model (possibly penalized to prevent overfitting with many regressors). The $R^2$ is usually computed in-sample; - the sign of the estimates $\hat{\beta}_k$, which indicates the
**direction**of the impact of each feature $x^k$ on $y$; - the $t$-statistics $t_{\hat{\beta_k}}$, which evaluate the
**magnitude**of this impact: regardless of its direction, large statistics in absolute value reveal prominent variables. Often, the $t$-statistics are translated into $p$-values which are computed under some suitable distributional assumptions.

**surrogate** models. The process is simple:

\begin{align}
\hat{f}(\textbf{X})=g(\textbf{X})+\textbf{error}
\end{align}

In [2]:

```
new_target = fit_RF.predict(X_short) # saving the predictions of tree as new target
decision_tree_model = tree.DecisionTreeRegressor(max_depth=3) # defining the global interpretable tree surrogate model
TreeSurrogate=decision_tree_model.fit(X_short,new_target) # fitting the surrogate
fig, ax = plt.subplots(figsize=(13, 8)) # setting the chart parameters
tree.plot_tree(TreeSurrogate,feature_names=features_short, ax=ax)
plt.show() # Plot!
```

FIGURE 13.1: Example of surrogate tree.

\begin{align}
I(k)=\sum_{n\in \mathcal{N}_k}G(n),
\end{align}

In [5]:

```
tree_VI = pd.DataFrame(data=fit_tree.feature_importances_,index=features_short,columns=['Tree']) # VI from tree model
RF_VI = pd.DataFrame(data=fit_RF.feature_importances_,index=features_short,columns=['RF']) # VI from random forest
XGB_VI = pd.DataFrame(data=fit_xgb.feature_importances_,index=features_short,columns=['XGB']) # VI from boosted trees
VI_trees=pd.concat([tree_VI,RF_VI,XGB_VI],axis=1) # Aggregate the VIs
VI_trees=VI_trees.loc[['Mkt_Cap_12M_Usd','Pb','Vol1Y_Usd']]/np.sum(VI_trees.loc[['Mkt_Cap_12M_Usd','Pb','Vol1Y_Usd']])
VI_trees.plot.bar(figsize=[10,6]) # Plotting sequence
```

Out[5]:

<AxesSubplot:>

FIGURE 13.2: Variable importance for tree-based models.

The baseline method to assess feature importance in the general case is the following:

In [6]:

```
from sklearn.linear_model import Ridge
import random
y_penalized = data_ml['R1M_Usd'].values # Dependent variable
X_penalized = data_ml[features].values # Predictors
fit_ridge_0 = Ridge(alpha=0.01) # Trained model
fit_ridge_0.fit(X_penalized_train, y_penalized_train) # Fit model
l_star= np.mean(np.square(fit_ridge_0.predict(X_penalized_train)-y_penalized_train)) # Loss
```

In [12]:

```
from collections import Counter
res = [] # Initialize
feature_random = random.sample(list((Counter(features)-Counter(features_short)).elements()), 12) # selecting fewer features for computation time sake
for feat in (features_short+feature_random): # Loop on the features
temp_data=training_sample[features].copy() # temp dataframe for all features
temp_data.loc[:,feat] = np.random.permutation(training_sample[feat]) # shuffling for feat[i]
fit_ridge_0.fit(temp_data[features], training_sample['R1M_Usd']) # fitting the model
result_VI=pd.DataFrame([feat],columns=['feat'])
result_VI['loss']=[np.mean(np.square(fit_ridge_0.predict(temp_data[features])-training_sample['R1M_Usd'])) - l_star] # Loss
res.append(result_VI) # appending through features loop
res = pd.concat(res)
res.set_index('feat',inplace=True)
```

Finally, we plot the results.

In [14]:

```
res.plot.bar(figsize=[10,6])
```

Out[14]:

<AxesSubplot:xlabel='feat'>

FIGURE 13.3: Variable importance for a ridge regression model.

**average impact** of $k$ on the predictions of the trained model $\hat{f}$. In order to do so, we assume that the feature space is random and we split it in two: $k$ versus $−k$, which stands for all features except for $k$. The partial dependence plot is defined as

\begin{equation}
\tag{13.1}
\bar{f}_k(x_k)=\mathbb{E}[\hat{f}(\textbf{x}_{-k},x_k)]=\int \hat{f}(\textbf{x}_{-k},x_k)d\mathbb{P}_{-k}(\textbf{x}_{-k}),
\end{equation}

d
P
−
k
(
⋅
)
is the (multivariate) distribution of the non-
k
features

x
−
k
. The above function takes the feature values

x
k
as argument and keeps all other features frozen via their sample distributions: this shows the impact of feature

k
solely. In practice, the average is evaluated using Monte-Carlo simulations:

\begin{equation}
\tag{13.2}
\bar{f}_k(x_k)\approx \frac{1}{M}\sum_{m=1}^M\hat{f}\left(x_k,\textbf{x}_{-k}^{(m)}\right),
\end{equation}

where

x
(
m
)
−
k
are independent samples of the non-
k
features.

In [15]:

```
from sklearn.inspection import PartialDependenceDisplay
PartialDependenceDisplay.from_estimator(fit_RF,training_sample[features_short], ['Pb'],kind='average')
```

Out[15]:

<sklearn.inspection._plot.partial_dependence.PartialDependenceDisplay at 0x2abb7384df0>

FIGURE 13.4: Partial dependence plot for the price-to-book ratio on the random forest model.

**simple interpretability**, which implies a limited number of variables with visual or textual representation. This is to make sure any human can easily understand the outcome of the tool;**local faithfulness**: the explanation holds for the vicinity of the instance.

\begin{align}
\xi(x)=\underset{g \in G}{\text{argmin}} \, \mathcal{L}(f,g,\pi_x)+\Omega(g),
\end{align}

\begin{align}
\mathcal{L}(f,g,\pi_x)=\sum_z \pi_x(z)(f(z)-g(z))^2
\end{align}

In Figure 13.5, we provide a simplified diagram of how LIME works.

We proceed with an example of implementation. There are several steps:

- Fit a model on some training data.
- Wrap everything using the lime() function.
- Focus on a few predictors and see their impact over a few particular instances (via the explain() function).

We start with the first step. This time, we work with a boosted tree model.

In [16]:

```
import lime # Package for LIME interpretation
import lime.lime_tabular
xgb_model = xgboost.XGBRegressor( # Parameters of the boosted tree
max_depth=5, # Max depth of each tree
learning_rate=0.5, # Learning rate
objective='reg:squarederror', # booster type of objective function
subsample=1, # Proportion of instance to be sampled (1 = all)
colsample_bytree=1, # Proportion of predictors to be sampled (1 = all)
gamma=0.1, # Penalization
n_estimators=10, # Number of trees
min_child_weight=10) # Min number of instances in each node
xgb_model.fit(train_features_xgb, train_label_xgb) # Training of the model
```

Out[16]:

XGBRegressor(base_score=0.5, booster='gbtree', callbacks=None, colsample_bylevel=1, colsample_bynode=1, colsample_bytree=1, early_stopping_rounds=None, enable_categorical=False, eval_metric=None, feature_types=None, gamma=0.1, gpu_id=-1, grow_policy='depthwise', importance_type=None, interaction_constraints='', learning_rate=0.5, max_bin=256, max_cat_threshold=64, max_cat_to_onehot=4, max_delta_step=0, max_depth=5, max_leaves=0, min_child_weight=10, missing=nan, monotone_constraints='()', n_estimators=10, n_jobs=0, num_parallel_tree=1, predictor='auto', random_state=0, ...)

Then, we head on to steps two and three.

In [19]:

```
explainer = lime.lime_tabular.LimeTabularExplainer(train_features_xgb.values, # values in tabular i.e. Matrix
mode='regression', # “classification” or “regression”
feature_names=train_features_xgb.columns,
verbose=1) # if true, print local prediction values from linear model
exp = explainer.explain_instance(train_features_xgb.iloc[0,:].values, # First instance in train_sample
predict_fn=xgb_model.predict, # prediction function
labels=train_label_xgb.iloc[0].values, # iterable with labels to be explained
distance_metric='euclidean', # Dist.func. "gower" is one alternative
num_samples=900, # Nb samples for loss function
num_features=6) # Nb of features shown (important ones)
exp.show_in_notebook(show_table=True) # Visual display
```

Intercept 0.0133735464008345 Prediction_local [0.01957318] Right: 0.020951485