diff --git a/NAMESPACE b/NAMESPACE
index 8f3872a..d1acc24 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -3,8 +3,8 @@
 S3method(print,or_ci)
 S3method(print,prop_ci)
 S3method(print,s_lsmeans)
-export("%>%")
 export(derive_bor)
+export(h_pairwise_survdiff)
 export(h_prep_prop)
 export(rrPostProb)
 export(s_get_coxph)
@@ -19,16 +19,20 @@ importFrom(dplyr,case_when)
 importFrom(dplyr,count)
 importFrom(dplyr,distinct)
 importFrom(dplyr,filter)
+importFrom(dplyr,full_join)
 importFrom(dplyr,group_by)
 importFrom(dplyr,left_join)
 importFrom(dplyr,mutate)
 importFrom(dplyr,row_number)
+importFrom(dplyr,rowwise)
 importFrom(dplyr,select)
+importFrom(dplyr,summarise)
 importFrom(dplyr,ungroup)
 importFrom(lifecycle,deprecated)
 importFrom(lubridate,days)
 importFrom(lubridate,ymd)
 importFrom(magrittr,"%>%")
+importFrom(magrittr,set_colnames)
 importFrom(rlang,":=")
 importFrom(rlang,.data)
 importFrom(rlang,sym)
diff --git a/R/package.R b/R/package.R
index d62a6aa..e6e589d 100644
--- a/R/package.R
+++ b/R/package.R
@@ -5,10 +5,11 @@
 "_PACKAGE"
 
 #' @import checkmate
+#' @importFrom magrittr %>% set_colnames
 #' @importFrom lifecycle deprecated
 #' @importFrom stats pbeta rbinom confint as.formula setNames coef quantile
-#' @importFrom dplyr add_count arrange case_when count distinct filter group_by
-#'  left_join mutate row_number select ungroup
+#' @importFrom dplyr add_count arrange case_when count distinct filter full_join
+#'  group_by left_join mutate row_number rowwise select summarise ungroup
 #' @importFrom rlang sym := .data
 #' @importFrom survival coxph Surv strata
 #' @importFrom lubridate ymd days
diff --git a/R/summarize-survival.R b/R/summarize-survival.R
index 673c0b8..54766af 100644
--- a/R/summarize-survival.R
+++ b/R/summarize-survival.R
@@ -21,8 +21,9 @@
 #' @param conf_level (`numeric`)\cr significance level for a two-side confidence
 #'  interval on survival curve, default is 0.95.
 #' @param strata (`character`)\cr a character vector used for stratification.
-#' @param pval_method (`string`)\cr method of comparing survival curves. Default is
-#'  "log-rank", others options can see [survminer::surv_pvalue()].
+#' @param rho (`number`)\cr a scalar parameter that controls the method
+#'  of comparing survival curves. Default is 0 that is log-rank , others options
+#'  can see [survival::survdiff()].
 #' @param pairwise (`logical`)\cr whether to conduct the pairwise comparison.
 #' @param ... other arguments to be passed to [survival::survfit()].
 #'
@@ -37,7 +38,7 @@
 #'
 #' @note
 #' This function mainly wraps `survival::survfit()` to compute the survival
-#' estimates, and `survminer::surv_pvalue()` to compute the p-value of comparing
+#' estimates, and `survival::survdiff()` to compute the p-value of comparing
 #' survival curves. The SE and p-value of difference rate are simply computed
 #' using normal distribution and Z statistic.
 #'
@@ -101,7 +102,7 @@ s_get_survfit <- function(data,
                           conf_type = c("log-log", "log", "plain"),
                           conf_level = 0.95,
                           strata = NULL,
-                          pval_method = "log-rank",
+                          rho = 0,
                           pairwise = FALSE,
                           ...) {
   assert_class(data, "data.frame")
@@ -110,7 +111,7 @@ s_get_survfit <- function(data,
   assert_numeric(time_point, null.ok = TRUE)
   assert_number(conf_level, lower = 0, upper = 1)
   assert_subset(strata, names(data))
-  assert_string(pval_method)
+  assert_number(rho)
   assert_logical(pairwise)
   conf_type <- match.arg(conf_type, c("log-log", "log", "plain"), several.ok = FALSE)
 
@@ -155,38 +156,6 @@ s_get_survfit <- function(data,
       select(group, n = n.start, events, median, lower = `0.95LCL`, upper = `0.95UCL`)
   }
 
-  # test survival curves
-  surv_test <- if (!is.null(km_fit$strata)) {
-    if (is.null(strata)) {
-      km_no_strata <- do.call(
-        survival::survfit,
-        args = list(
-          data = data,
-          formula = formula,
-          conf.int = conf_level,
-          conf.type = conf_type
-        )
-      )
-      survminer::surv_pvalue(km_no_strata, data = data, method = pval_method) %>%
-        tibble::as_tibble()
-    } else {
-      formula_strata <- as.formula(
-        paste0(format(formula), " + strata(", paste(strata, collapse = ", "), ")")
-      )
-      km_strata <- do.call(
-        survival::survfit,
-        args = list(
-          data = data,
-          formula = formula_strata,
-          conf.int = conf_level,
-          conf.type = conf_type
-        )
-      )
-      survminer::surv_pvalue(km_strata, data = data, method = pval_method) %>%
-        tibble::as_tibble()
-    }
-  }
-
   # overall survival rate
   cols <- c("time", "n.risk", "n.event", "n.censor", "surv", "std.err", "lower", "upper")
   overall <- summary(km_fit)
@@ -194,7 +163,7 @@ s_get_survfit <- function(data,
     tibble::as_tibble(overall[cols]) %>%
       mutate(group = grps)
   } else {
-    evt_ind <- which(overall$time < lag(overall$time)) - 1
+    evt_ind <- which(overall$time < dplyr::lag(overall$time)) - 1
     evt_ind <- c(evt_ind[1], diff(c(evt_ind, length(overall$time))))
     tibble::as_tibble(overall[cols]) %>%
       mutate(group = rep(grps, times = evt_ind))
@@ -243,6 +212,28 @@ s_get_survfit <- function(data,
       purrr::list_rbind()
   }
 
+  # test survival curves
+  surv_test <- if (!is.null(km_fit$strata)) {
+    survdiff <- h_pairwise_survdiff(
+      formula = formula,
+      strata = strata,
+      data = data,
+    )
+    if (!pairwise) {
+      tibble::tibble(
+        comparsion = paste(bylist[,1], bylist[,2], sep = " vs. "),
+        method = survdiff$method,
+        pval = as.numeric(na.omit(survdiff$p.value))
+      )
+    } else {
+      tibble::tibble(
+        comparsion = paste(bylist[,1], bylist[,2], sep = " vs. "),
+        method = survdiff$method,
+        pval = as.numeric(survdiff$p.value)[!is.na(as.numeric(survdiff$p.value))]
+      )
+    }
+  }
+
   # survival range in overall time
   range_evt <- surv_overall_rate %>%
     filter(n.event != 0) %>%
@@ -265,7 +256,8 @@ s_get_survfit <- function(data,
     mutate(
       min = min(event_min, censor_min, na.rm = TRUE),
       max = max(event_max, censor_max, na.rm = TRUE)
-    )
+    ) %>%
+    ungroup()
 
   structure(
     list(
@@ -281,13 +273,13 @@ s_get_survfit <- function(data,
         rate = surv_rate_diff
       ),
       params = list(
-        formula = format(c(formula, formula_strata)),
+        formula = format(formula),
         quantile = quantile,
         time_point = time_point,
         conf_type = conf_type,
         conf_level = conf_level,
         strata = strata,
-        pval_method = pval_method,
+        rho = rho,
         pairwise = pairwise
       )
     ),
diff --git a/R/utils.R b/R/utils.R
index fd0b1d1..be16abb 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -1,14 +1,59 @@
-#' Pipe operator
+#' @describeIn s_get_survfit a modified version of [survminer::pairwise_survdiff()]
+#' that can calculate the pairwise comparisons of survival curves regardless of
+#' whether stratification is given.
 #'
-#' See \code{magrittr::\link[magrittr:pipe]{\%>\%}} for details.
-#'
-#' @name %>%
-#' @rdname pipe
-#' @keywords internal
 #' @export
-#' @importFrom magrittr %>%
-#' @usage lhs \%>\% rhs
-#' @param lhs A value or the magrittr placeholder.
-#' @param rhs A function call using the magrittr semantics.
-#' @return The result of calling `rhs(lhs)`.
-NULL
+h_pairwise_survdiff <- function(formula,
+                                data,
+                                strata = NULL,
+                                rho = 0) {
+  group_var <- attr(stats::terms(formula), "term.labels")
+  formula <- if (!is.null(strata)) {
+    as.formula(
+      paste0(format(formula), " + strata(", paste(strata, collapse = ", "), ")")
+    )
+  } else {
+    formula
+  }
+
+  method <- if (rho == 0) {
+    "Log-Rank"
+  } else if (rho == 1) {
+    "Peto & Peto"
+  }
+
+  vardf <- unlist(data[, group_var])
+  group <- if (!is.factor(vardf)) {
+    droplevels(as.factor(vardf))
+  } else {
+    droplevels(vardf)
+  }
+
+  compare.levels <- function(i, j) {
+    subdf <- data %>% filter(!!sym(group_var) %in% c(levels(group)[c(i, j)]))
+    sdif <- survival::survdiff(formula, data = subdf, rho = rho)
+    stats::pchisq(sdif$chisq, length(sdif$n) - 1, lower.tail = FALSE)
+  }
+
+  level.names <- levels(group)
+  ix <- setNames(seq_along(level.names), paste0(group_var, "=", level.names))
+  pval <- outer(
+    ix[-1L], ix[-length(ix)],
+    function(ivec, jvec) {
+      vapply(seq_along(ivec), function(k) {
+        i <- ivec[k]
+        j <- jvec[k]
+        if (i > j) {
+          compare.levels(i, j)
+        } else {
+          NA_real_
+        }
+      }, FUN.VALUE = 0.05)
+    }
+  )
+
+  list(
+    method = ifelse(!is.null(strata), paste("Stratified", method), method),
+    p.value = pval
+  )
+}
diff --git a/man/pipe.Rd b/man/pipe.Rd
deleted file mode 100644
index 5fa90fe..0000000
--- a/man/pipe.Rd
+++ /dev/null
@@ -1,20 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/utils.R
-\name{\%>\%}
-\alias{\%>\%}
-\title{Pipe operator}
-\usage{
-lhs \%>\% rhs
-}
-\arguments{
-\item{lhs}{A value or the magrittr placeholder.}
-
-\item{rhs}{A function call using the magrittr semantics.}
-}
-\value{
-The result of calling \code{rhs(lhs)}.
-}
-\description{
-See \code{magrittr::\link[magrittr:pipe]{\%>\%}} for details.
-}
-\keyword{internal}
diff --git a/man/s_get_survfit.Rd b/man/s_get_survfit.Rd
index e9c58e9..2f2f5bd 100644
--- a/man/s_get_survfit.Rd
+++ b/man/s_get_survfit.Rd
@@ -1,7 +1,8 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/summarize-survival.R
+% Please edit documentation in R/summarize-survival.R, R/utils.R
 \name{s_get_survfit}
 \alias{s_get_survfit}
+\alias{h_pairwise_survdiff}
 \title{Summarize Survival Curves Analysis}
 \usage{
 s_get_survfit(
@@ -12,10 +13,12 @@ s_get_survfit(
   conf_type = c("log-log", "log", "plain"),
   conf_level = 0.95,
   strata = NULL,
-  pval_method = "log-rank",
+  rho = 0,
   pairwise = FALSE,
   ...
 )
+
+h_pairwise_survdiff(formula, data, strata = NULL, rho = 0)
 }
 \arguments{
 \item{data}{(\code{data.frame})\cr a data frame as input.}
@@ -37,8 +40,9 @@ interval on survival curve, default is 0.95.}
 
 \item{strata}{(\code{character})\cr a character vector used for stratification.}
 
-\item{pval_method}{(\code{string})\cr method of comparing survival curves. Default is
-"log-rank", others options can see \code{\link[survminer:surv_pvalue]{survminer::surv_pvalue()}}.}
+\item{rho}{(\code{number})\cr a scalar parameter that controls the method
+of comparing survival curves. Default is 0 that is log-rank , others options
+can see \code{\link[survival:survdiff]{survival::survdiff()}}.}
 
 \item{pairwise}{(\code{logical})\cr whether to conduct the pairwise comparison.}
 
@@ -62,9 +66,16 @@ statistical outputs, such as median and percentile survival time, survival
 time range, testing survival curve differences with or without stratification,
 time-point survival rates and difference of survival rates.
 }
+\section{Functions}{
+\itemize{
+\item \code{h_pairwise_survdiff()}: a modified version of \code{\link[survminer:pairwise_survdiff]{survminer::pairwise_survdiff()}}
+that calculates the pairwise comparisons of survival curves when the stratification
+is specified.
+
+}}
 \note{
 This function mainly wraps \code{survival::survfit()} to compute the survival
-estimates, and \code{survminer::surv_pvalue()} to compute the p-value of comparing
+estimates, and \code{survival::survdiff()} to compute the p-value of comparing
 survival curves. The SE and p-value of difference rate are simply computed
 using normal distribution and Z statistic.