From db1ad9c918ff507d5ab1cd5247686e36d66e5bb5 Mon Sep 17 00:00:00 2001
From: "Zhao, Yujie" <yujie.zhao@merck.com>
Date: Wed, 15 Jan 2025 13:51:49 -0500
Subject: [PATCH 1/7] add a vignette of fixed design simulations using
 `sim_fixed_n`

---
 vignettes/sim_fixed_design_simple.Rmd | 102 ++++++++++++++++++++++++++
 1 file changed, 102 insertions(+)
 create mode 100644 vignettes/sim_fixed_design_simple.Rmd

diff --git a/vignettes/sim_fixed_design_simple.Rmd b/vignettes/sim_fixed_design_simple.Rmd
new file mode 100644
index 00000000..a5d767af
--- /dev/null
+++ b/vignettes/sim_fixed_design_simple.Rmd
@@ -0,0 +1,102 @@
+---
+title: "Simulate Fixed Designs with Ease via sim_fixed_n"
+author: "Yujie Zhao and Keaven Anderson"
+output: rmarkdown::html_vignette
+bibliography: simtrial.bib
+vignette: >
+  %\VignetteIndexEntry{Simulate Fixed Designs with Ease via sim_fixed_n}
+  %\VignetteEngine{knitr::rmarkdown}
+---
+
+```{r, message=FALSE, warning=FALSE}
+library(gsDesign2)
+library(simtrial)
+library(dplyr)
+library(gt)
+
+set.seed(2025)
+```
+
+The `sim_fixed_n()` function simulates a two-arm trial with a single endpoint, accounting for time-varying enrollment, hazard ratios, and failure and dropout rates. 
+
+There are many advantages of calling `sim_fixed_n()` directly.
+
+- It is simple, which allows for a single function call to perform thousands of simulations.
+- It automatically implements a parallel computation backend, allowing users to achieve shorter running times.
+- It offers up to 5 options for determining the data cutoff through the `timing_type` parameter, allowing for flexibility in cuttings
+
+If people are interested in more complicated simulations, please refer to the vignette [Custom Fixed Design Simulations: A Tutorial on Writing Code from the Ground Up](https://merck.github.io/simtrial/articles/sim_fixed_design_custom.html).
+
+The process for simulating via `sim_fixed_n()` is outlined in Steps 1 to 3 below.
+
+# Step 1: Define design paramaters
+
+To run simulations for a fixed design, several design characteristics are required. The following lines of code create an unstratified 2-arm trial with equal randomization. The total sample size is set to 500, with a target event count of 300 and a total study duration of 36 months. The simulation is repeated 100 times. Enrollment will last for 12 months at a consistent enrollment rate. The median for the control arm is 10 months, with a delayed effect during the first 3 months followed by a hazard ratio of 0.7 thereafter. Additionally, there is a dropout rate of 0.001 over time.
+
+```{r}
+n_sim <- 100
+sample_size <- 500
+target_event <- 300
+total_duration <- 36
+stratum <- data.frame(stratum = "All", p = 1)
+block <- rep(c("experimental", "control"), 2)
+
+enroll_rate <- data.frame(stratum = "All", rate = 12, duration = 500 / 12)
+fail_rate <- data.frame(stratum = "All",
+                        duration = c(3, Inf), fail_rate = log(2) / 10, 
+                        hr = c(1, 0.6), dropout_rate = 0.001)
+```
+
+Besides specifying the design characteristics mentioned above, users also have the option to set these characteristics based on an asymptotic design (`x`). The following design computes the sample size for 85\% power. In this approach, users can obtain the sample size and targeted events from the output of `x`, specifically by using `sample_size <- x$analysis$n` and `target_event <- x$analysis$event`.
+```{r}
+x <- fixed_design_ahr(enroll_rate = enroll_rate, fail_rate = fail_rate, 
+                      alpha = 0.025, power = 0.85, ratio = 1, 
+                      study_duration = total_duration) |> to_integer()
+
+sample_size <- x$analysis$n
+target_event <- x$analysis$event
+```
+
+
+# Step 2: Run `sim_fixed_n()`
+
+Now that we have set up the design characteristics in Step 1, we can proceed to run `sim_fix_n()` for our simulations. This function automatically utilizes a parallel computing backend, which helps reduce the running time.
+
+The `timing_type` argument can take values of `1`, `2`, `3`, `4`, `5`, or any combination of these as a vector. Each number corresponds to a specific cutoff:
+  + `timing_type = 1`: Uses the planned study duration.
+  + `timing_type = 2`: The planned minimum follow-up after enrollment is complete.
+  + `timing_type = 3`: Reflects the planned minimum follow-up period after enrollment completion.
+  + `timing_type = 4`: The maximum of planned study duration and targeted event count cuts (i.e., using `timing_type = 1` and `timing_type = 2` together).
+  + `timing_type = 5`: The maximum of targeted event count and minimum follow-up cuts  (i.e., using `timing_type = 2` and `timing_type = 3` together).
+
+The `rho_gamma` argument is a data frame containing the variables `rho` and `gamma`, both of which should be greater than or equal to zero, to specify one Fleming-Harrington weighted log-rank test per row. For instance, setting `rho = 0 and gamma = 0` gives you the standard unweighted log-rank test, while `rho = 0 and gamma = 0.5` provides the weighted Fleming-Harrington (0, 0.5) log-rank test. If you're interested in tests other than the Fleming-Harrington weighted log-rank test, please refer to the vignette [Custom Fixed Design Simulations: A Tutorial on Writing Code from the Ground Up](https://merck.github.io/simtrial/articles/sim_fixed_design_custom.html).
+
+```{r, message=FALSE}
+sim_res <- sim_fixed_n(
+  n_sim = n_sim,
+  sample_size = sample_size, block = block, stratum = stratum,
+  target_event = target_event, total_duration = total_duration,
+  enroll_rate = enroll_rate, fail_rate = fail_rate,
+  timing_type = 1:5, rho_gamma = data.frame(rho = 0, gamma = 0))
+```
+
+The output of `sim_fixed_n` is a data frame with one row per simulated dataset per cutoff specified in `timing_type`, per test statistic specified in `rho_gamma`. 
+```{r}
+sim_res |> head() |> gt() |> tab_header("Overview Each Simulation Results")
+```
+
+# Step 3: Summarize simulations
+
+With the `r n_sim` simulations provided, users can summarize the simulated power and compare it to the target power of 85\% using some `dplyr` data manipulation.
+```{r}
+sim_res |>
+  group_by(cut) |>
+  summarize(`Simulated Power` = mean(z > qnorm(1 - 0.025))) |>
+  mutate(`Sample size` = sample_size,
+         `Targeted event` = target_event) |>
+  gt() |>
+  tab_header(title = "Summary of 100 simulations by 5 different cuts",
+             subtitle = "Tested by logrank")
+```
+
+## References

From 40818984dd2345ad17b033156dee26f6b9d0e9f2 Mon Sep 17 00:00:00 2001
From: "Zhao, Yujie" <yujie.zhao@merck.com>
Date: Wed, 15 Jan 2025 13:52:11 -0500
Subject: [PATCH 2/7] add a vignette of fixed design simulations by writing
 from the group up

---
 vignettes/sim_fixed_design_custom.Rmd | 349 ++++++++++++++++++++++++++
 1 file changed, 349 insertions(+)
 create mode 100644 vignettes/sim_fixed_design_custom.Rmd

diff --git a/vignettes/sim_fixed_design_custom.Rmd b/vignettes/sim_fixed_design_custom.Rmd
new file mode 100644
index 00000000..cc9cc371
--- /dev/null
+++ b/vignettes/sim_fixed_design_custom.Rmd
@@ -0,0 +1,349 @@
+---
+title: "Custom Fixed Design Simulations: A Tutorial on Writing Code from the Ground Up"
+author: "Yujie Zhao and Keaven Anderson"
+output: rmarkdown::html_vignette
+bibliography: simtrial.bib
+vignette: >
+  %\VignetteIndexEntry{Custom Fixed Design Simulations: A Tutorial on Writing Code from the Ground Up}
+  %\VignetteEngine{knitr::rmarkdown}
+---
+
+```{r, message=FALSE, warning=FALSE}
+library(gsDesign2)
+library(simtrial)
+library(dplyr)
+library(tibble)
+library(gt)
+library(doFuture)
+library(doRNG)
+library(tictoc)
+
+set.seed(2025)
+```
+The vignette [Simulate Fixed Designs with Ease via sim_fixed_n](https://merck.github.io/simtrial/articles/sim_fixed_design_simple.html) presents fixed design simulations using a single function call, `sim_fixed_n()`. It offers a simple and straightforward process for running simulations quickly.
+
+If users are interested in the following aspects, we recommend simulating from scratch rather than directly using `sim_fixed_n()`:
+
+- Tests beyond the logrank test or Fleming-Harrington weighted logrank tests, such as modestly weighted logrank tests, RMST tests, and milestone tests.
+- More complex cutoffs, such as analyzing data after 12 months of follow-up when at least 80% of the patient population is enrolled.
+- Different dropout rates in the control and experimental groups.
+
+The process for simulating from scratch is outlined in Steps 1 to 5 below.
+
+# Step 1: Simulate time-to-event data
+
+The `sim_pw_surv()` function allows for the simulation of a clinical trial with essentially arbitrary patterns of enrollment, failure rates, and censoring. To implement `sim_pw_surv()`, you need to specify 5 design characteristics to simulate time-to-event data:
+
+1. Sample Size (input as `n`).
+1. Stratified or Non-Stratified Designs (input as `stratum`).
+1. Randomization Ratio (input as `block`). The `sim_pw_surv()` function uses fixed block randomization.
+1. Enrollment Rate (input as `enroll_rate`). The `sim_pw_surv()` function supports piecewise enrollment, allowing the enrollment rate to be piecewise constant.
+1. Failure Rate (input as `fail_rate`) or time-to-event rate. The `sim_pw_surv()` function uses a piecewise exponential distribution for the failure rate, which makes it easy to define a distribution with changing failure rates over time. Specifically, in the $j$-th interval, the rate is denoted as $\lambda_j \geq 0$. We require that at least one interval has $\lambda_j > 0$. There are two methods for defining the failure rate:
+  + Specify the failure rate by treatment group, stratum, and time period. An example can be found in Scenario b).
+  + Create a `fail_rate` using `gsDesign2::define_fail_rate`, and then convert it to the required format using `to_sim_pw_surv()`. An example is provided in Scenario a).
+1. Dropout Rate (input as `dropout_Rate`). The `sim_pw_surv()` function accepts piecewise constant dropout rates, which may vary by treatment group. The configuration for dropout should be specified by treatment group, stratum, and time period, and setting up the dropout rate follows the same approach as the failure rate.
+
+
+## Scenario a) The simplest scenario
+
+We begin with the simplest implementation of `sim_pw_surv()`. The following lines of code will generate 500 subjects using equal randomization and an unstratified design.
+
+```{r}
+n_sim <- 100
+n <- 500
+stratum <- data.frame(stratum = "All", p = 1)
+block <- rep(c("experimental", "control"), 2)
+
+enroll_rate <- define_enroll_rate(rate = 12, duration = n / 12)
+
+fail_rate <- define_fail_rate(duration = c(6, Inf), fail_rate = log(2) / 10, 
+                              hr = c(1, 0.7), dropout_rate = 0.0001)
+
+uncut_data_a <- sim_pw_surv(n = n, stratum = stratum, block = block,
+                            enroll_rate = enroll_rate,
+                            fail_rate = to_sim_pw_surv(fail_rate)$fail_rate,
+                            dropout_rate = to_sim_pw_surv(fail_rate)$dropout_rate)
+```
+
+The output of `sim_pw_surv()` is subject-level observations, including stratum, enrollment time for the observation, treatment group the observation is randomized to, failure time, dropout time, calendar time of enrollment plot the minimum of failure time and dropout time (  `cte`), and an failure and dropout indicator (`fail = 1` is a failure, `fail = 0` is a dropout).
+```{r}
+uncut_data_a |> head() |> gt() |> tab_header("An Overview of Simulated TTE data")
+```
+
+## Scenario b) Differential dropout rates
+
+The dropout rate can differ between groups. For instance, in open-label studies, the control group may experience a higher dropout rate. The follow lines of code assumes the control group has a dropout rate of 0.002 for the first 10 months, which then decreases to 0.001 thereafter. In contrast, the experimental group has a constant dropout rate of 0.001 throughout the study. 
+
+```{r}
+differential_dropout_rate <- data.frame(
+  stratum = rep("All", 3), 
+  period = c(1, 2, 1), 
+  treatment = c("control", "control", "experimental"), 
+  duration = c(10, Inf, Inf), 
+  rate = c(0.002, 0.001, 0.001))
+
+uncut_data_b <- sim_pw_surv(n = n, stratum = stratum, block = block,
+                            enroll_rate = enroll_rate,
+                            fail_rate = to_sim_pw_surv(fail_rate)$fail_rate,
+                            dropout_rate = differential_dropout_rate)
+```
+
+## Scenario c) Stratified designs
+
+The following code assumes there are two strata (biomarker-positive and biomarker-negative) with equal prevalence of 0.5 for each. In the control arm, the median survival time is 10 months for biomarker-positive subjects and 8 months for biomarker-negative subjects. For both strata, the hazard ratio is 1 for the first 3 months, after which it decreases to 0.6 for biomarker-positive subjects and 0.8 for biomarker-negative subjects. The dropout rate is contently 0.001 for both strata over time.
+
+```{r}
+stratified_enroll_rate <- data.frame(
+  stratum = c("Biomarker positive", "Biomarker negative"),
+  rate = c(12, 12), 
+  duration = c(1, 1))
+
+stratified_fail_rate <- data.frame(
+  stratum = c(rep("Biomarker positive", 3), rep("Biomarker negative", 3)), 
+  period = c(1, 1, 2, 1, 1, 2), 
+  treatment = rep(c("control", "experimental", "experimental"), 2),
+  duration = c(Inf, 3, Inf, Inf, 3, Inf), 
+  rate = c(# failure rate of biomarker positive subjects: control arm, exp arm period 1, exp arm period 2
+           log(2) / 10, log(2) /10, log(2) / 10 * 0.6,
+           # failure rate of biomarker negative subjects: control arm, exp arm period 1, exp arm period 2
+           log(2) / 8, log(2) /8, log(2) / 8 * 0.8)
+  )
+
+stratified_dropout_rate <- data.frame(
+  stratum = rep(c("Biomarker positive", "Biomarker negative"), each = 2),
+  period = c(1, 1, 1, 1), 
+  treatment = c("control", "experimental", "control", "experimental"),
+  duration = rep(Inf, 4), 
+  rate = rep(0.001, 4)
+  )
+
+uncut_data_c <- sim_pw_surv(n = n,
+                            stratum = data.frame(stratum = c("Biomarker positive", "Biomarker negative"), 
+                                                 p = c(0.5, 0.5)),
+                            block = block,
+                            enroll_rate = stratified_enroll_rate,
+                            fail_rate = stratified_fail_rate,
+                            dropout_rate = stratified_dropout_rate
+                            )
+```
+
+
+For illustration purposes, we will focus on scenario b) for the following discussion.
+
+```{r}
+uncut_data <- uncut_data_b
+```
+
+# Step 2: Cut data
+
+The `get_analysis_date()` derives analysis date for interim/final analysis given multiple conditions, see [the help page of `get_analysis_date()` at the pkgdown website](https://merck.github.io/simtrial/reference/get_analysis_date.html). 
+
+Users can cut for analysis at the 24th month and there are 300 events, whichever arrives later. This is equivalent to `timing_type = 4` in `sim_fixed_n()`.
+```{r}
+cut_date_a <- get_analysis_date(data = uncut_data,
+                                planned_calendar_time = 24,
+                                target_event_overall = 300)
+```
+
+Users can also cut by the maximum of targeted 300 event and minimum follow-up 12 months. This is equivalent to `timing_type = 5` in `sim_fixed_n()`.
+```{r}
+cut_date_b <- get_analysis_date(data = uncut_data,
+                                min_followup = 12,
+                                target_event_overall = 300)
+```
+
+Users can cut data when there are 300 events, with maximum time extension to reach targeted events of 24 months. This is not enabled in `timing_type` of `sim_fixed_n()`.
+```{r}
+cut_date_c <- get_analysis_date(data = uncut_data,
+                                max_extension_for_target_event = 12,
+                                target_event_overall = 300)
+```
+
+Users can cut data after 12 months followup when 80\% of the patients are enrolled in the overall population as below. This is not enabled in `timing_type` of `sim_fixed_n()`.
+```{r}
+cut_date_d <- get_analysis_date(data = uncut_data,
+                                min_n_overall = 100 * 0.8,
+                                min_followup = 12)
+```
+
+More examples are available in the [reference page of `get_analysis_date()`](https://merck.github.io/simtrial/reference/get_analysis_date.html) For illustration purposes, we will focus on scenario d) for the following discussion.
+
+```{r}
+cut_date <- cut_date_d
+cat("The cutoff date is ", round(cut_date, 2))
+
+cut_data <- uncut_data |> cut_data_by_date(cut_date)
+cut_data |> head() |> gt() |> tab_header(paste0("An Overview of TTE data Cut at ", round(cut_date, 2), "Months"))
+```
+
+# Step 3: Run tests
+
+The `simtrial` package provides many options for testing methods, including (weighted) logrank tests, RMST test, milestone test, and MaxComboi test, see the [Section "Compute p-values/test statistics" at the pkgdown reference page] (https://merck.github.io/simtrial/reference/index.html#compute-p-values-test-statistics).
+
+The following code lists all possible tests available in `simtrial`. Users can select one of the tests listed above or combine the testing results to make comparisons across tests. For demonstration purposes, we will aggregate all tests together.
+```{r}
+# Logrank test
+sim_res_lr <- cut_data |> wlr(weight = fh(rho = 0, gamma = 0))
+
+# weighted logrank test by Fleming-Harrington weights
+sim_res_fh <- cut_data |> wlr(weight = fh(rho = 0, gamma = 0.5))
+
+# Modestly weighted logrank test
+sim_res_mb <- cut_data |> wlr(weight = mb(delay = Inf, w_max = 2))
+
+# Weighted logrank test by Xu 2017's early zero weights
+sim_res_xu <- cut_data |> wlr(weight = early_zero(early_period = 3))
+
+# RMST test
+sim_res_rmst <- cut_data |> rmst(tau = 10)
+
+# Milestone test
+sim_res_ms <- cut_data |> milestone(ms_time = 10)
+
+# Maxcombo tests comboing multiple weighted logrank test with Fleming-Harrington weights
+sim_res_mc <- cut_data |> maxcombo(rho = c(0, 0), gamma = c(0, 0.5))
+```
+
+The output of the tests mentioned above are lists including:
+
+- The testing method employed (WLR, RMST, milestone, or MaxCombo), which can be accessed using `sim_res_rmst$method`.
+- The parameters associated with the testing method. For instance, the RMST test parameter is 10, indicating that the RMST is evaluated at month 10. You can find this information using `sim_res_rmst$parameter`.
+- The point estimate and standard error for the testing method used. For example, the point estimate for RMST represents the survival difference between the experimental group and the control group. This estimate can be retrieved with `sim_res_rmst$estimate` and `sim_res_rmst$se`.
+- The Z-score for the testing method, accessible via `sim_res_rmst$z`. Please note that the Z-score is not provided for the MaxCombo test; instead, the p-value is reported (`sim_res_mc$p_value`).
+
+```{r}
+sim_res <- tribble(
+  ~Method, ~Parameter, ~Z, ~Estimate, ~SE, ~`P value`,
+  sim_res_lr$method, sim_res_lr$parameter, sim_res_lr$z, sim_res_lr$estimate, sim_res_lr$se, pnorm(-sim_res_lr$z),
+  sim_res_fh$method, sim_res_fh$parameter, sim_res_fh$z, sim_res_fh$estimate, sim_res_fh$se, pnorm(-sim_res_fh$z),
+  sim_res_mb$method, sim_res_mb$parameter, sim_res_mb$z, sim_res_mb$estimate, sim_res_mb$se, pnorm(-sim_res_mb$z),
+  sim_res_xu$method, sim_res_xu$parameter, sim_res_xu$z, sim_res_xu$estimate, sim_res_xu$se, pnorm(-sim_res_xu$z),
+  sim_res_rmst$method, sim_res_rmst$parameter|> as.character(), sim_res_rmst$z, sim_res_rmst$estimate, sim_res_rmst$se, pnorm(-sim_res_rmst$z),
+  sim_res_ms$method, sim_res_ms$parameter |> as.character(), sim_res_ms$z, sim_res_ms$estimate, sim_res_ms$se, pnorm(-sim_res_ms$z),
+  sim_res_mc$method, sim_res_mc$parameter, NA, NA, NA, sim_res_mc$p_value
+  ) 
+
+sim_res |> gt() |> tab_header("One Simulation Results")
+```
+
+
+# Step 4: Perform the above single simulation repeatedly
+
+We will now merge Steps 1 to 3 into a single function named `one_sim()`, which facilitates a single simulation run. The construction of `one_sim()` involves copying all the lines of code from Steps 1 to 3.
+
+```{r}
+one_sim <- function(sim_id = 1, 
+                    # arguments from Step 1: design characteristic
+                    n, stratum, enroll_rate, fail_rate, dropout_rate, block, 
+                    # arguments from Step 2； cutting method
+                    min_n_overall, min_followup,
+                    # arguments from Step 3； testing method
+                    fh, mb, xu, rmst, ms, mc
+                    ) {
+    # Step 1: simulate a time-to-event data
+    uncut_data <- sim_pw_surv(
+      n = n,
+      stratum = stratum,
+      block = block,
+      enroll_rate = enroll_rate,
+      fail_rate = fail_rate,
+      dropout_rate = dropout_rate) 
+    
+    ## Step 2: Cut data
+    cut_date <- get_analysis_date(min_n_overall = min_n_overall, min_followup = min_followup, data = uncut_data)
+    cut_data <- uncut_data |> cut_data_by_date(cut_date)
+    
+    # Step 3: Run tests
+    sim_res_lr <- cut_data |> wlr(weight = fh(rho = 0, gamma = 0))
+    sim_res_fh <- cut_data |> wlr(weight = fh(rho = fh$rho, gamma = fh$gamma))
+    sim_res_mb <- cut_data |> wlr(weight = mb(delay = mb$delay, w_max = mb$w_max))
+    sim_res_xu <- cut_data |> wlr(weight = early_zero(early_period = xu$early_period))
+    sim_res_rmst <- cut_data |> rmst(tau = rmst$tau)
+    sim_res_ms <- cut_data |> milestone(ms_time = ms$ms_time)
+    sim_res_mc <- cut_data |> maxcombo(rho = mc$rho, gamma = mc$gamma)
+    
+    sim_res <- tribble(
+      ~`Sim ID`, ~Method, ~Parameter, ~Z, ~Estimate, ~SE, ~`P value`,
+      sim_id, sim_res_lr$method, sim_res_lr$parameter, sim_res_lr$z, sim_res_lr$estimate, sim_res_lr$se, pnorm(-sim_res_lr$z),
+      sim_id, sim_res_fh$method, sim_res_fh$parameter, sim_res_fh$z, sim_res_fh$estimate, sim_res_fh$se, pnorm(-sim_res_fh$z),
+      sim_id, sim_res_mb$method, sim_res_mb$parameter, sim_res_mb$z, sim_res_mb$estimate, sim_res_mb$se, pnorm(-sim_res_mb$z),
+      sim_id, sim_res_xu$method, sim_res_xu$parameter, sim_res_xu$z, sim_res_xu$estimate, sim_res_xu$se, pnorm(-sim_res_xu$z),
+      sim_id, sim_res_rmst$method, sim_res_rmst$parameter|> as.character(), sim_res_rmst$z, sim_res_rmst$estimate, sim_res_rmst$se, pnorm(-sim_res_rmst$z),
+      sim_id, sim_res_ms$method, sim_res_ms$parameter |> as.character(), sim_res_ms$z, sim_res_ms$estimate, sim_res_ms$se, pnorm(-sim_res_ms$z),
+      sim_id, sim_res_mc$method, sim_res_mc$parameter, NA, NA, NA, sim_res_mc$p_value
+  ) 
+      
+    return(sim_res)
+}
+```
+
+After that, we will execute `one_sim()` multiple times using parallel computation. The following lines of code takes 24 works to run 100 simulations. 
+```{r}
+tic()
+registerDoFuture()
+registerDoRNG()
+plan("multisession", workers = 24)
+set.seed(2025)
+
+ans <- foreach(
+  sim_id = seq_len(n_sim),
+  .combine = "rbind",
+  .errorhandling = "stop"
+  ) %dorng% {
+    ans_new <- one_sim(
+      sim_id = sim_id, 
+      # arguments from Step 1: design characteristic
+      n = n, 
+      stratum = stratum, 
+      enroll_rate = enroll_rate, 
+      fail_rate = to_sim_pw_surv(fail_rate)$fail_rate, 
+      dropout_rate = differential_dropout_rate, 
+      block = block, 
+      # arguments from Step 2； cutting method
+      min_n_overall = 500 * 0.8,
+      min_followup = 12,
+      # arguments from Step 3； testing method
+      fh = list(rho = 0, gamma = 0.5), 
+      mb = list(delay = Inf, w_max = 2), 
+      xu = list(early_period = 3), 
+      rmst = list(tau = 10), 
+      ms = list(ms_time = 10), 
+      mc = list(rho = c(0, 0), gamma = c(0, 0.5))
+      )
+                              
+    ans_new
+  }
+
+plan("sequential")
+toc()
+```
+
+The output from the parallel computation resembles the output of `sim_fix_n()` described in the vignette [Simulate Fixed Designs with Ease via sim_fixed_n](https://merck.github.io/simtrial/articles/sim_fixed_design_simple.html). Each row in the output corresponds to the simulation results for each testing method per each repeation. 
+
+```{r}
+ans |> head() |> gt() |> tab_header("Overview Each Simulation results")
+```
+
+# Step 5: Summarize simulations
+
+Using the 100 parallel simulations provided above, users can summarize the simulated power and compare it across different testing methods with some data manipulation using `dplyr`. Please note that the power calculation for the MaxCombo test differs from the other tests, as it does not report a Z-score.
+
+```{r, message=FALSE}
+ans_non_mc <- ans |>
+  filter(Method != "MaxCombo") |>
+  group_by(Method, Parameter) %>% 
+  summarise(Power = mean(Z > -qnorm(0.025))) |>
+  ungroup()
+
+ans_mc <- ans |>
+  filter(Method == "MaxCombo") |>
+  summarize(Power = mean(`P value` < 0.025), Method = "MaxCombo", Parameter = "FH(0, 0) + FH(0, 0.5)") 
+
+ans_non_mc |>
+  union(ans_mc) |>
+  gt() |>
+  tab_header("Summary from 100 simulations")
+```
+
+
+## References

From 755b2293ce77ec0177503b6bdb91f7c512b94c74 Mon Sep 17 00:00:00 2001
From: "Zhao, Yujie" <yujie.zhao@merck.com>
Date: Wed, 15 Jan 2025 13:52:27 -0500
Subject: [PATCH 3/7] add a vignette of group sequential design simulations
 using `sim_gs_n`

---
 vignettes/sim_gs_design_simple.Rmd | 136 +++++++++++++++++++++++++++++
 1 file changed, 136 insertions(+)
 create mode 100644 vignettes/sim_gs_design_simple.Rmd

diff --git a/vignettes/sim_gs_design_simple.Rmd b/vignettes/sim_gs_design_simple.Rmd
new file mode 100644
index 00000000..da4882ae
--- /dev/null
+++ b/vignettes/sim_gs_design_simple.Rmd
@@ -0,0 +1,136 @@
+---
+title: "Simulate Group Sequential Designs with Ease via sim_gs_n"
+author: "Yujie Zhao and Keaven Anderson"
+output: rmarkdown::html_vignette
+bibliography: simtrial.bib
+vignette: >
+  %\VignetteIndexEntry{Simulate Group Sequential Designs with Ease via sim_gs_n}
+  %\VignetteEngine{knitr::rmarkdown}
+---
+
+```{r, message=FALSE, warning=FALSE}
+library(gsDesign2)
+library(simtrial)
+library(dplyr)
+library(gt)
+
+set.seed(2025)
+```
+
+The `sim_gs_n()` function simulates group sequential designs with fixed sample size and multiple analyses. There are many advantages of calling `sim_gs_n()` directly.
+
+- It is simple, which allows for a single function call to perform thousands of simulations.
+- It allows for a variety of testing methods, such as (weighted) logrank test, RMST test, and milestone test, via the `test = ...` argument.
+- It automatically implements a parallel computation backend, allowing users to achieve shorter running times.
+- It enables flexible data cutting methods for analyses, specifically: (1) planned calendar time, (2) targeted events, (3) maximum time extension to reach targeted events, (4) planned minimum time after the previous analysis, (5) minimal follow-up time after specified enrollment fraction, along with various combinations, using the `cut = ...` argument.
+
+The process for simulating via `sim_gs_n()` is outlined in Steps 1 to 3 below.
+
+If people are interested in more complicated simulations, please refer to the vignette [Custom Group Sequential Design Simulations: Crafting from Scratch](https://merck.github.io/simtrial/articles/sim_gs_design_custom.html).
+
+# Step 1: Define design paramaters
+
+To run simulations for a group sequential design, several design characteristics are required. The following lines of code create an unstratified 2-arm trial with equal randomization. Enrollment will last for 12 months at a consistent enrollment rate. The median for the control arm is 10 months, with a delayed effect during the first 3 months followed by a hazard ratio of 0.6 thereafter. Additionally, there is a dropout rate of 0.001 over time. The set up of these parameters follows a similar way as in the vignette [Simulate Fixed Designs with Ease via sim_fixed_n](https://merck.github.io/simtrial/articles/sim_fixed_design_simple.html).The total sample size is calculated for 90\% power and the simulation is repeated 100 times. 
+
+```{r}
+n_sim <- 100
+stratum <- data.frame(stratum = "All", p = 1)
+block <- rep(c("experimental", "control"), 2)
+enroll_rate <- data.frame(stratum = "All", rate = 1, duration = 12)
+fail_rate <- data.frame(stratum = "All",
+                        duration = c(3, Inf), fail_rate = log(2) / 10, 
+                        hr = c(1, 0.6), dropout_rate = 0.001)
+
+x <- gs_design_ahr(enroll_rate = enroll_rate, fail_rate = fail_rate,
+                   analysis_time = c(12, 24, 36), alpha = 0.025, beta = 0.1,
+                   upper = gs_spending_bound, lower = gs_b,
+                   upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
+                   lpar = rep(-Inf, 3)) |> to_integer()
+
+sample_size <- x$analysis$n |> max()
+event <- x$analysis$event
+eff_bound <- x$bound$z[x$bound$bound == "upper"]
+```
+
+```{r}
+cat("The total sample size is ", sample_size, ". \n")
+cat("The number of events at IA1, IA2 and FA are ", event, ". \n")
+cat("The efficacy bounds at IA1, IA2 and FA are", eff_bound, ". \n")
+```
+
+In addition to the above parameters, there are couple of more parameters required for a group sequential design. One required parameter is the testing method. The following lines of code generate a modestly weighted logrank test. Users can change these to other tests of interest. For example, Fleming-Harrington(0, 0.5) weighted logrank test by  `test <- wlr` and `weight <- fh(0, 0.5)`. More testing methods are available at [reference page of `simtrial`](https://merck.github.io/simtrial/reference/index.html#compute-p-values-test-statistics).
+```{r}
+test <- wlr
+weight <- mb(delay = Inf, w_max = 2)
+```
+
+The final step in constructing a group sequential design is setting the cutting method for each analysis. The `create_cut()` function includes 5 conditions for the cutoff: (1) planned calendar time, (2) targeted events, (3) maximum time extension to reach targeted events, (4) planned minimum time after the previous analysis, (5) minimal follow-up time after specified enrollment fraction. More details and examples are available at the [help page](https://merck.github.io/simtrial/reference/get_analysis_date.html).
+
+A straightforward method for cutting analyses is based on events. For instance, the following lines of code assume there will be 2 interim analyses and 1 final analysis cut when `r event` events occur. **In this event-driven approach, there is no need to update the efficacy boundary.**
+
+```{r}
+ia1_cut <- create_cut(target_event_overall = event[1])
+ia2_cut <- create_cut(target_event_overall = event[2])
+fa_cut <- create_cut(target_event_overall = event[3])
+
+cut <- list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut)
+```
+
+Users can set more complex cutting. For example, 
+
+- The first interim analysis occurs at 20 months, provided there are at least 100 events, which arrives later. However, if the target number of events is not reached, we will wait a maximum of 24 months.
+- The second interim analysis takes place at 32 months and requires at least 200 events, which arrives later. Additionally, this interim analysis will be scheduled at least 10 months after the first interim analysis.
+- The final analysis is set for 45 months and needs to have 350 events, which arrives later.
+
+**Please keep in mind that if the cut is not event-driven, boundary updates are necessary.** You can find more information about [the boundary updates in the boundary update vignette](https://merck.github.io/simtrial/articles/sim_fixed_design_simple.html). In this vignette, we use the event-driven cut for illustrative purposes.
+
+```{r, eval=FALSE}
+ia1_cut <- create_cut(
+  planned_calendar_time = 20, target_event_overall = 100,
+  max_extension_for_target_event = 24)
+
+ia2_cut <- create_cut(
+  planned_calendar_time = 32,
+  target_event_overall = 200,
+  min_time_after_previous_analysis = 10)
+
+fa_cut <- create_cut(
+  planned_calendar_time = 45,
+  target_event_overall = 350)
+
+cut <- list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut)
+```
+
+# Step 2: Run `sim_gs_n()`
+
+Now that we have set up the design characteristics in Step 1, we can proceed to run `sim_gs_n()` for `r n_sim` simulations. This function automatically utilizes a parallel computing backend, which helps reduce the running time.
+
+```{r, message=FALSE}
+sim_res <- sim_gs_n(
+  n_sim = n_sim,
+  sample_size = sample_size, stratum = stratum, block = block,
+  enroll_rate = enroll_rate, fail_rate = fail_rate,
+  test = test, weight = weight, cut = cut)
+```
+
+The output of `sim_gs_n` is a data frame with one row per simulation per analysis. 
+```{r}
+sim_res |> head() |> gt() |> tab_header("Overview Each Simulation results")
+```
+
+# Step 3: Summarize simulations
+
+With the `r n_sim` simulations provided, users can summarize the simulated power and compare it to the target power of 90\% using some `dplyr` data manipulation.
+```{r, message=FALSE}
+sim_res |>
+  left_join(data.frame(analysis = 1:3, eff_bound = eff_bound)) |>
+  group_by(analysis) |>
+  summarize(`Simulated power` = mean(z >= eff_bound)) |>
+  ungroup() |>
+  mutate(`Asymptotic power` = x$bound$probability[x$bound$bound == "upper"]) |>
+  gt() |>
+  tab_header("Summary of 100 simulations") |> 
+  fmt_number(columns = 2:3, decimals = 4)
+```
+
+## References

From 768555eecd12684cdfe2150c4173c644b287b2a4 Mon Sep 17 00:00:00 2001
From: "Zhao, Yujie" <yujie.zhao@merck.com>
Date: Wed, 15 Jan 2025 13:52:48 -0500
Subject: [PATCH 4/7] add a vignette of group sequential design simulations
 from scratch

---
 vignettes/sim_gs_design_custom.Rmd | 203 +++++++++++++++++++++++++++++
 1 file changed, 203 insertions(+)
 create mode 100644 vignettes/sim_gs_design_custom.Rmd

diff --git a/vignettes/sim_gs_design_custom.Rmd b/vignettes/sim_gs_design_custom.Rmd
new file mode 100644
index 00000000..3244421e
--- /dev/null
+++ b/vignettes/sim_gs_design_custom.Rmd
@@ -0,0 +1,203 @@
+---
+title: "Custom Group Sequential Design Simulations: Crafting from Scratch"
+author: "Yujie Zhao and Keaven Anderson"
+output: rmarkdown::html_vignette
+bibliography: simtrial.bib
+vignette: >
+  %\VignetteIndexEntry{Custom Group Sequential Design Simulations: Crafting from Scratch}
+  %\VignetteEngine{knitr::rmarkdown}
+---
+
+```{r, message=FALSE, warning=FALSE}
+library(gsDesign2)
+library(simtrial)
+library(dplyr)
+library(tibble)
+library(gt)
+library(doFuture)
+library(doRNG)
+library(tictoc)
+
+set.seed(2025)
+```
+
+
+The vignette [Simulate Group Sequential Designs with Ease via sim_gs_n](https://merck.github.io/simtrial/articles/sim_gs_design_simple.html) introduces the simulation of group sequential designs using `sim_gs_n()`. This function offers a simple and straightforward method for conducting rapid simulations with just a single function call.
+
+If users are interested in more complex scenarios such as the ones listed below, we recommend simulating from scratch rather than directly using `sim_gs_n()`:
+
+- Comparing different cutoffs.
+- Evaluating various testing methods.
+- Conducting distinct analyses for each test.
+- Analyzing different dropout rates in the control and experimental groups.
+- Testing by the MaxCombo method.
+
+The parameters setup to scratch a group sequential simulation is very similar to Step 1 of the vignette [Simulate Group Sequential Designs with Ease via sim_gs_n](https://merck.github.io/simtrial/articles/sim_gs_design_simple.html). To shorten the length of this vignette, we will use the same design characteristics and cutting method.
+
+```{r}
+n_sim <- 1e4
+stratum <- data.frame(stratum = "All", p = 1)
+block <- rep(c("experimental", "control"), 2)
+enroll_rate <- data.frame(stratum = "All", rate = 1, duration = 12)
+fail_rate <- data.frame(stratum = "All",
+                        duration = c(3, Inf), fail_rate = log(2) / 10, 
+                        hr = c(1, 0.6), dropout_rate = 0.001)
+
+x <- gs_design_ahr(enroll_rate = enroll_rate, fail_rate = fail_rate,
+                   analysis_time = c(12, 24, 36), alpha = 0.025, beta = 0.1,
+                   upper = gs_spending_bound, lower = gs_b,
+                   upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
+                   lpar = rep(-Inf, 3)) |> to_integer()
+
+sample_size <- x$analysis$n |> max()
+event <- x$analysis$event
+eff_bound <- x$bound$z[x$bound$bound == "upper"]
+
+ia1_cut <- create_cut(target_event_overall = event[1])
+ia2_cut <- create_cut(target_event_overall = event[2])
+fa_cut <- create_cut(target_event_overall = event[3])
+
+cut <- list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut)
+```
+
+```{r}
+cat("The total sample size is ", sample_size, ". \n")
+cat("The number of events at IA1, IA2 and FA are ", event, ". \n")
+cat("The efficacy bounds at IA1, IA2 and FA are", eff_bound, ". \n")
+```
+
+The process of simulating group sequential designs from scratch is very similar to that of fixed designs. The key difference is that group sequential designs require multiple analyses. For each analysis, the data cutting and testing adhere to the same procedures as Steps 1 to 3 described in the [fixed design vignette](https://merck.github.io/simtrial/articles/sim_fixed_design_custom.html). Therefore, we will omit these steps and jump to Step 4 of building a function `one_sim()` for 1 single simulation, which includes data generation, multiple cuttings and testings. 
+Instead of running just a single test, we conduct multiple tests for comparison.
+
+```{r}
+one_sim <- function(sim_id = 1, 
+                    # arguments from Step 1: design characteristic
+                    n, stratum, enroll_rate, fail_rate, dropout_rate, block, 
+                    # arguments from Step 2: cutting method
+                    cut,
+                    # arguments from Step 3: testing method
+                    fh, mb, xu, rmst, ms, mc
+                    ) {
+  
+    # Step 1: simulate time-to-event data
+    uncut_data <- sim_pw_surv(
+      n = n,
+      stratum = stratum,
+      block = block,
+      enroll_rate = enroll_rate,
+      fail_rate = fail_rate,
+      dropout_rate = dropout_rate) 
+    
+    # Step 2: Cut data -- !! this is different from the fixed design due to multiple analyses
+    n_analysis <- length(cut)
+    cut_date <- rep(-1, n_analysis)
+    cut_data <- list()
+    for (i in 1:n_analysis) {
+      cut_date[i] <- cut[[i]](uncut_data)
+      cut_data[[i]] <- uncut_data |> cut_data_by_date(cut_date[i])
+    }
+    
+    # Step 3: Run multiple tests -- !! this is different from the fixed design due to multiple analyses
+    sim_res <- NULL
+    for (i in 1:n_analysis) {
+      sim_res_lr <- cut_data[[i]] |> wlr(weight = fh(rho = 0, gamma = 0))
+      sim_res_fh <- cut_data[[i]] |> wlr(weight = fh(rho = fh$rho, gamma = fh$gamma))
+      sim_res_mb <- cut_data[[i]] |> wlr(weight = mb(delay = mb$delay, w_max = mb$w_max))
+      sim_res_xu <- cut_data[[i]] |> wlr(weight = early_zero(early_period = xu$early_period))
+      sim_res_rmst <- cut_data[[i]] |> rmst(tau = rmst$tau)
+      sim_res_ms <- cut_data[[i]] |> milestone(ms_time = ms$ms_time)
+      sim_res_mc <- cut_data[[i]] |> maxcombo(rho = mc$rho, gamma = mc$gamma)
+      
+      sim_res_new <- tribble(
+        ~`Sim ID`, ~Analysis, ~Method, ~Parameter, ~Z, ~Estimate, ~SE, ~`P value`,
+        sim_id, i, sim_res_lr$method, sim_res_lr$parameter, sim_res_lr$z, sim_res_lr$estimate, sim_res_lr$se, pnorm(-sim_res_lr$z),
+        sim_id, i, sim_res_fh$method, sim_res_fh$parameter, sim_res_fh$z, sim_res_fh$estimate, sim_res_fh$se, pnorm(-sim_res_fh$z),
+        sim_id, i, sim_res_mb$method, sim_res_mb$parameter, sim_res_mb$z, sim_res_mb$estimate, sim_res_mb$se, pnorm(-sim_res_mb$z),
+        sim_id, i, sim_res_xu$method, sim_res_xu$parameter, sim_res_xu$z, sim_res_xu$estimate, sim_res_xu$se, pnorm(-sim_res_xu$z),
+        sim_id, i, sim_res_rmst$method, sim_res_rmst$parameter|> as.character(), sim_res_rmst$z, sim_res_rmst$estimate, sim_res_rmst$se, pnorm(-sim_res_rmst$z),
+        sim_id, i, sim_res_ms$method, sim_res_ms$parameter |> as.character(), sim_res_ms$z, sim_res_ms$estimate, sim_res_ms$se, pnorm(-sim_res_ms$z),
+        sim_id, i, sim_res_mc$method, sim_res_mc$parameter, NA, NA, NA, sim_res_mc$p_value) 
+      
+      sim_res <- rbind(sim_res, sim_res_new)
+    }
+    
+    return(sim_res)
+}
+```
+
+Then we run `one_sim()` `r n_sim` times via parallel computation. 
+
+```{r}
+tic()
+registerDoFuture()
+registerDoRNG()
+plan("multisession", workers = 24)
+set.seed(2025)
+
+ans <- foreach(
+  sim_id = seq_len(n_sim),
+  .combine = "rbind",
+  .errorhandling = "stop"
+  ) %dorng% {
+    ans_new <- one_sim(
+      sim_id = sim_id, 
+      # arguments from Step 1: design characteristic
+      n = sample_size, stratum = stratum, enroll_rate = enroll_rate, 
+      fail_rate = to_sim_pw_surv(fail_rate)$fail_rate, 
+      dropout_rate = to_sim_pw_surv(fail_rate)$dropout_rate, 
+      block = block, 
+      # arguments from Step 2； cutting method
+      cut = cut,
+      # arguments from Step 3； testing method
+      fh = list(rho = 0, gamma = 0.5), 
+      mb = list(delay = Inf, w_max = 2), 
+      xu = list(early_period = 3), 
+      rmst = list(tau = 7), 
+      ms = list(ms_time = 10), 
+      mc = list(rho = c(0, 0), gamma = c(0, 0.5))
+      )
+                              
+    ans_new
+  }
+
+plan("sequential")
+toc()
+```
+
+The output of `sim_gs_n` is a data frame with one row per simulation per analysis per testing method. 
+```{r}
+ans |> head() |> gt() |> tab_header("Overview Each Simulation results")
+```
+
+Finally, using the `r n_sim` parallel simulations provided above, users can summarize the simulated power and compare it across different testing methods with some data manipulation using `dplyr`. 
+
+- For weighted logrank test, RMST test, and milestone test, we compare the Z-score of each analysis with their asymptotic efficacy boundaries obtained from the asymptotic design object `x`, i.e., `x$bound$z[x$bound$bound == "upper"]`.
+- The MaxCombo test is different from the tests mentioned above because it does not provide a Z-score. Instead, we compare its p-values with the planned alpha spending, i.e., `gsDesign::sfLDOF(alpha = 0.025, t = x$analysis$info_frac0)$spend`.
+
+```{r, message=FALSE}
+ans_non_mc <- ans |>
+  filter(Method != "MaxCombo") |>
+  left_join(x$bound |> select(analysis, bound, z) |> rename(eff_bound = z, Analysis = analysis)) |>
+  group_by(Analysis, Method, Parameter) %>% 
+  summarise(`Simulated power` = mean(Z > eff_bound)) |>
+  ungroup()
+
+ans_mc <- ans |>
+  filter(Method == "MaxCombo") |>
+  left_join(data.frame(analysis = 1:3, alpha_spend = gsDesign::sfLDOF(alpha = 0.025, t = x$analysis$info_frac0)$spend) |> rename(Analysis = analysis)) |>
+  group_by(Analysis) |>
+  summarize(`Simulated power` = mean(`P value` < alpha_spend), Method = "MaxCombo", Parameter = "FH(0, 0) + FH(0, 0.5)") |>
+  ungroup()
+
+ans_non_mc |>
+  union(ans_mc) |>
+  left_join(tibble(Analysis = 1:3, 
+                   `Asymptotic power of logrank` = x$bound$probability[x$bound$bound == "upper"])) |>
+  arrange(Analysis, Method) |>
+  gt() |>
+  tab_header(paste0("Summary from ", n_sim, " simulations")) |>
+  fmt_number(columns = 4:5, decimals = 4)
+```
+
+
+## References

From 04a233f53ccb1c309821726d849946a076355361 Mon Sep 17 00:00:00 2001
From: "Zhao, Yujie" <yujie.zhao@merck.com>
Date: Wed, 15 Jan 2025 13:53:09 -0500
Subject: [PATCH 5/7] update the pkgdown.yml

---
 _pkgdown.yml | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/_pkgdown.yml b/_pkgdown.yml
index f714373d..f06e55e2 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -86,11 +86,17 @@ articles:
   contents:
   - workflow
   - routines
-- title: "Simulations with NPH tests"
+- title: "Testing methods"
   contents:
   - modest-wlrt
   - maxcombo
   - rmst
+- title: "Simulate fixed/group sequential designs"
+  contents:
+  - sim_fixed_design_simple
+  - sim_gs_design_simple
+  - sim_fixed_design_custom
+  - sim_gs_design_custom
   - parallel
 - title: "NPH distribution approximations"
   contents:

From 6c13add6cbca35c0fe36b6b1d0b4191da1b856ea Mon Sep 17 00:00:00 2001
From: "Zhao, Yujie" <yujie.zhao@merck.com>
Date: Wed, 15 Jan 2025 14:01:58 -0500
Subject: [PATCH 6/7] add doRNG and tictoc as suggests in DESCRIPTION

---
 DESCRIPTION | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 410115a6..4372f7d4 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -60,7 +60,9 @@ Suggests:
     survMisc,
     survRM2,
     testthat,
-    tidyr
+    tidyr,
+    doRNG,
+    tictoc
 LinkingTo:
     Rcpp
 Roxygen: list(markdown = TRUE)

From eb802f74f90f565d68ace1a59c116fafe91ee8db Mon Sep 17 00:00:00 2001
From: "Zhao, Yujie" <yujie.zhao@merck.com>
Date: Wed, 15 Jan 2025 14:06:37 -0500
Subject: [PATCH 7/7] change 24 workers to 4 workers

---
 vignettes/sim_fixed_design_custom.Rmd | 4 ++--
 vignettes/sim_gs_design_custom.Rmd    | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/vignettes/sim_fixed_design_custom.Rmd b/vignettes/sim_fixed_design_custom.Rmd
index cc9cc371..45048bb6 100644
--- a/vignettes/sim_fixed_design_custom.Rmd
+++ b/vignettes/sim_fixed_design_custom.Rmd
@@ -277,12 +277,12 @@ one_sim <- function(sim_id = 1,
 }
 ```
 
-After that, we will execute `one_sim()` multiple times using parallel computation. The following lines of code takes 24 works to run 100 simulations. 
+After that, we will execute `one_sim()` multiple times using parallel computation. The following lines of code takes 4 works to run 100 simulations. 
 ```{r}
 tic()
 registerDoFuture()
 registerDoRNG()
-plan("multisession", workers = 24)
+plan("multisession", workers = 4)
 set.seed(2025)
 
 ans <- foreach(
diff --git a/vignettes/sim_gs_design_custom.Rmd b/vignettes/sim_gs_design_custom.Rmd
index 3244421e..14425be6 100644
--- a/vignettes/sim_gs_design_custom.Rmd
+++ b/vignettes/sim_gs_design_custom.Rmd
@@ -131,7 +131,7 @@ Then we run `one_sim()` `r n_sim` times via parallel computation.
 tic()
 registerDoFuture()
 registerDoRNG()
-plan("multisession", workers = 24)
+plan("multisession", workers = 4)
 set.seed(2025)
 
 ans <- foreach(