# Partial dependence plot, GBM multinomial

I'm using gbm package to create a gradient boosting model with multinomial family. I have problems in plotting and interpreting partial dependence plots with this type of models.

an example:

library(gbm)
data(iris)
GBM.model = gbm(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, data = iris,distribution = "multinomial")
plot.gbm(GBM.model, "Petal.Length")


I get the following plot: Could you help me to improve the plot and interpret it? I think it should have a line for each category of the response variable but I'm not sure. Is there another way of getting partial dependence plots for multinomial gbm models? I would appreciate any help with this

• Welcome to CV.SE. Your intuition that we should have a line for each category of the response variable is correct. (+1) Please see my answer below for more details. Sep 27, 2020 at 12:46

I think it is worth mentioning here that currently (Sep 2020) the gbm clear warns: Setting 'distribution = "multinomial"' is ill-advised as it is currently broken. It exists only for backwards compatibility. Use at your own risk.
In any case, we can get partial dependency plot through the package pdp. For example: partial(GBM.model, "Petal.Length", n.trees=100, which.class = 2, plot = TRUE, prob = TRUE)). Please note the which.class argument. If we do not supply it, we get only for class 1 (setosa).
To that extent, based on the graph you have attached we can tell you what is going on (wrong) with the gbm.plot. The function appears to have plotted all three class-specific partial dependency plots at the same graph and having connected the end of each class-specific PDP with the start of the next class' PDP. That is why we have these long straight lines from (7,-10) to (1,0) and then from (-7,-2) to (-3,1); we connect the PDPs. Finally please note that when plotting PDPs for classification it is usually more informative to plot the PDPs on the probability (i.e. the response) scale. The (-10,10) scale of the y-axis corresponds to the value on the logistic scale and while informative it is less obvious what is the impact on the final classification.
• If we want to have more customised plot I think the easier path is to pass the output of partial to an object without printing it directly (i.e. have plot = FALSE) and then using standard R functionality through legend to put the legend we would need, etc. If this answer is helpful please consider upvoting it and if it resolves your question marking it as the accepted answer. Sep 27, 2020 at 17:06