Documentation for gratia::posterior_samples() behaviour with offsets

[jonathonmellor]

Hi Gavin,

Really appreciate the work in this package - hoping to move our code to using it as much as possible.

Would it be possible to document (or happy to raise a PR if it's explained here) how gratia::posterior_samples behaves with predictions from models using an offset?

I'm unclear on whether the offset is ignored, or whether the offset is used to convert back to, for example, counts, before sampling from the error distribution?

[gavinsimpson]

Thanks for the question @jonathonmellor. What happens will depend on how you include the offset term in the model. This is because, as mentioned in ?mgcv::gam, whether or not predict.gam() ignores the offset depends on whether you included it via the offset() function in the model formula or via the offset argument to gam() etc.

posterior_samples() works by first calling fitted_samples() to get posterior expected values and then the relevant random number generator to generate new response values for those posterior expectations. fitted_samples() relies on predict.gam() (or predict.bam(), etc) and hence it and predicted_samples() and posterior_samples() will all be sensitive to how the offset term was included in the model.

If you included the offset through the formula, predict.gam() will require the offset variable to be present in the data at which you are trying to predict, so if you pass something to data it must contain the variable for the offset.

Basically, posterior_samples() is just a series of wrappers around predict.gam() and co, so any behaviour those functions have is inherited by posterior_samples().

[jonathonmellor]

for a negbinom example, if Xp doesn't have the offset included then it would be something like

sims <- Xp %*% t(betas + log(denom)) though ideally this wouldn't need to be specified.

[gavinsimpson]

Ah, right; I was forgetting that when type = "lpmatrix", which is how fitted_samples() uses predict.gam() it doesn't quite work how I described. In such instances, predict.gam() returns the Xp matrix with an attribute model.offset containing the vector of values for the offset term:

> Xp <- predict(mod, type = "lpmatrix")
> head(attr(Xp, "model.offset"))                                              
[1] 6.907755 6.906755 6.907755 6.908755 6.907755 6.908755

So, my initial thought is that the correct behaviour here would be to add this attribute to the matrix Xp and add a constant coefficient of 1 to betas for this column, or modify the line

sims <- Xp %*% t(betas)

to be

sims <- (Xp %*% t(betas)) + attr(Xp, "model.offset")

[jonathonmellor]

Oh interesting!

Option 2 feels more explicit in behaviour rather than fixing coefficients? I suppose it depends if you want to change how post_draws() and it's samplers work.

Happy to raise any PRs etc.

[gavinsimpson]

I believe this is now implemented via dd93df2. Would you mind checking that it works as per your expectations - it seemed to do the right thing on the example here but I haven't tested it widely.

[jonathonmellor]

Looks like it's working on both the toy example and on our real world data models as well. Thank you so much, really appreciate the quick update!

[gavinsimpson]

Great! Thanks again for raising the issue here.