---
title: 'Sequentially updating the likelihood of success of a Phase 3 pivotal time-to-event trial based on interim analyses or external information'
author: 
- name: Kaspar Rufibach, Paul Jordan, and Markus Abt
  affiliation: Department of Biostatistics, Roche Basel
date: '`r Sys.setlocale("LC_TIME", "C"); format(Sys.Date(), "%d %b %Y")`'
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Sequentially updating the likelihood of success of a Phase 3 pivotal time-to-event trial based on interim analyses or external information}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
bibliography: biblio.bib
---

```{r setup, include = FALSE}
now <- as.character(as.POSIXlt(Sys.time()))
today <- as.Date(substr(now, 1, 10))
now <- paste(today, " at ", substr(now, 12, 19), sep = "")

## set some knitr options
knitr::opts_chunk$set(echo = TRUE, cache = TRUE)
```


# Purpose
This document reproduces the computations in @rufibach_14 ([doi](http://dx.doi.org/10.1080/10543406.2014.972508)). All the references below to sections, figures, or tables are with respect to that paper. The functions used are all implemented in the `R` package @bpp and all the code below is also available in the help file of the function ```bpp_1interim```.

Further discussions around Bayesian Predictive Power (BPP), mainly about the choice of a prior, are provided in @rufibach_15 ([doi](http://dx.doi.org/10.1002/pst.1764)).


# Setup

We start with defining the parameters of the pivotal trial and the results from the two external studies: 

```{r, echo = TRUE, message = FALSE}
# load the package:
library(bpp)

# get the code for all the computations below:
?bpp_1interim

# specifications of the the pivotal trial
propA <- 0.5   # proportion of patients randomized to arm A
fac <- (propA * (1 - propA)) ^ (-1)
nevents <- c(0.5, 1) * 1600
finalSE <- sqrt(fac / nevents[2])
alphas <- c(0.001, 0.049)
za <- qnorm(1 - alphas / 2)
hrMDD <- exp(- za * sqrt(fac / nevents))
successmean <- log(hrMDD[2])

# efficacy and futility interim boundary
effi <- log(hrMDD[1])
futi <- log(1.025)

# specifications of the first external study
hr1 <- 0.396
sd1 <- 0.837

# specifications of the second external study
hr2 <- 0.287
sd2 <- 0.658

# grid to compute densities on
thetas <- seq(-0.65, 0.3, by = 0.01)
```


# Define the prior distributions
Then we define the parameters of the prior distributions, see Sections 4.1 and 4.2.

```{r, echo = TRUE, message = FALSE}
# ----------------------------------
# mean and sd of Normal prior:
# ----------------------------------
hr0 <- 0.85
sd0 <- 0.11

# all computations are done on the log(hazard ratio) scale, in order to 
# be able to approximate the distribution of the log(hazard ratio) via a
# Normal distribution:
priormean <- log(hr0)

# ----------------------------------
# parameters of pessimistic, or flat, prior:
# ----------------------------------
hr0flat <- 0.866
priormeanflat <- log(hr0flat)
width1 <- 0.21
height1 <- 2.48
```

Let us compare these two prior distributions:

```{r, echo = TRUE, results = 'asis', message = FALSE, fig.cap = "", fig.align = "center", fig.width = 7, fig.height = 5.5}
par(las = 1, mar = c(9, 5, 2, 1), mfrow = c(1, 1))
plot(0, 0, type = "n", xlim = c(-0.6, 0.3), ylim = c(-0.1, 5), xlab = "", ylab = "density", main = "")
title(expression("Normal and pessimistic prior density for "*theta), line = 0.7)
basicPlot(leg = FALSE, IntEffBoundary = NA, IntFutBoundary = NA, successmean = NA, priormean = NA)
lines(thetas, dnorm(thetas, mean = priormean, sd = sd0), col = 2, lwd = 2)
lines(thetas, dUniformNormalTails(thetas, mu = priormeanflat, width = width1, height = height1), lwd = 2, col = 3)
```


# Computation of all quantities in Table 1, and all the figures

In what follows, we compute all quantities that appear in Table 1 of @rufibach_14. In addition, the densities shown in Figures 1 and 2 are generated as well. All these computations are done for both, the Normal and the pessimistic prior.


## Interesting probabilities 

For the Normal prior:
```{r, echo = TRUE, message = FALSE}
# Normal prior probabilities to be below 0.7 or above 1:
lims <- c(0.7, 1)
pnorm1 <- plnorm(lims[1], meanlog = priormean, sdlog = sd0, lower.tail = TRUE, log.p = FALSE)   
pnorm2 <- plnorm(lims[2], meanlog = priormean, sdlog = sd0, lower.tail = FALSE, log.p = FALSE)
```

And the same for the pessimistic prior:
```{r, echo = TRUE, message = FALSE}
# pessimistic prior probabilities to be below 0.7 or above 1:
flat1 <- pUniformNormalTails(x = log(lims[1]), mu = priormeanflat, width = width1, height = height1)
flat2 <- 1 - pUniformNormalTails(x = log(lims[2]), mu = priormeanflat, width = width1, 
                                 height = height1)
```


## BPP at the beginning of the trial

Compute BPP based on the prior distributions and the specifications of the pivotal trial.

```{r, echo = TRUE, message = FALSE}
# ----------------------------------
# Normal prior:
# ----------------------------------
bpp0 <- bpp(prior = "normal", successmean = successmean, finalSE = finalSE, 
            priormean = priormean, priorsigma = sd0)

# ----------------------------------
# pessimistic prior:
# ----------------------------------
bpp0_1 <- bpp(prior = "flat", successmean = successmean, finalSE = finalSE, 
              priormean = priormeanflat, width = width1, height = height1)
```


## Update the prior distribution with external information

We have two external studies that define our likelihood with which we update the two prior distributions. 

```{r, echo = TRUE, message = FALSE}
# ----------------------------------
# Normal prior:
# ----------------------------------
up1 <- NormalNormalPosterior(datamean = log(hr1), sigma = sd1, n = 1, nu = priormean, 
                             tau = sd0)
bpp1 <- bpp(prior = "normal", successmean = successmean, finalSE = finalSE, 
            priormean = up1$postmean, priorsigma = up1$postsigma)

# update prior with second external study (result derived from pooled analysis: 
# Cox regression on patient level, stratified by study):
up2 <- NormalNormalPosterior(datamean = log(hr2), sigma = sd2, n = 1, nu = priormean, 
                             tau = sd0)
bpp2 <- bpp(prior = "normal", successmean = successmean, finalSE = finalSE, 
            priormean = up2$postmean, priorsigma = up2$postsigma)

# ----------------------------------
# pessimistic prior
# ----------------------------------
bpp1_1 <- integrate(FlatNormalPosterior, lower = -Inf, upper = Inf, successmean = successmean, 
                     finalSE = finalSE, interimmean = log(hr1), interimSE = sd1, 
                     priormean = priormeanflat, width = width1, height = height1,
                     subdivisions = 300)$value

bpp2_1 <- integrate(FlatNormalPosterior, -Inf, Inf, successmean = successmean, 
                     finalSE = finalSE, interimmean = log(hr2), 
                     interimSE = sd2, priormean = priormeanflat, 
                     width = width1, height = height1, subdivisions = 300)$value
```


## Update the prior distribution after not stopping at an interim analysis

The next quantities we are interested in is what happens to BPP if we do not stop at an interim analysis. We assess various scenarios, as in Table 1: interim is set up with futility and efficacy boundary, only one of the two, and we also provide updates using conditional power for the two extreme cases, namely that the estimate at the interim was known to be either equal to the futility or the efficacy boundary. Let us start with the Normal prior:

```{r, echo = TRUE, message = FALSE}
# ----------------------------------
# compute bpp after not stopping at interim, for Normal prior and various 
# assumptions on the amount of information we learn at the interim
# ----------------------------------

# assuming both boundaries:
bpp3.tmp <- bpp_1interim(prior = "normal", interimSE = sqrt(fac / nevents[1]), 
                         finalSE = finalSE, successmean = successmean, 
                         IntEffBoundary = effi, IntFutBoundary = futi, IntFix = 1, 
                         priormean = up2$postmean, propA = 0.5, thetas = thetas, 
                         priorsigma = up2$postsigma)
bpp3 <- bpp3.tmp$"BPP after not stopping at interim interval"
post3 <- bpp3.tmp$"posterior density interval"

# assuming only efficacy boundary:
bpp3_effi_only <- bpp_1interim(prior = "normal", interimSE = sqrt(fac / nevents[1]), 
                               finalSE = finalSE, successmean = successmean, 
                               IntEffBoundary = effi, IntFutBoundary = log(Inf), IntFix = 1, 
                               priormean = up2$postmean, propA = 0.5, thetas = thetas, 
                               priorsigma = 
                               up2$postsigma)$"BPP after not stopping at interim interval"

# assuming only futility boundary:
bpp3_futi_only <- bpp_1interim(prior = "normal", interimSE = sqrt(fac / nevents[1]), 
                               finalSE = finalSE, successmean = successmean, 
                               IntEffBoundary = log(0), IntFutBoundary = futi, IntFix = 1, 
                               priormean = up2$postmean, propA = 0.5, thetas = thetas, 
                               priorsigma = 
                               up2$postsigma)$"BPP after not stopping at interim interval"

# assuming interim efficacy boundary: 
bpp4.tmp <- bpp_1interim(prior = "normal", interimSE = sqrt(fac / nevents[1]), 
                         finalSE = finalSE, successmean = successmean, 
                         IntEffBoundary = effi, IntFutBoundary = Inf, IntFix = c(effi, futi), 
                         priormean = up2$postmean, propA = 0.5, thetas = thetas, 
                         priorsigma = up2$postsigma)
bpp4 <- bpp4.tmp$"BPP after not stopping at interim exact"[2, 1]
post4 <- bpp4.tmp$"posterior density exact"[, 1]

# assuming interim futility boundary: 
bpp5.tmp <- bpp_1interim(prior = "normal", interimSE = sqrt(fac / nevents[1]), 
                         finalSE = finalSE, successmean = successmean, 
                         IntEffBoundary = effi, IntFutBoundary = Inf, IntFix = futi, 
                         priormean = up2$postmean, propA = 0.5, thetas = thetas, 
                         priorsigma = up2$postsigma)
bpp5 <- bpp5.tmp$"BPP after not stopping at interim exact"[2, 1]
post5 <- bpp5.tmp$"posterior density exact"     # same as post4[, 2]
```


And the same for the pessimistic prior:

```{r, echo = TRUE, message = FALSE}
# ----------------------------------
# compute bpp after not stopping at interim, for pessimistic prior and various 
# assumptions on the amount of information we learn at the interim
# ----------------------------------

# assuming both boundaries:
bpp3.tmp_1 <- bpp_1interim(prior = "flat", interimSE = sqrt(fac / nevents[1]), 
                         finalSE = finalSE, successmean = successmean, 
                         IntEffBoundary = effi, IntFutBoundary = futi, IntFix = 1, 
                         priormean = up2$postmean, propA = 0.5, thetas = thetas, 
                         width = width1, height = height1)
bpp3_1 <- bpp3.tmp_1$"BPP after not stopping at interim interval"
post3_1 <- bpp3.tmp_1$"posterior density interval"

# assuming only efficacy boundary:
bpp3_1_effi_only <- bpp_1interim(prior = "flat", interimSE = sqrt(fac / nevents[1]), 
                               finalSE = finalSE, successmean = successmean, 
                               IntEffBoundary = effi, IntFutBoundary = log(Inf), IntFix = 1, 
                               priormean = up2$postmean, propA = 0.5, thetas = thetas, 
                               width = width1, 
                               height = height1)$"BPP after not stopping at interim interval"

# assuming only futility boundary:
bpp3_1_futi_only <- bpp_1interim(prior = "flat", interimSE = sqrt(fac / nevents[1]), 
                               finalSE = finalSE, successmean = successmean, 
                               IntEffBoundary = log(0), IntFutBoundary = futi, IntFix = 1, 
                               priormean = up2$postmean, propA = 0.5, thetas = thetas, 
                               width = width1, 
                               height = height1)$"BPP after not stopping at interim interval"

# assuming interim efficacy boundary: 
bpp4_1.tmp <- bpp_1interim(prior = "flat", interimSE = sqrt(fac / nevents[1]), 
                               finalSE = finalSE, successmean = successmean, 
                               IntEffBoundary = log(0), IntFutBoundary = effi, IntFix = effi, 
                               priormean = up2$postmean, propA = 0.5, thetas = thetas, 
                               width = width1, height = height1)
bpp4_1 <- bpp4_1.tmp$"BPP after not stopping at interim exact"[2, 1]
post4_1 <- bpp4_1.tmp$"posterior density exact"

# assuming interim futility boundary: 
bpp5_1 <- integrate(Vectorize(estimate_toIntegrate), lower = -Inf, upper = Inf, prior = "flat",
                    successmean = successmean, finalSE = finalSE, interimmean = futi, 
                    interimSE = sqrt(fac / nevents[1]), priormean = up2$postmean, 
                    width = width1, height = height1, subdivisions = 300)$value


bpp5_1.tmp <- bpp_1interim(prior = "flat", interimSE = sqrt(fac / nevents[1]), 
                           finalSE = finalSE, successmean = successmean, 
                           IntEffBoundary = log(0), IntFutBoundary = effi, IntFix = futi, 
                               priormean = up2$postmean, propA = 0.5, thetas = thetas, 
                               width = width1, height = height1)
bpp5_1 <- bpp5_1.tmp$"BPP after not stopping at interim exact"[2, 1]
post5_1 <- bpp5_1.tmp$"posterior density exact"
```


## Posterior densities

The resulting posteriors after the updates with the external studies are depicted below, compare Figure 1 in @rufibach_14. 

```{r, echo = TRUE, results = 'asis', message = FALSE, fig.cap = "", fig.align = "center", fig.width = 7, fig.height = 5}
par(las = 1, mar = c(9, 5, 2, 1), mfrow = c(1, 2), cex = 0.8)

# ----------------------------------
# Normal prior:
# ----------------------------------
plot(0, 0, type = "n", xlim = c(-0.6, 0.3), ylim = c(-0.1, 5), xlab = "", ylab = "density", 
     main = "")
title(expression("Normal prior density and corresponding posteriors for "*theta), line = 0.7)
basicPlot(leg = FALSE, IntEffBoundary = effi, IntFutBoundary = futi, successmean = successmean, 
          priormean = priormean)
lines(thetas, dnorm(thetas, mean = priormean, sd = sd0), col = 2, lwd = 2)
lines(thetas, dnorm(thetas, mean = up1$postmean, sd = up1$postsigma), col = 3, lwd = 2)
lines(thetas, dnorm(thetas, mean = up2$postmean, sd = up2$postsigma), col = 4, lwd = 2)
lines(thetas, post3, col = 1, lwd = 2)
legend(-0.64, 5.2, c("prior", "posterior after Sub1", "posterior after Sub1 & Sub2", 
       "posterior after Sub1 & Sub2 and not stopping at interim"), lty = 1, col = c(2:4, 1), 
       bty = "n", lwd = 2)

# ----------------------------------
# pessimistic prior:
# ----------------------------------

# first we have to compute the posteriors after the external updates:
flatpost1 <- rep(NA, length(thetas))
flatpost2 <- flatpost1
for (i in 1:length(thetas)){
  flatpost1[i] <- estimate_posterior(x = thetas[i], prior = "flat", interimmean = log(hr1), 
                                     interimSE = sd1, priormean = priormeanflat, width = width1, 
                                     height = height1)
  flatpost2[i] <- estimate_posterior(x = thetas[i], prior = "flat", interimmean = log(hr2), 
                                     interimSE = sd2, priormean = priormeanflat, width = width1, 
                                     height = height1)
}

plot(0, 0, type = "n", xlim = c(-0.6, 0.3), ylim = c(-0.10, 5), xlab = "", ylab = "density", main = "")
title(expression("Flat prior density and corresponding posteriors for "*theta), line = 0.7)
basicPlot(leg = FALSE, IntEffBoundary = effi, IntFutBoundary = futi, successmean = successmean, 
          priormean = priormeanflat)
lines(thetas, dUniformNormalTails(thetas, mu = priormeanflat, width = width1, height = height1), 
      lwd = 2, col = 2)
lines(thetas, flatpost1, col = 3, lwd = 2)
lines(thetas, flatpost2, col = 4, lwd = 2)
lines(thetas, post3_1, col = 1, lwd = 2)

legend(-0.64, 5.2, c("prior", "posterior after Sub1", "posterior after Sub1 & Sub2", 
       "posterior after Sub1 & Sub2 and not stopping at interim"), lty = 1, col = c(2:4, 1), 
       bty = "n", lwd = 2)
```


Next, the posteriors after updating with the interim information, as in Figure 2:

```{r, echo = TRUE, results = 'asis', message = FALSE, fig.cap = "", fig.align = "center", fig.width = 7, fig.height = 5.5}
par(las = 1, mar = c(9, 5, 2, 1), mfrow = c(1, 2), cex = 0.8)

# ----------------------------------
# Normal prior:
# ----------------------------------
plot(0, 0, type = "n", xlim = c(-0.6, 0.3), ylim = c(-0.1, 8), xlab = "", ylab = "density", 
     main = "")
title("Posteriors for after not stopping at interim, Normal prior", line = 0.7)
basicPlot(leg = FALSE, IntEffBoundary = effi, IntFutBoundary = futi, successmean = successmean, 
          priormean = priormean)
lines(thetas, post3, col = 1, lwd = 2)
lines(thetas, post4, col = 2, lwd = 2)
lines(thetas, post5, col = 3, lwd = 2)
leg2 <- c("interval knowledge", expression(hat(theta)*" = efficacy boundary"), 
          expression(hat(theta)*" = futility boundary"))
legend(-0.62, 8.2, leg2, lty = 1, col = 1:3, lwd = 2, bty = "n", 
       title = "posterior after not stopping at interim,")

# ----------------------------------
# pessimistic prior:
# ----------------------------------
plot(0, 0, type = "n", xlim = c(-0.6, 0.3), ylim = c(-0.10, 8), xlab = "", ylab = "density", 
     main = "")
title("Posteriors after not stopping at interim, pessimistic prior", line = 0.7)
basicPlot(leg = FALSE, IntEffBoundary = effi, IntFutBoundary = futi, successmean = successmean, 
          priormean = priormeanflat)
lines(thetas, post3_1, col = 1, lwd = 2)
lines(thetas, post4_1, col = 2, lwd = 2)
lines(thetas, post5_1, col = 3, lwd = 2)
leg.flat <- c("interval knowledge", expression(hat(theta)*" = efficacy boundary"), 
              expression(hat(theta)*" = futility boundary"))
legend(-0.62, 8.2, leg.flat, lty = 1, col = 1:3, lwd = 2, bty = "n", 
       title = "posterior after not stopping at interim,")
```


## Table 1 in @rufibach_14

Finally, we collect all the results computed above in one table, reproducing Table 1 in @rufibach_14.

```{r, echo = TRUE, results = 'asis', message = FALSE}
mat <- matrix(NA, ncol = 2, nrow = 10)
mat[, 1] <- c(pnorm1, pnorm2, bpp0, bpp1, bpp2, bpp3, bpp3_futi_only, bpp3_effi_only, bpp4, bpp5)
mat[, 2] <- c(flat1, flat2, bpp0_1, bpp1_1, bpp2_1, bpp3_1, bpp3_1_futi_only, bpp3_1_effi_only, 
              bpp4_1, bpp5_1)
colnames(mat) <- c("Normal prior", "Flat prior")
rownames(mat) <- c(paste("Probability for hazard ratio to be $\\le$ ", lims[1], sep = ""), 
paste("Probability for hazard ratio to be $\\ge$ ", lims[2], sep = ""), 
"PoS based on prior distribution", "PoS after Sub1", "PoS after Sub1 and Sub2", 
"PoS after not stopping at interim, assuming $\\hat \\theta \\in [\\theta_{futi}, \\theta_{effi}]$", 
"PoS after not stopping at interim, assuming $\\hat \\theta \\in [-\\infty, \\theta_{futi}]$", 
"PoS after not stopping at interim, assuming $\\hat \\theta \\in [\\theta_{effi}, \\infty]$", 
"PoS after not stopping at interim, assuming $\\hat \\theta = \\theta_{effi}$", 
"PoS after not stopping at interim, assuming $\\hat \\theta = \\theta_{futi}$")
mat <- as.data.frame(format(mat, digits = 2))
library(knitr); kable(mat)
```

# Update of BPP when not stopping the trial in two blinded interim analyses

The theory developed in @rufibach_14 can straightforwardly be extended to accomodate more than one interim analysis. Below, we illustrate the BPP update assuming a trial that did not stop after first a futility and second an efficacy interim analysis. Note that the function ```bpp_2interim``` in @bpp so far only allows specification of a Normal prior.

```{r, echo = TRUE, message = FALSE}
# ------------------------------------------------------------------------------------------
# Illustrate the update after two passed interims using the Gallium clinical trial
# ------------------------------------------------------------------------------------------

# mean and sd of Normal prior:
hr0 <- 0.9288563
priormean <- log(hr0)
priorsigma <- sqrt(4 / 12)

# specifications for pivotal study:
propA <- 0.5   
fac <- (propA * (1 - propA)) ^ (-1)
nevents <- c(111, 248, 370)
interimSE <- sqrt(fac / nevents[1:2])
finalSE <- sqrt(fac / nevents[3])
za <- c(3.9285726330559, 2.5028231888636, 1.9936294555664)
alphas <- 2 * (1 - pnorm(za))
hrMDD <- exp(- za * sqrt(fac / nevents))
successmean <- log(hrMDD[3])

# 2-d vector of efficacy and futility interim boundaries:
effi <- log(c(0, hrMDD[2]))  # first interim is for futility only
futi <- log(c(1, Inf))       # second interim is for efficacy only

bpp2 <- bpp_2interim(prior = "normal", interimSE = interimSE, finalSE = finalSE, 
             successmean = successmean, IntEffBoundary = effi, IntFutBoundary = futi, 
             priormean = priormean, thetas = thetas, priorsigma = priorsigma)

bpp2$"initial BPP"
bpp2$"BPP after not stopping at interim interval"
```


# Info
R version and packages used to generate this report:
 
R version: `r sessionInfo()$R.version$version.string`
 
Base packages: `r paste(sessionInfo()$basePkgs, collapse = " / ")`
 
Other packages: `r paste(names(sessionInfo()$otherPkgs), collapse = " / ")`
 
This document was generated on `r now`.


# References


