Add percentage of divergent transitions to convergence flag

[nehill197]

Should add divergent transitions as indicator of convergence.

Stan manual suggests a small number of divergent transitions is "often good enough to move forward" if you have good rhat and ess values.

Look into how divergent transitions are tracked within cmdstanr - if any are flagged (regardless of whether or not that sample is one of the thinned samples) it makes more sense to look at as a percentage.

Initial suggestion is < 1% but we should also research to see if there is an accepted %.

[StefanoMezzini]

the Stan page on common warnings has some helpful info. Specifically:

Even a small number of divergences after warmup cannot be safely ignored if completely reliable inference is desired. But if you get only few divergences and you get good Rhat and ESS values, the resulting posterior is often good enough to move forward. (I added the emphasis)

[StefanoMezzini]

There are though some people who are skeptical of any model with one or more divergences, and they say that divergences are only ok in early stages (if even) but should not be in final models. It should be possible to track all divergences with Stan because {mvgam} gives a rug plot of them with the trace plots

[StefanoMezzini]

another helpful resource: https://link.springer.com/article/10.1007/s11222-026-10825-9

For example, it is known that there are also false positive divergence warnings, and in such cases posterior SBC can be used to measure the amount of potential bias.

[StefanoMezzini]

It's hard to find a clear reference for what Stan considers a divergence. This is the best I've found so far, but it's not clear to me what the scales are or why the numbers are what they are:

The current default threshold is a factor of $10^3$, whereas when the leapfrog integrator is working properly, the divergences will be around $10^{−7}$ and do not compound due to the symplectic nature of the leapfrog integrator.

[nehill197]

Hey @StefanoMezzini thanks for starting to look into this!
Are you able to figure out whether the output from cmdstanr that eventually ends up in our glance object includes the divergent transitions only in the samples retained (i.e., excludes divergent transitions in the samples that were thinned out)?

[StefanoMezzini]

Yep! working on that now. I'm also looking at {mvgam} to see how they work with them

[StefanoMezzini]

the model always keeps all iterations -- this is visible in the glance object, since there are 284 divergences for the model with a total of 3 * 100 = 300 samples despite nthin = 10.

      n     K nchains niters nthin   ess  rhat converged num_divergent max_treedepth ebfmi
  <int> <int>   <int>  <int> <int> <int> <dbl> <lgl>             <dbl>         <int> <dbl>
1    40     5       3    100    10    11  675. FALSE               284             0  1.96

[StefanoMezzini]

I've added a column of the percentage of divergence transitions in the glance() output. The warnings (actually messages) were already there

> glance(x)
Warning: 71 of 3000 (2.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.

Warning: 2853 of 3000 (95.0%) transitions hit the maximum treedepth limit of 10.
See https://mc-stan.org/misc/warnings for details.
# A tibble: 1 × 13
      n     K nchains niters nthin   ess  rhat converged num_divergent perc_divergent max_treedepth perc_max_treedepth ebfmi
  <int> <int>   <int>  <int> <int> <int> <dbl> <lgl>             <dbl>          <dbl>         <int>              <dbl> <dbl>
1    40     3       3   1000     1    12  4.31 FALSE                71           2.37          2853               95.1  1.92

did we want to edit converged to force it to be FALSE if num_divergent > 0?

[nehill197]

That's great, thanks for looking into that.
I think we want to pick a perc_divergent value to be the cutoff, not a number.
I also think we may want to just keep the perc_divergent column in the table and exclude the num_divergent, because this one is getting quite wide and will be hard to squeeze in a single row in a report.
Did you find anything in your research that suggested a value that would be acceptable? I was thinking 1% originally but open to other options.

[StefanoMezzini]

I think it depends on the model complexity and desired accuracy/precision, although I am increasingly convinced that models should not have any divergences. It seems like we should be able to avoid false positives by increasing adapt_delta to decrease the step size and use a better stepwise linear approximation of Hamiltonian dynamics (see this page). Currently, I am leaning towards setting converged = FALSE whenever num_divergence > 1.

The divergences could be in part due to our frequent use of priors with modes near a boundary (e.g., exp, half normal)when likelihoods are unable to move the posterior mode away from the boundary. Conceptually this makes sense because a posterior with greatest density near 0 is sensible for some regression coefficients but nonsensical for many variance parameters which we assume to be > 0.

[StefanoMezzini]

here's an interesting example: https://mc-stan.org/learn-stan/case-studies/divergences_and_bias.html

[StefanoMezzini]

A problem with using both number of or % of divergent transitions as a metric to assess convergence is that the occurrence of rare events (divergences) is very sensitive to chain length. Thus, conveniently short chains are more likely to converge than long chains, and small percentages and numbers should instead be inflated if the chains are short. the example above has a good illustration of this issue

[sebdalgarno]

thanks @StefanoMezzini - one suggestion is to make convergence in this case dependent on whether analysis mode is 'paper' or 'report'. in case of 'report', we could use some threshold like @nehill197 1% and in 'paper' we could enforce no divergences? I think in many cases the divergences can be fixed (or mostly fixed) by adjusting adapt_delta - so this is a good way to enforce people to try that.
1% might be too lenient though - that is 20 divergences in normal model fit of 500 iters x 4 chains. 0.1%, or 0.2% might be more appropriate?