Optimality of testing procedures for survival data

TL;DR

This paper proposes an optimal testing procedure for survival data that accounts for non-proportional hazards, improving power prediction and serving as a benchmark for alternative tests in clinical trials.

Contribution

It introduces a new test that maximizes Bayesian predicted power conditioned on early data, addressing non-proportional hazards in survival analysis.

Findings

The proposed test outperforms traditional methods under non-proportional hazards.

Simulation studies based on cancer trial data demonstrate the test's effectiveness.

The method provides a benchmark for evaluating other survival analysis tests.

Abstract

Most statistical tests for treatment effects used in randomized clinical trials with survival outcomes are based on the proportional hazards assumption, which often fails in practice. Data from early exploratory studies may provide evidence of non-proportional hazards which can guide the choice of alternative tests in the design of practice-changing confirmatory trials. We study a test to detect treatment effects in a late-stage trial which accounts for the deviations from proportional hazards suggested by early-stage data. Conditional on early-stage data, among all tests which control the frequentist Type I error rate at a fixed $\alpha$ level, our testing procedure maximizes the Bayesian prediction of the finite-sample power. Hence, the proposed test provides a useful benchmark for other tests commonly used in presence of non-proportional hazards, for example weighted log-rank tests. We illustrate the approach in a simulations based on data from a published cancer immunotherapy phase III trial.

Figures

a of Ferris et al. ([2016](#bib.bib13)) by means of the DigitizeIt (TM) software (version 2.2) and the data extraction method of Guyot et al. ([2012](#bib.bib22)). Figure [1(a)](#S2.F1.sf1) shows the resulting Kaplan-Meier curves, which compare survival probabilities between the two study arms. These do not clearly separate in the initial 3-4 months of follow-up, a signal of delayed treatment effects.

a of Ferris et al. ([2016](#bib.bib13)) by means of the DigitizeIt (TM) software (version 2.2) and the data extraction method of Guyot et al. ([2012](#bib.bib22)). Figure [1(a)](#S2.F1.sf1) shows the resulting Kaplan-Meier curves, which compare survival probabilities between the two study arms. These do not clearly separate in the initial 3-4 months of follow-up, a signal of delayed treatment effects.

a of Ferris et al. ([2016](#bib.bib13)) by means of the DigitizeIt (TM) software (version 2.2) and the data extraction method of Guyot et al. ([2012](#bib.bib22)). Figure [1(a)](#S2.F1.sf1) shows the resulting Kaplan-Meier curves, which compare survival probabilities between the two study arms. These do not clearly separate in the initial 3-4 months of follow-up, a signal of delayed treatment effects.

Equations

$$p_{\theta }(x)=\prod _{i=1}^n{r}^{a_i}(1-r{)}^{1-a_i}{h}_{a_i}(t_i;\theta {)}^{1-d_i}{S}_{a_i}(t_i;\theta )g_i(t_i{)}^{d_i}{G}_i(t_i{)}^{1-d_i},$$

$$L(\theta ;x_e)=\prod _{i=1}^{n_e}{h}_{a_{e,i}}(t_{e,i};\theta {)}^{1-d_{e,i}}{S}_{a_{e,i}}(t_{e,i};\theta ),$$

$$BE{P}_{\phi }={\int }_{{\Theta }_1}{\pi }_{\phi }(\theta )p(\theta |x_e)d\theta ,$$

$$BE{P}_{\phi }=\int \int u(\theta ,\phi ,x)p_{\theta }(x)p(\theta |x_e)dxd\theta .$$

$$m(x)={\int }_{{\Theta }_1}L(\theta ;x)p(\theta |x_e)d\theta ,$$

$$0\le BE{P}_{\phi }-BE{P}_{{\phi }^{\prime }}\le f(\alpha ,r,n)=\frac{(1-r{)}^n}{\alpha }\sum _{s=0}^n{\left(\frac{r}{1-r}\right)}^s.$$

$$m(x)=\prod _{a=0}^1\prod _{j=1}^{k+1}{\left(\frac{v_{a,j}+s_{e,a,j}}{v_{a,j}+s_{e,a,j}+s_{a,j}}\right)}^{u_{a,j}+y_{e,a,j}+y_{a,j}}\frac{\Gamma (u_{a,j}+y_{e,a,j}+y_{a,j})}{\Gamma (u_{a,j}+y_{e,a,j})},$$

$$q(x)={\int }_{{\Theta }_1}{p}_{\theta }(x)\frac{p(\theta |x_e)}{P(H_1|x_e)}d\theta ,$$

$$Q(E)=\int \frac{{\sum }_{\sigma }I\left\{m(t,d,a_{\sigma })=m^{(k_{\alpha })}(x)\right\}q(t,d,a_{\sigma })}{{\sum }_{\sigma }q(t,d,a_{\sigma })}dQ(x),$$

$$\frac{{\sum }_{\sigma }I\left\{m(t,d,a_{\sigma })=m^{(k_{\alpha })}(x)\right\}q(t,d,a_{\sigma })}{{\sum }_{\sigma }q(t,d,a_{\sigma })}=\frac{{\sum }_{\sigma }I\left\{m(t,d,a_{\sigma })=m^{(k_{\alpha })}(x)\right\}m(t,d,a_{\sigma })}{{\sum }_{\sigma }m(t,d,a_{\sigma })}$$

$$\le \frac{{\sum }_{\sigma }I\left\{m(t,d,a_{\sigma })=m^{(k_{\alpha })}(x)\right\}m(t,d,a_{\sigma })}{{\sum }_{\sigma }I\left\{m(t,d,a_{\sigma })\ge m^{(k_{\alpha })}(x)\right\}m(t,d,a_{\sigma })}\le \frac{{\sum }_{\sigma }I\left\{m(t,d,a_{\sigma })=m^{(k_{\alpha })}(x)\right\}{m}^{(k_{\alpha })}(x)}{{\sum }_{\sigma }I\left\{m(t,d,a_{\sigma })\ge m^{(k_{\alpha })}(x)\right\}{m}^{(k_{\alpha })}(x)}$$

$$=\frac{{\sum }_{\sigma }I\left\{m(t,d,a_{\sigma })=m^{(k_{\alpha })}(x)\right\}}{{\sum }_{\sigma }I\left\{m(t,d,a_{\sigma })\ge m^{(k_{\alpha })}(x)\right\}}=\frac{\#\{j:m^{(j)}(x)=m^{(k_{\alpha })}(x)\}}{\#\{j:m^{(j)}(x)\ge m^{(k_{\alpha })}(x)\}},$$

$$Q(E)\le E_Q\left[\frac{1}{\left(\genfrac{}{}{0ex}{}{n}{{\sum }_{i=1}^n{a}_i}\right)\alpha }\right]=\sum _{s=0}^n\frac{1}{\left(\genfrac{}{}{0ex}{}{n}{s}\right)\alpha }\left(\genfrac{}{}{0ex}{}{n}{s}\right)r^s(1-r{)}^{n-s}=f(\alpha ,r,n).$$

Full text

Keywords: censored data; decision theory; design of clinical trials; hypothesis testing; proportional hazards.

1 Introduction

Researchers often use data generated by early-phase clinical studies to specify the protocol of randomized confirmatory phase III trials. Data predictive of confirmatory trial outcomes, including early estimates of treatment effects, are used to choose the primary endpoints (Gómez et al., 2014), the sample size (Lindley, 1997), the target populations (Lee and Wason, 2018), and other aspects of the study design (Brody, 2016). Still, in most cases prior information is not used to specify in the protocol, as mandated by regulatory agencies, which hypothesis testing procedure will be used in the final analyses to provide evidence of treatment effects. Agencies such as the U.S. Food and Drug Administration require the control of Type I and II errors at pre-specified rates (US Food and Drug Administration, 1998).

In Phase III trials, standard tests, such as Mantel’s log-rank, are often selected even for studies where prior data suggest their assumptions will be violated (Royston and Parmar, 2013; Alexander et al., 2018). For survival endpoints, methods related to the log-rank test are prevalent. Asymptotically, this is the most powerful test with a proportional hazards alternative (Fleming and Harrington, 2011). However, the proportional hazards assumption is often violated in practice, contributing to false-negative findings (Royston and Parmar, 2013), invalidating sample size calculations (Barthel et al., 2006), and affecting interim analyses (van Houwelingen et al., 2005).

Data from early-stage studies can inform about deviations from the proportional hazards assumption, suggesting the use of alternative methods (Royston and Parmar, 2013). Several extensions and alternatives are available to replace Mantel’s test, such as weighted (Fleming and Harrington, 2011) or adaptive log-rank tests (Yang and Prentice, 2010), and restricted mean survival tests (Royston and Parmar, 2013). Some of these procedures identify the most powerful test against specific alternatives, which may be discordant with early estimates from previous trials. Moreover, their optimality typically holds in a large-sample sense (e.g. in the local limit, for weighted log-rank tests; Fleming and Harrington, 2011).

We develop a statistical test to detect treatment effects in late-stage trials, accounting for deviations from the proportional hazards assumption indicated by early-phase studies (e.g. phase II trials). The proposed test does not belong to the weighted log-rank family or other common classes of tests. Starting from decision theory principles (Berger, 2013), we derive this test as the solution of a constrained decision problem (Ventz and Trippa, 2015): conditional on early-stage data, the test maximizes the predicted finite-sample power among all tests which control the frequentist Type I error rate of the late-stage study at a fixed $\alpha$ level. More precisely, the test maximizes the Bayesian predictive probability that the null hypothesis will be correctly rejected at the end of the confirmatory trial. The test is therefore a useful benchmark for other procedures applicable in the presence of non-proportional hazards.

As a motivating example, we consider the analysis of a randomized trial with delayed treatment effects on survival outcomes. This characteristic occurs when the treatment requires an induction period before it starts to provide therapeutic effects. When treatment effects are delayed, the hazard functions are not proportional and they separate across arms only later during follow-up (Fine, 2007). Initially overlapping survival curves (c.f. Figure 1(a)) are well documented in trials of cancer immunotherapies (Chen, 2013; Alexander et al., 2018). They can also be observed in other settings, such as in studies of breast cancer (Mehta et al., 2012) and melanoma (Robert et al., 2015) chemotherapies.

2 Example

We consider the survival times of the 361 patients with head and neck carcinomas that participated in CheckMate 141 study (Ferris et al., 2016), a Phase III trial that randomized patients to receive nivolumab, a novel cancer immunotherapy, or standard of care (SOC) in a 2:1 ratio. We reconstructed the individual-level data of this trial from Figure 1a of Ferris et al. (2016) by means of the DigitizeIt (TM) software (version 2.2) and the data extraction method of Guyot et al. (2012). Figure 1(a) shows the resulting Kaplan-Meier curves, which compare survival probabilities between the two study arms. These do not clearly separate in the initial 3-4 months of follow-up, a signal of delayed treatment effects.

3 Planning a late-stage trial

We plan a late-stage randomized trial with a survival endpoint and a sample size of $n$ patients. This will generate data $x=(t,d,a)$ to test if the treatment has positive effects on the primary outcome. Here, $t=(t_{1},\ldots,t_{n})$ are the observed follow-up times, $d=(d_{1},\ldots,d_{n})$ are the corresponding censoring indicators ($d_{i}=1$ if $t_{i}$ is censored, while $d_{i}=0$ if an event was observed), and $a=(a_{1},\ldots,a_{n})$ are the treatment arm indicators ($a_{i}=0$ or $a_{i}=1$ if the $i$-th patient is randomized to the control or treatment arm). Patients are assigned to arms with a fixed randomization probability. We assume that censoring times are non informative in the sense of Heitjan and Rubin (1991) and independent of treatment assignment.

For design purposes, we specify a model for the distribution that will generate the data $x$. This is described by a density $p_{\theta}(x)$ that depends on a vector of parameters $\theta\in\Theta$. Here $\theta$ may be infinite-dimensional if the model is semi- or non-parametric. Typically, $p_{\theta}(x)$ will have the form

[TABLE]

where $r\in(0,1)$ is the probability of assignment to arm $a=1$, $h_{a}(t;\theta)>0$ is the hazard function of arm $a=0,1$ (for example, in the exponential model, $h_{a}(t;\theta)=\theta_{a}$, $\theta=(\theta_{0},\theta_{1})\in\Theta=(0,+\infty)^{2}$), $S_{a}(t;\theta)=\exp\left(-\int_{0}^{t}h_{a}(s;\theta)ds\right)$ is the corresponding survival function, and $g_{i}(t)$ and $G_{i}(t)$ are the density and (left-continuous) survival function of the $i$-th patient’s censoring time. Here, the censoring mechanism is taken as known, a common assumption when planning new experiments (Chow et al., 2007). We will later explain that this assumption is not used in the development of the proposed testing procedure.

We consider the non-parametric null hypothesis $H_{0}:P\in\mathcal{P}_{0}$, where $P$ is the true data-generating distribution of $x$ (i.e. $P(A)$ is the probability that $x\in A$) and $\mathcal{P}_{0}$ is the class of all distributions which are invariant with respect to permutations of the treatment arm assignments. Hence, $P\in\mathcal{P}_{0}$ if $P(t,d,a)=P(t,d,a^{\prime})$ for all $a^{\prime}$ obtained by permuting the elements of $a$.

The alternative hypothesis is defined from model (1) as $H_{1}$: $P$ has density $p_{\theta}(x)$ for some $\theta\in\Theta_{1}$, where $\Theta_{1}$ is a subset of $\Theta$. For example, $\Theta_{1}$ may include all $\theta$ such that $h_{0}(t;\theta)\neq h_{1}(t;\theta)$, such that the median of $S_{1}(t;\theta)$ is greater than that of $S_{0}(t;\theta)$, or such that the restricted mean survival time is greater in arm $a=1$ (Royston and Parmar, 2013).

According to the definition of the null hypothesis, when the experimental treatment has no effects – regardless of whether the model $p_{\theta}(x)$ is correct or not – the treatment assignments $a_{1}$, $\ldots$, $a_{n}$ do not provide information about the follow-up times $t$ and censoring indicators $d$. Hence, the distribution of the data does not change if these are arbitrarily permuted (Fisher, 1935; Dawid, 1988; Good, 2006; Pesarin and Salmaso, 2010).

This definition of $H_{0}$ covers distributions $P(x)$ for which $(t_{1},d_{1},a_{1}),\ldots,(t_{n},d_{n},a_{n})$ - the observations from individual patients - are not independent and identically distributed. This may happen when investigators selectively enroll patients based on interim analyses or results from other studies published during the enrollment period. Additional examples include amendments of inclusion-exclusion criteria, or improvements of adjuvant therapies (Tamm and Hilgers, 2014). In these cases, if the experimental treatment has no effects, $P(x)$ remains invariant when the indicators of the assigned arms are permuted.

It is now necessary to choose which $\alpha$-level test $\varphi(x)$ should be used in the late-stage trial. A (randomized) test of $H_{0}$ is a function $\varphi(x)\in[0,1]$. After the data $x$ have been collected, the hypothesis $H_{0}$ is rejected with probability $\varphi(x)$ (Lehmann and Romano, 2006). A test $\varphi(x)$ is non-randomized if it can only attain the values 0 and 1. Only non-randomized tests are used in practice, but here we consider randomized tests because of their analytic advantages. The expected value $E_{P}[\varphi(x)]=\int_{\mathcal{X}_{n}}\varphi(x)dP(x)$ is equal to the probability of rejecting $H_{0}$ with data generated from the distribution $P$. If $\alpha\in(0,1)$ and $E_{P}[\varphi(x)]\leq\alpha$ for all $P\in\mathcal{P}_{0}$, then $\varphi(x)$ is said to have level $\alpha$.

4 Bayesian expected power

Different $\alpha$-level tests are usually compared with respect to their power functions $\pi_{\varphi}(\theta)=\int\varphi(x)p_{\theta}(x)dx$. If $\varphi_{1}(x)$ and $\varphi_{2}(x)$ are two $\alpha$-level tests for $H_{0}$ versus the simple alternative $H_{1}:\theta=\theta_{1}$, for some fixed $\theta_{1}\in\Theta_{1}$, then $\varphi_{1}(x)$ is preferred to $\varphi_{2}(x)$ if $\pi_{\varphi_{1}}(\theta_{1})\geq\pi_{\varphi_{2}}(\theta_{1})$. Such comparisons are difficult for composite alternative hypotheses. In fact, uniformly most powerful $\alpha$-level tests, i.e. tests achieving the maximum power across all alternative models $\theta_{1}\in\Theta_{1}$, may not exist (Lehmann and Romano, 2006).

To address this problem, some authors proposed to compare tests with respect to their average power. Specifically, the average power of a test $\varphi(x)$ is $\int_{\Theta_{1}}\pi_{\varphi}(\theta)p(\theta)d\theta$, where $p(\theta)$ is a distribution weighting each value of $\theta\in\Theta$ based on pre-experimental information (Spiegelhalter and Freedman, 1986; O’Hagan et al., 2005). With this metric, two tests are always comparable. Additionally, $\alpha$-level tests maximizing the average power always exist, although these may be randomized (Chen et al., 2007).

To allow data $x_{e}=(t_{e},d_{e},a_{e})$ from an early-stage trial to inform comparisons between tests, we consider a data-dependent prior $p(\theta|x_{e})$. We assume that the prior probability of the alternative hypothesis $H_{1}$ is positive, i.e. $P(H_{1}|x_{e})=\int_{\Theta_{1}}p(\theta|x_{e})d\theta>0$ (in general it may be $P(H_{1}|x_{e})<1$ because the support of the prior may be larger than $\Theta_{1}$).

Several approaches have been proposed to incorporate historical data in a prior distribution, e.g. power priors (Ibrahim et al., 2015), meta-analytic priors (Schmid et al., 2016), and commensurate priors (Hobbs et al., 2011). The effective sample size metric of Morita et al. (2008) can guide the specification of $p(\theta|x_{e})$ and help avoid overly informative choices.

For simplicity, we define $p(\theta|x_{e})$ as the posterior distribution $p(\theta|x_{e})\propto L(\theta;x_{e})p(\theta)$, where, letting $n_{e}$ be the early-stage trial sample size,

[TABLE]

while $p(\theta)$ is a prior distribution on $\Theta$ whose choice depends on the specific application context. In this definition of $p(\theta|x_{e})$ we implicitly assume homogeneous treatment effects and survival distributions across the early- and late-stage trials.

Extending the average power approach, the Bayesian expected power (BEP) of $\varphi(x)$ is

[TABLE]

a concept first introduced by Brown et al. (1987) and “rediscovered” by several authors (Liu, 2018). It is simple to observe that $BEP_{\varphi}=\textrm{Pr}(\varphi(x)\textrm{ rejects }H_{0}\textrm{ and }\theta\in\Theta_{1}|x_{e})$, the probability, conditional on the early-stage data, that $\varphi(x)$ will correctly reject $H_{0}$ at the end of the late-stage trial. This is often called the probability of success of the trial (Liu, 2018).

From the point of view of decision theory (Berger, 2013), the BEP is the expected value of the utility function $u(\theta,\varphi,x)=I\!\left\{\theta\in\Theta_{1}\right\}\varphi(x)$ (if $H_{1}$ holds, then the utility increases with the probability $\varphi(x)$ of rejecting $H_{0}$). Indeed,

[TABLE]

The problem of choosing which test to apply in the late-stage trial can thus be stated as a constrained maximization problem (Ventz and Trippa, 2015): among $\alpha$-level tests we optimize the BEP.

5 Tests maximizing the expected power

We construct an $\alpha$-level test with maximum Bayesian expected power for the null hypothesis $H_{0}:P\in\mathcal{P}_{0}$, which includes all distributions invariant with respect to permutations of the treatment assignments $a_{1}$, $\ldots$, $a_{n}$. We show that the optimal test can be obtained in the form of a permutation test. Any permutation test is obtained by computing or approximating the distribution of some real-valued test statistic $T(x)$ across all permutations of the treatment assignments, while the values of the follow-up times $t$ and censoring indicators $d$ are kept fixed at the observed values.

To be more formal, for each permutation $\sigma$ of $(1,\ldots,n)$, we denote with $a_{\sigma}=(a_{\sigma(1)}$, $\ldots$, $a_{\sigma(n)})$ the corresponding vector obtained by re-ordering the elements of $a=(a_{1},\ldots,a_{n})$. Moreover, if $T(x)$ is any real-valued statistics, for each $x=(t,d,a)$ we let $T^{(1)}(x)\leq\cdots\leq T^{(n!)}(x)$ be the ordered values of $T(t,d,a_{\sigma})$ as $\sigma$ varies across all $n!$ permutations.

The $\alpha$-level (randomized) permutation test $\varphi(x)$ of $H_{0}$ based on the test statistic $T(x)$ can now be defined as follows. First, let $k_{\alpha}=n!-\lfloor\alpha n!\rfloor$, so that, for each $x$, $T^{(k_{\alpha})}(x)$ is the $(1-\alpha)$-level quantile of $T^{(j)}(x)$ for $j=1,\ldots,n!$. Second, let $M^{+}(x)=\sum_{j=1}^{n!}I\{T^{(j)}(x)>T^{(k_{\alpha})}(x)\}$ and $M^{0}(x)=\sum_{j=1}^{n!}I\{T^{(j)}(x)=T^{(k_{\alpha})}(x)\}$ be the number of $T^{(j)}(x)$’s greater or equal to $T^{(k_{\alpha})}(x)$, respectively. Then, the permutation test $\varphi(x)$ is defined by letting $\varphi(x)=1$ when $T(x)>T^{(k_{\alpha})}(x)$, $\varphi(x)=0$ when $T(x)<T^{(k_{\alpha})}(x)$, and $\varphi(x)=(\alpha n!-M^{+}(x))/M^{0}(x)<1$ when $T(x)=T^{(k_{\alpha})}(x)$. This satisfies $E_{P}[\varphi(x)]=\alpha$ for all $P\in\mathcal{P}_{0}$ (Lehmann and Romano, 2006, Theorem 15.2.1).

In order to construct a maximum-BEP test, we define

[TABLE]

which is the marginal likelihood of the data $x$ given the early-stage data $x_{e}$. Without loss of generality, we assume that the density $p_{\theta}(x)\propto L(\theta;x)$ defined by Equation 1 is determined by a dominating measure invariant with respect to permutations of the treatment assignments.

Proposition 5.1**.**

Given the early-stage data $x_{e}$, the $\alpha$-level permutation test with test statistic $T(x)=m(x)$ maximizes the BEP in Equation (4) among all $\alpha$-level tests of $H_{0}$.

Proof.

The proof is provided in Appendix A.1. ∎

The maximum-BEP test $\varphi(x)$ requires to compute the marginal likelihood $m(x)$ of $x$ and of all permuted versions of the same dataset. Since $m(x)$ does not depend on the censoring distribution functions $G_{i}(t)$ that appear in Equation 1, the censoring mechanism is irrelevant to identify the optimal test. Note, however, that the censoring mechanism still determines the value of the BEP of $\varphi(x)$.

Since randomized tests are not used in practice, we consider the non-randomized version of the test $\varphi(x)$ from Theorem 5.1, i.e. $\varphi^{\prime}(x)=I\!\left\{m(x)>m^{(k_{\alpha})}(x)\right\}$ (this does not depend on $M^{+}(x)$ or $M^{0}(x))$. The non-randomized test $\varphi^{\prime}(x)$ rejects $H_{0}$ when the permutation p-value $\textrm{ppv}(x)=\sum_{\pi}I\!\left\{m(t,d,a_{\pi})\geq m(t,d,a)\right\}/n!$ is less or equal than $\alpha$, where the sum extends over all $n!$ permutations of $a_{1},\ldots,a_{n}$ (Lehmann and Romano, 2006, Section 15.2.1). Typically, $n!$ will be too large to compute $\textrm{ppv}(x)$ exactly. Hence, we implement the test $\varphi^{\prime}(x)$ by means of a Monte Carlo approximation. Accordingly, given $x$, a large random sample of permutations $\pi_{1},\ldots,\pi_{B}$ ($B=10^{3}$, say) is used to compute the estimate $\widehat{\textrm{ppv}}(x)=\sum_{i=1}^{B}I\!\left\{m(t,d,a_{\pi_{i}})\geq m(t,d,a)\right\}/B.$ The hypothesis $H_{0}$ is then rejected if $\widehat{\textrm{ppv}}(x)\leq\alpha$ (Pesarin and Salmaso, 2010, Section 1.9.3).

Since $\varphi^{\prime}(x)\leq\varphi(x)$, the non-randomized test $\varphi^{\prime}(x)$ is $\alpha$-level for $H_{0}$, but it may not achieve the maximum BEP. Still, $\varphi^{\prime}(x)$ is a useful benchmark for other tests of $H_{0}$, as its BEP will be close to that of the maximum-BEP test $\varphi(x)$ for large $n$. In fact, the following Proposition shows that, under mild conditions, the difference between $BEP_{\varphi}$ and $BEP_{\varphi^{\prime}}$ is bounded above by a known function $f(\alpha,r,n)$ such that $f(\alpha,r,n)\rightarrow 0$ as $n\rightarrow+\infty$ for all fixed levels $\alpha$ and randomization probabilities $r$.

Proposition 5.2**.**

Suppose that for all $x=(t,d,a)$ with $q(x)>0$ it is $m(x)\neq m(t,d,a_{\pi})$ for all permutations $\pi$ such that $a\neq a_{\pi}$. Then

[TABLE]

Proof.

The proof is provided in Appendix A.2. ∎

In general, $f(\alpha,r,n)$ will be fairly small. For example, if $\alpha=0.05$ and $r=1/2$, $f(\alpha,r,n)=2^{-n}(n+1)/\alpha<10^{-3}$ for all $n\geq 15$. For $r\neq 1/2$ it is $f(\alpha,r,n)=[(1-r)^{n+1}-r^{n+1}]/\alpha(1-2r)$; for $r=2/3$ (as in Section 2), $f(\alpha,r,n)<10^{-3}$ for all $n\geq 25$.

6 The piecewise exponential model

Because of its flexibility and tractability, to implement our maximum-BEP test we use a piecewise exponential model (Benichou and Gail, 1990). In this model, the hazard function $h_{a}(t;\theta)$ is constant over the intervals of a fixed partition $\tau_{0}=0<\tau_{1}<\cdots<\tau_{k}<+\infty=\tau_{k+1}$ of the time axis. In particular, $h_{a}(t;\theta)=\theta_{a,j}$ if $t\in[\tau_{j-1},\tau_{j})$, with $j=1,\ldots,k+1$, $t\in\mathbb{R}_{+}$, arms $a=0,1$, and $\theta=(\theta_{0,1},\ldots,\theta_{0,k+1},\theta_{1,1},\ldots,\theta_{1,k+1})\in\Theta=(0,+\infty)^{2(k+1)}$.

The likelihood function $L(\theta;x)$ of the piecewise exponential model depends on a simple set of sufficient statistics. Given data $x=(t,d,a)$, let $s_{a,j}=\sum_{i=1}^{n}\max(0,\min(\tau_{j}-\tau_{j-1},t_{i}-\tau_{j-1}))I\{a_{i}=a\}$ be the total time at risk spent in the interval $[\tau_{j},\tau_{j+1})$ by patients in arm $a$. Additionally, let $y_{a,j}=\sum_{i=1}^{n}(1-d_{i})I\{a_{i}=a,\ \tau_{j-1}\leq t_{i}<\tau_{j}\}$ be the number of events observed during $[\tau_{j-1},\tau_{j})$ in arm $a$. Then, $L(\theta;x)=\prod_{a=0}^{1}\prod_{j=1}^{k+1}\theta_{a,j}^{y_{a,j}}\exp(-\theta_{a,j}s_{a,j}).$

For convenience, we use a conjugate prior $p(\theta)$. This is obtained by letting all $\theta_{a,j}$ be independent and distributed as a gamma random variable with shape parameter $u_{a,j}$ and rate parameter $v_{a,j}$. With this choice, the distribution $p(\theta|x_{e})$ presents independent $\theta_{a,j}$ components which are gamma distributed with shape parameter $u_{a,j}+y_{e,a,j}$ and rate parameter $v_{a,j}+s_{e,a,j}$, where the $y_{e,a,j}$ and $s_{e,a,j}$ are the sufficient statistics of $x_{e}$. The marginal likelihood $m(x)$ needed to implement the maximum-BEP test can thus be obtained explicitly from Equation 5:

[TABLE]

where $\Gamma(z)$ is the gamma function.

As an example, Figure 1(a) shows the posterior means of the survival probabilities in the nivolumab or SOC arm of CheckMate 141 obtained from the piecewise exponential model. For all $j=1,\ldots,k=4$, we conveniently defined $\tau_{j}$ to be the $j$-th quintile of the distribution of follow-up times in the SOC arm ($\tau_{1}=2.2$, $\tau_{2}=4.0$, $\tau_{3}=5.3$, and $\tau_{4}=7.9$ months). In other words, the prior model is chosen by peeking at the early stage trial. Additionally, we specify gamma priors on the $\theta_{a,j}$ with $u_{a,j}=v_{a,j}=10^{-3}$ for all $a$ and $j$. The posterior estimates (Figure 1(a)) reflect the delayed separation in the Kaplan-Meier curves, as the estimated survival probabilities diverge only after 4 months of follow-up.

7 Application: trials with delayed treatment effects

7.1 Simulation study

As an illustration, we use CheckMate 141 data to simulate a large number of phase II and III trials with delayed treatment effects. In these simulations, we compare different tests with respect to their probability of rejecting the hypothesis of no treatment effects at the end of the phase III trial. We consider Mantel’s log-rank test and several others that account for delayed treatment effects: i) a lagged log-rank test that ignores the first $10\%$ of observed follow-up times (Zucker and Lakatos, 1990); ii) the Fleming-Harrington $G^{0,1}$ test, which gives more weight to late events (Fine, 2007); iii) the adaptive log-rank of Yang and Prentice (2010), which weighs events according to a preliminary estimate of the hazard functions; and iv) a test of the difference in Restricted Mean Survival Times (RMSTs) across study arms (we follow Huang and Kuan, 2018, and compute RMSTs up to the smallest of the two maximum event times observed in each study arm). We also implement the maximum-BEP test (using the Monte Carlo approach of Section 5) based on phase II data. For all tests, we consider $\alpha=0.05$ and a two-sided alternative hypothesis.

To simulate a trial of size $n$, we first sample with replacement $n$ patients from the CheckMate 141 data. Then, depending on the individual treatment assignment, we generate the corresponding survival times from the Kaplan-Meier curves of Figure 1(a). Assuming a maximum follow-up of 15 months, we generate patient’s censoring times by sampling independently from the empirical censoring distribution (Efron, 1981).

Using this approach, we iterate the following steps 10,000 times: i) we simulate a phase II trial of about half the size of CheckMate 141 ($n_{e}=180$); ii) using the simulated phase II data $x_{e}$, we compute the sufficient statistics of the piecewise exponential model (i.e. all $s_{e,a,j}$ and $y_{e,a,j}$; c.f. Section 6); we fix the $\tau_{j}$s at the quintiles of the follow-up times in the SOC arm from $x_{e}$; iii) we simulate a subsequent phase III study, generating phase III data $x$ of size $n=361$; iv) we compute the marginal likelihood $m(x)$ of the phase III data $x$ (c.f. Equation 6; we used gamma parameters $u_{a,j}=v_{a,j}=10^{-3}$); v) we apply our permutation test, as well as the others described above. The proportion of rejections across iterations is the Monte Carlo estimate of a test’s rejection probability.

Figure 1(b) reports the estimated rejection probabilities for each testing procedure. The maximum-BEP permutation test based on phase II data had the highest estimated probability of rejecting the null hypothesis (about 0.90). The $G^{0,1}$ test and the lagged log-rank test have both estimated rejection probabilities of about 0.87. The adaptive log-rank and RMST tests have lower rejection probabilities, 0.77 and 0.66 respectively. Mantel’s log-rank test has the lowest estimated rejection probability, i.e. 0.60, a third less than the one achieved by our test. This finding is consistent with previous studies, which highlighted how the log-rank test may suffer a severe loss of power when treatment effects are delayed (Fine, 2007; Chen, 2013; Alexander et al., 2018).

7.2 Robustness analysis

We consider four additional simulation scenarios in which the outcome distributions in phase II and III are not identical. In all scenarios, the distribution of the phase II data $x_{e}$ are the same as in Section 7.1, while the distribution of the phase III data $x$ are different. In Scenario 1, the data $x$ are generated from the predictive distribution $q(x)$ (see Equation 7): a value $\theta^{\prime}$ is first sampled from $p(\theta|x_{e})$, then $x$ is generated from the distribution $p_{\theta^{\prime}}(x)$. Here we assume $r=2/3$ as in CheckMate 141 and that censoring can only occur after 15 months of follow-up. Proposition 5.1 indicates that, in this scenario, our permutation test has the highest expected power. Scenarios 2-4 instead represent two settings in which our test may suffer from a loss of power. In Scenario 2, $x$ is generated by a different piecewise exponential model than the one used to construct the benchmark test. The phase III delay in treatment effects is shorter than in phase II data. Specifically, $x$ is generated by a model with only one cut-point, fixed at $\tau_{1}=2$ months, and whose parameters are set equal to the maximum likelihood estimates obtained from CheckMate 141 data. Scenario 3 is similar, but the cut-point is fixed at $\tau_{1}=8$ months to represent longer phase III delays than in phase II. Finally, in Scenario 4 there are no delays in treatment effects (i.e. no cut-points) in phase III and the standard proportional hazards assumption holds. Supplementary Figure 4 shows the phase III survival curves used in each scenario.

Figure 2 shows the results of this robustness analysis. As expected, in Scenario 1 our permutation test has a much higher rejection probability than all other tests (0.98). Instead, its performance is sub-optimal in Scenarios 2-4. Although in Scenario 3 our permutation test may be considered comparable with the others (rejection probability equal to 0.29), in Scenario 2 it has the lowest rejection probability (0.84, compared to 0.90 for the Mantel’s log-rank). The results from Scenario 4 highlight, as expected, that our test based on phase II data can have a lower power than Mantel’s log-rank if the arm-specific hazards are proportional in the phase III study. All these findings support the intuition that the power of the maximum-BEP test depends on how accurately the phase III survival distributions can be predicted from prior phase II data. If the phase II and III trial populations are markedly different, then a test specified using phase II data may perform poorly in the phase III study.

In the Supporting Information, we report on additional robustness analyses. First, we considered an additional simulation scenario where the survival curves cross (Supplementary Figure 5(b), panel a). In this scenario, the power of our permutation test was one of the highest among the compared tests (Supplementary Figure 5(b), panel b). In addition, we assessed how the power of the proposed test varied with the size of the early-phase trial. To do so, we repeated the simulations of Section 7.1 for different values of $n_{e}$ between 40 to 220 patients. The test’s power was $\geq 80\%$ for all considered values of $n_{e}$, and $\geq 90\%$ for $n_{e}\geq 160$ (Supplementary Figure 6).

8 Generalization to stratified designs

Treatment effects are often expected to vary across patients’ groups defined, for example, by gender or biomarkers. In such cases one can stratify patients with respect to covariates measured before randomization. We focus on the primary goal of testing whether the experimental treatment has no effects across all strata or if it is effective at least in some of the strata (alternative hypothesis), for example in one or multiple subgroups defined by a relevant biomarker (Freidlin et al., 2010).

Our approach can be easily generalized to this setting. For simplicity, we suppose each patient $i=1,\ldots,n$ is categorized by a binary covariate $z_{i}=0,1$, presence ($z_{i}=1$) or absence ($z_{i}=0$) of some marker. Data $x$ become $x=(t,d,a,z)$, where $z=(z_{1},\ldots,z_{n})\in\{0,1\}^{n}$. We assume that censoring is non-informative and independent of treatment assignments conditionally on $z_{1},\ldots,z_{n}$ (Heitjan and Rubin, 1991).

To illustrate, we specify a piecewise-exponential model $h_{a}(t;\theta,z)$ for the hazard function in arm $a=0,1$ for patients with marker level $z=0,1$ (Freidlin et al., 2010): $h_{a}(t;\theta,z)=\theta_{a,z,j}$ for all $t\in[\tau_{j-1},\tau_{j})$. The prior remains nearly identical to the previous sections. In particular, the marginal likelihood $m(t,d,a,z)$, similar to Equation 5, has a closed form expression.

We specify the null hypothesis $H_{0}^{\prime}:P\in\mathcal{P}^{\prime}_{0}$, where $\mathcal{P}^{\prime}_{0}$ is the class of all distributions which are invariant with respect to permutations of the treatment assignment $a$ within strata of $z$. More precisely, $P\in\mathcal{P}^{\prime}_{0}$ if and only if $p$, the density of $P$, satisfies $p(t,d,a,z)=p(t,d,a_{\sigma},z)$ for all permutations $\sigma$ of $(1,\ldots,n)$ such that $z_{\sigma(i)}=z_{i}$ for all $i=1,\ldots,n$.

With a simple modification, Proposition 5.1 still holds with this new definition of the null hypothesis. In Section 5, the maximum-BEP permutation test computes the distribution of $m(x)$ under $H_{0}$ by considering all permutations of the treatment arm indicators $a_{1},\ldots,a_{n}$. In the stratified case, only permutations $\sigma$ of $(1,\ldots,n)$ such that $z_{\sigma(i)}=z_{i}$ for all $i=1,\ldots,n$ are considered. If $\Sigma(z)$ is the set of all such $\sigma$, then the permutation p-value associated to the maximum-BEP test is given by $\textrm{ppv}(x)=\sum_{\sigma\in\Sigma(z)}I\{m(t,d,a_{\sigma},z)\geq m(t,d,a,z)\}/|\Sigma(z)|$, where $|\Sigma(z)|=(\sum_{i=1}^{n}z_{i})!(n-\sum_{i=1}^{n}z_{i})!$.

8.1 Simulation study of biomarker-stratified designs

Similarly as in Section 7.1, we use CheckMate 141 data to simulate phase II and III trials where 50% of patients express ($z_{i}=1$) a marker predictive of treatment effects (Patel and Kurzrock, 2015). In the simulations, we compare different tests for the final analysis of the phase III trial, each accounting for patients’ maker values. Specifically, we consider four (5%-level) tests of $H_{0}^{\prime}$: i) a test based on a stratified Cox proportional hazards model (a common approach in this setting; Mehrotra et al., 2012); ii) the procedure obtained by first performing separate log-rank tests within marker strata and then applying the Bonferroni correction (another common approach; Freidlin et al., 2014); iii) the test obtained by estimating separate log-normal accelerated failure time (AFT) models (Kalbfleisch and Prentice, 2002, Chapter 2) within marker strata, testing the statistical significance of their regression coefficients, and then applying the Bonferroni correction; and, lastly, iv) our maximum-BEP test, tailored to the simulated phase II data.

In detail, we simulate 10,000 phase II ($n_{e}=180$) and III ($n=361$) trials from CheckMate 141 data for each of 3 scenarios. In Scenario 1, the survival distributions in phase II and III trials are the same. Here, the survival time of a patient in arm $a=0$ or with marker $z=0$ (respectively, in arm $a=1$ with marker $z=1$) is generated from the SOC (nivolumab) Kaplan-Meier curve in Figure 1(a). Censoring times are generated as in Section 7.1.

In Scenarios 2 and 3, phase II data are generated as before, but the outcome distribution in phase III is different. In Scenario 2, phase III data are generated assuming proportional hazards within levels of $z$. Specifically, the survival time of a patient in arm $a=0$ or with marker $z=0$ (respectively, in arm $a=1$ with $z=1$) is simulated from an exponential distribution, fixed at the corresponding estimate from the SOC (respectively, nivolumab) arm of CheckMate 141. Scenario 3 is similar, but we do not assume proportional hazards in the strata of $z$. Rather, phase III data are generated as in Scenario 2, but using stratum-specific log-normal distributions (Kalbfleisch and Prentice, 2002).

Figure 3 shows the results of the simulations. For Scenario 1, the estimated rejection probabilities are 0.18 for the stratified Cox model, 0.26 for the Bonferroni-based log-rank test, 0.08 for the Bonferroni-based AFT test, and 0.49 for our test tailored to early-stage data. As in Section 7.2, in Scenario 2 the test based on a stratified Cox proportional hazards model and the Bonferroni-based log-rank test outperforms our tailored test, as their estimated rejection probabilities are 0.21, 0.29, and 0.17, respectively. Instead, in Scenario 3, the proposed test has an estimated rejection probability (0.32) higher than that of the other tests, although the power advantage is reduced compared to Scenario 1. We obtained qualitatively the same results when we repeated the simulations assuming a phase III sample size ($n=1,000$) sufficient to achieve 90% power with our test (Supplementary Figure 7). These results confirm that, in late-stage trials, a benefit can be attained when prior data are used to optimize hypothesis testing, but this benefit can be reduced when the early- and late-stage survival distributions are discordant.

9 Discussion

Data from previous studies should be routinely used to design late-stage clinical trials. This is especially relevant when standard assumptions, such as the proportional hazards assumption, might not hold. Our approach allows to specify a test for final analyses that achieves two goals: i) it accounts for the specific type of deviations from proportional hazards suggested by prior data - delayed treatment effects, crossing survival curves, or others (Uno et al., 2014); ii) it satisfies the requirements of regulatory agencies (Ventz and Trippa, 2015). The test maximizes a decision-theoretic criteria, it leverages prior data, and it is of $\alpha$-level for an interpretable null hypothesis.

Our contribution is distinct from the growing body of literature on trial designs that incorporate external data in the study analyses (Psioda and Ibrahim, 2018; Ventz et al., 2019). In contrast, we leverage early-phase data only during the design stage of the late-phase trial. We do not directly use any external data in the final analysis of the late-phase trial. As a consequence, our testing procedure controls the phase III false-positive error rate at the chosen $\alpha$ level, even when early-phase data may overestimate the treatment effect.

The use of appropriate metrics of treatment effects is important to interpret trial results. Rejecting the null hypothesis of equality of survival curves is often not sufficient to understand which treatment is clinically preferable, especially when survival curves cross. To do so, study arms must be compared using interpretable measures of the treatment effects (Saad et al., 2018). For example, in the context of Section 7 it could be relevant to consider the one-sided alternative hypothesis that the 12-months restricted mean survival time in the nivolumab arm is greater than in the standard of care arm (Royston and Parmar, 2013). As discussed in Section 3, our approach could be used to construct testing procedures that maximize the expected power to detect improvements in restricted mean survival time.

To plan a phase III trial using our approach, different designs can be compared in simulation studies consistent with early-phase data Liu (2018). For example, the phase III sample size can be chosen as follows. First, multiple hypothetical phase III trials, each based on a different sample size, are simulated on the basis of early-phase information (e.g. as done in Section 7). Then, the probabilities of detecting treatment effects associated to each sample size are estimated using the simulations. The smallest sample size whose estimated power achieves adequate levels (e.g. $>80\%$) is chosen for the design of the future phase III trial.

To implement our approach, we used a piece-wise exponential model with fixed cut-points chosen on the basis of prior data (Benichou and Gail, 1990). In general, however, it may be necessary to choose between models with a different number of cut-points - e.g. by comparing their fit to early-phase data using Akaike’s information criterion or similar metrics (Gelman et al., 2013, Chapter 7). Models that regularize and shrink the difference between $\theta$ components between adjacent time points could be advantageous, especially with small early-phase datasets. This is because the regularization could reduce the impact of the choice of cut-points (Ibrahim et al., 2005, Chapter 3; Murray et al., 2016). Additionally, investigators can add a prior distribution on the position and the number of cut-points (Demarqui et al., 2008). Other parametric and or non-parametric models different than the piece-wise exponential could also be considered in our approach (Ibrahim et al., 2005). Regardless of the chosen model, to implement our test it is necessary to compute the marginal likelihood of late-stage data. This may be complicated for some non-conjugate models, but several methods have been developed to approximate it (Pajor et al., 2017).

Although we derived our test assuming a single early-stage dataset, the use of multiple prior data sources may provide better outcome predictions for late-stage trials. Our approach can incorporate multiple datasets using power priors (Ibrahim et al., 2000) or hierarchical models (Spiegelhalter et al., 2004).

Our simulations, based on data from the CheckMate 141 trial, confirm that weighted log-rank tests can outperform other tests in the presence of delayed treatment effects. However, these tests depend on a set of tuning parameters, such as the duration of the lag time for lagged log-rank tests or the $\rho$ and $\delta$ coefficients of the $G^{\rho,\delta}$ family, which may be difficult to tune. Instead, our approach directly translates early-stage data into a test procedure for the late-stage trial. Robustness analyses highlight how the performance of our approach depends on the consistency of outcome data and the similarity of enrolled populations between phase II and III trials. Ensuring the transportability of results to subsequent trials remains a major concern in the design of exploratory clinical trials (Wang et al., 2006).

Acknowledgements

We thank G. Fell and T. Chen (Dana-Faber Cancer Institute) for helping with data collection. We also thank Massimiliano Russo (Harvard Medical School) for useful suggestions and Sarah Craver for help in manuscript editing.

Data Availability Statement

The data that supports the findings of this study are available in the supplementary material of this article.

Appendix: technical proofs

A.1 Proof of Proposition 5.1

We will use two technical lemmas. The first shows that for any probability law $Q\not\in\mathcal{P}_{0}$ (c.f. Section 3) it is possible to construct a permutation test which is $\alpha$-level for $H_{0}:P\in\mathcal{P}_{0}$ and has maximum power against the simple alternative $H_{1}:P=Q$. We assume that the densities $p_{\theta}(x)$ are taken with respect to the same dominating measure $\mu$, and that this is invariant with respect to permutations of the treatment assignments.

Lemma A.1**.**

Let $\varphi(x)$ be the $\alpha$-level permutation test of $H_{0}:P\in\mathcal{P}_{0}$ based on $T(x)=q(x)$, where $q(x)$ is the density of $Q\not\in\mathcal{P}_{0}$ with respect to $\mu$. If $\varphi^{\prime}(x)$ is another $\alpha$-level test of $H_{0}$, then $E_{Q}[\varphi^{\prime}(x)]\leq E_{Q}[\varphi(x)]$, i.e. $\varphi(x)$ has higher power under the alternative $H_{1}:P=Q$.

Proof.

Let $\mathcal{P}_{\mu}\subseteq\mathcal{P}_{0}$ be the set of all distributions dominated by $\mu$ that are invariant with respect to permutations of treatment assignment (a non-empty set, since it includes $q^{\prime}(t,d,a)=\sum_{\sigma}q(t,d,a_{\sigma})/n!$). By Theorem 2 of Lehmann and Stein (1949), $E_{Q}[\varphi^{\prime}(x)]\leq E_{Q}[\varphi(x)]$ for every test $\varphi^{\prime}(x)$ such that $E_{P}[\varphi^{\prime}(x)]\leq\alpha$ for all $P\in\mathcal{P}_{\mu}$. Now, if $\varphi^{\prime}(x)$ is an $\alpha$-level test of $H_{0}$, then $E_{P}[\varphi^{\prime}(x)]\leq\alpha$ for all $P\in\mathcal{P}_{\mu}$ and therefore $E_{Q}[\varphi^{\prime}(x)]\leq E_{Q}[\varphi(x)]$. ∎

The second lemma shows that to construct an $\alpha$-level maximum-BEP test of $H_{0}$ it is sufficient to construct an $\alpha$-level test of $H_{0}$ that has maximum power against a simple alternative $H_{1}:P=Q$ defined by a specific $Q\not\in\mathcal{P}_{0}$. More precisely, here $Q$ is the the predictive distribution of $x$ conditional on $\theta\in\Theta_{1}$, which is defined by

[TABLE]

Lemma A.2**.**

A test $\varphi(x)$ of $H_{0}$ maximizes the BEP (c.f. Equation 3) among all $\alpha$-level tests if and only if it maximizes the power $E_{Q}[\varphi(x)]$ among all $\alpha$-level tests.

Proof.

By Fubini’s theorem, the BEP of a test $\varphi(x)$ can be written as $BEP_{\varphi}=E_{Q}[\varphi(x)]\cdot P(H_{1}|x_{e})$. Since $P(H_{1}|x_{e})$ is positive and does not depend on $\varphi(x)$, $BEP_{\varphi}$ is maximal if and only if $E_{Q}[\varphi(x)]$ is also maximal. ∎

We are finally ready to prove Proposition 5.1.

Proof of Proposition 5.1.

By Lemma A.1, the $\alpha$-level permutation test $\varphi^{\prime}(x)$ based on $T^{\prime}(x)=q(x)$ maximizes the power $E_{Q}[\varphi(x)]$ among all $\alpha$-level tests of $H_{0}$. By Lemma A.2, $\varphi^{\prime}(x)$ thus has maximum BEP among all $\alpha$-level tests of $H_{0}$. It now suffices to show that $\varphi^{\prime}(x)=\varphi(x)$ for all $x$ such that $q(x)>0$, where $\varphi(x)$ is the $\alpha$-level permutation test based on $T(x)=m(x)$. To do so, note that, by Equation 1, if $q(x)>0$, then $m(x)>0$ as well, and the ratio $q(x)/m(x)$ is invariant with respect to permutations of the treatment arm assignments. The thesis now follows since $q(t,d,a_{\sigma})\propto m(t,d,a_{\sigma})$ for all permutations $\sigma$. ∎

A.2 Proof of Proposition 5.2

Denote with $\Pi(x)$ the set of all $\binom{n}{\sum_{i=1}^{n}a_{i}}$ distinct datasets obtainable from $x=(t,d,a)$ by permuting the elements of $a$. Let $Q$ be defined by the density $q(x)$ in Equation 7. By assumption, if $q(x)>0$ then $m(x_{1})\neq m(x_{2})$ for all $x_{1},x_{2}\in\Pi(x)$, $x_{1}\neq x_{2}$. By Proposition A.2, $0\leq BEP_{\varphi}-BEP_{\varphi^{\prime}}=E_{Q}[\varphi(x)-\varphi^{\prime}(x)]\leq Q(E)$, where $E$ is the set of all $x$ such that $m(x)=m^{(k_{\alpha})}(x)$. From Section 5.9 of Lehmann and Romano (2006),

[TABLE]

where both sums extend over all $n!$ permutations $\sigma$ of $(1,\ldots,n)$. If $q(x)>0$, then

[TABLE]

where the first equality follow because the ratio $q(x)/m(x)$ is invariant with respect to permutations $\sigma$ of $a$. Now, by the definitions of $m^{(k_{\alpha})}(x)$ and $k_{\alpha}$, the denumerator of the last fraction is greater or equal than $\alpha n!$. Instead, the numerator is equal to $n!/\binom{n}{\sum_{i=1}^{n}a_{i}}$, as i) $m^{(k_{\alpha})}(x)=m(x_{\alpha})$ for some $x_{\alpha}\in\Pi(x)$, ii) for each $x^{\prime}=(t,d,a^{\prime})\in\Pi(x)$ there are exactly $n!/\binom{n}{\sum_{i=1}^{n}a_{i}}$ permutations $\sigma$ such that $x^{\prime}=(t,d,a_{\sigma})$, and iii) $m(x)$ assumes distinct values on distinct points of $\Pi(x)$, by assumption. Thus, by Equation 8,

[TABLE]

Supporting Information