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 consider 4 approaches:

Lachin and Foulkes (1986): default recommendation.
Schoenfeld (1981): commonly used method.
Freedman (1982): may be useful for cases with large hazard reduction.
Bernstein and Lagakos (1978): while this is not recommended, we include it for completeness.

2.3.1 Lachin and Foulkes (1986)

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.3.2 Other methods

All fixed-design sample size routines now accept a method argument with default "LachinFoulkes" and alternatives "Schoenfeld", "Freedman", and "BernsteinLagakos". The same notation as above applies, and the implementations extend naturally to multiple strata, unequal randomization (ratio), and non-inferiority designs (hr0). Note that "Schoenfeld" and "Freedman" methods only support superiority testing (hr0 = 1). Additionally, "Freedman" does not support stratified populations.

The methods of Schoenfeld (1981), Freedman (1982), and Bernstein and Lagakos (1978) are all based on the number of events. Thus, to calculate sample size, we first calculate the number of events required and then find the sample size that yields this number of events. Schoenfeld (1981) and Freedman (1982) methods are generally stated in terms of the total number of events \(D\) required to achieve power \(1-\beta\) with type I error \(\alpha\) (1-sided). Both of these methods use only the null hypothesis variance in computing sample size. The Bernstein and Lagakos (1978) method is based on the same asymptotic approximation as Lachin and Foulkes (1986), but uses a different variance approximation for the log hazard ratio under the null hypothesis.

2.3.2.1 Schoenfeld (1981)

Schoenfeld (1981) derived the asymptotic distribution of the log-rank statistic under the assumption of proportional hazards and local alternatives (i.e., \(\ln(HR)\) is small). The total number of events \(D\) required is \[\begin{equation} D = \frac{(Z_\alpha + Z_\beta)^2}{\xi(1-\xi)[\ln(HR)]^2}. \label{eq:schoenfeld} \end{equation}\] This formula effectively assumes that the variance of the log hazard ratio is the same under the null and alternative hypotheses, specifically equal to \(1/[\xi(1-\xi)D]\). This corresponds to the variance under the null hypothesis where the event rate is the same in both groups. This method does not explicitly account for the different event rates in the control and experimental groups under the alternative hypothesis, nor does it account for the potential difference in censoring distributions or enrollment patterns beyond their effect on the total number of events.

Restrictions: The "Schoenfeld" method only supports superiority testing (hr0 = 1).

2.3.2.2 Freedman (1982)

Freedman (1982) provided a similar formula, but phrased in terms of the survival probabilities or hazard rates. Freedman also denoted the randomization ratio \(\phi\) and hazard ratio \(\theta\) as control:experimental which is one divided by our usage here (experimental:control). The Freedman result with the control:experimental convention is: \[\begin{equation} D = (Z_\alpha + Z_\beta)^2\frac{(1 + \theta\phi)^2}{\phi(1-\theta)^2}. \label{eq:freedman} \end{equation}\] Replacing \(\phi\) with \(\xi=1/(1+\phi)\) and \(\theta\) with \(1/HR\), this can be rewritten as \[\begin{equation} D = (Z_\alpha + Z_\beta)^2 \frac{(HR + (1-\xi)/\xi)^2} {\frac{1-\xi}{\xi}(HR-1)^{2}}. \label{eq:freedman2} \end{equation}\]

This is very similar to the Schoenfeld formula when there is equal randomization (\(\phi =1, \xi = 1/2\)); for \(HR\) close to 1, the Freedman contrast approximates \(\ln(HR)\) since \(2\frac{HR-1}{HR+1} \approx \ln(HR)\). Like the Schoenfeld method, it assumes equal variance under the null and alternative hypotheses.

Restrictions: The "Freedman" method only supports superiority testing (hr0 = 1) and does not support stratified populations.

2.3.2.3 Bernstein and Lagakos (1978)

Bernstein and Lagakos (1978) proposed a method that, like Lachin and Foulkes (1986), accounts for the difference in variance under the null and alternative hypotheses. However, they use a specific formulation for the variance under the alternative hypothesis. Using the notation from the previous section, the variance of the log hazard ratio estimator under the alternative hypothesis is estimated as: \[\begin{equation} V_1 = \frac{1}{E\left\{\delta|\lambda_{11},\eta_1\right\}} + \frac{1}{E\left\{\delta|\lambda_{10},\eta_0\right\}}. \end{equation}\] They then assume this same variance applies under the null hypothesis (effectively using the alternative hypothesis event counts to estimate the null variance). The sample size formula then becomes: \[\begin{equation} \sqrt{N}\ln(HR) = (Z_\alpha + Z_\beta) \sqrt{V_1 N}. \end{equation}\] This can be rewritten to solve for \(N\) (or effectively \(D\) if we consider the event counts): \[\begin{equation} D = \frac{(Z_\alpha + Z_\beta)^2 (1/p_C + 1/p_E)}{[\ln(HR)]^2} \end{equation}\] where \(p_C\) and \(p_E\) are the proportions of events in the control and experimental groups, respectively. Compared to Lachin and Foulkes (1986), this method tends to underestimate the variance under the null hypothesis (since event rates are typically lower/variance higher under H0 for superiority trials), leading to slightly smaller sample sizes. For stratified designs, Bernstein and Lagakos (1978) sum the variances across strata, similar to the method described above.

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,
                  method = c("LachinFoulkes", "Schoenfeld", "Freedman", "BernsteinLagakos")) {
  method <- match.arg(method)
  # set up parameters
  zalpha <- -qnorm(alpha / sided)
  zbeta <- if (is.null(beta)) NULL else -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}
  # ratio is validated in nSurv/gsSurv/gsSurvCalendar to be a single positive scalar
  Qe <- ratio / (1 + ratio) # Proportion of subjects in the experimental group
  Qc <- 1 - Qe
  # allocation of accrual by arm
  gammaC <- gamma * Qc
  gammaE <- gamma * Qe

  # For all methods, we need expected events under H1
  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
  )
  # Need expected events under control hazards for BernsteinLagakos and LF null variance
  eDC0_list <- eEvents(
    lambda = lambdaC, eta = etaC, gamma = gammaC,
    R = R, S = S, T = T, minfup = minfup
  )
  eDE0_list <- eEvents(
    lambda = lambdaC, eta = etaE, gamma = gammaE,
    R = R, S = S, T = T, minfup = minfup
  )
  eDC0 <- eDC0_list$d
  eDE0 <- eDE0_list$d
  # Compute total n for each  treatment group
  nC <- sum(eDC$n)
  nE <- sum(eDE$n)
  # For a subject enrolled in  a treatment group, what is the probability
  # for each stratum that they are in that stratum and have an event?
  pC1 <- eDC$d / nC
  pE1 <- eDE$d / nE
  pC0 <- eDC0 / sum(eDC0_list$n)
  pE0 <- eDE0 / sum(eDE0_list$n)
  # Inverse-variance weighting across strata
  v1_strata <- (1 / pC1) + (1 / pE1)
  v1_strat <- 1 / sum(1 / v1_strata)
  v0_strata <- (1 / pC0) + (1 / pE0)
  v0_strat <- 1 / sum(1 / v0_strata)
  # Total event proportion under H1 (for scaling events to sample size)
  total_d1 <- sum(eDC$d + eDE$d) / sum(eDC$n + eDE$n)
  # Expected enrollment per unit accrual rate (needed for event-based methods)
  mx <- sum(eDC$n + eDE$n)
  # Initialize power_val for methods that compute power directly
  power_val <- NULL
  if (identical(method, "LachinFoulkes")) {
    # compute H0 failure rates as average of control, experimental
    # This is used for LachinFoulkes only
    lambdaC0 <- (1 + hr * ratio) / (1 + hr0 * ratio) * lambdaC
    # 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
    # Compute variances
    var0_lf <- 1 / sum((1 / eDC0 + 1 / eDE0)^(-1))
    var1_lf <- 1 / sum((1 / eDC$d + 1 / eDE$d)^(-1))
    delta_log <- abs(log(hr / hr0))
    if (is.null(beta)) {
      # Power calculation matching vignette formula:
      # power = pnorm(-(qnorm(1-alpha)*sqrt(var0) - |delta|)/sqrt(var1))
      # Note: variances are already computed for the given gamma (sample size)
      power_val <- pnorm(-(zalpha * sqrt(var0_lf) - delta_log) / sqrt(var1_lf))
      n <- 1  # Multiplier is 1 for power calculation
    } else {
      # Sample size calculation
      n <- ((zalpha * sqrt(var0_lf) + zbeta * sqrt(var1_lf)) / delta_log)^2
    }
  } else if (identical(method, "Schoenfeld")) {
    # Schoenfeld formula: D = (zalpha + zbeta)^2 / (Qe*(1-Qe) * (log(HR))^2)
    # Check for non-inferiority/super-superiority (hr0 must equal 1)
    if (hr0 != 1) {
      stop("Schoenfeld method only supports superiority testing (hr0 = 1)")
    }
    # For stratified: use inverse-variance weighted variance
    # Check for single vs multiple strata
    if (!is.matrix(lambdaC) || ncol(lambdaC) == 1) {
      # Single stratum: use standard formula
      var_schoen <- 1 / (Qe * (1 - Qe))
      if (is.null(beta)) {
        # Power calculation: solve for zbeta from sample size
        # D = (zalpha + zbeta)^2 * var / delta^2
        # sqrt(D)*delta = (zalpha + zbeta)*sqrt(var)
        # zbeta = sqrt(D)*delta/sqrt(var) - zalpha
        n_current <- mx  # Current sample size from input gamma
        d_current <- n_current * total_d1  # Current events
        delta_log <- abs(log(hr))
        zbeta_calc <- sqrt(d_current) * delta_log / sqrt(var_schoen) - zalpha
        power_val <- pnorm(zbeta_calc)
        n <- 1  # Multiplier is 1 for power calculation
      } else {
        d_req <- (zalpha + zbeta)^2 * var_schoen / (log(hr)^2)
        n_sample <- d_req / total_d1
        n <- n_sample / mx
      }
    } else {
      # Multi-stratum: use stratified variance (Schoenfeld uses same var for H0 and H1)
      # For H0 variance with ratio != 1: take total events under H1 and divide by ratio
      # Experimental gets ratio/(ratio+1) proportion, control gets 1/(ratio+1)
      total_events_H1 <- eDC$d + eDE$d
      eDC0_sch <- total_events_H1 * Qc  # Control events under H0
      eDE0_sch <- total_events_H1 * Qe  # Experimental events under H0
      # var0sch then inverse-variance weighted
      var0sch <- 1 / sum((1 / eDC0_sch + 1 / eDE0_sch)^(-1))
      # H1 variance uses actual H1 events
      var1sch <- var0sch
      delta_log <- abs(log(hr / hr0))
      if (is.null(beta)) {
        # Power calculation matching vignette formula:
        # power = pnorm(-(qnorm(1-alpha)*sqrt(var0) - |delta|)/sqrt(var1))
        # Note: variances are already computed for the given gamma (sample size)
        power_val <- pnorm(-(zalpha * sqrt(var0sch) - delta_log) / sqrt(var1sch))
        n <- 1  # Multiplier is 1 for power calculation
      } else {
        # Sample size calculation for stratified Schoenfeld
        # Use the same formula as other methods: c = ((zalpha*sqrt(var0) + zbeta*sqrt(var1)) / |delta|)^2
        # This gives the multiplier for the accrual rate (same as BernsteinLagakos)
        c_mult <- ((zalpha * sqrt(var0sch) + zbeta * sqrt(var1sch)) / delta_log)^2
        n <- c_mult
      }
    }
  } else if (identical(method, "Freedman")) {
    # Freedman formula with contrast delta
    # Check for stratified populations
    if (is.matrix(lambdaC) && ncol(lambdaC) > 1) {
      stop("Stratified method not available for Freedman method")
    }
    # Check for non-inferiority/super-superiority (hr0 must equal 1)
    if (hr0 != 1) {
      stop("Freedman method only supports superiority testing (hr0 = 1)")
    }
    # Freedman uses control:experimental allocation ratio; convert to exp/control
    # Equivalent to D = (zalpha + zbeta)^2 * ratio * (hr + 1/ratio)^2 / (hr - 1)^2
    delta <- (hr - 1) / (hr + 1 / ratio)
    if (is.null(beta)) {
      # Power calculation: solve for zbeta from current expected events
      # D = ((zalpha + zbeta)^2 * ratio) / delta^2
      d_current <- mx * total_d1
      zbeta_calc <- sqrt(d_current / ratio) * abs(delta) - zalpha
      power_val <- pnorm(zbeta_calc)
      n <- 1
    } else {
      d_req <- ((zalpha + zbeta) / delta)^2 * ratio
      # scale to sample size, then convert to accrual rate multiplier
      n_sample <- d_req / total_d1
      n <- n_sample / mx
    }
  } else if (identical(method, "BernsteinLagakos")) {
    # BernsteinLagakos: 
    # Variance computations use 1/c_events + 1/e_events per stratum.
    # Both H1 and H0 variances use input control event rates for  1/c_events.
    # For H1 variance we compute 1/e_events using the input failure rates 
    # times the input hr; for H0 variance we set experimental rate to control
    # when computing e_events.
    # Both H0 and H1 are inverse-variance weighted across strata
    # Compute eDE0 for BernsteinLagakos H0 variance
    # For superiority (hr0 = 1): use lambdaC for experimental group
    # For non-inferiority/super-superiority (hr0 != 1): use lambdaC * hr0
    if (hr0 == 1) {
      # Superiority: experimental rate = control rate
      eDE0_bl_list <- eEvents(
        lambda = lambdaC, eta = etaE, gamma = gammaE,
        R = R, S = S, T = T, minfup = minfup
      )
    } else {
      # Non-inferiority/super-superiority: experimental rate = control * hr0
      eDE0_bl_list <- eEvents(
        lambda = lambdaC * hr0, eta = etaE, gamma = gammaE,
        R = R, S = S, T = T, minfup = minfup
      )
    }
    eDE0_bl <- eDE0_bl_list$d
    # H1 variance per stratum: 1/c_events + 1/e_events (using H1 events)
    var1_bl_strata <- 1 / eDC$d + 1 / eDE$d
    var1_bl <- 1 / sum(1 / var1_bl_strata)  # Inverse-variance weighted
    # H0 variance per stratum: 1/c_events + 1/e_events (using H0 events)
    # where experimental uses control rates (superiority) or control * hr0
    var0_bl_strata <- 1 / eDC0 + 1 / eDE0_bl
    var0_bl <- 1 / sum(1 / var0_bl_strata)  # Inverse-variance weighted
    delta_log <- abs(log(hr / hr0))
    power_val <- NULL
    if (is.null(beta)) {
      # Power calculation matching vignette formula:
      # power = pnorm(-(qnorm(1-alpha)*sqrt(var0) - |delta|)/sqrt(var1))
      power_val <- pnorm(
        -(zalpha * sqrt(var0_bl) - delta_log) / sqrt(var1_bl)
      )
      # For power calculation, n is the multiplier (1 = no chg from input gamma)
      n <- 1
    } else {
      # n calculation: c = ((zalpha*sqrt(var0) + zbeta*sqrt(var1)) / |delta|)^2
      # This gives the multiplier for the accrual rate
      c_mult <- ((zalpha * sqrt(var0_bl) + zbeta * sqrt(var1_bl)) /
        delta_log)^2
      n <- c_mult
    }
  } else {
    stop("Unsupported method")
  }

  # Determine variable type based on what was computed
  if (is.null(beta)) {
    variable_type <- "Power"
  } else {
    variable_type <- "Accrual rate"
  }
  
  rval <- list(
    alpha = alpha, sided = sided, beta = beta, 
    power = if (is.null(power_val)) 1 - beta else power_val,
    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 = variable_type
  )
  class(rval) <- "nSurv"
  return(rval)
}

# print.nSurv roxy [sinew] ----
#' @rdname nSurv
#' @export
# print.nSurv function [sinew] ----
print.nSurv <- function(x, digits = 3, show_strata = TRUE, ...) {
  if (!inherits(x, "nSurv")) {
    stop("Primary argument must have class nSurv")
  }

  # helpers to preserve input values and compact values
  input_value <- function(nm) {
    if (!is.null(x$inputs) && !is.null(x$inputs[[nm]])) {
      return(compact(x$inputs[[nm]]))
    }
    if (!is.null(x$call) && !is.null(x$call[[nm]])) {
      return(paste(deparse(x$call[[nm]]), collapse = ""))
    }
    nm
  }
  compact <- function(v) {
    if (is.null(v)) return("NULL")
    if (length(v) <= 3) return(paste(round(v, digits), collapse = ", "))
    paste0(
      paste(round(v[1:3], digits), collapse = ", "),
      " ... (len=", length(v), ")"
    )
  }

  cat(
    "nSurv fixed-design summary ",
    "(method=", x$method, "; target=", x$variable, ")\n",
    sep = ""
  )
  power_pct <- if (!is.null(x$beta)) {
    (1 - x$beta) * 100
  } else {
    x$power * 100
  }
  cat(
    sprintf(
      "HR=%.3f vs HR0=%.3f | alpha=%.3f (sided=%d) | power=%.1f%%\n",
      x$hr, x$hr0, x$alpha, x$sided, power_pct
    )
  )
  cat(
    sprintf(
      paste(
        "N=%.1f subjects | D=%.1f events |",
        "T=%.1f study duration |",
        "accrual=%.1f Accrual duration |",
        "minfup=%.1f minimum follow-up |",
        "ratio=%s randomization ratio (experimental/control)\n"
      ),
      x$n, x$d, x$T, x$T - x$minfup, x$minfup,
      paste(x$ratio, collapse = "/")
    )
  )

  if (show_strata && !is.null(x$eDC) && length(x$eDC) > 1) {
    cat("\nExpected events by stratum (H1):\n")
    evt <- data.frame(
      control = x$eDC,
      experimental = x$eDE,
      total = x$eDC + x$eDE
    )
    print(round(evt, digits))
  }

  cat("\nKey inputs (names preserved):\n")
  items <- c("gamma", "R", "lambdaC", "eta", "etaE", "S")
  desc_map <- list(
    gamma = "Accrual rate(s)",
    R = "Accrual rate duration(s)",
    lambdaC = "Control hazard rate(s)",
    eta = "Control dropout rate(s)",
    etaE = "Experimental dropout rate(s)",
    S = "Event and dropout rate duration(s)"
  )
  tab <- data.frame(
    desc = unlist(desc_map[items], use.names = FALSE),
    item = items,
    value = vapply(
      items,
      function(nm) {
        switch(nm,
          lambdaC = compact(x$lambdaC),
          gamma = compact(x$gamma),
          eta = compact(x$etaC),
          etaE = compact(x$etaE),
          R = compact(x$R),
          S = compact(x$S),
          compact(NULL)
        )
      },
      character(1)
    ),
    input = vapply(items, input_value, character(1)),
    stringsAsFactors = FALSE
  )
  print(tab, row.names = FALSE)

  invisible(x)
}

2.5.2 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
)

nSurv fixed-design summary (method=; target=Accrual rate)
HR=0.500 vs HR0=1.000 | alpha=0.025 (sided=1) | power=90.0%
N=429.6 subjects | D=90.1 events | T=2.0 study duration | accrual=0.5 Accrual duration | minfup=1.5 minimum follow-up | ratio=1 randomization ratio (experimental/control)

Key inputs (names preserved):
                               desc    item   value   input
                    Accrual rate(s)   gamma 859.238   gamma
           Accrual rate duration(s)       R     0.5       R
             Control hazard rate(s) lambdaC     0.2 lambdaC
            Control dropout rate(s)     eta     0.1     eta
       Experimental dropout rate(s)    etaE     0.1    etaE
 Event and dropout rate duration(s)       S    NULL       S

x <- getFromNamespace("nSurvival", "gsDesign")(
  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

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.4726
Expected sample size (total):       178.797
Accrual rates:
        [,1]    [,2]    [,3]
[1,] 35.7594 35.7594 17.8797
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):        179.1183
Expected sample size (total):       299.8657
Accrual rates:
       [,1]   [,2]   [,3]
[1,] 1.4808 5.9233 4.4425
[2,] 2.9616 2.2212 7.4041
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 <- -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 + qnorm(beta))
    } else {
      return(pnorm(-zb))
    }
  }
  # compute power
  power <- 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(paste(
        "With minfup = NULL, trial is over-powered for any follow-up duration.",
        "Reduce accrual rates (gamma), increase beta, or adjust assumptions."
      ))
    }
    if (right < 0) {
      stop(paste(
        "With minfup = NULL, trial is under-powered for any follow-up duration.",
        "Increase accrual rates (gamma), decrease beta, or adjust assumptions."
      ))
    }
    y <- 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 {
    left <- KTZ(minfup + .01,
      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
    )
    right <- KTZ(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
    )
    if (is.finite(left) && is.finite(right)) {
      if (left > 0 && right > 0) {
        stop(paste(
          "With T = NULL, trial is over-powered for any accrual duration.",
          "Reduce accrual rates (gamma), increase beta, or adjust assumptions."
        ))
      }
      if (left < 0 && right < 0) {
        stop(paste(
          "With T = NULL, trial is under-powered for any accrual duration.",
          "Increase accrual rates (gamma), decrease beta, or adjust assumptions."
        ))
      }
    }
    y <- 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#' Advanced time-to-event sample size calculation
#'
#' \code{nSurv()} is used to calculate the sample size for a two-arm clinical trial
#' with a time-to-event endpoint under the assumption of proportional hazards.
#' The default method assumes a fixed enrollment duration and fixed trial duration;
#' in this case the required sample size is achieved by increasing enrollment rates.
#' \code{nSurv()} implements the \code{Lachin and Foulkes (1986)} method as default.
#' Schoenfeld (1981), Freedman (1982), and Bernstein and Lagakos (1989) methods
#' are also supported; see Details.
#' \code{gsSurv()} combines \code{nSurv()} with \code{gsDesign()} to derive a
#' group sequential design for a study with a time-to-event endpoint.
#'
#' @details
#' The Lachin and Foulkes method uses both null and alternate hypothesis
#' variances to derive sample size and is extended here to support
#' non-inferiority, super-superiority, and stratified designs.
#' As an alternative, the Kim and Tsiatis (1990) method can be used with fixed
#' enrollment rates and either fixed enrollment duration or fixed minimum
#' follow-up.
#' The Schoenfeld approach uses the asymptotic distribution of the log-rank
#' statistic under the assumption of proportional hazards and local alternatives
#' (i.e., \eqn{\log(HR)} is small). The Freedman approach uses the same
#' asymptotic distribution and, like the Schoenfeld approach, uses just the
#' null hypothesis variance to derive sample size.
#' The Bernstein and Lagakos (1989) approach was derived to compute sample size
#' for a stratified model with a common proportional hazards assumption across
#' strata. Like the Lachin and Foulkes (1986) method, it uses both null and
#' alternate hypothesis variances to compute sample size; however, the null
#' hypothesis variance is computed differently. The Bernstein and Lagakos (1989)
#' approach uses the alternate hypothesis failure rate assumptions for both the
#' control and experimental groups, while the Lachin and Foulkes method uses
#' null hypothesis rates that average the alternate hypothesis failure rates to
#' get similar numbers of expected events under the null and alternate
#' hypotheses. Since the Lachin and Foulkes method has fewer events under the
#' null hypothesis (less statistical information), it calculates less power
#' than the Bernstein and Lagakos method.
#' 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{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 gsSurv nEventsIA tEventsIA 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. 
#' Note that rates corresponding the final time period are extended indefinitely.
#' @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; otherwise enrollment after \code{sum(R)} is 0.
#' @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.
#' The final time period is extended indefinitely for each stratum.
#' @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 method One of
#' \code{"LachinFoulkes"} (default), \code{"Schoenfeld"}, \code{"Freedman"},
#' or \code{"BernsteinLagakos"}. Note: \code{"Schoenfeld"} and \code{"Freedman"}
#' methods only support superiority testing (\code{hr0 = 1}).
#' \code{"Freedman"} does not support stratified populations.
#' @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 ofexpected 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}},
#' \code{\link{gsSurvCalendar}}, \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, survival, sample size, power, stratified, proportional hazards,
#' non-inferiority, super-superiority, Lachin and Foulkes, Schoenfeld, Freedman, Bernstein and Lagakos.
#' @examples
#'
#' # Vary accrual rate gamma to obtain power
#' # T, minfup and R all specified, although R will be adjusted on output
#' # gamma as input will be multiplied in output to achieve desired power
#' # Default method is Lachin and Foulkes
#' x_nsurv <- nSurv(lambdaC = log(2) / 6, R = 10, hr = .5, eta = .001, gamma = 1, 
#'       alpha = 0.02, beta = .15, T = 36, minfup = 12, method = "LachinFoulkes")
#' # Demonstrate print method
#' print(x_nsurv)
#' # Same assumptions for group sequential design
#' x_gs <- gsSurv(k = 4, sfl = gsDesign::sfPower, sflpar = .5, lambdaC = log(2) / 6, hr = .5,
#'       eta = .001, gamma = 1, T = 36, minfup = 12, method = "LachinFoulkes")
#' print(x_gs)
#' # Demonstrate xtable method for gsSurv
#' print(xtable::xtable(x_gs,
#'   footnote = "This is a footnote; note that it can be wide.",
#'   caption = "Caption example for xtable output."
#' ))
#' # Demonstrate nEventsIA method
#' # find expected number of events at time 12 in the above trial
#' nEventsIA(x = x_gs, tIA = 10)
#'
#' # find time at which 1/4 of events are expected
#' tEventsIA(x = x_gs, timing = .25)
#'
#' # Adjust accrual duration R to achieve desired power
#' # Trial duration T input as NULL and will be computed based on 
#' # accrual duration R and minimum follow-up duration minfup
#' # Minimum follow-up duration minfup is specified
#' # We use the Schoenfeld method to compute accrual duration R
#' # Control median survival time is 6
#' nSurv(lambdaC = log(2) / 6, hr = .5, eta = .001, gamma = 6, 
#'       alpha = .025, beta = .1, minfup = 12, T = NULL, method = "Schoenfeld")
#' # Same assumptions for group sequential design
#' gsSurv(k = 4, sfu = gsDesign::sfHSD, sfupar = -4, sfl = gsDesign::sfPower, sflpar = .5, 
#'       lambdaC = log(2) / 6, hr = .5, eta = .001, gamma = 6, 
#'       T = 36, minfup = 12, method = "Schoenfeld") |>
#' print()
#'
#' # Vary minimum follow-up duration minfup to obtain power
#' # Accrual duration R rate gamma are fixed and will not change on output. 
#' # Trial duration T and minimum follow-up minfup are input as NULL 
#' # and will be computed on output. 
#' # We will use the Freedman method to compute sample size
#' # Control median survival time is 6
#' # Often this method will fail as the accrual duration and rate provide too 
#' # little or too much sample size.
#' nSurv(lambdaC = log(2) / 6, hr = .5, eta = .001, gamma = 6, R = 25, 
#' alpha = .025, beta = .1, minfup = NULL, T = NULL, method = "Freedman")
#' # Same assumptions for group sequential design
#' gsSurv(k = 4, sfu = gsDesign::sfHSD, sfupar = -4, sfl = gsDesign::sfPower, sflpar = .5, 
#'       lambdaC = log(2) / 6, hr = .5, eta = .001, gamma = 6, 
#'       T = 36, minfup = 12, method = "Freedman") |>
#' print()
#'
#' # piecewise constant enrollment rates (vary accrual rate to achieve power)
#' # also piecewise constant failure rates
#' # will specify annualized enrollment and failure rates
#' nSurv(
#'   lambdaC = -log(c(.95, .97, .98)), # 5%, 3% and 2% annual failure rates
#'   S = c(1, 1), # 1 year duration for first 2 failure rates, 3rd continues indefinitely
#'   R = c(.25, .25, 1.5), # 2-year enrollment with ramp-up over first 1/2 year
#'   gamma = c(1, 3, 6), # 1, 3 and 6 annualized enrollment rates will be 
#'                       # multiplied by ratio to achieve desired power
#'   hr = .5, eta = -log(1 - .01), # 1% annual censoring rate
#'   minfup = 3, T = 5, # 5-year trial duration with 3-year minimum follow-up
#'   alpha = .025, beta = .1, method = "LachinFoulkes")
#' # Same assumptions for group sequential design
#' gsSurv(k = 4, sfu = gsDesign::sfHSD, sfupar = -4, sfl = gsDesign::sfPower, sflpar = .5, 
#'       lambdaC = -log(c(.95, .97, .98)), # 5%, 3% and 2% annual failure rates
#'       S = c(1, 1), # 1 year duration for first 2 failure rates, 3rd continues indefinitely
#'       R = c(.25, .25, 1.5), # 2-year enrollment with ramp-up over first 1/2 year
#'       gamma = c(1, 3, 6), # 1, 3 and 6 annualized enrollment rates will be 
#'                       # multiplied by ratio to achieve desired power
#'       hr = .5, eta = -log(1 - .01), # 1% annual censoring rate
#'       minfup = 3, T = 5, # 5-year trial duration with 3-year minimum follow-up
#'       alpha = .025, beta = .1, method = "LachinFoulkes") |>
#' print()
#'
#' # combine it all: 2 strata, 2 failure rate periods
#' # Note that method = "LachinFoulkes" may be preferred here
#' 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,
#'   alpha = .025, beta = .1, method = "BernsteinLagakos"
#' )
#' # Same assumptions for group sequential design
#' gsSurv(k = 4, sfu = gsDesign::sfHSD, sfupar = -4, sfl = gsDesign::sfPower, sflpar = .5, 
#'       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,
#'       alpha = .025, beta = .1, method = "BernsteinLagakos") |>
#' print()
#'
#' # Example to compute power for a fixed design.
#' # Trial duration T, minimum follow-up minfup and accrual duration R are all
#' # specified and will not change on output. 
#' # beta=NULL will compute power and output will be the same as if beta were specified.
#' # This option is not available for group sequential designs.
#' nSurv(lambdaC = log(2) / 6, hr = .5, eta = .001, gamma = 6, R = 25, 
#'       alpha = .025, beta = NULL, minfup = 12, T = 36, method = "LachinFoulkes") |>
#' print()
#' @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,
                  method = c("LachinFoulkes", "Schoenfeld", "Freedman", "BernsteinLagakos")) {
  method <- match.arg(method)
  input_vals <- list(
    gamma = gamma,
    R = R,
    lambdaC = lambdaC,
    eta = eta,
    etaE = etaE,
    S = S
  )
  # Validate ratio is a single positive scalar
  if (!is.numeric(ratio) || length(ratio) != 1 || ratio <= 0) {
    stop("ratio must be a single positive scalar")
  }
  if (!is.null(S)) {
    if (!is.numeric(S) || any(is.na(S)) || any(!is.finite(S)) || any(S <= 0)) {
      stop("S must be a numeric vector of positive values")
    }
  }
  if (is.null(R)) {
    stop("R must be specified and cannot be NULL")
  }
  if (is.null(beta) && (is.null(T) || is.null(minfup))) {
    stop("When beta is NULL, R, T, and minfup must all be specified")
  }
  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)
  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 {
    # Use LFPWE when T and minfup are both provided
    # This supports method argument for both power (beta = NULL) and 
    # sample size calculations
    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, method = method
    )
  }

  nameR <- nameperiod(cumsum(xx$R))
  stratnames <- paste("Stratum", seq_len(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
  xx$method <- method
  xx$call <- match.call()
  xx$inputs <- input_vals
  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.
method: fixed-design variance/event formulation (default "LachinFoulkes"; alternatives "Schoenfeld", "Freedman", "BernsteinLagakos"); applies per stratum, unequal randomization, and NI designs.

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.
method: As input (variance/event formulation).
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.