diff --git a/DESCRIPTION b/DESCRIPTION
index 30b7768fb2..47c631beca 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -133,7 +133,6 @@ Collate:
     'prop_diff.R'
     'prop_diff_test.R'
     'prune_occurrences.R'
-    'reflast.R'
     'response_biomarkers_subgroups.R'
     'response_subgroups.R'
     'riskdiff.R'
@@ -160,3 +159,4 @@ Collate:
     'utils_factor.R'
     'utils_grid.R'
     'utils_rtables.R'
+    'utils_split_funs.R'
diff --git a/NAMESPACE b/NAMESPACE
index 3161d23cd3..d21ec2a24f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -217,6 +217,7 @@ export(has_fraction_in_cols)
 export(has_fractions_difference)
 export(imputation_rule)
 export(keep_content_rows)
+export(keep_level_order)
 export(keep_rows)
 export(logistic_regression_cols)
 export(logistic_summary_by_flag)
diff --git a/R/reflast.R b/R/reflast.R
deleted file mode 100644
index 0257fb3f2d..0000000000
--- a/R/reflast.R
+++ /dev/null
@@ -1,69 +0,0 @@
-#' Split Function to Place Reference Group Facet Last
-#'
-#' Place reference group facet last during post-processing stage in a custom split function.
-#'
-#' @param splret result of the core split
-#' @param spl split object
-#' @param fulldf data.frame of incoming data to be split
-#'
-#' @export
-#'
-#' @seealso [rtables::make_split_fun()]
-#'
-#' @examples
-#' library(dplyr)
-#'
-#' # Define custom split function
-#' ref_last <- make_split_fun(post = list(ref_group_last))
-#'
-#' dat <- data.frame(
-#'   x = factor(letters[1:5], levels = letters[5:1]),
-#'   y = 1:5
-#' )
-#'
-#' # With rtables layout functions
-#' basic_table() %>%
-#'   split_cols_by("x", ref_group = "c", split_fun = ref_last) %>%
-#'   analyze("y") %>%
-#'   build_table(dat)
-#'
-#' # With tern layout funcitons
-#' adtte_f <- tern_ex_adtte %>%
-#'   filter(PARAMCD == "OS") %>%
-#'   mutate(
-#'     AVAL = day2month(AVAL),
-#'     is_event = CNSR == 0
-#'   )
-#'
-#' basic_table() %>%
-#'   split_cols_by(var = "ARMCD", ref_group = "ARM B", split_fun = ref_last) %>%
-#'   add_colcounts() %>%
-#'   surv_time(
-#'     vars = "AVAL",
-#'     var_labels = "Survival Time (Months)",
-#'     is_event = "is_event",
-#'   ) %>%
-#'   build_table(df = adtte_f)
-#'
-ref_group_last <- function(splret, spl, fulldf) {
-  if (!"ref_group_value" %in% slotNames(spl)) {
-    stop("Reference group is undefined.")
-  }
-
-  spl_var <- spl@payload
-  init_lvls <- levels(fulldf[[spl_var]])
-
-  if (!all(names(splret$values) %in% init_lvls)) {
-    stop("This split function does not work with combination facets.")
-  }
-
-  ref_group_pos <- which(init_lvls == spl@ref_group_value)
-  reord_lvls <- c(init_lvls[-ref_group_pos], init_lvls[ref_group_pos])
-  ord <- match(reord_lvls, names(splret$values))
-
-  make_split_result(
-    splret$values[ord],
-    splret$datasplit[ord],
-    splret$labels[ord]
-  )
-}
diff --git a/R/utils_split_funs.R b/R/utils_split_funs.R
new file mode 100644
index 0000000000..c57ad1a8a8
--- /dev/null
+++ b/R/utils_split_funs.R
@@ -0,0 +1,108 @@
+#' Custom Split Functions
+#'
+#' @description
+#' Collection of useful functions that are expanding on the core list of functions
+#' provided by `rtables`. See `?rtables::custom_split_funs` and [rtables::make_split_fun()]
+#' for more information on how to make a custom split function. All these functions
+#' work with [split_rows_by()] argument `split_fun` to modify the way the split
+#' happens.
+#'
+#' @seealso [rtables::make_split_fun()]
+#'
+#' @name utils_split_funs
+NULL
+
+#' @describeIn utils_split_funs split function to place reference group facet last
+#'  during post-processing stage.
+#'
+#' @return
+#' * `ref_group_last` returns an utility function that puts the reference group
+#'  as last and needs to be assigned to `split_fun`.
+#'
+#' @examples
+#' library(dplyr)
+#'
+#' dat <- data.frame(
+#'   x = factor(letters[1:5], levels = letters[5:1]),
+#'   y = 1:5
+#' )
+#' # ref_group_last
+#'
+#' # With rtables layout functions
+#' basic_table() %>%
+#'   split_cols_by("x", ref_group = "c", split_fun = ref_group_last) %>%
+#'   analyze("y") %>%
+#'   build_table(dat)
+#'
+#' # With tern layout funcitons
+#' adtte_f <- tern_ex_adtte %>%
+#'   filter(PARAMCD == "OS") %>%
+#'   mutate(
+#'     AVAL = day2month(AVAL),
+#'     is_event = CNSR == 0
+#'   )
+#'
+#' basic_table() %>%
+#'   split_cols_by(var = "ARMCD", ref_group = "ARM B", split_fun = ref_group_last) %>%
+#'   add_colcounts() %>%
+#'   surv_time(
+#'     vars = "AVAL",
+#'     var_labels = "Survival Time (Months)",
+#'     is_event = "is_event",
+#'   ) %>%
+#'   build_table(df = adtte_f)
+#'
+#' @export
+ref_group_last <- make_split_fun(
+  post = list(
+    function(splret, spl, fulldf) {
+      if (!"ref_group_value" %in% slotNames(spl)) {
+        stop("Reference group is undefined.")
+      }
+
+      spl_var <- spl@payload # can be substituted by splret$labels
+      init_lvls <- levels(fulldf[[spl_var]])
+
+      if (!all(names(splret$values) %in% init_lvls)) {
+        stop("This split function does not work with combination facets.")
+      }
+
+      ref_group_pos <- which(init_lvls == spl@ref_group_value)
+      reord_lvls <- c(init_lvls[-ref_group_pos], init_lvls[ref_group_pos])
+      ord <- match(reord_lvls, names(splret$values))
+
+      make_split_result(
+        splret$values[ord],
+        splret$datasplit[ord],
+        splret$labels[ord]
+      )
+    }
+  )
+)
+#' @describeIn utils_split_funs split function to keep original order of factor
+#'  levels in the split.
+#'
+#' @return
+#' * `keep_level_order` returns an utility function that keeps the original levels'.
+#'  It needs to be assigned to `split_fun`.
+#'
+#' @examples
+#' # keep_level_order --------
+#' # Even if default would bring ref_group first, the original order puts it last
+#' basic_table() %>%
+#'   split_cols_by("Species", ref_group = "virginica", split_fun = keep_level_order) %>%
+#'   analyze("Sepal.Length") %>%
+#'   build_table(iris)
+#'
+#' @export
+keep_level_order <- make_split_fun(
+  post = list(
+    function(splret, spl, fulldf, ...) {
+      # browser() if you enter here the order of splret seems already correct
+      ord <- order(names(splret$values))
+      make_split_result(splret$values[ord],
+                        splret$datasplit[ord],
+                        splret$labels[ord])
+    }
+  )
+)
diff --git a/man/ref_group_last.Rd b/man/ref_group_last.Rd
deleted file mode 100644
index 8e49244cd7..0000000000
--- a/man/ref_group_last.Rd
+++ /dev/null
@@ -1,57 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/reflast.R
-\name{ref_group_last}
-\alias{ref_group_last}
-\title{Split Function to Place Reference Group Facet Last}
-\usage{
-ref_group_last(splret, spl, fulldf)
-}
-\arguments{
-\item{splret}{result of the core split}
-
-\item{spl}{split object}
-
-\item{fulldf}{data.frame of incoming data to be split}
-}
-\description{
-Place reference group facet last during post-processing stage in a custom split function.
-}
-\examples{
-library(dplyr)
-
-# Define custom split function
-ref_last <- make_split_fun(post = list(ref_group_last))
-
-dat <- data.frame(
-  x = factor(letters[1:5], levels = letters[5:1]),
-  y = 1:5
-)
-
-# With rtables layout functions
-basic_table() \%>\%
-  split_cols_by("x", ref_group = "c", split_fun = ref_last) \%>\%
-  analyze("y") \%>\%
-  build_table(dat)
-
-# With tern layout funcitons
-adtte_f <- tern_ex_adtte \%>\%
-  filter(PARAMCD == "OS") \%>\%
-  mutate(
-    AVAL = day2month(AVAL),
-    is_event = CNSR == 0
-  )
-
-basic_table() \%>\%
-  split_cols_by(var = "ARMCD", ref_group = "ARM B", split_fun = ref_last) \%>\%
-  add_colcounts() \%>\%
-  surv_time(
-    vars = "AVAL",
-    var_labels = "Survival Time (Months)",
-    is_event = "is_event",
-  ) \%>\%
-  build_table(df = adtte_f)
-
-}
-\seealso{
-\code{\link[rtables:make_split_fun]{rtables::make_split_fun()}}
-}
diff --git a/man/utils_split_funs.Rd b/man/utils_split_funs.Rd
new file mode 100644
index 0000000000..ad50b17489
--- /dev/null
+++ b/man/utils_split_funs.Rd
@@ -0,0 +1,85 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils_split_funs.R
+\name{utils_split_funs}
+\alias{utils_split_funs}
+\alias{ref_group_last}
+\alias{keep_level_order}
+\title{Custom Split Functions}
+\usage{
+ref_group_last(df, spl, vals = NULL, labels = NULL, trim = FALSE, .spl_context)
+
+keep_level_order(
+  df,
+  spl,
+  vals = NULL,
+  labels = NULL,
+  trim = FALSE,
+  .spl_context
+)
+}
+\value{
+\itemize{
+\item \code{ref_group_last} returns an utility function that puts the reference group
+as last and needs to be assigned to \code{split_fun}.
+}
+
+\itemize{
+\item \code{keep_level_order} returns an utility function that keeps the original levels'.
+It needs to be assigned to \code{split_fun}.
+}
+}
+\description{
+Collection of useful functions that are expanding on the core list of functions
+provided by \code{rtables}. See \code{?rtables::custom_split_funs} and \code{\link[rtables:make_split_fun]{rtables::make_split_fun()}}
+for more information on how to make a custom split function. All these functions
+work with \code{\link[=split_rows_by]{split_rows_by()}} argument \code{split_fun} to modify the way the split
+happens.
+}
+\section{Functions}{
+\itemize{
+\item \code{ref_group_last()}: split function to place reference group facet last
+during post-processing stage.
+
+\item \code{keep_level_order()}: split function to keep original order of factor
+levels in the split.
+
+}}
+\examples{
+library(dplyr)
+
+dat <- data.frame(
+  x = factor(letters[1:5], levels = letters[5:1]),
+  y = 1:5
+)
+# ref_group_last
+
+# With rtables layout functions
+basic_table() \%>\%
+  split_cols_by("x", ref_group = "c", split_fun = ref_group_last) \%>\%
+  analyze("y") \%>\%
+  build_table(dat)
+
+# With tern layout funcitons
+adtte_f <- tern_ex_adtte \%>\%
+  filter(PARAMCD == "OS") \%>\%
+  mutate(
+    AVAL = day2month(AVAL),
+    is_event = CNSR == 0
+  )
+
+basic_table() \%>\%
+  split_cols_by(var = "ARMCD", ref_group = "ARM B", split_fun = ref_group_last) \%>\%
+  add_colcounts() \%>\%
+  surv_time(
+    vars = "AVAL",
+    var_labels = "Survival Time (Months)",
+    is_event = "is_event",
+  ) \%>\%
+  build_table(df = adtte_f)
+
+# keep_level_order
+
+}
+\seealso{
+\code{\link[rtables:make_split_fun]{rtables::make_split_fun()}}
+}
diff --git a/tests/testthat/test-reflast.R b/tests/testthat/test-utils_split_fun.R
similarity index 66%
rename from tests/testthat/test-reflast.R
rename to tests/testthat/test-utils_split_fun.R
index b613cc3fcf..04e95212f3 100644
--- a/tests/testthat/test-reflast.R
+++ b/tests/testthat/test-utils_split_fun.R
@@ -1,53 +1,75 @@
-split_fun <- make_split_fun(post = list(ref_group_last))
 
 testthat::test_that("ref_group_last split fun gives error when ref group is undefined", {
   lyt <- basic_table() %>%
-    split_cols_by(var = "ARMCD", split_fun = split_fun) %>%
+    split_cols_by(var = "ARMCD", split_fun = ref_group_last) %>%
     analyze("AGE")
 
   testthat::expect_error(build_table(lyt, df = tern_ex_adsl))
 })
 
-testthat::test_that("ref_group_last split fun gives error when used with combo facets", {
-  custom_splitfun <- make_split_fun(
-    post = list(
-      add_combo_facet("A_C", "Arms A+C", c("A: Drug X", "C: Combination")),
-      ref_group_last
-    )
-  )
+# testthat::test_that("ref_group_last split fun gives error when used with combo facets", {
+#   custom_splitfun <- make_split_fun(
+#     post = list(
+#       add_combo_facet("A_C", "Arms A+C", c("A: Drug X", "C: Combination")),
+#       ref_group_last
+#     )
+#   )
+#
+#   lyt <- basic_table() %>%
+#     split_cols_by("ARM", ref_group = "B: Placebo", split_fun = custom_splitfun) %>%
+#     analyze("AGE")
+#
+#   testthat::expect_error(build_table(lyt, df = tern_ex_adsl))
+# })
 
-  lyt <- basic_table() %>%
-    split_cols_by("ARM", ref_group = "B: Placebo", split_fun = custom_splitfun) %>%
-    analyze("AGE")
+testthat::test_that("analyze_vars works as expected with ref_group_last split fun", {
+  # Default behavior
+  result <- basic_table() %>%
+    split_cols_by(var = "ARMCD", ref_group = "ARM C") %>%
+    add_colcounts() %>%
+    analyze_vars(c("AGE", "STRATA2")) %>%
+    build_table(df = tern_ex_adsl)
 
-  testthat::expect_error(build_table(lyt, df = tern_ex_adsl))
-})
+  res <- testthat::expect_silent(result)
+  testthat::expect_identical(names(res), c("ARM C", "ARM A", "ARM B"))
 
-testthat::test_that("analyze_vars works as expected with ref_group_last split fun", {
+  # ref_group_last
   result <- basic_table() %>%
-    split_cols_by(var = "ARMCD", ref_group = "ARM B", split_fun = split_fun) %>%
+    split_cols_by(var = "ARMCD", ref_group = "ARM C", split_fun = ref_group_last) %>%
     add_colcounts() %>%
     analyze_vars(c("AGE", "STRATA2")) %>%
     build_table(df = tern_ex_adsl)
 
   res <- testthat::expect_silent(result)
-  testthat::expect_snapshot(res)
+  testthat::expect_snapshot(res[3:4, ])
+
+  # keep_level_order
+  result_ordered <- basic_table() %>%
+    split_cols_by(var = "ARMCD", ref_group = "ARM B", split_fun = keep_level_order) %>%
+    add_colcounts() %>%
+    analyze_vars(c("AGE", "STRATA2")) %>%
+    build_table(df = tern_ex_adsl)
+
+  res_ordered <- testthat::expect_silent(result_ordered)
+
+  testthat::expect_identical(toString(res), toString(res_ordered))
 })
 
 testthat::test_that("compare_vars works as expected with ref_group_last split fun", {
   result <- basic_table() %>%
-    split_cols_by(var = "ARMCD", ref_group = "ARM B", split_fun = split_fun) %>%
+    split_cols_by(var = "ARMCD", ref_group = "ARM B", split_fun = ref_group_last) %>%
     add_colcounts() %>%
     compare_vars("AGE") %>%
     build_table(df = tern_ex_adsl)
 
   res <- testthat::expect_silent(result)
-  testthat::expect_snapshot(res)
+  testthat::expect_identical(names(res), c("ARM A", "ARM C", "ARM B"))
+  testthat::expect_snapshot(res[1:2, ])
 })
 
 testthat::test_that("summarize_ancova works as expected with ref_group_last split fun", {
   result <- basic_table() %>%
-    split_cols_by(var = "ARMCD", ref_group = "ARM B", split_fun = split_fun) %>%
+    split_cols_by(var = "ARMCD", ref_group = "ARM B", split_fun = ref_group_last) %>%
     add_colcounts() %>%
     summarize_ancova(
       vars = "BMRKR1",
@@ -58,7 +80,8 @@ testthat::test_that("summarize_ancova works as expected with ref_group_last spli
     build_table(tern_ex_adsl)
 
   res <- testthat::expect_silent(result)
-  testthat::expect_snapshot(res)
+  testthat::expect_identical(names(res), c("ARM A", "ARM C", "ARM B"))
+  testthat::expect_snapshot(res[1:2, ])
 })
 
 testthat::test_that("binary endpoint layouts work as expected with ref_group_last split fun", {
@@ -67,7 +90,7 @@ testthat::test_that("binary endpoint layouts work as expected with ref_group_las
     dplyr::mutate(is_rsp = AVALC %in% c("CR", "PR"))
 
   result <- basic_table() %>%
-    split_cols_by(var = "ARM", ref_group = "B: Placebo", split_fun = split_fun) %>%
+    split_cols_by(var = "ARM", ref_group = "B: Placebo", split_fun = ref_group_last) %>%
     add_colcounts() %>%
     estimate_odds_ratio(vars = "is_rsp") %>%
     estimate_proportion_diff(vars = "is_rsp", table_names = "prop_diff") %>%
@@ -75,6 +98,7 @@ testthat::test_that("binary endpoint layouts work as expected with ref_group_las
     build_table(adrs_f)
 
   res <- testthat::expect_silent(result)
+  testthat::expect_identical(names(res), c("A: Drug X", "C: Combination", "B: Placebo"))
   testthat::expect_snapshot(res)
 })
 
@@ -87,7 +111,7 @@ testthat::test_that("time to event layouts works as expected with ref_group_last
     )
 
   result <- basic_table() %>%
-    split_cols_by(var = "ARMCD", ref_group = "ARM B", split_fun = split_fun) %>%
+    split_cols_by(var = "ARMCD", ref_group = "ARM B", split_fun = ref_group_last) %>%
     add_colcounts() %>%
     coxph_pairwise(
       vars = "AVAL",