---
title: "Probability of Success with Co-Data"
subtitle: "Use of Co-Data in Clinical Trails with RBesT"
author: "Sebastian Weber"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
  html_vignette:
    toc: true
  html_document:
    toc: true
    number_sections: true
    toc_float:
      collapsed: false
      smooth_scroll: false
  pdf_document:
    toc: true
  word_document:
    toc: true
vignette: >
  %\VignetteIndexEntry{Probability of Success with Co-Data}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
---

```{r, include=FALSE}
if(!require(RBesT)) {
    ## in case RBesT is not yet installed, please
    ## install RBesT from CRAN
    stop("RBesT R package is not installed. Please use install.packages to install.")
    ## example call to install.packages
    ## install.packages("RBesT")
    ## the following packages are recommended as well
    ## install.packages(c("bayesplot", "knitr", "dplyr", "purrr", "ggplot2"))
}
library(RBesT)
library(knitr)
library(bayesplot)
library(dplyr)
library(purrr)
library(ggplot2)
theme_set(theme_bw())
knitr::opts_chunk$set(
    fig.width = 1.62*4,
    fig.height = 4,
    cache=TRUE
)
## setup up fast sampling when run on CRAN
is_CRAN <- !identical(Sys.getenv("NOT_CRAN"), "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=50, RBesT.MC.iter=100, RBesT.MC.chains=2, RBesT.MC.thin=1)
}
```

<!-- To produce a html report, please execute in R -->
<!-- library(rmarkdown) -->
<!-- render("PoS_codata.Rmd") -->


The probability of success is a very useful concept to assess the
chances of success of a trial while taking uncertainties into
account. In this document we briefly introduce the needed formal
details, describe an example data set and then demonstrate various
ways of increasing complexity how the probability of success concept
can be applied to increasingly complex data situations.

# Introduction

The co-data concept has been introduced in [1]. It differs from the
use of historical data in that the approach makes use of contemporary
data. A meta-analytic-predictive (MAP) analysis assumes that
historical data is known at the time-point of specyfing the analysis
and is as such a retrospective summary of available data. The MAP
prior is then combined with the current trial data. A co-data approach
extends this sequential procedure to a meta-analytic-combined (MAC)
analysis. In the MAC approach all available data is analyzed in a
single step - that is, historical and concurrent data is combined in a
single inference step. Both approaches MAP and MAC yield exactly the
same results as is demonstrated in the appendix at the bottom. An
example for a co-data scenario in drug development is the simultaneous
execution of twin phase III trails for registriation. In such a
setting, a futility analysis at an interim analysis may take
historical an all contemporary data into account in a co-data
approach. This example has been discussed in [1] using the probability
of success (PoS) as metric to assess futility at an interim analysis
and is discussed here in detail.

The key property of the probability of success metric is the
consideration of uncertainty in parameters conditional on available
data. In contrast, the conditional power (CP) calculates the frequency
a given experemintal design will be successful for a known value of
the parameters. For example, in a 1-sample experiment with a one-sided
success criterion the trial is successful if the collected data $y_N$
of sample size $N$ exceeds some critical value $y_c$; recall that the
critical value $y_c$ is determined by the success criterion, prior and
sample size when evaluated. Assuming that the sampling model of the
data is $p(y|\theta)$, then

$$ CP_N(\theta) = \int I(y_N  > y_c) \, p(y_N|\theta) \, dy_N.$$

The integration over the data $y_N$ comprises all possible outcomes of
the trial. Note that before the start of the trial, $CP_N(\theta_0)$
is the type I error rate under the conventional null hypothesis
($\theta=\theta_0$) and $CP_N(\theta_a)$ the power of the trial under
the alternative ($\theta=\theta_a$). At an interim analysis at sample
size $n_I$, the conditional power is then evaluated conditional on the
observed data so far (the $n_I$ measurements) while the remaining
sample size ($N-n_I$) is random and distributed according to
$p(y|\theta)$ with $\theta$ set to some known value,

$$ CP_{N-n_I}(\theta|y_{n_I}) = \int I(y_{n_I} + y_{N-n_I} = y_N > y_c|y_{n_I}) \, p(y_{N-n_I}|\theta) \, dy_{N-n_I}.$$

The known value can be set equal to the observed point estimate at the
interim, $CP_{N-n_I}(\hat{\theta}_{I}|y_{n_I})$, or to the assumed
true alternative, $CP_{N-n_I}(\theta_a|y_{n_I})$, used to plan the
trial.

The probability of success in contrast assigns $\theta$ a distribution
and marginalizes the conditional power over this distribution. In
absence of additional trial external information this distribution is
the posterior for $p(\theta|y_{n_I})$ obtained from the prior for
$p(\theta)$ and the data collected up to the interim,

$$ PoS_I = \int CP_{N-n_I}(\theta|y_{n_I}) \, p(\theta|y_{n_I}) \, d\theta.$$

However, our knowledge about $\theta$ can be refined if other
data-sources like completed (historical) or concurrent trials are
available,

$$ PoS_{I,H,...} = \int CP_{N-n_I}(\theta|y_{n_I}) \, p(\theta|y_{I},y_{H},...) \, d\theta.$$

It is important to note that the conditional power is *always*
evaluated with respect to the trial data (and prior) only. Thus,
additional data-sources are not part of the analysis of the trial. In
practice this means that the probability of success is usually
calculated for a trial which uses non-informative priors, but at
interim we may use additional data-sources to refine our knowledge on
$\theta$ which will not be part of the trial analysis.

# Example Data Scenario

In the following the hypothetical example as in [1] is discussed. The
assumed endpoint is time-to-event, which is analyzed using the normal
approximation of the log-rank statistic for comparing two
groups. Under a 1:1 randomization the sampling distribution of the
log-hazard ratio is a normal distribution with sampling standard
deviation $2$ and a standard error of $\sqrt{4/N_{events}}$ for a
sample size of $N_{events}$. The historical data considered is a proof
of concept and a phase II trial. The twin phase III studies are event
driven. Each trial stops whenever a total of $379$ events is reached
and an interim is planned whenever at least $150$ events have
occured. The assumed true hazard ratio used for the design of the
trial is $0.8$.

Example data:
```{r}
trials <- data.frame(study  = c("PoC", "PhII", "PhIII_A", "PhIII_B"),
                     deaths = c(    8,     85,       162,       150),
                     HR     = c(  0.7,   0.75,      0.83,      0.78),
                     stringsAsFactors = FALSE
                     )
## under the normal approximation of the log-HR, the sampling sd is 2
## such that the standard errors are sqrt(4/events)
trials <- trials %>%
    mutate(logHR=log(HR), sem=sqrt(4/deaths))
kable(trials, digits=2)
```

The remaining outline of the vignette is to first evaluate the design
properties of the trials, then calculate the probability of success at
interim for the A trial only (using only trial A data). Next, the
probability of success is calculated using in addition the historical
data. Subsequently, the probability of success for trial A is
calculated also using the historical data *and* the concurrent phase
III data of trial B. Finally, the overall probability of success is
calculated, which is defined by the joint success of both trials. In
the appendix the equivalence of the MAP and MAC approach is
demonstrated.


# Statistical Trial Design Considerations

Key design choices:

- time-to-event endpoint
- phase III trials stop at target # of events 379
- null hypothesis of no difference in HR, $\theta_0 = 1.0$
- one-sided $\alpha = 0.025$
- alternative hypothesis assumes true HR of $\theta_a=0.75$
- interim when at least 150 events reached

Historical data:

- promising internal PoC
- promising phase II

Co-data:

- two phase III trials run in parallel $\Rightarrow$ each phase III
  trial is *concurrent* with the other

Define design choices
```{r}
Nev <- 379

alt_HR <- 0.75
alt_logHR <- log(alt_HR)

alpha <- 0.025
```

## Power calculation


Here we use the unit information prior as non-informative prior and
define it using the mean & effective sample size (ESS) specification:

```{r}
unit_inf <- mixnorm(c(1, 0, 1), sigma=2, param="mn")
unit_inf
```

Define conditional power for the overall trial:
```{r}
success_crit <- decision1S(1-alpha, 0)
## let's print the defined criterion
success_crit
design <- oc1S(unit_inf, Nev, success_crit, sigma=2)
```

Under the alternative these design choices result in 80% power

```{r}
design(alt_logHR)
```

The impact of the unit-information prior is minimal which can be seen
by comparing to the frequentist calculation:

```{r}
power.t.test(n=Nev, delta=-1*alt_logHR, sd=2, type="one.sample", sig.level=0.025, alternative="one.sided")
```

With RBesT we can explore the conditional power for a range of
alternatives:

```{r}
ggplot(data.frame(HR=c(0.5, 1.2)), aes(HR)) +
    stat_function(fun=compose(design, log)) +
    vline_at(c(alt_HR, 1.0), linetype=I(2)) +
    scale_y_continuous(breaks=seq(0,1,by=0.1)) +
    scale_x_continuous(breaks=c(alt_HR, seq(0,1.2,by=0.2))) +
    ylab(NULL) + xlab(expression(theta[a])) +
    ggtitle(paste("Power for N =", Nev, "events"))
```


## Critical value

The critical value determines at which observed logHR we *just*
conclude that the success criterion is fulfilled.

```{r}
design_crit <- decision1S_boundary(unit_inf, Nev, success_crit, sigma=2)

design_crit
exp(design_crit)
```

We can check this:

```{r}
success_crit(postmix(unit_inf, m=design_crit, n=379))
```

Ok, when observing the critical value, we get a success.

Now, what if we observe a 1% worse result?

```{r}
success_crit(postmix(unit_inf, m=design_crit+log(1.01), n=379))
```

No success then $\Rightarrow$ this is the critical boundary value.

# PoS at interim for phase III trial A only

## No use of historical information

Posterior of treatment effect at interim. The trial uses a
non-informative prior for the treatment effect:

```{r}
interim_A <- postmix(unit_inf, m=trials$logHR[3], se=trials$sem[3])
interim_A
```

Now we are interested in the PoS at trial completion. The prior to
use for the analysis of the second half is given by the data
collected so far.

```{r}
interim_pos_A <- pos1S(interim_A, Nev-trials$deaths[3], success_crit, sigma=2)
```

The returned function can now calculate the PoS assuming any
distribution on the treatment effect. In case we do not use any
historical information, then this is just the interim posterior:

```{r}
interim_pos_A(interim_A)
```

The above command integrates the *conditional power* over the
uncertainty which we have about the treatment effect as defined above
for $PoS_I$.

The conditional power and the operating characteristics of a trial
coincide whenever we do not condition on any observed data. The key
difference of the conditional power as compared to the probability of
success is that it assumes a known value for the parameter of
interest. This can be seen as follows: First define the conditional
power which is conditional on the observed data,
$CP_{N-n_I}(\theta|y_{n_I})$:

```{r}
interim_oc_A <- oc1S(interim_A, Nev-trials$deaths[3], success_crit, sigma=2)
```

The conditional power assuming the alternative is true (a HR of 0.75):

```{r}
interim_oc_A(alt_logHR)
```

In case there is no uncertainty of the treatment effect (here
$se=10^-4$), then this result agrees with the probability of success
calculation:

```{r}
interim_pos_A(mixnorm(c(1, alt_logHR, 1E-4)))
```

For trial B the calculation is:
```{r}    
interim_B <- postmix(unit_inf, m=trials$logHR[4], se=trials$sem[4])
interim_pos_B <- pos1S(interim_B, Nev-trials$deaths[4], success_crit, sigma=2)
interim_pos_B(interim_B)
```

## Use of historical information - MAP approach

However, we have historical information of which we can take advantage
at the interim for a better informed decision.

Our data before the phase III trials includes the PoC and the phase II
trial. We now derive from these a MAP prior; recall that the MAP prior
is the prediction of the log-hazard ratio of a future trial:

```{r}
base <- trials[1:2,]

set.seed(342345)
base_map_mc <- gMAP(cbind(logHR, sem) ~ 1 | study,
                    family=gaussian,
                    data=base,
                    weights=deaths,
                    tau.dist="HalfNormal", tau.prior=0.5,
                    beta.prior=cbind(0, 2))

forest_plot(base_map_mc, est="MAP")

base_map <- automixfit(base_map_mc)

plot(base_map)$mix + xlab(expression(log(theta)))

base_map
```

At the interim we have even more knowledge available on the treatment
effect through the interim data itself which we can include into the
MAP prior:

```{r}    
interim_A_combined <- postmix(base_map, m=trials$logHR[3], se=trials$sem[3])
```

The PoS for this posterior at interim (representing historical
*and* interim data collected) is:
```{r}
interim_pos_A(interim_A_combined)
```

Note that we have not redefined `interim_pos_A`, such that this
calculates the PoS for the phase III A trial taking into account that
the final analysis will use a non-informative prior.

For trial B the calculation is:
```{r}    
interim_B_combined <- postmix(base_map, m=trials$logHR[4], se=trials$sem[4])
interim_pos_B(interim_B_combined)
```

## Use of historical information - MAC approach

However, there is even more information which can be used here,
since the phase III result of trial B is also available:

```{r}
interim_map_mc <- update(base_map_mc, data=trials)
```

Now the trial B specific posterior at interim is
```{r}
kable(fitted(interim_map_mc), digits=3)
```

which we can extract as:


1. obtain posterior (which we restrict to the first 4 columns)

```{r}
interim_map_post <- as.matrix(interim_map_mc)[,1:4]

dim(interim_map_post) # posterior is given as matrix: iteration x parameter

head(interim_map_post, n=3)
```

2. turn MCMC posterior sample into parametric mixture

```{r}
interim_A_allcombined <- automixfit(interim_map_post[,"theta[3]"])
```

3. and finally evaluate the PoS

```{r}
interim_pos_A(interim_A_allcombined)
```

which aligns with the published result under the assumption of full
exchangeability.


For trial B computations are:

```{r}
interim_B_allcombined <- automixfit(interim_map_post[,"theta[4]"])
interim_pos_B(interim_B_allcombined)
```

# Differential discounting

Differential discounting allows to weight different data-sources
differently. For example, we may assume greater heterogeneity for
the historical data in comparison to the twin phase III trials.

Assign data to historical (2) and concurrent data strata (1):

```{r}
trials <- trials %>% mutate(stratum=c(2, 2, 1, 1))

kable(trials, digits=2)

set.seed(435345)
interim_diff_map_mc <- gMAP(cbind(logHR, sem) ~ 1 | study,
                            tau.strata=stratum,
                            family=gaussian,
                            data=trials,
                            weights=deaths,
                            tau.dist="HalfNormal", tau.prior=c(0.5, 1),
                            beta.prior=cbind(0, 2))

interim_diff_map_post <- as.matrix(interim_diff_map_mc)[,1:4]

interim_A_diff_allcombined <- automixfit(interim_diff_map_post[,"theta[3]"])
interim_B_diff_allcombined <- automixfit(interim_diff_map_post[,"theta[4]"])

interim_pos_A(interim_A_diff_allcombined)
interim_pos_B(interim_B_diff_allcombined)
```



# PoS for both phase III trials being successful

So far we have only calculated the individual PoS per trial, but
more interesting is the *overall* PoS for both trials being
successful.

Recall, the PoS is the conditional power integrated over an assumed
true effect distribution. Hence, we had for trial A:

```{r}
interim_pos_A(interim_A)
```

As explained, the conditional power is the operating characerstic
of a design when conditioning on the already observed data:
```{r}
interim_oc_A <- oc1S(interim_A, Nev-trials$deaths[3], success_crit, sigma=2)
```

The PoS is then the integral of the conditional power over the
parameter space $\theta$ representing our knowledge. This integral
can be evaluated in a Monte-Carlo (MC) approach as

$$ PoS_I = \int CP_{N-n_I}(\theta|y_{n_I}) \, p(\theta|y_{n_I}) \,
d\theta \approx \frac{1}{S} \sum_{i=1}^S CP(\theta_i),$$

whenever we have a sample of $p(\theta|y_{n_I})$ of size $S$... which
we have:

```{r}
interim_A_samp <- rmix(interim_A, 1E4)
mean(interim_oc_A(interim_A_samp))
```

This is an MC approach to calculating the PoS.

When now considering the probability for both trials being
successful we have to perform an MC integration over the joint
density $p(\theta_A,\theta_B|y_{n_{I_A}},y_{n_{I_B}})$

$$
\begin{aligned}
PoS &= \iint CP_{N-n_{I_A}}(\theta_A|y_{n_{I_A}}) \,
CP_{N-n_{I_B}}(\theta_B|y_{n_{I_B}})\,
p(\theta_A,\theta_B|y_{n_{I_A}},y_{n_{I_B}}) \, d\theta_A d\theta_B \\
& \approx
\frac{1}{S} \sum_{i=1}^S CP_{N-n_{I_A}}(\theta_{A,i}|y_{n_{I_A}}) \,
CP_{N-n_{I_B}}(\theta_{B,i}|y_{n_{I_B}}).
\end{aligned}
$$

Thus we need to also get the conditional power for trial B at interim...

```{r}
interim_oc_B <- oc1S(interim_B, Nev-trials$deaths[4], success_crit, sigma=2)
```

...and integrate over the posterior samples (differential discounting case)

```{r}
mean(interim_oc_A(interim_diff_map_post[,"theta[3]"]) * interim_oc_B(interim_diff_map_post[,"theta[4]"]))
```

which is slightly larger than assuming independence:
```{r}
interim_pos_A(interim_A) * interim_pos_B(interim_B)
```

This is due to dependence of the posteriors
```{r}
cor(interim_diff_map_post[,c("theta[3]", "theta[4]")])
```

For the full exchangeability case we have
```{r}
mean(interim_oc_A(interim_map_post[,"theta[3]"]) * interim_oc_B(interim_map_post[,"theta[4]"]))
```


# Summary

We have now calculated with increasing complexity the probability of
success for various data constellations. As new trials are only
conducted whenever previous trial results were positive, it is
important to take note of the potential selection bias. Moreover,
adding more historical data sources in this situation will likely
increase the probability of success as illustrated by this summary of
our preceding calculations.

Phase III trial A:
```{r}
## only interim data of trial A 
interim_pos_A(interim_A)
## in addition with prior historical data PoC & phase II data
interim_pos_A(interim_A_combined)
## finally with the interim data of the phase III B
interim_pos_A(interim_A_allcombined)
```

Phase III trial B:
```{r}
## only interim data of trial B
interim_pos_B(interim_B)
## in addition with prior historical data PoC & phase II data
interim_pos_B(interim_B_combined)
## finally with the interim data of the phase III A
interim_pos_B(interim_B_allcombined)
```

# Appendix: MAP and MAC equivalence

In the preceeding sections we have used MAP and MAC equivalence
already. The proof for the equivalence is presented reference in
[2]. The formal deriavtion is shown at the end of this section

While MAP and MAC provide the exact same results, the difference is a
sequential vs a joint analysis as (see also [2]):

1. MAP: Summarize historical information as MAP and then update the
   MAP with the trial result (MCMC, then `postmix`).
2. MAC: Directly summarize historical information and trial result in
   a single step (only MCMC on all data).

The two results above using MAP and MAC did not line up. The reason
here is that the MAP approach used the historical data and the phase
III trial A interim data only. In contrast, the MAC approach used the
historical data and interim phase III data of both trials. To show the
equivalence we need to align this mismatch of used data.

Run `gMAP` with base data and produce a large MCMC sample (10
chains) to get a very high precision.

```{r}
base_map_mc_2 <- gMAP(cbind(logHR, sem) ~ 1 | study,
                      family=gaussian,
                      data=base,
                      weights=deaths,
                      tau.dist="HalfNormal", tau.prior=0.5,
                      beta.prior=cbind(0, 2),
                      chains=ifelse(is_CRAN, 2, 20))
```


Force an accurate fit with 5 components:

```{r}
base_map_2 <- mixfit(base_map_mc_2, Nc=5)
base_map_2
```

Now, combine the MAP prior (representing historical knowledge) with
the interim data of trial A:

```{r}
interim_A_combined_2 <- postmix(base_map_2, m=trials$logHR[3], se=trials$sem[3])
```


1. Run the respective MAC analysis (thus we need historical data +
   phase III A trial, but excluding the phase III B data):

```{r}
interim_map_mc_2 <- update(base_map_mc_2, data=trials[-4,])

interim_map_post_2 <- as.matrix(interim_map_mc_2)[,1:3]
```


2. turn MCMC posterior sample into parametric mixture
```{r}
interim_A_allcombined_2 <- mixfit(interim_map_post_2[,"theta[3]"], Nc=5)

interim_A_allcombined_2
```


Now let's overlay the two posterior's

```{r}
ggplot(data.frame(logHR=c(-0.8,0.25)), aes(logHR)) +
    stat_function(fun=dmix, args=list(mix=interim_A_combined_2), aes(linetype="MAP")) +
    stat_function(fun=dmix, args=list(mix=interim_A_allcombined_2), aes(linetype="MAC")) +
    scale_linetype_discrete("Analysis") +
    ggtitle("Posterior log hazard of phase III A trial at interim")
```


The PoS is essentially the same

```{r}
interim_pos_A(interim_A_combined_2)
interim_pos_A(interim_A_allcombined_2)
```

## Formal MAP and MAC equivalence

The stated equivalence requires that the posterior of a trial
specific parameter

$$p(\theta_\star|y_\star,y_H),$$

which is conditional on the trial specific data $y_\star$ **and** the
historical data $y_H$ (MAC approach, joint use of $y_H,y_\star$), is
equivalent to obtaining the MAP prior $p(\theta_\star|y_H)$ based on
the historical data and then analyzing the new trial with this prior.

$$
\begin{aligned}
p(\theta_\star|y_\star,y_H) &\propto p(\theta_\star,\theta_H|y_\star,y_H) \\
 &\propto p(y_\star,y_H|\theta_\star,\theta_H) \,
 p(\theta_\star,\theta_H) \\ 
 &= p(y_\star|\theta_\star) \, p(y_H|\theta_H) \,
 p(\theta_\star,\theta_H) \\ 
&\propto p(y_\star|\theta_\star) \, p(\theta_\star,\theta_H|y_H) \\
&\propto p(y_\star|\theta_\star) \, p(\theta_\star|y_H)
\end{aligned}
$$

The equivalence holds under the use of the meta-analytic model.

# References

[1] Neuenschwander, B., Roychoudhury, S., & Schmidli, H. (2016). On the Use of Co-Data in Clinical Trials. Statistics in Biopharmaceutical Research, 8(3), 345-354. \doi{10.1080/19466315.2016.1174149}

[2] 1. Schmidli H, Gsteiger S, Roychoudhury S, O'Hagan A, Spiegelhalter D, Neuenschwander B. Robust meta-analytic-predictive priors in clinical trials with historical control information. Biometrics. 2014;70(4):1023-1032. \doi{10.1111/biom.12242}

# R Session Info

```{r}
sessionInfo()
```

```{r,include=FALSE}
options(.user_mc_options)
```
