# The Brier Score under Administrative Censoring:
Problems and Solutions

Håvard Kvamme<EMAIL_ADDRESS>
Ørnulf Borgan<EMAIL_ADDRESS>
Department of Mathematics
University of Oslo
P.O. Box 1053 Blindern
0316 Oslo, Norway

###### Abstract

The Brier score is commonly used for evaluating probability predictions. In
survival analysis, with right-censored observations of the event times, this
score can be weighted by the inverse probability of censoring (IPCW) to retain
its original interpretation. It is common practice to estimate the censoring
distribution with the Kaplan-Meier estimator, even though it assumes that the
censoring distribution is independent of the covariates. This paper discusses
the general impact of the censoring estimates on the Brier score and shows
that the estimation of the censoring distribution can be problematic. In
particular, when the censoring times can be identified from the covariates,
the IPCW score is no longer valid. For administratively censored data, where
the potential censoring times are known for all individuals, we propose an
alternative version of the Brier score. This administrative Brier score does
not require estimation of the censoring distribution and is valid even if the
censoring times can be identified from the covariates.

Keywords: survival analysis, time-to-event-prediction, customer churn, inverse
probability weighting, progressive type I censoring

##  1 Introduction

Recently, there has been an increasing interest in combining machine learning
methodology with survival analysis for improved time-to-event prediction. Some
methods extend the well known Cox regression with neural networks (Katzman et
al., 2018; Ching et al., 2018; Yousefi et al., 2017; Luck et al., 2017; Kvamme
et al., 2019), while others consider a more direct approach for optimizing the
likelihood for right-censored time-to-event data (Biganzoli et al., 1998; Lee
et al., 2018; Fotso, 2018; Gensheimer and Narasimhan, 2019; Kvamme and Borgan,
2019). Also worth mentioning is the Random Survival Forest (Ishwaran et al.,
2008) which makes decision trees based on the log-rank test and estimates the
cumulative hazards with the Nelson-Aalen estimator.

Although these methods are available for right-censored event times, a
substantial part of the machine learning community is not familiar with
survival analysis and might find it reasonable to instead apply binary
classifiers for time-to-event prediction. In short, a binary classifier
estimates the probability that an individual experience the event by time $t$,
and can be fitted by disregarding individuals censored before that time. As
empirical evaluation of predictions is central to machine learning, the best
way to convince this audience to use survival methodology would likely be to
show that the survival methods give improved evaluation scores.

Arguably, the two most common evaluation criteria for survival predictions are
the inverse probability of censoring weighted (IPCW) Brier score (Graf et al.,
1999; Gerds and Schumacher, 2006) and different versions of the concordance
index (Harrell Jr et al., 1982; Antolini et al., 2005; Uno et al., 2011; Gerds
et al., 2013). We will, in this paper, show that the IPCW Brier score can be
biased under administrative censoring, i.e., when the right-censoring times
are known for _all_ individuals. In fact, we will show that this bias benefits
the binary classifiers, meaning that the binary classifiers can get better
scores than corresponding survival methods, even though the estimates of the
binary classifiers are, in reality, much worse than those of the survival
methods. Furthermore, we will show that the IPCW Brier score is heavily
dependent on the estimates of the censoring distribution, making it
unattractive as an evaluation metrics as one can question the validity of the
results. We, therefore, propose the _administrative Brier score_ , which is
created for handling both these issues. In particular, it does not have the
potential bias of the IPCW score under administrative censoring, and it does
not require estimation of the censoring distribution.

To give the reader an understanding of the potential problems of the IPCW
Brier score, we will start with an illustrative example in Section 2. Then, in
Section 3, we will introduce the Brier score in detail and discuss more
carefully the issues of the IPCW scheme. We present our proposed alternative,
the administrative Brier score, in Section 3.2. The binary classifiers are
investigated in Section 4, where we will show their relationship to the
potential bias of the IPCW Brier score under administrative censoring. A
simulation study is conducted in Section 5 to empirically illustrate our
findings, and in Section 6, we investigate a real data set with administrative
censoring. A summary and concluding remarks are made in Section 7.

The code for the evaluation metrics, the survival methods, the simulations,
and the data sets is available at github.com/havakv/pycox.

##  2 A Real-World Example

To illustrate the potential issues with the IPCW Brier score, we consider an
example encountered while researching the KKBox Churn data set in Kvamme et
al. (2019). The task is to predict whether or not customers continue to
subscribe to the KKBox music streaming service $t$ days after their first
subscription. If customers leave their subscription they have _churned_ , and
these are the events we want to model. The follow-up times of the customers
depend on when they first subscribed, as we only follow customers up to
January 29, 2017. A substantial part of the customers do not experience the
event before this date and are instead _right-censored_. The type of right-
censoring encountered here is called _administrative censoring_ or
_progressive type I censoring_ , meaning that we actually know the censoring
times of _all_ the customers. So even for customers for whom we have observed
the churn time, we know the time they would have been censored.

We approach the modeling of the event-time distribution in two ways. The first
is with the Logistic-Hazard method (Brown, 1975; Biganzoli et al., 1998;
Gensheimer and Narasimhan, 2019; Kvamme and Borgan, 2019) which accounts for
censored observations by considering the likelihood for right-censored event
times. We use the version of the method described by Kvamme and Borgan (2019),
meaning that we parameterize the discrete hazards with a neural network in the
form of a multilayer perceptron (MLP).

The second approach is to fit a binary classifier for each time $t$ and remove
all customers censored before this time. The responses (labels) given to the
classifiers are indicators of whether or not each customer has churned. We
denote this as the _BCE_ method, as it minimizes the binary cross-entropy of
the survival estimates, where survival means that a customer has not yet
churned. The BCE method is an MLP with equivalent network structure to that of
the Logistic-Hazard, with each output node corresponding to a binary
classifier at time $t$. The method will be described in detail in Section 4.
As censored individuals are removed from the data set at their time of
censoring, the BCE method has a bias towards higher churn probabilities. We
would, therefore, expect the Logistic-Hazard to perform better for the KKBox
data set.

Figure 1: Brier scores of survival estimates from the Logistic-Hazard method
and the BCE method on the KKBox data set. The Brier scores are weighted by the
inverse probability of censoring estimated with Kaplan-Meier (left) and with a
Logistic-Hazard model (right).

The IPCW Brier scores with Kaplan-Meier censoring estimates (Graf et al.,
1999) are calculated for the survival estimates of the two methods, and a plot
of the results is shown in the left part of Figure 1. For higher times, we see
that the BCE method has better scores than the Logistic-Hazard method. This is
unexpected, as the censoring proportion increases with time, meaning the bias
of the BCE estimates increases with time. However, the Kaplan-Meier censoring
estimates are not covariate dependent, so the results are likely explained by
poor estimation of the censoring distribution. This was addressed by Gerds and
Schumacher (2006), but the Kaplan-Meier estimator continues to be the most
common choice for estimating the weights in the IPCW Brier score.

In the right part of Figure 1, we have plotted the IPCW Brier scores where the
censoring distribution is estimated with a Logistic-Hazard model similar to
that of the survival estimates. We now see that the BCE generally performs
worse than the Logistic-Hazard, but we still find that it gets better scores
than the Logistic-Hazard for the largest times. We do, however, believe that
the BCE actually performs worse than the Logistic-Hazard, also for the largest
times, and that the reason for these results can be explained by some
artifacts of the IPCW Brier scores applied to a data set with administrative
censoring.

In the following sections, we will present what we believe to be the cause of
these results, and we will use simulated data to reproduce the artifacts. We
will, also, propose the _administrative Brier score_ , which does not suffer
from the same vulnerabilities as the IPCW Brier score, and does not even
require estimation of the censoring distribution. In Section 6, we will
revisit the KKBox data set for a more in-depth analysis.

##  3 The Brier Score

In the following, we present the Brier score for evaluating time-to-event
predictions in the form of survival estimates. We will then introduce the
topic of right-censoring in survival analysis, and show how the IPCW Brier
score accounts for censored observations. This is followed by a discussion of
administratively censored observations and how they affect the IPCW Brier
score. We end the section by introducing a new Brier score for handling
administratively right-censored observations.

Let $T^{*}_{i}$ denote the event time of individual $i$ and let $f_{i}(t)$
denote the density function of this event time. The survival function of
individual $i$ is then given by

$\displaystyle
S_{i}(t)=\textnormal{P}(T^{*}_{i}>t)=\int_{t}^{\infty}f_{i}(u)\,du.$

In time-to-event prediction, we want to estimate $S_{i}(t)$ for all
individuals, and we will let $\pi_{i}(t)$ denote these estimates. To simplify
the presentation, we will disregard the estimation uncertainty, and consider
the $\pi_{i}(t)$’s as known (non-random) functions.

A reasonable metric for evaluating the predictive performance of the
$\pi_{i}(t)$’s, is to calculate the mean squared error of the estimates to the
true survival functions. For a data set with $n$ individuals, this score is

$\displaystyle\text{MSE}(t)=\frac{1}{n}\sum^{n}_{i=1}{\left[S_{i}(t)-\pi_{i}(t)\right]}^{2}.$
(1)

However, $S_{i}(t)$ is not known outside of simulations, and we instead
observe event times $T^{*}_{i}$ drawn from the event time distribution. The
Brier score for uncensored data approximates the true survival functions with
step-functions with jumps at the event times, giving

$\displaystyle\text{BS}(t)$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}{\left[\mathbbm{1}\\{T^{*}_{i}>t\\}-\pi_{i}(t)\right]}^{2}$
(2)
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\left[{\pi_{i}(t)}^{2}\mathbbm{1}\\{T_{i}^{*}\leq
t\\}+{[1-\pi_{i}(t)]}^{2}\mathbbm{1}\\{T_{i}^{*}>t\\}\right]$

The expectation of the Brier score is

$\displaystyle\mathbb{E}\left[\text{BS}(t)\right]$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\mathbb{E}\left[{\pi_{i}(t)}^{2}\mathbbm{1}\\{T_{i}^{*}\leq
t\\}+{[1-\pi_{i}(t)]}^{2}\mathbbm{1}\\{T_{i}^{*}>t\\}\right]$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\left[{\pi_{i}(t)}^{2}\,\textnormal{P}(T_{i}^{*}\leq
t)+{[1-\pi_{i}(t)]}^{2}\,\textnormal{P}(T_{i}^{*}>t)\right]$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\left[{\pi_{i}(t)}^{2}\,[1-S_{i}(t)]+{[1-\pi_{i}(t)]}^{2}\,S_{i}(t)\right]$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\left({\left[S_{i}(t)-\pi_{i}(t)\right]}^{2}+S_{i}(t)\left[1-S_{i}(t)\right]\right)$
$\displaystyle=\text{MSE}(t)+\frac{1}{n}\sum^{n}_{i=1}S_{i}(t)\left[1-S_{i}(t)\right].$

So the expected Brier score is the sum of the MSE in (1) and a constant that
does not depend on our survival estimates $\pi_{i}(t)$. This constant is the
irreducible error of approximating the true survival functions $S_{i}(t)$ with
the step-functions $\mathbbm{1}\\{T^{*}_{i}>t\\}$. So minimizing the expected
Brier score is equivalent to minimizing the MSE, and the minimum is obtained
for the true survival functions, i.e., $\pi_{i}(t)=S_{i}(t)$.

###  3.1 The IPCW Brier Score

In many applications, only a subset of the event times $T^{*}_{i}$ is
observed. For some individuals, we instead only know that the event time
occurs after some censoring time $C^{*}_{i}$. For modeling these data sets, we
consider the right-censored event time $T_{i}=\min\\{T^{*}_{i},\,C^{*}_{i}\\}$
and the event indicator $D_{i}=\mathbbm{1}\\{T^{*}_{i}\leq C^{*}_{i}\\}$. We
follow the common convention in survival analysis that when an event and
censoring time coincide, we observe the occurrence of the event.

As we only have partial information, the Brier score in (2) cannot be
calculated. We can, however, approximate it by weighting the scores of the
observed event times by the inverse probability of censoring (Graf et al.,
1999; Gerds and Schumacher, 2006). This is called _inverse probability of
censoring weighting_ (IPCW), and the IPCW Brier score is given by

$\displaystyle\text{BS}_{\text{IPCW}}(t)=\frac{1}{n}\sum_{i=1}^{n}\left[\frac{\pi_{i}{(t)}^{2}\mathbbm{1}\\{T_{i}\leq
t,D_{i}=1\\}}{G_{i}(T_{i}-)}+\frac{{[1-\pi_{i}(t)]}^{2}\mathbbm{1}\\{T_{i}>t\\}}{G_{i}(t)}\right],$
(3)

where $G_{i}(t)=\textnormal{P}(C^{*}_{i}>t)>0$, is the survival function of
the censoring distribution for individual $i$. Assuming, for simplicity, that
the $G_{i}(t)$’s are known functions, the expectation of the IPCW Brier score
is

$\displaystyle\mathbb{E}\left[\text{BS}_{\text{IPCW}}(t)\right]$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\mathbb{E}\left[\frac{{\pi_{i}(t)}^{2}\mathbbm{1}\\{T_{i}\leq
t,D_{i}=1\\}}{G_{i}(T_{i}-)}+\frac{{[1-\pi_{i}(t)]}^{2}\mathbbm{1}\\{T_{i}>t\\}}{G_{i}(t)}\right]$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}{\pi_{i}(t)}^{2}\,\mathbb{E}\left[\frac{\mathbbm{1}\\{T_{i}^{*}\leq
t,T_{i}^{*}\leq
C_{i}^{*}\\}}{G_{i}(T_{i}^{*}-)}\right]+{[1-\pi_{i}(t)]}^{2}\,\frac{\textnormal{P}(T_{i}^{*}>t,C_{i}^{*}>t)}{G_{i}(t)}$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}{\pi_{i}(t)}^{2}\,\int_{0}^{t}\frac{G_{i}(u-)\,f_{i}(u)}{G_{i}(u-)}\,du+{[1-\pi_{i}(t)]}^{2}\,\frac{G_{i}(t)\,S_{i}(t)}{G_{i}(t)}$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}{\pi_{i}(t)}^{2}\,[1-S_{i}(t)]+{[1-\pi_{i}(t)]}^{2}\,S_{i}(t),$
$\displaystyle=\text{MSE}(t)+\frac{1}{n}\sum^{n}_{i=1}S_{i}(t)\left[1-S_{i}(t)\right].$

We see that this is identical to the expectation of the uncensored Brier score
(2), so the IPCW Brier score is a reasonable approximation of the uncensored
Brier Score.

####  3.1.1 IPCW with Administrative Censoring

For administrative censoring, _all_ the censoring times $C^{*}_{i}$ are known,
regardless of whether individual $i$ experienced an event or was censored.
Although $C^{*}_{i}$, in this case, is not random, we will continue to
capitalize $C^{*}_{i}$ for simplicity. As we know the censoring times, the
censoring distribution is given by the step-function
$G_{i}(t)=\mathbbm{1}\\{C^{*}_{i}>t\\}$, and the IPCW Brier score (3) can be
written as

$\displaystyle\text{BS}_{\text{IPCW}}(t)=\frac{1}{n}\sum_{i=1}^{n}\left[{\pi_{i}(t)}^{2}\mathbbm{1}\\{T_{i}\leq
t,D_{i}=1\\}+{[1-\pi_{i}(t)]}^{2}\mathbbm{1}\\{T_{i}>t\\}\right].$ (4)

This is equivalent to the unweighted Brier score (2) on the subset of
individuals that are not censored, meaning we simply disregard the set of
censored individuals $\\{i:\,T_{i}\leq t,D_{i}=0\\}$. The step-function
$G_{i}(t)=\mathbbm{1}\\{C^{*}_{i}>t\\}$ violates the assumptions of the IPCW
Brier score which, as stated in Section 3.1, requires $G_{i}(t)>0$. As a
result, the expectation of the score is no longer equal to that of the
uncensored Brier score,

$\displaystyle\mathbb{E}\left[\text{BS}_{\text{IPCW}}(t)\right]$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\mathbb{E}\left[{\pi_{i}(t)}^{2}\mathbbm{1}\\{T_{i}\leq
t,D_{i}=1\\}+{[1-\pi_{i}(t)]}^{2}\mathbbm{1}\\{T_{i}>t\\}\right]$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\left[{\pi_{i}(t)}^{2}\,\textnormal{P}(T_{i}\leq
t,D_{i}=1)+{[1-\pi_{i}(t)]}^{2}\,\textnormal{P}(T_{i}>t)\right]$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\left[{\pi_{i}(t)}^{2}\,\textnormal{P}(T_{i}^{*}\leq\min\\{t,C_{i}^{*}\\})+{[1-\pi_{i}(t)]}^{2}\,\textnormal{P}(T_{i}^{*}>t)\mathbbm{1}\\{C_{i}^{*}>t\\}\right]$

$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\left[{\pi_{i}(t)}^{2}\,\left[1-S_{i}(\min\\{t,C_{i}^{*}\\})\right]+{[1-\pi_{i}(t)]}^{2}\,S_{i}(t)\mathbbm{1}\\{C_{i}^{*}>t\\}\right]$
$\displaystyle=\frac{1}{n}\sum^{n}_{i=1}\left[{\pi_{i}(t)}^{2}\,[1-S_{i}(t)]+{[1-\pi_{i}(t)]}^{2}\,S_{i}(t)\right]\mathbbm{1}\\{C^{*}_{i}>t\\}$
$\displaystyle\quad+\frac{1}{n}\sum^{n}_{i=1}{\pi_{i}(t)}^{2}\,[1-S_{i}(C_{i}^{*})]\,\mathbbm{1}\\{C^{*}_{i}\leq
t\\}.$

We see that the expectation is only equal to that of the uncensored Brier
score if all individuals have censoring times $C^{*}_{i}>t$. For this group,
with $C^{*}_{i}>t$, we still have that the minimizers of the score are
$\pi_{i}(t)=S_{i}(t)$. However, for the other group, with $C^{*}_{i}\leq t$,
we see that the minimizers are $\pi_{i}(t)=0$. Hence, according to this score,
the optimal survival estimates are
$\pi_{i}(t)=S_{i}(t)\,\mathbbm{1}\\{C^{*}_{i}>t\\}$. This is problematic, as
we want an evaluation metric that decreases as our estimates approach the true
survival functions $S_{i}(t)$.

Due to these biases, one would, in practice, not use the Brier score (4) with
$G_{i}(t)=\mathbbm{1}\\{C^{*}_{i}>t\\}$, but instead use the regular IPCW
Brier score (3) with estimates of the censoring distributions. We will,
however, see that these issues can still be present. Specifically, if the
covariates hold sufficient information to identify the administrative
censoring times, then the estimated censoring distributions can approach the
step functions $G_{i}(t)=\mathbbm{1}\\{C^{*}_{i}>t\\}$, and we get a Brier
score close to that of (4).

####  3.1.2 Estimation of the Censoring Distribution

In practice, we need to estimate the censoring distributions $G_{i}(t)$ to use
the IPCW Brier score. This estimation is, in itself, a time-to-event
prediction problem, and can be addressed in the same manner as the original
time-to-event problem. Graf et al. (1999) proposed to use the Kaplan-Meier
estimates of the censoring distribution, and this is still the most common
approach. However, the Kaplan-Meier estimator disregards the covariates,
meaning all individuals are assumed to have the same censoring distribution.
This can lead to biased censoring estimates, as addressed by Gerds and
Schumacher (2006).

In predictive modeling, we typically split our data set into a training set
used to fit the models, a validation set used for hyperparameter tuning, and a
test set used for evaluating the models’ predictions. We, therefore, only
consider the Brier score calculated on a test set in this paper. Both Graf et
al. (1999) and Gerds and Schumacher (2006) use the test set to estimate the
censoring distribution, which is reasonable for simple models such as the
Kaplan-Meier estimator and Cox regression. If we, however, want to use more
flexible models, such as the Logistic-Hazard with neural networks, fitting to
the test set is likely to results in overfitted censoring estimates. To the
best of our knowledge, this topic has not been addressed in the literature, so
there are no “best practices” for how to approach such estimation problems. In
this paper, we will, therefore, treat the censoring distribution in the same
manner as the event-time distribution, meaning we fit the censoring model to
the training set, and use a validation set to verify that the estimates are
reasonable.

When the censoring distribution is obtained from the Kaplan-Meier estimator on
the test set, the weights of the IPCW Brier score (3) sum to $n$. But for more
flexible methods, and methods fitted to the training set, this is not
necessarily the case. A more general version of the IPCW Brier score (3), that
works for any reasonable censoring estimates, is obtained by replacing $n$
with the sum of the weights

$\displaystyle\tilde{n}(t)=\sum_{i=1}^{n}\left[\frac{\mathbbm{1}\\{T_{i}\leq
t,D_{i}=1\\}}{G_{i}(T_{i}-)}+\frac{\mathbbm{1}\\{T_{i}>t\\}}{G_{i}(t)}\right].$

This ensures that the Brier score is between 0 and 1,

$\displaystyle\text{BS}_{\text{IPCW}}(t)=\frac{1}{\tilde{n}(t)}\sum_{i=1}^{n}\left[\frac{\pi_{i}{(t)}^{2}\mathbbm{1}\\{T_{i}\leq
t,D_{i}=1\\}}{G_{i}(T_{i}-)}+\frac{{[1-\pi_{i}(t)]}^{2}\mathbbm{1}\\{T_{i}>t\\}}{G_{i}(t)}\right].$
(5)

This is the version of the IPCW Brier score we use in all our experiments.

The censoring estimates in (5) are actually functions of the covariates,
$G_{i}(t)=G(t\,|\,\mathbf{x}_{i})$. If there is enough information in the
covariates $\mathbf{x}_{i}$ to identify the censoring times $C^{*}_{i}$, the
censoring estimates will approach the step-functions
$G_{i}(t)=\mathbbm{1}\\{C^{*}_{i}>t\\}$, meaning the scores approach the
biased Brier score (4). As an example of this, consider a study where all
censoring is due to administrative censoring at the closure of the study, and
one includes the start date for each individual as a covariate, then the
$C^{*}_{i}$’s can be identified. A more realistic example would be that a
combination of certain covariates can identify the start date of some
individuals and, consequently, a subset of the censoring times can be
identified.

If we estimate the censoring distribution with a flexible method, we might
experience that some of the estimates of $G_{i}(t)$ become very small. This
corresponds to very large weights, meaning that a single individual can
potentially dominate the score. To prevent this from occurring, we set a
maximum weight allowed. As an example, if we have a maximum weight of 100, we
do not allow estimates of $G_{i}(t)<0.01$, giving the interpretation that a
single individual can maximally represent 100 individuals. In practice, this
is ensured by setting weights larger than 100 to be 100. By decreasing the
maximum allowed weight, we reduce the variance of the IPCW Brier score at the
expense of introducing some bias.

It should now be clear that the IPCW Brier score’s dependence on the censoring
estimates is undesirable, as we do not want an evaluation criterion to be
heavily dependent on the choices made by the researcher. If the scores are
substantially affected by the censoring estimates, one can question the
validity of the results.

###  3.2 The Administrative Brier Score

We have shown that the IPCW Brier score may have undesirable behavior under
administrative censoring, and that estimation of the censoring distribution
can be problematic in general. We propose an alternative that alleviates both
of these problems. Our approach is to use the uncensored Brier score (2) and
simply remove individuals from evaluation after their administrative censoring
time. This is possible as we know the censoring times $C^{*}_{i}$ for all
individuals. However, this also means that the score is not applicable in
studies where we do not know all the censoring times.

To define the administrative Brier score, first note that for $C^{*}_{i}\geq
t$ we know whether $T^{*}_{i}\leq t$ or $T^{*}_{i}>t$ (remember that if
$T^{*}_{i}=C^{*}_{i}$ we assume that we observe the occurrence of the event).
The administrative Brier score is then

$\displaystyle\text{BS}_{\text{A}}(t)$
$\displaystyle=\frac{1}{\tilde{n}(t)}\sum_{i=1}^{n}{\left[\mathbbm{1}\\{T^{*}_{i}>t\\}-\pi_{i}(t)\right]}^{2}\,\mathbbm{1}\\{C^{*}_{i}\geq
t\\},$ $\displaystyle=\frac{1}{\tilde{n}(t)}\sum_{i:\,C^{*}_{i}\geq
t}{\left[\mathbbm{1}\\{T^{*}_{i}>t\\}-\pi_{i}(t)\right]}^{2},$ (6)

where we scale by

$\displaystyle\tilde{n}(t)=\sum^{n}_{i=1}\mathbbm{1}\\{C^{*}_{i}\geq t\\},$

instead of $n$. The expectation of the administrative Brier score is

$\displaystyle\mathbb{E}\left[\text{BS}_{\text{A}}(t)\right]$
$\displaystyle=\frac{1}{\tilde{n}(t)}\sum_{i=1}^{n}\mathbb{E}\left[{\left[\mathbbm{1}\\{T^{*}_{i}>t\\}-\pi_{i}(t)\right]}^{2}\right]\,\mathbbm{1}\\{C^{*}_{i}\geq
t\\}$
$\displaystyle=\frac{1}{\tilde{n}(t)}\sum^{n}_{i=1}\left[{\pi_{i}(t)}^{2}\,[1-S_{i}(t)]+{[1-\pi_{i}(t)]}^{2}\,S_{i}(t)\right]\mathbbm{1}\\{C^{*}_{i}\geq
t\\}$ $\displaystyle=\frac{1}{\tilde{n}(t)}\sum_{i:\,C^{*}_{i}\geq
t}\left[{\pi_{i}(t)}^{2}\,[1-S_{i}(t)]+{[1-\pi_{i}(t)]}^{2}\,S_{i}(t)\right],$

which, for the subset of individuals with $C^{*}_{i}\geq t$, is identical to
the uncensored Brier score. Individuals with $C^{*}_{i}<t$ give no
contribution to the score. As we have no information about event times
$T^{*}_{i}>C^{*}_{i}$, it is reasonable to only consider this subset.

In practice, it is probably reasonable that a subset of the administrative
censoring times can be identified by the covariates. We saw in Sections 3.1.1
and 3.1.2 that the IPCW Brier score can be biased for this type of data. As
the administrative Brier score is still minimized by the true survival
functions $\pi_{i}(t)=S_{i}(t)$, it is more robust to these issues than the
IPCW Brier score.

##  4 Binary Classifiers for Time-to-Event Prediction

In machine learning, it is quite common to approach time-to-event prediction
with binary classifiers rather than relying on survival methodology.
Especially for someone without experience in survival analysis, a binary
classifier may seem like a reasonable approach to time-to-event modeling.

A binary classifier can be constructed for a given time $t$ by minimizing the
binary cross-entropy (negative log-likelihood for Bernoulli data) and
disregarding individuals that were censored before time $t$. This gives the
loss function,

$\displaystyle\text{loss}_{\text{BCE}}(t)$
$\displaystyle=-\sum^{n}_{i=1}\Big{(}\mathbbm{1}\\{T_{i}>t\\}\log[\pi_{i}(t)]+\mathbbm{1}\\{T_{i}\leq
t,\,D_{i}=1\\}\log[1-\pi_{i}(t)]\Big{)}$ (7)
$\displaystyle=-\sum^{n}_{i=1}\Big{(}y_{i}\log[\pi_{i}(t)]+(1-y_{i})\log[1-\pi_{i}(t)]\Big{)}\,\Big{(}1-\mathbbm{1}\\{T_{i}\leq
t,\,D_{i}=0\\}\Big{)}.$

were the labels $y_{i}=\mathbbm{1}\\{T_{i}>t\\}$ denote if the events happen
after time $t$. If we want survival estimates for a range of times
$\tau_{1}<\tau_{2}<\cdots<\tau_{m}$, we can fit a model for each $\tau_{j}$.
Alternatively, if we use a model for $\pi_{i}(\tau_{j})$ that can be estimated
for multiple $\tau_{j}$’s simultaneously, such as a neural network with $m$
output nodes, we can fit the model to the sum of the $m$ losses

$\displaystyle\text{loss}_{\text{BCE}}=-\sum_{j=1}^{m}\sum^{n}_{i=1}\Big{(}\mathbbm{1}\\{T_{i}>\tau_{j}\\}\log[\pi_{i}(\tau_{j})]+\mathbbm{1}\\{T_{i}\leq\tau_{j},\,D_{i}=1\\}\log[1-\pi_{i}(\tau_{j})]\Big{)}.$
(8)

We refer to this approach as the BCE method.

The binary classifiers and BCE method are clearly biased, as the removal of
censored individuals decreases the survival estimates. In fact, if there is
sufficient information in the covariates to identify the censoring times
$C^{*}_{i}$, the survival estimates of the binary classifiers will approach

$\displaystyle\pi_{i}(t)=S_{i}(t)\mathbbm{1}\\{C^{*}_{i}>t\\}.$ (9)

The derivations leading to (9) are given in Appendix A. We recognize these
estimates as the minimizers of the IPCW Brier score in Section 3.1.1. So given
that the $C^{*}_{i}$’s can be identified from the covariates and we have a
sufficiently large data set, the binary classifiers are essentially optimal
for minimizing the IPCW Brier score. If the covariates only identify a subset
of the $C^{*}_{i}$’s, it is not clear whether or not a binary classifier will
give better IPCW Brier scores than a corresponding survival method. We believe
this might explain why the BCE method got a lower Brier score on the KKBox
data set in Section 2.

If a researcher ignores censored individuals in the test set in the same
manner as for the binary classifiers (7), the uncensored Brier score (2) will
take the form of the IPCW Brier score in (4) with
$G_{i}(t)=\mathbbm{1}\\{C^{*}_{i}>t\\}$,

$\displaystyle\text{BS}(t)=\frac{1}{n}\sum_{i=1}^{n}\left[{\pi_{i}(t)}^{2}\mathbbm{1}\\{T_{i}\leq
t,D_{i}=1\\}+{[1-\pi_{i}(t)]}^{2}\mathbbm{1}\\{T_{i}>t\\}\right].$

We know from Section 3.1.1 that the expectation of this Brier score is
minimized by the estimates (9).

###  4.1 Simulation with BCE

To illustrate the concerns addressed in the previous sections, we conduct a
simple simulation study. We draw event times from a discrete event-time
distribution with 1000 discrete time points and constant hazard
$h_{i}(t)=0.00084$. The censoring times are drawn uniformly over the full
duration span from 0 to 100. We will use the censoring times in two distinct
ways. The first is to consider the censoring times as random, which
corresponds to a linear survival function $G_{i}(t)$ for the censoring,
identical for all individuals. The second approach is to consider the
censoring times identifiable from the covariates, meaning we have step-
functions $G_{i}(t)=\mathbbm{1}\\{C^{*}_{i}>t\\}$. In this case, we will use a
monotone function of $C^{*}_{i}$ as a covariate, which represents the
inclusion of covariates that can be used to identify the censoring time. For
brevity, we will denote this covariate as $C^{*}_{i}$, as it is equivalent to
using $C^{*}_{i}$ as a covariate. In summary, the only distinction between the
two approaches is whether or not the covariates contain information about the
censoring time $C^{*}_{i}$.

We fit two neural networks with the BCE loss (8) to 10,000 simulated samples,
the first without the covariate identifying $C^{*}_{i}$ and the second with
this covariate. In Figure 2 we have plotted the estimated survival function
for six distinct censoring times $C^{*}_{i}$, marked by the vertical red
dotted lines. The plots also contain the true survival function (in blue) for
comparison. From the orange lines, we see a clear bias of the BCE method
(binary classifier) applied to right-censored data, as it always
underestimates the survival. On the other hand, the survival estimates by the
BCE method with censoring information, represented by the green lines, follow
the true survival function reasonably well up till the censoring times and
then fall to zero. This agrees well with the optimal estimates in (9).

Figure 2: Survival estimates from the BCE method with and without a monotone
function of the censoring time $C^{*}_{i}$ as a covariate. The vertical dotted
red line gives the censoring time $C^{*}_{i}$. The true survival function is
plotted in blue. All event times are drawn from the same distribution.

Next, we investigate the resulting Brier scores of these survival estimates.
In Figure 3 we have plotted the uncensored Brier score (top left), the IPCW
Brier scores (5) with Kaplan-Meier estimates of $G_{i}(t)$ (top right), the
IPCW scores with the true censoring distribution
$G_{i}(t)=\mathbbm{1}\\{C^{*}_{i}>t\\}$ (bottom left), and the proposed
administrative Brier score (bottom right) given by (3.2). Recall that the IPCW
scores with the true censoring distribution (bottom left) correspond to
removing censored individuals in the same manner as that of a binary
classifier. The uncensored Brier score is computed on an uncensored test set,
while the three other scores are computed on a censored test set with the same
censoring distributions as in the training set.

Figure 3: Brier scores from simulations with administrative censoring. The top
left plot uses an uncensored test set, while the three other use the right-
censored test set. The orange lines represent a BCE method without any
information about the censoring times, while the green lines represent a BCE
method with sufficient covariate information to identify the censoring times.

From the uncensored Brier score, it is clear that the true survival function
does better than the two BCE methods. However, for both the IPCW Brier scores
we see that the BCE method with censoring information performs the best. In
fact, when we use the true censoring step-function as weights, even the BCE
without knowledge of the censoring time has better scores than the true
survival function. These figures, therefore, clearly illustrate a potential
weakness of the IPCW applied to data sets with identifiable censoring times.

The administrative Brier score (bottom right), rightfully identifies the true
survival function as better than the estimates of the BCE methods. For this
score, the BCE with censoring information performs very closely to the true
survival function, but this is expected as the survival estimates $\pi_{i}(t)$
are close to the true survival function for $C^{*}_{i}>t$.

The simulations were made to represent a worst-case scenario for the IPCW
Brier score, in the sense that the covariates contain information that makes
the censoring times $C^{*}_{i}$ easily identifiable for all individuals. It is
unlikely to find such extreme results in real data sets, but it serves the
purpose of illustrating the potential issues.

##  5 Simulations

The simulations in Section 4.1 illustrated the potential issues of using the
IPCW Brier score on a data set with administrative censoring. The simulations
were, however, tailored to emphasize such issues by having a very simple
survival function, and censoring times $C^{*}_{i}$ that were very easy to
identify from the covariates. In reality, the scores depend on how well the
BCE method is able to estimate both the survival and censoring distribution.
Hence, we conduct a more reasonable simulation study to investigate the extent
of these issues.

We will use the framework presented by Kvamme and Borgan (2019) to create
simulated data sets. This means that we draw event times by sequentially
sampling from discrete-time hazards on the time grid
$\\{0.1,0.2,\ldots,100\\}$. The hazards are specified through the logit-hazard
function $g(t\,|\,\mathbf{x})\in\mathbb{R}$. Note that $g(t\,|\,\mathbf{x})$
is just the notation used by Kvamme and Borgan (2019), and is not related to
$G_{i}(t)$. The discrete-time hazards are given by the sigmoid

$\displaystyle
h(t\,|\,\mathbf{x})=\frac{1}{1+\exp\left[-g(t\,|\,\mathbf{x})\right]}.$

The logit-hazards $g(t\,|\,\mathbf{x})$ are made up of a weighted sum of three
functions $g_{\text{sin}}(t\,|\,\mathbf{x})$,
$g_{\text{con}}(t\,|\,\mathbf{x})$, and $g_{\text{acc}}(t\,|\,\mathbf{x})$,

$\displaystyle g(t\,|\,\mathbf{x})$
$\displaystyle=\alpha_{1}\,g_{\text{sin}}(t\,|\,\mathbf{x})+\alpha_{2}\,g_{\text{con}}(t\,|\,\mathbf{x})+\alpha_{3}\,g_{\text{acc}}(t\,|\,\mathbf{x}),$
(10) $\displaystyle g_{\text{sin}}(t\,|\,\mathbf{x})$
$\displaystyle=\gamma_{1}\sin(\gamma_{2}[t+\gamma_{3}])+\gamma_{4},$
$\displaystyle g_{\text{con}}(t\,|\,\mathbf{x})$ $\displaystyle=\gamma_{5},$
(11) $\displaystyle g_{\text{acc}}(t\,|\,\mathbf{x})$
$\displaystyle=\gamma_{6}\cdot t-10,$ $\displaystyle\alpha_{i}$
$\displaystyle=\frac{\exp(\gamma_{i+6})}{\sum_{j=1}^{3}\exp(\gamma_{j+6})},\quad\text{for
}i=1,2,3.$

Here each of the nine $\gamma_{k}$’s is a function of the covariates
$\mathbf{x}$, meaning $\gamma_{k}$ is just a simplified notation of
$\gamma_{k}(\mathbf{x})$. The exact definitions of these
$\gamma_{k}(\mathbf{x})$’s are given in Kvamme and Borgan (2019, Appendix A.1)
which also explains the scheme used to draw the covariates. With $\tau_{j}$
denoting the $j$’th time point in $\\{0.1,0.2,\ldots,100\\}$, the survival
function is given by

$\displaystyle
S_{i}(\tau_{j})=S(\tau_{j}\,|\,\mathbf{x}_{i})=\prod_{k=1}^{j}\left[1-h(\tau_{k}\,|\,\mathbf{x}_{i})\right].$
(12)

To incorporate administrative censoring in the simulations, we consider a
monotonically decreasing function $Q(t\,|\,\mathbf{x}_{i})$, and let the
censoring time $C^{*}_{i}$ be defined by a threshold $\epsilon$ such that
$Q(C^{*}_{i}\,|\,\mathbf{x}_{i})=\epsilon$. This gives the survival function
of the censoring distribution

$\displaystyle G_{i}(t)=\mathbbm{1}\\{Q(t\,|\,\mathbf{x}_{i})>\epsilon\\}$
(13)

Hence, the censoring is deterministic, while the complexity of
$Q(t\,|\,\mathbf{x}_{i})$ controls the complexity of the relationship between
the covariates $\mathbf{x}_{i}$ and the censoring time $C^{*}_{i}$. We will
let $Q(t\,|\,\mathbf{x}_{i})$ have the same functional form as a survival
function, meaning it is defined in the same manner as (12), but with its own
set of covariates independent of the covariates of the event-time
distribution. In all the experiments we set $\epsilon=0.2$.

In the experiments, we compare the BCE method with the Logistic-Hazard method.
We use the Logistic-Hazard because of its similarity to the BCE method, in
that it also minimizes the binary cross-entropy, but use the discrete hazards
instead of the survival estimates. We use the implementation of the Logistic-
Hazard by Kvamme and Borgan (2019), but the method has been described by
multiple authors (Brown, 1975; Biganzoli et al., 1998; Gensheimer and
Narasimhan, 2019).

###  5.1 Complicated Censoring Distribution

In the first study, we draw event times using (10) with two covariates per
$\gamma_{k}$. The function $Q(t\,|\,\mathbf{x})$ is also defined using (10)
and (12), but with an independent set of covariates, ensuring that the event
times are independent of the censoring times. Note that although the
covariates contain enough information to identify the censoring times, the
complexity of the censoring distribution makes this a somewhat hard task.

We fit models using all the covariates from both the event-time distribution
and the censoring distribution, giving a total of 36 covariates. We draw
10,000 individuals for training and testing, and 4,000 for a validation set
used for early stopping of our training procedure. The networks are ReLU-nets
with 4 layers and 32 nodes in each layer. Batch normalization and dropout with
a rate of 0.1 are applied between each layer. The BCE and Logistic-Hazard both
give estimates for 50 equidistant times points between 0 and 100, but we
perform constant density interpolation (linear interpolation of survival
estimates, see Kvamme and Borgan, 2019) to obtain predictions for all 1,000
time points.

Figure 4: Brier scores from simulations with complicated administrative
censoring.

In Figure 4, we have plotted the Brier score of an uncensored test set, the
administrative Brier score (3.2), and four IPCW Brier scores using (5). The
IPCW scores are calculated with Kaplan-Meier estimates for the censoring
distribution, the true censoring function in (13), and censoring distributions
estimated with Logistic-Hazard and max weights of 100 and 1000. Recall from
Section 3.1.2 that when we estimate the censoring distribution with other
methods than Kaplan-Meier, the weights can become very large, resulting in
unstable results. Hence we set a max weight, which is given above the
respective plots.

From the figure, we see that both methods perform rather poorly compared to
the true survival function, and the Logistic-Hazard performs better than the
BCE for all scores except for the IPCW with the true censoring distribution.
The problems with IPCW do not appear here because the functional form of the
censoring distribution is quite complicated. Hence, the BCE method is not able
to identify the $C^{*}_{i}$’s to the extent that it can take advantage of the
bias of the IPCW scores.

Figure 5: Brier scores from simulations with simple administrative censoring.

###  5.2 Simple Censoring Distribution

In the second simulations study, we let $Q(t\,|\,\mathbf{x})$ be defined by
(12) with $g(t\,|\,\mathbf{x})=\gamma_{5}$ from (11) and let $\gamma_{5}$ be a
function of 5 covariates. This gives a very simple censoring distribution,
meaning that the censoring times $C^{*}_{i}$ are quite easily identifiable. We
repeat the simulations with the event-time distribution unchanged, and the
results are displayed in Figure 5. First, we note that the performance of the
BCE method on the uncensored test set is worse than before. This is expected
because it is now easier for the BCE method to identify the censoring time
$C^{*}_{i}$, meaning the survival estimates fall to zero right after
$C^{*}_{i}$ (as in Figure 2). Nevertheless, the BCE method now gets better
IPCW scores than the Logistic-Hazard. This shows that when the administrative
censoring function is simple to learn, the BCE method exhibit more of the step
behavior of the optimal estimates in (9). We observe, however, that the
administrative Brier score rightfully considers the BCE method worse than the
Logistic-Hazard. In this regard, we argue that even though the IPCW Brier
score might work well for administrative censoring, the administrative Brier
score is the safer choice.

In the simulations, we have considered a censoring distribution with
covariates that are independent of the survival covariates, so the event-time
distribution is independent of the censoring distribution. This might explain
why the Kaplan-Meier estimates work so well. If one would encounter this
independence for real data, one could simply remove the censoring covariates,
as they do not affect the event-time distribution, and use IPCW with Kaplan-
Meier weights.

##  6 KKBox Churn Prediction

Finally, we revisit the KKBox data set discussed in Section 2. We fit the BCE
method, corresponding to multiple binary classifiers, and the Logistic-Hazard
method. The training set is of size 100,000 and we use a validation set of
size 20,000 for early stopping. We use 6-layer ReLU networks with 128 nodes in
each layer, and with batch normalization and a dropout rate of 0.1 between
each layer. Entity embeddings are used for the categorical covariates (Guo and
Berkhahn, 2016) with embedding sizes that are the square root of the number of
categories. The outputs of the networks are of size 150 representing an
equidistant grid of the full time-scale of the training set. Constant density
interpolation (Kvamme and Borgan, 2019) is applied to the survival estimates
to obtain predictions for all the time-points in the test set. We use the
AdamWR (Loshchilov and Hutter, 2019) optimizer with a cycle length of 1 epoch,
but we double the cycle length and multiply the learning rate by 0.8 after
each cycle. Also, we do not include weight decay. The data set is,
essentially, the data set presented by Kvamme et al. (2019), but including all
censoring times and an extra categorical covariate stating the payment method.
The code for obtaining the data set is available at github.com/havakv/pycox.

Figure 6: Brier scores on the KKBox data set.

The censoring distribution is estimated in two ways. The first is with the
Kaplan-Meier estimator and the second is with a Logistic-Hazard with the same
hyperparameters as the Logistic-Hazard used to estimate the churn
distribution. The Brier scores are computed on a test set of size 100,000 and
displayed in Figure 6. The figure contains the two plots from Figure 1, the
administrative Brier scores (3.2), and three other IPCW scores (5) with
different max weights. Note that a maximum weight of 1 corresponds to the
unweighted score (4) where the set of censored individuals $\\{i:\,T_{i}\leq
t,D_{i}=0\\}$ are removed. This is the same way the BCE method (7) handles
censored individuals.

As we recall from Section 2, the IPCW Brier score with Kaplan-Meier estimates
considers the BCE method the best (top left Figure 6). We now see from the
administrative scores (top right) that this is probably a wrong conclusion, as
the BCE method has worse administrative Brier scores than the Logistic-Hazard.
We also see that the IPCW scores obtained with Logistic-Hazard are very
dependent on the maximum weight allowed in the scores. For higher allowed
weights, we would consider the BCE method to perform worse than the Logistic-
Hazard, but note that for the highest times this is not the case.

The administrative scores are simpler to obtain than the IPCW scores, as they
do not require estimation of the censoring distribution. Also, from this case
study, we see that the administrative scores are probably a safer chose of
evaluation metric than the IPCW weighted Brier score.

##  7 Discussion

In this paper, we have addressed potential issues of the inverse probability
of censoring weighted (IPCW) Brier score, in particular for administrative
censoring. If the covariates have sufficient information to determine the
administrative censoring times, the IPCW Brier score will not be minimized by
the true survival functions, but instead by a function that falls to zero
after the censoring time. We have also shown that a binary classifier that
disregards censored individuals will approach this minimizer. As a
consequence, the binary classifier might actually get better IPCW Brier scores
than the true survival function. Due to this bias, we argue that the IPCW
Brier score needs to be applied with care.

The IPCW Brier score can be substantially affected by the estimated censoring
distribution. If this is the case, the validity of the scores can be
questioned. In regards to both issues, we propose the _administrative Brier
score_ that works for administrative right-censored event times and does not
require estimation of the censoring distribution. This means that it is
simpler to use than the IPCW scores.

We simulate examples where the IPCW score fails to identify the best survival
estimates, but the administrative Brier score still provides reasonable
scores. We also investigate a real-world churn data set with administrative
censoring and find that it exhibits some of the same behavior as our
simulations. This shows that the proposed administrative Brier score can be a
very useful evaluation metric.

In this paper, we have only investigated the IPCW Brier score, but note that
there are multiple other IPCW scores that might suffer from the same drawbacks
as the Brier score. For example the IPCW Binomial log-likelihood (Graf et al.,
1999) and the IPCW concordance index (Uno et al., 2011; Gerds et al., 2013).
Preliminary investigations of the IPCW Binomial log-likelihood (Graf et al.,
1999) suggest that it has the same issues as the IPCW Brier score while an
administrative version of the Binomial log-likelihood behaves in the same
manner as the administrative Brier score.

The issues discussed in this paper are mostly relevant for machine learning
methods, such as neural networks, as the issues are only notable for quite
precise estimates of the event-time distribution and the censoring times
$C^{*}_{i}$. Hence, we are unlikely to encounter such issues with classical
statistical models.

Large parts of the machine learning literature rely on empirical evaluation of
predictive methodology. We, therefore, believe that more research on the
evaluation metrics for right-censored survival data is needed.

Acknowledgments

This work was supported by The Norwegian Research Council 237718 through the
Big Insight Center for research-driven innovation.

##  A The BCE Survival Estimates

To better understand the survival estimates of the binary classifiers in
Section 4, we investigate the minimizers of the expected loss. Again, we
stress that the BCE method corresponds to a set of binary classifiers in the
manner given by the loss (7) and (8).

The expected loss of the binary classifier (7) is

$\displaystyle\mathbb{E}\left[\text{loss}_{\text{BCE}}(t)\right]$
$\displaystyle=-\sum^{n}_{i=1}\Big{(}\textnormal{P}(T_{i}>t)\log\left[\pi_{i}(t)\right]+\textnormal{P}(T_{i}\leq
t,\,D_{i}=1)\log\left[1-\pi_{i}(t)\right]\Big{)}$
$\displaystyle=-\sum^{n}_{i=1}\Big{(}\textnormal{P}(T^{*}_{i}>t,\,C^{*}_{i}>t)\log\left[\pi_{i}(t)\right]+\textnormal{P}(T^{*}_{i}\leq
t,\,C^{*}_{i}\geq T^{*}_{i})\log\left[1-\pi_{i}(t)\right]\Big{)}$
$\displaystyle=-\sum^{n}_{i=1}\Big{(}S_{i}(t)G_{i}(t)\log\left[\pi_{i}(t)\right]+\left[\int_{0}^{t}G_{i}(u-)f_{i}(u)du\right]\,\log\left[1-\pi_{i}(t)\right]\Big{)}.$

The minimizers of this expectation with respect to $\pi_{i}(t)$ can be found
by equating the partial derivatives with zero,

$\displaystyle\frac{\partial\mathbb{E}\left[\text{loss}_{\text{BCE}}(t)\right]}{\partial\pi_{i}(t)}$
$\displaystyle=-\frac{S_{i}(t)\,G_{i}(t)}{\pi_{i}(t)}+\frac{\int_{0}^{t}G_{i}(u-)f_{i}(u)du}{1-\pi_{i}(t)}=0,$

This gives the minimizers

$\displaystyle\pi_{i}(t)$
$\displaystyle=\frac{S_{i}(t)G_{i}(t)}{S_{i}(t)G_{i}(t)+\int_{0}^{t}G_{i}(u-)f_{i}(u)du}$
(A.1)
$\displaystyle\leq\frac{S_{i}(t)G_{i}(t)}{S_{i}(t)G_{i}(t)+\int_{0}^{t}G_{i}(t)f_{i}(u)du}$
$\displaystyle\leq\frac{S_{i}(t)}{S_{i}(t)+\int_{0}^{t}f_{i}(u)du}$
$\displaystyle=S_{i}(t).$

We see that the $\pi_{i}(t)$’s are underestimating the true survival as long
as there is censoring present. This is expected as the binary classifiers
remove censored individuals, decreasing the proportion of survived
individuals.

If the censoring times $C^{*}_{i}$ can be identified by the covariates, the
censoring distribution is given by the step function
$G_{i}(t)=\mathbbm{1}\\{C^{*}_{i}>t\\}$. The minimizer (A.1) can then be
written as

$\displaystyle\pi_{i}(t)$
$\displaystyle=\frac{S_{i}(t)\mathbbm{1}\\{C^{*}_{i}>t\\}}{S_{i}(t)\mathbbm{1}\\{C^{*}_{i}>t\\}+\int_{0}^{t}\mathbbm{1}\\{C^{*}_{i}\geq
u\\}f_{i}(u)du}.$

If $C^{*}_{i}>t$, we have

$\displaystyle\pi_{i}(t)=\frac{S_{i}(t)}{S_{i}(t)+\int_{0}^{t}f_{i}(u)du}=S_{i}(t),$

and if $C^{*}_{i}\leq t$, we have $\pi_{i}(t)=0$. The minimizer can,
therefore, be written as

$\displaystyle\pi_{i}(t)=S_{i}(t)\mathbbm{1}\\{C^{*}_{i}>t\\}.$

We recognize this as the minimizer of the IPCW Brier score with administrative
censoring from Section 3.1.1. So if there is sufficient information in the
covariates to identify the administrative censoring times $C^{*}_{i}$, the
binary classifiers will approach the minimizers of the IPCW Brier score.

## References

  * Antolini et al. (2005) Laura Antolini, Patrizia Boracchi, and Elia Biganzoli.  A time-dependent discrimination index for survival data.  _Statistics in Medicine_ , 24(24):3927–3944, 2005.
  * Biganzoli et al. (1998) Elia Biganzoli, Patrizia Boracchi, Luigi Mariani, and Ettore Marubini.  Feed forward neural networks for the analysis of censored survival data: a partial logistic regression approach.  _Statistics in Medicine_ , 17(10):1169–1186, 1998.
  * Brown (1975) Charles C. Brown.  On the use of indicator variables for studying the time-dependence of parameters in a response-time model.  _Biometrics_ , 31(4):863–872, 1975.
  * Ching et al. (2018) Travers Ching, Xun Zhu, and Lana X. Garmire.  Cox-nnet: An artificial neural network method for prognosis prediction of high-throughput omics data.  _PLoS Computational Biology_ , 14(4):e1006076, 2018.
  * Fotso (2018) Stephane Fotso.  Deep neural networks for survival analysis based on a multi-task framework.  _arXiv preprint arXiv:1801.05512_ , 2018.
  * Gensheimer and Narasimhan (2019) Michael F. Gensheimer and Balasubramanian Narasimhan.  A scalable discrete-time survival model for neural networks.  _PeerJ_ , 7:e6257, 2019.
  * Gerds and Schumacher (2006) Thomas A. Gerds and Martin 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) Thomas A. Gerds, Michael W. Kattan, Martin Schumacher, and Changhong Yu.  Estimating a time-dependent concordance index for survival prediction models with covariate dependent censoring.  _Statistics in Medicine_ , 32(13):2173–2184, 2013\.
  * Graf et al. (1999) Erika Graf, Claudia Schmoor, Willi Sauerbrei, and Martin Schumacher.  Assessment and comparison of prognostic classification schemes for survival data.  _Statistics in Medicine_ , 18(17-18):2529–2545, 1999.
  * Guo and Berkhahn (2016) Cheng Guo and Felix Berkhahn.  Entity embeddings of categorical variables.  _arXiv preprint arXiv:1604.06737_ , 2016.
  * Harrell Jr et al. (1982) Frank E. Harrell Jr, Robert M. Califf, David B. Pryor, Kerry L. Lee, and Robert A. Rosati.  Evaluating the yield of medical tests.  _Journal of the American Medical Association_ , 247(18):2543–2546, 1982.
  * Ishwaran et al. (2008) Hemant Ishwaran, Udaya B. Kogalur, Eugene H. Blackstone, and Michael S. Lauer.  Random survival forests.  _Annals of Applied Statistics_ , 2(3):841–860, 2008.
  * Katzman et al. (2018) Jared L. Katzman, Uri Shaham, Alexander Cloninger, Jonathan Bates, Tingting Jiang, and Yuval Kluger.  Deepsurv: personalized treatment recommender system using a Cox proportional hazards deep neural network.  _BMC Medical Research Methodology_ , 18(1), 2018.
  * Kvamme et al. (2019) Håvard Kvamme, Ørnulf Borgan, and Ida Scheel.  Time-to-event prediction with neural networks and Cox regression.  _Journal of Machine Learning Research_ , 20(129):1–30, 2019.
  * Kvamme and Borgan (2019) Håvard Kvamme and Ørnulf Borgan.  Continuous and discrete-time survival prediction with neural networks.  _arXiv preprint arXiv:1910.06724_ , 2019.
  * Lee et al. (2018) Changhee Lee, William R Zame, Jinsung Yoon, and Mihaela van der Schaar.  Deephit: A deep learning approach to survival analysis with competing risks.  In _Thirty-Second AAAI Conference on Artificial Intelligence_ , 2018\.
  * Loshchilov and Hutter (2019) Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization.  In _International Conference on Learning Representations_ , 2019.
  * Luck et al. (2017) Margaux Luck, Tristan Sylvain, Héloïse Cardinal, Andrea Lodi, and Yoshua Bengio.  Deep learning for patient-specific kidney graft survival analysis.  _arXiv preprint arXiv:1705.10245_ , 2017.
  * Uno et al. (2011) Hajime Uno, Tianxi Cai, Michael J. Pencina, Ralph B. D’Agostino, and L. J. Wei.  On the C-statistics for evaluating overall adequacy of risk prediction procedures with censored survival data.  _Statistics in Medicine_ , 30(10):1105–1117, 2011.
  * Yousefi et al. (2017) Safoora Yousefi, Fatemeh Amrollahi, Mohamed Amgad, Chengliang Dong, Joshua E. Lewis, Congzheng Song, David A. Gutman, Sameer H. Halani, Jose Enrique Velazquez Vega, Daniel J. Brat, et al.  Predicting clinical outcomes from large scale cancer genomic profiles with deep survival models.  _Scientific Reports_ , 7(11707), 2017.


# Universal Learned Image Compression with Low Computational Cost

###### Abstract

Recently, learned image compression methods have developed rapidly and
exhibited excellent rate-distortion performance when compared to traditional
standards, such as JPEG, JPEG2000 and BPG. However, the learning-based methods
suffer from high computational costs, which is not beneficial for deployment
on devices with limited resources. To this end, we propose shift-addition
parallel modules (SAPMs), including SAPM-E for the encoder and SAPM-D for the
decoder, to largely reduce the energy consumption. To be specific, they can be
taken as plug-and-play components to upgrade existing CNN-based architectures,
where the shift branch is used to extract large-grained features as compared
to small-grained features learned by the addition branch. Furthermore, we
thoroughly analyze the probability distribution of latent representations and
propose to use Laplace Mixture Likelihoods for more accurate entropy
estimation. Experimental results demonstrate that the proposed methods can
achieve comparable or even better performance on both PSNR and MS-SSIM metrics
to that of the convolutional counterpart with an about 2$\times$ energy
reduction.

Index Terms—  Image compression, computational costs, energy consumption

##  1 Introduction

Image compression is one of the most fundamental fields in signal processing,
which aims to reach a trade-off between bitrate and distortion. With the rapid
development of deep learning, learned image compression methods [1, 2, 3, 4,
5] have drawn much attention and exhibited their outstanding rate-distortion
performance when compared to the traditional standards, such as JPEG [6],
JPEG2000 [7] and BPG [8]. Although numerous progress has been made, these
learning-based methods suffer from large computational complexity, which is
not conducive to deployment on devices with limited resources. Therefore, it
is necessary to propose a novel image compression framework with low
computational costs while maintaining comparable rate-distortion performance.

Fig. 1: Visualization of reconstructed images Kodim04 from Kodak dataset and
unit energy cost comparison.

To solve the problem, several methods have been developed. One of the most
popular methods is pruning [9], which aims to remove redundant or less
important weights to compress and accelerate the original network. Another
popular method is knowledge distillation [10], and the performance of
lightweight student networks can be improved through the knowledge conveyed by
the corresponding teacher networks. However, these compressed models still
contain massive multiplication operations and consume a lot of energy, and do
not benefit deployment. Although there exist model quantization methods [11]
that can reduce energy consumption, they cannot reach comparable performance
when compared to their baselines. Recently, Chen et al. [12] pioneered
AdderNet, which utilizes the $L1$ norm to calculate the similarities between
inputs and filters. This novel operation can avoid multiplication and reduce
massive energy consumption. In the same year, Elhoushi et al. [13] proposed
DeepShift, which quantizes each weight of convolutional networks to a power of
2. In the aspect of hardware, applying bit-shift operation on an element is
equivalent to multiplying it by a power of 2, and it means that DeepShift can
realize bit-shift operation through special quantization. Thus, both AdderNet
and DeepShift are suitable for deployment to resource-constrained devices,
such as mobile phones and embedded devices.

Recently, Li et al. [14] have proposed AdderIC, which utilizes AdderNet to
construct an image compression framework. However, it still exists a
performance gap between AdderIC and its CNN counterpart. To this end, in this
paper, we try to make full use of both AdderNet and DeepShift to propose a
novel image compression framework with low computational costs and comparable
rate-distortion performance. Our main contributions are as follows:

$\bullet$ We propose shift-addition parallel modules (SAPMs), including SPAM-E
for the encoder and SAPM-D for the decoder, to largely reduce energy
consumption and facilitate deployment on those computing-constrained devices.
Besides, the proposed SAPMs are plug-and-play components, which means that
they have strong generalization ability and can be applied in other CNN-based
architectures.

$\bullet$ We thoroughly analyze the probability distribution of latent
representations and propose a Laplace Mixture Model for more accurate entropy
estimation.

$\bullet$ Experimental results show that the proposed methods can bring about
2$\times$ energy reduction while maintaining comparable or even better rate-
distortion performance on both PSNR and MS-SSIM metrics as shown in Fig. 1.

##  2 Proposed Method

###  2.1 Motivation

From [15, 16, 17], we can get energy consumption among different operations
and formats when implemented in a 45nm CMOS technology as follow:

Table 1: Energy consumption comparison Operation | Format | Energy Cost (pJ) | Improv.
---|---|---|---
# Mult. | FP32 | 3.70 | -
# Add. | FP32 | 0.90 | 4.1$\times$
FIX32 | 0.10 | 37$\times$
# Shift | FIX32 | 0.13 | 28.5$\times$

From Tab. 1, both shift and addition can bring a cheap operation, while
multiplication suffers from extremely large energy consumption. Besides,
through experimental observation, neither shift nor addition can achieve
similar performance as multiplication-based networks (e.g., CNN). To this end,
we try to combine these two weak players and develop low-power shift-addition
parallel modules (SAPMs) for image compression, where the shift branch can be
used to extract large-grained features and the addition branch can learn
small-grained features.

###  2.2 Shift-Addition Parallel Modules

####  2.2.1 Shift Branch

DeepShift was first proposed in [13], and it adopts bitwise shift operation to
construct networks. Specifically, there are two methods to train the DeepShift
model, including DeepShift-Q and DeepShift-PS. DeepShift-Q quantizes the
weight $W$ in networks to $W_{q}$ by rounding it to the nearest power of 2,
while DeepShift-PS trains shift and sign parameters directly. In this paper,
DeepShift-Q is adopted to train our shift branch and the forward pass as
follow:

$S=\operatorname{sign}(W),P=\log_{2}(\operatorname{abs}(W)),W_{q}=S\cdot
2^{p}\vspace{-0.1cm}$ (1)

where $S$ denotes sign matrix and $P$ represents shift matrix. Once quantized,
it can realize shift operation since applying bitwise shift on an element is
equivalent to multiplying it by a power of 2 in hardware.

Although bitwise shift operation exhibits its cheap cost in terms of energy
consumption, it cannot achieve the same or similar expressive capacity as its
original counterpart. This is because that shift branch will ignore some
important information during quantization, which means it can merely be used
to extract coarse-grained features.

Fig. 2: Probability distribution of latent representation $y$.

####  2.2.2 Addition Branch

AdderNet was first proposed by Chen et al. [12], which utilizes $L1$ norm
rather than cross-correlation to calculate the similarities between inputs and
filters. Due to the properties of the $L1$ norm, AdderNet can avoid expensive
multiplication and only involves addition or subtraction. What is more, the
vanilla AdderNet adopts a batch normalization (BN) after each adder layer to
prevent gradients from exploding. However, for the image compression task, it
is not helpful to improve performance once introducing the BN layer since it
cannot reduce pixel-wise redundancies. Thus, it is necessary to develop a new
scheme that can not only prevent gradients from exploding, but also help to
reduce spatial redundancies.

Fig. 3: General architecture. The detailed structures of the adder blocks are
shown on the right, while the shift block is a shift layer in the encoder and
a shift layer follow by a PixelShuffle layer [18] in the decoder. Context
model is composed of a 5$\times$5 mask convolution [19], while the structures
of hyper encoder, hyper decoder and factorized entropy model are the same as
[4].

Motivated by [20], we adopt 1$\times$1 implicitly deconvolution (ID) to
replace the BN after each adder layer. The core of it is to calculate the
covariance matrix $Cov$ of input $X$ and then get the deconvolution operation
$D$:

$\operatorname{\textit{Cov}}=\frac{1}{N}(X-\mu)^{T}(X-\mu),D=\operatorname{\textit{Cov}}^{-\frac{1}{2}}\vspace{-0.1cm}$
(2)

where $N$ and $\mu$ represent the number of samples and the mean of input,
respectively. Once $D$ is obtained, we can apply it after each centered input
$(X-\mu).$ Furthermore, according to the association rule of matrix
multiplication, we can implicitly obtain the output $Y$ by changing the
weights $W$ of networks as follow:

$Y=(X-\mu)\cdot D\cdot W=X\cdot(D\cdot W)-\mu\cdot D\cdot W\vspace{-0.1cm}$
(3)

Although the computation cost of 1$\times$1 ID is slightly higher than the BN,
it is still significantly lower than that of regular convolution layers (e.g.,
5$\times$5 Conv).

Besides, two kinds of shortcut connections are developed for powerful
expressive capacity, and the details will be discussed in Sec. 2.4. Similar to
CNN, the addition branch can extract fine-grained features and capture details
in images.

###  2.3 Laplace Mixture Model

Entropy models utilize parameterized distribution to fit the marginal
distribution of quantized latent representation $\hat{y}$. In other words, a
more accurate entropy model can save more bits, which can help to improve
rate-distortion performance significantly. Ballé et al. [2] firstly proposed a
univariate Gaussian distribution model for the hyperprior, and later works [3]
extended it to a mean and scale Gaussian distribution for accurate entropy
estimation. Furthermore, Cheng et al. [4] developed a Gaussian Mixture Model,
which utilizes $k$ Gaussian distributions to fit the marginal distribution and
demonstrates strong fitting abilities. In a word, these works are based on an
assumption that the prior distribution of latent representation is a Gaussian
distribution. However, whether the prior distribution is also Gaussian in the
proposed frameworks is unclear, and it should be explored.

We conduct several experiments and randomly select one of the channels to
visualize the probability distribution of latent representation in Fig. 2.
From the figure, we can see that the probability distribution of latent
representation is more Laplace than Gaussian, which means that the entropy
models used in CNN frameworks cannot fit the marginal distribution well. To
this end, we propose Laplace Mixture Model (LMM) to estimate the marginal
distribution more accurately.

$p_{\hat{y}\mid
z}(\hat{y}\mid\hat{z})\sim\sum_{k=1}^{K}w^{(k)}\mathcal{L}\left(\mu^{(k)},\sigma^{(k)}\right)\vspace{-0.2cm}$
(4)

where $\mu^{(k)}$, $\sigma^{(k)}$, $w^{(k)}$denote the mean, scale and weight
of the $k$-th mixture (Laplace distribution) respectively. Once LMM is
developed, it can help to save more bits during arithmetic encoding and
decoding as shown in Fig. 3

###  2.4 Overall Structure

Fig. 3 shows the general architecture. The backbone network is based on [1, 2]
and there are $n$ levels of transformation in the encoder or decoder, where
the convolution layer and transposed convolution layer are replaced by SAPM-E
and SAPM-D, respectively. In SAPM-E, an adder block is composed of an adder
layer followed by a 1$\times$1 ID. At the same time, we add an average pooling
to improve the expressive capacity of the addition branch. In SAPM-D,
1$\times$1 ID is used to increase the number of output channels, and the
upsampling can be realized by the PixelShuffle layer [18]. Similar to the role
of pooling in the encoder, we also add a 1$\times$1 Conv as a shortcut
connection to enhance the expressive capacity.

##  3 Experiments

###  3.1 Experimental Setup

Fig. 4: Performance evaluation on Kodak dataset.

We adopt 3 levels of transformation in the encoder or decoder in our
experiments. Specifically, for better performance, we employ a convolutional
layer followed by two SAPM-Es in the encoder, and two SAPM-Ds followed by a
transposed convolutional layer in the decoder, which is similar to the
settings of vanilla AdderNet. Besides, we choose generalized divisive
normalization (GDN) [1] as our nonlinear module. Other settings such as hyper
encoder, hyper decoder, and context model are consistent with [4]. After that,
the model is trained and optimized for mean squared error (MSE) using two sets
of $\lambda$ values (i.e., low bitrate: $\\{16,32,64,128,256,512\\}$, high
bitrate: $\\{1024,2048,4096\\}$) on the CLIC2020 training dataset [21] which
consists of 61894 images with 256$\times$256 pixels, and tested on the
standard Kodak dataset [22] with 24 images of 512$\times$768 or 768$\times$512
pixels. To be specific, in the first stage, we train the models with the
largest $\lambda$ value in each set with a batch size of 8, and then apply
Adam optimizer with the learning rate of $1\times 10^{-4}$ in the first
900,000 iterations and $1\times 10^{-5}$ in the remaining 100,000 iterations.
For other bitrates ($\lambda$ values), we adopt the models trained on bitrate
$\lambda=4096$ and $\lambda=512$ as pre-trained models and then fine-tune the
rest models for 320,000 iterations with a learning rate of $1\times 10^{-5}$.
At last, we adopt average bit-per-pixel (bpp), average peak signal-to-noise
ratio (PSNR), average multi-scale structural similarity (MS-SSIM) as well as
energy consumption as the metrics to evaluate the performance.

###  3.2 Performance Comparison

The rate-distortion performance on the Kodak dataset is shown in Fig. 4. We
compare our method with the corresponding CNN baseline and several traditional
standards, including JPEG [6], JPEG2000 [7] and BPG [8]. Regarding PSNR, the
proposed method largely outperforms JPEG and JPEG2000, while showing extremely
comparable performance with the CNN baseline and BPG. When comes to MS-SSIM,
our model greatly surpasses JPEG, JPEG2000 and BPG, and even exhibits better
performance than the CNN baseline at some high bitrate points. Moreover,
similar to the CNN counterpart, our method also produces visually pleasant
reconstructed images as shown in Fig. 1.

At last, we calculate the energy consumption of the proposed SAPMs according
to Tab. 1. Due to the low computation cost of shortcut connections and
1$\times$1 ID, we omit them and then make a comparison with the convolutional
layer. From Tab. 2, the proposed SAPMs can bring more than 2$\times$ energy
reduction compared to the convolutional layer, which facilitates deployment on
devices with limited resources.

Table 2: Unit energy comparison between CNN and SAPM. Operation | Format | Convolution | SAPM
---|---|---|---
# Mult. | FP32 | $3.70\times 1$ | 0
# Add. | FP32 | $0.90\times 1$ | $0.90\times 2$
FIX32 | 0 | $0.10\times 1$
# Shift | FIX32 | 0 | $0.13\times 1$
Total Energy Costs (pJ) | 4.60 | 2.03 (2.2$\times$$\downarrow$)

###  3.3 Ablation Studies

To verify the effectiveness of the proposed methods, several ablation
experiments are carried out as shown in Fig. 4. To begin with, either shift
branch or addition branch is less capable compared to their CNN baseline,
while integrating these two weak players can bring strong expressive capacity.
Furthermore, the proposed LMM achieves a more accurate entropy estimation for
the marginal distribution, which can save more bits and further improve the
rate-distortion performance. What is more, it is worth noting that the vanilla
AdderNet cannot be directly used in image compression because it would cause
great performance degradation and its rate-distortion (RD) curve cannot be
displayed in Fig. 4.

##  4 Conclusion

In this paper, we proposed SAPMs to largely reduce energy consumption in
learned image compression. Besides, we thoroughly analyze the probability
distribution of the latent representation and then develop a novel Laplace
Mixture Model for more accurate entropy estimation. At last, several
experiments are carried out to demonstrate that the proposed model can achieve
comparable or even better performance on both PSNR and MS-SSIM metrics to that
of the CNN counterpart while reducing more than 2$\times$ energy consumption.

## References

  * [1] Johannes Ballé, Valero Laparra, and Eero P. Simoncelli,  “End-to-end optimized image compression,”  in 5th International Conference on Learning Representations, ICLR 2017, 2017.
  * [2] Johannes Ballé, David Minnen, Saurabh Singh, Sung Jin Hwang, and Nick Johnston,  “Variational image compression with a scale hyperprior,”  in International Conference on Learning Representations, 2018.
  * [3] David Minnen, Johannes Ballé, and George Toderici,  “Joint autoregressive and hierarchical priors for learned image compression,”  Advances in Neural Information Processing Systems, vol. 2018-Decem, no. Nips, pp. 10771–10780, 2018.
  * [4] Zhengxue Cheng, Heming Sun, Masaru Takeuchi, and Jiro Katto,  “Learned image compression with discretized gaussian mixture likelihoods and attention modules,”  Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, pp. 7936–7945, 2020.
  * [5] Yueqi Xie, Ka Leong Cheng, and Qifeng Chen,  Enhanced Invertible Encoding for Learned Image Compression, vol. 1,  Association for Computing Machinery, 2021.
  * [6] Gregory K Wallace,  “The jpeg still picture compression standard,”  IEEE transactions on consumer electronics, vol. 38, no. 1, pp. xviii–xxxiv, 1992.
  * [7] Majid Rabbani and Rajan Joshi,  “An overview of the JPEG 2000 still image compression standard,”  Signal Processing: Image Communication, vol. 17, no. 1, pp. 3–48, jan 2002.
  * [8] F. Bellard,  “BPG image format,” Apr. 2018. [Online]. Available: https://bellard.org/bpg/.
  * [9] Huizi Mao, Song Han, Jeff Pool, Wenshuo Li, Xingyu Liu, Yu Wang, and William J Dally,  “Exploring the granularity of sparsity in convolutional neural networks,”  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops, 2017, pp. 13–20.
  * [10] Geoffrey Hinton, Oriol Vinyals, and Jeff Dean,  “Distilling the knowledge in a neural network,”  arXiv preprint arXiv:1503.02531, 2015.
  * [11] Matthieu Courbariaux, Itay Hubara, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio,  “Binarized Neural Networks: Training Deep Neural Networks with Weights and Activations Constrained to +1 or -1,”  arXiv preprint arXiv:1602.02830, 2016.
  * [12] Hanting Chen, Yunhe Wang, Chunjing Xu, Boxin Shi, Chao Xu, Qi Tian, and Chang Xu,  “Addernet: Do we really need multiplications in deep learning?,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2020, pp. 1468–1477.
  * [13] Mostafa Elhoushi, Zihao Chen, Farhan Shafiq, Ye Henry Tian, and Joey Yiwei Li,  “Deepshift: Towards multiplication-less neural networks,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2021, pp. 2359–2368.
  * [14] Bowen Li, Yao Xin, Chao Li, Youneng Bao, Fanyang Meng, and Yongsheng Liang,  “Adderic: Towards low computation cost image compression,”  in ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2022, pp. 2030–2034.
  * [15] Mark Horowitz,  “1.1 computing’s energy problem (and what we can do about it),”  in 2014 IEEE International Solid-State Circuits Conference Digest of Technical Papers (ISSCC). IEEE, 2014, pp. 10–14.
  * [16] William Dally,  “High-performance hardware for machine learning,”  NIPS Tutorial, 2015.
  * [17] Haoran You, Xiaohan Chen, Yongan Zhang, Chaojian Li, Sicheng Li, Zihao Liu, Zhangyang Wang, and Yingyan Lin,  “Shiftaddnet: A hardware-inspired deep network,”  arXiv preprint arXiv:2010.12785, 2020.
  * [18] Wenzhe Shi, Jose Caballero, Ferenc Huszar, Johannes Totz, Andrew P. Aitken, Rob Bishop, Daniel Rueckert, and Zehan Wang,  “Real-Time Single Image and Video Super-Resolution Using an Efficient Sub-Pixel Convolutional Neural Network,”  Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, vol. 2016-Decem, pp. 1874–1883, 2016.
  * [19] Aaron Van den Oord, Nal Kalchbrenner, Lasse Espeholt, Oriol Vinyals, Alex Graves, et al.,  “Conditional image generation with pixelcnn decoders,”  Advances in neural information processing systems, vol. 29, 2016\.
  * [20] C Ye,  “Network deconvolution,”  in International Conference on Learning Representations, 2020.
  * [21] Toderici George, Shi Wenzhe, Timofte Radu, Theis Lucas, Balle Johannes, Agustsson Eirikur, Nick Johnston, and Mentzer Fabian,  “Workshop and challenge on learned image compression (clic2020),” 2020\.
  * [22] Rich Franzen,  “Kodak lossless true color image suite (photocd pcd0992),” http://r0k.us/graphics/kodak/,  1999\.


11institutetext:  Department of Computer Science, Open University of the
Netherlands
22institutetext: Department of Multimedia Engineering, Osaka University, Japan
33institutetext: ICIS, Radboud University, The Netherlands

# What is understandable in Bayesian network explanations?

Raphaela Butz 11 Renée Schulz 22 Arjen Hommersom 1133 Marko van Eekelen 1133

###### Abstract

Explaining predictions from Bayesian networks, for example to physicians, is
non-trivial. Various explanation methods for Bayesian network inference have
appeared in literature, focusing on different aspects of the underlying
reasoning. While there has been a lot of technical research, there is very
little known about how well humans actually understand these explanations. In
this paper, we present ongoing research in which four different explanation
approaches were compared through a survey by asking a group of human
participants to interpret the explanations.

##  1 Introduction

Bayesian networks (BNs) gained rising interest in AI in medicine in the last
30 years. Part of the interest in Bayesian networks in this particular domain
lies in the apparent transparency and interpretability of models [1]. However,
it has been widely recognised that explaining the inference of Bayesian
networks is difficult as they do not imitate human reasoning [6], which may
hinder their applicability in the health care domain. Therefore, there have
been many proposals for explaining Bayesian networks in the field of
healthcare, e.g. [2, 5, 10] as well as in other fields [7],[8],[9],[11]. This
shows that explaining BNs is not a trivial task and raises the question which
methods or aspects of methods are actually useful in practice.

The goal of this paper is to compare and analyse the human understandability
of BNs using recent explanation methods from Vlek [8], Kyrimi [5], Timmer [7]
and Butz [2], which are based on [10] [9] [11]. We first briefly introduce the
explanation methods and our approach in Section 2. We then discuss the survey
that we used in Section 3. In Section 4 the results are evaluated and the
conclusion of our findings can be found in Section 5.

##  2 Methodology

For the survey, a mixed methods (qualitative and quantitative) approach was
evaluated to be necessary. We require qualitative (e.g. opinions, questions,
thoughts, and misunderstandings) as well as quantitative data (demographics,
self evaluation of understanding, and comparison of models on a scale).
Interviews as well as an online version of the interview in a survey format
were conducted. The interview and survey feature the same set of questions and
both methods are based on open-ended answer possibilities including a few
scale-based and multiple choice questions.

Even though this research is done to improve the use of explained BNs in the
medical and healthcare field, an approach was needed to broaden the range of
possible participants. Therefore, we did not use a BN with medical data to
minimize the necessity of background knowledge of the participants.
Furthermore, the survey was designed using a gamified approach to ease the
cognitive load on the participants’ end. Gamification is defined as the use of
game design elements used in non-game context [3]. In this survey, we chose
the element of storytelling as well as contextualised (to the story)
visualizations. The goal with this approach is to gain a generalisable outcome
to improve the explanation of BNs in all contexts, including healthcare.

For the analysis of the free text answers, the qualitative content analysis
process [4] was chosen. The coding of the data was done with a deductive
approach, meaning we chose analysing categories based on previous research.
There were two sets of categories that were used for the analysis. The first
set of categories features the ability of participants to come to a conclusion
based on the different explanation methods, with the subcategories of ’having
doubts’ or ’no doubts’ about the conclusion. It should be noted that only the
participants’ perception and not their actual understanding of the methods
were measured. The second set of categories describes the participants’
expressions (positive and negative) about specific aspects of the explanation
methods (e.g. certain visualization forms, structure or type of language).
Furthermore, the coding of the data was done by two experts individually for
comparable sets of coded data to ensure quality and further iteration of the
categories as needed. The categories as stated in this paper were refined in 3
iterations, developed through the analysis of the content and deeper
understanding of user behavior.

##  3 Case Study and Survey Design

For the design of the case study, several explanation methods were chosen
which we labeled with the names _Tree, Scenario, Table, Language_. The _Tree_
method by Timmer et al. [7] is using reasoning chains to explain the
derivation of probabilities resulting from the BN. This method requires two
phases. In the first phase, a support graph is constructed from the identified
chains of reasoning featuring the target variable as a root node. In the
second phase, arguments are subsequently constructed and presented in a
(visual) tree structure for the users to interpret.
Vlek et al. [8] _Scenario_ method uses Scenario nodes to evaluate the
coherence of a scenario in a BN. Scenarios, as modeled in this method, add a
narrative component to allow experts to gain overview over the entirety of the
statistical information within a BN while being presented with smaller
scenarios highlighting certain aspects of the BN’s information. The scenario
nodes in the BN allow for the calculation of how the evidence supports or
refutes each scenario node.
Kyrimi and Marsh [5] method’s structure consists of three successive levels of
detail featuring an incremental explanation. These levels are _1: Significant
evidence variables_ , _2: Information flow_ and _3: Significant evidence
impact on the intermediate variables_. Level 1 determines in which way the
evidence affects the variable of interest. Level 2 shows the intermediate
variables through which the information passes in the BN. Level 3 shows in
which way the evidence affects the variables through which the information
passes.
The last method, _Language_ , is describing the most probable explanation in
natural language and has been introduced in Butz, Hommersom and Van Eekelen
[2]. A argumentative probability tree is built by searching for the strongest
dependency starting from the evidence. To generate natural language from the
argumentative probability tree, boiler plates were introduced and generated in
a uniform structure due to the argumentative probability tree.
For all presented explanation methods, the BN fom Vlek [8] has been used as a
basis to have a comparable set of explanations. This means that the same data
was presented using the different explanation methods. During the design of
the final representations of each method, we found that in some cases,
creating such examples is not straightforward.
Some of the explanation methods feature technical terminology of its own and
therefore are not generally understandable without further explanation. One
example for this can be found in _Table_ [5], as the final representation
features a certain structure of ’levels’, however the potential reader of the
results does not actually know what that means. This means that some methods
generate another structure that is in need of explanation. To be able to
include these methods in the comparison, we included a base set of
explanations on how to read and understand certain parts of the
representation, e.g. ’levels’ in _Table_ and the indicator called ’infinite’
in _Tree_ [7].
For the survey, we used Vlek’s BN [8] on a homicide case to provide an equal
data set in order to prevent the results from being biased by data sets of
different quality. This was done with the intention that a non-specialized
data set can be understood by everyone, without specialized knowledge, as it
would have been necessary in the medical domain. An example for an explanation
method based on a medical condition, that requires specialized knowledge to
understand the explanation, can be found in [5]. In the survey, we asked the
participants to explain each method in a free text form or took notes while
interviewing. At the end, we gave them the opportunity to reflect on the
methods and to assess for themselves how informative they perceived them.

##  4 Survey Results

In total, 20 people from different professions participated in the survey. Out
of these 20 participants, 3 claimed having basic knowledge about BNs, 9 heard
of them, and 8 have never heard of them before. Participants indicated that
the BN itself and _Scenario_ were the easiest to understand with 11 votes,
while _Table_ was voted the most confusing with 11 votes followed by
_Language_ with 8 votes. Likewise, 11 people stated that the BN gives the most
insight, while 4 voted for _Scenario_. _Language_ came in third with 3 votes.

Some issues were frequently raised in some methods, such as missing
explanations about the meaning of arrows and how visualisations were generated
was asked 7 times. These issues were often mentioned in combination with
missing explanations for colour coding (6 times) and the percentage indicators
(9 times). It was stated 5 times that in _Tree_ , the labeling with ’infinite’
was confusing. Furthermore, the difference between the terms ”infinite” and
”certain” was not clear. That _Table_ is confusing was stated 10 times, 5
times with a focus on the levels and 4 times the average scenario was not
understood. That the language used in _Language_ was not natural and thus
confusing was mentioned 6 times, whilst 4 times it was voted the easiest to
understand. 2 times, however, it was stated that the text was too long and
that participants preferred information presented more to the point and
shorter.

The results of the qualitative content analysis can be found in Figure 1. The
most prominent results are that the participants had the least doubts about
their conclusion in _Scenario_ (12 people), whereas 13 people could not deduct
a conclusion from _Table_. Another interesting result is that the explanatory
ability of the BN hardly differs from that of _Language_ in the presented
example.

Figure 1: Results of the qualitative content analysis by categories.

##  5 Conclusion

A conclusion that can be drawn for all methods is that it is not only about
developing a method to explain a BN, but it is also essential to state what
prior knowledge is required to understand the explanation. The evaluation of
_Table_ was particularly difficult for the authors because the paper did not
specify what prior knowledge was necessary to understand the results of the
method. Furthermore, it was noticeable that many people spoke of scenarios (13
times), not only in _Scenario_ , which suggests that we like to think in
scenarios and want to see them confirmed or refuted.

Most participants stated _Scenario_ was the easiest to understand after the
BN. This result may be biased because the example that was used in the case
study was designed for the _Scenario_ method. Other types of data may have
been more suitable for other explanation methods.

In addition, it seems that a combination of visualization and simple, but
precise text helps most participants to make decisions. Natural language (and
longer notes or sentences) seems to be in a disadvantage for grasping
information quick and precisely in a written format. However, as this research
continues, we are aware that switching the output medium to e.g. audio might
have further impact on the usefulness of a representation.

## References

  * [1] Arrieta, A.B., Díaz-Rodríguez, N., Del Ser, J., Bennetot, A., Tabik, S., Barbado, A., García, S., Gil-López, S., Molina, D., Benjamins, R., et al.: Explainable artificial intelligence (XAI): Concepts, taxonomies, opportunities and challenges toward responsible ai. Information Fusion 58, 82–115 (2020)
  * [2] Butz, R., Hommersom, A., van Eekelen, M.: Explaining the most probable explanation. In: Scalable Uncertainty Management - 12th International Conference, SUM 2018, Milan, Italy, October 3-5, 2018, Proceedings. pp. 50–63 (2018), https://doi.org/10.1007/978-3-030-00461-3_4
  * [3] Deterding, S., Dixon, D., Khaled, R., Nacke, L.: From game design elements to gamefulness: defining” gamification”. In: Proceedings of the 15th international academic MindTrek conference. pp. 9–15 (2011)
  * [4] Elo, S., Kyngäs, H.: The qualitative content analysis process. Journal of Advanced Nursing 62(1), 107–115 (2008), https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1365-2648.2007.04569.x
  * [5] Kyrimi, E., Marsh, W.: A progressive explanation of inference in ’hybrid’ Bayesian networks for supporting clinical decision making. Proceedings of the Eighth International Conference on Probabilistic Graphical Models (Jan 2016)
  * [6] Lacave, C., Díez, F.J.: A review of explanation methods for Bayesian networks. The Knowledge Engineering Review 17(2), 107–127 (2002)
  * [7] Timmer, S.T., Meyer, J.J.C., Prakken, H., Renooij, S., Verheij, B.: A two-phase method for extracting explanatory arguments from Bayesian networks. International Journal of Approximate Reasoning 80, 475 – 494 (2017)
  * [8] Vlek, C.S., Prakken, H., Renooij, S., Verheij, B.: A method for explaining Bayesian networks for legal evidence with scenarios. Artificial Intelligence and Law 24(3), 285–324 (Sep 2016)
  * [9] Vreeswijk, G.A.W.: Argumentation in Bayesian Belief Networks, pp. 111–129. Springer Berlin Heidelberg, Berlin, Heidelberg (2005)
  * [10] Williams, M., Williamson, J.: Combining argumentation and Bayesian nets for breast cancer prognosis. Journal of Logic, Language and Information 15(1), 155–178 (2006)
  * [11] Yap, G.E., Tan, A.H., Pang, H.H.: Explaining inferences in Bayesian networks. Applied Intelligence 29(3), 263–278 (Dec 2008), http://dx.doi.org/10.1007/s10489-007-0093-8