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).
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):
find_beta_params(
theta = 0.25, # Event ratio
boundary = "lower",
boundary_target = 0.15,
interval_width = 0.95
)
#> alpha beta p2.5 p50.0 p97.5
#> 1 15 45 0.1498208 0.2472077 0.3659499We 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.
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.
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 you do not used them in your function):
arms: character vector of the currently activearmsin the trial.allocs: character vector of the allocations (trialarms) of all patients randomised in the trial, including those randomised toarmsthat are no longer active.ys: numeric vector of the outcomes of all patients randomised in the trial, including those randomised toarmsthat are no longer active.control: single character, the currentcontrolarm;NULLin trials without a commoncontrol.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 numeric matrix 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:
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:
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
)
#> Control Experimental arm A Experimental arm B
#> [1,] 30.96555 29.34973 29.26143
#> [2,] 30.47382 23.22668 25.08249
#> [3,] 31.04807 31.76577 19.81416
#> [4,] 17.00712 24.30809 16.36256
#> [5,] 21.31251 27.74615 22.63147
#> [6,] 25.50944 24.16283 30.29049
#> [7,] 16.60420 29.49526 28.75436
#> [8,] 25.17899 33.29374 30.87149
#> [9,] 23.72043 27.78537 29.89836
#> [10,] 30.50004 28.43694 26.62115Importantly, 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()).