DENNIS DOBLER AND KATHRIN MÖLLENHOFF

*Corresponding author, Vrije Universiteit Amsterdam, De Boelelaan 1111, 1081 HV, The Netherlands.

Vrije Universiteit Amsterdam, De Boelelaan 1111, 1081 HV, The Netherlands

# A nonparametric relative treatment effect for direct comparisons of censored
paired survival outcomes

Dennis Dobler  Kathrin Möllenhoff  Department of Mathematics, Faculty of
Science, Vrije Universiteit Amsterdam, North-Holland, The Netherlands
Mathematical Institute, Heinrich-Heine-University Düsseldorf, North-Rhine
Westphalia, Germany<EMAIL_ADDRESS>
(26 April 2016; 6 June 2016; 6 June 2016)

###### Abstract

[Summary]A very classical problem in statistics is to test the stochastic
superiority of one distribution to another. However, many existing approaches
are developed for independent samples and, moreover, do not take censored data
into account. We develop a new estimand-driven method to compare the
effectiveness of two treatments in the context of right-censored survival data
with matched pairs. With the help of competing risks techniques, the so-called
relative treatment effect is estimated. It quantifies the probability that the
individual undergoing the first treatment survives the matched individual
undergoing the second treatment. Hypothesis tests and confidence intervals are
based on a studentized version of the estimator, where resampling-based
inference is established by means of a randomization method. In a simulation
study, we found that the developed test exhibits good power, when compared to
competitors which are actually testing the simpler null hypothesis of the
equality of both marginal survival functions. Finally, we apply the
methodology to a well-known benchmark data set from a trial with patients
suffering from with diabetic retinopathy.

###### keywords:

Estimand, matched pairs, nonparametric statistics, randomization, survival
analysis.

††articletype: Article Type

##  1 Introduction

Testing the stochastic superiority of one distribution to another is a very
classical problem in statistics. For two independent and fully observable
samples, the $t$-test, the median test, and the Mann-Whitney test are three
well-known solutions to it. Likewise, in the case of two dependent samples,
i.e., if data consist of (matched) pairs, the $t$-test, the sign test, and the
Wilcoxon signed rank test based on the pair-wise differences could be used.
Extensions to the case of two dependent and _censored_ samples are not obvious
but some have been developed in the literature for a few decades. The
developed test statistics are based on ranked absolute within-pair differences
of the possibly censored survival times 41, a difference of counting processes
37, differences of efficient scores 31, ranking all censored observations
separately from the uncensored ones 3, a combination of terms for different
censoring and pair-wise ordering patterns 9, integrals of scores with respect
to differences of the sample-specific Nelson-Aalen estimators for cumulative
hazards 26, a combination of frailty-based log-rank tests 30, using further
prioritized outcomes, i.e., additional data, if the primary survival endpoint
does not offer decisive pair-wise comparisons 33. Extensions for covariates
were also developed 23, 25. Reviews of methods for paired survival data
including additional references to other approaches and discussions are also
available42, 28.

While these existing approaches undoubtedly offer many good approaches for
powerful statistical inference in the two-sample problem, an additional,
easily interpretable quantification of the discrepancy between both samples is
usually not available. One notable exception is the popular win ratio method
of Pocock et al.33 which recently had been exploited through win odds in order
to take ties into account in the inference method8. The perhaps most
straightforward approach for such a quantification is to compare the survival
chances for both groups at a fixed time point. This would provide a very
limited, yet easy-to-interpret summary. A more global impression of the
difference between both samples could be obtained by integration over time,
leading to the (restricted) mean survival time. We will treat this topic in
another forthcoming paper.

Instead, we will pursue an approach which is motivated by another estimand
with a very clear interpretation: the _relative treatment effect_. In brief,
it describes the probability that the lifetime under Treatment 1 is bigger
than the lifetime under Treatment 2. If this is (significantly) different from
0.5, a solid statistical conclusion can be drawn about the treatment
efficacies. At the same time, it is a simple probability which is easy to
communicate.

In the present paper, we will develop a nonparametric methodology with an
emphasis on the following quality criteria:

  1. 1.

begin the research with a clear formulation of an estimand of interest;

  2. 2.

make only very few and weak assumptions for the method to work;

  3. 3.

in particular, no continuity of survival functions is needed, i.e.,
instantaneous hazard rates need not exist;

  4. 4.

guaranteed large sample properties;

  5. 5.

a good statistical reliability even for small samples, i.e., good control of
the type-I error rate and confidence level, as well as a good power and narrow
confidence intervals, respectively.

All of these points are of crucial importance, in particular in the light of
the ICH E9(R1)17 guidelines on estimands in trial analysis. At the beginning
of Section A.5.1 therein, it is written that: “An estimand for the effect of
treatment relative to a control will be estimated by comparing the outcomes in
a group of subjects on the treatment to those in a similar group of subjects
on the control. For a given estimand, an aligned method of analysis, or
estimator, should be implemented that is able to provide an estimate on which
reliable interpretation can be based. The method of analysis will also support
calculation of confidence intervals and tests for statistical significance. An
important consideration for whether an interpretable estimate will be
available is the extent of assumptions that need to be made in the analysis.”

These statements clarify that even the most powerful inference method might
not be the preferable one if other criteria are not met, e.g., if no
intelligible estimand is available. Also, we wish to point out that most of
the methods rely on strong assumptions such as the equality of censoring times
for both members of a pair41, 3 or the continuity of survival distributions41,
3, 9, 30. In this sense, the power of the test we will develop in this paper
is not criterion of greatest importance although a powerful method is of
course welcome.

This paper is structured as follows. First, we will introduce the relative
treatment effect and explain its estimation in Section 2. At the end of that
section, we will relate the present approach to some others from the
literature. Second, in Section 3, we will present a method to make statistical
inference based on a data re-randomization technique13. Third, we will
investigate the large sample properties of the new method and explore its
small sample performance by means of a simulation study in Section 4. Therein,
also the power of the test will be assessed in a comparison to some competitor
methods. Next, in Section 5, we will apply the methodology to a well-known
benchmark data set from a trial with patients suffering from with diabetic
retinopathy, described by Huster et al.24. We will conclude with a discussion
in Section 6. The Supplementary Material contains all proofs, additional
technical details, and additional simulation results.

##  2 Paired survival and the relative treatment effect

###  2.1 Model and notation

We denote by $(T_{1},T_{2})$ a bivariate random vector on a probability space
$(\Omega,\mathcal{A},P)$. Each $T_{j}$ $(j=1,2)$ stands for a survival time of
a patient who was randomized to receive treatment $j$. For instance, the pair
results from a matching of two individuals with a similar physiology. As a
consequence, we generally assume $T_{1}$ and $T_{2}$ to be dependent. There
exist many approaches for the estimation of the bivariate survival
distribution of $(T_{1},T_{2})$ in the literature; see the paper by Pruitt34
for a comparison of six methods that are able to handle bivariate right-
censored data; Dai et al.10 developed an estimator under the more general
assumption of bivariate left-truncated and right-censored data. We do not mean
to give an exhaustive list of references in that direction and point to the
references included in the two just mentioned papers. Instead of estimating
the bivariate survival function, we aim at estimating a meaningful summary of
it, i.e., a treatment effect measure.

Let us now prepare the introduction of our estimand of interest. Let $\tau>0$
denote the maximum follow-up time of a study in which the superiority of
Treatment 1 to Treatment 2 shall be analyzed. Here, superiority means that
Treatment 1 prolongs the survival times compared to Treatment 2. Hence, it
seems constructive to consider the following probability:

$\displaystyle\tilde{\theta}=P(T_{1}>T_{2})+\tfrac{1}{2}P(T_{1}=T_{2}).$ (1)

The second term is important to give equal credit to both treatments in the
case of equal outcomes. Furthermore, we say that Treatment 1 is preferable if
$\tilde{\theta}>0.5$. However, due to the maximum follow-up time $\tau$,
$\tilde{\theta}$ is not always estimable; instead, we will focus on the
estimand

$\displaystyle\theta=P(\min(T_{1},\tau)>\min(T_{2},\tau))+\tfrac{1}{2}P(\min(T_{1},\tau)=\min(T_{2},\tau)),$
(2)

which we will we call the _relative treatment effect_ from now on. Brunner and
Munzel7 introduced this concept for the case of two independent samples. Many
more research papers in this context emerged afterwards. Just to mention two,
Munzel and Brunner29 and Konietschke and Pauly27 considered to the case of two
dependent but fully observable samples. Finally, for technical reasons, we
assume that $P(T_{1}>\tau,T_{2}>\tau)>0$ and that
$P(T_{1}=\tau)=P(T_{2}=\tau)=0$. The latter is always achievable by
artificially increasing $\tau$ by a very small number. Throughout the paper,
we assume that $\theta\in(0,1)$, i.e., no perfect superiority of one treatment
over the other.

As mentioned before, censoring is omnipresent in many medical studies. We thus
assume that survival times are independently right-censored, i.e., it is only
possible to observe an event if it occurred before a so-called censoring time,
say $C_{1}$ and $C_{2}$, respectively. These are allowed to be dependent, but
$(C_{1},C_{2})$ and $(T_{1},T_{2})$ are assumed to be independent. As a
consequence, the actually observable data are of the form
$(X_{1},\delta_{1},X_{2},\delta_{2})$, where $X_{j}=\min(T_{j},C_{j},\tau)$
and $\delta_{j}=1\\{\min(T_{j},\tau)\leq C_{j}\\}$, $j=1,2$; here,
$1\\{\cdot\\}$ denotes the indicator function. For estimation of the second
probability in the relative treatment effect, we consider the case of
$X_{j}=\tau$ as uncensored. On a side note, Efron14 also assumes the largest
observation as uncensored in order to achieve a so-called “self-consistency”
property for the Kaplan-Meier estimator; see Section 7 therein.

In the following, we will assume without loss of generality that the censoring
times are continuously distributed prior to $\tau$. In the case of discrete
components in their distribution, ties can be broken by adding very small
positive random numbers to them. These random numbers can be chosen small
enough so that the order of all event times among the censoring times is not
altered. Also, all statistical procedures considered below are not affected by
these small modifications of the censoring times.

###  2.2 Transformation of paired survival data into competing risks data, and
estimation

Let us now describe our novel estimation approach for the relative treatment
effect. We assume that our data set consists of $n$ independently and
identically distributed data points of the kind described above, i.e., of
independent paired right-censored data. For a facilitated estimation, we
transform the paired survival data into a competing risks data set. Note that
our transformation is similar in spirit to the one used by Scheike, Holst, and
Hjelmborg36: therein, a transformation of a paired _competing risks_ data set
into a univariate one facilitated the estimation of a concordance function.
Our transformation works as follows: if for a pair

  * 1.

the first entry is _observed_ to fail before the second, an event of type 1
occurred;

  * 2.

the second entry is _observed_ to fail before the first, an event of type 2
occurred;

  * 3.

both entries are _observed_ to fail simultaneously, an event of type 3
occurred.

All other cases are labelled right-censored. In each case, the (censored)
event time is set to be the minimum of all four event and censoring times. We
summarize the thus obtained competing risks data set as
$(Z_{i},\varepsilon_{i})=(\min(\check{T}_{i},\check{C}_{i}),\check{\varepsilon}_{i}\cdot
1\\{\check{T}_{i}\leq\check{C}_{i}\\}),i=1,\dots,n$. Here, $\check{T}_{i}$ is
the minimum of both (potentially unobservable) event times of pair $i$ and
$\check{C}_{i}$ is the minimum of both (potentially unobservable) censoring
times of pair $i$. Note that the type of event $\check{\varepsilon}_{i}$ is
unobservable in the case of a censoring. More technical details about the
transformation are given in Section A in the Supplementary Material.

One important consequence of the subsequent Proposition 2.1 is that the Aalen-
Johansen estimators 1 of $F_{j}$, say $\widehat{F}_{j,n}$, $j=1,2,3$, are
available for estimating $\theta$:

$\displaystyle\widehat{\theta}_{n}=\widehat{F}_{2,n}(\tau)+\tfrac{1}{2}\widehat{F}_{3,n}(\tau).$
(3)

The specific structure of these Aalen-Johansen estimators are given in the
Supplementary Material. Many statistical properties of the Aalen-Johansen
estimator are well-known; cf. 5, Section IV.4.

######  Proposition 2.1.

  1. (a)

(Representation) The relative treatment effect can be written as
$\theta=F_{2}(\tau)+\tfrac{1}{2}F_{3}(\tau)$, where $F_{j}(t)=P(\check{T}\leq
t,\check{\varepsilon}=j)$, $t\in[0,\tau]$, denotes the $j$-th cumulative
incidence function, $j=1,2,3$.

  2. (b)

(Sufficiency) The above-described data reduction is sufficient for $\theta$.

  3. (c)

(Efficiency) $\widehat{\theta}_{n}$ is the nonparametric maximum likelihood
estimator (NPMLE) of $\theta$.

Let us now translate two of the most important statistical results about
Aalen-Johansen estimators to the relative treatment effect estimator, that is,
consistency and asymptotic normality. To this end, we denote convergences in
probability and in distribution by $\stackrel{{\scriptstyle p}}{{\to}}$ and
$\stackrel{{\scriptstyle d}}{{\to}}$, respectively.

######  Theorem 2.2.

As $n\to\infty$, we have $\widehat{\theta}_{n}\stackrel{{\scriptstyle
p}}{{\to}}\theta$ and
$\sqrt{n}(\widehat{\theta}_{n}-\theta)\stackrel{{\scriptstyle
d}}{{\to}}N(0,\sigma^{2}_{\theta})$ with asymptotic variance
$\sigma^{2}_{\theta}=\sigma_{2}^{2}+\sigma_{23}+\tfrac{1}{4}\sigma_{3}^{2}\in(0,\infty)$.
Here, $\sigma_{j}^{2}$ and $\sigma_{jk}$, respectively, denote the asymptotic
variance of $\widehat{F}_{j,n}$ and the asymptotic covariance of
$\widehat{F}_{j,n}$ and $\widehat{F}_{k,n}$, $j,k=1,2,3$.

A more detailed formula for the asymptotic variance $\sigma^{2}_{\theta}$ is
offered in the Supplementary Material. The results of Theorem 2.2 could be
combined with a consistent variance estimator to construct Wald tests and
related confidence intervals. However, such inference procedures typically
exhibit a suboptimal control of the type-I error rate and the confidence
level, respectively, especially for small sample sizes. That is why we propose
a resampling-based approach in Section 3 below.

###  2.3 Discussion of related approaches in the literature

Let us review the present approach in the light of existing approaches and
suggestions from the literature. Seigel and Podgor37 proposed to compare the
counting processes for the competing risks of type 1 and 2 at time $\tau$, say
$N_{1}(\tau)-N_{2}(\tau)$, i.e., McNemar’s statistic. A statistical analysis
of this difference would require taking the censoring rate into account. In
contrast, the relative treatment effect estimator (2) relates to the
difference of cumulative incidence functions through a one-to-one mapping:
$F_{1}(\tau)-F_{2}(\tau)=(1-\theta)-\theta=1-2\theta$. Similarly,
$1-2\widehat{\theta}_{n}=\widehat{F}_{1,n}(\tau)-\widehat{F}_{2,n}(\tau)$. The
latter may be written as
$\int_{0}^{\tau}(N_{1}-N_{2})(du)/\widehat{G}_{n}(u-)$, where the denominator
is the left-continuous version of the Kaplan-Meier estimator for being
uncensored (in the competing risks data set). In that sense, our approach may
be called an inverse-probability-of-censoring-weighting (IPCW) version of the
suggestion by Seigel and Podgor37.

The previous representation also illustrates the difference from the class of
test statistics suggested by Dabrowska9: she proposed $T=\int
K_{u}(s)(\tilde{N}_{1}-\tilde{N}_{2})(ds)+\int
K_{c}(s)(\tilde{N}_{3}-\tilde{N}_{4})(ds)$. Here $\tilde{N}_{1}$ and
$\tilde{N}_{2}$ are defined like $N_{1}$, $N_{2}$, just based on the
completely uncensored pairs, $\tilde{N}_{3}$ and $\tilde{N}_{4}$ correspond to
the singly censored data points, and $K_{u}$ and $K_{c}$ are some scoring
processes.

Let us also compare the present method to the Kaplan-Meier estimator-based
approach in Dobler13. Therein, a different variant of the relative treatment
effect is analyzed, say
$\bar{\theta}=P(T_{11}>T_{22})+\tfrac{1}{2}P(T_{11}=T_{22})$. That is,
outcomes from different treatments of different individuals are compared.
Other names for this parameter are _D-value_ or _Mann-Whitney parameter_18,
11, and also the _C-index_20 is related. It is known 13 that $\bar{\theta}$
can be represented by an intergral that involves only the marginal survival
functions. Consequently, estimation could be based on two marginal and
dependent Kaplan-Meier estimators. In contrast, the method developed in the
present paper is based on direct comparison within each pair. In this sense,
confounding is avoided or at least reduced. The novel competing risks-based
approach will thus make better use of the available information on the
dependence structure within the paired observations. It has been argued18, 19
that it is challenging to draw causal conclusions for $\bar{\theta}$ if not
paired but only sample-specific measurements are available. In addition,
Example 2.3 illustrates that $\bar{\theta}$ might exhibit some undesirable
proporties which do not occur for $\theta$. Finally, the parameter
$\tilde{\theta}$ from (1) is related to the area under the ROC curve (AUC)
which is sometimes written in a similar way.32

We wish to illustrate some differences between the estimands $\theta$ and
$\bar{\theta}$ in addition to examples from the literature18, 19. In
particular, we will point out two cases with different implications for
subpopulations. Let us suppose that the following data set is fully
observable; a variant of the following examples has been kindly provided by
Katharina Kramer (University of Augsburg):

pair | $T_{1}$ | $T_{2}$ | $\textnormal{sign}(T_{1}-T_{2})$ | subgroup
---|---|---|---|---
1 | 2 | 1 | +1 | 1
2 | 4 | 3 | +1 | 1
3 | 6 | 5 | +1 | 2
4 | 8 | 7 | +1 | 2

Here, we also suppose that additional information (a subgroup) is available,
e.g., males and females. For simplicity, let $\tau=\infty$. In the example
above, the above-mentioned estimators of $\theta$ and $\bar{\theta}$ yield the
estimates $\widehat{\theta}_{n}=1$ and
$\widehat{\bar{\theta}}_{n}=\tfrac{10}{16}$, respectively. These are the
estimates for the whole population, i.e., both subgroups combined. We would
like to point out that in this completely observable case the estimators
simplify to the empirical fractions: $\widehat{\theta}_{n}$ is equivalent to
the sign test statistic and $\widehat{\bar{\theta}}_{n}$ is equivalent to the
Mann-Whitney U test statistic. Within the subsamples, however, the estimates
are $\widehat{\theta}_{j,n}=1$ and
$\widehat{\bar{\theta}}_{j,n}=\tfrac{3}{4}$, $j=1,2$. Extending this example
reveals that $\widehat{\bar{\theta}}_{n}$ can get arbitrarily close to
$\tfrac{1}{2}$, i.e., basically no Treatment 1 benefit, whereas
$\widehat{\bar{\theta}}_{j,n}=\tfrac{3}{4}>\tfrac{1}{2}$. This cannot happen
with $\hat{\theta}_{n}$ which is a convex combination of the subgroup-specific
$\hat{\theta}_{1,n}$ and $\hat{\theta}_{2,n}$ (in the fully observable case).

In contrast, the following example is a case where
$\widehat{\bar{\theta}}_{1,n}=\tfrac{3}{4}>\tfrac{1}{2}$ and
$\widehat{\bar{\theta}}_{2,n}=\tfrac{1}{2}$ for the subsamples but
$\widehat{\bar{\theta}}_{n}=\tfrac{7}{16}<\tfrac{1}{2}$ for the whole sample:

pair | $T_{1}$ | $T_{2}$ | $\textnormal{sign}(T_{1}-T_{2})$ | subgroup
---|---|---|---|---
1 | 2 | 1 | +1 | 1
2 | 4 | 3 | +1 | 1
3 | 5 | 6 | -1 | 2
4 | 8 | 7 | +1 | 2

That is, for each subgroup Treatment 1 seems beneficial or at least not
harmful in view of $\bar{\theta}$, whereas is seems harmful for the whole
population. In contrast, $\widehat{\theta}_{n}=\tfrac{3}{4}$ which is in line
with $\widehat{\theta}_{1,n}=1$ and $\widehat{\theta}_{2,n}=\tfrac{1}{2}$.
This illustrates that, in some special cases, drawing conclusions for
subpopulations can be difficult based on $\bar{\theta}$. As we saw, these
phenomenons do not seem to occur for $\theta$. It should be pointed out though
that such comparisons with subgroups are more challenging in the censored case
because there is in general no direct connection between
$\widehat{\theta}_{n}$ and $(\widehat{\theta}_{1,n},\widehat{\theta}_{2,n})$
any more. It might also be useful to complement the estimand $\theta$ with an
average treatment effect, e.g., differences of the restricted mean survival
times (RMST) $E(\min(T_{1},\tau)-\min(T_{2},\tau))$ to get additional insight
in the effectiveness of a treatment with respect to another one.

##  3 Inference on the relative treatment effect

One resampling option is a variant of the classical bootstrap 15, 16, i.e.,
draw $n$ times independently with replacement from the competing risks data
pairs $(Z_{i},\varepsilon_{i}),i=1,\dots,n$, and recompute
$\widehat{\theta}_{n}$ based on the drawn bootstrap sample. By repeating this
procedure a large number $B$ of times, the collection of the normalized
bootstrapped relative treatment effect estimators, say
$W_{n,b}^{*}=\sqrt{n}(\widehat{\theta}^{*}_{n,b}-\widehat{\theta}_{n}),b=1,\dots,B$,
can be used to estimate different aspects of the distribution of
$W_{n}=\sqrt{n}(\widehat{\theta}_{n}-\theta)$, e.g., the $(1-\alpha)$-quantile
which is necessary for a right-tailed test.

Another, perhaps less well-known resampling option, is given by data re-
randomization. The general procedure is similar to the bootstrap, except that
data points are not drawn with replacement but instead some other random
variation is introduced which is related to an algebraic group structure. One
popular example is random permutation of the data points in a two independent
samples setting which leads to permutation tests. Randomization tests can be
shown to be finitely exact if the re-randomization procedure reflects the data
generation process. Recent works revisited the finite exactness of
randomization tests.21, 22 Similarly, finite exactness of confidence intervals
can be shown under randomization-invariance (up to the discreteness of the
randomization distribution). In addition, Dobler13 argued the asymptotic
exactness of randomization tests (with finite exactness in special cases) even
if the data generation process does not exactly match the re-randomization
method.

We will apply one such randomization approach is the present paper. To
motivate it, consider for a moment the strong null hypothesis that both
treatments are completely exchangeable in every respect. In that case, we
would have for the relative treatment effect $\theta=0.5$. However, under the
weak null hypothesis $H_{0}:\theta=0.5$, it is well possible that the
treatments are not exchangeable; for example, if the survival functions
related to $T_{1}$ and $T_{2}$ are allowed to be different. A treatment re-
randomization, i.e., a random re-labeling of the assigned treatment within
each data pair, would lead to an artificial situation in which both treatments
are exchangeable, which in turn implies $H_{0}$. That is, the relative
treatment effect estimator based on the re-randomized data set targets the
value $0.5$. Thus, randomization tests and confidence intervals that are based
on critical values obtained through the described randomization procedure will
be finitely exact under the sharp null hypothesis of treatment
exchangeability. Additionally, it will be asymptotically exact under the weak
null hypothesis $H_{0}$ if the normalized relative treatment effect estimators
are suitably studentized so that the limit distribution is pivotal. Note that
this re-randomization is different from the approach used in classical
permutation tests for two independent samples, as the present data consist of
paired data, and the randomization is done within each pair.

Let us make the randomization approach more explicit. In terms of a fixed
competing risks data set $(Z_{i},\varepsilon_{i}),i=1,\dots,n$, a re-labeling
of both treatments would mean that the times $Z_{i}$ remain unchanged but
every occurrence of a type-1 and type-2 event will be randomly re-labeled a
type-1 or type-2 event, each with probability 50%. Denote a thus obtained
randomized data set by $(\tilde{Z}_{i},\tilde{\varepsilon}_{i}),i=1,\dots,n$,
and the resulting randomized relative treatment effect by
$\tilde{\theta}_{n}$, where the randomization procedure is considered to be
random. For realizations of the re-labeling, while the original data are kept
fixed, we will again use the index $b=1,\dots,B$ in the subscript to indicate
the randomization iteration. It remains to justify the conditional convergence
in distribution of the randomized relative treatment effect estimator. To this
end, we introduce the notation $\mathcal{L}$ for the distribution of a random
variable and let $d$ denote a distance which metrizes the space of
distributions on $\mathbb{R}$, e.g., the Prokhorov distance6 pp. 72-73.
Finally, let $\stackrel{{\scriptstyle P}}{{\to}}$ denote convergence in
probability. Let $\widehat{\sigma}_{\theta,n}^{2}$ and
$\tilde{\sigma}_{\theta,n}^{2}$ be suitable consistent estimators of the
asymptotic variances of $\widehat{\theta}_{n}$ and $\tilde{\theta}_{n}$,
respectively; see Section B in the Supplementary Material for this paper and
Subsection 4.1 for details. Furthermore, we wish to point out that
$P\text{-}\lim_{n\to\infty}\widehat{\sigma}_{\theta,n}^{2}=\sigma^{2}_{\theta}\neq
P\text{-}\lim_{n\to\infty}\tilde{\sigma}_{\theta,n}^{2}=\tilde{\sigma}^{2}_{\theta}$
in general if the two treatments are not exchangeable. This underlines the
necessity to studentize the (randomized) estimator.

######  Theorem 3.1.

As $n\to\infty$,

$d\Big{(}\mathcal{L}\Big{(}\frac{\sqrt{n}(\tilde{\theta}_{n}-0.5)}{\tilde{\sigma}_{\theta,n}}\
\Big{|}\
Z_{1},\varepsilon_{1},Z_{2},\varepsilon_{2},\dots\Big{)},\mathcal{L}\Big{(}\frac{\sqrt{n}(\widehat{\theta}_{n}-\theta)}{\widehat{\sigma}_{\theta,n}}\Big{)}\Big{)}\stackrel{{\scriptstyle
P}}{{\to}}0.$

Based on Theorem 3.1, consistent one-sided hypothesis tests for
$H_{0}:\theta\leq 0.5$ against $H_{a}:\theta>0.5$ of asymptotic level
$\alpha\in(0,1)$ are given by rejecting $H_{0}$ if and only if
$\widehat{\theta}_{n}$ exceeds $0.5+\widehat{\sigma}_{\theta,n}\cdot
c^{*}_{n}(1-\alpha)$. Here, $c^{*}_{n}(1-\alpha)$ denotes the conditional
$(1-\alpha)$-quantile of $(\tilde{\theta}_{n}-0.5)/\tilde{\sigma}_{\theta,n}$
given the data. Such tests control the significance level exactly under the
sharp null hypothesis of the exchangeability of both treatments (and censoring
distributions). Yet, they might still be slightly conservative for very small
sample sizes if the non-randomized version***By a randomized version of a
test, we mean tests which reject the null hypothesis with a probability
$\pi\in(0,1)$ if the test statistic is equal to the critical value. This is
not to be confused with the data re-randomization procedure which is used to
compute the critical or $p$-value. of the randomization tests are used.
Equivalently, $[\widehat{\theta}_{n}+\widehat{\sigma}_{\theta,n}\cdot
c^{*}_{n}(1-\alpha),1]$ are one-sided asymptotic level $(1-\alpha)$ confidence
intervals. Similarly, two-sided tests and confidence intervals can be obtained
if additionally the lower randomization-based quantiles are used. As mentioned
before, these confidence intervals will also be exact under treatment-
exchangeability, up to the discreteness of the randomization distribution.

Another option is to apply a differentiable and one-to-one transformation
$\phi:(0,1)\to\mathbb{R}$ to the relative treatment effect, e.g.,
$\phi(\theta)=\log(-\log(\theta))$. Confidence intervals based on
transformation-based statistics
$\\{\phi(\widehat{\theta}_{n})-\phi(\theta)\\}/\\{\phi^{\prime}(\widehat{\theta}_{n})\cdot\widehat{\sigma}_{\theta,n}\\}$
and the quantiles of their randomization version are ensured to be contained
in $[0,1]$. The delta-method in combination with Slutzky’s theorem justifies
the asymptotic correctness of these adjusted inference procedures.

##  4 Simulation studies

###  4.1 General remarks

The small sample properties are analyzed with the help of simulation studies.
Next to simulations for assessing the size of the proposed right-tailed tests
under the null hypothesis $H_{0}:\theta=0.5$, we also conducted simulations
regarding the power of these tests. All simulations were conducted under R
version 4.1.035. We also used the R package _etm_4 in which the Aalen-Johansen
estimators and all related variance and covariance estimators of Greenwood-
type are implemented.

###  4.2 Size under the null hypothesis

In our simulation study, we only included non-exchangeably distributed data in
order to have a fair comparison between the non-randomized asymptotic,
bootstrap, and randomization tests. We considered the following combinations
of simulation scenarios:

  * •

sample sizes: $n\in\\{25,50,75,100,125,150\\}$;

  * •

significance levels: $\alpha\in\\{1\%,5\%,10\%\\}$; only the results for
$\alpha=5\%$ are included in the main body of this paper, the others in
Section D.1 of the Supplementary Material;

  * •

copulas: Gumbel-Hougaard with parameter equal to 5;
Clayton with parameter equal to -0.6;

  * •

marginal distributions: $Exp(2)$ versus a $50/50$
$Exp(3)$-$Exp(\lambda)$-mixture;
$Gompertz(0.6,b)$ versus $Exp(3)$;

  * •

censoring distributions (the same in both treatment groups): $U(0,a)$ with
$a\in\\{1.1,1.6,2.7\\}$ for the first combination of marginal distributions,
and
$a\in\\{0.7,1,1.75\\}$ for the second combination;

  * •

time end point $\tau=1$ for the first combination of marginal distributions,
$\tau=0.6$ for the second combination;

  * •

5,000 iterations of each tests based on 2,000 bootstrap and randomization
iterations, respectively.

The above-indicated rate parameters $\lambda$ and $b$ of the marginal
distributions were found through numerous generations of large data sets and
they were chosen such that $\theta\approx 0.5$, i.e., the null hypothesis is
considered true. The parameters of the censoring distributions resulted in
censoring rates of $38\%$ to $42\%$ (strong), $27\%$ to $34\%$ (medium), and
$17\%$ to $27\%$ (light); these rates were found through simulations, where
truncations at $\tau$ were also considered censorings.

To elaborate a bit more on the simulation steps, we would like to point out
that the bivariate copula data are generated first. Next, the quantile
functions of the marginal distributions are used to transform the copula data
into bivariate data with the pre-selected dependence structure and marginal
distributions. Finally, censoring is introduced by taking the minimum of the
event and the simulated censoring times. For the purpose of estimating the
relative treatment effect, the synthetic data are next transformed into a
competing risks data set as described in Section 2.2. It should be pointed out
that the cumulative incidence functions underlying the transformed competing
risks data set neither have to be specified, nor do they play an important
role, except for at time $\tau$.

The results displayed in Table 1 can be summarized as follows. Among the tests
based on the untransformed relative treatment effect (“lin.” in the table),
the asymptotical test tends to be a bit liberal, while the bootstrap test is
somewhat conservative. This is more pronounced for the smaller sample sizes up
to $n=100$. The randomization test is generally closest to the selected
significance level. For larger sample sizes ($n\in\\{125,150\\}$), the sizes
of all tests approach the $5\%$ level quite accurately.

Table 1 about here.

Table 1: Simulated sizes of the right-tailed tests (in %) with nominal significance level $\alpha=5\%$. Abbreviations: Copula: GH = Gumbel-Hougaard; critical values: asy. = asymptotical normal, bs. = bootstrap, rand. = randomization; tests: lin. = linear, tra. = $\log$-$\log$-transformed. |  |  | light censoring | medium censoring | strong censoring
---|---|---|---|---|---
|  |  | asy. | bs. | rand. | asy. | bs. | rand. | asy. | bs. | rand.
copula | distribution | $n$ | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra.
GH | Exp vs. | $25$ | $7.0$ | $4.4$ | $3.7$ | $5.8$ | $5.6$ | $5.5$ | $7.1$ | $4.8$ | $4.3$ | $6.5$ | $5.7$ | $5.6$ | $6.4$ | $4.5$ | $3.6$ | $5.1$ | $6.7$ | $6.8$
| Exp mix | $50$ | $6.0$ | $4.4$ | $4.3$ | $5.2$ | $5.2$ | $5.2$ | $6.2$ | $4.7$ | $4.4$ | $5.4$ | $5.4$ | $5.4$ | $5.2$ | $4.2$ | $3.3$ | $3.9$ | $6.2$ | $6.2$
|  | $75$ | $5.6$ | $4.3$ | $4.3$ | $5.0$ | $5.0$ | $5.0$ | $6.7$ | $5.2$ | $4.9$ | $5.8$ | $5.9$ | $5.9$ | $4.5$ | $3.8$ | $3.0$ | $3.5$ | $5.9$ | $5.8$
|  | $100$ | $5.5$ | $4.6$ | $4.7$ | $5.2$ | $5.2$ | $5.2$ | $5.9$ | $4.7$ | $4.6$ | $5.2$ | $5.1$ | $5.1$ | $4.9$ | $4.2$ | $3.3$ | $3.8$ | $6.8$ | $6.8$
|  | $125$ | $5.1$ | $4.3$ | $4.5$ | $4.9$ | $4.9$ | $4.9$ | $5.4$ | $4.5$ | $4.5$ | $4.8$ | $5.1$ | $5.1$ | $4.7$ | $4.3$ | $3.5$ | $3.7$ | $7.2$ | $7.2$
|  | $150$ | $5.8$ | $5.1$ | $5.3$ | $5.5$ | $5.5$ | $5.5$ | $5.7$ | $4.9$ | $5.0$ | $5.4$ | $5.3$ | $5.3$ | $4.5$ | $3.9$ | $3.5$ | $3.7$ | $6.7$ | $6.7$
GH | Gompertz | $25$ | $8.1$ | $5.0$ | $2.2$ | $4.5$ | $6.2$ | $6.2$ | $9.7$ | $7.0$ | $3.0$ | $5.8$ | $9.3$ | $9.5$ | $11.5$ | $8.8$ | $4.2$ | $5.1$ | $15.9$ | $15.5$
| vs. Exp | $50$ | $6.0$ | $4.3$ | $3.7$ | $4.6$ | $5.2$ | $5.2$ | $6.2$ | $4.5$ | $2.6$ | $3.6$ | $5.9$ | $5.9$ | $8.4$ | $6.6$ | $2.8$ | $4.4$ | $12.7$ | $12.7$
|  | $75$ | $5.0$ | $3.7$ | $3.6$ | $4.4$ | $4.2$ | $4.2$ | $5.3$ | $4.0$ | $3.0$ | $3.8$ | $5.0$ | $5.1$ | $7.1$ | $5.8$ | $2.8$ | $3.9$ | $12.2$ | $12.2$
|  | $100$ | $5.3$ | $4.3$ | $4.5$ | $5.0$ | $5.0$ | $5.0$ | $5.8$ | $4.5$ | $4.0$ | $4.6$ | $5.4$ | $5.4$ | $5.7$ | $4.4$ | $2.3$ | $2.8$ | $10.1$ | $10.1$
|  | $125$ | $5.4$ | $4.5$ | $4.6$ | $5.1$ | $5.2$ | $5.2$ | $5.5$ | $4.6$ | $4.6$ | $5.0$ | $5.2$ | $5.2$ | $5.0$ | $4.1$ | $2.4$ | $2.9$ | $9.8$ | $9.8$
|  | $150$ | $5.1$ | $4.4$ | $4.6$ | $4.9$ | $4.8$ | $4.8$ | $5.3$ | $4.3$ | $4.5$ | $4.9$ | $4.8$ | $4.8$ | $4.4$ | $3.9$ | $2.6$ | $2.9$ | $9.1$ | $9.1$
Clayton | Exp vs. | $25$ | $6.2$ | $4.2$ | $3.2$ | $5.2$ | $5.0$ | $5.0$ | $6.4$ | $3.9$ | $3.1$ | $5.4$ | $5.2$ | $5.2$ | $6.4$ | $4.0$ | $3.1$ | $5.6$ | $4.7$ | $4.7$
| Exp mix | $50$ | $5.5$ | $4.0$ | $3.8$ | $4.7$ | $4.8$ | $4.8$ | $5.9$ | $4.3$ | $3.7$ | $4.7$ | $4.9$ | $4.9$ | $5.6$ | $4.0$ | $3.4$ | $4.6$ | $5.4$ | $5.4$
|  | $75$ | $5.6$ | $4.5$ | $4.3$ | $4.9$ | $5.0$ | $5.0$ | $5.3$ | $4.2$ | $3.9$ | $4.5$ | $4.8$ | $4.8$ | $5.3$ | $4.4$ | $4.0$ | $4.6$ | $5.4$ | $5.4$
|  | $100$ | $5.5$ | $4.6$ | $4.6$ | $4.9$ | $5.4$ | $5.4$ | $5.3$ | $4.3$ | $4.2$ | $4.6$ | $4.8$ | $4.8$ | $5.0$ | $4.1$ | $3.7$ | $4.1$ | $5.0$ | $4.9$
|  | $125$ | $4.7$ | $3.9$ | $4.0$ | $4.3$ | $4.7$ | $4.7$ | $5.3$ | $4.7$ | $4.7$ | $5.0$ | $5.3$ | $5.3$ | $5.2$ | $4.4$ | $3.8$ | $4.2$ | $5.3$ | $5.3$
|  | $150$ | $4.9$ | $4.1$ | $4.3$ | $4.5$ | $4.6$ | $4.6$ | $4.9$ | $4.1$ | $4.3$ | $4.5$ | $4.9$ | $4.9$ | $5.1$ | $4.6$ | $4.0$ | $4.2$ | $5.4$ | $5.4$
Clayton | Gompertz | $25$ | $6.8$ | $4.5$ | $3.0$ | $5.5$ | $5.4$ | $5.4$ | $7.3$ | $4.5$ | $3.3$ | $5.9$ | $5.3$ | $5.3$ | $8.3$ | $5.1$ | $3.9$ | $6.7$ | $5.8$ | $5.8$
| vs. Exp | $50$ | $5.5$ | $4.0$ | $3.6$ | $4.6$ | $5.0$ | $5.0$ | $6.7$ | $5.0$ | $4.1$ | $5.2$ | $5.8$ | $5.8$ | $6.7$ | $4.9$ | $3.9$ | $5.2$ | $5.6$ | $5.6$
|  | $75$ | $5.4$ | $4.2$ | $4.0$ | $4.8$ | $4.9$ | $4.9$ | $6.0$ | $4.3$ | $4.0$ | $4.9$ | $5.4$ | $5.4$ | $6.0$ | $4.6$ | $3.8$ | $4.7$ | $5.5$ | $5.5$
|  | $100$ | $5.8$ | $4.8$ | $4.8$ | $5.2$ | $5.4$ | $5.4$ | $5.7$ | $4.7$ | $4.3$ | $4.7$ | $5.3$ | $5.3$ | $6.3$ | $5.0$ | $4.2$ | $4.8$ | $5.8$ | $5.8$
|  | $125$ | $4.8$ | $4.0$ | $4.0$ | $4.4$ | $4.6$ | $4.6$ | $5.1$ | $4.5$ | $4.2$ | $4.6$ | $4.9$ | $4.9$ | $5.9$ | $4.9$ | $4.3$ | $4.8$ | $5.8$ | $5.8$
|  | $150$ | $5.4$ | $4.7$ | $4.7$ | $5.1$ | $5.2$ | $5.2$ | $5.9$ | $4.8$ | $4.7$ | $5.0$ | $5.5$ | $5.5$ | $5.7$ | $4.9$ | $4.4$ | $4.8$ | $5.6$ | $5.6$

One notable peculiarity is the combination of strong censoring, the Gumbel-
Hougaard copula, and, in particular, Gompertz versus exponential marginals:
here, the bootstrap test stays very conservative, and the randomization test
is very anti-conservative. There is a slight improvement when the sample size
increases. In this scenario, the asymptotic test is also anti-conservative but
it surprisingly performs better than the randomization test. For smaller
sample sizes ($n\in\\{25,50\\}$), similar observations about the simulation
results can be made for the related medium censoring case. Figure 1
illustrates that the above-discussed scenario is indeed a very challenging
one: many data points are converted into censorings on the competing risks
scale; the censoring rate amounts to about 52%. At the same time, there is a
very strong correlation between the survival times of a pair, and two very
different marginal distributions. For smaller sample sizes, most of these
characteristics are hardly visible. That is why this is by far the most
challenging simulation setting.

Figure 1: Scatterplots of $n=200$ simulated data points according to a Gumbel-
Hougaard copula and, respectively, Gumbel and exponentially distributed
marginals. The red line is the diagonal $y=x$. The dashed lines illustrate the
upper end of the censoring support (0.7) and $\tau=0.6$.
Upper panels: complete data set (no censoring); lower panels: censored data
set.
Left panels: raw data; censorings in both coordinates are denoted by "+";
censorings only in the horizontal coordinate are denoted "X"; censorings only
in the vertical coordinate are denoted "Y"; completely uncensored data points
are denoted by red circles.
Right panels: competing risks data after transformation. The event time is the
minimum of both coordinates. The symbols "1", "2", and "C" represent whether
the data point corresponds to an observed event of type 1, type 2, or to a
censoring.

The $\log$-$\log$-transformation (“tra.” in the table) seems to rectify the
liberality of the asymptotic test and the conservativeness of the bootstrap
test in most scenarios. However, the transformation has nearly no effect on
the randomization test.

###  4.3 Power simulations

In addition to the size simulations under $H_{0}:\theta=0.5$, we have also
conducted a simulation study to assess the power of the one-tailed version of
the developed test. Since the results for the transformed and untransformed
test statistics are very much alike, we have solely focused on the latter. We
also considered two competitor tests, also in their one-tailed versions: the
paired Prentice-Wilcoxon test31 which was found to be very powerful in the
comparative simulation study by 42; the stratified log-rank test30 which was
more closely analyzed under correlated frailty models. In order to ensure a
fair comparison, the randomization versions of all tests were used, such that
all of them control the significance level for finite sample sizes under
exchangeability. In this subsection, we chose the significance level
$\alpha=5\%$.

We considered the same two copulas as in the previous subsection, the sample
sizes $n\in\\{25,50,100\\}$, 1,000 test replications, 1,000 randomization
iterations, and the following three marginal distribution scenarios:

  1. 1.

Mixture of the $Exp(2)$-exponential and the $U(0,2)$-uniform distribution
against the $Exp(2)$-exponential distribution; the censoring times were
independently
$U(0,2.5)$-distributed and $\tau=1.9$. This scenario departs from the sharp
null hypothesis of exchangeability into alternatives with crossing hazard
rates at late time points close to $\tau$ as the mixing parameter puts more
and more weight on the uniform distribution.

  2. 2.

Mixture of the $Exp(2)$-exponential and the Gompertz distribution with shape
parameter $0.1$ and rate parameter $2$ against the $Exp(2)$-exponential
distribution; the censoring times were independently $U(0,2.5)$-distributed
and $\tau=1.8$. This scenario departs from the sharp null hypothesis of
exchangeability into alternatives with crossing hazard rates at central time
points as the mixing parameter puts more and more weight on the Gompertz
distribution.

  3. 3.

The $Exp(2/k)$-exponential distribution against the $Exp(2)$-exponential
distribution; the censoring times were independently $U(0,2)$-distributed and
$\tau=1.3$. This scenario departs from the sharp null hypothesis of
exchangeability into alternatives with parallel hazard rates as the scale
parameter $k$ increases from $1$ to $2$.

Figure 2: Simulated power of three selected right-tailed randomization-based
tests with significance level $\alpha=5\%$.

For now, we only focus on the results for Scenario 1 graphically presented in
Figure 2; the results for Scenarios 2 and 3 are presented in Section D.2 of
the Supplementary Material and they are similar to those presented here. We
can see from Figure 2 by comparing each combination of left and right panel
that the copula that connects the lifetimes apparently has only little
influence on the performance of the tests. Not surprisingly, the power of all
tests increases when the sample size increases (top to bottom in the figure)
and when we depart from the null hypothesis (from left to right within each
panel). The paired Prentice-Wilcoxon test31 is always the most powerful one.
In most cases, the proposed test has the next higher power but its performance
is generally very similar to that of the stratified log-rank test30.

Multiple comments are in order. First, in this simulation study we could
confirm the earlier findings42 that the paired Prentice-Wilcoxon test is
indeed quite powerful. This is also not surprising because it is based on an
efficient score approach. Next, the power of the stratified log-rank test
could possibly be greater if the optimal combination with the unstratified
log-rank test was used.30 However, an implementation of the combination would
be beyond the scope of the present paper. In addition, it should not be
forgotten that the considered competitor tests were proposed for the sharp
null hypothesis of equal survival distributions; it is only natural that such
tests potentially have a greater power than the proposed test which was
designed for the weak null hypothesis $H_{0}:\theta=0.5$. Another reason for
the relatively high power of the log-rank test and the paired Prentice-
Wilcoxon test is the fact that the generated lifetimes were not truncated at
$\tau$, whereas smaller values of $\tau$ let the relative treatment effect get
closer to $0.5$, i.e., closer to the null hypothesis. In the light of these
facts, the power of the proposed test procedure is very competitive. In
addition, as initially mentioned, our focus was on the development of an
easily interpretable estimand-based inference procedure, rather than
developing the most powerful test for comparing two survival functions.

##  5 Data example

We illustrate our methodology by re-analyzing a well known benchmark data set
which has been published in the R package survival38. The data set diabetic,
in detail described and analyzed by Huster et al.24, contains $394$
observations from a trial including $197$ patients with “high-risk" diabetic
retinopathy, a complication associated with diabetes mellitus that frequently
leads to blindness. In this trial each patient acts as its own control: one
eye was randomized to a laser photocoagulation, while the other eye received
no treatment.

Apart from suffering from diabetic retinopathy, the inclusion criterion of the
trial was a visual acuity of at least 20/100 in both eyes. The aim of the
study was to investigate the effect of the laser treatment on delaying the
onset of blindness, defined by a visual acuity of less than 5/200 at two
consecutive visits after four months. Thus, the survival times are the times
until blindness for the eyes. Censoring was caused by death, dropout, or the
end of the study. Consequently, we consider the censored paired survival
outcomes and their relative treatment effect.

The data set consists of two subgroups, defined by the type of diabetes, i.e.,
patients with juvenile onset diabetes (diagnosis before an age of 20, 114
patients) and adult onset diabetes (83 patients). In these subgroups,
respectively $78.9\%$ and $85.1\%$ of the patients had at least one censoring.
Further, there are several other covariates, e.g., the laser type and a risk
score.

For a better overview of the data, Figure 3(a) visualizes the time until
blindness for the total sample, regardless of the age at diagnosis. From the
pattern of the observations displayed in the figure, we conclude that the
laser treatment seems in general to delay blindness. For a separate
investigation of the juvenile and adult sample, respectively, Figure 3(b)
displays the estimated Kaplan-Meier curves for each eye, i.e. the one
achieving a laser treatment and the other one acting as the control. For both
subsamples, we observe a visible difference between the Kaplan-Meier curves
and again we conclude that the laser treatment seems to delay the onset of
blindness compared to the control, which becomes even more visible in the
adult sample.

Figure 3 about here.

(a) Time until blindness for the diabetic retinopathy data. The values on the
diagonal (denoted by "+") correspond to patients where no blindness occurred
throughout the observational period. Red circles indicate blindness of both
eyes, values below the diagonal (denoted by "Y") indicate blindness of the
control eye only, whereas values above the diagonal (denoted by "X") indicate
blindness of the treated eye only.

(b) Kaplan-Meier-curves for the two eyes (treated and control) for the
diabetic retinopathy data, fitted separately for the juvenile and the adult
sample.

Figure 3: Graphical summaries of the data set.

Table 2 about here.

Table 2: Comparison of $95\%$\- confidence intervals and p-values for the diabetic retinopathy data. Abbreviations: asy. = asymptotical normal, bs. = bootstrap, rand. = randomization; lin. = linear, tra. = $\log$-$\log$-transformed Method | Juvenile sample | Adult sample
---|---|---
$95\%$\- confidence interval | p-value |  $95\%$\- confidence interval | p-value
asy. | lin. |  [0.517, 0.678] | 0.017 |  [0.655, 0.807] | $<0.001$
tra. |  [0.513, 0.673] | 0.025 |  [0.646, 0.798] | $<0.001$
bs. | lin. |  [0.514, 0.680] | 0.014 |  [0.652, 0.802] | $<0.001$
tra. |  [0.517, 0.677] | 0.012 |  [0.655, 0.800] | $<0.001$
rand. | lin. |  [0.515, 0.680] | 0.025 |  [0.654, 0.809] | $<0.001$
tra. |  [0.515, 0.676] | 0.025 |  [0.651, 0.801] | $<0.001$

In order to confirm our visual findings, we will estimate the relative
treatment effect of the laser photocoagulation and the corresponding
confidence intervals, as well as performing the corresponding two-sided
hypothesis tests in order to assess the effectiveness of this therapy. We fix
the maximum follow-up time as $\tau=60$ (indicated by the dashed box in Figure
3(a)) and perform the analysis for both samples, that is juvenile and adult
onset diabetes, separately. For the juvenile sample, we obtain a relative
treatment effect of $\hat{\theta}_{J}=0.598$, for the adult sample we have
$\hat{\theta}_{A}=0.731$. The corresponding $95\%-$ confidence intervals and
the results of the (two-sided) hypothesis test for the different approaches
described in Section 3 are summarized in Table 2. The transformation used for
the analysis is given by $\phi(\theta)=\log(-\log(\theta))$ and all results
based on bootstrap were achieved by using $B=$2,000 bootstrap repetitions. For
each subsample, all confidence intervals and also all p-values are very
similar. In general, the effect of an additional transformation of the test
statistic is rather small and yields very similar confidence intervals and
test results. As already indicated by Figure 3(b), there is a notable
difference of treatment and control eye, which is even larger in the adult
sample. This is confirmed by p-values below $0.001$ for the adult sample for
all tests under consideration. For the juvenile sample, p-values lie between
$0.012$ and $0.025$. Hence, we conclude a significant treatment effect for
both samples at the significance level $\alpha=5\%$.

Our results are similar to the findings from Oakes and Feng30, who
investigated the same data set regarding the treatment effect. In their paper,
they propose three different test approaches and conclude that, for the
example at hand, all resulting treatment effects are significant.

##  6 Discussion

In this paper, we developed a new estimand in the context of paired, right-
censored survival data, the so-called relative treatment effect, to compare
the effectiveness of two treatments. Such data occur for instance in matched
pairs studies. We derived confidence intervals and hypothesis tests and could
demonstrate all desirable propoerties by means of a simulation study.

The relative treatment effect $\theta$ quantifies the stochastic ordering of
treatment outcomes but not how much bigger one survival times is than the
other. In this sense, it is a global measure for the superiority of the first
treatment, although other measures, e.g., about the actual size of the
differences, might also be of major importance in some applications. One
extension of the present method could be the incorporation of additional
patient covariates, e.g., those that are used for the matching of individuals
into pairs, which could be used to tackle the classification problem of who
specifically should receive which treatment. This could be achieved in terms
of semiparametric regression models or by involving relative treatment effects
in a machine learning algorithm.

There are multiple possibilities for other extensions of the present approach.
One open question is how to incorporate additional patients that could not be
matched with others or if multiple patients of one treatment group could be
matched with just one patient of the other group. The latter problem could
potentially be approached by means of an appropriate re-weighting of the
within-pair comparisons. However, there is the risk that the celebrated easy
interpretability of the relative treatment effect could be lost. Another way
to extend the present approach is the incorporation of additional covariates.
Due to the favourable competing risks approach, such an extension could be
achieved rather straightforwardly, e.g., by means of cause-specific hazard
models or subdistribution hazard models.

Finally, we note that the hypothesis test presented in this paper investigates
the significance of the treatment effect. However, there might also occur
situations where one is rather interested in testing whether the deviation of
the treatment effect of $0.5$ is not larger than pre-specified values
$\epsilon_{1}$ and $\epsilon_{2}$, respectively. In other words, this requires
an equivalence test for
$H_{0}^{eq}:\theta\notin(0.5-\varepsilon_{1},0.5+\varepsilon_{2})$ against
$H_{a}^{eq}:\theta\in(0.5-\varepsilon_{1},0.5+\varepsilon_{2})$ (see, for
example, Wellek40). Such an approach could provide a very flexible framework
for statistical inference, address numerous other research questions, and
consequently provide a useful addition to the test proposed in this paper. We
leave the development of such a procedure for future research.

## Acknowledgements

The authors wish to thank Arthur Allignol and Mark Clements for providing a
suitable update of the _etm_ R package, Marialuisa Restaino, Nan van Geloven,
Liesbeth de Wreede, Hein Putter, Giuliana Cortese, Thomas Scheike, Katharina
Kramer, and Paul Blanche for useful discussions.

## Conflict of interest

The authors declare that they have no conflict of interest.

## Data Availability Statement

The R code used for the simulation studies, for generating synthetical data
sets, and for the real data application is available at
https://github.com/dennis-dobler/relative_treatment_effect_paired_survival.

## References

  * Aalen and Johansen 1978 Aalen, O. O. and Johansen, S. (1978).  An empirical transition matrix for non-homogeneous Markov chains based on censored observations.  Scand. J. Statist. 5, 141–150.
  * Aalen 1978 Aalen, O. O. (1978).  Nonparametric inference for a family of counting processes.  Ann. Stat. 6, 701–726.
  * Albers 1988 Albers, W. (1988).  Combined rank tests for randomly censored paired data.  J. Am. Statist. Assoc. 83, 1159–1162.
  * Allignol et al. 2011 Allignol, A., Schumacher, M., and Beyersmann, J. (2011).  Empirical transition matrix of multi-state models: the etm package.  J. Statist. Softw. 38, 1–15.
  * Andersen et al. 1993 Andersen, P. K., Borgan, Ø., Gill, R. D., and Keiding, N. (1993).  Statistical models based on counting processes.  Springer, New York.
  * Billingsley 1999 Billingsley, P. (1999).  Convergence of probability measures.  John Wiley & Sons Inc., New York.
  * Brunner and Munzel 2000 Brunner, E. and Munzel, U. (2000).  The nonparametric Behrens-Fisher problem: asymptotic theory and a small-sample approximation.  Biometrical J. 42, 17–25.
  * Brunner, Vandemeulebroecke, and Mütze 2021 Brunner, E., Vandemeulebroecke, M., and Mütze, T. (2021)  Win odds: An adaptation of the win ratio to include ties.  Statist. Med. 40, 3367–3384.
  * Dabrowska 1990 Dabrowska, D. (1990).  Signed-rank tests for censored matched pairs.  J. Am. Statist. Assoc. 85, 478–485.
  * Dai et al. 2016 Dai, H., Restaino, M., and Wang, H. (2016).  A class of nonparametric bivariate survival function estimators for randomly censored and truncated data.  J. Nonparametr. Statist. 28, 736–751.
  * Dobler and Pauly 2018 Dobler, D. and Pauly, M. (2018).  Bootstrap-and permutation-based inference for the Mann–Whitney effect for right-censored and tied data.  Test 27, 639–658.
  * Dobler 2017 Dobler, D. (2017).  A discontinuity adjustment for subdistribution function confidence bands applied to right-censored competing risks data.  Electron. J. Statist. 11, 3673–3702.
  * Dobler 2023 Dobler, D. (2023).  Randomized Empirical Processes by Algebraic Groups, and Tests for Weak Null Hypotheses.  Bernoulli 29, 1109–1136.
  * Efron 1967 Efron, B. (1967).  The two sample problem with censored data.  In Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, volume 4, pages 831–853.
  * Efron 1979 Efron, B. (1979).  Bootstrap methods: another look at the jackknife.  Ann. Statist. 7, 1–26.
  * Efron 1981 Efron, B. (1981).  Censored data and the bootstrap.  J. Am. Statist. Assoc. 76, 312–319.
  * European Medicines Agency 2017 European Medicines Agency (2017).  Addendum on Estimands and Sensitivity Analysis in Clinical Trials to the Guideline on Statistical Principles for Clinical Trials.
  * Fay et al. 2018 Fay, M. P., Brittain, E. H., Shih, J. H., Follmann, D. A., and Gabriel, E. E. (2018).  Causal estimands and confidence intervals associated with Wilcoxon-Mann-Whitney tests in randomized experiments.  Statist. Med. 37, 2923–2937.
  * Greenland et al. 2020 Greenland, S., Fay, M. P., Brittain, E. H., Shih, J. H., Follmann, D. A., Gabriel, E. E., and Robins, J. M. (2020).  On causal inferences for personalized medicine: how hidden causal assumptions led to erroneous causal claims about the D-value.  Am. Statist. 74, 243–248.
  * Hartman et al. 2023 Hartman, N., Kim, S., He, K., and Kalbfleisch, J. D. (2023).  Pitfalls of the concordance index for survival outcomes.  Statist. Med. 42, 2179–2190.
  * Hemerik and Goeman 2018 Hemerik, J. and Goeman, J. J. (2018).  Exact testing with random permutations.  Test 27, 811–825.
  * Hemerik and Goeman 2021 Hemerik, J. and Goeman, J. J. (2021).  Another look at the lady tasting tea and differences between permutation tests and randomisation tests.  Int. Statist. Rev. 89, 367–381.
  * Holt and Prentice 1974 Holt, J. and Prentice, R. (1974).  Survival analyses in twin studies and matched pair experiments.  Biometrika 61, 17–30.
  * Huster et al. 1989 Huster, W. J., Brookmeyer, R., and Self, S. G. (1989).  Modelling paired survival data with covariates.  Biometrics 45, 145–156.
  * Jeong 2002 Jeong, G.-J. (2002).  A covariate-adjusted logrank test for paired survival data.  Comm. Statist. Appl. Meth. 9, 533–542.
  * Jung 1999 Jung, S.-H. (1999).  Rank tests for matched survival data.  Lifetime Data Anal. 5, 67–79.
  * Konietschke and Pauly 2012 Konietschke, F. and Pauly, M. (2012).  A studentized permutation test for the non-parametric Behrens-Fisher problem in paired data.  Electron. J. Statist. 6, 1358–1372.
  * Le-Rademacher and Brazauskas 2014 Le-Rademacher, J. and Brazauskas, R. (2014).  Inference for paired survival data.  Handbook of survival analysis, 615–632.  Chapman and Hall/CRC.
  * Munzel and Brunner 2002 Munzel, U. and Brunner, E. (2002).  An exact paired rank test.  Biometrical J. 44, 584–593.
  * Oakes and Feng 2010 Oakes, D. and Feng, C. (2010).  Combining stratified and unstratified log-rank tests in paired survival data.  Statist. Med. 29, 1735–1745.
  * O’Brien and Fleming 1987 O’Brien, P. and Fleming, T. (1987).  A paired Prentice-Wilcoxon test for censored paired data.  Biometrics 43, 169–180.
  * Pauly, Asendorf, and Konietschke 2016 Pauly, M., Asendorf, T., and Konietschke, F. (2016).  Permutation-based inference for the AUC: A unified approach for continuous and discontinuous data.  Biometrical J. 58, 1319–1337.
  * Pocock et al. 2012 Pocock, S., Ariti, C., Collier, T., and Wang, D. (2012).  The win ratio: a new approach to the analysis of composite endpoints in clinical trials based on clinical priorities.  Eur. Heart J. 33, 176–182.
  * Pruitt 1993 Pruitt, R. (1993).  Small sample comparison of six bivariate survival curve estimators.  J. Statist. Comput. Sim. 45, 147–167.
  * R Core Team 2021 R Core Team (2021).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Scheike, Holst, and Hjelmborg 2014 Scheike, T. H., Holst, K. K., and Hjelmborg, J. B. (2014).  Estimating twin concordance for bivariate competing risks twin data.  Statist. Med. 33, 1193–1204.
  * Seigel and Podgor 1982 Seigel, M. and Podgor, M. (1982).  A sign test for significance of differences in survivorship curves from paired truncated data.  Contr. Clin. Trials 3, 69–71.
  * Therneau 2020 Therneau, T. (2020).  A package for survival analysis in rR package, available at https://CRAN.R-project.org/package=survival.
  * van der Vaart and Wellner 1996 van der Vaart, A. W. and Wellner, J. A. (1996).  Weak Convergence and Empirical Processes. With Applications to Statistics. Springer, New York.
  * Wellek 2002 Wellek, S. (2002).  Testing statistical hypotheses of equivalence.  Chapman and Hall/CRC.
  * Woolson and Lachenbruch 1980 Woolson, R. and Lachenbruch, P. (1980).  Rank tests for censored matched pairs.  Biometrika 67, 597–606.
  * Woolson and O’Gorman 1992 Woolson, R. and O’Gorman, T. (1992).  A comparison of several tests for censored paired data.  Statist. Med. 11, 193–208.

##  Appendix A Technical details for the competing risks-based estimation of
the relative treatment effect

In this section, we will explain in detail how the transformation of the
paired, right-censored data to competing risks data can be achieved and why
these competing risks data are usable for valid estimation of the relative
treatment effect.

Let us thus assume that a data set consists of independently and identically
distributed data points $(X_{i1},\delta_{i1},X_{i2},\delta_{i2}),i=1,\dots,n$,
as described in Section 2 in the main manuscript. As explained there, we wish
to estimate the relative treatment effect $\theta$ with the help of a
competing risks approach. Our key strategy is to transform the data into a
competing risks data set:

$(Z_{i},\varepsilon_{i})=(\min(\check{T}_{i},\check{C}_{i}),\check{\varepsilon}_{i}\cdot
1\\{\check{T}_{i}\leq\check{C}_{i}\\}),\quad i=1,\dots,n,$

where $\check{T}_{i}=\min(T_{i1},T_{i2},\tau)$ with survival function
$S(t)=P(\check{T}_{i}>t)\ (0\leq t\leq\tau)$ denote the event times,
$\check{C}_{i}=\min(C_{i1},C_{i2})$ with survival function
$G(t)=P(\check{C}_{i}>t)\ (0\leq t\leq\tau)$ denote the censoring times, and
$\check{\varepsilon}_{i}\in\\{1,2,3\\}$ are the event indicators,
respectively. Again, $\check{T}_{i},\check{C}_{i}$, and
$\check{\varepsilon}_{i}$ are not fully observable but
$(Z_{i},\varepsilon_{i})$ is. The three different (artificial) events are
defined as follows:

  * 1.

If the first treated subject experienced an _observable_ event before the
second, i.e., $X_{i1}<X_{i2},\delta_{i1}=1$ or if
$X_{i1}=X_{i2},(\delta_{i1},\delta_{i2})=(1,0)$, we say the first type of
event ($\varepsilon_{i}=1$) took place at time $X_{i1}=T_{i1}$. Note though
that the second just mentioned case is impossible due to the independent and
continuous censoring assumptions:

$P(X_{i1}=X_{i2},\delta_{i1}=1,\delta_{i2}=0)\leq P(T_{i1}=C_{i2})=\int
P(C_{i2}=t)dP^{T_{i1}}(t)=0.$

  * 2.

Similarly, the second type of event ($\varepsilon_{i}=2$) took place at time
$X_{i2}=T_{i2}$ if $X_{i2}<X_{i1}$ and $\delta_{i2}=1$.

  * 3.

If both events were _observed_ simultaneously, i.e., $X_{i1}=X_{i2}$ and
$\delta_{i1}=\delta_{i2}=1$, the third type of event ($\varepsilon_{i}=3$)
took place at $X_{i1}=X_{i2}=T_{i1}=T_{i2}$.

In all other cases, the observation is censored ($\varepsilon_{i}=0$) at

$\min(X_{i1},X_{i2})=\min(C_{i1},C_{i2}).$

In the following, we suppress the index $i$ for ease of presentation when
there is no need to specify the precise subject. We denote the cause-specific
cumulative hazard functions in the competing risks framework by
$A_{j},j=1,2,3$. The following lemma is the crucial step towards estimating
the relative treatment effect $\theta$; in fact, it shows that no relevant
information is lost by the above-described conversion to a competing risks
data set.

######  Lemma A.1.

The cause-specific Nelson-Aalen estimators

$\widehat{A}_{j,n}(t)=\sum_{u\leq
t}\tfrac{\sum_{i=1}^{n}1\\{Z_{i}=u,\varepsilon_{i}=j\\}}{\sum_{i=1}^{n}1\\{Z_{i}\geq
u\\}}$

are uniformly consistent for $A_{j}(t),j=1,2,3$ in $t\in[0,\tau]$ as
$n\to\infty$. In addition, $\\{n^{1/2}(\widehat{A}_{j,n}-A_{j})\\}_{j=1}^{3}$
converges in distribution to a three-dimensional zero-mean Gaussian process as
$n\to\infty$.

The proof of Lemma A.1 (provided in Section C below) reveals that the
transformation of the paired survival data into competing risks data preserves
the underlying intensities. As a consequence, other properties of the Nelson-
Aalen estimators are retrieved, such as their interpretation as nonparametric
maximum likelihood estimators. Another consequence of Lemma A.1 is that the
Aalen-Johansen estimators1 for the cumulative incidence functions for all
event times are similarly estimable; cf. Dobler12 for the general case with
both continuous and discrete components in the event time distribution. To be
specific, estimation of the relative treatment effect is achievable as
follows:

$\widehat{\theta}_{n}=\widehat{F}_{2,n}(\tau)+\tfrac{1}{2}\widehat{F}_{3,n}(\tau)=\int_{0}^{\tau}\widehat{S}_{n}(u-)d(\widehat{A}_{2,n}+\tfrac{1}{2}\widehat{A}_{3,n})(u),$

where $\widehat{S}_{n}(t)=\prod_{u\leq
t}\\{1-d(\widehat{A}_{1,n}+\widehat{A}_{2,n}+\widehat{A}_{3,n})(u)\\}$ denotes
the Kaplan-Meier estimator of $S(t)=P(\check{T}>t)$ and $\widehat{F}_{j,n}$
are the Aalen-Johansen estimators of $F_{j}$, $j=2,3$. The minus sign in an
argument indicates the left-continuous version of a function.

##  Appendix B Asymptotic variances of the relative treatment effect
estimators, and consistent variance estimators

###  B.1 Asymptotic variances of $\sqrt{n}(\widehat{\theta}_{n}-\theta)$ and
$\sqrt{n}(\tilde{\theta}_{n}-\tfrac{1}{2})$

Tedious but straightforward calculations revealed the following expression for
the asymptotic variance of the relative treatment effect estimator,
$\sqrt{n}(\widehat{\theta}_{n}-\theta)$:

$\displaystyle\sigma^{2}_{\theta}$
$\displaystyle=\int_{0}^{\tau}\int_{0}^{\tau}S(u-)S(v-)\Big{[}\int_{0}^{\min(u,v)-}\frac{\sigma^{2}_{\bullet}(dw)}{1-\Delta
A_{\bullet}(w)}(A_{2}+\tfrac{1}{2}A_{3})(du)(A_{2}+\tfrac{1}{2}A_{3})(dv)$
$\displaystyle-2\int_{0}^{\min(u,v)-}\frac{(\sigma_{12}+\tfrac{1}{2}\sigma_{13}+\sigma_{2}^{2}+\tfrac{3}{2}\sigma_{23}+\tfrac{1}{2}\sigma_{3}^{2})(dw)}{1-\Delta
A_{\bullet}(w)}(A_{2}+\tfrac{1}{2}A_{3})(du)A_{\bullet}(dv)$
$\displaystyle+(\sigma_{2}^{2}+\sigma_{23}+\tfrac{1}{4}\sigma_{3}^{2})(\min(u,v))A_{\bullet}(du)A_{\bullet}(dv)\Big{]},$

where $A_{\bullet}=A_{1}+A_{2}+A_{3}$ is the all-cause cumulative hazard
function,

$\sigma_{j}^{2}(t)=\int_{0}^{t}\frac{1-\Delta A_{j}(u)}{S(u-)G(u-)}A_{j}(du),\
j=1,2,3,$

and

$\sigma_{j\ell}(t)=-\sum_{u\leq t}\frac{\Delta A_{j}(u)\Delta
A_{\ell}(u)}{S(u-)G(u-)},\ j\neq\ell,$

are the asymptotic variance and covariance functions of the normalized cause-
specific Nelson-Aalen estimators, respectively, and

$\sigma_{\bullet}^{2}(t)=\int_{0}^{t}\frac{1-\Delta
A_{\bullet}(u)}{S(u-)G(u-)}A_{\bullet}(du)=\sum_{j\neq\ell}(\sigma^{2}_{j}(t)+\sigma_{j\ell}(t)),\quad
t\in[0,\tau],$

is the asymptotic variance function of the normalized all-cause Nelson-Aalen
estimator. We propose to use the consistent Greenwood-type variance and
estimators of $\sigma_{j}^{2}$ and $\sigma_{j\ell}$; see, e.g., formulas
(4.4.17) and (4.4.18) in Andersen et al.5.

For the randomization version of the relative treatment effect estimator,
$\sqrt{n}(\tilde{\theta}_{n}-\tfrac{1}{2})$, we discovered a structure similar
to $\sigma^{2}_{\theta}$, except that all quantities
$\sigma_{j}^{2},\sigma_{j\ell},A_{j}$ are replaced by their randomization
counterparts; see the subsequent subsection for details.

###  B.2 Variance estimators

As motivated in the previous subsection, a consistent estimator for the
variance of the normalized randomized relative treatment effect,
$\sqrt{n}(\tilde{\theta}_{n}-\tfrac{1}{2})$, is given by

$\displaystyle\tilde{\sigma}^{2}_{\theta,n}=\int_{0}^{\tau}\int_{0}^{\tau}\tilde{S}_{n}(u-)\tilde{S}_{n}(v-)\Big{[}\int_{0}^{\min(u,v)-}\frac{\tilde{\sigma}^{2}_{\bullet,n}(dw)}{1-\Delta\tilde{A}_{\bullet,n}(w)}(\tilde{A}_{2,n}+\tfrac{1}{2}\tilde{A}_{3,n})(du)$
$\displaystyle\cdot(\tilde{A}_{2,n}+\tfrac{1}{2}\tilde{A}_{3,n})(dv)-2\int_{0}^{\min(u,v)-}\frac{(\tilde{\sigma}_{12,n}+\tfrac{1}{2}\tilde{\sigma}_{13,n}+\tilde{\sigma}_{2,n}^{2}+\tfrac{3}{2}\tilde{\sigma}_{23,n}+\tfrac{1}{2}\tilde{\sigma}_{3,n}^{2})(dw)}{1-\Delta\tilde{A}_{\bullet,n}(w)}$
$\displaystyle\cdot(\tilde{A}_{2,n}+\tfrac{1}{2}\tilde{A}_{3,n})(du)\tilde{A}_{\bullet,n}(dv)+(\tilde{\sigma}_{2,n}^{2}+\tilde{\sigma}_{23,n}+\tfrac{1}{4}\tilde{\sigma}_{3,n}^{2})(\min(u,v))\tilde{A}_{\bullet,n}(du)\tilde{A}_{\bullet,n}(dv)\Big{]}$

with $\tilde{S}_{n}$ and $\tilde{A}_{j,n}$ being the Kaplan-Meier and the
cause-specific Nelson-Aalen estimators, respectively, based on the randomized
sample $(\tilde{Z}_{i},\tilde{\varepsilon}_{i}),i=1,\dots,n$, and similarly
for the estimators $\tilde{\sigma}_{j\ell,n}$ and $\tilde{\sigma}^{2}_{j,n}$
of $\sigma_{j\ell}$ and $\sigma^{2}_{j}$, respectively, $j=1,2,3$, $\ell\neq
j$. The estimator $\widehat{\sigma}_{\theta,n}^{2}$ of $\sigma_{\theta}^{2}$
is similarly obtained, just based on the original sample
$(Z_{i},\varepsilon_{i}),i=1,\dots,n$.

##  Appendix C Proofs

Proof of Lemma A.1. We are going to show that the censored competing risks
data set exhibits the correct underlying hazard rates. The consistency and
asymptotic normality statements then follow from well-known results in the
literature, e.g., Andersen et al.5. We refer to Dobler12 for detailed
derivations regarding the general case of event times with both discrete and
continuous components.

We first consider continuous components of the hazard rate. Thus, at points of
continuity $t<\tau$ of the distribution of $\check{T}$, we have

$\displaystyle\lim_{dt\downarrow 0}P(Z\in[t,t+dt],\varepsilon=1\ |\ Z\geq
t)/dt$ $\displaystyle=\lim_{dt\downarrow
0}P(\\{X_{1}<X_{2},\delta_{1}=1,X_{1}\in[t,t+dt]\\}$
$\displaystyle\cup\\{X_{1}=X_{2}\in[t,t+dt],\delta_{1}=1,\delta_{2}=0\\}|\
X_{1}\geq t,X_{2}\geq t)/dt$ $\displaystyle=\lim_{dt\downarrow
0}\frac{P(T_{1}\in[t,t+dt],T_{2}\geq t+dt,C_{1}\geq t+dt,C_{2}\geq
t)}{P(T_{1}\geq t,T_{2}\geq t,C_{1}\geq t,C_{2}\geq t)dt}$
$\displaystyle=\lim_{dt\downarrow 0}\frac{P(T_{1}\in[t,t+dt],T_{2}\geq
t+dt)P(C_{1}\geq t+dt,C_{2}\geq t)}{P(T_{1}\geq t,T_{2}\geq t)P(C_{1}\geq
t,C_{2}\geq t)dt}$ $\displaystyle=\lim_{dt\downarrow
0}P(\check{T}\in[t,t+dt],\check{\varepsilon}=1\ |\ \check{T}\geq
t)/dt=\lim_{dt\downarrow 0}dA_{1}(t)/dt.$

A similar derivation holds for time points $t\leq\tau$ of discontinuity:

$\displaystyle P(Z=t,\varepsilon=1\ |\ Z\geq t)$
$\displaystyle=P(X_{1}=t,X_{2}>t,\delta_{1}=1\ |\ X_{1}\geq t,X_{2}\geq t)$
$\displaystyle=\frac{P(\min(T_{1},\tau)=t,T_{2}>t,C_{1}\geq
t,C_{2}>t)}{P(T_{1}\geq t,T_{2}\geq t,C_{1}\geq t,C_{2}\geq t)}$
$\displaystyle=P(\check{T}=t,\check{\varepsilon}=1\ |\ \check{T}\geq
t)=A_{1}(t)-A_{1}(t-).$

Here, $A_{1}(t-)=\lim_{u\uparrow t}A_{1}(u)$ denotes the left-continuous
version of $A_{1}$. The very same arguments can be used for the second
competing risk, i.e., $\check{\varepsilon}=2$.

For the third competing risk, we similarly have for any point $t<\tau$ of
continuity that

$\displaystyle\lim_{dt\downarrow 0}P(Z\in[t,t+dt],\varepsilon=3\ |\ Z\geq
t)/dt$ $\displaystyle=\lim_{dt\downarrow
0}P(X_{1}=X_{2}\in[t,t+dt],\delta_{1}=\delta_{2}=1\ |\ X_{1}\geq t,X_{2}\geq
t)/dt$ $\displaystyle=\lim_{dt\downarrow
0}\frac{P(T_{1}=T_{2}\in[t,t+dt],C_{1}\geq t+dt,C_{2}\geq t+dt)}{P(T_{1}\geq
t,T_{2}\geq 2,C_{1}\geq t,C_{2}\geq t)dt}$ $\displaystyle=\lim_{dt\downarrow
0}\frac{P(T_{1}=T_{2}\in[t,t+dt])P(C_{1}\geq t+dt,C_{2}\geq t+dt)}{P(T_{1}\geq
t,T_{2}\geq t)P(C_{1}\geq t,C_{2}\geq t)dt}$ $\displaystyle=\lim_{dt\downarrow
0}P(T_{1}=T_{2}\in[t,t+dt]\ |\ T\geq t)/dt$ $\displaystyle=\lim_{dt\downarrow
0}P(\check{T}\in[t,t+dt]\ |\ \check{T}\geq t)/dt=\lim_{dt\downarrow
0}dA_{3}(t)/dt.$

Again, a similar derivation holds for times $t\leq\tau$ of discontinuity:

$\displaystyle P(Z=t,\varepsilon=3\ |\ Z\geq t)$
$\displaystyle=P(X_{1}=X_{2}=t,\delta_{1}=\delta_{2}=1\ |\ X_{1}\geq
t,X_{2}\geq t)$
$\displaystyle=\frac{P(\min(T_{1},\tau)=\min(T_{2},\tau)=t,C_{1}\geq
t,C_{2}\geq t)}{P(T_{1}\geq t,T_{2}\geq 2,C_{1}\geq t,C_{2}\geq t)}$
$\displaystyle=\frac{P(\min(T_{1},\tau)=\min(T_{2},\tau)=t)P(C_{1}\geq
t,C_{2}\geq t)}{P(T_{1}\geq t,T_{2}\geq t)P(C_{1}\geq t,C_{2}\geq t)}$
$\displaystyle=P(\min(T_{1},\tau)=\min(T_{2},\tau)=t\ |\ T_{1}\geq t,T_{2}\geq
t)$ $\displaystyle=P(\check{T}=t,\check{\varepsilon}=3\ |\ \check{T}\geq
t)=A_{3}(t)-A_{3}(t-).$

This verifies that the Nelson-Aalen estimators estimate the correct
quantities. $\QED$
Proof of Proposition 2.1. Part (a) is obvious.
(b) Lemma A.1 revealed that the underlying intensity processes are preserved
under the data transformation. In addition, counting and at-risk process are
known to be sufficient statistics for the nonparametric intensity processes2.
It is thus apparent that the $\sigma$-algebra generated by the competing risks
data set is sufficient for $\theta$.
(c) It is well-known that Aalen-Johansen estimators have an interpretation as
NPMLEs; we refer to Sections IV.4.1.5 and IV.1.5 in Andersen et al.5 for a
detailed discussion about the concept of NPMLEs, and the NPMLE properties of
the Aalen-Johansen and the more fundamental Nelson-Aalen estimators,
respectively. Likewise, it clear that $\widehat{\theta}_{n}$ is the NPMLE of
$\theta$. $\QED$
Proof of Theorem 2.2. The stated convergence in distribution is a consequence
of the continuous mapping theorem in combination with the functional delta-
method because the relative treatment effect estimator depends on the cause-
specific Nelson-Aalen estimators through a combination of the following
functionals:

  * •

Wilcoxon functional (Section 3.9.4.1 in van der Vaart and Wellner39). For
càdlàg functions $f$ and functions $g$ of bounded variation (bounded by a
fixed constant), the functional $(f,g)\mapsto\int_{0}^{\tau}f(u)dg(u)$ is
Hadamard-differentiable with derivative
$(h,k)\mapsto\int_{0}^{\tau}h(u)dg(u)+\int_{0}^{\tau}f(u)dk(u)$, where the
latter integral is defined via integration by parts:
$\int_{0}^{\tau}f(u)dk(u)=f(\tau)k(\tau)-f(0)k(0)-\int_{0}^{\tau}k(u-)df(u)$
if $k$ is not of bounded variation.

  * •

Product integral (Section 3.9.4.5 in van der Vaart and Wellner39). For càdlàg
functions $A$ of bounded variation (bounded by a fixed constant), the
functional $A\mapsto\prod_{0<t\leq\tau}(1+A(dt))$ is Hadamard-differentiable
with derivative

$\alpha\mapsto\int_{0}^{\tau}\Big{(}\prod_{0<s<t}(1+A(ds))\Big{)}\alpha(dt)\Big{(}\prod_{t<v\leq\tau}(1+A(dv))\Big{)}$

which is again defined by integration by parts. Here and below, $\prod$ with a
time-continuous indexing stands for the product integral.

Due to the chain rule (Lemma 3.9.3 in van der Vaart and Wellner39), the
composition of both functionals is also Hadamard-differentiable. Thus, since
the Kaplan-Meier estimator satisfies $\widehat{S}_{n}(u)=\prod_{0<s\leq
u}(1+(-\widehat{A}_{1,n}-\widehat{A}_{2,n}-\widehat{A}_{3,n})(du))$, we obtain
the following asymptotic representation of the relative treatment effect
estimator:

$\displaystyle\sqrt{n}(\widehat{\theta}_{n}-\theta)=$
$\displaystyle\int_{0}^{\tau}\sqrt{n}(\widehat{S}_{n}(u-)-S(u-))(A_{2}+\tfrac{1}{2}A_{3})(du)$
$\displaystyle+\int_{0}^{\tau}S(u-)\sqrt{n}\big{[}(\widehat{A}_{2,n}+\tfrac{1}{2}\widehat{A}_{3,n})-(A_{2}+\tfrac{1}{2}A_{3})\big{]}(du)+o_{p}(1),$

where the first integral can be rewritten as

$\displaystyle-\int_{0}^{\tau}S(u-)\int_{0}^{u-}\frac{\sqrt{n}(\widehat{A}_{\bullet,n}-A_{\bullet})(dw)}{1-\Delta
A_{\bullet}(w)}(A_{2}+\tfrac{1}{2}A_{3})(du)+o_{p}(1).$

Here, $\widehat{A}_{\bullet,n}=\sum_{j=1}^{3}\widehat{A}_{j,n}$ is the all-
cause Nelson-Aalen estimator. Simplifying the expression in the previous
display once again, we get

$\displaystyle-\int_{0}^{\tau}S(u-)$
$\displaystyle\Big{[}\sqrt{n}(\widehat{A}_{\bullet,n}-A_{\bullet})(u-)+\sum_{w\in
D_{u-}}\frac{\Delta A_{\bullet}(w)}{1-\Delta
A_{\bullet}(w)}\cdot\Delta\sqrt{n}(\widehat{A}_{\bullet,n}-A_{\bullet})(w)\Big{]}$
$\displaystyle\times(A_{2}+\tfrac{1}{2}A_{3})(du)+o_{p}(1),$

where $D_{u-}\subset[0,u)$ denotes the set of all discontinuities of
$A_{\bullet}$ in $[0,u)$. In combination with the well-known limit behaviour
of the Nelson-Aalen estimators on function spaces, our previous reasoning and
the continuous mapping theorem establishes the asymptotic normality of the
normalized relative treatment effect estimator.

It remains to verify the asymptotic variance $\sigma^{2}_{\theta}$. First note
that, for deterministic càdlàg functions of bounded variation $a$ and $b$ and
for stochastic processes $W$ and $Y$ that satisfy the required integrability
conditions, Fubini’s theorem implies that

$cov\Big{(}\int_{0}^{\tau}W(u)da(u),\int_{0}^{\tau}Y(v)db(v)\Big{)}=\int_{0}^{\tau}\int_{0}^{\tau}cov\big{(}W(u),Y(v)\big{)}da(u)db(v).$

This, in combination with the asymptotic variances $\sigma^{2}_{j}$ and
covariances $\sigma_{j\ell}$ of the normalized cause-specific Nelson-Aalen
estimators and straightforward but tedious computations reveal the claimed
structure of asymptotic variance $\sigma^{2}_{\theta}$. $\QED$
Proof of Theorem 3.1. At first, we verify the conditions of Theorem 2 in
Dobler13 to establish the asymptotic normality of the randomized relative
treatment effect. For this, we notice that the (randomized) Nelson-Aalen
estimators $\tilde{A}_{j,n}$, $j=1,2,3$, are retrieved as functionals of the
empiricial process of $(\tilde{Z}_{i},\tilde{\varepsilon}_{i}),i=1,\dots,n$,
indexed by the class

$\mathcal{F}=\\{(z,e)\mapsto 1\\{z\leq t,e=j\\},\ (z,e)\mapsto 1\\{z\geq
t\\}:\ 0\leq t\leq\tau;\ j=1,2,3\\}.$

It is easy to see that $\mathcal{F}$ is a Vapnik-C̆ervonenkis class and it
hence satisfies the uniform entropy condition (2.5.1) in 39. As a consequence,
$\mathcal{F}$ is both $\mathbb{P}$ and $\tilde{\mathbb{P}}$-Donsker and both
distributions have bounded supremum norms, where
$\mathbb{P}=P^{(Z,\varepsilon)}$ and
$\tilde{\mathbb{P}}=P^{(\tilde{Z},\tilde{\varepsilon})}$. Similarly, the class

$\tilde{\mathcal{F}}=\\{(z,e)\mapsto\tfrac{1}{2}[f(z,j)+f(z,4-j)]:\
f\in\mathcal{F}\\}$

is $\mathbb{P}$-Donsker and $\mathbb{P}$ again has a bounded supremum norm
with respect to $\tilde{\mathcal{F}}$. Thus, Theorem 2 in 13 yields that the
conditional distribution of $\sqrt{n}(\tilde{\theta}_{n}-0.5)$ converges
weakly to a normal distribution in probability as $n\to\infty$.

The limit variance, say $\tilde{\sigma}_{\theta}^{2}$, is in general different
from $\sigma_{\theta}^{2}$. That is why the studentization based on
$\tilde{\sigma}_{\theta,n}^{2}$ is essential. To justify the consistency of
$\tilde{\sigma}_{\theta,n}^{2}$ for $\tilde{\sigma}_{\theta}^{2}$, note that
this variance estimator is a continuous functional of the randomization
empirical process, i.e., the empirical process based on
$(\tilde{Z}_{i},\tilde{\varepsilon}_{i}),i=1,\dots,n$. Since the class
$\mathcal{F}$ is Glivenko-Cantelli, the continuous mapping theorem yields the
consistency of the variance estimator. Finally, an application of a
conditional version of Slutkzy’s theorem concludes the proof. $\QED$

##  Appendix D Additional simulation results

###  D.1 Sizes under $H_{0}:\theta=0.5$ and nominal significance levels
$\alpha\in\\{1\%,10\%\\}$

Tables 3 and 4 summarize the simulation results for the size of the proposed
one-sided tests under $H_{0}:\theta=0.5$ with nominal significance levels
$\alpha=1\%$ and $10\%$, respectively. The remaining simulation scenarios are
the same as those which led to the results displayed in Table 1 for
$\alpha=5\%$ in the main manuscript. As the overall impression is similar to
the findings from that Table 1, no additional comments on the outcomes are
given here.

###  D.2 Additional power simulation results

Figures 4 and 5 present the power simulation results of the proposed
randomization-based one-sided test in comparison to the paired Prentice-
Wilcoxon test and the stratified log-rank test. Here, the marginal
distributions are different competing exponential distributions and
exponential-Gompertz mixtures versus an exponential distribution,
respectively. Since also here the overall findings are similar to those given
in Figure 2 in the main manuscript, no further comments are needed.

Figure 4: Simulated power of three selected right-tailed randomization-based
tests with significance level $\alpha=5\%$. Figure 5: Simulated power of three
selected right-tailed randomization-based tests with significance level
$\alpha=5\%$.

Table 3: Simulated sizes of the right-tailed tests (in %) with nominal
significance level $\alpha=1\%$.
Abbreviations: Copula: GH = Gumbel-Hougaard; critical values: asy. = asymptotical normal, bs. = bootstrap, rand. = randomization; tests: lin. = linear, tra. = $\log-\log$-transformed. |  |  | light censoring | medium censoring | strong censoring
---|---|---|---|---|---
|  |  | asy. | bs. | rand. | asy. | bs. | rand. | asy. | bs. | rand.
copula | distribution | $n$ | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra.
GH | Exm vs. | $25$ | $2.5$ | $0.4$ | $0.3$ | $2.2$ | $1.2$ | $1.2$ | $3.2$ | $0.6$ | $1.2$ | $2.9$ | $1.5$ | $1.4$ | $3.2$ | $0.6$ | $1.5$ | $2.5$ | $1.3$ | $1.3$
| Exp mix | $50$ | $1.6$ | $0.6$ | $0.5$ | $1.4$ | $1.1$ | $1.1$ | $2.1$ | $0.9$ | $0.9$ | $1.9$ | $1.4$ | $1.4$ | $2.7$ | $1.3$ | $1.5$ | $2.6$ | $1.4$ | $1.4$
|  | $75$ | $1.3$ | $0.6$ | $0.5$ | $1.0$ | $0.9$ | $0.9$ | $1.9$ | $0.7$ | $0.8$ | $1.6$ | $1.3$ | $1.3$ | $2.1$ | $1.3$ | $1.3$ | $1.8$ | $1.5$ | $1.5$
|  | $100$ | $1.5$ | $0.9$ | $1.0$ | $1.3$ | $1.3$ | $1.3$ | $1.4$ | $0.6$ | $0.7$ | $1.1$ | $0.9$ | $0.9$ | $2.4$ | $1.5$ | $1.2$ | $1.7$ | $1.6$ | $1.6$
|  | $125$ | $1.1$ | $0.6$ | $0.7$ | $1.0$ | $1.0$ | $1.0$ | $1.3$ | $0.8$ | $0.9$ | $1.2$ | $1.0$ | $1.0$ | $2.4$ | $1.8$ | $1.4$ | $1.8$ | $1.8$ | $1.8$
|  | $150$ | $1.4$ | $0.9$ | $1.1$ | $1.2$ | $1.3$ | $1.3$ | $1.2$ | $0.9$ | $0.9$ | $1.0$ | $1.1$ | $1.1$ | $2.1$ | $1.5$ | $1.4$ | $1.7$ | $1.3$ | $1.2$
GH | Gompertz | $25$ | $2.7$ | $0.9$ | $0.2$ | $0.7$ | $1.6$ | $1.7$ | $4.6$ | $1.3$ | $0.9$ | $1.3$ | $2.9$ | $3.1$ | $6.7$ | $2.1$ | $2.4$ | $1.3$ | $5.9$ | $5.9$
| vs. Exp | $50$ | $1.5$ | $0.6$ | $0.1$ | $0.9$ | $1.1$ | $1.1$ | $2.0$ | $1.0$ | $0.5$ | $1.1$ | $1.9$ | $1.8$ | $3.8$ | $2.0$ | $0.9$ | $2.3$ | $5.1$ | $5.0$
|  | $75$ | $1.3$ | $0.6$ | $0.5$ | $1.0$ | $1.0$ | $1.0$ | $1.4$ | $0.7$ | $0.4$ | $0.9$ | $1.4$ | $1.4$ | $3.1$ | $1.8$ | $0.8$ | $1.7$ | $4.7$ | $4.6$
|  | $100$ | $1.3$ | $0.6$ | $0.7$ | $1.1$ | $1.0$ | $1.0$ | $1.6$ | $0.8$ | $0.5$ | $0.9$ | $1.3$ | $1.4$ | $1.7$ | $1.0$ | $0.4$ | $0.9$ | $3.7$ | $3.7$
|  | $125$ | $1.4$ | $0.6$ | $0.9$ | $1.2$ | $1.1$ | $1.1$ | $1.2$ | $0.6$ | $0.6$ | $0.8$ | $1.0$ | $1.0$ | $1.7$ | $0.8$ | $0.4$ | $0.7$ | $3.9$ | $3.8$
|  | $150$ | $1.4$ | $0.9$ | $1.1$ | $1.4$ | $1.4$ | $1.4$ | $1.2$ | $0.6$ | $0.7$ | $0.9$ | $1.0$ | $1.0$ | $1.5$ | $0.8$ | $0.5$ | $0.7$ | $3.0$ | $3.0$
Clayton | Exp vs. | $25$ | $2.2$ | $0.5$ | $0.1$ | $2.4$ | $1.1$ | $1.2$ | $2.2$ | $0.4$ | $0.4$ | $2.6$ | $1.0$ | $0.9$ | $2.5$ | $0.4$ | $0.6$ | $2.6$ | $1.0$ | $1.0$
| Exp mix | $50$ | $1.6$ | $0.7$ | $0.4$ | $1.6$ | $1.2$ | $1.2$ | $1.7$ | $0.4$ | $0.5$ | $1.4$ | $0.9$ | $0.9$ | $1.9$ | $0.8$ | $0.8$ | $1.9$ | $1.1$ | $1.1$
|  | $75$ | $1.1$ | $0.5$ | $0.5$ | $0.9$ | $1.0$ | $1.0$ | $1.4$ | $0.7$ | $0.6$ | $1.3$ | $1.1$ | $1.1$ | $1.5$ | $0.8$ | $0.9$ | $1.3$ | $1.2$ | $1.2$
|  | $100$ | $1.6$ | $0.9$ | $1.0$ | $1.4$ | $1.4$ | $1.4$ | $1.4$ | $0.9$ | $0.7$ | $1.1$ | $1.2$ | $1.2$ | $1.3$ | $0.7$ | $0.7$ | $1.1$ | $1.1$ | $1.1$
|  | $125$ | $1.0$ | $0.6$ | $0.8$ | $0.9$ | $0.8$ | $0.9$ | $1.5$ | $0.9$ | $0.9$ | $1.3$ | $1.4$ | $1.4$ | $1.3$ | $0.8$ | $0.8$ | $1.1$ | $1.1$ | $1.2$
|  | $150$ | $1.0$ | $0.6$ | $0.7$ | $1.0$ | $0.9$ | $0.9$ | $1.2$ | $0.7$ | $0.8$ | $1.0$ | $1.0$ | $1.0$ | $1.2$ | $0.8$ | $0.8$ | $1.0$ | $1.1$ | $1.1$
Clayton | Gompertz | $25$ | $2.3$ | $0.6$ | $0.1$ | $2.3$ | $1.4$ | $1.5$ | $2.5$ | $0.4$ | $0.4$ | $2.6$ | $1.1$ | $1.0$ | $3.3$ | $0.3$ | $0.6$ | $3.1$ | $1.4$ | $1.4$
| vs. Exp | $50$ | $1.3$ | $0.6$ | $0.5$ | $1.1$ | $0.9$ | $0.9$ | $2.0$ | $0.8$ | $0.6$ | $1.7$ | $1.2$ | $1.1$ | $1.9$ | $0.7$ | $0.6$ | $1.8$ | $1.1$ | $1.1$
|  | $75$ | $1.0$ | $0.6$ | $0.4$ | $0.7$ | $0.8$ | $0.8$ | $1.6$ | $0.8$ | $0.8$ | $1.3$ | $1.1$ | $1.1$ | $1.5$ | $0.8$ | $0.7$ | $1.2$ | $1.1$ | $1.1$
|  | $100$ | $1.4$ | $0.8$ | $0.8$ | $1.2$ | $1.1$ | $1.1$ | $1.4$ | $0.8$ | $0.8$ | $1.0$ | $1.1$ | $1.1$ | $1.7$ | $1.0$ | $0.9$ | $1.3$ | $1.4$ | $1.4$
|  | $125$ | $1.0$ | $0.7$ | $0.7$ | $1.0$ | $1.0$ | $1.0$ | $1.4$ | $0.8$ | $0.8$ | $1.1$ | $1.2$ | $1.2$ | $1.7$ | $1.1$ | $1.0$ | $1.3$ | $1.5$ | $1.4$
|  | $150$ | $1.3$ | $0.8$ | $0.8$ | $1.0$ | $1.0$ | $1.0$ | $1.4$ | $0.8$ | $0.8$ | $1.0$ | $1.2$ | $1.2$ | $1.5$ | $1.2$ | $1.0$ | $1.3$ | $1.4$ | $1.4$

Table 4: Simulated sizes of the right-tailed tests (in %) with nominal
significance level $\alpha=10\%$.
Abbreviations: Copula: GH = Gumbel-Hougaard; critical values: asy. = asymptotical normal, bs. = bootstrap, rand. = randomization; tests: lin. = linear, tra. = $\log-\log$-transformed. |  |  | light censoring | medium censoring | strong censoring
---|---|---|---|---|---
|  |  | asy. | bs. | rand. | asy. | bs. | rand. | asy. | bs. | rand.
copula | distribution | $n$ | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra. | lin. | tra.
GH | Exp vs. | $25$ | $12.3$ | $10.1$ | $8.8$ | $10.4$ | $10.7$ | $10.6$ | $11.1$ | $9.2$ | $7.8$ | $9.4$ | $10.5$ | $10.5$ | $9.1$ | $7.9$ | $5.4$ | $6.7$ | $12.8$ | $12.7$
| Exp mix | $50$ | $10.9$ | $9.6$ | $9.4$ | $10.0$ | $10.4$ | $10.4$ | $10.4$ | $9.2$ | $8.9$ | $9.6$ | $10.3$ | $10.2$ | $7.2$ | $6.8$ | $4.7$ | $5.2$ | $12.4$ | $12.4$
|  | $75$ | $10.5$ | $9.4$ | $9.5$ | $10.2$ | $10.1$ | $10.1$ | $11.0$ | $9.8$ | $10.0$ | $10.6$ | $10.7$ | $10.7$ | $6.8$ | $6.1$ | $4.9$ | $5.1$ | $12.3$ | $12.3$
|  | $100$ | $10.7$ | $9.9$ | $10.1$ | $10.5$ | $10.4$ | $10.4$ | $10.6$ | $9.6$ | $9.7$ | $10.1$ | $10.1$ | $10.1$ | $7.0$ | $6.6$ | $5.3$ | $5.5$ | $12.3$ | $12.3$
|  | $125$ | $9.8$ | $8.9$ | $9.1$ | $9.4$ | $9.6$ | $9.6$ | $10.1$ | $9.3$ | $9.6$ | $9.8$ | $9.6$ | $9.6$ | $6.8$ | $6.3$ | $5.4$ | $5.6$ | $12.8$ | $12.8$
|  | $150$ | $10.7$ | $9.9$ | $10.4$ | $10.6$ | $10.5$ | $10.5$ | $10.4$ | $9.6$ | $9.8$ | $10.1$ | $9.9$ | $9.9$ | $6.7$ | $6.2$ | $5.7$ | $5.7$ | $12.7$ | $12.7$
GH | Gompertz | $25$ | $13.3$ | $11.0$ | $7.1$ | $9.4$ | $12.0$ | $12.0$ | $14.2$ | $12.0$ | $6.7$ | $9.4$ | $14.6$ | $14.8$ | $16.0$ | $13.9$ | $6.6$ | $9.3$ | $22.9$ | $22.7$
| vs. Exp | $50$ | $11.2$ | $9.7$ | $8.9$ | $9.8$ | $10.5$ | $10.5$ | $10.8$ | $9.2$ | $6.6$ | $7.4$ | $10.9$ | $10.9$ | $12.2$ | $10.9$ | $5.4$ | $7.0$ | $19.3$ | $19.3$
|  | $75$ | $9.7$ | $8.4$ | $8.5$ | $8.9$ | $9.1$ | $9.2$ | $10.2$ | $8.9$ | $7.7$ | $8.6$ | $10.1$ | $10.1$ | $10.5$ | $9.5$ | $5.8$ | $6.5$ | $18.6$ | $18.6$
|  | $100$ | $10.5$ | $9.4$ | $9.6$ | $10.0$ | $10.1$ | $10.2$ | $10.1$ | $9.0$ | $8.6$ | $9.1$ | $10.3$ | $10.3$ | $9.0$ | $8.3$ | $5.1$ | $5.7$ | $15.1$ | $15.2$
|  | $125$ | $10.8$ | $9.7$ | $10.1$ | $10.4$ | $10.5$ | $10.5$ | $10.4$ | $9.2$ | $9.5$ | $9.9$ | $10.1$ | $10.2$ | $8.7$ | $8.1$ | $5.5$ | $5.8$ | $15.4$ | $15.4$
|  | $150$ | $10.0$ | $9.4$ | $9.5$ | $9.8$ | $9.8$ | $9.8$ | $9.8$ | $9.0$ | $9.4$ | $9.7$ | $9.7$ | $9.7$ | $8.1$ | $7.3$ | $5.3$ | $5.7$ | $14.9$ | $14.8$
Clayton | Exp vs. | $25$ | $11.1$ | $8.9$ | $7.9$ | $9.3$ | $9.8$ | $9.8$ | $11.2$ | $8.9$ | $7.1$ | $8.9$ | $9.9$ | $9.9$ | $10.9$ | $8.7$ | $7.1$ | $9.0$ | $9.9$ | $9.9$
| Exp mix | $50$ | $9.9$ | $8.6$ | $8.3$ | $9.0$ | $9.6$ | $9.6$ | $10.6$ | $9.1$ | $8.6$ | $9.4$ | $10.0$ | $10.0$ | $9.6$ | $8.6$ | $7.0$ | $7.7$ | $9.9$ | $9.9$
|  | $75$ | $10.6$ | $9.3$ | $9.4$ | $9.9$ | $10.2$ | $10.2$ | $9.2$ | $8.1$ | $7.6$ | $8.2$ | $9.2$ | $9.3$ | $9.7$ | $8.6$ | $7.5$ | $8.1$ | $10.3$ | $10.3$
|  | $100$ | $10.8$ | $10.0$ | $10.1$ | $10.5$ | $10.7$ | $10.7$ | $9.4$ | $8.5$ | $8.6$ | $8.9$ | $9.7$ | $9.7$ | $9.0$ | $8.1$ | $7.4$ | $7.8$ | $10.1$ | $10.1$
|  | $125$ | $9.5$ | $8.8$ | $9.1$ | $9.5$ | $9.6$ | $9.6$ | $10.0$ | $9.2$ | $9.0$ | $9.3$ | $10.3$ | $10.3$ | $9.7$ | $8.9$ | $8.2$ | $8.4$ | $11.1$ | $11.1$
|  | $150$ | $10.2$ | $9.4$ | $9.7$ | $9.9$ | $10.2$ | $10.2$ | $9.5$ | $8.9$ | $8.8$ | $9.0$ | $10.3$ | $10.3$ | $9.1$ | $8.5$ | $7.8$ | $8.1$ | $10.6$ | $10.6$
Clayton | Gompertz | $25$ | $12.1$ | $9.6$ | $7.9$ | $9.5$ | $10.1$ | $10.2$ | $12.2$ | $9.7$ | $7.4$ | $9.6$ | $10.8$ | $10.8$ | $12.6$ | $10.3$ | $7.7$ | $10.1$ | $11.4$ | $11.4$
| vs. Exp | $50$ | $11.0$ | $9.5$ | $8.7$ | $9.5$ | $10.0$ | $10.0$ | $11.9$ | $10.4$ | $9.3$ | $10.2$ | $11.0$ | $11.0$ | $11.3$ | $9.6$ | $8.5$ | $9.2$ | $11.2$ | $11.3$
|  | $75$ | $10.1$ | $8.9$ | $8.7$ | $9.2$ | $9.4$ | $9.4$ | $11.2$ | $10.1$ | $9.3$ | $9.9$ | $11.1$ | $11.1$ | $10.6$ | $9.5$ | $8.4$ | $9.1$ | $11.4$ | $11.5$
|  | $100$ | $10.7$ | $9.7$ | $9.7$ | $10.1$ | $10.4$ | $10.4$ | $10.9$ | $9.9$ | $9.3$ | $9.8$ | $10.5$ | $10.5$ | $10.8$ | $9.7$ | $8.9$ | $9.3$ | $11.3$ | $11.3$
|  | $125$ | $9.6$ | $8.7$ | $8.9$ | $9.1$ | $9.4$ | $9.4$ | $9.3$ | $8.6$ | $8.2$ | $8.5$ | $9.4$ | $9.4$ | $10.4$ | $9.6$ | $8.5$ | $8.9$ | $10.8$ | $10.8$
|  | $150$ | $10.4$ | $9.6$ | $9.5$ | $9.8$ | $10.1$ | $10.1$ | $11.0$ | $10.2$ | $9.8$ | $10.0$ | $10.8$ | $10.8$ | $10.1$ | $9.6$ | $8.7$ | $9.0$ | $11.0$ | $11.0$