Speed comparison stan vs MCMCglmm for hierarchical regression

[Ryan King]

I mentioned on Andrew's blog that I had not had good results with stan v1.0 in high dimension, and he suggested I try again and post here.

The setup is a hierarchical regression model with two groups of random effects, one with 100 members and one with 150 members with 900 total individuals. The numbers from from a real dataset that I don't have permission to share, but I simulate data with similar characteristics. I try out both a continuous and binary outcome. I haven't done any real tuning, but have done the obvious things recommended in the manual.  The design matrices are somewhat sparse, and I know that MCMCglmm uses sparse matrix techniques. I didn't see an option for those in stan, or if eigen does so automatically.

The competitor is MCMCglmm, which uses a specialized block-gibbs sampler for regression problems. I dedicate 10% of the iterations to warmup and burnin for each. The results are that (r)stan 1.3 is quite a bit better than I recall from 1.0. The engines are similar on this small test, but I haven't dedicated the CPU time to do numerous replications. 

On continuous outcomes (a LMM) MCMCglmm is a bit faster (6.4 effective samples / sec vs 4.1 ES/s), and that the variation of stan is higher (on three trials the minimum ES/s was 1.8; watching the counter go the iteration number has a propensity to go quickly in streaks).

On binary outcomes (a GLMM), stan does better than MCMCglmm with the slice sampling option (2.6 ES/s vs 1.4 ES/s), which is superior to not slice sampling. Again, fairly variable ES/s.

Attached are scripts for both tests and a session info.

Is this a problem likely to benefit from tuning in the warmup or mass matrix options?

Thanks,

Ryan King

University of Chicago Pritzker School of Medicine

[Andrew Gelman]

Ryan:

Thanks for posting.  Just to be clear:  are you using the same effective sample size formula for MCMCglmm and Stan?

A

--
You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
 
 
<sessionInfo.txt><new.stan.mixed.modeltest.binary.R><new.stan.mixed.modeltest.R>

[Bob Carpenter]

Thanks so much for following up. The main reason for
the speedup from Stan 1.0 to 1.3 is that we've vectorized
the distributions.

Your models look about as efficient as they can
be coded in Stan. (There is an extra N-vector declaration
in the models for mu, but that shouldn't slow things
down noticeably.)

We've also noticed the high variation among chains for timing.
Part of it is initialization --- we use very diffuse initializations.
You might want to try it with init=0, which centers the initializations.
Even though they're centered to start (on the unconstrained scale),
they'll immediately get a random kick for the initial momentum.

We are working on adaptation with the goal of bringing down
the variation among chain timings.

- Bob

[Ryan King]

No, thanks for reminding me! I was using stan's n_eff printout and coda:effectiveSize. Using coda:effectiveSize for stan is slightly more generous for the binary outcome (3.1 ES/s).

Oddly enough, it made the continuous outcome look worse (2.7 ES/s). I ran that sampler again to make sure it wasn't a fluke. 

It is impressive that stan does about the same on gaussian outcomes and binary ones; binary outcomes are generally a much harder computation problem. Do you think it's a shortcoming of coda:effectiveSize? It's not clear how to access the function rstan is using for n_eff. I also noticed that longer warmup (1000) and run times (10k) tended to produce higher ESS than short ones (500 and 5000).

This is the correct way to access the ESS?
(effectiveSize(fit1@sim$samples[[1]]$effect2sd) +  effectiveSize(fit1@sim$samples[[1]]$effect1sd)) /2

Thanks,

Ryan King

[Jiqiang Guo]

I am not sure that I understood your question about the function for ESS in rstan.  rstan does not expose a function that calculates only the ESS.  But summary function for the fitted object will compute the ESS.  For the first example in https://code.google.com/p/stan/wiki/RStanGettingStarted#How_to_Install_RStan

for example, 

> summary(fit)$summary[,"n_eff"]

       mu       tau    eta[1]    eta[2]    eta[3]    eta[4]    eta[5]    eta[6]    eta[7]    eta[8]  theta[1]  theta[2]  theta[3] 

 160.3627  109.7633  949.8854 1139.1143  385.3387  802.7463 1044.9465 1104.5331  815.4161 1056.2776  397.1705  629.7662  232.4378 

 theta[4]  theta[5]  theta[6]  theta[7]  theta[8]      lp__ 

 416.2000  436.4769  345.8459  523.2209  450.7623  195.4303 

--

Jiqiang 

On Sun, Jun 2, 2013 at 12:50 PM, Ryan King <c.rya...@gmail.com> wrote:

It's not clear how to access the function rstan is using for n_eff.

[Andrew Gelman]

Or you can use the monitor() function to get ESS and Rhat from a 3-way array of simulations that are produced some other way, for example from Jags or MCMCglmm.

A

[Bob Carpenter]

To do apples-to-apples comparisons, you could use Coda.

It's possible using Rstan's read_stan_csv function or
Stan command-line bin/print function, but I don't see doc
for either of their input formats. The R function
doc says

This function reads the samples and using the
comments (beginning with "#") to create a stanfit object.

and the Stan doc says nothing about the input format (chapter
5 is about bin/print).

- Bob

>> stan-users+...@googlegroups.com <mailto:stan-users+...@googlegroups.com>.

[Andrew Gelman]

No, please don't use Coda. We want to use the ESS formula that we like best. To compare using coda could lead to problems if the chains are not mixing well.

[Bob Carpenter]

On 6/3/13 4:52 PM, Andrew Gelman wrote:
> No, please don't use Coda. We want to use the ESS formula that we like best. To compare using coda could lead to problems if the chains are not mixing well.

For someone to use RStan or bin/print, we'd need to document
what the input format looks like. I just took that up on
the stan-dev list (the question of should we do it being
now answered in the affirmative).

Until then, the choice seems to be to use Coda (with its more
generous R-hat and n_eff calcs) or to code up your own definitions.

A simple example or two might suffice until we get around to the
doc.

- Bob

[Andrew Gelman]

It's easy to use monitor() which is in rstan. It goes like this:
1. If the number of iterations per chain is "iter", and the number of chains is "chains", and the number of parameters saved your mcmc output is "d", store the simulations as a 3-way array of dimensions iter, chains, d (in that order).
2. Suppose those chains have been saved into a 3-way array called sims, and let warmup be the number of warmup iterations. Then just do monitor (sims, warmup).

It's all obtainable by doing the following in R:
library ("rstan")
help (monitor)

I really really don't want people to be using coda for comparisons!
A

[Bob Carpenter]

On 6/3/13 5:10 PM, Andrew Gelman wrote:
> It's easy to use monitor() which is in rstan. It goes like this:
> 1. If the number of iterations per chain is "iter", and the number of chains is "chains", and the number of parameters saved your mcmc output is "d", store the simulations as a 3-way array of dimensions iter, chains, d (in that order).
> 2. Suppose those chains have been saved into a 3-way array called sims, and let warmup be the number of warmup iterations. Then just do monitor (sims, warmup).
>
> It's all obtainable by doing the following in R:
> library ("rstan")
> help (monitor)
>
> I really really don't want people to be using coda for comparisons!

I didn't realize that there was a function that could do this
in RStan, but indeed monitor() is just what's needed to compute
summaries of samples and it's very clearly documented. I just
hadn't seen it.

I don't think there's really a need to do this kind of
thing in command-line Stan.

- Bob

[Ryan King]

Using rstan:monitor jitters the answers a bit

Continuous:

MCMCglmm 8.3 n_e/s

stan 2.7 n_e/s

Binary:

MCMCglmm (expanded) 1.4 n_e/s

MCMCglmm (slice) 0.7 n_e/s

stan 0.6 n_e/s

Thanks,

Ryan King