Return NA for Rhat, ESS, and MCSE if any of the chains is constant

[avehtari]

Currently, we return NA for Rhat, ESS, and MCSE if all draws have the same value as we can't know whether the case is structural constant or stuck chains. However, if only some chains are constant we don't return NA. This leads to overconfident ESS and MCSE estimates for extreme quantiles.

4 "chains" with different supports

q <- cbind(runif(n=1000,min=-1,max=1),runif(n=1000,min=-1,max=1),runif(n=1000,min=-1,max=1),runif(n=1000,min=-1,max=1)*1.25) 

produce following quantile ESSs. The dotted line shows the quantile 0.02375. Left from that quantile, 3 of the chains don't have any values and thus the indicator function is constant, and ESS starts to increase.

I propose (and @paul-buerkner agrees) that we should change the Rhat/ESS computation to return NA if any of the chains is constant. Then quantile plot would look like this

and no MCSE would be given for quantiles with alpha<0.02375

It would be good to add to the documentation explanation of this. T

[paul-buerkner]

Should now be fixed.

[mike-lawrence]

The ess_tail() function is now yielding NA for even it's own example.

mu <- extract_variable_matrix(example_draws(), "mu")
ess_tail(mu)
#> [1] NA

> sessionInfo()
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3

locale:
 [1] LC_CTYPE=en_CA.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_CA.UTF-8        LC_COLLATE=en_CA.UTF-8    
 [5] LC_MONETARY=en_CA.UTF-8    LC_MESSAGES=en_CA.UTF-8   
 [7] LC_PAPER=en_CA.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods  
[7] base     

other attached packages:
[1] posterior_1.0.0 aria_0.1.0     

loaded via a namespace (and not attached):
 [1] matrixStats_0.58.0   fs_1.5.0             usethis_2.0.1       
 [4] lubridate_1.7.10     devtools_2.4.1       httr_1.4.2          
 [7] rprojroot_2.0.2      rstan_2.21.2         tensorA_0.36.2      
[10] tools_4.1.0          backports_1.2.1      utf8_1.2.1          
[13] R6_2.5.0             DBI_1.1.1            colorspace_2.0-1    
[16] withr_2.4.2          tidyselect_1.1.1     gridExtra_2.3       
[19] prettyunits_1.1.1    processx_3.5.2       curl_4.3.1          
[22] compiler_4.1.0       rvest_1.0.0          cli_2.5.0           
[25] RNetCDF_2.4-2        xml2_1.3.2           desc_1.3.0          
[28] stringfish_0.15.1    scales_1.1.1         checkmate_2.0.0     
[31] readr_1.4.0          callr_3.7.0          stringr_1.4.0       
[34] digest_0.6.27        StanHeaders_2.21.0-7 pkgconfig_2.0.3     
[37] sessioninfo_1.1.1    dbplyr_2.1.1         fastmap_1.1.0       
[40] readxl_1.3.1         rlang_0.4.11         rstudioapi_0.13     
[43] generics_0.1.0       farver_2.1.0         RApiSerialize_0.1.0 
[46] jsonlite_1.7.2       dplyr_1.0.6          distributional_0.2.2
[49] inline_0.3.18        magrittr_2.0.1       loo_2.4.1           
[52] Rcpp_1.0.6           munsell_0.5.0        fansi_0.5.0         
[55] abind_1.4-5          lifecycle_1.0.0      stringi_1.6.2       
[58] pkgbuild_1.2.0       grid_4.1.0           parallel_4.1.0      
[61] forcats_0.5.1        crayon_1.4.1         haven_2.4.1         
[64] hms_1.1.0            knitr_1.33           ps_1.6.0            
[67] pillar_1.6.1         codetools_0.2-18     stats4_4.1.0        
[70] pkgload_1.2.1        reprex_2.0.0         glue_1.4.2          
[73] beepr_1.3            V8_3.4.2             modelr_0.1.8        
[76] data.table_1.14.0    remotes_2.3.0        renv_0.13.2         
[79] RcppParallel_5.1.4   vctrs_0.3.8          cellranger_1.1.0    
[82] testthat_3.0.2       gtable_0.3.0         purrr_0.3.4         
[85] tidyr_1.1.3          qs_0.24.1            assertthat_0.2.1    
[88] cachem_1.0.5         ggplot2_3.3.3        xfun_0.23           
[91] broom_0.7.6          tidyverse_1.3.1      roxygen2_7.1.1      
[94] audio_0.1-7          tibble_3.1.2         memoise_2.0.0       
[97] ellipsis_0.3.2       cmdstanr_0.4.0.9000

posterior GithubSHA1: ee38b68

[paul-buerkner]

I agree this behavior is unfortunate but consistent with the above described notion of how to handle variables that are constant within a chain. I have also come across such examples recently, and I don't know yet, how to best move forward with this. Any ideas and suggestions are very welcome.

[mike-lawrence]

how to handle variables that are constant within a chain.

But mu isn't constant in any of the chains

[mike-lawrence]

Ah, delving into the computations that ess_tail implies, it turns the samples into binary indicators of whether the sample is below X (where X is the 5%ile value of all samples, ignoring chain membership), and then each chain is split in two, so some of these half-chains can get a constant FALSE if they have few values below X.

And even for a single chain, if there's a difference in dispersion between the halves the less dispersed half will have a greater chance of having all FALSE values. So I guess NA for ess_tail is a good signal of non-convergence.

[paul-buerkner]

Yes this is the logic. So with relatively few iterations per chain and lots of chains, ess_tail will probably be NA more often than not. I agree it is not ideal thought. NA does not necessarily indicate non-convergence but rather an unclear state that could be caused by multiple things that are hard to differentiate, not all of which are a reason for concern (e.g., a variable that is constant by design).

[mjskay]

Could this be turned into a potentially useful teaching example in the docs? Something along the lines of "this is why you shouldn't calculate 90% intervals on a sample this small" followed by some simple examples, like ess_tail(rnorm(50)); ess_tail(rnorm(100)); ess_tail(rnorm(1000))?

I'm not sure how you can provide more useful information to users without raising errors or warnings, which seem likely to cause lots of false positives. We could create a subtype of numeric vectors that lets you provide some metadata on NAs, but it could get messy and I think would probably not be worth it.