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

TL;DR

This paper proposes a modified wild bootstrap method and variance estimator for more accurate confidence bands in survival analysis with tied event times, especially in competing risks data.

Contribution

It introduces a discontinuity adjustment technique for bootstrap and variance estimation, improving inference accuracy in tied survival data.

Findings

Enhanced bootstrap method with discontinuity adjustment improves coverage accuracy.

Simulation studies demonstrate better finite-sample performance.

Application to ICU pneumonia data confirms practical utility.

Abstract

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.

Tables

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

set-ups		$N(0,1)$		$Poi\left(1\right)-1$		weird
$p$	$n$	old	new	old	new	old	new
	50	87.31	88.47	89.44	90.61	88.93	90.36
	75	89.48	90.39	91.55	92.24	91.13	92.01
	100	90.64	91.26	92.32	92.96	92.21	92.64
	125	91.21	91.96	92.36	93.07	93.43	93.76
25%	150	91.12	91.66	92.48	92.98	93.37	93.75
	175	91.69	92.14	93.10	93.21	93.20	93.63
	200	91.99	92.46	93.31	93.80	93.63	94.05
	225	92.51	92.98	93.29	93.67	93.58	94.03
	250	92.48	92.88	93.65	93.91	93.41	93.84
	50	87.41	89.31	89.09	91.28	89.76	91.29
	75	89.35	90.99	91.28	92.71	91.03	92.59
	100	90.01	91.53	92.39	93.64	91.88	93.29
	125	90.56	91.92	92.95	93.90	92.29	93.43
50%	150	90.59	91.93	92.44	93.38	93.04	94.05
	175	91.21	92.22	92.58	93.49	92.66	93.73
	200	91.22	92.28	92.30	93.31	93.16	94.19
	225	91.44	92.61	92.70	93.73	93.25	94.16
	250	91.72	92.82	93.13	94.27	93.07	94.05
	50	87.23	90.60	89.12	91.56	89.47	92.11
	75	90.09	92.56	90.96	93.22	91.29	93.77
	100	90.58	92.86	92.13	94.23	91.59	93.83
	125	90.83	93.03	92.13	94.26	91.75	94.08
75%	150	90.84	93.11	91.96	93.94	92.29	94.35
	175	90.83	93.13	91.66	93.87	92.25	94.03
	200	90.91	93.11	92.45	94.46	92.50	94.49
	225	91.15	93.43	92.52	94.49	92.22	94.12
	250	91.76	93.76	92.38	94.36	92.73	94.53
	50	87.71	91.98	89.08	92.91	89.56	93.19
	75	88.79	92.70	90.26	93.74	90.21	93.64
	100	89.89	93.25	90.90	94.06	90.47	93.84
	125	89.76	93.14	90.92	93.88	90.44	93.64
100%	150	89.64	93.21	90.84	93.97	91.00	93.87
	175	90.10	93.24	90.52	93.52	91.25	93.99
	200	90.27	93.55	91.06	94.00	91.09	94.12
	225	90.01	93.07	91.30	94.19	91.08	94.26
	250	89.71	93.18	90.92	93.89	90.76	93.86

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

set-ups		$N(0,1)$		$Poi\left(1\right)-1$		weird
$p$	$n$	old	new	old	new	old	new
	50	87.72	88.64	89.37	90.69	89.87	90.84
	75	89.22	90.16	91.86	92.55	91.75	92.35
	100	90.68	91.22	92.43	93.07	93.06	93.31
	125	91.13	91.61	93.05	93.50	93.42	93.79
25%	150	91.63	92.16	93.5	93.77	93.48	93.79
	175	91.79	92.06	93.87	94.19	93.75	93.93
	200	91.93	92.37	93.84	93.94	93.70	93.88
	225	92.40	92.77	93.54	93.89	93.73	93.97
	250	92.17	92.61	94.03	94.26	94.24	94.56
	50	87.15	88.59	89.20	90.73	89.18	90.50
	75	89.78	90.80	91.46	92.47	91.18	92.22
	100	90.24	91.32	92.58	93.25	92.53	93.47
	125	91.17	92.11	92.99	93.70	93.06	93.76
50%	150	91.63	92.52	93.25	93.98	93.17	93.95
	175	92.65	93.29	92.67	93.34	93.06	93.71
	200	92.21	92.87	93.36	94.03	93.51	94.23
	225	92.20	92.98	94.07	94.80	93.09	93.96
	250	92.96	93.52	93.71	94.23	93.37	93.93
	50	87.63	89.82	89.19	91.41	88.95	90.91
	75	89.52	91.46	91.97	93.42	91.36	92.96
	100	90.30	91.90	92.45	93.76	92.52	93.95
	125	91.44	92.88	92.93	94.03	92.83	94.14
75%	150	91.49	92.76	93.20	94.50	93.22	94.23
	175	91.50	92.76	92.98	94.18	93.66	94.63
	200	92.00	93.08	93.92	94.94	93.75	94.85
	225	92.20	93.12	93.57	94.56	93.17	94.13
	250	92.49	93.78	93.51	94.55	93.69	94.81
	50	88.28	91.10	90.21	92.75	90.31	92.68
	75	89.81	92.2	91.70	93.90	91.77	93.57
	100	90.54	92.66	92.10	93.98	92.09	93.96
	125	91.12	93.19	92.58	94.36	92.74	94.47
100%	150	91.41	93.44	92.83	94.61	92.77	94.77
	175	91.60	93.59	92.84	94.57	92.84	94.79
	200	91.49	93.46	93.03	94.76	93.25	94.97
	225	92.32	93.98	92.92	94.70	93.38	94.95
	250	91.81	93.81	92.79	94.58	92.98	94.73

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

set-ups		$N(0,1)$		$Poi\left(1\right)-1$		weird
$p$	$n$	old	new	old	new	old	new
	50	87.21	88.03	89.93	90.94	90.15	90.69
	75	89.75	90.33	92.18	92.50	91.78	92.38
	100	91.01	91.61	93.12	93.52	92.71	93.08
	125	92.09	92.29	93.99	94.18	93.34	93.67
25%	150	91.87	92.39	93.63	93.78	93.69	93.85
	175	92.44	92.71	94.03	94.40	94.19	94.35
	200	92.76	93.03	94.18	94.35	93.76	94.07
	225	93.23	93.44	94.46	94.77	94.02	94.23
	250	92.84	93.11	93.88	94.10	94.16	94.41
	50	88.02	89.22	89.57	90.86	89.59	90.81
	75	89.88	90.74	92.37	93.28	92.34	93.00
	100	90.95	91.69	93.19	93.74	93.02	93.71
	125	91.98	92.46	93.44	93.83	93.91	94.2
50%	150	91.90	92.54	93.60	94.01	93.97	94.27
	175	92.95	93.40	94.00	94.42	94.11	94.45
	200	92.79	93.11	94.05	94.39	93.68	94.00
	225	92.78	93.25	94.16	94.53	94.28	94.67
	250	93.01	93.28	93.92	94.30	93.77	94.03
	50	88.42	89.91	90.81	92.12	90.51	91.96
	75	90.28	91.50	92.73	93.94	92.85	93.67
	100	91.43	92.45	93.30	94.21	93.02	93.85
	125	91.41	92.30	93.84	94.58	93.57	94.08
75%	150	92.73	93.31	94.01	94.77	93.54	94.25
	175	92.44	93.14	94.44	95.04	94.09	94.82
	200	93.04	93.84	94.13	94.65	94.02	94.58
	225	93.32	93.94	94.41	94.95	94.08	94.69
	250	93.03	93.59	94.27	94.85	94.29	94.77
	50	88.87	90.73	91.04	92.79	91.16	92.69
	75	90.95	92.58	93.21	94.28	93.02	94.24
	100	91.98	93.28	93.74	95.01	93.74	94.67
	125	92.91	93.94	93.69	94.88	94.13	95.13
100%	150	92.82	93.88	94.18	95.10	94.58	95.49
	175	93.16	94.25	94.88	95.59	94.85	95.74
	200	93.53	94.51	94.38	95.30	94.59	95.29
	225	94.06	94.76	94.63	95.59	94.69	95.39
	250	93.12	94.31	94.95	95.79	94.79	95.60

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

set-ups		$N(0,1)$		$Poi\left(1\right)-1$		weird
$p$	$n$	old	new	old	new	old	new
	50	88.65	89.93	89.33	90.79	88.85	90.18
	75	89.82	90.60	90.59	91.54	90.39	91.32
	100	90.51	91.11	90.93	91.82	91.13	91.88
	125	90.92	91.69	91.11	91.82	91.07	91.92
25%	150	91.41	91.70	91.53	92.20	92.21	92.74
	175	91.79	92.17	92.75	93.17	92.65	93.16
	200	92.24	92.67	92.18	92.59	92.26	92.73
	225	92.22	92.70	92.75	93.20	92.64	93.14
	250	92.38	93.03	92.53	93.05	92.44	92.68
	50	87.72	90.36	88.52	90.87	88.77	91.14
	75	89.17	91.06	89.54	91.43	89.92	91.87
	100	89.86	91.62	90.26	91.87	90.60	92.18
	125	90.86	92.28	91.27	92.69	91.00	92.39
50%	150	90.99	92.28	91.34	92.77	91.15	92.57
	175	91.41	92.67	91.76	92.84	91.12	92.42
	200	91.77	92.99	92.20	93.57	91.78	92.86
	225	91.70	92.87	91.72	92.95	91.94	93.14
	250	91.44	92.84	92.00	93.20	91.65	93.04
	50	87.35	90.83	88.74	92.14	88.32	91.82
	75	89.97	92.41	90.57	93.39	90.09	93.16
	100	89.89	92.60	90.32	93.29	90.19	93.01
	125	90.52	92.98	90.95	93.38	91.03	93.55
75%	150	90.61	93.03	91.15	93.80	91.16	93.67
	175	91.03	93.41	91.54	93.95	90.88	93.37
	200	91.29	93.53	90.98	93.49	91.30	93.41
	225	91.63	93.76	91.70	93.52	91.11	93.61
	250	91.44	93.55	91.55	93.78	91.95	94.24
	50	88.99	92.88	89.60	93.67	89.68	93.41
	75	90.32	93.74	90.71	94.31	90.23	93.99
	100	90.59	94.04	91.17	94.37	90.60	94.14
	125	91.27	94.59	91.79	94.81	91.56	94.55
100%	150	91.63	94.59	91.42	94.50	91.20	94.60
	175	91.53	94.56	91.95	94.80	90.97	94.22
	200	91.78	94.65	91.82	94.68	91.76	94.68
	225	91.94	94.92	91.87	94.81	91.59	94.56
	250	92.40	94.92	91.86	94.76	91.35	94.46

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

set-ups		$N(0,1)$		$Poi\left(1\right)-1$		weird
$p$	$n$	old	new	old	new	old	new
	50	87.21	88.70	89.27	90.49	89.00	90.36
	75	89.23	90.08	90.91	91.57	90.65	91.51
	100	90.54	90.86	91.86	92.55	91.40	91.89
	125	91.49	91.96	92.45	92.86	91.59	92.25
25%	150	92.52	92.80	92.31	92.94	92.18	92.65
	175	92.61	92.88	92.86	93.06	92.71	93.23
	200	92.25	92.76	92.66	92.88	92.80	93.20
	225	92.92	93.41	93.02	93.32	93.00	93.38
	250	92.37	92.78	92.90	93.29	93.19	93.46
	50	87.32	89.05	88.85	90.49	89.17	90.77
	75	89.94	91.20	91.08	92.27	90.65	91.81
	100	91.24	92.37	91.08	92.32	91.51	92.59
	125	91.05	92.16	91.61	92.49	91.88	92.77
50%	150	91.99	92.82	91.92	92.79	92.00	93.01
	175	92.37	93.10	92.53	93.21	92.66	93.38
	200	92.24	92.79	92.42	93.37	92.59	93.25
	225	92.29	93.22	92.75	93.44	92.87	93.64
	250	92.05	92.73	92.70	93.38	92.05	92.74
	50	88.15	90.60	88.32	91.01	89.49	91.71
	75	90.41	92.20	91.34	92.98	90.39	92.28
	100	91.25	92.74	90.97	92.66	91.19	93.00
	125	91.76	93.21	91.99	93.34	91.57	93.09
75%	150	91.67	92.85	91.73	93.04	91.90	93.31
	175	92.30	93.57	92.50	93.75	92.08	93.34
	200	92.25	93.46	92.11	93.50	92.05	93.46
	225	92.40	93.63	92.43	93.88	92.89	94.09
	250	92.70	93.79	92.38	93.66	92.56	93.67
	50	88.99	92.00	89.50	92.73	89.59	92.83
	75	91.10	93.57	91.14	93.59	90.91	93.76
	100	90.84	93.12	92.10	94.27	91.66	93.77
	125	92.06	93.80	92.06	94.18	92.25	94.27
100%	150	91.95	93.89	92.56	94.65	92.60	94.59
	175	92.41	94.23	92.87	94.70	92.27	94.35
	200	92.34	94.21	93.11	94.88	92.99	94.88
	225	92.73	94.31	93.25	94.85	93.06	94.89
	250	93.32	94.90	92.74	94.41	92.90	94.61

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

set-ups		$N(0,1)$		$Poi\left(1\right)-1$		weird
$p$	$n$	old	new	old	new	old	new
	50	88.39	89.53	89.93	91.16	89.14	90.05
	75	90.39	90.99	90.67	91.39	90.74	91.52
	100	91.39	91.89	91.75	92.26	91.92	92.30
	125	91.69	92.17	92.54	93.03	92.52	92.90
25%	150	92.07	92.51	92.57	92.83	92.92	93.28
	175	92.24	92.55	92.76	92.97	92.8	93.07
	200	93.01	93.27	92.58	93.04	92.96	93.21
	225	93.11	93.32	92.99	93.25	93.16	93.42
	250	92.88	93.03	93.60	93.69	93.05	93.31
	50	88.77	90.02	89.28	90.66	89.53	90.89
	75	90.43	91.37	91.27	92.29	90.87	91.92
	100	91.74	92.45	91.99	92.74	91.94	92.62
	125	92.12	92.82	92.49	93.11	92.08	92.61
50%	150	92.12	92.59	93.05	93.60	92.56	93.18
	175	92.43	93.00	92.96	93.25	93.00	93.45
	200	92.94	93.22	93.08	93.50	92.69	93.28
	225	92.67	93.01	93.20	93.92	93.51	93.90
	250	93.41	93.58	93.03	93.47	92.74	93.25
	50	89.17	90.83	89.70	91.53	89.65	91.36
	75	90.99	92.23	91.29	92.41	91.21	92.32
	100	91.15	92.12	92.36	93.24	92.26	93.37
	125	92.21	93.15	92.57	93.44	92.49	93.50
75%	150	92.64	93.39	93.08	93.98	92.68	93.46
	175	92.50	93.21	93.25	93.91	92.67	93.44
	200	92.73	93.39	93.34	94.15	93.10	93.78
	225	92.91	93.51	93.18	93.96	93.61	94.11
	250	93.34	94.01	93.84	94.39	92.92	93.67
	50	90.11	92.25	90.74	92.78	90.52	92.90
	75	91.22	92.70	92.57	94.08	91.94	93.44
	100	92.01	93.62	93.09	94.37	92.49	93.84
	125	92.39	93.45	93.28	94.64	93.59	94.56
100%	150	92.77	94.10	93.47	94.59	93.28	94.41
	175	93.22	94.35	93.66	94.84	94.00	94.93
	200	93.88	94.73	93.79	94.78	93.90	94.84
	225	93.74	94.62	94.25	95.07	93.69	94.70
	250	93.61	94.56	93.94	94.95	94.17	95.12

Table 7. Table 7 : Simulated coverage probabilities of confidence bands for F 1 subscript 𝐹 1 F_{1} in per cent where p = 0 𝑝 0 p=0

	Hall-Wellner						Equal Precision
	$N(0,1)$		$Poi\left(1\right)-1$		weird		$N(0,1)$		$Poi\left(1\right)-1$		weird
$n$	old	new	old	new	old	new	old	new	old	new	old	new
50	87.61	88.64	88.80	90.03	88.72	89.75	87.08	88.06	89.55	90.36	89.43	90.48
75	90.15	90.75	90.71	91.46	90.74	91.48	89.36	90.21	91.79	92.33	91.62	92.14
100	90.78	91.31	91.13	91.58	91.77	92.23	90.74	91.24	93.20	93.43	92.43	92.72
125	91.72	92.07	92.02	92.31	91.58	92.23	91.35	91.77	92.96	93.37	93.19	93.54
150	92.29	92.54	92.60	92.85	92.90	93.21	91.33	91.52	93.64	93.85	93.51	93.76
175	92.27	92.53	92.48	92.83	92.86	93.17	92.45	92.64	93.91	94.14	94.06	94.39
200	92.67	92.85	92.89	93.02	93.18	93.34	92.87	93.04	93.68	93.80	93.47	93.83
225	92.98	93.05	92.60	92.71	93.26	93.68	92.69	92.93	93.85	94.12	94.16	94.22
250	92.78	92.99	92.93	93.14	93.14	93.32	93.20	93.43	93.68	93.81	94.25	94.30

Equations

$$\alpha (t)=\underset{\Delta t\downarrow 0}{\lim }\frac{1}{\Delta t}P(T\in [t,t+\Delta t]|T\ge t)$$

$$N_{ji}(t)$$

$$\text{and}{Y}_i(t)$$

$${\Lambda }_1(t)=\sum _{i=1}^n{\Lambda }_{1i}(t)={\int }_0^{t}Y(u)\mathrm{d}{A}_1(u)=\sum _{i=1}^n{\int }_0^t{Y}_i(u)\mathrm{d}{A}_1(u),$$

$$⟨M_{1i}⟩(t)$$

$$\text{and}[M_{1i}](t)$$

$$\displaystyle[\sqrt{n}(\widehat{A}_{1}-A_{1}^{*})](t)=n\Big{[}\int_{0}^{\cdot}\frac{J\mathrm{d}M_{1}}{Y}\Big{]}(t)=n\int_{0}^{t}\frac{J}{Y^{2}}\mathrm{d}[M_{1}]=n\sum_{i=1}^{n}\int_{0}^{t}\frac{J}{Y^{2}}\mathrm{d}[M_{1i}],$$

$$[\sqrt{n}({\hat{A}}_1-A_1^{\ast })](t)$$

$$=n{\int }_0^t\frac{J}{Y^2}\mathrm{d}{N}_1-2n{\int }_0^t\frac{J}{Y^2}\Delta A_1\mathrm{d}{N}_1+n{\int }_0^t\frac{J}{Y}\Delta A_1\mathrm{d}{A}_1.$$

$${\hat{\sigma }}_1^2(t)=n{\int }_0^t\frac{J}{Y^2}\mathrm{d}{N}_1-2n{\int }_0^t\frac{J}{Y^2}\frac{\Delta N_1}{Y}\mathrm{d}{N}_1+n{\int }_0^t\frac{J}{Y}\frac{\Delta N_1}{Y}\frac{\mathrm{d}{N}_1}{Y}=n{\int }_0^t\frac{Y-\Delta N_1}{Y^3}\mathrm{d}{N}_1.$$

$${\hat{\tilde{W}}}_1(t)=\sqrt{n}\sum _{i=1}^n{\xi }_i{\int }_0^t\frac{\mathrm{d}{N}_{1i}}{Y}$$

$$t⟼⟨{\hat{\tilde{W}}}_1⟩(t)$$

$$\text{and}t⟼[{\hat{\tilde{W}}}_1](t)$$

$$\mathrm{d}{M}_{1i}(u)\text{by}{\xi }_i\sqrt{\frac{Y(u)-\Delta N_1(u)}{Y(u)}}\mathrm{d}{N}_{1i}(u),$$

$${\hat{W}}_1(t)=\sqrt{n}\sum _{i=1}^n{\xi }_i{\int }_0^t\sqrt{\frac{Y(u)-\Delta N_1(u)}{Y(u)}}\frac{\mathrm{d}{N}_{1i}(u)}{Y(u)}.$$

$${\mathcal{F}}_t=\sigma ({\xi }_i{N}_{1i}(u),N_{1i}(v),Y_i(v):u\in [0,t],v\in [0,K],i=1,\dots ,n).$$

$${\hat{\sigma }}_1^2:t⟼⟨{\hat{W}}_1⟩(t)=$$

$$\text{and}{\hat{\hat{\sigma }}}_1^2:t⟼[{\hat{W}}_1](t)=$$

$${\hat{W}}_1\stackrel{d}{⟶}{U}_1\sim \text{Gauss}(0,{\sigma }_1^2)\text{in outer probability}$$

$$\hat{S}(t)=\prod _{0\le u\le t}(1-\hat{A}(\mathrm{d}u)),t\in [0,K].$$

$$W_S(t)=\sqrt{n}(\hat{S}(t)-S(t))=S(t){\int }_0^t\frac{W(\mathrm{d}u)}{1-\Delta A(u)}+o_p(1).$$

$${\hat{W}}_S=\hat{S}{\int }_0^{\cdot }\frac{\hat{W}(\mathrm{d}u)}{1-\Delta \hat{A}(u)}$$

$$(s,t)\mapsto S(s)S(t){\int }_0^t\frac{\mathrm{d}A(u)}{(1-\Delta A(u))\bar{H}(u)};$$

$$⟨M_{1i},M_{2i}⟩(t)$$

$$[M_{1i},M_{2i}](t)$$

$$(W_1,W_2,\dots ,W_k)\stackrel{d}{⟶}(U_1,U_2,\dots ,U_k),$$

$$t\mapsto {\sigma }_j^2(t)={\int }_0^t\frac{1-\Delta A_j(u)}{\bar{H}(u)}\mathrm{d}{A}_j(u),j=1,2,\dots ,k$$

$$(s,t)\mapsto cov(U_j(s),U_{\ell }(t))=-{\int }_0^{s\wedge t}\frac{\Delta A_{\ell }(u)}{\bar{H}(u)}\mathrm{d}{A}_j(u)=:{\sigma }_{j\ell }(s\wedge t).$$

$${\hat{W}}_j(t)$$

$$\displaystyle+\frac{1}{\sqrt{2}}\sum_{\ell=1}^{k}sign(\ell-j)\Big{[}\sqrt{n}\sum_{i=1}^{n}\xi_{j\ell i}\int_{0}^{t}\sqrt{\frac{\Delta N_{j}(u)}{Y(u)}}\frac{\mathrm{d}N_{\ell i}(u)}{Y(u)}+\sqrt{n}\sum_{i=1}^{n}\xi_{\ell ji}\int_{0}^{t}\sqrt{\frac{\Delta N_{\ell}(u)}{Y(u)}}\frac{\mathrm{d}N_{ji}(u)}{Y(u)}\Big{]},$$

$$({\hat{W}}_1,{\hat{W}}_2,\dots ,{\hat{W}}_k)\stackrel{d}{⟶}(U_1,U_2,\dots ,U_k)\text{in outer probability},$$

$${\hat{\sigma }}_j^2(t)$$

$${\hat{\sigma }}_{j\ell }(t)$$

$${\hat{\hat{\sigma }}}_j^2(\cdot )$$

Full text

A Discontinuity Adjustment

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

March 8, 2024

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; Nelson-Aalen 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

[TABLE]

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\in\mathbb{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,

[TABLE]

respectively, with $j=1,\dots,k$ and $i=1,\dots,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}(t)=E(N_{11}(t))$ and $\bar{H}(t)=E(Y_{1}(t))$. Now, the sum $N_{1}(t)=\sum_{i=1}^{n}N_{1i}(t)$ has the compensator

[TABLE]

where $A_{1}(t)=\int_{0}^{t}\bar{H}(u)^{-1}\mathrm{d}H^{uc}_{1}(u)$ is the cumulative hazard function for a type 1 event. Therefore, $(M_{1}(t)=N_{1}(t)-\Lambda_{1}(t))_{t}$ is a square-integrable martingale. The Nelson-Aalen estimator for $A_{1}(t)$ is defined as the counting process integral $\widehat{A}_{1}(t)=\int_{0}^{t}Y^{-1}\mathrm{d}N_{1}$, where we let $0/0=0$. The limit covariance function of the normalized process $W_{1}=\sqrt{n}(\widehat{A}_{1}-A_{1})$ is determined by the predictable and optional variation processes of the square-integrable martingales $M_{1i}:t\mapsto N_{1i}(t)-\Lambda_{1i}(t)$. In the possibly discontinuous case, these are given by

[TABLE]

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\mapsto\int_{0}^{t}J(u)A_{1}(\mathrm{d}u)/Y(u)$ denote the compensator of the Nelson-Aalen estimator $\widehat{A}_{1}$, where $J(u)=1\{Y(u)>0\}$. By the usual martingale theory, we have the optional covariation process

[TABLE]

where the last equality follows from the orthogonality of all involved single martingales. Inserting the above representation for the optional variation process yields

[TABLE]

In a last step, we replace the unknown increments of the cumulative hazard function $\Delta A_{1}(u)$ with Nelson-Aalen increments $\Delta\widehat{A}_{1}(u)=Y(u)^{-1}\Delta N_{1}(u)$ and obtain the Greenwood-type variance estimator

[TABLE]

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}(K)<\infty$. A simple application of Glivenko-Cantelli theorems in combination with the continuous mapping theorem show that this Greenwood-type estimator is uniformly consistent for $\sigma_{1}^{2}(t)=\int_{0}^{t}(1-\Delta A_{1}(u))\bar{H}(u)^{-1}\mathrm{d}A_{1}(u)$ on $[0,K]$. 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 $\Delta A_{1}$ in the limit variance is crucially different from the continuous case in which the limit variance function reduces to $\tilde{\sigma}^{2}_{1}:t\mapsto\int_{0}^{t}\bar{H}^{-1}\mathrm{d}A_{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. 2013 or Bluhmki et al. 2016 for resampling more general counting process-based Nelson-Aalen estimators): First, the (independent) martingale increments $\mathrm{d}M_{1i}(u)$ in the martingale representation of $\sqrt{n}(\widehat{A}_{1}-A_{1})$ are replaced by independently weighted counting process increments, i.e., by $\xi_{i}\mathrm{d}N_{1i}(u)$. Here, the wild bootstrap weights $\xi_{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

[TABLE]

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

[TABLE]

When utilized as variance estimators, both variation processes in the previous display are uniformly consistent for $\tilde{\sigma}^{2}$. Unfortunately, this is even true in the discontinuous case where $\sigma^{2}\neq\tilde{\sigma}^{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\rightarrow\infty$ 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 $\mathrm{d}M_{1i}(u)$ should not be replaced by $\xi_{i}\mathrm{d}N_{1i}(u)$ but rather by something which – considered again as a martingale – reproduces the correct covariance structure. Therefore, we replace

[TABLE]

which results in the following wild bootstrap resampling version of the normalized Nelson-Aalen estimator:

[TABLE]

In a way similar to Bluhmki et al. (2016), one can show that $(\widehat{W}_{1}(t))_{t\in[0,K]}$ is a martingale with respect to the filtration $(\mathcal{F}_{t})_{t\in[0,K]}$ given by

[TABLE]

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

[TABLE]

It is well-known that the Greenwood-type variance estimator $\widehat{\sigma}_{1}^{2}$ is uniformly consistent for $\sigma^{2}_{1}$. Assuming the existence of the fourth moments of $\xi_{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}=\{t\in[0,K]:\Delta A_{1}(t)>0\}$ the subset of time points for which ties among the type $1$ events are possible. Similarly, we define $D_{2},\dots,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.

Condition 2.1**.**

$D=\bigcup_{j=1}^{k}D_{j}$* 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 $\mathcal{F}_{0}$ and as $n\rightarrow\infty$, we have the following conditional weak convergence

[TABLE]

on the càdlàg function space $D[0,K]$ equipped with the supremum distance topology, where $U_{1}$ is a Gaussian zero-mean martingale with variance function $t\mapsto\sigma^{2}_{1}(t)$. 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.

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 $S(t)=P(T>t)$ is

[TABLE]

It exhibits the martingale representation

[TABLE]

Thus, the discontinuity-adjusted wild bootstrapped normalized Kaplan-Meier estimator is

[TABLE]

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

[TABLE]

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\in\mathbb{N}$ competing risks, i.e. to $M_{1i},\dots,M_{ki}$, $i=1,\dots,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

[TABLE]

cf. the derivations in Section II.4 in Andersen et al. (1993). Define $W_{j}=\sqrt{n}(\widehat{A}_{j}-A_{j})$, $j=1,\dots,k$. The above variation processes are strong evidence that the normalized Nelson-Aalen estimators $W_{1},W_{2},\dots,W_{k}$ are not asymptotically independent anymore in the presence of ties. This is indeed the case:

Theorem 3.1**.**

As $n\rightarrow\infty$, we have on the product càdlàg function space $D^{k}[0,K]$, equipped with the $\max$-$\sup$-norm, that

[TABLE]

where $U_{1},U_{2},\dots,U_{k}$ are zero-mean Gaussian martingales with variance functions

[TABLE]

and covariance functions (for $j\neq\ell$)

[TABLE]

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 $\xi_{j\ell i},j,\ell=1,\dots,k,\ i=1,\dots,n$, be i.i.d. random variables with $E(\xi_{111})=0$ and $E(\xi^{2}_{111})=1$, which are also independent of the data. Denote by $N=\sum_{j=1}^{k}N_{j}$ the number of all kinds of events. Then, the single components of the wild bootstrap version of the multivariate Nelson-Aalen estimator $(W_{1},\dots,W_{k})$ are given by

[TABLE]

where $sign(x)=1\{x>0\}-1\{x<0\}$ 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 $\mathcal{F}_{0}$ and as $n\rightarrow\infty$, we have the following conditional weak convergence on the product càdlàg function space $D^{k}[0,K]$, equipped with the $\max$-$\sup$-norm:

[TABLE]

where $(U_{1},U_{2},\dots,U_{k})$ 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:

[TABLE]

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

[TABLE]

are the optional process-type (co)variance estimators motivated from the wild bootstrap martingale properties. Assume that all $\xi_{111}$ have finite fourth moments. By applications of Glivenko-Cantelli theorems in combination with the continuous mapping theorem, the Greenwood-type (co)variance estimators $\widehat{\sigma}_{j}$ and $\widehat{\sigma}_{j\ell}$ are shown to be uniformly consistent for $\sigma_{j}^{2}$ and $\sigma_{j\ell}$, 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 $E(\xi_{111}^{4})<\infty$: In points of continuity of all cumulative hazard functions, i.e. on $[0,K]\setminus D$, 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 $\widehat{\widehat{\sigma}}^{2}_{j}$ by $\widehat{\sigma}^{2}_{j}$ and apply the conditional Chebyshev inequality (given $\mathcal{F}_{0}$) in order to show the negligibility of the differences $\widehat{\widehat{\sigma}}^{2}_{j}-{\widehat{\sigma}}^{2}_{j}$ in probability. The last argument can be repeated for the covariance estimators. A final application of the continuous mapping theorem yields

[TABLE]

for all $j\neq\ell$ as $n\rightarrow\infty$.

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

Denote the cumulative incidence functions by $F_{j}(t),j=1,\dots,k$, which specify the probabilities to die due to cause $j$ during the time interval $[0,t]$. 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-U_{1}$, $A-A_{1}$, and $1-S-F_{1}$, respectively. Here, $A$ denotes the all-cause 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 $\widehat{F}_{1}(t)=\int_{0}^{t}\widehat{S}(u-)\mathrm{d}\widehat{A}_{1}(u)$ for the cumulative incidence function $F_{1}(t)=\int_{0}^{t}S(u-)\mathrm{d}A_{1}(u)$:

Theorem 4.1**.**

As $n\rightarrow\infty$,

[TABLE]

where $U_{F_{1}}$ is a zero-mean Gaussian process with covariance function

[TABLE]

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 $\sqrt{n}(\widehat{F}_{1}-F_{1})$ is given by

[TABLE]

where $\widehat{W}_{1}$ and $\widehat{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 $\mathcal{F}_{0}$ and as $n\rightarrow\infty$, we have the following weak convergence on the càdlàg function space $D[0,K]$, equipped with the $\sup$-norm:

[TABLE]

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 $\xi_{j\ell i}+1\sim Bin(Y(X_{i}),\max(1,Y(X_{i}))^{-1})$, 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 $\sigma^{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:

[TABLE]

Similarly, $\widehat{\widehat{\sigma}}^{2}_{F_{1}}$ is obtained by replacing all estimators $\widehat{\sigma}_{j}$ and $\widehat{\sigma}_{j\ell}$, $j\neq\ell$, by their wild bootstrap counterparts $\widehat{\widehat{\sigma}}_{j}$ and $\widehat{\widehat{\sigma}}_{j\ell}$, respectively. Their uniform (conditional) consistencies for $\sigma^{2}_{F_{1}}$ on $[0,K]$ 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), time-simultaneous confidence bands for $F_{1}$ can be deduced. In particular, let $\phi(s)=\log(-\log(1-s))$ be a transformation applied to $F_{1}$ in order to ensure band boundaries between 0 and 1 and let $g_{1}(s)=\log(1-\widehat{F}_{1}(s))/\widehat{\rho}(s)$ and $g_{2}(s)=\log(1-\widehat{F}_{1}(s))/(1+\widehat{\rho}^{2}(s))$ be weight functions leading to the usual equal precision and Hall-Wellner bands, respectively; see Andersen et al. (1993). Here $\widehat{\rho}^{2}(s)=\widehat{\sigma}_{F_{1}}^{2}(s)/(1-\widehat{F}_{1}(s))^{2}$. Let $\widehat{g}_{1}$ and $\widehat{g}_{2}$ be their wild bootstrap counterparts, i.e. the variance estimates $\widehat{\sigma}_{F_{1}}^{2}$ are replaced by $\widehat{\widehat{\sigma}}_{F_{1}}^{2}$. The confidence bands for $F_{1}$ are then derived from the asymptotics of the supremum distance $Z_{1j}=\sup_{u\in[t_{1},t_{2}]}|\sqrt{n}g_{j}(u)(\phi(\widehat{F}_{1}(u))-\phi(F_{1}(u)))|$ and its wild bootstrap counterpart $\widehat{Z}_{1j}=\sup_{u\in[t_{1},t_{2}]}|\widehat{g}_{j}(u)\phi^{\prime}(\widehat{F}_{1}(u))\widehat{W}_{F_{1}}(u)|$, where $[t_{1},t_{2}]\subset[0,K]$ and $j\in\{1,2\}$. Let $q_{0.95,j}$ be the conditional $95\%$ quantile of $\widehat{Z}_{1j}$ given the data. The resulting asymptotic $95\%$ confidence bands are $1-(1-\widehat{F}_{1}(s))^{\exp(\pm n^{-1/2}q_{0.95,j}/g_{j}(s))},\ s\in[t_{1},t_{2}],\ 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 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 $\alpha_{1}(t)=\exp(-t)$ and $\alpha_{2}(t)=1-\exp(-t)$ which yield the cumulative function of the first risk $F_{1}(t)=0.5(1-\exp(-2t))$. 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 $\{0,\frac{1}{k},\frac{2}{k},\dots\}$, where $k\in\{5,10,20\}$, and the discretization probabilities to be $p\in\{0,0.25,0.5,0.75,1\}$. The resulting theoretic cumulative incidence functions

[TABLE]

are presented in Figure 1, with the exception of the continuous function $F_{1}$. Here $[s]$ denotes the integer closest to $s\in\mathbb{R}$. For simulating data, which have the desired cumulative incidence function $F_{1}^{p,k}$, it is mandatory to first round the event times $T_{i}$, and then generate the event types $\varepsilon_{i}$ in a second step, according to the formula

[TABLE]

where $S:t\mapsto\exp(-t)$ 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 $i$th survival time is discretized, then we discretize the $i$th 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 $[0.25,0.75]$. 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\in\{5,10\}$, 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\in\{5,10\}$, 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_{1}^{p,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 variate-based 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.

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 $[5,55]$, 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\in\{50,75\}$ 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.

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).

Appendix A Detailed Derivation of the Asymptotic Variance of the Nelson-Aalen Estimator

Define $H^{uc}(t)=\sum_{j=1}^{k}EN_{j1}(t)$ as the probability that an uncensored event due to any cause occurs until time $t$. According to van der Vaart and Wellner (1996), p. 383f, we have $\sqrt{n}(\widehat{A}-A)\stackrel{{\scriptstyle d}}{{\longrightarrow}}\int_{0}^{\cdot}M^{uc}(\mathrm{d}u)/\bar{H}(u)$, where $M^{uc}(t)=G^{uc}(t)-\int_{0}^{t}\bar{G}(u)\mathrm{d}A(u)$ is a zero-mean Gaussian martingale. Its variance function is determined by

[TABLE]

[TABLE]

[TABLE]

Note that $A(t)=\int_{0}^{t}H^{uc}(\mathrm{d}u)/\bar{H}(u)$. Thus, for $s\leq t$, the covariance function of $M^{uc}$ at $(s,t)$ is

[TABLE]

We conclude, as in van der Vaart and Wellner (1996), that

[TABLE]

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 cause-specific cumulative hazard functions that

[TABLE]

Appendix B Consistency of the Wild Bootstrap for the Univariate Nelson-Aalen Estimator

Proof.

Without loss of generality, assume that $0,K\in D_{1}$ for simplifying notation. Write $D_{1}=\{d_{0},d_{1},\dots,d_{J}\}$ with the natural ordering $d_{j}<d_{j+1}$ for all $j=1,\dots,J$. Then $[0,K]\setminus D_{1}=\bigcup\limits_{j=1}^{J}(d_{j-1},d_{j})$. As argued in Bluhmki et al. (2016), it is now straightforward to show that each process $(\widehat{W}_{1}(t)-\widehat{W}_{1}(d_{j-1}))_{t}$ on each interval $[d_{j-1},d_{j})$, $j=1,\dots,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 $\widehat{W}_{1}(d_{j})-\widehat{W}_{1}(d_{j-1}):=\lim_{t\uparrow d_{j}}\widehat{W}_{1}(t)-\widehat{W}_{1}(d_{j-1})$, and this procedure introduces square-integrable martingales on the whole $[d_{j-1},d_{j}]$.

First, we notice that the conditional weak convergence in probability of the processes $(\widehat{W}_{1}(t)-\widehat{W}_{1}(d_{j-1}))_{t}$ on each interval $[d_{j-1},d_{j}]$, $j=1,\dots,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 $(\tilde{U}_{1j}(t))_{t\in[d_{j-1},d_{j}]}$, $j=1,\dots,J$. Due to the martingale extension above, Rebolledo’s limit theorem implies the almost sure continuity of $\tilde{U}_{1j}$ on each time interval. Furthermore, these are zero-mean processes with variance function $t\mapsto\sigma^{2}_{1}(t)-\sigma^{2}_{1}(d_{j-1})$.

Due to the continuity of the limit processes $\tilde{U}_{1j}$ on the intervals $[d_{j-1},d_{j}]$, 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 $\Delta\widehat{W}_{1}(d_{j})$ holds in probability by the already argued convergence of all finite-dimensional conditional distributions. Therefore, the independence of the (bootstrapped) Nelson-Aalen increments imply that, as $n\rightarrow\infty$, the conditional distribution of

[TABLE]

given $\mathcal{F}_{0}$ converges weakly in probability to the distribution of

[TABLE]

on the product Space $\mathbb{R}\times D[d_{0},d_{1}]\times\mathbb{R}\times\dots\times D[d_{J-1},d_{J}]\times\mathbb{R}$ equipped with the $\sup$-$\max$-norm. Here, all components are independent, and the normally distributed random variables $\tilde{U}_{1}(d_{j})$ have mean zero and variance $\Delta\sigma^{2}_{1}(d_{j})$, $j=1,\dots,J$.

Applying the functional, which is continuous with respect to the $\max$-$\sup$-norm,

[TABLE]

[TABLE]

to the previous limit theorem, the continuous mapping theorem implies that: Given $\mathcal{F}_{0}$ and as $n\rightarrow\infty$, the conditional distribution of

[TABLE]

converges weakly (on the product-function space $D[d_{0},d_{1}]\times\dots\times D[d_{J-1},d_{J}]\times\mathbb{R}$ equipped with the $\max$-$\sup$-norm) in probability to the distribution of

[TABLE]

Here the right boundary values are again considered as the left-hand limits given by the martingale extension theorem. The process $(U_{1}(t))_{t\in[0,K]}$ is a zero-mean Gaussian martingale with variance function $t\mapsto\sigma^{2}_{1}(t)$ with, in general, discontinuous sample paths.

Finally, we apply the continuous functional

[TABLE]

in order to obtain the desired conditional weak convergence for $(\widehat{W}_{1}(t))_{t\in[0,K]}$ in probability. ∎

Appendix C Derivation of the Asymptotic Covariance of Multiple Nelson-Aalen 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\rightarrow\infty$,

[TABLE]

by the continuous mapping theorem. The covariance function of $U$ is given by

[TABLE]

but on the other hand

[TABLE]

Solving for the unknown covariances on the right-hand side of the previous display, we obtain

[TABLE]

Due to symmetry and inductively, it follows that

[TABLE]

for $j\neq\ell,\ j,\ell=1,\dots,k,$ even in the situation of $k\in\mathbb{N}$ competing risks.

Appendix 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

[TABLE]

Further, consider arbitrary points of time $0\leq t_{1}\leq\dots\leq t_{m}\leq K$ and the vector

[TABLE]

Due to analogy, it is enough to calculate the entries corresponding to $(\widehat{W}_{1}(t_{1}),\widehat{W}_{1}(t_{2}))$ and $(\widehat{W}_{1}(t_{1}),\widehat{W}_{2}(t_{2}))$ of the limit of the matrix

[TABLE]

We start by calculating the entry for $a=1,b=2,j=\tilde{j}=1$, that is,

[TABLE]

Here, the last equality follows from $\Delta N_{1}\mathrm{d}N_{\ell}=\Delta N_{\ell}\mathrm{d}N_{1}$ for all $\ell$. By the Glivenko-Cantelli theorem in combination with the continuous mapping theorem, it follows that the quantity in the previous display converges to $\sigma_{1}^{2}(t_{1})$ in probability as $n\rightarrow\infty$.

Now, consider the entry of $\widehat{\Gamma}$ for $a=1,b=2,j=\tilde{j}=2$:

[TABLE]

By the same arguments as before, this is a consistent estimator for $\sigma_{12}(t_{1})$ as $n\rightarrow\infty$.

Appendix E Derivation of the Asymptotic Covariance Function of an Aalen-Johansen Estimator in the Competing Risks Set-Up

.

In order to derive the the covariance function of

[TABLE]

at any $(s,t)\in[0,K]^{2}$, we exemplarily calculate the covariance function of the first integral and the covariance function between both integrals. Hence, as the covariance function $(s,t)\mapsto\sigma^{2}_{1}(s\wedge t)$ of $U_{1}$ only increases along the diagonal,

[TABLE]

Furthermore, we similarly have for the covariance between both integrals that

[TABLE]

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:

[TABLE]