Predict or parametric bootstrap for a generalized mixed model

I want to fit a generalized mixed model to a dataset where the dependent variable is velocity (continuous positive numbers), the random variable is the subject ID, and the predictor is distance to some point (also a continuous positive number).
Here’s a contrived MWE:

using MixedModels, DataFrames, Distributions

n = 1000 # number of obs
x = 100rand(n) # predictor
y = 5rand(Gamma(), n) .+ rand(n).*x/100 # dependent 
data = DataFrame(x = x, y = y, id = rand(1:5, n)) # the data frame
categorical!(data, :id) # make the random variable a categorical one

m = fit(MixedModel, @formula(y ~ 1 + x + (1|id)), data, Gamma()) # fit the model

Apart from the P value, I’d also like to plot the confidence interval area around the fit, something like this (image arbitrarily stolen from here):

How do I get those bounds from a Generalized Mixed Model like m?

I think @dmbates has implemented parametric bootstrapping, at least for LinearMixedModel but maybe not for GeneralizedLinearMixedModel:

1 Like

Yes (docs are here), but I want it for a Gamma link function. So my only options (after some discussion on Slack) are to just use a Gaussian link function instead (thus using a LinearMixedModel) and/or removing the random effect by calculating the mean within subject ID.