diff --git a/DESCRIPTION b/DESCRIPTION
index 61f07a8c..9b95638e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -58,9 +58,10 @@ Suggests:
     knitr,
     rmarkdown,
     survMisc,
+    survRM2,
     testthat,
     tidyr
 LinkingTo:
     Rcpp
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.1
diff --git a/NAMESPACE b/NAMESPACE
index c2eb91d5..e4fbdc43 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -11,6 +11,7 @@ export(get_cut_date_by_event)
 export(mb_weight)
 export(pvalue_maxcombo)
 export(randomize_by_fixed_block)
+export(rmst)
 export(rpwexp)
 export(rpwexp_enroll)
 export(sim_fixed_n)
@@ -37,4 +38,5 @@ importFrom(mvtnorm,GenzBretz)
 importFrom(mvtnorm,pmvnorm)
 importFrom(survival,Surv)
 importFrom(survival,is.Surv)
+importFrom(survival,survfit)
 useDynLib(simtrial, .registration = TRUE)
diff --git a/R/check_args.R b/R/check_args.R
new file mode 100644
index 00000000..c1e3ef5c
--- /dev/null
+++ b/R/check_args.R
@@ -0,0 +1,83 @@
+#  Copyright (c) 2024 Merck & Co., Inc., Rahway, NJ, USA and its affiliates.
+#  All rights reserved.
+#
+#  This file is part of the simtrial program.
+#
+#  simtrial is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#' Check argument types, length, or dimension
+#'
+#' @param arg An argument to be checked.
+#' @param type A character vector of candidate argument type.
+#' @param length A numeric value of argument length or `NULL`.
+#' @param dim A numeric vector of argument dimension or `NULL`.
+#'
+#' @details If `type`, `length` or `dim` is `NULL`, the corresponding check will not be executed.
+#'
+#' @section Specification:
+#' \if{latex}{
+#'  \itemize{
+#'    \item Check if arg is NULL.
+#'    \item Extract the type, length and dim information from arg.
+#'    \item Compare with target values and report error message if it does not match.
+#'  }
+#'  }
+#' \if{html}{The contents of this section are shown in PDF user manual only.}
+#'
+#' @return Check failure detailed error message.
+#'
+#' @keywords internal
+#'
+#' @examples
+#' \dontrun{
+#' tbl <- as.data.frame(matrix(1:9, nrow = 3))
+#' simtrial:::check_args(arg = tbl, type = c("data.frame"))
+#'
+#' vec <- c("a", "b", "c")
+#' simtrial:::check_args(arg = vec, type = c("character"), length = 3)
+#' }
+check_args <- function(arg, type, length = NULL, dim = NULL) {
+  if (is.null(arg)) {
+    return(NULL)
+  }
+
+  arg <- as.vector(arg)
+
+  check <- list()
+  message <- list()
+
+  if (!is.null(type)) {
+    check[["type"]] <- any(class(arg) %in% type) & (!is.null(class(arg)))
+    message[["type"]] <- paste("The argument type did not match:", paste(type, collapse = "/"))
+  }
+
+  if (!is.null(length)) {
+    check[["length"]] <- all(length(arg) == length) & (!is.null(length(arg)))
+    message[["length"]] <- paste("The argument length is not", length)
+  }
+
+  if (!is.null(dim)) {
+    check[["dim"]] <- all(dim(arg) == dim) & (!is.null(dim(arg)))
+    message[["dim"]] <- paste("The argument dimension is not", paste(dim, collapse = ","))
+  }
+
+  check <- unlist(check)
+  message <- unlist(message)
+
+  if (!all(unlist(check))) {
+    stop(paste(message[!check], collapse = "\n"))
+  } else {
+    return(NULL)
+  }
+}
diff --git a/R/rmst.R b/R/rmst.R
new file mode 100644
index 00000000..f01bde29
--- /dev/null
+++ b/R/rmst.R
@@ -0,0 +1,278 @@
+#  Copyright (c) 2024 Merck & Co., Inc., Rahway, NJ, USA and its affiliates.
+#  All rights reserved.
+#
+#  This file is part of the simtrial program.
+#
+#  simtrial is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#' RMST difference of 2 arms
+#'
+#' @param data A time-to-event dataset with a column `tte` indicating the
+#'   survival time and a column of `event` indicating whether it is
+#'   event or censor.
+#' @param tau RMST analysis time.
+#' @param var_label_tte Column name of the TTE variable.
+#' @param var_label_event Column name of the event variable.
+#' @param var_label_group Column name of the grouping variable.
+#' @param reference A group label indicating the reference group.
+#' @param alpha Type I error.
+#'
+#' @return The z statistics.
+#'
+#' @export
+#'
+#' @examples
+#' data(ex1_delayed_effect)
+#' rmst(
+#'   data = ex1_delayed_effect,
+#'   var_label_tte = "month",
+#'   var_label_event = "evntd",
+#'   var_label_group = "trt",
+#'   tau = 10,
+#'   reference = "0"
+#' )
+rmst <- function(
+    data,
+    tau = 10,
+    var_label_tte = "tte",
+    var_label_event = "event",
+    var_label_group = "treatment",
+    reference = "control",
+    alpha = 0.05) {
+  res <- rmst_two_arm(
+    time_var = data[[var_label_tte]],
+    event_var = data[[var_label_event]],
+    group_var = data[[var_label_group]],
+    trunc_time = tau,
+    reference = reference,
+    alpha = alpha
+  )
+
+  ans <- data.frame(
+    rmst_arm1 = res$rmst_per_arm$rmst[res$rmst_per_arm$group != reference],
+    rmst_arm0 = res$rmst_per_arm$rmst[res$rmst_per_arm$group == reference],
+    rmst_diff = res$rmst_diff$rmst_diff,
+    z = res$rmst_diff$rmst_diff / res$rmst_diff$std
+  )
+
+  return(ans)
+}
+
+#' Calculate RMST difference
+#'
+#' @inheritParams rmst_single_arm
+#' @param group_var A vector of treatment groups.
+#' @param trunc_time A numeric vector of pre-defined cut-off time point(s).
+#' @param reference Group name of reference group for RMST comparison.
+#'   Default is the first group name by alphabetical order.
+#'
+#' @return
+#' A list of 2 data frames of RMST calculations:
+#' - `rmst_per_arm`: the calculation results per group.
+#' - `rmst_diff`: the calculation results of RMST differences.
+#'
+#' @keywords internal
+#'
+#' @examples
+#' data(ex1_delayed_effect)
+#' with(
+#'   ex1_delayed_effect,
+#'   simtrial:::rmst_two_arm(
+#'     time_var = month,
+#'     event_var = evntd,
+#'     group_var = trt,
+#'     trunc_time = 6,
+#'     reference = "0",
+#'     alpha = 0.05
+#'   )
+#' )
+rmst_two_arm <- function(
+    time_var,
+    event_var,
+    group_var,
+    trunc_time,
+    reference = sort(unique(group_var))[1],
+    alpha = 0.05) {
+  # Input type check
+  check_args(time_var, type = c("integer", "numeric"))
+  check_args(event_var, type = c("integer", "numeric"))
+  check_args(trunc_time, type = c("integer", "numeric"))
+  check_args(reference, type = c("character"))
+  check_args(alpha, type = c("integer", "numeric"))
+
+  # Input value check
+  stopifnot(time_var >= 0)
+  stopifnot(event_var %in% c(0, 1))
+  stopifnot(trunc_time >= 0)
+  stopifnot(0 <= alpha & alpha <= 1)
+
+  # Check truncation time
+  if (any(min(tapply(time_var, group_var, max)) < trunc_time)) {
+    stop(paste0(
+      "The truncation time must be shorter than the minimum of the largest observed time in each group: ",
+      sprintf("%.3f", min(tapply(time_var, group_var, max)))
+    ))
+  }
+
+  g_label <- sort(unique(as.character(group_var)))
+
+  # Calculate RMST for each group by rmst_single()
+  one_rmst <- function(x) {
+    indx <- group_var == x
+    rmst_single_arm(
+      time_var = time_var[indx],
+      event_var = event_var[indx],
+      tau = trunc_time,
+      group_label = x,
+      alpha = alpha
+    )
+  }
+
+  op_single <- do.call(rbind, lapply(g_label, one_rmst))
+
+  # Calculate RMST difference and corresponding confidence intervals between each group with reference group
+  diff_rmst <- function(x) {
+    df_rf <- op_single[op_single$group == reference, ]
+    df2 <- op_single[op_single$group == x, ]
+    cutoff_time <- trunc_time
+    group <- paste(unique(df2$group), "-", unique(df_rf$group))
+    rmst_diff <- df2$rmst - df_rf$rmst
+    variance <- df2$variance + df_rf$variance
+    std <- sqrt(variance)
+    lcl <- rmst_diff - stats::qnorm(1 - alpha / 2) * std
+    ucl <- rmst_diff + stats::qnorm(1 - alpha / 2) * std
+    data.frame(cutoff_time, group, rmst_diff, variance, std, lcl, ucl, stringsAsFactors = FALSE)
+  }
+
+  op_diff <- do.call(rbind, lapply(setdiff(g_label, reference), diff_rmst))
+
+  ans <- list(rmst_per_arm = op_single, rmst_diff = op_diff)
+
+  return(ans)
+}
+
+#' Calculate RMST for a single cut-off time point
+#'
+#' @param time_var A numeric vector of follow up time.
+#' @param event_var A numeric or integer vector of the status indicator;
+#'   0=alive 1=event.
+#' @param tau A value of pre-defined cut-off time point.
+#' @param group_label A character of customized treatment group name.
+#' @param alpha A numeric value of the significant level for RMST
+#'   confidence interval. Default is 0.05.
+#'
+#' @return
+#' A data frame of
+#' - Cutoff time: same as \code{tau};
+#' - Group label: same as \code{group_label};
+#' - Estimated RMST;
+#' - Variance, standard error, and CIs of the estimated RMST;
+#' - Number of events.
+#'
+#' @importFrom survival survfit Surv
+#'
+#' @keywords internal
+#'
+#' @examples
+#' data(ex1_delayed_effect)
+#' data_single_arm <- ex1_delayed_effect[ex1_delayed_effect$trt == 1, ]
+#' simtrial:::rmst_single_arm(
+#'   time_var = data_single_arm$month,
+#'   event_var = data_single_arm$evntd,
+#'   tau = 10,
+#'   group_label = "Treatment 1",
+#'   alpha = 0.05
+#' )
+rmst_single_arm <- function(
+    time_var,
+    event_var,
+    tau,
+    group_label = "Single Group",
+    alpha = 0.05) {
+  # Input type check
+  check_args(time_var, type = c("integer", "numeric"))
+  check_args(event_var, type = c("integer", "numeric"))
+  check_args(tau, type = c("integer", "numeric"), length = 1)
+  check_args(group_label, type = c("character", "factor"), length = 1)
+  check_args(alpha, type = c("integer", "numeric"))
+
+  # Input value check
+  stopifnot(time_var >= 0)
+  stopifnot(event_var %in% c(0, 1))
+  stopifnot(tau >= 0)
+  stopifnot(0 <= alpha & alpha <= 1)
+
+  # Fit a single Kaplan-Meier curve
+  fit <- survival::survfit(survival::Surv(time_var, event_var) ~ 1)
+
+  # Extract survival probability, number of event, number at risk,
+  # and number of censored along with observed time from the fitted model
+  # as a new data frame.
+  df <- data.frame(
+    time = fit$time,
+    n_risk = fit$n.risk,
+    n_event = fit$n.event,
+    n_censor = fit$n.censor,
+    surv = fit$surv,
+    stringsAsFactors = FALSE
+  )
+
+  # Filter df by survival time less or equal to the pre-specified cut-off time point tau
+  df_fit1 <- df[df$time <= tau, ]
+
+  # Add initial value of (time, survival) = (0,1) for calculating time difference
+  df_fit2 <- rbind(df_fit1, c(0, NA, NA, NA, 1))
+
+  # Add cut-off time if no records observed on the pre-specified time point
+  if (max(df_fit1$time) != tau) df_fit2 <- rbind(df_fit2, c(tau, NA, NA, NA, NA))
+
+  # Sort the data frame by time
+  df_fit2 <- df_fit2[order(df_fit2$time), ]
+  n_event <- df_fit2$n_event
+  n_risk <- df_fit2$n_risk
+
+  # Calculate the time difference and set the last value as NA
+  time_diff <- c(diff((sort(df_fit2$time))), NA)
+
+  # Calculate the area under the curve per time interval
+  area <- time_diff * df_fit2$surv
+
+  # Calculate the inverse cumulated area under the curve per observed time point A_i
+  big_a <- rev(c(0, cumsum(rev(area)[-1])))
+
+  # Calculation of dev refers to di / Yi * (Yi - di)
+  dev <- (n_event / (n_risk * (n_risk - n_event))) * (big_a^2)
+
+  # Based on the calculation, create a data frame with below items:
+  # cutoff_time is the input of pre-defined cut-off time point
+  cutoff_time <- tau
+  # group is the input group name
+  group <- group_label
+  # rmst is the estimated RMST
+  rmst <- sum(area, na.rm = TRUE)
+  # std is the standard error of the estimated RMST
+  variance <- sum(dev, na.rm = TRUE) * sum(n_event, na.rm = TRUE) / (sum(n_event, na.rm = TRUE) - 1)
+  std <- sqrt(variance)
+  # lcl and ucl are lower/upper control limit of CIs for RMST
+  lcl <- rmst - stats::qnorm(1 - alpha / 2) * std
+  ucl <- rmst + stats::qnorm(1 - alpha / 2) * std
+  event <- sum(n_event, na.rm = TRUE)
+
+  ans <- data.frame(
+    cutoff_time, group, rmst, variance, std, lcl, ucl, event,
+    stringsAsFactors = FALSE
+  )
+
+  return(ans)
+}
diff --git a/_pkgdown.yml b/_pkgdown.yml
index c7e828f1..ebc94ccf 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -47,6 +47,7 @@ reference:
   contents:
     - counting_process
     - pvalue_maxcombo
+    - rmst
 
 - title: "Randomization algorithms"
   contents:
@@ -78,6 +79,7 @@ articles:
   - arbitrary-hazard
   - modest-wlrt
   - pvalue-maxcombo
+  - rmst
 - title: "Best practices"
   contents:
   - parallel
diff --git a/man/check_args.Rd b/man/check_args.Rd
new file mode 100644
index 00000000..3d466f3d
--- /dev/null
+++ b/man/check_args.Rd
@@ -0,0 +1,48 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/check_args.R
+\name{check_args}
+\alias{check_args}
+\title{Check argument types, length, or dimension}
+\usage{
+check_args(arg, type, length = NULL, dim = NULL)
+}
+\arguments{
+\item{arg}{An argument to be checked.}
+
+\item{type}{A character vector of candidate argument type.}
+
+\item{length}{A numeric value of argument length or \code{NULL}.}
+
+\item{dim}{A numeric vector of argument dimension or \code{NULL}.}
+}
+\value{
+Check failure detailed error message.
+}
+\description{
+Check argument types, length, or dimension
+}
+\details{
+If \code{type}, \code{length} or \code{dim} is \code{NULL}, the corresponding check will not be executed.
+}
+\section{Specification}{
+
+\if{latex}{
+ \itemize{
+   \item Check if arg is NULL.
+   \item Extract the type, length and dim information from arg.
+   \item Compare with target values and report error message if it does not match.
+ }
+ }
+\if{html}{The contents of this section are shown in PDF user manual only.}
+}
+
+\examples{
+\dontrun{
+tbl <- as.data.frame(matrix(1:9, nrow = 3))
+simtrial:::check_args(arg = tbl, type = c("data.frame"))
+
+vec <- c("a", "b", "c")
+simtrial:::check_args(arg = vec, type = c("character"), length = 3)
+}
+}
+\keyword{internal}
diff --git a/man/rmst.Rd b/man/rmst.Rd
new file mode 100644
index 00000000..b6822fdd
--- /dev/null
+++ b/man/rmst.Rd
@@ -0,0 +1,50 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rmst.R
+\name{rmst}
+\alias{rmst}
+\title{RMST difference of 2 arms}
+\usage{
+rmst(
+  data,
+  tau = 10,
+  var_label_tte = "tte",
+  var_label_event = "event",
+  var_label_group = "treatment",
+  reference = "control",
+  alpha = 0.05
+)
+}
+\arguments{
+\item{data}{A time-to-event dataset with a column \code{tte} indicating the
+survival time and a column of \code{event} indicating whether it is
+event or censor.}
+
+\item{tau}{RMST analysis time.}
+
+\item{var_label_tte}{Column name of the TTE variable.}
+
+\item{var_label_event}{Column name of the event variable.}
+
+\item{var_label_group}{Column name of the grouping variable.}
+
+\item{reference}{A group label indicating the reference group.}
+
+\item{alpha}{Type I error.}
+}
+\value{
+The z statistics.
+}
+\description{
+RMST difference of 2 arms
+}
+\examples{
+data(ex1_delayed_effect)
+rmst(
+  data = ex1_delayed_effect,
+  var_label_tte = "month",
+  var_label_event = "evntd",
+  var_label_group = "trt",
+  tau = 10,
+  reference = "0"
+)
+}
diff --git a/man/rmst_single_arm.Rd b/man/rmst_single_arm.Rd
new file mode 100644
index 00000000..4c9e841a
--- /dev/null
+++ b/man/rmst_single_arm.Rd
@@ -0,0 +1,52 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rmst.R
+\name{rmst_single_arm}
+\alias{rmst_single_arm}
+\title{Calculate RMST for a single cut-off time point}
+\usage{
+rmst_single_arm(
+  time_var,
+  event_var,
+  tau,
+  group_label = "Single Group",
+  alpha = 0.05
+)
+}
+\arguments{
+\item{time_var}{A numeric vector of follow up time.}
+
+\item{event_var}{A numeric or integer vector of the status indicator;
+0=alive 1=event.}
+
+\item{tau}{A value of pre-defined cut-off time point.}
+
+\item{group_label}{A character of customized treatment group name.}
+
+\item{alpha}{A numeric value of the significant level for RMST
+confidence interval. Default is 0.05.}
+}
+\value{
+A data frame of
+\itemize{
+\item Cutoff time: same as \code{tau};
+\item Group label: same as \code{group_label};
+\item Estimated RMST;
+\item Variance, standard error, and CIs of the estimated RMST;
+\item Number of events.
+}
+}
+\description{
+Calculate RMST for a single cut-off time point
+}
+\examples{
+data(ex1_delayed_effect)
+data_single_arm <- ex1_delayed_effect[ex1_delayed_effect$trt == 1, ]
+simtrial:::rmst_single_arm(
+  time_var = data_single_arm$month,
+  event_var = data_single_arm$evntd,
+  tau = 10,
+  group_label = "Treatment 1",
+  alpha = 0.05
+)
+}
+\keyword{internal}
diff --git a/man/rmst_two_arm.Rd b/man/rmst_two_arm.Rd
new file mode 100644
index 00000000..ac71d75d
--- /dev/null
+++ b/man/rmst_two_arm.Rd
@@ -0,0 +1,56 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/rmst.R
+\name{rmst_two_arm}
+\alias{rmst_two_arm}
+\title{Calculate RMST difference}
+\usage{
+rmst_two_arm(
+  time_var,
+  event_var,
+  group_var,
+  trunc_time,
+  reference = sort(unique(group_var))[1],
+  alpha = 0.05
+)
+}
+\arguments{
+\item{time_var}{A numeric vector of follow up time.}
+
+\item{event_var}{A numeric or integer vector of the status indicator;
+0=alive 1=event.}
+
+\item{group_var}{A vector of treatment groups.}
+
+\item{trunc_time}{A numeric vector of pre-defined cut-off time point(s).}
+
+\item{reference}{Group name of reference group for RMST comparison.
+Default is the first group name by alphabetical order.}
+
+\item{alpha}{A numeric value of the significant level for RMST
+confidence interval. Default is 0.05.}
+}
+\value{
+A list of 2 data frames of RMST calculations:
+\itemize{
+\item \code{rmst_per_arm}: the calculation results per group.
+\item \code{rmst_diff}: the calculation results of RMST differences.
+}
+}
+\description{
+Calculate RMST difference
+}
+\examples{
+data(ex1_delayed_effect)
+with(
+  ex1_delayed_effect,
+  simtrial:::rmst_two_arm(
+    time_var = month,
+    event_var = evntd,
+    group_var = trt,
+    trunc_time = 6,
+    reference = "0",
+    alpha = 0.05
+  )
+)
+}
+\keyword{internal}
diff --git a/vignettes/rmst.Rmd b/vignettes/rmst.Rmd
new file mode 100644
index 00000000..7d98e1c9
--- /dev/null
+++ b/vignettes/rmst.Rmd
@@ -0,0 +1,155 @@
+---
+title: "Restricted mean survival time (RMST)"
+output: rmarkdown::html_vignette
+bibliography: simtrial.bib
+vignette: >
+  %\VignetteIndexEntry{Restricted mean survival time (RMST)}
+  %\VignetteEngine{knitr::rmarkdown}
+---
+
+```{r, include=FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+# Introduction
+
+Restricted mean survival time (RMST) is defined as the area under the survival curve up to a specific time point.  It can be interpreted as the average survival time during a defined time period ranging from time 0 to a specific follow-up time point, which is a straightforward and clinically meaningful way to interpret the contrast in survival between groups.
+
+The RMST may provide valuable information for comparing two survival curves when the proportional hazards assumption is not met, such as in cases of *crossing* or *delayed separation* of survival curves.
+
+# RMST vs. logrank
+
+The log-rank test calculates the test statistics using the survival rate at each time point, and then summarizes them to test the equality of the survival curves as a whole for the entire follow-up period.
+
+A comparison of the RMST between two survival curves provides an estimate of the duration of time gained or lost associated with the exposure.
+
+Although RMST has an advantage over the hazard ratio, a previous study showed that the difference in RMST often has operating characteristics similar to the log-rank test under the proportional hazards assumption [@royston2013restricted]. However, in the case of crossing survival curves, the efficacy of an intervention may be demonstrated by showing a difference in RMST between the two curves, although the log-rank test may fail to detect a difference because of the occurrence of non-proportional hazards.
+
+# Estimation of RMST in a single arm at a single time point
+
+Assume the event time as $T$, and the survival function as $S(t) = Pr(T>t)$. The restricted mean survival time at a pre-specified cutoff time point $\tau$ is
+$$
+  \text{RMST}(\tau) = E[\min (T, \tau)] = \int_{0}^{\tau} S(u) d u.
+$$
+Suppose there are $D$ events, with distinct observed event times at $t_1 < t_2 < \ldots <t_D$. For any $i = 1, \ldots, D$, let $Y_i$ be the number of risk just prior to $t_i$, and let $d_i$ be the number of subjects that fail at $t_i$. The Kaplan-Meier (product-limit) estimate of the survival function at $t_i$ is
+
+$$
+  \hat{S}(t_i)
+  =
+  \prod_{j=1}^{i}
+  \left(
+    1-\frac{d_{j}}{Y_{j}}
+  \right)
+$$
+Based on the definition and formula above, $\text{RMST}(\tau)$ can be estimated by
+$$
+  \widehat{\text{RMST}}(\tau)
+  =
+  \int_{0}^{\tau} \hat{S}(t) d t
+  =
+  \sum_{i=1}^{L_{\tau}} \hat{S}\left(t_{i-1}\right)\left(t_{i}-t_{i-1}\right)
+  +
+  \hat{S}\left(t_{L_{\tau}}\right)\left(\tau-t_{L_{\tau}}\right),
+$$
+where $L_{\tau}$ is the number of $t_i$ values that are less than $\tau$.
+
+The standard error of $\widehat{\text{RMST}}(\tau)$ can be estimated by
+$$
+  \hat{\sigma}
+  =
+  \widehat{\text{Var}}(\widehat{\text{RMST}}(\tau))
+  =
+  \sqrt{\sum_{i=1}^{L_\tau}
+  \frac{d_{i} A_{i}^{2}}{Y_{i}\left(Y_{i}-d_{i}\right)}}
+$$
+where
+$$
+  A_{i}
+  =
+  \int_{t_i}^{\tau} \hat{S}(t) d t
+  =
+  \sum_{j=i}^{L_\tau}
+  \hat{S}(t_{j}) (t_{j+1}-t_{j})
+  +
+  \hat{S}(t_{L_\tau})(\tau-t_{L_\tau})
+$$
+and $m=\sum_{j=1}^{L_\tau} d_{j}$. The $(1-\alpha)$ confidence interval for $\text{RMST}$ can be calculated as
+$$
+  \left[
+    \widehat{\operatorname{RMST}}(\tau) - z_{\alpha/2}\hat{\sigma},
+    \;\;
+    \widehat{\operatorname{RMST}}(\tau) + z_{\alpha/2}\hat{\sigma}
+  \right]
+$$
+where $\alpha$ is a predefined significant level, and $z_{\alpha/2}$ is the upper $1-\alpha/2$ critical value for the standard normal distribution.
+
+```{r, message=FALSE, warning=FALSE}
+# Simulate NPH data from the piecewise model
+library(simtrial)
+# Table display
+library(gt)
+```
+
+```{r}
+data(ex1_delayed_effect)
+data_single_arm <- ex1_delayed_effect[ex1_delayed_effect$trt == 1, ]
+simtrial:::rmst_single_arm(
+  time_var = data_single_arm$month,
+  event_var = data_single_arm$evntd,
+  tau = 10
+) |> gt()
+```
+
+## Estimation of RMST differences in 2 arms at a single time point
+
+Let $\text{RMST}_{1}(\tau)$ and $\text{RMST}_{2}(\tau)$ be the RMST of treatment group 1 and 2 at predefined time $\tau$, the RMST difference between the 2 treatment groups ($\theta$) can be defined as
+$$
+  \theta = \text{RMST}_1(\tau) - \text{RMST}_2(\tau).
+$$
+The expected value of $\theta$ is $E(\theta) = E[RMST_{1}(\tau)] - E[RMST_{2}(\tau)]$. If two treatment groups are independent, the variance of $\theta$ is
+$$
+  \text{Var}(\theta) = \sigma_{1}^{2} + \sigma_{2}^{2}
+$$
+
+Similarly, the $(1-\alpha)$ confidence interval for $\text{RMST}$ difference between 2 groups can be calculated as:
+$$
+  \left[
+    \hat{\theta} - z_{\alpha/2}\sqrt{\hat{\sigma}_1^2 + \hat{\sigma}_2^2},
+    \;\;
+    \hat{\theta} + z_{\alpha/2}\sqrt{\hat{\sigma}_1^2 + \hat{\sigma}_2^2}
+  \right].
+$$
+
+```{r}
+tau <- 10
+
+data(ex1_delayed_effect)
+
+ex1_delayed_effect |>
+  rmst(
+    var_label_tte = "month",
+    var_label_event = "evntd",
+    var_label_group = "trt",
+    tau = 10,
+    reference = "0"
+  )
+```
+
+The R package survRM2 [@uno2022survrm2] performs two-sample comparisons using RMST as a summary measure of the survival time distribution. Three kinds of between-group contrast metrics (i.e., the difference in RMST, the ratio of RMST and the ratio of the restricted mean time lost (RMTL)) are computed. It performs an ANCOVA-type covariate adjustment as well as unadjusted analyses for those measures. We use this R package as validation of `simtrial::rmst()`.
+
+```{r}
+verify <- survRM2::rmst2(
+  time = ex1_delayed_effect$month,
+  status = ex1_delayed_effect$evntd,
+  arm = ex1_delayed_effect$trt,
+  tau = tau,
+  alpha = 0.05
+)
+
+verify$RMST.arm1$rmst[1] - verify$RMST.arm0$rmst[1]
+```
+
+## References
diff --git a/vignettes/simtrial.bib b/vignettes/simtrial.bib
index 2eba68c4..d6afb79d 100644
--- a/vignettes/simtrial.bib
+++ b/vignettes/simtrial.bib
@@ -187,6 +187,17 @@ @article{NPHWGSimulation
   publisher = {Taylor \& Francis}
 }
 
+@article{royston2013restricted,
+  title     = {Restricted mean survival time: an alternative to the hazard ratio for the design and analysis of randomized trials with a time-to-event outcome},
+  author    = {Royston, Patrick and Parmar, Mahesh KB},
+  journal   = {BMC Medical Research Methodology},
+  volume    = {13},
+  number    = {1},
+  pages     = {1--15},
+  year      = {2013},
+  publisher = {BioMed Central}
+}
+
 @article{Schemper2009,
   title   = {The estimation of average hazard ratios by weighted Cox regression},
   author  = {Schemper, Michael and Wakounig, Samo and Heinze, Georg},
@@ -206,6 +217,14 @@ @article{Tsiatis
   year    = {1982}
 }
 
+@article{uno2022survrm2,
+  title  = {{survRM2}: Comparing Restricted Mean Survival Time},
+  author = {Hajime Uno and Lu Tian and Miki Horiguchi and Angel Cronin and Chakib Battioui and James Bell},
+  year   = {2022},
+  note   = {R package version 1.0-4},
+  url    = {https://CRAN.R-project.org/package=survRM2}
+}
+
 @article{Xi2017,
   title     = {A unified framework for weighted parametric multiple test procedures},
   author    = {Xi, Dong and Glimm, Ekkehard and Maurer, Willi and Bretz, Frank},
@@ -215,4 +234,4 @@ @article{Xi2017
   pages     = {918--931},
   year      = {2017},
   publisher = {Wiley Online Library}
-}
\ No newline at end of file
+}