from __future__ import division import numpy as np import warnings try: import matplotlib.pyplot as pl import matplotlib except ImportError: warnings.warn("matplotlib could not be loaded!") pass from . import labels from . import colors from ..common import convert_name, approximate_interactions def dependence_plot(ind, shap_values, features, feature_names=None, display_features=None, interaction_index="auto", color="#1E88E5", axis_color="#333333", cmap=None, dot_size=16, x_jitter=0, alpha=1, title=None, xmin=None, xmax=None, ax=None, show=True): """ Create a SHAP dependence plot, colored by an interaction feature. Plots the value of the feature on the x-axis and the SHAP value of the same feature on the y-axis. This shows how the model depends on the given feature, and is like a richer extenstion of the classical parital dependence plots. Vertical dispersion of the data points represents interaction effects. Grey ticks along the y-axis are data points where the feature's value was NaN. Parameters ---------- ind : int or string If this is an int it is the index of the feature to plot. If this is a string it is either the name of the feature to plot, or it can have the form "rank(int)" to specify the feature with that rank (ordered by mean absolute SHAP value over all the samples). shap_values : numpy.array Matrix of SHAP values (# samples x # features). features : numpy.array or pandas.DataFrame Matrix of feature values (# samples x # features). feature_names : list Names of the features (length # features). display_features : numpy.array or pandas.DataFrame Matrix of feature values for visual display (such as strings instead of coded values). interaction_index : "auto", None, int, or string The index of the feature used to color the plot. The name of a feature can also be passed as a string. If "auto" then shap.common.approximate_interactions is used to pick what seems to be the strongest interaction (note that to find to true stongest interaction you need to compute the SHAP interaction values). x_jitter : float (0 - 1) Adds random jitter to feature values. May increase plot readability when feature is discrete. alpha : float The transparency of the data points (between 0 and 1). This can be useful to the show density of the data points when using a large dataset. xmin : float or string Represents the lower bound of the plot's x-axis. It can be a string of the format "percentile(float)" to denote that percentile of the feature's value used on the x-axis. xmax : float or string Represents the upper bound of the plot's x-axis. It can be a string of the format "percentile(float)" to denote that percentile of the feature's value used on the x-axis. ax : matplotlib Axes object Optionally specify an existing matplotlib Axes object, into which the plot will be placed. In this case we do not create a Figure, otherwise we do. """ if cmap is None: cmap = colors.red_blue # create a matplotlib figure, if `ax` hasn't been specified. if not ax: figsize = (7.5, 5) if interaction_index != ind else (6, 5) fig = pl.figure(figsize=figsize) ax = fig.gca() else: fig = ax.get_figure() # convert from DataFrames if we got any if str(type(features)).endswith("'pandas.core.frame.DataFrame'>"): if feature_names is None: feature_names = features.columns features = features.values if str(type(display_features)).endswith("'pandas.core.frame.DataFrame'>"): if feature_names is None: feature_names = display_features.columns display_features = display_features.values elif display_features is None: display_features = features if feature_names is None: feature_names = [labels['FEATURE'] % str(i) for i in range(shap_values.shape[1])] # allow vectors to be passed if len(shap_values.shape) == 1: shap_values = np.reshape(shap_values, len(shap_values), 1) if len(features.shape) == 1: features = np.reshape(features, len(features), 1) ind = convert_name(ind, shap_values, feature_names) # plotting SHAP interaction values if len(shap_values.shape) == 3 and len(ind) == 2: ind1 = convert_name(ind[0], shap_values, feature_names) ind2 = convert_name(ind[1], shap_values, feature_names) if ind1 == ind2: proj_shap_values = shap_values[:, ind2, :] else: proj_shap_values = shap_values[:, ind2, :] * 2 # off-diag values are split in half # TODO: remove recursion; generally the functions should be shorter for more maintainable code dependence_plot( ind1, proj_shap_values, features, feature_names=feature_names, interaction_index=ind2, display_features=display_features, ax=ax, show=False, xmin=xmin, xmax=xmax ) if ind1 == ind2: ax.set_ylabel(labels['MAIN_EFFECT'] % feature_names[ind1]) else: ax.set_ylabel(labels['INTERACTION_EFFECT'] % (feature_names[ind1], feature_names[ind2])) if show: pl.show() return assert shap_values.shape[0] == features.shape[0], \ "'shap_values' and 'features' values must have the same number of rows!" assert shap_values.shape[1] == features.shape[1], \ "'shap_values' must have the same number of columns as 'features'!" # get both the raw and display feature values oinds = np.arange(shap_values.shape[0]) # we randomize the ordering so plotting overlaps are not related to data ordering np.random.shuffle(oinds) xv = features[oinds, ind].astype(np.float64) xd = display_features[oinds, ind] s = shap_values[oinds, ind] if type(xd[0]) == str: name_map = {} for i in range(len(xv)): name_map[xd[i]] = xv[i] xnames = list(name_map.keys()) # allow a single feature name to be passed alone if type(feature_names) == str: feature_names = [feature_names] name = feature_names[ind] # guess what other feature as the stongest interaction with the plotted feature if interaction_index == "auto": interaction_index = approximate_interactions(ind, shap_values, features)[0] interaction_index = convert_name(interaction_index, shap_values, feature_names) categorical_interaction = False # get both the raw and display color values color_norm = None if interaction_index is not None: cv = features[:, interaction_index] cd = display_features[:, interaction_index] clow = np.nanpercentile(cv.astype(np.float), 5) chigh = np.nanpercentile(cv.astype(np.float), 95) if type(cd[0]) == str: cname_map = {} for i in range(len(cv)): cname_map[cd[i]] = cv[i] cnames = list(cname_map.keys()) categorical_interaction = True elif clow % 1 == 0 and chigh % 1 == 0 and chigh - clow < 10: categorical_interaction = True # discritize colors for categorical features if categorical_interaction and clow != chigh: clow = np.nanmin(cv.astype(np.float)) chigh = np.nanmax(cv.astype(np.float)) bounds = np.linspace(clow, chigh, int(chigh - clow + 2)) color_norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N-1) # optionally add jitter to feature values if x_jitter > 0: if x_jitter > 1: x_jitter = 1 xvals = xv.copy() if isinstance(xvals[0], float): xvals = xvals.astype(np.float) xvals = xvals[~np.isnan(xvals)] xvals = np.unique(xvals) if len(xvals) >= 2: smallest_diff = np.min(np.diff(np.sort(xvals))) jitter_amount = x_jitter * smallest_diff xv += (np.random.ranf(size = len(xv))*jitter_amount) - (jitter_amount/2) # the actual scatter plot, TODO: adapt the dot_size to the number of data points? xv_nan = np.isnan(xv) xv_notnan = np.invert(xv_nan) if interaction_index is not None: # plot the nan values in the interaction feature as grey cvals = features[oinds, interaction_index].astype(np.float64) cvals_imp = cvals.copy() cvals_imp[np.isnan(cvals)] = (clow + chigh) / 2.0 cvals[cvals_imp > chigh] = chigh cvals[cvals_imp < clow] = clow p = ax.scatter( xv[xv_notnan], s[xv_notnan], s=dot_size, linewidth=0, c=cvals[xv_notnan], cmap=cmap, alpha=alpha, vmin=clow, vmax=chigh, norm=color_norm, rasterized=len(xv) > 500 ) p.set_array(cvals[xv_notnan]) else: p = ax.scatter(xv, s, s=dot_size, linewidth=0, color=color, alpha=alpha, rasterized=len(xv) > 500) if interaction_index != ind and interaction_index is not None: # draw the color bar if type(cd[0]) == str: tick_positions = [cname_map[n] for n in cnames] if len(tick_positions) == 2: tick_positions[0] -= 0.25 tick_positions[1] += 0.25 cb = pl.colorbar(p, ticks=tick_positions) cb.set_ticklabels(cnames) else: cb = pl.colorbar(p) cb.set_label(feature_names[interaction_index], size=13) cb.ax.tick_params(labelsize=11) if categorical_interaction: cb.ax.tick_params(length=0) cb.set_alpha(1) cb.outline.set_visible(False) bbox = cb.ax.get_window_extent().transformed(fig.dpi_scale_trans.inverted()) cb.ax.set_aspect((bbox.height - 0.7) * 20) # handles any setting of xmax and xmin # note that we handle None,float, or "percentile(float)" formats if xmin is not None or xmax is not None: if type(xmin) == str and xmin.startswith("percentile"): xmin = np.nanpercentile(xv, float(xmin[11:-1])) if type(xmax) == str and xmax.startswith("percentile"): xmax = np.nanpercentile(xv, float(xmax[11:-1])) if xmin is None or xmin == np.nanmin(xv): xmin = np.nanmin(xv) - (xmax - np.nanmin(xv))/20 if xmax is None or xmax == np.nanmax(xv): xmax = np.nanmax(xv) + (np.nanmax(xv) - xmin)/20 ax.set_xlim(xmin, xmax) # plot any nan feature values as tick marks along the y-axis xlim = ax.get_xlim() if interaction_index is not None: p = ax.scatter( xlim[0] * np.ones(xv_nan.sum()), s[xv_nan], marker=1, linewidth=2, c=cvals_imp[xv_nan], cmap=cmap, alpha=alpha, vmin=clow, vmax=chigh ) p.set_array(cvals[xv_nan]) else: ax.scatter( xlim[0] * np.ones(xv_nan.sum()), s[xv_nan], marker=1, linewidth=2, color=color, alpha=alpha ) ax.set_xlim(xlim) # make the plot more readable ax.set_xlabel(name, color=axis_color, fontsize=13) ax.set_ylabel(labels['VALUE_FOR'] % name, color=axis_color, fontsize=13) if title is not None: ax.set_title(title, color=axis_color, fontsize=13) ax.xaxis.set_ticks_position('bottom') ax.yaxis.set_ticks_position('left') ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax.tick_params(color=axis_color, labelcolor=axis_color, labelsize=11) for spine in ax.spines.values(): spine.set_edgecolor(axis_color) if type(xd[0]) == str: ax.set_xticks([name_map[n] for n in xnames], xnames, rotation='vertical', fontsize=11) if show: with warnings.catch_warnings(): # ignore expected matplotlib warnings warnings.simplefilter("ignore", RuntimeWarning) pl.show()