# Intelligent Vehicles and Transportation

Hüseyin Abut, Gerhard Schmidt, Kazuya Takeda, Jacob Lambert, and John H.L.
Hansen  Clément Courcelle  Northern Robotics Laboratory, Université Laval,
Québec, QC, Canada G1V 0A6  Dominic Baril  Northern Robotics Laboratory,
Université Laval, Québec, QC, Canada G1V 0A6  François Pomerleau and Johann
Laconte  Northern Robotics Laboratory, Université Laval, Québec, QC, Canada
G1V 0A6

###### Abstract

In the context of autonomous driving, vehicles are inherently bound to
encounter more extreme weather during which public safety must be ensured. As
climate is quickly changing, the frequency of heavy snowstorms is expected to
increase and become a major threat to safe navigation. While there is much
literature aiming to improve navigation resiliency to winter conditions, there
is a lack of standard metrics to quantify the loss of visibility of lidar
sensors related to precipitation. This chapter proposes a novel metric to
quantify the lidar visibility loss in real time, relying on the notion of
visibility from the meteorology research field. We evaluate this metric on the
Canadian Adverse Driving Conditions (CADC) dataset, correlate it with the
performance of a state-of-the-art lidar-based localization algorithm, and
evaluate the benefit of filtering point clouds before the localization
process. We show that the Iterative Closest Point (ICP) algorithm is
surprisingly robust against snowfalls, but abrupt events, such as snow gusts,
can greatly hinder its accuracy. We discuss such events and demonstrate the
need for better datasets focusing on these extreme events to quantify their
effect.

author#1 location#1 figurec

Clément Courcelle et al. On the Importance of Quantifying Visibility for
Autonomous Vehicles under Extreme Precipitation

On the Importance of Quantifying Visibility for Autonomous Vehicles under
Extreme PrecipitationOn the Importance of Quantifying Visibility for
Autonomous Vehicles under Extreme Precipitation

Keywords: Localization, Snow, Precipitation, Visibility, Lidar

##  1 Introduction

CBPL_Sec:Introduction

As robotics becomes more and more part of our everyday life, autonomous
vehicles start to roam the world alongside human drivers. In this regard,
localization is a crucial component of autonomous navigation. Lidar
lidarsensors have emerged as a standard for autonomous vehicle sensor suites,
providing dense information about the environment[1]. However, lidars can be
sensitive to the noise created by adverse weather like rain, fog or snow
because of the reflection of the laser beams by raindrops or snowflakes[2].
LABEL:CBPL_fig:intro shows an example of heavy snowfall that happened during
the winter in the campus of Quebec City, Canada. Such a level of precipitation
yields consequent noise on lidar data, where most of the lidar beams were
stopped by a snowflake, thus returning erroneous data. To this effect, various
snow removal filtersfilters have recently been proposed in the literature for
lidar scans[3, 4]. Datasets have also been released with the aim of developing
and evaluating localization algorithms robust to precipitation, such as the
Canadian Adverse Driving Conditions (CADC) dataset[5]. However, most work in
the literature uses qualitative metrics and there is no standard protocol to
quantify the snow precipitation in the vicinity of the vehicle.

Fig. 1: Robot driving in a snowstorm in Quebec City, Canada. Harsh snowstorms
yield considerable noise on lidar data, especially during snow gusts,
potentially leading to localization failure.CBPL_fig:intro

In this chapter, we propose a novel method to physically model the noise
encountered by lidars in snowstormssnowstorms. Using this method, we are able
to quantify the level of precipitation using only lidar information, relying
on the notion of visibility from the meteorology research field. We then
evaluate this metric on all sequences from the CADC dataset and correlate it
with the localization performance of a state-of-the-art lidar-based
localization algorithm. We demonstrate that using the data available in the
literature, lidar-based localization is currently robust to significant
snowfall. Nevertheless, we observe that localization algorithms will
potentially fail under extreme snowfall. This observation motivates additional
datasets that better feature extreme weather, effectively pushing localization
algorithms to their limit. Therefore, the contributions of this chapter are 1)
a novel metric to estimate the visibility in snowstorms; and 2) an in-depth
analysis of the impact of snowstorms on lidar-based localization.

##  2 Related Work

CBPL_Sec:related_work

Measurements of precipitationprecipitation type and intensity are often done
through the work of hydrologists and meteorologists through the use of a
disdrometer device[6]. These devices provide precipitation rates based on the
measures of the size and the velocity of precipitation particles. [7] measured
the rain intensity by using the rates provided by a meteorological station to
quantify the performance of lidar sensors under rain. For a similar objective
in rain and fog, [8] characterized weather by measuring meteorological
visibility, particle size distribution, and rainfall intensity. The
evaluations were performed indoor in fog chamber facilities, allowing the use
of several weather sensors. To provide outdoor field measurements, various
robotic deployments and datasets under snow precipitation have recently been
documented in the literature, measuring either the intensity or the density of
precipitation. [9] have deployed a vision-based teach-and-repeat navigation
system for 100 days at the University of Toronto Institute for Aerospace
Studies, conducting multi-seasonal, off-road navigation, for a total of
over$29\text{\,}\mathrm{km}$. Unfortunately, the intensity of snowfalls is
only qualitatively described in this work. Later, [10] deployed a similar
teach-and-repeat off-road navigation framework, but this time based on lidars.
The authors cumulated$18.8\text{\,}\mathrm{km}$ of autonomous navigation data
in boreal forest and classified the precipitation using the disdrometer
measures of a nearby meteorological station. [5] have released the CADC multi-
modal dataset including over$20\text{\,}\mathrm{km}$ of driving through two
years in winter under harsh weather in Waterloo, Canada. For this dataset, the
snow density is estimated by the number of points that are removed by the
Dynamic Radius Outlier Removal (DROR) filter[3]. This metric is dependent on
the lidar sensor used for the dataset, which is a Velodyne VLP-32C.
Alternatively, [4] have released the Winter Adverse Driving dataSet (WADS),
recorded in the Keweenaw peninsula in Michigan. In this dataset, the authors
quantify the intensity of snowfalls by measuring the snowfall rate in inches
per hour. Lastly,[11] have proposed the Boreas dataset, featuring
over$350\text{\,}\mathrm{km}$ of urban driving data in Canadian winter.
Precipitation characterization is limited to rain or snow in the Boreal
dataset, indicating neither its intensity nor its density. In this work, we
propose a novel metric to estimate the visibility in precipitation that is
adapted for real-time computation and mobile sensors measurements, while
correlating with the field of meteorology. Our goal is that this metric allows
standardizing how precipitation is quantified in various datasets while being
computable in real time for autonomous platforms.

Driving in adverse weatheradverse weather has revealed to be a major challenge
in autonomous navigation[12], leading to a lot of research aiming at analyzing
and reducing the impact of precipitation on different aspects of autonomous
navigation, such as localization and object detection. Moreover, [13] showed
that lidar-based mapping and localization algorithms perform surprisingly well
in snowstorms and even exceed radar-based algorithms. This result was also
observed in deployment by [10] who noted that it is not the snowfalls but
rather snow accumulated on the ground that disrupts the mapping and
localization. Various filters have also been proposed in the literature to
remove the precipitation noise on lidar scans. Firstly, the Radius Outlier
Removal (ROR) filter is included in the Point Cloud Library[14]. It is not
designed for snow specifically and therefore does not take into account that
the snow-induced noise is stronger close to the sensor, as shown by [15]. To
address this limitation, [3] developed the DROR filter. This filter is based
on the ROR, but the searching radius depends on the distance between the
current point and the sensor. This modification greatly improves the filter
efficiency for snow points further than five meters away from the sensor.
Secondly, the Statistical Outlier Removal (SOR) filter is also implemented in
the Point Cloud Library[14]. Similar to the ROR filter, it was not designed to
remove snow-induced noise specifically. Thus, [4] designed the Dynamic
Statistical Outlier Removal (DSOR) filter, using a threshold proportional to
the distance between the point and the sensor. The DSOR filter achieves better
recall but a slightly lower precision than the DROR, meaning that DSOR removes
more snow but also more environmental points. These works do not evaluate the
improvement of localization after the filter application, which was first done
by [16], who present the Dynamic Distance-Intensity Outlier Removal (DDIOR)
filter, relying on point intensity and neighborhood. The evaluation of the
localization accuracy was performed on one sequence of the WADS dataset. In
this work, we correlate the accuracy of a state-of-the-art localization
algorithm with snowstorm intensity. This snowstorm intensity is computed using
the novel visibility metric that we present in the next section.

##  3 Theory

CBPL_Sec:Theory

In this section, we look at the development of a density measure based on the
work of [17] who developed a method to measure density fieldsdensity fields in
the case of obstacle avoidance for safe navigation. We propose to extend their
method to our work, supplanting the obstacles field by the density field of
the snowstorm. This allows us to define the notion of visibilityvisibility in
the context of lidar measurements, with a strong connection with the notion of
visibility in the field of meteorologymeteorology [18, 19]. In the following,
we propose a method to estimate the snow density field around the robot. We
assume that the density of the snow can be modeled using a heterogeneous
Poisson point process, and that the density of the snow is constant in the $z$
direction for a small height difference. As such, we only take into account
points that fall within a $1\text{\,}\mathrm{m}$ strip in the $z$ direction,
centered on the position of the lidar. The probability
$\mathbb{P}\left(\mathtt{object}\right)$ that a beam goes through a space
$\mathcal{S}\subset\mathbb{R}^{2}$ and measures an object without colliding
with a snowflake can be computed as

$\displaystyle\mathbb{P}\left(\mathtt{object}\right)$
$\displaystyle=\exp\left(-\Lambda(\mathcal{S})\right),\quad\text{with
}\Lambda(\mathcal{S})=\int_{\mathcal{S}}\lambda(\bm{s})\operatorname{d\\!}\bm{s},$
(1)

where $\lambda(\bm{s})\in\mathbb{R}_{>0}$ is the density of snowflakes at the
position $\bm{s}\in\mathbb{R}^{2}$.

To estimate and store the density field, the space is tessellated into cells.
As such, each cell $c_{i}$ contains a density $\lambda_{i}$ corresponding to
the associated density at this position. Moreover, the lidar sensor is modeled
as having a collision area $A_{c}$, meaning that for each measurement, the
snowflake that caused the beam to be reflected was in a region of area $A_{c}$
around the measurement. As a recall, the density field can be estimated with

$\lambda_{i}=\frac{1}{A_{c}}\ln\left(1+\frac{h_{i}}{m_{i}}\right),{eq\mathrel{\mathop{\mathchar
58\relax}}lambda}$ (2)

where $h_{i}$ is the number of times a collision occurred in the cell $c_{i}$
and $m_{i}$ the number of times a beam went through the cell $c_{i}$ without
collision[17]. For extensive details on the mathematical proofs, please refer
to[20].

Finally, one can note that the estimated field also changes in time. As this
additional parameter greatly complexifies the estimation of the field, we
assume that the field is constant for a given time interval
$[t-\tau/2,t+\tau/2]$ and use all measurements in this interval to estimate
the field at time $t$. We used an interval size of
$\tau=$1\text{\,}\mathrm{s}$$ in our experiments, and each measurement has an
area of $A_{c}=$40\times 40\text{\,}{\mathrm{cm}}^{2}$$. Additionally, the
field is tessellated into cells of size $10\times
10\text{\,}{\mathrm{cm}}^{2}$. LABEL:CBPL_fig:density_field_example shows an
example of the resulting density field for different snowstorms. The left
field corresponds to a run of the CADC when the snowfall was calm and without
wind. As a result, the density around the robot is equally distributed and low
enough that one can distinguish structural elements in the point cloud. In the
second case, the field is estimated during a heavy gust of wind and snow,
yielding to a higher density of snow around the robot. This event was recorded
in the campus of Quebec City, in the same snowstorm represented on
LABEL:CBPL_fig:intro. Note that the major issue is that because of this higher
density, the range of the lidar is greatly reduced, and almost no structured
elements are visible in this case.

Fig. 2:  Examples of snow density fields with their associated point clouds.
Point clouds are colored according to the point's height for easier reading.
In the zoomed rectangles, gray areas correspond to unobserved space by the
lidar beams. Left: The robot is navigating while a weak snowstorm is in
effect, resulting in a density field of low values. Most of the environment is
still visible, where walls and other landmarks are easily recognizable. Right:
The robot is in a heavy snowstorm with a snow gust in front of it. Most of the
points correspond to snowflakes, and no clear object can be identified in the
point cloud. CBPL_fig:density_field_example

Then, we demonstrate the close relationship between our method and the concept
of visibility that was introduced in the textbook from [18] in Middleton1957,
then adapt this concept in the case of lidar sensors. In the field of
meteorology, visibility is defined as the distance at which a black object
becomes invisible. More precisely, an object is visible if the apparent
contrast between it and the background is high enough for a human observer to
distinguish it. The contrast $C$ is defined as

$C=\exp\left(-\sigma r\right),{CBPL_{e}q\mathrel{\mathop{\mathchar
58\relax}}apparent_{c}ontrast}$ (3)

where $r$ is the distance between the observer and the object, and $\sigma$ is
the atmospheric extinction coefficient average over the distance $r$ and over
the visible spectrum. The visibility $V$ is then defined as the distance at
which the contrast $C$ drops below a threshold $\epsilon$. In the following,
we propose to mimic this metric in the case of lidar measurements instead of a
human observer. We define the notion of p-visibility, which is the distance at
which the probability of perceiving an object with a lidar is at most $p$. For
a given field $\lambda(\cdot)$, the probability that the lidar will be able to
measure an object without being interrupted by a snowflake is

$\begin{aligned}
\mathbb{P}\left(\mathtt{object}\right)&=\exp\left(-\Lambda(\mathcal{S})\right)\\\
&=\exp\left(-\overline{\lambda}A_{\mathcal{S}}\right),\quad\overline{\lambda}=\frac{\Lambda(\mathcal{S})}{A_{\mathcal{S}}},\end{aligned}{CBPL_{e}q\mathrel{\mathop{\mathchar
58\relax}}p_{m}easure}$ (4)

where $\mathcal{S}$ is the space that the beam needs to cross to reach the
object to measure, $A_{\mathcal{S}}$ its corresponding area and
$\overline{\lambda}$ the average density of the field. Note that this equation
is strictly equivalent to LABEL:CBPL_eq:apparent_contrast with
$\sigma=\overline{\lambda}$ and $r=A_{\mathcal{S}}$. Assuming the area the
beam is traveling in a sector shape, its area is equal to
$A_{\mathcal{S}}=\frac{1}{2}d^{2}\alpha,$ where $d$ is the traveled distance
of the beam, and $\alpha$ is the aperture angle of the lidar sensor [21].
Using LABEL:CBPL_eq:p_measure, we can find the distance $V_{p}$ at which the
probability to perceive the object drops below the probability $p$:

$\begin{aligned}
&\mathbb{P}\left(\mathtt{object}\right)=\exp\left(-\overline{\lambda}\cdot\frac{V_{p}^{2}\alpha}{2}\right)=p\\\
\Leftrightarrow\quad&V_{p}=\sqrt{\frac{-2\ln{p}}{\overline{\lambda}\alpha}}.\end{aligned}{CBPL_{e}q\mathrel{\mathop{\mathchar
58\relax}}p_{v}isibility}$ (5)

Using this metric, one can estimate the p-visibility at any time in a
snowstorm. Furthermore, note that the whole pipeline can easily run online,
and thus the p-visibility can be estimated in real time.

One can note strong connections between the contrast visibility in
LABEL:CBPL_eq:apparent_contrast and ours from LABEL:CBPL_eq:p_visibility.
Indeed, the atmospheric extinction $\sigma$ corresponds in our case to the
local density of the snowfall $\overline{\lambda}$, that both serve the same
purpose. The higher these coefficients, the less visibility the current
observer has. In the first case, the coefficient depicts the reflectivity of
the atmosphere, whereas in our case the reflectivity is replaced by the
reflection of the snowflakes. As such, the density can also be seen as a
measure of reflectivity. The final difference between the equations comes from
the fact that a lidar beam travel in a cone shape [21], contrary to a single
light ray, leading to more and more space being occupied by the beam as it
travels in space. Using LABEL:CBPL_eq:p_visibility, we are able to compute the
p-visibility for any given field in real time. The average density of the
field $\overline{\lambda}$ is estimated using data that are at most
$5\text{\,}\mathrm{m}$ away from the robot. In the following, we fix the
probability at $p=0.5$ and shorten the term 0.5-visibility as simply
visibility. Thus, visibility is the distance at which a lidar beam has half a
chance of reaching an object and not being deflected by a snowflake.

As an example, the two density fields depicted in
LABEL:CBPL_fig:density_field_example lead to very different visibilities.
Using a standard lidar sensor, such as the Velodyne HDL-32E or Robosense RS-32
with an aperture angle of $0.085\text{\,}\mathrm{\SIUnitSymbolDegree}$, we can
compute the visibility of the lidar for a calm snowfall and a snow gust during
a heavier storm. For the snowfall (LABEL:CBPL_fig:density_field_example –
Left), the visibility is equal to approximately $40\text{\,}\mathrm{m}$. In
the case of the snow gust (LABEL:CBPL_fig:density_field_example – Right), the
visibility is reduced to only $5\text{\,}\mathrm{m}$, meaning that half of the
time the lidar will not be able to measure an obstacle past
$5\text{\,}\mathrm{m}$, and will instead return noise.

##  4 Results

CBPL_Sec:Results In this section, we evaluate our visibility metric on the
CADC dataset and correlate it with the performance of a state-of-the-art
localizationlocalization algorithm. A thorough description of the sensors used
to record the CADC dataset can be found in[5]. Firstly, we describe the
algorithm, as well as the localization performance metric that is used
throughout this section. Secondly, we present an analysis of the impact of
snowfall on localization. Then, we discuss the improvements of the
localization accuracy that come from filtering the snow from the point cloud.

To evaluate the correlation between our visibility metric and localization
accuracy, we have used a lidar-based framework similar to the one described by
[10]. The framework is based on the Iterative Closest Point Iterative Closest
Point (ICP) registration algorithm. This algorithm computes the rigid
transformation that minimizes the registration error between two subsequent
scans, with the help of a prior based on Inertial Measurement Unit (IMU) and
wheel encoder measurements. As demonstrated by [10], real-time systems cannot
compute localization and mapping by taking all lidar points into account due
to the high computational complexity. As such, a random sub-sampling keeping
$70\text{\,}\mathrm{\char 37\relax}$ of the input point cloud has been used.
Additionally, it should be noted that we use a map filter to remove dynamic
points after scan registration[22], helping to preserve the consistency of the
map. To evaluate the quality of the localization, we compute the Relative Pose
Error (RPE) as described by [23]. The estimated trajectory is compared to
Global Navigation Satellite System (GNSS) Real Time Kinematic (RTK) data,
considered as ground truth. We set the window duration to
$1\text{\,}\mathrm{s}$, as we seek to pinpoint sharp, abrupt localization
errors that would lead to dangerous behaviors of the vehicle. The evaluations
have been performed using the evo evaluation
library.111github.com/MichaelGrupp/evo

First, we analyze the impact of snowfall on localization.
LABEL:CBPL_fig:rpe_density shows the relative pose error as a function of the
visibility for all the sequences of the CADC. As we randomly subsample the
input point cloud, the dataset have been fully evaluated 50 times, leading to
$46\text{\,}\mathrm{h}$ of data. One can see that there is no clear
correlation between visibility and localization error. Surprisingly, the ICP
algorithm is resilient to a great level of snow precipitation, even when half
of the lidar beams cannot see past $8\text{\,}\mathrm{m}$. However, for low
visibility, some outliers tends to lengthen the error distribution up to
$40\text{\,}\mathrm{\char 37\relax}$.

Fig. 3:  Relative Pose Error (RPE) correlated with the visibility, estimated
as described in LABEL:CBPL_Sec:Theory. A logarithmic scale is used on the left
axis, and data are gathered into bins of $2.2\text{\,}\mathrm{m}$. The number
of points constituting each bin is displayed on the right axis. As the dataset
is focused on harsh weather conditions, most measurements happen in low
visibility. No clear correlation between the visibility and the localization
error can be established. CBPL_fig:rpe_density

To better highlight these outliers, LABEL:CBPL_fig:snowgust_time depicts the
localization error for two given runs of the CADC. Run 70 contains heavy snow
precipitation with an average visibility of ten meters, but maintains a RPE
below $3\text{\,}\mathrm{\char 37\relax}$. The environment features several
buildings, allowing the ICP to rely on many planes to provide an accurate
localization. Run 78 is recorded in a more rural area, meaning that there are
very few landmarks on the side of the road and particularly no buildings. The
environment is therefore much less constrained, resulting in a higher
localization error. In this specific run, the error is relatively constant,
but for two distinct events. The first one corresponds to a heavy gust caused
by a passing truck lifting powder snow in front of the lidar sensor. During
such an event, a local, very dense cloud of snowflakes can greatly obstruct
the field of view, as depicted in LABEL:CBPL_fig:density_field_example –
Right. As the density of snow reaches an extreme level, a large part of the
lidar scan consists of beams that hit snowflakes, instead of the static
surrounding environment. The second event simply corresponds to dynamic
obstacles (i.e., cars) obstructing the field of view of the lidar, leading to
a decreased localization accuracy[22]. One can observe that in the case of
wind gusts, the third quartile of the RPE error for this experiment's
distribution is very high, going above $10\text{\,}\mathrm{\char 37\relax}$.
Such an error can easily disrupt the entire navigation stack, leading to
hazardous behaviors in case of autonomous driving. This observation leads to
theorize that ICP-based localization is sensible to which points are randomly
filtered in the scan, showing potential localization failures if the system is
subject to denser, heterogeneous snowfall. However, capturing such events is
challenging as they are quite rare, and thus more data is needed before being
able to draw quantitative evaluations about these events.

Fig. 4: Comparison of the Relative Pose Error (RPE) of three sequences. Run 70
presents heavy falling snow in an urban environment. Run 78 is recorded on a
more rural road than most of the other sequences, and includes events like a
snow gust that disrupts the localization. The DSOR filter was applied to Run
78 to measure the improvement on localization accuracy.CBPL_fig:snowgust_time

We also look at the localization performance gain related to filtering the
snow with the DSOR filter[4]. LABEL:CBPL_fig:snowgust_time shows that an
advantage of using a filter is that it greatly reduces the interquartile range
of the error during events like a snow gust, meaning that the chance of a
total failure of the localization process is also decreased. Once again, a
generalization of this result would require more snow gust data, as these
events are scarce and only one snow gust has been recorded in the CADC.

LABEL:CBPL_fig:filters depicts the localization error given different
parameters of the DSOR filter, which is a state-of-the-art snow removal filter
at the time of writing. The filter has two parameters $s,r\in\mathbb{R}_{>0}$.
The parameter $s$ is used in the SOR filter[14] to determine the proportion of
points that are classified as outliers, while $r$ is added for DSOR to scale
the filter threshold with the distance between the point and the sensor. The
lower the parameters $s,r$ are, the more aggressive the filter becomes. Note
that setting $s,r\rightarrow+\infty$ is equivalent to not applying any filter
on the point cloud. Overall, we can see that applying the filter at its best
configuration yields an improvement of about $1\text{\,}\mathrm{\char
37\relax}$ on the median of the RPE. This aggressive tuning explains that the
noise of the snow gust can be filtered, even though the points are denser than
in usual snow precipitation.

Fig. 5: DSOR filter applied to Run 78 containing the snow gust event. Multiple
couples of parameters ${r,s}$ are tested to find the best configuration,
indicated by the red dot, which reduces the median of the Relative Pose Error
(RPE) in the sequence.CBPL_fig:filters

We have evaluated our visibility metric, described inLABEL:CBPL_Sec:Theory on
all runs of the CADC dataset and shown its correlation with localization
performance inLABEL:CBPL_fig:rpe_density. Through this result, we observe that
extreme precipitation conditions can potentially lead to localization failure,
however the data is too sparse to confirm. We argue that our metric could help
standardize the evaluation of the loss of visibility related to precipitation
on all datasets. Our results also motivate the need for more datasets that
feature extreme precipitations and visibility loss to challenge state-of-the-
art localization algorithms.

##  5 Conclusion

CBPL_Sec:Conclusion In this chapter, we introduced a novel visibility metric
for lidars in snowstorm scenarios, relying on the notion of visibility in the
meteorological field. We evaluated the CADC dataset with this metric, and
showed that there is no clear correlation between the visibility in snowstorms
and the localization performance of the ICP algorithm. However, events such as
snow gusts can lead to localization failure, but available datasets contain
too few recordings of these events to quantify their true impact. As future
work, we aim to use our metric to develop precipitation-aware filters and
vehicle control. The visibility metric can also be used to characterize and
classify snowstorm datasets from the robotics community. Furthermore, research
on new datasets targeting specifically extreme and short events, such as snow
gusts, would greatly improve the resilience of localization algorithms.

## Acknowledgements

This research was supported by the Natural Sciences and Engineering Research
Council of Canada (NSERC) through the grant CRDPJ 527642-18 SNOW (Self-driving
Navigation Optimized for Winter).

## References

  * [1] Ricardo Roriz, Jorge Cabral and Tiago Gomes  ``Automotive LiDAR Technology: A Survey''  In _IEEE Transactions on Intelligent Transportation Systems_ 23.7 IEEE, 2021, pp. 6282–6297  DOI: 10.1109/TITS.2021.3086804
  * [2] R. H. Rasshofer, M. Spies and H. Spies  ``Influences of weather phenomena on automotive laser radar systems''  In _Advances in Radio Science_ 9.B.2, 2011, pp. 49–60  DOI: 10.5194/ars-9-49-2011
  * [3] Nicholas Charron, Stephen Phillips and Steven L. Waslander  ``De-noising of lidar point clouds corrupted by snowfall''  In _Proceedings - 15th Conference on Computer and Robot Vision_ Toronto, Canada: IEEE, 2018, pp. 254–261  DOI: 10.1109/CRV.2018.00043
  * [4] Akhil Kurup and Jeremy Bos  ``DSOR: A Scalable Statistical Filter for Removing Falling Snow from LiDAR Point Clouds in Severe Winter Weather''  In _arXiv preprint arXiv:2109.07078_ , 2021  arXiv: http://arxiv.org/abs/2109.07078
  * [5] Matthew Pitropov et al.  ``Canadian Adverse Driving Conditions dataset''  In _International Journal of Robotics Research_ 40.4-5, 2021, pp. 681–690  DOI: 10.1177/0278364920979368
  * [6] Ali Tokay, David B. Wolff and Walter A. Petersen  ``Evaluation of the New Version of the Laser-Optical Disdrometer, OTT Parsivel2''  In _Journal of Atmospheric and Oceanic Technology_ 31.6, 2014, pp. 1276–1288  DOI: 10.1175/JTECH-D-13-00174.1
  * [7] A. Filgueira et al.  ``Quantifying the influence of rain in LiDAR performance''  In _Measurement: Journal of the International Measurement Confederation_ 95 Elsevier Ltd, 2017, pp. 143–148  DOI: 10.1016/j.measurement.2016.10.009
  * [8] M. Kutila et al.  ``Automotive LiDAR performance verification in fog and rain''  In _21st Conference on Intelligent Transportation Systems_ Maui, USA: IEEE, 2018, pp. 1695–1701  DOI: 10.1109/ITSC.2018.8569624
  * [9] Kirk MacTavish, Michael Paton and Timothy D. Barfoot  ``Selective memory: Recalling relevant experience for long-term visual localization''  In _Journal of Field Robotics_ 35.8, 2018, pp. 1265–1292  DOI: 10.1002/rob.21838
  * [10] Dominic Baril et al.  ``Kilometer-scale autonomous navigation in subarctic forests: challenges and lessons learned''  In _Field Robotics_ 2.50, 2022, pp. 1628–1660  arXiv: http://arxiv.org/abs/2111.13981
  * [11] Keenan Burnett et al.  ``Boreas: A Multi-Season Autonomous Driving Dataset''  In _arXiv preprint arXiv:2203.10168_ , 2022  arXiv: http://arxiv.org/abs/2203.10168
  * [12] Jessica Van Brummelen, Marie O'Brien, Dominique Gruyer and Homayoun Najjaran  ``Autonomous vehicle perception: The technology of today and tomorrow''  In _Transportation Research Part C: Emerging Technologies_ 89 Elsevier, 2018, pp. 384–406  DOI: 10.1016/j.trc.2018.02.012
  * [13] Keenan Burnett et al.  ``Are We Ready for Radar to Replace Lidar in All-Weather Mapping and Localization?''  In _IEEE Robotics and Automation Letters_ 7.4, 2022, pp. 10328 –10335  arXiv: http://arxiv.org/abs/2203.10174
  * [14] Radu Bogdan Rusu and Steve Cousins  ``3D is here: Point Cloud Library (PCL)''  In _International Conference on Robotics and Automation_ Shangai, China: IEEE, 2011  DOI: 10.1109/ICRA.2011.5980567
  * [15] Sébastien Michaud, Jean-François Lalonde and Philippe Giguère  ``Towards Characterizing the Behavior of LiDARs in Snowy Conditions''  In _IROS_ Hamburg, Germany: IEEE, 2015, pp. 1–6
  * [16] Weiqi Wang et al.  ``A Scalable and Accurate De-Snowing Algorithm for LiDAR Point Clouds in Winter''  In _Remote Sensing_ 14.6, 2022, pp. 1468  DOI: 10.3390/rs14061468
  * [17] J. Laconte et al.  ``Lambda-Field: A continuous counterpart of the Bayesian occupancy grid for risk assessment''  In _IROS_ Macau, China: IEEE, 2019, pp. 167–172
  * [18] William Edgar Knowles Middleton  ``Vision through the Atmosphere''  In _University of Toronto Press_ , 1952
  * [19] Roy M. Rasmussen et al.  ``The estimation of snowfall rate using visibility''  In _Journal of Applied Meteorology_ 38.10, 1999, pp. 1542–1563  DOI: 10.1175/1520-0450(1999)038<1542:TEOSRU>2.0.CO;2
  * [20] Johann Laconte et al.  ``A novel occupancy mapping framework for risk-aware path planning in unstructured environments''  In _Sensors_ 21.22, 2021, pp. 1–30  DOI: 10.3390/s21227562
  * [21] Johann Laconte, Simon-Pierre Deschenes, Mathieu Labussiere and Francois Pomerleau  ``Lidar Measurement Bias Estimation via Return Waveform Modelling in a Context of 3D Mapping''  In _ICRA_ Montreal, Canada: IEEE, 2019, pp. 8100–8106  DOI: 10.1109/ICRA.2019.8793671
  * [22] Francois Pomerleau et al.  ``Long-term 3D map maintenance in dynamic environments''  In _ICRA_ Hong Kong: IEEE, 2014, pp. 3712–3719  DOI: 10.1109/ICRA.2014.6907397
  * [23] Jrgen Sturm et al.  ``A benchmark for the evaluation of RGB-D SLAM systems''  In _IROS_ Vilamoura, Portugal: IEEE, 2012, pp. 573–580  DOI: 10.1109/IROS.2012.6385773

  *[CADC]: Canadian Adverse Driving Conditions
  *[ICP]: Iterative Closest Point
  *[DROR]: Dynamic Radius Outlier Removal
  *[WADS]: Winter Adverse Driving dataSet
  *[ROR]: Radius Outlier Removal
  *[SOR]: Statistical Outlier Removal
  *[DSOR]: Dynamic Statistical Outlier Removal
  *[DDIOR]: Dynamic Distance-Intensity Outlier Removal
  *[IMU]: Inertial Measurement Unit
  *[RPE]: Relative Pose Error
  *[GNSS]: Global Navigation Satellite System
  *[RTK]: Real Time Kinematic


# Comparison of Baseline Covariate Adjustment Methods for Restricted Mean
Survival Time

Keisuke Hanada  Kyowa Kirin Co., Ltd  Junji Moriya  Kyowa Kirin Co., Ltd
Masahiro Kojima111Address: Biometrics Department, R&D Division, Kyowa Kirin
Co., Ltd. Otemachi Financial City Grand Cube, 1-9-2 Otemachi, Chiyoda-ku,
Tokyo, 100-004, Japan. Tel: +81-3-5205-7200 E-Mail:
<EMAIL_ADDRESS>Kyowa Kirin Co., Ltd  The Institute of
Statistical Mathematics

###### Abstract

The restricted mean survival time is a clinically easy-to-interpret measure
that does not require any assumption of proportional hazards. We focus on two
ways to directly model the survival time and adjust the covariates. One is to
calculate the pseudo-survival time for each subject using leave-one-out, and
then perform a model analysis using all pseudo-values to adjust for
covariates. The pseudo-survival time is used to reflect information of
censored subjects in the model analysis. The other method adjusts for
covariates using subjects for whom the time-to-event was observed while
adjusting for the censored subjects using the inverse probability of censoring
weighting (IPCW). This paper evaluates the performance of these two methods in
terms of the power to detect group differences through a simple example
dataset and computer simulations. The simple example illustrates the intuitive
behavior of the two methods. With the method using pseudo-survival times, it
is difficult to interpret the pseudo-values. We confirm that the pseudo-
survival times are different from the actual data obtained in a primary
biliary cholangitis clinical trial because of the many censored data. In the
simulations, the method using IPCW is found to be more powerful. Even in the
case of group differences with respect to the censor incidence rates and
covariates, the method using IPCW maintains a nominal significance level for
the type-1 error rate. We conclude that the IPCW method should be used to
estimate the restricted mean survival time when adjusting the covariates.

Key words and phrases: restricted mean survival time, inverse probability of
censoring weighting

##  1 Introduction

The proportional hazards (PH) Cox model is a well-known method for adjusting
the baseline covariates when using time-to-event data. The PH Cox model
evaluates the influence of baseline covariate parameters through hazard
ratios. However, a statistically valid assessment of the PH Cox model requires
an assumption of PH. In fact, in anticancer drugs such as immune checkpoint
inhibitors, the survival curves of the treatment and control groups coincide
for some time after the start of treatment, and then separate [1, 2]. For such
survival curves, the PH assumption may not hold. Spruance et al. [3] reported
that physicians have difficulty in helping patients to understand the
treatment effect when assessed in terms of hazard ratios. Therefore, in recent
years, the restricted mean survival time (RMST), which does not require the PH
assumption and uses familiar averages, has attracted considerable attention.
Because the RMST uses the mean survival time as a summary measure, differences
between treatment groups are easy to interpret. The mean survival time is
calculated using the area under the survival percentage curve, as estimated by
the Kaplan-Meier method, which is adjusted for censored subjects. In other
words, the mean survival time is calculated as
${\widehat{\mu}}=\int^{\infty}_{0}{\widehat{S}}(t)dt$, where
${\widehat{S}}(t)$ is the Kaplan-Meier estimator at time $t$. The Kaplan-Meier
estimator, however, cannot be defined when the largest observation is
censored. To avoid this, it is necessary to restrict time to a certain time
point $\tau$.

A number of baseline covariate adjustment methods for RMST have been proposed.
The adjustment methods fall into indirect and direct categories. The indirect
method uses a PH Cox model with estimated cumulative hazards to derive the
RMST [4, 5]. Although the advantage of RMST is that it can be calculated
without assuming PH, the indirect adjustment method requires the assumption of
PH, which may not be appropriate. Additionally, it requires the estimation of
cumulative hazards. To overcome these problems, direct methods, which model
survival and adjust the covariates directly without the assumption of PH, have
been developed by Andersen et al. [6] and Tian et al. [7]. An RMST that is
adjusted for time-dependent covariates after initial treatment has been
proposed by Zhang et al. [8].

In the course of a controlled clinical trial with a primary endpoint that can
be characterized as a time-to-event, we had planned to analyze the difference
in RMSTs between two groups by adjusting the baseline covariates. However, we
were faced with the question of whether Andersen’s or Tian’s method should be
applied. Ambrogi et al. [9] examined the accuracy of parameter estimation for
Andersen’s and Tian’s methods, but the difference in statistical power between
these methods has not yet been confirmed. Karrison and Kocherginsky [10]
examined the effects of adjustment for covariates, but did not identify the
differences among adjustment methods.

In this research, we compare the performance of Andersen’s and Tian’s methods
in terms of the power to detect group differences through a simple example
dataset and simulations. We illustrate the intuitive behavior of the two
methods in the simple example, and evaluate the detection power through
computer simulations with Cox’s model as a reference. Simulations are
performed under various scenarios. Because Andersen’s and Tian’s methods can
be easily performed with the RMSTREG PROCEDURE in SAS, we use the SAS
environment to conduct the simulations. The SAS program code used in the
simulations is summarized in the Supplemental Material. Additionally, we
confirm the difference in RMSTs adjusted by covariates using data from two
actual trials.

The remainder of this paper is organized as follows. Section 2 introduces
Andersen’s and Tian’s methods. Section 3 describes the setting and results of
the computer simulations, before Section 4 describes the results of the two
actual clinical trials. We conclude this paper with a discussion of our
results in Section 5. The program code used for the SAS simulations, as well
as the SAS analysis program code for the two actual clinical trials, is
included in the Supplemental Material.

##  2 Methods

This section first presents a definition of RMST. Next, we describe the two
methods of adjustment for the baseline covariates, one using the pseudo-
survival times (PSTs) and the other based on inverse probability censoring
weighting (IPCW).

###  2.1 Definition of RMST

Suppose that $T$ is the survival time subject to right censoring by a
censoring time $C$, which is independent of $T$. The observable survival time
is $R=\min(T,C)$ with the censoring indicator $I(T\leq C)$. For a sample size
of $n$, the data of the $i$-th subject is $R_{i}$ for $i=1,2,\ldots,n$. For
any restricted time $\tau>0$, the RMST is defined as

$\displaystyle\mu(\tau)$ $\displaystyle=E[\min(T,\tau)]=\int_{0}^{\tau}S(t)dt$
(1)

where $S(t)=\Pr(T>t)$ is the survival function. The estimator of RMST
$\mu(\tau)$ is $\hat{\mu}(\tau)=E[\min(R,\tau)]=\int_{0}^{\tau}\hat{S}(t)dt$.
In a two-arm clinical trial with survival functions $S_{0}(t)$ and $S_{1}(t)$
in the control and treatment groups, respectively, the difference $d$ in RMST
between the groups is given by

$\displaystyle d$
$\displaystyle=\int_{0}^{\tau}S_{1}(t)dt-\int_{0}^{\tau}S_{0}(t)dt=\int_{0}^{\tau}[S_{1}(t)-S_{0}(t)]dt$
(2)

i.e., $d$ is the area between the survival curves. In the standard case, we
can test the null hypothesis $d=0$ by comparing
$z=\hat{d}/\mathrm{SE}(\hat{d})$ with a student-$t$ or (in large samples)
normal reference distribution [11].

For the covariate-adjusted survival function, a method using Cox’s model has
been proposed [4, 5], but this is computationally complex because the
nonparametric hazard function must be calculated. Therefore, to resolve the
problems of the indirect method, direct methods for modeling survival and
adjusting covariates without the assumption of PH were developed by Andersen
et al. [6] and Tian et al. [7].

###  2.2 Model-based estimation using PSTs

Andersen et al. [6] proposed a covariate adjustment method using PSTs to
consider the censored subjects. Let $X_{i}=$min$(R_{i},\tau)$, $i=1,\dots,n$,
be random variables. We assume that the RMST is $\mu(\tau)=E[X_{i}]$ and the
estimator of $\mu(\tau)$ is ${\widehat{\mu}}(\tau)$. The PST of the $i$-th
subject is defined as

$\displaystyle{\widehat{\mu}}_{i}$
$\displaystyle=n{\widehat{\mu}}(\tau)-(n-1){\widehat{\mu}}^{-i}(\tau)$ (3)

where ${\widehat{\mu}}(\tau)$ is calculated as the RMST using all subjects and
${\widehat{\mu}}^{-i}(\tau)$ is the “leave-one-out” estimator for $\mu(\tau)$
based on $X_{j},j\neq i$.

With ${\widehat{\mu}}_{i}$ as the objective variable, regression analysis can
be performed with a generalized linear model for the $q$ covariates
${\text{\boldmath$Z$}}_{i}$. The generalized linear model is defined as
follows:

$\displaystyle g({\widehat{\mu}}_{i})$
$\displaystyle={\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i}$ (4)

with the link function $g(\cdot)$, where $\beta$ is a $q$-dimensional
parameter vector. For the survival time, the link function is set to the
identity function or the log function. The regression parameters $\beta$ for
PST ${\widehat{\mu}}_{i}$ can be estimated using generalized estimating
equations as

$\displaystyle{\text{\boldmath$U$}}_{PST}({\text{\boldmath$\beta$}})$
$\displaystyle=\sum_{i=1}^{n}{\text{\boldmath$U$}}_{i}({\text{\boldmath$\beta$}})$
$\displaystyle=\sum_{i=1}^{n}\left(\\{{\widehat{\mu}}_{i}-g^{-1}({\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i})\\}\frac{\partial}{\partial{\text{\boldmath$\beta$}}}g^{-1}({\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i})\right)$
$\displaystyle=0.$ (5)

We state the generalized estimating equation in the appendix for the cases in
which the identity and log functions are applied as the link function.

Andersen et al. [12] showed that consistent estimates of $\beta$ could be
obtained from (2.2) and that variance estimates for the solution
${\widehat{\text{\boldmath$\beta$}}}$ could be obtained from the standard
sandwich estimator

$\displaystyle\hat{V}$
$\displaystyle=I(\hat{{\text{\boldmath$\beta$}}})^{-1}\hat{\mbox{var}}\\{U_{PV}({\text{\boldmath$\beta$}})\\}I(\hat{{\text{\boldmath$\beta$}}})^{-1},$
(6)

with

$\displaystyle
I({\text{\boldmath$\beta$}})=\frac{1}{{\hat{v}}^{2}}\sum_{i=1}^{n}\left(\frac{\partial
g^{-1}({\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i})}{\partial{\text{\boldmath$\beta$}}}\right)^{T}\left(\frac{\partial
g^{-1}({\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i})}{\partial{\text{\boldmath$\beta$}}}\right)$
(7)
$\displaystyle\hat{\mbox{var}}\\{U_{PV}({\text{\boldmath$\beta$}})\\}=\frac{1}{{\hat{v}}^{4}}\sum_{i=1}^{n}U_{i}(\hat{{\text{\boldmath$\beta$}}})U_{i}(\hat{{\text{\boldmath$\beta$}}})^{T},$
(8)

where ${\hat{v}}^{2}$ is the variance of ${\widehat{\mu}}_{i}$.

###  2.3 Model-based estimation using IPCW

Tian et al. [7] proposed a covariate adjustment method weighted by the
survival function of censoring. The influence of the censored data can be
adjusted by weighting the inverse of the survival function based on IPCW. The
IPCW estimating function for $\beta$ in model (4) is

$\displaystyle U_{IPCW}({\text{\boldmath$\beta$}})$
$\displaystyle=\sum_{i=1}^{n}\frac{I(X_{i}\leq
C_{i})\\{X_{i}-g^{-1}({\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i})\\}}{\hat{S}(X_{i})}{\text{\boldmath$Z$}}_{i}=0,$
(9)

where $\hat{S}(\cdot)$ is the Kaplan-Meier estimator of the censoring time $C$
based on $R_{i}$ and the censoring indicator $I(X_{i}\leq C_{i})$ for
$i=1,\dots,n$. We state the generalized estimating equation in the appendix
for the cases in which the identity and log functions are applied as the link
function.

Let $\hat{{\text{\boldmath$\beta$}}}$ be the unique root of
$U_{IPCW}({\text{\boldmath$\beta$}})=0$. An important property for building a
prediction model is that $\hat{{\text{\boldmath$\beta$}}}$ converges to a
constant $\bar{{\text{\boldmath$\beta$}}}$ in probability, even when model (4)
is misspecified. Additionally,
$n^{1/2}(\hat{{\text{\boldmath$\beta$}}}-\bar{{\text{\boldmath$\beta$}}})$
converges weakly to a normal distribution with mean zero as $n\to\infty$ [7].

###  2.4 Examples of Andersen’s and Tian’s methods

We now illustrate the intuitive behavior of Andersen’s and Tian’s methods with
a simple example dataset. The example dataset is presented in Table 1. This
dataset assumes that 12 subjects are followed and their survival times
observed, and consists of the survival time (weeks), censoring (yes/no), and
age (covariate). The cut-off is assumed to be 100 weeks. We add the PSTs
derived by leave-one-out in Andersen’s method. The PSTs are almost identical
to the original survival times when the event occurs. However, when the data
are censored, the time at which the event is supposed to have occurred is
calculated. The link function is set to identity to understand the
interpretation of the analysis results easily.

Table 1: Example dataset

Treatment group | Control group
---|---
ID | ST | Censor | Age | PST | ID | ST | Censor | Age | PST
1 | 20 | Yes | 60 | 78.4 | 7 | 20 | No | 70 | 20.0
2 | 40 | No | 80 | 30.4 | 8 | 30 | Yes | 60 | 78.4
3 | 60 | Yes | 70 | 100.4 | 9 | 40 | No | 60 | 30.4
4 | 80 | No | 70 | 75.4 | 10 | 50 | No | 80 | 42.9
5 | 100 | Yes | 60 | 106.6 | 11 | 80 | Yes | 70 | 106.6
6 | 100 | Yes | 60 | 106.6 | 12 | 100 | Yes | 60 | 106.6

ST: Survival Time, PST: Pseudo-Survival Time

The survival curves of the original data are shown in Figure 1, and the
survival curves of the PST are presented in Figure 2. From these figures, note
that the mean survival time for both groups has increased through the use of
PST.

Figure 1: Survival curve of original survival data of example dataset

\+ means censored.

Figure 2: Survival curve of pseudo-survival data of example dataset

The difference in RMSTs between the two groups (treatment group versus control
group) until the cut-off point of $100$ weeks is summarized in Table 2. The
difference in the two RMSTs does not change under Andersen’s method, even if
age (covariate) is included, because the average age is equal in both groups
when censored patients are included. However, the average age of patients for
whom the event was observed is 75 (years) in the actual group and 70 (years)
in the placebo group. Tian’s method reduces the difference in the two RMSTs by
adjusting the covariate. The difference between the groups can be reduced by
considering the censored information compared with that in the naive RMSTs.

Table 2: Differences in RMSTs in the example dataset

Parameter | Naive | Andersen | Tian | Andersen (age) | Tian (age)
---|---|---|---|---|---
Diff | 29.1 | 18.8 | 25.0 | 18.8 | 24.5
Age | - | - | - | -2.1 | -2.3

Andersen/Tian: statistical model includes only the treatment group as a fixed
effect. Andersen (age)/Tian (age): statistical model includes the treatment
group as a fixed effect and age as a covariate.

We also consider two additional cases in which the age data are changed to
confirm the influence caused by age. In the first case, there is no difference
in average age between groups for all patients, including those censored, but
there is a large difference in average age between groups including only
patients for whom the event occurs. In the second case, there is a difference
in average age for all patients, including those censored. Details of
additional analyses are provided in the appendix. For the first case,
Andersen’s method does not change the RMST difference results, while Tian’s
method increases the difference. In the second case, Andersen’s method reduces
the difference, whereas Tian’s method does not change the difference.

The example data and analysis program code are presented in the Supplemental
Material.

##  3 Simulation Study

We now present the results of simulations to evaluate the methods for deriving
the RMST with covariate adjustments.

###  3.1 Simulation setup

We consider 16 simulation scenarios to evaluate the performance of Andersen’s
and Tian’s methods in terms of the statistical power for identifying the
difference in RMSTs between treatment and control groups. Eight scenarios are
based on the case where the PH assumption between the two groups holds. In
scenario 1 (PH-S1), the true survival curves are consistent across the two
groups, allowing us to confirm the type-1 error rate. In PH-S2, the percentage
of censored data in the treatment group is higher than in the control group.
This enables us to confirm the performance of Andersen’s and Tian’s methods,
which handle censored data differently. In PH-S3, the survival curve of the
treatment group is better than in the placebo group. In PH-S4, the percentage
of censored data in the treatment group is higher than in the control group.
Even in randomized trials, the covariates may differ between the two groups.
Thus, in PH-S5–PH-S8, the averages of the true covariates differ between the
two groups, allowing us to confirm that the covariates can be adjusted when
differences occur. The remaining eight scenarios consider settings in which
the slope of the survival function changes at a certain time, called a change
point (CP). In CP-S1–CP-S4, the two survival functions intersect. In
CP-S5–CP-S8, the control group has many events up to the CP, and few events
after the CP.

For PH-S1–PH-S8, we use the regression model below to generate the simulation
data. Let the survival times $T$ be given by

$\displaystyle\log
T=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+\log(\varepsilon),$ (10)

where $\beta_{0}$ is the intercept parameter, $\beta_{1}$ is a parameter
reflecting the difference between treatment ($x_{1}=1$) and control
($x_{1}=0)$, $\beta_{2}$ is a parameter reflecting the covariate, $x_{2}$ is
normally distributed with variance $20$, and $\varepsilon$ is a random
variable obeying a Weibull distribution.

For CP-S1–CP-S8, we use the regression model with CPs to generate the
simulation data. Let the survival times $T$ be given by

$\displaystyle\log
T=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+\beta_{3}x_{1}x_{3}+\log(\varepsilon),\,\,\,\,\,x_{3}=\left\\{\begin{array}[]{cc}a&(x_{2}\leq
50\mbox{ }\mathrm{(change\,point)})\\\ b&(x_{2}>50\mbox{
}\mathrm{(change\,point)})\end{array}\right.$ (13)

where $\beta_{3}$ is a covariate parameter with CP = $50$. The parameter
settings are presented in Table 3 and the related survival curves are shown in
Figure 3. We generate 10,000 datasets with $n=10,20,\dots,90,100$ subjects for
each scenario. We determine the power to detect the differences in RMST
between groups without adjusting for covariates, known as the naive approach.
For Andersen’s and Tian’s methods, the statistical model for detecting group
differences is as follows:

$\displaystyle\exp({\gamma}_{1}+x_{1}{\gamma}_{2}+x_{2}{\gamma}_{3}).$ (14)

When the p-value of ${\hat{\gamma}}_{2}$ for the null hypothesis $H_{0}:$
${\gamma}_{2}=0$ is less than $0.05$, we consider the group differences to be
detected. For reference, we also determine the power of Cox’s PH model to
detect differences between the groups. The parametric function of Cox’s PH
model is given by (14).

Table 3: Parameter settings of simulation scenarios

|  |  |  |  |  | True mean of $x_{2}$ |  |  | Percentage of censored data
---|---|---|---|---|---|---|---|---|---
Scenario | $\tau$ | $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | Treatment | Control | $a$ | $b$ | Treatment | Control
PH-S1 | 0.2 | 2 | 0 | -0.1 | - | 50 | 50 | - | - | 0.1 | 0.1
PH-S2 | 0.2 | 2 | 0 | -0.1 | - | 50 | 50 | - | - | 0.4 | 0.1
PH-S3 | 0.2 | 2 | 0.2 | -0.1 | - | 50 | 50 | - | - | 0.1 | 0.1
PH-S4 | 0.2 | 2 | 0.2 | -0.1 | - | 50 | 50 | - | - | 0.4 | 0.1
PH-S5 | 0.2 | 2 | 0 | -0.1 | - | 48 | 52 | - | - | 0.1 | 0.1
PH-S6 | 0.2 | 2 | 0 | -0.1 | - | 48 | 52 | - | - | 0.4 | 0.1
PH-S7 | 0.2 | 2 | 0.2 | -0.1 | - | 48 | 52 | - | - | 0.1 | 0.1
PH-S8 | 0.2 | 2 | 0.2 | -0.1 | - | 48 | 52 | - | - | 0.4 | 0.1
CP-S1 | 0.5 | 2 | -0.5 | -0.1 | 2 | 50 | 50 | 1 | -0.5 | 0.1 | 0.1
CP-S2 | 0.5 | 2 | -0.5 | -0.1 | 2 | 50 | 50 | 1 | -0.5 | 0.4 | 0.1
CP-S3 | 0.5 | 2 | -0.5 | -0.1 | 2 | 48 | 52 | 1 | -0.5 | 0.1 | 0.1
CP-S4 | 0.5 | 2 | -0.5 | -0.1 | 2 | 48 | 52 | 1 | -0.5 | 0.4 | 0.1
CP-S5 | 0.6 | 4 | -0.5 | -0.1 | 3 | 50 | 50 | 1 | -0.5 | 0.1 | 0.1
CP-S6 | 0.6 | 4 | -0.5 | -0.1 | 3 | 50 | 50 | 1 | -0.5 | 0.4 | 0.1
CP-S7 | 0.6 | 4 | -0.5 | -0.1 | 3 | 48 | 52 | 1 | -0.5 | 0.1 | 0.1
CP-S8 | 0.6 | 4 | -0.5 | -0.1 | 3 | 48 | 52 | 1 | -0.5 | 0.4 | 0.1

For PH-S2–PH-S8, changes from PH-S1 are shown in red. For CP-S2–CP-S4, changes
from CP-S1 are shown in red. For CP-S6–CP-S8, changes from CP-S5 are shown in
red.

Figure 3: Survival curves of simulation scenarios

###  3.2 Simulation results

We present the results of the 16 simulation scenarios in Figure 4. For PH-S1,
the type-1 error rate can be controlled at the nominal significance level by
all methods. For PH-S2, Tian can control the type-1 error rate. However, the
type-1 error rates of Andersen, Cox and naive are inflated as the number of
subjects increases. For PH-S3, the powers of Andersen, Tian, and Cox adjusting
for covariates are higher than naive. For PH-S4, Andersen, Cox and naive have
the highest power but Tian had slightly lower power. For PH-S5, the Andersen,
Tian, and Cox models give a nominal significance level, whereas the naive
approach produces an inflated type-1 error rate as the number of subjects
increases because of the lack of adjustment in the covariate. For PH-S6, Tian
can control the type 1 error rate but other methods inflated the type-1 error
rate as the number of subjects increases. For PH-S7, the powers of Andersen,
Tian, and Cox adjusting for covariates are higher than naive. For the PH-S8,
Andersen and Cox have the highest power, Tian had slightly lower power, and
naive does not increase in power even when the number of subjects was
increased.

For the CP scenarios, the power of all methods tends to rise as the number of
subjects increases. Notably, all methods have very high power in CP-S4, and
similarly for CP-S5, except for Cox’s model, in which the power increase more
slowly.

Figure 4: Power of simulation scenarios

##  4 Re-analysis of Data from Two Clinical Trials

To evaluate Andersen’s and Tian’s methods, we re-analyzed the data of two
actual clinical trials using the RMSTREG PROCEDURE in SAS. The SAS analysis
program code is included in the Supplemental Material.

###  4.1 Re-analysis of clinical trial data in primary biliary cholangitis

We evaluated the performance of Andersen’s and Tian’s methods using the
primary biliary cholangitis (PBC) dataset in the survival package of R [13].
This dataset comes from a randomized placebo-controlled trial of primary
biliary cholangitis conducted between 1974 and 1984. The original PBC data
include non-randomized patients. Considering only randomized patients, we
selected 134 patients (72: D-penicillamine group, 62: placebo group) without
hepatomegaly or enlarged livers. The percentage of censoring in the
D-penicillamine group is 72.22% (52/72) and the percentage of censoring in the
placebo group is 74.19% (46/62). The survival curves are shown in Figure 5.

Figure 5: Survival curves of PBC data

P: Placebo group, D: D-penicillamine group, + denotes censored.

The covariates are edema (0: no edema, 0.5: untreated or successfully treated,
1: edema despite diuretic therapy), serum bilirubin (mg/dl), serum albumin
(mg/dl), standardized blood clotting time (protime), and age (years). The
restricted time was 12.34 yrs. The link function is set to identity to
understand the interpretation of the analysis results easily.

In the re-analysis results, the naive RMST $\pm$ standard error is $8.88\pm
0.62$ in the placebo group and $7.96\pm 0.61$ in the D-penicillamine group,
with a $-0.92\pm 0.54$ difference between the two groups (D-penicillamine
group $-$ placebo group). The confidence interval (p-value) is $[-1.97,0.10]$
$(0.086$). The difference in RMST between groups adjusted by Andersen’s method
is $-0.56\pm 0.67$ $[-1.88,0.76]$ $(0.405)$, and that between groups adjusted
by Tian’s method is $-0.93\pm 0.43$ $[-1.78,-0.09]$ $(0.030)$. Table 4
presents the estimators of all model parameters. The survival curves for the
PST are shown in Figure 7 in the appendix. We can confirm that the survival
curves of the original data and the pseudo-data are different.

Table 4: Analysis results of Andersen’s and Tian’s methods in PBC data

| Andersen | Tian
---|---|---
Parameter | Mean | SE | [CI] (p-value) | Mean | SE | [CI] (p-value)
Intercept | 20.20 | 6.49 | [7.47, 32.93] (0.002) | -3.39 | 3.70 | [-10.63, 3.86] (0.360)
Treatment group | -0.56 | 0.67 | [-1.88, 0.76] (0.405) | -0.93 | 0.43 | [-1.78, -0.09] (0.030)
Edema (1) | -1.29 | 2.58 | [-6.35, 3.77] (0.618) | -5.70 | 1.70 | [-9.04, -2.35] (0.001)
Edema (0.5) | 0.39 | 1.68 | [-2.90, 3.68] (0.814) | -4.40 | 0.52 | [-5.42, -3.38] (0.000)
Bilirubin | -0.49 | 0.14 | [-0.76, -0.23] (0.000) | -0.11 | 0.11 | [-0.32, 0.11] (0.332)
Albumin | 0.24 | 1.08 | [-1.87, 2.35] (0.823) | 0.89 | 0.57 | [-0.24, 2.01] (0.122)
Protime | -0.48 | 0.51 | [-1.48, 0.53] (0.351) | 0.50 | 0.23 | [0.04, 0.96] (0.033)
Age | -0.11 | 0.03 | [-0.17, -0.04] (0.002) | 0.05 | 0.02 | [0.00, 0.09] (0.056)

SE: Standard Error. CI: $95\%$ Confidence Interval. Treatment group and edema
are fixed effects.

###  4.2 Re-analysis of clinical trial data in malignant glioma

We also evaluated the performance of Andersen’s and Tian’s methods using the
randomized placebo-controlled clinical trial data for malignant glioma
reported by Brem et al. [14]. The analysis data include 110 patients in the
treatment group and 112 patients in the placebo group. The percentage of
censoring in the treatment group is 6.36% (7/110) and the percentage of
censoring in the placebo group is 7.14% (8/112). The survival curves are shown
in Figure 6.

Figure 6: Survival curves of malignant glioma data

Pl: Placebo group, Po: Polymer group, + denotes censored.

The covariates are age (years), years from diagnosis, Karnofsky performance
score (0: ¡ 70, 1: $\leq$ 70), race (0: other, 1: white), sex (0: female, 1:
male), nitro (0: no, 1: yes), tumor histopathology at implantation (path) (1:
glioblastoma, 2: anaplastic astrocytoma, 3: oligodendroglioma, 4: other), and
grade (0: quiescent, 1: active). The restricted time point was 150 week. The
link function is set to identity to understand the interpretation of the
analysis results easily.

The naive RMST $\pm$ standard error is $40.51\pm 2.03$ in the placebo group
and $47.40\pm 2.18$ in the polymer group, with a $6.88\pm 6.29$ difference
between the two groups (polymer group $-$ placebo group). The confidence
interval (p-value) is $[-5.45,19.21]$ $(0.274$). The difference in RMST
between groups adjusted by Andersen’s method is $-0.56\pm 0.67$ $[-1.88,0.76]$
$(0.405)$ and that between groups adjusted by Tian’s method is $-0.93\pm 0.43$
$[-1.78,-0.09]$ $(0.030)$. Table 5 gives the estimators of all model
parameters. The survival curves for the PST are shown in Figure 8 in the
appendix. We can confirm that the survival curves of the original data and the
pseudo-data are similar.

Table 5: Analysis results of Andersen’s and Tian’s methods in malignant glioma
data

| Andersen | Tian
---|---|---
Parameter | Mean | SE | [CI] (p-value) | Mean | SE | [CI] (p-value)
Intercept | 75.12 | 21.12 | [33.74, 116.51] (0.000) | 47.18 | 20.57 | [6.87, 87.49] (0.022)
Treatment group | 7.41 | 5.77 | [-3.89, 18.71] (0.199) | 12.10 | 5.33 | [1.64, 22.55] (0.023)
Age | -0.57 | 0.23 | [-1.01, -0.13] (0.012) | -0.13 | 0.19 | [-0.49, 0.24] (0.499)
Dx | 1.70 | 1.48 | [-1.21, 4.61] (0.252) | 1.71 | 1.10 | [-0.44, 3.87] (0.120)
PS | 12.69 | 5.25 | [2.41, 22.97] (0.016) | 8.75 | 5.41 | [-1.85, 19.36] (0.106)
Race | -17.54 | 13.46 | [-43.93, 8.84] (0.193) | -8.47 | 9.76 | [-27.61, 10.66] (0.385)
Sex | 11.60 | 5.35 | [1.12, 22.08] (0.030) | -0.29 | 5.2 | [-10.48, 9.9] (0.955)
Nitro | -9.98 | 5.50 | [-20.76, 0.8] (0.070) | -15.6 | 4.01 | [-23.47, -7.73] (0.000)
Path (4) | 47.74 | 20.07 | [8.4, 87.07] (0.017) | 22.46 | 11.49 | [-0.06, 44.98] (0.051)
Path (3) | 37.19 | 13.23 | [11.26, 63.12] (0.005) | 22.45 | 13.69 | [-4.37, 49.28] (0.101)
Path (2) | 19.30 | 9.22 | [1.22, 37.37] (0.036) | 15.14 | 7.57 | [0.31, 29.98] (0.045)
Grade | -15.68 | 11.74 | [-38.7, 7.34] (0.182) | -9.23 | 9.05 | [-26.96, 8.5] (0.307)

SE: Standard Error. CI: $95\%$ Confidence Interval. Treatment group and edema
are fixed effects.

##  5 Discussion

We evaluated the performance of Andersen’s and Tian’s baseline covariate
adjustment methods because we could not decide which method to include in a
statistical analysis plan for a clinical trial. Andersen’s method calculates
the PST using a leave-one-out approach for all subjects, including censored
subjects, and performs regression analysis using all pseudo-values. Tian’s
method uses IPCW to estimate the parameters of the regression model, and also
considers censoring. Therefore, there is a significant difference between the
two methods in the process of adjusting the covariates. In Andersen’s method,
the interpretation of pseudo-values is difficult. We consider Tian’s method to
be more intuitive than Andersen’s method. In this paper, we confirmed the
intuitive behavior of Andersen’s and Tian’s methods using a simple example
dataset.

Through a simulation study, we confirmed that Tian’s method has higher power
than Andersen’s method. Tian’s method is able to control the type-1 error rate
independently of censor bias, whereas Andersen’s method and Cox’s PH model
inflate the type-1 error rate. For scenarios PH-S3 and PH-S6, Cox’s PH model
only has slightly higher power than the RMST. Royston and Parmar [11], Tian et
al. [15], and Huang and Kuan [16] also found that the log-rank test had
slightly higher power than the RMST under the PH assumption. In the CP
scenarios, where the PH assumption does not hold, the RMST is more powerful.
In PH-S4, PH-S8, CP-S2, and CP-S4, when the treatment group has a higher
proportion of censored data, Andersen’s method and Cox’s PH model exhibit the
highest power. However, because the censors are biased toward the treatment
group, Andersen’s method and Cox’s PH model give an inflated type-1 error
rate. Hence, we consider Tian’s method to be more appropriate in terms of
retaining statistical power in certain situations.

We re-analyzed clinical PBC trial data provided by the R survival package and
malignant glioma trial data provided by Piantadosi [17]. For the PBC data,
Andersen’s method implies that bilirubin and age affect the RMST. In contrast,
Tian’s method implies that the treatment group, edema (1, 0.5), and protime
affect the RMST. These differences in results can be attributed to differences
in methodology. The percentage of censoring in the placebo group is 74.19%
(46/62), and the percentage of censoring in the D-penicillamine group is
72.22% (52/72). Thus, the high incidence of censoring may have contributed to
the difference in results. For Andersen’s method, there is no difference in
bilirubin between groups and an age difference of 2.9 (years) (47.4 in the
placebo group and 50.3 in the D-penicillamine group). Because the
D-penicillamine group is older and the estimate of the age parameter is
negative, the difference is smaller than the adjustment for covariates. For
Tian’s method, there is no difference in protime between groups. In the
placebo group, edema (0) is 87.5% (=14/16) and edema (1) is 12.5% (=2/16),
while in the D-penicillamine group, edema (0) is 75% (=15/20), edema (0.5) is
10% (=2/20), and edema (1) is 15% (=3/20). Although there is a small
difference in the edema distribution between groups, little adjustment is made
for differences in RMST because of the small number of observed events. We
discuss the adjusted difference in RMSTs in more detail in the appendix.

For the malignant glioma data, Andersen’s method implies that age, PS, sex,
and path (4, 3, and 2) affect the RMST. Tian’s method implies that the
treatment group, nitro, and path (4 and 2) affect the RMST. The reason for the
difference in estimation results is the different estimation methods. For
Andersen’s method, the difference in RMST is almost the same because there is
no significant difference in the background information. For Tian’s method,
there is no difference in nitro between groups; for the placebo group, path
(1) is 68.3%, path (2) is 14.4%, path (3) is 16.3%, and path (4) is 1.0%,
while in the polymer group, path (1) is 72.8%, path (2) is 12.6%, path (3) is
9.7%, and path (4) is 4.9%. The proportion of path (2) and path (3) in the
polymer group is lower than in the placebo group. The estimates of the path
(2) and path (3) parameters are positive, which widens the difference in RMST
between the groups. We discuss the adjusted difference in RMSTs in more detail
in the appendix.

In summary, the use of PSTs makes Andersen’s method difficult to interpret,
while Tian’s method retains its power in various situations. Even when group
differences occur with respect to the censor incidence rates and covariates,
Tian’s method maintains a nominal significance level for the type-1 error
rate. Using actual trial data, we have confirmed that the PSTs are different
from the actual data in the PBC case because many data are censored.
Andersen’s method should be used when the purpose is to adjust for group
differences in RMST using the background of all subjects, including the
censored subjects, rather than to adjust the backgrounds of only those
subjects for whom the event occurred. We conclude that Tian’s method can be
recommended.

Author Contributions. KH developed the method and prepared the simulation
program code, JM reviewed the overall study design, and MK wrote the bulk of
the paper. MK also analyzed the example data, simulation datasets, and
clinical trial data.

Acknowledgements. MK would like to thank Associate Professor Hisashi Noma for
his encouragement and helpful suggestions.

## References

  * [1] Hossein Borghaei, Luis Paz-Ares, Leora Horn, David R Spigel, Martin Steins, Neal E Ready, Laura Q Chow, Everett E Vokes, Enriqueta Felip, Esther Holgado, et al.  Nivolumab versus docetaxel in advanced nonsquamous non–small-cell lung cancer.  New England Journal of Medicine, 373(17):1627–1639, 2015.
  * [2] Caroline Robert, Georgina V Long, Benjamin Brady, Caroline Dutriaux, Michele Maio, Laurent Mortier, Jessica C Hassel, Piotr Rutkowski, Catriona McNeil, Ewa Kalinka-Warzocha, et al.  Nivolumab in previously untreated melanoma without braf mutation.  New England journal of medicine, 372(4):320–330, 2015.
  * [3] Spotswood L Spruance, Julia E Reid, Michael Grace, and Matthew Samore.  Hazard ratio in clinical trials.  Antimicrobial agents and chemotherapy, 48(8):2787–2792, 2004.
  * [4] Theodore Karrison.  Restricted mean life with adjustment for covariates.  Journal of the American Statistical Association, 82(400):1169–1176, 1987.
  * [5] David M Zucker.  Restricted mean life with covariates: modification and extension of a useful survival analysis method.  Journal of the American Statistical Association, 93(442):702–709, 1998.
  * [6] Per Kragh Andersen, Mette Gerster Hansen, and John P Klein.  Regression analysis of restricted mean survival time based on pseudo-observations.  Lifetime Data Anal., 10:335––350, 2004.
  * [7] Lu Tian, Lihui Zhao, and LJ Wei.  Predicting the restricted mean event time with the subject’s baseline covariates in survival analysis.  Biostatistics, 15(2):222–233, 2014.
  * [8] Chengfeng Zhang, Hongji Wu, Baoyi Huang, Hao Yuan, Yawen Hou, and Zheng Chen.  Restricted mean survival time regression model with time-dependent covariates.  Statistics in Medicine, 41:4081–4090, 2022.
  * [9] Federico Ambrogi, Simona Iacobelli, and Per Kragh Andersen.  Analyzing differences between restricted mean survival time curves using pseudo-values.  BMC medical research methodology, 22(1):1–12, 2022.
  * [10] Theodore Karrison and Masha Kocherginsky.  Restricted mean survival time: Does covariate adjustment improve precision in randomized clinical trials?  Clinical Trials, 15(2):178–188, 2018.
  * [11] Patrick Royston and Mahesh KB Parmar.  Restricted mean survival time: an alternative to the hazard ratio for the design and analysis of randomized trials with a time-to-event outcome.  BMC medical research methodology, 13(1):1–15, 2013.
  * [12] Per Kragh Andersen, John P. Klein, and Susanne Rosthøj.  Generalised linear models for correlated pseudo-observations, with applications to multi-state models.  Biometrika, 90(1):15–27, 2003.
  * [13] Terry M Therneau and Thomas Lumley.  Package ‘survival’.  R Top Doc, 128(10):28–33, 2015.
  * [14] Henry Brem, S Piantadosi, PC Burger, M Walker, R Selker, NA Vick, K Black, M Sisti, S Brem, G Mohr, et al.  Placebo-controlled trial of safety and efficacy of intraoperative controlled delivery by biodegradable polymers of chemotherapy for recurrent gliomas.  The Lancet, 345(8956):1008–1012, 1995.
  * [15] Lu Tian, Haoda Fu, Stephen J Ruberg, Hajime Uno, and Lee-Jen Wei.  Efficiency of two sample tests via the restricted mean survival time for analyzing event time observations.  Biometrics, 74(2):694–702, 2018.
  * [16] Bo Huang and Pei-Fen Kuan.  Comparison of the restricted mean survival time with the hazard ratio in superiority trials with a time-to-event end point.  Pharmaceutical statistics, 17(3):202–213, 2018.
  * [17] Steven Piantadosi.  Clinical trials: a methodologic perspective.  John Wiley & Sons, 2017.

##  Appendix A Appendix

###  A.1 Generalized estimating equations of Andersen’s method

[Case where the link function is the identity function]

The linear model is expressed as

$\displaystyle{\widehat{\mu}}_{i}={\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i}$
(15)

$\displaystyle{\text{\boldmath$U$}}_{PV}({\text{\boldmath$\beta$}})$
$\displaystyle=\sum_{i=1}^{n}{\text{\boldmath$U$}}_{i}({\text{\boldmath$\beta$}})$
$\displaystyle=\sum_{i=1}^{n}\left(\\{{\widehat{\mu}}_{i}-{\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i}\\}{\text{\boldmath$Z$}}_{i}\right)$
$\displaystyle=0.$ (16)

The estimator ${\widehat{\text{\boldmath$\beta$}}}$ can be obtained as

$\displaystyle{\widehat{\text{\boldmath$\beta$}}}=\sum_{i=1}^{n}{\widehat{\mu}}_{i}\left(\sum_{j=1}^{n}{\text{\boldmath$Z$}}_{j}{\text{\boldmath$Z$}}_{j}^{T}\right)^{-1}{\text{\boldmath$Z$}}_{i}.$
(17)

[Case where the link function is the log function]

The log-linear model is expressed as

$\displaystyle\log({\widehat{\mu}}_{i})={\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i}$
(18)

The estimator ${\widehat{\text{\boldmath$\beta$}}}$ can be obtained from the
following equation. The equations below are not in closed-form and need to be
estimated using numerical analysis.

$\displaystyle{\text{\boldmath$U$}}_{PV}({\text{\boldmath$\beta$}})$
$\displaystyle=\sum_{i=1}^{n}{\text{\boldmath$U$}}_{i}({\text{\boldmath$\beta$}})$
$\displaystyle=\sum_{i=1}^{n}\left(\\{{\widehat{\mu}}_{i}-\exp({\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i})\\}\exp({\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i}){\text{\boldmath$Z$}}_{i}\right)$
$\displaystyle=0.$ (19)

###  A.2 Estimating equations of Tian’s method

[Case where the link function is the identity function]

The linear model is expressed as

$\displaystyle X_{i}={\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i}$
(20)

$\displaystyle{\text{\boldmath$U$}}_{IPCW}({\text{\boldmath$\beta$}})=\sum_{i=1}^{n}\frac{I(X_{i}\leq
C_{i})\\{X_{i}-{\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i}\\}}{\hat{S}(X_{i})}{\text{\boldmath$Z$}}_{i}=0.$
(21)

The estimator ${\widehat{\text{\boldmath$\beta$}}}$ can be obtained as

$\displaystyle{\widehat{\text{\boldmath$\beta$}}}=\sum_{i=1}^{n}\frac{I(X_{i}\leq
C_{i})\ X_{i}}{\hat{S}(X_{i})}\left(\sum_{j=1}^{n}\frac{I(X_{j}\leq
C_{j})}{\hat{S}(X_{j})}{\text{\boldmath$Z$}}_{j}{\text{\boldmath$Z$}}_{j}^{T}\right)^{-1}{\text{\boldmath$Z$}}_{i}.$
(22)

[Case where the link function is the log function]

The log-linear model is expressed as

$\displaystyle\log(X_{i})={\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i}.$
(23)

The estimator ${\widehat{\text{\boldmath$\beta$}}}$ can be obtained from the
following equation. This equation is not in closed-form and needs to be
estimated using numerical analysis.

$\displaystyle{\text{\boldmath$U$}}_{IPCW}({\text{\boldmath$\beta$}})=\sum_{i=1}^{n}\frac{I(X_{i}\leq
C_{i})\\{X_{i}-\exp({\text{\boldmath$\beta$}}^{T}{\text{\boldmath$Z$}}_{i})\\}}{\hat{S}(X_{i})}{\text{\boldmath$Z$}}_{i}=0.$
(24)

###  A.3 Additional analyses of modified example dataset

We considered two additional cases in which the age data were changed to
confirm the influence caused by age. In the first case, there is no difference
in average age between groups for all patients, including those censored, but
there is a large difference in average age between groups of patients for whom
events occur. The second case considers a difference in average age for all
patients, including those censored.

In the first case, the average age of patients for whom the event occurred is
75 for the treatment group and 67 for the control group. Detailed data are
presented in Table 6. The red font indicates changed data.

Table 6: Example dataset

Treatment group | Control group
---|---
ID | ST | Censor | Age | PST | ID | ST | Censor | Age | PST
1 | 20 | Yes | 60 | 78.4 | 7 | 20 | No | 60 | 20.0
2 | 40 | No | 80 | 30.4 | 8 | 30 | Yes | 70 | 78.4
3 | 60 | Yes | 70 | 100.4 | 9 | 40 | No | 60 | 30.4
4 | 80 | No | 70 | 75.4 | 10 | 50 | No | 80 | 42.9
5 | 100 | Yes | 60 | 106.6 | 11 | 80 | Yes | 70 | 106.6
6 | 100 | Yes | 60 | 106.6 | 12 | 100 | Yes | 60 | 106.6

ST: Survival Time, PST: Pseudo-Survival Time

The analysis results are given in Table 7. For Tian’s method, the difference
in RMST is increased by adjusting for the age because of the higher average
age in the treatment group.

Table 7: Differences in RMSTs of example dataset

Parameter | Naive | Andersen | Tian | Andersen (age) | Tian (age)
---|---|---|---|---|---
Diff | 29.1 | 18.8 | 25.0 | 18.8 | 28.1
Age | - | - | - | -1.2 | -1.7

Andersen/Tian: statistical model includes only the treatment group as a fixed
effect. Andersen (age)/Tian (age): statistical model includes the treatment
group as a fixed effect and age as a covariate.

In the second case, the average age of all patients is 67 for the treatment
group and 68 for the control group. Detailed data are presented in Table 8.
The red font indicates the changed data.

Table 8: Example dataset

Treatment group | Control group
---|---
ID | ST | Censor | Age | PST | ID | ST | Censor | Age | PST
1 | 20 | Yes | 60 | 78.4 | 7 | 20 | No | 70 | 20.0
2 | 40 | No | 80 | 30.4 | 8 | 30 | Yes | 70 | 78.4
3 | 60 | Yes | 70 | 100.4 | 9 | 40 | No | 60 | 30.4
4 | 80 | No | 70 | 75.4 | 10 | 50 | No | 80 | 42.9
5 | 100 | Yes | 60 | 106.6 | 11 | 80 | Yes | 70 | 106.6
6 | 100 | Yes | 60 | 106.6 | 12 | 100 | Yes | 60 | 106.6

ST: Survival Time, PST: Pseudo-Survival Time

The analysis results are given in Table 9. For Andersen’s method, the
difference in RMST is reduced by adjusting for age because of the higher
average age in the control group.

Table 9: Differences in RMSTs of example dataset

Parameter | Naive | Andersen | Tian | Andersen (age) | Tian (age)
---|---|---|---|---|---
Diff | 29.1 | 18.8 | 25.0 | 15.4 | 24.5
Age | - | - | - | -2.0 | -2.3

Andersen/Tian: statistical model includes only the treatment group as a fixed
effect. Andersen (age)/Tian (age): statistical model includes the treatment
group as a fixed effect and age as a covariate.

###  A.4 Survival curves of PSTs in two actual clinical trials

The survival curves of the pseudo-values for PBC clinical trial data are shown
in Figure 7.

Figure 7: Survival curves of pseudo-values of PBC clinical trial data

The survival curves of PSTs of malignant glioma clinical trial data are shown
in Figure 8.

Figure 8: Survival curve of pseudo-values of malignant glioma clinical trial
data

###  A.5 Detailed analysis of data from two actual clinical trials

For the PBC data, Andersen’s method implied that bilirubin and age affect the
RMST, whereas Tian’s method implied that the treatment group, edema (1, 0.5),
and protime affect the RMST. These differences in results were attributed to
differences in methodology. The percentage of censoring in the placebo group
is 74.19% (46/62) and the percentage of censoring in the D-penicillamine group
is 72.22% (52/72). Thus, the high incidence of censoring may have contributed
to the difference in results. For Andersen’s method, the difference between
groups for bilirubin and age, for which the confidence intervals did not
include 0, was examined for all subjects, including the censored subjects.
Bilirubin was measured at 2.1 in the placebo group and 1.9 in the
D-penicillamine group, a difference of 0.2. The average age was 47.4 in the
placebo group and 50.3 in the D-penicillamine group, a difference of 2.9.
Because the D-penicillamine group was older and the estimate of the age
parameter was negative, the difference was smaller than the adjustment for
covariates. For Tian’s method, the difference between groups for protime and
the distribution of edema, for which the confidence intervals did not include
0, was examined for all subjects for whom the event occurred. Protime was
equal to 11.1 in both the placebo and D-penicillamine groups. For the
distribution of edema, edema (0) was 87.5% (=14/16) and edema (1) was 12.5%
(=2/16) in the placebo group, and edema (0) was 75% (=15/20), edema (0.5) was
10% (=2/20), and edema (1) was 15% (=3/20) in the D-penicillamine group.
Although there is a small difference in the edema distribution between groups,
little adjustment was made for differences in RMST because of the small number
of observed events.

For the malignant glioma data, Andersen’s method implied that age, PS, sex,
and path (4, 3, and 2) affect the RMST, whereas Tian’s method implied that the
treatment group, nitro, and path (4 and 2) affect the RMST. The reason for the
difference in estimation results was the different estimation methods. For
Andersen’s method, the difference between groups for age, PS, sex, and path,
for which the confidence intervals did not include 0, was examined for all
subjects, including the censored subjects. The average age was 47.6 in the
placebo group and 48.1 in the polymer group, a difference of 0.5. PS (1) was
50.0% in the placebo group and 55.5% in the polymer group, a difference of
5.5%. The proportion of males was 61.6% in the placebo group and 67.3% in the
polymer group, a difference of 5.7%. In the placebo group, Path (1) was 65.2%,
path (2) was 14.3%, path (3) was 17.9%, and path (4) was 2.7%, whereas in the
polymer group, path (1) was 69.1%, path (2) was 12.7%, path (3) was 13.6%, and
path (4) was 4.5%. The difference in RMST was almost the same because there
was no significant difference in the background information. For Tian’s
method, the difference between groups for nitro and the proportion of path
(2), for which the confidence intervals did not include 0, was examined for
all subjects for whom the event occurred. Nitro was measured at 45.2% in the
placebo group and 49.5% in the polymer group, a difference of 4.3%. For the
placebo group, path (1) was 68.3%, path (2) was 14.4%, path (3) was 16.3%, and
path (4) was 1.0%, whereas in the polymer group, path (1) was 72.8%, path (2)
was 12.6%, path (3) was 9.7%, and path (4) was 4.9%. The proportions of path
(2) and path (3) in the polymer group were lower than in the placebo group.
The estimates of the path (2) and path (3) parameters were positive, which
widened the difference in RMST between the groups.