Terrible sampling w/gamma regression

Modeling
,
1

Hi all,

Trying to run the following model:

data{
vector[126] LPS_concentration;
int Plate[126];
int Time2[126];
int Time[126];
int Subject[126];
int LPS[126];
}
parameters{
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;
}
model{
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!

2

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)?

3

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
4

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.

5

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.

6

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
7

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.

2 Likes
8

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.

3 Likes