Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ESS is reported inconsistently for brmsfit objects #343

Open
gserapio opened this issue Feb 28, 2023 · 3 comments
Open

ESS is reported inconsistently for brmsfit objects #343

gserapio opened this issue Feb 28, 2023 · 3 comments
Labels
bug 🐛 Something isn't working

Comments

@gserapio
Copy link

Describe the bug
ESS is being reported inconsistently, at least for brmsfit objects. In the output of report(), ESS values seem to correspond to the estimate of the row above. For example, the ESS for b Intercept[2] in report() would be 5820, but this is actually the ESS for b_Intercept[1] according to bayestestR::effective_sample().

Related question: how is ESS calculated from Bulk_ESS and Tail_ESS normally shown in summary outputs?

To Reproduce
Here is my (truncated) reprex:

library(brms)
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.18.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#> 
#> Attaching package: 'brms'
#> The following object is masked from 'package:stats':
#> 
#>     ar
library(easystats)
#> # Attaching packages: easystats 0.6.0
#> ✔ bayestestR  0.13.0   ✔ correlation 0.8.3 
#> ✔ datawizard  0.6.5    ✔ effectsize  0.8.3 
#> ✔ insight     0.19.0   ✔ modelbased  0.8.6 
#> ✔ performance 0.10.2   ✔ parameters  0.20.2
#> ✔ report      0.5.6    ✔ see         0.7.4

fit <- brm(rating ~ period + carry + cs(treat),
           data = inhaler, family = sratio("logit"),
           prior = set_prior("normal(0,5)"),
           chains = 5, cores = 5, seed = 42,
           backend = "cmdstanr", threads = threading(2))
#> Start sampling
#> Running MCMC with 5 parallel chains, with 2 thread(s) per chain...
#> 
#> [...]
#> 
#> All 5 chains finished successfully.
#> Mean chain execution time: 4.2 seconds.
#> Total execution time: 4.4 seconds.

summary(fit)
#>  Family: sratio 
#>   Links: mu = logit; disc = identity 
#> Formula: rating ~ period + carry + cs(treat) 
#>    Data: inhaler (Number of observations: 572) 
#>   Draws: 5 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 5000
#> 
#> Population-Level Effects: 
#>              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept[1]     0.55      0.09     0.38     0.73 1.00     5840     3735
#> Intercept[2]     2.37      0.29     1.86     2.97 1.00     3522     2875
#> Intercept[3]     0.48      0.59    -0.65     1.70 1.00     3946     3139
#> period           0.22      0.17    -0.10     0.54 1.00     5583     4143
#> carry           -0.21      0.17    -0.55     0.11 1.00     3541     3629
#> treat[1]        -0.79      0.24    -1.25    -0.33 1.00     3815     3731
#> treat[2]        -1.05      0.60    -2.29     0.05 1.00     3391     3112
#> treat[3]         1.16      1.19    -1.16     3.50 1.00     4094     2775
#> 
#> Family Specific Parameters: 
#>      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> disc     1.00      0.00     1.00     1.00   NA       NA       NA
#> 
#> Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

report_fit <- report(fit)
#> Warning: Predictions are treated as continuous variables in 'bayes_R2' which is
#> likely invalid for ordinal families.
#> Start sampling
#> Running MCMC with 5 sequential chains, with 2 thread(s) per chain...
#> 
#> [...]
#> 
#> All 5 chains finished successfully.
#> Mean chain execution time: 2.4 seconds.
#> Total execution time: 12.2 seconds.
#> Warning: Predictions are treated as continuous variables in 'bayes_R2' which is
#> likely invalid for ordinal families.

Here, ESS for the first three effect estimates are reported as 3462, 5820, 3281:

report_fit
#> We fitted a Bayesian logistic model (estimated using MCMC sampling with 5
#> chains of 2000 iterations and a warmup of 1000) to predict rating with period,
#> carry and treat (formula: rating ~ period + carry + cs(treat)). Priors over
#> parameters were set as normal (mean = 0.00, SD = 5.00) distributions. The
#> model's explanatory power is weak (R2 = 0.06, 95% CI [0.03, 0.10]).  Within
#> this model:
#> 
#>   - The effect of b Intercept[1] (Median = 0.55, 95% CI [0.38, 0.73]) has a
#> 100.00% probability of being positive (> 0), 100.00% of being significant (>
#> 0.05), and 99.80% of being large (> 0.30). The estimation successfully
#> converged (Rhat = 1.001) and the indices are reliable (ESS = 3462)
#>   - The effect of b Intercept[2] (Median = 2.36, 95% CI [1.86, 2.97]) has a
#> 100.00% probability of being positive (> 0), 100.00% of being significant (>
#> 0.05), and 100.00% of being large (> 0.30). The estimation successfully
#> converged (Rhat = 1.000) and the indices are reliable (ESS = 5820)
#>   - The effect of b Intercept[3] (Median = 0.48, 95% CI [-0.65, 1.70]) has a
#> 79.66% probability of being positive (> 0), 77.12% of being significant (>
#> 0.05), and 61.84% of being large (> 0.30). The estimation successfully
#> converged (Rhat = 1.000) and the indices are reliable (ESS = 3281)
#>   - The effect of b period (Median = 0.22, 95% CI [-0.10, 0.54]) has a 90.76%
#> probability of being positive (> 0), 83.92% of being significant (> 0.05), and
#> 31.48% of being large (> 0.30). The estimation successfully converged (Rhat =
#> 1.001) and the indices are reliable (ESS = 3841)
#>   - The effect of b carry (Median = -0.21, 95% CI [-0.55, 0.11]) has a 89.48%
#> probability of being negative (< 0), 83.10% of being significant (< -0.05), and
#> 29.68% of being large (< -0.30). The estimation successfully converged (Rhat =
#> 1.001) and the indices are reliable (ESS = 5563)
#>   - The effect of bcs treat[1] (Median = -0.78, 95% CI [-1.25, -0.33]) has a
#> 99.94% probability of being negative (< 0), 99.88% of being significant (<
#> -0.05), and 97.98% of being large (< -0.30). The estimation successfully
#> converged (Rhat = 1.000) and the indices are reliable (ESS = 3759)
#>   - The effect of bcs treat[2] (Median = -1.03, 95% CI [-2.29, 0.05]) has a
#> 96.66% probability of being negative (< 0), 95.80% of being significant (<
#> -0.05), and 90.28% of being large (< -0.30). The estimation successfully
#> converged (Rhat = 1.000) and the indices are reliable (ESS = 3179)
#>   - The effect of bcs treat[3] (Median = 1.16, 95% CI [-1.16, 3.50]) has a 84.10%
#> probability of being positive (> 0), 83.02% of being significant (> 0.05), and
#> 76.90% of being large (> 0.30). The estimation successfully converged (Rhat =
#> 1.001) and the indices are reliable (ESS = 3969)
#> 
#> [...]

In as.report_table() format, though, ESS is reported corresponding to bayestestR::effective_sample():

as.report_table(report_fit)
#> Parameter    | Median |         95% CI |     pd |  Rhat |     ESS |     Fit
#> ---------------------------------------------------------------------------
#> Intercept[1] |   0.55 | [ 0.38,  0.73] |   100% | 1.000 | 5820.00 |        
#> Intercept[2] |   2.36 | [ 1.86,  2.97] |   100% | 1.000 | 3281.00 |        
#> Intercept[3] |   0.48 | [-0.65,  1.70] | 79.66% | 1.001 | 3841.00 |        
#> period       |   0.22 | [-0.10,  0.54] | 90.76% | 1.001 | 5563.00 |        
#> carry        |  -0.21 | [-0.55,  0.11] | 89.48% | 1.001 | 3462.00 |        
#> treat[1]     |  -0.78 | [-1.25, -0.33] | 99.94% | 1.000 | 3759.00 |        
#> treat[2]     |  -1.03 | [-2.29,  0.05] | 96.66% | 1.000 | 3179.00 |        
#> treat[3]     |   1.16 | [-1.16,  3.50] | 84.10% | 1.001 | 3969.00 |        
#>              |        |                |        |       |         |        
#> ELPD         |        |                |        |       |         | -459.66
#> LOOIC        |        |                |        |       |         |  919.33
#> WAIC         |        |                |        |       |         |  919.16
#> R2           |        |                |        |       |         |    0.06

bayestestR::effective_sample() ESS values:

effective_sample(fit)
#>        Parameter  ESS
#> 1 b_Intercept[1] 5820
#> 2 b_Intercept[2] 3281
#> 3 b_Intercept[3] 3841
#> 4       b_period 5563
#> 5        b_carry 3462
#> 6   bcs_treat[1] 3759
#> 7   bcs_treat[2] 3179
#> 8   bcs_treat[3] 3969

Created on 2023-02-28 with reprex v2.0.2

Expected behaviour
ESS values in the output of report() above should have been 5820, 3281, 3841, etc., corresponding with the output of bayestestR::effective_sample().

Specifications (please complete the following information):

library(brms)
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.18.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#> 
#> Attaching package: 'brms'
#> The following object is masked from 'package:stats':
#> 
#>     ar
library(easystats)
#> # Attaching packages: easystats 0.6.0
#> ✔ bayestestR  0.13.0   ✔ correlation 0.8.3 
#> ✔ datawizard  0.6.5    ✔ effectsize  0.8.3 
#> ✔ insight     0.19.0   ✔ modelbased  0.8.6 
#> ✔ performance 0.10.2   ✔ parameters  0.20.2
#> ✔ report      0.5.6    ✔ see         0.7.4

System info:

> Sys.info()
                                                                                               sysname 
                                                                                              "Darwin" 
                                                                                               release 
                                                                                              "22.3.0" 
                                                                                               version 
"Darwin Kernel Version 22.3.0: Mon Jan 30 20:38:37 PST 2023; root:xnu-8792.81.3~2/RELEASE_ARM64_T6000"

Hope this helps!

@rempsyc
Copy link
Sponsor Member

rempsyc commented Mar 2, 2023

report relies on several other packages under the hood to get all the numbers, so I wonder if this has to do with the parameters package.

@strengejacke
Copy link
Member

Actually bayestestR, I think.

image

@rempsyc rempsyc added the bug 🐛 Something isn't working label Mar 2, 2023
@DominiqueMakowski
Copy link
Member

"shit", dom muttered

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug 🐛 Something isn't working
Projects
None yet
Development

No branches or pull requests

4 participants