# Monitoring machine learning (ML)-based risk prediction algorithms in the
presence of confounding medical interventions

Jean Feng1, Alexej Gossmann2, Gene Pennello2, Nicholas Petrick2,
Berkman Sahiner2, Romain Pirracchio1

( 1University of California, San Francisco
2U.S. Food and Drug Administration, Center for Devices and Radiological Health
)

###### Abstract

Performance monitoring of machine learning (ML)-based risk prediction models
in healthcare is complicated by the issue of confounding medical interventions
(CMI): when an algorithm predicts a patient to be at high risk for an adverse
event, clinicians are more likely to administer prophylactic treatment and
alter the very target that the algorithm aims to predict. A simple approach is
to ignore CMI and monitor only the untreated patients, whose outcomes remain
unaltered. In general, ignoring CMI may inflate Type I error because (i)
untreated patients disproportionally represent those with low predicted risk
and (ii) evolution in both the model and clinician trust in the model can
induce complex dependencies that violate standard assumptions. Nevertheless,
we show that valid inference is still possible if one monitors conditional
performance and if either conditional exchangeability or time-constant
selection bias hold. Specifically, we develop a new score-based cumulative sum
(CUSUM) monitoring procedure with dynamic control limits. Through simulations,
we demonstrate the benefits of combining model updating with monitoring and
investigate how over-trust in a prediction model may delay detection of
performance deterioration. Finally, we illustrate how these monitoring methods
can be used to detect calibration decay of an ML-based risk calculator for
postoperative nausea and vomiting during the COVID-19 pandemic.

##  1 Introduction

After a machine learning (ML)-based system is deployed in clinical practice,
real-world monitoring of the algorithm for potential performance degradation
is necessary for mitigating risk and is an important aspect of good machine
learning practice (GMLP) (U.S. Food and Drug Administration and Health Canada,
2021). Various methods for performance monitoring are available (Feng et al.,
2022b; Kahn et al., 1996), such as those from statistical process control. All
of these procedures assume an ideal data setting in which the prediction
target is observed. However, the data available for monitoring a ML-based risk
prediction algorithm are often subject to confounding medical interventions
(CMI), because predictions from the algorithm can alter the very outcome that
it aims to predict (Paxton et al., 2013; Lenert et al., 2019; Perdomo et al.,
2020).

As a motivating example, consider the Targeted Real-time Early Warning System
(TREWS) sepsis risk prediction algorithm (Adams et al., 2022), which estimates
the probability of a patient developing septic shock if they only receive
standard of care (SOC) and no additional interventions. This algorithm was
recently shown to reduce in-hospital mortality and organ failure rates. The
alert system increased the propensity and speed at which clinicians submitted
their first antibiotic order (Henry et al., 2022). Moreover, the authors found
that clinicians’ likelihood to interact with the TREWS alert system depended
on their previous interactions, and hypothesize that clinician trust will
continue to evolve with increased exposure to ML-based systems. Monitoring the
performance of the TREWS algorithm is especially important because (i) it
depends on electronic health record (EHR) data, which is prone to distribution
shifts, and (ii) overuse of antibiotics has major negative consequences. To
evaluate the model’s predictions for patient outcomes under SOC, a potential
solution is to restrict our attention to only those patients who received SOC,
as the counterfactual outcomes for the patients who received antibiotics are
unknown. However, in the likely scenario where high-risk patients are
preferentially selected to receive an intervention, procedures that monitor
marginal performance measures (e.g., misclassification rate or AUC) without
adjusting for CMI are biased, due to dependent censoring of the outcome
(Figure 9 of the Appendix). This can lead to inflated false alarm rates and/or
unnecessarily long detection delays.

In the offline setting, one can try to address the mismatch between the “SOC-
only” and general target population by reweighting the data by the inverse of
their treatment propensities. In the sequential setting, the closest work
along these lines is Sun et al. (2014), which combines inverse censoring
weights with a cumulative sum (CUSUM) algorithm to monitor survival outcomes.
However, proper error rate control is contingent on knowing the exact weights,
which is unlikely to hold in our setting. Moreover, we may not even be able to
accurately estimate the propensity weights, because it is difficult to
anticipate how clinician trust in an algorithm will evolve over time. Finally,
even with access to exact weights, near-violations of the positivity condition
(i.e. weights close to zero) can drastically slow down asymptotic convergence
of these procedures, thereby inflating the false alarm rate.

Given the difficulties of monitoring marginal performance measures using a
propensity-based approach, we propose to monitor conditional measures of
performance, which we show are independent of treatment propensities under one
of two ignorability conditions. The first is the well-known conditional
exchangeability assumption (Rubin, 1976). The second is the assumption of
time-constant selection bias, which has not yet been discussed in the
literature to our knowledge. If either condition holds, we can ignore CMI, use
a “standard” monitoring procedure to analyze the SOC-only data, and avoid
estimating treatment propensities altogether. We pay special attention to the
model calibration, a popular conditional performance measure and one of the
most common types of performance deterioration in practice (Hickey et al.,
2013; Davis et al., 2017).

In addition, we address two challenges in analyzing performance of ML-based
risk prediction models. First, the sequence of predictors among the SOC-only
patients can be highly nonstationary, because the ML algorithm and/or the
clinician’s interactions with the ML algorithm can evolve over time. Second,
the exact performance characteristics of an ML algorithm are often unknown
upfront and must be estimated, whereas many monitoring algorithms were
originally designed for settings where the pre-change data distribution is
known exactly (e.g. industrial manufacturing). Although there are procedures
that partially address these challenges (Dette and Gösmann, 2020; Zeileis and
Hornik, 2007; Gombay, 2017), we are not aware of a frequentist monitoring
procedure that adequately addresses both. To this end, we introduce a new
nonanticipating score-based CUSUM chart statistic and a computationally
efficient procedure for generating dynamic control limits. We note that
Bayesian approaches naturally handle nonstationarity and sources of
uncertainty, but posterior inference is computationally challenging over long
periods of time and for complex models of performance decay (Shiryaev, 1963;
West, 1986; Bhattacharya, 1994). As such, we leave the study of Bayesian
monitoring methods to future work.

Through simulation studies, we explore how various factors such as clinician
trust and model retraining impact our ability to detect performance decay. We
also apply score-based CUSUM monitoring to detect calibration decay of a ML-
based risk calculator for postoperative nausea and vomiting (PONV) on data
from the Multicenter Perioperative Outcomes Group (MPOG). For a locked PONV
risk calculator, the procedure fires an alarm during the COVID-19 pandemic. In
contrast, when we continually retrain the model, we find that the model can
steadily adapt to temporal shifts while remaining well-calibrated. Code is
available at https://github.com/jjfeng/monitoring_ML_CMI.

##  2 Two monitoring problems

Here we introduce the sequential monitoring problem in the “standard” setting
and then extend this to the setting with CMI. In the following section, we
discuss conditions under which a monitoring procedure designed for the former
setting provides valid inference for the latter. For convenience, Table 1
summarizes the mathematical notation used in this paper.

###  2.1 The standard monitoring problem

Following the framework set forth in Chu et al. (1996), suppose one observes
covariates $Z_{t}$ and predictors $Y_{t}$ for times $t=1,\cdots,m$, during
which the conditional distribution $Y_{t}|Z_{t}$ is constant. For some integer
$K>1$, we are interested in detecting structural change in this conditional
distribution for sequentially arriving observations $(Z_{t},Y_{t})$ for times
$t=m+1,\cdots,mK$. We assume there is at most one changepoint $\kappa$ during
this monitoring period, where $\kappa=\lfloor m\kappa^{\text{rel}}\rfloor$ for
some $\kappa^{\text{rel}}\in(1,K)$. If such a changepoint exists, the
conditional distribution $Y_{t}|Z_{t}$ is assumed to follow some model with
parameters $(\theta,\delta\mathbbm{1}\\{t\geq\kappa\\})$, where
$\theta\in\mathbb{R}^{p}$ describes the pre-change distribution and
$\delta\in\mathbb{R}^{d}$ describes the shift. The outcome $Y_{t}$ is assumed
to be conditionally independent of historical data given $Z_{t}$, so that the
distribution can be factorized into

$\displaystyle\begin{split}\Pr\left(Y_{1},\cdots,Y_{t},Z_{1},\cdots,Z_{t}\right)=&\prod_{i=1}^{t}\Pr\left(Y_{i}\mid
Z_{i};\theta,\delta\mathbbm{1}\\{i\geq\kappa\\}\right)\Pr\left(Z_{i}\mid
Z_{1},Y_{1},\cdots,Z_{i-1},Y_{i-1};\eta\right)\end{split}$ (1)

for $t=1,\cdots,mK$, where $\eta$ describes the conditional distribution of
$Z_{i}$ with respect to past monitoring data. The hypothesis test of interest
is

$\displaystyle\begin{split}H_{0}&:Y_{i}|Z_{i}\sim(\theta,0)\quad\forall
i=1,\cdots,mK\\\ H_{1}&:\exists\delta,\kappa\text{ s.t.
}Y_{i}|Z_{i}\sim(\theta,\delta\mathbbm{1}\\{i\geq\kappa\rfloor\\})\quad\forall
i=1,\cdots,mK,\end{split}$ (2)

where $\theta$ is a nuisance parameter.

We define a monitoring procedure as one with chart statistic $C_{m}(t)$ and
dynamic control limit (DCL) $h_{m}(t)$ at times $t=m+1,\cdots,mK$. The
procedure fires an alarm when the chart statistic first exceeds the control
limit, i.e. $\hat{T}_{m}=\inf\left\\{t:C_{m}(t)>h_{m}(t)\right\\}.$ Assuming
the operating characteristics for a monitoring procedure have been established
for the standard setting, our interest is to understand when these properties
transfer to the CMI setting.

###  2.2 The CMI monitoring problem

Data are generated in the CMI setting as follows. Upon arrival of a patient
with covariates $(X_{t},\tilde{X}_{t})$ at time $t$, the ML algorithm
$\hat{f}_{t}:\mathcal{X}\mapsto\mathcal{Q}$ outputs prediction
$\hat{f}_{t}(X_{t})$. The clinician then takes into account various
factors—such as the prediction, patient covariates, and prior experience with
the algorithm—to decide treatment $A_{t}$, where $A_{t}=0$ indicates SOC and
$A_{t}=1$ indicates additional intervention. Finally patient outcome $Y_{t}$
is observed. Following the potential outcomes framework, let $Y_{t}(a)$
indicate the patient outcome if treatment $a$ is administered. We assume the
realized outcome $Y_{t}$ is equal to the potential outcome $Y_{t}(A_{t})$
(consistency). Note that $t$ denotes patient ordering, not repeated
observations per patient.

This paper will focus on monitoring ML algorithms that predict the outcome
under SOC $Y_{t}(0)$, though our results can be extended to monitor algorithms
that predict treatment-specific outcomes. The central question is whether we
can simply monitor the subsequence of patients who were assigned SOC, as their
outcomes remain unaltered. To this end, let $\tau_{i}$ be the time of the
$i$th patient who received SOC, which can be viewed as a random stopping time.
We refer to the data observed at time points $\\{\tau_{i}:i=1,\cdots,mK\\}$ as
SOC-only data.

In simple cases, we aim to detect structural change in
$Y_{t}(0)\mid\hat{f}_{t}(X_{t})$. This corresponds to monitoring the
calibration curve if $\hat{f}_{t}$ is risk prediction algorithm and
positive/negative predictive values if $\hat{f}_{t}$ is a binary classifier.
For static algorithms, a structural change can only occur if there is a shift
in the joint distribution $(X_{t},Y_{t}(0))$. For evolving ML algorithms, a
structural change can occur either due to a shift in $(X_{t},Y_{t}(0))$ and/or
the updating procedure fails to maintain this conditional relationship. Both
failure modes are reasons for alarm.

When we directly apply a standard monitoring procedure to SOC-only data, this
corresponds to replacing $Z_{i}$ with $\hat{f}_{\tau_{i}}(X_{\tau_{i}})$ and
$Y_{i}$ with $Y_{\tau_{i}}$ to test the hypothesis

$\displaystyle\begin{split}H_{0}&:Y_{t}(0)|\hat{f}_{t}(X_{t}),t=\tau_{i}\sim(\theta,0)\quad\forall
i=1,\cdots,mK\\\ H_{1}&:\exists\delta,\kappa\text{ s.t.
}Y_{t}(0)|\hat{f}_{t}(X_{t}),t=\tau_{i}\sim(\theta,\delta\mathbbm{1}\\{i\geq\kappa\\})\quad\forall
i=1,\cdots,mK.\end{split}$ (3)

Note that this assumes the changepoint corresponds to the random stopping time
$\tau_{\lfloor m\kappa^{\text{rel}}\rfloor}$. In more complex settings where
clinical decision-making also depends on $\tilde{X}_{t}$, it is useful to
monitor for shifts in $Y_{t}(0)|\hat{f}_{t}(X_{t}),\tilde{X}_{t}$. This
corresponds to replacing $Z_{i}$ with
$\left(\hat{f}_{\tau_{i}}(X_{\tau_{i}}),\tilde{X}_{\tau_{i}}\right)$ and
$Y_{i}$ with $Y_{\tau_{i}}$ in the standard monitoring procedure to test

$\displaystyle\begin{split}H_{0}&:Y_{t}(0)|\hat{f}_{t}(X_{t}),\tilde{X}_{t},t=\tau_{i}\sim(\theta,0)\quad\forall
i=1,\cdots,mK\\\ H_{1}&:\exists\delta,\kappa\text{ s.t.
}Y_{t}(0)|\hat{f}_{t}(X_{t}),\tilde{X}_{t},t=\tau_{i}\sim(\theta,\delta\mathbbm{1}\\{i\geq\kappa\\})\quad\forall
i=1,\cdots,mK\end{split}.$ (4)

One can view (4) as monitoring stratified calibration curves—a stricter notion
of calibration in the hierarchy defined by Van Calster et al. (2016)—or
stratified predictive values.

In general, the operating characteristics of a monitoring procedure designed
for the standard setting may not transfer to the CMI setting because (i) the
distribution of the SOC-only data may not factorize per (1) and (ii) the
parameters may be biased. Nevertheless, the following section highlights
ignorability conditions under which valid statistical inference does transfer.
We note that no positivity assumptions are required to establish control of
the false alarm rate, because we only monitor for structural change with
respect to the population with non-zero probability of receiving SOC, even
though the null hypothesis applies to the broader population. Nevertheless,
statistical power and consistency of the monitoring procedure does depend on
positivity, as we show in our numerical experiments.

##  3 Ignorability assumptions

We now describe two conditions in which standard monitoring procedures applied
directly to SOC-only data, ignoring issues of CMI, provide valid inference.
The first ignorability condition is the well-known conditional exchangeability
assumption, also known as no unmeasured confounding. The second ignorability
condition allows treatment decisions to depend on unmeasured confounders but
requires selection bias to remain constant over time.

###  3.1 Conditional exchangeability

The simplest version of the conditional exchangeability assumption states that
treatment assignment $A_{t}$ is conditionally independent of potential outcome
$Y_{t}(0)$ given prediction $\hat{f}_{t}(X_{t})$, i.e.

$\displaystyle Y_{t}(0)$ $\displaystyle\perp
A_{t}\mid\hat{f}_{t}(X_{t})\quad\forall t=1,2,\cdots.$ (5)

In addition, let us assume the potential outcome $Y_{t}(0)$ is conditionally
independent of data at prior time points per

$\displaystyle Y_{t}(0)$
$\displaystyle\perp\left(A_{1},\hat{f}_{1}(X_{1}),Y_{1}(0),\cdots,A_{t-1},\hat{f}_{t-1}(X_{t-1}),Y_{t-1}(0)\right)\mid\hat{f_{t}}(X_{t}),A_{t}=0\quad\forall
t=1,2,\cdots.$ (6)

These conditions hold, for instance, if the ML algorithm is static (i.e.
$\hat{f}_{t}(X_{t})\equiv\hat{f}(X_{t})$ for all $t$) and the treatment
propensity at time $t$ solely depends on the ML algorithm’s predictions.
Figure 1 top shows a more complex example, in which both the treatment
propensities and ML algorithm evolve. To ensure that the conditional
independence assumption (6) holds, there cannot be an arrow from
$Y_{t-1}(a_{t-1})$ to $\hat{f}_{t}(X_{t})$, because this would introduce
collider bias. As such, we must keep the data for model updating and
monitoring separate, through mechanisms such as online sample-splitting.

Figure 1:  The top single world intervention graph (SWIG) (Richardson and
Robins, 2013) is an example satisfying conditional exchangeability (5) and
conditional independence (6). Black and gray variables are observed and
unobserved, respectively. The bottom SWIG satisfies time-constant selection
bias (9) and conditional independence (6) per the assumptions listed in
Example 1. The dotted ovals indicate determinative causes where $A_{t}=1$
whenever $A_{t}^{\prime}=1$ for all $t$.

Given (6), the distribution of SOC-only data can be factorized into

$\displaystyle\begin{split}&\Pr\left(Y_{\tau_{1}},\cdots,Y_{\tau_{t}},\hat{f}_{\tau_{1}}(X_{\tau_{1}}),\cdots,\hat{f}_{\tau_{t}}(X_{\tau_{t}})\right)\\\
=&\prod_{i=1}^{t}\Pr\left(Y_{\tau_{i}}|\hat{f}_{\tau_{i}}(X_{\tau_{i}})\right)\Pr\left(\hat{f}_{\tau_{i}}(X_{\tau_{i}})\mid\hat{f}_{\tau_{1}}(X_{\tau_{1}}),\cdots,\hat{f}_{\tau_{i-1}}(X_{\tau_{i-1}}),Y_{\tau_{1}},\cdots,Y_{\tau_{i-1}};\eta\right).\end{split}$
(7)

Then per the conditional exchangeability assumption (5), we have that the
conditional distribution in the SOC-only data coincides with that of the
general population, i.e.

$\displaystyle\Pr(Y_{\tau_{i}}|\hat{f}_{\tau_{i}}(X_{\tau_{i}}))=$
$\displaystyle\Pr(Y_{t}(0)|\hat{f}_{t}(X_{t}),t=\tau_{i}),$ (8)

which is parameterized by $(\theta,\delta\mathbbm{1}\\{i\geq\kappa\\})$. So
the distribution of SOC-only data factorizes into the same form as (1) with
monitoring target $Y_{t}(0)|\hat{f}_{t}(X_{t}),t=\tau_{i}$ and a standard
monitoring procedure that ignores CMI would provide valid inference for (3).

As is typically the case, the propensity to treat may also depend on other
predictors $\tilde{X}_{t}$. We can extend (5) and (6) by conditioning on
$\tilde{X}_{t}$ as well. Using similar arguments as above, we can show that
ignoring CMI still provides valid inference for testing (4).

###  3.2 Time-constant selection bias

When unmeasured confounders $U_{t}$ exist, conditional exchangeability no
longer holds and conditioning on $A_{t}=0$ results in selection bias (Figure 1
bottom). That is, the observed relationship between the outcome and the ML
predictions in the SOC-only data does not necessarily correspond to that of
the general population. Nevertheless, we can show that the treatment
propensities are still ignorable if selection bias remains constant over time.

The simplest version of the time-constant selection bias assumption states
that there is some function $h:\mathcal{Q}\mapsto\mathbb{R}$ such that

$\displaystyle
E\left[Y_{t}(0)\mid\hat{f}_{t}(X_{t})=q\right]-E\left[Y_{t}(0)\mid\hat{f}_{t}(X_{t})=q,A_{t}=0\right]=h\left(q\right)\qquad\forall
q\in\mathcal{Q},\forall t=1,2,\cdots.$ (9)

Nevertheless, we can still recover shifts in the conditional risk in the
general population from the SOC-only data because this bias cancels out, i.e.

$\displaystyle\begin{split}&E\left[Y_{t}(0)\mid\hat{f}_{t}(X_{t})=q,A_{t}=0\right]-E\left[Y_{1}(0)\mid\hat{f}_{1}(X_{1})=q,A_{1}=0\right]\\\
=&E\left[Y_{t}(0)\mid\hat{f}_{t}(X_{t})=q\right]-E\left[Y_{1}(0)\mid\hat{f}_{1}(X_{1})=q\right]\qquad\forall
q\in\mathcal{Q},\forall t=1,2,\cdots.\end{split}$ (10)

So combining time-constant selection bias with the conditional independence
assumption (6), we can plug into the factorization (7)

$\displaystyle\Pr(Y_{\tau_{i}}|\hat{f}_{\tau_{i}}(X_{\tau_{i}})=q)=\Pr(Y_{t}(0)|\hat{f}_{t}(X_{t})=q,A_{t}=0,t=\tau_{i};\theta^{\prime},\delta\mathbbm{1}\\{i\geq\kappa\\})$
(11)

where $\theta^{\prime}$ describes the conditional distribution
$Y_{t}(0)|\hat{f}_{t}(X_{t}),A=0$ prior to the changepoint. Importantly, the
value of $\delta$ remains unbiased, even though $\theta^{\prime}$ may not
equal $\theta$. Since $\theta^{\prime}$ is treated as a nuisance parameter, a
standard monitoring procedure that ignores CMI can still provide valid
inference for testing (3). A straightforward extension is to condition on
additional variables $\tilde{X}_{t}$ in assumptions (6) and (9), which implies
that a standard monitoring procedure that ignores CMI provides valid inference
for testing (4).

When does time-constant selection bias hold? We discuss one set of conditions
in Example 1 of the Appendix for the DAG (Figure 1 bottom). Briefly, suppose
the clinician makes an initial treatment decision $A_{t}^{\prime}\in\\{0,1\\}$
that she modifies based on some unmeasured confounder $U_{t}$ (e.g. a
biomarker) One of the conditions is that $A_{t}=1$ whenever
$A_{t}^{\prime}=1$, also known as determinative causation (Hernán et al.,
2004; VanderWeele and Robins, 2007, 2009).

##  4 The score-based CUSUM

In this section, we describe a new score-based CUSUM procedure with DCLs to
detect performance deterioration in the CMI setting. An example control chart
for this procedure is shown in Figure 2. The procedure is specially designed
to address two key challenges that arise when monitoring ML-based risk
prediction algorithms. First, the predictor sequence may be nonstationary due
to changes in the model and/or the clinician’s interactions with the model
over time. Second, the initial performance of the algorithm may be unknown and
must be estimated. DCLs have been previously used to address nonstationarity,
but only for likelihood-based test statistics and in settings where the pre-
change and shift parameters are known (Zhang et al., 2017), as the
distribution of such test statistics can be calculated under such assumptions.
Instead, we propose to use a nonanticipative score-based chart statistic, in
that the score for the $i$-th observation is calculated with respect to
parameters estimated using only historical data (Lorden and Pollak, 2005).
This preserves the martingale structure of the chart statistic even when the
nuisance parameter $\theta$ for the pre-change distribution is continually re-
estimated. We can then derive its asymptotic distribution and efficiently
construct DCLs using a simple parametric bootstrap procedure.

Figure 2:  Example control chart for the score-based CUSUM. The chart
statistic and control limits are shown in the blue and orange lines,
respectively. An alarm is fired when the chart statistic exceeds the control
limit.

For ease of notation, we will present the monitoring procedures in the context
of the standard setting. The procedure can be applied directly to the CMI
setting if the aforementioned ignorability conditions hold.

We quantify the operating characteristics of our procedure in terms of the
false alarm rate and its statistical power. That is, the probability of firing
an alarm before the changepoint, $\Pr\left(\hat{T}_{m}<\kappa\right)$, should
be controlled at level $\alpha>0$. In addition, the procedure should be
consistent, in that $\lim_{m\rightarrow\infty}\Pr\left(\hat{T}_{m}\leq
mK\right)=1$ under $H_{1}$. All proofs are in the Appendix.

For the theoretical analyses, we suppose the parametric model for
$Y_{t}|Z_{t}$ is correctly specified. For concreteness, consider the following
examples, which we will use in our empirical analyses. The first model
describes both the pre-change distribution and the structural change on the
log odds scale using logistic regression, i.e.

$\displaystyle\Pr\left(Y_{t}=1\mid
Z_{t};\theta,\delta\mathbbm{1}\\{t\geq\kappa\\}\right)$
$\displaystyle=\frac{1}{1+\exp\left(-(\theta+\delta\mathbbm{1}\\{t\geq\kappa\\})^{\top}Z_{t}\right)}.$
(12)

The second model describes the pre-change distribution on the log odds scale
but the structural change on the risk scale using

$\displaystyle\Pr\left(Y_{t}=1\mid
Z_{t};\theta,\delta\mathbbm{1}\\{t\geq\kappa\\}\right)$
$\displaystyle=\left[\frac{1}{1+\exp\left(-\theta^{\top}Z_{t}\right)}+(\delta\mathbbm{1}\\{t\geq\kappa\\})^{\top}Z_{t}\right]_{[0,1]},$
(13)

where $[x]_{[0,1]}=\min(1,\max(0,x))$. When conditional exchangeability holds,
we can monitor for structural change on any scale and use either (12) or (13).
If only time-constant selection bias holds, we are restricted to modeling
shifts on the risk scale and may only use (13).

###  4.1 Known pre-change parameter

We begin with the simplest setting in which the true value of the nuisance
parameter $\theta$, denoted $\theta_{0}$, is known. This holds, for instance,
when the model is known to be well-calibrated and one aims to test (3). Note
that the noncontamination dataset not needed in this setting.

For observation $(Z_{t},Y_{t})$, the score vector with respect to $\delta$
under the null hypothesis is $\left.\nabla_{\delta}\log p\left(Y_{t}\mid
Z_{t};\theta_{0},\delta\right)\right|_{\delta=0}$. Because the conditional
mean of the score is zero prior to the changepoint and nonzero after, we will
monitor for shifts in the average score using the cumulative sum (CUSUM)
(Page, 1954). In particular, for candidate changepoint $t^{\prime}$, define
the cumulative score up to time $t$ as

$\displaystyle\psi^{(\text{known})}_{m}(t^{\prime},t)=\sum_{i=t^{\prime}}^{t}\left.\nabla_{\delta}\log
p\left(Y_{i}\mid Z_{i};\theta_{0},\delta\right)\right|_{\delta=0}.$ (14)

Since the true changepoint time is unknown, the score-based CUSUM with respect
to norm $\|\cdot\|$ is defined as

$\displaystyle
C_{m}^{(\text{known})}(t)=\max_{t^{\prime}=m+1,\cdots,t}\left\|\psi^{(\text{known})}_{m}(t^{\prime},t)\right\|.$
(15)

In our empirical analyses, we use $\|\cdot\|_{1}$, though one can consider
other norms (for instance, using $\|\cdot\|_{2}$ would be more similar to
Rao’s score statistic).

We define DCLs $h_{m}(t)$ recursively using an alpha spending approach. Let
$\alpha^{\text{rel}}:[1,K]\mapsto[0,1]$ be the alpha-spending function, where
$\alpha^{\text{rel}}$ is assumed to be continuous and monotonically non-
decreasing. Then $h_{m}(t)$ is the minimal threshold at which the conditional
false alarm rate up to time $t$ matches the prespecified alpha-spending rate,
i.e.

$\displaystyle\Pr\left(\exists t^{\prime}\in\\{m+1,2,...,t\\}\text{ such that
}C_{m}^{(\text{known})}(t^{\prime})>h_{m}(t^{\prime})\mid\\{Z_{t^{\prime}}:t^{\prime}=1,...,t\\}\right)\leq\alpha^{\text{rel}}(t/m)$
(16)

under the null. To calculate the DCLs, we resample outcomes $Y^{*}_{t}$ given
$Z_{t}$ and the known pre-change parameters for all time points $t$,
constructing sequences $\\{(Z_{t},Y^{*(b)}_{t}):t=m+1,\cdots,mK\\}$ for
$b=1,\cdots,B$. By choosing a sufficiently large value for $B$, we then
estimate the distribution of the chart statistic and compute the DCLs.

We can show that the monitoring procedure is consistent if the average score
after the changepoint is bounded away from zero under the alternative, i.e.
there is some $c>0$ and $K^{\prime}\in(\kappa^{\text{rel}},K]$ such that

$\displaystyle\lim_{m\rightarrow\infty}\left\|\frac{1}{m}\sum_{t=\lfloor
m\kappa^{\text{rel}}\rfloor}^{\lfloor
mK^{\prime}\rfloor}\mathbb{E}\left[\left.\nabla_{\delta}\log
p\left(Y_{t}|Z_{t};\theta_{0},\delta\right)\right|_{\delta=0}\right]\right\|\geq
c,$ (17)

and the martingale

$\displaystyle\sum_{t=\lfloor m\kappa^{\text{rel}}\rfloor}^{\lfloor
mK^{\prime}\rfloor}\left.\nabla_{\delta}\log
p\left(Y_{t}|Z_{t};\theta_{0},\delta\right)\right|_{\delta=0}-\mathbb{E}\left[\left.\nabla_{\delta}\log
p\left(Y_{t}^{*}|Z_{t};\theta_{0},\delta\right)\right|_{\delta=0}\mid
Z_{t}\right]$ (18)

satisfies the martingale central limit theorem.

###  4.2 Unknown pre-change parameter

When the pre-change parameter $\theta$ is unknown, we need to estimate its
value and adjust the DCLs to reflect this additional source of uncertainty.
This situation is likely to occur when calibration data is limited (e.g.
deploying a model at a new site) or when monitoring the more complex
conditional distribution $Y_{t}(0)|\hat{f}_{t}(X_{t}),\tilde{X}_{t}$. To this
end, let the maximum likelihood estimator (MLE) for $\theta$ up to time $t$,
denoted $\hat{\theta}_{m,t}$, be the solution to the estimating equation
$\sum_{i=1}^{t}\nabla_{\theta}\log p\left(Y_{i}\mid Z_{i};\theta,0\right)=0$.
We define the score-based CUSUM chart statistic

$\displaystyle
C_{m}^{(\text{plugin})}(t)=\max_{t^{\prime}=m+1,\cdots,t}\left\|\psi^{(\text{plugin})}_{m}(t^{\prime},t)\right\|,$
(19)

where

$\displaystyle\psi^{(\text{plugin})}_{m}(t^{\prime},t)=\sum_{i=t^{\prime}}^{t}\left.\nabla_{\delta}\log
p\left(Y_{i}\mid Z_{i};\hat{\theta}_{m,i-1},\delta\right)\right|_{\delta=0}.$
(20)

The score for observation $(Z_{i},Y_{i})$ uses $\hat{\theta}_{m,i-1}$ rather
than $\hat{\theta}_{m,i}$, so (20) is nonanticipative.

To determine the operating characteristics of $C_{m}^{(\text{plugin})}$ under
the null hypothesis, we require the following assumptions. We use
$\delta_{0}=0$ to denote the value of $\delta$ under the null.

######  Assumption 4.1.

Under the null, there is a zero-mean $(p+d)$-dimensional non-degenerate
gaussian process $U$ such that

$\max_{m+1\leq i\leq
mK}\left\|\left[\frac{1}{\sqrt{m}}\sum_{j=1}^{i}\left(\begin{array}[]{c}\nabla_{\theta}\log
p\left(Y_{j}\mid Z_{j};\theta_{0},\delta_{0}\right)\\\ \nabla_{\delta}\log
p\left(Y_{j}\mid
Z_{j};\theta_{0},\delta_{0}\right)\end{array}\right)\right]-\left(\begin{array}[]{c}U_{\theta}(i/m)\\\
U_{\delta}(i/m)\end{array}\right)\right\|=o_{p}\left(1\right).$

where $U_{\theta}$ and $U_{\delta}$ are $p$\- and $d$-dimensional,
respectively.

######  Assumption 4.2.

Under the null, $\hat{\theta}_{m,i}$ is asymptotically linear with a remainder
term that converges uniformly to zero, i.e.

$\displaystyle\max_{m<i\leq
mK}\sqrt{m}\left\|\left(\hat{\theta}_{m,i}-\theta_{0}\right)-\mathbb{E}\left[-\sum_{j=1}^{i}\nabla_{\theta}^{2}\log
p\left(Y_{j}\mid
Z_{j};\theta_{0},\delta_{0}\right)\right]^{-1}\sum_{j=1}^{i}\nabla_{\theta}\log
p\left(Y_{j}\mid Z_{j};\theta_{0},\delta_{0}\right)\right\|$
$\displaystyle=o_{p}\left(1\right).$

######  Assumption 4.3.

Under the null, there exist functions
$\Lambda_{0}:[1,K]\mapsto\mathbb{R}^{p\times p}$ and
$\bar{V}_{0}:[1,K]\mapsto\mathbb{R}^{d\times p}$ such that

$\displaystyle\max_{m<i\leq
mK}\left\|\Lambda_{0}^{-1}\left(\frac{i}{m}\right)-m\mathbb{E}\left[-\sum_{j=1}^{i}\nabla_{\theta}^{2}\log
p\left(Y_{j}\mid
Z_{j};\theta_{0},\delta_{0}\right)\right]^{-1}\right\|=o_{p}(1)$
$\displaystyle\bar{V}_{0}(t)=\mathbb{E}\left[\nabla_{\theta}\nabla_{\delta}\log
p\left(Y_{\lfloor mt\rfloor}\mid Z_{\lfloor
mt\rfloor};\theta_{0},\delta_{0}\right)\right]\quad\forall t\in[1,K].$

These assumptions hold, for instance, under piecewise local stationarity (Wu
and Zhou, 2018; Horváth et al., 2021). We can then prove that
$\psi_{m}^{(\text{plugin})}$ is well-approximated by the process

$\displaystyle\phi_{m}(t_{1},t_{2})=\sum_{i=t_{1}}^{t_{2}}\nabla_{\delta}\log
p\left(Y_{i}\mid
Z_{i};\theta_{0},\delta_{0}\right)+\sum_{i=t_{1}}^{t_{2}}V_{0}\left(\frac{i}{m}\right)\Lambda_{0}^{-1}\left(\frac{i-1}{m}\right)\sum_{j=1}^{i-1}\nabla_{\theta}\log
p\left(Y_{j}\mid Z_{j};\theta_{0},\delta_{0}\right)$ (21)

under the null and that the latter converges weakly, as formalized below.

######  Theorem 4.1.

Suppose the null hypothesis is true and that Assumptions 4.1, 4.2, and 4.3
hold. In addition, suppose the second and third derivatives are bounded (see
Appendix for the formal assumption). Then

$\displaystyle\max_{m<t_{1},t_{2}\leq
mK}\frac{1}{\sqrt{m}}\left\|\psi_{m}^{(\text{plugin})}(t_{1},t_{2})-\phi_{m}(t_{1},t_{2})\right\|=o_{p}(1).$

and

$\displaystyle\left\\{(\nu_{1},\nu_{2})\mapsto\frac{1}{\sqrt{m}}\psi_{m}^{(\text{plugin})}(\lfloor
m\nu_{1}\rfloor,\lfloor
m\nu_{2}\rfloor)\right\\}_{(\nu_{1},\nu_{2})\in\Delta}\Rightarrow\left\\{(\nu_{1},\nu_{2})\mapsto
U_{\delta}(\nu_{2})-U_{\delta}(\nu_{1})+\int_{\nu_{1}}^{\nu_{2}}\bar{V}_{0}(v)\Lambda_{0}^{-1}(v)U_{\theta}(v)dv\right\\}_{(\nu_{1},\nu_{2})\in\Delta}$

where
$\Delta=\left\\{(\nu_{1},\nu_{2}):\nu_{1}<\nu_{2},\nu_{1}\in[1,K],\nu_{2}\in[1,K]\right\\}$.

In addition, Theorem B.1 in the Appendix proves that the procedure is
consistent if analogous assumptions hold under the alternative.

Given Theorem 4.1, we can determine the DCLs by analyzing the distribution of
$\phi_{m}$ among the resampled sequences
$\\{(Z_{t},Y^{*(b)}_{t}):t=m+1,\cdots,mK\\}$. One caveat is that this
technically requires sampling from the true $\theta$, but we only have an
estimate for its value. As such, we perform the parametric bootstrap and
resample the $t$-th outcome using the MLE estimated up to time $t-1$. The full
monitoring procedure is outlined in Algorithm 1. Additional implementation
details are provided in the Appendix.

Algorithm 1 Pseudocode for score-based CUSUM procedure with dynamic control
limits

Select time factor $K$, alpha spending function $\alpha^{\text{rel}}$, and
number of bootstrap sequences $B$.

Let $\mathcal{B}_{m}=\\{1,...,B\\}$ represent the bootstrap sequences that
have not been rejected at time $m$.

Observe non-contaminated data $\\{(Z_{t},Y_{t}):t=1,\cdots,m\\}$.

Calculate MLE $\hat{\theta}_{m,m}$

for $b=1,...,B$ do

Resample outcome $Y_{t}^{*(b)}$ given $Z_{t}$ for $t=1,\cdots,m$, with
$\theta=\hat{\theta}_{m,m}$ and $\delta=0$.

end for

for $t=m+1,...,mK$ do

Observe $(Z_{t},Y_{t})$.

Calculate chart statistic $C_{m}^{(plugin)}(t)$ and MLE $\hat{\theta}_{m,t}$.

for $b\in\mathcal{B}_{t-1}$ do

Resample outcome $Y_{t}^{*(b)}$ given $Z_{t}$ with
$\theta=\hat{\theta}_{m,t-1}$ and $\delta=0$.

Compute $\phi_{m}(t^{\prime},t)$ for $t^{\prime}=m+1,\cdots,t-1$ for the
$b$-th bootstrap sequence.

Calculate
$C^{*,(b)}_{m}(t)=\max_{t^{\prime}\in\\{m+1,\cdots,t\\}}\phi_{m}(t^{\prime},t)$.

end for

Set $h_{m}(t)$ such that the proportion of bootstrap chart statistics
exceeding the DCL is
$\left|\left\\{b:b\in\mathcal{B}_{t-1},C^{*,(b)}_{m}(t)>h_{m}(t)\right\\}\right|/B=\alpha^{\text{rel}}\left(t/m\right)-\alpha^{\text{rel}}\left((t-1)/m\right).$

Define $\mathcal{B}_{t}=\\{b:b\in\mathcal{B}_{t-1},C^{*,(b)}_{m}(t)\leq
h_{m}(t)\\}$.

if $C_{m}^{(\text{plugin})}(t)>h_{m}(t)$ then

Fire an alarm. Break.

end if

end for

##  5 Simulation studies

We now investigate how the score-based CUSUM is influenced by various factors,
including continual retraining of the ML algorithm, magnitude of the
structural change, and clinician trust. We vary each factor individually to
isolate its impact, though many factors can co-occur in real-world settings.

Details for all empirical experiments are in the Appendix. In all simulations,
the outcome $Y_{t}(0)$ given predictors $X_{t}\in\mathbb{R}^{p^{\prime}}$ and
$\tilde{X}_{t}\in\mathbb{R}$ is generated from a logistic regression model.
Unless specified otherwise, we set $p^{\prime}=8$ and train the risk
prediction algorithm using logistic regression. The nominal false alarm rate
for the score-based CUSUM is set to $\alpha=0.1$.

For comparison, we also implement Bayesian monitoring using Stan (Carpenter et
al., 2017), which performs Hamiltonian Monte Carlo (HMC). Additional
implementation details are in the Appendix. While our findings demonstrate
that Bayesian monitoring can yield similar results for appropriately chosen
priors, a number of major limitations must be addressed before it can be
recommended for practical use. First, there are currently no ready-to-use
software packages for long-term monitoring of the model (12). HMC is not
designed for sequential monitoring because it runs posterior inference from
scratch for every observation. In our experiments, HMC ran on the order of
hours whereas the score-based CUSUM ran on the order of seconds. Second,
posterior inference is generally difficult for models with structural
constraints such as (13), so we resorted to approximate inference when using
this model. Finally, Bayesian inference can be sensitive to misspecification
of the model and/or prior, so more concrete recommendations are needed on
their selection.

###  5.1 False alarm rate control

We begin with evaluating false alarm rate control of score-based CUSUM
monitoring in finite samples. The data are simulated under the null and
satisfy either the assumption of conditional exchangeability or time-constant
selection bias. A shift in the treatment propensities is introduced halfway
through the monitoring period.

Here we consider a locked ML algorithm. In the Conditional Exchangeability
simulation, the treatment propensities are generated according to a logistic
regression model with only $\hat{f}(X_{t})$ and $\tilde{X}_{t}$ as inputs; the
outcome is generated according to (12). In the Time-constant selection bias
simulation, the treatment propensities and outcome are generated per Example 1
of the Appendix. We consider two versions of both assumptions: one that only
conditions on $\hat{f}(X_{t})$ and another that conditions on
$(\hat{f}(X_{t}),\tilde{X}_{t})$. The former leads us to test (3) and the
latter leads us to test (4).

Conditional exchangeability

Time-constant selection bias

Figure 3:  Cumulative distribution of alarm times for score-based CUSUM
monitoring of a locked model $\hat{f}$ under the null. The assumption of
conditional exchangeability and time-constant selection bias are satisfied in
the top and bottom rows, respectively. In the left column, the assumptions
hold when conditioning on $\hat{f}(X_{t})$; in the right column, the
assumptions hold when conditioning on $(\hat{f}(X_{t}),\tilde{X}_{t})$. The
target false alarm rate is 0.1, which is achieved as the size of the non-
contaminated dataset $m$ increases.

We vary the size of the noncontamination dataset size $m$ and monitor for $mK$
time points, with $K$ set to $4$. As shown in Figure 3, Type I error is
inflated for small values of $m$, but converges to the nominal rate once $m$
is sufficiently large.

###  5.2 Monitoring a continually retrained model in stationary settings

An important safety requirement for continually retrained models is that they
remain well-calibrated over time. By wrapping such models within a monitoring
procedure, we can continually check that this requirement is met while
allowing model discrimination to change time. This approach is particularly
beneficial for online learning procedures that only provide weak (or no)
performance guarantees, such as black-box models. Here we investigate the
alarm rates for continually retrained ridge-penalized logistic regression
models and gradient boosted trees (GBT) in the stationary setting where data
is independently and identically distributed. This simulation also
simultaneously investigates the sensitivity of our monitoring procedure to
misspecification of the conditional distribution, since the monitoring model
is unlikely to hold uniformly over time.

We simulate a higher-dimensional setting with $p^{\prime}=50$ and a propensity
model that satisfies conditional exchageability. Both risk prediction
algorithms were continually retrained on all prior SOC data for a fixed set of
hyperparameters. Because GBT can be poorly calibrated, we recalibrated model
updates using Platt scaling. The monitoring procedures were implemented to
detect shifts on the logit and risk scales using (12) and (13), respectively.

Ridge-penalized logistic regression

Gradient boosted trees with Platt scaling

Figure 4:  Monitoring continually retrained ML risk prediction algorithms
using ridge-penalized logistic regression (top) and gradient boosted trees
with Platt scaling (bottom). Calibration curves of model updates from a single
replicate are plotted on the risk scale and logit scale in the first and
second columns, respectively. The average AUC of the model updates are shown
in the third column. The cumulative distribution of the alarm times are shown
for score-based CUSUM monitoring (fourth column) and Bayesian changepoint
monitoring (fifth column), which were implemented to detect for shifts on
either the risk scale or the logit scale.

For ridge-penalized logistic regression, the CUSUM alarm rate is close to the
nominal rate of $\alpha=0.1$ when monitoring on the risk scale and increases
to 0.16 when monitoring on the logit scale (Figure 4). Similar trends were
observed when using Bayesian monitoring. The difference in alarm rates is
explained by the calibration curves of the model updates: these curves
coincide nearly perfectly with the ideal diagonal line on the risk scale but
deviate slightly on the logit scale. As such, we recommend monitoring shifts
on the risk scale since (i) the risk scale is more meaningful in practice and
(ii) it is less sensitive to model misspecification. These results are highly
promising, since it shows how the false alarm rate for a continually updated
model can actually match that of a locked model and allow for steady
improvement in model discrimination over time.

For GBT, we found that Platt scaling improved model calibration but was unable
to maintain perfect calibration uniformly over time. Using the CUSUM, alarm
rates were 0.19 and 0.24 when monitoring on the risk and logit scales,
respectively; results from Bayesian inference were similar. This highlights
how one-time model recalibration is insufficient for maintaining model
reliability. Instead, one may need online calibration methods such as (Feng et
al., 2022a).

###  5.3 Big shifts. Little shifts. Locked models. Evolving models.

We now investigate alarm rates in the presence of a structural change. We vary
both the magnitude of the structural change and examine how model retraining
may affect alarm rates. In simulations large-shift and little-shift, we shrink
the coefficients of the outcome model at $t=50$ by 80% and 50%, respectively.
To retrain the ML model in nonstationary settings, we use an exponentially
weighted average forecaster (EWAF), a classic online learning algorithm that
can adapt to adversarial dataset shifts (Cesa-Bianchi and Lugosi, 2006; Feng,
2021).

Big Shift

Small Shift

Figure 5:  Monitoring locked versus continually retrained models in the
presence of big and small distribution shifts (top and bottom rows,
respectively). The calibration curves of the locked and continually retrained
models (first and second columns, respectively) are plotted over time,
including before and after the changepoint. The third and fourth columns show
the cumulative distribution of alarm times for score-based CUSUM monitoring
and Bayesian changepoint monitoring, respectively. The changepoint is
indicated in the right two columns by a dashed vertical line.

Across all scenarios, we find that the time to an alarm is longer if the
magnitude of the structural change is smaller (Figure 5). Perhaps more
interestingly, we consistently find that model retraining increases the time
to an alarm, depending on the quality of the model updates. In the large shift
simulation, the median time to alarm is delayed by 80 observations, because
the number of observations needed to detect a shift is small whereas that
needed to restore model calibration is much larger. In contrast, the alarm
rate drops to below 50% with model retraining in the small-shift simulation,
because detection requires many more observations and the EWAF was often able
to recalibrate the model during this time. Similar patterns can be seen when
using Bayesian monitoring. In summary, the alarm time can be viewed as the
result of a competition between model monitoring and updating. By designing a
sufficiently fast and adaptive model updating procedure, one can significantly
extend the total lifetime of a risk prediction algorithm.

###  5.4 Clinician trust

Intuitively, clinician trust can interfere with our ability to detect
performance decay. To investigate this, we simulate three levels of clinician
trust. In the No-trust simulation, the treatment propensity is uncorrelated
with the model’s prediction. In Calibrated-trust, the rate of treatment is the
same as the predicted risk. In Over-trust, the decision to treat essentially
thresholds on the model’s predictions, so that nearly everyone above 50%
predicted risk is treated and nearly everyone below 50% is untreated. We also
simulate two types of structural change: one in which shifts in the risk are
symmetric (symmetric-shift) and one in which risks shift the most among
patients with the highest initial risk (high-risk).

Shifts are largest in the high- and low-risk populations

Shifts are largest in only the high-risk population

Figure 6:  We simulate no trust, calibrated trust, and over trust in the ML
algorithm to understand its impact on detection delay. We simulate two types
of structural changes in the outcome model: shifts in the conditional risk
that is symmetric and largest among patients with the lowest and highest risks
(top row); and shifts that are largest only among patients with the highest
risks (bottom row). The left column shows how shifts in the conditional risk
vary with respect to the probability of a patient being assigned SOC, and thus
their probability of being sampled for monitoring. The middle and right
columns show the cumulative distribution of alarm times using score-based
CUSUM and Bayesian monitoring, respectively.

In the symmetric-shift simulation, increasing clinician trust had little
impact on alarm times, and may have even decreased detection delay (Figure 6
top). On the other hand, increasing clinician trust substantially increased
detection delay in the high-risk simulation (Figure 6 bottom). This difference
is explained by the distribution of risk shifts in the monitored patient
population. In the symmetric-shift simulation, increasing trust tended to
increase the representation of patients with lower predicted risks, which
corresponds to those experiencing bigger shifts in their risk. Whereas in the
high-risk simulation, increasing trust increased the representation of
subjects experiencing smaller shifts.

These results have an important practical implication: when designing
monitoring strategies, we should incorporate any prior knowledge regarding
which populations are likely to experience distribution shifts. If changes are
likely to be concentrated among subjects with high treatment propensities,
passive monitoring of SOC-only data may not sufficient. Instead, one may need
to actively augment the amount of information in the monitoring data. Prior
work has suggested educating healthcare providers on the appropriate use of ML
and warning against over-reliance (Finlayson et al., 2021) or randomizing
patients to receive no ML-based recommendations when it is deemed ethical
(Harris et al., 2022).

###  5.5 Sensitivity analysis of the time-constant selection bias assumption

Finally, we explore how violations of the time-constant selection bias
assumption can impact detection delay of a structural change. We first
simulate data that satisfies this time-constancy assumption based on Example
1. Then we introduce violations of this assumption by adding an edge from the
unmeasured confounder $U_{t}$ to the final treatment decision $A_{t}$ in the
DAG and setting a non-zero edge weight at times $t=100$ or $300$. Such shifts
in the propensity model could occur if, say, a clinician suspects performance
of the ML algorithm has decayed and decides to place more weight in the
unmeasured risk factor $U_{t}$.

Because the simulated violation is designed to dampen the shift observed in
the data, we find that detection delay increases as the violation occurs
earlier in time (Figure 7). In the worst case scenario, the structural change
and the shift in the treatment propensities occurs at the same time ($t=100$)
and power drops by 30%. Nevertheless, such a scenario is unlikely to happen in
practice since it assumes clinicians know exactly when performance decays.

In this simulation, we find that the power of Bayesian monitoring is much
lower than that for the score-based CUSUM. This is likely due to the
sensitivity of Bayesian inference to model misspecification: the monitoring
model assumes a single changepoint, whereas the observed data distribution
shifts at two time points. Consequently, its power drops by over 40%. In fact,
even without violations of the time-constancy assumption, the power of the
Bayesian procedure is much lower than that of the score-based CUSUM. This may
be due to difficulties in performing posterior inference for (13), which is
only partially differentiable.

Figure 7:  Monitoring a structural change at time $t=100$ (dashed vertical
line) when the assumption of time-constant selection bias is violated. We
introduce a shift in the treatment propensities that violates the time-
constant selection bias assumption at times $t=100$ and $t=300$. We also
simulate no violation of the time-constant selection bias by never introducing
this shift in the treatment propensities. For different risk prediction
values, we plot the conditional risks among the SOC-only population over time,
which are biased for those among the general population. We also plot the
oracle conditional risk for comparison. The middle and right columns show the
cumulative distribution of alarm times for score-based CUSUM monitoring and
Bayesian monitoring, respectively.

##  6 Monitoring a ML-based PONV risk calculator

Postoperative nausea and/or vomiting (PONV) is one of the most common side
effects of anesthesia. To reduce PONV rates, various risk prediction models
have been developed to guide the use of antiemetics. Here we simulate
monitoring a ML-based PONV risk calculator based on retrospective data from
the UCSF MPOG registry ($n=2434$). We define a patient as receiving additional
care if they received at least two antiemetics.

Using data from January 2018 to May 2019, we trained a random forest (RF) to
predict risk of PONV based on preoperative variables, including biological
sex, smoking status, age, ASA score, and blood test results. We then locked
the model, used the first 200 patients to initialize the monitoring
procedures, and started monitoring mid-December 2019. We suppose the
conditional exchangeability assumption holds with respect to
$\hat{f}_{t}(X_{t})$ and modeled the data using (12). Control limits were set
so that the false alarm rate/probability is 20%.

Monitoring a locked model

Monitoring a continually updated model

Figure 8:  Control charts for monitoring a ML-based risk prediction model for
Post-operative Nausea and Vomiting (PONV) using score-based monitoring (left)
and Bayesian inference (right). The top control charts monitor a locked model,
starting before the COVID pandemic. The control charts in the second row
monitor the continually retrained model, starting midway through the pandemic.
The ROC and calibration curves of the averaged models before and after January
2021 are shown in the bottom row.

The score-based and Bayesian monitoring procedures fired alarms in late 2020
(Figure 8 top). Both control charts suggest the shift started in May 2020, as
their chart statistics started to increase at this time. The detection of a
performance shift during the COVID-19 pandemic is not unexpected and may be
explained by many causes: (i) the anesthesia department had implemented
changes in antiemetic medication administration; (ii) there was a big shift in
the type of patients who received surgery during this time because of the
pandemic; and (iii) exposure to the SARS-Cov2 virus affected the overall
health of many patients. So if this ML algorithm was indeed deployed in the
hospital, it would be important to have a dedicated quality
assurange/improvement (QA/QI) team to investigate what the root causes are and
determine which corrective actions to take.

We then considered a scenario where the RF was retrained every 10
observations, where we split the data stream with 40% dedicated for model
retraining and 60% for model monitoring. We initialized the monitoring
procedure using observations from 200 patients starting July 2019. Due to data
splitting, initialization spanned a longer time period and monitoring began
October 2020. Although the AUC of the RF dropped to 0.59 in the beginning of
the pandemic due to shifts in the data, continual model retraining improved
the AUC to 0.63 (Figure 8 bottom). Moreover, no shifts in the calibration
curve were detected during this time.

##  7 Related work

The problem of CMI is closely related to the problem of verification bias
(Begg and Greenes, 1983; Zhou, 1994; Alonzo and Pepe, 2005), which arises when
disease status is only verified for a subset of patients. Verification bias is
typically presented using a missing data framework. However, given the
fundamental link missing data and causal inference, the assumptions and
techniques used to address CMI and verification bias share many similarities.
For instance, our assumption of conditional exchangeability can be viewed as
generalization of the missing at random assumption.

Given this connection with the missing data literature, our work is also
closely related to monitoring time-to-event data (Sego et al., 2009; Gandy et
al., 2010; Sun et al., 2014), many of which make similar conditional
independence assumptions. However, prior work assume that the pre-change
distribution is known, that the structural change has a very specific form
that is summarized by a single parameter (e.g. a proportional increase in the
hazard), and that the sequence of predictors are stationary. None of these
necessarily hold in our setting. Recently, Xue et al. (2021) proposed a batch-
sequential procedure for monitoring time-to-event data that relaxes many of
these assumptions; however, it relies on large batch sizes to justify
asymptotic approximations. In contrast, the score-based CUSUM procedure can be
used in a fully sequential manner.

Finally, we note that this work is part of the growing literature on
sequential testing of causal quantities using observational data (Li et al.,
2011; Cook et al., 2015; Waudby-Smith et al., 2021). Unlike these works, our
goal is test for structural change rather than the value of a time-constant
parameter.

##  8 Discussion

Although CMI can complicate performance monitoring of risk prediction models,
we have shown that CMI is ignorable when monitoring conditional performance
measures under either the assumption of conditional exchangeability or time-
constant selection bias. We introduce a new score-based CUSUM procedure with
DCLs that provides false alarm rate control, even if the ML algorithm and/or
the clinician’s interactions with the algorithm evolve over time.

Whereas prior work on performance monitoring has been restricted to locked
algorithms (Feng et al., 2022b), this work extends these ideas to monitor
continually updated algorithms. By wrapping online learning methods within a
monitoring framework, we greatly strengthen the performance guarantees of
these algorithms. In particular, many online learning methods only control
average performance over long time periods and do not protect against sudden
performance drops (Cesa-Bianchi and Lugosi, 2006; Feng, 2021). Using
monitoring procedures, we can address the latter.

Many issues warrant further investigation. Although we found the score-based
CUSUM to be somewhat robust to model misspecification, it will be important to
further improve robustness of the method. We have also assumed that the
treatment option is always clearly delineated and that patient outcomes are
observed immediately. This may not hold in certain settings, so future work
should consider more complex treatment variables and account for delays in
reporting. Finally, monitoring using SOC-only data is most powerful when
treatment decisions do not mask areas of major performance decay. Future work
should explore ways to address extreme treatment propensities, such as the
inclusion of other data sources and randomization.

## Acknowledgments

This work was greatly improved by helpful suggestions and feedback from Alan
Hubbard, Fan Xia, and Si Wen. We are grateful to Daniel Lazzareschi for
sharing the UCSF MPOG data.

## Funding

This work was supported by the Food and Drug Administration (FDA) of the U.S.
Department of Health and Human Services (HHS) as part of a financial
assistance award Center of Excellence in Regulatory Science and Innovation
grant to University of California, San Francisco (UCSF) and Stanford
University, (U01FD005978). The contents are those of the author(s) and do not
necessarily represent the official views of, nor an endorsement, by FDA/HHS,
or the U.S. Government.

## References

  * Adams et al. [2022] Roy Adams, Katharine E Henry, Anirudh Sridharan, and Hossein et al. Soleimani.  Prospective, multi-site study of patient outcomes after implementation of the TREWS machine learning-based early warning system for sepsis.  _Nat. Med._ , 28(7):1455–1460, July 2022.  URL http://dx.doi.org/10.1038/s41591-022-01894-0.
  * Alonzo and Pepe [2005] Todd A Alonzo and Margaret Sullivan Pepe.  Assessing accuracy of a continuous screening test in the presence of verification bias.  _J. R. Stat. Soc. Ser. C Appl. Stat._ , 54(1):173–190, January 2005.  URL https://onlinelibrary.wiley.com/doi/10.1111/j.1467-9876.2005.00477.x.
  * Begg and Greenes [1983] Colin B Begg and Robert A Greenes.  Assessment of diagnostic tests when disease verification is subject to selection bias.  _Biometrics_ , 39(1):207–215, 1983.  URL http://www.jstor.org/stable/2530820.
  * Bhattacharya [1994] P K Bhattacharya.  Some aspects of change-point analysis.  In _Institute of Mathematical Statistics Lecture Notes - Monograph Series_ , Lecture notes-monograph series, pages 28–56. Institute of Mathematical Statistics, Hayward, CA, 1994.  URL https://projecteuclid.org/ebook/download?urlId=10.1214/lnms/1215463112&isFullBook=false.
  * Carpenter et al. [2017] Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell.  Stan: A probabilistic programming language.  _J. Stat. Softw._ , 76(1), 2017.  URL https://www.osti.gov/biblio/1430202.
  * Cesa-Bianchi and Lugosi [2006] Nicolo Cesa-Bianchi and Gabor Lugosi.  _Prediction, Learning, and Games_.  Cambridge University Press, March 2006.  URL https://play.google.com/store/books/details?id=zDnRBlazhfYC.
  * Chu et al. [1996] Chia-Shang James Chu, Maxwell Stinchcombe, and Halbert White.  Monitoring structural change.  _Econometrica_ , 64(5):1045–1065, 1996.  URL http://www.jstor.org/stable/2171955.
  * Cook et al. [2015] Andrea J Cook, Robert D Wellman, Jennifer C Nelson, Lisa A Jackson, and Ram C Tiwari.  Group sequential method for observational data by using generalized estimating equations: application to vaccine safety datalink.  _J. R. Stat. Soc. Ser. C Appl. Stat._ , 64(2):319–338, 2015.  URL http://www.jstor.org/stable/24771896.
  * Davis et al. [2017] Sharon E Davis, Thomas A Lasko, Guanhua Chen, Edward D Siew, and Michael E Matheny.  Calibration drift in regression and machine learning models for acute kidney injury.  _J. Am. Med. Inform. Assoc._ , 24(6):1052–1061, November 2017.  URL http://dx.doi.org/10.1093/jamia/ocx030.
  * Dette and Gösmann [2020] Holger Dette and Josua Gösmann.  A likelihood ratio approach to sequential change point detection for a general class of parameters.  _J. Am. Stat. Assoc._ , 115(531):1361–1377, July 2020.  URL https://doi.org/10.1080/01621459.2019.1630562.
  * Feng [2021] Jean Feng.  Learning to safely approve updates to machine learning algorithms.  In _Proceedings of the Conference on Health, Inference, and Learning_ , CHIL ’21, pages 164–173, New York, NY, USA, April 2021. Association for Computing Machinery.  URL https://doi.org/10.1145/3450439.3451864.
  * Feng et al. [2022a] Jean Feng, Alexej Gossmann, Berkman Sahiner, and Romain Pirracchio.  Bayesian logistic regression for online recalibration and revision of risk prediction models with performance guarantees.  _Journal of the American Medical Informatics Association_ , 2022a.  URL http://dx.doi.org/10.1093/jamia/ocab280.
  * Feng et al. [2022b] Jean Feng, Rachael V Phillips, Ivana Malenica, Andrew Bishara, Alan E Hubbard, Leo A Celi, and Romain Pirracchio.  Clinical artificial intelligence quality improvement: towards continual monitoring and updating of AI algorithms in healthcare.  _npj Digital Medicine_ , 5(1):1–9, May 2022b.  URL https://www.nature.com/articles/s41746-022-00611-y.
  * Finlayson et al. [2021] Samuel G Finlayson, Adarsh Subbaswamy, Karandeep Singh, John Bowers, Annabel Kupke, Jonathan Zittrain, Isaac S Kohane, and Suchi Saria.  The clinician and dataset shift in artificial intelligence.  _N. Engl. J. Med._ , 385(3):283–286, July 2021\.  URL http://dx.doi.org/10.1056/NEJMc2104626.
  * Gandy et al. [2010] A Gandy, J T Kvaløy, A Bottle, and F Zhou.  Risk-adjusted monitoring of time to event.  _Biometrika_ , 97(2):375–388, June 2010.  URL https://academic-oup-com.ucsf.idm.oclc.org/biomet/article-pdf/97/2/375/583479/asq004.pdf.
  * Gombay [2017] Edit Gombay.  Editor’s special invited paper: On the efficient score vector in sequential monitoring.  _Sequential Analysis_ , 36(4):435–466, October 2017.  URL https://doi.org/10.1080/07474946.2017.1394728.
  * Harris et al. [2022] Steve Harris, Tim Bonnici, Thomas Keen, Watjana Lilaonitkul, Mark J White, and Nel Swanepoel.  Clinical deployment environments: Five pillars of translational machine learning for health.  _Frontiers in Digital Health_ , 4, 2022.  URL https://www.frontiersin.org/articles/10.3389/fdgth.2022.939292.
  * Henry et al. [2022] Katharine E Henry, Roy Adams, and Cassandra et al. Parent.  Factors driving provider adoption of the TREWS machine learning-based early warning system and its effects on sepsis treatment timing.  _Nat. Med._ , 28(7):1447–1454, July 2022.  URL http://dx.doi.org/10.1038/s41591-022-01895-z.
  * Hernán et al. [2004] Miguel A Hernán, Sonia Hernández-Díaz, and James M Robins.  A structural approach to selection bias.  _Epidemiology_ , 15(5):615–625, September 2004\.  URL http://dx.doi.org/10.1097/01.ede.0000135174.63482.43.
  * Hickey et al. [2013] Graeme L Hickey, Stuart W Grant, Gavin J Murphy, Moninder Bhabra, Domenico Pagano, Katherine McAllister, Iain Buchan, and Ben Bridgewater.  Dynamic trends in cardiac surgery: why the logistic EuroSCORE is no longer suitable for contemporary cardiac surgery and implications for future risk models.  _Eur. J. Cardiothorac. Surg._ , 43(6):1146–1152, June 2013.  URL http://dx.doi.org/10.1093/ejcts/ezs584.
  * Horváth et al. [2021] Lajos Horváth, Curtis Miller, and Gregory Rice.  Detecting early or late changes in linear models with heteroscedastic errors.  _Scand. Stat. Theory Appl._ , 48(2):577–609, June 2021.  URL https://onlinelibrary.wiley.com/doi/10.1111/sjos.12507.
  * Kahn et al. [1996] M G Kahn, T C Bailey, S A Steib, V J Fraser, and W C Dunagan.  Statistical process control methods for expert system performance monitoring.  _J. Am. Med. Inform. Assoc._ , 3(4):258–269, July 1996.  URL http://dx.doi.org/10.1136/jamia.1996.96413133.
  * Lenert et al. [2019] Matthew C Lenert, Michael E Matheny, and Colin G Walsh.  Prognostic models will be victims of their own success, unless…  _J. Am. Med. Inform. Assoc._ , 26(12):1645–1650, December 2019.  URL http://dx.doi.org/10.1093/jamia/ocz145.
  * Li et al. [2011] Lingling Li, Martin Kulldorff, Jennifer C Nelson, and Andrea J Cook.  A propensity Score-Enhanced sequential analytic method for comparative drug safety surveillance.  _Statistics in Biosciences_ , 2011.  URL https://link.springer.com/content/pdf/10.1007/s12561-011-9034-5.pdf.
  * Lorden and Pollak [2005] Gary Lorden and Moshe Pollak.  Nonanticipating estimation applied to sequential analysis and changepoint detection.  _Annals of Statistics_ , 33(3):1422–1454, June 2005.  URL https://projecteuclid.org/journals/annals-of-statistics/volume-33/issue-3/Nonanticipating-estimation-applied-to-sequential-analysis-and-changepoint-detection/10.1214/009053605000000183.full.
  * McCormick et al. [2012] Tyler H McCormick, Adrian E Raftery, David Madigan, and Randall S Burd.  Dynamic logistic regression and dynamic model averaging for binary classification.  _Biometrics_ , 68(1):23–30, March 2012.  URL http://dx.doi.org/10.1111/j.1541-0420.2011.01645.x.
  * Page [1954] E S Page.  Continuous inspection schemes.  _Biometrika_ , 41(1/2):100–115, 1954.  URL http://www.jstor.org/stable/2333009.
  * Paxton et al. [2013] Chris Paxton, Alexandru Niculescu-Mizil, and Suchi Saria.  Developing predictive models using electronic medical records: challenges and pitfalls.  _AMIA Annu. Symp. Proc._ , 2013:1109–1115, November 2013\.  URL https://www.ncbi.nlm.nih.gov/pubmed/24551396.
  * Perdomo et al. [2020] Juan Perdomo, Tijana Zrnic, Celestine Mendler-Dünner, and Moritz Hardt.  Performative prediction.  In Hal Daumé Iii and Aarti Singh, editors, _Proceedings of the 37th International Conference on Machine Learning_ , volume 119 of _Proceedings of Machine Learning Research_ , pages 7599–7609. PMLR, 2020\.  URL https://proceedings.mlr.press/v119/perdomo20a.html.
  * Richardson and Robins [2013] T S Richardson and J M Robins.  Single world intervention graphs (SWIGs): A unification of the counterfactual and graphical approaches to causality.  _Center for the Statistics and the Social Sciences_ , 2013.  URL https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.294.7647&rep=rep1&type=pdf.
  * Rubin [1976] Donald B Rubin.  Inference and missing data.  _Biometrika_ , 63(3):581–592, December 1976.  URL https://academic-oup-com.ucsf.idm.oclc.org/biomet/article-pdf/63/3/581/756166/63-3-581.pdf.
  * Sego et al. [2009] Landon H Sego, Marion R Reynolds, Jr, and William H Woodall.  Risk-adjusted monitoring of survival times.  _Stat. Med._ , 28(9):1386–1401, April 2009.  URL http://dx.doi.org/10.1002/sim.3546.
  * Shiryaev [1963] A N Shiryaev.  On optimum methods in quickest detection problems.  _Theory Probab. Appl._ , 8(1):22–46, January 1963\.  URL https://doi.org/10.1137/1108002.
  * Sun et al. [2014] Rena Jie Sun, John D Kalbfleisch, and Douglas E Schaubel.  A weighted cumulative sum (WCUSUM) to monitor medical outcomes with dependent censoring.  _Stat. Med._ , 33(18):3114–3129, August 2014\.  URL http://dx.doi.org/10.1002/sim.6139.
  * Tartakovsky et al. [2014] Alexander Tartakovsky, Igor Nikiforov, and Michele Basseville.  _Sequential Analysis: Hypothesis Testing and Changepoint Detection_.  CRC Press, August 2014.  URL https://play.google.com/store/books/details?id=zhsbBAAAQBAJ.
  * U.S. Food and Drug Administration and Health Canada [2021] U.S. Food and Drug Administration and Health Canada.  Good machine learning practice for medical device development, October 2021.  URL https://www.fda.gov/medical-devices/software-medical-device-samd/good-machine-learning-practice-medical-device-development-guiding-principles.
  * Van Calster et al. [2016] Ben Van Calster, Daan Nieboer, Yvonne Vergouwe, Bavo De Cock, Michael J Pencina, and Ewout W Steyerberg.  A calibration hierarchy for risk models was defined: from utopia to empirical data.  _J. Clin. Epidemiol._ , 74:167–176, June 2016.  URL http://dx.doi.org/10.1016/j.jclinepi.2015.12.005.
  * VanderWeele and Robins [2007] Tyler J VanderWeele and James M Robins.  Directed acyclic graphs, sufficient causes, and the properties of conditioning on a common effect.  _Am. J. Epidemiol._ , 166(9):1096–1104, November 2007.  URL http://dx.doi.org/10.1093/aje/kwm179.
  * VanderWeele and Robins [2009] Tyler J VanderWeele and James M Robins.  Minimal sufficient causation and directed acyclic graphs.  _Annals of Statistics_ , 37(3):1437–1465, June 2009.  URL https://projecteuclid.org/journals/annals-of-statistics/volume-37/issue-3/Minimal-sufficient-causation-and-directed-acyclic-graphs/10.1214/08-AOS613.full.
  * Waudby-Smith et al. [2021] Ian Waudby-Smith, David Arbour, Ritwik Sinha, Edward H Kennedy, and Aaditya Ramdas.  Doubly robust confidence sequences for sequential causal inference.  March 2021.  URL http://arxiv.org/abs/2103.06476.
  * West [1986] Mike West.  Bayesian model monitoring.  _J. R. Stat. Soc._ , 48(1):70–78, September 1986\.  URL https://onlinelibrary.wiley.com/doi/10.1111/j.2517-6161.1986.tb01391.x.
  * Wu and Zhou [2018] Weichi Wu and Zhou Zhou.  Gradient-based structural change detection for nonstationary time series M-estimation.  _Ann. Stat._ , 46(3):1197–1224, 2018.  URL https://www.jstor.org/stable/26542822.
  * Xue et al. [2021] Yishu Xue, Jun Yan, and Elizabeth D Schifano.  Simultaneous monitoring for regression coefficients and baseline hazard profile in cox modeling of time-to-event data.  _Biostatistics_ , 22(4):756–771, October 2021\.  URL http://dx.doi.org/10.1093/biostatistics/kxz064.
  * Zeileis and Hornik [2007] Achim Zeileis and Kurt Hornik.  Generalized m-fluctuation tests for parameter instability.  _Stat. Neerl._ , 61(4):488–508, November 2007\.  URL https://onlinelibrary.wiley.com/doi/10.1111/j.1467-9574.2007.00371.x.
  * Zhang et al. [2017] Xiang Zhang, Justin B Loda, and William H Woodall.  Dynamic probability control limits for risk-adjusted CUSUM charts based on multiresponses.  _Stat. Med._ , 36(16):2547–2558, July 2017.  URL http://dx.doi.org/10.1002/sim.7312.
  * Zhou [1994] X H Zhou.  Effect of verification bias on positive and negative predictive values.  _Stat. Med._ , 13(17):1737–1745, September 1994\.  URL http://dx.doi.org/10.1002/sim.4780131705.

Table 1: Mathematical symbols Symbol |  Meaning
---|---
$X_{t}$ |  Patient covariates that go into the ML algorithm
$\tilde{X}_{t}$ |  Additional patient covariates used to make treatment decisions
$Y_{t}$ |  Patient outcome
$A_{t}$ |  Treatment assignment
$\tau_{1},\tau_{2},\cdots$ |  Indices for the subsequence of times at which the patient was assigned standard-of-care (SOC)
$Z_{t}$ |  Predictors in the standard monitoring setting
$\hat{f}_{t}$ |  The ML algorithm at time $t$
$\kappa^{\text{rel}}\in(1,K)$ |  Position of changepoint in relative time
$\kappa$ |  Changepoint in absolute time, equal to $\lfloor m\kappa^{\text{rel}}\rfloor$
$\theta$ |  Parameter indexing the pre-change distribution
$\delta$ |  Parameter indexing the structural change
$m$ |  Size of dataset needed for initialization of monitoring procedures (also known as non-contaminated data)
$m+1,\cdots,mK$ |  The time period for monitoring structural change
$C_{m}(t)$ |  The chart statistic of a monitoring procedure at time $t$
$h_{m}(t)$ |  The control limit of a monitoring procedure at time $t$
$\hat{T}_{m}$ |  Alarm time of a monitoring procedure, i.e. when the chart statistic first exceeds the control limit

##  Appendix A Example satisfying the time-constant selection bias assumption

######  Example 1.

Consider the bottom single world intervention graph (SWIG) in Figure 1, where
$U_{t}$ is an unmeasured confounder. By the rules of D-separation, we have
that (6) holds. Suppose the distribution $(U_{t},A_{t}^{\prime})$ is constant
over time; as such, we will drop the time indices when denoting their marginal
and conditional distributions. For the conditional risk model, assume there is
no-additive interaction with respect to time and $U_{t}$, i.e.

$\operatorname{\mathbb{E}}\left[Y_{t}(0)\mid
X_{t}=x,U_{t}=u\right]=g_{0}(x,u;\theta_{0})+g_{1}(x;\delta_{1})\mathbbm{1}\\{t>\kappa\\}.$

Because $A_{t}=0$ implies that $A_{t}^{\prime}=0$, we have $U_{t}\perp
X_{t}\mid f_{t}(X_{t})=q,A_{t}=0$. Then for all times $t$ and
$q\in\mathcal{Q}$, we have that

$\displaystyle\operatorname{\mathbb{E}}\left[Y_{t}(0)\mid\hat{f}_{t}(X_{t})=q\right]-\operatorname{\mathbb{E}}\left[Y_{t}(0)\mid\hat{f}_{t}(X_{t})=q,A_{t}=0\right]$
(22) $\displaystyle=$
$\displaystyle\int\left(g_{0}(x_{t},u;\theta_{0})+g_{1}(x_{t};\delta_{1})\mathbbm{1}\\{t>\kappa\\}\right)p(x_{t}\mid\hat{f}_{t}(x_{t})=q)\left[p(u)-p(u\mid
a^{\prime}=0)\right]dx_{t}du$ (23) $\displaystyle=$ $\displaystyle\int
g_{0}(x_{t},u;\theta_{0})p(x_{t}\mid\hat{f}_{t}(x_{t})=q)\left[p(u)-p(u\mid
a^{\prime}=0)\right]dx_{t}du.$ (24)

There are various conditions under which (24) is time-constant. One
requirement is that $g_{0}(x,u;\theta_{0})$ is additive, in that
$g_{0}(x,u;\theta_{0})=g_{0,0}(x;\theta_{0})+g_{0,1}(u;\theta_{0})$.
Alternatively, we require the ML algorithm to be locked ($\hat{f}_{t}=\hat{f}$
for all $t$) and the distribution of $X_{t}$ to not vary over time.

##  Appendix B Proofs for the score-based CUSUM

Let $\mathcal{Z}$ and $\mathcal{Y}$ be the domains for the predictors and
outcomes. Let $(\theta_{0},\delta_{0})$ and $(\theta_{0},\delta_{1})$
parameterize the pre-change and post-change distribution, where $\delta_{0}=0$
and $\delta_{1}\neq 0$. We assume that $p(y|z;\theta,\delta)$ is 3-times
continuously differentiable with respect to $(\theta,\delta)$. For
convenience, denote

$\displaystyle\Lambda_{m}(i)$
$\displaystyle=\mathbb{E}\left[-\sum_{j=1}^{i}\nabla_{\theta}^{2}\log
p\left(Y_{j}\mid Z_{j};\theta_{0},\delta_{0}\right)\right]$ $\displaystyle
V_{m}(i)$ $\displaystyle=\mathbb{E}\left[\nabla_{\theta}\nabla_{\delta}\log
p\left(Y_{i}\mid Z_{i};\theta_{0},\delta_{0}\right)\middle|Z_{i}\right].$

For $v\in[1,K]$, define the limit of the MLEs as
$\bar{\hat{\theta}}(v)=\lim_{m\rightarrow\infty}\hat{\theta}_{m,\lfloor
mv\rfloor}(v)$.

We use the symbol $\Rightarrow$ to mean weak convergence in the space under
consideration. Throughout, we will use $c$ (sometimes with subscripts) to
denote constants, which may vary across contexts. When we write
$Z_{m}\leq_{p}c$, this means that asymptotically as $m\rightarrow\infty$, the
random variable $Z_{m}$ is bounded by some constant $c$ with probability 1.

###  B.1 Asymptotics under the null

Here we prove asymptotic convergence of the chart statistic under the null. In
addition to the assumptions listed in the main manuscript, we will require the
second and third derivatives of the likelihood to be bounded as follows.

######  Assumption B.1.

There exist constants $c_{1},c_{2}>0$ and some neighborhood $B(\theta_{0},r)$
centered at $\theta_{0}$ with radius $r>0$ such that

$\displaystyle\sup_{\tilde{\theta}\in
B(\theta_{0},r)}\sup_{y\in\mathcal{Y},z\in\mathcal{Z}}\left\|\nabla_{\theta}^{2}\nabla_{\delta}\log
p\left(y\mid z;\tilde{\theta},\delta_{0}\right)\right\|_{\infty}$
$\displaystyle\leq c_{1}$ (25)
$\displaystyle\sup_{y\in\mathcal{Y},z\in\mathcal{Z}}\left\|\nabla_{\theta}\nabla_{\delta}\log
p\left(y\mid z;\theta_{0},\delta_{0}\right)\right\|_{\infty}$
$\displaystyle\leq_{p}c_{2}.$ (26)

######  Lemma B.1.

Suppose Assumptions 4.1, 4.2, 4.3, and B.1 hold. Define

$\tilde{\phi}_{m}(t_{1},t_{2})=\sum_{i=t_{1}}^{t_{2}}\nabla_{\delta}\log
p\left(Y_{i}\mid
Z_{i};\theta_{0},\delta_{0}\right)+\sum_{i=t_{1}}^{t_{2}}V_{m}(i)\Lambda_{m}^{-1}(i-1)\sum_{j=1}^{i-1}\nabla_{\theta}\log
p\left(Y_{j}\mid Z_{j};\theta_{0},\delta_{0}\right).$

Under the null, we have

$\max_{m<t_{1},t_{2}\leq
mK}\frac{1}{\sqrt{m}}\left\|\psi_{m}^{(\text{plugin})}(t_{1},t_{2})-\tilde{\phi}_{m}(t_{1},t_{2})\right\|_{2}=o_{p}(1).$

###### Proof.

Consider the decomposition

$\displaystyle\frac{1}{\sqrt{m}}\left(\psi_{m}^{(\text{plugin})}(t_{1},t_{2})-\tilde{\phi}_{m}(t_{1},t_{2})\right)=$
$\displaystyle
R_{m}^{(1)}(t_{1},t_{2})+R_{m}^{(2)}(t_{1},t_{2})+R_{m}^{(3)}(t_{1},t_{2})$

where

$\displaystyle R_{m}^{(1)}(t_{1},t_{2})$
$\displaystyle=\frac{1}{\sqrt{m}}\sum_{i=t_{1}}^{t_{2}}\nabla_{\delta}\log
p\left(Y_{i}\mid
Z_{i};\hat{\theta}_{m,i-1},\delta_{0}\right)-\nabla_{\delta}\log
p\left(Y_{i}\mid
Z_{i};\theta_{0},\delta_{0}\right)-\nabla_{\theta}\nabla_{\delta}\log
p\left(Y_{i}\mid
Z_{i};\theta_{0},\delta_{0}\right)\left(\hat{\theta}_{m,i-1}-\theta_{0}\right)$
$\displaystyle R_{m}^{(2)}(t_{1},t_{2})$
$\displaystyle=\frac{1}{\sqrt{m}}\sum_{i=t_{1}}^{t_{2}}V_{m}(i)\left[\hat{\theta}_{m,i-1}-\theta_{0}-\Lambda_{m}^{-1}(i-1)\sum_{j=1}^{i-1}\nabla_{\theta}\log
p\left(Y_{j}\mid Z_{j};\theta_{0},\delta_{0}\right)\right]$ $\displaystyle
R_{m}^{(3)}(t_{1},t_{2})$
$\displaystyle=\frac{1}{\sqrt{m}}\sum_{i=t_{1}}^{t_{2}}\left[\nabla_{\theta}\nabla_{\delta}\log
p\left(Y_{i}\mid
Z_{i};\theta_{0},\delta_{0}\right)-V_{m}(i)\right]\left(\hat{\theta}_{m,i-1}-\theta_{0}\right).$

We will prove that each term in this decomposition is negligible, i.e.

$\max_{m<t_{1}<t_{2}\leq
mK}\left\|R_{m}^{(j)}(t_{1},t_{2})\right\|_{2}=o_{p}(1)\quad\forall j=1,2,3.$
(27)

#### First remainder term.

For any $\epsilon>0$, we have that

$\displaystyle\begin{split}\Pr\left(\max_{m<t_{1}<t_{2}\leq
mK}\left\|R_{m}^{(1)}(t_{1},t_{2})\right\|_{2}>\epsilon\right)\leq&\Pr\left(\max_{m<i\leq
mK}\left\|\hat{\theta}_{m,i}-\theta_{0}\right\|_{2}>c\frac{\log
m}{\sqrt{m}}\right)\\\ &+\Pr\left(\max_{m<t_{1}<t_{2}\leq
mK}\left\|R_{m}^{(1)}(t_{1},t_{2})\right\|_{2}>\epsilon,\max_{m<i\leq
mK}\left\|\hat{\theta}_{m,i}-\theta_{0}\right\|_{2}\leq c\frac{\log
m}{\sqrt{m}}\right).\end{split}$ (28)

The first summand on the RHS of (28) goes to zero because Assumptions 4.1,
4.2, and 4.3 imply

$\max_{m<i\leq
mK}\left\|\hat{\theta}_{m,i}-\theta_{0}\right\|_{2}=o_{p}\left(\frac{\log
m}{\sqrt{m}}\right).$ (29)

To bound the second summand, we have by Taylor’s theorem and Assumption B.1
that for sufficiently large $m$

$\displaystyle\left\|R_{m}^{(1)}(t_{1},t_{2})\right\|_{2}$
$\displaystyle\leq\frac{1}{2\sqrt{m}}\sum_{i=t_{1}}^{t_{2}}\left(\max_{\tilde{\theta}\in
B(\theta_{0},r)}\max_{z\in\mathcal{Z},y\in\mathcal{Y}}\left\|\nabla_{\theta}^{2}\nabla_{\delta}\log
p\left(y\mid
z;\tilde{\theta},\delta_{0}\right)\right\|_{\infty}\right)\left\|\hat{\theta}_{m,i-1}-\theta_{0}\right\|_{2}^{2}$
$\displaystyle\leq\frac{c_{1}}{2\sqrt{m}}\sum_{i=t_{1}}^{t_{2}}\left\|\hat{\theta}_{m,i-1}-\theta_{0}\right\|_{2}^{2}$
$\displaystyle\leq\frac{c_{1}}{2}\sqrt{m}\left(K-1\right)\max_{m<i\leq
mK}\left\|\hat{\theta}_{m,i}-\theta_{0}\right\|_{2}^{2}$
$\displaystyle=o_{p}\left(\left(K-1\right)\left(\log
m\right)^{2}/\sqrt{m}\right)$

for all $(t_{1},t_{2})$ where $m<t_{1}\leq t_{2}\leq mK$. So (27) holds for
$j=1$.

#### Second remainder term.

By Assumption B.1 and the Cauchy-Schwarz inequality, we have that

$\displaystyle\left\|R_{m}^{(2)}(t_{1},t_{2})\right\|_{2}$
$\displaystyle\leq\frac{c_{2}}{\sqrt{m}}\sum_{i=t_{1}}^{t_{2}}\max_{m<i\leq
mK}\left\|\hat{\theta}_{m,i-1}-\theta_{0}-\Lambda_{m}^{-1}(i-1)\sum_{j=1}^{i-1}\nabla_{\theta}\log
p\left(Y_{j}\mid Z_{j};\theta_{0},\delta_{0}\right)\right\|_{2}$
$\displaystyle\leq c\left(K-1\right)\sqrt{m}\max_{m<i\leq
mK}\left\|\hat{\theta}_{m,i}-\theta_{0}-\Lambda_{m}^{-1}(i)\sum_{j=1}^{i}\nabla_{\theta}\log
p\left(Y_{j}\mid Z_{j};\theta_{0},\delta_{0}\right)\right\|_{2}.$

Then by Assumption 4.2, this term is $o_{p}\left(1\right)$. So (27) holds for
$j=2$.

#### Third remainder term.

For any $\epsilon>0$, we have that

$\displaystyle\begin{split}&\Pr\left(\max_{t_{1},t_{2}}\left\|\sum_{i=t_{1}}^{t_{2}}\left(\nabla_{\theta}\nabla_{\delta}\log
p\left(Y_{i}\mid
Z_{i};\theta_{0},\delta_{0}\right)-V_{m}(i)\right)\left(\hat{\theta}_{m,i-1}-\theta_{0}\right)\right\|_{2}\geq\epsilon\sqrt{m}\right)\\\
\leq&\Pr\left(\max_{m<i\leq
mK}\left\|\hat{\theta}_{m,i}-\theta_{0}\right\|_{2}>c\frac{\log
m}{\sqrt{m}}\right)\\\
&+\Pr\left(\max_{t_{1},t_{2}}\left\|\sum_{i=t_{1}}^{t_{2}}\left(\nabla_{\theta}\nabla_{\delta}\log
p\left(Y_{i}\mid
Z_{i};\theta_{0},\delta_{0}\right)-V_{m}(i)\right)\left(\hat{\theta}_{m,i-1}-\theta_{0}\right)\right\|_{2}\geq\epsilon\sqrt{m},\max_{m<i\leq
mK}\left\|\hat{\theta}_{m,i}-\theta_{0}\right\|_{2}\leq c\frac{\log
m}{\sqrt{m}}\right).\end{split}$ (30)

Per (29), the first summand on the RHS of (30) goes to zero. The second
summand is bounded by

$\Pr\left(\max_{t_{1},t_{2}}\left\|\sum_{i=t_{1}}^{t_{2}}\left(\nabla_{\theta}\nabla_{\delta}\log
p\left(Y_{i}\mid
Z_{i};\theta_{0},\delta_{0}\right)-V_{m}(i)\right)\left(\hat{\theta}_{m,i-1}-\theta_{0}\right)\mathbbm{1}\left\\{\|\hat{\theta}_{m,i-1}-\theta_{0}\|_{2}\leq
c\frac{\log m}{\sqrt{m}}\right\\}\right\|_{2}\geq\epsilon\sqrt{m}\right).$
(31)

Because the outcome $Y_{i}$ is conditionally independent of past data given
$Z_{i}$, the elements in this summation form a martingale difference sequence,
i.e.

$\operatorname{\mathbb{E}}\left[G_{m}(i)\left(\hat{\theta}_{m,i-1}-\theta_{0}\right)\mathbbm{1}\left\\{\left\|\hat{\theta}_{m,i-1}-\theta_{0}\right\|_{2}\leq
c\frac{\log m}{\sqrt{m}}\right\\}\middle|\mathcal{F}_{i}\right]=0$

where we use the notational shorthand

$G_{m}(i)=\nabla_{\theta}\nabla_{\delta}\log p\left(Y_{i}\mid
Z_{i};\theta_{0},\delta_{0}\right)-V_{m}(i)$

and

$\mathcal{F}_{i}=\left(Z_{1},Y_{1},\cdots,Z_{i-1},Y_{i-1},Z_{i}\right).$

Moreover, by Assumption B.1, $G_{m}(i)$ is sub-Gaussian. That is, there is
some $\sigma^{2}>0$ such that for all $\lambda>0$, we have for all unit
vectors $u,v$ that

$\displaystyle\operatorname{\mathbb{E}}\left[\exp\left(\lambda
v^{\top}G_{m}(i)u\right)\middle|\mathcal{F}_{i-1}\right]$
$\displaystyle\leq\operatorname{\mathbb{E}}\left[\exp\left(\lambda^{2}\sigma^{2}\lambda_{\max}\left(G_{m}(i)\right)\right)\middle|\mathcal{F}_{i-1}\right].$

By the law of total expectations, we then have for any unit vector $v$ that

$\displaystyle\operatorname{\mathbb{E}}\left[\exp\left(\lambda
v^{\top}\sum_{i=t_{1}}^{t_{2}}G_{m}(i)\left(\hat{\theta}_{m,i-1}-\theta_{0}\right)\mathbbm{1}\left\\{\left\|\hat{\theta}_{m,i-1}-\theta_{0}\right\|_{2}\leq
c\frac{\log m}{\sqrt{m}}\right\\}\right)\right]$ $\displaystyle=$
$\displaystyle\operatorname{\mathbb{E}}\Biggr{[}\operatorname{\mathbb{E}}\left[\exp\left(\lambda
v^{\top}G_{m}(t_{2})\left(\hat{\theta}_{m,t_{2}-1}-\theta_{0}\right)\mathbbm{1}\left\\{\left\|\hat{\theta}_{m,t_{2}-1}-\theta_{0}\right\|_{2}\leq
c\frac{\log m}{\sqrt{m}}\right\\}\right)\middle|\mathcal{F}_{t_{2}-1}\right]$
$\displaystyle\qquad\times\exp\left(\lambda
v^{\top}\sum_{i=t_{1}}^{t_{2}-1}G_{m}(i)\left(\hat{\theta}_{m,i-1}-\theta_{0}\right)\mathbbm{1}\left\\{\left\|\hat{\theta}_{m,i-1}-\theta_{0}\right\|_{2}\leq
c\frac{\log m}{\sqrt{m}}\right\\}\right)\Biggr{]}$ $\displaystyle\leq$
$\displaystyle\exp\left(\lambda^{2}c^{2}\sigma^{2}\frac{\left(\log
m\right)^{2}}{m}\right)\operatorname{\mathbb{E}}\left[\exp\left(\lambda
v^{\top}\sum_{i=t_{1}}^{t_{2}-1}G_{m}(i)\left(\hat{\theta}_{m,i-1}-\theta_{0}\right)\mathbbm{1}\left\\{\left\|\hat{\theta}_{m,i-1}-\theta_{0}\right\|_{2}\leq
c\frac{\log m}{\sqrt{m}}\right\\}\right)\right]$ $\displaystyle\leq$
$\displaystyle\exp\left(\lambda^{2}c^{2}\sigma^{2}\left(K-1\right)\left(\log
m\right)^{2}\right).$

Using the Chernoff bound, we have that (31) is bounded by

$\displaystyle\sum_{m<t_{1}\leq t_{2}\leq
mK}\Pr\left(\left\|\sum_{i=t_{1}}^{t_{2}}\left(\nabla_{\theta}\nabla_{\delta}\log
p\left(Y_{i}\mid
Z_{i};\theta_{0},\delta_{0}\right)-V_{m}(i)\right)\left(\hat{\theta}_{m,i-1}-\theta_{0}\right)\mathbbm{1}\left\\{\left\|\hat{\theta}_{m,i-1}-\theta_{0}\right\|_{2}\leq
c\frac{\log m}{\sqrt{m}}\right\\}\right\|_{2}\geq\epsilon\sqrt{m}\right)$
$\displaystyle\leq$ $\displaystyle
m^{2}\left(K-1\right)^{2}\exp\left(\lambda^{2}c^{2}\sigma^{2}\left(K-1\right)\left(\log
m\right)^{2}-\epsilon^{2}m\right).$

This converges to zero as $m\rightarrow\infty$, so (27) holds for $j=3$. ∎

Using Lemma B.1 above, we are now ready to prove Theorem 4.1.

###### Proof of Theorem 4.1.

Consider the decomposition

$\displaystyle\tilde{\phi}_{m}(t_{1},t_{2})$
$\displaystyle=\phi_{m}(t_{1},t_{2})+R_{m}^{(1)}(t_{1},t_{2})+R_{m}^{(2)}(t_{1},t_{2})$
(32)

with remainder terms defined as

$\displaystyle R_{m}^{(1)}(t_{1},t_{2})=$
$\displaystyle\frac{1}{\sqrt{m}}\sum_{i=t_{1}}^{t_{2}}\left(V_{m}(i)-\bar{V}_{0}\left(i/m\right)\right)\Lambda_{m}^{-1}(i-1)\sum_{j=1}^{i-1}\nabla_{\theta}\log
p\left(Y_{j}\mid Z_{j};\theta_{0},\delta_{0}\right)$ $\displaystyle
R_{m}^{(2)}(t_{1},t_{2})=$
$\displaystyle\frac{1}{\sqrt{m}}\sum_{i=t_{1}}^{t_{2}}\bar{V}_{0}\left(i/m\right)\left(\Lambda_{m}^{-1}(i-1)-\frac{1}{m}\Lambda_{0}^{-1}\left(\frac{i-1}{m}\right)\right)\sum_{j=1}^{i-1}\nabla_{\theta}\log
p\left(Y_{j}\mid Z_{j};\theta_{0},\delta_{0}\right).$

We will first show the remainder terms are negligible, i.e.

$\max_{m<t_{1}<t_{2}\leq mK}\|R_{m}^{(j)}(t_{1},t_{2})\|=o_{p}(1)\quad\forall
j=1,2.$ (33)

#### First remainder term.

To bound the first remainder, note that the partial sums form a martingale due
to Assumption 4.3, i.e.

$\operatorname{\mathbb{E}}\left[\left(V_{m}(i)-\bar{V}_{0}\left(i/m\right)\right)\Lambda_{m}^{-1}(i-1)\sum_{j=1}^{i-1}\nabla_{\theta}\log
p\left(Y_{j}\mid
Z_{j};\theta_{0},\delta_{0}\right)\middle|\mathcal{F}_{i}\right]=0.$

Moreover, $\left(V_{m}(i)-\bar{V}_{0}(i/m)\right)\Lambda_{m}^{-1}(i)$ is sub-
Gaussian and

$\displaystyle\max_{m<i\leq mK}\left\|\sum_{j=1}^{i-1}\nabla_{\theta}\log
p\left(Y_{j}\mid Z_{j};\theta_{0},\delta_{0}\right)\right\|_{2}$
$\displaystyle=o_{p}\left(\frac{\log m}{\sqrt{m}}\right)$

per Assumptions 4.1 and 4.3. As such, we can use a similar martingale argument
as the previous lemma to prove that (33) is satisfied for $j=1$.

#### Second remainder term.

By the Cauchy-Schwarz inequality, we have that

$\displaystyle\frac{1}{\sqrt{m}}\sum_{i=t_{1}}^{t_{2}}\bar{V}_{0}\left(i/m\right)\left(\Lambda_{m}^{-1}(i-1)-\frac{1}{m}\Lambda_{0}^{-1}\left(\frac{i-1}{m}\right)\right)\sum_{j=1}^{i-1}\nabla_{\theta}\log
p\left(Y_{j}\mid Z_{j};\theta_{0},\delta_{0}\right)$ $\displaystyle\leq
c\sum_{i=t_{1}}^{t_{2}}\left\|\Lambda_{m}^{-1}(i)-\frac{1}{m}\Lambda_{0}^{-1}(i/m)\right\|_{2}\left\|\frac{1}{\sqrt{m}}\sum_{j=1}^{i-1}\nabla_{\theta}\log
p\left(Y_{j}\mid Z_{j};\theta_{0},\delta_{0}\right)\right\|_{2}$
$\displaystyle\leq
c\left(K-1\right)\left(\max_{i=m+1,\cdots,mK}\left\|m\Lambda_{m}^{-1}(i)-\Lambda_{0}^{-1}(i/m)\right\|_{2}\right)\left(\max_{m<i\leq
mK}\left\|\frac{1}{\sqrt{m}}\sum_{j=1}^{i-1}\nabla_{\theta}\log
p\left(Y_{j}\mid Z_{j};\theta_{0},\delta_{0}\right)\right\|_{2}\right)$

By Assumptions 4.1 and 4.3, it follows that (33) for $j=2$.

In addition, by Assumption 4.1, we have that

$\left\\{\nu\mapsto\frac{1}{\sqrt{m}}\sum_{j=1}^{\lfloor
m\nu\rfloor}\left(\begin{array}[]{c}\nabla_{\theta}\log p\left(Y_{j}\mid
Z_{j};\theta_{0},\delta_{0}\right)\\\
V_{0}(\nu)\Lambda_{0}^{-1}(\nu)\nabla_{\delta}\log p\left(Y_{j}\mid
Z_{j};\theta_{0},\delta_{0}\right)\end{array}\right)\right\\}\Rightarrow\left\\{\nu\mapsto\left(\begin{array}[]{c}U_{\theta}(i/m)\\\
V_{0}(\nu)\Lambda_{0}^{-1}(\nu)U_{\delta}(i/m)\end{array}\right)\right\\}.$

So by Slutsky’s theorem and the continuous mapping theorem, we have weak
convergence of the process $\phi_{m}$ with respect to the space of bounded
functions $f:\Delta\mapsto\mathbb{R}^{d}$ as follows

$\displaystyle\left\\{(\nu_{1},\nu_{2})\mapsto\phi_{m}(\nu_{1},\nu_{2})\right\\}_{(\nu_{1},\nu_{2})\in\Delta}$
$\displaystyle\Rightarrow\left\\{(\nu_{1},\nu_{2})\mapsto
U_{\delta}(\nu_{2})-U_{\delta}(\nu_{1})+\int_{\nu_{1}}^{\nu_{2}}\bar{V}_{0}(v)\Lambda_{0}^{-1}(v)U_{\theta}(v)dv\right\\}_{(\nu_{1},\nu_{2})\in\Delta}.$
(34)

Combining this result with Lemma B.1 and (33), the process
$\psi_{m}^{(\text{plugin})}$ converges weakly to the same limit as $\phi_{m}$.
∎

###  B.2 Asymptotics under the alternative

Suppose there is some $K^{\prime}\in(\kappa,K]$ that satisfies the following
assumptions. Assumptions B.2, B.3, and B.4 can be viewed as analogous but
simplified versions of the assumptions in Section B.1. Assumption B.5 assumes
the cumulative score process is characterized by some non-zero drift under the
alternative for some time period, even when we continually update the plugin
estimator $\hat{\theta}_{m,i}$. For example, this is likely to hold for values
of $K^{\prime}$ that are slightly larger than $\kappa$, since the plugin
estimators up to time $\lfloor mK^{\prime}\rfloor$ will not have strayed too
far from its actual value of $\theta_{0}$ prior to the changepoint.

######  Assumption B.2.

Under the alternative, suppose that

$\frac{1}{\sqrt{m}}\sum_{j=\lfloor m\kappa\rfloor}^{\lfloor
mK^{\prime}\rfloor}\nabla_{\delta}\log p\left(Y_{j}\mid
Z_{j};\bar{\hat{\theta}}\left(\frac{j-1}{m}\right),\delta_{0}\right)-\mathbb{E}_{\theta_{0},\delta_{1}\mathbbm{1}\left\\{i\geq
m\kappa\right\\}}\left[\nabla_{\delta}\log p\left(Y_{j}\mid
Z_{j};\bar{\hat{\theta}}\left(\frac{j-1}{m}\right),\delta_{0}\right)\middle|Z_{j}\right]=O_{P}(1).$

######  Assumption B.3.

Under the alternative, suppose that

$\max_{i=\lfloor m\kappa\rfloor,\cdots,\lfloor
mK^{\prime}\rfloor}\left\|\sqrt{m}\left(\hat{\theta}_{m,i}-\bar{\hat{\theta}}(i/m)\right)\right\|_{2}=O_{p}(1).$

######  Assumption B.4.

There is some $c>0$ and neighborhood $B$ that includes the set
$\left\\{\bar{\hat{\theta}}\left(v\right):v\in[\tau,K^{\prime}]\right\\}$ with
nonzero radius such that

$\sup_{\theta\in
B}\sup_{y\in\mathcal{Y},z\in\mathcal{Z}}\left\|\nabla_{\theta}\nabla_{\delta}\log
p\left(y\mid z;\theta,\delta_{0}\right)\right\|_{\infty}\leq_{p}c.$

######  Assumption B.5.

There is some $c>0$ such that

$\lim_{m\rightarrow\infty}\left\|\frac{1}{m}\sum_{j=\lfloor
m\kappa\rfloor}^{\lfloor
mK^{\prime}\rfloor}\mathbb{E}_{\theta_{0},\delta_{1}}\left[\nabla_{\delta}\log
p\left(Y_{j}\mid
Z_{j};\bar{\hat{\theta}}(t),\delta_{0}\right)\right]\right\|_{2}\geq c.$

######  Theorem B.1.

Suppose Assumptions B.2 to B.5 hold. Then under the alternative hypothesis, we
have

$\lim_{m\rightarrow\infty}\Pr\left(\exists t\in\\{m+1,\cdots,mK\\}\text{ such
that }C_{m}^{(\text{plugin})}(t)>h_{m}(t)\right)=1.$

###### Proof.

By the definition of $C_{m}^{(\text{plugin})}(t)$, suffices to prove that

$\frac{1}{\sqrt{m}}\psi_{m}^{(\text{plugin})}(\lfloor m\kappa\rfloor,\lfloor
mK^{\prime}\rfloor)=\frac{1}{\sqrt{m}}\sum_{i=\lfloor m\kappa\rfloor}^{\lfloor
mK^{\prime}\rfloor}\nabla_{\delta}\log p\left(Y_{i}\mid
Z_{i};\hat{\theta}_{m,i-1},\delta_{0}\right)$

goes to infinity. Consider the following decomposition

$\displaystyle\begin{split}\frac{1}{\sqrt{m}}\psi_{m}^{(\text{plugin})}(\lfloor
m\kappa\rfloor,\lfloor mK^{\prime}\rfloor)=&\frac{1}{\sqrt{m}}\sum_{i=\lfloor
m\kappa\rfloor}^{\lfloor
mK^{\prime}\rfloor}\operatorname{\mathbb{E}}_{\theta_{0},\delta_{1}}\left[\nabla_{\delta}\log
p\left(Y_{i}\mid
Z_{i};\bar{\hat{\theta}}\left(\frac{i-1}{m}\right),\delta_{0}\right)\right]\\\
&+\frac{1}{\sqrt{m}}\sum_{i=\lfloor m\kappa\rfloor}^{\lfloor
mK\rfloor}\left\\{\nabla_{\delta}\log p\left(Y_{i}\mid
Z_{i};\bar{\hat{\theta}}\left(\frac{i-1}{m}\right),\delta_{0}\right)-\mathbb{E}_{\theta_{0},\delta_{1}}\left[\nabla_{\delta}\log
p\left(Y_{i}\mid
Z_{i};\bar{\hat{\theta}}\left(\frac{i-1}{m}\right),\delta_{0}\right)\right]\right\\}\\\
&+R_{m}\end{split}$ (35)

The first term on the right hand side must diverge to infinity by Assumption
B.5. Per Assumption B.2, the second term is $O_{p}(1)$. Per Assumption B.3 and
Taylor’s theorem, the remainder $R_{m}$ satisfies

$\displaystyle\lim_{m\rightarrow\infty}\left\|R_{m}\right\|_{2}\leq\lim_{m\rightarrow\infty}\frac{c_{1}}{\sqrt{m}}\sup_{\\{\tilde{\theta}_{m,i-1}:i=\lfloor
m\kappa\rfloor,\cdots\lfloor mK^{\prime}\rfloor\\}\in B}\left\|\sum_{i=\lfloor
m\kappa\rfloor}^{\lfloor mK^{\prime}\rfloor}\nabla_{\theta}\nabla_{\delta}\log
p\left(Y_{i}\mid
Z_{i};\tilde{\theta}_{m,i-1},\delta_{0}\right)\left(\hat{\theta}_{m,i-1}-\bar{\hat{\theta}}\left(\frac{i-1}{m}\right)\right)\right\|_{2}.$
(36)

Combined with Assumption B.4, we have that

$\displaystyle\lim_{m\rightarrow\infty}\left\|R_{m}\right\|_{2}\leq\lim_{m\rightarrow\infty}c_{2}\sqrt{m}\max_{i=\lfloor
m\kappa\rfloor,\cdots,\lfloor
mK^{\prime}\rfloor}\left\|\hat{\theta}_{m,i-1}-\bar{\hat{\theta}}\left(\frac{i-1}{m}\right)\right\|_{2}=O_{p}(1).$
(37)

As such, the right hand side of (35) diverges to infinity, which means its
left hand side must also diverge to infinity. Thus we have our desired result.
∎

##  Appendix C Implementation details for the score-based CUSUM

There are a number of implementation decisions to make. First, the number of
sequences $B$ should be set to a value large enough such that the estimated
DCLs converge. In our simulations, we chose $B$ so that the chart statistic
for five or more bootstrapped sequences exceeded the DCL at each time step.
Next, one can maximize statistical power by tuning the shape of the alpha-
spending function. Here we simply use a linear alpha-spending function, but
future work may explore nonlinear functions instead. Finally, our theoretical
results allow for monitoring in a fully sequential manner or in batches of
observations. We found that batching had a negligible impact on detection
delay and improved both computational efficiency and convergence to the
asymptotic distribution. As such, we recommend setting the batch size so a
significant change is unlikely to occur within a batch. In our experiments,
the batch size is set to 10.

##  Appendix D Bayesian changepoint monitoring

In this section, we briefly review Bayesian changepoint monitoring and discuss
its implementation.

To determine if a structural change has occurred, the chart statistic in
Bayesian monitoring is the posterior probability of there having been a
change, i.e. $\hat{C}^{\text{bayes}}(t)=\Pr\left(\kappa\leq t\mid
Y_{1},\cdots,Y_{t},Z_{1},\cdots,Z_{t}\right)$. The operating characteristics
of a Bayesian procedure is defined by its probability of firing a false alarm
$\Pr\left(\kappa>\hat{T}^{\text{bayes}}\right)$, where
$\hat{T}^{\text{bayes}}$ is the alarm time and the probability is marginalized
over the prior. One can show that the static control limit of $1-\alpha$
indeed controls the false alarm probability at level $\alpha$, albeit
conservatively (see [Tartakovsky et al., 2014], page 318).

The performance of the Bayesian monitoring procedue is sensitive to the choice
of the prior. Given the minimax optimality of the Shiryaev-Roberts procedure
[Shiryaev, 1963], we use a modified geometric distribution for the prior of
$\kappa$, in which

$\displaystyle\pi(\kappa=t)\propto p(1-p)^{t-1}\quad\forall t=m+1,\cdots,mK$
(38)

with $p=1/mK$ and the probability of there being no changepoint is set to
$0.5$. We assume a normal prior for $\theta$ with the mean and covariance
matrix set to the results from maximum likelihood estimation on the non-
contaminated data. We also place a normal prior for $\delta$ with mean zero
and a diagonal covariance matrix. In the simulations, we set the diagonal
matrix so that the mean norm of $\delta$ in the prior is close to that of the
actual shift. In practice, such information is not known and one must rely on
prior knowledge.

In this paper, we use Hamiltonian Monte Carlo (HMC) implemented using Stan
[Carpenter et al., 2017]. While this was sufficient for running experiments,
there are a number of caveats when running HMC. First, HMC is not recommended
for long-term monitoring, because it runs posterior inference from scratch for
every new observation. A potential solution is to perform Bayesian filtering
with a Laplace approximation [McCormick et al., 2012, Feng et al., 2022b], but
existing methods have yet to be specialized for detecting a single
changepoint. Second, HMC requires the model to be fully differentiable, so it
is not guaranteed to provide valid inference for the risk shift model (13).

##  Appendix E Simulation details

For each simulation, we generate a $p^{\prime}$-random vector $X_{t}$ and
random variables $\tilde{X}_{t}$ and $U_{t}$. All these variables are drawn
independently from the uniform distribution from -1 to 1. The outcome
$Y_{t}(0)$ was generated using either (12) or (13) with
$Z_{t}=(X_{t},\tilde{X}_{t},U_{t},1)$. Treatments were assigned using one of
two models. The first is a logistic regression model with
$(\hat{f}_{t}(X_{t}),X_{t},\tilde{X}_{t},U_{t})$ with coefficients and
intercept denoted by $\gamma_{\text{LR}}$. The second sets treatment to
$A_{t}=\max(A_{t}^{(1)},A_{t}^{(2)})$, where $A_{t}^{(1)}$ and $A_{t}^{(2)}$
are generated using logistic regression models with inputs
$(\hat{f}_{t}(X_{t}),X_{t},\tilde{X}_{t},U_{t})$, with coefficients and
intercept denoted by $\gamma^{(1)}$ and $\gamma^{(2)}$, respectively. The
model parameters used to generate outcomes and treatment assignments are given
in Table 2.

Section |  Experiment(s) |  Outcome model |  Treatment model
---|---|---|---
5.1 |  CE with respect to $\hat{f}_{t}(X_{t})$ |  (12) with $\theta=(2,1,1,1,\vec{0}_{4},0,0,0)$ and $\delta=\vec{0}$ |  $\gamma_{\text{LR}}=(0.3,\vec{0}_{8},0,0,0)$ up to $t=mK/2$.
$\gamma_{\text{LR}}=(0.6,\vec{0}_{8},0,0,0)$ after $t=mK/2$.
5.1 |  CE with respect to $(\hat{f}_{t}(X_{t}),\tilde{X}_{t})$ |  (12) with $\theta=(2,1,1,1,\vec{0}_{4},1,0,0)$ and $\delta=\vec{0}$ |  $\gamma_{\text{LR}}=(0.3,\vec{0}_{8},0.1,0,0)$ prior to $t=mK/2$.
$\gamma_{\text{LR}}=(0.6,\vec{0}_{8},0.2,0,0)$ after $t=mK/2$.
5.1 |  TC with respect to $\hat{f}_{t}(X_{t})$ |  (13) with $\theta=(2,1,1,1,\vec{0}_{4},0,1,0)$ and $\delta=\vec{0}$ |  $\gamma^{(1)}=(0,\vec{0}_{8},0,1,-2)$ at all time points.
$\gamma^{(2)}=(0.2,\vec{0}_{8},0,0,0)$ up to $t=mK/2$.
$\gamma^{(2)}=(0.4,\vec{0}_{8},0,0,0)$ after $t=mK/2$.
5.1 |  TC with respect to $(\hat{f}_{t}(X_{t}),\tilde{X}_{t})$ |  (13) with $\theta=(2,1,1,1,\vec{0}_{4},1,1,0)$ and $\delta=\vec{0}$ |  $\gamma^{(1)}=(0,\vec{0}_{8},0,1,-2)$ at all time points.
$\gamma^{(2)}=(0.2,\vec{0}_{8},0.3,0,0)$ up to $t=mK/2$.
$\gamma^{(2)}=(0.4,\vec{0}_{8},0.6,0,0)$ after $t=mK/2$.
5.2 |  Model retrained using ridge-penalized logistic regression or gradient boosted trees |  (12) with $\theta=(2,1,1,\vec{0}_{47},0,0,0)$ and $\delta=\vec{0}$ |  $\gamma_{\text{LR}}=(0.5,\vec{0}_{50},0,0,0)$
5.3 |  Big and small shifts |  (12) with $\theta=(2,1,1,1,\vec{0}_{4},0,0,0)$ and $\delta=(-1.6,-0.8,-0.8,-0.8,\vec{0}_{4},0,0,0)$ or $\delta=(-1,-0.5,-0.5,-0.5,\vec{0}_{4},0,0,0)$ for the big or small shifts, respectively. |  $\gamma_{\text{LR}}=(0.15,\vec{0}_{8},0,0,0)$
5.4 |  Risks either shifted symmetrically or only among those with high-risk. |  (12) with $\theta=(2,1,1,1,\vec{0}_{4},0,0,0)$ and $\delta=(-1,-0.5,-0.5,-0.5,\vec{0}_{4},0,0,0)$ or $\delta=(-1,-0.5,-0.5,-0.5,\vec{0}_{4},0,0,-0.75)$ for symmetric- or high-risk shifts, respectively. |  No-trust: $\gamma_{\text{LR}}=(0.01,\vec{0}_{8},0,0,0)$
Calibrated-trust: $\gamma_{\text{LR}}=(1,\vec{0}_{8},0,0,0)$
Over-trust: $\gamma_{\text{LR}}=(5,\vec{0}_{8},0,0,0)$

5.5 |  Violations of the time-constant selection bias assumption are introduced at time $t^{\prime}$ |  (13) with $\theta=(2,1,1,1,\vec{0}_{4},0,0,0)$ and $\delta=(-0.1,-0.02,\vec{0}_{6},0,0,0)$ |  $\gamma^{(1)}=(0,\vec{0}_{8},0,1,-1)$ up to $t^{\prime}$.
$\gamma^{(1)}=(-0.5,\vec{0}_{8},0,1,-1)$ after $t^{\prime}$.
$\gamma^{(2)}=(0.8,\vec{0}_{8},0,0,0)$ at all time points.
Table 2:  Model parameters used to generate outcomes and treatment assignments
in the simulations.

##  Appendix F Naïve monitoring of the misclassification rate

Here we present an additional simulation that compares the false alarm rate of
score-based monitoring with a naïve CUSUM procedure that monitors the overall
misclassification rate to our proposed score-based CUSUM procedure that
monitors the conditional distribution of $Y_{t}(0)|\hat{f}_{t}(X_{t})$. We
assume conditional exchangeability with respect to $\hat{f}_{t}(X_{t})$. The
outcome is generated per (12) with $Z_{t}=X_{t}$ and
$\delta=(2,1,1,1,\vec{0}_{4},0)$. Treatment is assigned using a logistic
regression model with input as $\hat{f}_{t}(X_{t})$ with parameter
$\gamma_{\text{LR}}=(1,-0.5)$ up to time $t=200$ and
$\gamma_{\text{LR}}=(5,-2.5)$ thereafter. Rather than fitting a risk
prediction model, we fit a locked binary classifier for simplicity and
classify any observation to be positive if the predicted risk exceeds 0.7.

As shown in Figure 9, the naïve CUSUM has a substantially inflated false alarm
rate because the overall misclassification rate, without further adjustment,
is sensitive to shifts in clinician trust. On the other hand, the proposed
score-based CUSUM procedure controls the false alarm rate at the desired
level.

Figure 9: Comparison of naïve CUSUM procedure that monitors the unadjusted
overall misclassification rate versus score-based CUSUM monitoring of the
conditional distribution $Y_{t}(0)|\hat{f}_{t}(X_{t})$