R/gsSurv-nSurv.R

In gsDesign: Group Sequential Design

# nSurv roxy [sinew] ----
#' 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,
#'   \code{vignette("SpendingFunctionOverview")}, 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.
#'   Note that \code{sfupar} can be specified as a positive scalar for
#'   \code{sfLDOF} for a generalized O'Brien-Fleming spending function.
#' @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 (>= 1 and <= 80) controlling the number of numerical
#'   integration grid points. Default is 18, as recommended by Jennison and
#'   Turnbull (2000). Grid points are spread out in the tails for accurate
#'   probability calculations. Larger values provide more grid points and greater
#'   accuracy but slow down computation. Jennison and Turnbull (p. 350) note an
#'   accuracy of \eqn{10^{-6}} with \code{r = 16}. This parameter is normally
#'   not changed by users.
#' @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 show_gsDesign Logical; for \code{print.gsSurv()}, include gsDesign details.
#' @param show_strata Logical; for \code{print.gsSurv()}, include strata summaries.
#' @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 the same structure as \code{nEventsIA(..., simple=TRUE)} when
#'
#' @export
#'
#' @author Keaven Anderson \email{keaven_anderson@@merck.com}
#'
#' @seealso \code{\link{gsBoundSummary}}, \code{\link{xprint}},
#'   \code{\link{gsSurvCalendar}}, \link{gsDesign-package},
#'   \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
#' @keywords survival
#' @keywords sample size
#' @keywords power
#' @keywords stratified
#' @keywords proportional hazards
#' @keywords non-inferiority
#' @keywords super-superiority
#' @keywords Lachin and Foulkes
#' @keywords Schoenfeld
#' @keywords Freedman
#' @keywords 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()
# nSurv function [sinew] ----
nSurv <- function(
  lambdaC = log(2) / 6, hr = .6, hr0 = 1, eta = 0, etaE = NULL,
  gamma = 1, R = 12, S = NULL, T = 18, minfup = 6, 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)
}

# 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)
}