Skip to contents

Group sequential design using weighted log-rank test under non-proportional hazards

Usage

gs_design_wlr(
  enroll_rate = define_enroll_rate(duration = c(2, 2, 10), rate = c(3, 6, 9)),
  fail_rate = tibble(stratum = "All", duration = c(3, 100), fail_rate = log(2)/c(9, 18),
    hr = c(0.9, 0.6), dropout_rate = rep(0.001, 2)),
  weight = wlr_weight_fh,
  approx = "asymptotic",
  alpha = 0.025,
  beta = 0.1,
  ratio = 1,
  info_frac = NULL,
  info_scale = c("h0_h1_info", "h0_info", "h1_info"),
  analysis_time = 36,
  binding = FALSE,
  upper = gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = alpha),
  lower = gs_spending_bound,
  lpar = list(sf = gsDesign::sfLDOF, total_spend = beta),
  test_upper = TRUE,
  test_lower = TRUE,
  h1_spending = TRUE,
  r = 18,
  tol = 1e-06,
  interval = c(0.01, 1000)
)

Arguments

enroll_rate

Enrollment rates.

fail_rate

Failure and dropout rates.

weight

Weight of weighted log rank test:

  • "1" = unweighted.

  • "n" = Gehan-Breslow.

  • "sqrtN" = Tarone-Ware.

  • "FH_p[a]_q[b]" = Fleming-Harrington with p=a and q=b.

approx

Approximate estimation method for Z statistics.

  • "event_driven" = only work under proportional hazard model with log rank test.

  • "asymptotic".

alpha

One-sided Type I error.

beta

Type II error.

ratio

Experimental:Control randomization ratio (not yet implemented).

info_frac

Targeted information fraction at each analysis.

info_scale

Information scale for calculation. Options are:

  • "h0_h1_info" (default): variance under both null and alternative hypotheses is used.

  • "h0_info": variance under null hypothesis is used.

  • "h1_info": variance under alternative hypothesis is used.

analysis_time

Minimum time of analysis.

binding

Indicator of whether futility bound is binding; default of FALSE is recommended.

upper

Function to compute upper bound.

upar

Parameters passed to upper.

lower

Function to compute lower bound.

lpar

Parameters passed to lower.

test_upper

Indicator of which analyses should include an upper (efficacy) bound; single value of TRUE (default) indicates all analyses; otherwise, a logical vector of the same length as info should indicate which analyses will have an efficacy bound.

test_lower

Indicator of which analyses should include an lower bound; single value of TRUE (default) indicates all analyses; single value FALSE indicated no lower bound; otherwise, a logical vector of the same length as info should indicate which analyses will have a lower bound.

h1_spending

Indicator that lower bound to be set by spending under alternate hypothesis (input fail_rate) if spending is used for lower 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, r will not be changed by the user.

tol

Tolerance parameter for boundary convergence (on Z-scale).

interval

An interval that is presumed to include the time at which expected event count is equal to targeted event.

Value

A list with input parameters, enrollment rate, analysis, and bound.

Specification

The contents of this section are shown in PDF user manual only.

Examples

library(dplyr)
library(mvtnorm)
library(gsDesign)
library(gsDesign2)

# set enrollment rates
enroll_rate <- define_enroll_rate(duration = 12, rate = 500 / 12)

# set failure rates
fail_rate <- define_fail_rate(
  duration = c(4, 100),
  fail_rate = log(2) / 15, # median survival 15 month
  hr = c(1, .6),
  dropout_rate = 0.001
)

# Example 1 ----
# Boundary is fixed
x <- gsSurv(
  k = 3,
  test.type = 4,
  alpha = 0.025, beta = 0.2,
  astar = 0, timing = 1,
  sfu = sfLDOF, sfupar = 0,
  sfl = sfLDOF, sflpar = 0,
  lambdaC = 0.1,
  hr = 0.6, hr0 = 1,
  eta = 0.01, gamma = 10,
  R = 12, S = NULL,
  T = 36, minfup = 24,
  ratio = 1
)

gs_design_wlr(
  enroll_rate = enroll_rate,
  fail_rate = fail_rate,
  ratio = 1,
  alpha = 0.025, beta = 0.2,
  weight = function(x, arm0, arm1) {
    wlr_weight_fh(x, arm0, arm1, rho = 0, gamma = 0.5)
  },
  upper = gs_b,
  upar = x$upper$bound,
  lower = gs_b,
  lpar = x$lower$bound,
  analysis_time = c(12, 24, 36)
)
#> $input
#> $input$enroll_rate
#> # A tibble: 1 × 3
#>   stratum duration  rate
#>   <chr>      <dbl> <dbl>
#> 1 All           12  41.7
#> 
#> $input$fail_rate
#> # A tibble: 2 × 5
#>   stratum duration fail_rate dropout_rate    hr
#>   <chr>      <dbl>     <dbl>        <dbl> <dbl>
#> 1 All            4    0.0462        0.001   1  
#> 2 All          100    0.0462        0.001   0.6
#> 
#> $input$weight
#> function(x, arm0, arm1) {
#>     wlr_weight_fh(x, arm0, arm1, rho = 0, gamma = 0.5)
#>   }
#> <environment: 0x5650a44bcc38>
#> 
#> $input$approx
#> [1] "asymptotic"
#> 
#> $input$alpha
#> [1] 0.025
#> 
#> $input$beta
#> [1] 0.2
#> 
#> $input$ratio
#> [1] 1
#> 
#> $input$info_frac
#> [1] 0.3241690 0.7434051 1.0000000
#> 
#> $input$analysis_time
#> [1] 12 24 36
#> 
#> $input$upper
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x5650a3d7c118>
#> <environment: namespace:gsDesign2>
#> 
#> $input$upar
#> [1] 3.710303 2.511407 1.992970
#> 
#> $input$lower
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x5650a3d7c118>
#> <environment: namespace:gsDesign2>
#> 
#> $input$lpar
#> [1] -0.2361874  1.1703638  1.9929702
#> 
#> $input$test_upper
#> [1] TRUE
#> 
#> $input$test_lower
#> [1] TRUE
#> 
#> $input$h1_spending
#> [1] TRUE
#> 
#> $input$binding
#> [1] FALSE
#> 
#> $input$info_scale
#> [1] "h0_h1_info"
#> 
#> $input$r
#> [1] 18
#> 
#> $input$tol
#> [1] 1e-06
#> 
#> 
#> $enroll_rate
#> # A tibble: 1 × 3
#>   stratum duration  rate
#>   <chr>      <dbl> <dbl>
#> 1 All           12  30.6
#> 
#> $fail_rate
#> # A tibble: 2 × 5
#>   stratum duration fail_rate dropout_rate    hr
#>   <chr>      <dbl>     <dbl>        <dbl> <dbl>
#> 1 All            4    0.0462        0.001   1  
#> 2 All          100    0.0462        0.001   0.6
#> 
#> $bounds
#> # 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.00366     0.000104  3.71           0.434    0.000104
#> 2        1 lower     0.105       0.407    -0.236          1.05     0.593   
#> 3        2 upper     0.504       0.00587   2.51           0.688    0.00601 
#> 4        2 lower     0.158       0.894     1.17           0.840    0.121   
#> 5        3 upper     0.800       0.0193    1.99           0.775    0.0231  
#> 6        3 lower     0.200       0.981     1.99           0.775    0.0231  
#> 
#> $analysis
#> # A tibble: 3 × 9
#>   analysis  time     n event   ahr theta  info info0 info_frac
#>      <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>     <dbl>
#> 1        1    12  368.  78.9 0.781 0.626  2.65  2.66     0.133
#> 2        2    24  368. 181.  0.666 0.765 11.3  11.6      0.565
#> 3        3    36  368. 244.  0.639 0.732 20.0  20.9      1    
#> 
#> attr(,"class")
#> [1] "non_binding" "wlr"         "gs_design"   "list"       

# Example 2 ----
# Boundary derived by spending function
gs_design_wlr(
  enroll_rate = enroll_rate,
  fail_rate = fail_rate,
  ratio = 1,
  alpha = 0.025, beta = 0.2,
  weight = function(x, arm0, arm1) {
    wlr_weight_fh(x, arm0, arm1, rho = 0, gamma = 0.5)
  },
  upper = gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
  lower = gs_spending_bound,
  lpar = list(sf = gsDesign::sfLDOF, total_spend = 0.2),
  analysis_time = c(12, 24, 36)
)
#> $input
#> $input$enroll_rate
#> # A tibble: 1 × 3
#>   stratum duration  rate
#>   <chr>      <dbl> <dbl>
#> 1 All           12  41.7
#> 
#> $input$fail_rate
#> # A tibble: 2 × 5
#>   stratum duration fail_rate dropout_rate    hr
#>   <chr>      <dbl>     <dbl>        <dbl> <dbl>
#> 1 All            4    0.0462        0.001   1  
#> 2 All          100    0.0462        0.001   0.6
#> 
#> $input$weight
#> function(x, arm0, arm1) {
#>     wlr_weight_fh(x, arm0, arm1, rho = 0, gamma = 0.5)
#>   }
#> <environment: 0x5650a44bcc38>
#> 
#> $input$approx
#> [1] "asymptotic"
#> 
#> $input$alpha
#> [1] 0.025
#> 
#> $input$beta
#> [1] 0.2
#> 
#> $input$ratio
#> [1] 1
#> 
#> $input$info_frac
#> [1] 0.3241690 0.7434051 1.0000000
#> 
#> $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: 0x5650a3d7b238>
#> <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: 0x5650a1da1c50>
#> <environment: namespace:gsDesign>
#> 
#> $input$upar$total_spend
#> [1] 0.025
#> 
#> 
#> $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: 0x5650a3d7b238>
#> <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: 0x5650a1da1c50>
#> <environment: namespace:gsDesign>
#> 
#> $input$lpar$total_spend
#> [1] 0.2
#> 
#> 
#> $input$test_upper
#> [1] TRUE
#> 
#> $input$test_lower
#> [1] TRUE
#> 
#> $input$h1_spending
#> [1] TRUE
#> 
#> $input$binding
#> [1] FALSE
#> 
#> $input$info_scale
#> [1] "h0_h1_info"
#> 
#> $input$r
#> [1] 18
#> 
#> $input$tol
#> [1] 1e-06
#> 
#> 
#> $enroll_rate
#> # A tibble: 1 × 3
#>   stratum duration  rate
#>   <chr>      <dbl> <dbl>
#> 1 All           12  24.0
#> 
#> $fail_rate
#> # A tibble: 2 × 5
#>   stratum duration fail_rate dropout_rate    hr
#>   <chr>      <dbl>     <dbl>        <dbl> <dbl>
#> 1 All            4    0.0462        0.001   1  
#> 2 All          100    0.0462        0.001   0.6
#> 
#> $bounds
#> # 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.0000000732     3.30e-10  6.18           0.208    3.30e-10
#> 2        1 lower 0.000441         7.55e- 3 -2.43           1.85     9.92e- 1
#> 3        2 upper 0.301            2.57e- 3  2.80           0.625    2.57e- 3
#> 4        2 lower 0.0882           8.23e- 1  0.925          0.856    1.77e- 1
#> 5        3 upper 0.800            2.20e- 2  1.97           0.751    2.42e- 2
#> 6        3 lower 0.200            9.78e- 1  1.97           0.751    2.42e- 2
#> 
#> $analysis
#> # A tibble: 3 × 9
#>   analysis  time     n event   ahr theta  info info0 info_frac
#>      <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>     <dbl>
#> 1        1    12  288.  61.9 0.781 0.626  2.08  2.09     0.133
#> 2        2    24  288. 142.  0.666 0.765  8.86  9.07     0.565
#> 3        3    36  288. 191.  0.639 0.732 15.7  16.4      1    
#> 
#> attr(,"class")
#> [1] "non_binding" "wlr"         "gs_design"   "list"