---
title: "Visualisation with CRM"
author: "Kristian Brock"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Visualisation with CRM}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

In the general [introductory CRM vignette](CRM.html), we introduced the different flavours of the Continual Reassessment Method (CRM) implmented in `trialr`.
In this vignette, we demonstrate some visualisation methods that will be useful for conveying inferences to investigators and the wider research community.

The following visualisation methods are applicable in all `trialr` variants of the CRM because they rely only on the family of `prob_tox` parameters that are estimated by all CRM models.

For illustration, let us assume that we have treated 6 patients at 3 dose-levels:

| Patient | Cohort | Dose-level | DLT |
|:-------:|:------:|:----------:|:---:|
|    1    |   1    |      2     |  0  |
|    2    |   1    |      2     |  0  |
|    3    |   2    |      3     |  0  |
|    4    |   2    |      3     |  0  |
|    5    |   3    |      4     |  1  |
|    6    |   3    |      4     |  1  |

```{r}
outcomes <- '2NN 3NN 4TT'
```

and that we are using the empiric CRM with the following prior to target the dose with Pr(DLT) closest to 25%:

```{r}
skeleton <- c(0.05, 0.15, 0.25, 0.4, 0.6)
target <- 0.25
```

To access the Stan implementations of the CRM, we must load `trialr`:
```{r, message=FALSE, warning=FALSE}
library(trialr)
```

and to fit the data to the model, we run:
```{r, results = "hide", warning=FALSE, message=FALSE}
fit <- stan_crm(outcomes, skeleton = skeleton, target = target, 
                beta_sd = sqrt(1.34), seed = 123)
fit
```

Recall that we set the random number generator seed in demonstrations for reproducibility.
You may or may not choose to do this in a real analysis.

We saw [before](CRM.html) that dose-level 2 is closest to our target toxicity rate of 25%.
However, we will probably like to convey this information to the research world in a visually-appealing manner.
This is simple with access to the posterior samples.

For ease of plotting with `ggplot2`, we recompose the posterior samples to a tall "tidy" format using `tidyr` and `dplyr`:

```{r, warning=FALSE, message=FALSE}
library(tidyr)
library(dplyr)

prob_tox_samp <- as.data.frame(fit, 'prob_tox')
prob_tox_samp_tall <- prob_tox_samp %>%
  gather(Label, ProbTox) %>%
  mutate(
    DoseLevel = rep(1:ncol(prob_tox_samp), each = nrow(prob_tox_samp)),
    Draw = rep(1:nrow(prob_tox_samp), times = ncol(prob_tox_samp))
  )
```

The tidy data looks like this:

```{r}
prob_tox_samp_tall %>% head(10)
```

Boxplots would be a traditional way of visualising the distributions of the probability of toxicity at each dose:

```{r, fig.width=7, fig.height=7}
library(ggplot2)

prob_tox_samp_tall %>% 
  ggplot(aes(x = DoseLevel, y = ProbTox, group = DoseLevel)) +
  geom_boxplot() + 
  ylim(0, 1) + 
  labs(title = 'boxplot of Pr(DLT) under CRM')
```

However, boxplots give only limited information on the distributions.
For instance, it might be tempting to assume that the probability of toxicity is normally distributed at each dose-level.
The boxplots suggest some wide tails.
This inference is much more clear, however, using a violin-plot:

```{r, fig.width=7, fig.height=7}
prob_tox_samp_tall %>% 
  ggplot(aes(x = DoseLevel, y = ProbTox, group = DoseLevel)) +
  geom_violin(fill = 'orange') + 
  ylim(0, 1) + 
  labs(title = 'violin plot of Pr(DLT) under CRM')
```

If you are a fan of post-punk UK music (and you have installed the [`ggridges`](https://CRAN.R-project.org/package=ggridges) package), you may however prefer to show this information using a ridge plot, aka a [joyplot](https://en.wikipedia.org/wiki/Unknown_Pleasures)

```{r, fig.width=7, fig.height=7, message=FALSE}
library(ggridges)

prob_tox_samp_tall %>% 
  mutate(DoseLevel = factor(DoseLevel)) %>% 
  ggplot(aes(x = ProbTox, y = DoseLevel, fill = DoseLevel)) +
  geom_density_ridges() + 
  theme(legend.position = 'none') +
  labs(title = 'joyplot of Pr(DLT) under CRM') + 
  theme(legend.position = 'bottom')
```

Hopefully none of us would try to claim these posterior probabilities of toxicity are normally distributed under this model.
Assuming normality has been one method for performing posterior inference with CRM models in the non-MCMC setting.
With the posterior samples provided by `rstan`, we do not need to assume.

We will naturally want to visualise quantities beyond just the probability of toxicity.
We learned in the [introductory CRM vignette](CRM.html) that with the full Bayesian CRM provided by `trialr` and `rstan`, we can calculate the probability that each dose is the maximum tolerable dose (MTD).

We can visualise the MCMC candidates for the dose-toxicity curve on one plot.
Colouring them by the MTD candidate they propose (i.e. using a single colour for all the curves that suggest dose-level 1 is the maximum tolerable dose, etc), we get an idea of uncertainty still in this trial:

```{r, fig.width=7, fig.height=7, message=FALSE}
prob_tox_samp_tall %>% 
  filter(Draw <= 200) -> prob_tox_subset

prob_tox_subset %>% 
  left_join(
    prob_tox_subset %>% 
      group_by(Draw) %>% 
      summarise(MTD = which.min(abs(ProbTox - target)) %>% factor()),
    by = 'Draw') %>% 
  ggplot(aes(x = DoseLevel, y = ProbTox, group = Draw)) +
  geom_line(aes(col = MTD), alpha = 0.5) + 
  geom_hline(yintercept = target, col = 'red', linetype = 'dashed') + 
  labs(title = 'Identify of MTD is shrouded in mystery') +
  theme(legend.position = 'bottom')
```

We used just 200 curves above to avoid saturating the plot.
We can visualise that data rather more bluntly:

```{r, fig.width=7, fig.height=7, message=FALSE, warning=FALSE}
apply(prob_tox_samp, 1, function(x) which.min(abs(x - target))) %>% 
  data.frame(MTD = factor(.)) %>% 
  count(MTD) %>% 
  mutate(ProbMTD = n / sum(n)) %>% 
  ggplot(aes(x = MTD, y = ProbMTD, fill = MTD)) + 
  geom_col() +
  theme(legend.position = 'bottom')
```

In this interim stage, each of the first four doses could plausibly be the MTD, but the top dose looks unlikely.
This information was not readily available from some of the above plots of the probabilities of toxicity.

We might also like to visualise the probability that the toxicity rate at each dose exceeds our target toxicity rate.

```{r, fig.width=7, fig.height=7}
colMeans(prob_tox_samp > target) %>% 
  data.frame(ProbToxic = ., DoseLevel = 1:length(skeleton)) %>%
  ggplot(aes(x = DoseLevel, y = ProbToxic, fill = ProbToxic)) + 
  geom_col() + 
  scale_fill_gradient(low="green", high="red") + 
  labs(title = 'Posterior probability that each dose is too toxic',
       y = 'Prob(Tox > target)', fill = 'Prob too toxic') +
  theme(legend.position = 'bottom')
```

Based on our prior and the data assembled thus far, dose-levels 4 and 5 look quite likely to be overdoses.

## Other CRM vignettes
There are many vignettes illuminating the CRM in `trialr`:

- [An introduction to CRM in 'trialr'](CRM.html)
- [Visualisation in CRM](CRM-visualisation.html)
- [TITE-CRM](TITE-CRM.html)
- Stopping in CRM
- Case study using Levy, et al. (2006)
- Simulation


# trialr
`trialr` is available at https://github.com/brockk/trialr and https://CRAN.R-project.org/package=trialr 


# References
