# Robust Rational Polynomial Camera Modelling
for SAR and Pushbroom Imaging

Roland Akiki1,2 Roger Marí 1 Carlo de Franchis1,2
Jean-Michel Morel1 Gabriele Facciolo1

(1 Université Paris-Saclay, CNRS, ENS Paris-Saclay, Centre Borelli, France
2 Kayrros SAS)

###### Abstract

The Rational Polynomial Camera (RPC) model can be used to describe a variety
of image acquisition systems in remote sensing, notably optical and Synthetic
Aperture Radar (SAR) sensors. RPC functions relate 3D to 2D coordinates and
vice versa, regardless of physical sensor specificities, which has made them
an essential tool to harness satellite images in a generic way. This article
describes a terrain-independent algorithm to accurately derive a RPC model
from a set of 3D-2D point correspondences based on a regularized least squares
fit. The performance of the method is assessed by varying the point
correspondences and the size of the area that they cover. We test the
algorithm on SAR and optical data, to derive RPCs from physical sensor models
or from other RPC models after composition with corrective functions.

##  1 Introduction

Developing a remote sensing application requires a set of tools, one of which
is geolocation. Geolocation relates the 3D world coordinates to the 2D image.
This is represented by means of a projection function
$\mathcal{P}:\mathbf{R}^{3}\rightarrow\mathbf{R}^{2}$, that maps 3D points to
the image plane, and its inverse, the localization function
$\mathcal{L}:\mathbf{R}^{2}\times\mathbf{R}\rightarrow\mathbf{R}^{3}$. When
all the physical phenomena and components involved in the acquisition process
are known, the geolocation functions can be defined by a chain of operations
that model such factors, in what is known as a physical or rigorous sensor
model.

Pushbroom scanners are the most common optical satellite image acquisition
system, typically consisting of a single line of pixel sensors mounted on a
platform that captures each line of the image at a different moment in time.
As a result, the exterior orientation parameters, i.e. the perspective center
and the attitude angles, change from line to line. The intrinsic parameters
(e.g. pixel size, focal length, lens distortion), related to the physical
design of the sensor, are constant across the image [1]. A detailed
description of a simplified physical sensor model for pushbroom scanners can
be found in [2].

Figure 1: The RPC model is derived using a grid of 3D CNPs (Control Points)
and its projection onto the satellite image.

In the case of Synthetic Aperture Radar (SAR) images, the most used physical
sensor model is the Range-Doppler model detailed in [3]. SAR satellites send
an electromagnetic wave that is reflected on the ground. The image is acquired
line by line (similar to a pushbroom system), and the position of a ground
patch in the image is related to its distance to the sensor, known as the
range. The Range-Doppler model is constructed based on ephemeris data (time,
position and velocity samples along the orbit) and the acquisition timing
information. The ephemeris data needs to be interpolated to obtain continuous
geolocation functions along the orbit.

Image vendors have adopted the generic Rational Polynomial Camera111Also
referred to as Rational Polynomial Coefficients camera model or Rational
Function Model in the literature. (RPC) model to save customers from having to
deal with the complex specificities of rigorous sensor models. The RPC model
is independent of physical properties and offers flexibility to work with
different coordinate systems. RPCs have become essential metadata to process
satellite images in a generic way, from different sources and for multiple
tasks, e.g. photogrammetry and radargrammetry based 3D reconstruction or image
ortho-rectification and coregistration.

In this article, we describe a terrain-independent algorithm to fit a RPC
model from a physical sensor model or any other geolocation model. Our
contributions are:

  * -

An open-source implementation of the method as an easy-to-use Python package,
which is available at https://github.com/cmla/rpcfit.

  * -

An evaluation of the algorithm’s precision and robustness based on real
scenarios. We test our method using Sentinel-1 and WorldView-3 images, to fit
a SAR physical sensor model or correct an existing RPC model by composing it
with a complementary transformation.

###  1.1 Mathematical formulation of the RPC model

The RPC model defines the projection function $\mathcal{P}$ as

$r_{n}=\frac{a\left(X_{n},Y_{n},Z_{n}\right)}{b\left(X_{n},Y_{n},Z_{n}\right)}\quad
c_{n}=\frac{e\left(X_{n},Y_{n},Z_{n}\right)}{f\left(X_{n},Y_{n},Z_{n}\right)},$
(1)

where $a,b,e,f$ are cubic polynomials. $X,Y,Z$ represent the longitude,
latitude and height of a 3D point; and $r,c$ are the row and column of its
projection on the image plane.

Equation 1 uses normalized coordinates for better numerical stability, hence
the subscript $n$. Normalized values are in the range [-1, 1] and they are
obtained using two scalars, an offset and a scale factor:
$X_{n}=(X-X_{\text{offset}})/X_{\text{scale}}$, where $X$ could be $r,c,X,Y$
or $Z$ from Equation 1.

Each RPC polynomial $p$ is defined by 20 coefficients as

$\displaystyle\textstyle\begin{aligned}
p\left(X,Y,Z\right)&=\,\,p_{0}+p_{1}Z+p_{2}Y+p_{3}X+p_{4}ZY+p_{5}ZX\\\
+&p_{6}YX+p_{7}X^{2}+p_{8}Y^{2}+p_{9}Z^{2}+p_{10}ZYX\\\
+&p_{11}Z^{2}Y+p_{12}Z^{2}X+p_{13}Y^{2}Z+p_{14}Y^{2}X\\\
+&p_{15}ZX^{2}+p_{16}YX^{2}+p_{17}Z^{3}+p_{18}Y^{3}+p_{19}X^{3},\end{aligned}$

(2)

where $p_{i}$ is the $i$-th coefficient of $p$. Since we set $p_{0}=1$ for the
RPC denominator polynomials, a total of 78 coefficients need to be determined
to define $a$, $b$, $e$ and $f$ in Equation 1.

##  2 Related work

RPCs have been used for high-resolution optical satellite imaging since the
launch of Ikonos in 1999 [1, 4, 5, 6]. In the last decade, they have been
proven to be extremely accurate for SAR acquisition systems as well [7, 8].

The RPC model of a satellite image can be constructed using a set of
correspondences between image and object space coordinates. Depending on the
nature of these correspondences, the literature can be classified into
terrain-dependent or independent methods (or a combination of both). Terrain-
dependent strategies use Ground Control Points (GCPs), whose object and image
coordinates are known in advance relying on manual labeling or on-site
measurements. Oppositely, terrain-independent methods derive virtual sets of
2D-3D point correspondences from other geolocation functions, usually a
physical sensor model. Once the point correspondences are available, least
squares algorithms are typically used to estimate the RPC coefficients that
minimize the error between the projected 3D points and their image locations.

Several works have underlined the importance of using uniformly distributed
points in sufficient amount, covering the different parts of the image and the
whole altitude range of the scene [4, 6]. As a result, regularized least
squares methods have become widely used to gain robustness to different
configurations and enforce well-conditioned normal equations [4, 6, 8, 9].
Additionally, terrain-dependent strategies have explored the selection of
optimal and balanced subsets of GCPs, e.g. [4] propose a bucketing strategy or
[9] study the benefits of encouraging correspondences located at building
edges in urban scenarios. In contrast, terrain-independent methods can
arbitrarily generate regular sets of points, but require special care to the
boundaries and density of the structure, e.g. [8] investigate the impact of
different number of elevation layers for flat and mountainous areas.

##  3 Method

We follow a terrain-independent approach similar to [4] to fit an RPC model to
another input geolocation model. The data used to fit the model consists of a
3D grid of uniformly distributed Control Points (CNPs) within some longitude,
latitude and altitude boundaries (Fig. 1). The 2D image point of each CNP can
be obtained by projecting it with the input geolocation model, so that each
sample results in 5 normalized values, i.e. $(X_{i},Y_{i},Z_{i},r_{i},c_{i})$.
For simplicity, we drop the subscript $n$ of normalized coordinates and
replace it by $i$ to refer to the sample index.

Using $N$ CNPs, Equation 1 can be rewritten as a system of equations, in
matrix form, following the derivation of [4]:

$WTI-WG=0,$ (3)

where

$\displaystyle W=$ diag

$\displaystyle\left[\frac{1}{b(\bm{X}_{1})},...,\frac{1}{b(\bm{X}_{N})},\frac{1}{f(\bm{X}_{1})},...,\frac{1}{f(\bm{X}_{N})}\right]$

$\displaystyle b(\bm{X}_{i})=b(X_{i},Y_{i},Z_{i}),$ $\displaystyle T=$
$\displaystyle\,\text{block diag}\left[M_{r},M_{c}\right]\quad
M_{r},M_{c}\in\mathbf{R}^{N\times 39}$

$\displaystyle M_{r_{i}}=$
$\displaystyle\left[1,Z_{i},Y_{i},...,X^{3}_{i},-r_{i}Z_{i},-r_{i}Y_{i},...,-r_{i}X^{3}_{i}\right]$
$\displaystyle M_{c_{i}}=$
$\displaystyle\left[1,Z_{i},Y_{i},...,X^{3}_{i},-c_{i}Z_{i},-c_{i}Y_{i},...,-c_{i}X^{3}_{i}\right],$
$\displaystyle I=$
$\displaystyle\left[a_{0},...,a_{19},b_{1},...,b_{19},e_{0},...,e_{19},f_{1},...,f_{19}\right]^{T},$
$\displaystyle G=$
$\displaystyle\left[r_{0},...,r_{N},c_{0},...,c_{N}\right]^{T}.$

In Equation 3, $W$ is a weight matrix with shape $2N\times 2N$,
$b(\bm{X}_{i})$ denotes the RPC polynomial $b$ evaluated with the 3D
coordinates of the $i$-th CNP; $T$ is the design matrix with shape $2N\times
78$; $I$ is the solution vector with the 78 RPC coefficients necessary to
determine $a$, $b$, $e$ and $f$ in Equation 1; and $G$ is a vector of length
$2N$ containing the CNPs image coordinates.

Equation 3 can be solved by least squares minimization to estimate $I$, using
the normal equation

$T^{T}W^{2}TI-T^{T}W^{2}G=0.$ (4)

To increase numerical stability, ridge estimation regularization [9, 8] is
often added so that the normal equation becomes

$(T^{T}W^{2}T+h^{2}E)I-T^{T}W^{2}G=0,$ (5)

where $E$ is the identity matrix and $h$ is a scalar controlling the
regularization that is applied. To choose the best regularization factor $h$,
the L-curve criterion was introduced in [8]. This heuristic computes the log
norm of the solution $\left(\operatorname{log}\left\|I\right\|_{h}\right)$
versus the log norm of the residual $\left(\operatorname{log}\left\|WTI-
WG\right\|_{h}\right)$ across different values of $h$ that extend from the
minimal to the maximal singular value of T. This curve usually has a L-shape,
in which the optimum corresponds to the maximum regularization parameter that
achieves a small residual. The value corresponding to the corner of the curve,
at the position of maximal curvature, is taken to set $h$ automatically. For
non-weighted regularized least squares (i.e. weights are set to the identity,
$W=E$), the L-curve criterion is fast since the curvature can be computed with
closed form expressions [10].

Therefore, we first set $W^{(0)}=E$ (where the superscript denotes the
iteration number) and use the L-curve criterion to determine the optimal $h$
and an initial solution $I^{(0)}$. Then, for $i\geq 1$ , $W^{(i)}$ is
determined from $I^{(i-1)}$ and is plugged in Equation 5 to solve for
$I^{(i)}$ iteratively (the SVD least squares solver is used for stability).
The iterations stop when the change in terms of RMSE between the RPC projected
CNPs and their image coordinates becomes lower than a tolerance value. After
convergence some final ICCV (Iteration by Correcting Characteristic Value [8])
iterations are computed to remove possible biases introduced by the
regularization. The same stopping criterion based on the RMSE improvement is
used. Each ICCV iteration $k$ can be expressed as

$(T^{T}(W^{(k)})^{2}T+E)I^{(k)}=T^{T}(W^{(k)})^{2}G+I^{(k-1)}.$ (6)

##  4 Experiments

###  4.1 Data and use cases description

  * •

SAR. 32 Sentinel-1 SAR images (Table 1). The data is in Interferometric Wide
Swath (IW) mode, each product contains 3 subswaths, and each subswath contains
multiple bursts that need to be stitched together to get a continuous image.
We construct the Range-Doppler physical sensor model and use the method from
Section 3 to fit an equivalent RPC model for each image of the dataset.

  * •

Optical. 47 WorldView-3 panchromatic images (Table 1), from the 2016 IARPA
Multi-View Stereo 3D Mapping Challenge [11]. The original RPCs of the images
exhibit small inaccuracies, mainly due to inexact knowledge of the satellite
attitude angles, which cause a 3D point to be projected to non corresponding
pixels across different images. Bundle adjustment (BA) algorithms are a well-
known approach to correct RPC errors [5, 12, 13]. We apply a BA similar to
[12] to correc the projection function $\mathcal{P}$ of each RPC into a new
$\mathcal{P}_{\text{BA}}$, expressed as

$\mathcal{P}_{\text{BA}}(X)=\mathcal{P}(R(X-T-C)+C),$ (7)

where $X$ is a 3D point. That is, each RPC is corrected by applying a
translation $T$ followed by a rotation $R$ around an approximate camera center
$C$, before applying the original projection $\mathcal{P}$. $C$ is derived by
regressing a projective model from each RPC model. Our method from Section 3
is used to fit $\mathcal{P}_{\text{BA}}$ from the composition of $\mathcal{P}$
with $T,R$, $C$.

| SAR dataset | Optical dataset
---|---|---
platform | S1 A/B | WorldView-3
number of images | 32 | 47
geographic area | Albania | Argentina
(lon, lat) center | (18.82,41.02) | (-58.61,-34.47)
first acquisition date | 2019-08-03 | 2014-11-15
last acquisition date | 2020-02-05 | 2016-01-13
altitude range (m) | $\left[-533,2969\right]$ | $\left[-513,548\right]$
Table 1: SAR and optical data used in the experiments. The altitude range of the area covered by each collection of images is defined using the [min, max] values from the corresponding SRTM digital elevation model $\pm 500$ m to consider tall buildings or fine irregularities beyond bare ground level. (a) SAR |  (b)
---|---
(c) Optical |  (d)
Figure 2: RPC fitting error varying grid length and surface area on the SAR
(a-b) and optical data (c-d) described in Section 4.1. Each vertical bar
corresponds to the $[-\sigma/2,\,\mu,\,\sigma/2]$ values of the root-mean-
square error (RMSE) evaluated across the different images of the datasets, in
each dimension of the image plane, where $\mu$ corresponds to the mean and
$\sigma$ to the standard deviation.

###  4.2 Performance assessment

To apply the method described in Section 3, the altitude limits of the grid of
CNPs are set using the altitude ranges in Table 1. To assess the RPC fitting,
we use a grid of Check Points (CKPs), which are located in the middle of each
pair of consecutive CNPs. The RMSE between the image coordinates obtained by
projecting the CKPs with the output RPC and the image coordinates obtained
using the input geolocation model, measured in pixels, is used as evaluation
metric. We set a tolerance of $10^{-10}$ for the stopping criterion based on
the RMSE improvement, and a maximum of $20$ iterations for the weighted least
squares and the ICCV iterations.

Two type of experiments were conducted to assess the performance and
robustness of the method:

  * -

Varying surface area. For each image, we fit different RPCs by gradually
increasing the longitude and latitude limits of the grid of CNPs from a small
square centered at the image center to a larger square including the entire
image. The number of CNPs is fixed, with $50\times 50$ samples in the
longitude and latitude dimensions and 10 elevation layers (25000 points in
total).

  * -

Varying grid length. For each image, we fit different RPCs, by increasing the
number of CNPs in each elevation layer, i.e. $n\times n$ where $n$ is the grid
length. The longitude and latitude boundaries of the grid are fixed using the
equivalent limits of the image plane.

Overall, the results presented in Fig. 2 show that the RPCs constructed with
our method approximate the geolocation models with very high accuracy for the
two datasets and the two scenarios outlined in Section 4.1. The RMSE is in the
order of $10^{-4}$ pixels or less in both dimensions of the image plane for
the majority of configurations that were tested, which emphasizes the
robustness of the method.

The experiments with different grid lengths (Fig. 2b and 2d) show that 10
samples in the longitude and latitude directions is already a good choice and
increasing this number beyond 20 does not result in significant improvements.
The experiments with varying surface area (Fig. 2a and 2c) show that both the
overall RMSE values and its variation across the different images increase
with the size of the area being fitted. This is probably due to the fact that
the SAR physical sensor model is less smooth for large neighborhoods. A
similar behavior is obtained with the optical dataset, where the original RPCs
are known to behave locally as an affine camera [5, 13].

##  5 Conclusion

This article described an automatic algorithm to fit the RPC model of a
satellite image in a terrain-independent manner. The inputs of the method are
a regular grid of 3D points (CNPs), with multiple elevation layers, and the 2D
locations of the points on the image plane. We evaluated the method on real
scenarios using collections of SAR and optical satellite images, and assessed
its performance by varying the CNPs configuration. Finally, we release an
open-source implementation of the algorithm as an easy-to-use Python package.

## Acknowledgements

Work partly financed by IDEX Paris-Saclay IDI 2016, ANR-11-IDEX-0003-02,
Office of Naval research grant N00014-17-1-2552 and N00014-20-S-B001, DGA
Astrid project « filmer la Terre » no ANR-17-ASTR-0013-01, MENRT, and by a
grant from Région Île-de-France.

## References

  * [1] J. Grodecki,  “IKONOS stereo feature extraction - RPC approach,”  in ASPRS Annual Conference, 2001.
  * [2] C. de Franchis, E. Meinhardt-Llopis, D. Greslou, and G. Facciolo,  “Attitude refinement for orbiting pushbroom cameras: a simple polynomial fitting method,”  Image Processing On Line, vol. 2015, pp. 328–361, 2015.
  * [3] J.C. Curlander,  “Location of spaceborne SAR imagery,”  TGRS, vol. GE-20, no. 3, pp. 359–364, 1982.
  * [4] C.V. Tao and Y. Hu,  “A comprehensive study of the rational function model for photogrammetric processing,”  Photogrammetric Engineering & Remote Sensing, vol. 67-12, 2001\.
  * [5] C.S. Fraser, G. Dial, and J. Grodecki,  “Sensor orientation via RPCs,”  ISPRS, vol. 60, no. 3, pp. 182–194, 2006.
  * [6] T. Long, W. Jiao, and G. He,  “RPC estimation via $\ell_{1}$-norm regularized least squares (L1LS),”  TGRS, vol. 53, no. 8, 2015.
  * [7] G. Zhang, W. Fei, Z. Li, X. Zhu, and D. Li,  “Evaluation of the RPC model for spaceborne SAR imagery,”  Photogrammetric Engineering & Remote Sensing, vol. 76-6, pp. 727–733, 2010.
  * [8] L. Zhang, X. He, T. Balz, X. Wei, and M. Liao,  “Rational function modeling for spaceborne SAR datasets,”  ISPRS, vol. 66, no. 1, pp. 133–145, 2011.
  * [9] Y. Wang, Y. Zhang, and N. Su,  “RPC estimation via feature points for urban areas,”  in IGARSS, 2016, pp. 6684–6687.
  * [10] P.C. Hansen and D.P. O’Leary,  “The use of the L-curve in the regularization of discrete ill-posed problems,”  SIAM Journal on Scientific Computing, vol. 14, no. 6, pp. 1487–1503, 1993.
  * [11] M. Bosch, Z. Kurtz, S. Hagstrom, and M. Brown,  “A multiple view stereo benchmark for satellite imagery,”  in AIPR Workshop, 2016, pp. 1–9.
  * [12] R.A. Beyer, O. Alexandrov, and S. McMichael,  “The Ames Stereo Pipeline: NASA’s open source software for deriving and processing terrain data,”  Earth and Space Science, vol. 5, no. 9, pp. 537–548, 2018.
  * [13] R. Marí, C. de Franchis, E. Meinhardt-Llopis, and G. Facciolo,  “To bundle adjust or not: A comparison of relative geolocation correction strategies for satellite multi-view stereo,”  in ICCV Workshops, 2019.


# Self-Supervised Time-to-Event Modeling
with Structured Medical Records

Ethan Steinberg, Yizhe Xu, Jason Fries, Nigam Shah

(January 8, 2023)

###### Abstract

Time-to-event models (also known as survival models) are used in medicine and
other fields for estimating the probability distribution of the time until a
particular event occurs. While providing many advantages over traditional
classification models, such as naturally handling censoring, time-to-event
models require more parameters and are challenging to learn in settings with
limited labeled training data. High censoring rates, common in events with
long time horizons, further limit available training data and exacerbate the
risk of overfitting. Existing methods, such as proportional hazard or
accelerated failure time-based approaches, employ distributional assumptions
to reduce parameter size, but they are vulnerable to model misspecification.
In this work, we address these challenges with MOTOR a self-supervised model
that leverages temporal structure found in large-scale collections of
timestamped, but largely unlabeled events, typical of electronic health record
data. MOTOR defines a time-to-event pretraining task that naturally captures
the probability distribution of event times, making it well-suited to
applications in medicine. After pretraining on 8,192 tasks auto-generated from
2.7M patients (2.4B clinical events), we evaluate the performance of our
pretrained model after fine-tuning to unseen time-to-event tasks. MOTOR-
derived models improve upon current state-of-the-art C statistic performance
by 6.6% and decrease training time (in wall time) by up to 8.2 times. We
further improve sample efficiency, with adapted models matching current state-
of-the-art performance using 95% less training data.

Keywords— time-to-event, survival, deep learning, self-supervised learning,
transfer learning, foundation models

##  1 Introduction

Time-to-event models (which are also called survival models ) are used for
estimating the probability distribution of the time elapsed until a particular
event occurs given a set of features. Time-to-event models have two primary
advantages over traditional classification approaches. First, they can
estimate quantities of interest, such as the hazard rate at a particular time
or the median time until an event, that are difficult or impossible to compute
in other approaches. Second, time-to-event models account for various types of
censoring induced by different data collection mechanisms. Adjusting for
censoring is critical for obtaining unbiased estimates, which in turn is
important for creating fair models (Chang et al., 2022). As a healthcare
example, the 2013 ACC/AHA pooled cohort equations are used to identify
patients at high risk of atherosclerotic cardiovascular disease in 10 years
and determine the appropriate statin prescription (Goff et al., 2014). In non-
clinical settings, time-to-event models are used for churn prediction
(Perianez et al., 2016), real time advertisement bidding (Wu et al., 2015),
recommendation systems (Jing and Smola, 2017), and other business analytics.

Despite the benefits of time-to-event models, they come with a significant
downside of requiring large training sets to estimate the entire time-to-event
probability distribution, i.e., the probability of the event happening at
every time point rather than at a single time point. A variety of methods have
been proposed to mitigate this issue, generally by adding assumptions about
the time-to-event distribution. These assumptions, such as proportional
hazards (Cox, 1972) or parametric distributions (Martinsson, 2016; Avati et
al., 2020), significantly reduce the number of parameters necessary to learn a
time-to-event model, but may also degrade predictive performance when the
assumptions do not hold. Regardless of these techniques, learning higher-order
interactions between features in high-dimensional datasets is challenging with
limited uncensored labels.

Recent work in deep learning has approached the challenge of limited training
data by relying on transfer learning via increasingly powerful foundation
models (Bommasani et al., 2021). Foundation models leverage advances in self-
supervised learning to train models using large collections of unlabeled data.
The resulting pretrained models learn broadly useful features that can then be
rapidly adapted to new tasks via fine-tuning, substantially lowering the
amount of labeled data and time required to train new models. Self-supervised
learning requires designing a pretraining task that captures generally useful
patterns in unlabeled data, typically by predicting intrinsic properties of
the data itself. For example, autoregressive pretraining in NLP predicts the
next word in a sequence conditioned on prior words (Radford et al., 2019). In
structured electronic health records (EHRs), which consist of text-like
sequences of timestamped clinical events, autoregressive pretraining has also
resulted in strong performance and improved sample efficiency in fine-tuned
models (Steinberg et al., 2021). However large-scale pretraining is
underexplored for time-to-event modeling, with its benefits remaining unknown.
Moreover, standard autoregressive pretraining disregards much of the temporal
structure present in EHRs and it is unclear how these autoregressive
objectives perform when used for more complex time-to-event modeling.

In this work, we introduce MOTOR (Many Outcome Time Oriented Representations),
an approach for learning highly effective time-to-event models using transfer
learning. We build on existing transfer learning techniques by developing a
new self-supervised pretraining task designed specifically for time-to-event
modeling. We then train a 100M parameter Transformer (Vaswani et al., 2017)
model using a flexible piecewise exponential time-to-event modeling with our
specialized pretraining task. The resulting Transformer model is then fine-
tuned on arbitrary time-to-event target tasks of interest. Our main
contributions are the following:

  1. 1.

We introduce a self-supervised pretraining task for learning time-to-event
models and demonstrate that the resulting model helps us train significantly
better time-to-event models than state-of-the-art approaches with an average
improvement in C statistic of 6.6% across several tasks.

  2. 2.

We explain how to implement models that learn our time-to-event pretraining
task through large-scale piecewise exponential models in a parameter efficient
manner, by recognizing that these models can share time-dependent risk factors
between tasks.

  3. 3.

We evaluate several properties of MOTOR, including the choice of pretraining
task, sample efficiency, and compute requirements. We find that time-to-event
pretraining outperforms a more commonly used autoregressive pretraining
objective in EHR data. We also find that fine-tuned MOTOR matches performance
of state-of-the-art methods while using 95% as less training data. Finally,
fine-tuning reduces training time with large datasets by up to 8.2 times (in
wall time) over Random Survival Forest.

##  2 Related Work

###  2.1 Time-to-event Models

Accelerated failure time (AFT) models are a fully parametric approach for
modeling time-to-event outcomes where the time-to-event distribution is
parameterized with a specified probability distribution. AFT models assume
that the joint effect of features is to accelerate or decelerate the time to
an event by some constant, in which the effect of features is learned by
optimizing a loss function, e.g., negative log-likelihood function. Common
probability distributions include exponential, gamma, log-normal, log-
logistic, and Weibull. To allow for additional flexibility, piecewise
distributions that enable the use of different distribution parameters across
different regions of time can be used, with the piecewise exponential
distribution (Friedman, 1982) being a distinctly popular choice. More
recently, AFT models have been combined with deep learning by using a neural
network to map between input features and the parameters of specified
distributions (Fornili et al., 2014; Martinsson, 2016; Ren et al., 2018;
Kopper et al., 2020; Avati et al., 2020).

Cox PH models (Cox, 1972) are another common approach for time-to-event
modeling that reduce the complexity of the model space by assuming a constant
hazard ratio between instances over time. In situations where the proportional
hazards (PH) assumption does not hold, Cox PH models may perform poorly.
Similarly to AFT models, Cox models can also be combined with deep learning by
using a neural network to map between the input features and the learned
constant hazard ratio (Katzman et al., 2018).

Random Survival Forests (RSFs) (Ishwaran et al., 2008a) are a nonparametric
approach for time-to-event modeling that works by extending random forests
(Breiman, 2001) to time-to-event tasks. RSFs estimate cumulative hazard
functions using the Nelson-Aalen estimator. The best splitting variable at
each node is selected to maximize the survival differences between the two
daughter nodes. Similar to random forests, RSFs apply bagging, random
subsampling, and random variable selection techniques to improve model
performance. Moreover, RSFs automatically identify nonlinear relationships
among variables via the recursively partitioning scheme. Athey et al. (2019)
further introduces “honest” trees for building RSFs to mitigate overfitting in
which separate data sets are used for tree construction and prediction. One
major limitation of RSFs is that they are computationally inefficient when
applied to large-scale data (Rasmy et al., 2021a).

There has also been some work exploring how to use multi-task learning to
better fit time-to-event models in cases with limited data. Li et al. (2016)
proposed treating the time bins in piecewise time-to-event models as instances
of related tasks and explicitly regularizing models to share weights across
time bins. Gao and Cui (2022) explored using explicitly chosen “compatible”
clinical tasks (i.e. clinical tasks with similar risk structure) to better
learn neural network featurized Cox models. Wang and Sun (2021) evaluated
learning mortality and length of stay time-to-event models simultaneously
while training cancer time-to-event models. MOTOR differs from these prior
approaches by focusing on scaling up the number of auxiliary tasks
dramatically (up to 8,192 tasks) as part of its pretraining process.

###  2.2 Deep Learning For Structured EHR Data

There is a wide variety of work applying deep learning to structured
Electronic Health Records (EHR) data to train models for answering complex
medical questions. Early approaches, best exemplified by DeepPatient (Miotto
et al., 2016) and the feedforward model in Rajkomar et al. (2018), were
analogous to bag-of-words featurization from NLP and required transforming
patient EHR data into sparse, count-based feature matrices where each column
represents the number of times a particular medical concept occurred within a
patient’s record. While this approach is parameter efficient, it requires
making additional, often non-obvious, feature engineering choices to capture
temporal information about events.

Recent deep learning architectures have focused on directly incorporating the
temporal structure of EHR data by using either recurrent neural networks
(RNNs) or attention based models such as Transformers. These time-oriented
neural network models take sequential EHR data as input and model interactions
between events over time. Both standard neural network architectures and EHR
specific neural network architectures have been proposed for the EHR setting
such as GRU-D (Che et al., 2016), RETAIN (Choi et al., 2016a), and the
Transformer-based BEHRT (Li et al., 2019).

While older deep learning methods for EHR data focused on training models end-
to-end using small patient cohorts, recent methods increasingly use self-
supervised learning to pretrain models using millions of patient records. The
resulting foundation models (Bommasani et al., 2021) demonstrate strong
transfer learning properties and better performance in limited data settings
via fine-tuning and other adaptation techniques. Various self-supervised
objectives have been proposed for learning meaningful patterns in EHR data.
DeepPatient (Miotto et al., 2016) used a denoising autoencoder objective to
learn effective feedforward models. Recent approaches use objectives informed
by language modeling in NLP. CLMBR (Steinberg et al., 2021) proposed an
autoregressive, “next day” code prediction task to train RNN models. Many
authors have presented work using masked language modeling objectives to
predict masked or “corrupted” tokens in an input stream (Li et al., 2019;
Rasmy et al., 2021b; Pang et al., 2021; Zeng et al., 2022).

Prior EHR deep learning models have largely focused on traditional
classification tasks, however, there have been several efforts in combining
deep learning with time-to-event modeling on EHR data. Avati et al. (2020)
uses both RNNs and feedforward neural networks to power AFT models for
estimating each patient’s time until death. Rasmy et al. (2021a) uses RNNs as
a backbone for training Cox models that predict hospital mortality and
ventilation needs. Wen et al. (2022) uses feedforward networks to parameterize
Cox and DeepHIT (Ryu et al., 2020) models for predicting time until discharge.
MOTOR differs from this prior work in two respects. First, we generally use
larger and more flexible Transformer based piecewise exponential models that
can handle more complex interactions. Second, in order to handle the
corresponding increase in parameters, we focus on the benefits of pretraining
at a massive scale.

##  3 Methods

###  3.1 Problem Setup

Figure 1: An example electronic health record (EHR) structured using the
Observational Medical Outcomes Partnership (OMOP) Common Data Model. Each
circle represents a day containing 1 or more clinical events (black) or no
events (grey). Events correspond to different medical codes, color-coded by
OMOP category, with examples labeled in the bottom boxes.

We consider a set of $N$ patients and their EHRs $X$, where $X_{i}$ denotes
all information for patient $i$. Patient $X_{i}$ consists of a sequence of
timestamped clinical events, broadly defined as any action recorded in an EHR
using a medical code. Codes are symbols drawn from a finite vocabulary defined
by 26 different medical ontologies specified by our EHR data provider. For
example, many EHRs use ICD10 codes (Organization, 2004) to indicate diseases,
CPT codes (Dotson, 2013) to indicate procedures, RxNorm codes (Nelson et al.,
2011) to indicate drugs and so forth. Let $X_{ij}$ denote the $j$-th event in
the sequence, where each event consists of a code for that event and the
corresponding timestamp for when that event occurs. Timestamps may occur at
minute or day-level resolution and can co-occur with other events. In this
work, we only consider coded data and discard clinical notes and other
unstructured data.

Figure 1 gives an example EHR record to demonstrate the general temporal
structure of patient data. The x-axis shows the patient timeline that is
offset by date of birth. On day 0, we record the patient characteristics
(e.g., sex, race, and ethnicity). As time progresses, on day 29834, we observe
lab measurements (e.g., pulse rate and blood pressure), medical procedures
(e.g., CPT4/99214 Office/Other Outpatient Visit), and diagnosis (e.g.,
ICD10/I25.10 Coronary Arteriosclerosis); on day 31626, we obtain information
on clinical drugs (e.g., RxNorm /197361 Amlodipine 5 MG).

We are interested in learning models to predict time-to-event tasks for
patients in our dataset. For each patient $i$, let $T_{i}$ and $C_{i}$ denote
the time to an event of interest and censoring time, respectively.
$U_{i}=\mathrm{min}(T_{i},C_{i})$ is the observed follow-up time, and
$\delta_{i}=\mathbb{1}\\{T_{i}\leq C_{i}\\}$ is the event indicator, which is
1 if the patient experiences the event at $U_{i}$ or 0 if the patient is
censored at $U_{i}$. Our estimand of interest is

$\mathbb{P}(t|X_{i})=\mathbb{P}(T_{i}=t),$ (1)

where $t$ is an arbitrary prediction time.

###  3.2 MOTOR

We introduce MOTOR, a Transformer-based neural network model that is
pretrained using a self-supervised, time-to-event objective defined using
structured EHR data. Figure 2 provides an overview of the overall approach,
with patient records being fed into a Transformer backbone to convert them
into latent representations that are then used for time-to-event pretraining
predictions. The resulting pretrained MOTOR model is then combined with a
linear head to train downstream, time-to-event models.

Figure 2: An overview of the MOTOR neural network architecture and the time-
to-event task used for pretraining.

####  3.2.1 Transformer

Following the previous work of BEHRT (Li et al., 2019), we use a Transformer
based neural network for processing our EHR data. $X_{i}$, the sequence of
codes and timestamps for a patient $i$, is first transformed into a
combination of code embeddings and time embeddings. Code embeddings are
generated with a standard embedding layer with a vocab size of 65,536 (the
maximum number of token ids that can fit in a 16 bit integer). Time embeddings
are generated using a standard sinusoidal position embedding with time (in
days since birth) instead of the position index being fed to the sine and
cosine functions. These embeddings are then input into a Transformer to
generate representations for every patient and every event in each patient’s
record. The overall Transformer backbone contains 100M trainable parameters.

Our Transformer has four major changes that distinguish it from Transformer
architectures commonly used with structured EHR data. First, unlike BEHRT
which uses a masked language modeling objective that incorporates
bidirectional information from patient sequences, we use causally masked
attention (Dai and Le, 2015; Peters et al., 2018; Raffel et al., 2022). This
attention pattern ensures information only flows forwards in the sequence and
not backwards to guarantee that the model does not “cheat” by using
information from the future to predict the past. If fully connected attention
is used instead, models quickly learn shortcuts based on future information
which causes issues when applying the model in the standard case when future
data is not available. Second, we use rotary position embeddings (Su et al.,
2021) for combining the time and code embeddings at every attention subunit in
the network. Third, due to the long sequence lengths in our data, we use a
local attention mechanism, where attention is limited to 512 token context
windows, in order to avoid quadratic complexity (Roy et al., 2021). This
enables us to use sequence lengths up to 16,384 events (the maximum we can fit
in 32 GB of GPU memory). Fourth, we do not the explicit visit parts of BEHRT
(the segment and position features) as our source EHR data contains
overlapping visits, which is not compatible with BEHRT’s approach. Instead, we
represent visits as additional clinical events within the timeline.

The output of our model is one matrix per patient $R_{i}$, where each row
$R_{ij}$ corresponds to the patient $i$’s representation for the event $j$
that contains cumulative information up to and including event $j$. The width
of the row $R_{ij}$ is the size of the patient representation, which is a
tunable hyperparameter.

####  3.2.2 Piecewise Exponential Time-to-Event Model

We choose a modeling strategy that is amenable to large-scale deep learning in
a setting with thousands of tasks and provides a differentiable objective that
can support high censoring rates while being computationally efficient.
Piecewise Exponential Artificial Neural Networks (PEANN) (Fornili et al.,
2014), which combines neural networks with piecewise exponential models, is
one such approach. Other approaches, like Cox based approaches, require at
least one uncensored example in each batch, which introduces additional
training complexity by requiring batch sampling strategies. Thus, we use a
PEANN-style approach to construct piecewise exponential models for each task,
indexed by $k$, with the hazards parameterized by neural networks.

Piecewise exponential time-to-event models require splitting the observed
follow-up times into a fixed number of time bins $B$ that cover the full
possible range of event times from zero to positive infinity. Piecewise
exponential models assume that the distribution of event times within each bin
follows an exponential distribution (which is equivalent to assuming that the
hazard rate is constant within each bin). Let $\lambda_{ijkb}$ be the constant
hazard rate for patient $i$ after event $j$ for task $k$ within the time bin
$b$. Note that the bin construction can be arbitrary, the only restriction is
that it bins must start from zero and end at positive infinity. In our
experiments, we select our bins by having the time points of each bin based on
percentiles of the event times. This ensures that each bin contains an equal
fraction (1/$B$) of the observed events. We tune the number of bins $B$ as a
hyperparameter.

With hazard $\lambda_{ijkb}$, it is then possible to define the standard
maximum likelihood time-to-event loss function. Let $\delta_{ijkb}$ denote
whether or not the $i$th subject at the time of the $j$th event has an event
for the $k$th task within the $b$th bin, and let
$U_{ijkb}=\mathrm{max}\\{0,\mathrm{min}(T_{ijkb}-T_{ijk(b-1)},U_{ijkb}-T_{ijk(b-1)})\\}$
indicate the observed time in that bin. Recall the density function for
exponential distribution is $f(U;\lambda)=\lambda\,\mathrm{exp}(-\lambda\,U)$,
then the likelihood for uncensored subjects is $f(U;\lambda)$ and for censored
subjects is $F^{c}(U;\lambda)=1-F(U;\lambda)$, where
$F(U;\lambda)=1-\mathrm{exp}(-\lambda\,U)$ is the cumulative distribution
function. Thus, the likelihood function is

$L(U_{ijkb}|\lambda_{ijkb})=\\{\lambda_{ijkb}\,\mathrm{exp}(-\lambda_{ijkb}\,U_{ijkb})\\}^{\delta_{ijkb}}\\{\mathrm{exp}(-\lambda_{ijkb}\,U_{ijkb})\\}^{1-\delta_{ijkb}}$
(2)

Note that the likelihood function is specific to a particular patient index,
day index, task index, and bin index. In order to obtain the likelihood
function on the entire dataset, it is necessary to take the product over every
patient, event, task and bin.

$L(U|\lambda)=\prod_{i,j,k,b}L(U_{ijkb}|\lambda_{ijkb})$ (3)

####  3.2.3 Estimating the Hazard Given A Patient Representation

Given a particular patient representation, we need a method for computing the
vector of hazard rates $\lambda_{ijkb}$ in Equation (2). We assume that the
transformation $R_{ij}\mapsto{\log(\lambda_{ijkb})}$ is a linear mapping of
low rank, and compute it with multiple smaller matrix multiplications to save
on compute, memory, and parameter counts.

First, we transform our patient representation $R_{ij}$ into a feature matrix
$X_{ijb}$ for modeling the hazard rates in each time bin $b$ with a learned
linear transformation. The feature matrix has an inner width of size $p$,
which is a tuned parameter that defines the dimensionality of our per-time bin
features. We can then compute $\log(\lambda_{ijbk})=X_{ijb}\hat{\beta}_{k}$,
where $\hat{\beta}_{k}$ is a vector of learned parameters with length $p$ for
a specific task $k$.

####  3.2.4 Time-to-Event Pretraining Task

Given our Transformer encoder, algorithm to compute hazard, and piecewise
exponential time-to-event model specification, we define a pretraining task of
predicting the time to various events defined by codes in our EHR data. To
train a general time-to-event model that captures as much signal as possible,
we train on as many tasks (codes) as we can comfortably fit into GPU memory,
$8,192$ in our setup. We select the $8,192$ most informative codes as ranked
by entropy, defined by Shannon’s formula (Shannon, 1948):

$H(code)=-\mathbb{P}(code)\,\mathrm{log}\,\mathbb{P}(code)-(1-\mathbb{P}(code))\,\mathrm{log}\,\mathbb{P}(code),$
(4)

where $code$ is a code and $\mathbb{P}(code)$ is the probability of code
$code$ occurring within a particular day.

In order to provide a more conservative measure of performance on novel time-
to-event tasks, we take an additional step of explicitly removing codes for
our evaluation target tasks from the set of pretraining tasks. Note that
removing codes for evaluation is not strictly necessary for correctness as we
are using a patient split to ensure no leakage between pretraining and
evaluation.

We pretrain our model end-to-end using the negative log-likelihood derived
based on Equation (2) averaging over all $8,192$ codes. We train for 10 epochs
using the Adam optimizer with the standard GPT2 (Radford et al., 2019)
learning rate schedule. The learning rate and other hyperparameters are set
through a grid search on the validation set. See Appendix Appendix A.
Hyperparameter Search Grids for the full details of our hyperparameter grid
search.

####  3.2.5 Adapting MOTOR to Target Tasks

Given the pretrained model above, we can then train models for new target
tasks of interest through the process of fine-tuning. Specifically, we
estimate a new vector of parameters $\hat{\beta}_{k}^{\mathrm{\,new}}$ by
minimizing the new loss function defined by the new task. In order to increase
the runtime efficiency and avoid overfitting, we only tune those new
$\hat{\beta}_{k}^{\mathrm{\,new}}$ parameters for each task and do not update
any other weights. This is commonly done with other pretrained models such as
BERT. The simpler setup of only training the new
$\hat{\beta}_{k}^{\mathrm{\,new}}$ parameters, as compared to full fine-tuning
where all the weights are updated, also allows us to use a more sophisticated
optimizer with faster convergence rates. To take advantage of this, we use a
conjugate gradient optimizer when fine-tuning. In order to improve
performance, we augment our conjugate gradient optimizer with an L2 penalty
term and use the warm start optimization. We sweep the L2 regularization
strength for each target task independently.

##  4 Experiments

We evaluate the effectiveness of the MOTOR method by performing experiments
where we compare its overall performance to several alternative time-to-event
modeling approaches and evaluate the impact of pretraining data size and fine-
tuning data size on model performance.

###  4.1 Data and Compute

We use the 2022 release of the de-identified STAnford Research data repository
(STARR) Observational Medical Outcomes Partnership (OMOP) (Observational
Health Data Sciences and Informatics, 2022) clinical records dataset. STARR
OMOP contains longitudinal EHR data for patients at Stanford Hospital and
Stanford Children’s Hospital and contains 2,730,411 patients (2,368,208,275
clinical events) with most data occurring between 2014 and 2022. Each patient
record consists of demographic information (age, sex, ethnicity) and coded
clinical information including diagnosis codes, lab test orders and results,
medication orders, procedures, and visits.

We process our dataset by converting all of the tables within the OMOP schema
to a uniform timeline format of timestamped events, where each event is
annotated with both the time and the clinical code for that event. We
additionally perform some post-processing to correct several timestamp-related
issues in the source data. In particular, we move billing codes to the end of
visits to better reflect when those billing codes are clinically available and
move all events annotated at 12:00 AM to 11:59 PM to better account for how
12:00 AM timestamps frequently represent day-level timestamps. We also adjust
for events before birth by moving them to the date of birth.

In order to obtain unbiased evaluation metrics for the performance of methods,
we split our data into training, validation, and test sets, containing 70%,
15% and 15% of the patients’ records, respectively. We perform our splitting
by patient id such that each patient is assigned a split. The training set
(1.9 million patients) is used for pretraining and fine-tuning target task
time-to-event models, the validation set (0.4 million patients) is used for
hyperparameter selection, and the test set (0.4 million patients) is used for
out-of-sample evaluation. Every task (including both pretraining tasks and
target tasks) uses the exact same data splits to ensure that there is no data
leakage.

Experiments are performed in Stanford’s NERO compute environment using 24
Intel Xeon 2.70GHz CPU cores and 8 Nvidia V100 GPUs. Each pretraining and/or
fine-tuning job is only assigned one V100 at a time.

###  4.2 Target Task Setup

From the STARR dataset, we derive six time-to-event target tasks, Celiac
Disease, Heart Attack, Lupus, Nonalcoholic Fatty Liver Disease (NAFLD),
Pancreatic Cancer, and Stroke. These tasks were chosen as they cover a wide
spectrum of clinically interesting acute and chronic health conditions. Each
time-to-event modeling task aims to predict the time to the clinical event of
interest from a random visit in the patient record after the patient’s first
year at Stanford Hospital. The clinical events are defined by a set of ICD
codes. As noted in Section 3.2.4, we remove these ICD codes from our
pretraining setup to better evaluate model performance on unseen tasks.

For each task, we define a set of cases and controls, where cases have the
event in question and controls do not. As most patients do not have any of the
diseases in our set of tasks, the number of controls vastly outnumbers the
number of cases. This excess of controls creates issues evaluating some of the
less computationally efficient baselines (e.g., Random Survival Forests), so
we employ a naive subsampling strategy where 4 controls are randomly selected
for every case. Note that this naive subsampling introduces a certain degree
of censoring bias in that hazard rates are artificially scaled up due to the
removal of censored examples. Appendix Appendix B. Subsampling Bias provides
the formula for the exact hazard rate scaling implied by our subsampling
strategy. However, this does not affect our primary results as the resulting
hazard rate scaling bias is shared across all methods and the primary goal of
this study is to evaluate the relative performance of different modeling
approaches. Since we evaluate and train all models using the same subsampling,
the censoring bias affects them in a similar manner and our ranking of models’
performance should not be impacted. We considered correcting for this bias
with a weighting strategy (Kohler et al., 2002; Van der Laan and Robins,
2003), but extreme weights induce unstable training. The set of ICD codes and
the number of cases and controls are detailed in Table 1.

Table 1: Task setup details for six target tasks. Task Name | ICD 10 Codes | Number of Cases | Number of Controls
---|---|---|---
Celiac Disease | K90.0 | 3453 | 13812
Heart Attack | I21.* | 2815 | 11260
Lupus | M32.* | 3274 | 13812
Pancreatic Cancer | C25.* | 2163 | 8652
Nonalcoholic Fatty Liver Disease | K76.0 | 22013 | 88052
Stroke | I63.* | 11855 | 47420

###  4.3 Methods under Comparison

We compare MOTOR with three state-of-the-art time-to-event modeling approaches
that make different types of modeling assumptions: Cox PH (Cox, 1972),
DeepSurv (Katzman et al., 2018), and Random Survival Forests (Ishwaran et al.,
2008b).

Cox PH and Random Survival Forests cannot naively handle our sequential event
format as they require input features to be in the form of a feature matrix.
In order to run those methods on our dataset, we implement a featurization
strategy to transform our timestamped event data into rows and columns that
can be used for those methods. Following prior standard practice in the
medical domain (OHDSI, 2019), we implement a count featurization strategy
where a column is created for each event code in the source data and each
value in that column corresponds to the number of times that event code
appears in the patient record. We augment this featurization strategy by
adding a hyperparameter for ontology expansion where you inject additional
columns for higher level concepts, such as “Cancer” for “Lung Cancer”. Higher
level concepts are determined using the SNOMED ontology (El-Sappagh et al.,
2018), which provides relationships and groupings for coded information. This
ontology expansion makes it easier to learn models in low data settings by
increasing the overlap between the features available for each patient (Choi
et al., 2016b).

In principle, DeepSurv can accept the same timestamped input format as MOTOR,
however we observed poor performance due to overfitting of our larger
parameterized model on our small time-to-event datasets. Instead, we train
DeepSurv on the same feature matrices used for Random Survival Forests and Cox
PH. These models are less prone to overfitting because they do not require
estimating as many parameters as we can use simpler non-sequence oriented
neural network models. We parameterize DeepSurv with a three layer feedforward
neural network with a rectified linear unit (ReLU) activation function and
batch-normalization. As in MOTOR, we use dropout and early stopping for
regularization.

Appendix Appendix A. Hyperparameter Search Grids contains the full
hyperparameter search grids for all approaches as well as the software
versions used when applicable.

###  4.4 Experiment Schema

We focus on evaluating the model performance of MOTOR, as well as the methods
in Section 4.3, in the following aspects:

  1. a

Overall Performance: We compare MOTOR with three baseline time-to-event
methods using all available target task data.

  2. b

Performance of Pretraining Task: We examine the performance of MOTOR on the
pretraining task.

  3. c

Choice of Pretraining Task: We compare our time-to-event pretraining task to
an auto-regressive, next code pretraining task.

  4. d

Sample Size Used in Pretraining: We examine the impact of pretraining on model
performance by conducting pretraining on a decreasing proportion of the
training data. We evaluate using 5%, 10%, and 25% of our training split for
pretraining.

  5. e

Sample Size Used in Fine-tuning: We examine the impact of the fine-tuning
sample size on model performance. We use all the training samples for
pretraining but randomly sample 5%, 10% and 25% of the training split for
fine-tuning.

  6. f

Time-to-Train: We compare the implementation wall time of MOTOR to that of the
second most performant predictive method, Random Survival Forests.

###  4.5 Evaluation Metrics

We evaluate model performance using two metrics, the time-dependent C
statistic and the Nam-D’Agostino (ND) calibration metric (D’Agostino and Nam,
2003). These metrics were chosen because they can be used to evaluate time-
dependent risk (which cannot be done with time-independent metrics like
Harrell’s C statistic).

The time-dependent C statistic summarizes the ability of a predictive model to
rank patients by risk based on their instantaneous hazard functions.
Specifically, we apply the incident and dynamic definitions for the time-
dependent sensitivity and specificity, respectively, following the choices in
Heagerty and Zheng (2005):

$\displaystyle\mathrm{sen}^{\mathbb{I}}(c,t)$
$\displaystyle=\mathbb{P}(S_{i}>c|\mathrm{d}N_{i}^{*}(t)=1)$
$\displaystyle\mathrm{spe}^{\mathbb{D}}(c,t)$
$\displaystyle=\mathbb{P}(S_{i}\leq c|N_{i}^{*}(t)=0),$

where $S_{i}$ is the predicted risk score for subject $i$.
$N_{i}^{*}(t)=\mathbb{1}\\{T_{i}\leq t\\}$ and
$\mathrm{d}N_{i}^{*}(t)=N_{i}^{*}(t)-N_{i}^{*}(t-)$ follow the standard
definitions of a counting process. Then, the concordance can be expressed as a
weighted average of the area under the time-dependent receiver operating curve
(ROC)

$C=\frac{\int_{t}\mathrm{AUC}(t)\,w(t)\,\mathrm{d}t}{\int_{t}w(t)\,\mathrm{d}t}$
(5)

where $w(t)=f(t)\cdot S(t)$. $f(t)$ is the event rate at a particular time t
and $S(t)$ is the survival probability at that time. Both $f(t)$ and $S(t)$
are estimated using the Kaplan-Meier estimator. $\mathrm{AUC}(t)$ is the time-
specific area under ROC and can be computed by integration, and ROC can be
calculated using $\mathrm{sen}^{\mathbb{I}}$ and $\mathrm{spe}^{\mathbb{D}}$
with the standard formulation.

Note that the formula in Equation (5) is slightly modified from Heagerty and
Zheng (2005) because our dataset contains many events that happened at the
exact same time due to some events having only day-level time resolution.
Heagerty and Zheng (2005) assumes that all event times are distinct thus the
sum of the weights (the term in the denominator) is $\frac{1}{2}$. As we have
many tied event times, we cannot make the same assumption and have to
explicitly compute the denominator as $\int_{t}w(t)\,\mathrm{d}t$.

The ND calibration metric (D’Agostino and Nam, 2003) is constructed based on
the idea that the survival probability estimates, $\hat{S}(t)=P(T_{i}>t)$, for
a particular time $t$ should be calibrated with respect to the number of
observed survivors past that time point, which is not available in censored
data, thus we estimate it using the Kaplan-Meier estimator per D’Agostino and
Nam (2003). First, the data is sorted by $\hat{S}(t)$ given by a model and
split into $M$ equal size risk bins, where $\overline{p(t)}_{m}$ is the
average value of $\hat{S}(t)$ within bin $m$. Second, a KM estimator is fit
within each risk bin for computing actual survival probabilities. Finally, the
squared difference between the expected and actual values is computed and
summed across bins. The resulting statistic is $\chi^{2}$ with $M-1$ degrees
of freedom. In principle, you can run hypotheses tests using this statistic,
but we instead directly report and compare the raw statistic

$\chi^{2}(t)=\sum_{m=1}^{M}\frac{[KM(t)_{m}-\overline{p(t)}_{m}]^{2}\,n_{m}}{\overline{p(t)}_{m}(1-\overline{p(t)}_{m})},$
(6)

where $n_{m}$ is the number of observations in bin $m$. For all evaluations,
we set $t$ to the median event time. In addition, we make one minor
modification to this metric in that we remove the $n_{m}$ term in the
numerator as we are primarily concerned about relative performance. It’s a
constant factor that only depends on the number of patients for each task, and
removing it makes it easier to compare results across target tasks.

The time-dependent C statistic and ND test become increasingly unstable as the
censoring rate increases since they rely on Kaplan Meier estimators that
themselves become unstable at higher censoring rates. As such, it is generally
recommended to only compute these statistics on a subset of the possible times
to reduce the variance of estimates. As suggested by Heagerty and Zheng
(2005), we use the time-restricted variants of these evaluation metrics by
setting a time limit of the 90% percentile of observed event times and only
perform evaluations up to that time point.

In order to accurately quantify statistical significance, we use the paired
test-set bootstrap approach (Konietschke and Pauly, 2013) to estimate the 95%
confidence interval for the difference in model performance between MOTOR and
every other method. For every target task, we resample our test set with
replacement to obtain 1000 bootstrap samples. We then compute all metrics for
all of our models on each bootstrap sample, subtracting MOTOR’s performance
each time to obtain an estimate of how each method differs from MOTOR. We
finally compute confidence intervals by extracting the 2.5% and 97.5%
percentiles of the bootstrapped samples. All of our confidence intervals for
the differences in performance between methods are reported in Appendix
Appendix C. Confidence Intervals.

##  5 Results

###  5.1 Overall Performance

Table 2 provides the time-dependent concordance performance of the various
methods trained on the full dataset. MOTOR performs substantially better than
all of our baselines, providing statistically significant higher C statistics
in all of our six tasks (Table S2). In addition, we see approaches that make
stronger assumptions show worse ranking performance under complicated data
structures. The Cox model results in the lowest C statistics since it assumes
both proportional hazard and a linear relationship between the log hazard and
covariates. Both DeepSurv (which does not make this linear assumption) and
Random Survival Forests (which is completely nonparametric) show improved C
statistics.

We also report a particularly important ablation of MOTOR vs MOTOR without any
pretraining (MOTOR-WP). As MOTOR combines a sophisticated neural network model
with pretraining, it’s important to demonstrate that pretraining itself
provides value. MOTOR-WP uses the exact same piecewise exponential setup,
Transformer backbone, and the hyperparameter grid as MOTOR, with the only
change being that MOTOR-WP is only trained on the target tasks. MOTOR-WP
performs worse than MOTOR, demonstrating the importance of pretraining for
modeling time-to-event tasks.

Table 2: Time-dependent C statistic for the various methods on the full dataset. Larger values indicate better ranking performance. Bolding indicates the best performing model. RSF = Random Survival Forests. MOTOR-WP = MOTOR Without Pretraining Method | Celiac | Heart Attack | Lupus | NAFLD | Pancretic Cancer | Stroke
---|---|---|---|---|---|---
Cox PH | 0.689 | 0.761 | 0.770 | 0.726 | 0.793 | 0.779
DeepSurv | 0.704 | 0.823 | 0.790 | 0.800 | 0.811 | 0.830
RSF | 0.729 | 0.836 | 0.787 | 0.802 | 0.824 | 0.840
MOTOR-WP | 0.696 | 0.795 | 0.803 | 0.821 | 0.777 | 0.831
MOTOR | 0.802 | 0.884 | 0.850 | 0.859 | 0.865 | 0.874

Table 3 provides the corresponding overall ND calibration results. MOTOR
yields the smallest (modified) ND statistics for most of the target tasks;
however, the improvements are not statistically significant given the
overlapping confidence intervals in Table S3. Compared to the Cox PH model,
the nonlinear DeepSurv results in smaller ND statistics but only shows
statistically significant better performance for the Lupus task. Of note is
that the MOTOR-WP baseline has substantially worse calibration, which is a
common failure mode for neural network models (Guo et al., 2017).

Table 3: ND calibration for the various methods on the full dataset. Smaller values indicate smaller calibration errors. Bolding indicates the best performing model. RSF = Random Survival Forests. MOTOR-WP = Motor Without Pretraining. Method | Celiac | Heart Attack | Lupus | NAFLD | Pancretic Cancer | Stroke
---|---|---|---|---|---|---
Cox PH | 0.167 | 0.214 | 0.373 | 0.121 | 0.141 | 0.031
DeepSurv | 0.041 | 0.089 | 0.102 | 2.911 | 0.081 | 0.023
RSF | 0.299 | 0.225 | 0.095 | 0.091 | 0.345 | 0.112
MOTOR-WP | 4.970 | 3.990 | 23.882 | 0.077 | 7.104 | 0.083
MOTOR | 0.045 | 0.029 | 0.059 | 0.019 | 0.046 | 0.032

###  5.2 Performance of Pretraining Task

When MOTOR is being pretrained, we are implicitly training 8,192 piecewise
exponential time-to-event models across a variety of medical codes. Assessing
pretraining performance is important because it allows us to evaluate MOTOR’s
performance on a wider variety of tasks with vastly different censoring rates
than our subsampled target tasks. Figure 3 provides the C statistic as a
function of the amount of uncensored data for each task, with lines indicating
the C statistics of the first quartile (0.783), median (0.840), and third
quartile (0.898).

Performance on the pretraining tasks is strong across the vast majority of
tasks, as indicated by the relatively high first quartile C statistic of
0.783. It’s noteworthy that performance is strong even for highly censored
tasks, even when only one in a thousand patients have an uncensored event.

Figure 3: Pretraining task C statistic as a function of the fraction of
uncensored data. The dotted lines indicate the quartiles for the C statistic
across all fractions of uncensored data.

###  5.3 Choice of Pretraining Task

One of the main novel components of MOTOR is our time-to-event pretraining
task. In order to investigate how much this pretraining task contributes to
MOTOR’s performance, we perform a comparison with another common EHR
pretraining task: next code prediction. This autoregressive pretraining task
aims to predict the next code in an EHR sequence (Steinberg et al., 2021). One
limitation of the next code prediction setup is that we are unable to perform
a straightforward linear head fine-tuning step because the next code
prediction pretraining task does not contain the necessary linear
transformation from $R_{ij}$ to $X_{ijb}$ in our piecewise exponential model.

Instead, for this experiment, we perform full weight updates while fine-
tuning, updating every weight in the model using gradient descent. In order to
ensure a fair comparison, we apply the same tweak to MOTOR and then compare
performance on our target tasks. We thus ensure that both pretraining tasks
have a matched setup, with the same architecture, hyperparameter search grid,
and fine-tuning strategy.

Table 4 compares the next code pretraining task to our time-to-event
pretraining task. We see that the time-to-event pretraining task outperforms
the next code pretraining task, indicated by the statistically significant
higher C statistics on most of the target tasks (Table S4). This implies that
utilizing a time-to-event pretraining task is beneficial for predicting time-
to-event outcomes.

Table 4: Impact of the pretraining objective on the performance of MOTOR in terms of the time-dependent C statistic. We compare a “next code” baseline with our proposed time-to-event objective. Note that these experiments are conducted with full fine-tuning of the entire model. Objective | Celiac | Heart Attack | Lupus | NAFLD | Pancretic Cancer | Stroke
---|---|---|---|---|---|---
Next Code | 0.774 | 0.862 | 0.842 | 0.860 | 0.860 | 0.857
Time-to-Event | 0.800 | 0.881 | 0.858 | 0.864 | 0.861 | 0.874

###  5.4 Sample Size Used in Pretraining

As MOTOR relies on pretraining to achieve good performance, an important
question is how MOTOR’s performance changes as a function of available
pretraining data. We evaluate 5% to 25% uniform random samples of our
pretraining dataset and report performance for all target tasks. Table 5 shows
that performance degrades quite rapidly as the pretraining data size
decreases, quickly approaching the performance of the baselines for smaller
pretraining sizes. The decreased performance of MOTOR on these subsampled
datasets likely reflects overfitting that occurs when training our high-
parameter model on relatively small time-to-event datasets.

Table 5: Impact of pretraining on the performance of MOTOR in terms of the time-dependent C statistic. We utilize a subset of the full training data to conduct pretraining. K = thousand patients. M = million patients. # Patients For Pretraining | Celiac | Heart Attack | Lupus | NAFLD | Pancretic Cancer | Stroke
---|---|---|---|---|---|---
95K (5%) | 0.694 | 0.843 | 0.740 | 0.768 | 0.825 | 0.824
191K (10%) | 0.714 | 0.849 | 0.749 | 0.786 | 0.831 | 0.828
478K (25%) | 0.742 | 0.853 | 0.788 | 0.799 | 0.837 | 0.833
1.91M (100%) | 0.802 | 0.884 | 0.850 | 0.859 | 0.865 | 0.874

###  5.5 Sample Size Used for Fine-tuning

A key benefit of foundation models is improved sample efficiency, enabling the
adaptation of pretrained models to novel tasks while requiring fewer labeled
training examples. We evaluate MOTOR’s sample efficiency by randomly
subsampling the datasets for our target tasks and fine-tuning MOTOR on these
smaller datasets. We evaluate on 5%, 10% and 25% and 100% of the original
fine-tuning datasets.

Table 6 contains the time-dependent C statistics for each of those reduced
training sets. Even though decreasing fine-tuning sample sizes induces lower C
statistics (Table S6), MOTOR still has strong ranking performance even with
only 5% of the original training data for fine-tuning. This probably reflects
that there are many similar tasks within the pretraining set, enabling
efficient fine-tuning even given limited available task data for training.

Table 6: Impact of sample sizes for fine-tuning MOTOR on the time-dependent C statistic. We utilize a subset of the full training data to conduct pretraining. 100% corresponds to the entire training dataset. % Used | Celiac | Heart Attack | Lupus | NAFLD | Pancretic Cancer | Stroke
---|---|---|---|---|---|---
5% | 0.785 | 0.878 | 0.840 | 0.851 | 0.849 | 0.868
10% | 0.781 | 0.874 | 0.844 | 0.854 | 0.854 | 0.868
25% | 0.790 | 0.880 | 0.845 | 0.856 | 0.859 | 0.869
100% | 0.802 | 0.884 | 0.850 | 0.859 | 0.865 | 0.874

###  5.6 Time-to-Train

The ability to fit time-to-event models quickly is essential for a variety of
scenarios such as interactive model training and limited resource
environments. One of the main advantages of MOTOR is that almost all of the
weights can be frozen during fine-tuning, which means that one can learn
target task models very quickly by only computing gradients for the linear
heads.

In order to evaluate model training time, we perform timing experiments on all
of our six target tasks, comparing MOTOR to the current state-of-the-art,
Random Survival Forests. To simplify our experimental setup, we ignore the
cost of hyperparameter tuning and focus solely on the time to fit the final
target task models with the optimally determined hyperparameters. We provide
16 CPU cores, one V100 GPU, and 100 GB of RAM on isolated machines to both
Random Survival Forests and MOTOR.

Table 7 contains the clock time, in minutes for both methods on each task. The
main variable which affects the runtime is the number of patients per task,
which is described in Table 1. Note that Random Survival Forests is more
efficient on tasks with fewer patients, where the overhead of MOTOR is more
significant. However, this completely flips for the tasks with many patients,
in particular NAFLD and Stroke, where the increased efficiency of MOTOR allows
us to train target task models more quickly than Random Survival Forests, with
an 8.2x increase in speed for the NAFLD task. We note this is a somewhat
conservative time comparison that ignores the cost and additional workflow
complexity of materializing Random Survival Forests’ feature matrix.

Table 7: The time it takes to train a target task model, in minutes. Each run is given 16 CPU cores, one V100 GPU, and 100 GB of RAM. The fastest method for each target task is bolded. RSF = Random Survival Forests. Method | Celiac | Heart Attack | Lupus | NAFLD | Pancretic Cancer | Stroke
---|---|---|---|---|---|---
RSF | 1.82 | 1.21 | 5.80 | 93.31 | 2.08 | 22.22
MOTOR | 2.35 | 2.47 | 2.40 | 11.45 | 2.16 | 7.50

##  6 Discussion

We proposed a novel and effective approach, which we call MOTOR, for learning
time-to-event models using transfer learning. The use of transfer learning
allows MOTOR to mitigate the key issue of overfitting in standard time-to-
event analysis approaches, which allows MOTOR to use more sophisticated high-
parameter models even in limited data settings. From conducting a
comprehensive experiment, we show that MOTOR outperforms several state-of-the-
art approaches and provides estimators with strong ranking and calibration
performance across a range of clinical tasks. Such success is attributed to
the novel time-to-event pretraining process that enables MOTOR to learn
sophisticated interactions and event distributions over 8,192 informative
events.

In addition, MOTOR is scalable and computationally efficient. We only need to
pretrain MOTOR once on the entire dataset, and then the same pretrained model
can be reused many times on an arbitrary number of target tasks. As this fine-
tuning only requires updating a small number of weights, it is very
computationally efficient, allowing us to train task-specific models up to
eight times faster than Random Survival Forests. The decreased computation
costs of MOTOR make time-to-event models much more practical in clinical
settings with limited computational resources.

This work has several limitations. First, pretraining MOTOR is quite
expensive, often requiring several GPU days in our experiments. Even though
this is a one-time cost, it might be prohibitive in some severely resource-
constrained settings. Second, our experimental protocol performed evaluations
on subsampled target task datasets to allow us to evaluate more methods. Doing
so introduces a certain degree of censoring bias, which might affect our
results. Third, we only evaluated MOTOR on ICD code-derived tasks to simplify
the labeling effort. In some cases, assigning labels to specific medical
conditions or outcomes requires developing complicated phenotyping rules and
machine learning-based methods (Banda et al., 2018). We leave evaluating MOTOR
in these settings for future work.

##  7 Conclusion

In conclusion, we validated that MOTOR is an effective approach for learning
reliable time-to-event models for a variety of tasks, even in the presence of
limited data. The use of pretraining, and in particular the specialized time-
to-event pretraining task, allows us to learn flexible deep neural network
time-to-event models in low sample settings and low compute resource
situations where they would otherwise not be viable.

## Acknowledgements

This work was supported by R01 HL144555 from the National Heart, Lung, and
Blood Institute (NHLBI).

## References

  * Athey et al. [2019] S. Athey, J. Tibshirani, and S. Wager.  Generalized random forests.  _Ann. Statist._ , pages 47(2): 1148–1178, 2019.
  * Avati et al. [2020] A. Avati, T. Duan, S. Zhou, K. Jung, N. H. Shah, and A. Y. Ng.  Countdown regression: sharp and calibrated survival predictions.  In _Uncertainty in Artificial Intelligence_ , pages 145–155. PMLR, 2020.
  * Banda et al. [2018] J. M. Banda, M. Seneviratne, T. Hernandez-Boussard, and N. H. Shah.  Advances in electronic phenotyping: from rule-based definitions to machine learning models.  _Annual review of biomedical data science_ , 1:53, 2018\.
  * Bommasani et al. [2021] R. Bommasani, D. A. Hudson, E. Adeli, R. Altman, S. Arora, S. von Arx, M. S. Bernstein, J. Bohg, A. Bosselut, E. Brunskill, et al.  On the opportunities and risks of foundation models.  _arXiv preprint arXiv:2108.07258_ , 2021.
  * Breiman [2001] L. Breiman.  Random forests.  _Machine Learning_ , pages 45: 5–32, 2001.
  * Chang et al. [2022] T. Chang, M. W. Sjoding, and J. Wiens.  Disparate censorship & undertesting: A source of label bias in clinical machine learning, 2022.  URL https://arxiv.org/abs/2208.01127.
  * Che et al. [2016] Z. Che, S. Purushotham, K. Cho, D. Sontag, and Y. Liu.  Recurrent neural networks for multivariate time series with missing values, 2016.  URL https://arxiv.org/abs/1606.01865.
  * Choi et al. [2016a] E. Choi, M. T. Bahadori, J. A. Kulas, A. Schuetz, W. F. Stewart, and J. Sun.  Retain: An interpretable predictive model for healthcare using reverse time attention mechanism, 2016a.  URL https://arxiv.org/abs/1608.05745.
  * Choi et al. [2016b] E. Choi, M. T. Bahadori, L. Song, W. F. Stewart, and J. Sun.  Gram: Graph-based attention model for healthcare representation learning, 2016b.  URL https://arxiv.org/abs/1611.07012.
  * Cox [1972] D. R. Cox.  Regression models and life-tables.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 34(2):187–202, 1972.
  * D’Agostino and Nam [2003] R. D’Agostino and B.-H. Nam.  Evaluation of the performance of survival analysis models: Discrimination and calibration measures.  In _Advances in Survival Analysis_ , volume 23 of _Handbook of Statistics_ , pages 1–25. Elsevier, 2003.  doi: https://doi.org/10.1016/S0169-7161(03)23001-7.  URL https://www.sciencedirect.com/science/article/pii/S0169716103230017.
  * Dai and Le [2015] A. M. Dai and Q. V. Le.  Semi-supervised sequence learning.  In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 28. Curran Associates, Inc., 2015.  URL https://proceedings.neurips.cc/paper/2015/file/7137debd45ae4d0ab9aa953017286b20-Paper.pdf.
  * Dotson [2013] P. Dotson.  Codes: What Are They, Why Are They Necessary, and How Are They Developed?  _Adv Wound Care (New Rochelle)_ , 2(10):583–587, Dec 2013.
  * El-Sappagh et al. [2018] S. El-Sappagh, F. Franda, F. Ali, and K. S. Kwak.  SNOMED CT standard ontology based on the ontology for general medical science.  _BMC Med Inform Decis Mak_ , 18(1):76, Aug 2018\.
  * Fornili et al. [2014] M. Fornili, F. Ambrogi, P. Boracchi, and E. Biganzoli.  Piecewise exponential artificial neural networks (peann) for modeling hazard function with right censored data.  In E. Formenti, R. Tagliaferri, and E. Wit, editors, _Computational Intelligence Methods for Bioinformatics and Biostatistics_ , pages 125–136, Cham, 2014. Springer International Publishing.  ISBN 978-3-319-09042-9.
  * Friedman [1982] M. Friedman.  Piecewise exponential models for survival data with covariates.  _The Annals of Statistics_ , 10(1):101–113, 1982\.
  * Gao and Cui [2022] Y. Gao and Y. Cui.  Clinical time-to-event prediction enhanced by incorporating compatible related outcomes.  _PLOS Digit Health_ , 1(5), 2022.
  * Goff et al. [2014] D. Goff, D. Lloyd-Jones, G. Bennett, S. Coady, R. D’Agostino, R. Gibbons, P. Greenland, D. Lackland, D. Levy, C. O’Donnell, et al.  American college of cardiology/american heart association task force on practice guidelines. 2013 acc/aha guideline on the assessment of cardiovascular risk: a report of the american college of cardiology/american heart association task force on practice guidelines.  _Circulation_ , 129(25 suppl 2):49–73, 2014.
  * Guo et al. [2017] C. Guo, G. Pleiss, Y. Sun, and K. Q. Weinberger.  On calibration of modern neural networks.  _CoRR_ , abs/1706.04599, 2017.  URL http://arxiv.org/abs/1706.04599.
  * Heagerty and Zheng [2005] P. J. Heagerty and Y. Zheng.  Survival model predictive accuracy and roc curves.  _Biometrics_ , 61(1):92–105, 2005.
  * Ishwaran et al. [2008a] H. Ishwaran, U. Kogalur, E. Blackstone, and et al.  Random survival forests.  _Ann Appl Stat_ , pages 2(3):841–860, 2008a.
  * Ishwaran et al. [2008b] H. Ishwaran, U. B. Kogalur, E. H. Blackstone, and M. S. Lauer.  Random survival forests.  _The Annals of Applied Statistics_ , 2(3), sep 2008b.  doi: 10.1214/08-aoas169.  URL https://doi.org/10.1214%2F08-aoas169.
  * Jing and Smola [2017] H. Jing and A. J. Smola.  Neural survival recommender.  In _Proceedings of the Tenth ACM International Conference on Web Search and Data Mining_ , WSDM ’17, page 515–524, New York, NY, USA, 2017. Association for Computing Machinery.  ISBN 9781450346757.  doi: 10.1145/3018661.3018719.  URL https://doi.org/10.1145/3018661.3018719.
  * Katzman et al. [2018] J. L. Katzman, U. Shaham, A. Cloninger, J. Bates, T. Jiang, and Y. Kluger.  DeepSurv: personalized treatment recommender system using a cox proportional hazards deep neural network.  _BMC Medical Research Methodology_ , 18(1), feb 2018.  doi: 10.1186/s12874-018-0482-1.  URL https://doi.org/10.1186%2Fs12874-018-0482-1.
  * Kohler et al. [2002] M. Kohler, K. Máthé, and M. Pintér.  Prediction from randomly right censored data.  _Journal of Multivariate Analysis_ , 80(1):73–100, 2002.
  * Konietschke and Pauly [2013] F. Konietschke and M. Pauly.  Bootstrapping and permuting paired t-test type statistics.  _Statistics and Computing_ , 24(3):283–296, Jan. 2013.  doi: 10.1007/s11222-012-9370-4.  URL https://doi.org/10.1007/s11222-012-9370-4.
  * Kopper et al. [2020] P. Kopper, S. Pölsterl, C. Wachinger, B. Bischl, A. Bender, and D. Rügamer.  Semi-structured deep piecewise exponential models.  _CoRR_ , abs/2011.05824, 2020.  URL https://arxiv.org/abs/2011.05824.
  * Li et al. [2016] Y. Li, J. Wang, J. Ye, and C. K. Reddy.  A multi-task learning formulation for survival analysis.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, page 1715–1724, New York, NY, USA, 2016. Association for Computing Machinery.  ISBN 9781450342322.  doi: 10.1145/2939672.2939857.  URL https://doi.org/10.1145/2939672.2939857.
  * Li et al. [2019] Y. Li, S. Rao, J. R. A. Solares, A. Hassaine, D. Canoy, Y. Zhu, K. Rahimi, and G. Salimi-Khorshidi.  Behrt: Transformer for electronic health records, 2019.  URL https://arxiv.org/abs/1907.09538.
  * Martinsson [2016] E. Martinsson.  _Wtte-rnn: Weibull time to event recurrent neural network_.  PhD thesis, Chalmers University of Technology & University of Gothenburg, 2016.
  * Miotto et al. [2016] R. Miotto, L. Li, B. A. Kidd, and J. T. Dudley.  Deep Patient: An Unsupervised Representation to Predict the Future of Patients from the Electronic Health Records.  _Sci Rep_ , 6:26094, May 2016.
  * Nelson et al. [2011] S. J. Nelson, K. Zeng, J. Kilbourne, T. Powell, and R. Moore.  Normalized names for clinical drugs: RxNorm at 6 years.  _J Am Med Inform Assoc_ , 18(4):441–448, 2011\.
  * Observational Health Data Sciences and Informatics [2022] Observational Health Data Sciences and Informatics.  Omop common data model.  https://ohdsi.github.io/CommonDataModel/index.html, 2022.
  * OHDSI [2019] OHDSI.  _The Book of OHDSI: Observational Health Data Sciences and Informatics_.  OHDSI, 2019.  ISBN 9781088855195.  URL https://books.google.com/books?id=JxpnzQEACAAJ.
  * Organization [2004] W. H. Organization.  Icd-10 : international statistical classification of diseases and related health problems : tenth revision, 2004.
  * Pang et al. [2021] C. Pang, X. Jiang, K. S. Kalluri, M. Spotnitz, R. Chen, A. Perotte, and K. Natarajan.  Cehr-bert: Incorporating temporal information from structured ehr data to improve prediction tasks.  In S. Roy, S. Pfohl, E. Rocheteau, G. A. Tadesse, L. Oala, F. Falck, Y. Zhou, L. Shen, G. Zamzmi, P. Mugambi, A. Zirikly, M. B. A. McDermott, and E. Alsentzer, editors, _Proceedings of Machine Learning for Health_ , volume 158 of _Proceedings of Machine Learning Research_ , pages 239–260. PMLR, 04 Dec 2021.  URL https://proceedings.mlr.press/v158/pang21a.html.
  * Perianez et al. [2016] A. Perianez, A. Saas, A. Guitart, and C. Magne.  Churn prediction in mobile social games: Towards a complete assessment using survival ensembles.  In _2016 IEEE International Conference on Data Science and Advanced Analytics (DSAA)_. IEEE, oct 2016.  doi: 10.1109/dsaa.2016.84.  URL https://doi.org/10.1109%2Fdsaa.2016.84.
  * Peters et al. [2018] M. E. Peters, M. Neumann, M. Iyyer, M. Gardner, C. Clark, K. Lee, and L. Zettlemoyer.  Deep contextualized word representations.  In _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers)_ , pages 2227–2237, New Orleans, Louisiana, June 2018. Association for Computational Linguistics.  doi: 10.18653/v1/N18-1202.  URL https://aclanthology.org/N18-1202.
  * Radford et al. [2019] A. Radford, J. Wu, R. Child, D. Luan, D. Amodei, and I. Sutskever.  Language models are unsupervised multitask learners.  2019\.
  * Raffel et al. [2022] C. Raffel, N. Shazeer, A. Roberts, K. Lee, S. Narang, M. Matena, Y. Zhou, W. Li, and P. J. Liu.  Exploring the limits of transfer learning with a unified text-to-text transformer.  _J. Mach. Learn. Res._ , 21(1), jun 2022.  ISSN 1532-4435.
  * Rajkomar et al. [2018] A. R. Rajkomar, E. Oren, K. Chen, A. Dai, N. Hajaj, M. Hardt, P. J. Liu, X. Liu, J. Marcus, M. Sun, P. P. Sundberg, H. Yee, K. Zhang, Y. Zhang, G. Flores, G. Duggan, J. Irvine, Q. Le, K. Litsch, A. Mossin, J. J. Tansuwan, D. Wang, J. Wexler, J. Wilson, D. Ludwig, S. Volchenboum, K. Chou, M. Pearson, S. Madabushi, N. Shah, A. Butte, M. Howell, C. Cui, G. Corrado, and J. Dean.  Scalable and accurate deep learning for electronic health records.  _npj Digital Medicine_ , 2018.  URL https://www.nature.com/articles/s41746-018-0029-1.
  * Rasmy et al. [2021a] L. Rasmy, M. Nigo, B. S. Kannadath, Z. Xie, B. Mao, K. Patel, Y. Zhou, W. Zhang, A. Ross, H. Xu, and D. Zhi.  Covrnn—a recurrent neural network model for predicting outcomes of covid-19 patients: model development and validation using ehr data.  _medRxiv_ , 2021a.  doi: 10.1101/2021.09.27.21264121.  URL https://www.medrxiv.org/content/early/2021/09/29/2021.09.27.21264121.
  * Rasmy et al. [2021b] L. Rasmy, Y. Xiang, Z. Xie, C. Tao, and D. Zhi.  Med-BERT: pretrained contextualized embeddings on large-scale structured electronic health records for disease prediction.  _npj Digital Medicine_ , 4(1), May 2021b.  doi: 10.1038/s41746-021-00455-y.  URL https://doi.org/10.1038/s41746-021-00455-y.
  * Ren et al. [2018] K. Ren, J. Qin, L. Zheng, Z. Yang, W. Zhang, L. Qiu, and Y. Yu.  Deep recurrent survival analysis, 2018.  URL https://arxiv.org/abs/1809.02403.
  * Roy et al. [2021] A. Roy, M. Saffar, A. Vaswani, and D. Grangier.  Efficient content-based sparse attention with routing transformers.  _Transactions of the Association for Computational Linguistics_ , 9:53–68, 2021.  doi: 10.1162/tacl˙a˙00353.  URL https://aclanthology.org/2021.tacl-1.4.
  * Ryu et al. [2020] J. Y. Ryu, M. Y. Lee, J. H. Lee, B. H. Lee, and K. S. Oh.  DeepHIT: a deep learning framework for prediction of hERG-induced cardiotoxicity.  _Bioinformatics_ , 36(10):3049–3055, May 2020\.
  * Shannon [1948] C. E. Shannon.  A mathematical theory of communication.  _The Bell System Technical Journal_ , 27:379–423, 1948\.  URL http://plan9.bell-labs.com/cm/ms/what/shannonday/shannon1948.pdf.
  * Steinberg et al. [2021] E. Steinberg, K. Jung, J. A. Fries, C. K. Corbin, S. R. Pfohl, and N. H. Shah.  Language models are an effective representation learning technique for electronic health record data.  _Journal of Biomedical Informatics_ , 113:103637, 2021.
  * Su et al. [2021] J. Su, Y. Lu, S. Pan, A. Murtadha, B. Wen, and Y. Liu.  Roformer: Enhanced transformer with rotary position embedding, 2021.  URL https://arxiv.org/abs/2104.09864.
  * Van der Laan and Robins [2003] M. J. Van der Laan and J. M. Robins.  _Unified methods for censored longitudinal data and causality_ , volume 5.  Springer, 2003.
  * Vaswani et al. [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. u. Kaiser, and I. Polosukhin.  Attention is all you need.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 30. Curran Associates, Inc., 2017.  URL https://proceedings.neurips.cc/paper/2017/file/3f5ee243547dee91fbd053c1c4a845aa-Paper.pdf.
  * Wang and Sun [2021] Z. Wang and J. Sun.  Survtrace: Transformers for survival analysis with competing events.  _CoRR_ , abs/2110.00855, 2021.  URL https://arxiv.org/abs/2110.00855.
  * Wen et al. [2022] Y. Wen, M. F. Rahman, Y. Zhuang, M. Pokojovy, H. Xu, P. McCaffrey, A. Vo, E. Walser, S. Moen, and T. B. Tseng.  Time-to-event modeling for hospital length of stay prediction for COVID-19 patients.  _Mach Learn Appl_ , 9:100365, Sep 2022.
  * Wu et al. [2015] W. C.-H. Wu, M.-Y. Yeh, and M.-S. Chen.  Predicting winning price in real time bidding with censored data.  In _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’15, page 1305–1314, New York, NY, USA, 2015. Association for Computing Machinery.  ISBN 9781450336642.  doi: 10.1145/2783258.2783276.  URL https://doi.org/10.1145/2783258.2783276.
  * Zeng et al. [2022] X. Zeng, S. L. Linwood, and C. Liu.  Pretrained transformer framework on pediatric claims data for population specific tasks.  _Scientific Reports_ , 12(1):3651, Mar 2022.  ISSN 2045-2322.  doi: 10.1038/s41598-022-07545-1.  URL https://doi.org/10.1038/s41598-022-07545-1.

## Appendix A. Hyperparameter Search Grids

Hyperparameters | Values
---|---
Count Based Representation
min_patients | 100, 1000
ontology_expansion | false, true
Cox PH
lambda | automatic
software_version | glmnet 4.1-6
DeepSurv
dropout | 0, 0.1
inner_layer_size | 32, 128, 256, 512
learning_rate |  $10^{-3}$, $10^{-4}$, $10^{-5}$
software_version | pycox 0.2.3
Random Survival Forests
node_size | 5 , 15, 50, 100
min_time_bins | 8, 16, 32
software_version | rfsrc 3.1.0
MOTOR
dropout | 0, 0.1
learning_rate |  $10^{-3}$, $10^{-4}$, $10^{-5}$
num_time_bins | 8, 16
survival_dim | 256, 512
inner_dim | 768
layers | 6
max_sequence_length | 16,384
vocabulary_size | 65,536
software_version | piton d9d22c6
Table S1: Hyperparmeter search grids of methods under comparison in our
experiments. The software version for implementing each method is also shown.

## Appendix B. Subsampling Bias

For our experiments, we use a simple subsampling strategy where we drop a
fraction $d$ of the censored examples in our dataset before modeling and
evaluation. This causes a non-uniform scaling of the hazard rates. Let
$H_{i}(t)$ be the hazard rate for a particular time t and let
$\widetilde{H_{i}}(t)$ be the hazard rate for a particular time $t$ induced by
our subsampling strategy. Equation 7 provides the ratio of
$\widetilde{H_{i}}(t)$ and $H_{i}(t)$. Proof Appendix B. Subsampling Bias is
the derivation of that equation.

$\tfrac{\widetilde{H_{i}}(t)}{H_{i}(t)}=\frac{1}{1-d\,\mathbb{P}(C_{i}<T_{i}|C_{i}>t,T_{i}>t)}$
(7)

###### Proof.

The hazard in our original dataset before sampling is

$\displaystyle H_{i}(t)$
$\displaystyle=\frac{\mathbb{P}(T_{i}=t)}{\mathbb{P}(T_{i}>t)}$
$\displaystyle=\frac{\mathbb{P}(T_{i}=t)\,\mathbb{P}(C_{i}>t)}{\mathbb{P}(T_{i}>t)\,\mathbb{P}(C_{i}>t)}$
$\displaystyle=\frac{\mathbb{P}(T_{i}=t,C_{i}>t)}{\mathbb{P}(T_{i}>t,C_{i}>t)}\hskip
28.45274pt\mathrm{assume}\,\,C_{i}\perp\perp T_{i}$
$\displaystyle=\frac{\mathbb{P}(T_{i}=t,C_{i}>t)}{\mathbb{P}(C_{i}<T_{i},C_{i}>t,T_{i}>t)+\mathbb{P}(C_{i}\geq
T_{i},C_{i}>t,T_{i}>t)}$

The hazard in our dataset after sampling is

$\displaystyle\widetilde{H_{i}}(t)$
$\displaystyle=\frac{\mathbb{P}(T_{i}=t,C_{i}>t)}{(1-d)\,\mathbb{P}(C_{i}<T_{i},T_{i}>t,C_{i}>t)+\mathbb{P}(C_{i}\geq
T_{i},T_{i}>t,C_{i}>t)}$
$\displaystyle=\frac{\mathbb{P}(T_{i}=t,C_{i}>t)}{((1-d)\,\mathbb{P}(C_{i}<T_{i}|T_{i}>t,C_{i}>t)+\mathbb{P}(C_{i}\geq
T_{i}|T_{i}>t,C_{i}>t))\,\mathbb{P}(T_{i}>t,C_{i}>t)}$
$\displaystyle=\frac{\mathbb{P}(T_{i}=t,C_{i}>t)}{((1-d)\,\mathbb{P}(C_{i}<T_{i}|T_{i}>t,C_{i}>t)+(1-\mathbb{P}(C_{i}<T_{i}|T_{i}>t,C_{i}>t)))\,\mathbb{P}(T_{i}>t,C_{i}>t)}$
$\displaystyle=\frac{\mathbb{P}(T_{i}=t,C_{i}>t)}{(1-d\,\mathbb{P}(C_{i}<T_{i}|T_{i}>t,C_{i}>t))\,\mathbb{P}(T_{i}>t,C_{i}>t)}$
$\displaystyle=\frac{H_{i}(t)}{1-d\,\mathbb{P}(C_{i}<T_{i}|T_{i}>t,C_{i}>t)}$

Thus,

$\frac{\widetilde{H_{i}}(t)}{H_{i}(t)}=\frac{1}{1-d\,\mathbb{P}(C_{i}<T_{i}|C_{i}>t,T_{i}>t)}$

∎

## Appendix C. Confidence Intervals

Table S2: 95% confidence intervals of difference in time-dependent C statistic
between MOTOR and the other methods on the full dataset. Higher numbers are
considered better for this metric. Asterisks indicate statistically
significant entries at p = 0.05. RSF = Random Survival Forest. MOTOR-WP =
MOTOR Without Pretraining.

Method | Celiac | Heart Attack | Lupus | NAFLD | Pancretic Cancer | Stroke
---|---|---|---|---|---|---
Cox PH | [-0.139, -0.086]* | [-0.143, -0.103]* | [-0.099, -0.060]* | [-0.142, -0.125]* | [-0.092, -0.049]* | [-0.104, -0.086]*
DeepSurv | [-0.122, -0.072]* | [-0.075, -0.045]* | [-0.078, -0.042]* | [-0.066, -0.053]* | [-0.075, -0.034]* | [-0.051, -0.037]*
RSF | [-0.093, -0.050]* | [-0.063, -0.033]* | [-0.081, -0.045]* | [-0.064, -0.051]* | [-0.060, -0.023]* | [-0.041, -0.028]*
MOTOR-WP | [-0.132, -0.079]* | [-0.110, -0.069]* | [-0.068, -0.026]* | [-0.044, -0.032]* | [-0.113, -0.063]* | [-0.050, -0.036]*
MOTOR | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000]

Table S3: 95% confidence intervals for the difference in ND calibration
between MOTOR and the methods on the full dataset. Lower numbers are better
for this metric. Asterisks indicate statistically significant entries at p =
0.05. RSF = Random Survival Forest. MOTOR-WP = Motor Without Pretraining.

Method | Celiac | Heart Attack | Lupus | NAFLD | Pancretic Cancer | Stroke
---|---|---|---|---|---|---
Cox PH | [0.011, 0.239]* | [0.068, 0.321]* | [0.173, 0.522]* | [0.072, 0.135]* | [-0.094, 0.192] | [-0.029, 0.030]
DeepSurv | [-0.081, 0.088] | [-0.028, 0.162] | [-0.035, 0.155] | [2.592, 3.267]* | [-0.120, 0.193] | [-0.037, 0.042]
RSF | [0.120, 0.426]* | [0.049, 0.323]* | [-0.050, 0.131] | [0.040, 0.108]* | [0.099, 0.490]* | [0.043, 0.125]*
MOTOR-WP | [2.556, 10.321]* | [2.073, 10.747]* | [10.942, 51.232]* | [0.023, 0.099]* | [2.859, 19.977]* | [-0.001, 0.156]
MOTOR | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000]

Table S4: 95% confidence intervals for the difference in MOTOR’s time-
dependent C statistic when the pretraining task is changed. We evaluate both a
“next code” baseline as well as our proposed time-to-event objective. Note
that these experiments are done with full fine-tuning of the entire model.
Asterisks indicate statistically significant entries at p = 0.05.

Objective | Celiac | Heart Attack | Lupus | NAFLD | Pancretic Cancer | Stroke
---|---|---|---|---|---|---
Next Code | [-0.043, -0.007]* | [-0.031, -0.007]* | [-0.029, -0.004]* | [-0.008, 0.000] | [-0.018, 0.014] | [-0.022, -0.012]*
Time-to-Event | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000]

Table S5: 95% confidence intervals for the difference in time-dependent C
statistic between pretraining on 100% of the training data vs more limited
subsets. Asterisks indicate statistically significant entries at p = 0.05.

% Used | Celiac | Heart Attack | Lupus | NAFLD | Pancretic Cancer | Stroke
---|---|---|---|---|---|---
5% | [-0.130, -0.084]* | [-0.052, -0.028]* | [-0.132, -0.089]* | [-0.099, -0.084]* | [-0.058, -0.022]* | [-0.057, -0.044]*
10% | [-0.109, -0.066]* | [-0.046, -0.023]* | [-0.124, -0.081]* | [-0.080, -0.067]* | [-0.052, -0.016]* | [-0.053, -0.040]*
25% | [-0.078, -0.041]* | [-0.041, -0.020]* | [-0.081, -0.044]* | [-0.066, -0.054]* | [-0.041, -0.014]* | [-0.047, -0.035]*
100% | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000]

Table S6: 95% confidence intervals for the difference in time-dependent C
statistic between fine-tuning on 100% of the training data vs subsets of the
data. Asterisks indicate statistically significant entries at p = 0.05.

% Used | Celiac | Heart Attack | Lupus | NAFLD | Pancretic Cancer | Stroke
---|---|---|---|---|---|---
5% | [-0.029, -0.005]* | [-0.010, -0.001]* | [-0.018, -0.003]* | [-0.010, -0.007]* | [-0.025, -0.004]* | [-0.009, -0.004]*
10% | [-0.033, -0.008]* | [-0.015, -0.004]* | [-0.012, 0.000] | [-0.006, -0.004]* | [-0.018, -0.003]* | [-0.008, -0.004]*
25% | [-0.018, -0.005]* | [-0.007, -0.000]* | [-0.010, 0.000] | [-0.005, -0.003]* | [-0.010, -0.002]* | [-0.006, -0.004]*
100% | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000] | [0.000, 0.000]