Rhat and ESS from the print(stan_fit) different from the functions

General
1

Hello, I estimated a model using Rstan and with function stan(), then checked the results with function print(stan_fit, pars="…"). It returned back very small ESS and large Rhat. However, when I used the function Rhat(), ess_bulk(), and ess_tail() to compute these indices again, the results seemed much more normal. I don’t understand why I got different values and which one I should believe. Could anyone help me with it? Thanks very much!

1 Like
2

It’s likely that you have different versions of these statistics called. Can you specify how did you call Rhat and ESS functions, so I can tell the specific versions. See Search 2.5 million pages of mathematics and statistics articles for description of different versions.

Can you tell the values that you consider small, large and normal?

3

Thanks for your reply! Here is how I called Rhat and ESS functions:
I first used stan_est = extract(stan_fit) to get the parameter estimates in all iterations. Then used theta_est = stan_est$theta to get estimates of a specific parameter, theta.
Because I had 2 chains and each with 1000 iterations, so I formatted theta_est as: theta_est2 = matrix(theta_est, ncol=2). I guess this may let me have iterations from 1 chain in each column, but I can be wrong. Lastly, I called Rhat and ESS functions as: Rhat(theta_est2); ess_bulk(theta_est2); ess_tail(theta_est2).

4

Ah this explains it. The bad design choice in extract strikes again. If look at the help of extract

     extract(object, pars, permuted = TRUE, inc_warmup = FALSE, 
       include = TRUE)

you see that by default it permutes the iterations which discards the needed information about chains and autocorrelations for Rhat and ESS. Use
stan_est = extract(stan_fit, permuted=FALSE)
or
stan_est = as.matrix(stan_fit)
or even better, use posterior package

We also recommend to run at least 4 chains.

6

Thank you very much!