adoptr_jss.R
In adoptr: Adaptive Optimal Two-Stage Designs

## ----setup, include=FALSE-------------------------------------------
library(adoptr)
library(dplyr)
library(tidyr)
library(ggplot2)

buildvignette <- as.logical(Sys.getenv("NOTCRAN", unset = FALSE))

knitr::opts_chunk$set(
    engine='R', 
    highlight=FALSE, 
    prompt=TRUE, 
    tidy=FALSE,
    eval=buildvignette
)
backup_options <- options()
options(prompt = "R> ", continue = "+  ", width = 70, useFancyQuotes = FALSE)

## ----plot-opts, include=FALSE---------------------------------------
#  add_opts <- function(p) {
#    p + theme_bw() +
#    theme(
#      panel.grid   = element_blank(),
#      legend.title = element_blank(),
#      legend.position = "bottom"
#    )
#  }
#  
#  
#  # extract legend
#  g_legend <- function(a.gplot){
#    tmp <- ggplot_gtable(ggplot_build(a.gplot))
#    leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
#    legend <- tmp$grobs[[leg]]
#    return(legend)
#  }
#  

## ----class-diagram, out.width='100%', echo=FALSE, fig.cap='Overview of the most important classes and methods (in italic) in the R-package adoptr. A subclassing relationship is indicated by a connecting line to the  corresponding super class above it. The most important methods for each class are listed under the respective class name in italic font.'----
#  knitr::include_graphics("structure.png")

## ----case-1-data-distribution---------------------------------------
#  datadist <- Normal(two_armed = TRUE)

## ----case-1-hypotheses----------------------------------------------
#  null        <- PointMassPrior(theta = .0, mass = 1.0)
#  alternative <- PointMassPrior(theta = .3, mass = 1.0)
#  power       <- Power(dist = datadist, prior = alternative)
#  toer        <- Power(dist = datadist, prior = null)
#  mss         <- MaximumSampleSize()

## ----case-1-ess-----------------------------------------------------
#  ess <- ExpectedSampleSize(dist = datadist, prior = alternative)

## -------------------------------------------------------------------
#  initial_design <- get_initial_design(theta = 0.3, alpha = 0.025,
#                                       beta = 0.1, type = "two-stage",
#                                       dist = datadist, order = 7)

## ----collapse=FALSE-------------------------------------------------
#  evaluate(toer, initial_design)
#  evaluate(power, initial_design)

## ----collapse=FALSE-------------------------------------------------
#  evaluate(toer  <= .025, initial_design)
#  evaluate(power >= .9, initial_design)

## ----case-1-optimization--------------------------------------------
#  opt1 <- minimize(ess, subject_to(power >= 0.9, toer  <= 0.025),
#                   initial_design)

## -------------------------------------------------------------------
#  cp <- ConditionalPower(dist = datadist, prior = alternative)
#  summary(opt1$design, "Power" = power, "ESS" = ess, "CP" = cp)

## ----standard-case, fig.height=2.25*1.25, fig.width=6*1.25, out.width='100%', fig.cap="Optimal sample size, critical value, and conditional power plotted against the interim test statistic (built-in plot method)."----
#  plot(opt1$design, `Conditional power` = cp)

## -------------------------------------------------------------------
#  prior <- ContinuousPrior(
#    pdf     = function(theta) dnorm(theta, mean = .3, sd = .1),
#    support = c(-1, 1),
#    order   = 25)

## ----case-2-ess-----------------------------------------------------
#  ess <- ExpectedSampleSize(dist = datadist, prior = prior)

## ----case-2-pwr-mcr-------------------------------------------------
#  epower <- Power(dist = datadist, prior = condition(prior, c(.1, 1)))

## ----case-2-evaluate-epwr-opt1--------------------------------------
#  evaluate(epower, opt1$design)

## ----case-2-optimization, warning=FALSE-----------------------------
#  opt2 <- minimize(ess, subject_to(epower >= 0.9, toer <= 0.025),
#                   initial_design,
#                   opts = list(algorithm = "NLOPT_LN_COBYLA",
#                               xtol_rel = 1e-5, maxeval = 20000))

## -------------------------------------------------------------------
#  `n(X_1)`      <- ConditionalSampleSize()
#  `E[n(X_1)^2]` <- expected(composite({`n(X_1)`^2}),
#                            data_distribution = datadist,
#                            prior = prior)

## ----case-3-optimization--------------------------------------------
#  opt3 <- minimize(composite({`E[n(X_1)^2]` - 200000*epower}),
#                   subject_to(toer <= 0.025), initial_design)

## ----power-utility--------------------------------------------------
#  evaluate(epower, opt3$design)

## ----comparison, echo=FALSE, warning=FALSE, fig.height=2.25*1.25, fig.width=6*1.25, out.width='100%', fig.cap="Comparison of optimal designs under a point prior, a continuous prior, and a utility maximization approach."----
#  
#  x1    <- seq(0, 3, by = .01)
#  theta <- seq(0, .4, by = .01)
#  
#  plot_data_1 <- tibble(
#      type   = c("Point Prior", "Continuous Prior", "Utility Maximization"),
#      design = c(opt1$design, opt2$design, opt3$design)
#    ) %>%
#    group_by(type) %>%
#    do(
#      x1 = x1,
#      n  = adoptr::n(.$design[[1]], x1),
#      c2 = c2(.$design[[1]], x1)
#    )  %>%
#    unnest(., cols = c(x1, n, c2)) %>%
#    mutate(
#      section = ifelse(
#        is.finite(c2),
#        "continuation",
#        ifelse(c2 == -Inf, "efficacy", "futility")
#      )
#    ) %>%
#    gather(variable, value, n, c2)
#  
#  plot_data_2 <- tibble(
#      type   = c("Point Prior", "Continuous Prior", "Utility Maximization"),
#      design = c(opt1$design, opt2$design, opt3$design)
#    ) %>%
#    group_by(type) %>%
#    do(
#      theta = theta,
#      pow   = sapply(theta,
#                     function(d) {
#                         evaluate(Power(datadist, PointMassPrior(d, 1)), .$design[[1]])
#                       })
#    )  %>%
#    unnest(., cols = c(theta, pow)) %>%
#    gather(variable, value, pow)
#  
#  
#  
#  
#  
#  p1 <- {
#    ggplot(filter(plot_data_1, variable == "n"), aes(x1, value, color = type)) +
#      xlab(expression(x[1])) +
#      geom_line(aes(group = interaction(section, type))) +
#      scale_y_continuous("n", limits = c(0, 550), breaks = seq(0, 500, 100))
#    } %>%
#    add_opts
#  
#  p2 <- {
#    ggplot(filter(plot_data_1, variable == "c2"), aes(x1, value, color = type)) +
#      xlab(expression(x[1])) +
#      geom_line(aes(group = interaction(section, type))) +
#      scale_y_continuous(expression(c[2]), breaks = seq(0, 5, by = .5))
#  } %>%
#    add_opts
#  
#  
#  p3 <- {
#    ggplot(filter(plot_data_2, variable == "pow"), aes(theta, value, color = type)) +
#      xlab(expression(theta)) +
#      geom_line(aes(group = type)) +
#      scale_y_continuous("Power", breaks = seq(0, 1, .1))
#  } %>%
#    add_opts
#  
#  
#  gridExtra::grid.arrange(
#    gridExtra::arrangeGrob(
#      p1 + theme(legend.position = "none"),
#      p2 + theme(legend.position = "none"),
#      p3 + theme(legend.position = "none"),
#      nrow = 1
#    ),
#    g_legend(p1),
#    nrow = 2,
#    heights = c(10, 1)
#  )

## ----case-4-hypotheses----------------------------------------------
#  prior <- PointMassPrior(theta = .3, mass = 1.0)
#  ess   <- ExpectedSampleSize(dist = datadist, prior = prior)
#  cp    <- ConditionalPower(dist = datadist, prior = prior)
#  power <- expected(cp, data_distribution = datadist, prior = prior)

## -------------------------------------------------------------------
#  opt4 <- minimize(ess, subject_to(toer <= 0.025, power >= 0.9, cp >= 0.8),
#                   initial_design)

## ----cp-constraint, echo=FALSE, warning=FALSE, fig.height=2.25*1.25, fig.width=6*1.25, out.width="100%", fig.cap="Optimal designs with and without conditional power constraint."----
#  x1 <- seq(0, 3, by = .01)
#  
#  plot_data_2 <- tibble(
#    type   = c("No CP constraint", "With CP constraint"),
#    design = c(opt1$design, opt4$design)
#  ) %>%
#      group_by(type) %>%
#      do(
#          x1 = x1,
#          n  = adoptr::n(.$design[[1]], x1),
#          c2 = c2(.$design[[1]], x1),
#          cp = evaluate(cp, .$design[[1]], x1)
#      )  %>%
#      unnest(., c(x1, n, c2, cp)) %>%
#      mutate(
#          section = ifelse(
#              is.finite(c2),
#              "continuation",
#              ifelse(c2 == -Inf, "efficacy", "futility")
#          )
#      ) %>%
#      gather(variable, value, n, c2, cp)
#  
#  add_opts <- function(p) {
#    p + xlab(expression(x[1])) +
#    theme_bw() +
#    theme(
#      panel.grid   = element_blank(),
#      legend.title = element_blank(),
#      legend.position = "bottom"
#    )
#  }
#  
#  
#  p1_2 <- {
#    ggplot(filter(plot_data_2, variable == "n"), aes(x1, value, color = type)) +
#        geom_line(aes(group = interaction(section, type))) +
#        scale_y_continuous("n", limits = c(0, 500), breaks = seq(0, 500, 100))
#  } %>%
#    add_opts
#  
#  p2_2 <- {
#    ggplot(filter(plot_data_2, variable == "c2"), aes(x1, value, color = type)) +
#        geom_line(aes(group = interaction(section, type))) +
#        scale_y_continuous(expression(c[2]), breaks = seq(0, 5, by = .5))
#  } %>%
#    add_opts
#  
#  p3_2 <- {
#    ggplot(filter(plot_data_2, variable == "cp"), aes(x1, value, color = type)) +
#        geom_line(aes(group = interaction(section, type))) +
#        scale_y_continuous("Conditional Power", breaks = seq(0, 1, by = .1))
#  } %>%
#    add_opts
#  
#  
#  
#  gridExtra::grid.arrange(
#    gridExtra::arrangeGrob(
#      p1_2 + theme(legend.position = "none"),
#      p2_2 + theme(legend.position = "none"),
#      p3_2 + theme(legend.position = "none"),
#      nrow = 1
#    ),
#    g_legend(p1_2),
#    nrow = 2,
#    heights = c(10, 1)
#  )

## -------------------------------------------------------------------
#  initial_design@n1  <- 80
#  initial_design@c1f <- 0
#  initial_design     <- make_fixed(initial_design, n1, c1f)

## ----message=FALSE--------------------------------------------------
#  opt5 <- minimize(ess, subject_to(toer <= 0.025, power >= 0.9),
#                   initial_design)

## ----tunable, echo=FALSE, warning=FALSE, fig.height=2.25*1.25, fig.width=6*1.25, out.width="100%", fig.cap="Comparison of fully optimal design and optimal design with fixed first-stage sample size."----
#  
#  x1    <- seq(-.5, 3, by = .01)
#  theta <- seq(-.2, .5, by = .01)
#  
#  plot_data_1 <- tibble(
#      type   = c("Optimal", "Optimal with fixed n1"),
#      design = c(opt1$design, opt5$design)
#    ) %>%
#    group_by(type) %>%
#    do(
#        x1 = x1,
#        n  = adoptr::n(.$design[[1]], x1),
#        c2 = c2(.$design[[1]], x1)
#    )  %>%
#    unnest(., c(x1, n, c2)) %>%
#    mutate(
#        section = ifelse(
#            is.finite(c2),
#            "continuation",
#            ifelse(c2 == -Inf, "efficacy", "futility")
#        )
#    ) %>%
#    gather(variable, value, n, c2)
#  
#  
#  
#  plot_data_2 <- tibble(
#      type   = c("Optimal", "Optimal with fixed n1"),
#      design = c(opt1$design, opt5$design)
#    ) %>%
#    group_by(type) %>%
#    do(
#        theta = theta,
#        pow   = sapply(theta,
#                       function(d) {
#                         evaluate(ExpectedSampleSize(datadist, PointMassPrior(d, 1)), .$design[[1]])
#                       })
#    )  %>%
#    unnest(., cols  = c(theta, pow)) %>%
#    gather(variable, value, pow)
#  
#  
#  
#  p1 <- {
#    ggplot(filter(plot_data_1, variable == "n"), aes(x1, value, color = type)) +
#        xlab(expression(x[1])) +
#        geom_line(aes(group = interaction(section, type))) +
#        scale_y_continuous("n", limits = c(0, 500), breaks = seq(0, 500, 100))
#  } %>%
#    add_opts
#  
#  p2 <- {
#    ggplot(filter(plot_data_1, variable == "c2"), aes(x1, value, color = type)) +
#        xlab(expression(x[1])) +
#        geom_line(aes(group = interaction(section, type))) +
#        scale_y_continuous(expression(c[2]), breaks = seq(0, 5, by = .5))
#  } %>%
#    add_opts
#  
#  
#  p3 <- {
#    ggplot(filter(plot_data_2, variable == "pow"), aes(theta, value, color = type)) +
#        xlab(expression(theta)) +
#        geom_line(aes(group = type)) +
#        scale_y_continuous("Expected Sample Size") +
#        geom_vline(xintercept = 0.3, colour = "grey", size=0.3)
#  } %>%
#    add_opts
#  
#  
#  gridExtra::grid.arrange(
#    gridExtra::arrangeGrob(
#      p1 + theme(legend.position = "none"),
#      p2 + theme(legend.position = "none"),
#      p3 + theme(legend.position = "none"),
#      nrow = 1
#    ),
#    g_legend(p1),
#    nrow = 2,
#    heights = c(10, 1)
#  )

## ----reset-options, include=FALSE-----------------------------------
#  options(backup_options)