customizing_plots.R
In RBesT: R Bayesian Evidence Synthesis Tools

## ---- SETTINGS-knitr, include=FALSE-------------------------------------------
## knitr settings used to build vignettes
library(RBesT)
library(knitr)
library(ggplot2)
theme_set(theme_bw())
knitr::knit_hooks$set(pngquant = knitr::hook_pngquant)
knitr::opts_chunk$set(
  dev = "ragg_png",
  dpi = 72,
  fig.retina = 2,
  fig.width = 1.62*4,
  fig.height = 4,
  fig.align = "center",
  out.width = "100%",
  pngquant = "--speed=1 --quality=50"
  )

## ---- SETTINGS-sampling, include=FALSE----------------------------------------
## sampling settings used to build vignettes
## setup up fast sampling when run on CRAN
is_CRAN <- Sys.getenv("NOT_CRAN", "true") != "true"
## NOTE: for running this vignette locally, please uncomment the
## following line:
## is_CRAN <- FALSE
.user_mc_options <- list()
if (is_CRAN) {
    .user_mc_options <- options(RBesT.MC.warmup=250, RBesT.MC.iter=500, RBesT.MC.chains=2, RBesT.MC.thin=1, RBesT.MC.control=list(adapt_delta=0.9))
}
set.seed(6475863)

## -----------------------------------------------------------------------------
# Load required libraries
library(RBesT)
library(ggplot2)
library(dplyr)
library(tidyr)
library(bayesplot)

# Default settings for bayesplot
color_scheme_set("blue")
theme_set(theme_default(base_size=12))

# Load example gMAP object
set.seed(546346)
map_crohn <- gMAP(cbind(y, y.se) ~ 1 | study,
                  family=gaussian,
                  data=transform(crohn, y.se=88/sqrt(n)),
                  weights=n,
                  tau.dist="HalfNormal", tau.prior=44,
                  beta.prior=cbind(0,88))

print(map_crohn)

## -----------------------------------------------------------------------------
forest_plot(map_crohn)

## -----------------------------------------------------------------------------
forest_plot(map_crohn, model="both") + legend_move("right")

## -----------------------------------------------------------------------------
# preview a color scheme
color_scheme_view("mix-blue-red")
# and now let's use it
color_scheme_set("mix-blue-red") 
forest_plot(map_crohn)
color_scheme_set("gray")
forest_plot(map_crohn)

## -----------------------------------------------------------------------------
color_scheme_set("blue")
forest_plot(map_crohn, size=0.5, alpha=0)

## ---- fig.width=10, fig.height=6----------------------------------------------
# adjust the base font size
theme_set(theme_default(base_size=16))
forest_plot(map_crohn, model="both", est="MAP", size=1) + legend_move("right") +
  labs(title="Forest plot", subtitle="Results of Meta-Analytic-Predictive (MAP) analysis", 
       caption="Plot shows point estimates (posterior medians) with 95% intervals")

## ---- eval=FALSE--------------------------------------------------------------
#  ggsave("plot1.png", last_plot(), width=1.62*2.78, height=2.78, unit="in") # too small for the chosen font size
#  ggsave("plot2.png", last_plot(), width=1.62*5.56, height=5.56, unit="in") # fits a single ppt slide quite well

## ---- echo=FALSE--------------------------------------------------------------
theme_set(theme_default(base_size=12))
forest_plot(map_crohn) + labs(title="Original forest plot")

## -----------------------------------------------------------------------------
# Extract the data from the returned object
fp_data <- forest_plot(map_crohn)$data
print(fp_data, digits=2)
# Use a two-component map mixture to compute the vertical line location
map_mix <- mixfit(map_crohn, Nc=2)
# Finally compose a ggplot call for the desired graph
ggplot(fp_data, aes(x=study, y=median, ymin=low, ymax=up, linetype=model, shape=model)) + 
  geom_pointrange(size=0.7, position=position_dodge(width=0.5)) + 
  geom_hline(yintercept=qmix(map_mix, 0.5), linetype=3, alpha=0.5) +
  coord_flip() + theme_bw(base_size=12) + theme(legend.position="None") + 
  labs(x="", y="Response", title="Modified forest plot")

## -----------------------------------------------------------------------------
# Define prior distributions
prior_pbo <- mixbeta(inf1=c(0.60, 19, 29), inf2=c(0.30, 4, 5), rob=c(0.10, 1, 1))
prior_trt <- mixbeta(c(1, 1/3, 1/3))

# Study sample size
n_trt <- 50
n_pbo <- 20

# Create decision rules and designs to represent success and futility
success <- decision2S(pc=c(0.90, 0.50), qc=c(log(1), log(0.50)), lower.tail=TRUE, link="log")
futility <- decision2S(pc=c(0.10, 0.50), qc=c(log(1), log(0.50)), lower.tail=FALSE, link="log")
design_suc <- oc2S(prior_trt, prior_pbo, n_trt, n_pbo, success)
design_fut <- oc2S(prior_trt, prior_pbo, n_trt, n_pbo, futility)
crit_suc <- decision2S_boundary(prior_trt, prior_pbo, n_trt, n_pbo, success)
crit_fut <- decision2S_boundary(prior_trt, prior_pbo, n_trt, n_pbo, futility)

# Create a data frame that holds the outcomes for y1 (treatment) that define success and futility, 
# conditional on the number of events on y2 (placebo)
outcomes <- data.frame(y2=c(0:n_pbo), suc=crit_suc(0:n_pbo), fut=crit_fut(0:n_pbo), max=n_trt)
outcomes$suc <- with(outcomes, ifelse(suc<0, 0, suc)) # don't allow negative number of events

# Finally put it all together in a plot. 
o <- 0.5 # offset
ggplot(outcomes, aes(x=y2, ymin=-o, ymax=suc+o)) + geom_linerange(size=4, colour="green", alpha=0.6) + 
  geom_linerange(aes(ymin=suc+o, ymax=fut+o), colour="orange", size=4, alpha=0.6) + 
  geom_linerange(aes(ymin=fut+o, ymax=max+o), colour="red", size=4, alpha=0.6) + 
  annotate("text", x=c(2,14), y=c(36,8), label=c("STOP", "GO"), size=10) + 
  scale_x_continuous(breaks=seq(0,n_pbo,by=2)) + 
  scale_y_continuous(breaks=seq(0,n_trt,by=4)) + 
  labs(x="Events on placebo", y="Events on treatment", title="Study outcomes") + 
  coord_flip() + theme_bw(base_size=12)

## ---- fig.width=8.5, fig.height=5.5-------------------------------------------
# Define the grid of true event rates for which to evaluate OC
p_trt <- seq(0, 0.5, length=200)
p_pbo <- c(0.35, 0.40, 0.45, 0.50)

# Loop through the values for placebo and compute outcome probabilities 
oc_list <- lapply(p_pbo, function(x) {
  p_go <- design_suc(p_trt, x)
  p_stop <- design_fut(p_trt, x)
  data.frame(p_trt, p_pbo=x, Go=p_go, Stop=p_stop, Indeterminate=1-p_go-p_stop)
})
# The above returns a list, so we bind the elements together into one data frame
oc <- bind_rows(oc_list) 
# And convert from wide to long format
oc <- gather(oc, "Outcome", "Probability", 3:5)
oc$facet_text <- as.factor(paste("True placebo rate = ", oc$p_pbo, sep=""))

# Finally we are ready to plot
ggplot(oc, aes(x=p_trt, y=Probability, colour=Outcome, linetype=Outcome)) + 
  facet_wrap(~facet_text) + geom_line(size=1) + 
  scale_colour_manual(values=c("green", "orange", "red"), name="Outcome") + 
  scale_linetype(guide=FALSE) + 
  geom_hline(yintercept=c(0.1, 0.8), linetype=3) + 
  scale_y_continuous(breaks=seq(0, 1, by=0.2)) + 
  labs(x="True event rate for treatment", y="Probability", title="Operating Characteristics") + 
  theme_bw(base_size=12) 

## -----------------------------------------------------------------------------
sessionInfo()