2

I'm trying to do a bayesian gamma regression with stan. I know the correct link function is the inverse canonical link, but if i dont use a log link parameters can be negative, and enter in a gamma distribution with a negative value, that obviously can't be possible. how can i deal with it?

parameters {
vector[K] betas; //the regression parameters
real beta0;
real<lower=0.001, upper=100 > shape; //the variance parameter
}
transformed parameters {
vector<lower=0>[N] eta; //the expected values (linear predictor)
vector[N] alpha; //shape parameter for the gamma distribution
vector[N] beta; //rate parameter for the gamma distribution

eta <- beta0 + X*betas; //using the log link 
}
model {  
beta0 ~ normal( 0 , 2^2 ); 

for(i in 2:K)
betas[i] ~ normal( 0 , 2^2 ); 

y ~ gamma(shape,shape * eta);
}
Jacopo Soppelsa
  • 31
  • 3

1 Answers1

0

I was struggling with this a couple weeks ago, and while I don't consider this a definitive answer, hopefully it is still helpful. For what it's worth, McCullagh and Nelder directly acknowledge this inappropriate support of the canonical link function. They advise that one must constrain the betas to properly match the support. Here's the relevant passage

The canonical link function yields sufficient statistics which are linear functions of the data and it is given by η = 1/μ. Unlike the canonical links for the Poisson and binomial distributions, the reciprocal transformation, which is often interpretable as the rate of a process, does not map the range of μ onto the whole real line. Thus the requirement that η > 0 implies restrictions on the βs in any linear model. Suitable precautions must be taken in computing β_hat so that negative values of μ_hat are avoided.

-- McCullagh and Nelder (1989). Generalized Linear Models. p. 291

It depends on your X values, but as far as I can tell (please correct me someone!) in an MCMC-based Bayesian case, you can achieve this by either using a truncated prior on the betas or a strong enough prior on your intercept to make the inappropriate regions numerically impossible to reach.

In my case, I ultimately used an identity link with a strong positive prior intercept and that was sufficient and yielded reasonable results.

Also, the choice of link really depends on your X. As the passage above implies, the use of the canonical link assumes that your linear model is in rate space. Using log or identity link functions appear to be also very common, and ultimately it's about providing a space that offers a sufficient span for the linear function to capture the response.

merv
  • 67,214
  • 13
  • 180
  • 245
  • Thanks a lot, it worked using a more informative prior for beta0, the problem using a log link is that my professor said it isn't the true link for a gamma distribution because the link is the inverse, and she doesn't want to use a log. I don't know if I could use a log link when an inverse is required. By the way using a prior that force betas to be positive is very restrictive. My x are all positive values and dummies – Jacopo Soppelsa Dec 13 '18 at 15:48
  • @JacopoSoppelsa glad that worked. I agree, one should not blindly force `beta`s to be positive, and specifically not if all your `X` values are nonnegative. To clarify, I was only stating that you could truncate them (or a strong prior would also work), but did not specify to what range since again that depends on `X`. If your professor wants a specific direct-interpretability of the beta values in rate space, then that would make sense to insist on using that link. I wouldn't go so far as to call it "the true link" though. – merv Dec 13 '18 at 16:16
  • thanks a lot, it helped me so much. Last thing, are you saying that using a canonical link let me interpret betas at a linear regression but I can use a not canonical link but there are interpretabily problem? – Jacopo Soppelsa Dec 13 '18 at 16:21
  • @JacopoSoppelsa only that the interpretation is more direct. For example in purely categorical predictor space (with 0,1 dummies), then the betas values are directly interpretable as rate effects. That is, the sizes of the betas have meaning in rate space. If you're doing survival time analysis that could be very informative. If all you care about are whether the betas are significant and would be satisfied with *relative* effect sizes (e.g., ranking your predictors), then I think you are little more free to use other links. I am not an expert though (just a PhD student). – merv Dec 13 '18 at 16:39
  • @JacopoSoppelsa You may want to also try asking questions like this on [CrossValidated](https://stats.stackexchange.com), as there will be many more experienced statisticians there to offer insights. – merv Dec 13 '18 at 16:41