From 24147650dae78e56f8a1bf90fa6187e6133b53d9 Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Wed, 7 Feb 2024 15:55:37 -0500 Subject: [PATCH 01/30] initial draft --- vignettes/articles/story-update-boundary.Rmd | 251 +++++++++++++++++++ 1 file changed, 251 insertions(+) create mode 100644 vignettes/articles/story-update-boundary.Rmd diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd new file mode 100644 index 00000000..b8a87583 --- /dev/null +++ b/vignettes/articles/story-update-boundary.Rmd @@ -0,0 +1,251 @@ +--- +title: "Efficacy and Futility Boundary Update" +author: "Yujie Zhao and Keaven M. Anderson" +output: + rmarkdown::html_document: + toc: true + toc_float: true + toc_depth: 2 + number_sections: true + highlight: "textmate" + css: "custom.css" + code_folding: "hide" +bibliography: "gsDesign2.bib" +vignette: > + %\VignetteIndexEntry{Efficacy and Futility Boundary Update} + %\VignetteEncoding{UTF-8} + %\VignetteEngine{knitr::rmarkdown} +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) + +library(gsDesign) +library(gsDesign2) +library(gt) +``` + +# Design assumptions +```{r} +alpha <- 0.025 +beta <- 0.1 +ratio <- 1 + +# enrollment +enroll_rate <- define_enroll_rate( + duration = c(2, 2, 10), + rate = c(3, 6, 9)) + +# failure and dropout +fail_rate <- define_fail_rate( + duration = c(3, 100), + fail_rate = log(2) / 9, + hr = c(1, 0.6), + dropout_rate = .0001) + +# upper and lower boundary +# take the 2-sided symmetric design as an example +upper <- gs_spending_bound +upar <- list(sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL, timing = NULL) + +lower <- gs_spending_bound +lpar <- list(sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL, timing = NULL) +``` + +# Derive group sequential design + +We assume there are totally 2 analyses, including FA. The IA happens at month 20, and FA occurs at month 36. +```{r} +x <- gs_design_ahr( + enroll_rate = enroll_rate, + fail_rate = fail_rate, + alpha = alpha, + beta = beta, + info_frac = NULL, + analysis_time = c(20, 36), + ratio = ratio, + upper = gs_spending_bound, + upar = upar, + lower = lower, + lpar = lpar, + binding = FALSE) |> to_integer() +``` + +In the planed design, we have + +- planned events: `r round(x$analysis$event, 0)` +- planned information fraction (timing): `r round(x$analysis$info_frac, 4)` +- planned alpha spending: `r sfLDOF(0.025, x$analysis$info_frac)$spend` + +```{r} +x |> + summary() |> + gt() |> + tab_header(title = "Original design") |> + tab_style( + style = list( + cell_fill(color = "lightcyan"), + cell_text(weight = "bold") + ), + locations = cells_body( + columns = `Alternate hypothesis`, + rows = Bound == "Efficacy" + )) |> + tab_style( + style = list( + cell_fill(color = "#F9E3D6"), + cell_text(weight = "bold") + ), + locations = cells_body( + columns = `Alternate hypothesis`, + rows = Bound == "Futility" + ) + ) +``` + +# Update bounds at time of analysis {.tabset} + +Assume in practice, there are 180 and 300 events observed in IA and FA, respectively. +```{r} +# observed vs. planned events +event_observed <- c(180, 300) +event_planned <- x$analysis$event +``` + +Now we calculate the updated timing based on the following 2 approaches. + +## Event fraction based + +The updated timing is `r round(c(event_observed[1] / max(event_planned), 1), 4)`, if we divide the observed IA events by the planed FA events. We will spend alpha at IA based on its actual information fraction. The rest of the alpha will be all spent to FA. + +```{r} +# updated spending time +upper_spending_time <- c(event_observed[1] / max(event_planned), 1) +lower_spending_time <- upper_spending_time + +# updated upar & lpar +upar_updated <- x$input$upar +upar_updated$timing <- upper_spending_time +lpar_updated <- x$input$lpar +lpar_updated$timing <- lower_spending_time +``` + +```{r} +# observed statistical information +info_observed <- gs_info_ahr( + enroll_rate = x$enroll_rate, + fail_rate = x$fail_rate, + ratio = ratio, + event = event_observed, + analysis_time = x$analysis$time) +``` + +```{r} +x_updated <- gs_power_npe( + theta = info_observed$theta, + theta0 = 0, + theta1 = x$analysis$theta, + info = info_observed$info, + info0 = x$analysis$info0, + info1 = x$analysis$info, + binding = x$input$binding, + upper = x$input$upper, + upar = upar_updated, + lower = x$input$lower, + lpar = lpar_updated, + test_upper = TRUE, + test_lower = TRUE) + +x_updated |> + gt() |> + tab_header(title = "Updated design", subtitle = "by event fraction based timing") |> + tab_style( + style = list( + cell_fill(color = "lightcyan"), + cell_text(weight = "bold") + ), + locations = cells_body( + columns = probability, + rows = bound == "upper" + ) + ) |> + tab_style( + style = list( + cell_fill(color = "#F9E3D6"), + cell_text(weight = "bold") + ), + locations = cells_body( + columns = probability, + rows = bound == "lower" + ) + ) +``` + +## Information fraction based + +The updated timing is `r round(info_observed$info / max(info_observed$info), 4)`, which is the information fraction derived by the `gs_info_ahr()` by inputing the observed events. + +```{r} +# updated spending time +upper_spending_time <- info_observed$info / max(info_observed$info) +lower_spending_time <- upper_spending_time + +# updated upar & lpar +upar_updated <- x$input$upar +upar_updated$timing <- upper_spending_time +lpar_updated <- x$input$lpar +lpar_updated$timing <- lower_spending_time +``` + +```{r} +# observed statistical information +info_observed <- gs_info_ahr( + enroll_rate = x$enroll_rate, + fail_rate = x$fail_rate, + ratio = ratio, + event = event_observed, + analysis_time = x$analysis$time) +``` + +```{r} +x_updated <- gs_power_npe( + theta = info_observed$theta, + theta0 = 0, + theta1 = x$analysis$theta, + info = info_observed$info, + info0 = x$analysis$info0, + info1 = x$analysis$info, + binding = x$input$binding, + upper = x$input$upper, + upar = upar_updated, + lower = x$input$lower, + lpar = lpar_updated, + test_upper = TRUE, + test_lower = TRUE) + +x_updated |> + gt() |> + tab_header(title = "Updated design", subtitle = "by information fraction based timing") |> + tab_style( + style = list( + cell_fill(color = "lightcyan"), + cell_text(weight = "bold") + ), + locations = cells_body( + columns = probability, + rows = bound == "upper" + ) + ) |> + tab_style( + style = list( + cell_fill(color = "#F9E3D6"), + cell_text(weight = "bold") + ), + locations = cells_body( + columns = probability, + rows = bound == "lower" + ) + ) +``` + +# References From a3e7fce92b85bda6811143b4aaa898ad731328d3 Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Tue, 13 Feb 2024 14:30:18 -0500 Subject: [PATCH 02/30] 2nd draft after discussion with Keaven --- vignettes/articles/story-update-boundary.Rmd | 282 ++++++++++++------- 1 file changed, 186 insertions(+), 96 deletions(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index b8a87583..bd49265a 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -9,7 +9,6 @@ output: number_sections: true highlight: "textmate" css: "custom.css" - code_folding: "hide" bibliography: "gsDesign2.bib" vignette: > %\VignetteIndexEntry{Efficacy and Futility Boundary Update} @@ -23,10 +22,15 @@ knitr::opts_chunk$set(echo = TRUE) library(gsDesign) library(gsDesign2) library(gt) +library(tibble) +library(dplyr) ``` # Design assumptions -```{r} + +We assume a total of two analyses, including a final Analysis (FA). The initial Analysis (IA) occurs at month 20, followed by the FA at month 36. The enrollment period spans 14 months, with the first 2 months having an enrollment rate of 3, the next 2 months with a rate of 6, and the final 10 months with a rate of 9. The planned enrollment rate to achieve the 90\% power is a multiptier of 3, 6, and 9. Additionally, there is a delayed treatment effect in the first 3 months (HR = 1), which transitions to an HR of 0.6 thereafter. The median of the control arm is 9 months, and the dropout rate is 0.0001. + +```{r, echo=FALSE} alpha <- 0.025 beta <- 0.1 ratio <- 1 @@ -42,20 +46,18 @@ fail_rate <- define_fail_rate( fail_rate = log(2) / 9, hr = c(1, 0.6), dropout_rate = .0001) +``` -# upper and lower boundary -# take the 2-sided symmetric design as an example -upper <- gs_spending_bound -upar <- list(sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL, timing = NULL) +# One-sided design {.tabset} -lower <- gs_spending_bound -lpar <- list(sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL, timing = NULL) -``` +In this example, we have efficacy bounds at both the IA and FA. The spending function used is the @lan1983discrete spending function with a total alpha of `r alpha`, which approximates an O'Brien-Fleming bound. -# Derive group sequential design +## Planed design -We assume there are totally 2 analyses, including FA. The IA happens at month 20, and FA occurs at month 36. ```{r} +upper <- gs_spending_bound +upar <- list(sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL) + x <- gs_design_ahr( enroll_rate = enroll_rate, fail_rate = fail_rate, @@ -66,9 +68,10 @@ x <- gs_design_ahr( ratio = ratio, upper = gs_spending_bound, upar = upar, - lower = lower, - lpar = lpar, - binding = FALSE) |> to_integer() + test_upper = TRUE, + lower = gs_b, + lpar = rep(-Inf, 2), + test_lower = FALSE) |> to_integer() ``` In the planed design, we have @@ -76,8 +79,9 @@ In the planed design, we have - planned events: `r round(x$analysis$event, 0)` - planned information fraction (timing): `r round(x$analysis$info_frac, 4)` - planned alpha spending: `r sfLDOF(0.025, x$analysis$info_frac)$spend` +- planned efficacy bounds: `r round(x$bound$z[x$bound$bound == "upper"], 4)` -```{r} +```{r, echo=FALSE} x |> summary() |> gt() |> @@ -88,7 +92,7 @@ x |> cell_text(weight = "bold") ), locations = cells_body( - columns = `Alternate hypothesis`, + columns = Z, rows = Bound == "Efficacy" )) |> tab_style( @@ -97,135 +101,216 @@ x |> cell_text(weight = "bold") ), locations = cells_body( - columns = `Alternate hypothesis`, + columns = Z, rows = Bound == "Futility" ) ) ``` -# Update bounds at time of analysis {.tabset} +## Update bounds at time of analysis -Assume in practice, there are 180 and 300 events observed in IA and FA, respectively. -```{r} +Assume in practice, there are 180 and 280 events observed in IA and FA, respectively. + +```{r, echo=FALSE} # observed vs. planned events -event_observed <- c(180, 300) +event_observed <- c(180, 280) event_planned <- x$analysis$event + +tibble(Analysis = c("IA", "FA"), + `Planned number of events` = event_planned, + `Observed number of events` = event_observed) |> + gt() |> + tab_header("Planned vs. Observed events") ``` -Now we calculate the updated timing based on the following 2 approaches. +Considering the observed number of events mentioned above, we utilize the `gs_power_npe()` function to generate updated efficacy bounds. As `gs_power_npe()` is a lower-level function that may be unfamiliar to some users, we will start by providing an introduction to its arguments as preliminaries. -## Event fraction based +- *3 theta-related arguments*: + + `theta`: natural parameter for group sequential design representing expected incremental drift at all analyses, i.e., $-\log(\text{AHR})$. It is used for crossing probability calculation. For example, if we set `theta = 0`, then the crossing probability is the type I error. Default is `0.1`. + + `theta0`: natural parameter for H0 and it decides the upper bounds. Default is `NULL`. + + `theta1`: natural parameter for H1 and it decides the lower bounds. Default is `NULL`. + +- *3 statistical information-related arguments*: + + `info`: statistical information at all analyses for input `theta`. Default is `1`. + + `info0`: statistical information under H0, if different than `info`. It impacts null hypothesis bound calculation. Default is `NULL`. + + `info1`: statistical information under hypothesis used for futility bound calculation if different from `info`. It impacts futility hypothesis bound calculation. Default is `NULL`. + +- *Efficacy/Futility boundary-related arguments*: + + `upper`, `upar` and `test_upper`: efficacy bounds related, same as `gs_design_ahr()`. + + `lower`, `lpar` and `test_lower`: futility bounds related, same as `gs_design_ahr()`. + + `binding`: `TRUE` or `FALSE` indicating whether it is binding or non-binding designs. + +When implementing `gs_power_npe()`, we initially set up `theta = 0`, which provides the crossing probability under the null hypothesis. Users have the flexibility to modify the value of `theta` if they are interested in the crossing probability under alternative hypotheses or other scenarios. In this implementation, we leave the setup of `theta0` with its default imputed value of `0`, indicating the null hypothesis. Additionally, as we are solely dealing with efficacy bounds, we do not adjust `theta1`, which influences the futility bound. -The updated timing is `r round(c(event_observed[1] / max(event_planned), 1), 4)`, if we divide the observed IA events by the planed FA events. We will spend alpha at IA based on its actual information fraction. The rest of the alpha will be all spent to FA. +Moving forward, we proceed to set up the values of the `info` arguments for the input `theta = 0`. Since the statistical information is event-based, and `theta = 0` (null hypothesis), the observed statistical information under null hypothesis is +$$ + \text{observed number of events} \times \frac{r}{1 + r} \times \frac{1}{1 + r}, +$$ +where $r$ is the randomization ratio (experimental : control). In this example, the observed statistical information is calculated as `info = event_observed / 4`, and thus, we input `r event_observed / 4` for the `info` argument. In this case, we retain the default imputed value for `info0`, which is set as `info`. The `info` represents the statistical information under the null hypothesis (`theta = 0`). Furthermore, we leave `info1` as it is, considering our focus on efficacy bounds, while noting that `info1` does affect the futility bounds. -```{r} -# updated spending time -upper_spending_time <- c(event_observed[1] / max(event_planned), 1) -lower_spending_time <- upper_spending_time - -# updated upar & lpar -upar_updated <- x$input$upar -upar_updated$timing <- upper_spending_time -lpar_updated <- x$input$lpar -lpar_updated$timing <- lower_spending_time -``` - -```{r} -# observed statistical information -info_observed <- gs_info_ahr( - enroll_rate = x$enroll_rate, - fail_rate = x$fail_rate, - ratio = ratio, - event = event_observed, - analysis_time = x$analysis$time) -``` +Lastly, we establish the parameters for the upper and lower bounds. The setup of `upper`, `upar`, `test_upper`, `lower`, `lpar`, and `test_lower` closely resembles the original design, with the exception of the `upar` setup. Here, we have introduced the `timing = ...` parameter. At the IA, the observed number of events is `r event_observed[1]`. If we employ the interim analysis alpha spending strategy, the spending timing at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. Here, `r round(event_planned[2], 0)` represents the planned number of events at FA. The remaining alpha will be allocated to the FA stage. ```{r} x_updated <- gs_power_npe( - theta = info_observed$theta, - theta0 = 0, - theta1 = x$analysis$theta, - info = info_observed$info, - info0 = x$analysis$info0, - info1 = x$analysis$info, - binding = x$input$binding, - upper = x$input$upper, - upar = upar_updated, - lower = x$input$lower, - lpar = lpar_updated, + # theta + theta = 0, + # info + info = event_observed / 4, + # upper bound + upper = gs_spending_bound, + upar = list(sf = gsDesign::sfLDOF, + total_spend = alpha, + param = NULL, + timing = c(event_observed[1] / max(event_planned), 1)), test_upper = TRUE, - test_lower = TRUE) - + # lower bound + lower = gs_b, + lpar = rep(-Inf, 2), + test_lower = FALSE, + # binding + binding = FALSE) +``` +So the updated efficacy bounds are `r x_updated$z[x_updated$bound == "upper"]`. +```{r, echo=FALSE} x_updated |> + filter(bound == "upper") |> gt() |> - tab_header(title = "Updated design", subtitle = "by event fraction based timing") |> + tab_header(title = "Updated design", + subtitle = paste0("with observed ", paste(event_observed, collapse = ", "), " events")) |> tab_style( style = list( cell_fill(color = "lightcyan"), cell_text(weight = "bold") ), locations = cells_body( - columns = probability, + columns = z, rows = bound == "upper" ) - ) |> + ) |> + tab_footnote(footnote = "Crossing pbability under H0.", locations = cells_column_labels(columns = probability)) +``` + +# Two-sided asymmetric design, beta-spending with non-binding lower bound {.tabset} + +In this section, we investigate a 2 sided asymmetric design, with the non-binding beta-spending futility bounds. Beta-spending refers to error spending for the lower bound crossing probabilities under the alternative hypothesis. Non-binding assumes the trial continues if the lower bound is crossed for Type I and Type II error computation. + +## Planed design + +In the original designs, we employ the Lan DeMets O'Brien-Fleming spending function [@lan1983discrete] for both efficacy and futility bounds. The total spending for efficacy is `r alpha`, and for futility is `r beta`. Besides, we assume the futility test only happens at IA. + +```{r} +upper <- gs_spending_bound +upar <- list(sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL) +lower <- gs_spending_bound +lpar <- list(sf = gsDesign::sfLDOF, total_spend = beta, param = NULL) + +x <- gs_design_ahr( + enroll_rate = enroll_rate, + fail_rate = fail_rate, + alpha = alpha, + beta = beta, + info_frac = NULL, + analysis_time = c(20, 36), + ratio = ratio, + upper = gs_spending_bound, + upar = upar, + test_upper = TRUE, + lower = lower, + lpar = lpar, + test_lower = c(TRUE, FALSE), + binding = FALSE) |> to_integer() +``` + +In the planed design, we have + +- planned events: `r round(x$analysis$event, 0)` +- planned information fraction (timing): `r round(x$analysis$info_frac, 4)` +- planned alpha spending: `r sfLDOF(0.025, x$analysis$info_frac)$spend` +- planned efficacy bounds: `r round(x$bound$z[x$bound$bound == "upper"], 4)` +- planned futility bounds: `r round(x$bound$z[x$bound$bound == "lower"], 4)` + +Since we added futility bounds, the sample size and number of events are larger than what we have in the 1-sided example. +```{r, echo=FALSE} +x |> + summary() |> + gt() |> + tab_header(title = "Original design") |> + tab_style( + style = list( + cell_fill(color = "lightcyan"), + cell_text(weight = "bold") + ), + locations = cells_body( + columns = Z, + rows = Bound == "Efficacy" + )) |> tab_style( style = list( cell_fill(color = "#F9E3D6"), cell_text(weight = "bold") ), locations = cells_body( - columns = probability, - rows = bound == "lower" + columns = Z, + rows = Bound == "Futility" ) ) ``` -## Information fraction based +## Update bounds at time of analysis -The updated timing is `r round(info_observed$info / max(info_observed$info), 4)`, which is the information fraction derived by the `gs_info_ahr()` by inputing the observed events. +In practice, let us assume that there were 180 events observed at IA stage and 280 events at the FA stage. -```{r} -# updated spending time -upper_spending_time <- info_observed$info / max(info_observed$info) -lower_spending_time <- upper_spending_time - -# updated upar & lpar -upar_updated <- x$input$upar -upar_updated$timing <- upper_spending_time -lpar_updated <- x$input$lpar -lpar_updated$timing <- lower_spending_time -``` +```{r, echo=FALSE} +# observed vs. planned events +event_observed <- c(180, 280) +event_planned <- x$analysis$event -```{r} -# observed statistical information -info_observed <- gs_info_ahr( - enroll_rate = x$enroll_rate, - fail_rate = x$fail_rate, - ratio = ratio, - event = event_observed, - analysis_time = x$analysis$time) +tibble(Analysis = c("IA", "FA"), + `Planned number of events` = event_planned, + `Observed number of events` = event_observed) |> + gt() |> + tab_header("Planned vs. Observed events") ``` +Again, we use `gs_power_npe()` to calculate the updated efficacy and futility bounds. We initially set up `theta = 0`, which provides the crossing probability under the null hypothesis. Users have the flexibility to modify the value of `theta` if they are interested in the crossing probability under alternative hypotheses or other scenarios. We set up `theta0 = 0` for null hypothesis. In this implementation, we set up `theta1` as the weighted average of the piecewise HR, i.e., +$$ + -\log(\text{AHR}) = -\log(\sum_{i=1}^m \text{HR}_m), +$$ +where the weight is decided by the observed number of events. For example, assume there are 50 events observed at the first 3 months (HR = 1), and 150 events observed after 3 months (HR = 0.6), we can derive `theta1` as $-\left[\frac{50}{280} \times \log(1) + \frac{230}{280} \times \log(0.6)\right]$. + +Moving forward, we proceed to set up the values of the `info` arguments for the input `theta = 0`. Similar to the examples in the 1-sided design above, we set `info` as `r event_observed`/4, where `r event_observed` is the observed number of events. Because `info0` is under null hypothesis, we leave it and use its default imput values same as `info`. Furthermore, we leave `info1` as its default imput values same as `info`. + +Lastly, we establish the parameters for the upper and lower bounds. The setup of `upper`, `upar`, `test_upper`, `lower`, `lpar`, and `test_lower` closely resembles the original design, with the exception of the `upar` and `lpar` setup. Here, we have introduced the `timing = ...` parameter. At the IA, the observed number of events is `r event_observed[1]`. If we employ the interim analysis alpha spending strategy, the spending timing at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. Here, `r round(event_planned[2], 0)` represents the planned number of events at FA. The remaining alpha will be allocated to the FA stage. + ```{r} x_updated <- gs_power_npe( - theta = info_observed$theta, + # theta + theta = 0, theta0 = 0, - theta1 = x$analysis$theta, - info = info_observed$info, - info0 = x$analysis$info0, - info1 = x$analysis$info, - binding = x$input$binding, + theta1 = -(50 * log(1) + 230 * log(0.6))/280, + # info, # ?? shall we set up info1? + info = event_observed / 4, + # upper bound upper = x$input$upper, - upar = upar_updated, - lower = x$input$lower, - lpar = lpar_updated, + upar = list(sf = gsDesign::sfLDOF, + total_spend = alpha, + param = NULL, + timing = c(event_observed[1] / max(event_planned), 1)), test_upper = TRUE, - test_lower = TRUE) + # lower bound + lower = x$input$lower, + lpar = list(sf = gsDesign::sfLDOF, + total_spend = beta, + param = NULL, + timing = c(event_observed[1] / max(event_planned), 1)), + test_lower = c(TRUE, FALSE), + binding = x$input$binding) +``` +```{r, echo=FALSE} x_updated |> gt() |> - tab_header(title = "Updated design", subtitle = "by information fraction based timing") |> + tab_header(title = "Updated design", subtitle = "by event fraction based timing") |> tab_style( style = list( cell_fill(color = "lightcyan"), @@ -235,7 +320,7 @@ x_updated |> columns = probability, rows = bound == "upper" ) - ) |> + ) |> tab_style( style = list( cell_fill(color = "#F9E3D6"), @@ -248,4 +333,9 @@ x_updated |> ) ``` + + + + + # References From ea720163bf8534fe4f75914da355496d588c6d21 Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Tue, 13 Feb 2024 16:26:40 -0500 Subject: [PATCH 03/30] add few more text --- vignettes/articles/story-update-boundary.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index bd49265a..e66fedda 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -278,7 +278,7 @@ $$ $$ where the weight is decided by the observed number of events. For example, assume there are 50 events observed at the first 3 months (HR = 1), and 150 events observed after 3 months (HR = 0.6), we can derive `theta1` as $-\left[\frac{50}{280} \times \log(1) + \frac{230}{280} \times \log(0.6)\right]$. -Moving forward, we proceed to set up the values of the `info` arguments for the input `theta = 0`. Similar to the examples in the 1-sided design above, we set `info` as `r event_observed`/4, where `r event_observed` is the observed number of events. Because `info0` is under null hypothesis, we leave it and use its default imput values same as `info`. Furthermore, we leave `info1` as its default imput values same as `info`. +Moving forward, we proceed to set up the values of the `info` arguments for the input `theta = 0`. Similar to the examples in the 1-sided design above, we set `info` as `r event_observed`/4, where `r event_observed` is the observed number of events. Because `info0` is under null hypothesis, we leave it and use its default imput values same as `info`. Furthermore, we leave `info1` as its default imput values same as `info` with the local null hypothesis approach. Though people can explictly write the formula of `info1` in terms of the observed events, but it would be easiest with unblinded data which we probably don't want to get into using. Lastly, we establish the parameters for the upper and lower bounds. The setup of `upper`, `upar`, `test_upper`, `lower`, `lpar`, and `test_lower` closely resembles the original design, with the exception of the `upar` and `lpar` setup. Here, we have introduced the `timing = ...` parameter. At the IA, the observed number of events is `r event_observed[1]`. If we employ the interim analysis alpha spending strategy, the spending timing at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. Here, `r round(event_planned[2], 0)` represents the planned number of events at FA. The remaining alpha will be allocated to the FA stage. From 2dd0f2eed19cad7cc5e39dc738d23cf6d6ee0c4f Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Tue, 13 Feb 2024 16:31:38 -0500 Subject: [PATCH 04/30] rewording and correct typo --- vignettes/articles/story-update-boundary.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index e66fedda..62e50991 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -278,7 +278,7 @@ $$ $$ where the weight is decided by the observed number of events. For example, assume there are 50 events observed at the first 3 months (HR = 1), and 150 events observed after 3 months (HR = 0.6), we can derive `theta1` as $-\left[\frac{50}{280} \times \log(1) + \frac{230}{280} \times \log(0.6)\right]$. -Moving forward, we proceed to set up the values of the `info` arguments for the input `theta = 0`. Similar to the examples in the 1-sided design above, we set `info` as `r event_observed`/4, where `r event_observed` is the observed number of events. Because `info0` is under null hypothesis, we leave it and use its default imput values same as `info`. Furthermore, we leave `info1` as its default imput values same as `info` with the local null hypothesis approach. Though people can explictly write the formula of `info1` in terms of the observed events, but it would be easiest with unblinded data which we probably don't want to get into using. +Moving forward, we proceed to set up the values of the `info` arguments for the input `theta = 0`. Similar to the examples in the 1-sided design above, we set `info` as `r event_observed`/4, where `r event_observed` is the observed number of events. Because `info0` is under null hypothesis, we leave it and use its default imputed values same as `info`. Furthermore, we leave `info1` as its default imputed values same as `info` with the local null hypothesis approach. Though people can explicitly write the formula of `info1` in terms of the observed events, the difference is very minor, and it would be easiest with unblinded data. Lastly, we establish the parameters for the upper and lower bounds. The setup of `upper`, `upar`, `test_upper`, `lower`, `lpar`, and `test_lower` closely resembles the original design, with the exception of the `upar` and `lpar` setup. Here, we have introduced the `timing = ...` parameter. At the IA, the observed number of events is `r event_observed[1]`. If we employ the interim analysis alpha spending strategy, the spending timing at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. Here, `r round(event_planned[2], 0)` represents the planned number of events at FA. The remaining alpha will be allocated to the FA stage. From 14d05cbebbfe4ec8a309f28ab56dc4690b856466 Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Tue, 13 Feb 2024 16:46:09 -0500 Subject: [PATCH 05/30] fix lintr issues --- vignettes/articles/story-update-boundary.Rmd | 66 ++++++++------------ 1 file changed, 26 insertions(+), 40 deletions(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 62e50991..311b7ff9 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -159,8 +159,8 @@ x_updated <- gs_power_npe( # upper bound upper = gs_spending_bound, upar = list(sf = gsDesign::sfLDOF, - total_spend = alpha, - param = NULL, + total_spend = alpha, + param = NULL, timing = c(event_observed[1] / max(event_planned), 1)), test_upper = TRUE, # lower bound @@ -172,21 +172,18 @@ x_updated <- gs_power_npe( ``` So the updated efficacy bounds are `r x_updated$z[x_updated$bound == "upper"]`. ```{r, echo=FALSE} -x_updated |> +x_updated |> filter(bound == "upper") |> gt() |> - tab_header(title = "Updated design", + tab_header(title = "Updated design", subtitle = paste0("with observed ", paste(event_observed, collapse = ", "), " events")) |> tab_style( style = list( cell_fill(color = "lightcyan"), - cell_text(weight = "bold") - ), + cell_text(weight = "bold")), locations = cells_body( columns = z, - rows = bound == "upper" - ) - ) |> + rows = bound == "upper")) |> tab_footnote(footnote = "Crossing pbability under H0.", locations = cells_column_labels(columns = probability)) ``` @@ -207,11 +204,11 @@ lpar <- list(sf = gsDesign::sfLDOF, total_spend = beta, param = NULL) x <- gs_design_ahr( enroll_rate = enroll_rate, fail_rate = fail_rate, - alpha = alpha, + alpha = alpha, beta = beta, - info_frac = NULL, + info_frac = NULL, analysis_time = c(20, 36), - ratio = ratio, + ratio = ratio, upper = gs_spending_bound, upar = upar, test_upper = TRUE, @@ -238,22 +235,17 @@ x |> tab_style( style = list( cell_fill(color = "lightcyan"), - cell_text(weight = "bold") - ), + cell_text(weight = "bold")), locations = cells_body( columns = Z, - rows = Bound == "Efficacy" - )) |> + rows = Bound == "Efficacy")) |> tab_style( style = list( cell_fill(color = "#F9E3D6"), - cell_text(weight = "bold") - ), + cell_text(weight = "bold")), locations = cells_body( columns = Z, - rows = Bound == "Futility" - ) - ) + rows = Bound == "Futility")) ``` ## Update bounds at time of analysis @@ -267,7 +259,7 @@ event_planned <- x$analysis$event tibble(Analysis = c("IA", "FA"), `Planned number of events` = event_planned, - `Observed number of events` = event_observed) |> + `Observed number of events` = event_observed) |> gt() |> tab_header("Planned vs. Observed events") ``` @@ -289,48 +281,42 @@ x_updated <- gs_power_npe( theta0 = 0, theta1 = -(50 * log(1) + 230 * log(0.6))/280, # info, # ?? shall we set up info1? - info = event_observed / 4, + info = event_observed / 4, # upper bound upper = x$input$upper, - upar = list(sf = gsDesign::sfLDOF, - total_spend = alpha, - param = NULL, + upar = list(sf = gsDesign::sfLDOF, + total_spend = alpha, + param = NULL, timing = c(event_observed[1] / max(event_planned), 1)), test_upper = TRUE, # lower bound lower = x$input$lower, - lpar = list(sf = gsDesign::sfLDOF, - total_spend = beta, - param = NULL, + lpar = list(sf = gsDesign::sfLDOF, + total_spend = beta, + param = NULL, timing = c(event_observed[1] / max(event_planned), 1)), test_lower = c(TRUE, FALSE), binding = x$input$binding) ``` ```{r, echo=FALSE} -x_updated |> +x_updated |> gt() |> tab_header(title = "Updated design", subtitle = "by event fraction based timing") |> tab_style( style = list( cell_fill(color = "lightcyan"), - cell_text(weight = "bold") - ), + cell_text(weight = "bold")), locations = cells_body( columns = probability, - rows = bound == "upper" - ) - ) |> + rows = bound == "upper")) |> tab_style( style = list( cell_fill(color = "#F9E3D6"), - cell_text(weight = "bold") - ), + cell_text(weight = "bold")), locations = cells_body( columns = probability, - rows = bound == "lower" - ) - ) + rows = bound == "lower")) ``` From 35ba0ff54451fb476a7801235c2049970f2bce22 Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Tue, 13 Feb 2024 16:47:27 -0500 Subject: [PATCH 06/30] update pkgdown.yml --- _pkgdown.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/_pkgdown.yml b/_pkgdown.yml index db61dae5..3c409c34 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -139,6 +139,7 @@ articles: - articles/story-npe-integration - articles/story-nph-futility - articles/story-arbitrary-distribution + - articles/story-update-boundary - title: "Spending bound" desc: > From 21d7a36f3e93cd229de7ff06dc487a409a1fc942 Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Tue, 13 Feb 2024 16:58:34 -0500 Subject: [PATCH 07/30] fix lintr issues --- vignettes/articles/story-update-boundary.Rmd | 236 ++++++++----------- 1 file changed, 100 insertions(+), 136 deletions(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 311b7ff9..3365df62 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -58,20 +58,19 @@ In this example, we have efficacy bounds at both the IA and FA. The spending fun upper <- gs_spending_bound upar <- list(sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL) -x <- gs_design_ahr( - enroll_rate = enroll_rate, - fail_rate = fail_rate, - alpha = alpha, - beta = beta, - info_frac = NULL, - analysis_time = c(20, 36), - ratio = ratio, - upper = gs_spending_bound, - upar = upar, - test_upper = TRUE, - lower = gs_b, - lpar = rep(-Inf, 2), - test_lower = FALSE) |> to_integer() +x <- gs_design_ahr(enroll_rate = enroll_rate, + fail_rate = fail_rate, + alpha = alpha, + beta = beta, + info_frac = NULL, + analysis_time = c(20, 36), + ratio = ratio, + upper = gs_spending_bound, + upar = upar, + test_upper = TRUE, + lower = gs_b, + lpar = rep(-Inf, 2), + test_lower = FALSE) |> to_integer() ``` In the planed design, we have @@ -82,29 +81,18 @@ In the planed design, we have - planned efficacy bounds: `r round(x$bound$z[x$bound$bound == "upper"], 4)` ```{r, echo=FALSE} -x |> - summary() |> +x |> + summary() |> gt() |> tab_header(title = "Original design") |> - tab_style( - style = list( - cell_fill(color = "lightcyan"), - cell_text(weight = "bold") - ), - locations = cells_body( - columns = Z, - rows = Bound == "Efficacy" - )) |> - tab_style( - style = list( - cell_fill(color = "#F9E3D6"), - cell_text(weight = "bold") - ), - locations = cells_body( - columns = Z, - rows = Bound == "Futility" - ) - ) + tab_style(style = list(cell_fill(color = "lightcyan"), + cell_text(weight = "bold")), + locations = cells_body(columns = Z, + rows = Bound == "Efficacy")) |> + tab_style(style = list(cell_fill(color = "#F9E3D6"), + cell_text(weight = "bold")), + locations = cells_body(columns = Z, + rows = Bound == "Futility")) ``` ## Update bounds at time of analysis @@ -119,7 +107,7 @@ event_planned <- x$analysis$event tibble(Analysis = c("IA", "FA"), `Planned number of events` = event_planned, `Observed number of events` = event_observed) |> - gt() |> + gt() |> tab_header("Planned vs. Observed events") ``` @@ -151,24 +139,22 @@ where $r$ is the randomization ratio (experimental : control). In this example, Lastly, we establish the parameters for the upper and lower bounds. The setup of `upper`, `upar`, `test_upper`, `lower`, `lpar`, and `test_lower` closely resembles the original design, with the exception of the `upar` setup. Here, we have introduced the `timing = ...` parameter. At the IA, the observed number of events is `r event_observed[1]`. If we employ the interim analysis alpha spending strategy, the spending timing at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. Here, `r round(event_planned[2], 0)` represents the planned number of events at FA. The remaining alpha will be allocated to the FA stage. ```{r} -x_updated <- gs_power_npe( - # theta - theta = 0, - # info - info = event_observed / 4, - # upper bound - upper = gs_spending_bound, - upar = list(sf = gsDesign::sfLDOF, - total_spend = alpha, - param = NULL, - timing = c(event_observed[1] / max(event_planned), 1)), - test_upper = TRUE, - # lower bound - lower = gs_b, - lpar = rep(-Inf, 2), - test_lower = FALSE, - # binding - binding = FALSE) +x_updated <- gs_power_npe(theta = 0, + # info + info = event_observed / 4, + # upper bound + upper = gs_spending_bound, + upar = list(sf = gsDesign::sfLDOF, + total_spend = alpha, + param = NULL, + timing = c(event_observed[1] / max(event_planned), 1)), + test_upper = TRUE, + # lower bound + lower = gs_b, + lpar = rep(-Inf, 2), + test_lower = FALSE, + # binding + binding = FALSE) ``` So the updated efficacy bounds are `r x_updated$z[x_updated$bound == "upper"]`. ```{r, echo=FALSE} @@ -177,14 +163,12 @@ x_updated |> gt() |> tab_header(title = "Updated design", subtitle = paste0("with observed ", paste(event_observed, collapse = ", "), " events")) |> - tab_style( - style = list( - cell_fill(color = "lightcyan"), - cell_text(weight = "bold")), - locations = cells_body( - columns = z, - rows = bound == "upper")) |> - tab_footnote(footnote = "Crossing pbability under H0.", locations = cells_column_labels(columns = probability)) + tab_style(style = list(cell_fill(color = "lightcyan"), + cell_text(weight = "bold")), + locations = cells_body(columns = z, + rows = bound == "upper")) |> + tab_footnote(footnote = "Crossing pbability under H0.", + locations = cells_column_labels(columns = probability)) ``` # Two-sided asymmetric design, beta-spending with non-binding lower bound {.tabset} @@ -199,23 +183,22 @@ In the original designs, we employ the Lan DeMets O'Brien-Fleming spending funct upper <- gs_spending_bound upar <- list(sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL) lower <- gs_spending_bound -lpar <- list(sf = gsDesign::sfLDOF, total_spend = beta, param = NULL) - -x <- gs_design_ahr( - enroll_rate = enroll_rate, - fail_rate = fail_rate, - alpha = alpha, - beta = beta, - info_frac = NULL, - analysis_time = c(20, 36), - ratio = ratio, - upper = gs_spending_bound, - upar = upar, - test_upper = TRUE, - lower = lower, - lpar = lpar, - test_lower = c(TRUE, FALSE), - binding = FALSE) |> to_integer() +lpar <- list(sf = gsDesign::sfLDOF, total_spend = beta, param = NULL) + +x <- gs_design_ahr(enroll_rate = enroll_rate, + fail_rate = fail_rate, + alpha = alpha, + beta = beta, + info_frac = NULL, + analysis_time = c(20, 36), + ratio = ratio, + upper = gs_spending_bound, + upar = upar, + test_upper = TRUE, + lower = lower, + lpar = lpar, + test_lower = c(TRUE, FALSE), + binding = FALSE) |> to_integer() ``` In the planed design, we have @@ -228,24 +211,18 @@ In the planed design, we have Since we added futility bounds, the sample size and number of events are larger than what we have in the 1-sided example. ```{r, echo=FALSE} -x |> - summary() |> +x |> + summary() |> gt() |> tab_header(title = "Original design") |> - tab_style( - style = list( - cell_fill(color = "lightcyan"), - cell_text(weight = "bold")), - locations = cells_body( - columns = Z, - rows = Bound == "Efficacy")) |> - tab_style( - style = list( - cell_fill(color = "#F9E3D6"), - cell_text(weight = "bold")), - locations = cells_body( - columns = Z, - rows = Bound == "Futility")) + tab_style(style = list(cell_fill(color = "lightcyan"), + cell_text(weight = "bold")), + locations = cells_body(columns = Z, + rows = Bound == "Efficacy")) |> + tab_style(style = list(cell_fill(color = "#F9E3D6"), + cell_text(weight = "bold")), + locations = cells_body(columns = Z, + rows = Bound == "Futility")) ``` ## Update bounds at time of analysis @@ -275,53 +252,40 @@ Moving forward, we proceed to set up the values of the `info` arguments for the Lastly, we establish the parameters for the upper and lower bounds. The setup of `upper`, `upar`, `test_upper`, `lower`, `lpar`, and `test_lower` closely resembles the original design, with the exception of the `upar` and `lpar` setup. Here, we have introduced the `timing = ...` parameter. At the IA, the observed number of events is `r event_observed[1]`. If we employ the interim analysis alpha spending strategy, the spending timing at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. Here, `r round(event_planned[2], 0)` represents the planned number of events at FA. The remaining alpha will be allocated to the FA stage. ```{r} -x_updated <- gs_power_npe( - # theta - theta = 0, - theta0 = 0, - theta1 = -(50 * log(1) + 230 * log(0.6))/280, - # info, # ?? shall we set up info1? - info = event_observed / 4, - # upper bound - upper = x$input$upper, - upar = list(sf = gsDesign::sfLDOF, - total_spend = alpha, - param = NULL, - timing = c(event_observed[1] / max(event_planned), 1)), - test_upper = TRUE, - # lower bound - lower = x$input$lower, - lpar = list(sf = gsDesign::sfLDOF, - total_spend = beta, - param = NULL, - timing = c(event_observed[1] / max(event_planned), 1)), - test_lower = c(TRUE, FALSE), - binding = x$input$binding) +x_updated <- gs_power_npe(theta = 0, + theta0 = 0, + theta1 = -(50 * log(1) + 230 * log(0.6))/280, + # info + info = event_observed / 4, + # upper bound + upper = x$input$upper, + upar = list(sf = gsDesign::sfLDOF, + total_spend = alpha, + param = NULL, + timing = c(event_observed[1] / max(event_planned), 1)), + test_upper = TRUE, + # lower bound + lower = x$input$lower, + lpar = list(sf = gsDesign::sfLDOF, + total_spend = beta, + param = NULL, + timing = c(event_observed[1] / max(event_planned), 1)), + test_lower = c(TRUE, FALSE), + binding = x$input$binding) ``` ```{r, echo=FALSE} x_updated |> gt() |> tab_header(title = "Updated design", subtitle = "by event fraction based timing") |> - tab_style( - style = list( - cell_fill(color = "lightcyan"), - cell_text(weight = "bold")), - locations = cells_body( - columns = probability, - rows = bound == "upper")) |> - tab_style( - style = list( - cell_fill(color = "#F9E3D6"), - cell_text(weight = "bold")), - locations = cells_body( - columns = probability, - rows = bound == "lower")) + tab_style(style = list(cell_fill(color = "lightcyan"), + cell_text(weight = "bold")), + locations = cells_body(columns = probability, + rows = bound == "upper")) |> + tab_style(style = list(cell_fill(color = "#F9E3D6"), + cell_text(weight = "bold")), + locations = cells_body(columns = probability, + rows = bound == "lower")) ``` - - - - - # References From 4cd384dbb7ee96023e5088387d60a8f66da8acfc Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Wed, 14 Feb 2024 11:03:20 -0500 Subject: [PATCH 08/30] fix lintr --- vignettes/articles/story-update-boundary.Rmd | 22 +++++++++----------- 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 3365df62..65f021a7 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -36,16 +36,14 @@ beta <- 0.1 ratio <- 1 # enrollment -enroll_rate <- define_enroll_rate( - duration = c(2, 2, 10), - rate = c(3, 6, 9)) +enroll_rate <- define_enroll_rate(duration = c(2, 2, 10), + rate = c(3, 6, 9)) # failure and dropout -fail_rate <- define_fail_rate( - duration = c(3, 100), - fail_rate = log(2) / 9, - hr = c(1, 0.6), - dropout_rate = .0001) +fail_rate <- define_fail_rate(duration = c(3, 100), + fail_rate = log(2) / 9, + hr = c(1, 0.6), + dropout_rate = .0001) ``` # One-sided design {.tabset} @@ -106,7 +104,7 @@ event_planned <- x$analysis$event tibble(Analysis = c("IA", "FA"), `Planned number of events` = event_planned, - `Observed number of events` = event_observed) |> + `Observed number of events` = event_observed) |> gt() |> tab_header("Planned vs. Observed events") ``` @@ -141,7 +139,7 @@ Lastly, we establish the parameters for the upper and lower bounds. The setup of ```{r} x_updated <- gs_power_npe(theta = 0, # info - info = event_observed / 4, + info = event_observed / 4, # upper bound upper = gs_spending_bound, upar = list(sf = gsDesign::sfLDOF, @@ -167,7 +165,7 @@ x_updated |> cell_text(weight = "bold")), locations = cells_body(columns = z, rows = bound == "upper")) |> - tab_footnote(footnote = "Crossing pbability under H0.", + tab_footnote(footnote = "Crossing pbability under H0.", locations = cells_column_labels(columns = probability)) ``` @@ -254,7 +252,7 @@ Lastly, we establish the parameters for the upper and lower bounds. The setup of ```{r} x_updated <- gs_power_npe(theta = 0, theta0 = 0, - theta1 = -(50 * log(1) + 230 * log(0.6))/280, + theta1 = -(50 * log(1) + 230 * log(0.6)) / 280, # info info = event_observed / 4, # upper bound From 6afa9866fe385e96bb78e61fa7a2bb9837519907 Mon Sep 17 00:00:00 2001 From: Nan Xiao Date: Wed, 14 Feb 2024 19:19:01 -0500 Subject: [PATCH 09/30] Improve code style for boundary update vignette --- vignettes/articles/story-update-boundary.Rmd | 515 ++++++++++++------- 1 file changed, 333 insertions(+), 182 deletions(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 65f021a7..521aebea 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -1,5 +1,5 @@ --- -title: "Efficacy and Futility Boundary Update" +title: "Efficacy and futility boundary update" author: "Yujie Zhao and Keaven M. Anderson" output: rmarkdown::html_document: @@ -11,279 +11,430 @@ output: css: "custom.css" bibliography: "gsDesign2.bib" vignette: > - %\VignetteIndexEntry{Efficacy and Futility Boundary Update} - %\VignetteEncoding{UTF-8} + %\VignetteIndexEntry{Efficacy and futility boundary update} %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} -knitr::opts_chunk$set(echo = TRUE) +knitr::opts_chunk$set(echo = FALSE) +``` -library(gsDesign) +```{r, message=FALSE, warning=FALSE} library(gsDesign2) -library(gt) -library(tibble) -library(dplyr) ``` # Design assumptions -We assume a total of two analyses, including a final Analysis (FA). The initial Analysis (IA) occurs at month 20, followed by the FA at month 36. The enrollment period spans 14 months, with the first 2 months having an enrollment rate of 3, the next 2 months with a rate of 6, and the final 10 months with a rate of 9. The planned enrollment rate to achieve the 90\% power is a multiptier of 3, 6, and 9. Additionally, there is a delayed treatment effect in the first 3 months (HR = 1), which transitions to an HR of 0.6 thereafter. The median of the control arm is 9 months, and the dropout rate is 0.0001. +We assume a total of two analyses, including a final Analysis (FA). +The initial Analysis (IA) occurs at month 20, followed by the FA at month 36. +The enrollment period spans 14 months, with the first 2 months having an +enrollment rate of 3, the next 2 months with a rate of 6, and the final +10 months with a rate of 9. +The planned enrollment rate to achieve the 90\% power is a multiplier of +3, 6, and 9. Additionally, there is a delayed treatment effect in the first +3 months (HR = 1), which transitions to an HR of 0.6 thereafter. +The median of the control arm is 9 months, and the dropout rate is 0.0001. -```{r, echo=FALSE} +```{r} alpha <- 0.025 beta <- 0.1 ratio <- 1 -# enrollment -enroll_rate <- define_enroll_rate(duration = c(2, 2, 10), - rate = c(3, 6, 9)) - -# failure and dropout -fail_rate <- define_fail_rate(duration = c(3, 100), - fail_rate = log(2) / 9, - hr = c(1, 0.6), - dropout_rate = .0001) +# Enrollment +enroll_rate <- define_enroll_rate( + duration = c(2, 2, 10), + rate = c(3, 6, 9) +) + +# Failure and dropout +fail_rate <- define_fail_rate( + duration = c(3, 100), + fail_rate = log(2) / 9, + hr = c(1, 0.6), + dropout_rate = .0001 +) ``` # One-sided design {.tabset} -In this example, we have efficacy bounds at both the IA and FA. The spending function used is the @lan1983discrete spending function with a total alpha of `r alpha`, which approximates an O'Brien-Fleming bound. +In this example, we have efficacy bounds at both the IA and FA. The spending +function used is the @lan1983discrete spending function with a total alpha of +`r alpha`, which approximates an O'Brien-Fleming bound. ## Planed design -```{r} +```{r, echo=TRUE} upper <- gs_spending_bound upar <- list(sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL) -x <- gs_design_ahr(enroll_rate = enroll_rate, - fail_rate = fail_rate, - alpha = alpha, - beta = beta, - info_frac = NULL, - analysis_time = c(20, 36), - ratio = ratio, - upper = gs_spending_bound, - upar = upar, - test_upper = TRUE, - lower = gs_b, - lpar = rep(-Inf, 2), - test_lower = FALSE) |> to_integer() +x <- gs_design_ahr( + enroll_rate = enroll_rate, + fail_rate = fail_rate, + alpha = alpha, + beta = beta, + info_frac = NULL, + analysis_time = c(20, 36), + ratio = ratio, + upper = gs_spending_bound, + upar = upar, + test_upper = TRUE, + lower = gs_b, + lpar = rep(-Inf, 2), + test_lower = FALSE +) |> to_integer() ``` -In the planed design, we have +In the planed design, we have - planned events: `r round(x$analysis$event, 0)` - planned information fraction (timing): `r round(x$analysis$info_frac, 4)` -- planned alpha spending: `r sfLDOF(0.025, x$analysis$info_frac)$spend` +- planned alpha spending: `r gsDesign::sfLDOF(0.025, x$analysis$info_frac)$spend` - planned efficacy bounds: `r round(x$bound$z[x$bound$bound == "upper"], 4)` -```{r, echo=FALSE} +```{r} x |> summary() |> - gt() |> - tab_header(title = "Original design") |> - tab_style(style = list(cell_fill(color = "lightcyan"), - cell_text(weight = "bold")), - locations = cells_body(columns = Z, - rows = Bound == "Efficacy")) |> - tab_style(style = list(cell_fill(color = "#F9E3D6"), - cell_text(weight = "bold")), - locations = cells_body(columns = Z, - rows = Bound == "Futility")) + gt::gt() |> + gt::tab_header(title = "Original design") |> + gt::tab_style( + style = list( + gt::cell_fill(color = "lightcyan"), + gt::cell_text(weight = "bold") + ), + locations = gt::cells_body( + columns = Z, + rows = Bound == "Efficacy" + ) + ) |> + gt::tab_style( + style = list( + gt::cell_fill(color = "#F9E3D6"), + gt::cell_text(weight = "bold") + ), + locations = gt::cells_body( + columns = Z, + rows = Bound == "Futility" + ) + ) ``` -## Update bounds at time of analysis +## Update bounds at time of analysis -Assume in practice, there are 180 and 280 events observed in IA and FA, respectively. +Assume in practice, there are 180 and 280 events observed in IA and FA, respectively. -```{r, echo=FALSE} -# observed vs. planned events +```{r} +# Observed vs. planned events event_observed <- c(180, 280) event_planned <- x$analysis$event -tibble(Analysis = c("IA", "FA"), - `Planned number of events` = event_planned, - `Observed number of events` = event_observed) |> - gt() |> - tab_header("Planned vs. Observed events") +tibble::tibble( + Analysis = c("IA", "FA"), + `Planned number of events` = event_planned, + `Observed number of events` = event_observed +) |> + gt::gt() |> + gt::tab_header("Planned vs. Observed events") ``` -Considering the observed number of events mentioned above, we utilize the `gs_power_npe()` function to generate updated efficacy bounds. As `gs_power_npe()` is a lower-level function that may be unfamiliar to some users, we will start by providing an introduction to its arguments as preliminaries. +Considering the observed number of events mentioned above, we utilize the +`gs_power_npe()` function to generate updated efficacy bounds. +As `gs_power_npe()` is a lower-level function that may be unfamiliar to some +users, we will start by providing an introduction to its arguments as preliminaries. - *3 theta-related arguments*: - + `theta`: natural parameter for group sequential design representing expected incremental drift at all analyses, i.e., $-\log(\text{AHR})$. It is used for crossing probability calculation. For example, if we set `theta = 0`, then the crossing probability is the type I error. Default is `0.1`. - + `theta0`: natural parameter for H0 and it decides the upper bounds. Default is `NULL`. - + `theta1`: natural parameter for H1 and it decides the lower bounds. Default is `NULL`. - + - `theta`: natural parameter for group sequential design representing expected + incremental drift at all analyses, i.e., $-\log(\text{AHR})$. + It is used for crossing probability calculation. + For example, if we set `theta = 0`, then the crossing probability is the + type I error. Default is `0.1`. + - `theta0`: natural parameter for H0 and it decides the upper bounds. + Default is `NULL`. + - `theta1`: natural parameter for H1 and it decides the lower bounds. + Default is `NULL`. + - *3 statistical information-related arguments*: - + `info`: statistical information at all analyses for input `theta`. Default is `1`. - + `info0`: statistical information under H0, if different than `info`. It impacts null hypothesis bound calculation. Default is `NULL`. - + `info1`: statistical information under hypothesis used for futility bound calculation if different from `info`. It impacts futility hypothesis bound calculation. Default is `NULL`. - -- *Efficacy/Futility boundary-related arguments*: - + `upper`, `upar` and `test_upper`: efficacy bounds related, same as `gs_design_ahr()`. - + `lower`, `lpar` and `test_lower`: futility bounds related, same as `gs_design_ahr()`. - + `binding`: `TRUE` or `FALSE` indicating whether it is binding or non-binding designs. - -When implementing `gs_power_npe()`, we initially set up `theta = 0`, which provides the crossing probability under the null hypothesis. Users have the flexibility to modify the value of `theta` if they are interested in the crossing probability under alternative hypotheses or other scenarios. In this implementation, we leave the setup of `theta0` with its default imputed value of `0`, indicating the null hypothesis. Additionally, as we are solely dealing with efficacy bounds, we do not adjust `theta1`, which influences the futility bound. + - `info`: statistical information at all analyses for input `theta`. + Default is `1`. + - `info0`: statistical information under H0, if different than `info`. + It impacts null hypothesis bound calculation. Default is `NULL`. + - `info1`: statistical information under hypothesis used for futility bound + calculation if different from `info`. + It impacts futility hypothesis bound calculation. Default is `NULL`. -Moving forward, we proceed to set up the values of the `info` arguments for the input `theta = 0`. Since the statistical information is event-based, and `theta = 0` (null hypothesis), the observed statistical information under null hypothesis is +- *Efficacy/Futility boundary-related arguments*: + - `upper`, `upar` and `test_upper`: efficacy bounds related, same as `gs_design_ahr()`. + - `lower`, `lpar` and `test_lower`: futility bounds related, same as `gs_design_ahr()`. + - `binding`: `TRUE` or `FALSE` indicating whether it is binding or non-binding designs. + +When implementing `gs_power_npe()`, we initially set up `theta = 0`, +which provides the crossing probability under the null hypothesis. +Users have the flexibility to modify the value of `theta` if they are +interested in the crossing probability under alternative hypotheses or +other scenarios. In this implementation, we leave the setup of `theta0` +with its default imputed value of `0`, indicating the null hypothesis. +Additionally, as we are solely dealing with efficacy bounds, +we do not adjust `theta1`, which influences the futility bound. + +Moving forward, we proceed to set up the values of the `info` arguments +for the input `theta = 0`. Since the statistical information is event-based, +and `theta = 0` (null hypothesis), the observed statistical information +under null hypothesis is $$ \text{observed number of events} \times \frac{r}{1 + r} \times \frac{1}{1 + r}, $$ -where $r$ is the randomization ratio (experimental : control). In this example, the observed statistical information is calculated as `info = event_observed / 4`, and thus, we input `r event_observed / 4` for the `info` argument. In this case, we retain the default imputed value for `info0`, which is set as `info`. The `info` represents the statistical information under the null hypothesis (`theta = 0`). Furthermore, we leave `info1` as it is, considering our focus on efficacy bounds, while noting that `info1` does affect the futility bounds. +where $r$ is the randomization ratio (experimental : control). +In this example, the observed statistical information is calculated as +`info = event_observed / 4`, and thus, we input `r event_observed / 4` +for the `info` argument. In this case, we retain the default imputed value +for `info0`, which is set as `info`. The `info` represents the +statistical information under the null hypothesis (`theta = 0`). +Furthermore, we leave `info1` as it is, considering our focus on efficacy bounds, +while noting that `info1` does affect the futility bounds. + +Lastly, we establish the parameters for the upper and lower bounds. +The setup of `upper`, `upar`, `test_upper`, `lower`, `lpar`, and `test_lower` +closely resembles the original design, with the exception of the `upar` setup. +Here, we have introduced the `timing = ...` parameter. +At the IA, the observed number of events is `r event_observed[1]`. +If we employ the interim analysis alpha spending strategy, the spending timing +at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. +Here, `r round(event_planned[2], 0)` represents the planned number of events +at FA. The remaining alpha will be allocated to the FA stage. + +```{r, echo=TRUE} +x_updated <- gs_power_npe( + theta = 0, + # Info + info = event_observed / 4, + # Upper bound + upper = gs_spending_bound, + upar = list( + sf = gsDesign::sfLDOF, + total_spend = alpha, + param = NULL, + timing = c(event_observed[1] / max(event_planned), 1) + ), + test_upper = TRUE, + # Lower bound + lower = gs_b, + lpar = rep(-Inf, 2), + test_lower = FALSE, + # Binding + binding = FALSE +) +``` -Lastly, we establish the parameters for the upper and lower bounds. The setup of `upper`, `upar`, `test_upper`, `lower`, `lpar`, and `test_lower` closely resembles the original design, with the exception of the `upar` setup. Here, we have introduced the `timing = ...` parameter. At the IA, the observed number of events is `r event_observed[1]`. If we employ the interim analysis alpha spending strategy, the spending timing at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. Here, `r round(event_planned[2], 0)` represents the planned number of events at FA. The remaining alpha will be allocated to the FA stage. +So the updated efficacy bounds are `r x_updated$z[x_updated$bound == "upper"]`. ```{r} -x_updated <- gs_power_npe(theta = 0, - # info - info = event_observed / 4, - # upper bound - upper = gs_spending_bound, - upar = list(sf = gsDesign::sfLDOF, - total_spend = alpha, - param = NULL, - timing = c(event_observed[1] / max(event_planned), 1)), - test_upper = TRUE, - # lower bound - lower = gs_b, - lpar = rep(-Inf, 2), - test_lower = FALSE, - # binding - binding = FALSE) -``` -So the updated efficacy bounds are `r x_updated$z[x_updated$bound == "upper"]`. -```{r, echo=FALSE} x_updated |> - filter(bound == "upper") |> - gt() |> - tab_header(title = "Updated design", - subtitle = paste0("with observed ", paste(event_observed, collapse = ", "), " events")) |> - tab_style(style = list(cell_fill(color = "lightcyan"), - cell_text(weight = "bold")), - locations = cells_body(columns = z, - rows = bound == "upper")) |> - tab_footnote(footnote = "Crossing pbability under H0.", - locations = cells_column_labels(columns = probability)) + dplyr::filter(bound == "upper") |> + gt::gt() |> + gt::tab_header( + title = "Updated design", + subtitle = paste0("with observed ", paste(event_observed, collapse = ", "), " events") + ) |> + gt::tab_style( + style = list( + gt::cell_fill(color = "lightcyan"), + gt::cell_text(weight = "bold") + ), + locations = gt::cells_body( + columns = z, + rows = bound == "upper" + ) + ) |> + gt::tab_footnote( + footnote = "Crossing pbability under H0.", + locations = gt::cells_column_labels(columns = probability) + ) ``` # Two-sided asymmetric design, beta-spending with non-binding lower bound {.tabset} -In this section, we investigate a 2 sided asymmetric design, with the non-binding beta-spending futility bounds. Beta-spending refers to error spending for the lower bound crossing probabilities under the alternative hypothesis. Non-binding assumes the trial continues if the lower bound is crossed for Type I and Type II error computation. +In this section, we investigate a 2 sided asymmetric design, with the +non-binding beta-spending futility bounds. Beta-spending refers to +error spending for the lower bound crossing probabilities under the +alternative hypothesis. Non-binding assumes the trial continues if the +lower bound is crossed for Type I and Type II error computation. ## Planed design -In the original designs, we employ the Lan DeMets O'Brien-Fleming spending function [@lan1983discrete] for both efficacy and futility bounds. The total spending for efficacy is `r alpha`, and for futility is `r beta`. Besides, we assume the futility test only happens at IA. +In the original designs, we employ the Lan DeMets O'Brien-Fleming +spending function [@lan1983discrete] for both efficacy and futility bounds. +The total spending for efficacy is `r alpha`, and for futility is `r beta`. +Besides, we assume the futility test only happens at IA. -```{r} +```{r, echo=TRUE} upper <- gs_spending_bound upar <- list(sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL) lower <- gs_spending_bound lpar <- list(sf = gsDesign::sfLDOF, total_spend = beta, param = NULL) -x <- gs_design_ahr(enroll_rate = enroll_rate, - fail_rate = fail_rate, - alpha = alpha, - beta = beta, - info_frac = NULL, - analysis_time = c(20, 36), - ratio = ratio, - upper = gs_spending_bound, - upar = upar, - test_upper = TRUE, - lower = lower, - lpar = lpar, - test_lower = c(TRUE, FALSE), - binding = FALSE) |> to_integer() +x <- gs_design_ahr( + enroll_rate = enroll_rate, + fail_rate = fail_rate, + alpha = alpha, + beta = beta, + info_frac = NULL, + analysis_time = c(20, 36), + ratio = ratio, + upper = gs_spending_bound, + upar = upar, + test_upper = TRUE, + lower = lower, + lpar = lpar, + test_lower = c(TRUE, FALSE), + binding = FALSE +) |> to_integer() ``` -In the planed design, we have +In the planed design, we have - planned events: `r round(x$analysis$event, 0)` - planned information fraction (timing): `r round(x$analysis$info_frac, 4)` -- planned alpha spending: `r sfLDOF(0.025, x$analysis$info_frac)$spend` +- planned alpha spending: `r gsDesign::sfLDOF(0.025, x$analysis$info_frac)$spend` - planned efficacy bounds: `r round(x$bound$z[x$bound$bound == "upper"], 4)` - planned futility bounds: `r round(x$bound$z[x$bound$bound == "lower"], 4)` -Since we added futility bounds, the sample size and number of events are larger than what we have in the 1-sided example. -```{r, echo=FALSE} +Since we added futility bounds, the sample size and number of events are +larger than what we have in the 1-sided example. + +```{r} x |> summary() |> - gt() |> - tab_header(title = "Original design") |> - tab_style(style = list(cell_fill(color = "lightcyan"), - cell_text(weight = "bold")), - locations = cells_body(columns = Z, - rows = Bound == "Efficacy")) |> - tab_style(style = list(cell_fill(color = "#F9E3D6"), - cell_text(weight = "bold")), - locations = cells_body(columns = Z, - rows = Bound == "Futility")) + gt::gt() |> + gt::tab_header(title = "Original design") |> + gt::tab_style( + style = list( + gt::cell_fill(color = "lightcyan"), + gt::cell_text(weight = "bold") + ), + locations = gt::cells_body( + columns = Z, + rows = Bound == "Efficacy" + ) + ) |> + gt::tab_style( + style = list( + gt::cell_fill(color = "#F9E3D6"), + gt::cell_text(weight = "bold") + ), + locations = gt::cells_body( + columns = Z, + rows = Bound == "Futility" + ) + ) ``` -## Update bounds at time of analysis +## Update bounds at time of analysis -In practice, let us assume that there were 180 events observed at IA stage and 280 events at the FA stage. +In practice, let us assume that there were 180 events observed at IA stage +and 280 events at the FA stage. -```{r, echo=FALSE} -# observed vs. planned events +```{r} +# Observed vs. planned events event_observed <- c(180, 280) event_planned <- x$analysis$event -tibble(Analysis = c("IA", "FA"), - `Planned number of events` = event_planned, - `Observed number of events` = event_observed) |> - gt() |> - tab_header("Planned vs. Observed events") +tibble::tibble( + Analysis = c("IA", "FA"), + `Planned number of events` = event_planned, + `Observed number of events` = event_observed +) |> + gt::gt() |> + gt::tab_header("Planned vs. Observed events") ``` -Again, we use `gs_power_npe()` to calculate the updated efficacy and futility bounds. We initially set up `theta = 0`, which provides the crossing probability under the null hypothesis. Users have the flexibility to modify the value of `theta` if they are interested in the crossing probability under alternative hypotheses or other scenarios. We set up `theta0 = 0` for null hypothesis. In this implementation, we set up `theta1` as the weighted average of the piecewise HR, i.e., +Again, we use `gs_power_npe()` to calculate the updated efficacy and +futility bounds. We initially set up `theta = 0`, which provides the +crossing probability under the null hypothesis. Users have the flexibility to +modify the value of `theta` if they are interested in the crossing probability +under alternative hypotheses or other scenarios. +We set up `theta0 = 0` for null hypothesis. In this implementation, we set up +`theta1` as the weighted average of the piecewise HR, i.e., $$ -\log(\text{AHR}) = -\log(\sum_{i=1}^m \text{HR}_m), $$ -where the weight is decided by the observed number of events. For example, assume there are 50 events observed at the first 3 months (HR = 1), and 150 events observed after 3 months (HR = 0.6), we can derive `theta1` as $-\left[\frac{50}{280} \times \log(1) + \frac{230}{280} \times \log(0.6)\right]$. - -Moving forward, we proceed to set up the values of the `info` arguments for the input `theta = 0`. Similar to the examples in the 1-sided design above, we set `info` as `r event_observed`/4, where `r event_observed` is the observed number of events. Because `info0` is under null hypothesis, we leave it and use its default imputed values same as `info`. Furthermore, we leave `info1` as its default imputed values same as `info` with the local null hypothesis approach. Though people can explicitly write the formula of `info1` in terms of the observed events, the difference is very minor, and it would be easiest with unblinded data. - -Lastly, we establish the parameters for the upper and lower bounds. The setup of `upper`, `upar`, `test_upper`, `lower`, `lpar`, and `test_lower` closely resembles the original design, with the exception of the `upar` and `lpar` setup. Here, we have introduced the `timing = ...` parameter. At the IA, the observed number of events is `r event_observed[1]`. If we employ the interim analysis alpha spending strategy, the spending timing at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. Here, `r round(event_planned[2], 0)` represents the planned number of events at FA. The remaining alpha will be allocated to the FA stage. - -```{r} -x_updated <- gs_power_npe(theta = 0, - theta0 = 0, - theta1 = -(50 * log(1) + 230 * log(0.6)) / 280, - # info - info = event_observed / 4, - # upper bound - upper = x$input$upper, - upar = list(sf = gsDesign::sfLDOF, - total_spend = alpha, - param = NULL, - timing = c(event_observed[1] / max(event_planned), 1)), - test_upper = TRUE, - # lower bound - lower = x$input$lower, - lpar = list(sf = gsDesign::sfLDOF, - total_spend = beta, - param = NULL, - timing = c(event_observed[1] / max(event_planned), 1)), - test_lower = c(TRUE, FALSE), - binding = x$input$binding) +where the weight is decided by the observed number of events. For example, +assume there are 50 events observed at the first 3 months (HR = 1), +and 150 events observed after 3 months (HR = 0.6), we can derive `theta1` as +$-\left[\frac{50}{280} \times \log(1) + \frac{230}{280} \times \log(0.6)\right]$. + +Moving forward, we proceed to set up the values of the `info` arguments for +the input `theta = 0`. Similar to the examples in the 1-sided design above, +we set `info` as `r event_observed`/4, where `r event_observed` is the +observed number of events. Because `info0` is under null hypothesis, +we leave it and use its default imputed values same as `info`. +Furthermore, we leave `info1` as its default imputed values same as `info` +with the local null hypothesis approach. Though people can explicitly write +the formula of `info1` in terms of the observed events, the difference is +very minor, and it would be easiest with unblinded data. + +Lastly, we establish the parameters for the upper and lower bounds. +The setup of `upper`, `upar`, `test_upper`, `lower`, `lpar`, and `test_lower` +closely resembles the original design, with the exception of the `upar` and +`lpar` setup. Here, we have introduced the `timing = ...` parameter. +At the IA, the observed number of events is `r event_observed[1]`. +If we employ the interim analysis alpha spending strategy, the spending timing +at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. +Here, `r round(event_planned[2], 0)` represents the planned number of events +at FA. The remaining alpha will be allocated to the FA stage. + +```{r, echo=TRUE} +x_updated <- gs_power_npe( + theta = 0, + theta0 = 0, + theta1 = -(50 * log(1) + 230 * log(0.6)) / 280, + # Info + info = event_observed / 4, + # Upper bound + upper = x$input$upper, + upar = list( + sf = gsDesign::sfLDOF, + total_spend = alpha, + param = NULL, + timing = c(event_observed[1] / max(event_planned), 1) + ), + test_upper = TRUE, + # Lower bound + lower = x$input$lower, + lpar = list( + sf = gsDesign::sfLDOF, + total_spend = beta, + param = NULL, + timing = c(event_observed[1] / max(event_planned), 1) + ), + test_lower = c(TRUE, FALSE), + binding = x$input$binding +) ``` -```{r, echo=FALSE} +```{r} x_updated |> - gt() |> - tab_header(title = "Updated design", subtitle = "by event fraction based timing") |> - tab_style(style = list(cell_fill(color = "lightcyan"), - cell_text(weight = "bold")), - locations = cells_body(columns = probability, - rows = bound == "upper")) |> - tab_style(style = list(cell_fill(color = "#F9E3D6"), - cell_text(weight = "bold")), - locations = cells_body(columns = probability, - rows = bound == "lower")) + gt::gt() |> + gt::tab_header(title = "Updated design", subtitle = "by event fraction based timing") |> + gt::tab_style( + style = list( + gt::cell_fill(color = "lightcyan"), + gt::cell_text(weight = "bold") + ), + locations = gt::cells_body( + columns = probability, + rows = bound == "upper" + ) + ) |> + gt::tab_style( + style = list( + gt::cell_fill(color = "#F9E3D6"), + gt::cell_text(weight = "bold") + ), + locations = gt::cells_body( + columns = probability, + rows = bound == "lower" + ) + ) ``` # References From 7cab218f3c777dde4e74c619d8c8ad25f0d17cc9 Mon Sep 17 00:00:00 2001 From: keaven Date: Wed, 14 Feb 2024 20:47:21 -0500 Subject: [PATCH 10/30] ongoing edits --- vignettes/articles/story-update-boundary.Rmd | 44 +++++++++++++------- 1 file changed, 28 insertions(+), 16 deletions(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 65f021a7..3cd8eeb7 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -28,7 +28,12 @@ library(dplyr) # Design assumptions -We assume a total of two analyses, including a final Analysis (FA). The initial Analysis (IA) occurs at month 20, followed by the FA at month 36. The enrollment period spans 14 months, with the first 2 months having an enrollment rate of 3, the next 2 months with a rate of 6, and the final 10 months with a rate of 9. The planned enrollment rate to achieve the 90\% power is a multiptier of 3, 6, and 9. Additionally, there is a delayed treatment effect in the first 3 months (HR = 1), which transitions to an HR of 0.6 thereafter. The median of the control arm is 9 months, and the dropout rate is 0.0001. +We begin with design assumptions. +We assume two analyses: an interim analysis (IA) final Analysis (FA). The IA is planned for month 20 after opening enrollment, followed by the FA at month 36. +The planned enrollment period spans 14 months, with the first 2 months having an enrollment rate of 1/3 the final rate, the next 2 months with a rate of 2/3 of the final rate, and the final rate for the remaining 10 months. +To obtain the targeted 90\% power, these rates will be multiplied by a constant. +The control arm is assumed to follow an exponential distribution with a median of 9 months and the dropout rate is 0.0001 per month regardless of treatment group. +Finally, the experimental treatment group is piecewise exponential with a 3-month delayed treatment effect; i.e. in the first 3 months HR = 1 and the HR is 0.6 thereafter. ```{r, echo=FALSE} alpha <- 0.025 @@ -37,10 +42,10 @@ ratio <- 1 # enrollment enroll_rate <- define_enroll_rate(duration = c(2, 2, 10), - rate = c(3, 6, 9)) + rate = (1:3)/3) # failure and dropout -fail_rate <- define_fail_rate(duration = c(3, 100), +fail_rate <- define_fail_rate(duration = c(3, Inf), fail_rate = log(2) / 9, hr = c(1, 0.6), dropout_rate = .0001) @@ -48,9 +53,9 @@ fail_rate <- define_fail_rate(duration = c(3, 100), # One-sided design {.tabset} -In this example, we have efficacy bounds at both the IA and FA. The spending function used is the @lan1983discrete spending function with a total alpha of `r alpha`, which approximates an O'Brien-Fleming bound. +For the design, we have efficacy bounds at both the IA and FA. The spending function used is the @lan1983discrete spending function with a total alpha of `r alpha`, which approximates an O'Brien-Fleming bound. -## Planed design +## Planned design ```{r} upper <- gs_spending_bound @@ -71,12 +76,14 @@ x <- gs_design_ahr(enroll_rate = enroll_rate, test_lower = FALSE) |> to_integer() ``` -In the planed design, we have +The planned design targets: - planned events: `r round(x$analysis$event, 0)` -- planned information fraction (timing): `r round(x$analysis$info_frac, 4)` +- planned information fraction at interim and final analysis: `r round(x$analysis$info_frac, 4)` - planned alpha spending: `r sfLDOF(0.025, x$analysis$info_frac)$spend` - planned efficacy bounds: `r round(x$bound$z[x$bound$bound == "upper"], 4)` + +We note that rounding up the final targeted events increases power slightly over the targeted 90%. ```{r, echo=FALSE} x |> @@ -95,7 +102,9 @@ x |> ## Update bounds at time of analysis -Assume in practice, there are 180 and 280 events observed in IA and FA, respectively. +We assume 180 and 280 events observed at the IA and FA, respectively. +We will assume the differences from planned are due to logistical considerations. +We also assume the protocol specifies that the full $\alpha$ will be spent at the final analysis even in a case like this where there is a shortfall of events versus the design plan. ```{r, echo=FALSE} # observed vs. planned events @@ -103,18 +112,21 @@ event_observed <- c(180, 280) event_planned <- x$analysis$event tibble(Analysis = c("IA", "FA"), - `Planned number of events` = event_planned, - `Observed number of events` = event_observed) |> + `Planned events` = event_planned, + `Observed events` = event_observed) |> gt() |> tab_header("Planned vs. Observed events") ``` -Considering the observed number of events mentioned above, we utilize the `gs_power_npe()` function to generate updated efficacy bounds. As `gs_power_npe()` is a lower-level function that may be unfamiliar to some users, we will start by providing an introduction to its arguments as preliminaries. +We will utilize the `gs_power_npe()` function to update efficacy bounds based on the observed events. +We start by providing an introduction to its arguments. - *3 theta-related arguments*: - + `theta`: natural parameter for group sequential design representing expected incremental drift at all analyses, i.e., $-\log(\text{AHR})$. It is used for crossing probability calculation. For example, if we set `theta = 0`, then the crossing probability is the type I error. Default is `0.1`. - + `theta0`: natural parameter for H0 and it decides the upper bounds. Default is `NULL`. - + `theta1`: natural parameter for H1 and it decides the lower bounds. Default is `NULL`. + + `theta`: a vector with the natural parameter for the group sequential design; represents expected ; $-\log(\text{AHR})$ in this case. + This is used for boundary crossing probability calculation. + For example, if we set `theta = 0`, then the crossing probability is the Type I error. + + `theta0`: natural parameter for H0; used for determining the upper (efficacy) spending and bounds. If the default of `NULL` is given, then this is replaced by 0. + + `theta1`: natural parameter for H1; used for determining the lower spending and bounds. The default is `NULL`, in which case the value is replaced with the input `theta`. - *3 statistical information-related arguments*: + `info`: statistical information at all analyses for input `theta`. Default is `1`. @@ -173,7 +185,7 @@ x_updated |> In this section, we investigate a 2 sided asymmetric design, with the non-binding beta-spending futility bounds. Beta-spending refers to error spending for the lower bound crossing probabilities under the alternative hypothesis. Non-binding assumes the trial continues if the lower bound is crossed for Type I and Type II error computation. -## Planed design +## Planned design In the original designs, we employ the Lan DeMets O'Brien-Fleming spending function [@lan1983discrete] for both efficacy and futility bounds. The total spending for efficacy is `r alpha`, and for futility is `r beta`. Besides, we assume the futility test only happens at IA. @@ -199,7 +211,7 @@ x <- gs_design_ahr(enroll_rate = enroll_rate, binding = FALSE) |> to_integer() ``` -In the planed design, we have +In the planned design, we have - planned events: `r round(x$analysis$event, 0)` - planned information fraction (timing): `r round(x$analysis$info_frac, 4)` From eb45801b178d3682a124318ad7da62c05665f553 Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Thu, 15 Feb 2024 10:20:43 -0500 Subject: [PATCH 11/30] fix lintr issues --- vignettes/articles/story-update-boundary.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 5f16756a..46c6a967 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -43,7 +43,7 @@ ratio <- 1 # enrollment enroll_rate <- define_enroll_rate(duration = c(2, 2, 10), - rate = (1:3)/3) + rate = (1:3) / 3) # failure and dropout fail_rate <- define_fail_rate(duration = c(3, Inf), From 96db2d2afb9bf9f165fc8e2c101168a6ac594251 Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Thu, 15 Feb 2024 18:15:45 -0500 Subject: [PATCH 12/30] move `ahr_blinded` from inst to R --- NAMESPACE | 2 + {inst/helper_function => R}/ahr_blinded.R | 56 +++++--- man/ahr_blinded.Rd | 74 ++++++++++ vignettes/articles/story-update-boundary.Rmd | 142 +++++++++++++++++-- 4 files changed, 241 insertions(+), 33 deletions(-) rename {inst/helper_function => R}/ahr_blinded.R (67%) create mode 100644 man/ahr_blinded.Rd diff --git a/NAMESPACE b/NAMESPACE index 8a5a39ab..eee1a1c9 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -9,6 +9,7 @@ S3method(summary,gs_design) S3method(to_integer,fixed_design) S3method(to_integer,gs_design) export(ahr) +export(ahr_blinded) export(as_gt) export(as_rtf) export(define_enroll_rate) @@ -81,6 +82,7 @@ importFrom(stats,pnorm) importFrom(stats,qnorm) importFrom(stats,stepfun) importFrom(stats,uniroot) +importFrom(survival,Surv) importFrom(tibble,tibble) importFrom(utils,tail) useDynLib(gsDesign2, .registration = TRUE) diff --git a/inst/helper_function/ahr_blinded.R b/R/ahr_blinded.R similarity index 67% rename from inst/helper_function/ahr_blinded.R rename to R/ahr_blinded.R index 0473f0ef..596c01d7 100644 --- a/inst/helper_function/ahr_blinded.R +++ b/R/ahr_blinded.R @@ -48,12 +48,9 @@ #' \if{html}{The contents of this section are shown in PDF user manual only.} #' #' @return A \code{tibble} with one row containing -#' `AHR` blinded average hazard ratio based on assumed period-specific hazard ratios input in `failRates` +#' `ahr` blinded average hazard ratio based on assumed period-specific hazard ratios input in `fail_rate` #' and observed events in the corresponding intervals -#' `Events` total observed number of events, `info` statistical information based on Schoenfeld approximation, -#' and info0 (information under related null hypothesis) for each value of `totalDuration` input; -#' if `simple=FALSE`, `Stratum` and `t` (beginning of each constant HR period) are also returned -#' and `HR` is returned instead of `AHR` +#' `event` total observed number of events, and `info0` (information under related null hypothesis). #' #' @examples #' \donttest{ @@ -61,8 +58,8 @@ #' library(survival) #' ahr_blinded( #' Srv = Surv( -#' time = simtrial::Ex2delayedEffect$month, -#' event = simtrial::Ex2delayedEffect$evntd +#' time = simtrial::ex2_delayed_effect$month, +#' event = simtrial::ex2_delayed_effect$evntd #' ), #' intervals = c(4, 100), #' hr = c(1, .55), @@ -71,32 +68,45 @@ #' } #' #' @export -ahr_blinded <- function(Srv = Surv( - time = simtrial::Ex1delayedEffect$month, - event = simtrial::Ex1delayedEffect$evntd - ), - intervals = array(3, 3), +ahr_blinded <- function(Srv = survival::Surv(time = simtrial::ex1_delayed_effect$month, + event = simtrial::ex1_delayed_effect$evntd), + intervals = c(3, Inf), hr = c(1, .6), ratio = 1) { - msg <- "hr must be a vector of positive numbers" - if (!is.vector(hr, mode = "numeric")) stop(msg) - if (min(hr) <= 0) stop(msg) + # Input checking + if (!is.vector(hr, mode = "numeric") | min(hr) <= 0) { + stop("ahr_blinded: hr must be a vector of positive numbers.") + } - events <- simtrial::pwexpfit(Srv, intervals)[, 3] + tte_data <- data.frame(time = Srv[, "time"], status = Srv[, "status"]) + if (nrow(subset(tte_data, time > sum(intervals) & status > 0)) > 0) { + intervals_imputed <- c(intervals, Inf) + } else { + intervals_imputed <- intervals + } + if (length(intervals_imputed) != length(hr)) { + stop("ahr_blinded: the piecewise model specified hr and intervals are not aligned.") + } + + # Fit the survival data into piecewise exponential model + event <- simtrial::fit_pwexp(Srv, intervals)[, 3] nhr <- length(hr) - nx <- length(events) + nx <- length(event) + # Add to hr if length shorter than intervals - if (length(hr) < length(events)) hr <- c(hr, rep(hr[nhr], nx - nhr)) + if (length(hr) < length(event)) { + hr <- c(hr, rep(hr[nhr], nx - nhr)) + } # Compute blinded AHR - theta <- sum(log(hr[1:nx]) * events) / sum(events) + theta <- - sum(log(hr[1:nx]) * event) / sum(event) # Compute adjustment for information - Qe <- ratio / (1 + ratio) + q_e <- ratio / (1 + ratio) ans <- tibble( - Events = sum(events), AHR = exp(theta), - theta = theta, info0 = sum(events) * (1 - Qe) * Qe + event = sum(event), ahr = exp(-theta), + theta = theta, info0 = sum(event) * (1 - q_e) * q_e ) return(ans) -} +} \ No newline at end of file diff --git a/man/ahr_blinded.Rd b/man/ahr_blinded.Rd new file mode 100644 index 00000000..a8628ccd --- /dev/null +++ b/man/ahr_blinded.Rd @@ -0,0 +1,74 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ahr_blinded.R +\name{ahr_blinded} +\alias{ahr_blinded} +\title{Blinded estimation of average hazard ratio} +\usage{ +ahr_blinded( + Srv = survival::Surv(time = simtrial::ex1_delayed_effect$month, event = + simtrial::ex1_delayed_effect$evntd), + intervals = c(3, Inf), + hr = c(1, 0.6), + ratio = 1 +) +} +\arguments{ +\item{Srv}{input survival object (see \code{Surv}); note that only 0=censored, 1=event for \code{Surv}} + +\item{intervals}{Vector containing positive values indicating interval lengths where the +exponential rates are assumed. +Note that a final infinite interval is added if any events occur after the final interval +specified.} + +\item{hr}{vector of hazard ratios assumed for each interval} + +\item{ratio}{ratio of experimental to control randomization.} +} +\value{ +A \code{tibble} with one row containing +\code{ahr} blinded average hazard ratio based on assumed period-specific hazard ratios input in \code{fail_rate} +and observed events in the corresponding intervals +\code{event} total observed number of events, and \code{info0} (information under related null hypothesis). +} +\description{ +Based on blinded data and assumed hazard ratios in different intervals, compute +a blinded estimate of average hazard ratio (AHR) and corresponding estimate of statistical information. +This function is intended for use in computing futility bounds based on spending assuming +the input hazard ratio (hr) values for intervals specified here. +} +\section{Specification}{ + +\if{latex}{ + \itemize{ + \item Validate if input hr is a numeric vector. + \item Validate if input hr is non-negative. + \item Simulate piece-wise exponential survival estimation with the inputs survival object Srv + and intervals. + \item Save the length of hr and events to an object, and if the length of hr is shorter than + the intervals, add replicates of the last element of hr and the corresponding numbers of events + to hr. + \item Compute the blinded estimation of average hazard ratio. + \item Compute adjustment for information. + \item Return a tibble of the sum of events, average hazard raito, blinded average hazard + ratio, and the information. + } +} +\if{html}{The contents of this section are shown in PDF user manual only.} +} + +\examples{ +\donttest{ +library(simtrial) +library(survival) +ahr_blinded( + Srv = Surv( + time = simtrial::ex2_delayed_effect$month, + event = simtrial::ex2_delayed_effect$evntd + ), + intervals = c(4, 100), + hr = c(1, .55), + ratio = 1 +) +} + +} diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 46c6a967..4fafb87d 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -50,7 +50,8 @@ fail_rate <- define_fail_rate(duration = c(3, Inf), fail_rate = log(2) / 9, hr = c(1, 0.6), dropout_rate = .0001) - +# IA and FA analysis time +analysis_time <- c(20, 36) ``` # One-sided design {.tabset} @@ -69,7 +70,7 @@ x <- gs_design_ahr( alpha = alpha, beta = beta, info_frac = NULL, - analysis_time = c(20, 36), + analysis_time = analysis_time, ratio = ratio, upper = gs_spending_bound, upar = upar, @@ -331,14 +332,27 @@ x |> ) ``` -## Update bounds at time of analysis +## Update bounds at time of analysis {.tabset} In practice, let us assume that there were 180 events observed at IA stage and 280 events at the FA stage. ```{r} +set.seed(2024) + +observed_data <- simtrial::sim_pw_surv(n = x$analysis$n[x$analysis$analysis == 2], + stratum = data.frame(stratum = "All", p = 1), + block = c(rep("control", 2), rep("experimental", 2)), + enroll_rate = x$enroll_rate, + fail_rate = (fail_rate |> simtrial::to_sim_pw_surv())$fail_rate, + dropout_rate = (fail_rate |> simtrial::to_sim_pw_surv())$dropout_rate) + +observed_ia_data <- observed_data |> simtrial::cut_data_by_date(analysis_time[1]) +observed_fa_data <- observed_data |> simtrial::cut_data_by_date(analysis_time[2]) + # Observed vs. planned events -event_observed <- c(180, 280) +# event_observed <- c(180, 280) +event_observed <- c(sum(observed_ia_data$event), sum(observed_fa_data$event)) event_planned <- x$analysis$event tibble::tibble( @@ -355,15 +369,67 @@ futility bounds. We initially set up `theta = 0`, which provides the crossing probability under the null hypothesis. Users have the flexibility to modify the value of `theta` if they are interested in the crossing probability under alternative hypotheses or other scenarios. -We set up `theta0 = 0` for null hypothesis. In this implementation, we set up -`theta1` as the weighted average of the piecewise HR, i.e., +We set up `theta0 = 0` for null hypothesis. As for `theta1`, there are 2 options to set it up. + +- **Option 1:** `theta1` is the weighted average of the piecewise HR, i.e., $$ -\log(\text{AHR}) = -\log(\sum_{i=1}^m \text{HR}_m), $$ -where the weight is decided by the observed number of events. For example, -assume there are 50 events observed at the first 3 months (HR = 1), -and 150 events observed after 3 months (HR = 0.6), we can derive `theta1` as -$-\left[\frac{50}{280} \times \log(1) + \frac{230}{280} \times \log(0.6)\right]$. +where the weight is decided by the observed number of events. +```{r} +observed_first_3_month <- observed_data |> simtrial::cut_data_by_date(3) +``` +In this example, there are `r sum(observed_first_3_month$event)` events at the first 3 months (HR = 1), and `r event_observed[2] - sum(observed_first_3_month$event)` events observed after 3 months (HR = 0.6), we can derive `theta1` as +`r sum(observed_first_3_month$event)` / `r event_observed[2]` $\times \log(1) +$ `r event_observed[2] - sum(observed_first_3_month$event)` $\times \log(0.6)$ = `r -(sum(observed_first_3_month$event) * log(1) + (event_observed[2] - sum(observed_first_3_month$event)) * log(0.6)) / event_observed[2]`. + +- **Option 2:** we use the blinded data to estimate the AHR and `theta1` is the negative logarithm of the estimated AHR. +```{r} +# source(file.path(system.file("helper_function", package = "gsDesign2"), "ahr_blinded.R")) + +ahr_blinded <- function(Srv = survival::Surv( + time = simtrial::Ex1delayedEffect$month, + event = simtrial::Ex1delayedEffect$evntd + ), + intervals = array(3, 3), + hr = c(1, .6), + ratio = 1) { + msg <- "hr must be a vector of positive numbers" + if (!is.vector(hr, mode = "numeric")) stop(msg) + if (min(hr) <= 0) stop(msg) + + event <- simtrial::fit_pwexp(Srv, intervals)[, 3] + nhr <- length(hr) + nx <- length(event) + # Add to hr if length shorter than intervals + if (length(hr) < length(event)) hr <- c(hr, rep(hr[nhr], nx - nhr)) + + # Compute blinded AHR + theta <- - sum(log(hr[1:nx]) * event) / sum(event) + + # Compute adjustment for information + q_e <- ratio / (1 + ratio) + + ans <- tibble( + event = sum(event), ahr = exp(-theta), + theta = theta, info0 = sum(event) * (1 - q_e) * q_e + ) + return(ans) +} +``` + +```{r} +blinded_ahr <- ahr_blinded(Srv = survival::Surv(time = observed_fa_data$tte, + event = observed_fa_data$event), + intervals = c(3, Inf), + hr = c(1, 0.6), + ratio = 1) + +blinded_ahr |> + gt::gt() |> + gt::tab_header("Blinded estimation of average hazard ratio") +``` + +In this example, the estimated AHR is `r blinded_ahr$ahr`, and accordingly, the `theta1` values is `r blinded_ahr$theta`. Moving forward, we proceed to set up the values of the `info` arguments for the input `theta = 0`. Similar to the examples in the 1-sided design above, @@ -385,6 +451,7 @@ at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. Here, `r round(event_planned[2], 0)` represents the planned number of events at FA. The remaining alpha will be allocated to the FA stage. +### Option 1: `theta1` is the weighted average ```{r, echo=TRUE} x_updated <- gs_power_npe( theta = 0, @@ -440,4 +507,59 @@ x_updated |> ) ``` +### Option 2: `theta1` is based on blinded ahr estimate +```{r, echo=TRUE} +x_updated <- gs_power_npe( + theta = 0, + theta0 = 0, + theta1 = blinded_ahr$theta, + # Info + info = event_observed / 4, + # Upper bound + upper = x$input$upper, + upar = list( + sf = gsDesign::sfLDOF, + total_spend = alpha, + param = NULL, + timing = c(event_observed[1] / max(event_planned), 1) + ), + test_upper = TRUE, + # Lower bound + lower = x$input$lower, + lpar = list( + sf = gsDesign::sfLDOF, + total_spend = beta, + param = NULL, + timing = c(event_observed[1] / max(event_planned), 1) + ), + test_lower = c(TRUE, FALSE), + binding = x$input$binding +) +``` + +```{r} +x_updated |> + gt::gt() |> + gt::tab_header(title = "Updated design", subtitle = "by event fraction based timing") |> + gt::tab_style( + style = list( + gt::cell_fill(color = "lightcyan"), + gt::cell_text(weight = "bold") + ), + locations = gt::cells_body( + columns = probability, + rows = bound == "upper" + ) + ) |> + gt::tab_style( + style = list( + gt::cell_fill(color = "#F9E3D6"), + gt::cell_text(weight = "bold") + ), + locations = gt::cells_body( + columns = probability, + rows = bound == "lower" + ) + ) +``` # References From 96aef191fbdfaffe886489ae187f5b9e6d3c136d Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Thu, 15 Feb 2024 18:15:56 -0500 Subject: [PATCH 13/30] update the vignette --- vignettes/articles/story-update-boundary.Rmd | 152 +++++++------------ 1 file changed, 58 insertions(+), 94 deletions(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 4fafb87d..1ee8e65e 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -135,43 +135,9 @@ tibble(Analysis = c("IA", "FA"), tab_header("Planned vs. Observed events") ``` -We will utilize the `gs_power_npe()` function to update efficacy bounds based on the observed events. -We start by providing an introduction to its arguments. - -- *3 theta-related arguments*: - + `theta`: a vector with the natural parameter for the group sequential design; represents expected ; $-\log(\text{AHR})$ in this case. - This is used for boundary crossing probability calculation. - For example, if we set `theta = 0`, then the crossing probability is the Type I error. - + `theta0`: natural parameter for H0; used for determining the upper (efficacy) spending and bounds. If the default of `NULL` is given, then this is replaced by 0. - + `theta1`: natural parameter for H1; used for determining the lower spending and bounds. The default is `NULL`, in which case the value is replaced with the input `theta`. - -- *3 statistical information-related arguments*: - + `info`: statistical information at all analyses for input `theta`. Default is `1`. - + `info0`: statistical information under H0, if different than `info`. It impacts null hypothesis bound calculation. Default is `NULL`. - + `info1`: statistical information under hypothesis used for futility bound calculation if different from `info`. It impacts futility hypothesis bound calculation. Default is `NULL`. - -- *Efficacy/Futility boundary-related arguments*: - + `upper`, `upar` and `test_upper`: efficacy bounds related, same as `gs_design_ahr()`. - + `lower`, `lpar` and `test_lower`: futility bounds related, same as `gs_design_ahr()`. - + `binding`: `TRUE` or `FALSE` indicating whether it is binding or non-binding designs. - -When implementing `gs_power_npe()`, we initially set up `theta = 0`, which provides the crossing probability under the null hypothesis. Users have the flexibility to modify the value of `theta` if they are interested in the crossing probability under alternative hypotheses or other scenarios. In this implementation, we leave the setup of `theta0` with its default imputed value of `0`, indicating the null hypothesis. Additionally, as we are solely dealing with efficacy bounds, we do not adjust `theta1`, which influences the futility bound. - -- *3 statistical information-related arguments*: - - `info`: statistical information at all analyses for input `theta`. - Default is `1`. - - `info0`: statistical information under H0, if different than `info`. - It impacts null hypothesis bound calculation. Default is `NULL`. - - `info1`: statistical information under hypothesis used for futility bound - calculation if different from `info`. - It impacts futility hypothesis bound calculation. Default is `NULL`. - -- *Efficacy/Futility boundary-related arguments*: - - `upper`, `upar` and `test_upper`: efficacy bounds related, same as `gs_design_ahr()`. - - `lower`, `lpar` and `test_lower`: futility bounds related, same as `gs_design_ahr()`. - - `binding`: `TRUE` or `FALSE` indicating whether it is binding or non-binding designs. - -When implementing `gs_power_npe()`, we initially set up `theta = 0`, +We will utilize the `gs_power_npe()` function to update efficacy bounds based on the observed events. +The deatils of its arguments are explained in the Appendix. +We initially set up `theta = 0`, which provides the crossing probability under the null hypothesis. Users have the flexibility to modify the value of `theta` if they are interested in the crossing probability under alternative hypotheses or @@ -334,10 +300,9 @@ x |> ## Update bounds at time of analysis {.tabset} -In practice, let us assume that there were 180 events observed at IA stage -and 280 events at the FA stage. +In practice, let us assume the observed data is generaed by `simtrial::sim_pw_surv`. -```{r} +```{r, echo=TRUE} set.seed(2024) observed_data <- simtrial::sim_pw_surv(n = x$analysis$n[x$analysis$analysis == 2], @@ -347,14 +312,16 @@ observed_data <- simtrial::sim_pw_surv(n = x$analysis$n[x$analysis$analysis == 2 fail_rate = (fail_rate |> simtrial::to_sim_pw_surv())$fail_rate, dropout_rate = (fail_rate |> simtrial::to_sim_pw_surv())$dropout_rate) +observed_first_3_month <- observed_data |> simtrial::cut_data_by_date(3) observed_ia_data <- observed_data |> simtrial::cut_data_by_date(analysis_time[1]) observed_fa_data <- observed_data |> simtrial::cut_data_by_date(analysis_time[2]) # Observed vs. planned events -# event_observed <- c(180, 280) event_observed <- c(sum(observed_ia_data$event), sum(observed_fa_data$event)) event_planned <- x$analysis$event +``` +```{r} tibble::tibble( Analysis = c("IA", "FA"), `Planned number of events` = event_planned, @@ -369,55 +336,30 @@ futility bounds. We initially set up `theta = 0`, which provides the crossing probability under the null hypothesis. Users have the flexibility to modify the value of `theta` if they are interested in the crossing probability under alternative hypotheses or other scenarios. -We set up `theta0 = 0` for null hypothesis. As for `theta1`, there are 2 options to set it up. +We set up `theta0 = 0` for null hypothesis. As for `theta1`, there are 2 options to set it up. +These 2 options arrive at the same value of `theta1`. - **Option 1:** `theta1` is the weighted average of the piecewise HR, i.e., $$ - -\log(\text{AHR}) = -\log(\sum_{i=1}^m \text{HR}_m), + -\log(\text{AHR}) = -\log(\sum_{i=1}^m w_i \; \text{HR}_m), $$ -where the weight is decided by the observed number of events. -```{r} -observed_first_3_month <- observed_data |> simtrial::cut_data_by_date(3) +where the weight is decided by the ratio of observed number of events and the observed final events. +In this example, the number of observed events at the first 3 months (HR = 1) can be derived as +```{r, echo=TRUE} +event_first_3_month <- sum(observed_data$fail_time < 3) +event_first_3_month +``` +and the number of observed events after 3 months (HR = 0.6) are +```{r, echo=TRUE} +event_after_3_month <- event_observed[2] - event_first_3_month +event_after_3_month ``` -In this example, there are `r sum(observed_first_3_month$event)` events at the first 3 months (HR = 1), and `r event_observed[2] - sum(observed_first_3_month$event)` events observed after 3 months (HR = 0.6), we can derive `theta1` as -`r sum(observed_first_3_month$event)` / `r event_observed[2]` $\times \log(1) +$ `r event_observed[2] - sum(observed_first_3_month$event)` $\times \log(0.6)$ = `r -(sum(observed_first_3_month$event) * log(1) + (event_observed[2] - sum(observed_first_3_month$event)) * log(0.6)) / event_observed[2]`. +We can derive `theta1` as +-[`r event_first_3_month` / `r event_observed[2]` $\times \log(1) +$ `r event_after_3_month` / `r event_observed[2]` $\times \log(0.6)$] = `r round(-(event_first_3_month * log(1) + event_after_3_month * log(0.6)) / event_observed[2], 4)`. - **Option 2:** we use the blinded data to estimate the AHR and `theta1` is the negative logarithm of the estimated AHR. -```{r} -# source(file.path(system.file("helper_function", package = "gsDesign2"), "ahr_blinded.R")) - -ahr_blinded <- function(Srv = survival::Surv( - time = simtrial::Ex1delayedEffect$month, - event = simtrial::Ex1delayedEffect$evntd - ), - intervals = array(3, 3), - hr = c(1, .6), - ratio = 1) { - msg <- "hr must be a vector of positive numbers" - if (!is.vector(hr, mode = "numeric")) stop(msg) - if (min(hr) <= 0) stop(msg) - - event <- simtrial::fit_pwexp(Srv, intervals)[, 3] - nhr <- length(hr) - nx <- length(event) - # Add to hr if length shorter than intervals - if (length(hr) < length(event)) hr <- c(hr, rep(hr[nhr], nx - nhr)) - - # Compute blinded AHR - theta <- - sum(log(hr[1:nx]) * event) / sum(event) - - # Compute adjustment for information - q_e <- ratio / (1 + ratio) - - ans <- tibble( - event = sum(event), ahr = exp(-theta), - theta = theta, info0 = sum(event) * (1 - q_e) * q_e - ) - return(ans) -} -``` -```{r} +```{r, echo=TRUE} blinded_ahr <- ahr_blinded(Srv = survival::Surv(time = observed_fa_data$tte, event = observed_fa_data$event), intervals = c(3, Inf), @@ -429,7 +371,11 @@ blinded_ahr |> gt::tab_header("Blinded estimation of average hazard ratio") ``` -In this example, the estimated AHR is `r blinded_ahr$ahr`, and accordingly, the `theta1` values is `r blinded_ahr$theta`. +In this example, the estimated AHR is `r round(blinded_ahr$ahr, 4)`, and accordingly, the `theta1` values is `r round(blinded_ahr$theta, 4)`. + +Both Option 1 and Option 2 yield the same value of theta1. Option 1 is applicable +when the event counts are available, while Option 2 is suitable for scenarios where +time-to-event data is available, without the corresponding event counts. Moving forward, we proceed to set up the values of the `info` arguments for the input `theta = 0`. Similar to the examples in the 1-sided design above, @@ -441,22 +387,15 @@ with the local null hypothesis approach. Though people can explicitly write the formula of `info1` in terms of the observed events, the difference is very minor, and it would be easiest with unblinded data. -Lastly, we establish the parameters for the upper and lower bounds. -The setup of `upper`, `upar`, `test_upper`, `lower`, `lpar`, and `test_lower` -closely resembles the original design, with the exception of the `upar` and -`lpar` setup. Here, we have introduced the `timing = ...` parameter. -At the IA, the observed number of events is `r event_observed[1]`. -If we employ the interim analysis alpha spending strategy, the spending timing -at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. -Here, `r round(event_planned[2], 0)` represents the planned number of events -at FA. The remaining alpha will be allocated to the FA stage. +Lastly, we establish the parameters for the upper and lower bounds, +following the similar procedure as we have in the 1-sided design. ### Option 1: `theta1` is the weighted average ```{r, echo=TRUE} x_updated <- gs_power_npe( theta = 0, theta0 = 0, - theta1 = -(50 * log(1) + 230 * log(0.6)) / 280, + theta1 = -(event_first_3_month * log(1) + event_after_3_month * log(0.6)) / event_observed[2], # Info info = event_observed / 4, # Upper bound @@ -507,7 +446,7 @@ x_updated |> ) ``` -### Option 2: `theta1` is based on blinded ahr estimate +### Option 2: `theta1` is based on blinded ahr estimation ```{r, echo=TRUE} x_updated <- gs_power_npe( theta = 0, @@ -562,4 +501,29 @@ x_updated |> ) ) ``` + +# Appendix + +## Arguments of `gs_power_npe()` + +We provide an introduction to its arguments. + +- *3 theta-related arguments*: + + `theta`: a vector with the natural parameter for the group sequential design; represents expected ; $-\log(\text{AHR})$ in this case. + This is used for boundary crossing probability calculation. + For example, if we set `theta = 0`, then the crossing probability is the Type I error. + + `theta0`: natural parameter for H0; used for determining the upper (efficacy) spending and bounds. If the default of `NULL` is given, then this is replaced by 0. + + `theta1`: natural parameter for H1; used for determining the lower spending and bounds. The default is `NULL`, in which case the value is replaced with the input `theta`. + +- *3 statistical information-related arguments*: + + `info`: statistical information at all analyses for input `theta`. Default is `1`. + + `info0`: statistical information under H0, if different than `info`. It impacts null hypothesis bound calculation. Default is `NULL`. + + `info1`: statistical information under hypothesis used for futility bound calculation if different from `info`. It impacts futility hypothesis bound calculation. Default is `NULL`. + +- *Efficacy/Futility boundary-related arguments*: + + `upper`, `upar` and `test_upper`: efficacy bounds related, same as `gs_design_ahr()`. + + `lower`, `lpar` and `test_lower`: futility bounds related, same as `gs_design_ahr()`. + + `binding`: `TRUE` or `FALSE` indicating whether it is binding or non-binding designs. + + # References From 22d9b35ca0eac4524e89297958587b709da27d92 Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Fri, 16 Feb 2024 10:40:28 -0500 Subject: [PATCH 14/30] fix cmd check --- DESCRIPTION | 3 ++- R/ahr_blinded.R | 16 ++++++++-------- _pkgdown.yml | 1 + man/ahr_blinded.Rd | 8 ++++---- vignettes/articles/story-update-boundary.Rmd | 10 +++++----- 5 files changed, 20 insertions(+), 18 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 38cc6839..aaf59bc0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -48,7 +48,8 @@ Imports: tibble, tidyr, utils, - Rcpp + Rcpp, + survival Suggests: covr, ggplot2, diff --git a/R/ahr_blinded.R b/R/ahr_blinded.R index 596c01d7..892b9a32 100644 --- a/R/ahr_blinded.R +++ b/R/ahr_blinded.R @@ -21,7 +21,7 @@ #' the input hazard ratio (hr) values for intervals specified here. #' @importFrom tibble tibble #' @importFrom survival Surv -#' @param Srv input survival object (see \code{Surv}); note that only 0=censored, 1=event for \code{Surv} +#' @param surv input survival object (see \code{Surv}); note that only 0=censored, 1=event for \code{Surv} #' @param intervals Vector containing positive values indicating interval lengths where the #' exponential rates are assumed. #' Note that a final infinite interval is added if any events occur after the final interval @@ -34,7 +34,7 @@ #' \itemize{ #' \item Validate if input hr is a numeric vector. #' \item Validate if input hr is non-negative. -#' \item Simulate piece-wise exponential survival estimation with the inputs survival object Srv +#' \item Simulate piece-wise exponential survival estimation with the inputs survival object surv #' and intervals. #' \item Save the length of hr and events to an object, and if the length of hr is shorter than #' the intervals, add replicates of the last element of hr and the corresponding numbers of events @@ -57,7 +57,7 @@ #' library(simtrial) #' library(survival) #' ahr_blinded( -#' Srv = Surv( +#' surv = Surv( #' time = simtrial::ex2_delayed_effect$month, #' event = simtrial::ex2_delayed_effect$evntd #' ), @@ -68,17 +68,17 @@ #' } #' #' @export -ahr_blinded <- function(Srv = survival::Surv(time = simtrial::ex1_delayed_effect$month, +ahr_blinded <- function(surv = survival::Surv(time = simtrial::ex1_delayed_effect$month, event = simtrial::ex1_delayed_effect$evntd), intervals = c(3, Inf), hr = c(1, .6), ratio = 1) { # Input checking - if (!is.vector(hr, mode = "numeric") | min(hr) <= 0) { + if (!is.vector(hr, mode = "numeric") || min(hr) <= 0) { stop("ahr_blinded: hr must be a vector of positive numbers.") } - tte_data <- data.frame(time = Srv[, "time"], status = Srv[, "status"]) + tte_data <- data.frame(time = surv[, "time"], status = surv[, "status"]) if (nrow(subset(tte_data, time > sum(intervals) & status > 0)) > 0) { intervals_imputed <- c(intervals, Inf) } else { @@ -89,7 +89,7 @@ ahr_blinded <- function(Srv = survival::Surv(time = simtrial::ex1_delayed_effect } # Fit the survival data into piecewise exponential model - event <- simtrial::fit_pwexp(Srv, intervals)[, 3] + event <- simtrial::fit_pwexp(surv, intervals)[, 3] nhr <- length(hr) nx <- length(event) @@ -109,4 +109,4 @@ ahr_blinded <- function(Srv = survival::Surv(time = simtrial::ex1_delayed_effect theta = theta, info0 = sum(event) * (1 - q_e) * q_e ) return(ans) -} \ No newline at end of file +} diff --git a/_pkgdown.yml b/_pkgdown.yml index 3c409c34..3d1f1dbd 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -46,6 +46,7 @@ reference: - gs_info_ahr - gs_power_ahr - gs_design_ahr + - ahr_blinded - title: "Weighted logrank" desc: > Functions for the weighted logrank test (WLR) method. diff --git a/man/ahr_blinded.Rd b/man/ahr_blinded.Rd index a8628ccd..419e6987 100644 --- a/man/ahr_blinded.Rd +++ b/man/ahr_blinded.Rd @@ -5,7 +5,7 @@ \title{Blinded estimation of average hazard ratio} \usage{ ahr_blinded( - Srv = survival::Surv(time = simtrial::ex1_delayed_effect$month, event = + surv = survival::Surv(time = simtrial::ex1_delayed_effect$month, event = simtrial::ex1_delayed_effect$evntd), intervals = c(3, Inf), hr = c(1, 0.6), @@ -13,7 +13,7 @@ ahr_blinded( ) } \arguments{ -\item{Srv}{input survival object (see \code{Surv}); note that only 0=censored, 1=event for \code{Surv}} +\item{surv}{input survival object (see \code{Surv}); note that only 0=censored, 1=event for \code{Surv}} \item{intervals}{Vector containing positive values indicating interval lengths where the exponential rates are assumed. @@ -42,7 +42,7 @@ the input hazard ratio (hr) values for intervals specified here. \itemize{ \item Validate if input hr is a numeric vector. \item Validate if input hr is non-negative. - \item Simulate piece-wise exponential survival estimation with the inputs survival object Srv + \item Simulate piece-wise exponential survival estimation with the inputs survival object surv and intervals. \item Save the length of hr and events to an object, and if the length of hr is shorter than the intervals, add replicates of the last element of hr and the corresponding numbers of events @@ -61,7 +61,7 @@ the input hazard ratio (hr) values for intervals specified here. library(simtrial) library(survival) ahr_blinded( - Srv = Surv( + surv = Surv( time = simtrial::ex2_delayed_effect$month, event = simtrial::ex2_delayed_effect$evntd ), diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 1ee8e65e..81eaaf9c 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -336,7 +336,7 @@ futility bounds. We initially set up `theta = 0`, which provides the crossing probability under the null hypothesis. Users have the flexibility to modify the value of `theta` if they are interested in the crossing probability under alternative hypotheses or other scenarios. -We set up `theta0 = 0` for null hypothesis. As for `theta1`, there are 2 options to set it up. +We set up `theta0 = 0` for null hypothesis. As for `theta1`, there are 2 options to set it up. These 2 options arrive at the same value of `theta1`. - **Option 1:** `theta1` is the weighted average of the piecewise HR, i.e., @@ -360,14 +360,14 @@ We can derive `theta1` as - **Option 2:** we use the blinded data to estimate the AHR and `theta1` is the negative logarithm of the estimated AHR. ```{r, echo=TRUE} -blinded_ahr <- ahr_blinded(Srv = survival::Surv(time = observed_fa_data$tte, - event = observed_fa_data$event), +blinded_ahr <- ahr_blinded(surv = survival::Surv(time = observed_fa_data$tte, + event = observed_fa_data$event), intervals = c(3, Inf), hr = c(1, 0.6), ratio = 1) -blinded_ahr |> - gt::gt() |> +blinded_ahr |> + gt::gt() |> gt::tab_header("Blinded estimation of average hazard ratio") ``` From bdf77e782bfe00942b23aca75a236e497e41f5d4 Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Fri, 16 Feb 2024 11:12:07 -0500 Subject: [PATCH 15/30] fix cmd check --- DESCRIPTION | 3 ++- R/ahr_blinded.R | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index aaf59bc0..58f2ca78 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -56,7 +56,8 @@ Suggests: kableExtra, knitr, rmarkdown, - testthat + testthat, + simtrial VignetteBuilder: knitr LinkingTo: diff --git a/R/ahr_blinded.R b/R/ahr_blinded.R index 892b9a32..7f601da5 100644 --- a/R/ahr_blinded.R +++ b/R/ahr_blinded.R @@ -69,7 +69,7 @@ #' #' @export ahr_blinded <- function(surv = survival::Surv(time = simtrial::ex1_delayed_effect$month, - event = simtrial::ex1_delayed_effect$evntd), + event = simtrial::ex1_delayed_effect$evntd), intervals = c(3, Inf), hr = c(1, .6), ratio = 1) { From 18a26d563750daffb4e4d5960843109497c75cec Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Fri, 16 Feb 2024 11:52:30 -0500 Subject: [PATCH 16/30] move some text to appendix --- vignettes/articles/story-update-boundary.Rmd | 154 +++++++++++-------- 1 file changed, 91 insertions(+), 63 deletions(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 81eaaf9c..89d0391a 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -52,6 +52,9 @@ fail_rate <- define_fail_rate(duration = c(3, Inf), dropout_rate = .0001) # IA and FA analysis time analysis_time <- c(20, 36) + +# randomization ratio +ratio <- 1 ``` # One-sided design {.tabset} @@ -123,11 +126,13 @@ We assume 180 and 280 events observed at the IA and FA, respectively. We will assume the differences from planned are due to logistical considerations. We also assume the protocol specifies that the full $\alpha$ will be spent at the final analysis even in a case like this where there is a shortfall of events versus the design plan. -```{r} +```{r, echo=TRUE} # Observed vs. planned events event_observed <- c(180, 280) event_planned <- x$analysis$event +``` +```{r} tibble(Analysis = c("IA", "FA"), `Planned events` = event_planned, `Observed events` = event_observed) |> @@ -136,53 +141,23 @@ tibble(Analysis = c("IA", "FA"), ``` We will utilize the `gs_power_npe()` function to update efficacy bounds based on the observed events. -The deatils of its arguments are explained in the Appendix. -We initially set up `theta = 0`, -which provides the crossing probability under the null hypothesis. -Users have the flexibility to modify the value of `theta` if they are -interested in the crossing probability under alternative hypotheses or -other scenarios. In this implementation, we leave the setup of `theta0` -with its default imputed value of `0`, indicating the null hypothesis. -Additionally, as we are solely dealing with efficacy bounds, -we do not adjust `theta1`, which influences the futility bound. - -Moving forward, we proceed to set up the values of the `info` arguments -for the input `theta = 0`. Since the statistical information is event-based, -and `theta = 0` (null hypothesis), the observed statistical information -under null hypothesis is -$$ - \text{observed number of events} \times \frac{r}{1 + r} \times \frac{1}{1 + r}, -$$ -where $r$ is the randomization ratio (experimental : control). -In this example, the observed statistical information is calculated as -`info = event_observed / 4`, and thus, we input `r event_observed / 4` -for the `info` argument. In this case, we retain the default imputed value -for `info0`, which is set as `info`. The `info` represents the -statistical information under the null hypothesis (`theta = 0`). -Furthermore, we leave `info1` as it is, considering our focus on efficacy bounds, -while noting that `info1` does affect the futility bounds. - -Lastly, we establish the parameters for the upper and lower bounds. -The setup of `upper`, `upar`, `test_upper`, `lower`, `lpar`, and `test_lower` -closely resembles the original design, with the exception of the `upar` setup. -Here, we have introduced the `timing = ...` parameter. -At the IA, the observed number of events is `r event_observed[1]`. -If we employ the interim analysis alpha spending strategy, the spending timing -at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. -Here, `r round(event_planned[2], 0)` represents the planned number of events -at FA. The remaining alpha will be allocated to the FA stage. +The details of its arguments and implementations are explained in the Appendix. ```{r, echo=TRUE} x_updated <- gs_power_npe( + # `theta = 0` provides the crossing probability under H0. + # Users have the flexibility to modify the value of `theta, + # if they are interested in the crossing probability under H1 or other scenarios. theta = 0, - # Info - info = event_observed / 4, + # Observed statistical information under H0 with equal randomization + info = event_observed * ratio / (1 + ratio)^2, # Upper bound upper = gs_spending_bound, upar = list( sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL, + # The remaining alpha will be allocated to the FA stage. timing = c(event_observed[1] / max(event_planned), 1) ), test_upper = TRUE, @@ -300,7 +275,7 @@ x |> ## Update bounds at time of analysis {.tabset} -In practice, let us assume the observed data is generaed by `simtrial::sim_pw_surv`. +In practice, let us assume the observed data is generated by `simtrial::sim_pw_surv()`. ```{r, echo=TRUE} set.seed(2024) @@ -312,7 +287,6 @@ observed_data <- simtrial::sim_pw_surv(n = x$analysis$n[x$analysis$analysis == 2 fail_rate = (fail_rate |> simtrial::to_sim_pw_surv())$fail_rate, dropout_rate = (fail_rate |> simtrial::to_sim_pw_surv())$dropout_rate) -observed_first_3_month <- observed_data |> simtrial::cut_data_by_date(3) observed_ia_data <- observed_data |> simtrial::cut_data_by_date(analysis_time[1]) observed_fa_data <- observed_data |> simtrial::cut_data_by_date(analysis_time[2]) @@ -331,12 +305,11 @@ tibble::tibble( gt::tab_header("Planned vs. Observed events") ``` -Again, we use `gs_power_npe()` to calculate the updated efficacy and -futility bounds. We initially set up `theta = 0`, which provides the -crossing probability under the null hypothesis. Users have the flexibility to -modify the value of `theta` if they are interested in the crossing probability -under alternative hypotheses or other scenarios. -We set up `theta0 = 0` for null hypothesis. As for `theta1`, there are 2 options to set it up. +Again, we use `gs_power_npe()` to calculate the updated efficacy and futility bounds. +The details of its arguments and implementations are explained in the Appendix. +We initially set up `theta = 0` for crossing probability under the null hypothesis +and `theta0 = 0` for null hypothesis. +As for `theta1`, there are 2 options to set it up. These 2 options arrive at the same value of `theta1`. - **Option 1:** `theta1` is the weighted average of the piecewise HR, i.e., @@ -365,7 +338,9 @@ blinded_ahr <- ahr_blinded(surv = survival::Surv(time = observed_fa_data$tte, intervals = c(3, Inf), hr = c(1, 0.6), ratio = 1) +``` +```{r} blinded_ahr |> gt::gt() |> gt::tab_header("Blinded estimation of average hazard ratio") @@ -377,33 +352,27 @@ Both Option 1 and Option 2 yield the same value of theta1. Option 1 is applicabl when the event counts are available, while Option 2 is suitable for scenarios where time-to-event data is available, without the corresponding event counts. -Moving forward, we proceed to set up the values of the `info` arguments for -the input `theta = 0`. Similar to the examples in the 1-sided design above, -we set `info` as `r event_observed`/4, where `r event_observed` is the -observed number of events. Because `info0` is under null hypothesis, -we leave it and use its default imputed values same as `info`. -Furthermore, we leave `info1` as its default imputed values same as `info` -with the local null hypothesis approach. Though people can explicitly write -the formula of `info1` in terms of the observed events, the difference is -very minor, and it would be easiest with unblinded data. - -Lastly, we establish the parameters for the upper and lower bounds, -following the similar procedure as we have in the 1-sided design. - ### Option 1: `theta1` is the weighted average ```{r, echo=TRUE} x_updated <- gs_power_npe( + # `theta = 0` provides the crossing probability under H0. + # Users have the flexibility to modify the value of `theta, + # if they are interested in the crossing probability under H1 or other scenarios. theta = 0, + # -log(AHR) for H0 theta0 = 0, + # -log(AHR) for H1 + # It is used for determining the lower spending and bounds. theta1 = -(event_first_3_month * log(1) + event_after_3_month * log(0.6)) / event_observed[2], - # Info - info = event_observed / 4, + # Observed statistical information under H0 with equal randomization + info = event_observed * ratio / (1 + ratio)^2, # Upper bound upper = x$input$upper, upar = list( sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL, + # The remaining alpha will be allocated to the FA stage. timing = c(event_observed[1] / max(event_planned), 1) ), test_upper = TRUE, @@ -413,6 +382,7 @@ x_updated <- gs_power_npe( sf = gsDesign::sfLDOF, total_spend = beta, param = NULL, + # The remaining alpha will be allocated to the FA stage. timing = c(event_observed[1] / max(event_planned), 1) ), test_lower = c(TRUE, FALSE), @@ -449,17 +419,24 @@ x_updated |> ### Option 2: `theta1` is based on blinded ahr estimation ```{r, echo=TRUE} x_updated <- gs_power_npe( + # `theta = 0` provides the crossing probability under H0. + # Users have the flexibility to modify the value of `theta, + # if they are interested in the crossing probability under H1 or other scenarios. theta = 0, + # -log(AHR) for H0 theta0 = 0, + # -log(AHR) for H1 + # It is used for determining the lower spending and bounds. theta1 = blinded_ahr$theta, - # Info - info = event_observed / 4, + # Observed statistical information under H0 with equal randomization + info = event_observed * ratio / (1 + ratio)^2, # Upper bound upper = x$input$upper, upar = list( sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL, + # The remaining alpha will be allocated to the FA stage. timing = c(event_observed[1] / max(event_planned), 1) ), test_upper = TRUE, @@ -469,6 +446,7 @@ x_updated <- gs_power_npe( sf = gsDesign::sfLDOF, total_spend = beta, param = NULL, + # The remaining alpha will be allocated to the FA stage. timing = c(event_observed[1] / max(event_planned), 1) ), test_lower = c(TRUE, FALSE), @@ -525,5 +503,55 @@ We provide an introduction to its arguments. + `lower`, `lpar` and `test_lower`: futility bounds related, same as `gs_design_ahr()`. + `binding`: `TRUE` or `FALSE` indicating whether it is binding or non-binding designs. +## Explanations of `gs_power_npe()` arguments set up in the one-sided design +We initially set up `theta = 0`, +which provides the crossing probability under the null hypothesis. +Users have the flexibility to modify the value of `theta` if they are +interested in the crossing probability under alternative hypotheses or +other scenarios. In this implementation, we leave the setup of `theta0` +with its default imputed value of `0`, indicating the null hypothesis. +Additionally, as we are solely dealing with efficacy bounds, +we do not adjust `theta1`, which influences the futility bound. + +Moving forward, we proceed to set up the values of the `info` arguments +for the input `theta = 0`. Since the statistical information is event-based, +and `theta = 0` (null hypothesis), the observed statistical information +under null hypothesis is +$$ + \text{observed number of events} \times \frac{r}{1 + r} \times \frac{1}{1 + r}, +$$ +where $r$ is the randomization ratio (experimental : control). +In this example, the observed statistical information is calculated as +`info = event_observed / 4`, and thus, we input `r event_observed / 4` +for the `info` argument. In this case, we retain the default imputed value +for `info0`, which is set as `info`. The `info` represents the +statistical information under the null hypothesis (`theta = 0`). +Furthermore, we leave `info1` as it is, considering our focus on efficacy bounds, +while noting that `info1` does affect the futility bounds. + +Lastly, we establish the parameters for the upper and lower bounds. +The setup of `upper`, `upar`, `test_upper`, `lower`, `lpar`, and `test_lower` +closely resembles the original design, with the exception of the `upar` setup. +Here, we have introduced the `timing = ...` parameter. +At the IA, the observed number of events is `r event_observed[1]`. +If we employ the interim analysis alpha spending strategy, the spending timing +at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. +Here, `r round(event_planned[2], 0)` represents the planned number of events +at FA. The remaining alpha will be allocated to the FA stage. + +## Explanations of `gs_power_npe()` arguments set up in the two-sided asymmetric design, beta-spending with non-binding lower bound + +Moving forward, we proceed to set up the values of the `info` arguments for +the input `theta = 0`. Similar to the examples in the 1-sided design above, +we set `info` as `r event_observed`/4, where `r event_observed` is the +observed number of events. Because `info0` is under null hypothesis, +we leave it and use its default imputed values same as `info`. +Furthermore, we leave `info1` as its default imputed values same as `info` +with the local null hypothesis approach. Though people can explicitly write +the formula of `info1` in terms of the observed events, the difference is +very minor, and it would be easiest with unblinded data. + +Lastly, we establish the parameters for the upper and lower bounds, +following the similar procedure as we have in the 1-sided design. # References From 3524fa285efc1f0504c18961489108560e6fa7a0 Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Fri, 16 Feb 2024 12:43:43 -0500 Subject: [PATCH 17/30] fix lintr issues --- vignettes/articles/story-update-boundary.Rmd | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 89d0391a..54dc81a7 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -145,8 +145,8 @@ The details of its arguments and implementations are explained in the Appendix. ```{r, echo=TRUE} x_updated <- gs_power_npe( - # `theta = 0` provides the crossing probability under H0. - # Users have the flexibility to modify the value of `theta, + # `theta = 0` provides the crossing probability under H0. + # Users have the flexibility to modify the value of `theta, # if they are interested in the crossing probability under H1 or other scenarios. theta = 0, # Observed statistical information under H0 with equal randomization @@ -355,14 +355,14 @@ time-to-event data is available, without the corresponding event counts. ### Option 1: `theta1` is the weighted average ```{r, echo=TRUE} x_updated <- gs_power_npe( - # `theta = 0` provides the crossing probability under H0. - # Users have the flexibility to modify the value of `theta, + # `theta = 0` provides the crossing probability under H0. + # Users have the flexibility to modify the value of `theta, # if they are interested in the crossing probability under H1 or other scenarios. theta = 0, # -log(AHR) for H0 theta0 = 0, # -log(AHR) for H1 - # It is used for determining the lower spending and bounds. + # It is used for determining the lower spending and bounds. theta1 = -(event_first_3_month * log(1) + event_after_3_month * log(0.6)) / event_observed[2], # Observed statistical information under H0 with equal randomization info = event_observed * ratio / (1 + ratio)^2, @@ -419,14 +419,14 @@ x_updated |> ### Option 2: `theta1` is based on blinded ahr estimation ```{r, echo=TRUE} x_updated <- gs_power_npe( - # `theta = 0` provides the crossing probability under H0. - # Users have the flexibility to modify the value of `theta, + # `theta = 0` provides the crossing probability under H0. + # Users have the flexibility to modify the value of `theta, # if they are interested in the crossing probability under H1 or other scenarios. theta = 0, # -log(AHR) for H0 theta0 = 0, # -log(AHR) for H1 - # It is used for determining the lower spending and bounds. + # It is used for determining the lower spending and bounds. theta1 = blinded_ahr$theta, # Observed statistical information under H0 with equal randomization info = event_observed * ratio / (1 + ratio)^2, From ab3ffb8ad6d19b96b51d1120c24875aca7463a4e Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Thu, 22 Feb 2024 15:52:38 -0500 Subject: [PATCH 18/30] add references --- vignettes/articles/story-update-boundary.Rmd | 1 + 1 file changed, 1 insertion(+) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 54dc81a7..33005bea 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -554,4 +554,5 @@ very minor, and it would be easiest with unblinded data. Lastly, we establish the parameters for the upper and lower bounds, following the similar procedure as we have in the 1-sided design. + # References From c75361a63aa6f738a9b001039ddc090fc70b2dfa Mon Sep 17 00:00:00 2001 From: keaven Date: Sat, 24 Feb 2024 12:54:13 -0500 Subject: [PATCH 19/30] Copyright update --- R/gs_info_wlr.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/gs_info_wlr.R b/R/gs_info_wlr.R index 776c4f2f..89c1b7a5 100644 --- a/R/gs_info_wlr.R +++ b/R/gs_info_wlr.R @@ -1,9 +1,9 @@ # Copyright (c) 2024 Merck & Co., Inc., Rahway, NJ, USA and its affiliates. # All rights reserved. # -# This file is part of the gsdmvn program. +# This file is part of the gsDesign2 program. # -# gsdmvn is free software: you can redistribute it and/or modify +# gsDesign2 is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. From a1a4f9d5a9760e0ad7d12a19eaf9588ed1b3a591 Mon Sep 17 00:00:00 2001 From: keaven Date: Sat, 2 Mar 2024 09:19:09 -0500 Subject: [PATCH 20/30] Reuse design input at time of design update --- vignettes/articles/story-update-boundary.Rmd | 39 ++++++++++---------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 33005bea..25836b4d 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -144,29 +144,28 @@ We will utilize the `gs_power_npe()` function to update efficacy bounds based on The details of its arguments and implementations are explained in the Appendix. ```{r, echo=TRUE} +# Take spending setup from original design +upar <- x$upper$upar +# Now update timing parameter for the interim analysis. +# Interim spending based on observed events divided by final planned events. +# The remaining alpha will be allocated to FA. +upar$timing <- c(event_observed[1] / max(event_planned), 1) + x_updated <- gs_power_npe( - # `theta = 0` provides the crossing probability under H0. - # Users have the flexibility to modify the value of `theta, - # if they are interested in the crossing probability under H1 or other scenarios. + # theta = 0 provides the crossing probability under H0. theta = 0, - # Observed statistical information under H0 with equal randomization - info = event_observed * ratio / (1 + ratio)^2, - # Upper bound - upper = gs_spending_bound, - upar = list( - sf = gsDesign::sfLDOF, - total_spend = alpha, - param = NULL, - # The remaining alpha will be allocated to the FA stage. - timing = c(event_observed[1] / max(event_planned), 1) - ), - test_upper = TRUE, - # Lower bound - lower = gs_b, - lpar = rep(-Inf, 2), - test_lower = FALSE, + # Observed statistical information under H0 + info = event_observed * x$input$ratio / (1 + x$input$ratio)^2, + # Upper bound uses spending function from planned design + upper = x$input$upper, + upar = upar, + test_upper = x$input$test_upper, + # No lower bound, but copy this from input + lower = x$input$lower, + lpar = x$input$lpar, + test_lower = x$input$test_lower, # Binding - binding = FALSE + binding = x$input$binding ) ``` From 576c5ff7d99c66548553329ff312ee033cc5e300 Mon Sep 17 00:00:00 2001 From: keaven Date: Thu, 14 Mar 2024 11:43:52 -0400 Subject: [PATCH 21/30] ongoing edits --- .../feature-comparison-gsDesign-gsDesign2.Rmd | 94 +++++++++++++++++++ 1 file changed, 94 insertions(+) create mode 100644 vignettes/articles/feature-comparison-gsDesign-gsDesign2.Rmd diff --git a/vignettes/articles/feature-comparison-gsDesign-gsDesign2.Rmd b/vignettes/articles/feature-comparison-gsDesign-gsDesign2.Rmd new file mode 100644 index 00000000..149efcc0 --- /dev/null +++ b/vignettes/articles/feature-comparison-gsDesign-gsDesign2.Rmd @@ -0,0 +1,94 @@ +--- +title: "gsDesign vs. gsDesign2 Feature Comparison" +output: + rmarkdown::html_document: + toc: true + toc_float: true + toc_depth: 2 + number_sections: true + highlight: "textmate" + css: "custom.css" +bibliography: "gsDesign2.bib" +vignette: > + %\VignetteEngine{knitr::rmarkdown} + %\VignetteIndexEntry{gsDesign vs. gsDesign2 Feature Comparison} +--- + + +\newcommand{\boxedcheckmark} + {{\ooalign{$\Box$\cr\hidewidth$\checkmark$\hidewidth}}} + + +$\boxedcheckmark$ + +$\Box$ + + + +# Heading 1 + +## Heading 2 + +### Heading 3 + +# Markdown syntax + +This is an **R Markdown** HTML vignette. Markdown is a simple _formatting syntax_ +for authoring HTML, PDF, and Word documents. + +$$x^n + y^n = z^n$$ + +When you click the _Knit_ button a document will be generated that includes both +content as well as the output of any embedded R code chunks within the document. + +Here is an inline equation $A = \pi r^2$. + +# Syntax highlighting + +This is inline `code`. + +This is a code block: + +```{r} +seq_len(10) +``` + +This is another code block, without evaluation: + + +# Tables + +```{r, echo=FALSE} +knitr::kable(head(mtcars)) +``` + + + +# Markdown tables + +You can use any online Markdown table generator to create these tables. + +| Feature | gsDesign | gsDesign2 | +| ------------- |:-------------: | -----: | +| Nonproportional hazards | ❌ | ✔️ | +| Allow skipping bound at an analysis | ❌ | ✔️ | +| Alternates to logrank for survival analysis | ❌ | ✔️ | +| Alternate spending approaches | ❌ | ✔️ | +| HR bounds for futility | ❌ | ✔️ | +| Stratified design for binomial | ❌ | ✔️ | + +# Markdown figures + +Without caption: + +```markdown +#![](example.png){width=50% class="center-block"} +``` + +With caption: + +```markdown +#![Insert caption here.](example.png){width=50%} +``` + +# References From 8441052b3883184991281c3235a556fd15ada2b1 Mon Sep 17 00:00:00 2001 From: keaven Date: Thu, 14 Mar 2024 17:06:30 -0400 Subject: [PATCH 22/30] ongoing updates --- vignettes/articles/story-update-boundary.Rmd | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 25836b4d..6f075c21 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -28,8 +28,8 @@ library(tibble) # Design assumptions -We begin with design assumptions. -We assume two analyses: an interim analysis (IA) final Analysis (FA). The IA is planned for month 20 after opening enrollment, followed by the FA at month 36. +We assume two analyses: an interim analysis (IA) and a final Analysis (FA). +The IA is planned 20 months after opening enrollment, followed by the FA at month 36. The planned enrollment period spans 14 months, with the first 2 months having an enrollment rate of 1/3 the final rate, the next 2 months with a rate of 2/3 of the final rate, and the final rate for the remaining 10 months. To obtain the targeted 90\% power, these rates will be multiplied by a constant. The control arm is assumed to follow an exponential distribution with a median of 9 months and the dropout rate is 0.0001 per month regardless of treatment group. @@ -63,6 +63,8 @@ For the design, we have efficacy bounds at both the IA and FA. We use the @lan19 ## Planned design +We use the + ```{r, echo=TRUE} upper <- gs_spending_bound upar <- list(sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL) @@ -73,6 +75,7 @@ x <- gs_design_ahr( alpha = alpha, beta = beta, info_frac = NULL, + info_scale = "h0_info", analysis_time = analysis_time, ratio = ratio, upper = gs_spending_bound, @@ -87,8 +90,8 @@ x <- gs_design_ahr( The planned design targets: - planned events: `r round(x$analysis$event, 0)` -- planned information fraction at interim and final analysis: `r round(x$analysis$info_frac, 4)` -- planned alpha spending: `r gsDesign::sfLDOF(0.025, x$analysis$info_frac)$spend` +- planned information fraction for interim and final analysis: `r round(x$analysis$info_frac, 4)` +- planned alpha spending: `r round(gsDesign::sfLDOF(0.025, x$analysis$info_frac)$spend, 4)` - planned efficacy bounds: `r round(x$bound$z[x$bound$bound == "upper"], 4)` We note that rounding up the final targeted events increases power slightly over the targeted 90%. @@ -96,7 +99,7 @@ We note that rounding up the final targeted events increases power slightly over ```{r} x |> summary() |> - gt::gt() |> + as_gt() |> gt::tab_header(title = "Original design") |> gt::tab_style( style = list( @@ -158,7 +161,7 @@ x_updated <- gs_power_npe( info = event_observed * x$input$ratio / (1 + x$input$ratio)^2, # Upper bound uses spending function from planned design upper = x$input$upper, - upar = upar, + upar = x$input$upar, test_upper = x$input$test_upper, # No lower bound, but copy this from input lower = x$input$lower, @@ -222,6 +225,7 @@ x <- gs_design_ahr( alpha = alpha, beta = beta, info_frac = NULL, + info_scale = "h0_info", analysis_time = c(20, 36), ratio = ratio, upper = gs_spending_bound, @@ -248,7 +252,7 @@ larger than what we have in the 1-sided example. ```{r} x |> summary() |> - gt::gt() |> + as_gt() |> gt::tab_header(title = "Original design") |> gt::tab_style( style = list( @@ -481,6 +485,7 @@ x_updated |> # Appendix + ## Arguments of `gs_power_npe()` We provide an introduction to its arguments. From 0181f59acd1b166de99eba992d41e4959c7968e9 Mon Sep 17 00:00:00 2001 From: keaven Date: Fri, 15 Mar 2024 19:02:38 -0400 Subject: [PATCH 23/30] near final? --- vignettes/articles/story-update-boundary.Rmd | 125 +++++++------------ 1 file changed, 43 insertions(+), 82 deletions(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 6f075c21..667f44a5 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -24,6 +24,7 @@ knitr::opts_chunk$set(echo = FALSE) library(gsDesign2) library(gt) library(tibble) +library(dplyr) ``` # Design assumptions @@ -57,13 +58,21 @@ analysis_time <- c(20, 36) ratio <- 1 ``` +We use the null hypothesis information for boundary crossing probability calculations under both the null and alternate hypotheses. +This will also imply the null hypothesis information will be used for the information fraction used in spending functions to derive the design. + +```{r} +info_scale = "h0_info" +``` + # One-sided design {.tabset} For the design, we have efficacy bounds at both the IA and FA. We use the @lan1983discrete spending function with a total alpha of `r alpha`, which approximates an O'Brien-Fleming bound. ## Planned design -We use the + + ```{r, echo=TRUE} upper <- gs_spending_bound @@ -159,6 +168,7 @@ x_updated <- gs_power_npe( theta = 0, # Observed statistical information under H0 info = event_observed * x$input$ratio / (1 + x$input$ratio)^2, + info_scale = "h0_info", # Upper bound uses spending function from planned design upper = x$input$upper, upar = x$input$upar, @@ -204,12 +214,12 @@ In this section, we investigate a 2 sided asymmetric design, with the non-binding beta-spending futility bounds. Beta-spending refers to error spending for the lower bound crossing probabilities under the alternative hypothesis. Non-binding assumes the trial continues if the -lower bound is crossed for Type I and Type II error computation. +lower bound is crossed for Type I, but not Type II error computation. ## Planned design -In the original designs, we employ the Lan DeMets O'Brien-Fleming -spending function [@lan1983discrete] for both efficacy and futility bounds. +In the original designs, we employ the Lan-DeMets spending function used to approximate O'Brien-Fleming +bounds [@lan1983discrete] for both efficacy and futility bounds. The total spending for efficacy is `r alpha`, and for futility is `r beta`. Besides, we assume the futility test only happens at IA. @@ -276,19 +286,21 @@ x |> ) ``` -## Update bounds at time of analysis {.tabset} +## Update bounds at time of analysis In practice, let us assume the observed data is generated by `simtrial::sim_pw_surv()`. ```{r, echo=TRUE} set.seed(2024) -observed_data <- simtrial::sim_pw_surv(n = x$analysis$n[x$analysis$analysis == 2], - stratum = data.frame(stratum = "All", p = 1), - block = c(rep("control", 2), rep("experimental", 2)), - enroll_rate = x$enroll_rate, - fail_rate = (fail_rate |> simtrial::to_sim_pw_surv())$fail_rate, - dropout_rate = (fail_rate |> simtrial::to_sim_pw_surv())$dropout_rate) +observed_data <- simtrial::sim_pw_surv( + n = x$analysis$n[x$analysis$analysis == 2], + stratum = data.frame(stratum = "All", p = 1), + block = c(rep("control", 2), rep("experimental", 2)), + enroll_rate = x$enroll_rate, + fail_rate = (fail_rate |> simtrial::to_sim_pw_surv())$fail_rate, + dropout_rate = (fail_rate |> simtrial::to_sim_pw_surv())$dropout_rate +) observed_ia_data <- observed_data |> simtrial::cut_data_by_date(analysis_time[1]) observed_fa_data <- observed_data |> simtrial::cut_data_by_date(analysis_time[2]) @@ -333,6 +345,13 @@ event_after_3_month We can derive `theta1` as -[`r event_first_3_month` / `r event_observed[2]` $\times \log(1) +$ `r event_after_3_month` / `r event_observed[2]` $\times \log(0.6)$] = `r round(-(event_first_3_month * log(1) + event_after_3_month * log(0.6)) / event_observed[2], 4)`. +```{r} +theta1 <- -(event_first_3_month * log(1) + event_after_3_month * log(0.6)) / event_observed[2] +theta1 +``` + + + - **Option 2:** we use the blinded data to estimate the AHR and `theta1` is the negative logarithm of the estimated AHR. ```{r, echo=TRUE} @@ -349,24 +368,27 @@ blinded_ahr |> gt::tab_header("Blinded estimation of average hazard ratio") ``` -In this example, the estimated AHR is `r round(blinded_ahr$ahr, 4)`, and accordingly, the `theta1` values is `r round(blinded_ahr$theta, 4)`. +```{r} +theta1 <- blinded_ahr$theta +theta1 +``` + -Both Option 1 and Option 2 yield the same value of theta1. Option 1 is applicable +Both Option 1 and Option 2 yield the same value of `theta1`. Option 1 is applicable when the event counts are available, while Option 2 is suitable for scenarios where time-to-event data is available, without the corresponding event counts. +A value of `theta1` is computed with one of the above methods and incorporated below. -### Option 1: `theta1` is the weighted average ```{r, echo=TRUE} x_updated <- gs_power_npe( # `theta = 0` provides the crossing probability under H0. - # Users have the flexibility to modify the value of `theta, + # Users have the flexibility to modify the value of `theta`, # if they are interested in the crossing probability under H1 or other scenarios. theta = 0, - # -log(AHR) for H0 + # -log(AHR) = 0 for H0 is used for determining upper spending theta0 = 0, - # -log(AHR) for H1 - # It is used for determining the lower spending and bounds. - theta1 = -(event_first_3_month * log(1) + event_after_3_month * log(0.6)) / event_observed[2], + # -log(AHR) for H1 is used for determining the lower spending and bounds. + theta1 = theta1, # Observed statistical information under H0 with equal randomization info = event_observed * ratio / (1 + ratio)^2, # Upper bound @@ -393,72 +415,11 @@ x_updated <- gs_power_npe( ) ``` -```{r} -x_updated |> - gt::gt() |> - gt::tab_header(title = "Updated design", subtitle = "by event fraction based timing") |> - gt::tab_style( - style = list( - gt::cell_fill(color = "lightcyan"), - gt::cell_text(weight = "bold") - ), - locations = gt::cells_body( - columns = probability, - rows = bound == "upper" - ) - ) |> - gt::tab_style( - style = list( - gt::cell_fill(color = "#F9E3D6"), - gt::cell_text(weight = "bold") - ), - locations = gt::cells_body( - columns = probability, - rows = bound == "lower" - ) - ) -``` -### Option 2: `theta1` is based on blinded ahr estimation -```{r, echo=TRUE} -x_updated <- gs_power_npe( - # `theta = 0` provides the crossing probability under H0. - # Users have the flexibility to modify the value of `theta, - # if they are interested in the crossing probability under H1 or other scenarios. - theta = 0, - # -log(AHR) for H0 - theta0 = 0, - # -log(AHR) for H1 - # It is used for determining the lower spending and bounds. - theta1 = blinded_ahr$theta, - # Observed statistical information under H0 with equal randomization - info = event_observed * ratio / (1 + ratio)^2, - # Upper bound - upper = x$input$upper, - upar = list( - sf = gsDesign::sfLDOF, - total_spend = alpha, - param = NULL, - # The remaining alpha will be allocated to the FA stage. - timing = c(event_observed[1] / max(event_planned), 1) - ), - test_upper = TRUE, - # Lower bound - lower = x$input$lower, - lpar = list( - sf = gsDesign::sfLDOF, - total_spend = beta, - param = NULL, - # The remaining alpha will be allocated to the FA stage. - timing = c(event_observed[1] / max(event_planned), 1) - ), - test_lower = c(TRUE, FALSE), - binding = x$input$binding -) -``` ```{r} -x_updated |> +x_updated |> + dplyr::filter(abs(z) < Inf) |> gt::gt() |> gt::tab_header(title = "Updated design", subtitle = "by event fraction based timing") |> gt::tab_style( From dc85a73cb662aa5fb0cf8094d37a038008444cb0 Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Mon, 18 Mar 2024 11:09:30 -0400 Subject: [PATCH 24/30] 2 lintr issues --- vignettes/articles/story-update-boundary.Rmd | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 667f44a5..23fa98ba 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -62,7 +62,7 @@ We use the null hypothesis information for boundary crossing probability calcula This will also imply the null hypothesis information will be used for the information fraction used in spending functions to derive the design. ```{r} -info_scale = "h0_info" +info_scale <- "h0_info" ``` # One-sided design {.tabset} @@ -418,7 +418,7 @@ x_updated <- gs_power_npe( ```{r} -x_updated |> +x_updated |> dplyr::filter(abs(z) < Inf) |> gt::gt() |> gt::tab_header(title = "Updated design", subtitle = "by event fraction based timing") |> From 3e5ece8744c0726083d0d5fc6c2dab3352c6860d Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Mon, 18 Mar 2024 11:16:44 -0400 Subject: [PATCH 25/30] delete a redundant vignette --- .../feature-comparison-gsDesign-gsDesign2.Rmd | 94 ------------------- 1 file changed, 94 deletions(-) delete mode 100644 vignettes/articles/feature-comparison-gsDesign-gsDesign2.Rmd diff --git a/vignettes/articles/feature-comparison-gsDesign-gsDesign2.Rmd b/vignettes/articles/feature-comparison-gsDesign-gsDesign2.Rmd deleted file mode 100644 index 149efcc0..00000000 --- a/vignettes/articles/feature-comparison-gsDesign-gsDesign2.Rmd +++ /dev/null @@ -1,94 +0,0 @@ ---- -title: "gsDesign vs. gsDesign2 Feature Comparison" -output: - rmarkdown::html_document: - toc: true - toc_float: true - toc_depth: 2 - number_sections: true - highlight: "textmate" - css: "custom.css" -bibliography: "gsDesign2.bib" -vignette: > - %\VignetteEngine{knitr::rmarkdown} - %\VignetteIndexEntry{gsDesign vs. gsDesign2 Feature Comparison} ---- - - -\newcommand{\boxedcheckmark} - {{\ooalign{$\Box$\cr\hidewidth$\checkmark$\hidewidth}}} - - -$\boxedcheckmark$ - -$\Box$ - - - -# Heading 1 - -## Heading 2 - -### Heading 3 - -# Markdown syntax - -This is an **R Markdown** HTML vignette. Markdown is a simple _formatting syntax_ -for authoring HTML, PDF, and Word documents. - -$$x^n + y^n = z^n$$ - -When you click the _Knit_ button a document will be generated that includes both -content as well as the output of any embedded R code chunks within the document. - -Here is an inline equation $A = \pi r^2$. - -# Syntax highlighting - -This is inline `code`. - -This is a code block: - -```{r} -seq_len(10) -``` - -This is another code block, without evaluation: - - -# Tables - -```{r, echo=FALSE} -knitr::kable(head(mtcars)) -``` - - - -# Markdown tables - -You can use any online Markdown table generator to create these tables. - -| Feature | gsDesign | gsDesign2 | -| ------------- |:-------------: | -----: | -| Nonproportional hazards | ❌ | ✔️ | -| Allow skipping bound at an analysis | ❌ | ✔️ | -| Alternates to logrank for survival analysis | ❌ | ✔️ | -| Alternate spending approaches | ❌ | ✔️ | -| HR bounds for futility | ❌ | ✔️ | -| Stratified design for binomial | ❌ | ✔️ | - -# Markdown figures - -Without caption: - -```markdown -#![](example.png){width=50% class="center-block"} -``` - -With caption: - -```markdown -#![Insert caption here.](example.png){width=50%} -``` - -# References From 9806efa18b855ae457c318281fc6cbfc6a1b487e Mon Sep 17 00:00:00 2001 From: Nan Xiao Date: Mon, 18 Mar 2024 12:17:30 -0400 Subject: [PATCH 26/30] Style `ahr_blinded()` documentation and code --- R/ahr_blinded.R | 93 +++++++++++++++++++++++++--------------------- man/ahr_blinded.Rd | 56 ++++++++++++++-------------- 2 files changed, 80 insertions(+), 69 deletions(-) diff --git a/R/ahr_blinded.R b/R/ahr_blinded.R index 7f601da5..b4cf7780 100644 --- a/R/ahr_blinded.R +++ b/R/ahr_blinded.R @@ -1,3 +1,6 @@ +# Copyright (c) 2024 Merck & Co., Inc., Rahway, NJ, USA and its affiliates. +# All rights reserved. +# # This file is part of the gsDesign2 program. # # gsDesign2 is free software: you can redistribute it and/or modify @@ -15,49 +18,55 @@ #' Blinded estimation of average hazard ratio #' -#' Based on blinded data and assumed hazard ratios in different intervals, compute -#' a blinded estimate of average hazard ratio (AHR) and corresponding estimate of statistical information. -#' This function is intended for use in computing futility bounds based on spending assuming -#' the input hazard ratio (hr) values for intervals specified here. -#' @importFrom tibble tibble -#' @importFrom survival Surv -#' @param surv input survival object (see \code{Surv}); note that only 0=censored, 1=event for \code{Surv} -#' @param intervals Vector containing positive values indicating interval lengths where the -#' exponential rates are assumed. -#' Note that a final infinite interval is added if any events occur after the final interval -#' specified. -#' @param hr vector of hazard ratios assumed for each interval -#' @param ratio ratio of experimental to control randomization. +#' Based on blinded data and assumed hazard ratios in different intervals, +#' compute a blinded estimate of average hazard ratio (AHR) and corresponding +#' estimate of statistical information. +#' This function is intended for use in computing futility bounds based on +#' spending assuming the input hazard ratio (hr) values for intervals +#' specified here. +#' +#' @param surv Input survival object (see [survival::Surv()]); +#' note that only 0 = censored, 1 = event for [survival::Surv()]. +#' @param intervals Vector containing positive values indicating +#' interval lengths where the exponential rates are assumed. +#' Note that a final infinite interval is added if any events occur +#' after the final interval specified. +#' @param hr Vector of hazard ratios assumed for each interval. +#' @param ratio Ratio of experimental to control randomization. +#' +#' @return A `tibble` with one row containing +#' - `ahr` - Blinded average hazard ratio based on assumed period-specific +#' hazard ratios input in `fail_rate` and observed events in the +#' corresponding intervals. +#' - `event` - Total observed number of events. +#' - `info0` - Information under related null hypothesis. #' #' @section Specification: #' \if{latex}{ #' \itemize{ #' \item Validate if input hr is a numeric vector. #' \item Validate if input hr is non-negative. -#' \item Simulate piece-wise exponential survival estimation with the inputs survival object surv -#' and intervals. -#' \item Save the length of hr and events to an object, and if the length of hr is shorter than -#' the intervals, add replicates of the last element of hr and the corresponding numbers of events -#' to hr. +#' \item Simulate piece-wise exponential survival estimation with the inputs +#' survival object surv and intervals. +#' \item Save the length of hr and events to an object, and if the length of +#' hr is shorter than the intervals, add replicates of the last element of hr +#' and the corresponding numbers of events to hr. #' \item Compute the blinded estimation of average hazard ratio. #' \item Compute adjustment for information. -#' \item Return a tibble of the sum of events, average hazard raito, blinded average hazard -#' ratio, and the information. +#' \item Return a tibble of the sum of events, average hazard ratio, +#' blinded average hazard ratio, and the information. #' } #' } #' \if{html}{The contents of this section are shown in PDF user manual only.} #' -#' @return A \code{tibble} with one row containing -#' `ahr` blinded average hazard ratio based on assumed period-specific hazard ratios input in `fail_rate` -#' and observed events in the corresponding intervals -#' `event` total observed number of events, and `info0` (information under related null hypothesis). +#' @importFrom tibble tibble +#' @importFrom survival Surv +#' +#' @export #' #' @examples -#' \donttest{ -#' library(simtrial) -#' library(survival) #' ahr_blinded( -#' surv = Surv( +#' surv = survival::Surv( #' time = simtrial::ex2_delayed_effect$month, #' event = simtrial::ex2_delayed_effect$evntd #' ), @@ -65,14 +74,14 @@ #' hr = c(1, .55), #' ratio = 1 #' ) -#' } -#' -#' @export -ahr_blinded <- function(surv = survival::Surv(time = simtrial::ex1_delayed_effect$month, - event = simtrial::ex1_delayed_effect$evntd), - intervals = c(3, Inf), - hr = c(1, .6), - ratio = 1) { +ahr_blinded <- function( + surv = survival::Surv( + time = simtrial::ex1_delayed_effect$month, + event = simtrial::ex1_delayed_effect$evntd + ), + intervals = c(3, Inf), + hr = c(1, .6), + ratio = 1) { # Input checking if (!is.vector(hr, mode = "numeric") || min(hr) <= 0) { stop("ahr_blinded: hr must be a vector of positive numbers.") @@ -94,19 +103,19 @@ ahr_blinded <- function(surv = survival::Surv(time = simtrial::ex1_delayed_effec nx <- length(event) # Add to hr if length shorter than intervals - if (length(hr) < length(event)) { - hr <- c(hr, rep(hr[nhr], nx - nhr)) - } + if (length(hr) < length(event)) hr <- c(hr, rep(hr[nhr], nx - nhr)) # Compute blinded AHR - theta <- - sum(log(hr[1:nx]) * event) / sum(event) + theta <- -sum(log(hr[1:nx]) * event) / sum(event) # Compute adjustment for information q_e <- ratio / (1 + ratio) ans <- tibble( - event = sum(event), ahr = exp(-theta), - theta = theta, info0 = sum(event) * (1 - q_e) * q_e + event = sum(event), + ahr = exp(-theta), + theta = theta, + info0 = sum(event) * (1 - q_e) * q_e ) return(ans) } diff --git a/man/ahr_blinded.Rd b/man/ahr_blinded.Rd index 419e6987..a23c01f1 100644 --- a/man/ahr_blinded.Rd +++ b/man/ahr_blinded.Rd @@ -13,28 +13,35 @@ ahr_blinded( ) } \arguments{ -\item{surv}{input survival object (see \code{Surv}); note that only 0=censored, 1=event for \code{Surv}} +\item{surv}{Input survival object (see \code{\link[survival:Surv]{survival::Surv()}}); +note that only 0 = censored, 1 = event for \code{\link[survival:Surv]{survival::Surv()}}.} -\item{intervals}{Vector containing positive values indicating interval lengths where the -exponential rates are assumed. -Note that a final infinite interval is added if any events occur after the final interval -specified.} +\item{intervals}{Vector containing positive values indicating +interval lengths where the exponential rates are assumed. +Note that a final infinite interval is added if any events occur +after the final interval specified.} -\item{hr}{vector of hazard ratios assumed for each interval} +\item{hr}{Vector of hazard ratios assumed for each interval.} -\item{ratio}{ratio of experimental to control randomization.} +\item{ratio}{Ratio of experimental to control randomization.} } \value{ A \code{tibble} with one row containing -\code{ahr} blinded average hazard ratio based on assumed period-specific hazard ratios input in \code{fail_rate} -and observed events in the corresponding intervals -\code{event} total observed number of events, and \code{info0} (information under related null hypothesis). +\itemize{ +\item \code{ahr} - Blinded average hazard ratio based on assumed period-specific +hazard ratios input in \code{fail_rate} and observed events in the +corresponding intervals. +\item \code{event} - Total observed number of events. +\item \code{info0} - Information under related null hypothesis. +} } \description{ -Based on blinded data and assumed hazard ratios in different intervals, compute -a blinded estimate of average hazard ratio (AHR) and corresponding estimate of statistical information. -This function is intended for use in computing futility bounds based on spending assuming -the input hazard ratio (hr) values for intervals specified here. +Based on blinded data and assumed hazard ratios in different intervals, +compute a blinded estimate of average hazard ratio (AHR) and corresponding +estimate of statistical information. +This function is intended for use in computing futility bounds based on +spending assuming the input hazard ratio (hr) values for intervals +specified here. } \section{Specification}{ @@ -42,26 +49,23 @@ the input hazard ratio (hr) values for intervals specified here. \itemize{ \item Validate if input hr is a numeric vector. \item Validate if input hr is non-negative. - \item Simulate piece-wise exponential survival estimation with the inputs survival object surv - and intervals. - \item Save the length of hr and events to an object, and if the length of hr is shorter than - the intervals, add replicates of the last element of hr and the corresponding numbers of events - to hr. + \item Simulate piece-wise exponential survival estimation with the inputs + survival object surv and intervals. + \item Save the length of hr and events to an object, and if the length of + hr is shorter than the intervals, add replicates of the last element of hr + and the corresponding numbers of events to hr. \item Compute the blinded estimation of average hazard ratio. \item Compute adjustment for information. - \item Return a tibble of the sum of events, average hazard raito, blinded average hazard - ratio, and the information. + \item Return a tibble of the sum of events, average hazard ratio, + blinded average hazard ratio, and the information. } } \if{html}{The contents of this section are shown in PDF user manual only.} } \examples{ -\donttest{ -library(simtrial) -library(survival) ahr_blinded( - surv = Surv( + surv = survival::Surv( time = simtrial::ex2_delayed_effect$month, event = simtrial::ex2_delayed_effect$evntd ), @@ -70,5 +74,3 @@ ahr_blinded( ratio = 1 ) } - -} From f8e9dd4cb3d27262f8b4bc092fcbc03c27cb38ef Mon Sep 17 00:00:00 2001 From: Nan Xiao Date: Mon, 18 Mar 2024 12:20:49 -0400 Subject: [PATCH 27/30] Fix check note on no visible binding for global variable `status` --- R/globals.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/globals.R b/R/globals.R index 0e1e6b44..30fb57a3 100644 --- a/R/globals.R +++ b/R/globals.R @@ -23,6 +23,8 @@ utils::globalVariables( ".", ".SD", # From `ahr()` c("stratum", "rate", "hr", "treatment", "time", "info0", "info"), + # From `ahr_blinded()` + c("status"), # From `as_gt.gs_design()` c("Bound", "Alternate hypothesis", "Null hypothesis", "Analysis"), # From `expected_accrual()` From 121485ae8c8a205d1740f6c83924fcb170629348 Mon Sep 17 00:00:00 2001 From: Nan Xiao Date: Mon, 18 Mar 2024 12:24:51 -0400 Subject: [PATCH 28/30] Sort dependencies alphabetically --- DESCRIPTION | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 1966dbd4..2ad68f4a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -46,11 +46,11 @@ Imports: npsurvSS, r2rtf, stats, + survival, tibble, tidyr, utils, - Rcpp, - survival + Rcpp Suggests: covr, ggplot2, From d68e4de2cb4ea0e3dbf17aed46453eccc9c31ddd Mon Sep 17 00:00:00 2001 From: Nan Xiao Date: Mon, 18 Mar 2024 12:58:24 -0400 Subject: [PATCH 29/30] Remove unnecessary `library()` calls to avoid masking and improve text style --- vignettes/articles/story-update-boundary.Rmd | 213 +++++++++++-------- 1 file changed, 123 insertions(+), 90 deletions(-) diff --git a/vignettes/articles/story-update-boundary.Rmd b/vignettes/articles/story-update-boundary.Rmd index 23fa98ba..cda953af 100644 --- a/vignettes/articles/story-update-boundary.Rmd +++ b/vignettes/articles/story-update-boundary.Rmd @@ -22,44 +22,52 @@ knitr::opts_chunk$set(echo = FALSE) ```{r, message=FALSE, warning=FALSE} library(gsDesign2) -library(gt) -library(tibble) -library(dplyr) ``` # Design assumptions -We assume two analyses: an interim analysis (IA) and a final Analysis (FA). -The IA is planned 20 months after opening enrollment, followed by the FA at month 36. -The planned enrollment period spans 14 months, with the first 2 months having an enrollment rate of 1/3 the final rate, the next 2 months with a rate of 2/3 of the final rate, and the final rate for the remaining 10 months. -To obtain the targeted 90\% power, these rates will be multiplied by a constant. -The control arm is assumed to follow an exponential distribution with a median of 9 months and the dropout rate is 0.0001 per month regardless of treatment group. -Finally, the experimental treatment group is piecewise exponential with a 3-month delayed treatment effect; i.e. in the first 3 months HR = 1 and the HR is 0.6 thereafter. - +We assume two analyses: an interim analysis (IA) and a final analysis (FA). +The IA is planned 20 months after opening enrollment, followed by the FA at +month 36. +The planned enrollment period spans 14 months, with the first 2 months having +an enrollment rate of 1/3 the final rate, the next 2 months with a rate of 2/3 +of the final rate, and the final rate for the remaining 10 months. +To obtain the targeted 90\% power, these rates will be multiplied by a constant. +The control arm is assumed to follow an exponential distribution with a median +of 9 months and the dropout rate is 0.0001 per month regardless of treatment group. +Finally, the experimental treatment group is piecewise exponential with a +3-month delayed treatment effect; that is, in the first 3 months HR = 1 and +the HR is 0.6 thereafter. ```{r} alpha <- 0.025 beta <- 0.1 ratio <- 1 -# enrollment -enroll_rate <- define_enroll_rate(duration = c(2, 2, 10), - rate = (1:3) / 3) +# Enrollment +enroll_rate <- define_enroll_rate( + duration = c(2, 2, 10), + rate = (1:3) / 3 +) -# failure and dropout -fail_rate <- define_fail_rate(duration = c(3, Inf), - fail_rate = log(2) / 9, - hr = c(1, 0.6), - dropout_rate = .0001) +# Failure and dropout +fail_rate <- define_fail_rate( + duration = c(3, Inf), + fail_rate = log(2) / 9, + hr = c(1, 0.6), + dropout_rate = .0001 +) # IA and FA analysis time analysis_time <- c(20, 36) -# randomization ratio +# Randomization ratio ratio <- 1 ``` -We use the null hypothesis information for boundary crossing probability calculations under both the null and alternate hypotheses. -This will also imply the null hypothesis information will be used for the information fraction used in spending functions to derive the design. +We use the null hypothesis information for boundary crossing probability +calculations under both the null and alternate hypotheses. +This will also imply the null hypothesis information will be used for the +information fraction used in spending functions to derive the design. ```{r} info_scale <- "h0_info" @@ -67,13 +75,12 @@ info_scale <- "h0_info" # One-sided design {.tabset} -For the design, we have efficacy bounds at both the IA and FA. We use the @lan1983discrete spending function with a total alpha of `r alpha`, which approximates an O'Brien-Fleming bound. +For the design, we have efficacy bounds at both the IA and FA. +We use the @lan1983discrete spending function with a total alpha of `r alpha`, +which approximates an O'Brien-Fleming bound. ## Planned design - - - ```{r, echo=TRUE} upper <- gs_spending_bound upar <- list(sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL) @@ -96,14 +103,15 @@ x <- gs_design_ahr( ) |> to_integer() ``` -The planned design targets: +The planned design targets: -- planned events: `r round(x$analysis$event, 0)` -- planned information fraction for interim and final analysis: `r round(x$analysis$info_frac, 4)` -- planned alpha spending: `r round(gsDesign::sfLDOF(0.025, x$analysis$info_frac)$spend, 4)` -- planned efficacy bounds: `r round(x$bound$z[x$bound$bound == "upper"], 4)` - -We note that rounding up the final targeted events increases power slightly over the targeted 90%. +- Planned events: `r round(x$analysis$event, 0)` +- Planned information fraction for interim and final analysis: `r round(x$analysis$info_frac, 4)` +- Planned alpha spending: `r round(gsDesign::sfLDOF(0.025, x$analysis$info_frac)$spend, 4)` +- Planned efficacy bounds: `r round(x$bound$z[x$bound$bound == "upper"], 4)` + +We note that rounding up the final targeted events increases power slightly +over the targeted 90\%. ```{r} x |> @@ -134,9 +142,11 @@ x |> ## Update bounds at time of analysis -We assume 180 and 280 events observed at the IA and FA, respectively. +We assume 180 and 280 events observed at the IA and FA, respectively. We will assume the differences from planned are due to logistical considerations. -We also assume the protocol specifies that the full $\alpha$ will be spent at the final analysis even in a case like this where there is a shortfall of events versus the design plan. +We also assume the protocol specifies that the full $\alpha$ will be spent at +the final analysis even in a case like this where there is a shortfall of events +versus the design plan. ```{r, echo=TRUE} # Observed vs. planned events @@ -145,14 +155,17 @@ event_planned <- x$analysis$event ``` ```{r} -tibble(Analysis = c("IA", "FA"), - `Planned events` = event_planned, - `Observed events` = event_observed) |> - gt() |> - tab_header("Planned vs. Observed events") +tibble::tibble( + Analysis = c("IA", "FA"), + `Planned events` = event_planned, + `Observed events` = event_observed +) |> + gt::gt() |> + gt::tab_header("Planned vs. Observed events") ``` -We will utilize the `gs_power_npe()` function to update efficacy bounds based on the observed events. +We will utilize the `gs_power_npe()` function to update efficacy bounds +based on the observed events. The details of its arguments and implementations are explained in the Appendix. ```{r, echo=TRUE} @@ -164,7 +177,7 @@ upar <- x$upper$upar upar$timing <- c(event_observed[1] / max(event_planned), 1) x_updated <- gs_power_npe( - # theta = 0 provides the crossing probability under H0. + # `theta = 0` provides the crossing probability under H0 theta = 0, # Observed statistical information under H0 info = event_observed * x$input$ratio / (1 + x$input$ratio)^2, @@ -218,8 +231,9 @@ lower bound is crossed for Type I, but not Type II error computation. ## Planned design -In the original designs, we employ the Lan-DeMets spending function used to approximate O'Brien-Fleming -bounds [@lan1983discrete] for both efficacy and futility bounds. +In the original designs, we employ the Lan-DeMets spending function used to +approximate O'Brien-Fleming bounds [@lan1983discrete] for both efficacy and +futility bounds. The total spending for efficacy is `r alpha`, and for futility is `r beta`. Besides, we assume the futility test only happens at IA. @@ -248,13 +262,13 @@ x <- gs_design_ahr( ) |> to_integer() ``` -In the planned design, we have +In the planned design, we have -- planned events: `r round(x$analysis$event, 0)` -- planned information fraction (timing): `r round(x$analysis$info_frac, 4)` -- planned alpha spending: `r gsDesign::sfLDOF(0.025, x$analysis$info_frac)$spend` -- planned efficacy bounds: `r round(x$bound$z[x$bound$bound == "upper"], 4)` -- planned futility bounds: `r round(x$bound$z[x$bound$bound == "lower"], 4)` +- Planned events: `r round(x$analysis$event, 0)` +- Planned information fraction (timing): `r round(x$analysis$info_frac, 4)` +- Planned alpha spending: `r gsDesign::sfLDOF(0.025, x$analysis$info_frac)$spend` +- Planned efficacy bounds: `r round(x$bound$z[x$bound$bound == "upper"], 4)` +- Planned futility bounds: `r round(x$bound$z[x$bound$bound == "lower"], 4)` Since we added futility bounds, the sample size and number of events are larger than what we have in the 1-sided example. @@ -288,10 +302,10 @@ x |> ## Update bounds at time of analysis -In practice, let us assume the observed data is generated by `simtrial::sim_pw_surv()`. +In practice, let us assume the observed data is generated by `simtrial::sim_pw_surv()`. ```{r, echo=TRUE} -set.seed(2024) +set.seed(42) observed_data <- simtrial::sim_pw_surv( n = x$analysis$n[x$analysis$analysis == 2], @@ -320,28 +334,36 @@ tibble::tibble( gt::tab_header("Planned vs. Observed events") ``` -Again, we use `gs_power_npe()` to calculate the updated efficacy and futility bounds. +Again, we use `gs_power_npe()` to calculate the updated efficacy and futility bounds. The details of its arguments and implementations are explained in the Appendix. -We initially set up `theta = 0` for crossing probability under the null hypothesis -and `theta0 = 0` for null hypothesis. +We initially set up `theta = 0` for crossing probability under the null hypothesis +and `theta0 = 0` for null hypothesis. As for `theta1`, there are 2 options to set it up. -These 2 options arrive at the same value of `theta1`. +These 2 options arrive at the same value of `theta1`. + +- **Option 1:** `theta1` is the weighted average of the piecewise HR, that is, -- **Option 1:** `theta1` is the weighted average of the piecewise HR, i.e., $$ -\log(\text{AHR}) = -\log(\sum_{i=1}^m w_i \; \text{HR}_m), $$ -where the weight is decided by the ratio of observed number of events and the observed final events. -In this example, the number of observed events at the first 3 months (HR = 1) can be derived as + +where the weight is decided by the ratio of observed number of events and the +observed final events. +In this example, the number of observed events at the first 3 months (HR = 1) +can be derived as + ```{r, echo=TRUE} event_first_3_month <- sum(observed_data$fail_time < 3) event_first_3_month ``` + and the number of observed events after 3 months (HR = 0.6) are + ```{r, echo=TRUE} event_after_3_month <- event_observed[2] - event_first_3_month event_after_3_month ``` + We can derive `theta1` as -[`r event_first_3_month` / `r event_observed[2]` $\times \log(1) +$ `r event_after_3_month` / `r event_observed[2]` $\times \log(0.6)$] = `r round(-(event_first_3_month * log(1) + event_after_3_month * log(0.6)) / event_observed[2], 4)`. @@ -350,16 +372,19 @@ theta1 <- -(event_first_3_month * log(1) + event_after_3_month * log(0.6)) / eve theta1 ``` - - -- **Option 2:** we use the blinded data to estimate the AHR and `theta1` is the negative logarithm of the estimated AHR. +- **Option 2:** we use the blinded data to estimate the AHR and `theta1` is the + negative logarithm of the estimated AHR. ```{r, echo=TRUE} -blinded_ahr <- ahr_blinded(surv = survival::Surv(time = observed_fa_data$tte, - event = observed_fa_data$event), - intervals = c(3, Inf), - hr = c(1, 0.6), - ratio = 1) +blinded_ahr <- ahr_blinded( + surv = survival::Surv( + time = observed_fa_data$tte, + event = observed_fa_data$event + ), + intervals = c(3, Inf), + hr = c(1, 0.6), + ratio = 1 +) ``` ```{r} @@ -373,9 +398,8 @@ theta1 <- blinded_ahr$theta theta1 ``` - -Both Option 1 and Option 2 yield the same value of `theta1`. Option 1 is applicable -when the event counts are available, while Option 2 is suitable for scenarios where +Both Option 1 and Option 2 yield the same value of `theta1`. Option 1 is applicable +when the event counts are available, while Option 2 is suitable for scenarios where time-to-event data is available, without the corresponding event counts. A value of `theta1` is computed with one of the above methods and incorporated below. @@ -387,7 +411,7 @@ x_updated <- gs_power_npe( theta = 0, # -log(AHR) = 0 for H0 is used for determining upper spending theta0 = 0, - # -log(AHR) for H1 is used for determining the lower spending and bounds. + # -log(AHR) for H1 is used for determining the lower spending and bounds theta1 = theta1, # Observed statistical information under H0 with equal randomization info = event_observed * ratio / (1 + ratio)^2, @@ -397,7 +421,7 @@ x_updated <- gs_power_npe( sf = gsDesign::sfLDOF, total_spend = alpha, param = NULL, - # The remaining alpha will be allocated to the FA stage. + # The remaining alpha will be allocated to the FA stage timing = c(event_observed[1] / max(event_planned), 1) ), test_upper = TRUE, @@ -407,7 +431,7 @@ x_updated <- gs_power_npe( sf = gsDesign::sfLDOF, total_spend = beta, param = NULL, - # The remaining alpha will be allocated to the FA stage. + # The remaining alpha will be allocated to the FA stage timing = c(event_observed[1] / max(event_planned), 1) ), test_lower = c(TRUE, FALSE), @@ -415,8 +439,6 @@ x_updated <- gs_power_npe( ) ``` - - ```{r} x_updated |> dplyr::filter(abs(z) < Inf) |> @@ -446,27 +468,36 @@ x_updated |> # Appendix - ## Arguments of `gs_power_npe()` We provide an introduction to its arguments. - *3 theta-related arguments*: - + `theta`: a vector with the natural parameter for the group sequential design; represents expected ; $-\log(\text{AHR})$ in this case. - This is used for boundary crossing probability calculation. - For example, if we set `theta = 0`, then the crossing probability is the Type I error. - + `theta0`: natural parameter for H0; used for determining the upper (efficacy) spending and bounds. If the default of `NULL` is given, then this is replaced by 0. - + `theta1`: natural parameter for H1; used for determining the lower spending and bounds. The default is `NULL`, in which case the value is replaced with the input `theta`. - + - `theta`: a vector with the natural parameter for the group sequential design; + represents expected ; $-\log(\text{AHR})$ in this case. + This is used for boundary crossing probability calculation. + For example, if we set `theta = 0`, then the crossing probability is the + Type I error. + - `theta0`: natural parameter for H0; used for determining the upper + (efficacy) spending and bounds. If the default of `NULL` is given, + then this is replaced by 0. + - `theta1`: natural parameter for H1; used for determining the lower spending + and bounds. The default is `NULL`, in which case the value is replaced with + the input `theta`. + - *3 statistical information-related arguments*: - + `info`: statistical information at all analyses for input `theta`. Default is `1`. - + `info0`: statistical information under H0, if different than `info`. It impacts null hypothesis bound calculation. Default is `NULL`. - + `info1`: statistical information under hypothesis used for futility bound calculation if different from `info`. It impacts futility hypothesis bound calculation. Default is `NULL`. - + - `info`: statistical information at all analyses for input `theta`. + Default is `1`. + - `info0`: statistical information under H0, if different than `info`. + It impacts null hypothesis bound calculation. Default is `NULL`. + - `info1`: statistical information under hypothesis used for futility bound + calculation if different from `info`. It impacts futility hypothesis bound + calculation. Default is `NULL`. + - *Efficacy/Futility boundary-related arguments*: - + `upper`, `upar` and `test_upper`: efficacy bounds related, same as `gs_design_ahr()`. - + `lower`, `lpar` and `test_lower`: futility bounds related, same as `gs_design_ahr()`. - + `binding`: `TRUE` or `FALSE` indicating whether it is binding or non-binding designs. + - `upper`, `upar` and `test_upper`: efficacy bounds related, same as `gs_design_ahr()`. + - `lower`, `lpar` and `test_lower`: futility bounds related, same as `gs_design_ahr()`. + - `binding`: `TRUE` or `FALSE` indicating whether it is binding or non-binding designs. ## Explanations of `gs_power_npe()` arguments set up in the one-sided design @@ -483,9 +514,11 @@ Moving forward, we proceed to set up the values of the `info` arguments for the input `theta = 0`. Since the statistical information is event-based, and `theta = 0` (null hypothesis), the observed statistical information under null hypothesis is + $$ \text{observed number of events} \times \frac{r}{1 + r} \times \frac{1}{1 + r}, $$ + where $r$ is the randomization ratio (experimental : control). In this example, the observed statistical information is calculated as `info = event_observed / 4`, and thus, we input `r event_observed / 4` @@ -498,7 +531,7 @@ while noting that `info1` does affect the futility bounds. Lastly, we establish the parameters for the upper and lower bounds. The setup of `upper`, `upar`, `test_upper`, `lower`, `lpar`, and `test_lower` closely resembles the original design, with the exception of the `upar` setup. -Here, we have introduced the `timing = ...` parameter. +Here, we have introduced the `timing` parameter. At the IA, the observed number of events is `r event_observed[1]`. If we employ the interim analysis alpha spending strategy, the spending timing at IA is calculated as `r event_observed[1]` / `r round(event_planned[2], 0)`. @@ -517,7 +550,7 @@ with the local null hypothesis approach. Though people can explicitly write the formula of `info1` in terms of the observed events, the difference is very minor, and it would be easiest with unblinded data. -Lastly, we establish the parameters for the upper and lower bounds, +Lastly, we establish the parameters for the upper and lower bounds, following the similar procedure as we have in the 1-sided design. # References From 00ca543af328057ec90f6fb53ba80a92993b55f3 Mon Sep 17 00:00:00 2001 From: "Zhao, Yujie" Date: Mon, 18 Mar 2024 13:10:05 -0400 Subject: [PATCH 30/30] add theta's return details --- R/ahr_blinded.R | 2 ++ man/ahr_blinded.Rd | 2 ++ 2 files changed, 4 insertions(+) diff --git a/R/ahr_blinded.R b/R/ahr_blinded.R index b4cf7780..11ba4dbe 100644 --- a/R/ahr_blinded.R +++ b/R/ahr_blinded.R @@ -40,6 +40,8 @@ #' corresponding intervals. #' - `event` - Total observed number of events. #' - `info0` - Information under related null hypothesis. +#' - `theta` - Natural parameter for group sequential design representing +#' expected incremental drift at all analyses. #' #' @section Specification: #' \if{latex}{ diff --git a/man/ahr_blinded.Rd b/man/ahr_blinded.Rd index a23c01f1..77391653 100644 --- a/man/ahr_blinded.Rd +++ b/man/ahr_blinded.Rd @@ -33,6 +33,8 @@ hazard ratios input in \code{fail_rate} and observed events in the corresponding intervals. \item \code{event} - Total observed number of events. \item \code{info0} - Information under related null hypothesis. +\item \code{theta} - Natural parameter for group sequential design representing +expected incremental drift at all analyses. } } \description{