# Simulation based calibration for RBesT

#### Sebastian Weber

#### Tue Nov 23 23:08:09 2021

This report documents the results of a simulation based calibration (SBC) run for `RBesT`. The calibration data will be generated whenever relevant changes to the `gMAP` function were made. The calibration runs are performed for typical use cases of `gMAP`. These include the three likelihoods (binomial, gaussian & Poisson), a sparse (\(2\) trials) and dense (\(10\) trials) data situation and finally a run with a very/less conservative prior choice for between-trial heterogeniety parameter.

The calibration data presented here has been generated at and with the `RBesT` git version as:

    ## Created:  2021-09-03 08:43:14 UTC
    ## git hash: 3b254d41091496843b68a39477b149202796a5f5
    ## MD5:      8f363c657046968a77eee164222f45dd

The MD5 hash of the calibration data file presented here must match the above listed MD5:

    ##                    calibration.rds
    ## "8f363c657046968a77eee164222f45dd"

# Introduction

Simulation based calibration (SBC) is a necessary condition which must be met for any Bayesian analysis with proper priors. The details are presented in Talts, et. al (see https://arxiv.org/abs/1804.06788).

Self-consistency of any Bayesian analysis with a proper prior:

\[ p(\theta) = \iint \mbox{d}\tilde{y} \, \mbox{d}\tilde{\theta} \, p(\theta|\tilde{y}) \, p(\tilde{y}|\tilde{\theta}) \, p(\tilde{\theta}) \] \[ \Leftrightarrow p(\theta) = \iint \mbox{d}\tilde{y} \, \mbox{d}\tilde{\theta} \, p(\theta,\tilde{y},\tilde{\theta}) \]

SBC procedure:

Repeat \(s=1, ..., S\) times:

1. Sample from the prior \[\tilde{\theta} \sim p(\theta)\]

2. Sample fake data \[\tilde{y} \sim p(y|\tilde{\theta})\]

3. Obtain \(L\) posterior samples \[\{\theta_1, ..., \theta_L\} \sim p(\tilde{\theta}|\tilde{y})\]

4. Calculate the rank \(r_s\) of the prior draw \(\tilde{\theta}\) wrt to the posterior sample \(\{\theta_1, ..., \theta_L\} \sim p(\tilde{\theta}|\tilde{y})\) which falls into the range \([0,L]\) out of the possible \(L+1\) ranks. The rank is calculated as \[r_s = \sum_{l=1}^L \mathbb{I}[ \theta_l < \tilde{\theta}]\]

The \(S\) ranks then form a uniform \(0-1\) density and the count in each bin has a binomial distribution with probability of \[p(r \in \mbox{Any Bin}) =\frac{(L+1)}{S}.\]

## Hierarchical intercept only (random-effects intercept) model for binomial, gaussian and Poisson likelihood

Likelihood:

- Binary \[y_i|\theta_{j} \sim \mbox{Bernoulli}(\theta_j), \] \[g(\theta) = \mbox{logit}(\theta)\]
- Normal \[y_i|\theta_{j} \sim \mbox{Normal}(\theta_j, \sigma^2), \] \[g(\theta) = \theta\]
- Poisson \[y_i|\theta_{j} \sim \mbox{Poisson}(\theta_j), \] \[g(\theta) = \log(\theta)\]

Hierarchical prior:

\[ g(\theta_j)|\mu,\tau \sim \mbox{Normal}(\mu, \tau^2)\]

\[\mu \sim \mbox{Normal}(m_\mu, s^2_\mu)\] \[\tau \sim \mbox{Normal}^+(0, s^2_\tau)\]

The fake data simulation function returns for binomial and Poisson data the sum of the responses while for normal the mean summary is used. Please refer to the `sbc_tools.R` and `make_reference_rankhist.R` R programs for the implementation details.

The reference runs are created with \(L=1023\) posterior draws for each replication and a total of \(S=10^4\) replications are run per case. For the evaluation here the results are reduced to \(B=L'+1=64\) bins to ensure a sufficiently large sample size per bin.

# SBC results

## Sampler Diagnostics Overview

family | problem | sd_tau | N | total_divergent | min_ess | max_Rhat | total_large_Rhat | min_lp_ess_bulk | min_lp_ess_tail
--- | --- | --- | --- | --- | --- | --- | --- | --- | ---
binomial | dense | 0.5 | 10000 | 0 | 917 | 1.012595 | 0 | 407 | 244
binomial | dense | 1.0 | 10000 | 0 | 571 | 1.012849 | 0 | 388 | 391
binomial | sparse | 0.5 | 10000 | 11 | 465 | 1.009765 | 0 | 390 | 261
binomial | sparse | 1.0 | 10000 | 11 | 162 | 1.012056 | 0 | 325 | 281
gaussian | dense | 0.5 | 10000 | 0 | 395 | 1.015950 | 0 | 343 | 319
gaussian | dense | 1.0 | 10000 | 3 | 436 | 1.019491 | 0 | 317 | 271
gaussian | sparse | 0.5 | 10000 | 1532 | 16 | 1.167858 | 0 | 41 | 31
gaussian | sparse | 1.0 | 10000 | 2831 | 2 | 2.434608 | 1 | 3 | 23
poisson | dense | 0.5 | 10000 | 1 | 451 | 1.011987 | 0 | 265 | 327
poisson | dense | 1.0 | 10000 | 0 | 384 | 1.012953 | 0 | 266 | 340
poisson | sparse | 0.5 | 10000 | 3018 | 5 | 1.714576 | 1 | 10 | 50
poisson | sparse | 1.0 | 10000 | 3326 | 6 | 1.531513 | 1 | 8 | 110

Note: Large Rhat is defined as exceeding 1.2.

# Summary Statistics

## \(\chi^2\) Statistic, \(\mu\)

| problem | likelihood | sd_tau | parameter | statistic | df | p.value
--- | --- | --- | --- | --- | --- | --- | ---
1 | dense | binomial | 0.5 | mu | 37.747 | 63 | 0.995
13 | dense | binomial | 1 | mu | 64.026 | 63 | 0.440
25 | dense | gaussian | 0.5 | mu | 65.011 | 63 | 0.406
37 | dense | gaussian | 1 | mu | 46.362 | 63 | 0.943
49 | dense | poisson | 0.5 | mu | 62.490 | 63 | 0.494
61 | dense | poisson | 1 | mu | 59.021 | 63 | 0.619
73 | sparse | binomial | 0.5 | mu | 62.298 | 63 | 0.501
77 | sparse | binomial | 1 | mu | 73.485 | 63 | 0.172
81 | sparse | gaussian | 0.5 | mu | 59.904 | 63 | 0.587
85 | sparse | gaussian | 1 | mu | 71.680 | 63 | 0.212
89 | sparse | poisson | 0.5 | mu | 57.024 | 63 | 0.688
93 | sparse | poisson | 1 | mu | 66.240 | 63 | 0.366

## \(\chi^2\) Statistic, \(\tau\)

| problem | likelihood | sd_tau | parameter | statistic | df | p.value
--- | --- | --- | --- | --- | --- | --- | ---
2 | dense | binomial | 0.5 | tau | 58.624 | 63 | 0.633
14 | dense | binomial | 1 | tau | 84.032 | 63 | 0.040
26 | dense | gaussian | 0.5 | tau | 69.683 | 63 | 0.263
38 | dense | gaussian | 1 | tau | 53.555 | 63 | 0.796
50 | dense | poisson | 0.5 | tau | 67.430 | 63 | 0.328
62 | dense | poisson | 1 | tau | 79.309 | 63 | 0.080
74 | sparse | binomial | 0.5 | tau | 51.635 | 63 | 0.846
78 | sparse | binomial | 1 | tau | 61.517 | 63 | 0.529
82 | sparse | gaussian | 0.5 | tau | 63.424 | 63 | 0.461
86 | sparse | gaussian | 1 | tau | 63.757 | 63 | 0.450
90 | sparse | poisson | 0.5 | tau | 51.251 | 63 | 0.855
94 | sparse | poisson | 1 | tau | 61.773 | 63 | 0.520

## \(\chi^2\) Statistic, group estimates \(\theta\)

| problem | likelihood | sd_tau | parameter | statistic | df | p.value
--- | --- | --- | --- | --- | --- | --- | ---
3 | dense | binomial | 0.5 | theta[1] | 69.709 | 63 | 0.262
4 | dense | binomial | 0.5 | theta[10] | 62.925 | 63 | 0.479
5 | dense | binomial | 0.5 | theta[2] | 49.830 | 63 | 0.886
6 | dense | binomial | 0.5 | theta[3] | 56.115 | 63 | 0.718
7 | dense | binomial | 0.5 | theta[4] | 59.789 | 63 | 0.591
8 | dense | binomial | 0.5 | theta[5] | 65.203 | 63 | 0.400
9 | dense | binomial | 0.5 | theta[6] | 73.779 | 63 | 0.166
10 | dense | binomial | 0.5 | theta[7] | 75.750 | 63 | 0.130
11 | dense | binomial | 0.5 | theta[8] | 59.763 | 63 | 0.592
12 | dense | binomial | 0.5 | theta[9] | 71.859 | 63 | 0.208
15 | dense | binomial | 1 | theta[1] | 59.648 | 63 | 0.597
16 | dense | binomial | 1 | theta[10] | 79.834 | 63 | 0.075
17 | dense | binomial | 1 | theta[2] | 62.016 | 63 | 0.511
18 | dense | binomial | 1 | theta[3] | 48.115 | 63 | 0.917
19 | dense | binomial | 1 | theta[4] | 52.698 | 63 | 0.819
20 | dense | binomial | 1 | theta[5] | 57.254 | 63 | 0.680
21 | dense | binomial | 1 | theta[6] | 66.970 | 63 | 0.343
22 | dense | binomial | 1 | theta[7] | 88.947 | 63 | 0.017
23 | dense | binomial | 1 | theta[8] | 53.261 | 63 | 0.804
24 | dense | binomial | 1 | theta[9] | 71.706 | 63 | 0.211
27 | dense | gaussian | 0.5 | theta[1] | 54.208 | 63 | 0.777
28 | dense | gaussian | 0.5 | theta[10] | 66.317 | 63 | 0.363
29 | dense | gaussian | 0.5 | theta[2] | 75.597 | 63 | 0.133
30 | dense | gaussian | 0.5 | theta[3] | 44.134 | 63 | 0.966
31 | dense | gaussian | 0.5 | theta[4] | 54.784 | 63 | 0.760
32 | dense | gaussian | 0.5 | theta[5] | 66.854 | 63 | 0.346
33 | dense | gaussian | 0.5 | theta[6] | 63.027 | 63 | 0.475
34 | dense | gaussian | 0.5 | theta[7] | 51.046 | 63 | 0.860
35 | dense | gaussian | 0.5 | theta[8] | 64.166 | 63 | 0.435
36 | dense | gaussian | 0.5 | theta[9] | 50.829 | 63 | 0.865
39 | dense | gaussian | 1 | theta[1] | 51.917 | 63 | 0.839
40 | dense | gaussian | 1 | theta[10] | 79.616 | 63 | 0.077
41 | dense | gaussian | 1 | theta[2] | 68.928 | 63 | 0.284
42 | dense | gaussian | 1 | theta[3] | 69.274 | 63 | 0.274
43 | dense | gaussian | 1 | theta[4] | 64.243 | 63 | 0.433
44 | dense | gaussian | 1 | theta[5] | 44.442 | 63 | 0.963
45 | dense | gaussian | 1 | theta[6] | 80.691 | 63 | 0.066
46 | dense | gaussian | 1 | theta[7] | 52.326 | 63 | 0.829
47 | dense | gaussian | 1 | theta[8] | 59.123 | 63 | 0.615
48 | dense | gaussian | 1 | theta[9] | 71.437 | 63 | 0.218
51 | dense | poisson | 0.5 | theta[1] | 61.760 | 63 | 0.521
52 | dense | poisson | 0.5 | theta[10] | 66.995 | 63 | 0.342
53 | dense | poisson | 0.5 | theta[2] | 63.642 | 63 | 0.454
54 | dense | poisson | 0.5 | theta[3] | 65.382 | 63 | 0.394
55 | dense | poisson | 0.5 | theta[4] | 72.870 | 63 | 0.185
56 | dense | poisson | 0.5 | theta[5] | 69.760 | 63 | 0.261
57 | dense | poisson | 0.5 | theta[6] | 63.898 | 63 | 0.445
58 | dense | poisson | 0.5 | theta[7] | 61.120 | 63 | 0.544
59 | dense | poisson | 0.5 | theta[8] | 81.203 | 63 | 0.061
60 | dense | poisson | 0.5 | theta[9] | 62.797 | 63 | 0.484
63 | dense | poisson | 1 | theta[1] | 54.938 | 63 | 0.755
64 | dense | poisson | 1 | theta[10] | 57.626 | 63 | 0.668
65 | dense | poisson | 1 | theta[2] | 57.114 | 63 | 0.685
66 | dense | poisson | 1 | theta[3] | 81.946 | 63 | 0.055
67 | dense | poisson | 1 | theta[4] | 68.224 | 63 | 0.304
68 | dense | poisson | 1 | theta[5] | 82.010 | 63 | 0.054
69 | dense | poisson | 1 | theta[6] | 53.222 | 63 | 0.805
70 | dense | poisson | 1 | theta[7] | 87.322 | 63 | 0.023
71 | dense | poisson | 1 | theta[8] | 59.776 | 63 | 0.592
72 | dense | poisson | 1 | theta[9] | 80.909 | 63 | 0.064
75 | sparse | binomial | 0.5 | theta[1] | 73.523 | 63 | 0.171
76 | sparse | binomial | 0.5 | theta[2] | 54.541 | 63 | 0.767
79 | sparse | binomial | 1 | theta[1] | 58.957 | 63 | 0.621
80 | sparse | binomial | 1 | theta[2] | 60.083 | 63 | 0.581
83 | sparse | gaussian | 0.5 | theta[1] | 66.061 | 63 | 0.372
84 | sparse | gaussian | 0.5 | theta[2] | 54.656 | 63 | 0.764
87 | sparse | gaussian | 1 | theta[1] | 51.494 | 63 | 0.850
88 | sparse | gaussian | 1 | theta[2] | 58.010 | 63 | 0.654
91 | sparse | poisson | 0.5 | theta[1] | 71.987 | 63 | 0.205
92 | sparse | poisson | 0.5 | theta[2] | 51.085 | 63 | 0.859
95 | sparse | poisson | 1 | theta[1] | 68.352 | 63 | 0.301
96 | sparse | poisson | 1 | theta[2] | 50.701 | 63 | 0.868

## Session Info

    ## R version 3.6.3 (2020-02-29)
    ## Platform: x86_64-pc-linux-gnu (64-bit)
    ## Running under: Ubuntu 18.04.5 LTS
    ##
    ## Matrix products: default
    ## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
    ## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
    ##
    ## locale:
    ##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
    ##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
    ##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
    ##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
    ##  [9] LC_ADDRESS=C               LC_TELEPHONE=C
    ## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
    ##
    ## attached base packages:
    ## [1] tools     stats     graphics  grDevices utils     datasets  methods
    ## [8] base
    ##
    ## other attached packages:
    ##  [1] purrr_0.3.4        rstan_2.19.3       StanHeaders_2.19.2 RBesT_1.6-3
    ##  [5] testthat_2.3.2     devtools_2.3.0     usethis_1.6.0      ggplot2_3.3.0
    ##  [9] broom_0.5.5        tidyr_1.0.2        dplyr_0.8.5        assertthat_0.2.1
    ## [13] knitr_1.28         rmarkdown_2.1
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] Rcpp_1.0.4.6       mvtnorm_1.1-0      lattice_0.20-44    prettyunits_1.1.1
    ##  [5] ps_1.3.2           rprojroot_1.3-2    digest_0.6.25      plyr_1.8.6
    ##  [9] R6_2.4.1           ggridges_0.5.2     backports_1.1.6    stats4_3.6.3
    ## [13] evaluate_0.14      highr_0.8          pillar_1.4.3       rlang_0.4.5
    ## [17] rstudioapi_0.11    callr_3.4.3        checkmate_2.0.0    desc_1.2.0
    ## [21] stringr_1.4.0      loo_2.2.0          munsell_0.5.0      compiler_3.6.3
    ## [25] xfun_0.13          pkgconfig_2.0.3    pkgbuild_1.0.6     rstantools_2.1.1
    ## [29] htmltools_0.4.0    tidyselect_1.0.0   tibble_3.0.0       gridExtra_2.3
    ## [33] codetools_0.2-18   matrixStats_0.56.0 fansi_0.4.1        crayon_1.3.4
    ## [37] withr_2.1.2        grid_3.6.3         nlme_3.1-152       gtable_0.3.0
    ## [41] lifecycle_0.2.0    magrittr_1.5       scales_1.1.0       RcppParallel_5.1.4
    ## [45] cli_2.0.2          stringi_1.4.6      fs_1.4.1           remotes_2.1.1
    ## [49] ellipsis_0.3.0     generics_0.0.2     vctrs_0.2.4        Formula_1.2-3
    ## [53] glue_1.4.0         processx_3.4.2     pkgload_1.0.2      parallel_3.6.3
    ## [57] yaml_2.2.1         inline_0.3.15      colorspace_1.4-1   sessioninfo_1.1.1
    ## [61] bayesplot_1.7.1    memoise_1.1.0