diff --git a/.github/scripts/benchmark.R b/.github/scripts/benchmark.R
new file mode 100644
index 000000000..caf2962e9
--- /dev/null
+++ b/.github/scripts/benchmark.R
@@ -0,0 +1,50 @@
+suppressPackageStartupMessages({
+  library(bench)
+  library(dplyr)
+  library(pkgload)
+})
+
+# Load the local cards package
+pkgload::load_all(".")
+
+# We use the trial dataset from gtsummary to test
+library(gtsummary)
+data_big <- trial[rep(seq_len(nrow(trial)), 20), ]
+
+cat("Running benchmarks for ard_summary and ard_tabulate...\n")
+
+res_summary <- bench::mark(
+  ard_summary(
+    data_big,
+    variables = c(age, marker),
+    by = trt
+  ),
+  check = FALSE,
+  iterations = 10,
+  filter_gc = FALSE
+)
+
+res_tabulate <- bench::mark(
+  ard_tabulate(
+    data_big,
+    variables = c(grade, trt),
+    by = response
+  ),
+  check = FALSE,
+  iterations = 10,
+  filter_gc = FALSE
+)
+
+# Format the results into a markdown table
+cat("### Performance Benchmark Results\n\n", file = "bench_report.md")
+cat("These benchmarks run on 20x replicated `gtsummary::trial` dataset.\n\n", file = "bench_report.md", append = TRUE)
+
+cat("#### `ard_summary`\n", file = "bench_report.md", append = TRUE)
+knitr::kable(summary(res_summary)[, c("expression", "min", "median", "itr/sec", "mem_alloc", "gc/sec", "n_itr", "n_gc", "total_time")]) |>
+  cat(file = "bench_report.md", append = TRUE, sep = "\n")
+
+cat("\n#### `ard_tabulate`\n", file = "bench_report.md", append = TRUE)
+knitr::kable(summary(res_tabulate)[, c("expression", "min", "median", "itr/sec", "mem_alloc", "gc/sec", "n_itr", "n_gc", "total_time")]) |>
+  cat(file = "bench_report.md", append = TRUE, sep = "\n")
+
+cat("\nBenchmark report generated successfully.\n")
diff --git a/.github/workflows/benchmark.yaml b/.github/workflows/benchmark.yaml
new file mode 100644
index 000000000..0d291a024
--- /dev/null
+++ b/.github/workflows/benchmark.yaml
@@ -0,0 +1,40 @@
+name: Performance Benchmark
+on:
+  pull_request:
+    branches: [main, master]
+
+jobs:
+  benchmark:
+    if: startsWith(github.event.pull_request.title, 'perf')
+    runs-on: ubuntu-latest
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+    permissions:
+      pull-requests: write
+    steps:
+      - name: Checkout PR
+        uses: actions/checkout@v4
+
+      - name: Setup R
+        uses: r-lib/actions/setup-r@v2
+
+      - name: Install Dependencies & Cache
+        uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          dependencies: '"hard"'
+          extra-packages: |
+            any::bench
+            any::knitr
+            any::gtsummary
+            any::pkgload
+          cache: true
+          cache-version: 3
+
+      - name: Run Benchmarks
+        run: Rscript .github/scripts/benchmark.R
+
+      - name: Post Comment
+        uses: mshick/add-pr-comment@v3
+        with:
+          message-path: bench_report.md
+          message-id: 'add-benchmark-pr-comment'
diff --git a/R/ard_summary.R b/R/ard_summary.R
index 78a590b48..53ffab5ae 100644
--- a/R/ard_summary.R
+++ b/R/ard_summary.R
@@ -262,33 +262,36 @@ ard_summary.data.frame <- function(data,
                                     by, strata, data,
                                     new_col_name = "...ard_all_stats...") {
   df_nested[[new_col_name]] <-
-    map(
+    lapply(
       df_nested[["...ard_nested_data..."]],
       function(nested_data) {
-        map(
-          variables,
-          function(variable) {
-            map2(
-              statistic[[variable]], names(statistic[[variable]]),
-              function(fun, fun_name) {
-                .lst_results_as_df(
-                  x = # calculate results, and place in tibble
-                    eval_capture_conditions(
-                      getOption(
-                        "cards.calculate_stats_as_ard.eval_fun",
-                        default = expr(do.call(fun, args = list(stats::na.omit(nested_data[[variable]]))))
-                      )
-                    ),
-                  variable = variable,
-                  fun_name = fun_name,
-                  fun = fun
-                )
-              }
-            ) |>
-              unname()
-          }
-        ) |>
-          dplyr::bind_rows()
+        res_list <- unlist(
+          lapply(
+            variables,
+            function(variable) {
+              lapply(
+                names(statistic[[variable]]),
+                function(fun_name) {
+                  fun <- statistic[[variable]][[fun_name]]
+                  .lst_results_as_df(
+                    x = # calculate results, and place in tibble
+                      eval_capture_conditions(
+                        getOption(
+                          "cards.calculate_stats_as_ard.eval_fun",
+                          default = expr(do.call(fun, args = list(stats::na.omit(nested_data[[variable]]))))
+                        )
+                      ),
+                    variable = variable,
+                    fun_name = fun_name,
+                    fun = fun
+                  )
+                }
+              )
+            }
+          ),
+          recursive = FALSE
+        )
+        dplyr::bind_rows(res_list)
       }
     )
 
@@ -324,33 +327,30 @@ ard_summary.data.frame <- function(data,
   # unnesting results if needed
   if (.is_named_list(x$result, allow_df = TRUE)) {
     if (is.data.frame(x$result)) x$result <- unclass(x$result)
-    df_ard <-
-      dplyr::tibble(
-        stat_name = names(x$result),
-        result = unname(x$result),
-        warning = list(x$warning),
-        error = list(x$error)
-      )
-  }
-  # if result is not a nested list, return a single row tibble
-  else {
-    df_ard <-
-      map(x, list) |>
-      dplyr::as_tibble() |>
-      dplyr::mutate(
-        stat_name =
-        # if the function is a "cards_fn" AND the result is missing, use the provided placeholder stat names
-          case_switch(
-            is_empty(.env$x$result) && is_cards_fn(.env$fun) ~ list(get_cards_fn_stat_names(.env$fun)),
-            .default = .env$fun_name
-          )
-      ) |>
-      tidyr::unnest("stat_name")
+    stat_name <- names(x$result)
+    res_list <- unname(x$result)
+  } else {
+    # if result is not a nested list, return a single row tibble
+    stat_name <- if (is_empty(x$result) && is_cards_fn(fun)) {
+      get_cards_fn_stat_names(fun)
+    } else {
+      fun_name
+    }
+    res_list <- list(x$result)
   }
 
-  df_ard |>
-    dplyr::mutate(variable = .env$variable) |>
-    dplyr::rename(stat = "result")
+  len <- length(stat_name)
+
+  out <- data.frame(
+    stat_name = stat_name,
+    stringsAsFactors = FALSE
+  )
+  out$stat <- res_list
+  out$warning <- rep_len(list(x$warning), len)
+  out$error <- rep_len(list(x$error), len)
+  out$variable <- rep_len(variable, len)
+
+  out
 }
 
 
diff --git a/R/ard_tabulate.R b/R/ard_tabulate.R
index b359d1025..d403b676b 100644
--- a/R/ard_tabulate.R
+++ b/R/ard_tabulate.R
@@ -325,23 +325,46 @@ ard_tabulate.data.frame <- function(data,
             tab_stats = tab_stats
           )
 
-        df_result_tabulation |>
-          .nesting_rename_ard_columns(variable = variable, by = by, strata = strata) |>
-          dplyr::mutate(
-            across(any_of(c("...ard_n...", "...ard_N...", "...ard_p...", "...ard_n_cum...", "...ard_p_cum...")), as.list),
-            across(c(matches("^group[0-9]+_level$"), any_of("variable_level")), as.list)
-          ) |>
-          tidyr::pivot_longer(
-            cols = any_of(c("...ard_n...", "...ard_N...", "...ard_p...", "...ard_n_cum...", "...ard_p_cum...")),
-            names_to = "stat_name",
-            values_to = "stat"
-          ) |>
-          dplyr::mutate(
-            stat_name =
-              gsub(pattern = "^...ard_", replacement = "", x = .data$stat_name) %>%
-                gsub(pattern = "...$", replacement = "", x = .)
-          ) |>
-          dplyr::filter(.data$stat_name %in% tab_stats[["tabulation"]])
+        df_res <- .nesting_rename_ard_columns(df_result_tabulation, variable = variable, by = by, strata = strata)
+
+        # Convert grouping columns to list
+        for (col in c(grep("^group[0-9]+_level$", names(df_res), value = TRUE), intersect(names(df_res), "variable_level"))) {
+          if (!is.list(df_res[[col]])) {
+            df_res[[col]] <- as.list(df_res[[col]])
+          }
+        }
+
+        all_pivot_cols <- c("...ard_n...", "...ard_N...", "...ard_p...", "...ard_n_cum...", "...ard_p_cum...")
+        cols_in_df <- all_pivot_cols[all_pivot_cols %in% names(df_res)]
+
+        stat_names_clean <- gsub("...$", "", gsub("^...ard_", "", cols_in_df))
+
+        keep_idx <- stat_names_clean %in% tab_stats[["tabulation"]]
+        cols_to_pivot <- cols_in_df[keep_idx]
+        stat_names_clean <- stat_names_clean[keep_idx]
+
+        if (length(cols_to_pivot) == 0) {
+          return(dplyr::tibble())
+        }
+
+        n_rows <- nrow(df_res)
+        n_cols <- length(cols_to_pivot)
+
+        # We must drop all potential pivot columns that weren't selected to match pivot_longer behavior
+        fixed_cols <- setdiff(names(df_res), cols_in_df)
+        df_out <- df_res[fixed_cols]
+        df_out <- df_out[rep(seq_len(n_rows), each = n_cols), , drop = FALSE]
+
+        df_out$stat_name <- rep(stat_names_clean, times = n_rows)
+
+        stat_list <- vector("list", n_rows * n_cols)
+        for (i in seq_along(cols_to_pivot)) {
+          indices <- seq(i, length(stat_list), by = n_cols)
+          stat_list[indices] <- as.list(df_res[[cols_to_pivot[i]]])
+        }
+        df_out$stat <- stat_list
+
+        dplyr::as_tibble(df_out)
       }
     ) |>
     dplyr::bind_rows()
diff --git a/R/ard_tabulate_value.R b/R/ard_tabulate_value.R
index d0ed786f0..0293ad9ab 100644
--- a/R/ard_tabulate_value.R
+++ b/R/ard_tabulate_value.R
@@ -84,15 +84,22 @@ ard_tabulate_value.data.frame <- function(data,
     fmt_fun = fmt_fun,
     stat_label = stat_label
   ) |>
-    dplyr::filter(
-      pmap(
-        list(.data$variable, .data$variable_level),
-        function(variable, variable_level) {
-          variable_level %in% .env$value[[variable]]
+    dplyr::filter({
+      keep <- logical(dplyr::n())
+      val <- .env$value
+      for (v in names(val)) {
+        idx <- which(.data$variable == v)
+        if (length(idx) > 0) {
+          v_levels <- .data$variable_level[idx]
+          if (all(lengths(v_levels) == 1L)) {
+            keep[idx] <- unlist(v_levels, recursive = FALSE, use.names = FALSE) %in% val[[v]]
+          } else {
+            keep[idx] <- vapply(v_levels, function(x) x %in% val[[v]], logical(1))
+          }
         }
-      ) |>
-        unlist()
-    ) |>
+      }
+      keep
+    }) |>
     dplyr::mutate(context = "tabulate_value")
 
   # append attributes ----------------------------------------------------------
diff --git a/R/ard_total_n.R b/R/ard_total_n.R
index ca180be85..7ec1b3b20 100644
--- a/R/ard_total_n.R
+++ b/R/ard_total_n.R
@@ -29,13 +29,20 @@ ard_total_n.data.frame <- function(data, ...) {
   check_data_frame(data)
 
   # calculate total N ----------------------------------------------------------
-  data |>
-    dplyr::mutate(..ard_total_n.. = TRUE) |>
-    ard_tabulate_value(
-      variables = "..ard_total_n..",
-      statistic = list(..ard_total_n.. = "N"),
-      value = list(..ard_total_n.. = TRUE)
+  ans <-
+    dplyr::tibble(
+      variable = "..ard_total_n..",
+      context = "total_n",
+      stat_name = "N",
+      stat_label = "N",
+      stat = list(nrow(data)),
+      fmt_fun = list(0L),
+      warning = list(NULL),
+      error = list(NULL)
     ) |>
-    dplyr::mutate(context = "total_n") |>
-    dplyr::select(-all_ard_variables("levels"))
+    as_card()
+
+  attr(ans, "args") <- list(variables = "..ard_total_n..", by = character(0), strata = character(0))
+
+  ans
 }
diff --git a/R/tidy_ard_order.R b/R/tidy_ard_order.R
index bcdb81cf8..133167098 100644
--- a/R/tidy_ard_order.R
+++ b/R/tidy_ard_order.R
@@ -37,30 +37,38 @@ tidy_ard_column_order <- function(x, group_order = c("ascending", "descending"))
   set_cli_abort_call()
   group_order <- arg_match(group_order)
 
-  # specify the ordering the grouping variables
-  group_cols <-
-    data.frame(colname = dplyr::select(x, all_ard_groups()) |> names()) |>
-    dplyr::arrange(
-      case_switch(
-        group_order == "ascending" ~ as.integer(unlist(str_extract_all(.data$colname, "\\d+"))),
-        group_order == "descending" ~ dplyr::desc(as.integer(unlist(str_extract_all(.data$colname, "\\d+"))))
-      ),
-      .data$colname
-    ) |>
-    dplyr::pull("colname")
+  nms <- names(x)
+  group_cols <- grep("^group[0-9]+$", nms, value = TRUE)
+  group_level_cols <- grep("^group[0-9]+_level$", nms, value = TRUE)
+  all_group_cols <- c(group_cols, group_level_cols)
 
-  # selecting the columns in the tidy order
-  dplyr::select(
-    x,
-    all_of(group_cols),
-    all_ard_variables(),
-    any_of(c(
-      "context",
-      "stat_name", "stat_label", "stat", "stat_fmt", "fmt_fun",
-      "warning", "error"
-    )),
-    dplyr::everything()
-  )
+  if (length(all_group_cols) > 0) {
+    nums <- as.integer(gsub("[^0-9]", "", all_group_cols))
+    if (group_order == "ascending") {
+      o <- order(nums, all_group_cols)
+    } else {
+      o <- order(-nums, all_group_cols)
+    }
+    ordered_groups <- all_group_cols[o]
+  } else {
+    ordered_groups <- character(0)
+  }
+
+  var_cols <- grep("^variable[0-9]*$", nms, value = TRUE)
+  var_lvl_cols <- grep("^variable[0-9]*_level$", nms, value = TRUE)
+  all_var_cols <- c(var_cols, var_lvl_cols)
+
+  std_cols <- intersect(c(
+    "context",
+    "stat_name", "stat_label", "stat", "stat_fmt", "fmt_fun",
+    "warning", "error"
+  ), nms)
+
+  other_cols <- setdiff(nms, c(ordered_groups, all_var_cols, std_cols))
+
+  new_cols <- c(ordered_groups, all_var_cols, std_cols, other_cols)
+
+  x[, new_cols, drop = FALSE]
 }
 
 
@@ -69,18 +77,28 @@ tidy_ard_column_order <- function(x, group_order = c("ascending", "descending"))
 tidy_ard_row_order <- function(x) {
   set_cli_abort_call()
 
-  # get columns that dictate ordering
-  cols <- x |>
-    dplyr::select(all_ard_groups(c("names", "levels"))) |>
-    names()
-  if (!is_empty(cols)) {
-    max_group_n <- as.integer(unlist(str_extract_all(cols, "\\d+"))) |> max()
-    cols <-
-      map(seq_len(max_group_n), ~ c(paste0("group", .x), paste0("group", .x, "_level"))) |>
-      unlist() |>
-      intersect(cols)
+  nms <- names(x)
+  group_cols <- grep("^group[0-9]+$", nms, value = TRUE)
+  group_lvl_cols <- grep("^group[0-9]+_level$", nms, value = TRUE)
+  cols <- c(group_cols, group_lvl_cols)
+
+  if (length(cols) > 0) {
+    nums <- as.integer(gsub("[^0-9]", "", cols))
+    max_group_n <- max(nums)
+
+    target_cols <- character(0)
+    for (i in seq_len(max_group_n)) {
+      target_cols <- c(target_cols, paste0("group", i), paste0("group", i, "_level"))
+    }
+    cols <- intersect(target_cols, cols)
+
+    if (length(cols) > 0) {
+      order_list <- lapply(cols, function(col) match(x[[col]], unique(x[[col]])))
+      o <- do.call(order, order_list)
+      x <- x[o, , drop = FALSE]
+      rownames(x) <- NULL
+    }
   }
 
-  # perform the ordering
-  x |> dplyr::arrange(across(all_of(cols), .fns = function(x) match(x, unique(x))))
+  x
 }
diff --git a/benchmarks/bench_ard_summary.R b/benchmarks/bench_ard_summary.R
new file mode 100644
index 000000000..53f19b23d
--- /dev/null
+++ b/benchmarks/bench_ard_summary.R
@@ -0,0 +1,25 @@
+suppressPackageStartupMessages({
+  library(bench)
+  library(dplyr)
+  library(pkgload)
+})
+
+pkgload::load_all("/home/kpagacz/jnj/gtsummary-perf/cards")
+
+# We can use the trial dataset from gtsummary to test
+library(gtsummary)
+data_big <- trial[rep(seq_len(nrow(trial)), 20), ]
+
+cat("Benchmarking ard_summary...\n")
+
+res <- bench::mark(
+  ard_summary(
+    data_big,
+    variables = c(age, marker),
+    by = trt
+  ),
+  check = FALSE,
+  iterations = 10,
+  filter_gc = FALSE
+)
+print(res)
diff --git a/benchmarks/bench_ard_tabulate.R b/benchmarks/bench_ard_tabulate.R
new file mode 100644
index 000000000..b790512a1
--- /dev/null
+++ b/benchmarks/bench_ard_tabulate.R
@@ -0,0 +1,23 @@
+suppressPackageStartupMessages({
+  library(bench)
+  library(dplyr)
+  library(pkgload)
+})
+
+pkgload::load_all("/home/kpagacz/jnj/gtsummary-perf/cards")
+library(gtsummary)
+data_big <- trial[rep(seq_len(nrow(trial)), 20), ]
+
+cat("Benchmarking ard_tabulate...\n")
+
+res <- bench::mark(
+  ard_tabulate(
+    data_big,
+    variables = c(grade, trt),
+    by = response
+  ),
+  check = FALSE,
+  iterations = 10,
+  filter_gc = FALSE
+)
+print(res)
diff --git a/benchmarks/test_nest.R b/benchmarks/test_nest.R
new file mode 100644
index 000000000..b6f6780d3
--- /dev/null
+++ b/benchmarks/test_nest.R
@@ -0,0 +1,122 @@
+library(cards)
+library(dplyr)
+data_big <- gtsummary::trial[rep(seq_len(nrow(gtsummary::trial)), 2), ]
+
+res1 <- cards::nest_for_ard(
+  data_big,
+  by = "trt",
+  strata = "grade"
+)
+
+# Implementation of the optimized version
+nest_for_ard2 <- function(data, by = NULL, strata = NULL, key = "data",
+                          rename_columns = TRUE, list_columns = TRUE,
+                          include_data = TRUE,
+                          include_by_and_strata = FALSE) {
+  if (is_empty(by) && is_empty(strata)) {
+    return((dplyr::tibble("{key}" := list(data))))
+  }
+
+  n_missing <- nrow(data) - nrow(tidyr::drop_na(data, all_of(by), all_of(strata)))
+  if (n_missing > 0L) {
+    cli::cli_inform("{n_missing} missing observation{?s} in the {.val {c(by, strata)}} column{?s} have been removed.")
+  }
+
+  # create nested strata data --------------------------------------------------
+  if (!is_empty(strata)) {
+    df_strata <-
+      data[strata] |>
+      tidyr::drop_na() |>
+      dplyr::distinct() |>
+      dplyr::arrange(across(all_of(strata)))
+  }
+
+  # create nested by data --------------------------------------------------
+  if (!is_empty(by)) {
+    lst_unique_vals <-
+      by |>
+      lapply(FUN = function(x) data[[x]] |> cards:::.unique_and_sorted()) |>
+      stats::setNames(nm = by)
+
+    df_by <- tidyr::expand_grid(!!!lst_unique_vals)
+  }
+
+  # combining by and strata data sets into one, as needed ----------------------
+  if (!is_empty(by) && is_empty(strata)) {
+    df_return <- df_by
+  } else if (is_empty(by) && !is_empty(strata)) {
+    df_return <- df_strata
+  } else if (!is_empty(by) && !is_empty(strata)) {
+    df_return <-
+      df_strata |>
+      dplyr::mutate(
+        "{key}" := list(df_by),
+        .before = 0L
+      ) |>
+      tidyr::unnest(cols = all_of(key))
+  }
+
+  if (isTRUE(include_data)) {
+    # Optimized subsetting logic using inner_join
+    cols_to_keep <- if (!include_by_and_strata) {
+      setdiff(names(data), c(by, strata))
+    } else {
+      names(data)
+    }
+
+    # Left join to find group index for each row
+    df_return_idx <- df_return
+    df_return_idx[["..group_idx.."]] <- seq_len(nrow(df_return))
+
+    data_mapped <- dplyr::inner_join(
+      data |> dplyr::mutate("..row_idx.." = seq_len(nrow(data))),
+      df_return_idx[c(by, strata, "..group_idx..")],
+      by = c(by, strata)
+    )
+
+    row_indices <- split(
+      data_mapped[["..row_idx.."]],
+      factor(data_mapped[["..group_idx.."]], levels = seq_len(nrow(df_return)))
+    )
+
+    df_return[[key]] <- lapply(
+      row_indices,
+      function(idx) {
+        # drop = FALSE handles the case when there's only 1 column
+        data[idx, cols_to_keep, drop = FALSE]
+      }
+    )
+  }
+
+  # put variable levels in list to preserve types when stacked -----------------
+  if (isTRUE(list_columns)) {
+    df_return <-
+      df_return |>
+      dplyr::mutate(across(.cols = -any_of(key), .fns = as.list))
+  }
+
+  # rename by and strata columns to group## and group##_level ------------------
+  if (isTRUE(rename_columns)) {
+    df_return <-
+      df_return |>
+      cards:::.nesting_rename_ard_columns(by = by, strata = strata)
+  }
+
+  ard_final <- df_return |>
+    dplyr::as_tibble()
+
+  attr(ard_final, "args") <- list(
+    by = by,
+    strata = strata
+  )
+
+  ard_final
+}
+
+res2 <- nest_for_ard2(
+  data_big,
+  by = "trt",
+  strata = "grade"
+)
+
+cat("Identical? ", identical(res1, res2), "\n")