A Discontinuity Adjustment

arXiv:1702.01081v1 [] 3 Feb 2017

for Subdistribution Function Confidence Bands
Applied to Right-Censored Competing Risks Data
Dennis Dobler
Ulm University, Institute of Statistics,
Helmholtzstr. 20, 89081 Ulm, Germany
email: dennis.dobler@uni-ulm.de

February 6, 2017
Summary
The wild bootstrap is the resampling method of choice in survival analytic applications. Theoretic justifications rely on the assumption of existing intensity functions which is equivalent to an
exclusion of ties among the event times. However, such ties are omnipresent in practical studies. It turns out that the wild bootstrap should only be applied in a modified manner that corrects
for altered limit variances and emerging dependencies. This again ensures the asymptotic exactness of inferential procedures. An analogous necessity is the use of the Greenwood-type variance
estimator for Nelson-Aalen estimators which is particularly preferred in tied data regimes. All
theoretic arguments are transferred to bootstrapping Aalen-Johansen estimators for cumulative incidence functions in competing risks. An extensive simulation study as well as an application to
real competing risks data of male intensive care unit patients suffering from pneumonia illustrate
the practicability of the proposed technique.
Keywords: Aalen-Johansen estimator; Counting process; Discontinuous cumulative hazard functions; Discontinuous cumulative incidence functions; Greenwood-type variance estimator; NelsonAalen estimator; Survival analysis; Tied event times; Wild bootstrap.

1

Introduction

In survival analysis, the wild bootstrap is a frequently utilized resampling technique; cf. Lin
(1997), Martinussen and Scheike (2006), Beyersmann et al. (2013). Typically, large sample properties are verified by relying on the assumption of existing transition intensities. In the context of
a survival time T this implies the existence of a function
αptq “ lim

∆tÓ0

1
P pT P rt, t ` ∆ts | T ě tq
∆t

or, similarly, the existence of a Lebesgue density of the distribution P T . Without any doubt, this
assumption is too strict in practical applications where events are usually recorded on a discrete
time lattice, e.g. on a daily basis. For estimating unknown limit variances of Nelson-Aalen estimators it is well-known (e.g. Andersen et al. 1993) that Greenwood-type estimators should be
utilized in the presence of ties. Allignol et al. (2010) even found a general preference (also under
left-truncation) for this kind of estimator.
From now on, we assume that there are discrete components in the event time distribution and
that the event, if observed, can be classified to one out of k different causes, i.e. competing risks.
Suppose that n P N i.i.d. individuals participate in our study, but that their observability may
be independently right-censored by i.i.d. censoring variables. Based on these observations, all
available collected information is contained in the counting and at risk processes,
Nji ptq “ 1t“ individual i is observed to fail in r0, ts due to risk j” u
and

Yi ptq “ 1t“ individual i is under observation at time t ´ ” u

respectively, with j “ 1, . . . , k and i “ 1, . . . , n. This notation may be used to extend the arguments below to also incorporate independent left-truncation in the sense of Andersen et al.
(1993),
Chapter III. Define H1uc ptq “ EpN11 ptqq and H̄ptq “ EpY1 ptqq. Now, the sum N1 ptq “
řn
i“1 N1i ptq has the compensator
Λ1 ptq “

n
ÿ
i“1

żt
Λ1i ptq “

Y puqdA1 puq “
0

n żt
ÿ

Yi puqdA1 puq,

i“1 0

şt
where A1 ptq “ 0 H̄puq´1 dH1uc puq is the cumulative hazard function for a type 1 event. Therefore,
pM1 ptq “ N1 ptq ´ Λ1 ptqqt is a square-integrable martingale. The Nelson-Aalen estimator for A1 ptq
ş
p1 ptq “ t Y ´1 dN1 , where we let 0{0 “ 0. The limit
is defined as the counting process integral A
0
?
p1 ´A1 q is determined by the predictable
covariance function of the normalized process W1 “ npA
and optional variation processes of the square-integrable martingales M1i : t ÞÑ N1i ptq ´ Λ1i ptq.
1

In the possibly discontinuous case, these are given by
żt
∆Λ1i puqdΛ1i puq
żt
żt
rM1i sptq “ N1i ptq ´ 2 ∆Λ1i puqdN1i puq ` ∆Λ1i puqdΛ1i puq,

xM1i yptq “ Λ1i ptq ´

0

and

0

0

respectively, cf. Andersen et al. (1993), Section II.4. First, we would like to heuristically discuss
again why Greenwood-type variances estimators are a good choice if ties in theşevent times are post
sible. This is as well seen with the help of covariation processes: Let A˚1 : t ÞÑ 0 JpuqA1 pduq{Y puq
p1 , where Jpuq “ 1tY puq ą 0u. By the
denote the compensator of the Nelson-Aalen estimator A
usual martingale theory, we have the optional covariation process
”
?
˚
p
r npA1 ´ A1 qsptq “ n

ż¨
0

JdM1 ı
ptq “ n
Y

żt
0

n żt
ÿ
J
J
drM1 s “ n
drM1i s,
2
Y
Y2
i“1 0

where the last equality follows from the orthogonality of all involved single martingales. Inserting
the above representation for the optional variation process yields
n
ÿ
?
˚
p
r npA1 ´ A1 qsptq “ n

n żt
n żt
ÿ
ÿ
J
J
J
dN1i ´ 2n
∆Λ1i dN1i ` n
∆Λ1i dΛ1i
2
2
Y
Y
Y2
i“1 0
i“1 0
i“1 0
żt
żt
żt
J
J
J
dN1 ´ 2n
∆A1 dN1 ` n
∆A1 dA1 .
“n
2
2
0 Y
0 Y
0 Y

żt

In a last step, we replace the unknown increments of the cumulative hazard function ∆A1 puq with
p1 puq “ Y puq´1 ∆N1 puq and obtain the Greenwood-type variance
Nelson-Aalen increments ∆A
estimator
żt
p12 ptq
σ

“n
0

J
dN1 ´ 2n
Y2

żt
0

J ∆N1
dN1 ` n
Y2 Y

żt
0

J ∆N1 dN1
“n
Y Y Y

żt
0

Y ´ ∆N1
dN1 .
Y3

One could as well proceed similarly by starting with the predictable variation process instead.
Now, for the large sample properties, let K ą 0 be a terminal time such that A1 pKq ă 8.
A simple application of Glivenko-Cantelli theorems in combination with the continuous mapping
şt
theorem show that this Greenwood-type estimator is uniformly consistent for σ12 ptq “ 0 p1 ´
∆A1 puqqH̄puq´1 dA1 puq on r0, Ks. In fact, this limit is the asymptotic variance of the normalized
Nelson-Aalen estimator in the general set up where ties are possible; see e.g. van der Vaart and
Wellner (1996), Example 3.9.19, for such a result in the classical survival set up and Appendix A
for a detailed calculation for the present competing risks situation. The subtraction of ∆A1 in
the limit variance is crucially different from the continuous case in which the limit variance func2

şt
tion reduces to σ̃12 : t ÞÑ 0 H̄ ´1 dA1 . It is precisely this difference which is the cause for the
inconsistency of the usual wild bootstrap procedure applied to Nelson-Aalen estimators and, as a
consequence, also when applied to Aalen-Johansen estimators of cumulative incidence functions.
The present article is organized as follows: In Section 2, we first discuss the implications of
the usual wild bootstrap procedure, propose a discontinuity adjustment for the wild bootstrap,
and finally present a conditional central limit theorem for this new technique. In the following
Section 3, we show that the Nelson-Aalen estimators for different competing risks are in general
asymptotically dependent. Therefore, we present in a next step an extension of the first proposal
for the wild bootstrap adjustment that guarantees the correct limit dependence structure between
the components for different risks. This technique has some direct implications on resampling the
Aalen-Johansen estimator for cumulative incidence functions as these depend on all cause-specific
hazard functions and, therefore, also on their dependencies. We present conditional central limit
theorems corresponding to this set-up in Section 4, where also variance estimators for these AalenJohansen estimators and time-simultaneous confidence bands for cumulative incidence functions
are deduced. The performance of these bands in terms of coverage probabilities is analyzed in a
simulation study in Section 5 and there it is compared to the behaviour of confidence bands based
on the usual, unadjusted wild bootstrap. In this connection, we consider different variations of
discretization coarseness and discretization probabilities. All considered resampling techniques
are applied to a real data example with competing risks in Section 6, where confidence bands for
the probability of an alive discharge of male patients with pneumonia from intensive care units are
constructed. We conclude with a small discussion in Section 7. All proofs and various detailed
derivations are presented in Appendices A–E.

2

Proposed Wild Bootstrap Adjustment for a Univariate NelsonAalen Estimator

The wild bootstrap is typically applied in the following way (see e.g. Lin 1997, Martinussen and
Scheike 2006, Beyersmann et al. 2013 or Bluhmki et al. 2016 for resampling more general counting
process-based Nelson-Aalen estimators): First, the (independent) martingale increments dM1i puq
? p
in the martingale representation of npA
1 ´ A1 q are replaced by independently weighted counting
process increments, i.e., by ξi dN1i puq. Here, the wild bootstrap weights ξi are i.i.d. with zero mean
and variance 1. Bluhmki et al. (2016) argued that the resulting wild bootstrapped Nelson-Aalen
estimator given by
żt
n
? ÿ
dN1i
x
W̃1 ptq “ n ξi
Y
0
i“1
3

constitutes a square-integrable martingale with predictable and optional variation processes
żt

dN1 p2
“ σ̃1 ptq
2
0 Y
żt
n
ÿ
dN1i p
x
2
p12 ptq.
t ÞÝÑ rW̃1 sptq “ n ξi
“ σ̃
2
0 Y
i“1

x yptq “ n
t ÞÝÑ xW̃
1
and

When utilized as variance estimators, both variation processes in the previous display are uniformly
consistent for σ̃ 2 . Unfortunately, this is even true in the discontinuous case where σ 2 ‰ σ̃ 2 . Therefore, the wild bootstrap procedure, as it is presently applied in the literature, fails to be consistent
for the correct limit distribution of the Nelson-Aalen estimator whenever ties are an inevitable
phenomenon in the analyzed data set. This is exactly the reason why this wild bootstrap approach
consistently overestimates the true variance in case of ties. This phenomenon does not vanish as
n Ñ 8 and it is stronger pronounced for coarser lattices of discrete event times; see the simulation
study in Section 5. This problem occurs while resampling Aalen-Johansen estimators for cumulative incidence functions in competing risks data and it occurs as well while resampling general
Nelson-Aalen estimators if ties are present.
This non-trivial problem calls for a general solution. In the present article, we exemplify the
subsequent solution in the right-censored competing risks set-up. Extensions and modifications to
wild bootstrap versions of more general Nelson-Aalen estimators or of Aalen-Johansen estimators
in general Markovian situations may be obtained in a similar manner, but the limit variances will
be much more complicated. The crucial defect in the wild bootstrap resampling scheme described
above is that the martingale increments dM1i puq should not be replaced by ξi dN1i puq but rather by
something which – considered again as a martingale – reproduces the correct covariance structure.
Therefore, we replace
d
dM1i puq by

ξi

Y puq ´ ∆N1 puq
dN1i puq,
Y puq

which results in the following wild bootstrap resampling version of the normalized Nelson-Aalen
estimator:
żtd
n
? ÿ
Y puq ´ ∆N1 puq dN1i puq
x1 ptq “ n ξi
W
.
Y puq
Y puq
0
i“1
x1 ptqqtPr0,Ks is a martingale with
In a way similar to Bluhmki et al. (2016), one can show that pW
respect to the filtration pFt qtPr0,Ks given by
Ft “ σpξi N1i puq, N1i pvq, Yi pvq : u P r0, ts, v P r0, Ks, i “ 1, . . . , nq.
4

Analogously, it is easy to see that its predictable and optional variation processes are given by
żt

Y puq ´ ∆N1 puq
dN1 puq
Y 3 puq
0
żt
n
ÿ
Y puq ´ ∆N1 puq
2
2
p
x
p1 : t ÞÝÑ rW1 sptq “n ξi
dN1i puq.
and σ
Y 3 puq
0
i“1
p12
σ

x1 yptq “n
: t ÞÝÑ xW

p12 is uniformly consistent for σ12 . AsIt is well-known that the Greenwood-type variance estimator σ
suming the existence of the fourth moments of ξi , a simple application of Chebyshev’s inequality
shows the (conditional) consistency of the second estimator; the uniformity in the conditional convergence in probability follows from a Pólya-type argument. The conditional weak convergence of
the finite-dimensional marginal distributions of the wild bootstrapped Nelson-Aalen process follows easily by an application of Theorem A.1 in Beyersmann et al. (2013). This also shows that
the proposed wild bootstrap approach succeeds in maintaining the correct asymptotic covariance
function which had been our aim in the first place.
Denote by D1 “ tt P r0, Ks : ∆A1 ptq ą 0u the subset of time points for which ties among
the type 1 events are possible. Similarly, we define D2 , . . . , Dk for the other risks. Throughout the
article, we assume the following technical condition in order to conclude the conditional tightness
of the wild bootstrapped Nelson-Aalen process.
Condition 2.1. D “

Ťk
j“1

Dj has a finite cardinality.

In practical applications, this assumption is naturally satisfied: A finite end-of-study time and
measurements on a daily or weekly basis result in a finite lattice. A proof of conditional tightness
finally yields the following conditional central limit theorem for the Nelson-Aalen process:
Theorem 2.2. Assume Condition 2.1. Given F0 and as n Ñ 8, we have the following conditional
weak convergence
d
x1 ÝÑ
W
U1 „ Gaussp0, σ12 q

in outer probability

on the càdlàg function space Dr0, Ks equipped with the supremum distance topology, where U1
is a Gaussian zero-mean martingale with variance function t ÞÑ σ12 ptq. That is, the modified
wild bootstrap succeeds in reproducing the same limit process of the Nelson-Aalen process, in
particular, if ties are present.
The proof is given in Appendix B. When resampling a functional of a multivariate NelsonAalen estimator such as the Aalen-Johansen estimator, it is mandatory to take the covariance structure between all cause-specific Nelson-Aalen estimators into account. In order to reflect this in the
5

resampling scheme, a further adjustment needs to be done as shown in the following section. We
conclude this section with an application of the present approach to the Kaplan-Meier estimator.
Remark 2.3. Consider the case of only one risk, i.e. k “ 1 and W “ W1 . The Kaplan-Meier
estimator for the survival function Sptq “ P pT ą tq is
p “
Sptq

ź

p
p1 ´ Apduqq,

t P r0, Ks.

0ďuďt

It exhibits the martingale representation
?
p ´ Sptqq “ Sptq
WS ptq “ npSptq

żt
0

W pduq
` op p1q.
1 ´ ∆Apuq

Thus, the discontinuity-adjusted wild bootstrapped normalized Kaplan-Meier estimator is
ż¨
xS “ Sp
W
0

x pduq
W
p
1 ´ ∆Apuq

and Theorem 2.2 in combination with the continuous mapping theorem yields the correct limit
process distribution, i.e. a zero-mean Gaussian process with covariance function given by
żt
ps, tq ÞÑ SpsqSptq
0

dApuq
;
p1 ´ ∆ApuqqH̄puq

cf. e.g. Example 3.9.31 in van der Vaart and Wellner (1996).

3

Extension to the Joint Convergence of Multiple Nelson-Aalen
Estimators in Competing Risks Models

We now extend the above martingale notation to all k P N competing risks, i.e. to M1i , . . . , Mki ,
i “ 1, . . . , n. For analyses involving two different cumulative cause-specific hazard functions
in a competing risks set-up, it is also important to take the asymptotic covariance structure of
both Nelson-Aalen estimators into account. In the absolutely continuous case, this asymptotic
covariance function vanishes as all Nelson-Aalen estimators are asymptotically independent; cf.
Andersen et al. (1993), Theorem IV.1.2. In the presence of ties, however, the situation is quite
6

different: Here we have for the martingales M1i and M2i of Section 1 that
żt
Yi puq∆A1 puqdA2 puq,

xM1i , M2i yptq “ ´
ż0t

0

0

0

Yi puq∆A1 puqdA2 puq,

Yi puq∆A2 puqdM1i puq `

Yi puq∆A1 puqdM2i puq ´

rM1i , M2i sptq “ ´

żt

żt

? p
cf. the derivations in Section II.4 in Andersen et al. (1993). Define Wj “ npA
j ´ Aj q, j “
1, . . . , k. The above variation processes are strong evidence that the normalized Nelson-Aalen
estimators W1 , W2 , . . . , Wk are not asymptotically independent anymore in the presence of ties.
This is indeed the case:
Theorem 3.1. As n Ñ 8, we have on the product càdlàg function space Dk r0, Ks, equipped with
the max-sup-norm, that
d

pW1 , W2 , . . . , Wk q ÝÑ pU1 , U2 , . . . , Uk q,
where U1 , U2 , . . . , Uk are zero-mean Gaussian martingales with variance functions
żt
t ÞÑ

σj2 ptq

“
0

1 ´ ∆Aj puq
dAj puq,
H̄puq

j “ 1, 2, . . . , k

and covariance functions (for j ‰ `)
ż s^t
ps, tq ÞÑ covpUj psq, U` ptqq “ ´
0

∆A` puq
dAj puq “: σj` ps ^ tq.
H̄puq

See Appendix C for a derivation of this asymptotic covariance function. In order to account for
this dependence structure in a joint convergence consideration, the wild bootstrap of the previous
section needs to be adjusted once more. Therefore, let ξj`i , j, ` “ 1, . . . , k, i “ 1, . . . , n, be i.i.d.
2
random variablesřwith Epξ111 q “ 0 and Epξ111
q “ 1, which are also independent of the data.
k
Denote by N “ j“1 Nj the number of all kinds of events. Then, the single components of the
wild bootstrap version of the multivariate Nelson-Aalen estimator pW1 , . . . , Wk q are given by
n
? ÿ
x
Wj ptq “ n ξjji

żtd

Y puq ´ ∆N puq dNji puq
Y puq
Y puq
0
i“1
d
żt
żtd
k
n
n
”? ÿ
1 ÿ
∆Nj puq dN`i puq ? ÿ
∆N` puq dNji puq ı
`?
signp` ´ jq n ξj`i
` n ξ`ji
,
Y puq Y puq
Y puq Y puq
2 `“1
0
0
i“1
i“1
7

where signpxq “ 1tx ą 0u ´ 1tx ă 0u is the signum function. This signum function is important
in order to insure the required negative covariance between all components. The following large
sample properties hold:
Theorem 3.2. Assume Condition 2.1. Given F0 and as n Ñ 8, we have the following conditional
weak convergence on the product càdlàg function space Dk r0, Ks, equipped with the max-supnorm:
d
x1 , W
x2 , . . . , W
xk q ÝÑ
pW
pU1 , U2 , . . . , Uk q

in outer probability,

where pU1 , U2 , . . . , Uk q is the same Gaussian martingale as in Theorem 3.1.
The proof is given in Appendix D. Note that, if we are interested in just a single univariate
Nelson-Aalen estimator, the present approach yields the same limit distribution as the wild bootstrap technique proposed in Section 2. Hence, it does – asymptotically – not matter which of both
techniques is applied to the univariate Nelson-Aalen estimator.
Variance and covariance estimators (also for the wild bootstrap versions) are again motivated
by the predictable and optional covariation processes of the involved martingales. The resulting
estimators turn out to be same as those obtained by the plug-in method:
żt

Y puq ´ ∆Nj puq
dNj puq,
j “ 1, . . . , k,
Y 3 puq
0
żt
∆Nj puq
pj` ptq “ ´n
σ
dN` puq,
j, ` “ 1, . . . , k; j ‰ `,
3
0 Y puq
pj2 ptq
σ

“n

are the usual Greenwood-type (co)variance estimators and
n
ÿ

ż¨

Y puq ´ ∆Nj puq
dNji puq
Y 3 puq
0
i“1
ż
ż¨
n
n
ı
ÿ
n ÿ ” ÿ 2 ¨ ∆Nj puq
∆N` puq
2
ξj`i
dN
puq
`
ξ
dN
puq
,
`
`i
ji
`ji
3
3
2 `‰j i“1
0 Y puq
0 Y puq
i“1
ż
ż¨
n
n
ı
ÿ
∆N` puq
n ” ÿ 2 ¨ ∆Nj puq
2
p
pj` p¨q “ ´
σ
ξj`i
dN
puq
`
ξ
dN
puq
,
`i
ji
`ji
3
3
2 i“1
0 Y puq
0 Y puq
i“1
p
pj2 p¨q “ n
σ

2
ξjji

j “ 1, . . . , k,
j, ` “ 1, . . . , k; j ‰ `,

are the optional process-type (co)variance estimators motivated from the wild bootstrap martingale
properties. Assume that all ξ111 have finite fourth moments. By applications of Glivenko-Cantelli
theorems in combination with the continuous mapping theorem, the Greenwood-type (co)variance
pj and σ
pj` are shown to be uniformly consistent for σj2 and σj` , respectively. For the wild
estimators σ
8

bootstrap-type (co)variance estimators, we can parallel the arguments in the proof of Theorem 2.2,
4
after first assuming the existence of fourth moments Epξ111
q ă 8: In points of continuity of
all cumulative hazard functions, i.e. on r0, KszD, Rebolledo’s martingale central limit theorem
applies and it also implies the uniform consistency of the optional variation process increments.
p
pj2 and
pj2 by σ
In points of discontinuity, which are finitely many by assumption, we approximate σ
apply the conditional Chebyshev inequality (given F0 ) in order to show the negligibility of the
p
pj2 in probability. The last argument can be repeated for the covariance estimators.
pj2 ´ σ
differences σ
A final application of the continuous mapping theorem yields
p
p
pj2 , σ
pj2 ÝÑ σj
σ

p
p
pj` , σ
pj` ÝÑ σj`
σ

and

uniformly on r0, Ks in (conditional outer) probability,

for all j ‰ ` as n Ñ 8.

4

Extension to Aalen-Johansen Estimators for Cumulative Incidence Functions in Competing Risks

Denote the cumulative incidence functions by Fj ptq, j “ 1, . . . , k, which specify the probabilities
to die due to cause j during the time interval r0, ts. For ease of presentation, we consider the
situation of k “ 2 competing risks which is achieved by aggregating all but the first risk to be
the second competing risk. The general results are obtained by replacing U2 , A2 , and F2 in the
representations below by U ´ U1 , A ´ A1 , and 1 ´ S ´ F1 , respectively. Here, A denotes the allcause cumulative hazard function. Utilizing the functional delta-method in combination with the
weak convergence results for the Nelson-Aalen estimator,
we as well obtain a weak convergence
şt
p
p
p1 puq for the cumulative incidence
theorem for the Aalen-Johansen estimator F1 ptq “ 0 Spu´qdA
şt
function F1 ptq “ 0 Spu´qdA1 puq:
Theorem 4.1. As n Ñ 8,
WF1

?
d
“ npFp1 ´ F1 q ÝÑ UF1 “

ż¨
0

1 ´ F2 pu´q ´ F1 p¨q
dU1 puq `
1 ´ ∆Apuq

ż¨
0

F1 pu´q ´ F1 p¨q
dU2 puq,
1 ´ ∆Apuq

where UF1 is a zero-mean Gaussian process with covariance function
ż s^t
σF2 1

: ps, tq ÞÑ

p1 ´ F2 puq ´ F1 psqqp1 ´ F2 puq ´ F1 ptqq 1 ´ ∆A1 puq
dA1 puq
p1 ´ ∆Apuqq2
H̄puq
0
ż s^t
pF1 puq ´ F1 psqqpF1 puq ´ F1 ptqq 1 ´ ∆A2 puq
`
dA2 puq
p1 ´ ∆Apuqq2
H̄puq
0
9

ż s^t
´
0

ż s^t
´
0

p1 ´ F2 puq ´ F1 psqqpF1 puq ´ F1 ptqq ∆A1 puq
dA2 puq
p1 ´ ∆Apuqq2
H̄puq
p1 ´ F2 puq ´ F1 ptqqpF1 puq ´ F1 psqq ∆A2 puq
dA1 puq.
p1 ´ ∆Apuqq2
H̄puq

For the application of the functional delta-method, note that the Aalen-Johansen estimator in
the present competing risks framework is a combination of the Wilcoxon and the product integral functional applied to the multivariate Nelson-Aalen estimator. Both of these functionals are
Hadamard-differentiable as shown for example in Section 3.9 of van der Vaart and Wellner (1996).
A derivation of the above asymptotic covariance function is presented in Appendix E.
?
Now, an appropriate wild bootstrap version of npFp1 ´ F1 q is given by
ż¨
xF1 ptq “
W
0

ż¨ p
1 ´ Fp2 pu´q ´ Fp1 p¨q x
F1 pu´q ´ Fp1 p¨q x
dW1 puq `
dW2 puq,
p
p
1 ´ ∆Apuq
1 ´ ∆Apuq
0

x1 and W
x2 are again the wild bootstrap versions of the Nelson-Aalen estimator as presented
where W
in Section 3. Using similar martingale arguments as in Appendix B, we obtain the following
conditional central limit theorem for the wild bootstrap version of the Aalen-Johansen estimator:
Theorem 4.2. Assume Condition 2.1. Given F0 and as n Ñ 8, we have the following weak
convergence on the càdlàg function space Dr0, Ks, equipped with the sup-norm:
d
xF1 ÝÑ
W
UF1

in probability,

where UF1 is the same Gaussian process as in Theorem 4.1.
Remark 4.3 (The weird bootstrap). Note that the very same proofs may be applied to verify that
the above conditional central limit theorems hold for the weird bootstrap as well. This resampling scheme corresponds to choosing ξj`i ` 1 „ BinpY pXi q, maxp1, Y pXi qq´1 q, where Xi is the
censoring or event time of individual i, whichever comes first. This is a particular choice of the
data-dependent multiplier bootstrap of Dobler et al. (2015). In their article, heuristic arguments
for the second order correctness under absolute continuity of the data have shown that centered
unit Poisson variates and weird bootstrap multipliers perform favorably in comparison to standard
normal wild bootstrap weights. In order to also check the preference of either of the first two resampling procedures in the present set-up, where ties are allowed, we included the weird bootstrap
yielding competing inference methods into the subsequent simulation study.
10

Estimators for σF2 1 and its wild bootstrap variant are obtained similarly as such estimators for
the Nelson-Aalen (co)variances, i.e. via plug-in:
ż s^t
pF2 1
σ

: ps, tq ÞÑ
0

p1 ´ Fp2 pu´q ´ Fp1 psqqp1 ´ Fp2 pu´q ´ Fp1 ptqq 2
dp
σ1 puq
2
p
p1 ´ ∆Apuqq

ż s^t p
pF1 pu´q ´ Fp1 psqqpFp1 pu´q ´ Fp1 ptqq 2
`
dp
σ2 puq
2
p
p1 ´ ∆Apuqq
0
ż s^t
p1 ´ Fp2 pu´q ´ Fp1 psqqpFp1 pu´q ´ Fp1 ptqq
`
dp
σ12 puq
2
p
p1 ´ ∆Apuqq
0
ż s^t
p1 ´ Fp2 pu´q ´ Fp1 ptqqpFp1 pu´q ´ Fp1 psqq
σ21 puq.
dp
`
2
p
p1 ´ ∆Apuqq
0
p
pF2 1 is obtained by replacing all estimators σ
pj and σ
pj` , j ‰ `, by their wild bootstrap
Similarly, σ
p
p
pj` , respectively. Their uniform (conditional) consistencies for σF2 1 on r0, Ks
pj and σ
counterparts σ
follow immediately by the uniform consistency of the Nelson-Aalen (co)variance estimators and
the continuous mapping theorem.
Remark 4.4 (Deduced confidence bands). Following the lines of Beyersmann et al. (2013), timesimultaneous confidence bands for F1 can be deduced. In particular, let φpsq “ logp´ logp1 ´ sqq
be a transformation applied to F1 in order to ensure band boundaries between 0 and 1 and let
g1 psq “ logp1 ´ Fp1 psqq{p
ρpsq and g2 psq “ logp1 ´ Fp1 psqq{p1 ` ρp2 psqq be weight functions leading
to the usual equal precision and Hall-Wellner bands, respectively; see Andersen et al. (1993).
pF2 1 psq{p1 ´ Fp1 psqq2 . Let gp1 and gp2 be their wild bootstrap counterparts, i.e. the
Here ρp2 psq “ σ
p
pF2 1 are replaced by σ
pF2 1 . The confidence bands for F1 are then derived from
variance estimates σ
?
the asymptotics of the supremum distance Z1j “ supuPrt1 ,t2 s | ngj puqpφpFp1 puqq ´ φpF1 puqqq| and
xF1 puq|, where rt1 , t2 s Ă r0, Ks
its wild bootstrap counterpart Zp1j “ supuPrt ,t s |p
gj puqφ1 pFp1 puqqW
1 2

and j P t1, 2u. Let q0.95,j be the conditional 95% quantile of Zp1j given the data. The resulting
´1{2
asymptotic 95% confidence bands are 1 ´ p1 ´ Fp1 psqqexpp˘n q0.95,j {gj psqq , s P rt1 , t2 s, j “ 1, 2.

5

Small Sample Behaviour

We empirically assess the difference between the common wild bootstrap approach and the adjusted wild bootstrap proposed in this article via simulation studies. We simulated the wild bootstrap procedures based on standard normal and centered unit Poisson multipliers as well as the
11

weird bootstrap of Remark 4.3. These methods are compared in terms of the simulated coverage probabilities of the confidence bands described in Remark 4.4. We consider a simulation set-up motivated by Dobler and Pauly (2014), i.e. we chose the cause-specific hazard rates
α1 ptq “ expp´tq and α2 ptq “ 1 ´ expp´tq which yield the cumulative function of the first risk
F1 ptq “ 0.5p1 ´ expp´2tqq. In order to allow for tied data, we pre-specify different discretization lattices and round different proportions of the population to the nearest discretization point.
In particular, we choose the discretization lattices to be t0, k1 , k2 , . . . u, where k P t5, 10, 20u, and
the discretization probabilities to be p P t0, 0.25, 0.5, 0.75, 1u. The resulting theoretic cumulative
incidence functions
F1p,k ptq “ pF1

´ rkt ´ 0.5s
k

`

0.5 ¯
` p1 ´ pqF1 ptq
k

are presented in Figure 1, with the exception of the continuous function F1 . Here rss denotes
the integer closest to s P R. For simulating data, which have the desired cumulative incidence
function F1p,k , it is mandatory to first round the event times Ti , and then generate the event types εi
in a second step, according to the formula
ˇ rT ks
¯ F pu ` 1 q ´ F pmaxpu ´ 1 , 0qq
1
1
ˇ i
2k
2k
“u “
,
P εi “ 1ˇ
1
1
k
Spmaxpu ´ 2k , 0qq ´ Spu ` 2k
q
´

where S : t ÞÑ expp´tq denotes the survival function of the continuous random variables Ti .
Censoring is introduced by i.i.d. standard exponentially distributed random variables. If the
ith survival time is discretized, then we discretize the ith censoring time as well. Finally, we take
the minimum out of each such pair and mark an individual as censored whenever the (discretized)
censoring time precedes the (discretized) survival time. The sample size increases from n “ 50
to n “ 250 in steps of 25. We choose the time interval, along which asymptotic 95% confidence
bands shall be constructed, to be r0.25, 0.75s. The simulations have been conducted using R version
3.2.3 (R Development Core Team, 2016) using 10,000 outer Monte Carlo iterations and 999 wild
bootstrap replicates.
Tables 1 to 6 contain the simulated coverage probabilities of equal precision and Hall-Wellner
bands for simulation set-ups with a discrete component in the cumulative incidence function, i.e.
p ą 0. The columns of simulation results corresponding to the common wild / weird bootstrap
procedures are entitled old, whereas the columns showing the results of the respective adjusted
wild / weird bootstrap are entitled new.
At first, we start with a discussion on the choice of multipliers. For equal precision bands and
in almost all set-ups, there is a pronounced superiority of the wild bootstrap with centered unit
Poisson multipliers and the weird bootstrap over the respective coverage probabilities of the bands
based on standard normal weights. This is true for the common resampling procedures as well
as for the proposed adjusted bootstraps. For Hall-Wellner bands, this superiority is not as much
pronounced and sometimes even the confidence bands based on standard normal multipliers yield
12

the most accurate coverage probabilities. But in cases, where this is so, the deviance is only very
small. The phenomenon, that standard normal multipliers yield a worse performance than those
with skewness equal to one, is in line with the findings in a revised version of Dobler et al. (2015)
where also heuristic theoretic arguments for a second-order correctness of both superior resampling
procedures are provided. As there is, all in all, not much of a difference between the simulated
coverage probabilities of the centered unit Poisson wild bootstrap and the weird bootstrap, we only
focus on the results of the Poisson choice. In general, the equal precision bands are more accurate
than the Hall-Wellner bands.
Here, the discretization-adjusted wild bootstrap almost always yielded coverage probabilities
closer to the nominal level in comparison to the unadjusted wild bootstrap. The deviances between
these coverage probabilities of each of those two resampling procedures appears to be larger the
higher the discretization probability p and the coarser the discretization lattice is. For instance, this
difference even amounts to 4.03 percentage points in case of the Hall-Wellner bands, n “ 50, and
p “ 1 and to 3.83 percentage points in case of equal precision bands and the same p and n.
In case of k P t5, 10u, the coverage probabilities of the common wild bootstrap do not appear
to converge at all towards 95% as the sample size increases. Instead, the simulated probabilities
fluctuate around 93% or even 92%. On the other hand, the discretization-adjusted equal precision
wild bootstrap bands yield much better coverage probabilities which are greater than 94% or at
least in the high 93%-region for larger sample sizes. In contrast to the unadjusted procedure, we
observe for small samples and for the adjusted confidence bands coverage probabilities closer to
the nominal level for higher discretization probabilities p. This is only reasonable, as p “ 100%
corresponds to a multivariate, but not an infinite-dimensional statistical problem. We do not see
this tendency for the unadjusted procedure in case of k P t5, 10u, which again stresses that it is not
suitable for these kinds of tied data regimes.
Finally, Table 7 shows the corresponding results for the scenario in which the continuous F1
is the true cumulative incidence function and the usual wild bootstrap technique yields asymptotically exact inference procedures. However, it is surprising to see that even here the adjusted
wild bootstrap yields more accurate confidence bands in comparison to the unadjusted procedure.
Therefore, there is apparently no loss at all in utilizing the discretization adjustment – with or
without ties in the data.
All in all, we conclude that the proposed discontinuity adjustment should always be applied
in order to greatly improve the coverage probabilities of confidence bands for F1p,k . The present
simulation results show this improvement, which amounts to up to two or three percentage points
for smaller samples, in many conducted simulation scenarios. As the standard normal variatebased wild bootstrap disappoints in general, our final advice is to combine the present discontinuity
adjustment with the wild bootstrap based on the Poisson-distributed random variables or the weird
bootstrap. Additionally, equal precision bands should be preferred to Hall-Wellner confidence
bands due to the slight but frequent difference in coverage probabilities.

13

Figure 1: Cumulative incidents functions F1p,k underlying the present simulations.
14

set-ups
p
n
50
75
100
125
25% 150
175
200
225
250
50
75
100
125
50% 150
175
200
225
250
50
75
100
125
75% 150
175
200
225
250
50
75
100
125
100% 150
175
200
225
250

N p0, 1q
old
new
87.31 88.47
89.48 90.39
90.64 91.26
91.21 91.96
91.12 91.66
91.69 92.14
91.99 92.46
92.51 92.98
92.48 92.88
87.41 89.31
89.35 90.99
90.01 91.53
90.56 91.92
90.59 91.93
91.21 92.22
91.22 92.28
91.44 92.61
91.72 92.82
87.23 90.60
90.09 92.56
90.58 92.86
90.83 93.03
90.84 93.11
90.83 93.13
90.91 93.11
91.15 93.43
91.76 93.76
87.71 91.98
88.79 92.70
89.89 93.25
89.76 93.14
89.64 93.21
90.10 93.24
90.27 93.55
90.01 93.07
89.71 93.18

P oip1q ´ 1
old
new
89.44 90.61
91.55 92.24
92.32 92.96
92.36 93.07
92.48 92.98
93.10 93.21
93.31 93.80
93.29 93.67
93.65 93.91
89.09 91.28
91.28 92.71
92.39 93.64
92.95 93.90
92.44 93.38
92.58 93.49
92.30 93.31
92.70 93.73
93.13 94.27
89.12 91.56
90.96 93.22
92.13 94.23
92.13 94.26
91.96 93.94
91.66 93.87
92.45 94.46
92.52 94.49
92.38 94.36
89.08 92.91
90.26 93.74
90.90 94.06
90.92 93.88
90.84 93.97
90.52 93.52
91.06 94.00
91.30 94.19
90.92 93.89

weird
old
new
88.93 90.36
91.13 92.01
92.21 92.64
93.43 93.76
93.37 93.75
93.20 93.63
93.63 94.05
93.58 94.03
93.41 93.84
89.76 91.29
91.03 92.59
91.88 93.29
92.29 93.43
93.04 94.05
92.66 93.73
93.16 94.19
93.25 94.16
93.07 94.05
89.47 92.11
91.29 93.77
91.59 93.83
91.75 94.08
92.29 94.35
92.25 94.03
92.50 94.49
92.22 94.12
92.73 94.53
89.56 93.19
90.21 93.64
90.47 93.84
90.44 93.64
91.00 93.87
91.25 93.99
91.09 94.12
91.08 94.26
90.76 93.86

Table 1: Simulated coverage probabilities of equal precision bands in per cent where k “ 5

15

set-ups
p
n
50
75
100
125
25% 150
175
200
225
250
50
75
100
125
50% 150
175
200
225
250
50
75
100
125
75% 150
175
200
225
250
50
75
100
125
100% 150
175
200
225
250

N p0, 1q
old
new
87.72 88.64
89.22 90.16
90.68 91.22
91.13 91.61
91.63 92.16
91.79 92.06
91.93 92.37
92.40 92.77
92.17 92.61
87.15 88.59
89.78 90.80
90.24 91.32
91.17 92.11
91.63 92.52
92.65 93.29
92.21 92.87
92.20 92.98
92.96 93.52
87.63 89.82
89.52 91.46
90.30 91.90
91.44 92.88
91.49 92.76
91.50 92.76
92.00 93.08
92.20 93.12
92.49 93.78
88.28 91.10
89.81 92.2
90.54 92.66
91.12 93.19
91.41 93.44
91.60 93.59
91.49 93.46
92.32 93.98
91.81 93.81

P oip1q ´ 1
old
new
89.37 90.69
91.86 92.55
92.43 93.07
93.05 93.50
93.5 93.77
93.87 94.19
93.84 93.94
93.54 93.89
94.03 94.26
89.20 90.73
91.46 92.47
92.58 93.25
92.99 93.70
93.25 93.98
92.67 93.34
93.36 94.03
94.07 94.80
93.71 94.23
89.19 91.41
91.97 93.42
92.45 93.76
92.93 94.03
93.20 94.50
92.98 94.18
93.92 94.94
93.57 94.56
93.51 94.55
90.21 92.75
91.70 93.90
92.10 93.98
92.58 94.36
92.83 94.61
92.84 94.57
93.03 94.76
92.92 94.70
92.79 94.58

weird
old
new
89.87 90.84
91.75 92.35
93.06 93.31
93.42 93.79
93.48 93.79
93.75 93.93
93.70 93.88
93.73 93.97
94.24 94.56
89.18 90.50
91.18 92.22
92.53 93.47
93.06 93.76
93.17 93.95
93.06 93.71
93.51 94.23
93.09 93.96
93.37 93.93
88.95 90.91
91.36 92.96
92.52 93.95
92.83 94.14
93.22 94.23
93.66 94.63
93.75 94.85
93.17 94.13
93.69 94.81
90.31 92.68
91.77 93.57
92.09 93.96
92.74 94.47
92.77 94.77
92.84 94.79
93.25 94.97
93.38 94.95
92.98 94.73

Table 2: Simulated coverage probabilities of equal precision bands in per cent where k “ 10

16

set-ups
p
n
50
75
100
125
25% 150
175
200
225
250
50
75
100
125
50% 150
175
200
225
250
50
75
100
125
75% 150
175
200
225
250
50
75
100
125
100% 150
175
200
225
250

N p0, 1q
old
new
87.21 88.03
89.75 90.33
91.01 91.61
92.09 92.29
91.87 92.39
92.44 92.71
92.76 93.03
93.23 93.44
92.84 93.11
88.02 89.22
89.88 90.74
90.95 91.69
91.98 92.46
91.90 92.54
92.95 93.40
92.79 93.11
92.78 93.25
93.01 93.28
88.42 89.91
90.28 91.50
91.43 92.45
91.41 92.30
92.73 93.31
92.44 93.14
93.04 93.84
93.32 93.94
93.03 93.59
88.87 90.73
90.95 92.58
91.98 93.28
92.91 93.94
92.82 93.88
93.16 94.25
93.53 94.51
94.06 94.76
93.12 94.31

P oip1q ´ 1
old
new
89.93 90.94
92.18 92.50
93.12 93.52
93.99 94.18
93.63 93.78
94.03 94.40
94.18 94.35
94.46 94.77
93.88 94.10
89.57 90.86
92.37 93.28
93.19 93.74
93.44 93.83
93.60 94.01
94.00 94.42
94.05 94.39
94.16 94.53
93.92 94.30
90.81 92.12
92.73 93.94
93.30 94.21
93.84 94.58
94.01 94.77
94.44 95.04
94.13 94.65
94.41 94.95
94.27 94.85
91.04 92.79
93.21 94.28
93.74 95.01
93.69 94.88
94.18 95.10
94.88 95.59
94.38 95.30
94.63 95.59
94.95 95.79

weird
old
new
90.15 90.69
91.78 92.38
92.71 93.08
93.34 93.67
93.69 93.85
94.19 94.35
93.76 94.07
94.02 94.23
94.16 94.41
89.59 90.81
92.34 93.00
93.02 93.71
93.91 94.2
93.97 94.27
94.11 94.45
93.68 94.00
94.28 94.67
93.77 94.03
90.51 91.96
92.85 93.67
93.02 93.85
93.57 94.08
93.54 94.25
94.09 94.82
94.02 94.58
94.08 94.69
94.29 94.77
91.16 92.69
93.02 94.24
93.74 94.67
94.13 95.13
94.58 95.49
94.85 95.74
94.59 95.29
94.69 95.39
94.79 95.60

Table 3: Simulated coverage probabilities of equal precision bands in per cent where k “ 20

17

set-ups
p
n
50
75
100
125
25% 150
175
200
225
250
50
75
100
125
50% 150
175
200
225
250
50
75
100
125
75% 150
175
200
225
250
50
75
100
125
100% 150
175
200
225
250

N p0, 1q
old
new
88.65 89.93
89.82 90.60
90.51 91.11
90.92 91.69
91.41 91.70
91.79 92.17
92.24 92.67
92.22 92.70
92.38 93.03
87.72 90.36
89.17 91.06
89.86 91.62
90.86 92.28
90.99 92.28
91.41 92.67
91.77 92.99
91.70 92.87
91.44 92.84
87.35 90.83
89.97 92.41
89.89 92.60
90.52 92.98
90.61 93.03
91.03 93.41
91.29 93.53
91.63 93.76
91.44 93.55
88.99 92.88
90.32 93.74
90.59 94.04
91.27 94.59
91.63 94.59
91.53 94.56
91.78 94.65
91.94 94.92
92.40 94.92

P oip1q ´ 1
old
new
89.33 90.79
90.59 91.54
90.93 91.82
91.11 91.82
91.53 92.20
92.75 93.17
92.18 92.59
92.75 93.20
92.53 93.05
88.52 90.87
89.54 91.43
90.26 91.87
91.27 92.69
91.34 92.77
91.76 92.84
92.20 93.57
91.72 92.95
92.00 93.20
88.74 92.14
90.57 93.39
90.32 93.29
90.95 93.38
91.15 93.80
91.54 93.95
90.98 93.49
91.70 93.52
91.55 93.78
89.60 93.67
90.71 94.31
91.17 94.37
91.79 94.81
91.42 94.50
91.95 94.80
91.82 94.68
91.87 94.81
91.86 94.76

weird
old
new
88.85 90.18
90.39 91.32
91.13 91.88
91.07 91.92
92.21 92.74
92.65 93.16
92.26 92.73
92.64 93.14
92.44 92.68
88.77 91.14
89.92 91.87
90.60 92.18
91.00 92.39
91.15 92.57
91.12 92.42
91.78 92.86
91.94 93.14
91.65 93.04
88.32 91.82
90.09 93.16
90.19 93.01
91.03 93.55
91.16 93.67
90.88 93.37
91.30 93.41
91.11 93.61
91.95 94.24
89.68 93.41
90.23 93.99
90.60 94.14
91.56 94.55
91.20 94.60
90.97 94.22
91.76 94.68
91.59 94.56
91.35 94.46

Table 4: Simulated coverage probabilities of Hall-Wellner bands in per cent where k “ 5

18

set-ups
p
n
50
75
100
125
25% 150
175
200
225
250
50
75
100
125
50% 150
175
200
225
250
50
75
100
125
75% 150
175
200
225
250
50
75
100
125
100% 150
175
200
225
250

N p0, 1q
old
new
87.21 88.70
89.23 90.08
90.54 90.86
91.49 91.96
92.52 92.80
92.61 92.88
92.25 92.76
92.92 93.41
92.37 92.78
87.32 89.05
89.94 91.20
91.24 92.37
91.05 92.16
91.99 92.82
92.37 93.10
92.24 92.79
92.29 93.22
92.05 92.73
88.15 90.60
90.41 92.20
91.25 92.74
91.76 93.21
91.67 92.85
92.30 93.57
92.25 93.46
92.40 93.63
92.70 93.79
88.99 92.00
91.10 93.57
90.84 93.12
92.06 93.80
91.95 93.89
92.41 94.23
92.34 94.21
92.73 94.31
93.32 94.90

P oip1q ´ 1
old
new
89.27 90.49
90.91 91.57
91.86 92.55
92.45 92.86
92.31 92.94
92.86 93.06
92.66 92.88
93.02 93.32
92.90 93.29
88.85 90.49
91.08 92.27
91.08 92.32
91.61 92.49
91.92 92.79
92.53 93.21
92.42 93.37
92.75 93.44
92.70 93.38
88.32 91.01
91.34 92.98
90.97 92.66
91.99 93.34
91.73 93.04
92.50 93.75
92.11 93.50
92.43 93.88
92.38 93.66
89.50 92.73
91.14 93.59
92.10 94.27
92.06 94.18
92.56 94.65
92.87 94.70
93.11 94.88
93.25 94.85
92.74 94.41

weird
old
new
89.00 90.36
90.65 91.51
91.40 91.89
91.59 92.25
92.18 92.65
92.71 93.23
92.80 93.20
93.00 93.38
93.19 93.46
89.17 90.77
90.65 91.81
91.51 92.59
91.88 92.77
92.00 93.01
92.66 93.38
92.59 93.25
92.87 93.64
92.05 92.74
89.49 91.71
90.39 92.28
91.19 93.00
91.57 93.09
91.90 93.31
92.08 93.34
92.05 93.46
92.89 94.09
92.56 93.67
89.59 92.83
90.91 93.76
91.66 93.77
92.25 94.27
92.60 94.59
92.27 94.35
92.99 94.88
93.06 94.89
92.90 94.61

Table 5: Simulated coverage probabilities of Hall-Wellner bands in per cent where k “ 10

19

set-ups
p
n
50
75
100
125
25% 150
175
200
225
250
50
75
100
125
50% 150
175
200
225
250
50
75
100
125
75% 150
175
200
225
250
50
75
100
125
100% 150
175
200
225
250

N p0, 1q
old
new
88.39 89.53
90.39 90.99
91.39 91.89
91.69 92.17
92.07 92.51
92.24 92.55
93.01 93.27
93.11 93.32
92.88 93.03
88.77 90.02
90.43 91.37
91.74 92.45
92.12 92.82
92.12 92.59
92.43 93.00
92.94 93.22
92.67 93.01
93.41 93.58
89.17 90.83
90.99 92.23
91.15 92.12
92.21 93.15
92.64 93.39
92.50 93.21
92.73 93.39
92.91 93.51
93.34 94.01
90.11 92.25
91.22 92.70
92.01 93.62
92.39 93.45
92.77 94.10
93.22 94.35
93.88 94.73
93.74 94.62
93.61 94.56

P oip1q ´ 1
old
new
89.93 91.16
90.67 91.39
91.75 92.26
92.54 93.03
92.57 92.83
92.76 92.97
92.58 93.04
92.99 93.25
93.60 93.69
89.28 90.66
91.27 92.29
91.99 92.74
92.49 93.11
93.05 93.60
92.96 93.25
93.08 93.50
93.20 93.92
93.03 93.47
89.70 91.53
91.29 92.41
92.36 93.24
92.57 93.44
93.08 93.98
93.25 93.91
93.34 94.15
93.18 93.96
93.84 94.39
90.74 92.78
92.57 94.08
93.09 94.37
93.28 94.64
93.47 94.59
93.66 94.84
93.79 94.78
94.25 95.07
93.94 94.95

weird
old
new
89.14 90.05
90.74 91.52
91.92 92.30
92.52 92.90
92.92 93.28
92.8 93.07
92.96 93.21
93.16 93.42
93.05 93.31
89.53 90.89
90.87 91.92
91.94 92.62
92.08 92.61
92.56 93.18
93.00 93.45
92.69 93.28
93.51 93.90
92.74 93.25
89.65 91.36
91.21 92.32
92.26 93.37
92.49 93.50
92.68 93.46
92.67 93.44
93.10 93.78
93.61 94.11
92.92 93.67
90.52 92.90
91.94 93.44
92.49 93.84
93.59 94.56
93.28 94.41
94.00 94.93
93.90 94.84
93.69 94.70
94.17 95.12

Table 6: Simulated coverage probabilities of Hall-Wellner bands in per cent where k “ 20

20

n
50
75
100
125
150
175
200
225
250

N p0, 1q
old
new
87.61 88.64
90.15 90.75
90.78 91.31
91.72 92.07
92.29 92.54
92.27 92.53
92.67 92.85
92.98 93.05
92.78 92.99

Hall-Wellner
P oip1q ´ 1
old
new
88.80 90.03
90.71 91.46
91.13 91.58
92.02 92.31
92.60 92.85
92.48 92.83
92.89 93.02
92.60 92.71
92.93 93.14

weird
old
new
88.72 89.75
90.74 91.48
91.77 92.23
91.58 92.23
92.90 93.21
92.86 93.17
93.18 93.34
93.26 93.68
93.14 93.32

N p0, 1q
old
new
87.08 88.06
89.36 90.21
90.74 91.24
91.35 91.77
91.33 91.52
92.45 92.64
92.87 93.04
92.69 92.93
93.20 93.43

Equal Precision
P oip1q ´ 1
old
new
89.55 90.36
91.79 92.33
93.20 93.43
92.96 93.37
93.64 93.85
93.91 94.14
93.68 93.80
93.85 94.12
93.68 93.81

weird
old
new
89.43 90.48
91.62 92.14
92.43 92.72
93.19 93.54
93.51 93.76
94.06 94.39
93.47 93.83
94.16 94.22
94.25 94.30

Table 7: Simulated coverage probabilities of confidence bands for F1 in per cent where p “ 0

6

Real Data Example

We applied the present discretization adjustment to the sir.adm data-set of the R package mvna.
It consists of competing risks data of patients who are in an intensive care unit (ICU), where
the event of primary interest, “alive discharge out of ICU”, competes against the secondary event
“death in ICU”. For seeing the difference between the common and the new approach more clearly,
we analyzed the subset of all male patients suffering from pneumonia. Out of these n “ 63
individuals, five have been right-censored and 41 out of all 44 type 1 events fell into the time
interval r5, 55s, which we chose for confidence band construction. Due to the worse performance
of the wild bootstrap based on standard normal multipliers as seen in Section 5, we derived these
bands by only using centered unit Poisson variates. In order to minimize the computational error in
the quantile-finding process, 99,999 wild bootstrap iterations have been conducted. The confidence
bands resulting from the weird bootstrap almost coincide with those just described. Therefore, they
are not shown.
The resulting confidence bands are presented in Figure 2. For both kinds of bands, (equal
precision bands in the upper panel, Hall-Wellner bands in the lower panel), we see that the discretization adjustment leads to a widening in comparison to the unadjusted bands. This is in line
with the results of the simulation study in Section 5, where the unadjusted bands appeared to be the
most liberal, i.e. the smallest. In particular, the adjusted equal precision bands cover an additional
area of 2.1 percentage points at the terminal time point t “ 55, whereas this deviance even amounts
to 3.3 percentage points for the Hall-Wellner bands. This might not appear to be much at a first
glance at the plots in Figure 2. But in fact, it may be the cause for a formidable improvement of the
bands’ coverage probability: The simulation results of Section 5 for k “ 20, discretization probability p “ 100%, and sample sizes n P t50, 75u suggest that the adjusted wild bootstrap procedure
21

might improve the coverage probabilities of both kinds of bands by approximately two percentage
points. With a view towards the liberal behaviour of the unadjusted bands, these enhancements of
the coverage probabilities are highly worthwhile.

Figure 2: Asymptotic 95% equal precision (upper) and Hall-Wellner bands (lower panel) for the
cumulative incidence function of the competing risk “alive discharge out of ICU” for male patients
suffering from pneumonia.
22

7

Discussion and Future Research

In this article, we analyzed a discontinuity adjustment of the common wild bootstrap applied to
right-censored competing risks data. This adjustment is absolutely necessary, as ties in the data introduce an asymptotic dependence between multiple cause-specific Nelson-Aalen estimators. The
common wild bootstrap fails in reproducing this effect since it establishes independence for all
sample sizes. The problem is even more involved for Aalen-Johansen estimators of cumulative
incidence functions, which are non-linear functionals of all cause-specific hazards. Simulation
results reported the striking liberality of the unadjusted bands which also fail to keep the nominal level asymptotically. Instead, the discretization-adjusted wild bootstrap greatly improves the
coverage probability. This effect is more pronounced the more discrete the event times are. But
even in the absolutely continuous case, the suggested procedure appears to perform preferably.
Therefore, we advise to always use the adjustment when right-censored competing risks data shall
be analyzed. The real data example reveals that the discontinuity adjustment does only lead to
slight widening of the common wild bootstrap-based bands which is already enough to improve
the coverage accuracy greatly.
The presented wild bootstrap approach may be extended to more general Markovian multi-state
models since the martingale arguments of Appendix B still apply. A still open question is, whether
the common wild bootstrap also fails in case of tied survival data which are assumed to follow the
Cox proportional hazards model (Cox, 1972). See Lin et al. (1993) for the wild bootstrap applied
to absolutely continuous survival data following the Cox model. And if it fails, does the method
proposed in this article require further modification?

Acknowledgements
The author likes to thank Markus Pauly for helpful discussions. He also appreciates the support by
the German Research Foundation (DFG).

A

Detailed Derivation of the Asymptotic Variance of the NelsonAalen Estimator

Define H uc ptq “

řk

ENj1 ptq as the probability that an uncensored event due to any cause occurs
? p
d
´ Aq ÝÑ
until
time t. According to van der Vaart and Wellner
(1996), p. 383f, we have npA
ş¨ uc
ş
t
M pduq{H̄puq, where M uc ptq “ Guc ptq ´ 0 ḠpuqdApuq is a zero-mean Gaussian martingale.
0
j“1

23

Its variance function is determined by
EGuc psqGuc ptq “ H uc ps ^ tq ´ H uc psqH uc ptq,
E ḠpsqḠptq “ H̄ps _ tq ´ H̄psqH̄ptq,
EGuc psqḠptq “ pH uc psq ´ H uc pt´qq1tt ď su ´ H uc psqH̄ptq.
şt
Note that Aptq “ 0 H uc pduq{H̄puq. Thus, for s ď t, the covariance function of M uc at ps, tq is
EpM uc psqM uc ptqq
żsżt

H̄pu _ vq ´ H̄puqH̄pvq
dH uc pvqdH uc puq
“ H psq ´ H psqH ptq `
H̄puq
H̄pvq
0 0
żt
dH uc puq
´ rpH uc psq ´ H uc pu´qq1tu ď su ´ H uc psqH̄puqs
H̄puq
ż0s
dH uc pvq
´ rpH uc ptq ´ H uc pv´qq1tv ď tu ´ H uc ptqH̄pvqs
H̄pvq
0
żsżu”
żsżt”
ı
ı
1
1
uc
uc
uc
uc
uc
´ 1 dH pvqdH puq `
´ 1 dH uc pvqdH uc puq
“ H psq ´ H psqH ptq `
H̄puq
H̄pvq
0 0
ż0 s u
´ pH uc psq ` H uc ptqqApsq ` 2 H uc pu´qdApuq ` 2H uc psqH uc ptq
ż s0
żs
” 1
ı
uc
uc
uc
uc
uc
uc
“ H psq ` H psqH ptq ` pH ptq ´ H puqq
´ 1 dH puq ` pApuq ´ H uc puqqdH uc puq
H̄puq
0
ż0 s
´ pH uc psq ` H uc ptqqApsq ` 2 H uc pu´qdApuq
0
żs
żs
uc
uc
uc
uc
uc
uc
uc
H uc puqdH uc puq
H puqdApuq `
“ H psq ` H psqH ptq ` H ptqApsq ´ H ptqH psq ´
0
0
żs
żs
uc
uc
uc
uc
` pApuq ´ H puqqdH puq ´ pH psq ` H ptqqApsq ` 2 H uc pu´qdApuq
0
0
żs
żs
żs
uc
uc
uc
uc
“ H psq ´ H psqApsq ´
H puqdApuq `
ApuqdH puq ` 2 H uc pu´qdApuq
0
0
0
żs
żs
“ H uc psq ´
H uc puqdApuq `
H uc pu´qdApuq
0
0
żs
żs
żs
uc
“
H̄puqdApuq ´
∆H puqdApuq “
H̄puqp1 ´ ∆ApuqqdApuq
uc

0

uc

uc

0

0

24

We conclude, as in van der Vaart and Wellner (1996), that
?
d
p ´ Aq ÝÑ
npA

ż¨
0

´ ż ¨ 1 ´ ∆Apuq
¯
1
uc
dM puq „ Gauss 0,
dApuq ,
H̄puq
H̄puq
0

where Gauss again indicates that the limit process is a Gaussian martingale.
The very same calculations hold true if each H uc is replaced with H1uc for the subdistribution
function of an uncensored type 1 event. Therefore, we have for Nelson-Aalen estimators for causespecific cumulative hazard functions that
?
d
p1 ´ A1 q ÝÑ
npA

ż¨
0

B

¯
´
1
2
p¨q
.
puq
„
Gauss
0,
σ
dMuc
1
1
H̄puq

Consistency of the Wild Bootstrap for the Univariate NelsonAalen Estimator

Proof. Without loss of generality, assume that 0, K P D1 for simplifying notation. Write D1 “
td0 , d1 , . . . , dJ u with the natural ordering dj ă dj`1 for all j “ 1, . . . , J. Then r0, KszD1 “
J
Ť
pdj´1 , dj q. As argued in Bluhmki et al. (2016), it is now straightforward to show that each
j“1

x1 ptq ´ W
x1 pdj´1 qqt on each interval rdj´1 , dj q, j “ 1, . . . , J, defines a square-integrable
process pW
martingale. Since such martingales can be extended to the right boundary of the time interval,
x1 pdj q ´ W
x1 pdj´1 q :“ limtÒd W
x1 ptq ´ W
x1 pdj´1 q, and this
we may define the boundary values W
j

procedure introduces square-integrable martingales on the whole rdj´1 , dj s.
x1 ptq ´
First, we notice that the conditional weak convergence in probability of the processes pW
x1 pdj´1 qqt on each interval rdj´1 , dj s, j “ 1, . . . , J, is already implied by exactly the same ReW
bolledo’s martingale central limit theorem arguments as in Bluhmki et al. (2016). Denote the limit
Gaussian martingale processes as pŨ1j ptqqtPrdj´1 ,dj s , j “ 1, . . . , J. Due to the martingale extension
above, Rebolledo’s limit theorem implies the almost sure continuity of Ũ1j on each time interval.
Furthermore, these are zero-mean processes with variance function t ÞÑ σ12 ptq ´ σ12 pdj´1 q.
Due to the continuity of the limit processes Ũ1j on the intervals rdj´1 , dj s, we are able to switch
from the Skorohod topology to the more convenient sup-norm metrization; see the discussion
in Section II.8 in Andersen et al. (1993). At each t “ dj , the weak conditional convergence
25

x1 pdj q holds in probability by the already argued convergence of all finitein distribution of ∆W
dimensional conditional distributions. Therefore, the independence of the (bootstrapped) NelsonAalen increments imply that, as n Ñ 8, the conditional distribution of
´
¯
x
x
x
x
x
x
x
∆W1 pd0 q, W1 pt0 q´W1 pd0 q, ∆W1 pd1 q, . . . , W1 ptJ´1 q´W1 pdJ´1 q, ∆W1 pdJ q

t0 Prd0 ,d1 s,...,tJ´1 PrdJ´1 ,dJ s

given F0 converges weakly in probability to the distribution of
pŨ1 pd0 q, Ũ11 pt1 q, Ũ1 pd1 q, . . . , Ũ1J ptJ q, Ũ1 pdJ qqt1 Prd0 ,d1 s,...,tJ PrdJ´1 ,dJ s
on the product Space R ˆ Drd0 , d1 s ˆ R ˆ ¨ ¨ ¨ ˆ DrdJ´1 , dJ s ˆ R equipped with the sup-max-norm.
Here, all components are independent, and the normally distributed random variables Ũ1 pdj q have
mean zero and variance ∆σ12 pdj q, j “ 1, . . . , J.
Applying the functional, which is continuous with respect to the max-sup-norm,
ψ : R ˆ Drd0 , d1 s ˆ R ˆ ¨ ¨ ¨ ˆ DrdJ´1 , dJ s ˆ R ÝÑ Drd0 , d1 s ˆ ¨ ¨ ¨ ˆ DrdJ´1 , dJ s ˆ R,

´
px0 , y1 pt1 q, x1 , . . . , yJ ptJ q, xJ qt1 ,...,tJ ÞÝÑ x0 ` y1 pt1 q, x0 ` y1 pd1 q ` x1 ` y2 pt2 q,
. . . , x0 `

J´1
ÿ

pyj pdj q ` xj q ` yJ ptJ q, x0 `

j“1

J
ÿ

¯
pyj pdj q ` xj q

j“1

t1 ,...,tJ

to the previous limit theorem, the continuous mapping theorem implies that: Given F0 and as
n Ñ 8, the conditional distribution of
´
¯
x
x
x
W1 pt1 q, . . . , W1 ptJ q, W1 pdJ q

t1 Prd0 ,d1 s,...,tJ PrdJ´1 ,dJ s

converges weakly (on the product-function space Drd0 , d1 s ˆ ¨ ¨ ¨ ˆ DrdJ´1 , dJ s ˆ R equipped with
the max-sup-norm) in probability to the distribution of
´

¯
U1 pt1 q, . . . , U1 ptJ q, U1 pdJ q

.

t1 Prd0 ,d1 s,...,tJ PrdJ´1 ,dJ s

Here the right boundary values are again considered as the left-hand limits given by the martingale
26

extension theorem. The process pU1 ptqqtPr0,Ks is a zero-mean Gaussian martingale with variance
function t ÞÑ σ12 ptq with, in general, discontinuous sample paths.
Finally, we apply the continuous functional
φ : Drd0 , d1 s ˆ ¨ ¨ ¨ ˆ DrdJ´1 , dJ s ˆ R ÝÑ Dr0, Ks,
J
´ÿ
¯
py1 pt1 q, . . . , yJ ptJ q, xJ qt1 ,...,tJ ÞÝÑ
yj ptq ¨ 1rdj´1 ,dj q ptq ` xJ ¨ 1tKu ptq

tPr0,Ks

j“1

x1 ptqqtPr0,Ks in probability.
in order to obtain the desired conditional weak convergence for pW

C

Derivation of the Asymptotic Covariance of Multiple NelsonAalen Estimators for Cumulative Cause-Specific Hazards

Consider the situation of two competing risks. In order to derive the asymptotic covariance of two
cause-specific Nelson-Aalen estimators, we first note that, as n Ñ 8,
2
ÿ
?
?
d
p
pj ´ Aj q ÝÑ
npA ´ Aq “
npA
U1 ` U2 “ U
j“1

by the continuous mapping theorem. The covariance function of U is given by
ż s^t
2

covpU psq, U ptqq “ σ ps ^ tq “
0

1 ´ ∆Apuq
dApuq
H̄puq

but on the other hand
covpU psq, U ptqq “

2
ÿ

covpUj psq, Uj ptqq ` covpU1 psq, U2 ptqq ` covpU1 ptq, U2 psqq.

j“1

Solving for the unknown covariances on the right-hand side of the previous display, we obtain
2 ż s^t
ÿ
1 ´ ∆Apuq
1 ´ ∆Aj puq
covpU1 psq, U2 ptqq ` covpU1 ptq, U2 psqq “
dApuq ´
dAj puq
H̄puq
H̄puq
0
j“1 0
ż s^t
ż s^t
∆A1 puq
∆A2 puq
“´
dA2 puq ´
dA1 puq
H̄puq
H̄puq
0
0

ż s^t

27

Due to symmetry and inductively, it follows that
ż s^t
σj` ps ^ tq “ covpUj psq, U` ptqq “ ´
0

∆Aj puq
dA` puq
H̄puq

for j ‰ `, j, ` “ 1, . . . , k, even in the situation of k P N competing risks.

D

Consistency of the Discretization-Adjusted Wild Bootstrap
for the Multivariate Nelson-Aalen Estimator

The proof of tightness follows along the same lines as that of Theorem 2.2. It only remains to
calculate the finite-dimensional marginal limit distributions. These are calculated with the help of
Theorem A.1 in Beyersmann et al. (2013): Therefore, we abbreviate
xj ptq “
W

n
ÿ

ξjji Zn,jji `

i“1

k
ÿ

signp` ´ jq

n
ÿ

ξj`i Z̃n,j`i `

i“1

`“1

k
ÿ

signp` ´ jq

n
ÿ

ξ`ji Z̃n,`ji .

i“1

`“1

Further, consider arbitrary points of time 0 ď t1 ď ¨ ¨ ¨ ď tm ď K and the vector
x1 pt1 q, . . . , W
x1 ptm q, W
x2 pt1 q, . . . , W
x2 ptm q . . . , W
xk pt1 q, . . . , W
xk ptm qq.
pW
x1 pt1 q, W
x1 pt2 qq and pW
x1 pt1 q, W
x2 pt2 qq
Due to analogy, it is enough to calculate the entries corresponding to pW
of the limit of the matrix
n ”
´
ı
ÿ
ÿ
ÿ
p :“ 1tj “ j̃u
Γ
Zn,jji pta qZn,jji ptb q `
Z̃n,j`i pta qZ̃n,j`i ptb q `
Z̃n,`ji pta qZ̃n,`ji ptb q
i“1
n ”
ÿ

´ 1tj ‰ j̃u

`‰j

`‰j

ı¯
Z̃n,j j̃i pta qZ̃n,j j̃i ptb q ` Z̃n,j̃ji pta qZ̃n,j̃ji ptb q

.

a,b“1,...,m; j,j̃“1,...,k

i“1

We start by calculating the entry for a “ 1, b “ 2, j “ j̃ “ 1, that is,
n ”
ÿ

Zn,11i pt1 qZn,11i pt2 q `

i“1

“n

k
ÿ

Z̃n,1`i pt1 qZ̃n,1`i pt2 q `

`“2
n ż t1
ÿ
i“1 0

k
ÿ

ı
Z̃n,`1i pt1 qZ̃n,`1i pt2 q

`“2
k
ÿ

n ż t1
n ż t1
” ÿ
ÿ
∆N1 puq dN`i puq
∆N` puq dN1i puq ı
Y puq ´ ∆N puq dN1i puq 1
`
n
`
n
Y puq
Y 2 puq
2 `“2 i“1 0 Y puq Y 2 puq
Y puq Y 2 puq
i“1 0

28

ż t1
“n
0

ż t1
“n
0

ż t1
ż t1
k
Y puq ´ ∆N puq dN1 puq 1 ÿ ”
∆N` puq dN1 puq ı
∆N1 puq dN` puq
`
`
n
n
Y 2 puq
Y puq
2 `“2
Y 2 puq Y puq
Y 2 puq Y puq
0
0
Y puq ´ ∆N1 puq dN1 puq
.
Y 2 puq
Y puq

Here, the last equality follows from ∆N1 dN` “ ∆N` dN1 for all `. By the Glivenko-Cantelli
theorem in combination with the continuous mapping theorem, it follows that the quantity in the
previous display converges to σ12 pt1 q in probability as n Ñ 8.
p for a “ 1, b “ 2, j “ j̃ “ 2:
Now, consider the entry of Γ
ż t1
ż t1
n ”
ı
ÿ
1”
∆N1 puq dN2 puq
∆N2 puq dN1 puq ı
Z̃n,12i pt1 qZ̃n,12i pt2 q ` Z̃n,21i pt1 qZ̃n,21i pt2 q “ ´ n
`
n
.
´
2
Y 2 puq Y puq
Y 2 puq Y puq
0
0
i“1
By the same arguments as before, this is a consistent estimator for σ12 pt1 q as n Ñ 8.

E

Derivation of the Asymptotic Covariance Function of an AalenJohansen Estimator in the Competing Risks Set-Up

.
In order to derive the the covariance function of
ż¨
ż¨
1 ´ F2 pu´q ´ F1 p¨q
F1 pu´q ´ F1 p¨q
dU1 puq `
dU2 puq,
1 ´ ∆Apuq
1 ´ ∆Apuq
0
0
at any ps, tq P r0, Ks2 , we exemplarily calculate the covariance function of the first integral and the
covariance function between both integrals. Hence, as the covariance function ps, tq ÞÑ σ12 ps ^ tq
of U1 only increases along the diagonal,
żt
´ ż s 1 ´ F pu´q ´ F psq
¯
1 ´ F2 pu´q ´ F1 ptq
2
1
dU1 puq,
dU1 puq
cov
1 ´ ∆Apuq
1 ´ ∆Apuq
0
0
ż s^t
p1 ´ F2 pu´q ´ F1 psqqp1 ´ F2 pu´q ´ F1 ptqq 2
“
dσ1 puq
p1 ´ ∆Apuqq2
0
ż s^t
p1 ´ F2 pu´q ´ F1 psqqp1 ´ F2 pu´q ´ F1 ptqq 1 ´ ∆A1 puq
“
dA1 puq.
p1 ´ ∆Apuqq2
H̄puq
0
29

Furthermore, we similarly have for the covariance between both integrals that
żt
´ ż s 1 ´ F pu´q ´ F psq
¯
F1 pu´q ´ F1 ptq
2
1
cov
dU1 puq,
dU2 puq
1 ´ ∆Apuq
1 ´ ∆Apuq
0
0
ż s^t
p1 ´ F2 pu´q ´ F1 psqqpF1 pu´q ´ F1 ptqq
“
dσ12 puq
p1 ´ ∆Apuqq2
0
ż s^t
p1 ´ F2 pu´q ´ F1 psqqpF1 pu´q ´ F1 ptqq ∆A1 puq
dA2 puq.
“´
p1 ´ ∆Apuqq2
H̄puq
0
Finally, including also the remaining two analogous terms, we obtain the following asymptotic
covariance function of the Aalen-Johansen estimator for the first cumulative incidence function as
the sum of all four covariance functions:
ż s^t
ps, tq ÞÑ

p1 ´ F2 pu´q ´ F1 psqqp1 ´ F2 pu´q ´ F1 ptqq 1 ´ ∆A1 puq
dA1 puq
p1 ´ ∆Apuqq2
H̄puq
0
ż s^t
pF1 pu´q ´ F1 psqqpF1 pu´q ´ F1 ptqq 1 ´ ∆A2 puq
`
dA2 puq
p1 ´ ∆Apuqq2
H̄puq
0
ż s^t
p1 ´ F2 pu´q ´ F1 psqqpF1 pu´q ´ F1 ptqq ∆A1 puq
´
dA2 puq
p1 ´ ∆Apuqq2
H̄puq
0
ż s^t
p1 ´ F2 pu´q ´ F1 ptqqpF1 pu´q ´ F1 psqq ∆A2 puq
dA1 puq.
´
p1 ´ ∆Apuqq2
H̄puq
0

References
A. Allignol, M. Schumacher, and J. Beyersmann. A note on variance estimation of the AalenJohansen estimator of the cumulative incidence function in competing risks, with a view towards
left-truncated data. Biom. J., 52(1):126–137, 2010.
P. K. Andersen, Ø. Borgan, R. D. Gill, and N. Keiding. Statistical Models Based on Counting
Processes. Springer, New York, 1993.
J. Beyersmann, M. Pauly, and S. Di Termini. Weak Convergence of the Wild Bootstrap for the
Aalen-Johansen Estimator of the Cumulative Incidence Function of a Competing Risk. Scandinavian Journal of Statistics, 40(3), 2013.
T. Bluhmki, D. Dobler, J. Beyersmann, and M. Pauly. The Wild Bootstrap for Multivariate NelsonAalen Estimators. Submitted article. arXiv preprint:1602.02071v2, 2016.
D. R. Cox. Regression Models and Life-Tables. Journal of the Royal Statistical Society, Series B
(Methodological), 34(2):187–220, 1972.
30

D. Dobler and M. Pauly. Bootstrapping Aalen-Johansen processes for competing risks: Handicaps,
solutions, and limitations. Electronic Journal of Statistics, 8(2):2779–2803, 2014.
D. Dobler, J. Beyersmann, and M Pauly. Non-strange weird resampling for complex survival data.
Preprint arXiv:1507.02838, 2015.
D. Y. Lin. Non-parametric inference for cumulative incidence functions in competing risks studies.
Statistics and Medicine, 16:901–910, 1997.
D. Y. Lin, L. J. Wei, and Z. Ying. Checking the Cox model with cumulative sums of martingalebased residuals. Biometrika, 80(3):557–572, 1993.
T. Martinussen and T. H. Scheike. Dynamic regression models for survival data. Statistics for
Biology and Health. Springer, New York, 2006.
R Development Core Team.
R: A Language and Environment for Statistical Computing.
R Foundation for Statistical Computing, Vienna, Austria, 2016.
URL
http://www.R-project.org.
A. W. van der Vaart and J. Wellner. Weak Convergence and Empirical Processes. Springer, New
York, 1996.

31

