comparison of hand-coded model to rstanarm

[Travis Riddle]

Hi all,

I'm trying to figure out why I'm getting a slightly different set of results from a simple model coded in rstanarm vs. one that I wrote myself. I have data on GPA for about 1200 students. The median number of measurements per student is 8, with a maximum of 12, and a handful for whom we have fewer than 8 because they left the district or dropped out or something. I've estimated a simple hierarchical, random intercept model like so:

data {

  int<lower=1> n; //number of obs

  int<lower=1> p; //number of participants

  vector[n] y; // observed vals

  int s[p]; //number of obs per participant

  int sub_obs[n]; //subject observation indicator

}

parameters {

  vector[p] b_mu; //mean for student p

  vector<lower=0>[p] sigma_i; //variance for student p

  real<lower=0> sigma; //one grand sigma 

  real mu;  //sampled grand mean

}

model {

  int pos;

  mu ~ normal(0,1);

  sigma_i ~ normal(0,1);

  b_mu ~ normal(mu, sigma);

  pos <- 1;

  for (i in 1:p){

    y[pos:pos+s[i]] ~ normal(b_mu[i], sigma_i[i]);

    pos <- pos+s[i];

  }

}

This appears to work really well - all three chains run in under a minute, there are no troublesome divergence warnings like a few other models I've tried, the chains look good, effective sample sizes, Rhats, etc don't raise any red flags. Further, the estimates are in line with what I'd expect - the median of the posterior for `mu` is 2.76, [2.71-2.80].

When I run what I assume is the same model in rstanarm:

stan_glmer(overallgpa ~ 1 + (1|ID), data=df.mod, prior_intercept=normal(0,1))

things are a bit different. It runs quite a bit more slowly (around 6 minutes per chain), and the estimates across multiple runs are consistently lower (median of `mu` is 2.69, [2.64-2.75]), and there are some troublesome diagnostics - high autocorrelation in the chains, low n_eff & Rhat > 1.01 for mu, and so on. My best guess is that rstanarm uses non-centered parameterizations by default, and this particular model happens to fall on the left-hand side of figures 8 and 9 in Betancourt & Girolami. I tried looking at the rstanarm code produced from the model, but it's a bit too hairy to make much sense of without some serious effort on my part.

Any guidance here would be much appreciated!

-Travis

[Ben Goodrich]

On Tuesday, May 31, 2016 at 3:33:05 PM UTC-4, Travis Riddle wrote:

When I run what I assume is the same model in rstanarm:

stan_glmer(overallgpa ~ 1 + (1|ID), data=df.mod, prior_intercept=normal(0,1))

things are a bit different. It runs quite a bit more slowly (around 6 minutes per chain), and the estimates across multiple runs are consistently lower (median of `mu` is 2.69, [2.64-2.75]), and there are some troublesome diagnostics - high autocorrelation in the chains, low n_eff & Rhat > 1.01 for mu, and so on. My best guess is that rstanarm uses non-centered parameterizations by default, and this particular model happens to fall on the left-hand side of figures 8 and 9 in Betancourt & Girolami. I tried looking at the rstanarm code produced from the model, but it's a bit too hairy to make much sense of without some serious effort on my part.

stan_glmer does use a non-centered parameterization. There is also a gamma prior on the between group standard deviations.

Ben

[Travis Riddle]

Thanks for the sanity check, Ben. Is the gamma prior responsible for the slightly lower values of mu?

[Ben Goodrich]

On Tuesday, May 31, 2016 at 4:10:28 PM UTC-4, Travis Riddle wrote:

Thanks for the sanity check, Ben. Is the gamma prior responsible for the slightly lower values of mu?

Maybe, but it is all interrelated. The gamma prior does have a longer tail than the half-normal you used.

Ben