Hard coding a dataset and parameter relationship in Stan which works in Jags

[Daniel Padfield]

Hi Folks

I am just starting out on my Bayesian adventures and am just trying to get my head round some code my supervisor has given me from another collaborator.

We have a bunch of community metabolic rate data that is dependent on temperature and mass. We have 1 temperature point per rate point but variable numbers of mass points (as there are > 1 individuals that make each community).

The relationship of metabolic rate with mass is a power law at the individual level and my collaborator hard-coded it into the jags model from the raw mass points. A variable for each rate variable is created by log(mass_1^a + .... mass_i^a). A simplified example for 3 points of rate and temperature would be:

'model{
  # priors
  bT ~ dnorm(0, 0.001)        # effect of temperature
  a ~ dnorm(0, 0.001)         # size scaling exponent
  ln_r0 ~ dnorm(0, 0.001)     # community normalisation constant
  tau_R ~ dgamma(1.0E-3,1.0E-3)
  var_R <- 1/tau_R


  # size corrected biomass
  ln.biomass[1] <- log(
  pow(2.54, a) + pow(3.28, a) + pow(3.54, a) + pow(2.64, a))
  ln.biomass[2] <- log(
  pow(2.64, a) + pow(2.61, a) + pow(2.67, a) + pow(3.36, a) + pow(2.86, a))
  ln.biomass[3] <- log(
  pow(2.65, a) + pow(3.60, a) + pow(3.39, a))


  for(i in 1:length(ln_rate)){
    ln_rate[i] ~ dnorm(mu[i], tau_R)
    mu[i] <- log(exp(ln_r0) * exp(ln.biomass[i])) + bT*temp[i]
  }
}
'

I was wondering how I could go about hard coding this relationship in Stan or whether there was another way of doing it in Stan using arrays or matrices.

So far I have tried to alter the jags code to Stan by putting the ln.biomass[] code into the model{} specification but I somehow doubt this is the best way to do such things.

'data{
  int N;                              // number of samples
  vector[N] ln_rate;                  // rate
  vector[N] temp;                     // temperature
}
parameters{  
  real<lower=0> sigma;      // standard deviation
  real ln_r0;               // community normalisation constant
  real bT;                  // effect of temperature
  real a;                   // size scaling exponent
}
model{
  vector[N] mu;
  vector[N] ln_biomass;
  ln_r0 ~ normal(0,10);
  bT ~ normal(0,10);
  a ~ normal(0,10);
  sigma ~ cauchy(0,2);
 
  ln_biomass[1] = log(
    pow(2.54, a) + pow(3.28, a) + pow(3.54, a) + pow(2.64, a));
  ln_biomass[2] = log(
    pow(2.64, a) + pow(2.61, a) + pow(2.67, a) + pow(3.36, a) + pow(2.86, a));
  ln_biomass[3] = log(
    pow(2.65, a) + pow(3.60, a) + pow(3.39, a));
 
  # likelihood regression for size dependence of rate
  for(i in 1:N) {
    mu[i]  = log(exp(ln_r0) * exp(ln_biomass[i])) + bT*temp;
  }
  ln_rate ~  normal(mu, sigma);
}'

Any advice on how best to code this would be much appreciated.

Many thanks

Dan

[Bob Carpenter]

Sorry for not getting back to you sooner---I've been swamped
with grant applications and StanCon.

The short answer is that you can translate that model line for line,
but it includes a lot of bad modeling practices inherited from BUGS
(too diffuse priors) and things we just don't need to do in Stan (work
with precision matrices).

See below for some more details.

> On Jan 18, 2017, at 5:09 AM, Daniel Padfield <d.padfi...@gmail.com> wrote:
>
> Hi Folks
>
> I am just starting out on my Bayesian adventures and am just trying to get my head round some code my supervisor has given me from another collaborator.
>
> We have a bunch of community metabolic rate data that is dependent on temperature and mass. We have 1 temperature point per rate point but variable numbers of mass points (as there are > 1 individuals that make each community).
>
> The relationship of metabolic rate with mass is a power law at the individual level and my collaborator hard-coded it into the jags model from the raw mass points. A variable for each rate variable is created by log(mass_1^a + .... mass_i^a). A simplified example for 3 points of rate and temperature would be:
>
> 'model{
> # priors
> bT ~ dnorm(0, 0.001) # effect of temperature
> a ~ dnorm(0, 0.001) # size scaling exponent
> ln_r0 ~ dnorm(0, 0.001) # community normalisation constant
> tau_R ~ dgamma(1.0E-3,1.0E-3)
> var_R <- 1/tau_R
>
>
> # size corrected biomass
> ln.biomass[1] <- log(
> pow(2.54, a) + pow(3.28, a) + pow(3.54, a) + pow(2.64, a))
> ln.biomass[2] <- log(
> pow(2.64, a) + pow(2.61, a) + pow(2.67, a) + pow(3.36, a) + pow(2.86, a))
> ln.biomass[3] <- log(
> pow(2.65, a) + pow(3.60, a) + pow(3.39, a))
>
>
> for(i in 1:length(ln_rate)){
> ln_rate[i] ~ dnorm(mu[i], tau_R)
> mu[i] <- log(exp(ln_r0) * exp(ln.biomass[i])) + bT*temp[i]
> }
> }
> '
>
> I was wondering how I could go about hard coding this relationship in Stan or whether there was another way of doing it in Stan using arrays or matrices.
>
> So far I have tried to alter the jags code to Stan by putting the ln.biomass[] code into the model{} specification but I somehow doubt this is the best way to do such things.

That's actually fine. Given that it depends on a parameter (a), you can't
precompute it. Whether you make it a transformed parameter or a local
variable depends on whether you want to save it (transformed parameters
get saved).

>
> 'data{
> int N; // number of samples
> vector[N] ln_rate; // rate
> vector[N] temp; // temperature
> }
> parameters{
> real<lower=0> sigma; // standard deviation
> real ln_r0; // community normalisation constant
> real bT; // effect of temperature
> real a; // size scaling exponent
> }
> model{
> vector[N] mu;
> vector[N] ln_biomass;
> ln_r0 ~ normal(0,10);
> bT ~ normal(0,10);
> a ~ normal(0,10);
> sigma ~ cauchy(0,2);

These aren't the same priors, but should be better if you know the scales
of the variables are roughly in (-20, 20) and if you think sigma might
be a lot higher than 2 (otherwise, you can use a half normal instead of
half Cauchy).

> ln_biomass[1] = log(
> pow(2.54, a) + pow(3.28, a) + pow(3.54, a) + pow(2.64, a));
> ln_biomass[2] = log(
> pow(2.64, a) + pow(2.61, a) + pow(2.67, a) + pow(3.36, a) + pow(2.86, a));
> ln_biomass[3] = log(
> pow(2.65, a) + pow(3.60, a) + pow(3.39, a));

This would be more arithmetically stable as log_sum_exp calculation:

ln_biomass[1] = log_sum_exp({a * log(2.64), a * log(2.67), ...});

And then you can put it all together as a one-liner if you'd prefer:

vector[N] ln_biomass
= to_vector({ log_sum_exp({a * log(2.64), a * log(2.67), ...}),
log_sum_exp({...}),
log_sum_exp({...}) });

but it's probably clearer the way you have it. Mitzi's working on
vector expressions but we won't have that until Stan 2.15 at the earliest.

> # likelihood regression for size dependence of rate
> for(i in 1:N) {
> mu[i] = log(exp(ln_r0) * exp(ln_biomass[i])) + bT*temp;
> }

log and exp round trips are a red flag for numerical stability. Instead,
use this:

mu[i] = ln_r0 + ln_biomass[i] + bT * temp;

which then lets you vectorize the whole loop as:

mu = ln_r0 + bT * temp + ln_biomass;

I reordered so the scalar additions happen first. At that point, this is
simple enough that you don't need a local, and could just drop the expression
for mu into the normal itself.

> ln_rate ~ normal(mu, sigma);
> }'
>
> Any advice on how best to code this would be much appreciated.

If you don't mind, please let us know how it works compared to the JAGS
version (you need to use Stan's n_eff calculations for a fair comparison
as our n_eff estimates are much more conservative than JAGS to account
for potential non-convergence---they converge to the same answer as
R-hat goes to 1).

- Bob

[Daniel Padfield]

Hi Bob

Thanks so much for this worked a charm on the example. Will at some point get round to applying it to my supervisors old data to compare the two methods for you.

Until next time

Cheers

Dan