1  Computing expected number of events over time

1.1 General formulation

We will consider a single treatment group and stratum and assume subjects enroll according to a Poisson process with entry rate \(g(t)\geq 0\) for \(t>0\). Denote the number of subjects enrolled at or before time \(t\) by \(N(t)\) for \(t\geq 0\). The expected number of subjects enrolled by time \(t\) is simply \[\begin{equation} E\{N(t)\}=G(t)=\int_0^t g(s)ds. \end{equation}\] We assume the study runs from time 0 to some time \(t_m>0\) where \(m\) is a positive integer we will further define shortly. We assume the time to event is independent and identically distributed for all subjects enrolled. We also assume the censoring time is independent and identically distributed for all subjects enrolled. For an individual, let \(X>0\) denote the duration from entry into the trial until an event and \(Y>0\) denote the duration from entry into the trial until loss-to-follow-up. The duration of time a subject is followed is \(T=\min(X,Y)\). We assume the pair \(X\), \(Y\) is independent of entry time.

Denote the number of events observed by time \(t\) by \(D(t)\) for \(t\geq 0\). The expected number of events observed by time \(t_m\) is \[\begin{equation} E\{D(t_m)\}=\int_0^{t_m} P\{X\leq t_m-u, X\leq Y\} g(u)du. \end{equation}\] With a change of variables from \(u\) to \(t=t_m-u\) we have \[\begin{equation} E\{D(t_m)\} =\int_0^{t_m} P\{X\leq t, X\leq Y\} g(U\leq t_m-t)dt.\label{eq:Edelta} \end{equation}\]

1.2 Piecewise failure and enrollment distibutions

Here we define a piecewise random variable and key related probabilities and follow with a corresponding general formulation for calculating \(P\{\delta=1\}\).

1.2.1 General definitions

Recall that \(m\) is a positive integer. For \(i=1,2,\ldots,m\) we assume We assume \(0=t_0<t_1<...<t_m\) and for \(i=1,2,\ldots,m\) let \(s_i=t_i-t_{i-1}\). We define random variables that will be used to define \(X\) and \(Y\) on a piecewise basis over the intervals \((t_{i-1}, t_{i}\)], \(i=1,2,\ldots,m\) as follows.

  • For \(i=1,2,\ldots,m\), we assume \(X_i>0\), \(Y_i>0\) are random variables that are independent of the study entry time \(U\).
  • We let \(X_i\) and \(Y_i\) define \(X\) and \(Y\), respectively, on the interval \((t_{i-1},t_i]\), \(i=1,2,\ldots,m\), as follows.
    • If \(X_1\leq t_1=s_1\), then \(X=X_1\).
    • For \(i=2,3,\ldots,m-1\) recursively define \(X=t_{i-1}+X_i\) if \(X\) has not been yet defined and \(X_i\leq s_i\).
    • If \(X\) has not been otherwise defined above, let \(X=t_{m-1}+X_m\).

\(Y\) is defined analogously by \(Y_1,Y_2,\ldots,Y_m\).

For \(i=1,2,\ldots,m\) we assume \[\begin{equation} q_i=P\{X_i>s_i, Y_i>s_i\}\label{eq:qi} \end{equation}\] \[\begin{equation} Q_i\equiv P\{X>t_i,Y>t_i\}=\prod_{j=1}^iP\{X_i>s_i, Y_i>s_i\}=\prod_{j=1}^i q_j. \end{equation}\] \[\begin{equation} p_i=P\{t_{i-1}<X\leq t_i,X\leq Y\}=Q_{i-1}P\{X_i\leq s_i,X_i\leq Y_i\} \end{equation}\] and let \[\begin{equation} P_i\equiv P\{T\leq t_i, X\leq Y\} =\sum_{j=1}^i p_j. \end{equation}\] As an example, if \(X_i\), \(Y_i\), \(i=1,2,\ldots,m\) are all independent these probability assumptions are automatically satisified.

We will break the calculation of interest into intervals. Denote \[\begin{equation} E\{D(t_m)\}= \sum_{i=1}^m A_i \end{equation}\] where for \(i=1,2,\ldots,m\) \[\begin{eqnarray} A_i&=&\int_{t_{i-1}}^{t_i}P\{X\leq t_i, X\leq Y\}g(t_m-t)dt \nonumber\\ &=&P_{i-1}(G(t_m-t_{i-1})-G(t_m-t_i))+Q_{i-1}B_i\label{eq:genEdelta} \end{eqnarray}\] where \[\begin{equation} B_i=\int_0^{s_i}P\{X_i\leq s,X_i\leq Y_i\}g(t_m-(t_{i-1}+s))ds \label{eq:Bi}\end{equation}\] We are now left to calculate \(p_i\), \(q_i\), and \(B_i\), \(i=1,2,\ldots,m\) given specific distributional assumptions. Applying equations (\(\ref{eq:qi}\))-(\(\ref{eq:Bi}\)) we can then compute \(E\{D(t_m)\}.\)

1.2.2 Definition of a piecewise exponential failure time

For \(i=1,2,\ldots,m\) we assume \(X_i\) and \(Y_i\) are independent and exponentially distributed with failure rates \(\lambda_i\) and \(\eta_i\), respectively, where \(X_i\) will be referred to as an event time and \(Y_i\) a censoring or lost-to-follow-up time. The random variables \(X_i\) and \(Y_i\) will be used to define \(T\) and \(\delta\) over the interval \((t_{i-1},t_i]\). We define \(\delta_i=1\) if \(\min(X_i, Y_i)\leq s_i\) and \(\delta_i=0\) otherwise, \(i=1,2,\ldots,m\). We will let \(J\) be the first interval \(i\) in which \(\min(X_i,Y_i)\leq s_i\) and \(J=m+1\) if \(\min(X_i,Y_i)>s_i\), \(i=1,2,\ldots,m\). We define the time to event \(T\equiv t_m\) if \(J=m+1\) and \(T\equiv t_{J-1}+\min(X_J,Y_J)\), otherwise. We define the indicator that \(T\) is an event time rather than a lost-to-follow-up time by the random variable \(\delta^\prime=1\) if \(J<m+1\) and \(T=t_{J-1}+X_J\). The random variable \(T\) will be referred to as a piecewise exponential failure time with event indicator \(\delta^\prime\).

1.2.3 Exponential failure and censoring

We assume \(m=1\), \(X\sim\text{Exponential}(\lambda)\) and \(Y\sim\text{Exponential}(\eta)\). We assume uniform entry over a time \((0,r]\) where \(0 < r \leq t_1\) which implies \(g(u) = 1/r\) on \((0, r]\) and is 0 elsewhere. It is straightforward to calculate that \[ P\{X\leq t,X\leq Y\}=\frac{\lambda}{\lambda+\eta}(1-e^{-(\lambda+\eta)t}) \]

From (\(\ref{eq:genEdelta}\)) we only need calculate \[\begin{eqnarray} E\{D(t_m)\}&=&B_1\nonumber\\ &=&\int_0^{t_1}\frac{\lambda}{\lambda+\eta}\left(1-e^{-(\lambda+\eta)t}\right) g(t_1-t)dt \nonumber\\ &=&\frac{\lambda}{(\lambda+\eta)r} \int_{t_1-r}^{t_1}\left(1-e^{-(\lambda+\eta)t}\right)dt \nonumber\\ &=&\frac{\lambda}{\lambda+\eta} \left( 1-\frac{e^{-(\lambda+\eta)(t_1-r)}-e^{-(\lambda+\eta)t_1}}{(\lambda+\eta)r} \right).\label{eq:exppdel} \end{eqnarray}\] which corresponds with Lachin and Foulkes (1986).

With entry following a truncated exponential distribution for \(0<u\leq r<t_1\) and \(\gamma\neq 0\) and \(\gamma\neq \lambda+\eta\) \[ g(u)=\frac{e^{-\gamma u}}{1-e^{-\gamma r}} \] it is straightforward to show \[\begin{eqnarray} E\{D(t_1)\}&=&\int_{t_1-r}^{t_1}\frac{\lambda}{\lambda+\eta} (1-e^{-(\lambda+\eta)t}) \gamma e^{-\gamma(t_1-t)}dt\label{eq:gamEdelta}\\ &=&\frac{\lambda}{\lambda+\eta}\left(1+ \frac{\gamma e^{-(\lambda+\eta)t_1}}{\lambda+\eta-\gamma}\left( \frac{1-e^{(\lambda+\eta-\gamma)r}}{1-e^{-\gamma r}}\right)\right) \label{eq:gamEdelta2} \end{eqnarray}\] which again is consistent with Lachin and Foulkes (1986). If \(\gamma=\lambda+\eta\), it is straightforward to show either directly by integrating (\(\ref{eq:gamEdelta}\)) or by applying l’Hôpital’s rule to (\(\ref{eq:gamEdelta2}\)) that \[\begin{equation} E\{D(t_1)\}=\lambda\left(\frac{1}{\gamma}- \frac{re^{-\gamma t_1}}{1-e^{-\gamma r}}\right). \end{equation}\]

1.2.4 Piecewise constant enrollment

We now assume subjects enroll at a constant rate in each of \(m\) intervals. Without loss of generality, we assume that for \(i=0,1,2,\ldots,m\), that the enrollment rate is [g(t)=_{m+1-i}] for \(t\) in the interval \((t_m-t_i,t_m-t_{i-1})\). If subject recruiting does not occur throughout the course of the trial, some of the \(\gamma_i\) values may be 0.

1.2.5 Expected number of events

We now apply equations (\(\ref{eq:qi}\))-(\(\ref{eq:Bi}\)) to compute \(P\{\delta=1\}\) using the following. \[\begin{equation} q_i=e^{-(\lambda_i+\eta_i)s_i}\label{eq:pwqi} \end{equation}\] \[\begin{equation} p_i=Q_{i-1}\frac{\lambda_i}{\lambda_i+\eta_i}(1-e^{-(\lambda_i+\eta_i)s_i}) \end{equation}\] For \(0\leq s<s_i\) we have \(t_m-t_i<t_m-(t_{i-1}-s)\leq t_m-t_i\) which implies \(g(t_m-(t_{i-1}-s))=\gamma_{m+1-i}\) which implies \[\begin{eqnarray} B_i&=&\int_0^{s_i}\frac{\gamma_{m+1-i}\lambda_i}{\lambda_i+\eta_i} (1-e^{-(\lambda_i+\eta_i)s})ds\nonumber\\ &=&\frac{\gamma_{m+1-i}\lambda_i}{\lambda_i+\eta_i} \left(s_i- \frac{1-e^{-(\lambda_i+\eta_i)s_i}}{\lambda_i+\eta_i}\right).\label{eq:pwpdel} \end{eqnarray}\] We note that when \(\lambda_i+\eta_i=0\) that \(p_i=0\) and \(B_i=0\) by definition; this is accounted for in eEvents1(). As an exercise, show that (\(\ref{eq:pwqi}\))-(\(\ref{eq:pwpdel}\)) can be used along with (\(\ref{eq:qi}\))-(\(\ref{eq:Bi}\)) to produce (\(\ref{eq:exppdel}\)).

1.2.6 Testing event count calculations

We assume a study duration of \(T=20\) months. Enrollment is assumed to occur over the course of 5 months at a rate of 5 per month for the first 2 months, 10 per month for 1 month, and 20 per month for months 4 and 5. This results in an expected enrollment of 60. The uniform rates to reflect enrollment probability over time for a given study subject are 1/12 from 0 to 2, 2/12=1/6 from 2 to 3, and 4/12=1/3 from 3 to 5; enrollment rates are 0 thereafter. The exponential event rate is assumed to be high initially with a rate of .05 in the first month, .02 for the second month, and .01 per month thereafter. The rate of censoring is assumed to be exponential with a rate of .03 per month from time of enrollment through 20 months. Based on this, the \(t_i\) values are 0, 1, 2, 3, 5 and corresponding timepoints at the end of the study 20, 19, 18, 17, 15. Summing \(A_i\) from 1 to 6 yields \(E\{D(20)\}=11.023\). The calculations in the following table were computed in an Excel spreadsheet and can be used for checking the following calculations in R.

\(i\) \(t_i\) \(\gamma_{m+1-i}\) \(\lambda_i\) \(\eta_i\) \(q_i\) \(Q_i\) \(p_i\) \(P_i\) \(B_i\) \(A_i\)
1 1 0 .05 .01 0.9418 0.9418 0.0485 0.0485 0.0000 0.0000
2 2 0 .02 .01 0.9704 0.9139 0.0186 0.0671 0.0000 0.0000
3 15 0 .01 .01 0.7711 0.7047 0.1046 0.1717 0.0000 0.0000
4 17 20 .01 .01 0.9608 0.6771 0.0138 0.1855 0.3947 7.1464
5 18 10 .01 .01 0.9802 0.6637 0.0067 0.1922 0.0497 1.8889
6 20 5 .01 .01 0.9608 0.6376 0.0130 0.2052 0.0987 1.9877

1.2.7 R routine for computing expected number of events for a single stratum

Calculations for the expected number of events as just outlined are performed using the function eEvents1(). This is not accessible when the gsDesign package is loaded since input variable checks are not performed for this lower-level function. The more general function eEvents() will be presented in the next section. It performs the same calculations for a stratified population.

Following are the arguments for eEvents1():

  • lambda: a vector of event rates
  • eta: a vector of dropout rates; if length is 1, this will be transformed to the length of lambda
  • gamma: a vector of enrollment rates
  • R: a vector of durations for enrollment rates; should have the same length as gamma; note that if sum(R) > Tfinal - minfup then the actual accrual is limited to Tfinal - minfup but that the output R is still as input
  • S: a vector of durations of failure rate periods. If lambda has length 1, this should be NULL; otherwise, this should have length 1 less than lambda. The final failure rate in lambda is assumed to continue indefinitely, and thus is not specified in S
  • T: a scalar indicating the time at which the expected number of events and enrollment is to be computed
  • Tfinal: a scalar indicating the time at which the trial is to be completed. If NULL, this is assumed the same as T
  • minfup: a scalar indicating the minimum follow-up for each observation at the end of the trial. Note that min(Tfinal - minfup, sum(R)) is when enrollment ends and thus we must have Tfinal > minfup.

The reader may skip the next section as examples later demonstrate application of the routing.

1.2.8 The eEvents1 source code

# Do not edit source code in gsSurv.R
# This should be modified in https://github.com/keaven/gsSurv
# eEvents1 function [sinew] ----
eEvents1 <- function(lambda = 1, eta = 0, gamma = 1, R = 1, S = NULL,
                     T = 2, Tfinal = NULL, minfup = 0, simple = TRUE) {
  if (is.null(Tfinal)) {
    Tfinal <- T
    minfupia <- minfup
  } else {
    minfupia <- max(0, minfup - (Tfinal - T))
  }

  nlambda <- length(lambda)
  if (length(eta) == 1 & nlambda > 1) {
    eta <- rep(eta, nlambda)
  }
  T1 <- cumsum(S)
  T1 <- c(T1[T1 < T], T)
  T2 <- T - cumsum(R)
  T2[T2 < minfupia] <- minfupia
  i <- 1:length(gamma)
  gamma[i > length(unique(T2))] <- 0
  T2 <- unique(c(T, T2[T2 > 0]))
  T3 <- sort(unique(c(T1, T2)))
  if (sum(R) >= T) T2 <- c(T2, 0)
  nperiod <- length(T3)
  s <- T3 - c(0, T3[1:(nperiod - 1)])

  lam <- rep(lambda[nlambda], nperiod)
  et <- rep(eta[nlambda], nperiod)
  gam <- rep(0, nperiod)

  for (i in length(T1):1)
  {
    indx <- T3 <= T1[i]
    lam[indx] <- lambda[i]
    et[indx] <- eta[i]
  }
  for (i in min(length(gamma) + 1, length(T2)):2) {
    gam[T3 > T2[i]] <- gamma[i - 1]
  }
  q <- exp(-(c(lam) + c(et)) * s)
  Q <- cumprod(q)
  indx <- 1:(nperiod - 1)
  Qm1 <- c(1, Q[indx])
  p <- c(lam) / (c(lam) + c(et)) * c(Qm1) * (1 - c(q))
  p[is.nan(p)] <- 0
  P <- cumsum(p)
  B <- c(gam) / (c(lam) + c(et)) * c(lam) * (c(s) - (1 - c(q)) / (c(lam) + c(et)))
  B[is.nan(B)] <- 0
  A <- c(0, P[indx]) * c(gam) * c(s) + c(Qm1) * c(B)
  if (!simple) {
    return(list(
      lambda = lambda, eta = eta, gamma = gamma, R = R, S = S,
      T = T, Tfinal = Tfinal, minfup = minfup, d = sum(A),
      n = sum(gam * s), q = q, Q = Q, p = p, P = P, B = B, A = A, T1 = T1,
      T2 = T2, T3 = T3, lam = lam, et = et, gam = gam
    ))
  } else {
    return(list(
      lambda = lambda, eta = eta, gamma = gamma, R = R, S = S,
      T = T, Tfinal = Tfinal, minfup = minfup, d = sum(A),
      n = sum(c(gam) * c(s))
    ))
  }
}
eEvents1(
  lambda = c(1, 2, 3, 4), eta = 0, gamma = c(1, 2), R = c(1, 5), S = c(1.5, 2.5, 3.5),
  T = 10, Tfinal = NULL, minfup = 0, simple = FALSE
)
$lambda
[1] 1 2 3 4

$eta
[1] 0 0 0 0

$gamma
[1] 1 2

$R
[1] 1 5

$S
[1] 1.5 2.5 3.5

$T
[1] 10

$Tfinal
[1] 10

$minfup
[1] 0

$d
[1] 10.999

$n
[1] 11

$q
[1] 2.231302e-01 6.737947e-03 2.753645e-05 2.478752e-03 1.831564e-02

$Q
[1] 2.231302e-01 1.503439e-03 4.139938e-08 1.026188e-10 1.879529e-12

$p
[1] 7.768698e-01 2.216267e-01 1.503398e-03 4.129676e-08 1.007393e-10

$P
[1] 0.7768698 0.9984966 1.0000000 1.0000000 1.0000000

$B
[1] 0.0000000 0.0000000 6.3333517 2.5012394 0.7545789

$A
[1] 0.000000 0.000000 6.998998 3.000000 1.000000

$T1
[1]  1.5  4.0  7.5 10.0

$T2
[1] 10  9  4

$T3
[1]  1.5  4.0  7.5  9.0 10.0

$lam
[1] 1 2 3 4 4

$et
[1] 0 0 0 0 0

$gam
[1] 0 0 2 2 1

We provide an example where both \(\lambda_i=0\) and \(\eta_i=0\) for some \(i\). This option will be used later when follow-up is cut off at some point in time for a study.

eEvents1(
  lambda = c(1, 2, 3, 0), eta = 0, gamma = c(1, 2), R = c(1, 5), S = c(1.5, 2.5, 3.5),
  T = 10, Tfinal = NULL, minfup = 0, simple = FALSE
)
$lambda
[1] 1 2 3 0

$eta
[1] 0 0 0 0

$gamma
[1] 1 2

$R
[1] 1 5

$S
[1] 1.5 2.5 3.5

$T
[1] 10

$Tfinal
[1] 10

$minfup
[1] 0

$d
[1] 10.999

$n
[1] 11

$q
[1] 2.231302e-01 6.737947e-03 2.753645e-05 1.000000e+00 1.000000e+00

$Q
[1] 2.231302e-01 1.503439e-03 4.139938e-08 4.139938e-08 4.139938e-08

$p
[1] 0.776869840 0.221626721 0.001503398 0.000000000 0.000000000

$P
[1] 0.7768698 0.9984966 1.0000000 1.0000000 1.0000000

$B
[1] 0.000000 0.000000 6.333352 0.000000 0.000000

$A
[1] 0.000000 0.000000 6.998998 3.000000 1.000000

$T1
[1]  1.5  4.0  7.5 10.0

$T2
[1] 10  9  4

$T3
[1]  1.5  4.0  7.5  9.0 10.0

$lam
[1] 1 2 3 0 0

$et
[1] 0 0 0 0 0

$gam
[1] 0 0 2 2 1

1.2.9 Extending expected number of events calculations to a stratified population

We extend our notation to describe a stratified population here, to show the calculations we wish to perform.

The stratified case is made reasonably general, allowing enrollment, event rates and dropout rates to change over time in whatever piecewise fashion is required. Without loss of generality, we assume the same time intervals for constant enrollment and failure rates as before. If these differ by strata, then any time where a change is made for any stratum is included for all strata in the notation. Whereas we had vectors representing enrollment rates, failure rates and censoring rates previously we now have matrices with the column of each matrix representing a vector for each stratum. In an effort to make this explicit, we provide notation. We assume there are \(s\geq 1\) strata in the population. For \(1\leq i\leq m\) we previously let \(\lambda_i\) represent a failure rate in the \(i\)-th time period. We now let \(\boldsymbol{\lambda}\) represent a \(m\times s\) matrix with failure rates for a stratum in each column and for \(1\leq i\leq m\), \(1\leq j \leq s\) we let \(\lambda_{ij}\) of \(\boldsymbol{\lambda}\) represent the failure rate in the \(i\)-th time period and \(j\)-th stratum. We let the \(m\times s\) matrix \(\boldsymbol{\eta}\) represent the analogous set of dropout rates. We assume \(n\) enrollment rate periods and we let \(\boldsymbol{\gamma}\) represent an \(n\times s\) matrix of enrollment rates for each time period.

1.2.10 eEvents source code

The routine eEvents() works for a stratified as well as an unstratified population. It is accessible to users when the gsDesign package is loaded as checking is performed on input variables. The generality implied in the matrix notation just defined is available through eEvents().

# eEvents roxy [sinew] ----
#' Expected number of events for a time-to-event study
#'
#' `eEvents()` is used to calculate the expected number of events for a
#' population with a time-to-event endpoint.  It is based on calculations
#' demonstrated in Lachin and Foulkes (1986) and is fundamental in computations
#' for the sample size method they propose. Piecewise exponential survival and
#' dropout rates are supported as well as piecewise uniform enrollment. A
#' stratified population is allowed. Output is the expected number of events
#' observed given a trial duration and the above rate parameters.
#'
#' \code{eEvents()} produces an object of class \code{eEvents} with the number
#' of subjects and events for a set of pre-specified trial parameters, such as
#' accrual duration and follow-up period. The underlying power calculation is
#' based on Lachin and Foulkes (1986) method for proportional hazards assuming
#' a fixed underlying hazard ratio between 2 treatment groups. The method has
#' been extended here to enable designs to test non-inferiority. Piecewise
#' constant enrollment and failure rates are assumed and a stratified
#' population is allowed. See also \code{\link{nSurvival}} for other Lachin and
#' Foulkes (1986) methods assuming a constant hazard difference or exponential
#' enrollment rate.
#'
#' \code{print.eEvents()} formats the output for an object of class
#' \code{eEvents} and returns the input value.
#'
#' @param lambda scalar, vector or matrix of event hazard rates; rows represent
#' time periods while columns represent strata; a vector implies a single
#' stratum.
#' @param eta scalar, vector or matrix of dropout hazard rates; rows represent
#' time periods while columns represent strata; if entered as a scalar, rate is
#' constant across strata and time periods; if entered as a vector, rates are
#' constant across strata.
#' @param gamma a scalar, vector or matrix of rates of entry by time period
#' (rows) and strata (columns); if entered as a scalar, rate is constant
#' across strata and time periods; if entered as a vector, rates are constant
#' across strata.
#' @param R a scalar or vector of durations of time periods for recruitment
#' rates specified in rows of \code{gamma}. Length is the same as number of
#' rows in \code{gamma}. Note that the final enrollment period is extended as
#' long as needed.
#' @param S a scalar or vector of durations of piecewise constant event rates
#' specified in rows of \code{lambda}, \code{eta} and \code{etaE}; this is NULL
#' if there is a single event rate per stratum (exponential failure) or length
#' of the number of rows in \code{lambda} minus 1, otherwise.
#' @param T time of analysis; if \code{Tfinal=NULL}, this is also the study
#' duration.
#' @param Tfinal Study duration; if \code{NULL}, this will be replaced with
#' \code{T} on output.
#' @param minfup time from end of planned enrollment (\code{sum(R)} from output
#' value of \code{R}) until \code{Tfinal}.
#' @param x an object of class \code{eEvents} returned from \code{eEvents()}.
#' @param digits which controls number of digits for printing.
#' @param ... Other arguments that may be passed to the generic print function.
#' @return \code{eEvents()} and \code{print.eEvents()} return an object of
#' class \code{eEvents} which contains the following items: \item{lambda}{as
#' input; converted to a matrix on output.} \item{eta}{as input; converted to a
#' matrix on output.} \item{gamma}{as input.} \item{R}{as input.} \item{S}{as
#' input.} \item{T}{as input.} \item{Tfinal}{planned duration of study.}
#' \item{minfup}{as input.} \item{d}{expected number of events.}
#' \item{n}{expected sample size.} \item{digits}{as input.}
#' @examples
#'
#' # 3 enrollment periods, 3 piecewise exponential failure rates
#' str(eEvents(
#'   lambda = c(.05, .02, .01), eta = .01, gamma = c(5, 10, 20),
#'   R = c(2, 1, 2), S = c(1, 1), T = 20
#' ))
#'
#' # control group for example from Bernstein and Lagakos (1978)
#' lamC <- c(1, .8, .5)
#' n <- eEvents(
#'   lambda = matrix(c(lamC, lamC * 2 / 3), ncol = 6), eta = 0,
#'   gamma = matrix(.5, ncol = 6), R = 2, T = 4
#' )
#'
#' @aliases print.eEvents
#' @author Keaven Anderson \email{keaven_anderson@@merck.com}
#' @seealso \link{gsDesign package overview}, \link{plot.gsDesign},
#' \code{\link{gsDesign}}, \code{\link{gsHR}},
#' \code{\link{nSurvival}}
#' @references Lachin JM and Foulkes MA (1986), Evaluation of Sample Size and
#' Power for Analyses of Survival with Allowance for Nonuniform Patient Entry,
#' Losses to Follow-Up, Noncompliance, and Stratification. \emph{Biometrics},
#' 42, 507-519.
#'
#' Bernstein D and Lagakos S (1978), Sample size and power determination for
#' stratified clinical trials. \emph{Journal of Statistical Computation and
#' Simulation}, 8:65-73.
#' @keywords design
#' @rdname eEvents
#' @export
# eEvents function [sinew] ----
eEvents <- function(lambda = 1, eta = 0, gamma = 1, R = 1, S = NULL, T = 2,
                    Tfinal = NULL, minfup = 0, digits = 4) {
  if (is.null(Tfinal)) {
    if (minfup >= T) {
      stop("Minimum follow-up greater than study duration.")
    }
    Tfinal <- T
    minfupia <- minfup
  } else {
    minfupia <- max(0, minfup - (Tfinal - T))
  }

  if (!is.matrix(lambda)) {
    lambda <- matrix(lambda, nrow = length(lambda))
  }
  if (!is.matrix(eta)) {
    eta <- matrix(eta, nrow = nrow(lambda), ncol = ncol(lambda))
  }
  if (!is.matrix(gamma)) {
    gamma <- matrix(gamma, nrow = length(R), ncol = ncol(lambda))
  }
  n <- rep(0, ncol(lambda))
  d <- n
  for (i in 1:ncol(lambda))
  {
    # KA: updated following line with as.vector statements 10/16/2017
    a <- eEvents1(
      lambda = as.vector(lambda[, i]), eta = as.vector(eta[, i]),
      gamma = as.vector(gamma[, i]), R = R, S = S, T = T,
      Tfinal = Tfinal, minfup = minfup
    )
    n[i] <- a$n
    d[i] <- a$d
  }
  T1 <- cumsum(S)
  T1 <- unique(c(0, T1[T1 < T], T))
  nper <- length(T1) - 1
  names1 <- round(T1[1:nper], digits)
  namesper <- paste("-", round(T1[2:(nper + 1)], digits), sep = "")
  namesper <- paste(names1, namesper, sep = "")
  if (nper < dim(lambda)[1]) {
    lambda <- matrix(lambda[1:nper, ], nrow = nper)
  }
  if (nper < dim(eta)[1]) {
    eta <- matrix(eta[1:nper, ], nrow = nper)
  }
  rownames(lambda) <- namesper
  rownames(eta) <- namesper
  colnames(lambda) <- paste("Stratum", 1:ncol(lambda))
  colnames(eta) <- paste("Stratum", 1:ncol(eta))
  T2 <- cumsum(R)
  T2[T - T2 < minfupia] <- T - minfupia
  T2 <- unique(c(0, T2))
  nper <- length(T2) - 1
  names1 <- round(c(T2[1:nper]), digits)
  namesper <- paste("-", round(T2[2:(nper + 1)], digits), sep = "")
  namesper <- paste(names1, namesper, sep = "")
  if (nper < length(gamma)) {
    gamma <- matrix(gamma[1:nper, ], nrow = nper)
  }
  rownames(gamma) <- namesper
  colnames(gamma) <- paste("Stratum", 1:ncol(gamma))
  x <- list(
    lambda = lambda, eta = eta, gamma = gamma, R = R,
    S = S, T = T, Tfinal = Tfinal,
    minfup = minfup, d = d, n = n, digits = digits
  )
  class(x) <- "eEvents"
  return(x)
}

1.2.12 Example calculations in R

We now perform the computation from Section 1.2.5 using eEvents(), the function available for users from the gsDesign package. We begin with a single stratum with piecewise exponential failure and piecewise continuous enrollment.

eEvents(
  lambda = c(.05, .02, .01), eta = .01, gamma = c(5, 10, 20),
  R = c(2, 1, 2), S = c(1, 1), T = 20
)
Study duration:              Tfinal=20
Analysis time:                    T=20
Accrual duration:                   5
Min. end-of-study follow-up: minfup=0
Expected events (total):            11.023
Expected sample size (total):       60
Number of strata:                   1
Accrual rates:
    Stratum 1
0-2         5
2-3        10
3-5        20
Event rates:
     Stratum 1
0-1       0.05
1-2       0.02
2-20      0.01
Censoring rates:
     Stratum 1
0-1       0.01
1-2       0.01
2-20      0.01

Following is an example with multiple strata.

x <- eEvents(
  lambda = matrix(c(.05, .02, .01, .1, .04, .02), nrow = 3),
  eta = .01, gamma = matrix(c(5, 10, 20, 5, 10, 20), nrow = 3),
  R = c(2, 1, 2), S = c(1, 1), T = 20
)
x$n
[1] 60 60
x$d
[1] 11.02302 19.95135

The variables lambda, eta, and gamma provide the failure rates, censoring rates and enrollment rates, respectively. R is a scalar or vector with the duration of each enrollment rate period. The variable S specifies the duration of initial failure rate periods; should be NULL if only one failure rate period or 1 less than the number of failure rate periods, otherwise. After the initial two failure rate periods, the final failure rate is assumed to continue indefinitely. Each of these can be a scalar, vector or matrix; see the eEvents() help file for further clarification.

An additional capability for these routines is to incorporate a minimum follow-up requirement at the end of the study. We will want to use the routine for both final and interim timepoints. To appropriately incorporate the minimum follow-up at final analysis, the minimum end-of-study follow-up (minfup which defaults to 0) and the time at the end of the study (Tfinal; default is NULL) may be specified. If not specified, the final analysis time Tfinal is assumed to be the same as the input variable T and will be output as such. T is always the study time (time since beginning of enrollment) we are considering for the expected number of events. The minimum follow-up at an interim analysis may be 0 or, if the interim timing is close enough to the final, minfup - (Tfinal - T). The routine will stop accrual according to the minimum of Tfinal - minfup and the sum of values in R; this is done without changing the input values of either on return from the routine. As a simple example, we consider the following modification of the above.

eEvents(
  lambda = c(.05, .02, .01), eta = .01, gamma = c(5, 10, 20),
  R = c(2, 1, 20), S = c(1, 1), T = 18, Tfinal = 22, minfup = 6
)
Study duration:              Tfinal=22
Analysis time:                    T=18
Accrual duration:                   16
Min. end-of-study follow-up: minfup=6
Expected events (total):            35.2387
Expected sample size (total):       280
Number of strata:                   1
Accrual rates:
     Stratum 1
0-2          5
2-3         10
3-16        20
Event rates:
     Stratum 1
0-1       0.05
1-2       0.02
2-18      0.01
Censoring rates:
     Stratum 1
0-1       0.01
1-2       0.01
2-18      0.01

1.2.13 Literature example: stratified with exponential failure, no dropouts

Bernstein and Lagakos (1978) provided FORTRAN code and examples for a stratified sample size, but did not provide for random dropouts and only accounted for exponential failure rates. This is equivalent to \(\eta=0\) for the more general case here. Their example assumed 2 years of accrual, 2 years of follow-up, 3 strata with 40%, 40% and 20% of enrollment, equal randomization between 2 arms (control and experimental), control group failure rates in the three strata of 1, .8 and .5, a hazard ratio of 2/3 (experimental to control). The expected proportion of deaths by treatment group and stratum are noted as .941, .899, .767 (control) and .854, .788, .625 (experimental). To have an expected number of enrollees of 1 in each group, we require an enrollment rate of .5 over the 2 years, or \(\gamma=.5\). We calculate these numbers as follows

lamC <- c(1, .8, .5)
n <- eEvents(
  lambda = matrix(c(lamC, lamC * 2 / 3), ncol = 6), eta = 0,
  gamma = matrix(.5, ncol = 6), R = 2, T = 4
)
n$d
[1] 0.9414902 0.8992911 0.7674558 0.8544147 0.7883950 0.6252700

Now for each patient enrolled in the study overall, we compute the expected number events per treatment group given the above enrollment proportions and event rates. First, we do this calculations using the above computation. Then, we follow with direct computation from eEvents() by setting the rates so there is an expected enrollment of 1 patient.

sum(n$d[1:3] * c(.4, .4, .2)) / 2
[1] 0.4449018
sum(n$d[4:6] * c(.4, .4, .2)) / 2
[1] 0.391089
n <- eEvents(
  lambda = matrix(c(lamC, lamC * 2 / 3), ncol = 6), eta = 0,
  gamma = matrix(.25 * c(.4, .4, .2, .4, .4, .2), ncol = 6), R = 2, T = 4
)
c(sum(n$d[1:3]), sum(n$d[4:6]))
[1] 0.4449018 0.3910890