diff --git a/NAMESPACE b/NAMESPACE
index 16ddbdd..f7ffbaa 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -19,6 +19,7 @@ export(s_propci)
 import(checkmate)
 importFrom(dplyr,add_count)
 importFrom(dplyr,arrange)
+importFrom(dplyr,bind_rows)
 importFrom(dplyr,case_when)
 importFrom(dplyr,count)
 importFrom(dplyr,distinct)
diff --git a/R/package.R b/R/package.R
index b37d535..32ce4bd 100644
--- a/R/package.R
+++ b/R/package.R
@@ -8,8 +8,8 @@
 #' @importFrom magrittr set_colnames set_rownames
 #' @importFrom lifecycle deprecated
 #' @importFrom stats pbeta rbinom confint as.formula setNames coef quantile
-#' @importFrom dplyr add_count arrange case_when count distinct filter full_join
-#'  group_by left_join mutate row_number rowwise select summarise ungroup
+#' @importFrom dplyr add_count arrange bind_rows case_when count distinct filter
+#'  full_join group_by left_join mutate row_number rowwise select summarise ungroup
 #' @importFrom rlang sym := .data
 #' @importFrom rtables analyze basic_table build_table in_rows rcell non_ref_rcell
 #'  split_cols_by
diff --git a/R/pkg-methods.R b/R/pkg-methods.R
index 8f65410..b6a822c 100644
--- a/R/pkg-methods.R
+++ b/R/pkg-methods.R
@@ -5,32 +5,83 @@ print.s_lsmeans <- function(x, ...) {
   cat("Model Call: ", append = FALSE)
   print(x$model$call)
   cat("\n")
-  varc <- paste0(x$model$xlev[[x$params$var]], collapse = ", ")
-  cat("Predictor/Treatment: ", x$params$var, " (", varc, ")\n", sep = "")
-  if (!is.null(x$params$by)) {
-    byc <- paste0(x$model$xlev[[x$params$by]], collapse = ", ")
-    cat("Group by: ", x$params$by, " (", byc, ")\n", sep = "")
+
+  grp_var <- x$params$var
+  by_var <- x$params$by
+  df <- x$data
+  ref_col <- levels(x$lsm_est[[grp_var]])[1]
+  side <- if (x$params$alternative == "greater" & x$params$null >= 0) {
+    "p-value (Superiority Test)"
+  } else if (x$params$alternative == "greater" & x$params$null < 0) {
+    "p-value (Non-inferiority Test)"
+  } else if (x$params$alternative == "less" & x$params$null >= 0) {
+    "p-value (Non-Superiority Test)"
+  } else if (x$params$alternative == "less" & x$params$null < 0) {
+    "p-value (inferiority Test)"
+  } else {
+    "p-value"
   }
 
-  cat("\n")
-  cat("Least-squares Means Estimates:\n")
-  print(x$lsm_est)
+  df <- bind_rows(
+    x$lsm_est,
+    x$lsm_contr
+  ) %>%
+    rowwise() %>%
+    mutate(
+      !!sym(grp_var) := case_when(
+        is.na(!!sym(grp_var)) ~ strsplit(as.character(.data$contrast), " - ")[[1]][1],
+        TRUE ~ !!sym(grp_var)
+      )
+    )
 
-  if (x$params$contrast) {
-    cat("\n")
-    cat("Contrast Estimates of Least-squares Means:\n")
-    if (x$params$alternative == "two.sided") {
-      cat(sprintf("Null hypothesis is \u03b8 equal to %s.\n", x$params$null))
-    } else if (x$params$alternative == "greater" & x$params$null >= 0) {
-      cat(sprintf("Null hypothesis is \u03b8 non-superiority to %s.\n", x$params$null))
-    } else if (x$params$alternative == "greater" & x$params$null < 0) {
-      cat(sprintf("Null hypothesis is \u03b8 inferiority to %s.\n", x$params$null))
-    } else if (x$params$alternative == "less" & x$params$null >= 0) {
-      cat(sprintf("Null hypothesis is \u03b8 superiority to %s.\n", x$params$null))
-    }
-    print(x$lsm_contr)
+  a_lsmeam_func <- function(df, .var, .in_ref_col, est, contr, vst) {
+    curgrp <- df[[.var]][1]
+    df_est <- df %>%
+      filter(!!sym(grp_var) == curgrp & is.na(.data$contrast))
+    df_contr <- df %>%
+      filter(!!sym(grp_var) == curgrp & .data$contrast == paste(curgrp, "-", ref_col))
+    in_rows(
+      rcell(unlist(df_est[, c("estimate", "SE"), drop = TRUE]), format = "xx.xx (xx.xx)"),
+      rcell(unlist(df_est[, c("lower.CL", "upper.CL"), drop = TRUE]),
+        format = "(xx.xx, xx.xx)",
+        indent_mod = 1L
+      ),
+      non_ref_rcell(
+        unlist(df_contr[, c("estimate", "SE"), drop = TRUE]),
+        .in_ref_col,
+        format = "xx.xx (xx.xx)"
+      ),
+      non_ref_rcell(
+        unlist(df_contr[, c("lower.CL", "upper.CL"), drop = TRUE]),
+        .in_ref_col,
+        format = "(xx.xx, xx.xx)",
+        indent_mod = 1L
+      ),
+      non_ref_rcell(
+        df_contr[, "p.value", drop = TRUE],
+        .in_ref_col,
+        format = "x.xxxx | (<0.0001)",
+        indent_mod = 1L
+      ),
+      .names = c(
+        "LS Mean (SE)", "95% CI",
+        "Difference LS Mean (SE)", "95% CI",
+        side
+      )
+    )
   }
 
+  tbl <- basic_table() %>%
+    split_cols_by(grp_var, ref_group = ref_col)
+  if (!is.null(by_var)) {
+    tbl <- tbl %>% rtables::split_rows_by(by_var, label_pos = "topleft")
+  }
+  result <- tbl %>%
+    analyze(grp_var, a_lsmeam_func) %>%
+    rtables::append_topleft("  Statistics") %>%
+    build_table(df = df)
+  print(result)
+
   invisible(x)
 }
 
diff --git a/R/summarize-lsmeans.R b/R/summarize-lsmeans.R
index 45b68ea..3eed71c 100644
--- a/R/summarize-lsmeans.R
+++ b/R/summarize-lsmeans.R
@@ -15,7 +15,8 @@
 #' @param object (`fitted model`)\cr any fitted model that can be accepted by
 #'  `emmeans::emmeans()`, such as object from `lm` for ANCOVA and `mmrm` for MMRM.
 #' @param var (`string`)\cr string specifying the name of the predictor, such as
-#'  it would be treatment group variable (TRT01PN) in the clinical trials.
+#'  it would be treatment group variable (TRT01PN) in the clinical trials. Default
+#'  the first level is defined as reference group.
 #' @param by (`string`)\cr string specifying the name of the grouped variable.
 #'  The estimates and contrasts will be evaluated separately for each elements
 #'  from the grouped variable.
@@ -74,25 +75,20 @@
 #' @examples
 #' # refactor the level order:
 #' data(fev_data, package = "mmrm")
-#' fev_data$ARMCD <- factor(fev_data$ARMCD, level = c("TRT", "PBO"))
 #'
-#' # fit mmrm model:
+#' # fit MMRM model:
 #' fit <- mmrm::mmrm(
 #'   formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
 #'   reml = TRUE, method = "Kenward-Roger", vcov = "Kenward-Roger-Linear",
 #'   data = fev_data
 #' )
 #'
-#' # estimate ARMCD within whole visits:
-#' s_get_lsmeans(fit, "ARMCD")
-#'
 #' # estimate ARMCD by visits:
 #' s_get_lsmeans(fit, "ARMCD", by = "AVISIT")
 #'
 #' # estimate ARMCD by visits for superiority testing with null hypothesis of 2:
 #' s_get_lsmeans(fit, "ARMCD", by = "AVISIT", null = 2, alternative = "greater")
 #'
-#'
 #' # fit ANCOVA model:
 #' fit2 <- fev_data |>
 #'   dplyr::filter(VISITN == 4 & !is.na(FEV1)) |>
@@ -132,7 +128,7 @@ s_get_lsmeans <- function(object,
     )
 
   contr <- if (contrast) {
-    emmeans::contrast(ems, adjust = "none", method = "pairwise")
+    emmeans::contrast(ems, adjust = "none", method = "revpairwise")
   }
   contr_ci <- tibble::as_tibble(
     data.frame(confint(contr, level = conf.level))
@@ -150,6 +146,7 @@ s_get_lsmeans <- function(object,
   structure(
     list(
       model = ems@model.info,
+      data = object$data,
       lsm_est = lsm_est_res,
       lsm_contr = contr_est_res,
       params = list(
diff --git a/man/s_get_lsmeans.Rd b/man/s_get_lsmeans.Rd
index 5e27294..14c7e66 100644
--- a/man/s_get_lsmeans.Rd
+++ b/man/s_get_lsmeans.Rd
@@ -23,7 +23,8 @@ s_get_lsmeans(
 \code{emmeans::emmeans()}, such as object from \code{lm} for ANCOVA and \code{mmrm} for MMRM.}
 
 \item{var}{(\code{string})\cr string specifying the name of the predictor, such as
-it would be treatment group variable (TRT01PN) in the clinical trials.}
+it would be treatment group variable (TRT01PN) in the clinical trials. Default
+the first level is defined as reference group.}
 
 \item{by}{(\code{string})\cr string specifying the name of the grouped variable.
 The estimates and contrasts will be evaluated separately for each elements
@@ -83,25 +84,20 @@ There is no any p value adjusted method involved.
 \examples{
 # refactor the level order:
 data(fev_data, package = "mmrm")
-fev_data$ARMCD <- factor(fev_data$ARMCD, level = c("TRT", "PBO"))
 
-# fit mmrm model:
+# fit MMRM model:
 fit <- mmrm::mmrm(
   formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
   reml = TRUE, method = "Kenward-Roger", vcov = "Kenward-Roger-Linear",
   data = fev_data
 )
 
-# estimate ARMCD within whole visits:
-s_get_lsmeans(fit, "ARMCD")
-
 # estimate ARMCD by visits:
 s_get_lsmeans(fit, "ARMCD", by = "AVISIT")
 
 # estimate ARMCD by visits for superiority testing with null hypothesis of 2:
 s_get_lsmeans(fit, "ARMCD", by = "AVISIT", null = 2, alternative = "greater")
 
-
 # fit ANCOVA model:
 fit2 <- fev_data |>
   dplyr::filter(VISITN == 4 & !is.na(FEV1)) |>
diff --git a/tests/spelling.R b/tests/spelling.R
index 6713838..13f77d9 100644
--- a/tests/spelling.R
+++ b/tests/spelling.R
@@ -1,3 +1,6 @@
-if(requireNamespace('spelling', quietly = TRUE))
-  spelling::spell_check_test(vignettes = TRUE, error = FALSE,
-                             skip_on_cran = TRUE)
+if (requireNamespace("spelling", quietly = TRUE)) {
+  spelling::spell_check_test(
+    vignettes = TRUE, error = FALSE,
+    skip_on_cran = TRUE
+  )
+}
diff --git a/tests/testthat/test-pkg-methods.R b/tests/testthat/test-pkg-methods.R
index d93c4df..3a6bfed 100644
--- a/tests/testthat/test-pkg-methods.R
+++ b/tests/testthat/test-pkg-methods.R
@@ -7,16 +7,14 @@ test_that("print.s_lsmeans works as expected for mmrm", {
   )
   res <- capture_output(print(s_get_lsmeans(fit, "ARMCD", by = "AVISIT")))
   expect_match(res, "Model Call: mmrm::mmrm(formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT", fixed = TRUE)
-  expect_match(res, "Predictor/Treatment: ARMCD (PBO, TRT)", fixed = TRUE)
-  expect_match(res, "Group by: AVISIT (VIS1, VIS2, VIS3, VIS4)", fixed = TRUE)
-  expect_match(res, "Least-squares Means Estimates", fixed = TRUE)
-  expect_match(res, "Null hypothesis is θ equal to 0", fixed = TRUE)
+  expect_match(res, "LS Mean (SE)", fixed = TRUE)
+  expect_match(res, "Difference LS Mean (SE)", fixed = TRUE)
 
   res2 <- capture_output(print(s_get_lsmeans(fit, "ARMCD", by = "AVISIT", null = 2, alternative = "greater")))
-  expect_match(res2, "Null hypothesis is θ non-superiority to 2", fixed = TRUE)
+  expect_match(res2, "p-value (Superiority Test)", fixed = TRUE)
 
   res3 <- capture_output(print(s_get_lsmeans(fit, "ARMCD", by = "AVISIT", null = -2, alternative = "greater")))
-  expect_match(res3, "Null hypothesis is θ inferiority to -2", fixed = TRUE)
+  expect_match(res3, "p-value (Non-inferiority Test)", fixed = TRUE)
 })
 
 test_that("print.prop_ci works as expected", {
diff --git a/tests/testthat/test-summarize-lsmeans.R b/tests/testthat/test-summarize-lsmeans.R
index a4107d3..92ed704 100644
--- a/tests/testthat/test-summarize-lsmeans.R
+++ b/tests/testthat/test-summarize-lsmeans.R
@@ -28,14 +28,14 @@ test_that("s_get_lsmeans works as expected by visit with default arguments", {
   expect_equal(
     res$lsm_contr[1:2, ],
     tibble::tibble(
-      contrast = factor(c("PBO - TRT", "PBO - TRT")),
+      contrast = factor(c("TRT - PBO", "TRT - PBO")),
       AVISIT = factor(c("VIS1", "VIS2"), levels = c("VIS1", "VIS2", "VIS3", "VIS4")),
-      estimate = c(-3.774230, -3.732304),
+      estimate = c(3.774230, 3.732304),
       SE = c(1.0817635, 0.8633334),
       df = c(145.5520, 145.2771),
-      lower.CL = c(-5.912224, -5.438620),
-      upper.CL = c(-1.636236, -2.025988),
-      t.ratio = c(-3.488960, -4.323131),
+      lower.CL = c(1.636236, 2.025988),
+      upper.CL = c(5.912224, 5.438620),
+      t.ratio = c(3.488960, 4.323131),
       p.value = c(6.417790e-04, 2.840654e-05)
     ),
     tolerance = 0.0001
@@ -44,7 +44,6 @@ test_that("s_get_lsmeans works as expected by visit with default arguments", {
 
 test_that("s_get_lsmeans works as expected for superiority testing with null hypothesis of 2", {
   data(fev_data, package = "mmrm")
-  fev_data$ARMCD <- factor(fev_data$ARMCD, level = c("TRT", "PBO"))
   fit <- mmrm::mmrm(
     formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
     reml = TRUE, method = "Kenward-Roger", vcov = "Kenward-Roger-Linear",
@@ -97,13 +96,13 @@ test_that("s_get_lsmeans works as expected for ancova", {
   expect_equal(
     res$lsm_contr,
     tibble::tibble(
-      contrast = c("PBO - TRT"),
-      estimate = c(-3.89943),
+      contrast = c("TRT - PBO"),
+      estimate = c(3.89943),
       SE = c(1.686347),
       df = c(128),
-      lower.CL = c(-7.236155),
-      upper.CL = c(-0.562705),
-      t.ratio = c(-2.312354),
+      lower.CL = c(0.562705),
+      upper.CL = c(7.236155),
+      t.ratio = c(2.312354),
       p.value = c(0.02235472)
     ),
     tolerance = 0.0001