8
$\begingroup$

I noticed there's a difference in partial dependence calculated by R package gbm and Python's scikit-learn.

Here's gbm's partial dependence of median value on median income of the California housing dataset:

enter image description here

And here's scikit-learn': enter image description here

It's easy to see that R's partial dependence ranges from 1.5 to 4.5, whereas scikit-learn's from -0.5 to 1.5, but the shape of the lines is nearly the same. I don't get why it's like that.

Relevant code:

R

library(oem)
library(gbm)
data(calHousing)
X <- calHousing[ ,!(colnames(calHousing) == "medianValue")]
y <- calHousing$medianValue / 100000
gbm.model <- gbm.fit(X, y, distribution="gaussian", n.trees=100, interaction.depth=12, shrinkage=0.15)
plot(gbm.model, i.var="medianIncome")

Python code is a copy-paste from scikit-learn' example page.

print(__doc__)

import numpy as np
import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

from sklearn.cross_validation import train_test_split
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.ensemble.partial_dependence import plot_partial_dependence
from sklearn.ensemble.partial_dependence import partial_dependence
from sklearn.datasets.california_housing import fetch_california_housing


def main():
    cal_housing = fetch_california_housing()
    # import ipdb; ipdb.set_trace()
    # split 80/20 train-test
    X_train, X_test, y_train, y_test = train_test_split(cal_housing.data,
                                                        cal_housing.target,
                                                        test_size=0.2,
                                                        random_state=1)
    names = cal_housing.feature_names

    print('_' * 80)
    print("Training GBRT...")
    clf = GradientBoostingRegressor(n_estimators=100, max_depth=12, min_samples_split=10,
                                    learning_rate=0.15, loss='ls', subsample=0.5,
                                    random_state=1)
    clf.fit(X_train, y_train)
    print("done.")

    print('_' * 80)
    print('Convenience plot with ``partial_dependence_plots``')
    print

    features = [0, 5, 1, 2, (5, 1)]
    fig, axs = plot_partial_dependence(clf, X_train, features,
                                       feature_names=names,
                                       n_jobs=3, grid_resolution=100)
    fig.suptitle('Partial dependence of house value on nonlocation features\n'
                 'for the California housing dataset')
    plt.subplots_adjust(top=0.9)  # tight_layout causes overlap with suptitle
    plt.show()


# Needed on Windows because plot_partial_dependence uses multiprocessing
if __name__ == '__main__':
    main()
$\endgroup$
2
  • 1
    $\begingroup$ My guess is some sort of standardization took place for the outcome variable for scikit-learn (either in the underlying data or for the plotting function) that didn't take place in the R example. $\endgroup$
    – dmartin
    Commented Apr 2, 2016 at 20:40
  • $\begingroup$ Max depth and interaction depth do not mean the same thing. $\endgroup$ Commented Aug 27, 2016 at 18:54

1 Answer 1

9
$\begingroup$

Scikit-learn center the partial dependence with the mean of the target value, R does not.

Here is an example using diabetes dataset.

R

data(diabetes, package="lars")

y        <- diabetes$y
x        <- diabetes$x
class(x) <- "matrix"
data     <- data.frame(y, as.data.frame(x))

model <- gbm::gbm(formula = y ~ . , data = data, distribution = "gaussian", 
                  shrinkage = 1, bag.fraction = 1, n.trees = 100,
                  interaction.depth = 2, verbose = T, keep.data = F)

partial <- plot.gbm(dgbm, i.var = 1, return.grid = T)
plot(partial[, 2] - mean(y), type = "l")

R partial plot

Python

import numpy as np
import sklearn
import matplotlib.pyplot as plt
import sklearn.datasets
import sklearn.ensemble
from sklearn.ensemble.partial_dependence import partial_dependence

diabetes = sklearn.datasets.load_diabetes()
X= diabetes.data
y= diabetes.target



gbm = sklearn.ensemble.GradientBoostingRegressor(loss='ls', learning_rate=1, max_leaf_nodes=3, min_samples_leaf=10,
                                             n_estimators=100, verbose=True)
model_gbm = gbm.fit(X, y)

partial, axe = partial_dependence(gbrt=model_gbm, X=X, target_variables=(0))

plt.plot(partial.T)

enter image description here

$\endgroup$
2
  • $\begingroup$ Great answer, I wish I could upvote it twice. Any idea why Scikit-learn would subtract off the mean of the target variable? Isn't is nice to see, for a given independent variable, what the predicted target variable is after averaging over your data? $\endgroup$
    – Chris A.
    Commented Aug 25, 2016 at 15:49
  • $\begingroup$ @EtienneKintzler, how would the scale on y-axis get affected in case of classification. $\endgroup$ Commented Apr 21, 2022 at 10:12

Your Answer

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

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