ggplot2 error when using stan_diag

[Diego Barneche]

Hello everyone,

I am interested in learning more about posteriors visualisation using stan model fits, but I'm quite new to stan. For example, I do not really understand why the function stan_diag is failing on my machine using a very simple linear model fit (traceplot runs fine though). I have just installed the latest R and rstan versions (and all associated packages). Looking at the error message, it seems that the error could be related to the function grid.arrange from package gridExtra (potential explanations here and here)? Or am I missing something very simple? 

Thanks very much for your assistance!

Best, Diego

Please see below my .stan and .R codes, sessionInfo(), and the error message (in red) generated when using stan_diag:

test.stan

data {

int<lower=1> N; //number of observations

real y[N]; //response variable

real x[N]; //predictor variable

}

parameters {

real b1; //fixed effect - intercept

real b2; //fixed effect - slope

real<lower=0> sdModel;

}

model {

b1       ~  normal(0, 1/1.0E-6);

b2       ~  normal(0, 1/1.0E-6);

sdModel  ~  cauchy(0, 5);

    for(i in 1:N) {

    y[i] ~  normal(b1 + b2*x[i], sdModel);

    }

}

generated quantities {

    real logLik[N];

    for(i in 1:N) {

    logLik[i]   <-  normal_log(y[i] , b1 + b2*x[i], sdModel);

    }

}

R code:

library(rstan)

rstan_options(auto_write = TRUE)

options(mc.cores = parallel::detectCores())

sessionInfo()

x  <-  seq(-10, 10, length.out=100)

y  <-  2 - 3 * x + rnorm(100)

stanDat   <-  list(N  =  length(x),

              x  =  x,

              y  =  y)

stanFit  <-  stan(file = '~/Desktop/test.stan', data=stanDat, iter=1e4, thin=5, chains=3)

print(stanFit, pars=c('b1', 'b2', 'sdModel'), probs=c(0.025, 0.5, 0.975)) # works fine

traceplot(stanFit, pars=c('b1', 'b2', 'sdModel')) # works fine

stan_diag(stanFit, pars=c('b1', 'b2', 'sdModel')) # fails

sessionInfo():

> sessionInfo()

R version 3.3.1 (2016-06-21)

Platform: x86_64-apple-darwin13.4.0 (64-bit)

Running under: OS X 10.11.4 (El Capitan)

locale:

[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8

attached base packages:

[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:

[1] rstan_2.10.1         StanHeaders_2.10.0-2 ggplot2_2.1.0       

loaded via a namespace (and not attached):

 [1] colorspace_1.2-6 scales_0.4.0     plyr_1.8.4       parallel_3.3.1   inline_0.3.14    gtable_0.2.0     gridExtra_2.2.1  Rcpp_0.12.5      grid_3.3.1      

[10] stats4_3.3.1     munsell_0.4.3   

Error (related to ggplot2?):

> print(stanFit, pars=c('b1', 'b2', 'sdModel'), probs=c(0.025, 0.5, 0.975)) # works fine

Inference for Stan model: test.

3 chains, each with iter=10000; warmup=5000; thin=5; 

post-warmup draws per chain=1000, total post-warmup draws=3000.

         mean se_mean   sd  2.5%   50% 97.5% n_eff Rhat

b1       2.15       0 0.11  1.93  2.15  2.37  2798    1

b2      -3.01       0 0.02 -3.04 -3.01 -2.97  2777    1

sdModel  1.06       0 0.08  0.92  1.06  1.23  2997    1

Samples were drawn using NUTS(diag_e) at Thu Jun 30 08:58:42 2016.

For each parameter, n_eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor on split chains (at 

convergence, Rhat=1).

 The estimated Bayesian Fraction of Missing Information is a measure of

 the efficiency of the sampler with values close to 1 being ideal.

 For each chain, these estimates are

 0.8 0.9 1

> traceplot(stanFit, pars=c('b1', 'b2', 'sdModel')) # works fine

> stan_diag(stanFit, pars=c('b1', 'b2', 'sdModel')) # fails

Error in gList(lp_hist = list(grobs = list(list(x = 0.5, y = 0.5, width = 1,  : 

  only 'grobs' allowed in "gList"

In addition: Warning messages:

1: In grob$wrapvp <- vp : Coercing LHS to a list

2: In grob$wrapvp <- vp : Coercing LHS to a list

3: In grob$wrapvp <- vp : Coercing LHS to a list

> 

[Jonah Gabry]

Do you have the most recent versions of the grid and gridExtra packages installed? It's possible that we're depending on a minimum version of those without realizing it.

While we sort this out, you can get similar plots to stan_diag with our shinystan package/GUI.

Jonah

[Diego Barneche]

Thanks for your reply Jonah.

Yes, I started my R installation from scratch yesterday including all packages installations. They should be the newest (as per version on CRAN).

I'll try shinystan now, thanks for the suggestion.

[Jonah Gabry]

Hi Diego, there was something I didn't notice in the code you posted: there's no "pars" argument to the stan_diag function, which is why it's not working for you. But it's not giving you an informative error message either, which is a problem. I'll add an informative error message in the next release. 

Jonah

[Jonah Gabry]

I just fixed this on the develop branch of rstan on github: 

https://github.com/stan-dev/rstan/issues/308

[Diego Barneche]

Oh, silly mistake! Sorry about that, and thanks very much for your answer. It works now!

Also, thanks for fixing the error message, that's very helpful.

Best,

Diego