---
title: "Advanced example"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Advanced example}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6, 
  fig.align = "center"
)
```

This vignette provides an **advanced** example on how to:

1. Specify a trial design with `setup_trial()` using user-written functions for
generating outcomes and drawing samples from the posterior distributions
2. Transform the raw outcome estimates in a way appropriate for your needs
3. Derive and specify prior probability distributions

The manual for `setup_trial()` contains another example (more succinctly
annotated); see the **Examples** sections in `?setup_trial()`. As an
alternative to using custom functions, the convenience functions
`setup_trial_binom()` and `setup_trial_norm()` may be used for simulating trials
with binary, binomially distributed or continuous, normally distributed
outcomes, using sensible defaults and flat priors.

This vignette is intended for advanced users and assume general knowledge of
writing `R` functions and specifying prior distributions for (simple, conjugate)
Bayesian models. Before reading this vignette, we highly encourage you to read
the Basic examples vignette (see `vignette("Basic-examples", "adaptr")`) with
more basic examples using the simplified `setup_trial_binom()` function and a
thorough discussion of the settings applicable to all trial designs.

## Preamble

In this example, we set up a trial with three `arms`, one of which is a common
`control`, and an undesirable binary outcome (e.g., mortality).

This examples creates a custom version of the `setup_trial_binom()` function
using non-flat priors for the event rates in each arm (`setup_trial_binom()`
uses flat priors), and returning event probabilities as percentages (instead of
fractions), to illustrate the use of a custom function to summarise the raw
outcome data.

While `setup_trial()` attempts to validate the custom functions by
assessing their output during trial specification, some edge cases might elude
this validation. We, therefore, urge users specifying custom functions to
carefully test more complex functions themselves before actual use. 

If you go through the trouble of writing a nice set of functions for generating
outcomes and sampling from posterior distributions, please consider adding these
to the package. This way, others can benefit from your work and it helps
validate them. See how on the
[GitHub](https://github.com/INCEPTdk/adaptr#contributing) page under
**Contributing**.

Although the user-written custom functions below do not depend on the `adaptr`
package, as the first thing we load the package:

```{r setup}
library(adaptr)
```

--and then set the global seed to ensure reproducible results in the vignette:
```{r}
set.seed(89)
```

We define the functions and (for illustration purposes and as a sanity check)
print their outputs. They are, then, invoked by `setup_trial()` (in the final
code chunk of this vignette).

## Functions for generating outcomes

This function should take a single argument (`allocs`), a character vector
containing the allocations (names of trial `arms`) of all patients included
since the last adaptive analysis. The function must return a *numeric* vector,
regardless of actual outcome type (so, e.g., categorical outcomes must be
encoded as numeric). The returned numeric vector must be of the same length, and
the values in the same order as `allocs`. That is, the third element of `allocs`
specifies the allocation of the third patient randomised since the last adaptive
analysis, and (correspondingly) the third element of the returned vector is that
patient's outcome.

It sounds complicated, but it becomes clearer when we actually specify the 
function (essentially a re-implementation of the built-in function used by
`setup_trial_binom()`):

```{r}
get_ys_binom_custom <- function(allocs) {
  # Binary outcome coded as 0/1 - prepare returned vector of appropriate length
  y <- integer(length(allocs))
  
  # Specify trial arms and true event probabilities for each arm
  # These values should exactly match those supplied to setup_trial
  # NB! This is not validated, so this is the user's responsibility
  arms <- c("Control", "Experimental arm A", "Experimental arm B")
  true_ys <- c(0.25, 0.27, 0.20)
  
  # Loop through arms and generate outcomes
  for (i in seq_along(arms)) {
    # Indices of patients allocated to the current arm
    ii <- which(allocs == arms[i])
    # Generate outcomes for all patients allocated to current arm
    y[ii] <- rbinom(length(ii), 1, true_ys[i])
  }
  
  # Return outcome vector
  y
}
```

To illustrate how the function works, we first generate random allocations for
50 patients using equal allocation probabilities, the default behaviour of
`sample()`. By enclosing the call in parentheses, the resulting
allocations are printed:

```{r}
(allocs <- sample(c("Control", "Experimental arm A", "Experimental arm  B"),
                  size = 50, replace = TRUE))
```

Next, we generate random outcomes for these patients: 
```{r}
(ys <- get_ys_binom_custom(allocs))
```

## Functions for drawing posterior samples

The `setup_trial_binom()` function uses *beta-binomial conjugate prior models*
in each arm, with `beta(1, 1)` priors. These priors are uniform (≈
"non-informative") on the probability scale and corresponds to the same amount
of information as provided by 2 patients (1 with and 1 without an event),
described in greater detail by, e.g., *Ryan et al, 2019*
([10.1136/bmjopen-2018-024256](https://doi.org/10.1136/bmjopen-2018-024256)).

Our custom function for generating posterior draws also uses beta-binomial
conjugate prior models, but with informative priors. Informative
priors may prevent undue influence of random, early fluctuations in a trial by
pulling posterior estimates closer to the prior when limited data are available.

We seek relatively weakly informative priors centred on previous knowledge (or
beliefs), and so before we can actually define a function for generating
posterior draws based on informative priors, we need to derive such a prior.

### Informative priors

We assume that we have prior knowledge corresponding to a belief that the best
estimate for the true event probability in the `control` arm is 0.25 (25%),
and that the true event probability is between 0.15 and 0.35 (15-35%) with
95% probability. The mean of a beta distribution is simply
`[number of events]/[number of patients]`.

To derive a beta distribution that reflects this prior belief, we use
`find_beta_params()`, a helper function included in the `adaptr` (see
`?find_beta_params` for details):

```{r}
find_beta_params(
  theta = 0.25, # Event probability
  boundary = "lower",
  boundary_target = 0.15,
  interval_width = 0.95
)
```

We thus see that our prior belief prior roughly corresponds to previous
randomisation of 60 patients of whom 15 (`alpha`) experienced the event and 45
(`beta`) did not.

Even though we may expect event probabilities to differ in the non-control arms,
for this example we consider this prior appropriate for all arms as we consider
event probabilities smaller/larger than those represented by the prior unlikely.

Below, we illustrate the effects of this prior compared to the default
`beta(1, 1)` prior used in `setup_trial_binom()` for a single trial arm with 20
patients randomised, 12 events and 8 non-events. This corresponds to an
estimated event probability of 0.6 (60%), far more than the expected 0.25 (25%).
This could come about by random fluctuations when few patients have been
randomised, even if our prior beliefs are correct.

```{r, echo = FALSE}
# Helper function for plotting priors and posteriors together in the vignette
plot_prior_posterior <- function(events, non_events, prior_events, prior_non_events) {
  beta_curve <- function(prior_alpha, prior_beta, events = 0, non_events = 0, type = "prior") {
    ggplot2::geom_function(
      ggplot2::aes(colour = sprintf("beta(%i, %i) prior", prior_alpha, prior_beta), linetype = type), 
      fun = dbeta, args = list(shape1 = prior_alpha + events, shape2 = prior_beta + non_events)
    )
  }
  
  ggplot2::ggplot(data.frame(x = seq(from = 0, to = 1, by = 0.01)), ggplot2::aes(x = x)) +
    beta_curve(1, 1) +
    beta_curve(1, 1, events, non_events, "posterior") +
    beta_curve(prior_events, prior_non_events) +
    beta_curve(prior_events, prior_non_events, events, non_events, "posterior") +
    ggplot2::theme_minimal() +
    ggplot2::theme(legend.title = ggplot2::element_blank()) +
    ggplot2::scale_colour_brewer(palette = "Set1") +
    ggplot2::scale_x_continuous(breaks = 0.1 * 0:10, minor_breaks = 0.05 * 0:20, 
                                labels = paste0(0:10 * 10, "%"), expand = c(0.01, 0.01)) +
    ggplot2::scale_y_continuous(breaks = NULL, expand = c(0, 0), name = NULL) + 
    ggplot2::labs(x = "Event probability",
                  title = sprintf("Single trial arm - %i patients (%i events, %i non-events)", events + non_events, events, non_events)
  )
}
```

```{r prior_posterior_20, echo=FALSE}
plot_prior_posterior(12, 8, 15, 45)
```

Next, we illustrate the effects on the same prior after 200 patients have been
randomised to the same arm, with 56 events and 144 non-events, which
corresponds to an estimated event probability of 0.28 (28%), which
is more similar to the expected event probability.

```{r prior_posterior_200, echo=FALSE}
plot_prior_posterior(56, 144, 15, 45)
```

When comparing to the previous plot, we clearly see that once more patients have
been randomised, the larger sample of observed data starts to dominate the
posterior, so that the prior exerts less influence on the posterior
distribution (both posterior distributions are very alike despite the very
different prior distributions).

### Defining the function to generate posterior draws

There are a number of important things to be aware of when specifying this
function. First, it must accept all the following arguments (with the exact same
names, *even* if they are not used in the function):

- `arms`: character vector of the *currently active* `arms` in the trial.
- `allocs`: character vector of the allocations (trial `arms`) of *all*
patients randomised in the trial, *including* those randomised to `arms`
that are no longer active.
- `ys`: numeric vector of the outcomes of *all* patients randomised in the
trial, *including* those randomised to `arms` that are no longer active.
- `control`: single character, the *current* `control` arm; `NULL` in trials
without a common `control`.
- `n_draws`: single integer, the number of posterior draws to generate for
*each* arm.

Alternatively, unused arguments can be left out and the ellipsis (`...`)
included as the final argument in your function.

Second, the order of `allocs` and `ys` must match: the fifth element of `allocs`
represents the allocation of the fifth patient, while the fifth element of `ys`
represent the outcome of the same patient.

Third, `allocs` and `ys` are provided for **all** patients, including
those randomised to `arms` that are no longer active. This is done as users in
some situations may want to use these data when generating posterior draws for
the other (and currently **active**) `arms`.

Fourth, `adaptr` does not restrict how posterior samples are drawn.
Consequently, Markov chain Monte Carlo- or variational inference-based methods
may be used, and other packages supplying such functionality may be called by
user-provided functions. However, using more complex methods than simple
conjugate models substantially increases simulation run time. Consequently,
simpler models are well-suited for use in simulations.

Fifth, the function must return a `matrix` of numeric values with `length(arms)`
columns and `n_draws` rows, with the **currently active** `arms` as column
names. That is, each row must contain one posterior draw for each arm. `NA`'s
are not allowed, so even if no patients have been randomised to an arm yet,
valid numeric values should be returned (e.g., drawn from the prior or from
another very diffuse *posterior* distribution). Even if the outcome is not truly
numeric, both the vector of outcomes provided to the function (`ys`) and the
returned `matrix` with posterior draws must be encoded as numeric.

With this in mind, we are ready to specify the function:

```{r}
get_draws_binom_custom <- function(arms, allocs, ys, control, n_draws) {
  # Setup a list to store the posterior draws for each arm
  draws <- list()
  
  # Loop through the ACTIVE arms and generate posterior draws
  for (a in arms) {
    # Indices of patients allocated to the current arm
    ii <- which(allocs == a)
    # Sum the number of events in the current arm
    n_events <- sum(ys[ii])
    # Compute the number of patients in the current arm
    n_patients <- length(ii)
    # Generate draws using the number of events, the number of patients
    # and the prior specified above: beta(15, 45)
    # Saved using the current arms' name in the list, ensuring that the
    # resulting matrix has column names corresponding to the ACTIVE arms
    draws[[a]] <- rbeta(n_draws, 15 + n_events, 45 + n_patients - n_events)
  }
  
  # Bind all elements of the list column-wise to form a matrix with
  # 1 named column per ACTIVE arm and 1 row per posterior draw.
  # Multiply result with 100, as we're using percentages and not proportions
  # in this example (just to correspond to the illustrated custom function to
  # generate RAW outcome estimates below)
  do.call(cbind, draws) * 100
}
```

We now call the function using the previously generated `allocs` and `ys`.

To avoid cluttering, we only generate 10 posterior draws from each arm in this
example:

```{r}
get_draws_binom_custom(
  # Only currently ACTIVE arms, but all are considered active at this time
  arms = c("Control", "Experimental arm A", "Experimental arm B"),
  allocs = allocs, # Generated above
  ys = ys, # Generated above
  # Input control arm, argument is supplied even if not used in the function
  control = "Control",
  # Input number of draws (for brevity, only 10 draws here)
  n_draws = 10
)
```

**Importantly**, less than 100 posterior draws from each arm is not allowed when
setting up the trial specification, to avoid unstable results (see
`setup_trial_binom()`).

## Specifying the function to calculate raw outcome estimates

Finally, a custom function may be specified to calculate **raw** summary
estimates in each arm; these raw estimates are *not* posterior estimates, but
can be considered the maximum likelihood point estimates in this example. This
function must take a numeric vector (all outcomes in an arm) and return a single
numeric value. This function is called separately for each arm. Because we
express results as *percentages* and not *proportions* in this example, this
function simply calculates the outcome percentage in each arm:

```{r}
fun_raw_est_custom <- function(ys) {
  mean(ys) * 100
}
```

We now call this function on the outcomes in the `"Control"` arm, as an example:

```{r}
cat(sprintf(
  "Raw outcome percentage estimate in the 'Control' group: %.1f%%", 
  fun_raw_est_custom(ys[allocs == "Control"])
))
```

## Setup the trial specification

With all the functions defined, we can now setup the trial specification. As
stated above, some validation of all the custom functions is carried out when
the trial is setup:

```{r}
setup_trial(
  arms = c("Control", "Experimental arm A", "Experimental arm B"),
  
  # true_ys, true outcome percentages (since posterior draws and raw estimates
  # are returned as percentages, this must be a percentage as well, even if
  # probabilities are specified as proportions internally in the outcome
  # generating function specified above
  true_ys = c(25, 27, 20),
  
  # Supply the functions to generate outcomes and posterior draws
  fun_y_gen = get_ys_binom_custom,
  fun_draws = get_draws_binom_custom,
  
  # Define looks
  max_n = 2000,
  look_after_every = 100,
  
  # Define control and allocation strategy
  control = "Control",
  control_prob_fixed = "sqrt-based",
  
  # Define equivalence assessment - drop non-control arms at > 90% probability
  # of equivalence, defined as an absolute difference of 10 %-points
  # (specified on the percentage-point scale as the rest of probabilities in
  # the example)
  equivalence_prob = 0.9,
  equivalence_diff = 10,
  equivalence_only_first = TRUE,
  
  # Input the function used to calculate raw outcome estimates
  fun_raw_est = fun_raw_est_custom,
  
  # Description and additional information
  description = "custom trial [binary outcome, weak priors]",
  add_info = "Trial using beta-binomial conjugate prior models and beta(15, 45) priors in each arm."
)
```

As `setup_trial()` runs with no errors or warnings, the custom trial has been
successfully specified and may be run by `run_trial()` or `run_trials()` or
calibrated by `calibrate_trial()`. If the custom functions provided to
`setup_trial()` calls other custom functions (or uses objects defined by the
user outside these functions) or if functions loaded from non-base `R` packages
are used, please be aware that exporting these objects/functions or prefixing
them with their namespace is necessary if simulations are conducted using
multiple cores. See `setup_cluster()` and `run_trial()` for additional details
on how to export necessary functions and objects.
