Using RBesT to reproduce Schmidli et. al “Robust MAP Priors”

Sebastian Weber

2017-07-12

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

  1. use mixture priors to approximate accuratley numerical MCMC MAP priors

  2. 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).

Operating Characteristics, Table 1

Type I error and power
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
Sample size
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

Additional power Figure under varying pc

Bias and rMSE, Figure 1

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.

Ulcerative colitis example

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: 0
map_auto <- automixfit(map)

## advanced: look at AIC of all fitted models
sapply(attr(map_auto, "models"), AIC)
##          4          3          2          1 
## -10976.675 -10962.891 -10877.326  -9672.574
print(map_auto)
## EM for Beta Mixture Model
## Log-Likelihood = 5499.337
## 
## Univariate beta mixture
## Mixture Components:
##   comp1       comp2       comp3       comp4      
## w  0.42795398  0.41885398  0.12079611  0.03239593
## a 12.66068764  2.63622438  2.45105356  0.68401854
## b 98.73362970 21.96275310  8.90036038  1.44030153
## 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)

qplot(r, sd, data=post, colour=prior, shape=prior) + coord_cartesian(ylim=c(0,0.17))

sessionInfo()
## R version 3.3.3 (2017-03-06)
## Platform: i386-w64-mingw32/i386 (32-bit)
## Running under: Windows 7 (build 7601) Service Pack 1
## 
## locale:
## [1] LC_COLLATE=C                        LC_CTYPE=German_Switzerland.1252   
## [3] LC_MONETARY=German_Switzerland.1252 LC_NUMERIC=C                       
## [5] LC_TIME=German_Switzerland.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] foreach_1.4.3   bayesplot_1.2.0 RBesT_1.2-2     Rcpp_0.12.11   
##  [5] knitr_1.16      dplyr_0.7.1     purrr_0.2.2.2   readr_1.1.1    
##  [9] tidyr_0.6.3     tibble_1.3.3    ggplot2_2.2.1   tidyverse_1.1.1
## 
## loaded via a namespace (and not attached):
##  [1] checkmate_1.8.3      jsonlite_1.5         reshape2_1.4.2      
##  [4] magrittr_1.5         gtable_0.2.0         rmarkdown_1.6       
##  [7] hms_0.3              xml2_1.1.1           htmltools_0.3.6     
## [10] forcats_0.2.0        stringr_1.2.0        haven_1.1.0         
## [13] broom_0.4.2          cellranger_1.1.0     Formula_1.2-1       
## [16] lattice_0.20-34      StanHeaders_2.16.0-1 plyr_1.8.4          
## [19] lubridate_1.6.0      gridExtra_2.2.1      iterators_1.0.8     
## [22] pkgconfig_2.0.1      R6_2.2.2             digest_0.6.12       
## [25] stats4_3.3.3         colorspace_1.3-2     bindrcpp_0.2        
## [28] rprojroot_1.2        rstan_2.16.2         stringi_1.1.5       
## [31] yaml_2.1.14          lazyeval_0.2.0       codetools_0.2-15    
## [34] evaluate_0.10.1      labeling_0.3         compiler_3.3.3      
## [37] httr_1.2.1           bindr_0.1            backports_1.1.0     
## [40] inline_0.3.14        munsell_0.4.3        psych_1.7.5         
## [43] modelr_0.1.0         readxl_1.0.0         highr_0.6           
## [46] parallel_3.3.3       assertthat_0.2.0     tools_3.3.3         
## [49] foreign_0.8-67       mvtnorm_1.0-6        scales_0.4.1        
## [52] glue_1.1.1           rlang_0.1.1          mnormt_1.5-5        
## [55] nlme_3.1-131         rvest_0.3.2          grid_3.3.3