---
title: "Introduction to the 'tipmap' package"
subtitle: "Tipping point analysis for clinical trials that employ Bayesian dynamic borrowing" 
package: tipmap
author: "Christian Stock"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_document:
          theme: "default"
          highlight: "default"
          toc: true
          toc_float: true
bibliography: references.bib
csl: jrss_style.csl
#bibliography: '`r system.file("references.bib", package="tipmap")`'
#csl: '`r system.file("jrss_style.csl", package="tipmap")`'
vignette: >
  %\VignetteIndexEntry{Introduction to the 'tipmap' package}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options:
  chunk_output_type: console
---


```{r setup, include = F}
knitr::opts_chunk$set(
  echo = T, collapse = T, warning = F, message = F, 
  prompt = T, comment = "#",
  out.width = "100%"
)
library(ggplot2)
ggplot2::theme_set(theme_bw())
```


## Purpose of the package

The R package **tipmap** implements a tipping point analysis for clinical trials that employ Bayesian dynamic borrowing of a treatment effect from external evidence via robust meta-analytic predictive (MAP) priors. A tipping point analysis allows to assess how much weight on the informative component of a robust MAP prior is needed to conclude that the investigated treatment is efficacious, based on the total evidence. The package mainly provides an implementation of a graphical approach proposed by @Best2021 for different one-sided evidence levels (80%, 90%, 95%, 97.5%).

Tipping point analyses can be useful both at the planning and the analysis stage of a clinical trial that uses external information. At the _planning stage_, they can help to determe (pre-specify) a weight of the informative component of the MAP prior for a primary analysis. 
Various possible results of the planned trial in the target population and implications for the treatment effect estimate and statistical inferences based on the total evidence may be explored for a range of weights. 
Through this exercise, in addition to other criteria, decision-makers can develop a sense under which circumstances they would still feel comfortable to establish efficacy with a specific level of certainty. A preferred primary weight will typically be a compromise between the belief in the applicability of the data and operating characteristics of the resulting design specifications.
At the _analysis stage_, tipping point analyses can be used as a sensitivity analysis to assess the dependency of the treatment effect estimate and statistical inferences on the weight of the informative component of the MAP prior. This can also be understood in the sense of a reverse-Bayes analysis (@Held2022).

This vignette shows an exemplary application of the tipping point analysis with hypothetical data.

Further functions in this package (not illustrated in this vignette) facilitate the specification of a robust MAP prior via expert elicitation, specifically the choice of a primary weight (using the roulette method).

Intended use of the **tipmap**-package is the planning, analysis and interpretation of (small) clinical trials in pediatric drug development, where extrapolation of efficacy, often through Bayesian methods, has become increasingly common (@Gamalo2022; @ICHE11A; @Ionan2023; @Travis2023). For an application, see @Maher2024.

For the implementation of the MAP prior approach, including computation of the posterior distribution, the **RBesT**-package is used (@Weber2021).


## Generating data for an exemplary tipping point analysis

In this vignette, we assume that results from three clinical trials conducted in adult patients (the _source_ population) are available, which share key features with a new trial among pediatric patients (the _target_ population).
For example, they had been conducted in the same indication, studied the same drug and provided results on an endpoint of interest for the target population. This means a certain degree of _exchangeability_ between the trials in the source and target population can be assumed. 
The similarity in disease and response to treatment between source and target population always need to be carefully considered in any setting, usually by clinical experts in the disease area. 

We assume that it is supported by medical evidence and now planned to consider the trials in adult patients in a Bayesian dynamic borrowing approach, and we would like to create a robust MAP prior (@Schmidli2014). The treatment effect measure of interest is assumed to be a mean difference between a treated group and a control group with respect to a continuous endpoint.

### Derivation of MAP prior based on trials in the source population

We start by specifying an object that contains the prior data.

The function `create_prior_data()` takes vectors of total sample sizes, treatment effect estimates and their standard errors as arguments and generates a data frame. A study label is optional.
```{r, eval=T, echo=T}
library(tipmap)
prior_data <- create_prior_data(
  n_total = c(160, 240, 320),
  est = c(1.16, 1.43, 1.59),
  se = c(0.46, 0.35, 0.28)
)
```

```{r, eval=F, echo=F}
# compute standard deviation of change in each arm (assumed equal);
# for two-sample data:
sd <- prior_data$se / sqrt(1/(prior_data$n_total/2) + 1/(prior_data$n_total/2)) 
sigma1 <- mean(sd)
sigma1 # = 2.70826
nt <- 15; nc <- 15
f <- (nt+nc)^2 / (nt*nc)
sigma2 <- sqrt(f)*sigma1
sigma2 # = 5.41652
```

```{r, eval=T, echo=T}
print(prior_data)
```

We then generate a MAP prior from our prior data using the **RBesT**-package (@Weber2021).
```{r, eval=T, echo=T}
set.seed(123)
uisd <- 5.42
map_mcmc <- RBesT::gMAP(
  formula = cbind(est, se) ~ 1 | study_label,
  data = prior_data,
  family = gaussian,
  weights = n_total,
  tau.dist = "HalfNormal",
  tau.prior = cbind(0, uisd / 16),
  beta.prior = cbind(0, uisd)
  )
```

A few additional specifications are needed to be made to fit the MAP prior model; for details see @Neuenschwander2020 or @Weber2021. The variable `uisd` here represents an assumed unit-information standard deviation and the specification of the prior on between-trial heterogeneity parameter tau follows recommendations to consider moderate heterogeneity for a two-group parameter, such as the mean difference (@Neuenschwander2020).

This is a summary of the fitted model based on samples from the posterior distribution:
```{r, eval=T, echo=T}
summary(map_mcmc)
```

A forest plot of the Bayesian meta-analysis is shown in Figure 1. It is augmented with meta-analytic shrinkage estimates per trial. The figure shows the per-trial point estimates
(light point) and the 95% frequentist confidence intervals (dashed line) and the model derived
median (dark point) and the 95% credible interval of the meta-analytic model. 

```{r forest_plot, eval=T, echo=T, fig.width=6, fig.height=3, dev=c('png'), out.width="70%", fig.cap='Figure 1: Forest plot.'}
plot(map_mcmc)$forest_model
```

Subsequently, the MAP prior is approximated by a mixture of conjugate normal distributions. The parametric form facilitates the computation of posteriors when the MAP prior is combined with results from the trial in the target population.
```{r, eval=T, echo=T}
map_prior <- RBesT::automixfit(
  sample = map_mcmc,
  Nc = seq(1, 4),
  k = 6,
  thresh = -Inf
  )
```

The approximation yields a mixture of two normals:
```{r, eval=T, echo=T}
print(map_prior)
```

The density of the parametric mixture together with a histogram of MCMC samples from the `map_mcmc` object is shown in Figure 2.

```{r map_prior_dens, eval=T, echo=T, fig.width=6, fig.height=3, dev=c('png'), out.width="70%", fig.cap='Figure 2: Overlay of the MCMC histogram of the MAP prior and the fitted parametric mixture approximation.'}
plot(map_prior)$mix
```

The derivation of the MAP prior is now complete. For normal likelihoods the parametric representation by a mixture of normals can be used to calculate posterior distributions analytically.


### Trial results in the target population

We now create a numeric vector with data on pediatric trial (the total sample size, the treatment effect estimate and its standard error). In the planning stage, this may be an expected result.

```{r, eval=T, echo=T}
pediatric_trial <- create_new_trial_data(n_total = 30, est = 1.02, se = 1.4)
```

```{r, eval=T, echo=T}
print(pediatric_trial)
```

The function `create_new_trial_data()` computes quantiles, assuming normally distributed errors. This is merely used to plot a confidence interval for the treatment effect estimate obtained in the target trial in the tipping point plot.


## Performing the tipping point analysis

### Computation of posteriors for a range of weights

We can now compute posterior distributions for a range of weights
using the function `create_posterior_data()`.
```{r, eval=T, echo=T}
posterior <- create_posterior_data(
  map_prior = map_prior,
  new_trial_data = pediatric_trial,
  sigma = uisd)
```

```{r, eval=T, echo=T}
head(posterior, 4)
```
```{r, eval=T, echo=T}
tail(posterior, 4)
```

The resulting data frame has `r dim(posterior)[1]` rows and `r dim(posterior)[2]` columns. 
The weights increase incrementally in steps of 0.005 from 0 to 1, i.e.\ posterior quantiles for `r length(posterior$weight)` weights are computed. For each weight the data frame contains the following `r dim(posterior)[2]-1` posterior quantiles.
```{r, eval=F, echo=F}
length(posterior$weight)
dim(posterior)[1]
class(posterior)
colnames(posterior)[-1]
```

```{r, eval=T, echo=T}
colnames(posterior)[-1]
```
These posterior quantiles can be directly used for inferences based on the total evidence (new data and prior combined). They reflect one-sided 99%, 97.5%, 95%, 90%, 80%, and 50% evidence levels for a given weight, respectively.


### Creating the tipping point plot

The function to produce the tipping point plot is called `tipmap_plot()`, it requires a dataframe with data on all components generated by the function `create_tipmap_data()`.

```{r, eval=T, echo=T}
tipmap_data <- create_tipmap_data(
  new_trial_data = pediatric_trial,
  posterior = posterior,
  map_prior = map_prior)
```

```{r tipmap_plot, eval=T, echo=T, fig.width=8, fig.height=5, dev=c('png'), out.width="95%", fig.cap='Figure 3: Tipping point plot.'}
(p1 <- tipmap_plot(tipmap_data = tipmap_data))
```

At the center of the plot, quantiles of the posterior distribution are shown, defining
credible intervals of the effect estimate and also one-sided evidence-levels for given weights 
of the informative component of the MAP prior. 
The intersections between the lines connecting the respective quantiles and the horizontal line at 0 (the null effect) are referred to as tipping points (indicated by vertical lines in red color). They indictae the minimum weight that is required to conclude that the treatment is efficacious for a given one-sided evidence level (@Best2021).
On the left and right side of the plot, the treatment effect estimate obtained in the trial in the (pediatric) target population (with 95% confidence interval) and the MAP prior (with 95% credible interval) are shown, respectively.

The plot is a `ggplot`-object that can be modified accordingly.
For example, if we had chosen a primary weight of 0.38, we could add a vertical reference line at this position. There are additional features to customize the plot in the `tipmap_plot()` function, see `help(tipmap_plot)`.

```{r tipmap_plot_refline, eval=T, echo=T, fig.width=8, fig.height=5, dev=c('png'), out.width="95%", fig.cap='Figure 4: Tipping point plot with reference line.'}
primary_weight <- 0.38
(p2 <- p1 + ggplot2::geom_vline(xintercept = primary_weight, col="green4"))
```

We see from Figure 4 that, for a weight of 0.38, there is a probability of larger than 90% but less than 95% based on the posterior distribution that the treatment effect is larger than 0, i.e. the treatment is efficacious.

### Extracting quantities of interest

The data frame with posteriors for all weights can be filtered to obtain posterior quantiles for weights of specific interest by the function `get_posterior_by_weight()`:

```{r, eval=T, echo=T}
get_posterior_by_weight(
  posterior = posterior, 
  weight = c(primary_weight)
  )
```
The function `get_tipping_points()` extracts tipping points for one-sided 80%, 90%, 95% and 97.5% evidence levels, respectively.

```{r, eval=T, echo=T}
tipp_points <- get_tipping_points(
  tipmap_data = tipmap_data,  
  quantile = c(0.2, 0.1, 0.05, 0.025)
)
tipp_points
```

Calculating the precise posterior probability that treatment effect exceeds a threshold value is possible via functions in the **RBesT**-package.

```{r, eval=T, echo=T}
prior_primary <- RBesT::robustify(
  priormix = map_prior,
  weight = (1 - primary_weight),
  m = 0,
  n = 1,
  sigma = uisd
  )
```

```{r, eval=T, echo=T}
posterior_primary <- RBesT::postmix(
  priormix = prior_primary,
  m = pediatric_trial["mean"],
  se = pediatric_trial["se"]
  )
```

```{r, eval=F, echo=F}
summary(posterior_primary)
```

The posterior probability that the treatment effect is larger than 0, 0.5 and 1, respectively, can be assessed through the cumulative distribution function of the posterior.

```{r, eval=T, echo=T}
round(1 - RBesT::pmix(posterior_primary, q = 0), 3)
round(1 - RBesT::pmix(posterior_primary, q = 0.5), 3)
round(1 - RBesT::pmix(posterior_primary, q = 1), 3)
```

This is illustrated by a cumulative density curve of the posterior.

```{r cumulative_dens, eval=T, echo=T, fig.width=7, fig.height=4.5, dev=c('png'), out.width="80%", fig.cap='Figure 5: Cumulative density of posterior with weight w=0.38.'}
library(ggplot2)
plot(posterior_primary, fun = RBesT::pmix) +
  scale_x_continuous(breaks = seq(-1, 2, 0.5)) +
  scale_y_continuous(breaks = 1-c(1, 0.927, 0.879, 0.782, 0.5, 0),
                     limits = c(0,1),
                     expand = c(0,0)
                     ) +
  ylab("Cumulative density of posterior with w=0.38") +
  xlab("Quantile") +
  geom_segment(aes(x = 0,
                   y = RBesT::pmix(mix = posterior_primary, q = 0), 
                   xend = 0, 
                   yend = 1), 
               col="red") +
  geom_segment(aes(x = 0.5,
                   y = RBesT::pmix(mix = posterior_primary, q = 0.5), 
                   xend = 0.5, 
                   yend = 1), 
               col="red") + 
  geom_segment(aes(x = 1,
                   y = RBesT::pmix(mix = posterior_primary, q = 1), 
                   xend = 1, 
                   yend = 1), 
               col="red") + 
  theme_bw()
```


As a further example, for the weight corresponding to the tipping point of the one-sided evidence-level of 95% (=0.51), we would obtain a posterior probability of 95% that the treatment effect is larger than 0.

```{r, eval=T, echo=T}
tipp_points[3]
```

```{r, eval=T, echo=T}
prior_95p <- RBesT::robustify(
  priormix = map_prior,
  weight = (1 - tipp_points[3]),
  m = 0,
  n = 1,
  sigma = uisd
  )
```

```{r, eval=T, echo=T}
posterior_95p <- RBesT::postmix(
  priormix = prior_95p,
  m = pediatric_trial["mean"],
  se = pediatric_trial["se"]
  )
```

```{r, eval=T, echo=T}
round(1 - RBesT::pmix(posterior_95p, q = 0), 3)
```


## References