Why don't partial dependence plots match model predictions?

Background

My training is in statistics, but I'm interested in machine learning. My models involve nonlinear relationships between 2-5 predictors and a single response (all variables continuous). There are clear, direct causal links between each of the predictors and the response, which are well-known - but the interactions between the predictors are complex and not well-understood. I am interested in fitting machine learning models and then exploring the relationship between interacting predictors and the response (I am also exploring this theoretically to develop process-based models).

Problem

I can see two ways to to characterise how a predictor relates to a response in a fitted machine learning model

1. Partial dependence plots. Often these are calculated with just one predictor at a time, and if I understand correctly, calculating the relationship between the single predictor and the response implicitly involves 'setting' the remaining predictors at their median value. However, there is no reason why this is necessary - one can calculate the partial dependence across a grid of all predictor values. This is interesting to me because it lets me explore how the predictors interact to influence the response.

2. The second way is to simply use the fitted model to make predictions on a new dataset across a gradient in the predictor. For multiple predictors (my case), one can use the fitted model to make predictions across a grid of all predictor values.

My question is: why do the two approaches above produce different results? I would expect them to be identical and would like to understand which is a better approach for my purposes.

Here is an example with 2 predictors showing that the two approaches produce different results. I use a random forest for the example but I don't think this is specific to that type of model and would be interested in a more general answer.

# Load libraries and set graphical defaults
library(randomForest)
library(pdp)
library(ggplot2)
theme_set(theme_bw(base_size = 15))

# Set seed
set.seed(6253)

# Generate dataset with x1 and x2 as predictors and y as response
dat <- data.frame(x1 = rnorm(10000, 0, 1),
x2 = rnorm(10000, 0, 1))
dat$$y <- (dat$$x1 * 8) + (dat$$x2 * 5) + (dat$$x1 * dat$x2 * 3) + rnorm(10000, 0, 2) # Checking that the predictors are strongly related to the response summary(lm(y ~ x1*x2, data = dat)) plot(dat$y ~ predict(lm(y ~ x1*x2, data = dat)))

# Fitting random forest
rf1 <- randomForest(y ~ x1 + x2, data = dat)

# Generate new dataset with grid of x1 and x2 values, to fill with model predictions
newdat <- expand.grid(x1 = seq(-1, 1, 0.2),
x2 = seq(-1, 1, 0.2))

# Make predictions using partial dependence (plot) approach
# Uses pdp package but results do not seem to depend on the package
newdat$$pdp_preds <- partial(rf1, pred.var = c("x1", "x2"), pred.grid = newdat[,1:2])$$yhat

# Make predictions by just using the fitted model and predict()
newdat$rf_preds <- predict(rf1, newdat) # Plot the two types of predictions against each other # They are correlated but far from identical ggplot(newdat, aes(rf_preds, pdp_preds)) + geom_point(size = 3, alpha = 0.6) + geom_abline(slope = 1) + xlab('Using predict() approach') + ylab('Using partial dependence plot approach') # Visualise the x1*x2 surface from the partial dependence calculations ggplot(newdat, aes(x1, x2, z = pdp_preds)) + geom_contour_filled() + ggtitle('Using partial dependence plot approach') # Visualise the x1*x2 surface from predict() ggplot(newdat, aes(x1, x2, z = rf_preds)) + geom_contour_filled() + ggtitle('Using predict() approach')  Plots Comparing the two approaches: x1*x2 surface from the partial dependence plot approach: x1*x2 surface from the predict() approach: 1 Answer It is an ordering issue, partial does not return the data in the order you think it does. newdat is ordered first by x1 and then x2, while the results from partial are ordered first by x2 and then x1. To that extent, from the currently produced PDP surfaces, we can actually see that one is the transposed version of the other. Being more particular, if we define: QQ = partial(rf, pred.var = c("x1", "x2"), pred.grid = newdat) as its output and then plot the PDP surface as ggplot(newdat, aes(QQ$x1, QQ$x2, z = QQ$yhat)) + geom_contour_filled() + ggtitle('Using partial dependence plot approach') we get the same result as with the direct approach.

• Oof, thanks! This is definitely the issue. Quick question to confirm my understanding: do you agree that the two approaches (using model to predict and using partial dependence calculations) should in principle always produce identical output?
– mkt
Aug 13 at 10:52
• Almost... they should here but there are some caveats. PDPs do not always capture reasonable interactions between features, and especially tail-behaviour insights can be misleading. In this example x1 and x2 are essentially uncorrelated so we are fine but especially with more realistic sample doing some sort of marginalisation helps us get more generalisable and less variable insights. Aug 13 at 20:37
• Thanks, this is useful! Would you be able to point me to any useful resources to understand these caveats better? I'm happy to post this as a new Q too.
– mkt
Aug 13 at 22:05
• Glad I could help! :) I haven't find many good resources, IML by Molnar is very good. Read the chapter on PDPs (8.1 at the time of writing) and then the chapter on ALEPs (8.2 at the time of writing this), they cover these issues very nicely. The (few) references within should be enough to carry you forward. In general, PDPs have been almost an afterthought of a long time. The original GBM paper from Friedman (2001) gives them about 3 pages and then the matter was considered a done deal for almost 15 years. Only post 2015 have they been re-visited. Aug 13 at 23:02