Behaviour of draw() with a model offset

[DEHewitt]

Hi @gavinsimpson,

Not really an issue per se rather just trying to clarify something that I haven't been able to find an answer to. If I run the following model:

m <- gam(y ~ s(x0) + s(x1), family = nb(), data = data)
gratia::draw(m, select = "s(x0)")

My understanding is that the resultant plot provides the marginal effects of x0 on y - that is the effect of x0 averaged over the values of x1. Is the same true if I run the following:

m <- gam(y ~ s(x0) + s(x1) + offset(log(z)), family = nb(link = "log"), data = data)
gratia::draw(m, select = "s(x0)")

Except now it is the effect of x0 averaged over the values of x1 and z? Is it possible to set z to some reference level? In my case z refers to fishing effort and I am interested in presenting plots that show the estimates of x0, x1, etc. on harvest given 1 unit of effort (i.e., z = 1).

Thanks for your help, and for all the effort you go to answering GAM-related questions across the web (e.g., StackOverflow). Between that and your blog you've taught me most of what I know!

[gavinsimpson]

I would describe the plots for the first model as showing the partial effect of x0 say on y conditional upon the effect of x1 being set equal to 0, on the scale of the linear predictor (link scale; in this case the log scale).

For the second model, I would say the plots show the same thing as for the first model, except that now the response is not in units of y but in units of y per z. We interpret the expected value of y, $\mathbb{E}(y_i) = \lambda_i$, as a rate, with $lambda_i$ giving the expected count for observation $i$. With an offset, this becomes the expected rate of y per unit z. If z was population in each region and y was the number of deaths recorded per region, the model with z as an offset would be for the number of deaths per head of population in each region. As such, the plot of s(x0) would show the partial effect of x0 on the rate pf y per z.

I'm a bit rusty in the terminology but IIRC these are also marginal effects. The plots certainly seem to correspond to the slopes produced by the marginaleffects package (which is what it calls "marginal effects")

There is such a thing as an average marginal effect, which is the average of the marginal effects (slope) of each data point.

[gavinsimpson]

As for

Is it possible to set z to some reference level? In my case z refers to fishing effort and I am interested in presenting plots that show the estimates of x0, x1, etc. on harvest given 1 unit of effort (i.e., z = 1).

no, it is not possible to do this with gratia; if you want that, use data_slice to generate data over x0 holding x1 at some value (just don't include it in the slice and the function will do this for you), while providing a value for z. Then pass the model and the data slice to fitted_values(). If you were to also add exclude = "s(x1)" to the fitted_values() call and scale = "link" you'll get something close to what you want, and if you used terms = "s(x0)", scale = "link" you should get exactly what you want.