Package 'gsDesign'

Convert a summary table object to a gt object

Description

as_gt() is deprecated in favor of lt(), which produces a lightweight HTML table without the heavy gt dependency. as_gt() is kept for one release so existing code that customizes the output with gt functions keeps working; it still returns a gt_tbl object and requires gt to be installed. New code should use lt(); see lt-methods for the available arguments, which mirror those of as_gt().

Usage

as_gt(x, ...)

## S3 method for class 'gsBinomialExactTable'
as_gt(
  x,
  ...,
  title = "Operating Characteristics for the Truncated SPRT Design",
  subtitle = "Assumes trial evaluated sequentially after each response",
  theta_label = gt::html("Underlying<br>response rate"),
  bound_label = c("Futility bound", "Efficacy bound"),
  prob_decimals = 2,
  en_decimals = 1,
  rr_decimals = 0
)

Arguments

x	Object to be converted.
...	Other parameters that may be specific to the object.
title	Table title.
subtitle	Table subtitle.
theta_label	Label for theta.
bound_label	Label for bounds.
prob_decimals	Number of decimal places for probability of crossing.
en_decimals	Number of decimal places for expected number of observations when bound is crossed or when trial ends without crossing.
rr_decimals	Number of decimal places for response rates.

Value

A gt_tbl object that may be extended by overloaded versions of as_gt().

See Also

lt(), lt-methods

---

Save a summary table object as an RTF file

Description

Convert and save the summary table object created with as_table to an RTF file using r2rtf; currently only implemented for gsBinomialExact.

Usage

as_rtf(x, ...)

## S3 method for class 'gsBinomialExactTable'
as_rtf(
  x,
  file,
  ...,
  title = paste("Operating Characteristics by Underlying Response Rate for",
    "Exact Binomial Group Sequential Design"),
  theta_label = "Underlying Response Rate",
  response_outcome = TRUE,
  bound_label = if (response_outcome) c("Futility Bound", "Efficacy Bound") else
    c("Efficacy Bound", "Futility Bound"),
  en_label = "Expected Sample Sizes",
  prob_decimals = 2,
  en_decimals = 1,
  rr_decimals = 0
)

## S3 method for class 'gsBoundSummary'
as_rtf(
  x,
  file,
  ...,
  title = "Boundary Characteristics for Group Sequential Design",
  footnote_p_onesided = "one-side p-value for experimental better than control",
  footnote_appx_effect_at_bound = NULL,
  footnote_p_cross_h0 = "Cumulative type I error assuming binding futility bound",
  footnote_p_cross_h1 = "Cumulative power under the alternate hypothesis",
  footnote_specify = NULL,
  footnote_text = NULL
)

Arguments

x	Object to be saved as RTF file.
...	Other parameters that may be specific to the object.
file	File path for the output.
title	Title of the report.
theta_label	Label for theta.
response_outcome	Logical values indicating if the outcome is response rate (TRUE) or failure rate (FALSE). The default value is TRUE.
bound_label	Label for bounds. If the outcome is response rate, then the label is "Futility bound" and "Efficacy bound". If the outcome is failure rate, then the label is "Efficacy bound" and "Futility bound".
en_label	Label for expected number.
prob_decimals	Number of decimal places for probability of crossing.
en_decimals	Number of decimal places for expected number of observations when bound is crossed or when trial ends without crossing.
rr_decimals	Number of decimal places for response rates.
footnote_p_onesided	Footnote for one-side p-value.
footnote_appx_effect_at_bound	Footnote for approximate effect treatment at bound.
footnote_p_cross_h0	Footnote for cumulative type I error.
footnote_p_cross_h1	Footnote for cumulative power under the alternate hypothesis.
footnote_specify	Vector of string to put footnote super script.
footnote_text	Vector of string of footnote text.

Details

Currently only implemented for gsBinomialExact objects. Creates a table to summarize an object. For gsBinomialExact, this summarized operating characteristics across a range of effect sizes.

Value

as_rtf() returns the input x invisibly.

See Also

vignette("binomialSPRTExample")

Examples

# as_rtf for gsBinomialExact
zz <- gsBinomialExact(
  k = 3, theta = seq(0.1, 0.9, 0.1), n.I = c(12, 24, 36),
  a = c(-1, 0, 11), b = c(5, 9, 12)
)
zz |>
  as_table() |>
  as_rtf(
    file = tempfile(fileext = ".rtf"),
    title = "Power/Type I Error and Expected Sample Size for a Group Sequential Design"
  )

safety_design <- binomialSPRT(
  p0 = .04, p1 = .1, alpha = .04, beta = .2, minn = 4, maxn = 75
)
safety_power <- gsBinomialExact(
  k = length(safety_design$n.I),
  theta = seq(.02, .16, .02),
  n.I = safety_design$n.I,
  a = safety_design$lower$bound,
  b = safety_design$upper$bound
)
safety_power |>
  as_table() |>
  as_rtf(
    file = tempfile(fileext = ".rtf"),
    theta_label = "Underlying\nAE Rate",
    prob_decimals = 3,
    bound_label = c("Low Rate", "High Rate")
  )
# as_rtf for gsBoundSummary
xgs <- gsSurv(lambdaC = .2, hr = .5, eta = .1, T = 2, minfup = 1.5)
gsBoundSummary(xgs, timename = "Year", tdigits = 1) |> as_rtf(file = tempfile(fileext = ".rtf"))

ss <- nSurvival(
  lambda1 = .2, lambda2 = .1, eta = .1, Ts = 2, Tr = .5,
  sided = 1, alpha = .025, ratio = 2
)
xs <- gsDesign(nFixSurv = ss$n, n.fix = ss$nEvents, delta1 = log(ss$lambda2 / ss$lambda1))
gsBoundSummary(xs, logdelta = TRUE, ratio = ss$ratio) |> as_rtf(file = tempfile(fileext = ".rtf"))

xs <- gsDesign(nFixSurv = ss$n, n.fix = ss$nEvents, delta1 = log(ss$lambda2 / ss$lambda1))
gsBoundSummary(xs, logdelta = TRUE, ratio = ss$ratio) |>
  as_rtf(file = tempfile(fileext = ".rtf"),
  footnote_specify = "Z",
  footnote_text = "Z-Score")

---

Create a summary table

Description

Create a tibble to summarize an object; currently only implemented for gsBinomialExact.

Usage

as_table(x, ...)

## S3 method for class 'gsBinomialExact'
as_table(x, ...)

Arguments

x	Object to be summarized.
...	Other parameters that may be specific to the object.

Details

Currently only implemented for gsBinomialExact objects. Creates a table to summarize an object. For gsBinomialExact, this summarized operating characteristics across a range of effect sizes.

Value

A tibble which may have an extended class to enable further output processing.

See Also

vignette("binomialSPRTExample")

Examples

b <- binomialSPRT(p0 = .1, p1 = .35, alpha = .08, beta = .2, minn = 10, maxn = 25)
b_power <- gsBinomialExact(
  k = length(b$n.I), theta = seq(.1, .45, .05), n.I = b$n.I,
  a = b$lower$bound, b = b$upper$bound
)
b_power |>
  as_table() |>
  lt()

---

Power Table for Binomial Tests

Description

Creates a power table for binomial tests with various control group response rates and treatment effects. The function can compute power and Type I error either analytically or through simulation. With large simulations, the function is still fast and can produce exact power values to within simulation error.

Usage

binomialPowerTable(
  pC = c(0.8, 0.9, 0.95),
  delta = seq(-0.05, 0.05, 0.025),
  n = 70,
  ratio = 1,
  alpha = 0.025,
  delta0 = 0,
  scale = "Difference",
  failureEndpoint = TRUE,
  simulation = FALSE,
  nsim = 1e+06,
  adj = 0,
  chisq = 0
)

Arguments

pC	Vector of control group response rates.
delta	Vector of treatment effects (differences in response rates).
n	Total sample size.
ratio	Ratio of experimental to control sample size.
alpha	Type I error rate.
delta0	Non-inferiority margin.
scale	Scale for the test ("Difference", "RR", or "OR").
failureEndpoint	Logical indicating if the endpoint is a failure (TRUE) or success (FALSE).
simulation	Logical indicating whether to use simulation (TRUE) or analytical (FALSE) power calculation.
nsim	Number of simulations to run when simulation = TRUE.
adj	Use continuity correction for the testing (default is 0; only used if simulation = TRUE).
chisq	Chi-squared value for the test (default is 0; only used if simulation = TRUE).

Details

The function binomialPowerTable() creates a grid of all combinations of control group response rates and treatment effects. All out of range values (i.e., where the experimental group response rate is not between 0 and 1) are filtered out. For each combination, it computes the power either analytically using nBinomial() or through simulation using simBinomial(). When using simulation, the simPowerBinomial() function (not exported) is called internally to perform the simulations. Assuming $p$ is the true probability of a positive test, the simulation standard error is

$\text{SE} = \sqrt{p(1 - p) / \text{nsim}}.$

For example, when approximating an underlying Type I error rate of 0.025, the simulation standard error is 0.000156 with 1000000 simulations and the approximated power 95 is 0.025 +/- 1.96 * SE = 0.025 +/- 0.000306.

Value

A data frame containing:

pC

Control group response or failure rate.

delta

Treatment effect.

pE

Experimental group response or failure rate.

Power

Power for the test (asymptotic or simulated).

See Also

nBinomial, simBinomial

Examples

# Create a power table with analytical power calculation
power_table <- binomialPowerTable(
  pC = c(0.8, 0.9),
  delta = seq(-0.05, 0.05, 0.025),
  n = 70
)

# Create a power table with simulation
power_table_sim <- binomialPowerTable(
  pC = c(0.8, 0.9),
  delta = seq(-0.05, 0.05, 0.025),
  n = 70,
  simulation = TRUE,
  nsim = 10000
)

---

Utility functions to verify variable properties

Description

Utility functions to verify an objects's properties including whether it is a scalar or vector, the class, the length, and (if numeric) whether the range of values is on a specified interval. Additionally, the checkLengths function can be used to ensure that all the supplied inputs have equal lengths.

isInteger is similar to is.integer except that isInteger(1) returns TRUE whereas is.integer(1) returns FALSE.

checkScalar is used to verify that the input object is a scalar as well as the other properties specified above.

checkVector is used to verify that the input object is an atomic vector as well as the other properties as defined above.

checkRange is used to check whether the numeric input object's values reside on the specified interval. If any of the values are outside the specified interval, a FALSE is returned.

checkLength is used to check whether all of the supplied inputs have equal lengths.

Usage

checkLengths(..., allowSingle = FALSE)

checkRange(
  x,
  interval = 0:1,
  inclusion = c(TRUE, TRUE),
  varname = deparse(substitute(x)),
  tol = 0
)

checkScalar(x, isType = "numeric", ...)

checkVector(x, isType = "numeric", ..., length = NULL)

isInteger(x)

Arguments

...	For the checkScalar and checkVector functions, this input represents additional arguments sent directly to the checkRange function. For the checkLengths function, this input represents the arguments to check for equal lengths.
allowSingle	logical flag. If TRUE, arguments that are vectors comprised of a single element are not included in the comparative length test in the checkLengths function. Partial matching on the name of this argument is not performed so you must specify 'allowSingle' in its entirety in the call.
x	any object.
interval	two-element numeric vector defining the interval over which the input object is expected to be contained. Use the inclusion argument to define the boundary behavior.
inclusion	two-element logical vector defining the boundary behavior of the specified interval. A TRUE value denotes inclusion of the corresponding boundary. For example, if interval=c(3,6) and inclusion=c(FALSE,TRUE), then all the values of the input object are verified to be on the interval (3,6].
varname	character string defining the name of the input variable as sent into the function by the caller. This is used primarily as a mechanism to specify the name of the variable being tested when checkRange is being called within a function.
tol	numeric scalar defining the tolerance to use in testing the intervals of the checkRange function.
isType	character string defining the class that the input object is expected to be.
length	integer specifying the expected length of the object in the case it is a vector. If length=NULL, the default, then no length check is performed.

Value

isInteger: Boolean value as checking result Other functions have no return value, called for side effects

Examples

# check whether input is an integer
isInteger(1)
isInteger(1:5)
try(isInteger("abc")) # expect error

# check whether input is an integer scalar
checkScalar(3, "integer")

# check whether input is an integer scalar that resides
# on the interval on [3, 6]. Then test for interval (3, 6].
checkScalar(3, "integer", c(3, 6))
try(checkScalar(3, "integer", c(3, 6), c(FALSE, TRUE))) # expect error

# check whether the input is an atomic vector of class numeric,
# of length 3, and whose value all reside on the interval [1, 10)
x <- c(3, pi, exp(1))
checkVector(x, "numeric", c(1, 10), c(TRUE, FALSE), length = 3)

# do the same but change the expected length; expect error
try(checkVector(x, "numeric", c(1, 10), c(TRUE, FALSE), length = 2))

# create faux function to check input variable
foo <- function(moo) checkVector(moo, "character")
foo(letters)
try(foo(1:5)) # expect error with function and argument name in message

# check for equal lengths of various inputs
checkLengths(1:2, 2:3, 3:4)
try(checkLengths(1, 2, 3, 4:5)) # expect error

# check for equal length inputs but ignore single element vectors
checkLengths(1, 2, 3, 4:5, 7:8, allowSingle = TRUE)

---

Testing, Confidence Intervals, Sample Size and Power for Comparing Two Binomial Rates

Description

Support is provided for sample size estimation, power, testing, confidence intervals and simulation for fixed sample size trials (that is, not group sequential or adaptive) with two arms and binary outcomes. Both superiority and non-inferiority trials are considered. While all routines default to comparisons of risk-difference, options to base computations on risk-ratio and odds-ratio are also included.

nBinomial() computes sample size or power using the method of Farrington and Manning (1990) for a trial to test the difference between two binomial event rates. The routine can be used for a test of superiority or non-inferiority. For a design that tests for superiority nBinomial() is consistent with the method of Fleiss, Tytun, and Ury (but without the continuity correction) to test for differences between event rates. This routine is consistent with the Hmisc package routines bsamsize and bpower for superiority designs. Vector arguments allow computing sample sizes for multiple scenarios for comparative purposes.

testBinomial() computes a Z- or Chi-square-statistic that compares two binomial event rates using the method of Miettinen and Nurminen (1980). This can be used for superiority or non-inferiority testing. Vector arguments allow easy incorporation into simulation routines for fixed, group sequential and adaptive designs.

ciBinomial() computes confidence intervals for 1) the difference between two rates, 2) the risk-ratio for two rates or 3) the odds-ratio for two rates. This procedure provides inference that is consistent with testBinomial() in that the confidence intervals are produced by inverting the testing procedures in testBinomial(). The Type I error alpha input to ciBinomial is always interpreted as 2-sided.

simBinomial() performs simulations to estimate the power for a Miettinen and Nurminen (1985) test comparing two binomial rates for superiority or non-inferiority. As noted in documentation for bpower.sim() in the HMisc package, by using testBinomial() you can see that the formulas without any continuity correction are quite accurate. In fact, Type I error for a continuity-corrected test is significantly lower (Gordon and Watson, 1996) than the nominal rate. Thus, as a default no continuity corrections are performed.

varBinomial computes blinded estimates of the variance of the estimate of 1) event rate differences, 2) logarithm of the risk ratio, or 3) logarithm of the odds ratio. This is intended for blinded sample size re-estimation for comparative trials with a binary outcome.

Testing is 2-sided when a Chi-square statistic is used and 1-sided when a Z-statistic is used. Thus, these 2 options will produce substantially different results, in general. For non-inferiority, 1-sided testing is appropriate.

You may wish to round sample sizes up using ceiling().

Farrington and Manning (1990) begin with event rates p1 and p2 under the alternative hypothesis and a difference between these rates under the null hypothesis, delta0. From these values, actual rates under the null hypothesis are computed, which are labeled p10 and p20 when outtype=3. The rates p1 and p2 are used to compute a variance for a Z-test comparing rates under the alternative hypothesis, while p10 and p20 are used under the null hypothesis. This computational method is also used to estimate variances in varBinomial() based on the overall event rate observed and the input treatment difference specified in delta0.

Sample size with scale="Difference" produces an error if p1-p2=delta0. Normally, the alternative hypothesis under consideration would be p1-p2-delta0$>0$. However, the alternative can have p1-p2-delta0$<0$.

Usage

ciBinomial(x1, x2, n1, n2, alpha = 0.05, adj = 0, scale = "Difference")

nBinomial(
  p1,
  p2,
  alpha = 0.025,
  beta = 0.1,
  delta0 = 0,
  ratio = 1,
  sided = 1,
  outtype = 1,
  scale = "Difference",
  n = NULL
)

simBinomial(
  p1,
  p2,
  n1,
  n2,
  delta0 = 0,
  nsim = 10000,
  chisq = 0,
  adj = 0,
  scale = "Difference"
)

testBinomial(
  x1,
  x2,
  n1,
  n2,
  delta0 = 0,
  chisq = 0,
  adj = 0,
  scale = "Difference",
  tol = 1e-11
)

varBinomial(x, n, delta0 = 0, ratio = 1, scale = "Difference")

Arguments

x1	Number of “successes” in the control group
x2	Number of “successes” in the experimental group
n1	Number of observations in the control group
n2	Number of observations in the experimental group
alpha	type I error; see sided below to distinguish between 1- and 2-sided tests
adj	With adj=1, the standard variance with a continuity correction is used for a Miettinen and Nurminen test statistic This includes a factor of $n / (n - 1)$ where $n$ is the total sample size. If adj is not 1, this factor is not applied. The default is adj=0 since nominal Type I error is generally conservative with adj=1 (Gordon and Watson, 1996).
scale	“Difference”, “RR”, “OR”; see the scale parameter documentation above and Details. This is a scalar argument.
p1	event rate in group 1 under the alternative hypothesis
p2	event rate in group 2 under the alternative hypothesis
beta	type II error
delta0	A value of 0 (the default) always represents no difference between treatment groups under the null hypothesis. delta0 is interpreted differently depending on the value of the parameter scale. If scale="Difference" (the default), delta0 is the difference in event rates under the null hypothesis (p10 - p20). If scale="RR", delta0 is the logarithm of the relative risk of event rates (p10 / p20) under the null hypothesis. If scale="LNOR", delta0 is the difference in natural logarithm of the odds-ratio under the null hypothesis log(p10 / (1 - p10)) - log(p20 / (1 - p20)).
ratio	sample size ratio for group 2 divided by group 1
sided	2 for 2-sided test, 1 for 1-sided test
outtype	nBinomial only; 1 (default) returns total sample size; 2 returns a data frame with sample size for each group (n1, n2; if n is not input as NULL, power is returned in Power; 3 returns a data frame with total sample size (n), sample size in each group (n1, n2), Type I error (alpha), 1 or 2 (sided, as input), Type II error (beta), power (Power), null and alternate hypothesis standard deviations (sigma0, sigma1), input event rates (p1, p2), null hypothesis difference in treatment group means (delta0) and null hypothesis event rates (p10, p20).
n	If power is to be computed in nBinomial(), input total trial sample size in n; this may be a vector. This is also the sample size in varBinomial, in which case the argument must be a scalar.
nsim	The number of simulations to be performed in simBinomial()
chisq	An indicator of whether or not a chi-square (as opposed to Z) statistic is to be computed. If delta0=0 (default), the difference in event rates divided by its standard error under the null hypothesis is used. Otherwise, a Miettinen and Nurminen chi-square statistic for a 2 x 2 table is used.
tol	Default should probably be used; this is used to deal with a rounding issue in interim calculations
x	Number of “successes” in the combined control and experimental groups.

Value

testBinomial() and simBinomial() each return a vector of either Chi-square or Z test statistics. These may be compared to an appropriate cutoff point (e.g., qnorm(.975) for normal or qchisq(.95,1) for chi-square).

ciBinomial() returns a data frame with 1 row with a confidence interval; variable names are lower and upper.

varBinomial() returns a vector of (blinded) variance estimates of the difference of event rates (scale="Difference"), logarithm of the odds-ratio (scale="OR") or logarithm of the risk-ratio (scale="RR").

With the default outtype=1, nBinomial() returns a vector of total sample sizes is returned. With outtype=2, nBinomial() returns a data frame containing two vectors n1 and n2 containing sample sizes for groups 1 and 2, respectively; if n is input, this option also returns the power in a third vector, Power. With outtype=3, nBinomial() returns a data frame with the following columns:

n	A vector with total samples size required for each event rate comparison specified
n1	A vector of sample sizes for group 1 for each event rate comparison specified
n2	A vector of sample sizes for group 2 for each event rate comparison specified
alpha	As input
sided	As input
beta	As input; if n is input, this is computed
Power	If n=NULL on input, this is 1-beta; otherwise, the power is computed for each sample size input
sigma0	A vector containing the standard deviation of the treatment effect difference under the null hypothesis times sqrt(n) when scale="Difference" or scale="OR"; when scale="RR", this is the standard deviation time sqrt(n) for the numerator of the Farrington-Manning test statistic x1-exp(delta0)*x2.
sigma1	A vector containing the values as sigma0, in this case estimated under the alternative hypothesis.
p1	As input
p2	As input
p10	group 1 event rate used for null hypothesis
p20	group 2 event rate used for null hypothesis

Author(s)

Keaven Anderson [email protected]

References

Farrington, CP and Manning, G (1990), Test statistics and sample size formulae for comparative binomial trials with null hypothesis of non-zero risk difference or non-unity relative risk. Statistics in Medicine; 9: 1447-1454.

Fleiss, JL, Tytun, A and Ury (1980), A simple approximation for calculating sample sizes for comparing independent proportions. Biometrics;36:343-346.

Gordon, I and Watson R (1985), The myth of continuity-corrected sample size formulae. Biometrics; 52: 71-76.

Miettinen, O and Nurminen, M (1985), Comparative analysis of two rates. Statistics in Medicine; 4 : 213-226.

See Also

Normal,uniroot

Examples

# Compute z-test test statistic comparing 39/500 to 13/500
# use continuity correction in variance
x <- testBinomial(x1 = 39, x2 = 13, n1 = 500, n2 = 500, adj = 1)
x
pnorm(x, lower.tail = FALSE)

# Compute with unadjusted variance
x0 <- testBinomial(x1 = 39, x2 = 23, n1 = 500, n2 = 500)
x0
pnorm(x0, lower.tail = FALSE)

# Perform 50k simulations to test validity of the above
# asymptotic p-values
# (you may want to perform more to reduce standard error of estimate)
sum(as.double(x0) <=
  simBinomial(p1 = .078, p2 = .078, n1 = 500, n2 = 500, nsim = 10000)) / 10000
sum(as.double(x0) <=
  simBinomial(p1 = .052, p2 = .052, n1 = 500, n2 = 500, nsim = 10000)) / 10000

# Perform a non-inferiority test to see if p2=400 / 500 is within 5% of
# p1=410 / 500 use a z-statistic with unadjusted variance
x <- testBinomial(x1 = 410, x2 = 400, n1 = 500, n2 = 500, delta0 = -.05)
x
pnorm(x, lower.tail = FALSE)

# since chi-square tests equivalence (a 2-sided test) rather than
# non-inferiority (a 1-sided test),
# the result is quite different
pchisq(testBinomial(
  x1 = 410, x2 = 400, n1 = 500, n2 = 500, delta0 = -.05,
  chisq = 1, adj = 1
), 1, lower.tail = FALSE)

# now simulate the z-statistic witthout continuity corrected variance
sum(qnorm(.975) <=
  simBinomial(p1 = .8, p2 = .8, n1 = 500, n2 = 500, nsim = 100000)) / 100000

# compute a sample size to show non-inferiority
# with 5% margin, 90% power
nBinomial(p1 = .2, p2 = .2, delta0 = .05, alpha = .025, sided = 1, beta = .1)

# assuming a slight advantage in the experimental group lowers
# sample size requirement
nBinomial(p1 = .2, p2 = .19, delta0 = .05, alpha = .025, sided = 1, beta = .1)

# compute a sample size for comparing 15% vs 10% event rates
# with 1 to 2 randomization
nBinomial(p1 = .15, p2 = .1, beta = .2, ratio = 2, alpha = .05)

# now look at total sample size using 1-1 randomization
n <- nBinomial(p1 = .15, p2 = .1, beta = .2, alpha = .05)
n
# check if inputing sample size returns the desired power
nBinomial(p1 = .15, p2 = .1, beta = .2, alpha = .05, n = n)

# re-do with alternate output types
nBinomial(p1 = .15, p2 = .1, beta = .2, alpha = .05, outtype = 2)
nBinomial(p1 = .15, p2 = .1, beta = .2, alpha = .05, outtype = 3)

# look at power plot under different control event rate and
# relative risk reductions
library(dplyr)
library(ggplot2)
p1 <- seq(.075, .2, .000625)
len <- length(p1)
p2 <- c(p1 * .75, p1 * 2/3, p1 * .6, p1 * .5)
Reduction <- c(rep("25 percent", len), rep("33 percent", len), 
               rep("40 percent", len), rep("50 percent", len))
df <- tibble(p1 = rep(p1, 4), p2, Reduction) |>
  mutate(`Sample size` = nBinomial(p1, p2, beta = .2, alpha = .025, sided = 1))
ggplot(df, aes(x = p1, y = `Sample size`, col = Reduction)) + 
  geom_line() + 
  xlab("Control group event rate") +
  ylim(0,6000) +
  ggtitle("Binomial sample size computation for 80 pct power")

# compute blinded estimate of treatment effect difference
x1 <- rbinom(n = 1, size = 100, p = .2)
x2 <- rbinom(n = 1, size = 200, p = .1)
# blinded estimate of risk difference variance
varBinomial(x = x1 + x2, n = 300, ratio = 2, delta0 = 0)
# blinded estimate of log-risk-ratio variance
varBinomial(x = x1 + x2, n = 300, ratio = 2, delta0 = 0, scale = "RR")
# blinded estimate of log-odds-ratio variance
varBinomial(x = x1 + x2, n = 300, ratio = 2, delta0 = 0, scale = "OR")

---

Exact confidence intervals for vaccine or prevention efficacy

Description

Computes a fixed look Clopper–Pearson confidence interval for vaccine or prevention efficacy after conditioning on the total number of events.

Usage

ciBinomialExact(
  x,
  n,
  ratio = 1,
  conf.level = 0.95,
  alternative = c("two.sided", "lower", "upper")
)

Arguments

x	Number of events in the experimental group.
n	Total number of events.
ratio	Experimental-to-control randomization ratio.
conf.level	Confidence level.
alternative	Character string specifying a two-sided interval, a lower confidence bound, or an upper confidence bound for efficacy.

Details

Conditional on n, x is binomial with event probability p in the experimental group. Clopper–Pearson limits for p are transformed using efficacy = 1 - p / (ratio * (1 - p)). Because this transformation is decreasing, the probability limits are reversed on the efficacy scale. The same calculation applies to vaccine efficacy (VE) and prevention efficacy (PE).

Value

A one-row data frame containing the observed counts, randomization ratio, efficacy estimate, confidence limits, and confidence level.

See Also

repeatedCIBinomialExact, sequentialCIBinomialExact

Examples

ciBinomialExact(x = 16, n = 78, ratio = 3)

---

Sample size re-estimation based on conditional power

Description

ssrCP() adapts 2-stage group sequential designs to 2-stage sample size re-estimation designs based on interim analysis conditional power. This is a simple case of designs developed by Lehmacher and Wassmer, Biometrics (1999). The conditional power designs of Bauer and Kohne (1994), Proschan and Hunsberger (1995), Cui, Hung and Wang (1999) and Liu and Chi (2001), Gao, Ware and Mehta (2008), and Mehta and Pocock (2011). Either the estimated treatment effect at the interim analysis or any chosen effect size can be used for sample size re-estimation. If not done carefully, these designs can be very inefficient. It is probably a good idea to compare any design to a simpler group sequential design; see, for example, Jennison and Turnbull (2003). The a assumes a small Type I error is included with the interim analysis and that the design is an adaptation from a 2-stage group sequential design Related functions include 3 pre-defined combination test functions (z2NC, z2Z, z2Fisher) that represent the inverse normal combination test (Lehmacher and Wassmer, 1999), the sufficient statistic for the complete data, and Fisher's combination test. Power.ssrCP computes unconditional power for a conditional power design derived by ssrCP.

condPower is a supportive routine that also is interesting in its own right; it computes conditional power of a combination test given an interim test statistic, stage 2 sample size and combination test statistic. While the returned data frame should make general plotting easy, the function plot.ssrCP() prints a plot of study sample size by stage 1 outcome with multiple x-axis labels for stage 1 z-value, conditional power, and stage 1 effect size relative to the effect size for which the underlying group sequential design was powered.

Sample size re-estimation using conditional power and an interim estimate of treatment effect was proposed by several authors in the 1990's (see references below). Statistical testing for these original methods was based on combination tests since Type I error could be inflated by using a sufficient statistic for testing at the end of the trial. Since 2000, more efficient variations of these conditional power designs have been developed. Fully optimized designs have also been derived (Posch et all, 2003, Lokhnygina and Tsiatis, 2008). Some of the later conditional power methods allow use of sufficient statistics for testing (Chen, DeMets and Lan, 2004, Gao, Ware and Mehta, 2008, and Mehta and Pocock, 2011).

The methods considered here are extensions of 2-stage group sequential designs that include both an efficacy and a futility bound at the planned interim analysis. A maximum fold-increase in sample size (maxinc)from the supplied group sequential design (x) is specified by the user, as well as a range of conditional power (cpadj) where sample size should be re-estimated at the interim analysis, 1-the targeted conditional power to be used for sample size re-estimation (beta) and a combination test statistic (z2) to be used for testing at the end of the trial. The input value overrun represents incremental enrollment not included in the interim analysis that is not included in the analysis; this is used for calculating the required number of patients enrolled to complete the trial.

Whereas most of the methods proposed have been based on using the interim estimated treatment effect size (default for ssrCP), the variable theta allows the user to specify an alternative; Liu and Chi (2001) suggest that using the parameter value for which the trial was originally powered is a good choice.

Usage

condPower(
  z1,
  n2,
  z2 = z2NC,
  theta = NULL,
  x = gsDesign(k = 2, timing = 0.5, beta = beta),
  ...
)

ssrCP(
  z1,
  theta = NULL,
  maxinc = 2,
  overrun = 0,
  beta = x$beta,
  cpadj = c(0.5, 1 - beta),
  x = gsDesign(k = 2, timing = 0.5),
  z2 = z2NC,
  ...
)

## S3 method for class 'ssrCP'
plot(
  x,
  z1ticks = NULL,
  mar = c(7, 4, 4, 4) + 0.1,
  ylab = "Adapted sample size",
  xlaboffset = -0.2,
  lty = 1,
  col = 1,
  ...
)

z2NC(z1, x, ...)

z2Z(z1, x, n2 = x$n.I[2] - x$n.I[1], ...)

z2Fisher(z1, x, ...)

Power.ssrCP(x, theta = NULL, delta = NULL, r = 18)

Arguments

z1	Scalar or vector with interim standardized Z-value(s). Input of multiple values makes it easy to plot the revised sample size as a function of the interim test statistic.
n2	stage 2 sample size to be used in computing sufficient statistic when combining stage 1 and 2 test statistics.
z2	a combination function that returns a test statistic cutoff for the stage 2 test based in the interim test statistic supplied in z1, the design x and the stage 2 sample size.
theta	If NULL (default), conditional power calculation will be based on estimated interim treatment effect. Otherwise, theta is the standardized effect size used for conditional power calculation. Using the alternate hypothesis treatment effect can be more efficient than the estimated effect size; see Liu and Chi, Biometrics (2001).
x	A group sequential design with 2 stages (k=2) generated by gsDesign. For plot.ssrCP, x is a design returned by ssrCP().
...	Allows passing of arguments that may be needed by the user-supplied function, codez2. In the case of plot.ssrCP(), allows passing more graphical parameters.
maxinc	Maximum fold-increase from planned maximum sample size in underlying group sequential design provided in x.
overrun	The number of patients enrolled before the interim analysis is completed who are not included in the interim analysis.
beta	Targeted Type II error (1 - targeted conditional power); used for conditional power in sample size reestimation.
cpadj	Range of values strictly between 0 and 1 specifying the range of interim conditional power for which sample size re-estimation is to be performed. Outside of this range, the sample size supplied in x is used.
z1ticks	Test statistic values at which tick marks are to be made on x-axis; automatically calculated under default of NULL
mar	Plot margins; see help for par
ylab	y-axis label
xlaboffset	offset on x-axis for printing x-axis labels
lty	line type for stage 2 sample size
col	line color for stage 2 sample size
delta	Natural parameter values for power calculation; see gsDesign for a description of how this is related to theta.
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 $10^{-6}$ with r = 16. This parameter is normally not changed by users.

Value

ssrCP returns a list with the following items:

x	As input.
z2fn	As input in z2.
theta	standardize effect size used for conditional power; if NULL is input, this is computed as z1/sqrt(n1) where n1 is the stage 1 sample size.
maxinc	As input.
overrun	As input.
beta	As input.
cpadj	As input.
dat	A data frame containing the input z1 values, computed cutoffs for the standard normal test statistic based solely on stage 2 data (z2), stage 2 sample sizes (n2), stage 2 conditional power (CP), standardize effect size used for conditional power calculation (theta), and the natural parameter value corresponding to theta (delta). The relation between theta and delta is determined by the delta0 and delta1 values from x: delta = delta0 + theta(delta1-delta0).

Author(s)

Keaven Anderson [email protected]

References

Bauer, Peter and Kohne, F., Evaluation of experiments with adaptive interim analyses, Biometrics, 50:1029-1041, 1994.

Chen, YHJ, DeMets, DL and Lan, KKG. Increasing the sample size when the unblinded interim result is promising, Statistics in Medicine, 23:1023-1038, 2004.

Gao, P, Ware, JH and Mehta, C, Sample size re-estimation for adaptive sequential design in clinical trials, Journal of Biopharmaceutical Statistics", 18:1184-1196, 2008.

Jennison, C and Turnbull, BW. Mid-course sample size modification in clinical trials based on the observed treatment effect. Statistics in Medicine, 22:971-993", 2003.

Lehmacher, W and Wassmer, G. Adaptive sample size calculations in group sequential trials, Biometrics, 55:1286-1290, 1999.

Liu, Q and Chi, GY., On sample size and inference for two-stage adaptive designs, Biometrics, 57:172-177, 2001.

Mehta, C and Pocock, S. Adaptive increase in sample size when interim results are promising: A practical guide with examples, Statistics in Medicine, 30:3267-3284, 2011.

See Also

gsDesign

Examples

library(ggplot2)
# quick trick for simple conditional power based on interim z-value, stage 1 and 2 sample size
# assumed treatment effect and final alpha level
# and observed treatment effect
alpha <- .01 # set final nominal significance level
timing <- .6 # set proportion of sample size, events or statistical information at IA
n2 <- 40 # set stage 2 sample size events or statistical information
hr <- .6 # for this example we will derive conditional power based on hazard ratios
n.fix <- nEvents(hr=hr,alpha=alpha) # you could otherwise make n.fix an arbitrary positive value
# this just derives a group sequential design that should not change sample size from n.fix
# due to stringent IA bound
x <- gsDesign(k=2,n.fix=n.fix,alpha=alpha,test.type=1,sfu=sfHSD,
sfupar=-20,timing=timing,delta1=log(hr))
# derive effect sizes for which you wish to compute conditional power
hrpostIA = seq(.4,1,.05)
# in the following, we convert HR into standardized effect size based on the design in x
powr <- condPower(x=x,z1=1,n2=x$n.I[2]-x$n.I[1],theta=log(hrpostIA)/x$delta1*x$theta[2])
ggplot(
  data.frame(
    x = hrpostIA,
    y = condPower(
      x = x,
      z1 = 1,
      n2 = x$n.I[2] - x$n.I[1],
      theta = log(hrpostIA) / x$delta1 * x$theta[2]
    )
  ),
  aes(x = x, y = y)
) +
  geom_line() +
  labs(
    x = "HR post IA",
    y = "Conditional power",
    title = "Conditional power as a function of assumed HR"
  )

# Following is a template for entering parameters for ssrCP
# Natural parameter value null and alternate hypothesis values
delta0 <- 0
delta1 <- 1
# timing of interim analysis for underlying group sequential design
timing <- .5
# upper spending function
sfu <- sfHSD
# upper spending function paramater
sfupar <- -12
# maximum sample size inflation
maxinflation <- 2
# assumed enrollment overrrun at IA
overrun <- 25
# interim z-values for plotting
z <- seq(0, 4, .025)
# Type I error (1-sided)
alpha <- .025
# Type II error for design
beta <- .1
# Fixed design sample size
n.fix <- 100
# conditional power interval where sample
# size is to be adjusted
cpadj <- c(.3, .9)
# targeted Type II error when adapting sample size
betastar <- beta
# combination test (built-in options are: z2Z, z2NC, z2Fisher)
z2 <- z2NC

# use the above parameters to generate design
# generate a 2-stage group sequential design with
x <- gsDesign(
  k = 2, n.fix = n.fix, timing = timing, sfu = sfu, sfupar = sfupar,
  alpha = alpha, beta = beta, delta0 = delta0, delta1 = delta1
)
# extend this to a conditional power design
xx <- ssrCP(
  x = x, z1 = z, overrun = overrun, beta = betastar, cpadj = cpadj,
  maxinc = maxinflation, z2 = z2
)
# plot the stage 2 sample size
plot(xx)
# demonstrate overlays on this plot
# overlay with densities for z1 under different hypotheses
lines(z, 80 + 240 * dnorm(z, mean = 0), col = 2)
lines(z, 80 + 240 * dnorm(z, mean = sqrt(x$n.I[1]) * x$theta[2]), col = 3)
lines(z, 80 + 240 * dnorm(z, mean = sqrt(x$n.I[1]) * x$theta[2] / 2), col = 4)
lines(z, 80 + 240 * dnorm(z, mean = sqrt(x$n.I[1]) * x$theta[2] * .75), col = 5)
axis(side = 4, at = 80 + 240 * seq(0, .4, .1), labels = as.character(seq(0, .4, .1)))
mtext(side = 4, expression(paste("Density for ", z[1])), line = 2)
text(x = 1.5, y = 90, col = 2, labels = expression(paste("Density for ", theta, "=0")))
text(x = 3.00, y = 180, col = 3, labels = expression(paste("Density for ", theta, "=",
 theta[1])))
text(x = 1.00, y = 180, col = 4, labels = expression(paste("Density for ", theta, "=",
 theta[1], "/2")))
text(x = 2.5, y = 140, col = 5, labels = expression(paste("Density for ", theta, "=",
 theta[1], "*.75")))
# overall line for max sample size
nalt <- xx$maxinc * x$n.I[2]
lines(x = par("usr")[1:2], y = c(nalt, nalt), lty = 2)

# compare above design with different combination tests
# use sufficient statistic for final testing
xxZ <- ssrCP(
  x = x, z1 = z, overrun = overrun, beta = betastar, cpadj = cpadj,
  maxinc = maxinflation, z2 = z2Z
)
# use Fisher combination test for final testing
xxFisher <- ssrCP(
  x = x, z1 = z, overrun = overrun, beta = betastar, cpadj = cpadj,
  maxinc = maxinflation, z2 = z2Fisher
)
# combine data frames from these designs
y <- rbind(
  data.frame(xx$dat, Test = "Normal combination"),
  data.frame(xxZ$dat, Test = "Sufficient statistic"),
  data.frame(xxFisher$dat, Test = "Fisher combination")
)
# plot stage 2 statistic required for positive combination test
ggplot2::ggplot(data = y, ggplot2::aes(x = z1, y = z2, col = Test)) + 
ggplot2::geom_line()
# plot total sample size versus stage 1 test statistic
ggplot2::ggplot(data = y, ggplot2::aes(x = z1, y = n2, col = Test)) + 
ggplot2::geom_line()
# check achieved Type I error for sufficient statistic design
Power.ssrCP(x = xxZ, theta = 0)

# compare designs using observed vs planned theta for conditional power
xxtheta1 <- ssrCP(
  x = x, z1 = z, overrun = overrun, beta = betastar, cpadj = cpadj,
  maxinc = maxinflation, z2 = z2, theta = x$delta
)
# combine data frames for the 2 designs
y <- rbind(
  data.frame(xx$dat, "CP effect size" = "Obs. at IA"),
  data.frame(xxtheta1$dat, "CP effect size" = "Alt. hypothesis")
)
# plot stage 2 sample size by design
ggplot2::ggplot(data = y, ggplot2::aes(x = z1, y = n2, col = CP.effect.size)) + 
ggplot2::geom_line()
# compare power and expected sample size
y1 <- Power.ssrCP(x = xx)
y2 <- Power.ssrCP(x = xxtheta1)
# combine data frames for the 2 designs
y3 <- rbind(
  data.frame(y1, "CP effect size" = "Obs. at IA"),
  data.frame(y2, "CP effect size" = "Alt. hypothesis")
)
# plot expected sample size by design and effect size
ggplot2::ggplot(data = y3, ggplot2::aes(x = delta, y = en, col = CP.effect.size)) + 
ggplot2::geom_line() +
ggplot2::xlab(expression(delta)) + 
ggplot2::ylab("Expected sample size")
# plot power by design and effect size
ggplot2::ggplot(data = y3, ggplot2::aes(x = delta, y = Power, col = CP.effect.size)) + 
ggplot2::geom_line() + 
ggplot2::xlab(expression(delta))

---

Expected number of events for a time-to-event study

Description

eEvents() is used to calculate the expected number of events for a population with a time-to-event endpoint. It is based on calculations demonstrated in Lachin and Foulkes (1986) and is fundamental in computations for the sample size method they propose. Piecewise exponential survival and dropout rates are supported as well as piecewise uniform enrollment. A stratified population is allowed. Output is the expected number of events observed given a trial duration and the above rate parameters.

Usage

eEvents(
  lambda = 1,
  eta = 0,
  gamma = 1,
  R = 1,
  S = NULL,
  T = 2,
  Tfinal = NULL,
  minfup = 0,
  digits = 4
)

## S3 method for class 'eEvents'
print(x, digits = 4, ...)

Arguments

lambda	Scalar, vector or matrix of event hazard rates; rows represent time periods while columns represent strata; a vector implies a single stratum.
eta	Scalar, vector or matrix of dropout hazard rates; rows represent time periods while columns represent strata; if entered as a scalar, rate is constant across strata and time periods; if entered as a vector, rates are constant across strata.
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.
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. In eEvents(), recruitment after the final period is assumed to be 0 if T exceeds sum(R); the final period is only extended in designs where T is solved for (e.g., nSurv with T = NULL).
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 length of the number of rows in lambda minus 1, otherwise.
T	Time of analysis; if Tfinal=NULL, this is also the study duration.
Tfinal	Study duration; if NULL, this will be replaced with T on output.
minfup	Time from end of planned enrollment (sum(R) from output value of R) until Tfinal.
digits	Which controls number of digits for printing.
x	An object of class eEvents returned from eEvents().
...	Other arguments that may be passed to the generic print function.

Details

eEvents() produces an object of class eEvents with the number of subjects and events for a set of pre-specified trial parameters, such as accrual duration and follow-up period. The underlying power calculation is based on Lachin and Foulkes (1986) method for proportional hazards assuming a fixed underlying hazard ratio between 2 treatment groups. The method has been extended here to enable designs to test non-inferiority. Piecewise constant enrollment and failure rates are assumed and a stratified population is allowed. See also nSurvival for other Lachin and Foulkes (1986) methods assuming a constant hazard difference or exponential enrollment rate.

print.eEvents() formats the output for an object of class eEvents and returns the input value.

Value

eEvents() and print.eEvents() return an object of class eEvents which contains the following items:

lambda	As input; converted to a matrix on output.
eta	As input; converted to a matrix on output.
gamma	As input.
R	As input.
S	As input.
T	As input.
Tfinal	Planned duration of study.
minfup	As input.
d	Expected number of events.
n	Expected sample size.
digits	As input.

Author(s)

Keaven Anderson [email protected]

References

Lachin JM and Foulkes MA (1986), Evaluation of Sample Size and Power for Analyses of Survival with Allowance for Nonuniform Patient Entry, Losses to Follow-Up, Noncompliance, and Stratification. Biometrics, 42, 507-519.

Bernstein D and Lagakos S (1978), Sample size and power determination for stratified clinical trials. Journal of Statistical Computation and Simulation, 8:65-73.

See Also

vignette("gsDesignPackageOverview"), plot.gsDesign, gsDesign, gsHR, nSurvival

Examples

# 3 enrollment periods, 3 piecewise exponential failure rates
str(eEvents(
  lambda = c(.05, .02, .01), eta = .01, gamma = c(5, 10, 20),
  R = c(2, 1, 2), S = c(1, 1), T = 20
))

# Control group for example from Bernstein and Lagakos (1978)
lamC <- c(1, .8, .5)
n <- eEvents(
  lambda = matrix(c(lamC, lamC * 2 / 3), ncol = 6), eta = 0,
  gamma = matrix(.5, ncol = 6), R = 2, T = 4
)

---

One-Sample Binomial Routines

Description

gsBinomialExact computes power/Type I error and expected sample size for a group sequential design in a single-arm trial with a binary outcome. This can also be used to compare event rates in two-arm studies. The print function has been extended using print.gsBinomialExact to print gsBinomialExact objects. Similarly, a plot function has been extended using plot.gsBinomialExact to plot gsBinomialExact objects.

binomialSPRT computes a truncated binomial sequential probability ratio test (SPRT) which is a specific instance of an exact binomial group sequential design for a single arm trial with a binary outcome.

gsBinomialPP computes a truncated binomial (group) sequential design based on predictive probability.

nBinomial1Sample uses exact binomial calculations to compute power and sample size for single arm binomial experiments.

gsBinomialExact is based on the book "Group Sequential Methods with Applications to Clinical Trials," Christopher Jennison and Bruce W. Turnbull, Chapter 12, Section 12.1.2 Exact Calculations for Binary Data. This computation is often used as an approximation for the distribution of the number of events in one treatment group out of all events when the probability of an event is small and sample size is large.

An object of class gsBinomialExact is returned. On output, the values of theta input to gsBinomialExact will be the parameter values for which the boundary crossing probabilities and expected sample sizes are computed.

Note that a[1] equal to -1 lower bound at n.I[1] means 0 successes continues at interim 1; a[2]==0 at interim 2 means 0 successes stops trial for futility at 2nd analysis. For final analysis, set a[k] equal to b[k]-1 to incorporate all possibilities into non-positive trial; see example.

The sequential probability ratio test (SPRT) is a sequential testing scheme allowing testing after each observation. This likelihood ratio is used to determine upper and lower cutoffs which are linear and parallel in the number of responses as a function of sample size. binomialSPRT produces a variation the the SPRT that tests only within a range of sample sizes. While the linear SPRT bounds are continuous, actual bounds are the integer number of response at or beyond each linear bound for each sample size where testing is performed. Because of the truncation and discretization of the bounds, power and Type I error achieve will be lower than the nominal levels specified by alpha and beta which can be altered to produce desired values that are achieved by the planned sample size. See also example that shows computation of Type I error when futility bound is considered non-binding.

Note that if the objective of a design is to demonstrate that a rate (e.g., failure rate) is lower than a certain level, two approaches can be taken. First, 1 minus the failure rate is the success rate and this can be used for planning. Second, the role of beta becomes to express Type I error and alpha is used to express Type II error.

Plots produced include boundary plots, expected sample size, response rate at the boundary and power.

gsBinomial1Sample uses exact binomial computations based on the base R functions qbinom() and pbinom(). The tabular output may be convenient for plotting. Note that input variables are largely not checked, so the user is largely responsible for results; it is a good idea to do a run with outtype=3 to check that you have done things appropriately. If n is not ordered (a bad idea) or not sequential (maybe OK), be aware of possible consequences.

nBinomial1Sample is based on code from Marc Schwartz [email protected]. The possible sample size vector n needs to be selected in such a fashion that it covers the possible range of values that include the true minimum. NOTE: the one-sided evaluation of significance is more conservative than using the 2-sided exact test in binom.test.

Usage

gsBinomialExact(
  k = 2,
  theta = c(0.1, 0.2),
  n.I = c(50, 100),
  a = c(3, 7),
  b = c(20, 30)
)

binomialSPRT(
  p0 = 0.05,
  p1 = 0.25,
  alpha = 0.1,
  beta = 0.15,
  minn = 10,
  maxn = 35
)

## S3 method for class 'gsBinomialExact'
plot(x, plottype = 1, ...)

## S3 method for class 'binomialSPRT'
plot(x, plottype = 1, ...)

nBinomial1Sample(
  p0 = 0.9,
  p1 = 0.95,
  alpha = 0.025,
  beta = NULL,
  n = 200:250,
  outtype = 1,
  conservative = FALSE
)

Arguments

k	Number of analyses planned, including interim and final. For gsBinomialExact(), k = 1 supports a fixed design.
theta	Vector of possible underling binomial probabilities for a single binomial sample.
n.I	Sample size at analyses (increasing positive integers); vector of length k.
a	Number of "successes" required to cross lower bound cutoffs to reject p1 in favor of p0 at each analysis; vector of length k; -1 (minimum allowed) means no lower bound.
b	Number of "successes" required to cross upper bound cutoffs for rejecting p0 in favor of p1 at each analysis; vector of length k; value > n.I means no upper bound.
p0	Lower of the two response (event) rates hypothesized.
p1	Higher of the two response (event) rates hypothesized.
alpha	Nominal probability of rejecting response (event) rate p0 when it is true.
beta	Nominal probability of rejecting response (event) rate p1 when it is true. If NULL, Type II error will be computed for all input values of n and output will be in a data frame.
minn	Minimum sample size at which sequential testing begins.
maxn	Maximum sample size.
x	Item of class gsBinomialExact or binomialSPRT for print.gsBinomialExact. Item of class gsBinomialExact for plot.gsBinomialExact. Item of class binomialSPRT for item of class plot.binomialSPRT.
plottype	1 produces a plot with counts of response at bounds (for binomialSPRT, also produces linear SPRT bounds); 2 produces a plot with power to reject null and alternate response rates as well as the probability of not crossing a bound by the maximum sample size; 3 produces a plot with the response rate at the boundary as a function of sample size when the boundary is crossed; 6 produces a plot of the expected sample size by the underlying event rate (this assumes there is no enrollment beyond the sample size where the boundary is crossed).
...	arguments passed through to ggplot.
n	sample sizes to be considered for nBinomial1Sample. These should be ordered from smallest to largest and be > 0.
outtype	Operative when beta != NULL. 1 means routine will return a single integer sample size while for output=2a data frame is returned (see Value); note that this not operative is beta is NULL in which case a data table is returned with Type II error and power for each input value of n.
conservative	operative when outtype=1 or 2 and beta != NULL. Default FALSE selects minimum sample size for which power is at least 1-beta. When conservative=TRUE, the minimum sample sample size for which power is at least 1-beta and there is no larger sample size in the input n where power is less than 1-beta.

Value

gsBinomialExact() returns a list of class gsBinomialExact and gsProbability (see example); when displaying one of these objects, the default function to print is print.gsProbability(). The object returned from gsBinomialExact() contains the following elements:

k	As input.
theta	As input.
n.I	As input.
lower	A list containing two elements: bound is as input in a and prob is a matrix of boundary crossing probabilities. Element i,j contains the boundary crossing probability at analysis i for the j-th element of theta input. All boundary crossing is assumed to be binding for this computation; that is, the trial must stop if a boundary is crossed.
upper	A list of the same form as lower containing the upper bound and upper boundary crossing probabilities.
en	A vector of the same length as theta containing expected sample sizes for the trial design corresponding to each value in the vector theta.

binomialSPRT produces an object of class binomialSPRT that is an extension of the gsBinomialExact class. The values returned in addition to those returned by gsBinomialExact are:

intercept	A vector of length 2 with the intercepts for the two SPRT bounds.
slope	A scalar with the common slope of the SPRT bounds.
alpha	As input. Note that this will exceed the actual Type I error achieved by the design returned.
beta	As input. Note that this will exceed the actual Type II error achieved by the design returned.
p0	As input.
p1	As input.

nBinomial1Sample produces a data frame with power for each input value in n if beta=NULL. Otherwise, a sample size achieving the desired power is returned unless the minimum power for the values input in n is greater than or equal to the target or the maximum yields power less than the target, in which case an error message is shown. The input variable outtype has no effect if beta=NULL. Otherwise, outtype=1 results in return of an integer sample size and outtype=2 results in a data frame with one record which includes the desired sample size. When a data frame is returned, the variables include:

p0	Input null hypothesis event (or response) rate.
p1	Input alternative hypothesis (or response) rate; must be > p0.
alpha	Input Type I error.
beta	Input Type II error except when input is NULL in which case realized Type II error is computed.
n	sample size.
b	cutoff given n to control Type I error; value is NULL if no such value exists.
alphaR	Type I error achieved for each output value of n; less than or equal to the input value alpha.
Power	Power achieved for each output value of n.

Note

The gsDesign technical manual is available at https://keaven.github.io/gsd-tech-manual/.

Author(s)

Jon Hartzel, Yevgen Tymofyeyev and Keaven Anderson [email protected]

References

Jennison C and Turnbull BW (2000), Group Sequential Methods with Applications to Clinical Trials. Boca Raton: Chapman and Hall.

Code for nBinomial1Sample was based on code developed by [email protected].

See Also

gsProbability

Examples

library(ggplot2)

zz <- gsBinomialExact(
  k = 3, theta = seq(0.1, 0.9, 0.1), n.I = c(12, 24, 36),
  a = c(-1, 0, 11), b = c(5, 9, 12)
)

# let's see what class this is
class(zz)

# because of "gsProbability" class above, following is equivalent to
# \code{print.gsProbability(zz)}
zz

# also plot (see also plots below for \code{binomialSPRT})
# add lines using geom_line()
plot(zz) + 
ggplot2::geom_line()

# now for SPRT examples
x <- binomialSPRT(p0 = .05, p1 = .25, alpha = .1, beta = .2)
# boundary plot
plot(x)
# power plot
plot(x, plottype = 2)
# Response (event) rate at boundary
plot(x, plottype = 3)
# Expected sample size at boundary crossing or end of trial
plot(x, plottype = 6)

# sample size for single arm exact binomial

# plot of table of power by sample size
# note that outtype need not be specified if beta is NULL
nb1 <- nBinomial1Sample(p0 = 0.05, p1=0.2,alpha = 0.025, beta=NULL, n = 25:40)
nb1
library(scales)
ggplot2::ggplot(nb1, ggplot2::aes(x = n, y = Power)) + 
ggplot2::geom_line() + 
ggplot2::geom_point() + 
ggplot2::scale_y_continuous(labels = percent)

# simple call with same parameters to get minimum sample size yielding desired power
nBinomial1Sample(p0 = 0.05, p1 = 0.2, alpha = 0.025, beta = .2, n = 25:40)

# change to 'conservative' if you want all larger sample
# sizes to also provide adequate power
nBinomial1Sample(p0 = 0.05, p1 = 0.2, alpha = 0.025, beta = .2, n = 25:40, conservative = TRUE)

# print out more information for the selected derived sample size
nBinomial1Sample(p0 = 0.05, p1 = 0.2, alpha = 0.025, beta = .2, n = 25:40, conservative = TRUE,
 outtype = 2)
# Reproduce realized Type I error alphaR
stats::pbinom(q = 5, size = 39, prob = .05, lower.tail = FALSE)
# Reproduce realized power
stats::pbinom(q = 5, size = 39, prob = 0.2, lower.tail = FALSE)
# Reproduce critical value for positive finding
stats::qbinom(p = 1 - .025, size = 39, prob = .05) + 1
# Compute p-value for 7 successes
stats::pbinom(q = 6, size = 39, prob = .05, lower.tail = FALSE)
# what happens if input sample sizes not sufficient?
## Not run:  
  nBinomial1Sample(p0 = 0.05, p1 = 0.2, alpha = 0.025, beta = .2, n = 25:30)

## End(Not run)

---

Boundary derivation - low level

Description

gsBound() and gsBound1() are lower-level functions used to find boundaries for a group sequential design. They are not recommended (especially gsBound1()) for casual users. These functions do not adjust sample size as gsDesign() does to ensure appropriate power for a design.

gsBound() computes upper and lower bounds given boundary crossing probabilities assuming a mean of 0, the usual null hypothesis. gsBound1() computes the upper bound given a lower boundary, upper boundary crossing probabilities and an arbitrary mean (theta).

The function gsBound1() requires special attention to detail and knowledge of behavior when a design corresponding to the input parameters does not exist.

Usage

gsBound(I, trueneg, falsepos, tol = 1e-06, r = 18, printerr = 0)

gsBound1(theta, I, a, probhi, tol = 1e-06, r = 18, printerr = 0)

Arguments

I	Vector containing statistical information planned at each analysis.
trueneg	Vector of desired probabilities for crossing upper bound assuming mean of 0.
falsepos	Vector of desired probabilities for crossing lower bound assuming mean of 0.
tol	Tolerance for error (scalar; default is 0.000001). Normally this will not be changed by the user. This does not translate directly to number of digits of accuracy, so use extra decimal places.
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 $10^{-6}$ with r = 16. This parameter is normally not changed by users.
printerr	If this scalar argument set to 1, this will print messages from underlying C program. Mainly intended to notify user when an output solution does not match input specifications. This is not intended to stop execution as this often occurs when deriving a design in gsDesign that uses beta-spending.
theta	Scalar containing mean (drift) per unit of statistical information.
a	Vector containing lower bound that is fixed for use in gsBound1.
probhi	Vector of desired probabilities for crossing upper bound assuming mean of theta.

Value

Both routines return a list. Common items returned by the two routines are:

k	The length of vectors input; a scalar.
theta	As input in gsBound1(); 0 for gsBound().
I	As input.
a	For gsbound1, this is as input. For gsbound this is the derived lower boundary required to yield the input boundary crossing probabilities under the null hypothesis.
b	The derived upper boundary required to yield the input boundary crossing probabilities under the null hypothesis.
tol	As input.
r	As input.
error	Error code. 0 if no error; greater than 0 otherwise.

gsBound() also returns the following items:

rates	a list containing two items:
falsepos	vector of upper boundary crossing probabilities as input.
trueneg	vector of lower boundary crossing probabilities as input.

gsBound1() also returns the following items:

problo	vector of lower boundary crossing probabilities; computed using input lower bound and derived upper bound.
probhi	vector of upper boundary crossing probabilities as input.

Note

The gsDesign technical manual is available at https://keaven.github.io/gsd-tech-manual/.

Author(s)

Keaven Anderson [email protected]

References

Jennison C and Turnbull BW (2000), Group Sequential Methods with Applications to Clinical Trials. Boca Raton: Chapman and Hall.

See Also

vignette("gsDesignPackageOverview"), gsDesign, gsProbability

Examples

# set boundaries so that probability is .01 of first crossing
# each upper boundary and .02 of crossing each lower boundary
# under the null hypothesis
x <- gsBound(
  I = c(1, 2, 3) / 3, trueneg = rep(.02, 3),
  falsepos = rep(.01, 3)
)
x

#  use gsBound1 to set up boundary for a 1-sided test
x <- gsBound1(
  theta = 0, I = c(1, 2, 3) / 3, a = rep(-Inf, 3),
  probhi = c(.001, .009, .015)
)
x$b

# check boundary crossing probabilities with gsProbability
y <- gsProbability(k = 3, theta = 0, n.I = x$I, a = x$a, b = x$b)$upper$prob

#  Note that gsBound1 only computes upper bound
#  To get a lower bound under a parameter value theta:
#      use minus the upper bound as a lower bound
#      replace theta with -theta
#      set probhi as desired lower boundary crossing probabilities
#  Here we let set lower boundary crossing at 0.05 at each analysis
#  assuming theta=2.2
y <- gsBound1(
  theta = -2.2, I = c(1, 2, 3) / 3, a = -x$b,
  probhi = rep(.05, 3)
)
y$b

#  Now use gsProbability to look at design
#  Note that lower boundary crossing probabilities are as
#  specified for theta=2.2, but for theta=0 the upper boundary
#  crossing probabilities are smaller than originally specified
#  above after first interim analysis
gsProbability(k = length(x$b), theta = c(0, 2.2), n.I = x$I, b = x$b, a = -y$b)

---

Conditional Power at Interim Boundaries

Description

gsBoundCP() computes the total probability of crossing future upper bounds given an interim test statistic at an interim bound. For each interim boundary, assumes an interim test statistic at the boundary and computes the probability of crossing any of the later upper boundaries.

See Conditional power section of manual for further clarification. See also Muller and Schaffer (2001) for background theory.

Usage

gsBoundCP(x, theta = "thetahat", r = 18)

Arguments

x	An object of type gsDesign or gsProbability
theta	if "thetahat" and class(x)!="gsDesign", conditional power computations for each boundary value are computed using estimated treatment effect assuming a test statistic at that boundary (zi/sqrt(x$n.I[i]) at analysis i, interim test statistic zi and interim sample size/statistical information of x$n.I[i]). Otherwise, conditional power is computed assuming the input scalar value theta.
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 $10^{-6}$ with r = 16. This parameter is normally not changed by users.

Value

For one-sided designs, a vector of length x$k-1 containing conditional power at each interim upper bound. For two-sided designs, a matrix with x$k-1 rows and columns CPlo and CPhi, containing conditional power at the lower and upper bounds, respectively. Designs with harm bounds (test.type 7 or 8) also have a CPharm column. An absent bound has conditional power NA.

Note

The gsDesign technical manual is available at https://keaven.github.io/gsd-tech-manual/.

Author(s)

Keaven Anderson [email protected]

References

Jennison C and Turnbull BW (2000), Group Sequential Methods with Applications to Clinical Trials. Boca Raton: Chapman and Hall.

Muller, Hans-Helge and Schaffer, Helmut (2001), Adaptive group sequential designs for clinical trials: combining the advantages of adaptive and classical group sequential approaches. Biometrics;57:886-891.

See Also

gsDesign, gsProbability, gsCP

Examples

# set up a group sequential design
x <- gsDesign(k = 5)
x

# compute conditional power based on interim treatment effects
gsBoundCP(x)

# compute conditional power based on original x$delta
gsBoundCP(x, theta = x$delta)

---

Calibrate Futility Spending at Fixed Information

Description

Compatibility wrapper for gsCPOSFutilitySpending(mode = "fixed_information"), retaining the legacy argument, diagnostic and error names. Fit beta spending parameters to conditional assurance targets while holding the reference information and efficacy boundaries fixed. Overall power may change. Total beta is solved internally so that the spending rule and the terminal decision at the fixed efficacy boundary are consistent.

Usage

gsCAFutilitySpending(
  x,
  target_ca,
  i = seq_along(target_ca),
  sfl = x$lower$sf,
  prior,
  control = list()
)

Arguments

x	A gsDesign, gsSurv, gsSurvCalendar, or gsSurvPower design with test.type 3 or 4. The analysis information fractions in x$timing are held fixed. These fractions are distinct from the times supplied to the spending functions; see Information fractions and spending times below.
target_ca	Conditional assurance targets strictly between zero and one.
i	Unique active interim futility indices, one per target.
sfl	Supported lower spending function or its name. Defaults to x$lower$sf, the reference futility spending function. Supply sfl to override it. See Details.
prior	List with finite numeric vectors z (standardized effects on the gsCPOS() theta scale) and wgts (nonnegative prior masses or density-weighted quadrature weights). Weights must have positive total and are normalized internally. The prior must be supplied explicitly and is held fixed throughout fitting. For the default gsBoundSummary() prior, use normalGrid(mu = x$delta / 2, sigma = 10 / sqrt(x$n.fix)).
control	Named numerical controls. Use ca_tol for the maximum absolute target residual (default 1e-4, finite and in (0, 0.1)). All other controls and defaults are as in gsCPOSFutilitySpending: start, lower, upper, maxit, reltol, backward, and trace. Unknown or invalid controls are errors.

Details

This is the fixed information counterpart of gsCPOSFutilitySpending. The prior, information, efficacy boundaries, spending times and testing indicators are held fixed. For each candidate spending parameter set, gsBound1() derives interim futility bounds, and gsProbability() evaluates total beta when final lower and upper bounds coincide. An internal scalar solve makes that beta agree with the beta used by the spending function.

The probability target conditions on continuation through an analysis, not on an observed statistic at a boundary. An externally calculated fixed design gsPOS() benchmark can be supplied as a target. This implements a spending-based fixed information conditional assurance rule, not a full sample size re-estimation procedure.

With nonbinding futility (test.type = 4), the efficacy-only type I error specification is unchanged. With binding futility (test.type = 3), changing futility while freezing efficacy can change actual type I error, including increasing it above the reference alpha. Inspect achieved_type1; this function does not promise preservation of alpha or power in that case. The reference nominal alpha remains in x$alpha.

Replay by calling this function with the original reference, the fitted spending family, prior and targets, and control$start set to the fitted free parameters. A usual power-preserving gsDesign() call is not an equivalent reconstruction. The internally fitted beta and spending parameters together specify the lower spending rule at the fixed information.

Value

A c("gsCAFutilitySpending", "gsDesign") object, with caFutilitySpending diagnostics analogous to those of gsCPOSFutilitySpending, using target_ca, achieved_ca and ca_tol. Additional fields include achieved_beta, achieved_type1, and fixed_information. The returned beta is achieved overall beta, not the reference beta. Error classes use the prefix gsCAFutilitySpending.

Information fractions and spending times

All six spending calibrators keep the reference analysis information fractions x$timing fixed. These are the cumulative information fractions x$n.I / x$n.I[x$k], ending at 1. For example, x$timing = c(.5, .75, 1) keeps the analyses at 50%, 75%, and 100% of the final information. Power-preserving calibration may change the maximum information and thus the absolute information n.I at every analysis while retaining these fractions. Fixed information calibration (gsCAFutilitySpending() and gsCPOSFutilitySpending(mode = "fixed_information")) also holds n.I and the efficacy boundaries fixed.

Spending times are the inputs to the spending functions, stored in x$upper$sTime and x$lower$sTime. These are also retained during calibration but may differ from the information fractions; for example, calendar-based spending uses fractions of calendar time. Thus fixed information fractions do not mean that spending must use information time, or that absolute calendar analysis dates must be fixed. See Survival designs for how survival calendar times are handled.

Survival designs

Survival inputs retain their survival classes and endpoint assumptions. Candidate probabilities are evaluated on the statistical event count scale. Power-preserving probability calibration of fixed duration, rate-scaled designs rebuilds the survival plan only for the selected fit and checks all targets again on the returned object. If the deferred search or that check fails, calibration retries once with full survival reconstruction, retaining the best available internal parameters as starting values. Effect calibration and accrual or follow-up duration solves reconstruct the plan for every candidate. For gsSurv() and gsSurvCalendar() inputs, information fractions, spending times, and the enrollment/follow-up constraint are retained; enrollment rates or durations are recalculated as required. Calendar designs with fixed enrollment and follow-up retain their calendar schedule up to numerical tolerance. Stored calls are not evaluated.

For gsSurvPower() inputs, power-preserving calibration fixes the realized calendar times and enrollment periods and rescales enrollment rates to attain the fitted event counts. The evaluated alternative x$hr and its achieved power are used, even if the original design alternative x$hr1 differed. Original event trigger and calendar cap rules are not reapplied: the realized schedule becomes the new plan. Fixed information conditional POS/CA calibration instead retains the survival plan, event counts and efficacy bounds while updating futility and achieved power.

Priors and explicit theta remain standardized drifts per square root event, not hazard ratios. Rounding with toInteger() after calibration can change the target; calibration of an already rounded reference may return noninteger event counts. The final analysis is not a valid target index for interim calibration: i identifies the interim bound or continuation event at which the target is evaluated.

See Also

gsCPOSFutilitySpending, gsCPOS, gsPOS, gsBound1

Examples

x <- gsDesign(k = 2, test.type = 4, sflpar = 0)
prior <- list(z = c(0, x$delta), wgts = c(.2, .8))
target <- gsCPOS(1, x, prior$z, prior$wgts)
fit <- gsCAFutilitySpending(x, target, prior = prior,
                           control = list(start = 1))
fit$caFutilitySpending[c("sflpar", "achieved_ca", "achieved_beta")]
stopifnot(identical(fit$n.I, x$n.I),
          identical(fit$upper$bound, x$upper$bound))

---

Conditional and Predictive Power, Overall and Conditional Probability of Success

Description

gsCP() computes conditional boundary crossing probabilities at future planned analyses for a given group sequential design assuming an interim z-statistic at a specified interim analysis. While gsCP() is designed toward computing conditional power for a variety of underlying parameter values, condPower is built to compute conditional power for a variety of interim test statistic values which is useful for sample size adaptation (see ssrCP). gsPP() averages conditional power across a posterior distribution to compute predictive power. gsPI() computes Bayesian prediction intervals for future analyses corresponding to results produced by gsPP(). gsPosterior() computes the posterior density for the group sequential design parameter of interest given a prior density and an interim outcome that is exact or in an interval. gsPOS() computes the probability of success for a trial using a prior distribution to average power over a set of theta values of interest. gsCPOS() assumes no boundary has been crossed before and including an interim analysis of interest, and computes the probability of success based on this event. Note that gsCP() and gsPP() take only the interim test statistic into account in computing conditional probabilities, while gsCPOS() conditions on not crossing any bound through a specified interim analysis.

See Conditional power section of manual for further clarification. See also Muller and Schaffer (2001) for background theory.

For gsPP(), gsPI(), gsPOS() and gsCPOS(), the prior distribution for the standardized parameter theta () for a group sequential design specified through a gsDesign object is specified through the arguments theta and wgts. This can be a discrete or a continuous probability density function. For a discrete function, generally all weights would be 1. For a continuous density, the wgts would contain integration grid weights, such as those provided by normalGrid.

For gsPosterior, a prior distribution in prior must be composed of the vectors z density. The vector z contains points where the prior is evaluated and density the corresponding density or, for a discrete distribution, the probabilities of each point in z. Densities may be supplied as from normalGrid() where grid weights for numerical integration are supplied in gridwgts. If gridwgts are not supplied, they are defaulted to 1 (equal weighting). To ensure a proper prior distribution, you must have sum(gridwgts * density) equal to 1; this is NOT checked, however.

Usage

gsCP(x, theta = NULL, i = 1, zi = 0, r = 18)

gsPP(
  x,
  i = 1,
  zi = 0,
  theta = c(0, 3),
  wgts = c(0.5, 0.5),
  r = 18,
  total = TRUE
)

gsPI(
  x,
  i = 1,
  zi = 0,
  j = 2,
  level = 0.95,
  theta = c(0, 3),
  wgts = c(0.5, 0.5)
)

gsPosterior(x = gsDesign(), i = 1, zi = NULL, prior = normalGrid(), r = 18)

gsPOS(x, theta, wgts)

gsCPOS(i, x, theta, wgts)

Arguments

x	An object of type gsDesign or gsProbability
theta	a vector with $\theta$ value(s) at which conditional power is to be computed; for gsCP() if NULL, an estimated value of $\theta$ based on the interim test statistic (zi/sqrt(x$n.I[i])) as well as at x$theta is computed. For gsPosterior, this may be a scalar or an interval; for gsPP and gsCP, this must be a scalar.
i	analysis at which interim z-value is given; must be from 1 to x$k-1
zi	interim z-value at analysis i (scalar)
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 $10^{-6}$ with r = 16. This parameter is normally not changed by users.
wgts	Weights to be used with grid points in theta. Length can be one if weights are equal, otherwise should be the same length as theta. Values should be positive, but do not need to sum to 1.
total	The default of total=TRUE produces the combined probability for all planned analyses after the interim analysis specified in i. Otherwise, information on each analysis is provided separately.
j	specific analysis for which prediction is being made; must be >i and no more than x$k
level	The level to be used for Bayes credible intervals (which approach confidence intervals for vague priors). The default level=.95 corresponds to a 95% credible interval. level=0 provides a point estimate rather than an interval.
prior	provides a prior distribution in the form produced by normalGrid

Value

gsCP() returns an object of the class gsProbability. Based on the input design and the interim test statistic, the output gsDesign object has bounds for test statistics computed based on solely on observations after interim i. Boundary crossing probabilities are computed for the input $\theta$ values. See manual and examples.

gsPP() if total==TRUE, returns a real value indicating the predictive power of the trial conditional on the interim test statistic zi at analysis i; otherwise returns vector with predictive power for each future planned analysis.

gsPI() returns an interval (or point estimate if level=0) indicating 100level% credible interval for the z-statistic at analysis j conditional on the z-statistic at analysis i<j. The interval does not consider intervending interim analyses. The probability estimate is based on the predictive distribution used for gsPP() and requires a prior distribution for the group sequential parameter theta specified in theta and wgts.

gsPosterior() returns a posterior distribution containing the the vector z input in prior$z, the posterior density in density, grid weights for integrating the posterior density as input in prior$gridwgts or defaulted to a vector of ones, and the product of the output values in density and gridwgts in wgts.

gsPOS() returns a real value indicating the probability of a positive study weighted by the prior distribution input for theta.

gsCPOS() returns a real value indicating the probability of a positive study weighted by the posterior distribution derived from the interim test statistic and the prior distribution input for theta conditional on an interim test statistic.

Note

The gsDesign technical manual is available at https://keaven.github.io/gsd-tech-manual/.

Author(s)

Keaven Anderson [email protected]

References

Jennison C and Turnbull BW (2000), Group Sequential Methods with Applications to Clinical Trials. Boca Raton: Chapman and Hall.

Proschan, Michael A., Lan, KK Gordon and Wittes, Janet Turk (2006), Statistical Monitoring of Clinical Trials. NY: Springer.

Muller, Hans-Helge and Schaffer, Helmut (2001), Adaptive group sequential designs for clinical trials: combining the advantages of adaptive and classical group sequential approaches. Biometrics;57:886-891.

See Also

normalGrid, gsDesign, gsProbability, gsBoundCP, ssrCP, condPower

Examples

library(ggplot2)
# set up a group sequential design
x <- gsDesign(k = 5)
x

# set up a prior distribution for the treatment effect
# that is normal with mean .75*x$delta and standard deviation x$delta/2
mu0 <- .75 * x$delta
sigma0 <- x$delta / 2
prior <- normalGrid(mu = mu0, sigma = sigma0)

# compute POS for the design given the above prior distribution for theta
gsPOS(x = x, theta = prior$z, wgts = prior$wgts)

# assume POS should only count cases in prior where theta >= x$delta/2
gsPOS(x = x, theta = prior$z, wgts = prior$wgts * (prior$z >= x$delta / 2))

# assuming a z-value at lower bound at analysis 2, what are conditional
# boundary crossing probabilities for future analyses
# assuming theta values from x as well as a value based on the interim
# observed z
CP <- gsCP(x, i = 2, zi = x$lower$bound[2])
CP

# summing values for crossing future upper bounds gives overall
# conditional power for each theta value
CP$theta
t(CP$upper$prob) %*% c(1, 1, 1)

# compute predictive probability based on above assumptions
gsPP(x, i = 2, zi = x$lower$bound[2], theta = prior$z, wgts = prior$wgts)

# if it is known that boundary not crossed at interim 2, use
# gsCPOS to compute conditional POS based on this
gsCPOS(x = x, i = 2, theta = prior$z, wgts = prior$wgts)

# 2-stage example to compare results to direct computation
x <- gsDesign(k = 2)
z1 <- 0.5
n1 <- x$n.I[1]
n2 <- x$n.I[2] - x$n.I[1]
thetahat <- z1 / sqrt(n1)
theta <- c(thetahat, 0, x$delta)

# conditional power direct computation - comparison w gsCP
pnorm((n2 * theta + z1 * sqrt(n1) - x$upper$bound[2] * sqrt(n1 + n2)) / sqrt(n2))

gsCP(x = x, zi = z1, i = 1)$upper$prob

# predictive power direct computation - comparison w gsPP
# use same prior as above
mu0 <- .75 * x$delta * sqrt(x$n.I[2])
sigma2 <- (.5 * x$delta)^2 * x$n.I[2]
prior <- normalGrid(mu = .75 * x$delta, sigma = x$delta / 2)
gsPP(x = x, zi = z1, i = 1, theta = prior$z, wgts = prior$wgts)
t <- .5
z1 <- .5
b <- z1 * sqrt(t)
# direct from Proschan, Lan and Wittes eqn 3.10
# adjusted drift at n.I[2]
pnorm(((b - x$upper$bound[2]) * (1 + t * sigma2) +
  (1 - t) * (mu0 + b * sigma2)) /
  sqrt((1 - t) * (1 + sigma2) * (1 + t * sigma2)))

# plot prior then posterior distribution for unblinded analysis with i=1, zi=1
xp <- gsPosterior(x = x, i = 1, zi = 1, prior = prior)
plot(x = xp$z, y = xp$density, type = "l", col = 2, xlab = expression(theta), ylab = "Density")
points(x = x$z, y = x$density, col = 1)

# add posterior plot assuming only knowlede that interim bound has
# not been crossed at interim 1
xpb <- gsPosterior(x = x, i = 1, zi = 1, prior = prior)
lines(x = xpb$z, y = xpb$density, col = 4)

# prediction interval based in interim 1 results
# start with point estimate, followed by 90% prediction interval
gsPI(x = x, i = 1, zi = z1, j = 2, theta = prior$z, wgts = prior$wgts, level = 0)
gsPI(x = x, i = 1, zi = z1, j = 2, theta = prior$z, wgts = prior$wgts, level = .9)

---

Exact conditional power for a group sequential binomial design

Description

Computes exact conditional probabilities of crossing future boundaries, given the cumulative experimental arm event count at an interim analysis.

Usage

gsCPBinomialExact(
  x,
  i = 1,
  x.i,
  theta = NULL,
  ve = NULL,
  ratio = x$ratio,
  binding = TRUE
)

Arguments

x	An exact binomial spending design returned by [toBinomialExact()].
i	Interim analysis at which conditioning occurs; must be less than 'x$k'.
x.i	Cumulative experimental arm events observed at analysis 'i'.
theta	Conditional probabilities that a future event occurs in the experimental arm. If 'NULL', the observed probability and the values in 'x$theta' are used.
ve	Vaccine or prevention efficacy assumptions. These are converted to conditional event probabilities using 'ratio'. Specify at most one of 'theta' and 've'.
ratio	Experimental-to-control randomization ratio. This defaults to the ratio retained by [toBinomialExact()]. It is required when 've' is supplied and is also used to report efficacy corresponding to 'theta'.
binding	Logical indicating whether future futility or harm boundaries stop the trial. With 'FALSE', conditional efficacy power ignores these future non-binding boundaries. This argument does not restrict the observed result at analysis 'i'.

Details

Conditional on the interim count, future experimental arm event increments are independent binomial random variables. Their distribution is propagated through the remaining integer boundaries. Thus the calculation is exact under the same conditional binomial assumptions used by [gsBinomialExact()].

The calculation conditions only on the statistic at analysis 'i', as [gsCP()] does. It is computed regardless of whether the observed count has crossed a current efficacy, futility, or harm boundary. This gives a hypothetical projection if follow-up were to continue; it does not reverse a stopping decision. The 'binding' argument controls only whether futility or harm boundaries at future analyses are enforced. Thus, the default 'binding = TRUE' includes future non-efficacy stopping, whereas 'binding = FALSE' ignores it.

Value

An object of class 'gsBinomialExactCP'. Components include future absolute event counts, assumed conditional event probabilities and efficacies, per-analysis efficacy and futility crossing probabilities, continuation probabilities, and total conditional power and futility.

See Also

[gsCP()], [toBinomialExact()], [VEtable()]

Examples

design <- gsSurv(
  k = 2, test.type = 4, timing = .5, ratio = 3,
  hr = .3, hr0 = .7
)
exact_design <- toBinomialExact(design)
gsCPBinomialExact(exact_design, i = 1, x.i = 20, ve = c(.5, .7))

---

Calibrate Futility Spending to Conditional Power Targets

Description

gsCPFutilitySpending() selects parameters for a beta spending futility boundary so that conditional power at one or more interim futility bounds matches specified targets. Each candidate design is reconstructed with gsDesign(), which recalculates the statistical information required to retain the reference design's unconditional power.

Usage

gsCPFutilitySpending(
  x,
  target_cp,
  i = seq_along(target_cp),
  sfl = x$lower$sf,
  theta = NULL,
  control = list()
)

Arguments

x	A gsDesign, gsSurv, gsSurvCalendar, or gsSurvPower design with test.type 3, 4, 7, or 8. The analysis information fractions in x$timing are held fixed. These fractions are distinct from the times supplied to the spending functions; see Information fractions and spending times below.
target_cp	Numeric vector of conditional power targets strictly between zero and one, one per interim index in i. Each target is the probability of crossing any future efficacy bound, from analysis i + 1 through the final analysis x$k, conditional on the statistic at analysis i equaling its candidate futility bound. With theta = NULL, the future effect is the observed effect implied by that bound, lower$bound[i] / sqrt(n.I[i]); otherwise theta specifies the future effect. Future efficacy and futility boundaries remain in force. With multiple targets, each is a total future conditional power at its corresponding interim, not a next analysis conditional power.
i	Interim analysis indices corresponding to target_cp. Values must identify active futility bounds, be unique, and be in 1:(x$k - 1). Results are ordered by analysis.
sfl	Futility spending function, supplied as a supported function or its character name. By default, use x$lower$sf, the reference design's futility spending function. Supply sfl to override it.
theta	Optional future effect for conditional power. A scalar is recycled; otherwise its length must equal target_cp. When NULL, the observed effect at each candidate lower bound is used.
control	Optional named list of numerical solver settings. Unspecified components retain their defaults; ordinary use does not require this argument. These settings affect the search, not the effect assumption or CP target. Unknown names and invalid values produce input errors. start Initial spending parameters (default NULL). When omitted, use the reference design's parameters if its spending function and parameter count match; otherwise use the family defaults below. Custom spending functions require an explicit starting vector. lower, upper Search limits for the spending parameters, not Z-boundaries (both default NULL, selecting the family defaults below). Supplied start and limit vectors must be finite numeric vectors matching the number of free parameters. Limits must satisfy lower < upper and contain the starting values. cp_tol Maximum absolute difference between achieved and targeted CP at every selected analysis (default 1e-4). Must be a finite scalar strictly between 0 and 0.1. maxit Positive integer iteration limit for the joint optimizer (default 500), not a global limit on design evaluations or the one-parameter root search. reltol Positive finite scalar controlling internal numerical convergence (default 1e-10). This does not replace the final cp_tol acceptance check. backward Allow a latest-to-earliest coordinate fallback if the initial joint fit misses the targets (default TRUE). trace Display joint optimizer progress (default FALSE). Both backward and trace must be nonmissing scalar logical values.

Details

Conditional power is the probability of crossing any efficacy bound after the selected interim, through the final analysis, conditional on the interim statistic equaling the candidate lower bound. When theta is NULL, the future effect is the observed effect implied by that bound, lower$bound[i] / sqrt(n.I[i]), following the default convention in gsCP(). The calculation conditions on the interim statistic even though it is at a stopping boundary; future futility bounds remain part of the conditional power calculation.

For a three-analysis design, target_cp = .3, i = 1 targets probability 0.3 of efficacy at analysis 2 or 3, conditional on $Z_1 = a_1$, where $a_1$ is the fitted futility bound. With a two-parameter spending family, target_cp = c(.3, .5), i = c(1, 2) additionally targets probability 0.5 of efficacy at analysis 3 conditional on $Z_2 = a_2$. Each target uses its own interim conditioning value and, when theta = NULL, its own bound-implied future effect. These are total future conditional powers, not probabilities of efficacy only at the next analysis or probabilities that are added together across targeted interims.

One-target calibration supports the one-parameter families sfHSD, sfPower, sfExponential, and sfLDOF. Two-target calibration supports sfLogistic, sfBetaDist, sfCauchy, sfNormal, sfExtremeValue, and sfExtremeValue2. sfLinear may be used with any number of targets and is the intended choice for more than two. Its knot times are fixed at the lower spending times for the targeted analyses; the fitted cumulative spending proportions are constrained to be strictly increasing and between zero and one.

With multiple targets, joint constrained optimization first uses the reference or supplied starting parameters. If needed, a latest-to-earliest coordinate solve supplies alternative starting values for another joint fit. A result is returned only when every conditional power residual is within control$cp_tol. Searches may stop early once all residuals are at most min(control$cp_tol / 10, 1e-7); computations retain double precision.

The fitted lower spending parameters depend on the complete design, including efficacy spending. For test.type 7 and 8 they may also depend on harm spending. Changing any of those specifications requires recalibration. toInteger() carries the fitted spending function and parameters forward, but rounding information can change the achieved conditional power and does not trigger recalibration.

Value

A calibrated design inheriting from gsCPFutilitySpending and gsDesign, retaining gsSurv and gsSurvPower when applicable. The cpFutilitySpending component contains targets, achieved conditional powers, effects, fitted spending metadata, information, reference efficacy and harm specifications, and solver diagnostics.

Invalid inputs raise a gsCPFutilitySpending_input_error. A well-formed target that cannot be attained raises a gsCPFutilitySpending_infeasible_error; optimizer failure raises a gsCPFutilitySpending_convergence_error.

Information fractions and spending times

All six spending calibrators keep the reference analysis information fractions x$timing fixed. These are the cumulative information fractions x$n.I / x$n.I[x$k], ending at 1. For example, x$timing = c(.5, .75, 1) keeps the analyses at 50%, 75%, and 100% of the final information. Power-preserving calibration may change the maximum information and thus the absolute information n.I at every analysis while retaining these fractions. Fixed information calibration (gsCAFutilitySpending() and gsCPOSFutilitySpending(mode = "fixed_information")) also holds n.I and the efficacy boundaries fixed.

Spending times are the inputs to the spending functions, stored in x$upper$sTime and x$lower$sTime. These are also retained during calibration but may differ from the information fractions; for example, calendar-based spending uses fractions of calendar time. Thus fixed information fractions do not mean that spending must use information time, or that absolute calendar analysis dates must be fixed. See Survival designs for how survival calendar times are handled.

Spending parameter search defaults

When not supplied in control, parameter limits and fallback starting values are selected by family:

Family	Start	Lower	Upper
sfHSD	-2	-40	40
sfPower	1	1e-4	50
sfExponential	0.5	1e-4	1.5
sfLDOF	1	0.005	20
sfBetaDist	c(1, 1)	c(1e-3, 1e-3)	c(50, 50)
Other supported two-parameter families	c(0, 1)	c(-20, 1e-3)	c(20, 50)
Custom function	Required	-20 per parameter	20 per parameter

For sfLinear, start instead contains one cumulative spending proportion per target, strictly increasing and in (0, 1). When omitted, starting proportions are derived from the reference design's cumulative lower spending divided by beta and adjusted to satisfy these constraints. User-supplied lower and upper are not supported for its constrained parameterization.

Survival designs

Survival inputs retain their survival classes and endpoint assumptions. Candidate probabilities are evaluated on the statistical event count scale. Power-preserving probability calibration of fixed duration, rate-scaled designs rebuilds the survival plan only for the selected fit and checks all targets again on the returned object. If the deferred search or that check fails, calibration retries once with full survival reconstruction, retaining the best available internal parameters as starting values. Effect calibration and accrual or follow-up duration solves reconstruct the plan for every candidate. For gsSurv() and gsSurvCalendar() inputs, information fractions, spending times, and the enrollment/follow-up constraint are retained; enrollment rates or durations are recalculated as required. Calendar designs with fixed enrollment and follow-up retain their calendar schedule up to numerical tolerance. Stored calls are not evaluated.

For gsSurvPower() inputs, power-preserving calibration fixes the realized calendar times and enrollment periods and rescales enrollment rates to attain the fitted event counts. The evaluated alternative x$hr and its achieved power are used, even if the original design alternative x$hr1 differed. Original event trigger and calendar cap rules are not reapplied: the realized schedule becomes the new plan. Fixed information conditional POS/CA calibration instead retains the survival plan, event counts and efficacy bounds while updating futility and achieved power.

Priors and explicit theta remain standardized drifts per square root event, not hazard ratios. Rounding with toInteger() after calibration can change the target; calibration of an already rounded reference may return noninteger event counts. The final analysis is not a valid target index for interim calibration: i identifies the interim bound or continuation event at which the target is evaluated.

See Also

gsDesign, gsSurv, gsCP, sfLinear, toInteger, gsPPFutilitySpending

Examples

# Lan-DeMets O'Brien-Fleming efficacy spending with futility only at IA 1.
x <- gsDesign(
  k = 3, test.type = 4, timing = c(.5, .75),
  sfu = sfLDOF,
  sfl = sfHSD, sflpar = 1,
  testLower = c(TRUE, FALSE, FALSE)
)
target_cp <- .3
# With theta = NULL (the default), CP uses the observed effect implied
# by the interim futility bound.
fit <- gsCPFutilitySpending(x, target_cp = target_cp, i = 1)
fit$cpFutilitySpending[c("target_cp", "achieved_cp", "sflpar")]

# Use the fitted spending parameter in the final gsDesign.
final_design <- gsDesign(
  k = 3, test.type = 4, timing = c(.5, .75),
  sfu = sfLDOF,
  sfl = sfHSD, sflpar = fit$cpFutilitySpending$sflpar,
  testLower = c(TRUE, FALSE, FALSE)
)

# The final design has CP 0.3000 at the first interim futility bound.
gsBoundSummary(
  final_design,
  exclude = "B-value"
)

# Use the same test type, timing, and spending in a survival design.
# Other survival inputs use gsSurv() defaults. The IA 1 futility CP row
# again shows 0.3000, agreeing with target_cp up to numerical tolerance.
surv_design <- gsSurv(
  k = 3, test.type = 4, timing = c(.5, .75),
  sfu = sfLDOF,
  sfl = sfHSD, sflpar = fit$cpFutilitySpending$sflpar,
  testLower = c(TRUE, FALSE, FALSE)
)
gsBoundSummary(surv_design, exclude = "B-value")

# Survival designs can also be calibrated directly.
surv_fit <- gsCPFutilitySpending(surv_design, target_cp = .7, i = 1)
gsBoundSummary(surv_fit, exclude = "B-value")

# Optionally require an absolute CP residual no larger than 0.000001.
fit_tight <- gsCPFutilitySpending(
  x, target_cp = target_cp, i = 1,
  control = list(cp_tol = 1e-6)
)
abs(fit_tight$cpFutilitySpending$achieved_cp - target_cp) <= 1e-6

# Two targets use a two-parameter family and active futility at both interims.
x_two <- gsDesign(
  k = 3, test.type = 4, timing = c(.5, .75),
  sfl = sfLogistic, sflpar = c(0, 1)
)
fit_two <- gsCPFutilitySpending(
  x_two, target_cp = c(.3, .5), i = c(1, 2), theta = NULL
)
# IA 1: efficacy at IA 2 or FA; IA 2: efficacy at FA.
# Independently recompute total future CP using each fitted bound's effect.
achieved <- vapply(1:2, function(j) {
  cp <- gsCP(fit_two, i = j, zi = fit_two$lower$bound[j],
             theta = fit_two$lower$bound[j] / sqrt(fit_two$n.I[j]))
  sum(cp$upper$prob[, 1])
}, numeric(1))
stopifnot(max(abs(achieved - c(.3, .5))) <= 1e-4,
          isTRUE(all.equal(fit_two$timing, x_two$timing)))

---

Calibrate Futility Spending to Conditional Probability of Success

Description

Select beta spending parameters to match gsCPOS() targets at selected interims, either preserving reference power or fixing reference information.

Usage

gsCPOSFutilitySpending(
  x,
  target_cpos,
  i = seq_along(target_cpos),
  sfl = x$lower$sf,
  prior,
  control = list(),
  mode = c("preserve_power", "fixed_information")
)

Arguments

x	A gsDesign, gsSurv, gsSurvCalendar, or gsSurvPower design with test.type 3 or 4. The analysis information fractions in x$timing are held fixed. These fractions are distinct from the times supplied to the spending functions; see Information fractions and spending times below.
target_cpos	Numeric conditional assurance targets strictly between zero and one, one per selected interim.
i	Unique active interim futility indices, defaulting to seq_along(target_cpos). Results are ordered by analysis.
sfl	Supported lower spending function or its name. Defaults to x$lower$sf, the reference futility spending function. Supply sfl to override it. See Details.
prior	List with finite numeric vectors z (standardized effects on the gsCPOS() theta scale) and wgts (nonnegative prior masses or density-weighted quadrature weights). Weights must have positive total and are normalized internally. The prior must be supplied explicitly and is held fixed throughout fitting. For the default gsBoundSummary() prior, use normalGrid(mu = x$delta / 2, sigma = 10 / sqrt(x$n.fix)).
control	Named list of solver controls. cpos_tol is the maximum absolute target residual (default 1e-4, finite and in (0, 0.1)). The remaining controls have the same definitions and defaults as in gsPPFutilitySpending: start, lower, and upper (all NULL); maxit (500); reltol (1e-10); backward (TRUE); and trace (FALSE). Unknown, unnamed, duplicate or invalid controls are errors. pp_tol and cp_tol are not accepted here.
mode	Design constraint: "preserve_power" (default) rebuilds candidates with gsDesign(), allowing maximum information to change. "fixed_information" holds information and efficacy boundaries fixed, solving total beta internally and allowing overall power to change.

Details

Both modes target the existing gsCPOS() calculation. In fixed information mode, gsBound1() derives lower bounds and an internal beta solve makes spending consistent with the final decision. For binding futility, changing lower bounds while fixing efficacy can increase actual type I error above nominal alpha; inspect achieved_type1. Replay this mode with the original reference and fitted parameters as control$start; an ordinary power-preserving gsDesign() call is not equivalent. The legacy gsCAFutilitySpending interface is a compatibility wrapper.

Conditional assurance is the prior-averaged probability of future efficacy rejection given that neither stopping boundary has been crossed through interim i. It includes the entire previous stopping history. Conditioning on continuation reweights the effect distribution; this differs from gsPP(), which conditions on an exact interim statistic, and from unconditional gsPOS(). A point prior gives fixed effect success conditional on continuation, not conditional power at the futility bound.

A fixed design gsPOS() value at a separately chosen feasible sample size can be computed once and supplied as a benchmark target. The benchmark must not be recomputed from candidates. Preserving power while adjusting information is an extension of a fixed information conditional assurance rule, not an implementation of sample size re-estimation. Inspect sample size inflation, overall operating characteristics and effect sizes at all bounds.

The supported spending families and requirement for one target per free parameter are the same as for gsPPFutilitySpending, including two-parameter families and sfLinear. The shared solver accepts only fits meeting every target tolerance; failure is not a proof of global infeasibility.

Candidate continuation probabilities at or below sqrt(.Machine$double.eps) are rejected to avoid unstable conditioning. Harm bound designs are unsupported because gsCPOS() does not account for their harm stopping probability in its denominator. Replay parameters with the complete reference design and the same prior. Spending functions retain their usual timing flexibility, but changed timing, testing indicators or rounding need not retain exact target values.

Value

A c("gsCPOSFutilitySpending", "gsDesign") object, retaining survival classes when applicable. Component cposFutilitySpending contains targets (target_cpos), achieved values (achieved_cpos), residuals, indices, normalized prior, fitted parameters, information and power, reference settings and solver diagnostics. It also records continuation_probability, joint_future_efficacy and unconditional_pos. Errors inherit from gsCPOSFutilitySpending_error, with suffixes _input_error, _infeasible_error or _convergence_error. Both modes also record mode, fixed_information, achieved_beta and achieved_type1.

Spending parameter search defaults

When not supplied in control, parameter limits and fallback starting values are selected by family:

Family	Start	Lower	Upper
sfHSD	-2	-40	40
sfPower	1	1e-4	50
sfExponential	0.5	1e-4	1.5
sfLDOF	1	0.005	20
sfBetaDist	c(1, 1)	c(1e-3, 1e-3)	c(50, 50)
Other supported two-parameter families	c(0, 1)	c(-20, 1e-3)	c(20, 50)
Custom function	Required	-20 per parameter	20 per parameter

For sfLinear, start instead contains one cumulative spending proportion per target, strictly increasing and in (0, 1). When omitted, starting proportions are derived from the reference design's cumulative lower spending divided by beta and adjusted to satisfy these constraints. User-supplied lower and upper are not supported for its constrained parameterization.

Information fractions and spending times

All six spending calibrators keep the reference analysis information fractions x$timing fixed. These are the cumulative information fractions x$n.I / x$n.I[x$k], ending at 1. For example, x$timing = c(.5, .75, 1) keeps the analyses at 50%, 75%, and 100% of the final information. Power-preserving calibration may change the maximum information and thus the absolute information n.I at every analysis while retaining these fractions. Fixed information calibration (gsCAFutilitySpending() and gsCPOSFutilitySpending(mode = "fixed_information")) also holds n.I and the efficacy boundaries fixed.

Spending times are the inputs to the spending functions, stored in x$upper$sTime and x$lower$sTime. These are also retained during calibration but may differ from the information fractions; for example, calendar-based spending uses fractions of calendar time. Thus fixed information fractions do not mean that spending must use information time, or that absolute calendar analysis dates must be fixed. See Survival designs for how survival calendar times are handled.

Survival designs

Survival inputs retain their survival classes and endpoint assumptions. Candidate probabilities are evaluated on the statistical event count scale. Power-preserving probability calibration of fixed duration, rate-scaled designs rebuilds the survival plan only for the selected fit and checks all targets again on the returned object. If the deferred search or that check fails, calibration retries once with full survival reconstruction, retaining the best available internal parameters as starting values. Effect calibration and accrual or follow-up duration solves reconstruct the plan for every candidate. For gsSurv() and gsSurvCalendar() inputs, information fractions, spending times, and the enrollment/follow-up constraint are retained; enrollment rates or durations are recalculated as required. Calendar designs with fixed enrollment and follow-up retain their calendar schedule up to numerical tolerance. Stored calls are not evaluated.

For gsSurvPower() inputs, power-preserving calibration fixes the realized calendar times and enrollment periods and rescales enrollment rates to attain the fitted event counts. The evaluated alternative x$hr and its achieved power are used, even if the original design alternative x$hr1 differed. Original event trigger and calendar cap rules are not reapplied: the realized schedule becomes the new plan. Fixed information conditional POS/CA calibration instead retains the survival plan, event counts and efficacy bounds while updating futility and achieved power.

Priors and explicit theta remain standardized drifts per square root event, not hazard ratios. Rounding with toInteger() after calibration can change the target; calibration of an already rounded reference may return noninteger event counts. The final analysis is not a valid target index for interim calibration: i identifies the interim bound or continuation event at which the target is evaluated.

See Also

gsCPOS, gsPOS, gsPPFutilitySpending, gsCPFutilitySpending

Examples

x <- gsDesign(k = 3, test.type = 4, timing = c(.5, .75), sflpar = 1)
prior <- list(z = c(0, x$delta), wgts = c(.2, .8))
target <- gsCPOS(i = 1, x = x, theta = prior$z, wgts = prior$wgts)
fit <- gsCPOSFutilitySpending(
  x, target_cpos = target, i = 1, prior = prior,
  control = list(start = 0)
)
fit$cposFutilitySpending$sflpar
gsCPOS(i = 1, x = fit, theta = prior$z, wgts = prior$wgts)
fixed <- gsCPOSFutilitySpending(
  x, target_cpos = target, prior = prior, mode = "fixed_information"
)
fixed$cposFutilitySpending[c("mode", "achieved_cpos", "achieved_beta")]
stopifnot(identical(fixed$n.I, x$n.I))

---

Group sequential design interim density function

Description

Given an interim analysis i of a group sequential design and a vector of real values zi, gsDensity() computes an interim density function at analysis i at the values in zi. For each value in zi, this interim density is the derivative of the probability that the group sequential trial does not cross a boundary prior to the i-th analysis and at the i-th analysis the interim Z-statistic is less than that value. When integrated over the real line, this density computes the probability of not crossing a bound at a previous analysis. It corresponds to the subdistribution function at analysis i that excludes the probability of crossing a bound at an earlier analysis.

The initial purpose of this routine was as a component needed to compute the predictive power for a trial given an interim result; see gsPP.

See Jennison and Turnbull (2000) for details on how these computations are performed.

Usage

gsDensity(x, theta = 0, i = 1, zi = 0, r = 18)

Arguments

x	An object of type gsDesign or gsProbability
theta	a vector with $\theta$ value(s) at which the interim density function is to be computed.
i	analysis at which interim z-values are given; must be from 1 to x$k
zi	interim z-value at analysis i (scalar)
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 $10^{-6}$ with r = 16. This parameter is normally not changed by users.

Value

zi	The input vector zi.
theta	The input vector theta.
density	A matrix with length(zi) rows and length(theta) columns. The subdensity function for z[j], theta[m] at analysis i is returned in density[j,m].

Note

The gsDesign technical manual is available at https://keaven.github.io/gsd-tech-manual/.

Author(s)

Keaven Anderson [email protected]

References

Jennison C and Turnbull BW (2000), Group Sequential Methods with Applications to Clinical Trials. Boca Raton: Chapman and Hall.

See Also

gsDesign, gsProbability, gsBoundCP

Examples

library(ggplot2)
# set up a group sequential design
x <- gsDesign()

# set theta values where density is to be evaluated
theta <- x$theta[2] * c(0, .5, 1, 1.5)

# set zi values from -1 to 7 where density is to be evaluated
zi <- seq(-3, 7, .05)

# compute subdensity values at analysis 2
y <- gsDensity(x, theta = theta, i = 2, zi = zi)

# plot sub-density function for each theta value
plot(y$zi, y$density[, 3],
  type = "l", xlab = "Z",
  ylab = "Interim 2 density", lty = 3, lwd = 2
)
lines(y$zi, y$density[, 2], lty = 2, lwd = 2)
lines(y$zi, y$density[, 1], lwd = 2)
lines(y$zi, y$density[, 4], lty = 4, lwd = 2)
title("Sub-density functions at interim analysis 2")
legend(
  x = c(3.85, 7.2), y = c(.27, .385), lty = 1:5, lwd = 2, cex = 1.5,
  legend = c(
    expression(paste(theta, "=0.0")),
    expression(paste(theta, "=0.5", delta)),
    expression(paste(theta, "=1.0", delta)),
    expression(paste(theta, "=1.5", delta))
  )
)

# add vertical lines with lower and upper bounds at analysis 2
# to demonstrate how likely it is to continue, stop for futility
# or stop for efficacy at analysis 2 by treatment effect
lines(rep(x$upper$bound[2], 2), c(0, .4), col = 2)
lines(rep(x$lower$bound[2], 2), c(0, .4), lty = 2, col = 2)

# Replicate part of figures 8.1 and 8.2 of Jennison and Turnbull text book
# O'Brien-Fleming design with four analyses

x <- gsDesign(k = 4, test.type = 2, sfu = "OF", alpha = .1, beta = .2)

z <- seq(-4.2, 4.2, .05)
d <- gsDensity(x = x, theta = x$theta, i = 4, zi = z)

plot(z, d$density[, 1], type = "l", lwd = 2, ylab = expression(paste(p[4], "(z,", theta, ")")))
lines(z, d$density[, 2], lty = 2, lwd = 2)
u <- x$upper$bound[4]
text(expression(paste(theta, "=", delta)), x = 2.2, y = .2, cex = 1.5)
text(expression(paste(theta, "=0")), x = .55, y = .4, cex = 1.5)

---

Design Derivation

Description

gsDesign() is used to find boundaries and trial size required for a group sequential design.

Many parameters normally take on default values and thus do not require explicit specification. One- and two-sided designs are supported. Two-sided designs may be symmetric or asymmetric. Wang-Tsiatis designs, including O'Brien-Fleming and Pocock designs can be generated. Designs with common spending functions as well as other built-in and user-specified functions for Type I error and futility are supported. Type I error computations for asymmetric designs may assume binding or non-binding lower bounds. The print function has been extended using print.gsDesign() to print gsDesign objects; see examples.

The user may ignore the structure of the value returned by gsDesign() if the standard printing and plotting suffice; see examples.

delta and n.fix are used together to determine what sample size output options the user seeks. The default, delta=0 and n.fix=1, results in a ‘generic’ design that may be used with any sampling situation. Sample size ratios are provided and the user multiplies these times the sample size for a fixed design to obtain the corresponding group sequential analysis times. If delta>0, n.fix is ignored, and delta is taken as the standardized effect size - the signal to noise ratio for a single observation; for example, the mean divided by the standard deviation for a one-sample normal problem. In this case, the sample size at each analysis is computed. When delta=0 and n.fix>1, n.fix is assumed to be the sample size for a fixed design with no interim analyses. See examples below.

Following are further comments on the input argument test.type which is used to control what type of error measurements are used in trial design. The manual may also be worth some review in order to see actual formulas for boundary crossing probabilities for the various options. Options 3 and 5 assume the trial stops if the lower bound is crossed for Type I and Type II error computation (binding lower bound). For the purpose of computing Type I error, options 4 and 6 assume the trial continues if the lower bound is crossed (non-binding lower bound); that is a Type I error can be made by crossing an upper bound after crossing a previous lower bound. Beta-spending refers to error spending for the lower bound crossing probabilities under the alternative hypothesis (options 3 and 4). In this case, the final analysis lower and upper boundaries are assumed to be the same. The appropriate total beta spending (power) is determined by adjusting the maximum sample size through an iterative process for all options. Since options 3 and 4 must compute boundary crossing probabilities under both the null and alternative hypotheses, deriving these designs can take longer than other options. Options 5 and 6 compute lower bound spending under the null hypothesis. Options 7 and 8 add a harm bound to beta spending futility: option 7 treats both futility and harm bounds as binding for Type I error, while option 8 treats both as non-binding. Futility and harm are mutually exclusive lower tail stopping outcomes with separately specified spending functions.

Usage

gsDesign(
  k = 3,
  test.type = 4,
  alpha = 0.025,
  beta = 0.1,
  astar = 0,
  delta = 0,
  n.fix = 1,
  timing = 1,
  sfu = sfHSD,
  sfupar = -4,
  sfl = sfHSD,
  sflpar = -2,
  sfharm = sfLDPocock,
  sfharmparam = -2,
  tol = 1e-06,
  r = 18,
  n.I = 0,
  maxn.IPlan = 0,
  nFixSurv = 0,
  endpoint = NULL,
  delta1 = 1,
  delta0 = 0,
  overrun = 0,
  usTime = NULL,
  lsTime = NULL,
  testUpper = TRUE,
  testLower = TRUE,
  testHarm = TRUE
)

## S3 method for class 'gsDesign'
xtable(
  x,
  caption = NULL,
  label = NULL,
  align = NULL,
  digits = NULL,
  display = NULL,
  ...
)

Arguments

k	Number of analyses planned, including interim and final.
test.type	1=one-sided 2=two-sided symmetric 3=two-sided, asymmetric, beta-spending with binding lower bound 4=two-sided, asymmetric, beta-spending with non-binding lower bound 5=two-sided, asymmetric, lower bound spending under the null hypothesis with binding lower bound 6=two-sided, asymmetric, lower bound spending under the null hypothesis with non-binding lower bound 7=two-sided, asymmetric, with binding futility and binding harm bounds 8=two-sided, asymmetric, with non-binding futility and non-binding harm bounds. See details, examples and manual.
alpha	Type I error, always one-sided. Default value is 0.025.
beta	Type II error, default value is 0.1 (90% power).
astar	Total spending for the lower (test.type 5 or 6) or harm (test.type 7 or 8) bound under the null hypothesis. Default is 0. For test.type 5 or 6, astar specifies the total probability of crossing a lower bound at all analyses combined. For test.type 7 or 8, astar specifies total harm spending under the null, calibrated ignoring futility stopping but retaining efficacy stopping. The default input astar = 0 selects 0.1 for harm designs and $1 -$alpha for types 5 and 6. Actual harm stopping probabilities account for futility and can be smaller than the spending targets. Capping harm at an active futility bound can also reduce attained spending. See the harm monitoring section in gsDesign.
delta	Effect size for theta under alternative hypothesis. This can be set to the standardized effect size to generate a sample size if n.fix=NULL. See details and examples.
n.fix	Sample size for fixed design with no interim; used to find maximum group sequential sample size. For a time-to-event outcome, input number of events required for a fixed design rather than sample size and enter fixed design sample size (optional) in nFixSurv. See details and examples.
timing	Sets relative timing of interim analyses. Default of 1 produces equally spaced analyses. Otherwise, this is a vector of length k or k-1. The values should satisfy 0 < timing[1] < timing[2] < ... < timing[k-1] < timing[k]=1.
sfu	A spending function or a character string indicating a boundary type (that is, “WT” for Wang-Tsiatis bounds, “OF” for O'Brien-Fleming bounds and “Pocock” for Pocock bounds). For one-sided and symmetric two-sided testing is used to completely specify spending (test.type=1, 2), sfu. The default value is sfHSD which is a Hwang-Shih-DeCani spending function. See details, vignette("SpendingFunctionOverview"), manual and examples.
sfupar	Real value, default is $-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 sfupar specifies any parameters needed for the spending function specified by sfu; this is not needed for spending functions (sfLDOF, sfLDPocock) or bound types (“OF”, “Pocock”) that do not require parameters. Note that sfupar can be specified as a positive scalar for sfLDOF for a generalized O'Brien-Fleming spending function.
sfl	Specifies the spending function for lower boundary crossing probabilities when asymmetric, two-sided testing is performed (test.type = 3, 4, 5, or 6). Unlike the upper bound, only spending functions are used to specify the lower bound. The default value is sfHSD which is a Hwang-Shih-DeCani spending function. The parameter sfl is ignored for one-sided testing (test.type=1) or symmetric 2-sided testing (test.type=2). See details, spending functions, manual and examples.
sflpar	Real value, default is $-2$, which, with the default Hwang-Shih-DeCani spending function, specifies a less conservative spending rate than the default for the upper bound.
sfharm	A spending function for the harm bound, used with test.type = 7 or test.type = 8. Default is sfLDPocock. See spendingFunction for details.
sfharmparam	Real value, default is $-2$. Parameter for the harm bound spending function sfharm; ignored by sfLDPocock.
tol	Tolerance for error (default is 0.000001). Normally this will not be changed by the user. This does not translate directly to number of digits of accuracy, so use extra decimal places.
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 $10^{-6}$ with r = 16. This parameter is normally not changed by users. With options(gsDesign.quadrature = "gl") (see the section on numerical integration below) r scales the number of Gauss-Legendre nodes relative to its default.
n.I	Used for re-setting bounds when timing of analyses changes from initial design; see examples.
maxn.IPlan	Used for re-setting bounds when timing of analyses changes from initial design; see examples.
nFixSurv	If a time-to-event variable is used, nFixSurv computed as the sample size from nSurvival may be entered to have gsDesign compute the total sample size required as well as the number of events at each analysis that will be returned in n.fix; this is rounded up to an even number.
endpoint	An optional character string that should represent the type of endpoint used for the study. This may be used by output functions. Types most likely to be recognized initially are "TTE" for time-to-event outcomes with fixed design sample size generated by nSurvival() and "Binomial" for 2-sample binomial outcomes with fixed design sample size generated by nBinomial().
delta1	delta1 and delta0 may be used to store information about the natural parameter scale compared to delta that is a standardized effect size. delta1 is the alternative hypothesis parameter value on the natural parameter scale (e.g., the difference in two binomial rates).
delta0	delta0 is the null hypothesis parameter value on the natural parameter scale.
overrun	Scalar or vector of length k-1 with patients enrolled that are not included in each interim analysis.
usTime	Default is NULL in which case upper bound spending time is determined by timing. Otherwise, this should be a vector of length k with the spending time at each analysis (see Details section of gsDesign).
lsTime	Default is NULL in which case lower bound spending time is determined by timing. Otherwise, this should be a vector of length k with the spending time at each analysis (see Details section of gsDesign).
testUpper	Indicator of which analyses should include an upper (efficacy) bound. A single value of TRUE (default) indicates all analyses have an efficacy bound. Otherwise, a logical vector of length k indicating which analyses will have an efficacy bound. Overridden to all TRUE for test.type 1 and 2. Must be TRUE at the final analysis to achieve targeted power. At each analysis, at least one of testUpper, testLower, or testHarm must be TRUE. Where testUpper is FALSE, the upper bound is set to +Inf and displayed as NA in output.
testLower	Indicator of which analyses should include a lower (futility) bound. A single value of TRUE (default) indicates all analyses have a lower bound; FALSE indicates none. Otherwise, a logical vector of length k. Ignored for test.type 1 (one-sided, no lower bound). Overridden to all TRUE for test.type 2 (symmetric). For test.type 3–8, at least one analysis must be TRUE. Where testLower is FALSE, the lower bound is set to -Inf and displayed as NA in output.
testHarm	Indicator of which analyses should include a harm bound. A single value of TRUE (default) indicates all analyses have a harm bound; FALSE indicates none. Otherwise, a logical vector of length k. Only used for test.type 7 or 8; at least one analysis must be TRUE for those types. Where testHarm is FALSE, the harm bound is set to -Inf and the bound is displayed as NA in output. Cumulative harm crossing probability from earlier analyses is still displayed.
x	An R object of class found among methods(xtable). See below on how to write additional method functions for xtable.
caption	Character vector of length 1 or 2 containing the table's caption or title. If length is 2, the second item is the "short caption" used when LaTeX generates a "List of Tables". Set to NULL to suppress the caption. Default value is NULL.
label	Character vector of length 1 containing the LaTeX label or HTML anchor. Set to NULL to suppress the label. Default value is NULL.
align	Character vector of length equal to the number of columns of the resulting table, indicating the alignment of the corresponding columns. Also, "|" may be used to produce vertical lines between columns in LaTeX tables, but these are effectively ignored when considering the required length of the supplied vector. If a character vector of length one is supplied, it is split as strsplit(align, "")[[1]] before processing. Since the row names are printed in the first column, the length of align is one greater than ncol(x) if x is a data.frame. Use "l", "r", and "c" to denote left, right, and center alignment, respectively. Use "p{3cm}" etc. for a LaTeX column of the specified width. For HTML output the "p" alignment is interpreted as "l", ignoring the width request. Default depends on the class of x.
digits	Numeric vector of length equal to one (in which case it will be replicated as necessary) or to the number of columns of the resulting table or matrix of the same size as the resulting table, indicating the number of digits to display in the corresponding columns. Since the row names are printed in the first column, the length of the vector digits or the number of columns of the matrix digits is one greater than ncol(x) if x is a data.frame. Default depends on the class of x. If values of digits are negative, the corresponding values of x are displayed in scientific format with abs(digits) digits.
display	Character vector of length equal to the number of columns of the resulting table, indicating the format for the corresponding columns. Since the row names are printed in the first column, the length of display is one greater than ncol(x) if x is a data.frame. These values are passed to the formatC function. Use "d" (for integers), "f", "e", "E", "g", "G", "fg" (for reals), or "s" (for strings). "f" gives numbers in the usual xxx.xxx format; "e" and "E" give n.ddde+nn or n.dddE+nn (scientific format); "g" and "G" put x[i] into scientific format only if it saves space to do so. "fg" uses fixed format as "f", but digits as number of significant digits. Note that this can lead to quite long result strings. Default depends on the class of x.
...	Additional arguments. (Currently ignored.)

Value

An object of the class gsDesign. This class has the following elements and upon return from gsDesign() contains:

k	As input.
test.type	As input.
alpha	As input.
beta	As input.
astar	As input, except when test.type=5 or 6 and astar is input as 0; in this case astar is changed to 1-alpha; for types 7 and 8, input 0 is changed to 0.1.
delta	The standardized effect size for which the design is powered. Will be as input to gsDesign() unless it was input as 0; in that case, value will be computed to give desired power for fixed design with input sample size n.fix.
n.fix	Sample size required to obtain desired power when effect size is delta.
timing	A vector of length k containing the portion of the total planned information or sample size at each analysis.
tol	As input.
r	As input.
n.I	Vector of length k. If values are input, same values are output. Otherwise, n.I will contain the sample size required at each analysis to achieve desired timing and beta for the output value of delta. If delta=0 was input, then this is the sample size required for the specified group sequential design when a fixed design requires a sample size of n.fix. If delta=0 and n.fix=1 then this is the relative sample size compared to a fixed design; see details and examples.
maxn.IPlan	As input.
nFixSurv	As input.
nSurv	Sample size for Lachin and Foulkes method when nSurvival is used for fixed design input. If nSurvival is used to compute n.fix, then nFixSurv is inflated by the same amount as n.fix and stored in nSurv. Note that if you use gsSurv for time-to-event sample size, this is not needed and a more complete output summary is given.
endpoint	As input.
delta1	As input.
delta0	As input.
overrun	As input.
usTime	As input.
lsTime	As input.
testUpper	Logical vector of length k indicating which analyses have an efficacy (upper) bound.
testLower	Logical vector of length k indicating which analyses have a futility (lower) bound.
testHarm	Logical vector of length k indicating which analyses have a harm bound (only for test.type 7 or 8).
upper	Upper bound spending function, boundary and boundary crossing probabilities under the NULL and alternate hypotheses. See vignette("SpendingFunctionOverview") and manual for further details.
lower	Lower bound spending function, boundary and boundary crossing probabilities at each analysis. Lower spending is under alternative hypothesis (beta spending) for test.type=3 or 4. For test.type=2, 5 or 6, lower spending is under the null hypothesis. For test.type=1, output value is NULL. See vignette("SpendingFunctionOverview") and manual. For test.type=7 or 8, crossing probabilities exclude harm crossings.
harm	Harm bound spending function, boundary, and crossing probabilities at each analysis for test.type=7 or 8. Harm, lower, and upper crossing probabilities are mutually exclusive and include futility stopping, unlike the harm spending calibration.
theta	Standarized effect size under null (0) and alternate hypothesis. If delta is input, theta[1]=delta. If n.fix is input, theta[1] is computed using a standard sample size formula (pseudocode): ((Zalpha+Zbeta)/theta[1])^2=n.fix.
falseprobnb	For test.type=4 or 6, this contains false positive probabilities under the null hypothesis assuming that crossing a futility bound does not stop the trial.
en	Expected sample size accounting for early stopping. For time-to-event outcomes, this would be the expected number of events (although gsSurv will give expected sample size). For information-based-design, this would give the expected information when the trial stops. If overrun is specified, the expected sample size includes the overrun at each interim.

An object of class "xtable" with attributes specifying formatting options for a table

Numerical integration

Boundary crossing probabilities are computed by the recursive numerical integration of Jennison and Turnbull (2000, Chapter 19). By default the integration grid of that chapter is used (Simpson's rule on 12 r - 3 points concentrated around the mean of the test statistic), which gives probabilities accurate to about $10^{-7}$ with the default r = 18. Setting options(gsDesign.quadrature = "gl") switches all computations (gsDesign(), gsProbability(), gsBound(), gsBound1(), gsDensity() and everything built on them) to Gauss-Legendre quadrature on the continuation region of each analysis, with a number of nodes that adapts to the width of the region relative to the neighboring information increments. In benchmarked designs this is several times faster and accurate to about $10^{-12}$, with differences from the default grid within the accuracy of that grid. The rule uses at most 992 nodes, so extremely small information increments can reduce accuracy. normalGrid() always returns the Jennison and Turnbull grid.

Harm monitoring

For types 7 and 8, harm spending is calibrated under the null as if futility stopping were ignored. Earlier harm and efficacy stops are retained. This convention is separate from binding versus non-binding efficacy calibration: type 7 protects alpha assuming both lower stopping rules are followed, whereas type 8 protects alpha even if they are ignored. Reported harm$prob, power, expected information, and gsBoundSummary() crossing probabilities assume all active stopping rules are followed. Thus actual harm probabilities may be smaller than cumsum(harm$spend). Where both lower bounds are active, harm is capped at futility to retain their ordering; this can also prevent full spending.

The defaults astar = 0.1 (selected by input 0) and sfharm = sfLDPocock provide relatively early spending and aim to flag moderately small one-sided p-values favoring harm. At harm boundary $h$, that nominal p-value is $\Phi(h)$, not the cumulative spending and not the efficacy direction p-value displayed by gsBoundSummary(). There is no fixed nominal p-value cutoff guaranteed at every analysis. With more analyses, a larger total harm spending astar may be desired to retain similarly permissive nominal harm thresholds, at the cost of more false harm signals under the null. This changes harm spending, not futility beta spending. Inspect the resulting nominal thresholds for the actual analysis schedule. See vignette("HarmBound") for examples.

Note

The gsDesign technical manual is available at https://keaven.github.io/gsd-tech-manual/.

Author(s)

Keaven Anderson [email protected]

References

Jennison C and Turnbull BW (2000), Group Sequential Methods with Applications to Clinical Trials. Boca Raton: Chapman and Hall. Lan KK, DeMets DL (1989). Group sequential procedures: calendar versus information time. Statistics in medicine 8(10):1191-8. Liu, Q, Lim, P, Nuamah, I, and Li, Y (2012), On adaptive error spending approach for group sequential trials with random information levels. Journal of biopharmaceutical statistics; 22(4), 687-699.

See Also

vignette("gsDesignPackageOverview"), gsBoundSummary, plot.gsDesign, gsProbability, vignette("SpendingFunctionOverview"),

Normal xtable

Examples

library(ggplot2)
#  symmetric, 2-sided design with O'Brien-Fleming-like boundaries
#  lower bound is non-binding (ignored in Type I error computation)
#  sample size is computed based on a fixed design requiring n=800
x <- gsDesign(k = 5, test.type = 2, n.fix = 800)

# note that "x" below is equivalent to print(x) and print.gsDesign(x)
x
plot(x)
plot(x, plottype = 2)

# Assuming after trial was designed actual analyses occurred after
# 300, 600, and 860 patients, reset bounds
y <- gsDesign(
  k = 3, test.type = 2, n.fix = 800, n.I = c(300, 600, 860),
  maxn.IPlan = x$n.I[x$k]
)
y

#  asymmetric design with user-specified spending that is non-binding
#  sample size is computed relative to a fixed design with n=1000
sfup <- c(.033333, .063367, .1)
sflp <- c(.25, .5, .75)
timing <- c(.1, .4, .7)
x <- gsDesign(
  k = 4, timing = timing, sfu = sfPoints, sfupar = sfup, sfl = sfPoints,
  sflpar = sflp, n.fix = 1000
)
x
plot(x)
plot(x, plottype = 2)

# same design, but with relative sample sizes
gsDesign(
  k = 4, timing = timing, sfu = sfPoints, sfupar = sfup, sfl = sfPoints,
  sflpar = sflp
)

---