Confidence bands for multiplicative hazards models: Flexible resampling approaches.

We propose new resampling-based approaches to construct asymptotically valid time-simultaneous confidence bands for cumulative hazard functions in multistate Cox models. In particular, we exemplify the methodology in detail for the simple Cox model with time-dependent covariates, where the data may be subject to independent right-censoring or left-truncation. We use simulations to investigate their finite sample behavior. Finally, the methods are utilized to analyze two empirical examples with survival and competing risks data.

INTRODUCTION

Multiplier resampling (also known as the wild bootstrap) has evolved as a state‐of‐the‐art choice for inference on cumulative incidences or hazards in nonparametric and semiparametric survival and multistate models in event history analysis; see Lin et al. (1993) and Lin et al. (1994) for Cox models and Lin (1997) for competing risks set‐ups. The basic idea is to mimic martingale representations of the estimators by replacing nonobservable martingale residuals with randomly weighted counting processes. This approach has been extended in various directions, allowing for a greater flexibility (Beyersmann et al., 2013; Dobler & Pauly, 2014; Dobler et al., 2017) and multiple, possibly recurrent, states (Bluhmki et al., 2018, 2019). We will here focus on semiparametric regression models. The most used regression model in survival analysis is Cox's proportional hazard model (Cox, 1972). It is highly useful to predict survival probabilities based on covariates. Such survival functions are typically provided with point‐wise confidence intervals, and this is implemented in all major software packages. However, when interest is in the survival function as a whole, it is preferable to report it together with uniform confidence intervals. These so‐called confidence bands describe the uncertainty of the whole survival function estimation. But this is infrequently done in practice because there are few programs that construct such uniform bands. In addition, apart from only few exceptions such as Lin (1997), systematic evaluations of finite sample results with demonstrations of the bands’ performance are rarely available in the literature. In the present article we provide such results in Cox models and we investigate various new resampling bands that exhibit a reliable performance even for small sample sizes.

Our main achievements are the introduction of resampling strategies that jointly mimic the unknown distribution of baseline and parameter estimators in Cox models and corresponding multistate versions. We provide theoretical justification for the resampling approaches based on martingale arguments, thereby allowing for different mechanisms of incomplete observations and simplifying many arguments.

There exist many resampling approaches for the Cox model with independent right‐censoring. One approach is based on the martingale representations of the Breslow estimator for the cumulative hazard and the parameter estimator, where all unobservable martingale residuals are replaced with reweighted counting processes (e.g., Lin et al., 1994). Another possibilty that works more generally, for example for rate estimation, is to replace the martingale instead with (e.g., Spiekerman & Lin, 1998). This exploits the model's semiparametric structure more actively.

We will also use a more natural approach, which starts one step earlier by rewriting the score equations for the baseline function and the Euclidean parameter: we apply the multiplier techniques from above to martingale representations of the score equations. This leads to new equations, which are solved by quantities depending on the multipliers. Hence, paralleling the same steps as for the original estimators, their resampling counterparts are obtained in an elementary way. Throughout the article, we follow the flexible approach of Beyersmann et al. (2013) and allow for general multiplier distributions.

For ease of presentation, we exemplify the new methodology mainly for the Cox model in a survival setting but we also discuss extensions to more general multistate and other regression models. We are able to base our theoretical derivations for the resampling approaches on martingale arguments, which are novel in the semiparametric context. Thus, intricate derivations for verifying conditional tightness, that is, tightness of the resampling process conditioned on the observed data, are no longer required and are handled conveniently by the martingale central limit theorems. The martingale arguments for the multiplier bootstrap are motivated by those for resampling nonparametric Aalen‐Johansen estimators (Bluhmki et al., 2018). In particular, we prove similar martingale properties of the resampling counterparts as for the original estimators. Moreover, mirroring the martingale structure in the bootstrap world allows for a simple interpretation and easy incorporation of missingness mechanisms, for example, independent right‐censoring or left‐truncation. To conclude, our findings allow for a wide range of applications some of which will be discussed in more detail in future papers.

The paper is organized as follows: Section 2 sketches how estimators in Cox models are found and it lists the technical conditions required for all of their large sample properties. Section 3 contains a derivation of all considered multiplier bootstrap procedures and also theoretical statements about their validity. In addition, we discuss several important extensions to more general multistate or other regression models. In Section 4, we compare the resampling procedures in an extensive simulation study. Section 5 has a brief demonstration of the methodology in a survival setting, where interest is on constructing confidence bands for the survival function for patients with acute myocardial infarction, and another demonstration of confidence bands for cumulative incidence functions for patients after a bone marrow transplantation. Finally, we discuss the findings in Section 6.

The online supporting information contains all proofs whose novelty lies in the considerable simplification of the technical arguments by transferring martingale methods to the resampling procedures in the present semiparametric framework. It also contains a second small simulation study, algorithmic descriptions of the multiplier bootstraps, and links to GitHub pages with example R code for user‐friendly applicability of our procedures.

JOINT LARGE SAMPLE PROPERTIES IN THE COX MODEL

We consider multiplicative Cox (or proportional hazards) models given by the intensity of the counting process of subject , possibly subject to independent left‐truncation and/or right‐censoring, given its possibly time‐dependent ‐dimensional vector of predictable covariates . That is, . Here, counts the occurence of a particular event for subject not later than treatment time is the corresponding at‐risk indicator at that time, is the baseline hazard function, and is an unknown regression parameter (Cox, 1972). Let be a terminal evaluation time on the treatment time‐scale. Throughout we assume all to be contained in a bounded set and that the cumulative baseline hazard is finite at . Standard arguments yield the Breslow estimator of and the maximum partial likelihood parameter estimator of . To illustrate this, we simplify derivations of Scheike and Zhang (2002) for Cox‐Aalen models to the Cox model Equation (1): the score equation for is and solved by . We will make use of the definitions where and for any vector , and is the indicator that any individual is under risk shortly before . For simplicity, the notion of is usually suppressed. Replacing for in the score equation for , that is, and defining , we obtain a solvable score equation for : Through its solution we obtain the Breslow estimator of . The central limit theorem for stated below requires the negative Jacobian of with respect to , that is, , where and . Due to the assumed boundedness of the covariates, it follows from Theorems VII.2.2 (p. 498) and VII.2.3 (p. 503f.) in Andersen et al. (1993) that and are asymptotically Gaussian provided the following regularity conditions hold; see Condition VII.2.1 in Andersen et al. (1993, p. 497). Let denote convergence in probability as and the gradient operator with respect to .

There exists a neighborhood of and functions , and such that for each :

, where denotes the Euclidean norm;

is a continuous function of uniformly in and bounded on ;

is bounded away from zero on ;

and for ;

is positive definite, where .

Throughout, we assume Condition 1 to be satisfied. Note that (a) and (b) immediately imply for each as long as . The proofs of the aforementioned theorems in Andersen et al. (1993) make use of the following asymptotic representations of the estimators: These representations motivate the first bootstrap approach in Section 3. Here, is the probability limit of , and is a square‐integrable martingale in ; cf. Andersen et al. (1993, Section VII.2.2; p. 496ff) .

MULTIPLIER RESAMPLING APPROACHES AND MAIN THEOREMS

While inference about can be based on the asymptotic normality of its estimator (e.g., Martinussen & Scheike, 2006; p. 184ff.), the complicated limit process of the normalized Breslow estimator hinders straightforward time‐simultaneous inference about or functionals thereof, for example, the survival function. As a remedy we will investigate two resampling strategies, both of which make use of independent i.i.d. zero mean and unit variance random variables , which are henceforth called multipliers. Since and are dependent, we have to ensure that the bootstrap procedures respect their dependence structure.

The ‘Classical’ multiplier bootstrap

The first method is inspired by the resampling procedures of Lin et al. (1994) and Spiekerman and Lin (1998). They use standard normal multipliers, which is motivated from the asymptotic martingale representation of the normalized estimators. The idea is to replace the martingale increments with (Lin et al., 1994) or with , using estimated martingale increments (Spiekerman & Lin, 1998). We extend these approaches in several ways: our approaches are not restricted to normal multipliers ; cf. Beyersmann et al. (2013) where centered Poisson multipliers exhibited a better finite sample behavior in a nonparametric setting. Also, to be in line with the usual guidelines for bootstrap‐based inference (Hall & Wilson, 1991), we provide and analyze suitable bootstrap versions of the covariance estimators. Finally, in Section 3.4 an extension to more general models is outlined.

We exemplify the general resampling idea for by introducing multiplier bootstrap versions of the score equation defining vector and its negative Jacobian . Here and below, all objects with asterisks refer to bootstrap quantities. Following the asymptotic representation Equations (5) and (6), we find their following multiplier bootstrap counterparts:

Note that Equations (9) and (10) are essentially definitions of and , respectively, and that the terms in Equations (5) and (6) have been dropped. Alternatively to this multiplier bootstrap choice, the Spiekerman and Lin (1998)‐type martingale increment estimates may replace in Equations (7) and (10). A bootstrap‐type covariance estimate similar to Equation (8) has been suggested by Dobler and Pauly (2014) in a nonparametric competing risks context. Here, it is additionally motivated from martingale arguments: defining and as in Equations (7) and (8) but replacing , we found that is the optional variation process of the square‐integrable martingale ; see the supporting information for details.

The next subsection offers an alternative resampling strategy, which does not ignore the terms in Equations (5) and (6). Instead, they are naturally incorporated into the resampling step.

Bootstrapping the score equations

The idea of the second approach is to replace martingale representations of score equations with their multiplier counterparts. To this end, we parallel the classic approach of jointly solving two score equations and first expand the score equation in Equation (2) to . A multiplier bootstrap counterpart thereof is now given by replacing with with , and with : Keeping fixed, its “solution” for is Next, consider a martingale representation of the score Equation (3) for : This motivates a multiplier resampling version given by Inserting for eventually yields the final multiplier bootstrap score equation The last equality is due to . Hence, coincides with formula Equation (7). Define as the solution to Equation (13). Similarly as in the proof of Theorem VII.2.1 in Andersen et al. (1993, p. 497f.), it can be shown that the probability of the existence of tends to one and that converges to zero in probability as ; see also the proof of Theorem 1 below for similar arguments.

Finally, we combine Equation (12) with to find a resampling version of the Breslow estimator, denoted . After centering it at and multiplying it with , it equals A Taylor expansion around of the first term on the right‐hand side reveals the striking similarity to decomposition Equation (10). However, the current multiplier bootstrap approach does not ignore the term resulting from the Taylor expansion. Another nice property of this “estimating equation” approach is the similar treatment for bootstrap and original estimator, which is in line with general recommendations for constructing resampling algorithms (Efron & Tibshirani, 1994). We are going to use the mean value theorem (cf. Feng et al., 2013) to analyze the asymptotic behavior of . To this end, we introduce a modified Jacobian operator. For a continuously differentiable function and a matrix we define The mean value theorem yields that , where every column of is on the line segment between and .

Consistency and confidence bands for the cumulative hazard

The validity proof of both resampling strategies from Sections 3.1 and 3.2 is based on

Under Conditions 1(a)‐(e) it holds that, as if the resampling is done via method Equation (9) or (13).

,

if each column of converges in probability to ,

given all observations, in probability, i.e. it is asymptotically multivariately normally distributed,

The next theorem shows that both bootstraps based on the approach (Lin et al., 1994) are asymptotically valid. Therein, denotes a distance that metrizes weak convergence on , e.g. the Prohorov distance (Dudley, 2002, p. 309ff.), and and are, respectively, the unconditional and conditional law of a random variable .

Under Condition 1 it holds for both resampling strategies Equation (9) or (13) that the asymptotic distributions of and coincide, that is, as , where and .

The asymptotic variance function of (and thus also of ) can be found in Andersen et al. (1993), Corollary VII.2.4, p. 505, where a consistent estimator is given. In our simulation study in Section 4, we chose the multiplier bootstrap counterpart of to be the empirical variance of the bootstrap realizations of .

We use the theorem to construct time‐simultaneous confidence bands for on fixed intervals . In particular, we obtain results similar to those of Lin et al. (1994): denoting by a continuously differentiable and strictly increasing or decreasing function, asymptotic level confidence bands for on are given by where is a possibly random weight function. Typical choices are in case of the transformation , and for the transformation . The resulting confidence bands correspond to the so‐called equal precision (for or ) and Hall‐Wellner bands (for or ), respectively. Let and be the multiplier bootstrap analogs of and , respectively, , that is, we utilized instead of and the corresponding empirical variances instead of . We chose to be the quantile of the conditional distribution of , and the naïve choice for is the corresponding quantile of . However, the latter resulted in numerical instabilities, and we preferred the asymptotically equivalent choice .

It follows from Theorem 1 that all confidence bands are valid for large sample sizes. To additionally assess their small sample properties, we compare them in Monte‐Carlo simulations in Section 4. There, we also analyze the analogous behavior of the resampling approaches based on . It should be noted that the confidence bands obtained from the ‘classical’ multiplier bootstrap from Section 3.1 with standard normal multipliers are not equal to those proposed in Lin et al. (1994). The key difference is that we additionally use a resampling version of the weight functions (particularly of the variance estimator) to be in line with the usual guidelines for bootstrap‐based inference (Hall & Wilson, 1991). We additionally refer to Allignol et al. (2018), where a similar recommendation was made.

Extensions to more general models and more on inference

A related multiplier bootstrap approach carries over to more general models in multistate set‐ups. In particular, as long as the counting process martingale can be mimicked with the help of bootstrap multipliers, the asymptotics of the resampled estimators result in almost the same way as for the original estimators. Thus, the above methodology can straightforwardly be extended to multistate models with states and multiplicative intensity processes for each transition , where . Different to above, this model allows for an arbitrary number of transitions between different states. Following Bluhmki et al. (2018), the above multiplier bootstrap approach remains applicable: instead of the previously used multipliers , more general white noise processes with zero mean and unit variance are required. Randomly weighting the counting process increments leads to . The evaluation occurs at a finite number of times, so no existence problems arise for the white noise processes, and the implementation is straightforward.

The martingale arguments still hold, so that the multiplier bootstrap mimics the joint stochastic behavior of the parameter and multivariate hazard function estimators. Indeed, according to Bluhmki et al. (2018), independent white noise processes and for transitions give rise to orthogonal square‐integrable martingales in and , with respect to the filtration Here, and are ‐predictable random functions. The predictable variation process of the first martingale is , which nicely reflects the properties of the original estimators: the corresponding orthogonal counting process martingales and have predictable variation , and similarly for . In this sense, the bootstrap‐based predictable variation processes estimate the original predictable variation processes. Extending the proofs in the supporting information based on these findings, large‐sample properties of estimators in multistate problems transfer to their multiplier bootstrap versions in an obvious way, as long as the original estimators admit martingale representations. Similar arguments extend to Cox‐Aalen multiplicative‐additive intensities (Scheike & Zhang, 2002) or other models with martingale structures.

The incorporation of certain filtered (e.g., right‐censored) observations is again allowed and important inferential applications follow: confidence bands for cumulative transition hazards or incidence functions, and tests for null hypotheses formulated in terms of the parameters can be constructed. Here, bootstrap‐based versions of score or Wald‐type test statistics (Martinussen & Scheike, 2006, p. 185f.) might ensure a proper finite sample behavior. However, detailed evaluations of these applications require additional extensive simulations and further elaborations in future research. For simplicity, we continue with our focus on the Cox model Equation (1) and assess the impact of the proposed methods in simulations.

SIMULATION STUDY

To compare the performances of the various resampling approaches described in Section 3.3, we conducted a simulation study in which we covered situations of small to large sample sizes: . The generated data follow the Cox survival model with baseline hazard rate , normally distributed covariates with standard deviation , and regression parameter . Censoring times are standard exponentially distributed, truncated at . The time interval, along which 95% confidence bands for the cumulative baseline hazard function shall be constructed, was . Here we chose the start time of because “the approximations tend to be poor for close to 0” (Lin et al., 1994, p. 77). Another simulation study that evaluates the multiplier bootstrap for confidence bands for baseline cumulative incidence functions is given in Suporting Information.

As bootstrap multipliers , we considered the common choice , centered unit Poisson variables with unit skewness, and centered unit exponential variables with a skewness of 2. We simulated all confidence bands for that were introduced in Section 3.3, i.e. ‐ and nontransformed Hall‐Wellner and equal precision bands. We considered replacing the martingale increments with or , denoted in Table 1 as “” and “”, respectively, and also both kinds of resampling algorithms, the direct resampling method of Section 3.1 and the method of Section 3.2 in which the estimating equations were bootstrapped. We note that, apart from the specific variance estimation scheme, the bands resulting from the cox.aalen function in the R package timereg are essentially the same as the multiplier bootstrap from Section 3.1 with no transformation and standard normal multipliers and the results for these are therefore not shown. For each set‐up and type of band, we constructed 10,000 confidence bands, where each was based on 1,000 resampling iterations. The obtained empirical coverage probabilities given in Table 1 are the observed frequencies of constructed bands containing completely.

Table 1

Simulated coverage probabilities (in %, rounded) of various 95% confidence bands for the baseline cumulative hazard function and sample sizes

			Hall‐Wellner				Equal precision
Distribution of G_i		Resampling approach	Estimating equation		Direct resampling		Estimating equation		Direct resampling
	n		id	log	id	log	id	log	id	log
N(0,1)	100	dN	86.8	93.4	87.2	93.4	85.6	92.7	85.5	92.5
		dM	86.9	93.2	87.2	93.6	85.7	92.6	85.2	92.4
	200	dN	90.8	94.2	90.9	94.2	89.8	93.9	89.8	93.8
		dM	90.6	94.0	90.7	94.2	89.4	93.3	89.3	93.2
	400	dN	93.2	95.1	93.3	95.0	92.4	94.1	92.4	94.1
		dM	93.2	94.9	93.2	95.1	92.3	94.1	92.4	94.0
Exp(1)−1	100	dN	88.4	94.4	88.5	94.5	88.8	95.3	89.0	95.3
		dM	88.1	94.3	88.4	94.5	88.7	95.4	87.9	94.6
	200	dN	91.3	94.6	91.3	94.7	91.3	95.4	91.4	95.3
		dM	91.2	94.5	91.2	94.6	91.0	94.9	90.8	94.6
	400	dN	93.6	95.3	93.5	95.2	93.1	94.9	93.1	95.0
		dM	93.9	95.4	94.1	95.5	93.5	95.4	93.4	95.2
Poi(1)−1	100	dN	87.0	93.3	87.1	93.6	86.1	93.3	86.2	93.0
		dM	87.0	93.6	87.4	93.7	86.0	93.2	86.0	92.8
	200	dN	90.6	94.1	90.7	94.2	89.6	93.6	89.7	93.7
		dM	91.0	94.2	91.2	94.3	90.0	94.1	90.0	94.0
	400	dN	93.3	95.0	93.3	95.1	92.4	94.2	92.6	94.2
		dM	93.3	95.1	93.1	95.1	92.4	94.3	92.3	94.2

Distributions of the multipliers : standard normal (top panel), centered unit exponential (middle panel), centered unit Poisson (bottom panel)

When the sample size is , all methods give a reasonable performance. For smaller samples there are notable differences, and it seems that the log‐transform improves the bands’ performances considerably. Whether the bootstrap is based on or does not seem important, and in terms of computations it is considerably easier and faster to use the multipliers based on . Despite theoretical and practical advantages of Poisson variables over standard normal multipliers in nonparametric competing risks models (Dobler et al., 2017), this choice of bootstrap multipliers does not seem highly important here. However, centered exponential multipliers tend to perform the best. Also, the particular resampling method, be it the direct approach of Section 3.1 or the estimating equation approach of Section 3.2, does not seem to have a clearly positive or negative impact on the outcomes.

Apart from the coverage probabilities, we also computed the bands’ median widths for . Because of the marginal differences among multiple resampling approaches, Figure 1 only shows the median widths for the method with standard normal multipliers, based on resampling the estimating equation. Again, we let . In all panels the trend is the same: at early time points the equal precision bands are slightly wider than the Hall‐Wellner bands. But for later time points it is the opposite; and the widths of the Hall‐Wellner bands grow rapidly due to quickly increasing variances. Also, the log‐transformed bands are slightly wider than their untransformed counterparts. All in all, we recommend to use log‐transformed equal precision bands for because of their good coverage probabilities (log‐transformation) and their reasonable slimness.

Figure 1

Median widths of untransformed and log‐transformed confidence bands for sample sizes (top), (middle), (bottom panel). Solid and dashed: Hall‐Wellner, pointed and point‐dashed: equal precision. This figure appears in color in the electronic version of this article

Finally, we remark that our simulation results are only partially comparable to those of Lin et al. (1994): we construct confidence bands for , that is, for an individual with ; their bands are for survival curves for multiple covariates using also different transformantions.

Overall, however, the bands here and in Lin et al. (1994) seem to work well.

DATA EXAMPLES

In this section, we illustrate how the resampling approach can be used in two practical settings, based on data that are available in the timereg R‐package. The R‐code is given in the supporting information. First, we show how confidence bands can be used and constructed in a standard survival setting. The key point is that a simultaneous confidence band is the type of interval estimate most often of interest unless focus is on a particular survival probability at a specific time such as, for example, 5‐year survival. In addition, we find confidence intervals for derived measures, in particular, for the restricted mean life. The second example illustrates in a competing risks setting how simple it is to combine the multiplier bootstrap based on Cox models for cause‐specific hazards to construct confidence bands for the cumulative incidence.

Survival estimation

We consider the TRACE study (Jensen et al., 1997), where interest is on survival after acute myocardial infarction for 1878 consecutive patients included in the study. The data‐set is available in the timereg R‐package. Here for sake of illustration we focus interest on the covariates diabetes (1/0), sex and age. We depict the survival predictions with confidence bands for a male with average age (66.9 years) and with or without diabetes, as well as the standard 95% point‐wise confidence intervals for comparison.

We note that the hazard ratio related to diabetes is 1.82 with 95% confidence interval (1.50, 2.18). Thus reflecting that diabetes is a factor that leads to increased mortality. More interestingly, in regard to absolute mortality rates, the higher risk of diabetics is reflected in our estimated survival curves for males with average age and with diabetes (lower set of curves, with confidence bands (region) and point‐wise intervals (dotted lines) or without diabetes (upper set of curves). The regions are standard timereg‐bands based on direct resampling and via , and for comparison we also show the resampling‐based equal precision bands based on a log‐transform (broken lines). We note a slight difference between the different confidence bands. In addition we note that the bands are quite a bit wider than the point‐wise intervals. As the latter do not provide simultaneous coverage, the bands should be used to provide uncertainty about the entire survival curve as shown in Figure 2.

Figure 2

Survival function estimates for males with average age, with (fat broken line) or without diabetes (fat solid line), with 95 % confidence intervals (dotted lines), multiplier bootstrapped log‐transform 95 % confidence bands (broken lines) and 95 % confidence bands of timereg (regions). This figure appears in color in the electronic version of this article

We finally illustrate how the joint asymptotic distribution of the baseline and parameter estimators can be used with other functionals. To this end, consider the restricted mean life of an individual with lifetime and covariate vector , that is, with estimator . We refer to Karrison (1987) for the concept of adjusting the restricted mean life to covariates. To get a description of the uncertainty of based on multiplier bootstrap constructions we can simply apply the functional to the obtained bootstrap samples. It follows that has the same asymptotic distribution as due to Hadamard differentiability. Thus, we can easily construct symmetric 95% confidence intervals for the restricted mean and their differences based on the bootstrap. The key point being that these are very easy to get at when the bootstrap estimates are at hand. In addition, still due to Hadamard differentiability, differences in restricted mean life and uncertainties thereof are also estimable using the obtained bootstrap estimates.

For example, the direct multiplier bootstrap approach based on and standard normal multipliers yields that males with diabetes have a restricted mean life within the first 5 years at 3.87 (3.74, 4.00) for males without diabetes and 3.15 (2.91, 3.41) with diabetes. Males with diabetes thus lose 0.71 (0.49, 0.93) years within the first 5 years. Confidence intervals for other functionals are similarly found by means of the functional Δ method.

Cumulative incidence for competing risks

We here consider data on bone marrow transplantation from patients that are treated for myelodysplasia (Sierra et al., 2002). This data is also available in the the timereg R‐package. Of 408 patients 161 died from treatment‐related causes, and 87 had relapse before dying, and the remaining were censored. We consider these as two competing risks of death.

We wish to estimate the risk of dying from treatment‐related causes by a certain point in time, and do this by modeling the two cause‐specific hazards using Cox regression models based on the covariates platelet count (High 1/Low 0), tcell‐depletion (yes 1/no 0), and standardized age ([age‐35]/15). We then combine the estimates for a particular covariate combination (Low platelets, no tcell‐depletion, and average age (35 years) to estimate the cumulative incidence , where is the probability of surviving both causes, and is the intensity for a treatment‐related death.

We denote the baseline cumulative hazards and the parameters of both Cox models by . Note that the Cox models for each cause‐specific hazard are resampled independently. Their combination yields inference on the cumulative incidence functional

Again, has the same asymptotic distribution as due to Hadamard differentiability of the functional. We thus simply need to apply the multiplier bootstrap to the two Cox models and then combine estimates. Subsequently, we can construct a confidence band for based on the realizations at hand. Figure 3 shows the cumulative incidence curve for a subject with covariate vector with 95% equal precision confidence bands with (dotted lines) and without log‐transform (broken lines).

Figure 3

Cumulative incidence curve for a subject with covariate with 95% equal precision confidence bands with (dotted lines) and without log‐transform (broken lines)

DISCUSSION AND FURTHER RESEARCH

Despite their importance and considerable interest in making survival predictions based on semiparametric regression models such as the Cox model, confidence bands are not used much in practice. This is probably because the most commonly‐used Cox model routines do not implement simultaneous confidence bands. The aim of this work is to investigate some natural and simple multiplier bootstrap approaches for filling this gap. In particular, we show in the supporting information that the proposed bootstrap solutions do asymptotically have the desired properties. A key point in our proofs is that major characteristics of the procedure can be deduced from martingale arguments. This enormously facilitates the transfer of classical proofs for the estimators to their multiplier bootstrap counterparts. The approach evidently generalizes to more complex models as long as martingale structures for the relevant counting processes are preserved, that is, Cox models in multistate models or Fine‐Gray regression models for subdistribution functions (Fine & Gray, 1999) with administrative censoring. Future work will focus on the adaptation to more complex designs.

In addition, we consider the finite sample performance of various confidence bands. It turns out that the choice of transformation and weight function considerably influence widths and coverage of the bands. Among the bootstrap weights, centered exponential multipliers exhibit the best coverage probability. In particular, their combination with the log‐transformation leads to very accurate results, even for a small sample size (). Furthermore, the equal precision bands exhibit much smaller widths in comparison to the Hall‐Wellner bands.

Another nice feature of our approach is its simple extensibility to find confidence intervals for functionals. We illustrated this with the restricted mean life based on a Cox model.