# Backward Joint Model for Dynamic Prediction using Multivariate Longitudinal
and Competing Risk Data

###### Abstract

Joint modeling is a useful approach to dynamic prediction of clinical outcomes
using longitudinally measured predictors. When the outcomes are competing risk
events, fitting the conventional shared random effects joint model often
involves intensive computation, especially when multiple longitudinal
biomarkers are be used as predictors, as is often desired in prediction
problems. Motivated by a longitudinal cohort study of chronic kidney disease,
this paper proposes a new joint model for the dynamic prediction of end-stage
renal disease with the competing risk of death. The model factorizes the
likelihood into the distribution of the competing risks data and the
distribution of longitudinal data given the competing risks data. The
estimation with the EM algorithm is efficient, stable and fast, with a one-
dimensional integral in the E-step and convex optimization for most parameters
in the M-step, regardless of the number of longitudinal predictors. The model
also comes with a consistent albeit less efficient estimation method that can
be quickly implemented with standard software, ideal for model building and
diagnotics. This model enables the prediction of future longitudinal data
trajectories conditional on being at risk at a future time, a practically
significant problem that has not been studied in the statistical literature.
We study the properties of the proposed method using simulations and a real
dataset and compare its performance with the shared random effects joint
model.

Wenhao Li1, Liang Li1,∗, Brad C. Astor2, Wei Yang3, Tom H. Greene4

1 Department of Biostatistics, The University of Texas MD Anderson Cancer
Center, Houston, TX, U.S.A.

2 School of Medicine and Public Health, University of Wisconsin-Madison,
Madison, WI, U.S.A.

3 Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA,
U.S.A.

4 School of Medicine, University of Utah, Salt Lake City, UT, U.S.A.

∗ _email_<EMAIL_ADDRESS>
Keywords: Competing risks; Dynamic prediction; Joint model; Multivariate
longitudinal data; Survival analysis.

##  1 Introduction

Chronic kidney disease (CKD) is a major public health burden, affecting
approximately 14% of the adult population in the US. It is characterized by a
gradual and often irreversible deterioration of kidney function that may occur
over years (Appel et al., 2023). The stages of CKD are defined based on kidney
function, usually quantified by the glomerular filtration rate (GFR, in
mL/min/1.73 m2). GFR estimates how well the kidneys filter waste and excess
fluid from the blood. Stages 1, 2, 3A, 3B, 4 and 5 correspond to GFR in
$\geqslant 90$, $60-89$, $45-59$, $30-44$, $15-29$, and $<15$ mL/min/1.73 m2,
respectively. At stage 5, patients need to be prepared for entering End-Stage
Renal Disease (ESRD), when the kidneys stop functioning adequately and either
dialysis or kidney transplant is required to sustain life. More severe CKD is
progressively associated with many adverse outcomes, including cardiovascular
disease, stroke, heart failure and peripheral vascular disease (Appel et al.,
2023). Major causes of CKD in the US include diabetes mellitus and
hypertension. CKD, aging, and these comorbidities may result in death, though
the exact cause can be difficult to ascertain and is often multifactorial.

Predicting ESRD is important for CKD research, patient counseling, and
clinical management of the disease. Various prediction models have been
proposed (Lim et al., 2022), but they leave one or more of statistical
challenges unaddressed. First, the prediction should ideally adapt to the
ever-changing at-risk population and the potentially time-varying association
between the predictor and outcome (Li et al., 2016). The at-risk population
changes because of both aging and systematic differences between those who
progress to ESRD early and those who do not. The relative importance of
predictors may vary at various stages of the disease. Second, it is important
to have the capacity to incorporate prior clinical history into prediction, as
opposed to only the current health status. Third, the prediction of ESRD
should adjust for the competing risk by death to avoid bias and produce proper
interpretation. Fourth, not all patients are equally concerned about the near-
future risk of ESRD. Those at earlier stages have low ESRD risk and are more
interested in the future progression trajectory of GFR if ESRD or death does
not occur (GFR is undefined when a patient reaches ESRD or death, obviously).

The statistical challenges above can be handled in the framework of dynamic
prediction (Taylor et al., 2005; Yu et al., 2008; Rizopoulos, 2011). A typical
approach is to estimate the conditional joint distribution of the longitudinal
predictor variables (e.g., GFR) and survival data (e.g., the times of ESRD and
death), given the time-invariant predictors. Once this joint distribution is
estimated from the training dataset, the prediction for any subject in the
validation dataset can be calculated straightforwardly from the conditional
distribution given the longitudinal history before the prediction time.
Therefore, the key to dynamic prediction is to build an appropriate joint
model of longitudinal and survival data. Comprehensive reviews on this topic
are available (Rizopoulos, 2012; Elashoff et al., 2016; Andrinopoulou et al.,
2021), but the published statistical methods are insufficient to study the
aforementioned CKD problem.

First, most published joint modeling research presented data applications and
simulations with a very small number of longitudinal variables. The main
reason is that the likelihood function of these models is complicated and
contains an analytically intractable numerical integration for each of the $n$
subjects in the training dataset. Since the typical joint model uses shared
random effects to quantify the association between longitudinal and survival
data, the dimension of the integral increases linearly with the number of
longitudinal variables, and the number of unknown parameters in the random
effects distribution increases quadratically. In addition, the marginal log-
likelihood function is usually not concave with multiple local maxima, which
may cause non-convergence or computational instability in numerical
optimization, especially when the number of unknown parameters is large. These
are widely recognized difficulties in the implementation of such models.
(Rizopoulos, 2012; Elashoff et al., 2016) They become particularly troublesome
when we use joint models for prediction problems, which may require a large
number of longitudinal predictors for improved predictive accuracy. Second,
the aforementioned computational problem exacerbates the joint modeling with
competing risk outcomes, whose multiple survival endpoints further increase
the number of shared random effects. The limited literature on this topic
presented numerical studies with a small number of random effects (or
dimension of integrals), typically coming from one or two longitudinal
variables, and with “brute force” maximum likelihood or Bayesian computation
to tackle the model fitting problem (Elashoff et al., 2007a, b; Li et al.,
2009, 2012; Andrinopoulou et al., 2014, 2015; Rajeswaran et al., 2018). When
the number of longitudinal variables increases, it becomes increasingly
difficult to use this strategy to maintain computational stability and avoid
excessive computing time. Third, while there is previous literature on the
prediction of future longitudinal biomarker trajectories (Taylor et al., 2005;
Yu et al., 2008; Rizopoulos, 2012), this was done without conditioning on the
fact that the future survival endpoint does not occur and hence is not
suitable for our research problem (e.g., the predicted future GFR at year 3 is
meaningless if the patient dies before that time).

In this paper, we propose a new type of joint model that avoids these
limitations and is therefore suitable for the CKD research problem and beyond.
It extends the basic idea of backward joint model (BJM) (Shen and Li, 2021) to
the competing risk setting and also enables the prediction of future
longitudinal renal biomarkers in the absence of death or ESRD. The model and
estimation can be easily scaled for a large number of longitudinal variables
because the dimension of the numerical integral is always capped at 1. For
most model parameters, the log-likelihood has a concave shape, which
facilitates convergence. The proposed method also comes with a simpler version
that can be implemented quickly with standard statistical software to obtain
consistent albeit slightly less efficient estimation. The proposed BJM is also
flexible in the sense that it can use various competing risk models or
longitudinal data models as its building blocks for the prediction.

##  2 Data, Estimands and Models

###  2.1 Data and notation

The training dataset has $n$ subjects, indexed by $i=1,2,...,n$. Let
$\tilde{T}_{i}$ denote the true time to a terminal event and $C_{i}$ denote
the censoring time. All the time variables in this paper are expressed as the
time gap since baseline unless stated otherwise. The terminal event can be one
of $J$ event types, denoted by $\tilde{D}_{i}\in\\{1,2,...,J\\}$. In our
motivating data application in CKD, $J=2$, and the terminal event is either
ESRD ($j=1$) or death ($j=2$). The observed survival time is
$T_{i}=\min(\tilde{T}_{i},C_{i})$. The observed event type is
$D_{i}=\delta_{i}\tilde{D}_{i}\in\\{0,1,...,J\\}$, where the censoring
indicator $\delta_{i}=\bm{1}(\tilde{T}_{i}<C_{i})$ equals $1$ when the
terminal event is observed and $0$ when censored. Let $\bm{V}_{i}$ denote the
vector of baseline covariates. Suppose that there are $M$ longitudinally
measured continuous biomarkers, denoted by $m=1,2,...,M$. We denote the
repeated measures of the $m$-th biomarker of subject $i$ by $\bm{Y}_{mi}$,
with the corresponding vector of time points $\bm{t}_{mi}$. The number of
longitudinal biomarker measurements is $n_{mi}$, which varies by biomarker and
subject. The measurement times of a biomarker may be irregularly spaced
between baseline and $T_{i}$, and different biomarkers may be measured at
different time points. Let
$\bm{Y}_{i}=(\bm{Y}_{1i}^{T},\bm{Y}_{2i}^{T},...,\bm{Y}_{Mi}^{T})^{T}$ and
$\bm{t}_{i}=(\bm{t}_{1i}^{T},\bm{t}_{2i}^{T},...,\bm{t}_{Mi}^{T})^{T}$ collect
all the longitudinal biomarker data and their measurement times in vector
forms. We follow two conventional assumptions from the joint modeling
literature. First, the biomarker measurement times are non-informative in the
sense that the distribution of $\bm{t}_{i}$ in $(0,T_{i})$ is independent of
the rest of the data. Second, the censoring is independent in the sense that
$C_{i}$ is conditionally independent of $\bm{Y}_{i}$ and $\tilde{T}_{i}$ given
$\bm{V}_{i}$.

###  2.2 Prediction estimands

The models in this paper are motivated by the research objective, estimating
the predicted distribution of the survival and longitudinal outcomes given the
longitudinal history prior to the time of prediction. Therefore, we must
define the estimands of this prediction problem before describing the models.
We use the subscript $o$ to index a generic subject in the validation dataset.
Let $s$ be the time of prediction and $\overline{\bm{Y}_{o}(s)}$ represent the
longitudinal biomarker history prior to $s$ (i.e., all the biomarker
measurements before $s$). $\Delta$ denotes a pre-specified prediction horizon.
The prediction estimand of scientific interest is defined and derived in (1).
This is the probability of having the terminal event of type $j$
($j=1,2,...,J$) in the next $\Delta$ years (or other time units), given that
the subject is still at-risk at the prediction time $s$ ($T_{o}>s$) and given
the observed history ($\bm{V}_{o}$, $\overline{\bm{Y}_{o}(s)}$). The
derivation used the independent censoring assumption.

$\displaystyle P(s<\tilde{T}_{o}\leqslant
s+\Delta,\tilde{D}_{o}=j|\overline{\bm{Y}_{o}(s)},T_{o}>s,\bm{V}_{o})=\dfrac{P(\overline{\bm{Y}_{o}(s)},s<\tilde{T}_{o}\leqslant
s+\Delta,\tilde{D}_{o}=j|\bm{V}_{o})}{P(\overline{\bm{Y}_{o}(s)},\tilde{T}_{o}>s|\bm{V}_{o})},~{}\forall
j$ (1)

We use $P(.)$ to denote probability or probability density, when there is no
ambiguity. We are also interested in the dynamic prediction of future
longitudinal biomarker distribution. This is formulated as the estimand in (2)
for a future time $t>s$. The notation $\bm{Y}_{o}(t)$ denotes the generic
subject $o$’s longitudinal biomarker value at $t$.

$\displaystyle
P(\bm{Y}_{o}(t)=\bm{y}|\tilde{T}_{o}>t,\overline{\bm{Y}_{o}(s)},T_{o}>s,\bm{V}_{o})=\frac{P(\bm{Y}_{o}(t)=\bm{y},\overline{\bm{Y}_{o}(s)},\tilde{T}_{o}>t|\bm{V}_{o})}{P(\overline{\bm{Y}_{o}(s)},\tilde{T}_{o}>t|\bm{V}_{o})},$
(2)

This is the probability density function of the longitudinal biomarker at a
future time $t$, given that the subject is at-risk by the time of prediction
($T_{o}>s$), and the observed history
$\\{\overline{\bm{Y}_{o}(s)},\bm{V}_{o}\\}$, and conditional on the absence of
the terminal event by time $t$. A key novel aspect of this definition is that
the predicted longitudinal biomarker at a future time $t$ can only be defined
when the terminal event does not occur by that time. In our data application,
for example, the GFR does not exist if the patient has reached ESRD or death.
We may choose $t=s+\Delta$ to obtain the biomarker prediction at the horizon,
or choose $t$ to be any time between $s$ and $s+\Delta$ to obtain the
predicted longitudinal trajectory in $(s,s+\Delta)$. The notation in (2) is
designed for the most general case that we predict the multivariate density of
all the $M$ biomarkers. In our data application, we are interested in the
single most important longitudinal biomarker and similar expression can be
formulated for that situation. Given the estimated density function, various
interpretable features of the density, such as quantiles, mean and mode, can
be subsequently obtained, as illustrated in the data application.

###  2.3 Backward joint models with extrapolation or two parts

If the conditional joint distribution
$P(\bm{Y}_{o},\tilde{T}_{o},\tilde{D}_{o}|\bm{V}_{o})$ can be estimated from
the training dataset, then the prediction estimands in (1) and (2) can be
easily calculated. Since
$P(\bm{Y}_{o},\tilde{T}_{o},\tilde{D}_{o}|\bm{V}_{o})=P(\tilde{T}_{o},\tilde{D}_{o}|\bm{V}_{o})P(\bm{Y}_{o}|\tilde{T}_{o},\tilde{D}_{o},\bm{V}_{o})$,
we propose sub-models for the two conditionals. The first is a competing risk
regression with baseline covariates. The second sub-model can be a
multivariate linear mixed model when $\bm{Y}_{o}$ is continuous. Since
$\bm{Y}_{o}$ occurs prior to $\tilde{T}_{o}$, we call this approach a
“backward” joint model (BJM). This likelihood decomposition resembles the
pattern mixture model in the missing data literature but there are fundamental
differences in model formulation, estimation and scope of application, as
elucidated previously (Shen and Li, 2021). From (1) and (2), all that is
needed for prediction purposes is the conditional joint distribution
$P(\bm{Y}_{o},\tilde{T}_{o},\tilde{D}_{o}|\bm{V}_{o})$. In fact, any model
that leads to the estimation of this joint distribution can be used here. The
model does not have to reflect the data-generating mechanism in the sense that
it consists of the conditional distributions of later data given earlier data.
The data-generating mechanism of real data is usually far more complicated
than a statistical model can describe.

Both prediction estimands involve
$P(\overline{\bm{Y}_{o}(s)},\tilde{T}_{o}>t|\bm{V}_{o})=\int_{t}^{\infty}P(\overline{\bm{Y}_{o}(s)}|\tilde{T}_{o}=u,\bm{V}_{o})P(\tilde{T}_{o}=u|\bm{V}_{o})du,~{}(t\geqslant
s)$. In typical biomedical studies, $\tilde{T}_{o}$ and $C_{o}$ have finite
supports. We consider two study designs. First, there is a time $\tau_{max}$
such that $P(C_{o}>\tau_{max}|\bm{V}_{o})>0$ and
$P(\tilde{T}_{o}>\tau_{max}|\bm{V}_{o})=0$, i.e., the support of $C_{o}$ is
longer than the support of $\tilde{T}_{o}$. Both
$P(\overline{\bm{Y}_{o}(s)}|\tilde{T}_{o}=u,\bm{V}_{o})$ and
$P(\tilde{T}_{o}=u|\bm{V}_{o})$ are identifiable at every $u$ on the support
of $\tilde{T}_{o}$, and we just need a sub-model for each conditional
distribution. This is the design where the follow-up is generally longer than
the survival time of study subjects. Second, there is a time $\tau_{max}$ such
that $P(\tilde{T}_{o}>\tau_{max}|\bm{V}_{o})>0$ and
$P(C_{o}>\tau_{max}|\bm{V}_{o})=0$ for at least some subjects in the
population, i.e., the support of $C_{o}$ is shorter than the support of
$\tilde{T}_{o}$. $P(\overline{\bm{Y}_{o}(s)}|\tilde{T}_{o}=u,\bm{V}_{o})$ and
$P(\tilde{T}_{o}=u|\bm{V}_{o})$ are not always identifiable from the data when
$u>\tau_{max}$. This is the design where the follow-up is not long enough to
estimate the entire survival function. To solve this problem, we write
$P(\overline{\bm{Y}_{o}(s)},\tilde{T}_{o}>t|\bm{V}_{o})$ as

$\int_{t}^{\tau_{max}}P(\overline{\bm{Y}_{o}(s)}|\tilde{T}_{o}=u,\bm{V}_{o})P(\tilde{T}_{o}=u|\bm{V}_{o})du+P(\overline{\bm{Y}_{o}(s)}|\tilde{T}_{o}>\tau_{max},\bm{V}_{o})P(\tilde{T}_{o}>\tau_{max}|\bm{V}_{o}).$
(3)

A sub-model for survival data can identify both
$P(\tilde{T}_{o}=u|\bm{V}_{o})$, $t<u\leqslant\tau_{max}$, and
$P(\tilde{T}_{o}>\tau_{max}|\bm{V}_{o})$. Two sub-models for longitudinal data
are needed to identify
$P(\overline{\bm{Y}_{o}(s)}|\tilde{T}_{o}=u,\bm{V}_{o})$,
$t<u\leqslant\tau_{max}$ and
$P(\overline{\bm{Y}_{o}(s)}|\tilde{T}_{o}>\tau_{max},\bm{V}_{o})$ for those
with $\tilde{T}>\tau_{max}$. These subjects with $\tilde{T}>\tau_{max}$ are
referred to as long-term survivors (LTS; Wang et al, 2020). In real data
analysis, we can empirically use the estimated marginal distribution of
$\tilde{T}$ to decide whether the extra longitudinal sub-model for the LTS is
needed. If the Kaplan-Meier curve of $\tilde{T}$ descends to near $0$, that
suggests that the LTS group is likely very small. We call it the BJM with
extrapolation (BJM-EX). The extrapolation means that the estimated
distribution $P(\tilde{T}|\bm{V})$ extrapolates beyond the support of the data
sample. Alternatively, the LTS sub-model is needed and we call it the two-part
BJM (BJM-TP). For BJM-TP, the prediction time $s<\tau_{max}-\Delta$ because
there is no observed terminal event beyond $\tau_{max}$ and the model does not
extrapolate beyond the data range. For BJM-EX, there is no such restriction
because the model extrapolates.

###  2.4 Sub-models for competing risk outcomes

Any competing risk regression model with time-invariant covariates can be
fitted to $\\{T_{i},D_{i},\bm{V}_{i};i=1,2,...,n\\}$ and estimate
$P(\tilde{T}_{o},\tilde{D}_{o}|\bm{V}_{o})$. Common choices include the cause-
specific hazard models, the Fine-Gray sub-distribution hazard model, or the
vertical modeling (Nicolaie et al., 2010), among others. We used vertical
modeling in the numerical studies because it directly models
$P(\tilde{T}_{i}|\bm{V}_{i})$, which is convenient for predicting longitudinal
data in (2). We compared the BJMs with vertical modeling and with cause-
specific hazard models, and the predictive performances were similar in our
data application.

The likelihood of a vertical model for
$P(\tilde{T}_{i},\tilde{D}_{i}|\bm{V}_{i})$ can be factorized as the
likelihood product of $P(\tilde{T}_{i}|\bm{V}_{i})$ and
$P(\tilde{D}_{i}|\tilde{T}_{i},\bm{V}_{i})$. The former can be formulated as a
semi-parametric Cox model, a parametric Weibull model, or any other survival
regression model that fits the data well. The latter is a multinomial
regression model, according to Nicolaie et al. (2010). Specifically,
$P(\tilde{D}_{i}=j|\tilde{T}_{i}=t,\bm{V}_{i})=\dfrac{\exp(\bm{B}(t)^{T}\bm{\xi}^{(\tilde{T})}_{j}+\bm{V}_{i}^{T}\bm{\xi}^{(\bm{V})}_{j})}{\sum_{j=1}^{J}\exp(\bm{B}(t)^{T}\bm{\xi}^{(\tilde{T})}_{j}+\bm{V}_{i}^{T}\bm{\xi}^{(\bm{V})}_{j})},j=1,...,J.$.
The $\bm{B}(t)$ is a vector of pre-specified regression spline basis to
accommodate the nonlinear effect of $\tilde{T}_{i}$ on the event types. When
there are two event types, this is a logistic regression with partially linear
specification for the covariates $\tilde{T}_{i}$ and $\bm{V}_{i}$. The
estimation of vertical modeling follows Nicolaie et al. (2010).

###  2.5 Sub-models for longitudinal outcomes

For the BJM-EX, $P(\bm{Y}_{o}|\tilde{T}_{o},\tilde{D}_{o},\bm{V}_{o})$ is
modeled by a multivariate linear mixed model, with $T_{i}$, $D_{i}$,
$\bm{V}_{i}$, $\bm{t}_{i}$ and their interactions as fixed-effect covariates
and with random effects (intercept, or intercept and slope) to account for the
longitudinal correlation and the correlation among different biomarkers. The
sub-model of the $m$-th longitudinal biomarker ($m=1,2,...,M$) can be
formulated in the following general form.

$Y_{mi}(t)=\sum_{l=0}^{L}\vartheta_{ml}(\bm{V}_{i},\phi(\tilde{T}_{i}),\tilde{D}_{i})g_{l}(t)+\bm{u}_{mi}(t)^{T}\bm{b}_{mi}+\epsilon_{mi}(t)~{}~{},~{}~{}t\in[0,\tilde{T}_{i}).$
(4)

$Y_{mi}(t)$ is subject $i$’s $m$-th biomarker, measured at time $t$. In the
training data, $t\in\bm{t}_{mi}$. The first term describes the conditional
mean trajectory (CMT), given baseline covariates ($\bm{V}$), terminal event
time ($\tilde{T}$) and type ($\tilde{D}$). The CMT is defined for subjects
with the same trajectory lengths ($\tilde{T}$), because it does not lead to
proper interpretation to average the longitudinal biomarker trajectories of
different lengths when the follow-up duration is correlated with longitudinal
data, a defining feature of joint modeling problems that is sometimes called
informative dropout. The second term includes random effects that describe the
subject-specific, mean-zero deviation from the CMT. The random intercept
($\bm{u}_{mi}(t)=1$), random intercept plus slope
($\bm{u}_{mi}(t)=(1,t)^{T}$), or other nonlinear random effect terms can be
used. The $g_{l}(t),l=0,1,...,L$ are basis functions that characterize the
shape of the CMT. For example, the linear CMT has $L=1$, $g_{0}(t)=1$, and
$g_{1}(t)=t$. The coefficients of the basis functions,
$\vartheta_{ml}(\bm{V}_{i},g_{1}(\tilde{T}_{i}),\tilde{D}_{i})$, are
interpreted as the intercept, slope, etc., allowing the CMT to vary with
$\bm{V}_{i}$, $\tilde{T}_{i}$ and $\tilde{D}_{i}$. For example,
$\vartheta_{m0}(\bm{V}_{i},\phi(\tilde{T}_{i}),\tilde{D}_{i})=\bm{\beta}_{m1}^{T}\bm{V}_{i}+\nu_{m0}\tilde{D}_{i}+\nu_{m2}\phi(\tilde{T}_{i})+\nu_{m4}\phi(\tilde{T}_{i})\tilde{D}_{i}$,
$\vartheta_{m1}(\bm{V}_{i},\tilde{T}_{i},\tilde{D}_{i})=\bm{\beta}_{m2}^{T}\bm{V}_{i}+\nu_{m1}\tilde{D}_{i}+\nu_{m3}\phi(\tilde{T}_{i})+\nu_{m5}\phi(\tilde{T}_{i})\tilde{D}_{i}$.
The pre-specified transformation $\phi(.)$, e.g. $\log$, may be used to reduce
the skewness of the time variable $\tilde{T}_{i}$.

The residuals $\varepsilon_{mi}(t)$ are independent noises with mean zero and
variance $\sigma_{e,m}^{2}$. Extensions to residuals with auto-correlation are
straightforward. The random effects from all $M$ biomarkers
$\bm{b}_{i}=(\bm{b}_{1i}^{T},\bm{b}_{2i}^{T},...,\bm{b}_{Mi}^{T})^{T}$ follows
a multivariate normal distribution with mean zero and covariance matrix
$\bm{\Omega}$. When model (4) has random intercept and slope, $\bm{\Omega}$ is
a $2M\times 2M$ matrix that is not restricted to be block diagonal, allowing
the biomarkers to be correlated.

The BJM-TP has an extra longitudinal sub-model for the LTS:

$Y_{mi}(t)=\sum_{l=0}^{L}\vartheta_{ml}^{e}(\bm{V}_{i})g_{l}(t)+\bm{u}_{mi}^{e}(t)^{T}\bm{b}^{e}_{mi}+\epsilon_{mi}^{e}(t)~{}~{},~{}~{}t\in[0,\tau_{max})~{}\mathrm{and}~{}\tilde{T}_{i}>\tau_{max}.$
(5)

The superscript $e$ symbolizes quantities in the “extra” LTS model. The model
formulation and interpretation are similar to (4) but with a few differences.
First, model (5) does not involve $\tilde{T}_{i}$ and $\tilde{D}_{i}$ in the
covariates because the LTS includes all subjects who survive beyond
$\tau_{max}$ regardless of their actual $\tilde{T}_{i}$ and $\tilde{D}_{i}$.
Second, since the LTS is a heterogeneous group, models (5) and (4) cannot have
the same random effect variance or residual variance. Presumably, the
variances are larger in model (5).

##  3 Estimation

The BJM’s likelihood factorization
$P(\bm{Y}_{o},\tilde{T}_{o},\tilde{D}_{o}|\bm{V}_{o})=P(\tilde{T}_{o},\tilde{D}_{o}|\bm{V}_{o})P(\bm{Y}_{o}|\tilde{T}_{o},\tilde{D}_{o},\bm{V}_{o})$
suggests a two-step estimation procedure. In the first step, we fit the
competing risk regression model (Section 2.4). Many standard software packages
are available. In the second step, we estimate
$P(\bm{Y}_{o}|\tilde{T}_{o},\tilde{D}_{o},\bm{V}_{o})$ from the longitudinal
sub-models in Section 2.5 and this is the focus of this section. We propose
two approaches. Section 3.1 presents a simple method that can be conveniently
implemented with standard linear mixed model software. This method only uses
training data subjects whose terminal event times are observed and those still
at-risk by $\tau_{max}$. Since many censored data are not used, we call this
approach the complete-case analysis (CCA). A remarkably useful feature of the
BJM is that the CCA produces consistent estimation, albeit less efficient. In
contrast, typical joint models of longitudinal and survival data would require
specialized software and they have limited functionality and flexibility. The
CCA of BJM enables quick and convenient exploratory analysis, model building,
model assumption diagnosis, and preliminary prediction accuracy assessment.
Once the model specification is completed, we can use the fully likelihood-
based analysis in Section 3.2 to obtain efficient estimators and the final
results of the prediction model and prediction accuracy.

###  3.1 Complete-case analysis (CCA)

To see why the parameters of longitudinal sub-model can be consistently
estimated with the complete cases in the training data (i.e., $D_{i}\neq 0$),
we observe that, for any $j=1,2,...,J$,
$P(\bm{Y}|\tilde{T}=t,\tilde{D}=j,\bm{V})=P(\bm{Y}|\tilde{T}=t,C>t,\tilde{D}=j,\bm{V})=P(\bm{Y}|T=t,D=j,\bm{V})$.
The simple derivation is justified by the conditional independence of
censoring. It holds at any $t\in(0,\infty)$ for the BJM-EX, and any
$t\in(0,\tau_{max})$ for the BJM-TP. $P(\bm{Y}|T=t,D=j,\bm{V})$ can be
estimated via a linear mixed model according to the specification in (4).
Additionally,
$P(\bm{Y}|\tilde{T}>\tau_{max},\bm{V})=P(\bm{Y}|\tilde{T}>\tau_{max},C>\tau_{max},\bm{V})=P(\bm{Y}|T>\tau_{max},\bm{V})$
in the BJM-TP. Hence, model (5) can be estimated from the subset of training
data with $T_{i}>\tau_{max}$. Consistent estimation of CCA has been
demonstrated in previous research on models with a similar likelihood
factorization (Li et al., 2018; Wang et al., 2020; Shen and Li, 2021; Wang et
al., 2023).

###  3.2 EM algorithm for the analysis of all data

This section describes the EM algorithm that calculates the maximum likelihood
estimators (MLE) of the parameters in the longitudinal sub-models. Since it
uses both censored and non-censored data, the estimation is more efficient
than the CCA (Li et al., 2018; Wang et al., 2020; Shen and Li, 2021). The
censored $\tilde{T}_{i}$ and $\tilde{D}_{i}$ are treated as missing data in
the complete data log-likelihood. Since the parameters in the competing risk
sub-model are estimated consistently in Section 2.4 and are fixed at their
estimated values during the EM iterations, the method in this section is a
type of pseudo maximum likelihood estimation (Gong and Samaniego, 1981).

Let $\bm{\Theta}$ denote the parameters of the longitudinal sub-model. With
BJM-TP, $\bm{\Theta}$ includes additionally the parameters $\bm{\Theta}^{e}$
from the LTS sub-model (5). The estimated parameters from the survival sub-
model are denoted by $\hat{\bm{\Theta}}_{surv}$. Since
$\hat{\bm{\Theta}}_{surv}$ is fixed throughout the EM iterations, we omit it
for the formula below. In the following, we mainly describe the EM algorithm
of the BJM-EX, but will point out the changes needed when the BJM-TP is used.
The complete data log-likelihood is $\sum_{i=1}^{n}\left\\{\log
P(\bm{Y}_{i}|\tilde{T}_{i},\tilde{D}_{i},\bm{V}_{i};\bm{\Theta})+\log
P(\tilde{T}_{i},\tilde{D}_{i}|\bm{V}_{i})\right\\}$. The observed data of all
subjects in the training dataset is collected in
$\\{\bm{Y},\bm{T},\bm{\delta},\bm{D},\bm{V}\\}$. The EM algorithm uses the
consistent estimators from the CCA as the initial values in the first EM
iteration. At the beginning of the $q$-th iteration ($q=1,2,...$), the
estimated parameter is denoted by $\hat{\bm{\Theta}}^{(q)}$.

In the E-step of the $q$-th iteration, we calculate the conditional
expectation of the complete data log-likelihood as in (6). This expectation is
taken over the conditional distribution of the missing data (censored
$\tilde{T}_{i}$ and $\tilde{D}_{i}$) given the observed data, with parameters
in this distribution set at $\hat{\bm{\Theta}}^{(q)}$.

$\displaystyle~{}\sum_{i=1}^{n}E_{(\tilde{T}_{i},\tilde{D}_{i})}\left\\{\log
P(\bm{Y}_{i}|\tilde{T}_{i},\tilde{D}_{i},\bm{V}_{i};\bm{\Theta}^{(q+1)})+\log
P(\tilde{T}_{i},\tilde{D}_{i}|\bm{V}_{i})\Big{|}\bm{Y},\bm{T},\bm{\delta},\bm{D};\hat{\bm{\Theta}}^{(q)}\right\\}$
(6) $\displaystyle\propto$
$\displaystyle~{}\sum_{i=1}^{n}\sum_{j=1}^{J}\int_{\tilde{T}_{i}=T_{i}}^{\infty}\left\\{\log
P(\bm{Y}_{i}|\tilde{T}_{i},\tilde{D}_{i}=j,\bm{V}_{i};\bm{\Theta}^{(q+1)})\right\\}P(\tilde{T}_{i},\tilde{D}_{i}=j|\bm{Y},\bm{T},\bm{\delta},\bm{D};\hat{\bm{\Theta}}^{(q)})d\tilde{T}_{i}$

There is always a one-dimensional integral in (6), regardless of how many
longitudinal biomarker variables are involved. This integral is needed only
for censored subjects. In contrast, the widely used shared random effects
joint model (Rizopoulos, 2012) has $n$ integrals in its likelihood, and the
dimension of the integral increases with the number of longitudinal variables.
This feature is the key reason that BJM is computationally much easier. The
integration in (6) goes from $T_{i}$ to $\infty$. With BJM-TP, this
calculation is split into two parts, analogously to (3). The first integrates
from $T_{i}$ to $\tau_{max}$ using model (4) and the second uses the LTS model
(5).

The integral in (6) can be calculated with the trapezoidal rule. If $T_{i}$ is
a censoring time, we assume that $\tilde{T}_{i}$ falls in one of a series of
pre-specified connected intervals, denoted by
$\\{\mathbb{L}_{ik}=[a_{ik},b_{ik});k=1,...,K_{i}\\}$. The intervals are
connected (e.g., $b_{ik}=a_{i(k+1)}$) and starts from $T_{i}$
($a_{i1}=T_{i}$). With BJM-EX, we can set the ending time of the last interval
at a very large survival time (e.g., 100 year old). With BJM-TP, the last
interval is $[\tau_{max},\infty)$ and the LTS model is used on this interval.
A moderate $K_{i}$ can produce accurate results when the residual survival
beyond $T_{i}$ is not very large or when the integrand is not excessively
wiggly (Li et al., 2018). The target function (6) of BJM-EX then becomes:

$\displaystyle\sum_{i=1}^{n}\sum_{k=1}^{K_{i}}\sum_{j=1}^{J}\left\\{\log
P(\bm{Y}_{i}|\tilde{T}_{i}=\dfrac{a_{ik}+b_{ik}}{2},\tilde{D}_{i}=j,\bm{V}_{i};\bm{\Theta}^{(q+1)})\right\\}P(\tilde{T}_{i}\in\mathbb{L}_{ik},\tilde{D}_{i}=j|\bm{Y},\bm{T},\bm{\delta},\bm{D};\hat{\bm{\Theta}}^{(q)}).$
(7)

With BJM-TP, the target function is the same except that the terms
corresponding to the last interval ($\mathbb{L}_{iK_{i}}$) are replaced by
$\left\\{\log
P(\bm{Y}_{i}|\tilde{T}_{i}\in\mathbb{L}_{iK_{i}},\bm{V}_{i};\bm{\Theta}^{(q+1)})\right\\}P(\tilde{T}_{i}\in\mathbb{L}_{iK_{i}}|\bm{Y},\bm{T},\bm{\delta},\bm{D};\hat{\bm{\Theta}}^{(q)})$.

In the M-step, (7) is maximized with respect to $\bm{\Theta}^{(q+1)}$,
producing an updated estimator $\hat{\bm{\Theta}}^{(q+1)}$ that will enter the
next ($q+1$-th) iteration. The only terms containing the unknown parameters
are $\log
P(\bm{Y}_{i}|\tilde{T}_{i}=(a_{ik}+b_{ik})/2,\tilde{D}_{i}=j,\bm{V}_{i};\bm{\Theta}^{(q+1)})$.
These are log of multivariate normal density functions, which are quadratic
forms of the regression coefficients in the longitudinal sub-models.
$P(\tilde{T}_{i}\in\mathbb{L}_{ik},\tilde{D}_{i}=j|\bm{Y},\bm{T},\bm{\delta},\bm{D};\hat{\bm{\Theta}}^{(q)})$
is calculated in the E-step and treated as fixed in the maximization of the
M-step. Therefore, the regression coefficients have closed-form solutions (the
variance-covariance matrices of the random effects do not; see below). This
feature, in conjunction with using the consistent estimators from the CCA as
the initial value of the EM algorithm, speeds up the computation considerably
and improves the robustness of the iterations. We have observed a 100%
convergence rate with the EM and the computing time is much shorter compared
to the shared random effects joint model.

The number of regression coefficients in the BJM increases linearly with the
number of longitudinal biomarkers. However, the number of parameters in the
variance-covariance matrix of the random effects goes up quadratically, which
poses a challenge for any joint models that aim to accommodate multivariate
longitudinal data. However, the BJM has a special feature that helps avoid
this computational burden. During the EM algorithm, we fix the variance-
covariance matrix of the random effects and residuals at their consistently
estimated values from the CCA. This approach is justified by the pseudo
maximum likelihood estimation (Gong and Samaniego, 1981).

###  3.3 Variance estimation

We use the bootstrap for variance estimation. Specifically, we randomly sample
the subjects in the training data with replacement, creating a bootstrap
dataset. The bootstrap dataset is put through the BJM estimation procedure to
produce point estimators of the model parameters. This process is repeated 200
times, and the variances of the estimated model parameters can be calculated
as the sample variance of the 200 point estimators. For a subject in the
validation dataset, the prediction (1) or (2) is a function of the estimated
BJM model parameters and the subject’s data prior to the prediction time. The
variance of the prediction can be calculated by Delta-method to account for
the sampling variability of the estimated BJM model. As an alternative
approach that avoids the mathematical derivation of the Delta-method, we can
also calculate the prediction using each of the 200 estimated model parameters
from the bootstrap. The resulting sample variance quantifies the variance of
the prediction. Since the estimated model parameters are approximately normal,
a Wald-type confidence interval can be obtained from the estimated variance.

###  3.4 Prediction accuracy measures

We follow the previous literature in choosing the prediction accuracy
measures. For predicting competing risk outcomes, we use the time-dependent
Brier score (BS) and the time-dependent AUC (Definition A in Wu and Li, 2018),
both designed for competing risk data. For predicting the future longitudinal
biomarkers, we use RMSE (root mean squared error). Since our data application
is CKD research, we also use two commonly used measures in that field for
predicting GFR, the P30 and P50 (Levey, 1999). They are the proportions of
subjects for whom the predicted GFR falls within 30% or 50% of the observed
GFR at the horizon. Higher AUC, lower BS, lower RMSE, higher P30 and P50
suggest better predictive accuracy. When the GFR measurements are irregularly
spaced and there is no GFR data at the horizon, the closest measurement prior
to the horizon is used.

##  4 Data Application

We illustrate the proposed BJM methodology by analyzing data from the African
American Study of Kidney Disease and Hypertension (AASK) Cohort Study. AASK is
a longitudinal cohort study of 1,094 CKD subjects. Detailed descriptions of
AASK data were reported previously (Li et al., 2016; Wu et al., 2020). By the
end of the follow-up, 318 subjects (29%) reached ESRD and 176 subjects (16%)
died before reaching ESRD. We did not consider survival after ESRD in a semi-
competing risk analysis because the patient experience and treatments were
substantially different after reaching ESRD. We used two baseline predictors,
age and gender. We considered three longitudinally measured predictive
biomarkers, the estimated glomerular filtration rate (eGFR) using the MDRD
equation, 24-hour urine protein to creatinine ratio (UPCr), and serum albumin.
These biomarkers have known clinical relationships with the survival outcomes
and hence are ideal choices as predictors.

Web Figure S1 displays the Kaplan-Meier curve of the composite endpoint of
ESRD or death, whichever occurs first, and the cumulative incidence functions
of ESRD and death separately. The distribution of the composite event in this
population can be estimated up to 12 years, with more than 40% of the events
expected beyond that time. This supported the use of BJM-TP. For comparison
purposes, we also presented results from BJM-EX. We set $\tau_{max}=12$ years,
the maximum observed event time in the data. We experimented with both
vertical modeling and cause-specific hazard modeling as the survival sub-
model, and they produced similar prediction accuracy. The reported results
here came from BJM with vertical modeling. Figure 3 displays the CMTs of eGFR,
which are approximately linear and suggest a clear association with the time
of the composite endpoint and event types, with an interaction between the
two.

Including more predictive longitudinal biomarkers improves the prediction
accuracy. The BJM remains computationally convenient in this situation. To
illustrate the benefit of using more longitudinal predictors, we considered
BJMs with one (BJM-1, with eGFR), two (BJM-2, with eGFR and UPCr), or three
(BJM-3, with eGFR, UPCr, and albumin) longitudinal biomarkers. Our EM
algorithm achieved 100% convergence rates in these analyses. Since this paper
is mainly about prediction, we focus on prediction results here, and defer the
extensive details on model specification, parameter estimation and
interpretation, and sensitivity analysis to the Web Appendix. With BJM-3, the
extrapolation model (EX) has 48 parameters and the two-part model (TP) has 72
parameters. The BJM includes a competing risk regression model and a linear
mixed model for each longitudinal parameter. These are conventional
statistical models in which all the parameters are well identified. The number
of parameters characterizing the correlation of multiple longitudinal
biomarkers increases quadratically with the number of biomarkers, but these
parameters are inevitable in any joint models with multivariate longitudinal
biomarkers. The TP model has 50% more parameters because it uses two parts to
improve flexibility without extrapolation, but all the parameters in its LTS
sub-model are identified.

We used five-fold cross-validation to evaluate the prediction accuracy at a
horizon of 3 years (Table 1). Additional results under 1-year horizon are
reported in the Web Appendix. In predicting ESRD and death, Table 1 shows that
BJMs with more longitudinal biomarkers generally performed better than those
with less, with BJM-3-TP performing similarly to or better than all other
methods in both AUC and BS. The TP methods generally performed better than the
EX methods, likely because it is less prone to model misspecification when the
tail distribution of the composite event cannot be estimated. The performance
of BJM-3-EX rose on par with their counterpart TP model (BJM-3-TP), while
BJM-2-EX and BJM-1-EX did not. The performance of EX models is more difficult
to speculate about probably because the accuracy of the extrapolation is
unknown and cannot be verified from data. Nonetheless, BJM-3-EX achieved
similar predictive accuracy for ESRD and death as BJM-3-TP did but with fewer
parameters, making it an empirically useful model for prediction purposes. In
predicting the eGFR, TP methods substantially outperformed EX methods. The BS
appears to be more sensitive to between-method differences than the AUC,
possibly because the latter is rank-based. The incidence of ESRD can be
predicted with very high AUC because the treatment decision of dialysis or
renal transplantation is strongly influenced by the patient’s eGFR.

For comparison purposes, we fitted the widely-used shared random effects joint
model (SJM) using the JMBayes2 package in R. The model included all three
longitudinal biomarkers (i.e., denoted by SJM-3 in Table 1), whose
trajectories were modeled by linear mixed models with subject-specific random
intercepts and random slopes, and these random effects were shared with the
cause-specific hazard models as covariates. The BJMs outperformed the SJM in
predicting ESRD and death. The results of eGFR prediction are mixed. However,
it seems that the SJM does not condition on the absence of the composite event
when predicting eGFR (Section 7.2 of Rizopoulos, 2012), which makes this
comparison difficult. The BJM is computationally much faster than the SJM. To
complete the five-fold cross-validation on a personal computer with 2.9 GHz
CPU and 32 GB RAM, SJM-3 took 3 hours whereas BJM-3-TP with EM algorithm took
13 minutes. The CCA of BJM can implemented with standard statistical software
and the computing time is almost negligible.

Figure 4 depicts individual dynamic predictions of ESRD and death for two AASK
patients. The patient on the left (A) had stable eGFR before year 3 but went
into a steady decline from year 3 to year 8. The UPCr also increased over
time, indicating increasing kidney damage. The predicted future probability of
ESRD increases sharply in response. This patient eventually reached ESRD. The
patient on the right had stable biomarkers and the predicted ESRD risk is low.
The mortality risk, however, increases over time, presumably due to both aging
and low ESRD risk (a competing risk phenomenon). Figure 5 illustrates two
useful graphical presentations of the joint prediction of competing risk and
longitudinal outcomes. The first uses the quantile curves to visualize how the
predicted eGFR distribution conditional on absence of ESRD or death varies
over time, in addition to the predicted cumulative incidence functions of ESRD
and death. This patient’s eGFR was within a stable range between the two
prediction times, and the predicted ESRD was similar as a result. However, we
can observe a subtle decrease in the predicted eGFR in the second prediction,
perhaps due to the lower observed eGFR just prior to the prediction and
increased UPCr. The cumulative incidence of death increases due to time
elapsing and aging. The second graphical presentation is clinically more
interpretable (Hu et al., 2012), where the future probabilities of being in
each CKD stage and death are visualized. From year 2 to year 4, the eGFR
stabilized around 30, which triggered the model to lower the probability of
ESRD but enlarge the probability area corresponding to CKD stage 4 (eGFR
between 15 and 29). The probability of death increased slightly due to time
elapsing and aging.

We performed a sensitivity analysis to study whether the prediction accuracy
is affected by the specifications for random effects in (4), random intercept
only or both random intercept and slope, and the choices of $\phi(.)$ in that
model, $\phi(t)=\log(t)$ or $\phi(t)=t$. The random slope adds flexibility to
the individual-level longitudinal trajectories within each survival stratum
(i.e., given $\tilde{T}$) but also adds more parameters to the variance matrix
for the random effects. A log transformation may fit the data better when
$\tilde{T}$ is skewed. Web Table S4 shows that the prediction accuracy is not
sensitive to these alternative formulations. We speculate that the survival
time in the AASK data is not very skewed (Web Figure S1). For this reason, we
used the simpler model specification of random intercept and no log
transformation.

##  5 Simulation

We conducted simulations to study the finite sample performance of the
proposed estimation procedure for both TP and EX models, and to compare the
relative efficiency of CCA vs. EM algorithm. The CCA and EM methods both
produce consistent estimation, but CCA is easier to use, making it worthwhile
to assess its efficiency cost. We already compared the predictive performance
of BJM and SJM in the AASK data. We do not compare them in the simulations
because they operate under different modeling assumptions and it is impossible
to simulate data so that both are correctly specified. Within the likelihood
framework, consistent and efficient estimation of the BJM guarantees optimal
predictive accuracy.

We conducted and reported simulation studies for the BJM-EX and BJM-TP
separately. In each study, we simulated 300 independent Monte Carlo datasets
of $n$ subjects from the BJM-EX model (or BJM-TP model), $n=300$ or $1,000$,
analyzed each dataset using both the CCA and EM algorithm, and calculated the
bias of estimated model parameters and the relative efficiency between CCA and
EM. We simulated one baseline covariate and three longitudinal biomarkers. The
longitudinal sub-models had a linear CMT and random intercept. The competing
risk model is a vertical model with a Weibull regression for the time to
composite event. The censoring rate is 40% in the BJM-EX data. For BJM-TP, the
overall censoring rate is 40% and 20% subjects were censored at $\tau_{max}$.
The BJM-EX and BJM-TP were fit to their respective datasets so that we can
evaluate the estimation methods under the correct model specification. The CCA
was used to obtain the initial values for the EM. All EM algorithms converged
under the criterion that all longitudinal sub-model parameters differed by no
more than $10^{-4}$ in consecutive iterations. Further details of the data-
generating schemes are in the Web Appendix. Figure 1 shows that the finite-
sample biases of the EM estimators are generally small and drop to nearly zero
when the sample size increases from 300 to 1,000, supporting consistent
estimation. The EM algorithm has notable efficiency improvement over the CCA,
especially among the LTS sub-model parameters, probably because the CCA only
uses the subjects at-risk by $\tau_{max}$ to estimate the LTS, which could
have a small sample size.

##  6 Discussion

While the idea of backward joint modeling has been proposed for dynamic
prediction (Shen and Li, 2021) and non-dynamic prediction (Li et al., 2018)
problems, this paper is novel in several aspects: predicting competing risk
outcomes, predicting the future distribution of longitudinal outcome
conditional on being at risk, dual approaches (EX and TP) designed for
different scenarios. The estimation algorithm of BJM ensures that model
fitting remains computationally tractable with a large number of longitudinal
variables. This is a common issue in prediction problems, but there is
currently no research on joint modeling that guarantees computationally
feasible algorithms. We compared the prediction accuracy of BJM and the widely
used SJM in our data application, and BJM performed notably better. While a
generalizable conclusion cannot usually be drawn from one data analysis, and
it is difficult to compare the prediction accuracy in simulations because BJM
and SJM operate under different modeling assumptions, there are two features
in the BJM that help improve its predictive performance. The first is model
formulation. The SJM was proposed decades ago for explaining the association
between longitudinal and survival data, or for informative dropout problems.
It was recently repurposed for use in dynamic prediction problems (Rizopoulos,
2011). In contrast, the BJM was designed from the ground up for dynamic
prediction. As shown in Section 2, we first proposed the predictive estimand
in (1) and (2) and then proposed a sub-model for each component in that
estimand. Any appropriate sub-models for those components can be used in
combination to improve prediction accuracy in various contexts, an interesting
topic for future research. The second is model estimation. The SJM must be
estimated by specialized software, which often offers limited flexibility. The
lengthy computing time may also limit the scope of exploration. The BJM is
much faster to estimate, and the model formulation can even be explored in CCA
with standard software. We believe that the ease of estimation accelerates the
development of BJM variants for various research problems.

Another approach for dynamic prediction is the landmark modeling (van
Houwelingen and Putter, 2011; Zheng and Heagerty, 2005). The joint modeling
has many advantages over landmark modeling, including efficient likelihood-
based estimation and easily incorporating the longitudinal history into the
prediction. However, the landmark modeling is easy to implement, often with
standard software, and it can handle large number of predictor variables (Li
et al., 2023). With BJM, we have made joint modeling as computationally
efficient, scalable, and easy to develop as landmark modeling. In addition, we
are not aware of any landmark modeling approach that can predict future
longitudinal data distribution. It is not straightforward with landmark
modeling because the prediction of longitudinal data must be made conditional
on the absence of a future clinical event, which requires modeling their joint
distribution.

A limitation of this paper is that the longitudinal sub-models are parametric.
The BJM models the CMT, not the individual-level trajectories. The CMT usually
has less curvature and hence a parametric model may suffice. Model checking
can also be done conveniently using CCA. Nonetheless, it is difficult to
ensure that many longitudinal biomarkers are all modeled adequately with
parametric models. Future work is needed to improve the flexibility of the
fixed effect part with semi-parametric or non-parametric specification. The
random effect part can also be improved by using penalized splines for
subject-specific trajectory deviation from the CMT, though that would cause a
substantial increase in the number of parameters in the covariance matrix of
spline coefficients. However, this difficulty challenges all joint models with
multivariate nonlinear longitudinal trajectories. Recent work suggested that,
even if the random effect part is misspecified, the rest of the BJM model can
be estimated consistently (Wang et al., 2023). In the sensitivity analysis
(Web Table S4), we found that adding a random slope causes little change in
the prediction accuracy. Since the BJM characterizes the association between
survival outcome and longitudinal biomarkers through the fixed effects (as
opposed to using random effects, such as in the SJM; Shen and Li, 2021), it
remains to be studied whether the BJM is robust to misspecification of the
random effect part. Another limitation of the BJM is that it sometimes
requires an extra LTS sub-model. SJM does not have this need because it uses a
time-varying hazard function. While the LTS can be easily estimated, it
increases the number of parameters by about 50%. At a chosen $\tau_{max}$, the
percent of subjects contributing to the LTS estimation is
$P(\tilde{T}>\tau_{max})P(C>\tau_{max})$. The censoring distribution may
depend on staggered study entry, administrative censoring, and dropout
patterns. When $P(C>\tau_{max})$ is low, we may need to decrease $\tau_{max}$
to ensure enough sample size for the LTS estimation. The longitudinal and
survival data beyond $\tau_{max}$ are artificially censored as a result.

## Acknowledgements

This research is supported in part by NIH grants P30CA016672 and R01DK118079.

## Supporting Information

Web Appendices, Tables, and Figures are available with this paper.

## References

  * Andrinopoulou et al. (2021) Andrinopoulou, E.-R., Harhay, M. O., Ratcliffe, S. J., and Rizopoulos, D. (2021).  Reflection on modern methods: Dynamic prediction using joint models of longitudinal and time-to-event data.  International Journal of Epidemiology 50, 1731–1743.
  * Andrinopoulou et al. (2015) Andrinopoulou, E.-R., Rizopoulos, D., Takkenberg, J. J., and Lesaffre, E. (2015).  Combined dynamic predictions using joint models of two longitudinal outcomes and competing risk data.  Statistical Methods in Medical Research 26, 1787–1801.
  * Andrinopoulou et al. (2014) Andrinopoulou, E.-R., Rizopoulos, D., Takkenberg, J. J. M., and Lesaffre, E. (2014).  Joint modeling of two longitudinal outcomes and competing risk data.  Statistics in Medicine 33, 3167–3178.
  * Appel et al. (2023) Appel, L. J., Grams, M., Woodward, M., Harris, K., Arima, H., Chalmers, J., Yatsuya, H., Tamakoshi, K., Li, Y., Coresh, J., et al. (2023).  Estimated glomerular filtration rate, albuminuria, and adverse outcomes: An individual-participant data meta-analysis.  JAMA 330, 1266–1277.
  * Elashoff et al. (2016) Elashoff, R., li, G., and Li, N. (2016).  Joint Modeling of Longitudinal and Time-to-Event Data.  Chapman and Hall/CRC.
  * Elashoff et al. (2007a) Elashoff, R. M., Li, G., and Li, N. (2007a).  An approach to joint analysis of longitudinal measurements and competing risks failure time data.  Statistics in Medicine 26, 2813–2835.
  * Elashoff et al. (2007b) Elashoff, R. M., Li, G., and Li, N. (2007b).  A joint model for longitudinal measurements and survival data in the presence of multiple failure types.  Biometrics 64, 762–771.
  * Gong and Samaniego (1981) Gong, G. and Samaniego, F. J. (1981).  Pseudo maximum likelihood estimation: Theory and applications.  The Annals of Statistics 9, 861–869.
  * Hu et al. (2012) Hu, B., Li, L., Wang, X., and Greene, T. (2012).  Nonparametric multistate representations of survival and longitudinal data with measurement error.  Statistics in Medicine 31, 2303–2317.
  * Levey (1999) Levey, A. S. (1999).  A more accurate method to estimate glomerular filtration rate from serum creatinine: A new prediction equation.  Annals of Internal Medicine 130, 461\.
  * Li et al. (2016) Li, L., Luo, S., Hu, B., and Greene, T. (2016).  Dynamic prediction of renal failure using longitudinal biomarkers in a cohort study of chronic kidney disease.  Statistics in Biosciences 9, 357–378.
  * Li et al. (2018) Li, L., Wu, C.-H., Ning, J., Huang, X., Shih, Y.-C. T., and Shen, Y. (2018).  Semiparametric estimation of longitudinal medical cost trajectory.  Journal of the American Statistical Association 113, 582–592.
  * Li et al. (2009) Li, N., Elashoff, R. M., Li, G., and Saver, J. (2009).  Joint modeling of longitudinal ordinal data and competing risks survival times and analysis of the NINDS rt-PA stroke trial.  Statistics in Medicine 29, 546–557.
  * Li et al. (2012) Li, N., Elashoff, R. M., Li, G., and Tseng, C. (2012).  Joint analysis of bivariate longitudinal ordinal outcomes and competing risks survival times with nonparametric distributions for random effects.  Statistics in Medicine 31, 1707–1721.
  * Li et al. (2023) Li, W., Li, L., and Astor, B. C. (2023).  A comparison of two approaches to dynamic prediction: Joint modeling and landmark modeling.  Statistics in Medicine 42, 2101–2115.
  * Lim et al. (2022) Lim, D. K. E., Boyd, J. H., Thomas, E., Chakera, A., Tippaya, S., Irish, A., Manuel, J., Betts, K., and Robinson, S. (2022).  Prediction models used in the progression of chronic kidney disease: A scoping review.  PLOS ONE 17, e0271619.
  * Nicolaie et al. (2010) Nicolaie, M. A., van Houwelingen, H. C., and Putter, H. (2010).  Vertical modeling: A pattern mixture approach for competing risks modeling.  Statistics in Medicine 29, 1190–1205.
  * Rajeswaran et al. (2018) Rajeswaran, J., Blackstone, E. H., and Barnard, J. (2018).  Joint modeling of multivariate longitudinal data and competing risks using multiphase sub-models.  Statistics in Biosciences 10, 651–685.
  * Rizopoulos (2011) Rizopoulos, D. (2011).  Dynamic predictions and prospective accuracy in joint models for longitudinal and time-to-event data.  Biometrics 67, 819–829.
  * Rizopoulos (2012) Rizopoulos, D. (2012).  Joint Models for Longitudinal and Time-to-Event Data.  Chapman and Hall/CRC.
  * Shen and Li (2021) Shen, F. and Li, L. (2021).  Backward joint model and dynamic prediction of survival with multivariate longitudinal data.  Statistics in Medicine 40, 4395–4409.
  * Taylor et al. (2005) Taylor, J. M., Yu, M., and Sandler, H. M. (2005).  Individualized predictions of disease progression following radiation therapy for prostate cancer.  Journal of Clinical Oncology 23, 816–825.
  * van Houwelingen and Putter (2011) van Houwelingen, H. and Putter, H. (2011).  Dynamic Prediction in Clinical Survival Analysis.  CRC Press.
  * Wang et al. (2023) Wang, S., Ning, J., Xu, Y., Shih, T. Y., Shen, Y., and Li, L. (2023).  Statistical modeling of longitudinal medical cost trajectory with insurance claims data.  Annals of Applied Statistics page In press.
  * Wang et al. (2020) Wang, S., Shen, Y., Shih, Y.-C. T., Xu, Y., and Li, L. (2020).  Statistical modeling of longitudinal medical cost trajectory: renal cell cancer care cost analyses.  Biostatistics 24, 244–261.
  * Wang et al. (2023) Wang, Y., Nan, B., and Kalbfleisch, J. D. (2023).  Kernel estimation of bivariate time-varying coefficient model for longitudinal data with terminal event.  Journal of the American Statistical Association pages 1–10.
  * Wu and Li (2018) Wu, C. and Li, L. (2018).  Quantifying and estimating the predictive accuracy for censored time-to-event data with competing risks.  Statistics in Medicine 37, 3106–3124.
  * Wu et al. (2020) Wu, C., Li, L., and Li, R. (2020).  Dynamic prediction of competing risk events using landmark sub-distribution hazard model with multiple longitudinal biomarkers.  Statistical Methods in Medical Research 29, 3179–3191.
  * Yu et al. (2008) Yu, M., Taylor, J. M. G., and Sandler, H. M. (2008).  Individual prediction in prostate cancer studies using a joint longitudinal survival–cure model.  Journal of the American Statistical Association 103, 178–187.
  * Zheng and Heagerty (2005) Zheng, Y. and Heagerty, P. J. (2005).  Partly conditional survival models for longitudinal data.  Biometrics 61, 379–391.

Figure 1:  Bias and efficiency of the BJM estimation in simulations (left: BJM-3-EX; right: BJM-3-TP). The horizontal axis is the index of parameters, grouped by their roles in the model. For better visualization, we used different background colors for different groups of parameters. For the EX model, the parameter groups from left to the right are longitudinal sub-models for biomarkers 1-3, survival sub-model, and variance-covariance matrix. For the TP model, the parameter groups from left to the right are longitudinal sub-models for biomarkers 1-3, longitudinal LTS sub-models for biomarkers 1-3, survival sub-model, and variance-covariance matrix. The number of parameters in each parameter group can be read from the horizontal axis. The 4 curves demonstrate the percent bias (drawn to the vertical axis on the left) and relative efficiency (drawn to the vertical axis on the right) of each estimated parameter from the EM algorithm, under the sample size of $n=300$ or $n=1,000$. The relative efficiency is measured by the ratio of the empirical standard deviations of the point estimators from the EM and CCA. Values less than 1 indicate that EM is more efficient than CCA; smaller values indicate more efficiency improvement. Figure 2: The Kaplan-Meier survival curve of the composite event of ESRD and death in the AASK study on the left and the cumulative incidence functions of ESRD (red) and death (black) on the right.  Figure 3:  The conditional mean trajectory (CMT) of eGFR at 4 selected composite event times in the AASK study. In each plot, we aggregated the longitudinal biomarker data from all subjects who had ESRD (red) or death (black) in that year. The dots show the mean biomarker by month, averaged across all subjects who had a measurement at that month. When the measurement times are non-informative, these biomarker means are unbiased estimators of the CMT at each month. A linear model is fitted to the monthly mean biomarker data for those with ESRD and death separately, with pointwise confidence intervals shown as gray bands.  Figure 4:  Illustrating the dynamic prediction of ESRD and death with two AASK patients (patient A on the left and patient B on the right). The horizontal axis is the year since baseline. The longitudinal measurements of eGFR (green dots; connected by segments to form a continuous trajectory curve) are drawn to the vertical axis on the left; the longitudinal measurements of UPCr (blue dots and curve) and albumin (orange dots and curve) are drawn to the vertical axis on the right. The six horizontal dotted lines in the background mark the probability of 0, 0.2, 0.4, 0.6, 0.8, and 1.0, from bottom to the top. The predicted probabilities of ESRD (red dots; connected by linear segments to form a continuous curve) and death (black dots and curve) within a horizon of 3 years can be read from their positions among these six lines. The vertical red line marks the ESRD time of patient A; the vertical black line marks the death time of patient B. This dynamic prediction plot visualizes how the personalized, real-time risks of ESRD and death vary with the three longitudinal biomarkers. The prediction also incorporated baseline predictors. Figure 5: Illustrating the joint prediction of the future ESRD, death and eGFR distributions in two graphical presentations. The top plots show the prediction of AASK patient A at two prediction times (vertical brown line), based on the accumulating biomarker history to the left. The bottom plots show patient B. These are the same two patients in Figure 4 and the plotting symbols are similar. The horizon is 3 years. With patient A, we plotted the predicted probabilities at 8 equally spaced future time points ($t=0,0.5,1.0,...,3.0$). At each time, we calculated the density function of eGFR and its 9 quantiles (10%, …, 90%). These quantiles are plotted to the vertical axis on the left. For better visualization, we colored the areas between adjacent quantile curves into 8 green shaded areas. The thick green curve in the middle indicates the predicted mode of eGFR. With patient B, we plotted the multistate representation of all predicted outcomes (Hu et al., 2012). The probabilities of death, ESRD and CKD stages (defined by eGFR; see Section 1) were calculated up to the horizon. These probabilities were drawn to the dotted horizontal lines in the background, and they sum up to 1 at each time point. Table 1: Prediction accuracy results of BJM and SJM in the AASK data analysis. The prediction accuracy was evaluated among at-risk subjects at Year 1-4 after baseline. The prediction horizon is 3 years. SJM: shared random effects joint model.  Accuracy | Type | Model | Year1 | Year2 | Year3 | Year4
---|---|---|---|---|---|---
AUC | ESRD | BJM-3-TP | 0.9068 | 0.9378 | 0.9267 | 0.9008
AUC | ESRD | BJM-2-TP | 0.9068 | 0.9358 | 0.9259 | 0.8987
AUC | ESRD | BJM-1-TP | 0.8727 | 0.9135 | 0.9322 | 0.9135
AUC | ESRD | BJM-3-EX | 0.9057 | 0.9398 | 0.9342 | 0.9182
AUC | ESRD | BJM-2-EX | 0.9096 | 0.9492 | 0.9527 | 0.9587
AUC | ESRD | BJM-1-EX | 0.8796 | 0.9314 | 0.9358 | 0.9277
AUC | ESRD | SJM-3 | 0.8941 | 0.8945 | 0.9506 | 0.9449
AUC | Death | BJM-3-TP | 0.7354 | 0.6430 | 0.6736 | 0.6725
AUC | Death | BJM-2-TP | 0.7083 | 0.6099 | 0.6569 | 0.6620
AUC | Death | BJM-1-TP | 0.6586 | 0.6387 | 0.6162 | 0.5888
AUC | Death | BJM-3-EX | 0.5999 | 0.6554 | 0.6403 | 0.6295
AUC | Death | BJM-2-EX | 0.5873 | 0.6114 | 0.6098 | 0.6163
AUC | Death | BJM-1-EX | 0.6242 | 0.6497 | 0.6605 | 0.6566
AUC | Death | SJM-3 | 0.6709 | 0.6316 | 0.6243 | 0.6874
BS | ESRD | BJM-3-TP | 0.0672 | 0.0585 | 0.0677 | 0.0714
BS | ESRD | BJM-2-TP | 0.0659 | 0.0586 | 0.0677 | 0.0710
BS | ESRD | BJM-1-TP | 0.0776 | 0.0690 | 0.0652 | 0.0672
BS | ESRD | BJM-3-EX | 0.0686 | 0.0582 | 0.0673 | 0.0693
BS | ESRD | BJM-2-EX | 0.0677 | 0.0640 | 0.0739 | 0.0809
BS | ESRD | BJM-1-EX | 0.0809 | 0.0786 | 0.0812 | 0.0921
BS | ESRD | SJM-3 | 0.1177 | 0.1151 | 0.1078 | 0.0718
BS | Death | BJM-3-TP | 0.0422 | 0.0512 | 0.0642 | 0.0748
BS | Death | BJM-2-TP | 0.0425 | 0.0514 | 0.0646 | 0.0741
BS | Death | BJM-1-TP | 0.0422 | 0.0486 | 0.0588 | 0.0674
BS | Death | BJM-3-EX | 0.0442 | 0.0504 | 0.0614 | 0.0710
BS | Death | BJM-2-EX | 0.0550 | 0.0699 | 0.0878 | 0.1020
BS | Death | BJM-1-EX | 0.0547 | 0.0694 | 0.0868 | 0.1009
BS | Death | SJM-3 | 0.0867 | 0.1081 | 0.1279 | 0.1306
RMSE | eGFR | BJM-3-TP | 8.1559 | 8.7276 | 8.4907 | 8.2903
RMSE | eGFR | BJM-3-EX | 10.6981 | 10.6853 | 9.5622 | 8.4871
RMSE | eGFR | SJM-3 | 8.8914 | 8.9091 | 8.0943 | 6.8213
P30 | eGFR | BJM-3-TP | 0.9156 | 0.9241 | 0.9469 | 0.9308
P30 | eGFR | BJM-3-EX | 0.7891 | 0.8254 | 0.8607 | 0.9000
P30 | eGFR | SJM-3 | 0.9213 | 0.9195 | 0.9110 | 0.9466
P50 | eGFR | BJM-3-TP | 0.9600 | 0.9709 | 0.9812 | 0.9853
P50 | eGFR | BJM-3-EX | 0.9553 | 0.9520 | 0.9614 | 0.9742
P50 | eGFR | SJM-3 | 0.9683 | 0.9831 | 0.9874 | 0.9871
# at-risk patients in AASK data | 1065 | 1005 | 955 | 875


[a]Rajnandini Mukherjee

# The $I=1/2$ and $3/2$ K-$\pi$ scattering length with domain wall fermions at
physical pion mass with all-to-all propagators

Nils Asmussen  Felix Erben  Jonathan Flynn  Andreas Jüttner  Christopher T.
Sachrajda

###### Abstract

We present our calculations for the $I=1/2$ and $3/2$ K-$\pi$ s-wave
scattering length with physical quark masses, extracted from the interaction
energy of Euclidean two-point functions. We use the domain wall fermion action
with physical quark masses at a single lattice spacing. We are specifically
interested in the systematic effects due to around-the-world terms on the
overall determination of the scattering length. We present our progress and
discuss the various systematic effects in our preliminary results.

##  1 Introduction

To further our understanding of strong interactions in the standard model, it
is of interest to look at particle resonances at low energy, such as kaon-pion
scattering. While this can be formulated using an effective-theory approach
such as chiral perturbation theory [1], simulations of lattice QCD can provide
both qualitative and quantitative insight from first principles. This is one
of the cleanest scattering processes with regard to systematics and is a
stepping stone towards high precision calculations that are particulary
sensitive to certain systematic effects, such as in [2] and [3].

In this project, we compute the s-wave $K\pi$ scattering length in the two
isospin channels: $I=1/2$ and $3/2$. There are previous determinations for
these two quantities [4, 5, 6] at unphysically heavy pion and/or kaon masses
which use chiral perturbation theory to extrapolate the results to the
physical pion and kaon mass. In our work, we calculate the scattering length
directly at physical masses, using 23 configurations on a $48^{3}\times 96$
lattice with inverse spacing $a^{-1}=1.730(4)$ GeV [7]. We use the Möbius
domain wall fermion action and an all-to-all setup with meson-field
contractions for computing propagators. We compute a number of low modes of
the Dirac operator exactly, and the high modes are stochastically sampled by
inverting the Dirac operator on $Z_{2}\times Z_{2}$ noise sources.

##  2 Methodology

We are interested in the s-wave scattering length, $a_{0}$, described by
Lüscher’s formula [8] for interactions in a Euclidean volume of spatial extent
$L$ as

$\displaystyle\Delta E_{K\pi}$ $\displaystyle=-\frac{2\pi
a_{0}}{\mu_{K\pi}L^{3}}\left(1+c_{1}\frac{a_{0}}{L}+c_{2}\frac{a_{0}^{2}}{L^{2}}\right)+\mathcal{O}\Big{(}\frac{1}{L^{6}}\Big{)},$
(1)

where $c_{1}=-2.837297$, $c_{2}=6.375183$ and
$\mu_{K\pi}=m_{\pi}m_{K}/(m_{\pi}+m_{K})$ is the reduced mass. Eqn. (1)
relates the scattering length to $\Delta E_{K\pi}$, the binding energy of the
$K\pi$ composite state. This in turn can be extracted by studying the ground
state of the 2-point function of some $K\pi$ operator. In our study, we choose
the following $K\pi$ operator

$\displaystyle O_{K\pi}(t)=K(t+\delta)\pi(t),$ (2)

where the kaon and pion are separated in time by $\delta=a$ in order to avoid
unwanted contributions from additional Wick contractions of the correlated
$Z_{2}$ noise sources. The corresponding 2-point function

$\displaystyle C_{K\pi}(t)=$ $\displaystyle\langle
O^{\text{snk}}_{K\pi}(t)O^{\dagger\text{src}}_{K\pi}(0)\rangle=\langle
K(t+\delta)\pi(t)(K(\delta)\pi(0))^{\dagger}\rangle,$ (3)

has a straightforward spectral decomposition in infinite volume with ground
state energy $E_{K\pi}=m_{K}+m_{\pi}+\Delta E_{K\pi}$. However, periodic
boundary conditions introduce ‘around-the-world’ effects, giving rise to
additional terms

$\displaystyle C_{K\pi}(t)=$ $\displaystyle|\langle
K\pi|K(\delta)\pi(0)|0\rangle|^{2}\left(e^{-E_{K\pi}t}+e^{-E_{K\pi}(T-t)}\right)$
(4)
$\displaystyle+|\langle\pi|K(\delta)\pi(0)|K\rangle|^{2}e^{-m_{K}t}e^{-m_{\pi}(T-t)}$
(5) $\displaystyle+|\langle
K|K(\delta)\pi(0)|\pi\rangle|^{2}e^{-m_{\pi}t}e^{-m_{K}(T-t)}$ (6)
$\displaystyle+\cdots$

where $T$ is the time extent of the lattice, and further terms are suppressed
by multiples of pion and kaon masses. The 2-point function has a cosh-like
term (4) and two around-the-world (ATW) terms (5) and (6) corresponding to
single meson states going around the periodic boundary. The binding energy,
$\Delta E_{K\pi}$, appears in the cosh-like term and to extract it precisely
from $C_{K\pi}$, we must account for the two ATW terms. Note that this is the
case even with the modest statistical precision of our study.

Using all-to-all propagators should allow for low statistical error per
configuration in $C_{K\pi}(t)$ by averaging multiple source planes. However,
it is possible to have further noise cancellations on each configuration by
instead studying the ratio

$\displaystyle
R_{K\pi}(t)=\frac{C_{K\pi}(t)}{C_{K}(t)C_{\pi}(t)}=\frac{\langle
K(t+\delta)\pi(t)(K(\delta)\pi(0))^{\dagger}\rangle}{\langle\pi(t)\pi(0)^{\dagger}\rangle\langle
K(t+\delta)K(\delta)^{\dagger}\rangle},$ (7)

which also separates into a cosh-like piece, $R_{K\pi}^{\text{cosh}}(t)$, and
an ATW piece $R_{K\pi}^{\text{ATW}}(t)$. Here we use the pion and kaon 2-point
functions, which we separately fit to also extract the values of the pion and
kaon masses and their amplitudes.

The ATW part of this ratio is given by

$\displaystyle R^{\text{ATW}}_{K\pi}(t)=$
$\displaystyle\left|\frac{\langle\pi|K(\delta)\pi(0)|K\rangle}{\langle
0|K|K\rangle\langle\pi|\pi|0\rangle}\right|^{2}\frac{e^{-m_{K}t}e^{-m_{\pi}(T-t)}}{\left(e^{-m_{\pi}t}+e^{-m_{\pi}(T-t)}\right)\left(e^{-m_{K}t}+e^{-m_{K}(T-t)}\right)}$
(8) $\displaystyle+\left|\frac{\langle K|K(\delta)\pi(0)|\pi\rangle}{\langle
0|K|K\rangle\langle\pi|\pi|0\rangle}\right|^{2}\frac{e^{-m_{\pi}t}e^{-m_{K}(T-t)}}{\left(e^{-m_{\pi}t}+e^{-m_{\pi}(T-t)}\right)\left(e^{-m_{K}t}+e^{-m_{K}(T-t)}\right)}.$
(9)

All terms in this equation are fully determined by the information from the
pion and kaon 2-point functions except for the matrix elements
$\langle\pi|K(\delta)\pi(0)|K\rangle$ and $\langle
K|K(\delta)\pi(0)|\pi\rangle$. In order to compute them, we separately
construct correlation functions — called $C_{\pi K\pi K}$ and $C_{KK\pi\pi}$
respectively — that carry these matrix elements in their ground state
contributions. In particular, we use the following 3-point functions

$\displaystyle C_{\pi K\pi K}$
$\displaystyle=\langle\pi(\Delta)K(t+\delta)\pi(t)K(0)\rangle$ (10)
$\displaystyle=\langle
0|\pi|\pi\rangle\langle\pi|K(\delta)\pi(0)|K\rangle\langle K|K|0\rangle
e^{-m_{\pi}(\Delta-t)}e^{-m_{K}t}+\cdots,$ $\displaystyle C_{KK\pi\pi}$
$\displaystyle=\langle K(\Delta)K(t+\delta)\pi(t)\pi(0)\rangle$ (11)
$\displaystyle=\langle 0|K|K\rangle\langle
K|K(\delta)\pi(0)|\pi\rangle\langle\pi|\pi|0\rangle
e^{-m_{\pi}t}e^{-m_{K}(\Delta-t)}+\cdots.$

The unknown matrix elements $\langle\pi|K(\delta)\pi(0)|K\rangle$ and $\langle
K|K(\delta)\pi(0)|\pi\rangle$ can then be extracted by fitting these 3-point
functions to their analytical forms. In this way we can now completely account
for the ATW contributions and what remains is to extract the signal. As
before, the signal lies in the cosh-like term

$\displaystyle
R_{K\pi}^{\text{cosh}}(t)=A\frac{\left(e^{-({m_{\pi}+m_{K}}+\Delta
E_{K\pi})t}+e^{-(m_{\pi}+m_{K}+\Delta
E_{K\pi})(T-t)}\right)}{\left(e^{-m_{\pi}t}+e^{-m_{\pi}(T-t)}\right)\left(e^{-m_{K}t}+e^{-m_{K}(T-t)}\right)},$
(12)

where all the terms are known except the amplitude, $A$, and $\Delta E_{K\pi}$
which we fit for. Using the fit value of $\Delta E_{K\pi}$ in Lüscher’s
formula (1), we are able to determine the value of the scattering length. The
lattice quantities required for this analysis are the pion and kaon 2-point
functions, the $K\pi$ 2-point function, and the two 3-point functions for
determining the unknown ATW matrix elements.

Now that the procedure for calculating the scattering length is well defined,
recall that we must do this for the two isospin channels $I=1/2$ and $3/2$.
The above described procedure remains the same for the two isospin channels,
however, in each case a given correlator is made up of weighted diagrammatic
contributions corresponding to the isospin-dependent Wick contractions for
that channel. The source operator in eqn. (3) is fixed to

$O_{K\pi}^{\text{src}}(t)=T\left[(\overline{u}\gamma_{5}s)(t+\delta)(\overline{u}\gamma_{5}d)(t)\right],$
(13)

while different sink operators are used to project onto the isospin $1/2$ or
$3/2$ channels

$\displaystyle O_{K\pi}^{\text{snk}\,I=1/2}(t)$
$\displaystyle=T\left[(\overline{d}\gamma_{5}s)(t+\delta)\;(\overline{u}\gamma_{5}d)(t)-\tfrac{1}{2}(\overline{u}\gamma_{5}s)(t+\delta)\;(\overline{d}\gamma_{5}d)(t)\right.$
$\displaystyle\qquad+\left.\tfrac{1}{2}(\overline{u}\gamma_{5}s)(t+\delta)\;(\overline{u}\gamma_{5}u)(t)\right],$
(14) $\displaystyle O_{K\pi}^{\text{snk}\,I=3/2}(t)$
$\displaystyle=T\left[(\overline{u}\gamma_{5}s)(t+\delta)\;(\overline{u}\gamma_{5}d)(t)\right].$
(15)

The $C_{K\pi}$ correlation function is thus constructed from the combinations

$\displaystyle C_{K\pi}^{I=1/2}=D+\frac{1}{2}C-\frac{3}{2}R,\quad
C_{K\pi}^{I=3/2}=D-C,$ (16)

where $D$, $C$, and $R$ are the direct, crossed and rectangular contractions,
as illustrated in figure 1. The same distinction is made when constructing the
3-point functions (10) and (11) for the two isospin channels.

(a) D (b) C (c) R

Figure 1: Diagrams corresponding to various Wick contractions that contribute
to $C_{K\pi}$. Operators with subscripts $1$ and $2$ correspond to the source
and sink operators respectively.

##  3 Results

Previous results for the $K\pi$ scattering length in the two isospin channels
are listed in table 1. We find sensible fits to $R_{K\pi}(t)$ in both
channels, as shown in figure 2. However, we find a systematic drift in the fit
results upon varying the $t_{min}$ and $t_{max}$ of the fit interval, shown in
figure 3. Thus the fit result is not stable. To explain this trend in the fit
value of $\Delta E_{K\pi}$ we consider the following possible explanations:
effect of ATW contributions, contamination by excited states and availability
of statistics.

|  $m_{\pi}/$MeV | $m_{\pi}a_{0}^{I=1/2}$ | $m_{\pi}a_{0}^{I=3/2}$
---|---|---|---
Sasaki (2014) [4] | extrap | 0.142(14)(27) | -0.0469(24)(20)
166 | 0.158(36) | -0.108(12)
ETMC (2018) [5] | extrapI=I/2/$\chi_{PT}^{I=3/2}$ | 0.163(3) | -0.059(2)
hadspec (2019) [6] | 239 | 0.46(3) |
284 | 0.79(13) |
Table 1: Previous computations of $K\pi$ scattering length. The first error is
statistical and the second is systematic error due to chiral extrapolation.

(a)  (b)

Figure 2: Fits to $R_{K\pi}(t)$ data in the two isospin channels while
accounting for ATW contributions (the orange points correspond to the fit
interval). Fit values of $\Delta E_{K\pi}$ correspond to scattering lengths
$m_{\pi}a_{0}^{I=1/2}=0.160(3)$ and $m_{\pi}a_{0}^{I=3/2}=-0.0377(2)$ using
Lüscher’s formula in eqn. (1).

(a)  (b)

Figure 3: Variation in fit value of $\Delta E_{K\pi}$ with fit intervals,
extracted from $R_{K\pi}(t)$ for each isospin channel. This shows a strong
systematic dependence of the results on the choice of both $t_{min}$
($x$-axis) and $t_{max}$ (color in legend).

###  3.1 ‘Around-the-world’ effects

(a)  (b)

Figure 4: Effect of the ATW terms in fitting $R_{K\pi}(t)$. These plots show
the time-dependence of the ratio of the cosh-like component of fits using two
different ansatz: $R_{K\pi}(t)=R_{K\pi}^{cosh}(t)$ and
$R_{K\pi}(t)=R_{K\pi}^{cosh}(t)+R_{K\pi}^{ATW}(t)$.

(a)  (b)

Figure 5: Variation in fit value of $\Delta E_{K\pi}$ with fit intervals,
extracted from $R_{K\pi}(t)$ from two ansatz, one accounting for ATW
contributions, and the other not. It can be seen here that the ATW terms have
a constant additive contribution to the value of $\Delta E_{K\pi}$.

By looking at the functional form of $R_{K\pi}^{\text{ATW}}(t)$, eqns. (10)
and (11), it is evident that ATW terms have time-dependent contributions. As
the systematic trend we see in the fits to $R_{K\pi}(t)$ is also time-
dependent, it is worthwhile to study the effect of these ATW terms in the
overall determination of $\Delta E_{K\pi}$. We investigate this effect by
studying the cosh-like component of fits done with and without accounting for
the ATW terms. This is to quantify how the ATW terms are absorbed into the
cosh only ansatz when the ATW terms are unaccounted for, as shown in figure 4.
This tells us that the ATW terms have a time-dependent linear contribution at
all times. To understand their effect on the value of $\Delta E_{K\pi}$, we
look at eqns. (10-12) at small $t$, where

$\displaystyle R_{K\pi}(t)\approx A+B\Delta
E_{K\pi}t+C_{ATW}t+\mathcal{O}(t^{2})$ (17)

which implies that the ATW terms shift the value of $\Delta E_{K\pi}$ when
they are not accounted for. This shift is shown in figure 5 for both isospin
channels. The effect of the ATW terms is particularly important in the $I=3/2$
channel where the shift is comparable in magnitude to $\Delta E_{K\pi}$.

ATW terms are substantial and present at all times, which is a particular
feature of our study owing to the (low) physical pion mass which allows the
ATW terms to be more dominant. It is therefore important to account for these
terms in any further studies made at the physical point for a reliable
determination of $\Delta E_{K\pi}$. However, as shown in figure 5, they lead
to a constant shift and do not account for the time-dependence.

###  3.2 Excited states

To investigate further the source of the time-dependent trend in the fit
values of $\Delta E_{K\pi}$, we next study the effect of excited states in our
$R_{K\pi}$ data. We do this by looking at the effective mass plot for
$R_{K\pi}$ using our data, generated using all-to-all propagators, and data
generated using a distillation procedure [2] where the ground state of a
$K\pi$ operator is determined to a greater accuracy. This is shown in figure
6, which reveals excited state contamination at low times. However these
effects die out at larger times, and are negligible in the range where we
perform fits (the grey area in the plots). While the distillation data has
greater precision in comparison to the all-to-all data, this comparison shows
us that the systematic trend we see in the fit value of $\Delta E_{K\pi}$ is
unlikely a consequence of excited state effects.

(a)  (b)

Figure 6: Effective mass plots for $R_{K\pi}(t)$ in the two isospin channels,
comparing the contamination of excited states in the all-to-all data compared
to data generated via distillation. The grey region corresponds to the region
where fits are made. $N_{c}$ refers to the number of configurations used.

###  3.3 Statistics

$\Delta E_{K\pi}$ is an extremely delicate quantity in comparison to the other
observables we compute in this analysis. For example, the pion and kaon masses
on the lattice are $am_{\pi}=0.0803(2)$ and $am_{K}=0.2887(2)$, and therefore
$a\Delta E_{K\pi}$ is of the order of magnitude of the error in these masses.
In our study we attempt to extract this using 23 configurations. While the
all-to-all setup allows us to compute propagators with relatively low
statistical error due to availability of a large number of source planes per
configuration, the overall number of configurations is a limiting factor.
Figure 6 also indicates this effect of statistics as the distillation study
has more configurations ($N_{c}=35$). However, comparison of the effective
masses suggests it to be unlikely that the time-dependent trend in the fit for
$\Delta E_{K\pi}$ is a direct effect of statistics, especially in the $I=3/2$
channel where the all-to-all data is in good agreement with the distillation
data. We are therefore in need of larger statistics to pin down the
systematics effects we see.

##  4 Conclusion

In this study we have attempted to calculate the $K\pi$ scattering length at
physical pion mass. As this is one of the first calculations done directly at
the physical point for both isospin channels, we encounter new features and
new issues. The most interesting feature of our study is the effect of around-
the-world contributions, which are substantial due to the light pion. This
effect is significant at all times on the lattice especially when extracting a
small and sensitive quantity such as the binding energy, $\Delta E_{K\pi}$.
The time-dependent trend in the value of this quantity remains unexplained and
may be a combination of effects considered above or an indicator of physics
that has not been taken into account. We look forward with interest towards
other analyses carried out at the physical point, where we expect around-the-
world effects to play an important role in the study of scattering lengths.

##  5 Acknowledgements

The distillation data is shared by Nelson Lachini, Antonin Portelli, and F.E.
This work used the DiRAC Extreme Scaling service at the University of
Edinburgh, operated by the Edinburgh Parallel Computing Centre on behalf of
the STFC DiRAC HPC Facility (www.dirac.ac.uk) under project codes dp008 and
dp207 using Grid and Hadrons. This equipment was funded by BEIS capital
funding via STFC capital grant ST/R00238X/1 and STFC DiRAC Operations grant
ST/R001006/1. DiRAC is part of the National e-Infrastructure. R.M. is
supported by the Presidential Scholarship from the University of Southampton.
C.T.S. is partially supported by an Emeritus Fellowship from the Leverhulme
Trust and by STFC (UK) grants ST/P000711/1 and ST/T000775/1. A.J. and J.M.F.
also acknowledge funding by ST/P000711/1 2017-2021 and ST/T000775/1. F.E. is
supported in part by UK STFC grant ST/P000630/1. F.E. also received funding
from the European Research Council (ERC) under the European Union’s Horizon
2020131 research and innovation programme under grant agreement No 757646.

## References

  * [1] M. Jamin et al., _S-wave Kpi scattering in chiral perturbation theory with resonances,Nucl. Phys. B 587, 331-362 (2000)_
  * [2] N. Lachini et al., _$K\pi$ scattering at physical pion mass using distillation, PoS LATTICE2021 435 (2022)_
  * [3] F. Erben et al., _Exploring distillation at the SU(3) flavour symmetric point,PoS LATTICE2022_
  * [4] K. Sasaki et al., _Scattering lengths for two pseudoscalar meson systems,Phys. Rev. D 89, 054502 (2014)_
  * [5] C. Helmes et al., _Hadron-Hadron Interactions from $N_{f}=2+1+1$ Lattice QCD: $I=3/2$ $\pi K$ Scattering Length, Phys. Rev. D 98, 114511 (2018)_
  * [6] D. Wilson et al., _The quark-mass dependence of elastic $\pi K$ scattering from QCD, Phys. Rev. Lett. 123, 042002 (2019)_
  * [7] T. Blum et al., _Domain wall QCD with physical quark masses,Phys. Rev. D 93, 074505 (2016)_
  * [8] M Lüscher, _Volume dependence of the energy spectrum in massive quantum field theories,Commun. Math. Phys. 104 (1986) 177._


# On a question of Matt Baker regarding the dollar game

Marine Cases

(April 2023)

###### Abstract

In [1], Matt Baker wrote the following about the dollar game played on a graph
: “The total number of borrowing moves required to win the game when playing
the “borrowing binge strategy” is independent of which borrowing moves you do
in which order! Note, however, that it is usually possible to win in fewer
moves by employing lending moves in combination with borrowing moves. The
optimal strategy when one uses both kinds of moves is not yet understood.”

In this article we give a lower bound on the minimum number $M_{\min}$ of such
moves of an optimal algorithm in term of the number of moves $M_{0}$ of the
borrowing binge strategy. Concretely, we have :
$M_{\min}\geq\dfrac{M_{0}}{n-1}$ where $n$ is the number of vertices of the
graph. This bound is tight.

###### Contents

  1. 1 Introduction
    1. 1.1 Definitions
  2. 2 Main Theorem
    1. 2.1 Examples where the bound in Theorem 1 is optimal
  3. 3 Acknowledgements

##  1 Introduction

In this note, we consider the dollar game on a connected graph $G$ with
initial divisor $C$, and focus on the minimum number of borrowing and lending
moves needed to win the game (i.e, reach an all non-negative divisor). If the
game is winnable, then it can always be won using only borrowing moves. More
precisely, any winnable game can be won by simply using the borrowing binge
strategy : any time there are vertices in debt, pick one of them and do a
borrowing move. Repeat until everyone is out of debt!
This strategy is optimal if you are constrained to only use borrowing moves,
but not if you can do both lending and borrowing moves, as seen below:

-1020233-11-113borrow on red2-10212borrow on red12-1211borrow on red012111borrow on red002004lend on blueEffective divisor reached in $M_{\min}={\color[rgb]{0.1953125,0.1953125,0.78515625}\definecolor[named]{pgfstrokecolor}{rgb}{0.1953125,0.1953125,0.78515625}1}$ moveEffective divisor reached in $M_{0}={\color[rgb]{0.78515625,0.1953125,0.1953125}\definecolor[named]{pgfstrokecolor}{rgb}{0.78515625,0.1953125,0.1953125}4}$ movesBinge Borrowing algorithm

In this article we will prove a bound $M_{\min}\geq\dfrac{M_{0}}{n-1}$ on how
far from optimal the Borrowing Binge strategy can be, where n is the number of
vertices of G.

We start by reinterpreting the question in terms of the chip firing game
through a well-known bijection $v_{0}\rightarrow K-v_{0}$ (See [3], p 40),
where $K=(d_{i}-1)_{i}$ is the maximal stable divisor.

In this context, the borrowing binge strategy becomes the standard greedy
chip-firing game : any time there are vertices with at least as many chips as
their degree, pick one of them and do a lending move; repeat until the divisor
is stable and no move can be made. The effective divisor then becomes a stable
divisor. Then the question becomes what is the minimal number of lending and
borrowing moves to reach a stable divisor.

The proof of our main theorem $M_{\min}\geq\dfrac{M_{0}}{n-1}$, with
$M_{\min}$ the minimum number of moves of an optimal algorithm and $M_{0}$ the
number of moves of the borrowing binge strategy, relies on two lemmas :

  1. 1.

The principle of least action which says that the greedy chip firing game
strategy is optimal to reach a stable divisor if we restrict to only lending
moves.

  2. 2.

A lemma which tells us by performing only lending moves, we can at most cut
the length by a multiplicative factor $\dfrac{1}{n-1}$ if we allow borrowing
moves as well.

Combining the two lemmas give the theorem.

###  1.1 Definitions

First, let us introduce some basic definitions. We will consider only simple
connected graphs throughout this article.

######  Definition 1.

  * •

The Laplacian matrix of G, denoted $\Delta$, is a matrix that captures the
structure of a graph. Its elements are given by:

$\Delta_{ij}=\begin{cases}d_{i}&\text{if }i=j.\\\ -1&\text{if }i\neq j\text{
and vertices }i\text{ and }j\text{ are adjacent.}\\\
0&\text{otherwise.}\end{cases}$

where $d_{i}$ is the degree of the vertex $i$.

  * •

A divisor is a vector in $\mathbf{Z}^{n}$ that represents the number of chips
on every vertex of the graph.

  * •

Given two divisors $C$ and $C^{\prime}$, and $v\in\mathbf{Z}^{n}$, the
notation $C\xrightarrow{v}C^{\prime}$ means that $C^{\prime}=C-\Delta\cdot v$.

  * •

Let $C$ be a divisor:

    * –

A lending move on vertex $i$ corresponds to $C\xrightarrow{e_{i}}C^{\prime}$
and means that the vertex $i$ lends one chip to each of its neighbors.

    * –

A borrowing move on vertex $i$ corresponds to
$C\xrightarrow{-e_{i}}C^{\prime}$ and means that the vertex $i$ borrows one
chip from each of its neighbors.

  * •

An effective divisor $C$ is a divisor with $C\geq 0$.

  * •

A stable divisor $C$ is a divisor with $C\leq K=(d_{i}-1)_{i}$, meaning that
one cannot perform a lending move without becoming indebted, i.e
$(C^{\prime})_{i}<0$ for some i .

  * •

A closest stable divisor $C^{*}$ to a given divisor $C$ is defined as a stable
divisor such that the number of moves (either borrowing or lending) required
to reach any other stable divisor from $C$ is no less than the number of moves
required to reach $C^{*}$ from $C$.

  * •

The distance between two divisors is the minimal number of moves required to
reach one from the other.

Let $C$ be a divisor and consider $C^{\prime}=K-C$. Then a borrowing move if a
vertex $i$ has a non positive number of chips in $C$ corresponds to a lending
move if a vertex $i$ has at least $d_{i}$ chips.

We have then the following equivalences:

Binge strategy for dollar game |  Chip firing game
---|---
Effective divisor |  Stable divisor
Borrow on negative vertices |  Lend on vertices with at least as many chips as the degree
Borrowing binge strategy to reach an effective divisor |  Standard greedy chip firing game (without a sink) to reach a stable divisor
Riemann-Roch formula (see [2]): number of chips that can be removed before no longer having an effective divisor in the class (or, equivalently, the borrowing binge is no longer terminating) |  Riemann-Roch formula: number of chips that can be added before no longer having stable divisor in the class (or, equivalently, the game is no longer terminating)

##  2 Main Theorem

######  Theorem 1.

Let $C\xrightarrow{v_{0}}X_{0}$, where $v_{0}$ is the firing sequence found
with the greedy algorithm. If $\left|v_{0}\right|_{1}=M_{0}$ then the closest
stable divisor must be at distance at least $\dfrac{M_{0}}{n-1}$.

The two examples afterwards prove that the bound in Theorem 1 is tight.

The following is a classic lemma from chip-firing game theory :

######  Lemma 1.

(Least action principle.) Let $C\xrightarrow{v_{0}}X_{0}$ be the sequence from
the greedy algorithm of the classical chip-firing game. Then, if
$C\xrightarrow{v_{1}}X_{1}$ is another sequence to another stable
configuration, we have that if $v_{1}\geq 0$, then $v_{1}\geq v_{0}$. (For a
proof of this lemma, see [3], pp.39-40)

######  Lemma 2.

Let $v\geq 0$ be a vector with at least one zero coordinate such that
$C\xrightarrow{v}X$. Then the shortest path $v^{\prime}$ from $C$ to $X$
satisfies $\left|v^{\prime}\right|_{1}\geq\dfrac{\left|v\right|_{1}}{n-1}$.

###### Proof.

In this proof, $k/2$ mean the integer $\lfloor k/2\rfloor$, for every integer
$k$.

Since $\mathbf{1}=\begin{pmatrix}1\\\ \vdots\\\ 1\end{pmatrix}$ is a basis for
$\text{ker }\Delta$, we are looking for
$v_{min}=v+k\cdot\mathbf{1},k\in\mathbf{Z}$ with minimal $L^{1}$-norm.

First, note that if $v^{\prime}$ has $m$ positive coordinates, $l$ zero
coordinates, and $r=n-m-l$ negative coordinates, then

$\displaystyle\left|v^{\prime}-1\cdot\mathbf{1}\right|_{1}$
$\displaystyle=\left|v^{\prime}\right|_{1}+r+l-m$
$\displaystyle=\left|v^{\prime}\right|_{1}+n-2m$
$\displaystyle<\left|v^{\prime}\right|_{1}$ $\displaystyle\text{if
}m>\dfrac{n}{2}$

and

$\displaystyle\left|v^{\prime}+1\cdot\mathbf{1}\right|_{1}$
$\displaystyle=\left|v^{\prime}\right|_{1}-r+l+m$
$\displaystyle=\left|v^{\prime}\right|_{1}+n-2r$
$\displaystyle<\left|v^{\prime}\right|_{1}$ $\displaystyle\text{if
}r>\dfrac{n}{2}.$

Thus the minimum can be obtained only when the number of positive and negative
components are both at most $n/2$. Therefore, we must remove
$k\cdot\mathbf{1}$ from $v$, with $k$ chosen so that the new vector
$v_{min}=v-k\cdot\mathbf{1}$ has exactly $n/2$ positive or zero components and
$n/2$ negative or zero components. In this case, the $L^{1}$ norm of this
vector will be minimal. Concretely :

  * •

For $n$ even, the vector with $\frac{n}{2}-1$ positive components,
$\frac{n}{2}$ negative components and 1 zero component is a minimum.

  * •

For $n$ odd, the minimum is unique. The vector with $\frac{n}{2}$ positive and
negative components, and 1 zero component is the minimun.

Both cases can thus be dealt with simultaneously by writing: $\frac{n-1}{2}$
positive components, $\frac{n}{2}$ negative components and 1 zero component.

Let us now compute this minimal norm. We want to show that
$\left|v_{\min}\right|_{1}\geq\dfrac{\left|v\right|_{1}}{n-1}$.

Without loss of generality, let us assume that $v$ satifies

$v=\begin{pmatrix}x_{1}\\\ \vdots\\\ x_{n}\end{pmatrix}\quad
0=x_{1}\leq\ldots\leq x_{n}.$

(We sort its components for simplicity.)

Then $k=x_{\frac{n-1}{2}}$ and

$v_{min}=\begin{pmatrix}0\\\ x_{2}\\\ \vdots\\\
x_{n}\end{pmatrix}-x_{\frac{n-1}{2}}\cdot\begin{pmatrix}1\\\ \vdots\\\
1\end{pmatrix}=\begin{pmatrix}-x_{\frac{n-1}{2}}\\\ x_{2}-x_{\frac{n-1}{2}}\\\
\vdots\\\ x_{n}-x_{\frac{n-1}{2}}\end{pmatrix}.$

Thus

$\displaystyle\left|\begin{pmatrix}-x_{\frac{n-1}{2}}\\\
x_{2}-x_{\frac{n-1}{2}}\\\ \vdots\\\
x_{n}-x_{\frac{n-1}{2}}\end{pmatrix}\right|_{1}$
$\displaystyle=x_{\frac{n-1}{2}}+\sum^{\frac{n-1}{2}-1}_{i=2}x_{\frac{n-1}{2}}-x_{i}+\sum^{n}_{i=\frac{n-1}{2}+1}x_{i}-x_{\frac{n}{2}}$
$\displaystyle=x_{\frac{n-1}{2}}+\underbrace{\sum^{\frac{n-1}{2}-1}_{i=2}x_{\frac{n-1}{2}}-x_{i}}_{\geq
0}+\underbrace{\sum^{n-1}_{i=\frac{n-1}{2}+1}x_{i}-x_{\frac{n-1}{2}}}_{\geq
0}+(x_{n-1}-x_{\frac{n-1}{2}})$ $\displaystyle\geq x_{n}$
$\displaystyle\geq\dfrac{\sum^{n}_{i=2}x_{i}}{n-1}=\dfrac{\sum^{n-1}_{i=1}x_{i}}{n-1}\text{
since $x_{1}=0$}$ $\displaystyle=\dfrac{\left|v\right|_{1}}{n-1}.$

∎

Using these two lemmas, we can show our theorem:

###### Proof.

Let $C\xrightarrow{v_{0}}X_{0}$ be the sequence from the greedy algorithm. We
want to find a closest stable divisor. Let $v_{1}$ be a sequence from $C$ to a
closest stable divisor. Since $\begin{pmatrix}1\\\ \vdots\\\ 1\end{pmatrix}$
is a basis for $\text{ker }\Delta$, we can always choose $v_{1}\geq 0$ with
one zero coordinate. By Lemma 2, we have $v_{1}\geq v_{0}\geq 0$.

Now, by Lemma 1, the shortest path from C to the stable divisor that $v_{1}$
reaches is $v_{1}^{{}^{\prime}}$, with
$\left|v_{1}^{\prime}\right|_{1}\geq\dfrac{\left|v_{1}\right|_{1}}{n-1}$.

Moreover, since $v_{1}\geq v_{0}\geq 0$ , we have
$M_{min}=\left|v_{1}^{\prime}\right|_{1}\geq\dfrac{\left|v_{1}\right|_{1}}{n-1}\geq\dfrac{\left|v_{0}\right|_{1}}{n-1}=\dfrac{M_{0}}{n-1}$.
∎

###  2.1 Examples where the bound in Theorem 1 is optimal

In this section, we will construct two examples to show that the bound we
obtained is tight.
Example 1 :

$-nk$$k$$k$$k$$k$$\cdots$$-k$$0$$0$$0$$0$$\cdots$$C$$\Delta(v_{0})$$X_{0}$

Let $C\xrightarrow{v_{0}}X_{0}$ with $v_{0}=\begin{pmatrix}k\\\ \vdots\\\ k\\\
0\\\ k\\\ \vdots\\\ k\end{pmatrix}.$ We can choose
$v_{min}=v_{0}-k\cdot\begin{pmatrix}1\\\ 1\\\ \vdots\\\ 1\\\ 1\\\
1\end{pmatrix}=\begin{pmatrix}0\\\ \vdots\\\ 0\\\ -k\\\ 0\\\ \vdots\\\ 0\\\
\end{pmatrix}.$

Note here that we didn’t even need to find another stable divisor : the same
divisor with another (shorter but equivalent) vector of firings give us the
best possible result.

Example 2:  The preceding example offers an optimal, linear (in $n$)
improvement. One might think that this is the only way to achieve the optimal
bound, through an equivalent and much shorter path to the same divisor, via
the following algorithm : compute the greedy path, identify the vector with
the smallest norm in $\mathbf{coker}(\Delta)$, and use this ‘renormalization’
to obtain a linear improvement. This suggests that, once we take into account
the ‘renormalized’ path to the greedy divisor, we might no longer obtain a
linear improvement by finding a path to a closer stable divisor. However, our
second example demonstrates that this is not the case. More precisely, we
prove that we can achieve an improvement by a factor of $n/2$ by reaching a
stable divisor entirely different from the greedy one, and this factor of
$n/2$ holds even when we account for the ‘renormalized’ path to the greedy
divisor.

$n$ will denote an even integer in the following.

$-nk$$k$$k$$k$$k$$\cdots$$k$$0$$0$$0$$0$$0$$0$$-nk$$-nkkk$$\cdots$Complete
graph $K_{\frac{n}{2}}$$\frac{n}{2}$ blue verticesFire every blue vertex $k$
times Fire the pink vertex $k$ times
$\frac{n}{2}k$$0$$0$$0$$0$$\cdots$0$0$$0$$0$$0$$0$$0$$-nkkk$$\cdots$$-\frac{n}{2}k$$-\frac{n}{2}k$$0$$0$$0$$0$$0$$\cdots$0${\scriptstyle-}k$${\scriptstyle-}k$${\scriptstyle-}k$${\scriptstyle-}k$${\scriptstyle-}k$${\scriptstyle-}k$0$-nkkk$$\cdots$

Let $C\xrightarrow{\Delta{(v_{0})}}X_{0}$ with $v_{0}=\begin{pmatrix}k\\\
\vdots\\\ k\\\ 0\\\ 0\\\ \vdots\\\ 0\end{pmatrix}.$ Here $v_{0}$ is the best
possible among all vectors in $v_{0}+ker\Delta$, since there are at most $n/2$
positive and negative components (cf lemma 2). Thus, to improve, we must find
a sequence to another stable divisor. Note that the vector
$v_{min}=\begin{pmatrix}0\\\ \vdots\\\ 0\\\ -k\\\ 0\\\ \vdots\\\ 0\\\
\end{pmatrix}$ satisfies $C\xrightarrow{\Delta v_{\text{min}}}X_{1}$ and we
have $\left|V_{\min}\right|_{1}=\dfrac{\left|v_{0}\right|_{1}}{\frac{n}{2}}$.

##  3 Acknowledgements

I would like to thank Matt Baker for his advice (both on a mathematical and a
non-mathematical level) and for suggesting this problem to me.

## References

  * [1] Matthew Baker.  The dollar game on a graph (with 9 surprises).  2010\.  URL: https://drive.google.com/file/d/1RMmfJ_0E6Gy-R6a59mkR_donKzSNHxJq/view.
  * [2] Matthew Baker and Serguei Norine.  Riemann-Roch and Abel-Jacobi theory on a finite graph.  Advances in Mathematics, 215(2):766–788, 2007.  URL: https://arxiv.org/abs/math/0608360.
  * [3] Caroline Klivans.  The Mathematics of Chip-firing.  CRC Press, 11 2018.  URL: https://www.dam.brown.edu/people/cklivans/Chip-Firing.pdf.