Hierarchical Model - Hyper priors for a Gamma Distribution

Modeling
1

I am working on building a hierarchical RL model for a parameter drawn from a gamma distribution (beta_1 variable). I am attempting to re-parameterize the gamma distribution in terms of its mean and variance and then set the hyperpriors from a uniform (0,20) but am getting an insane number of divergences so am definitely just doing something wrong. Could anyone point me in the right direction?

I am new to heirarchical modeling so any help will be greatly appreciated!

data {
 int Trial; 
 int outcome[Trial,5];
 int choice[Trial,5];
 int stim1[Trial,5];
 int stim2[Trial,5];
 int feed[Trial,5]; 
 int  times[Trial,5]; 
 int<lower=0, upper=1> current_choice[Trial,5] ;  
 
 
}

transformed data {
 
 matrix[Trial,5] time_passed;
 for(j in 1:5){
   for(i in 2:Trial){
     time_passed[i,j]=(times[i,j]-times[i-1,j])/(24*3600*1000);
   }}
}

parameters {
 // hyperparameters for learning rate  
 real<lower=0,upper=1>  w; 
 real<lower=0> kap; 
 // and beta 
 real<lower=0> m_b; 
 real<lower=0> s_b;
 
 //per subject parameter 
 vector<lower=0,upper=1>[5] alpha_1;
 vector<lower=0>[5] beta_1;
 real  lambda_1; 
 
}

transformed parameters {
 real v_b= square(s_b);
 real shape_b=square(m_b)/v_b; 
 real rate_b= m_b/v_b; 
}

model {
 vector[2] current;
 real a; 
 real b; 
 real k; 
 // hyper prior for the learning rate - i.e. the overall learning rate 
 w ~ beta(1,9); // highest level priors 
 kap ~ gamma(.01,.01); //highest level priors - alpha didnt have a problem 
 k= kap+2; 
 a=w*(k-2)+1; 
 b=(1-w)*(k-2)+1; 
 //uniform higher level priors for beta_1 variable  
 m_b ~ uniform(1,20); 
 s_b ~ uniform(1,20);  
 
 for(n in 1:5){
   vector [Trial] PE; 
   vector[105] ev;
   alpha_1[n] ~ beta(a,b);
   PE=rep_vector(0,Trial);
   ev=rep_vector(0,105); 
   beta_1[n] ~ gamma(shape_b,rate_b); 
   lambda_1[n] ~ gamma(3,1);
   
   for (t in 1:Trial){
     if(t>1){
       
       ev = (.5-ev)* exp(-lambda_1[n]*time_passed[t,n])+.5; 
       
     }
     current[1] = ev[stim1[t,n]+1];
     current[2] = ev[stim2[t,n]+1];
     
     current_choice[t,n] ~ bernoulli_logit(beta_1[n]*current[2] - beta_1[n]*current[1]);
     
     if(feed[t,n]==1){
       PE[t] = outcome[t,n] - ev[choice[t,n]+1];
       ev[choice[t,n]+1] +=  alpha_1[n] * PE[t]; 
     }
     
   }
 }
}
1 Like
2

Before we start debugging the whole model, I’d like to make two requests:

i) Can you write a simpler version of your model that still displays the pathological behaviour?
ii) What do Gamma random variables drawn from this prior specification look like?
In other words, if I write
\mu_b \sim \operatorname{Uniform}(1, 20),
v_b \sim \operatorname{Uniform}(1, 20),
Y \sim \operatorname{Gamma}(\alpha(\mu_b, v_b), \beta(\mu_b, v_b)),
what does Y look like?
Finally, it’s important to note that the same way reparametrisation can make geometry nicer, it can also make it tougher. If you’re particularly unlucky, the reparametrisation with the most interpretable parameters will also be the one with worse geometry.

2 Likes
3

One more thing, maybe @maxbiostat can help to clear out.

The gamma distribution is not a location scale model, what I mean is that alpha hyper-parameter doesn’t behave like a mean.

So maybe re-parametrize a gamma in terms of a mean and variance is not a good idea

1 Like
4

Hi all, would try reparameterising a gamma in terms of its mean and shape (not variance) as the variance scales with the mean (given shape).

1 Like
5

Hi, have you made any progress with this issue? I’m curious! And good for posterity to report back ;)