3
$\begingroup$

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:

Comparing the two approaches

x1*x2 surface from the partial dependence plot approach:

Using the partial dependence plot approach

x1*x2 surface from the predict() approach:

Using the predict() approach

$\endgroup$

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$
    – mkt
    Aug 13 at 10:52
  • 1
    $\begingroup$ 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. $\endgroup$
    – usεr11852
    Aug 13 at 20:37
  • $\begingroup$ 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. $\endgroup$
    – mkt
    Aug 13 at 22:05
  • 1
    $\begingroup$ 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. $\endgroup$
    – usεr11852
    Aug 13 at 23:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.