# Inverse-Weighted Survival Games

Xintian Han
NYU
<EMAIL_ADDRESS>
Mark Goldstein 11footnotemark: 1
NYU
<EMAIL_ADDRESS>
Aahlad Puli
NYU
<EMAIL_ADDRESS>Thomas Wies
NYU
<EMAIL_ADDRESS>
Adler J. Perotte
Columbia University
<EMAIL_ADDRESS>
Rajesh Ranganath
NYU
<EMAIL_ADDRESS>
Equal Contribution.

###### Abstract

Deep models trained through maximum likelihood have achieved state-of-the-art
results for survival analysis. Despite this training scheme, practitioners
evaluate models under other criteria, such as binary classification losses at
a chosen set of time horizons, e.g. Brier score (bs) and Bernoulli log
likelihood (bll). Models trained with maximum likelihood may have poor bs or
bll since maximum likelihood does not directly optimize these criteria.
Directly optimizing criteria like bs requires inverse-weighting by the
censoring distribution. However, estimating the censoring model under these
metrics requires inverse-weighting by the failure distribution. The objective
for each model requires the other, but neither are known. To resolve this
dilemma, we introduce Inverse-Weighted Survival Games. In these games,
objectives for each model are built from re-weighted estimates featuring the
other model, where the latter is held fixed during training. When the loss is
proper, we show that the games always have the true failure and censoring
distributions as a stationary point. This means models in the game do not
leave the correct distributions once reached. We construct one case where this
stationary point is unique. We show that these games optimize bs on
simulations and then apply these principles on real world cancer and
critically-ill patient data.

##  1 Introduction

Survival analysis is the modeling of time-to-event distributions and is widely
used in healthcare to predict time from diagnosis to death, risk of disease
recurrence, and changes in level of care. In survival data, events, known as
failures, are often right-censored, i.e., only a lower bound on the time is
observed, for instance, when a patient leaves a study before failing. Under
certain assumptions, maximum likelihood estimators are consistent for survival
modeling (Kalbfleisch and Prentice, 2002).

Recently, deep survival models have obtained state-of-the-art results
(Ranganath et al., 2016; Alaa and van der Schaar, 2017; Katzman et al., 2018;
Kvamme et al., 2019; Zhong and Tibshirani, 2019). Common among these are
discrete-time models (Yu et al., 2011; Lee et al., 2018; Fotso, 2018; Lee et
al., 2019; Ren et al., 2019; Kvamme and Borgan, 2019b; Kamran and Wiens, 2021;
Goldstein et al., 2020; Sloma et al., 2021) even when data are continuous
because they can borrow classification architectures and flexibly approximate
continuous densities (Miscouridou et al., 2018).

Though training is often based on maximum likelihood, criteria such as Brier
score (bs) and Bernoulli log likelihood (bll) have been used to evaluate
survival models (Haider et al., 2020). The bs and bll are classification
losses adapted for survival by treating the model as a binary classifier at
various time horizons (will the event occur before or after 5 years?) (Kvamme
and Borgan, 2019b; Lee et al., 2019; Steingrimsson and Morrison, 2020). Bs can
also be motivated by calibration (section 3) which is valuable because
survival probabilities are used to communicate risk (Sullivan et al., 2004).
However bs and bll are challenging to estimate because they require inverse
probability of censor-weighting (ipcw), which depends on the true censoring
distribution (Van der Laan et al., 2003).

Though consistent, due to finite data, maximum likelihood may lead to models
with poor bs and bll. But directly optimizing these criteria is challenging
because ipcw estimation requires solving an additional survival modeling
problem to estimate the unknown censoring distribution. This poses a re-
weighting dilemma: each model is required for training the other under these
criteria but neither are known.

To resolve the dilemma, we introduce Inverse-Weighted Survival Games for
training with respect to criteria such as bs and bll. We pose survival
analysis as a game with the failure and censoring models as players. Each
model’s loss is built from ipcw estimates featuring the other model. Inspired
by game theory (Neumann and Morgenstern, 2007; Letcher et al., 2019), we ask:
should the censoring model’s re-weighting role in the failure objective be
considered part of the censoring objective? We find the answer to be no. In
each step of training, each model follows gradients of its loss with the other
model held fixed to compute weights.

When the loss is proper (e.g. bs, bll) (Gneiting and Raftery, 2007), we show
that games have the true failure and censoring distributions as a stationary
point. This means the models in the game do not leave the correct
distributions once reached. We then describe one case where this stationary
point is unique. Finally, we show that inverse-weighted game training achieves
better bs and bll than maximum likelihood methods on simulations and real
world cancer and ill-patient data.111Code is available at
https://github.com/rajesh-lab/Inverse-Weighted-Survival-Games

##  2 Notation and background on ipcw

#### Notation.

Let $T$ be a failure time with cdf $F(t)=P(T\leq t)$ , density $f$, survival
function $\overline{F}=1-F$, and model $F_{\theta_{T}}$. Let $C$ be a
censoring time with cdf $G$, density $g$, $\overline{G}=1-G$, and model
$G_{\theta_{C}}$. This means $\overline{G}(t)=P(C>t)$. Let
$\overline{G}(t^{-})$ denote $P(C\geq t)$. We observe features $X$, time
$U=\min(T,C)$ and $\Delta=\mathbbm{1}\left[T\leq C\right]$. For discrete
models over $K$ times, let $\theta_{Tt}=P_{\theta}(T=t)$ and
$\theta_{Ct}=P_{\theta}(C=t)$.

#### Models.

We focus on deep discrete models like those studied in Lee et al. (2018);
Kvamme and Borgan (2019b). The model maps inputs $X$ to a categorical
distribution over times. When the observations are continuous, a
discretization scheme is necessary. Following Kvamme and Borgan (2019b);
Goldstein et al. (2020), we set bins to correspond to quantiles of observed
times. We represent all times by the lower boundary of their respective
interval.

#### Assumptions.

We assume i.i.d. data and random censoring: $T\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}C~{}|~{}X$ (Kalbfleisch and Prentice, 2002).
We also require the censoring positivity assumption (Gerds et al., 2013). Let
$f=dF$. Then:

$\displaystyle\exists\epsilon\quad\text{ s.t. }\quad\forall x\,\forall
t\in\\{t\leq
t_{\text{max}}~{}|~{}f(t|x)>0\\},\quad\overline{G}(t^{-}|x)\geq\epsilon>0,$
(1)

i.e. it is possible that censoring events occur late-enough for us to observe
failures up until a maximum time $t_{\text{max}}$. Truncating at a maximum
time is necessary in practice for continuous distributions because datasets
may have no samples in the tails, leading to practical positivity violations
(Gerds et al., 2013). This truncation happens implicitly for categorical
models by choosing the last bin.

In this work, we model the censoring distribution. This task is dual to the
original survival problem: the roles of censoring and failure times are
reversed. Therefore, to observe censoring events properly, we also require a
version of eq. 1 to hold with the roles of $F$ and $G$ reversed (appendix A).

####  Ipcw estimators.

Inverse probability of censor-weighting (ipcw) is a method for estimation
under censoring (Van der Laan et al., 2003; Bang and Robins, 2005). Consider
the marginal mean $\mathop{\mathbb{E}}[T]$. Ipcw reformulates such
expectations in terms of observed data. Using ipcw, we can show that:

$\displaystyle\mathop{\mathbb{E}}[T]=\mathop{\mathbb{E}}_{X}\mathop{\mathbb{E}}_{T|X}\Big{[}\frac{\mathop{\mathbb{E}}[\mathbbm{1}\left[T\leq
C|X\right]]}{\mathop{\mathbb{E}}[\mathbbm{1}\left[T\leq
C|X\right]]}T\Big{]}=\mathop{\mathbb{E}}_{X}\mathop{\mathbb{E}}_{T|X}\mathop{\mathbb{E}}_{C|X}\Bigg{[}\frac{\mathbbm{1}\left[T\leq
C\right]}{\overline{G}(T^{-}|X)}T\Bigg{]}=\mathop{\mathbb{E}}_{T,C,X}\Bigg{[}\frac{\Delta
U}{\overline{G}(U^{-}|X)}\Bigg{]}$

We derive this fully in appendix B. The second equality holds because
$\Delta=1\implies U=T$ and means we can identify $\mathop{\mathbb{E}}[T|X]$
provided that we know $G$ and that random censoring and positivity hold.

##  3 Time-dependent survival evaluations

Figure 1: Test nll and bs for 5 different models, each trained with a
different learning rate

Brier score (bs) (Brier and Allen, 1951) is proper for classification, meaning
that it has a minimum at the true data distribution (Gneiting and Raftery,
2007). The bs is often adapted for survival evaluations (Lee et al., 2019;
Kvamme et al., 2019; Haider et al., 2020). For time $t$, it computes
differences between the cdf and true event status at $t$, turning survival
analysis into a classification problem at a given time horizon:

$\displaystyle\text{BS}(t;\theta)$
$\displaystyle=\mathop{\mathbb{E}}\Big{[}\Big{(}F_{\theta_{T}}(t~{}|~{}X)-\mathbbm{1}\left[T\leq
t\right]\Big{)}^{2}\Big{]}$ (2)

Bs is often used as a proxy for marginal calibration error (Kumar et al.,
2018; Lee et al., 2019), which measures differences between cdf levels
$\alpha\in[0,1]$ and observed proportions of datapoints with
$F_{\theta}(T|X)<\alpha$ (Demler et al., 2015). This usage of bs stems from
its decomposition into calibration plus a refinement (discriminative) term
(DeGroot and Fienberg, 1983).

Unfortunately one cannot compute bs unmodified since $\mathbbm{1}\left[T\leq
t\right]$ is unobserved for a point censored before $t$. Ipcw bs (Graf et al.,
1999; Gerds and Schumacher, 2006) estimates bs $(t)$ under censoring:

$\displaystyle\text{BS}(t;\theta)$
$\displaystyle=\mathop{\mathbb{E}}\Big{[}\frac{\overline{F}_{\theta_{T}}(t~{}|~{}X)^{2}\Delta\mathbbm{1}\left[U\leq
t\right]}{\overline{G}(U^{-}~{}|~{}X)}+\frac{F_{\theta_{T}}(t~{}|~{}X)^{2}\mathbbm{1}\left[U>t\right]}{\overline{G}(t~{}|~{}X)}\Big{]}.$
(3)

eq. 3 is equivalent to eq. 2 (appendix C). Negative Bernoulli log likelihood
(bll) is similar, but with log loss (appendix D). bs and bll are proper for
classification at each time $t$, so their sum or integral over $t$ is still
proper (appendix H).

#### Proper objectives differ.

Though negative log likelihood (nll), bs and bll all have the same true
distribution at optimum with infinite data, they may yield significantly
different solutions in practice. For example, nll-trained models may not
achieve good bs (Kvamme and Borgan, 2019a). In fig. 1, we show test set nll
and bs for 5 models trained with nll at different learning rates on Gamma-
simulated data (described in section 5.1). Nll does not align with bs: models
that have low nll may not have low bs. Model 4 has the lowest nll but not the
lowest bs. When a practitioner requires good performance under bs or bll, they
should optimize directly for those metrics.

#### Re-weighting dilemma.

Censoring introduces challenges because we must use ipcw to estimate bs and
bll. Crucially, the $G$ in eq. 3 is the true censoring distribution rather
than a model, but during training, we only have access to models. This poses a
dilemma: can the models be used in re-weighting estimates during training to
successfully optimize these criteria under censoring?

##  4 Inverse-Weighted Survival Games

A reasonable attempt to solve the dilemma is to jointly optimize the sum of
$F_{\theta}$ and $G_{\theta}$’s losses where each model re-weights the other’s
loss. The expectation is that both models will improve over training and yield
reliable ipcw estimates for each other. Concretely, consider this for eq. 3
plus the same objective with the roles of $F_{\theta}$ and $G_{\theta}$
reversed. Unfortunately, there exist solutions to this optimization problem
with smaller loss than for the distributions from which the data was
generated, making this summed objective improper for the pair of
distributions. In fig. 2, we plot this for ipcw bs $(t=1)$ for models over two
timesteps222$\lx@glossaries@gls@link{acronym}{bs}{\leavevmode\textsc{bs}}(t=1)$
is proper for distributions with support over two timesteps because
$\lx@glossaries@gls@link{acronym}{bs}{\leavevmode\textsc{bs}}(t=K)$ for a
model with support over $K$ timesteps is always $0$, so the summed bs equals
$\lx@glossaries@gls@link{acronym}{bs}{\leavevmode\textsc{bs}}(1)$. as a
function of each model’s single parameter.

(a) Contours of summed objective

(b) Gradients of the game

Figure 2:  Figure 2(a): the sum of $F$ and $G$’s ipcw bs (1) scores is an
improper joint objective for the failure and censoring models. Figure 2(b): in
contrast, as shown in the gradient field, the one timestep game has a unique
stationary point at the true data generating distributions.

To address this phenomenon, we introduce Inverse-Weighted Survival Games. In
these games, a failure player and censor player simultaneously minimize their
own loss function. The failure and censoring model are featured in both loss
functions and playing the game results in a trained failure and censoring
model. We show in experiments that these games produce models with good bs,
bll, and concordance relative to those trained with maximum likelihood.

For simplicity, we present the games for marginal categorical models. The
analysis can be extended to conditional parameterizations with the usual
caveats shared by maximum likelihood. Our experiments explore the conditional
setting.

###  4.1 Basic definition of game

We follow the setup in Letcher et al. (2019). A differentiable $n$-player game
consists of $n$ players each with loss $\ell_{i}$ and parameter (or state)
$\theta_{i}$. Player $i$ controls only $\theta_{i}$ and aims to minimize
$\ell_{i}$. However, each $\ell_{i}$ is a function of the whole state
$\theta=(\theta_{i},\theta_{-i})$333For two players, when $\ell_{1}=-\ell_{2}$
the game is called zero-sum and can be written as a minimax game, sometimes
referred to as adversarial (e.g. as in Generative Adversarial Networks (gans)
(Goodfellow et al., 2014)). Games with $\ell_{1}\neq-\ell_{2}$ are called non-
zero-sum. . The simultaneous gradient is the gradient of the losses w.r.t.
each players’ respective parameters:

$\displaystyle\xi(\theta)=[\nabla_{\theta_{1}}\ell_{1},\ldots,,\nabla_{\theta_{n}}\ell_{n}]$

The dynamics of the game refers to following $-\xi$. The solution concepts in
games are equilibria (the game analog of optima) and stationary points. One
necessary condition for equilibria is finding a stationary point
$\theta^{\star}$ such that $\xi(\theta^{\star})=0$. The simplest algorithm
follows the dynamics to find stationary points and is called simultaneous
gradient descent. With learning rate $\eta$,

$\displaystyle\theta\leftarrow\theta-\eta\xi(\theta).$

This can be interpreted as each player taking their best possible move at each
instant.

###  4.2 Constructing survival games

We specify an Inverse-Weighted Survival Game as follows. First, choose a loss
$L$ used to construct the losses for each player. Next, derive the ipcw form
$L_{I}$ that can be used to compute $L$ under censoring: for true failure and
censoring distributions $F^{\star}$ and $G^{\star}$, the losses $L$ and
$L_{I}$ are related through
$L_{I}(F_{\theta_{T}};G^{\star})=L(F_{\theta_{T}})$ and
$L_{I}(G_{\theta_{C}};F^{\star})=L(G_{\theta_{C}})$. The loss functions for
the two players are defined as:

$\displaystyle\begin{split}\ell_{F}(\theta)\triangleq
L_{I}(F_{\theta_{T}};G_{\theta_{C}}),\quad\ell_{G}(\theta)\triangleq
L_{I}(G_{\theta_{C}};F_{\theta_{T}})\end{split}$ (4)

Compared to eq. 3, we have replaced the true re-weighting distributions with
models. Finally, the failure player and censor player minimize their
respective loss functions w.r.t. only their own parameters:

$\displaystyle\texttt{failure
player:}\min_{\theta_{T}}\ell_{F},\quad\quad\texttt{censor
player:}\min_{\theta_{C}}\ell_{G}$

One example of these games is the ipcw bs $(t)$ game, derived in appendix C.
With $\overline{\Delta}=1-\Delta$,

$\displaystyle\begin{split}\ell_{F}^{t}(\theta)&=\mathop{\mathbb{E}}\Big{[}\frac{\overline{F}_{\theta_{T}}(t)^{2}\Delta\mathbbm{1}\left[U\leq
t\right]}{\overline{G}_{\theta_{C}}(U^{-})}+\frac{F_{\theta_{T}}(t)^{2}\mathbbm{1}\left[U>t\right]}{\overline{G}_{\theta_{C}}(t)}\Big{]}\\\
\ell_{G}^{t}(\theta)&=\mathop{\mathbb{E}}\Big{[}\frac{\overline{G}_{\theta_{C}}(t)^{2}\overline{\Delta}\mathbbm{1}\left[U\leq
t\right]}{\overline{F}_{\theta_{T}}(U)}+\frac{G_{\theta_{C}}(t)^{2}\mathbbm{1}\left[U>t\right]}{\overline{F}_{\theta_{T}}(t)}\Big{]}.\end{split}$
(5)

In section 4.3, we show that this formulation (fig. 2(b)) has formal
advantages over the optimization in fig. 2(a) for particular choices of $L$.

#### Multiple Timesteps.

The example is specified for a given $t$, but the games can be designed for
multiple timesteps. We use bs for a $K$ timestep model to demonstrate.
$\lx@glossaries@gls@link{acronym}{bs}{\leavevmode\textsc{bs}}(K)$ is $0$ for
any model: the left terms contain $\overline{F}_{\theta_{T}}$ and
$\overline{G}_{\theta_{C}}$, which are are both $0$ when evaluated at $K$; in
the right terms, $\mathbbm{1}\left[U>K\right]$ is always $0$. One option is to
define summed games with:

$\displaystyle\ell_{F}=\sum_{t=1}^{K-1}\ell_{F}^{t},\quad\ell_{G}=\sum_{t=1}^{K-1}\ell_{G}^{t}$

The summed game is shown in algorithm 1. Alternatively, instead of the sum, it
is possible to find solutions for all timesteps with respect to one pair of
models $(F_{\theta},G_{\theta})$. For $K\text{-}1$ timesteps this can be
formalized as a $2(K\text{-}1)$-player game: there is a failure player and
censor player for the loss at each $t$:

$\displaystyle t^{th}\text{-}\texttt{failure
player:}\min_{\theta_{Tt}}\ell_{F}^{t},\quad\quad t^{th}\text{-}\texttt{censor
player:}\min_{\theta_{Ct}}\ell_{G}^{t}$

We study theory that applies to both approaches in section 4.3, namely that
the true failure and censoring distribution are stationary points in both
types of games. We prove additional properties about uniqueness of the
stationary point for a special case of the multi-player game in section 4.4.
Summed games are more stable to optimize in practice because they optimize
objectives at all time steps w.r.t all parameters, while multiplayer games can
only improve each time step’s loss w.r.t. one parameter. We study the summed
games empirically in section 5.

Input: Choice of losses $\ell_{F},\ell_{G}$, learning rate $\gamma$

Initialize: $T$ model parameters $\theta_{T}$ and $C$ model parameters
$\theta_{C}$ randomly

repeat

Set $g_{T}=0$ and $g_{C}=0$

for $t=1$ to $K-1$ do

$g_{T}=g_{T}+d\ell_{F}^{t}/d\theta_{T}$

$g_{C}=g_{C}+d\ell_{G}^{t}/d\theta_{C}$

end for

$\theta_{T}\leftarrow\theta_{T}-\gamma g_{T}$ and
$\theta_{C}\leftarrow\theta_{C}-\gamma g_{C}$

until convergence

Output: $\theta_{T},\theta_{C}$

Algorithm 1 Following Gradients in Summed Games

###  4.3 Ipcw games have a stationary point at data distributions

Among a game’s stationary points should be the true failure and censoring
distributions.

######  Proposition 1.

Assume
$\exists\theta_{T}^{\star}\in\Theta_{T},\exists\theta_{C}^{\star}\in\Theta_{C}$
such that $F^{\star}=F_{\theta_{T}^{\star}}$ and
$G^{\star}=G_{\theta_{C}^{\star}}$. Assume the game losses $\ell_{F},\ell_{G}$
are based on proper losses $L$ and that the games are only computed at times
for which positivity holds. Then $(\theta_{T}^{\star},\theta_{C}^{\star})$ is
a stationary point of the game eq. 4.

The proof is in appendix I. The result holds for summed and multi-player games
using bs, bll, or other proper scoring rules such as auc.444Though often
reported, the time-dependent concordance$(t)$ is not proper (Blanche et al.,
2019). When games are built from such objectives, the set of solutions
includes $(\theta_{T}^{\star},\theta_{C}^{\star})$ and models do not leave
this correct solution when reached. Under the stated assumptions, this result
holds for discrete and continuous distributions. However, as mentioned in
section 2, in practice a truncation time must be picked to ensure the
assumptions are met for continuous distributions.

###  4.4 Uniqueness for Discrete Brier Games

We provide a stronger result for the bs game in eq. 5 when solving all
timesteps with multi-player games: its only stationary point is located at the
true failure and censoring distributions.

######  Proposition 2.

For discrete models over $K$ timesteps, assuming that $\theta^{\star}_{Tt}>0$
and $\theta^{\star}_{Ct}>0$, the solution
$(\theta_{T}^{\star},\theta_{C}^{\star})$ is the only stationary point for the
multi-player bs game shown in algorithm 2 for times $t\in\\{1,\ldots,K-1\\}$

The proof is in appendix J. To illustrate this, fig. 2(b) shows that, unlike
the minimization in fig. 2(a), the ipcw bs game moves to the correct solution
at its unique stationary point.

##  5 Experiments

We run experiments on a simulation with conditionally Gamma times, a semi-
simulated survival dataset based on mnist, several sources of cancer data, and
data on critically-ill hospital patients.

#### Losses.

We build categorical models in 3 ways: the standard nll method (eq. 8), the
ipcw bs game and the negative ipcw bll game.

#### Metrics.

For these models we report bs (uncensored for simulations and Kaplan-Meier
(km)-weighted for real data), bll (also uncensored or weighted), concordance
which measures the proportion of pairs whose predicted risks are ordered
correctly (Harrell Jr et al., 1996), and nll. We report mean and standard
deviation of the metrics over 5 different seeds. In all plots, the middle
solid line represents the mean and the shaded band represents the standard
deviation.

#### Model Description.

In all experiments except for mnist, we use a 3-hidden-layer ReLU network that
outputs $20$ categorical bins (more bin choices in section G.1). For mnist we
first use a small convolutional network and follow with the same fully-
connected network.

Model and training details including learning rate and model selection can be
found in appendix F.

###  5.1 Simulation Studies

#### Data.

We draw $X\in\mathbb{R}^{32}\sim\mathcal{N}(0,10I)$ and
$T\sim\text{Gamma}(\text{mean}=\mu_{t})$ where $\mu_{t}$ a log-linear function
of $X$. The censoring times are also gamma with mean $0.9*\mu_{t}$. Both
distributions have constant variance $0.05$. It holds that
$T\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}C~{}|~{}X$. Each random seed draws a new
dataset.

#### Results.

fig. 3 demonstrates that the games optimize the true uncensored bs, and,
though more slowly with respect to training size, log-likelihood does too. The
games have better test-set performance on all metrics for small training size.
All methods converge on similar performance when there is enough data (though
enough is highly-dependent on dimensionality and model class).

(a) Uncensored bs (b) Uncensored Neg bll (c) Concordance (d) Categorical nll

Figure 3: Test set evaluation metrics (y-axis) on the Gamma simulation versus
number of training points (x-axis) for three methods. Each point in the plot
represents the evaluation metric value of a fully trained model with that
number of training points. Lower is better for all the metrics except for
concordance.

#### Calibration.

We include a qualitative investigation of model calibration on the gamma
simulation trained with 2000 datapoints. fig. 4 shows that the bs game
achieves near-perfect calibration while the two likelihood-based methods
suffer some error. This is expected since likelihood-based methods do not
directly optimize calibration while bs does (section 3).

Figure 4: Calibration curves (Avati et al., 2019) comparing game-training and
nll-training. .

###  5.2 Semi-simulated studies

#### Data.

Survival-mnist (Gensheimer, 2019; Pölsterl, 2019) draws times conditionally on
mnist label $Y$. This means digits define risk groups and
$T\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X~{}|~{}Y$. Times within a digit are i.i.d.
The model only sees the image pixels $X$ as covariates so it must learn to
classify digits (risk groups) to model times. We follow Goldstein et al.
(2020) and use Gamma times. $T\sim\text{Gamma}(\text{mean}=10*(Y+1))$. We set
the variance constant to $0.05$. Lower digit labels $Y$ yield earlier event
times. $C$ is drawn similarly but with mean $9.9*(Y+1)$. Each random seed
draws a new dataset.

#### Results.

This experiment demonstrates that better nll does not correspond to better
performance on bs, bll, and concordance. Similarly to the previous experiment,
fig. 5 shows that game methods attain better uncensored test-set bs and bll on
survival-mnist than likelihood-based training does. The games likewise attain
higher concordance. Nll training performs better at the metric it directly
optimizes. This experiment also establishes that it is possible to optimize
through deep convolutional models with batch norm and pooling using the game
training methods.

(a) Uncensored bs (b) Uncensored Neg bll (c) Concordance (d) Categorical nll

Figure 5: Test set evaluation metrics (y-axis) on survival-mnist versus number
of training points (x-axis) for three methods. Each point in the plot
represents the evaluation metric value of a fully trained model with that
number of training points. Lower is better for all the metrics except for
concordance.

###  5.3 Real Datasets

#### Data.

We use several datasets used in recent papers (Chen, 2020; Kvamme et al.,
2019) and available in the python packages DeepSurv (Katzman et al., 2018) and
PyCox (Kvamme et al., 2019), and the R Survival (Therneau, 2021). The datasets
are:

  * •

Molecular Taxonomy of Breast Cancer International Consortium (metabric)
(Curtis et al., 2012)

  * •

Rotterdam Tumor Bank (rott) (Foekens et al., 2000) and German Breast Cancer
Study Group (gbsg) (Schumacher et al., 1994) combined into one dataset (rott.
& gbsg)

  * •

Study to Understand Prognoses Preferences Outcomes and Risks of Treatment
(support) (Knaus et al., 1995) which includes severely ill hospital patients

For more description see appendix F. For real data, there is no known ground
truth for the censoring distribution, which means evaluation requires
assumptions. Following the experiments in Kvamme et al. (2019), we assume that
censoring is marginal estimate with km to evaluate models.555This is also the
route taken in the R packages Survival (Therneau, 2021), PEC (Mogensen et al.,
2012), and riskRegression (Gerds et al., 2020).

#### Results.

On metabric, games attain lower (better) km-weighted bs and bll than nll-
training when the number of datapoints is small, and have better concordance
and nll though they do not directly optimize them. As data size increases, all
methods converge to similar performance. On rott. & gbsg, the trend is
similar: games optimize the bs and bll more rapidly as a function of training
set size than nll-training does. Again, all methods converge to similar
performance in all metrics when the number of datapoints is large enough. All
methods perform similarly on support.

#### Caveats.

First, though popular, these survival datasets are low-dimensional, so any of
the objectives can perform well on the metrics with just several hundred
points. We see that this is distinct from mnist, where thousands of points
were required to improve performance. Second, though possible, it may not be
true that censoring is marginal on these datasets, which would mean that the
bs and bll results only have their interpretation conditional on a particular
set of covariates. Our method is also correct when the censoring is marginal
though. Lastly, no method is stable for all metrics, for all training sizes,
on all seeds for all datasets.

(a) km-weighted bs (b) km-weighted Neg bll (c) concordance (d) categorical nll

Figure 6:  Test set evaluation metrics (y-axis) on metabric versus number of
training points (x-axis) for three methods. Each point in the plot represents
the evaluation metric value of a fully trained model with that number of
training points. Lower is better for all the metrics except for concordance.

(a) km-weighted bs (b) km-weighted Neg bll (c) concordance (d) categorical nll

Figure 7:  Test set evaluation metrics (y-axis) on rott. & gbsg versus number
of training points (x-axis) for three methods. Each point in the plot
represents the evaluation metric value of a fully trained model with that
number of training points. Lower is better for all the metrics except for
concordance.

(a) km-weighted bs (b) km-weighted Neg bll (c) concordance (d) categorical nll

Figure 8:  Test set evaluation metrics (y-axis) on support versus number of
training points (x-axis) for three methods. Each point in the plot represents
the evaluation metric value of a fully trained model with that number of
training points. Lower is better for all the metrics except for concordance.

##  6 Related Work

#### Nuisance parameters.

Under non-informative censoring, the censoring distribution is unrelated to
the failure distribution, but estimating it can help improve learning the
failure distribution; here, the censoring distribution is a nuisance
parameter. Existing causal estimation methods propose two-stage procedures
where the first stage estimates the nuisance-parameter (e.g. propensity score)
and the second stage uses the learned nuisance-parameter as-is to define an
estimator or loss function for the target parameter (causal effect). (Van Der
Laan and Rubin, 2006; Van der Laan and Rose, 2011; Chernozhukov et al., 2018;
Foster and Syrgkanis, 2019). In this work, we instead show that estimating the
target (failure model objective or failure model itself) can benefit from a
coupled estimation procedure where the nuisance parameter (censoring model) is
also trained simultaneously. The failure model needs the censoring
distribution to compute bs but censoring estimation needs the failure model,
and despite this circular dependence, we characterize a case where the game
training leads to the true data generating distributions.

#### Double Robust Censoring Unbiased Transformations.

For functions $h$, Rubin and van der Laan (2007) estimate conditional mean
$\mathop{\mathbb{E}}[h(T,X)|X]$ under censoring using a double-robust
estimator: given estimates of the conditional failure and censoring cumulative
distribution functions (cdfs) $\hat{F}(t|X)$ and $\hat{G}(c|X)$, the estimator
of $\mathop{\mathbb{E}}[h(T,X)|X]$ is unbiased when either nuisance cdf is
correct. However, here we are concerned with estimating a quantity to be used
as a loss for learning $\hat{F}$. We therefore presumably do not already have
an estimate of $\hat{F}$ to be used in a doubly-robust estimator.

#### Censoring Unbiased Losses for Deep Learning.

Steingrimsson and Morrison (2020) build failure model loss functions based on
the estimators from Rubin and van der Laan (2007). Their bs loss extends ipcw
bs estimation to the doubly-robust case and to our knowledge is the first
instance of ipcw-based estimation procedures being used in a general purpose
way to define loss functions for deep survival analysis.

However, their censoring distribution is estimated once before training and
held fixed rather than incorporated into a joint training procedure as in the
games in this work. The fixed censoring estimate is implemented by km, which
assumes a marginal censoring distribution. Making the marginal assumption when
censoring is truly conditional should not yield a performant model under the
bs criteria since the training objective does not directly estimate or
optimize the true bs that would be measured under no censoring. When marginal
censoring does hold, km estimation, which is non-parametric, may be a simpler
and stable choice versus the game, depending on sample size, data variance,
and conditional parameterization assumptions. But since it is in general
unknown if censoring is marginal, we use conditional models which are also
correct under marginal censoring.

##  7 Discussion

In this work, we propose a new training method for survival models under
censored data. We argue that on finite data, it is important to close the gap
between training methodology and the desired evaluation criteria. We showed in
the experiments that better nll does not correspond to better performance on
bs, bll, and concordance, all evaluations of interest in survival analysis.

The main trend in our experimental results was that data size matters: smaller
meant the game methods performed better than nll and enough data meant that
they perform similarly, which is expected since all objectives are proper.
However enough data is hard to define: it depends on dimensionality and on the
data generating distribution and model class. It is a great direction to build
more precise understanding on how objectives behave differently even when they
have the same optimum on infinite data: though likelihood is known to be
asymptotically efficient for survival analysis, more analysis is necessary for
comparing likelihood and Brier score’s trade-offs on small sample sizes.

In the experiments, we focus on categorical models. On the other hand,
proposition 1 applies to continuous distributions as well, provided that
positivity can be satisfied. However, this is rare in practice because most
survival data has a final follow-up time, and even before this time there may
be very few samples with late times (Gerds et al., 2013). For this reason,
working with continuous distributions requires picking a truncation time and
playing games only up to that time.

Evaluation on real data under censoring requires assumptions. It is important
to further consider how to better assess test-set performance on metrics such
as bs, bll, and concordance. Because concordance is not proper (Blanche et
al., 2019), we do not build objectives from it here, but it too is not
invariant to censoring. Regarding games, we showed properties about stationary
points. More analysis is necessary to describe important convergence
properties of optimizing these games.

#### Social Impact.

Survival models are deployed in hospital settings and have high impact on
public health. In this work, we saw benefits of a new training approach for
these models, but no training method is a panacea. Practitioners of survival
analysis must take great care to consider various training and validation
approaches, as well as consider possible test distribution shifts, prior to
deployment.

## Acknowledgments and Disclosure of Funding

This work was supported by:

  * •

NIH/NHLBI Award R01HL148248

  * •

NSF Award 1922658 NRT-HDR: FUTURE Foundations, Translation, and Responsibility
for Data Science.

  * •

NSF Award 1514422 TWC: Medium: Scaling proof-based verifiable computation

  * •

NSF Award 1815633 SHF

## References

  * Alaa and van der Schaar (2017) A. M. Alaa and M. van der Schaar.  Deep multi-task gaussian processes for survival analysis with competing risks.  In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , pages 2326–2334. Curran Associates Inc., 2017\.
  * Andersen et al. (2012) P. K. Andersen, O. Borgan, R. D. Gill, and N. Keiding.  _Statistical models based on counting processes_.  Springer Science & Business Media, 2012.
  * Avati et al. (2019) A. Avati, T. Duan, S. Zhou, K. Jung, N. H. Shah, and A. Y. Ng.  Countdown regression: Sharp and calibrated survival predictions.  In A. Globerson and R. Silva, editors, _Proceedings of the Thirty-Fifth Conference on Uncertainty in Artificial Intelligence, UAI 2019, Tel Aviv, Israel, July 22-25, 2019_ , page 28. AUAI Press, 2019.  URL http://auai.org/uai2019/proceedings/papers/28.pdf.
  * Bang and Robins (2005) H. Bang and J. M. Robins.  Doubly robust estimation in missing data and causal inference models.  _Biometrics_ , 61(4):962–973, 2005.
  * Blanche et al. (2013) P. Blanche, J.-F. Dartigues, and H. Jacqmin-Gadda.  Review and comparison of roc curve estimators for a time-dependent outcome with marker-dependent censoring.  _Biometrical Journal_ , 55(5):687–704, 2013.
  * Blanche et al. (2019) P. Blanche, M. W. Kattan, and T. A. Gerds.  The c-index is not proper for the evaluation of-year predicted risks.  _Biostatistics_ , 20(2):347–357, 2019.
  * Brier and Allen (1951) G. W. Brier and R. A. Allen.  Verification of weather forecasts.  In _Compendium of meteorology_ , pages 841–848. Springer, 1951.
  * Chen (2020) G. H. Chen.  Deep kernel survival analysis and subject-specific survival time prediction intervals.  In _Machine Learning for Healthcare Conference_ , pages 537–565. PMLR, 2020.
  * Chernozhukov et al. (2018) V. Chernozhukov, D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, and J. Robins.  Double/debiased machine learning for treatment and structural parameters, 2018.
  * Curtis et al. (2012) C. Curtis, S. P. Shah, S.-F. Chin, G. Turashvili, O. M. Rueda, M. J. Dunning, D. Speed, A. G. Lynch, S. Samarajiwa, Y. Yuan, et al.  The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups.  _Nature_ , 486(7403):346–352, 2012.
  * DeGroot and Fienberg (1983) M. H. DeGroot and S. E. Fienberg.  The comparison and evaluation of forecasters.  _Journal of the Royal Statistical Society: Series D (The Statistician)_ , 32(1-2):12–22, 1983.
  * Demler et al. (2015) O. V. Demler, N. P. Paynter, and N. R. Cook.  Tests of calibration and goodness-of-fit in the survival setting.  _Statistics in medicine_ , 34(10):1659–1680, 2015\.
  * Foekens et al. (2000) J. A. Foekens, H. A. Peters, M. P. Look, H. Portengen, M. Schmitt, M. D. Kramer, N. Brünner, F. Jänicke, M. E. Meijer-van Gelder, S. C. Henzen-Logmans, et al.  The urokinase system of plasminogen activation and prognosis in 2780 breast cancer patients.  _Cancer research_ , 60(3):636–643, 2000.
  * Foster and Syrgkanis (2019) D. J. Foster and V. Syrgkanis.  Orthogonal statistical learning.  _arXiv preprint arXiv:1901.09036_ , 2019.
  * Fotso (2018) S. Fotso.  Deep neural networks for survival analysis based on a multi-task framework.  _arXiv preprint arXiv:1801.05512_ , 2018.
  * Gensheimer (2019) B. Gensheimer, Michael F.and Narasimhan.  A scalable discrete-time survival model for neural networks.  _PeerJ 7:e6257_ , 2019.
  * Gerds and Schumacher (2006) T. A. Gerds and M. Schumacher.  Consistent estimation of the expected brier score in general survival models with right-censored event times.  _Biometrical Journal_ , 48(6):1029–1040, 2006\.
  * Gerds et al. (2013) T. A. Gerds, M. W. Kattan, M. Schumacher, and C. Yu.  Estimating a time-dependent concordance index for survival prediction models with covariate dependent censoring.  _Statistics in Medicine_ , 32(13):2173–2184, 2013\.
  * Gerds et al. (2020) T. A. Gerds, P. Blanche, T. H. Scheike, R. Mortensen, M. Wright, N. Tollenaar, J. Muschelli, U. B. Mogensen, and B. Ozenne.  Package ‘riskregression’, 2020.
  * Gneiting and Raftery (2007) T. Gneiting and A. E. Raftery.  Strictly proper scoring rules, prediction, and estimation.  _Journal of the American statistical Association_ , 102(477):359–378, 2007.
  * Goldstein et al. (2020) M. Goldstein, X. Han, A. M. Puli, A. Perotte, and R. Ranganath.  X-cal: Explicit calibration for survival analysis.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Goodfellow et al. (2014) I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio.  Generative adversarial networks.  _arXiv preprint arXiv:1406.2661_ , 2014.
  * Graf et al. (1999) E. Graf, C. Schmoor, W. Sauerbrei, and M. Schumacher.  Assessment and comparison of prognostic classification schemes for survival data.  _Statistics in medicine_ , 18(17-18):2529–2545, 1999.
  * Haider et al. (2020) H. Haider, B. Hoehn, S. Davis, and R. Greiner.  Effective ways to build and evaluate individual survival distributions.  _Journal of Machine Learning Research_ , 21(85):1–63, 2020.
  * Harrell Jr et al. (1996) F. E. Harrell Jr, K. L. Lee, and D. B. Mark.  Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors.  _Statistics in medicine_ , 15(4):361–387, 1996\.
  * Hung and Chiang (2010a) H. Hung and C.-T. Chiang.  Estimation methods for time-dependent auc models with survival data.  _Canadian Journal of Statistics_ , 38(1):8–26, 2010a.
  * Hung and Chiang (2010b) H. Hung and C.-t. Chiang.  Optimal composite markers for time-dependent receiver operating characteristic curves with censored survival data.  _Scandinavian journal of statistics_ , 37(4):664–679, 2010b.
  * Kalbfleisch and Prentice (2002) J. D. Kalbfleisch and R. L. Prentice.  _The Statistical Analysis of Failure Time Data_.  Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., 2 edition, 2002.
  * Kamran and Wiens (2021) F. Kamran and J. Wiens.  Estimating calibrated individualized survival curves with deep learning.  2021\.
  * Katzman et al. (2018) J. L. Katzman, U. Shaham, A. Cloninger, J. Bates, T. Jiang, and Y. Kluger.  Deepsurv: personalized treatment recommender system using a cox proportional hazards deep neural network.  _BMC medical research methodology_ , 18(1):24, 2018.
  * Knaus et al. (1995) W. A. Knaus, F. E. Harrell, J. Lynn, L. Goldman, R. S. Phillips, A. F. Connors, N. V. Dawson, W. J. Fulkerson, R. M. Califf, N. Desbiens, et al.  The support prognostic model: Objective estimates of survival for seriously ill hospitalized adults.  _Annals of internal medicine_ , 122(3):191–203, 1995.
  * Kumar et al. (2018) A. Kumar, S. Sarawagi, and U. Jain.  Trainable calibration measures for neural networks from kernel mean embeddings.  In _International Conference on Machine Learning_ , pages 2805–2814, 2018.
  * Kvamme and Borgan (2019a) H. Kvamme and Ø. Borgan.  The brier score under administrative censoring: Problems and solutions.  _arXiv preprint arXiv:1912.08581_ , 2019a.
  * Kvamme and Borgan (2019b) H. Kvamme and Ø. Borgan.  Continuous and discrete-time survival prediction with neural networks.  _arXiv preprint arXiv:1910.06724_ , 2019b.
  * Kvamme et al. (2019) H. Kvamme, Ørnulf Borgan, and I. Scheel.  Time-to-event prediction with neural networks and cox regression.  _Journal of Machine Learning Research_ , 20(129):1–30, 2019.  URL http://jmlr.org/papers/v20/18-424.html.
  * Lee et al. (2018) C. Lee, W. R. Zame, J. Yoon, and M. van der Schaar.  Deephit: A deep learning approach to survival analysis with competing risks.  In _Thirty-Second AAAI Conference on Artificial Intelligence_ , 2018\.
  * Lee et al. (2019) C. Lee, W. Zame, A. Alaa, and M. Schaar.  Temporal quilting for survival analysis.  In _The 22nd international conference on artificial intelligence and statistics_ , pages 596–605. PMLR, 2019.
  * Letcher et al. (2019) A. Letcher, D. Balduzzi, S. Racaniere, J. Martens, J. Foerster, K. Tuyls, and T. Graepel.  Differentiable game mechanics.  _The Journal of Machine Learning Research_ , 20(1):3032–3071, 2019.
  * Miscouridou et al. (2018) X. Miscouridou, A. Perotte, N. Elhadad, and R. Ranganath.  Deep survival analysis: Nonparametrics and missingness.  In _Machine Learning for Healthcare Conference_ , pages 244–256, 2018\.
  * Mogensen et al. (2012) U. B. Mogensen, H. Ishwaran, and T. A. Gerds.  Evaluating random forests for survival analysis using prediction error curves.  _Journal of Statistical Software_ , 50(11):1–23, 2012.  URL https://www.jstatsoft.org/v50/i11.
  * Neumann and Morgenstern (2007) J. v. Neumann and O. Morgenstern.  _Theory of Games and Economic Behavior_.  Princeton University Press, 2007.
  * Pölsterl (2019) S. Pölsterl.  Sebastian pölsterl, Jul 2019.  URL https://k-d-w.org/blog/2019/07/survival-analysis-for-deep-learning/.
  * Ranganath et al. (2016) R. Ranganath, A. Perotte, N. Elhadad, and D. Blei.  Deep survival analysis.  _arXiv preprint arXiv:1608.02158_ , 2016.
  * Ren et al. (2019) K. Ren, J. Qin, L. Zheng, Z. Yang, W. Zhang, L. Qiu, and Y. Yu.  Deep recurrent survival analysis.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pages 4798–4805, 2019.
  * Rubin and van der Laan (2007) D. Rubin and M. J. van der Laan.  A doubly robust censoring unbiased transformation.  _The international journal of biostatistics_ , 3(1), 2007\.
  * Schumacher et al. (1994) M. Schumacher, G. Bastert, H. Bojar, K. Hübner, M. Olschewski, W. Sauerbrei, C. Schmoor, C. Beyerle, R. Neumann, and H. Rauschecker.  Randomized 2 x 2 trial evaluating hormonal treatment and the duration of chemotherapy in node-positive breast cancer patients. german breast cancer study group.  _Journal of Clinical Oncology_ , 12(10):2086–2093, 1994.
  * Sloma et al. (2021) M. Sloma, F. Syed, M. Nemati, and K. S. Xu.  Empirical comparison of continuous and discrete-time representations for survival prediction.  In _Survival Prediction-Algorithms, Challenges and Applications_ , pages 118–131. PMLR, 2021.
  * Steingrimsson and Morrison (2020) J. A. Steingrimsson and S. Morrison.  Deep learning for survival outcomes.  _Statistics in medicine_ , 39(17):2339–2349, 2020\.
  * Sullivan et al. (2004) L. M. Sullivan, J. M. Massaro, and R. B. D’Agostino Sr.  Presentation of multivariate data for clinical use: The framingham study risk score functions.  _Statistics in medicine_ , 23(10):1631–1660, 2004\.
  * Therneau (2021) T. M. Therneau.  _A Package for Survival Analysis in R_ , 2021.  URL https://CRAN.R-project.org/package=survival.  R package version 3.2-11.
  * Uno et al. (2007) H. Uno, T. Cai, L. Tian, and L.-J. Wei.  Evaluating prediction rules for t-year survivors with censored regression models.  _Journal of the American Statistical Association_ , 102(478):527–537, 2007.
  * Van der Laan and Rose (2011) M. J. Van der Laan and S. Rose.  _Targeted learning: causal inference for observational and experimental data_.  Springer Science & Business Media, 2011.
  * Van Der Laan and Rubin (2006) M. J. Van Der Laan and D. Rubin.  Targeted maximum likelihood learning.  _The international journal of biostatistics_ , 2(1), 2006\.
  * Van der Laan et al. (2003) M. J. Van der Laan, M. Laan, and J. M. Robins.  _Unified methods for censored longitudinal data and causality_.  Springer Science & Business Media, 2003.
  * Wolbers et al. (2014) M. Wolbers, P. Blanche, M. T. Koller, J. C. Witteman, and T. A. Gerds.  Concordance for prognostic models with competing risks.  _Biostatistics_ , 15(3):526–539, 2014.
  * Yadlowsky et al. (2019) S. Yadlowsky, S. Basu, and L. Tian.  A calibration metric for risk scores with survival data.  In _Machine Learning for Healthcare Conference_ , pages 424–450, 2019\.
  * Yu et al. (2011) C.-N. Yu, R. Greiner, H.-C. Lin, and V. Baracos.  Learning patient-specific cancer survival distributions as a sequence of dependent regressors.  In _Advances in Neural Information Processing Systems_ , pages 1845–1853, 2011.
  * Zhong and Tibshirani (2019) C. Zhong and R. Tibshirani.  Survival analysis as a classification problem.  _arXiv preprint arXiv:1909.11171_ , 2019.

##  Appendix A Notation, Assumptions, and Likelihoods in More Detail

###  A.1 Notation

Let $T$ be a failure time with cdf $F$. $T$’s survival function is defined by
$\overline{F}=1-F$. We denote failure models by $F_{\theta_{T}}$. Let $C$ be a
censoring time with cdf $G$, survival function $\overline{G}$, and model
$G_{\theta_{C}}$. Under right-censoring, define $U=\min(T,C)$,
$\Delta=\mathbbm{1}\left[T\leq C\right]$ and we observe
$(X_{i},U_{i},\Delta_{i})$. We use $\overline{G}(t^{-})$ to denote $P(C\geq
t)$.

###  A.2 Assumptions

We assume i.i.d. data and random censoring: $T\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}C~{}|~{}X$ [Kalbfleisch and Prentice, 2002].
Derivations in this work also require the censoring positivity assumption
[Gerds et al., 2013]. Let $f=dF$ (a failure density) and $g=dG$ (a censoring
density). Then we assume

$\displaystyle\exists\epsilon\quad\text{ s.t. }\quad\forall x\,\forall
t\in\\{t\leq
t_{\text{max}}~{}|~{}f(t|x)>0\\},\quad\overline{G}(t^{-}|x)\geq\epsilon>0,$
(6)

for some truncation time $t_{\text{max}}$. Truncating at a maximum time is
necessary in practice for continuous distributions because datasets may have
no samples in the tails, leading to practical positivity violations [Gerds et
al., 2013]. This truncation happens implicitly for categorical models by
choosing the bins.

To observe censoring events properly, we also require a version of eq. 1 to
hold with the roles of $F$ and $G$ reversed:

$\displaystyle\exists\epsilon\quad\text{ s.t. }\quad\forall x\,\forall
t\in\\{t\leq
t_{\text{max}}~{}|~{}g(t|x)>0\\},\quad\overline{F}(t|x)\geq\epsilon>0.$ (7)

$t_{\text{max}}$ should be chosen so that these two conditions hold.

###  A.3 Likelihoods

As mentioned, we assume data are i.i.d. and censoring is random
$T\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}C~{}|~{}X$. Under these assumptions, the
likelihood, by definition [Andersen et al., 2012], is:

$\displaystyle L(\theta_{T},\theta_{C})$
$\displaystyle=\prod_{i}\Big{[}f_{\theta_{T}}(U_{i})\overline{G}_{\theta_{C}}(U_{i}^{-})\Big{]}^{\Delta_{i}}\Big{[}g_{\theta_{C}}(U_{i})\overline{F}_{\theta_{T}}(U_{i})\Big{]}^{1-\Delta_{i}},$
(8)

When a failure is observed, $\Delta_{i}=\mathbbm{1}\left[T_{i}\leq
C_{i}\right]=1$ so we compute the failure density or mass $f$ at the observed
time $U_{i}=T_{i}$. In this case, the only thing we know about the censoring
time is $C_{i}\geq T_{i}=U_{i}$. We therefore compute $P(C_{i}\geq
T_{i})=P(C_{i}\geq
U_{i})=1-G_{\theta_{C}}(U_{i}^{-})=\overline{G}_{\theta_{C}}(U_{i}^{-})$.
Likewise, when a censoring time is observed, $\Delta_{i}=0$ so we compute the
censoring density or mass $g$ at the observed censoring time $U_{i}=C_{i}$. In
this case, the only thing we know about the failure time is that
$T_{i}>C_{i}$. We therefore compute
$P(T_{i}>C_{i})=P(T_{i}>U_{i})=1-F(U_{i})=\overline{F}(U_{i})$.

Under the additional assumption of non-informativeness -that $F$ and $G$ don’t
share parameters and therefore $\theta_{T},\theta_{C}$ are distinct- the $g/G$
terms are constant wrt $\theta_{T}$ and the $f/F$ terms are constant wrt
$\theta_{C}$. In this case, when one is modeling failures, they can use the
partial failure likelihood:

$\displaystyle L(\theta_{T})^{\text{partial}}$
$\displaystyle=\prod_{i}\Big{[}f_{\theta_{T}}(U_{i})\Big{]}^{\Delta_{i}}\Big{[}\overline{F}_{\theta_{T}}(U_{i})\Big{]}^{1-\Delta_{i}}$

And when one is modeling censoring they can use the partial censoring
likelihood:

$\displaystyle L(\theta_{C})^{\text{partial}}$
$\displaystyle=\prod_{i}\Big{[}\overline{G}_{\theta_{C}}(U_{i}^{-})\Big{]}^{\Delta_{i}}\Big{[}g_{\theta_{C}}(U_{i})\Big{]}^{1-\Delta_{i}}$

###  A.4 Failure partial likelihood depends on true censoring distribution

We now show that the failure partial likelihood’s scale depends on the true
sampling distribution of censoring times, even if the censoring model has
dropped as a constant in the objective. The expected likelihood is:

$\displaystyle\mathop{\mathbb{E}}_{\begin{subarray}{c}T\sim
F_{\theta_{T}^{*}},C\sim G_{\theta_{C}^{*}}\\\
U=\min(T,C),\Delta=\mathbbm{1}\left[T\leq
C\right]\end{subarray}}\Big{[}f_{\theta_{T}}(U)^{\mathbbm{1}\left[\Delta=1\right]}\overline{F}_{\theta_{T}}(U)^{\mathbbm{1}\left[\Delta=0\right]}\Big{]}$

The reason is that $\Delta$ and $U$ depend on T and C (therefore on
$F_{\theta_{T}^{*}}$ and $G_{\theta_{C}^{*}}$). We now constructively show
that the failure model’s nll can vary with the true censoring distribution.
Let us consider a marginal survival analysis problem (no features) and random
censoring. The log nll is:

$\displaystyle\mathop{\mathbb{E}}_{F_{\theta_{T}^{*}},G_{\theta_{C}^{*}}}[\Delta\log
f_{\theta_{T}}(U)]+\mathop{\mathbb{E}}_{F_{\theta_{T}^{*}},G_{\theta_{C}^{*}}}[(1-\Delta)\log\overline{F}_{\theta_{T}}(U)]$

Now consider an $F_{\theta_{T}^{*}}$ whose support starts at time $1$ (e.g.
uniform over 1,2,3) and $G_{\theta_{C}^{*}}$ such that there is probability
$\rho$ that $C=0$ and probability $1-\rho$ that $C$ take a value above the
support of $T$ (e.g. >3). Points are therefore only censored at time 0 or
uncensored.

$\displaystyle\mathop{\mathbb{E}}_{F_{\theta_{T}^{*}},G_{\theta_{C}^{*}}}[\Delta\log
f_{\theta_{T}^{*}}(U)]+\mathop{\mathbb{E}}_{F_{\theta_{T}^{*}},G_{\theta_{C}^{*}}}[(1-\Delta)\log\overline{F}_{\theta_{T}^{*}}(U)]$
$\displaystyle=(1-\rho)\mathop{\mathbb{E}}_{F_{\theta_{T}^{*}}}[\log
f_{\theta_{T}^{*}}(T)]+\rho\mathop{\mathbb{E}}_{G_{\theta_{C}^{*}}}[\log\overline{F}_{\theta_{T}^{*}}(C)]$
$\displaystyle=(1-\rho)\mathop{\mathbb{E}}_{F_{\theta_{T}^{*}}}[\log
f_{\theta_{T}^{*}}(T)]+\rho\mathop{\mathbb{E}}_{G_{\theta_{C}^{*}}}[\log\overline{F}_{\theta_{T}^{*}}(0)]$
$\displaystyle=(1-\rho)\mathop{\mathbb{E}}_{F_{\theta_{T}^{*}}}[\log
f_{\theta_{T}^{*}}(T)]+\rho\mathop{\mathbb{E}}_{G_{\theta_{C}^{*}}}[\log 1]$
$\displaystyle=(1-\rho)\mathop{\mathbb{E}}_{F_{\theta_{T}^{*}}}[\log
f_{\theta_{T}^{*}}(T)]+\rho\mathop{\mathbb{E}}_{G_{\theta_{C}^{*}}}[0]$
$\displaystyle=(1-\rho)\mathop{\mathbb{E}}_{F_{\theta_{T}^{*}}}[\log
f_{\theta_{T}^{*}}(T)]$

This quantity depends on $\rho$. This shows that the failure model’s nll
depends on the true sampling distribution of censoring times.

##  Appendix B Ipcw Primer

Ipcw is a technique for estimation under censoring [Gerds and Schumacher,
2006]. Consider estimating the marginal mean of
$T:\mathop{\mathbb{E}}[T]=\mathop{\mathbb{E}}_{X}\mathop{\mathbb{E}}_{T|X}[T]$.
$T$ is not observed for all datapoints. Instead, we observe $U=\min(T,C)$ and
$\Delta=\mathbbm{1}\left[T\leq C\right]$. Ipcw reformulates such expectations
in terms of observed data. Using this method, we can show that:

$\displaystyle\begin{split}\mathop{\mathbb{E}}_{X}\mathop{\mathbb{E}}_{T|X}[T]&=\mathop{\mathbb{E}}_{X}\mathop{\mathbb{E}}_{T|X}\Bigg{[}\frac{\mathop{\mathbb{E}}_{C|X}\mathbbm{1}\left[T\leq
C\right]}{\mathop{\mathbb{E}}_{C^{\prime}|X}\mathbbm{1}\left[T\leq
C^{\prime}\right]}T\Bigg{]}\\\
&=\mathop{\mathbb{E}}_{X}\mathop{\mathbb{E}}_{T|X}\mathop{\mathbb{E}}_{C|X}\Bigg{[}\frac{\mathbbm{1}\left[T\leq
C\right]}{\mathop{\mathbb{E}}_{C^{\prime}|X}\mathbbm{1}\left[T\leq
C^{\prime}\right]}T\Bigg{]}\\\
&=\mathop{\mathbb{E}}_{T,C,X}\Bigg{[}\frac{\mathbbm{1}\left[T\leq
C\right]}{\mathop{\mathbb{E}}_{C^{\prime}|X}\mathbbm{1}\left[T\leq
C^{\prime}\right]}T\Bigg{]}\\\
&=\mathop{\mathbb{E}}_{T,C,X}\Bigg{[}\frac{\mathbbm{1}\left[T\leq
C\right]}{\mathbb{P}(C^{\prime}\geq T|X)}T\Bigg{]}\\\
&=\mathop{\mathbb{E}}_{T,C,X}\Bigg{[}\frac{\mathbbm{1}\left[T\leq
C\right]}{\overline{G}(T^{-}|X)}T\Bigg{]}\\\
&=\mathop{\mathbb{E}}_{U,\Delta,X}\Bigg{[}\frac{\Delta
U}{\overline{G}(U^{-}|X)}\Bigg{]}\end{split}$

We have used $C^{\prime}$ in the denominator to emphasize that it is not a
function of $C$ in the integral over the numerator indicator once that
expectation is moved out. We have used random censoring to go from
$\mathop{\mathbb{E}}_{T|X}\mathop{\mathbb{E}}_{C|X}$ to the joint
$\mathop{\mathbb{E}}_{T,C|X}$. The last equality changes from the complete
data distribution to the observed distribution and holds because
$\Delta=1\implies U=T$. This means we can estimate the expectation, provided
that we know $G$ and that random censoring and positivity (eq. 1) hold. In
practice, we must learn the censoring distribution, a challenging task as it
is also censored.

Graf et al. [1999] develop the ipcw bs. Gerds and Schumacher [2006] extend it
to conditional censoring and Kvamme and Borgan [2019a] specialize to
administrative censoring. Gerds et al. [2013], Wolbers et al. [2014] develop
the ipcw concordance. ipcw estimators for several forms of area under curve
(auc) have been studied in Hung and Chiang [2010a, b], Blanche et al. [2013,
2019], Uno et al. [2007]. Yadlowsky et al. [2019] derive an ipcw estimator for
binary survival calibration.

##  Appendix C Deriving ipcw Brier Scores

We derive the ipcw bs introduced by Graf et al. [1999], Gerds and Schumacher
[2006]. In the below let F-BS be the F model’s BS and let F-BS-CW be its
censor-weighted version. The censor-weighted failure bs:

$\displaystyle\text{F-BS-
CW}(t)=\mathop{\mathbb{E}}_{T,C}\Big{[}\frac{(1-F_{\theta}(t))^{2}\mathbbm{1}\left[T\leq
C\right]\mathbbm{1}\left[U\leq t\right]}{P_{\theta}(C^{\prime}\geq
U)}+\frac{F_{\theta}(t)^{2}\mathbbm{1}\left[U>t\right]}{P_{\theta}(C^{\prime}>t)}\Big{]}$

where $U=\min(T,C)$ and $F_{\theta}=P_{\theta}(T\leq\cdot)$, It’s relationship
to the regular bs is:

$\displaystyle\text{F-BS}(t)$
$\displaystyle=\mathop{\mathbb{E}}_{T}\Big{[}\Big{(}F_{\theta}(t)-\mathbbm{1}\left[T\leq
t\right]\Big{)}^{2}\Big{]}$
$\displaystyle=\mathop{\mathbb{E}}_{T}\Big{[}(1-F_{\theta}(t))^{2}\mathbbm{1}\left[T\leq
t\right]+F_{\theta}(t)^{2}\mathbbm{1}\left[T>t\right]\Big{]}$
$\displaystyle=\mathop{\mathbb{E}}_{T}\Big{[}\frac{\mathop{\mathbb{E}}_{C}\mathbbm{1}\left[T\leq
C\right]}{\mathop{\mathbb{E}}_{C^{\prime}}\mathbbm{1}\left[T\leq
C^{\prime}\right]}(1-F_{\theta}(t))^{2}\mathbbm{1}\left[T\leq
t\right]+\frac{\mathop{\mathbb{E}}_{C}\mathbbm{1}\left[C>t\right]}{\mathop{\mathbb{E}}_{C^{\prime}}\mathbbm{1}\left[C^{\prime}>t\right]}F_{\theta}(t)^{2}\mathbbm{1}\left[T>t\right]\Big{]}$
$\displaystyle=\mathop{\mathbb{E}}_{T,C}\Big{[}\frac{(1-F_{\theta}(t))^{2}\mathbbm{1}\left[T\leq
C\right]\mathbbm{1}\left[T\leq
t\right]}{\mathop{\mathbb{E}}_{C^{\prime}}\mathbbm{1}\left[T\leq
C^{\prime}\right]}+\frac{F_{\theta}(t)^{2}\mathbbm{1}\left[T>t\right]\mathbbm{1}\left[C>t\right]}{\mathop{\mathbb{E}}_{C^{\prime}}\mathbbm{1}\left[C^{\prime}>t\right]}\Big{]}$
$\displaystyle=\mathop{\mathbb{E}}_{T,C}\Big{[}\frac{(1-F_{\theta}(t))^{2}\mathbbm{1}\left[T\leq
C\right]\mathbbm{1}\left[T\leq t\right]}{P_{\theta}(C^{\prime}\geq
T)}+\frac{F_{\theta}(t)^{2}\mathbbm{1}\left[T>t\right]\mathbbm{1}\left[C>t\right]}{P_{\theta}(C^{\prime}>t)}\Big{]}$
$\displaystyle=\mathop{\mathbb{E}}_{T,C}\Big{[}\frac{(1-F_{\theta}(t))^{2}\mathbbm{1}\left[T\leq
C\right]\mathbbm{1}\left[U\leq t\right]}{P_{\theta}(C^{\prime}\geq
U)}+\frac{F_{\theta}(t)^{2}\mathbbm{1}\left[U>t\right]}{P_{\theta}(C^{\prime}>t)}\Big{]}$
$\displaystyle=\text{F-BS-CW}(t)$

The expectation comes out due to $T\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}C$. The last line follows from $T\leq
C\implies U=T$ (in the left term) and
$\mathbbm{1}\left[T>t\right]\mathbbm{1}\left[C>t\right]=\mathbbm{1}\left[U>t\right]$
(in the right term). Define likewise the failure-weighted censor bs

$\displaystyle\text{G-BS-
CW}(t)=\mathop{\mathbb{E}}_{T,C}\Big{[}\frac{(1-G_{\theta}(t))^{2}\mathbbm{1}\left[C<T\right]\mathbbm{1}\left[U\leq
t\right]}{P_{\theta}(T^{\prime}>U)}+\frac{G_{\theta}(t)^{2}\mathbbm{1}\left[U>t\right]}{P_{\theta}(T^{\prime}>t)}\Big{]}$

where $G_{\theta}=P_{\theta}(C\leq\cdot)$. The relationship to the censoring
distribution’s bs is:

$\displaystyle\text{G-BS}(t)$
$\displaystyle=\mathop{\mathbb{E}}_{C}\Big{[}\Big{(}G_{\theta}(t)-\mathbbm{1}\left[C\leq
t\right]\Big{)}^{2}\Big{]}$
$\displaystyle=\mathop{\mathbb{E}}_{C}\Big{[}(1-G_{\theta}(t))^{2}\mathbbm{1}\left[C\leq
t\right]+G_{\theta}(t)^{2}\mathbbm{1}\left[C>t\right]\Big{]}$
$\displaystyle=\mathop{\mathbb{E}}_{C}\Big{[}\frac{\mathop{\mathbb{E}}_{T}\mathbbm{1}\left[C<T\right]}{\mathop{\mathbb{E}}_{T^{\prime}}\mathbbm{1}\left[C<T^{\prime}\right]}(1-G_{\theta}(t))^{2}\mathbbm{1}\left[C\leq
t\right]+\frac{\mathop{\mathbb{E}}_{T}\mathbbm{1}\left[T>t\right]}{\mathop{\mathbb{E}}_{T^{\prime}}\mathbbm{1}\left[T^{\prime}>t\right]}G_{\theta}(t)^{2}\mathbbm{1}\left[C>t\right]\Big{]}$
$\displaystyle=\mathop{\mathbb{E}}_{T,C}\Big{[}\frac{(1-G_{\theta}(t))^{2}\mathbbm{1}\left[C<T\right]\mathbbm{1}\left[C\leq
t\right]}{\mathop{\mathbb{E}}_{T^{\prime}}\mathbbm{1}\left[C<T^{\prime}\right]}+\frac{G_{\theta}(t)^{2}\mathbbm{1}\left[T>t\right]\mathbbm{1}\left[C>t\right]}{\mathop{\mathbb{E}}_{T^{\prime}}\mathbbm{1}\left[T^{\prime}>t\right]}\Big{]}$
$\displaystyle=\mathop{\mathbb{E}}_{T,C}\Big{[}\frac{(1-G_{\theta}(t))^{2}\mathbbm{1}\left[C<T\right]\mathbbm{1}\left[C\leq
t\right]}{P_{\theta}(T^{\prime}>C)}+\frac{G_{\theta}(t)^{2}\mathbbm{1}\left[T>t\right]\mathbbm{1}\left[C>t\right]}{P_{\theta}(T^{\prime}>t)}\Big{]}$
$\displaystyle=\mathop{\mathbb{E}}_{T,C}\Big{[}\frac{(1-G_{\theta}(t))^{2}\mathbbm{1}\left[C<T\right]\mathbbm{1}\left[U\leq
t\right]}{P_{\theta}(T^{\prime}>U)}+\frac{G_{\theta}(t)^{2}\mathbbm{1}\left[U>t\right]}{P_{\theta}(T^{\prime}>t)}\Big{]}$
$\displaystyle=\text{G-BS-CW}(t)$

The expectation comes out due to $T\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}C$. The last line follows from $C<T\implies
U=C$ (in the left term) and
$\mathbbm{1}\left[T>t\right]\mathbbm{1}\left[C>t\right]=\mathbbm{1}\left[U>t\right]$
(in the right term).

##  Appendix D Negative Bernoulli Log Likelihood

Negative bll is similar to bs, but replaces the squared error with negated log
loss:

$\displaystyle\text{NBLL}(t;\theta)$
$\displaystyle=\mathop{\mathbb{E}}_{T,C,X}\Big{[}-\log({F_{\theta_{T}}}(t~{}|~{}X))\mathbbm{1}\left[T\leq
t\right]-\log(\overline{F}_{\theta_{T}}(t~{}|~{}X))\mathbbm{1}\left[T>t\right]\Big{]}$

Ipcw bll can likewise be written as [Kvamme et al., 2019]:

$\displaystyle\text{F-NBLL-CW}(t;\theta)$
$\displaystyle=\mathop{\mathbb{E}}_{T,C,X}\Big{[}\frac{-\log({F_{\theta_{T}}}(t~{}|~{}X))\Delta\mathbbm{1}\left[U\leq
t\right]}{G(U^{-}~{}|~{}X)}+\frac{-\log(\overline{F}_{\theta_{T}}(t~{}|~{}X))\mathbbm{1}\left[U>t\right]}{G(t~{}|~{}X)}\Big{]}$

##  Appendix E Game Algorithm

Input: Choice of losses $\ell_{F},\ell_{G}$, learning rate $\gamma$

Initialize $\theta_{Tt}$ and $\theta_{Ct}$ randomly for $t=1,\dots,K-1$

repeat

//  for each parameter of each player

for $t=1$ to $K-1$ do

$g_{Tt}\leftarrow d\ell_{F}^{t}/d\theta_{Tt}$

$g_{Ct}\leftarrow d\ell_{G}^{t}/d\theta_{Ct}$

end for

//  for each parameter of each player

for $t=1$ to $K-1$ do

$\theta_{Tt}\leftarrow\theta_{Tt}-\gamma g_{Tt}$

$\theta_{Ct}\leftarrow\theta_{Ct}-\gamma g_{Ct}$

end for

until convergence

Algorithm 2 Following Gradients in Multi-Player Games

##  Appendix F Experiments

###  F.1 Data

#### Gamma Simulation

We draw $x$ from a 32 dimensional multivariate normal $\mathcal{N}(0,10I)$. We
simulate conditionally gamma failure times with mean $\mu_{t}$ a log-linear
function of $x$ with coefficients for each feature drawn $\text{Unif}(0,0.1)$.
The censoring times are also conditionally gamma with mean $0.9*\mu_{t}$. Both
distributions have constant variance $0.05$. $\alpha,\beta$ parameterization
of the gamma is recovered from mean, variance by $\alpha=\mu^{2}/\sigma^{2}$
and $\beta=\mu/\sigma^{2}$. $T$ and $C$ are conditionally independent given
$X$. Each random seed draws a new dataset.

We report metrics as a function of training size. We use training sizes
[200,400,600,800,1000]. We use validation size 1024 and testing size 2048.

#### Survival mnist

Survival-mnist [Gensheimer, 2019, Pölsterl, 2019] draws times conditionally on
mnist label $Y$. This means digits define risk groups and
$T\mathchoice{\mathrel{\hbox to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X~{}|~{}Y$. Times within a digit are i.i.d.
The model only sees the image pixels $X$ as covariates so it must learn to
classify digits (risk groups) to model times. The PyCox package [Kvamme et
al., 2019] uses Exponential times. We follow Goldstein et al. [2020] and use
Gamma times. $T$’s mean is $10*(Y+1)$ so that lower labels $Y$ mean sooner
event times. We set the variance constant to $0.05$. $C$ is drawn similarly
but with $9.9*(Y+1)$. Each random seed draws a new dataset.

We report metrics as a function of training size. We use training sizes [512,
1024, 2048, 4096, 8192, 10240]. We use validation size 1024 and testing size
2048.

#### Real Data

We report results on

  * •

support [Knaus et al., 1995] which includes severely ill hospital patients.
There are 14 features. we split into 5,323 for training, 1774 for validation,
and 1776 for testing.

  * •

metabric [Curtis et al., 2012]. There are 9 features. We split into 1,142 for
training, 380 for validation, and 382 for testing.

  * •

rott [Foekens et al., 2000] and gbsg [Schumacher et al., 1994] combined into
one dataset (rott. & gbsg). There are 7 features. We split into 1,339 for
training, 446 for validation, and 447 for testing.

For more description see Therneau [2021], Katzman et al. [2018], Chen [2020].

In the main text, we report results on a subset of these datasets with metrics
as a function of training size. We use training sizes [10, 20, 30, 40, 50, 60,
70, 80, 90, 100, 110, 120, 130, 140, 150, 175, 200]. We use validation size
300 and always use the entire testing set. We standardize all real data with
the training set mean and standard deviation.

###  F.2 Models

In all experiments except for mnist, we use a 3-hidden-layer ReLU network. The
hidden sizes are [128, 64, 64] for the Gamma simulation and [128,256,64] for
the real data. We output $20$ categorical bins. See section G.1 for different
choices of number of bins, which did not show any significant differences in
results. For mnist we first use a small convolutional network and follow with
the same fully-connected network, but using hidden sizes [512,256,64].

###  F.3 Training

We use learning rate $0.001$ in all experiments for all losses using the Adam
optimizer. We train for 300 epochs for the simulated data and 200 for the real
data. For all data and all losses, this was enough to overfit on the training
data. We use no weight decay or dropout.

###  F.4 Model Selection

We select the best model on the validation set using the following approach:

  1. 1.

Save the $F$ and $G$ models from all the epochs in $F$-set and $G$-set.

  2. 2.

Randomly choose a model $\tilde{F}$ in the $F$-set.

  3. 3.

Use $\tilde{F}$ as the weight for $\ell_{G}$. Find the model $\tilde{G}$ from
$G$-set to minimize $\ell_{G}$ weighted by $\tilde{F}$.

  4. 4.

Use $\tilde{G}$ as the weight for $\ell_{F}$. Find the model $\tilde{F}$ from
$F$-set to minimize $\ell_{F}$ weighted by $\tilde{G}$.

  5. 5.

Repeat steps 3 and 4 until convergence.

Once converged, we use $\tilde{F}$ and $\tilde{G}$ as our best model to
evaluate at the test set. The above approach plays as similar role as the
game. Instead of gradient descent, this time we select a model from a set to
play the game. We first fix $F$ to find the best $G$ based on $\ell_{G}$ and
then fix $G$ to find the best $F$ based on $\ell_{F}$.

##  Appendix G Ablations

###  G.1 Changing number of bins on MNIST

Changing number of categorical bins (K) in [10,20,30,40,50]. Cannot directly
compare between two choices of K due to changing meaning of
likelihood/BS/Concordance but can compare nll and bs-Game at each K. Trends
similar across all choices of K.

(a) Uncensored bs (b) Uncensored Neg bll (c) Concordance (d) Categorical nll

Figure 9: 10 bins. nll (Blue). bs-Game (Orange).

(a) Uncensored bs (b) Uncensored Neg bll (c) Concordance (d) Categorical nll

Figure 10: 20 bins. nll (Blue). bs-Game (Orange).

(a) Uncensored bs (b) Uncensored Neg bll (c) Concordance (d) Categorical nll

Figure 11: 30 bins. nll (Blue). bs-Game (Orange).

(a) Uncensored bs (b) Uncensored Neg bll (c) Concordance (d) Categorical nll

Figure 12: 40 bins. nll (Blue). bs-Game (Orange).

(a) Uncensored bs (b) Uncensored Neg bll (c) Concordance (d) Categorical nll

Figure 13: 50 bins. nll (Blue). bs-Game (Orange).

##  Appendix H  Proof of Summed or Integrated Brier Score to be proper

######  Proposition 3.

Assume we have a list of time $t_{1},\dots,t_{K}$. Assume the true
distribution for $T$ is $F^{*}=F_{\theta_{T}^{*}}$ in eq. 2. We have:

  * •

The summed bs $\sum_{i=1}^{K}BS(t_{i};\theta)$ is proper, i.e., it has one
minimizer at the true parameters $\theta_{T}^{*}$.

  * •

The integrated bs $\int_{t_{1}}^{t_{K}}BS(t;\theta)dt$ is proper, i.e., it has
one minimizer at the true parameters $\theta_{T}^{*}$.

###### Proof.

Since $\lx@glossaries@gls@link{acronym}{bs}{\leavevmode\textsc{bs}}(t)$ is
proper, it has one minimizer at $\theta_{T}^{\star}$, i.e., for
$\theta_{T}\neq\theta_{T}^{*}$,
$\lx@glossaries@gls@link{acronym}{bs}{\leavevmode\textsc{bs}}(t;\theta_{T}^{\star})\leq\lx@glossaries@gls@link{acronym}{bs}{\leavevmode\textsc{bs}}(t;\theta_{T})$
for all $t$. Since this holds for all $t$, we then have:

$\sum_{i=1}^{K}BS(t_{i};\theta_{T}^{\star})\leq\sum_{i=1}^{K}BS(t_{i};\theta_{T}).$

This means that the summed Brier Score at $\theta_{T}^{*}$ is smaller than at
any other $\theta_{T}$. The summed bs has one minimizer at the true parameters
$\theta_{T}^{*}$, i.e., it is proper. Since the BS inequality holds for all
$t$, we also have

$\int_{t_{1}}^{t_{K}}BS(t;\theta_{T}^{\star})dt\leq\int_{t_{1}}^{t_{K}}BS(t;\theta_{T})dt$

This means that the integrated Brier Score at $\theta_{T}^{*}$ is smaller than
at any other $\theta_{T}$. The integrated bs has one minimizer at the true
parameters $\theta_{T}^{*}$, i.e., it is proper. ∎

##  Appendix I Proof of proposition 1

Here we prove that the true solution is a stationary point of the game. We
restate the proposition here.

######  Proposition.

Assume
$\exists\theta_{T}^{\star}\in\Theta_{T},\exists\theta_{C}^{\star}\in\Theta_{C}$
such that $F^{\star}=F_{\theta_{T}^{\star}}$ and
$G^{\star}=G_{\theta_{C}^{\star}}$. Assume the game losses $\ell_{F},\ell_{G}$
are based on proper losses $L$ and that the games are only computed at times
for which positivity holds. Then $(\theta_{T}^{\star},\theta_{C}^{\star})$ is
a stationary point of the game eq. 4.

$\displaystyle\begin{split}\ell_{F}(\theta)=L_{I}(F_{\theta_{T}};G_{\theta_{C}}),\quad\ell_{G}(\theta)&=L_{I}(G_{\theta_{C}};F_{\theta_{T}})\end{split}$
(4)

###### Proof.

In $\ell_{F}(\theta)$, by the definition of the ipcw estimator, when
$\theta_{C}=\theta_{C}^{*}$,
$L_{I}(F_{\theta_{T}};G_{\theta_{C}})=L(F_{\theta_{T}})$. Due to the fact that
$L$ is proper, $\theta_{T}^{*}$ is a minimizer for $L(F_{\theta_{T}})$. Then
at $(\theta_{T},\theta_{C})=(\theta_{T}^{*},\theta_{C}^{*})$, we have

$\left.\frac{d\ell_{F}(\theta)}{d\theta_{T}}\right|_{\begin{subarray}{c}\theta_{T}=\theta_{T}^{*}\\\
\theta_{C}=\theta_{C}^{*}\end{subarray}}=\left.\frac{dL_{I}(F_{\theta_{T}};G_{\theta_{C}^{*}})}{d\theta_{T}}\right|_{\theta_{T}=\theta_{T}^{*}}=\left.\frac{dL(F_{\theta_{T}})}{d\theta_{T}}\right|_{\theta_{T}=\theta_{T}^{*}}=0$

Similarly for $\ell_{G}(\theta)$, we have

$\left.\frac{d\ell_{G}(\theta)}{d\theta_{C}}\right|_{\begin{subarray}{c}\theta_{C}=\theta_{C}^{*}\\\
\theta_{T}=\theta_{T}^{*}\end{subarray}}=\left.\frac{dL_{I}(G_{\theta_{C}};F_{\theta_{T}^{*}})}{d\theta_{C}}\right|_{\theta_{C}=\theta_{C}^{*}}=\left.\frac{dL(G_{\theta_{C}})}{d\theta_{C}}\right|_{\theta_{C}=\theta_{C}^{*}}=0$

Since the two gradients are zero, the game will stay at the true parameters.
Therefore, $(\theta_{T}^{\star},\theta_{C}^{\star})$ is a stationary point of
the game eq. 4. ∎

##  Appendix J  Proof of proposition 2

Here we prove that under one construction of the game in algorithm 2, the true
solution is the unique stationary point of the game. We restate the
proposition here.

######  Proposition.

Consider discrete distributions over $K$ times. Let
$\theta_{T}=\\{\theta_{T1},\cdots,\theta_{T(K-1)}\\}$,
$\theta_{Tt}=P_{\theta}(T=t)$, $F_{\theta_{T}}(t)=\sum_{k=1}^{t}\theta_{Tk}$,
and likewise for $C,\theta_{C}$. Assuming that $\theta^{\star}_{Tt}>0$ and
$\theta^{\star}_{Ct}>0$, the solution
$(\theta_{T}^{\star},\theta_{C}^{\star})$ is the only stationary point for the
multi-player bs game shown in algorithm 2 for times $t\in\\{1,\ldots,K-1\\}$

###### Proof.

We show by induction on the time $t$ of the ipcw bs game that the simultaneous
gradient equations are only satisfied at $\hat{\theta_{T}}=\theta_{T}^{\star}$
and $\hat{\theta}_{C}=\theta_{C}^{\star}$. There is a lot of arithmetic but
eventually it comes down to (1) substitution of one variable for another (2)
assuming all previous timestep parameters are correct (induction) (3) finding
the zeros of a quadratic (4) showing that one of the two solutions is the
correct parameter and the other is invalid.

Note: this proof uses the notation that $\hat{\theta}$ is a model parameter
and $\theta^{\star}$ is the correct one.

###  J.1 BS(1) (base case)

We can compute the expectations defining $\text{F-BS-CW}(1)$ and $\text{G-BS-
CW}(1)$ in closed form. That gives us:

$\displaystyle\text{F-BS-CW}(1)$
$\displaystyle=\theta^{\star}_{T1}(1-\hat{\theta}_{T1})^{2}+(1-\theta^{\star}_{T1})(1-\theta^{\star}_{C1})\frac{\hat{\theta}_{T1}^{2}}{1-\hat{\theta}_{C1}}$
$\displaystyle\text{G-BS-CW}(1)$
$\displaystyle=\frac{\theta^{\star}_{C1}(1-\theta^{\star}_{T1})(1-\hat{\theta}_{C1})^{2}}{1-\hat{\theta}_{T1}}+(1-\theta^{\star}_{T1})(1-\theta^{\star}_{C1})\frac{\hat{\theta}_{C1}^{2}}{1-\hat{\theta}_{T1}}$

The derivatives are

$\displaystyle\frac{d\text{F-BS-CW}(1)}{d\hat{\theta}_{T1}}$
$\displaystyle=2\frac{(1-\theta^{\star}_{T1})(1-\theta^{\star}_{C1})}{1-\hat{\theta}_{C1}}\hat{\theta}_{T1}-2(1-\hat{\theta}_{T1})\theta^{\star}_{T1}=0$
$\displaystyle\frac{d\text{G-BS-CW}(1)}{d\hat{\theta}_{C1}}$
$\displaystyle=2\frac{(1-\theta^{\star}_{T1})(1-\theta^{\star}_{C1})}{1-\hat{\theta}_{T1}}\hat{\theta}_{C1}-2\frac{(1-\theta^{\star}_{T1})(1-\hat{\theta}_{C1})\theta^{\star}_{C1}}{1-\hat{\theta}_{T1}}=0$

We can take each derivative equation and write one variable in terms of the
other. First, taking $d\text{F-BS-CW}/d\hat{\theta}_{T1}$ and writing
$\hat{\theta}_{T1}$ in terms of $\hat{\theta}_{C1}$:

$\displaystyle\frac{d\text{F-BS-
CW}(1)}{d\hat{\theta}_{T1}}=2\frac{(1-\theta^{\star}_{T1})(1-\theta^{\star}_{C1})}{1-\hat{\theta}_{C1}}\hat{\theta}_{T1}-2(1-\hat{\theta}_{T1})\theta^{\star}_{T1}=0$

implies

$\displaystyle\frac{(1-\theta^{\star}_{T1})(1-\theta^{\star}_{C1})}{1-\hat{\theta}_{C1}}\hat{\theta}_{T1}=(1-\hat{\theta}_{T1})\theta^{\star}_{T1}$
$\displaystyle\frac{(1-\theta^{\star}_{T1})(1-\theta^{\star}_{C1})}{1-\hat{\theta}_{C1}}\hat{\theta}_{T1}+\theta^{\star}_{T1}\hat{\theta}_{T1}=\theta^{\star}_{T1}$
$\displaystyle\Big{(}\frac{(1-\theta^{\star}_{T1})(1-\theta^{\star}_{C1})}{1-\hat{\theta}_{C1}}+\theta^{\star}_{T1}\Big{)}\hat{\theta}_{T1}=\theta^{\star}_{T1}$
$\displaystyle\hat{\theta}_{T1}=\frac{\theta^{\star}_{T1}}{\Big{(}\frac{(1-\theta^{\star}_{T1})(1-\theta^{\star}_{C1})}{1-\hat{\theta}_{C1}}+\theta^{\star}_{T1}\Big{)}}$

Now solving for $\hat{\theta}_{C1}$ in the G-BS-CS derivative:

$\displaystyle\frac{d\text{G-BS-
CW}(1)}{d\hat{\theta}_{C1}}=2\frac{(1-\theta^{\star}_{T1})(1-\theta^{\star}_{C1})}{1-\hat{\theta}_{T1}}\hat{\theta}_{C1}-2\frac{(1-\theta^{\star}_{T1})(1-\hat{\theta}_{C1})\theta^{\star}_{C1}}{1-\hat{\theta}_{T1}}=0$

implies

$\displaystyle\frac{(1-\theta^{\star}_{T1})(1-\theta^{\star}_{C1})}{1-\hat{\theta}_{T1}}\hat{\theta}_{C1}=\frac{(1-\theta^{\star}_{T1})(1-\hat{\theta}_{C1})}{1-\hat{\theta}_{T1}}\theta^{\star}_{C1}$

Given $1-\theta^{\star}_{T1}\neq 0$ and $1-\hat{\theta}_{T1}\neq 0$, we have

$(1-\theta^{\star}_{C1})\hat{\theta}_{C1}=(1-\hat{\theta}_{C1})\theta^{\star}_{C1}$

which gives us $\hat{\theta}_{C1}=\theta^{\star}_{C1}$. Given
$1-\theta^{\star}_{T1}\neq 0$ and $1-\hat{\theta}_{T1}\neq 0$, the above
derivative equations jointly imply

$\displaystyle\hat{\theta}_{T1}=\Big{(}\theta^{\star}_{T1}\Big{)}\Bigg{(}\frac{(1-\theta^{\star}_{T1})(1-\theta^{\star}_{C1})}{1-\hat{\theta}_{C1}}+\theta^{\star}_{T1}\Bigg{)}^{-1},\quad\hat{\theta}_{C1}=\theta^{\star}_{C1}$

Substituting $\hat{\theta}_{C1}=\theta^{\star}_{C1}$ in the formula for
$\hat{\theta}_{T1}$ in terms of $\hat{\theta}_{C1}$, we have

$\hat{\theta}_{T1}=\Big{(}\theta^{\star}_{T1}\Big{)}\Bigg{(}\frac{(1-\theta^{\star}_{T1})(1-\theta^{\star}_{C1})}{1-\theta^{\star}_{C1}}+\theta^{\star}_{T1}\Bigg{)}^{-1}=\frac{\theta^{\star}_{T1}}{(1-\theta^{\star}_{T1})+\theta^{\star}_{T1}}=\theta^{\star}_{T1}$

Therefore, under the assumptions, for the BS(1) case, we have the only
stationary point at the two true 1st-timestep parameters:
$\hat{\theta}_{T1}=\theta_{T1}^{\star}$ and
$\hat{\theta}_{C1}=\theta_{C1}^{\star}$.

###  J.2 Induction step

We can proceed by induction over timesteps. Claim: given $P_{\theta}(T\leq
a)=P^{\star}(T\leq a)$ and $P_{\theta}(C\leq a)=P^{\star}(C\leq a),\quad
a=1,\dots,k$, the stationary point of the game BS(k+1) has to satisfy
$P_{\theta}(T=k+1)=P^{\star}(T=k+1)$ and $P_{\theta}(C=k+1)=P^{\star}(C=k+1)$
i.e. $\hat{\theta}_{T,k+1}=\theta^{\star}_{T,k+1}$ and
$\hat{\theta}_{C,k+1}=\theta^{\star}_{C,k+1}$.We first simplify F-BS-CW.

$\displaystyle\text{F-BS-
CW}(k+1)=\mathop{\mathbb{E}}_{T,C}\Big{[}\frac{(1-F_{\theta}(k+1))^{2}\mathbbm{1}\left[T\leq
C\right]\mathbbm{1}\left[U\leq k+1\right]}{P_{\theta}(C^{\prime}\geq
U)}+\frac{F_{\theta}(k+1)^{2}\mathbbm{1}\left[U>k+1\right]}{P_{\theta}(C^{\prime}>k+1)}\Big{]}$

We simplify each term of F-BS-CW separately. The left term of F-BS-CW is

$\displaystyle\mathop{\mathbb{E}}_{T,C}\frac{(1-F_{\theta}(k+1))^{2}\mathbbm{1}\left[T\leq
C\right]\mathbbm{1}\left[U\leq k+1\right]}{P_{\theta}(C^{\prime}\geq U)}$
$\displaystyle=$ $\displaystyle
P_{\theta}(T>k+1)^{2}\mathop{\mathbb{E}}_{T,C}\frac{\mathbbm{1}\left[T\leq
C\right]\mathbbm{1}\left[U\leq k+1\right]}{P_{\theta}(C^{\prime}\geq U)}$
$\displaystyle=$ $\displaystyle
P_{\theta}(T>k+1)^{2}\sum_{a=1}^{K}\sum_{b=1}^{K}P^{\star}(T=a)P^{\star}(C=b)\frac{\mathbbm{1}\left[a\leq
b\right]\mathbbm{1}\left[\min(a,b)\leq
k+1\right]}{P_{\theta}(C^{\prime}\geq\min(a,b))}\quad$
$\displaystyle\quad\quad\quad\Big{[}\text{condition }\mathbbm{1}\left[a\leq
b\right]\text{moves from indicator to sum limits and }\min(a,b)=a\Big{]}$
$\displaystyle=$ $\displaystyle
P_{\theta}(T>k+1)^{2}\sum_{a=1}^{K}\sum_{b=a}^{K}\frac{P^{\star}(T=a)P^{\star}(C=b)\mathbbm{1}\left[a\leq
k+1\right]}{P_{\theta}(C^{\prime}\geq a)}$
$\displaystyle\quad\quad\quad\Big{[}\text{condition }\mathbbm{1}\left[a\leq
k+1\right]\text{ moves from indicator to sum limit}\Big{]}$ $\displaystyle=$
$\displaystyle
P_{\theta}(T>k+1)^{2}\sum_{a=1}^{k+1}\sum_{b=a}^{K}\frac{P^{\star}(T=a)P^{\star}(C=b)}{P_{\theta}(C^{\prime}\geq
a)}$ $\displaystyle=$ $\displaystyle
P_{\theta}(T>k+1)^{2}\sum_{a=1}^{k+1}P^{\star}(T=a)\sum_{b=a}^{K}\frac{P^{\star}(C=b)}{P_{\theta}(C^{\prime}\geq
a)}$ $\displaystyle=$ $\displaystyle
P_{\theta}(T>k+1)^{2}\sum_{a=1}^{k+1}P^{\star}(T=a)\frac{P^{\star}(C\geq
a)}{P_{\theta}(C^{\prime}\geq a)}$ $\displaystyle\Bigg{[}\text{induction
hypothesis: }P_{\theta}(C\leq a)=P^{\star}(C\leq a),\quad a=1,\dots,k\implies
P_{\theta}(C>a)=P^{\star}(C>a),\quad a=1,\dots,k\Bigg{]}$ $\displaystyle=$
$\displaystyle P_{\theta}(T>k+1)^{2}\sum_{a=1}^{k+1}P^{\star}(T=a)\cdot 1$
$\displaystyle=$ $\displaystyle P_{\theta}(T>k+1)^{2}P^{\star}(T\leq k+1)$
$\displaystyle=$
$\displaystyle(1-\sum_{i=1}^{k}\hat{\theta}_{Ti}-\hat{\theta}_{T(k+1)})^{2}\sum_{i=1}^{k+1}\theta^{\star}_{Ti}$
$\displaystyle\quad\quad\quad\Big{[}\text{induction hypothesis:
}P_{\theta}(T\leq a)=P^{\star}(T\leq a),\quad a=1,\dots,k\Big{]}$
$\displaystyle=$
$\displaystyle(1-\sum_{i=1}^{k}\theta^{\star}_{Ti}-\hat{\theta}_{T(k+1)})^{2}\sum_{i=1}^{k+1}\theta^{\star}_{Ti}$
$\displaystyle=$ $\displaystyle(1-p-x)^{2}(p+t)$
$\displaystyle\overset{\Delta}{=}$ $\displaystyle A,\quad\quad\text{where
}p=\sum_{i=1}^{k}\theta^{\star}_{Ti},q=\sum_{i=1}^{k}\theta^{\star}_{Ci},x=\hat{\theta}_{T(k+1)},y=\hat{\theta}_{C(k+1)}.,t=\theta^{\star}_{T(k+1)}c=\theta^{\star}_{C(k+1)}.$

The right term of F-BS-CW is

$\displaystyle\mathop{\mathbb{E}}_{T,C}\frac{F_{\theta}(k+1)^{2}\mathbbm{1}\left[U>k+1\right]}{P_{\theta}(C^{\prime}>k+1)}$
$\displaystyle=\frac{F_{\theta}(k+1)^{2}}{P_{\theta}(C^{\prime}>k+1)}\mathop{\mathbb{E}}_{T,C}\mathbbm{1}\left[U>k+1\right]$
$\displaystyle\quad\quad\quad\Bigg{[}T~{}\text{and}~{}C~{}\text{are
independent means
}\mathbbm{1}\left[U>z\right]=\mathbbm{1}\left[T>z\right]\mathbbm{1}\left[C>z\right]\Bigg{]}$
$\displaystyle=\frac{F_{\theta}(k+1)^{2}}{P_{\theta}(C^{\prime}>k+1)}P^{\star}(T>k+1)P^{\star}(C>k+1)$
$\displaystyle=\frac{(\sum_{i=1}^{k+1}\hat{\theta}_{Ti})^{2}}{1-\sum_{i=1}^{k+1}\hat{\theta}_{Ci}}(1-\sum_{i=1}^{k+1}\theta^{\star}_{Ti})(1-\sum_{i=1}^{k+1}\theta^{\star}_{Ci})$
$\displaystyle\quad\Bigg{[}\text{induction hypothesis: }P_{\theta}(T\leq
a)=P^{\star}(T\leq a)\text{ and }P_{\theta}(C\leq a)=P^{\star}(C\leq a),\quad
a=1,\dots,k\Bigg{]}$
$\displaystyle=\frac{(\sum_{i=1}^{k}\theta^{\star}_{Ti}+\hat{\theta}_{T(k+1)})^{2}}{1-\sum_{i=1}^{k}\theta^{\star}_{Ci}-\hat{\theta}_{C(k+1)}}(1-\sum_{i=1}^{k+1}\theta^{\star}_{Ti})(1-\sum_{i=1}^{k+1}\theta^{\star}_{Ci})$
$\displaystyle=\frac{(p+x)^{2}}{1-q-y}(1-p-t)(1-q-c)\triangleq B$

where again
$p=\sum_{i=1}^{k}\theta^{\star}_{Ti},q=\sum_{i=1}^{k}\theta^{\star}_{Ci},x=\hat{\theta}_{T(k+1)},y=\hat{\theta}_{C(k+1)},t=\theta^{\star}_{T(k+1)},c=\theta^{\star}_{C(k+1)}$.
To summarize, F-BS-CW$(k+1)=A+B$:

$\displaystyle\text{F-BS-
CW}(k+1)=(1-p-x)^{2}(p+t)+\frac{(p+x)^{2}}{1-q-y}(1-p-t)(1-q-c)$

Then we simplify G-BS-CW.

$\displaystyle\text{G-BS-
CW}(k+1)=\mathop{\mathbb{E}}_{T,C}\Big{[}\frac{(1-G_{\theta}(k+1))^{2}\mathbbm{1}\left[C<T\right]\mathbbm{1}\left[U\leq
k+1\right]}{P_{\theta}(T^{\prime}>U)}+\frac{G_{\theta}(k+1)^{2}\mathbbm{1}\left[U>k+1\right]}{P_{\theta}(T^{\prime}>k+1)}\Big{]}$

The left term of G-BS-CW

$\displaystyle\mathop{\mathbb{E}}_{T,C}\frac{(1-G_{\theta}(k+1))^{2}\mathbbm{1}\left[C<T\right]\mathbbm{1}\left[U\leq
k+1\right]}{P_{\theta}(T^{\prime}>U)}$ $\displaystyle=$
$\displaystyle(1-G_{\theta}(k+1))^{2}\mathop{\mathbb{E}}_{T,C}\frac{\mathbbm{1}\left[C<T\right]\mathbbm{1}\left[U\leq
k+1\right]}{P_{\theta}(T^{\prime}>U)}$ $\displaystyle=$
$\displaystyle(1-G_{\theta}(k+1))^{2}\sum_{a=1}^{K}\sum_{b=1}^{K}P^{\star}(C=a)P^{\star}(T=b)\frac{\mathbbm{1}\left[a<b\right]\mathbbm{1}\left[\min(a,b)\leq
k+1\right]}{P_{\theta}(T^{\prime}>\min(a,b))}$
$\displaystyle\quad\quad\quad\text{condition
}\mathbbm{1}\left[a<b\right]\text{moves from indicator to sum limits
and}\min(a,b)=a$ $\displaystyle=$
$\displaystyle(1-G_{\theta}(k+1))^{2}\sum_{a=1}^{K}\sum_{b=a+1}^{K}\frac{P^{\star}(C=a)P^{\star}(T=b)\mathbbm{1}\left[a\leq
k+1\right]}{P_{\theta}(T^{\prime}>a)}$
$\displaystyle\quad\quad\quad\text{condition }\mathbbm{1}\left[a\leq
k+1\right]\text{ moves from indicator to sum limits}$ $\displaystyle=$
$\displaystyle(1-G_{\theta}(k+1))^{2}\sum_{a=1}^{k+1}\sum_{b=a+1}^{K}\frac{P^{\star}(C=a)P^{\star}(T=b)}{P_{\theta}(T^{\prime}>a)}$
$\displaystyle\quad\quad\quad\Bigg{[}\text{split sum over a into two terms: 1
through k, and k+1, recall b starts at a+1}\Bigg{]}$ $\displaystyle=$
$\displaystyle(1-G_{\theta}(k+1))^{2}\left(\sum_{a=1}^{k}\sum_{b=a+1}^{K}\frac{P^{\star}(C=a)P^{\star}(T=b)}{P_{\theta}(T^{\prime}>a)}+\sum_{b=k+2}^{K}\frac{P^{\star}(C=k+1)P^{\star}(T=b)}{P_{\theta}(T^{\prime}>k+1)}\right)$
$\displaystyle=$
$\displaystyle(1-G_{\theta}(k+1))^{2}\left(\sum_{a=1}^{k}P^{\star}(C=a)\sum_{b=a+1}^{K}\frac{P^{\star}(T=b)}{P_{\theta}(T^{\prime}>a)}+P^{\star}(C=k+1)\sum_{b=k+2}^{K}\frac{P^{\star}(T=b)}{P_{\theta}(T^{\prime}>k+1)}\right)$
$\displaystyle=$
$\displaystyle(1-G_{\theta}(k+1))^{2}\left(\sum_{a=1}^{k}\frac{P^{\star}(C=a)P^{\star}(T\geq
a+1)}{P_{\theta}(T^{\prime}>a)}+\frac{P^{\star}(C=k+1)P^{\star}(T>k+1)}{P_{\theta}(T^{\prime}>k+1)}\right)$
$\displaystyle=$
$\displaystyle(1-G_{\theta}(k+1))^{2}\left(\sum_{a=1}^{k}\frac{P^{\star}(C=a)P^{\star}(T>a)}{P_{\theta}(T^{\prime}>a)}+\frac{P^{\star}(C=k+1)P^{\star}(T>k+1)}{P_{\theta}(T^{\prime}>k+1)}\right)$
$\displaystyle\Bigg{[}\text{induction hypothesis: }P_{\theta}(T\leq
a)=P^{\star}(T\leq a),\quad a=1,\dots,k\implies
P_{\theta}(T>a)=P^{\star}(T>a),\quad a=1,\dots,k\Bigg{]}$ $\displaystyle=$
$\displaystyle(1-G_{\theta}(k+1))^{2}\left(\sum_{a=1}^{k}P^{\star}(C=a)+\frac{P^{\star}(C=k+1)P^{\star}(T>k+1)}{P_{\theta}(T^{\prime}>k+1)}\right)$
$\displaystyle=$
$\displaystyle(1-\sum_{i=1}^{k}\hat{\theta}_{Ci}-\hat{\theta}_{C(k+1)})^{2}\left(\sum_{i=1}^{k}\theta^{\star}_{Ci}+\frac{\theta^{\star}_{C(k+1)}(1-\theta^{\star}_{T(k+1)}-\sum_{i=1}^{k}\theta^{\star}_{Ti})}{1-\sum_{i=1}^{k}\hat{\theta}_{Ti}-\hat{\theta}_{T(k+1)}}\right)$
$\displaystyle\quad\Bigg{[}\text{induction hypothesis: }P_{\theta}(T\leq
a)=P^{\star}(T\leq a)\quad\text{ and }\quad P_{\theta}(C\leq
a)=P^{\star}(C\leq a),\quad a=1,\dots,k\Bigg{]}$ $\displaystyle=$
$\displaystyle(1-\sum_{i=1}^{k}\theta^{\star}_{Ci}-\hat{\theta}_{C(k+1)})^{2}\left(\sum_{i=1}^{k}\theta^{\star}_{Ci}+\frac{\theta^{\star}_{C(k+1)}(1-\theta^{\star}_{T(k+1)}-\sum_{i=1}^{k}\theta^{\star}_{Ti})}{1-\sum_{i=1}^{k}\theta^{\star}_{Ti}-\hat{\theta}_{T(k+1)}}\right)$
$\displaystyle=$
$\displaystyle(1-q-y)^{2}(q+\frac{c(1-t-p)}{1-p-x})\overset{\Delta}{=}C$

By symmetry with $B$, the right term is

$\displaystyle\mathop{\mathbb{E}}_{T,C}\frac{G_{\theta}(k+1)^{2}\mathbbm{1}\left[U>k+1\right]}{P_{\theta}(T^{\prime}>k+1)}=\frac{(q+y)^{2}}{1-p-x}(1-q-c)(1-p-t)\overset{\Delta}{=}D$

Again using
$p=\sum_{i=1}^{k}\theta^{\star}_{Ti},q=\sum_{i=1}^{k}\theta^{\star}_{Ci},x=\hat{\theta}_{T(k+1)},y=\hat{\theta}_{C(k+1)},t=\theta^{\star}_{T(k+1)},c=\theta^{\star}_{C(k+1)}$,
we have

$\displaystyle\text{G-BS-CW}(k+1)$ $\displaystyle=C+D$
$\displaystyle=(1-q-y)^{2}(q+\frac{c(1-t-p)}{1-p-x})+\frac{(q+y)^{2}}{1-p-x}(1-q-c)(1-p-t)$

The stationary point satisfies

$\displaystyle\frac{\partial\text{G-wt-FBS(k+1)}}{\partial x}$
$\displaystyle=\frac{\partial A}{\partial x}+\frac{\partial B}{\partial x}$
$\displaystyle=-2(1-p-x)(p+t)+2\frac{(p+x)}{1-q-y}(1-p-t)(1-q-c)$
$\displaystyle=0$ $\displaystyle\frac{\partial\text{F-wt-GBS(k+1)}}{\partial
y}$ $\displaystyle=\frac{\partial C}{\partial y}+\frac{\partial D}{\partial
y}$
$\displaystyle=-2(1-q-y)(q+\frac{c(1-t-p)}{1-p-x})+2\frac{(q+y)}{1-p-x}(1-q-c)(1-p-t)=0$

It’s a system of quadratic equations with two unknowns. The system has
analytical solutions. Solving the above equations for $x,y$ by Mathematica (it
is quite a long derivation manually), the solutions are

$x=t,y=c$

or

$\displaystyle x=$ $\displaystyle(1/(-q+q^{2}+qc))(cp-qcp-qt+q^{2}t+ct$
$\displaystyle-(p(-1+q+c+qp-q^{2}p-cp+qt-q^{2}t-ct))/((-1+q)(p+t))$
$\displaystyle+(qp(-1+q+c+qp-q^{2}p-cp+qt-q^{2}t-ct))/((-1+q)(p+t))$
$\displaystyle-(t(-1+q+c+qp-q^{2}p-cp+qt-q^{2}t-ct))/((-1+q)(p+t))$
$\displaystyle+(qt(-1+q+c+qp-q^{2}p-cp+qt-q^{2}t-ct))/((-1+q)(p+t)))$
$\displaystyle y=$
$\displaystyle(-1+q+c+qp-q^{2}p-cp+qt-q^{2}t-ct)/((-1+q)(p+t))$

To check if this second solution is valid, it would need to be the case that
$q+y<1$ because we only consider $k+1<K$. If we ask mathematica to simplify
q+y that satisfies the above solution, then this holds:

$\displaystyle q+y=\frac{-1+q-c(-1+p+t)}{(-1+q)(p+t)}$

The numerator and the denominator are both negative. If $k+1<K$ (we know bs at
K is 0 and also we only have K-1 parameters), the numerator minus denominator
=

$\displaystyle-1+q-c(-1+p+t)-(-1+q)(p+t)$
$\displaystyle=(-1+q)(1-p-t)-c(-1+p+t)$ $\displaystyle=(-1+q+c)(1-p-t)$
$\displaystyle<0$

Therefore,

$\sum_{i=1}^{k}\theta^{\star}_{Ci}+\hat{\theta}_{C(k+1)}=q+y>1$

This is invalid. So

$x=t,y=c$

is the only solution, i.e.,
$\hat{\theta}_{T(k+1)}=\theta^{\star}_{T(k+1)},\hat{\theta}_{C(k+1)}=\theta^{\star}_{C(k+1)}$.
By induction, we conclude that

$\hat{\theta}_{Ti}=\theta^{\star}_{Ti},\hat{\theta}_{Ci}=\theta^{\star}_{Ci},i=1,\dots,K-1$

By $\hat{\theta}_{TK}=1-\sum_{i=1}^{K-1}\hat{\theta}_{Ti}$ and
$\hat{\theta}_{CK}=1-\sum_{i=1}^{K-1}\hat{\theta}_{Ci}$, we have

$\hat{\theta}_{TK}=\theta^{\star}_{TK},\hat{\theta}_{CK}=\theta^{\star}_{CK}$

Therefore,

$\hat{\theta}_{Ti}=\theta^{\star}_{Ti},\hat{\theta}_{Ci}=\theta^{\star}_{Ci},i=1,\dots,K$

is the only stationary point for the game. ∎