easystats / report

:scroll: :tada: Automated reporting of objects in R
https://easystats.github.io/report/
Other
694 stars 70 forks source link

ESS is reported inconsistently for `brmsfit` objects #343

Open gserapio opened 1 year ago

gserapio commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

Actually bayestestR, I think.

image

DominiqueMakowski commented 1 year ago

"shit", dom muttered