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`

?