# Translating predictive distributions into informative priors

Andrew A. Manderson111MRC Biostatistics Unit, the University of Cambridge and
The Alan Turing Institute<EMAIL_ADDRESS>Robert J. B.
Goudie222MRC Biostatistics Unit, the University of Cambridge.
<EMAIL_ADDRESS>
(June, 2023)

###### Abstract

When complex Bayesian models exhibit implausible behaviour, one solution is to
assemble available information into an informative prior. Challenges arise as
prior information is often only available for the observable quantity, or some
model-derived marginal quantity, rather than directly pertaining to the
natural parameters in our model. We propose a method for translating available
prior information, in the form of an elicited distribution for the observable
or model-derived marginal quantity, into an informative joint prior. Our
approach proceeds given a parametric class of prior distributions with as yet
undetermined hyperparameters, and minimises the difference between the
supplied elicited distribution and corresponding prior predictive
distribution. We employ a global, multi-stage Bayesian optimisation procedure
to locate optimal values for the hyperparameters. Three examples illustrate
our approach: a cure-fraction survival model, where censoring implies that the
observable quantity is _a priori_ a mixed discrete/continuous quantity; a
setting in which prior information pertains to $R^{2}$ – a model-derived
quantity; and a nonlinear regression model.

##  1 Introduction

A key asset to the Bayesian paradigm is the conceptual ease with which prior
information is incorporated into models. For complex, nonlinear,
overparameterised, or otherwise partially identified (Gustafson, 2015) models,
including such prior information is essential to exclude model behaviours that
conflict with reality, and/or known qualities of the phenomena being modelled.
Prior information can also improve computation of estimates of the posterior,
making otherwise unusable models suitable for inference. However, it is for
precisely the models for which prior information is so important that setting
appropriate priors for parameters is hardest.

We consider in this paper the task of forming such appropriate informative
priors. We distinguish two tasks: predictive elicitation; and the subsequent
translation into a prior for a given model. Predictive elicitation is an
approach in which elicitation (O’Hagan _et al._ , 2006; Falconer _et al._ ,
2021; Low Choy, 2012) proceeds via predictive distributions, often for
observable quantities, and is thought to be the most reliable and available
form of information (Kadane and Wolfson, 1998). Predictive elicitation is also
model-agnostic, meaning the complex, time-consuming process of elicitation
does not need to be repeated for each variant of a model. A recent,
comprehensive review of elicitation is undertaken by Mikkola _et al._ (2021).

Translation is the process of using information from predictive elicitation to
set the corresponding informative prior for a model. Translation, as a
distinct step in the prior specification process, has received less attention
than elicitation. A simple predictive approach is to directly model the
elicited, observable quantity. This direct approach requires no translation.
For example, the Bayesian quantile-parameterised likelihood (Hadlock, 2017;
Keelin and Powley, 2011) approach of Perepolkin _et al._ (2021) involves
updating directly elicited information in light of observations. Such direct
approaches are currently only feasible for simple models with no latent
structure. For models with simple latent structure, it is sometimes possible
to elicit information about an invertible function of the parameters (e.g.
Chaloner _et al._ , 1993). In these cases it is possible to analytically
translate the elicited information into a prior for the parameters.
Translation is also clear for conjugate models (Percy, 2002), as if we specify
the target prior predictive distribution using the conjugate distribution,
then the prior predictive distribution determines the values for the
hyperparameters of the prior (related to the idea of “reverse Bayes” in Held
_et al._ , 2022). Translation, however, is unclear in general for nonconjugate
models (Gribok _et al._ , 2004), although techniques for specific models with
specific latent structure are numerous, and include linear regression
(Ibrahim, 1997), logistic regression (Bedrick _et al._ , 1997; Chen _et al._ ,
1999), Cox models (Chaloner _et al._ , 1993), contingency table analyses
(Good, 1967), hierarchical models (Hem, 2021, noting that the space of
possible hierarchical models is vast), and autoregressive time-series models
(Jarociński and Marcet, 2019). Nevertheless, a model-agnostic approach with a
corresponding generic implementation would be preferable, as noted by Gelman
_et al._ (2017) and Mikkola _et al._ (2021).

Our approach to translation builds on the idea of predictive checks Gabry _et
al._ , 2019; Gelman _et al._ , 2017; Box, 1980; Winkler, 1967, the
“hypothetical future samples” of; van Zundert _et al._ , 2022, which are an
often recommended (Gelman _et al._ , 2020; van Zundert _et al._ , 2022) tool
in assessing the concordance of the prior predictive distribution and the
elicited predictive information or elicited data distribution. If concordance
between these two distributions is low, then the Bayesian workflow (Gelman _et
al._ , 2020) proceeds by adjusting the prior to better match the prior
predictive distribution to the elicited information. However, for many classes
of models, manually adjusting the prior in this manner is infeasible; the
complexity required to describe the phenomena of interest muddies the
relationship between the prior and the distribution of the observables, and so
a more automated method is required. Wang _et al._ (2018) and Thomas _et al._
(2020) have proposed approaches in which either regions of observable space or
specific realisations are labelled as plausible or implausible by experts, and
then an suitable prior formed that accounts for this information, via history
matching and a “human in the loop” processes driven by a Gaussian process
model respectively. Albert _et al._ (2012) propose a supra-Bayesian approach
intended for multiple experts, in which a Bayesian model is formed for
quantiles or probabilities elicited from the experts. Another approach, and
the closest in motivation and methodology to ours, is Hartmann _et al._ 2020;
da Silva _et al._ , 2019, which is partly inspired by, which employs a
Dirichlet likelihood for elicited predictive quantiles to handle both
elicitation and translation.

In this paper we develop a method, and software package, for translation: that
is, for constructing an informative prior distribution for model parameters
that results in a desired _target_ prior predictive distribution. Using as the
starting point for our method a target prior predictive distribution means
that our approach is applicable regardless of how the target is elicited, and
means that uncertainty is directly and intuitively represented as a
distribution. We define a suitable loss function between the prior predictive
distribution, given specific values for the hyperparameters of the prior, and
this target distribution. The loss function is intentionally generic and
permits data that are discrete, continuous, or a mixture thereof. We minimise
this loss via a generic, simulation-based, global optimisation process to
locate optimal hyperparameters. We highlight throughout that solutions to this
optimisation problem are rarely unique, so to regularise the problem we adopt
a multiple objective approach. The global optimisation approach is also
selected with generality in mind, rendering our method applicable to models
where derivative information is unavailable. We make our method available in
an R package pbbo333The release used in this paper is available at
https://doi.org/10.5281/zenodo.7736707. (R Core Team, 2022). Our method is
illustrated in three challenging, moderate dimension problems; a nonlinear
regression model, a model using predictive information on a model-derived
quantity – a situation not explicitly covered by prior predictive elicitation;
and a cure fraction survival model.

##  2 Translating elicited prior predictive distributions

In this section we propose a set of desirable properties for any translation
method, and our mathematical framework and optimisation strategy that together
constitute our methodology. We aim to minimise the difference between the
prior predictive distribution and our elicited target distribution, whilst
also promoting the marginal variance of the model’s parameters. Satisfying the
first of these goals produces priors faithful to the supplied information; the
second promotes uniqueness and replicability – we elaborate on the precise
meaning of these properties momentarily. We adopt a multi-objective
optimisation approach to locate priors satisfying these requirements.

###  2.1 Desiderata

We postulate three key properties that we would like a translation method to
satisfy: faithfulness, uniqueness, and replicability.

_Faithfulness_ We consider a prior faithful if it accurately encodes the
target data distribution provided by the elicitation subject. Faithfulness is
a property of both the model, as it must be possible for the model to
represent the information, and the procedure employed to obtain the prior.
Especially with simple models and prior structures, not all target prior
predictive distributions can be encoded.

_Uniqueness_ In a complex model there may be many prior distributions that
imply the same prior predictive distribution. These prior distributions will
all be equally faithful. Should uniqueness be desired – and it seems a
reasonable enough desiderata most of the time – we must distinguish between
priors based on other properties. In Section 2.4 we propose to distinguish
priors based on their marginal standard deviations, but other properties could
be easily incorporated into our approach.

_Replicability_ We call a procedure/method consistent if it obtains the same,
or very similar, prior across independent replications, given the same target.
This property is particularly important to assess for methods, like ours, that
make use of simulation-based or otherwise stochastic estimates. Global,
gradient-free optimisers, when applied to noisy loss functions, offer no
guarantee of finding the global minimum in finite time (Liberti, 2008; Mullen,
2014). We assess replicability empirically in all our examples.

These properties are partly inspired by other works in the prior elicitation
and specification literature. Faithfulness is closely related to Johnson _et
al._ (2010a)’s definition of _validity_ (see also Johnson _et al._ , 2010b)
and O’Hagan _et al._ (2006)’s use of _faithful_ in Chapter 8 (and throughout
the book). However, their concerns are specific to the elicitation process –
do the quantities elicited represent what the experts believe? – and not to
the subsequent translational step. Our conception of uniqueness and the need
for regularisation is noted by da Silva _et al._ (2019) and Stefan _et al._
(2022), and is similar to the notion of model _sloppiness_ of Gutenkunst _et
al._ (2007).

We have introduced our desiderata in what we believe to be their order of
importance. Firstly, without faithfulness the procedure has not achieved the
main aim of translating our knowledge into the prior distribution.
Subsequently, given a suite of faithful priors, regularising the problem until
it is unique allows us to select one in a clear and replicable way. Such
uniqueness inducing regularisation schemes often improve a procedure’s
replicability. Replicability ultimately also relies on the empirical behaviour
of the procedure when applied to the model of interest. We note that the
desiderata are not binary, and at times we may wish to sacrifice some amount
of the latter properties for improved faithfulness. We also envisage settings
where sacrificing some model flexibility, and thus faithfulness, for a marked
increase in replicability increases the usefulness or persuasiveness of a
model.

###  2.2 The target predictive distribution $\text{T}(Y)$

Our methodology assumes a target predictive distribution (a cumulative
distribution function, CDF) for the observable quantity, $\text{T}(Y)$, has
been chosen. The target distribution can be chosen in any manner, but we
recommend using predictive elicitation (Kadane and Wolfson, 1998). In brief,
such an elicitation proceeds by querying experts about the observable quantity
at a small number quantiles, then fitting an appropriate parametric
distribution to the elicited values (see Chapter 6 of O’Hagan _et al._ ,
2006). We assume a (mixture of) standard distributions can describe the target
predictive distribution function $\text{T}(Y)$, and that we can draw samples
from this distribution.

We often wish to elicit information about the observable quantity $Y$
conditional on some known values of a covariate. For example, when using the
linear model $Y=X\beta+\varepsilon$ we may elicit information about $Y$ at a
fixed set of values for $X$. Further suppose $X$ is an experimental design
specified before collecting observations of $Y$, or comprises observational
covariates whose values are known prior to model construction. In such
settings we can elicit $r=1,\ldots,R$ conditional target distributions
$\text{T}(Y\mid X_{r})$.

We elect to describe our methodology in this covariate-specific setting, as it
readily reduces to the covariate-independent case.

###  2.3 Total predictive discrepancy (primary objective)

Consider a joint model for observables $Y\in\mathcal{Y}\subseteq{}\mathbb{R}$
and parameters $\theta\in\Theta\subseteq\mathbb{R}^{Q}$, given hyperparameters
$\lambda\in\Lambda\subset\mathbb{R}^{L}$ and covariates
$X\in\mathcal{X}\subseteq\mathbb{R}^{C}$. This joint model has CDF
$\text{P}(Y,\theta\mid\lambda,X)$ and prior predictive CDF
$\text{P}(Y\mid\lambda,X)$ for $Y$. Choosing the CDF as the basis for our
methodology enables us to be agnostic to whether the observable is continuous,
discrete, or a mixture thereof. We will use _distribution_ to refer to the CDF
of a stochastic quantity, and _density_ to refer to the corresponding
probability density function (where it exists).

Further suppose the target CDF $\text{T}(Y\mid X_{r})$ has been elicited at
$R$ values of the covariate vector denoted $\\{X_{r}\\}_{r=1}^{R}$, which we
stack in the covariate matrix $\boldsymbol{X}=\left[X_{1}^{\top}\cdots
X_{R}^{\top}\right]\in\boldsymbol{\mathcal{X}}\subseteq\mathbb{R}^{R\times
C}$. We assume that each target $\text{T}(Y\mid X_{r})$ has identical support
to $\text{P}(Y\mid\lambda,X_{r})$. Lastly, it will be convenient to denote
$\text{T}(Y\mid\boldsymbol{X})=\prod_{r=1}^{R}\text{T}(Y\mid X_{r})$, with
$\text{P}(Y\mid\lambda,\boldsymbol{X})$ and
$\text{P}(\theta\mid\lambda,\boldsymbol{X})$ defined analogously.

We quantify the difference between the prior predictive and target by the
_covariate-specific predictive discrepancy_ , which we define to be

$\tilde{D}(\lambda\mid\boldsymbol{X})=\frac{1}{R}\sum_{r=1}^{R}\int
d(\text{P}(Y\mid\lambda,X_{r}),\text{T}(Y\mid X_{r}))\text{d}\text{T}(Y\mid
X_{r}),$ (1)

for some discrepancy function $d(\cdot,\cdot)$. The Riemann-Stieltjes integral
in Equation 1 means this definition applies when $\mathcal{Y}$ is continuous,
discrete, or a mixture thereof. Minimising Equation (1) admits the optimal
hyperparameter
$\lambda^{*}=\min_{\lambda\in\Lambda}\tilde{D}(\lambda\mid\boldsymbol{X})$.
The covariate-independent equivalent $\tilde{D}(\lambda)$ is obtained by
setting $R=1$ and ignoring all conditioning on $X_{r}$ in Equation (1).

The discrepancy function $d(\cdot,\cdot)$ takes two CDFs as its arguments.
Inspired by the Cramér-von Mises (von Mises, 1947) and Anderson-Darling
(Anderson and Darling, 1952) distributional tests we define, for arbitrary
CDFs $\text{M}(Y)$ and $\text{P}(Y)$, two options for the discrepancy
function,

$d^{\text{CvM}}(\text{M}(Y),\text{P}(Y))=\left(\text{M}(Y)-\text{P}(Y\right))^{2},\quad
d^{\text{AD}}(\text{M}(Y),\text{P}(Y))=\frac{\left(\text{M}(Y)-\text{P}(Y\right))^{2}}{\text{P}(Y)(1-\text{P}(Y))}.$
(2)

Both discrepancies are proper scoring rules (Gneiting and Raftery, 2007) as
they are minimised iff $\text{M}(Y)=\text{P}(Y)$ for all $Y\in\mathcal{Y}$.
Supposing $\text{P}(Y\mid\lambda,X_{r})$ is flexible enough to exactly match
$\text{T}(Y\mid X_{r})$ for some $\lambda^{*}$, then both discrepancies will
yield the same $\lambda^{*}$. Differences arise when
$\text{P}(Y\mid\lambda,X_{r})$ is insufficiently flexible, in which case
Anderson-Darling discrepancy $d^{\text{AD}}$ places more emphasis on matching
the tails of two CDFs under consideration. Furthermore, we will have to resort
to a finite-sample approximation to Equation (1) (which we detail
momentarily), and the Anderson-Darling discrepancy is more challenging to
accurately compute.

###  2.4 Regularising estimates of $\lambda^{*}$ by promoting the marginal
standard deviation (secondary objective)

The optimisation problem of minimising Equation (1) is often underspecified.
Specifically, there are many optimal values $\lambda^{*}$ that yield values of
$\tilde{D}(\lambda^{*}\mid\boldsymbol{X})$ that are practically
indistinguishable (noted by da Silva _et al._ , 2019) but with immensly
differing prior distributions $\text{P}(\theta\mid\lambda^{*},\boldsymbol{X})$
and corresponding marginals for components of $\theta$. In terms of our
desiderata, there are many equally faithful priors (which immediately implies
a lack of uniqueness), thus we have an optimisation problem with solutions
that are difficult to replicate due to nonuniqueness. This is not surprising
because we are providing information only on $Y$, which is typically of lower
dimension than $\theta$. A particularly challenging case for uniqueness is in
models with additive noise forms (such as Equation (17)). In this case it is
challenging to avoid attributing all the variability to the noise term if the
prior for the noise is to be specified; in this case it will generally be
necessary to fix a prior for the variance using knowledge of the measurement
process.

To handle more general cases of lack of uniqueness, we seek to encode the
following principle into our methodology: given two estimates of $\lambda^{*}$
which have equivalent values of $\tilde{D}(\lambda^{*}\mid\boldsymbol{X})$, we
prefer the one with the larger variance for
$\text{P}(\theta\mid\lambda^{*},\boldsymbol{X})$. This preference induces less
of a restriction on the possible values of $\theta$ in the posterior.

We make use of this principle by adopting a multi-objective approach to prior
construction and, therefore, now derive a suitable mathematical quantity
measuring the variability of $\theta$. There are numerous suitable functions
measuring such variability, and our methodology is agnostic to the particular
functional form. Most generally, we define the secondary objective
$\tilde{N}(\lambda\mid\boldsymbol{X})$ as comprising any such suitable
function $n(\theta)$ with

$\tilde{N}(\lambda\mid\boldsymbol{X})=\int
n(\theta)\,\text{d}\text{P}(\theta\mid\lambda,\boldsymbol{X}).$ (3)

In this paper we consider only one form for $n(\theta)$, and so hereafter the
second objective, which we also seek to minimise, is always

$\tilde{N}(\lambda\mid\boldsymbol{X})=-\frac{1}{Q}\sum_{q=1}^{Q}\log\left(\text{SD}_{\,\text{P}(\theta_{q}\mid\lambda,\boldsymbol{X})}\left[\theta_{q}\right]\right),$
(4)

where $\text{SD}_{\text{P}(Z)}[Z]$ is the standard deviation of $Z$ under
distribution $\text{P}(Z)$. This quantity is the mean of the marginal log
standard deviations of each of the $Q$ components of
$\theta\in\Theta\subseteq\mathbb{R}^{Q}$, which we negate so as to promote
marginal variability when performing minimisation. We work with the standard
deviation (instead of the variance) and take the logarithm thereof to minimise
the contribution of any particularly extreme marginal (i.e. marginal with
relatively low or high variance). Equations (3) and (4) make explicit the
dependence on $\text{P}(\theta\mid\lambda,\boldsymbol{X})$, and thus
$\lambda$, for clarity. We often have analytic expressions for
$\text{SD}_{\text{P}(\theta\mid\lambda,\boldsymbol{X})}[\theta_{q}]$, but
precise estimates are also simple to obtain using Monte Carlo. Note that
Equation (4) assumes
$\text{SD}_{\,\text{P}(\theta_{q}\mid\lambda,\boldsymbol{X})}[\theta_{q}]$
exists, and is nonzero and finite for all $q$ and $\lambda\in\Lambda$. Should
this not be true, for example if one of the marginals of
$\text{P}(\theta_{q}\mid\lambda,\boldsymbol{X})$ is a Cauchy distribution, we
can instead employ alternative, robust estimators of scale (Kravchuk and
Pollett, 2012; Rousseeuw and Croux, 1993).

###  2.5 Post optimisation decision step

We jointly minimise Equations (1) and (4) using a multi-objective optimisation
algorithm, which we will cover in Section 2.6. By adopting a multiple
objective approach to the translation problem, we obtain a set of possible
$\lambda$ values which comprise the Pareto frontier
$\mathcal{P}=\\{\lambda_{l}\\}_{l=1}^{\lvert\mathcal{P}\rvert}$, the set of
all “non-dominated” choices for $\lambda$, meaning that no point in
$\mathcal{P}$ is preferable in _both_ objectives to any of the remaining
points in $\mathcal{P}$ (for an introduction to multi-objective optimisation
problems see Chapter 2 of Deb, 2001). For each $\lambda$ in $\mathcal{P}$ we
compute the loss

$\tilde{L}(\lambda)=\log(\tilde{D}(\lambda\mid\boldsymbol{X}))+\kappa\,\tilde{N}(\lambda\mid\boldsymbol{X}),$
(5)

where the value of $\kappa>0$ expresses our relative belief in the importance
of the secondary objective. We take the log of
$\tilde{D}(\lambda\mid\boldsymbol{X})$ to be on a similar scale to our
definition of $\tilde{N}(\lambda\mid\boldsymbol{X})$ in Equation (4), but
stress that this not necessary should there be a more appropriate scale on
which to define $\tilde{L}(\lambda)$. The optimal value is then chosen such
that $\lambda^{*}:=\min\limits_{\lambda\in\mathcal{P}}\tilde{L}(\lambda)$.

This optimum is clearly sensitive to the choice of $\kappa$. If the relative
belief between the two objectives were known, $\kappa$ could be set directly,
but this will not usually be straightforward to assess. However, the multi-
objective optimisation approach makes it computationally inexpensive to test
many values of $\kappa$ (the set of which we denote with $\mathcal{K}$), and
thus plot the Pareto frontier coloured by loss, with the minimum loss point
indicated. These can guide our choice of $\kappa$: we can seek a value of
$\kappa$ with the minimum loss point not on the extreme of the Pareto
frontier, since we would like to balance the two objectives. This approach is
particularly useful in settings where the scales of the two optima differ
markedly, which we further discuss in Appendix C. Where it is feasible to
replicate the optimisation procedure, we can additionally seek a choice of
$\kappa$ that leads to Pareto frontiers with minimal variability across
replicates, since this suggests the optimal solution can be estimated
reliably.

###  2.6 Optimisation strategy

With the mathematical framework defined, we now turn to discuss the many
practicalities of optimisation. We use a two-stage global optimisation process
to construct a prior with the desiderata listed in Section 2.1. The first
stage focuses entirely on faithfulness by minimising only
$\tilde{D}(\lambda\mid\boldsymbol{X})$, using the variant of controlled random
search 2 (CRS2, Price, 1983) proposed by Kaelo and Ali (2006). Stage one
output is then used to initialise stage two, which additionally focuses on
uniqueness and replicability by employing multi-objective Bayesian
optimisation (Frazier, 2018; Zaefferer _et al._ , 2012) to jointly minimise
$\tilde{D}(\lambda\mid\boldsymbol{X})$ and
$\tilde{N}(\lambda\mid\boldsymbol{X})$. We focus on faithfulness, and thus
$D(\lambda\mid\boldsymbol{X})$, as a separate stage because minimising
$D(\lambda\mid\boldsymbol{X})$ is considerably more challenging than
minimising $N(\lambda\mid\boldsymbol{X})$. By initialising stage two with
faithful estimates for $\lambda$ we can, in the second stage, spend more
computational resources on finding points along the Pareto frontier. The
resulting optimal prior should suitably encode the information in the target
predictive distribution, without being overly confident for model parameters.
An idealised form of this process is illustrated in Figure 1.

Note that almost all global optimisation methods require, in the absence of
other constraints, $\Lambda$ to be compact subset of $\mathbb{R}^{L}$. We
require that the practitioner specify the upper/lower limits for each
dimension of $\Lambda$.

Figure 1: An idealised depiction of the methodology we introduce in this paper
when the corresponding densities are available. The starting point, or "stage
zero" (S0) depicts the elicited target prior predictive distribution for the
observable quantity $Y$ depicted and denoted by its density $\text{t}(Y)$ (red
line). Stage one starts (S1 – start) from an initial prior (blue line), which
does not match the target and is uninformative for the model parameter
$\theta$. Optimisation proceeds by using controlled random search to minimise
the predictive discrepancy (grey shaded area) and produces optimal
hyperparameters $\lambda^{*}$. Stage one produces (S1 – end) a very faithful
prior predictive distribution, but is overly confident for the model
parameter. Stage two (S2) uses multi-objective Bayesian optimisation and
corrects this overconfidence for $\theta$ with only a small increase in
predictive discrepancy.

####  2.6.1 Evaluating the objectives

As noted previously, evaluating $\tilde{N}(\lambda\mid\boldsymbol{X})$ for
models where analytic results or simple Monte Carlo estimates are available is
straightforward, and we denote the corresponding estimate (or, if available,
the analytic form) of $\tilde{N}(\lambda\mid\boldsymbol{X})$ with
$N(\lambda\mid\boldsymbol{X})$. However, there are two immediate challenges to
evaluating $\tilde{D}(\lambda\mid\boldsymbol{X})$:

  1. 1.

the prior predictive CDF $\text{P}(Y\mid\lambda,\boldsymbol{X})$ is often
analytically unavailable;

  2. 2.

the integral in Equation (1) is almost always intractable.

We address the former with a Monte Carlo based empirical CDF (ECDF), and the
latter with importance sampling. Specifically, given a particular value of
$\lambda$ and $X_{r}$, we draw $S_{r}$ samples
$\boldsymbol{y}_{r}^{(\text{P})}=(y_{s,r})_{s=1}^{S_{r}}$ with
$\boldsymbol{y}_{r}^{(\text{P})}\sim\text{P}(Y\mid\lambda,X_{r})$ to form the
ECDF $\hat{\text{P}}(Y\mid\lambda,X_{r},\boldsymbol{y}_{r}^{(\text{P})})$. To
apply importance sampling we rewrite the integral in Equation (1) with respect
to importance distribution $\text{Q}(Y\mid X_{r})$ and Radon-Nikodym
derivative $\text{d}\text{T}(Y\mid X_{r})\mathop{/}\text{d}\text{Q}(Y\mid
X_{r})$, such that

$\tilde{D}(\lambda\mid\boldsymbol{X})=\frac{1}{R}\sum_{r=1}^{R}\int
d(\text{P}(Y\mid\lambda,X_{r}),\text{T}(Y\mid
X_{r}))\frac{\text{d}\text{T}(Y\mid X_{r})}{\text{d}\text{Q}(Y\mid
X_{r})}\text{d}\text{Q}(Y\mid X_{r}).$ (6)

When $Y$ is discrete or of mixed type, $\text{q}(Y\mid X_{r})$ is instead a
probability mass function or an appropriate mixture of discrete and continuous
densities. Supposing we draw $I_{r}$ importance samples
$(y_{i,r})_{i=1}^{I_{r}}\sim\text{Q}(Y\mid X_{r})$, we denote the importance
sampling approximation to $\tilde{D}(\lambda\mid\boldsymbol{X})$ with
$D(\lambda\mid\boldsymbol{X})$ such that

$D(\lambda\mid\boldsymbol{X})=\frac{1}{R}\sum_{r=1}^{R}\frac{1}{I_{r}}\sum_{i=1}^{I_{r}}d(\text{P}(y_{i,r}\mid\lambda,X_{r}),\text{T}(y_{i,r}\mid
X_{r}))\frac{\text{d}\text{T}(y_{i,r}\mid X_{r})}{\text{d}\text{Q}(y_{i,r}\mid
X_{r})}.$ (7)

Note that we write $\text{Q}(Y\mid X_{r})$, to make clear that the importance
distribution could be covariate-specific, but in straightforward settings a
common importance distribution for all $R$ covariate values will be
appropriate.

We select $\text{Q}(Y\mid X_{r})$ using information about the support
$\mathcal{Y}$, and samples from $\text{P}(Y\mid\lambda,X_{r})$ and
$\text{T}(Y\mid X_{r})$. For more details see Appendix A. Finally, for
numerical stability we evaluate $D(\lambda\mid\boldsymbol{X})$ on the log
scale, with details available in Appendix B. This process is summarised in
Algorithm 1.

Algorithm 1 Importance sampling and ECDF approximation
$D(\lambda\mid\boldsymbol{X})$

1Targets $\text{T}(Y\mid X_{r})$ for $r=1,\ldots,R$; samplers for generating
points from $\text{T}(Y\mid X_{r})$ and $\text{P}(Y\mid\lambda,X_{r})$;
discrepancy $d(\cdot,\cdot)$; number of samples to draw $S_{r}$; number of
importance samples $I_{r}$; observable support $\mathcal{Y}$

2

3function Evaluate log $(D$($\lambda\mid\boldsymbol{X})$)

4 for $r$ in $1\ldots R$ do

5 Sample prior predictive
$\boldsymbol{y}^{(\text{P})}_{r}=(y_{s,r}^{(\text{P})})_{s=1}^{S_{r}}\sim\text{P}(Y\mid\lambda,X_{r})$

6 Use $\boldsymbol{y}^{(\text{P})}_{r}$ to form the ECDF
$\hat{\text{P}}(Y\mid\lambda,X_{r},\boldsymbol{y}^{(\text{P})}_{r})$

7 Sample target
$\boldsymbol{y}^{(\text{T})}_{r}=(y_{s,r}^{(\text{T})})_{s=1}^{S_{r}}\sim\text{T}(Y\mid
X_{r})$

8 Choose importance distribution $\text{Q}(Y\mid X_{r})$ via Appendix A

9 Sample importance points $(y_{i,r})_{i=1}^{I_{r}}\sim\text{Q}(Y\mid X_{r})$

10 end for

11 Compute $\log(D(\lambda\mid\boldsymbol{X}))$ using Equations (19) – (22) in
Appendix B

12 return: Value of $\log(D(\lambda\mid\boldsymbol{X}))$

13end function

####  2.6.2 Optimisation, stage 1

In this stage we focus solely on faithfulness by minimising
$D(\lambda\mid\boldsymbol{X})$. We do so using CRS2 (Price, 1983) with local
mutation (Kaelo and Ali, 2006), which we run for $N_{\text{CRS2}}$ iterations.
We make use of the final optimum value $\lambda^{*}$, as well as all of the
$N_{\text{CRS2}}$ trial points to obtain a design $\mathcal{D}$ for the next
stage. The design comprises values of $\lambda$, and their corresponding
values of $\log(D(\lambda\mid\boldsymbol{X}))$. A (small) number of padding
points $N_{\text{pad}}$ are added to $\mathcal{D}$ for numerical robustness in
stage 2. The result is the design
$\mathcal{D}=\left\\{\lambda_{i},\log(D(\lambda_{i}\mid\boldsymbol{X}))\right\\}_{i=1}^{N_{\text{design}}+N_{\text{pad}}}$,
whose construction is detailed in Algorithm 3 in Appendix D.

Whilst CRS2 was not designed to minimise noisy functions, it appears
empirically robust to small quantities of noise. We can make the noise in
$D(\lambda\mid\boldsymbol{X})$ arbitrarily small (by increasing $S_{r}$ and
$I_{r}$), but doing so usually incurs an enormous computational cost.
Carefully balancing the noise in the objective, and thus the quality of the
stage one solution, against the cost of evaluation yields a faithful optimum
$\lambda^{*}$ and useful design $\mathcal{D}$ in an acceptable amount of time.

####  2.6.3 Optimisation, stage 2

Stage two accounts for the secondary objective $N(\lambda\mid\boldsymbol{X})$
in addition to the primary objective. We adjust our optimisation technique to
affect this change in emphasis, and use multi-objective Bayesian optimisation,
via MSPOT (Zaefferer _et al._ , 2012), to jointly minimise
$D(\lambda\mid\boldsymbol{X})$ and $N(\lambda\mid\boldsymbol{X})$. MSPOT uses
a separate Gaussian process (GP) approximation to each of the objectives, and
evaluates these approximations at many points from Latin hypercube designs
(Stein, 1987). In each of the $N_{\text{BO}}$ iterations, the best points
under the current GP approximations are evaluated using the actual objectives.
These evaluations accumulate and thus iteratively improve the GP
approximations. After $N_{\text{BO}}$ iterations, the evaluated points are
reduced to their Pareto frontier (Kung _et al._ , 1975). Algorithm 4 in
Appendix D describes in detail the MSPOT algorithm as applied to two
objectives.

The noisy and computationally expensive nature of our objectives, particularly
$D(\lambda\mid\boldsymbol{X})$, necessitates an approach such as MSPOT.
Employing approximate GP models for the objectives allows us to screen
potential values of $\lambda\in\Lambda$ inexpensively, and avoid evaluating
the actual objectives at values of $\lambda$ far from optimal. Moreover, the
GP is a flexible yet data efficient model to use as an approximation and can,
through appropriate choice of kernel, capture correlation or other complex
relationships between components of $\lambda$ and the objective.

We use an optional batching technique in stage two because the computational
cost of evaluating the GP growing cubically in the number of points
$N_{\text{BO}}$ used in its construction. Batching avoids this problem by
subsampling to remove similar, or otherwise uninformative, points from the
collection used to form the surrogate GP model. Specifically, we run MSPOT for
a fixed number of iterations, then select a subsample of the points evaluated
in the completed batch as initialisation points for the following batch. In
our experience this approach can produce equivalent or better priors in less
time than a single batch of many iterations. The exact subsampling strategy is
detailed in Algorithm 5 in Appendix D.

###  2.7 Benchmarking and other empirical considerations

We will show results for both the multi-objective approach and the single-
objective approach, which does not regularise via the second objective (the
optimisation process is identical except stage 2 considers only
$D(\lambda\mid\boldsymbol{X})$).

Having selected $\kappa$ we compute the optimal $\lambda$ by minimising
$L(\lambda)$. Given $\lambda^{*}$, we will empirically assess the performance
of our method against the desiderata (Section 2.1). Faithfulness is
straightforward to assess empirically by comparing the target distribution
$\text{t}(Y\mid X_{r})$ and the estimated optimal prior predictive
distribution $\text{p}(Y\mid\lambda^{*},X_{r})$. Replicability and uniqueness
are more challenging to disentangle empirically: in the absense of
replicability, we are not able to conclude whether uniqueness holds or not. We
will first assess replicability by examining the stability of the components
of the loss in Equation 5 across independent replications of the optimisation
procedure. When the loss is stable across replicates, we will assess
uniqueness by examining whether the optimal prior predictive distribution
$\text{p}(Y\mid\lambda^{*},X_{r})$ and prior
$\text{P}(\theta\mid\lambda^{*},\boldsymbol{X})$ are stable across replicates;
the stability of both provides good evidence that of uniqueness.

###  2.8 Summary

Our method for specifying a prior given predictive information requires:

  1. 1.

a method for sampling $\text{P}(Y\mid\lambda,\boldsymbol{X})$;

  2. 2.

upper and lower limits that render $\Lambda$ a compact subset of
$\mathbb{R}^{L}$;

  3. 3.

a target CDF $\text{T}(Y\mid\boldsymbol{X})$ which, for numerical stability,
must be implemented on the log-scale;

  4. 4.

a method for generating samples from $\text{T}(Y\mid\boldsymbol{X})$, as the
importance sampler depends on information contained in such samples;

  5. 5.

a choice of $\kappa$ (and we will present a diagnostic plot–see e.g. Figure
10–which assists in making this decision).

Algorithm 2 describes, in pseudocode, our proposed methodology.

Algorithm 2 Methodology to translate prior predictive information into a prior
for the parameters in a complex model

1$\log(D(\lambda\mid\boldsymbol{X}))$ (evaluable using Algorithm 1); secondary
objective $N(\lambda\mid\boldsymbol{X})$; $\kappa$; number of Bayesian
optimisation iterations $N_{\text{BO}}$; number of batches $N_{\text{batch}}$;
number of CRS2 iterations $N_{\text{CRS2}}$; number of importance samples per-
covariate $I_{r}$; number of prior predictive samples per-covariate $S_{r}$.

2

3function pbbo($\kappa,N_{\text{BO}},N_{\text{batch}}$)

4 Minimising $\log(D(\lambda\mid\boldsymbol{X}))$ alone, compute the initial
design $\mathcal{D}$ using CRS2 via Algorithm 3

5 for $b$ in $1\ldots N_{\text{batch}}$ do

6 Jointly minimising $\log(D(\lambda\mid\boldsymbol{X}))$ and
$N(\lambda\mid\boldsymbol{X})$, compute the $b$th Pareto Frontier
$\mathcal{P}_{b}$ and complete design $\mathcal{D}_{b}$ using Algorithm 4,
initialising with design $\mathcal{D}$

7 Update design $\mathcal{D}$ using $\mathcal{P}_{b}$ and $\mathcal{D}_{b}$
via Algorithm 5

8 end for

9 With final Pareto frontier $\mathcal{P}_{N_{\text{batch}}}$, compute
$\lambda^{*}=\min
L(\lambda)=\min\limits_{\lambda\in\mathcal{P}_{N_{\text{batch}}}}\log(D(\lambda\mid\boldsymbol{X}))+\kappa
N(\lambda\mid\boldsymbol{X})$

10 return: $\lambda^{*}$

11end function

###  2.9 The pbbo R package

We implement our methodology in an R package (R Core Team, 2022) called pbbo,
available from https://github.com/hhau/pbbo. pbbo444The release associated
with this paper is available at https://doi.org/10.5281/zenodo.6817350. builds
on top of mlrMBO (Bischl _et al._ , 2018) for multi-objective Bayesian
optimisation, nlopt and nloptr (Ypma _et al._ , 2022; Johnson, 2014) for
global optimisation using CRS2 (Kaelo and Ali, 2006), and other packages for
internal functionality and logging (Wickham _et al._ , 2019; Rowe, 2016;
Maechler _et al._ , 2021). The code implementing the examples we consider in
Sections 3 to 5, which further illustrate pbbo, can be found at
https://gitlab.com/andrew-manderson/pbbo-paper.

##  3 Calibrating a cure fraction survival model

Cure models (Peng and Taylor, 2014; Amico and Van Keilegom, 2018) for survival
data are useful when a cure mechanism is physically plausible _a priori_ , and
when individuals are followed up for long enough to be certain all censored
individuals in our data are “cured”. Such lengthy follow ups are not always
possible, but a cure model remains plausible when a large fraction of the
censored observations occur after the last observed event time. However, we
cannot distinguish in the right tail of the survival time distribution between
censored uncured individuals and genuinely cured individuals.

We suppose here that we possess prior knowledge on the fraction of individuals
likely to be cured, and the distribution of event times amongst the uncured,
and seek to translate this information into a reasonable prior for the
parameters in a cure model. Consider individuals $n=1,\ldots,N$ with event
times $Y_{n}$ and censoring times $C_{n}$, such that $Y_{n}\in(0,C_{n}]$. We
suppose censoring times are distributed such that $C_{n}\sim
20+\text{Exp}(1)$. We simulate data for $N=50$ individuals, each with $B=4$
correlated covariates. We sample a single correlation matrix
$\boldsymbol{Q}\sim\text{LKJ}(5)$ (Lewandowski _et al._ , 2009) and
subsequently covariates
$\tilde{\mathbf{x}}_{n}\sim\text{MultiNormal}(\boldsymbol{0},\boldsymbol{Q})$.
This results in marginally-standardised yet correlated covariates.

There are several properties of this model that make it an interesting subject
for prior specification methodology. The observed quantity, and thus the
target distribution, is of mixed discrete/continuous type due to censoring.
Additionally, we specify a model with a nontrivial correlation structure,
about which we wish to specify an informative prior. Eliciting informative
priors for correlation structures is known to be challenging. Finally,
identifiability is known to be challenging in cure models (Peng and Taylor,
2014), and so the model is a demanding test of our regularisation procedure.

###  3.1 Target survival time distribution and covariate generation

Suppose that individuals are followed up for an average (but arbitrary) of 21
units of time, with those who experience the event doing so a long time before
the end of follow up. Furthermore, suppose we believe that, _a priori_ , $5\%$
of the patients will be cured, with $0.2\%$ of events censored due to
insufficient follow up.

A target distribution that is consistent with our beliefs comprises a point
mass of $0.05$ at $C_{n}$, and a lognormal distribution with location
$\mu^{\text{LN}}=\log(3)$ and scale $\sigma^{\text{LN}}=2\mathop{/}3$ for
$Y_{n}<C_{n}$. This choice of lognormal has $99.8\%$ of its mass residing
below 21, and thus produces event times that are “well separated” from the
censoring time, as required by a cure fraction model. Denoting the lognormal
CDF with $\text{LogNormal}(Y;\mu,\sigma^{2})$, we define the target CDF

$\text{T}(Y_{n}\mid
C_{n})=0.95\frac{\text{F}^{\text{LN}}(Y_{n};\mu^{\text{LN}},(\sigma^{\text{LN}})^{2})}{Z_{n}}+0.05\mathbbm{1}_{\left\\{Y_{n}=C_{n}\right\\}},\qquad
Y_{n}\in(0,C_{n}],$ (8)

where
$Z_{n}=\text{LogNormal}(C_{n};\mu^{\text{LN}},\left(\sigma^{\text{LN}}\right)^{2})$
is the required normalising constant. This individual-specific construction of
$\text{T}(Y_{n}\mid C_{n})$ implies that the number of covariates $R=N$, since
the censoring time $C_{n}$ functions as a covariate.

We simulate data for this example with $N=50$ individuals, each with $B=4$
correlated covariates. In line with our target distribution, simulated
censoring times are distributed such that $C_{n}\sim 20+\text{Exp}(1)$. We
sample a single correlation matrix $\boldsymbol{Q}\sim\text{LKJ}(5)$
(Lewandowski _et al._ , 2009) and subsequently covariates
$\tilde{\mathbf{x}}_{n}\sim\text{MultiNormal}(\boldsymbol{0},\boldsymbol{Q})$.
This results in marginally-standardised yet correlated covariates.

###  3.2 Model

A cure model for survival data, expressed in terms of its survival function,
is

$S(Y\mid
X,\theta)=\pi+(1-\pi)\tilde{S}(Y\mid\tilde{\boldsymbol{X}},\tilde{\theta}),$
(9)

where a proportion $\pi\in(0,1)$ of the population are _cured_ and never
experience the event of interest. The survival times for the remaining $1-\pi$
proportion of the population are distributed according to the _uncured_
survival function $\tilde{S}(Y\mid\tilde{\boldsymbol{X}},\tilde{\theta})$. We
use the tilde in $\tilde{\boldsymbol{X}}$ and $\tilde{\theta}$ to denote
quantities specific to the uncured survival distribution, and denote
$\theta=(\pi,\tilde{\theta})$ to align with our general notation.

Right censoring results in $Y_{n}=C_{n}$. The censoring indicator
$\delta_{n}=\mathbbm{1}_{\left\\{Y_{n}<C_{n}\right\\}}$ is 0 for right
censored events, and is 1 otherwise. We denote with $\tilde{\mathbf{x}}_{n}$
the $n$th row of the $N\times B$ covariate matrix $\tilde{\boldsymbol{X}}$.
Our model supposes that the uncured event times are distributed according to a
Weibull regression model, with survival function
$\tilde{S}(Y_{n}\mid\tilde{\theta},\tilde{\mathbf{x}}_{n},C_{n})$ and hazard
$\tilde{h}(Y_{n}\mid\tilde{\theta},\tilde{\mathbf{x}}_{n},C_{n})$ such that

$\begin{gathered}\tilde{S}(Y_{n}\mid\tilde{\theta},\tilde{\mathbf{x}}_{n},C_{n})=\exp\left\\{-Y_{n}^{\gamma}\exp\left\\{\beta_{0}+\tilde{\mathbf{x}}_{n}\boldsymbol{\beta}\right\\}\right\\},\quad
Y_{n}\in(0,C_{n}]\\\
\tilde{h}(Y_{n}\mid\tilde{\theta},\tilde{\mathbf{x}}_{n},C_{n})=\gamma
Y_{n}^{\gamma-1}\exp\left\\{\beta_{0}+\tilde{\mathbf{x}}_{n}\boldsymbol{\beta}\right\\},\end{gathered}$
(10)

with $\tilde{\theta}=(\gamma,\beta_{0},\boldsymbol{\beta})$. The likelihood
for the $n$th individual is

$\displaystyle\text{p}(Y_{n}\mid\theta,\tilde{\mathbf{x}}_{n},C_{n})=$
$\displaystyle\left((1-\pi)\tilde{S}(Y_{n}\mid\tilde{\theta},\tilde{\mathbf{x}}_{n},C_{n})\tilde{h}(Y_{n}\mid\tilde{\theta},\tilde{\mathbf{x}}_{n},C_{n})\right)^{\delta_{n}}$
(11)
$\displaystyle\times\left(\pi+(1-\pi)\tilde{S}(Y_{n}\mid\tilde{\theta},\tilde{\mathbf{x}}_{n},C_{n})\right)^{1-\delta_{n}}.$

To align the notation in this example with that introduced in Section 2 we
denote $Y=(Y_{n})_{n=1}^{N}$ and $X=(C_{n},\tilde{\mathbf{x}}_{n})_{n=1}^{N}$.
Note that we are using $X$ to represent _all_ information that must be
conditioned on, including the censoring times. This is necessary because, in
our more general notation, the support of $Y\mid X_{r}$ is truncated to an
interval that depends on $X_{r}$, making the censoring times necessary to
fully-specify $Y\mid X_{r}$.

We will seek to identify optimal values of the hyperparamers
$\lambda=(\alpha,\beta,\mu_{0},\sigma^{2}_{0},s_{\beta},\boldsymbol{\omega},\boldsymbol{\eta},a_{\pi},b_{\pi})^{\top}$,
where:

$\begin{gathered}\pi\sim\text{Beta}(a_{\pi},b_{\pi}),\quad\gamma\sim\text{Gamma}(\alpha,\beta),\\\
\beta_{0}\sim\text{Normal}(\mu_{0},\sigma_{0}^{2}),\quad\boldsymbol{\beta}\sim\text{MVSkewNormal}(\boldsymbol{0},\boldsymbol{S},\boldsymbol{\eta}).\end{gathered}$
(12)

The upper and lower limits that define $\Lambda$ are specified in Table 2 in
Appendix E. Our use of the multivariate skew-normal prior has two motivations.
The skewness is necessary to incorporate the nonlinear relationship between
the hazard and the effect of the covariates, and a covariance structure is
used to account for fact that not all the elements of $\boldsymbol{\beta}$ can
be large simultaneously. We assume that $\tilde{\boldsymbol{X}}$ is marginally
(column-wise) standardised, and can thus decompose
$\boldsymbol{S}=\text{diag}(s_{\beta})\,\,\boldsymbol{\Omega}\,\,\text{diag}(s_{\beta})$
where $s_{\beta}$ is the scale of the prior marginals of $\boldsymbol{\beta}$.
The standardisation allows us to use only one $s_{\beta}$ instead of one per
covariate. We elect to parameterise $\boldsymbol{\Omega}$ using the $(B-1)!=6$
elements that uniquely determine its Cholesky factor, which we denote
$\boldsymbol{\omega}=(\omega_{1},\ldots,\omega_{6})^{\top}\in[-1,1]^{6}$.
These elements are transformed into $\boldsymbol{\Omega}$ using the partial
correlation method of Lewandowski _et al._ (2009), also employed by the Stan
math library (Stan Development Team, 2022). The $B$-vector $\boldsymbol{\eta}$
controls, but is not equal to, the marginal skewness for each element of
$\boldsymbol{\beta}$ using the multivariate skew-normal definition of Azzalini
and Valle (1996), as implemented in the sn package (Azzalini, 2022).

###  3.3 Tuning parameters and further details

In this example we again employ the multi-objective approach with $N(\lambda)$
as defined in Equation 4. We also run single-objective optimisation to
facilitate assessing uniqueness. Total predictive discrepancy is computed
using both the Anderson-Darling and Cramér-von Mises discrepancy functions.

We use $N_{\text{CRS2}}=2000$ CRS2 iterations, followed by
$N_{\text{batch}}=3$ batches of multi-objective Bayesian optimisation using
$N_{\text{BO}}=200$ iterations per batch, carrying forward
$N_{\text{design}}=60$ points between batches. The predictive discrepancy
function, $D(\lambda)$, is estimated using $S_{r}=10^{4}$ samples from the
prior predictive, and evaluated using $I_{r}=5\times 10^{3}$ importance
samples from an appropriate mixed discrete/continuous importance density.

###  3.4 Results

####  3.4.1 Choosing $\kappa$ for the multi-objective approach

We inspect the Pareto frontiers and compute the minimum loss points for 6
values of $\kappa\in\\{0.1,0.2,0.3,0.5,1,2\\}$. These are displayed, for the
Anderson-Darling discrepancy, in Figure 2, and Appendix E.2 Figure 17 for the
Cramér-von Mises discrepancy. To reiterate, the assumption is that nearby
points on the Pareto frontier are “similar” priors for
$\text{p}(\theta\mid\lambda)$.

The maximum and minimum values for $\kappa$ yield a number of minimum loss
points on the extremes of the Pareto frontier. For the maximum, $\kappa=2$,
these points correspond to unsuitable values for $\lambda$ due to a lack of
faithfulness, which we will demonstrate momentarily. The remaining options for
$\kappa$ result in similar minimum loss points, illustrating that the minimum
loss point is constant, or very similar, across a range of $\kappa$ values.
Such insensitivity is ideal as we can obtain sensible solutions for a wide
variety of $\kappa$ values. We select $\kappa=0.3$ as this simultaneously
minimises the variability in loss and both objective functions.

Figure 2: Pareto frontiers for the survival example using the Anderson-Darling
discrepancy, for the values of $\kappa$ we consider. Note that the range
associated with the colour scale differ between panels. The red crosses
($\color[rgb]{0.70703125,0.05859375,0.125}\definecolor[named]{pgfstrokecolor}{rgb}{0.70703125,0.05859375,0.125}{+}$)
indicate the minimum loss point on each frontier, for each value of $\kappa$.

####  3.4.2 Replicability of loss at the optima

Figure 3 displays $\log(D(\lambda\mid\boldsymbol{X}))$, or $\log(D(\lambda))$
in the single objective setting, $N(\lambda\mid\boldsymbol{X})$, and
$L(\lambda\mid\boldsymbol{X})$ at the replicate optima $\lambda^{*}$, using
$\kappa=0.3$ in the multi-objective setting. We observe a tight distribution
of objective values, particularly for $\log(D(\lambda^{*}\mid\boldsymbol{X}))$
and $\log(D(\lambda^{*}))$, suggesting our results are replicable in this
example.

Figure 3: Value of the objectives $D(\lambda^{*}\mid\boldsymbol{X})$,
$N(\lambda^{*}\mid\boldsymbol{X})$, and the loss
$L(\lambda\mid\boldsymbol{X})$ at the replicated optima $\lambda^{*}$ for
$\kappa=0.3$. The single objective optima and predictive discrepancy
$\log(D(\lambda^{*}))$ are displayed in the right panel. The Anderson-Darling
discrepancy function is displayed in red, with the Cramér-von Mises in blue.
The points have been jittered along the x-axis for readability.

####  3.4.3 Faithfulness

We assess faithfulness by inspecting Figure 4, which displays the prior
predictive distributions at the optima and target for the randomly selected
individual $n=9$ in our simulated population (other individuals are not
visually distinguishable). In addition to the possible values of $\kappa$, the
optimal prior predictive distributions are also displayed using the single
objective approach for reference. Across all replicates and different
$\kappa$, the estimated optimal prior predictive distribution is highly
faithful to the target. The multi-objective results are close to those
obtained using the single objective approach, except for the two largest
values of $\kappa$. By overly valuing the secondary objective, the negative
mean log standard deviation, a fraction of the maximum $\kappa$ fits are
considerably worse. The poor fits correspond to the minimum loss points in the
bottom-right corner of the $\kappa=2$ panel in Figure 2.

Figure 4: Estimated optimal prior predictive densities
$\text{p}(Y_{n}\mid\lambda^{*})$ (red/blue lines and dots) and target
densities $\text{t}(Y_{n}\mid C_{N})$ (black lines and crosses) for randomly
selected individual $n=9$. The continuous portions of the densities are
displayed as lines, with the discrete/censored portion displayed as a point at
$Y_{n}=C_{n}$. The top rows (red) correspond to the multi-objective approach
using the values of $\kappa$ in the row titles, with the bottom row (blue)
displaying the single objective results. The left and right columns use the
Anderson-Darling and Cramér-Von Mises discrepancies respectively. Densities
are truncated to $[0,0.75]$ for readability.

####  3.4.4 Uniqueness

Figure 5: Estimated optimal prior marginal densities of
$\text{p}(\theta\mid\lambda^{*})$ for each component of $\theta$ in the
survival example. The top four rows (red) are obtained using the multi-
objective approach, under both discrepancy functions (AD – Anderson-Darling,
CvM – Cramér-von Mises). The selected value of $\kappa=0.3$ (top two rows), as
well as the extreme value $\kappa=2$ (middle two rows) are displayed. The
bottom two rows (blue) use the single, predictive discrepancy objective and
both discrepancy functions. Both axes are manually truncated for readability.
As such, many of the replicates for $\theta=\gamma$ and $\kappa=2$ are not
visible because the samples used to estimate $\text{p}(\gamma\mid\lambda^{*})$
lie between $10^{3}$ and $5\times 10^{5}$.

We now evaluate the uniqueness of the optimisation problem. Figure 5 displays
the marginals of $\theta$ under both the multi-objective ($\kappa=0.3,2$) and
the single objective approach, for each independent replicate.

The single objective approach consistently locates the degenerate, non-unique
solution where all the variation in the uncensored event times is attributed
to the baseline hazard shape $\gamma$ and the intercept $\beta_{0}$: i.e. all
the mass for $\boldsymbol{\beta}$ (the regression coefficients) is close to 0.
The combination of $\gamma$ and $\beta_{0}$ is evidently poorly identified,
and further calculation reveals that only the derived product
$\gamma\exp{\beta_{0}}$ is uniquely determined. Given the inter-replicate
consistency previously observed in Figure 3 we infer that the optimisation
process is locating equally good, in terms of predictive discrepancy, optima
that correspond to different $\text{P}(\theta\mid\lambda^{*})$, which is
precisely our definition of nonuniquness. Furthermore, the concentration
around $\boldsymbol{\beta}=0$ means these priors are unlikely to be efficient
for inferring the (possible) relationship between covariates and outcome – a
very strong signal in the data would be needed to overcome this prior.

Our multi-objective approach produces a more reasonable and disperse prior,
but still does not admit a unique optimal prior. The typical marginal prior
for $\boldsymbol{\beta}$ has been widened, and for $\kappa=0.3$ there is a
preference for the optima surrounding $\gamma\approx 7.5,\beta_{0}\approx-10$;
an improvement in uniqueness over the single objective approach, but
imperfect. When $\kappa=2$ the process regularly produces a marginal prior for
the cure fraction $\pi$ that gives considerable mass to values of $\pi>0.3$,
which is incongruent with our target. This is also visible in the censored
portions of the prior predictive estimates
$\text{p}(Y_{n}\mid\lambda^{*},C_{n})$ displayed in Figure 4.

Figure 6: Contours of the log prior density
$\log(\text{p}(\beta_{3},\beta_{4}\mid\lambda^{*}))$ at the optima. The left
column corresponds to the single objective approach, with the multi-objective
approach using $\kappa=0.3$ displayed in the right column. The top and bottom
row use the Cramér-von Mises and Anderson-Darling discrepancies respectively.
Note that, for clarity, we only plot the final 12 of 30 replicates, with
unique colouring for each replicate.

Eliciting covariance structures is challenging, and in this example we
included a covariance matrix in the hyperparameters for our model. In Figure 6
we display the bivariate prior marginal densities for two elements of
$\boldsymbol{\beta}$, specifically $\beta_{3}$ and $\beta_{4}$, for both the
multi-objective approach with $\kappa=0.3$ and the single objective approach.
Nonuniqueness is visible in both sets of estimates. There are marginal
densities of both positive and negative marginal skewness, and pairwise
correlation. The wider typical marginal for $(\beta_{3},\beta_{4})$ obtained
using the multi-objective approach is again visible in the right panel of
Figure 6. Lastly, there is a visible difference between the discrepancy
functions in the single objective approach, with the Cramer-von Mises
discrepancy typically producing a slightly wider marginal.

###  3.5 Example summary

Our procedure estimates priors that faithfully represent provided information
about the survival distribution. Replicability is reasonable in this case,
particularly for the loss, although the noise in Figure 2 would, ideally, be
smaller. Several priors appear to be consistent with the target distribution
we specify, and so the prior our procedure selects is not unique. The multi-
objective approach induces more concentration than the single objective
approach in the choice of prior for $\gamma$, but it is still not unique. A
particular challenge for uniqueness is the covariance structure for the vector
of covariate coefficients.

In the information we supply via $\text{T}(Y_{n}\mid C_{n})$, we are
completely certain of the fraction of cured patients _a priori_. Such
certainty is unlikely to be uncovered when eliciting information from experts.
A more elaborate construction of $\text{T}(Y_{n}\mid C_{n})$, or elaborate
methodology, may be able to represent such uncertainty, but the example
remains challenging without this additional complication.

##  4 Priors from model-derived quantities

Consider the linear model $Y=\boldsymbol{X}\boldsymbol{\beta}+\varepsilon$ for
$n\times p$ design matrix $\boldsymbol{X}$ and $p$-vector of coefficients
$\boldsymbol{\beta}$ indexed by $j=1,\ldots,p$, and where the noise
$\varepsilon$ has zero mean and variance $\sigma^{2}$. Suppose information
about the fraction of variance explained by the model is available – from
previous similar experiments, or from knowledge of the measurement process –
in the form of a plausible distribution for the coefficient of determination,
$R^{2}$, which for this model can be computed as

$R^{2}=1-\frac{\sigma^{2}}{n^{-1}\boldsymbol{\beta}^{\top}\boldsymbol{X}^{\top}\boldsymbol{X}\boldsymbol{\beta}+\sigma^{2}}.$
(13)

assuming that the columns of $\boldsymbol{X}$ have been centred. Our aim is to
use our knowledge of $R^{2}$ to set suitable priors for the regression
coefficients $\beta$. This idea was the inspiration for a class of shrinkage
priors (Zhang and Bondell, 2018; Zhang _et al._ , 2022), but we would like to
make this idea applicable to a wider selection of prior structures.

To illustrate the effect of including knowledge of $R^{2}$ on increasingly
complicated priors, we investigate the selection of appropriate
hyperparameters for three priors for the regression coefficients: two
shrinkage priors, and a simple Gaussian prior. We simultaneously vary the
covariate-independent target distribution $\text{T}(R^{2})$ to assess:

  * •

each prior’s ability to faithfully encode the information present across a
wide variety of target distributions;

  * •

uniqueness of the optimisation problem, and replicability of the single-
objective variant of our optimisation algorithm, for each prior/target pair.

Note that we assume throughout that the noise $\varepsilon$ is distributed
according to a Gaussian distribution with zero mean and variance $\sigma^{2}$,
with an $\text{InverseGamma}(a_{1},b_{1})$ prior on $\sigma^{2}$. To
demonstrate the challenge that this poses for uniqueness, we will seek to
select suitable $a_{1}$ and $b_{1}$ using our methodology.

##### Gaussian prior

The Gaussian prior has only one hyperparameter $\gamma$, which controls the
ratio of prior variability due to $\boldsymbol{\beta}$ to that of
$\varepsilon$, and is

$\beta_{j}\sim\text{N}\left(0,\frac{\sigma^{2}}{\gamma}\right).$ (14)

Hence, we denote parameters
$\boldsymbol{\theta}_{\text{GA}}=(\boldsymbol{\beta},\sigma^{2})$, and seek
optimum values for hyperparameters
$\boldsymbol{\lambda}_{\text{GA}}=(\gamma,a_{1},b_{1})$.

##### Dirichlet-Laplace prior (Dir. Lap.)

Bhattacharya _et al._ (2015) introduce the Dirichlet-Laplace shrinkage prior,
which is defined for the $j$th coefficient such that

$\beta_{j}\sim\text{Laplace}\left(0,\sigma\phi_{j}\tau\right),\quad(\phi_{1},\ldots,\phi_{p})\sim\text{Dirichlet}(\alpha,\ldots,\alpha),\quad\tau\sim\text{Gamma}(p\alpha,1\mathop{/}2).$
(15)

There is a single hyperparameter $\alpha$, with smaller values of $\alpha$
yielding more sparsity in $\boldsymbol{\beta}$. Thus we denote
$\boldsymbol{\lambda}_{\text{DL}}=(\alpha,a_{1},b_{1})$ and
$\boldsymbol{\theta}_{\text{DL}}=(\boldsymbol{\beta},\sigma^{2},\phi_{1},\ldots,\phi_{p},\tau)$.

##### Regularised horseshoe prior (Reg. Horse.)

The regularised horseshoe of Piironen and Vehtari (2017) is the most complex
of the priors. With more intermediary stochastic quantities between the
hyperparameters and $R^{2}$, as well as less linearity in the relationship
between the aforementioned, it is the most flexible of the priors. These
properties make finding optimal values of the hyperparameters more
challenging. The prior is

$\begin{gathered}c^{2}\sim\text{InvGamma}\left(\frac{\nu}{2},\frac{\nu
s^{2}}{2}\right),\quad\omega\sim\text{Cauchy}^{+}\left(0,\frac{p_{0}}{p-p_{0}}\sqrt{\frac{\sigma^{2}}{n}}\right),\quad\delta_{j}\sim\text{Cauchy}^{+}(0,1),\\\
\tilde{\delta}^{2}_{j}=\frac{c^{2}\delta^{2}_{j}}{c^{2}+\omega^{2}\delta^{2}_{j}},\quad\beta_{j}\sim\text{N}(0,\omega^{2}\tilde{\delta}^{2}_{j}),\end{gathered}$
(16)

where $\text{Cauchy}^{+}$ denotes the Cauchy distribution truncated to
$[0,\infty)$. Equation 16 leaves us free to choose three prior-specific
hyperparameters, $(p_{0},\nu,s^{2})$. Thus
$\boldsymbol{\lambda}_{\text{HS}}=(p_{0},\nu,s^{2},a_{1},b_{1})$ and
$\boldsymbol{\theta}_{\text{HS}}=(\boldsymbol{\beta},\sigma^{2},c^{2},\omega,\delta_{1},\ldots\delta_{p})$.
Whilst the regularised horseshoe is carefully designed to make
$(p_{0},\nu,s^{2})$ interpretable and easy to choose, here we aim to see if
values of these hyperparameters can be chosen to match an informative prior
for $R^{2}$.

###  4.1 Target predictive distributions

We consider sixteen different $\text{Beta}(s_{1},s_{2})$ distributions as our
target $\text{T}(R^{2})$, with
$\\{s_{1},s_{2}\\}\in\mathcal{S}\times\mathcal{S}$ and $\mathcal{S}$ chosen to
be four exponentially-spaced values between and including $1\mathop{/}3$ and
$3$ (i.e. equally-spaced between $\log(1\mathop{/}3)$ and $\log(3)$). These
values represent a variety of potential forms of the supplied target
predictive distribution for $R^{2}$.

###  4.2 Evaluation setup and tuning parameters

We fix $n=50$ and $p=80$ with entries in $\boldsymbol{X}$ drawn from a
standard Gaussian distribution, and assess replicability using 10 independent
runs for each prior and target. The support $\Lambda$ for the hyperparameters
is defined in Table 3 in Appendix F. We run pbbo with $S=10^{4}$ samples from
the prior predictive distribution, use both $d^{\text{AD}}$ and
$d^{\text{CvM}}$ as discrepancy functions, and evaluate the log predictive
discrepancy using $I=5\times 10^{3}$ samples from a $\text{Uniform}(0,1)$
importance distribution. We run the first stage of the optimiser for
$N_{\text{CRS2}}=1000$ iterations, and subsequently perform both single and
multi-objective Bayesian optimisation for $N_{\text{batch}}=1$ batch of
$N_{\text{BO}}=150$ iterations, using $N_{\text{design}}=50$ points from the
first stage. The single objective approach illustrates that differences in
flexibility between priors also induce differences in uniqueness, and
highlights issues in choosing a prior for the additive noise parameter
$\sigma^{2}$. Choosing $\kappa$ is challenging in the multi-objective
approach, as its value should depend on the target, the discrepancy function,
and the prior. These dependencies result in 96 possible choices of $\kappa$,
which is and infeasible number of choices to make in this example. Instead we
fix $\kappa=0.5$ for all multi-objective settings. Lastly, we compute the
secondary objective $N(\lambda)$ using the marginal standard deviation for all
elements in $\theta$ as in Equation (4), except for those quantities where,
for some $\lambda\in\Lambda$, the standard deviation is undefined. In such
cases we use the robust scale estimator from Rousseeuw and Croux (1993).

###  4.3 Results

We first assess the replicability of the discrepancy by plotting the value of
the objective at the optima $\lambda^{*}$ across replications. In Figure 7
each value of $\log(D(\lambda^{*}))$ is the mean of 10 evaluations of
$\log(D(\lambda))$ at each $\lambda^{*}$ to minimise residual noise in the
objective. Generally, it appears that the discrepancy is replicable for the
both Gaussian and Dirichlet-Laplace, for both Anderson-Darling and Cramér-von
Mises discrepancies. In contrast, the results for the Regularised Horseshoe
prior appear to replicate poorly under the Anderson-Darling discrepancy, but
reasonably under the Cramér-von Mises discrepancy.

Figure 7: Total log predictive discrepancy at the optima
$\log(D(\lambda^{*}))$. The target densities (denoted in the column titles)
correspond to the bottom row of Figure 8 (i.e. the same as Figure 9). Each
point corresponds to one of the 10 distinct replicates, and its value is the
mean of 10 evaluations of $\log(D(\lambda^{*}))$ for the same $\lambda^{*}$.
The top row displays Anderson-Darling discrepancy function, with the bottom
row displaying the Cramér-von Mises discrepancy function. Each prior is
displayed in the same colour as Figure 9, with the point shape displaying the
single or multi-objective approach.

We evaluate the faithfulness of the resulting prior distributions by
inspecting the densities $\text{p}(R^{2}\mid\lambda^{*})$ and
$\text{t}(R^{2})$ for the various targets (all distributions in this example
have corresponding densities). A selected subset of the pairs of
$(s_{1},s_{2})$ values are displayed in Figure 8 (complete results are in
Appendix F Figure 21). The faithfulness of the Gaussian prior is universally
poor, which we investigate further in Appendix F. Both shrinkage priors
perform better in cases where one of $s_{1}$ or $s_{2}$ is less than 1, with
the regularised horseshoe performing better for the $s_{1}=s_{2}>1$ cases.
Interestingly, the results are not symmetric in $s_{1}$ and $s_{2}$; the
Dirichlet-Laplace prior is able to match the $s_{1}=3,s_{2}=0.69$ target well,
with many of regularised horseshoe replicates performing poorly; whilst the
relative performance is reversed for $s_{1}=0.69,s_{2}=3$ (see Figure 21).
There is also perceptibly more variability in the regularised horseshoe
replicates, which suggests the optimisation problem is more challenging and
the predictive discrepancy objective is noisier. The multi-objective generally
produces more variable sets of optima, which is expected, as it is a more
difficult optimisation problem but we do not allow it additional computational
resources. There is little visible difference between the two discrepancy
functions. Finally, as the values of $s_{1}$ and $s_{2}$ increase, the
faithfulness of the shrinkage priors generally decreases. Across the full set
of simulations, the regularised horseshoe is evidently the most flexible
(Appendix F Figure 21).

Figure 8: Optimal prior predictive densities $\text{p}(R^{2}\mid\lambda^{*})$
for the three priors considered, for selected target densities. The column-
titles denote the target, which is also plotted as a black dashed line. The
row-titles denote which discrepancy is used to compute the total predictive
discrepancy, and whether the single or multi-objective approach to
optimisation is used. Each replicate of the Gaussian (’Gaussian’ – green),
Dirichlet-Laplace (’Dir. Lap.’ – red), and regularised horseshoe (’Reg.
Horse.’ – blue) is drawn in their respective colours. Density values are
trimmed to $[0,10]$ for readability.

To assess uniqueness, we consider estimated optimal hyperparameter values
$\lambda^{*}$ in each replicate. Figure 9 displays the estimates for $s_{1}=3$
and $s_{2}\in\\{0.33,0.69,1.44,3\\}$, which corresponds to the targets in
Figure 8. The estimates for $\gamma$ and $\alpha$, for the Gaussian and
Dirichlet-Laplace priors respectively, are consistent across replicates, which
suggests the optima may be unique. This remains true even for targets where
the prior is not faithful to the target, e.g. the $\text{Beta}(3,3)$ target.
There is more variability in the hyperparameters of the regularised horseshoe
prior. There does appear to be unique solution for $u$ for the
$\text{Beta}(3,0.33)$ and $\text{Beta}(3,0.69)$ targets, whereas $p_{0}$ and
$s^{2}$ are highly variable across replicates, which may reflect nonuniqueness
or may be due to the lack of replicability (discussed above) of this
optimisation for the regularised horseshoe.

The hyperparameters $(a_{1},b_{1})$ for the additive noise variance
$\sigma^{2}$ are highly variable across replications for almost all
prior/target combinations. This reflects the anticipated lack of uniqueness
when incorporating such hyperparameters. It is particularly striking for the
Dirichlet-Laplace prior when $s_{2}\in\\{0.33,0.69\\}$, where we observe
consistent and excellent fits/faithfulness to the target, but these do not
correspond to replicable estimates for $(a_{1},b_{1})$. These settings are
also interesting as the choice of single or multi-objective approach greatly
impacts the optimum values of $a_{1}$ and $b_{1}$. Faithfulness of the multi-
objective optima, illustrated in Figure 8, are not appreciably worse than the
single objective approach, but the inclusion of $\sigma^{2}$ into the
secondary objective has resulted in optimal values of $a_{1}$ and $b_{1}$ that
maximise the dispersion of the marginal prior (i.e. small $a_{1}$ and large
$b_{1}$). Asymptotic results are also known for the Gaussian prior, and in
Appendix F we further assess replicability by benchmarking our optimisation
process against suitable ‘true’ (asymptotically) values.

Figure 9: Optimal values $\lambda^{*}$ for each of the three priors
considered. Columns contain (possibly prior-specific) hyperparameters, with
the point colour corresponding to a specific prior. Each point’s shape
corresponds to the combination of discrepancy function and single or multi-
objective approach. The target beta densities (denoted by the row panel
titles) correspond to Figure 8.

###  4.4 Comparison of discrepancy functions

Our optimisation procedure has minimised $\log(D(\lambda))$ using both the
Anderson-Darling and Cramér-von Mises discrepancy functions. The former places
extra emphasis on matching the tails of the target, and thus the Regularised
Horseshoe values in the top row of Figure 7 differ from our expectations given
the results in the bottom row of Figure 8. Take, for example, the
$s_{1}=3,s_{2}=0.69$ case. It is plainly evident from Figure 8 that the
regularised horseshoe prior provides a better fit to the target distribution
at $\boldsymbol{\lambda}^{*}_{\text{HS}}$, and yet the corresponding
$\log(D(\boldsymbol{\lambda}_{\text{HS}}^{*}))$ values in the top row of
Figure 7 suggest that it is considerably worse that the Gaussian prior at
$\boldsymbol{\lambda}_{\text{GA}}^{*}$. To reconcile this apparent
contradiction, we inspect $\log(D(\lambda))$ at the optima computed using the
Cramér-Von Mises discrepancy function. These values are displayed in the
bottom row of Figure 7, whose values closely match our expectations given
Figure 8. Given the range of behaviours of $\text{p}(R^{2}\mid\lambda^{*})$
for all the optima, we can conclude that the Anderson-Darling more heavily
penalises over-estimation of the tails of $\text{p}(R^{2}\mid\lambda^{*})$
than under-estimation. This does not discount it as an optimisation objective,
but does complicate comparisons between competing priors.

###  4.5 Example summary

This example illustrates the use of a model-derived, nonobservable quantity
about which we have prior information as the basis for an informative prior.
The most flexible shrinkage model (the regularised horseshoe prior) was
clearly the most faithful to the supplied information in almost all cases.
Conversely, the Gaussian prior is the most replicable and unique, but the lack
of faithfulness means it is unsuitable to use as a prior when seeking to use a
Beta prior on $R^{2}$. The example also illustrates the difficulty of learning
about the prior for the additive noise term. Attempts to ameliorate this
difficulty by adopting the multi-objective approach merely hide the issue; for
the Gaussian and Dirichlet-Laplace priors, where replicability is evident, we
always select the inverse gamma prior that maximises the standard deviation of
$\sigma^{2}$ by being on the boundary of $\Lambda$. Such a prior for
$\sigma^{2}$ is unlikely to prove appropriate nor represent our prior
knowledge.

##  5 A human-aware prior for a human growth model

We now consider a nonlinear regression model for human growth. There are a
number of properties of this example that make it an interesting test for our
methodology. First, we find it difficult to specify priors congruent with
desired prior predictive distributions for such models; both the nonlinearity
and the obvious dependence of height on age complicate prior specification.
Data for human growth are also readily available, so we can assess the impact
of the prior on many distinct data and posteriors. Second, the model we
consider is also poorly behaved under the flat prior, so some prior
information is required to stabilise and/or regularise the estimate of the
posterior. Finally, this example is also considered by Hartmann _et al._
(2020), and so there is a suitable comparator for our results.

Suppose an individual has their height measured at age $t_{m}$ (in years) for
$m=1,\ldots,M$, with corresponding measurement $y_{m}$ (in centimetres). The
first Preece-Baines model (Preece and Baines, 1978) for human height is,

$\displaystyle y_{m}$ $\displaystyle=h(t_{m};\theta)+\varepsilon_{m}$ (17)
$\displaystyle=h_{1}-\frac{2(h_{1}-h_{0})}{\exp\\{s_{0}(t_{m}-\gamma)\\}+\exp\\{s_{1}(t_{m}-\gamma)\\}}+\varepsilon_{m},$
(18)

with $\varepsilon_{m}\sim\text{N}(0,\sigma^{2}_{y})$. Some constraints are
required to identify this model and ensure its physical plausibility:
specifically, we require $0<h_{0}<h_{1}$ and $0<s_{0}<s_{1}$. A
parameterisation that respects these constraints and is easier to work with
uses $\delta_{h}=h_{1}-h_{0}$ instead of $h_{1}$, and $\delta_{s}=s_{1}-s_{0}$
in place of $s_{1}$. All of $(h_{0},\delta_{h},s_{0},\delta_{s})$ thus have
the same positivity constraint. Finally, we also constrain $\gamma$ such that
$\gamma\in(\min_{m}(t_{m}),\max_{m}(t_{m}))$. However, these constraints are
not sufficient to make the model plausible for all permissible parameter
values – the denominator of the fraction can be very small, yielding negative
heights.

To align with the notation introduced in Section 2 we denote the parameters by
$\theta=(h_{0},\delta_{h},s_{0},\delta_{s},\gamma)$. As in Hartmann _et al._
(2020), we choose for each of the $q=1,\ldots,5$ elements of $\theta$ an
independent $\text{LogNormal}(\mu_{q},s^{2}_{q})$ prior. We will seek to
identify the optimal values of
$\lambda=\left(\mu_{q},s^{2}_{q}\right)_{q=1}^{5}$. Table 5 in Appendix 5
lists the upper and lower limits we choose for each component of $\lambda$.

We do not consider the measurement error variance $\sigma_{y}^{2}$ as part of
$\theta$. Doing so introduces a degenerate solution for $\lambda$, where all
variability in $Y$ is singularly attributable to $\varepsilon_{m}$ and thus
$\sigma^{2}_{y}$. Such a prior seems undesirable, so instead we fix the prior
for $\sigma_{y}^{2}$ to reflect the measurement process for human height;
measurement errors are unlikely to be more than one or two centimetres, so
values of $\sigma_{y}^{2}\approx 1$ seem reasonable. Thus we set
$\sigma_{y}\sim\text{LogNormal}(0,0.2^{2})$.

###  5.1 Target predictive distributions

Our data are assumed to originate from a sample of adolescent humans,
uniformly distributed between ages 2 and 18, and evenly split between sexes.
We consider two different target prior predictive distributions. We first
consider a _covariate-independent_ prior predictive density $\text{t}(Y)$ with
corresponding CDF $\text{T}(Y)$ for human heights across the entire age-range,
derived by summarising external data. This target (Figure 12) is a mixture of
3 gamma densities specified to approximate the external data, which is
multimodal due to the fact that humans grow in spurts. We also consider a
_covariate-specific_ $\text{T}(Y\mid X_{r})$ in which we specify the target
predictive distribution $\text{T}(Y\mid X_{r})$ of human heights at ages
$X_{r}\in(2,8,13,18)$ with $r=1,\ldots,4$. Each $\text{T}(Y\mid X_{r})$ is
normal (see Figure 13, and Appendix G.1 for details).

###  5.2 Tuning parameters and comparison with Hartmann et. al.

For both targets we obtain $\lambda^{*}$ using both the single objective and
multi-objective optimisation processes (see Sections 2.6 and 2.7). We use only
the Cramér-Von Mises discrepancy in this example, because we were not able to
reliably compute the Anderson-Darling discrepancy due to numerical
instabilities. We run pbbo using $S=5\times 10^{4}$ samples from
$\text{p}(Y\mid\lambda)$ and likewise $S_{r}=5\times 10^{4}$ samples from
$\text{p}(Y\mid\lambda,X_{r})$ for each of the 4 values of $X_{r}$. We use
$I=5\times 10^{3}$ and $I_{r}=5\times 10^{3}$ importance samples, with
$N_{\text{CRS2}}=2000$ CRS2 iterations, $N_{\text{batch}}=5$ Bayesian
optimisation batches each of $N_{\text{BO}}=250$ iterations, and carry forward
$N_{\text{design}}=50$ points per batch. We assess replicability using 30
independent runs of each objective/target pair.

For this example, we also assess the stability of the posterior: our prior
ideally admits a posterior amenable to sampling (i.e. removes spurious modes
and eliminates computational difficulties present when using a noninformative
prior). It also seems desirable that similar priors should yield similar
posteriors. We assess stability by exploit the robust sampling diagnostics in
Stan. Should any diagnostic flag an issue with the sampling of the posterior
we can be confident that something is amiss with the model, although the
converse is not immediately true; a lack of warnings does not imply the model
is behaving appropriately.

We compare the priors and posteriors produced by our methodology with the
posteriors produced using a flat, improper prior as a benchmark. We consider,
separately, each of the 93 individuals in the growth data (Tuddenham and
Snyder, 1954) provided by the fda package (Ramsay _et al._ , 2022) in R (R
Core Team, 2022). This is a form of prior sensitivity analysis, but distinct
from the ideas of Roos _et al._ (2015) which consider only one particular
realisation of the data. By considering each individual in the growth data
separately, as opposed to jointly in a hierarchical model, we heighten the
importance of including appropriate prior information. We sample each
posterior using Stan (Stan Development Team, 2021), setting adapt_delta = 0.95
and max_treedepth = 12 to minimise false positive warning messages.

Hartmann _et al._ (2020) also considered this example, so we compare also to
these results. However, there are key differences between our approaches that
must be kept in mind when comparing results. Hartmann _et al._ (2020) elicit 6
predictive quantiles at ages $t=(0,2.5,10,17.5)$, as opposed to entire
predictive distributions at ages $t=(2,8,13,18)$ which underpin the covariate-
specific version of our method. We use different ages because the model of
Preece and Baines (1978) is stated to be accurate and robust for ages greater
than 2. Hartmann _et al._ (2020) include a noise parameter in their definition
of $\theta$. The exact interpretation of this parameter is complicated by
their choice of Weibull likelihood, rendering the distribution of the
measurement errors sensitive to conditional mean of the model (this is still
the case despite their choice of Weibull parameterisation). Finally, Hartmann
_et al._ (2020) elicit quantiles from 5 different users and report an
estimated $\lambda^{*}$ for each user. These estimates (reproduced in Appendix
G.3) allow us to compare optimal the selected priors
$\text{p}(\theta\mid\lambda^{*})$. They do not report whether each user’s
estimate is consistent over repeated runs of their optimisation algorithm, and
do not discuss the issue of estimate replicability.

###  5.3 Results

####  5.3.1 Choosing $\kappa$

The optimal choice of $\kappa$ is target specific, and so we separately choose
an appropriate $\kappa\in\mathcal{K}=\\{0.05,0.1,\ldots,0.5\\}$ for the
covariate-independent and covariate-specific targets. As a heuristic, we
choose the value of $\kappa$ that yields the minimum variability of
$L(\lambda\mid\boldsymbol{X})$ for $\lambda\in\mathcal{P}$ amongst the
replicates, though this heuristic is unsuitable if, as we expect to be more
commonly the case, only one run of the optimiser is made. In such settings we
recommend plotting the Pareto frontiers as in Figure 10 and visually checking
that the minimum loss point is not at either extrema of the frontier.

Figure 10 displays Pareto frontiers for a selection of values
$\kappa\in\\{0.1,0.2,0.3,0.4\\}$, and the associated minimum loss points.
There is notable inter-replicate variability, with the Pareto frontier for
some replicates being totally dominated by other replicates. This is due to
the stochastic properties of the global optimisers we employ.

We select $\kappa=0.2$ for the covariate-specific target given our minimum
variability heuristic. Visible in most replicates is an inflection point at
values of $N(\lambda)\approx 1.5$ (the Y-axis in Figure 10) around which the
minimum loss points cluster, suggesting any of these points likely admits a
reasonable value for $\lambda^{*}$. The results for the covariate-independent
target are similar (See Appendix G Figure 22) and there we select
$\kappa=0.15$.

Figure 10: Pareto frontiers for each $\kappa\in\mathcal{K}$ for the covariate-
specific example. The minimum loss point for each replicate is plotted with
$\color[rgb]{0.70703125,0.05859375,0.125}\definecolor[named]{pgfstrokecolor}{rgb}{0.70703125,0.05859375,0.125}{+}$.

####  5.3.2 Replicability of discrepancy

Figure ~11 indicates that, across replicates, the optimisation procedure
identifies optimal $\lambda^{*}$ values with reasonably, but not entirely,
consistent predictive discrepancies.

Figure 11: Final predictive discrepancy $\log(D(\lambda^{*}))$, or
$\log(D(\lambda^{*}\mid\boldsymbol{X}))$ for the covariate-specific target.
The multiple objective optimisation approach uses the previously selected
values of $\kappa=0.15$ for the covariate-independent target, and $\kappa=0.2$
for the covariate-specific target. Horizontal jitter has been applied for
readability.

####  5.3.3 Faithfulness

Figure 12 displays the target and prior predictive density estimates in the
covariate-independent case. In the covariate-independent target case, we see
that introducing the secondary objective (right panel) produces estimates of
$\lambda^{*}$ that are congruent with estimates from the single objective
case, though the estimates are more variable. Both single and multi-objective
approaches result in reasonably, but not entirely, faithful densities for
$\text{p}(Y\mid\lambda^{*})$. However, most optimum priors seem to accumulates
additional probability surrounding $Y=h_{1}\approx 155$, resulting in
individual trajectories attaining their adult height $h_{1}$ for younger than
expected ages $t$ (which we will later confirm in Figure 14).

For the covariate-specific target, displayed in Figure 13, the secondary
objective introduces in some replicates outlying estimates for
$\text{p}(Y\mid\lambda^{*},X_{r})$, most clearly visible for $X_{1}=2$ and
$X_{2}=8$. Both the single and multi-objective approaches struggle to match
the prior predictive distribution at all ages, with consistently poorer
faithfulness for $X_{1}=2$. Empirically, it does not seem possible to match
all four supplied target prior predictive distributions simultaneously, given
the mathematical structure of the model. Lastly, because $\text{t}(Y\mid
X_{1}=2)$ is substantially narrower than the other targets, it is optimal,
under the Cramér-Von Mises discrepancy, to select wider priors better matching
the older age target distributions.

Figure 12: The covariate-independent marginal target density $\text{t}(Y)$
(red) and prior predictive densities $\text{p}(Y\mid\lambda^{*})$ for each of
the 30 replicates (blue lines). The replicates in the right panel are obtained
after $\kappa=0.15$ is chosen. Figure 13: Covariate-specific target densities
$\text{t}(Y\mid X_{r})$ (red lines) and prior predictive densities
$\text{p}(Y\mid\lambda^{*},X_{r})$ for each of the 30 replicates (blue lines).
The replicates in the right column are obtained after $\kappa=0.2$ is chosen.

####  5.3.4 Mean prior growth trajectories

We now explore whether the priors we estimate produce reasonable and
appropriately uncertain mean growth trajectories _a priori_. This is similar
to our requirement that priors are faithful to the supplied information, but
is more general as it considers all possible values of $t$. Figure 14 suggests
that both the covariate-independent and covariate-specific targets yield
plausible mean growth trajectories. However, the covariate-independent priors
are significantly more uncertain, resulting in implausible heights having _a
priori_ support. This contrasts with the covariate-specific priors, which
interpolate between the supplied targets with an acceptable degree of
uncertainty. It is interesting to note that the covariate-specific target has
similar levels of uncertainty across all ages, further suggesting that the
model may not be flexible enough to simultaneously match all the targets when
they have varying variance, as here. All 5 of the priors from Hartmann _et
al._ (2020), for a narrower uncertainty interval, are implausible in both
shape and width when viewed on this scale. It also seems unlikely that these
priors accurately reflect the information provided by the experts in Hartmann
_et al._ (2020), but this information is not reported.

Figure 14: Prior predictive for the model without noise
$\text{p}(h(t;\theta)\mid\lambda^{*})$ for each replicate/user from Hartmann
et al. (top right panel), the covariate-independent target (middle row) and
covariate-specific target (bottom row) for the multi-objective and single
objective settings (left column and right column respectively). Note that this
quantity does not include measurement error. Solid lines depict the mean, with
the grey regions representing the 95% prior predictive intervals, except for
the Hartmann panel, where the intervals are only 75% wide for visualisation
purposes. The y-axis is truncated to $(70,200)$ and uncertainty intervals are
also truncated to this range. The red lines in the covariate row correspond to
our supplied $\text{t}(Y\mid X_{r})$ densities, and represent the same
information as in Figure 13.

####  5.3.5 Posterior stability

Figure 15 displays for each of the 93 individuals in the growth data whether
Stan emits a warning message when estimating the posterior. The flat prior
consistently produces posteriors that emit warnings, with some individuals
particularly prone to warning messages (i.e. warnings are very correlated
within individual columns), suggesting that their data are less informative
than other individuals. Warnings are correlated within rows for the Hartmann
et. al. priors, indicating that some of the priors are more suitable
(replications 1 and 5) for the individuals in the growth data. Amongst our
results we note that the covariate-specific approach produces fewer warnings
than the covariate-independent approach in both the single- or multi-objective
cases. This reflects the additional information available in the covariate-
specific setting, and that this information results in more plausible priors.
Warnings are particularly correlated within specific priors (i.e. across rows)
for the covariate-independent approach, suggesting that these priors are
inappropriate for many individuals. The multi-objective approach (third row of
Figure 15) produces a small number of additional warnings above the equivalent
single-objective approach (fifth row of Figure 15), illustrating the trade-off
between single-objective priors that are as informative as possible, and
multiple-objective priors that are slightly less informative and thus perform
less well in settings where strong prior information is required.

Figure 15: Presence/absence of Stan warnings for all individuals (columns) in
the FDA package growth data and replicate prior estimates (rows). Each
replicate corresponds to a run of the optimisation process and thus a
different prior, except the flat, improper prior which is identical for each
replicate.

####  5.3.6 Uniqueness of marginal priors and posteriors

Figure 16 shows the priors for $(h_{0},\delta_{s})\in\theta$ for 30 replicates
with the covariate-specific target. The priors for our covariate-specific
target exhibit substantial variability. There appear to be two distinct
unimodal priors for $h_{0}$ with similar loss, suggesting that
$\text{T}(Y\mid\boldsymbol{X})$ does not provided enough information to
uniquely determine a prior distribution. However both priors are significantly
broader than the Hartmann et. al. priors. The marginal priors are slightly
more consistent for $\delta_{s}$ but variability persists.

Figure 16 also shows the posteriors for these parameters when using the
(uninformative and thus challenging) data from individual $n=26$. The flat
prior produces a multimodal posterior555We run Stan with the default 4 chains
to detect convergence warnings for Figure 15, but in Figure 16 we plot only
one chain per replicate to better highlight multi-modal posteriors. for our
parameters of interest, demonstrating the lack of stability when computing the
posterior and the need for some prior information. As a result, the posteriors
are under all methods strongly depend on the prior distribution used, which as
noted is not stable under any method here. However, all the possibility of a
posterior for $\delta_{s}$ with significant mass above 2 is removed by both of
our priors, as is desirable, because such solutions correspond to extremely
fast growth spurts that are physiologically implausible and are unsupported by
the data from individual $n=26$.

The priors and posteriors for the complete set of parameters $\theta$ are
displayed in Appendix G [Figures 23 and 24].

Figure 16: A comparison of the priors (blue) produced for $h_{0}$ by our
method using the covariate-specific target (top right) using the multiple
objective function ($\kappa=0.2$) and the single objective function
($\kappa=\text{NA}$, bottom right); Hartmann et al. (2020) (bottom left), and
with no prior displayed for the flat prior scenario (top left). The
corresponding posteriors for individual $n=26$ under each of these priors are
displayed in (red).

###  5.4 Example summary

The priors estimated by our procedure in this example are broadly faithful to
the supplied information, but appear to be constrained by the model’s
inflexibility that makes it difficult to simultaneously match the $t=2$ and
$t=18$ targets in the covariate-specific case. They regularise the posterior
sufficiently enabling accurate posterior sampling (model inadequacy
notwithstanding), with the covariate-specific, multi-objective method proving
most useful, but are arguably over concentrated and occasionally prevent the
model from fitting the data well. Uniqueness is improved by our secondary
objective, but perfect uniqueness across all replicates remains illusive and
may not be possible with only the information provided in
$\text{T}(Y\mid\boldsymbol{X})$. We observe a small improvement in uniqueness
attributable to the secondary objective (see Appendix G [Figures 23 and 24]);
some of the variability may be a result of lack of perfect replicability of
the optimisation rather than non-uniqueness.

##  6 Conclusion

In this paper we develop methodology, and software, for specifying priors
given predictive information about an observable or model-derived quantity. We
employ a CDF-based, multi-objective global optimisation approach to this
translation problem to make our approach widely applicable. Adopting a global
optimisation approach allows any kind of model to be specified and optimised,
not just those for which we can compute reparameterisation gradients. The
global optimisation approach also allows us to provide our functionality
directly in R, with which we envisage the majority of the users of our method
will be familiar. Our CDF-based predictive discrepancy is also generic as it
permits continuous, discrete, and mixed observable types. We apply our
methodology in three challenging example models, each of which we interrogate
for faithfulness, identifiability, and uniqueness. Each example is of moderate
dimension (3–17), with each representing a difficult structural elicitation
problem. Finally, our delineation between elicitation and translation, with
our emphasis on the latter, is a contribution to an under-explored area of
prior specification.

Our inspiration, for both methodology and examples, arises from applications
and models we have encountered in applied work. The Preece-Baines model is a
typical complex, nonlinear regression model for an intuitive and well
understood observable quantity, but for which the model parameters are not
easily understood. Setting a prior for models congruent with our knowledge is
difficult without a method for translation such as we have proposed. We
previously considered a survival model similar to the cure fraction model,
where we knew _a priori_ the fraction of cured/censored observations and a
distribution of likely survival times, in our earlier work (Manderson and
Goudie, 2022). Setting an appropriate prior for this model proved challenging,
and would have benefited greatly from the translation methodology introduced
in this paper. Finally, prior knowledge on $R^{2}$ has proved a valuable
mathematical basis for the R2-D2 shrinkage prior (Zhang _et al._ , 2022), and
more generally there are numerous model-derived quantities about which
practitioners possess prior information. Methods for translation are valuable
in this setting, as information about such model-derived quantities is often
difficult to express. An envisaged future example of this type considers
clustering models. In that setting we elicit an informative prior for number
of clusters, or the typical size of a cluster, which are derived from the
clustering model. Such quantities are readily reasoned about by experts, as
opposed to the parameters governing each cluster. Including this type of
information in complex models seems critical for stable and reliable Bayesian
inference.

One limitation of the current work is that we only partly address non-
uniqueness for flexible models and uninformative target distributions. In
these settings we emphases that our methodology remains valuable as a means to
assess the implications of a particular target distribution on a specific
model. For specific model-target pairs where uniqueness remains challenging,
our methodology can still provide useful insight into consequences of certain
$\text{T}(Y\mid\boldsymbol{X})$. We can delineate between components of
$\lambda$ that are well identified and thus consistently estimated, and those
that are not. Furthermore, we can directly inspect $\text{T}(Y\mid X_{r})$
different significantly from $\text{P}(Y\mid\lambda^{*},X_{r})$, and consider
if such differences are attributable to model inflexibility or implausible
targets. Finally, we are compelled to assess our choice of which components
should make up $\lambda$, and whether we have other information that we could
employ to fix certain components within $\lambda$ (e.g. the fixed prior for
the noise in the human height example). Another limitation that should be
noted is that global optimisation methods lack guarantees of finding the
global optimum in finite time, so the the performance of this optimisation
process to other problems is unclear.

## Acknowledgments and data availability

We thank Daniela De Angelis and Mevin Hooten for their feedback on an earlier
version of this manuscript.

This work was supported by The Alan Turing Institute under the UK Engineering
and Physical Sciences Research Council (EPSRC) [EP/N510129/1] and the UK
Medical Research Council [programme codes MC_UU_00002/2 and MC_UU_00002/20].
No original data were generated as part of this study, the growth data used in
Section 5 are available as part of the fda package for R (Ramsay _et al._ ,
2022) available on CRAN (https://cran.r-project.org/). For the purpose of open
access, the author has applied a Creative Commons Attribution (CC BY) licence
to any Author Accepted Manuscript version arising.

## Appendices

###  A Importance sampling

Appropriate importance distributions are crucial to obtaining an accurate and
low variance estimate of $D(\lambda\mid\boldsymbol{X})$. For values of
$\lambda$ far from optimal, $\text{P}(Y\mid\lambda,\boldsymbol{X})$ can differ
considerably from $\text{T}(Y\mid\boldsymbol{X})$. Given a specific $X_{r}$ we
require an importance distribution $\text{Q}(Y\mid X_{r})$ that places
substantial mass in the high probability regions of both $\text{T}(Y\mid
X_{r})$ and $\text{P}(Y\mid\lambda,X_{r})$, as it is in these regions that
$d(\cdot,\cdot)$ is largest. But we cannot exert too much effort on finding
these densities as they are specific to each value of $\lambda$, and must be
found anew for each $\lambda$.

We use three quantities to guide our choice of $\text{Q}(Y\mid X_{r})$, these
being the support $\mathcal{Y}$, the samples
$\boldsymbol{y}_{r}^{(\text{P})}\sim\text{P}(Y\mid\lambda,X_{r})$, and the
samples $\boldsymbol{y}_{r}^{(\text{T})}\sim\text{T}(Y\mid X_{r})$. Of primary
concern is the support. If $\mathcal{Y}=\mathbb{R}$ then we use a mixture of
Student-$t_{5}$ distributions; for $\mathcal{Y}=\mathbb{R}=(0,\infty)$ we
employ a mixture of gamma distributions; and for $\mathcal{Y}=(0,a]$ with
known $a$, we opt for a mixture of Beta distributions with a discrete
component at $Y=a$. The parameters of the mixture components are estimated
using the method of moments. Specifically, denoting the empirical mean of
$\boldsymbol{y}_{r}^{(\text{P})}$ as $\hat{\mu}^{(\text{P})}$ and the
empirical variance by $\hat{v}^{(\text{P})}$, with $\hat{\mu}^{(\text{T})}$
and $\hat{v}^{(\text{T})}$ defined correspondingly for
$\boldsymbol{y}_{r}^{(\text{T})}$, Table 1 details our method of moments
estimators for the mixture components.

In this paper we limit ourselves to one dimensional $\mathcal{Y}$, where
importance sampling is mostly well behaved or can be tamed using a reasonable
amount of computation. This covers many models, and with the covariate-
specific target it includes regression models. It is harder to elicit
$\text{T}(Y\mid\boldsymbol{X})$ for higher dimensional data spaces, and the
difficulties with higher dimensional importance sampling are well known.

$\mathcal{Y}$ | $\text{Q}_{r}(Y)$ | Parameter estimates | Mixture weights |  Notes
---|---|---|---|---
$\mathbb{R}$ | $\begin{aligned} &\pi_{1}\text{Student-}t_{5}(Y;\hat{\mu}_{1},\hat{s}_{1})+\\\ &\pi_{2}\text{Student-}t_{5}(Y;\hat{\mu}_{2},\hat{s}_{2})\end{aligned}$ | $\begin{aligned} \hat{\mu}_{1}=\hat{\mu}^{(\text{P})},&\,\hat{s}_{1}=c\sqrt{\hat{v}^{(\text{P})}}\\\ \hat{\mu}_{2}=\hat{\mu}^{(\text{T})},&\,\hat{s}_{2}=c\sqrt{\hat{v}^{(\text{T})}}\\\ \end{aligned}$ | $\pi_{1}=\pi_{2}=0.5$ |  $c$ defaults to $1.05$
$(0,\infty)$ | $\begin{aligned} &\pi_{1}\text{Gamma}(Y;\hat{\alpha}_{1},\hat{\beta}_{1})+\\\ &\pi_{2}\text{Gamma}(Y;\hat{\alpha}_{2},\hat{\beta}_{2})\end{aligned}$ | $\begin{aligned} \hat{\alpha}_{1}=\frac{(\hat{\mu}^{(\text{P})})^{2}}{\tilde{\omega}^{(\text{P})}},&\,\hat{\beta}_{1}=\frac{\hat{\mu}^{(\text{P})}}{\tilde{\omega}^{(\text{P})}}\\\ \hat{\alpha}_{2}=\frac{(\hat{\mu}^{(\text{T})})^{2}}{\tilde{\omega}^{(\text{T})}},&\,\hat{\beta}_{2}=\frac{\hat{\mu}^{(\text{T})}}{\tilde{\omega}^{(\text{T})}}\\\ \end{aligned}$ | $\pi_{1}=\pi_{2}=0.5$ |  $\tilde{\omega}=\min(c^{2}\hat{v},10^{5})$,
$c$ defaults to $1.05$
$(0,a)$ | $\begin{aligned} &\frac{\pi_{1}}{a}\text{Beta}\left(\frac{Y}{a};\hat{a}_{1},\hat{b}_{1}\right)+\\\ &\frac{\pi_{2}}{a}\text{Beta}\left(\frac{Y}{a};\hat{a}_{2},\hat{b}_{2}\right)+\\\ &\pi_{3}\mathbbm{1}_{\\{Y=a\\}}\end{aligned}$ | $\begin{aligned} \hat{a}_{1}&=\hat{\mu}^{(\text{P})}\left[\frac{\hat{\mu}^{(\text{P})}}{\tilde{\omega}^{(\text{P})}}(1-\hat{\mu}^{(\text{P})})-1\right]\\\ \hat{b}_{1}&=\frac{(1-\hat{\mu}^{(\text{P})})}{\hat{\mu}^{(\text{P})}}\hat{a}_{1}\\\ \hat{a}_{2}&=\hat{\mu}^{(\text{T})}\left[\frac{\hat{\mu}^{(\text{T})}}{\tilde{\omega}^{(\text{T})}}(1-\hat{\mu}^{(\text{T})})-1\right]\\\ \hat{b}_{2}&=\frac{(1-\hat{\mu}^{(\text{T})})}{\hat{\mu}^{(\text{T})}}\hat{a}_{2}\end{aligned}$ | $\begin{gathered}\pi_{1}=\pi_{2}=0.45\\\ \pi_{3}=0.05\end{gathered}$ |  $\tilde{\omega}=\max(c^{2}\hat{v},10^{-6})$,
$c$ defaults to $1.05$
Table 1: Importance distributions and method of moments estimators for their
constituent parametric distributions. Note that $c$ is a user-selected tuning
parameter to enable the construction of wider importance distributions.

###  B Evaluating $D(\lambda\mid\boldsymbol{X})$

For both numerical stability and optimisation performance (Eriksson and
Poloczek, 2021; Snoek _et al._ , 2014) we evaluate
$D(\lambda\mid\boldsymbol{X})$ on the log scale. This is because far from
optimal values of $\lambda$ have corresponding $D(\lambda\mid\boldsymbol{X})$
many orders of magnitude larger than near optimal values of $\lambda$.
Furthermore, the Gaussian process approximation that underlies Bayesian
optimisation assumes constant variance, necessitating a log or log-like
transformation.

Suppose again that we sample
$\boldsymbol{y}_{r}^{(\text{P})}\sim\text{P}(Y\mid\lambda,X_{r})$, from which
we form the ECDF
$\hat{\text{P}}(Y\mid\lambda,X_{r},\boldsymbol{y}_{r}^{(\text{P})})$. We also
select an appropriate importance distribution $\text{Q}(Y\mid X_{r})$ and
density $\text{q}(Y\mid X_{r})$ using Appendix A, and sample importance points
$(y_{i,r})_{i=1}^{I_{r}}\sim\text{Q}(Y\mid X_{r})$. Define the intermediary
quantity $z(y_{i,r})$ as

$z(y_{i,r})=\log\left(d\left(\hat{\text{P}}(y_{i,r}\mid\lambda,X_{r},\boldsymbol{y}_{r}^{(\text{P})}),\text{T}(y_{i,r}\mid
X_{r})\right)\right)+\log\left(\text{t}(y_{i,r}\mid
X_{r})\right)-\log\left(\text{q}(y_{i,r}\mid X_{r})\right),$ (19)

and then rewrite Equation (7) to read

$\log(D(\lambda\mid\boldsymbol{X}))=-\log(R)+\log\left(\sum_{r=1}^{R}\exp\left\\{-\log(I_{r})+\log\left(\sum_{i=1}^{I_{r}}\exp\left\\{z(y_{i,r})\right\\}\right)\right\\}\right).$
(20)

All $\log(\sum\exp\\{\cdot\\})$ terms are computed using the numerically
stable form (Blanchard _et al._ , 2021).

Accurately evaluating $\log(d(\cdot,\cdot))$ in Equation (19) involves
managing the discrete nature of the ECDF (that it returns exactly zero or one
for some inputs), and using specialised functions for each discrepancy to
avoid issues with floating point arithmetic. We compute
$\log(d^{\text{CvM}}(\cdot,\cdot))$ using

$\log\left(d^{\text{CvM}}\left(\hat{\text{P}}(y_{i,r}\mid\lambda,X_{r},\boldsymbol{y}_{r}^{(\text{P})}),\text{T}(y_{i,r})\right)\right)=2\log\left(\left\lvert\hat{\text{P}}(y_{i,r}\mid\lambda,X_{r},\boldsymbol{y}_{r}^{(\text{P})})-\exp\\{\mathcal{T}(y_{i,r})\\}\right\rvert\right),$
(21)

where $\mathcal{T}(y_{i,r})=\log(\text{T}(y_{i,r}))$. The log-CDF (LCDF) is
often more numerically accurate for improbable values of $y_{i,r}$, and so our
methodology assumes that it is this LCDF form in which the target distribution
is supplied. However, because the ECDF can return exact zero/one values there
is no way to perform this computation on the log scale. We thus employ high
precision floating point numbers when exponentiating the LCDF values, using
Rmpfr (Maechler _et al._ , 2021), to avoid evaluating $\log(0)$.

For $\log(d^{\text{AD}}(\cdot,\cdot))$, additional care must be taken as the
denominator of $d^{\text{AD}}$ in Equation (2) tends to underflow to zero.
Thus we evaluate it using

$\begin{multlined}\log\left(d^{\text{AD}}\left(\hat{\text{P}}(y_{i,r}\mid\lambda,X_{r},\boldsymbol{y}_{r}^{(\text{P})}),\text{T}(y_{i,r})\right)\right)=\\\
2\log\left(\left\lvert\hat{\text{P}}(y_{i,r}\mid\lambda,X_{r},\boldsymbol{y}_{r}^{(\text{P})})-\exp\\{\mathcal{T}(y_{i,r})\\}\right\rvert\right)-\mathcal{T}(y_{i,r})-\texttt{log1mexp}(-\mathcal{T}(y_{i,r})),\end{multlined}\log\left(d^{\text{AD}}\left(\hat{\text{P}}(y_{i,r}\mid\lambda,X_{r},\boldsymbol{y}_{r}^{(\text{P})}),\text{T}(y_{i,r})\right)\right)=\\\
2\log\left(\left\lvert\hat{\text{P}}(y_{i,r}\mid\lambda,X_{r},\boldsymbol{y}_{r}^{(\text{P})})-\exp\\{\mathcal{T}(y_{i,r})\\}\right\rvert\right)-\mathcal{T}(y_{i,r})-\texttt{log1mexp}(-\mathcal{T}(y_{i,r})),$
(22)

where $\texttt{log1mexp}(x)=\log(1-\exp\\{-x\\})$ is implemented by the Rmpfr
package (Maechler, 2012). Such precision is necessary for improbably large
values of $y_{i,r}$ under T, as the CDF/LCDF often rounds to 1/0
(respectively). It is not always feasible to evaluate Equation (22) with
sufficient accuracy to avoid under/over-flow issues – it requires a high-
precision implementation of $\mathcal{T}(y_{i,r})$ for extreme $y_{i,r}$ and
many additional bits of precision for both $y_{i,r}$ and the result. In these
settings we revert to $\log(d^{\text{CvM}}(\cdot,\cdot))$.

###  C Further notes on choosing $\kappa$

Advantages of multi-objective optimisation are most immediately apparent when
the scales of our objectives differ markedly. Consider the equivalent
linearised approach, where we select $\kappa$ _before_ optimisation and
directly optimise $\tilde{L}(\lambda\mid\boldsymbol{X})$. It is generally not
possible to know the range of the values of
$\tilde{D}(\lambda\mid\boldsymbol{X})$ and
$\tilde{N}(\lambda\mid\boldsymbol{X})$ before optimisation. Selecting an
appropriate $\kappa$ without this knowledge is prohibitively difficult,
leaving only the computationally expensive trial-and-error approach – where we
re-run the optimiser for each new possible value of $\kappa$ – as a plausible
strategy for choosing $\kappa$. In contrast, given $\mathcal{P}$ it is
computationally trivial to recompute $\lambda^{*}$ for many possible values of
$\kappa$ _after_ optimisation (e.g. each panel of Figure 10 is trivial to
compute). We can thus select $\kappa$ in a problem-specific manner for
practically no additional computational cost to that of the multi-objective
optimiser. Note that the multi-objective optimisation approach is more
expensive that the linearised approach, but this additional cost is dwarfed by
the number of re-runs of the latter typically required to select $\kappa$.

###  D Algorithmic descriptions of the optimisation process

####  D.1 CRS2 as an initialiser for Bayesian optimisation

Algorithm 3 describes our use of CRS2 (Kaelo and Ali, 2006) to obtain a
suitable design to initialise the Bayesian multi-objective optimisation
approach in step 2.

Algorithm 3 Using CRS2 to find an initial design for Bayesian optimisation

1Log total predictive discrepancy $\log(D(\lambda\mid\boldsymbol{X}))$
(evaluable using Algorithm 1), number of CRS2 iterations to run
$N_{\text{CRS2}}$, number of points in final design $N_{\text{design}}$,
number of additional padding points to add for numerical stability
$N_{\text{pad}}$, hyperparameter support $\Lambda$

2

3function Initial design($N_{\text{CRS2}},N_{\text{design}},N_{\text{pad}}$)

4 Initialise $\mathcal{S}=\\{\\}$, an empty set to hold possible design points

5 for $i$ in $1\ldots N_{\text{CRS2}}$ do

6 Minimising $\log(D(\lambda\mid\boldsymbol{X}))$, get the $i$th trial point
$\tilde{\lambda}_{i}$ and value $\log(D(\lambda_{i}\mid\boldsymbol{X}))$ from
CRS2 with local mutation (Kaelo and Ali, 2006)

7 Compute
$\tilde{w}_{i}=-\exp\left\\{\log\left(D(\tilde{\lambda}_{i}\mid\boldsymbol{X})\right)\right\\}$

8 Concatenate
$\mathcal{S}=\mathcal{S}\cup\left\\{\tilde{\lambda}_{i},\log(D(\tilde{\lambda}_{i}\mid\boldsymbol{X})),\tilde{w}_{i}\right\\}$

9 end for

10 Normalise weights such that
$w_{i}=\exp\left\\{\tilde{w}_{i}-\log\left(\sum_{i=1}^{N_{\text{CRS2}}}\exp\left\\{\tilde{w}_{i}\right\\}\right)\right\\}$

11 Subsample without replacement $N_{\text{design}}$ values from $\mathcal{S}$
according to the normalised weights, and store in
$\mathcal{D}=\\{\lambda_{i},\log(D(\lambda_{i}\mid\boldsymbol{X}))\\}_{i=1}^{N_{\text{design}}}$

12 Sample $N_{\text{pad}}$ points from a Latin hypercube design spanning
$\Lambda$ (Stein, 1987), evaluate $\log(D(\lambda\mid\boldsymbol{X}))$ at
these points, and add them to $\mathcal{D}$

13 return:
$\mathcal{D}=\\{\lambda_{i},\log(D(\lambda_{i}\mid\boldsymbol{X}))\\}_{i=1}^{N_{\text{design}}+N_{\text{pad}}}$

14end function

####  D.2 MSPOT

Algorithm 4 describes, in our notation, the MSPOT (Zaefferer _et al._ , 2012)
algorithm for two objectives. Note that within the algorithm we suppress each
objective’s dependence on $\boldsymbol{X}$ for brevity.

Algorithm 4 Global two-objective Bayesian optimisation using MSPOT (Zaefferer
_et al._ , 2012)

1Primary objective $D(\lambda)$, secondary objective $N(\lambda)$, initial
design
$\mathcal{D}=\left\\{\lambda_{i},D(\lambda_{i}),N(\lambda_{i})\right\\}_{i=1}^{N_{\text{design}}+N_{\text{pad}}}$,
number of iterations $N_{\text{BO}}$, number of new points to evaluate the
surrogate models at $N_{\text{new}}$, number of evaluations to add to the
design within an iteration $N_{\text{eval}}$, hyperparameter support $\Lambda$

2

3function Bayesian optimisation using MSPOT($N_{\text{BO}}$)

4 for $i$ in $1\ldots N_{\text{BO}}$ do

5 Form Gaussian process (GP) approximations to $D(\lambda)$ and $N(\lambda)$
using $\mathcal{D}$

6 Generate a new Latin hypercube design $\mathcal{N}$ of size $N_{\text{new}}$
covering $\Lambda$, such that $N_{\text{new}}\gg N_{\text{design}}$

7 for $k$ in $1\ldots N_{\text{new}}$ do

8 Use the GPs to estimate $\hat{D}(\lambda_{k})$ and $\hat{N}(\lambda_{k})$

9 Add these to $\mathcal{N}$ so that
$\mathcal{N}_{k}=\left\\{\lambda_{k},\hat{D}(\lambda_{k}),\hat{N}(\lambda_{k})\right\\}$

10 end for

11 Truncate $\mathcal{N}$ to $N_{\text{eval}}$ points according to the non-
dominated sorting rank and hypervolume contribution (Beume _et al._ , 2007;
Deb, 2001; Deb _et al._ , 2002; Beume _et al._ , 2009) of each point in
$\left\\{D(\lambda_{k}),N(\lambda_{k})\right\\}_{k=1}^{N_{\text{new}}}$ with
$N_{\text{eval}}\ll N_{\text{new}}$

12 for $j$ in $1\ldots N_{\text{eval}}$ do

13 Evaluate the objectives $D(\lambda_{j})$ and $N(\lambda_{j})$ for
$\lambda_{j}\in\mathcal{N}$

14 Add these evaluations to
$\mathcal{D}=\mathcal{D}\cup\left\\{\lambda_{j},D(\lambda_{j}),N(\lambda_{j})\right\\}$

15 end for

16 end for

17 Compute the Pareto frontier
$\mathcal{P}=\left\\{\lambda_{i},D(\lambda_{i}),N(\lambda_{i})\right\\}_{i=1}^{\lvert\mathcal{P}\rvert}$
from
$\mathcal{D}=\left\\{\lambda_{i},D(\lambda_{i}),N(\lambda_{i})\right\\}_{i=1}^{N_{\text{design}}+N_{\text{pad}}+N_{\text{BO}}N_{\text{eval}}}$
(Kung _et al._ , 1975, see)

18 return: $\mathcal{P}$ and $\mathcal{D}$

19end function

####  D.3 Inter batch resampling

Algorithm 5 describes our inter-batch resampling algorithm that we
occasionally adopt in stage two of our optimisation process.

Algorithm 5 Resample the outputs from a previous batch to obtain a design for
the current one.

1Pareto frontier
$\mathcal{P}=\left\\{\lambda_{i},\log(D(\lambda_{i}\mid\boldsymbol{X})),N(\lambda_{i}\mid\boldsymbol{X})\right\\}_{i=1}^{\lvert\mathcal{P}\rvert}$
and all evaluated points
$\mathcal{E}=\left\\{\lambda_{i},\log(D(\lambda_{i}\mid\boldsymbol{X})),N(\lambda_{i}\mid\boldsymbol{X})\right\\}_{i=1}^{\lvert\mathcal{E}\rvert}$
from previous batch (with
$\lvert\mathcal{P}\rvert\ll\lvert\mathcal{E}\rvert$), number of design points
$N_{\text{design}}$, number of padding points $N_{\text{pad}}$, hyperparameter
support $\Lambda$

2

3function Next batch design($N_{\text{design}},N_{\text{pad}}$)

4 Initialise $\mathcal{D}=\mathcal{P}$

5 Compute the weights $w_{i}$ for all points in $\mathcal{E}$ in the same
manner as Algorithm 3 so that
$\mathcal{E}=\left\\{\lambda_{i},\log(D(\lambda_{i}\mid\boldsymbol{X})),N(\lambda_{i}\mid\boldsymbol{X}),w_{i}\right\\}_{i=1}^{\lvert\mathcal{E}\rvert}$

6 Sample without replacement
$\max\left(N_{\text{design}}-\lvert\mathcal{P}\rvert,0\right)$ points from
$\mathcal{E}$ according to the weights and add these points to $\mathcal{D}$

7 Sample $N_{\text{pad}}$ points from a Latin hypercube design covering
$\Lambda$ and add these to $\mathcal{D}$

8 return: $\mathcal{D}$ such that
$\lvert\mathcal{D}\rvert=\max(N_{\text{design}},\lvert\mathcal{P}\rvert)+N_{\text{pad}}$

9end function

###  E Additional information for the cure fraction survival example

####  E.1 Hyperparameter support $\Lambda$

See Table 2

Hyperparameter | Lower | Upper | # Elements
---|---|---|---
$\alpha$ | $\epsilon$ | 20 | 1
$\beta$ | $\epsilon$ | 20 | 1
$\mu_{0}$ | -10 | 10 | 1
$\sigma_{0}$ | $\epsilon$ | 10 | 1
$s_{\beta}$ | $\epsilon$ | 10 | 1
$\boldsymbol{\omega}$ | -1 + $\epsilon$ | 1 - $\epsilon$ | 6
$\boldsymbol{\eta}$ | -5 | 5 | 4
$a_{\pi}$ | 1 | 50 | 1
$b_{\pi}$ | 1 | 50 | 1
Table 2: Hyperparameters $\lambda$ for the cure fraction model, their upper
and lower limits that define $\Lambda$, and the number of elements in the
hyperparameter (which is 1 for all scalar quantities). Note that
$\epsilon=10^{-4}$ is added or subtracted to the limits to avoid degenerate
estimates for $\lambda$.

####  E.2 Pareto Frontiers for Cramér-von Mises discrepancy

See Figure 17.

Figure 17: Pareto frontiers for the survival example using the Cramér-von
Mises discrepancy, for the values of $\kappa$ we consider. Note that the
colour scale displaying loss is panel-specific. The red crosses
($\color[rgb]{0.70703125,0.05859375,0.125}\definecolor[named]{pgfstrokecolor}{rgb}{0.70703125,0.05859375,0.125}{+}$)
indicate the minimum loss point on each frontier, for each value of $\kappa$.

###  F Additional information for the $R^{2}$ example

####  F.1 Hyperparameter support $\Lambda$ – faithfulness experiment

See Table 3. Note that for the Dirichlet-Laplace prior, Zhang and Bondell
(2018) suggest bounding $\alpha\in[(\max(n,p))^{-1},1\mathop{/}2]$. In our
experiments we regularly encountered optimal values of $\alpha$ on the lower
boundary, so we use instead $1\mathop{/}(3\max(n,p))$ as a lower bound.

Prior | Hyperparameter | Lower | Upper
---|---|---|---
Gaussian | $a_{1}$ | 2 | 500
Gaussian | $b_{1}$ | 0.2 | 500
Gaussian | $\gamma$ | 1 | 500
Dir. Lap. | $a_{1}$ | 2 | 500
Dir. Lap. | $b_{1}$ | 0.2 | 500
Dir. Lap. | $\alpha$ | $1\mathop{/}(3\max(n,p))$ | $1\mathop{/}2$
Reg. Horse. | $a_{1}$ | 2 | 500
Reg. Horse. | $b_{1}$ | 0.2 | 500
Reg. Horse. | $p_{0}$ | 1 | $p\mathop{/}2$
Reg. Horse. | $\nu$ | 1 | 80
Reg. Horse. | $s^{2}$ | $10^{-5}$ | 100
Table 3: Hyperparameters $\lambda$ for the $R^{2}$ example and their
upper/lower limits that define $\Lambda$.

####  F.2 A comparison to an asymptotic result

The poor fit for the Gaussian prior observed in Figure 8 could be attributed
to issues in the optimisation process, or to the lack of flexibility in the
prior. To investigate, we compare the results for $\lambda_{\text{GA}}$ to
Theorem 5 of Zhang and Bondell (2018), which is an asymptotic result regarding
the optimal value of $\lambda_{GA}$ for a target $\text{Beta}(s_{1},s_{2})$
density for $R^{2}$. We compare pairs of $(n_{k},p_{k})$ for $k=1,\ldots,5$,
noting that assumption (A4) of Zhang and Bondell requires that
$p_{k}=\text{o}(n_{k})$ as $k\rightarrow\infty$ (for strictly increasing
sequences $p_{k}$ and $n_{k}$). Thus we consider values of $p$ such that
$p_{1}=80$ with $p_{k}=2p_{k-1}$ and $n$ with $n_{1}=50$ and
$n_{k}=n_{k-1}^{1.2}$, both for $k=2,\ldots,5$. Each $(n_{k},p_{k})$ pair is
replicated 20 times, and for each replicate we generate a different
$\boldsymbol{X}$ matrix with standard normal entries. As the target density we
choose $s_{1}=5,s_{2}=10$ – a “more Gaussian” target than previously
considered and thus, we speculate, possibly more amenable to translation with
a Gaussian prior for $\beta$. We also use this example as an opportunity to
assess if there are notable differences between the Cramér-Von Mises
discrepancy and the Anderson-Darling discrepancy as defined in Equation (2).
The support $\Lambda$ for $\lambda_{\text{GA}}$ differs slightly from the
example in the main text, and is defined in Table 4, as matching our target
with larger design matrices requires considerably larger values of $\gamma$.

The computation of $R^{2}$ becomes increasingly expensive as $n_{k}$ and
$p_{k}$ increase, which limits the value of some of our method’s tuning
parameters. The approximate discrepancy function uses $S=2000$ samples from
the prior predictive and is evaluated using $I=500$ importance samples. We run
CRS2 for $N_{\text{CRS2}}=500$ iterations, using $N_{\text{design}}=50$ in the
initial design for the subsequent single batch of Bayesian optimisation, which
uses $N_{\text{BO}}=100$ iterations.

Prior | Hyperparameter | Lower | Upper
---|---|---|---
Gaussian | $a_{1}$ | $2+10^{-6}$ | 50
Gaussian | $b_{1}$ | 0.2 | 50
Gaussian | $\gamma$ | 1 | 5000
Table 4: Hyperparameters $\lambda$ for the asymptotic example, and their
upper/lower limits that define $\Lambda$.

##### Results

Figure 18 displays the results in terms of the normalised difference between
the $\gamma$ we estimate $\gamma_{\text{pbbo}}^{*}$, and the asymptotic result
of Zhang and Bondell $\gamma_{\text{asym}}^{*}$. Our typical finite sample
estimate is slightly larger than the asymptotic result, and the difference
increases with $n_{k}$ and $p_{k}$. The variability of the normalised
difference remains roughly constant, and thus reduces on an absolute scale,
though extrema seem to occur more frequently for larger $n_{k}$ and $p_{k}$.
These simulations suggest that the asymptotic regime has not been reached even
at the largest $n_{k}$ and $p_{k}$ values we assessed.

Figure 18: Relative difference between the value of $\gamma$ obtained using
our methodology ($\gamma_{\text{pbbo}}^{*}$) and Theorem 5 of Zhang and
Bondell (2018) ($\gamma_{\text{asym}}^{*}$).

The estimates of $\gamma$ are not themselves particularly illuminating: we
should instead look for differences in the distribution of $R^{2}$ at the
optima, which is to say on the “data” scale. Figure 19 displays the target
distribution and the prior predictive distribution at the optima
$\text{p}(R^{2}\mid\lambda^{*}_{GA})$. The fit is increasingly poor as $n$ and
$p$ increase, and there is little difference both between the two
discrepancies and with each discrepancies replications. The lack of difference
implies that the optimisation process is consistently locating the same minima
for $D(\lambda)$. We conclude that either 1) the ability of the model to match
the target depends on there being additional structure in $\boldsymbol{X}$, or
2) it is not possible to encode the information in a $\text{Beta}(5,10)$ prior
for $R^{2}$ into the Gaussian prior.

Figure 19: The target density $\text{t}(R^{2})$ and optimal prior predictive
densities $\text{p}(R^{2}\mid\lambda^{*})$ under both the Cramér-von Mises
(red, left column) and Anderson-Darling (blue, right column) discrepancies.
There are 20 replicates of each discrepancy in this plot.

This example also further illustrates the difficulties inherent in acquiring a
prior for additive noise terms. Specifically, in this example it is difficult
to learn $(a_{1},b_{1})$, despite the fact that the contribution of
$\sigma^{2}$ to Equation (13) is not purely additive. However, as we see in
Figure 20, estimates are uniformly distributed across the permissible space,
except for bunching at the upper and lower bounds of $\Lambda$. Note that for
numerical and computational stability, we constrain $a_{1}\in(2,50]$ and
$b_{1}\in(0.2,50]$ in this example. This contrasts with similarity between
replicates visible in Figure 19, and is thus evidence that
$(\hat{a}_{1},\hat{b}_{1})$ have no apparent effect on the value of
$D(\lambda^{*})$. We should instead set the prior for $\sigma^{2}$ based on
external knowledge of the measurement process for $Y$.

Figure 20: Histograms of scaled estimates of $(a_{1}^{*},b_{1}^{*})$ for the
settings considered in Section F.2. Estimates have been scaled to $[0,1]$ for
visualisation purposes using the upper and lower limits defined in Table 3.

The regularisation method we employ in the two other examples in the main text
is unlikely to assist in estimating $(a_{1},b_{1})$. Promoting a larger mean
log marginal standard deviation, with the knowledge $D(\lambda)$ is
insensitive to the value of $(a_{1},b_{1})$, would simply pick the largest
possible value for $b_{1}^{2}\mathop{/}\left((a_{1}-1)^{2}(a_{1}-2)\right)$,
which occurs when $a_{1}$ is at its minimum allowable value and $b_{1}$ its
corresponding maximum.

####  F.3 Full faithfulness results

The complete results from the faithfulness experiment are displayed in Figure
21.

Figure 21: As in Figure 8 but for all values of $(s_{1},s_{2})$ denoted in the
facet panels titles. The performance of the regularised horseshoe is superior
to the Dirichlet-Laplace, both of which are vast improvements over the
Gaussian.

###  G Additional information for the Preece-Baines example

####  G.1 Details for $\text{T}(Y)$ and $\text{T}(Y\mid X_{r})$

Denote with $\text{Gamma}(Y;\alpha,\beta)$ the CDF of the gamma distribution
with shape parameter $\alpha$ and rate $\beta$; and
$\text{Normal}(Y;\xi,\omega^{2})$ the CDF of the normal distribution with mean
$\xi$ and standard deviation $\omega$. We define the covariate-independent
target

$\text{T}(Y)=0.38\,\text{Gamma}(Y;45.49,0.44)+0.36\,\text{Gamma}(Y;115.41,0.81)+0.27\,\text{Gamma}(Y;277.51,1.64),$
(23)

and the covariate-specific target

$\begin{gathered}\text{T}(Y\mid
X_{1}=2)=\text{Normal}(Y;88,3.5^{2}),\quad\text{T}(Y\mid
X_{2}=8)=\text{Normal}(Y;130,5.5^{2}),\\\ \text{T}(Y\mid
X_{3}=13)=\text{Normal}(Y;160,8^{2}),\quad\text{T}(Y\mid
X_{4}=18)=\text{Normal}(Y;172,9.5^{2}).\end{gathered}$ (24)

####  G.2 Hyperparameter support $\Lambda$

Table 5 contains the upper and lower limits for each hyperparameter, thus
defining the feasible region $\Lambda$.

Parameter ($\theta_{q}$) |  $\mu_{q}$ – Lower |  $\mu_{q}$ – Upper |  $\sigma_{q}$ – Lower |  $\sigma_{q}$ – Upper
---|---|---|---|---
$\theta_{1}=h_{0}$ | 130 | 185 | $\epsilon$ | 30
$\theta_{2}=\delta_{h}$ | $\epsilon$ | 30 | $\epsilon$ | 2
$\theta_{3}=s_{0}$ | $\epsilon$ | 0.2 | $\epsilon$ | 0.1
$\theta_{4}=\delta_{s}$ | $\epsilon$ | 1.5 | $\epsilon$ | 0.2
$\theta_{5}=\gamma$ | 9 | 15 | $\epsilon$ | 1
Table 5: Parameter vector $\theta$ and associated model specific parameter.
The rightmost four columns of the table define the upper and lower limits for
the hyperparameters $(\mu_{q},\sigma_{q})$, thus defining $\Lambda$.
Informative bounds are required for numerical stability of the data generating
process, and an $\epsilon=10^{-6}$ is required to avoid nonsensical optimal
values of $\lambda$.

####  G.3 Hartmann _et al._ (2020) priors

Table 6 contains the priors elicited by Hartmann _et al._ (2020) (extracted
from the supplementary material of that paper) for the parameters in the
Preece-Baines example. To generate the prior predictive samples displayed in
Figure 14 we draw, for each user, $\theta$ from the corresponding lognormal
distribution then compute $h(t;\theta)$ using Equation (18) (without the error
term) and 250 values of $t$ spaced evenly between ages $2$ and $18$.

User | Parameter | Expectation | Variance | Lognormal $\mu$ | Lognormal $\sigma$
---|---|---|---|---|---
1 | $h_{0}$ | 162.80 | 4.20 | 5.09 | 0.01
1 | $h_{1}$ | 174.50 | 0.80 | 5.16 | 0.01
1 | $s_{0}$ | 0.10 | 0.10 | -3.50 | 1.55
1 | $s_{1}$ | 3.30 | 0.21 | 1.18 | 0.14
1 | $\theta$ | 13.40 | 0.01 | 2.60 | 0.01
2 | $h_{0}$ | 153.73 | 1.60 | 5.04 | 0.01
2 | $h_{1}$ | 191.74 | 4.32 | 5.26 | 0.01
2 | $s_{0}$ | 0.04 | 0.01 | -4.21 | 1.41
2 | $s_{1}$ | 2.00 | 4.30 | 0.33 | 0.85
2 | $\theta$ | 15.90 | 0.70 | 2.76 | 0.05
3 | $h_{0}$ | 148.80 | 1.86 | 5.00 | 0.01
3 | $h_{1}$ | 177.14 | 3.68 | 5.18 | 0.01
3 | $s_{0}$ | 0.07 | 0.00 | -2.75 | 0.43
3 | $s_{1}$ | 4.54 | 37.83 | 0.99 | 1.02
3 | $\theta$ | 11.31 | 0.21 | 2.42 | 0.04
4 | $h_{0}$ | 162.80 | 0.02 | 5.09 | 0.00
4 | $h_{1}$ | 174.50 | 0.01 | 5.16 | 0.00
4 | $s_{0}$ | 0.10 | 0.01 | -2.65 | 0.83
4 | $s_{1}$ | 1.60 | 1.70 | 0.22 | 0.71
4 | $\theta$ | 14.70 | 0.90 | 2.69 | 0.06
5 | $h_{0}$ | 162.60 | 0.85 | 5.09 | 0.01
5 | $h_{1}$ | 174.40 | 0.90 | 5.16 | 0.01
5 | $s_{0}$ | 0.10 | 0.01 | -2.65 | 0.83
5 | $s_{1}$ | 3.40 | 0.01 | 1.22 | 0.03
5 | $\theta$ | 14.60 | 0.02 | 2.68 | 0.01
Table 6: Priors elicited by Hartmann _et al._ (2020) for each of the 5 users
they study. Hartmann _et al._ provide their results in the form of expected
values and variances for the parameters of the model, we compute the
corresponding lognormal location $\mu$ and scale $\sigma$ parameters from this
information. Values are rounded to two digits of precision.

####  G.4 Pareto frontiers for the covariate-independent target

The Pareto frontier for the covariate-independent target and all values of
$\kappa\in\mathcal{K}$ is displayed in Figure 22.

Figure 22: Pareto frontiers for each $\kappa\in\mathcal{K}$ for the covariate-
independent example. The minimum loss point for each replicate is plotted with
$\color[rgb]{0.70703125,0.05859375,0.125}\definecolor[named]{pgfstrokecolor}{rgb}{0.70703125,0.05859375,0.125}{+}$.

####  G.5 Full marginal prior and posterior comparison plots

Figures 23 and 24 are extended versions of Figure 16, and display the prior
and posterior estimates for all the parameters in $\theta$. Consistency and
uniqueness remain, evidently, challenging and as yet unobtainable.

Figure 23: A comparison of the priors (blue) produced by our method using the
covariate-independent marginal target (bottom two rows); and Hartmann et al.
(2020) (second row), with no prior displayed for the flat prior scenario. The
corresponding posteriors (red) for individual $n=26$ under each of these
priors are displayed as dashed lines. Note that y-axes change within columns
and are limited to values that clip some of the priors/posteriors for
readability. Figure 24: Otherwise identical to Figure 23 but the bottom two
rows display the results obtained using the covariate-specific target.

## References

  * Albert _et al._ (2012) Albert, I., Donnet, S., Guihenneuc-Jouyaux, C., Low-Choy, S., Mengersen, K., and Rousseau, J. (2012). “Combining expert opinions in prior elicitation,” Bayesian Analysis 7(3), 503–532, doi: 10.1214/12-BA717.
  * Amico and Van Keilegom (2018) Amico, M., and Van Keilegom, I. (2018). “Cure Models in Survival Analysis,” Annual Review of Statistics and Its Application 5(1), 311–342, doi: 10.1146/annurev-statistics-031017-100101.
  * Anderson and Darling (1952) Anderson, T. W., and Darling, D. A. (1952). “Asymptotic theory of certain "goodness of fit" criteria based on stochastic processes,” The Annals of Mathematical Statistics 23(2), 193–212, doi: 10.1214/aoms/1177729437.
  * Azzalini (2022) Azzalini, A. (2022). “sn: The skew-normal and related distributions such as the skew-t and the SUN” .
  * Azzalini and Valle (1996) Azzalini, A., and Valle, A. D. (1996). “The multivariate skew-normal distribution,” Biometrika 83(4), 715–726, doi: 10.1093/biomet/83.4.715.
  * Bedrick _et al._ (1997) Bedrick, E. J., Christensen, R., and Johnson, W. (1997). “Bayesian binomial regression: Predicting survival at a trauma center,” The American Statistician 51(3), 211–218, doi: 10.2307/2684890.
  * Beume _et al._ (2009) Beume, N., Fonseca, C. M., Lopez-Ibanez, M., Paquete, L., and Vahrenhold, J. (2009). “On the complexity of computing the hypervolume indicator,” IEEE Transactions on Evolutionary Computation 13(5), 1075–1082, doi: 10.1109/TEVC.2009.2015575.