Group-sequential logrank methods for trial designs using bivariate non-competing event-time outcomes

Abstract

We discuss the multivariate (2L-variate) correlation structure and the asymptotic distribution for the group-sequential weighted logrank statistics formulated when monitoring two correlated event-time outcomes in clinical trials. The asymptotic distribution and the variance–covariance for the 2L-variate weighted logrank statistic are derived as available in various group-sequential trial designs. These methods are used to determine a group-sequential testing procedure based on calendar times or information fractions. We apply the theoretical results to a group-sequential method for monitoring a clinical trial with early stopping for efficacy when the trial is designed to evaluate the joint effect on two correlated event-time outcomes. We illustrate the method with application to a clinical trial and describe how to calculate the required sample sizes and numbers of events.

Appendices

Appendix

Proof of Theorem 1

Let ${M}_{ik}^{({\ell })}(t)={N}_{ik}^{({\ell })}(t)-\int _{0}^{t}{Y}_{ik}^{({\ell })}(x)d{\Lambda }_{g_i k}(x)$ and $\{{\mathscr {F}}_{k,t}^{({\ell })}:t\ge 0 \}$ be a standard filtration generated from the history through time t for the kth outcome and the ${\ell }$th analysis (${\mathscr {F}}_{k,t}^{({\ell })}$ is the smallest $\sigma$-algebra generated by $\{{N}_{ik}^{({\ell })}(x),{N}_{ik}^{C({\ell })}(x): 0\le x\le t,i=1,\cdots ,n_{\ell }\}$, where ${N}_{ik}^{C({\ell })}(t)={\mathbb {1}}\{T_{ik}^{({\ell })}\le t,$${\Delta }_{ik}^{({\ell })}=0\}$ is a censoring counting process). As is well-known, ${M}_{ik}^{({\ell })}(t)$ has the ${\mathscr {F}}_{k,t}^{({\ell })}$-martingale property. We discuss the asymptotic behavior using the decomposition of the weighted logrank process $U_k^{({\ell })}(t)=\hat{m}_k^{({\ell })}(t)+n_{{\ell }}^{-\frac{1}{2}}\mathcal {M}_{k}^{({\ell })}(t)$ from the definition of $U_k^{({\ell })}$, where

$$\begin{aligned}&\mathcal {M}_{k}^{({\ell })}(t) =\int _{0}^{t}\hat{H}_k^{({\ell })}(x)n_{{\ell }}^{\frac{1}{2}}{\textstyle \sum _{i=1}^{n_L}\mathrm{d}{\widetilde{M}}_{ik}^{({\ell })}(x)} =\int _{0}^{t}\hat{H}_k^{({\ell })}(x)n_{{\ell }}^{\frac{1}{2}}\mathrm{d}{\widetilde{M}}_{\bullet k}^{({\ell })}(x),\\&\mathrm{d}{\widetilde{M}}_{\bullet k}^{({\ell })}(x) = \frac{\mathrm{d}\overline{M}_{2k}^{({\ell })}(x)}{\overline{{Y}}_{2k}^{({\ell })}(x)}-\frac{\mathrm{d}\overline{M}_{1k}^{({\ell })}(x)}{\overline{{Y}}_{1k}^{({\ell })}(x)},\,\,\,\, \mathrm{d}\overline{M}_{jk}^{({\ell })}(x)=\textstyle \sum _{i=1}^{n_{\ell }}{\mathbb {1}}\{g_i=j\}\mathrm{d}{{M}}_{ik}^{({\ell })}(x), \\&\mathrm{d}{\widetilde{M}}_{ik}^{({\ell })}(x) ={\mathbb {1}}{\{i\le n_{{\ell }}\}} \left\{ \frac{{\mathbb {1}}{\{g_{i}=2\}}}{\overline{{Y}}_{2k}^{({\ell })}(x)}-\frac{{\mathbb {1}}{\{g_{i}=1\}}}{\overline{{Y}}_{1k}^{({\ell })}(x)} \right\} \mathrm{d}{{M}}_{ik}^{({\ell })}(x), \end{aligned}$$

and $\mathcal {M}_{k}^{({\ell })}(t)$ is ${\mathscr {F}}_{k,t}^{({\ell })}$-martingale because $\hat{H}_k^{({\ell })}(t)$ is ${\mathscr {F}}_{k,t}^{({\ell })}$-predictable.

Let $\hat{\varvec{Z}}^*=(\hat{Z}^*_{1}({\tau }_1),\ldots ,\hat{Z}^*_1({\tau }_L),\hat{Z}^*_{2}({\tau }_1),\ldots ,\hat{Z}^*_{2}({\tau }_L))^\mathrm{T}$ and let $\hat{Z}^*_k(\tau _{\ell })$ be $\hat{Z}_k(\tau _{\ell })$ whose denominator is replaced by the limit version,

$$\begin{aligned} \hat{Z}^*_k(\tau _{\ell }) =n_{\ell }^{\frac{1}{2}}\displaystyle \frac{U_k^{({\ell })}(\tau _{\ell })}{\sqrt{{V}_{kk}^{0({\ell })}(\tau _{\ell })}} =n_{\ell }^{\frac{1}{2}}\hat{\mu }_{k{\ell }}+\xi _{k{\ell }}\mathcal {M}_{k}^{({\ell })}(\tau _{\ell }), \end{aligned}$$

where we write $\xi _{k{\ell }}=1/\sqrt{{V}_{kk}^{0({\ell })}(\tau _{{\ell }})}$ for simplicity. The distribution of $\hat{\varvec{Z}}-D_{\scriptstyle {\varvec{n}}}\hat{\varvec{\mu }}$ is asymptotically equivalent to

$$\begin{aligned} \hat{\varvec{Z}}^* - D_{\scriptstyle {\varvec{n}}}\hat{\varvec{\mu }} =\left( \xi _{11}\mathcal {M}_{1}^{(1)}(\tau _1),\ldots ,\xi _{1L}\mathcal {M}_{1}^{(L)}(\tau _L), \xi _{21}\mathcal {M}_{2}^{(1)}(\tau _1),\ldots ,\xi _{2L}\mathcal {M}_{2}^{(L)}(\tau _L) \right) ^\mathrm{T} \end{aligned}$$

because the dominated convergence theorem works by the convergence of $\hat{V}_{kk}^{0({\ell })}(\tau _{\ell })$${\mathop {\rightarrow }\limits ^{P}} {V}_{kk}^{0({\ell })}(\tau _{\ell })$ uniformly on ${\ell }=1,\ldots ,L$ as $n_L\ge \cdots \ge n_1\rightarrow \infty$. We find it necessary to study the covariance of $\mathcal {M}_{k}^{({\ell })}$’s for characterizing the distribution of $\hat{\varvec{Z}}^*-D_{\scriptstyle {\varvec{n}}}\hat{\varvec{\mu }}$.

In the proof hereafter, it is sufficient to consider the case of $L=2$. As a function related to the characteristic function of ${\mathcal {M}}_{k}^{({\ell })}(t)$, define

$$\begin{aligned} G_{k}^{({\ell })}(t)=\exp \left( {\varvec{i} z_{k{\ell }}{\mathcal {M}}_{k}^{({\ell })}(t)+{\textstyle \frac{z_{k{\ell }}^2}{2}}{\langle }{\mathcal {M}}_{k}^{({\ell })},{\mathcal {M}}_{k}^{({\ell })}{\rangle }(t)} \right) \end{aligned}$$

for a real non-zero $z_{k{\ell }}$ and $\varvec{i}=\sqrt{-1}$, where ${\langle }m_1,m_2{\rangle }$ denotes a predictable covariance process for two martingales $m_1$ and $m_2$. In this case we have

$$\begin{aligned} {\langle }{\mathcal {M}}_{k}^{({\ell })},{\mathcal {M}}_{k}^{({\ell }')}{\rangle }(t) = n_{{\ell }}^{\frac{1}{2}} n_{{\ell }'}^{\frac{1}{2}} \int _{0}^{t}\hat{H}_k^{({\ell })}(x)\hat{H}_k^{({\ell }')}(x) \left\{ \frac{\mathrm{d}{\Lambda }_{1k}(x)}{\overline{Y}_{1k}^{({\ell }\vee {\ell }')}(x)} +\frac{\mathrm{d}{\Lambda }_{2k}(x)}{\overline{Y}_{2k}^{({\ell }\vee {\ell }')}(x)} \right\} , \end{aligned}$$

following the standard martingale theory of survival analysis (see Fleming and Harrington (1991)). The consistency of $\hat{S}^{({\ell })}_{jk}$, the Glivenko-Cantelli theorem, and Conditions 1 and 3 imply $\sup _{0\le x\le \tau _{{\ell }}}$$|\hat{H}_k^{({\ell })}(x)-{H}_k^{({\ell })}(x)|{\mathop {\rightarrow }\limits ^{P}}0$ and

$$\begin{aligned} \sup _{0\le x\le \tau _{{\ell }}} \left| {\hat{H}_k^{({\ell })}(x)}\bigl /n_{{\ell }}^{-1}\bar{Y}_{jk}^{({\ell })}(x) -h_{jk}^{({\ell })}(x) \right| {\mathop {\rightarrow }\limits ^{P}}0\quad \mathrm{as}\quad n_{{\ell }}\rightarrow \infty , \end{aligned}$$

(9)

where

$$\begin{aligned} h_{jk}^{({\ell })}(x) =\frac{{H}_{k}^{({\ell })}(x)}{a_{j{\ell }} y_{jk}^{({\ell })}(x)} =W_k^{({\ell })}(x) \frac{a_{j'{\ell }} S_{j'k}(x_-)}{S_{\bullet k}^{({\ell })}(x_-)},\,\,\,\,j'=3-j, \end{aligned}$$

and note that $0\le {H}_k^{({\ell })}(x)<\infty$ for $x\in [0,\tau _{{\ell }}]$, ${H}_k^{({\ell })}(x)=0$ for $\tau _{{\ell }}<x$ and $0\le h_{jk}^{({\ell })}(x)<\infty$ for all x. The univariate asymptotic result provides $E(\mathrm{e}^{\scriptstyle {\varvec{i}} z_{k{\ell }}{\mathcal {M}}_{k}^{({\ell })}(t)})\rightarrow \exp (-\frac{z_{k{\ell }}^2}{2}V_{kk}(t,t\mid \tau _{{\ell }},\tau _{{\ell }}))$ as $n_{{\ell }}\rightarrow \infty$, which corresponds to the following convergences,

$$\begin{aligned} E(G_{k}^{({\ell })}(t))\rightarrow 1 \quad \mathrm{and}\quad {\langle }{\mathcal {M}}_{k}^{({\ell })},{\mathcal {M}}_{k}^{({\ell }')}{\rangle }(t){\mathop {\rightarrow }\limits ^{P}} V_{kk}(t,t\mid \tau _{{\ell }},\tau _{{\ell }'}) \end{aligned}$$

(Nishiyama 2011). For different $k\ne k'$, it is difficult to show joint normality with correlation between ${\mathcal {M}}_{k}^{({\ell })}$ and ${\mathcal {M}}_{k'}^{({\ell })}$ with standard martingale theory of counting processes (Fleming and Harrington 1991; Andersen et al. 1993). However, we overcome the challenge applying Ito’s formula. The discrete Ito’s formula (Jacod and Shiryaev 2003; Huang and Strawderman 2006) provides the decomposition of $G_{k}^{({\ell })}(t)$,

$$\begin{aligned} G_{k}^{({\ell })}(t)-1= & {} \displaystyle \sum _{j=1,2} \int _{0}^{t}G_{k}^{({\ell })}(x-) \widetilde{H}_{jk}^{\mathrm{a}({\ell })}(x) \mathrm{d}{\overline{M}}_{jk}^{({\ell })}(x) \nonumber \\&+\displaystyle \sum _{j=1,2} \int _{0}^{t}G_{k}^{({\ell })}(x-)\widetilde{H}_{jk}^{\underline{(}{\ell })}(x)\overline{Y}_{jk}^{({\ell })}(x)\mathrm{d}{{\Lambda }}_{jk}(x), \end{aligned}$$

(10)

where, with $\varvec{i}_1=-\varvec{i}$ and $\varvec{i}_2=\varvec{i}$,

$$\begin{aligned} \widetilde{H}_{jk}^{\mathrm{a}({\ell })}(x)= & {} \exp \left( \varvec{i}_j z_{k{\ell }} \textstyle \frac{\sqrt{n_{{\ell }}}\hat{H}_{k}^{({\ell })}(x)}{\overline{Y}_{jk}^{({\ell })}(x)}\right) -1,\\ \widetilde{H}_{jk}^{\mathrm {b}({\ell })}(x)= & {} \exp \left( \varvec{i}_j z_{k{\ell }} \textstyle \frac{\sqrt{n_{{\ell }}}\hat{H}_{k}^{({\ell })}(x)}{\overline{Y}_{jk}^{({\ell })}(x)}\right) -1-\varvec{i}_j z_{k{\ell }} \displaystyle \frac{\sqrt{n_{{\ell }}}\hat{H}_{k}^{({\ell })}(x)}{\overline{Y}_{jk}^{({\ell })}(x)}\\&+\frac{z_{k{\ell }}^2}{2}\left( \displaystyle \frac{\sqrt{n_{{\ell }}}\hat{H}_{k}^{({\ell })}(x)}{\overline{Y}_{jk}^{({\ell })}(x)}\right) ^2. \end{aligned}$$

The expectation of the right-hand side of (10) converges to zero as $n_{{\ell }}\rightarrow \infty$, because

$$\begin{aligned} \begin{array}{l} E\left( \textstyle \int _{0}^{t}G_{k}^{({\ell })}(x-)\widetilde{H}_{jk}^{\mathrm{a}({\ell })}(x)\mathrm{d}{\overline{M}}_{jk}^{({\ell })}(x)\right) =0\\ \mathrm{and}\quad \textstyle E\left( \int _{0}^{t}G_{k}^{({\ell })}(x-)\widetilde{H}_{jk}^{\mathrm {b}({\ell })}(x)\overline{Y}_{jk}^{({\ell })}(x)\mathrm{d}{{\Lambda }}_{jk}(x)\right) \rightarrow 0\\ \end{array} \end{aligned}$$

(11)

by the martingale property of $\overline{M}_{jk}^{({\ell })}$ and the Lindeberg condition, respectively. In fact, using the integrable martingale property of $G_{k}^{({\ell })}(x_{-})$ and the well-known inequality

$$\begin{aligned} \left| \exp (\varvec{i}c)-1-\varvec{i}c+\textstyle \frac{1}{2}c^{2} \right| \le {\mathbb {1}}\{|c|\le \varepsilon \}|c|^{3}+{\mathbb {1}}\{|c|>\varepsilon \}|c|^{2} \end{aligned}$$

for any real c, the latter result of (11) is obtained as

$$\begin{aligned}&E\left( \textstyle \int _0^t \left| G_{k}^{({\ell })}(x_{-})\widetilde{H}_{jk}^{\mathrm {b}({\ell })}(x) \right| \bar{Y}_{jk}^{({\ell })}(x)\mathrm{d}{{\Lambda }}_{jk}(x) \right) \\&\quad \le \exp \left( \textstyle \frac{z_{k{\ell }}^{2}}{2} {\langle }{\mathcal {M}}_{k}^{({\ell })},{\mathcal {M}}_{k}^{({\ell })}{\rangle }(t) \right) \times \\&\qquad \Biggl \{ E\left( \textstyle \int _0^t |c_{jk{\ell }}(x)|^{3} {\mathbb {1}}\{|c_{jk{\ell }}(x)|\le \varepsilon \} \bar{Y}_{jk}^{({\ell })}(x) \mathrm{d}{\Lambda }_{jk}(x) \right) \\&\qquad + E\left( \textstyle \int _0^t {|c_{jk{\ell }}(x)|^{2}{\mathbb {1}}\{|c_{jk{\ell }}(x)|>\varepsilon \}\bar{Y}_{jk}^{({\ell })}(x)} \mathrm{d}{\Lambda }_{jk}(x) \right) \Biggr \} {\rightarrow } 0 \end{aligned}$$

as $n_{{\ell }}\rightarrow \infty$, where $\varepsilon$ is an arbitrary positive number, $c_{jk{\ell }}(x)=z_{k{\ell }}\sqrt{n_{{\ell }}}\hat{H}_{k}^{({\ell })}(x)/\bar{Y}_{jk}^{({\ell })}(x)$ and we have $\sqrt{n_{{\ell }}}c_{jk{\ell }}(x){\mathop {\rightarrow }\limits ^{P}} {z_{k{\ell }}h_{jk}^{({\ell })}(x)}$ uniformly on $(0,\tau _{{\ell }}]$ from (9). Hence, we have

$$\begin{aligned} E((G_{1}^{({\ell })}(t)-1)(G_{2}^{({\ell }')}(s)-1))\rightarrow E(G_{1}^{({\ell })}(t)G_{2}^{({\ell }')}(s))-1 \end{aligned}$$

(12)

as $n_{{\ell }},n_{{\ell }'}\rightarrow \infty$ by the univariate results of $E(G_{k}^{({\ell })}(t))\rightarrow 1$, while using the formula (10) we can also find

$$\begin{aligned}&E((G_{1}^{({\ell })}(t)-1)(G_{2}^{({\ell }')}(s)-1)) \rightarrow \nonumber \\&\quad \displaystyle \sum _{j=1,2} \displaystyle \int _{0}^{t}\int _{0}^{s} E\left( G_{1}^{({\ell })}(x_-)G_{2}^{({\ell }')}(y_-) \widetilde{H}_{j1}^{\mathrm{a}({\ell })}(x)\widetilde{H}_{j2}^{\mathrm{a}({\ell }')}(y) \mathrm{d}{\overline{M}}_{j1}^{({\ell })}(x)\mathrm{d}{\overline{M}}_{j2}^{({\ell }')}(y) \right) \qquad \end{aligned}$$

(13)

as $n_{{\ell }},n_{{\ell }'}\rightarrow \infty$. Similarly to showing the latter result of (11), with asymptotic equality, we can replace the terms $\mathrm{e}^{\scriptstyle {\varvec{i}}_{j}(c_{j1{\ell }}(x)+c_{j2{\ell }'}(y))}$ and $\mathrm{e}^{\scriptstyle {\varvec{i}}_{j}c_{jk{\ell }}(\cdot )}$ included in (13) by

$$\begin{aligned} 1+\varvec{i}_{j}\{c_{j1{\ell }}(x)+c_{j2{\ell }'}(y)\}-\textstyle \frac{1}{2}\{c_{j1{\ell }}(x)+c_{j2{\ell }'}(y)\}^{2} \quad \mathrm{and}\quad 1+\varvec{i}_{j}c_{jk{\ell }}(\cdot )-\textstyle \frac{1}{2}c_{jk{\ell }}(\cdot )^{2}, \end{aligned}$$

respectively. In fact, we can show that

$$\begin{aligned} \tilde{H}_{j1}^{\mathrm{a}({\ell })}(x)\tilde{H}_{j2}^{\mathrm{a}({\ell }')}(y)= & {} \mathrm{e}^{\scriptstyle {\varvec{i}}_{j}(c_{j1{\ell }}(x)+c_{j2{\ell }'}(y))} -\mathrm{e}^{\scriptstyle {\varvec{i}}_{j} c_{j1{\ell }}(x)}-\mathrm{e}^{\scriptstyle {\varvec{i}}_{j} c_{j2{\ell }'}(y)}+1\\= & {} -c_{j1{\ell }}(x)c_{j2{\ell }'}(y)+o_{P}(1/\sqrt{n_{{\ell }}n_{{\ell }'}}) \end{aligned}$$

from the convergence result of $\sqrt{n_{{\ell }}}c_{jk{\ell }}(x)$. Hence, we have

$$\begin{aligned} \sqrt{n_{{\ell }}n_{{\ell }'}} \tilde{H}_{j1}^{\mathrm{a}({\ell })}(x)\tilde{H}_{j2}^{\mathrm{a}({\ell }')}(y) {\mathop {\rightarrow }\limits ^{P}} -z_{1{\ell }}z_{2{\ell }'}{h_{j1}^{({\ell })}(x)h_{j2}^{({\ell }')}(y)} \end{aligned}$$

(14)

as $n_{{\ell }},n_{{\ell }'}\rightarrow \infty$, so that we can apply this result to (13). Also, similar to Prentice and Cai (1992) and Sugimoto et al. (2013, 2017), we can show

$$\begin{aligned} \frac{1}{\hat{a}_{j{\ell }\wedge {\ell }'}n_{{\ell }\wedge {\ell }'}} E\left( {\iint } \mathrm{d}\overline{M}_{j1}^{({\ell })}(x) \mathrm{d}\overline{M}_{j2}^{({\ell }')}(y) \right)= & {} E\left( {\iint } {\mathrm{d}}M_{i1}^{({\ell })}(x) {\mathrm{d}}M_{i2}^{({\ell }')}(y) \mid g_i=j\right) \\= & {} {\iint } C_{{\ell }\wedge {\ell }'}(x\wedge y)A_{j}({\mathrm{d}}x,{\mathrm{d}}y). \end{aligned}$$

For simplicity, let $\phi (t,s)=E(G_{1}^{({\ell })}(t)G_{2}^{({\ell }')}(s))$. From (12), (13), (14), $\hat{\gamma }_{{\ell }}{\mathop {\rightarrow }\limits ^{P}}{\gamma }_{{\ell }}$, $\hat{a}_{j{\ell }}{\mathop {\rightarrow }\limits ^{P}}{a}_{j{\ell }}$ (Conditions 1–2) and the dominated convergence theorem, we have the integral equation for $\phi (t,s)$ under $n_{{\ell }},n_{{\ell }'}\rightarrow \infty$,

$$\begin{aligned}&\phi (t,s)-1=-z_{1{\ell }}z_{2{\ell }'}\sqrt{\displaystyle \frac{\gamma _{{\ell }\wedge {\ell }'}}{\gamma _{{\ell }\vee {\ell }'}}}\nonumber \\&\quad \times \int _{0}^{t}\int _{0}^{s} \phi (x_-,y_-) \sum _{j=1}^{2} a_{j{\ell }\wedge {\ell }'} h_{j1}^{({\ell })}(x) h_{j2}^{({\ell }')}(y) C_{{\ell }\wedge {\ell }'}(x\wedge y)A_j({\mathrm{d}}x,{\mathrm{d}}y). \end{aligned}$$

(15)

Similarly to bivariate survival function (Dabrowska 1988), the two-dimensional Volterra integral equation

$$\begin{aligned} \phi (t,s)=1 +\textstyle \int _{0}^{t}\int _{0}^{s}\phi (x_-,y_-)b_{12}({\mathrm{d}}x,{\mathrm{d}}y) \quad \mathrm{with}\quad \phi (t,0)=\phi (0,s)=1 \end{aligned}$$

is solved as $\phi (t,s)=\exp [\int _{0}^{t}\int _{0}^{s}\{b_{12}({\mathrm{d}}x,{\mathrm{d}}y)-b_{1}({\mathrm{d}}x,y)b_{2}(x,{\mathrm{d}}y)\}]$, where

$$\begin{aligned} b_{1}({\mathrm{d}}x,y)=\phi ({\mathrm{d}}x,y)/\phi (x_-,y_-) \quad \mathrm{and}\quad b_{2}(x,{\mathrm{d}}y)=\phi (x,{\mathrm{d}}y)/\phi (x_-,y_-). \end{aligned}$$

However, note that it is difficult to obtain $b_{k}(x,y)$, $k=1,2$ by directly differentiating (15) because of including the expectation of non-differentiable $M_{i1}^{({\ell })}(x)$ and $M_{i2}^{({\ell }')}(y)$. Alternatively, we can use the formula (10) again for the purpose, so that by the discussion similar to obtaining (15), as $n_{{\ell }},n_{{\ell }'}\rightarrow \infty$, we have

$$\begin{aligned} \int \phi ({\mathrm{d}}x,y)= & {} \int \left\{ E\left( {\mathrm{d}}G_{1}^{({\ell })}(x){\mathrm{d}}G_{2}^{({\ell }')}(y_-)\right) +E\left( {\mathrm{d}}G_{1}^{({\ell })}(x)G_{2}^{({\ell }')}(y_-)\right) \right\} \\\rightarrow & {} \int \phi (x_-,y_-) E\left( \textstyle \sum _{j} \widetilde{H}_{j1}^{\mathrm{a}({\ell })}(x) \mathrm{d}{\overline{M}}_{j1}^{({\ell })}(x) \right) =0. \end{aligned}$$

This yields $\iint b_{1}({\mathrm{d}}x,y)b_{2}(x,{\mathrm{d}}y)=0$. Hence, the solution of (15) is

$$\begin{aligned} \phi (t,s)= \exp \left( -z_{1{\ell }}z_{2{\ell }'} V_{12}(t,s\mid \tau _{{\ell }},\tau _{{\ell }'}) \right) . \end{aligned}$$

Therefore, if $E(\iint {\mathrm{d}}{\overline{M}}_{j1}^{({\ell })}(x){\mathrm{d}}{\overline{M}}_{j2}^{({\ell }')}(y))\ne 0$, the correlation between the two martingales works, which results in $E(G_{1}^{({\ell })}(t)G_{2}^{({\ell }')}(s))\ne 1$ but concludes

$$\begin{aligned} E(G_{1}^{({\ell })}(t)G_{2}^{({\ell }')}(s))\phi (t,s)^{-1}\rightarrow 1 \quad \mathrm{as}\quad n_L\ge \cdots \ge n_1\rightarrow \infty . \end{aligned}$$

In summary, these results provide that the characteristic function of marginal martingale vector $({\mathcal {M}}_{k}^{({\ell })}(t),{\mathcal {M}}_{k'}^{({\ell }')}(s))^\mathrm{T}$ converges to that of bivariate normal distribution as

$$\begin{aligned}&E\left( \mathrm{e}^{ \scriptstyle {\varvec{i}}z_{k{\ell }}{\mathcal {M}}_{k}^{({\ell })}(t) +\scriptstyle {\varvec{i}}z_{k'{\ell }'}{\mathcal {M}}_{k'}^{({\ell }')}(s) }\right) \\&\rightarrow \exp \left( -{\textstyle \frac{1}{2}}z_{k{\ell }}^2 V_{kk}(t,s\mid \tau _{\ell },\tau _{{\ell }})\right. \\&\qquad \left. -z_{k{\ell }}z_{k'{\ell }'} V_{kk'}(t,s\mid \tau _{\ell },\tau _{{\ell }'}) -{\textstyle \frac{1}{2}}z_{k'{\ell }'}^2 V_{k'k'}(t,s\mid \tau _{{\ell }'},\tau _{{\ell }'}) \right) \\&=\left\{ \begin{array}{ll} \exp \left( -2 z_{k{\ell }}^{2} V_{kk}(t,s\mid \tau _{{\ell }},\tau _{{\ell }})\right) &{} \quad \mathrm{if}\quad k=k',{\ell }={\ell }', \\ \exp \left( -\frac{1}{2}\{z_{k{\ell }} V_{kk}(t,s \mid \tau _{{\ell }},\tau _{{\ell }})^{1/2}+z_{k{\ell }'}V_{kk}(t,s\mid \tau _{{\ell }'},\tau _{{\ell }'})^{1/2}\}^{2}\right) &{}\quad \mathrm{if}\quad k=k',{\ell }\ne {\ell }', \\ \hbox {same as the above form} &{}\quad \mathrm{otherwise}.\\ \end{array} \right. \end{aligned}$$

A replication of the similar discussion provides that $({\mathcal {M}}_{1}^{(1)}(t),{\mathcal {M}}_{1}^{(2)}(t),{\mathcal {M}}_{2}^{(1)}(t),$${\mathcal {M}}_{2}^{(2)}(t))$ converges in distribution to a multivariate normal distribution with zero mean vector and covarince matrix

$$\begin{aligned} \left( \begin{array}{cccc} V_{11}(t,s\mid \tau _1,\tau _{1}), &{} &{} &{} \\ V_{11}(t,s\mid \tau _2,\tau _{1}), &{}\quad V_{11}(t,s\mid \tau _2,\tau _{2}), &{} &{} \\ V_{21}(t,s\mid \tau _1,\tau _{1}), &{}\quad V_{21}(t,s\mid \tau _1,\tau _{2}), &{}\quad V_{22}(t,s\mid \tau _1,\tau _{1}), &{}\\ V_{21}(t,s\mid \tau _2,\tau _{1}), &{}\quad V_{12}(t,s\mid \tau _2,\tau _{2}), &{}\quad V_{22}(t,s\mid \tau _2,\tau _{1}), &{}\quad V_{22}(t,s\mid \tau _2,\tau _{2}) \\ \end{array} \right) . \end{aligned}$$

These results lead imidiately to the convergence of $\hat{\varvec{Z}}^* - D_{\scriptstyle {\varvec{n}}}\hat{\varvec{\mu }}$ in distibution to ${\varvec{Z}}^* - D_{\scriptstyle {\varvec{n}}}{\varvec{\mu }}$, as summarized in Theorem 1. $\square$

Some additional results

Table 2 of Sect. 4 displays the results obtained under the assumption of a late time-dependent association (Clayton copula) for the joint survival distribution of the two event-time outcomes. The users may be interested in how the results change if the other types of dependency between two outcomes are assumed. In Table 3, we provide results from the design stage calculated under the same assumptions as Table 2 except that the joint survival distribution is replaced by an early time-dependent association (Gumbel copula). The pattern of the results of MSS, MEN and AEN under Gumbel copula are quite similar to Table 2, but, as the correlation is higher, their reduction rates from the values at zero correlation are slightly larger than those under Clayton copula.

Table 3 Sample sizes, number of events, and empirical powers in a group-sequential trial with two co-primary outcomes under an early time-dependent association (Gumbel copula)

As indicated by one referee, an important matter of concern is how the Type I error rates are controlled or not. In fact, the proposed design method is based on asymptotic results. To answer such a problem, we evaluate the behavior of the actual Type I error rates under sample sizes calculated by the proposed methods. Using ARDENT study, we consider three settings of $(\psi _1,\psi _2)=$ (1.0, 1.0), (0.565, 1.0) and (1.0, 0.721) (both null hypotheses and the two marginals) under the same configurations as Sect. 4, and we confirm the behavior via Monte-Carlo simulation with 1,000,000 runs. For the simulation, a trial ended at the planned follow-up duration. When the observed numbers were larger than the planned ones, the critical value at the final analysis was recalculated based on

$$\begin{aligned} 1-P\left( Z_{k1}<c_{k1},\ldots ,Z_{kL}<\tilde{c}_{kL}\mid \mathrm{H}_{0k}\right) =\alpha _k, \end{aligned}$$

where $\tilde{c}_{kL}$ is the critical value at the final analysis, recalculated such that the above equation is satisfied to control the Type I error adequately if the planned numbers are different from the observed ones.

Table 4 Simulation assessment: probability of rejecting null hypothesis under Clayton copula

Table 5 Simulation assessment: probability of rejecting null hypothesis under Gumbel copula

Table 6 Variance, calendar time and information fraction corresponding to the other endpoint’s information fraction

Tables 4 and 5 show the results of the actual Type I error rates, which are corresponding to the situations under null hypotheses of Tables 2 and 3 under Clayton and Gumbel copulas, respectively. Where the columns “Both” and “ALO” give the probabilities to reject two null hypotheses of OC1 and OC2 jointly (Both) and at least one (ALO), respectively, and “OC1” and “OC2” provide the probabilities to reject two single hypotheses of OC1 and OC2, respectively. We observe that the results of “Joint” are well controlled at the nominal error rate 2.5% in the three cases. Those of “ALO” are less than $2\times 2.5\%$ only at both null hypotheses and reflect the effect of multiplicity using two times testing. Also, the results of “OC1” and “OC2’ are well controlled at the nominal Type I error rate in three cases. Therefore, our method works well in controlling the nominal Type I error rate under the calculated sample size.

Table 1 of Sect. 4 displays the planning information for a group-sequential design at the fixed analysis time points (48 and 96 weeks) considered in ARDENT trial. Other group-sequential designs based on selected information fractions can be constructed. Table 6 displays the planning information for a group-sequential design for information fractions of 0.5 and 1.0.