Derive spending bound for group sequential boundary
Source:R/gs_spending_bound.R
gs_spending_bound.RdComputes one bound at a time based on spending under given distributional
assumptions. While user specifies gs_spending_bound() for use with other
functions, it is not intended for use on its own.
Most important user specifications are made through a list provided to
functions using gs_spending_bound().
Function uses numerical integration and Newton-Raphson iteration to derive
an individual bound for a group sequential design that satisfies a
targeted boundary crossing probability. Algorithm is a simple extension of
that in Chapter 19 of Jennison and Turnbull (2000).
Arguments
- k
Analysis for which bound is to be computed.
- par
A list with the following items:
sf(class spending function).total_spend(total spend).param(any parameters needed by the spending functionsf()).timing(a vector containing values at which spending function is to be evaluated orNULLif information-based spending is used).max_info(whentimingisNULL, this can be input as positive number to be used withinfofor information fraction at each analysis).
- hgm1
Subdensity grid from
h1()(k=2) orhupdate()(k>2) for analysis k-1; if k=1, this is not used and may beNULL.- theta
Natural parameter used for lower bound only spending; represents average drift at each time of analysis at least up to analysis k; upper bound spending is always set under null hypothesis (theta = 0).
- info
Statistical information at all analyses, at least up to analysis k.
- efficacy
TRUE(default) for efficacy bound,FALSEotherwise.- test_bound
A logical vector of the same length as
infoshould indicate which analyses will have a bound.- r
Integer value controlling grid for numerical integration as in Jennison and Turnbull (2000); default is 18, range is 1 to 80. Larger values provide larger number of grid points and greater accuracy. Normally
rwill not be changed by the user.- tol
Tolerance parameter for convergence (on Z-scale).
References
Jennison C and Turnbull BW (2000), Group Sequential Methods with Applications to Clinical Trials. Boca Raton: Chapman and Hall.
Author
Keaven Anderson keaven_anderson@merck.com
Examples
gs_power_ahr(
analysis_time = c(12, 24, 36),
event = c(30, 40, 50),
binding = TRUE,
upper = gs_spending_bound,
upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
lower = gs_spending_bound,
lpar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL)
)
#> $input
#> $input$enroll_rate
#> # A tibble: 3 × 3
#> stratum duration rate
#> <chr> <dbl> <dbl>
#> 1 All 2 3
#> 2 All 2 6
#> 3 All 10 9
#>
#> $input$fail_rate
#> # A tibble: 2 × 5
#> stratum duration fail_rate dropout_rate hr
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 All 3 0.0770 0.001 0.9
#> 2 All 100 0.0385 0.001 0.6
#>
#> $input$event
#> [1] 30 40 50
#>
#> $input$analysis_time
#> [1] 12 24 36
#>
#> $input$upper
#> function (k = 1, par = list(sf = gsDesign::sfLDOF, total_spend = 0.025,
#> param = NULL, timing = NULL, max_info = NULL), hgm1 = NULL,
#> theta = 0.1, info = 1:3, efficacy = TRUE, test_bound = TRUE,
#> r = 18, tol = 1e-06)
#> {
#> if (length(test_bound) == 1 && k > 1) {
#> test_bound <- rep(test_bound, k)
#> }
#> if (!is.null(par$timing)) {
#> timing <- par$timing
#> }
#> else {
#> if (is.null(par$max_info)) {
#> timing <- info/max(info)
#> }
#> else {
#> timing <- info/par$max_info
#> }
#> }
#> spend <- par$sf(alpha = par$total_spend, t = timing, param = par$param)$spend
#> old_spend <- 0
#> for (i in 1:k) {
#> if (test_bound[i]) {
#> xx <- spend[i] - old_spend
#> old_spend <- spend[i]
#> spend[i] <- xx
#> }
#> else {
#> spend[i] <- 0
#> }
#> }
#> spend <- spend[k]
#> if (!efficacy) {
#> if (spend <= 0) {
#> return(-Inf)
#> }
#> if (length(theta) == 1)
#> theta <- rep(theta, length(info))
#> a <- qnorm(spend) + sqrt(info[k]) * theta[k]
#> if (k == 1) {
#> return(a)
#> }
#> mu <- theta[k] * sqrt(info[k])
#> extreme_low <- mu - 3 - 4 * log(r)
#> extreme_high <- mu + 3 + 4 * log(r)
#> adelta <- 1
#> j <- 0
#> while (abs(adelta) > tol) {
#> hg <- hupdate(theta = theta[k], info = info[k], a = -Inf,
#> b = a, thetam1 = theta[k - 1], im1 = info[k -
#> 1], gm1 = hgm1, r = r)
#> i <- length(hg$h)
#> pik <- sum(hg$h)
#> adelta <- spend - pik
#> dplo <- hg$h[i]/hg$w[i]
#> if (adelta > dplo) {
#> adelta <- 1
#> }
#> else if (adelta < -dplo) {
#> adelta <- -1
#> }
#> else {
#> adelta <- adelta/dplo
#> }
#> a <- a + adelta
#> if (a > extreme_high) {
#> a <- extreme_high
#> }
#> else if (a < extreme_low) {
#> a <- extreme_low
#> }
#> if (abs(adelta) < tol) {
#> return(a)
#> }
#> j <- j + 1
#> if (j > 20) {
#> stop(paste("gs_spending_bound(): bound_update did not converge for lower bound calculation, analysis",
#> k, " !"))
#> }
#> }
#> }
#> else {
#> if (spend <= 0) {
#> return(Inf)
#> }
#> if (length(theta) == 1)
#> theta <- rep(theta, length(info))
#> b <- qnorm(spend, lower.tail = FALSE)
#> if (k == 1) {
#> return(b)
#> }
#> mu <- theta[k] * sqrt(info[k])
#> extreme_low <- mu - 3 - 4 * log(r)
#> extreme_high <- mu + 3 + 4 * log(r)
#> bdelta <- 1
#> j <- 1
#> while (abs(bdelta) > tol) {
#> hg <- hupdate(theta = 0, info = info[k], a = b, b = Inf,
#> thetam1 = 0, im1 = info[k - 1], gm1 = hgm1, r = r)
#> pik <- sum(hg$h)
#> bdelta <- spend - pik
#> dpikdb <- hg$h[1]/hg$w[1]
#> if (bdelta > dpikdb) {
#> bdelta <- 1
#> }
#> else if (bdelta < -dpikdb) {
#> bdelta <- -1
#> }
#> else {
#> bdelta <- bdelta/dpikdb
#> }
#> b <- b - bdelta
#> if (b > extreme_high) {
#> b <- extreme_high
#> }
#> else if (b < extreme_low) {
#> b <- extreme_low
#> }
#> if (abs(bdelta) < tol) {
#> return(b)
#> }
#> j <- j + 1
#> if (j > 20) {
#> stop(paste("gs_spending_bound(): bound_update did not converge for lower bound calculation, analysis",
#> k, " !"))
#> }
#> }
#> }
#> }
#> <bytecode: 0x56342efe2a08>
#> <environment: namespace:gsDesign2>
#>
#> $input$upar
#> $input$upar$sf
#> function (alpha, t, param = NULL)
#> {
#> checkScalar(alpha, "numeric", c(0, Inf), c(FALSE, FALSE))
#> checkVector(t, "numeric", c(0, Inf), c(TRUE, FALSE))
#> if (is.null(param) || param < 0.005 || param > 20)
#> param <- 1
#> checkScalar(param, "numeric", c(0.005, 20), c(TRUE, TRUE))
#> t[t > 1] <- 1
#> if (param == 1) {
#> rho <- 1
#> txt <- "Lan-DeMets O'Brien-Fleming approximation"
#> parname <- "none"
#> }
#> else {
#> rho <- param
#> txt <- "Generalized Lan-DeMets O'Brien-Fleming"
#> parname <- "rho"
#> }
#> z <- -qnorm(alpha/2)
#> x <- list(name = txt, param = param, parname = parname, sf = sfLDOF,
#> spend = 2 * (1 - pnorm(z/t^(rho/2))), bound = NULL, prob = NULL)
#> class(x) <- "spendfn"
#> x
#> }
#> <bytecode: 0x56342d0cba30>
#> <environment: namespace:gsDesign>
#>
#> $input$upar$total_spend
#> [1] 0.025
#>
#> $input$upar$param
#> NULL
#>
#> $input$upar$timing
#> NULL
#>
#>
#> $input$lower
#> function (k = 1, par = list(sf = gsDesign::sfLDOF, total_spend = 0.025,
#> param = NULL, timing = NULL, max_info = NULL), hgm1 = NULL,
#> theta = 0.1, info = 1:3, efficacy = TRUE, test_bound = TRUE,
#> r = 18, tol = 1e-06)
#> {
#> if (length(test_bound) == 1 && k > 1) {
#> test_bound <- rep(test_bound, k)
#> }
#> if (!is.null(par$timing)) {
#> timing <- par$timing
#> }
#> else {
#> if (is.null(par$max_info)) {
#> timing <- info/max(info)
#> }
#> else {
#> timing <- info/par$max_info
#> }
#> }
#> spend <- par$sf(alpha = par$total_spend, t = timing, param = par$param)$spend
#> old_spend <- 0
#> for (i in 1:k) {
#> if (test_bound[i]) {
#> xx <- spend[i] - old_spend
#> old_spend <- spend[i]
#> spend[i] <- xx
#> }
#> else {
#> spend[i] <- 0
#> }
#> }
#> spend <- spend[k]
#> if (!efficacy) {
#> if (spend <= 0) {
#> return(-Inf)
#> }
#> if (length(theta) == 1)
#> theta <- rep(theta, length(info))
#> a <- qnorm(spend) + sqrt(info[k]) * theta[k]
#> if (k == 1) {
#> return(a)
#> }
#> mu <- theta[k] * sqrt(info[k])
#> extreme_low <- mu - 3 - 4 * log(r)
#> extreme_high <- mu + 3 + 4 * log(r)
#> adelta <- 1
#> j <- 0
#> while (abs(adelta) > tol) {
#> hg <- hupdate(theta = theta[k], info = info[k], a = -Inf,
#> b = a, thetam1 = theta[k - 1], im1 = info[k -
#> 1], gm1 = hgm1, r = r)
#> i <- length(hg$h)
#> pik <- sum(hg$h)
#> adelta <- spend - pik
#> dplo <- hg$h[i]/hg$w[i]
#> if (adelta > dplo) {
#> adelta <- 1
#> }
#> else if (adelta < -dplo) {
#> adelta <- -1
#> }
#> else {
#> adelta <- adelta/dplo
#> }
#> a <- a + adelta
#> if (a > extreme_high) {
#> a <- extreme_high
#> }
#> else if (a < extreme_low) {
#> a <- extreme_low
#> }
#> if (abs(adelta) < tol) {
#> return(a)
#> }
#> j <- j + 1
#> if (j > 20) {
#> stop(paste("gs_spending_bound(): bound_update did not converge for lower bound calculation, analysis",
#> k, " !"))
#> }
#> }
#> }
#> else {
#> if (spend <= 0) {
#> return(Inf)
#> }
#> if (length(theta) == 1)
#> theta <- rep(theta, length(info))
#> b <- qnorm(spend, lower.tail = FALSE)
#> if (k == 1) {
#> return(b)
#> }
#> mu <- theta[k] * sqrt(info[k])
#> extreme_low <- mu - 3 - 4 * log(r)
#> extreme_high <- mu + 3 + 4 * log(r)
#> bdelta <- 1
#> j <- 1
#> while (abs(bdelta) > tol) {
#> hg <- hupdate(theta = 0, info = info[k], a = b, b = Inf,
#> thetam1 = 0, im1 = info[k - 1], gm1 = hgm1, r = r)
#> pik <- sum(hg$h)
#> bdelta <- spend - pik
#> dpikdb <- hg$h[1]/hg$w[1]
#> if (bdelta > dpikdb) {
#> bdelta <- 1
#> }
#> else if (bdelta < -dpikdb) {
#> bdelta <- -1
#> }
#> else {
#> bdelta <- bdelta/dpikdb
#> }
#> b <- b - bdelta
#> if (b > extreme_high) {
#> b <- extreme_high
#> }
#> else if (b < extreme_low) {
#> b <- extreme_low
#> }
#> if (abs(bdelta) < tol) {
#> return(b)
#> }
#> j <- j + 1
#> if (j > 20) {
#> stop(paste("gs_spending_bound(): bound_update did not converge for lower bound calculation, analysis",
#> k, " !"))
#> }
#> }
#> }
#> }
#> <bytecode: 0x56342efe2a08>
#> <environment: namespace:gsDesign2>
#>
#> $input$lpar
#> $input$lpar$sf
#> function (alpha, t, param = NULL)
#> {
#> checkScalar(alpha, "numeric", c(0, Inf), c(FALSE, FALSE))
#> checkVector(t, "numeric", c(0, Inf), c(TRUE, FALSE))
#> if (is.null(param) || param < 0.005 || param > 20)
#> param <- 1
#> checkScalar(param, "numeric", c(0.005, 20), c(TRUE, TRUE))
#> t[t > 1] <- 1
#> if (param == 1) {
#> rho <- 1
#> txt <- "Lan-DeMets O'Brien-Fleming approximation"
#> parname <- "none"
#> }
#> else {
#> rho <- param
#> txt <- "Generalized Lan-DeMets O'Brien-Fleming"
#> parname <- "rho"
#> }
#> z <- -qnorm(alpha/2)
#> x <- list(name = txt, param = param, parname = parname, sf = sfLDOF,
#> spend = 2 * (1 - pnorm(z/t^(rho/2))), bound = NULL, prob = NULL)
#> class(x) <- "spendfn"
#> x
#> }
#> <bytecode: 0x56342d0cba30>
#> <environment: namespace:gsDesign>
#>
#> $input$lpar$total_spend
#> [1] 0.025
#>
#> $input$lpar$param
#> NULL
#>
#> $input$lpar$timing
#> NULL
#>
#>
#> $input$test_lower
#> [1] TRUE
#>
#> $input$test_upper
#> [1] TRUE
#>
#> $input$ratio
#> [1] 1
#>
#> $input$binding
#> [1] TRUE
#>
#> $input$info_scale
#> [1] "h0_h1_info"
#>
#> $input$r
#> [1] 18
#>
#> $input$tol
#> [1] 1e-06
#>
#>
#> $enroll_rate
#> # A tibble: 3 × 3
#> stratum duration rate
#> <chr> <dbl> <dbl>
#> 1 All 2 3
#> 2 All 2 6
#> 3 All 10 9
#>
#> $fail_rate
#> # A tibble: 2 × 5
#> stratum duration fail_rate dropout_rate hr
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 All 3 0.0770 0.001 0.9
#> 2 All 100 0.0385 0.001 0.6
#>
#> $bound
#> # A tibble: 6 × 7
#> analysis bound probability probability0 z `~hr at bound` `nominal p`
#> <int> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 upper 0.00706 0.000867 3.13 0.316 0.000867
#> 2 1 lower 0.000935 0.00658 -2.48 2.49 0.993
#> 3 2 upper 0.115 0.00921 2.37 0.505 0.00892
#> 4 2 lower 0.00912 0.113 -1.21 1.42 0.888
#> 5 3 upper 0.324 0.0250 2.01 0.607 0.0222
#> 6 3 lower 0.0251 0.323 -0.474 1.12 0.682
#>
#> $analysis
#> analysis time n event ahr theta info info0
#> 1 1 14.90817 108 30.00008 0.7865726 0.2400702 7.373433 7.50002
#> 2 2 24.00000 108 49.06966 0.7151566 0.3352538 11.999266 12.26741
#> 3 3 36.00000 108 66.23948 0.6833395 0.3807634 16.267921 16.55987
#> info_frac info_frac0
#> 1 0.4532499 0.4529033
#> 2 0.7376029 0.7407917
#> 3 1.0000000 1.0000000
#>
#> attr(,"class")
#> [1] "ahr" "gs_design" "list"