# Long-Term Lithium Abundance Signatures following Planetary Engulfment

Jason Sevilla11affiliationmark:  , Aida Behmard22affiliationmark:  , Jim
Fuller11affiliationmark:  1Astronomy Department, California Institute of
Technology, Pasadena, CA 91125, USA  2Division of Geological and Planetary
Sciences, California Institute of Technology, Pasadena, CA 91125, USA

###### Abstract

Planetary engulfment events can occur while host stars are on the main
sequence. The addition of rocky planetary material during engulfment will lead
to refractory abundance enhancements in the host star photosphere, but the
level of enrichment and its duration will depend on mixing processes that
occur within the stellar interior, such as convection, diffusion, and
thermohaline mixing. We examine engulfment signatures by modeling the
evolution of photospheric lithium abundances. Because lithium can be burned
before or after the engulfment event, it produces unique signatures that vary
with time and host star type. Using MESA stellar models, we quantify the
strength and duration of these signatures following the engulfment of a 1, 10,
or 100 $M_{\oplus}$ planetary companion with bulk Earth composition, for
solar-metallicity host stars with masses ranging from 0.5$-$1.4 $M_{\odot}$.
We find that lithium is quickly depleted via burning in low-mass host stars
($\lesssim 0.7\,M_{\odot}$) on a time scale of a few hundred Myrs, but
significant lithium enrichment signatures can last for Gyrs in G-type stars
($\sim\\!0.9\,M_{\odot}$). For more massive stars (1.3$-$1.4 $M_{\odot}$),
engulfment can enhance internal mixing and diffusion processes, potentially
decreasing the surface lithium abundance. Our predicted signatures from
exoplanet engulfment are consistent with observed lithium-rich solar-type
stars and abundance enhancements in chemically inhomogeneous binary stars.

###### Subject headings:

diffusion — instabilities — planet-star interactions — stars: abundances —
stars:interiors

##  1\. Introduction

Refractory element abundances in stellar photospheres can be used as tracers
for astrophysical events. Lithium is a particularly useful tracer because it
is consumed in thermonuclear reactions at relatively low temperatures
($T\approx 3\times 10^{6}\,{\rm K})$ compared to other refractory species, and
is thus depleted over stellar lifetimes. This implies that lithium abundances
can shed light on recent events that altered stellar chemistry beyond birth
compositions.

Photospheric lithium is depleted over time by convection and other stellar
interior processes that can mix it down into the lithium burning region.
Observations (see summary in somers2016) show that lithium abundances for such
stars evolve and become depleted on timescales of millions to billions of
years, and are thus affected by stellar interior processes acting throughout
main sequence (MS) lifetimes. dumont2021 found that replicating such MS
lithium depletion requires consideration of non-standard mixing processes.
Thus, predictions of how long certain lithium abundance patterns will persist
are strongly dependent on which mixing processes are involved, e.g.,
convection, thermohaline mixing, gravitational settling, element diffusion,
mixing via shear instabilities, etc.

Observed lithium abundance patterns may also be affected by planet formation
and evolution processes, such as planetary engulfment. However, the evolution
of lithium enrichment signatures resulting from engulfment events is not
totally understood. sandquist2002 conjectured that lithium enrichments may
indicate stellar pollution resulting from planetary engulfment.
soaresfurtado2021 found significant lithium enrichment after modeling
planetary engulfment events, but they did not consider effects from diffusion,
overshoot mixing, and thermohaline mixing. On the other hand, theado2012
included thermohaline mixing, and found that engulfment could actually deplete
surface lithium abundances below pre-engulfment levels. As we shall see, these
processes can alter surface abundances by orders of magnitude and deserve
careful consideration.

If planetary engulfment does affect photospheric lithium abundances, we might
expect to see a difference in lithium observations of planet host stars
compared to stars lacking known planets. israelian2009 observed lithium
depletion in solar-type stars that harbor exoplanets compared to field stars.
These lithium patterns may have resulted from prior engulfment of inner
planets, or from sequestration of lithium within observed planets. Conversely,
baumann2010 performed statistical tests to examine lithium abundances between
two stellar populations, namely metal-rich solar analogues with and without
observed exoplanets. They found no significant difference between the two
samples, casting doubt on possible correlations between lithium abundances and
planetary engulfment.

One particular type of mixing, thermohaline instability, is expected to be
especially strong after engulfment events have occurred (vauclair2004;
garaud2011; bauer2018; bauer2019). This is because thermohaline mixing is
driven by an inverse mean molecular weight gradient, which is present after
engulfment when heavy planetary material is deposited within outer layers of
the engulfing star. Thermohaline mixing could attenuate lithium enrichment
over time, or even cause lithium depletion below the primordial level. This
mixing process was the cause of the lithium depletion signatures observed by
theado2012. While these results are valuable, they cannot be extrapolated to a
wide range of stellar types because theado2012 focused primarily on solar-type
stars.

To elucidate the connection between lithium enrichment and planetary
engulfment while considering relevant mixing processes, we ran stellar models
with the MESA stellar evolution code (paxton2011; paxton2013; paxton2015;
paxton2018; paxton2019). In Section 2, we discuss our MESA stellar models and
implementation of non-standard mixing processes such as thermohaline
instabilities. In Section LABEL:sec:results, we present analysis of the
lithium abundances following planetary engulfment when different mixing
processes are present, with varying amounts of accreted mass, and in
comparison to models without engulfment. We summarize observable engulfment
signature timescales in Section LABEL:sec:timescale, and discuss their
implications for the evolution of planet engulfment signatures in Section
LABEL:sec:discussion.

##  2\. Stellar Models

We computed our stellar models using the open-source 1D stellar evolution code
MESA (paxton2011; paxton2013; paxton2015; paxton2018; paxton2019). This allows
us to simulate mixing processes in the stars after planetary engulfment and to
monitor surface lithium abundance over time. We ran non-rotating stellar
models with zero-age main sequence (ZAMS) masses of 0.5–1.4 $M_{\odot}$, with
solar metallicities of $Z$ = 0.017. Models that included accretion were run in
three stages. In the first stage, we evolved the star up to its ZAMS phase.
Planetary engulfment was simulated via accretion of bulk Earth composition
material in the second stage (see Section LABEL:sec:accrete). In the third
stage, we evolved the star up to the end of its MS lifetime. The surface
lithium abundance in this final stage is the data of primary interest. We also
ran models without accretion for comparison. In the first stage, convection is
the only mixing process active. In the second and third stages, each model
utilized a selection of mixing settings that included convective overshoot
mixing, thermohaline mixing, atomic diffusion, and a minimum D_mix coefficient
to model additional relevant (but poorly understood) mixing processes.

###  2.1. Input Physics

Convective overshoot represents mixing that occurs near the convective
boundaries of the star. We used an exponential scheme, which MESA takes from
herwig2000. The formula for the overshoot mixing coefficient is given in
paxton2011:

$D_{\rm ov}=D_{\rm conv,0}\hskip
1.99168pt\exp\left(-\frac{2z}{f\lambda_{P,0}}\right),\hskip 2.84526pt$ (1)

where $D_{\rm conv,0}$ the diffusion coefficent taken from a point determined
by the user, $\lambda_{P,0}$ is the pressure scale height at this point, $z$
is the distance away from this point, and $f$ determines the characteristic
size of the overshooting region in terms of $\lambda_{P,0}$. The settings in
the MESA inlist allow us to input $f$ as well as $f_{0}$, which is how many
scale heights into the convective zone the aforementioned point will be placed
to calculate $D_{\rm conv,0}$. We chose $f=0.02$ and $f_{0}=0.005$ in all of
our runs.

Thermohaline mixing is a double diffusive instability that occurs in the
presence of an inverse mean molecular weight gradient and a stabilizing
entropy gradient, i.e., material of higher mean molecular weight and entropy
lying above material of lower mean molecular weight and entropy. The material
with high mean molecular weight sinks downwards in long strands, and material
with lower mean molecular weight rises up. This type of mixing is especially
relevant after planetary engulfment, where heavy planetary material is
deposited near the stellar surface. The diffusive mixing from thermohaline has
been calibrated from numerical simulations (e.g., denissenkov2010;
traxler2011; harrington2019).

We utilize MESA’s thermohaline mixing prescription based on brown2013 in our
MESA models. brown2013 provides a more accurate prescription of thermohaline
mixing compared to previous implementation (e.g., kippenhahn1980) via an
improved fingering convection model that is supported by 3D numerical
simulations (zemskova2014). The thermohaline instability is only active when
the composition gradient $\nabla_{\mu}=d\ln\mu/d\ln P$ is sufficiently
negative, and the effective thermohaline diffusivity $D_{\rm th}$ approaches
zero when $\nabla_{\mu}$ approaches the stability boundary (see discussion in
bauer2019).

The prescription for how MESA handles elemental diffusion is provided in
paxton2018, using Burgers’ diffusion equations. All isotopes are assigned to a
diffusion class that covers a range of atomic masses. Each isotope in a class
is treated identically to that class’s representative isotope. These diffusion
classes are summarized in Table 2.1.

Table 1List of representative isotopes and atomic mass ranges for each diffusion class.  Representative Isotope | Atomic Mass Range
---|---
1H | $1-2$
4He | $3-4$
7Li | $5-7$
12C | $8-12$
16O | $13-16$
24Mg | $17-24$
28Si | $25-28$
56Fe | $>28$


# Bayesian design and analysis of external pilot trials for complex
interventions

Duncan T. Wilson1 James M. S. Wason2,3 Julia Brown1 Amanda J. Farrin1 Rebecca
E. A. Walwyn1

(1 - Leeds Institute of Clinical Trials Research, University of Leeds, Leeds,
UK
2 - Institute of Health and Society, Newcastle University, Newcastle, UK
3 - MRC Biostatistics Unit, University of Cambridge, Cambridge, UK)

###### Abstract

External pilot trials of complex interventions are used to help determine if
and how a confirmatory trial should be undertaken, providing estimates of
parameters such as recruitment, retention and adherence rates. The decision to
progress to the confirmatory trial is typically made by comparing these
estimates to pre-specified thresholds known as progression criteria, although
the statistical properties of such decision rules are rarely assessed. Such
assessment is complicated by several methodological challenges, including the
simultaneous evaluation of multiple endpoints, complex multi-level models,
small sample sizes, and uncertainty in nuisance parameters. In response to
these challenges, we describe a Bayesian approach to the design and analysis
of external pilot trials. We show how progression decisions can be made by
minimising the expected value of a loss function, defined over the whole
parameter space to allow for preferences and trade-offs between multiple
parameters to be articulated and used in the decision making process. The
assessment of preferences is kept feasible by using a piecewise constant
parametrisation of the loss function, the parameters of which are chosen at
the design stage to lead to desirable operating characteristics. We describe a
flexible, yet computationally intensive, nested Monte Carlo algorithm for
estimating operating characteristics. The method is used to revisit the design
of an external pilot trial of a complex intervention designed to increase the
physical activity of care home residents.

##  1 Introduction

Complex interventions, defined as those comprised of several interacting
components [1], can be challenging to evaluate in randomised controlled trials
(RCTs) due to factors such as slow patient recruitment, poor levels of
adherence to the intervention, and low completeness of follow-up data. To
identify these problems prior to the main RCT we often conduct small trials
[1] known as pilots. These typically take the same form as the planned RCT but
with a considerably lower sample size [2]. If there is a seamless transition
between the pilot and the main RCT, with all data being pooled and used in the
final analysis, they are known as internal pilots. External pilots, in
contrast, are carried out separately to the main RCT with a clear gap between
the two trials. Pilot trials, which aim to inform the feasibility and optimal
design of a subsequent definitive trial [3], are distinct from phase II
trials, which focus instead on assessing potential efficacy and safety.

The data generated by an external pilot trial is used to help decide if the
main RCT should go ahead, and if so, whether the intervention or the trial
design should be adjusted to ensure success. In the UK, the National Institute
for Health Research ask that these _progression criteria_ are pre-specified
and included in the research plan [4], and the recent CONSORT extension to
randomised pilot trials requires their reporting [5]. A single pilot trial can
collect data on several progression criteria, often focused on the
aforementioned areas of recruitment, protocol adherence, and data collection
[6]. Although they may take the form of single threshold values leading to
binary stop/go decision rules, investigators are increasingly using two
thresholds to accommodate an intermediate decision between stopping altogether
and progressing straight to the main trial, which would allow progression but
only after some adjustments have been made [5]. The need for appropriate
progression criteria is clear when we consider the consequences of poor post-
pilot progression decisions. If the criteria are too lax, there is a greater
risk that the main trial will go ahead but found to be infeasible and thus a
waste of resources; if the criteria are too strict, a promising intervention
may be discarded under the mistaken belief that the main trial would be
infeasible. Despite this, there is little published guidance about how they
should be determined[6, 7].

In addition to pre-specifying progression criteria, another key design
decision is the choice of pilot sample size. Conventional methods of sample
size determination, which focus on ensuring the trial will have sufficient
power to detect a target difference in the primary outcome, are rarely used
since they would lead to a pilot sample size comparable with the main trial
sample size. Several methods for pilot sample size determination instead aim
to provide a sufficiently precise estimate of the variance in the primary
outcome measure to inform the sample size of the main trial [8, 9, 10, 11, 12,
13]. Others have suggested a simple rule of thumb for when the goal is to
identify unforeseen problems [14]. While some have noted that the low sample
size in pilots may lead to a considerable probability that a certain
progression criterion will be met (or missed) due to random sampling variation
[12, 15], and despite the consequences of making the wrong progression
decision, the statistical properties of pilot decision rules are rarely used
to inform the choice of sample size. This may be due to the methodological
challenges commonly found in pilot trials of complex interventions, including
the simultaneous evaluation of multiple endpoints, complex multi-level models,
small sample sizes, and prior uncertainty in nuisance parameters [16].

In this paper we will describe a method for designing and analysing external
pilot trials which addresses these challenges. We take a Bayesian view, where
progression decisions are made to minimise the expected value of a loss
function. We propose a loss function with three parameters whose values can be
determined either through direct elicitation of preferences or by considering
the pilot trial operating characteristics they lead to. The operating
characteristics we propose are all unconditional probabilities (with respect
to a prior distribution) of making incorrect decisions, also known as
assurances [17]. Using assurances rather than the analogous frequentist error
rates brings several benefits, including the ability to make use of existing
knowledge whilst allowing for any uncertainty, and a more natural
interpretation [18]. As we will show, assurances are also useful when our
preferences for different end-of-trial decisions are based on several
attributes in a complex way that involves trading off some against others.

The remainder of this paper is organised as follows. In Section 2 we describe
the general framework for pilot design and analysis, some operating
characteristics used for evaluation, and a routine for optimising the design.
Two illustrative examples are then described in Sections 3 and 4. Finally, we
discuss implications and limitations in Section 5.

##  2 Methods

###  2.1 Prior specification

Consider a pilot trial which will produce data $x$ according to model
$p(x|\theta)$. We decompose the parameters into $\theta=(\phi,\psi)$, where
$\phi$ denotes the parameters of substantive interest and $\psi$ the nuisance
parameters. We follow [19] and assume that two joint prior distributions of
$\theta$ have been specified. First, the _analysis_ prior $p_{A}(\theta)$ is
that which will be used when fitting the model once the pilot data is
obtained. It has been argued that regulators are unlikely to accept the prior
beliefs of the trial sponsor for analysis of the data [17, 20], and as such a
weakly or non-informative prior should be used for $p_{A}(\theta)$ in order to
“let the data drive the inference” [19]. The choice of such a prior will be
depend on the specific model being used, although methodological guidance for
various specific cases such as logistic regression [21] and hierarchical
models [22] is available. It should be emphasised, however, that the typically
small sample size of a pilot trial can mean the effect of the analysis prior
is non-negligible. As such, the analysis prior should provide a credible and
justifiable representation of prior ignorance, avoiding extreme default
choices which may place too much prior weight on infeasible regions of the
parameter space.

The _design_ prior $p_{D}(\theta)$ will be used when evaluating the
statistical performance of a proposed pilot trial design. It may be considered
as purely hypothetical in the spirit of a ‘what-if’ analysis [19], in which
case several candidate design priors may be suggested and performance
evaluated under each of these. Alternatively, and as we will assume in the
remainder of this paper, $p_{D}(\theta)$ can be a completely subjective prior
which fully expresses our knowledge and uncertainty in the parameters at the
design stage. Although eliciting such a prior is potentially challenging, many
examples describing successful practical applications of expert elicitation
for clinical trial design are available [20, 18, 23], as are tools for its
conduct such as the Sheffield Elicitation Framework (SHELF) [24].

###  2.2 Analysis and progression decisions

After observing the pilot data $x$, we must decide whether or not to progress
to the main RCT. We consider three possible actions following the
aforementioned ‘traffic light’ system commonly used in pilot trials:

  * •

red - discard the intervention and stop all future development or evaluation;

  * •

amber - proceed to the main RCT, but only after some modifications to the
intervention, the planned trial design, or both; or

  * •

green - proceed immediately to the main RCT.

In what follows we will denote these decisions by $r,a$ and $g$ respectively.
We assume that our preferences between the three possible decisions are
influenced by $\phi$ but independent of $\psi$, formalising the separation of
$\theta$ into substantive and nuisance components. We partition the
substantive parameter space $\Phi$ into three subspaces $\Phi_{I}$, for
$I=R,A,G$. Each subspace label corresponds to the decision we would make if we
knew the true value of $\phi$. For example, if $\phi\in\Phi_{R}$ then the
optimal decision is $r$(ed) - halt development and do not proceed to a
definitive trial. We will henceforth refer to these three subsets as
_hypotheses_. Throughout, we will distinguish hypothesis $I$ from the
corresponding optimal decision $i$ by using upper and lower case letters
respectively.

When $\phi\in\Phi_{I}$ and we choose a decision $j\neq i$, there will be
negative consequences. In particular, we may make three kinds of mistakes:
proceed to an infeasible main RCT; discard a promising intervention; or make
unnecessary adjustments to the intervention or trial design. We denote these
errors as $E_{1}$, $E_{2}$, $E_{3}$ respectively. The occurrence of error $j$
will be denoted by $E_{j}=1$, otherwise $E_{j}=0$. An error’s occurrence will
be a function of the decision made $d$ and the true parameter value $\phi$,
i.e. $E_{j}(d,\phi):\\{r,a,g\\}\times\Phi\rightarrow\\{0,1\\}$ for $j=1,2,3$.
We then use a loss function to express the preferences of the decision
maker(s) on the space of possible events $E_{1}\times E_{2}\times E_{3}$ under
uncertainty, defined as

$L(d,\phi)=c_{1}E_{1}(d,\phi)+c_{2}E_{2}(d,\phi)+c_{3}E_{3}(d,\phi).$

Note that the additive form of the loss function implies that the our
preferences for any one of the attributes $E_{1},E_{2},E_{3}$ are independent
of the values taken by the others [25].

To determine appropriate values of the parameters $c_{1},c_{2},c_{3}$ we first
scale the loss function by setting $c_{1}+c_{2}+c_{3}=1$. Thus, a loss of 0 is
obtained if no errors occur, and a loss of 1 is obtained if all errors occur
(although note that this is not possible in this setting). We then follow the
procedure described by French and Rios Insua (page 99) [25], eliciting some
judgements from the decision maker(s) and using these to determine the values
of $c_{1},c_{2},c_{3}$. One such judgement involves a simple gamble of
obtaining the event $(E_{1}=0,E_{2}=0,E_{3}=0)$ with probability $1-p_{1}$ and
the event $(E_{1}=1,E_{2}=0,E_{3}=1)$ with probability $p_{1}$. The decision
maker is asked to compare this gamble against an alternative of obtaining the
event $(E_{1}=1,E_{2}=0,E_{3}=0)$ for certain, and to adjust the value of
$p_{1}$ until they feel indifferent between the two options. Since this
indifference implies the expected losses of the two options are equal, we will
then have

$p_{1}(c_{1}+c_{3})=c_{1}.$

Similarly, we can ask the decision maker(s) to consider a gamble between the
event $(E_{1}=0,E_{2}=0,E_{3}=0)$ with probability $1-p_{2}$ and the event
$(E_{1}=1,E_{2}=1,E_{3}=0)$ with probability $p_{2}$, and compare this against
the option of obtaining $(E_{1}=1,E_{2}=0,E_{3}=0)$ for certain. Again, by
determining the value of $p_{2}$ which corresponds to indifference and thus
equal expected loss, we deduce that

$p_{2}(c_{1}+c_{2})=c_{1}.$

This gives three equations which can be solved to obtain

$c_{1}=\frac{-p_{1}p_{2}}{p_{1}p_{2}-p_{1}-p_{2}},~{}c_{2}=\frac{p_{1}p_{2}-p_{1}}{p_{1}p_{2}-p_{1}-p_{2}}~{}c_{3}=\frac{p_{1}p_{2}-p_{2}}{p_{1}p_{2}-p_{1}-p_{2}}.$

Note that the two specific judgements suggested here are only two of many
possible similar questions which could be posed to the decision maker(s). It
is recommended that more indifferences are elicited in order to seek out any
inconsistencies and further clarify their true preferences.

The loss function will then take values as given in Table 1. For example,
suppose we make a ‘green’ decision under the ‘amber’ hypothesis. The
subsequent trial will be infeasible because the necessary adjustments will not
have been made; but we have also discarded a promising intervention, since it
would have been redeemed had the adjustments been made. The overall loss is
therefore $c_{1}+c_{2}$.

Table 1: Losses associated with each decision under each hypothesis. |  | Hypothesis
---|---|---
|  | $\phi\in\Phi_{R}$ | $\phi\in\Phi_{A}$ | $\phi\in\Phi_{G}$
Decision | $r$ | 0 | $c_{2}$ | $c_{2}$
$a$ | $c_{1}+c_{3}$ | 0 | $c_{3}$
$g$ | $c_{1}$ | $c_{1}+c_{2}$ | 0

Given a loss function with parameters $\mathbf{c}=(c_{1},c_{2},c_{3})$ we
follow the principle of maximising expected utility (or in our case,
minimising the expected loss) when making a progression decision. We first use
the pilot data in conjugation with the analysis prior $p_{A}(\theta)$ to
obtain a posterior $p(\phi~{}|~{}x)$, and then choose the decision $i^{*}$
such that

$\displaystyle i^{*}$
$\displaystyle=\operatorname*{arg\,min}_{i\in\\{r,a,g\\}}\mathbb{E}_{\phi|x}[L(i,\phi)]$
(1) $\displaystyle=\operatorname*{arg\,min}_{i\in\\{r,a,g\\}}\int
L(i,\phi)p(\phi|x)d\phi.$ (2)

We can simplify this expression by noting that, given the piecewise constant
nature of the loss function, the expected loss of each decision depends only
on the posterior probabilities $p_{I}=Pr[\phi\in\Phi_{I}~{}|~{}x]$ for
$I=R,A,G$. We then have

$\displaystyle\mathbb{E}_{\phi|x}[L(r,\phi)]$
$\displaystyle=p_{A}c_{3}+p_{G}c_{3},$ (3)
$\displaystyle\mathbb{E}_{\phi|x}[L(a,\phi)]$
$\displaystyle=p_{R}c_{1}+p_{R}c_{2}+p_{G}c_{2},$ (4)
$\displaystyle\mathbb{E}_{\phi|x}[L(g,\phi)]$
$\displaystyle=p_{R}c_{1}+p_{A}c_{1}+p_{A}c_{3}.$ (5)

For some simple models which admit a conjugate analysis, the posterior
probabilities $p_{I}$ can be obtained exactly. Otherwise, Monte Carlo
estimates can be computed based on the samples from the joint posterior
distribution generated by an MCMC analysis of the pilot data. Specifically,
given $M$ samples $\phi^{(1)},\phi^{(2)},\ldots,\phi^{(M)}\sim
p(\phi~{}|~{}x)$,

$p_{I}\approx\frac{1}{M}\sum_{k=1}^{M}\mathbb{I}(\phi^{(k)}\in\Phi_{I}),$ (6)

where $\mathbb{I}(.)$ is the indicator function.

###  2.3 Operating characteristics

Defining a loss function and following the steps of the preceding section
effectively prescribes a decision rule mapping the pilot data sample space
$\mathcal{X}$ to the decision space $\\{r,a,g\\}$. To gain some insight at the
design stage into the properties of this rule, we propose to calculate some
trial operating characteristics. These take the form of unconditional
probabilities of making an error when following the rule, calculated with
respect to the design prior $p_{D}(\theta)$. We consider the following:

  * •

$OC_{1}=Pr[a~{}\&~{}\phi\in\Phi_{R}]+Pr[g~{}\&~{}\phi\in\Phi_{R}\cup\Phi_{A}]$
\- probability of proceeding to an infeasible main RCT;

  * •

$OC_{2}=Pr[r~{}\&~{}\phi\in\Phi_{A}\cup\Phi_{G}]+Pr[g~{}\&~{}\phi\in\Phi_{A}]$
\- probability of discarding a promising intervention;

  * •

$OC_{3}=Pr[a~{}\&~{}\phi\in\Phi_{R}\cup\Phi_{G}]$ \- probability of making
unnecessary adjustments to the intervention or the trial design.

These operating characteristics can be estimated using simulation. First, we
draw $N$ samples
$(\theta^{(1)},x^{(1)}),(\theta^{(2)},x^{(2)}),\ldots,(\theta^{(N)},x^{(N)})$
from the joint distribution $p(\theta,x)=p(x|\theta)p_{D}(\theta)$. For each
data set we then apply the analysis and decision making procedure described in
Section 2.2, using some vector $\mathbf{c}$ to parametrise the loss function.
This results in $N$ decisions $i^{(k)}$ which can be contrasted with the
corresponding true parameter value $\theta^{(k)}$ and in which hypothesis it
resides, noting if any of the three types of errors have been made. MC
estimates of the operating characteristics can then be calculated as the
proportion of occurrences of each type of error in the $N$ simulated cases.
Assuming that $N$ is large, the unbiased MC estimate of an operating
characteristic with true probability $p$ will be approximately normally
distributed with variance $p(1-p)/N$.111Note that in the case of complex
models which do no admit a conjugate analysis, the posterior probabilities
obtained using an MCMC analysis will themselves be approximate and as such the
optimal decision will be subject to error, which may increase the variance of
the operating characteristic estimates. However, this issue can be sidestepped
by assuming that, for each data set, the analysis that is simulated
corresponds exactly to the analysis that would be carried out in practice. In
particular, we assume that exactly $M$ posterior samples will be generated by
the same MCMC algorithm, using the same seed in the random number generator.

###  2.4 Optimisation

Elicitation of the loss function parameters $\mathbf{c}=(c_{1},c_{2},c_{3})$
in the manner described in Section 2.2 may be challenging, particularly when
multiple decision-makers are involved [26]. An alternative way to determine
$\mathbf{c}$ is through examining the operating characteristics it leads to
(for some fixed pilot design). As $\mathbf{c}$ is adjusted, the balance
between the conflicting objectives of minimising each OC will change, and the
task is then to find the $\mathbf{c}$ which returns the best balance from the
perspective of the decision-maker. Formally, and thinking of operating
characteristics as functions of $\mathbf{c}$, we wish to solve the multi-
objective optimisation problem

$\min_{\mathbf{c}\in\mathcal{C}}~{}\left(OC_{1}(\mathbf{c}),~{}OC_{2}(\mathbf{c}),~{}OC_{3}(\mathbf{c})\right)$
(7)

where $\mathcal{C}=\\{c_{1},c_{2}\in[0,1]~{}|~{}c_{1}+c_{2}\leq 1\\}$.

Since the three objectives are in conflict there will be no single solution
which simultaneously minimises each one. We would instead like to find a set
$\mathcal{C}^{*}=\\{\mathbf{c}^{(1)},\mathbf{c}^{(2)},\ldots,\mathbf{c}^{(K)}\\}$
such that each member provides a different balance between minimising the
three operating characteristics. If there exist
$\mathbf{c},\mathbf{c}^{\prime}\in\mathcal{C}^{*}$ such that
$OC_{i}(\mathbf{c}^{\prime})\leq OC_{i}(\mathbf{c})$ for all $i\in\\{1,2,3\\}$
and $OC_{i}(\mathbf{c}^{\prime})<OC_{i}(\mathbf{c})$ for some
$i\in\\{1,2,3\\}$, we say that $\mathbf{c}^{\prime}$ dominates $\mathbf{c}$.
In this case, because $\mathbf{c}$ leads to worse (or at least no better)
values of all three operating characteristics when compared to
$\mathbf{c}^{\prime}$, we have no reason to include it in our set
$\mathcal{C}^{*}$. Because the search space $\mathcal{C}$ has only two
dimensions, problem (7) can be approximately solved by generating a uniform
random sample of $\mathbf{c}$’s and estimating the operating characteristics
for each. Any parameters which are dominated in this set can then be
discarded, and the operating characteristics of those which remain can be
illustrated graphically. The decision maker(s) can then view the range of
available options, all providing different trade-offs amongst the three
operating characteristics, and choose from amongst them.

To solve the problem in a timely manner we must be able to estimate operating
characteristics quickly. Noting from equation (3) that the expected loss of
each decision depends only on $\mathbf{c}$ and the posterior probabilities
$p_{R},p_{A}$ and $p_{G}$, we first generate $N$ samples of these posterior
probabilities and then use this same set of samples for every evaluation. This
approach not only ensures that optimisation is computationally feasible, but
also means that differences in operating characteristics are entirely due to
differences in costs, as opposed to differences in the random posterior
probability samples.

##  3 Illustrative example - Child psychotherapy (TIGA-CUB)

TIGA-CUB (Trial on Improving Inter-Generational Attachment for Children
Undergoing Behaviour problems) was a two-arm, individually-randomised,
controlled pilot trial informing the feasibility and design of a confirmatory
RCT comparing Child Psychotherapy (CP) to Treatment as Usual (TaU), for
children with treatment resistant conduct disorders. The trial aimed to
recruit $60$ primary carer-child dyads, to be randomised equally to each arm.
This sample size was chosen to give desired levels of precision in the
estimates of the common standard deviation of the primary outcome, the follow-
up rate, and the adherence rate. Here, we focus on the latter two parameters
and consider how our proposed method could have informed the design of TIGA-
CUB.

We model the number of participants successfully followed-up (denoted $f$)
using a binomial distribution with parameter $p_{f}$, and similarly the number
successfully adhering to the intervention (denoted $a$) with a binomial
distribution with parameter $p_{a}$. For a fixed pilot trial per-arm sample
size $n$, the parameters of the model are $\phi=(p_{f},p_{a})$, with no
nuisance parameters. Further assuming that the numbers followed-up and
adhering are independent, the likelihood is then

$p(f,a|p_{f},p_{a})=\left[{2n\choose
f}p_{f}^{f}(1-p_{f})^{2n-f}\right]\times\left[{n\choose
a}p_{a}^{a}(1-p_{a})^{n-a}\right].$

At the design stage, the follow-up rate $p_{f}$ was thought to be somewhere in
the range 62% to 92%, while the adherence rate $p_{a}$ was thought to lie
between 40% and 95%. We reflect these ranges of uncertainty in our design
priors by using beta distributions $p_{f}\sim Beta(40,10)$ (thus giving a
prior mean of 0.8), and $p_{a}\sim Beta(11.2,4.8)$ (giving a prior mean of
0.7). We assume that a uniform ‘non-informative’ prior $Beta(1,1)$ will be
used for each parameter in the analysis.

TIGA-CUB’s progression criteria included only simple stop/go thresholds, with
no intermediate ‘amber’ decisions. As such, in this example we partition the
parameter space into two hypotheses, $\Phi_{G}$ and $\Phi_{R}$. For the
purposes of illustration we define the hypothesis $\Phi_{G}$ as the subset of
the parameter space where $p_{f}>=0.8$ and $p_{a}>=0.7$, hypothesis $\Phi_{R}$
being its complement. Thus, in this example we do not consider there to be a
trade-off between the two parameters of interest. For the main trial to be
feasible, both must be above their respective thresholds. The prior
distributions on parameters $p_{f}$ and $p_{a}$ imply an _a priori_
probability of 0.28 that $\phi\in\Phi_{G}$, i.e. that both follow-up and
adherence are sufficiently high.

In this special case, the loss function is

$L(d,\phi)=c_{1}E_{1}(d,\phi)+c_{2}E_{2}(d,\phi)$

and the expected losses of decisions $g$ and $r$ will be
$\mathbb{E}_{\phi|x}[L(g,\phi)]=c_{1}p_{R}$ and
$\mathbb{E}_{\phi|x}[L(r,\phi)]=c_{2}p_{G}$, where $p_{R}+p_{G}=1$ and
$c_{1}+c_{2}=1$. Decision $g$ is therefore optimal whenever $p_{G}>c_{1}$. The
posterior probability $p_{G}$ can be easily calculated given the pilot data
due to the beta prior distributions being conjugate. Specifically, given a
total sample size $n$ and observing $x_{f}$ participants with follow-up and
$x_{a}$ participants with adherence, the posterior probability
$Pr[\phi\in\Phi_{G}~{}|~{}x]$ is given by

$p_{G}=[1-F(0.8;1+x_{f},1+n-x_{f})]\times[1-F(0.7;1+x_{a},1+n/2-x_{a})],$ (8)

where $F(y;\alpha,\beta)$ denotes the cumulative probability function of the
beta distribution with parameters $\alpha,\beta$.

At the design stage we can calculate the probability of an infeasible trial
($OC_{1}$),

$\displaystyle Pr[g,\phi\in\Phi_{R}]$
$\displaystyle=\int_{\Phi_{R}}Pr[g~{}|~{}\phi]p(\phi)d\phi$ (9)
$\displaystyle=\int_{\Phi_{R}}\left(\sum_{x_{f}=0}^{n}\left[\sum_{x_{a}=0}^{n/2}\mathbb{I}(p_{G}<c_{1}~{}|~{}x_{f},x_{a},n)p(x_{a}~{}|~{}\phi)\right]p(x_{f}~{}|~{}\phi)\right)p(\phi)d\phi,$
(10)

and similarly for the probability of discarding a promising intervention. As
these calculations can be computationally expensive for moderate $n$ due to
the nested summation term, we use Monte Carlo approximations as described in
Section 2.

Keeping the sample size fixed at $n=30$ per arm, we estimated the operating
characteristics using a range of cost parameters values
$c_{1}=0,0.02,0.04,\ldots,1$ using $M=10^{6}$ Monte Carlo samples. The results
are plotted in Figure 1, with some specific values of $c_{1}$ highlighted. The
decision-maker can decide which point on the operating characteristic curve
best reflects their own priorities in terms of the two types of error. For
example, if the consequences of running an infeasible main RCT are considered
less important than those of needlessly discarding a potentially effective
intervention, the decision-maker may choose to set $c_{1}=0.2$ and would
obtain $OC_{1}=0.19,OC_{2}=0.05$.

Figure 1: Probabilities of an infeasible main trial ($OC_{1}$) and of
discarding a promising intervention ($OC_{2}$) for a range of loss parameters
$c_{1}$ when sample size is fixed at $n=30$.

To examine the effect of adjusting the sample size, we evaluated the operating
characteristics obtained for $n=10,12,14,\ldots,50$ per arm whilst setting
$c_{1}=0.2,0.36,0.5$. The results are shown in Figure 2. Each line includes a
shaded area denoting the 95% Monte Carlo error intervals, although these are
so small as to be illegible given the high number ($M=10^{6}$) of MC samples
used for each calculation. Although operating characteristics generally
improve as the sample size is increased, we see that for $c_{1}=0.36$ and 0.5
the probability of an infeasible main trial, $OC_{1}$, remains flat whilst
$OC_{2}$ has a downward trend. As we would expect, the the expected loss
reduces smoothly as $n$ increases in all cases. In contrast, there is some
variability beyond that explained by MC error in the OCs. This can be
explained by the discrete nature of simulated adherence and follow-up data.
Our results show that, for the design priors and hypotheses used in this
example, the chosen sample size in TIGA-CUB of $n=30$ can provide error rates
broadly in line with conventional type I and II error rates under the usual
hypothesis testing framework.

Figure 2: Probabilities of an infeasible main trial ($OC_{1}$) and of
discarding a promising intervention ($OC_{2}$) for a range of per-arm sample
sizes and different values of the loss parameter $c_{1}$.

##  4 Illustrative example - Physical activity in care homes (REACH)

The REACH (Research Exploring Physical Activity in Care Homes) trial aimed to
inform the feasibility and design of a future definitive RCT assessing a
complex intervention designed to increase the physical activity of care home
residents [27]. The trial was cluster randomised at the care home level, with
twelve care homes in total randomised equally between treatment as usual (TaU)
and the intervention plus TaU.

Data on several feasibility outcomes were collected. Here, we focus on four:
recruitment (measured in terms of the average number of residents in each care
home who participate in the trial, or average cluster size); adherence (a
binary indicator at the care home level indicating if the intervention was
fully implemented); data completion (a binary indicator for each resident of
successful follow-up at the planned primary outcome time of 12 months); and
potential efficacy (a continuous measure of physical activity at the resident
level). Progression criteria using the traffic light system were pre-specified
for all of these outcomes except potential efficacy, as detailed in Table 2

Table 2: Pre-specified progression criteria used in the original REACH design. Outcome | Red | Amber | Green
---|---|---|---
Recruitment (avg. per care home) | Less than 8 | Between 8 and 10 | At least 10
Adherence | Less than 50% | Between 50 and 75% | At least 75%
Follow-up | Less than 65% | Between 65 and 75% | At least 75%

Denoting the size of the $j$th cluster by $m_{j}$, we assume that cluster
sizes are normally distributed, $m_{j}\sim
N(\mu_{c},\sigma^{2}),j=1,\ldots,2k$. We further assume that the probability
of a participant being followed-up is constant across clusters and arms, and
that the total number follows a binomial distribution $f\sim
Bin(\sum_{j=1}^{2k}m_{j},p_{f})$. The number of care homes which successfully
adhere to the intervention is assumed to binomially distributed, $a\sim
Bin(k,p_{a})$.

The continuous measure of physical activity is expected to be correlated
within care homes. We model this using a random intercept, where the outcome
$y_{ij}$ of resident $i$ in care home $j$ is

$y_{ij}=X_{j}\times Y_{j}\times\mu+u_{j}+\varepsilon_{i}.$ (11)

Here, $X_{j}$ is a binary indicator of care home $j$ being randomised to the
intervention arm, $Y_{j}$ is a binary indicator of care home $j$ successfully
adhering to the intervention, $\mu$ is the mean treatment effect,
$u_{j}\sim\mathcal{N}(0,\sigma_{B}^{2})$ is the random effect for care home
$j$, and $\varepsilon_{i}\sim\mathcal{N}(0,\sigma_{W}^{2})$ is the residual
for resident $i$. We parametrise the model using the intracluster correlation
coefficient, $\rho=\sigma_{B}^{2}/(\sigma_{B}^{2}+\sigma_{W}^{2})$.

The parameters describing average cluster size, follow-up and adherence rates,
and mean treatment effect are of substantive interest when making progression
decisions, giving $\phi=(\mu_{c},p_{f},p_{a},\mu)$. The remainder are nuisance
parameters, $\psi=(\sigma^{2},\rho,\sigma_{W}^{2})$.

###  4.1 Model specification

To begin specifying a model for the REACH trial, we first note that the four
substantive parameters can be divided into two pairs. Firstly, mean cluster
size and follow-up rate relate to the amount of information which a
confirmatory trial will gather. Secondly, potential efficacy and adherence
relate to the effectiveness of the intervention, where effectiveness is
thought of as the effect which will be obtained in practice when the effect of
non-adherence is accounted for. We expect that a degree of trade-off between
adherence and potential efficacy will be acceptable, with a decrease in one
being compensated by an increase in the other. Likewise, low mean cluster size
could be compensated to some extent by higher follow-up rate, and vice versa.

While there may be trade-offs within these pairs of parameters, we do not
expect trade-offs between them. A trial with no effectiveness will be futile
regardless of the amount of information collected, and so should not be
conducted. Similarly, a confirmatory trial should not be conducted if it is
highly unlikely to produce enough information for the research question to be
adequately answered. We therefore consider the sub-spaces of $\Phi$ formed by
these parameter pairs, partition these into hypotheses, and combine these
together. Constructing hypotheses in these two-dimensional spaces is
cognitively simpler than working in the original four dimensional space, not
least because they can be easily illustrated graphically.

Formally, let $\Phi^{i}$ be the sub-space of mean cluster size and follow-up
rate, and $\Phi^{e}$ be that of adherence and potential efficacy. Having
specified hypotheses $\Phi^{i}_{I},\Phi^{e}_{I}$ for $I=R,A,G$, we then have

$\displaystyle\phi\in\begin{cases}\Phi_{R}\text{ if
}\phi^{i}\in\Phi^{i}_{R}\text{ or }\phi^{e}\in\Phi^{e}_{R}\\\ \Phi_{G}\text{
if }\phi^{i}\in\Phi^{i}_{G}\text{ and }\phi^{e}\in\Phi^{e}_{G}\\\
\Phi_{A}\text{ otherwise}.\\\ \end{cases}$ (12)

####  4.1.1 Follow-up and cluster size

Recall that cluster sizes are assumed to be normally distributed with mean
$\mu_{c}$ and variance $\sigma^{2})$. A normal-inverse-gamma prior

$\sigma^{2}\sim\Gamma^{-1}(\alpha_{0},\beta_{0}),~{}\mu_{c}\sim
N(\mu_{0},\sigma^{2}/\nu_{0})$ (13)

is placed on the mean and variance to allow for prior uncertainty in both
parameters. It was anticipated that an average of 8 - 12 residents would be
recruited in each care home. To reflect this prior belief we set the hyper-
parameters to $\mu_{0}=10,\nu_{0}=6,\alpha_{0}=20,\beta_{0}=39$, giving a
prior cluster size of 10 with mean variance 2.05.

For the probability of successful follow-up, $p_{f}$, we take a Beta
distribution with hyper-parameters $\alpha_{0}=22.4,\beta_{0}=9.6$ as the
prior. This gives a prior with a mean of 0.7 and a standard deviation of 0.08.

To partition the parameter space into hypotheses, we first consider the case
where follow-up is perfect, i.e. $p_{f}=1$. Conditional on this, we reason
that a mean cluster size of below 5 should lead to a red decision (stop
development), whereas a size of above 7 should lead to a green decision
(proceed to the main trial). As the probability of successful follow-up
decreases, we suppose that this can be compensated by an increase in mean
cluster size. We assume the nature of this trade-off is linear and decide that
if $p_{f}$ were reduced to 0.8, we would want to have a mean cluster size of
at least 8 to consider decisions $a$ or $g$. We further decide that a follow-
up rate of less than $p_{f}=0.6$ would be critically low, regardless of the
mean cluster size, and should always lead to decision $r$. Similarly, a
follow-up rate of $0.6\leq p_{f}<0.66$ should lead to modification of the
intervention or trial design. Together, these conditions lead to the following
partitioning of the parameter space:

$(p_{f},\mu_{c})\in\begin{cases}\Phi^{i}_{R}\text{ if }p_{f}<0.6\text{ or
}20-15p_{f}>\mu_{c}\\\ \Phi^{i}_{G}\text{ if }p_{f}>0.66\text{ and
}22-15p_{f}<\mu_{c}\\\ \Phi^{i}_{A}\text{ otherwise.}\end{cases}$ (14)

The hypotheses are illustrated in Figure 3 (a). Having specified both the
hypotheses and the prior distribution for these two parameters, we can obtain
prior probabilities of each hypothesis by sampling from the prior and
calculating the proportion of these samples falling into the regions
$\Phi^{i}_{R},\Phi^{i}_{A}$ and $\Phi^{i}_{G}$. We have plotted 1000 samples
from the prior in Figure 3 (a), falling into hypotheses
$\Phi^{i}_{R},\Phi^{i}_{A}$ and $\Phi^{i}_{G}$ in proportions 0.354, 0.517,
0.129 respectively. This demonstrates that there is significant prior
uncertainty regarding the optimal decision, indicating the potential value of
the pilot trial.

(a)  (b)

Figure 3: Marginal hypotheses over parameters for (a) follow-up rate $p_{f}$
and mean cluster size $\mu_{c}$; and (b) adherence rate $p_{a}$ and potential
efficacy $\mu$. Each point is a sample from the joint prior distribution.

####  4.1.2 Adherence and potential efficacy

Having defined priors and hypotheses with respect to cluster size and follow-
up, we now consider adherence and potential efficacy. Recall that the number
of care homes which successfully adhere to the intervention delivery plan is
assumed to be binomially distributed with probability $p_{a}$. We assume that
adherence is absolute in the sense that all residents in a care home which
does not successfully deliver the intervention will not receive any of the
treatment effect. We place a Beta prior on $p_{a}$, with hyper-parameters
$\alpha=28.8$ and $\beta=3.2$ giving a prior mean of 0.9 and a standard
deviation of 0.05.

For the continuous measure of physical activity, we place priors on the mean
effect $\mu$, the intracluster correlation coefficient $\rho$, and the within-
cluster variance $\sigma_{W}^{2}$ in the manner suggested in [22].
Specifically, we choose

$\displaystyle\mu$ $\displaystyle\sim N(0.2,0.25^{2})$ (15)
$\displaystyle\sigma_{W}^{2}$ $\displaystyle\sim\Gamma^{-1}(50,45)$ (16)
$\displaystyle\rho$ $\displaystyle\sim Beta(1.6,30.4).$ (17)

To reflect prior expectation of an ICC around 0.05 but possibly as large as
0.1, the hyperparameters give a prior mean of 0.05 for the ICC with a prior
probability of 0.104 that it will exceed 0.1.

While there is potential for adherence to be improved after the pilot, we
assume there will be little opportunity to improve the potential efficacy of
the intervention. Moreover, we suppose an absolute improvement in adherence of
up to around 0.1 is feasible. To define the hypotheses in this subspace we
first set a minimal level of potential efficacy to be 0.1, and decide that we
would be happy to make decision $g$ at this point if and only if adherence is
perfect. As $p_{a}$ reduces from 1, a corresponding linear increase in
potential efficacy is considered to maintain the overall effectiveness of the
intervention. The rate of substitution for this trade-off is determined to be
approximately 0.57 units of potential efficacy per unit of adherence
probability. We consider an absolute lower limit in adherence of $p_{a}=0.5$,
below which we will always consider decision $r$ to be optimal. Taking these
considerations together, the marginal hypotheses are defined as

$(p_{a},\mu)\in\begin{cases}\Phi^{e}_{R}\text{ if }p_{a}<0.5\text{ or
}0.96-0.57\mu>p_{a}\\\ \Phi^{e}_{G}\text{ if }p_{a}>0.6\text{ and
}1.06-0.57\mu<p_{a}\\\ \Phi^{e}_{A}\text{ otherwise.}\end{cases}$ (18)

The hypotheses are illustrated in Figure 3 (b). Again, a sample of size 1000
from the joint marginal prior distribution $p(p_{a},\mu)$ is also plotted,
falling into hypotheses $\Phi^{e}_{R},\Phi^{e}_{A}$ and $\Phi^{e}_{G}$ in
proportions 0.234, 0.470, 0.296 respectively. As before, this indicates
substantial prior uncertainty regarding the optimal decision and thus supports
the use of a pilot study.

The marginal hypotheses are combined together using equation (12). Considering
the same 1000 samples from the design prior plotted in Figure 3, these now
fall into the regions $\Phi_{R},\Phi_{A}$ and $\Phi_{G}$ in proportions 0.507,
0.458, and 0.035 respectively. Note that the prior probabilities of these
overall hypotheses are quite different to those of the marginal hypotheses. In
particular, there is a considerable increase in the probability that decision
$r$ will be optimal, and a considerable decrease that decision $g$ will be.

###  4.2 Evaluation

####  4.2.1 Weakly informative analysis

We applied the proposed method assuming that a weakly informative joint prior
distribution will be used at the analysis stage222Full details of the weakly
informative prior are given in the appendix.. We took the sample size of the
trial to be $k=6$ clusters per arm. For calculating operating characteristics
we generated $N=10^{4}$ samples from the joint distribution
$p(\theta,x)=p(x|\theta)p_{D}(\theta)$. We analysed each simulated data set
using Stan via the R package rstan [28], in each case generating 5000 samples
in four chains and discarding the first 2500 samples in each to allow for
burn-in, leading to $M=10^{4}$ posterior samples in total. This gave a maximum
Monte Carlo error of approximately 0.005 when estimating a posterior
probability $Pr[\phi\in\Phi_{I}~{}|~{}x]$, which we considered sufficient.
These posterior samples were then used to find the posterior probabilities of
each hypothesis, for each simulated data set.

We evaluated the operating characteristics for a sample of parameters
$(c_{1},c_{2},c_{3})$ as described in Section 2.4. A total of 254 parameter
vectors were evaluated, of which 62 led to operating characteristics which
were worse in every respect than some other vector (i.e. dominated) and were
discarded. The operating characteristics of the non-dominated parameters are
shown in Figure 4. The three operating characteristics are found to be highly
correlated. In particular, changing the parameters to give a lower probability
of discarding a promising intervention ($OC_{2}$) tends to lead to a reduction
in the probability of making an unnecessary adjustment ($OC_{3}$). When
selecting $(c_{1},c_{2})$, the key decision appears to be trading off the
probability of an infeasible trial, ($OC_{1}$), against $OC_{2}$. There is a
very limited opportunity to minimise $OC_{3}$ at the expense of these. For
example, compare points $b$ and $c$ in Figure 4, details of which are given in
Table 3. We see that point $c$ reduces $OC_{3}$ by 0.078 in comparison to
point $b$, but only at the expense of increase in $OC_{1}$ and $OC_{2}$ of
0.13 and 0.145 respectively.

Figure 4: Operating characteristics of the example pilot trial for a range of loss parameter vectors, when a weakly informative analysis prior is used. Table 3: Estimated operating characteristics (with standard errors) of the REACH trial for the three loss parameter vectors highlighted in Figure 4, when a weakly informative analysis prior is used. Costs have been rounded to 2 decimal places; operating characteristics and their errors to 3. Label | $(c_{1},c_{2},c_{3})$ | $OC_{1}$ | $OC_{2}$ | $OC_{3}$
---|---|---|---|---
a | (0.07, 0.9, 0.03) | 0.107 (0.003) | 0.108 (0.003) | 0.232 (0.004)
b | (0.18, 0.58, 0.24) | 0.021 (0.001) | 0.394 (0.005) | 0.08 (0.003)
c | (0.01, 0.29, 0.7) | 0.151 (0.004) | 0.539 (0.005) | 0.002 (0)

We would expect to see a clear relationship between the value of parameters
$c_{1},c_{2},c_{3}$ and the operating characteristics they relate to. We
explore this in Figure 5 with scatter plots of each parameter against each
operating characteristic. The results show that there is indeed a strong
relationship between the loss assigned to discarding a promising intervention,
$c_{2}$, and the probability that this event will occur, $OC_{2}$ (see centre
plot). Moreover, $c_{2}$ also seems to be the main determinant of operating
characteristics $OC_{1}$ and $OC_{3}$. The implication is that once the
$c_{2}\in[0,1]$ has been chosen, the operating characteristics of the trial
depend only weakly on the way in which the remaining $1-c_{3}$ is allocated to
$c_{1}$ and $c_{3}$. This appears to be due to the fact that, regardless of
how errors are weighted, the way we have defined our prior distributions and
hypotheses means we are much more likely to make the error of discarding a
promising intervention than the other types of error. The cost we assign to
this error is therefore more influential on the overall operating
characteristics than the other costs.

Figure 5: Relationships between the three loss parameters ($x$ axes) and
resulting operating characteristics ($y$ axes).

To illustrate the effect of varying sample size in the REACH trial, we set the
loss function parameters to that of point $a$ in Figure 4 and Table 3,
$(c_{1},c_{2},c_{3})=(0.07,0.9,0.03)$. We then estimated the operating
characteristics obtained for $k=6,12,18$ clusters per arm. Note that we
considered only three choices of sample size due to the significant
computational burden of each evaluation. The results are plotted in Figure 6.
Increasing the sample size appears to have little effect on $OC_{1}$ and
$OC_{3}$, while leading to a decrease in $OC_{2}$, the probability of
discarding a promising intervention. This behaviour reflects the priorities
encoded by the costs parameter, where $c_{2}=0.9$.

Figure 6: Operating characteristics of the REACH trial for per-arm sample
sizes $k=6,12,18$ and setting $(c_{1},c_{2},c_{3})=(0.069,0.116,0.815)$. Error
bars denote 95% confidence intervals. All points have been adjusted
horizontally to avoid overlap.

####  4.2.2 Incorporating subjective priors

Rather than use weakly or non-informative priors when analysing the pilot
data, we may instead want to make use of the (subjective) elicited knowledge
of parameter values described in the design prior $p_{D}(\theta)$.
Anticipating criticisms of a fully subjective analysis, we can envisage two
particular cases where this might be appropriate. Firstly, using the
components of the design prior which describe the nuisance parameters $\psi$
while maintaining weakly informative priors on substantive parameters $\phi$.
Secondly, when very little data on a specific substantive parameter is going
to be collected in the pilot, using the informative design prior for that
parameter could substantially improve operating characteristics.

We replicated the above analysis for these two scenarios. For the second, we
used informative priors for all nuisance parameters and for the probability of
adherence, $p_{a}$. Recall that this is informed by a binary indicator at the
care home level and only in the intervention arm, and will therefore have very
little pilot data bearing on it. For each case we used the same $N$ samples of
parameters and pilot data which were used in the weakly informative case,
repeating the Bayesian analysis using the appropriate analysis prior and
obtaining estimated posterior probabilities $p_{R},p_{A}$ and $p_{G}$ as
before. These were used in conjunction with the same set of loss parameter
vectors $\mathcal{C}$ to obtain corresponding operating characteristics.

For brevity we will refer to the three cases as Weakly Informative (WI),
Informative Nuisance (IN), and Informative Nuisance and Adherence (INA).
Comparing the operating characteristics of cases WI and IN, we found very
little difference (further details are provided in the appendix). When we
contrast cases WI and INA, however, there is a clear distinction. Using the
INA analysis prior will lead to larger probabilities of an infeasible trial
($OC_{1}$) and of unnecessary adjustment ($OC_{2}$), while reducing the
probability of discarding a promising intervention ($OC_{3}$), for almost all
loss parameters. The expected loss is always lower for the INA analysis than
for WI, as we would expect.

Figure 7: Operating characteristics and expected utilities for weakly (WI) and
partially informative (INA)

##  5 Discussion

When deciding if and how a definitive RCT of a complex intervention should be
conducted, and basing this decision on an analysis of data from a small pilot
trial, there is a risk we will inadvertently make the wrong choice. A Bayesian
analysis of pilot data followed by decision making based on a loss function
can help ensure this risk is minimised. The expected results of such a pilot
can be evaluated through simulation at the design stage, producing operating
characteristics which help us understand the potential for the pilot to lead
to better decision making. These evaluations can in turn be used to find the
loss function which leads to the most desirable operating characteristics, and
to inform the choice of sample size.

Our proposal has been motivated by some salient characteristics of complex
intervention pilot trials, and offers several potential benefits over standard
pilot trial design and analysis techniques. The Bayesian approach to analysis
means that complex multi-level models can be used to describe the data, even
when the sample size is small. In contrast to the usual application of
independent progression criteria for several parameters of interest, we
provide a way for preferential relationships between parameters to be
articulated and used when making decisions. Using a subjective prior
distribution on unknown parameters at the design stage allows both our
knowledge and our uncertainty to be fully expressed, meaning we can leverage
external information whilst also avoiding decisions which are highly sensitive
to imprecise point estimates.

Our proposed design is related to the literature on assurance calculations for
clinical trials [17], applying the idea of using unconditional event
probabilities as operating characteristics to the pilot trial setting. In
doing so we have shown how assurances can be defined for multiple substantive
parameters with trade-offs between them, and with respect to the ‘traffic
light’ red/amber/green decision structure commonly found in pilot trials. The
multi-objective optimisation framework we have used to inform trial design
allows the decision-maker to explicitly consider the different trade-offs
between operating characteristics which are available, and select that which
best reflects their own preferences. A similar approach has been taken in the
context of phase II trials using the statistical concept of admissible designs
[29, 30]. This can be contrasted with the conventional and much criticised
approach common in the frequentist context, where arbitrary constraints are
placed on type I and II error rates in order to define a single optimal design
[31].

The benefits brought by the Bayesian approach must be set against the
challenges it brings, particularly in terms of computation time and
implementation. In terms of the latter, we are required to specify a joint
prior distribution over the parameters $\theta$ and a partitioning of the
parameter space into the three hypotheses. The specification of the prior
distribution may be a challenging and time-consuming task. Although some
relevant data relating to similar contexts may be available, for example in
systematic reviews or observational studies, expert opinion may still be
required to articulate the relevance of such data to the problem at hand. When
no data are available, which is not unlikely given the early phase nature of
pilot studies, expert opinion will be the only source of information. Although
potentially challenging, many examples describing successful practical
applications of elicitation for clinical trial design are available [20, 18,
23], as are tools for its conduct such as the Sheffield Elicitation Framework
(SHELF) [24]. Dividing the parameter space into three hypotheses may also
prove challenging in practice, particularly when trade-offs between more than
two parameters are to be elicited. There is a need for methodological research
investigating how methods for multi-attribute preference elicitation, such as
those set out in [26], can be applied in this context.

The computational burden of the proposed method is significant, particularly
when the model is too complex to allow a conjugate analysis to be used when
sampling from the posterior distribution. We have used a nested Monte Carlo
sampling scheme to estimate operating characteristics, as seen elsewhere [19,
17, 32]. One potential approach to improve efficiency is to use non-parametric
regression to predict the expected losses of Equation (3) based on some
simulated data, thus bypassing the need to undertake a full MCMC analysis for
each of the $N$ samples in the outer loop. This approach has been shown to be
successful in the context of expected value of information calculations [33,
34]. The computational difficulties will be particularly pertinent when using
our approach to determine sample size, as several evaluations of different
sample size choices will be required. If the choice of sample size can be
framed as an optimisation problem, methods for efficient global optimisation
of computationally expensive functions such as those described in [35, 36] may
be useful [16]. Alternatively, one of several rules-of-thumb for choosing
pilot sample size [3, 9, 11, 13] could be used, with the resulting operating
characteristics evaluated using the proposed method.

We have defined our procedure in terms of a loss function, where the decision
making following the pilot will minimise the expected loss. However, the
piecewise constant loss function we have proposed may not adequately represent
the preferences of the decision maker. For example, we may object to the loss
associated with discarding a promising intervention being independent of
exactly how effective the intervention is. An alternative is to try to define
a richer representation of the loss function through direct elicitation of the
decision makers preferences under uncertainty [25], leading to a fully
decision-theoretic approach to design and analysis [37]. However, as
previously noted by others [38, 39, 40], implementation of these approaches
has been limited in practice and this may be indicative of their feasibility.

The proposed method could be extended in several ways. More operating
characteristics could be defined and used in design optimisation, more
complicated trade-off relationships between multiple parameters could be
addressed, or the hypotheses could be expanded to include nuisance parameters
which would be used as part of the sample size calculation in the main RCT. A
particularly interesting avenue for future research is to consider how to
model post-pilot trial actions in more detail. For example, while we allow for
the possibility of making an ‘amber’ decision, indicating that modifications
to the intervention or trial design should be made, we do not model what that
decision will actually look like and how it should relate to the observed
pilot data. Methodology for jointly modelling a pilot and subsequent main RCT
in this manner could be informed by developments for designing phase II/III
programs in the drug setting [41, 42, 43, 44].

### Acknowledgements

We would like to thank Alex Wright-Hughes, Robert Cicero, and the TIGA-CUB and
REACH trial teams for discussions which helped shape the scope of this paper.

### Data availability statement

All simulated data used in this manuscript, together with the code used to
generate it, is available at https://github.com/DTWilson/Bayesian_pilot.

### Funding

This work was supported by the Medical Research Council under Grant
MR/N015444/1 to D.T.W. and Grant MC_UU_00002/6 to J.M.S.W.

## References

  * [1] Peter Craig, Paul Dieppe, Sally Macintyre, Susan Michie, Irwin Nazareth, and Mark Petticrew.  Developing and evaluating complex interventions: the new medical research council guidance.  BMJ: British Medical Journal, 337, 9 2008.
  * [2] Sandra M. Eldridge, Gillian A. Lancaster, Michael J. Campbell, Lehana Thabane, Sally Hopewell, Claire L. Coleman, and Christine M. Bond.  Defining feasibility and pilot studies in preparation for randomised controlled trials: Development of a conceptual framework.  PLOS ONE, 11(3):e0150205, mar 2016.
  * [3] Gillian A. Lancaster, Susanna Dodd, and Paula R. Williamson.  Design and analysis of pilot studies: recommendations for good practice.  Journal of Evaluation in Clinical Practice, 10(2):307–312, 2004\.
  * [4] National Institute for Health Research.  Research for patient benefit (rfpb) programme guidance on applying for feasibility studies, 2017.
  * [5] Sandra M Eldridge, Claire L Chan, Michael J Campbell, Christine M Bond, Sally Hopewell, Lehana Thabane, and Gillian A Lancaster.  CONSORT 2010 statement: extension to randomised pilot and feasibility trials.  BMJ, page i5239, oct 2016.
  * [6] Kerry N L Avery, Paula R Williamson, Carrol Gamble, Elaine O’Connell Francischetto, Chris Metcalfe, Peter Davidson, Hywel Williams, and Jane M Blazeby.  Informing efficient randomised controlled trials: exploration of challenges in developing progression criteria for internal pilot studies.  BMJ Open, 7(2):e013537, feb 2017.
  * [7] Lisa V Hampson, Paula R Williamson, Martin J Wilby, and Thomas Jaki.  A framework for prospectively defining progression rules for internal pilot studies monitoring recruitment.  Statistical Methods in Medical Research, 0(0):0962280217708906, 2017\.  PMID: 28589752.
  * [8] Richard H. Browne.  On the use of a pilot sample for sample size determination.  Statistics in Medicine, 14(17):1933–1940, 1995.
  * [9] Steven A. Julious.  Sample size of 12 per group rule of thumb for a pilot study.  Pharmaceutical Statistics, 4(4):287–291, 2005.
  * [10] Julius Sim and Martyn Lewis.  The size of a pilot study for a clinical trial should be calculated in relation to considerations of precision and efficiency.  Journal of Clinical Epidemiology, 65(3):301–308, mar 2012.
  * [11] M Teare, Munyaradzi Dimairo, Neil Shephard, Alex Hayman, Amy Whitehead, and Stephen Walters.  Sample size requirements to estimate key design parameters from external pilot randomised controlled trials: a simulation study.  Trials, 15(1):264, 2014.
  * [12] Sandra M Eldridge, Ceire E Costelloe, Brennan C Kahan, Gillian A Lancaster, and Sally M Kerry.  How big should the pilot study for my cluster randomised trial be?  Statistical Methods in Medical Research, 2015.
  * [13] Amy L Whitehead, Steven A Julious, Cindy L Cooper, and Michael J Campbell.  Estimating the sample size for a pilot randomised trial to minimise the overall trial sample size for the external pilot and main trial for a continuous outcome variable.  Statistical Methods in Medical Research, 2015.
  * [14] Wolfgang Viechtbauer, Luc Smits, Daniel Kotz, Luc Budé, Mark Spigt, Jan Serroyen, and Rik Crutzen.  A simple formula for the calculation of sample size in pilot studies.  Journal of Clinical Epidemiology, 68(11):1375–1379, nov 2015.
  * [15] Cindy L Cooper, Amy Whitehead, Edward Pottrill, Steven A Julious, and Stephen J Walters.  Are pilot trials useful for predicting randomisation and attrition rates in definitive studies: A review of publicly funded trials.  Clinical Trials, 0(0):1740774517752113, 2018.  PMID: 29361833.
  * [16] D. T. Wilson, R. E. Walwyn, J. Brown, A. J. Farrin, and S. R. Brown.  Statistical challenges in assessing potential efficacy of complex interventions in pilot or feasibility studies.  Statistical Methods in Medical Research, 25(3):997–1009, jun 2015\.
  * [17] Anthony O’Hagan, John W. Stevens, and Michael J. Campbell.  Assurance in clinical trial design.  Pharmaceutical Statistics, 4(3):187–201, 2005.
  * [18] Adam Crisp, Sam Miller, Douglas Thompson, and Nicky Best.  Practical experiences of adopting assurance as a quantitative framework to support decision making in drug development.  Pharmaceutical Statistics, 0(0), 2018.
  * [19] Fei Wang and Alan E. Gelfand.  A simulation-based approach to bayesian sample size determination for performance under a given model and for separating models.  Statistical Science, 17(2):pp. 193–208, 2002.
  * [20] Rosalind J. Walley, Claire L. Smith, Jeremy D. Gale, and Phil Woodward.  Advantages of a wholly bayesian approach to assessing efficacy in early drug development: a case study.  Pharmaceutical Statistics, 14(3):205–215, apr 2015.
  * [21] Andrew Gelman, Aleks Jakulin, Maria Grazia Pittau, and Yu-Sung Su.  A weakly informative default prior distribution for logistic and other regression models.  The Annals of Applied Statistics, 2(4):1360–1383, dec 2008.
  * [22] David J. Spiegelhalter.  Bayesian methods for cluster randomized trials with continuous responses.  Statistics in Medicine, 20(3):435–452, 2001.
  * [23] Nigel Dallow, Nicky Best, and Timothy H Montague.  Better decision making in drug development through adoption of formal prior elicitation.  Pharmaceutical Statistics, 0(0), 2018.
  * [24] Anthony O’Hagan, Caitlin E. Buck, Alireza Daneshkhah, J. Richard Eiser, Paul H. Garthwaite, David J. Jenkinson, Jeremy E. Oakley, and Tim Rakow.  Uncertain Judgements: Eliciting Experts’ Probabilities.  John Wiley and Sons, 2006.
  * [25] Simon French and David Rios Insua.  Statistical Decision Theory.  Number 9 in Kendall’s Library of Statistics. Oxford University Press, 2000\.
  * [26] Ralph L. Keeney and Howard Raiffa.  Decisions with multiple objectives: preferences and value tradeoffs.  John Wiley & Sons, 1976.
  * [27] Anne Forster, , Jennifer Airlie, Karen Birch, Robert Cicero, Bonnie Cundill, Alison Ellwood, Mary Godfrey, Liz Graham, John Green, Claire Hulme, Rebecca Lawton, Vicki McLellan, Nicola McMaster, and Amanda Farrin.  Research exploring physical activity in care homes (REACH): study protocol for a randomised controlled trial.  Trials, 18(1), apr 2017.
  * [28] Stan Development Team.  RStan: the R interface to Stan, 2016.  R package version 2.14.1.
  * [29] Sin-Ho Jung, Taiyeong Lee, Kyung Mann Kim, and Stephen L. George.  Admissible two-stage designs for phase II cancer clinical trials.  Statistics in Medicine, 23(4):561–569, 2004.
  * [30] Adrian P. Mander, James M.S. Wason, Michael J. Sweeting, and Simon G. Thompson.  Admissible two-stage designs for phase II cancer clinical trials that incorporate the expected sample size under the alternative hypothesis.  Pharmaceutical Statistics, 11(2):91–96, 2012.
  * [31] Peter Bacchetti.  Current sample size conventions: Flaws, harms, and alternatives.  BMC Medicine, 8(1):17, Mar 2010.
  * [32] Alexander J. Sutton, Nicola J. Cooper, David R. Jones, Paul C. Lambert, John R. Thompson, and Keith R. Abrams.  Evidence-based sample size calculations based upon updated meta-analysis.  Statistics in Medicine, 26(12):2479–2500, 2007.
  * [33] Mark Strong, Jeremy E. Oakley, and Alan Brennan.  Estimating multiparameter partial expected value of perfect information from a probabilistic sensitivity analysis sample.  Medical Decision Making, 34(3):311–326, apr 2014.
  * [34] Mark Strong, Jeremy E Oakley, Alan Brennan, and Penny Breeze.  Estimating the expected value of sample information using the probabilistic sensitivity analysis sample: a fast, nonparametric regression-based method.  Medical Decision Making, 35(5):570–583, 2015.
  * [35] Donald R. Jones.  A taxonomy of global optimization methods based on response surfaces.  Journal of Global Optimization, 21(4):345–383, 2001.
  * [36] Olivier Roustant, David Ginsbourger, and Yves Deville.  DiceKriging, DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization.  Journal of Statistical Software, 51(1):1–55, 2012.
  * [37] Dennis V. Lindley.  The choice of sample size.  Journal of the Royal Statistical Society: Series D (The Statistician), 46(2):129–138, 1997.
  * [38] Lawrence Joseph and David B. Wolfson.  Interval-based versus decision theoretic criteria for the choice of sample size.  Journal of the Royal Statistical Society: Series D (The Statistician), 46(2):145–149, 1997.
  * [39] Peter Bacchetti, Charles E. McCulloch, and Mark R. Segal.  Simple, defensible sample sizes based on cost efficiency.  Biometrics, 64(2):577–585, jun 2008.
  * [40] John Whitehead, Elsa Valdés-Márquez, Patrick Johnson, and Gordon Graham.  Bayesian sample size for exploratory clinical trials incorporating historical data.  Statistics in Medicine, 27(13):2307–2327, 2008.
  * [41] Nigel Stallard.  Optimal sample sizes for phase II clinical trials and pilot studies.  Statistics in Medicine, 31(11-12):1031–1042, 2012.
  * [42] James M. S. Wason, Thomas Jaki, and Nigel Stallard.  Planning multi-arm screening studies within the context of a drug developmentâ€‰program.  Statistics in Medicine, 32(20):3424–3435, 2013.
  * [43] Heiko Götte, Armin Schüler, Marietta Kirchner, and Meinhard Kieser.  Sample size planning for phase II trials based on success probabilities for phase III.  Pharmaceutical Statistics, 14(6):515–524, sep 2015.
  * [44] Marietta Kirchner, Meinhard Kieser, Heiko Götte, and Armin Schüler.  Utility-based optimization of phase II/III programs.  Statist. Med., 35(2):305–316, aug 2015.