Feature Ranking

Feature Ranking#

Michael J. Pyrcz, Professor, The University of Texas at Austin

Chapter of e-book “Applied Machine Learning in Python: a Hands-on Guide with Code”.

Cite this e-Book as:

Pyrcz, M.J., 2024, Applied Machine Learning in Python: A Hands-on Guide with Code [e-book]. Zenodo. doi:10.5281/zenodo.15169139

The workflows in this book and more are available here:

Cite the MachineLearningDemos GitHub Repository as:

Pyrcz, M.J., 2024, MachineLearningDemos: Python Machine Learning Demonstration Workflows Repository (0.0.3) [Software]. Zenodo. DOI: 10.5281/zenodo.13835312. GitHub repository: GeostatsGuy/MachineLearningDemos

By Michael J. Pyrcz
© Copyright 2024.

This chapter is a tutorial for / demonstration of Feature Ranking.

YouTube Lecture: check out my lectures on:

These lectures are all part of my Machine Learning Course on YouTube with linked well-documented Python workflows and interactive dashboards. My goal is to share accessible, actionable, and repeatable educational content. If you want to know about my motivation, check out Michael’s Story.

Motivation for Feature Ranking#

There are often many predictor features (input variables) available for us to work with when building our prediction models.

there are good reasons to be selective; throwing in every possible feature is not a good idea!

In general, for the best prediction model, careful selection of the fewest features that provide the most information is a good practice.

Here’s why:

blunders – more predictor features result in more complicated workflows that require more professional time and provide more opportunities for mistakes in the workflow
difficult to visualize – higher-dimensional models, i.e., models with a larger number of predictor features, are more difficult to visualize
model checking – more complicated models may be more difficult to interrogate, interpret, and QC
predictor feature redundancy – more predictor features increase the likelihood of redundant features. The inclusion of highly redundant, collinear, or multicollinear features can increase model variance, increase model instability, and decrease prediction accuracy on testing data
computational time – more predictor features generally increase the computational time required to train the model and the computational storage required; i.e., the model may be less compact and portable
model overfit – the risk of overfit increases with the number of features due to increased model complexity
model extrapolation – many predictor features result in a high-dimensional model space with less data coverage and a greater likelihood of model extrapolation, which may be inaccurate

The primary concern with many predictor features is the curse of dimensionality. Let’s summarize the curse!

Curse of Dimensionality#

Data and Model Visualization - we can visualize projections/slices of higher-dimensional data, but we cannot directly visualize the full geometry beyond 3D, i.e., assess the model fit to the data, evaluate interpolation vs. extrapolation.

consider a 5D example shown as a matrix scatter plot, even in this case there is an extreme marginalization to 2D for each plot,

Example 5D data as a matrix scatter plot.

Sampling - the number of samples sufficient to infer statistics like the joint probability, $P(x_1,\ldots,x_m)$.

recall the calculation of a histogram or normalized histogram: we establish bins and calculate frequencies or probabilities in each bin.
we require a nominal number of data samples for each bin, so we require $𝑛=𝑛_{𝑠/𝑏𝑖𝑛} \cdot 𝑛_{𝑏𝑖𝑛𝑠}$ samples in 1D
but in mD we required $n$ samples to calculate the discretized joint probability,

\[ n=n_{s/bin} \cdot n_{bins}^m$ \]

for example, 10 samples per bin with 35 bins requires 12,250 samples in 2D, and 428,750 samples in 3D

Example 2D data with 35 bins for each feature.

Sample Coverage - the range of the sample values cover the predictor feature space.

fraction of the possible solution space that is sampled, for 1 feature we assume 80% coverage
remember, we usually, directly sample only $\frac{1}{10^7}$ of the volume of the subsurface
yes, the concept of coverage is subjective, how much data to cover? What about gaps? etc.

now if there is 80% coverage for 2 features the 2D coverage is 64%

coverage is,

\[ c = c_1^m \]

Distorted Space - high dimensional space is distorted.

take the ratio of the volume of an inscribed hypersphere in a hypercube,

\[ \frac{\pi^{\frac{m}{2}}}{m 2^{m-1} \Gamma\left(\frac{m}{2}\right)} \to 0 \quad \text{as} \quad m \to \infty \]

note, $\Gamma(𝑛)=(𝑛−1)!$.
high dimensional space is all corners and no middle and most of high dimensional space is far from the middle (all corners!).
as a result distances in high dimensional space lose sensitivity, i.e., for any random points in the space the expected pairwise distances all become the same,

\[ \lim_{m \to \infty} \left( \mathbb{E}\left[\text{dist}_{\text{max}}(m) - \text{dist}_{\text{min}}(m)\right] \right) \to 0 \]

the limit of the expectation of the range of pairwise distances over random points in hyper-dimensional space tends to zero. If distances are almost all the same, Euclidian distance is no longer meaningful!

The ratio of the volume of a hypersphere within a hypercube.

here’s the severity of the distortion for various dimensionalities,

m	nD / 2D
2	1.0
5	0.28
10	0.003
20	0.00000003

Multicollinearity - higher dimensional datasets are more likely to have collinearity or multicollinearity.

Feature linearly described by other features resulting in high model variance.

What is a Good Predictor Feature?#

Before we start feature ranking, let’s define a good predictor feature. There are two important components:

High Relevancy — a strong relationship with the response feature; i.e., the predictor feature is related to the feature that we are attempting to predict.
Low Redundancy — weak relationships with the other predictor features; i.e., the predictor feature provides unique information about the feature that we are attempting to predict.

What is Feature Ranking?#

Feature ranking is a set of methods that assign relative importance or value to each predictor feature with respect to the information it provides for inference and its importance for predicting a response feature. There are a wide variety of possible methods to accomplish this.

For feature ranking and selection,

my recommendation is a ‘wide-array’ approach with multiple analyses and metrics, while understanding the assumptions and limitations of each method.

Here are the general types of metrics that we consider for feature ranking.

Visual Inspection of Data Distributions and Scatter Plots
Statistical Summaries
Model-based
Recursive Feature Elimination

Also, we should not neglect expert knowledge, i.e.,

if additional information is known about physical processes, causation, reliability, and availability of predictor features, this should be integrated into assigning feature ranking.

Load the Required Libraries#

The following code loads the required libraries.

import geostatspy.GSLIB as GSLIB                              # GSLIB utilities, visualization and wrapper
import geostatspy.geostats as geostats                        # GSLIB methods convert to Python  
import geostatspy
print('GeostatsPy version: ' + str(geostatspy.__version__))

GeostatsPy version: 0.0.79

We will also need some standard packages. These should have been installed with Anaconda 3.

ignore_warnings = True                                        # ignore warnings?
import numpy as np                                            # ndarrays for gridded data
import pandas as pd                                           # DataFrames for tabular data
from sklearn import preprocessing                             # remove encoding error
from sklearn.feature_selection import RFE                     # for recursive feature selection
from sklearn.feature_selection import mutual_info_regression  # mutual information
from sklearn.linear_model import LinearRegression             # linear regression model
from sklearn.ensemble import RandomForestRegressor            # model-based feature importance
from sklearn.model_selection import GridSearchCV              # hyperparameter tuning
from sklearn import metrics                                   # measures to check our models
from statsmodels.stats.outliers_influence import variance_inflation_factor # variance inflation factor
import os                                                     # set working directory, run executables
import math                                                   # basic math operations
import random                                                 # for random numbers
import matplotlib.pyplot as plt                               # for plotting
from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks
from matplotlib.colors import ListedColormap                  # custom color maps
from matplotlib.patches import Patch                          # custom bar chart fills
import matplotlib.ticker as mtick                             # control tick label formatting
import seaborn as sns                                         # for matrix scatter plots
from scipy import stats                                       # summary statistics
import numpy.linalg as linalg                                 # for linear algebra
import scipy.spatial as sp                                    # for fast nearest neighbor search
import scipy.signal as signal                                 # kernel for moving window calculation
from numba import jit                                         # for numerical speed up
from statsmodels.stats.weightstats import DescrStatsW
plt.rc('axes', axisbelow=True)                                # plot all grids below the plot elements
if ignore_warnings == True:                                   
    import warnings
    warnings.filterwarnings('ignore')
cmap = plt.cm.inferno                                         # color map

For the Shapley value approach to feature ranking, we need an additional package and must initialize JavaScript support.

after running this block, you should see a hexagon with the text js, indicating that JavaScript support is ready
if you do not have the shap package installed, uncomment the installation line:

!{sys.executable} -m pip install shap

import sys
#!{sys.executable} -m pip install shap
import shap
shap.initjs()

If you get a package import error, you may first need to install the package. This can usually be accomplished by opening a command window on Windows and typing python -m pip install [package-name]. More assistance is available in the documentation for the respective package.

Design Custom Color Map#

Account for significance by masking nonsignificant values.

this is for demonstration only and could be updated for each plot based on the confidence and uncertainty of the results

my_colormap = plt.cm.get_cmap('RdBu_r', 256)                  # make a custom colormap
newcolors = my_colormap(np.linspace(0, 1, 256))               # define colormap space
white = np.array([250/256, 250/256, 250/256, 1])              # define white color (4 channel)
#newcolors[26:230, :] = white                                 # mask all correlations less than abs(0.8)
#newcolors[56:200, :] = white                                 # mask all correlations less than abs(0.6)
newcolors[76:180, :] = white                                  # mask all correlations less than abs(0.4)
signif = ListedColormap(newcolors)                            # assign as listed colormap
         
my_colormap = plt.cm.get_cmap('inferno', 256)                 # make a custom colormap
newcolors = my_colormap(np.linspace(0, 1, 256))               # define colormap space
white = np.array([250/256, 250/256, 250/256, 1])              # define white color (4 channel)
#newcolors[26:230, :] = white                                 # mask all correlations less than abs(0.8)
newcolors[0:12, :] = white                                    # mask all correlations less than abs(0.6)
#newcolors[86:170, :] = white                                 # mask all correlations less than abs(0.4)
sign1 = ListedColormap(newcolors)                             # assign as listed colormap

Declare Functions#

Here are a few functions to assist with calculating metrics for feature ranking and creating other plots:

plot_corr — plot a correlation matrix
partial_corr — partial correlation coefficient
semipar_corr — semipartial correlation coefficient
mutual_matrix — mutual information matrix containing all pairwise mutual information
mutual_information_objective — my modified version of the MRMR loss function, $(I_{xy} - \mathrm{average}(I_{xx}))$, for feature ranking, using all other predictor features
delta_mutual_information_quotient — change in the mutual information quotient when adding or removing a specific feature, using all other predictor features for comparison
plot_prediction_check - cross validation plot with outlier highlighting and error residual statistics
plot_ranking_summary - custom meat map with feature ranks over various feature ranking methods as a visual summary
weighted_avg_and_std — average and standard deviation accounting for data weights
weighted_percentile — percentile accounting for data weights
histogram_bounds — add confidence intervals to histograms
add_grid — convenience function to add major and minor gridlines to improve plot interpretability

Here are the functions:

def feature_rank_plot(pred,metric,mmin,mmax,nominal,title,ylabel,mask): # feature ranking plot
    mpred = len(pred); mask_low = nominal-mask*(nominal-mmin); mask_high = nominal+mask*(mmax-nominal)
    plt.plot(pred,metric,color='black',zorder=20)
    plt.scatter(pred,metric,marker='o',s=10,color='black',zorder=100)
    plt.plot([-0.5,mpred-0.5],[0.0,0.0],color='black',ls='--',linewidth = 1.0,zorder=1)
    plt.fill_between(np.arange(0,mpred,1),np.zeros(mpred),metric,where=(metric < nominal),interpolate=True,color='dodgerblue',alpha=0.3)
    plt.fill_between(np.arange(0,mpred,1),np.zeros(mpred),metric,where=(metric > nominal),interpolate=True,color='lightcoral',alpha=0.3)
    plt.fill_between(np.arange(0,mpred,1),np.full(mpred,mask_low),metric,where=(metric < mask_low),interpolate=True,color='blue',alpha=0.8,zorder=10)
    plt.fill_between(np.arange(0,mpred,1),np.full(mpred,mask_high),metric,where=(metric > mask_high),interpolate=True,color='red',alpha=0.8,zorder=10)  
    plt.xlabel('Predictor Features'); plt.ylabel(ylabel); plt.title(title)
    plt.ylim(mmin,mmax); plt.xlim([-0.5,mpred-0.5]); add_grid();
    plt.xticks(rotation=270.0)
    return

def plot_corr(corr_matrix,title,limits,mask):                 # plots a graphical correlation matrix 
    my_colormap = plt.cm.get_cmap('RdBu_r', 256)          
    newcolors = my_colormap(np.linspace(0, 1, 256))
    white = np.array([256/256, 256/256, 256/256, 1])
    white_low = int(128 - mask*128); white_high = int(128+mask*128)
    newcolors[white_low:white_high, :] = white                # mask all correlations less than abs(0.8)
    newcmp = ListedColormap(newcolors)
    m = corr_matrix.shape[0]
    im = plt.matshow(corr_matrix,fignum=0,vmin = -1.0*limits, vmax = limits,cmap = newcmp)
    plt.xticks(range(len(corr_matrix.columns)), corr_matrix.columns); ax = plt.gca()
    ax.xaxis.set_label_position('bottom'); ax.xaxis.tick_bottom()
    plt.yticks(range(len(corr_matrix.columns)), corr_matrix.columns)
    plt.colorbar(im, orientation = 'vertical')
    plt.title(title)
    for i in range(0,m):
        plt.plot([i-0.5,i-0.5],[-0.5,m-0.5],color='black')
        plt.plot([-0.5,m-0.5],[i-0.5,i-0.5],color='black')
    plt.ylim([-0.5,m-0.5]); plt.xlim([-0.5,m-0.5])
    plt.xticks(rotation=270.0)

def partial_corr(C):                                          # partial correlation by Fabian Pedregosa-Izquierdo, f@bianp.net
    C = np.asarray(C)
    p = C.shape[1]
    P_corr = np.zeros((p, p), dtype=float)
    for i in range(p):
        P_corr[i, i] = 1
        for j in range(i+1, p):
            idx = np.ones(p, dtype=bool)
            idx[i] = False
            idx[j] = False
            beta_i = linalg.lstsq(C[:, idx], C[:, j])[0]
            beta_j = linalg.lstsq(C[:, idx], C[:, i])[0]
            res_j = C[:, j] - C[:, idx].dot( beta_i)
            res_i = C[:, i] - C[:, idx].dot(beta_j)
            corr = stats.pearsonr(res_i, res_j)[0]
            P_corr[i, j] = corr
            P_corr[j, i] = corr
    return P_corr

def semipartial_corr(C):                                      # Michael Pyrcz modified the function above by Fabian Pedregosa-Izquierdo, f@bianp.net for semipartial correlation
    C = np.asarray(C)
    p = C.shape[1]
    P_corr = np.zeros((p, p), dtype=float)
    for i in range(p):
        P_corr[i, i] = 1
        for j in range(i+1, p):
            idx = np.ones(p, dtype=bool)
            idx[i] = False
            idx[j] = False
            beta_i = linalg.lstsq(C[:, idx], C[:, j])[0]
            res_j = C[:, j] - C[:, idx].dot( beta_i)
            res_i = C[:, i] 
            corr = stats.pearsonr(res_i, res_j)[0]
            P_corr[i, j] = corr
            P_corr[j, i] = corr
    return P_corr

def mutual_matrix(df,features):                               # calculate mutual information matrix
    mutual = np.zeros([len(features),len(features)])
    for i, ifeature in enumerate(features):
        for j, jfeature in enumerate(features):
            if i != j:
                mutual[i,j] = mutual_info_regression(df.iloc[:,i].values.reshape(-1, 1),np.ravel(df.iloc[:,j].values))[0]
    mutual /= np.max(mutual) 
    for i, ifeature in enumerate(features):
        mutual[i,i] = 1.0
    return mutual

def mutual_information_objective(x,y):                        # modified from MRMR loss function, Ixy - average(Ixx)
    mutual_information_quotient = []
    for i, icol in enumerate(x.columns):
        Vx = mutual_info_regression(x.iloc[:,i].values.reshape(-1, 1),np.ravel(y.values.reshape(-1, 1)))
        Ixx_mat = []
        for m, mcol in enumerate(x.columns):
            if i != m:
                Ixx_mat.append(mutual_info_regression(x.iloc[:,m].values.reshape(-1, 1),np.ravel(x.iloc[:,i].values.reshape(-1, 1))))
        Wx = np.average(Ixx_mat)
        mutual_information_quotient.append(Vx/Wx)
    mutual_information_quotient  = np.asarray(mutual_information_quotient).reshape(-1)
    return mutual_information_quotient

def delta_mutual_information_quotient(x,y):                   # standard mutual information quotient
    delta_mutual_information_quotient = []               
    
    Ixy = []
    for m, mcol in enumerate(x.columns):
        Ixy.append(mutual_info_regression(x.iloc[:,m].values.reshape(-1, 1),np.ravel(y.values.reshape(-1, 1))))
    Vs = np.average(Ixy)
    Ixx = []
    for m, mcol in enumerate(x.columns):
        for n, ncol in enumerate(x.columns):
            Ixx.append(mutual_info_regression(x.iloc[:,m].values.reshape(-1, 1),np.ravel(x.iloc[:,n].values.reshape(-1, 1))))
    Ws = np.average(Ixx) 
    
    for i, icol in enumerate(x.columns):          
        Ixy_s = []                                          
        for m, mcol in enumerate(x.columns):
            if m != i:
                Ixy_s.append(mutual_info_regression(x.iloc[:,m].values.reshape(-1, 1),np.ravel(y.values.reshape(-1, 1))))
        Vs_s = np.average(Ixy_s)
        Ixx_s = []
        for m, mcol in enumerate(x.columns):
            if m != i:
                for n, ncol in enumerate(x.columns):
                    if n != i:
                        Ixx_s.append(mutual_info_regression(x.iloc[:,m].values.reshape(-1, 1),np.ravel(x.iloc[:,n].values.reshape(-1, 1))))                  
        Ws_s = np.average(Ixx_s)
        delta_mutual_information_quotient.append((Vs/Ws)-(Vs_s/Ws_s))
    
    delta_mutual_information_quotient  = np.asarray(delta_mutual_information_quotient).reshape(-1)  
    return delta_mutual_information_quotient

def plot_prediction_check(y, y_hat, threshold=0.0, ymin=None, ymax=None,xlabel=None,ylabel=None,title=None): # cross validation plot
    y = np.asarray(y); y_hat = np.asarray(y_hat)

    if ymin is None:
        ymin = min(y.min(), y_hat.min())
    if ymax is None:
        ymax = max(y.max(), y_hat.max())
    if xlabel is None:
        xlabel = 'Actual'
    if ylabel is None:
        ylabel = 'Estimated'

    plt.scatter(y, y_hat,s=40, alpha=0.8,linewidths=0.3,color = 'black', edgecolors="black")
    plt.plot([ymin, ymax], [ymin, ymax],color='black', linestyle='-', linewidth=1.2) # perfect prediction: y_hat = y
    plt.plot([ymin, ymax],[ymin + threshold, ymax + threshold],color='black', linestyle='--', linewidth=0.8) # residual threshold: y_hat = y ± threshold
    plt.plot([ymin, ymax],[ymin - threshold, ymax - threshold],color='black', linestyle='--', linewidth=0.8)

    residual = y_hat - y                                      # identify samples outside the threshold
    residual_mean = np.mean(residual)                         # residual statistics
    residual_sigma = np.std(residual)
    residual_p10 = np.percentile(residual, 10)
    residual_p90 = np.percentile(residual, 90)
    
    outlier_idx = np.where(np.abs(residual) > threshold)[0]

    plt.scatter(y[outlier_idx], y_hat[outlier_idx],s=180,facecolors='white',edgecolors='black',linewidths=1.2,zorder=3) # highlight outliers

    for i in outlier_idx:                                     # add sample index
        plt.text(y[i], y_hat[i], str(i),ha='center', va='center',fontsize=6, zorder=4)

    legend_text = (                                           # residual statistics legend
        'Residual Statistics\n'
        f'Mean = {residual_mean:.1f}\n'
        f'σ = {residual_sigma:.1f}\n'
        f'P10 = {residual_p10:.1f}\n'
        f'P90 = {residual_p90:.1f}'
    )

    plt.text(ymax - 0.03 * (ymax - ymin),ymin + 0.03 * (ymax - ymin),legend_text,ha='right',va='bottom',fontsize=9,
        color='black',zorder=100,bbox=dict(facecolor='white',edgecolor='black',alpha=1.0,boxstyle='round,pad=0.4'))
    plt.xlim(ymin, ymax); plt.ylim(ymin, ymax); plt.xlabel(xlabel); plt.ylabel(ylabel); plt.gca().set_aspect('equal', adjustable='box')

    if title is not None:
        plt.title(title)

def plot_ranking_summary(rankings, features, ranking_names, title='Feature Ranking Summary'):
    rankings = np.asarray(rankings)

    if rankings.ndim != 2:
        raise ValueError("rankings must be a list of 1D ranking arrays.")

    if rankings.shape[1] != len(features):
        raise ValueError("Each ranking array must have the same length as features.")

    if rankings.shape[0] != len(ranking_names):
        raise ValueError("ranking_names must have one name for each ranking array.")

    rank_matrix = np.zeros_like(rankings, dtype=int)          # Convert each ranking metric to rank order, rank 1 = best, larger rank = worse

    for i in range(rankings.shape[0]):
        rank_matrix[i, np.argsort(-rankings[i])] = np.arange(1, rankings.shape[1] + 1)

    mean_rank = np.mean(rank_matrix, axis=0)                  # average rank across all methods for ordering
    indices = np.argsort(mean_rank)
    rank_plot = rank_matrix[:, indices].T
    
    x_edges = np.arange(len(ranking_names) + 1) - 0.5
    y_edges = np.arange(len(features) + 1) - 0.5
    
    ax = plt.gca()
    mesh = ax.pcolormesh(x_edges,y_edges,rank_plot,cmap='inferno_r',shading='flat',edgecolors='black',linewidth=0.5)
    
    ax.set_xticks(range(len(ranking_names))); ax.set_xticklabels(ranking_names, rotation=90)
    ax.set_yticks(range(len(features))); ax.set_yticklabels(np.asarray(features)[indices])
    ax.set_xlim(-0.5, len(ranking_names) - 0.5); ax.set_ylim(len(features) - 0.5, -0.5)
    
    cbar = plt.colorbar(mesh)
    cbar.set_label('Rank (1 = Best)'); cbar.ax.invert_yaxis()
            
    plt.tight_layout()

def weighted_avg_and_std(values, weights):                    # calculate weighted statistics (Eric O Lebigot, stack overflow)
    average = np.average(values, weights=weights)
    variance = np.average((values-average)**2, weights=weights)
    return (average, math.sqrt(variance))

def weighted_percentile(data, weights, perc):                 # calculate weighted percentile (iambr on StackOverflow @ https://stackoverflow.com/questions/21844024/weighted-percentile-using-numpy/32216049) 
    ix = np.argsort(data)
    data = data[ix] 
    weights = weights[ix] 
    cdf = (np.cumsum(weights) - 0.5 * weights) / np.sum(weights) 
    return np.interp(perc, cdf, data)

def histogram_bounds(values,weights,color):                   # add uncertainty bounds to a histogram          
    p10 = weighted_percentile(values,weights,0.1); avg = np.average(values,weights=weights); p90 = weighted_percentile(values,weights,0.9)
    plt.plot([p10,p10],[0.0,45],color = color,linestyle='dashed')
    plt.plot([avg,avg],[0.0,45],color = color)
    plt.plot([p90,p90],[0.0,45],color = color,linestyle='dashed')

def add_grid():                                               # add major and minor gridlines
    plt.gca().grid(True, which='major',linewidth = 1.0); plt.gca().grid(True, which='minor',linewidth = 0.2) # add y grids
    plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)
    plt.gca().xaxis.set_minor_locator(AutoMinorLocator()); plt.gca().yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks

Set the Working Directory#

I always like to do this because:

I don’t lose track of files
it simplifies subsequent reads and writes by avoiding the need to include the full file path each time

When I work on projects, I like to use separate folders for:

inputs
intermediate files
result tables
result plots
documentation and papers

To remain organized, I find it helpful to set and change the working directory throughout my workflows.

#os.chdir("d:/PGE383")                                   # set the working directory

Loading Tabular Data#

Here’s the command to load our comma-delimited data file into a Pandas DataFrame object.

Let’s load the provided multivariate, spatial dataset, unconv_MV_v4.csv, from my GeoDataSets repository.

This dataset contains variables from 1,000 unconventional wells, including:

Feature	Description	Units
`Por`	Average well porosity	$\%$
`Perm`	Log-transformed permeability to linearize relationships with other variables	$\log_{10}(\mathrm{mD})$
`AI`	Acoustic impedance	$\mathrm{kg/m^3 \cdot m/s \cdot 10^6}$
`Brittle`	Brittleness ratio	$\%$
`TOC`	Total organic carbon	$\%$
`VR`	Vitrinite reflectance	$\%$
`Prod`	Initial production, 90-day average	$\mathrm{MCFPD}$

Note that the dataset is synthetic, so it can be freely used with citation to the source:

Pyrcz, Michael J. (2021). GeoDataSets: Synthetic Subsurface Data Repository (0.0.1). Zenodo. https://doi.org/10.5281/zenodo.5564874

We load the data with the Pandas read_csv function into a DataFrame called my_data and then preview it to make sure it loaded correctly.

Note, you will have to update the path in quotes below to your own working directory. The path format is different on a Mac; for example, you could use "~/PGE".

# Configure SSL certificate verification for GitHub data loading
import ssl
import certifi
import urllib.request

ssl._create_default_https_context = lambda: ssl.create_default_context(
    cafile=certifi.where()
)

idata = 0
if idata == 0:
    df = pd.read_csv('https://raw.githubusercontent.com/GeostatsGuy/GeoDataSets/master/unconv_MV_v4.csv') # load data from Dr. Pyrcz's GitHub repository  
    
    response = 'Prod'                                         # specify the response feature
    x = df.copy(deep = True); x = x.drop(['Well',response],axis='columns') # make predictor and response DataFrames
    Y = df.loc[:,response]

    features = x.columns.values.tolist() + [Y.name]           # store the names of the features
    pred = x.columns.values.tolist()
    resp = Y.name
    
    xmin = [6.0,0.0,1.0,10.0,0.0,0.9]; xmax = [24.0,10.0,5.0,85.0,2.2,2.9] # set the minimum and maximum values for plotting
    Ymin = 500.0; Ymax = 9000.0
    
    predlabel = ['Porosity (%)','Permeability (mD)','Acoustic Impedance (kg/m2s*10^6)','Brittleness Ratio (%)', # set the names for plotting
                 'Total Organic Carbon (%)','Vitrinite Reflectance (%)']
    resplabel = 'Normalized Initial Production (MCFPD)'
    
    predtitle = ['Porosity','Permeability','Acoustic Impedance','Brittleness Ratio', # set the units for plotting
                 'Total Organic Carbon','Vitrinite Reflectance']
    resptitle = 'Normalized Initial Production'

    featurelabel = predlabel + [resplabel]                    # make feature labels and titles for concise code
    featuretitle = predtitle + [resptitle]
    
    m = len(pred) + 1
    mpred = len(pred)

# elif idata == 1:
#     names = {'Porosity':'Por'}
    
#     df = pd.read_csv('https://raw.githubusercontent.com/GeostatsGuy/GeoDataSets/master/12_sample_data.csv') # load data from Dr. Pyrcz's GitHub repository  
#     df = df.rename(columns=names)
#     df['Por'] = df['Por'] * 100.0; df['AI'] = df['AI'] / 1000.0; 
#     df.drop('Unnamed: 0',axis=1,inplace=True) 
    
#     features = df.columns.values.tolist()                          # store the names of the features

#     xmin = [0.0,0.0,0.0,4.0,0.0,6.5,1.4,1600.0,10.0,1300.0,1.6]; xmax = [10000.0,10000.0,1.0,19.0,500.0,8.3,3.6,6200.0,50.0,2000.0,12.0] # set the minimum and maximum values for plotting
    
#     flabel = ['Well (ID)','X (m)','Y (m)','Depth (m)','Porosity (fraction)','Permeability (mD)','Acoustic Impedance (kg/m2s*10^6)','Facies (categorical)',
#               'Density (g/cm^3)','Compressible velocity (m/s)','Youngs modulus (GPa)', 'Shear velocity (m/s)', 'Shear modulus (GPa)'] # set the names for plotting

#     ftitle = ['Well','X','Y','Depth','Porosity','Permeability','Acoustic Impedance','Facies',
#               'Density','Compressible velocity','Youngs modulus', 'Shear velocity', 'Shear modulus']

elif idata == 2:  
    df = pd.read_csv('https://raw.githubusercontent.com/GeostatsGuy/GeoDataSets/master/res21_2D_wells.csv') # load data from Dr. Pyrcz's GitHub repository  
    
    response = 'CumulativeOil'                                    # specify the response feature
    x = df.copy(deep = True); x = x.drop(['Well_ID','X','Y',response],axis='columns') # make predictor and response DataFrames
    Y = df.loc[:,response]
    
    features = x.columns.values.tolist() + [Y.name]               # store the names of the features
    pred = x.columns.values.tolist()
    resp = Y.name

    xmin = [1.0,0.0,0.0,4.0,0.0,6.5,1.4,1600.0,10.0,1300.0,1.6]; xmax = [75.0,10000.0,10000.0,19.0,500.0,8.3,3.6,6200.0,50.0,2000.0,12.0] # set the minimum and maximum values for plotting
    Ymin = 0.0; Ymax = 3000.0
    
    predlabel = ['Well (ID)','X (m)','Y (m)','Porosity (fraction)','Permeability (mD)','Acoustic Impedance (kg/m2s*10^6)',
              'Density (g/cm^3)','Compressible velocity (m/s)','Youngs modulus (GPa)', 'Shear velocity (m/s)', 'Shear modulus (GPa)'] 
    resplabel = 'Cumulative Production (MSTB)'
    
    predtitle = ['Well','X','Y','Porosity','Permeability','Acoustic Impedance',
              'Density (g/cm^3)','Compressible velocity','Youngs modulus', 'Shear velocity', 'Shear modulus'] 
    resptitle = 'Cumulative Production'

    featurelabel = predlabel + [resplabel]                    # make feature labels and titles for concise code
    featuretitle = predtitle + [resptitle]
    
    m = len(pred) + 1
    mpred = len(pred)

Above we specify the ranges for each feature as lists containing:

the minimum and maximum values
we could also calculate the feature ranges directly from the data with code like this:

Pormin = np.min(df['Por'].values)                             # extract ndarray of data table column
Pormax = np.max(df['Por'].values)                             # and calculate min and max

However, this would not necessarily result in easy-to-understand color bars and axis scales.

it is better to carefully select the ranges to maximize plot and gridline clarity
feature labels for the axes are included as lists for ease of plotting

Visualize the DataFrame#

Visualizing the DataFrame is a useful first check of the data.

many things can go wrong, e.g., we loaded the wrong data, all the features did not load, etc.

We can preview the data by utilizing the pandas head method of the DataFrame class.

provides a fast and clean format for checking the data
add the parameter n=13 to see the first 13 rows of the dataset

df.head(n=13)                                                 # we could also use this command for a table preview 

	Well	Por	Perm	AI	Brittle	TOC	VR	Prod
0	1	12.08	2.92	2.80	81.40	1.16	2.31	1695.360819
1	2	12.38	3.53	3.22	46.17	0.89	1.88	3007.096063
2	3	14.02	2.59	4.01	72.80	0.89	2.72	2531.938259
3	4	17.67	6.75	2.63	39.81	1.08	1.88	5288.514854
4	5	17.52	4.57	3.18	10.94	1.51	1.90	2859.469624
5	6	14.53	4.81	2.69	53.60	0.94	1.67	4017.374438
6	7	13.49	3.60	2.93	63.71	0.80	1.85	2952.812773
7	8	11.58	3.03	3.25	53.00	0.69	1.93	2670.933846
8	9	12.52	2.72	2.43	65.77	0.95	1.98	2474.048178
9	10	13.25	3.94	3.71	66.20	1.14	2.65	2722.893266
10	11	15.04	4.39	2.22	61.11	1.08	1.77	3828.247174
11	12	16.19	6.30	2.29	49.10	1.53	1.86	5095.810104
12	13	16.82	5.42	2.80	66.65	1.17	1.98	4091.637316

Summary Statistics for Tabular Data#

There are many efficient methods to calculate summary statistics from tabular data in DataFrames.

the pandas describe method of the DataFrame class provides the count, mean, standard deviation, minimum, maximum, and quartiles in a convenient data table
use transpose to flip the table so that features are on the rows and the statistics are on the columns
the reported percentiles may be modified by setting the percentiles parameter. For example, to report the commonly used P10, P50, and P90:

df.describe(percentiles=[0.1,0.5,0.9]).transpose()

df.describe().transpose()                                     # calculate summary statistics for the data

	count	mean	std	min	25%	50%	75%	max
Well	200.0	100.500000	57.879185	1.000000	50.750000	100.500000	150.250000	200.000000
Por	200.0	14.991150	2.971176	6.550000	12.912500	15.070000	17.402500	23.550000
Perm	200.0	4.330750	1.731014	1.130000	3.122500	4.035000	5.287500	9.870000
AI	200.0	2.968850	0.566885	1.280000	2.547500	2.955000	3.345000	4.630000
Brittle	200.0	48.161950	14.129455	10.940000	37.755000	49.510000	58.262500	84.330000
TOC	200.0	0.990450	0.481588	-0.190000	0.617500	1.030000	1.350000	2.180000
VR	200.0	1.964300	0.300827	0.930000	1.770000	1.960000	2.142500	2.870000
Prod	200.0	3864.407081	1553.277558	839.822063	2686.227611	3604.303506	4752.637555	8590.384044

Ranking features is really an effort to understand the features and their relationships with each other. We will start with basic data visualization and move to more complicated methods such as partial correlation and recursive feature elimination.

Coverage#

Let’s start with the concept of feature coverage.

if a feature is available for only a small proportion of the samples, then we may not want to include it, as it will result in more issues with feature imputation, i.e., estimation of missing data
by removing a couple of features with poor coverage, we may improve our model because there are limitations with feature imputation. Feature imputation can introduce bias into statistics and additional error into our prediction models
if listwise deletion is applied to deal with missing values, features with low coverage can result in a lot of removed data!

Let’s start with a data completeness bar chart showing the proportion of missing records.

plt.subplot(111)
(df.isnull().sum()/len(df)).plot(kind = 'bar')                # calculate DataFrame with percentage missing by feature
plt.xlabel('Feature'); plt.ylabel('Percentage of Missing Values'); plt.title('Data Completeness'); plt.ylim([0.0,1.0])
plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=0.8, wspace=0.2, hspace=0.2); add_grid(); plt.show()

_images/aee47f1261ac842a4b27dd2dda43319d24880ba6a87fd412ec1aa3a8c3044301.png

For the provided example dataset, the plot should show zero missing values. There are no missing data, so the Proportion of Missing Records is 0.0 for all features.

If you wanted to test this plot with some missing data, run this code first:

proportion_NaN = 0.1                                    # proportion of values in DataFrame to remove

remove = np.random.random(df.shape) < proportion_NaN    # make the boolean array for removal
print('Fraction of removed values in mask ndarray = ' + str(round(remove.sum()/remove.size,3)) + '.')

df_mask = df.mask(remove)                               # make a new DataFrame with specified proportion removed

Then remove the above code and reload the data to continue and obtain consistent results with the discussions below.

Warning: the data completeness plot does not tell the whole story. For example, if 20% of feature A is missing and 20% of feature B is missing, are these missing values in the same or different samples? This can have a large impact if you perform listwise deletion.

if there is not too much data, then we can actually visualize data coverage across all samples and features in a boolean table or heatmap (see below)
this method may identify specific samples with many missing features that could be removed to improve overall coverage, or other trends or structures in the missing data that may result in sampling bias

df_temp = df.copy(deep=True)                                  # make a deep copy of the DataFrame
df_bool = df_temp.isnull()                                    # true is value, false if NaN
#df_bool = df_bool.set_index(df_temp.pop('UWI'))              # set the index / feature for the heat map y column
heat = sns.heatmap(df_bool, cmap=['r','w'], annot=False, fmt='.0f',cbar=False,linecolor='black',linewidth=0.1) # make the binary heat map, no bins
heat.set_xticklabels(heat.get_xticklabels(), rotation=90, fontsize=8)
heat.set_yticklabels(heat.get_yticklabels(), rotation=0, fontsize=8)
heat.set_title('Data Completeness Heatmap',fontsize=16); heat.set_xlabel('Feature',fontsize=12); heat.set_ylabel('Sample (Index)',fontsize=12)
plt.subplots_adjust(left=0.0, bottom=0.0, right=1.8, top=1.6, wspace=0.2, hspace=0.2); plt.show()

_images/7ef6b30311ccd1ae428f2984c9bdd9999ef6f0a3557694854edd67192b8859a0.png

Once again, this plot should be quite boring for the provided dataset with perfect coverage; every cell should be filled in red.

add the code above to remove some records and test this plot. White cells represent missing records.

Feature Imputation#

See the chapter on feature imputation to learn what to do about missing data.

For now, a concise treatment here: we will only apply listwise deletion and move on.

we remove all records with any missing feature values. While this is quite simple, it is a sledgehammer approach to ensure the perfect coverage required by the feature ranking methods that we are about to demonstrate.

df.dropna(axis=0,how='any',inplace=True)                      # likewise deletion

Summary Statistics#

In any multivariate work, we should start with the univariate analysis: summary statistics of one variable at a time.

The summary statistics ranking method is qualitative; we are asking:

Are there data issues?
Do we trust the features? Do we trust all the features equally?
Are there issues that need to be addressed before we develop any multivariate workflows?

There are many efficient methods to calculate summary statistics from tabular data in DataFrames, including:

the describe() method provides the count, mean, standard deviation, minimum, maximum, and quartiles in a compact data table. We use the transpose() method to flip the table so that features are on the rows and the statistics are on the columns.

df.describe().transpose()                                     # DataFrame summary statistics

	count	mean	std	min	25%	50%	75%	max
Well	200.0	100.500000	57.879185	1.000000	50.750000	100.500000	150.250000	200.000000
Por	200.0	14.991150	2.971176	6.550000	12.912500	15.070000	17.402500	23.550000
Perm	200.0	4.330750	1.731014	1.130000	3.122500	4.035000	5.287500	9.870000
AI	200.0	2.968850	0.566885	1.280000	2.547500	2.955000	3.345000	4.630000
Brittle	200.0	48.161950	14.129455	10.940000	37.755000	49.510000	58.262500	84.330000
TOC	200.0	0.990450	0.481588	-0.190000	0.617500	1.030000	1.350000	2.180000
VR	200.0	1.964300	0.300827	0.930000	1.770000	1.960000	2.142500	2.870000
Prod	200.0	3864.407081	1553.277558	839.822063	2686.227611	3604.303506	4752.637555	8590.384044

Summary statistics are a critical first step in data checking.

this includes the number of valid (non-null) values for each feature. The count statistic reports the number of non-null observations for each variable and excludes np.NaN values.
we can see the general behavior of each feature, including central tendency, such as the mean, and dispersion, such as the standard deviation.
we can identify issues with negative values, extreme values, and values that are outside the range of plausible values for each property.
the data looks to be in pretty good shape, and for brevity, we will skip outlier detection. Let’s look at the univariate distributions.

Univariate Distributions#

As with summary statistics, this feature ranking method is:

a qualitative check for issues with the data and to assess our confidence in each feature.

It is better not to include a feature with low confidence in its quality, as it may be misleading while adding to model complexity, as discussed previously.

nbins = 20                                                    # number of histogram bins
for i, feature in enumerate(features):                        # plot histograms with central tendency and P10 and P90 labeled
    plt.subplot(4,2,i+1)
    y,_,_ = plt.hist(x=df[feature],weights=None,bins=nbins,alpha = 0.8,edgecolor='black',color='darkorange',density=True)
    plt.xlabel(feature); plt.ylabel('Frequency'); plt.ylim([0.0,y.max()*1.10]); plt.title(featuretitle[i]); add_grid() 

plt.subplots_adjust(left=0.0, bottom=0.0, right=2., top=4.1, wspace=0.15, hspace=0.3); plt.show()

_images/1e269b4ec33bbd75745fb58f1c3a85376aa1b955f9315ecab362983711a13ae2.png

The univariate distributions look good:

there are no obvious outliers.
the permeability is positively skewed, as is often observed.
the corrected TOC has a small spike, but it appears reasonable.

Bivariate Distributions#

Matrix scatter plots are a very efficient way to observe the full set of combinatorial bivariate relationships between the features.

this provides another opportunity to use data visualization to identify potential data issues.

We can assess the presence of:

collinearity — linear relationships between pairs of features.
constraints — regions of the cross plot with few or no samples.
nonlinearity — nonlinear relationships between features.
heteroscedasticity — a relationship between the variance and the mean, or a change in conditional variance with conditional mean.

pairgrid = sns.PairGrid(df)                                   # matrix scatter plots
pairgrid = pairgrid.map_upper(plt.scatter, color = 'darkorange', edgecolor = 'black', alpha = 0.8, s = 10)
pairgrid = pairgrid.map_diag(plt.hist, bins = 20, color = 'darkorange',alpha = 0.8, edgecolor = 'k') # map a density plot to the lower triangle
pairgrid = pairgrid.map_lower(sns.kdeplot, cmap = plt.cm.inferno, 
                              shade = False, shade_lowest = False, alpha = 1.0, n_levels = 10)
pairgrid.add_legend()
plt.subplots_adjust(left=0.0, bottom=0.0, right=0.9, top=0.9, wspace=0.2, hspace=0.2); plt.show()

_images/edcc4a5fd6dd90cf3d3aa3fa9d635952e7cfa3feab7a30fc97df684fcd801ab3.png

This matrix scatter plot communicates a lot of information. How could we use this plot for feature ranking?

we can identify features that are closely related to each other. For example, if two features have an almost perfect monotonic, linear, or near-linear relationship, we should consider removing one. This is a simple case of collinearity that can result in model instability, as discussed above.
we can check for linear versus nonlinear relationships. If we observe nonlinear bivariate relationships, this may impact the choice of feature-ranking methods or the quality of results from methods that assume linear relationships.
we can identify constraints, multiple modes, and heteroscedasticity between variables. These patterns may indicate the presence of multiple populations that should be analyzed and modeled separately.

Warning: Bivariate visualization and analysis are not sufficient to understand all of the multivariate relationships in the data.

bivariate visualization involves a substantial collapse of information through marginalization. Relationships that depend on three or more features may not be apparent in any individual bivariate plot.
multicollinearity includes strong linear relationships involving two or more predictor features. Some forms of multicollinearity can be difficult to identify from bivariate plots alone.

Covariance#

Covariance provides a measure of the strength and direction of the linear relationship between each predictor feature and the response feature.

Now we specify that the goal of this study is to predict production, our response variable, from the other available predictor features. We are thinking predictively, not inferentially; we want to estimate the function, $\hat{f}$, to accomplish this:

\[ Y = \hat{f}(X_1,\ldots,X_n) \]

where $Y$ is our response feature and $X_1,\ldots,X_n$ are our predictor features. If we retained all of our predictor features to predict the response, we would have:

\[ Prod = \hat{f}(Por,Perm,AI,Brittle,TOC,VR) \]

Now back to covariance. The covariance is defined as:

\[ C_{xy} = \frac{\sum_{i=1}^{n} (x_i - \overline{x})(y_i - \overline{y})}{(n-1)} \]

Important points about covariance:

measures the strength and direction of the linear relationship between two variables.
is sensitive to the dispersion / variance of both the predictor and response.
is scale-dependent, so the magnitude of covariance depends on the units and dispersion of the predictor and response.

We can use the following command to build a covariance matrix:

df.iloc[:,1:8].cov()                                    # covariance matrix sliced predictors vs. response

The output is a new Pandas DataFrame, so we can slice the last column to get a Pandas Series (ndarray with names) containing the covariances between all predictor features and the response.

covariance = df.iloc[:,df.columns.get_indexer(features)].cov().iloc[len(features)-1,:len(features)] # calculate covariance matrix and slice for only pred - resp
cov_matrix = df.iloc[:,df.columns.get_indexer(features)].cov()
plt.subplot(121)
plot_corr(cov_matrix,'Covariance Matrix',4000.0,0.1)          # using our correlation matrix visualization function
plt.xlabel('Features'); plt.ylabel('Features') 

plt.subplot(122)
feature_rank_plot(features[:-1],covariance[:-1],-10000.0,10000.0,0.0,'Feature Ranking, Covariance with ' + resp,'Covariance',0.1)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=0.8, wspace=0.2, hspace=0.3); plt.show()

_images/e1c88bdca0409c9c9b79a21bb648a24636a4bcf45fd8bbb479d3c3b21bd8e1b1.png

The covariance is not very useful for feature ranking because its magnitude depends on the variance of each feature, and the variance depends on the units used to measure the feature.

For example:

what is the variance of porosity when it is measured as a fraction versus a percentage?
what is the variance of permeability when it is measured in Darcy versus milliDarcy?

We can show that if we apply a constant multiplier, $c$, to a feature, $X$, the variance changes according to this relationship. The proof follows from the expectation formulation of variance:

\[ \sigma_{cX}^2 = c^2 \cdot \sigma_X^2 \]

For example, converting porosity from a fraction to a percentage multiplies the feature by 100. Therefore, the variance increases by a factor of $100^2 = 10{,}000$.

Conversely, moving from percentage to fraction decreases the variance of porosity by a factor of 10,000.

Therefore:

the variance of each feature is potentially arbitrary based on the units used to measure it, except when all features are measured in the same units.

Warning: covariance magnitudes are consequently not directly comparable across features with different units or scales.

Correlation coefficients are standardized covariances; therefore, they remove this arbitrary magnitude issue and provide a common scale for comparing the strength of linear relationships.

Correlation Coefficient#

The correlation coefficient, also known as the Pearson product moment correlation coefficient, provides a measure of the strength and direction of the linear relationship between each predictor feature and the response feature.

\[ \rho_{xy} = \frac{\sum_{i=1}^{n} (x_i - \overline{x})(y_i - \overline{y})}{(n-1)\sigma_x \sigma_y}, \, -1.0 \le \rho_{xy} \le 1.0 \]

The correlation coefficient:

measures the strength and direction of the linear relationship between two variables.
removes the sensitivity to the dispersion / variance of both the predictor and response feature by normalizing the covariance by the product of their standard deviations.
provides a standardized measure, from $-1$ to $1$, that allows the linear relationships between features to be compared regardless of their units.

We can use the following command to calculate a correlation matrix with the pandas corr method of the DataFrame class:

df.iloc[:,1:8].corr()

The output is a new Pandas DataFrame, so we can slice the last column to get a Pandas Series (ndarray with names) containing the correlations between all predictor features and the response.

correlation = df.iloc[:,df.columns.get_indexer(features)].corr().iloc[len(features)-1,:len(features)] # calculate covariance matrix and slice for only pred - resp
corr_matrix = df.iloc[:,df.columns.get_indexer(features)].corr()

plt.subplot(121)
plot_corr(corr_matrix,'Correlation Matrix',1.0,0.5)           # using our correlation matrix visualization function
plt.xlabel('Features'); plt.ylabel('Features')

plt.subplot(122)
feature_rank_plot(features[:-1],correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Correlation with ' + resp,'Correlation',0.5)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=0.8, wspace=0.2, hspace=0.3); plt.show()

_images/2f3d46d48328553f2918b2dd58bb21a987c70d35550eaff3e6b9a626eb4376ec.png

From the correlation matrix heat map and predictor feature to response feature correlation line plot, we can observe:

porosity, permeability, and total organic carbon have the strongest linear relationships with production
acoustic impedance has a weak negative relationship with production
brittleness’ correlation with production is very close to 0.0, likely due to the nonlinear relationship. There is a brittleness ratio sweet spot for production (rock that is not too soft nor too hard)!

We could also look at the full correlation matrix to evaluate the potential for redundancy between predictor features.

strong correlation between porosity and permeability, and between porosity and TOC
strong negative correlation between TOC and acoustic impedance

We are still limited to a strict linear relationship. The rank correlation allows us to relax this assumption.

Rank Correlation Coefficient#

The rank correlation coefficient, also known as the Spearman correlation coefficient, applies the rank transform to the data prior to calculating the correlation coefficient. To calculate the rank transform, simply replace the data values with their ranks, $R_x = 1,\dots,n$, where $n$ is the number of observations, with $1$ assigned to the minimum value and $n$ assigned to the maximum value.

\[ \rho_{R_x R_y} = \frac{\sum_{i=1}^{n} (R_{x_i} - \overline{R_x})(R_{y_i} - \overline{R_y})}{(n-1)\sigma_{R_x} \sigma_{R_y}}, \, -1.0 \le \rho_{xy} \le 1.0 \]

For example, the rank ordering of the $x$ data can be represented as:

\[ x_\alpha, \, \forall \alpha = 1,\dots, n, \, | \, x_i \ge x_j \, \forall \, i \gt j \]

and the corresponding ranks are:

\[ R_{x_i} = i \]

The rank correlation:

measures a monotonic relationship, relaxing the assumption of linearity
removes the sensitivity to the dispersion / variance of the predictor and response by replacing the original data values with their ranks

We can use the following command to build a rank correlation matrix and calculate the p-value:

stats.spearmanr(df.iloc[:,1:8])

The output includes a correlation matrix and a p-value matrix. We can slice the last column of the correlation matrix to get a Pandas Series (ndarray with names) containing the correlations between all predictor features and the response.

We also get a very convenient pval 2D ndarray containing the two-sided (two-tail, summing the probabilities over both tails) p-value for a hypothesis test with:

\[ H_o: \rho_{R_x, R_y} = 0 \]

\[ H_1: \rho_{R_x, R_y} \ne 0 \]

Recall: if the p-value is less than the assigned $\alpha$-value, we reject the null hypothesis that the correlation is zero; i.e., we can say that the rank correlation is statistically significant.

Let’s keep the p-values between all the predictor features and our response feature.

rank_correlation, rank_correlation_pval = stats.spearmanr(df.iloc[:,df.columns.get_indexer(features)]) # calculate the rank correlation coefficient
rank_matrix = pd.DataFrame(rank_correlation,columns=corr_matrix.columns)
rank_correlation = rank_correlation[:,len(features)-1][:len(features)]

print("\nRank Correlation Coefficient p-values:")
print("".ljust(15) + "".join(f"{feature:>12}" for feature in features))
for i, feature in reversed(list(enumerate(features))):
    print(f"{feature:<15}" + "".join(f"{rank_correlation_pval[i,j]:12.5f}" for j in range(len(features))))
print("\n")

plt.subplot(121)
plot_corr(rank_matrix,'Rank Correlation Matrix',1.0,0.5)      # using our correlation matrix visualization function
plt.xlabel('Features'); plt.ylabel('Features')

plt.subplot(122)
feature_rank_plot(features[:-1],rank_correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Rank Correlation with ' + resp,'Rank Correlation',0.5)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=0.8, wspace=0.2, hspace=0.3); plt.show()

Rank Correlation Coefficient p-values:
                        Por        Perm          AI     Brittle         TOC          VR        Prod
Prod                0.00000     0.00000     0.00000     0.41461     0.00000     0.00033     0.00000
VR                  0.02433     0.13414     0.00000     0.00027     0.00000     0.00000     0.00033
TOC                 0.00000     0.00000     0.00000     0.00291     0.00000     0.00000     0.00000
Brittle             0.00286     0.19281     0.17951     0.00000     0.00291     0.00027     0.41461
AI                  0.00000     0.00037     0.00000     0.17951     0.00000     0.00000     0.00000
Perm                0.00000     0.00000     0.00037     0.19281     0.00000     0.13414     0.00000
Por                 0.00000     0.00000     0.00000     0.00286     0.00000     0.02433     0.00000

_images/e29d83946a542c70e4a67ad7babef63a6e4cac82cc31a8594856c51800b1a265.png

From the rank correlation matrix heat map and the predictor feature-to-response feature rank correlation line plot, we can observe:

that the rank correlation coefficients are similar to the correlation coefficients, suggesting that nonlinearity and outliers are not substantially impacting the correlation-based feature ranking.

With regard to the rank correlation p-values:

at a typical $\alpha$-value of 0.05, only the rank correlation between brittleness and production fails to reject the null hypothesis; therefore, it is not significantly different from 0.0.

It is useful to look at the difference between the correlation coefficient and rank correlation coefficient.

plt.subplot(121)                                              # plot correlation matrix with significance colormap
diff = corr_matrix.values - rank_matrix.values
diff_matrix = pd.DataFrame(diff,columns=corr_matrix.columns)
plot_corr(diff_matrix,'Correlation - Rank Correlation',0.1,0.1) # using our correlation matrix visualization function
plt.xlabel('Features'); plt.ylabel('Features')

corr_diff = correlation - rank_correlation

plt.subplot(122)
feature_rank_plot(features[:-1],corr_diff[:-1],-0.10,0.10,0.0,'Correlation Coefficient - Rank Correlation Coefficient','Correlation Diffference',0.1)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=0.8, wspace=0.2, hspace=0.3); plt.show()

_images/3cd6fe8aba079ab6d451cf0340bcc46968b3bca396869b6b2c9b79cc91e77205.png

From the correlation and rank correlation difference matrix heat map and the predictor feature-to-response feature correlation and rank correlation difference line plot, we can observe:

the correlations of porosity and vitrinite reflectance with production increase with the rank transform, which may be caused by a reduced impact of nonlinearity and outliers
the correlation of brittleness with production decreases when we reduce the impact of nonlinearity and outliers

All of these methods up to now have considered one feature at a time. We can also consider methods that consider all features jointly to ‘isolate’ the relationship of each predictor feature with the response feature.

Warning — covariance, correlation coefficient, and rank correlation coefficient only account for feature relevancy and not feature redundancy.

Partial Correlation Coefficient#

This is a linear correlation coefficient that controls for the effects of all the remaining variables.

$\rho_{XY.Z}$ and $\rho_{YX.Z}$ are the partial correlations between $X$ and $Y$, and $Y$ and $X$, respectively, after controlling for $Z$ feature(s).

To calculate the partial correlation coefficient between $X$ and $Y$ given $Z_i, \forall \quad i = 1,\ldots,m-2$, the remaining features, we use the following steps:

perform linear, least-squares regression to predict $X$ from $Z_i, \forall \quad i = 1,\ldots,m-2$. $X$ is regressed on the predictors to calculate the estimate, $X^*$.
calculate the residuals in Step #1, $X-X^*$, where $X^* = f(Z_{1,\ldots,m-2})$ is the linear regression model estimate.
perform linear, least-squares regression to predict $Y$ from $Z_i, \forall \quad i = 1,\ldots,m-2$. $Y$ is regressed on the predictors to calculate the estimate, $Y^*$.
calculate the residuals in Step #3, $Y-Y^*$, where $Y^* = f(Z_{1,\ldots,m-2})$ is the linear regression model estimate.
calculate the correlation coefficient between the residuals from Steps #2 and #4, $\rho_{X-X^*,Y-Y^*}$.

The partial correlation provides a measure of the linear relationship between $X$ and $Y$ while controlling for the effects of the other features on both $X$ and $Y$.

We use the function, partial_corr, declared previously, taken from,

Fabian Pedregosa-Izquierdo, f@bianp.net. The original code is on GitHub.

To use this method, we must specify:

two variables to compare, $X$ and $Y$
other variables to control, $Z_{1,\ldots,m-2}$
linear relationships between the variables
no significant outliers

For interpretation and statistical inference, it is also useful to consider whether the variables are approximately bivariate normal.

Warning: Partial correlation only accounts for linear relationships for relevancy and redundancy. Strong nonlinear relationships may therefore be missed. The interpretation of the correlation coefficient can also be affected by outliers and departures from the assumptions used for statistical inference.

From inspection of the matrix scatter plot,

the data looks somewhat ok
but we have some departures from bivariate normality.

We could consider Gaussian univariate transforms to improve this. This option is provided later.

partial_correlation = partial_corr(df.iloc[:,df.columns.get_indexer(features)]) # calculate the partial correlation coefficients
partial_matrix = pd.DataFrame(partial_correlation,columns=corr_matrix.columns)
partial_correlation = partial_correlation[:,len(features)-1][:len(features)] # extract a single row and remove production with itself        

plt.subplot(121)
plot_corr(partial_matrix,'Partial Correlation Matrix',1.0,0.5) # using our correlation matrix visualization function
plt.xlabel('Features'); plt.ylabel('Features')

plt.subplot(122)
feature_rank_plot(features[:-1],partial_correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Partial Correlation with ' + resp,'Partial Correlation',0.5)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=0.8, wspace=0.2, hspace=0.3); plt.show()

_images/acd3ab801095340914c7d41b39381ba574719a483edeb9d031ba3a61865b6632.png

Now we see a lot of new information about the unique contributions of each predictor feature! We can observe,

porosity and permeability are strongly correlated with each other, so their partial correlations are reduced substantially
the absolute correlations of acoustic impedance and vitrinite reflectance with production increase, reflecting their unique linear associations with production after controlling for the other predictor features
total organic carbon flipped signs! When we control for all other variables, it has a negative relationship with production.

With the partial correlation coefficients, we have controlled for the influence of all other predictor features on both the specific predictor feature and the response feature. The semipartial correlation filters out the influence of all other predictor features on the specific predictor feature, while retaining the raw response variable.

Semipartial Correlation Coefficient#

The semipartial correlation coefficient controls for the effects of all the remaining features, $Z$, on $X$, and then calculates the correlation between the residual $X-X^*$ and $Y$.

note: we do not control for the influence of the $Z$ features on the response feature, $Y$, as we did with the partial correlation coefficient.

To calculate the semipartial correlation coefficient between $X$ and $Y$ given $Z_i, \forall \quad i = 1,\ldots,m-2$, the remaining features, we use the following steps:

perform linear, least-squares regression to predict $X$ from $Z_i, \forall \quad i = 1,\ldots,m-2$. $X$ is regressed on the remaining predictor features to calculate the estimate, $X^*$
calculate the residuals in Step #1, $X-X^*$, where $X^* = f(Z_{1,\ldots,m-2})$ is the linear regression model estimate
calculate the correlation coefficient between the residuals from Step #2, $X-X^*$, and the $Y$ response feature, $\rho_{X-X^*,Y}$

The semipartial correlation coefficient provides a measure of the linear relationship between $X$ and $Y$ while controlling for the effect of the other $Z$ predictor features on the predictor feature, $X$, to get,

the unique linear association of $X$ with $Y$.

Warning: Semipartial correlation only accounts for linear relationships for relevancy and redundancy. Strong nonlinear relationships may therefore be missed. The interpretation of the correlation coefficient can also be affected by outliers and departures from the assumptions used for statistical inference.

Below we use the function, semipartial_corr, modified from spartial_corr, taken from,

Fabian Pedregosa-Izquierdo, f@bianp.net. The original code is on GitHub.

semipartial_correlation = semipartial_corr(df.iloc[:,df.columns.get_indexer(features)]) # calculate the semi-partial correlation coefficients
semipartial_matrix = pd.DataFrame(semipartial_correlation,columns=corr_matrix.columns)
semipartial_correlation = semipartial_correlation[:,len(features)-1][:len(features)] # extract a single row and remove production with itself

plt.subplot(121)
plot_corr(semipartial_matrix,'Semi-partial Correlation Matrix',1.0,0.5) # using our correlation matrix visualization function
plt.xlabel('Features'); plt.ylabel('Features')

plt.subplot(122)
feature_rank_plot(features[:-1],semipartial_correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Semipartial Correlation with ' + resp,'Semipartial Correlation',0.5)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=0.8, wspace=0.2, hspace=0.3); plt.show()

_images/e458ca94c5e55703a0d2ea41f4eaaa552252064a3f9dfb4bac57db86c6bf4d0a.png

For the semipartial correlation coefficients, we can observe,

porosity, permeability and vitrinite reflectance are the most important by semipartial correlation, similar to partial correlation
all other predictor features have quite low semipartial correlations compared to partial correlation

This is a good moment to visualize, interpret and compare of all the results from the correlation analysis for feature ranking with,

a plot of correlation, rank correlation and partial correlction.

plt.subplot(131)
feature_rank_plot(features[:-1],correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Correlation with ' + resp,'Correlation',0.5)

plt.subplot(132)
feature_rank_plot(features[:-1],rank_correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Rank Correlation with ' + resp,'Rank Correlation',0.5)

plt.subplot(133)
feature_rank_plot(features[:-1],partial_correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Partial Correlation with ' + resp,'Partial Correlation',0.5)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.3, top=0.8, wspace=0.3, hspace=0.15); plt.show()

_images/de56e10de16d7d23bc7e2ad9dc2737ab012cebfd72d2cd0b6bdf088414766df6.png

Based on the correlation analysis:

porosity and permeability both have high relevancy, but are also redundancy with eachother. Nevertheless, both remain highly ranked relative to the other predictor features, and domain expertise would favor retaining both
TOC has high relevancy, but also appears to have very high redundancy. Its partial correlation with production has flipped to negative and is small, suggesting that much of its apparent relationship with production may be redundant with information provided by the other predictor features
AI has moderate relevancy, but very little apparent redundancy and therefore brings new information to the feature set
VR has low relevancy and low redundancy, and ultimately ranks just below permeability

Feature Ranking with Feature Transformations#

There are many reasons to perform feature transformations,

see the feature transformation chapter for more information
also, as mentioned above for partial and semipartial correlation, a distribution transformation may assist with compliance with metric assumptions.

As an exercise and check for our previous covariance and correlation ranking, let’s transform all the features and repeat the previous calculations and visualizations, then compare the results.

we know this will have an impact on covariance, but what about the correlation and rank correlation coefficients?

For the workflow below, we:

make a new DataFrame with all variables standardized
apply the selected transformation to all features and store the results

You can choose between:

Standardization — affine correction to scale the distributions to have $\overline{x} = 0$ and $\sigma_x = 1.0$, without changing the shape of the distribution
Normal Score Transform — distribution transformation of each feature to a standard normal distribution, i.e., Gaussian distribution with $\overline{x} = 0$ and $\sigma_x = 1.0$.

Use this block to select between standardization and Gaussian transformation:

transformation = 0                                            # 0 - affine, standardization, 1 - Gaussian transformation

For standardization, also known as affine correction, the features are transformed to have a mean of 0.0 and a standard deviation of 1.0 without changing the shape of their distributions.

For Gaussian transformation, also known as the Normal Score Transform or Gaussian anamorphosis, the full distribution of each feature is transformed to a Gaussian shape (bell curve) with a mean of 0.0 and a standard deviation of 1.0:

transformation = 1                                            # 0 - affine, standardization, 1 - Gaussian transformation

transformation = 1                                            # 0 - affine, standardization, 1 - Gaussian transformation

dfS = pd.DataFrame()
dfS['Well'] = df['Well'].values                               # add well index

if transformation == 0:                                       # affine correction, standardization to a mean of 0 and variance of 1  
    for feature in features:  
        dfS[feature] = GSLIB.affine(df[feature].values,0.0,1.0)

elif transformation == 1:                                     # nscor transformation, Gaussian with a mean of 0 and variance of 1
    for feature in features:
        dfS[feature],d1,d2 = geostats.nscore(df,feature)   
    
dfS.head()

	Well	Por	Perm	AI	Brittle	TOC	VR	Prod
0	1	-0.964092	-0.780664	-0.285841	2.432379	0.312053	1.114651	-1.780464
1	2	-0.832725	-0.378580	0.446827	-0.195502	-0.272809	-0.325239	-0.392079
2	3	-0.312053	-1.069155	1.722384	2.004654	-0.272809	2.241403	-0.832725
3	4	0.730638	1.325516	-0.531604	-0.590284	0.131980	-0.325239	0.815126
4	5	0.698283	0.298921	0.365149	-2.870033	1.047216	-0.259823	-0.531604

Regardless of the transformation that you choose, it is best practice to check the summary statistics:

ensure that the target mean and standard deviation are achieved

dfS.describe()                                                # check the summary statistics

	Well	Por	Perm	AI	Brittle	TOC	VR	Prod
count	200.000000	200.000000	200.000000	2.000000e+02	2.000000e+02	200.000000	200.000000	2.000000e+02
mean	100.500000	-0.009700	0.010306	9.732356e-03	8.028717e-05	0.014152	0.017360	1.617223e-03
std	57.879185	1.040456	1.005488	1.000221e+00	1.000278e+00	0.989223	1.000401	9.949811e-01
min	1.000000	-4.991462	-3.355431	-2.782502e+00	-2.870033e+00	-2.336891	-2.899210	-2.483589e+00
25%	50.750000	-0.670577	-0.647337	-6.588985e-01	-6.705770e-01	-0.670577	-0.651072	-6.705770e-01
50%	100.500000	0.006267	0.006267	8.881784e-16	8.881784e-16	0.018807	0.006267	8.881784e-16
75%	150.250000	0.670577	0.678574	6.705770e-01	6.705770e-01	0.682378	0.682642	6.705770e-01
max	200.000000	2.807034	2.807034	2.807034e+00	2.807034e+00	2.807034	2.807034	2.807034e+00

We should also check the matrix scatter plot again.

If you performed a normal score transform, you have standardized the mean and variance and corrected the univariate shape of each distribution, but the bivariate relationships may still depart from a Gaussian distribution.

Warning: univariate Gaussian anamorphosis does not guarantee bivariate Gaussianity.

pairgrid = sns.PairGrid(dfS)                                  # matrix scatter plots
pairgrid = pairgrid.map_upper(plt.scatter, color = 'darkorange', edgecolor = 'black', alpha = 0.8, s = 10)
pairgrid = pairgrid.map_diag(plt.hist, bins = 20, color = 'darkorange',alpha = 0.8, edgecolor = 'k') # map a density plot to the lower triangle
pairgrid = pairgrid.map_lower(sns.kdeplot, cmap = plt.cm.inferno, 
                              shade = False, shade_lowest = False, alpha = 1.0, n_levels = 10)
pairgrid.add_legend()
plt.subplots_adjust(left=0.0, bottom=0.0, right=0.9, top=0.9, wspace=0.2, hspace=0.2); plt.show()

_images/49dcc0cdd1ffe6d82887553705ed24f69e8093a67644888796db83fc89968b52.png

This is the new DataFrame with standardized variables. Now we repeat the previous calculations:

we can be more efficient this time and use relatively compact code to calculate and then visualize the results.

Here’s the compact calculation.

stand_covariance = dfS.iloc[:,dfS.columns.get_indexer(features)].cov().iloc[len(features)-1,:len(features)]
stand_correlation = dfS.iloc[:,dfS.columns.get_indexer(features)].corr().iloc[len(features)-1,:len(features)]
stand_rank_correlation, stand_rank_correlation_pval = stats.spearmanr(dfS.iloc[:,dfS.columns.get_indexer(features)])
stand_rank_correlation = stand_rank_correlation[:,len(features)-1][:len(features)]
stand_partial_correlation = partial_corr(dfS.iloc[:,dfS.columns.get_indexer(features)]) # calculate the partial correlation coefficients
stand_partial_correlation = stand_partial_correlation[:,len(features)-1][:len(features)]
stand_semipartial_correlation = semipartial_corr(dfS.iloc[:,dfS.columns.get_indexer(features)]) # calculate the semi-partial correlation coefficients
stand_semipartial_correlation = stand_semipartial_correlation[:,len(features)-1][:len(features)] 

Now the compact visualziation code.

plt.subplot(2,3,1)
feature_rank_plot(features[:-1],correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Correlation with ' + resp,'Correlation',0.5)

plt.subplot(2,3,2)
feature_rank_plot(features[:-1],rank_correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Rank Correlation with ' + resp,'Rank Correlation',0.5)

plt.subplot(2,3,3)
feature_rank_plot(features[:-1],partial_correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Partial Correlation with ' + resp,'Partial Correlation',0.5)

plt.subplot(2,3,4)
feature_rank_plot(features[:-1],stand_correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Correlation with ' + resp,'Correlation of Standardized',0.5)

plt.subplot(2,3,5)
feature_rank_plot(features[:-1],stand_rank_correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Rank Correlation with ' + resp,'Rank Correlation of Standardized',0.5)

plt.subplot(2,3,6)
feature_rank_plot(features[:-1],stand_partial_correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Partial Correlation with ' + resp,'Partial Correlation of Standardized',0.5)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.3, top=1.6, wspace=0.3, hspace=0.4); plt.show()

_images/46b8f282715840779b17f6faa885b1b619e448ce807d3f1da7d6d342f9e63086.png

What is the impact of our feature transformation? We can observe:

porosity now stands out as the strongest predictor feature, with high relevancy and low redundancy
permeability and AI remain highly relevant, but both now have high redundancy

Conditional Statistics and Conditional Distributions#

We separate the wells into low- and high-production groups and assess the differences in the conditional statistics and conditional distributions.

this provides a flexible, non-parametric method without distributional model assumptions to compare the relationship between each predictor feature and the response feature, in our example production
if the conditional statistics or conditional distributions differ substantially between the low- and high-production groups, then that feature is potentially informative

The steps for visualizing conditional statistics and distributions are:

apply a binary transformation to the response feature to form a new categorical feature. In this example, we apply a production threshold to transform the production response feature into “low” and “high” categories:

df['tProd'] = np.where(df['Prod']>=4000, 'High', 'Low') 

standardize all of our features to observe their relative differences together, i.e., so they all fit on the same y-axis:

x = df[['Por','Perm','AI','Brittle','TOC','VR']]
x_stand = (x - x.mean()) / (x.std())      

note, this code extracts the features into a new DataFrame x, then applies the standardization operation independently to each column (feature)

form a new, convenient DataFrame with the new response categorical feature and the standardized predictor features:

x = pd.concat([df['tProd'],x_stand.iloc[:,0:6]],axis=1)

apply the melt command to unpivot the DataFrame:

x = pd.melt(x,id_vars="tProd",var_name="Predictors",value_name='Standardized_Value')

now the DataFrame (6 features x 200 samples) is represented with 12000 rows, with:
- production: Low or High
- features: Por, Perm, AI, Brittle, TOC or VR
- standardized feature value

now pass this DataFrame to the violinplot or boxplot function from seaborn:

x is our predictor features
y is the standardized value for each predictor feature (all now in one column)
hue is the production level, High or Low
split is True so the violins are split in half
inner is "quart" so the 25th, 50th, and 75th percentiles are plotted

threshold = 2000.0                                            # select low and high threshold for production

df['tProd'] = np.where(df[resp] >= threshold, 'High', 'Low')  # temporary feature with low and high labels

x_temp = df[pred]                                             # temporary dataframe with feature, low / high label and standardized features
x_temp_stand = (x_temp - x_temp.mean()) / x_temp.std()        # standardize features
x_temp = pd.concat([df['tProd'], x_temp_stand.iloc[:, 0:len(pred)]],axis=1)
x_temp = pd.melt(x_temp,id_vars="tProd",var_name="Predictor Feature",value_name='Standardized Predictor Feature') # unpivot the DataFrame

plt.subplot(111)
palette = {"Low": "#0057B7", "High": "#FFD700"}               # fill colors

sns.violinplot(x="Predictor Feature",y="Standardized Predictor Feature",hue="tProd",data=x_temp,split=True,inner="quart",width=0.4,palette=palette)

plt.legend(title=None,loc="lower left",ncol=2,frameon=True)   # horizonal stacked legend 

for i in range(len(pred) - 1):                                # vertical separators between predictor features
    plt.axvline(i + 0.5,color='black',linestyle='-',linewidth=0.5,alpha=0.25)

plt.xticks(rotation=90)
plt.title('Conditional Standardized Predictor Feature Distributions Given Response Feature Threshold = ' + str(threshold))

ax = plt.gca()
ax.yaxis.set_major_locator(MultipleLocator(2.0))               # major horizontal grid lines every 2.0
ax.yaxis.set_minor_locator(MultipleLocator(0.2))               # minor horizontal grid lines every 0.2
ax.grid(axis='y', which='major', linewidth=0.8, alpha=0.5)     # grid lines
ax.grid(axis='y', which='minor', linewidth=0.4, alpha=0.2)

plt.subplots_adjust(left=0.0, bottom=0.0, right=1.5, top=1.0,wspace=0.2, hspace=0.2); plt.show()

df = df.drop(['tProd'], axis=1)

_images/0a8048b24ac08e91a80ff1f6e70fe7b463a00a0dfbbbe86cb6f7487db40b6624.png

From the violin plot we can observe:

the conditional distributions of porosity, permeability, and TOC show the largest differences between low- and high-production wells, including differences in their quartiles and modes.

We can replace the seaborn violinplot with the seaborn boxplot to display box-and-whisker plots of the conditional distributions.

box-and-whisker plots provide a clear display of the conditional P25, P50, and P75, along with the whiskers extending to the most extreme values within the Tukey outlier fences for low- and high-production wells.

threshold = 2000.0                                            # select low and high threshold for production

df['tProd'] = np.where(df[resp] >= threshold, 'High', 'Low')  # temporary feature with low and high labels

x_temp = df[pred]                                             # temporary dataframe with feature, low / high label and standardized features
x_temp_stand = (x_temp - x_temp.mean()) / x_temp.std()        # standardize features
x_temp = pd.concat([df['tProd'], x_temp_stand.iloc[:, 0:len(pred)]],axis=1)
x_temp = pd.melt(x_temp,id_vars="tProd",var_name="Predictor Feature",value_name='Standardized Predictor Feature') # unpivot the DataFrame

plt.subplot(111)
palette = {"Low": "#0057B7", "High": "#FFD700"}               # fill colors

# sns.violinplot(x="Predictor Feature",y="Standardized Predictor Feature",hue="tProd",data=x_temp,split=True,inner="quart",palette=palette)
sns.boxplot(x="Predictor Feature", y="Standardized Predictor Feature", hue="tProd",data=x_temp,width=0.4,palette=palette)

plt.legend(title=None,loc="lower left",ncol=2,frameon=True)   # horizonal stacked legend 

for i in range(len(pred) - 1):                                # vertical separators between predictor features
    plt.axvline(i + 0.5,color='black',linestyle='-',linewidth=0.5,alpha=0.25)

plt.xticks(rotation=90)
plt.title('Conditional Standardized Predictor Feature Distributions Given Response Feature Threshold = ' + str(threshold))

ax = plt.gca()
ax.yaxis.set_major_locator(MultipleLocator(2.0))               # major horizontal grid lines every 2.0
ax.yaxis.set_minor_locator(MultipleLocator(0.2))               # minor horizontal grid lines every 0.2
ax.grid(axis='y', which='major', linewidth=0.8, alpha=0.5)     # grid lines
ax.grid(axis='y', which='minor', linewidth=0.4, alpha=0.2)

plt.subplots_adjust(left=0.0, bottom=0.0, right=1.5, top=1.0,wspace=0.2, hspace=0.2); plt.show()

df = df.drop(['tProd'], axis=1)

_images/89fc9c5d81f4b53ff5cd0a371d19a360c9c535c237c439499643d19c47529ac8.png

From the conditional boxplot we can observe that the conditional distributions of porosity, permeability, and total organic carbon have the largest differences between low- and high-production wells.

we can also specifically observe outliers in permeability and total organic carbon above the upper Tukey outlier bound
vitrinite reflectance has outliers both below the lower Tukey outlier bound and above the upper Tukey outlier bound

Variance Inflation Factor#

A measure of linear multicollinearity between a predictor feature ($X_i$) and all other predictor features ($X_j, \forall j \ne i$).

Steps to calculate the variance inflation factor:

calculate a linear regression for a predictor feature given all the other predictor features:

\[\begin{split} X_i = \beta_0 + \sum_{\substack{j=1\\j\ne i}}^m \beta_j X_j + b_0 \end{split}\]

from the above model, determine the coefficient of determination, $R^2$, which measures the proportion of variance in $X_i$ explained by the other predictor features. Then, calculate the Variance Inflation Factor as:

\[ VIF = \frac{1}{1 - R^2} \]

repeat steps 1 and 2 for all predictor features.

How do we interpret variance inflation factor? Here’s a table that may be helpful

VIF	$R_X^2$	Interpretation	Guidance
5	0.80	Substantial redundancy	Warning, investigate feature
10	0.90	Strong redundancy	Candidate for removal
20	0.95	Severe redundancy	Strong candidate for removal
40	0.98	Very severe redundancy	Consider removing the feature unless there is a strong reason to retain it

Note that the variance inflation factor only accounts for redundancy among predictor features and does not account for the relevance of a predictor feature to the response feature.

commonly applied as a feature screening tool rather than for comprehensive feature ranking and selection
features with very high variance inflation factors may be dropped prior to feature ranking with other methods

vif_values = []
for i in range(df[pred].values.shape[1]):
    vif_values.append(variance_inflation_factor(df[pred].values, i))

vif_values = np.asarray(vif_values); indices = np.argsort(vif_values)[::-1] # find indices for descending order               

bar_colors = []                                               # define colors by VIF level
for vif in vif_values[indices]:
    if vif < 5:
        bar_colors.append("seagreen")
    elif vif < 20:
        bar_colors.append("gold")
    elif vif < 40:
        bar_colors.append("darkorange")
    else:
        bar_colors.append("firebrick")

plt.subplot(111)
plt.title("Variance Inflation Factor")
plt.bar(range(len(vif_values)),vif_values[indices],edgecolor="black",color=bar_colors,alpha=0.6,align="center")
plt.axhline(5, linestyle="--",linewidth=1.5,color='black',label="5: Warning") # add VIF thresholds
plt.axhline(20, linestyle="--",linewidth=1.5,color='black',label="10: Removal candidate")
plt.axhline(40, linestyle="--",linewidth=1.5,color='black',label="20: Remove")

plt.xticks(np.arange(len(pred)),np.array(pred)[indices].tolist(),rotation=90)

plt.gca().yaxis.grid(True, which="major", linewidth=1.0)
plt.gca().yaxis.grid(True, which="minor", linewidth=0.2)
plt.gca().tick_params(which="major", length=7)
plt.gca().tick_params(which="minor", length=4)

plt.gca().xaxis.set_minor_locator(AutoMinorLocator())
plt.gca().yaxis.set_minor_locator(AutoMinorLocator())

plt.xlim([-0.5, len(pred)-0.5])
plt.yscale("log")

plt.xlabel("Predictor Feature")
plt.ylabel("Variance Inflation Factor")

legend_elements = [                                           # legend for VIF categories
    Patch(facecolor="seagreen", edgecolor="black", alpha=0.6, label="< 5: OK"),
    Patch(facecolor="gold", edgecolor="black", alpha=0.6, label="5–20: Warning"),
    Patch(facecolor="darkorange", edgecolor="black", alpha=0.6, label="20–40: Removal candidate"),
    Patch(facecolor="firebrick", edgecolor="black", alpha=0.6, label="> 40: Remove")
]

plt.legend(handles=legend_elements, loc="upper right")

plt.subplots_adjust(left=0.0, bottom=0.0, right=2., top=1.,wspace=0.2, hspace=0.5); plt.show()

_images/b87d2f718e23a53bb042d783fe45935169c944107c6111ef70e95e1b5b74dfd0.png

From our variance inflation factor reverse-order-sorted bar chart, we can observe:

vitrinite reflectance has the most linear redundancy, while permeability has the least linear redundancy with the other predictor features.
remember, a high variance inflation factor indicates high linear redundancy!

Warning: the thresholds and sorted feature categories are a demonstration of a possible scheme and communication tool for feature screening with variance inflation factor. VIF values of 5, 10, and 40 are used here as thresholds for moderate redundancy, high redundancy and very high redundancy, respectively. The appropriate thresholds depend on the specific application and should not be interpreted as strict rules for feature removal.

Warning: variance inflation factor only evaluates linear redundancy, i.e., relationships with the other predictor features, and does not evaluate relevancy, i.e., the linear relationship with the response feature. Also, variance inflation factor is used as a screening tool to identify features with very high redundancy, and not as a comprehensive feature ranking tool.

Now let’s cover model-based feature ranking methods.

Model-based Feature Ranking Methods#

Some predictive machine learning models provide their own, built-in measures of feature importance. For example:

linear regression — the slope coefficients provide the gradient, or rate of change, of the response feature given a change in each predictor feature. The coefficients can be compared as measures of feature importance when the predictor features have been standardized to a common scale.
decision tree — the summation of the reduction in regression error (or classification impurity) from splitting on each predictor feature, then transformed by an L1 normalization to sum to 1.0, provides an impurity-based measure of the relative importance of each predictor feature in the prediction model.
Shapley value — a model-agnostic method based on summarizing marginal contributions, i.e., the impact on the prediction of including a feature versus excluding the feature, across different combinations of the other features and over a variety of prediction examples.

Note, there are workflows for neural network feature importance based on network weights, but for reasonably complicated networks, interpreting feature importance directly from the network weights is quite impractical.

Warning: model-based feature ranking is always dependent on the quality of the prediction model. Always first check the prediction model before using feature importance from the model.

B Coefficients#

Our first model-based feature importance measure is known as $B$ coefficients. These are:

linear regression model coefficients
without standardization of the variables
sensitive to the units, scale, and variance of the predictor features

Linear regression model parameters can be used for feature ranking by accounting for,

relevance, the linear relationship between the predictor feature and the response feature
redundancy, the linear relationships between the predictor feature and the other predictor features

Let’s use the linear regression method that is available in the SciPy Python package.

The estimator for $Y$ is simply the linear equation:

\begin{equation} Y^* = \sum_{i=1}^{m} b_i X_i + b_0 \end{equation}

The $b_i$ coefficients are solved to minimize the squared error between the estimates, $Y^*$, and the values in the training dataset, $Y$.

Warning: due to the sensitivity to the units, scale, and variance of the predictor features, $\beta$ coefficients (beta coefficients) instead of $B$ coefficients are commonly used for feature importance. Only use $B$ coefficients when the features have consistent units, or when differences in the variance and scale of the predictor features are meaningful.

reg = LinearRegression()                                      # instantiate a linear regression model                   
reg.fit(df[pred],df[resp])                                    # train the model
b = reg.coef_
Y_hat = reg.predict(df[pred])

plt.subplot(121)                                              # check model over training data, for brevity not a holdout cross validation
plot_prediction_check(df[resp], Y_hat, threshold=1500.0, ymin=0.0, ymax=9000,xlabel='Actual Production (MCFPD)',ylabel='Estimated Production (MCFPD)',
                     title='Model Cross Validation Plot - Training Data Only')
add_grid()

plt.subplot(122)
feature_rank_plot(pred,b,-1000.0,1000.0,0.0,'Feature Ranking, Linear Regression B Coefficients with ' + resp,r'Linear Regression Slope, $b_1$',0.5)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2., top=1.,wspace=0.1, hspace=0.3); plt.show()

_images/90ba93d73b541f31e7e8b36794b6b05ae91c5351962d98edf84c5ee96dda0385.png

Prior to using any model-based feature ranking measure, it is critical to first check the model.

the model cross validation plot includes residual summary statistics and highlights samples with large residuals for further investigation
note, this is a simple cross validation over the training data, i.e., the data used to train the model, to check the model fit to the training data
a check of model generalization, i.e., prediction for cases not used to train the model, requires holdout cross validation or K-fold cross validation

Here is the $b$ coefficients bar chart, ordered over our features from $b_i$, $i = 1,\ldots,n$. We can observe,

we see the negative contribution of AI and TOC

Warning: $b$ coefficients are sensitive to the magnitude, variance, and units of the predictor features and therefore should not be used directly for feature ranking when the predictor features are on different scales. We can remove this sensitivity by working with standardized features.

The $b$ coefficient feature ranking results are very sensitive to the scale and variance of the predictor features, and the units of the predictor features are not consistent. Therefore, we ignore these results and continue to $\beta$ coefficients.

Beta Weights#

Our second model-based feature importance measure is $\beta$ coefficients. Buliding on $b$ coefficients, these are:

also, linear regression model coefficients
but with standardization of the variables

Due to prior feature standardization to a mean of 0.0 and variance of 1.0, these are,

not sensitive to the absolute and relative magnitude and variance of the predictor features

Like $b$ coefficients, $\beta$ coefficients from linear regression model parameters can be used for feature ranking by accounting for,

relevance, the linear relationship between the predictor feature and the response feature
redundancy, the linear relationships between the predictor feature and the other predictor features

Let’s use the linear regression method that is available in the SciPy Python package.

The estimator for $Y^s$ standardized, with predictor features standardized, $X^s_i$, is simply the linear equation:

\begin{equation} \hat{Y}^{s} = \sum_{i=1}^{m} \beta_i X^s_i \end{equation}

Notice that the constant term, usually denoted as $b_0$ or $c$ is missing? Given the predictor and response features are standardized,

\[ \bar{y}_z=0,\qquad \bar{x}_{i,z}=0 \]

and for multiple regression, the intercept is,

\[ b_0=\bar{y}-\sum_{i=1}^{p}b_i\bar{x}_i \]

Therefore, after standardization of both predictor and response features, the intercept must be,

\[ b_0=0-\sum_{i=1}^{p}b_i(0)=0 \]

Warning: Assuming the multilinear regression constant term is 0.0 requires that both the response feature, $Y$, and all predictor features, $X$ have been standardized. If only the predictor features are standardized, the intercept is generally not zero.

We are ready to calculate $\beta$ coefficients since we have just standardized all our variables, note,

standardization $\rightarrow$ $\beta$ coefficients, same data relationships but with a common scale.

If you selected Gaussian anamorphosis above,

nscore transform $\rightarrow$ linear regression coefficients in transformed Gaussian space, with common scale and transformed distribution.

Warning: while Gaussian anamorphosis does set the mean and variance to 0.0 and 1.0 respectively and the resulting multilinear coefficients may be useful, they are not technically $\beta$ coefficients.

reg = LinearRegression()                                      # instantiate linear regression model
reg.fit(dfS[pred],dfS[resp])                                  # train model
beta = reg.coef_                                              # extract model parameters
Y_hat_beta = reg.predict(dfS[pred])                           # predict at training data

plt.subplot(121)                                              # check model over training data, for brevity not a holdout cross validation
plot_prediction_check(dfS[resp], Y_hat_beta, threshold=1.0, ymin=-3.0, ymax=3.0,
                    xlabel='Actual Production (Standardized MCFPD)',ylabel='Estimated Production (Standardized MCFPD)',
                    title='Model Cross Validation Plot - Training Data Only')
add_grid()

plt.subplot(122)
feature_rank_plot(pred,beta,-1.0,1.0,0.0,'Feature Ranking, Linear Regression B Coefficients with ' + resp,r'Linear Regression Slope, $b_1$',0.5)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2., top=1.,wspace=0.1, hspace=0.3); plt.show()

_images/7ce904c5b70b0250f93f667d1e3f91c3db475a5f766d6ca46f1b2d6de9dba67e.png

We can observe:

the change between the $b$ and $\beta$ coefficients is not simply a constant scaling of the ranking metrics, because the linear model coefficients are sensitive to the range and magnitude of the current predictor feature and the other predictor features
with the $\beta$ coefficients, porosity, acoustic impedance, and vitrinite reflectance are indicated as the most important features
the rankings of acoustic impedance and TOC have decreased significantly compared with the $b$ coefficient results

Feature Importance#

A variety of machine learning methods provide measures of feature importance, for example,

decision trees summarize the reduction in mean square error through inclusion of each feature as:

\[ FI(x) = \sum_{t \in T_f} \frac{N_t}{N} \Delta_{MSE_t} \]

where $T_f$ are all nodes with feature $x$ as the split, $N_t$ is the number of training samples reaching node $t$, $N$ is the total number of samples in the dataset and $\Delta_{MSE_t}$ is the reduction in MSE with the $t$ split.

Note, feature importance can be calculated in a similar manner to MSE reduction summarization above for the case of classification trees with Gini Impurity.

Let’s look at the feature importance from a random forest regression model fit to our data.

We instantiate a random forest with default hyperparameters. This results in unlimited complexity, over-trained trees in our forest. The averaging of these trees takes care of the overfit issue.
Then we train our random forest and extract the feature importances, calculated as the expectated feature importance over all the trees in the forest.
we can also extract the feature importances over all the trees in the forest and summarize with the standard deviation to access the robustness, uncertainty of our feature importance measure

For more information check out my lecture on random forest predictive machine learning.

random_forest = RandomForestRegressor(n_estimators=10,max_leaf_nodes=20) # instantiate the random forest  
random_forest = random_forest.fit(df[pred],df[resp])          # fit the random forest
importance_rank = random_forest.feature_importances_          # extract the expected feature importances
importance_rank_stand = importance_rank/np.max(importance_rank) # calculate relative mutual information

std = np.std([tree.feature_importances_ for tree in random_forest.estimators_],axis=0) # calculate stdev over trees
indices = np.argsort(importance_rank)[::-1]                   # find indices for descending order

Y_hat_tree = random_forest.predict(df[pred])

plt.subplot(121)                                              # check model over training data, for brevity not a holdout cross validation
plot_prediction_check(df[resp], Y_hat_tree, threshold=1500.0, ymin=0.0, ymax=9000.0,
                    xlabel='Actual Production (MCFPD)',ylabel='Estimated Production (MCFPD)',
                    title='Model Cross Validation Plot - Training Data Only')
add_grid()

plt.subplot(122)                                              # plot the feature importance 
baseline = 1.0 / x.shape[1]                                   # equal-share reference

bar_colors = ['darkorange' if v >= baseline else 'white'
              for v in importance_rank[indices]]

plt.title("Random Forest-based Feature Importances")
plt.bar(range(x.shape[1]), importance_rank[indices],
        edgecolor='black', color=bar_colors, alpha=0.8,
        yerr=std[indices], align='center')

plt.xticks(range(x.shape[1]), x.columns[indices], rotation=90)
plt.xlim([-0.5, x.shape[1]-0.5]); plt.ylim([0.,1.0])

plt.axhline(baseline, color='black', linestyle='--', linewidth=1.0) # equal-share baseline

plt.gca().yaxis.grid(True, which='major', linewidth=1.0)
plt.gca().yaxis.grid(True, which='minor', linewidth=0.2)
plt.gca().tick_params(which='major', length=7)
plt.gca().tick_params(which='minor', length=4)
plt.gca().xaxis.set_minor_locator(AutoMinorLocator())
plt.gca().yaxis.set_minor_locator(AutoMinorLocator())

plt.xlabel('Predictor Feature'); plt.ylabel('Feature Importance')

legend_elements = [                                           # legend
    Patch(facecolor='darkorange', edgecolor='black', alpha=0.8, label='Strong'),
    Patch(facecolor='white', edgecolor='black', label='Weak')
]
plt.legend(handles=legend_elements, title='Feature Importance',loc='upper right')

plt.subplots_adjust(left=0.0, bottom=0.0, right=2., top=1., wspace=0.1, hspace=0.5); plt.show()

_images/f8841e38d31a4749bdedaf8c5ed0d5dc779df73eb01b36d5d67a4b88ee527fd0.png

We can observe:

this is our best-performing model for accuracy over the training data, giving us increased confidence in the model-based feature ranking results

Tree-based feature importance should generally be interpreted relative to the other predictor features rather than using an absolute threshold,

tree-based feature importance is L1-normalized to sum to 1.0

As a rule of thumb, features with importance substantially greater than the typical feature importance can be considered strong contributors, while features with importance near zero can be considered weak contributors,

the equal-share baseline, $1/n$, provides a simple reference when the feature importances are normalized to sum to one
we suggest features above this baseline as strong features and features below this baseline as weak features as a useful communication tool

Here’s the tree-based feature importance from the random forest, appended to the previous feature ranking results from,

correlation analysis
$B$ and $\beta$ coefficients

plt.subplot(321)
feature_rank_plot(features[:-1],rank_correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Rank Correlation with ' + resp,'Rank Correlation',0.5)

plt.subplot(322)
feature_rank_plot(features[:-1],partial_correlation[:-1],-1.0,1.0,0.0,'Feature Ranking, Partial Correlation with ' + resp,'Partial Correlation',0.5)

plt.subplot(323)
feature_rank_plot(pred,b[0:len(pred)],-1000.0,1000.0,0.0,'Feature Ranking, B Coefficients with ' + resp,'B Coefficients',0.5)

plt.subplot(324)
feature_rank_plot(pred,beta[0:len(pred)],-1.0,1.0,0.0,'Feature Ranking, Beta Coefficients with ' + resp,'Beta Coefficients',0.5)

plt.subplot(325)
feature_rank_plot(pred,importance_rank_stand,0.0,1.0,0.0,'Feature Ranking, Feature Importance with ' + resp,'Standardized Feature Importance',0.5)

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.2, top=3.2, wspace=0.2, hspace=0.3); plt.show()

_images/0791ec61375338a592786ce11115c433e7aaa792b3ff80540b16f21cbcf2aaae.png

Mutual Information#

Mutual information is a generalized approach that quantifies the mutual dependence between two features.

quantifies the amount of information gained from observing one feature about the other
avoids any assumption about the form of the relationship (e.g. no assumption of linear relationship)
compares the joint probabilities to the product of the marginal probabilities

For discrete or binned continuous features $X$ and $Y$, mutual information is calculated as:

\[ I(X;Y) = \sum_{y \in Y} \sum_{x \in X}P_{X,Y}(x,y) log \left( \frac{P_{X,Y}(x,y)}{P_X(x) \cdot P_Y(y)} \right) \]

Note, given independence between $X$ and $Y$:

\[ P_{X,Y}(x,y) = P_X(x) \cdot P_Y(y) \]

therefore if the two features are independent then the $log \left( \frac{P_{X,Y}(x,y)}{P_X(x) \cdot P_Y(y)} \right) = 0$

The joint probability $P_{X,Y}(x,y)$ is a weighting term on the sum and enforces closure.

parts of the joint distribution with greater density have greater impact on the mutual information metric

For continuous (and nonbinned) features we can applied the integral form.

\[ I(X;Y) = \int_{Y} \int_{X}P_{X,Y}(x,y) log \left( \frac{P_{X,Y}(x,y)}{P_X(x) \cdot P_Y(y)} \right) dx dy \]

We get a sorted list of the indices in decreasing order of importance with the command

indices = np.argsort(importances)[::-1]

the slice reverses the order, for descending order of feature importance.

x_df = df.loc[:,pred]                                         # separate DataFrames for predictor and response features
y_df = df.loc[:,resp]

mi = mutual_info_regression(x_df,np.ravel(y_df))              # calculate mutual information
mi /= np.max(mi)                                              # calculate relative mutual information

indices = np.argsort(mi)[::-1]                                # find indices for descending order

print("Feature ranking:")                                     # write out the feature importances
for f in range(x.shape[1]):
    print("%d. feature %s = %f" % (f + 1, x.columns[indices][f], mi[indices[f]]))

plt.subplot(111)                                              # plot the relative mutual information 
plt.title("Mutual Information")
plt.bar(range(x.shape[1]), mi[indices],edgecolor = 'black',
       color="darkorange",alpha=0.6,align="center")
plt.xticks(range(x.shape[1]), x.columns[indices],rotation=90)
plt.xlim([-0.5, x.shape[1]-0.5]); plt.ylim([0,1.3])
plt.gca().yaxis.grid(True, which='major',linewidth = 1.0); plt.gca().yaxis.grid(True, which='minor',linewidth = 0.2) # add y grids
plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)
plt.gca().xaxis.set_minor_locator(AutoMinorLocator()); plt.gca().yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks
plt.xlabel('Predictor Feature'); plt.ylabel('Mutual Information')

plt.subplots_adjust(left=0.0, bottom=0.0, right=2., top=1., wspace=0.2, hspace=0.5); plt.show()

Feature ranking:
feature Por = 1.000000
feature Perm = 0.350161
feature TOC = 0.272141
feature Brittle = 0.073684
feature AI = 0.055586
feature VR = 0.002482

_images/8381acb8df622d0316449db3bf8fd4f6681f7e8238dcd74093cd4578a95e2a6d.png

Mutual information feature importance should generally be interpreted relative to the other predictor features rather than using an absolute threshold,

mutual information is normalized by the maximum mutual information to give a maximum value of 1.0
the feature with the maximum mutual information therefore has a normalized importance of 1.0, while features with lower mutual information have values between 0 and 1

As a rule of thumb, features with normalized mutual information substantially closer to 1.0 can be considered strong contributors, while features with values near 0.0 can be considered weak contributors.

unlike tree-based feature importance, there is no simple equal-share baseline such as $1/n$ for maximum-normalized mutual information
the relative distribution of the normalized mutual information values provides a useful communication tool for identifying strong and weak features

Mutual Information Accounting For Relevance and Redundancy#

Mutual information for feature ranking only accounts for,

feature relevance, through mutual dependence between the predictor feature and the response feature

but does not account for,

feature redundancy, through mutual dependence between the predictor feature and the other predictor features

Maximum Relevance Minimum Redundancy (MRMR) is a feature ranking approach that applies mutual information and accounts for both relevance and redundancy. For a subset of features, an MRMR measure is calculated as,

average mutual information between the predictor features in the subset and the response feature, minus the average mutual information between the predictor features in the subset,

\begin{equation} MID = \frac{1}{|S|}{\sum_{\alpha \in S} I(X_{\alpha},Y) } - \frac{1}{|S|^2} {\sum_{\alpha \in S} \sum_{\beta \in S} I(X_{\alpha},X_{\beta})} \end{equation}

This measure is posed as $relevance - redundancy$. An alternative MRMR measure is,

average mutual information between the predictor features in the subset and the response feature, divided by the average mutual information between the predictor features in the subset,

\begin{equation} MIQ = \frac{ \frac{1}{|S|}{\sum_{\alpha \in S} I(X_{\alpha},Y) } }{ \frac{1}{|S|^2} {\sum_{\alpha \in S} \sum_{\beta \in S} I(X_{\alpha},X_{\beta})} } \end{equation}

as a measure of $\frac{relevance}{redundancy}$.

The MRMR workflow includes evaluating the MRMR measure over various combinations of predictor features and selecting the combination that maximizes the MRMR measure.

Mutual Information Accounting For Relevance and Redundancy OFAT Variants#

Merzoug et al. (in press) propose an MRMR-based method with one-feature-at-a-time (OFAT) predictor feature ranking. For OFAT ranking, each predictor feature is considered individually as the selected feature, while the remaining predictor features are used to evaluate redundancy. We modify this approach to the following calculation:

relevance - the mutual information between the selected predictor feature, $X_i$, and the response feature, $Y$
redundancy - the average mutual information between the selected predictor feature, $X_i$, and the remaining predictor features, $X_{\alpha}, \alpha \ne i$
with the quotient form of the calculation from Gulgezen, Cataltepe and Yu (2009)

This modified version of the Maximum Relevance - Minimum Redundancy (MRMR) objective function for OFAT ranking scores the selected predictor feature $X_i$’s relevance as its mutual information with the response feature,

\begin{equation} I(X_i,Y) \end{equation}

and redundancy as the average mutual information between the selected predictor feature, $X_i$, and the remaining predictor features,

\begin{equation} \frac{1}{m-1} \sum_{\substack{\alpha=1 \ \alpha \ne i}}^m I(X_i,X_{\alpha}) \end{equation}

where $X$ are predictor features, $Y$ is the response feature, $X_i$ is the specific predictor feature being scored, $m$ is the number of predictor features, and $I()$ is mutual information between the indicated features.

One formulation is a simple difference, relevance minus redundancy,

\[\begin{split} \Phi_{\Delta}(X_i,Y) = I(X_i,Y) - \frac{1}{m-1} \sum_{\substack{\alpha=1 \\ \alpha \ne i}}^m I(X_i,X_{\alpha}) \end{split}\]

An alternative is a ratio,

\[\begin{split} \Phi_{r}(X_i,Y) = \frac{ I(X_i,Y) }{ \frac{1}{m-1} \sum_{\substack{\alpha=1 \\ \alpha \ne i}}^m I(X_i,X_{\alpha})} \end{split}\]

Here are the feature ranks for the mutual information relevance minus redundancy, $\Phi_{\Delta}(X_i,Y)$, approach.

obj_mutual = mutual_information_objective(x_df,y_df)
indices_obj = np.argsort(obj_mutual)[::-1]                    # find indices for descending order

plt.subplot(111)                                              # plot the relative mutual information 
plt.title("OFAT MRMR Objective Function for Mutual Information-based Feature Selection")
plt.bar(range(x.shape[1]), obj_mutual[indices_obj],
       color="darkorange",alpha = 0.6, align="center",edgecolor="black")
plt.xticks(range(x.shape[1]), x.columns[indices_obj],rotation=90)
plt.gca().yaxis.grid(True, which='major',linewidth = 1.0); plt.gca().yaxis.grid(True, which='minor',linewidth = 0.2) # add y grids
plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)
plt.gca().xaxis.set_minor_locator(AutoMinorLocator()); plt.gca().yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks
plt.xlim([-0.5, x.shape[1]-0.5]); plt.xlabel('Predictor Feature'); plt.ylabel('Feature Importance')
plt.subplots_adjust(left=0.0, bottom=0.0, right=2., top=1., wspace=0.2, hspace=0.5)
plt.show()

_images/ff8cd5b9084d119834424e20e0b668b7a461a549b02c6fa819c3bc001896f1fd.png

Delta Mutual Information Quotient Accounting for Relevance and Redundancy#

Gulgezen, Cataltepe and Yu (2009) develop the mutual information quotient as an OFAT ranking metric.

Akmal et al. (in press) suggest the use of a delta measure with the mutual information quotient

The change in $MIQ$ from the inclusion and removal of a specific predictor feature, $X_i$, may be calculated as,

\begin{equation} \Delta MIQ_i = \frac{ \frac{1}{|S|}{\sum_{\alpha=1}^m I(X_{\alpha},Y) } }{ \frac{1}{|S|^2} {\sum_{\alpha=1}^m \sum_{\beta=1}^m I(X_{\alpha},X_{\beta})} } - \frac{ \frac{1}{|S|}{\sum_{\alpha=1,\alpha \ne i}^m I(X_{\alpha},Y) } }{ \frac{1}{|S|^2} {\sum_{\alpha=1,\alpha \ne i}^m \sum_{\beta=1,\beta \ne i}^m I(X_{\alpha},X_{\beta})} } \end{equation}

delta_mutual_information = delta_mutual_information_quotient(x_df,y_df)

indices_delta_mutual_information = np.argsort(delta_mutual_information)[::-1] # find indices for descending order

plt.subplot(111)                                              # plot the relative mutual information 
plt.title("Delta Mutual Information Quotient")
plt.bar(range(x.shape[1]), delta_mutual_information[indices_delta_mutual_information],
       color="darkorange",alpha = 0.6,align="center",edgecolor = 'black')
plt.xticks(range(x.shape[1]), x.columns[indices_delta_mutual_information],rotation=90)
plt.xlim([-0.5, x.shape[1]-0.5])
plt.gca().yaxis.grid(True, which='major',linewidth = 1.0); plt.gca().yaxis.grid(True, which='minor',linewidth = 0.2) # add y grids
plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)
plt.gca().xaxis.set_minor_locator(AutoMinorLocator()); plt.gca().yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks
plt.plot([-0.5,x.shape[1]-0.5],[0,0],color='black',lw=3); plt.xlabel('Predictor Feature'); plt.ylabel('Feature Importance')
plt.subplots_adjust(left=0.0, bottom=0.0, right=2., top=1., wspace=0.2, hspace=0.5)
plt.show()

_images/4cb4bc8342ea26d8c1890f88a411dda2c9d33b0528fe3a4fdd7dca583cfa15f2.png

It is instructive to compare delta mutual information and variance inflation factor (VIF) ranking. Both of these methods account for predictor feature redundancy.

but VIF assumes linear relationships between predictor features and does not account for feature relevance to the response feature

The plot below may be applied to compare and contrast any feature ranks,

the axes are the ranks in each ranking metric
the intervals identify inconsistency in feature ranks between ranking metrics

threshold = 1.0                                               # rank threshold for significant difference
ymin = 0; ymax = len(pred) + 1

plt.scatter(stats.rankdata(delta_mutual_information),stats.rankdata(-vif_values),c='black',edgecolor='black')
for i, feature in enumerate(x.columns):
    plt.annotate(feature, (stats.rankdata(delta_mutual_information)[i]-0.2,stats.rankdata(-vif_values)[i]+0.1))
plt.xlabel('Delta Mutual Information Rank'); plt.ylabel('Variance Inflation Factor Rank')
plt.title('Variance Inflation Factor vs. Delta Mutual Information Feature Ranking')
plt.xlim(0,len(pred)+1.0); plt.ylim(0,len(pred)+1.0)

plt.plot([ymin, ymax],[ymin + threshold, ymax + threshold],color='black', linestyle='--', linewidth=0.8,zorder=100) # residual threshold: y_hat = y ± threshold
plt.plot([ymin, ymax],[ymin - threshold, ymax - threshold],color='black', linestyle='--', linewidth=0.8,zorder=100)

plt.fill_between([ymin, ymax],[ymin + threshold, ymax + threshold], [ymax, ymax], color='coral',alpha=0.2,zorder=1)
plt.fill_between([ymin, ymax],[ymin - threshold, ymax - threshold], [ymin, ymin], color='dodgerblue',alpha=0.2,zorder=1)

plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.0, wspace=0.2, hspace=0.2); plt.show()

_images/4ded8275579ad2186bce108c62c58e1b7760bd2d846b97accab47e8df440fb8b.png

After comparing VIF and delta mutual information, we can observe:

delta mutual information ranks porosity significantly higher than VIF, likely because VIF accounts for predictor feature redundancy but does not account for the strong relevance of porosity to the response feature.

Summary of All Ranking Methods#

Now we have a wide array of criteria to rank our features.

the $B$ coefficient have the same issue as covariance, sensitivity to the univariate variance

Here’s a reminder on salient points for each feature ranking metric,

Feature Ranking Metric	Relevancy	Redundancy	Non-linearity	Heteroscedasticity
Correlation	Yes	No	No	No
Rank correlation	Yes	No	Limited Monotonic	No
$b$ coefficients	Yes	Yes	No	No
$\beta$ coefficients	Yes	Yes	No	No
Random forest feature importance	Yes	Partially	Yes	Yes
Mutual information	Yes	No	Yes	Yes

Additional points and reminders:

correlation coefficient - measures linear relevance between a predictor and response and is sensitive to outliers.
rank correlation coefficient - measures monotonic relevance and therefore captures some nonlinear, monotonic relationships and is less sensitive to outliers.
$b$ coefficients - measure linear relevance and redundancy through the multivariate linear model, but are sensitive to the scale and variance of the features.
$\beta$ coefficients - measure linear relevance and redundancy through the multivariate linear model and are not sensitive to the scale or variance of the features, since the features are standardized before fitting the multivariate linear model.
Random forest feature importance - can capture nonlinear relationships and interactions, but its treatment of redundancy is indirect and can distribute importance among correlated predictors; it is not sensitive to feature scale or variance.
Mutual information - measures general statistical, nonparametric dependence and therefore can capture nonlinear relationships, but pairwise mutual information does not account for predictor redundancy.

rankings = [correlation[:-1],rank_correlation[:-1],b[0:len(pred)],beta[0:len(pred)],importance_rank,mi] # list of all metrics

ranking_names = ['Correlation','Rank Correlation','$b$ Coefficients',r'$\beta$ Coefficients','Random Forest','Mutual Information']

plot_ranking_summary(rankings, pred, ranking_names)

_images/e76d114d47d5de37c88f205188655e61bf1a27cbb4e3da72a661638b861f8516.png

Given all of these methods, I would rank the variables as:

Porosity - consistently ranks highly, with only the $b$ coefficients disagreeing. However, $b$ coefficients are sensitive to the scale and variance of the features, while the $\beta$ coefficients are more reliable for comparing features on different scales and agree with porosity as the top-ranked feature.
Permeability - consistently ranks as the second or third feature, with a second rank for metrics based on relevance and a third rank for metrics that account for both relevance and redundancy, likely due to its strong relationship with porosity.
VR - ranks highly with metrics that account for feature redundancy and could provide useful new information to support our prediction model.

I have assigned these ranks by observing the general trend across these metrics. Of course, there are,

trade-offs between feature redundancy and feature relevance, and between model complexity and model interpretability
we should not neglect expert knowledge. If additional information is known about physical processes, causation, and the reliability and availability of variables, this should be integrated when assigning ranks.

There is an empirical approach to sorting predictor features known as recursive feature elimination.

Recursive Feature Elimination#

Recursive Feature Elimination (RFE) works by recursively removing features and building a model with the remaining features. The steps are:

build a model and use the model-based feature ranking to rank the features, e.g., feature importance or the $\beta$ coefficient
remove the least important feature, also referred to as pruning the least important feature
repeat until there is only one feature remaining

The feature ranks are,

first feature removed is the worst
second feature removed is the second worst

$\vdots$

second last feature remaining in the second best
the last feature remaining is the best

The ‘scikit-learn’ Python package has a recursive feature elimination function, RFE in the feature_selection module,

the function accepts any prediction model with information about feature importance

Below, multilinear regression with $\beta$ coefficients is applied in the recursive feature elimination function.

rfe_linear = RFE(LinearRegression(), n_features_to_select=1, verbose=0)
rfe_linear = rfe_linear.fit(dfS[pred].values, np.ravel(dfS[resp]))

rfe_order = np.argsort(rfe_linear.ranking_)                   # RFE ranking: 1 = best, larger numbers = removed earlier

print('Recursive Feature Elimination: Multilinear Regression')

for i, index in enumerate(rfe_order):
    print('Rank #' + str(i + 1) + ' ' + pred[index])

Recursive Feature Elimination: Multilinear Regression
Rank #1 Por
Rank #2 VR
Rank #3 AI
Rank #4 TOC
Rank #5 Perm
Rank #6 Brittle

The advantages with the recursive elimination method:

the actual model can be used in assessing feature ranks
the ranking is based on accuracy of the estimate

but this method is sensitive to:

choice of model
training dataset

The feature ranks are quite different from our previous methods. Many have moved from the previous assessment. Perhaps we should use a more flexible modeling method.

Warning: recursive feature elimination is model-based and the goodness of the feature ranks depends on a good model. It would still be useful to check the model with a cross validation plot.

Let’s check the model prediction performance for the training data over the recursive feature elimination steps.

n_features = len(pred)

plt.figure(figsize=(15, 12))

for i in range(n_features):                                   # loop over the recursive feature elimination steps
    selected_indices = rfe_order[:n_features-i]               # features retained at this RFE step
    selected_features = [pred[j] for j in selected_indices]

    model = LinearRegression()                                # fit model with the remaining features
    model.fit(dfS[selected_features].values, np.ravel(dfS[resp]))

    y_hat = model.predict(dfS[selected_features].values)      # predictions

    plt.subplot(3, 3, i + 1)
    plot_prediction_check(dfS[resp].values,y_hat,threshold=10.0,ymin=dfS[resp].min(),ymax=dfS[resp].max(),
        xlabel='Actual',ylabel='Estimated',title=f'Training Prediction Check - {len(selected_features)} Features')

    feature_text = ', '.join(selected_features)               # list the features in the model
    plt.text(0.03, 0.97,feature_text,transform=plt.gca().transAxes,ha='left',va='top', fontsize=7,
        bbox=dict(facecolor='white', edgecolor='black', alpha=0.85))

plt.tight_layout()

_images/f8609dad3f14a2bf2cb0b3ef146d9b3c14f0b73e1511b7ca8c69281a4fab323f.png

From checking the model performance at training data, we can observe:

the models are ok, without evidence of underfit nor extreme error
a feature can be ranked as important by RFE, but removing it may have little effect on model performance.

Conversely:

a feature that is removed early may still be part of a model that predicts reasonably well because its information is redundant with other features.

Some more details about hyperparameter tuning for recursive feature elimination:

Repeated hyperparameter tuning — Hyperparameters could be tuned independently for each iteration of the feature set. This is often unnecessary and answers a different question: “What is the best model for each feature subset?” This may distract from the original question: “What are the feature ranks for a fixed model configuration?”
Tune once, then hold hyperparameters constant — Tune the model over a reasonable set of hyperparameters and then hold the selected hyperparameters constant across all iterations of feature elimination. This provides a more consistent basis for comparing feature subsets and their resulting feature rankings.

To demonstration this, we demonstrate recursive feature elimination with random forest,

let’s first perform hyperparameter tuning and then we will hold all hyperparameters constant over all the recursive feature elimination iterations.
hyperparameter tuning is performed by grid search, with k-fold cross validation

rf_tuning = RandomForestRegressor(random_state=42)            # instantiate random forest model

params = {                                                    # hyperparameter values for grid search tuning
    'max_leaf_nodes': [5, 10, 20, 40],
    'n_estimators': [100, 300, 500],
    'max_features': [0.5, 0.75, 1.0]
}

grid = GridSearchCV(rf_tuning, params, cv=5, scoring='neg_mean_squared_error', n_jobs=-1) # grid search hyperparameter tuning
grid.fit(dfS[pred].values, np.ravel(dfS[resp]))

print('Random Forest: tuned hyperparameters: '); print(grid.best_params_)

Random Forest: tuned hyperparameters: 
{'max_features': 1.0, 'max_leaf_nodes': 40, 'n_estimators': 500}

Now we pass these tuned hyperparameters to the random forest model within the recursive feature elimination function,

rfe_rf = RFE(RandomForestRegressor(**grid.best_params_, random_state=42), n_features_to_select=1, verbose=0) # instantiate tuned model
rfe_rf = rfe_linear.fit(dfS[pred].values, np.ravel(dfS[resp])) # fit tuned model to training data

rfe_rf_order = np.argsort(rfe_rf.ranking_)                    # RFE ranking: 1 = best, larger numbers = removed earlier

print('Recursive Feature Elimination: Random Forest')

for i, index in enumerate(rfe_rf_order):
    print('Rank #' + str(i + 1) + ' ' + pred[index])

Recursive Feature Elimination: Random Forest
Rank #1 Por
Rank #2 VR
Rank #3 AI
Rank #4 TOC
Rank #5 Perm
Rank #6 Brittle

Once again, check the model predictions over the training data.

note random forest requires a minimum of 2 predictor features, so we stop with the 2 predictor feature model and then assign first and second place based on feature importance with that model

n_features = len(pred)
min_features = 2

plt.figure(figsize=(15, 12))

for i in range(n_features - min_features + 1):                # loop from all features down to 2
    n_selected = n_features - i
    selected_indices = rfe_order[:n_selected]                 # features retained at this RFE step
    selected_features = [pred[j] for j in selected_indices]

    rfe_rf = RandomForestRegressor(**grid.best_params_, random_state=42) 
    rfe_rf.fit(dfS[selected_features].values, np.ravel(dfS[resp])) # fit Random Forest using the features retained at this RFE step

    y_hat = rfe_rf.predict(dfS[selected_features].values)     # predictions

    plt.subplot(3, 3, i + 1)
    plot_prediction_check(dfS[resp].values,y_hat,threshold=10.0,ymin=dfS[resp].min(),ymax=dfS[resp].max(),xlabel='Actual',
      ylabel='Estimated',title=f'Training Prediction Check - {len(selected_features)} Features'
    )

    feature_text = ', '.join(selected_features)
    plt.text(0.03, 0.97, feature_text,transform=plt.gca().transAxes,ha='left',va='top',fontsize=7,
       bbox=dict(facecolor='white', edgecolor='black', alpha=0.85))

plt.tight_layout()

_images/caccbf373831d92073aa62bb13c82840a0ba8e322006c1de8059858ba3d52943.png

From checking the model performance on the training data, we can observe:

the models reproduce the training data reasonably well, with no obvious evidence of underfitting or extreme residual errors;
a feature can be ranked as important by RFE, but removing it may have little effect on model performance.

Warning: Since model performance is evaluated only on the training data, these results do not assess the model’s ability to generalize to new data, and the model could still be overfit. Random Forest models are generally relatively resistant to overfitting, particularly when compared with individual decision trees. Overfitting remains possible, but the risk is reduced by initial k-fold cross-validation for hyperparameter tuning using all features. The selected hyperparameters are then held constant throughout the RFE process.

This is a justifiable workflow because it separates two important aspects of the modeling process:

Hyperparameter tuning: k-fold cross-validation is used once, using all features, to select a reasonable Random Forest configuration.
Recursive feature elimination: the selected hyperparameters are held constant while features are recursively removed. Training-data prediction checks are then used to demonstrate how model performance changes as features are eliminated.

Shapley Values#

Shapley value for feature ranking provides a very useful and flexible, model-based local and global feature importance by learning contribution of each feature to the prediction. Let’s start with the concepts of local and global feature importance.

All our previous feature ranking methods are,

global feature importance - a single metric that represents the importance of the feature, based on a statistic or model predictions, that is assumed representative for all future predictions of the model.

Shapley values also provide,

local feature importance - metrics that represent the importance of the features, for each specific model prediction. With local feature importance, the feature rankings may change for each prediction.

To demonstrate local and global feature importance, consider this visualization of a trained and tuned prediction model given predictor features porosity and brittleness to predict production.

Given a model we can calculate global feature importance over various testing, withheld data (left) or local feature importance for one testing datum (right).

Global feature importance summarizes the impact of porosity and brittleness over all predictions (white circles on the left), while local feature importance provides the impact of porosity and brittleness for an single prediction (white circle on the right).

This ability to calculate the impact of each predictor feature for each prediction of the repsonse feature, results in an important second use case for Shapely values,

Explainable Machine Learning: complicated models are often required but have low interpretability.

Two choices to improve model interpretability:

reduce the complexity of the models, but may also reduce model accuracy
develop improved, agnostic (for any model) model diagnostics

Shapely values are model agnostic, as shown in the figure below, we put any model in the box and simpley pass a variety of prediction cases, known as background, and observe the model predictions to learn the behavoir of the model.

Shapley is a general method for feature importance, based on the model as a black box (model agnostic).

Now let’s backup and talk about the original application of Shapley values, from co-operative game theory. We will talk about games and players and not features and models, but don’t worry - we will get back there.

Shapley Values in Cooperative Game Theory#

In cooperative game theory, we must calculate the allocation of earnings over the game participants. For example,

player 1 and 2 play a game together and earn $1,000,000.00. How do we divide these earnings?

This is accomplished through the concept of,

marginal contribution – the impact of one player on the result in a single case. For example, we can calculate,

Marginal Contribution Player 1 = Outcome Player 1 – Outcome No Players
Marginal Contribution Player 1 = Outcome Player 1 & 2 – Outcome Player 2 Only

Both of these are marginal contributions, but with different coalition sizes. We could continue up to the number of players, for example,

Marginal Contribution Player 1 = Outcome Player 1, 2 and 3 – Outcome Player 2 and 3
Marginal Contribution Player 1 = Outcome Player 1, 2, 3 and 4 – Outcome Player 2, 3 and 4

But we are missing a lot of other combinations over which we can calculate marginal contribution. To generalize, let’s introduce the idea of coalition size, where the coalition size indicates the number of other players present when Player 1 is added:

First level - Player 1 subtract no one plays (zero earnings)

Second level - Player 1 and another player subtract that other player by themselves

Third level - Player 1 and two other players subtract those two other players together

Fourth level - Player 1 and three other players subtract those three other players together

To clarify, here is a table, with all marginal contributions for four players, $P_1,P_2,P_3,P_4$,

Coalition Size	Coalition $S$	Add $P_1$	Add $P_2$	Add $P_3$	Add $P_4$
0	$\emptyset$	$v(P_1)-v(\emptyset)$	$v(P_2)-v(\emptyset)$	$v(P_3)-v(\emptyset)$	$v(P_4)-v(\emptyset)$
1	$P_1$	—	$v(P_1,P_2)-v(P_1)$	$v(P_1,P_3)-v(P_1)$	$v(P_1,P_4)-v(P_1)$
1	$P_2$	$v(P_1,P_2)-v(P_2)$	—	$v(P_2,P_3)-v(P_2)$	$v(P_2,P_4)-v(P_2)$
1	$P_3$	$v(P_1,P_3)-v(P_3)$	$v(P_2,P_3)-v(P_3)$	—	$v(P_3,P_4)-v(P_4)$
1	$P_4$	$v(P_1,P_4)-v(P_4)$	$v(P_2,P_4)-v(P_4)$	$v(P_3,P_4)-v(P_4)$	—
2	$P_1,P_2$	—	—	$v(P_1,P_2,P_3)-v(P_1,P_2)$	$v(P_1,P_2,P_4)-v(P_1,P_2)$
2	$P_1,P_3$	—	$v(P_1,P_2,P_3)-v(P_1,P_3)$	—	$v(P_1,P_3,P_4)-v(P_1,P_3)$
2	$P_1,P_4$	—	$v(P_1,P_2,P_4)-v(P_1,P_4)$	$v(P_1,P_3,P_4)-v(P_1,P_4)$	—
2	$P_2,P_3$	$v(P_1,P_2,P_3)-v(P_2,P_3)$	—	—	$v(P_2,P_3,P_4)-v(P_2,P_3)$
2	$P_2,P_4$	$v(P_1,P_2,P_4)-v(P_2,P_4)$	—	$v(P_2,P_3,P_4)-v(P_2,P_4)$	—
2	$P_3,P_4$	$v(P_1,P_3,P_4)-v(P_3,P_4)$	$v(P_2,P_3,P_4)-v(P_3,P_4)$	—	—
3	$P_1,P_2,P_3$	—	—	—	$v(P_1,P_2,P_3,P_4)-v(P_1,P_2,P_3)$
3	$P_1,P_2,P_4$	—	—	$v(P_1,P_2,P_3,P_4)-v(P_1,P_2,P_4)$	—
3	$P_1,P_3,P_4$	—	$v(P_1,P_2,P_3,P_4)-v(P_1,P_3,P_4)$	—	—
3	$P_2,P_3,P_4$	$v(P_1,P_2,P_3,P_4)-v(P_2,P_3,P_4)$	—	—	—

Recall, all of these marginal contributions are calculated from the same game. How do we summarize all of these marginal contributions into a single value for each player that will inform splitting the earnings?

we do a weighted average over the coalition sizes
since there are 4 coalition sizes, each coalition size gets $\frac{1}{4}$ weight
since coalition sizes 0 and 3 have only one case they both get the coalition size weight $\frac{1}{4}$
since coalition sizes 1 and 2 have 3 cases each, they share the coalition size weight $\frac{1}{4} \times \frac{1}{3} = \frac{1}{12}$

Coalition Size	Number of other players	Number of coalitions	Weight per coalition
0	0	1	$1/4$
1	1	3	$1/12$
2	2	3	$1/12$
3	3	1	$1/4$

Now we can formulate the equation for the Shapley value of $P_1$ of the four players $P_1,P_2,P_3,P_4$, as a weighted average with each coalition size on a separate line,

\[\begin{split} \begin{aligned} \phi_1 =\;& \frac{1}{4}\left[v(P_1)-v(\emptyset)\right] \\[4pt] &+\frac{1}{12}\left[ v(P_1,P_2)-v(P_2) + v(P_1,P_3)-v(P_3) + v(P_1,P_4)-v(P_4) \right] \\[4pt] &+\frac{1}{12}\left[ v(P_1,P_2,P_3)-v(P_2,P_3) + v(P_1,P_2,P_4)-v(P_2,P_4) + v(P_1,P_3,P_4)-v(P_3,P_4) \right] \\[4pt] &+\frac{1}{4}\left[ v(P_1,P_2,P_3,P_4)-v(P_2,P_3,P_4) \right]. \end{aligned} \end{split}\]

This illustrates an important concept about Shapley values,

it is not simply $v(P_1)$, nor a single marginal contribution, for example, $v(P_1,P_2)-v(P_2)$, but it measures how much the model value changes when $P_1$ is added across all possible combinations of other features, from a coalition size of 0 to a coalition size of 3.

While illustrative and comprehensive, the equation notation is unwieldy,

we can introduce the notation for marginal contribution, $\Delta_i(S)$, for the marginal contribution of player $P_i$ when added to coalition $S$:

\[ \Delta_i(S) = v(S \cup \{P_i\}) - v(S) \]

Now we can generalize the above equation for any case as,

\[ \phi_i = \sum_{S\subseteq N\setminus\{i\}} \frac{|S|!(n-|S|-1)!}{n!} \left[v(S\cup\{P_i\})-v(S) \right]. \]

where:

$\phi_i$ is the Shapley value for player $P_i$, representing the average marginal contribution of $P_i$ across all possible coalitions of the other players.
$P_i$ is the player or feature being evaluated.
$N$ is the set of all players or features.
$S$ is a subset (coalition) of $N$ that does not contain $P_i$.
$|S|$ is the number of players or features in coalition $S$.
$n$ is the total number of players or features, such that $n=|N|$.
$v(S)$ is the value of the model or coalition $S$.
$v(S\cup{P_i})$ is the value of the model or coalition after adding player $P_i$ to $S$.
$v(S\cup{P_i})-v(S)$ is the marginal contribution of player $P_i$ when added to coalition $S$.
$\frac{|S|!(n-|S|-1)!}{n!}$ is the Shapley weighting factor for coalition $S$, which accounts for the number of possible player orderings that place $P_i$ immediately after the players in $S$.
$\sum_{S\subseteq N\setminus{i}}$ indicates that the marginal contribution is summed over all possible coalitions $S$ that can be formed from the players

If the notation and factorials are confusing remember,

Shapley value is just the weighted average over orders of marginal contributions, earnings with subtract without a specific player.

Also, here’s an illutration of Shapley value for The Beatles to allocate records sold to John Lennon, given other band members Paul McCartney, George Harrison, and Ringo Starr.

Shapley value illustration for The Beatles for records sold.

I like to use this example in class,

helps explain the concept of Shapley values
I am surprised by how many of my students still like The Beatles.

Simple Shapley Value Examples#

Here’s three very simple demonstrations for Shapley values.

Shapley Additive Partnership#

Two people work together, person 1, $P_1$, and person 2, $P_2$, and earn $125,000. How should the profits be split?

\[ f(\mathrm{P}_1)=\$50,000,\qquad f(\mathrm{P}_2)=\$75,000,\qquad f(\mathrm{person}_1,\mathrm{person}_2)=\$125,000 \]

We average over the 0 and 1 coalition sizes,

\[ \phi_{P_1} = \frac{1}{2}\left[ f(\mathrm{P}_1) - \emptyset \right] + \frac{1}{2} \left[f(\mathrm{P}_1,\mathrm{P}_2)-f(\mathrm{P}_2)\right] \]

\[ \phi_{P_2} = \frac{1}{2}\left[ f(\mathrm{P}_2) - \emptyset \right] + \frac{1}{2} \left[f(\mathrm{P}_1,\mathrm{P}_2)-f(\mathrm{P}_1)\right] \]

Substituting into this equation we get,

\[ \phi_{P_1} = \frac{1}{2}\left[ \$50,000 - 0 \right] + \frac{1}{2} \left[ \$125,000 - \$75,000 \right] = \$50,000 \]

\[ \phi_{P_2} = \frac{1}{2}\left[ \$75,000 - 0 \right] + \frac{1}{2} \left[ \$125,000 - \$50,000 \right] = \$75,000 \]

Note, the two allocations sum to the total earnings,

\[ $50,000+$75,000=$125,000. \]

this is an important property of Shapley values, known as efficiency.

Shapley Synergistic Partnership#

The numbers above may seem “cooked” as the earnings for both players 1 and 2 were equal to the earning of each player on their own. Here’s is the same case, but with synergy, i.e., the two players together produce more earnings than the sum of their individual efforts.

Now how should the profits be split?

\[ f(\mathrm{P}_1)=\$50,000,\qquad f(\mathrm{P}_2)=\$70,000,\qquad f(\mathrm{person}_1,\mathrm{person}_2)=\$150,000 \]

We average over the 0 and 1 coalition sizes,

\[ \phi_{P_1} = \frac{1}{2}\left[ f(\mathrm{P}_1) - \emptyset \right] + \frac{1}{2} \left[f(\mathrm{P}_1,\mathrm{P}_2)-f(\mathrm{P}_2)\right] \]

\[ \phi_{P_2} = \frac{1}{2}\left[ f(\mathrm{P}_2) - \emptyset \right] + \frac{1}{2} \left[f(\mathrm{P}_1,\mathrm{P}_2)-f(\mathrm{P}_1)\right] \]

Substituting into this equation we get,

\[ \phi_{P_1} = \frac{1}{2}\left[ \$50,000 - 0 \right] + \frac{1}{2} \left[ \$150,000 - \$70,000 \right] = \$65,000 \]

\[ \phi_{P_2} = \frac{1}{2}\left[ \$70,000 - 0 \right] + \frac{1}{2} \left[ \$150,000 - \$50,000 \right] = \$85,000 \]

Note, the two allocations still sum to the total earnings and the Shapley values demonstrate efficiency again,

\[ $65,000+$85,000=$150,000. \]

Shapley Dysnergistic Partnership#

Now, for completeness, consider the partnership that should have never happenned, a case with dysnergy, i.e., the two players together produce less earnings than the sum of their individual efforts, i.e., they do better working alone.

Now how should the profits be split?

\[ f(\mathrm{P}_1)=\$50,000,\qquad f(\mathrm{P}_2)=\$70,000,\qquad f(\mathrm{person}_1,\mathrm{person}_2)=\$90,000 \]

We average over the 0 and 1 coalition sizes,

\[ \phi_{P_1} = \frac{1}{2}\left[ f(\mathrm{P}_1) - \emptyset \right] + \frac{1}{2} \left[f(\mathrm{P}_1,\mathrm{P}_2)-f(\mathrm{P}_2)\right] \]

\[ \phi_{P_2} = \frac{1}{2}\left[ f(\mathrm{P}_2) - \emptyset \right] + \frac{1}{2} \left[f(\mathrm{P}_1,\mathrm{P}_2)-f(\mathrm{P}_1)\right] \]

Substituting into this equation we get,

\[ \phi_{P_1} = \frac{1}{2}\left[ \$50,000 - 0 \right] + \frac{1}{2} \left[ \$90,000 - \$70,000 \right] = \$35,000 \]

\[ \phi_{P_2} = \frac{1}{2}\left[ \$70,000 - 0 \right] + \frac{1}{2} \left[ \$90,000 - \$50,000 \right] = \$55,000 \]

Note, the two allocations still sum to the total earnings and the Shapley values demonstrate efficiency again,

\[ $35,000+$55,000=$90,000. \]

Now we are ready to return to feature ranking.

Shapley Values for Feature Ranking#

We have demonstrated Shapley values for cooperative game theory the contribution of each player through summarization over marginal contributions.

Now for machine learning applications, we change,

player $\rightarrow$ feature, $X_i$
earnings $rightarrow$ model prediction, \hat{y} = f(X)

if we return to our previous example with 4 players, we can update our table for 4 predictor features, $X_1, X_2, X_3, X_4$ and denote the model predictions with $f(\cdot)$.

Coalition Size	Coalition $S$	Add $X_1$	Add $X_2$	Add $X_3$	Add $X_4$
0	$\emptyset$	$f(X_1)-f(\emptyset)$	$f(X_2)-f(\emptyset)$	$f(X_3)-f(\emptyset)$	$f(X_4)-f(\emptyset)$
1	$X_1$	—	$f(X_1,X_2)-f(X_1)$	$f(X_1,X_3)-f(X_1)$	$f(X_1,X_4)-f(X_1)$
1	$X_2$	$f(X_1,X_2)-f(X_2)$	—	$f(X_2,X_3)-f(X_2)$	$f(X_2,X_4)-f(X_2)$
1	$X_3$	$f(X_1,X_3)-f(X_3)$	$f(X_2,X_3)-f(X_3)$	—	$f(X_3,X_4)-f(X_3)$
1	$X_4$	$f(X_1,X_4)-f(X_4)$	$f(X_2,X_4)-f(X_4)$	$f(X_3,X_4)-f(X_4)$	—
2	$X_1,X_2$	—	—	$f(X_1,X_2,X_3)-f(X_1,X_2)$	$f(X_1,X_2,X_4)-f(X_1,X_2)$
2	$X_1,X_3$	—	$f(X_1,X_2,X_3)-f(X_1,X_3)$	—	$f(X_1,X_3,X_4)-f(X_1,X_3)$
2	$X_1,X_4$	—	$f(X_1,X_2,X_4)-f(X_1,X_4)$	$f(X_1,X_3,X_4)-f(X_1,X_4)$	—
2	$X_2,X_3$	$f(X_1,X_2,X_3)-f(X_2,X_3)$	—	—	$f(X_2,X_3,X_4)-f(X_2,X_3)$
2	$X_2,X_4$	$f(X_1,X_2,X_4)-f(X_2,X_4)$	—	$f(X_2,X_3,X_4)-f(X_2,X_4)$	—
2	$X_3,X_4$	$f(X_1,X_3,X_4)-f(X_3,X_4)$	$f(X_2,X_3,X_4)-f(X_3,X_4)$	—	—
3	$X_1,X_2,X_3$	—	—	—	$f(X_1,X_2,X_3,X_4)-f(X_1,X_2,X_3)$
3	$X_1,X_2,X_4$	—	—	$f(X_1,X_2,X_3,X_4)-f(X_1,X_2,X_4)$	—
3	$X_1,X_3,X_4$	—	$f(X_1,X_2,X_3,X_4)-f(X_1,X_3,X_4)$	—	—
3	$X_2,X_3,X_4$	$f(X_1,X_2,X_3,X_4)-f(X_2,X_3,X_4)$	—	—	—

Also, we can update our general equation for calculating Shapley value for feature, $X_i$, as $\phi_{X_i}$,

\[ \phi_{X_i} = \sum_{S\subseteq N\setminus\{i\}} \frac{|S|!(n-|S|-1)!}{n!} \left[f(S\cup\{X_i\})-f(S) \right]. \]

where:

$\phi_i$ is the Shapley value for feature $X_i$, representing the average marginal contribution of $X_i$ to the model prediction across all possible coalitions of the other features.
$X_i$ is the feature being evaluated.
$N$ is the set of all features in the model.
$S$ is a subset (coalition) of $N$ that does not contain $X_i$.
$|S|$ is the number of features in coalition $S$.
$n$ is the total number of features, such that $n=|N|$.
$f(S)$ is the model prediction using the features in coalition $S$.
$f(S\cup{X_i})$ is the model prediction after adding feature $X_i$ to coalition $S$.
$f(S\cup{X_i})-f(S)$ is the marginal contribution of feature $X_i$ when added to coalition $S$.
$\frac{|S|!(n-|S|-1)!}{n!}$ is the Shapley weighting factor for coalition $S$, which accounts for the number of possible feature orderings that place $X_i$ immediately after the features in $S$.
$\sum_{S\subseteq N\setminus{i}}$ indicates that the marginal contribution is summed over all possible coalitions $S$ that can be formed from the features other than $X_i$.

Now we have Shapley value for feature inportance and explainable machine learning artificial intelligence models. There are just a couple of additional details that we need to cover,

removing features from the model
summarizing local Shapley value to calculate a global Shapley value as a global measure of feature inportance

Removing Features from the Model#

To calculate Shapley values for feature ranking and explainable AI, we need the model predictions for all possible coalitions of predictor features. For example, for a model with four predictor features, $X_1$, $X_2$, $X_3$, and $X_4$, the required model predictions include:

Coalition size 0: $f(\emptyset)$
Coalition size 1: $f(X_1)$, $f(X_2)$, $f(X_3)$, $f(X_4)$
Coalition size 2: $f(X_1,X_2)$, $f(X_1,X_3)$, $f(X_1,X_4)$, $f(X_2,X_3)$, $f(X_2,X_4)$, $f(X_3,X_4)$
Coalition size 3: $f(X_1,X_2,X_3)$, $f(X_1,X_2,X_4)$, $f(X_1,X_3,X_4)$, $f(X_2,X_3,X_4)$
Coalition size 4: $f(X_1,X_2,X_3,X_4)$

For four features, this gives $2^4=16$ possible coalitions, including the empty coalition. The Shapley value for each feature summarizes its marginal contribution across these possible coalitions, with each marginal contribution assigned an appropriate Shapley weighting factor.

But there is a problem:

we only have one model, $f(X_1,X_2,X_3,X_4)$.
in general, we cannot ask our model for a prediction while simply omitting one feature, such as $f(X_1,X_2,X_3)$, or two features, such as $f(X_1,X_2)$, or three features, such as $f(X_1)$.

Naive Shapley Approach for Removing Features#

The naïve approach to solving this problem is to train the full combination of models required for all possible feature coalitions.

However, we do not want to do this if our goal is feature importance or explainable AI for diagnosing our original model, $f(X_1,X_2,X_3,X_4)$.

In general, this approach is not recommended.

Shapley with Feature Imputation for Removed Features#

A more practical method is to apply feature imputation to replace a removed feature with a representative value. There are a variety of approaches, similar to those used for feature imputation:

impute the expected value of the removed feature,

\[ f(X_1, X_2, X_3, X_4) = f(X_1, X_2, X_3, X_4 = E[X_4]) \]

impute the median value of the removed feature,

\[ f(X_1, X_2, X_3, X_4) = f(X_1, X_2, X_3, X_4 = P50_{X_4}) \]

Shapley with Tree-based Prediction Models#

Tree-based prediction models provide a practical approach to removing features without retraining the model. When a removed feature is encountered at a decision node, the prediction can be calculated by averaging over the branches, weighted by the number of training samples following each branch.

Here’s a simple decision tree and a single prediction case with feature $X_2$ removed and $X_1 = 25$.

Removing feature X2 from a simple decision tree for a single prediction.

Since $X_2$ is not encountered along any of the tree branches for this prediction, no action is required.

\[ f(X_1 = 25) = 20 \]

Now we remove feature $X_1$ for a single prediction case with $X_2=60$.

Removing feature X1 from a simple decision tree for a single prediction.

Now we encounter the removed feature $X_1$ at the first decision node. Since $X_1$ is unavailable, we cannot determine which branch to follow. Instead, we go down both branches and weight each path by the number of training samples following that path.

the left path leads to 2 leaf nodes with 15 and 45 samples
the right path leads to 2 leaf nodes with 30 and 10 samples

We weight the paths as:

left path as $\frac{15+45}{15+45+30+10}=\frac{60}{100}$
right path as $\frac{30+10}{15+45+30+10}=\frac{40}{100}$

At the second-level decision node on the left path from the first level, we encounter $X_1$ again. Since $X_1$ is still removed, we again go down both paths and weight each path to the leaf nodes as:

left path as $\frac{15}{15+45}=\frac{15}{60}$
right path as $\frac{45}{15+45}=\frac{45}{60}$

At the second-level decision node on the right path from the first level, we encounter $X_2$. Since $X_2$ is available, we can follow the appropriate branch based on $X_2=60$ and proceed to the leaf node.

Therefore, the prediction with $X_1$ removed is:

\[ f(X_2 = 60) = \frac{60}{100} \left[ \frac{15}{60}\times20 + \frac{45}{60}\times70 \right] + \frac{40}{100}[130] = 86.5 \]

Global Shapley by Summarizing Local Shapley#

A single measure of feature importance for each predictor can be calculated by summarizing local feature importance over a diverse set of predictions, called the background.

See the plot on the right below, which shows the Shapley value for each feature across the background predictions.

the mean absolute Shapley value for each feature provides a global measure of Shapley-based feature importance, summarizing the magnitude of the impact of each predictor feature on the model predictions across the background.

The plot on the left is the sorted bar chart of the mean absolute Shapley values, from the most important feature to the least important feature.

Shapley-based feature ranking by mean absolute Shapley value.

Shapley Value for Feature Ranking Demonstration#

Let’s take a random subset of the data as background values to evaluate our model.

we subset the data for faster calculation
we should ensure efficient coverage of the predictor feature space

Since Shapley values are model-based, we must first build a model.

Build a Random Forest Model#

Let’s start with a good random forest model and examine both local and global Shapley values.

Three different sets of model hyperparameters may be selected to demonstrate the impact of model performance on Shapley-based feature ranking:

underfit model
overfit model
tuned model

seed = 73093                                                  # set the random forest hyperparameters

selected_model = 3                                            # select: 1 - underfit model, 2 - overfit model, 3 - tuned model

if selected_model == 1:                                       # underfit random forest
    max_leaf_nodes = 2                                          
    num_tree = 10
    max_features = 2
elif selected_model == 2:                                     # overfit random forest
    max_leaf_nodes = 50                                           
    num_tree = 1
    max_features = 6
elif selected_model == 3:                                     # tuned random forest - see above for tuning code
    max_leaf_nodes = 40
    num_tree = 500
    max_features = 1

rfr = RandomForestRegressor(max_leaf_nodes=max_leaf_nodes, random_state=seed,n_estimators=num_tree, max_features=max_features)
rfr.fit(X = x, y = Y)

Y_hat = predict_train = rfr.predict(x)

MSE = metrics.mean_squared_error(Y,Y_hat)
Var_Explained = metrics.explained_variance_score(Y,Y_hat)
print('Mean Squared Error on Training = ', round(MSE,2),', Variance Explained =', round(Var_Explained,2))

importances = rfr.feature_importances_                        # expected (global) importance over the forest fore each predictor feature
std = np.std([rfr.feature_importances_ for tree in rfr.estimators_],axis=0)
indices = np.argsort(importances)[::-1].tolist()

plt.subplot(121)
plt.scatter(Y,Y_hat,s=None, c='darkorange',marker=None, cmap=None, norm=None, vmin=None, vmax=None, alpha=0.8, linewidths=0.3, edgecolors="black")
plt.title('Random Forest Model'); plt.xlabel('Actual Production (MCFPD)'); plt.ylabel('Estimated Production (MCFPD)')
plt.xlim(0,7000); plt.ylim(0,7000)
plt.arrow(0,0,7000,7000,width=0.02,color='black',head_length=0.0,head_width=0.0)

plt.subplot(122)
plt.title("Feature Importances")
plt.bar([pred[i] for i in indices],rfr.feature_importances_[indices],color="darkorange", alpha = 0.8, edgecolor = 'black', yerr=std[indices], align="center")
plt.ylim(0,1), plt.xlabel('Predictor Features'); plt.ylabel('Feature Importance')
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.2, top=0.8, wspace=0.2, hspace=0.2); plt.show()

Mean Squared Error on Training =  88405.22 , Variance Explained = 0.96

_images/331310169e7b3f9815b08cda9de56191260399ffecb7a95147f55b505923ae3a.png

Calculate Shapley Values#

Let’s select some background data at random to calculate local Shapley values and then summarize them with global Shapley measures.

Background Samples are selected as a random subset from the full dataset. Why not just use all the data as background?

Shapley values can be computationally expensive to calculate, as we need predictions for all combinations of features to calculate the marginal contributions that are summarized as Shapley values.
The background data should be representative, so we want to sample the original data in a manner that provides good coverage of the prediction cases of interest and reduces bias in our feature importance assessment.
Generalization versus specific prediction cases, using all the data as background provides an overall assessment of feature importance. However, we may instead want to carefully select background data to explore specific types of prediction cases.

For simplicity, here we randomly select $n$ data samples as background.

background = shap.sample(x,nsamples=50,random_state=73073) 
model_explainer = shap.TreeExplainer(rfr)
shap_values = model_explainer.shap_values(background)         # global Shapley Measures

Local Shapley Values#

Let’s start by looking at the local Shapley values to demonstrate the concept of efficiency.

first, let’s confirm that the output from the shap function is an $\left[n_{background},m\right]$ NumPy array, where $n_{background}$ is the number of background samples and $m$ is the number of predictor features.

shap_values.shape                                             # shape of the shapley value output, n background samples x m features

(50, 6)

Custom Local Force Plot#

We have the local Shapley values for each prediction in the background cases. Let’s visualize one to demonstrate this concept.

I coded this custom visualization to clearly communicate local Shapley values and the concept of efficiency.
We start at the average of the training response feature and add the local Shapley value for each predictor feature to reach the model prediction.

nback = 7                                                     # select a background sample index (from 0 - 49)

resp_avg  = np.average(Y_hat)
yhat = rfr.predict(background.iloc[[nback]])

current = resp_avg

plt.subplot(111)

plt.plot([current,current],[0,0.3],color='black',lw=2,zorder=1)
plt.plot([current-2,current],[0.2,0.3],color='black',lw=2,zorder=1)
plt.plot([current,current+2],[0.3,0.2],color='black',lw=2,zorder=1)
for i in range(len(pred)+1):
    plt.scatter(current,i+0.5,color='grey',edgecolor='black',zorder=10)
    if i < len(pred):
        if shap_values[nback,i] > 0.0:
            color = 'red'
        else:
            color = 'blue'
        plt.plot([current,current + shap_values[nback,i]],[i+1,i+1],color=color,lw=2,zorder=1)
        plt.plot([current,current],[i+0.6,i+1],color=color,lw=2,zorder=1)
        plt.plot([current + shap_values[nback,i],current + shap_values[nback,i]],[i+1,i+1.3],color=color,lw=2,zorder=1)
        plt.plot([current + shap_values[nback,i]-2,current + shap_values[nback,i]],[i+1.2,i+1.3],color=color,lw=2,zorder=1)
        plt.plot([current + shap_values[nback,i],current + shap_values[nback,i]+2],[i+1.3,i+1.2],color=color,lw=2,zorder=1)
        if shap_values[nback,i] > 0.0:
            plt.annotate('+ ' + str(np.round(shap_values[nback,i],0)),[current + shap_values[nback,i]*0.5,i+1.1],ha='center')
        else:
            plt.annotate('- ' + str(np.round(abs(shap_values[nback,i]),0)),[current + shap_values[nback,i]*0.5,i+1.1],ha='center')
        current = current + shap_values[nback,i]

plt.plot([current,current],[i+0.7,i+1],color='black',lw=2,zorder=1)
plt.plot([current-2,current],[i+0.9,i+1],color='black',lw=2,zorder=1)
plt.plot([current,current+2],[i+1,i+0.9],color='black',lw=2,zorder=1)

plt.plot([resp_avg,resp_avg],[-0.5,len(pred)+1.5],color='black',ls='--',zorder=1)
plt.plot([yhat,yhat],[-0.5,len(pred)+1.5],color='black',ls='--',zorder=1)
plt.annotate('Response Feature, Training Average',[resp_avg-8,1.0],rotation=90.0)
plt.annotate('Model Prediction',[yhat-8,1.0],rotation=90.0)

plt.yticks(ticks=np.arange(len(pred)+2), labels=[r'None / $\overline{y}$'] + pred + [r'$\hat{y}=f(X)$'])
add_grid(); plt.ylim([-0.5,len(pred)+1.5])
plt.xlabel('Production (MCFPD)'); plt.ylabel('Feature'); plt.title('Local Shapley Values, Background Index: ' + str(nback))
plt.subplots_adjust(left=0.0, bottom=0.0, right=2.2, top=0.8, wspace=0.2, hspace=0.2); plt.show()

_images/3bfaa72c54be939fe051bf42e2706736708cec1b5bbe1730bafc66a02c70c3c2.png

I like this custom local force plot because it shows:

for a single prediction
the impact of each feature
how the prediction moves from the global mean of the response feature to the actual prediction

Now, let’s review the plots available from the Shap Python package:

force_plot — visualization of the impact of all features across all background sample predictions
local_force_plot — visualization of the impact of all features on a single prediction

Shapley Force Plot#

We can simultaneously visualize all of the Shapley values for all of the sample data, in the order of the background dataset.

blue indicates a reduction in the predicted production, while red indicates an increase in predicted production

We are visualizing all background sample data at once. Reorder by the original sample ordering and select the $n_{back}$ index to compare a specific prediction to the custom local force plot above.

shap.force_plot(model_explainer.expected_value,shap_values,background,out_names = ['Production'],feature_names=pred,)  

Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.

Local Force Plot#

We pick a specific sample from the background data and visualize its force plot.

We can see the basis of the plot above: the Shapley values for all features given the local set of values for sample $i$, $(x_i)$.

Compare this result to the custom plot that I made above, and you will see that it communicates the same information.

shap.force_plot(model_explainer.expected_value,shap_values[nback],background.iloc[[nback]],show=False,feature_names = pred)

Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.

Appreciation to Xuesong Ma for the suggestion to improve the above local Shapley value content and visualizations.

Global Shapley Values#

Let’s review the global Shapley measures.

sorted bar chart of the arithmetic average of the absolute SHAP value over the background data
sorted plot of the SHAP value over the background data
plot of the SHAP value over the background data as a violin plot

Note: all of these methods use the global average ($E[X_i]$) for each feature to impute the feature value for cases where feature $i$ is not included.

plt.subplot(131)
shap.summary_plot(show=False,feature_names = pred, shap_values = shap_values, features = background, plot_type="bar",color = "darkorange",cmap = plt.cm.inferno)
plt.ylabel('Predictor Features')

plt.subplot(132)
shap.summary_plot(show=False,feature_names = pred, shap_values = shap_values, features = background,cmap = plt.cm.inferno)

plt.subplot(133)
shap.summary_plot(show=False,feature_names = pred, shap_values = shap_values, features = background,plot_type = "violin")

plt.subplots_adjust(left=0.0, bottom=0.0, right=2.2, top=1.2, wspace=0.2, hspace=0.2)
plt.show()

_images/ed1857e816934df0eced6ec51a500c097a7df7c6b050e0a56835d4721538f075.png

The center and right plots show the Shapley values for each feature over all of the randomly selected background samples, while the plot on the left shows the bar chart of the mean absolute Shapley values. We can observe from the global Shapley values that:

Porosity, Permeability, and TOC are the top features.

Comments#

This was a basic treatment of feature ranking. Much more could be done and discussed, I have many more resources. Check out my shared resource inventory and the YouTube lecture links at the start of this chapter with resource links in the videos’ descriptions.

I hope this is helpful,

Michael

About the Author#

Michael Pyrcz is a professor in the Cockrell School of Engineering, and the Jackson School of Geosciences, at The University of Texas at Austin, where he researches and teaches subsurface, spatial data analytics, geostatistics, and machine learning. Michael is also,

the principal investigator of the Energy Analytics freshmen research initiative and a core faculty in the Machine Learn Laboratory in the College of Natural Sciences, The University of Texas at Austin
an associate editor for Computers and Geosciences, and a board member for Mathematical Geosciences, the International Association for Mathematical Geosciences.

Michael has written over 100 peer-reviewed publications, a Python package for spatial data analytics, co-authored a textbook on spatial data analytics, Geostatistical Reservoir Modeling and author of two recently released e-books, Applied Geostatistics in Python: a Hands-on Guide with GeostatsPy and Applied Machine Learning in Python: a Hands-on Guide with Code.

All of Michael’s university lectures are available on his YouTube Channel with links to 100s of Python interactive dashboards and well-documented workflows in over 40 repositories on his GitHub account, to support any interested students and working professionals with evergreen content. To find out more about Michael’s work and shared educational resources visit his Website.

Want to Work Together?#

I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.

Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I’d be happy to drop by and work with you!
Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PI is Professor John Foster)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!
I can be reached at mpyrcz@austin.utexas.edu.

I’m always happy to discuss,

Michael

Michael Pyrcz, Ph.D., P.Eng. Professor, Cockrell School of Engineering and The Jackson School of Geosciences, The University of Texas at Austin

Coalition Size	Coalition \(S\)	Add \(P_1\)	Add \(P_2\)	Add \(P_3\)	Add \(P_4\)
0	\(\emptyset\)	\(v(P_1)-v(\emptyset)\)	\(v(P_2)-v(\emptyset)\)	\(v(P_3)-v(\emptyset)\)	\(v(P_4)-v(\emptyset)\)
1	\(P_1\)	—	\(v(P_1,P_2)-v(P_1)\)	\(v(P_1,P_3)-v(P_1)\)	\(v(P_1,P_4)-v(P_1)\)
1	\(P_2\)	\(v(P_1,P_2)-v(P_2)\)	—	\(v(P_2,P_3)-v(P_2)\)	\(v(P_2,P_4)-v(P_2)\)
1	\(P_3\)	\(v(P_1,P_3)-v(P_3)\)	\(v(P_2,P_3)-v(P_3)\)	—	\(v(P_3,P_4)-v(P_4)\)
1	\(P_4\)	\(v(P_1,P_4)-v(P_4)\)	\(v(P_2,P_4)-v(P_4)\)	\(v(P_3,P_4)-v(P_4)\)	—
2	\(P_1,P_2\)	—	—	\(v(P_1,P_2,P_3)-v(P_1,P_2)\)	\(v(P_1,P_2,P_4)-v(P_1,P_2)\)
2	\(P_1,P_3\)	—	\(v(P_1,P_2,P_3)-v(P_1,P_3)\)	—	\(v(P_1,P_3,P_4)-v(P_1,P_3)\)
2	\(P_1,P_4\)	—	\(v(P_1,P_2,P_4)-v(P_1,P_4)\)	\(v(P_1,P_3,P_4)-v(P_1,P_4)\)	—
2	\(P_2,P_3\)	\(v(P_1,P_2,P_3)-v(P_2,P_3)\)	—	—	\(v(P_2,P_3,P_4)-v(P_2,P_3)\)
2	\(P_2,P_4\)	\(v(P_1,P_2,P_4)-v(P_2,P_4)\)	—	\(v(P_2,P_3,P_4)-v(P_2,P_4)\)	—
2	\(P_3,P_4\)	\(v(P_1,P_3,P_4)-v(P_3,P_4)\)	\(v(P_2,P_3,P_4)-v(P_3,P_4)\)	—	—
3	\(P_1,P_2,P_3\)	—	—	—	\(v(P_1,P_2,P_3,P_4)-v(P_1,P_2,P_3)\)
3	\(P_1,P_2,P_4\)	—	—	\(v(P_1,P_2,P_3,P_4)-v(P_1,P_2,P_4)\)	—
3	\(P_1,P_3,P_4\)	—	\(v(P_1,P_2,P_3,P_4)-v(P_1,P_3,P_4)\)	—	—
3	\(P_2,P_3,P_4\)	\(v(P_1,P_2,P_3,P_4)-v(P_2,P_3,P_4)\)	—	—	—

Coalition Size	Coalition \(S\)	Add \(X_1\)	Add \(X_2\)	Add \(X_3\)	Add \(X_4\)
0	\(\emptyset\)	\(f(X_1)-f(\emptyset)\)	\(f(X_2)-f(\emptyset)\)	\(f(X_3)-f(\emptyset)\)	\(f(X_4)-f(\emptyset)\)
1	\(X_1\)	—	\(f(X_1,X_2)-f(X_1)\)	\(f(X_1,X_3)-f(X_1)\)	\(f(X_1,X_4)-f(X_1)\)
1	\(X_2\)	\(f(X_1,X_2)-f(X_2)\)	—	\(f(X_2,X_3)-f(X_2)\)	\(f(X_2,X_4)-f(X_2)\)
1	\(X_3\)	\(f(X_1,X_3)-f(X_3)\)	\(f(X_2,X_3)-f(X_3)\)	—	\(f(X_3,X_4)-f(X_3)\)
1	\(X_4\)	\(f(X_1,X_4)-f(X_4)\)	\(f(X_2,X_4)-f(X_4)\)	\(f(X_3,X_4)-f(X_4)\)	—
2	\(X_1,X_2\)	—	—	\(f(X_1,X_2,X_3)-f(X_1,X_2)\)	\(f(X_1,X_2,X_4)-f(X_1,X_2)\)
2	\(X_1,X_3\)	—	\(f(X_1,X_2,X_3)-f(X_1,X_3)\)	—	\(f(X_1,X_3,X_4)-f(X_1,X_3)\)
2	\(X_1,X_4\)	—	\(f(X_1,X_2,X_4)-f(X_1,X_4)\)	\(f(X_1,X_3,X_4)-f(X_1,X_4)\)	—
2	\(X_2,X_3\)	\(f(X_1,X_2,X_3)-f(X_2,X_3)\)	—	—	\(f(X_2,X_3,X_4)-f(X_2,X_3)\)
2	\(X_2,X_4\)	\(f(X_1,X_2,X_4)-f(X_2,X_4)\)	—	\(f(X_2,X_3,X_4)-f(X_2,X_4)\)	—
2	\(X_3,X_4\)	\(f(X_1,X_3,X_4)-f(X_3,X_4)\)	\(f(X_2,X_3,X_4)-f(X_3,X_4)\)	—	—
3	\(X_1,X_2,X_3\)	—	—	—	\(f(X_1,X_2,X_3,X_4)-f(X_1,X_2,X_3)\)
3	\(X_1,X_2,X_4\)	—	—	\(f(X_1,X_2,X_3,X_4)-f(X_1,X_2,X_4)\)	—
3	\(X_1,X_3,X_4\)	—	\(f(X_1,X_2,X_3,X_4)-f(X_1,X_3,X_4)\)	—	—
3	\(X_2,X_3,X_4\)	\(f(X_1,X_2,X_3,X_4)-f(X_2,X_3,X_4)\)	—	—	—

Feature	Description	Units
`Por`	Average well porosity	\(\%\)
`Perm`	Log-transformed permeability to linearize relationships with other variables	\(\log_{10}(\mathrm{mD})\)
`AI`	Acoustic impedance	\(\mathrm{kg/m^3 \cdot m/s \cdot 10^6}\)
`Brittle`	Brittleness ratio	\(\%\)
`TOC`	Total organic carbon	\(\%\)
`VR`	Vitrinite reflectance	\(\%\)
`Prod`	Initial production, 90-day average	\(\mathrm{MCFPD}\)

Feature Ranking

Contents

Feature Ranking#

Motivation for Feature Ranking#

Curse of Dimensionality#

What is a Good Predictor Feature?#

What is Feature Ranking?#

Load the Required Libraries#

Design Custom Color Map#

Declare Functions#

Set the Working Directory#

Loading Tabular Data#

Visualize the DataFrame#

Summary Statistics for Tabular Data#

Coverage#

Feature Imputation#

Summary Statistics#

Univariate Distributions#

Bivariate Distributions#

Covariance#

Correlation Coefficient#

Rank Correlation Coefficient#

Partial Correlation Coefficient#

Semipartial Correlation Coefficient#

Feature Ranking with Feature Transformations#

Conditional Statistics and Conditional Distributions#

Variance Inflation Factor#

Model-based Feature Ranking Methods#

B Coefficients#

Beta Weights#

Feature Importance#

Mutual Information#

Mutual Information Accounting For Relevance and Redundancy#

Mutual Information Accounting For Relevance and Redundancy OFAT Variants#

Delta Mutual Information Quotient Accounting for Relevance and Redundancy#

Summary of All Ranking Methods#

Recursive Feature Elimination#

Shapley Values#

Shapley Values in Cooperative Game Theory#

Simple Shapley Value Examples#

Shapley Additive Partnership#

Shapley Synergistic Partnership#

Shapley Dysnergistic Partnership#

Shapley Values for Feature Ranking#

Removing Features from the Model#

Naive Shapley Approach for Removing Features#

Shapley with Feature Imputation for Removed Features#

Shapley with Tree-based Prediction Models#

Global Shapley by Summarizing Local Shapley#

Shapley Value for Feature Ranking Demonstration#

Build a Random Forest Model#

Calculate Shapley Values#

Local Shapley Values#

Custom Local Force Plot#

Shapley Force Plot#

Local Force Plot#

Global Shapley Values#

Comments#

About the Author#

Want to Work Together?#