A class of two-sample nonparametric statistics for binary and time-to-event outcomes

Introduction

In many clinical studies, two or more endpoints are investigated aiming to provide a comprehensive picture of the treatment’s benefits and harms. Survival analysis has often been the sharp focus of clinical trial research. However, when there is more than one event of interest, the time until the appearance of the event is not always the unique center of attention; often the occurrence of an event over a fixed time period is as well an outcome of interest.

In the context of cancer immunotherapies trials, short-term binary endpoints based on the tumor size, such as objective response, are common in early-phase trials, whereas overall survival remains the gold standard in late-phase trials.^1,2 Since traditional oncology endpoints may not capture the clinical benefit of cancer immunotherapies, the idea of looking at both tumor response and survival has grown from the belief that together may achieve a better characterization of the clinical response. ³

Several authors have considered both objective response and overall survival as primary endpoints in cancer trials. Lai and Zee ⁴ proposed a single-arm phase II trial design with a tumor response rate and a time-to-event outcome, such as overall survival or progression progression-free survival. In their design, the dependence between the binary response and the time-to-event outcome is modeled through a Gaussian copula. Lai et al. ⁵ proposed a two-step sequential design in which the response rate and the time to the event are jointly modeled. Their approach relates the response rate and the time to the event by means of a mixture model build on the basis of the Cox proportional hazards model assumption. Chen and Wang ⁶ presented a joint model for binary marker responses and survival outcomes for clustered data. They based the statistical inference on a multivariate penalized likelihood method and estimated the standard errors using a jackknife resampling method.

An additional challenge in immunotherapy trials lies in the fact that delayed effects are likely to be found, bringing the need of alternative methods accounting for the non-proportionality of the hazards. ⁷ Statistics that look at differences between integrated weighted survival curves, such as those defined by Pepe and Fleming^8,9 and extended by Gu et al. ¹⁰ , are better suited to detect early or late survival differences and do not depend on the proportional hazards assumption. In this paper, we aim to propose a class of two-sample statistics that could be used in phase II/III design to jointly evaluate the efficacy on binary and survival endpoints, even in the presence of delayed treatment effects.

The problem of how to analyze multiple outcomes has been widely discussed in the literature.^11,12 Depending on the inferential goals, this problem can be addressed in different ways. When the trial’s objective is to test the null hypothesis of no effect on any of the endpoints, one can use the univariate tests together with a multiplicity adjustment to control the inflation of the type I error. The most popular multiple testing procedures (e.g., Bonferroni procedure ¹³ ) do not make any assumption regarding the joint distribution of the univariate tests, and therefore may lead to conservative designs when the endpoints are correlated. When aiming at demonstrating that the treatment has an effect across the endpoints without necessarily specifying the minimum effect on any of the endpoints, one can consider methods based on a combination of the univariate statistics, which generally require assumptions about the joint test distribution.

Several approaches have been developed allowing for the joint distribution of test statistics. O’Brien ¹⁴ and Pocock et al. ¹⁵ proposed global test statistics through the sum of individual statistics. O’Brien ¹⁴ developed a generalized least squares method by combining multiple statistics into a single hypothesis test when variables are normally distributed; whereas Pocock et al. ¹⁵ extended O’Brien’s work to asymptotically normal test statistics. Hothorn et al. ¹⁶ and Pipper et al. ¹⁷ approached the problem of testing multiple hypotheses using parametric and semi-parametric models. Hothorn et al. ¹⁶ used the limiting distribution of the parameter estimators to build upon the corresponding test statistics and their joint distribution. Based on that, their approach corrects the significance level by means of the simultaneous asymptotic normality of the test statistics. Pipper et al. ¹⁷ proposed a procedure for evaluating the efficacy in trials with multiple endpoints of different types. Their procedure is based on simultaneous asymptotic normality of the effect estimators from the single models for each endpoint together with multiple testing adjustments.

Extensive research has been done on joint modeling of longitudinal measurements and survival data (comprehensive overviews can be found in Tsiatis and Davidian, ¹⁸ Rizopoulos ¹⁹ and Papageorgiou et al. ²⁰ ). In most cases, the primary focus is on characterizing the association between the longitudinal and event time processes. The common framework is to relate the time-to-event and longitudinal outcomes through the proportional hazard model. Nevertheless, the relationship between binary response at a specific time point and survival outcome has received less attention. ⁶

In this paper, we have followed the idea launched by Pocock et al. ¹⁵ of combining multiple test statistics into a single hypothesis test. Specifically, we propose a class of statistics based on a weighted sum of a difference in proportions test and a weighted Kaplan–Meier test-based on the difference of survival functions to evaluate an overall effect on the binary and survival endpoints. Our proposal adds versatility to the study design by enabling different follow-up periods for each endpoint, and flexibility by incorporating weights. We define these weights to specify unequal priorities to the different endpoints and to anticipate the type of time-to-event difference to be detected.

This article is organized as follows. In Section ‘A general class of binary and survival test statistics’, we present the class of statistics for binary and time-to-event outcomes. In Section ‘Large sample results’, we set out the assumptions and present the large sample distribution theory for the proposed statistics. In Section ‘On the choice of weights’, we introduce different weights and discuss their choice. We give an overview of our R package SurvBin in Section ‘Implementation’ and illustrate our proposal with a recent immunotherapy trial in Section ‘Example’. In Section ‘Simulation study’, we evaluate the performance of these statistics in terms of the significance level and the statistical power with a simulation study. We conclude with a discussion.

All the required functions to use these statistics have been implemented in R and have been made available at: https://github.com/MartaBofillRoig/SurvBin.

A general class of binary and survival test statistics

Consider a randomized controlled trial comparing two treatment groups, control group ( i = 0 ) and intervention group ( i = 1 ), each composed of n ( i ) individuals, and denote by n = n ( 0 ) + n ( 1 ) the total sample size. Suppose that both groups are followed over the time interval [ 0 , τ ] and are compared on the basis of the following two endpoints: the occurrence of an event ε b before τ b ( 0 < τ b ≤ τ ), and the time to a different event ε s within the interval [ τ 0 , τ ] ( 0 ≤ τ 0 < τ ). For the i -th group ( i = 0 , 1 ), let p ( i ) ( τ b ) be the probability of having the event ε b before τ b , and S ( i ) ( ⋅ ) be the survival function of the time to the event ε s .

We consider the problem of testing simultaneously H b , 0 : p ( 0 ) ( τ b ) = p ( 1 ) ( τ b ) and H s , 0 : S ( 0 ) ( t ) = S ( 1 ) ( t ) , ∀ t ∈ [ τ 0 , τ ] , aiming to demonstrate an overall effect across the binary and survival endpoints, either by a higher probability of the occurrence of ε b , H b , 1 : p ( 0 ) ( τ b ) < p ( 1 ) ( τ b ) , or an improved survival with respect to ε s in the intervention group, H s , 1 : S ( 0 ) ( t ) ≤ S ( 1 ) ( t ) , ∀ t ∈ [ τ 0 , τ ] , S ( 0 ) ( ⋅ ) ≠ S ( 1 ) ( ⋅ ) . The hypothesis problem can then be formalized as follows:

{ H 0 : p ( 0 ) ( τ b ) = p ( 1 ) ( τ b ) and S ( 0 ) ( t ) = S ( 1 ) ( t ) , ∀ t ∈ [ τ 0 , τ ] H 1 : p ( 0 ) ( τ b ) < p ( 1 ) ( τ b ) or S ( 0 ) ( t ) ≤ S ( 1 ) ( t ) , ∀ t ∈ [ τ 0 , τ ] , ∃ t * ∈ [ τ 0 , τ ] , S ( 0 ) ( t * ) < S ( 1 ) ( t * )

(1)

U n ω ( τ 0 , τ b , τ ; Q ^ ) = ω b ⋅ U b , n ( τ b ) σ ^ b + ω s ⋅ U s , n ( τ 0 , τ ; Q ^ ) σ ^ s

(2)

U b , n ( τ b ) = n ( 0 ) n ( 1 ) n ( p ^ ( 1 ) ( τ b ) − p ^ ( 0 ) ( τ b ) )

(3)

U s , n ( τ 0 , τ ; Q ^ ) = n ( 0 ) n ( 1 ) n ( ∫ τ 0 τ Q ^ ( t ) ⋅ ( S ^ ( 1 ) ( t ) − S ^ ( 0 ) ( t ) ) d t )

(4)

The weights ω control the relative relevance of each outcome—if any—and the random weight function Q ^ ( ⋅ ) serves two purposes: to specify the type of survival difference that may exist between groups and to stabilize the variance of the difference of the two Kaplan–Meier functions. Some well-known special cases of Q ^ ( ⋅ ) are follows:

Q ^ ( t ) = G ^ ( t − ) , where G ^ ( t − ) is the pooled Kaplan–Meier estimator for the censoring distribution. This choice of Q ^ ( t ) down-weights the contributions on those times where the censoring is heavy.

The statistics in the L -class are defined for possible different follow-up configurations based on different choices of: the overall follow-up period, τ ; the time where the binary event is evaluated, τ b ; and the origin time for the survival outcome, τ 0 ; taking into account that 0 < max { τ 0 , τ b } < τ . There are however no restrictions on whether or not these periods overlap and, if they do, how much and when. We illustrate two different situations with different configurations for τ 0 , τ b , τ in Figure 1. The first case is exemplified by an HIV therapeutic vaccination study where safety-tolerability response (binary outcome) and time-to-viral rebound (survival outcome) are outcomes of interest. Whereas the safety–tolerability is evaluated at week 6 ( τ b = 6 ), the time-to-viral rebound is evaluated from week 6 to week 18 ( τ 0 = 6 and τ = 18 ). ²² The second example in the area of immunotherapy trials includes a binary outcome (objective response), evaluated at month 6, and overall survival, evaluated from randomization until year 4 ( τ 0 = 0 , τ b = 0.5 and τ = 4 ). ²³

OPEN IN VIEWER

Figure 1.

Illustration of two different follow-up configurations, the red and blue arrows represent the time-frame for binary and time-to-event outcomes, respectively. The red line goes from the start of the study (at time-point 0 ) until the binary outcome is evaluated at time τ b . The blue (dashed) line goes from when the time-to-event information begins to be collected ( τ 0 ) to the end of the study ( τ ).

The L -class statistics includes several statistical tests. If τ 0 = 0 , τ b = τ and ω b = ω s , then, U n ω ( Q ^ ) corresponds to the global test statistic proposed by Pocock et al. ¹⁵ . If ε b = ε s , τ 0 = τ b , and ω b = ω s , the statistic U n ω ( Q ^ ) is the equivalent of the linear combination test of Logan et al. ²⁴ when there is no censorship until τ b for testing for differences in survival curves after a pre-specified time-point.

Large sample results

In this section, we derive the asymptotic distribution of the L -class of statistics given in (2) under the null hypothesis and under contiguous alternatives, present an estimator of their asymptotic variance, and discuss the consistency of the L -statistics against any alternative hypothesis of the form of H 1 in (1). We start the section with the conditions we require for the L -class of statistics. To make the paper more concise and more readable, proofs and technical details are in the Supplemental material. Moreover, stratified L -tests are also presented and discussed in the Supplemental material.

Asymptotic distribution

To derive the asymptotic distribution of the statistic U n ω ( Q ^ ) , we use that U n ω ( Q ^ ) can be approximated by U n ω ( Q ) , the same statistic with the weights replaced by its deterministic function (see Lemma 1 in the Supplemental material). Roughly speaking, thanks to this approximation, we can ignore the randomness of Q ^ ( ⋅ ) and use U n ω ( Q ) to obtain the limiting distribution of U n ω ( Q ^ ) . In what follows, we state the asymptotic distributions under the null hypothesis in Theorem 1 and under a sequence of contiguous alternatives in Theorem 2.

Theorem 1

Let U n ω ( Q ^ ) be the statistic defined in (2). Under the conditions outlined in ‘Further Notation and Assumptions’, if the null hypothesis H 0 : H s , 0 ∩ H b , 0 holds, U n ω ( Q ^ ) converges in distribution, as n → + ∞ , to a normal distribution as follows:

U n ω ( Q ^ ) → N ( 0 , ω b 2 + ω s 2 + 2 ω b ω s ⋅ σ b s σ b ⋅ σ s )

where σ b 2 , σ s 2 stand for the asymptotic variances of U b , n and U s , n ( Q ) , respectively, and σ b s is the covariance between U b , n and U s , n ( Q ) . Their corresponding expressions are given by:

σ b 2 = ∑ i = 0 , 1 ( 1 − π ( i ) ) p ( i ) ( τ b ) ( 1 − p ( i ) ( τ b ) )

(5)

σ s 2 = − ∑ i = 0 , 1 ( 1 − π ( i ) ) ∫ τ 0 τ ( K τ ( i ) ( t ) ) 2 ( S ( i ) ( t ) ) 2 G ( i ) ( t ) d S ( i ) ( t )

(6)

σ b s = ∑ i = 0 , 1 ( 1 − π ( i ) ) ⋅ ( I { τ max = τ b } ⋅ ∫ τ 0 τ b K τ b ( i ) ( t ) S ( i ) ( t ) ⋅ ( p N ( i ) ( t ) − p ( i ) ( τ b ) ) d S ( i ) ( t ) + ∫ τ max τ K τ ( i ) ( t ) S ( i ) ( t ) ⋅ p ( i ) ( τ b ) ( d S X ( i ) ( t ) − − d S ( i ) ( t ) ) )

(7)

where τ max = max ( τ 0 , τ b ) , K τ * ( i ) ( t ) = ∫ t τ * Q ( u ) S ( i ) ( u ) d u ( τ * = τ or τ b ), p N ( i ) ( t ) = P ( X i j = 1 | d N i j ( t ) = 1 ) , and S X ( i ) ( t ) = P ( T i j > t | X i j = 1 ) for i = 0 , 1 .

Recall that σ b 2 , σ s 2 , and σ b s depend on τ 0 , τ b , τ , but we omit them for notational simplicity.

Theorem 2

Let U n ω ( Q ^ ) be the statistic defined in (2). Under the conditions outlined in the “Further notation and assumptions” section, consider the following sequences of contiguous alternatives for both binary and time-to-event hypotheses satisfying, as n → + ∞ :

n ( p n ( 1 ) − p ( 0 ) ) → g

and

n ( S n ( 1 ) ( t ) − S ( 0 ) ( t ) ) → G ( t )

for some constant g ∈ R + and bounded function G ( ⋅ ) ∈ R + , and ∀ t ∈ [ τ 0 , τ ] . Then, under contiguous alternatives of the form:

H 1 , n : n ( p n ( 1 ) − p ( 0 ) ) = g and n ( S n ( 1 ) ( t ) − S ( 0 ) ( t ) ) = G ( t ) , ∀ t ∈ [ τ 0 , τ ]

(8)

we have that:

U n ω ( Q ^ ) → N ( ω b g + ω s ∫ τ 0 τ Q ( t ) G ( t ) d t , ω b 2 + ω s 2 + 2 ω b ω s σ b s σ b ⋅ σ s )

in distribution as n → + ∞ , where σ b 2 , σ s 2 and σ b s are given in (5)–(7), respectively.

The covariance in (7) involves the conditional probabilities S X ( i ) ( t ) and p N ( i ) ( t ) , while S X ( i ) ( t ) represents the survival function for responders –individuals that have had the binary event ε b –, p N ( i ) ( t ) stands for the probability of being a responder among individuals experiencing ε s at t . Also note that, if τ b < τ 0 , the survival experience starts after the binary event has been evaluated and only involves the second integral in (7).

We notice that the efficiency of the L -statistics, U n ω ( Q ^ ) , under contiguous alternatives is driven by the non-centrality parameter μ c = ω b g + ω s ∫ τ 0 τ Q ( t ) G ( t ) d t , that is, by the sum of the weighted non-centrality parameters of U b , n and U s , n ( Q ^ ) .

Variance estimation and consistency

We now describe how to use the L -statistics to test H 0 versus H 1 given in (1). In particular, we propose a consistent estimator of the asymptotic variance of U n ω ( Q ^ ) , and present the standardized L -statistics to test H 0 : H s , 0 ∩ H b , 0 .

The asymptotic variance of U n ω ( Q ^ ) , given in Theorem 1, can be consistently estimated by:

Var ^ ( U n ω ( Q ^ ) ) = ω b 2 + ω s 2 + 2 ω b ω s σ ^ b s σ ^ b ⋅ σ ^ s

(9)

where σ ^ b , σ ^ s , and σ ^ b s denote the estimates of σ b , σ s and σ b s , and are given by:

σ ^ b 2 = p ^ ( τ b ) ( 1 − p ^ ( τ b ) )

(10)

σ ^ s 2 = − ∫ τ 0 τ ( K ^ τ ( t ) ) 2 S ^ ( t ) S ^ ( t − ) ⋅ n ( 0 ) G ^ ( 0 ) ( t − ) + n ( 1 ) G ^ ( 1 ) ( t − ) G ^ ( 0 ) ( t − ) G ^ ( 1 ) ( t − ) d S ^ ( t )

(11)

σ ^ b s = − ∫ τ 0 τ b K ^ τ b ( t ) ( ∑ i = 0 , 1 n − n ( i ) n ⋅ λ ^ X , T ( i ) ( t ) d t + p ^ ( τ b ) ⋅ d S ^ ( t ) S ^ ( t ) ) + ∫ τ b τ K ^ τ ( t ) ⋅ p ^ ( τ b ) S ^ ( t − ) ( − S ^ ( t − ) ⋅ d S ^ ( t ) S ^ ( t ) + ∑ i = 0 , 1 n − n ( i ) n ⋅ S ^ X ( i ) ( t − ) ⋅ d S ^ X ( i ) ( t ) S ^ X ( i ) ( t ) )

(12)

where K ^ τ * ( t ) = ∫ t τ * Q ^ ( u ) S ^ ( u ) d u ( τ * = τ or τ b ), S ^ ( t ) is the pooled Kaplan–Meier estimator of the survival functions, p ^ ( τ b ) is the pooled estimator of the probabilities p ( i ) ( τ b ) , S ^ X ( i ) ( t ) is the Kaplan–Meier estimator of S X ( i ) ( t ) , and λ ^ X , T ( i ) ( t ) is the estimator of λ X , T ( i ) ( t ) = lim d t → 0 P ( X i j = 1 , t ≤ T i j < t + d t | T i j > t ) / d t .

The variance estimator presented in (9) is obtained assuming that the variances of the two groups are equal (pooled estimator). An unpooled variance estimator is proposed in the Supplemental material. For both pooled and unpooled estimators, smoothing techniques are used to estimate λ X , T ( i ) ( t ) over the time period [ τ 0 , τ b ] . In this work, we have chosen kernel smoothing methods. Note that resampling methods can also be used to get an estimator of the variance of U n ω ( Q ^ ) . In the simulation section, we will discuss the results using the pooled, unpooled, and bootstrap variance estimators.

To test the global null hypothesis H 0 : H s , 0 ∩ H b , 0 in (1), we consider the normalized statistic of U n ω ( Q ^ ) :

U n ω ( Q ^ ) / Var ^ ( U n ω ( Q ^ ) )

(13)

Because this statistic (13) converges in distribution to a standard normal distribution, it can be used to test H 0 : H s , 0 ∩ H b , 0 by comparing U n ω ( Q ^ ) / Var ^ ( U n ω ( Q ^ ) ) to a standard normal distribution. Moreover, for positive Q ( ⋅ ) , the statistic is consistent against any alternative hypothesis of the form of H 1 in (1).

On the choice of weights

An important consideration when applying the statistics proposed in this paper is the choice of the weight functions. The L -class of statistics involves the already mentioned random weight function Q ^ ( t ) and deterministic weights ω = ( ω b , ω s ) . These weights are defined according to different purposes and have different roles in the statistic U n ω ( Q ) . In this section, we include different weights and discuss some of their strengths as well as shortcomings. The list provided is not exhaustive, other weights are possible and might be useful in specific circumstances.

Choice of Q ^ ( ⋅ )

The choice of Q ^ ( ⋅ ) might be very general as long as Q ^ ( ⋅ ) converges in probability to a function Q ( ⋅ ) , and both Q ^ ( ⋅ ) and Q ( ⋅ ) satisfy the conditions outlined in the ‘Further notation and assumptions’ section. In this section, we center our attention on a family of Q ^ ( ⋅ ) weights of the form:

Q ^ ( t ) = f ^ ( t ) ⋅ v ^ ( t ) ,

where: (i) f ^ ( ⋅ ) is a data-dependent function that converges, in probability to f ( ⋅ ) , a nonnegative piecewise continuous function with bounded variation on [ 0 , 1 ] . The term f ^ ( t ) takes care of the expected differences between survival functions and can be used as well to emphasize some parts of the follow-up according to the time-points ( τ 0 , τ b , τ s ); (ii) the weights v ^ ( ⋅ ) converge in probability to a deterministic positive bounded weight function v ( ⋅ ) . The main purpose of the weight v ^ ( t ) is to ensure the stability of the variance of the difference of the two Kaplan–Meier functions. To do so, we make the additional assumption that:

| v ( t ) | ≤ Γ ⋅ G ( i ) ( t ) 1 / 2 + δ and | v ^ ( t ) | ≤ Γ ⋅ G ^ ( i ) ( t ) 1 / 2 + δ

for all t ∈ [ τ 0 , τ ] , i = 0 , 1 and for some constants Γ , δ > 0 .

Different choices of f ^ ( t ) yield other known statistics. For instance, if f ( ⋅ ) = 1 , U s , n ( Q ^ ) corresponds to the Weighted Kaplan–Meier statistics.^8,9 Whenever f ^ and v ^ correspond to the weights (15) and (14), respectively, introduced below, we have the statistic proposed by Shen and Cai. ²⁵ Furthermore, note that the weight functions of the form Q ^ ( t ) = f ^ ( t ) ⋅ v ^ ( t ) are similar to those proposed by Shen and Cai ²⁵ ; while they assume that f ^ is a bounded continuous function, we assume that f ^ ( ⋅ ) is a nonnegative piecewise continuous function with bounded variation on [ 0 , 1 ] , and instead of only considering the Pepe and Fleming weight function corresponding to (15), we also allow for different weight functions v ^ ( t ) . Finally, if we do not consider any weight, that is, if Q ^ ( t ) = 1 , ∀ t , U s , n ( Q ^ ) corresponds to the difference of restricted mean survival times from τ 0 to τ .

In what follows, we outline different choices of v ^ ( t ) and f ^ ( t ) , together with a brief discussion for each one:

We require v ^ ( t ) to be small towards the end of the observation period if censoring is heavy. The usual weight functions v ^ ( t ) involve Kaplan–Meier estimators of the censoring survival functions. The most common weight functions are as follows:

v ^ c ( t ) = n G ^ ( 0 ) ( t − ) G ^ ( 1 ) ( t − ) n ( 0 ) G ^ ( 0 ) ( t − ) + n ( 1 ) G ^ ( 1 ) ( t − )

(14)

and v ^ ( t ) = v ^ c ( t ) , both proposed by Pepe and Fleming. ⁸ Among other properties, v ^ c ( ⋅ ) has been proved to be a competitor to the log-rank test for the proportional hazards alternative. ⁸ Note that if the censoring survival functions are equal for both groups and the sampling design is balanced ( n ( 0 ) = n ( 1 ) ), then, the differences in Kaplan–Meier estimators are weighted by the censoring survival function, that is, w ( t ) = C ( t ) = C ( i ) ( t ) for i = 0 , 1 . Also note that w ( t ) = 1 for uncensored data.

Analogously to Fleming and Harrington ²¹ statistics, f ^ ( t ) could be used to specify the type of expected differences between survival functions. That is, if we set:

f ( S ^ ( t − ) ) = S ^ ( t − ) ρ ( 1 − S ^ ( t − ) ) γ , ρ , γ ≥ 0

(15)

the choice ρ > 0 , γ = 0 leads to a test to detect early differences, while ρ = 0 , γ > 0 leads to a test to detect late differences; and ρ = γ = 0 leads to a test evenly distributed over time and corresponds to the weight function of the log-rank.

To put more emphasis on those times after the binary follow-up period we might consider:

f ( t ) = { a , t < τ b 1 − a , t ≥ τ b

for a < 0.5 .

Implementation

We have developed the SurvBin package to facilitate the use of the L -statistics and is now available on GitHub (https://github.com/MartaBofillRoig/SurvBin). The SurvBin package contains two key functions: lstats to compute the standardized L -statistic, U n ω ( Q ^ ) / Var ^ ( U n ω ( Q ^ ) ) , using the variance estimator given in Section ‘Large sample results’; and lstats_boots to compute the standardized L -statistic by using a bootstrap procedure to estimate the variance.

The SurvBin package also provides the functions survbinCov to calculate σ ^ b s ; and bintest and survtest to compute the univariate binary and survival statistics (3) and (4), U b , n ( τ b ) / σ ^ b and U s , n ( τ 0 , τ ; Q ^ ) / σ ^ s , respectively. In addition, the SurvBin package includes the function simsurvbin that can be used to simulate bivariate binary and survival data in a variety of situations.

The main function lstats can be called by: (time, status, binary, treat,tau0, tau, taub, rho, gam, eta, wb, ws, var_est) where time, status, binary and treat are vectors of the right-censored data, the status indicator, the binary data and the treatment group indicator, respectively; tau0, tau, taub denote the follow-up configuration; wb, ws are the weights ω ; rho, gam, eta are scalar parameters that controls the weight Q ^ ( t ) , which is given by Q ^ ( t ) = G ^ ( t − ) η ⋅ S ^ ( t − ) ρ ⋅ ( 1 − S ^ ( t − ) ) γ ; and var_est indicates the variance estimate to use (pooled or unpooled).

In this work, we estimate λ X , T ( i ) ( t ) by means of the Epanechnikov kernel function, and the local bandwidth selection and the boundary correction described by Muller and Wang ²⁶ by using the muhaz package. ²⁷

Example

Melanoma has been considered a good target for immunotherapy and its treatment has been a key goal in recent years. Here we consider a randomized, double-blind, phase III trial whose primary objective was to determine the safety and efficacy of the combination of a melanoma immunotherapy (gp100) together with an antibody vaccine (ipilimumab) in patients with previously treated metastatic melanoma. ²³ Despite the original endpoint was objective response rate at week 12, it was amended to overall survival and then considered a secondary endpoint. A total of 676 patients were randomly assigned to receive ipilimumab plus gp100, ipilimumab alone, or gp100 alone. The study was designed to have at least 90 % power to detect a difference in overall survival between the ipilimumab-plus-gp100 and gp100-alone groups at a two-sided α level of 0.05 , using a log-rank test. Cox proportional-hazards models were used to estimate hazard ratios and to test their significance. The results showed that ipilimumab with gp100 improved overall survival as compared with gp100 alone in patients with metastatic melanoma. However, the treatment had a delayed effect and an overlap between the Kaplan–Meier curves was observed during the first six months. Hence, the proportional hazards assumption appeared to be no longer valid, and a different approach would have been advisable.

To illustrate our proposal, we consider the comparison between the ipilimumab-plus-gp100 and gp100-alone groups based on the overall survival and objective response as multiple primary endpoints of the study. For this purpose, we have reconstructed individual observed times by scanning the overall survival Kaplan–Meier curves reported in Figure 1A of Hodi et al. ²³ using the reconstructKM package ²⁸ (see Figure 2), and, afterwards, we have simulated the binary response to mimic the percentage of responses obtained in the study.

OPEN IN VIEWER

Figure 2.

Kaplan–Meier curves for overall survival for ipilimumab-plus-gp100 and gp100-alone groups (arms 1 and 0, respectively).

Using the data obtained, we employ the L -statistic by means of the function lstats in the SurvBin package. To do so, we need to specify the weights ( Q ^ , ω ) to be used, and the time-points ( τ 0 , τ b , τ ). In our particular case, we take τ 0 = 0 , τ b = 0.5 , τ = 4 according to the trial design, choose Q ^ ( t ) = G ^ ( t − ) ⋅ ( 1 − S ^ ( t − ) ) to account for censoring and delayed effects in late times, and ( ω b , ω s ) = ( 0.25 , 0.75 ) to emphasize the importance of overall survival over objective response.

As shown below, the function lstats returns the standardized L -statistic, together with the L -statistic and its standard deviation, and the individual statistics. The value of the L -statistic, U n ω ( Q ^ ) in (2), is 3.24 and is obtained by using the values of U b , n ( τ b ) and U s , n ( τ 0 , τ ; Q ^ ) ( 0.45 and 2.44 , respectively), and σ ^ b and σ ^ s ( 0.24 and 0.66 ). The statistic U n ω ( Q ^ ) / Var ^ ( U n ω ( Q ^ ) ) equals 4.10 and is computed by using the variance estimator in (9) and then by means of σ ^ b and σ ^ b s together with the estimated covariance σ ^ b s ( − 0.0002 ).

Since we obtained U n ω ( Q ^ ) / Var ^ ( U n ω ( Q ^ ) ) = 4.10 > z α = 0.05 , we have a basis to reject H 0 and conclude that there is an overall effect of ipilimumab across the binary endpoint of tumor reduction and overall survival in patients with metastatic melanoma. Note that, in this example, we have been using the pooled variance estimator for the L -statistic. We have also calculated the statistic using the unpooled and bootstrap variance estimators (see Supplemental material) and notice that the results were not substantially different.

We also explored how much the results might vary if we change the weights assigned. Table 1 shows the values of the L -statistic when we modify the weights ω . Moreover, we compared the values of the statistics for the survival endpoint when using the weighted Kaplan–Meier statistics and weighted log-rank statistics. We note no major differences when using weighted Kaplan–Meier tests with different weights, nor in the weighted log-rank tests. This is because the survivals are the same at the beginning of the curves, and therefore putting a higher weight at the end of the curve does not have as much influence as if there was a small effect at the beginning of the trial. On the other hand, we observe that the weighted Kaplan–Meier statistics take higher values than the weighted log-rank tests for all combinations of weights.

OPEN IN VIEWER

Table 1.

Comparison of statistics (at 4 years) using different weights: L -statistics with different combinations of ( ω b , ω s ) ; and weighted Kaplan–Meier (WKM) statistics and weighted log-rank (WLR) statistics with different ( ρ , γ ) .

Statistics	Weights	Standardized test
L -statistics	( ω b , ω s ) = ( 0.25 , 0.75 )	4.11
L -statistics	( ω b , ω s ) = ( 0.5 , 0.5 )	3.93
L -statistics	( ω b , ω s ) = ( 0.75 , 0.25 )	2.94
WKM statistics	( ρ , γ , η ) = ( 0 , 0 , 0 )	3.70
WKM statistics	( ρ , γ , η ) = ( 0 , 1 , 0 )	3.69
WLR statistics	( ρ , γ ) = ( 0 , 0 )	3.13
WLR statistics	( ρ , γ ) = ( 0 , 1 )	2.75

Simulation study

Design

We conducted a simulation study to evaluate our proposal in terms of the statistical power and the type-I error with small sample sizes. We generated bivariate binary and time-to-event data through a copula-based framework and used conditional sampling as described in Trivedi and Zimmer. ²⁹

The parameters used for the simulation (summarized in Table 2) have been the following:

Clayton’s copula with association parameter between the marginal distributions of the binary and time-to-event outcomes equal to θ = 0.001 , 0.51 , 0.91 . These values correspond, respectively, to Spearman’s rank correlation values equal to 0.0002 , 0.32 , 0.45 , which represent increasing associations between the binary and time-to-event outcomes. We have not considered higher values of θ as they do not fulfill the condition that S X ( i ) ( τ ) > 0 ( i = 0 , 1 ).

Weibull survival functions, S b , a ( 0 ) ( t ) = e − ( t / b ) a , with a = 0.5 , 1 , 2 and b = 1 .

Probability of having the binary endpoint p ( 0 ) = 0.1 , 0.3 .; and

Sample size per arm n ( i ) = 250 for i = 0 , 1 and total sample size n = 500 .

The censoring distributions between groups were assumed equal and uniform U ( 0 , c ) with c = 3 .

Two different follow-up configurations were considered for τ 0 < τ b ≤ τ : (i) τ 0 = 0 , τ b = 0.5 , τ = 1 ; and (ii) τ 0 = 0 , τ b = τ = 1 .

We have considered the weights: Q ^ ( t ) = G ^ ( t − ) η ⋅ S ^ ( t − ) ρ ⋅ ( 1 − S ^ ( t − ) ) γ with η = 1 and ρ , γ = 0 , 1 . When simulating under the null hypothesis, we considered ( ω b , ω s ) equal to ( 0.5 , 0.5 ) ; whereas when simulating under the alternative hypotheses, we considered ( ω b , ω s ) equal to ( 0.25 , 0.75 ) , ( 0.5 , 0.5 ) , and ( 0.75 , 0.25 ) .

OPEN IN VIEWER

Table 2.

Scenarios used in the simulation study.

Parameter	Value	Parameter	Value
p ( 0 )	0.1 , 0.3	a	0.5 , 1 , 2
b	1	c	3
θ	0.001 , 0.51 , 0.91	n ( i ) ( i = 0 , 1 )	250
τ b	0.5 , 1	τ	1
ρ , γ	0 , 1	η	1
d	0 , 0.075	HR	0.75 , 1
( ω b , ω s )	( 0.25 , 0.75 ),( 0.5 , 0.5 ),( 0.75 , 0.25 )	t *	0 , 0.5

The simulations under the alternative hypothesis considered four different situations depending on whether there is a treatment effect on both endpoints and the type of difference between the survival curves. Specifically, the following cases were considered:

Effect on both binary and survival endpoints. The effect on the survival endpoint satisfies the proportional hazards assumption, that is, the hazard ratio (HR) between treatment groups is constant over the study duration;

Effect on the binary endpoint and non-effect on the survival endpoint ( H s , 0 : S ( 0 ) ( t ) = S ( 1 ) ( t ) , ∀ t );

Non-effect on the binary endpoint ( H b , 0 : p ( 0 ) ( τ b ) = p ( 1 ) ( τ b ) ) and effect on the survival endpoint with HR constant over the study duration;

Effect on both binary and survival endpoints. The treatment differences on the survival endpoint have a delayed effect, that is, the survival functions are assumed to be equal until time t * , and there is a constant hazard ratio (HR) between treatment groups from t * to τ .

We used d = p ( 1 ) ( τ b ) − p ( 0 ) ( τ b ) = 0.075 to simulate the effects on the binary endpoint. For the survival endpoint, we considered HR = 0.75 under proportional hazards, and HR = 0.70 and t * = 0.5 under delayed effects.

We evaluated the empirical significance level and the statistical power using the L -statistics with pooled, unpooled, and bootstrap variance estimators, and for the sake of the comparison, using the Bonferroni procedure. In addition, we presented the empirical results for testing the individual hypothesis H b , 0 and H s , 0 by using the statistics (3) and (4).

The total number of scenarios was 1504 ( 144 under the null hypothesis and 1360 under the alternative hypothesis). We ran 100 , 000 replicates and estimated the significance level ( α = 0.05 ) for each scenario under the null hypothesis. We ran 1000 replicates and estimated the statistical power for each scenario under alternative hypotheses. We performed all computations using the R software (version 4.0.2).

Power properties

When there is a treatment effect on both endpoints and the proportional hazards assumption is fulfilled (case 1), we obtained empirical powers with medians 0.85 , 0.85 , 0.83 using the L -statistics with pooled, unpooled, and bootstrap variance estimators, respectively; whereas the median of the empirical powers using Bonferroni was 0.80 . When there is a treatment effect on both endpoints and there are delayed effects (case 4), the empirical powers for the L -statistics have medians 0.73 , 0.75 , 0.70 using the pooled, unpooled, and bootstrap variance estimators, respectively; whereas the median of the empirical powers using Bonferroni was 0.63 .

Table 3 summarizes the simulation results on the power across different parameters in case 1. We compared the performance of the different variance estimators and noticed that the empirical powers do not substantially differ between them. We also observed that the power is not affected by the different weight functions in the case of proportional hazards. We obtained higher powers when emphasizing late-differences between the survival curves ( γ = 1 ) in the case of delayed effects (median powers of 0.85 , 0.85 , 0.84 for unpooled, pooled, and bootstrap variance estimators, respectively, with γ = 1 against 0.83 , 0.83 , 0.82 for unpooled, pooled, and bootstrap variance estimators with γ = 0 ).

OPEN IN VIEWER

Table 3.

Median empirical size and median empirical power from 100 , 000 and 1000 replications, respectively. The empirical size and powers are calculated using: the L -statistics (in (2)) according to the pooled, unpooled, bootstrap variance estimators (labeled as Pooled, Unpooled, and Boots.); and the Bonferroni procedure (Bonf.). Under the null hypothesis there is no effect on any of the endpoints ( d = 0 , HR = 1 ). Under the alternative hypothesis there is effect on both endpoints (Case 1: d = 0.075 , HR = 0.75 ) and the effect on the survival endpoint satisfies the proportional hazards assumptions ( t * = 0 ).

		Empirical size				Empirical powers (Case 1)
		Pooled	Unpooled	Boots.	Bonf.	Pooled	Unpooled	Boots.	Bonf.
( τ b , τ )	( 0.5 , 1 )	0.053	0.053	0.050	0.050	0.84	0.85	0.83	0.80
	( 1 , 1 )	0.056	0.055	0.050	0.050	0.85	0.85	0.83	0.79
θ	0.001	0.051	0.051	0.052	0.051	0.87	0.88	0.87	0.82
	0.510	0.054	0.055	0.049	0.050	0.84	0.83	0.81	0.79
	0.910	0.056	0.056	0.048	0.049	0.82	0.82	0.79	0.78
p ( 0 )	0.1	0.053	0.054	0.050	0.051	0.89	0.89	0.88	0.84
	0.3	0.053	0.055	0.050	0.049	0.78	0.78	0.76	0.73
a	0.5	0.053	0.054	0.051	0.050	0.87	0.87	0.86	0.82
	1	0.053	0.054	0.050	0.050	0.84	0.85	0.83	0.80
	2	0.054	0.055	0.048	0.050	0.82	0.81	0.79	0.78
( ρ , γ , η )	(0,0,1)	0.054	0.053	0.051	0.050	0.84	0.85	0.83	0.80
	(0,1,1)	0.054	0.055	0.051	0.050	0.84	0.85	0.83	0.80
	(1,0,1)	0.053	0.055	0.050	0.050	0.83	0.82	0.80	0.79
	(1,1,1)	0.054	0.053	0.050	0.050	0.86	0.86	0.84	0.81

Figure 3 shows boxplots for the empirical powers using the pooled, unpooled, and bootstrap variance estimators and the Bonferroni procedure. These simulations show the superiority of the L -statistics over the Bonferroni procedure, in terms of power, both under proportional hazards and under delayed effects and regardless of the choice of the weights ( ω b , ω s ).

OPEN IN VIEWER

Figure 3.

Boxplot of empirical powers based on scenarios in Table 2. The empirical powers are calculated using: the L -statistics (in (2)) according to the pooled, unpooled, bootstrap variance estimators; the Bonferroni procedure; and the individual statistics (3) and (4). The individual statistics for the binary and survival endpoints are labeled, respectively, as BE and SE. The color indicates which combination of weights ( ω b , ω s ) were used: red for ( ω b = 0.25 , ω s = 0.75 ); blue for ( ω b = 0.5 , ω s = 0.5 ); and green for ( ω b = 0.75 , ω s = 0.25 ).

When there is a treatment effect on only one of the endpoints (cases 2 and 3), the behavior of the power mainly relies on the pre-specified weights ( ω b , ω s ) (see Figure 3). If the survival endpoint is considered clinically more important than the binary endpoint and we use the weights ( ω b = 0.25 , ω s = 0.75 ), then the median of the empirical powers is around 0.42 in case 2 (i.e., when there is treatment effect on the survival endpoint) and around 0.40 in case 3 (i.e., when there is no effect on the survival endpoint). We found a similar behavior when the binary endpoint is considered more important and ( ω b = 0.75 , ω s = 0.25 ).

If both endpoints are equally important and there is treatment effect in only one of them, the empirical powers using the L -statistics take values between the power would have had used the two individual statistics. Given that the Bonferroni procedure assigns more importance to the more highly significant of the endpoints, ¹⁵ the powers are in this case higher using Bonferroni than using L -statistics.

We have also evaluated the empirical powers if we have had an endpoint with a small effect instead of no effect. We observe that the difference in powers between the L -statistics and Bonferroni procedure is smaller, and that the powers using L -statistics could be even higher than the ones using Bonferroni (see Supplemental material). Furthermore, we have assessed how the association between the binary and survival endpoints might affect the results by considering Frank’s copula instead of Clayton’s copula. We notice that there is no substantial difference between the results using Frank’s and Clayton’s copulas. Results using Frank’s copula can be found in the Supplemental material.

Discussion

We have proposed a class of statistics for a two-sample comparison based on two different outcomes: one dichotomous taking care, in most occasions, of short short-term effects, and a second one addressed to detect long long-term differences in a survival endpoint. Such statistics test the equality of proportions and the equality of survival functions. The approach combines a score test for the difference in proportions and a weighted Kaplan-Meier test-based for the difference of survival functions. The statistics are fully non-parametric and α level for testing the null hypothesis of no effect on any of these two outcomes. The statistics in the L -class are appealing in situations when both outcomes are relevant, regardless of how the follow-up periods of each outcome are, and even when the hazards are not proportional with respect to the time-to-event outcome or in the presence of delayed treatment effects, albeit the survival curves are supposed not to cross.

In this work, we focus mainly on clinical trials where the groups to be compared are supposed to only differ with respect to the treatment they receive. However, in observational studies and also sometimes in clinical trials, there is a need to adjust for other covariates than treatment. Whenever the number of covariates is not too large, stratified tests can be useful. The stratified version of the L -test is deferred to the Supplemental material.

We have incorporated weighted functions in the L -statistics to control the relative relevance of each outcome and to specify the type of survival differences that may exist between groups. In our proposed statistics, the weights ( ω b , ω s ) have been defined with the goal of incorporating the potential difference in clinical importance between the binary and survival endpoint, and therefore they must be fixed in the planning stage. As shown in the simulation study, the power of the trial will depend on the trial objectives and then on the relevance of each of the endpoints by means of ( ω b , ω s ) . A sensitivity analysis could be carried out calculating the L -statistics based on a pre-specified set of weights, as it is often done with the weighted log-rank statistics. ³⁰ The extension of these statistics incorporating data-driven weights to maximize the power will be considered in future works.

The testing procedure using the L -class of statistics satisfies a property called coherence that says that the nonrejection of an intersection hypothesis implies the nonrejection of any sub-hypothesis it implies, that is, H s , 0 and H b , 0 . ³¹ However, the testing procedure based on the L -class of statistics does not fulfil the consonant property that states that the rejection of the global null hypothesis implies the rejection of at least one of its sub-hypothesis. Bittman et al. ³² faced the problem of how to combine tests into a multiple testing procedure for obtaining a procedure that satisfies the coherence and consonance principles. An extension of this work to obtain a testing procedure that satisfies both properties could be an important research line to consider. If, however, it is desired to conclude an efficacy claim for the individual endpoints, alternative methods could also be considered, such as an adjusted Bonferroni (taking into account the joint distribution of the univariate tests stated in Section ‘Large sample results’), or a two-stage method in which efficacy is assessed in the univariate tests after detecting an overall effect. ³³

In this paper, we have considered weighted Kaplan–Meier tests-based for comparing the survival curves motivated by recent immunotherapy trials in which the proportional hazards assumption is in doubt. However, if the survival endpoint satisfies the proportional hazards assumption, the combination of a weighted log-rank test and difference in proportion test could also be considered. Previous works by Wei and Lachin^34,35 could then be used for this problem. The comparison of the relative efficiency between tests combining binary and survival data using weighted log-rank tests or weighted Kaplan–Meier tests and differences in proportions under a range of alternative hypotheses is open for future research.

This work has been restricted to those cases in which censoring and survival do not prevent assessing the binary endpoint response. Two questions remain open for future research concerning this point. The first is the possibility of censorship occurring earlier and making it impossible to assess the binary endpoint. The second is the fact that the survival event may preclude observing the binary endpoint. While the second case refers to a competing risk problem, the first refers to a censoring regarding binary endpoint. Regarding the competing-risk problem, it should also be noted that the time-to-event is usually more relevant in practice. And, in this situation, it is generally considered that if the survival endpoint happens before the binary endpoint, the binary endpoint would have taken the value representing the worst outcome. Further research is needed to study the implications of competing risks or censorship on this type of design. We plan to extend the L -class of statistics to more general censoring schemes where the binary endpoint could be censored or where there may be a case of competing risk. Last but not least, extensions to sequential and adaptive procedures in which the binary outcome could be tested at more than one time-point remain open for future research.

Supplemental Material