A Discontinuity Adjustment for Subdistribution Function Confidence Bands Applied to Right-Censored Competing Risks Data

,


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 ∆t P pT P rt, t `∆ts | T ě tq 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, N ji ptq " 1t" individual i is observed to fail in r0, ts due to risk j" u and Y i 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 H uc 1 ptq " EpN 11 ptqq and Hptq " EpY 1 ptqq.Now, the sum N 1 ptq " ř n i"1 N 1i ptq has the compensator where A 1 ptq " ş t 0 Hpuq ´1dH uc 1 puq is the cumulative hazard function for a type 1 event.Therefore, pM 1 ptq " N 1 ptq ´Λ1 ptqq t is a square-integrable martingale.The Nelson-Aalen estimator for A 1 ptq is defined as the counting process integral p A 1 ptq " ş t 0 Y ´1dN 1 , where we let 0{0 " 0. The limit covariance function of the normalized process W 1 " ?np p A 1 ´A1 q is determined by the predictable and optional variation processes of the square-integrable martingales M 1i : t Þ Ñ N 1i ptq ´Λ1i ptq.In the possibly discontinuous case, these are given by xM 1i yptq " Λ 1i ptq ´ż t 0 ∆Λ 1i puqdΛ 1i puq and rM 1i sptq " N 1i ptq ´2 ż t 0 ∆Λ 1i puqdN 1i puq `ż t 0 ∆Λ 1i puqdΛ 1i puq, 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 possible.This is as well seen with the help of covariation processes: Let A 1 : t Þ Ñ ş t 0 JpuqA 1 pduq{Y puq denote the compensator of the Nelson-Aalen estimator p A 1 , where Jpuq " 1tY puq ą 0u.By the usual martingale theory, we have the optional covariation process where the last equality follows from the orthogonality of all involved single martingales.Inserting the above representation for the optional variation process yields In a last step, we replace the unknown increments of the cumulative hazard function ∆A 1 puq with Nelson-Aalen increments ∆ p A 1 puq " Y puq ´1∆N 1 puq and obtain the Greenwood-type variance estimator p σ 2 1 ptq " n 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 A 1 pKq ă 8.A simple application of Glivenko-Cantelli theorems in combination with the continuous mapping theorem show that this Greenwood-type estimator is uniformly consistent for σ 2 1 ptq " ş t 0 p1 ∆A 1 puqq Hpuq ´1dA 1 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 ∆A 1 in the limit variance is crucially different from the continuous case in which the limit variance function reduces to σ2 1 : t Þ Ñ ş t 0 H´1 dA 1 .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 Aalen-Johansen 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 Nelson-

Aalen Estimator
The wild bootstrap is typically applied in the following way (see e.g.Lin 1997, Martinussen and Scheike 2006, Beyersmann et al. 2013or Bluhmki et al. 2016 for resampling more general counting process-based Nelson-Aalen estimators): First, the (independent) martingale increments dM 1i puq in the martingale representation of ?np p A 1 ´A1 q are replaced by independently weighted counting process increments, i.e., by ξ i dN 1i 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 constitutes a square-integrable martingale with predictable and optional variation processes 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 dM 1i puq should not be replaced by ξ i dN 1i puq but rather by something which -considered again as a martingale -reproduces the correct covariance structure.Therefore, we replace which results in the following wild bootstrap resampling version of the normalized Nelson-Aalen estimator: x W 1 ptq " ?n In a way similar to Bluhmki et al. (2016), one can show that p x W 1 ptqq tPr0,Ks is a martingale with respect to the filtration pF t q tPr0,Ks given by F t " σpξ i N 1i puq, N 1i pvq, Y i pvq : u P r0, ts, v P r0, Ks, i " 1, . . ., nq.
Analogously, it is easy to see that its predictable and optional variation processes are given by It is well-known that the Greenwood-type variance estimator p σ 2 1 is uniformly consistent for σ 2 1 .Assuming 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 Pó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 D 1 " tt P r0, Ks : ∆A 1 ptq ą 0u the subset of time points for which ties among the type 1 events are possible.Similarly, we define D 2 , . . ., D k 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.
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 F 0 and as n Ñ 8, we have the following conditional weak convergence 1 q in outer probability on the càdlàg function space Dr0, Ks equipped with the supremum distance topology, where U 1 is a Gaussian zero-mean martingale with variance function t Þ Ñ σ 2 1 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 Nelson-Aalen 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 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.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 M 1i , . . ., M ki , 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 different: Here we have for the martingales M 1i and M 2i of Section 1 that cf. the derivations in Section II.4 in Andersen et al. (1993).Define W j " ?np p A j ´Aj q, j " 1, . . ., k.The above variation processes are strong evidence that the normalized Nelson-Aalen estimators W 1 , W 2 , . . ., W k 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 càdlàg function space D k r0, Ks, equipped with the max-sup-norm, that where U 1 , U 2 , . . ., U k are zero-mean Gaussian martingales with variance functions t Þ Ñ σ 2 j ptq " ż t 0 1 ´∆A j puq Hpuq dA j puq, j " 1, 2, . . ., k and covariance functions (for j ‰ ℓ) ps, tq Þ Ñ covpU j psq, U ℓ ptqq " ´ż s^t 0 ∆A ℓ puq Hpuq dA j puq ": σ jℓ ps ^tq.
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.random variables with Epξ 111 q " 0 and Epξ 2 111 q " 1, which are also independent of the data.Denote by N " ř k j"1 N j the number of all kinds of events.Then, the single components of the wild bootstrap version of the multivariate Nelson-Aalen estimator pW 1 , . . ., W k q are given by 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 F 0 and as n Ñ 8, we have the following conditional weak convergence on the product càdlàg function space D k r0, Ks, equipped with the max-supnorm: where pU 1 , U 2 , . . ., U k 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: are the usual Greenwood-type (co)variance estimators and 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 estimators p σ j and p σ jℓ are shown to be uniformly consistent for σ 2 j and σ jℓ , respectively.For the wild bootstrap-type (co)variance estimators, we can parallel the arguments in the proof of Theorem 2.2, after first assuming the existence of fourth moments Epξ 4 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.In points of discontinuity, which are finitely many by assumption, we approximate p p σ 2 j by p σ 2 j and apply the conditional Chebyshev inequality (given F 0 ) in order to show the negligibility of the differences p p σ 2 j ´p σ 2 j in probability.The last argument can be repeated for the covariance estimators.A final application of the continuous mapping theorem yields p σ 2 j , p p σ 2 j p ÝÑ σ j and p σ jℓ , p p σ jℓ p ÝÑ σ jℓ uniformly on r0, Ks in (conditional outer) probability, for all j ‰ ℓ as n Ñ 8.

Extension to Aalen-Johansen Estimators for Cumulative Incidence Functions in Competing Risks
Denote the cumulative incidence functions by F j 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 U 2 , A 2 , and F 2 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 theorem for the Aalen-Johansen estimator p F 1 ptq " ş t 0 p Spu´qd p A 1 puq for the cumulative incidence function F 1 ptq " ş t 0 Spu´qdA 1 puq: where U F 1 is a zero-mean Gaussian process with covariance function 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 ?np p F 1 ´F1 q is given by where x W 1 and x W 2 are again the wild bootstrap versions of the Nelson-Aalen estimator as presented 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 F 0 and as n Ñ 8, we have the following weak convergence on the càdlàg function space Dr0, Ks, equipped with the sup-norm: where U F 1 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 pX i q, maxp1, Y pX i qq ´1q, where X i 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.
Estimators for σ 2 F 1 and its wild bootstrap variant are obtained similarly as such estimators for the Nelson-Aalen (co)variances, i.e. via plug-in: Similarly, p p σ 2 F 1 is obtained by replacing all estimators p σ j and p σ jℓ , j ‰ ℓ, by their wild bootstrap counterparts p p σ j and p p σ jℓ , respectively.Their uniform (conditional) consistencies for σ 2 F 1 on r0, Ks 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 F 1 can be deduced.In particular, let ϕpsq " logp´logp1 ´sqq be a transformation applied to F 1 in order to ensure band boundaries between 0 and 1 and let g 1 psq " logp1 ´p F 1 psqq{p ρpsq and g 2 psq " logp1 ´p F 1 psqq{p1 `p ρ 2 psqq be weight functions leading to the usual equal precision and Hall-Wellner bands, respectively; see Andersen et al. (1993).

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 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 F 1 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, 1 k , 2 k , . . .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 are presented in Figure 1, with the exception of the continuous function F 1 .Here rss denotes the integer closest to s P R. For simulating data, which have the desired cumulative incidence function F p,k 1 , it is mandatory to first round the event times T i , and then generate the event types ε i in a second step, according to the formula where S : t Þ Ñ expp´tq denotes the survival function of the continuous random variables T i .
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 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 F 1 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 F p,k 1 .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.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 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.

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)  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 H uc 1 for the subdistribution function of an uncensored type 1 event.Therefore, we have for Nelson-Aalen estimators for causespecific cumulative hazard functions that B Consistency of the Wild Bootstrap for the Univariate Nelson-

Aalen Estimator
Proof.Without loss of generality, assume that 0, K P D 1 for simplifying notation.Write D 1 " td 0 , d 1 , . . ., d J u with the natural ordering d j ă d j`1 for all j " 1, . . ., J. Then r0, KszD 1 " J Ť j"1 pd j´1 , d j q.As argued in Bluhmki et al. (2016), it is now straightforward to show that each process p x W 1 ptq ´x W 1 pd j´1 qq t on each interval rd j´1 , d j q, j " 1, . . ., J, defines a square-integrable martingale.Since such martingales can be extended to the right boundary of the time interval, we may define the boundary values x W 1 pd j q ´x W 1 pd j´1 q :" lim tÒd j x W 1 ptq ´x W 1 pd j´1 q, and this procedure introduces square-integrable martingales on the whole rd j´1 , d j s.
First, we notice that the conditional weak convergence in probability of the processes p x W 1 ptq x W 1 pd j´1 qq t on each interval rd j´1 , d j s, j " 1, . . ., J, is already implied by exactly the same Rebolledo's martingale central limit theorem arguments as in Bluhmki et al. (2016).Denote the limit Gaussian martingale processes as p Ũ1j ptqq tPrd j´1 ,d j s , j " 1, . . ., J. Due to the martingale extension above, Rebolledo's limit theorem implies the almost sure continuity of Ũ1j on each time interval.
Due to the continuity of the limit processes Ũ1j on the intervals rd j´1 , d j 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 " d j , the weak conditional convergence in distribution of ∆ x W 1 pd j q holds in probability by the already argued convergence of all finitedimensional conditional distributions.Therefore, the independence of the (bootstrapped) Nelson-Aalen increments imply that, as n Ñ 8, the conditional distribution of given F 0 converges weakly in probability to the distribution of p Ũ1 pd 0 q, Ũ11 pt 1 q, Ũ1 pd 1 q, . . ., Ũ1J pt J q, Ũ1 pd J qq t 1 Prd 0 ,d 1 s,...,t J Prd J´1 ,d J s on the product Space R ˆDrd 0 , d 1 s ˆR ˆ¨¨¨ˆDrd J´1 , d J s ˆR equipped with the sup-max-norm.
Here, all components are independent, and the normally distributed random variables Ũ1 pd j q have mean zero and variance ∆σ 2 1 pd j q, j " 1, . . ., J.
by the continuous mapping theorem.The covariance function of U is given by covpU psq, U ptqq " σ 2 ps ^tq " 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 Further, consider arbitrary points of time 0 ď t 1 ď ¨¨¨ď t m ď K and the vector Due to analogy, it is enough to calculate the entries corresponding to p x W 1 pt 1 q, x W 1 pt 2 qq and p x W 1 pt 1 q, x W 2 pt 2 qq of the limit of the matrix p Γ :" ´1tj " ju n ÿ i"1 " Z n,jji pt a qZ n,jji pt b q `ÿ ℓ‰j Zn,jℓi pt a q Zn,jℓi pt b q `ÿ ℓ‰j Zn,ℓji pt a q Zn,ℓji pt b q ı ´1tj ‰ ju n ÿ i"1 " Zn,j ji pt a q Zn,j ji pt b q `Z n, jji pt a q Zn, jji pt b q ı¯a ,b"1,...,m; j, j"1,...,k .
We start by calculating the entry for a " 1, b " 2, j " j " 1, that is, Here, the last equality follows from ∆N 1 dN ℓ " ∆N ℓ dN 1 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 σ 2 1 pt 1 q in probability as n Ñ 8. Now, consider the entry of p Γ for a " 1, b " 2, j " j " 2: By the same arguments as before, this is a consistent estimator for σ 12 pt 1 q as n Ñ 8.  We consider the same competing risks setup as above, i.e., we assume that there are k P N competing risks and n P N i.i.d.random event times T 1 , ..., T n , which are independently rightcensored and distributed as a random variable T " S. Here, S denotes the survival function, i.e., Sptq " P pT ą tq for all t ě 0; S need not be continuous.Then, we denote the probability that an individual is under observation at time t´, that is, just before time t, by Hptq " P pminpT, Cq ě tq " Spt´qGpt´q for all t ě 0. Here, C " G with survival function Gptq " P pC ą tq denotes a generic censoring time which is assumed to be independent of T .Also, t Þ Ñ f pt´q denotes the left-continuous version of a right-continuous function f : r0, 8q Ñ R. Furthermore, let p A j denote the cause-specific Nelson-Aalen estimator for the cumulative hazard function A j of type j events, p S the Kaplan-Meier estimator for the Survival function S, and p F j " ş .0 p Spu´qd p A j puq the Aalen-Johansen estimator for the cumulative incidence function F j " ş .0 Spu´qdA j puq for all j P t1, ..., ku.In addition to the assumptions in the main manuscript, it is actually required that HpKq ą 0 for K ě 0 to ensure finite variances σ 2 j pKq, j P t1, . . ., ku, in Theorem 3.1 above.Theorem 4.1 above states for k " 2 competing risks that as n Ñ 8 on the càdlàg space Dr0, Ks equipped with the sup-norm, where U F 1 is a zero-mean Gaussian process with covariance function where A " ř k j"1 A j and D " tt P r0, Ks : ∆Aptq ą 0u is the set of discontinuity time points.However, we found that the right-continuous versions F 1 , F 2 , S must appear in the covariance function above in all occurrences of F 1 pu´q, F 2 pu´q, Spu´q.
In order to prove this, we go one step back: By Theorem 3.1 above,

?
n ´p A 1 ´A1 , ..., p A k ´Ak ¯d Ý Ñ pU 1 , ..., U k q holds as n Ñ 8 on the product space D k r0, Ks equipped with the max-sup norm, where U 1 , ..., U k are zero-mean Gaussian-martingales with covpU j ptq, U j psqq " ż t^s 0 1 ´∆A j puq Hpuq dA j puq ": σ 2 j pt ^sq, covpU j ptq, U ℓ psqq " ´ż t^s 0 ∆A ℓ puq Hpuq dA j puq ": σ jℓ pt ^sq for all t, s P r0, Ks, j, ℓ P t1, ..., ku, j ‰ ℓ.We further note that the limit pU 1 , ..., U k q is separable since G uc 1 , ..., G uc k and G in Appendix A are tight, which follows by the main empirical central limit theorems in Vaart and Wellner (2023) The covariance function can be calculated analogously to Appendix E. Here, the last sum may be simplified to ř uPD,uďs,t

Figure 1 :
Figure 1: Cumulative incidents functions F p,k 1 underlying the present simulations.

Figure 2 :
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.

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: ´F1 ptqqpF 1 pu´q ´F1 psqq Hpuq ∆A 2 puq p1 ´∆Apuqq 2 dA 1 puq.
Remark 2.3.Consider the case of only one risk, i.e. k " 1 and W " W 1 .The Kaplan-Meier estimator for the survival function Sptq " P pT ą tq is

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

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

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

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

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

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

Table 7 :
Simulated coverage probabilities of confidence bands for F 1 in per cent where p " 0 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?Note that Aptq " ş t 0 H uc pduq{ Hpuq.Thus, for s ď t, the covariance function of M uc at ps, tq is covpU j psq, U j ptqq `covpU 1 psq, U 2 ptqq `covpU 1 ptq, U 2 psqq.Solving for the unknown covariances on the right-hand side of the previous display, we obtain covpU 1 psq, U 2 ptqq `covpU 1 ptq, U 2 psqq " Due to symmetry and inductively, it follows that σ jℓ ps ^tq " covpU j psq, U ℓ ptqq " ´ż s^t 0 ∆A j puq Hpuq dA ℓ puq for j ‰ ℓ, j, ℓ " 1, ..., k, even in the situation of k P N competing risks.D Consistency of the Discretization-Adjusted Wild Bootstrapfor the Multivariate Nelson-Aalen Estimator Ks 2 , we exemplarily calculate the covariance function of the first integral and the covariance function between both integrals.Hence, as the covariance function ps, tq Þ Ñ σ 2 1 ps ^tq of U 1 only increases along the diagonal, , as in Example 3.10.20.BV 1 r0, Ks, g P Dr0, Ks, where BV 1 r0, Ks denote the set of all càdlàg functions Dr0, Ks of total variation bounded by 1.As in Example 3.10.33 in Vaart and Wellner (2023), the functional delta method yields on D 2 r0, Ks, where the integral is defined by integration by parts since U 1 `U2 is not of bounded variation.Hence, we get as n Ñ 8 on D 2 r0, Ks ˆBV 1 r0, Ks by Slutsky's lemma.Note that the mapψ : D 2 r0, Ks ˆBV 1 r0, Ks Ñ Dr0, Ks, p Ã, B, Cq Þ Ñ Ãp.q Cp.q ´ż .Theorem E.1 (Corrected Theorem 4.1).As n Ñ 8, we have on the càdlàg space Dr0, Ks where U F 1 is a zero-mean Gaussian process with covariance functionσ 2 F 1 : ps, tq Þ Ñ