next release

[weberse2]

The main point about the upcoming release is planned to be the sum to zero parametrization trick for hierarchical models. This removes a parameter from the over-parametrized hierarchical model and has shown in benchmarks to be easier to the sampler. Essentially I can drop the default adapt_delta setting for NUTS from 0.99 to 0.95 for common models fitted with RBesT. Below is are instructions for running the attached benchmark, which is a fake data simulation fitted with priors which are rather weak, but used in practice. This benchmark used to be a challenge for RBesT which is now not any longer the case as one can see from a medium sized job - the new implementation has essentially no more divergent transitions while the old one struggles:

=====================================================
divergent-MAP stress test: legacy vs s2z
=====================================================
trials: 1000 | chains: 4 | J = 3 groups

Table: divergences and cost, arms side by side

setting      trials   delta legacy   delta s2z   div% legacy   div% s2z   div/fit legacy   div/fit s2z   ESS/kgrad legacy   ESS/kgrad s2z   sec legacy   sec s2z
----------  -------  -------------  ----------  ------------  ---------  ---------------  ------------  -----------------  --------------  -----------  --------
shipped        1000          0.990       0.950          39.0        0.3             1.01             0               7.00           37.19          1.3       0.9
matched99      1000          0.990       0.990          37.4        0.0             0.85             0               7.04           23.43          1.3       1.3
matched95      1000          0.950       0.950          92.7        0.5             8.45             0              10.27           37.50          0.8       0.9
tuned          1000          0.990       0.990          50.1        0.0             1.30             0               8.46           27.71          1.2       1.2
hard           1000          0.999       0.999           1.4        0.0             0.05             0               4.93           17.03          4.0       4.4

Instructions and code here (Opus vibe coded):

benchmark-divergent-map_R.txt

Reproducing the s2z vs. legacy benchmark

The benchmark-divergent-map.R script in the package root compares the new
sum-to-zero (s2z) parametrization of gMAP() against the legacy centred one on
the "divergent MAP" fake data simulation: a two-arm trial with a
covariate-by-treatment interaction, analysed as a random-effects meta-analysis
of the within-level treatment contrasts. With J = 2 groups tau is informed by
two numbers only, which leaves the long funnel that makes this a fair stress
test rather than a favourable one.

Both arms come from the same installed model — the legacy arm sets the user
facing opt-out (use_s2z = 0, i.e. options(RBesT.MC.s2z = FALSE)), the s2z arm is
the default — driven with identical data, seeds and control arguments. So one
install of this branch is all that is needed, and nothing has to be compiled per
arm.

Run it

# 1. install this branch into build/installed/
make dev-install

# 2. quick local check (~1 min): 4 trials, each arm at its shipped adapt_delta
Rscript benchmark-divergent-map.R --trials 4 --settings shipped --lib build/installed

# 3. all options, settings and cluster recipes
Rscript benchmark-divergent-map.R --help

Run it from the package root. --lib build/installed puts the dev install on
the library path; drop it if that library is already on R_LIBS_USER.

Useful variations

# use all cores of a workstation
Rscript benchmark-divergent-map.R --trials 50 \
  --backend mirai --daemons 12 --lib build/installed

# geometry comparison at a matched target (0.99 is the value legacy was tuned at)
Rscript benchmark-divergent-map.R --trials 50 --settings matched99 --lib build/installed

# does the advantage survive better informed heterogeneity?
Rscript benchmark-divergent-map.R --trials 50 --groups 5 \
  --backend mirai --daemons 12 --lib build/installed

# where should the default adapt_delta of the new geometry sit?
Rscript benchmark-divergent-map.R --trials 500 --arms s2z \
  --settings target95,target975,target99 --backend mirai --daemons 12

# cluster: LSF out of the box, mirai/clustermq/SLURM recipes at the end of the script
Rscript benchmark-divergent-map.R --trials 2000 --backend lsf --jobs 400 \
  --modules 'module load R/4.4.1'

Settings ladder: shipped (each arm at its own default, 0.99 legacy / 0.95
s2z), matched99, matched95, tuned, hard. The target* settings vary
only adapt_delta and are meant for --arms s2z alone.

A handful of trials cannot separate a 2% from a 5% divergence rate, so keep
local runs as smoke tests and use a worker pool for the numbers you want to
quote. Results are bit-identical across backends: trial data and Stan seed are
determined by the grid row.

Output

The script prints, per setting and arm, divergence rates, ESS per 1000
gradients, runtime and a paired per-trial contrast, and writes the raw results
to benchmark-divergent-map.rds (override with --out). Re-report a saved run
without refitting:

source("benchmark-divergent-map.R")
dmap_report(readRDS("benchmark-divergent-map.rds"))

The headline claim to check is that the s2z row at a mild adapt_delta matches
or beats the legacy row at a harsher one.

[weberse2]

I also now tried the centering and auto centering options as implemented in the respective brms branch... and the auto centering technique seems to buy further improvements. These are not yet in this PR, but I will add these once I have them documented, tested and have my head around it - which is worth it as it looks to me.

=====================================================
divergent-MAP stress test: legacy vs s2z (s2z, center, auto)
=====================================================
trials: 2000 | chains: 4 | J = 4 groups

Table: divergences, arms side by side

setting      trials   div% auto   div% center   div% legacy   div% s2z   div/fit auto   div/fit center   div/fit legacy   div/fit s2z
----------  -------  ----------  ------------  ------------  ---------  -------------  ---------------  ---------------  ------------
shipped        2000           0          99.6          19.3        0.5              0           100.25             0.39             0
matched99      2000           0          99.0          18.5        0.0              0            37.47             0.30             0
matched95      2000           0          99.4          75.6        0.4              0            95.92             3.64             0
tuned          2000           0          99.0          29.8        0.0              0            41.95             0.56             0
hard           2000           0          94.7           1.0        0.0              0            19.08             0.01             0

Table: cost, arms side by side

setting      trials   ESS/kgrad auto   ESS/kgrad center   ESS/kgrad legacy   ESS/kgrad s2z   sec auto   sec center   sec legacy   sec s2z
----------  -------  ---------------  -----------------  -----------------  --------------  ---------  -----------  -----------  --------
shipped        2000           106.04              42.41               9.56           42.83        0.8          1.4          1.3       1.0
matched99      2000            49.34              20.05               9.46           26.18        1.4          2.9          1.3       1.5
matched95      2000           105.22              42.53              14.51           42.88        0.8          1.4          0.9       1.0
tuned          2000            68.99              23.52              11.63           31.12        1.3          2.6          1.2       1.4
hard           2000            31.62               9.76               6.48           19.43        6.3          9.3          4.5       5.3