Additional power Figure under varying pc
The following script demonstrates the RBesT library to reproduce the main results from Schmidli et al., Biometrics 70, 1024, 2014.
The two main ideas of the paper are
use mixture priors to approximate accuratley numerical MCMC MAP priors
robustify informative MAP priors by adding a suitable non-informative component to the informative MAP prior
As example an adaptive design for a binomial endpoint is considered:
Stage 1: mI in test treatment and nI in control (e.g., mI = 20, nI = 15);
Stage 2: (m - mI ) in test treatment and max(n - ESSI, nmin) in control (e.g., nmin = 5).
| pc | 0.3_beta | 0.3_mix50 | 0.3_mix90 | 0.3_unif | 0_beta | 0_mix50 | 0_mix90 | 0_unif |
|---|---|---|---|---|---|---|---|---|
| 0.1 | 81.5 | 83.4 | 82.1 | 88.9 | 0.1 | 0.3 | 0.1 | 1.8 |
| 0.2 | 87.2 | 81.7 | 85.2 | 79.8 | 1.7 | 1.9 | 1.7 | 2.4 |
| 0.3 | 93.2 | 79.6 | 86.6 | 77.3 | 6.1 | 4.1 | 5.6 | 2.3 |
| 0.4 | 97.9 | 78.3 | 84.0 | 78.6 | 13.4 | 4.9 | 9.6 | 2.3 |
| 0.5 | 99.7 | 81.1 | 81.1 | 81.2 | 26.1 | 4.1 | 10.5 | 2.7 |
| 0.6 | 100.0 | 88.8 | 86.7 | 88.4 | 44.9 | 3.1 | 6.7 | 2.7 |
| pc | 0.3_beta | 0.3_mix50 | 0.3_mix90 | 0.3_unif | 0_beta | 0_mix50 | 0_mix90 | 0_unif |
|---|---|---|---|---|---|---|---|---|
| 0.1 | 20 | 25.1 | 20.7 | 38 | 20 | 25.1 | 20.7 | 38 |
| 0.2 | 20 | 25.5 | 20.9 | 38 | 20 | 25.5 | 20.9 | 38 |
| 0.3 | 20 | 28.9 | 21.7 | 38 | 20 | 28.9 | 21.7 | 38 |
| 0.4 | 20 | 33.7 | 23.8 | 38 | 20 | 33.7 | 23.8 | 38 |
| 0.5 | 20 | 37.5 | 27.5 | 38 | 20 | 37.5 | 27.5 | 38 |
| 0.6 | 20 | 38.9 | 32.3 | 38 | 20 | 38.9 | 32.3 | 38 |
Reproduction of Fig. 1 in Robust MAP Prior paper.
The bias and rMSE calculations are slightly involved as the sample size depends on the first stage.
Clinical example to exemplify the methodology.
## set seed to guarantee exact reproducible results
set.seed(25445)
map <- gMAP(cbind(r, n-r) ~ 1 | study,
family=binomial,
data=colitis,
tau.dist="HalfNormal",
beta.prior=2,
tau.prior=1)## Assuming default prior location for beta: 0map_auto <- automixfit(map)
## advanced: look at AIC of all fitted models
sapply(attr(map_auto, "models"), AIC)## 4 2 3 1
## -11069.796 -10999.555 -11007.529 -9905.178## EM for Beta Mixture Model
## Log-Likelihood = 5545.898
##
## Univariate beta mixture
## Mixture Components:
## comp1 comp2 comp3 comp4
## w 0.43136888 0.39133817 0.09149862 0.08579434
## a 13.65638421 3.44347219 7.78461437 0.91188780
## b 105.20424486 32.07393852 25.44873451 2.78278722## use best fitting mixture model as prior
prior <- map_auto
pl <- plot(prior)
pl$mix + ggtitle("MAP prior for ulcerative colitis")
Colitis MAPs from paper for further figures.
mapCol <- list(
one = mixbeta(c(1,2.3,16)),
two = mixbeta(c(0.77, 6.2, 50.8), c(1-0.77, 1.0, 4.7)),
three = mixbeta(c(0.53, 2.5, 19.1), c(0.38, 14.6, 120.2), c(0.08, 0.9, 2.9))
)## Warning in mixdist3(...): Weights do not sum to 1. Rescaling accordingly.mapCol <- c(mapCol, list(twoRob=robustify(mapCol$two, weight=0.1, mean=1/2),
threeRob=robustify(mapCol$three, weight=0.1, mean=1/2)
)
)Posterior for different remission rates, Figure 3
N <- 20
post <- foreach(prior=names(mapCol), .combine=rbind) %do% {
res <- data.frame(mean=rep(NA, N+1), sd=0, r=0:N)
for(r in 0:N) {
res[r+1,1:2] <- summary(postmix(mapCol[[prior]], r=r, n=N))[c("mean", "sd")]
}
res$prior <- prior
res
}
qplot(r, mean, data=post, colour=prior, shape=prior) + geom_abline(slope=1/20)
## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] C/UTF-8/C/C/C/C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] foreach_1.4.7 forcats_0.4.0 stringr_1.4.0 readr_1.3.1
## [5] tidyr_1.0.0 tibble_2.1.3 tidyverse_1.2.1 scales_1.0.0
## [9] ggplot2_3.2.1 purrr_0.3.3 dplyr_0.8.3 bayesplot_1.7.0
## [13] knitr_1.25 RBesT_1.5-4 Rcpp_1.0.2
##
## loaded via a namespace (and not attached):
## [1] httr_1.4.1 jsonlite_1.6 modelr_0.1.5
## [4] StanHeaders_2.19.0 Formula_1.2-3 assertthat_0.2.1
## [7] highr_0.8 stats4_3.5.2 cellranger_1.1.0
## [10] yaml_2.2.0 pillar_1.4.2 backports_1.1.5
## [13] lattice_0.20-38 glue_1.3.1 digest_0.6.21
## [16] checkmate_1.9.4 rvest_0.3.4 colorspace_1.4-1
## [19] htmltools_0.4.0 plyr_1.8.4 pkgconfig_2.0.3
## [22] rstan_2.19.2 broom_0.5.2 haven_2.1.1
## [25] mvtnorm_1.0-11 processx_3.4.1 generics_0.0.2
## [28] ellipsis_0.3.0 withr_2.1.2 lazyeval_0.2.2
## [31] cli_1.1.0 magrittr_1.5 crayon_1.3.4
## [34] readxl_1.3.1 evaluate_0.14 ps_1.3.0
## [37] nlme_3.1-137 xml2_1.2.2 pkgbuild_1.0.6
## [40] tools_3.5.2 loo_2.1.0 prettyunits_1.0.2
## [43] hms_0.5.1 lifecycle_0.1.0 matrixStats_0.55.0
## [46] munsell_0.5.0 callr_3.3.2 compiler_3.5.2
## [49] rlang_0.4.0 grid_3.5.2 ggridges_0.5.1
## [52] iterators_1.0.12 rstudioapi_0.10 labeling_0.3
## [55] rmarkdown_1.16 gtable_0.3.0 codetools_0.2-15
## [58] inline_0.3.15 reshape2_1.4.3 R6_2.4.0
## [61] gridExtra_2.3 lubridate_1.7.4 zeallot_0.1.0
## [64] stringi_1.4.3 parallel_3.5.2 vctrs_0.2.0
## [67] tidyselect_0.2.5 xfun_0.10