Trouble recreating the stan_glm for Gamma with Inverse Link

Modeling
1

Hi all! I’m relatively new to rstan, and for an Epidemiology project I need to run a Bayesian Gamma regression, with inverse link. I have gotten it to work already within the rstanarm library, but I tried to recreate it myself to understand it better and to get used to the stan language; and I’ve had some issues.

Just to clarify what the inverse link entails, we have y \sim Gamma(\alpha, \beta), with \alpha being the shape parameter, and \beta the rate. The inverse link is \mu^{-1} = \beta_0 + \beta_1 x_1 + \beta_2 x_2 = \eta. If we reparametrize with \mu = \alpha / \beta, we can write y \sim Gamma(\alpha, \eta\alpha).

My issue is as follows: because both shape and rate need to be strictly positive, the linear predictor \eta should also be positive. Despite running the below code (constraining it to be positive), I run into lots of divergent transitions (like, hundreds) and chains that are all over the place.

In r:

# Simulate some data
x1 <- 1:150
x2 <- x1^2
b <- c(0, 0.05, -0.00001)
eta <- (b[1] + b[2] * x1 + b[3] * x2)
y <- rgamma(150, shape = 0.7, rate = 0.7 * eta)

With stan:

data {
  int<lower=0> N;
  vector<lower=0>[N] y;
  vector[N] x1;
  vector[N] x2;
}

parameters {
  real b0;
  real b1;
  real b2;
  real<lower=0.001,upper=100 > shape;
}

transformed parameters{
  vector<lower=0>[N] eta; // linear predictor
  eta = b0 + b1 * x1 + b2 * x2;
}

model {
  b0 ~ normal(0, 50);
  b1 ~ normal(0, 50);
  b2 ~ normal(0, 50);
  shape ~ uniform(0.01, 50);
  for(i in 1:N){
      y[i] ~ gamma(shape, shape * eta[i]); 
  }
}

My viewer is full messages such as:
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: validate transformed params: eta[1] is -3.50457, but must be greater than or equal to 0 (in ‘model865c2a234c78_gamm’ at line 18)

Any ideas on what I am doing wrong here? What does rstanarm do differently? (I have looked at some of the source code, but I am lost to be honest) Thanks in advance!

NB: Any code I have found so far deals with the log link, which is substantially nicer since \eta is by default positive!
(See: http://seananderson.ca/2014/04/08/gamma-glms/)

2

The inverse link is real bad. In rstanarm, we have to work around this by defining b0_raw as positive, and then defining b0 in transformed parameters as

transformed parameters {
  vector<lower=0>[N] eta = b1 * x1 + b2 * x2;
  real b0 = b0_raw - min(eta);
  eta = eta + b0;
}

Even still, you need to either get a lot of divergences or have to take really small steps. And the priors become basically uninterpretable.

3

Thanks for the quick answer! I was curious since I’ve also been trying to implement it in JAGS, also to no avail.

And the main motivation for pursuing Bayesian here was that the inverse link is also a pain in the normal glm, where with non-simulated data it doesn’t converge very easily!