Skip to contents

Group sequential design using average hazard ratio under non-proportional hazards

Usage

gs_design_ahr(
  enroll_rate = define_enroll_rate(duration = c(2, 2, 10), rate = c(3, 6, 9)),
  fail_rate = define_fail_rate(duration = c(3, 100), fail_rate = log(2)/c(9, 18), hr =
    c(0.9, 0.6), dropout_rate = 0.001),
  alpha = 0.025,
  beta = 0.1,
  info_frac = NULL,
  analysis_time = 36,
  ratio = 1,
  binding = FALSE,
  upper = gs_b,
  upar = gsDesign::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),
  h1_spending = TRUE,
  test_upper = TRUE,
  test_lower = TRUE,
  info_scale = c("h0_h1_info", "h0_info", "h1_info"),
  r = 18,
  tol = 1e-06,
  interval = c(0.01, 1000)
)

Arguments

enroll_rate

Enrollment rates.

fail_rate

Failure and dropout rates.

alpha

One-sided Type I error.

beta

Type II error.

info_frac

Targeted information fraction at each analysis.

analysis_time

Minimum time of analysis.

ratio

Experimental:Control randomization ratio (not yet implemented).

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.

h1_spending

Indicator that lower bound to be set by spending under alternate hypothesis (input fail_rate) if spending is used for lower bound.

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.

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.

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.

Details

To be added.

Specification

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

Examples

library(gsDesign)
#> 
#> Attaching package: ‘gsDesign’
#> The following objects are masked from ‘package:gsDesign2’:
#> 
#>     as_gt, as_rtf
library(gsDesign2)
library(dplyr)

# Example 1 ----
# call with defaults
gs_design_ahr()
#> $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$alpha
#> [1] 0.025
#> 
#> $input$beta
#> [1] 0.1
#> 
#> $input$ratio
#> [1] 1
#> 
#> $input$info_frac
#> [1] 1
#> 
#> $input$analysis_time
#> [1] 36
#> 
#> $input$upper
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$upar
#> [1] 4.332634 2.339816 2.011793
#> 
#> $input$lower
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$lpar
#> [1] -1.281552      -Inf      -Inf
#> 
#> $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: 3 × 3
#>   stratum duration  rate
#>   <chr>      <dbl> <dbl>
#> 1 All            2  13.2
#> 2 All            2  26.4
#> 3 All           10  39.7
#> 
#> $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: 1 × 7
#>   analysis bound probability probability0     z `~hr at bound` `nominal p`
#>      <dbl> <chr>       <dbl>        <dbl> <dbl>          <dbl>       <dbl>
#> 1        1 upper         0.9        0.025  1.96          0.795      0.0250
#> 
#> $analysis
#> # A tibble: 1 × 9
#>   analysis  time     n event   ahr theta  info info0 info_frac
#>      <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>     <dbl>
#> 1        1    36  476.  292. 0.683 0.381  71.7  73.0         1
#> 
#> attr(,"class")
#> [1] "non_binding" "ahr"         "gs_design"   "list"       

# Example 2 ----
# Single analysis
gs_design_ahr(analysis_time = 40)
#> $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$alpha
#> [1] 0.025
#> 
#> $input$beta
#> [1] 0.1
#> 
#> $input$ratio
#> [1] 1
#> 
#> $input$info_frac
#> [1] 1
#> 
#> $input$analysis_time
#> [1] 40
#> 
#> $input$upper
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$upar
#> [1] 4.332634 2.339816 2.011793
#> 
#> $input$lower
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$lpar
#> [1] -1.281552      -Inf      -Inf
#> 
#> $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: 3 × 3
#>   stratum duration  rate
#>   <chr>      <dbl> <dbl>
#> 1 All            2  11.9
#> 2 All            2  23.8
#> 3 All           10  35.6
#> 
#> $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: 1 × 7
#>   analysis bound probability probability0     z `~hr at bound` `nominal p`
#>      <dbl> <chr>       <dbl>        <dbl> <dbl>          <dbl>       <dbl>
#> 1        1 upper         0.9        0.025  1.96          0.791      0.0250
#> 
#> $analysis
#> # A tibble: 1 × 9
#>   analysis  time     n event   ahr theta  info info0 info_frac
#>      <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>     <dbl>
#> 1        1    40  428.  280. 0.678 0.389  68.8  69.9         1
#> 
#> attr(,"class")
#> [1] "non_binding" "ahr"         "gs_design"   "list"       

# Example 3 ----
# Multiple analysis_time
gs_design_ahr(analysis_time = c(12, 24, 36))
#> $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$alpha
#> [1] 0.025
#> 
#> $input$beta
#> [1] 0.1
#> 
#> $input$ratio
#> [1] 1
#> 
#> $input$info_frac
#> [1] 0.3080415 0.7407917 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: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$upar
#> [1] 4.332634 2.339816 2.011793
#> 
#> $input$lower
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$lpar
#> [1] -1.281552      -Inf      -Inf
#> 
#> $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: 3 × 3
#>   stratum duration  rate
#>   <chr>      <dbl> <dbl>
#> 1 All            2  14.0
#> 2 All            2  27.9
#> 3 All           10  41.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: 4 × 7
#>   analysis bound probability probability0     z `~hr at bound` `nominal p`
#>      <int> <chr>       <dbl>        <dbl> <dbl>          <dbl>       <dbl>
#> 1        1 upper    0.000507   0.00000737  4.33          0.411  0.00000737
#> 2        1 lower    0.0111     0.100      -1.28          1.30   0.9       
#> 3        2 upper    0.566      0.00965     2.34          0.734  0.00965   
#> 4        3 upper    0.900      0.0251      2.01          0.795  0.0221    
#> 
#> $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  419.  95.0 0.811 0.210  23.4  23.8     0.309
#> 2        2    24  503. 228.  0.715 0.335  55.9  57.1     0.738
#> 3        3    36  503. 308.  0.683 0.381  75.8  77.1     1    
#> 
#> attr(,"class")
#> [1] "non_binding" "ahr"         "gs_design"   "list"       

# Example 4 ----
# Specified information fraction
# \donttest{
gs_design_ahr(info_frac = c(.25, .75, 1), analysis_time = 36)
#> $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$alpha
#> [1] 0.025
#> 
#> $input$beta
#> [1] 0.1
#> 
#> $input$ratio
#> [1] 1
#> 
#> $input$info_frac
#> [1] 0.25 0.75 1.00
#> 
#> $input$analysis_time
#> [1] 36
#> 
#> $input$upper
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$upar
#> [1] 4.332634 2.339816 2.011793
#> 
#> $input$lower
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$lpar
#> [1] -1.281552      -Inf      -Inf
#> 
#> $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: 3 × 3
#>   stratum duration  rate
#>   <chr>      <dbl> <dbl>
#> 1 All            2  14.2
#> 2 All            2  28.4
#> 3 All           10  42.5
#> 
#> $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: 4 × 7
#>   analysis bound probability probability0     z `~hr at bound` `nominal p`
#>      <int> <chr>       <dbl>        <dbl> <dbl>          <dbl>       <dbl>
#> 1        1 upper    0.000282   0.00000737  4.33          0.376  0.00000737
#> 2        1 lower    0.0166     0.100      -1.28          1.34   0.9       
#> 3        2 upper    0.585      0.00965     2.34          0.737  0.00965   
#> 4        3 upper    0.900      0.0249      2.01          0.797  0.0221    
#> 
#> $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  10.7  371.  78.3 0.823 0.195  19.3  19.6     0.251
#> 2        2  24.4  510. 235.  0.714 0.337  57.4  58.7     0.747
#> 3        3  36    510. 313.  0.683 0.381  76.9  78.3     1    
#> 
#> attr(,"class")
#> [1] "non_binding" "ahr"         "gs_design"   "list"       
# }

# Example 5 ----
# multiple analysis times & info_frac
# driven by times
gs_design_ahr(info_frac = c(.25, .75, 1), analysis_time = c(12, 25, 36))
#> $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$alpha
#> [1] 0.025
#> 
#> $input$beta
#> [1] 0.1
#> 
#> $input$ratio
#> [1] 1
#> 
#> $input$info_frac
#> [1] 0.25 0.75 1.00
#> 
#> $input$analysis_time
#> [1] 12 25 36
#> 
#> $input$upper
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$upar
#> [1] 4.332634 2.339816 2.011793
#> 
#> $input$lower
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$lpar
#> [1] -1.281552      -Inf      -Inf
#> 
#> $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: 3 × 3
#>   stratum duration  rate
#>   <chr>      <dbl> <dbl>
#> 1 All            2  14.0
#> 2 All            2  27.9
#> 3 All           10  41.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: 4 × 7
#>   analysis bound probability probability0     z `~hr at bound` `nominal p`
#>      <int> <chr>       <dbl>        <dbl> <dbl>          <dbl>       <dbl>
#> 1        1 upper    0.000507   0.00000737  4.33          0.411  0.00000737
#> 2        1 lower    0.0111     0.100      -1.28          1.30   0.9       
#> 3        2 upper    0.600      0.00965     2.34          0.738  0.00965   
#> 4        3 upper    0.900      0.0248      2.01          0.795  0.0221    
#> 
#> $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  419.  95.0 0.811 0.210  23.4  23.7     0.309
#> 2        2    25  503. 236.  0.711 0.341  57.8  59.1     0.763
#> 3        3    36  503. 308.  0.683 0.381  75.7  77.1     1    
#> 
#> attr(,"class")
#> [1] "non_binding" "ahr"         "gs_design"   "list"       
# driven by info_frac
# \donttest{
gs_design_ahr(info_frac = c(1 / 3, .8, 1), analysis_time = c(12, 25, 36))
#> $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$alpha
#> [1] 0.025
#> 
#> $input$beta
#> [1] 0.1
#> 
#> $input$ratio
#> [1] 1
#> 
#> $input$info_frac
#> [1] 0.3333333 0.8000000 1.0000000
#> 
#> $input$analysis_time
#> [1] 12 25 36
#> 
#> $input$upper
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$upar
#> [1] 4.332634 2.339816 2.011793
#> 
#> $input$lower
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$lpar
#> [1] -1.281552      -Inf      -Inf
#> 
#> $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: 3 × 3
#>   stratum duration  rate
#>   <chr>      <dbl> <dbl>
#> 1 All            2  13.9
#> 2 All            2  27.8
#> 3 All           10  41.7
#> 
#> $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: 4 × 7
#>   analysis bound probability probability0     z `~hr at bound` `nominal p`
#>      <int> <chr>       <dbl>        <dbl> <dbl>          <dbl>       <dbl>
#> 1        1 upper    0.000645   0.00000737  4.33          0.425  0.00000737
#> 2        1 lower    0.00928    0.100      -1.28          1.29   0.9       
#> 3        2 upper    0.640      0.00965     2.34          0.742  0.00965   
#> 4        3 upper    0.900      0.0244      2.01          0.795  0.0221    
#> 
#> $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.5  439.  102. 0.806 0.216  25.2  25.6     0.334
#> 2        2  26.4  501.  246. 0.706 0.348  60.1  61.4     0.797
#> 3        3  36    501.  307. 0.683 0.381  75.4  76.8     1    
#> 
#> attr(,"class")
#> [1] "non_binding" "ahr"         "gs_design"   "list"       
# }

# Example 6 ----
# 2-sided symmetric design with O'Brien-Fleming spending
# \donttest{
gs_design_ahr(
  analysis_time = c(12, 24, 36),
  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),
  h1_spending = FALSE
)
#> $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$alpha
#> [1] 0.025
#> 
#> $input$beta
#> [1] 0.1
#> 
#> $input$ratio
#> [1] 1
#> 
#> $input$info_frac
#> [1] 0.3080415 0.7407917 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: 0x557dd1cabf28>
#> <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: 0x557dd108f3d8>
#> <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: 0x557dd1cabf28>
#> <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: 0x557dd108f3d8>
#> <environment: namespace:gsDesign>
#> 
#> $input$lpar$total_spend
#> [1] 0.025
#> 
#> $input$lpar$param
#> NULL
#> 
#> $input$lpar$timing
#> NULL
#> 
#> 
#> $input$test_upper
#> [1] TRUE
#> 
#> $input$test_lower
#> [1] TRUE
#> 
#> $input$h1_spending
#> [1] FALSE
#> 
#> $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  13.7
#> 2 All            2  27.5
#> 3 All           10  41.2
#> 
#> $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.00226        0.0000538  3.87          0.449   0.0000538
#> 2        1 lower 0.000000613    0.0000538 -3.87          2.23    1.00     
#> 3        2 upper 0.550          0.00921    2.36          0.730   0.00919  
#> 4        2 lower 0.00000125     0.00921   -2.36          1.37    0.991    
#> 5        3 upper 0.900          0.0250     2.01          0.794   0.0222   
#> 6        3 lower 0.00000128     0.0250    -2.01          1.26    0.978    
#> 
#> $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  412.  93.4 0.811 0.210  23.0  23.3     0.309
#> 2        2    24  494. 224.  0.715 0.335  54.9  56.1     0.738
#> 3        3    36  494. 303.  0.683 0.381  74.4  75.8     1    
#> 
#> attr(,"class")
#> [1] "ahr"       "gs_design" "list"     
# }
# 2-sided asymmetric design with O'Brien-Fleming upper spending
# Pocock lower spending under H1 (NPH)
# \donttest{
gs_design_ahr(
  analysis_time = c(12, 24, 36),
  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::sfLDPocock, total_spend = 0.1, param = NULL, timing = NULL),
  h1_spending = TRUE
)
#> $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$alpha
#> [1] 0.025
#> 
#> $input$beta
#> [1] 0.1
#> 
#> $input$ratio
#> [1] 1
#> 
#> $input$info_frac
#> [1] 0.3080415 0.7407917 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: 0x557dd1cabf28>
#> <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: 0x557dd108f3d8>
#> <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: 0x557dd1cabf28>
#> <environment: namespace:gsDesign2>
#> 
#> $input$lpar
#> $input$lpar$sf
#> function (alpha, t, param) 
#> {
#>     checkScalar(alpha, "numeric", c(0, Inf), c(FALSE, FALSE))
#>     checkVector(t, "numeric", c(0, Inf), c(TRUE, FALSE))
#>     t[t > 1] <- 1
#>     x <- list(name = "Lan-DeMets Pocock approximation", param = NULL, 
#>         parname = "none", sf = sfLDPocock, spend = alpha * log(1 + 
#>             (exp(1) - 1) * t), bound = NULL, prob = NULL)
#>     class(x) <- "spendfn"
#>     x
#> }
#> <bytecode: 0x557dcb1bda40>
#> <environment: namespace:gsDesign>
#> 
#> $input$lpar$total_spend
#> [1] 0.1
#> 
#> $input$lpar$param
#> NULL
#> 
#> $input$lpar$timing
#> NULL
#> 
#> 
#> $input$test_upper
#> [1] TRUE
#> 
#> $input$test_lower
#> [1] TRUE
#> 
#> $input$h1_spending
#> [1] TRUE
#> 
#> $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  16.5
#> 2 All            2  32.9
#> 3 All           10  49.4
#> 
#> $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.00305    0.0000538  3.87           0.481   0.0000538
#> 2        1 lower     0.0430     0.268     -0.619          1.12    0.732    
#> 3        2 upper     0.638      0.00921    2.36           0.750   0.00920  
#> 4        2 lower     0.0823     0.874      1.13           0.871   0.129    
#> 5        3 upper     0.900      0.0250     1.98           0.813   0.0240   
#> 6        3 lower     0.100      0.975      1.97           0.813   0.0243   
#> 
#> $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  494.  112. 0.811 0.210  27.6  28.0     0.309
#> 2        2    24  593.  269. 0.715 0.335  65.9  67.3     0.738
#> 3        3    36  593.  364. 0.683 0.381  89.3  90.9     1    
#> 
#> attr(,"class")
#> [1] "ahr"       "gs_design" "list"     
# }

# Example 7 ----
# \donttest{
gs_design_ahr(
  alpha = 0.0125,
  analysis_time = c(12, 24, 36),
  upper = gs_spending_bound,
  upar = list(sf = gsDesign::sfLDOF, total_spend = 0.0125, param = NULL, timing = NULL),
  lower = gs_b,
  lpar = rep(-Inf, 3)
)
#> $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$alpha
#> [1] 0.0125
#> 
#> $input$beta
#> [1] 0.1
#> 
#> $input$ratio
#> [1] 1
#> 
#> $input$info_frac
#> [1] 0.3080415 0.7407917 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: 0x557dd1cabf28>
#> <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: 0x557dd108f3d8>
#> <environment: namespace:gsDesign>
#> 
#> $input$upar$total_spend
#> [1] 0.0125
#> 
#> $input$upar$param
#> NULL
#> 
#> $input$upar$timing
#> NULL
#> 
#> 
#> $input$lower
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$lpar
#> [1] -Inf -Inf -Inf
#> 
#> $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: 3 × 3
#>   stratum duration  rate
#>   <chr>      <dbl> <dbl>
#> 1 All            2  16.1
#> 2 All            2  32.2
#> 3 All           10  48.3
#> 
#> $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: 3 × 7
#>   analysis bound probability probability0     z `~hr at bound` `nominal p`
#>      <int> <chr>       <dbl>        <dbl> <dbl>          <dbl>       <dbl>
#> 1        1 upper    0.000619   0.00000679  4.35          0.435  0.00000679
#> 2        2 upper    0.505      0.00371     2.68          0.719  0.00371   
#> 3        3 upper    0.900      0.0125      2.28          0.785  0.0114    
#> 
#> $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  483.  109. 0.811 0.210  27.0  27.4     0.309
#> 2        2    24  579.  263. 0.715 0.335  64.3  65.8     0.738
#> 3        3    36  579.  355. 0.683 0.381  87.2  88.8     1    
#> 
#> attr(,"class")
#> [1] "non_binding" "ahr"         "gs_design"   "list"       

gs_design_ahr(
  alpha = 0.0125,
  analysis_time = c(12, 24, 36),
  upper = gs_b,
  upar = gsDesign::gsDesign(
    k = 3, test.type = 1, n.I = c(.25, .75, 1),
    sfu = sfLDOF, sfupar = NULL, alpha = 0.0125
  )$upper$bound,
  lower = gs_b,
  lpar = rep(-Inf, 3)
)
#> $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$alpha
#> [1] 0.0125
#> 
#> $input$beta
#> [1] 0.1
#> 
#> $input$ratio
#> [1] 1
#> 
#> $input$info_frac
#> [1] 0.3080415 0.7407917 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: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$upar
#> [1] 4.859940 2.658446 2.280095
#> 
#> $input$lower
#> function (par = NULL, k = NULL, ...) 
#> {
#>     if (is.null(k)) {
#>         return(par)
#>     }
#>     else {
#>         return(par[k])
#>     }
#> }
#> <bytecode: 0x557dd1ca8fd8>
#> <environment: namespace:gsDesign2>
#> 
#> $input$lpar
#> [1] -Inf -Inf -Inf
#> 
#> $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: 3 × 3
#>   stratum duration  rate
#>   <chr>      <dbl> <dbl>
#> 1 All            2  16.1
#> 2 All            2  32.2
#> 3 All           10  48.3
#> 
#> $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: 3 × 7
#>   analysis bound probability probability0     z `~hr at bound` `nominal p`
#>      <int> <chr>       <dbl>        <dbl> <dbl>          <dbl>       <dbl>
#> 1        1 upper   0.0000938  0.000000587  4.86          0.395 0.000000587
#> 2        2 upper   0.513      0.00393      2.66          0.721 0.00393    
#> 3        3 upper   0.900      0.0125       2.28          0.785 0.0113     
#> 
#> $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  483.  110. 0.811 0.210  27.0  27.4     0.309
#> 2        2    24  580.  263. 0.715 0.335  64.4  65.9     0.738
#> 3        3    36  580.  356. 0.683 0.381  87.3  88.9     1    
#> 
#> attr(,"class")
#> [1] "non_binding" "ahr"         "gs_design"   "list"       
# }