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
 A: 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.
