I'm using the gbm.step package in R to look at the influence of three continuous variables on my continuous response variable. I have 234 observations. The model:

poa.tc2.lr005.bg0.5 <- gbm.step(data=poa, 
                          gbm.x = 8:10, 
                          gbm.y = 7,    
                          tree.complexity = 2, 
                          family = "gaussian", 
                          #n.trees = 50,
                          #n.folds = 10,       
                          #step.size = 25, 
                          max.trees = 10000,    
                          prev.stratify = FALSE, 
                          learning.rate = 0.005, 
                          bag.fraction = 0.5)

After settling on some initial parameters of tc, lr and bag fraction, I would like to produce partial dependency plots for my predictor variables using gbm.plot:

gbm.plot(poa.tc2.lr005.bg0.5, n.plots=3,
         write.title = F, 
         y.label="Marginal effect on gs")

The general trend of the dependencies makes sense given the data, however the y-axis scale is throwing me off.

The y-axis range is between -0.04 to 0.06, however, my response variable range is between 0.02 and 0.38. With the new range, it's hard to accurately interpret the results.

My questions:

  1. Why doesn't the y-axis in the dependency plot reflect the range of my dependent variable? Are the values normalized to something?

  2. How do I extract the values used to construct these graphs? I would like to reconstruct the dependency plots in a different program. I have tried

    poa.fitted <- (poa.tc2.lr005.bg0.5$fitted)

but those values are not the same used in the dependency plots generated by the gbm.plot code above. Is there a different output for these values that I should be looking for?


3 Answers 3


First of all, I think you're using the gbm.step function from the dismo package.

Even though your y ranges from 0.02 to 0.38, the model can still decide that certain variables (or certain ranges of a variable) have a negative contribution to y's value. If this is the case, the marginal effects plot will include negative values.

Finally, use the plot.gbm function from the gbm package to get the values used for the marginal dependency plots:

model <- gbm(Sepal.Width  ~ ., distribution = "gaussian", data=iris)
plot(model, 'Species', return.grid=TRUE)
     Species        y
1     setosa 3.067434
2 versicolor 3.052529
3  virginica 3.052915

I think this plot uses a slightly different approach than dismo::gbm.plot, but the results should be similar.

  • $\begingroup$ thanks for catching my error in the first line. The plot.gbm function returned a list that is different from the values produced by my original gbm.plot function. From my understanding, these new values are the modeled values of the response variable plotted against an even distribution of the predictor, after accounting for the average effects from the other predictors. Am I on the right track? Any idea how to retrieve the original values from the gbm.plot function? Thanks for your help. $\endgroup$ Commented Nov 5, 2014 at 18:10
  • $\begingroup$ @woodland_creature Yes, the new values are the modeled values of the response variable, after accounting for the average effects of the other variables. I'm not super familiar with the gbm.plot function but it doesn't look like there's an option to return the plot data. It looks like the function's written entirely in R, so you could copy the source code and modify it to return the plot data. Type: gbm.plot into your console. $\endgroup$
    – Zach
    Commented Nov 5, 2014 at 19:06
  • $\begingroup$ Great, thanks! Using the modeled values, does this change the interpretation of the plot at all? In other words, when the y-axis contained negative values, it was easy to conclude that high predictor values had a negative effect on the response (I have a negative relationship). Now that I have the modeled values and the same negative relationship, can I safely interpret this plot to mean "high values of the predictor variable will result in lower values of the response, given the model"? Also, I'll try messing with the source code in gbm.plot, but my R skills will likely come up short. $\endgroup$ Commented Nov 5, 2014 at 19:43
  • $\begingroup$ @woodland_creature: The plot is basically showing you the model's view of the world, so yes, a negative region of the plot indicates the model thinks values in that region will decrease the response. If my answer is sufficient, please feel free to up-vote it and/or mark it as accepted. $\endgroup$
    – Zach
    Commented Nov 5, 2014 at 21:58
  • $\begingroup$ I know this is an old post but what does it actually mean to have a "negative contribution to y"? Does it mean that the values of x that made y negative actually tend to decrease values for dependent variable towards it's lower end? $\endgroup$ Commented Sep 16, 2017 at 3:11

The accepted answer is correct. However, I wanted to elaborate on some of the comments to the answer.

dismo::gbm.plot essentially calls gbm:plot.gbm in it's source code. The grid returned from plot.gbm is the same grid gbm.plot uses to make it's plots. dismo::gbm.plot just adds additional options to gbm:plot.gbm. If you want to recreate your own plots use gbm:plot.gbm with return.grid=TRUE to return a grid of 100 predictors and responses.

If you want to transform the responses like in dismo::gbm.plot or like in Elith et al. 2008 you can run functions on the grid.

#gbm.object= your gbm object
#j = variable number

k <- match(gbm.object$contributions$var[j],gbm.object$gbm.call$predictor.names)

response.matrix <- gbm::plot.gbm(gbm.object, k, return.grid = TRUE)
predictors[[j]] <- response.matrix[,1]   

#If factor variable
predictors[[j]] <- factor(predictors[[j]],levels = levels(gbm.object$gbm.call$dataframe[,gbm.object$gbm.call$gbm.x[k]]))

responses[[j]] <- response.matrix[,2] - mean(response.matrix[,2])

#Perform transformations
#Put it on a log scale

#Center the response to have zero mean over the data distribution
  responses[[j]]<-scale(responses[[j]], scale = FALSE)

#Add a smoothed line to an active plot
loess(responses[[j]] ~ predictors[[j]], span = 0.3)
  • $\begingroup$ "transform responses like dismo::gbm.plot" - gbm.plot's only modification to the plot.gbm data is the optional 'common.scale'. For Elith 2008 working guide, Fig 6: "Y axes are on the logit scale and are centred to have zero mean over the data distribution" so you may well be right. Though I'm not sure the benefit of this? Trying with dummy data, this centres values around 0, Inverts their position relative to 0 (top becomes bottom), Exponentiates them (midrange values push towards extremes), Then rescales from 0:1 to +/- values by subtracting the mean from each value. $\endgroup$
    – dez93_2000
    Commented Oct 18, 2021 at 18:16

A very simple way of getting the values used for plotting: 1. Type plot.gbm into your console. 2. Copy the function code to an empty script file. 3. On the first row type gbm.plot2 <- in front of function(..., this defines a new function gbm.plot2() that produces the same output as the original function. 4. Add print(X) before the very last }. 5. Call the new function temp = gbm.plot2(your.gbm.model, 1) and the data is saved as temp.

  • $\begingroup$ As I understand it this will fall foul of the central issue: gbm.plot produces the line files but doesn't easily allow extraction of values. Plot.gbm already allows easy value extraction (return.grid=TRUE) but the values are different. $\endgroup$
    – dez93_2000
    Commented Apr 28, 2018 at 14:36

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.