Survival of the Fittest Group: Factorial Analyses of Treatment Effects under Independent Right-Censoring
Dennis Dobler Markus Pauly




1*]Dennis Dobler
2]Markus Pauly


[1]Department of Mathematics, Vrije Universiteit Amsterdam, Netherlands.
[2]Institute of Statistics, Ulm University, Germany.

Survival of the Fittest Group: 
 Factorial Analyses of Treatment Effects 
 under Independent Right-Censoring
    [
    December 30, 2023
============================================================================================================




^* e-mail:  d.dobler@vu.nl

This paper introduces new effect parameters for factorial survival designs with possibly right-censored time-to-event data. In the special case of a two-sample design it coincides with the concordance or Wilcoxon parameter in survival analysis. More generally, the new parameters describe treatment or interaction effects and we develop estimates and tests to infer their presence. We rigorously study the asymptotic properties by means of empirical process techniques and additionally suggest wild bootstrapping for a consistent and distribution-free application of the inference procedures. 

The small sample performance is discussed based on simulation results.
The practical usefulness of the developed methodology is exemplified on a data example about patients with colon cancer by conducting one- and two-factorial analyses.






Keywords: Empirical process,  Factorial designs, Kaplan-Meier estimator, Quadratic forms, Wild bootstrap





 








§ MOTIVATION AND INTRODUCTION

Many biomedical and clinical trials are planned as factorial designs. Here, not only the (main) effects of separate factors but also interaction effects that are related to possibly complex factor combinations are of importance. 
Such interaction effects

may even alter the interpretation of main effects leading to the established comment by <cit.> that 
`it is desirable for reports of factorial trials to include estimates of the interaction between the treatments'.
 
On the other hand, nonparametric estimation and the inference of adequate effects in such designs can be rather involved. 
In particular, most existing inference procedures have focused on testing hypotheses formulated in terms of distribution functions <cit.> which cannot be inverted to obtain confidence intervals or regions for meaningful effects. 
Only recently, nonparametric methods for inferring adequate effects in general factorial designs with independent and dependent observations have been established <cit.>. These procedures are, however, only developed for completely observed data and not applicable for 

partially observed time-to-event data. 

Since many clinical studies are concerned with survival outcomes, adequate statistical inference methods for complex factorial time-to-event designs are of particular interest. 



To detect main effects, weighted logrank tests or their extensions may be applied in case of two or multiple samples  
<cit.>. However, these procedures only infer conclusions

in terms of cumulative hazard functions and cannot be applied to obtain 
concrete effect parameters with informative confidence intervals nor tests for the presence of interactions. 
In practice, interaction effects are usually modeled with the help of
Cox-, Aalen- or even Cox-Aalen regression models <cit.> with factors as covariates and incorporated interaction terms. However, although very flexible, these models are usually more driven towards hazards modeling by continuous covariates while the incorporation of several factor variables (e.g., via multiple dummy variables per factor) can become cumbersome; especially when interactions are incorporated, see also <cit.> and <cit.> for the uncensored case. 

The above problems directly motivate a nonparametric approach for estimating and inferring main and interaction effects in factorial designs with censored observations. 
So far, the only existing methods in this context are given by the nonparametric survival procedures of <cit.> and <cit.>.
They are based on a purely nonparametric model that does not require any multiplicative or additive structure of the hazards and can even be applied for arbitrary, possibly non-continuous survival distributions (i.e., it can be readily used for survival times rounded to days, weeks or months). 
Moreover, it leads to tests for

main and interaction effects in case of independent right-censored data. 
However, these tests suffer from several drawbacks: 
the procedure is based on a rather strong assumption on the underlying censorship distribution 
which is often hard to verify in practical situations. In addition, 
null hypotheses are only formulated in terms of distribution functions. 
As a result, there is no direct quantification and estimation of main and interaction effects in terms of confidence intervals as, e.g., required by regulatory authorities (ICH E9 Guideline, 1998, p. 25).

This is to be changed in the current paper. We develop and rigorously analyze nonparametric inference procedures, i.e. tests and confidence intervals, for meaningful effect sizes in factorial survival designs, where data may be subject to random right-censoring. 

Similar to the adaption of the <cit.> test to the two-sample survival set-up by <cit.>, we consider the recently proposed unweighted nonparametric effects of <cit.> and extend their ansatz to a general survival setting. In the special case of proportional hazards, these effects have a direct relationship to hazard ratios in two-sample settings <cit.> while they remain meaningful in case of non-proportional hazards. This fact makes the effect sizes even more appealing for practical purposes.

The paper is organized as follows. The statistical model and important results on the basic estimators are presented in Section <ref>. The resulting test statistic for the null hypotheses of interest is stated and mathematically analyzed in Section <ref>. Since the asymptotic distribution of the test statistic depends on unknown parameters, 
we propose a distribution-free multiplier resampling approach in Section <ref> and prove its consistency. 
In Section <ref> it is supplemented by a simulation study to assess the finite sample properties of the proposed procedure.
They are then exemplified on a colon cancer study in Section <ref>, where in the original study <cit.> the analysis was made in terms of Cox models.
Finally, the paper closes with concluding comments in Section <ref>.
All proofs are deferred to the technical Appendix.




§ THE SET-UP



To establish the general model, we consider sequences of mutually independent random variables

    T_ikind∼ S_i      C_ikind∼ G_i     ( i = 1, …, d,  k = 1, …, n_i),

where T_ik denotes the actual survival time of subject k in group i and C_ik the corresponding censoring variable. Moreover, to even allow for ties or survival times rounded to weeks or months, the survival functions S_i and G_i, i=1, …, d, defined on (0,∞) may be possibly discontinuous.
That is, the corresponding hazard rates may, but need not exist. 
The actually observable data consist of the right-censored survival times X_ik = T_ik∧ C_ik and the uncensoring indicators δ_ik = 1 { T_ik≤ C_ik}, 
 i = 1, …, d,  k = 1, …, n_i. In this set-up, a factorial structure can be incorporated by splitting up indices, see Section <ref> for details. 
 
 
 In the special case of d=2 groups with continuous survival times <cit.> introduced an estimator for the concordance probability
 
    w= P(T_11 > T_21) = - ∫  S_1 dS_2

 that a randomly chosen subject from the first group survives longer than someone from the second group. 
 If all subjects are completely observable, this effect size w reduces to the well-known Mann-Whitney effect underlying the <cit.> test.
 Inference procedures for w and related quantities in survival set-ups (such as the concordance parameter or the average hazard ratio)  have, e.g. been developed by <cit.>. However, an extension of the definition of w to the more general design (<ref>), allowing for an arbitrary factorial structure, is not straightforward. In particular, for the case of completely observed data, <cit.> and <cit.> point out several pitfalls that may lead to paradoxical results when working with a `wrong' extension of w. Adopting their solution to the present situation, we  introduce an additional `benchmark' survival time Z, independent of the above, with averaged survival function Z∼S̅ = 1/d∑_i=1^d S_i. This is used to extend w to 
 
 
    p̃_i = P(T_i1 > Z) + 1/2 P(T_i1 = Z) = 
     - ∫ S_i^±̣̅S,

 where the superscript ± denotes the average of a right-continuous function and its left-continuous version. 
 The use of such normalized survival functions adequately handles discrete components of the survival distribution, i.e. ties in the data are explicitly allowed.
  
 The choice of the effect parameter p̃_i is motivated by recent findings on nonparametric analyses of factorial designs with 
 complete observations in <cit.>. dobler2017nonparametric2 They stress that other choices, e.g.
 pairwise-comparisons of all concordance probabilities w or comparisons with the weighted survival function
∑_i=1^d n_i/NS_i instead of S̅, may easily result in paradoxical outcomes. This is no issue for the effects p̃_i which are sample size independent. 
 
 For later calculations, we emphasize that the effect parameters are balanced in the mean.
 In particular, we have 
 
 
 
    1/d∑_i=1^d p̃_i = - ∫S̅^±̣̅S = 1/2  and  p̃_i = - ∑_j = 1 
     j ≠ i^d ∫ S_i^±Ṣ_j + 1/2d.

 From a practical point of view, estimation of the p̃_i's would need 'arbitrarily' large survival times since the integral is defined on (0,∞).
 However,  every study ends at a certain point in time. For practical applicability, we therefore assume that the censoring times are bounded and we have to modify the p̃_i's accordingly:  
 denote by τ > 0 the largest possible censoring time. 
 In comparisons of survival times, which belong to different groups and which exceed τ,
 
 no group shall be favored.  In other words, the remaining mass has to be split up equally among the groups.  Technically, this is realized by setting the remaining mass of the survival functions to zero: S_i(τ) =0. 
 Redefining S_i and S̅ from now on as the survival functions of min(T_i1, τ) and min(Z, τ), respectively, this translates into the nonparametric concordance effects 
 
    p_i = P(min(T_i1, τ) > min(Z, τ)) + 1/2 P(min(T_i1, τ) = min(Z, τ)) = - ∫ S_i^±̣̅S.

 Obviously, all of the above-discussed positive properties of the effects parameter p̃_i also transfer to the nonparametric concordance effects p_i:
 it is a meaningful effect measure for ordinal and metric data, sample size independent, and allows for a suitable treatment of ties. 
 
  We aggregate all effects into the vector p̱ = (p_1, …, p_d)' and borrow a trick from <cit.> and <cit.> to express them as
  
    p̱  = (I̱_d ⊗1/d1̱_d') · (w̱_1', …,  w̱_d')' =: E̱_d ·w̱.

Here, w̱_i = (w_1i, …, w_di)' = - ∫ S_i^±̱̣S is the ^d-vector of effects for direct comparisons of group i with respect to all groups j=1,…,d,
 and S̱ = (S_1, …, S_d)' is the aggregation of all survival functions. Moreover, I̱_d denotes the identity matrix in ^d, 1̱_d  the d-dimensional vector of 1's and the symbol ⊗ denotes the Kronecker product. In this way the ith entry of w̱_i is w_ii = 1/2 which makes sense because equal groups should be valued equally high. 
 Anyhow, Equation (<ref>) shows that the problem of estimating p̱ reduces to the estimation of the pair-wise effects w_ji. But this can be achieved by substituting each involved survival function S_i by its Kaplan-Meier estimator S_i, i=1,…,d <cit.>. Proceeding in this way we denote by 
 w̱ and w̱_i these estimated counterparts of w̱ and w̱_i. Let N = ∑_i=1^d n_i be the total sample size. 
 Below we establish the asymptotic normality of 
 √(N)(w̱ - w̱) under the following framework:
 
    N^-1 n :=  (n_1/N, …, n_d/N)' →λ := (λ_1, …, λ_d)' ∈ (0,1)^d

 as min n →∞. To give a detailed description of the resulting asymptotic covariance structure, however, we first have to introduce some additional notation: 
  Let D[0,τ] be the space of all càdlàg-functions on [0,τ], equipped with the Skorokhod metric, and BV[0,τ] ⊂ D[0,τ] its subspace of càdlàg-functions with bounded variation.
 For the subsequent arguments it is essential that we can represent  w̱ = ϕ∘S̱  as a functional of S̱. In particular, the functional
 
    ϕ:   (BV[0,τ])^d →^d^2,      (f_1, …, f_d)' ⟼( - ∫ f_i^±f̣_j )_i,j = 1^d,

 with inner index j, is Hadamard-differentiable at S̱; see the proof of Lemma <ref> below for details.
 We denote its Hadamard-derivative at S̱ by ϕ̣_S̱, which is a continuous linear functional.
 
 
 
 
 
 
 
 
 For technical reasons, we assume throughout that P(T_i1 > τ) >0 for all groups i=1, …, d.
 We may now state the first preliminary but essential convergence result.
 
 Under the asymptotic regime (<ref>) we have
 
    √(N) ( w̱ - w̱ ) W̱,

 where W̱ has a centered multivariate normal distribution on ^d^2.
 
 In particular, we can write W̱ = ϕ̣_S̱· diag(λ)^-1/2U̱,
 where U̱ consists of independent, zero-mean Gaussian processes U_1, …, U_d with covariance functions 
 
    Γ_i(r,s) = S_i(r) S_i(s) ∫_0^r ∧ sΛ̣_i/S_i- G_i- (1 - ΔΛ_i),    i=1,…,d,

 where Λ_i denotes the cumulative hazard function corresponding to S_i =  (1 - Λ̣_i), i=1, …, d; 
 the symbol  denotes the product integral <cit.>. 
 Here, a minus sign in a subscript indicates the left-continuous version of a function and ΔΛ = Λ - Λ_- is the jump size function of Λ.
Note that the covariance matrix of W̱ is singular; in particular, ( ϕ̣_S̱·  diag(λ)^-1/2U̱ )_i,i = 0 for all i = 1, …, d.
 The other entries (i ≠ j) are distributed as follows:
 
    ( ϕ̣_S̱·  diag(λ)^-1/2U̱ )_i,j 
      =  1/√(λ_i)∫ U_i^±Ṣ_j 
        - 1/√(λ_j)∫ U_j^±Ṣ_i  
      ∼ N (0, 1/λ_i∫∫Γ_i^±±Ṣ_j Ṣ_j 
      + 1/λ_j∫∫Γ_j^±±Ṣ_i Ṣ_i ) .

 Here, the double appearance of ± signs means the average of all four combinations of left- and right- continuous versions in both arguments of a two-parameter function.
 
 
 Let us now turn to the estimation of the nonparametric concordance effects p̱.
 A matrix multiplication of w with E̱_d from the left is basically the same as taking the mean with respect to the inner index j. This immediately brings us to the first main result:
 
  
  Under the asymptotic regime (<ref>) we have
 
    √(N) ( p̱ - p̱) := √(N)E̱_d ( w̱ - w̱ ) E̱_d W̱ =( 1/d∑_i=1^d 1/√(λ_i)∫ U_i^±Ṣ_j - 1/√(λ_j)∫ U_j^±̣̅S  )_j=1^d,
 
 where E̱_d W̱ has the variance-covariance matrix V̱ with the following entries:
 
    V_ii = 1/λ_i∫∫Γ_i^±±̣̅S ( S̅ - 2/d S_i )
      + 1/d^2∑_j=1^d 1/λ_j∫∫Γ_j^±±Ṣ_i Ṣ_i

  in the ith diagonal entry, i = 1, …, d, and 
  
    V_ij = 1/d^2∑_j=1^d 1/λ_j∫∫Γ_j^±±Ṣ_i Ṣ_j
       - 1/d1/λ_i∫∫Γ_i^±±̣̅S Ṣ_j 
       - 1/d1/λ_j∫∫Γ_j^±±̣̅S Ṣ_i

   in the off-diagonal entries (i,j),  i ≠ j.
 
 A more compact form of the matrix V̱ is given in Appendix <ref>.
 
 
 
  
  

§ CHOICE OF TEST STATISTIC

  
  
  In order to develop hypothesis tests based on the estimator p̱, 
  we next need to find a consistent estimator V̱_N for V̱. 
 
  A natural choice  is to plug in estimators for all unknown quantities that are involved in V̱.
  In particular, we use the Kaplan-Meier estimators for all survival functions and
  Γ_i(s,t) =  S_i(s)  S_i(t) n_i ∫_0^s ∧ t [ Y_i (1 - ΔΛ_i)]^-1Λ_i for each covariance function Γ_i,
  where Y_i is the number at risk process and Λ_i is the Nelson-Aalen estimator for the cumulative hazard matrix in group i.
  Note that if ΔΛ_i(u) = 1, we also have S_i(u) = 0 in which case we let Γ_i(s,t) = 0 if s ≥ u or t ≥ u.
  We denote the resulting covariance matrix estimator by V̱_N.
  
  which results in
  
    V̱_N 
        
        
        
        
          =     N [ 1/d^2∫∫1̱_d' diag( n)^-1Γ^±±S̱S̱'
        + diag( n)^-1∫∫ diag(Γ^±±) S̅S̅
        - 1/d∫∫ diag( n)^-1Γ^±±S̅S̱'
        - ( 1/d∫∫ diag( n)^-1Γ^±±S̅S̱' )' ].

  
   
  Under the asymptotic regime (<ref>) we have the consistency
  V̱_N V̱. 
  

  
  All of the developed convergence results are now utilized to find the most natural test statistic.
  First note that the asymptotic covariance matrix V̱ is singular 
  since 1̱_d' √(N) (p̱ - p̱) ≡ 0, 
  whence r(V̱) ≤ d-1 follows.
  Furthermore, it is not at all obvious whether
  the ranks of the Moore-Penrose inverse r( (C̱V̱_N C̱' )^+ ) 
  converge in probability to the rank r( (C̱V̱C̱')^+)
  for a compatible contrast matrix C̱.
  Hence, the Wald-type statistic 
  N p̱' C̱' (C̱V̱_N C̱')^+ C̱p̱
  is not suitable for testing H_0^p(C̱): C̱p̱ = 0̱:
  Its asymptotic behaviour is unclear and, hence, there is no reasonable choice of critical values.
  
  
  Instead, we utilize a statistic that does not rely on the uncertain convergence of ranks of generalized inverses.
  This leads us to the survival version of the so-called ANOVA-rank-type statistic
  
    F_N(Ṯ) = N/tr(ṮV̱_N)p̱' Ṯp̱,

  where Ṯ = C̱' (C̱C̱')^+ C̱ is the unique projection matrix onto the column space of C̱. 
  
  Below we analyze both, its asymptotic behaviour under null hypotheses of the from H_0^p(C̱): C̱p̱ = 0̱ and under the corresponding alternative hypotheses H_a^p(C̱): C̱p̱≠0̱.
  
   
   Assume the asymptotic regime (<ref>) and that tr(ṮV̱)>0.
   
    
  a) Under  H_0^p(C̱) and as N →∞, we have
   F_N(Ṯ) χ = W̱' E̱_d' ṮE̱_d W̱/tr(ṮV̱)
   which is non-degenerate and non-negative with E(χ)=1.
    
  b) Under H_a^p(C̱) and as N →∞, we have
    F_N(Ṯ) ∞.
   
  
  
  
  
  
  
  Rather not:
  
  Another, commonly used approach for obtaining critical values in a multi-sample problem
  is the permutation technique applied to the ANOVA-type statistic
  
    F_N(C̱) = F_N(Ṯ) = N/tr(ṮV̱_N)p̱' Ṯp̱,

  where Ṯ = C̱' (C̱C̱')^+ C̱ is the unique projection matrix on the column space of C̱.
 
  As the distribution of χ depends on unknown quantities (cf. Theorem <ref>) the test statistic F_N(Ṯ) in (<ref>) is no asymptotic pivot. To nevertheless 
  obtain proper critical values which lead to asymptotically exact inference procedures we next propose and study a resampling approach.


  
  
  

§ RATHER NOT: PERMUTATION APPROACH

  
  Simulations (by Frank) have shown that the following approach does not yield a good approximation:
  
  
  
  A possible problem when trying to apply the permutation technique is that the eigenvalues involved in the asymptotic distribution of F_N(Ṯ) may be altered in the asymptotic distribution of the permutation version of F_N(Ṯ).
  Nevertheless, this permutation technique may be utilized for a new kind of approximation technique
  that will be competing with the ANOVA-Box-type-p tests of Brunner et al. (2016)
  or with related ANOVA-Pearson-type-p tests.
  The new approximation works in the following way:
  First we introduce a studentized test statistic, i.e.,
  
    F_N,stud(Ṯ) = N p̱' Ṯp̱ - tr(ṮV̱_N)/√(2 tr(ṮV̱_N ṮV̱_N)).

  This choice of studentization is justified by the asymptotic expectation of N p̱' Ṯp̱ under H_0^p being the sum of eigenvalues ∑_i=1^d ϱ_i(ṮV̱)
  which in turn is the same as the trace.
  The asymptotic variance of N p̱' Ṯp̱ under H_0^p equals 2 ∑_i=1^d ϱ_i^2(ṮV̱).
  But the sum of the squared eigenvalues is the same as the sum of eigenvalues of (ṮV̱)^2, i.e., ∑_i=1^d ϱ_i(ṮV̱ṮV̱).
  
  
  
  
  We now establish a similar studentization of the permutation version of the ANOVA-type statistic.
  To this end, denote by p̱^π and V̱_N^π the permutation versions 
  of p̱ and V̱, respectively.
  By the same arguments as above, similar estimators for the asymptotic expectation and variance are utilized, leading to
  
    F_N,stud^π(Ṯ) = N (p̱^π)' Ṯp̱^π - tr(ṮV̱_N^π)/√(2 tr(ṮV̱_N^πṮV̱_N^π)).


  It remains so be proven that the conditional distribution of N (p̱^π)' Ṯp̱^π is consistent for L(∑_i=1^d ϱ_i(ṮV̱^π) C_i^2) 
  where the C_i are i.i.d. N(0,1) random variables and V̱^π is the limit in probability of V̱_N^π.
  Furthermore, simulations will be conducted in order to compare the performance of the suggested permutation-based approximation procedure to the well-known Box-type approximation technique.
  
  
  
  
  

§ KOMMENTARE VON MARKUS

  
  Markus meinte, dass Frank einen aehnlichen Permutationsansatz wie oben bereits ausprobiert hatte... die Ergebnisse waren leider nicht gut.
  
  Generell sollten wir entweder wieder eine Momenten-Anpassungsmethode waehlen...
  oder direkt den Wild BS! Angewandt auf die Martingaldarstellung des Kaplan-Meiers!
  Dann sollte der richtige Limes herauskommen und alles muesste klappen!
  
  Wir wuerden wieder ab irgendeinem Punkt abschneiden und dahinterliegende Masse auf alle Gruppen verteilen (wie in unserem ersten Konkordanz-Papier).
  In einem Kommentar koennen wir festhalten, dass wir auch jeweils bis zu den Endzeitpunkten gehen koennten.
  Der Kaplan-Meier ist nach Gill auch in dieser Situation noch schwach konvergent.
  Der BS laesst sich nach meiner Diss auch darauf anwenden!
  Und der Wild BS ebenfalls, wenn man wieder mit Martingalargumenten an das ganze herangeht!!!

 
  
  
  

§ MOMENT-MATCHING APPROXIMATIONS

  
  streichen!?
  
  In other situations it was seen that approximation techniques aiming at matching moments of quadratic form statistics may yield similar or even better small sample results when compared to asymptotically exact inference procedures based on different kinds of resampling techniques;
  see ... and ... for such approximations in the context of approximating Cramér-von Mises-type statistics in Aalen-Johansen estimators for cumulative incidence functions under independent right-censoring and left-truncation.
  

  In order to be able to match as many moments as possible, we introduce a studentized version of the above ANOVA-type statistic, i.e.
  
    F_N,stud(Ṯ) = N p̱' Ṯp̱ - tr(ṮV̱_N)/√(2 tr(ṮV̱_N ṮV̱_N)).

  This choice of studentization is justified by the asymptotic expectation of N p̱' Ṯp̱ under H_0^p being the sum of eigenvalues ∑_i=1^d ϱ_i(ṮV̱)
  which in turn is the same as the trace.
  The asymptotic variance of N p̱' Ṯp̱ under H_0^p equals 2 ∑_i=1^d ϱ_i^2(ṮV̱).
  But the sum of the squared eigenvalues is the same as the sum of eigenvalues of (ṮV̱)^2, i.e. ∑_i=1^d ϱ_i(ṮV̱ṮV̱).
  Using a principle component decomposition, it is furthermore easy to see that the asymptotic null distribution of F_N(Ṯ) is L(∑_i=1^d ϱ_i(ṮV̱) B_i^2) 
  where the B_i are i.i.d. N(0,1) random variables.
  We would like to stress that the estimators involved in the studentization (<ref>) are obviously consistent under condition (<ref>).
   

§ INFERENCE VIA MULTIPLIER BOOTSTRAP

  
  
  In this section, we apply suitably tailored multiplier bootstrap techniques
  in order to approximate the small sample distribution of F_N( Ṯ ).
  To this end, we consider the situation under H_0^p(C̱) in which case we may expand
  
    F_N( Ṯ )     = N/tr(ṮV̱_N) ( p̱ - p̱ )' Ṯ ( p̱ - p̱ ) 
        
         = N/tr(ṮV̱_N) ( ϕ̣_S̱· ( S̱ -  S̱  ))' E̱_d'  Ṯ E̱_d ( ϕ̣_S̱· ( S̱ - S̱ )) + o_p(1) ,

  where S̱ is the vectorial aggregation of all Kaplan-Meier estimators S_1, …,  S_d.
  First, we replace the martingale residuals, that are attached to the Kaplan-Meier estimators, 
  with independent centered random variables which have approximately the same variance.
  In particular, we replace √(N) (S_i - S_i) with 
  
  
  
    S(t) ·√(N)∑_k=1^n_i G_ik∫_0^t [ (Y_i(u) - Δ N_i (u)) Y_i(u) ]^-1/2Ṇ_ik(u);

  cf. <cit.> for a similar wild bootstrap Greenwood-type correction for tied survival data.
  Here we utilized the usual counting process notation <cit.>: 
  N_ik indicates whether the event of interest already took place for individual k in group i.
  The wild bootstrap multipliers G_ik, i=1, …, n_i,  i=1, …, d, 
  are i.i.d. with zero mean and unit variance and also independent of the data.
  In <cit.> a similar multiplier resampling approach is
  applied to Nelson-Aalen and Aalen-Johansen estimators in one- and two-sample problems.
  
 In a next step toward the construction of a wild bootstrap statistic, we replace ϕ̣_S̱ with ϕ̣_S̱.
 Let us denote the thus obtained wild bootstrap version of 
  √(N)ϕ̣_S̱· ( S̱ - S̱ ) by
  W̱_N^*. Conditionally on the data, this d^2-variate random vector is for large N 
  approximately normally distributed and its limit distribution coincides with that of W̱; see the proof of Theorem <ref> below for details.
  
  
   In the case of continuously distributed survival data one might as well extend both central limit theorems
   to the unrestricted Skorohod space, i.e. to D[0,τ_j] instead of D[0,K].
   Here τ_j is the right end point of the support of S_j which may be different for each j=1, …, d.
   This convergence in distribution of the univariate Kaplan-Meier estimator was shown in <cit.> and <cit.>.
   That this convergence result holds as well for resampling counterparts was first shown in <cit.>,
   where Efron's (1981)efron81 bootstrap for the Kaplan-Meier estimator was extended to the whole support of the survival function with the help of martingale arguments.
   The same martingale arguments apply 
   when considering the martingale properties (cf. ) of the wild bootstrapped Nelson-Aalen estimators
   
    √(N)∑_i=1^n_j G_ij∫_0^t √(Y_j(u) - Δ N_j(u)/Y_j^3(u))Ṇ_ij(u)

    which contribute essentially to the Kaplan-Meier estimators. 
    
    
  
  Finally, a wild bootstrap version F_N^*(Ṯ) of F_N(Ṯ) requires 
 that we also use a consistent wild bootstrap-type estimator tr(ṮV̱_N^*) of tr(ṮV̱_N).
  It is found by replacing the estimators Γ_i with
  
    Γ^*_i(s,t) =  S_i(s)  S_i(t)  n_i ∑_k=1^n_i G_ik^2 ∫_0^s ∧ tṆ_ik/ (Y_i - Δ N_i) Y_i.

  Its conditional consistency was argued in <cit.> and a sufficient condition for this is E(G_11^4) < ∞.
  These wild bootstrap-type variance estimators also have the nice interpretation of optional variation processes of the wild bootstrapped Kaplan-Meier estimators (Dobler, 2017). 
  
  Hence, the resulting wild bootstrap version of F_N( C̱ ) is
  
    F_N^*( Ṯ )     = 1/tr(ṮV̱_N^*)W̱_N^*' E̱_d'  Ṯ E̱_d W̱_N^*.

  The following conditional central limit theorem ensures the consistency of this resampling approach. 
  
   
   Assume E(G_11^4) < ∞ and that the conditions of Theorem <ref> hold. 
   Conditionally on (X_ik, δ_ik), i =1, …, d,  k = 1, …, n_i,  we have for all underlying values of p̱
   
    F_N(Ṯ)^* χ = W̱' E̱_d' ṮE̱_d W̱/ tr(ṮV̱)

   in probability as N →∞.
  
 We would like to stress that the limit distribution coincides with that of F_N(Ṯ) under H_0^p(C̱). 
 For the wild bootstrap version F_N(Ṯ)^*, however, the convergence result holds under both, the null and the alternative hypothesis, 
 i.e. its conditional distribution always approximates the correct null distribution of the test statistic. 
  
  
   Another possible resampling technique for F_N(Ṯ) is the bootstrap by <cit.>,
   where, for calculating resampling versions of the Kaplan-Meier estimators,  
   d independent bootstrap samples with sample sizes n_1, …, n_d are drawn separately for each group.
   On the other hand, bootstrapping based on pooled bootstrap samples 
   would result in modified covariances
   and thus in an incorrect limit distribution of the quadratic form.
   This is why the pooled resampling approach of <cit.> does not generalize to the present set-up.
  


  We conclude the theoretical part of this article with a presentation of deduced inference procedures for the effect sizes p̱. 
  
  To this end, let c^*_N,α denote the (1-α)-quantile, α∈ (0,1), of the conditional distribution of F_N^*(Ṯ) given the data.
  In practice, this quantile is approximated via simulation by repeatedly generating sets of the wild bootstrap multipliers G_ki.
  
   
   Under the assumptions of Theorem <ref>, the test
   
    φ_N = 1{ F_N(Ṯ) > c_N,α^* }

   is asymptotically exact and consistent.
   That is, E( φ_N ) →α· 1_H_0^p(C̱) + 1_H_a^p(C̱) as N →∞.
  



Now, let r be the number of columns of C̱' and denote by c̱_1, …, c̱_r its column vectors.
The presentation of a simultaneous confidence region for the contrasts c̱_ℓ' p̱, ℓ = 1, …, r, in Corollary <ref> below will be done in an implicit manner. 

 
   
   Under the assumptions of Theorem <ref>, an asymptotically exact (1-α)-confidence ellipsoid for the contrasts c̱_ℓ' p̱, ℓ = 1, …, r,   
   is given by 
   
    CE = CE_N,1-α(C̱) = {v̱∈^r :
        (C̱p̱ - v̱)' (C̱C̱')^+ ( C̱p̱ - v̱) ≤tr(ṮV̱_N)/N c_N,α^* }.

   That is, P( C̱p̱∈ CE ) → 1 - α as N →∞.
   
  
  
  
  

§ SIMULATIONS


  In this section, we assess the small sample properties of the test φ_N as proposed in Corollary <ref>. 
  
  

 §.§ Behaviour under null hypotheses

  
    We first focus on its type I error control with respect to 
  
   
  * various kinds of contrast matrices
   
  * and different censoring intensities.
  
Design and Sample Sizes. For ease of presentation we restrict ourselves to a design with d=6 groups with different sample size layouts: 

  we considered small samples in a balanced design with ṉ_1 = (n_1,…,n_6)' = (10,10,10,10,10,10)' and two unbalanced designs with ṉ_2 = (n_1,…,n_6)' = (10,12,14,10,12,14)' and ṉ_3 = (10,12,14,14,10,12)', respectively. 
  To obtain designs with moderate to large sample sizes we increase these vectors component-wise by the factors K ∈{2,3,5,10}. Moreover, depending on the question of interest, we below distinguish between a one-way layout with six independent groups and a 2× 3 two-way design.

  Censoring Framework. We considered exponentially distributed censoring random variables C_i1ind∼Exp(λ_i) with the following vectors 
  λ = (λ_1,…,λ_6)' of rate parameters: 
   
    λ_1 = 0.4 · 1,  λ_2 = 0.5 · 1,  λ_3 = 2/3 · 1,   λ_4 = (0.4, 0.5, 2/3, 0.4, 0.5, 2/3)',  λ_5 = (0.4, 0.5, 2/3, 2/3, 0.5, 0.4)',
  where 1 ∈^6 is the vector consisting of 1s only. 
  Thus, the first three settings correspond to equal censoring mechanisms with increased censoring rate from λ_1 to λ_3. The other two (λ_4 and 
  λ_5) lead to unequal censoring. By considering all 75 possible combinations, many possible effects of censoring and sample size assignments are analyzed. 
  For example, in the set-up with ṉ_2, K=10 and λ_4, larger sample sizes are matched with a stronger censoring rate in an unbalanced design.

  Contrast Matrices and Null Hypotheses.
  We simulated the true significance level of the tests for the null hypotheses H_0^p(C̱): C̱p̱ = 0̱ for two designs and different contrast matrices of interest:

  In case of a one-way design with d=6 groups we were interested in the null hypotheses of `no group effect'  or        
  `equality of all treatment effects' H_0^p(C̱_1): {C̱_1 p̱ = 0̱} = {p_1=…=p_6}. This may be described by considering the matrix C̱_1 = P̱_6, where here and below P̱_d = I̱_d - J̱_d/d ≡I̱_d -  1_d  1_d' /d denotes the d-dimensional centering matrix.


    C̱_1 = [  1 -1         ;     1 -1      ;        ⋱  ⋱   ;           1 -1 ]∈^5 × 6.

 Next, we consider a 2× 3 two-way layout with two factors A (with two levels) and B 
  (with three levels). This is incorporated in Model (<ref>) by setting d=a· b= 2· 3 =6 and splitting up the index i into two indices i_1=1,2 (for the levels of factor A) and i_2=1,2,3 (for the levels of factor B). Thus, we obtain survival times T_i_1i_2k, k=1,…,n_i_1i_2, and corresponding 
 nonparametric concordance effects p_i_1i_2. More complex factorial designs can be incorporated similarly. In this 2× 3 set-up we are now interested in testing the null hypotheses of 
  
   
  (A) `No main effect of factor A': H_0^p(C̱_2,A): {C̱_2,Ap̱ = 0̱} = {p̅_1· = p̅_2·},
   
  (B) `No main effect of factor B': H_0^p(C̱_2,B): {C̱_2,Bp̱ = 0̱} = {p̅_· 1 = p̅_· 2 = p̅_· 3} and
   
  (AB) `No A× B interaction effect': H_0^p(C̱_2,AB): {C̱_2,ABp̱ = 0̱} 
   = {p_i_1i_2 - p̅_i_1· -p̅_· i_2 + p̅_·· = 0  for all  i_1,i_2},
  
  where p̅_i_1·, p̅_· i_2 and p̅_·· denote the means over the dotted indices. In particular, the corresponding contrast matrices are given by 
   C̱_2,A = P̱_2 ⊗1/3J̱_3,
 C̱_2,B = 1/2J̱_2  ⊗P̱_3, and C̱_2,AB = P̱_2 ⊗P̱_3,
where ⊗ indicates the Kronecker product.

Survival Distributions. For ease of presentation we only considered a rather challenging scenario, where the groups follow different survival distributions. In particular, we simulated 



 
  (G1) a lognormal distribution with meanlog parameter 0 and sdlog parameter 0.2726 for the 
 first group,
 
  (G2)  a Weibull distribution with scale parameter 1.412 and shape parameter 1.1 for the second group,
 
  (G3) a Gamma-distribution with scale parameter 0.4 and shape parameter 2.851 for the third group and
 
  (G4-G6)        mixing distributions of all pair combinations of the first three survival functions for the last three groups.

The first three survival functions are illustrated in Figure <ref>. We note that preliminary simulations for more crude scenarios with identical survival distributions in all groups exhibited a much better type-I-error control of our testing procedure (results not shown).




Anyhow, the parameters of the above distributions were chosen in such a way that the nonparametric concordance effects of all groups are equal, i.e. p_i = 0.5 for all i=1,…, 6 (one-way) and 
p_i_1i_2 = 0.5 for all i_1=1,2;  i_2=1,2,3 (two-way), respectively. Thus, all considered null hypotheses are true. 





We would like to stress that the case of continuously distributed survival times corresponds to an infinite-dimensional problem and is thus more difficult than the discrete case. For example, the simulation study in <cit.> confirms this observation in a related problem: the convergence rate of the actual coverage probabilities of confidence bands 
to the nominal confidence level is much faster the more discretely distributed the survival data is. Moreover, to make the simulation scenario even more challenging, we considered the situation with infinite τ to also get an indication of the functionality of the test in this case.




Simulations.
  We chose as wild bootstrap multipliers centered unit Poisson variables because a formal Edgeworth expansion and two simulation studies in <cit.> indicate
  that  those have theoretical and practical advantages over the common choice of standard normal multipliers.
We chose the nominal level α=5% and conducted each test 10,000 times for K=1,2,3 and 5,000 times for K=5,10 because of the massively increasing computational complexity for large samples.
Each test was based on critical values that were found using 1,999 wild bootstrap iterations. 
All simulations were conducted with the help of the R computing environment <cit.>.

Results.
The true type-I-error results for the four different null hypotheses are shown in Table <ref> (left panel: one-way for H_0^p(C̱_1) and right panel: two-way for 
H_0^p(C̱_2,A)) and Table <ref> (two-way for H_0^p(C̱_2,B) in the left and H_0^p(C̱_2,AB) in the right panel). 
It is apparent that all simulated  levels are elevated for the smallest sample sizes (K=1),
especially for the one-way test: here almost all type I error probabilities are between 13.0% and 17.7%.
For the two-way tests, these probabilities are mainly between 8.1% and 11.7% in this case (K=1).
On the one hand, this is due to the relatively strong censoring rates:
for λ = 0.4, the censoring probabilities across all sample groups are between 33% and 37% (found by simulating 100,000 censoring and survival time random variables each), for λ = 0.5, these probabilities range from 39.5% to 41.5%, and for λ = 2/3, they even reach values between 48.5% and 49%; resulting in only 5 to 7 uncensored observations per group. On the other hand, not to restrict the time horizon in inferential procedures about survival functions appears to slightly slow down the convergence of type I error probabilities to the nominal level as the sample size increases; 
see <cit.> for similar findings in the context of confidence bands for unrestricted survival functions.
However, the error probabilities recover for samples of double size (i.e. between 20 and 28) already:
in the one-way design, these error rates drop to mainly 8.2% - 9.9%, 
and in all two-way tests, we even achieve rates of mainly 6.1% - 8%.
If the sample sizes are tripled (i.e. between 30 and 42), most of the type I error probabilities are between 7%-8% (one-way) or 5.2%-6.9% (two-way).
In case of the sample size factor K=5, all results are only slightly liberal,
and for K=10 (i.e. sample sizes between 100 and 140), we see that the nominal level is well attained.
              









 §.§ Behaviour under shift alternatives


In addition to the simulations of the previous subsection, we also conducted a small power simulation of the above tests. For the alternative hypotheses, we considered a shift model: taking the same six basic survival and censoring functions as in the first set of simulations, we shift all survival and censoring times of the first sample group by δ∈{0.1, 0.2, …, 1}.
In this way, we maintain the same censoring rates as before and the distance to the null hypotheses is gradually increased: for growing δ >0, we obtain a growing relative effect p_1 > 0.5 (one-way) and 
p_11 > 0.5 (two-way), respectively. For each of the above considered contract matrices, C̱_1, C̱_2,A, C̱_2,B, C̱_2,AB,
we conducted one set of simulations with different unbalanced sample sizes and censoring rate combinations. 
For each set-up, we increased the sample sizes by the factors K=1,3,5. The results are displayed in Figure <ref>.

We see that, even for the smallest sample sizes (between 10 and 14), the power of the two-way testing procedures increase to 0.5 or 0.6 as the shift parameter approaches 1. For larger samples sizes the theoretically proven consistency is apparent. In comparison, the one-way test has a much higher power: 
For the undersized case (K=1) it already reaches a power of 0.8 while for moderate to larger sample sizes the power is almost 1 for shift parameters δ≥ 0.5. In comparison to the two-way procedure its superior power is, however, partially paid at the price of its pronounced liberality; especially for small sample sizes.


All in all, the simulations confirm that all tests have a satisfactory power with increasing sample size and/or shift parameter while maintaining a reasonable control of the nominal level for sample sizes of 30 to 42 already.


  


  

§ DATA EXAMPLE

  
  We illustrate the developed theory on a dataset from a colon cancer study <cit.>. 
  Considering the patients in Stage C, that is, there had been metastases to regional lymph nodes,
  the data consist of eligible 929 patients suffering from colon cancer. Survival (measured in days) was the primary endpoint of the study. We focus on the two factors `gender' and `treatment' (with three levels) to obtain a crossed 2× 3 survival design which is in line with a setting from our simulation study. In particular, there were
  315 patients in the observation group,
  310 others were treated with levamisole,
  and 304 received levamisole, combined with fluorouracil.
  Levamisole was originally used as an anthelmintic drug and fluorouracil (5-FU) is a medicine to treat various types of cancer. 
  The patients in the study had been randomized into one of these three treatment groups.
  Also, there were nearly as many women (445) as men (484) involved in the study.
  Figure <ref> depicts the Kaplan-Meier estimates of the survival probabilities for each treatment × sex subgroup.
  We refer to <cit.> for more details about the study.
  The dataset is freely accessible via the  command  after having loaded the package condSURV <cit.>. 
  
  The aim is now to investigate the presence of main or interaction effects of treatment and gender. As there are several ties in the data (roughly 16%; see Appendix <ref> for details) and we do not want to impose specific distributional assumptions, we focus on the nonparametric concordance effects. To this end, we first have to choose a proper τ. From our retrospective view, the most reasonable choice is found by determining for each group the minimal observed censoring time that exceeds all observed survival times in that group. We call these censoring times “terminal times”. Then,  τ is set to be the minimal terminal time.   In doing so, the group with that minimum terminal time does not benefit nor does it suffer from having the earliest terminal time when compared to the other groups.
  
  The first block in Table <ref> shows the sample sizes of all subgroups.
  In the present data example, the minimal terminal time is τ = 2173; see the second block of Table <ref>.
  In view of the sample sizes and the censoring rates given in the third block of Table <ref>,
  we compare the present dataset with the simulation set-ups in Section <ref>: a
  similarly strong censorship is obtained for λ_3 and comparable sample sizes n ∈ [100,140] for the choice K=10. Thus, judging from the rightmost columns of Tables <ref> and <ref>, we find it is safe to assume actual type I error probabilities of about 5.1% to 5.9% of the proposed nonparametric one- and two-way survival tests. 
  
  

  
  
  

  We tested the data in one- and two-factorial set-ups and chose α = 5% as the significance level.
  As in the simulation study, we used B=1,999 bootstrap iterations for each test.
  For the tests in the two-factorial model, we considered the null hypotheses corresponding to no main treatment effect, no main effect in sex, and no interaction effect between both.
  The test results, by means of p-values, are shown in Table <ref>.
  
  


  We found a significant indication against the equality of all d=6 groups (p-value < 0.001) but this difference between groups could not be inferred to result from a difference between the sexes (p-value = 0.331).
  However, we found a significant treatment effect (p-value < 0.001)
  as well as a significant interaction effect between treatment and sex (p-value < 0.001).
  The p-values in Table <ref> have not been adjusted for a type I error multiplicity but it is obvious that the results remain the same after an application of, say, the Bonferroni-Holm procedure.
  
  Indeed, looking at the rightmost block of Table <ref>,
  we agree with the findings of the hypothesis tests:
  the gender effect seems to be cancelled out if the treatment groups are combined,
  but within the male gender there seems to be a big difference in the concordance effects (p_1i_2∈ [0.459, 0.581]).
  Also, the interaction effect is apparent,
  as the female groups do not seem to strongly benefit from any treatment (p_2i_2∈ [0.483, 0.501]).
  On the other hand, the male groups exhibit a worse than average survival fitness in the observation and the levamisole treatment group (p_11=0.475, p_12=0.459)
  but a much better than average survival fitness for the combination treatment (p_13 = 0.581). Here the value 
  p_13 = 0.581 roughly means that a randomly generated observation from this specific group survives 
  a randomly generated observation from the mean distribution of all groups with probability 58.1%. Taking another look at the Kaplan-Meier curves in Figure <ref>,
  we immediately see that our concordance effects and the test outcomes make sense.
  We clearly see that there is a big difference in the male survival probabilities
  (the combination treatment group is superior to the levamisole treatment group which is in turn superior to the observation group)
  but there is not much of a difference between the female groups' survival curves.
  Indeed, comparing the Kaplan-Meier curve of the pooled males' survival times with that of the pooled females' times, 
  we graphically find no evident main gender effect.
  The plot of both Kaplan-Meier estimators is shown in Appendix <ref>.
  
  Finally, we relate our results to the original findings of <cit.>
  whose analyses involve the Cox proportional hazards model and logrank tests.
  The authors also detected that “Therapy with levamisole plus fluorouracil produced an unequivocal advantage over observation” and that levamisole alone did not produce a detectable effect.
  Furthermore, they concluded from an exploratory subset analysis that the “levamisole-fluorouracil treatment appeared to have the greatest advantage among male patients [...]”.
  This is exactly what we confirm in our analysis based on the nonparametric two-factorial tests.
  However, <cit.> neither account for the present ties in the data (note that the Cox regression postulates continuous outcomes) nor did they clearly stress the rather weak effect of the  levamisole-fluorouracil treatment for women. 
  They just state that their “results show [...] striking contradictions to those of subset analyses reported in the NCCTG study, in which levamisole plus fluorouracil was found to be most effective in reducing the risk of recurrence among female patients [...]” among other subgroups of patients.
  
  
  
  

§ DISCUSSION

  
We proposed novel nonparametric inference procedures for the analysis of factorial survival data that may be subject to independent random right-censoring. Critical values are obtained from a multiplier wild bootstrap approach leading to 
asymptotically valid tests and confidence regions for meaningful effect parameters. Thereby, the procedures do not require any multiplicative or additive hazard structure nor specific distributional survival and censoring assumptions.  In particular, different group distributions are allowed and 
ties are accounted for accordingly. Moreover, different to the nonparametric survival procedures of <cit.> and <cit.>, our methods are not only driven towards hypothesis testing but also to uncertainty  quantification of the underlying effect estimators. The latter can be used to comprehensibly describe and infer main and interaction effects in general nonparametric factorial survival designs with an arbitrary number of fixed factors. 
Together with a 1-1 connection with hazard ratios in proportional two-sample designs <cit.>,
this makes the new methods appealing for practical purposes. 

To investigate their theoretical properties, we rigorously proved central limit theorems of the underlying statistics and consistency of the corresponding procedures. In addition, 
extensive simulations were conducted for one- and two-way designs to also assess their finite sample properties in terms of power and type-I-error control. In case of small sample sizes with less than 10 completely observed subjects per group, they revealed a liberal behaviour; especially for the one-way testing procedure. However, for moderate to larger sample sizes the asymptotic results kicked in and the stated theoretical results were recovered.

Finally, the methods were used to exemplify the analysis of survival data in a study about treatments for colon cancer patient within a two-factorial survival design. As severe ties were present in the data, classical hazard based methods were not directly applicable. In comparison, our newly proposed nonparametric methods provided a very decent alternative for the analysis of such factorial survival designs without postulating any strict assumptions. 


To allow for a straightforward application, it is planned to implement the procedure into an easy to use R-package. 
In future research we will consider the case of stochastically ordered subgroups, for which a multiple testing algorithm could be developed with the aim to detect significantly different collections of all subgroups:
subgroups with no significant differences in the nonparametric concordance effects may be combined to facilitate the interpretation of the outcomes
and to ultimately serve for the development of different, more personalized medicines, one for each new subgroup combination.
Moreover, extensions of the current methodology to ordered alternatives or factorial designs obtained via stratified sampling will be part of a practically useful consecutive testing procedure.



  

§ ACKNOWLEDGEMENTS

   Markus Pauly likes to thank for the support from the German Research Foundation (Deutsche Forschungsgemeinschaft).

   plainnat
  
  
  
  

§ APPENDIX

  
  

§ PROOFS

  
  
  

 §.§ Proof of Lemma <ref>

  
  The proof of this lemma is based on an extension of the results of <cit.>.
  In that article it is shown that the functional
  ϕ: (BV[0,τ])^2 →,  (f,g) ↦ - ∫_0^τ f^±(u) g̣(u)
  is Hadamard-differentiable in functions of bounded variation; see the proof of their Theorem 1 in the supplement to that article for more details.
  Hence, a similar result transfers to the functional ϕ in the present article because ϕ is simply a multivariate aggregation of functionals of the same type as ϕ.
  
  For the sake of completeness, we point out that the ((i-1)d+j)th entry of continuous and linear Hadamard-derivative of ϕ at S̱ is given by
 
    (ϕ̣_S̱·ẖ)_i,j = ∫ h_j^±Ṣ_i - ∫ h_i^±Ṣ_j,   ẖ = (h_1, …, h_d)' ∈ (D[0,τ])^d,

 if h_1(0) = … = h_d(0) = 0 which indeed is the case in the present application. 
  
  The other fundamental result that we require is the convergence of all normalized Kaplan-Meier estimators in distribution.
  It is well-known that √(n_i) (  S_i - S_i)  U_i on the Skorokhod space D[0,τ] as n_i →∞ in the present situation of independent right-censoring.
  Hence, due to the independence of all sample groups, it also follows that
  diag(n_1, …, n_d) · (S̱ - S̱) U̱ on (D[0,τ])^d as min(n_1, …, n_d) →∞.
  
  Thus, after having added the multiplicative term √(n_i)/√(n_i) in each component i=1, …, d,
  it follows by means of the functional delta-method (cf. Theorem 3.9.4 in ) 
  that √(N) ( ϕ(S̱) - ϕ(S̱)) converges in distribution as asserted.
  

 §.§ Proof of Theorem <ref>

  
  This convergence result follows immediately from Lemma <ref> in combination with the continuous mapping theorem.
  

 §.§ Proof of Lemma <ref>


  It is well-known that each Kaplan-Meier estimator S_i
  and also each covariance function estimator Γ_i is uniformly consistent on [0,τ] for S_i and on [0,τ]^2 for Γ_i, respectively, i=1, …, d.
  Note here that due to the assumption P(T_i1 > τ) > 0 for all groups i the cumulative hazard functions are bounded on the interval of interest: Λ_i(τ) < ∞.
  
  Since V̱_N is a continuous functional of these estimators,
  it follows from the continuous mapping theorem that the matrix V̱_N is consistent for V̱.
  Likewise, we conclude that the estimator tr(ṮV̱_N) is consistent for the trace tr(ṮV̱) due to continuity of the trace and linear maps. 

 §.§ Proof of Theorem <ref>

  
  Again, the result essentially follows from the continuous mapping theorem as p̱↦p̱' Ṯp̱ is a continuous map.
  Thus, taking Lemma <ref> into account, an application of Slutsky's lemma reveals that
  
    F_N(C̱) d⟶χ

  under H_0^p(C̱) as N →∞.
  
  To show that χ has unit expectation,
  we note that the trace in the denominator of F_N(C̱) accounts for the asymptotic expectation of the quadratic form term:
  the expectation of a quadratic form X̱' A̱X̱ for a random vector X̱ with mean μ = E(X̱) and covariance matrix Σ = cov( X̱ ) equals
  E(X̱' A̱X̱ ) = tr( A̱Σ ) + μ' A̱μ;
  cf. Corollary 3.2b.1 in  <cit.>.
  In our situation, where A̱ = Ṯ and Ṯp̱ = 0̱ under H_0^p(C̱),
  the asymptotic expectation of N p̱' Ṯp̱ under H_0^p thus simplifies to tr( Ṯ V ).
  
  
  On the other hand, under the alternative hypothesis H_a^p(C̱), we have that 
  
    F_N(Ṯ) =  N/tr(ṮV̱_N)p̱' Ṯp̱ =  N/tr(ṮV̱_N) (p̱ - p̱)' Ṯ (p̱ - p̱) + N/tr(ṮV̱_N) (2 p̱ - p̱)' Ṯp̱.

  The first term on the right-hand side converges in distribution to χ as N →∞.
  Of the second term, (2 p̱ - p̱)' Ṯp̱ converges to the positive number p̱' Ṯp̱.
  Thus, an application of Slutzky's lemma yields that F_N(Ṯ) p→∞ as min n →∞.
  

 §.§ Proof of Theorem <ref>

  
  This conditional central limit theorem holds
  because for each wild bootstrapped Kaplan-Meier estimator a conditional central limit theorem holds separately; cf. <cit.> who suggested a wild bootstrap resampling technique for independently right-censored competing risks and survival data in the presence of ties.
  They are combined with the help of the continuous mapping theorem 
  and also with the consistency of the wild bootstrapped covariance estimators which were shown in that same paper.
  The consistency proofs of these methods, which were developped for competing risks data, 
  immediately translate into 
  
  the conditional weak convergence in probability of W̱_N in probability towards diag(λ)^-1/2U̱. Hence, the continuous mapping theorem implies that the conditional distribution of F_N^*(Ṯ) given the data converges weakly in probability to the distribution of χ.
  

§ ASYMPTOTIC COVARIANCE MATRIX V̱ OF √(N) ( P̱ - P̱)

  
  
  The covariance matrix V̱ in Theorem <ref> has the compact form
  
    V̱ =    1/d^2∑_j=1^d 1/λ_j∫∫Γ_j^±±̱̣S ̱̣S'
        + ∫∫ diag ( (1/λ_iΓ_i^±±)_i=1^d ) ̣̅S ̣̅S 
        - 1/d∫∫(1/λ_iΓ_i^±±)_i=1^d ̣̅S ̱̣S'
        - ( 1/d∫∫(1/λ_iΓ_i^±±)_i=1^d ̣̅S ̱̣S' )' 
    
       =    1/d^2∫∫1̱_d' diag(λ)^-1Γ^±±̱̣S ̱̣S'
        + diag(λ)^-1∫∫ diag(Γ^±±) ̣̅S ̣̅S 
        - 1/d∫∫ diag(λ)^-1Γ^±±̣̅S ̱̣S'
        - ( 1/d∫∫ diag(λ)^-1Γ^±±̣̅S ̱̣S' )' ,

  where Γ = (Γ_i)_i = 1, …, d 
  is the vectorial aggregation of asymptotic covariance functions of all Kaplan-Meier estimators.
  
  
  
  

§ ADDITIONAL INFORMATION ON THE DATA ANALYSIS

  
  
  
  
  As stated in Section 6 of the paper there are severe ties present in the dataset. The explicit numbers of ties are presented in Table <ref> below. They show that 130 different observations occurred at least twice in the data set. In particular, only 781 different realizations are encountered among all 929 individuals leading to 
a ties rate of roughly 148/929≈ 16%.              


After choosing the truncation time τ=2173 days as described in the paper, the numbers of ties further increase as illustrated in Table <ref>.
This leads to a ties rate of approximately 439/929≈ 47% in the final data analysis.



Finally, we present in Figure <ref> the Kaplan-Meier curves with respect to gender, thus after a combination of all three male subgroups and all three female subgroups.
Even though the test for the main effect in gender was not detected to be significant in the set-up with six subgroups, we found that gender became highly significant if the two-sample problem is considered (p-value < 0.001).
Indeed, in Figure <ref> the survival times of male patients appear to be slightly stochastically larger than those of the females.
We thus see that the tests further gain in power when different subgroups are combined as the involved sample sizes increase a lot. 
Of course, one should be careful with this new result because combining subgroups may not always make sense. In this case, especially in view of the significant interaction effect between treatment and gender, it could be reasonable to treat women differently than men.



  
