Power and Sample Size Calculations for the Restricted Mean Time Analysis of Prioritized Composite Endpoints

Abstract

As a new way of reporting treatment effect, the restricted mean time in favor (RMT-IF) of treatment measures the net average time the treated have had a less serious outcome than the untreated over a specified time window. With multiple outcomes of differing severity, this offers a more interpretable and data-efficient alternative to the prototypical restricted mean (event-free) survival time. To facilitate its adoption in actual trials, we develop simple approaches to power and sample size calculations and implement them in user-friendly R programs. In doing so we model the bivariate outcomes of death and a nonfatal event using a Gumbel–Hougaard copula with component-wise proportional hazards structures, under which the RMT-IF estimand is derived in closed form. In a standard set-up for censoring, the variance of the nonparametric effect-size estimator is simplified and computed via a hybrid of numerical and Monte Carlo integrations, allowing us to compute the power and sample size as functions of component-wise hazard ratios. Simulation studies show that these formulas provide accurate approximations in realistic settings. To illustrate our methods, we consider designing a new trial to evaluate treatment effect on the composite outcomes of death and cancer relapse in lymph node-positive breast cancer patients, with baseline parameters calculated from a previous study.

Keywords:

• Clinical trials
• Composite endpoints
• Copula models
• Estimands
• Restricted mean survival time

1 Introduction

Robust metrics like the restricted mean survival time (RMST) are increasingly welcomed in the design and analysis of clinical trials (Royston and Parmar 2011, 2013; Uno et al. 2014; McCaw et al. 2019). Unlike the traditionally used hazard ratio, the RMST, which measures the average survival time gained within a fixed time window, remains interpretable and nonparametrically estimable whether or not the hazards are proportional over time. When nonfatal events like hospitalization or cancer relapse are added to the survival endpoint, Mao (2021) considered an extension called the restricted mean time in favor (RMT-IF) of treatment, which measures the average time gained having had a less severe outcome. The RMT-IF draws on all events but also properly differentiates fatal versus nonfatal ones in their severity. As a result, it is statistically more efficient and clinically more interpretable than the standard restricted mean event-free survival time (RMEST), the RMST applied to the first composite event (Freemantle et al. 2003; Anker and McMurray 2012).

To adopt a new metric in real trials, we need to be able to calculate its power and sample size. While the specifics necessarily depend on the metric itself, such calculations always revolve around two central quantities: an effect-size estimand for the between-group difference and a variance measure for the within-group variation. Collectively the “signal” and “noise” give rise to the noncentrality parameter under the alternative hypothesis, upon which the power and sample size depend (see e.g., Ch. 14 of van der Vaart 1998). In this regard the RMST enjoys a natural advantage—its estimand is a simple integral of the survival function, and the corresponding plug-in Kaplan–Meier estimator can be easily tackled in the martingale framework (Fleming and Harrington 1991), resulting in a closed-form variance (Tian et al. 2018; Yung and Liu 2020; Eaton et al. 2020).

Unfortunately, the RMST’s ease with such derivations does not extend to composite endpoints. Although the RMT-IF estimand can be similarly expressed as an integral of the survival functions for the component events, the correlations between them make the resulting estimator resistant to standard martingale theory (Wei et al. 1989; Lin et al. 2000). In fact, Mao (2021) derived the variance of the RMT-IF estimator only through its robust influence functions (Tsiatis 2006), without an analytic form. Moreover, to describe the multiple outcomes realistically, one would need a joint model that accommodates not only between-component correlations but also, ideally, component-specific treatment effects.

Similar challenges exist in the power and sample size calculations for the win ratio (Pocock et al. 2012), a pairwise comparison scheme to analyze composite endpoints, and have recently been addressed by Mao et al. (2021) via a combination of copula models (Oakes 1989), U-statistic theory (Hoeffding 1948), and numerical computations. Although the RMT-IF is also pairwise defined, its estimand is reducible to individual-based quantities (namely, the component-wise survival functions). This makes variance calculation easier than it is for win-ratio-like U-statistics (Mao et al. 2021). More importantly, unlike the win ratio (Bebu and Lachin 2016; Mao 2019; Dong et al. 2020), the RMT-IF as an estimand is defined through the uncensored outcomes (as is the RMST) (Akacha et al. 2017; Akacha, Bretz and Ruberg 2017; ICH 2020). This means that we can calculate the effect size directly from the outcome distributions without numerical averaging over the censoring distribution (Mao et al. 2021). On the other hand, the RMT-IF involves an extra parameter, that is, the restricting time, which affects the interpretation of the effect size as well as the operating characteristics of the test. A flexible framework for trial design should allow the user to vary the choice of this important parameter.

We aim to develop tools that are fully aligned with the above features of the RMT-IF. Section 2 lays out the main steps. We start off with a bivariate copula model under which the survival and nonfatal event times each follow a proportional hazards model with constant baseline rate. We allow the treatment hazard ratio (HR) to be different on death and the nonfatal event. The component-specific HRs are translated to the effect size on the RMT-IF by an analytic function, yielding the “signal” of the test. For the “noise” part, we simplify the robust influence functions (Mao 2021) into an integral of martingale-like processes under a standard set-up for censoring. The restricting time appears in the upper limit of this integral. Thanks to this expression, we only need a one-off numerical routine to compute the variances for a spectrum of restricting times. These quantities not only provide us ready estimates of power and sample size, but can also be used to retrieve other design information such as the length of patient accrual, which simultaneously affects the power (via administrative censoring) and sample size. The proposed methods are vetted by simulations in Section 3 and are illustrated on the design of a future breast cancer trial in Section 4 using baseline parameters estimated from an existing study. Section 5 concludes the article with some practical guidelines and discussions of future work.

2 Methods

2.1 Set-Up

Let $D^{(a)}$ denote the survival time and $T^{(a)}$ the nonfatal event time for a generic patient in group a, where a = 1 and 0 indicate the treatment and control, respectively. To derive usable formulas for RMT-IF’s effect size and variance, we need a simple and realistic joint model for $(D^{(a)},T^{(a)})$. To that purpose, consider the Gumbel–Hougaard copula model, a popular choice for composite outcomes (Oakes 1989; Luo et al. 2015; Oakes 2016; Mao et al. 2021). Specifically, let(1) $$\begin{array}{c}\text{Pr}(D^{(a)}>s,T^{(a)}>t)\\ \hspace{1em}=\exp \left(-{\left[{\left\{\exp (a{\theta }_D){\lambda }_{D}s\right\}}^{\kappa }+{\left\{\exp (a{\theta }_H){\lambda }_{H}t\right\}}^{\kappa }\right]}^{1/\kappa }\right),\end{array}$$(1) where $({\lambda }_D,{\theta }_D)$ and $({\lambda }_H,{\theta }_H)$ are parameters related to the marginal distributions of $D^{(a)}$ and $T^{(a)}$, respectively, and $\kappa \ge 1$ determines their Kendall’s rank correlation $\mathcal{C}$ by $\mathcal{C}=1-{\kappa }^{-1}$ (Oakes 1989). Clearly, model (1) implies marginal exponential distributions for both death and the nonfatal event, the former with baseline hazard rate λ_D and treatment-to-control HR $\exp ({\theta }_D)$, and the latter with baseline hazard rate λ_H and treatment-to-control HR $\exp ({\theta }_H)$. Of note, we posit (1) only as a working model to derive the properties of the original nonparametric estimator for the RMT-IF (Mao 2021), not to construct a parametric one based on this model.

2.2 Expression of Effect Size

The RMT-IF was originally formulated on multistate outcomes $Y^{(a)}(t)$ from group a (a = 1, 0) with state space $\{0,1,\dots ,K,\infty \}$ for some $K\in {\mathbb{Z}}^{+}$. A higher number indicates a worse state and $\infty$ indicates death. Let $\tau >0$ be a prespecified restricting time. By generalized pairwise comparison (Buyse 2010), the average time patients in group a fare better than those in group $1-a$ over $[0,\tau ]$ is(2) $$\begin{array}{c}{W}_{a,1-a}(\tau )=E\left[{\int }_0^{\tau }I\{Y^{(a)}(t)<Y^{(1-a)}(t)\}\mathrm{d}t\right]\\ ={\int }_0^{\tau }\text{Pr}\{Y^{(a)}(t)<Y^{(1-a)}(t)\}\mathrm{d}t,\end{array}$$(2) where $I(·)$ is the indicator function. To measure the net favorability of the treatment against control, the RMT-IF considers the difference $\mu (\tau )=W_{1,0}(\tau )-W_{0,1}(\tau )$. We can interpret $\mu (\tau )$ as the net average time gained by the treatment in a more favorable state as compared to the control. For a two-state life-death model, we can show that $\mu (\tau )$ reduces to the group difference in the RMST. Like the standard RMST, $\mu (\tau )$ is defined through uncensored outcomes. This ensures that the estimand has a scientifically consistent meaning that does not depend on specific studies.

In our set-up with death prioritized over the nonfatal event, we may write $Y^{(a)}(t)=I(T^{(a)}\le t)+\infty ·I(D^{(a)}\le t)$. Hence, $Y^{(a)}(t)=0,1,$ and $\infty$ if the patient has been event-free, experienced the nonfatal event but remained alive, and died, respectively, by time t. As a result, $Y^{(a)}(t)<Y^{(1-a)}(t)$ means one of two mutually exclusive scenarios: (i) event-free versus having experienced the nonfatal event but alive, that is, ${\tilde{T}}^{(a)}>t$ and $T^{(1-a)}\le t<D^{(1-a)}$, where ${\tilde{T}}^{(a)}=D^{(a)}\wedge T^{(a)}$ (with $b\wedge c=\min (b,c)$); (ii) alive (regardless of the nonfatal event) versus dead, that is, $D^{(a)}>t$ and $D^{(1-a)}\le t$. By (2), we have that$$\begin{array}{c}{W}_{a,1-a}(\tau )={\int }_0^{\tau }\text{Pr}({\tilde{T}}^{(a)}>t,T^{(1-a)}\le t<D^{(1-a)})\mathrm{d}t\\ \hspace{1em}+{\int }_0^{\tau }\text{Pr}(D^{(a)}>t,D^{(1-a)}\le t)\mathrm{d}t\\ ={\int }_0^{\tau }{\tilde{S}}^{(a)}(t)\left\{S^{(1-a)}(t)-{\tilde{S}}^{(1-a)}(t)\right\}\mathrm{d}t\\ \hspace{1em}+{\int }_0^{\tau }{S}^{(a)}(t)\left\{1-S^{(1-a)}(t)\right\}\mathrm{d}t,\end{array}$$where $S^{(a)}(t)=\text{Pr}(D^{(a)}>t)$ and ${\tilde{S}}^{(a)}(t)=\text{Pr}({\tilde{T}}^{(a)}>t)$. Hence,(3) $$\begin{array}{c}\mu (\tau )=W_{1,0}(\tau )-W_{0,1}(\tau )\\ ={\int }_0^{\tau }\left\{{\tilde{S}}^{(1)}(t)S^{(0)}(t)-{\tilde{S}}^{(0)}(t)S^{(1)}(t)\right\}\mathrm{d}t\\ \hspace{1em}+{\int }_0^{\tau }\left\{S^{(1)}(t)-S^{(0)}(t)\right\}\mathrm{d}t.\end{array}$$(3)

The second term on the far right hand side of (3) is just an alternate expression of the group difference in the standard RMST, that is, $E(D^{(1)}\wedge \tau )-E(D^{(0)}\wedge \tau )$ (Royston and Parmar 2011).

We have already established that $S^{(a)}(t)=\exp \left\{-\exp (a{\theta }_D){\lambda }_{D}t\right\}$. It is also easy to see from (1) that ${\tilde{S}}^{(a)}(t)=\text{Pr}(D^{(a)}>t,T^{(a)}>t)=\exp \left\{-\exp (a\theta )\lambda t\right\}$, where $\lambda ={({\lambda }_D^{\kappa }+{\lambda }_H^{\kappa })}^{1/\kappa }$ and $\theta ={\kappa }^{-1}\log \left\{\exp (\kappa {\theta }_D){\lambda }_D^{\kappa }+\exp (\kappa {\theta }_H){\lambda }_H^{\kappa }\right\}-\log (\theta )$. Inserting these expressions into (3) and carrying out the integration, we find that(4) $$\begin{array}{c}\mu (\tau )=f\left\{\tau |\exp (\theta )\lambda +{\lambda }_D,\lambda +\exp ({\theta }_D){\lambda }_D\right\}\\ \hspace{1em}+f\left\{\tau |\exp ({\theta }_D){\lambda }_D,{\lambda }_D\right\},\end{array}$$(4) where $f(\tau |x,y)=x^{-1}\left\{1-\exp (-x\tau )\right\}-y^{-1}\left\{1-\exp (-y\tau )\right\}$. By similar derivation, we can express the group difference in RMST and RMEST as(5) $$\begin{array}{c}{\mu }_D(\tau )=E(D^{(1)}\wedge \tau )-E(D^{(0)}\wedge \tau )=f\left\{\tau |\exp ({\theta }_D){\lambda }_D,{\lambda }_D\right\}\\ \text{and}\\ {\mu }_E(\tau )=E({\tilde{T}}^{(1)}\wedge \tau )-E({\tilde{T}}^{(0)}\wedge \tau )=f\left(\tau |\exp (\theta )\lambda ,\lambda \right),\end{array}$$(5) respectively.

2.3 Test Statistic and Null Variance

In practice, the outcome events may be censored due to study termination or random loss to follow-up. Denote the censoring time of the generic patient in group a by $C^{(a)}$ (a = 1, 0), which is assumed to be independent of $(D^{(a)},T^{(a)})$. We then observe $({\tilde{X}}^{(a)},{\tilde{\delta }}^{(a)},X^{(a)},{\delta }^{(a)})$, where ${\tilde{X}}^{(a)}={\tilde{T}}^{(a)}\wedge C^{(a)},{\tilde{\delta }}^{(a)}=I({\tilde{T}}^{(a)}\le C^{(a)}),X^{(a)}=D^{(a)}\wedge C^{(a)}$, and ${\delta }^{(a)}=I(D^{(a)}\le C^{(a)})$. Given a random n_a -sample of $({\tilde{X}}^{(a)},{\tilde{\delta }}^{(a)},X^{(a)},{\delta }^{(a)})$ for group a (a = 1, 0), one can estimate ${\tilde{S}}^{(a)}(t)$ and $S^{(a)}(t)$ by the usual Kaplan–Meier estimator using the subsamples for $({\tilde{X}}^{(a)},{\tilde{\delta }}^{(a)})$ and $(X^{(a)},{\delta }^{(a)})$, respectively, and then insert the estimates in (3) to obtain a nonparametric estimator $\widehat{\mu }(\tau )$ for $\mu (\tau )$.

Let ${\widehat{\sigma }}^2$ denote the influence function-based variance estimator for $n^{1/2}\{\widehat{\mu }(\tau )-\mu (\tau )\}$ given in Proposition 1 of Mao (2021), where $n=n_1+n_0$ is the total sample size. Then, a two-sided level-α test $(0<\alpha <1)$ rejects the null hypothesis of no treatment effect (under which ${\theta }_D={\theta }_H=0$ so that $\mu (\tau )=0$) if $S_n:=n^{1/2}\widehat{\sigma }{(\tau )}^{-1}|\widehat{\mu }(\tau )|>z_{1-\alpha /2}$, where $z_{1-\alpha /2}={\Phi }^{-1}(1-\alpha /2)$ and $\Phi (·)$ is the cumulative distribution function of the standard normal distribution.

The power of this test depends on the distribution of S_n under the alternative hypothesis H_A , when θ_D and/or θ_H are nonzero. Under a standard asymptotic setting in which ${\theta }_D,{\theta }_H=O(n^{-1/2})$, we can use the normal approximation $S_n\stackrel{H_A}{\sim }N\left\{n^{1/2}{\sigma }_0{(\tau )}^{-1}\mu (\tau ),1\right\}$, where ${\sigma }_0^2(\tau )$ is the asymptotic variance of $n^{1/2}\widehat{\mu }(\tau )$ under the null (see e.g., sec. 14.1 of van der Vaart 1998). Assume without loss of generality that $\mu (\tau )>0$. Then the rejection rate under H_A is approximated by(6) $${\text{Pr}}_{H_A}(S_n>z_{1-\alpha /2})=\Phi \left\{n^{1/2}{\sigma }_0{(\tau )}^{-1}\mu (\tau )-z_{1-\alpha /2}\right\}.$$(6)

Conversely, we can use (6) to solve for the sample size n needed for a target power $1-\beta$:(7) $$n=\frac{{\sigma }_0^2(\tau ){(z_{1-\alpha /2}+z_{1-\beta })}^2}{\mu {(\tau )}^2}.$$(7)

In both (6) and (7), $\mu (\tau )$ can be calculated through (4) as a function of ${\theta }_D,{\theta }_H,{\lambda }_D,{\lambda }_H$, and κ. It thus remains to calculate ${\sigma }_0^2(\tau )$. Unlike the estimand, the variance depends on the censoring distribution, which is heavily influenced by trial design. As a standard set-up, suppose that patients are recruited uniformly over an initial period $[0,{\tau }_b]$ and are followed until time τ_c , where ${\tau }_c>{\tau }_b$. This results in administrative censoring times (those due to study termination) distributed as $\text{Un}[{\tau }_c-{\tau }_b,{\tau }_c]$. In addition, we assume a constant hazard λ_L for random loss to follow-up. Combining the two mechanisms together, we have that $C^{(a)}\sim \text{Un}[{\tau }_c-{\tau }_b,{\tau }_c]\wedge \text{Expn}({\lambda }_L)$ so that $\text{Pr}(C^{(a)}>t)=[\{{\tau }_b^{-1}({\tau }_c-t)\}\wedge 1]\exp (-{\lambda }_{L}t)$.

With that we can express in closed form the key quantities in the influence functions for variance estimation. Of particular importance are the “at-risk” probabilities$$\begin{array}{c}{\tilde{\pi }}^{(0)}(t)=\text{Pr}({\tilde{X}}^{(0)}>t)=\left(\frac{{\tau }_c-t}{{\tau }_b}\wedge 1\right)\exp \{-({\lambda }_L+\lambda )t\}\\ \text{and}\\ {\pi }^{(0)}(t)=\text{Pr}(X^{(0)}>t)=\left(\frac{{\tau }_c-t}{{\tau }_b}\wedge 1\right)\exp \{-({\lambda }_L+{\lambda }_D)t\}.\end{array}$$

With details relegated to the supplementary materials, we can reduce the null variance to ${\sigma }_0^2(\tau )=q^{-1}{(1-q)}^{-1}{\zeta }_0^2(\tau )$, where $q=n_1/n$ is the proportion of patients randomized to the treatment and(8) $$\begin{array}{c}{\zeta }_0^2(\tau )=E({\int }_0^{\tau }[\exp \{-(\lambda +{\lambda }_D)t\}\\ \left\{{\tilde{\delta }}^{(0)}I({\tilde{X}}^{(0)}\le t){\tilde{\pi }}^{(0)}{({\tilde{X}}^{(0)})}^{-1}-\lambda {\int }_0^{{\tilde{X}}^{(0)}\wedge t}{\tilde{\pi }}^{(0)}{(s)}^{-1}\mathrm{d}s\right\}\\ \hspace{1em}+\exp (-{\lambda }_{D}t)\left\{1-\exp (-\lambda t)\right\}\\ \{{\delta }^{(0)}I(X^{(0)}\le t){\pi }^{(0)}{(X^{(0)})}^{-1}\\ -{\lambda }_D{\int }_0^{X^{(0)}\wedge t}{\pi }^{(0)}{(s)}^{-1}\times \mathrm{d}s\left\}\right]\mathrm{d}t{)}^2.\end{array}$$(8)

The expectand on the right hand side of (8) can be computed via numerical integration, as the integrand consists only of analytic functions. The same strategy can be applied to the outer expectation using closed-form density functions of the data. But that is less convenient than the Monte Carlo approach—generate a large sample of $({\tilde{X}}^{(0)},{\tilde{\delta }}^{(0)},X^{(0)},{\delta }^{(0)})$ and approximate the expectation by the empirical average. Our numerical experiment shows that a Monte Carlo sample size of $N=10,000$ combined with M = 500 knots for numerical integration mostly controls the error within $\pm 1\%$. That level of precision is more than needed for subsequent use. Moreover, a single Monte Carlo sample can be used to evaluate ${\zeta }_0^2(\tau )$ for multiple τ (with a variable upper limit for the integral). In R, it takes only about 0.5 sec on a Dell XPS desktop with Intel Core i9-10900X (3.70 GHz) to evaluate a typical instance of ${\zeta }_0^2(\tau )$ for a whole spectrum of τ.

2.4 Practical Implementations

So far we have derived the “signal” $\mu (\tau )$ and “noise” ${\zeta }_0^2(\tau )$ needed to calculate the power and sample size. To highlight their dependence on the underlying parameters, we switch to the notation $\mu (\tau ;{\theta }_D,{\theta }_H,{\lambda }_D,{\lambda }_H,\kappa )$ and ${\zeta }_0^2(\tau ;{\lambda }_D,{\lambda }_H,\kappa ,{\tau }_b,{\tau }_c,{\lambda }_L)$, respectively. Thus, (6) and (7) become(9) $$\begin{array}{c}\text{Power:}\mathrm{ }1-\beta \\ =\Phi \left\{\frac{\sqrt{nq(1-q)}\mu (\tau ;{\theta }_D,{\theta }_H,{\lambda }_D,{\lambda }_H,\kappa )}{{\zeta }_0(\tau ;{\lambda }_D,{\lambda }_H,\kappa ,{\tau }_b,{\tau }_c,{\lambda }_L)}-z_{1-\alpha /2}\right\}\end{array}$$(9) and(10) $$\text{Sample}\mathrm{ }\text{size:}\mathrm{ }n=\frac{{\zeta }_0^2(\tau ;{\lambda }_D,{\lambda }_H,\kappa ,{\tau }_b,{\tau }_c,{\lambda }_L){(z_{1-\alpha /2}+z_{1-\beta })}^2}{q(1-q){\mu }^2(\tau ;{\theta }_D,{\theta }_H,{\lambda }_D,{\lambda }_H,\kappa )},$$(10) respectively. Usually patients are randomized 1:1 into the treatment and control groups, in which case q = 0.5. Nevertheless, the formulas can also handle unequal allocations should that be desirable in practice. The evaluation of ${\zeta }_0^2(\tau ;{\lambda }_D,{\lambda }_H,\kappa ,{\tau }_b,{\tau }_c,{\lambda }_L)$ is the rate-limiting step. But since it is independent of $({\theta }_D,{\theta }_H)$, we only need a one-off computation to get the power and sample size for different HRs.

Much hassle is saved by using our R functions rmtpower() and rmtss(), which carry out the calculations in (9) and (10), respectively. The rmtpower() function accepts vector-valued n, τ, $\exp ({\theta }_D)$, and $\exp ({\theta }_H)$ (the last two of which must be commensurate), along with scaler-valued ${\lambda }_D,{\lambda }_H,\kappa ,{\tau }_b,{\tau }_c,{\lambda }_L$, and α, and returns a three-dimensional array of power for every combination of the specified sample size, restricting time, and HRs. The rmtss() function works in a similar way except that a vector-valued $1-\beta$ should be supplied in place of n and a three-dimensional array of sample size will be returned for every combination of the specified target power, restricting time, and HRs. In both functions, the user can choose whether to apply the same to the RMST and RMEST, whose effect sizes are given in (5) and whose variances are obtained as a byproduct of the computation in Section 2.3. More details on this can be found in the supplementary materials.

Flexible use of the power and sample size programs can provide other useful information about trial design. For example, given a plausible estimate of patient accrual rate ρ, the investigator may need to know how long the accrual should last in order to have enough patients to achieve the target power. One can answer that question by plotting n as a function of ${\tau }_b\in [0,{\tau }_c]$, using rmtss(), with other parameters held constant, and then intersecting the curve with the straight line $n=\rho {\tau }_b$. The x-coordinate of the intersection point is the smallest τ_b needed to achieve the target power and the y-coordinate is the corresponding sample size. Nonintersecting curves mean that the target power is unattainable under the specified τ_c and ρ. To attain it, the investigator must either lengthen the trial or hasten recruitment.

Finally, the baseline parameters $({\lambda }_D,{\lambda }_H,\kappa )$ for the outcome distribution can be estimated from an external “pilot” study, that is, one conducted earlier with a similar patient population. We can do this using the simple estimators proposed in Mao et al. (2021). Specifically, let $({\tilde{X}}_i^{(0)},{\tilde{\delta }}_i^{(0)},X_i^{(0)},{\delta }_i^{(0)})$ $(i=1,\dots ,N_0)$ denote an N₀-random sample of $({\tilde{X}}^{(0)},{\tilde{\delta }}^{(0)},X^{(0)},{\delta }^{(0)})$ to be used as “pilot” data. We can then estimate λ_D , κ, and λ_H by$${\widehat{\lambda }}_D=\frac{\sum _{i=1}^{N_0}{\delta }_i^{(0)}}{\sum _{i=1}^{N_0}{X}_i^{(0)}},\hspace{1em}\widehat{\kappa }=\log (1-{\widehat{\lambda }}_{H\#}/\widehat{\lambda })/\log ({\widehat{\lambda }}_D/\widehat{\lambda }),$$and$${\widehat{\lambda }}_H={\widehat{\lambda }}_{H\#}{}^{1/\widehat{\kappa }}{\widehat{\lambda }}^{1-1/\widehat{\kappa }},$$ respectively, where$$\widehat{\lambda }=\frac{\sum _{i=1}^{N_0}{\tilde{\delta }}_i^{(0)}}{\sum _{i=1}^{N_0}{\tilde{X}}_i^{(0)}}\hspace{1em}\mathrm{ }\text{and}\hspace{1em}{\widehat{\lambda }}_{H\#}=\frac{\sum _{i=1}^{N_0}I({\tilde{X}}^{(0)}<X^{(0)})}{\sum _{i=1}^{N_0}{\tilde{X}}_i^{(0)}}.$$

These estimators are justified in sec. 3.2 of Mao et al. (2021) and implemented in the R function gumbel.est().

3 Simulation Studies

We used simulations to assess the accuracy of the power and sample size formulas. We first focused on the power function in (9) for the RMT-IF in comparison with those for the RMST and RMEST given in Section S1.2 of the supplementary materials. Let ${\tau }_b=3,{\tau }_c=4$, and ${\lambda }_L=0.01$. For the outcome model (1), we set ${\lambda }_D=0.05,0.4$, ${\lambda }_H=0.8$, and $\kappa =1.5$ (giving rise to Kendall’s rank correlation 33.3%). With ${\lambda }_D=0.05$, the observed mortality rate is about 10% and nonfatal event rate about 80%; with ${\lambda }_D=0.4$, the observed mortality rate is about 55% and nonfatal event rate about 65%.

Fixing $\exp ({\theta }_D)=0.8$, we plot the theoretical power curves of the RMT-IF, RMST, and RMEST as functions of $\exp ({\theta }_H)$ in Figure 1 under sample size n = 2000, Type I error rate $\alpha =0.05$, and treatment-to-control randomization ratio 1:1 (q = 0.5). Whether ${\lambda }_D=0.05$ or 0.4, the power curve of the RMST stays flat with regard to $\exp ({\theta }_H)$ as the test does not draw on the nonfatal event at all. When mortality rate is low, the RMST is easily outperformed by both RMT-IF and RMEST as soon as the nonfatal event takes on a moderate treatment effect. When mortality rate is high, death alone is sufficient to power the RMST. On the other hand, the RMEST can be severely underpowered when the treatment has minimal effect on the nonfatal event. This is not surprising as the composite first event is still dominated by the nonfatal type, which can mask the beneficial effect on overall survival. By contrast, the RMT-IF, with overall survival prioritized, stays largely robust to the varying effect size on the nonfatal component. The RMT-IF’s general power advantage over the two competitors also owes to its using more events per patient.

Fig. 1 Theoretical and empirical power of the RMT-IF, RMST, and RMEST under model (1) with ${\lambda }_H=0.8$, $\exp ({\theta }_D)=0.8$ and n = 2000 as a function of $\exp ({\theta }_H)$. Dashed lines, theoretical power curve; round dots, empirical power based on 5000 replicates.

These comparisons are based on theoretical formulas, albeit numerically evaluated. We also computed the empirical rejection rates of the three tests under different $\exp (\theta )$ based on 5000 replicates of simulated data, and then overlaid them on the theoretical curves in Figure 1. We can see that the empirical values roughly coincide with the theoretical predictions, showing the validity of the asymptotic formulas in finite samples.

Next, we evaluated the sample size formula (10) for the RMT-IF. We used the same set-up as in the first simulations except that we set a moderate death hazard ${\lambda }_D=0.2$ and varied κ over 1.0, 2.0, and 3.0 (corresponding to between-component Kendall’s rank correlations $0,50\%$, and 66.7%, respectively, Oakes 1989). The observed mortality rate under this set-up is about 35% (thus, censoring rate 65%) and nonfatal event rate about 75%. We considered randomization ratios 1:1 (q = 0.5) and 2:1 (q = 0.67). For each scenario, the sample size n required to achieve power $1-\beta =0.8$ or 0.9 is computed by (10), and then 5000 samples of size n are generated to compute the empirical power of the RMT-IF test. Table 1 shows strong agreement between the empirical and target power in all scenarios, especially for large n, indicating the overall accuracy of the formula. We also considered set-ups with lower and higher censoring rates in the supplementary materials and the results are comparable to those in Table 1.

Table 1 Simulation results on the empirical power of sample sizes computed from (10) under model (1) with HRs $\exp ({\theta }_D)$ and $\exp ({\theta }_H)$.

			Power $=0.8$				Power $=0.9$
			$\tau =1.5$		$\tau =3.0$		$\tau =1.5$		$\tau =3.0$
q	κ	$e^{{\theta }_D},e^{{\theta }_H}$	n	EP	n	EP	n	EP	n	EP
0.50	1.0	0.6, 0.6	214	0.82	157	0.82	287	0.92	210	0.91
		0.7, 0.8	695	0.81	469	0.81	931	0.88	627	0.90
		0.8, 0.8	998	0.79	781	0.80	1335	0.91	1046	0.90
	2.0	0.6, 0.6	278	0.82	197	0.82	372	0.92	264	0.90
		0.7, 0.8	989	0.81	635	0.79	1324	0.91	850	0.90
		0.8, 0.8	1264	0.79	960	0.80	1693	0.90	1285	0.89
	3.0	0.6, 0.6	295	0.82	209	0.82	395	0.93	279	0.92
		0.7, 0.8	1084	0.82	687	0.80	1452	0.92	919	0.89
		0.8, 0.8	1341	0.80	1013	0.80	1795	0.91	1355	0.91
0.67	1.0	0.6, 0.6	241	0.82	177	0.82	323	0.90	237	0.90
		0.7, 0.8	782	0.81	527	0.81	1047	0.91	706	0.91
		0.8, 0.8	1122	0.80	879	0.81	1502	0.90	1177	0.92
	2.0	0.6, 0.6	312	0.82	222	0.83	418	0.92	297	0.92
		0.7, 0.8	1113	0.82	714	0.81	1490	0.90	956	0.91
		0.8, 0.8	1422	0.80	1080	0.80	1904	0.91	1446	0.90
	3.0	0.6, 0.6	332	0.83	235	0.83	445	0.93	314	0.92
		0.7, 0.8	1220	0.80	773	0.82	1633	0.91	1034	0.90
		0.8, 0.8	1509	0.79	1139	0.81	2020	0.90	1525	0.91

Note: EP, empirical power; each scenario is based on 5000 replicates.

Finally, we turned to the problem of finding the minimum length of accrual τ_b given a fixed accrual rate ρ and total length of follow-up τ_c . We used the same set-up as that of Table 1 except that we considered a variable τ_b while fixing $\exp ({\theta }_D)=\exp ({\theta }_H)=0.8,\kappa =1.5$, and τ = 3 for simplicity. For target power $1-\beta =0.8$ and 0.9, we plotted the sample size n computed by (10) as a function of τ_b in Figure 2. As τ_b approaches τ_c , the administrative censoring time becomes stochastically shorter, driving up the censoring rate and thus the total number of patients needed. Then, as suggested in Section 2.4, we found the intersection points of the sample size curve with $n=\rho {\tau }_b$ for $\rho =400,800$ per unit time, as shown in Figure 2. Based on the coordinates $({\tau }_b,n)$ of each intersection point, we simulated 5000 samples to compute the empirical power of the RMT-IF test. As shown in Table 2, the empirical power is approximately equal to the target power under which τ_b and n are derived.

Fig. 2 Finding the minimum τ_b and n under fixed accrual rate ρ and total trial duration ${\tau }_c=4$. Dashed line, ρ = 400; dotted line, ρ = 800.

Table 2 Simulation results on the empirical power of n and τ_b computed from Figure 2.

		Power $=0.8$			Power $=0.9$
q	ρ	τb	n	EP	τb	n	EP
0.50	400	2.15	860	0.81	3.02	1208	0.90
	800	1.07	856	0.80	1.44	1152	0.89
0.67	400	2.44	976	0.79	3.79	1516	0.89
	800	1.21	968	0.81	1.62	1296	0.90

NOTE: see note to Table 1.

4 A Real Example

Lymph node-positive breast cancer occurs when cancerous cells spread from the primary location to nearby lymph nodes. Left untreated, the malignancy is likely to metastasize to other parts of the body, threatening the patient’s life. The standard treatment of node-positive breast cancer is chemotherapy, sometimes with hormonal supplements. In the late 1980s, the German Breast Cancer (GBC) study group evaluated the effectiveness of different cycles of chemotherapy and of hormonal treatment with tamoxifen on 720 node-positive breast cancer patients in a 2 × 2 design (Schmoor et al. 1996). The study showed moderate benefits of longer duration of chemotherapy and hormonal treatment on the patient’s relapse-free survival time. On the other hand, both chemotherapy and hormonal treatment are associated with severe side effects such as hair loss, nausea, diarrhea, and heart problems.

Suppose that an investigational treatment is to be tested against the standard care of node-positive breast cancer patients. Patient outcomes can be analyzed using the RMT-IF with cancer relapse as the nonfatal event along with death. We can use our power and sample size formulas to design the new trial, with baseline parameters estimated from the GBC study. To do so, we consider the same subset of 686 patients analyzed in Sauerbrei and Royston (1999). With median and maximum lengths of follow-up of 3.8 and 7.2 years, respectively, the observed mortality rate is about 25% and relapse rate about 38%. More formally, using the procedures described in Section 2.4, we find that ${\widehat{\lambda }}_D=0.069$ year^{– 1}, ${\widehat{\lambda }}_H=0.131$ year^{– 1}, and $\widehat{\kappa }=3.9$ (the last of which means that Kendall’s rank correlation between death and relapse is $1-{\widehat{\kappa }}^{-1}\approx 75\%$). Collectively they lead to a relapse-free survival hazard of $\widehat{\lambda }=0.134$ year^{– 1}. To check the validity of (1) as a working model, especially the assumption of constant hazard rates, we plot the Nelsen–Aalen estimates of the component-wise and overall cumulative hazards in Figure S1 of the supplementary materials. The nonparametric curves appear well approximated by the model-based straight lines, indicating a reasonable fit of the Gumbel–Hougaard model. Based on the hazard rate estimates, the 5-year overall and relapse-free survival rates are 71% and 51%, respectively; the 5-year RMST and RMEST are 4.2 and 3.6 years, respectively. These estimates are roughly in line with common epidemiologic evidence on node-positive breast cancer patients at large (Chu et al. 1996).

To size the new trial for a 5-year RMT-IF test with level 0.05, suppose that the trial will last ${\tau }_c=7$ years in total, with the first ${\tau }_b=3$ years devoted to patient enrollment. With randomization ratio 1:1 and 2:1, the sample size as a function of the HRs $\exp ({\theta }_D)$ and $\exp ({\theta }_H)$ is computed by (10) and displayed as surface plots in Figure 3. In particular, under randomization ratio 1:1 and HRs $\{\exp ({\theta }_D),\exp ({\theta }_H)\}=(0.6,0.6)$ and $(0.9,0.9)$, we need n = 471 and 8450 to achieve 80% power, and n = 630 and 11,312 to achieve 90% power, respectively.

Fig. 3 Sample size needed to power a 5-year RMT-IF test as a function of the HRs on death and relapse with baseline parameters estimated from the GBC data.

To explore the relationship between the rate and duration of patient accrual, we compute τ_b as a function of ρ under the fixed ${\tau }_c=7$ years using the method described in Section 2.4. For varying effect sizes, the corresponding τ_b and associated sample size n are tabulated in Table 3. Given a plausible estimate of the accrual rate, the investigator can use this table to look up the length of accrual needed to reach the target power.

Table 3 Duration of accrual τ_b (years) and sample size n under different accrual rate ρ (year^{– 1}) with baseline parameters estimated from the GBC data.

		Power $=0.8$				Power $=0.9$
		q = 0.50		q = 0.67		q = 0.50		q = 0.67
$e^{{\theta }_D},e^{{\theta }_H}$	ρ	τb	n	τb	n	τb	n	τb	n
0.80, 0.80	500	4.10	2050	4.67	2335	5.94	2970	–	–
	1000	2.04	2040	2.29	2290	2.73	2730	3.07	3070
	2000	1.02	2040	1.15	2300	1.37	2740	1.54	3080
0.80, 0.85	800	3.90	3120	4.42	3536	5.46	4368	–	–
	1200	2.59	3108	2.91	3492	3.46	4152	3.91	4692
	2400	1.30	3120	1.46	3504	1.73	4152	1.95	4680
0.80, 0.90	1300	4.08	5304	4.64	6032	5.87	7631	–	–
	2100	2.51	5271	2.82	5922	3.36	7056	3.79	7959
	3000	1.76	5280	1.98	5940	2.35	7050	2.64	7920

NOTE: indicates a scenario with target power unattainable under ${\tau }_c=$ 7 years.

Finally, it should be noted that in some studies a new line of treatment is initiated for patients upon relapse of cancer. The treatment switch as an intercurrent event (ICH 2020) is better accounted for as a competing risk (Conner and Trinquart 2021). In such cases, the RMT-IF can be estimated using the cumulative incidence functions of the component events as suggested in sec. 3.2 of Mao (2021). Meanwhile, sample size calculation needs to be redesigned accordingly.

5 Concluding Remarks

Under a bivariate copula model, we have derived power and sample size formulas for the RMT-IF, a new way to analyze prioritized composite endpoints of death and nonfatal events. Our approach allows for between-component correlation (quantifiable through data on previous studies) as well as component-specific treatment effects. The simplicity of the developed routines owes to the analytic expression of the estimand and the simplified form of the null variance, which together make it possible to compute the noncentrality parameter efficiently for varying restricting times. Investigators can now use the R-package rmt (https://cran.r-project.org/package=rmt) to design new trials based on the RMT-IF analysis.

The Gumbel–Hougaard copula in (1) is particularly convenient working model for the RMT-IF as it leads to a simple expression for the estimand. Other copulas may not be nearly as nice in that regard (Nelsen 2007). However, the model is also highly restrictive, especially concerning the time-constancy of hazard rates. In practice, it is recommended that the investigator perform model checking using data from similar studies, like we did in Section 4, before using the formulas to design new trials. In the future, it will be of interest to relax the working model to allow for time-dependent and non-proportional hazard structures, similarly to Yung and Liu (2020) in the univariate case, and to try out other joint models, for example, Clayton copulas (Clayton 1978), for sensitivity analysis.

We have considered only the case with a single occurrence of nonfatal event besides death. In some trials, especially those in chronic diseases, one patient can experience an adverse event repeatedly (Rogers et al. 2014; Mao and Kim 2021). In the recent INfluenza Vaccine to Effectively Stop CardioThoracic Events and Decompensated Heart Failure (INVESTED) trial (NCT02787044), for example, a high-risk cardiovascular patient could be hospitalized up to eight times per influenza season (Vardeny et al. 2021). Such recurrent events can be analyzed by the RMT-IF equally well, only with a slightly modified definition of “a better outcome” (i.e., being alive or having had fewer recurrent events). However, joint models like (1) are ill suited for recurrent events because they rarely hold when the event times are inherently ordered (Mao et al. 2022). Consequently, a new framework is needed for trial design.

We have focused on the power and sample size calculations under a fixed design. For ethical as well as financial considerations, real trials often prefer, or even mandate, interim analysis in the course of follow-up to assess whether the trial should be stopped early (Jennison and Turnbull 1999). With a model-free metric like the RMST and RMT-IF, group sequential analysis is complicated by the lack of an explicit correlation structure between the interim test statistics. Recently, Luo et al. (2019) obtained simple results for group sequential RMST tests assuming piecewise exponential distributions for the event and censoring times. Lu and Tian (2021) detailed a fully nonparametric approach which determines the rejection regions empirically using the data available at each interim look. The design and implementation of group sequential RMT-IF analysis remain an open problem.

Supplementary Materials

Supplementary Materials online include technical details and additional simulation results referenced in Sections 2 and 3, as well as the R-code and data to reproduce the analysis in Section 4. The methods proposed are implemented in the R-package rmt openly available on the Comprehensive R Archive Network (https://cran.r-project.org/package=rmt).

Additional information

Funding

This research was supported by the National Institutes of Health grant R01HL149875.