# DySurv: Dynamic Deep Learning Model for Survival Prediction in the ICU

Munib Mesinovic, Peter Watkinson, Tingting Zhu M. Mesinovic is with the
Department of Engineering Science, University of Oxford, Oxford, UK, e-mail:
<EMAIL_ADDRESS>Watkinson is with the Critical Care Research
Group, Nuffield Department of Clinical Neurosciences, University of Oxford,
Oxford, UKT. Zhu is with the Department of Engineering Science, University of
Oxford, Oxford, UK

###### Abstract

Survival analysis helps approximate underlying distributions of time-to-events
which in the case of critical care like in the ICU can be a powerful tool for
dynamic mortality risk prediction. Extending beyond the classical Cox model,
deep learning techniques have been leveraged over the last years relaxing the
many constraints of their counterparts from statistical methods. In this work,
we propose a novel conditional variational autoencoder-based method called
DySurv which uses a combination of static and time-series measurements from
patient electronic health records in estimating risk of death dynamically in
the ICU. DySurv has been tested on standard benchmarks where it outperforms
most existing methods including other deep learning methods and we evaluate it
on a real-world patient database from MIMIC-IV. The predictive capacity of
DySurv is consistent and the survival estimates remain disentangled across
different datasets supporting the idea that dynamic deep learning models based
on conditional variational inference in multi-task cases can be robust models
for survival analysis.

###### Index Terms:

deep learning, healthcare, personalized medicine, prognostication, survival
analysis, variational autoencoders

##  I Introduction

Survival analysis refers to statistical approaches at estimating distributions
of event times or times it takes for an event to happen as well as rates of
survival over time accounting for censoring. The events in question can be
machine failures in industry or the occurrence of specific diseases and death
[1]. In clinical practice, survival analysis can play a key role and provide
valuable insights into predicting patient outcomes and guiding treatment
decisions [2]. While most traditional applications of methods rely on
statistical models in epidemiology, with the rise of deep learning techniques,
personalised estimation of survival times for individual patients has become
possible [3]. Survival analysis can provide dynamic risk estimation for a
population or an individual patient over a period of time which helps track
the progression of risk and when combined with comparative analysis, variation
in survivability between treatment regimes. In settings like intensive care
units (ICUs), such robust prediction frameworks would be especially useful.
Besides, deep learning is more robust in detecting non-linear patterns in
data, estimating unknown distributions, and learning from complex data
modalities like time-series and images.

The primary objective of survival analysis is to approximate the underlying
distribution of events or hitting times known as the survival function taking
into account the influence of input features [4]. There are diverging
approaches to solving this problem such as assuming a specific distribution
like the Weibull for the survival function and then estimating the parameters
of this distribution using specific models from deep learning [5]. Other
approaches like the standard Cox Proportional Hazards (Cox) model are
considered semi-parametric as they assume a specific base distribution for the
hazard function (related to the survival function) which is then altered
through an exponential factor of the features [6]. In most cases, however, the
features are static and the models suffer under the limitations of the
assumptions made like the proportional hazards assumption in the case of Cox
and some of its derivatives. Using time-series deep learning models can allow
for learning long-term and sequential patterns in patient data without being
restricted to using just static features or just the most recent measurements
while also directly estimating the survival function (or alternatives like the
hazard function) through custom loss implementations.

Our work aims to address these limitations by using components like Long
Short-Term Memory (LSTM) to extract relevant patterns from time-series data
while capturing the long-term dependencies present in a patient’s stay in the
ICU. In addition, instead of assuming a prior distribution for the survival or
hazard function, we directly model the survival through an adapted negative
log-likelihood loss function accounting for censoring in the data. Other deep
learning implementations have used similar components like RNNs but we decide
to use the advantage of variational autoencoders (VAEs) and conditional
extensions which achieve superior performance in latent space generation while
also being useful for generative modelling [7, 8, 9]. We aim to show that
conditional variational inference and autoencoder reconstruction tasks can aid
in learning from complex time-series data by extracting latent features later
used in optimising the survival task. We validate our approach both in a
static and time-varying settings using a combination of benchmark datasets
from survival analysis as well as a recent public ICU database. The ICU
remains a key component of delivering critical care where decisions need to be
made urgently and whose effects can be seen within hours of the stay [10]. In
such a time-sensitive environment, providing dynamic risk estimation through
survival curves and time-dependent prediction can aid in prioritising patients
most urgently in need of care while considering the progression of
deterioration throughout the stay.

Using the fact that VAEs have been established as robust tools for
probabilistic modelling capable of capturing the underlying distribution of
complex high-dimensional data through learning salient low-dimensional
representations, we apply them to model the relationship between covariates
and survival times, enabling accurate estimation of survival probabilities and
event predictions. Since the model is not just multi-modal in input but also
in output with both reconstruction and survival tasks included, conditional
variational autoencoders (CVAE) provide a stable prediction performance [11].
By mapping high-dimensional covariates into a lower-dimensional space, CVAEs
can effectively capture the underlying structure and identify hidden factors
that significantly impact survival outcomes and predict for both binary labels
of mortality as well as times of death. This capability allows deeper insights
into time-resolved mortality risk prediction and disease progression in
critical care.

##  II Materials and Methods

###  II-A Data

Data for survival analysis contains three main sets of variables, the first is
the feature set $\mathcal{X}$ which can consist of static or time-series
features (the latter having measurements at potentially different sampling
frequencies), the time-to-event for the events in question or censoring
respectively, and the outcome label for the event like in standard machine
learning tasks [12]. In all of our dataset implementations, we standardise the
start time for the time-to-events to zero and for MIMIC IV we will detail
further pre-processing steps. In general, the time-to-event values can be left
to be continuous depending on the model being considered but we discretise the
time set into 10 equally spaced time periods for our model in the fashion of
DeepHit and Dynamic-DeepHit [9]. Furthermore, measurements are often right-
censored, meaning that patients can leave the study or be lost to follow-up
and not everyone needs to have experienced the event with their time-to-event
and time-varying features reflecting this. An assumption commonly made
elsewhere and here is that this censoring is unimportant and independent of
the outcome of the study itself [13]. Thus, the dataset can be represented as

$\mathcal{D}=\left\\{\left(\mathcal{X}^{i},\mathcal{t}^{i},y^{i}\right)\right\\}_{i=1}^{N}$
(1)

With $\mathcal{X}$ representing the feature matrix, $\mathcal{t}$ being the
time-to-event as the minimum of the event and censoring time, $y$ being the
label for the outcome, and $N$ samples included. As we will be working with
time-series data in the case of MIMIC-IV, a patient $\mathcal{X}^{i}$ can be
seen as

$\mathcal{X}^{i}=\left\\{\boldsymbol{x}_{1}^{i},\boldsymbol{x}_{2}^{i},\ldots,\boldsymbol{x}_{j}^{i}\right\\}$
(2)

Where $j$ is the length of the time-series for $1\leq j<J^{i}$ where $J$ is
the maximum time step and $\boldsymbol{x}_{j}^{i}$ contains $M$ features
$\left[x_{j,1}^{i},x_{j,2}^{i},\ldots,x_{j,M}^{i}\right]^{\top}$ with
timestamps of measurements
$\left[t_{1}^{i},t_{2}^{i},\ldots,t_{j}^{i}\right]$.

Standard benchmark datasets contain only static features but here we implement
survival analysis on the ICU dataset from MIMIC-IV which contains both static
and time-series features. To show performance across different datasets and
with different sizes, we will succinctly introduce these datasets. Across all
datasets, the event in question is death. The datasets were split into 60%
training, and 20% each for validation and testing. Quantile transformations
have been applied for standardisation and fit only on the training dataset.
Please see appendix for a full list of features for each dataset.

####  II-A1 SUPPORT

The Study to Understand Prognoses and Preferences for Outcomes and Risks of
Treatments contains data from five care academic centres in the United States
for general population survival for the next six months [14]. The result of
the study was a prognostic model to estimate survival for seriously ill
hospitalised patients. The dataset consists of 8,873 samples and 14 features.

####  II-A2 METABRIC

The Molecular Taxonomy of Breast Cancer International Consortium contains
genetic and clinical data from breast cancer patients with 1,904 samples and 9
features [15].

####  II-A3 GBSG

The Rotterdam & German Breast Cancer Study Group contains treatment and
clinical data on 2,232 breast cancer patients with 6 features [16].

####  II-A4 NWTCO

The National Wilm’s Tumor dataset contains staging and clinical data on 4,028
Wilms’ tumor patients with 6 features [17].

####  II-A5 sac3

The simulated dataset contains discrete time event-times with 44 features and
100,000 samples [18].

####  II-A6 sac_admin5

The simulated dataset contains discrete time event-times with 5 features and
50,000 samples [19].

####  II-A7 MIMIC IV

We also implement survival analysis on the de-identified real-world ICU
dataset Medical Information Mart for Intensive Care (MIMIC-IV v. 2.0, July
2022) which includes discharge information for over 15,000 additional ICU
patients compared to the previous release [20]. The dataset contains data from
the Beth Israel Deaconess Medical Center collected between 2008 and 2019. The
dataset contains 71,935 samples of ICU stays with 33 static features
(categorical features were one-hot encoded) and 65 time-varying features.
Static variables include age, sex, unit of admission, and others which did not
have missingness and for the time-series we decided to use forward filling as
clinicians in practice would similarly just look at the last recorded
measurement. If the first set of measurements is missing for some time-varying
features we backward fill for that feature for that patient so all time-series
features across all patients can later be aligned to start from the beginning
of admission and earliest record. Our processing of MIMIC-IV follows from our
previous work but is adapted for the survival scenario with the ICU length of
stay or the maximum time horizon for the event times defined as 10 days and
time-series features taken in 72-hour timesteps [21].

###  II-B Survival Function and Losses

In survival analysis, the main underlying goal is the estimation of the
survival function which represents the probability that no event occurs until
a time t and can be written as

$S(t)=\operatorname{P}(T>t)=\int_{t}^{+\infty}f(u)du=1-F(t)$ (3)

Where $f(t)$ is the probability density function of event time

$f(t)=\lim_{\Delta t\rightarrow 0}\frac{\operatorname{P}(t<T\leq t+\Delta
t)}{\Delta t}$

and $F(t)$ corresponds to the cumulative risk or incidence function
$F(t)=\operatorname{P}(T\leq t)$. $T$ represents the time of the event,
$\operatorname{P}$ is the probability, and $t$ is the specitic timestamp for
risk estimation. The key to training a deep learning model is to learn an
estimate of the cumulative incidence or risk function $\hat{F}(t)$ as the
joint distribution of the event time and outcome label given the observations.
As we discretise the time into intervals, we can estimate this event
probability across arbitrary periods and remain faithful to the original
survival analysis formulation instead of resorting to chained binary
classification [22]. Once we have an estimate of the risk, we can then simply
obtain the survival function and curves by taking its negation from one as
shown in equation (3). We also discretise the time scale in the style of other
popular deep learning approaches so we depend on the probability mass function
instead of the probability density function.

We undertake estimation of the underlying cumulative risk by primarily
optimising for the negative log-likelihood of the joint distribution of the
event time and outcome with right-censoring. For those patients who have
suffered the event, we capture both the outcome and the time at which it
occurs. For censored patients, we capture the censoring time conditioned on
the measurements recorded prior to the censoring. If we assume
$a_{t}=\hat{P}\left(T=t\mid\mathcal{X}\right)$ (the output of the last layer
node of the neural network module) represents the probability of experiencing
the event at time $t$, then the loss can be represented as

$\displaystyle\mathcal{L}_{1}=-$
$\displaystyle\sum_{i=1}^{N}\left[1\left(y^{i}\neq\varnothing\right)\cdot\log\left(\frac{a_{t^{i}}^{i}}{1-\sum_{y\neq\varnothing}\sum_{n\leq
t_{Ji}^{i}}a_{n}^{i}}\right)\right.$ (4)
$\displaystyle\left.+1\left(y^{i}=\varnothing\right)\cdot\log\left(1-\sum_{y\neq\varnothing}\hat{F}\left(t^{i}\mid\mathcal{X}^{i}\right)\right)\right]$

Where $1$ is the indicator function. By optimising for this loss, we estimate
the actual risk distribution for each patient and a prediction can be made for
arbitrary times across event times as outputs of the last layer.

As data are complex (i.e., static and time-varying), we seek to learn a latent
representation using a variational autoencoder that would improve learning for
the task of survival analysis. There are two loss terms, namely the
reconstruction (mean squared error) and the Kullback–Leibler (KL) divergence
in the variational autoencoder. Estimating the distribution of the underlying
latent factors relies on minimising the KL divergence between an approximation
of the true posterior and the true distribution both of which are assumed to
be multivariate Gaussians. As such, learning these distributions means
optimising for their parameters, the mean and standard deviation. The ability
of the decoder to successfully reconstruct the input is captured with a simple
mean squared error term between the reconstruction of the input and the input
itself. Thus, the loss for variational inference can be seen as

$\displaystyle\mathcal{L}_{2}=L_{vae}(E,D)=\left\|\mathcal{X}-\mathcal{X}_{\text{recon
}}\right\|_{2}$ (5) $\displaystyle+\frac{1}{2}\sum_{i=1}^{z_{\text{dim
}}}[\left(\mu_{i}^{2}+\sigma_{i}^{2}\right)-1-\log\left(\sigma_{i}^{2}\right))]$

Where $z$ is the sampled latent vector from the probabilistic encoder for the
learned Gaussian distribution with mean $\mu$ and standard deviation $\sigma$.
For training, however, since sampling is a stochastic process, we use the
reparameterization trick to backpropagate the gradient and represent the
latent vector as the sum of a deterministic variable and an auxiliary
independent random variable $\varepsilon$ [23]

$\varepsilon\sim N(0,1)\quad z=\mu+\varepsilon*\sigma\quad\rightarrow\quad
z\sim N(\mu,\sigma)$ (6)

This latent vector is then used as input to a neural network module in
learning the survival task. During training, the decoder uses the latent
vector and the condition vector (survival labels or times of death) as input
which helps the latent space capture other information instead of trying to
better reconstruct the input. At test time, the decoder is not used anymore,
and the latent space is used for prediction in the survival task. Details on
the network structure implementing the optimisation will be discussed later,
but the total loss can then be presented as

$\mathcal{L}=\alpha\mathcal{L}_{1}+(1-\alpha)\mathcal{L}_{2}$ (7)

Where $\alpha$ is the balancing coefficient between the two losses
($\mathcal{L}_{1}$ and $\mathcal{L}_{2}$) and $0\leq\alpha\leq 1$. $\alpha$ is
considered as a hyperparameter that is optimised during training according to
multi-objective optimization principles.

###  II-C Benchmark Models

To compare DySurv to other survival analysis methods, we implement a selection
of the most popular and consistently cited methods in the field and provide a
short description of each below. We will first introduce discrete-time methods
which rely on discretising the event times into specified durations, and then
follow with continuous-time methods.

####  II-C1 PMF

The parametrisation of the Probability Mass Function (PMF) of the event times
is another way of estimation without resorting to using discrete-time risk or
hazard in likelihood optimisation. We described the continuous probability
density function without using discrete time boundaries. It is the foundation
of other methods like DeepHit and Multi-Task Logistic Regression. It similarly
resorts to optimising a negative log-likelihood loss but instead of using the
cumulative risk function, it uses the approximations of the PMF and survival
functions [18]. Since we can establish a direct representational relation
between the risk and survival functions, the PMF loss can be seen as an
alternative to our own loss function. In our and others’ implementation, the
PMF method is a simple Multi-Layer Perceptron (MLP) optimised for this loss
and we use the same structure as we use for our survival module to be
described later.

####  II-C2 MTLR

Multi-Task Logistic Regression provides a generalization of the binomial log-
likelihood to jointly model the sequence of binary labels for each time
interval risk prediction. This method similarly minimises the negative log-
likelihood with the PMF and survival function terms but the network outputs
are cumulatively summed in reverse to no certain advantage and, in fact, just
adds computational complexity and numeric instability [24, 25].

####  II-C3 BCESurv

This Binary Cross-Entropy for Survival is a method consisting of a set of
binary classifiers that remove individuals as they are censored. The loss is
the binary cross entropy of the survival estimates at a set of discrete times,
with targets that are indicators of surviving each time. Each output node in
the last layer corresponds to a binary classifier evaluated at that time
point. As censored patients are removed, the method is biased towards those
with higher event probabilities [19].

####  II-C4 DeepHit

The single-risk version of DeepHit is a deep learning model whose output nodes
are softmaxed to jointly model the probabilities between the event times and
the time durations are discretised like in our case. The model depends on
optimising both the negative log-likelihood loss based on the cumulative
incidence function and a ranking loss built on the intuition of the
concordance. The ranking loss penalises incorrect ordering of patient pairs in
which the patient that remains longer in the study should have a lower risk at
the end point for the patient with the shorter stay. Including this loss
function allows the model to directly optimise for the concordance metric
which is also the only metric of evaluation used in their paper hence leading
to potentially biased and inflated results. Subsequently reproduced work has
shown that indeed this model is not calibrated well and the inclusion of this
ranking loss, while helping to show better performance as measured by
concordance, significantly lags across other metrics in survival analysis when
compared to simpler models [26, 8, 22, 27].

####  II-C5 Logistic Hazard

The Logistic Hazard method is a submodular implementation of our own deep
learning model using the loss in (4) with an MLP that similarly parametrises
the PMF of the survival times [28]. DySurv expands on this method to include
it as a component in the framework with the variational autoencoder to jointly
optimise for both tasks of reconstruction and latent space formation as well
as survival estimation. A key difference between the log-likelihood loss used
here (and in our model) and in DeepHit is that logistic hazards do not allow
for survival past the maximum time horizon.

####  II-C6 CoxTime

CoxTime is a relative risk model that extends Cox regression beyond the
proportional hazards and is the first of the continuous-time methods. The
standard Cox regression model which we will not spend space introducing here
consists of a baseline hazard term (defined cumulatively in the loss by a pre-
selected estimator such as Breslow) and a relative risk term which is an
exponential factor of the weighted linear combination of features. The basic
model assumes constant proportionality between the patients’ hazards over time
and is thus restrictive. In other words, the difference between survival
likelihoods for a given time is proportional to the difference in feature or
hazard values for patients. CoxTime goes around this assumption by
parametrising the relative risk term as a function of time and not just the
features, thus the non-proportional behaviour over time is modelled by
allowing for time to be considered alongside the features [29].

####  II-C7 CoxCC

CoxCC (Case Control) is just a proportional implementation of the CoxTime
model and is the closest to the standard Cox implementation where the
minimisation of the partial log-likelihood is done with stochastic gradient
descent by averaging over constrained risk sets for mini-batch learning
instead of the entire dataset like in classical survival analysis [29].

####  II-C8 DeepSurv

DeepSurv is a deep learning model that directly minimises the negative partial
log-likelihood as defined in the standard Cox model. It is one of the first
deep learning implementations for survival analysis and the output of the
model is the log-risk term of the Cox model which accounts for nonlinearity
[30]. There is no indirect estimation of cumulative risk or survival through
likelihood estimation like in the previous methods, thus DeepSurv similarly
suffers under limitations of the Cox such as the proportionality assumption
[31].

####  II-C9 PCHazard

The last continuous-time method we introduce and implement is PCHazard which
assumes that the continuous-time hazard rate (instantaneous value of risk) is
piecewise constant. The method relies on optimising for the likelihood
contribution which mimics the MTLR approach albeit in the continuous-time
setting with the hazards parametrised by a simple MLP. The piecewise constant
causes the likelihood to behave like a Poisson likelihood [32]. Despite the
method operating in continuous-time, the hazards are defined in time intervals
which rely on discretisation steps from the observed continuous event times
and censoring times, while we discretise the times to a predefined set of time
flagposts [18].

A summary of all the methods can be seen in Table I.

TABLE I: Survival analysis methods investigated Method | Time scale | Reference
---|---|---
PMF | discrete-time | [18]
MTLR | discrete-time | [24]
BCESurv | discrete-time | [19]
DeepHit | discrete-time | [26]
Logistic Hazard | discrete-time | [28]
CoxTime | continuous-time | [29]
CoxCC | continuous-time | [29]
DeepSurv | continuous-time | [30]
PCHazard | continuous-time | [18]

###  II-D DySurv

Our proposed deep learning model leverages the established structure of a
variational autoencoder and combines it with a survival analysis learning MPL
module to simultaneously optimise for both reconstruction and likelihood
losses. Figure 1 shows the structure of the proposed model. Using a simple
autoencoder has been shown to lead to overfitting and imbalanced learning of
the reconstruction task that could harm learning the survival task whereas a
VAE’s objective function is based on the reconstruction loss from a randomly
sampled vector allowing for more robustness [33]. We concatenate the feature
vectors from the static and time-series features together before feeding them
into an encoder equipped with a Long-Short Term Memory (LSTM) cell. As
mentioned earlier, we use 72 1-hour timesteps and static variables are
included by replicating them before concatenation with the time-series vector.
Indeed, the time-series components are compressed into a latent representation
which is then used, as is, as input to the estimation of the negative log-
likelihood loss function. Following the LSTM unit, the remaining part of the
encoder consists of an MLP module with 3 layers. The encoder and decoder are
mirrored in their structure with the hidden neurons in the MLP layers
consisting of 3 times feature length, then 5 times, then 3 times, before
passing onto the last output layer. For the encoder, the last layer is used in
approximating the mean and standard deviation of the latent Gaussian
distribution. As for the decoder, the last output layer is used for
reconstructing the input vector. Once the latent vector is sampled from the
Gaussian distribution defined by these parameters, the lower-dimensional
latent factors are used as input for an MLP module in optimising the survival
task. The survival module is similar to the MLP components of the encoder and
decoder. All components are jointly optimised through multi-loss optimisation.
All of the components have been investigated with and without dropout
included. The output of the survival module is 10 nodes softmaxed hence
properly jointly distributed, each of which gives the probability that the
patient has suffered the event (death) at that specific time interval. In the
MIMIC-IV case, each node then corresponds to 1-day risk prediction in the ICU
as the maximum time horizon is 10 days. The hyperparameters optimised in the
network through grid search using the training and validation set included
learning rate, batch size, $\alpha$, and dropout proportion. To minimise
overfitting, we employ early stopping techniques.

Figure 1: Proposed DySurv model for dynamic risk prediction based on a
combination of static and time-series EHR data

###  II-E Metrics

In this section, we will switch the notation of samples from superscript to
subscript, hence $\mathbf{x}^{i}$ is now $\mathbf{x}_{i}$ for sample $i$.
Since we are no longer making single risk predictions at specific times and
estimating distributions of event times for censored samples, different
evaluation metrics must apply than those used in classic machine learning
classification and prediction settings. The most common metric for evaluating
survival analysis models is the concordance index $C_{ind}$, which estimates
the probability that, for a random pair of samples, the predicted survival
times (risk probabilities) of the two samples have the same ordering as their
true survival times [29]. This explanation works perfectly for settings of
proportional hazards where the ordering does not change over time but for our
purposes, we will not be limited by such assumptions. Hence, we will rely on
using the time-dependent extension $C_{ind}^{\mathrm{td}}$ with some
modifications accounting for predictions independent of feature observations
having a concordance of 0.5. The metric can be represented as

$C_{ind}^{\mathrm{td}}=\mathrm{P}\left\\{\hat{S}\left(T_{i}\mid\mathbf{x}_{i}\right)<\hat{S}\left(T_{i}\mid\mathbf{x}_{j}\right)\mid
T_{i}<T_{j},D_{i}=1\right\\}$ (8)

Where $\hat{S}$ indicates the estimated survival probabilities are used and
$D_{i}=1$ that only those who experienced the event are considered in this
metric. A noted limitation of this metric is its obvious bias and dependence
on the censoring distribution as only non-censored samples are considered
making it affected by the length of stay and the censoring proportion that
increases over the length of stay. To this end, we decide to use additional
metrics for more holistic evaluation especially as previously proposed models
like DeepHit were found to be ungeneralisable when evaluated using other
metrics besides concordance. We also evaluate our model using the Integrated
Brier Score or IBS. The Brier Score is similar to the mean squared error as it
represents the average squared distances between the predicted and the true
survival probability (approximated with step functions with jumps at the event
times) and is always a number between 0 and 1, with 0 being the best possible
value [25]. In fact, the expectation of the Brier Score contains the mean
squared error as one of its additive terms, so minimising one is minimising
the other [19]. Since we need to know the event times for calculating IBS and
we do not have access to all the samples’ event times in right-censoring, an
adjusted metric called the inverse probability of censoring weights Brier
Score (IPCW) is used instead to approximate the times by weighting the scores
of the observed event times by the inverse probability of censoring. The
equation used is thus

$\displaystyle\operatorname{BS}_{\operatorname{IPCW}}(t)=\frac{1}{n}\sum_{i=1}^{n}[\frac{\hat{S}_{i}(t)^{2}1\left\\{T_{i}\leq
t,D_{i}=1\right\\}}{\hat{G}_{i}\left(T_{i}-\right)}+$ (9)
$\displaystyle\frac{\left[1-\hat{S}_{i}(t)\right]^{2}1\left\\{T_{i}>t\right\\}}{\hat{G}_{i}(t)}]$

Where $\hat{G}_{i}(t)=\mathrm{P}\left(C_{i}^{*}>t\right)>0$ is the Kaplan-
Meier estimate of the censoring distribution for sample $i$ and $C_{i}^{*}$ is
the censoring time. The expected value of this metric is the same as for the
uncensored Brier Score. As one notices, this metric is evaluated at specific
times whereas the Integrated Brier Score or IBS provides a general evaluation
of model performance at all times

$\mathrm{IBS}=\frac{1}{\max\left(T_{i}\right)}\int_{0}^{\max\left(T_{i}\right)}BS_{\operatorname{IPCW}}(t)dt$
(10)

But a limitation of this metric is the biased assumption of the censoring
distribution estimated by the nonparametric Kaplan-Meier estimator which
disregards the features, meaning all samples are assumed to have the same
censoring distribution. This can be addressed by using an administrative
extension of the metric that requires access to all the censoring times but a
discussion of this is left to the reader to peruse as they see fit [19].
Lastly, we introduce the IPCW (negative) binomial log-likelihood or NBLL from
classic binary classification which measures both discrimination and
calibration of the estimates and use its integrated extension for all times
INBLL

$\displaystyle\operatorname{BLL}(t)=\frac{1}{N}\sum_{i=1}^{N}[\frac{\log\left[1-\hat{S}\left(t\mid\mathbf{x}_{i}\right)\right]1\left\\{T_{i}\leq
t,D_{i}=1\right\\}}{\hat{G}\left(T_{i}\right)}+$ (11)
$\displaystyle\frac{\log\left[\hat{S}\left(t\mid\mathbf{x}_{i}\right)\right]1\left\\{T_{i}>t\right\\}}{\hat{G}(t)}]$

$\mathrm{IBLL}=\frac{1}{\max\left(T_{i}\right)}\int_{0}^{\max\left(T_{i}\right)}\operatorname{BLL}(t)dt$
(12)

For both of the last metrics, we approximate the integrals by numerical
integration (for 100 timesteps as based on previous literature), and the time
span is the duration of the test set as these metrics are only evaluated on
the test set [29].

##  III Results

To holistically evaluate DySurv, we present a set of experiments and
comparisons with other benchmark survival analysis models across multiple
datasets. We not only present the discriminative performance of the model as
measured by concordance, but also its calibration as measured by IBS and IBLL.
As the main aim of the model is to issue dynamic survival scores for patients
in the ICU, inclusion of the real-world MIMIC-IV electronic health records
dataset provides the paper with practical significance in addition to a
methodological contribution. The results consist of two major experiments, one
is the ability of the model to successfully learn from static data which is
present in all the datasets, and the other to learn from a combination of
static and time-varying data such as MIMIC-IV. For these purposes, Dynamic-
DeepHit is the only relevant comparison as other survival analysis models deal
with static data only. Tables II and III show these results across datasets
and for all models included. We implemented these methods in PyTorch (PyCox)
v. 1.10.1 using a MacBook Air M1 2021 laptop with data processing completed
using pandas and SQL.

TABLE II: Test results on different datasets for survival analysis models and DySurv as evaluated by three different metrics introduced in Materials and Methods. For concordance, higher is better, and for the other two metrics, lower is better. All of the results are an average of five random seeds. | $C_{ind}^{\mathrm{td}}$ | IBS | IBLL |  | $C_{ind}^{\mathrm{td}}$ | IBS | IBLL
---|---|---|---|---|---|---|---
SUPPORT |  |  |  | METABRIC |  |  |
PMF | 57.9 | 0.195 | 0.574 | PMF | 63.8 | 0.168 | 0.497
MTLR | 55.3 | 0.205 | 0.775 | MTLR | 56.8 | 0.172 | 0.527
BCESurv | 55.3 | 0.290 | 2.08 | BCESurv | 56.8 | 0.138 | 0.477
DeepHit | 57.3 | 0.273 | 0.678 | DeepHit | 65.5 | 0.123 | 0.415
Logistic Hazard | 53.5 | 0.206 | 0.762 | Logistic Hazard | 59.0 | 0.163 | 0.498
CoxTime | 59.5 | 0.193 | 0.565 | CoxTime | 65.4 | 0.114 | 0.361
CoxCC | 59.7 | 0.192 | 0.563 | CoxCC | 65.9 | 0.166 | 0.508
DeepSurv | 60.6 | 0.190 | 0.559 | DeepSurv | 62.4 | 0.176 | 0.541
PCHazard | 55.1 | 0.206 | 0.633 | PCHazard | 51.4 | 0.160 | 0.547
DySurv | 64.7 | 0.190 | 0.561 | DySurv | 64.5 | 0.120 | 0.387
GBSG |  |  |  | NWTCO |  |  |
PMF | 68.5 | 0.179 | 0.528 | PMF | 69.7 | 0.122 | 0.389
MTLR | 65.6 | 0.180 | 0.542 | MTLR | 66.8 | 0.109 | 0.403
BCESurv | 65.6 | 0.156 | 0.481 | BCESurv | 69.1 | 0.108 | 0.393
DeepHit | 68.1 | 0.174 | 0.514 | DeepHit | 71.1 | 0.118 | 0.348
Logistic Hazard | 67.4 | 0.179 | 0.537 | Logistic Hazard | 66.5 | 0.108 | 0.396
CoxTime | 68.4 | 0.171 | 0.510 | CoxTime | 70.7 | 0.110 | 0.343
CoxCC | 59.6 | 0.205 | 0.597 | CoxCC | 70.3 | 0.110 | 0.373
DeepSurv | 68.5 | 0.180 | 0.531 | DeepSurv | 68.3 | 0.115 | 0.391
PCHazard | 55.8 | 0.182 | 0.574 | PCHazard | 60.2 | 0.118 | 0.465
DySurv | 70.4 | 0.164 | 0.499 | DySurv | 70.3 | 0.111 | 0.347
sac3 |  |  |  | sac_admin5 |  |  |
PMF | 74.3 | 0.125 | 0.391 | PMF | 71.5 | 0.124 | 0.387
MTLR | 65.0 | 0.124 | 0.539 | MTLR | 65.7 | 0.122 | 0.520
BCESurv | 67.8 | 0.163 | 0.586 | BCESurv | 68.4 | 0.164 | 0.505
DeepHit | 74.2 | 0.184 | 0.527 | DeepHit | 71.6 | 0.186 | 0.396
Logistic Hazard | 72.0 | 0.120 | 0.492 | Logistic Hazard | 70.7 | 0.118 | 0.481
CoxTime | 78.7 | 0.117 | 0.362 | CoxTime | 78.5 | 0.117 | 0.362
CoxCC | 76.4 | 0.124 | 0.384 | CoxCC | 76.7 | 0.122 | 0.381
DeepSurv | 76.1 | 0.126 | 0.390 | DeepSurv | 77.4 | 0.119 | 0.371
PCHazard | 64.0 | 0.135 | 0.514 | PCHazard | 65.1 | 0.123 | 0.503
DySurv | 80.6 | 0.112 | 0.359 | DySurv | 79.6 | 0.116 | 0.361
TABLE III: Test results on MIMIC-IV dataset for survival analysis models and DySurv as evaluated by three different metrics introduced in Materials and Methods. For concordance, higher is better, and for the other two metrics, lower is better. All of the results are an average of five random seeds. | $C_{ind}^{\mathrm{td}}$ | IBS | IBLL
---|---|---|---
MIMIC-IV |  |  |
PMF | 50.9 | 0.126 | 0.389
MTLR | 52.4 | 0.126 | 0.389
BCESurv | 52.2 | 0.157 | 0.473
DeepHit | 54.4 | 0.137 | 0.421
Logistic Hazard | 52.6 | 0.122 | 0.396
CoxTime | 53.1 | 0.122 | 0.337
CoxCC | 52.9 | 0.123 | 0.393
DeepSurv | 54.2 | 0.128 | 0.403
PCHazard | 51.0 | 0.122 | 0.378
DySurv (static) | 55.7 | 0.111 | 0.360
Dynamic-DeepHit | 56.0 | 0.143 | 0.376
DySurv ($+$ time-series) | 57.9 | 0.122 | 0.320

In survival analysis, it is also vital to show the performance of the models
in creating survival curves, or estimates of the survival probabilities over
time for different patients. Whereas traditional survival analysis models from
statistics rely on risk sets and computing survival estimates for a
population, an advantage of deep learning models is that the task of survival
estimation can apply to an individual patient sample. We provide survival
curves for a random group of five samples/patients across all datasets that
show a clear separation of risk (entanglement indicates the model has not
successfully learned the risk trajectories) for different samples as generated
by DySurv. Figure 2 indicates these for all the datasets except MIMIC-IV
considered.

(a) Survival curve (estimate of survival probability over time) for five
random samples from SUPPORT

(b) Survival curve (estimate of survival probability over time) for five
random samples from METABRIC

(c) Survival curve (estimate of survival probability over time) for five
random samples from GBSG

(d) Survival curve (estimate of survival probability over time) for five
random samples from NWTCO

(e) Survival curve (estimate of survival probability over time) for five
random samples from sac3

(f) Survival curve (estimate of survival probability over time) for five
random samples from sac5

Figure 2: The set of survival curves for benchmark datasets as generated by
DySurv that shows clear extrapolation of risk across different samples. Since
DySurv provides discrete estimates over time-intervals, the curves were
generated with additional interpolation commonly used in other discrete
methods in survival analysis.

For MIMIC-IV, the data input consists of both static and time-series data.
Here we provide survival curve results for both cases as well as for Dynamic-
DeepHit in the time-series case and show how DySurv has greater capacity at
separating risk trajectories in both cases as Figure 3 shows.

(a) Survival curve (estimate of survival probability over time) for five
random samples from MIMIC-IV using only static features with DySurv

(b) Survival curve (estimate of survival probability over time) for five
random samples from MIMIC-IV using both static and time-series features with
DySurv

(c) Survival curve (estimate of survival probability over time) for five
random samples from MIMIC-IV using both static and time-series features with
Dynamic-DeepHit

Figure 3: The set of survival curves for MIMIC-IV ICU EHR dataset as generated
by DySurv that shows clear extrapolation of risk across different patients as
compared to both static, time-series, and Dynamic-DeepHit cases. The same set
of patients was used in generating the survival curves across cases.

Upon deployment of a trained model, DySurv can generate risk estimates through
time for each patient while using their history of observations. We do not
rely on using landmarking methods or a specific pre-defined time for risk
prediction as the score are issued across the entire time interval
simultaneously. For data pre-processing purposes, a time horizon is selected
corresponding to 24 hours with 72-hour timespan for the LSTM in the case of
MIMIC-IV since ICU risk assessments often use information over 72 hours for
the next 24-hour risk prediction [34].

##  IV Discussion

The first set of results relates to applying DySurv only on static data from
several benchmark datasets of varying sizes. We see that for the vast majority
of these benchmarks, DySurv outperforms both standard statistical as well as
deep learning alternatives across all metrics except for METABRIC and NWTCO on
the concordance where DeepHit tends to perform slightly better. This is
probably due to the implementation of the biased ranking loss mentioned
earlier that aids in having better discriminative performance as measured by
the concordance metric but that is not reflected as measured by the other two
metrics. Similar behaviour for DeepHit has been oberved in another study by
[29]. We also see that the non-VAE implementation of the logistic hazard
performs much worse than DySurv across all experiments, thereby strengthening
the idea that adding variational inference to the logistic hazard can aid in
learning the survival task. This improvement occurs despite having the
additional task of reconstruction now. The predictive advantage comes from the
identification of lower dimensional latent vectors used in the survival task
instead of the raw features directly. Furthermore, on very large synthetic
datasets, such as sac3 and sac_admin5, DySurv performs better due to having a
greater amount of data to learn from. A limitation of the other benchmark
datasets is their relatively small size may constrain DySurv from learning its
optimal parameters to provide netter survival prediction.

While there are previous models that have attempted to use autoencoders for
survival analysis such as [33] and [35], they have not explored variational
inference extensively. These models also rely on optimisation of the Cox
partial log-likelihood loss, hence being restricted by the proportionality
assumption and they do not account for dynamic time-series or time-varying
features in the input. Furthermore, work by [33] suggests that the VAE model’s
learned compact latent representation directly aids in the Cox model’s
improved performance. This intuition is precisely what we have also seen in
our results albeit on a larger, more flexible, and more complex scenario of
time-series ICU risk prediction with direct joint distribution estimation. The
concatenation autoencoder from [35] is not even compatible with static data,
by far the most common modality in survival analysis, thereby limiting its
relevance significantly. DySurv addresses all these limitations and provides a
flexible solution to dynamic survival analysis with deep learning.

As for MIMIC-IV results, we see that even with only static data in the input,
DySurv manages to outperform other survival analysis models. When time-series
data is added to the input in multi-modal fashion, it can help the model
improve its performance and outperform both CoxTime and Dynamic-DeepHit across
all metrics. Survival curves in Figure 3 show clearly that the static version
of DySurv has the clearest separation of survival trajectories for different
patients and that adding time-series data makes the task more complex. This is
probably due to the reconstruction task being a lot more difficult now that
time-series are involved and the model makes sacrifices on the survival task
front. Nonetheless, compared to Dynamic-DeepHit whose survival curves are
barely disentangled and hence unusable, the time-series version of DySurv
still provides relatively reliable survival estimation. We also see that, as
we would expect clinically, the survival of the patients significantly changes
in the last few days in the ICU, starting to drop a few days before death.
Previous work has shown that earlier times in the ICU correspond with higher
survival rates. This suggests that the identification of the time period when
survival rates drop dramatically in the ICU can help with targetting earlier
treatment for those most at risk [36]. Patient deterioration of this kind in
the ICU can help aid physicians in emergency medicine in identifying those
patients most at risk in the future and proactively addressing their needs.

In this paper, we present a novel dynamic risk prediction model for survival
in the ICU based on deep learning in survival analysis paradigms. Our method
builds on a combination of previous work including Dynamic-DeepHit by
leveraging direct learning of the joint distribution of the first event time
and the event through log-likelihood optimisation with logistic hazards.
Theoretically, this approach is an alternative to the risk log-likelihood loss
function of Dynamic-DeepHit itself that does not use a ranking loss for biased
inflation of concordance results. Our DySurv model is capable of learning from
complex EHR ICU time-series data and extracting lower-dimensional latent
representations that can be useful for learning the survival task while also
balancing reconstruction. As the model has been difficult to train due to loss
instabilities and sensitivity to hyperparameter selection, future work can
explore including a regularisation component to the loss terms. And finally,
using the underlying latent distribution to directly model an alternative of
the survival distribution like the Weibull distribution instead of using a
Gaussian intermediate.

##  V Contributions

MM and TZ conceived and designed the study. The data was curated by MM. Formal
analysis was undertaken by MM. Development of the statistical analysis and
machine learning methodologies was completed by MM. Supervision was provided
by PW and TZ. Visualisations, writing, and editing was done by MM. The
corresponding author and TZ had full access to all data.

##  VI Data sharing

The data used in the study is all publicly available. MIMIC-IV can be accessed
after some initial requests and approvals from the PhysioNet website. The
remaining datasets have their sources listed in the PyCox library repository.

##  VII Acknowledgements

## References

  * [1] Lee MLT, Whitmore GA.  Threshold regression for survival analysis: modeling event times by a stochastic process reaching a boundary.  Statistical Science. 2006.
  * [2] Yoon J, Alaa A, Cadeiras M, Van Der Schaar M.  Personalized donor-recipient matching for organ transplantation.  In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 31; 2017. .
  * [3] Luck M, Sylvain T, Cardinal H, Lodi A, Bengio Y.  Deep learning for patient-specific kidney graft survival analysis.  arXiv preprint arXiv:170510245. 2017.
  * [4] Singh R, Mukhopadhyay K.  Survival analysis in clinical trials: Basics and must know areas.  Perspectives in clinical research. 2011;2(4):145.
  * [5] Liu X.  Survival analysis: models and applications.  John Wiley & Sons; 2012.
  * [6] Lin H, Zelterman D. Modeling survival data: extending the Cox model. Taylor & Francis; 2002.
  * [7] Giunchiglia E, Nemchenko A, van der Schaar M.  Rnn-surv: A deep recurrent model for survival analysis.  In: Artificial Neural Networks and Machine Learning–ICANN 2018: 27th International Conference on Artificial Neural Networks, Rhodes, Greece, October 4-7, 2018, Proceedings, Part III 27. Springer; 2018. p. 23-32.
  * [8] Ren K, Qin J, Zheng L, Yang Z, Zhang W, Qiu L, et al.  Deep recurrent survival analysis.  In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 33; 2019. p. 4798-805.
  * [9] Lee C, Yoon J, Van Der Schaar M.  Dynamic-deephit: A deep learning approach for dynamic survival analysis with competing risks based on longitudinal data.  IEEE Transactions on Biomedical Engineering. 2019;67(1):122-33.
  * [10] Dummitt B, Zeringue A, Palagiri A, Veremakis C, Burch B, Yount B.  Using survival analysis to predict septic shock onset in ICU patients.  Journal of Critical Care. 2018;48:339-44.
  * [11] Itkina M, Ivanovic B, Senanayake R, Kochenderfer MJ, Pavone M.  Evidential sparsification of multimodal latent spaces in conditional variational autoencoders.  Advances in Neural Information Processing Systems. 2020;33:10235-46.
  * [12] Bradburn MJ, Clark TG, Love SB, Altman DG.  Survival analysis part II: multivariate data analysis–an introduction to concepts and methods.  British journal of cancer. 2003;89(3):431-6.
  * [13] Tsiatis AA, Davidian M.  Joint modeling of longitudinal and time-to-event data: an overview.  Statistica Sinica. 2004:809-34.
  * [14] Knaus WA, Harrell FE, Lynn J, Goldman L, Phillips RS, Connors AF, et al.  The SUPPORT prognostic model: Objective estimates of survival for seriously ill hospitalized adults.  Annals of internal medicine. 1995;122(3):191-203.
  * [15] Curtis C, Shah SP, Chin SF, Turashvili G, Rueda OM, Dunning MJ, et al.  The genomic and transcriptomic architecture of 2,000 breast tumours reveals novel subgroups.  Nature. 2012;486(7403):346-52.
  * [16] Schumacher M, Bastert G, Bojar H, Hübner K, Olschewski M, Sauerbrei W, et al.  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. 1994;12(10):2086-93.
  * [17] Breslow NE, Chatterjee N.  Design and analysis of two-phase studies with binary outcome applied to Wilms tumour prognosis.  Journal of the Royal Statistical Society: Series C (Applied Statistics). 1999;48(4):457-68.
  * [18] Kvamme H, Borgan Ø.  Continuous and discrete-time survival prediction with neural networks.  arXiv preprint arXiv:191006724. 2019.
  * [19] Kvamme H, Borgan Ø.  The brier score under administrative censoring: Problems and solutions.  arXiv preprint arXiv:191208581. 2019.
  * [20] Johnson A, Bulgarelli L, Pollard T, Horng S, Celi LA, Mark R.  Mimic-iv.  version 04) PhysioNet https://doi org/1013026/a3wn-hq05. 2020.
  * [21] Mesinovic M, Watkinson P, Zhu T.  XMI-ICU: Explainable Machine Learning Model for Pseudo-Dynamic Prediction of Mortality in the ICU for Heart Attack Patients.  arXiv preprint arXiv:230506109. 2023.
  * [22] Sun Z, Dong W, Shi J, He K, Huang Z.  Attention-based deep recurrent model for survival prediction.  ACM Transactions on Computing for Healthcare. 2021;2(4):1-18.
  * [23] Kingma DP, Welling M.  Auto-encoding variational bayes.  arXiv preprint arXiv:13126114. 2013.
  * [24] Yu CN, Greiner R, Lin HC, Baracos V.  Learning patient-specific cancer survival distributions as a sequence of dependent regressors.  Advances in neural information processing systems. 2011;24.
  * [25] Fotso S.  Deep neural networks for survival analysis based on a multi-task framework.  arXiv preprint arXiv:180105512. 2018.
  * [26] Lee C, Zame W, Yoon J, Van Der Schaar M.  Deephit: A deep learning approach to survival analysis with competing risks.  In: Proceedings of the AAAI conference on artificial intelligence. vol. 32; 2018. .
  * [27] Nagpal C, Li X, Dubrawski A.  Deep survival machines: Fully parametric survival regression and representation learning for censored data with competing risks.  IEEE Journal of Biomedical and Health Informatics. 2021;25(8):3163-75.
  * [28] Gensheimer MF, Narasimhan B.  A scalable discrete-time survival model for neural networks.  PeerJ. 2019;7:e6257.
  * [29] Kvamme H, Borgan Ø, Scheel I.  Time-to-event prediction with neural networks and Cox regression.  arXiv preprint arXiv:190700825. 2019.
  * [30] Katzman JL, Shaham U, Cloninger A, Bates J, Jiang T, Kluger Y.  DeepSurv: personalized treatment recommender system using a Cox proportional hazards deep neural network.  BMC medical research methodology. 2018;18(1):1-12.
  * [31] Thammasorn P, Schaub SK, Hippe DS, Spraker MB, Peeken JC, Wootton LS, et al.  Regularizing the Deepsurv network using projection loss for medical risk assessment.  IEEE Access. 2022;10:8005-20.
  * [32] Friedman M.  Piecewise exponential models for survival data with covariates.  The Annals of Statistics. 1982;10(1):101-13.
  * [33] Kim S, Kim K, Choe J, Lee I, Kang J.  Improved survival analysis by learning shared genomic information from pan-cancer data.  Bioinformatics. 2020;36(Supplement_1):i389-98.
  * [34] Yu S, Leung S, Heo M, Soto GJ, Shah RT, Gunda S, et al.  Comparison of risk prediction scoring systems for ward patients: a retrospective nested case-control study.  Critical Care. 2014;18:1-9.
  * [35] Tong L, Mitchel J, Chatlin K, Wang MD.  Deep learning based feature-level integration of multi-omics data for breast cancer patients survival analysis.  BMC medical informatics and decision making. 2020;20:1-12.
  * [36] Simchen E, Sprung CL, Galai N, Zitser-Gurevich Y, Bar-Lavi Y, Gurman G, et al.  Survival of critically ill patients hospitalized in and out of intensive care units under paucity of intensive care unit beds.  Critical care medicine. 2004;32(8):1654-61.

##  Appendix A Appendix

###  A-A Data Description

A detailed list of features used in the study from different datasets can be
seen in Tables IV and V.

TABLE IV: Features extracted from the bechmark survival datasets. Feature | Type | Feature | Type
---|---|---|---
Sex | binary | NWTCO |
Age | integer | Age | integer
Number of comorbidities | integer | Relapse | binary
SUPPORT Coma Score | continuous | Histology Local | binary
Years of Education | continuous | Histology Institute | binary
APS III Score | continuous | Tumor Stage | categorical
Diabetes | binary | Study Type | categorical
Dementia | binary | GBSG |
3rd Day MAP | continuous | Age | integer
3rd Day WBC | continuous | Hormone Treatment | binary
3rd Day HR | continuous | Menopause Treatment | binary
3rd Day RR | continuous | Tumor Grade | categorical
Temperature | continuous | Number of positive nodes | categorical
3rd Day Bilirubin | continuous | Progesterone (fmol) | continuous
METABRIC |  | Estrogen (fmol) | continuous
Age | integer |  |
Gene I | binary |  |
Gene II | binary |  |
Gene III | binary |  |
Gene IV | binary |  |
Hormone Treatment | binary |  |
Radiotherapy | binary |  |
Chemotherapy | binary |  |
ER-positive | binary |  |
TABLE V: Features extracted from the MIMIC-IV database. The features include demographic data collected for all patients, ICU unit-specific information like type of unit, hospital information like regional location, time since admission, vital signs including respiratory rate and blood pressure, and biochemical measurements including blood glucose and hemoglobin. Feature | Type | Feature | Type
---|---|---|---
Sex | binary | Braden Score | continuous
Age | integer | Strength L Arm | continuous
Height | continuous | Strength R Arm | continuous
Weight | continuous | Strength L Leg | continuous
Hour of Admission | integer | Strength R Leg | continuous
Time Since Admission | continuous | GCS - Eye | continuous
Eye Response | continuous | GCS - Motor | continuous
Motor Response | continuous | GCS - Verbal | continuous
Verbal Response | continuous | Daily Weight | continuous
Ethnicity | categorical | ALT | continuous
Unit Type | categorical | AST | continuous
Admission Location | categorical | HCO3 | continuous
Insurance | categorical | Hct | continuous
Time-series (summary features) |  |  |
ALT | continuous | Alkaline Phosphatase | continuous
Anion Gap | continuous | AST | continuous
Base Excess | continuous | Bicarbonate | continuous
Bilirubin | continuous | Calcium | continuous
Total CO2 | continuous | Chloride | continuous
Creatinine | continuous | Glucose | continuous
Hematocrit | continuous | Hemoglobin | continuous
INR(PT) | continuous | Lactate | continuous
MCH | continuous | MCHC | continuous
MCV | continuous | Magnesium | continuous
PT | continuous | PTT | continuous
Phosphate | continuous | Platelet Count | continuous
Potassium | continuous | RDW | continuous
Red Blood Cells | continuous | Sodium | continuous
Urea Nitrogen | continuous | White Blood Cells | continuous
pCO2 | continuous | pH | continuous
pO2 | continuous | JH-HLM | continuous
Dyspnea Assessment | continuous | Daily Weight | continuous
Glucose | continuous | Heart Rate | continuous
DBP | continuous | SBP | continuous
O2 Flow | continuous | O2 Sat (%) | continuous
Pain Level | continuous | Pain Level Response | continuous
Phosphorous | continuous | Respiratory Rate | continuous
Richmond-RAS Scale | continuous |  Temperature (∘F) | continuous


# Reward Shaping for Building Trustworthy Robots in Sequential Human-Robot
Interaction

Yaohui Guo1 and X. Jessie Yang1 and Cong Shi1 *This work was supported by the
Rackham pre-doctoral fellowship awarded to the first author and the Air Force
Office of Scientific Research under Grant No. FA9550-20-1-0406.1Yaohui Guo, X.
Jessie Yang and Cong Shi are with the Department of Industrial and Operations
Engineering, University of Michigan, 1209 Beal Avenue, Ann Arbor, Michigan
yaohuig, xijyang<EMAIL_ADDRESS>
###### Abstract

Trust-aware human-robot interaction (HRI) has received increasing research
attention, as trust has been shown to be a crucial factor for effective HRI.
Research in trust-aware HRI discovered a dilemma — maximizing task rewards
often leads to decreased human trust, while maximizing human trust would
compromise task performance. In this work, we address this dilemma by
formulating the HRI process as a two-player Markov game and utilizing the
reward-shaping technique to improve human trust while limiting performance
loss. Specifically, we show that when the shaping reward is potential-based,
the performance loss can be bounded by the potential functions evaluated at
the final states of the Markov game. We apply the proposed framework to the
experience-based trust model, resulting in a linear program that can be
efficiently solved and deployed in real-world applications. We evaluate the
proposed framework in a simulation scenario where a human-robot team performs
a search-and-rescue mission. The results demonstrate that the proposed
framework successfully modifies the robot’s optimal policy, enabling it to
increase human trust at a minimal task performance cost.

##  I Introduction

Human trust plays a crucial role in human-robot interaction (HRI) as it
mediates the human’s reliance on the robot, thus directly affecting the
effectivenes of the human-robot team [1, 2, 3]. As a result, researchers have
proposed trust-aware human-robot planning [4], which equips a robot with the
ability to estimate and anticipate human trust and enables it to strategically
plan its actions to foster better cooperation, improve teamwork, and
ultimately enhance the overall performance of the human-robot team.

Trust-aware HRI explicitly consider human trust in the robot’s decision-making
processes [4, 5, 6]. Chen et al. [4] modeled the sequential HRI process as a
partially observable Markov decision process (POMDP) and incorporated human
trust into the state transition function, allowing the robot to optimize its
objectives in the interaction process. However, the authors discovered a
dilemma — maximizing task rewards often leads to decreased human trust, while
maximizing human trust would compromise task performance. Guo et al. further
investigated the dilemma and showed that the robot could intentionally
“deceive” its human partner in order to maximize its rewards [5]. To address
this issue, they proposed adding a “trust-seeking” term to the reward function
to encourage the robot to increase human trust and avoid deception. However,
it remains unclear how to design such terms.

Figure 1: This figure illustrates the idea of improving human trust through
reward shaping. The horizontal axis is the total task reward while the
vertical axis is the total trust-related reward. The shaded area stands for
the collection of policies of our interest. A point
$\mathbf{V}^{\pi}\in\mathbb{R}^{2}$ indicates the value a policy $\pi$ can
earn during the interaction process. The policy $\pi^{*}$ is the optimal
policy for earning the task reward as its value attains the maximum along the
total task reward axis; policy $\pi^{t*}$ is the optimal policy for gaining
human trust. The bold line is the Pareto front, i.e., the set of policies that
yield the best possible trade-offs between human trust and total task rewards.
Our goal is to find a policy $\pi^{\prime}$ such that it seeks to improve
human trust at a small cost of task reward.

The problem of balancing the robot’s total task reward and human trust is
fundamentally a bi-objective optimization problem, as illustrated in figure 1.
It has been shown in previous studies that there does not always exist a
unique dominant policy for the robot – a policy that can maximize the robot’s
total reward and human trust at the same time [5, 4]. A prevalent concept in
bi-objective optimization involves determining the Pareto front, which
represents an optimal balance between conflicting objectives, thereby allowing
the system designer to select a suitable policy from this set. However,
computing the Pareto front in an MDP can be computationally heavy, and it
requires expert knowledge to choose an optimal one from the Pareto front.

Instead, we address the reward design problem in trust-aware sequential HRI by
reward-shaping, a method used in reinforcement learning to guide the learning
process of an agent by modifying its reward function. We model the sequential
HRI process as a two-player Markov game and aim to design a shaping reward
that can guide the robot to gain human trust while guaranteeing small
performance loss. Specifically, we prove that if the shaping reward is a
carefully designed potential-based function, the performance loss can be
bounded by the potential function evaluated at the final states of the Markov
game. We show that applying the proposed framework to the experience-based
trust model results in a linear program that can be efficiently solved. To
evaluate the proposed method, we simulated a scenario where a human-robot team
performs a search-and-rescue (SAR) mission. The results demonstrate that our
method successfully modifies the robot’s optimal policy, enabling it to
increase human trust while satisfying the constraint on performance loss.

The contributions of this work include:

  * •

Proposing a computational framework for balancing task performance and human
trust in trust-aware sequential HRI.

  * •

Developing a novel reward-shaping method for designing the reward term with a
theoretical guarantee on the performance loss.

The remainder of this paper is organized as follows: we first review relevant
literature in section II. In section III, we formulate the trust-aware human-
robot interaction problem as a two-player Markov game. In section IV, we
present the reward-shaping method for designing the trust-seeking term. In
section V, we present the simulation scenario used to evaluate the proposed
framework. Finally, in section VI, we discuss the results and summarize the
work.

##  II Related Work

In this work, we investigate the problem of reward design for trustworthy
robots in trust-aware sequential HRI. Our method is developed upon previous
research, including studies on human trust in robots and trust-aware decision-
making, as well as reward shaping. We review the relevant literature in this
section.

It is worth noting that trust has been widely studied in different areas and
given different definitions. In particular, in this work, we use the
definition given by Lee and See [1], which highlights the uncertainty in HRI:
“trust is the attitude that an agent will help achieve an individual’s goals
in a situation characterized by uncertainty and vulnerability”.

###  II-A Computational Trust Model in HRI

Previous literature attempted to understand the development and evolution of
human trust in an HRI process. Empirical studies have been conducted to
understand the temporal dynamics of trust when a person interacts with
autonomy repeatedly [7, 8, 9, 10, 11, 12, 13]. Three major properties that
characterize how a person’s trust in autonomy changes due to moment-to-moment
interactions with autonomy are summarized in [8, 14], namely continuity,
negativity bias, and stabilization.

In addition, several computational trust models have been developed, including
[8, 15, 16, 17, 18]. Xu and Dudek proposed the online probabilistic trust
inference model (OPTIMo) [15], which employs Bayesian networks to estimate
human trust based on the autonomous agent’s performance and human behavioral
signals. Soh et al. [16] proposed a Bayesian model that combines Gaussian
processes and recurrent neural networks to predict trust levels for different
tasks. Based on the three properties identified in empirical studies, Guo and
Yang [8] modeled trust as a Beta random variable, parameterized by positive
and negative interaction experiences a person has with a robot. In their
following work, they extended the model to the multi-human-multi-robot case by
introducing trust propagation between agents [18].

For a detailed review of the computational models, see [3].

###  II-B Trust-aware Planning

Endowed with a trust prediction model, a robot can predict how human trust
changes due to moment-to-moment interactions and in turn plan its actions
accordingly. Existing studies in trust-aware decision-making modeled HRI
processes as Markov decision processes (MDPs). Chen et al. [4] proposed the
trust-POMDP to let a robot actively estimate and exploit its human teammate’s
trust. Their human-subject study showed that purely optimizing task
performance in a human-robot team may lead to decreased human trust. Losey and
Sadigh [19] modeled human-robot interaction as a two-player POMDP where the
human does not know the robot’s objective. They proposed 4 ways for the robot
to formulate the human’s perception of the robot’s objective and showed the
robot will be more communicative if it assumes the human trusts the robot and
thus increases the human’s involvement.

###  II-C Reward Shaping in Trust-aware HRI

As pure task-driven rewards lead to low human trust, Chen et al. [4] showed
that pure trust-seeking rewards would correct such issues but lead to
suboptimal task performance. Guo et al. [5] further examined the robot’s
optimal policies under different human behavior assumptions and suggested that
adding a decaying trust-seeking term to the reward function would encourage
the robot to seek high human trust during the initial interaction and maintain
high human trust throughout the entire process. The idea of adding an
additional term to the task reward is called reward shaping, which is a
technique initially developed to accelerate the learning speed of an agent in
reinforcement learning. Ng et al. [20] demonstrated that in an infinite-
horizon MDP, the optimal policy with respect to the shaped reward remains
optimal in the original model, provided that the shaping reward can be
expressed as a particular combination of potential functions. In the following
section, we leverage this potential-based reward-shaping approach to devise a
trust reward that yields minimal performance degradation.

##  III Problem Formulation

In this section, we formulate the sequential HRI problem as a two-player
Markov game. In addition, we present the experience-based trust dynamics
model, which will be later used in the proposed framework.

###  III-A Sequential HRI as a Two-player Game

We consider the scenario where a human $h$ and a robot $r$ work as a team
collaboratively for $N$ rounds, as shown in figure 2. We formulate this
process as a finite-horizon two-player Markov game
$M=\langle\mathcal{S},\mathcal{A}^{h},\mathcal{A}^{r},T,R^{h},R^{r}\rangle$,
where $\mathcal{S}$ is the set of environment states, $\mathcal{A}^{h}$ and
$\mathcal{A}^{r}$ are the sets of actions available to the human and the
robot, $T$ is a Markov kernel from
$\mathcal{S}\times\mathcal{A}^{h}\times\mathcal{A}^{r}$ to $\mathcal{S}$
($\sigma$-algebra omitted) specifying the transitional probability of the
process, and $R^{h}$ and $R^{r}$ are the reward functions of the two agents.
At round $n$, these two agents observe the current state
${s}_{n}\in\mathcal{S}$ of the environment and then take actions. The robot
selects an action $a_{n}^{r}=\pi({s}_{n})\in\mathcal{A}^{r}$ according to its
policy $\pi$. Then, the human observes the robot’s action and takes action
$a_{n}^{h}=f\left(s_{n},a_{n}^{r}\right)\in\mathcal{A}^{h}$ according to his
or her policy $f$. Their actions transition the environment to a new state
${s}_{n+1}$ according to the probability
$T\left(\cdot\middle|{s}_{n},a_{n}^{h},a_{n}^{r}\right)$, and give the human
and the robot rewards
$x_{n}^{h}=R^{h}\left(a_{n}^{h},a_{n}^{r},s_{n},s_{n+1}\right)$ and
$x_{n}^{r}=R^{r}\left(a_{n}^{h},a_{n}^{r},s_{n},s_{n+1}\right)$ respectively.
The game starts from the initial state $s_{1}$ and terminates after $N$ steps
at state $s_{N+1}$, where the latter depends on the policies $\pi$ and $f$.
The goal of the robot is to maximize the expected discounted total payoff

$J_{M}(s_{1})=\mathbb{E}\left[\sum_{n=1}^{N}\gamma^{n-1}x_{n}^{r}\right].$ (1)

The optimal policy, possibly not unique, is a policy $\pi^{*}$ that maximizes
$J_{M}(s_{1})$, i.e., $\pi^{*}=\arg\max_{\pi\in\Pi}J_{M}(s_{1})$, where $\Pi$
is the set of available policies.

To analyze the proposed model, we define the value function. Given any policy
$\pi$ at round $n$, we define the value function
$V_{n}^{\pi}:\mathcal{S}\rightarrow\mathbb{R}$ as the expected discounted
reward by the end of round $N$, i.e.,

$V_{n}^{\pi}(s)=\mathbb{E}_{\pi}\left[\sum_{i=n}^{N}\gamma^{i-n}x_{i}^{r}\right],$
(2)

where $x_{i}$ is the reward received by the robot during round $i$ by
following the policy $\pi$ from state $s$ and thereafter, and the expectation
is over the state-transitions taken upon executing $\pi$.

((a)) Sequential HRI as a two player Markov game. At round $n$, these two
agents observe the current state of the environment $s_{n}\in\mathcal{S}$. The
robot selects an action $a_{n}^{r}=\pi({s}_{n})\in\mathcal{A}^{r}$ according
to its policy $\pi$. Then, the human observes the robot’s action and takes
action $a_{n}^{h}=f\left({s}_{n},a_{n}^{r}\right)\in\mathcal{A}^{h}$ according
to his or her policy $f$. The joint actions then act on the environment, which
transitions to a new state ${s}_{n+1}$ with probability
$T\left(\cdot\middle|{s}_{n},a_{n}^{h},a_{n}^{r}\right)$.

((b)) This information flow chart highlights the roles of the policies $\pi$
and $f$ in the decision process. The robot selects its action $a^{r}_{n}$
based on its policy $\pi$ and the human selects his or her action based on the
behavior model $f$.

Figure 2: A general sequential HRI framework.

Our approach assumes that the human agent plays a supervisory role by
observing the robot’s actions before deciding on their own action. This
mirrors real-world scenarios where humans can intervene in robots’ operations.
The cases where the human and robot agents act simultaneously is a special
setting where the human policy $f$ is independent of the robot’s current
action $a^{r}_{n}$, i.e., $f\left(s_{n},a_{n}^{r}\right)=f^{\prime}(s_{n})$
for some $f^{\prime}:\mathcal{S}\rightarrow\mathcal{A}^{h}$. Furthermore, we
assume the game is Markovian. Non-Markovian games can also be formulated in
this framework through Markovian embedding, i.e., we can include more
variables in the state space to make the transitional probability and the
policies Markovian. Finally, given our assumption that the policy $f$ of the
human is stationary, the game is deemed Markovian from the robot’s
perspective. This allows us to treat the game as an MDP for analysis purposes.

###  III-B Experience-based Trust Dynamics

To model human trust dynamics, we utilize the experience-based trust model
developed in [8]. This model is general enough to be applied to various human-
robot interaction (HRI) settings since it only requires the robot’s
performance as input. In the model, trust $t_{n}$ before the $n$th interaction
is defined as a random variable that follows a Beta distribution
($t_{n}\sim\operatorname{Beta}(\alpha_{n},\beta_{n})$). The two positive shape
parameters, $\alpha_{n}$ and $\beta_{n}$, both greater than or equal to 1,
represent the cumulative positive and negative interaction experience the
human has had with the robot, and they are updated by

$(\alpha_{n+1},\beta_{n+1})=(\alpha_{n}+w^{s}p_{n},\beta_{n}+w^{f}(1-p_{n})),$
(3)

where $w^{s}p_{n}$ and $w^{f}(1-p_{n})$ are the experience gains from the
robot’s success and failure, and parameters $w^{s}$ and $w^{f}$ determine the
unit gains. Here, $p_{n}\in[0,1]$ represents the robot’s performance measure
on the $n$th task.

##  IV Reward Shaping for Trustworthy Policy

Previous research has shown that the robot may exhibit manipulative behavior
to achieve better task rewards at the cost of losing human trust. To prevent
the robot from engaging in this manipulative behavior, we introduce a trust
reward function $R^{t}$ such that, at the end of round $n$, the robot receives
the composite reward $R:=R^{t}+R^{r}$ instead of $R^{r}$. Such reward $R^{t}$
shapes the behavior of the learning agent in a Markov process and thus is
named shaping reward in literature.

By introducing the shaping reward $R^{t}$, the original Markov game
$M=\langle\mathcal{S},\mathcal{A}^{h},\mathcal{A}^{r},T,R^{h},R^{r}\rangle$ is
transformed into a new game
$M^{\prime}=\langle\mathcal{S},\mathcal{A}^{h},\mathcal{A}^{r},T,R^{h},R\rangle$,
where $R$ is the composite reward. Our research question is how to design the
shaping reward $R^{t}$ to promote trust while minimizing the loss of task
performance. In particular, let
$\pi^{\prime*}=\arg\max_{\pi\in\Pi}J_{M^{\prime}}(s_{1})$ be an optimal policy
in $M^{\prime}$. Then the performance loss that the robot suffers amounts to
$V_{1}^{\pi^{*}}(s_{1})-V_{1}^{\pi^{\prime*}}(s_{1})$. Our objective is to
select a proper $R^{t}$ to increase human trust while limiting the performance
loss by some $\epsilon>0$.

###  IV-A Bounding the Performance Loss

Ng et al. showed that, for an infinite-horizon MDP, the optimal policy in
$M^{\prime}$ is also optimal in $M$ if $R^{t}$ is a potential-based shaping
function [20]. A shaping reward $R^{t}$ is said to be potential-based if there
exists some real-valued function $\Phi:\mathcal{S}\rightarrow\mathbb{R}$ such
that for all
$s\in\mathcal{S}-\\{s_{1}\\},a^{r}\in\mathcal{A}^{r},a^{h}\in\mathcal{A}^{h},s^{\prime}\in\mathcal{S}$,

$R^{t}\left(a^{r},a^{h},s,s^{\prime}\right)=\gamma\Phi\left(s^{\prime}\right)-\Phi(s).$

This result indicates which types of shaping rewards can ensure zero
performance loss, i.e.,
$V_{1}^{\pi^{*}}(s_{1})-V_{1}^{\pi^{\prime*}}(s_{1})=0$, for an infinite-
horizon MDP. However, our scenario is different in two aspects: first, we have
a finite horizon; second, we allow some sacrifice in task performance to
improve human trust. Despite these differences, we are inspired by the proof
technique in [20] and propose that a carefully-designed potential-based
shaping reward $R^{t}$ can guarantee a small performance loss in our setting.

######  Theorem 1.

Let
$M=\langle\mathcal{S},\mathcal{A}^{h},\mathcal{A}^{r},T,R^{h},R^{r}\rangle$
and
$M^{\prime}=\langle\mathcal{S},\mathcal{A}^{h},\mathcal{A}^{r},T,R^{h},R\rangle$
be two $N$-horizon Markov games with $R=R^{r}+R^{t}$. Let
$\pi^{*}=\arg\max_{\pi\in\Pi}J_{M}(s_{1})$ and
$\pi^{\prime*}=\arg\max_{\pi\in\Pi}J_{M^{\prime}}(s_{1})$ be two optimal
policies w.r.t. $M$ and $M^{\prime}$ respectively and $\epsilon$ be a positive
number. Suppose there exists some real-valued function
$\Phi:\mathcal{S}\rightarrow\mathbb{R}$ such that

$R^{t}\left(a^{r},a^{h},s,s^{\prime}\right)=\gamma\Phi\left(s^{\prime}\right)-\Phi(s)$
(4)

for all $s,s^{\prime}\in\mathcal{S}$. Then

$V_{1}^{\pi^{*}}(s_{1})-V_{1}^{\pi^{\prime*}}(s_{1})\leqslant\epsilon$ (5)

if

$\mathbb{E}_{\pi^{\prime*}}[\Phi(s_{N+1})]-\mathbb{E}_{\pi^{*}}[\Phi(s_{N+1})]\leqslant\gamma^{-N}\epsilon,$
(6)

where $s_{N+1}$ is the final state of the Markov process when starting from
state $s_{1}$ and following the corresponding policy for $N$ rounds.

###### Proof.

Let

$x_{n}^{t}=R^{t}\left(a_{n}^{h},a_{n}^{r},s_{n},s_{n+1}\right)$ (7)

and define the value function of $M^{\prime}$ similar to Eq. (2) as

$V_{n}^{\prime\pi}(s)=\mathbb{E}_{\pi}\left[\sum_{i=n}^{N}\gamma^{i-n}\left(x_{i}^{r}+x_{i}^{t}\right)\right].$
(8)

By Eq. (2) and Eq. (8), we have

${\displaystyle
V_{1}^{\prime\pi}(s_{1})-V_{1}^{\pi}(s_{1})=\mathbb{E}_{\pi}\left[\sum\limits_{n=1}^{N}\gamma^{n-1}x_{n}^{t}\right].}$
(9)

We can express the second term on the right-hand side by Eqs. (4) and (7) as

$\displaystyle{\displaystyle\mathbb{E}_{\pi}\left[\sum\limits_{n=1}^{N}\gamma^{n-1}x_{n}^{t}\right]}$
(10) $\displaystyle=$
$\displaystyle{\displaystyle\mathbb{E}_{\pi}\left[\sum\limits_{n=1}^{N}\gamma^{n-1}R^{t}\left(a_{n}^{r},a_{n}^{h},s_{n},s_{n+1}\right)\right]}$
$\displaystyle=$
$\displaystyle{\displaystyle\mathbb{E}_{\pi}\left[\sum_{n=1}^{N}\gamma^{n}\Phi(s_{n+1})-\sum_{n=1}^{N}\gamma^{n-1}\Phi(s_{n})\right]}$
$\displaystyle=$
$\displaystyle{\displaystyle\gamma^{N}\mathbb{E}_{\pi}[\Phi(s_{N+1})]-\Phi(s_{1})}.$

By combining Eqs. (9) and (10) we obtain

${\displaystyle
V_{1}^{\prime\pi}(s_{1})-V_{1}^{\pi}(s_{1})=\gamma^{N}\mathbb{E}_{\pi}[\Phi(s_{N+1})]-\Phi(s_{1}).}$
(11)

Now we can bound the performance loss when executing $\pi^{\prime*}$ instead
of $\pi^{*}$ on $M$:

$\displaystyle V_{1}^{\pi^{*}}(s_{1})-V_{1}^{\pi^{\prime*}}(s_{1})$
$\displaystyle=$ $\displaystyle
V_{1}^{\prime\pi^{*}}(s_{1})-V_{1}^{\prime\pi^{\prime*}}(s_{1})$
$\displaystyle+\gamma^{N}(\mathbb{E}_{\pi^{\prime*}}[\Phi(s_{N+1})]-\mathbb{E}_{\pi^{*}}[\Phi(s_{N+1})])$
$\displaystyle\leqslant$
$\displaystyle\gamma^{N}(\mathbb{E}_{\pi^{\prime*}}[\Phi(s_{N+1})]-\mathbb{E}_{\pi^{*}}[\Phi(s_{N+1})])$
$\displaystyle\leqslant$ $\displaystyle\epsilon,$

where the first equality follows from Eq. (11), the first inequality holds
because $\pi^{\prime*}$ is the optimal policy for $M^{\prime}$, and the last
inequality follows from the hypothesis Eq. (6). ∎

Theorem 1 allows us to bound the performance loss by the value of the
potential function on the terminal states of the game. It directly implies the
following corollary:

######  Corollary 1.

Let $M$ and $M^{\prime}$ be two Markov games that satisfy the conditions in
theorem 1, $\pi^{*}$ and $\pi^{\prime*}$ be two optimal policies respectively,
and $R^{t}$ be a potential-based shaping reward. Let
$\mathcal{S}_{s_{1}}^{\pi}(N+1)$ be the set of states that are reachable from
state $s_{1}$ when following policy $\pi$ after $N$ rounds. Then

$V_{1}^{\pi^{*}}(s_{1})-V_{1}^{\pi^{\prime*}}(s_{1})\leqslant\epsilon$ (12)

if

$\max\Phi\left(\mathcal{S}_{s_{1}}^{\pi^{\prime*}}(N+1)\right)-\min\Phi\left(\mathcal{S}_{s_{1}}^{\pi^{*}}(N+1)\right)\leqslant\gamma^{-N}\epsilon.$

###### Proof.

Apparently,

$\displaystyle\mathbb{E}_{\pi^{\prime*}}[\Phi(s_{N+1})]-\mathbb{E}_{\pi^{*}}[\Phi(s_{N+1})]$
$\displaystyle\leqslant$
$\displaystyle\max\Phi\left(\mathcal{S}_{s_{1}}^{\pi^{\prime*}}(N+1)\right)-\min\Phi\left(\mathcal{S}_{s_{1}}^{\pi^{*}}(N+1)\right)$
$\displaystyle\leqslant$ $\displaystyle\gamma^{-N}\epsilon.$

The result follows from theorem 1. ∎

Compared with theorem 1, the condition in corollary 1 is easier to verify.

###  IV-B Trust-seeking Shaping Reward via Experience-based Trust Dynamics

In this section, we apply corollary 1 and the experience-based trust model to
design a shaping reward for encouraging trust-seeking behavior. Recall that,
in the experience-based trust model, human trust is represented as an
experience tuple $(\alpha,\beta)$. Suppose we have a Markov game
$M=\langle\mathcal{S},\mathcal{A}^{h},\mathcal{A}^{r},T,R^{h},R^{r}\rangle$
where the human trust $(\alpha,\beta)$ can be extracted from the state
variable $s$ as $(\alpha,\beta)=g(s)$. We define a potential-based trust
reward as

$\displaystyle R^{t}\left(a^{r},a^{h},s,s^{\prime}\right)$
$\displaystyle=\gamma\phi\left(g\left(s^{\prime}\right)\right)-\phi(g(s))$
(13)
$\displaystyle=\gamma\phi\left(\alpha^{\prime},\beta^{\prime}\right)-\phi(\alpha,\beta),$

with $\phi$ be to determined. Here the actual potential function is
$\Phi:=\phi\circ g$.

Let $R=R^{r}+R^{t}$ be the new reward and
$M^{\prime}=\langle\mathcal{S},\mathcal{A}^{h},\mathcal{A}^{r},T,R^{h},R\rangle$
be the transformed Markov game. We will apply corollary 1 to bound the
performance loss. For any policy $\pi$,

$\Phi\left(\mathcal{S}_{s_{1}}^{\pi}(N+1)\right)=\phi\left(g\left(\mathcal{S}_{s_{1}}^{\pi}(N+1)\right)\right)=\phi\left(\mathcal{E}_{s_{1}}^{\pi}(N+1)\right),$
(14)

where $\mathcal{E}_{s_{1}}^{\pi}(N+1)$ is the collection of reachable trust
state at round $N+1$. Based on Eq. (3), trust $(\alpha_{N+1},\beta_{N+1})$ at
round $N+1$ is

${\textstyle\left(\alpha_{1}+w^{s}\sum_{n=1}^{N}p_{n},\beta_{1}+w^{f}\left(N-\sum_{n=1}^{N}p_{n}\right)\right)},$

where $p_{1},p_{2},\dotsc,p_{N}$ are the robot’s performance values during the
interaction. Since $\sum_{n=1}^{N}p_{n}\in[0,N]$, $(\alpha_{N+1},\beta_{N+1})$
lies on the line

$l_{N+1}=\\{(\alpha_{1}+w^{s}t,\beta_{1}+w^{f}(N-t))\ |\ t\in[0,N]\\},$ (15)

which indicates that

$\mathcal{E}_{s_{1}}^{\pi}(N+1)\subseteq l_{N+1}.$ (16)

Therefore, by Eqs. (14) and (16), it suffices to have

$\max\phi(l_{N+1})-\min\phi(l_{N+1})\leqslant\gamma^{-N}\epsilon$ (17)

for the condition in corollary 1 to hold, where $l_{N+1}$ is given in Eq.
(15).

Eq. (17) constrains the choices of $\phi$ such that the performance loss is
within $\epsilon$. In addition, we should design the function $\phi$ such that
the shaping reward $R^{t}$ optimizes human trust. For example, if we want to
increase human trust, we should reward the robot for if the future state has
higher trust and penalize the robot otherwise. By Eq. (13), one way to achieve
this is to add the following constraints:

$\displaystyle\gamma\phi\left(\alpha^{\prime},\beta^{\prime}\right)-\phi(\alpha,\beta)\geqslant
0,$ $\displaystyle\ \text{if}\
\frac{\alpha^{\prime}}{\alpha^{\prime}+\beta^{\prime}}\geqslant\frac{\alpha}{\alpha+\beta};$
(18)
$\displaystyle\gamma\phi\left(\alpha^{\prime},\beta^{\prime}\right)-\phi(\alpha,\beta)<0,$
$\displaystyle\ \text{otherwise}.$

Here
$\frac{\alpha^{\prime}}{\alpha^{\prime}+\beta^{\prime}}\geqslant\frac{\alpha}{\alpha+\beta}$
indicates trust state $\left(\alpha^{\prime},\beta^{\prime}\right)$ has higher
expected trust compared to $(\alpha,\beta)$, since we assume human trust
follows the Beta distribution $\operatorname{Beta}(\alpha,\beta)$. Another
example is trust calibration. If our goal is to calibrate human trust around a
point $t^{*}$, we can let $R^{t}$ to reward the robot if human trust is moving
towards $t^{*}$ by forcing
$\gamma\phi\left(\alpha^{\prime},\beta^{\prime}\right)-\phi(\alpha,\beta)\geqslant
0$ if
$\left|\frac{\alpha^{\prime}}{\beta^{\prime}+\alpha^{\prime}}-t^{*}\right|\leqslant\left|\frac{\alpha}{\beta+\alpha}-t^{*}\right|$
and
$\gamma\phi\left(\alpha^{\prime},\beta^{\prime}\right)-\phi(\alpha,\beta)<0$
otherwise.

##  V Case Study

To assess our framework’s effectiveness, we simulate a search-and-rescue (SAR)
mission, comparing the robot’s optimal policies, both with and without a
shaping reward.

###  V-A The SAR Mission

The SAR mission was inspired by the work of Wang et al. [21], where a human
and a robot work together to search multiple sites in a town for potential
hazards. At each site, the robot enters first to scan for threats and then
advises the human whether to wear protective gear before entering. However,
wearing the heavy gear is time-consuming, and if there is no threat, it wastes
valuable time. On the other hand, if the human chooses not to wear the
protective gear and there is a threat, they risk injury. The objective is to
complete the mission as quickly as possible while minimizing the human’s
health loss.

We assume that the human-robot team starts to search from site $1$ until site
$N$. At site $n$, a threat indicator $\eta_{n}$ is drawn from a Bernoulli
distribution $\operatorname{Bern}(d_{n})$. There is a threat in site $n$ if
$\eta_{n}=1$ and no threat otherwise. The danger level $d_{n}$ is drawn from
the uniform distribution $\operatorname{U}[0,1]$. The human-robot team does
not know $\eta_{n}$ or $d_{n}$. Instead, prior to the start of the mission,
the team is provided with an estimation $d_{n}^{h}$ of $d_{n}$. Before
entering site $n$, the robot will analyze the site based on its sensory input
and reach a more accurate estimation $d_{n}^{r}$ of $d_{n}$. $d_{k}^{h}$ and
$d_{k}^{r}$ follow Beta distribution
$\operatorname{Beta}(\kappa^{h}d_{n},\kappa^{h}(1-d_{n}))$ and
$\operatorname{Beta}(\kappa^{r}d_{n},\kappa^{r}(1-d_{n}))$ respectively. We
assume that $\kappa^{r}>\kappa^{h}\geqslant 1$ such that the robot has a more
accurate assessment of $d_{n}$ compared with the human. We summarize the
probability model as follows:

Threat indicator $\displaystyle\eta_{n}$
$\displaystyle\overset{\text{i.i.d.}}{\sim}\operatorname{Bern}(d_{n})$ Danger
level $\displaystyle d_{n}$
$\displaystyle\overset{\text{i.i.d.}}{\sim}\operatorname{U}[0,1]$
$\displaystyle\text{Human's estimation of }d_{n}$ $\displaystyle d_{n}^{h}$
$\displaystyle\overset{\text{i.i.d.}}{\sim}\operatorname{Beta}(\kappa^{h}d_{n},\kappa^{h}(1-d_{n}))$
$\displaystyle\text{Robot's estimation of }d_{n}$ $\displaystyle d_{n}^{r}$
$\displaystyle\overset{\text{i.i.d.}}{\sim}\operatorname{Beta}(\kappa^{r}d_{n},\kappa^{r}(1-d_{n}))$

We formulate the SAR mission as a two-player Markov game
$M=\langle\mathcal{S},\mathcal{A}^{h},\mathcal{A}^{r},T,R^{h},R^{r}\rangle$
introduced in section III-A. The state space $\mathcal{S}=[1,\infty)^{2}$
comprises all possible experience pairs $(\alpha,\beta)$ that represent the
human’s trust. The initial state is $s_{1}=(\alpha_{1},\beta_{1})$. The
robot’s action set, denoted as $\mathcal{A}^{r}=\\{0,1\\}$, includes two
options: recommending wearing and recommending not wearing the protective
gear, represented by $a^{r}=1$ and $a^{r}=0$, respectively. Similarly, the
human’s action set, denoted as $\mathcal{A}^{h}=\\{0,1\\}$, includes two
options: wearing or not wearing protective gear, represented by $a^{h}=1$ and
$a^{h}=0$, respectively. We define the robot’s performance at site $n$ as
$p_{n}=\mathbf{1}\left\\{a_{n}^{r}=\eta_{n}\right\\}$, which evaluates to 1 if
the robot’s recommendation agrees with the presence of the threat
($a_{n}^{r}=\eta_{n}$) and 0 otherwise. The state $(\alpha_{n},\beta_{n})$ at
site $n$ transitions to

$(\alpha_{n+1},\beta_{n+1})=(\alpha_{n}+w^{s}p_{n},\beta_{n}+w^{f}(1-p_{n}))$

at site $n+1$ as specified in the experience-based trust model. We assume the
robot has already learned the parameters $w^{s}$ and $w^{f}$ from its previous
interactions with the human. We also assume that the robot and the human share
the same task reward, i.e.,
$R^{h}=R^{r}=-w^{\text{H}}\Delta^{\text{H}}-w^{\text{T}}\Delta^{\text{T}}$,
where $\Delta^{\text{H}}$ and $\Delta^{\text{T}}$ are the time cost and the
health cost for the human-robot team and $w^{\text{H}}$ and $w^{\text{T}}$ are
the corresponding weights. The values of $\Delta^{\text{H}}$ and
$\Delta^{\text{T}}$ are given in table I, and $w^{\text{H}}$ and
$w^{\text{T}}$ are set to 1 and 0.2. The discount factor $\gamma$ is set to
0.9.

We assume that the human follows the reverse-psychology policy as introduced
in [5], where the human will likely comply with the robot’s recommendation
when human trust is high and will do the opposite when trust is low.
Specifically, we have

$\Pr\left(a_{n}^{h}=a_{n}^{r}\right)=\frac{\alpha_{n}}{\alpha_{n}+\beta_{n}}\text{
and }\Pr\left(a_{n}^{h}\neq
a_{n}^{r}\right)=\frac{\beta_{n}}{\alpha_{n}+\beta_{n}},$

where $\frac{\alpha_{n}}{\alpha_{n}+\beta_{n}}$ is the expected human trust
since the human trust follows the beta distribution
$\operatorname{Beta}(\alpha_{n},\beta_{n})$.

Given the model conditions above, we can calculate the probabilities of all
four cases listed in table I, for different actions $a^{r}$. This enables us
to determine the expected immediate reward for each $a^{r}$ and thus apply the
value iteration method to derive an optimal policy for the robot.

Table I: Value Table of $(\Delta^{\text{H}},\Delta^{\text{T}})$ | Protective gear
---|---
Yes ($a^{h}=1$) | No ($a^{h}=0$)
Threat existence | Yes ($\eta=1$) | $(1,300)$ | $(100,50)$
No ($\eta=0$) | $(0,250)$ | $(0,30)$

###  V-B Shaping the Reward

Suppose that we are interested in increasing human trust during the
interaction while limiting the performance loss by $\epsilon$. We apply the
reward-shaping technique developed in section IV to achieve this goal. We
define the shaping reward as
$R^{t}\left(\alpha,\beta,\alpha^{\prime},\beta^{\prime}\right)=\gamma\Phi\left(\alpha^{\prime},\beta^{\prime}\right)-\Phi(\alpha,\beta)$,
where the function $\Phi:\mathcal{S}\rightarrow\mathbb{R}$ is to be
determined. We observe that, from a state $(\alpha,\beta)$, the next state can
either be
$\left(\alpha^{\uparrow},\beta^{\uparrow}\right)=\left(\alpha+w^{s},\beta\right)$,
where the expectation of trust increases, or
$\left(\alpha^{\downarrow},\beta^{\downarrow}\right)=\left(\alpha,\beta+w^{f}\right)$,
where the expectation of trust decreases. Instead of using the generic method
in Eq. (18), we can incentivize the trust-seeking behavior by maximizing the
reward difference between the two future states, i.e., maximizing the value
$R^{t}\left(\alpha,\beta,\alpha^{\uparrow},\beta^{\uparrow}\right)-R^{t}\left(\alpha,\beta,\alpha^{\downarrow},\beta^{\downarrow}\right)$.
Combining with the performance loss constraint in Eq. (17), we obtain the
following optimization problem

$\displaystyle\max$ $\displaystyle
R^{t}\left(\alpha,\beta,\alpha^{\uparrow},\beta^{\uparrow}\right)-R^{t}\left(\alpha,\beta,\alpha^{\downarrow},\beta^{\downarrow}\right)$
(19) s.t.
$\displaystyle\max\Phi(l_{N+1})-\min\Phi(l_{N+1})\leqslant\gamma^{-N}\epsilon,$

where $l_{N+1}$ is defined in Eq. (15). We choose a linear potential function:
$\Phi(\alpha,\beta)=a\alpha+b\beta$ with $a,b\in\mathbb{R}$. With some
algebraic manipulation, (19) becomes a clean linear program:

$\displaystyle\max$ $\displaystyle aw^{s}-bw^{f}$ (20) s.t.
$\displaystyle\frac{-\gamma^{-N}\epsilon}{N}\leqslant
aw^{s}-bw^{f}\leqslant\frac{\gamma^{-N}\epsilon}{N}.$

As the above program is underdetermined, we enforce an extra constraint that
$b=0$ and then obtain the optimal solution
$(a,b)=(\frac{\gamma^{-N}\epsilon}{Nw^{s}},0)$. Therefore, the optimal
potential function is
$\Phi(\alpha,\beta)=\frac{\gamma^{-N}\epsilon}{Nw^{s}}\alpha$, and the shaping
reward is

$R^{t}\left(\alpha,\beta,\alpha^{\prime},\beta^{\prime}\right)=\frac{\gamma^{-N}\epsilon}{Nw^{s}}\left(\gamma\alpha^{\prime}-\alpha\right).$
(21)

###  V-C Simulation Results

Let $M$ be the original Markov game without the shaping reward and
$M^{\prime}$ be the one with shaping reward $R^{t}$. We solve $M^{\prime}$
with 4 different values of $\epsilon$ and, in figure 3, plot the optimal
action $\pi^{\prime*}(\alpha_{1},\beta_{1})$ in $M^{\prime}$, the value
function $V_{1}^{\prime\pi^{\prime*}}(\alpha_{1},\beta_{1})$ of
$\pi^{\prime*}$ in $M^{\prime}$ (defined in Eq. (8)), and the value function
$V_{1}^{\pi^{\prime*}}(\alpha_{1},\beta_{1})$ of $\pi^{\prime*}$ in $M$
(defined in Eq. (2)), for various values of $(\alpha_{1},\beta_{1})$. In the
optimal action plots, the black area corresponds to the states where the
optimal action of the robot is not recommending the human to wear protective
gear, i.e., $a^{r}_{1}=0$, while the white area corresponds to recommending to
wear the gear, i.e., $a^{r}_{1}=1$.

((a)) $\epsilon=0$. No reward shaping.

((b)) $\epsilon=30$.

((c)) $\epsilon=100$.

((d)) $\epsilon=300$.

Figure 3:  Plots of the optimal action $\pi^{\prime*}(\alpha_{1},\beta_{1})$,
value function $V_{1}^{\prime\pi^{\prime*}}(\alpha_{1},\beta_{1})$ in the
Markov game $M^{\prime}$ with shaping reward, and value function
$V_{1}^{\pi^{\prime*}}(\alpha_{1},\beta_{1})$ in the original game $M$, under
different values of $\epsilon$. In the optimal action plots, the black area
corresponds to the states where the optimal action of the robot is not
recommending the human to wear protective gear, i.e., $a^{r}_{1}=0$, and the
white area corresponds to recommending to wear the gear. The robot’s estimated
danger level $d^{r}_{1}$ at the first site is 0.06.

At the first site, the robot’s perceived danger level $d^{r}_{1}$ is 0.06,
which indicates that the “trustworthy” action is to recommend the human not to
wear the gear. In figure 3(a), we set $\epsilon=0$ so there is no shaping
reward. The optimal action is reversed around the 45-degree line, which means
the robot will reverse its action when the human trust crosses a threshold.
This manipulative behavior is consistent with the finding in [5]. In figure
3(b), we set $\epsilon=30$ to allow reward shaping at the cost of a moderate
amount of performance loss. We can observe that the black area grows larger
compared with that of figure 3(a), which implies the trust reward guided the
robot to gain human trust by recommending not to wear the gear. In figure
3(d), we set a large value of 300 for $\epsilon$, and the black area is the
largest among all the settings, indicating that the robot will choose the
“righteous” action in most cases. Moreover, a comparison across
$V_{1}^{\pi^{\prime*}}$ with different $\epsilon$ values shows that the task
reward loss caused by the shaping reward is within $\epsilon$, and initial
states with higher trust have a higher value in $M$ when $\epsilon$ is set to
be higher. This shows the effectiveness of the proposed reward-shaping method.

Finally, we notice an intriguing behavior wherein, as $\epsilon$ increases,
the value function $V_{1}^{\prime\pi^{\prime*}}$ exhibits higher values at
states characterized by lower trust. For instance, given $\epsilon=300$,
states denoted by low $\alpha_{1}$ and high $\beta_{1}$ (low-trust states),
display greater values in contrast to those states with high trust. This
observation may be attributed to two factors. Firstly, with an increment in
$\epsilon$, the shaping reward $R^{t}$ magnitude also escalates, as is
evidenced by Eq. (21). As a result, the pattern of
$V_{1}^{\prime\pi^{\prime*}}$ is largely dictated by the pattern of
$\sum_{n}R^{t}$. Secondly, it appears that the algorithm identifies low-trust
states as having a greater potential for earning $R^{t}$, particularly when
compared to higher-trust states.

In conclusion, the simulation results show that our reward-shaping method
successfully guides the robot to actively gain human trust, overcoming the
manipulative behavior in the pure performance-driven setting.

##  VI Conclusion

In this work, we proposed a framework to balance task reward and human trust.
We formulated the problem as a reward-shaping problem and proposed a novel
framework to solve it. We evaluated the proposed framework in a simulation
scenario where a human-robot team performs a search-and-rescue mission. The
results showed that the proposed framework successfully modifies the robot’s
optimal policy, enabling it to increase human trust with a minimal task
performance cost. However, the work should be viewed in light of the following
limitations. First, we only provide a sufficient condition to guarantee small
performance loss. A necessary condition is in need to complete the theory.
Second, we used linear potential functions for designing the shaping reward.
The effectiveness of other forms of potential functions can be investigated in
future research.

## References

  * [1] J. D. Lee and K. A. See, “Trust in automation: Designing for appropriate reliance,” _Human Factors_ , vol. 46, no. 1, pp. 50–80, 2004.
  * [2] T. B. Sheridan, “Human-Robot Interaction: Status and Challenges,” _Human Factors_ , vol. 58, no. 4, pp. 525–532, 2016.
  * [3] B. C. Kok and H. Soh, “Trust in robots: Challenges and opportunities,” _Current Robotics Reports_ , pp. 1–13, 2020.
  * [4] M. Chen, S. Nikolaidis, H. Soh, D. Hsu, and S. Srinivasa, “Trust-aware decision making for human-robot collaboration: Model learning and planning,” _ACM Transactions on Human-Robot Interaction (THRI)_ , vol. 9, no. 2, pp. 1–23, 2020.
  * [5] Y. Guo, C. Shi, and X. J. Yang, “Reverse psychology in trust-aware human-robot interaction,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 3, pp. 4851–4858, 2021.
  * [6] S. Bhat, J. B. Lyons, C. Shi, and X. J. Yang, “Clustering trust dynamics in a human-robot sequential decision-making task,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 4, pp. 8815–8822, 2022.
  * [7] X. J. Yang, V. V. Unhelkar, K. Li, and J. A. Shah, “Evaluating effects of user experience and system transparency on trust in automation,” in _Proceedings of the 12th ACM/IEEE International Conference on Human-Robot Interaction (HRI ’17)_. ACM, 2017, pp. 408–416.
  * [8] Y. Guo and X. J. Yang, “Modeling and Predicting Trust Dynamics in Human-Robot Teaming: A Bayesian Inference Approach,” _International Journal of Social Robotics_ , vol. 13, pp. 1899–1909, 2021.
  * [9] X. J. Yang, C. Schemanske, and C. Searle, “Toward Quantifying Trust Dynamics: How People Adjust Their Trust After Moment-to-Moment Interaction With Automation,” _Human Factors_ , p. 00187208211034716, 2021.
  * [10] J. D. Lee and N. Moray, “Trust control strategies and allocation of function in human-machine systems,” _Ergonomics_ , vol. 35, no. 10, pp. 1243–1270, 1992.
  * [11] D. Manzey, J. Reichenbach, and L. Onnasch, “Human performance consequences of automated decision aids: The impact of degree of automation and system experience,” _Journal of Cognitive Engineering and Decision Making_ , vol. 6, no. 1, pp. 57–87, 2012.
  * [12] E. J. de Visser, M. M. Peeters, M. F. Jung, S. Kohn, T. H. Shaw, R. Pak, and M. A. Neerincx, “Towards a theory of longitudinal trust calibration in human–robot teams,” _International Journal of Social Robotics_ , vol. 12, no. 2, pp. 459–478, 2020.
  * [13] B. F. Malle and D. Ullman, “A multidimensional conception and measure of human-robot trust,” in _Trust in human-robot interaction_. Elsevier, 2021, pp. 3–25.
  * [14] X. J. Yang, Y. Guo, and C. Schemanske, “From Trust to Trust Dynamics: Combining Empirical and Computational Approaches to Model and Predict Trust Dynamics in Human-Autonomy Interaction,” in _Human-Automation Interaction: Transportation_ , V. G. Duffy, S. J. Landry, J. D. Lee, and N. A. Stanton, Eds., 2023, pp. 253–265.
  * [15] A. Xu and G. Dudek, “Optimo: Online probabilistic trust inference model for asymmetric human-robot collaborations,” in _2015 10th ACM/IEEE International Conference on Human-Robot Interaction (HRI)_. IEEE, 2015, pp. 221–228.
  * [16] H. Soh, Y. Xie, M. Chen, and D. Hsu, “Multi-task trust transfer for human–robot interaction,” _The International Journal of Robotics Research_ , vol. 39, no. 2-3, pp. 233–249, 2020.
  * [17] W.-L. Hu, K. Akash, N. Jain, and T. Reid, “Real-time sensing of trust in human-machine interactions,” _IFAC-PapersOnLine_ , vol. 49, no. 32, pp. 48–53, 2016.
  * [18] Y. Guo, X. J. Yang, and C. Shi, “Enabling Team of Teams: A Trust Inference and Propagation (TIP) Model in Multi-Human Multi-Robot Teams,” in _Proceedings of Robotics: Science and Systems_ , Daegu, Republic of Korea, July 2023.
  * [19] D. P. Losey and D. Sadigh, “Robots that take advantage of human trust,” in _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2019, pp. 7001–7008.
  * [20] A. Y. Ng, D. Harada, and S. Russell, “Policy invariance under reward transformations: Theory and application to reward shaping,” in _ICML_ , vol. 99, 1999, pp. 278–287.
  * [21] N. Wang, D. V. Pynadath, and S. G. Hill, “The impact of pomdp-generated explanations on trust and performance in human-robot teams.” in _AAMAS_ , 2016, pp. 997–1005.