Simulation based calibration for OncoBayes2

Fri May 7 18:58:59 2021

This report documents the results of a simulation based calibration (SBC) run for OncoBayes2. TODO

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

## Created:  2021-05-07 09:46:51 UTC
## git hash: 20b6ac288f3269d4c5137140e97e454552bd75d0
## MD5:      d2152012b4939338ed237dc96e2df7c2

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

##                    calibration.rds 
## "d2152012b4939338ed237dc96e2df7c2"

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}.\]

Model description TODO

The fake data simulation function returns … TODO. 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

model problem N total_divergent min_ess max_Rhat total_large_Rhat min_lp_ess_bulk min_lp_ess_tail
combo2_EX base 9950 0 362 1.012 0 67 159
combo2_EX warmup_base 50 0 739 1.004 0 247 256
combo2_EXNEX base 9950 0 22 1.058 0 100 237
combo2_EXNEX warmup_base 50 0 146 1.011 0 226 245
combo3_EXNEX base 9950 0 39 1.048 0 140 177
combo3_EXNEX warmup_base 50 0 250 1.024 0 258 403
log2bayes_EXNEX base 9950 0 4 2.362 1 33 233
log2bayes_EXNEX warmup_base 50 0 54 1.027 0 280 391

Large Rhat is defined as exceeding \(1.2\).

\(\chi^2\) Statistic, Model 1: Single-agent logistic regression

param statistic df p.value
beta_group[A,I(log(drug_A/1)),intercept] 30.726 31 0.480
beta_group[A,I(log(drug_A/1)),log_slope] 34.733 31 0.295
beta_group[B,I(log(drug_A/1)),intercept] 24.205 31 0.802
beta_group[B,I(log(drug_A/1)),log_slope] 39.450 31 0.142
beta_group[C,I(log(drug_A/1)),intercept] 21.837 31 0.888
beta_group[C,I(log(drug_A/1)),log_slope] 30.310 31 0.501
mu_log_beta[I(log(drug_A/1)),intercept] 21.914 31 0.886
mu_log_beta[I(log(drug_A/1)),log_slope] 32.781 31 0.380
tau_log_beta[STRAT,I(log(drug_A/1)),intercept] 32.486 31 0.393
tau_log_beta[STRAT,I(log(drug_A/1)),log_slope] 42.918 31 0.075

\(\chi^2\) Statistic, Model 2: Double combination, fully exchangeable

param statistic df p.value
beta_group[A,I(log(drug_A/1)),intercept] 26.554 31 0.694
beta_group[A,I(log(drug_A/1)),log_slope] 30.336 31 0.500
beta_group[A,I(log(drug_B/1)),intercept] 27.507 31 0.646
beta_group[A,I(log(drug_B/1)),log_slope] 36.173 31 0.240
beta_group[B,I(log(drug_A/1)),intercept] 29.818 31 0.527
beta_group[B,I(log(drug_A/1)),log_slope] 23.411 31 0.834
beta_group[B,I(log(drug_B/1)),intercept] 44.211 31 0.058
beta_group[B,I(log(drug_B/1)),log_slope] 39.462 31 0.142
beta_group[C,I(log(drug_A/1)),intercept] 36.838 31 0.217
beta_group[C,I(log(drug_A/1)),log_slope] 37.517 31 0.195
beta_group[C,I(log(drug_B/1)),intercept] 37.069 31 0.209
beta_group[C,I(log(drug_B/1)),log_slope] 34.426 31 0.307
eta_group[A,I(drug_A/1 * drug_B/1)] 30.323 31 0.501
eta_group[B,I(drug_A/1 * drug_B/1)] 47.136 31 0.032
eta_group[C,I(drug_A/1 * drug_B/1)] 41.933 31 0.091
mu_eta[I(drug_A/1 * drug_B/1)] 28.998 31 0.569
mu_log_beta[I(log(drug_A/1)),intercept] 39.053 31 0.152
mu_log_beta[I(log(drug_A/1)),log_slope] 24.448 31 0.792
mu_log_beta[I(log(drug_B/1)),intercept] 40.915 31 0.110
mu_log_beta[I(log(drug_B/1)),log_slope] 30.803 31 0.476
tau_eta[STRAT,I(drug_A/1 * drug_B/1)] 34.784 31 0.292
tau_log_beta[STRAT,I(log(drug_A/1)),intercept] 34.643 31 0.298
tau_log_beta[STRAT,I(log(drug_A/1)),log_slope] 24.390 31 0.794
tau_log_beta[STRAT,I(log(drug_B/1)),intercept] 20.653 31 0.921
tau_log_beta[STRAT,I(log(drug_B/1)),log_slope] 24.314 31 0.798

\(\chi^2\) Statistic, Model 3: Double combination, EXchangeable/NonEXchangeable model

param statistic df p.value
beta_group[A,I(log(drug_A/1)),intercept] 27.680 31 0.638
beta_group[A,I(log(drug_A/1)),log_slope] 28.947 31 0.572
beta_group[A,I(log(drug_B/1)),intercept] 29.498 31 0.543
beta_group[A,I(log(drug_B/1)),log_slope] 30.586 31 0.487
beta_group[B,I(log(drug_A/1)),intercept] 46.650 31 0.035
beta_group[B,I(log(drug_A/1)),log_slope] 36.410 31 0.231
beta_group[B,I(log(drug_B/1)),intercept] 31.046 31 0.464
beta_group[B,I(log(drug_B/1)),log_slope] 40.256 31 0.123
beta_group[C,I(log(drug_A/1)),intercept] 35.629 31 0.260
beta_group[C,I(log(drug_A/1)),log_slope] 28.576 31 0.591
beta_group[C,I(log(drug_B/1)),intercept] 23.987 31 0.811
beta_group[C,I(log(drug_B/1)),log_slope] 36.141 31 0.241
eta_group[A,I(drug_A/1 * drug_B/1)] 31.757 31 0.429
eta_group[B,I(drug_A/1 * drug_B/1)] 32.915 31 0.373
eta_group[C,I(drug_A/1 * drug_B/1)] 34.701 31 0.296
mu_eta[I(drug_A/1 * drug_B/1)] 32.762 31 0.381
mu_log_beta[I(log(drug_A/1)),intercept] 37.466 31 0.197
mu_log_beta[I(log(drug_A/1)),log_slope] 30.150 31 0.510
mu_log_beta[I(log(drug_B/1)),intercept] 46.938 31 0.033
mu_log_beta[I(log(drug_B/1)),log_slope] 44.506 31 0.055
tau_eta[STRAT,I(drug_A/1 * drug_B/1)] 44.858 31 0.051
tau_log_beta[STRAT,I(log(drug_A/1)),intercept] 32.237 31 0.405
tau_log_beta[STRAT,I(log(drug_A/1)),log_slope] 28.909 31 0.574
tau_log_beta[STRAT,I(log(drug_B/1)),intercept] 26.406 31 0.702
tau_log_beta[STRAT,I(log(drug_B/1)),log_slope] 18.502 31 0.963

\(\chi^2\) Statistic, Model 4: Triple combination, EX/NEX model

param statistic df p.value
beta_group[A,I(log(drug_A/1)),intercept] 31.718 31 0.430
beta_group[A,I(log(drug_A/1)),log_slope] 24.653 31 0.783
beta_group[A,I(log(drug_B/1)),intercept] 25.587 31 0.741
beta_group[A,I(log(drug_B/1)),log_slope] 28.813 31 0.579
beta_group[A,I(log(drug_C/1)),intercept] 24.026 31 0.809
beta_group[A,I(log(drug_C/1)),log_slope] 21.894 31 0.886
beta_group[B,I(log(drug_A/1)),intercept] 29.402 31 0.548
beta_group[B,I(log(drug_A/1)),log_slope] 26.016 31 0.721
beta_group[B,I(log(drug_B/1)),intercept] 19.840 31 0.939
beta_group[B,I(log(drug_B/1)),log_slope] 37.350 31 0.200
beta_group[B,I(log(drug_C/1)),intercept] 55.750 31 0.004
beta_group[B,I(log(drug_C/1)),log_slope] 28.090 31 0.617
beta_group[C,I(log(drug_A/1)),intercept] 42.208 31 0.086
beta_group[C,I(log(drug_A/1)),log_slope] 28.710 31 0.584
beta_group[C,I(log(drug_B/1)),intercept] 27.578 31 0.643
beta_group[C,I(log(drug_B/1)),log_slope] 32.256 31 0.404
beta_group[C,I(log(drug_C/1)),intercept] 43.763 31 0.064
beta_group[C,I(log(drug_C/1)),log_slope] 31.162 31 0.458
eta_group[A,I(drug_A/1 * drug_B/1 * drug_C/1)] 30.643 31 0.484
eta_group[A,I(drug_A/1 * drug_B/1)] 35.334 31 0.271
eta_group[A,I(drug_A/1 * drug_C/1)] 40.979 31 0.108
eta_group[A,I(drug_B/1 * drug_C/1)] 35.968 31 0.247
eta_group[B,I(drug_A/1 * drug_B/1 * drug_C/1)] 36.218 31 0.238
eta_group[B,I(drug_A/1 * drug_B/1)] 30.816 31 0.476
eta_group[B,I(drug_A/1 * drug_C/1)] 40.653 31 0.115
eta_group[B,I(drug_B/1 * drug_C/1)] 28.666 31 0.587
eta_group[C,I(drug_A/1 * drug_B/1 * drug_C/1)] 29.926 31 0.521
eta_group[C,I(drug_A/1 * drug_B/1)] 48.019 31 0.026
eta_group[C,I(drug_A/1 * drug_C/1)] 22.592 31 0.864
eta_group[C,I(drug_B/1 * drug_C/1)] 39.027 31 0.152
mu_eta[I(drug_A/1 * drug_B/1 * drug_C/1)] 33.990 31 0.326
mu_eta[I(drug_A/1 * drug_B/1)] 40.090 31 0.127
mu_eta[I(drug_A/1 * drug_C/1)] 35.846 31 0.251
mu_eta[I(drug_B/1 * drug_C/1)] 42.989 31 0.074
mu_log_beta[I(log(drug_A/1)),intercept] 27.942 31 0.624
mu_log_beta[I(log(drug_A/1)),log_slope] 25.082 31 0.764
mu_log_beta[I(log(drug_B/1)),intercept] 37.408 31 0.198
mu_log_beta[I(log(drug_B/1)),log_slope] 39.002 31 0.153
mu_log_beta[I(log(drug_C/1)),intercept] 32.013 31 0.416
mu_log_beta[I(log(drug_C/1)),log_slope] 18.195 31 0.967
tau_eta[STRAT,I(drug_A/1 * drug_B/1 * drug_C/1)] 28.358 31 0.603
tau_eta[STRAT,I(drug_A/1 * drug_B/1)] 19.437 31 0.947
tau_eta[STRAT,I(drug_A/1 * drug_C/1)] 34.694 31 0.296
tau_eta[STRAT,I(drug_B/1 * drug_C/1)] 28.998 31 0.569
tau_log_beta[STRAT,I(log(drug_A/1)),intercept] 15.949 31 0.988
tau_log_beta[STRAT,I(log(drug_A/1)),log_slope] 25.453 31 0.747
tau_log_beta[STRAT,I(log(drug_B/1)),intercept] 29.843 31 0.525
tau_log_beta[STRAT,I(log(drug_B/1)),log_slope] 38.758 31 0.159
tau_log_beta[STRAT,I(log(drug_C/1)),intercept] 18.944 31 0.956
tau_log_beta[STRAT,I(log(drug_C/1)),log_slope] 36.691 31 0.222

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] mvtnorm_1.0-11     tibble_2.1.3       rstan_2.19.3      
##  [4] StanHeaders_2.19.2 abind_1.4-5        Formula_1.2-3     
##  [7] checkmate_1.9.4    OncoBayes2_0.7-0   testthat_2.2.1    
## [10] Rcpp_1.0.2         devtools_2.2.1     usethis_1.5.1     
## [13] ggplot2_3.2.1      broom_0.5.2        tidyr_1.0.0       
## [16] dplyr_0.8.3        assertthat_0.2.1   knitr_1.25        
## [19] rmarkdown_1.16    
## 
## loaded via a namespace (and not attached):
##  [1] lattice_0.20-38    prettyunits_1.0.2  ps_1.3.0          
##  [4] zeallot_0.1.0      rprojroot_1.3-2    digest_0.6.21     
##  [7] plyr_1.8.4         R6_2.4.0           ggridges_0.5.1    
## [10] backports_1.1.5    stats4_3.6.1       evaluate_0.14     
## [13] highr_0.8          pillar_1.4.2       rlang_0.4.0       
## [16] lazyeval_0.2.2     rstudioapi_0.10    callr_3.3.2       
## [19] desc_1.2.0         stringr_1.4.0      loo_2.1.0         
## [22] munsell_0.5.0      compiler_3.6.1     xfun_0.10         
## [25] pkgconfig_2.0.3    pkgbuild_1.0.6     rstantools_2.0.0  
## [28] htmltools_0.4.0    tidyselect_0.2.5   gridExtra_2.3     
## [31] codetools_0.2-16   matrixStats_0.55.0 crayon_1.3.4      
## [34] withr_2.1.2        grid_3.6.1         nlme_3.1-141      
## [37] gtable_0.3.0       lifecycle_0.1.0    magrittr_1.5      
## [40] scales_1.0.0       cli_1.1.0          stringi_1.4.3     
## [43] fs_1.3.1           remotes_2.1.0      ellipsis_0.3.0    
## [46] generics_0.0.2     vctrs_0.2.0        glue_1.3.1        
## [49] purrr_0.3.3        processx_3.4.1     pkgload_1.0.2     
## [52] parallel_3.6.1     yaml_2.2.0         inline_0.3.15     
## [55] colorspace_1.4-1   sessioninfo_1.1.1  bayesplot_1.7.0   
## [58] memoise_1.1.0