diff --git a/DESCRIPTION b/DESCRIPTION
index 7cfea7f5..fe53c949 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: simtrial
 Type: Package
 Title: Clinical Trial Simulation
-Version: 0.4.1.1
+Version: 0.4.1.2
 Authors@R: c(
     person("Keaven", "Anderson", email = "keaven_anderson@merck.com", role = c("aut")),
     person("Yujie", "Zhao", email = "yujie.zhao@merck.com", role =  c("ctb","cre")),
diff --git a/R/maxcombo.R b/R/maxcombo.R
index 1c29e626..0f4e1d12 100644
--- a/R/maxcombo.R
+++ b/R/maxcombo.R
@@ -35,7 +35,7 @@
 #'
 #' @return A list containing the test method (`method`),
 #' parameters of this test method (`parameter`),
-#' point estimation of the treatment effect (`estimation`),
+#' point estimate of the treatment effect (`estimate`),
 #' standardized error of the treatment effect (`se`),
 #' Z-score of each test of the MaxCombo (`z`),
 #' p-values (`p_value`)
@@ -107,10 +107,10 @@ maxcombo <- function(
 
     res$rho[i] <- rg_unique$rho[i]
     res$gamma[i] <- rg_unique$gamma[i]
-    res$estimation[i] <- weighted_o_minus_e_total
+    res$estimate[i] <- weighted_o_minus_e_total
     res$se[i] <- sqrt(weighted_var_total)
     res$var[i] <- weighted_var_total
-    res$z[i] <- res$estimation[i] / res$se[i]
+    res$z[i] <- res$estimate[i] / res$se[i]
   }
 
   # Merge back to full set of pairs ----
@@ -128,7 +128,7 @@ maxcombo <- function(
   temp <- data.frame(rho = rho, gamma = gamma)
   temp$x <- paste0("FH(", temp$rho, ", ", temp$gamma, ")")
   ans$parameter <- paste(temp$x, collapse = " + ")
-  ans$estimation <- NULL
+  ans$estimate <- NULL
   ans$se <- NULL
 
   # Get z statistics for input rho, gamma combinations
diff --git a/R/milestone.R b/R/milestone.R
index 512fdd68..21f145d1 100644
--- a/R/milestone.R
+++ b/R/milestone.R
@@ -32,7 +32,7 @@
 #' @return A list frame containing:
 #'   - `method` - The method, always `"milestone"`.
 #'   - `parameter` - Milestone time point.
-#'   - `estimation` - Survival difference between the experimental and control arm.
+#'   - `estimate` - Survival difference between the experimental and control arm.
 #'   - `se` - Standard error of the control and experimental arm.
 #'   - `z` - Test statistics.
 #'
@@ -53,7 +53,6 @@
 #' cut_data |>
 #'   milestone(10, test_type = "naive")
 milestone <- function(data, ms_time, test_type = c("log-log", "naive")) {
-
   test_type <- match.arg(test_type)
 
   # Fit into KM curves
@@ -83,7 +82,7 @@ milestone <- function(data, ms_time, test_type = c("log-log", "naive")) {
   if (na_ctrl + na_exp == 2) {
     z <- -Inf
   } else {
-    if (test_type == "naive"){
+    if (test_type == "naive") {
       z_numerator <- surv_diff
       z_denominator <- surv_exp * sqrt(sigma2_exp) + surv_ctrl * sqrt(sigma2_ctrl)
     } else if (test_type == "log-log") {
@@ -92,7 +91,7 @@ milestone <- function(data, ms_time, test_type = c("log-log", "naive")) {
     }
   }
 
-  ans$estimation <- z_numerator
+  ans$estimate <- z_numerator
   ans$se <- z_denominator
   ans$z <- z_numerator / z_denominator
 
diff --git a/R/rmst.R b/R/rmst.R
index 9ea813d9..6bc0ec9a 100644
--- a/R/rmst.R
+++ b/R/rmst.R
@@ -104,7 +104,7 @@ rmst <- function(
   ans <- list()
   ans$method <- "RMST"
   ans$parameter <- tau
-  ans$estimation <- res$rmst_diff$rmst_diff
+  ans$estimate <- res$rmst_diff$rmst_diff
   ans$se <- res$rmst_diff$std
   ans$z <- res$rmst_diff$rmst_diff / res$rmst_diff$std
 
diff --git a/R/sim_fixed_n.R b/R/sim_fixed_n.R
index 1d2b34d9..f8ca3a4b 100644
--- a/R/sim_fixed_n.R
+++ b/R/sim_fixed_n.R
@@ -377,7 +377,7 @@ sim_fixed_n <- function(
   return(results)
 }
 
-# Build a function to calculate test related statistics (e.g., z, estimation, se, etc.) and log-hr
+# Build a function to calculate test related statistics (e.g., z, estimate, se, etc.) and log-hr
 doAnalysis <- function(d, rho_gamma, n_stratum) {
   if (nrow(rho_gamma) == 1) {
     res <- d |>
@@ -386,7 +386,7 @@ doAnalysis <- function(d, rho_gamma, n_stratum) {
     ans <- data.frame(
       method = res$method,
       parameter = res$parameter,
-      estimation = res$estimation,
+      estimate = res$estimate,
       se = res$se,
       z = res$z
     )
@@ -397,7 +397,7 @@ doAnalysis <- function(d, rho_gamma, n_stratum) {
     ans <- data.frame(
       method = rep(res$method, nrow(rho_gamma)),
       parameter = rep(res$parameter, nrow(rho_gamma)),
-      estimation = rep("-", nrow(rho_gamma)),
+      estimate = rep("-", nrow(rho_gamma)),
       se = rep("-", nrow(rho_gamma)),
       z = res$z,
       p_value = rep(res$p_value, nrow(rho_gamma))
diff --git a/R/sim_gs_n.R b/R/sim_gs_n.R
index 1aff59cb..c385836d 100644
--- a/R/sim_gs_n.R
+++ b/R/sim_gs_n.R
@@ -22,15 +22,15 @@
 #' arguments will change as we add additional features.
 #'
 #' @inheritParams sim_fixed_n
-#' @param test One or more test functions such as [wlr()], [maxcombo()], or
-#'   [rmst()]. If a single test function is provided, it will be applied at each
-#'   cut. Alternatively a list of functions created by [create_test()]. The list
-#'   form is experimental and currently limited. It only accepts one test per
-#'   cutting (in the future multiple tests may be accepted), and all the tests
-#'   must consistently return the same exact results (again this may be more
-#'   flexible in the future). Importantly, note that the simulated data set is
-#'   always passed as the first positional argument to each test function
-#'   provided.
+#' @param test One or more test functions such as [wlr()], [rmst()], or
+#'   [milestone()] ([maxcombo()] can only be applied by itself). If a single
+#'   test function is provided, it will be applied at each cut. Alternatively a
+#'   list of functions created by [create_test()]. The list form is experimental
+#'   and currently limited. It only accepts one test per cutting (in the future
+#'   multiple tests may be accepted), and all the tests must consistently return
+#'   the same exact results (again this may be more flexible in the future).
+#'   Importantly, note that the simulated data set is always passed as the first
+#'   positional argument to each test function provided.
 #' @param cut A list of cutting functions created by [create_cut()], see
 #'   examples.
 #' @param seed Random seed.
@@ -81,7 +81,7 @@
 #' # - At least 20 months have elapsed after enrolling 200/400 subjects, with a
 #' #   minimum of 20 months follow-up.
 #' # However, if events accumulation is slow, we will wait for a maximum of 24 months.
-#' ia1 <- create_cut(
+#' ia1_cut <- create_cut(
 #'   planned_calendar_time = 20,
 #'   target_event_overall = 100,
 #'   max_extension_for_target_event = 24,
@@ -95,7 +95,7 @@
 #' # - At least 250 events have occurred.
 #' # - At least 10 months after IA1.
 #' # However, if events accumulation is slow, we will wait for a maximum of 34 months.
-#' ia2 <- create_cut(
+#' ia2_cut <- create_cut(
 #'   planned_calendar_time = 32,
 #'   target_event_overall = 200,
 #'   max_extension_for_target_event = 34,
@@ -106,133 +106,125 @@
 #' # The final analysis will occur at the later of the following 2 conditions:
 #' # - At least 45 months have passed since the start of the study.
 #' # - At least 300 events have occurred.
-#' fa <- create_cut(
+#' fa_cut <- create_cut(
 #'   planned_calendar_time = 45,
 #'   target_event_overall = 350
 #' )
 #'
-#' # Test 1: regular logrank test
+#' # Example 1: regular logrank test at all 3 analyses
 #' sim_gs_n(
 #'   n_sim = 3,
 #'   sample_size = 400,
 #'   enroll_rate = enroll_rate,
 #'   fail_rate = fail_rate,
 #'   test = wlr,
-#'   cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
+#'   cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
 #'   seed = 2024,
 #'   weight = fh(rho = 0, gamma = 0)
 #' )
 #'
-#' # Test 2: weighted logrank test by FH(0, 0.5)
+#' # Example 2: weighted logrank test by FH(0, 0.5) at all 3 analyses
 #' sim_gs_n(
 #'   n_sim = 3,
 #'   sample_size = 400,
 #'   enroll_rate = enroll_rate,
 #'   fail_rate = fail_rate,
 #'   test = wlr,
-#'   cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
+#'   cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
 #'   seed = 2024,
 #'   weight = fh(rho = 0, gamma = 0.5)
 #' )
 #'
-#' # Test 3: weighted logrank test by MB(3)
+#' # Example 3: weighted logrank test by MB(3) at all 3 analyses
 #' sim_gs_n(
 #'   n_sim = 3,
 #'   sample_size = 400,
 #'   enroll_rate = enroll_rate,
 #'   fail_rate = fail_rate,
 #'   test = wlr,
-#'   cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
+#'   cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
 #'   seed = 2024,
 #'   weight = mb(delay = 3)
 #' )
 #'
-#' # Test 4: weighted logrank test by early zero (6)
+#' # Example 4: weighted logrank test by early zero (6) at all 3 analyses
 #' sim_gs_n(
 #'   n_sim = 3,
 #'   sample_size = 400,
 #'   enroll_rate = enroll_rate,
 #'   fail_rate = fail_rate,
 #'   test = wlr,
-#'   cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
+#'   cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
 #'   seed = 2024,
 #'   weight = early_zero(6)
 #' )
 #'
-#' # Test 5: RMST
+#' # Example 5: RMST at all 3 analyses
 #' sim_gs_n(
 #'   n_sim = 3,
 #'   sample_size = 400,
 #'   enroll_rate = enroll_rate,
 #'   fail_rate = fail_rate,
 #'   test = rmst,
-#'   cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
+#'   cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
 #'   seed = 2024,
 #'   tau = 20
 #' )
 #'
-#' # Test 6: Milestone
+#' # Example 6: Milestone at all 3 analyses
 #' sim_gs_n(
 #'   n_sim = 3,
 #'   sample_size = 400,
 #'   enroll_rate = enroll_rate,
 #'   fail_rate = fail_rate,
 #'   test = milestone,
-#'   cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
+#'   cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
 #'   seed = 2024,
 #'   ms_time = 10
 #' )
 #'
-#' # Test 7: MaxCombo (WLR-FH(0,0) + WLR-FH(0, 0.5))
-#' # for all analyses
+#' # Example 7: WLR with fh(0, 0.5) test at IA1,
+#' # WLR with mb(6, Inf) at IA2, and milestone test at FA
+#' ia1_test <- create_test(wlr, weight = fh(rho = 0, gamma = 0.5))
+#' ia2_test <- create_test(wlr, weight = mb(delay = 6, w_max = Inf))
+#' fa_test <- create_test(milestone, ms_time = 10)
 #' sim_gs_n(
 #'   n_sim = 3,
 #'   sample_size = 400,
 #'   enroll_rate = enroll_rate,
 #'   fail_rate = fail_rate,
-#'   test = maxcombo,
-#'   cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
-#'   seed = 2024,
-#'   rho = c(0, 0),
-#'   gamma = c(0, 0.5)
+#'   test = list(ia1 = ia1_test, ia2 = ia2_test, fa = fa_test),
+#'   cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
+#'   seed = 2024
 #' )
 #'
-#' # Test 8: MaxCombo (WLR-FH(0,0.5) + milestone(10))
-#' # for all analyses
-#' \dontrun{
-#' sim_gs_n(
-#'   n_sim = 3,
-#'   sample_size = 400,
-#'   enroll_rate = enroll_rate,
-#'   fail_rate = fail_rate,
-#'   test = maxcombo(test1 = wlr, test2 = milestone),
-#'   cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
-#'   seed = 2024,
-#'   test1_par = list(weight = fh(rho = 0, gamma = 0.5)),
-#'   test2_par = list(ms_time = 10)
+#' # WARNING: Multiple tests per cut will be enabled in a future version.
+#' #          Currently does not work.
+#' # Example 8: At IA1, we conduct 3 tests, LR, WLR with fh(0, 0.5), and RMST test.
+#' # At IA2, we conduct 2 tests, LR and WLR with early zero (6).
+#' # At FA, we conduct 2 tests, LR and milestone test.
+#' ia1_test <- list(
+#'   test1 = create_test(wlr, weight = fh(rho = 0, gamma = 0)),
+#'   test2 = create_test(wlr, weight = fh(rho = 0, gamma = 0.5)),
+#'   test3 = create_test(rmst, tau = 20)
+#' )
+#' ia2_test <- list(
+#'   test1 = create_test(wlr, weight = fh(rho = 0, gamma = 0)),
+#'   test2 = create_test(wlr, weight = early_zero(6))
+#' )
+#' fa_test <- list(
+#'   test1 = create_test(wlr, weight = fh(rho = 0, gamma = 0)),
+#'   test3 = create_test(milestone, ms_time = 20)
 #' )
-#' }
-#'
-#' # Test 9: MaxCombo (WLR-FH(0,0) at IAs
-#' # and WLR-FH(0,0) + milestone(10) + WLR-MB(4,2) at FA)
 #' \dontrun{
 #' sim_gs_n(
 #'   n_sim = 3,
 #'   sample_size = 400,
 #'   enroll_rate = enroll_rate,
 #'   fail_rate = fail_rate,
-#'   test = list(ia1 = wlr, ia2 = wlr, fa = maxcombo),
-#'   cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
-#'   seed = 2024,
-#'   test_par = list(
-#'     ia1 = list(weight = fh(rho = 0, gamma = 0)),
-#'     ia2 = list(weight = fh(rho = 0, gamma = 0)),
-#'     ia3 = list(
-#'       test1_par = list(weight = fh(rho = 0, gamma = 0)),
-#'       test2_par = list(ms_time = 10),
-#'       test3_par = list(delay = 4, w_max = 2)
-#'     )
-#'   )
+#'   test = list(ia1 = ia1_test, ia2 = ia2_test, fa = fa_test),
+#'   cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
+#'   seed = 2024
 #' )
 #' }
 sim_gs_n <- function(
@@ -294,11 +286,18 @@ sim_gs_n <- function(
 
       # Test
       ans_1sim_new <- test[[i_analysis]](simu_data_cut, ...)
-      ans_1sim_new$analysis <- i_analysis
-      ans_1sim_new$cut_date <- cut_date[i_analysis]
-      ans_1sim_new$sim_id <- sim_id
-      ans_1sim_new$n <- nrow(simu_data_cut)
-      ans_1sim_new$event <- sum(simu_data_cut$event)
+      ans_1sim_new <- c(sim_id = sim_id, ans_1sim_new)
+      ans_1sim_new <- append(
+        x = ans_1sim_new,
+        values = c(
+          analysis = i_analysis,
+          cut_date = cut_date[i_analysis],
+          n = nrow(simu_data_cut),
+          event = sum(simu_data_cut$event)
+        ),
+        after = 3
+      )
+      ans_1sim_new <- convert_list_to_df_w_list_cols(ans_1sim_new)
 
       # rbind simulation results for all IA(s) and FA in 1 simulation
       ans_1sim <- rbind(ans_1sim, ans_1sim_new)
@@ -422,3 +421,25 @@ multitest <- function(data, ...) {
   }
   return(output)
 }
+
+# Convert a list to a one row data frame using list columns
+convert_list_to_df_w_list_cols <- function(x) {
+  stopifnot(is.list(x), !is.data.frame(x))
+
+  new_list <- vector(mode = "list", length = length(x))
+  names(new_list) <- names(x)
+
+  for (i in seq_along(x)) {
+    if (length(x[[i]]) > 1) {
+      new_list[[i]] <- I(list(x[[i]]))
+    } else {
+      new_list[i] <- x[i]
+    }
+  }
+
+  # Convert the list to a data frame with one row
+  df_w_list_cols <- do.call(data.frame, new_list)
+  stopifnot(nrow(df_w_list_cols) == 1)
+
+  return(df_w_list_cols)
+}
diff --git a/R/wlr.R b/R/wlr.R
index d897a350..ecc96698 100644
--- a/R/wlr.R
+++ b/R/wlr.R
@@ -27,7 +27,7 @@
 #'
 #' @return A list containing the test method (`method`),
 #' parameters of this test method (`parameter`),
-#' point estimation of the treatment effect (`estimation`),
+#' point estimate of the treatment effect (`estimate`),
 #' standardized error of the treatment effect (`se`),
 #' Z-score (`z`), p-values (`p_value`).
 #'
@@ -95,9 +95,9 @@ wlr <- function(data, weight, return_variance = FALSE) {
     x <- x |> fh_weight(rho = weight$rho, gamma = weight$gamma)
 
     ans$parameter <- paste0("FH(rho=", weight$rho, ", gamma=", weight$gamma, ")")
-    ans$estimation <- sum(x$weight * x$o_minus_e)
+    ans$estimate <- sum(x$weight * x$o_minus_e)
     ans$se <- sqrt(sum(x$weight^2 * x$var_o_minus_e))
-    ans$z <- ans$estimation / ans$se
+    ans$z <- ans$estimate / ans$se
   } else if (inherits(weight, "mb")) {
     x <- x |> mb_weight(delay = weight$delay, w_max = weight$w_max)
 
diff --git a/man/maxcombo.Rd b/man/maxcombo.Rd
index 40648b4b..ebc5fc05 100644
--- a/man/maxcombo.Rd
+++ b/man/maxcombo.Rd
@@ -32,7 +32,7 @@ see details; Default: \code{FALSE}.}
 \value{
 A list containing the test method (\code{method}),
 parameters of this test method (\code{parameter}),
-point estimation of the treatment effect (\code{estimation}),
+point estimate of the treatment effect (\code{estimate}),
 standardized error of the treatment effect (\code{se}),
 Z-score of each test of the MaxCombo (\code{z}),
 p-values (\code{p_value})
diff --git a/man/milestone.Rd b/man/milestone.Rd
index 87e8abfb..0552d8d7 100644
--- a/man/milestone.Rd
+++ b/man/milestone.Rd
@@ -30,7 +30,7 @@ A list frame containing:
 \itemize{
 \item \code{method} - The method, always \code{"milestone"}.
 \item \code{parameter} - Milestone time point.
-\item \code{estimation} - Survival difference between the experimental and control arm.
+\item \code{estimate} - Survival difference between the experimental and control arm.
 \item \code{se} - Standard error of the control and experimental arm.
 \item \code{z} - Test statistics.
 }
diff --git a/man/sim_gs_n.Rd b/man/sim_gs_n.Rd
index 103a39ad..9e93de22 100644
--- a/man/sim_gs_n.Rd
+++ b/man/sim_gs_n.Rd
@@ -36,15 +36,15 @@ for experimental vs. control, and dropout rates by stratum and time period.}
 \item{block}{As in \code{\link[=sim_pw_surv]{sim_pw_surv()}}. Vector of treatments to be included
 in each block.}
 
-\item{test}{One or more test functions such as \code{\link[=wlr]{wlr()}}, \code{\link[=maxcombo]{maxcombo()}}, or
-\code{\link[=rmst]{rmst()}}. If a single test function is provided, it will be applied at each
-cut. Alternatively a list of functions created by \code{\link[=create_test]{create_test()}}. The list
-form is experimental and currently limited. It only accepts one test per
-cutting (in the future multiple tests may be accepted), and all the tests
-must consistently return the same exact results (again this may be more
-flexible in the future). Importantly, note that the simulated data set is
-always passed as the first positional argument to each test function
-provided.}
+\item{test}{One or more test functions such as \code{\link[=wlr]{wlr()}}, \code{\link[=rmst]{rmst()}}, or
+\code{\link[=milestone]{milestone()}} (\code{\link[=maxcombo]{maxcombo()}} can only be applied by itself). If a single
+test function is provided, it will be applied at each cut. Alternatively a
+list of functions created by \code{\link[=create_test]{create_test()}}. The list form is experimental
+and currently limited. It only accepts one test per cutting (in the future
+multiple tests may be accepted), and all the tests must consistently return
+the same exact results (again this may be more flexible in the future).
+Importantly, note that the simulated data set is always passed as the first
+positional argument to each test function provided.}
 
 \item{cut}{A list of cutting functions created by \code{\link[=create_cut]{create_cut()}}, see
 examples.}
@@ -102,7 +102,7 @@ ratio <- 1 # Randomization ratio (experimental:control)
 # - At least 20 months have elapsed after enrolling 200/400 subjects, with a
 #   minimum of 20 months follow-up.
 # However, if events accumulation is slow, we will wait for a maximum of 24 months.
-ia1 <- create_cut(
+ia1_cut <- create_cut(
   planned_calendar_time = 20,
   target_event_overall = 100,
   max_extension_for_target_event = 24,
@@ -116,7 +116,7 @@ ia1 <- create_cut(
 # - At least 250 events have occurred.
 # - At least 10 months after IA1.
 # However, if events accumulation is slow, we will wait for a maximum of 34 months.
-ia2 <- create_cut(
+ia2_cut <- create_cut(
   planned_calendar_time = 32,
   target_event_overall = 200,
   max_extension_for_target_event = 34,
@@ -127,133 +127,119 @@ ia2 <- create_cut(
 # The final analysis will occur at the later of the following 2 conditions:
 # - At least 45 months have passed since the start of the study.
 # - At least 300 events have occurred.
-fa <- create_cut(
+fa_cut <- create_cut(
   planned_calendar_time = 45,
   target_event_overall = 350
 )
 
-# Test 1: regular logrank test
+# Example 1: regular logrank test at all 3 analyses
 sim_gs_n(
   n_sim = 3,
   sample_size = 400,
   enroll_rate = enroll_rate,
   fail_rate = fail_rate,
   test = wlr,
-  cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
+  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
   seed = 2024,
   weight = fh(rho = 0, gamma = 0)
 )
 
-# Test 2: weighted logrank test by FH(0, 0.5)
+# Example 2: weighted logrank test by FH(0, 0.5) at all 3 analyses
 sim_gs_n(
   n_sim = 3,
   sample_size = 400,
   enroll_rate = enroll_rate,
   fail_rate = fail_rate,
   test = wlr,
-  cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
+  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
   seed = 2024,
   weight = fh(rho = 0, gamma = 0.5)
 )
 
-# Test 3: weighted logrank test by MB(3)
+# Example 3: weighted logrank test by MB(3) at all 3 analyses
 sim_gs_n(
   n_sim = 3,
   sample_size = 400,
   enroll_rate = enroll_rate,
   fail_rate = fail_rate,
   test = wlr,
-  cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
+  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
   seed = 2024,
   weight = mb(delay = 3)
 )
 
-# Test 4: weighted logrank test by early zero (6)
+# Example 4: weighted logrank test by early zero (6) at all 3 analyses
 sim_gs_n(
   n_sim = 3,
   sample_size = 400,
   enroll_rate = enroll_rate,
   fail_rate = fail_rate,
   test = wlr,
-  cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
+  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
   seed = 2024,
   weight = early_zero(6)
 )
 
-# Test 5: RMST
+# Example 5: RMST at all 3 analyses
 sim_gs_n(
   n_sim = 3,
   sample_size = 400,
   enroll_rate = enroll_rate,
   fail_rate = fail_rate,
   test = rmst,
-  cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
+  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
   seed = 2024,
   tau = 20
 )
 
-# Test 6: Milestone
+# Example 6: Milestone at all 3 analyses
 sim_gs_n(
   n_sim = 3,
   sample_size = 400,
   enroll_rate = enroll_rate,
   fail_rate = fail_rate,
   test = milestone,
-  cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
+  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
   seed = 2024,
   ms_time = 10
 )
 
-# Test 7: MaxCombo (WLR-FH(0,0) + WLR-FH(0, 0.5))
-# for all analyses
+# Example 7: WLR with fh(0, 0.5) test at IA1,
+# WLR with mb(6, Inf) at IA2, and milestone test at FA
+ia1_test <- create_test(wlr, weight = fh(rho = 0, gamma = 0.5))
+ia2_test <- create_test(wlr, weight = mb(delay = 6, w_max = Inf))
+fa_test <- create_test(milestone, ms_time = 10)
 sim_gs_n(
   n_sim = 3,
   sample_size = 400,
   enroll_rate = enroll_rate,
   fail_rate = fail_rate,
-  test = maxcombo,
-  cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
-  seed = 2024,
-  rho = c(0, 0),
-  gamma = c(0, 0.5)
-)
-
-# Test 8: MaxCombo (WLR-FH(0,0.5) + milestone(10))
-# for all analyses
-\dontrun{
-sim_gs_n(
-  n_sim = 3,
-  sample_size = 400,
-  enroll_rate = enroll_rate,
-  fail_rate = fail_rate,
-  test = maxcombo(test1 = wlr, test2 = milestone),
-  cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
-  seed = 2024,
-  test1_par = list(weight = fh(rho = 0, gamma = 0.5)),
-  test2_par = list(ms_time = 10)
+  test = list(ia1 = ia1_test, ia2 = ia2_test, fa = fa_test),
+  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
+  seed = 2024
 )
-}
 
-# Test 9: MaxCombo (WLR-FH(0,0) at IAs
-# and WLR-FH(0,0) + milestone(10) + WLR-MB(4,2) at FA)
+# WARNING: Multiple tests per cut will be enabled in a future version.
+#          Currently does not work.
+# Example 8: At IA1, we conduct 3 tests, LR, WLR with fh(0, 0.5), and RMST test.
+# At IA2, we conduct 2 tests, LR and WLR with early zero (6).
+# At FA, we conduct 2 tests, LR and milestone test.
+ia1_test <- list(test1 = create_test(wlr, weight = fh(rho = 0, gamma = 0)),
+                 test2 = create_test(wlr, weight = fh(rho = 0, gamma = 0.5)),
+                 test3 = create_test(rmst, tau = 20))
+ia2_test <- list(test1 = create_test(wlr, weight = fh(rho = 0, gamma = 0)),
+                 test2 = create_test(wlr, weight = early_zero(6)))
+fa_test <- list(test1 = create_test(wlr, weight = fh(rho = 0, gamma = 0)),
+                 test3 = create_test(milestone, ms_time = 20))
 \dontrun{
 sim_gs_n(
   n_sim = 3,
   sample_size = 400,
   enroll_rate = enroll_rate,
   fail_rate = fail_rate,
-  test = list(ia1 = wlr, ia2 = wlr, fa = maxcombo),
-  cut = list(ia1 = ia1, ia2 = ia2, fa = fa),
-  seed = 2024,
-  test_par = list(
-    ia1 = list(weight = fh(rho = 0, gamma = 0)),
-    ia2 = list(weight = fh(rho = 0, gamma = 0)),
-    ia3 = list(
-      test1_par = list(weight = fh(rho = 0, gamma = 0)),
-      test2_par = list(ms_time = 10),
-      test3_par = list(delay = 4, w_max = 2)
-    )
-  )
+  test = list(ia1 = ia1_test, ia2 = ia2_test, fa = fa_test),
+  cut = list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut),
+  seed = 2024
 )
 }
 \dontshow{\}) # examplesIf}
diff --git a/man/wlr.Rd b/man/wlr.Rd
index 5b7129dc..d8f39b15 100644
--- a/man/wlr.Rd
+++ b/man/wlr.Rd
@@ -19,7 +19,7 @@ see details; Default: \code{FALSE}.}
 \value{
 A list containing the test method (\code{method}),
 parameters of this test method (\code{parameter}),
-point estimation of the treatment effect (\code{estimation}),
+point estimate of the treatment effect (\code{estimate}),
 standardized error of the treatment effect (\code{se}),
 Z-score (\code{z}), p-values (\code{p_value}).
 }
diff --git a/tests/testthat/helper-sim_gs_n.R b/tests/testthat/helper-sim_gs_n.R
index 70f9a8b2..9dddbfcf 100644
--- a/tests/testthat/helper-sim_gs_n.R
+++ b/tests/testthat/helper-sim_gs_n.R
@@ -42,7 +42,7 @@ test_cutting <- function() {
   # - At least 20 months have elapsed after enrolling 200/400 subjects, with a
   #   minimum of 20 months follow-up
   # However, if events accumulation is slow, we will wait for a maximum of 24 months.
-  ia1 <- create_cut(
+  ia1_cut <- create_cut(
     planned_calendar_time = 20,
     target_event_overall = 100,
     max_extension_for_target_event = 24,
@@ -55,7 +55,7 @@ test_cutting <- function() {
   # - At least 250 events have occurred
   # - At least 10 months after IA1
   # However, if events accumulation is slow, we will wait for a maximum of 34 months.
-  ia2 <- create_cut(
+  ia2_cut <- create_cut(
     planned_calendar_time = 32,
     target_event_overall = 200,
     max_extension_for_target_event = 34,
@@ -65,10 +65,10 @@ test_cutting <- function() {
   # The final analysis will occur at the later of the following 2 conditions:
   # - At least 45 months have passed since the start of the study
   # - At least 300 events have occurred
-  fa <- create_cut(
+  fa_cut <- create_cut(
     planned_calendar_time = 45,
     target_event_overall = 350
   )
 
-  return(list(ia1 = ia1, ia2 = ia2, fa = fa))
+  return(list(ia1 = ia1_cut, ia2 = ia2_cut, fa = fa_cut))
 }
diff --git a/tests/testthat/test-unvalidated-rmst.R b/tests/testthat/test-unvalidated-rmst.R
index 9db59684..5aa72afb 100644
--- a/tests/testthat/test-unvalidated-rmst.R
+++ b/tests/testthat/test-unvalidated-rmst.R
@@ -12,7 +12,7 @@ test_that("rmst() snapshot test", {
     rmst_diff = 0.8650492799679741,
     z = 2.2178796367487963
   )
-  expect_equal(observed$estimation, expected$rmst_diff)
+  expect_equal(observed$estimate, expected$rmst_diff)
   expect_equal(observed$z, expected$z)
 })
 
diff --git a/tests/testthat/test-unvalidated-sim_gs_n.R b/tests/testthat/test-unvalidated-sim_gs_n.R
index b56d23c9..fc9abe8f 100644
--- a/tests/testthat/test-unvalidated-sim_gs_n.R
+++ b/tests/testthat/test-unvalidated-sim_gs_n.R
@@ -3,272 +3,355 @@
 
 # See helper-sim_gs_n.R for helper functions
 
-test_that("Test 1: regular logrank test", {
-  # observed <- sim_gs_n(
-  #   n_sim = 3,
-  #   sample_size = 400,
-  #   enroll_rate = test_enroll_rate(),
-  #   fail_rate = test_fail_rate(),
-  #   test = wlr,
-  #   cut = test_cutting(),
-  #   seed = 2024,
-  #   weight = fh(rho = 0, gamma = 0)
-  # )
-  # expected <- data.frame(
-  #   rho = numeric(9),
-  #   gamma = numeric(9),
-  #   z = c(
-  #     -3.7486049782713247, -4.53034007934394, -4.316452743033609,
-  #     -3.4771440155825752, -3.8631501353780324, -3.2777779731288317,
-  #     -3.075862925191481, -3.619345457605645, -4.2225917786532925
-  #   ),
-  #   analysis = rep(1:3, 3),
-  #   cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
-  #   sim_id = rep(1:3, each = 3L),
-  #   n = rep(400L, 9L),
-  #   event = c(229, 295, 355, 241, 290, 350, 226, 282, 350)
-  # )
-  # expect_equal(observed, expected)
-  expect_equal(1 + 1, 2)
+test_that("regular logrank test", {
+  observed <- sim_gs_n(
+    n_sim = 3,
+    sample_size = 400,
+    enroll_rate = test_enroll_rate(),
+    fail_rate = test_fail_rate(),
+    test = wlr,
+    cut = test_cutting(),
+    seed = 2024,
+    weight = fh(rho = 0, gamma = 0)
+  )
+  expected <- data.frame(
+    sim_id = rep(1:3, each = 3L),
+    method = rep("WLR", 9L),
+    parameter = rep("FH(rho=0, gamma=0)", 9L),
+    analysis = rep(1:3, 3),
+    cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
+    n = rep(400L, 9L),
+    event = c(229, 295, 355, 241, 290, 350, 226, 282, 350),
+    estimate = c(
+      -28.194825408790173, -38.32580538077858, -39.49229553865729,
+      -26.84871111584948, -32.548237296118835, -30.06631062297029,
+      -23.063020152157016, -30.16329862679027, -38.75506042018556
+    ),
+    se = c(
+      7.521418120132856, 8.459807588292295, 9.149247748025164, 7.7214837796562,
+      8.425309955739992, 9.17277218574699, 7.498065002594769, 8.333909813280053,
+      9.178026778744327
+    ),
+    z = c(
+      -3.7486049782713247, -4.53034007934394, -4.316452743033609,
+      -3.4771440155825752, -3.8631501353780324, -3.2777779731288317,
+      -3.075862925191481, -3.619345457605645, -4.2225917786532925
+    )
+  )
+  expect_equal(observed, expected)
 })
 
-test_that("Test 2: weighted logrank test by FH(0, 0.5)", {
-  # observed <- sim_gs_n(
-  #   n_sim = 3,
-  #   sample_size = 400,
-  #   enroll_rate = test_enroll_rate(),
-  #   fail_rate = test_fail_rate(),
-  #   test = wlr,
-  #   cut = test_cutting(),
-  #   seed = 2024,
-  #   weight = fh(rho = 0, gamma = 0.5)
-  # )
-  # expected <- data.frame(
-  #   rho = numeric(9),
-  #   gamma = rep(0.5, 9L),
-  #   z = c(
-  #     -4.149161171743935, -4.778107819550277, -4.2607297587160256,
-  #     -3.605092910242299, -3.945081123231263, -2.919179640988388,
-  #     -3.1432278107909206, -3.640458610667732, -4.243289152457
-  #   ),
-  #   analysis = rep(1:3, 3),
-  #   cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
-  #   sim_id = rep(1:3, each = 3L),
-  #   n = rep(400L, 9L),
-  #   event = c(229, 295, 355, 241, 290, 350, 226, 282, 350)
-  # )
-  # expect_equal(observed, expected)
-  expect_equal(1 + 1, 2)
+test_that("weighted logrank test by FH(0, 0.5)", {
+  observed <- sim_gs_n(
+    n_sim = 3,
+    sample_size = 400,
+    enroll_rate = test_enroll_rate(),
+    fail_rate = test_fail_rate(),
+    test = wlr,
+    cut = test_cutting(),
+    seed = 2024,
+    weight = fh(rho = 0, gamma = 0.5)
+  )
+  expected <- data.frame(
+    sim_id = rep(1:3, each = 3L),
+    method = rep("WLR", 9L),
+    parameter = rep("FH(rho=0, gamma=0.5)", 9L),
+    analysis = rep(1:3, 3),
+    cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
+    n = rep(400L, 9L),
+    event = c(229, 295, 355, 241, 290, 350, 226, 282, 350),
+    estimate = c(
+      -16.934217242190208, -24.448179085866208, -25.51076208491462,
+      -15.500239367897708, -19.967690764549445, -17.5390556887186,
+      -12.664624037110103, -18.05051250570102, -25.59217169575864
+    ),
+    se = c(
+      4.081359229309616, 5.116707284384241, 5.987416130471112, 4.299539499761723,
+      5.061414490811455, 6.008213897648367, 4.02917790229253, 4.958307300296486,
+      6.031210878226377
+    ),
+    z = c(
+      -4.149161171743935, -4.778107819550277, -4.2607297587160256,
+      -3.605092910242299, -3.945081123231263, -2.919179640988388,
+      -3.1432278107909206, -3.640458610667732, -4.243289152457
+    )
+  )
+  expect_equal(observed, expected)
 })
 
-test_that("Test 3: weighted logrank test by MB(3)", {
-  # observed <- sim_gs_n(
-  #   n_sim = 3,
-  #   sample_size = 400,
-  #   enroll_rate = test_enroll_rate(),
-  #   fail_rate = test_fail_rate(),
-  #   test = wlr,
-  #   cut = test_cutting(),
-  #   seed = 2024,
-  #   weight = mb(delay = 3)
-  # )
-  # expected <- data.frame(
-  #   z = c(
-  #     -3.797133894694147, -4.581330588107247, -4.3496437937060906,
-  #     -3.5011312494121394, -3.886541892591609, -3.2792862684447983,
-  #     -3.114079263266195, -3.6587146250230145, -4.2632793831797855
-  #   ),
-  #   analysis = rep(1:3, 3),
-  #   cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
-  #   sim_id = rep(1:3, each = 3L),
-  #   n = rep(400L, 9L),
-  #   event = c(229, 295, 355, 241, 290, 350, 226, 282, 350)
-  # )
-  # expect_equal(observed, expected)
-  expect_equal(1 + 1, 2)
+test_that("weighted logrank test by MB(3)", {
+  observed <- sim_gs_n(
+    n_sim = 3,
+    sample_size = 400,
+    enroll_rate = test_enroll_rate(),
+    fail_rate = test_fail_rate(),
+    test = wlr,
+    cut = test_cutting(),
+    seed = 2024,
+    weight = mb(delay = 3)
+  )
+  expected <- data.frame(
+    sim_id = rep(1:3, each = 3L),
+    method = rep("WLR", 9L),
+    parameter = rep("MB(delay = 3, max_weight = Inf)", 9L),
+    analysis = rep(1:3, 3),
+    cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
+    n = rep(400L, 9L),
+    event = c(229, 295, 355, 241, 290, 350, 226, 282, 350),
+    estimate = c(
+      -34.1924345680359, -46.745479781695614, -48.190848712798775,
+      -32.192766832733724, -39.186116293163025, -36.14077883676622,
+      -27.895635278073794, -36.66945854377384, -47.28630966948101
+    ),
+    se = c(
+      9.004800861990685, 10.203472306286518, 11.079263268070516, 9.194961439431635,
+      10.082514836095871, 11.020928299103936, 8.957907914269184, 10.022497598741575,
+      11.091534337637594
+    ),
+    z = c(
+      -3.797133894694147, -4.581330588107247, -4.3496437937060906,
+      -3.5011312494121394, -3.886541892591609, -3.2792862684447983,
+      -3.114079263266195, -3.6587146250230145, -4.2632793831797855
+    )
+  )
+  expect_equal(observed, expected)
 })
 
-test_that("Test 4: weighted logrank test by early zero (6)", {
-  # observed <- sim_gs_n(
-  #   n_sim = 3,
-  #   sample_size = 400,
-  #   enroll_rate = test_enroll_rate(),
-  #   fail_rate = test_fail_rate(),
-  #   test = wlr,
-  #   cut = test_cutting(),
-  #   seed = 2024,
-  #   weight = early_zero(6)
-  # )
-  # expected <- data.frame(
-  #   z = c(
-  #     -4.552617167258777, -5.188572984743822, -4.686073828268738,
-  #     -3.185533497487861, -3.5975030245947046, -2.786930008687834,
-  #     -2.3673440974318556, -3.0630537456426414, -3.7816194091003705
-  #   ),
-  #   analysis = rep(1:3, 3),
-  #   cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
-  #   sim_id = rep(1:3, each = 3L),
-  #   n = rep(400L, 9L),
-  #   event = c(229, 295, 355, 241, 290, 350, 226, 282, 350)
-  # )
-  # expect_equal(observed, expected)
-  expect_equal(1 + 1, 2)
+test_that("weighted logrank test by early zero (6)", {
+  observed <- sim_gs_n(
+    n_sim = 3,
+    sample_size = 400,
+    enroll_rate = test_enroll_rate(),
+    fail_rate = test_fail_rate(),
+    test = wlr,
+    cut = test_cutting(),
+    seed = 2024,
+    weight = early_zero(6)
+  )
+  expected <- data.frame(
+    sim_id = rep(1:3, each = 3L),
+    method = rep("WLR", 9L),
+    parameter = rep("Xu 2017 with first 6 months of 0 weights", 9L),
+    analysis = rep(1:3, 3),
+    cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
+    n = rep(400L, 9L),
+    event = c(229, 295, 355, 241, 290, 350, 226, 282, 350),
+    estimate = c(
+      -21.998993527998245, -32.129973499986654, -33.29646365786535,
+      -17.199406900467533, -22.89893308073689, -20.417006407588342,
+      -11.776058510868394, -18.876336985501645, -27.468098778896934
+    ),
+    se = c(
+      4.8321641639910355, 6.192448982496682, 7.105407400328277, 5.399223368403168,
+      6.36522969520413, 7.325984629661098, 4.974375513742725, 6.162587585135993,
+      7.2635809708390155
+    ),
+    z = c(
+      -4.552617167258777, -5.188572984743822, -4.686073828268738,
+      -3.185533497487861, -3.5975030245947046, -2.786930008687834,
+      -2.3673440974318556, -3.0630537456426414, -3.7816194091003705
+    )
+  )
+  expect_equal(observed, expected)
 })
 
-test_that("Test 5: RMST", {
-  # observed <- sim_gs_n(
-  #   n_sim = 3,
-  #   sample_size = 400,
-  #   enroll_rate = test_enroll_rate(),
-  #   fail_rate = test_fail_rate(),
-  #   test = rmst,
-  #   cut = test_cutting(),
-  #   seed = 2024,
-  #   tau = 20
-  # )
-  # expected <- data.frame(
-  #   rmst_arm1 = c(
-  #     12.466259284156251, 12.444204897288326, 12.425100778728808,
-  #     12.392111715564337, 12.496963791557544, 12.479119007501355, 12.62769367846186,
-  #     12.737915554271744, 12.740241766667666
-  #   ),
-  #   rmst_arm0 = c(
-  #     9.585107633112955, 9.591073977478539, 9.590592780789704, 9.824721964671674,
-  #     10.097271436421035, 10.110783864663125, 10.340195893022198,
-  #     10.289798076615766, 10.261299533752227
-  #   ),
-  #   rmst_diff = c(
-  #     2.8811516510432966, 2.8531309198097876, 2.834507997939104, 2.567389750892662,
-  #     2.3996923551365086, 2.36833514283823, 2.287497785439662, 2.4481174776559786,
-  #     2.478942232915438
-  #   ),
-  #   z = c(
-  #     3.7899815357169184, 3.991862864282945, 3.980100861311682, 3.474868814723485,
-  #     3.2950209410683957, 3.2541151987300845, 2.9805344295194454,
-  #     3.3009521580248022, 3.3504301652133
-  #   ),
-  #   analysis = rep(1:3, 3),
-  #   cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
-  #   sim_id = rep(1:3, each = 3L),
-  #   n = rep(400L, 9L),
-  #   event = c(229, 295, 355, 241, 290, 350, 226, 282, 350)
-  # )
-  # expect_equal(observed, expected)
-  expect_equal(1 + 1, 2)
+test_that("RMST", {
+  observed <- sim_gs_n(
+    n_sim = 3,
+    sample_size = 400,
+    enroll_rate = test_enroll_rate(),
+    fail_rate = test_fail_rate(),
+    test = rmst,
+    cut = test_cutting(),
+    seed = 2024,
+    tau = 20
+  )
+  expected <- data.frame(
+    sim_id = rep(1:3, each = 3L),
+    method = rep("RMST", 9L),
+    parameter = rep(20, 9L),
+    analysis = rep(1:3, 3),
+    cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
+    n = rep(400L, 9L),
+    event = c(229, 295, 355, 241, 290, 350, 226, 282, 350),
+    estimate = c(
+      2.8811516510432966, 2.8531309198097876, 2.834507997939104, 2.567389750892662,
+      2.3996923551365086, 2.36833514283823, 2.287497785439662, 2.4481174776559786,
+      2.478942232915438
+    ),
+    se = c(
+      0.7602020283980866, 0.7147367073498636, 0.7121698913441514,
+      0.7388450867596175, 0.7282783320819864, 0.7277969580678861,
+      0.7674790677752639, 0.7416397937499536, 0.7398877489385366
+    ),
+    z = c(
+      3.7899815357169184, 3.991862864282945, 3.980100861311682, 3.474868814723485,
+      3.2950209410683957, 3.2541151987300845, 2.9805344295194454,
+      3.3009521580248022, 3.3504301652133
+    )
+  )
+  expect_equal(observed, expected)
 })
 
-test_that("Test 6: Milestone", {
-  # observed <- sim_gs_n(
-  #   n_sim = 3,
-  #   sample_size = 400,
-  #   enroll_rate = test_enroll_rate(),
-  #   fail_rate = test_fail_rate(),
-  #   test = milestone,
-  #   cut = test_cutting(),
-  #   seed = 2024,
-  #   ms_time = 10
-  # )
-  # expected <- data.frame(
-  #   method = rep("milestone", 9L),
-  #   z = c(
-  #     9.252619142383594, 12.078380683791904, 12.078380683791904, 5.565741269919053,
-  #     5.457930240636103, 5.457930240636103, 9.051772787302813, 9.054982526543846,
-  #     9.054982526543846
-  #   ),
-  #   ms_time = rep(10, 9L),
-  #   surv_ctrl = c(
-  #     0.40800409626773176, 0.40972689075630214, 0.40972689075630214,
-  #     0.4718268722892688, 0.46670065754089335, 0.46670065754089335,
-  #     0.46149611243704863, 0.46199999999999974, 0.46199999999999974
-  #   ),
-  #   surv_exp = c(
-  #     0.568975019886668, 0.5849999999999997, 0.5849999999999997, 0.5922853919588814,
-  #     0.5840900715499292, 0.5840900715499292, 0.6150543366195163,
-  #     0.6139773404060171, 0.6139773404060171
-  #   ),
-  #   surv_diff = c(
-  #     0.16097092361893622, 0.1752731092436976, 0.1752731092436976,
-  #     0.12045851966961263, 0.11738941400903585, 0.11738941400903585,
-  #     0.15355822418246762, 0.1519773404060174, 0.1519773404060174
-  #   ),
-  #   std_err_ctrl = c(
-  #     0.03693587681297664, 0.034952703615152854, 0.034952703615152854,
-  #     0.03614098127448581, 0.035432630739150366, 0.035432630739150366,
-  #     0.035815727559287504, 0.03540131462139614, 0.03540131462139614
-  #   ),
-  #   std_err_exp = c(
-  #     0.03662189834863626, 0.03484070894801079, 0.03484070894801079,
-  #     0.035312669921649095, 0.034912158581439694, 0.034912158581439694,
-  #     0.03505127094114008, 0.034738243333119145, 0.034738243333119145
-  #   ),
-  #   analysis = rep(1:3, 3),
-  #   cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
-  #   sim_id = rep(1:3, each = 3L),
-  #   n = rep(400L, 9L),
-  #   event = c(229, 295, 355, 241, 290, 350, 226, 282, 350)
-  # )
-  # expect_equal(observed, expected)
-  expect_equal(1 + 1, 2)
+test_that("Milestone", {
+  observed <- sim_gs_n(
+    n_sim = 3,
+    sample_size = 400,
+    enroll_rate = test_enroll_rate(),
+    fail_rate = test_fail_rate(),
+    test = milestone,
+    cut = test_cutting(),
+    seed = 2024,
+    ms_time = 10,
+    test_type = "naive"
+  )
+  expected <- data.frame(
+    sim_id = rep(1:3, each = 3L),
+    method = rep("milestone", 9L),
+    parameter = rep(10, 9L),
+    analysis = rep(1:3, 3),
+    cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
+    n = rep(400L, 9L),
+    event = c(229, 295, 355, 241, 290, 350, 226, 282, 350),
+    estimate = c(
+      0.16097092361893622, 0.1752731092436976, 0.1752731092436976,
+      0.12045851966961263, 0.11738941400903585, 0.11738941400903585,
+      0.15355822418246762, 0.1519773404060174, 0.1519773404060174
+    ),
+    se = c(
+      0.0735577751616129, 0.06979341256316365, 0.06979341256316365,
+      0.0714536511961349, 0.07034478932059006, 0.07034478932059006,
+      0.07086699850042757, 0.07013955795451529, 0.07013955795451529
+    ),
+    z = c(
+      2.188360418259918, 2.5113130710591616, 2.5113130710591616, 1.685827353160205,
+      1.6687719892662431, 1.6687719892662431, 2.1668509663428335,
+      2.1667849760982665, 2.1667849760982665
+    )
+  )
+  expect_equal(observed, expected)
 })
 
-test_that("Test 7: MaxCombo (WLR-FH(0,0) + WLR-FH(0, 0.5))", {
-  # observed <- sim_gs_n(
-  #   n_sim = 3,
-  #   sample_size = 400,
-  #   enroll_rate = test_enroll_rate(),
-  #   fail_rate = test_fail_rate(),
-  #   test = maxcombo,
-  #   cut = test_cutting(),
-  #   seed = 2024,
-  #   rho = c(0, 0),
-  #   gamma = c(0, 0.5)
-  # )
-  # expected <- data.frame(
-  #   p_value = c(
-  #     2.6155386454673746e-05, 1.4330486162172917e-06, 1.247801863046849e-05,
-  #     0.0002358380298724816, 6.130077643518028e-05, 0.0007667834024346343,
-  #     0.001216230102102256, 0.00020471863687732128, 1.7249355113824194e-05
-  #   ),
-  #   analysis = rep(1:3, 3),
-  #   cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
-  #   sim_id = rep(1:3, each = 3L),
-  #   n = rep(400L, 9L),
-  #   event = c(229, 295, 355, 241, 290, 350, 226, 282, 350)
-  # )
-  # expect_equal(observed, expected)
-  expect_equal(1 + 1, 2)
+test_that("WLR with fh(0, 0.5) test at IA1, WLR with mb(6, Inf) at IA2, and milestone test at FA", {
+  ia1_test <- create_test(wlr, weight = fh(rho = 0, gamma = 0.5))
+  ia2_test <- create_test(wlr, weight = mb(delay = 6, w_max = Inf))
+  fa_test <- create_test(milestone, ms_time = 10, test_type = "naive")
+
+  observed <- sim_gs_n(
+    n_sim = 3,
+    sample_size = 400,
+    enroll_rate = test_enroll_rate(),
+    fail_rate = test_fail_rate(),
+    test = list(ia1 = ia1_test, ia2 = ia2_test, fa = fa_test),
+    cut = test_cutting(),
+    seed = 2024
+  )
+  expected <- data.frame(
+    sim_id = rep(1:3, each = 3L),
+    method = rep(c("WLR", "WLR", "milestone"), 3),
+    parameter = rep(c("FH(rho=0, gamma=0.5)", "MB(delay = 6, max_weight = Inf)", "10"), 3),
+    analysis = rep(1:3, 3),
+    cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
+    n = rep(400L, 9L),
+    event = c(229, 295, 355, 241, 290, 350, 226, 282, 350),
+    estimate = c(
+      -16.934217242190208, -55.13395025199291, 0.1752731092436976,
+      -15.500239367897708, -44.437051762182506, 0.11738941400903585,
+      -12.664624037110103, -41.66249375734963, 0.1519773404060174
+    ),
+    se = c(
+      4.081359229309616, 11.636448672437368, 0.06979341256316365, 4.299539499761723,
+      11.272930548434758, 0.07034478932059006, 4.02917790229253, 11.25721499206715,
+      0.07013955795451529
+    ),
+    z = c(
+      -4.149161171743935, -4.738039225196407, 2.5113130710591616,
+      -3.605092910242299, -3.9419254444313574, 1.6687719892662431,
+      -3.1432278107909206, -3.7009592325196583, 2.1667849760982665
+    )
+  )
+  expect_equal(observed, expected)
+})
+
+test_that("MaxCombo (WLR-FH(0,0) + WLR-FH(0, 0.5))", {
+  observed <- sim_gs_n(
+    n_sim = 3,
+    sample_size = 400,
+    enroll_rate = test_enroll_rate(),
+    fail_rate = test_fail_rate(),
+    test = maxcombo,
+    cut = test_cutting(),
+    seed = 2024,
+    rho = c(0, 0),
+    gamma = c(0, 0.5)
+  )
+  expected <- data.frame(
+    sim_id = rep(1:3, each = 3L),
+    method = rep("MaxCombo", 9L),
+    parameter = rep("FH(0, 0) + FH(0, 0.5)", 9L),
+    analysis = rep(1:3, 3),
+    cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
+    n = rep(400L, 9L),
+    event = c(229, 295, 355, 241, 290, 350, 226, 282, 350),
+    z = I(list(
+      c(-3.7486049782713247, -4.149161171743935),
+      c(-4.53034007934394, -4.778107819550277),
+      c(-4.316452743033609, -4.2607297587160256),
+      c(-3.4771440155825752, -3.605092910242299),
+      c(-3.8631501353780324, -3.945081123231263),
+      c(-3.2777779731288317, -2.919179640988388),
+      c(-3.075862925191481, -3.1432278107909206),
+      c(-3.619345457605645, -3.640458610667732),
+      c(-4.2225917786532925, -4.243289152457)
+    )),
+    p_value = c(
+      2.6155386454673746e-05, 1.4330486162172917e-06, 1.247801863046849e-05,
+      0.0002358380298724816, 6.130077643518028e-05, 0.0007667834024346343,
+      0.001216230102102256, 0.00020471863687732128, 1.7249355113824194e-05
+    )
+  )
+  expect_equal(observed, expected)
 })
 
 test_that("sim_gs_n() accepts different tests per cutting", {
-  # wlr_cut1 <- create_test(wlr, weight = fh(rho = 0, gamma = 0))
-  # wlr_cut2 <- create_test(wlr, weight = fh(rho = 0, gamma = 0.5))
-  # wlr_cut3 <- create_test(wlr, weight = fh(rho = 0.5, gamma = 0))
-  #
-  # observed <- sim_gs_n(
-  #   n_sim = 3,
-  #   sample_size = 400,
-  #   enroll_rate = test_enroll_rate(),
-  #   fail_rate = test_fail_rate(),
-  #   test = list(wlr_cut1, wlr_cut2, wlr_cut3),
-  #   cut = test_cutting(),
-  #   seed = 2024
-  # )
-  # expected <- data.frame(
-  #   rho = rep(c(0, 0, 0.5), 3),
-  #   gamma = rep(c(0, 0.5, 0), 3),
-  #   z = c(
-  #     -3.7486049782713247, -4.778107819550277, -4.189693884801371,
-  #     -3.4771440155825752, -3.945081123231263, -3.438138809871842,
-  #     -3.075862925191481, -3.640458610667732, -3.9489173860678495
-  #   ),
-  #   analysis = rep(1:3, 3),
-  #   cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
-  #   sim_id = rep(1:3, each = 3L),
-  #   n = rep(400L, 9L),
-  #   event = c(229, 295, 355, 241, 290, 350, 226, 282, 350)
-  # )
-  # expect_equal(observed, expected)
-  expect_equal(1 + 1, 2)
+  wlr_cut1 <- create_test(wlr, weight = fh(rho = 0, gamma = 0))
+  wlr_cut2 <- create_test(wlr, weight = fh(rho = 0, gamma = 0.5))
+  wlr_cut3 <- create_test(wlr, weight = fh(rho = 0.5, gamma = 0))
+
+  observed <- sim_gs_n(
+    n_sim = 3,
+    sample_size = 400,
+    enroll_rate = test_enroll_rate(),
+    fail_rate = test_fail_rate(),
+    test = list(wlr_cut1, wlr_cut2, wlr_cut3),
+    cut = test_cutting(),
+    seed = 2024
+  )
+  expected <- data.frame(
+    sim_id = rep(1:3, each = 3L),
+    method = rep("WLR", 9L),
+    parameter = rep(c("FH(rho=0, gamma=0)", "FH(rho=0, gamma=0.5)", "FH(rho=0.5, gamma=0)"), 3),
+    analysis = rep(1:3, 3),
+    cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
+    n = rep(400L, 9L),
+    event = c(229, 295, 355, 241, 290, 350, 226, 282, 350),
+    estimate = c(
+      -28.194825408790173, -24.448179085866208, -28.98456223760244,
+      -26.84871111584948, -19.967690764549445, -23.830324019953483,
+      -23.063020152157016, -18.05051250570102, -27.319166131937404
+    ),
+    se = c(
+      7.521418120132856, 5.116707284384241, 6.918062043326721, 7.7214837796562,
+      5.061414490811455, 6.931169838614448, 7.498065002594769, 4.958307300296486,
+      6.9181407107482125
+    ),
+    z = c(
+      -3.7486049782713247, -4.778107819550277, -4.189693884801371,
+      -3.4771440155825752, -3.945081123231263, -3.438138809871842,
+      -3.075862925191481, -3.640458610667732, -3.9489173860678495
+    )
+  )
+  expect_equal(observed, expected)
 })
 
 test_that("sim_gs_n() requires a test for each cutting", {
@@ -290,3 +373,75 @@ test_that("sim_gs_n() requires a test for each cutting", {
     "If you want to run different tests at each cutting"
   )
 })
+
+test_that("sim_gs_n() can combine wlr(), rmst(), and milestone() tests", {
+  test_cut1 <- create_test(wlr, weight = fh(rho = 0, gamma = 0))
+  test_cut2 <- create_test(rmst, tau = 20)
+  test_cut3 <- create_test(milestone, ms_time = 10, test_type = "naive")
+
+  observed <- sim_gs_n(
+    n_sim = 3,
+    sample_size = 400,
+    enroll_rate = test_enroll_rate(),
+    fail_rate = test_fail_rate(),
+    test = list(test_cut1, test_cut2, test_cut3),
+    cut = test_cutting(),
+    seed = 2024
+  )
+  expected <- data.frame(
+    sim_id = rep(1:3, each = 3L),
+    method = rep(c("WLR", "RMST", "milestone"), 3),
+    parameter = rep(c("FH(rho=0, gamma=0)", "20", "10"), 3),
+    analysis = rep(1:3, 3),
+    cut_date = c(24, 32, 45, 24, 32, 46.219327415802894, 24, 32, 50.86585486314699),
+    n = rep(400L, 9L),
+    event = c(229, 295, 355, 241, 290, 350, 226, 282, 350),
+    estimate = c(
+      -28.194825408790173, 2.8531309198097876, 0.1752731092436976,
+      -26.84871111584948, 2.3996923551365086, 0.11738941400903585,
+      -23.063020152157016, 2.4481174776559786, 0.1519773404060174
+    ),
+    se = c(
+      7.521418120132856, 0.7147367073498636, 0.06979341256316365, 7.7214837796562,
+      0.7282783320819864, 0.07034478932059006, 7.498065002594769,
+      0.7416397937499536, 0.07013955795451529
+    ),
+    z = c(
+      -3.7486049782713247, 3.991862864282945, 2.5113130710591616,
+      -3.4771440155825752, 3.2950209410683957, 1.6687719892662431,
+      -3.075862925191481, 3.3009521580248022, 2.1667849760982665
+    )
+  )
+  expect_equal(observed, expected)
+})
+
+test_that("convert_list_to_df_w_list_cols() is robust to diverse input", {
+  x <- list(
+    num_single = 0.5,
+    num_multi = seq(0, 1, by = 0.1),
+    chr_single = "a",
+    chr_multi = letters,
+    int_1 = 1L,
+    int_multi = 1L:10L
+  )
+  observed <- convert_list_to_df_w_list_cols(x)
+  expected <- data.frame(
+    num_single = 0.5,
+    num_multi = I(list(
+      c(
+        0, 0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6000000000000001,
+        0.7000000000000001, 0.8, 0.9, 1
+      )
+    )),
+    chr_single = "a",
+    chr_multi = I(list(
+      c(
+        "a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l", "m", "n", "o",
+        "p", "q", "r", "s", "t", "u", "v", "w", "x", "y", "z"
+      )
+    )),
+    int_1 = 1L,
+    int_multi = I(list(1:10))
+  )
+  expect_equal(observed, expected)
+})