partial regression plots in natural scale (Poisson)

[nicolashazzi5]

Hi everyone!
I would like to ask about partial regression plots using poisson or nb distribution using the draw function.

I am currently using the draw() function as follows:
draw(m1, select = "s(month)", residuals = TRUE)
to plot the partial regression plots of my GAM model. However, I would like to visualize the Y-axis in its natural (response) scale, as my model uses a negative binomial distribution.

I was also wondering if it is possible to shift the scale so that the intercept is included, making the plot more interpretable. I know that in mgcv this can be done using:
plot(gam_model, seWithMean = TRUE, shift = coef(gam_model)[1])
Is there a similar option to achieve this in gratia?

Best regards,
Nicolás Hazzi

[gavinsimpson]

Note that adding just the model constant term is not sufficient to put the partial effect plot on to the natural response scale in the case of models, like yours, that involve a non-identity link function.

You can achieve what you want two ways with gratia:

1. Read ?draw.gam and note the arguments fun and constant. In that case, your call would be:

   draw(
     m1,
     select = "s(month)",
     residuals = TRUE,
     fun = inv_link(m1),           # extract the inverse of model's link function
     constant = model_constant(m1) # extract the model constant term (Intercept)
   )
   

2. A more general approach is to use conditional estimates from the model itself. These are predictions from the model at some values of the covariate(s) of interest and may include or exclude terms from the model as needed, or these can be held constant. For example:

   ds <- m1 |> data_slice(month = evenly(month))
   
   m1 |>
     conditional_values(condition = "month") |>
     draw()
   

   will produce a plot of estimates conditioned on 100 evenly spaced values of month while holding any other covariates at their median or modal values in the training data for continuous or factor variables respectively. You can add exclude = "s(foo)" to the call to exclude one or more model terms; this is useful if you want to exclude the effects of any random terms from the predictions, to give population or subject-average estimates depending on whether the model uses an identity link function or not, respectively. You can then add the data to the plot use a point layer with & geom_point(data = bar, aes(x = month, y = Y)) for example where bar is the data frame used to fit the model and Y is the response variable in the model.