sec0em2.75em subsec2.75em3.5em

# Optimal sizing of a holdout set for safe predictive model updating

Sami Haidar-Wehbe1,⟂, Samuel R. Emerson1,⟂, Louis J. M. Aslett1,2, James
Liley1
1: Department of Mathematical Sciences, Durham University, UK
2: Alan Turing Institute, London, UK
$\perp$: Equal contribution
Correspondence to {louis.aslett<EMAIL_ADDRESS>
###### Abstract

Predictive risk scores are increasingly used to guide clinical or other
interventions in complex settings, particularly healthcare. Directly updating
a risk score used to guide interventions leads to biased risk estimates. We
propose updating using a ‘holdout set’ – a subset of the population that does
not receive risk-score-guided interventions – to prevent this. Since samples
in the holdout set do not benefit from risk predictions, its size must trade
off performance of the updated risk score whilst minimising the number of held
out samples. We prove that this approach outperforms simple alternatives, and
by defining a general loss function describe conditions under which an optimal
holdout size (OHS) can be readily identified. We introduce parametric and
semi-parametric algorithms for OHS estimation and demonstrate their use on a
recent risk score for pre-eclampsia. Based on these results, we argue that a
holdout set is a safe, viable and easily implemented means to safely update
predictive risk scores. Holdout set; Machine learning in healthcare; Machine
learning safety; Model update; Predictive model; Treatment decision.

###  1 Introduction

Risk scores estimate the probability of an event $Y$ given predictors $X$.
Their use has become routine in medical practice (Topol, 2019), where $Y$
typically represents adverse event incidence and $X$ various clinical
observations. Once calculated, risk scores may be used to guide interventions,
perhaps modifying $X$, with the aim of decreasing the probability of an
adverse $Y$. For example, the QRISK3 score predicts thromboembolic risk given
predictors including age and hypertension (Hippisley-Cox et al., 2017), and a
high score may prompt prescription of antihypertensives.

Risk scores are typically developed by regressing observations of $Y$ on $X$.
Should the distribution of $\left(X,Y\right)$ subsequently change, or ‘drift’,
then risk estimates may become biased (Tsymbal, 2004). This can happen
naturally over time, meaning that risk scores typically need to be updated
periodically to maintain accuracy.

Updating of the risk score will involve obtaining new observations of $(X,Y)$.
Crucially the distribution of $(X,Y)$ may also change due to the effect of the
risk score itself: that is, high predicted risk of an adverse event may
trigger intervention to reduce that risk. The effect of such interventions may
be difficult or impossible to infer or measure, and indeed the fact that
intervention took place may be unrecorded. In the QRISK3 example above,
individuals prescribed antihypertensives in response to higher QRISK3 scores
should have lower thromboembolic risk than they would have if QRISK3 was not
used. Should a new risk score be fitted to observed $(X,Y)$, the effect of
hypertension on risk would be underestimated, and overall risk estimates
upwards-biased. This bias is worsened by heavier intervention resulting in
risk scores becoming ‘victims of their own success’ (Lenert et al., 2019).
This framework of directly updating a risk score on an ‘intervened’ population
has been termed ‘repeated risk minimisation’ (Perdomo et al., 2020) in the
context when such bias is accounted for or ‘naïve updating’ (Liley et al.,
2021b) when it is not.

In a previous paper we briefly noted that this bias could be avoided by
splitting the population on which the score can be used into an ‘intervention’
set and a ‘holdout’ set, with an updated model trained on the latter (Liley et
al., 2021b). In this work, we formally develop this proposal for practical use
in real-world predictive risk score updating and prove its suitability. We
address a vital tension in the choice of an optimal size (OHS) for the holdout
set: for the risk score to be accurate, the holdout set should not be too
small; but any samples in the holdout set will not benefit from risk scores,
so nor should it be too large.

Our work is methodological, and pre-empts a class of problems rather than
analysing an existing risk score. We show that the updating paradox above is
inevitable for risk scores on complex systems intended to guide interventions.
The holdout set approach provides a simple method to facilitate ongoing use of
such risk scores in the presence of both drift and intervention effects.

Contributions in this area are important due to rapidly evolving legislation.
Currently, the European Union treat each update of a risk model as a separate
risk score requiring re-approval, but in the USA a proactive approach is taken
with a ‘total-life cycle’ paradigm which allows practitioners to update risk
models as necessary without requesting approval (USFDA and others, 2019). This
approach could allow updating-induced biases to go undetected, and highlights
the need for safe updating methods in risk score deployment. The use of
holdout sets as examined in this work offers one potential solution.

Our paper is structured in the following way. In Section 2, we review relevant
literature and more precisely define the problem. We then give a motivating
example. In Section 4, we demonstrate dominance of the holdout-set approach
over simpler approaches, and construct an optimisation problem to find an OHS.
In Section 5, we describe two algorithms for OHS estimation, using explicit
parametrisation and using Bayesian emulation. In Section 6, we support our
findings with numerical demonstrations and resolve our motivating example by
applying our methods to a risk score for pre-eclampsia (PRE) to estimate an
OHS for updating it.

###  2 Review of related work

Widespread collection of electronic health records has spurred development of
new diagnostic and prognostic risk scores (Cook and Collins, 2015; Liley et
al., 2021a), which can allow detection of patterns too complex for humans to
discover. Examples of such scores in widespread use include: EuroSCORE II,
which predicts mortality risk at hospital discharge following cardiac surgery
(Nashef et al., 2012); and the STS risk score from the USA predicting risk of
postoperative mortality (Shahian et al., 2018). Many such scores have
demonstrable efficacy in clinical trials and in-vivo (Chalmers et al., 2013;
Wallace et al., 2014; Hippisley-Cox et al., 2017).

An important general concern with these scores is continued accuracy of
predictions. A 2011 review found that risk scores for hospital readmission
perform poorly and highlighted issues with design of their trials (Kansagara
et al., 2011). More recently, an analysis of a sepsis response score used
during the COVID pandemic found increasing risk overestimation over time
Finlayson et al. (2020). Various efforts have been made to standardise
procedures in risk score estimation to address these issues Collins et al.
(2015).

Several algorithms have been developed to update models with new data in the
presence of drift (Lu et al., 2018), which ideally leads to the best possible
model performance after every update. Adaptation of model updating to avoid
naïve updating-induced bias requires explicit causal reasoning (Sperrin et
al., 2019) and often further data collection (Liley, 2021). In a seminal
paper, Perdomo et al. (2020) analyse asymptotic behaviour of repeated naïve
updating, giving necessary and sufficient conditions for successive
predictions to converge to a stable setting where they ‘predict their own
effect’. Other approaches to optimise a general loss function by modulating
parameters of the risk score are developed in Mendler-Dünner et al. (2020);
Drusvyatskiy and Xiao (2020); Li and Wai (2021) and Izzo et al. (2021). These
approaches seek to minimise a ‘performative’ loss to the population in the
presence of an arbitrary risk score. Our alternative approach seeks to target
risk scores which reliably estimate the same quantity, namely $P(Y\mid X)$ in
a ‘native’ system prior to risk score deployment. Our approach is well-suited
to settings where the performative loss is essentially intractable, requiring
cost estimates of risk scores only in limited settings. We note ‘stability’ is
not necessarily desirable in terms of distribution of interventions: in the
QRISK3 setting, if an individual is at untreated risk of 50% and treated risk
of 10%, with treatment distributed proportionally to assessed risk, a ‘stable’
risk score would assess risk as e.g. 30%, prompting a mild intervention after
which true risk remained at 30%, regardless of treatment cost.

We found no literature directly addressing the focus in this paper:
determining how large a holdout set should be. Similar problems do arise in
clinical trial design: Stallard et al. (2017) estimate the optimal size of
clinical trial groups for a rare disease in which individuals not in the trial
stand to gain more than those in it, using a Bayesian decision-theoretic
approach accounting for benefit to future patients in the population.

OHS estimation requires quantification of expected material costs when using
risk scores trained to data of various sizes. Such costs will typically depend
on the error of risk predictions. The relation of predictive error to training
set size is known as the ‘learning curve’, which can sometimes be accurately
parametrised (Amari, 1993). A recent review paper suggests a power-law is
accurate for simple models (Viering and Loog, 2021).

###  3 Motivating example

In this example, we consider the ASPRE score (Akolekar et al., 2013) for
evaluating risk of pre-eclampsia (PRE), a hypertensive complication of
pregnancy, on the basis of predictors derived from ultrasound scans in early
pregnancy. Although treatable, PRE confers a serious risk to both the fetus
and the mother. The risk of PRE is lowered by treatment with aspirin through
the second and third trimesters (Rolnik et al., 2017b), but aspirin therapy
itself confers a slight risk, contraindicating universal treatment, and
suggesting prescription of aspirin only if the risk of PRE is sufficiently
high or other indications are present (ACOG, 2016). The ASPRE score was
developed to aid clinicians in estimating PRE risk and has been shown to be
useful in prioritising patients for aspirin therapy (Rolnik et al., 2017a). We
will not differentiate early- and late-stage PRE.

Due to changing population demographics, the influence of risk factors is
likely to change over time, and the ASPRE score will need to be periodically
updated. As discussed, a naïve re-fitting of a risk score on the basis of
($X:$) maternal assessment in early pregnancy and ($Y:$) eventual PRE
incidence could lead to dangerous underestimation of PRE risk, due to
individuals previously assessed as high risk being treated in response to the
assessment.

This could be avoided by maintaining a held-out set. For patients in this set,
no ASPRE score would be calculated at first scan, and treatment would be
according to best practice in the absence of a risk score. An updated ASPRE
score can then be fitted to data from these patients. Patients in this holdout
set go without the benefit of the ASPRE score, leading to a less accurate
allocation of prophylactic treatment (aspirin) and consequently a higher risk
of PRE (Rolnik et al., 2017a). However, an inappropriately small hold-out set
would lead to an inaccurate updated model, reducing the benefit of future use
of the score.

Our first question is whether a hold-out set is worth the cost, as opposed to
simply continuing to use the existing score or updating naively. We will show
that, as long as drift and intervention effects occur, the cost is generally
justified. A second question is, if implemented, how large should a holdout
set be? We will address this question in Section 6.3.

###  4 Theory

####  4.1 General setup

Our general strategy for safe model updating using a hold-out set is
illustrated as a causal graph in Figure 1. Each of the three columns is called
an ‘epoch’ (0,1,2 in subscripts), representing a period of time in which a
risk score is deployed and data gathered to update the risk score. Epochs
correspond to consecutive intervals of continuous time, where epoch 0
corresponds to interval $[0,1)$, epoch 1 to $[1,2)$ and so on. Ellipses
containing $X$ or $Y$ are covariates and outcomes (respectively) of
populations of samples. We use the shorthand $\\{X_{e}^{h}\\}$ and
$\\{X_{e}^{i}\\}$ (holdout and intervention respectively) to mean sets of
samples from a random process $X_{t}$, with $t\in[e,e+1)$ and
$X_{t}\sim\mu_{t}$, where $\mu_{t}$ may change with $t$ but has constant
support (e.g. $\mathds{R}^{p}$), and $\\{Y_{e}^{h}\\}$ and $\\{Y_{e}^{i}\\}$
to be corresponding set of observations of random processes $Y^{h}_{t}$,
$Y^{i}_{t}$ representing outcomes. We leave the distributions of
$Y^{h}_{t}|X_{t}$ and $Y^{i}_{t}|X_{t}$ unspecified for the moment. We presume
that the times $t$ corresponding to $\\{X_{e}^{i}\\}\cup\\{X_{e}^{h}\\}$ are
uniformly distributed on $[e,e+1)$ but will not presume this for
$\\{X_{e}^{i}\\}$ or $\\{X_{e}^{h}\\}$ individually.

Under a ‘native’ setting prior to deployment of a risk score, $\\{X_{0}\\}$
and $\\{Y_{0}\\}$ have a single causal link, modelled by risk score $\rho_{0}$
(leftmost epoch). Once $\rho_{0}$ is in use in the intervention set in epoch 1
(ellipses $\\{X_{1}^{i}\\}$, $\\{Y_{1}^{i}\\}$), a second causal pathway
through $\rho_{0}$ is established from $\\{X_{1}^{i}\\}$ to $\\{Y_{1}^{i}\\}$,
but there remains only one causal pathway from $\\{X_{1}^{h}\\}$ to
$\\{Y_{1}^{h}\\}$ in the holdout set (middle epoch). The updating process can
be continued rightwards ($\rho_{1}$, $\rho_{2}$, $\dots$).

Figure 1: Dynamics of a risk model across epochs when trained on a holdout
set.

To demonstrate why we opt to use a holdout set approach, we will define two
functions dependent on $t$. We define $f_{t}(x)$ as the probability of event
$Y_{t}=1$ given $X_{t}=x$ if no risk score is in place, and $g_{t}(x;\rho)$ as
the probability of event $Y_{t}=1$ given $X_{t}=x$ when a risk score $\rho$ is
used. In the holdout set schema in figure 1, we have
$P(Y_{t}^{h}=1)=f_{t}(X_{t}^{h})$ and
$P(Y_{t}^{i}=1)=g_{t}(X_{t}^{i},\rho_{e-1})$ for $t\in[e,e+1)$.

We generally wish to estimate $f_{t}$ rather than $g_{t}$, since $f_{t}(x)$
gives a risk of the event in question under standard practice. An agent may
opt to intervene in addition to standard practice if $f_{t}(x)$ is high. The
presence of drift is tantamount to an assumption that $f_{t}$ changes with
$t$, and indeed we will assume that it changes Lipschitz-continuously. Change
in $f_{t}$ with $t$ gives a motivation for updating a risk score at all, and
continuous change underlies the assumption that a risk score will still be
useful for a period of time after it is fitted. Drift will occur in many
realistic settings (Tsymbal, 2004). The presence of intervention effects
asserts that $f_{t}$ and $g_{t}$ differ; or rather, that some action may be
taken in response to the risk score. Since there is little reason to use a
risk score otherwise, we may assume intervention effects are typically
present.

Suppose we have a risk score $\rho_{0}$ fitted at $t\approx 0$ which
approximates $f_{0}$. We consider several options for what to do at $t=1$:
firstly, continue using our score $\rho_{0}$ (no update); secondly, refit a
risk score $\rho_{1}^{a}$ to observations $X_{1},Y_{1}$ with
$X_{1}\sim\mu_{1}$ and $Y_{1}|X_{1}\sim g_{t}(X_{1};\rho_{0})$ and use the new
risk score at $t=1$ (naive update) or thirdly, between $t=0$ and $t=1$, leave
out a ‘holdout’ set of samples and refit a new risk score $\rho_{1}^{b}$ to
these samples (holdout-update). If using the no-update option, performance
would be poor as $\rho_{0}\approx f_{0}\not\approx f_{1}$ (due to drift). If
using the naive-update option, performance would be poor as
$\rho_{1}^{a}\approx g_{1}(\cdot;\rho_{0})\not\approx f_{1}$ (due to
intervention). The holdout-update option can allow $\rho_{1}^{b}$ to
arbitrarily closely approximate $f_{1}$, by choosing the holdout set at times
in $[1-\epsilon,1)$ as $\epsilon\to 0$ (given continuous change in $f_{t}$).
Because of this, we can generally choose a finite holdout set size
$|\\{X_{1}^{h}\\}|$ above which the error in estimating $f_{2}$ will be less
than even the asymptotic error of the no-update or naive-update options. We
state and prove this observation formally in Supplement S2. We will briefly
demonstrate this by simulation in section 6.1.

####  4.2 Cost specification

We consider the aggregate sets of samples
$\left(\\{X_{e}^{h}\\},\\{Y_{e}^{h}\\}\right)\cup\left(\\{X_{e+1}^{i}\\},\\{Y_{e+1}^{i}\\}\right)$,
and denote the total number of samples as $N$. We denote the $n$ samples in
the ‘holdout’ set $\left(\\{X_{e}^{h}\\},\\{Y_{e}^{h}\\}\right)$ as $D_{n}$
(where ‘$D$’ indicates ‘data’). Since we are choosing
$\left(\\{X_{e}^{h}\\},\\{Y_{e}^{h}\\}\right)$ to be as close in time to
$\left(\\{X_{e+1}^{i}\\},\\{Y_{e+1}^{i}\\}\right)$ as possible, we will
presume that we have

$\displaystyle X_{e}^{h}$ $\displaystyle\sim\mu_{e+1}\hskip 20.0pt$
$\displaystyle\mathbb{E}(Y_{e}^{h}|X_{e}^{h})$
$\displaystyle=f_{e+1}(X_{e}^{h})$ $\displaystyle X_{e+1}^{i}$
$\displaystyle\sim\mu_{e+1}\hskip 20.0pt$
$\displaystyle\mathbb{E}(Y_{e+1}^{h}|X_{e+1}^{h})$
$\displaystyle=g_{e+1}(X_{e+1}^{h})$

A risk score is fitted to $D_{n}$ which approximates $f_{e+1}(x)$ and is used
in the intervention set $\left(\\{X_{e+1}^{i}\\},\\{Y_{e+1}^{i}\\}\right)$. We
will presume that samples in $D_{n}$ are pairwise independent, as are samples
in $\left(\\{X_{e+1}^{i}\\},\\{Y_{e+1}^{i}\\}\right)$, although samples in the
latter depend on $D_{n}$ through the fitted risk score.

We define $C_{1}(X)$ and $C_{2}(X;D_{n})$ as random variables associated with
the total ‘cost’ of an observation with covariates $X$ in the holdout set in
epoch $e$ and intervention set in epoch $e+1$ respectively, covering both the
cost of potentially managing the event $Y=1$ and any costs of intervention.
Function $C_{2}(X;D_{n})$ depends on $D_{n}$ only through the risk score
fitted to $D_{n}$. We define the expected cost per observation in the holdout
and intervention sets, respectively:

$k_{1}=\mathbb{E}_{X\sim\mu_{e+1},C_{1}}\left\\{C_{1}(X)\right\\}\hskip
20.0ptk_{2}(n)=\mathbb{E}_{X\sim\mu_{e+1},C_{2}}\left[\mathbb{E}_{D_{n}}\left\\{C_{2}(X;D_{n})\right\\}\right]$
(1)

Subscripted values $C_{1}$ and $C_{2}$ indicate variance in $C_{1}(X)$,
$C_{2}(X;D_{n})$ independent of $X,D_{n}$. We now express total cost $\ell$
across the sample population as

$\ell(n)=\underbrace{k_{1}n}_{\text{Holdout
set}}+\underbrace{k_{2}(n)(N-n)}_{\text{Intervention set}}$ (2)

The meaning of the $\ell$ is contextual; for instance, in QRISK3 it may mean
total number of deaths for a fixed healthcare budget.

####  4.3 Sufficient conditions for OHS

In this section, we consider conditions under which the cost $\ell(n)$ can be
readily optimised. We discuss estimation of $N$, $k_{1}$ and $k_{2}(\cdot)$ in
section 5. We begin with the following assumptions:

######  Assumption 1.

$k_{1}$ does not depend on $n$: in a medical context, treatment plans and
outcomes for patients without risk scores do not depend on the number of such
patients.

######  Assumption 2.

$k_{2}(n)$ is monotonically decreasing in $n$: the more data available to
train the risk score, the greater its clinical utility.

######  Assumption 3.

There exists $M\in(0,N)$ such that $n\geq M\Leftrightarrow k_{2}(n)\leq
k_{1}$: a good enough risk score will lead to better patient outcomes than
baseline treatment, and a poor enough risk score fitted to small amounts of
data leads to worse expected outcomes than baseline treatment.

######  Assumption 4.

$\mathbb{E}\left\\{k_{2}(i+1)-k_{2}(i)\right\\}>\mathbb{E}\left\\{k_{2}(j+1)-k_{2}(j)\right\\}$
for $1\leq i<j\leq N-1$, with expectations over training data: the ‘learning
curve’ for our risk score is convex; there are diminishing returns in the cost
per patient from adding more samples to the training data.

We may extend the domain of $k_{2}(\cdot)$, $\ell(\cdot)$ to the real interval
$[0,N)$ such that both functions are smooth; and $k_{2}^{\prime}(n)<0$, given
assumption 2, $k_{2}^{\prime\prime}(n)>0$, given assumption 4. This leads to
the result that there exists an optimal size for the holdout set minimizing
the expected total cost. The proof is given in Supplement S3.

######  Theorem 1.

Suppose assumptions 1-4 hold. Then there exists an OHS
$N_{*}\in\\{1,\dots,N-1\\}$ with $N\in\mathbb{N}$, such that:
$\ell(i)\geq\ell(j)\textrm{ for }0<i<j<N_{*}$ and $\ell(i)\leq\ell(j)\textrm{
for }N_{*}<i<j<N$

As an immediate consequence, we note that the OHS always exceeds the minimal
training sample size required to match baseline treatment.

######  Corollary 1.

The value of $N_{*}$ always exceeds the value of $M$ in assumption 3, since if
$N^{\prime}<M$ we have
$\ell(N^{\prime})=k_{1}N^{\prime}+k_{2}(N^{\prime})(N-N^{\prime})>k_{1}N^{\prime}+k_{1}(N-N^{\prime})=k_{1}N\geq\ell(N_{*})$

Consequently assumption 4 may be relaxed for $i,j<M$; we need only be
concerned with the behaviour of $k_{2}(n)$ at realistically large values of
$n$, rather than $n\in{1,\dots,M}$. We also have

$\lim_{N\to\infty}\lim_{n\to N}\ell^{\prime}(n)=\lim_{N\to\infty}\lim_{n\to
N}\left\\{k_{1}-k_{2}(n)+(N-n)k_{2}^{\prime}(n)\right\\}=k_{1}-\lim_{n\to\infty}k_{2}(n)$

and since $k_{1}>k_{2}(n)>0$ for large $n$, we have that expected total costs
$\ell(n)$ are increasing, but bounded by the per observation expected cost of
baseline treatment $k_{1}$.

####  4.4 Robustness to assumptions

The applicability of assumptions 1-4 in real world settings requires careful
consideration. We address violations of assumptions 2 and 4 in section 5.3 and
assumptions 1 and 4 here.

Assumption 1 is fundamental to the success of the holdout set concept. It may
be violated if, for instance, agents who can make interventions learn the
behaviour of a risk score and apply this to samples with no score. Such
violations are not of serious concern: if we presume that such changes in
agents endure over time, then they can be considered as simply contributing to
drift, which need not be independent of holdout set size (Supplement S2).

If ethically appropriate, assumption 1 could be assured by partitioning agents
to manage only samples in holdout sets or only in intervention sets (e.g.,
cluster randomisation). This requires assuming that changes in agent behaviour
as above _do not_ endure until the following epoch.

If assumption 3 fails due to $k_{2}(0)\leq k_{1}$, we may show a weaker result
(the presence of a non-trivial but potentially non-unique minimum loss) by
replacing assumptions 2, 3 and 4 with:

######  Assumption 5.

There exists an $0<M<N$ such that
$\frac{N-M}{N}(k_{1}-k_{2}(M))>k_{1}-k_{2}(0)$

This assumption is essentially stating that at some point the risk score will
greatly outperform a risk score built with no data. This leads to the result
that:

######  Theorem 2.

Suppose assumptions 1 and 5 hold, and $k_{2}(0)\leq k_{1}$. Then there exists
an $N_{*}\in\\{1,\dots,N-1\\}$ such that: $\ell(i)\geq\ell(N_{*})\textrm{ for
}i\in\\{1,\dots,N-1\\}$ and $\ell(i)>\ell(N_{*})\textrm{ for }i\in\\{0,N\\}$

In a setting in which $k_{1}>k_{2}(0)$ but one or more of assumptions 2-4 do
not hold, we have

######  Theorem 3.

Suppose assumption 1 holds, $k_{1}<k_{2}(0)$ and there exists $0<M<N$ such
that $k_{2}(M)<k_{1}$. Then there exists an $N_{*}\in\\{1,\dots,N-1\\}$ such
that: $\ell(i)\geq\ell(N_{*})\textrm{ for }i\in\\{1,\dots,N-1\\}$ and
$\ell(i)>\ell(N_{*})\textrm{ for }i\in\\{0,N\\}$.

Both results are proved in Supplement S3

In the setting where $k_{1}=k_{2}(0)$ and assumption 1 does not hold, it may
be reasonable in some settings to assign samples in the holdout set risk
scores based on no data (for example risk scores generated entirely from
expert opinion) and blind agents to holdout/intervention status. Under this
setting we may have greater assurance in assumption 1.

###  5 Estimation of OHS

####  5.1 Estimation of $k_{2}(n)$

We are aiming to find a holdout set size which minimises costs during epoch
$e\geq 1$ ($t\in[e,e+1)$) (noting that the holdout set will be used late in
the epoch when $t\approx e+1$) and must do this during epoch $e-1$ (when
$t<e$). We have the following:

  1. 1.

An approximate number of samples on which the model will be used or refitted

  2. 2.

A cohort of samples $(X,Y)$ with $X\sim\mu_{e-\epsilon}\approx\mu_{e}$,
$Y|X\sim f_{e-\epsilon}(X)\approx f_{e}(X)$, with $\epsilon$ small

where 2 is from a holdout set if $e>1$ or from initial training data if
$e-1=0$. We aim to estimate the cost function $\ell(n)$ at $t\in[e,e+1)$. Our
approach is to estimate $\ell(n)$ for $t=e-\epsilon$ and assume that the OHS
is approximately conserved from $t=e-\epsilon$ to $t=e+1$, though drift may
occur in $\ell(n)$.

At time $t=e-\epsilon$, we now need to estimate constants $N$, $k_{1}$, and
the function $k_{2}$. The constants $N$ and $k_{1}$ are straightforward: $N$,
the total number of samples on which a predictive score can be fitted or used,
will usually be known or specified (item 1); and $k_{1}$, the average cost per
sample under baseline behaviour without a score, can be estimated from
observed costs in the cohort in item 2. The function $k_{2}(\cdot)$ (equation
1) is more difficult to estimate, as it concerns quantifying costs of
hypothetical risk scores. We may tractably estimate $k_{2}(\cdot)$ by assuming
that

$\mathbb{E}_{X\sim\mu_{e},C_{2}}\left\\{C_{2}(X;D_{n})\right\\}=\mathscr{L}\\{\textrm{err}(\rho_{D_{n}})\\}$
(3)

where $\rho_{D_{n}}$ is a risk score fitted to $D_{n}$, $\textrm{err}(\cdot)$
is a measure of error, and $\mathscr{L}$ is some nondecreasing function. We
claim that in general circumstances we may take $\mathscr{L}(\cdot)$ to be
linear and $\textrm{err}(\cdot)$ to be expected mean-squared error (MSE) or a
similar general loss over $X\sim\mu_{t}$. We derive this heuristically in
Supplement S4 and derive expressions for $k_{2}(n)$ directly in a specific
case in section 6.3. Once $\mathscr{L}$ is known, this allows $k_{2}(n)$ to be
estimated readily by establishing the ‘learning curve’ of the risk score using
item 2.

Some direct estimates of $k_{2}(n)$ are necessary to determine $\mathscr{L}$.
A direct option is to designate subcohorts of the intervention set in epoch
$e-1$ to receive risk scores fitted to smaller subsamples of available
training data, allowing direct observation of the costs of such risk scores.
While simple, this approach is potentially ethically tenuous and expensive.
Simpler options include estimating the function $\mathscr{L}$ through expert
opinion or other outside information.

In summary, we recommend that $k_{2}(n)$ is estimated by jointly making a
small number of estimates during epoch $e-1$, either directly or indirectly,
to establish $\mathscr{L}$, and thereafter estimated by evaluating the error
of a risk score fitted to $n$ samples using the set in item 2 and transforming
it according to estimated $\mathscr{L}$.

####  5.2 Parametric estimation of OHS

A natural algorithm for estimating the OHS is immediately suggested by Theorem
1: assume $k_{2}$ is known up to parameters $\theta$, and estimate $N$,
$k_{1}$ and $\theta$ to estimate the OHS. Parameters $\theta$ of $k_{2}$ may
be estimated from observations of pairs $\\{n,k_{2}(n)\\}$, potentially with
error in $k_{2}(n)$. To minimize the number of estimates of $k_{2}(n)$ we
iteratively add observations $(n,k_{2}(n))$ to an existing set of observations
so as to greedily reduce expected error in the resultant OHS estimate.

We suggest a routine parametric algorithm (algorithm 1) with estimation of
asymptotic confidence intervals. Full details of theory, proofs and algorithms
are given in Supplement S5.

1 $\mathbf{n},\mathbf{k_{2}},\boldsymbol{\sigma}^{2}\leftarrow$ some initial
values $\mathbf{n}$ with $(\boldsymbol{k_{2}})_{i}\approx k_{2}(n_{i})$,
$(\boldsymbol{\sigma}^{2})_{i}=\textrm{var}(\hat{k_{2}}(n_{i}))$;

2 while  _$|\mathbf{n}| <$ total iterations_ do

3 Find best new value $\tilde{n}$ to add to $\mathbf{n}$ ;

4 Estimate $\hat{k_{2}}(\tilde{n})\approx k_{2}(\tilde{n})$ ;

5 $\mathbf{n}\leftarrow(\mathbf{n}\cup\tilde{n})$,
$\mathbf{k_{2}}\leftarrow\left\\{\mathbf{k_{2}}\cup\hat{k_{2}}(\tilde{n})\right\\}$,
$\boldsymbol{\sigma}^{2}\leftarrow\boldsymbol{\sigma}^{2}\cup\textrm{var}\left\\{\hat{k_{2}}(\tilde{n})\right\\}$
;

6

7 end while

8return Re-estimate OHS $n_{*}^{final}$ from
$\mathbf{n},\mathbf{k_{2}},\boldsymbol{\sigma}$ ;

Algorithm 1 Parametric OHS estimation overview

Consistency of algorithm 3 depends on whether $\mathbf{n}$ eventually contains
enough elements of sufficient multiplicity to estimate $\Theta_{0}$
consistently. Sampling some positive proportion of values of $\mathbf{n}$
randomly from $\\{1,\dots,N\\}$ guarantees that the multiplicity of all
$n\in\mathbf{n}$ almost surely eventually exceeds any finite value, readily
ensuring consistency. Finite-sample bias of $n_{*}^{final}$ depends on
$\nabla_{\Theta}n_{*}$ and the variance of $\Theta$. See Supplementary Figure
8 for typical forms of $\nabla_{\Theta}n_{*}$.

####  5.3 Semi-parametric (emulation) estimation of OHS

Parametrization of $k_{2}(n)$ may be inappropriate if the learning curve of
the risk score or the relation between the learning curve and $k_{2}(n)$ (from
section 5.1) are complex (Viering and Loog, 2021). We propose a second
algorithm which is less reliant on assuming a parametric form for $k_{2}(n)$,
using Bayesian optimisation (Brochu et al., 2010). We quantify the uncertainty
of the form for $k_{2}(n)$ through the construction of a Gaussian process
emulator of $\ell$ . The prior mean function for this emulator takes a
particular parametric form, but crucially can deviate from this prior function
with the addition of data.

We take the minimum of its posterior mean over $n$ to be our OHS estimate,
efficiently choosing values of $n$ at which to estimate $\ell(n)$ using an
‘expected improvement’ function $EI(\cdot)$, which also provides a natural
stopping criterion; if $EI(n)>\tau$ we roughly expect the minimum cost to
decrease by at least $\tau$ from adding another estimate of $\ell(n)$ to our
data. Our approximate procedure is given in algorithm 2. Detailed algorithm
and proofs of consistency are in Supplement S6.

1 $\mathbf{n},\mathbf{d},\boldsymbol{\sigma}^{2}\leftarrow$ some initial
values $\mathbf{n}$ with $(\boldsymbol{d})_{i}=d(n_{i})\approx\ell(n)$,
$(\boldsymbol{\sigma}^{2})_{i}=\textrm{var}\\{d(n_{i})\\}$ ;

2 Estimate mean and variance of Gaussian process $\ell(n)$ and function
$EI(n)$;

3 while  _$\max_{n\in\\{1,\dots,N\\}}\\{EI(n)\\} >\tau$_ do

4 $\tilde{n}\leftarrow\arg\max_{n\in\\{1,\dots,N\\}}EI(n)$ ;

5 Estimate $d(\tilde{n})\approx k_{2}(\tilde{n})$ ;

6 $\mathbf{n}\leftarrow(\mathbf{n}\cup\tilde{n})$;
$\mathbf{d}\leftarrow(\mathbf{d}\cup d\\{\tilde{n})\\}$;
$\boldsymbol{\sigma}^{2}\leftarrow\left[\boldsymbol{\sigma}^{2}\cup\textrm{var}\\{d(\tilde{n})\\}\right]$
;

7 Re-estimate mean and variance of Gaussian process $\ell(n)$ and function
$EI(n)$ ;

8

9 end while

return $n_{*}^{final}=\arg\min_{n\in\\{1,\dots,N\\}}\left\\{\textrm{Mean of
$\ell(n)$}\right\\}$

Algorithm 2 Emulation OHS estimation; minimum cost improvement $\tau$

###  6 Simulations

####  6.1 Simulation showing dominance of holdout set approach

We briefly demonstrate the theory described in section 4.1 using simulated
data similar to our motivating example (details in Supplement S2.1). Figure 2
shows total cost of no-update (‘none’), naive-update (‘naive’) and holdout-
update (‘H.S’) at two holdout sizes over time. As drift occurs in $f_{t}$, the
total costs of the holdout-set approaches are lowest. Choice of holdout set
size aims to balance increased costs due to non-intervention in the holdout
set (the ‘spikes’) against inaccuracy in fitted scores after drift. We
demonstrate the natural emergence of an OHS in a simulated context in
Supplement S7.

Figure 2: Loss over time of no-update, naive-update and holdout-update
strategies.

####  6.2 Comparison of parametric and emulation algorithms

In this section, we give circumstances in which one of algorithms 1, 2 may be
preferable to the other. We consider two versions of the function $k_{2}(n)$:

$k_{2}^{p}(n)=an^{-b}+c\hskip 10.0pt\mbox{, }\hskip
10.0ptk_{2}^{np}(n)=an^{-b}+c+\frac{10^{4}}{\sqrt{2\pi}}\exp\left(-\frac{1}{2}\left(\frac{n-4\times
10^{4}}{8\times 10^{3}}\right)^{2}\right)$

where: ‘p’/‘np’: ‘parametric assumptions satisfied/not satisfied’, and
$\theta=(a,b,c)=(10000,1.2,0.2)$. We assume $N$ and $k_{1}$ are known to be
$1\times 10^{5}$ and $0.4$ respectively. For emulation, we use a kernel width
$\zeta=5000$ and variance $\sigma_{u}^{2}$ of $1\times 10^{7}$.

The function $k_{2}^{np}(n)$ exhibits ‘double-descent’ behaviour
(Supplementary Figures 11(a), 11(b)), which is possible for learning curves
(Viering and Loog, 2021). This violates assumptions 2, 4.

We firstly show the distribution of estimates of OHS using our two algorithms
when $k_{2}$ takes either form above. To fit $k_{2}$, we use $\mathbf{n}$ as
200 values of $n$ randomly chosen from $\\{1,\dots,N\\}$, with values
$\mathbf{k_{2}}$ independently sampled as $(\mathbf{k_{2}})_{i}\sim
N\\{k_{2}(n_{i}),\sigma_{i}^{2}\\}$, where $\boldsymbol{\sigma}\sim
U(0.001,0.02)$. Supplementary Figure 11(c) shows the distributions and medians
of OHS estimates using the parametric and emulation algorithm in settings with
parametric assumptions either satisfied or unsatisfied.

The parametric OHS estimate is empirically unbiased and has less variance than
the emulation estimate when parametric assumptions are satisfied, but is
biased when they are not. Variance of OHS estimates using the emulation method
is lower when parametric assumptions are not satisfied because the true cost
function has a sharper minimum in that case (see Supplementary Figure 11(b)).
Because the cost function is ‘flat’ around the minimum in the setting where
parametric assumptions are satisfied (Supplementary Figure 11(b)), the
consequences of the high variance of the semi-parametric (emulation) estimator
are minimal, as the cost is similar across a range of values near the OHS.

Figure 4: Convergence rates with parametric (top) and emulation (bottom)
algorithms, using either a random (black) or greedy (gray) method to select
the next value of $n$ to add to $\mathbf{n}$, with parametric assumptions
satisfied (left) or unsatisfied (right). Simulations run for 200 datasets from
underlying model. In larger panels, horizontal lines show true optimal holdout
set (OHS) size; vertical lines indicate when $\geq 2.5\%$ of runs return that
OHS within a grid of size 1000. Smaller panels show root mean-square error
between total costs from simulations and minimal total cost. Note variable
axis scaling on left and right.

We next examine the consequences of sampling the ‘next’ value of $n$
$\tilde{n}$ greedily, using equation 31 as in algorithm 3 or $EI$ as in
algorithm 4, rather than randomly with $\tilde{n}$ chosen uniformly in
$\\{1,\dots,N\\}$, by comparing the rates of convergence of OHS estimates.
Figure 4 show medians and a rough discrete kernel estimate of OHS estimates at
various sizes of $|\mathbf{n}|$, where $\mathbf{n}$ is generated by adding
points $\tilde{n}$ either randomly or systematically, after choosing the
initial five values in $\mathbf{n}$ randomly.

Convergence is faster when ‘next points’ are picked systematically rather than
randomly and when using parametric estimates (though these are biased and
inconsistent for $k_{2}^{np}$). This is highlighted by the smaller panels
which show the root mean-square error between the total cost at the estimated
optimal sizes and the total cost at the true OHS. In particular, that the non-
parametric method shows bifurcation detecting both local minima whilst the
parametric method converges to a mid-point which is far from optimal in terms
of total costs.

####  6.3 Illustration in realistic setting

We now return to our motivating example. Supposing ASPRE is to be refitted
every five years, the intervention set should include all individuals in the
subsequent years before the model is refitted, and all individuals not used in
the next refitting procedure. Suppose we are refitting ASPRE to use in a
population of 5 million individuals, from which we have approximately 80,000
new pregnancies per year. Thus $N\approx 400,000$ (SE: 1500). See Supplement
S8 for further details.

We presume a simple clinical action in which a fixed proportion $\pi=10\%$ of
individuals at the highest assessed PRE risk are treated with aspirin. We
assume that if untreated with aspirin, a proportion $\pi_{0}$ of individuals
designated to be ‘low-risk’ (lowest 90%) will develop PRE, as will a
proportion $\pi_{1}$ of individuals designated high-risk. Under current pre-
ASPRE best-practice guidelines we have $\pi_{0}\approx 0.02$ (SE $0.0009$) and
$\pi_{1}\approx 0.08$ (SE $0.008$) (see Supplement S8). Aspirin reduces PRE
risk by approximately $1-\alpha=63\%$ (SE 0.09) (Rolnik et al., 2017b). We
denote ‘cost’ as simply the number of cases of PRE in a population, so total
expected cost per individual under ‘baseline’ treatment (clinical actions
without the aid of a risk model) is

$k_{1}=\pi_{0}(1-\pi)+\pi_{1}\pi\alpha\approx 0.022$ (4)

with standard error approximately 0.001. Note that this is not equal to the
untreated PRE risk in the population, since some proportion of individuals are
treated pre-emptively.

The data used to fit the initial ASPRE model can be used to estimate
$k_{2}(n)$ for potential model updates. We do not have access to this dataset,
but demonstrate estimation of a learning curve on synthetic data designed to
resemble it. Although $k_{2}(n)$ is easy and fast to estimate here, in order
to mimic a real example where such estimation is time consuming or costly we
restrict ourselves to use only $|\mathbf{n}|=120$ values of $n$, determined
using either algorithm 3 or 4. For both algorithms, we assumed a power-law
form $k_{2}\\{n;\theta=(a,b,c)\\}=an^{-b}+c$.

Using the parametric algorithm, we found an OHS of 10271 (90% CI 8103-12438),
with minimum cost (expected cases over five years) of 8172. Using the
emulation algorithm, we found an OHS of 13313 with an expected cost of 8164,
with holdout sizes of 9210-17619 having a probability $>0.1$ of cost $<8164$.
Figure 5 shows estimated cost functions, OHSs, and error using the two
algorithms.

(a)

(b)

(c)

Figure 5:  Estimation of cost functions (black lines), OHS (red dots) and
error using parametric (middle) and emulation (left; blue lines
$\mu(n)+3\sqrt{\Psi(n)}$) algorithms, and change in estimated OHS and error
with number of sample points $|\mathbf{n}|$ (right; red lines: parametric,
blue/dashed lines: emulation). Note that the ‘best’ points (black dots) to
optimize parametric estimation are spread-out to estimate $\theta$ well, but
for emulation they are clustered for accurate locally approximation. Error
measures for OHS in parametric and emulation algorithms (red/blue shaded
respectively) have different meanings and are not comparable.

###  7 Concluding remarks

In this work we propose the use of a holdout set to safely update predictive
models, and describe considerations in determining the optimal size of such a
set. We establish theoretical properties of this optimal size under reasonable
assumptions, and develop two algorithms for estimating it, evaluating their
use in both a toy simulation and a real-life motivated simulation. The holdout
set approach comprises a practical and simple approach to an important problem
in practical machine learning, which will become particularly important as
risk scores come to be used to prompt intervention in real-world applications.

Other solutions to the model updating problem have been proposed, such as
developing models that introduce missing causal connections (Alaa and van der
Schaar, 2018; Sperrin et al., 2019) and using several predictive scores
together (Liley, 2021). Such solutions tend to be difficult to implement: more
comprehensive modelling generally requires some observation of the
intervention, requiring more data over time, and parallel use of several risk
scores is highly complex. With this in mind, holdout sets could prove valuable
in updating strategies since, in principle, they can be applied in any
setting. Availability of direct risk-score-naive data through the holdout set
can also facilitates post-deployment maintenance and surveillance (Sperrin et
al., 2019).

We strongly suggest planning an updating strategy for a risk model _before_ it
is deployed. This work illustrates one strategy in this direction and we hope
stimulates both use of and extensions of such methods for safe predictive
score updating. For the interested reader, see further discussion in
Supplement S9.

### Acknowledgements

SH’s contributions arose from an MSc dissertation for the MISCADA programme at
Durham University. LJMA was partially supported by a Health Programme
Fellowship at the Alan Turing Institute. JL and LJMA were partially supported
by Wave 1 of The UKRI Strategic Priorities Fund under the EPSRC Grant
EP/T001569/1, particularly the ‘Health’ theme within that grant and The Alan
Turing Institute; and by Health Data Research UK, an initiative funded by
UKRI, Department of Health and Social Care (England), the devolved
administrations, and leading medical research charities; SRE is funded by the
EPSRC doctoral training partnership at Durham University, grant reference
EP/R513039/1. We thank Catalina Vallejos, Sebastian Vollmer and Bilal Mateen
for helpful discussion.

### Code availability

We have implemented all functionality described in this manuscript in an R
package on CRAN (OptHoldoutSize) and at
https://github.com/jamesliley/OptHoldoutSize. Pipelines to generate all
results in this manuscript are available at
https://github.com/jamesliley/OptHoldoutSize_pipelines.

### References

  * ACOG (2016) Practice advisory on low-dose aspirin and prevention of preeclampsia: updated recommendations.  American College of Obstetricians and Gynecologists (ACOG).  Cited by: §3.
  * R. Akolekar, A. Syngelaki, L. Poon, D. Wright, and K. H. Nicolaides (2013) Competing risks model in early screening for preeclampsia by biophysical and biochemical markers.  Fetal diagnosis and therapy 33 (1),  pp. 8–15.  Cited by: §3.
  * A. M. Alaa and M. van der Schaar (2018) Autoprognosis: automated clinical prognostic modeling via Bayesian optimization with structured kernel learning.  arXiv preprint arXiv:1802.07207.  Cited by: §7.
  * S. Amari (1993) A universal theorem on learning curves.  Neural networks 6 (2),  pp. 161–166.  Cited by: §2.
  * E. Brochu, V. M. Cora, and N. De Freitas (2010) A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning.  arXiv preprint arXiv:1012.2599.  Cited by: §5.3.
  * J. Chalmers, M. Pullan, B. Fabri, J. McShane, M. Shaw, N. Mediratta, and M. Poullis (2013) Validation of EuroSCORE II in a modern cohort of patients undergoing cardiac surgery.  European Journal of Cardio-Thoracic Surgery 43 (4),  pp. 688–694.  Cited by: §2.
  * G. S. Collins, J. B. Reitsma, D. G. Altman, and K. G. Moons (2015) Transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD): the TRIPOD statement.  Circulation 131 (2),  pp. 211–219.  Cited by: §2.
  * J. A. Cook and G. S. Collins (2015) The rise of big clinical databases.  British Journal of Surgery 102 (2),  pp. e93–e101.  External Links: Document, Link Cited by: §2.
  * D. Drusvyatskiy and L. Xiao (2020) Stochastic optimization with decision-dependent distributions.  arXiv preprint arXiv:2011.11173.  Cited by: §2.
  * S. G. Finlayson, A. Subbaswamy, K. Singh, J. Bowers, A. Kupke, J. Zittrain, I. S. Kohane, and S. Saria (2020) The clinician and dataset shift in artificial intelligence.  The New England Journal of Medicine,  pp. 283–286.  Cited by: §2.
  * J. Hippisley-Cox, C. Coupland, and P. Brindle (2017) Development and validation of QRISK3 risk prediction algorithms to estimate future risk of cardiovascular disease: prospective cohort study.  BMJ 357.  Cited by: §1, §2.
  * Z. Izzo, J. Zou, and L. Ying (2021) How to learn when data gradually reacts to your model.  arXiv preprint arXiv:2112.07042.  Cited by: §2.
  * D. Kansagara, H. Englander, A. Salanitro, D. Kagen, C. Theobald, M. Freeman, and S. Kripalani (2011) Risk prediction models for hospital readmission: a systematic review.  Jama 306 (15),  pp. 1688–1698.  Cited by: §2.
  * M. C. Lenert, M. E. Matheny, and C. G. Walsh (2019) Prognostic models will be victims of their own success, unless….  Journal of the American Medical Informatics Association 26 (12),  pp. 1645–1650.  Cited by: §1.
  * Q. Li and H. Wai (2021) State dependent performative prediction with stochastic approximation.  arXiv preprint arXiv:2110.00800.  Cited by: §2.
  * J. Liley, G. Bohner, S. R. Emerson, B. A. Mateen, K. Borland, D. Carr, S. Heald, S. D. Oduro, J. Ireland, K. Moffat, R. Porteous, S. Riddell, S. Rogers, N. Cunningham, C. Holmes, K. Payne, S. J. Vollmer, C. A. Vallejos, and L. J. Aslett (2021a) Development and assessment of a machine learning tool for predicting emergency admission in Scotland.  medRxiv.  Cited by: §2.
  * J. Liley, S. R. Emerson, B. A. Mateen, C. A. Vallejos, L. J. Aslett, and S. J. Vollmer (2021b) Model updating after interventions paradoxically introduces bias.  AISTATS proceedings.  Cited by: §1, §1.
  * J. Liley (2021) Stacking interventions for equitable outcomes.  arXiv preprint arXiv:2110.04163.  Cited by: §2, §7.
  * J. Lu, A. Liu, F. Dong, F. Gu, J. Gama, and G. Zhang (2018) Learning under concept drift: a review.  IEEE Transactions on Knowledge and Data Engineering 31 (12),  pp. 2346–2363.  Cited by: §2.
  * C. Mendler-Dünner, J. C. Perdomo, T. Zrnic, and M. Hardt (2020) Stochastic optimization for performative prediction.  arXiv preprint arXiv:2006.06887.  Cited by: §2.
  * S. A. Nashef, F. Roques, L. D. Sharples, J. Nilsson, C. Smith, A. R. Goldstone, and U. Lockowandt (2012) EuroSCORE II.  European journal of cardio-thoracic surgery 41 (4),  pp. 734–745.  Cited by: §2.
  * J. Perdomo, T. Zrnic, C. Mendler-Dünner, and M. Hardt (2020) Performative prediction.  In International Conference on Machine Learning,  pp. 7599–7609.  Cited by: §1, §2.
  * D. L. Rolnik, D. Wright, L. Poon, A. Syngelaki, N. O’Gorman, C. de Paco Matallana, R. Akolekar, S. Cicero, D. Janga, M. Singh, et al. (2017a) ASPRE trial: performance of screening for preterm pre-eclampsia.  Ultrasound in obstetrics & gynecology 50 (4),  pp. 492–495.  Cited by: §3, §3.
  * D. L. Rolnik, D. Wright, L. C. Poon, N. O’Gorman, A. Syngelaki, C. de Paco Matallana, R. Akolekar, S. Cicero, D. Janga, M. Singh, et al. (2017b) Aspirin versus placebo in pregnancies at high risk for preterm preeclampsia.  New England Journal of Medicine 377 (7),  pp. 613–622.  Cited by: §3, §6.3.
  * D. M. Shahian, J. P. Jacobs, V. Badhwar, P. A. Kurlansky, A. P. Furnary, J. C. Cleveland Jr, K. W. Lobdell, C. Vassileva, M. C. W. von Ballmoos, V. H. Thourani, et al. (2018) The Society of Thoracic Surgeons 2018 adult cardiac surgery risk models: part 1—background, design considerations, and model development.  The Annals of thoracic surgery 105 (5),  pp. 1411–1418.  Cited by: §2.
  * M. Sperrin, D. Jenkins, G. P. Martin, and N. Peek (2019) Explicit causal reasoning is needed to prevent prognostic models being victims of their own success.  Journal of the American Medical Informatics Association 26 (12),  pp. 1675–1676.  Cited by: §2, §7.
  * N. Stallard, F. Miller, S. Day, S. W. Hee, J. Madan, S. Zohar, and M. Posch (2017) Determination of the optimal sample size for a clinical trial accounting for the population size.  Biometrical Journal 59 (4),  pp. 609–625.  Cited by: §2.
  * E. J. Topol (2019) High-performance medicine: the convergence of human and artificial intelligence.  Nature medicine 25 (1),  pp. 44–56.  Cited by: §1.
  * A. Tsymbal (2004) The problem of concept drift: definitions and related work.  Computer Science Department, Trinity College Dublin 106 (2),  pp. 58.  Cited by: §1, §4.1.
  * USFDA et al. (2019) Proposed regulatory framework for modifications to artificial intelligence/machine learning (ai/ml)-based software as a medical device (samd).  Cited by: §1.
  * T. Viering and M. Loog (2021) The shape of learning curves: a review.  arXiv preprint arXiv:2103.10948.  Cited by: §2, §5.3, §6.2.
  * E. Wallace, E. Stuart, N. Vaughan, K. Bennett, T. Fahey, and S. M. Smith (2014) Risk prediction models to predict emergency hospital admission in community-dwelling adults: a systematic review.  Medical care 52 (8),  pp. 751.  Cited by: §2.

## Supplementary materials

###  S1 General notation

Throughout this document we will take $X$ to generically mean ‘covariates’ and
$Y$ to mean ‘outcomes’. The subscript $t$ will be taken to mean ‘time’ in a
continuous sense, and the subscript $e$ to mean ‘epoch’, referring to
consecutive episodes of time. The superscript $h$ will correspond to the
holdout set, and $i$ to the intervention set. The number $N$ will refer to the
total number of samples on which a risk score may be trained or used during an
epoch, and $n$ to denote a holdout set size, usually taken to be variable. We
denote the standard normal PDF and CDF by $\phi(\cdot)$, $\Phi(\cdot)$
respectively.

###  S2 Formal statement of dominance of holdout-set based updating

As in section 4, we presume a random process $X_{t}$, $t\in\mathds{R}^{+}$,
where $\\{X_{e}^{h}\\}\cup\\{X_{e}^{i}\\}$ is a set of instances of $X_{t}$
for values of $t$ in $[e,e+1)$. We presume $X_{t}$ has distribution $\mu_{t}$
where $\mu_{t}$ has constant support $\mathscr{X}$. We will assume $\mu_{t}$
varies continuously, in a manner to be specified.

We consider two functions $f_{t}:\mathscr{X}\to[0,1]$ and
$g_{t}(\cdot;\rho):\mathscr{X}\to[0,1]$, which correspond to risk given a
value of $X_{t}$ without and with an existing risk score $\rho$ in place
respectively, as described in section 4. For our purposes in this section,
$\rho$ will always be $\rho_{0}$ (fitted at time $t=0$) so we will omit the
dependence of $g_{t}$ on $\rho$. To avoid causal considerations, we will
define two further Bernoulli random variables: $Y_{f}$, with
$P(Y^{f}_{t}=1|X_{f}=x)=f_{t}(x)$ and $Y_{g}$, with
$P(Y^{g}_{t}=1|X_{f}=x)=g_{t}(x)$.

We aim to estimate risk scores $\rho_{\cdot}:\mathscr{X}\to[0,1]$ which
estimate $f_{t}$ for $t$ on some interval $[e,e+1)$. We will, for the moment,
aim to estimate $f_{t}$ for $t\in[1,2)$. We presume that we already have:

  1. 1.

A risk score $\rho_{0}$ which approximates $f_{0}$ to arbitrary accuracy

  2. 2.

Arbitrarily many samples of $X_{t},Y^{g}_{t}$ with $t\in[1-\epsilon_{1},1)$

  3. 3.

A total of $m$ samples of $X_{t},Y^{f}_{t}$ with $t\in[1-\epsilon_{1},1)$
(holdout set)

As in the main manuscript, considering the use of a risk score in time period
$t\in[1,2)$, we will call the use of $\rho_{0}$ a ‘no-update’ strategy, the
use of a risk score $\rho_{(n)}$ fitted to samples 2 a ‘naive-update’ strategy
and the use of a risk score $\rho_{(h)}$ fitted to samples 3 a ‘holdout-
update’ strategy.

We denote the mean absolute error and mean square error of a risk score $\rho$
at time $t$ as

$\displaystyle\ell_{t}(\rho)$ $\displaystyle=\int|\rho(x)-f_{t}(x)|d\mu_{t}$
$\displaystyle\ell_{t}^{2}(\rho)$
$\displaystyle=\int(\rho(x)-f_{t}(x))^{2}d\mu_{t}$

The following theorem establishes that under certain conditions both measures
of error are lower for $\rho_{(h)}$ than for $\rho_{0}$ or $\rho_{(n)}$ when
$t\in[1,1+\delta)$ for some $\delta>0$.

######  Theorem S4.

Suppose that:

  1. 1.

$|f_{0}(x)-f_{1}(x)|\geq\gamma_{1}$ for all
$x\in\mathscr{X}_{1}\subseteq\mathscr{X}$ with
$P_{\mu_{1}}(\mathscr{X}_{1})=\kappa_{1}>0$ (non-negligible drift)

  2. 2.

$|g_{1-\epsilon}(x)-f_{1}(x)|\geq\gamma_{2}$ for all
$x\in\mathscr{X}_{2}\subseteq\mathscr{X}$ and $0\leq\epsilon<\epsilon_{1}$
with $P_{\mu_{1}}(\mathscr{X}_{2})=\kappa_{2}>0$ (non-negligible intervention)

  3. 3.

$f_{t}(x)$ is uniformly $\alpha$-Lipschitz continuous in $t$ (continuous drift
in $Y^{f}|X$)

  4. 4.

$\mu_{t}$ is $\alpha_{2}$-Lipschitz continuous in $t$ with respect to total
variational distance distance (continuous drift in $X$)

  5. 5.

If a risk score $\rho(x)$ is fitted to data $D(n)$ consisting of $n$
identically distributed samples from $X_{t},Y^{f}_{t}$ or $X_{t},Y^{g}_{t}$
then for any $\epsilon>0$ we have $\lim_{n\to\infty}P_{D(n)}\left(\sup_{x\in
X}|f_{t}(x)-\rho(x)|<\epsilon\right)=1$ (risk score consistency)

Suppose $\rho_{1}$, $\rho_{2}$, $\rho_{3}$ are fitted as above. Then given any
$\epsilon_{1},\epsilon_{2}>0$ and some fixed $\delta$ satisfying

$\delta<\min\left(\frac{\gamma_{1}\kappa_{1}}{2\alpha+\gamma_{1}\alpha_{2}},\frac{\gamma_{2}\kappa_{2}}{2\alpha+\gamma_{2}\alpha_{2}}\right)$
(5)

we can find an $m_{1}$ such that for $m\geq m_{1}$ we have

$P_{D(m)}\left(\ell_{t}(\rho_{(h)}<\min\left(\ell_{t}(\rho_{0}),\ell_{t}(\rho_{(n)})\right)\right)>1-\epsilon_{2}$
(6)

for $1\leq t\leq 1+\delta$. If $\delta$ satifies:

$\displaystyle\delta<\min$
$\displaystyle\left(\frac{\sqrt{\left(2\alpha^{2}+\alpha_{2}\gamma_{1}^{2}\right)^{2}+4\alpha_{2}\gamma_{1}^{2}\kappa_{1}}-\left(2\alpha^{2}+\alpha_{2}\gamma_{1}^{2}\right)}{2\alpha^{2}},\right.$
$\displaystyle\phantom{=}\left.\frac{\sqrt{\left(2\alpha^{2}+\alpha_{2}\gamma_{2}^{2}\right)^{2}+4\alpha_{2}\gamma_{2}^{2}\kappa_{2}}-\left(2\alpha^{2}+\alpha_{2}\gamma_{2}^{2}\right)}{2\alpha^{2}}\right)$
(7)

then we may find an $m_{2}$ such that for $m\geq m_{2}$

$P_{D(m)}\left(\ell_{t}^{2}(\rho_{(h)})<\min\left(\ell_{t}^{2}(\rho_{0}),\ell_{t}^{2}(\rho_{(n)})\right)\right)>1-\epsilon_{2}$
(8)

for $1\leq t\leq 1+\delta$.

###### Proof.

From $\alpha$-uniform Lipschitz continuity of $f_{t}$ (assumption 3) we have,
for all $x$, $\delta$:

$f_{t}(x)-\alpha\delta\leq f_{t+\delta}(x)\leq f_{t}(x)+\alpha\delta$ (9)

and from assumption 4 we have

$\left|\int_{\mathscr{X}_{0}}(d\mu_{t}-d\mu_{t+\delta})\right|\leq\alpha_{2}\delta$
(10)

From assumption 5 we have

$P\left(\sup_{x\in\mathscr{X}}|\rho_{(h)}(x)-f_{1}(x)|<\epsilon\right)>1-q(\epsilon,m)$
(11)

where $q(\epsilon,m)$ is a function which converges to 0 as $m\to\infty$ for
any $\epsilon>0$.

Now for $\delta<\frac{\gamma_{1}\kappa_{1}}{2\alpha+\gamma_{1}\alpha_{2}}$ we
have

$\displaystyle\ell_{1+\delta}(\rho_{0})$
$\displaystyle=\int|f_{0}(x)-f_{1+\delta}(x)|d\mu_{t+\delta}$
$\displaystyle\geq\int|f_{0}(x)-f_{1}(x)|d\mu_{t+\delta}-\int|f_{1+\delta}(x)-f_{1}(x)|d\mu_{t+\delta}$
$\displaystyle\geq\int_{\mathscr{X}_{1}}\inf_{x\in\mathscr{X}_{1}}|f_{0}(x)-f_{1}(x)|d\mu_{t+\delta}-\int\sup_{x\in\mathscr{X}}|f_{1+\delta}(x)-f_{1}(x)|d\mu_{t+\delta}$
$\displaystyle\geq\int_{\mathscr{X}_{0}}\gamma_{1}d\mu_{t+\delta}-\alpha\delta$
$\displaystyle\geq\gamma_{1}\left(\int_{\mathscr{X}_{0}}d\mu_{t}-\alpha_{2}\delta\right)-\alpha\delta$
$\displaystyle\geq\gamma_{1}(\kappa_{1}-\alpha_{2}\delta)-\alpha\delta$
$\displaystyle>\alpha\delta$ (12)

For $\delta$ satisfying inequality 7 we have

$\displaystyle\ell^{2}_{1+\delta}(\rho_{0})$
$\displaystyle=\int\left(f_{0}(x)-f_{1+\delta}(x)\right)^{2}d\mu_{t+\delta}$
$\displaystyle=\int(f_{0}(x)-f_{1}(x))^{2}d\mu_{t+\delta}+\int(f_{1+\delta}(x)-f_{1}(x))^{2}d\mu_{t+\delta}$
$\displaystyle\phantom{=}-\int(f_{1+\delta}(x)-f_{1}(x))(f_{0}(x)-f_{1}(x))d\mu_{t+\delta}$
$\displaystyle\geq\int_{\mathscr{X}_{0}}\inf_{x\in\mathscr{X}_{0}}|f_{0}(x)-f_{1}(x)|^{2}d\mu_{t+\delta}+\int\inf_{x\in\mathscr{X}}|f_{1+\delta}(x)-f_{1}(x)|^{2}d\mu_{t+\delta}$
$\displaystyle\phantom{=}-2\left(\int(f_{1+\delta}(x)-f_{1}(x))^{2}d\mu_{t+\delta}\int(f_{0}(x)-f_{1}(x))^{2}d\mu_{t+\delta}\right)^{1/2}$
$\displaystyle\geq\gamma_{1}^{2}\int_{\mathscr{X}_{0}}d\mu_{t+\delta}+0$
$\displaystyle\phantom{=}-2\left(\int\sup_{x\in\mathscr{X}}|f_{1+\delta}(x)-f_{1}(x)|^{2}d\mu_{t+\delta}\int\sup_{x\in\mathscr{X}}|f_{0}(x)-f_{1}(x)|^{2}d\mu_{t+\delta}\right)^{1/2}$
$\displaystyle\geq\gamma_{1}^{2}\left(\int_{\mathscr{X}_{0}}d\mu_{t}-\alpha_{2}\delta\right)-2\sqrt{\alpha^{2}\delta^{2}\cdot\alpha^{2}}$
$\displaystyle=\gamma_{1}^{2}(\kappa_{1}-\alpha_{2}\delta)-2\alpha^{2}\delta$
$\displaystyle>\alpha^{2}\delta^{2}$ (13)

Similarly for
$\delta<\frac{\gamma_{2}\kappa_{2}}{2\alpha+\gamma_{1}\alpha_{2}}$ and
$0<\epsilon\leq\epsilon_{1}$ we have

$\ell_{2+\delta}(\rho_{2})\geq\gamma_{2}(\kappa_{2}-\alpha_{2}\delta)-\alpha\delta>\alpha\delta$
(14)

and for $\delta$ satisfying inequality 7 we have

$\displaystyle\ell^{2}_{1+\delta}(\rho_{2})\geq\gamma_{2}^{2}(\kappa_{2}-\alpha_{2}\delta)-2\alpha^{2}\delta>\alpha^{2}\delta^{2}$
(15)

Looking at $\rho_{(h)}$ we have

$\displaystyle\ell_{1+\delta}(\rho_{(h)})$
$\displaystyle=\int|\rho_{(h)}(x)-f_{1+\delta}(x)|d\mu_{t+\delta}$
$\displaystyle\leq\int|\rho_{(h)}(x)-f_{1}(x)|d\mu_{t+\delta}+\int|f_{1+\delta}(x)-f_{1}(x)|d\mu_{t+\delta}$
$\displaystyle\leq\int\sup_{x\in\mathscr{X}}|\rho_{(h)}(x)-f_{1}(x)|d\mu_{t+\delta}+\int\sup_{x\in\mathscr{X}}|f_{1+\delta}(x)-f_{1}(x)|d\mu_{t+\delta}$
$\displaystyle=\sup_{x\in\mathscr{X}}|\rho_{(h)}(x)-f_{1}(x)|+\alpha\delta$
(16)

and

$\displaystyle\ell^{2}_{1+\delta}(\rho_{(h)})$
$\displaystyle=\int(\rho_{(h)}(x)-f_{1+\delta}(x))^{2}d\mu_{t+\delta}$
$\displaystyle\leq\int\sup_{x\in\mathscr{X}}|\rho_{(h)}(x)-f_{1+\delta}(x)|^{2}d\mu_{t+\delta}$
$\displaystyle\leq\int\left(\sup_{x\in\mathscr{X}}|\rho_{(h)}(x)-f_{1}(x)|+\sup_{x\in\mathscr{X}}|f_{1+\delta}(x)-f_{1}(x)|\right)^{2}d\mu_{t+\delta}$
$\displaystyle\leq\left(\sup_{x\in\mathscr{X}}|\rho_{(h)}(x)-f_{1}(x)|+\alpha\delta\right)^{2}$
(17)

Given the a choice of $\delta$ satisfying inequality 5, we have

$\epsilon_{3}\stackrel{{\scriptstyle\text{def}}}{{=}}\min\left(\gamma_{1}(\kappa_{1}-\alpha_{2}\delta),\gamma_{2}(\kappa_{2}-\alpha_{2}\delta)\right)-2\alpha\delta>0$
(18)

so we may select $m_{1}$ such that for $m\geq m_{1}$ we have

$P\left(\sup_{x\in\mathscr{X}}|\rho_{(h)}(x)-f_{1}(x)|<\epsilon_{3}\delta\right)>1-\epsilon_{2}$

in which case, given inequalities 12, 14, and 16, expression 6 holds.

Likewise, given a choice of $\delta$ satisfying inequality 7, we have

$\epsilon_{4}\stackrel{{\scriptstyle\text{def}}}{{=}}\sqrt{\min\left(\gamma_{1}^{2}(\kappa_{1}-\alpha_{2}\delta)-2\alpha^{2}\delta,\gamma_{2}^{2}(\kappa_{2}-\alpha_{2}\delta)-2\alpha^{2}\delta\right)}-\alpha\delta>0$
(19)

so we may select $m_{2}$ such that for $m\geq m_{2}$ we have

$P\left(\sup_{x\in\mathscr{X}}|\rho_{(h)}(x)-f_{1}(x)|<\epsilon_{4}\right)>1-\epsilon_{2}$

in which case, given inequalities 13, 15, and 17, expression 8 holds ∎

######  Remark 1 (Robustness of assumptions).

The assumptions of theorem S4 all fundamentally underlie usability of risk
scores in general, with the exception of assumption 1, which is a well-
established effect for many risk scores Tsymbal (2004); Zliobaite (2010).

Assumptions 3,4 state that ‘drift’ must be continuous with time. This is
essentially because continuous drift implies that a risk score that is
consistent at time $t$ will still be useful at a time $t+\epsilon$. Without
such an assumption, risk scores are essentially unusable. We note that the
conditional distribution of intervened outcumes $Y^{g}|X$ need not be
continuous. Obvious extensions to the theorem will give similar results for
almost-continuous drift.

Assumption 2 essentially states that use of the risk score leads to some
change in the risk of an outcome given covariates. Except in special
circumstances, this is typically the very reason for a risk score to be fitted
and used, so this assumption is generally justified.

Assumption 1 essentially states that a non-negligible change occurs over time
in the distribution of the non-intervened outcome conditional on covariates;
that is, should no risk score be in place, the risk to an individual sample
with fixed covariates will vary non-negligibly over time for a positive
proportion of samples. This assumption underlies the need to update risk
scores; if essentially no drift occurs, a risk score need not be updated at
all.

Assumption 2 may be weakened to $|g_{1}(x)-f_{1}(x)|>\gamma_{2}$ if uniform
Lipschitz continuity of $g_{t}$ in $t$ can be assumed.

The processes $\mu_{t}$ and $f_{t}$, which govern drift processes, are not
considered to be random, but the theorem trivially holds if either has random
pertubations for which all elements of their respective domain remain
constrained by the theorem assumptions.

####  S2.1 Simulation

We simulated a population of $2\times 10^{5}$ samples at 50 timepoints, with
ten timepoints per epoch (time between updates). We considered a risk score on
22 ‘visible’ covariates similar to those of the ASPRE score Rolnik et al.
(2017a) with true risk also depending on a ‘latent’ covariate not included in
the risk score. We designated the true risk function $f_{t}$ as a logistic
model with coefficients varying continuously as a Gaussian process of $t$. At
each time point, we computed risk scores using each method, and made
interventions on the 10% of samples with highest predicted risk by reducing
values of visible and latent covariates. We defined total cost as the sum of
post-intervention risk across all samples. Hold-out sets were used in the
final time-point of each epoch.

###  S3 Proof of Theorems 1, 2, and 3

######  Theorem 1.

Suppose assumptions 1-4 in the main manuscript hold. Then there exists a
$N_{*}\in(0,N)$ with $N\in\mathbb{N}$, which we call the optimal holdout set
size, such that:

$\displaystyle\ell(i)\geq\ell(j)$ $\displaystyle\textrm{ for }0<i<j<N_{*}$
$\displaystyle\ell(i)\leq\ell(j)$ $\displaystyle\textrm{ for }N_{*}<i<j<N$

###### Proof.

As discussed in the main manuscript, we may impose that

$\displaystyle\frac{\partial}{\partial n}k_{2}(n)$ $\displaystyle<0$ (20)

Since both $k_{2}(n)$ and $(N-n)$ are positive and monotonically decreasing in
$n$, so is $k_{2}(n)(N-n)$. Now

$\displaystyle\ell^{\prime}(n)$ $\displaystyle=\frac{\partial}{\partial
n}\left(k_{1}n+k_{2}\left(n\right)\left(N-n\right)\right)$ (21)
$\displaystyle=k_{1}+k_{2}^{\prime}\left(n\right)\left(N-n\right)-k_{2}\left(n\right)$
(22)
$\displaystyle=\left(k_{1}-k_{2}\left(n\right)\right)+k_{2}^{\prime}\left(n\right)\left(N-n\right)$
(23)

By assumption 3 in the main manuscript, $k_{1}<k_{2}(0)$, and, from equation
20, $k_{2}^{\prime}(0)<0$, so both terms in equation 23 are negative when
$n=0$ and $\ell^{\prime}(0)<0$. When $n=N$, the second term vanishes while the
first one is positive, as $k_{1}>k_{2}(N)$ by assumption 3 in the main
manuscript. We thus have $\ell^{\prime}(N)>0$. By assumption, $\ell$ is
smooth, so by Bolzano’s Theorem, there must exist at least one point $n_{*}$
for which $\ell^{\prime}\left(n_{*}\right)=0$, which is an extremum of $\ell$.

We now prove that this extremum is unique and a minimum. First, by assumption
4 in the main manuscript, we may impose that

$\frac{\partial^{2}}{\partial n^{2}}k_{2}(n)>0$ (24)

Taking the second derivative of $\ell$

$\displaystyle\frac{\partial^{2}}{\partial n^{2}}\ell(n)$
$\displaystyle=\frac{\partial^{2}}{\partial
n^{2}}\left(k_{1}n+k_{2}(n)(N-n)\right)$ (25)
$\displaystyle=k_{2}^{\prime\prime}(n)(N-n)-2k_{2}^{\prime}(n)$ (26)

and using equations 20 and 24, we see that $\ell^{\prime\prime}(n)$ is
strictly positive, and, as a consequence, $\ell^{\prime}(n)$ is monotonically
increasing. Therefore, the extremum of $\ell(n)$ at $n_{*}$ we found earlier
is unique and, as $\ell^{\prime\prime}(n)>0$, it is a minimum.

If $n_{*}\in 1..(N-1)$, let $N_{*}=n_{*}$. If $n_{*}\notin\mathbb{N}$, let
$N_{*}$ be the closest natural number to either side of $n_{*}$. From
assumption 3 in the main manuscript, $N_{*}$ cannot be $0$ or $N$. In both
scenarios, this completes the proof.

∎

For the proofs of theorems 2 and 3, we will use the following lemma:

######  Lemma S1.

Suppose assumption 1 holds and there exists $0<M<N$ such that
$k_{2}(M)<k_{1}$. Then $\ell(M)<\ell(N)$.

###### Proof.

$\displaystyle k_{2}(M)<k_{1}$ $\displaystyle\implies(N-M)k_{2}(M)<(N-M)k_{1}$
$\displaystyle\implies(N-M)k_{2}(M)+Mk_{1}<Nk_{1}$
$\displaystyle\implies\ell(M)<\ell(N)$

∎

We now have

######  Theorem 2.

Suppose assumptions 1 and 5 hold, and $k_{2}(0)\leq k_{1}$. Then there exists
an $N_{*}\in\\{1,\dots,N-1\\}$ such that: $\ell(i)\geq\ell(N_{*})\textrm{ for
}i\in\\{1,\dots,N-1\\}$ and $\ell(i)>\ell(N_{*})\textrm{ for }i\in\\{0,N\\}$

###### Proof.

All that is needed to show that there exists a holdout set size $M$ where
$\ell(M)<\ell(0)$ and $\ell(M)<\ell(N)$. This will be the $M$ in assumption 5.

It immediately follows from assumption 5 that $k_{2}(M)<k_{1}$, so by Lemma S1
$\ell(M)<\ell(N)$. If $k_{1}=k_{2}(0)$ then we are done as $\ell(N)=\ell(0)$.
If $k_{1}>k_{2}(0)$ then from assumption 5 we have

$\displaystyle k_{1}-k_{2}(0)<\frac{N-M}{N}(k_{1}-k_{2}(M))$
$\displaystyle\implies N(k_{1}-k_{2}(0))<(N-M)(k_{1}-k_{2}(M))$
$\displaystyle\implies Nk_{1}-Nk_{2}(0)<Nk_{1}-(N-M)k_{2}(M)-Mk_{1}$
$\displaystyle\implies\ell(N)-\ell(0)<\ell(N)-\ell(M)$
$\displaystyle\implies\ell(M)<\ell(0)$

as needed ∎

######  Theorem 3.

Suppose assumption 1 holds, $k_{1}<k_{2}(0)$ and there exists $0<M<N$ such
that $k_{2}(M)<k_{1}$. Then there exists an $N_{*}\in\\{1,\dots,N-1\\}$ such
that: $\ell(i)\geq\ell(N_{*})\textrm{ for }i\in\\{1,\dots,N-1\\}$ and
$\ell(i)>\ell(N_{*})\textrm{ for }i\in\\{0,N\\}$.

###### Proof.

Note

$k_{1}<k_{2}(0)\implies Nk_{1}<Nk_{2}(0)\implies\ell(N)<\ell(0)$

so all that is needed is $\ell(M)<\ell(N)$, which is given by Lemma S1. ∎

###  S4 Estimation of $k_{2}(n)$

We argue in this section that $k_{2}(n)$ can generally be a modelled as a
linear function of expected mean squared error of the risk score. Suppose that
we are interested in making predictions at a time $t$. We consider a risk
score $\rho(X)$ which is an inexact approximation of $f_{t}(X)$ (recalling the
definition of $f_{t}$ from section S2 as the probability of $Y=1$ given $X$ at
time $t$ when no risk score is used). We will write
$c_{2}(x,\rho)=\mathbb{E}_{C_{2}}\\{C_{2}(x;d_{n})\\}$, where
$\mathbb{E}_{C_{2}}$ indicates expectation over randomness in actual cost for
a given sample with a given risk score, and $\rho$ a risk score fitted to
samples $d_{n}$.

If it is reasonable to assume that $c_{2}(x,\rho)$ has a straightforward form
in terms of $\rho$, $f_{t}$, then a corresponding form for $k_{2}$ may be
immediate. For instance, if one of

$\displaystyle c_{2}(x,\rho)$ $\displaystyle=c^{0}+c^{1}|\rho-f_{t}(x)|$
$\displaystyle c_{2}(x,\rho)$ $\displaystyle=c^{0}+c^{1}(\rho-f_{t}(x))^{2}$
$\displaystyle c_{2}(x,\rho)$ $\displaystyle=c^{0}+c^{1}f_{t}(x)(\rho-
f_{t}(x))^{2}$

holds, for some constants $c^{0}$, $c^{1}$, then $k_{2}$ will be linear in

$\displaystyle\mathbb{E}_{D_{n}}\left\\{\mathbb{E}_{X}\left[|\rho(X)-f_{t}(X)|\right]\right\\}$
$\displaystyle\mathbb{E}_{D_{n}}\left\\{\mathbb{E}_{X}\left[(\rho(X)-f_{t}(X))^{2}\right]\right\\}$
$\displaystyle\mathbb{E}_{D_{n}}\left\\{\mathbb{E}_{X}\left[f_{t}(X)(\rho(X)-f_{t}(X))^{2}\right]\right\\}$
(27)

respectively. As discussed in the main manuscript, this reduces the estimation
of $k_{2}(n)$ to estimating the ‘learning curve’ of a risk score, and
expectations of risk score accuracy measures such as those above over $X$ and
$D_{n}$ can be readily estimated for small $n$ given training samples $X,Y$.

In more general cases where simple forms of $c_{2}(x,\rho)$ cannot be assumed,
we claim that if $c_{2}(x,\rho)$ is smooth in $\rho$, we should generally
expect $k_{2}(n)$ to be approximately linear in the expected mean-square error
of the risk score.

We work from the following heuristic:

> For any given sample and a range of possible risk scores for that sample,
> the intervention taken will minimise the expected cost for the risk score
> which is unbiased.

This is equivalent to

$\arg\min_{\rho}c_{2}(x,\rho)=f_{t}(x)$ (28)

for all $x$ in the domain of $X$. We suppose firstly that $c_{2}(x,\rho)$ is
smooth in $\rho$, and write

$\displaystyle c_{2}(x,\rho)$
$\displaystyle=c_{2}\left(x,f_{t}(x)\right)+\frac{1}{2}\frac{\partial^{2}c_{2}}{\partial\rho^{2}}\left(x,f_{t}(x)\right)\left(\rho-
f_{t}(x)\right)^{2}+O\left((\rho-f_{t}(x))^{3}\right)$

noting that $\frac{\partial c_{2}}{\partial\rho}\left(x,f_{t}(x)\right)=0$ by
assumption.

The value $\frac{\partial^{2}c_{2}}{\partial\rho^{2}}\left(x,f_{t}(x)\right)$
represents the curvature with respect to $\rho$ of the function
$c_{2}(x,\rho)$ about $\rho=f_{t}(x)$. Practically, this corresponds to the
tolerance or robustness of the intervention: the amount of cost incurred due
to a given deviation of the risk score from $f_{t}(x)$. We claim that this
quantity will thus have relatively low variation across values of $x$, as the
degree of robustness should be roughly constant.

Given this, we have

$\displaystyle\mathbb{E}_{X}\left[c_{2}(X;\rho)\right]$
$\displaystyle=\mathbb{E}_{X}\bigg{[}c_{2}\left(x,f_{t}(x)\right)$
$\displaystyle\phantom{=}+\frac{\partial^{2}c_{2}}{\partial\rho^{2}}\left(x,f_{t}(x)\right)(\rho(X)-f_{t}(X))^{2}$
$\displaystyle\phantom{=}\left.+O\left(\sup_{x}(\rho(x)-f_{t}(x))^{3}\right)\right]$
$\displaystyle\approx\mathbb{E}_{X}\left[c_{2}\left(x,f_{t}(x)\right)\right]$
$\displaystyle\phantom{=}+\mathbb{E}_{X}\left[\frac{\partial^{2}c_{2}}{\partial\rho^{2}}\left(x,f_{t}(x)\right)\right]\mathbb{E}_{X}\left[(\rho(X)-f_{t}(X))^{2}\right]$
$\displaystyle\phantom{=}+O\left(\sup_{x}(\rho(x)-f_{t}(x))^{3}\right)$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}c_{0}+c_{2}\,\textrm{MSE}(\rho)+O\left(\sup_{x}(\rho(x)-f_{t}(x))^{3}\right)$

where $c_{0}$, $c_{2}$ are independent of $\rho$ and $D_{n}$, and
$\textrm{MSE}(\rho)$ is the standard mean-square error of $\rho$. Hence

$\displaystyle k_{2}(n)$
$\displaystyle=\mathbb{E}_{D_{n}}\left\\{\mathbb{E}_{X}\left[c_{2}(X;\rho)\right]\right\\}$
$\displaystyle\approx c_{0}+c_{2}\mathbb{E}_{D_{n}}\left\\{MSE(\rho)\right\\}$

so $k_{2}(n)$ is approximately linear in
$\mathbb{E}_{D_{n}}\left\\{MSE(\rho)\right\\}$.

Suppose that we have a simple setting where we have a single intervention
which we may use, which has a proportional effect on the risk of $Y=1$ (that
is, $g_{t}(x)=(1-\alpha)f_{t}(x)$, with $\alpha<1$). We intervene on a sample
if their risk score exceeds a particular threshold $\rho_{0}$. The cost
function $c_{2}(x,\rho)$ is now discontinuous, but we may simply derive the
form of $k_{2}$ in terms of risk score performance.

We assume that the intervention has a fixed cost $\gamma_{i}$, and that an
event $Y=1$ has a fixed cost $\gamma_{y}$. Then

$\displaystyle c_{2}(x,\rho)$
$\displaystyle=\mathbb{E}_{C_{2}}\left\\{(\textrm{Cost of an
event})+(\textrm{Cost of intervening})\right\\}$
$\displaystyle=\begin{cases}\gamma_{y}f_{t}(x)+0&\textrm{if }\rho<\rho_{0}\\\
\gamma_{y}\alpha f_{t}(x)+\gamma_{i}&\textrm{if }\rho\geq\rho_{0}\end{cases}$
(29)

disregarding potential baseline costs common to all samples. We may now apply
the heuristic above more directly, by presuming that the threshold $\rho_{0}$
is chosen so as to minimise the expectation of $c_{2}(X,\rho)$ over $X$ under
the assumption that $\rho(X)=f_{t}(X)$. In other words, we choose the
threshold that gives us the best outcome assuming the risk score is correct.

This implies that for $f_{t}(x)<\rho_{0}$, we have
$\gamma_{y}f_{t}(x)<\gamma_{y}\alpha f_{t}(x)+\gamma_{i}$ (if true risk is
below the threshold, it is cheaper not to intervene) and for
$f_{t}(x)\geq\rho_{0}$, we have $\gamma_{y}f_{t}(x)\geq\gamma_{y}\alpha
f_{t}(x)+\gamma_{i}$ (if true risk is above the threshold, it is cheaper to
intervene).

We now have:

$\displaystyle c_{2}(x,\rho)-c_{2}(x,f_{t}(x))$
$\displaystyle=\begin{cases}\gamma_{y}f_{t}(x)&\textrm{if
}\rho<\rho_{0},f_{t}(x)<\rho_{0}\\\ \gamma_{y}\alpha
f_{t}(x)+\gamma_{i}-\gamma_{y}f_{t}(x)&\textrm{if
}\rho\geq\rho_{0},f_{t}(x)<\rho_{0}\\\ \gamma_{y}f_{t}(x)-(\gamma_{y}\alpha
f_{t}(x)+\gamma_{i})&\textrm{if }\rho<\rho_{0},f_{t}(x)\geq\rho_{0}\\\
\gamma_{y}\alpha f_{t}(x)+\gamma_{i}&\textrm{if
}\rho\geq\rho_{0},f_{t}(x)\geq\rho_{0}\end{cases}$
$\displaystyle=1_{(\rho<\rho_{0})\,XOR\,(f_{t}(x)<\rho_{0})}|\gamma_{y}(\alpha-1)f_{t}(x)+\gamma_{i}|$

and, denoting $\delta_{\rho}(x)=(\rho(x)<\rho_{0})\,XOR\,(f_{t}(x)<\rho_{0})$,
we have

$\displaystyle k_{2}(n)$
$\displaystyle=k_{2}(n)-\lim_{n\to\infty}k_{2}(n)+\lim_{n\to\infty}k_{2}(n)$
$\displaystyle=\lim_{n\to\infty}k_{2}(n)+\mathbb{E}_{D_{n}}\left\\{\mathbb{E}_{X}\left[c_{2}(X,\rho)-c_{2}(X,f_{t})\right]\right\\}$
$\displaystyle=c^{0}+\mathbb{E}_{D_{n}}\left\\{\mathbb{E}_{X}\left[1_{\delta_{\rho}(X)}|c^{1}f_{t}(x)+c^{2}|\right]\right\\}$

where $c^{0}$, $c^{1}$, $c^{2}$ are constant, and readily estimated from
several observations of $k_{2}$. This form is unsurprising: if a risk score
$\rho(X)$ is such that the sign of $\rho(X)-\rho_{0}$ agrees with the sign of
$f_{t}(X)-\rho_{0}$, it will have identical cost to a risk score $\rho(X)$
which agrees with $f_{t}(X)$ everywhere.

###  S5 Parametric OHS estimation

In this section, we describe estimation of optimal holdout set sizes by
explicit parametrisation of the function $k_{2}(n)$. As in section 5.2 in the
main paper, we take $k_{2}(n)=k_{2}(n;\theta)$. We will take $k_{2}^{\prime}$,
$k_{2}^{\prime\prime}$, $\ell^{\prime}$ to mean partial derivatives with
respect to $n$, and the shorthand $\Theta=(N,k_{1},\theta)$ and
$\Theta_{0}=\mathbb{E}(\Theta)$. We will also write $n_{*}=n_{*}(\Theta)$,
$\ell(n_{*})=\ell\\{n_{*}(\Theta);\Theta\\}$, $n_{0}=n_{*}(\Theta_{0})$ and
$\ell(n_{0})=\ell\\{n_{*}(\Theta_{0}),\Theta_{0}\\}$ for brevity. We presume
that $\Theta$ is an unbiased estimate of $\Theta_{0}$, so $\Theta_{0}$
corresponds to ‘true’ parameter values.

We firstly develop asymptotic confidence intervals for parametric OHS
estimates to link error in parameter estimates to error in optimal size. The
sample-size $m$ used in the following denotes a proxy for effort expended in
estimating $\Theta_{0}$.

######  Theorem S5.

Assume that $k_{2}^{\prime\prime}(n;\theta)$, $k_{2}^{\prime}(n;\theta)$ and
$\nabla_{\theta}k_{2}(n;\theta)$ are continuous in $n$ and $\theta$ in some
neighbourhood of $(n_{0},\Theta_{0})$, and that $\Theta_{0}$ parametrizes a
setting satisfying assumptions 1-4. Suppose that $\Theta$ behaves as a mean of
of $m$ appropriately-distributed samples in satisfying
$\sqrt{m}(\Theta-\Theta_{0})\to N\left(0,\Sigma\right)$ in distribution where
$\Theta_{0}$ does not depend on $m$, that an estimate $\widehat{\Sigma}$ of
$\Sigma$ is available which is independent of $\Theta$ and satisfies
$||\widehat{\Sigma}-\Sigma||_{2}\to 0$ in distribution, and that $n_{0}$ is
finite and unique as above. Then denoting

$\beta_{\Theta}=\frac{\partial^{2}\ell}{\partial
n\partial\Theta_{i}}\bigg{/}\frac{\partial^{2}\ell}{\partial n^{2}},\hskip
30.0pt\gamma_{\Theta}=\frac{\partial\ell}{\partial\Theta_{i}}$

we may uniquely define
$n_{0}=\left\\{n:\ell^{\prime}(n;\Theta_{0})=0\right\\}$ and we have

$\sqrt{m}(n_{*}-n_{0})\to
N\left(0,\beta_{\Theta_{0}}^{t}\Sigma\beta_{\Theta_{0}}\right),\hskip
30.0pt\sqrt{m}\left\\{\ell(n_{*})-\ell(n_{0})\right\\}\to
N\left(0,\gamma_{\Theta_{0}}^{t}\Sigma\gamma_{\Theta_{0}}\right)$

in distribution, and denoting $z_{\alpha}=\Phi^{-1}\left(1-\alpha/2\right)$,
the confidence intervals

$I_{\alpha}(\Theta,\hat{\Sigma})=\left[n_{*}(\Theta)\pm
z_{\alpha}\sqrt{\frac{\beta_{\Theta}^{t}\widehat{\Sigma}\beta_{\Theta}}{m}}\hskip
5.0pt\right],\hskip 20.0ptJ_{\alpha}(\Theta,\hat{\Sigma})=\left[\ell(n_{*})\pm
z_{\alpha}\sqrt{\frac{\gamma_{\Theta}^{t}\widehat{\Sigma}\gamma_{\Theta}}{m}}\hskip
5.0pt\right]$

satisfy $P\left\\{n_{0}\in I_{\alpha}(\Theta,\hat{\Sigma})\right\\}\to
1-\alpha$ and $P\left\\{\ell(n_{0})\in
J_{\alpha}(\Theta,\hat{\Sigma})\right\\}\to 1-\alpha$ as $m\to\infty$.

The proof is given in Supplement S5.2 below. A consequence is that for
sufficiently accurately estimated costs, the OHS will be a non-trivial size:

######  Corollary S2.

Under assumptions of Theorems 1, S5, $P\\{1<n_{*}(\Theta)<N\\}\to 1$ as
$m\to\infty$.

In light of the proportionality assumption in Section 5.1, and the tendency of
the accuracy of a risk score with number of training samples (‘learning
curve’) to follow a power-law form (Viering and Loog, 2021), we recommend
considering such a parametric form for $k_{2}$ (i.e.
$k_{2}(n;\theta)=an^{-b}+c$ with $\theta=(a,b,c)$), and provide explicit
asymptotic confidence intervals for this setting in Supplement S5.1. Examples
of variation in $n_{*}$ and $\ell(n_{*})$ with a power-law form for $k_{2}$,
are shown in Supplementary Figures 8, 9.

Note that confidence intervals must be interpreted with with care: if the
sampling distributions for $k_{1}$ and $\theta$ admit the possibility that
assumptions of Theorem 1 are violated such that
$P\left[k_{1}<\lim\inf_{n\to\infty}\\{k_{2}(n,\theta)\\}\right]>0$ then the
standard error of $n_{*}$ does not exist, as $n_{*}$ can be undefined. Finite-
sample confidence intervals may be constructed by bootstrapping (see function
ci_ohs() in our R package OptHoldoutSize).

Our parametric algorithm assumes $\Theta$ is estimated from a multiset
$\mathbf{n}$ of values in $\\{1,\dots,N\\}$ and estimates $\mathbf{d}$ of
$k_{2}(n)$ for each $n\in\mathbf{n}$ with known finite sampling variances
$\boldsymbol{\sigma}^{2}$. For certain multisets $\mathbf{n}$, estimates of
$\Theta$ will not converge; for instance, if $\mathbf{n}$ contains only a
single value repeated. The value $m$ in Theorem S5 should be interpreted as an
‘effective’ population size, such that
$\sqrt{m}\left\\{\Theta(\mathbf{n})-\Theta_{0}\right\\}\to
N\left(0,\Sigma\right)$ in distribution.

Given that our eventual aim to estimate the OHS with minimal error, we suggest
the following way to iteratively select a new value $\tilde{n}$ at which an
estimate $\hat{k_{2}}(\tilde{n})$ of $k_{2}(\tilde{n})$ should be made, given
a set $\mathbf{n}$ of points at which estimates $\mathbf{k_{2}}$ of
$k_{2}(\mathbf{n})$ have been made already. We denote by
$\Theta(\mathbf{n},\mathbf{k_{2}},\boldsymbol{\sigma})$,
$\hat{\Sigma}(\mathbf{n},\mathbf{k_{2}},\boldsymbol{\sigma})$ and
$I_{\alpha}(\mathbf{n},\mathbf{k_{2}},\boldsymbol{\sigma})$ respectively the
estimates of $\Theta_{0}$,
$\lim_{m\to\infty}\textrm{var}\left[\sqrt{m}\left\\{\Theta(\mathbf{n},\mathbf{k_{2}},\boldsymbol{\sigma})-\Theta_{0}\right\\}\right]$
and the width of the confidence interval
$I_{\alpha}\left\\{\Theta(\mathbf{n},\mathbf{k_{2}},\boldsymbol{\sigma}),\hat{\Sigma}(\mathbf{n},\mathbf{k_{2}},\boldsymbol{\sigma})\right\\}$.
Suppose we have the option of estimating $d(n)$ for one value of
$n\in\\{1,\dots,N\\}$ with known variance
$\textrm{var}\left\\{d(n)\right\\}=\sigma^{2}$. We select $\tilde{n}$ as:

$\tilde{n}=\arg\min_{n}\mathbb{E}_{d(n)\sim
N\left[k_{2}\left\\{n,\Theta(\mathbf{n},\mathbf{k_{2}},\boldsymbol{\sigma})\right\\},\sigma^{2}\right]}\left[I_{\alpha}\left\\{\mathbf{n}\cup
n,\mathbf{k_{2}}\cup\hat{k_{2}}(n),\boldsymbol{\sigma}\cup\sigma\right\\}\right]$
(31)

that is, ‘select the $\tilde{n}$ which will minimize the expected OHS
confidence interval width if added to our set $\mathbf{n}$, with expectation
computed with respect to our current parameter estimates’. If no minimum
exists, $\tilde{n}$ is selected uniformly from ${1,\dots,N}$. Algorithm 3, an
expanded version of algorithm 1 in the main manuscript, shows our full
parametric estimation procedure.

1 $\mathbf{n},\mathbf{k_{2}},\boldsymbol{\sigma}^{2}\leftarrow$ some initial
values $\mathbf{n}$ with $(\boldsymbol{k_{2}})_{i}=\hat{k_{2}}(n_{i})\approx
k_{2}(n_{i})$,
$(\boldsymbol{\sigma}^{2})_{i}=\textrm{var}(\hat{k_{2}}(n_{i}))$;

2 while  _$|\mathbf{n}| <n_{add}$_ do

3 Find best new value $\tilde{n}$ to add to $\mathbf{n}$ as per formula 31 ;

4 Estimate $\hat{k_{2}}(\tilde{n})\approx k_{2}(\tilde{n})$ ;

5 $\mathbf{n}\leftarrow(\mathbf{n}\cup\tilde{n})$,
$\mathbf{k_{2}}\leftarrow\left\\{\mathbf{k_{2}}\cup\hat{k_{2}}(\tilde{n})\right\\}$,
$\boldsymbol{\sigma}^{2}\leftarrow\boldsymbol{\sigma}^{2}\cup\textrm{var}\left\\{\hat{k_{2}}(\tilde{n})\right\\}$
;

6

7 end while

8Re-estimate OHS
$n_{*}^{final}=n_{*}\left\\{\Theta\left(\mathbf{n},\mathbf{k_{2}},\boldsymbol{\sigma}\right)\right\\}$
;

return $n_{*}^{final}$

Algorithm 3 Parametric OHS estimation; with $n_{add}$ estimates of
$k_{2}(\cdot)$

####  S5.1 Explicit partial derivatives for $n^{*}$, $\ell$ with power-law
parametrisation

If we assume a power-law form of $k_{2}$, parametrised by
$\theta=(a,b,c,k_{1},N)$;

$k_{2}(n;\theta)=an^{-b}+c$ (32)

then we have

$\displaystyle\frac{\partial n_{*}}{\partial a}$
$\displaystyle=\frac{1}{a}\left(\frac{bNn_{*}-(b-1)n_{*}^{2}}{b(b+1)N-b(b-1)n_{*}}\right)$
$\displaystyle\frac{\partial n_{*}}{\partial b}$
$\displaystyle=\frac{Nn_{*}(b\log(n_{*})-1)-n_{*}^{2}\left((b-1)\log(n_{*})-1\right)}{b(b+1)N-b(b-1)n_{*}}$
$\displaystyle\frac{\partial n_{*}}{\partial c}$
$\displaystyle=\frac{1}{a}\left(\frac{n_{*}^{b+2}}{b(b+1)N-b(b-1)n_{*}}\right)$
$\displaystyle\frac{\partial n_{*}}{\partial k_{1}}$
$\displaystyle=\frac{1}{a}\left(\frac{-n_{*}^{b+2}}{b(b+1)N-b(b-1)n_{*}}\right)$
$\displaystyle\frac{\partial n_{*}}{\partial N}$
$\displaystyle=\frac{bn_{*}}{b(b+1)N-b(b-1)n_{*}}$

and, more simply

$\displaystyle\frac{\partial}{\partial a}\ell(n_{*};\theta)$
$\displaystyle=(N-n_{*})n_{*}^{-b}$ $\displaystyle\frac{\partial}{\partial
b}\ell(n_{*};\theta)$ $\displaystyle=-\log(n_{*})(N-n_{*})an_{*}^{-b}$
$\displaystyle\frac{\partial}{\partial c}\ell(n_{*};\theta)$
$\displaystyle=N-n_{*}$ $\displaystyle\frac{\partial}{\partial
k_{1}}\ell(n_{*};\theta)$ $\displaystyle=n_{*}$
$\displaystyle\frac{\partial}{\partial N}\ell(n_{*};\theta)$
$\displaystyle=an_{*}^{-b}+c$

####  S5.2 Proof of Theorem S5

######  Theorem S5.

Assume that $k_{2}^{\prime\prime}(n;\theta)$, $k_{2}^{\prime}(n;\theta)$ and
$\nabla_{\theta}k_{2}(n;\theta)$ are continuous in $n$ and $\theta$ in some
neighbourhood of $(n_{0},\Theta_{0})$, and that $\Theta_{0}$ parametrizes a
setting satisfying assumptions 1-4. Suppose that $\Theta$ behaves as a mean of
of $m$ appropriately-distributed samples in satisfying
$\sqrt{m}(\Theta-\Theta_{0})\to N\left(0,\Sigma\right)$ in distribution where
$\Theta_{0}$ does not depend on $m$, that an estimate $\widehat{\Sigma}$ of
$\Sigma$ is available which is independent of $\Theta$ and satisfies
$||\widehat{\Sigma}-\Sigma||_{2}\to 0$ in distribution, and that $n_{0}$ is
finite and unique as above. Then denoting

$\beta_{\Theta}=\frac{\partial^{2}\ell}{\partial
n\partial\Theta_{i}}\bigg{/}\frac{\partial^{2}\ell}{\partial n^{2}},\hskip
30.0pt\gamma_{\Theta}=\frac{\partial\ell}{\partial\Theta_{i}}$

we may uniquely define
$n_{0}=\left\\{n:\ell^{\prime}(n;\Theta_{0})=0\right\\}$ and we have

$\sqrt{m}(n_{*}-n_{0})\to
N\left(0,\beta_{\Theta_{0}}^{t}\Sigma\beta_{\Theta_{0}}\right),\hskip
30.0pt\sqrt{m}\left\\{\ell(n_{*})-\ell(n_{0})\right\\}\to
N\left(0,\gamma_{\Theta_{0}}^{t}\Sigma\gamma_{\Theta_{0}}\right)$ (33)

in distribution, and denoting $z_{\alpha}=\Phi^{-1}\left(1-\alpha/2\right)$,
the confidence intervals

$I_{\alpha}(\Theta,\hat{\Sigma})=\left[n_{*}(\Theta)\pm
z_{\alpha}\sqrt{\frac{\beta_{\Theta}^{t}\widehat{\Sigma}\beta_{\Theta}}{m}}\hskip
5.0pt\right],\hskip 20.0ptJ_{\alpha}(\Theta,\hat{\Sigma})=\left[\ell(n_{*})\pm
z_{\alpha}\sqrt{\frac{\gamma_{\Theta}^{t}\widehat{\Sigma}\gamma_{\Theta}}{m}}\hskip
5.0pt\right]$

satisfy $P\left\\{n_{0}\in I_{\alpha}(\Theta,\hat{\Sigma})\right\\}\to
1-\alpha$ and $P\left\\{\ell(n_{0})\in
J_{\alpha}(\Theta,\hat{\Sigma})\right\\}\to 1-\alpha$ as $m\to\infty$.

As above, we consider $n_{*}$ as a function of parameters
$\Theta=(N,k_{1},\theta)$ (where $k_{2}(\cdot)=k_{2}(\cdot;\theta)$), write
$n_{*}=n_{*}(\Theta)$, set $\Theta_{0}=E(\Theta)$, $n_{0}=n_{*}(\Theta_{0})$
and $\ell(n_{0})=\ell(n_{0};\Theta_{0})$ and
$\ell(n_{*})=\ell(n_{*}(\Theta),\Theta)$. As discussed above, $n_{*}$ and
$\ell(n_{*})$ do not generally have means or standard errors.

###### Proof.

From $\ell(n)=k_{1}n+k_{2}(n;\theta)(N-n)$ and
$n_{*}=\\{n:\ell^{\prime}(n;\Theta)=0\\}$, where such $n_{*}$ is unique, we
have (as per section 5.2)

$\displaystyle(\nabla n_{*})_{i}=\frac{\partial
n_{*}}{\partial\Theta_{i}}=\frac{\frac{\partial^{2}l}{\partial
n\partial\Theta_{i}}}{\frac{\partial^{2}\ell}{\partial
n^{2}}}=(\beta_{\Theta_{0}})_{i}$
$\displaystyle\left(\nabla\ell(n_{*})\right)_{i}=\frac{\partial\ell}{\partial\Theta_{i}}=(\gamma_{\Theta_{0}})_{i}$

for all components $\Theta_{i}$ of $\Theta$. Thus partial derivatives of
$n_{*}$ exist as long as

$\frac{\partial^{2}\ell}{\partial n^{2}}>0$ (34)

By assumption, $\ell(\cdot;\Theta_{0})$ has a minimum at $n_{0}$. Since

$\frac{\partial^{2}\ell}{\partial n^{2}}=\frac{\partial^{2}}{\partial
n^{2}}k_{2}(n;\theta)-2\frac{\partial}{\partial n}k_{2}(n;\theta)$ (35)

where both terms are continuous in a neighbourhood of $n_{0}$, $\Theta_{0}$ by
assumption, the value of $\frac{\partial^{2}\ell}{\partial n^{2}}$ must be
positive in some (possibly smaller) neighbourhood $R_{\delta}$ of
$(n_{0},\Theta_{0})$ of width $2\delta$, and hence all partial derivatives of
$n_{*}$ and $\ell(n_{*})$ are defined (and indeed continuous) in $R_{\delta}$.
Within $R_{\delta}$ we have

$\displaystyle n_{*}(\Theta)$ $\displaystyle=n_{*}(\Theta_{0})+(\nabla
n_{*}|_{\Theta=\Theta_{0}})\cdot\left(\Theta-\Theta_{0}\right)+O\left(||\Theta-\Theta_{0}||_{2}\right)$
$\displaystyle=n_{0}+\beta_{\Theta_{0}}^{t}\cdot\left(\Theta-\Theta_{0}\right)+O\left(||\Theta-\Theta_{0}||_{2}\right)$
(36) $\displaystyle\ell(n_{*})$
$\displaystyle=\ell(n_{*}(\Theta_{0});\Theta_{0})+(\nabla\ell(n_{*})|_{\Theta=\Theta_{0}})\cdot\left(\Theta-\Theta_{0}\right)+O\left(||\Theta-\Theta_{0}||_{2}\right)$
$\displaystyle=\ell(n_{0})+\gamma_{\Theta_{0}}^{t}\cdot\left(\Theta-\Theta_{0}\right)+O\left(||\Theta-\Theta_{0}||_{2}\right)$
(37)

from which, given the assumption of asymptotic normality of $\Theta$,
assertions 33 follow. We note that despite this convergence in distribution,
$n_{*}$ and $\ell(n_{*})$ do not generally have first or second moments for
finite $m$.

We now have

$\displaystyle P\left(n_{0}\geq
n_{*}(\Theta)+z_{\alpha}\sqrt{\frac{\beta_{\Theta}\widehat{\Sigma}\beta_{\Theta}^{t}}{m}}\right)$
$\displaystyle=P\left(\frac{\sqrt{m}}{z_{\alpha}}(n_{0}-n_{*}(\Theta))\geq\sqrt{\beta_{\Theta}^{t}\widehat{\Sigma}\beta_{\Theta}}\right)$
$\displaystyle=P\left(\frac{\sqrt{m}}{z_{\alpha}}(n_{0}-n_{*}(\Theta))\geq\left(\beta_{\Theta_{0}}^{t}\Sigma\beta_{\Theta_{0}}+\right.\right.$
$\displaystyle\phantom{=}\phantom{P((n_{0}-n_{*}(\Theta))\geq}\beta_{\Theta}^{t}\left(\widehat{\Sigma}-\Sigma\right)\beta_{\Theta}+$
$\displaystyle\phantom{=}\left.\left.\phantom{P((n_{0}-n_{*}(\Theta))\geq}(\beta_{\Theta}-\beta_{\Theta_{0}})^{t}\Sigma\left(\beta_{\Theta}+\beta_{\Theta_{0}}\right)\right)^{\frac{1}{2}}\right)$
$\displaystyle\to
P\left(\frac{\sqrt{m}}{z_{\alpha}}(n_{0}-n_{*}(\Theta))\geq\sqrt{\beta_{\Theta_{0}}^{t}\Sigma\beta_{\Theta_{0}}}\right)$
$\displaystyle=\frac{\alpha}{2}$ (38)

since, by the assumption of convergence of $\widehat{\Sigma}$

$\displaystyle\left|\beta_{\Theta}^{t}\left(\Sigma-\widehat{\Sigma}\right)\beta_{\Theta}\right|$
$\displaystyle\leq||\beta_{\Theta}||_{2}||\Sigma-\widehat{\Sigma}||_{2}$
$\displaystyle\to_{p}0$ (39)

and, since $P(\Theta\in R_{\delta})\to 1$ by the asymptotic normality of
$\Theta$, we have from 36

$\displaystyle\left|(\beta_{\Theta}-\beta_{\Theta_{0}})^{t}\Sigma\left(\beta_{\Theta}+\beta_{\Theta_{0}}\right)\right|$
$\displaystyle=O\left(||\beta_{\Theta}-\beta_{\Theta_{0}}||_{2}\right)$
$\displaystyle\to_{p}0$ (40)

Thus, combining with the corresponding limit for the lower end of
$I_{\alpha}(\Theta,\hat{\Sigma})$:

$P(n_{0}\in I_{\alpha}(\Theta,\hat{\Sigma}))\to 1-\alpha$ (41)

as required. An identical argument holds for
$J_{\alpha}(\Theta,\hat{\Sigma})$.

∎

###  S6 Estimation of OHS by Bayesian Emulation

Our second algorithm for estimation of optimal holdout sizes uses Bayesian
emulation (Brochu et al., 2010). In many cases, it may be difficult or
unrealistic to provide a precise parametric form for the function $k_{2}(n)$.
The function depends both on the ‘learning curve’ of the risk score, which may
be complex (Viering and Loog, 2021), and the relationship of the risk score
accuracy to the accrued cost, which may be nonlinear. Here, we propose a
second algorithm which is less reliant on assuming a particular parametric
form for $k_{2}(n)$.

As in the main manuscript, we approximate the cost function $\ell$ as an
‘emulator’ modelled as a Gaussian process, and take the minimum of its
posterior mean over $n$ as our OHS estimate. It is worth noting that whilst
gaining an accurate approximation of the cost function is important, the main
goal is to ascertain the minimum of this function, not provide a universally
effective approximation at all points. Therefore, we aim to choose the
location of design points $\mathbf{n}$ in order to efficiently obtain the
minimum of the cost function, and hence the OHS.

First we must construct an emulator which approximates the cost function. We
begin with an initial set of design points $\mathbf{n}$ and their
corresponding observed cost estimates $\mathbf{d}$ with sampling variances
$\boldsymbol{\sigma}^{2}$, noting that $\boldsymbol{\sigma}$ has a slightly
different meaning to that in section 5.2. The prior for our emulator is,
following Vernon et al. (2018), $\ell(n)=m(n,\Theta)+u(n)$ with mean function
$m(n,\Theta)=k_{1}n+k_{2}(n;\theta)(N-n)$, given some initial estimate of
$\Theta=(N,k_{1},\theta)$, and $u(n)$ a zero-mean Gaussian process

$u(n)\sim\mathcal{GP}\left\\{0,k(n,n^{\prime})\right\\}\hskip
20.0ptk(n,n^{\prime})=\sigma_{u}^{2}\exp{\left\\{-\left(\frac{n-n^{\prime}}{\zeta}\right)^{2}\right\\}}$

where $k$ is chosen to enforce smoothness in $\ell(n)$, though other
covariance functions having varying degrees of smoothness could be used. The
hyperparameters $\theta$, $\sigma_{u}$ and $\zeta$ are problem-specific and
must be specified; however, we will show that for sufficiently large
$|\mathbf{n}|$ mis-specification of $\theta$, $\sigma_{u}$ and $\theta$ is
overcome.

We denote $d_{i}\sim N\left\\{l(n_{i}),\sigma_{i}^{2}\right\\}$ where
$d_{i}=(\mathbf{d})_{i}$ is the $i$th element of $\mathbf{d}$, etc. Since
$\mathbf{n}$ may be a multiset, we take $\mathbf{n}^{1}$ as the set of unique
values in $\mathbf{n}$, with $\mathbf{d}^{1}$, $\boldsymbol{\sigma}^{1}$
defined correspondingly with $d^{1}_{i}$ as an inverse-variance weighted mean
of $\\{d_{j}:n_{j}=n^{1}_{i}\\}$ with sample variance $(\sigma^{1}_{i})^{2}$:

$d^{1}_{i}=\frac{\sum_{j:n_{j}=n^{1}_{i}}d_{j}\sigma_{j}^{-2}}{\sum_{j:n_{j}=n^{1}_{i}}\sigma^{-2}_{j}}\hskip
20.0pt(\sigma^{1}_{i})^{2}=\sum_{j:n_{j}=n^{1}_{i}}\sigma^{-2}_{j}$ (42)

noting that $\sigma^{1}_{i}$ may change with $i$. Alternatively, we may
account for the variation in $\mathbf{d}$ through ‘inactive’ variables and opt
to use a ’nugget’ term; this approach is described in detail in Supplement
S6.1.

Now with input $n$, with an unevaluated loss value, our emulator specifies
that the joint distribution of $\ell(n)$ and our observed output values
$\mathbf{d}^{1}$ is:

$\begin{bmatrix}\ell(n)\\\
\mathbf{d^{1}}\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}m(n,\Theta)\\\
m(\mathbf{n^{1}},\Theta)\end{bmatrix},\begin{bmatrix}k(n,n)&k(n,\mathbf{n^{1}})\\\
k(\mathbf{n^{1}},n)&k(\mathbf{n^{1}},\mathbf{n^{1}})+\textrm{diag}\\{(\boldsymbol{\sigma^{1}})^{2}\\}\end{bmatrix}\right)$

where $m(\mathbf{n^{1}},\Theta)^{T}_{i}=m(n^{1}_{i},\Theta)$,
$k(n,\mathbf{n^{1}})_{i}=k(\mathbf{n^{1}},n)^{T}_{i}=k(n,n^{1}_{i})$,
$k(\mathbf{n},\mathbf{n})_{ij}=k(n^{1}_{i},n^{1}_{j})$,
$\textrm{diag}\\{(\boldsymbol{\sigma^{1}})^{2}\\}_{ij}=(\sigma^{1}_{i})^{2}1_{i=j}$.
By obtaining the conditional posterior distribution
$\pi_{\mathbf{n}}=\pi\\{\ell(n)\mid
n,\mathbf{n^{1}},\mathbf{d^{1}},\boldsymbol{\sigma^{1}}\\}$ and taking the
expectation and variance we gain the Bayes linear update equations (Vernon et
al., 2018):

$\displaystyle\mu(n)$
$\displaystyle=\mathbb{E}_{\pi_{\mathbf{n}}}\\{\ell(n)\\}=m(n,\Theta)+k(n,\mathbf{n^{1}})\left[k(\mathbf{n^{1}},\mathbf{n^{1}})+\textrm{diag}\\{(\boldsymbol{\sigma^{1}})^{2}\\}\right]^{-1}\\{\mathbf{d^{1}}-m(\mathbf{n^{1}},\Theta)\\}$
$\displaystyle\Psi(n)$
$\displaystyle=\textrm{var}_{\pi_{\mathbf{n}}}\\{\ell(n)\\}=k(n,n)-k(n,\mathbf{n^{1}})[k(\mathbf{n^{1}},\mathbf{n^{1}})+\textrm{diag}\\{(\boldsymbol{\sigma^{1}})^{2}\\}]^{-1}k(\mathbf{n^{1}},n)$

In algorithm 3, selection of new design points should generally favour well-
spaced points across $\\{1,\dots,N\\}$ for both exploration and exploitation.
Here, since we wish both to estimate the OHS accurately but also locally
approximate $\ell$ well, we choose the next $n$ in a way which predominantly
but not completely favours exploitation. We use the ‘expected improvement’,
which measures discrepancy between the emulator at a certain design point and
the known minimum $EI(\cdot)$ (Brochu et al., 2010):

$EI(n)=\\{d^{-}-\mu(n)\\}\Phi\left(\frac{d^{-}-\mu(n)}{\sqrt{\Psi(n)}}\right)+\sqrt{\Psi(n)}\phi\left(\frac{d^{-}-\mu(n)}{\sqrt{\Psi(n)}}\right)$

where $d^{-}=\min_{i}\\{\mathbf{d^{1}}_{i}\\}$, and

$\tilde{n}=\arg\max_{n\in\\{1,\dots,N\\}}EI(n)=\arg\max\left(\mathbb{E}_{\pi_{\mathbf{n}}}\left[\max\\{0,d^{-}-\ell(n)\\}\right]\right)$

We see that by formulating the problem in terms of $EI(\cdot)$, there is a
natural stopping criterion on the size of $\mathbf{n}$: setting a threshold
$EI(\tilde{n})>\tau$ allows us to specify that for each iteration that we
expect total cost to improve by at least $\tau$ over our current known minimum
$d^{-}$. Examples of $\mu(\cdot),\Psi(\cdot),EI(\cdot)$ are shown in
Supplementary Figure 10. This leads to algorithm 4 for OHS estimation by
Bayesian Emulation (a more precise version of algorithm 2 in the main
manuscript).

1 $\mathbf{n},\mathbf{d},\boldsymbol{\sigma}^{2}\leftarrow$ some initial
values $\mathbf{n}$ with $(\boldsymbol{d})_{i}=d(n_{i})\approx\ell(n)$,
$(\boldsymbol{\sigma}^{2})_{i}=\textrm{var}\\{d(n_{i})\\}$ ;

2 Coalesce $\mathbf{n},\mathbf{d},\boldsymbol{\sigma}$ into
$\mathbf{n^{1}},\mathbf{d^{1}},\boldsymbol{\sigma^{1}}$ as above ;

3 Estimate functions $\mu(n)$, $\Psi(n)$, $EI(n)$, with
$\Theta=\Theta(\mathbf{n^{1}},\mathbf{d^{1}},\boldsymbol{\sigma^{1}})$ ;

4 while  _$\max_{n\in\\{1,\dots,N\\}}\\{EI(n)\\} >\tau$_ do

5 $\tilde{n}\leftarrow\arg\max_{n\in\\{1,\dots,N\\}}EI(n)$ ;

6 Estimate $d(\tilde{n})\approx k_{2}(\tilde{n})$ ;

7 $\mathbf{n}\leftarrow(\mathbf{n}\cup\tilde{n})$;
$\mathbf{d}\leftarrow(\mathbf{d}\cup d\\{\tilde{n})\\}$;
$\boldsymbol{\sigma}^{2}\leftarrow\left[\boldsymbol{\sigma}^{2}\cup\textrm{var}\\{d(\tilde{n})\\}\right]$
;

8 Coalesce $\mathbf{n},\mathbf{d},\boldsymbol{\sigma}$ into
$\mathbf{n^{1}},\mathbf{d^{1}},\boldsymbol{\sigma^{1}}$ ;

9 Re-estimate functions $\mu(n)$, $\Psi(n)$, $EI(n)$, with
$\Theta=\Theta(\mathbf{n^{1}},\mathbf{d^{1}},\boldsymbol{\sigma^{1}})$ ;

10

11 end while

return $n_{*}^{final}=\arg\min_{n\in\\{1,\dots,N\\}}\left\\{\mu(n)\right\\}$

Algorithm 4 Emulation OHS estimation; minimum cost improvement $\tau$

Various results on the consistency of the expected improvement algorithm have
been proved, albeit in differing settings; either with noiseless observations
$\mathbf{d}$ (Locatelli, 1997; Vazquez and Bect, 2010; Bull, 2011) or with
noisy observations with known variance (Ryzhov, 2016). We prove the following
consistency results specifically for the setting of this work in Supplement
S6.2.

######  Theorem S6.

If $\ell(n)$, $\boldsymbol{\sigma}$, and $m(n,\Theta)$ are almost surely
bounded and $d_{i}\sim N\left\\{l(n_{i}),\sigma_{i}^{2}\right\\}$ then for
every $n\in\\{1,\ldots,N\\}$, as the multiplicity of $n$ in $\mathbf{n}$ tends
to $\infty$ we have $\mu(n)\longrightarrow\ell(n)$ and $\Psi(n)\to 0$ almost
surely with respect to variation in $\mathbf{d}$.

also noting the following simple result, proved in Supplement S6.3:

######  Theorem S7.

Given the conditions of Theorem S6, for every $n\in\\{1,\ldots,N\\}$, as the
multiplicity of $n$ in $\mathbf{n}$ tends to $\infty$,

$EI(n)\longrightarrow 0$

almost surely with respect to randomness in $\mathbf{d}$

These results assert that $\mu(n)$ can eventually approximate any loss
function sufficiently well given enough estimates of $\ell$ at all values of
$n$. It is not obvious that this is guaranteed by algorithm 4, although we
show that this generally does occur in the following, the proof of which is
given in Supplement S6.4:

######  Theorem S8.

If $\ell(n)$, $\boldsymbol{\sigma}$, and $m(n,\Theta)$ are almost surely
bounded and $d_{i}\sim N\left\\{l(n_{i}),\sigma_{i}^{2}\right\\}$ then under
algorithm 4 with $\tau=0$, the value $\mu(\tilde{n})$ converges almost surely
to $\ell(\tilde{n})$ for every $\tilde{n}\in\\{1,\ldots,N\\}$.

We characterize the error in $n_{*}$ using ‘the number of values of $n$ for
which the probability of the true cost at holdout set size $n$ is less than
the estimated minimum cost exceeds $1-\alpha$’, or formally:
$\left\\{n:pr_{\textbf{n}_{\pi}}\left\\{\ell(n)<\mu(n_{*})\right\\}\geq
1-\alpha\right\\}$, although this should not be interpreted as a credible set
for $n_{*}$. This is implemented in our R package OptHoldoutSize, available on
CRAN.

####  S6.1 Emulation of cost function with nugget term

Rather than explaining the variation of values in $\mathbf{d}$ corresponding
to a design point in $\mathbf{n^{1}}$ as approximation error of a
deterministic loss function, we can explain this variation as the result of
not including active variables, being the data $(X,Y)$. Note that as a
consequence we are now not emulating a deterministic function $\ell(n)$ as we
are not generalising the loss through expectations, we are generalising the
loss through omission of the data which generated $\mathbf{d}$. To clarify
this distinction we replace the loss function $\ell(n)$ with the stochastic
function $\mathcal{E}(n)$.

Now we may specify variation in $\mathbf{d}$ using a ‘nugget’ term $w(n)$,
following Bower et al. (2010):

$\mathcal{E}(n)=m(n)+u(n)+w(n)$ (43)

where $m(n)$ and $u(n)$ are as before but now $w(n)$ represents our nugget
term, which we again specify as a Gaussian process:

$w(n)\sim\mathcal{GP}(0,\kappa(n,n^{\prime}))$ (44)

with

$\kappa(n,n^{\prime})=\begin{cases}\kappa(n)&\mbox{if }n=n^{\prime}\\\
0&\mbox{otherwise}\end{cases}$ (45)

Since there is less variance in risk scores fitted to larger datasets, we
expect less variance in $\mathcal{E}(n)$ for larger $n$, so we specify
$\kappa(n)$ as a monotonically decreasing function in $n$.

The joint distribution between $\mathcal{E}(n)$ and $\mathbf{d^{1}}$ is now:

$\begin{bmatrix}\mathcal{E}(n)\\\
\mathbf{d^{1}}\end{bmatrix}\sim\mathcal{N}\left(\begin{bmatrix}m(n)\\\
m(\mathbf{n^{1}})\end{bmatrix},\begin{bmatrix}k(n,n)+\kappa(n)&k(n,\mathbf{n^{1}})\\\
k(\mathbf{n^{1}},n)&k(\mathbf{n^{1}},\mathbf{n^{1}})+\textrm{diag}(\kappa(\mathbf{n^{1}}))\end{bmatrix}\right)$
(46)

This then gives our Bayes linear update equations in terms of
$\pi_{\mathbf{n}}=\pi(\mathcal{E}(n)|n,\mathbf{n^{1}},\mathbf{d^{1}})$ as

$\displaystyle\mu(n)$
$\displaystyle=\mathbb{E}_{\pi_{\mathbf{n^{1}}}}(\mathcal{E}(n))$
$\displaystyle=m(n)+k(n,\mathbf{n^{1}})[k(\mathbf{n^{1}},\mathbf{n^{1}})+\textrm{diag}(\kappa(\mathbf{n^{1}}))]^{-1}(\mathbf{d^{1}}-m(\mathbf{n^{1}}))$
(47) $\displaystyle\Psi(n)$
$\displaystyle=\textrm{var}_{\pi_{\mathbf{n^{1}}}}(\mathcal{E}(n))$
$\displaystyle=k(n,n)+\kappa(n)-k(n,\mathbf{n})[k(\mathbf{n^{1}},\mathbf{n^{1}})+\textrm{diag}(\kappa(\mathbf{n^{1}}))]^{-1}k(\mathbf{n^{1}},n)$
(48)

Note that this differs only slightly from the emulator constructed in section
5.3, with the main difference being we now attribute uncertainty in the loss
values as an inherent behaviour of our emulator and not in the procedure to
obtain these loss values. As a result $\kappa(n)$ does not decrease as the
multiplicity of elements of $\mathbf{n}$ increases, which represents a major
disadvantage to the uncertainty representation in section 5.3.

One may then be sceptical of the benefit of duplicating design points for this
method, and whilst it is possible to use this method without duplication (i.e
$\mathbf{n}=\mathbf{n^{1}}$), the consequence of this would be that we are
heavily reliant on a singular sample to locate the minimum which could be
misleading. Averaging various samples at the same design point mitigates this
potential problem, as does replacing $d^{-}$ with
$\mu^{-}=\min_{i}\\{\mu(\mathbf{n^{1}}_{i})\\}$ as detailed in Brochu et al.
(2010). Taking the median of samples instead of a weighted mean is more
appropriate here as we are not seeking to accurately approximate an
expectation, instead we only wish to avoid extreme samples misleading our
search for the minimum.

####  S6.2 Proof of Theorem S6

######  Theorem S6.

If $\ell(n)$, $\boldsymbol{\sigma}$, and $m(n,\Theta)$ are almost surely
bounded and $d_{i}\sim N\left\\{l(n_{i}),\sigma_{i}^{2}\right\\}$ then for
every $n\in\\{1,\ldots,N\\}$, as the multiplicity of $n$ in $\mathbf{n}$ tends
to $\infty$ we have $\mu(n)\longrightarrow\ell(n)$ and $\Psi(n)\to 0$ almost
surely with respect to variation in $\mathbf{d}$.

###### Proof.

Assume W.L.O.G that $(\mathbf{n}^{1})_{1}=n$. Since $\boldsymbol{\sigma}$ is
bounded, we have (from equation 42)
$\textrm{var}\left((\mathbf{d}^{1})_{1}\right)=(\boldsymbol{\sigma^{1}})^{2}_{1}\to
0$, so $(\mathbf{d}^{1})_{1}\longrightarrow\ell(n)$ almost surely. We now
prove that
$k(n,\mathbf{n^{1}})[k(\mathbf{n^{1}},\mathbf{n^{1}})+\textrm{diag}((\boldsymbol{\sigma^{1}})^{2})]^{-1}=(1,0,\ldots,0)$
when $(\boldsymbol{\sigma^{1}})_{1}=0$. Now:

$\displaystyle
k(n,\mathbf{n^{1}})[k(\mathbf{n^{1}},\mathbf{n^{1}})+\textrm{diag}((\boldsymbol{\sigma^{1}})^{2})]^{-1}$
$\displaystyle=(1,0,\ldots,0)$ $\displaystyle\Leftrightarrow\hskip
20.0ptk(n,\mathbf{n^{1}})$
$\displaystyle=(1,0,\ldots,0)*k(\mathbf{n^{1}},\mathbf{n^{1}}+\textrm{diag}((\boldsymbol{\sigma^{1}})^{2}))$
(49)

and
$k(n,\mathbf{n^{1}})=(1,0,\ldots,0)*k(\mathbf{n^{1}},\mathbf{n^{1}}+\textrm{diag}((\boldsymbol{\sigma^{1}})^{2})_{-1})$
is true by definition as the first row of
$k(\mathbf{n^{1}},\mathbf{n^{1}})+\textrm{diag}((\boldsymbol{\sigma^{1}})^{2})$
is $k(n,\mathbf{n^{1}})$. Therefore,

$\mu(n)=m(n,\Theta)+(1,0,\ldots,0)(\mathbf{d^{1}}-m(\mathbf{n^{1}},\Theta))=m(n,\Theta)+d(n)-m(n,\Theta)=d(n)=\ell(n)$
(50)

almost surely, and

$\Psi(n)=k(n,n)-(1,0,\ldots,0)k(\mathbf{n^{1}},n)=k(n,n)-k(n,n)=0$ (51)

in the limit.

∎

####  S6.3 Proof of Theorem S7

######  Theorem S7.

Given the conditions of Theorem S6, for every $n\in\\{1,\ldots,N\\}$, as the
multiplicity of $n$ in $\mathbf{n}$ tends to $\infty$,

$EI(n)\longrightarrow 0$

almost surely with respect to randomness in $\mathbf{d}$

###### Proof.

From Theorem S6 we have that $\mu(n)\longrightarrow\ell(n)<\infty$ and
$d^{1}_{i}\longrightarrow\ell(n^{1}_{i})$, so therefore in the limit we can
state $pr_{\mathbf{d}}(-\infty<d^{-}-\mu(n)\leq 0)=1$. Indeed, let $j$ be the
index such that $n^{1}_{j}=n$. If in the limit $d^{-}>\mu(n)=d(n)$ then this
implies that $d^{1}_{j}<\min_{i}\\{d^{1}_{i}\\}$ which is a contradiction.
Also note from Theorem S6 that $\Psi(n)\longrightarrow 0$ and that
$\Phi\left(\cdot\right)\in(0,1)$,
$\phi\left(\cdot\right)\in(0,(2\pi)^{-1/2}]$. As a result the following two
scenarios have joint probability 1:

  * •

$d^{-}-\mu(n)=0$ in the limit: As $\Phi\left(\cdot\right)$,
$\phi\left(\cdot\right)$ are bounded and $\Psi(n)=0$ in the limit, we also
have $EI(n)=0$ in the limit.

  * •

$\infty<d^{-}-\mu(n)<0$ in the limit: As $\Psi(n)=0$ in the limit,
$\Phi\left(\frac{d^{-}-\mu(n)}{\sqrt{\Psi(n)}}\right)=0$ in the limit. As
$\phi\left(\cdot\right)$ is bounded we have that $EI(n)=0$ in the limit.

which proves the corollary.

∎

####  S6.4 Proof of Theorem S8

######  Theorem S8.

If $\ell(n)$, $\boldsymbol{\sigma}$, and $m(n,\Theta)$ are almost surely
bounded and $d_{i}\sim N\left\\{l(n_{i}),\sigma_{i}^{2}\right\\}$ then under
algorithm 4 with $\tau=0$, the value $\mu(\tilde{n})$ converges almost surely
to $\ell(\tilde{n})$ for every $\tilde{n}\in\\{1,\ldots,N\\}$.

###### Proof.

Our overall argument is to show that algorithm 4 leads to the multiplicity of
$\tilde{n}$ in $\mathbf{n}$ tending to infinity, from which the result follows
from Theorem S6.

To do this, we begin with the following two lemmas, the second of which
describes the limiting behaviour of $EI(n)$ according to how often $n$ occurs
in $\mathbf{n}$: namely that if the multiplicity of $n$ in $\mathbf{n}$
diverges, the value of $EI(n)$ converges to 0; otherwize, it remains positive.
We introduce the index $EI_{\mathbf{n}}(n)$ to indicate the dependence of
$EI(n)$ on $\mathbf{n}$ and assume that the function $\ell(n)$ is fixed. For a
multiset $\mathbf{n}_{i}$, we denote $\textrm{mult}_{\mathbf{n_{i}}}(n)$ as
the multiplicity of $n$ in $\mathbf{n_{i}}$.

######  Lemma S2.

Suppose $m\times m$ matrix $A$ is symmetric. Denote by $I^{1}$ the $m\times m$
matrix with $I^{1}_{ij}=1_{i=j=1}$. Let $x$ be a vector of length $m$ and
denote by $A_{x}$ the matrix $A$ with its top row replaced by $x$. Then for
$p$ in any interval containing 0 on which $A+pI^{1}$ is invertible we have

$\frac{\partial}{\partial
p}\left(x^{T}(A+pI^{1})^{-1}x\right)=-\frac{\left|A_{x}\right|^{2}}{\left|A+pI^{1}\right|^{2}}$
(52)

###### Proof.

If $M(p)$ is invertible in a neighbourhood of $p$ we have $\frac{\partial
M^{-1}}{\partial p}=-M^{-1}\frac{\partial M}{\partial p}M^{-1}$, and if $M$ is
symmetric with dimensions $m\times m$ and first row $M_{1}$, then
$MI^{1}M=M_{1}M_{1}^{T}$. Since $(A+pI)$ and $A$ differ only in the top row,
we have $\textrm{adj}(A+pI)_{1}=\textrm{adj}(A)_{1}$, where
$\textrm{adj}(\cdot)$ indicates the adjugate matrix and $\cdot_{1}$ the top
row. We now have

$\displaystyle\frac{\partial}{\partial p}\left(x^{T}(A+pI^{1})^{-1}x\right)$
$\displaystyle=-x^{T}(A+pI^{1})^{-1}\frac{\partial(A+pI^{1})}{\partial
p}(A+pI^{1})^{-1}x$ $\displaystyle=-x^{T}(A+pI^{1})^{-1}I^{1}(A+pI^{1})^{-1}x$
$\displaystyle=\frac{x^{T}\textrm{adj}(A+pI^{1})I^{1}\textrm{adj}(A+pI^{1})x}{|A+pI^{1}|^{2}}$
$\displaystyle=-\frac{x^{T}\textrm{adj}(A+pI^{1})_{1}\textrm{adj}(A+pI^{1})_{1}^{T}x}{|A+pI^{1}|^{2}}$
$\displaystyle=-\frac{x^{T}\textrm{adj}(A)_{1}\textrm{adj}(A)_{1}^{T}x}{|A+pI^{1}|^{2}}$
$\displaystyle=-\frac{|A_{x}|^{2}}{|A+pI^{1}|^{2}}$

as required. ∎

######  Lemma S3.

Let $S_{1}$ and $S_{2}$ be disjoint subsets of $[N]=\\{1,\dots,N\\}$ with
$S_{1}\cup S_{2}=[N]$. For a multiset $\mathbf{n}$ denote

$\displaystyle q_{1}(\mathbf{n})$ $\displaystyle=\max_{n\in
S_{1}}\textrm{mult}_{\mathbf{n}}(n)$ $\displaystyle q_{2}(\mathbf{n})$
$\displaystyle=\min_{n\in S_{2}}\textrm{mult}_{\mathbf{n}}(n)$

Suppose we have infinite sequences $\mathbf{n}$, $\mathbf{d}$,
$\boldsymbol{\sigma}$, where $\mathbf{d}\sim
N(l(\mathbf{n}),\boldsymbol{\sigma}^{2})$ and $\boldsymbol{\sigma^{2}}$ is
upper-bounded, and let $\mathbf{n_{i}}$, $\mathbf{d_{i}}$,
$\boldsymbol{\sigma_{i}}$ denote the (multiset) first $i$ elements of each
sequence. Suppose that $q_{1}(\mathbf{n_{i}})\leq m_{1}$ for all $i$ and
$q_{2}(\mathbf{n_{i}})\to\infty$, and the set
$\left\\{k_{2}(n,\Theta_{i})=k_{2}\left(n,\Theta(\mathbf{n_{i}},\mathbf{d_{i}},\boldsymbol{\sigma_{i}})\right):n\in\\{1,\dots,N\\},i\in\mathbb{N}\right\\}$
is almost surely asymptotically bounded. Then for sufficiently large
$\sigma_{u}$:

$\lim\sup_{i\to\infty}EI_{\mathbf{n}_{i}}(n)=\begin{cases}e_{n}>0&\textrm{if
}n\in S_{1}\\\ 0&\textrm{if }n\in S_{2}\end{cases}$ (54)

almost surely.

###### Proof.

We will in fact show that even $\lim\inf EI_{\mathbf{n_{i}}}(n)>0$ for $n\in
S_{1}$, but $\lim\sup$ will suffice for our purposes. We note that

$EI_{\mathbf{n}_{i}}(n)>0\Leftrightarrow\sqrt{\Psi_{\mathbf{n}_{i}}(n)}\phi\left(\frac{d^{-}_{\mathbf{n}_{i}}-\mu_{\mathbf{n}_{i}}(n)}{\sqrt{\Psi_{\mathbf{n}_{i}}(n)}}\right)>(\mu_{\mathbf{n}_{i}}(n)-d^{-}_{\mathbf{n}_{i}})\Phi\left(\frac{d^{-}_{\mathbf{n}_{i}}-\mu_{\mathbf{n}_{i}}(n)}{\sqrt{\Psi_{\mathbf{n}_{i}}(n)}}\right)$
(55)

We will show that for all $n$, we have

$P\left(-\infty<\lim\inf_{i\to\infty}\left(d_{\mathbf{n_{i}}}^{-}-\mu_{\mathbf{n_{i}}}(n)\right)\right)=1$
(56)

By the argument in theorem S6 and corollary S7 we have for $n\in S_{2}$ that
$\lim_{i\to\infty}\Psi_{\mathbf{n_{i}}}(n)=0$, from which both sides of 55
converge to 0. For $n\in S_{1}$ we will show
$\lim_{i\to\infty}\Psi_{\mathbf{n_{i}}}(n)>0$, in which case we may define

$z_{\mathbf{n}_{i}}(n)=\frac{\mu_{\mathbf{n}_{i}}(n)-d^{-}_{\mathbf{n}_{i}}}{\sqrt{\Psi_{\mathbf{n}_{i}}(n)}}$
(57)

from which inequality 55 reduces to

$\phi(z_{\mathbf{n}_{i}}(n))>z_{\mathbf{n}_{i}}(n)\Phi(-z_{\mathbf{n}_{i}}(n))$
(58)

which holds for all $0\leq z_{\mathbf{n}_{i}}(n)<-\infty$. Since
$z_{\mathbf{n}_{i}}(n)$ is asymptotically bounded between positive values, the
result follows.

Beginning with $d_{\mathbf{n_{i}}}^{-}$, we note that $d_{\mathbf{n_{i}}}^{-}$
is the minimum of

  1. 1.

Values of $\mathbf{d^{1}_{i}}$ corresponding to values of $\mathbf{n^{1}_{i}}$
in $S_{1}$; and

  2. 2.

Values of $\mathbf{d^{1}_{i}}$ corresponding to values of $\mathbf{n^{1}_{i}}$
in $S_{2}$

For sufficiently large $s$, the sequence $\\{n_{j}=(\mathbf{n})_{j}:j>s\\}$
never contains any $n\in S_{1}$ again; hence, the minimum of item 1 is
determined after finitely many $i$ and its limit is finite. Since
$\boldsymbol{\sigma}$ is upper-bounded, all values of $\mathbf{d^{1}_{i}}$ in
item 2 converge to finite values in $\\{\ell(n):n\in S_{2}\\}$ almost surely.
Hence $d_{\mathbf{n_{i}}}^{-}$ converges almost surely to a finite value.

Since $\lim\sup_{i\to\infty}$ and $\lim\inf_{i\to\infty}$ of
$m\left(n;\Theta\left(\mathbf{n_{i}},\mathbf{d_{i}},\boldsymbol{\sigma_{i}}\right)\right)$
are almost surely finite, all terms in $\mu(n)$ are asymptotically finite,
from which equation 56 follows.

It remains to consider $\Psi_{\mathbf{n}_{i}}(n)$ for $n\in S_{1}$. Firstly
take $n\in\mathbf{n^{1}}$ and suppose W.L.O.G that
$\boldsymbol{n^{1}_{i}}_{1}=n$. Since $n\in S_{1}$ we have
$\lim_{i\to\infty}\textrm{mult}_{\mathbf{n_{i}}}(n)>0$ so
$\lim_{i\to\infty}(\boldsymbol{\sigma^{1}_{i}})_{1}$ exists and is positive.
Denoting $\boldsymbol{\sigma^{\prime}}$ as $\boldsymbol{\sigma^{1}_{i}}$ with
0 substituted for the first element, we have

$\displaystyle\frac{\partial}{\partial(\boldsymbol{\sigma^{1}_{i}})^{2}_{1}}\Psi_{\mathbf{n_{i}}}(n)$
$\displaystyle=\frac{\partial}{\partial(\boldsymbol{\sigma^{1}_{i}})^{2}_{1}}\left(k(n,n)-k(n,\mathbf{n^{1}_{i}})[k(\mathbf{n^{1}_{i}},\mathbf{n^{1}_{i}})+\textrm{diag}((\boldsymbol{\sigma^{1}_{i}})^{2})]^{-1}k(\mathbf{n^{1}_{i}},n)\right)$
$\displaystyle=\frac{\left|k(\mathbf{n^{1}_{i}},\mathbf{n^{1}_{i}})+\textrm{diag}((\boldsymbol{\sigma^{\prime}})^{2})\right|^{2}}{\left|k(\mathbf{n^{1}_{i}},\mathbf{n^{1}_{i}})+\textrm{diag}((\boldsymbol{\sigma^{1}_{i}})^{2})\right|^{2}}$
$\displaystyle>0$

by lemma S2; hence $\Psi_{\mathbf{n_{i}}}(n)$, considered as a function of
$(\boldsymbol{\sigma^{1}_{i}})^{2}_{1}$, is increasing. Given that
$\lim_{i\to\infty}(\boldsymbol{\sigma^{1}_{i}})_{j}$ is 0 for
$(\mathbf{n^{1}_{i}})_{j}\in S_{2}$ and is positive for
$(\mathbf{n^{1}_{i}})_{j}\in S_{1}$, we conclude that
$\lim_{i\to\infty}\Psi_{\mathbf{n_{i}}}(n)$ is positive when $n\in S_{1}$ and
$n\in\mathbf{n^{1}}$.

If $n\notin\mathbf{n^{1}}$, so $n$ never occurs in any $\mathbf{n_{i}}$, then
we firstly note that since $k(n,n)<k(n,m)$ for any $m\neq n$, we have:

$k(n,n)-k(n,\mathbf{n^{1}_{i}})[k(\mathbf{n^{1}_{i}},\mathbf{n^{1}_{i}}))]^{-1}k(\mathbf{n^{1}_{i}},n)>0$
(59)

This omits the term $\textrm{diag}((\boldsymbol{\sigma^{1}_{i}})^{2}$ from the
expression for $\Psi_{\mathbf{n_{i}}}(n)$. However, if we denote
$k^{\prime}_{j}$ the matrix
$k(\mathbf{n^{1}_{i}},\mathbf{n^{1}_{i}})+\textrm{diag}((\boldsymbol{\sigma^{1}_{i}})^{2})$
with the $j$th row replaced by $k(n,\mathbf{n^{1}_{i}})$, we have from lemma
S2:

$\frac{\partial}{\partial(\boldsymbol{\sigma^{1}_{i}})^{2}_{j}}\Psi_{\mathbf{n_{i}}}(n)=\frac{\left|k^{\prime}_{j}\right|^{2}}{\left|k(\mathbf{n^{1}_{i}},\mathbf{n^{1}_{i}})+\textrm{diag}((\boldsymbol{\sigma^{1}_{i}})^{2})\right|^{2}}>0$

for any element $(\boldsymbol{\sigma^{1}_{i}})^{2}_{j}$ of
$(\boldsymbol{\sigma^{1}_{i}})^{2}$; hence $\Psi_{\mathbf{n_{i}}}(n)$ is
increasing in any such element and its positivity follows. This completes the
proof of the lemma.

∎

Now suppose that some $n\in\\{1,\dots,N\\}$ occurs only finitely often in
$\mathbf{n^{1}}$. Then there must be some largest set $S_{1}$ of such $n$,
with complement $S_{2}=\\{1,\dots,N\\}\setminus S_{1}$. Since every element in
$S_{1}$ occurs in $\mathbf{n^{1}}$ with finite multiplicity there must be some
$j$ such that no $n\in S_{1}$ occurs amongst the values
$\\{(\mathbf{n^{1}})_{j+1},(\mathbf{n^{1}})_{j+2},\dots\\}$. But from lemma
S3, there will almost surely eventually be some $J>j$ for which some value in
$\\{EI_{\mathbf{n_{J}}}(n):n\in S_{1}\\}$ exceeds all values in
$\\{EI_{\mathbf{n_{J}}}(n):n\in S_{1}\\}$, and hence
$(\mathbf{n^{1}})_{J+1}\in S_{1}$ (as long as $\tau$ is sufficiently small),
contradicting the choice of $j$. So the event that an $n\in\\{1,\dots,N\\}$
occurs in $\mathbf{n^{1}}$ with finite multiplicity has probability 0. This
completes the proof.

∎

###  S7 Optimal holdout set size arising from a simulated example

In this section, we analyse the dynamics of a roughly realistic, binary
outcome system, subject to predictions from different families of risk models.
Our main aim is to demonstrate the natural emergence of an optimal holdout set
size from a reasonable setting.

We generated datasets with a population size $N=5000$ with seven standard
normally distributed covariates and outcomes $Y$ under a ground-truth logistic
model, either with interaction terms (i.e., non-linear) or without (linear).
We considered risk scores $\rho$ derived from either logistic regression
models (not including interaction terms) or random forests. We designated cost
functions $C_{1}$, $C_{2}$ to have value 0 for true-negatives, 0.5 for false-
or true- positives, and 1 for false-negatives.

Figure 7 shows simulation results using either linear or logistic prediction
models and linear or non-linear underlying models for $Y\mid X$. We can
observe that an optimal holdout set size can arise naturally from standard
predictive models, since empirical $k_{2}$ curves for both a random forest and
logistic regression satisfy assumptions 2 and 4 in the main manuscript. The
optimal holdout set size occurs at a value $n$ smaller than that at which
$k_{2}(n)$ is nearly ‘flat’, indicating that unnecessarily large training sets
are suboptimal. However, since $\ell(n)$ rises only linearly as $n$ increases,
it is generally less costly to slightly overestimate rather than underestimate
the optimal holdout set size. Finally, the rightmost panels illustrate that
the optimal holdout set size is not necessarily smaller for a more accurate
model: the random forest model (non-lin $\rho$) in the non-linear underlying
case (right panels) leads to uniformly lower expected costs $k_{2}(n)$ at all
potential holdout set sizes, although the optimal holdout set size is larger.

Figure 7: Examples of cost functions as per Theorem 1 arising naturally from a
basic risk score, with varying underlying model (und.), risk score type
($\rho$) and one point-wise standard deviation (shaded regions). The
contributions of terms $k_{1}n$ to $\ell(n)$ depend only on the underlying
model and are the same in each column. OHS: optimal holdout set size

###  S8 Estimation of parameters for optimal holdout size in ASPRE

####  S8.1 Estimation of $N$

As per our assumptions, we presume we are refitting ASPRE to use in a
population of 5 million individuals, from which we have approximately 80,000
new pregnancies per year. The incidence of pregnancy per year is now

$\frac{8\times 10^{4}}{5\times 10^{6}}=\frac{1}{125}$ (60)

so we have

$\displaystyle N$ $\displaystyle\approx 5\times 8\times 10^{4}$
$\displaystyle=400000$ (61)

with standard error

$\displaystyle SE(N)$ $\displaystyle\approx 5\sqrt{5\times
10^{6}\times\frac{1}{125}\left(1-\frac{1}{125}\right)}$ $\displaystyle\approx
1500$ (62)

####  S8.2 Estimation of $k_{1}$ and $k_{2}$

We assume $\pi=10\%\approx 2707/25797$, the proportion of individuals assigned
to the treatment group in Rolnik et al. (2017a) due to having an estimated
risk of PRE $>1\%$.

To estimate $k_{1}$, we considered the study reported in O’Gorman et al.
(2017) assessing sensitivity and specificity of NICE and ACOG guidelines in
assessing PRE risk. In this study, 8775 indivduals were assessed, amongst
which 239 developed PRE, for an overall incidence of $239/8875\approx 0.027$.
We estimated the performance of a ‘baseline’ estimator of PRE risk (that is,
in the absence of any ASPRE score) by linearly interpolating the points
corresponding to ‘ACOG aspirin’, ‘NICE’ and ‘ACOG’ on ROC curves in Figure 1.
On this basis, a baseline estimator identifying the 10% of individuals at
highest PRE risk (approximately 800) would correspond to the point $(x,y)$ on
the interpolated ROC curve with

$239x+(8775-239)y=0.1\times 8875$ (63)

which occurs at roughly a 20% detection (true positive) rate and a 10% false
positive rate, close to that of the NICE guidelines.

Since few women in the study were treated with aspirin, we assume that PRE
rates in the highest-10% and lowest-10% risk groups assessed by baseline risk
(NICE) are untreated risk (that is, if not treated with aspirin). At the
inferred true and false positive rates, we would expect that amongst the 10%
of women designated highest-risk by the NICE guidelines, we have a PRE rate of

$\displaystyle\pi_{1}$
$\displaystyle\approx\frac{\textrm{TPR}\times\textrm{(Num. PRE)}}{\textrm{Num.
positive}}$ $\displaystyle=\frac{0.2\times 239}{0.1\times 8875}$
$\displaystyle\approx 0.054$ (64)

and amongst the 90% designated lower risk, a PRE rate of

$\displaystyle\pi_{0}$
$\displaystyle\approx\frac{(1-\textrm{TPR})\times\textrm{(Num.
PRE)}}{\textrm{Num. negative}}$ $\displaystyle=\frac{0.8\times 239}{0.9\times
8875}$ $\displaystyle\approx 0.024$ (65)

Given that true positive rates of the NICE guidelines are computed as a
fraction with denominator 239, we presume standard errors of $\pi_{1}$ and
$\pi_{0}$ of

$\displaystyle SE(\pi_{1})$
$\displaystyle\approx\frac{\sqrt{\pi_{1}(1-\pi_{1})}}{8875\times 0.1}$
$\displaystyle\approx 0.0076$ $\displaystyle SE(\pi_{0})$
$\displaystyle\approx\frac{\sqrt{\pi_{0}(1-\pi_{0})}}{8875\times 0.9}$
$\displaystyle\approx 0.0017$ (66)

Now, treating errors in $\pi_{0}$, $\pi_{1}$ and $\alpha$ as pairwise
independent

$\displaystyle k_{1}$ $\displaystyle=\pi_{0}(1-\pi)+\pi_{1}\pi\alpha$
$\displaystyle\approx 0.0235$ $\displaystyle SE(k_{1})$
$\displaystyle=SE\left(\pi_{0}(1-\pi)+\pi_{1}\pi\alpha\right)$
$\displaystyle\approx 0.0016$ (67)

We estimate the population prevalence $\pi_{PRE}$ of untreated PE as the
frequency observed in the original ASPRE data:

$\pi_{PRE}=\frac{1426}{57974}\approx 2.4\%$ (68)

Note that, although this is approximately equal to $\pi_{0}$, they are
different quantities: $\pi_{0}$ is the population frequency of PRE amongst
individuals at the lowest 90% risk by NICE guidelines.

Denoting $\pi_{1}(n)$ as the untreated risk of PRE in the top 10% of
individuals according to an ASPRE score trained on $n$ individuals (and
$\pi_{0}(n)$ correspondingly), we note that it is equal to the sensitivity (or
TPR) of the risk score at the level where proportion $\pi$ of individuals are
designated high-risk. Thus for any training set size $n$

$\pi_{0}(n)=\frac{\pi_{PRE}-\pi\pi_{1}(n)}{1-\pi}$ (69)

so the average cost to an individual in the intervention set may be expressed
in terms of $\pi_{1}(n)$:

$\displaystyle k_{2}(n)$ $\displaystyle=\pi_{0}(n)(1-\pi)+\pi_{1}(n)\pi\alpha$
$\displaystyle=\pi_{PRE}-\pi\pi_{1}(n)(1-\alpha)$ (70)

####  S8.3 Implementation

We implemented the complete ASPRE model as described in Rolnik et al. (2017b).
We simulated a population of individuals with a similar distribution of ASPRE
model covariates. We computed the ASPRE scores for our simulated individuals,
and found a linear transformation of these scores such that, should the scores
exactly specify the probability of PRE, the expected population prevalence and
sensitivity of the score would match those reported in Rolnik et al. (2017a):
prevalence $\pi_{PRE}$, and sensitivity amongst 10% highest scores: 12.3%. We
then simulated PRE incidence according to these transformed scores.

We found that a generalised linear model with logistic link performed almost
as well as the ASPRE score on our simulated data, so we used this model type
to estimate the learning curve in the interests of simplicity.

To choose values $\mathbf{n}$ and $\mathbf{k_{2}}/\mathbf{d}$, we initially
chose a set $\mathbf{n}$ of 20 random values from $[500,30000]$. For each size
$n$ in $\mathbf{n}$, we took a random sample of our data of size $n$, fitted a
logistic model to that sample, and estimated corresponding expected costs per
individual $\mathbf{k_{2}}$ as above. We fitted values
$\theta=\theta(\mathbf{n},\mathbf{k_{2}})=(a,b,c)$ parametrising $k_{2}$ as
the maximum-likelihood estimator of $\theta$ under the model

$(\mathbf{k_{2}})_{i}\sim
N\left(k_{2}((\mathbf{n})_{i},\theta),\sigma^{2}\right)\sim
N(a(\mathbf{n})_{i}^{-b}+c,\sigma^{2})$ (71)

for a fixed values $\sigma$, noting that the estimate of $\theta$ is
independent of $\sigma$. For the parametric algorithm, we then set all values
of $\boldsymbol{\sigma}$ to the same value, chosen empirically as the sample
variance of

$\mathbf{k_{2}}-k_{2}(\mathbf{n},\theta(\mathbf{n},\mathbf{k_{2}}))$ (72)

For the emulation algorithm, we set values $\mathbf{d}$ as

$\mathbf{d}_{i}=k_{1}(\mathbf{n})_{i}+(\mathbf{k_{2}})_{i}(N-(\mathbf{n})_{i})$
(73)

transforming values $\boldsymbol{\sigma}$ correspondingly for use in the
emulation algorithm. We then sequentially chose 100 additional values
$\mathbf{n}$ using both algorithm 3 and 4, setting $\boldsymbol{\sigma}$ as
the same value found in 72. After choosing the 120 values of $\mathbf{n}$
using algorithm 3, we re-estimated $\mathbf{k_{2}}/\mathbf{d}$ for each of
these values before estimating the OHS and confidence interval to avoid any
potential regression-to-the mean effects from choosing next-values-of-$n$ so
as to minimise estimated confidence interval width.

Our complete pipeline is available at
https://github.com/jamesliley/OptHoldoutSize_pipelines, and a comprehensive
vignette is included in our R package OptHoldoutSize on CRAN and at
https://github.com/jamesliley/OptHoldoutSize.

###  S9 Supplementary discussion

In this section we discuss the use and consequences of several of our methods
and results in greater depth.

Theorem 1 establishes that straightforward conditions on the system to be
modelled lead to an optimal holdout set size which is straightforward to find,
and indeed that the cost function is convex (in a discrete sense). An obvious
limitation is the requirement for convexity of $k_{2}$; risk score learning
curves, on which $k_{2}$ depends, are not necessarily convex for complex
models (Viering and Loog, 2021). In practice, the cost function $\ell(n)$ is
still generally well-behaved even if convexity of $k_{2}$ is violated, and may
be approximated using our emulation algorithm (as demonstrated in
Supplementary Figure 11(b) and Figure 4 in the main manuscript).

The use of a Gaussian process emulator to approximate the true loss function
enables an automatic selection of the optimal holdout set size under fewer
assumptions, although the efficiency of this method heavily depends on the
quality of our emulator. Various extensions of the emulator may improve our
surrogate of the loss function, for example specifying priors on the
parameters $\theta,\sigma^{2}_{u},\zeta$ and using the likelihood provided by
the Gaussian process to marginalize out these parameters. An explicit approach
is given in Andrianakis and Challenor (2011), but under linearity assumptions
which do not hold in our case, so analytic tractability would be lost. If we
were able to cheaply estimate the derivative of the cost function at design
points, this could be incorporated into our emulator (Killeya, 2004), enabling
greater posterior accuracy around these points. Direct estimation of gradients
from only estimates of $\ell(n)$ usually requires double the number of
evaluations as estimation of $\ell(n)$ values, and so has the potential to
become a more costly procedure than the method presented in section 5.3.

Our simulations and theoretical findings illuminate several non-obvious
properties of the optimal holdout set size. Firstly, from section 6, we note
that the optimal holdout set size is often fairly small, and training a risk
score until improvement in $k_{2}$ is negligible is counterproductive. More
accurate risk scores do not necessarily lead to lower optimal holdout set
size, but generally lead to lower loss. Given corollary 1, it is generally
better to err on a higher side of the optimal holdout set size, since cost
increases at most linearly.

###  S10 Supplementary figures

Figure 8: Dependence of optimal holdout set size on parameters of estimated
learning curve ($a,b,c$, with $k_{2}(n;a,b,c)=an^{-b}+c$), cost in
intervention set $k_{1}$, and total number of samples $N$. Figures show change
in optimal holdout set size $n_{*}$ while varying one parameter and holding
others constant at $(a,b,c)=\left(\frac{3}{2},\frac{3}{2},\frac{1}{4}\right)$,
$k_{1}=\frac{4}{5}$, $N=10^{4}$.

Figure 9: Dependence of minimum total cost on parameters of estimated learning
curve ($a,b,c$, with $k_{2}(n;a,b,c)=an^{-b}+c$), cost in intervention set
$k_{1}$, and total number of samples $N$. Figures show change in minimal cost
$\ell(n_{*})$ while varying one parameter and holding others constant at
$(a,b,c)=\left(\frac{3}{2},\frac{3}{2},\frac{1}{4}\right)$,
$k_{1}=\frac{4}{5}$, $N=10^{4}$.

(a)

(b)

Figure 10:  Left panel shows emulator constructed using three $k_{2}()$ values
(see pipelines). Function $m(n,\Theta)$ is constructed using $\theta$ derived
from these three $k_{2}()$ estimates. Note reduced pointwise posterior
variance at sample points. Rightmost panel shows expected improvement plot for
the emulator constructed in panel10(a). Note local minima at existing sample
points.

(a) Versions of $k_{2}$

(b) Costs associated with $k_{2}$

(c) Distr/meds of est. optimal holdout set size

Figure 11: Parametric and emulation algorithms with parametric assumptions
satisfied or unsatisfied. OHS: optimal holdout set size

### References

  * Y. Andrianakis and P. Challenor (2011) Parameter estimation for Gaussian process emulators.  Technical report Technical report, Managing Uncertainty in Complex Models.  Cited by: §S9.
  * R. G. Bower, M. Goldstein, and I. Vernon (2010) Galaxy formation: a Bayesian uncertainty analysis.  Bayesian analysis 5 (4),  pp. 619–669.  Cited by: §S6.1.