# Simulation based calibration for RBesT

#### Sebastian Weber

#### Thu May 28 08:59:30 2020

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:  2020-05-27 11:26:59 UTC
    ## git hash: f5a086c54d19b4972e3896dbb8efa2022d7d8015
    ## MD5:      1221d359f47b231f1483fa4274d42113

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

    ##                    calibration.rds
    ## "1221d359f47b231f1483fa4274d42113"

# 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[10] | 62.925 | 63 | 0.479
4 | dense | binomial | 0.5 | theta[1] | 69.709 | 63 | 0.262
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[10] | 79.834 | 63 | 0.075
16 | dense | binomial | 1 | theta[1] | 59.648 | 63 | 0.597
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[10] | 66.317 | 63 | 0.363
28 | dense | gaussian | 0.5 | theta[1] | 54.208 | 63 | 0.777
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[10] | 79.616 | 63 | 0.077
40 | dense | gaussian | 1 | theta[1] | 51.917 | 63 | 0.839
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[10] | 66.995 | 63 | 0.342
52 | dense | poisson | 0.5 | theta[1] | 61.760 | 63 | 0.521
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[10] | 57.626 | 63 | 0.668
64 | dense | poisson | 1 | theta[1] | 54.938 | 63 | 0.755
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.1 (2019-07-05)
    ## Platform: x86_64-pc-linux-gnu (64-bit)
    ## Running under: Ubuntu 16.04.6 LTS
    ##
    ## Matrix products: default
    ## BLAS:   /usr/lib/libblas/libblas.so.3.6.0
    ## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
    ##
    ## locale:
    ## [1] C
    ##
    ## attached base packages:
    ## [1] tools     stats     graphics  grDevices utils     datasets  methods
    ## [8] base
    ##
    ## other attached packages:
    ##  [1] purrr_0.3.3        rstan_2.19.2       StanHeaders_2.19.0
    ##  [4] RBesT_1.6-1        testthat_2.2.1     Rcpp_1.0.2
    ##  [7] devtools_2.2.1     usethis_1.5.1      ggplot2_3.2.1
    ## [10] broom_0.5.2        tidyr_1.0.0        dplyr_0.8.3
    ## [13] assertthat_0.2.1   knitr_1.25         rmarkdown_1.16
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] mvtnorm_1.0-11     lattice_0.20-38    prettyunits_1.0.2
    ##  [4] ps_1.3.0           zeallot_0.1.0      rprojroot_1.3-2
    ##  [7] digest_0.6.21      plyr_1.8.4         R6_2.4.0
    ## [10] ggridges_0.5.1     backports_1.1.5    stats4_3.6.1
    ## [13] evaluate_0.14      highr_0.8          pillar_1.4.2
    ## [16] rlang_0.4.0        lazyeval_0.2.2     rstudioapi_0.10
    ## [19] callr_3.3.2        checkmate_1.9.4    desc_1.2.0
    ## [22] stringr_1.4.0      loo_2.1.0          munsell_0.5.0
    ## [25] compiler_3.6.1     xfun_0.10          pkgconfig_2.0.3
    ## [28] pkgbuild_1.0.6     htmltools_0.4.0    tidyselect_0.2.5
    ## [31] tibble_2.1.3       gridExtra_2.3      codetools_0.2-16
    ## [34] matrixStats_0.55.0 crayon_1.3.4       withr_2.1.2
    ## [37] grid_3.6.1         nlme_3.1-141       gtable_0.3.0
    ## [40] lifecycle_0.1.0    magrittr_1.5       scales_1.0.0
    ## [43] cli_1.1.0          stringi_1.4.3      fs_1.3.1
    ## [46] remotes_2.1.0      ellipsis_0.3.0     generics_0.0.2
    ## [49] vctrs_0.2.0        Formula_1.2-3      glue_1.3.1
    ## [52] processx_3.4.1     pkgload_1.0.2      parallel_3.6.1
    ## [55] yaml_2.2.0         inline_0.3.15      colorspace_1.4-1
    ## [58] sessioninfo_1.1.1  bayesplot_1.7.0    memoise_1.1.0