diff --git a/NAMESPACE b/NAMESPACE
index d5c5c5ed3..437a50106 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -24,6 +24,7 @@ S3method(process_selectors,data.frame)
 S3method(tibble::tbl_sum,card)
 export("%>%")
 export(add_calculated_row)
+export(add_hierarchical_zero_rows)
 export(alias_as_fmt_fn)
 export(alias_as_fmt_fun)
 export(all_ard_group_n)
diff --git a/NEWS.md b/NEWS.md
index e3d349f78..78c9eca8b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,9 @@
 # cards 0.9.0.9002
 
+## New Features and Functions
+
+* Added `add_hierarchical_zero_rows()` to add zero-count rows to `ard_stack_hierarchical()` and `ard_stack_hierarchical_count()` results for the level combinations that do not occur in the data, e.g. an SOC with no events. The expected combinations are passed as a data frame with columns named after the hierarchy variables. (#602, @Melkiades)
+
 ## Bug Fixes
 
 * Fixed a bug in `print.compare_ard()` where the rows present in `x` but not in `y` were announced and then not printed. (#605, @malanbos)
diff --git a/R/add_hierarchical_zero_rows.R b/R/add_hierarchical_zero_rows.R
new file mode 100644
index 000000000..a62b115ef
--- /dev/null
+++ b/R/add_hierarchical_zero_rows.R
@@ -0,0 +1,199 @@
+#' Add Zero Rows to Stacked Hierarchical ARDs
+#'
+#' @description `r lifecycle::badge('experimental')`\cr
+#'
+#' A stacked hierarchical ARD only has rows for the levels observed in the data:
+#' an SOC with no events, or an AE never reported under an SOC, has no rows.
+#'
+#' `add_hierarchical_zero_rows()` adds these rows. Pass the level
+#' combinations you expect in `levels`: each combination missing from `x` is
+#' added with `n = 0`, the denominator `N` of its `by` group, and `p = 0`. As in
+#' the rest of the ARD, `p` is `NaN` only when the `by` group is empty (`N = 0`).
+#'
+#' @param x (`card`)\cr
+#'   a stacked hierarchical ARD of class `'card'` created using [ard_stack_hierarchical()] or
+#'   [ard_stack_hierarchical_count()].
+#' @param levels (`data.frame`)\cr
+#'   the expected level combinations, one per row. The columns are named after the
+#'   hierarchy variables of `x` and start from the top of the hierarchy, e.g. `AESOC`
+#'   alone, or `AESOC` and `AEDECOD`. Their order does not matter.
+#'
+#' @details
+#' Rows are only added for the variables with summary rows in `x` (see the
+#' `include` argument of [ard_stack_hierarchical()]).
+#'
+#' The result is sorted alphanumerically, as returned by [ard_stack_hierarchical()],
+#' so that each added level sits under its parent. Use [sort_ard_hierarchical()]
+#' afterwards for a different order.
+#'
+#' @return an ARD data frame of class 'card'
+#' @seealso [sort_ard_hierarchical()]
+#' @name add_hierarchical_zero_rows
+#'
+#' @examplesIf (identical(Sys.getenv("NOT_CRAN"), "true") || identical(Sys.getenv("IN_PKGDOWN"), "true"))
+#' # AEs of two SOCs: no severe AE is reported for EYE DISORDERS,
+#' # and only some of the CARDIAC DISORDERS AEs are severe
+#' adae <- ADAE |>
+#'   dplyr::filter(AESOC %in% c("CARDIAC DISORDERS", "EYE DISORDERS"))
+#'
+#' ard <- adae |>
+#'   dplyr::filter(AESEV == "SEVERE") |>
+#'   ard_stack_hierarchical(
+#'     variables = c(AESOC, AEDECOD),
+#'     by = TRTA,
+#'     denominator = ADSL,
+#'     id = USUBJID
+#'   )
+#'
+#' # add the SOC with no severe AEs
+#' add_hierarchical_zero_rows(ard, levels = dplyr::distinct(adae, AESOC))
+#'
+#' # add every SOC/AE combination of the two SOCs
+#' add_hierarchical_zero_rows(ard, levels = dplyr::distinct(adae, AESOC, AEDECOD))
+NULL
+
+#' @rdname add_hierarchical_zero_rows
+#' @export
+add_hierarchical_zero_rows <- function(x, levels) {
+  set_cli_abort_call()
+
+  # check inputs ---------------------------------------------------------------
+  check_not_missing(x)
+  check_not_missing(levels)
+  check_class(x, c("ard_stack_hierarchical", "ard_stack_hierarchical_count"))
+  check_data_frame(levels)
+
+  ard_args <- attributes(x)$args
+  variables <- ard_args$variables
+
+  # the columns of `levels` are the top hierarchy variables, in any order
+  top_variables <- variables[seq_len(ncol(levels))]
+  if (ncol(levels) == 0L || !setequal(names(levels), top_variables)) {
+    cli::cli_abort(
+      c(
+        "The columns of {.arg levels} must be hierarchy variables of {.arg x}, starting from the top.",
+        i = "The hierarchy variables of {.arg x}, from the top, are {.val {variables}}."
+      ),
+      call = get_cli_abort_call()
+    )
+  }
+  if (anyNA(levels)) {
+    cli::cli_abort(
+      "The {.arg levels} argument cannot contain missing values.",
+      call = get_cli_abort_call()
+    )
+  }
+  levels <- dplyr::distinct(levels[top_variables])
+
+  # only the variables with summary rows in `x` get new rows, and only `n`, `N`
+  # and `p` have a known value for an unobserved level
+  complete_variables <- intersect(top_variables, ard_args$include)
+  is_summary_row <- x$context %in% c("hierarchical", "hierarchical_count") & x$variable %in% complete_variables
+  unexpected_stats <- setdiff(x$stat_name[is_summary_row], c("n", "N", "p"))
+  if (!is_empty(unexpected_stats)) {
+    cli::cli_abort(
+      c(
+        "Cannot add unobserved levels to an ARD with the {.val {unexpected_stats}} statistic{?s}.",
+        i = "Only the {.val {c('n', 'N', 'p')}} statistics can be set for an unobserved level."
+      ),
+      call = get_cli_abort_call()
+    )
+  }
+
+  # build the rows of the unobserved levels ------------------------------------
+  new_rows <-
+    lapply(
+      complete_variables,
+      \(variable) {
+        .unobserved_level_rows(
+          x = x[is_summary_row & x$variable == variable, ],
+          levels = levels,
+          variable = variable,
+          ancestors = variables[seq_len(match(variable, variables) - 1L)]
+        )
+      }
+    ) |>
+    dplyr::bind_rows()
+
+  if (nrow(new_rows) == 0L) {
+    return(x)
+  }
+
+  # sort as `ard_stack_hierarchical()` does, so each added level sits under its parent
+  dplyr::bind_rows(x, new_rows) |>
+    sort_ard_hierarchical(sort = "alphanumeric")
+}
+
+# rows for the combinations of `levels` missing from `x`, the summary rows of
+# `variable`. Each layout of `x` (one per set of group columns: with the `by`
+# variables, or overall without them) is completed from a template, the rows of
+# one observed level, which hold every `by` group with its denominator `N`
+.unobserved_level_rows <- function(x, levels, variable, ancestors) {
+  if (nrow(x) == 0L) {
+    cli::cli_abort(
+      "Cannot add levels for {.val {variable}}: {.arg x} has no rows for this variable to copy.",
+      call = get_cli_abort_call()
+    )
+  }
+  expected <- dplyr::distinct(levels[c(ancestors, variable)])
+  expected_chr <- dplyr::mutate(expected, across(everything(), as.character))
+  group_cols <- grep("^group[0-9]+$", names(x), value = TRUE)
+
+  lapply(
+    vctrs::vec_group_loc(x[group_cols])$loc,
+    \(loc) {
+      x_layout <- x[loc, ]
+      # the level columns of the ancestors and of `variable` in this layout
+      level_cols <- c(
+        sprintf("%s_level", group_cols[match(ancestors, unlist(x_layout[1L, group_cols]))]),
+        "variable_level"
+      )
+      observed <- x_layout[level_cols] |>
+        unlist_ard_columns(columns = everything()) |>
+        dplyr::mutate(across(everything(), as.character)) |>
+        stats::setNames(names(expected))
+      missing <- vctrs::vec_slice(expected, !vctrs::vec_in(expected_chr, observed))
+      if (nrow(missing) == 0L) {
+        return(NULL)
+      }
+
+      # copy the template once per missing combination, then set its levels
+      template <- x_layout[vctrs::vec_equal(observed, vctrs::vec_slice(observed, 1L)), ]
+      i_missing <- rep(seq_len(nrow(missing)), each = nrow(template))
+      new <- template[rep(seq_len(nrow(template)), times = nrow(missing)), ]
+      for (i in seq_along(level_cols)) {
+        new[[level_cols[i]]] <- .as_levels_like(
+          missing[[i]][i_missing],
+          like = template[[level_cols[i]]][[1L]],
+          arg = paste0("levels$", names(missing)[i])
+        )
+      }
+
+      # `n` is 0 and so is `p = n / N`, except in an empty `by` group (`N = 0`):
+      # there `p` is NaN, which is exactly where the template's own `p` is NaN
+      new$stat[new$stat_name == "n"] <- list(0L)
+      is_p <- new$stat_name == "p"
+      new$stat[is_p] <- lapply(new$stat[is_p], \(p) if (is.nan(p)) NaN else 0)
+      for (col in intersect(c("warning", "error"), names(new))) {
+        new[[col]] <- rep(list(NULL), nrow(new))
+      }
+      if ("stat_fmt" %in% names(new)) {
+        new <- apply_fmt_fun(new, replace = TRUE)
+      }
+      new
+    }
+  ) |>
+    dplyr::bind_rows()
+}
+
+# new level values in the type of the levels already in the ARD (`like`), as a
+# list column. A factor keeps its levels, extended by any new value
+.as_levels_like <- function(values, like, arg) {
+  values <-
+    if (is.factor(like)) {
+      factor(as.character(values), levels = union(levels(like), as.character(values)))
+    } else {
+      vctrs::vec_cast(values, like, x_arg = arg, call = get_cli_abort_call())
+    }
+  as.list(values)
+}
diff --git a/_pkgdown.yml b/_pkgdown.yml
index d7487c9bd..992455cd2 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -83,6 +83,7 @@ reference:
       - as_cards_fn
   - subtitle: "Wrangle ARD"
     contents:
+      - add_hierarchical_zero_rows
       - diff_ard_hierarchical
       - filter_ard_hierarchical
       - sort_ard_hierarchical
diff --git a/man/add_hierarchical_zero_rows.Rd b/man/add_hierarchical_zero_rows.Rd
new file mode 100644
index 000000000..9de9ffb42
--- /dev/null
+++ b/man/add_hierarchical_zero_rows.Rd
@@ -0,0 +1,66 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/add_hierarchical_zero_rows.R
+\name{add_hierarchical_zero_rows}
+\alias{add_hierarchical_zero_rows}
+\title{Add Zero Rows to Stacked Hierarchical ARDs}
+\usage{
+add_hierarchical_zero_rows(x, levels)
+}
+\arguments{
+\item{x}{(\code{card})\cr
+a stacked hierarchical ARD of class \code{'card'} created using \code{\link[=ard_stack_hierarchical]{ard_stack_hierarchical()}} or
+\code{\link[=ard_stack_hierarchical_count]{ard_stack_hierarchical_count()}}.}
+
+\item{levels}{(\code{data.frame})\cr
+the expected level combinations, one per row. The columns are named after the
+hierarchy variables of \code{x} and start from the top of the hierarchy, e.g. \code{AESOC}
+alone, or \code{AESOC} and \code{AEDECOD}. Their order does not matter.}
+}
+\value{
+an ARD data frame of class 'card'
+}
+\description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#experimental}{\figure{lifecycle-experimental.svg}{options: alt='[Experimental]'}}}{\strong{[Experimental]}}\cr
+
+A stacked hierarchical ARD only has rows for the levels observed in the data:
+an SOC with no events, or an AE never reported under an SOC, has no rows.
+
+\code{add_hierarchical_zero_rows()} adds these rows. Pass the level
+combinations you expect in \code{levels}: each combination missing from \code{x} is
+added with \code{n = 0}, the denominator \code{N} of its \code{by} group, and \code{p = 0}. As in
+the rest of the ARD, \code{p} is \code{NaN} only when the \code{by} group is empty (\code{N = 0}).
+}
+\details{
+Rows are only added for the variables with summary rows in \code{x} (see the
+\code{include} argument of \code{\link[=ard_stack_hierarchical]{ard_stack_hierarchical()}}).
+
+The result is sorted alphanumerically, as returned by \code{\link[=ard_stack_hierarchical]{ard_stack_hierarchical()}},
+so that each added level sits under its parent. Use \code{\link[=sort_ard_hierarchical]{sort_ard_hierarchical()}}
+afterwards for a different order.
+}
+\examples{
+\dontshow{if ((identical(Sys.getenv("NOT_CRAN"), "true") || identical(Sys.getenv("IN_PKGDOWN"), "true"))) withAutoprint(\{ # examplesIf}
+# AEs of two SOCs: no severe AE is reported for EYE DISORDERS,
+# and only some of the CARDIAC DISORDERS AEs are severe
+adae <- ADAE |>
+  dplyr::filter(AESOC \%in\% c("CARDIAC DISORDERS", "EYE DISORDERS"))
+
+ard <- adae |>
+  dplyr::filter(AESEV == "SEVERE") |>
+  ard_stack_hierarchical(
+    variables = c(AESOC, AEDECOD),
+    by = TRTA,
+    denominator = ADSL,
+    id = USUBJID
+  )
+
+# add the SOC with no severe AEs
+add_hierarchical_zero_rows(ard, levels = dplyr::distinct(adae, AESOC))
+
+# add every SOC/AE combination of the two SOCs
+add_hierarchical_zero_rows(ard, levels = dplyr::distinct(adae, AESOC, AEDECOD))
+\dontshow{\}) # examplesIf}
+}
+\seealso{
+\code{\link[=sort_ard_hierarchical]{sort_ard_hierarchical()}}
+}
diff --git a/tests/testthat/_snaps/add_hierarchical_zero_rows.md b/tests/testthat/_snaps/add_hierarchical_zero_rows.md
new file mode 100644
index 000000000..79c54ded9
--- /dev/null
+++ b/tests/testthat/_snaps/add_hierarchical_zero_rows.md
@@ -0,0 +1,68 @@
+# add_hierarchical_zero_rows() errors with incorrect inputs
+
+    Code
+      add_hierarchical_zero_rows(ard_tabulate(cards::ADSL, variables = AGEGR1),
+      levels = levels_all)
+    Condition
+      Error in `add_hierarchical_zero_rows()`:
+      ! The `x` argument must be class <ard_stack_hierarchical/ard_stack_hierarchical_count>, not a <card> object.
+
+---
+
+    Code
+      add_hierarchical_zero_rows(ard, levels = "EYE DISORDERS")
+    Condition
+      Error in `add_hierarchical_zero_rows()`:
+      ! The `levels` argument must be class <data.frame>, not a string.
+
+---
+
+    Code
+      add_hierarchical_zero_rows(ard, levels = data.frame(TRTA = "Placebo"))
+    Condition
+      Error in `add_hierarchical_zero_rows()`:
+      ! The columns of `levels` must be hierarchy variables of `x`, starting from the top.
+      i The hierarchy variables of `x`, from the top, are "AESOC" and "AEDECOD".
+
+---
+
+    Code
+      add_hierarchical_zero_rows(ard, levels = dplyr::distinct(levels_all, AEDECOD))
+    Condition
+      Error in `add_hierarchical_zero_rows()`:
+      ! The columns of `levels` must be hierarchy variables of `x`, starting from the top.
+      i The hierarchy variables of `x`, from the top, are "AESOC" and "AEDECOD".
+
+---
+
+    Code
+      add_hierarchical_zero_rows(ard, levels = data.frame(AESOC = NA_character_))
+    Condition
+      Error in `add_hierarchical_zero_rows()`:
+      ! The `levels` argument cannot contain missing values.
+
+---
+
+    Code
+      add_hierarchical_zero_rows(ard, levels = data.frame(AESOC = 1))
+    Condition
+      Error in `add_hierarchical_zero_rows()`:
+      ! Can't convert `levels$AESOC` <double> to <character>.
+
+---
+
+    Code
+      add_hierarchical_zero_rows(ard_cum, levels = levels_all)
+    Condition
+      Error in `add_hierarchical_zero_rows()`:
+      ! Cannot add unobserved levels to an ARD with the "n_cum" statistic.
+      i Only the "n", "N", and "p" statistics can be set for an unobserved level.
+
+---
+
+    Code
+      add_hierarchical_zero_rows(dplyr::filter(ard, variable != "AEDECOD"), levels = levels_all)
+    Condition
+      Error in `add_hierarchical_zero_rows()`:
+      ! Cannot add levels for "AEDECOD": `x` has no rows for this variable to copy.
+
diff --git a/tests/testthat/test-add_hierarchical_zero_rows.R b/tests/testthat/test-add_hierarchical_zero_rows.R
new file mode 100644
index 000000000..497ff5870
--- /dev/null
+++ b/tests/testthat/test-add_hierarchical_zero_rows.R
@@ -0,0 +1,274 @@
+skip_on_cran()
+
+ADAE_subset <- cards::ADAE |>
+  dplyr::filter(AESOC %in% c("CARDIAC DISORDERS", "EYE DISORDERS"))
+
+# two CARDIAC DISORDERS AEs, each reported in all three arms: EYE DISORDERS and
+# the other CARDIAC DISORDERS AEs are unobserved
+ADAE_observed <- ADAE_subset |>
+  dplyr::filter(AEDECOD %in% c("MYOCARDIAL INFARCTION", "SINUS BRADYCARDIA"))
+
+ard <- ard_stack_hierarchical(
+  data = ADAE_observed,
+  variables = c(AESOC, AEDECOD),
+  by = TRTA,
+  denominator = cards::ADSL,
+  id = USUBJID
+)
+
+# the 27 SOC/AE combinations of the two SOCs
+levels_all <- dplyr::distinct(ADAE_subset, AESOC, AEDECOD)
+
+# "SOC / AE" of each AE row with the `n` statistic, so once per arm
+ae_rows <- function(x) {
+  x <- x |>
+    dplyr::filter(variable == "AEDECOD", stat_name == "n") |>
+    unlist_ard_columns()
+  sort(paste(x$group2_level, x$variable_level, sep = " / "))
+}
+expected_ae_rows <- sort(rep(paste(levels_all$AESOC, levels_all$AEDECOD, sep = " / "), times = 3L))
+
+test_that("add_hierarchical_zero_rows() works", {
+  expect_silent(ard_full <- add_hierarchical_zero_rows(ard, levels = levels_all))
+
+  # each SOC/AE combination has a row per arm, under its SOC
+  expect_equal(ae_rows(ard_full), expected_ae_rows)
+
+  # every added level, SOC or AE, has n = 0, the N of its arm, and p = 0
+  ard_added <- dplyr::filter(
+    ard_full,
+    variable %in% c("AESOC", "AEDECOD"),
+    !variable_level %in% list("CARDIAC DISORDERS", "MYOCARDIAL INFARCTION", "SINUS BRADYCARDIA")
+  )
+  expect_equal(nrow(ard_added), 26L * 3L * 3L) # 1 SOC and 25 AEs, in 3 arms, with n, N and p
+  expect_identical(unique(ard_added$stat[ard_added$stat_name == "n"]), list(0L))
+  expect_identical(unique(ard_added$stat[ard_added$stat_name == "p"]), list(0))
+  expect_equal(
+    ard_added |> dplyr::filter(stat_name == "N") |> dplyr::distinct(group1_level, stat),
+    ard |> dplyr::filter(variable == "AESOC", stat_name == "N") |> dplyr::distinct(group1_level, stat)
+  )
+
+  # the rows of `x` are unchanged: removing the zero-count levels gives back `x`
+  expect_equal(filter_ard_hierarchical(ard_full, sum(n) > 0), ard)
+
+  # the result is sorted like the output of `ard_stack_hierarchical()`
+  expect_equal(sort_ard_hierarchical(ard_full, sort = "alphanumeric"), ard_full)
+
+  # the column order of `levels` does not matter
+  expect_equal(
+    add_hierarchical_zero_rows(ard, levels = levels_all[c("AEDECOD", "AESOC")]),
+    ard_full
+  )
+})
+
+test_that("add_hierarchical_zero_rows() adds the top level only", {
+  ard_soc <- add_hierarchical_zero_rows(ard, levels = dplyr::distinct(ADAE_subset, AESOC))
+
+  # EYE DISORDERS is added, a row per arm and statistic, with no AE rows under it
+  expect_equal(nrow(ard_soc), nrow(ard) + 9L)
+  expect_true(any(ard_soc$variable_level %in% list("EYE DISORDERS")))
+  expect_equal(ae_rows(ard_soc), ae_rows(ard))
+})
+
+test_that("add_hierarchical_zero_rows() returns `x` unchanged when no level is missing", {
+  expect_identical(
+    add_hierarchical_zero_rows(ard, levels = dplyr::distinct(ADAE_observed, AESOC, AEDECOD)),
+    ard
+  )
+})
+
+test_that("add_hierarchical_zero_rows() works with `overall`, `over_variables`, `attributes` and `total_n`", {
+  ard_opts <- ard_stack_hierarchical(
+    data = ADAE_observed,
+    variables = c(AESOC, AEDECOD),
+    by = TRTA,
+    denominator = cards::ADSL,
+    id = USUBJID,
+    overall = TRUE,
+    over_variables = TRUE,
+    attributes = TRUE,
+    total_n = TRUE
+  )
+  ard_full <- add_hierarchical_zero_rows(ard_opts, levels = levels_all)
+
+  # an added AE has a row per arm and an overall row, all under its SOC
+  expect_equal(
+    ard_full |>
+      dplyr::filter(variable == "AEDECOD", variable_level %in% list("GLAUCOMA"), stat_name == "n") |>
+      dplyr::select(all_ard_groups()) |>
+      unlist_ard_columns() |>
+      as.data.frame(),
+    data.frame(
+      group1 = c("TRTA", "TRTA", "TRTA", "AESOC"),
+      group1_level = c("Placebo", "Xanomeline High Dose", "Xanomeline Low Dose", "EYE DISORDERS"),
+      group2 = c("AESOC", "AESOC", "AESOC", NA),
+      group2_level = c("EYE DISORDERS", "EYE DISORDERS", "EYE DISORDERS", NA)
+    ),
+    ignore_attr = "args"
+  )
+  # with the N of each arm, and of all subjects in the overall row
+  expect_equal(
+    ard_full |>
+      dplyr::filter(variable == "AEDECOD", variable_level %in% list("GLAUCOMA"), stat_name == "N") |>
+      dplyr::pull(stat) |>
+      unlist(),
+    unname(c(table(cards::ADSL$TRTA), nrow(cards::ADSL)))
+  )
+
+  # the rows that do not depend on the levels are unchanged
+  other_rows <- \(x) dplyr::filter(x, context %in% c("attributes", "total_n") | variable == "..ard_hierarchical_overall..")
+  expect_equal(other_rows(ard_full), other_rows(ard_opts))
+})
+
+test_that("add_hierarchical_zero_rows() works with `include`", {
+  ard_ae <- ard_stack_hierarchical(
+    data = ADAE_observed,
+    variables = c(AESOC, AEDECOD),
+    by = TRTA,
+    denominator = cards::ADSL,
+    id = USUBJID,
+    include = AEDECOD
+  )
+  ard_full <- add_hierarchical_zero_rows(ard_ae, levels = levels_all)
+
+  # AE rows are added under their SOC, and no SOC summary rows are created
+  expect_equal(ae_rows(ard_full), expected_ae_rows)
+  expect_false("AESOC" %in% ard_full$variable)
+})
+
+test_that("add_hierarchical_zero_rows() works with ard_stack_hierarchical_count()", {
+  ard_count <- ard_stack_hierarchical_count(
+    data = ADAE_observed,
+    variables = c(AESOC, AEDECOD),
+    by = TRTA,
+    denominator = cards::ADSL
+  )
+  ard_full <- add_hierarchical_zero_rows(ard_count, levels = levels_all)
+
+  expect_equal(ae_rows(ard_full), expected_ae_rows)
+
+  # an added level only has the `n` statistic, set to 0
+  ard_eye <- dplyr::filter(ard_full, variable == "AESOC", variable_level %in% list("EYE DISORDERS"))
+  expect_identical(unique(ard_eye$stat_name), "n")
+  expect_identical(unique(ard_eye$stat), list(0L))
+})
+
+test_that("add_hierarchical_zero_rows() works with three hierarchy variables", {
+  ard_sev <- ard_stack_hierarchical(
+    data = ADAE_subset,
+    variables = c(AESOC, AEDECOD, AESEV),
+    by = TRTA,
+    denominator = cards::ADSL,
+    id = USUBJID
+  )
+  # no severe GLAUCOMA is reported
+  ard_full <- add_hierarchical_zero_rows(
+    ard_sev,
+    levels = data.frame(AESOC = "EYE DISORDERS", AEDECOD = "GLAUCOMA", AESEV = "SEVERE")
+  )
+
+  # only the severity is added, a row per arm and statistic, under its SOC and AE
+  expect_equal(nrow(ard_full), nrow(ard_sev) + 9L)
+  ard_severe <- dplyr::filter(
+    ard_full,
+    variable == "AESEV", variable_level %in% list("SEVERE"), group3_level %in% list("GLAUCOMA")
+  )
+  expect_equal(unique(unlist(ard_severe$group2_level)), "EYE DISORDERS")
+  expect_identical(unique(ard_severe$stat[ard_severe$stat_name == "n"]), list(0L))
+})
+
+test_that("add_hierarchical_zero_rows() keeps factor levels as factors", {
+  # the factor levels are the observed ones only, so the added levels are new
+  ard_fct <- ADAE_observed |>
+    dplyr::mutate(AESOC = factor(AESOC), AEDECOD = factor(AEDECOD)) |>
+    ard_stack_hierarchical(
+      variables = c(AESOC, AEDECOD),
+      by = TRTA,
+      denominator = cards::ADSL,
+      id = USUBJID
+    )
+  ard_full <- add_hierarchical_zero_rows(ard_fct, levels = levels_all)
+
+  # added levels are factors too, so the level columns still unlist to factors
+  expect_s3_class(unlist(ard_full$variable_level[ard_full$variable == "AESOC"]), "factor")
+  expect_s3_class(unlist(ard_full$group2_level[ard_full$variable == "AEDECOD"]), "factor")
+  expect_equal(ae_rows(ard_full), expected_ae_rows)
+})
+
+test_that("add_hierarchical_zero_rows() keeps `p` as NaN in an empty `by` group", {
+  # "Arm Z" has no subjects, so N = 0 and p = 0 / 0
+  arm_levels <- c(sort(unique(cards::ADSL$TRTA)), "Arm Z")
+  ard_z <- ard_stack_hierarchical(
+    data = dplyr::mutate(ADAE_observed, TRTA = factor(TRTA, levels = arm_levels)),
+    variables = c(AESOC, AEDECOD),
+    by = TRTA,
+    denominator = dplyr::mutate(cards::ADSL, TRTA = factor(TRTA, levels = arm_levels)),
+    id = USUBJID
+  )
+  ard_full <- add_hierarchical_zero_rows(ard_z, levels = levels_all)
+
+  ard_eye_p <- dplyr::filter(
+    ard_full,
+    variable == "AESOC", variable_level %in% list("EYE DISORDERS"), stat_name == "p"
+  )
+  expect_equal(unlist(ard_eye_p$stat), c(0, 0, 0, NaN))
+})
+
+test_that("add_hierarchical_zero_rows() formats the added levels of a formatted ARD", {
+  ard_full <- add_hierarchical_zero_rows(apply_fmt_fun(ard), levels = levels_all)
+
+  ard_eye <- dplyr::filter(ard_full, variable == "AESOC", variable_level %in% list("EYE DISORDERS"))
+  expect_equal(unlist(ard_eye$stat_fmt[ard_eye$stat_name %in% c("n", "p")]), rep(c("0", "0.0"), times = 3L))
+})
+
+test_that("add_hierarchical_zero_rows() errors with incorrect inputs", {
+  # `x` is not a stacked hierarchical ARD
+  expect_snapshot(
+    add_hierarchical_zero_rows(ard_tabulate(cards::ADSL, variables = AGEGR1), levels = levels_all),
+    error = TRUE
+  )
+
+  # `levels` is not a data frame
+  expect_snapshot(
+    add_hierarchical_zero_rows(ard, levels = "EYE DISORDERS"),
+    error = TRUE
+  )
+
+  # a column of `levels` is not a hierarchy variable
+  expect_snapshot(
+    add_hierarchical_zero_rows(ard, levels = data.frame(TRTA = "Placebo")),
+    error = TRUE
+  )
+
+  # the columns of `levels` do not start from the top of the hierarchy
+  expect_snapshot(
+    add_hierarchical_zero_rows(ard, levels = dplyr::distinct(levels_all, AEDECOD)),
+    error = TRUE
+  )
+
+  # `levels` has missing values
+  expect_snapshot(
+    add_hierarchical_zero_rows(ard, levels = data.frame(AESOC = NA_character_)),
+    error = TRUE
+  )
+
+  # a level of the wrong type
+  expect_snapshot(
+    add_hierarchical_zero_rows(ard, levels = data.frame(AESOC = 1)),
+    error = TRUE
+  )
+
+  # a statistic with no known value for an unobserved level
+  ard_cum <- ard
+  ard_cum$stat_name[ard_cum$variable == "AESOC" & ard_cum$stat_name == "N"] <- "n_cum"
+  expect_snapshot(
+    add_hierarchical_zero_rows(ard_cum, levels = levels_all),
+    error = TRUE
+  )
+
+  # `x` has no rows to copy for a variable
+  expect_snapshot(
+    add_hierarchical_zero_rows(dplyr::filter(ard, variable != "AEDECOD"), levels = levels_all),
+    error = TRUE
+  )
+})