Automated Dependence PlotsDavid I. Inouye
Purdue University
dinouye@purdue.edu

Liu Leqi, Joon Sik Kim
Carnegie Mellon University{leqil,joonsikk}@cs.cmu.edu

Abstract
In practical applications of machine learning,
it is necessary to look beyond standard metrics
such as test accuracy in order to validate various
qualitative properties of a model. Partial dependence plots (PDP), including instance-specific
PDPs (i.e., ICE plots), have been widely used as
a visual tool to understand or validate a model.
Yet, current PDPs suffer from two main drawbacks: (1) a user must manually sort or select interesting plots, and (2) PDPs are usually limited
to plots along a single feature. To address thesedrawbacks, we formalize a method for automating the selection of interesting PDPs and extend
PDPs beyond showing single features to show
the model response along arbitrary directions,for example in raw feature space or a latent
space arising from some generative model. We
demonstrate the usefulness of our automateddependence plots (ADP) across multiple usecases and datasets including model selection,
bias detection, understanding out-of-sample behavior, and exploring the latent space of a generative model. The code is available at https://github.com/davidinouye/adp.

1

INTRODUCTION

Modern applications of machine learning (ML) involve acomplex web of social, political, and regulatory issues aswell as standard technical issues such as covariate shift ortraining dataset bias. Although the most common validation metric is test set accuracy, the aforementioned issues
cannot be resolved by test set accuracy alone—and sometimes these issues directly oppose high test set accuracy(e.g., privacy concerns). Especially in certain applications
such as autonomous cars or automated loan approval,
practitioners have become concerned with unexpected
Proceedings of the 36th Conference on Uncertainty in Artificial
Intelligence (UAI), PMLR volume 124, 2020.

Bryon Aragam
University of Chicago
bryon@chicagobooth.edu

Pradeep Ravikumar
Carnegie Mellon University
pradeepr@cs.cmu.edu

and incorrect model behaviors, such as behaviors for data
points that were not seen during training or testing (e.g.,classifying a person as part of the road or approving alarge fraudulent loan). Thus, average-based validation
may be insufficient to validate a model. Rather, given
the difficulty of specifying expected model behavior apriori, qualitative methods for validating a model andhighlighting the most interesting or unusual model behaviors beyond average-based methods are needed.
One popular approach for qualitatively validating or understanding the effect of a particular feature on the model
response is a partial dependence plot (PDP) [8]. A
PDP plots the average model response across one feature marginalized over the other features and thus gives a
global view of the feature effect on the model response.
Because PDPs take an average over some features, they
may not be as effective at showing unusual behaviors that
may be important for certain applications as noted above.
To address this gap, instance-specific PDPs such as Individual Conditional Expectation (ICE) plots [10] could
be used to show the PDP plot for a single target instancerather than averaged over a set of instances. While PDPs
were invented many years ago, they have continued to be
widely used for understanding model behavior because
they are simple to interpret (e.g., [33]).
In order to investigate unusual behaviors for safety-criticalapplications, practitioners must manually inspect O(nd)instance-specific PDPs, where n is the number of samples
they want to investigate and d is the number of features.
This manual inspection is prohibitive, even for moderaten and d. Additionally, because PDPs only consider asingle feature (i.e., axis-aligned directions), they can missimportant interactions between features which may be
critical in certain applications. One could use 2D PDP
heatmaps, however, this would increase the number ofplots from O(nd) to O(nd2 ). Thus, despite widespreaduse, current PDPs suffer from two main drawbacks: (1) auser must manually select interesting plots, and (2) PDPs
are usually limited to plots along a single feature.

1. We formalize the concept of interestingness or unexpectedness of dependence plots by defining twoclasses of plot utility functions that have multiple
instantiations including utilities to compare two models and validate (or invalidate) certain properties suchas linearity, monotonicity, and Lipschitz continuity.
2. We generalize PDPs beyond axis-aligned directions
to consider the model response along sparse linear directions. In particular, we optimize a specified utility
measure over sparse linear directions using a greedy
coordinate pairs algorithm. For tabular data where
the features are inherently interpretable, we optimize
for a sparse linear direction in the raw feature space.
For rich data such as images or text, we propose tofind sparse linear directions in a latent representation
space (e.g., via a VAE) but show the corresponding images along this direction in the original raw
feature space.
3. We demonstrate the usefulness of the resulting automated dependence plots (ADPs) across multiple
use-cases and datasets including model selection,
bias detection, understanding out-of-sample behavior, and exploring the latent space of a generativemodel.
Related work. The importance of safety and robustness
in ML is now well-acknowledged [24, 31], especiallygiven the current dynamic regulatory environment surrounding the use of AI in production [13, 32]. A popular
approach to auditing and checking models before deployment is “explaining” a black-box model post-hoc. Bothearly and recent work in explainable ML rely on local
approximations (e.g., [11, 23]). Other recent papers havestudied the influence of training samples or features [5, 6].
These have been extended to subsets of informative features [2, 5, 19] to explain a prediction. Other approaches
employ influence functions [18] and prototype selection[34, 35]. A popular class of approaches take the perspective of local approximations, such as linear approximations [19, 23], and saliency maps [27, 28, 30]. A crucial

LeastMonotonicUtility=0.017
2.5
Score

To address these drawbacks, we formalize a method toautomate the selection of interesting PDP plots and toextend PDPs to directional dependence plots, which show
the model response in more general directions—either
in raw feature space for tabular data or in a latent space
for, e.g., visual or textual data. An illustrative example ofsuch a directional dependence plot is given in Fig. 1. Here,multiple features are being changed and the specific plotis optimized to show the least monotonic direction. This
optimization, which is one of our key contributions, is
how we automate the selection of “interesting” or relevantplots. We summarize our contributions as follows:

2.0

TargetModel
BestPartial

1.5
1659
19

1544
30

1430
1316
41
53
(amount, age)

1202
64

1087
75

Figure 1: The above example highlights the two key innovations of our proposed automated dependence plots: (1)Finding interesting plots by formalizing and optimizing
plot utility measures—in this case, the utility is basedon the non-monotonicity of the curve. (2) Optimizing
over directions that change multiple features rather than
only varying a single feature as in standard PDP plots (asindicated on the x-axis). For the target loan application(designated by the red point), the directional dependence
plot (solid line) shows the change in model scores (y-axis)
along a direction in which two numeric features (amountand age) are varied. The plot was optimized so that themodel response along the plot is the least monotonic (forcomparison, the best monotone regression line is shownvia the dotted line; see Sec 2 for more details).
caveat with such approximation-based methods is that thequality of the approximation is often unknown, and istypically only valid in small neighborhoods of an instance(although we note recent work on global approximations[12]). In contrast to these previous feature selection methods, our approach leverages a utility measure to selectfeatures or directions.

2

PLOT UTILITY MEASURESIn this section, we define two classes of utility measures
for generic dependence plots including PDP plots, ICE
plots, and our generalization to directional dependenceplots described in the next section. The intuition is thatthe plot utility measure quantifies the “interestingness” orrelevance of a particular feature or direction. As we will
discuss in more detail in Sec. 4, it is implicit here that wewish to maximize the utility (cf. (7)).
Notation. We assume that the input space is X = Rd —we could consider categorical variables by showing bar
charts as in categorical PDP plots but for simplicity,
we will focus on continuous features. For instancespecific plots, we will denote the target instance as
x0 ∈ Rd . Let f : Rd → R be a black-box function,i.e., we can query f to obtain pairs (x, f (x)). We do
not assume that f is differentiable or even continuous
in order to allow non-differentiable models such as ran-dom forests. We will denote a plot by it’s corresponding univariate function f˜(t) : R → R, where we visualize this function for some bounded interval, i.e., fort ∈ [a, b]—the bounds for standard PDP plots are usually based on the minimum and maximum values alongeach feature but we will generalize this for directionalplots in the next section. For other multivariate functions g(x), we will denote similarly the corresponding
univariate function as g̃(t). As an example, for PDP
f˜i (t) = E[f (. . . , xi−1 , t, xi+1 , · · · )] where the expectation is an empirical expectation with respect to somedataset (often a training dataset) and t ranges over theminimum and maximum of the i-th feature. For instancespecific PDP (ICE), f˜i (t) = f (. . . , xi−1 , t, xi+1 , · · · )where the input x is fixed to a target point x0 except
for the i-th feature—this can be seen as a PDP plot where
the dataset is a single point. We will denote U (f˜, a, b) as
a plot utility function where we usually suppress the dependence on the bounds a and b for simplicity and merely
denote the utility as U (f˜).
In the next subsections, we carefully develop and define
two general classes of plot utility measures. Both classes
of utilities compare to another plot, which we will show
as a dotted line as can be seen in Fig. 1. This helps
give the reader an interpretable reference for the utility
measure itself—e.g., showing the best bounded Lipschitz
approximation to the plot. Of course, we do not claim thatour proposed set of utility measures is exhaustive, andindeed, this paper will hopefully lead to more creativeand broader classes of utility measures.
2.1

MODEL CONTRAST UTILITY MEASURES

A natural way to measure the utility of a plot with respect to one model is to contrast it to the same plot basedon a different model. Given another multivariate model
gx0 (x) : Rd → R, which could be defined with respect toa target instance, we define the contrast utility measure:Rb
Uc (f˜) = a L( f˜(t), g̃x0 (t) )dt ,(1)
where L is a loss function, e.g., squared or absolute loss,
and g̃x0 (t) is the corresponding plot function with respect
to gx0 (x). By maximizing Uc (f˜), we can find the plotf˜ that differs the most from the baseline model g̃x0 (t).
The baseline could be another simple model like logistic regression or some local approximation of the modelaround the target instance such as explanations produced
by LIME [23]—this would provide a way to critique local
approximation methods and show where and how theydiffer from the true model.
Example 1 (Variance of Plot). While the PDP paper [8]did not propose ways for sorting or selecting interestingplots, a commonly used method to sort PDP curves isby the variance of the plot. The utility of a plot basedon variance can be seen as special case of the model
contrast utility where g̃(x) = c = E[f˜(s)] = 1/(b −Rba) a f˜(s)ds—i.e., an expectation of f˜ with respect toa uniform distribution between a and b—and the loss
function is squared error. We will use this as the default
sorting mechanism for comparison in our experiments
because it is the simplest and seems to be common inpractice.
Example 2 (Contrast with Constant Model). Sometimes
we may want to create a contrast model that is dependent
on the target point x0 . The simplest case is where the contrast model is a constant fixed to the original prediction
value, i.e., gx0 (x) = cx0 = f (x0 ) Note that in this casethe comparison function depends on x0 . This contrast toa constant model can find directions that deviate the mostfrom the prediction; this implicitly finds plots that are notflat and significantly affect the prediction value.
Example 3 (Contrast with Validated Linear Model). Suppose an organization has already deployed a carefully
validated linear model—i.e., the linear parameters were
carefully checked by domain experts to make sure themodel behaves reasonably with respect to all features.
The organization would like to improve the model’s performance by using a new model, but wants to see how the
new model compares to their carefully validated linear
model to see where it differs the most. In this case, the organization could let the contrast model be their previous
model, i.e., gx0 (x) = gLinear (x) where g does not dependon the target point x0 .
Example 4 (Contrast Random Forest and DNN). An
organization may want to compare two different model
classes such as random forests and deep neural networks(DNN) to diagnose if there are significant differences inthese model classes or if they are relatively similar. In
this case, the contrast model g(x) would be the random
forest or DNN.
Example 5 (Contrast with local approximations used for
explanations). We can also compare the true model with
explanation methods based on local approximation suchas LIME [23] or gradient-based explanation methods[26, 28, 30]. We can simply use the local approximation
to the model centered at the target point as the contrast
model,i.e., gx0 (x) = fˆx0 (x), where fˆx0 is the local approximation centered around x0 . Thus, the found diagnostic curve will show the maximum deviation of the truemodel from the local approximation model being used foran explanation. Importantly, this allows our diagnostic
method to assess the goodness of local approximation
explanation methods showing when they are reasonable
and when they may fail; see Fig. 2.

ModelComparisonUtility=0.018
Probability

0.7

TargetModel
Local Approx.

0.6
0.5
0.4
4

17.6

31.2
44.8
Loan Duration

58.4

72

Figure 2: This plot illustrates using the model contrast
utility where gx0 (dotted line) is an explanation model
based on the gradient similar to the local linear explanation models in LIME [23]. Notice how it shows where theapproximation may be appropriate (duration > 46) andwhere it might be far from the true model (duration < 46).
2.2

FUNCTIONAL PROPERTY(IN)VALIDATION UTILITY MEASURESIn many contexts, a user may be more interested in validating (or invalidating) certain functional properties of amodel, such as monotonicity or smoothness. For example,
if it is expected that a model should be increasing withrespect to a feature (e.g., income in a loan scoring model),then we’d like to check that this property holds true (atleast approximately). Let H be a class of univariate functions h̃ : R → R that represents a property that encodes
acceptable or expected behaviors. To measure deviationof a plot from this class of functions, take the minimum
expected loss over all h ∈ H:
Z b
(2)
Up (f˜) = minL(f˜(t), h̃(t))dth̃∈Hawhere as usual, L is a loss function. The minimization in(2) is a univariate regression problem that can be solved
using standard techniques. This utility will find dependence plots that maximally or minimally violate the functional properties encoded by H.
Example 6 (Linearity (in)validation via linear regression).
A user might want to view the plots that are the least linear
to see if there is some unusual plots that may need further
investigation. For this example, the class of functions
would merely be linear functions, i.e., H is the set of
univariate linear functions. This problem could be solved
easily using standard linear regression methods.
Example 7 (Monotonicity (in)validation via isotonic regression). In many applications, it may be known that themodel output should behave simply with respect to certain
features. For example, one might expect that the score ismonotonic in terms of income in a loan scoring model.
In this case, the class of functions should be the set ofmonotonic functions, i.e., H is the set of all monotonic
functions. The resulting problem can be efficiently solved
using isotonic regression [1]—and this is what we do in

our experiments; see Fig. 4 for an example of validating(or invalidating) the monotonic property.
Example 8 (Lipschitz-boundedness (in)validation via
constrained least squares). Another property that an organization might want to validate is whether the function
has a small Lipschitz constant along the curve. Formally,
they may wish to check if the following condition holds:f˜(t2 ) − f˜(t1 )≤ L,
t2 − t1∀t1 , t2 ∈ [a, b](3)

where L is a fixed Lipschitz constant. Thus, the corresponding class of functions HL is the set of Lipschitz
continuous functions with a Lipschitz constant less than
L. In practice, we can solve this problem via constrained
least squares optimization—similar to isotonic regression(details in supplement). An example of using Lipschitz
bounded functions for H can be seen in Fig. 3. This utility
will find the curve that maximally violates the bounded
Lipschitz condition; this curve may also be useful in finding where the model is particularly sensitive to changes
in the input since the derivative along the curve will belarger than the Lipschitz constant.
1.0
True Model
Best Approx. with L ≤ 5

0.5
0.0
−0.5
0.0

0.2

0.4

0.6

0.8

1.0

Figure 3: Example of contrast utility function when thefunction class is Lipschitz continuous with L ≤ 5.

3

DIRECTIONAL DEPENDENCE
PLOTS

In this section, we formally define a directional dependence plot, which generalizes classical PDPs and instancespecific PDPs. First, recall that a PDP f˜i varies thei-th feature and averages over the other features in adataset, i.e., f˜i (t) = Ex [f (x + tei )], where ei is thei-th standard basis vector. In this form of f˜i , the generalization to directions is quite straightforward as we
can replace ei by v where v is an arbitrary unit vector:f˜v (t) = Ex [f (x + tv)]. The instance-specific dependence plot with respect to x0 is:f˜v,x0 (t) = f (x0 + tv) .

(4)

In future sections, we will often suppress the dependenceon v or x0 if it is understood from the context. Because
this paper focuses on finding unusual or interesting plots

for the purposes of model validation especially in safetycritical applications, we will focus on the local instancespecific directional plots in our experiments but theseideas can clearly be applied to “global” directional plots.
We note that both for interpretability and computational
reasons, we will usually assume that v is sparse or lowdimensional in the following sections. In the next sections,
we will develop generalizations for directional plots inlatent spaces, parameter spaces of transformations, andother general considerations for directional plots.
Directions in a latent space. For rich data such as images, directions in the original input space (e.g., raw
pixels) may not be very informative. In these cases, it
may be more intuitive to move along directions in a latent space, such as one arising from a generative model¯
G : Rd 7→ Rd , that generates input vectors in Rd given¯
latent vectors in Rd , where d¯ is usually smaller than d(e.g., a low dimensional representation). Examples ofsuch models include VAEs [17, 22] and deep generativemodels via normalizing flows [7, 9, 14, 20, 21]. We then
optimize for directions in the latent space of the generative model G(·) rather than the raw input space itself.
This allows us to define directions in the latent space thatcan correspond to arbitrary curves in the latent space:f˜G,x0 (t) = f (G(z0 + tv)) ,(5)

where z0 = G−1 (x0 ) and G−1 (x0 ) is an (approximate)
input to the generative model that would have generated
x0 , i.e., G(G−1 (x0 )) ≈ x0 . For VAEs [17, 22], the decoder network acts as G and the encoder network acts as
an approximate G−1 . For normalizing flow-based models, both G(z) and G−1 (x0 ) can be computed exactly by
construction [7, 9, 14, 20, 21].
Directions in the parameter space of known transformations. In certain domains, there are natural classesof transformations that are semantically meaningful. Examples include adding white noise to an audio file or
blurring an image, removing semantic content like one
letter from a stop sign, or changing the style of an image
using style transfer [15]. We can consider directions in
the parameter space of a set of these known transformations, which could be arbitrarily complex and non-linear.
We will denote each transformation as λv : Rd → Rd
with parameter v ∈ [0, 1], where λ0 corresponds to theidentity transformation. For example, if λv a rotation operator, then v would represent the angle of rotation wherev = 0 is 0 degrees and v = 1 is 180 degrees. Givenan ordered set of ` different transformations denoted by
(1)
(2)
(`)
(λv , λv , . . . , λv ), we can define a plot function using

a composition of these simpler transformations:f˜v,x0 (t) = f (Λx0 (vt))(2) (1)Λx0 (v) , λ(`)
v` (. . . λv2 (λv1 (x0 )) . . . )(6)where v ∈ [0, 1]` is the optimization parameter. Thus,
directions in the parameter space correspond to non-linear
transformation curves in the original input space. We notethat the transformations can be arbitrarily non-linear andcomplex—even deep neural networks.
Realistic plot bounds. In some contexts, it is desirable
to explicitly audit the behavior of f off its training data
manifold. This is helpful for detecting unwanted bias
and checking robustness to rare, confusing, or adversarial examples—especially in safety critical applications
such as autonomous vehicles. In other contexts, we may
wish to ensure that the selected plots show realistic combinations of variables—that is, that they stay within thebounds of the training data distribution. For example, it
may not be useful to show an example of a high income
person who also receives social security benefits. Fortunately, it is not hard to enforce this constraint: First,
bound the endpoints of the plot using a box constraint
based on the minimum and maximum of each feature—
this is all that classical PDP plots consider. However, thissimple bounding may not be reasonable when we canmove along multiple features. Thus, we can also train adensity model on the training data to account for featuredependencies and further restrict the bounds of the plot toonly lie in regions with a density value above a threshold.
In our experiments, we use a simple multivariate Gaussian density estimate to create the boundaries of our plots
but more complex deep density estimators could be used
(e.g., [7, 14, 20, 21]).
Visualizing directional dependence plots. In typical
PDP plots, the feature values are shown on the horizontal
axis—yielding direct interpretability if the feature itself
is interpretable. For directional plots, we show the set offeature values that are changing that correspond to eachvalue of t (see the x-axis of Fig. 1). For example, if both
age and duration of a loan application are changing acrossa directional plot, we show pairs of age and duration along
the horizontal axis. For image-based examples, we show
the image in the original feature space corresponding todifferent t values to maintain interpretability.

4

OPTIMIZATION OF PLOT UTILITY

So far we have discussed: 1) how to measure the plot
utility and 2) how to generalize beyond axis-aligned directions to arbitrary directions. With these pieces in place,the objective is to find directional dependence plots thatare simultaneously interpretable and interesting. To dothis, we restrict the directions under consideration (usually to sparse directions) and optimize for the directionthat shows the best (either highest or lowest) utility:v ∗ = arg max U (f˜v,x0 ) ,(7)v

where v is restricted to some interpretable subset (e.g.,
sparse vectors). This optimization problem can also be
optimized over a set of target instances as follows:v ∗ , x∗0 = arg max U (f˜v,x0 ) ,

(8)

v,x0

where x is one instance in a specified dataset (e.g., a random sample of the training dataset). This finds both atarget data point and a direction where f exhibits unusual
or interesting behavior defined by the utility U , thus enabling users to see the worst case behavior of the model.
Note that the dataset used for optimization may not be
part of the training data because it may be new data where
no class labels are given.
Since f is assumed to be an arbitrary black-box, where
only function evaluations are possible, the optimization
problem even for a single target instance in (7) is usually nonconvex and nonsmooth. This model-agnostic
setup generalizes to many interesting settings includingones in which the model is private, distributed or nondifferentiable (e.g., boosted trees) such that exact gradient
information would be impossible or difficult to compute.
Given these restrictions, we must resort to zeroth-order
optimization. While we could use general zeroth-order
optimization techniques, we require that the directions
are sparse so that the resulting directional plots are interpretable. Thus for both computational and practicalreasons, we propose a greedy optimization scheme called
greedy coordinate pairs (GCP) that adds non-zero coordinates in a greedy fashion, as outlined in Algorithm 1. We
initialize this algorithm by finding the single feature withthe highest utility—this is the same as computing the bestaxis-aligned direction for standard PDP plots.
Computational complexity. The total complexity of
the GCP algorithm—which is easily parallelized—isO(dkM I), where M is the number of angles tested and I
is the number of iterations. In our experiments, GCP typically found an interesting direction within 10 iterations.
By comparison, the computational complexity of standard
PDP plots for all dimensions in terms of the number ofmodel evaluations is O(ndk). Thus, compared to PDP
plots, our method is at most O(M I) slower.
Optimization evaluation. To test the effectiveness of
GCP (Algorithm 1), we ran two tests. First, we compared
the optimal utility values returned by GCP to the utility

Algorithm 1 Greedy Coordinate Pairs (GCP)
Input: Plot utility U , target point x0 , max number offeatures D, grid size M
Output: Optimized direction v ∗G(i, j, θ) ≡ Givens rotation matrix for coordinates iand j with angle θπ 2π
Θ = {0, M, M , · · · , π}v ← arg max i∈{1,...,d} U (f˜ei ) {Coordinate-wise optimization}
while v not converged doi∗ , j ∗ , θ∗ ←arg max i,j,θ∈Θ U (f˜G(i,j,θ)v )Ps.t.
k I([G(i, j, θ)v]k 6= 0) ≤ D
v ← G(i∗ , j ∗ , θ∗ )v
end whileof 10,000 randomly selected directional plots based on arandom forest classifier. In all cases, GCP returned thehighest values. For example, GCP found a directional
plot f˜ with U (f˜) = 0.0016, compared with 0.0007 forthe best random curve. In the second experiment, we
generated random synthetic models in which the directional plot with the highst utility could be determined inadvance, making it possible to evaluate how often GCP
selects the optimal one. We evaluated its performanceon directions with at most one, two, or three nonzero
coordinates, and found that in 100%, 97%, and 98% of
simulations, GCP found the optimal directions. In the few
cases GCP did not find the optimal directions, this was dueto the randomness in generating examples whose optimal
directions are more difficult to identify (e.g.„ the optimal
curve was nearly constant). Details and further results onthese experiments can be found in the appendix. Thus,
we conclude that despite our algorithm being greedy, our
GCP algorithm is empirically reasonable.

5

EXPERIMENTS

We present five concrete use cases: 1) Selecting among
several models, 2) Bias detection, 3) Out-of-sample behavior in computer vision, 4) Discovering interesting model
properties via generative models, and 5) Evaluating robustness to covariate shift. We have put (5) and somedetails of the experiments in the appendix. A python module implementing the proposed framework is available athttps://github.com/davidinouye/adp.
5.1

SELECTING A MODEL FOR LOAN
PREDICTION

Suppose we have trained multiple models, each with similar test accuracies. Which of these models should be deployed? To answer this question, directional dependenceADP

RandomForest (train=0.81, test=0.63)
LeastMonotonicUtility=0.18
Target
Model
Mono.

2
0
−2

Target
Model
Mono.

0

4330 8410 12490 16570 20649
amountRandomForest (train=0.81, test=0.63)
Variance=0.121

PDP

2

PDP

0.6
0.4

3366 6481 9597 12712 15828
amount5000

10000
amount

PDP

0.6

15000

RandomForest (train=0.81, test=0.63)
Variance=1.201
ICE
Target

2

0

5000

10000
amount15000

2503514 6778 10042 13306 16570
amountDNN (train=0.83, test=0.62)
Variance=0.160
0.8

PDP

0.4
0.2
0

5000

10000
amount15000

20

GradientBoost (train=0.77, test=0.62)
Variance=0.988
2ICE
Target40
duration60

DNN (train=0.83, test=0.62)
Variance=3.698
5ICE
Target

0
−5
−10

−2

−2

Target
Model
Mono.

0.6

0

0

0

GradientBoost (train=0.77, test=0.62)
Variance=0.110

0.2
0

5

−10
250

0.4

0.2

DNN (train=0.83, test=0.62)
LeastMonotonicUtility=2.4

−5

−2
250

ICE

GradientBoost (train=0.77, test=0.62)
LeastMonotonicUtility=0.068

0

5000

10000
amount

15000

0

5000

10000
amount15000

Figure 4: Our dependence plots (top, gray) selected by finding the least monotonic direction elucidates potentiallyproblematic behavior of the various models whereas PDP (middle) and instance-specific ICE (bottom) plots do nothighlight this non-monotonic behavior. The PDP and ICE plots were selected according to the simple highest varianceutility, which is used in practice for standard PDP plots. All plots are based on 100 random samples, where the PDP
plots average over the samples and our plots and the ICE plots additionally select a specific target point. The dotted line
(labeled “Mono.”) is the best monotonic regression plot.
plots can be used to detect undesirable or unexpected behaviours. We explore this use-case of qualitative model
selection on a dataset of German loan application data inorder to find non-monotonic plots. For example, it may
be undesirable to deploy a model that decreases a candidate’s score in response to an increase in their income.
We train a random forest, a gradient boosted tree, and a
deep neural network. For this simple example, we consider axis-aligned directions with only one non-zero sothat we can compare to PDP and ICE plots; additionally,
we optimize the utility over 100 random target points x0
as in Eq. 8—thus providing an estimate of the worst-casebehavior of the model. The test accuracies of these models were all close to 62%—thus, a stakeholder would not
be able to select a model based on accuracy alone. The
plots for our method using the monotonicity utility compared to PDP and ICE using the simple variance utility
can be seen in Fig. 4. In addition to a single number thatquantifies non-monotonicity, our directional dependenceplots show the user both the location and severity of theworst-case non-monotonicity. On the other hand, PDP

and instance-specific PDP (ICE) plots selected by variance prefer models that have large ranges, but may haveexpected or uninteresting patterns within this range. By
contrast, our dependence plots can highlight more subtle
and interesting patterns such as non-monotonicity. Our
directional dependence plots suggest that random forests
or gradient boosted trees may be preferable since theirworst-case behaviors are nearly monotonic, whereas the
DNN model is far from monotonic.
5.2

BIAS DETECTION IN RECIDIVISM
MODEL PREDICTION

In many real-world applications, it is essential to checkmodels for bias. A contemporary example of significant
interest in recent years concerns recidivism prediction
instruments (RPIs), which use data to predict the likelihood of a convicted criminal to re-offend [3, 4]. Given
an instance x0 = (u, x2 , . . . , xd ), consider what the output of f would have been had the protected attributeu ∈ {0, 1} (e.g., race or gender) been flipped. In certain

SVM Score

ModelComparisonUtility=0.32
0.0
−0.5

Target (White)
Original Race
Other Race

−1.0
1

58.2

115.4
172.6
229.8
Time Served (Month)287

Figure 5: A directional plot using the model contrast
utility where the comparison model is the same model but
the race is flipped. Notice that bias between races is far
from uniform even switching bias direction. Expanded
figure with flipping both gender and race on two target
points can be found in the appendix.
cases, such protected attributes might not be explicitly
available, in which case, we could use proxy attributes orfactors, though we do not focus on that here. A modelthat ignores the protected attributes would see little orno change in f as a result of this change. In this simple
example, we explicitly ignore dependencies between theprotected attribute and other features though this would be
important to consider for any significant conclusions to be
made. Given this situation, we select the model contrast
utility (Sec 2) with a special definition for the comparison
model defined as follows: gu0 (u) = fu0 (σ(u)), where[σ(u)]i is 1−ui if the i-th feature is the protected attribute,and ui otherwise, essentially flipping only the protected
attribute and leaving all other features untouched. There
are two cases: Either (a) No such plot deviates too much,in which case this is evidence that f is not biased, or(b) there is a dimension along which f is significantly
biased. A directional plot for flipping race from white toblack based on data from [25] using a kernel SVM model
can be seen in Fig 5. One can clearly see that the model
behaves in biased ways with respect to race: The effectof time served on the risk score clearly depends on raceand even switches bias direction. For this high stakes
application, the use of such a model could be problematic,
and our directional plots highlight this fact easily evenfor non-technical users. Finally, these directional plots
avoid averaging the data into a single metric and offer
more insight into the location and form of the model biasthat may depend on the inputs in complex ways.
5.3

UNDERSTANDING OUT-OF-SAMPLE
BEHAVIOR FOR TRAFFIC SIGN
DETECTION

When the model is deployed in practice, understanding
how the model behaves outside of the training dataset isoften critical. For example, for a traffic sign detector on
autonomous vehicles, will the detector predict correctly if
a traffic sign is rotated even though the training data did
not contain any rotated signs? What other variations of thedata is the detector susceptible to? For this use-case, wefirst trained a convolution neural network on German Traffic Sign Recognition Benchmark (GTSRB) dataset [29]which achieves 97% test accuracy, that will act as the stop
sign detector. We consider transformation curves based on
five standard image transformations: rotate, blur, brighten,
desaturate, and contrast. Each of these transformations
creates images outside of the training data, and we areinterested in finding the interesting combinations of thesetransformations that influence the prediction score of thedetector in different ways. Fig. 6 depicts the resulting
plots we generated through optimizing the corresponding
utilities. The least constant direction (top) simultaneously
adjusts contrast and brightness, which is expected since
this transformation gradually removes almost all information about the stop sign. The most constant direction(middle) adjusts saturation and contrast, which may besurprising since it suggests that the model ignores color.
Finally, the least monotonic direction (bottom) is rotation,
which suggests that the model identifies a white horizontal region but ignores the actual letters “STOP” since it
still predicts correctly when the stop sign is flipped over.
5.4

DISCOVERING INTERESTING MODELPROPERTIES VIA GENERATIVE MODELSOur directional dependence plots can also search nonlinear directions (i.e., curves) in the input space, by discovering linear directions in some latent feature spacelearned by generative models. This can be particularly
useful when linear directions in the input feature space
are less semantic compared to non-linear directions, like
images. In this example, we demonstrate how directionalplots effectively select the most interesting non-linear
directions in the pixel space.
On the MNIST data, we train a VAE [17] for learninga 10-dimensional latent data representation, as well as a
CNN classifier. The directions in the latent space of deep
networks are found to be semantic in several ways [16],so we optimize in this space for the least constant, mostconstant, and least monotonic in terms of the prediction
score change of the classifier. In Figure 7, given a test
instance of digit 1, we show the varying prediction probability (left) along with the set of reconstructed images
along the direction (right). Least-constant utility (top row)
discovers directions that add distinctive curves and loops
to the image inducing a steep drop in the prediction probability, while the most-constant utility (middle row) finds
directions that keep the relative shape of the digit 1 con-

Figure 6: The dependence plots for traffic sign detection show which transformation is the most sensitive (top), theleast sensitive (middle), and the least monotonic (bottom). These curves highlight both expected model trends (top)and unexpected trends (middle and bottom), where the model seems to ignore color and fails when the stop is rotatedpartially but works again when the stop sign is almost flipped over.

Figure 7: Discovering various non-linear directions in images using VAE. Given the target instance of digit 1, theleast-constant utility highlights that adding more curves and loops to the image induces the prediction probability todrop the most (top row), while the most-constant utility identifies directions that preserve the relative shapes of the digit(middle row). The least-monotonic utility exposes the directions that make the prediction probability fluctuate the most– moving along these directions alter the digit to different numbers, e.g., 5 and 3.
stant throughout. Least-monotonic utility (bottom row)
successfully exposes directions that change the digit 1 to 5
and 3, inducing the most non-monotonic prediction probability changes. Note that such discovery is not limited
to images — as long as some semantically meaningful
feature embedding is available (e.g., word embedding),
the same process can be used to automatically discover interesting directions within the input feature space. Whilethese use cases demonstrate our method across a wide
range of scenarios, a user study for various contexts suchas exploratory data analysis, model checking, or modelcomparison would be an excellent area for future research.

Overall, we believe our framework opens up new possibilities and challenges for automating and extending PDP
plots.

Acknowledgements
The authors acknowledge the support of DARPA via
FA87501720152 and Accenture. DII acknowledges support from Northrop Grumman. JSK acknowledges support from Kwanjeong Educational Foundation. BA acknowledges support from the Robert H. Topel Faculty
Research Fund.

References[1] M. J. Best and N. Chakravarti. Active set algorithms forisotonic regression; a unifying framework. Mathematical
Programming, 47(1-3):425–439, 1990.
[2] J. Chen, L. Song, M. Wainwright, and M. Jordan. Learning
to explain: An information-theoretic perspective on modelinterpretation. In ICML, pages 883–892, 2018.
[3] A. Chouldechova. Fair prediction with disparate impact:
A study of bias in recidivism prediction instruments. Big
data, 5(2):153–163, 2017.
[4] A. Chouldechova and M. G’Sell. Fairer and more accurate,
but for whom? Workshop on Fairness, Accountability, andTransparency in Machine Learning, 2017.
[5] A. Datta, S. Sen, and Y. Zick. Algorithmic transparency
via quantitative input influence: Theory and experimentswith learning systems. In IEEE Symposium on Security
and Privacy (SP), pages 598–617. IEEE, 2016.

[18] P. W. Koh and P. Liang. Understanding black-box predictions via influence functions. In Proceedings of the 34th
International Conference on Machine Learning-Volume
70, pages 1885–1894. JMLR. org, 2017.
[19] S. M. Lundberg and S.-I. Lee. A unified approach tointerpreting model predictions. In NIPS, 2017.
[20] J. B. Oliva, A. Dubey, M. Zaheer, B. Póczos, R. Salakhutdinov, E. P. Xing, and J. Schneider. Transformation autoregressive networks. In ICML, 2018.
[21] G. Papamakarios, T. Pavlakou, and I. Murray. Masked
autoregressive flow for density estimation. In NIPS, 2017.
[22] D. J. Rezende, S. Mohamed, and D. Wierstra. Stochastic backpropagation and approximate inference in deep
generative models. In ICML, 2014.
[23] M. T. Ribeiro, S. Singh, and C. Guestrin. Why should itrust you?: Explaining the predictions of any classifier. In
KDD, 2016.

[6] A. Dhurandhar, P.-Y. Chen, R. Luss, C.-C. Tu, P. Ting,
K. Shanmugam, and P. Das. Explanations based on the
missing: towards contrastive explanations with pertinent
negatives. In NeurIPS, pages 590–601. Curran Associates
Inc., 2018.

[24] S. Saria and A. Subbaswamy. Tutorial: Safe and reliablemachine learning. arXiv preprint arXiv:1904.07204, 2019.

[7] L. Dinh, J. Sohl-Dickstein, and S. Bengio. Density estimation using real NVP. In ICLR, 2017.

[26] A. Shrikumar, P. Greenside, and A. Kundaje. Learning important features through propagating activation differences.
ICML, 2017.

[8] J. H. Friedman. Greedy function approximation: A gradient boosting machine. The Annals of Statistics, 29(5):1189–1232, 2001.
[9] M. Germain, K. Gregor, I. Murray, and H. Larochelle.
MADE: Masked autoencoder for distribution estimation.
In ICML, 2015.
[10] A. Goldstein, A. Kapelner, J. Bleich, and E. Pitkin. Peeking inside the black box: Visualizing statistical learningwith plots of individual conditional expectation. Journalof Computational and Graphical Statistics, 24(1):44–65,
jan 2015.
[11] R. Guidotti, A. Monreale, S. Ruggieri, D. Pedreschi,
F. Turini, and F. Giannotti. Local rule-based explanationsof black box decision systems, 2019.

[25] P. Schmidt and A. D. Witte. Predicting Recidivism in North
Carolina, 1978 and 1980. Inter-university Consortium for
Political and Social Research, 1988.

[27] K. Simonyan, A. Vedaldi, and A. Zisserman. Deep insideconvolutional networks: Visualising image classification
models and saliency maps. Workshop at ICLR, 2014.
[28] D. Smilkov, N. Thorat, B. Kim, F. Viégas, and M. Wattenberg. Smoothgrad: removing noise by adding noise.
Visualization for Deep Learning Workshop at ICML, 2017.
[29] J. Stallkamp, M. Schlipsing, J. Salmen, and C. Igel. The
German Traffic Sign Recognition Benchmark: A multiclass classification competition. In IEEE International
Joint Conference on Neural Networks, pages 1453–1460,
2011.
[30] M. Sundararajan, A. Taly, and Q. Yan. Axiomatic attribution for deep networks. ICML, 2017.

[12] W. Guo, S. Huang, Y. Tao, X. Xing, and L. Lin. Explaining
deep learning models–a bayesian non-parametric approach.
In NeurIPS, pages 4514–4524, 2018.

[31] K. R. Varshney and H. Alemzadeh. On the safety of machine learning: Cyber-physical systems, decision sciences,and data products. Big data, 5(3):246–255, 2017.

[13] G. Hadfield and J. Clark. Regulatory markets for AI safety.
Safe Machine Learning Workshop at ICLR, 2019.

[32] S. Wachter, B. Mittelstadt, and C. Russell. Counterfactual
explanations without opening the black box: Automateddecisions and the gdpr. SSRN Electronic Journal, 11 2017.
doi: 10.2139/ssrn.3063289.

[14] D. I. Inouye and P. Ravikumar. Deep density destructors.
In ICML, pages 2172–2180, jul 2018.
[15] J. Johnson, A. Alahi, and L. Fei-Fei. Perceptual losses forreal-time style transfer and super-resolution. In B. Leibe,J. Matas, N. Sebe, and M. Welling, editors, Computer Vision – ECCV 2016, pages 694–711, Cham, 2016. Springer
International Publishing. ISBN 978-3-319-46475-6.
[16] B. Kim, M. Wattenberg, J. Gilmer, C. Cai, J. Wexler,
F. Viegas, et al. Interpretability beyond feature attribution: Quantitative testing with concept activation vectors(tcav). In International Conference on Machine Learning,
pages 2668–2677, 2018.
[17] D. P. Kingma and M. Welling. Auto-encoding variational
bayes. In International Conference on Learning Representations (ICLR), 2014.

[33] J. Wexler, M. Pushkarna, T. Bolukbasi, M. Wattenberg,
F. Viégas, and J. Wilson. The what-if tool: Interactive
probing of machine learning models. IEEE transactions on
visualization and computer graphics, 26(1):56–65, 2019.
[34] C.-K. Yeh, J. Kim, I. E.-H. Yen, and P. K. Ravikumar.
Representer point selection for explaining deep neuralnetworks. In NeurIPS, pages 9291–9301, 2018.
[35] X. Zhang, A. Solar-Lezama, and R. Singh. Interpreting
neural network judgments via minimal, stable, and symbolic corrections. In NeurIPS, pages 4874–4885, 2018.