False warning about max Rhat? - rstan

General
1

With the new rstan I get a warning

Warning messages:
1: The largest R-hat is 1.09, indicating chains have not mixed.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#r-hat 
2: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help.

However when I sort my summary fit for Rhat the max Rhat does not match the warning

 par                 mean se_mean     sd  `2.5%`  `25%`  `50%`  `75%` `97.5%` n_eff  Rhat
   <chr>              <dbl>   <dbl>  <dbl>   <dbl>  <dbl>  <dbl>  <dbl>   <dbl> <dbl> <dbl>
 1 exposure_rate[29]  0.212 0.00706 0.0626  0.0934  0.166  0.211  0.253   0.352  78.6  1.04
 2 lambda_log[34]     7.65  0.00642 0.0624  7.54    7.60   7.65   7.68    7.78   94.5  1.04
 3 lambda_log[177]    6.07  0.00702 0.0954  5.88    6.00   6.07   6.13    6.26  185.   1.03
 4 lambda_log[1]      5.70  0.00624 0.0693  5.58    5.65   5.70   5.74    5.83  123.   1.02
 5 sigma[23]          0.659 0.00677 0.0937  0.481   0.594  0.655  0.718   0.861 191.   1.02
 6 lambda_log[68]     6.89  0.00190 0.0350  6.82    6.87   6.89   6.92    6.96  339.   1.02
 7 lambda_log[91]     6.10  0.00475 0.0639  5.95    6.06   6.10   6.14    6.23  181.   1.02
 8 lambda_log[12]     7.01  0.00333 0.0597  6.89    6.97   7.01   7.05    7.14  322.   1.02
 9 exposure_rate[13] -0.384 0.00370 0.0506 -0.473  -0.420 -0.383 -0.348  -0.287 187.   1.02
10 lambda_log[84]     7.16  0.00401 0.0660  7.04    7.12   7.16   7.21    7.28  271.   1.02

And if I plot the least mixed parameter

image
image.jpg955×929 102 KB

Model


data {

	// Reference matrix inference

	int<lower=0> G;
	int<lower=0> S;
	int CL;

 	int counts_linear[CL] ;
	int G_linear[CL] ;
	int S_linear[CL] ;

	real<upper=0> sigma_slope;
  real sigma_intercept;
  real<lower=0>sigma_sigma;


}
parameters {

  // Gene-wise properties of the data
  vector[G] lambda_log;
  vector[G] sigma;

  vector[S] exposure_rate;

	real<lower=0> lambda_mu;
  real<lower=0> lambda_sigma;
  real<upper=0> lambda_skew;

}

model {

  // Overall properties of the data
  lambda_mu ~ normal(0,2);
	lambda_sigma ~ normal(0,2);
	lambda_skew ~ normal(0,1);

exposure_rate ~ normal(0,1);
sum(exposure_rate) ~ normal(0, 0.001 * S);
lambda_log ~ skew_normal(lambda_mu, lambda_sigma, lambda_skew);
sigma ~ normal(sigma_slope * lambda_log + sigma_intercept, sigma_sigma);

counts_linear ~ neg_binomial_2_log(lambda_log[G_linear] + exposure_rate[S_linear], 1.0 ./ exp(sigma[G_linear]));

}

sessionInfo()

R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS:   /stornext/System/data/apps/R/R-3.6.0/lib64/R/lib/libRblas.so
LAPACK: /stornext/System/data/apps/R/R-3.6.0/lib64/R/lib/libRlapack.so

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

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

other attached packages:
 [1] data.tree_0.7.8       limma_3.40.2          foreach_1.4.4         rstan_2.19.2          StanHeaders_2.18.1-10 magrittr_1.5          forcats_0.4.0         stringr_1.4.0        
 [9] dplyr_0.8.3           purrr_0.3.2           readr_1.3.1           tidyr_0.8.3.9000      tibble_2.1.3          ggplot2_3.1.1         tidyverse_1.2.1      

loaded via a namespace (and not attached):
 [1] nlme_3.1-140              matrixStats_0.54.0        lubridate_1.7.4           RColorBrewer_1.1-2        httr_1.4.0                tools_3.6.0               backports_1.1.4          
 [8] utf8_1.1.4                R6_2.4.0                  lazyeval_0.2.2            colorspace_1.4-1          withr_2.1.2               tidyselect_0.2.5          gridExtra_2.3            
[15] prettyunits_1.0.2         processx_3.3.1            compiler_3.6.0            cli_1.1.0                 rvest_0.3.4               arrayhelpers_1.0-20160527 xml2_1.2.0               
[22] influenceR_0.1.0          plotly_4.9.0              labeling_0.3              diptest_0.75-7            scales_1.0.0              callr_3.2.0               digest_0.6.20            
[29] pkgconfig_2.0.2           htmltools_0.3.6           htmlwidgets_1.3           rlang_0.4.0               readxl_1.3.1              rstudioapi_0.10           svUnit_0.7-12            
[36] visNetwork_2.0.6          generics_0.0.2            jsonlite_1.6              rgexf_0.15.3              inline_0.3.15             tidyTranscriptomics_0.1.0 loo_2.1.0                
[43] Matrix_1.2-17             Rcpp_1.0.1                munsell_0.5.0             fansi_0.4.0               viridis_0.5.1             yaml_2.2.0                edgeR_3.26.3             
[50] stringi_1.4.3             ggstance_0.3.1            pkgbuild_1.0.3            plyr_1.8.4                grid_3.6.0                parallel_3.6.0            listenv_0.7.0            
[57] crayon_1.3.4              lattice_0.20-38           haven_2.1.0               hms_0.4.2                 locfit_1.5-9.1            zeallot_0.1.0             ps_1.3.0                 
[64] pillar_1.4.2              igraph_1.2.4.1            widyr_0.1.1               codetools_0.2-16          stats4_3.6.0              XML_3.98-1.19             glue_1.3.1               
[71] packrat_0.5.0             downloader_0.4            data.table_1.12.2         modelr_0.1.4              vctrs_0.2.0               cellranger_1.1.0          gtable_0.3.0             
[78] future_1.13.0             assertthat_0.2.1          broom_0.5.2               coda_0.19-2               viridisLite_0.3.0         iterators_1.0.10          Rook_1.1-1               
[85] tidybayes_1.0.4           DiagrammeR_1.0.1          globals_0.12.4            ellipsis_0.2.0.1          brew_1.0-6
2

What about the lp__ quantity?

3

That link says to “Look at Bulk- and Tail-ESS for further information.” which is provided by the rstan::rhat function rather than summary. You can also access it via monitor.

4

Just to share the opposite case (summary gives higher Rhat than warning)

Ordering summary by Rhat (1.17) warning 1.09.

   par                    mean se_mean     sd   `2.5%`    `25%`   `50%`   `75%` `97.5%` n_eff  Rhat
   <chr>                 <dbl>   <dbl>  <dbl>    <dbl>    <dbl>   <dbl>   <dbl>   <dbl> <dbl> <dbl>
 1 lambda_log[3202]     2.97   0.250   1.47   -1.64     2.74     3.29    3.86    4.60    34.4  1.17
 2 lambda_log[1597]     3.24   0.103   0.731   0.883    2.99     3.31    3.62    4.21    49.9  1.14
 3 prop_3[16,12]        0.0218 0.00267 0.0195  0.00268  0.00862  0.0156  0.0286  0.0763  53.2  1.12
 4 prop_2[16,7]         0.0301 0.00277 0.0206  0.00795  0.0157   0.0244  0.0379  0.0891  55.3  1.12
 5 prop_3[16,9]         0.949  0.00277 0.0213  0.892    0.941    0.955   0.964   0.976   59.2  1.11
 6 prop_2[16,5]         0.951  0.00277 0.0213  0.894    0.943    0.956   0.965   0.977   59.2  1.11
 7 prop_3[4,13]         0.0155 0.00108 0.0103  0.00276  0.00838  0.0130  0.0199  0.0440  90.2  1.06
 8 sigma_inv_log[3202]  2.50   0.0739  0.757   0.735    2.10     2.55    2.98    3.84   105.   1.05
 9 lambda_log[3118]     1.29   0.0474  0.725  -0.0755   0.848    1.27    1.69    2.82   234.   1.05
10 lambda_log[2387]    -3.21   0.0932  1.09   -5.33    -3.89    -3.21   -2.54   -1.08   136.   1.05
# … with 10,361 more rows
Warning messages:
1: The largest R-hat is 1.09, indicating chains have not mixed.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#r-hat 
2: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#bulk-ess 
3: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
http://mc-stan.org/misc/warnings.html#tail-ess

Just for completeness, I show the traceplots of the “worst” performing parameters

image

I know that 250 iterations are not enough, and the high Rhat depends on this as well. However, the inconsistency between summary and warning is a bit confusing.