# Gamma regression with known shape parameter - brms ignores offsets for shape in predictions?

Sometimes we know the shape parameter, e.g. when we look at the sampling distribution for the variance (of a normal r.v.) from a study.

My first attempt at doing this didn’t work, but I got this to work with distributional regression using ` bf( precision ~ 1 + (1|study) + offset(log(known_shape)), shape ~ 0 + offset(log(known_shape))`, where as far as I can tell no actual model is fit for shape and it’s just being treated as a constant (as per the Stan code I inspected). However, it looks like the offset for the shape is being ignored when I use the `predict` method.

Have I made some obvious mistake or is there an easier way to do this?

``````library(tidyverse)
library(brms)

varests = tibble(vest = c(1, 1.2, 1.1, 0.7),
n = c(50, 48, 100, 10),
study=factor(1:4)) %>%
mutate(precision = 1/vest,
nu = n-1,
tvest = 1/(vest*nu/2))

# Note for estimated variances, we have
#   vest * nu /2 ~ Gamma(shape = nu/2, rate = 1/sigma^2)
# So, to model variances across studies, we want to model something like
#   expected -log(variance) = 1 + (1|study)
#    or (because there's no negative-log-link function):
#   expected log(precision) = 1 + (1|study)
# and fix shape = nu/2.

# Try 1: This does not work, apparently a constant prior defined in terms of
# a variable is not possible like this.
brmfit1 = brm(tvest ~ 1 + (1|study),
data = varests,
prior=prior(normal(0,2),class="Intercept") +
prior(constant(nu/2),class="shape"),
backend="cmdstanr"
)

# Try 2: this samples and by looking at the Stan code via stancode(brmfit2)
# it seems to do what I want.
brmfit2 = brm(bf( precision ~ 1 + (1|study) + offset(log(nu/2)),
shape ~ 0 + offset(log(nu/2)) ),
data = varests,
prior=prior(normal(0,2),class="Intercept"),
backend="cmdstanr",

# However, the predictions I get, here seem wrong:
predict(brmfit2, newdata=tibble(nu=100, study=1L))
# Gets me about 2.2 (as median)
predict(brmfit2, newdata=tibble(nu=10, study=1L))
# Gets me about 0.22 (as median)
# This seems to be off by exactly what you'd expect when the offset in the shape
# parameter is ommitted.
``````
• Operating System: Windows 10
• brms Version: 2.16.1

I checked via debug and at least in brms 2.17+ the offset seems to be included correctly in the predictions.

@paul.buerkner , I hope I’ve got mixed up somehow, but if not, I think I found what the issue is. It’s in the following bits in the Stan code:

``````...
// initialize linear predictor term
vector[N] mu = Intercept + rep_vector(0.0, N) + offsets;
...
for (n in 1:N) {
// add more terms to the linear predictor
mu[n] += r_1_1[J_1[n]] * Z_1_1[n];
}
...
for (n in 1:N) {
// apply the inverse link function
mu[n] = shape[n] * exp(-(mu[n]));
}
target += gamma_lpdf(Y | shape, mu);
``````

It seems to me that `mu[n]` ends up containing the expected value (i.e. shape/rate), but the `gamma_lpdf` function expects shape and rate as its parameters, but instead we are giving it shape and shape/rate. This is then correctly reflected in `posterior_predict_gamma` (obviously, important that the two match up). Thus, if I re-predict the same data that I modelled (i.e. not changing the offsets, e.g. `predict(brmfit2, newdata=varests)`), I get sensible looking predictions back, but once I predict for new data with different offfsets, the answers end up being nonsense, I think.

1 Like

Ouch, did I get this totally wrong and everything is fine??? Ìf the regression model is for the mean and things should just work the way I want, if I use:

``````brmfit3 = brm(bf( vest ~ 1 + (1|study),
shape ~ 0 + offset(log(nu/2)) ),
data = varests,