Use log form in the parameters block and unlog in the transformed parameters block

Modeling
1

Hi all,

I am new to this forum. I was hoping someone could help me with my interpretation of the following paper of Baio and Blangiardo

Basically, the model has the following priors:

home ~ Normal(0, 10000)
att ~ Normal(mu_att, tau_att)
def ~ Normal(mu_def, tau_def)

and the corresponding hyperpriors:

mu_att ~ Normal(0, 10000)
tau_att ~ Gamma(0.1, 0.1)
mu_def ~ Normal(0, 10000)
tau_fef ~ Gamma(0.1, 0.1)

Visually it looks as follows

The Stan code:

data {
   int<lower=1> n_teams;
   int<lower=1> n_matches;
   int<lower=1, upper=n_teams> home_team[n_matches];
   int<lower=1, upper=n_teams> away_team[n_matches];
   int<lower=0> y_1[n_matches];
   int<lower=0> y_2[n_matches];
}

parameters {
   real home;
   real mu_att;
   real mu_def;
   real tau_att;
   real tau_def;

   vector[n_teams-1] att_free;
   vector[n_teams-1] def_free;
}

transformed parameters {
   vector[n_teams] att;
   vector[n_teams] def;

   // sum-to-zero constraint
   for (team in 1:(n_teams-1)) {
      att[team] = att_free[team];
      def[team] = def_free[team];
   }
 
   att[n_teams] = -sum(att_free);
   def[n_teams] = -sum(def_free); 
}

model {
   vector[n_matches] theta_1;
   vector[n_matches] theta_2;

   // priors
   home ~ normal(0, 10000);
   mu_att ~ normal(0, 10000);
   mu_def ~ normal(0, 10000);
   tau_att ~ gamma(0.1, 0.1);
   tau_def ~ gamma(0.1, 0.1);

   att_free ~ normal(mu_att, 1/tau_att);    
   def_free ~ normal(mu_def, 1/tau_def);    

   for (match in 1:n_matches) {
      theta_1[match] = att[home_team[match]] + def[away_team[match]] + home;
      theta_2[match] = att[away_team[match]] + def[home_team[match]];

      y_1 ~ poisson_log(theta_1[match]);
      y_2 ~ poisson_log(theta_2[match]);
   }
}

The above Stan code is heavily influence by the following topic on this forum.
One of the answers says:

“Except your original posterior is going to be messed up because you did not constrain tau_att and tau_def to be positive. It would probably be better to declare them in log form in the parameters block and then antilog them in the transformed parameters block.”

How can I achieve this in my Stan code? Please advice.
If you see other room of improvements, I am all ears!
Thanks in advance!!

2

In principle you can always specify the priors in any way you want and use an arbitrary transformation of the variables afterwards (in Stan that would be the transformed parameters or model block). Most times you have to do this anyway once your parameters go into the model, to the transformation of parameters from those with priors on them is just another function in the overall model.

Specifying parameters on log scale is an example when you want the parameter to be positive, or when you don’t have good prior information on the scale of the parameter (\theta may be 1, 10, or 1000, for instance). So you can just exponentiate the parameter with the same base (e.g. for ln take exp, for log_2 take 2^\theta).

You won’t actually write the log in your code, that’s implicit in exponentiating it afterwards, so it’d look something like this:

parameters {
   real log_tau_att; // the log part is only in the name as a reminder
   real log_tau_def;
}

transformed parameters {
   real tau_att = exp(log_tau_att);
   real tau_def = exp(log_tau_def);
}

model {
   vector[n_matches] theta_1;
   vector[n_matches] theta_2;

   // priors
   log_tau_att ~ normal(mu, sigma);
   log_tau_def ~ normal(mu, sigma); // here you will have to think what are reasonable priors for this parameterization
}

keep in mind that if you are using gamma priors on your tau_ parameters you are already constraining them to be positive since they have zero probability of being lower than zero. So if that is the best way of representing your knowledge about the parameters you can keep that. Alternatively, you can constrain them by simply specifying the parameter as real<lower=0> tau_att;, that way you can use a half-normal or half-cauchy prior on the parameter, for instance.

1 Like
3

Thanks for reaching out. It indeed makes sense that using a gamma prior on my tau_att and tau_def constraints them to be >0. However, when I run the code I get the following error message:

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: normal_lpdf: Scale parameter is -0.156121, but must be > 0! (in ‘baio-blangiardo.stan’ at line 43)
Exception: normal_lpdf: Scale parameter is -0.311934, but must be > 0! (in ‘baio-blangiardo.stan’ at line 44)

Line 43 and 44 are:

att_free ~ normal(mu_att, 1/tau_att);    
def_free ~ normal(mu_def, 1/tau_def);

So, whats going wrong? In addition, why do we need to do 1 / x?
Thanks for helping out

4

If the numbers were lower I’d think it’s some kind of numerical instability, where you get things like negative 1e-10 instead of zero and it will raise an error. Not sure if there’s an underlying issue where you are actually getting negative values where you shouldn’t, or if they haven’t been evaluated as zero probability and still go through as “valid” because the posterior probability hasn’t been computed yet.

If the latter is the case it will be rejected either way, it’s just that the normal will raise an exception for negative variance (invalid parameter) before actually computing the the likelihood but not the gamma (zero-probability observation). I think the poisson will complain about its parameter being zero, so it’s also surprising that the gamma won’t.

In any case, unless this is a symptom of something actually wrong it just means the proposal will be rejected, so it won’t affect your posterior and if it doesn’t happen often shouldn’t affect performance much.

1 Like
5

I think the problem here is that Stan doesn’t initialize the chains from the prior distribution, but rather on the interval [-2,2] on the unconstrained scale. To avoid bad inits, it would be good to specify <lower = 0> when you declare these parameters. Under the hood, this is equivalent to declaring the parameters of the log scale, exponentiating, doing the necessary Jacobian adjustment, and then sampling.

Some alternative options (each of which seems like an unnecessary hassle) would be:

  • Leaving the model untouched but passing your own, valid inits. This could still run into numerical pathologies near the boundary, however.
  • Declaring your constrained parameters on the unconstrained scale (the log scale), exponentiating, keeping your gamma prior for the exponentiated parameter, and doing the necessary Jacobian adjustment.
  • Declaring your constrained parameters on the log scale and exponentiating, but figuring out a good prior to use on the log-scale parameter (rather than the exponentiated parameter), which eliminates the need for a Jacobian adjustment.
1 Like