Sample size calculation for multi-arm parallel design with restricted mean survival time

Abstract

With the recent advances in oncology treatment, restricted mean survival time (RMST) is increasingly being used to replace the routine approach based on hazard ratios in randomized controlled trials for time-to-event outcomes. While RMST has been widely applied in single-arm and two-arm designs, challenges still exist in comparing RMST in multi-arm trials with three or more groups. In particular, it is unclear in the literature how to compare more than one intervention simultaneously or perform multiple testing based on RMST, and sample size determination is a major obstacle to its penetration to practice. In this paper, we propose a novel method of designing multi-arm clinical trials with right-censored survival endpoint based on RMST that can be applied in both phase II/III settings using a global ${\chi }^2$ test as well as a modeling-based multiple comparison procedure. The framework provides a closed-form sample size formula built upon a multi-arm global test and a sample size determination procedure based on multiple-comparison in the phase II dose-finding study. The proposed method enjoys strong robustness and flexibility as it requires less a priori set-up than conventional work, and obtains a smaller sample size while achieving the target power. In the assessment of sample size, we also incorporate practical considerations, including the presence of non-proportional hazards and staggered patient entry. We evaluate the validity of our method through simulation studies under various scenarios. Finally, we demonstrate the accuracy and stability of our method by implementing it in the design of two real clinical trial examples.

1. Introduction

Multi-arm trials with time-to-event outcomes are increasingly common today due to the limited number of patients and the presence of multiple promising treatment regimens. In a review of PubMed published randomized trials in December 2012 by Odutayo et al., ¹ 79% of them were parallel-group designs, of which 21% had more than two groups. These trials typically involve a combination of active treatments, different regimens or doses of an intervention, a placebo, or standard treatment. The primary objective of the trial will dictate the appropriate analysis strategy, such as using a single global test to examine intervention differences or adopting a dose–response model in phase II to detect trends.

In the design stage of a multi-arm trial, an important issue is that the sample size required is more complex than a two-arm trial. The complexity arises because of the limited information available for all arms and the need to account for multiple comparisons. In the past, Makuch and Simon ² and Liu and Dahlberg ³ proposed heuristic approaches to calculate the required sample size for two-group comparisons and then multiplied it by the number of groups. These methods do not give adequate power and have stringent survival curve settings. Later approaches built the sample size formula on the global test of no overall treatment effects. Ahnn and Anderson ⁴ proposed a sample size formula for comparing more than two survival curves using the log-rank test, and they ⁵ generalized it by allowing for flexibility in patient entry, non-compliance and more using the Tarone–Ware class of test statistics. Recent extensions proposed by Halabi and Singh ⁶ and Barthel et al. ⁷ have addressed issues such as unequal allocation ratios and non-proportional hazards (NPHs) in survival curves. Besides parallel design for multi-arm comparison, multi-arm design based on dose-finding studies such as the method proposed by Andrillon et al. ⁸ and Bornkamp et al. ⁹ have also emerged, as well as the sample size calculation procedure developed by Deng et al. ¹⁰

The methods discussed earlier mainly use hazard ratios (HRs) to quantify treatment effects analytically. However, with the recent evolution of oncology studies, where the proportional hazards (PHs) assumption is often violated, there is a growing trend of using supplementary measurements, such as mean survival times, to summarize the treatment benefit. NPHs are prevalent in immunotherapy studies, where the treatment effect may be delayed such as the case in Alexander et al. ¹¹ This delay makes it challenging to draw clinically meaningful conclusions from HRs, which could lead to misleading results due to model mis-specification.^12,13 Furthermore, the HR, averaged over the duration of the trial, does not provide any information about the absolute effect along with accumulated risks or timescale. As a result, it may be difficult for clinicians and patients to comprehend the results. In such cases, restricted mean survival time (RMST), by allowing the treatment effect to be communicated throughout different periods of time in life, can give clinicians a better understanding of treatment benefit and avoid misleading conclusions when NPHs occurs. ¹⁴

In light of the advantages of RMST, there is a need for the development of sample size formulae based on this measure in multi-arm trials. In fact, there is a surge of interest in developing study designs with RMST, such as the non-inferiority trial proposed by Weir and Trinquart , ¹⁵ equality trial by Eaton et al. ¹⁴ and group sequential design by Lu and Tian ¹⁶ and Murray and Tsiatis. ¹⁷ However, there has been little research on sample size calculation for comparing RMST in multi-arm trials. To fill this gap, we propose two main strategies for sample size determination and analysis with different aims. According to the Consolidated Standards of Reporting Trials (CONSORT) extension for reporting multi-arm parallel trials, ¹⁸ two of the main objectives and analysis strategies in multi-arm trials are to compare all groups at once with a global test and to examine the trend of different doses of intervention with dose–response models. The first research objective is commonly seen in phase III trials. The second objective becomes more popular these years in phase II dose-finding studies such as the chronic obstructive pulmonary disease (COPD) example (see Section 5.2), where four once-daily doses were considered along with a placebo and the five dosage levels were treated as five different arms and an optimal dose along with a dose–response model were finally selected. Given these two aims, we enriched the study for multi-arm treatments using RMST by providing innovative designs tailored to these two objectives. In particular, under the first objective of comparing all groups, we proposed a ${\chi }^2$ test for the overall RMST difference and derived a closed-form sample size formula. Under the second objective of finding the optimal dosage level, we extend the classical multiple comparison procedures combined with modeling techniques (MCP-Mod) proposed by Bretz et al. ¹⁹ to the context of RMST measure and develop a new sample size determination procedure in this context. Both strategies offer valuable insights for researchers working with the RMST endpoint in multi-arm trials. Notice that multi-arm designs may have several other features, such as crossover or factorial designs. However, due to the limitations of space in our manuscript, we primarily focused on the parallel design and its related strategies.

In addition to the above-mentioned benefits, another appealing aspect of the proposed RMST-based sample size formulae is their robustness to the misspecification of distribution assumption compared to conventional HR-based sample size formulae. This is because it requires less setting-up a priori before the trial on survival curves and the variance estimator of the RMST ratio is often smaller than that for HR. ²⁰ Furthermore, our approach is also flexible and cost-saving, as it adjusts for piece-wise survival distribution and staggered patient entry in the sample size formula, and it turns out in simulation studies that it achieves the targeted power with a smaller sample size than required by other methods. In addition, RMST benefits from its robust nonparametric inference when compared with the Cox model-based HR. ²¹ The asymptotic normality of RMST and its analytic form of the variance estimator facilitate the calculation of power and sample size. Given the above flexibility and robustness, it can be widely applicable in the practical development of cancertreatments.

The rest of this article is organized as follows. In Section 2, we propose a global test for equality of RMSTs in multi-arm parallel design and the underlying sample size determination method accommodating NPH and staggered patient entry. In Section 3, we provide the MCP-Mod technique built on RMST with its sample size determination procedure. The performance of the proposed methods is demonstrated with simulation results in Section 4 and two real trial examples in Section 5. Section 6 concludes this paper with some discussions.

2. Multi-arm parallel design for global comparison

Multi-arm parallel design with the aim of global comparison is common in phase II and phase III trials. For example, one of the primary objectives of phase II trials with multiple candidate dose levels is to identify the optimal biological dose to be used in phase III trials that produce the most desirable effect. It is appropriate to compare multiple treatment arms altogether using a single global test to detect the efficacy signal, which is critical to the go/no-go decision at the end of phase II studies. In addition, in phase III trials, multiple treatment comparisons are also very common when considering drug combinations in cancer treatments. In this section, we propose a sample size calculation procedure built upon RMST under the global test of examining whether the intervention effects differ instead of how they differ.

2.1. Global test

Suppose in a clinical trial the primary endpoint is a disease-related event, such as death or disease progression, and the variable of interest is the time to event, namely $T$ . Consider a trial setting with $K$ parallel study arms labeled $k=1,\dots ,K$ , in which $N$ individual subjects are randomized into one of them with $n_k(k=1,\dots ,K)$ individuals in the $k$ th arm. We let $k=1$ represent the control group. The aim is to compare the RMST of the $K$ arms globally at a specified restriction time $\tau$ . Let $R_k(\tau )$ be the RMST in group $k$ ( $k=1,\dots ,K$ ). It is of interest to test

$$H_0:R_1(\tau )=\cdots =R_K(\tau )$$

(1)

against

$$H_1:R_k(\tau )\ne R_l(\tau )\text{ for at least one }k\ne l(1\le k,l\le K)$$

(2)

Further, let $\mathit{R}(\tau )=(R_1(\tau ),\dots ,R_K(\tau ){)}^{\prime }$ and assume for simplicity $\mathit{R}\equiv \mathit{R}(\tau )$ .

In the cases with right-censored data, RMST is often estimated by the area under the Kaplan–Meier curve $\hat{R}(\tau )={\int }_0^{\tau }\hat{S}(t)dt$ , where $\hat{S}(t)$ is the Kaplan–Meier estimator of the survival function $S(t)=P(T>t)$ . Meier ²² proved that the asymptotic mean and variance of $\hat{R}(\tau )$ are ${\int }_0^{\tau }S(t)dt$ and ${\sigma }^2(\tau )=2{\int }_0^{\tau }{\int }_u^{\tau }S(u)S(v)C(u)dvdu$ , respectively, where $C(u)$ is the sum of squared coefficient of variation for $S(u)$ . Further, Klein and Moeschberger ²³ provided the Greenwood plug-in variance estimator in the following form:

$$V_k=\widehat{Var}({\hat{R}}_k(\tau ))=\sum _{j=1}^{n_k}{[{\int }_{t_j}^{\tau }{\hat{S}}_k(t)dt]}^2\frac{{\delta }_{kj}}{Y_{kj}(Y_{kj}-{\delta }_{kj})}$$

(3)

where ${\delta }_{kj},Y_{kj},k=1,\dots ,K,j=1,\dots ,n_k$ are the censoring indicators and number of at-risk samples at the observed time of $j$ th sample in group $k$ , and $D_k=\sum _{j=1}^{n_k}{\delta }_{kj}$ is the total number of observed events in group $k$ . Tian et al. ²⁴ demonstrated that under a mild condition on the censoring distribution, the asymptotic normality for ${\hat{R}}_k(\tau )$ is valid up to the time point $\tau$ , where $\tau$ can be less than or equal to the largest follow-up time in the study. Following their condition, we assume a sufficient condition on the censoring distribution that

$$\underset{t\to {\tau }_B}{lim}\frac{g_B(t)}{({\tau }_B-t{)}^{1+ϵ}}>0$$

for an $ϵ\in (0,1),$ where $B$ is the censoring time, $g_B(t)$ is the probability density function of $B$ and ${\tau }_B=inf\{\tau |P(B\ge \tau )=0\}$ . This condition is satisfied as long as the censoring distribution does not approach zero too fast. Given the sample data, ${\hat{\tau }}_B$ , the finite sample approximation of ${\tau }_B$ can be defined as the largest follow-up time (either censored or uncensored). If we assume a staggered patient entry under piece-wise uniform distribution with moderate rates, as in most clinical trials, this condition is still satisfied because censoring does not become sparse even close to ${\hat{\tau }}_B$ under the flat and moderate rates of patient accrual. Hence, the test proposed in this article is based on the asymptotic normality of ${\hat{R}}_k(\tau )$ for $\tau \in [0,{\hat{\tau }}_B]$ .

Notice that the global hypothesis in (1) can also be expressed as a form of linear combinations of RMST. We let

$${\mathit{C}}_{K-1,K}=\left(\begin{array}{ccccc}-1& 1& 0& \cdots & 0\\ -1& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -1& 0& 0& \cdots & 1\end{array}\right)$$

(4)

be the coefficient matrix with the $k$ th row being a vector of coefficients for the $k$ th linear combination of $\mathit{R}$ , the global test is equivalent to testing $H_0:\mathit{C}\mathit{R}=\mathbf{\text{0}}$ against $H_1:\mathit{C}\mathit{R}\ne \mathbf{\text{0}}$ . Similar to Wald’s statistic testing uncensored means, we propose a ${\chi }^2$ statistic for the RMST global test

$$Q=(\mathit{C}\hat{\mathit{R}}{)}^{\prime }(\mathit{C}\hat{\mathbf{\Sigma }}{\mathit{C}}^{\prime }{)}^{-1}(\mathit{C}\hat{\mathit{R}})$$

(5)

where $\hat{\mathit{R}}$ is the estimated RMST until $\tau$ and $\hat{\mathbf{\Sigma }}=\mathrm{diag}(V_1,\dots ,V_K)$ is the estimated variance covariance matrix of $\hat{\mathit{R}}$ , where $V_k=\widehat{\mathrm{Var}}({\hat{R}}_k),k=1,\dots ,K$ is the Greenwood plug-in variance estimator of RMST of the $k$ th arm as specified in (3). Due to the asymptotic properties that $\hat{\mathit{R}}$ follows a multivariate normal distribution $N(\mathit{R},\mathbf{\Sigma })$ , where $\mathbf{\Sigma }$ can be consistently estimated by $\hat{\mathbf{\Sigma }}$ , $Q$ has an asymptotic ${\chi }^2(K-1)$ distribution under $H_0$ .

2.2. Sample size calculations for the global test

In the study design stage, we propose to set a priori specification of restricted time $\tau$ in the protocol guided by clinical, practical, and financial considerations as no data is obtained before the trial. In addition, deriving the sample size formula for survival endpoints in multi-arm trials is a complex task that typically requires more specification of the survival curves and their corresponding variances for each treatment group. However, obtaining complete information about survival distributions is often challenging for investigators. To address this issue, we propose a new approach that only requires specification of the survival distribution of the treatment group with the maximum treatment effect, which we call the representative treatment group. For other treatment arms, we assume the same estimated variance of RMST as that in the representative treatment arm. By avoiding the need to specify all treatment survival distributions, this technique offers greater flexibility and reduces the risk of mis-specifying survival curves, making it a valuable tool in real-world situations. We will later confirm through empirical study that this assumption is reasonable, and even if violated, it does not significantly affect the power.

Consider a multi-arm clinical trial with $K$ arms, where the first arm represents the control group. We assume that the variance in each treatment group is equal and denoted as $\mathrm{Var}({\hat{R}}_2)=\cdots =\mathrm{Var}({\hat{R}}_K)={\sigma }^2(\tau )$ . However, the control arm may have a different variance, denoted as ${\sigma }_1^2(\tau )$ , and we define a constant $c={\sigma }^2(\tau )/{\sigma }_1^2(\tau )$ to account for this difference. To simplify the study design, we assume an equal allocation of $n$ subjects in each group. Under these assumptions, we can express the difference in RMST between each treatment group and the control group as ${\mathrm{\Delta }}_k=R_k(\tau )-R_1(\tau )$ , where $k=2,\dots ,K$ . The hypothesis in (1) can be formulated as $H_0:{\mathrm{\Delta }}_k=0$ for $k=2,\dots ,K$ . With this notation, the test statistic $Q$ can be simplified as follows:

$$Q=\frac{\sum _{k=2}^K{\hat{\mathrm{\Delta }}}_k^2}{V}-\frac{{\left(\sum _{k=2}^K{\hat{\mathrm{\Delta }}}_k\right)}^2}{(K-1+c)V}$$

(6)

where ${\hat{\mathrm{\Delta }}}_k={\hat{R}}_k(\tau )-{\hat{R}}_1(\tau )$ and $V={\hat{\sigma }}^2(\tau ).$ A detailed derivation can be found in the Appendix. It is important to note that the simplification discussed here is only intended to facilitate the sample size calculation procedure by linking the required sample size to the desired power through statistical information. The exact test still takes the form of (5). Under the alternative hypothesis $H_1$ , where ${\mathrm{\Delta }}_k\ne 0$ for at least one $k=2,\dots ,K$ , we have $\mathit{C}\hat{\mathit{R}}\sim N(\mathbf{\Delta },\mathit{C}\hat{\mathbf{\Sigma }}{\mathit{C}}^{\prime })$ . As a result, the test statistic $Q$ follows a non-central ${\chi }^2(K-1,\delta )$ distribution, where $K-1$ represents the degrees of freedom and $\delta$ is the non-centrality parameter.

On the one hand, the non-centrality parameter $\delta$ needs to be known to achieve a given power $1-\beta$ and significance level $\alpha$ that

$$\delta =inf\{\delta :1-\beta <P(Q>{\chi }_{1-\alpha }^2(K-1)|Q\sim {\chi }^2(K-1,\delta ))\}$$

(7)

where ${\chi }_{1-\alpha }^2(K-1)$ is the upper $\alpha$ quantile of ${\chi }^2(K-1)$ .

On the other hand, the non-centrality parameter $\delta$ can be derived from the known ${\mathrm{\Delta }}_k$ ’s and $V$ . By equating the first moment $E(Q|H_1)$ to $\delta +K-1$ , that is, the expectation of ${\chi }^2(K-1,\delta )$ , we have

$$\delta =\frac{\sum _{k=2}^K{\mathrm{\Delta }}_k^2}{{\sigma }^2(\tau )}-\frac{{\left(\sum _{k=2}^K{\mathrm{\Delta }}_k\right)}^2}{(K-1+c){\sigma }^2(\tau )}$$

(8)

A detailed derivation of $E(Q|H_1)$ can be found in the Appendix.

Then the required sample size can be solved by equating (7) and (8) after replacing ${\sigma }^2(\tau )$ by its asymptotic value. Specifically, Pepe and Fleming ²⁵ showed that given the known RMST margin $\mathbf{\Delta }$ and setting the survival distribution $S(t),S_1(t)$ , hazard function $\lambda (t),{\lambda }_1(t)$ , cumulative loss to follow-up distribution $G(t),G_1(t)$ of the representative treatment group and control group, respectively, as well as the restriction time $\tau$ , asymptotic variance in the representative treatment group can be expressed as

$$V={\int }_0^{\tau }\frac{{({\int }_t^{\tau }S(x)dx)}^2\lambda (t)}{nS(t)(1-G(t))}dt$$

(9)

where $n=n_2=\cdots =n_K$ . $V_1={\hat{\sigma }}_1^2(\tau )$ can be derived in a similar manner. Then, for a simple case of equal allocation of sample size, the sample size required in each arm is given by

$$n=\frac{\delta {\int }_0^{\tau }\left[\frac{{({\int }_t^{\tau }S(x)dx)}^2\lambda (t)}{S(t)(1-G(t))}\right]dt}{\sum _{k=2}^K{\mathrm{\Delta }}_k^2-{\left(\sum _{k=2}^K{\mathrm{\Delta }}_k\right)}^2(K-1+\hat{c})}$$

(10)

where

$$\hat{c}=V/V_1={\int }_0^{\tau }\frac{{({\int }_t^{\tau }S(x)dx)}^2\lambda (t)}{S(t)(1-G(t))}dt/{\int }_0^{\tau }\frac{{({\int }_t^{\tau }{S}_1(x)dx)}^2{\lambda }_1(t)}{S_1(t)(1-G_1(t))}dt$$

It is worth noting that the sample size calculation procedure presented in (10) assumes a simple case of equal allocation of sample size. However, in practice, it is often the case that more patients are assigned to the new treatment than those in the control group to learn more about the new drug. In this case, an unbalanced allocation ratio between treatment and control is used, denoted by $r=n_2/n_1=\cdots =n_K/n_1\ne 1$ , where $n_1$ is the sample size for the control group. The sample size required for the control arm can be obtained by modifying equation (10) as follows:

$$n_1=\frac{\delta {\int }_0^{\tau }\left[\frac{{({\int }_t^{\tau }S(x)dx)}^2\lambda (t)}{S(t)(1-G(t))}\right]dt}{\sum _{k=2}^K{\mathrm{\Delta }}_k^2-{\left(\sum _{k=2}^K{\mathrm{\Delta }}_k\right)}^2(K-1+\hat{c}r)}$$

(11)

Note that the above sample size formulae can also accommodate NPH by flexibly selecting the survival distribution such as a piece-wise exponential distribution.

Another concern is that in reality, patients are often not all recruited at the beginning of the trial, but enter the trial sequentially over an accrual period, and then are followed up for another fixed period with no new patients entering. At the end of the trial, administrative censoring occurs with some of them not experiencing the events yet. This is known as staggered patient entry, and it can be accommodated in the sample size calculation framework with a minor modification. Define a cumulative distribution function of the patient recruitment time $F^R(t)$ such as the piece-wise uniform distribution within an accrual period from calendar time 0 to $t_1$ , and the planned follow-up time $t_2$ , where $t_1<t_2$ , administrative censoring is brought into the trial. Let $s$ be the calendar time when the patient is accrued, due to the administrative censoring, the potential exposure time since patient entry becomes $T=min\{u,t_2-s\},$ where $u$ denotes the patient survival time. Then the survival function, hazard function, and censoring distribution adjusting for staggered patient entry are, respectively,

$$\begin{array}{rl}{S}^E(t)& =S(t)F^R(t_2-t),{\lambda }^E(t)=\lambda (t)+f^R(t_2-t)/F^R(t_2-t)\\ G^E(t)& =1-(1-G(t))F^R(t_2-t)\end{array}$$

Hence the required sample size in the control group under a certain allocation ratio $r$ becomes

$$n_1=\frac{\delta {\int }_0^{\tau }\left[\frac{{({\int }_t^{\tau }{S}^E(x)dx)}^2{\lambda }^E(t)}{S^E(t)(1-G^E(t))}\right]dt}{\sum _{k=2}^K{\mathrm{\Delta }}_k^2-{\left(\sum _{k=2}^K{\mathrm{\Delta }}_k\right)}^2/(K-1+\hat{c}r)}$$

(12)

3. Multi-arm dose-finding studies with modelling

Parallel design aims to detect an overall difference in intervention effects. In phase-II dose-finding studies where the main objective is to examine a trend of the dosage effects, the proposed global test is not sufficient and a dose–response model should be considered in the test. In this section, we propose a more sophisticated hybrid approach for multi-arm dose-finding study with RMST motivated by MCP-Mod ¹⁹ as well as the sample size calculation. The test still utilizes the linear combinations of RMST but in addition, provides flexibility in modeling the relationship between dosage levels and RMST. We will start this section by introducing MCP-Mod first.

3.1. Multiple comparison procedures combined with modeling techniques (MCP-Mod)

MCP-Mod is a unified approach in dose-finding study, which innovatively combines the proof-of-concept (PoC) with dose-finding into a single study. By combining multiple comparison procedures with modeling techniques into one step, this method leads to a seamless design. In their framework, the response $Y$ is an uncensored continuous variable observed for a given set of groups of patients corresponding to an increasing dose sequence $d_1<\cdots <d_K$ , where $d_1$ represents the dose in the placebo group. For the purpose of detecting a global dose–response trend, MCP-Mod assumes a linear model of the response

$$Y_{ij}={\mu }_i+{ϵ}_{ij},{ϵ}_{ij}\stackrel{ind}{\sim }N(0,{\sigma }^2),i=1,\dots ,K,j=1,\dots ,n_i$$

where ${\mu }_i=f(d_i;\mathit{\theta })$ is the mean response at dose $d_i$ and $f(\cdot )$ is the candidate dose–response model with parameter $\mathit{\theta }$ , while $i$ refers to the dose group and $j$ refers to the patient within group $i$ . Historical data and clinical experts’ opinions will provide initial knowledge or guess to determine the parameters of the dose–response model.

To achieve the goal of model selection, a set of $M$ parameterized candidate models $f_1,\dots ,f_M$ is assumed, corresponding to the mean response vectors for model $m$ under the $K$ doses ${\mu }_m=({\mu }_{m1},\dots ,{\mu }_{mK})$ , $m=1,\dots ,M$ . Under these models, contrast coefficient ${\mathit{c}}_m^{\prime }=(c_{m1},\dots ,c_{mK})$ subject to ${\mathit{c}}_m^{\prime }\mathbf{1}=0$ and $\sum _{i=1}^K{c}_{mi}^2=1$ is further determined for the purpose of constructing the test of PoC for each model (see Section 3.1** in Bretz et al. ¹⁹ ). Then each of the dose–response models is tested by a single contrast test with $H_0^m:{\mathit{c}}_m^{\prime }\mathit{\mu }=0$ against $H_1^m:{\mathit{c}}_m^{\prime }\mathit{\mu }>0$ . This leads to the contrast test statistic defined as

$$T_m=\frac{\sum _{i=1}^K{c}_{mi}{\overline{Y}}_i}{S\sqrt{\sum _{i=1}^K{c}_{mi}^2/n_i}},m=1,\dots ,M$$

where $S^2=\sum _{i=1}^K\sum _{j=1}^{n_i}(Y_{ij}-{\overline{Y}}_i{)}^2/(N-K)$ is the pooled variance and $N=\sum _{i=1}^K{n}_i$ . Note that although the true model is unknown, MCP-Mod can control the familywise error rate (FWER) when establishing PoC and selecting the best model out of the $M$ candidates. To achieve a single decision of the existence of the trend, that is, at least one $H_0^m,m=1,\dots ,M$ is rejected, the final test statistic for verifying PoC takes $T_{max}=\underset{m}{max}\{T_m\}$ . PoC is established if $T_{max}\ge q_{1-\alpha }$ , where $q_{1-\alpha }$ is a multiplicity adjusted critical value obtained by the equicoordinate upper $\alpha$ quantile of the joint distribution of the vector of contrast statistics ${\mathit{T}}^{\prime }=(T_1,\dots ,T_M)$ under $H_0$ in (1). It follows that ${\mathit{T}}^{\prime }\sim MV{T}_M(N-K;\mathbf{0},\mathcal{I})$ , a multivariate $t$ distribution of dimension $M$ with $N-K$ degrees of freedom and correlation matrix $\mathcal{I}=({\rho }_{ij})$ , where

$${\rho }_{ij}=(\sum _{l=1}^K{c}_{il}{c}_{jl}/n_l)/{(\sum _{l=1}^K{c}_{il}^2/n_l\sum _{l=1}^K{c}_{jl}^2/n_l)}^{1/2},1\le i,j\le K$$

Once PoC is established, the best model among all significant candidate models is selected for dose estimation in the last stage of MCP-Mod. A $P$ -value method or other criteria Akaike information criterion or Bayesian information criterion can be adopted to select the best model. After that, the target dose of interest can be estimated by the selected dose–response model, but this is beyond the scope of this paper.

3.2. MCP-Mod with endpoint measured by RMST

Now, we consider dose-finding clinical trials where the response of interest can be right-censored such as time-to-toxicity in oncology trial designs with an aim of detecting an overall trend as well as model selection, and completing the MCP-Mod for a right-censored response from the perspective of RMST. In this case, although the survival outcome $T$ no longer follows a normal distribution, we can still apply the candidate models on RMST thanks to the asymptotic properties of the RMST estimator. Suppose that patients are followed up for an observation window until time $\tau$ . We concentrate on the PoC stage to test the difference of RMST among dosage levels and model selection, while the followed dose selection procedure can be done naturally by any covariate-adjustment RMST method such as the method proposed by Tian et al., ²⁶ which we leave to the readers for further exploration.

Assume that with some clinical knowledge, a candidate model set $\mathcal{M}=\{M_m:m=1,\dots ,M\}$ and the model parameters at a specified time point $\tau$ are given, that is, for each dose $d_1,\dots ,d_K$ planned, the model generates ${\mathit{R}}_m=(R_{m1},\dots ,R_{mK})$ . ${\mathit{R}}_m$ can be viewed as ${\mathit{\mu }}_m$ in Section 3.1, and $R_{mi}=f_m(d_i,{\mathit{\theta }}_m)$ . Frequently used potential candidate models $f_m(\cdot )$ are listed non-exhaustively in Table 1. Further assuming the pre-determined known contrast coefficients ${\mathit{c}}_m^{\prime }=(c_{m1},\dots ,c_{mK})$ has been computed (Section 3.2.1) under the constraints ${\mathit{c}}_m^{\prime }\mathbf{1}=0$ and $\sum _{i=1}^K{c}_{mi}^2=1$ , we test on the null hypothesis $H_0^m:{\mathit{c}}_m^{\prime }\mathit{R}=0$ . Without loss of generality, we suppose that if a dose–response relationship exists, dosages higher than $d_1$ always result in a larger RMST. Then the alternative hypothesis corresponds to $H_1^m:{\mathit{c}}_m^{\prime }\mathit{R}>0,m=1,\dots ,K,$ while in the other way around a left-tailed test is adopted. Either way, the test statistic becomes

$$Z_m=\frac{\sum _{i=1}^K{c}_{mi}{\hat{R}}_i}{\sqrt{\sum _{i=1}^K{c}_{mi}^2{V}_i}}$$

(13)

where ${\hat{R}}_i$ is the estimated RMST and $V_i$ is the same Greenwood plug-in estimator in group $i$ as defined in equation (3). In this way, the unified ${\chi }^2$ statistic proposed in equation (5) is split into $M$ $Z$ -statistics to achieve the goal of model selection. Due to the asymptotic normality of ${\hat{R}}_i,i=1,\dots ,K$ , $Z_m$ follows $N(0,1)$ under $H_0^m$ and jointly ${\mathit{Z}}^{\prime }=(Z_1,\dots ,Z_M)$ follows $MVN(0,\mathcal{I})$ , where $\mathcal{I}=({\rho }_{ij})$ is the variance–covariance matrix with

$${\rho }_{ij}=Corr(Z_i,Z_j)=\frac{\sum _{l=1}^K{c}_{il}{c}_{jl}}{\sqrt{\sum _{l=1}^K{c}_{il}^2\sum _{l=1}^K{c}_{jl}^2}}$$

(14)

A detailed derivation is shown in the Appendix. Then, the critical value $q_{1-\alpha }$ can be derived by computing the equicoordinate quantile function of this distribution, or by solving

$$P(max\{Z_1,\dots ,Z_M\}\ge q_{1-\alpha }|H_0)=\alpha$$

(15)

In order to combine them into a single decision of PoC, we take $Z_{max}=max\{Z_1\cdots ,Z_M\}$ such that for a critical value $q_{1-\alpha }$ , a significant dose–response signal is established if $Z_{max}>q_{1-\alpha }$ .

Table 1.

Typical dose–response models used in candidate set.

Model	$f(d,\mathit{\theta })$
Linear	$E_0+\gamma d$
Exponential	$E_0+E_1\exp \{d/\gamma \}$
Logistic	$E_0+E_{max}/(1+\exp \{(E{D}_{50}-d)/\gamma \})$
Quadratic	$E_0+{\beta }_{1}d+{\beta }_2{d}^2$
Emax	$E_0+E_{max}d/(E{D}_{50}+d)$
Linear log-dose	$E_0+\gamma \log (d+1)$

Once at least one individual test is significant, PoC is established with a reference model set leading to significant contrast tests. Note that FWER is also strongly controlled at $\alpha$ for each single contrast test. In the second step after PoC is verified, we go through the $M$ tests one by one and decide upon the best model among those in the reference set. Good candidate models can be evaluated by simple criteria such as the minimum $P$ -value or absolute predicted error. ²⁶ It is shown in the simulation that under a sufficient sample size, the model can be selected accurately even if the candidate model parameters do not fit the true ones well.

An additional advantage of the proposed method is its emphasis on modeling. Once the multiple comparison procedure is completed with a selected optimized model, the final step in MCP-Mod consists of fitting the dose–response model. We propose to estimate the relationship between RMST and dosage levels by nonlinear least squares estimates given in Bates and Watts. ²⁷ Since the optimization relies heavily on the initial parameter estimates and sometimes the data is not enough to fit a model due to limited dosage levels, we propose to do data augmentation before modeling to increase the amount of data. Specifically, we add simulated copies of RMSTs at each dosage level $k\in \{1,\dots ,K\}$ from $N({\hat{R}}_k,V_k)$ to the original ${\hat{R}}_k$ . Then the estimated model can be used in dose selection, for example, the selection of the smallest dose that shows a clinically relevant and statistically significant effect. ²⁸

3.2.1. Optimum contrast coefficients

The implementation of the above method starts by assuming the parametric candidate models, doses, and sample size, followed by computing the optimum contrast coefficients. Candidate models fitting the true curves well, together with the contrast coefficients will boost the power of the test and reduce the required sample size. In this subsection, we discuss the determination of optimum contrast coefficients given a specified model $f_m(d;\mathit{\theta }),m\in \{1,\dots ,M\}$ . Under the same sense from Bretz et al., ¹⁹ we aim to maximize the probability of rejecting the null hypothesis when the model is correct.

Let $g(\mathit{c},{\mathit{R}}_m)$ be the objective function of rejecting $H_0$ when the model $f_m(d,\mathit{\theta })$ is correct, where ${\mathit{R}}_m=(R_{m1},\dots ,R_{mK})$ is the RMST vector under model $m$ and $\mathit{c}$ is the vector of contrast coefficients. Denote ${\mathit{c}}_{\mathrm{opt}}(f_m)$ the vector of optimum contrast coefficient under model $m\in \{1,\dots ,M\}$ with constraints of ${\mathit{c}}_{\mathrm{opt}}^{\prime }\mathbf{1}=0$ and $\sum _{i=1}^K{c}_{\mathrm{opt},i}^2=1$ . It can be shown that maximizing the chance of rejecting $H_0$ when the model is correct it is equivalent to

$$\begin{array}{rl}& {\mathit{c}}_{\mathrm{opt}}(f_m)=\underset{\mathit{c}}{\mathrm{argmax}}g(\mathit{c},{\mathit{R}}_m)\\ & g(\mathit{c},{\mathit{R}}_m)=\frac{{\mathit{c}}^{\prime }{\mathit{R}}_m}{\sqrt{{\mathit{c}}^{\prime 2}{\mathit{\sigma }}_m^2}}\end{array}$$

(16)

where the theoretical variance of RMST under model $m$ is ${\mathit{\sigma }}_m^2=({\sigma }_{m1}^2,\dots ,{\sigma }_{mK}^2{)}^{\prime }$ is derived from (9) by assuming exponential survival $\mathrm{Exp}({\lambda }_{mi})$ for group $i$ and a same loss to follow-up distribution. Specifically, we obtain ${\lambda }_{mi}$ by solving $R_{mi}={\int }_0^{\tau }\exp (-{\lambda }_{mi}t)dt$ . Notice that ${\mathit{c}}_{\mathrm{opt}}(f_m)$ cannot be expressed in closed form so we adopt the technique of spherical parameterization ²⁹ to reparameterize the coefficients. The detailed optimization procedures are presented in the Appendix. In this way, constrained optimization translates into a more robust and easier unconstrained one, and standard numerical optimization is available in any software to obtain ${\mathit{c}}_{\mathrm{opt}}(f_m)$ .

3.3. Power and sample size calculations for MCP-Mod of RMST

In practice, it is a crucial step at the study design stage to determine the sample size required for certain purposes. For dose-finding studies with time-to-event response variable, we aim to determine the sample size required to ensure a minimum power level $1-\beta$ for detecting the dose–response relationship under a certain dose–response model with RMST vector ${\mathit{R}}_m=(R_{m1},\dots ,R_{mK})$ . For simplicity, we restrict ourselves to the equal allocation of sample sizes in each dose group, while the general cases of unbalanced design are also straightforward. Assuming the candidate set of models $\mathcal{M}$ , the common cumulative loss to follow-up distribution $G(t)$ , hazard and survival functions ${\lambda }_{mi}(t),S_{mi}(t),m=1,\dots ,M,i=1,\dots ,K$ (e.g. the exponential survival distribution) that yields the corresponding $R_{mi}=f_m(d_i;\mathit{\theta })$ , restriction time $\tau$ , type-I error $\alpha$ , and power $1-\beta$ , under the above MCP-Mod method for survival endpoint, we first consider the power calculation under model $m\in \mathcal{M}$ . It follows from the distribution of the test statistic $Z_1,\dots ,Z_M$ in Section 3.2 that the power to detect PoC is

$$\begin{array}{rl}{\pi }_m(n)& =P(max\{Z_1,\dots ,Z_M\}\ge q_{1-\alpha }|\mathit{R}={\mathit{R}}_m)\\ & =1-P(Z_1<q_{1-\alpha },\dots ,Z_M<q_{1-\alpha }|\mathit{R}={\mathit{R}}_m)\end{array}$$

(17)

where $q_{1-\alpha }$ is the critical value. ${\pi }_m(n)$ can be computed by the “pmvnorm” function in R of a multivariate normal distribution with correlation matrix $\mathcal{I}$ and mean vector

$$(\frac{\sum _{i=1}^K{c}_{1i}{R}_{1i}}{\sqrt{\sum _{i=1}^K{c}_{1i}^2{\hat{\sigma }}_{1i}^2}},\dots ,\frac{\sum _{i=1}^K{c}_{Mi}{R}_{Mi}}{\sqrt{\sum _{i=1}^K{c}_{Mi}^2{\hat{\sigma }}_{Mi}^2}})$$

where

$${\hat{\sigma }}_{mi}^2={\int }_0^{\tau }\frac{{({\int }_t^{\tau }{S}_{mi}(x)dx)}^2{\lambda }_i(t)}{n_{mi}{S}_{mi}(t)(1-G(t))}dt$$

is the theoretical variance under model $m$ ³⁰ in terms of the design parameters. Then under certain sample size $n=n_1+\cdots +n_K$ , a power vector $\mathit{\pi }(n)=({\pi }_1(n),\dots ,{\pi }_M(n))$ can be obtained without any knowledge of the true model under $H_1$ . To synthesize the power vector into a single value, we can adopt different weight functions or any summary function, for example, we can simply use a conservative approach to condense it into ${\pi }^{*}(n)=min\{{\pi }_1(n),\dots ,{\pi }_M(n)\}$ as the unified power measure. Then the sample size can be determined by searching for the smallest integer $n$ from an initial value until the unified power reaches the target $1-\beta$ .

In the next section, we will investigate how the restriction time $\tau$ and RMST margins between each treatment and control group affect the sample size required. Basically, the required sample size relies mainly on the combination of effects of candidate models and parameters.

4. Simulation studies

We conducted simulation studies to investigate the performance of the proposed sample size determination methods for the global test and dose-finding studies. Specifically, we aimed to answer the following questions: Does the proposed sample size formula for global testing provide accurate estimates under PH and NPH scenarios? Secondly, how well do the sample size estimator and model selection procedure with MCP-Mod perform in dose-finding studies? Thirdly, how do factors such as restriction time, number of treatment arms, RMST margins, and candidate model settings in MCP-Mod affect the performance of sample size determination? Lastly, how does the performance differ between the two proposed sample size determination procedures and conventional methods in multi-arm designs? To answer these questions, we look into global tests for PH scenarios and NPH scenarios, as well as the dose-finding studies separately in the following subsections.

In all sets of simulations, we assumed patients would be accrued at a constant rate over two consecutive periods before $t_1=12$ months and would be followed up until $t_2=36$ months. Specifically, patient entry followed a piece-wise uniform distribution with a rate of 1 out of 16 before 8 months and then a rate of 1 out of 8 until 12 months. Independent right-censored time follows an exponential distribution with an annual dropout rate of 5%. Other design specifications included a 9-month median survival time for the control group with exponential survival, a goal of attaining 90% power, and a type-I error of $\alpha =0.05$ . We present all simulation results of the required sample size and empirical power based on 1000 simulated trials in Tables 2 to 4.

Table 2.

Simulation results for PHs in the global test: $n$ is the sample size required to achieve 90% power in each group, Power $\mathrm{\%}$ is the empirical power under sample size $n$ based on 1000 simulations. ${\mathrm{\Delta }}_2,{\mathrm{\Delta }}_3$ , and ${\mathrm{\Delta }}_4$ are the RMST differences between the treatment group and the control group, and numbers in the parentheses are the corresponding HRs.

	K = 3						K = 4
	Parameters		RMST		Barthel		Parameters			RMST		Barthel
$\tau$	${\mathrm{\Delta }}_2$	${\mathrm{\Delta }}_3$	$n$	Power%	$n$	Power%	${\mathrm{\Delta }}_2$	${\mathrm{\Delta }}_3$	${\mathrm{\Delta }}_4$	$n$	Power%	$n$	Power%
24	2.04 (.75)	2.04 (.75)	320	90.4	345	91.4	2.04 (.75)	2.04 (.75)	2.04 (.75)	315	90.8	339	90.3
	2.04 (.75)	0 (1)	336	92.7	363	93.0	2.04 (.75)	2.04 (.75)	0 (1)	247	92.2	264	92.7
	2.04 (.75)	0.75 (.9)	427	89.6	467	93.4	2.04 (.75)	0 (1)	0 (1)	337	90.0	366	91.3
	0.75 (.9)	0.75 (.9)	2303	90.0	2488	92.8	2.04 (.75)	0.75 (.9)	0 (1)	375	88.6	405	92.1
	3.56 (.6)	0 (1)	113	90.7	126	93.3	3.56 (.6)	2.04 (.75)	0.75 (.9)	143	90.6	158	93.3
	3.56 (.6)	0.75 (.9)	133	89.4	150	93.9	3.56 (.6)	3.56 (.6)	3.56 (.6)	105	89.3	111	90.7
	0 (1)	0 (1)	200	5.2	—	—	0 (1)	0 (1)	0 (1)	200	5.6	—	—-
18	1.45 (.75)	1.45 (.75)	378	89.9	449	95.7	1.45 (.75)	1.45 (.75)	1.45 (.75)	379	90.9	439	94.7
	1.45 (.75)	0 (1)	382	90.5	478	95.4	1.45 (.75)	1.45 (.75)	0 (1)	277	89.8	346	94.4
	1.45 (.75)	0.55 (.9)	495	90.3	612	95.0	1.45 (.75)	0 (1)	0 (1)	381	91.3	483	96.3
	0.55 (.9)	0.55 (.9)	2656	90.3	3212	95.2	2.04 (.75)	0.75 (.9)	0 (1)	370	89.4	533	95.0
	2.48 (.6)	0 (1)	139	92.0	169	96.0	2.48 (.6)	1.45 (.75)	0.55 (.9)	179	91.8	210	96.4
	2.48 (.6)	0.55 (.9)	165	91.1	200	95.8	2.48 (.6)	2.48 (.6)	2.48 (.6)	134	91.8	145	92.4
	0 (1)	0 (1)	200	4.7	–	–	0 (1)	0 (1)	0 (1)	200	5.2	–	–

PH: proportional hazard; RMST: restricted mean survival time; HR: hazard ratio.

Table 3.

Simulation results for NPHs in global test at $\tau =24$ : Power $\mathrm{\%}$ is the empirical power with sample size $n$ based on 1000 simulations. ${\mathrm{HR}}_1$ and ${\mathrm{HR}}_2$ refer to the hazard ratio in treatment Groups 2 and 3, respectively, with piece-wise constants in different periods, and the corresponding RMST difference at 24 months is listed in the column of ${\mathrm{\Delta }}_2$ and ${\mathrm{\Delta }}_3$ . $n$ : sample size required in each group when the design model is true under NPH and when the design model is wrong under PH.

	Parameters						Adjusted for NPH				PH
	${\mathrm{HR}}_1$			${\mathrm{HR}}_2$			RMST		Barthel		RMST		Barthel
Scenario			${\mathrm{\Delta }}_2$			${\mathrm{\Delta }}_3$	$n$	Power%	$n$	Power%	$n$	Power%	$n$	Power%
Early effect	$<12$	$\ge 12$		$<12$	$\ge 12$
	0.75	1	1.63	0.75	1	1.63	488	91.4	510	92.8	320	74.7	345	76.4
	0.75	1	1.63	1	1	0	494	90.9	538	94.7	336	75.6	363	78.2
Late effect	$<3$	$\ge 3$		$<3$	$\ge 3$
	1	0.75	1.40	1	0.75	1.40	695	88.8	728	89.7	320	57.4	345	58.9
	1	0.75	1.40	1	1	0	744	91.9	739	91.2	336	55.8	363	60.7
	1	0.71	1.63	1	0.71	1.63	501	90.0	524	91.1	228	54.0	246	61.0
	1	0.71	1.63	1	1	0	541	92.2	535	91.9	241	58.2	262	59.9
Middle effect	3–18	$<3$ , $>18$		3–18	$<3$ , $>18$
	0.6	1	2.23	0.6	1	2.23	275	91.1	275	91.1	106	49.1	114	49.6
	0.6	1	2.23	0.75	1	1.31	370	90.7	368	89.4	143	51.5	157	52.9
	0.695	1	1.63	0.695	1	1.63	508	89.7	509	89.4	203	53.3	219	54.3
	0.695	1	1.63	1	1	0	543	91.1	521	88.3	214	51.2	234	56.5
Crossing survival	$<3$	$\ge 3$		$<3$	$\ge 3$
	1.1	0.7	1.46	1.1	0.8	0.85	870	90.7	885	90.9	277	43.7	302	43.4
	1.2	0.5	2.59	1.2	0.5	2.59	213	89.0	223	91.2	60	36.9	65	42.9
	1.1	0.674	1.63	1.1	0.674	1.63	521	90.1	540	90.8	174	43.1	187	45.6
	1.1	0.674	1.63	1	1	0	570	91.2	544	87.9	184	42.7	201	46.0

NPH: non-proportional hazard; RMST: restricted mean survival time; HR: hazard ratio.

Table 4.

Simulation results for dose-finding study using multiple comparison procedures combined with modeling techniques (MCP-Mod) and the global test. $n$ : sample size per treatment arm. Selected models: the model that is selected most time in 1000 simulations.

		Parameters					Global test		MCP-Mod
Case	$\tau$	${\mathrm{\Delta }}_1$	${\mathrm{\Delta }}_2$	${\mathrm{\Delta }}_3$	${\mathrm{\Delta }}_4$	${\mathrm{\Delta }}_5$	$n$	Power%	$n$	Power%	Selected models	s.e.
a	24	0.75	1.42	2.14	2.72	3.56	142	90.4	137	99.3	Linear	0.0573
	18	0.55	1.02	1.52	1.92	2.48	165	89.3	142	98.4	Linear	0.0405
	12	0.32	0.58	0.86	1.08	1.38	216	89.5	141	93.4	Linear	0.0260
b	24	0.75	1.77	2.52	3.13	3.56	126	90.8	137	99.5	Logistic	0.0568
	18	0.55	1.26	1.78	2.20	2.48	157	90.2	142	98.9	Logistic	0.0397
	12	0.32	0. 72	1.00	1.23	1.38	193	91.0	141	95.0	Logistic	0.0221
c	24	0.75	0.91	1.08	1.59	3.56	166	89.2	137	96.7	Exponential	0.3456
	18	0.55	0.66	0.78	1.14	2.48	192	90.3	142	92.6	Exponential	0.2464
	12	0.32	0.38	0.45	0.65	1.38	248	90.5	141	86.8	Exponential	0.1403
d	24	0.75	2.04	3.56	2.52	1.59	150	90.5	137	92.4	Quadratic	0.1859
	18	0.55	1.45	2.48	1.78	1.14	175	89.9	142	87.1	Quadratic	0.4222
	12	0.32	0.82	1.38	1.00	0.65	229	87.7	141	78.8	Quadratic	0.2259

In particular, we also apply the sample size determined in multi-arm design by the log-rank global test proposed by Barthel et al. ⁷ as a comparison in Tables 2 and 3. Notice that many sample size determination procedures for two treatments’ comparison have been proposed for survival outcomes. For example, a parametric way (exponential survival model-based) from Lachine and Foulkes ³¹ and a Cox PH model-based method from Freedman. ³² The former accommodates staggered patient entry and loss to follow-up but the latter does not. To illustrate the better performance of the proposed multi-arm design with these conventional methods, we apply Bonferroni correction on these two methods to derive the sample size as well as the statistical power. In the study design, the postulated HR is based on Groups 1 and 2. The sample size is validated through the test for equivalence of exponential survival rate (Lachin ³³ ) and log-rank test for the two methods, respectively.

4.1. Global test under PHs

The simulation results for the PHs setting are presented in Table 2. We considered $K=3$ and $4$ with varying RMST margins $\mathbf{\Delta }$ (effect size) between each pair of treatment and control, while the restriction time was fixed at 24 and 18 months. We set the representative treatment arm to be the second arm and assumed exponential survival for all treatment groups. By varying the RMST and the corresponding HRs (listed in brackets in Table 2) in different arms, we investigated how the violation of the equal variance assumption affected the power and sample size. We determined the sample size based on (12), which adjusts for staggered patient entry, and then conducted the global test proposed in Section 2.1 to investigate the empirical power under the determined sample size. The simulation results for conventional methods with Bonferroni correction are shown in Table S1 in the Supplemental materials.

We can observe that when the treatment difference is moderate, the power can be well maintained at around 90% under the calculated sample size, which verifies the validity of the assumption of equal RMST variance across treatment groups. Basically, the smaller the variation among treatment effects, the closer the empirical power to the target. Compared with the method proposed by Barthel et al., ⁷ the RMST method is more accurate and reduces the cost in most cases, while the log-rank-based method is conservative, especially when the analysis time is short. The required sample size is influenced by a combined effect of the effect size $\mathbf{\Delta }$ and how many groups share the same RMST. Roughly speaking, a larger effect size (difference between treatment and control) along with a larger variation of the RMST among different treatment arms leads to a smaller sample size. Further, the restriction time $\tau$ has a large impact on the sample size required. Under the same survival function, there is a $15\mathrm{\%}$ to $20\mathrm{\%}$ increase in sample size when $\tau =18$ compared with $\tau =24$ . The last row in the table of each $\tau$ shows the type I error, which is well maintained at around $0.05$ . From Table S1 in the Supplemental Material, we can tell that the Bonferroni correction under Lachin and Foulkes’ method often necessitates a larger sample size due to its stringent model assumptions, but the power is well maintained because the sample size determination accommodates staggered patient entry and loss to follow-up. On the other hand, Freedman’s sample size calculation is less conservative than Lachin and Foulkes’ method but more conservative than the proposed method. Hence the empirical power is larger than it designed because it is not able to accommodate staggered patient entry as well as the loss to follow-up. As these situations often exist in real clinical trials, the proposed sample size determination based on the global RMST test is more appropriate than the multiple correction methods.

4.2. Global test under NPHs

The simulation results for a wide range of NPH scenarios are presented in Table 3. To simulate the NPH scenarios, we used a piece-wise exponential distribution. We considered a 3-arm design, and the hazards in the treatment arms varied over different time periods, leading to different overall $\mathbf{\Delta }$ . We calculated the sample size at $\tau =24$ . The column named by ‘Adjusted for Non-PH” presents the method of determining sample size using parameters specified in Table 3, which takes NPHs into account. The representative treatment group is still Group 2. To facilitate comparison, the last four columns in the table show the results when the PHs assumption are adopted without considering changes in HR over time, and the underlying assumed HR is fixed at the minimum HR in each group, thus resulting in a smaller sample size required and a very low power. With the parameters designed, we simulated four situations, including early effect, late effect, middle effect, and crossing survival functions. These situations are common in targeted therapy and cancer immunotherapy. For example, late divergence of the survival curves occurs when there is a bridging therapy in each arm before the action of a therapy, while early separation of survival curves occurs when the drug is effective at the beginning, but patients experience a leveling off thereafter. The simulation results for conventional methods with Bonferroni correction are shown in Table S2 in the Supplemental Materials.

In summary, the simulation results in Table 3 demonstrate the effectiveness of the proposed RMST method in accurately estimating the required sample size and achieving the targeted power under a wide range of NPHs scenarios. The results also highlight the importance of considering the NPH assumption in the design stage, as using the PH assumption can lead to a severe underestimation of sample size, especially when crossing survivals occurs at a relatively early time before $\tau$ . Furthermore, we can observe that under the same RMST margin ${\mathrm{\Delta }}_k,k=2,3$ , the early effect case requires fewer samples than the middle effect, while in the late effect situation, RMST design requires a larger sample size. This is reasonable since in a delayed effect situation, the potential differences of survival curves are more likely to be partially taken into account in the followed-up $\tau$ units of time. In this example, Barthel’s sample size can also maintain the power well but is slightly inferior to the proposed RMST test with regard to the robustness. From Table S2 in the Supplemental material, we can tell that Freedman’s sample size is the least conservative while Lachin and Foulkes’ sample size is the most conservative, but they both face severe power dilution when NPHs exist. The simulation emphasizes the importance of using RMST-based tests and sample size determination as it accommodates NPHs.

In summary, the sample size determination method is quite robust to the equal variance assumption in treatment arms, both in the constant hazard and non-constant hazard cases, and is more accurate and cost-saving than that of Barthel et al. ⁷ These findings highlight the importance of considering the potential violation of the proportional hazards assumption and using RMST as a tool when designing clinical trials and calculating sample size efficiently.

Table 5.

Prior dose–response shapes with $d\in \{0,1,\dots ,5\}$

Model	$f_{12}(d;\theta )$	$f_{18}(d;\theta )$	$f_{24}(d;\theta )$
Linear	$9+0.3d$	$10+0.5d$	$11+0.7d$
Exponential	$9+0.1\exp (d/2)$	$10+0.2\exp (d/2)$	$11+0.3\exp (d/2)$
Logistic	$9+2/(1+\exp ((2.2-d)/0.5))$	$10+2.8/(1+\exp ((2-d)/0.5))$	$11+3.6/(1+\exp ((1.8-d)/0.5))$
Quadratic	$9+0.8d-0.1d^2$	$10+1.4d-0.2d^2$	$11+2d-0.3d^2$
Emax	$9+2d/(2.2+d)$	$10+2.8d/(2+d)$	$11+3.6d/(1.8+d)$

4.3. MCP-Mod with endpoint measured by RMST under dose-finding studies

In this subsection, we study the performance of MCP-Mod for endpoint measured by RMST in a dose-finding study via simulation. We investigated six dose levels $d=0,1,2,3,4,5$ , with RMST measured at three restriction time points $\tau =24,18,12.$ Five different dose–response models for RMST were assumed prior to data generation and are given in Table 5. All the shapes are within a similar range. Four different cases with different HRs in each group are investigated: (a) $\mathrm{HR}=(0.9,0.82,0.74,0.68,0.6)$ ; (b) $\mathrm{HR}=(0.9,0.78,0.7,0.64,0.6)$ ; (c) $\mathrm{HR}=(0.9,0.88,0.86,0.8,0.6)$ ; and (d) $\mathrm{HR}=(0.9,0.75,0.6,0.7,0.8)$ . Figure 1 displays the true dose–response shapes for all scenarios regarding RMSTs, where the HRs on the optimal dose are all at 0.6. Under these settings, by adopting a conservative summary power function ${\pi }^{*}(n)=min\{{\pi }_1(n),\dots ,{\pi }_M(n)\}$ , we calculate the required sample size using the proposed method for the purpose of establishing the PoC in Table 4. The empirical power under the required sample size and model selection results via the minimum $p$ -value approach are shown in the table as well. Then given the selected model, we fit the corresponding nonlinear parametric model and summarize the standard errors for nonlinear least squares in the table. To compare the sample size under the proposed two methods with and without modeling, we also compute the sample size determined through a global test. Basically, by assuming candidate models of the dose–response shape in the design stage, the sample size needed can be dramatically reduced by at most $38\mathrm{\%}$ .

Figure 1.

True dose–response shapes in simulation studies.

The required sample size relies on the collective effect of model parameters, summary power function, restriction times, and the margin of the effects. Generally speaking, the restriction time alone has little effect on the sample size while other factors affect the sample size a lot. When the candidate model fits the data well, power can be better maintained at the target level of 0.9. Here, “shape” means the shape of the curve as well as the absolute difference between the true curve and candidate curves. For example, in Case d, though the candidate shapes are not exactly the same as the true dose–response shape shown, the absolute difference between the assumed RMST level and true RMST are quite close when $\tau =24$ , thus resulting in a well maintained power level. On the other hand, the good thing is that the model selection performance is good no matter how the candidate models fit the true dose–response relationship.

5. Trial examples

Two real clinical studies with right-censored data are analyzed in this section to illustrate the proposed sample size determination procedure.

5.1. CheckMate 451 trial

CheckMate 451 trial, a double-blind, three-arm, phase III study, ³⁴ with an aim to evaluate the efficacy of nivolumab plus ipilimumab and nivolumab monotherapy as maintenance therapy for extensive-disease small-cell lung cancer following first-line chemotherapy. In total, 834 patients were randomly assigned in a $1:1:1$ ratio to the nivolumab plus ipilimumab arm, nivolumab arm, and placebo arm. However, there was a three-month delay between the patients’ last chemotherapy dose and random drug allocation, resulting in NPHs of the KM curves. The survival function of the placebo arm was modeled using four hazard pieces, with assumed overall survival (OS) rates of 90%, 47%, 15%, and 9% at 3, 9, 18, and 26 weeks, respectively. For the combination of nivolumab and ipilimumab and nivolumab alone, a delayed effect versus placebo was assumed, with a HR of 1 for the first 3 months and an HR of 0.68 thereafter. The study also considered a 5% probability of dropout by month 6, and patients were uniformly accrued over 28 weeks. A HR of 0.72 was deemed clinically relevant, and the trial was designed with 90% power and 0.05 type I error, with a total follow-up duration of $t_2=36$ weeks.

To determine the sample size, we used the same setup and assumed that the nivolumab group shares the same RMST variance with the combined therapy group. Instead of using HRs, we first transformed the assumed HRs into corresponding RMSTs. We then determined the restriction time $\tau$ using the “percent maturity” method proposed by Royston and Parmar, ³⁵ to select the minimum $\tau$ such that the percent maturity reaches 1. In this trial, we chose $\tau =21$ . We also calculated the sample size per arm with varying $\tau$ from 7 to 30 months in Figure 2. The figure illustrates that when the restriction time is >18 months, the required sample size can be smaller than that in the original trial design. This means that it may be possible to achieve the same statistical power with a smaller sample size, which can help to reduce the cost and time of the trial. However, it is important to note that extending the analysis time may also increase the risk of confounding factors and other biases, which should be carefully considered when making such decisions.

Figure 2.

Sample size per arm under different $\tau$ in CheckMate 451 trial. Dashed line: $n=278$ in the original trial design.

Under these assumptions, our proposed global RMST test method predicts that a total sample size of 750 patients is sufficient to detect a clinically relevant difference vector in RMST of $\mathrm{\Delta }=(1.829,0)$ months, reducing the original sample size by 10%. Moreover, even though the treatment variances are not the same as assumed in the setup, the empirical power from 1000 simulations under $n=750$ of detecting a clinically relevant RMST at 21 months still reaches 91.1%. To validate the required sample size, we re-calculate the sample size using the same $\tau$ value, but with re-estimated design parameters based on their published results. These results correspond to the overall survival (OS) rates of 92%, 53%, 23%, and 15% at 3, 9, 18, and 26 weeks, respectively. Additionally, the HRs were found to be 0.92 and 0.84 for nivolumab versus placebo, and nivolumab plus ipilimumab versus placebo, respectively. The re-calculated sample size required was calculated to be 791, which does not deviate a lot from the sample size determined by the designed parameters.

It is important to note that the primary goal of the original trial is to detect the difference between the nivolumab plus ipilimumab arm and the control arm. Under the proposed sample size determination method, we found that with a size of 250 per arm, the empirical power of detecting the clinical RMST difference between the control and nivolumab plus ipilimumab arm in a pairwise manner at 21 months can also reach 93.2%. This justifies the validity of applying our proposed sample size method.

5.2. MCP-Mod with RMST in COPD example

In this example, we analyzed data from a COPD study to demonstrate the use of MCP-Mod for endpoint measured by RMST based on the time to the first COPD exacerbation. We adopted the same design as in Deng et al., ¹⁰ where four once-daily doses of the study medication were considered at 5, 25, 50, and 100 mg, along with a placebo, with an equal allocation of patients in each dose group. The desired type I error and targeted power were set to 0.05 and 85%, respectively. Since the median survival time in the placebo group was 0.5 year, we set the restriction time at $\tau =1$ .

To reconstruct the same six candidate models shown in Figure 3 between dose and RMST, we back-transformed from hazard to RMST under exponential survival, instead of log hazard used in Deng et al. ¹⁰ Using the proposed method, we found that a total of 405 patients were required with $n=81$ patients per arm to achieve 85% power at $\tau =1$ . This is slightly larger than the sample size of $n=75$ per arm used in Deng et al. ¹⁰ However, this is a reasonable tradeoff as the proposed RMST-based design allows for more flexibility and robustness to distribution assumption than HR-based designs.

Figure 3.

Candidate dose–response relationships for the chronic obstructive pulmonary disease (COPD) study: black line: linear model with $E_0=0.541,\gamma =0.0014$ ; red line: Emax model with $E_{50}=43.81,E_{max}=0.1986$ ; orange line: Emax model with $E_{50}=5.6830,E_{max}=0.1458$ ; blue line: exponential model with $E_1=0.002302,\gamma =24.38$ ; purple line: logistic model with $E_{max}=0.1381,E_{50}=39.64,\gamma =6.96$ ; green line: beta model with $a=0.71,b=0.99,D=120$ .

To assess the accuracy of our proposed sample size calculation procedure, we conducted 1000 simulations with 81 patients per arm and calculated the empirical and theoretical power for each candidate model, as well as the overall average power when $\tau =1$ . Table 6 displays all the results, comparing them to those from the MCP-Mod proposed by Deng et al. ¹⁰

Table 6.

Theoretical and empirical power of MCP-Mod for RMST endpoint comparing with that based on the hazard.

	Hazard		RMST
Method	Empirical power	Theoretical power	Empirical power	Theoretical power
Linear	0.827	0.827	0.848	0.844
Emax1	0.857	0.863	0.872	0.860
Emax2	0.881	0.872	0.855	0.850
Exponential	0.828	0.878	0.843	0.851
Logistic	0.920	0.917	0.911	0.927
Beta	0.804	0.851	0.781	0.781
Average	0.851	0.851	0.852	0.852

^a

MCP-Mod: multiple comparison procedures combined with modeling techniques; RMST: restricted mean survival time.

The empirical power under all models, except for exponential and logistic, was higher than the theoretical power in the proposed MCP-Mod based on RMST, resulting in an overall empirical power of 88.7%. This indicates that our proposed method performs well in practice. The difference between theoretical power and empirical power is within 5%, which could be attributed to the instability of the RMST estimator and variance estimator.

We also computed the empirical type I error when the RMST was the same for all treatment groups as that in the placebo group. The resulting size was 0.048, which is the same as that reported by Deng et al. ¹⁰ This example well illustrates the potential of using RMST instead of HR in a multi-arm dose-finding study design. With this approach, the designed sample size can still achieve the targeted power, while providing a more interpretable dose–response relationship.

Notice that one concern of using the RMST test in real-world examples is its long follow-up time. However, it is not specific to the RMST-based method but rather a characteristic of other dose-finding studies involving late-onset toxicities. To provide further empirical evidence and address the concern, we compared the proposed RMST-based method with the Cox model-based method. ¹⁰ Based on the study design, the trial would have been terminated at the restriction time of 1 year using the proposed method. While in Deng et al., ¹⁰ the trial was terminated when a total of 242 events were observed (Type II censoring design). To unify the termination time, we simulated the data 1000 times and observed that the average termination time using the method from Deng et al. ¹⁰ with 242 observed events was close to 0.97 year, which is in proximity to the restriction time of 1 year used in our proposed method. This finding highlights the similarity in the follow-up times required by different methodologies when dealing with right-censored survival data.

6. Discussion

In this article, we present the multi-arm design and sample size calculation for survival data with endpoint measured by RMST with the aims of detecting global effect and modeling dose–response relationship, respectively. The main contribution lies in its first attempt to apply RMST in multi-arm design for survival data. By utilizing RMST, the resulting sample size can achieve a more accurate power target. This is due to the stability of the RMST ratio variance, which outperforms the conventional HR approach. We propose a closed form and a simulation-based sample size calculation procedure in realistic settings, which is quite flexible allowing staggered patient entry as well as NPHs that commonly occur in modern drug development. The proposed method requires fewer a priori assumptions, as it only necessitates the survival curves of the control and representative treatment group, rather than assuming the distributions of all treatments. Complex trial settings in survival data have a combined effect on sample size, including the magnitude of the targeted RMST difference to be detected, the accrual distribution, the control group survival distribution, a representative treatment group distribution, loss to follow-up, length of the study and follow-up duration, restriction time, the allocation ratio between control and treatment, and Type I and Type II errors. A related problem is how we plan a trial with an appropriate restriction time $\tau$ in the study protocol. Despite the clinical and scientific perspective, several data-dependent methods can be directly applied in multi-arm design, such as the method developed in Eaton et al. ¹⁴ based on the probability of RMST to be estimable, or the probability that data is mature enough for final analysis proposed by Royston and Parmar. ³⁵ The selection of restriction time $\tau$ indeed affects the sample size as well as the performance of the RMST test a lot. In this paper, due to the complexity of multi-arm trial design, we propose to set a priori specification of $\tau$ in the trial protocol for the primary analysis, guided by clinical considerations. This is in line with the recommendations by Royston and Parmar. ³⁵ In practice, many issues need to be considered to select an appropriate $\tau$ . First, there is a tradeoff between the choice of $\tau$ and the sample size. When patient enrollment is easily achievable, a smaller $\tau$ can be selected; otherwise, the sample size may be small, and a relatively longer $\tau$ should be adopted. Secondly, the nature of the disease also influences the appropriate follow-up periods as a certain cancer may require a longer follow-up period. In addition, practical constraints, such as financial resources, must also be taken into account, meaning that $\tau$ cannot be excessively long with a limited budget. Of course, there are other data-driven approaches to select $\tau$ in the absence of this clinically relevant choice as well. For example, Royston and Parmar ³⁵ also suggested selecting a value of $\tau$ slightly lower than the maximum expected follow-up time.

Simulation results of various situations show that the proposed sample size determination method performs well in acquiring the targeted power and outperforms the conservative HR-based sample size formula. It also validates the emphasis on taking particular care of the parameter settings in the study design stage, otherwise, large bias and power loss may occur. Trial examples we performed show that the proposed RMST-based method with an appropriate choice of the restriction time $\tau$ requires a similar sample size as that in the HR-based method, but can accommodate much more complicated survival curves and better maintains the power level required than the HR-based method.

We have focused on the multi-arm design analyzed at a pre-specified restriction time $\tau$ or at the end of the trial. In practice, for the purpose of modifying the trial design and re-estimating the sample size, an interim analysis is often conducted in the ongoing trial. Specifically, the group sequential method in randomized clinical trials based on the difference in RMST was studied theoretically by Lu and Tian, ¹⁶ where the test statistic for two-arm comparison has an independent increments covariance structure. In the context of multi-arm global design, to develop sequential analysis at interim times $t_1$ and $t_2$ , the joint distribution is no longer a multivariate normal but a multivariate chi-square, and the theorem will not be directly applicable, which requires further exploration. In addition, other multi-arm designs, such as factorial design, multi-arm multi-stage design, or adaptive design have not been considered yet and our work opens opportunities for future development of these novel designs.

Supplemental Material

sj-pdf-1-smm-10.1177_09622802231219852 - Supplemental material for Sample size calculation for multi-arm parallel design with restricted mean survival time

Supplemental material, sj-pdf-1-smm-10.1177_09622802231219852 for Sample size calculation for multi-arm parallel design with restricted mean survival time by Yaxian Chen, Kwok Fai Lam and Jiajun Xu in Statistical Methods in Medical Research