2 Sample size for a fixed design

There are two approaches to sizing a study from the literature that are presented here. Both cases will be developed for cases where time-to-event distribution assumptions include a stratified population, proportional hazards and piecewise-exponential failure rates.

The Lachin and Foulkes (1986) approach where enrollment and study duration are fixed while the enrollment rates are changed proportionately to obtain sufficient enrollment and endpoints to ensure the desired power.
The Kim and Tsiatis (1990) method where enrollment rates are fixed and duration of enrollment is varied to ensure power.

2.1 Notation

For the remainder of this document we consider two treatment groups: an experimental group with parameters denoted with a subscript \(E\) and a control group with parameters denoted with a subscript \(C\). We will work just with piecewise exponential failure rates and piecewise continuous enrollment. The proportion of patients randomized to the experimental group will be denoted by \(\xi\). Recall that we have used \(g()\) to denote the enrollment rate over time. We will now assume the enrollment rate in the experimental group is \(\xi g()\), while \((1-\xi)g()\) is the enrollment rate in the control group. For piecewise constant enrollment we assumed \(g(t)=\gamma_i\) for \(t\) in \((t_{i-1}, t_i]\), for \(i\) from 1 to \(m\) and \(0=t_0<t_1\ldots<t_m\).

We will use matrix to notation to abbeviate many of the assumptions for each treatment group.

2.2 The basic equation

The two methods will be developed in separate functions initially. Then a single function nSurv() will be developed that is intended for to be made visible to users from the gsDesign package.

2.3 Asymptotic approximation for power and sample size

We will use the formulation of Lachin and Foulkes (1986) and their equation (2.4) for the asymptotic approximation relating power to sample size for a proportional hazards model. We generalize their results for the case when the null hypothesis does not specify equal hazards. Denote the null hypothesis failure rates for control and experimental treatment groups as \(\lambda_{00}\) and \(\lambda_{01}\), respectively. Denote the alternate hypothesis rates as \(\lambda_{10}\) and \(\lambda_{11}\). Further, denote the alternate hypothesis hazard ratio \(h_1=\lambda_{11}/\lambda_{10}\), and the null hypothesis hazard ratio \(h_0=\lambda_{01}/\lambda_{00}\). We let censoring rates be specific to the control (\(\eta_0\)) and experimental (\(\eta_1\)) groups; these values are only implicit in the equations below. Lachin and Foulkes assumed a null hypothesis with no difference between failure rates in the control and experimental rates and test for superiority. That is, \(\lambda_{00}=\lambda_{01}\) (\(h_0=1\)) and \(\lambda_{10}<\lambda_{01}\) (\(h_1<1\)). They set event rates under the null hypothesis so that the the weighted average event rate is the same under the null and alternate hypotheses: \[\begin{equation} \lambda_{00}=\lambda_{01}=\bar\lambda=(1-\xi)\lambda_{10}+\xi\lambda_{11}.\label{eq:lambdabar} \end{equation}\] The apparent intent of this is to equalize the variance for the log hazard ratio under null and alternative hypotheses. While this will not be exactly the case, we will check the expected number of events under the null and alternate hypotheses later to see that the approximation works reasonably well. The Lachin and Foulkes power equation for proportional hazards translates in our notation to: \[\begin{eqnarray} \sqrt{N}\ln(HR) &=& Z_\alpha \sqrt{E\left\{\delta|\bar{\lambda},\eta\right\}^{-1}(\xi^{-1}+(1-\xi)^{-1})} \nonumber\\ &&+Z_\beta\sqrt{ E\left\{\delta|\lambda_1,\eta\right\}^{-1}\xi^{-1}+ E\left\{\delta|\lambda_0,\eta\right\}^{-1}(1-\xi)^{-1}} \label{eq:npwr} \end{eqnarray}\] Lachin and Foulkes did not cover any cases other than equality under the null hypothesis; (i.e., the assumed \(h_0 \neq 1\) or, equivalently, \(\lambda_{00}\neq \lambda_{01}\)). Equation (\(\ref{eq:npwr}\)) generalizes in this case to \[\begin{eqnarray} \sqrt{N}\ln\left(\frac{h_1}{h_0}\right) &=& Z_\alpha \sqrt{E\left\{\delta|\lambda_{01},\eta_1\right\}^{-1}\xi^{-1}+ E\left\{\delta|\lambda_{00},\eta_0\right\}^{-1}(1-\xi)^{-1}} \nonumber\\ &&+Z_\beta\sqrt{ E\left\{\delta|\lambda_{11},\eta_1\right\}^{-1}\xi^{-1}+ E\left\{\delta|\lambda_{10},\eta_0\right\}^{-1}(1-\xi)^{-1}} \label{eq:npwra} \end{eqnarray}\] Note that the computations inside square root signs are variance approximations. When multiple strata are considered, we sum the variance approximations across strata; we will assume a common hazard ratio across strata - although it would appear that power for testing for a weighted average of the log hazard ratios across strata might be considered with the same approach. With version 3.x.x we changed the weighting across strata and the piecewise exponential failure windows to inverse-variance weighted sum. [IN DEVELOPMENT!!!]. To apply equation (\(\ref{eq:npwra}\)), we need specific rates under both the null and alternate hypotheses. For superiority above we started with alternate hypothesis rates and computed rates for the null hypothesis. This also seems appropriate when testing for non-inferiority. Since \[\lambda_{01}=h_0\lambda_{00},\] the average of the null hypothesis rates weighted by randomization ratio of experimental to control \(r\) is \[\bar\lambda=\frac{\lambda_{00}+r\lambda_{01}}{1+r}=\frac{\lambda_{00}(1+h_0r)}{1+r}.\] We have the analogous equation under the alternate hypothesis and we will assume the same \(\bar\lambda\): \[\bar\lambda=\frac{\lambda_{10}(1+h_1r)}{1+r}.\] Given the fixed values for \(h_0\), \(h_1\), \(r\) and \(\lambda_{10}\) we can solve for \[\lambda_{00}=\frac{1+h_1r}{1+h_0r}\lambda_{10}\] and \[\lambda_{10}=h_0\lambda_{00}.\]

2.4 A period naming utility

We will want to make it clear upon return from any sample size computation what enrollment periods and failure rate periods were assumed. A simple period naming function will make this possible. We define the function nameperiod() and provide an initial example:

# nameperiod function [sinew] ----
nameperiod <- function(R, digits = 2) {
  if (length(R) == 1) {
    return(paste("0-", round(R, digits), sep = ""))
  }
  R0 <- c(0, R[1:(length(R) - 1)])
  return(paste(round(R0, digits), "-", round(R, digits), sep = ""))
}

Now we consider enrollment period durations, convert these to cumulative durations and name them.

R <- c(1, 2, 3)
gamma <- matrix(c(4, 5, 6), nrow = 3)
rownames(gamma) <- nameperiod(cumsum(R))
gamma

2.5 Fixed enrollment and study duration

Equation (\(\ref{eq:npwra}\)) can be solved directly given the values of \(E\{\delta|\lambda,\eta\}\) for the various values of \(\lambda\) and \(\eta\). That is, these values do not change with proportional increases or decreases in enrollment rates over a fixed duration of time.

2.5.1 LFPWE source

Note that on output, the input enrollment duration R is altered to ensure T - minfup = sum(R). This may eliminate some enrollment periods and/or alter the final returned value of R. The returned value of gamma is also checked so that the length of the returned value matches the length of the returned R.

# LFPWE function [sinew] ----
LFPWE <- function(alpha = .025, sided = 1, beta = .1,
                  lambdaC = log(2) / 6, hr = .5, hr0 = 1, etaC = 0, etaE = 0,
                  gamma = 1, ratio = 1, R = 18, S = NULL, T = 24, minfup = NULL) {
  # set up parameters
  zalpha <- -stats::qnorm(alpha / sided)
  zbeta <- -stats::qnorm(beta)
  if (is.null(minfup)) minfup <- max(0, T - sum(R))
  if (length(R) == 1) {
    R <- T - minfup
  } else if (sum(R) != T - minfup) {
    cR <- cumsum(R)
    nR <- length(R)
    if (cR[length(cR)] < T - minfup) {
      cR[length(cR)] <- T - minfup
    } else {
      cR[cR > T - minfup] <- T - minfup
      cR <- unique(cR)
    }
    if (length(cR) > 1) {
      R <- cR - c(0, cR[1:(length(cR) - 1)])
    } else {
      R <- cR
    }
    if (nR != length(R)) {
      if (is.vector(gamma)) {
        gamma <- gamma[1:length(R)]
      } else {
        gamma <- gamma[1:length(R), ]
      }
    }
  }
  ngamma <- length(R)
  if (is.null(S)) {
    nlambda <- 1
  } else {
    nlambda <- length(S) + 1
  }
  Qe <- ratio / (1 + ratio)
  Qc <- 1 - Qe

  # compute H0 failure rates as average of control, experimental
  if (length(ratio) == 1) {
    lambdaC0 <- (1 + hr * ratio) / (1 + hr0 * ratio) * lambdaC
    gammaC <- gamma * Qc
    gammaE <- gamma * Qe
  } else {
    lambdaC0 <- lambdaC %*% diag((1 + hr * ratio) / (1 + hr0 * ratio))
    gammaC <- gamma %*% diag(Qc)
    gammaE <- gamma %*% diag(Qe)
  }
  # do computations
  eDC0 <- eEvents(
    lambda = lambdaC0, eta = etaC, gamma = gammaC,
    R = R, S = S, T = T, minfup = minfup
  )$d
  eDE0 <- eEvents(
    lambda = lambdaC0 * hr0, eta = etaE, gamma = gammaE,
    R = R, S = S, T = T, minfup = minfup
  )$d
  eDC <- eEvents(
    lambda = lambdaC, eta = etaC, gamma = gammaC,
    R = R, S = S, T = T, minfup = minfup
  )
  eDE <- eEvents(
    lambda = lambdaC * hr, eta = etaE, gamma = gammaE,
    R = R, S = S, T = T, minfup = minfup
  )

  n <- ((zalpha * sqrt(1 / sum(eDC0) + 1 / sum(eDE0)) +
    zbeta * sqrt(1 / sum(eDC$d) + 1 / sum(eDE$d))
  ) / log(hr / hr0))^2

  mx <- sum(eDC$n + eDE$n)
  rval <- list(
    alpha = alpha, sided = sided, beta = beta, power = 1 - beta,
    lambdaC = lambdaC, etaC = etaC, etaE = etaE, gamma = n * gamma,
    ratio = ratio, R = R, S = S, T = T, minfup = minfup,
    hr = hr, hr0 = hr0, n = n * mx, d = n * sum(eDC$d + eDE$d),
    eDC = eDC$d * n, eDE = eDE$d * n, eDC0 = eDC0 * n, eDE0 = eDE0 * n,
    eNC = eDC$n * n, eNE = eDE$n * n, variable = "Accrual rate"
  )
  class(rval) <- "nSurv"
  return(rval)
}

# print.nSurv roxy [sinew] ----
#' @rdname nSurv
#' @export
# print.nSurv function [sinew] ----
print.nSurv <- function(x, digits = 4, ...) {
  if (!inherits(x, "nSurv")) {
    stop("Primary argument must have class nSurv")
  }
  x$digits <- digits
  x$sided <- 1
  cat("Fixed design, two-arm trial with time-to-event\n")
  cat("outcome (Lachin and Foulkes, 1986).\n")
  cat("Solving for: ", x$variable, "\n")
  cat("Hazard ratio                  H1/H0=",
    round(x$hr, digits),
    "/", round(x$hr0, digits), "\n",
    sep = ""
  )
  cat("Study duration:                   T=",
    round(x$T, digits), "\n",
    sep = ""
  )
  cat("Accrual duration:                   ",
    round(x$T - x$minfup, digits), "\n",
    sep = ""
  )
  cat("Min. end-of-study follow-up: minfup=",
    round(x$minfup, digits), "\n",
    sep = ""
  )
  cat("Expected events (total, H1):        ",
    round(x$d, digits), "\n",
    sep = ""
  )
  cat("Expected sample size (total):       ",
    round(x$n, digits), "\n",
    sep = ""
  )
  enrollper <- periods(x$S, x$T, x$minfup, x$digits)
  cat("Accrual rates:\n")
  print(round(x$gamma, digits))
  cat("Control event rates (H1):\n")
  print(round(x$lambda, digits))
  if (max(abs(x$etaC - x$etaE)) == 0) {
    cat("Censoring rates:\n")
    print(round(x$etaC, digits))
  } else {
    cat("Control censoring rates:\n")
    print(round(x$etaC, digits))
    cat("Experimental censoring rates:\n")
    print(round(x$etaE, digits))
  }
  cat("Power:                 100*(1-beta)=",
    round((1 - x$beta) * 100, digits), "%\n",
    sep = ""
  )
  cat("Type I error (", x$sided,
    "-sided):   100*alpha=",
    100 * x$alpha, "%\n",
    sep = ""
  )
  if (min(x$ratio == 1) == 1) {
    cat("Equal randomization:          ratio=1\n")
  } else {
    cat(
      "Randomization (Exp/Control):  ratio=",
      x$ratio, "\n"
    )
  }
}

2.5.2 Experimental version

# LFPWEnew function [sinew] ----
LFPWEnew <- function(alpha = .025, sided = 1, beta = .1,
                     lambdaC = log(2) / 6, hr = .5, hr0 = 1, etaC = 0, etaE = 0,
                     gamma = 1, ratio = 1, R = 18, S = NULL, T = 24, minfup = NULL) {
  # set up parameters
  zalpha <- -stats::qnorm(alpha / sided)
  if (is.null(minfup)) minfup <- max(0, T - sum(R))
  if (length(R) == 1) {
    R <- T - minfup
  } else if (sum(R) != T - minfup) {
    cR <- cumsum(R)
    nR <- length(R)
    if (cR[length(cR)] < T - minfup) {
      cR[length(cR)] <- T - minfup
    } else {
      cR[cR > T - minfup] <- T - minfup
      cR <- unique(cR)
    }
    if (length(cR) > 1) {
      R <- cR - c(0, cR[1:(length(cR) - 1)])
    } else {
      R <- cR
    }
    if (nR != length(R)) {
      if (is.vector(gamma)) {
        gamma <- gamma[1:length(R)]
      } else {
        gamma <- gamma[1:length(R), ]
      }
    }
  }
  ngamma <- length(R)
  if (is.null(S)) {
    nlambda <- 1
  } else {
    nlambda <- length(S) + 1
  }
  Qe <- ratio / (1 + ratio)
  Qc <- 1 - Qe

  # compute H0 failure rates as average of control, experimental
  if (length(ratio) == 1) {
    lambdaC0 <- (1 + hr * ratio) / (1 + hr0 * ratio) * lambdaC
    gammaC <- gamma * Qc
    gammaE <- gamma * Qe
  } else {
    lambdaC0 <- lambdaC %*% diag((1 + hr * ratio) / (1 + hr0 * ratio))
    gammaC <- gamma %*% diag(Qc)
    gammaE <- gamma %*% diag(Qe)
  }
  # do computations
  eDC0 <- eEvents(
    lambda = lambdaC0, eta = etaC, gamma = gammaC,
    R = R, S = S, T = T, minfup = minfup
  )$d
  eDE0 <- eEvents(
    lambda = lambdaC0 * hr0, eta = etaE, gamma = gammaE,
    R = R, S = S, T = T, minfup = minfup
  )$d
  eDC <- eEvents(
    lambda = lambdaC, eta = etaC, gamma = gammaC,
    R = R, S = S, T = T, minfup = minfup
  )
  eDE <- eEvents(
    lambda = lambdaC * hr, eta = etaE, gamma = gammaE,
    R = R, S = S, T = T, minfup = minfup
  )
  # sample size
  if (!is.null(beta)) {
    zbeta <- -stats::qnorm(beta)
    ##### OLD
    n_old <- ((zalpha * sqrt(1 / sum(eDC0) + 1 / sum(eDE0)) +
      zbeta * sqrt(1 / sum(eDC$d) + 1 / sum(eDE$d))
    ) / log(hr / hr0))^2
    cat("Control events (H0): ")
    cat(eDC0)
    cat("\nControl events (H1): ")
    cat(eDC$d)
    cat("\nExperimental events (H0): ")
    cat(eDE0)
    cat("\nExperimental events (H1): ")
    cat(eDE$d)
    cat("\nOld n: ")
    cat(n_old)
    cat("\n")
    #####
    ##### NEW 20230327; NO CHANGE FOR UNSTRATIFIED
    n <- ((zalpha * sqrt(1 / sum((1 / eDC0 + 1 / eDE0)^(-1))) +
      zbeta * sqrt(1 / sum((1 / eDC$d + 1 / eDE$d)^(-1)))
    ) / log(hr / hr0))^2
    cat("New n: ")
    cat(n)
    cat("\n")
  } else {
    # power
  }
  mx <- sum(eDC$n + eDE$n)
  rval <- list(
    alpha = alpha, sided = sided, beta = beta, power = 1 - beta,
    lambdaC = lambdaC, etaC = etaC, etaE = etaE, gamma = n * gamma,
    ratio = ratio, R = R, S = S, T = T, minfup = minfup,
    hr = hr, hr0 = hr0, n = n * mx, d = n * sum(eDC$d + eDE$d),
    eDC = eDC$d * n, eDE = eDE$d * n, eDC0 = eDC0 * n, eDE0 = eDE0 * n,
    eNC = eDC$n * n, eNE = eDE$n * n, variable = "Accrual rate"
  )
  class(rval) <- "nSurv"
  return(rval)
}

2.5.3 Examples

In the Merck (ELSTIC) guidance on time-to-event sample size, an example with \(\alpha=.025\), \(\beta=.1\), \(\lambda_C=.2\), \(\lambda_E=.1\), \(\eta_C=\eta_E=.1\), enrollment period of 1/2 year and total study time of 2 years. The Lachin and Foulkes (1986) method is documented to produce a sample size of 430, which is the same as the following if rounded up to the nearest even number. This is reproduced here by both the nSurvival() routine, which has been validated, and by LFPWE().

LFPWE(
  alpha = .025, sided = 1, beta = .1, lambdaC = .2, hr = 1 / 2,
  etaC = .1, etaE = .1, gamma = 1, ratio = 1, R = .5, S = NULL, T = 2
)

Fixed design, two-arm trial with time-to-event
outcome (Lachin and Foulkes, 1986).
Solving for:  Accrual rate 
Hazard ratio                  H1/H0=0.5/1
Study duration:                   T=2
Accrual duration:                   0.5
Min. end-of-study follow-up: minfup=1.5
Expected events (total, H1):        90.0987
Expected sample size (total):       429.6189
Accrual rates:
[1] 859.2377
Control event rates (H1):
[1] 0.2
Censoring rates:
[1] 0.1
Power:                 100*(1-beta)=90%
Type I error (1-sided):   100*alpha=2.5%
Equal randomization:          ratio=1

x <- gsDesign:::nSurvival(lambda1 = .2, lambda2 = .1, eta = .1, Tr = .5, Ts = 2)
x$n

[1] 429.6189

x$nEvents

[1] 90.09875

Fixed design, two-arm trial with time-to-event
outcome (Lachin and Foulkes, 1986).
Study duration (fixed):          Ts=2
Accrual duration (fixed):        Tr=0.5
Uniform accrual:              entry="unif"
Control median:      log(2)/lambda1=3.5
Experimental median: log(2)/lambda2=6.9
Censoring median:        log(2)/eta=6.9
Control failure rate:       lambda1=0.2
Experimental failure rate:  lambda2=0.1
Censoring rate:                 eta=0.1
Power:                 100*(1-beta)=90%
Type I error (1-sided):   100*alpha=2.5%
Equal randomization:          ratio=1
Sample size based on hazard ratio=0.5 (type="rr")
Sample size (computed):           n=430
Events required (computed): nEvents=91

Same example with new version:

LFPWEnew(
  alpha = .025, sided = 1, beta = .1, lambdaC = .2, hr = 1 / 2,
  etaC = .1, etaE = .1, gamma = 1, ratio = 1, R = .5, S = NULL, T = 2
)

We return to the Bernstein and Lagakos (1978) example of a stratified trial with exponential failure times. Their methods incorporate only a variance estimate under the alternate hypothesis and yields a sample size of 172, slightly smaller than what we compute here.

x <- LFPWE(
  alpha = .05, beta = .2, hr = 2 / 3, ratio = 1,
  R = 2, S = NULL, T = 4,
  lambdaC = matrix(c(1, .8, .5), nrow = 1),
  etaC = matrix(0, nrow = 1, ncol = 3),
  etaE = matrix(0, nrow = 1, ncol = 3),
  gamma = matrix(c(.4, .4, .2), nrow = 1)
)
x

Fixed design, two-arm trial with time-to-event
outcome (Lachin and Foulkes, 1986).
Solving for:  Accrual rate 
Hazard ratio                  H1/H0=0.6667/1
Study duration:                   T=4
Accrual duration:                   2
Min. end-of-study follow-up: minfup=2
Expected events (total, H1):        149.455
Expected sample size (total):       178.7758
Accrual rates:
        [,1]    [,2]    [,3]
[1,] 35.7552 35.7552 17.8776
Control event rates (H1):
     [,1] [,2] [,3]
[1,]    1  0.8  0.5
Censoring rates:
     [,1] [,2] [,3]
[1,]    0    0    0
Power:                 100*(1-beta)=80%
Type I error (1-sided):   100*alpha=5%
Equal randomization:          ratio=1

Next, we consider the piecewise exponential case

x <- LFPWE(
  alpha = .025, beta = .1,
  lambdaC = matrix(c(.05, .02, .01), ncol = 1), hr = .6,
  etaC = matrix(.01, nrow = 3, ncol = 1), etaE = matrix(.01, nrow = 3, ncol = 1),
  gamma = matrix(c(5, 10, 20), ncol = 1),
  R = c(2, 1, 2), S = c(1, 1), T = 20
)
x

Fixed design, two-arm trial with time-to-event
outcome (Lachin and Foulkes, 1986).
Solving for:  Accrual rate 
Hazard ratio                  H1/H0=0.6/1
Study duration:                   T=20
Accrual duration:                   5
Min. end-of-study follow-up: minfup=15
Expected events (total, H1):        164.1408
Expected sample size (total):       1099.533
Accrual rates:
         [,1]
[1,]  91.6277
[2,] 183.2555
[3,] 366.5109
Control event rates (H1):
     [,1]
[1,] 0.05
[2,] 0.02
[3,] 0.01
Censoring rates:
     [,1]
[1,] 0.01
[2,] 0.01
[3,] 0.01
Power:                 100*(1-beta)=90%
Type I error (1-sided):   100*alpha=2.5%
Equal randomization:          ratio=1

Finally, we have all of: piecewise exponential failure, piecewise constant enrollment, and a stratified population with differing randomization ratio by stratum.

x <- LFPWE(
  lambdaC = log(2) / matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3),
  gamma = matrix(c(2, 4, 8, 3, 6, 10), nrow = 2),
  hr = .6, S = c(3, 6), ratio = 1:3, R = c(3, 3), minfup = 6, T = 30
)
x

Fixed design, two-arm trial with time-to-event
outcome (Lachin and Foulkes, 1986).
Solving for:  Accrual rate 
Hazard ratio                  H1/H0=0.6/1
Study duration:                   T=30
Accrual duration:                   24
Min. end-of-study follow-up: minfup=6
Expected events (total, H1):        171.4585
Expected sample size (total):       287.077
Accrual rates:
       [,1]   [,2]   [,3]
[1,] 1.4177 5.6707 4.2530
[2,] 2.8353 2.1265 7.0883
Control event rates (H1):
       [,1]   [,2]   [,3]
[1,] 0.2310 0.1155 0.0770
[2,] 0.1733 0.0866 0.0578
[3,] 0.1386 0.0693 0.0462
Censoring rates:
[1] 0
Power:                 100*(1-beta)=90%
Type I error (1-sided):   100*alpha=2.5%
Randomization (Exp/Control):  ratio= 1 2 3

2.6 Fixed enrollment rate methods

Now we assume fixed enrollment rates and adjust either the enrollment duration or follow-up duration to obtain a targeted power or number of events. As this is following the strategy of Kim and Tsiatis (1990), we name the function KT(). KT() is primarily a root-finding program based on the routine KTZ() which is used to calculate power based on (eq:npwr) or the expected number of events under the alternate hypothesis (using eEvents()) given enrollment rates, duration of a study and duration of follow-up.

2.6.1 KTZ source code and discussion

KTZ() can compute the difference from a targeted number of events or targeted power. It can also compute power for a trial with given design characteristics. Depending on the structure of other input variables, a non-NULL input value of x may be used to set

the minimum follow-up after a given enrollment duration (if minfup is input as NULL), or
the duration of enrollment.

If simple is input as TRUE, then the routine is used to return a difference from a target value for

a targeted total number of events in n1Target if n1Target is not NULL, or
a targeted power (target is actually set to the standard normal inverse of 1 - beta), if beta if both n1Target and beta are NULL on input.

If simple is input as FALSE, then an nSurv object is returned giving the power and type II error for the input parameter values for the trial.

# KTZ function [sinew] ----
KTZ <- function(x = NULL, minfup = NULL, n1Target = NULL,
                lambdaC = log(2) / 6, etaC = 0, etaE = 0,
                gamma = 1, ratio = 1, R = 18, S = NULL, beta = .1,
                alpha = .025, sided = 1, hr0 = 1, hr = .5, simple = TRUE) {
  zalpha <- -stats::qnorm(alpha / sided)
  Qc <- 1 / (1 + ratio)
  Qe <- 1 - Qc
  # set minimum follow-up to x if that is missing and x is given
  if (!is.null(x) && is.null(minfup)) {
    minfup <- x
    if (sum(R) == Inf) {
      stop("If minimum follow-up is sought, enrollment duration must be finite")
    }
    T <- sum(R) + minfup
    variable <- "Follow-up duration"
  } else if (!is.null(x) && !is.null(minfup)) { # otherwise, if x is given, set it to accrual duration
    T <- x + minfup
    R[length(R)] <- Inf
    variable <- "Accrual duration"
  } else { # otherwise, set follow-up time to accrual plus follow-up
    T <- sum(R) + minfup
    variable <- "Power"
  }
  # compute H0 failure rates as average of control, experimental
  if (length(ratio) == 1) {
    lambdaC0 <- (1 + hr * ratio) / (1 + hr0 * ratio) * lambdaC
    gammaC <- gamma * Qc
    gammaE <- gamma * Qe
  } else {
    lambdaC0 <- lambdaC %*% diag((1 + hr * ratio) / (1 + hr0 * ratio))
    gammaC <- gamma %*% diag(Qc)
    gammaE <- gamma %*% diag(Qe)
  }

  # do computations
  eDC <- eEvents(
    lambda = lambdaC, eta = etaC, gamma = gammaC,
    R = R, S = S, T = T, minfup = minfup
  )
  eDE <- eEvents(
    lambda = lambdaC * hr, eta = etaE, gamma = gammaE,
    R = R, S = S, T = T, minfup = minfup
  )
  # if this is all that is needed, return difference
  # from targeted number of events
  if (simple && !is.null(n1Target)) {
    return(sum(eDC$d + eDE$d) - n1Target)
  }
  eDC0 <- eEvents(
    lambda = lambdaC0, eta = etaC, gamma = gammaC,
    R = R, S = S, T = T, minfup = minfup
  )
  eDE0 <- eEvents(
    lambda = lambdaC0 * hr0, eta = etaE, gamma = gammaE,
    R = R, S = S, T = T, minfup = minfup
  )
  # compute Z-value related to power from power equation
  zb <- (log(hr0 / hr) -
    zalpha * sqrt(1 / sum(eDC0$d) + 1 / sum(eDE0$d))) /
    sqrt(1 / sum(eDC$d) + 1 / sum(eDE$d))
  # if that is all that is needed, return difference from
  # targeted value
  if (simple) {
    if (!is.null(beta)) {
      return(zb + stats::qnorm(beta))
    } else {
      return(stats::pnorm(-zb))
    }
  }
  # compute power
  power <- stats::pnorm(zb)
  beta <- 1 - power
  # set accrual period durations
  if (sum(R) != T - minfup) {
    if (length(R) == 1) {
      R <- T - minfup
    } else {
      nR <- length(R)
      cR <- cumsum(R)
      cR[cR > T - minfup] <- T - minfup
      cR <- unique(cR)
      cR[length(R)] <- T - minfup
      if (length(cR) == 1) {
        R <- cR
      } else {
        R <- cR - c(0, cR[1:(length(cR) - 1)])
      }
      if (length(R) != nR) {
        gamma <- matrix(gamma[1:length(R), ], nrow = length(R))
        gdim <- dim(gamma)
      }
    }
  }
  rval <- list(
    alpha = alpha, sided = sided, beta = beta, power = power,
    lambdaC = lambdaC, etaC = etaC, etaE = etaE,
    gamma = gamma, ratio = ratio, R = R, S = S, T = T,
    minfup = minfup, hr = hr, hr0 = hr0, n = sum(eDC$n + eDE$n),
    d = sum(eDC$d + eDE$d), tol = NULL, eDC = eDC$d, eDE = eDE$d,
    eDC0 = eDC0$d, eDE0 = eDE0$d, eNC = eDC$n, eNE = eDE$n,
    variable = variable
  )
  class(rval) <- "nSurv"
  return(rval)
}

# input minimum follow-up in x and consider R fixed
KTZ(
  x = 2, lambdaC = log(2) / matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3),
  gamma = matrix(c(2, 4, 8, 3, 6, 10), nrow = 2), n1Target = 0,
  hr = .6, S = c(3, 6), R = c(3, 15), minfup = NULL
)

[1] 148.8083

# set total time to x + minfup (R is variable)
KTZ(
  x = 10, lambdaC = log(2) / matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3),
  gamma = matrix(c(2, 4, 8, 3, 6, 10), nrow = 2), n1Target = 0,
  hr = .6, S = c(3, 6), R = c(3, 15), minfup = 2
)

[1] 65.46827

# same as 2nd example with differing randomization ratio by stratum
KTZ(
  x = 10, lambdaC = log(2) / matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3),
  gamma = matrix(c(2, 4, 8, 3, 6, 10), nrow = 2), n1Target = 0,
  hr = .6, S = c(3, 6), ratio = 1:3, R = c(3, 15), minfup = 2
)

[1] 61.67763

2.6.2 KT source code

# KT function [sinew] ----
KT <- function(alpha = .025, sided = 1, beta = .1,
               lambdaC = log(2) / 6, hr = .5, hr0 = 1, etaC = 0, etaE = 0,
               gamma = 1, ratio = 1, R = 18, S = NULL, minfup = NULL,
               n1Target = NULL, tol = .Machine$double.eps^0.25) {
  # set up parameters
  ngamma <- length(R)
  if (is.null(S)) {
    nlambda <- 1
  } else {
    nlambda <- length(S) + 1
  }
  Qe <- ratio / (1 + ratio)
  Qc <- 1 - Qe
  if (!is.matrix(lambdaC)) lambdaC <- matrix(lambdaC)
  ldim <- dim(lambdaC)
  nstrata <- ldim[2]
  nlambda <- ldim[1]
  etaC <- matrix(etaC, nrow = nlambda, ncol = nstrata)
  etaE <- matrix(etaE, nrow = nlambda, ncol = nstrata)
  if (!is.matrix(gamma)) gamma <- matrix(gamma)
  gdim <- dim(gamma)
  eCdim <- dim(etaC)
  eEdim <- dim(etaE)

  # search for trial duration needed to achieve desired power
  if (is.null(minfup)) {
    if (sum(R) == Inf) {
      stop("Enrollment duration must be specified as finite")
    }
    left <- KTZ(.01,
      lambdaC = lambdaC, n1Target = n1Target,
      etaC = etaC, etaE = etaE, gamma = gamma, ratio = ratio,
      R = R, S = S, beta = beta, alpha = alpha, sided = sided,
      hr0 = hr0, hr = hr
    )
    right <- KTZ(1000,
      lambdaC = lambdaC, n1Target = n1Target,
      etaC = etaC, etaE = etaE, gamma = gamma, ratio = ratio,
      R = R, S = S, beta = beta, alpha = alpha, sided = sided,
      hr0 = hr0, hr = hr
    )
    if (left > 0) stop("Enrollment duration over-powers trial")
    if (right < 0) stop("Enrollment duration insufficient to power trial")
    y <- stats::uniroot(
      f = KTZ, interval = c(.01, 10000), lambdaC = lambdaC,
      etaC = etaC, etaE = etaE, gamma = gamma, ratio = ratio,
      R = R, S = S, beta = beta, alpha = alpha, sided = sided,
      hr0 = hr0, hr = hr, tol = tol, n1Target = n1Target
    )
    minfup <- y$root
    xx <- KTZ(
      x = y$root, lambdaC = lambdaC,
      etaC = etaC, etaE = etaE, gamma = gamma, ratio = ratio,
      R = R, S = S, minfup = NULL, beta = beta, alpha = alpha,
      sided = sided, hr0 = hr0, hr = hr, simple = F
    )
    xx$tol <- tol
    return(xx)
  } else {
    y <- stats::uniroot(
      f = KTZ, interval = minfup + c(.01, 10000), lambdaC = lambdaC,
      n1Target = n1Target,
      etaC = etaC, etaE = etaE, gamma = gamma, ratio = ratio,
      R = R, S = S, minfup = minfup, beta = beta,
      alpha = alpha, sided = sided, hr0 = hr0, hr = hr, tol = tol
    )
    xx <- KTZ(
      x = y$root, lambdaC = lambdaC,
      etaC = etaC, etaE = etaE, gamma = gamma, ratio = ratio,
      R = R, S = S, minfup = minfup, beta = beta, alpha = alpha,
      sided = sided, hr0 = hr0, hr = hr, simple = F
    )
    xx$tol <- tol
    return(xx)
  }
}

# set minimum follow-up (allow R to vary to get power)
KT(
  lambdaC = log(2) / matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3),
  gamma = matrix(c(2, 4, 8, 3, 6, 10), nrow = 2), n1Target = 10,
  hr = .6, S = c(3, 6), ratio = 1:3, R = c(3, 15), minfup = 2
)

Fixed design, two-arm trial with time-to-event
outcome (Lachin and Foulkes, 1986).
Solving for:  Accrual duration 
Hazard ratio                  H1/H0=0.6/1
Study duration:                   T=4.6979
Accrual duration:                   2.6979
Min. end-of-study follow-up: minfup=2
Expected events (total, H1):        10
Expected sample size (total):       43.1662

Warning in max(periods): no non-missing arguments to max; returning -Inf

Accrual rates:
     [,1] [,2] [,3]
[1,]    2    8    6
[2,]    4    3   10
Control event rates (H1):
       [,1]   [,2]   [,3]
[1,] 0.2310 0.1155 0.0770
[2,] 0.1733 0.0866 0.0578
[3,] 0.1386 0.0693 0.0462
Censoring rates:
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0
Power:                 100*(1-beta)=10.9406%
Type I error (1-sided):   100*alpha=2.5%
Randomization (Exp/Control):  ratio= 1 2 3

# allow minimum follow-up to vary to obtain power
KT(
  lambdaC = log(2) / matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3),
  gamma = matrix(c(2, 4, 8, 3, 6, 10), nrow = 2), n1Target = 100,
  hr = .6, S = c(3, 6), ratio = 1:3, R = c(3, 10), minfup = NULL
)

Fixed design, two-arm trial with time-to-event
outcome (Lachin and Foulkes, 1986).
Solving for:  Follow-up duration 
Hazard ratio                  H1/H0=0.6/1
Study duration:                   T=16.7158
Accrual duration:                   13
Min. end-of-study follow-up: minfup=3.7158
Expected events (total, H1):        100.0001
Expected sample size (total):       218
Accrual rates:
     [,1] [,2] [,3]
[1,]    2    8    6
[2,]    4    3   10
Control event rates (H1):
       [,1]   [,2]   [,3]
[1,] 0.2310 0.1155 0.0770
[2,] 0.1733 0.0866 0.0578
[3,] 0.1386 0.0693 0.0462
Censoring rates:
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    0    0    0
[3,]    0    0    0
Power:                 100*(1-beta)=69.7639%
Type I error (1-sided):   100*alpha=2.5%
Randomization (Exp/Control):  ratio= 1 2 3

2.6.3 Checks

We will check same sample size and power calculations using all different methods we have using a common example. We do not print the output of each check due to volume, but allow the reader the option of running the code to make comparisons. Note that nSurvival() which was validated versus external results in a previous release of the gsDesign package and that LFPWE() was previously checked versus nSurvival(). Thus, we begin by deriving a sample size using LFPWE(). The input parameters are specifically kept separate at the beginning of the code so that the user can easily check using whatever parameters interest them. Note that nSurvival() will only work with a single stratum and failure rate period.

lambdaC <- log(2) / 20
hr <- .5
hr0 <- 1
S <- NULL
T <- 30
R <- 20
minfup <- 10
gamma <- 8
alpha <- .025
sided <- 1
alpha <- alpha
beta <- .1
etaC <- 0
etaE <- 0
ratio <- 1

gsDesign::nSurvival(
  lambda1 = lambdaC, lambda2 = lambdaC * hr, Ts = T,
  Tr = R, eta = etaC, alpha = alpha, sided = sided, beta = beta
)

Fixed design, two-arm trial with time-to-event
outcome (Lachin and Foulkes, 1986).
Study duration (fixed):          Ts=30
Accrual duration (fixed):        Tr=20
Uniform accrual:              entry="unif"
Control median:      log(2)/lambda1=20
Experimental median: log(2)/lambda2=40
Censoring only at study end (eta=0)
Control failure rate:       lambda1=0.035
Experimental failure rate:  lambda2=0.017
Censoring rate:                 eta=0
Power:                 100*(1-beta)=90%
Type I error (1-sided):   100*alpha=2.5%
Equal randomization:          ratio=1
Sample size based on hazard ratio=0.5 (type="rr")
Sample size (computed):           n=228
Events required (computed): nEvents=89

y <- LFPWE(
  gamma = gamma, hr = hr, lambdaC = lambdaC, R = R, T = T,
  minfup = minfup, alpha = alpha, sided = sided, beta = beta
)

Now we move on to derive the above result using KTZ() which computes power. With simple = TRUE (the default) and beta = NULL, the routine will should return the power to match the desired 90% input above. The first call derives R based on the input duration of accrual input in x. The second call sets the minimum follow-up duration to minfup to the input value of x.

KTZ(
  x = y$T - y$minfup, gamma = y$gamma, hr = y$hr,
  lambdaC = y$lambdaC, R = Inf, minfup = y$minfup
)
KTZ(
  x = y$minfup, gamma = y$gamma, hr = y$hr, lambdaC = y$lambdaC,
  minfup = NULL, R = y$T - y$minfup
)

Now we move on to deriving the sample size using KT(). The first call derives the duration of enrollment required and the second call the minimum follow-up required to obtain the desired power.

if (length(R) == 1) {
  temR <- Inf
} else {
  temR <- R
  temR[length(R)] <- Inf
}
KT(gamma = y, hr = hr, lambdaC = lambdaC, R = temR, minfup = minfup)
KT(gamma = y$gamma, hr = hr, lambdaC = lambdaC, R = R)

2.6.4 Examples

Using the output design as input, we compute the power both for hazard ratios of .6 and .75:

x <- KTZ(
  gamma = y$gamma, hr = .6, lambdaC = y$lambdaC,
  minfup = y$minfup, R = y$T - y$minfup, simple = F
)
x$power

[1] 0.69822

x <- KTZ(
  gamma = y$gamma, hr = .75, lambdaC = y$lambdaC,
  minfup = y$minfup, R = y$T - y$minfup, simple = F
)
x$power

[1] 0.3063416

2.7 nSurv: the general fixed design sample size and power function

The code for the general sample size and power routine nSurv() is short as it simply calls routines derived above. Following the code, we document the input and output for the routine.

2.7.1 Source code# nSurv roxy [sinew] ----
#' Advanced time-to-event sample size calculation
#'
#' \code{nSurv()} is used to calculate the sample size for a clinical trial
#' with a time-to-event endpoint and an assumption of proportional hazards.
#' This set of routines is new with version 2.7 and will continue to be
#' modified and refined to improve input error checking and output format with
#' subsequent versions. It allows both the Lachin and Foulkes (1986) method
#' (fixed trial duration) as well as the Kim and Tsiatis(1990) method (fixed
#' enrollment rates and either fixed enrollment duration or fixed minimum
#' follow-up). Piecewise exponential survival is supported as well as piecewise
#' constant enrollment and dropout rates. The methods are for a 2-arm trial
#' with treatment groups referred to as experimental and control. A stratified
#' population is allowed as in Lachin and Foulkes (1986); this method has been
#' extended to derive non-inferiority as well as superiority trials.
#' Stratification also allows power calculation for meta-analyses.
#' \code{gsSurv()} combines \code{nSurv()} with \code{gsDesign()} to derive a
#' group sequential design for a study with a time-to-event endpoint.
#'
#' \code{print()}, \code{xtable()} and \code{summary()} methods are provided to
#' operate on the returned value from \code{gsSurv()}, an object of class
#' \code{gsSurv}. \code{print()} is also extended to \code{nSurv} objects. The
#' functions \code{\link{gsBoundSummary}} (data frame for tabular output),
#' \code{\link{xprint}} (application of \code{xtable} for tabular output) and
#' \code{summary.gsSurv} (textual summary of \code{gsDesign} or \code{gsSurv}
#' object) may be preferred summary functions; see example in vignettes. See
#' also \link{gsBoundSummary} for output
#' of tabular summaries of bounds for designs produced by \code{gsSurv()}.
#'
#' Both \code{nEventsIA} and \code{tEventsIA} require a group sequential design
#' for a time-to-event endpoint of class \code{gsSurv} as input.
#' \code{nEventsIA} calculates the expected number of events under the
#' alternate hypothesis at a given interim time. \code{tEventsIA} calculates
#' the time that the expected number of events under the alternate hypothesis
#' is a given proportion of the total events planned for the final analysis.
#'
#' \code{nSurv()} produces an object of class \code{nSurv} 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.
#'
#' When study duration (\code{T}) and follow-up duration (\code{minfup}) are
#' fixed, \code{nSurv} applies exactly the Lachin and Foulkes (1986) method of
#' computing sample size under the proportional hazards assumption when For
#' this computation, enrollment rates are altered proportionately to those
#' input in \code{gamma} to achieve the power of interest.
#'
#' Given the specified enrollment rate(s) input in \code{gamma}, \code{nSurv}
#' may also be used to derive enrollment duration required for a trial to have
#' defined power if \code{T} is input as \code{NULL}; in this case, both
#' \code{R} (enrollment duration for each specified enrollment rate) and
#' \code{T} (study duration) will be computed on output.
#'
#' Alternatively and also using the fixed enrollment rate(s) in \code{gamma},
#' if minimum follow-up \code{minfup} is specified as \code{NULL}, then the
#' enrollment duration(s) specified in \code{R} are considered fixed and
#' \code{minfup} and \code{T} are computed to derive the desired power. This
#' method will fail if the specified enrollment rates and durations either
#' over-powers the trial with no additional follow-up or underpowers the trial
#' with infinite follow-up. This method produces a corresponding error message
#' in such cases.
#'
#' The input to \code{gsSurv} is a combination of the input to \code{nSurv()}
#' and \code{gsDesign()}.
#'
#' \code{nEventsIA()} is provided to compute the expected number of events at a
#' given point in time given enrollment, event and censoring rates. The routine
#' is used with a root finding routine to approximate the approximate timing of
#' an interim analysis. It is also used to extend enrollment or follow-up of a
#' fixed design to obtain a sufficient number of events to power a group
#' sequential design.
#'
#' @aliases nSurv print.nSurv print.gsSurv xtable.gsSurv
#' @param x An object of class \code{nSurv} or \code{gsSurv}.
#' \code{print.nSurv()} is used for an object of class \code{nSurv} which will
#' generally be output from \code{nSurv()}. For \code{print.gsSurv()} is used
#' for an object of class \code{gsSurv} which will generally be output from
#' \code{gsSurv()}. \code{nEventsIA} and \code{tEventsIA} operate on both the
#' \code{nSurv} and \code{gsSurv} class.
#' @param digits Number of digits past the decimal place to print
#' (\code{print.gsSurv.}); also a pass through to generic \code{xtable()} from
#' \code{xtable.gsSurv()}.
#' @param lambdaC scalar, vector or matrix of event hazard rates for the
#' control group; rows represent time periods while columns represent strata; a
#' vector implies a single stratum.
#' @param hr hazard ratio (experimental/control) under the alternate hypothesis
#' (scalar).
#' @param hr0 hazard ratio (experimental/control) under the null hypothesis
#' (scalar).
#' @param eta scalar, vector or matrix of dropout hazard rates for the control
#' group; 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 etaE matrix dropout hazard rates for the experimental group specified
#' in like form as \code{eta}; if NULL, this is set equal to \code{eta}.
#' @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 when variable enrollment duration is
#' specified (input \code{T=NULL}), 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 study duration; if \code{T} is input as \code{NULL}, this will be
#' computed on output; see details.
#' @param minfup follow-up of last patient enrolled; if \code{minfup} is input
#' as \code{NULL}, this will be computed on output; see details.
#' @param ratio randomization ratio of experimental treatment divided by
#' control; normally a scalar, but may be a vector with length equal to number
#' of strata.
#' @param sided 1 for 1-sided testing, 2 for 2-sided testing.
#' @param alpha type I error rate. Default is 0.025 since 1-sided testing is
#' default.
#' @param beta type II error rate. Default is 0.10 (90\% power); NULL if power
#' is to be computed based on other input values.
#' @param tol for cases when \code{T} or \code{minfup} values are derived
#' through root finding (\code{T} or \code{minfup} input as \code{NULL}),
#' \code{tol} provides the level of error input to the \code{uniroot()}
#' root-finding function. The default is the same as for \code{\link{uniroot}}.
#' @param k Number of analyses planned, including interim and final.
#' @param test.type \code{1=}one-sided \cr \code{2=}two-sided symmetric \cr
#' \code{3=}two-sided, asymmetric, beta-spending with binding lower bound \cr
#' \code{4=}two-sided, asymmetric, beta-spending with non-binding lower bound
#' \cr \code{5=}two-sided, asymmetric, lower bound spending under the null
#' hypothesis with binding lower bound \cr \code{6=}two-sided, asymmetric,
#' lower bound spending under the null hypothesis with non-binding lower bound.
#' \cr See details, examples and manual.
#' @param astar Normally not specified. If \code{test.type=5} or \code{6},
#' \code{astar} specifies the total probability of crossing a lower bound at
#' all analyses combined.  This will be changed to \eqn{1 - }\code{alpha} when
#' default value of 0 is used.  Since this is the expected usage, normally
#' \code{astar} is not specified by the user.
#' @param timing Sets relative timing of interim analyses in \code{gsSurv}.
#' Default of 1 produces equally spaced analyses.  Otherwise, this is a vector
#' of length \code{k} or \code{k-1}.  The values should satisfy \code{0 <
#' timing[1] < timing[2] < ... < timing[k-1] < timing[k]=1}. For
#' \code{tEventsIA}, this is a scalar strictly between 0 and 1 that indicates
#' the targeted proportion of final planned events available at an interim
#' analysis.
#' @param sfu A spending function or a character string indicating a boundary
#' type (that is, \dQuote{WT} for Wang-Tsiatis bounds, \dQuote{OF} for
#' O'Brien-Fleming bounds and \dQuote{Pocock} for Pocock bounds).  For
#' one-sided and symmetric two-sided testing is used to completely specify
#' spending (\code{test.type=1, 2}), \code{sfu}.  The default value is
#' \code{sfHSD} which is a Hwang-Shih-DeCani spending function.  See details,
#' \link{Spending_Function_Overview}, manual and examples.
#' @param sfupar Real value, default is \eqn{-4} which is an
#' O'Brien-Fleming-like conservative bound when used with the default
#' Hwang-Shih-DeCani spending function. This is a real-vector for many spending
#' functions.  The parameter \code{sfupar} specifies any parameters needed for
#' the spending function specified by \code{sfu}; this will be ignored for
#' spending functions (\code{sfLDOF}, \code{sfLDPocock}) or bound types
#' (\dQuote{OF}, \dQuote{Pocock}) that do not require parameters.
#' @param sfl Specifies the spending function for lower boundary crossing
#' probabilities when asymmetric, two-sided testing is performed
#' (\code{test.type = 3}, \code{4}, \code{5}, or \code{6}).  Unlike the upper
#' bound, only spending functions are used to specify the lower bound.  The
#' default value is \code{sfHSD} which is a Hwang-Shih-DeCani spending
#' function.  The parameter \code{sfl} is ignored for one-sided testing
#' (\code{test.type=1}) or symmetric 2-sided testing (\code{test.type=2}).  See
#' details, spending functions, manual and examples.
#' @param sflpar Real value, default is \eqn{-2}, which, with the default
#' Hwang-Shih-DeCani spending function, specifies a less conservative spending
#' rate than the default for the upper bound.
#' @param r Integer value controlling grid for numerical integration as in
#' Jennison and Turnbull (2000); default is 18, range is 1 to 80. Larger values
#' provide larger number of grid points and greater accuracy.  Normally
#' \code{r} will not be changed by the user.
#' @param usTime Default is NULL in which case upper bound spending time is
#' determined by \code{timing}. Otherwise, this should be a vector of length
#' code{k} with the spending time at each analysis (see Details in help for \code{gsDesign}).
#' @param lsTime Default is NULL in which case lower bound spending time is
#' determined by \code{timing}. Otherwise, this should be a vector of length
#' \code{k} with the spending time at each analysis (see Details in help for \code{gsDesign}).
#' @param tIA Timing of an interim analysis; should be between 0 and
#' \code{y$T}.
#' @param target The targeted proportion of events at an interim analysis. This
#' is used for root-finding will be 0 for normal use.
#' @param simple See output specification for \code{nEventsIA()}.
#' @param footnote footnote for xtable output; may be useful for describing
#' some of the design parameters.
#' @param fnwid a text string controlling the width of footnote text at the
#' bottom of the xtable output.
#' @param timename character string with plural of time units (e.g., "months")
#' @param caption passed through to generic \code{xtable()}.
#' @param label passed through to generic \code{xtable()}.
#' @param align passed through to generic \code{xtable()}.
#' @param display passed through to generic \code{xtable()}.
#' @param auto passed through to generic \code{xtable()}.
#' @param ... other arguments that may be passed to generic functions
#' underlying the methods here.
#' @return \code{nSurv()} returns an object of type \code{nSurv} with the
#' following components: \item{alpha}{As input.} \item{sided}{As input.}
#' \item{beta}{Type II error; if missing, this is computed.} \item{power}{Power
#' corresponding to input \code{beta} or computed if output \code{beta} is
#' computed.} \item{lambdaC}{As input.} \item{etaC}{As input.} \item{etaE}{As
#' input.} \item{gamma}{As input unless none of the following are \code{NULL}:
#' \code{T}, \code{minfup}, \code{beta}; otherwise, this is a constant times
#' the input value required to power the trial given the other input
#' variables.} \item{ratio}{As input.} \item{R}{As input unless \code{T} was
#' \code{NULL} on input.} \item{S}{As input.} \item{T}{As input.}
#' \item{minfup}{As input.} \item{hr}{As input.} \item{hr0}{As input.}
#' \item{n}{Total expected sample size corresponding to output accrual rates
#' and durations.} \item{d}{Total expected number of events under the alternate
#' hypothesis.} \item{tol}{As input, except when not used in computations in
#' which case this is returned as \code{NULL}.  This and the remaining output
#' below are not printed by the \code{print()} extension for the \code{nSurv}
#' class.} \item{eDC}{A vector of expected number of events by stratum in the
#' control group under the alternate hypothesis.} \item{eDE}{A vector of
#' expected number of events by stratum in the experimental group under the
#' alternate hypothesis.} \item{eDC0}{A vector of expected number of events by
#' stratum in the control group under the null hypothesis.} \item{eDE0}{A
#' vector of expected number of events by stratum in the experimental group
#' under the null hypothesis.} \item{eNC}{A vector of the expected accrual in
#' each stratum in the control group.} \item{eNE}{A vector of the expected
#' accrual in each stratum in the experimental group.} \item{variable}{A text
#' string equal to "Accrual rate" if a design was derived by varying the
#' accrual rate, "Accrual duration" if a design was derived by varying the
#' accrual duration, "Follow-up duration" if a design was derived by varying
#' follow-up duration, or "Power" if accrual rates and duration as well as
#' follow-up duration was specified and \code{beta=NULL} was input.}
#'
#' \code{gsSurv()} returns much of the above plus an object of class
#' \code{gsDesign} in a variable named \code{gs}; see \code{\link{gsDesign}}
#' for general documentation on what is returned in \code{gs}.  The value of
#' \code{gs$n.I} represents the number of endpoints required at each analysis
#' to adequately power the trial. Other items returned by \code{gsSurv()} are:
#' \item{gs}{A group sequential design (\code{gsDesign}) object.}
#' \item{lambdaC}{As input.} \item{etaC}{As input.} \item{etaE}{As input.}
#' \item{gamma}{As input unless none of the following are \code{NULL}:
#' \code{T}, \code{minfup}, \code{beta}; otherwise, this is a constant times
#' the input value required to power the trial given the other input
#' variables.} \item{ratio}{As input.} \item{R}{As input unless \code{T} was
#' \code{NULL} on input.} \item{S}{As input.} \item{T}{As input.}
#' \item{minfup}{As input.} \item{hr}{As input.} \item{hr0}{As input.}
#' \item{eNC}{Total expected sample size corresponding to output accrual rates
#' and durations.} \item{eNE}{Total expected sample size corresponding to
#' output accrual rates and durations.} \item{eDC}{Total expected number of
#' events under the alternate hypothesis.} \item{eDE}{Total expected number of
#' events under the alternate hypothesis.} \item{tol}{As input, except when not
#' used in computations in which case this is returned as \code{NULL}.  This
#' and the remaining output below are not printed by the \code{print()}
#' extension for the \code{nSurv} class.} \item{eDC}{A vector of expected
#' number of events by stratum in the control group under the alternate
#' hypothesis.} \item{eDE}{A vector of expected number of events by stratum in
#' the experimental group under the alternate hypothesis.} \item{eDC0}{A vector
#' of expected number of events by stratum in the control group under the null
#' hypothesis.} \item{eDE0}{A vector of expected number of events by stratum in
#' the experimental group under the null hypothesis.} \item{eNC}{A vector of
#' the expected accrual in each stratum in the control group.} \item{eNE}{A
#' vector of the expected accrual in each stratum in the experimental group.}
#' \item{variable}{A text string equal to "Accrual rate" if a design was
#' derived by varying the accrual rate, "Accrual duration" if a design was
#' derived by varying the accrual duration, "Follow-up duration" if a design
#' was derived by varying follow-up duration, or "Power" if accrual rates and
#' duration as well as follow-up duration was specified and \code{beta=NULL}
#' was input.}
#'
#' \code{nEventsIA()} returns the expected proportion of the final planned
#' events observed at the input analysis time minus \code{target} when
#' \code{simple=TRUE}. When \code{simple=FALSE}, \code{nEventsIA} returns a
#' list with following components: \item{T}{The input value \code{tIA}.}
#' \item{eDC}{The expected number of events in the control group at time the
#' output time \code{T}.} \item{eDE}{The expected number of events in the
#' experimental group at the output time \code{T}.} \item{eNC}{The expected
#' enrollment in the control group at the output time \code{T}.} \item{eNE}{The
#' expected enrollment in the experimental group at the output time \code{T}.}
#'
#' \code{tEventsIA()} returns
#' @author Keaven Anderson \email{keaven_anderson@@merck.com}
#' @seealso \code{\link{gsBoundSummary}}, \code{\link{xprint}},
#' \link{gsDesign package overview}, \link{plot.gsDesign},
#' \code{\link{gsDesign}}, \code{\link{gsHR}}, \code{\link{nSurvival}}
#' @references Kim KM and Tsiatis AA (1990), Study duration for clinical trials
#' with survival response and early stopping rule. \emph{Biometrics}, 46, 81-92
#'
#' 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.
#'
#' Schoenfeld D (1981), The Asymptotic Properties of Nonparametric Tests for
#' Comparing Survival Distributions. \emph{Biometrika}, 68, 316-319.
#' @keywords design
#' @examples
#'
#' # vary accrual rate to obtain power
#' nSurv(lambdaC = log(2) / 6, hr = .5, eta = log(2) / 40, gamma = 1, T = 36, minfup = 12)
#'
#' # vary accrual duration to obtain power
#' nSurv(lambdaC = log(2) / 6, hr = .5, eta = log(2) / 40, gamma = 6, minfup = 12)
#'
#' # vary follow-up duration to obtain power
#' nSurv(lambdaC = log(2) / 6, hr = .5, eta = log(2) / 40, gamma = 6, R = 25)
#'
#' # piecewise constant enrollment rates (vary accrual duration)
#' nSurv(
#'   lambdaC = log(2) / 6, hr = .5, eta = log(2) / 40, gamma = c(1, 3, 6),
#'   R = c(3, 6, 9), minfup = 12
#' )
#'
#' # stratified population (vary accrual duration)
#' nSurv(
#'   lambdaC = matrix(log(2) / c(6, 12), ncol = 2), hr = .5, eta = log(2) / 40,
#'   gamma = matrix(c(2, 4), ncol = 2), minfup = 12
#' )
#'
#' # piecewise exponential failure rates (vary accrual duration)
#' nSurv(lambdaC = log(2) / c(6, 12), hr = .5, eta = log(2) / 40, S = 3, gamma = 6, minfup = 12)
#'
#' # combine it all: 2 strata, 2 failure rate periods
#' nSurv(
#'   lambdaC = matrix(log(2) / c(6, 12, 18, 24), ncol = 2), hr = .5,
#'   eta = matrix(log(2) / c(40, 50, 45, 55), ncol = 2), S = 3,
#'   gamma = matrix(c(3, 6, 5, 7), ncol = 2), R = c(5, 10), minfup = 12
#' )
#'
#' # example where only 1 month of follow-up is desired
#' # set failure rate to 0 after 1 month using lambdaC and S
#' nSurv(lambdaC = c(.4, 0), hr = 2 / 3, S = 1, minfup = 1)
#'
#' # group sequential design (vary accrual rate to obtain power)
#' x <- gsSurv(
#'   k = 4, sfl = sfPower, sflpar = .5, lambdaC = log(2) / 6, hr = .5,
#'   eta = log(2) / 40, gamma = 1, T = 36, minfup = 12
#' )
#' x
#' print(xtable::xtable(x,
#'   footnote = "This is a footnote; note that it can be wide.",
#'   caption = "Caption example."
#' ))
#' # find expected number of events at time 12 in the above trial
#' nEventsIA(x = x, tIA = 10)
#'
#' # find time at which 1/4 of events are expected
#' tEventsIA(x = x, timing = .25)
#' @export
# nSurv function [sinew] ----
nSurv <- function(lambdaC = log(2) / 6, hr = .6, hr0 = 1, eta = 0, etaE = NULL,
                  gamma = 1, R = 12, S = NULL, T = NULL, minfup = NULL, ratio = 1,
                  alpha = 0.025, beta = 0.10, sided = 1, tol = .Machine$double.eps^0.25) {
  if (is.null(etaE)) etaE <- eta
  # set up rates as matrices with row and column names
  # default is 1 stratum if lambdaC not input as matrix
  if (is.vector(lambdaC)) lambdaC <- matrix(lambdaC)
  ldim <- dim(lambdaC)
  nstrata <- ldim[2]
  nlambda <- ldim[1]
  rownames(lambdaC) <- paste("Period", 1:nlambda)
  colnames(lambdaC) <- paste("Stratum", 1:nstrata)
  etaC <- matrix(eta, nrow = nlambda, ncol = nstrata)
  etaE <- matrix(etaE, nrow = nlambda, ncol = nstrata)
  if (!is.matrix(gamma)) gamma <- matrix(gamma)
  gdim <- dim(gamma)
  eCdim <- dim(etaC)
  eEdim <- dim(etaE)

  if (is.null(minfup) || is.null(T)) {
    xx <- KT(
      lambdaC = lambdaC, hr = hr, hr0 = hr0, etaC = etaC, etaE = etaE,
      gamma = gamma, R = R, S = S, minfup = minfup, ratio = ratio,
      alpha = alpha, sided = sided, beta = beta, tol = tol
    )
  } else if (is.null(beta)) {
    xx <- KTZ(
      lambdaC = lambdaC, hr = hr, hr0 = hr0, etaC = etaC, etaE = etaE,
      gamma = gamma, R = R, S = S, minfup = minfup, ratio = ratio,
      alpha = alpha, sided = sided, beta = beta, simple = F
    )
  } else {
    xx <- LFPWE(
      lambdaC = lambdaC, hr = hr, hr0 = hr0, etaC = etaC, etaE = etaE,
      gamma = gamma, R = R, S = S, T = T, minfup = minfup, ratio = ratio,
      alpha = alpha, sided = sided, beta = beta
    )
  }

  nameR <- nameperiod(cumsum(xx$R))
  stratnames <- paste("Stratum", 1:ncol(xx$lambdaC))
  if (is.null(xx$S)) {
    nameS <- "0-Inf"
  } else {
    nameS <- nameperiod(cumsum(c(xx$S, Inf)))
  }
  rownames(xx$lambdaC) <- nameS
  colnames(xx$lambdaC) <- stratnames
  rownames(xx$etaC) <- nameS
  colnames(xx$etaC) <- stratnames
  rownames(xx$etaE) <- nameS
  colnames(xx$etaE) <- stratnames
  rownames(xx$gamma) <- nameR
  colnames(xx$gamma) <- stratnames
  return(xx)
}

2.7.2 Input variables

The input variables for nSurv() are:

lambdaC: scalar, vector or matrix of event hazard rates for the control group; rows represent time periods while columns represent strata; a vector implies a single stratum.
hr: hazard ratio (experimental/control) under the alternate hypothesis (scalar).
hr0: hazard ratio (experimental/control) under the null hypothesis (scalar).
eta: scalar, vector or matrix of dropout hazard rates for the control group; rows represent time periods while columns represent strata; if entered as a scalar, rate is constant accross strata and time periods; if entered as a vector, rates are constant accross strata.
etaE: matrix dropout hazard rates for the experimental group specified in like form as eta; if NULL, this is set equal to eta.
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 accross strata and time periods; if entered as a vector, rates are constant accross strata.
R: a scalar or vector of durations of time periods for recruitment rates specified in rows of gamma. Length is the same as number of rows in gamma. Note that when variable enrollment duration is specified, final enrollment period is extended as long as needed.
S: a scalar or vector of durations of piecewise constant event rates specified in rows of lambda, eta and etaE; this is NULL if there is a single event rate per stratum (exponential failure) or of length 1 less than the number of rows in lambda otherwise.
T: study duration.
minfup: follow-up of last patient enrolled.
ratio: randomization ratio of experimental treatment divided by control; normally a scalar, but can be a vector with length equal to number of strata.
sided: 1 for 1-sided testing, 2 for 2-sided testing.
alpha: type I error rate. Default is 0.025 since 1-sided testing is default.
beta: type II error rate. Default is 0.10 (\(90\%\) power); NULL if power is to be computed based on other input values.
tol: for cases when values are derived through route finding, the level of error input to the uniroot() root-finding function.

2.7.3 Output structure

The output of nSurv() is of class nSurv, which is a list of the following items:

alpha: As input.
sided: As input.
beta: Type II error; if missing, this is computed.
power: Power corresponding to input beta or computed if output beta is computed.
lambdaC: As input.
etaC: As input.
etaE: As input.
gamma: As input unless none of the following are NULL: T, minfup, beta; otherwise, this is a constant times the input value required to power the trial given the other input variables.
ratio: As input.
R: As input unless T was NULL on input.
S: As input.
T: As input.
minfup: As input.
hr: As input.
hr0: As input.
n: Total expected sample size corresponding to output accrual rates and durations.
d: Total expected number of events under the alternate hypothesis.
tol: As input. This and the remaining output below are not printed by the print() extension for the nSurv class.
eDC: A vector of expected number of events by stratum in the control group under the alternate hypothesis.
eDE: A vector of expected number of events by stratum in the experimental group under the alternate hypothesis.
eDC0: A vector of expected number of events by stratum in the control group under the null hypothesis.
eDE0: A vector of expected number of events by stratum in the experimental group under the null hypothesis.
eNC: A vector of the expected accrual in each stratum in the control group.
eNE: A vector of the expected accrual in each stratum in the experimental group.

2.8 Documentation and primary examples for fixed design sample size and power

We present four examples for the advanced time-to-event sample size routine nSurv() corresponding to the four types of calculations peformed. All assume an input hazard ratio, dropout rates, Type I error rate. The first three assume an input Type II error rate, while the fourth computes power. The types of calculations are:

Based on a fixed enrollment duration and follow-up period, derive the enrollment rates required to obtain adequate power (Lachin and Foulkes 1986).
Based on fixed enrollment rates and a fixed minimum follow-up period, derive accrual duration required to obtain adequate power (Kim and Tsiatis 1990).
Based on fixed enrollment rates and accrual duration, derive study duration required to obtain adequate power (Kim and Tsiatis 1990).
Based on fixed accrual duration, accrual rates and follow-up, derive study power.

2.8.1 Derive accrual rates

We begin by specifying values for all three of 1) control event rates (lambdaC), 2) accrual duration (R), and 3) minimum follow-up (minfup). Relative accrual rates are specified for the time periods specified in R in gamma; these are modified on output to give absolute accrual rates required to achieve the desired power.

x <- nSurv(
  lambdaC = lambdaC, hr = hr, hr0 = hr0, eta = etaC, etaE = etaE,
  gamma = gamma, R = R, S = S, T = T, minfup = minfup, ratio = ratio,
  alpha = alpha, beta = beta
)

2.8.2 Derive accrual and study duration

Now we set the input study duration (T) to NULL so that it can be adjusted to achieve desired power; enrollment duration is altered to total study duration minus the input minimum follow-up.

x <- nSurv(
  lambdaC = lambdaC, hr = hr, hr0 = hr0, eta = etaC, etaE = etaE,
  gamma = gamma, R = R, S = S, T = NULL, minfup = minfup, ratio = ratio,
  alpha = alpha, beta = beta
)
x

Fixed design, two-arm trial with time-to-event
outcome (Lachin and Foulkes, 1986).
Solving for:  Accrual duration 
Hazard ratio                  H1/H0=0.5/1
Study duration:                   T=35.836
Accrual duration:                   25.836
Min. end-of-study follow-up: minfup=10
Expected events (total, H1):        88.3566
Expected sample size (total):       206.6883
Accrual rates:
        Stratum 1
0-25.84         8
Control event rates (H1):
      Stratum 1
0-Inf    0.0347
Censoring rates:
      Stratum 1
0-Inf         0
Power:                 100*(1-beta)=90%
Type I error (1-sided):   100*alpha=2.5%
Equal randomization:          ratio=1

KTZ(
  x = x$T - sum(x$R), lambdaC = lambdaC, hr = hr, hr0 = hr0, eta = etaC, etaE = etaE,
  gamma = gamma, R = x$R, S = x$S, minfup = NULL, ratio = ratio,
  alpha = alpha, beta = .1, simple = TRUE
)

[1] 6.233063e-08

2.8.3 Derive follow-up duration

Next we set the input follow-up (minfup) to NULL so that it can be derived to achieve the desired power. Note that this routine will fail if accrual is too little or too much in the specified time period. Thus, it is often easier to allow study and accrual duration to vary as just shown.

x <- nSurv(
  lambdaC = lambdaC, hr = hr, hr0 = hr0, eta = etaC, etaE = etaE,
  gamma = gamma, R = R, S = S, T = T, minfup = NULL, ratio = ratio,
  alpha = alpha, beta = beta
)

2.8.4 Derive power

Finally, we compute the power for a design by setting the input beta to NULL.

nSurv(
  lambdaC = lambdaC, hr = hr, hr0 = hr0, eta = etaC, etaE = etaE,
  gamma = gamma, R = R, S = S, T = T, minfup = minfup, ratio = ratio,
  alpha = alpha, beta = NULL
)

Fixed design, two-arm trial with time-to-event
outcome (Lachin and Foulkes, 1986).
Solving for:  Power 
Hazard ratio                  H1/H0=0.5/1
Study duration:                   T=30
Accrual duration:                   20
Min. end-of-study follow-up: minfup=10
Expected events (total, H1):        62.3423
Expected sample size (total):       160
Accrual rates:
     Stratum 1
0-20         8
Control event rates (H1):
      Stratum 1
0-Inf    0.0347
Censoring rates:
      Stratum 1
0-Inf         0
Power:                 100*(1-beta)=77.9917%
Type I error (1-sided):   100*alpha=2.5%
Equal randomization:          ratio=1

2.8.5 Examples with multiple strata, enrollment rates and failure rates

Fixed accrual duration and minimum follow-up duration — solve for accrual rate.

x <- nSurv(
  lambdaC = log(2) / matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3),
  gamma = matrix(c(2, 4, 8, 3, 6, 10), nrow = 2),
  hr = .6, S = c(3, 6), R = c(3, 3), minfup = 6, T = 30
)
x

Fixed design, two-arm trial with time-to-event
outcome (Lachin and Foulkes, 1986).
Solving for:  Accrual rate 
Hazard ratio                  H1/H0=0.6/1
Study duration:                   T=30
Accrual duration:                   24
Min. end-of-study follow-up: minfup=6
Expected events (total, H1):        160.4241
Expected sample size (total):       255.7091
Accrual rates:
     Stratum 1 Stratum 2 Stratum 3
0-3     1.2628    5.0510    3.7883
3-24    2.5255    1.8941    6.3138
Control event rates (H1):
      Stratum 1 Stratum 2 Stratum 3
0-3      0.2310    0.1155    0.0770
3-9      0.1733    0.0866    0.0578
9-Inf    0.1386    0.0693    0.0462
Censoring rates:
      Stratum 1 Stratum 2 Stratum 3
0-3           0         0         0
3-9           0         0         0
9-Inf         0         0         0
Power:                 100*(1-beta)=90%
Type I error (1-sided):   100*alpha=2.5%
Equal randomization:          ratio=1

Fixed accrual rate and minimum follow-up duration — solve for accrual duration.

x <- nSurv(
  lambdaC = log(2) / matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3),
  gamma = matrix(c(2, 4, 8, 3, 6, 10), nrow = 2),
  hr = .6, S = c(3, 6), R = c(3, 3), minfup = 6, T = NULL
)
x

Fixed design, two-arm trial with time-to-event
outcome (Lachin and Foulkes, 1986).
Solving for:  Accrual duration 
Hazard ratio                  H1/H0=0.6/1
Study duration:                   T=22.647
Accrual duration:                   16.647
Min. end-of-study follow-up: minfup=6
Expected events (total, H1):        160.7145
Expected sample size (total):       279.9995
Accrual rates:
        Stratum 1 Stratum 2 Stratum 3
0-3             2         8         6
3-16.65         4         3        10
Control event rates (H1):
      Stratum 1 Stratum 2 Stratum 3
0-3      0.2310    0.1155    0.0770
3-9      0.1733    0.0866    0.0578
9-Inf    0.1386    0.0693    0.0462
Censoring rates:
      Stratum 1 Stratum 2 Stratum 3
0-3           0         0         0
3-9           0         0         0
9-Inf         0         0         0
Power:                 100*(1-beta)=90%
Type I error (1-sided):   100*alpha=2.5%
Equal randomization:          ratio=1

Fixed accrual rate and accrual duration — solve for minimum follow-up. Note that this option is generally not recommended as the accrual rate or accrual duration are the more common variables to solve for.

x <- nSurv(
  lambdaC = log(2) / matrix(c(3, 4, 5, 6, 8, 10, 9, 12, 15), nrow = 3),
  gamma = matrix(c(2, 4, 8, 3, 6, 10), nrow = 2),
  hr = .6, S = c(3, 6), R = c(3, 15), minfup = NULL, T = 15
)
x

Fixed design, two-arm trial with time-to-event
outcome (Lachin and Foulkes, 1986).
Solving for:  Follow-up duration 
Hazard ratio                  H1/H0=0.6/1
Study duration:                   T=21.5363
Accrual duration:                   18
Min. end-of-study follow-up: minfup=3.5363
Expected events (total, H1):        160.8979
Expected sample size (total):       303
Accrual rates:
     Stratum 1 Stratum 2 Stratum 3
0-3          2         8         6
3-18         4         3        10
Control event rates (H1):
      Stratum 1 Stratum 2 Stratum 3
0-3      0.2310    0.1155    0.0770
3-9      0.1733    0.0866    0.0578
9-Inf    0.1386    0.0693    0.0462
Censoring rates:
      Stratum 1 Stratum 2 Stratum 3
0-3           0         0         0
3-9           0         0         0
9-Inf         0         0         0
Power:                 100*(1-beta)=90%
Type I error (1-sided):   100*alpha=2.5%
Equal randomization:          ratio=1