Gamma regression in Stan vs Frequentist approach

Modeling
11

First, the original post didn’t constrain phi to be > 0 so there’s tons of problems just with that. Add real<lower=0> phi. The cauchy priors probably don’t help and Gelman doesn’t recommend those anymore. Try something like normal(0, 4) that is reasonably diffuse.

An issue with the gamma distribution shows correlation with the parameters in the fisher information matrix. See Posterior estimates of rate and shape of gamma distribution are dependent.

The given parameterization from @Max_Mantei and the one in the post I linked to above give 0s for the off diagonal part of the Fisher information matrix which will help sampling.

Btw here’s what I get with the OP’s data and this slightly updated model with no additional reparameterization of the gamma distribution:

  variable    mean  median      sd     mad      q5     q95  rhat ess_bulk ess_tail
  <chr>      <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
1 betas[1]  0.178   0.178  0.00419 0.00411  0.170   0.184   1.00    1818.    2392.
2 betas[2]  0.0471  0.0471 0.00378 0.00379  0.0410  0.0534  1.00    4319.    2902.
3 betas[3] -0.0693 -0.0693 0.00465 0.00472 -0.0769 -0.0615  1.00    3408.    2427.
4 betas[4] -0.0228 -0.0228 0.00454 0.00457 -0.0304 -0.0155  1.00    3950.    3139.
5 betas[5] 13.8    13.8    0.0203  0.0200  13.8    13.8     1.00    1673.    2336.

data {
  int N; //the number of observations
  int K; //the number of columns in the model matrix
  vector[N] y; //the response
  matrix[N,K] X; //the model matrix
}
parameters {
  vector[K] betas; //the regression parameters
  real<lower=0> phi; //the variance parameter
}
transformed parameters {
  vector[N] mu = exp(X * betas); //the expected values (linear predictor)
  vector[N] alpha = mu .* mu / phi; //shape parameter for the gamma distribution
  vector[N] beta = mu / phi; //rate parameter for the gamma distribution
}
model {  
  betas ~ normal(0, 4);

  y ~ gamma(alpha,beta);
}