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

Interfaces brms
1

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,
              family=Gamma(link="log"),
              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,
              family=brmsfamily(family="Gamma", link="log", link_shape="log"),
              prior=prior(normal(0,2),class="Intercept"),
              backend="cmdstanr",
              control = list(adapt_delta=0.99))

# 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
2

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

3

@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
4

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,
              family=brmsfamily(family="Gamma", link="log", link_shape="log"),
              prior=prior(normal(0,2),class="Intercept"),
              backend="cmdstanr",
              control = list(adapt_delta=0.99))

I.e. we model the mean, then 1/mean is rate/shape and then the code multiplies by shape and we provide the correct thing to the likelihood function. A sanity check on whether I finally got it right would be appreciated, apologies for the confusion.

5

Yes, exactly. This is how brms intends to do it and from the code I am confident it does what is intended.