Precision of Rhat values for posterior means - stan_summary

[KelsonGs]

Hey there,

In the following discussion "Rhat" refers to the Rhat value for posterior means as presented in the summary output file. The issue I raise relates to how Rhat values are reported in the StanSample.jl summary output file (given by read_summary method). The models I am using have Rhat values in the interval (1.0,1.1), these values appear to be rounded to 1dp, hence for any Rhat < 1.05 the corresponding summary file presents the value as 1.0. The motivation for this choice of rounding Rhat to 1.0 appears to follow from the discussion on Rhat given in the CmdStan user guide (linked below) which suggests Rhat values less than 1.05 are desirable, therefore rounding to 1dp is sufficient for assessing Rhat under this criterion.

The issue I have is that I am going by the Rhat criterion presented in Vehtari et al. (2020), this criterion suggests Rhat values be < 1.01. Therefore, to assess whether my Rhat values fulfil this criterion I need more precise Rhat values, i.e. Rhat values rounded to >= 2dp. Alternatively, is there some existing method or workaround that would allow me to generate Rhat values with higher precision?

Thanks!

Refs:
Vehtari et al. (2020), https://arxiv.org/abs/1903.08008 - Accessed: http://www.stat.columbia.edu/~gelman/research/published/rhat.pdf

CmdStan user guide; stansummary discussion. Accessed: https://mc-stan.org/docs/cmdstan-guide/stansummary.html#model-parameters-and-quantities-of-interest

I am using CmdStan version 2.32.2 and Julia version 1.9.1.

[goedman]

Hi Kelson ( @KelsonGs ),

Thanks for your email. I'll take a look asap (I'm traveling at the moment and don't have very good connectivity). Or you can make a PR.

Best
Rob

[goedman]

Just to confirm we're talking about the same method, if I run:

In [2]: include("/Users/rob/.julia/dev/StanSample/example/bernoulli.jl");
[ Info: /Users/rob/.julia/dev/StanSample/example/tmp/bernoulli.stan updated.
4000×1 DataFrame
  Row │ theta     
      │ Float64   
──────┼───────────
    1 │ 0.573698
    2 │ 0.55348
    3 │ 0.429174
    4 │ 0.186664
    5 │ 0.125927
    6 │ 0.109919
    7 │ 0.265181
    8 │ 0.274059
    9 │ 0.241493
<snipped>

I get:

In [3]: read_summary(sm)
8×10 DataFrame
 Row │ parameters     mean       mcse          std        5%         50%        95%        ess         n_eff/s      r_hat        
     │ Symbol         Float64    Float64       Float64    Float64    Float64    Float64    Float64     Float64      Float64      
─────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ lp__           -8.16742     0.0187357   0.760978   -9.66123   -7.87367   -7.64004   1649.7       48520.6       1.00034
   2 │ accept_stat__   0.915255    0.00193465  0.127225    0.639421   0.970274   1.0       4324.56     127193.0       1.00107
   3 │ stepsize__      1.01558     0.0657878   0.0931778   0.85911    1.05703    1.10394      2.00602      59.0005    2.22175e13
   4 │ treedepth__     1.3655      0.0080534   0.48163     1.0        1.0        2.0       3576.59     105194.0       0.999351
   5 │ n_leapfrog__    2.3995      0.0191298   1.07642     1.0        3.0        3.0       3166.25      93125.0       1.0049
   6 │ divergent__     0.0       NaN           0.0         0.0        0.0        0.0        NaN           NaN       NaN
   7 │ energy__        8.66522     0.0263062   1.04349     7.69044    8.34701   10.7545    1573.47      46278.6       1.00074
   8 │ theta           0.333709    0.00355884  0.130675    0.131946   0.326891   0.561284  1348.25      39654.3       1.00231

[KelsonGs]

Hi Rob, Great to hear from you. The output you provided is correct in that it is the output from the stan_summary method. Curiously, your output gives theta with quite a few sig-figs, my output (picture below) has all the r_hat values for the parameters of interest given to 1dp. I've found that stan_summary always gives me param r_hat values to 1dp. Some of my models (with lower draws & warmups) give r_hat of 1.1, but the value never gets more precise than 1dp. Given that your output does not have the same problem, might this be a version issue? I am running CmdStan 2.3.2 using Miniconda, and the StanSample v7.4.2. Best, Kelson [image: image.png]

…

On Thu, 31 Aug 2023 at 17:11, Rob J Goedman ***@***.***> wrote: Just to confirm we're talking about the same method, if I run: include("/Users/rob/.julia/dev/StanSample/example/bernoulli.jl"); I get: In [3]: read_summary(sm) 8×10 DataFrame Row │ parameters mean mcse std 5% 50% 95% ess n_eff/s r_hat │ Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 ─────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── 1 │ lp__ -8.16742 0.0187357 0.760978 -9.66123 -7.87367 -7.64004 1649.7 48520.6 1.00034 2 │ accept_stat__ 0.915255 0.00193465 0.127225 0.639421 0.970274 1.0 4324.56 127193.0 1.00107 3 │ stepsize__ 1.01558 0.0657878 0.0931778 0.85911 1.05703 1.10394 2.00602 59.0005 2.22175e13 4 │ treedepth__ 1.3655 0.0080534 0.48163 1.0 1.0 2.0 3576.59 105194.0 0.999351 5 │ n_leapfrog__ 2.3995 0.0191298 1.07642 1.0 3.0 3.0 3166.25 93125.0 1.0049 6 │ divergent__ 0.0 NaN 0.0 0.0 0.0 0.0 NaN NaN NaN 7 │ energy__ 8.66522 0.0263062 1.04349 7.69044 8.34701 10.7545 1573.47 46278.6 1.00074 8 │ theta 0.333709 0.00355884 0.130675 0.131946 0.326891 0.561284 1348.25 39654.3 1.00231 — Reply to this email directly, view it on GitHub <#73 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AJG2IS2UQY2WQFNV7NLGPXDXYAMIJANCNFSM6AAAAAA4EF2RLE> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[KelsonGs]

1537×10 DataFrame
Row │ parameters mean mcse std 5% 50% 95% ess n_eff/s r_hat
│ String Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64
──────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
1 │ lp__ -2431.6 3.75133 129.663 -2621.72 -2442.95 -2204.98 1194.71 0.0885355 1.00222
2 │ accept_stat__ 0.841554 0.00304239 0.132944 0.58228 0.873387 0.995751 1909.45 0.141502 1.00385
3 │ stepsize__ 0.173684 0.00423944 0.00599574 0.167575 0.178405 0.180828 2.00017 0.000148225 2.08698e11
4 │ treedepth__ 4.88718 0.00855484 0.316376 4.0 5.0 5.0 1367.67 0.101353 1.01297
5 │ n_leapfrog__ 29.44 0.124608 4.75315 15.0 31.0 31.0 1455.03 0.107827 1.01142
6 │ divergent__ 0.0 NaN 0.0 0.0 0.0 0.0 NaN NaN NaN
7 │ energy__ 3191.58 3.79672 132.583 2961.03 3202.16 3387.69 1219.43 0.0903676 1.00213
8 │ beta[1.1] 0.26 0.001 0.44 -0.46 0.26 1.0 197282.0 15.0 1.0
9 │ beta[1.2] 2.8 0.011 0.88 1.4 2.7 4.3 5993.0 0.44 1.0
10 │ beta[1.3] -0.0042 0.00099 0.46 -0.77 -0.0035 0.76 218355.0 16.0 1.0
11 │ beta[1.4] -0.49 0.0021 0.47 -1.3 -0.48 0.27 50786.0 3.8 1.0
12 │ beta[1.5] -0.97 0.0043 0.51 -1.8 -0.95 -0.16 14243.0 1.1 1.0

[goedman]

@KelsonGs

Any chance you can send me your Stan model and (simulated?) data so I can reconstruct the issue?

I'm working and testing on Julia 1.9, 1.10 and 1.11_DEV and currently using cmdstan 2.32.0. You can see that in the tmp directory .csv files, e.g. in my example bernoulli_chain_1.csv. The summary file is in bernoulli_summary.csv. Can you send me your summary file?

Just copy and paste them between 2 "```" lines, e.g.

Here your copied stuff

And for StanSample I use v7.4.2

This looks more like a modeling issue, but I'm not sure.

Rob

[KelsonGs]

@goedman

I've emailed you the relevant files.

Best,
Kelson

[goedman]

@KelsonGs

Thanks, got them earlier while biking. Will try to generate simulated input and see how far I can get.

If you have the pre-computed data, that would be even better.

Thanks,
Rob

[goedman]

Hi Kelson ( @KelsonGs ),

Had a quick look at the summary result and HB_Stan.jl.

In HB_Stan.jl you use:

...
  df_summ = read_summary(sm)
  df_summ = DataFrame(df_summ)
...

The second DataFrame conversion is not doing anything I think, but also would not lead to rounding.

I am surprised by the estimation control settings for MCMC, these seem very high, and combined with the fact that the summary output seems to indicate that most estimated values are smaller than the standard deviation I wonder if there is a (much) smaller generative model for your data, something like (if that makes sense!):

X = data[1]
Y = data[2]
  
data= Dict("n_attribute"=> 3,
  "n_subject" => 5
  "n_question" => 10,
  "y_sim" => Y,
  "x_sim" => X
  )

to check the model.

But you might already have done all of that. In that case, I would probably try to post the model with needed data on the Stan mailing list. Or as an interim step I could try to run it. But my real experience with these kind of models is limited.

Best,
Rob

[KelsonGs]

Hi Rob,

Thanks for taking a look. I have emailed you one of the data files used in the estimation and the relevant pre-processing code to read, format and estimate (using HB_Stan.jl) needed to replicate the issue.

If you were interested in running the model, you may want to significantly decrease the samples and warmups in the MCMC controls for the HB_Stan.jl file I sent -- seeing as the actual results do not matter in regards to uncovering the issue of rhat precision in the summary output. I have done such tests myself using 5000 draws, 1000 warmups and found the same issue with the precision of rhat values in the summary output.

I will try running the Bernoulli example code you referenced in your first reply, the output from this will confirm whether the issue is with my specific model/codefiles or some other issue.

Finally, in regards to the estimation controls, these are directly taken from our runs on an HPC cluster, these runs were stress-testing a data-generating algorithm, hence the large number of samples/warmups. The number of questions, attributes and subjects is fixed based on the data generation, the actual MCMC controls are based on values that give us convergence - so unfortunately, no changing those.

Thank your for your continued help! I very much appreciate it.

Best,
Kelson

[goedman]

Hi @KelsonGs

As I also mentioned in my direct email, this morning I ran your model and was able to reproduce similar Rhat values for your model. Looking at some chain plots and the summary results I didn't see anything clearly wrong.

So I decided to experiment a bit with a slightly simpler model (see below pdf file) and what I found is that as soon as the Stan model combines coefficients, Rhat is reported as 1.0.

In below file this is demonstrated with 3 models estimating coefficients [a, b, sigma], [a, b[1], b[2], sigma] and [a, b, c, sigma]. In particular models 2 and 3 estimate similar values for b[1]==b and b[2]==c (as expected), but report Rhat differently.

🎈 kelson_02.pdf

The Rhat values are not visible, so I've copied them below:

ms_02.r_hat = [1.00044
1.00966
1.33853e13
1.01022
1.00864
NaN
1.00135
1.0011
1.0
1.0
1.0]

ms_03_r_hat = [1.00192
1.00349
1.94892e13
1.00508
1.00621
NaN
1.00263
0.999963
1.0002
1.00016
1.00165]

[goedman]

This is a script that can be used in the REPL:

using Pkg
using LinearAlgebra, HDF5
using StanSample

stan_02 = "
data {
	int<lower=1> N; 
	vector[N] x;
	vector[N] y;
	vector[N] z;
}
parameters {
	real a;
	vector[2] b;
	real<lower=0> sigma;
}
model {
	vector[N] mu;
	mu = a + b[1] * x + b[2] * y;
	a ~ normal(0, 2);
	b ~ normal(0, 1);
	sigma ~ exponential(1);
	z ~ normal(mu, sigma);
}
";

stan_03 = "
data {
	int<lower=1> N; 
	vector[N] x;
	vector[N] y;
	vector[N] z;
}
parameters {
	real a;
	real b;
	real c;
	real<lower=0> sigma;
}
model {
	vector[N] mu;
	mu = a + b * x + c * y;
	a ~ normal(0, 2);
	b ~ normal(0, 1);
	c ~ normal(0, 1);
	sigma ~ exponential(1);
	z ~ normal(mu, sigma);
}
";

N = 20
x = rand([-1, 1], N)
y = rand([-1, 1], N)
z = 1.0 .+ 5.0 .* x .+ 3.0 .* y .+ rand(Normal(1, 1), N)
data = (N=N, x=x, y=y, z=z)

path = expanduser("~/.julia/dev/StanSample/test/test_Kelson")
tmpdir = joinpath(path, "tmp")
isdir(tmpdir) && rm(tmpdir; recursive=true)

sm_02 = SampleModel("stan_02", stan_02, tmpdir)
rc_02 = stan_sample(sm_02; data)

if success(rc_02)
  ms_02 = read_summary(sm_02)
  ms_02.r_hat |> display
end

sm_03 = SampleModel("stan_03", stan_03, tmpdir)
rc_03 = stan_sample(sm_03; data)

if success(rc_03)
  ms_03 = read_summary(sm_03)
  ms_03.r_hat |> display
end

[goedman]

Another option would be to use MCMCChains.jl:

In [16]: chns = read_samples(sm_02, :mcmcchains)
Chains MCMC chain (1000×4×4 Array{Float64, 3}):

Iterations        = 1:1:1000
Number of chains  = 4
Samples per chain = 1000
parameters        = a, b.1, b.2, sigma
internals         = 

Summary Statistics
  parameters      mean       std      mcse    ess_bulk    ess_tail      rhat   ess_per_sec 
      Symbol   Float64   Float64   Float64     Float64     Float64   Float64       Missing 

           a    2.1674    0.2671    0.0047   3315.1546   2467.1093    1.0008       missing
         b.1    4.5323    0.2913    0.0064   2565.2265   1576.2727    1.0028       missing
         b.2    2.8664    0.2858    0.0062   2580.3754   1745.4285    1.0005       missing
       sigma    1.1244    0.2481    0.0063   2020.8906   1465.7429    1.0028       missing

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

           a    1.6480    1.9996    2.1620    2.3348    2.6989
         b.1    3.8605    4.3807    4.5569    4.7236    5.0248
         b.2    2.2347    2.6990    2.8931    3.0542    3.3602
       sigma    0.7713    0.9576    1.0832    1.2445    1.7235

This is after running above MWE.

[goedman]

Hi Kelson ( @KelsonGs ),

This is as far I think I can get. Not shown here, I tried the same with your model and get similar results (Rhat values close to 1). I'm still a bit concerned about the sigma estimates in your model (e.g. using plot_chains())

We could ask about the Rhat issue on the Stan mailing list.

Rob

[goedman]

Hi @KelsonGs

Can I close this discussion?

Thanks,
Rob