demo/robustMAP.R

In RBesT: R Bayesian Evidence Synthesis Tools

#' ---
#' title: "Using RBesT to reproduce Schmidli et al. \"Robust MAP Priors\""
#' author: "Sebastian Weber"
#' date: "`r Sys.Date()`"
#' output: rmarkdown::html_vignette
#' vignette: >
#'   %\VignetteIndexEntry{Using RBesT to reproduce Schmidli et al. "Robust MAP Priors"}
#'   %\VignetteEngine{knitr::rmarkdown}
#'   %\VignetteDepends{foreach}
#'   %\VignetteDepends{reshape2}
#' ---
#'
#' 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).
#'
#' To render a report with all the figures, please run:
#'
#' - ` rmarkdown::render("robustMAP.R")`
#'
#' - ` rmarkdown::render("robustMAP.R", "pdf_document")` (for a pdf version)
#'

library(foreach)
library(ggplot2)
library(dplyr)
library(tidyr)
library(RBesT)
library(knitr)
theme_set(theme_bw())
knitr::opts_chunk$set(
  fig.width = 7,
  fig.height = 4
)

#'
#' ## Operating Characteristics, Table 1
#'

#' the different priors
map <- list()
map$beta <- mixbeta(c(1.0, 4, 16))
map$mix90 <- mixbeta(c(0.9, 4, 16), c(0.1, 1, 1))
## map$mix70 <- mixbeta(c(0.7, 4, 16), c(0.3, 1, 1))
map$mix50 <- mixbeta(c(0.5, 4, 16), c(0.5, 1, 1))
## map$mix30 <- mixbeta(c(0.3, 4, 16), c(0.7, 1, 1))
## map$mix10 <- mixbeta(c(0.1, 4, 16), c(0.9, 1, 1))
map$unif <- mixbeta(c(1.0, 1, 1))

unif <- map$unif

## for the adaptive design the calculation is a bit involved as we
## have to calculate all possible ctl sample sizes which is determined
## by the ESS at the intermediate

OC_adaptBinary2S <- function(N1, Ntarget, Nmin, M, ctl.prior, treat.prior, pc, pt, decision) {
  ## calculate the different possible ESS values which we get after
  ## stage1
  r1 <- 0:N1
  ESSstage1 <- vector("double", N1 + 1)
  for (r in r1) {
    ESSstage1[r + 1] <- round(ess(postmix(ctl.prior, r = r, n = N1), method = "morita"))
  }

  ## number of patients enrolled in stage 2
  N2 <- pmax(Ntarget - ESSstage1, Nmin)

  ## total number of patients enrolled
  N <- N1 + N2

  P <- try(data.frame(pc = pc, pt = pt))
  if (inherits(P, "try-error")) {
    stop("pc and pt need to be of same size")
  }

  ## calculate for each scenario and sample size of the control the
  ## power
  power_all <- matrix(0, N1 + 1, nrow(P))
  for (r in r1) {
    ## power_all[r+1,] <- OC_binary2S(N[r+1], M, ctl.prior, treat.prior, P$pc, P$pt, crit)
    design_calc <- oc2S(treat.prior, ctl.prior, M, N[r + 1], decision)
    power_all[r + 1, ] <- design_calc(P$pt, P$pc)
  }

  ## finally take the mean with the respective weight which corresponds
  ## to the weight how the respective sample size occur
  w <- sapply(P$pc, function(p) dbinom(r1, N1, p))

  data.frame(power = colSums(power_all * w), samp = colSums(w * N))
}


## minimum of patients enrolled in stage 2
Nmin <- 5
## number of patients enrolled to control in stage 1
N1 <- 15
## target number of patients overall
Ntarget <- 40
## target number of patients in treatment group
M <- 40

## decision function: P(x1 - x2 > 0) > 0.975
dec <- decision2S(0.975, 0, lower.tail = FALSE)

cases <- expand.grid(prior = names(map), pc = seq(0.1, 0.6, by = 0.1), delta = c(0, 0.3))

## the mixture cases have a varying ess at the interim and need the
## adaptive function ...
cases.mix <- grep("mix", names(map), value = TRUE)
cases.fix <- grep("mix", names(map), value = TRUE, invert = TRUE)

resMix <- foreach(i = cases.mix, .combine = rbind) %do% {
  design <- subset(cases, prior == i)
  cbind(design, OC_adaptBinary2S(N1, Ntarget, Nmin, M, map[[i]], unif, design$pc, design$pc + design$delta, dec))
}

## ... for the non-mixture priors the ESS is fixed at the intermediate
## step such that the much faster oc2S can be used directly
resFix <- foreach(i = cases.fix, .combine = rbind) %do% {
  design <- subset(cases, prior == i)
  prior <- map[[i]]
  Nc <- Ntarget - round(ess(prior, method = "morita"))
  design_calc <- oc2S(unif, prior, M, Nc, dec)
  cbind(design, power = design_calc(design$pc + design$delta, design$pc), samp = Nc)
}

powerTable <- rbind(resMix, resFix)

P <- expand.grid(pc = c(0.2, 0.3, 0.4, 0.5), pt = seq(0.05, 0.95, by = 0.025))

powerMix <- foreach(i = cases.mix, .combine = rbind) %do% {
  cbind(P, prior = i, OC_adaptBinary2S(N1, Ntarget, Nmin, M, map[[i]], unif, P$pc, P$pt, dec))
}

powerFix <- foreach(i = cases.fix, .combine = rbind) %do% {
  prior <- map[[i]]
  Nc <- Ntarget - round(ess(prior, method = "morita"))
  design_calc <- oc2S(unif, prior, M, Nc, dec)
  cbind(P, prior = i, power = design_calc(P$pt, P$pc), samp = Nc)
  ## cbind(P, prior=i, power=OC_binary2S(Nc, M, prior, unif, P$pc, P$pt, dec), samp=Nc)
}

power <- rbind(powerMix, powerFix)


ocAdapt <- powerTable[, -ncol(powerTable)] %>%
  unite(case, delta, prior) %>%
  transform(power = 100 * power) %>%
  spread(case, power)
ocAdaptSamp <- powerTable[, -(ncol(powerTable) - 1)] %>%
  unite(case, delta, prior) %>%
  spread(case, samp)

kable(ocAdapt, digits = 1, caption = "Type I error and power")

kable(ocAdaptSamp, digits = 1, caption = "Sample size")

#'
#' ## Additional power Figure under varying pc
#'

ggplot(power, aes(pt - pc, power, colour = prior)) +
  geom_line() +
  facet_wrap(~pc) +
  scale_y_continuous(breaks = seq(0, 1, by = 0.2)) +
  scale_x_continuous(breaks = seq(-0.8, 0.8, by = 0.2)) +
  coord_cartesian(xlim = c(-0.15, 0.5)) +
  geom_hline(yintercept = 0.025, linetype = 2) +
  geom_hline(yintercept = 0.8, linetype = 2) +
  ggtitle("Prob. for alternative for different pc")

#'
#' ## Bias and rMSE, Figure 1
#'
#' Reproduction of Fig. 1 in Robust MAP Prior paper.
#'


plot(map$beta, prob = 1)
plot(map$mix50)

#' The bias and rMSE calculations are slightly involved as the sample
#' size depends on the first stage.

est <- foreach(case = names(map), .combine = rbind) %do% {
  ## prior to consider
  prior <- map[[case]]

  ## calculate the different possible ESS values which we get after
  ## stage1
  r1 <- 0:N1
  ESSstage1 <- c()
  for (r in r1) {
    ESSstage1 <- c(ESSstage1, round(ess(postmix(prior, r = r, n = N1), method = "morita", loc = "mode")))
  }

  ## number of patients enrolled in stage 2
  N2 <- pmax(Ntarget - ESSstage1, 5)

  ## total number of patients enrolled
  N <- N1 + N2

  ## we need the maximal possible number of patients
  Nmax <- max(N)

  ## calculate for each i = 0 to N1 possible responders in stage one
  ## the posterior when observing 0 to N[i] in total. Calculate for
  ## each scenario outcome E(p) and E(p^2)
  m <- matrix(0, N1 + 1, Nmax + 1)
  m2 <- matrix(0, N1 + 1, Nmax + 1)
  for (i in seq_along(N)) {
    n <- N[i]
    for (r in 0:n) {
      res <- summary(postmix(prior, r = r, n = n))[c("mean", "sd")]
      m[i, r + 1] <- res["mean"]
      m2[i, r + 1] <- res["sd"]^2 + m[i, r + 1]^2
    }
  }

  ## now collect the terms correctly weighted for each assumed true rate
  bias <- rMSE <- c()
  pt <- seq(0, 1, length = 101)
  for (p in pt) {
    ## weight for each possible N at stage 1
    wnp <- dbinom(0:N1, N1, p)

    ## E(p) and E(p^2) for each possible N at stage 1
    Mnm <- rep(0, N1 + 1)
    Mnm2 <- rep(0, N1 + 1)

    ## for a given weight at stage 1....
    for (i in seq(N1 + 1)) {
      n <- N[i]
      ## weights of possible outcomes when having n draws in
      ## stage1, we go up to Nmax+1 to get a vector of correct
      ## length; all entries above n are set to 0 from dbinom as
      ## expected as we can never observe more counts than the
      ## number of trials...
      wp <- dbinom(0:Nmax, n, p)

      Mnm[i] <- sum(m[i, ] * wp)
      Mnm2[i] <- sum(m2[i, ] * wp)
    }

    ## ... which we average over possible outcomes in stage 1
    Mm <- sum(wnp * Mnm)
    Mm2 <- sum(wnp * Mnm2)

    bias <- c(bias, (Mm - p))
    rMSE <- c(rMSE, sqrt(Mm2 - 2 * p * Mm + p^2))
  }
  data.frame(p = pt, bias = bias, rMSE = rMSE, prior = case)
}


ggplot(est, aes(p, 100 * bias, colour = prior)) +
  geom_line() +
  ggtitle("Bias")
ggplot(est, aes(p, 100 * rMSE, colour = prior)) +
  geom_line() +
  ggtitle("rMSE")


#'
#' ##  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
)

map_auto <- automixfit(map)

## advanced: look at AIC of all fitted models
sapply(attr(map_auto, "models"), AIC)

print(map_auto)

## 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))
)
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
}

ggplot(post, aes(r, mean, colour = prior, shape = prior)) +
  geom_point() +
  geom_abline(slope = 1 / 20)
ggplot(post, aes(r, sd, colour = prior, shape = prior)) +
  geom_point() +
  coord_cartesian(ylim = c(0, 0.17))

sessionInfo()