diff --git a/R/milestone.R b/R/milestone.R
index 20e3bc91..3cc3d91e 100644
--- a/R/milestone.R
+++ b/R/milestone.R
@@ -24,7 +24,15 @@
 #'   - `treatment` - Grouping variable.
 #' @param ms_time Milestone analysis time.
 #'
-#' @return A data frame containing the test statistics.
+#' @return A data frame containing:
+#' - `method` - The method, always `"milestone"`.
+#' - `z` - Test statistics.
+#' - `ms_time` - Milestone time point.
+#' - `surv0` - Survival rate of the control arm.
+#' - `surv1` - Survival rate of the experimental arm.
+#' - `surv_diff` - Survival difference between the experimental and control arm.
+#' - `std_err0` - Standard error of the control arm.
+#' - `std_err1` - Standard error of the experimental arm.
 #'
 #' @export
 #'
@@ -54,6 +62,11 @@ milestone <- function(data, ms_time) {
     z <- diff_survival / sqrt(var_survival)
   }
 
-  ans <- data.frame(z = z)
+  ans <- data.frame(
+    method = "milestone", z = z, ms_time = ms_time,
+    surv0 = fit_res$surv[1], surv1 = fit_res$surv[2],
+    surv_diff = diff_survival,
+    std_err0 = fit_res$std.err[1], std_err1 = fit_res$std.err[2]
+  )
   return(ans)
 }
diff --git a/R/sim_gs_n.R b/R/sim_gs_n.R
index 043ab23c..d3132974 100644
--- a/R/sim_gs_n.R
+++ b/R/sim_gs_n.R
@@ -22,50 +22,58 @@
 #' arguments will change as we add additional features.
 #'
 #' @inheritParams sim_fixed_n
-#' @param test a test function such as \code{\link{wlr}},
-#'   \code{\link{maxcombo}}, or \code{\link{rmst}}. The simulated data set is
+#' @param test A test function such as [wlr()],
+#'   [maxcombo()], or [rmst()]. The simulated data set is
 #'   passed as the first positional argument to the test function provided.
-#' @param cutting a list of cutting functions created by
-#'   \code{\link{create_cutting}}, see examples
-#' @param seed random seed
+#' @param cutting A list of cutting functions created by [create_cutting()],
+#'   see examples.
+#' @param seed Random seed.
 #' @param ... Arguments passed to the test function provided by the argument
-#'   \code{test}
+#'   `test`.
+#'
+#' @return A data frame summarizing the simulation ID, analysis date,
+#'   z statistics or p-values.
 #'
-#' @return a data frame summarizing the simulation ID, analysis date, z statistics or p-values
 #' @export
 #'
 #' @examples
 #' library(gsDesign2)
 #'
-#' # parameters for enrollment
-#' enroll_rampup_duration <- 4 # duration for enrollment ramp up
-#' enroll_duration <- 16       # total enrollment duration
-#' enroll_rate <- define_enroll_rate(duration = c(enroll_rampup_duration,
-#'                                                enroll_duration - enroll_rampup_duration),
-#'                                   rate = c(10, 30))
+#' # Parameters for enrollment
+#' enroll_rampup_duration <- 4 # Duration for enrollment ramp up
+#' enroll_duration <- 16 # Total enrollment duration
+#' enroll_rate <- define_enroll_rate(
+#'   duration = c(
+#'     enroll_rampup_duration,
+#'     enroll_duration - enroll_rampup_duration
+#'   ),
+#'   rate = c(10, 30)
+#' )
 #'
-#' # parameters for treatment effect
-#' delay_effect_duration <- 3  # delay treatment effect in months
-#' median_col <- 9             # survival median of the control arm
-#' median_exp <- c(9, 14)      # survival median of the experimental arm
+#' # Parameters for treatment effect
+#' delay_effect_duration <- 3 # Delay treatment effect in months
+#' median_col <- 9 # Survival median of the control arm
+#' median_exp <- c(9, 14) # Survival median of the experimental arm
 #' dropout_rate <- 0.001
-#' fail_rate <- define_fail_rate(duration = c(delay_effect_duration, 100),
-#'                               fail_rate = log(2) /  median_col,
-#'                               hr = median_col / median_exp,
-#'                               dropout_rate = dropout_rate)
+#' fail_rate <- define_fail_rate(
+#'   duration = c(delay_effect_duration, 100),
+#'   fail_rate = log(2) / median_col,
+#'   hr = median_col / median_exp,
+#'   dropout_rate = dropout_rate
+#' )
 #'
-#' # other related parameters
-#' alpha <- 0.025              # type I error
-#' beta <- 0.1                 # type II error
-#' ratio <- 1                  # randomization ratio (exp:col)
+#' # Other related parameters
+#' alpha <- 0.025 # Type I error
+#' beta <- 0.1 # Type II error
+#' ratio <- 1 # Randomization ratio (experimental:control)
 #'
 #' # Define cuttings of 2 IAs and 1 FA
 #' # IA1
 #' # The 1st interim analysis will occur at the later of the following 3 conditions:
-#' # - At least 20 months have passed since the start of the study
-#' # - At least 100 events have occurred
+#' # - At least 20 months have passed since the start of the study.
+#' # - At least 100 events have occurred.
 #' # - At least 20 months have elapsed after enrolling 200/400 subjects, with a
-#' #   minimum of 20 months follow-up
+#' #   minimum of 20 months follow-up.
 #' # However, if events accumulation is slow, we will wait for a maximum of 24 months.
 #' ia1 <- create_cutting(
 #'   planned_calendar_time = 20,
@@ -77,9 +85,9 @@
 #'
 #' # IA2
 #' # The 2nd interim analysis will occur at the later of the following 3 conditions:
-#' # - At least 32 months have passed since the start of the study
-#' # - At least 250 events have occurred
-#' # - At least 10 months after IA1
+#' # - At least 32 months have passed since the start of the study.
+#' # - 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_cutting(
 #'   planned_calendar_time = 32,
@@ -90,8 +98,8 @@
 #'
 #' # FA
 #' # 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
+#' # - At least 45 months have passed since the start of the study.
+#' # - At least 300 events have occurred.
 #' fa <- create_cutting(
 #'   planned_calendar_time = 45,
 #'   target_event_overall = 350
@@ -106,7 +114,8 @@
 #'   test = wlr,
 #'   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
 #'   seed = 2024,
-#'   weight = fh(rho = 0, gamma = 0))
+#'   weight = fh(rho = 0, gamma = 0)
+#' )
 #'
 #' # Test 2: weighted logrank test by FH(0, 0.5)
 #' sim_gs_n(
@@ -117,10 +126,10 @@
 #'   test = wlr,
 #'   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
 #'   seed = 2024,
-#'   weight = fh(rho = 0, gamma = 0.5))
-#'
+#'   weight = fh(rho = 0, gamma = 0.5)
+#' )
 #'
-#' # Test 3: weighted logrank test by MB(6)
+#' # Test 3: weighted logrank test by MB(3)
 #' sim_gs_n(
 #'   n_sim = 3,
 #'   sample_size = 400,
@@ -129,7 +138,8 @@
 #'   test = wlr,
 #'   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
 #'   seed = 2024,
-#'   weight = mb(delay = 3))
+#'   weight = mb(delay = 3)
+#' )
 #'
 #' # Test 4: weighted logrank test by early zero (6)
 #' sim_gs_n(
@@ -140,7 +150,8 @@
 #'   test = wlr,
 #'   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
 #'   seed = 2024,
-#'   weight = early_zero(6))
+#'   weight = early_zero(6)
+#' )
 #'
 #' # Test 5: RMST
 #' sim_gs_n(
@@ -151,9 +162,23 @@
 #'   test = rmst,
 #'   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
 #'   seed = 2024,
-#'   tau = 20)
+#'   tau = 20
+#' )
+#'
+#' # Test 6: Milestone
+#' sim_gs_n(
+#'   n_sim = 3,
+#'   sample_size = 400,
+#'   enroll_rate = enroll_rate,
+#'   fail_rate = fail_rate,
+#'   test = milestone,
+#'   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
+#'   seed = 2024,
+#'   ms_time = 10
+#' )
 #'
-#' # Test 6: maxcombo (FH(0,0) + FH(0, 0.5))
+#' # Test 7: MaxCombo (WLR-FH(0,0) + WLR-FH(0, 0.5))
+#' # for all analyses
 #' sim_gs_n(
 #'   n_sim = 3,
 #'   sample_size = 400,
@@ -163,66 +188,94 @@
 #'   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
 #'   seed = 2024,
 #'   rho = c(0, 0),
-#'   gamma = c(0, 0.5))
+#'   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),
+#'   cutting = 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 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),
+#'   cutting = 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)
+#'     )
+#'   )
+#' )
+#' }
 sim_gs_n <- function(
-    # number of simulations
-  n_sim = 1000,
-  # sample size
-  sample_size = 500,
-  # multinomial probability distribution for stratum enrollment
-  stratum = data.frame(stratum = "All", p = 1),
-  # enrollment rates
-  enroll_rate = data.frame(duration = c(2, 2, 10), rate = c(3, 6, 9)),
-  # failure rates
-  fail_rate = data.frame(
-    stratum = "All",
-    duration = c(3, 100),
-    fail_rate = log(2) / c(9, 18),
-    hr = c(.9, .6),
-    dropout_rate = rep(.001, 2)
-  ),
-  # fixed block randomization specification
-  block = rep(c("experimental", "control"), 2),
-  # default is to to logrank testing
-  # but alternative tests (such as rmst, maxcombo) can be specified
-  test = wlr,
-  # cutting for IA(s) and FA
-  cutting = NULL,
-  # random seed
-  seed = 2024,
-  # arguments passed to `test`
-  ...
-){
-  # input checking
+    n_sim = 1000,
+    sample_size = 500,
+    stratum = data.frame(stratum = "All", p = 1),
+    enroll_rate = data.frame(duration = c(2, 2, 10), rate = c(3, 6, 9)),
+    fail_rate = data.frame(
+      stratum = "All",
+      duration = c(3, 100),
+      fail_rate = log(2) / c(9, 18),
+      hr = c(.9, .6),
+      dropout_rate = rep(.001, 2)
+    ),
+    block = rep(c("experimental", "control"), 2),
+    test = wlr,
+    cutting = NULL,
+    seed = 2024,
+    ...) {
+  # Input checking
   # TODO
 
-  # simulate for n_sim times
+  # Simulate for `n_sim` times
   ans <- NULL
   for (sim_id in seq_len(n_sim)) {
     set.seed(seed + sim_id)
-    # generate data
+    # Generate data
     simu_data <- sim_pw_surv(
       n = sample_size,
       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)
+      dropout_rate = to_sim_pw_surv(fail_rate)$dropout_rate
+    )
 
-    # initialize the cut date of IA(s) and FA
+    # Initialize the cut date of IA(s) and FA
     n_analysis <- length(cutting)
     cut_date <- rep(-100, n_analysis)
     ans_1sim <- NULL
 
     for (i_analysis in seq_len(n_analysis)) {
-
-      # get cut date
+      # Get cut date
       cut_date[i_analysis] <- cutting[[i_analysis]](data = simu_data)
 
-      # cut the data
+      # Cut the data
       simu_data_cut <- simu_data |> cut_data_by_date(cut_date[i_analysis])
 
-      # test
+      # Test
       ans_1sim_new <- test(simu_data_cut, ...)
       ans_1sim_new$analysis <- i_analysis
       ans_1sim_new$cut_date <- cut_date[i_analysis]
diff --git a/man/milestone.Rd b/man/milestone.Rd
index c8168d8a..0f039c30 100644
--- a/man/milestone.Rd
+++ b/man/milestone.Rd
@@ -17,7 +17,17 @@ milestone(data, ms_time)
 \item{ms_time}{Milestone analysis time.}
 }
 \value{
-A data frame containing the test statistics.
+A data frame containing:
+\itemize{
+\item \code{method} - The method, always \code{"milestone"}.
+\item \code{z} - Test statistics.
+\item \code{ms_time} - Milestone time point.
+\item \code{surv0} - Survival rate of the control arm.
+\item \code{surv1} - Survival rate of the experimental arm.
+\item \code{surv_diff} - Survival difference between the experimental and control arm.
+\item \code{std_err0} - Standard error of the control arm.
+\item \code{std_err1} - Standard error of the experimental arm.
+}
 }
 \description{
 Milestone test for two survival curves
diff --git a/man/sim_gs_n.Rd b/man/sim_gs_n.Rd
index 4df396b6..ee612317 100644
--- a/man/sim_gs_n.Rd
+++ b/man/sim_gs_n.Rd
@@ -36,20 +36,21 @@ 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}{a test function such as \code{\link{wlr}},
-\code{\link{maxcombo}}, or \code{\link{rmst}}. The simulated data set is
+\item{test}{A test function such as \code{\link[=wlr]{wlr()}},
+\code{\link[=maxcombo]{maxcombo()}}, or \code{\link[=rmst]{rmst()}}. The simulated data set is
 passed as the first positional argument to the test function provided.}
 
-\item{cutting}{a list of cutting functions created by
-\code{\link{create_cutting}}, see examples}
+\item{cutting}{A list of cutting functions created by \code{\link[=create_cutting]{create_cutting()}},
+see examples.}
 
-\item{seed}{random seed}
+\item{seed}{Random seed.}
 
 \item{...}{Arguments passed to the test function provided by the argument
-\code{test}}
+\code{test}.}
 }
 \value{
-a data frame summarizing the simulation ID, analysis date, z statistics or p-values
+A data frame summarizing the simulation ID, analysis date,
+z statistics or p-values.
 }
 \description{
 WARNING: This experimental function is a work-in-progress. The function
@@ -58,35 +59,41 @@ arguments will change as we add additional features.
 \examples{
 library(gsDesign2)
 
-# parameters for enrollment
-enroll_rampup_duration <- 4 # duration for enrollment ramp up
-enroll_duration <- 16       # total enrollment duration
-enroll_rate <- define_enroll_rate(duration = c(enroll_rampup_duration,
-                                               enroll_duration - enroll_rampup_duration),
-                                  rate = c(10, 30))
-
-# parameters for treatment effect
-delay_effect_duration <- 3  # delay treatment effect in months
-median_col <- 9             # survival median of the control arm
-median_exp <- c(9, 14)      # survival median of the experimental arm
+# Parameters for enrollment
+enroll_rampup_duration <- 4 # Duration for enrollment ramp up
+enroll_duration <- 16 # Total enrollment duration
+enroll_rate <- define_enroll_rate(
+  duration = c(
+    enroll_rampup_duration,
+    enroll_duration - enroll_rampup_duration
+  ),
+  rate = c(10, 30)
+)
+
+# Parameters for treatment effect
+delay_effect_duration <- 3 # Delay treatment effect in months
+median_col <- 9 # Survival median of the control arm
+median_exp <- c(9, 14) # Survival median of the experimental arm
 dropout_rate <- 0.001
-fail_rate <- define_fail_rate(duration = c(delay_effect_duration, 100),
-                              fail_rate = log(2) /  median_col,
-                              hr = median_col / median_exp,
-                              dropout_rate = dropout_rate)
+fail_rate <- define_fail_rate(
+  duration = c(delay_effect_duration, 100),
+  fail_rate = log(2) / median_col,
+  hr = median_col / median_exp,
+  dropout_rate = dropout_rate
+)
 
-# other related parameters
-alpha <- 0.025              # type I error
-beta <- 0.1                 # type II error
-ratio <- 1                  # randomization ratio (exp:col)
+# Other related parameters
+alpha <- 0.025 # type I error
+beta <- 0.1 # type II error
+ratio <- 1 # randomization ratio (exp:col)
 
 # Define cuttings of 2 IAs and 1 FA
 # IA1
 # The 1st interim analysis will occur at the later of the following 3 conditions:
-# - At least 20 months have passed since the start of the study
-# - At least 100 events have occurred
+# - At least 20 months have passed since the start of the study.
+# - At least 100 events have occurred.
 # - At least 20 months have elapsed after enrolling 200/400 subjects, with a
-#   minimum of 20 months follow-up
+#   minimum of 20 months follow-up.
 # However, if events accumulation is slow, we will wait for a maximum of 24 months.
 ia1 <- create_cutting(
   planned_calendar_time = 20,
@@ -98,9 +105,9 @@ ia1 <- create_cutting(
 
 # IA2
 # The 2nd interim analysis will occur at the later of the following 3 conditions:
-# - At least 32 months have passed since the start of the study
-# - At least 250 events have occurred
-# - At least 10 months after IA1
+# - At least 32 months have passed since the start of the study.
+# - 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_cutting(
   planned_calendar_time = 32,
@@ -111,8 +118,8 @@ ia2 <- create_cutting(
 
 # FA
 # 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
+# - At least 45 months have passed since the start of the study.
+# - At least 300 events have occurred.
 fa <- create_cutting(
   planned_calendar_time = 45,
   target_event_overall = 350
@@ -127,7 +134,8 @@ sim_gs_n(
   test = wlr,
   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
   seed = 2024,
-  weight = fh(rho = 0, gamma = 0))
+  weight = fh(rho = 0, gamma = 0)
+)
 
 # Test 2: weighted logrank test by FH(0, 0.5)
 sim_gs_n(
@@ -138,10 +146,10 @@ sim_gs_n(
   test = wlr,
   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
   seed = 2024,
-  weight = fh(rho = 0, gamma = 0.5))
-
+  weight = fh(rho = 0, gamma = 0.5)
+)
 
-# Test 3: weighted logrank test by MB(6)
+# Test 3: weighted logrank test by MB(3)
 sim_gs_n(
   n_sim = 3,
   sample_size = 400,
@@ -150,7 +158,8 @@ sim_gs_n(
   test = wlr,
   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
   seed = 2024,
-  weight = mb(delay = 3))
+  weight = mb(delay = 3)
+)
 
 # Test 4: weighted logrank test by early zero (6)
 sim_gs_n(
@@ -161,7 +170,8 @@ sim_gs_n(
   test = wlr,
   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
   seed = 2024,
-  weight = early_zero(6))
+  weight = early_zero(6)
+)
 
 # Test 5: RMST
 sim_gs_n(
@@ -172,9 +182,23 @@ sim_gs_n(
   test = rmst,
   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
   seed = 2024,
-  tau = 20)
+  tau = 20
+)
+
+# Test 6: Milestone
+sim_gs_n(
+  n_sim = 3,
+  sample_size = 400,
+  enroll_rate = enroll_rate,
+  fail_rate = fail_rate,
+  test = milestone,
+  cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
+  seed = 2024,
+  ms_time = 10
+)
 
-# Test 6: maxcombo (FH(0,0) + FH(0, 0.5))
+# Test 7: MaxCombo (WLR-FH(0,0) + WLR-FH(0, 0.5))
+# for all analyses
 sim_gs_n(
   n_sim = 3,
   sample_size = 400,
@@ -184,5 +208,45 @@ sim_gs_n(
   cutting = list(ia1 = ia1, ia2 = ia2, fa = fa),
   seed = 2024,
   rho = c(0, 0),
-  gamma = c(0, 0.5))
+  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),
+  cutting = 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 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),
+  cutting = 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)
+    )
+  )
+)
+}
 }