Help converting hierarchical model from JAGS

[Noah]

Hello,

I have a hierarchical model in JAGS.  Appears to work fine, but takes forever on my current machine (latest core, i7, 16GB RAM, etc.)  So, I want to try things in stan.

Two questions:
1) Will it run faster/better from the command line than from within R with the rstan package?

2) Can someone help me with the format of the model file.  I'm having trouble understanding how to convert this jags code to stan code.  This is a hierarchical model with factors for both "group" and "program".  Looking at counts of students based on their group membership. (Pre-defined cluster) and the program they might enroll

Thanks!

Model below
================================

model{

    for(i in 1:N){
        eta[i] <-  ranef_group[group[i]] + ranef_prog[program[i]]+ b_var_1 * var_1[i] + b_var_2 * var_2[i]  + b_var_3 * var_3[i]
        log(lambda[i]) <- eta[i]
        count[i] ~ dpois(lambda[i])
    }
       
    b_var_1    ~ dnorm(0, tau_var_1)
    b_var_2 ~ dnorm(0, tau_var_2)
    b_var_3 ~ dnorm(0, tau_var_3)
   
    tau_var_1 <- pow(sigma_var_1, -2)
    tau_var_2 <- pow(sigma_var_2, -2)
    tau_var_3 <- pow(sigma_var_3, -2)
    tau_group <- pow(sigma_group, -2)
    tau_prog <- pow(sigma_prog, -2)
   
    sigma_var_1 ~ dunif(0,100)
    sigma_var_2 ~ dunif(0,100)
    sigma_var_3 ~ dunif(0,100)   
    sigma_group ~ dunif(0,100)
    sigma_prog ~ dunif(0,100)
   
   
    for (i in 1:n_group){ 
        ranef_group[i] ~ dnorm(0, tau_group)
    }
   
    for (i in 1:n_prog){ 
        ranef_prog[i] ~ dnorm(0, tau_prog)
    }       
}

[Bob Carpenter]

On 11/30/12 12:20 PM, Noah wrote:
> Hello,
>
> I have a hierarchical model in JAGS. Appears to work fine, but takes forever on my current machine (latest core, i7,
> 16GB RAM, etc.) So, I want to try things in stan.
>
> Two questions:
> 1) Will it run faster/better from the command line than from within R with the rstan package?

It runs the same compiled code. There's a bit of
overhead from using R, but that's only in getting
data in and out. It can take more size to store the
data in R, too.

> 2) Can someone help me with the format of the model file. I'm having trouble understanding how to convert this jags
> code to stan code. This is a hierarchical model with factors for both "group" and "program". Looking at counts of
> students based on their group membership. (Pre-defined cluster) and the program they might enroll
>
> Thanks!
>
> Model below
> ================================

You need to declare data and parameters. Look in
the manual to see how to do that, and make sure
you get the constraints on the deviation parameters sigma.

> model{
>
> for(i in 1:N){
> eta[i] <- ranef_group[group[i]] + ranef_prog[program[i]]+ b_var_1 * var_1[i] + b_var_2 * var_2[i] + b_var_3 *
> var_3[i]
> log(lambda[i]) <- eta[i]

No reason to make this separate, you just do

lambda[i] <- exp(eta[i]);

> count[i] ~ dpois(lambda[i])

You can just do this in Stan:

count[i] ~ poisson(exp(lambda[i]));

> }
>
> b_var_1 ~ dnorm(0, tau_var_1)

b_var_1 ~ normal(0,sigma_var_1);

> b_var_2 ~ dnorm(0, tau_var_2)
> b_var_3 ~ dnorm(0, tau_var_3)
>
> tau_var_1 <- pow(sigma_var_1, -2)
> tau_var_2 <- pow(sigma_var_2, -2)
> tau_var_3 <- pow(sigma_var_3, -2)
> tau_group <- pow(sigma_group, -2)
> tau_prog <- pow(sigma_prog, -2)
>
> sigma_var_1 ~ dunif(0,100)
> sigma_var_2 ~ dunif(0,100)
> sigma_var_3 ~ dunif(0,100)
> sigma_group ~ dunif(0,100)
> sigma_prog ~ dunif(0,100)
>
>
> for (i in 1:n_group){
> ranef_group[i] ~ dnorm(0, tau_group)
> }
>
> for (i in 1:n_prog){
> ranef_prog[i] ~ dnorm(0, tau_prog)
> }
> }

Also, things get executed in order, so that if
you did want to transform a variable for using it with <-
then you need to do it BEFORE using it.

- Bob

[Ben Goodrich]

On Friday, November 30, 2012 12:20:00 PM UTC-5, Noah wrote:

I have a hierarchical model in JAGS.  Appears to work fine, but takes forever on my current machine (latest core, i7, 16GB RAM, etc.)  So, I want to try things in stan.

For general remarks about speed, including the Matt trick used below, see the FAQ thread:

https://groups.google.com/d/msg/stan-users/4gv3fNCqSNk/J6ZItL2ZJ-IJ

Two questions:
1) Will it run faster/better from the command line than from within R with the rstan package?

The speed for 1 chain is essentially the same, but RStan makes it easier to execute chains in parallel and combine them for analysis. See

help(sflist2stanfit)

in rstan for an example of the parallel syntax, assuming you are not using Windows.
 

2) Can someone help me with the format of the model file.  I'm having trouble understanding how to convert this jags code to stan code.  This is a hierarchical model with factors for both "group" and "program".  Looking at counts of students based on their group membership. (Pre-defined cluster) and the program they might enroll

data {
  int<lower=1> N; // observations
  int<lower=1> J; // groups
  int<lower=1> P; // programs
  int<lower=1> K; // covariates

  int<lower=1> group[J];
  int<lower=1> program[P];
  matrix[N,K]; var;
  int<lower=0> count[N];
}
transformed data {
  real zeros[K];
  for (k in 1:K) { // needed for normal priors
    zeros[k] <- 0.0;
  }
}
parameters {
  vector[K] b;                          // coefficients
  vector[J] ranef_group_z;              // errors in random effects for groups
  vector[P] ranef_prog_z;               // errors in random effects for programs
  vector<lower=0,upper=100>[K] sigma_b; // can use something besides 100 for the upper-bound
  real<lower=0,upper=100> sigma_group;  // or no upper-bound at all ...
  real<lower=0,upper=100> sigma_prog;   // but then you probably need a proper prior on sigma_group and sigma_prog
}
transformed parameters {
  vector[J] ranef_group;
  vector[P] ranef_prog;
  // we call this the Matt trick, see the FAQ thread
  ranef_group <- sigma_prog * ranef_group_z;
  ranef_prog <- sigma_prog * ranef_prog_z;
}
model {
  real eta[N];

  for(i in 1:N) {
    eta[i] <-  ranef_group[group[i]] + ranef_prog[program[i]] + 

               // this part will lead to parser errors; see next message for a correction
               b_var_1 * var_1[i] + b_var_2 * var_2[i] + b_var_3 * var_3[i];
    count[i] ~ poisson(exp(eta[i]));
  }
  // In Stan >= v1.0.4, you will be able to do this, instead of the line directly above:
  // count ~ poisson_log(eta);

  // priors
  b ~ normal(zeros, sigma_var);
  ranef_group_z ~ normal(0.0, 1.0);
  ranef_prog_z ~ normal(0.0, 1.0);
  // if you use uniform priors on standard deviations then you don't need to specify them explicitly
}
generated quantities {
  vector[K] tau_b;
  real tau_group;
  real tau_prog;
  for (k in 1:K) {
    tau_b[k] <- 1.0 / pow(sigma_b[k], 2.0);
  }
  tau_group <- 1.0 / pow(sigma_group, 2.0);
  tau_prog <- 1.0 / pow(sigma_prog, 2.0);
}

That is not exactly the same model because putting a uniform prior on a standard deviation is not equivalent to putting a uniform prior on a precision. But that is basically how it could be done in Stan. The manual has a lot of similar code snippets.

Ben

[Ben Goodrich]

I said this part wrong:

On Friday, November 30, 2012 2:42:47 PM UTC-5, Ben Goodrich wrote:

model {
  real eta[N];

  for(i in 1:N) {
    eta[i] <-  ranef_group[group[i]] + ranef_prog[program[i]] +

               b_var_1 * var_1[i] + b_var_2 * var_2[i] + b_var_3 * var_3[i];

 
The way I set it up, there is no b_var_k but a vector b of length K because doing matrix algebra is faster. So, I should have said

model {
  vector[N] eta;
  eta <- var * b;
  for (i in 1:N) {
    eta[i] <- eta[i] + ranef_group[group[i]] + ranef_prog[program[i]];

and so on.

Ben

[Noah]

Ben,

Great advice.  Thanks for the help!