Group sequential design using weighted log-rank test under non-proportional hazards
Source:R/gs_design_wlr.R
gs_design_wlr.RdGroup 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_b,
upar = gsDesign(k = 3, test.type = 1, n.I = c(0.25, 0.75, 1), sfu = sfLDOF, sfupar =
NULL)$upper$bound,
lower = gs_b,
lpar = c(qnorm(0.1), -Inf, -Inf),
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
FALSEis 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 asinfoshould 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 valueFALSEindicated no lower bound; otherwise, a logical vector of the same length asinfoshould 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,
rwill 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.
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: 0x56342eb4ad48>
#>
#> $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: 0x56342efe38e8>
#> <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: 0x56342efe38e8>
#> <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: 0x56342eb4ad48>
#>
#> $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: 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$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.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"