Why don't partial dependence plots match model predictions?

Question

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:

Answer 1

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.