Terrible sampling w/gamma regression

Hi all,

Trying to run the following model:

vector[126] LPS_concentration;
int Plate[126];
int Time2[126];
int Time[126];
int Subject[126];
int LPS[126];
vector<lower=0>[2] a;
vector[18] a_Subject;
vector[2] b_Time;
vector[2] b_Time2;
vector[2] a_Plate;
real<lower=0> scale;
vector[126] mu;
scale ~ exponential( 20.41);
a_Plate ~ normal( 0 , 0.1 );
b_Time2 ~ normal( 0 , 0.001 );
b_Time ~ normal( -0.012 , 0.01 );
a_Subject ~ normal( 0 , 0.237 );
a ~ gamma( 1.248/0.029 , 1/0.029 );
for ( i in 1:126 ) {
mu[i] = a[LPS[i]] + a_Subject[Subject[i]] + b_Time[LPS[i]] * Time[i] + b_Time2[LPS[i]] * Time2[i] + a_Plate[Plate[i]];
mu[i] = exp(mu[i]);
LPS_concentration ~ gamma( mu/scale , 1/scale );

Have ran up to 100,000 iterations. These priors come from preliminary data, but I suspect they may be part of the problem…any help greatly appreciated!

I assume you’ve done prior predictive checks and also used fake data to see that you can recover the parameters using your model? If that’s the case then,

just a thought, in the stats pkg gamma uses inverse link as a link function, while in e.g. brms it’s the log link function. I don’t know why, but perhaps Paul used log link for a reason? @scholz and @Solomon could perhaps weigh in here (@Solomon you played with gamma recently for a another purpose)?

The inverse is the canonical link for the gamma but brms uses the log as default for all lower-bound distributions. That might explain the difference.

Could you expand on the problem? Does it take forever to sample or are the samples not useful / result in bad diagnostics.
So in addition to prior predictive checks and sbc, I would look at trace plots, and diagnostics like \widehat{R} and ESS to see if there is a specific part of your model that causes the problems.
In addition, some of the priors look rather narrow to me, which might be supported by preliminary data but relaxing them a little could help sampling if, for some reason, the preliminary data differs from your actual data.

1 Like

My thoughts will be very limited when it comes to raw Stan. But I have done some gamma regression with brm() and have found the log link samples much better than the inverse link.

It samples slowly, yes, but not terribly so. More concerning is that the Rhat and n_eff are also pretty bad. Some parameters okay, but many have Rhat between 1.1 and 1.3 with n_eff < 50, and the ‘scale’ parameter is consistently bad with a Rhat of 17 and n_eff of 2.

I have not ran it on fake data, but it does basically produce results that are consistent with simpler versions of the model that did sample well.

That makes me wonder if it’d help if you modeled log(scale) instead. Then instead of an exponential prior, you could use something like normal(log(1/20.41), 1).

1 Like

I agree with Solomons points. The prior for the scale parameter puts most weight below 0.1, is that consistent with what you expect? Maybe relax it a little to better explore the parameter space.


That is what I expect, but I nevertheless loosened up the prior w/ dexp(2). This sampled well, and indeed the scale is < 0.1.