Tesaro, 1000 Winter St, Waltham, MA 02451
E-mail<EMAIL_ADDRESS>
# Sample size calculation for the Andersen-Gill model comparing rates of
recurrent events

Yongqiang Tang  Ronan Fitzpatrick  Tesaro, 1000 Winter Street, Waltham, MA
02451, USA
Statistical Solutions Ltd., 4500 Avenue 4000, Cork Airport Business Park,
Cork, T12 NX7D, Ireland

(28 November 2018; 11 June 2019; 25 June 2019)

###### Abstract

[Summary] Recurrent events arise frequently in biomedical research, where the
subject may experience the same type of events more than once. The Andersen-
Gill (AG) model has become increasingly popular in the analysis of recurrent
events particularly when the event rate is not constant over time. We propose
a procedure for calculating the power and sample size for the robust Wald test
from the AG model in superiority, noninferiority and equivalence clinical
trials. Its performance is demonstrated by numerical examples. Sample SAS code
is provided in the supplementary material.

###### keywords:

Mixed Poisson process; Noninferiority and equivalence trials; Overdispersion;
Proportional rates/means model; Sandwich variance

††articletype: Research Article

##  1 Introduction

Recurrent events are frequently encountered in biomedical research, where the
subject may experience the same type of events more than once. Examples
include attacks in hereditary angioedema, exacerbations in chronic obstructive
pulmonary disease, bleeds in hemophilia, relapses in multiple sclerosis, and
infections in chronic granulomatous disease (CGD). In clinical trials, the
recurrent events are commonly analyzed by the negative binomial (NB)
regression 1, 2, 3. The NB regression assumes constant event rates over time,
which may fail to hold in some applications 4, 5, 6. The Andersen-Gill (AG)
model 7 provides a popular alternative tool for the analysis of recurrent
events, and it allows arbitrary event rate functions. The AG model often
yields similar treatment effect estimates (i.e. ratio of event rates between
groups) to the NB regression in empirical studies when the event rate is
roughly constant over time 1.

††The paper was published in Statistics in Medicine 2019 (Volume 38, Issue 24,
Pages 4819 - 4827). There was an error in Equation A4 in the appendix. It does
not affect design 1, but appears to slightly overestimate the sample size for
design 2 with staggered entry. The result becomes better after the correction
in the sense that the nominal power generally becomes closer to the simulated
power for design 2. The corrected contents were highlighted in red.

Sample size calculation is critical in designing a clinical trial to ensure
sufficient power to detect an important treatment effect. Sample size
methodology has been well developed for the NB regression; Please see Tang 8
and references therein. Matsui 4 and Song et al 9 derive sample size formulae
for the robust log-rank test 10, which is a nonparametric test suitable only
for superiority trials. In this paper, we propose a power and sample size
calculation procedure for the robust Wald test from the AG model 11. It is
applicable to superiority, noninferiority (NI) and equivalence trials. Two
designs are considered. In one design, the planned treatment duration is the
same for all subjects. In the other design, subjects are enrolled at different
calendar time, but administratively censored at the same calendar time. We
introduce the sample size procedure in Section $2$, and assess its performance
numerically in Section 3.

##  2 Power and sample size formulae

Andersen and Gill 7 provides a simple extension of the Cox proportional
hazards model to the analysis of recurrent events. Suppose $n$ subjects are
randomized to either the active ($x_{i}=1$) or control ($x_{i}=0$) treatment
in a clinical trial. Let $T_{i}$ be the follow-up time for subject $i$,
$Y_{i}(t)=I(T_{i}\geq t)$ the indicator function that subject $i$ is still
under observation at time $t$, and $N_{i}(t)$ the number of events experienced
by subject $i$ by time $t$. Inference for the event rate ratio $\exp(\beta)$
between treatment groups is based on the following partial likelihood

$PL(\beta)=\prod_{i=1}^{n}\prod_{\\{t:Y_{i}(t)=1\\}}\left[\frac{\exp(\beta
x_{i})}{\sum_{j=1}^{n}Y_{j}(t)\exp(\beta x_{j})}\right]^{dN_{i}(t)}.$ (1)

An attractive feature of the AG model is that the baseline event rate function
can be of arbitrary shape. We assume a constant event rate ratio over time,
but the AG model can handle time-varying treatment effects.

To obtain the maximum likelihood estimate (MLE) $\hat{\beta}$, we solve the
score function

$U(\beta)=\frac{\partial\log[PL(\beta)]}{\beta}=\sum_{i=1}^{n}\int_{0}^{\tau}[x_{i}-\bar{x}(\beta,t)]dN_{i}(t)=0,$

where $S^{(k)}(\beta,t)=n^{-1}\sum_{i=1}^{n}Y_{i}(t)x_{i}^{k}\exp(\beta
x_{i})$, $\bar{x}(\beta,t)=\frac{S^{(1)}(\beta,t)}{S^{(0)}(\beta,t)}$, and
$\tau$ is the maximum treatment duration in the trial.

If all covariates are time invariant (the covariates measured after
randomization are rarely used to assess the treatment effect in clinical
trials since the covariates may be affected by the treatment), the AG model
assumes that the time increments between events are independent according to a
Poisson process, but the recurrent events are generally dependent within a
subject 1. The Poisson-type assumption can be relaxed by using the sandwich
variance estimator, and the validity of this robust approach is justified by
Lin et al 11 for arbitrary dependence structures among recurrent events if the
proportional rate or mean assumption is met. For this reason, the robust
approach is also called the proportional rates/means model. The sandwich
variance estimate 11 for $\hat{\beta}$ is
$n^{-1}\hat{V}_{\beta}=n^{-1}\hat{I}_{\beta}^{-1}\hat{\Sigma}_{\beta}\hat{I}_{\beta}^{-1}$,
where
$\hat{\Lambda}_{0}(t)=\sum_{i=1}^{n}\int_{0}^{\tau}[nS^{(0)}(\hat{\beta},t)]^{-1}dN_{i}(t)$,
$d\hat{M}_{i}(t)=dN_{i}(t)-Y_{i}(t)\exp(\hat{\beta}x_{i})d\hat{\Lambda}_{0}(t)$,
$\hat{U}_{i}=\int_{0}^{\tau}[x_{i}-\bar{x}(\hat{\beta},t)]d\hat{M}_{i}(t)$,
$\hat{\Sigma}_{\beta}=n^{-1}\sum_{i=1}^{n}\hat{U}_{i}^{2}$, and
$\hat{I}_{\beta}=n^{-1}\sum_{i=1}^{n}\int_{0}^{\tau}[x_{i}-\bar{x}(\hat{\beta},t)]^{2}Y_{i}(t)\exp(\beta
x_{i})d\hat{\Lambda}_{0}(t)=n^{-1}\sum_{i=1}^{n}\int_{0}^{\tau}[\bar{x}(\hat{\beta},t)-\bar{x}^{2}(\hat{\beta},t)]dN_{i}(t)\,$.
The two-sided $100(1-\alpha)\%$ confidence interval (CI) for $\beta$ is

$[c_{l},c_{u}]=[\hat{\beta}-z_{1-\alpha/2}\sqrt{n^{-1}\hat{V}_{\beta}},\hat{\beta}+z_{1-\alpha/2}\sqrt{n^{-1}\hat{V}_{\beta}}],$

where $z_{p}$ is the $p$-th percentile of the standard normal distribution
$N(0,1)$.

In the sample size determination, we assume a mixed Poisson process (MPP)
model 4, 9, 12 for the event process. Let
$\Lambda_{g}(t)=\text{E}[N_{i}(t)|x_{i}=g]$ be the mean event function for
group $g$. The MPP introduces a random effect $\epsilon_{i}$ with mean $1$ and
variance $\kappa_{g}$ for each subject. Given $\epsilon_{i}$, the subject in
group $g$ follows a Poisson process with mean function
$\epsilon_{i}\Lambda_{g}(t)$. Subjects with $\epsilon_{i}>1$
($\epsilon_{i}<1$) tend to experience more (less) events than the average in
the population. The dispersion parameter $\kappa_{g}$ measures the between-
subject heterogeneity. Inclusion of important risk factors in the model may
reduce heterogeneity 12. The MPP provides a natural way to handle
overdispersion in recurrent events in that the variance of $N_{i}(t)$ is
larger than its mean 12. The mixing distribution for the random effect
$\epsilon_{i}$ is unspecified in the AG model. The NB regression uses a gamma
mixing distribution, and the event count $N_{i}(t)$ follows the NB
distribution 3.

In Appendix A.1, we show that $\hat{V}_{\beta}$ converges in probability to
$V_{\beta}$

$V_{\beta}=\frac{p_{1}[A_{1}+\kappa_{1}B_{1}]+p_{0}[A_{0}+\kappa_{0}B_{0}]}{\left(\int_{0}^{\tau}\frac{[p_{0}\pi_{0}(t)]\,[p_{1}\pi_{1}(t)\exp(\beta)]}{p_{0}\pi_{0}(t)+p_{1}\pi_{1}(t)\exp(\beta)}d\Lambda_{0}\right)^{2}},$
(2)

where $p_{g}$ is the proportion of subjects randomized to treatment group $g$,
$\pi_{g}(t)$ is the probability that a subject in group $g$ remains in the
study at time $t$,
$\omega_{0}(t)=\frac{p_{1}\pi_{1}(t)\exp(\beta)}{p_{1}\pi_{1}(t)\exp(\beta)+p_{0}\pi_{0}(t)}$,
$\omega_{1}(t)=1-\omega_{0}(t)$,
$A_{g}=\int_{t=0}^{\tau}\omega_{g}^{2}(t)\pi_{g}(t)d\Lambda_{g}(t)$, and
$B_{g}=2\int_{t=0}^{\tau}\left[\int_{s=0}^{t}\omega_{g}(s)d\Lambda_{g}(s)\right]\pi_{g}(t)\omega_{g}(t)d\Lambda_{g}(t)$.
We allow the loss to follow-up distribution $G_{g}(t)=1-\pi_{g}(t)$ and the
dispersion parameter $\kappa_{g}$ to differ between the two treatment groups.

At the design stage, it is often reasonable to assume the same dropout
distribution in the two treatment groups (i.e. $\pi_{1}(t)=\pi_{0}(t)$ for all
$t$), and $V_{\beta}$ reduces to

$V_{\beta}=\frac{1}{p_{1}E_{1}}+\frac{1}{p_{0}E_{0}}+2\left(\frac{\kappa_{1}}{p_{1}}\frac{F_{1}}{E_{1}^{2}}+\frac{\kappa_{0}}{p_{0}}\frac{F_{0}}{E_{0}^{2}}\right)=\left[\frac{1}{p_{1}\exp(\beta)}+\frac{1}{p_{0}}\right]\frac{1}{E_{0}}+\left[\frac{\kappa_{1}}{p_{1}}+\frac{\kappa_{0}}{p_{0}}\right]\frac{2F_{0}}{E_{0}^{2}},$
(3)

where $E_{g}=\int\pi_{g}(t)d\Lambda_{g}(t)$ and
$F_{g}=\int_{t=0}^{\tau}\pi_{g}(t)\Lambda_{g}(t)d\Lambda_{g}(t)$. In general,
formula (2) can be well approximated by the term between the two equal signs
in formula (3) if the dropout distribution differs between the two groups.

In Appendix A.2, we provide analytic expressions of $E_{g}$ and $F_{g}$ for
the Weibull and piecewise constant event rate functions when the dropout
pattern is identical in the two groups in two types of clinical trial designs.
In practical applications, almost any event rate function can be approximated
reasonably well by the piecewise constant function.

###  2.1 Superiority and NI trials

Suppose a lower event rate is desirable. In both superiority and NI trials,
the hypothesis can be written as

$H_{0}:\exp(\beta)\geq M_{0}\text{ or }\beta\geq\log(M_{0})\text{ \it versus
}H_{1}:\exp(\beta)<M_{0}\text{ or }\beta<\log(M_{0}).$ (4)

In a superiority trial, the objective is to demonstrate that the experimental
treatment can lower the event rate, and we set $M_{0}=1$. The NI trial aims to
show that the experimental treatment is not worse than the standard control
treatment by $M_{0}$, where $M_{0}>1$ is the prespecifed NI margin on the rate
ratio.

The power for test (4) is given by

$\displaystyle\begin{aligned}
\Pr(c_{u}<&\log(M_{0}))=\Pr\left[Z<\frac{-z_{1-\alpha/2}\sqrt{n^{-1}\hat{V}_{\beta}}-\beta+\log(M_{0})}{\sqrt{n^{-1}V_{\beta}}}\right]\approx\Phi\left[\frac{\sqrt{n}|\log(M_{0})-\beta|}{\sqrt{V_{\beta}}}-z_{1-\alpha/2}\right],\end{aligned}$
(5)

where $Z=(\hat{\beta}-\beta)/\sqrt{n^{-1}V_{\beta}}$ is asymptotically
distributed as $N(0,1)$. The required sample size is

$n=\frac{(z_{1-\alpha/2}+z_{P})^{2}V_{\beta}}{[\log(M_{0})-\beta]^{2}}.$ (6)

As mentioned in Tang 3, Equation (6) is identical to the upper size bound of
Tang 3, 13 for the NB regression (the dispersion parameter may differ between
the two groups in Tang 13) under the assumption of constant event rates if the
dropout pattern is the same in the two groups since
$F_{0}=\lambda_{0}^{2}\text{E}(T_{i}^{2})/2$,
$E_{0}=\lambda_{0}\text{E}(T_{i})$, and

$V_{\beta}=\left[\frac{1}{p_{1}\exp(\beta)}+\frac{1}{p_{0}}\right]\frac{1}{\lambda_{0}\text{E}(T_{i})}+\left[\frac{\kappa_{1}}{p_{1}}+\frac{\kappa_{0}}{p_{0}}\right]\frac{\text{E}(T_{i}^{2})}{\text{E}^{2}(T_{i})}.$

In this special situation, the AG model is almost as powerful as the NB
regression when the variation in the patients’ follow-up time $T_{i}$ is
small, and the two models yield the same power if all subjects have the same
follow-up time $T_{1}=\ldots=T_{n}$. However, the AG model does not require
specifying the mixing distribution.

The NI test is one-sided, and the actual type I error is $\alpha/2$. In
superiority trials, a two-sided test (i.e $H_{0}:\exp(\beta)=1$ vs
$H_{1}:\exp(\beta)\neq 1$) is often used in practice. Formulae (5) and (6) can
be used for the two-sided test since there is little chance that the observed
outcomes will be significantly better in the control group than in the
experimental group if the experimental treatment is truly more effective than
the control treatment 14. The power and sample size formulae (5) and (6)
remain the same if higher event rates indicate better health ($M_{0}\leq 1$)
and the experimental treatment is truly superior or clinically noninferior to
the control treatment in improving the event rate.

###  2.2 Equivalence trials

In an equivalence trial, the objective is to demonstrate that the experimental
treatment is neither superior nor inferior to the standard control treatment.
If the $100(1-\alpha)\%$ CI for $\exp(\beta)$ lies completely within the
interval $[M_{l},M_{u}]$, we can claim clinical equivalence of the two
treatments, where $M_{l}<1$ and $M_{u}>1$ are the prespecified margins. The
hypothesis is

$H_{0}:\exp(\beta)\geq M_{u}\text{ or }\exp(\beta)\leq M_{l}\text{ \it versus
}H_{1}:M_{l}<\exp(\beta)<M_{u}.$

The equivalence test can be viewed as the two one-sided tests and the type I
error is $\alpha/2$. The power is given by

$\displaystyle\begin{aligned}
P&=\Pr(\hat{\beta}+z_{1-\alpha/2}\sqrt{n^{-1}\hat{V}_{\beta}}<\log(M_{u})\text{
and }\hat{\beta}-z_{1-\alpha/2}\sqrt{n^{-1}\hat{V}_{\beta}}>\log(M_{l}))\\\
&\approx\Phi\left(\frac{\sqrt{n}[\log(M_{u})-\beta]}{\sqrt{V_{\beta}}}-z_{1-\alpha/2}\right)-\Phi\left(\frac{\sqrt{n}[\log(M_{l})-\beta]}{\sqrt{V_{\beta}}}+z_{1-\alpha/2}\right).\end{aligned}$
(7)

Formula (7) assumes that
$z_{1-\alpha/2}\sqrt{n^{-1}\hat{V}_{\beta}}+\log(M_{l})<\log(M_{u})-z_{1-\alpha/2}\sqrt{n^{-1}\hat{V}_{\beta}}$
and hence $2z_{1-\alpha/2}\sqrt{n^{-1}\hat{V}_{\beta}}<\log(M_{u}/M_{l})$,
which may not hold with a positive probability. Formula (7) works well in
large samples or when the estimated power is large, but generally
underestimates the power in small samples. The argument is the same as that
for continuous outcomes 15, 14.

The required sample size can be obtained by numerical inversion of the power
formula (7). In the special case when
$\Delta=\log(M_{u})-\beta=\beta-\log(M_{l})=\log(M_{u}/M_{l})/2$, the sample
size is given by

$n=\frac{(z_{1-\alpha/2}+z_{(1+P)/2})^{2}V_{\beta}}{\Delta^{2}}.$ (8)

##  3 Numerical examples

###  3.1 Example 1

We illustrate the sample size calculation for superiority trials by the
analysis of a CGD trial 11, 4. CGD is a rare immune system disorder
characterized by recurrent pyogenic infections. A total of $128$ patients were
randomized to gamma interferon or placebo. The trial was terminated early for
efficacy on basis of an interim analysis of the time to the first infection.
In the trial, $14$ ($22.2\%$) out of $63$ treated patients and $30$ ($46.2\%$)
out of $65$ patients on placebo had at least one infection. Furthermore, $9$
placebo patients and $4$ treated patients experienced at least $2$ infections.

One objective is to estimate the infection rate ratio between the two
treatments. The NB regression gives an estimate of $0.3566$ ($95\%$ CI:
$[0.1934,0.6575]$) while the AG model yields an estimate of $0.3338$ ($95\%$
CI: $[0.1814,0.6143]$). As evidenced by the exploratory analysis of Matsui 4,
the rate of infections may not be constant over time. For this reason, the AG
model is more appropriate for analyzing the CGD trial since it allows
arbitrary event rate function.

Suppose we want to design a new trial to assess the effect of a new
experimental product on the infection rate. We assume the event rate function
is of Weibull form $\lambda_{0}(t)=\psi\nu t^{\nu-1}$ in the placebo arm, and
the event rate ratio between the two treatments is constant
$\lambda_{1}(t)/\lambda_{0}(t)=\exp(\beta)=0.6$ over time. We get the MLE
$(\hat{\psi},\hat{\nu},\hat{\kappa})=(1.097^{1.221},1.221,0.871)$ by fitting a
NB process 12 model to the data using the SAS NLMIXED procedure on basis of
the likelihood function given in Equation (20) of Dean and Balshaw 16. Matsui
4 obtained similar point estimates based on the generalized estimating
equations (GEE) for the MPP 16.

To determine the sample size, we assume a common dispersion parameter and
identical dropout pattern in the two groups. We set $\psi=1.1$, $\nu=1.2$,
$\kappa=0.8$, which are close to the MLE. The treatment allocation ratio is
$p_{1}:p_{0}=1:1$ or $2:1$. We also perform sensitivity analyses to calculate
the sample sizes at alternative parameter values $\kappa=0.4,1.2$, $\psi=1.5$,
$\nu=0.9$. Both design $1$ (planned treatment duration $\tau_{c}=1$ year for
all patients) and design $2$ (accrual period $\tau_{a}=0.5$ year, additional
treatment duration $\tau_{c}=1$ year, constant enrollment rate $\eta=0$) are
considered (please refer to Appendix A.2 for details). In both designs, the
loss to follow-up distribution is exponential with mean $1/\delta=4$ years.

Table 1 reports the sample size and power estimates at the target $90\%$ power
and one-sided type I error $\alpha/2=0.025$. The empirical power is evaluated
based on $40,000$ trials. The data are simulated using Algorithm $2$ of Tang
12 and analyzed using the SAS PHREG procedure. There is more than $95\%$
chance that the simulated power lies within $2\sqrt{0.9*0.1/40000}=0.3\%$ of
the true power.  In both designs, In design $1$, the simulated power is within
$1\%$ of the nominal power in nearly all cases. The performance slightly
deteriorates in design $2$ possibly because of larger variation in the follow-
up time and higher overall dropout rate.

|  |  | Design $1$ | Design $2$
---|---|---|---|---
|  |  | balanced size | unbalanced size | balanced size | unbalanced size
|  |  | total | power ($\%$) | total | power ($\%$) | total | power ($\%$) | total | power ($\%$)
$\kappa$ | $\psi$ | $\nu$ | size | nominal | SIM | size | nominal | SIM | size | nominal | SIM | size | nominal | SIM
$0.4$ | $1.1$ | $0.9$ | $289$ | $90.05$ | $91.03$ | $304$ | $90.00$ | $89.05$ | ${\color[rgb]{1,0,0}256}$ | ${\color[rgb]{1,0,0}90.01}$ | ${\color[rgb]{1,0,0}90.67}$ | ${\color[rgb]{1,0,0}271}$ | ${\color[rgb]{1,0,0}90.04}$ | ${\color[rgb]{1,0,0}89.38}$
|  | $1.2$ | $294$ | $90.01$ | $90.68$ | $310$ | $90.03$ | $89.33$ | ${\color[rgb]{1,0,0}251}$ | ${\color[rgb]{1,0,0}90.10}$ | ${\color[rgb]{1,0,0}90.75}$ | ${\color[rgb]{1,0,0}265}$ | ${\color[rgb]{1,0,0}90.02}$ | ${\color[rgb]{1,0,0}89.34}$
| $1.5$ | $0.9$ | $231$ | $90.12$ | $90.78$ | $244$ | $90.03$ | $89.13$ | ${\color[rgb]{1,0,0}207}$ | ${\color[rgb]{1,0,0}90.05}$ | ${\color[rgb]{1,0,0}90.63}$ | ${\color[rgb]{1,0,0}220}$ | ${\color[rgb]{1,0,0}90.03}$ | ${\color[rgb]{1,0,0}89.11}$
|  | $1.2$ | $235$ | $90.07$ | $90.45$ | $249$ | $90.07$ | $88.90$ | ${\color[rgb]{1,0,0}204}$ | ${\color[rgb]{1,0,0}90.13}$ | ${\color[rgb]{1,0,0}90.54}$ | ${\color[rgb]{1,0,0}217}$ | ${\color[rgb]{1,0,0}90.09}$ | ${\color[rgb]{1,0,0}89.30}$
$0.8$ | $1.1$ | $0.9$ | $358$ | $90.02$ | $90.77$ | $382$ | $90.01$ | $89.15$ | ${\color[rgb]{1,0,0}328}$ | ${\color[rgb]{1,0,0}90.08}$ | ${\color[rgb]{1,0,0}90.67}$ | ${\color[rgb]{1,0,0}351}$ | ${\color[rgb]{1,0,0}90.03}$ | ${\color[rgb]{1,0,0}89.35}$
|  | $1.2$ | $365$ | $90.03$ | $90.77$ | $390$ | $90.05$ | $89.23$ | ${\color[rgb]{1,0,0}324}$ | ${\color[rgb]{1,0,0}90.03}$ | ${\color[rgb]{1,0,0}90.48}$ | ${\color[rgb]{1,0,0}348}$ | ${\color[rgb]{1,0,0}90.03}$ | ${\color[rgb]{1,0,0}89.34}$
| $1.5$ | $0.9$ | $300$ | $90.07$ | $90.39$ | $322$ | $90.03$ | $89.05$ | ${\color[rgb]{1,0,0}278}$ | ${\color[rgb]{1,0,0}90.02}$ | ${\color[rgb]{1,0,0}90.52}$ | ${\color[rgb]{1,0,0}300}$ | ${\color[rgb]{1,0,0}90.02}$ | ${\color[rgb]{1,0,0}89.59}$
|  | $1.2$ | $306$ | $90.08$ | $90.74$ | $328$ | $90.01$ | $89.34$ | ${\color[rgb]{1,0,0}277}$ | ${\color[rgb]{1,0,0}90.04}$ | ${\color[rgb]{1,0,0}90.60}$ | ${\color[rgb]{1,0,0}300}$ | ${\color[rgb]{1,0,0}90.09}$ | ${\color[rgb]{1,0,0}89.42}$
$1.2$ | $1.1$ | $0.9$ | $428$ | $90.06$ | $90.74$ | $460$ | $90.01$ | $89.49$ | ${\color[rgb]{1,0,0}399}$ | ${\color[rgb]{1,0,0}90.05}$ | ${\color[rgb]{1,0,0}90.38}$ | ${\color[rgb]{1,0,0}431}$ | ${\color[rgb]{1,0,0}90.02}$ | ${\color[rgb]{1,0,0}88.97}$
|  | $1.2$ | $436$ | $90.04$ | $90.59$ | $469$ | $90.01$ | $89.55$ | ${\color[rgb]{1,0,0}398}$ | ${\color[rgb]{1,0,0}90.06}$ | ${\color[rgb]{1,0,0}90.24}$ | ${\color[rgb]{1,0,0}431}$ | ${\color[rgb]{1,0,0}90.05}$ | ${\color[rgb]{1,0,0}89.38}$
| $1.5$ | $0.9$ | $369$ | $90.03$ | $90.13$ | $400$ | $90.03$ | $89.17$ | ${\color[rgb]{1,0,0}349}$ | ${\color[rgb]{1,0,0}90.00}$ | ${\color[rgb]{1,0,0}90.56}$ | ${\color[rgb]{1,0,0}380}$ | ${\color[rgb]{1,0,0}90.01}$ | ${\color[rgb]{1,0,0}89.51}$
|  | $1.2$ | $376$ | $90.02$ | $90.24$ | $408$ | $90.04$ | $89.28$ | ${\color[rgb]{1,0,0}351}$ | ${\color[rgb]{1,0,0}90.07}$ | ${\color[rgb]{1,0,0}90.32}$ | ${\color[rgb]{1,0,0}382}$ | ${\color[rgb]{1,0,0}90.02}$ | ${\color[rgb]{1,0,0}89.31}$
Table 1: Estimated sample size at the nominal $90\%$ power and simulated power
(SIM) at the calculated size in designing a new GCD superiority trial
[1] SIM is evaluated using $40,000$ simulated trials.
[2] Losses to follow-up are exponentially distributed with mean $1/\delta=4$ years (annual dropout rate $22.1\%$) in both arms. |  | Unequal dispersion(a) | Unequal dropout(b)
---|---|---|---
|  |  |  | design 1 | design 2 |  |  | design 1 | design 2
|  |  |  | total | power ($\%$) | total | power ($\%$) |  |  | total | power ($\%$) | total | power ($\%$)
$\psi$ | $\nu$ | $\kappa_{0}$ | $\kappa_{1}$ | size | nominal | SIM | size | nominal | SIM | $\kappa_{0}$ | $\kappa_{1}$ | size | nominal | SIM | size | nominal | SIM
1.1 | 0.9 | 0.4 | 0.8 | $324$ | $90.08$ | $91.12$ | ${\color[rgb]{1,0,0}292}$ | ${\color[rgb]{1,0,0}90.05}$ | ${\color[rgb]{1,0,0}90.95}$ | 0.4 | 0.4 | $287$ | $90.06$ | $90.89$ | ${\color[rgb]{1,0,0}254}$ | ${\color[rgb]{1,0,0}90.01}$ | ${\color[rgb]{1,0,0}90.57}$
|  | 0.4 | 1.2 | $358$ | $90.02$ | $91.07$ | ${\color[rgb]{1,0,0}328}$ | ${\color[rgb]{1,0,0}90.08}$ | ${\color[rgb]{1,0,0}91.10}$ | 0.8 | 0.8 | $356$ | $90.02$ | $90.55$ | ${\color[rgb]{1,0,0}326}$ | ${\color[rgb]{1,0,0}90.07}$ | ${\color[rgb]{1,0,0}90.65}$
|  | 0.8 | 1.2 | $393$ | $90.04$ | $90.84$ | ${\color[rgb]{1,0,0}363}$ | ${\color[rgb]{1,0,0}90.03}$ | ${\color[rgb]{1,0,0}90.80}$ | 1.2 | 1.2 | $426$ | $90.06$ | $90.46$ | ${\color[rgb]{1,0,0}397}$ | ${\color[rgb]{1,0,0}90.04}$ | ${\color[rgb]{1,0,0}90.26}$
1.1 | 1.2 | 0.4 | 0.8 | $330$ | $90.06$ | $90.79$ | ${\color[rgb]{1,0,0}287}$ | ${\color[rgb]{1,0,0}90.01}$ | ${\color[rgb]{1,0,0}90.69}$ | 0.4 | 0.4 | $292$ | $90.05$ | $90.37$ | ${\color[rgb]{1,0,0}248}$ | ${\color[rgb]{1,0,0}90.01}$ | ${\color[rgb]{1,0,0}90.36}$
|  | 0.4 | 1.2 | $365$ | $90.03$ | $91.02$ | ${\color[rgb]{1,0,0}324}$ | ${\color[rgb]{1,0,0}90.03}$ | ${\color[rgb]{1,0,0}90.99}$ | 0.8 | 0.8 | $363$ | $90.06$ | $90.58$ | ${\color[rgb]{1,0,0}322}$ | ${\color[rgb]{1,0,0}90.03}$ | ${\color[rgb]{1,0,0}90.66}$
|  | 0.8 | 1.2 | $400$ | $90.00$ | $90.75$ | ${\color[rgb]{1,0,0}361}$ | ${\color[rgb]{1,0,0}90.04}$ | ${\color[rgb]{1,0,0}90.69}$ | 1.2 | 1.2 | $434$ | $90.06$ | $90.29$ | ${\color[rgb]{1,0,0}396}$ | ${\color[rgb]{1,0,0}90.04}$ | ${\color[rgb]{1,0,0}90.33}$
1.5 | 0.9 | 0.4 | 0.8 | $265$ | $90.04$ | $90.82$ | ${\color[rgb]{1,0,0}243}$ | ${\color[rgb]{1,0,0}90.09}$ | ${\color[rgb]{1,0,0}90.94}$ | 0.4 | 0.4 | $229$ | $90.06$ | $90.55$ | ${\color[rgb]{1,0,0}206}$ | ${\color[rgb]{1,0,0}90.11}$ | ${\color[rgb]{1,0,0}90.59}$
|  | 0.4 | 1.2 | $300$ | $90.07$ | $91.26$ | ${\color[rgb]{1,0,0}278}$ | ${\color[rgb]{1,0,0}90.02}$ | ${\color[rgb]{1,0,0}91.12}$ | 0.8 | 0.8 | $298$ | $90.01$ | $90.36$ | ${\color[rgb]{1,0,0}277}$ | ${\color[rgb]{1,0,0}90.05}$ | ${\color[rgb]{1,0,0}90.47}$
|  | 0.8 | 1.2 | $334$ | $90.00$ | $90.70$ | ${\color[rgb]{1,0,0}314}$ | ${\color[rgb]{1,0,0}90.06}$ | ${\color[rgb]{1,0,0}90.89}$ | 1.2 | 1.2 | $368$ | $90.06$ | $90.33$ | ${\color[rgb]{1,0,0}348}$ | ${\color[rgb]{1,0,0}90.02}$ | ${\color[rgb]{1,0,0}90.26}$
1.5 | 1.2 | 0.4 | 0.8 | $270$ | $90.03$ | $90.63$ | ${\color[rgb]{1,0,0}240}$ | ${\color[rgb]{1,0,0}90.02}$ | ${\color[rgb]{1,0,0}90.62}$ | 0.4 | 0.4 | $233$ | $90.05$ | $90.51$ | ${\color[rgb]{1,0,0}202}$ | ${\color[rgb]{1,0,0}90.08}$ | ${\color[rgb]{1,0,0}90.07}$
|  | 0.4 | 1.2 | $306$ | $90.08$ | $91.13$ | ${\color[rgb]{1,0,0}277}$ | ${\color[rgb]{1,0,0}90.04}$ | ${\color[rgb]{1,0,0}90.72}$ | 0.8 | 0.8 | $304$ | $90.05$ | $90.63$ | ${\color[rgb]{1,0,0}276}$ | ${\color[rgb]{1,0,0}90.08}$ | ${\color[rgb]{1,0,0}90.48}$
|  | 0.8 | 1.2 | $341$ | $90.05$ | $90.67$ | ${\color[rgb]{1,0,0}314}$ | ${\color[rgb]{1,0,0}90.06}$ | ${\color[rgb]{1,0,0}90.79}$ | 1.2 | 1.2 | $375$ | $90.06$ | $90.37$ | ${\color[rgb]{1,0,0}349}$ | ${\color[rgb]{1,0,0}90.01}$ | ${\color[rgb]{1,0,0}90.45}$
Table 2: Estimated sample size and simulated power (SIM) at the nominal $90\%$
power in the presence of unequal dropout or dispersion
[1] SIM is evaluated using $10,000$ simulated trials
[2] The treatment allocation ratio is $1:1$
(a) Losses to follow-up are exponentially distributed with mean $1/\delta=4$
years (annual dropout rate $22.1\%$) in both arms.
(b) Losses to follow-up are exponentially distributed with $\delta_{1}=0.15$
and $\delta_{0}=0.35$ (annual dropout rates $13.9\%$ and $29.5\%$) in the two
arms.
[3]  The sample size and nominal power estimates are updated for design 2 with
staggered entry. The simulated power may be different from the previously
reported values after re-runing the simulation for design 1.

###  3.2 Example 2

We conduct simulations to assess the performance of the proposed method in the
presence of unequal dispersion or differential dropout. Two scenarios are
considered. In one scenario, the dispersion parameters in the two groups are
different. In the other scenario, we assume different loss to follow-up
distributions for the two groups. The setup is otherwise similar to that in
the example $1$. The parameter values and simulation results are presented in
Table 2. The performance of the power and sample size method is almost as good
as that in Example $1$.

###  3.3 Example 3

Simulation is conducted to assess the proposed sample size method for NI and
equivalence trials. For illustration purposes, we assume a piecewise constant
event rate function for the control arm $\lambda_{0}(t)=1.0I(0\leq
t<0.4)+1.25I(0.4\leq t<0.8)+1.5I(0.8\leq t\leq 1)$, the event rate ratio
between the active and control arm is
$\exp(\beta)=\lambda_{1}(t)/\lambda_{0}(t)=0.9$ or $1.0$, and the dispersion
parameter is $\kappa=0.8$ or $1.2$. Only design $1$ is considered, and the
planned treatment duration is $\tau_{c}=1$ year for all patients. The
treatment allocation ratio is $1:1$. The loss to follow-up is exponentially
distributed with mean $1/\delta=4$ years. The margin is $M_{0}=1.25$ in the NI
trials, and $(M_{l},M_{u})=(0.75,1.25)$ in the equivalence trials.

Table 3 reports the sample size and power estimates at the target $80\%$ power
and one-sided type I error $\alpha/2=0.025$. The empirical power is evaluated
based on $10,000$ simulated trials. There is more than $95\%$ chance that the
simulated power lies within $2\sqrt{0.8*0.2/10000}=0.8\%$ of the true power.
The simulated power at the calculated sample size is generally close to the
target $80\%$ power, indicating the accuracy of the proposed method.

NI trials (a) | Equivalence trials (b)
---|---
|  | total | power ($\%$) |  |  | total | power ($\%$)
$\kappa$ | $\exp(\beta)$ | size | nominal | SIM | $\kappa$ | $\exp(\beta)$ | size | nominal | SIM
$0.8$ | $0.9$ | $547$ | $80.00$ | $80.26$ | $0.8$ | $0.9$ | $1781$ | $80.02$ | $79.49$
| $1.0$ | $1153$ | $80.03$ | $79.96$ |  | $1.0$ | $1262$ | $80.02$ | $79.58$
$1.2$ | $0.9$ | $675$ | $80.04$ | $80.32$ | $1.2$ | $0.9$ | $2195$ | $80.01$ | $81.20$
| $1.0$ | $1429$ | $80.02$ | $80.31$ |  | $1.0$ | $1564$ | $80.01$ | $80.47$
Table 3:  Estimated sample size at the nominal $80\%$ power and simulated
power (SIM) at the calculated sample size based on $10,000$ NI or equivalence
trials
(a) NI margin is $M_{0}=1.25$
(b) Equivalence margin is $(M_{l},M_{u})=(0.75,1.25)$

##  4 Discussion

We derive the power and sample size formulae for comparing recurrent rates in
superiority, NI and equivalence trials using the robust Wald test from the AG
model. The method allows the dispersion parameter, dropout rate, and/or sample
size to differ between treatment groups. Numerical examples demonstrate the
accuracy of the proposed method in moderate-to-large samples. It is always
recommended to run simulation studies to verify the power calculation
particularly when the sample size is relatively small.

We calculate the variance $V_{\beta}$ and the sample size at given event rate
function, dispersion parameter and dropout rate. These parameters may be
estimated from the historical trials using parametric methods. It is flexible
to adjust the parameter values and conduct sensitivity analyses to examine how
the sample size estimates vary with these parameter values. Please see Example
$1$ for illustration. It is possible to estimate $V_{\beta}$ from the
historical trials by nonparametric methods 9. However, the nonparametric
approach may require that the new trial is sufficiently similar to the
historical trial in terms of the study population, treatment duration, drop
rates, etc.

The robust AG approach has several limitations. First, the AG model uses a
common baseline hazard function for all events, and assumes that the risk of
an event is unaffected by any early events that occur within the same subject.
Therefore, the AG model is not suitable if the occurrence of early events
increases the risk for subsequent ones 1, 17. The AG model provides a
convenient way to estimate an overall treatment effect, but it would be
difficult to estimate the event specific treatment effect 1, 17, which is
useful for studying whether the treatment effect reduces after the patients
experience one or more events. Second, when the sample size is small, the
sandwich variance estimator tends to underestimate the true variance and have
large sampling variability, leading to inflated type I error rate 18, 19. In
the GEE methodology, the bias corrected sandwich variance estimator has been
proposed for small sample inferences 18, 19, 20. It is possible to extend the
bias correction method to the analysis of recurrent events. An alternative
strategy for the analysis of small trials is to use the robust score test
instead of the robust Wald test 20.

##  Appendix A Appendix: Technical details

###  A.1 A brief proof of equations (2) and (3)

By Lin et al 11, $\hat{V}_{\beta}$ is a consistent estimate of $V_{\beta}$

$V_{\beta}=\frac{\text{E}[\text{E}(U_{i}^{2}|x_{i})]}{\text{E}^{2}(I_{\beta})}=\frac{p_{1}\Sigma_{1}+p_{0}\Sigma_{0}}{\left(\int_{0}^{\tau}\omega_{1}(t)\omega_{0}(t)[p_{1}\pi_{1}(t)d\Lambda_{1}+p_{0}\pi_{0}(t)d\Lambda_{0}]\right)^{2}},$
(9)

where $d\,M_{i}(t)=dN_{i}(t)-Y_{i}(t)\exp(\beta x_{i})d\Lambda_{0}(t)$,
$U_{i}=\int_{0}^{\tau}[x_{i}-\bar{x}(\beta,t)]d\,{M}_{i}(t)$,
$\Sigma_{g}=\text{E}(U_{i}^{2}|x_{i}=g)$, and
$I_{\beta}=n^{-1}\sum_{i=1}^{n}\int_{0}^{\tau}[\bar{x}(\beta,t)-\bar{x}^{2}(\beta,t)]dN_{i}(t)$.
By Lemma $1$ in the web-based supplementary material of Song et al 9, we get

$\Sigma_{g}=\text{E}\left[\int_{0}^{\tau}(g-\bar{x}(\beta,t))dM_{i}(t)\int_{0}^{\tau}(g-\bar{x}(\beta,s))dM_{i}(s)\right]=A_{g}+\kappa
B_{g}$ (10)

for subjects in group $g$, where

$\displaystyle\begin{aligned}
A_{g}&=\text{E}\left[\int_{0}^{\tau}Y_{i}(t)\omega_{g}^{2}(t)d\Lambda_{g}(t)\right]=\int\omega_{g}^{2}(t)\pi_{g}(s)d\Lambda_{g}(t),\\\
B_{g}&=\text{E}\left[\int_{0}^{\tau}\int_{0}^{\tau}Y_{i}(t)Y_{i}(s)\omega_{g}(t)\omega_{g}(s)d\Lambda_{g}(t)d\Lambda_{g}(s)\right]=2\int_{s=0}^{\infty}\left[\int_{t=0}^{s}\omega_{g}(t)d\Lambda(t)\right]\pi_{g}(s)\omega_{g}(s)d\Lambda_{g}(s).\end{aligned}$

Inserting Equation (10) into Equation (9) yields Equation (2).

Equation (3) holds under equal dropout since $\omega_{0}(t)\equiv
p_{1}\exp(\beta)/D$, $\omega_{1}(t)\equiv p_{0}/D$,
$I_{\beta}=p_{0}\omega_{0}(t)E_{0}=p_{1}\omega_{1}(t)E_{1}$,
$A_{i}=\omega_{i}^{2}(t)E_{i}$, $B_{i}=2\omega_{i}^{2}(t)F_{i}$,
$E_{1}=E_{0}\exp(\beta)$, $F_{1}=F_{0}\exp(2\beta)$ and
$F_{1}/E_{1}^{2}=F_{0}/E_{0}^{2}$, where $D=p_{0}+p_{1}\exp(\beta)$.

###  A.2 Asymptotic variance expressions in two designs under equal dropout

####  A.2.1 Design 1

The planned treatment duration is $\tau_{c}$ years for each subject (the
accrual period is irrelevant in the sample size calculation). The loss to
follow-up is exponentially distributed with mean $\delta^{-1}$. The
probability that a subject is in the trial at time $t$ after randomization is
$\pi(t)=\exp(-\delta t)$.

Weibull event rate function

Suppose the rate function is $\lambda_{0}(t)=\psi\nu t^{\nu-1}$ and the mean
function is $\Lambda_{0}(t)=\psi t^{\nu}$ for the recurrent event in the
control group, where $\psi$ is a scale parameter and $\nu$ is a shape
parameter. Let $\text{IG}(\nu,a)=\int_{t=0}^{a}t^{\nu-1}\exp(-t)dt$ be the
incomplete gamma function and
$\text{IG}(\nu,a,b)=\int_{t=b}^{a}t^{\nu-1}\exp(-t)dt$. We have

$E_{0}^{\text{I}}=\int_{0}^{\tau_{c}}\pi(t)d\Lambda_{0}(t)=\psi\nu\int_{0}^{\tau_{c}}\exp(-\delta
t)t^{\nu-1}dt=\begin{cases}\frac{\psi\nu}{\delta^{\nu}}\text{IG}(\nu,\delta\tau_{c})&\text{
if }\delta\neq 0\\\ \psi\tau_{c}^{\nu}&\text{ at }\delta=0,\end{cases}$
$F_{0}^{\text{I}}=\int_{0}^{\tau_{c}}\pi(t)\Lambda_{0}(t)d\Lambda_{0}(t)=\psi^{2}\nu\int_{0}^{\tau_{c}}\exp(-\delta
t)t^{2\nu-1}dt=\begin{cases}\frac{\psi^{2}\nu}{\delta^{2\nu}}\text{IG}(2\nu,\delta\tau_{c})&\text{
if }\delta\neq 0\\\ \psi^{2}\tau_{c}^{2\nu}/2&\text{ at }\delta=0.\end{cases}$

Piecewise constant event rate function

Let $\lambda_{0}(t)=\sum_{k=1}^{d}\tilde{\lambda}_{k}I(l_{k-1}\leq t<l_{k})$,
where $l_{0}=0$, $l_{d}=\tau_{c}$. Then
$\Lambda_{0}(t)=\Lambda_{0}(l_{k-1})+\tilde{\lambda}_{k}(t-l_{k-1})$ when
$l_{k-1}\leq t<l_{k}$. Let $\Delta_{k}=l_{k}-l_{k-1}$ and
$G_{km}=\int_{l_{k-1}}^{l_{k}}\exp[-\delta(t-l_{k-1})](t-l_{k-1})^{m}dt=\int_{0}^{\Delta_{k}}\exp[-\delta
t]t^{m}dt$ for $m=0,1,2$. Then

$\begin{cases}G_{k0}=\frac{1-\exp(-\delta\Delta_{k})}{\delta},\,\,G_{k1}=\frac{1-(1+\delta\Delta_{k})\exp(-\delta\Delta_{k})}{\delta^{2}},\,\,G_{k2}=\frac{2-(\delta^{2}\Delta_{k}^{2}+2\delta\Delta_{k}+2)\exp(-\delta\Delta_{k})}{\delta^{3}}&\text{
if }\delta>0\\\
G_{k0}=\Delta_{k},\,\,G_{k1}=\frac{\Delta_{k}^{2}}{2},\,\,G_{k2}=\frac{\Delta_{k}^{3}}{3}&\text{
at }\delta=0\end{cases}$ (11)

We have

$E_{0}^{\text{I}}=\int_{0}^{\tau_{c}}\pi(t)d\Lambda_{0}(t)=\sum_{k=1}^{d}\tilde{\lambda}_{k}\int_{l_{k-1}}^{l_{k}}\pi(t)dt=\begin{cases}\sum_{k=1}^{d}\tilde{\lambda}_{k}\exp[-\delta
l_{k-1}]G_{k0}&\text{ if }\delta\neq 0\\\ \Lambda(\tau_{c})&\text{ at
}\delta=0,\end{cases}$
$F_{0}^{\text{I}}=\int_{0}^{\tau_{c}}\pi(t)\Lambda_{0}(t)d\Lambda_{0}(t)=\sum_{k=1}^{d}\tilde{\lambda}_{k}\int_{l_{k-1}}^{l_{k}}\pi(t)\Lambda_{0}(t)dt=\begin{cases}\sum_{k=1}^{d}\tilde{\lambda}_{k}\exp[-\delta
l_{k-1}][\Lambda(l_{k-1})G_{k0}+\tilde{\lambda}_{k}G_{k1}]&\text{ if
}\delta\neq 0\\\ \Lambda^{2}(\tau_{c})/2&\text{ at }\delta=0.\end{cases}$

####  A.2.2 Design 2

Subjects are enrolled during an accrual period of $\tau_{a}$ years, and
followed for an additional $\tau_{c}$ years after the closure of recruitment.
The total study duration is $\tau=\tau_{a}+\tau_{c}$ years. Suppose the entry
time for a subject is distributed with density function given by

$f(e_{i})=\frac{\eta\exp(-\eta e_{i})}{1-\exp(-\eta\tau_{a})},\text{ where
}0\leq e_{i}\leq\tau_{a}.$

The entry distribution is convex (faster patient entry at the beginning) if
$\eta>0$, and concave (lagging patient entry) if $\eta<0$, and uniform
$f(e_{i})=1/\tau_{a}$ if $\eta\rightarrow 0$. In terms of the sample size
calculation, design $1$ can be viewed as a special case of design $2$ by
setting $\tau_{a}=0$.

Given the entry time $e_{i}$, the maximum follow-up for an individual is
$\tau-e_{i}$. We assume the loss to follow-up is exponentially distributed
with mean $\delta^{-1}$. The probability that a subject is still in the trial
at time $t$ after randomization is

$\displaystyle\begin{aligned}
\pi(t)&=\Pr(T_{i}>t)=\Pr(T_{i}>t|e_{i}+t\leq\tau)\Pr(e_{i}+t\leq\tau)+\Pr(T_{i}>t|e_{i}+t>\tau)\Pr(e_{i}+t>\tau)\\\
&=\begin{cases}\exp(-\delta t)\text{ if }t\leq\tau_{c}\\\ \exp(-\delta
t){\color[rgb]{1,0,0}\frac{1-\exp[-\eta(\tau-t)]}{1-\exp(-\eta\tau_{a})}}\text{
if }\tau_{c}<t\leq\tau.\end{cases}\end{aligned}$ (12)

When $\eta\rightarrow 0$, we shall replace
${\color[rgb]{1,0,0}\frac{1-\exp[-\eta(\tau-t)]}{1-\exp(-\eta\tau_{a})}}$ by
its limiting value ${\color[rgb]{1,0,0}(\tau-t)/\tau_{a}}$ in Equation (12).
In design 2, it is easy to see that

$\int_{0}^{\tau}\pi(t)d\Lambda_{0}(t)=E^{\text{I}}+\int_{\tau_{c}}^{\tau}\pi(t)d\Lambda_{0}(t)\text{
and
}\int_{0}^{\tau}\pi(t)\Lambda_{0}(t)d\Lambda_{0}(t)=F_{0}^{\text{I}}+\int_{\tau_{c}}^{\tau}\pi(t)\Lambda_{0}(t)d\Lambda_{0}(t),$

where $E_{0}^{\text{I}}$ and $F_{0}^{\text{I}}$ are defined in Appendix A.2.1.
Below we give analytic expression for
$\int_{\tau_{c}}^{\tau}\pi(t)d\Lambda_{0}(t)$ and
$\int_{\tau_{c}}^{\tau}\pi(t)\Lambda_{0}(t)d\Lambda_{0}(t)$ at $\eta=0$. The
expressions are omitted when $\eta\neq 0$ due to limited space.

Weibull event rate function

Suppose $\lambda_{0}(t)=\psi\nu t^{\nu-1}$. When $\eta=0$, we get

$\displaystyle\begin{aligned}
\int_{\tau_{c}}^{\tau}\pi(t)d\Lambda_{0}(t)&=\begin{cases}{\color[rgb]{1,0,0}\frac{\tau\psi\nu}{\tau_{a}\delta^{\nu}}\text{IG}(\nu,\delta\tau,\delta\tau_{c})-\frac{\psi\nu}{\tau_{a}\delta^{\nu+1}}\text{IG}(\nu+1,\delta\tau,\delta\tau_{c})}&\text{
if }\delta>0\\\
{\color[rgb]{1,0,0}\frac{\tau\psi}{\tau_{a}}[\tau^{\nu}-\tau_{c}^{\nu}]-\frac{\psi\nu}{\tau_{a}(\nu+1)}[\tau^{\nu+1}-\tau_{c}^{\nu+1}]}&\text{
at }\delta=0\end{cases}\\\
\int_{\tau_{c}}^{\tau}\pi(t)\Lambda_{0}(t)d\Lambda_{0}(t)&=\begin{cases}{\color[rgb]{1,0,0}\frac{\tau\psi^{2}\nu}{\tau_{a}\delta^{2\nu}}\text{IG}(2\nu,\delta\tau,\delta\tau_{c})-\frac{\psi^{2}\nu}{\tau_{a}\delta^{2\nu+1}}\text{IG}(2\nu+1,\delta\tau,\delta\tau_{c})}&\text{
if }\delta>0\\\
{\color[rgb]{1,0,0}\frac{\tau\psi^{2}}{2\tau_{a}}[\tau^{2\nu}-\tau_{c}^{2\nu}]-\frac{\psi^{2}\nu}{\tau_{a}(2\nu+1)}[\tau^{2\nu+1}-\tau_{c}^{2\nu+1}]}&\text{
at }\delta=0\end{cases}\\\ \end{aligned}$ (13)

Piecewise constant event rate function

Suppose $\lambda_{0}(t)=\sum_{k=1}^{d^{*}}\tilde{\lambda}_{k}I(l_{k-1}\leq
t<l_{k})$, where $l_{d^{*}}=\tau=\tau_{a}+\tau_{c}$ and $l_{d}=\tau_{c}$. For
notational convenience, if $\tau_{c}$ is not a knot, it can be added as a
knot. When $\eta=0$,

$\displaystyle\begin{aligned}
&\int_{\tau_{c}}^{\tau}\pi(t)d\Lambda_{0}(t)={\color[rgb]{1,0,0}\sum_{k=d+1}^{d^{*}}\int_{l_{k-1}}^{l_{k}}\lambda_{k}\exp(-\delta
t)\frac{\tau-t}{\tau_{a}}dt=\sum_{k=d+1}^{d^{*}}\frac{\lambda_{k}}{\tau_{a}}\exp(-\delta
l_{k-1})\left[(\tau-l_{k-1})G_{k0}-G_{k1}\right]},\\\
&\int_{\tau_{c}}^{\tau}\pi(t)\Lambda_{0}(t)d\Lambda_{0}(t)={\color[rgb]{1,0,0}\sum_{k=d+1}^{d^{*}}\frac{\lambda_{k}}{\tau_{a}}\exp(-\delta
l_{k-1})\left\\{\Lambda_{k-1}[(\tau-l_{k-1})G_{k0}-G_{k1}]+\lambda_{k}[(\tau-
l_{k-1})G_{k1}-G_{k2}]\right\\}},\end{aligned}$

where $G_{k0}$, $G_{k1}$ and $G_{k2}$ are defined in Equation (11).

## References

  * 1 Wang Y, Meyerson L, Tang Y, Qian N. Statistical methods for the analysis of relapse data in MS clinical trials. Journal of the Neurological Sciences 2009; 285: 206 - 11.
  * 2 Aban IB, Cutter GR, Mavinga N. Inferences and power analysis concerning two negative binomial distributions with an application to MRI lesion counts data. Computational Statistics & Data Analysis 2009; 53: 820 -33.
  * 3 Tang Y. Sample Size Estimation for Negative Binomial Regression Comparing Rates of Recurrent Events with Unequal Follow-Up Time. Journal of Biopharmaceutical Statistics 2015; 25: 1100 - 13.
  * 4 Matsui S. Sample size calculations for comparative clinical trials with over-dispersed Poisson process data. Statistics in Medicine 2005; 24: 1339 - 56.
  * 5 Nicholas R, Straube S, Schmidli H, Schneider S, Friede T. Trends in annualized relapse rates in relapsing-remitting multiple sclerosis and consequences for clinical trial design. Multiple sclerosis 2011; 17: 1211 \- 17.
  * 6 Inusah S, Sormani MP, Cofield SS, et al. Assessing changes in relapse rates in multiple sclerosis. Multiple sclerosis 2010; 16: 1414 - 21\.
  * 7 Andersen PK, Gill RD. Cox’s regression model counting process: a large sample study. Annual of Statistics 1982; 10: 1100 - 20.
  * 8 Tang Y. Negative Binomial Regression: Sample Size with Unequal Follow-Up Times. In: CRC Press; 2017.
  * 9 Song R, Kosorok MR, Cai J. Robust covariate-adjusted log-rank statistics and corresponding sample size formula for recurrent events data. Biometrics 2008; 92: 741 -50.
  * 10 Lawless JF, Nadeau C. Some Simple Robust Methods for the Analysis of Recurrent Events. Technometrics 1995; 37: 158 - 68.
  * 11 Lin DY, Wei LJ, Yang I, Ying Z. Semiparametric Regression for the Mean and Rate Functions of Recurrent Events. Journal of the Royal Statistical Society B 2000; 62: 711 - 30.
  * 12 Tang Y. Algorithms for imputing partially observed recurrent events with applications to multiple imputation in pattern mixture models. Journal of Biopharmaceutical Statistics 2018; 28: 518-533.
  * 13 Tang Y. Sample size for comparing negative binomial rates in noninferiority and equivalence trials with unequal follow-up time. Journal of Biopharmaceutical Statistics 2018; 28: 475-491.
  * 14 Tang Y. Exact and approximate power and sample size calculations for analysis of covariance in randomized clinical trials with or without stratification. Journal of Biopharmaceutical research 2018; 10: 274-286.
  * 15 Tang Y. A Noniterative Sample Size Procedure for Tests Based on t Distributions. Statistics in Medicine 2018; 37: 3197 - 3213\.
  * 16 Dean CB, Balshaw R. Efficiency Lost By Analyzing Counts Rather Than Event Times in Poisson and Overdispersed Poisson Regression Models. Journal of the American Statistical Association 1997; 92: 1387 - 98.
  * 17 Amorim LD, Cai J. Modelling recurrent events: a tutorial for analysis in epidemiology. International journal of epidemiology 2015; 44: 1 - 10.
  * 18 Kauermann G, Carroll RJ. A note on the efficiency of sandwich covariance matrix estimation. Journal of the American Statistical Association 2001; 96: 1387 - 96.
  * 19 Lu B, Preisser JS, Qaqish BF, Suchindran C, Bangdiwala S, Wolfson M. A comparison of two bias-corrected covariance estimators for generalized estimating equations. Biometrics 2007; 63: 935-41.
  * 20 Guo X, Pan W, Connett JE, Hannan PJ, French SA. Small-sample performance of the robust score test and its modifications in generalized estimating equations. Statistics in Medicine 2005; 24: 3479 - 95.