diff --git a/DESCRIPTION b/DESCRIPTION
index d03055cc..4fa12eb6 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: pammtools
 Title: Piece-Wise Exponential Additive Mixed Modeling Tools for Survival Analysis
-Version: 0.8.0
-Date: 2026-06-17
+Version: 0.8.1
+Date: 2026-07-08
 Authors@R: c(
     person("Andreas", "Bender", , "andreas.bender@stat.uni-muenchen.de", role = c("aut", "cre"), comment=c(ORCID = "0000-0001-5628-8611")),
     person("Fabian", "Scheipl", , "fabian.scheipl@stat.uni-muenchen.de", role = c("aut"), comment = c(ORCID = "0000-0001-8172-3603")),
diff --git a/NEWS.md b/NEWS.md
index be047e50..67a4129d 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,21 @@
+# pammtools 0.8.1
+
+## Enhancements
+* The cumulative post-processing functions (`add_cumu_hazard()`,
+  `add_surv_prob()`, `add_cif()` -- for both their default and `pamm_ic`
+  methods -- and `add_trans_prob()`) now **error** when `newdata` is not grouped
+  so that the time variable is unique within each group (typically a forgotten
+  `group_by()`), instead of silently accumulating cumulative quantities across
+  distinct covariate profiles / causes / transitions and returning wrong
+  curves. This generalises the guard already used for the interval-censoring
+  pooling path (and mirrors the `stopifnot()` guard in the RMST example). All
+  of these functions accept `check_grouping = FALSE` to opt out (e.g. for
+  advanced workflows that intentionally accumulate over rows with varying
+  covariates). The guard detects mis-grouping via repeated time values within
+  a group: grids built with `make_newdata()` are always caught, but hand-built
+  grids stacking profiles with *disjoint* time grids are indistinguishable
+  from a single profile with time-varying covariates and pass undetected.
+
 # pammtools 0.8.0
 
 This release collects all changes since the last CRAN version (0.7.4),
diff --git a/R/add-functions.R b/R/add-functions.R
index 2f267e10..91b420be 100644
--- a/R/add-functions.R
+++ b/R/add-functions.R
@@ -158,7 +158,9 @@ preproc_reference <- function(reference, cnames, n_rows) {
 assert_ci_supported <- function(object, ci, ci_type) {
   if (ci && ci_type %in% c("default", "delta") && !inherits(object, "lm")) {
     stop(
-      "ci_type = \"", ci_type, "\" needs analytic (coefficient-based) ",
+      "ci_type = \"",
+      ci_type,
+      "\" needs analytic (coefficient-based) ",
       "confidence intervals, which this model does not provide. ",
       "Use ci_type = \"sim\".",
       call. = FALSE
@@ -173,7 +175,8 @@ resolve_time_var <- function(time_var, object, newdata) {
     # prefer the interval end point "tend" whenever it is available (as in
     # make_newdata grids); this also covers alternative backends that are
     # neither gam nor scam. Fall back to the "interval" factor otherwise.
-    time_var <- if (is_gam || "tend" %in% colnames(newdata)) "tend" else "interval"
+    time_var <- if (is_gam || "tend" %in% colnames(newdata)) "tend" else
+      "interval"
   } else {
     assert_string(time_var)
   }
@@ -233,10 +236,17 @@ resolve_time_var <- function(time_var, object, newdata) {
 #' @details
 #' When computing cumulative hazards or survival probabilities across groups,
 #' the input data must be grouped via \code{group_by()} prior to calling
-#' \code{add_cumu_hazard()} or \code{add_surv_prob()}. Omitting
-#' \code{group_by()} will not produce an error or warning but will return
-#' silently incorrect results, as the cumulative hazard will be accumulated
-#' over the entire dataset rather than within each group.
+#' \code{add_cumu_hazard()} or \code{add_surv_prob()}, so that the cumulative
+#' quantity is accumulated within each covariate profile rather than across the
+#' whole dataset. If \code{newdata} still contains several profiles per group
+#' (i.e.\ repeated \code{time_var} values within a group, typically a forgotten
+#' \code{group_by()}), the functions now \strong{stop with an error} rather than
+#' returning silently incorrect results. Set \code{check_grouping = FALSE} to
+#' skip this safeguard. Note the check detects the \emph{common} mis-grouping
+#' cases -- in particular any grid built with \code{\link{make_newdata}}, where
+#' profiles share time values -- but cannot catch hand-built grids that stack
+#' profiles with disjoint time grids, as these are indistinguishable from a
+#' single profile with time-varying covariates.
 #' See the \href{https://adibender.github.io/pammtools/articles/convenience.html#cumulative-hazard}{workflow vignette}
 #' for a worked example.
 #'
@@ -292,20 +302,33 @@ add_hazard.default <- function(
   # make_X/get_coefs/get_Vp triplet.
   if (!is.null(reference) || type == "link") {
     if (!inherits(object, c("glm", "scam"))) {
-      stop("`reference` and `type = \"link\"` require a coefficient-based model.")
+      stop(
+        "`reference` and `type = \"link\"` require a coefficient-based model."
+      )
     }
     return(hazard_ci(
-      object, newdata,
-      reference = reference, ci = ci, type = type, ci_type = ci_type,
-      time_var = time_var, se_mult = se_mult, ...
+      object,
+      newdata,
+      reference = reference,
+      ci = ci,
+      type = type,
+      ci_type = ci_type,
+      time_var = time_var,
+      se_mult = se_mult,
+      ...
     ))
   }
   if (ci && ci_type %in% c("default", "delta")) {
     assert_ci_supported(object, ci, ci_type)
     return(hazard_ci(
-      object, newdata,
-      ci = ci, type = type, ci_type = ci_type,
-      time_var = time_var, se_mult = se_mult, ...
+      object,
+      newdata,
+      ci = ci,
+      type = type,
+      ci_type = ci_type,
+      time_var = time_var,
+      se_mult = se_mult,
+      ...
     ))
   }
 
@@ -366,6 +389,54 @@ hazard_ci <- function(
 }
 
 
+#' Guard against silently-wrong cumulative results from mis-grouped `newdata`
+#'
+#' The cumulative post-processing functions (`add_cumu_hazard`,
+#' `add_surv_prob`, `add_cif`, `add_trans_prob`) accumulate hazards over the
+#' rows of each group of `newdata`. If several covariate profiles (or causes /
+#' transitions) share a group -- typically a forgotten `group_by()` -- the
+#' accumulation silently runs across them and returns wrong curves. A correctly
+#' grouped prediction grid has a unique `time_var` value per group, so a
+#' duplicate within a group signals this situation. (The converse does not
+#' hold: profiles stacked with *disjoint* time grids have unique `time_var`
+#' values and pass undetected -- such input is indistinguishable from a
+#' legitimate single profile with time-varying covariates, so it cannot be
+#' guarded against. Grids built with `make_newdata()` always share time values
+#' across profiles and are therefore always caught.) This generalises
+#' the `stopifnot()` guard used for the RMST helper to the whole `add_*`
+#' family, turning a silent numerical error into an explicit one.
+#' @keywords internal
+stop_if_undergrouped_for_cumulation <- function(
+  newdata,
+  time_var,
+  fun = "add_*"
+) {
+  if (is.null(time_var) || !time_var %in% names(newdata)) {
+    return(invisible(newdata))
+  }
+  by_group <- split(newdata[[time_var]], group_indices(newdata))
+  if (any(vapply(by_group, anyDuplicated, integer(1)) > 0L)) {
+    stop(
+      "`",
+      fun,
+      "()` received `newdata` with repeated '",
+      time_var,
+      "' values within a group, so distinct covariate profiles (or ",
+      "causes/transitions) would be accumulated together, producing silently ",
+      "incorrect cumulative estimates. Group `newdata` so that '",
+      time_var,
+      "' is unique within each group before calling `",
+      fun,
+      "()`, e.g. ",
+      "`newdata |> group_by(<profile columns>) |> ",
+      fun,
+      "(...)`.",
+      call. = FALSE
+    )
+  }
+  invisible(newdata)
+}
+
 #' @rdname add_hazard
 #' @export
 add_cumu_hazard <- function(newdata, object, ...) {
@@ -381,6 +452,12 @@ add_cumu_hazard <- function(newdata, object, ...) {
 #'   so that cumulative hazards start at the natural origin (consistent with
 #'   \code{\link{add_surv_prob}}, \code{\link{add_cif}} and
 #'   \code{\link{add_trans_prob}}).
+#' @param check_grouping Logical. If \code{TRUE} (default), the function checks
+#'   that \code{newdata} is grouped so that the time variable is unique within
+#'   each group and \code{\link[base]{stop}}s otherwise, guarding against
+#'   silently accumulating cumulative quantities across distinct covariate
+#'   profiles (a forgotten \code{group_by()}). Set to \code{FALSE} to skip the
+#'   check (e.g.\ for internal calls on already-validated data).
 #' @importFrom dplyr bind_cols
 #' @seealso \code{\link[mgcv]{predict.gam}},
 #' \code{\link[pammtools]{add_surv_prob}}
@@ -394,10 +471,19 @@ add_cumu_hazard.default <- function(
   time_var = NULL,
   interval_length = "intlen",
   boundary = TRUE,
+  check_grouping = TRUE,
   ...
 ) {
   interval_length <- quo_name(enquo(interval_length))
 
+  if (check_grouping) {
+    stop_if_undergrouped_for_cumulation(
+      newdata,
+      resolve_time_var(time_var, object, newdata),
+      "add_cumu_hazard"
+    )
+  }
+
   if (!overwrite) {
     if ("cumu_hazard" %in% names(newdata)) {
       stop(
@@ -505,9 +591,14 @@ get_cumu_hazard <- function(
     # analytic CI on the hazard, then propagate to the cumulative hazard
     vars_exclude <- c(vars_exclude, "se", "ci_lower", "ci_upper")
     newdata <- hazard_ci(
-      object, newdata,
-      type = "response", ci = TRUE, ci_type = ci_type,
-      time_var = time_var, se_mult = se_mult, ...
+      object,
+      newdata,
+      type = "response",
+      ci = TRUE,
+      ci_type = ci_type,
+      time_var = time_var,
+      se_mult = se_mult,
+      ...
     )
     if (ci_type == "default") {
       mutate_args <- mutate_args %>%
@@ -560,10 +651,17 @@ get_cumu_hazard <- function(
 #' @details
 #' When computing cumulative hazards or survival probabilities across groups,
 #' the input data must be grouped via \code{group_by()} prior to calling
-#' \code{add_cumu_hazard()} or \code{add_surv_prob()}. Omitting
-#' \code{group_by()} will not produce an error or warning but will return
-#' silently incorrect results, as the cumulative hazard will be accumulated
-#' over the entire dataset rather than within each group.
+#' \code{add_cumu_hazard()} or \code{add_surv_prob()}, so that the cumulative
+#' quantity is accumulated within each covariate profile rather than across the
+#' whole dataset. If \code{newdata} still contains several profiles per group
+#' (i.e.\ repeated \code{time_var} values within a group, typically a forgotten
+#' \code{group_by()}), the functions now \strong{stop with an error} rather than
+#' returning silently incorrect results. Set \code{check_grouping = FALSE} to
+#' skip this safeguard. Note the check detects the \emph{common} mis-grouping
+#' cases -- in particular any grid built with \code{\link{make_newdata}}, where
+#' profiles share time values -- but cannot catch hand-built grids that stack
+#' profiles with disjoint time grids, as these are indistinguishable from a
+#' single profile with time-varying covariates.
 #' See the \href{https://adibender.github.io/pammtools/articles/convenience.html#cumulative-hazard}{workflow vignette}
 #' for a worked example.
 #'
@@ -594,10 +692,14 @@ add_surv_prob.default <- function(
   time_var = NULL,
   interval_length = "intlen",
   boundary = TRUE,
+  check_grouping = TRUE,
   ...
 ) {
   interval_length <- quo_name(enquo(interval_length))
   time_var <- resolve_time_var(time_var, object, newdata)
+  if (check_grouping) {
+    stop_if_undergrouped_for_cumulation(newdata, time_var, "add_surv_prob")
+  }
   # The boundary is a continuous-time row at time == 0 (survival 1); see
   # add_cumu_hazard(). Only added for models predicted on the continuous time
   # axis (gam/scam/pamm), not for interval-factor models (glm/PEM).
@@ -703,9 +805,14 @@ get_surv_prob <- function(
   if (ci && ci_type %in% c("default", "delta")) {
     vars_exclude <- c(vars_exclude, "se", "ci_lower", "ci_upper")
     newdata <- hazard_ci(
-      object, newdata,
-      type = "response", ci = TRUE, ci_type = ci_type,
-      time_var = time_var, se_mult = se_mult, ...
+      object,
+      newdata,
+      type = "response",
+      ci = TRUE,
+      ci_type = ci_type,
+      time_var = time_var,
+      se_mult = se_mult,
+      ...
     )
     if (ci_type == "default") {
       mutate_args <- mutate_args %>%
@@ -966,7 +1073,7 @@ get_sim_ci <- function(
 ) {
   H <- sim_hazard(object, newdata, nsim, ...)
 
-  newdata$se       <- apply(H, 1, sd)
+  newdata$se <- apply(H, 1, sd)
   newdata$ci_lower <- row_quantile(H, alpha / 2)
   newdata$ci_upper <- row_quantile(H, 1 - alpha / 2)
 
@@ -987,7 +1094,7 @@ sim_cumulative_ci <- function(
   ...
 ) {
   intlen <- newdata[[interval_length]]
-  grp    <- dplyr::group_indices(newdata)
+  grp <- dplyr::group_indices(newdata)
 
   H <- sim_hazard(object, newdata, nsim, ...)
   draws <- matrix(
@@ -1003,17 +1110,21 @@ sim_cumulative_ci <- function(
 
 get_sim_ci_cumu <- function(newdata, object, ...) {
   sim_cumulative_ci(
-    newdata, object,
+    newdata,
+    object,
     accumulate = function(w, g) ave(w, g, FUN = cumsum),
-    names = c("cumu_lower", "cumu_upper"), ...
+    names = c("cumu_lower", "cumu_upper"),
+    ...
   )
 }
 
 get_sim_ci_surv <- function(newdata, object, ...) {
   sim_cumulative_ci(
-    newdata, object,
+    newdata,
+    object,
     accumulate = function(w, g) exp(-ave(w, g, FUN = cumsum)),
-    names = c("surv_lower", "surv_upper"), ...
+    names = c("surv_lower", "surv_upper"),
+    ...
   )
 }
 
@@ -1031,13 +1142,22 @@ get_sim_ci_surv <- function(newdata, object, ...) {
 #' @param cause_var Character. Column name of the 'cause' variable.
 #' @param interval_length \code{Character}, defaults to \code{"intlen"}.
 #'   contains the interval length in `newdata`.
+#' @param check_grouping Logical. If \code{TRUE} (default), \code{stop} if
+#'   \code{newdata} is not grouped so that the time variable is unique within
+#'   each group, guarding against silently accumulating the cumulative
+#'   incidence across distinct covariate profiles or causes. Note that
+#'   \code{check_grouping = FALSE} only skips this profile-level safeguard;
+#'   the independent check that \code{newdata} is grouped by \code{cause}
+#'   (inside \code{get_cif()}) still applies.
 #'
 #' @details
 #' When computing cumulative incidence for multiple groups, the input data must
-#' be grouped via \code{group_by()} before calling this function. Omitting
-#' \code{group_by()} will not produce an error or warning but will return
-#' silently incorrect results, as the cumulative incidence will be accumulated
-#' over the entire dataset rather than within each group.
+#' be grouped via \code{group_by()} (by cause and any covariates) before calling
+#' this function. If \code{newdata} still contains several profiles per group
+#' (repeated \code{time_var} values within a group, typically a forgotten
+#' \code{group_by()}), the function now \strong{stops with an error} rather than
+#' returning silently incorrect results, as the cumulative incidence would
+#' otherwise be accumulated across profiles rather than within each group.
 #'
 #' The returned data contains one boundary row per group at \code{time_var = 0}
 #' for plotting cumulative incidence from the time origin. On this row,
@@ -1089,11 +1209,16 @@ add_cif.default <- function(
   cause_var = "cause",
   time_var = NULL,
   interval_length = "intlen",
+  check_grouping = TRUE,
   ...
 ) {
   interval_length <- quo_name(enquo(interval_length))
   time_var <- resolve_time_var(time_var, object, newdata)
 
+  if (check_grouping) {
+    stop_if_undergrouped_for_cumulation(newdata, time_var, "add_cif")
+  }
+
   if (!overwrite) {
     if ("cif" %in% names(newdata)) {
       stop(
@@ -1176,13 +1301,13 @@ get_cif.default <- function(
   assert_choice(interval_length, colnames(newdata))
 
   newdata <- arrange(newdata, .data[[time_var]], .by_group = TRUE)
-  causes     <- levels(newdata[[cause_var]])
+  causes <- levels(newdata[[cause_var]])
   cause_data <- as.character(unique(newdata[[cause_var]]))
   if (length(cause_data) > 1) {
     stop("Did you forget to group by cause?")
   }
   dt <- newdata[[interval_length]]
-  n  <- nrow(newdata)
+  n <- nrow(newdata)
 
   # Stack the rows once per cause (same covariates/time, only `cause` changes) so
   # that a single hazard prediction supplies all cause-specific hazards -- and,
@@ -1199,8 +1324,8 @@ get_cif.default <- function(
   #   dCIF_k = (h_k / h.) * S(t-) * (1 - exp(-h. * dt)),   S via all-cause hazard.
   cif_from_hazards <- function(haz_by_cause) {
     total <- Reduce(`+`, haz_by_cause)
-    surv  <- c(1, utils::head(exp(-cumsum(total * dt)), -1))
-    incr  <- ifelse(
+    surv <- c(1, utils::head(exp(-cumsum(total * dt)), -1))
+    incr <- ifelse(
       total > 0,
       (haz_by_cause[[cause_data]] / total) * surv * (1 - exp(-total * dt)),
       0
@@ -1213,8 +1338,13 @@ get_cif.default <- function(
   split_by_cause <- function(h) lapply(idx, function(ix) h[ix])
 
   # point: plug-in at the point hazards
-  newdata[["cif"]] <- pmin(pmax(
-    cif_from_hazards(split_by_cause(get_hazard(object, stacked))), 0), 1)
+  newdata[["cif"]] <- pmin(
+    pmax(
+      cif_from_hazards(split_by_cause(get_hazard(object, stacked))),
+      0
+    ),
+    1
+  )
 
   if (ci) {
     H <- sim_hazard(object, stacked, nsim)
@@ -1451,14 +1581,26 @@ transition_state_table <- function(transitions, sep = "->") {
 #'   contains the interval length in `newdata`.
 #' @param transition \code{Character}, defaults to \code{"transition"}.
 #'   contains the transition labels in `newdata`.
+#' @param check_grouping Logical. If \code{TRUE} (default), the function checks
+#'   that \code{newdata} is grouped so that the time variable is unique within
+#'   each group once \code{transition} is part of the grouping, and
+#'   \code{\link[base]{stop}}s otherwise, guarding against silently accumulating
+#'   transition probabilities across distinct covariate profiles (a forgotten
+#'   \code{group_by()}). Set to \code{FALSE} to skip the check. As for the
+#'   other cumulative \code{add_*} functions, hand-built grids stacking
+#'   profiles with disjoint time grids cannot be detected (see
+#'   \code{\link{add_cumu_hazard}}).
 #' @param ... Further arguments passed to underlying methods.
 #'
 #' @details
 #' When computing transition probabilities for multiple groups, the input data must
-#' be grouped via \code{group_by()} before calling this function. Omitting
-#' \code{group_by()} will not produce an error or warning but will return
-#' silently incorrect results, as the transition probability will be accumulated
-#' over the entire dataset rather than within each group.
+#' be grouped via \code{group_by()} before calling this function. If \code{newdata}
+#' still contains several covariate profiles per (group, transition) -- i.e.\
+#' repeated \code{time_var} values within a group once \code{transition} is added
+#' to the grouping, typically a forgotten \code{group_by()} -- the function now
+#' \strong{stops with an error} rather than returning silently incorrect results,
+#' as the transition probability would otherwise be accumulated across profiles
+#' rather than within each group.
 #'
 #' The returned data contains one boundary row per group and transition at
 #' \code{time_var = 0} for plotting transition probabilities from the time
@@ -1493,6 +1635,7 @@ add_trans_prob <- function(
   time_var = "tend",
   interval_length = "intlen",
   transition = "transition",
+  check_grouping = TRUE,
   ...
 ) {
   orig_names <- names(newdata)
@@ -1539,6 +1682,14 @@ add_trans_prob <- function(
     old_groups <- group_vars(newdata)
   }
 
+  # Transition probabilities accumulate per (covariate profile x transition).
+  # Once `transition` is part of the grouping, `time_var` must be unique within
+  # each group; a duplicate means covariate profiles were not separated (a
+  # forgotten `group_by()`), which would silently mix them.
+  if (check_grouping) {
+    stop_if_undergrouped_for_cumulation(newdata, time_var, "add_trans_prob")
+  }
+
   split_groups <- setdiff(old_groups, transition_var)
   ordering_vars <- c(split_groups, transition_var, time_var)
   newdata <- newdata[
@@ -1556,7 +1707,8 @@ add_trans_prob <- function(
         ci = FALSE,
         time_var = time_var,
         interval_length = interval_length,
-        boundary = FALSE
+        boundary = FALSE,
+        check_grouping = FALSE
       )
     )
   }
diff --git a/R/cumulative-coefficient.R b/R/cumulative-coefficient.R
index bf160c1c..48f2f1a8 100644
--- a/R/cumulative-coefficient.R
+++ b/R/cumulative-coefficient.R
@@ -244,7 +244,8 @@ get_cumu_coef_baseline <- function(
       model,
       time_var = time_var,
       interval_length = interval_length,
-      boundary = FALSE
+      boundary = FALSE,
+      check_grouping = FALSE
     ) %>%
     mutate(
       method = class(model)[1],
diff --git a/R/pool-ic.R b/R/pool-ic.R
index 235791bf..e2c3ec48 100644
--- a/R/pool-ic.R
+++ b/R/pool-ic.R
@@ -24,6 +24,9 @@ NULL
 
 # Average a per-fit point estimate (the value column produced by an `add_*`
 # default with ci = FALSE) across the imputation fits -> MI point estimate.
+# Cumulative adders pass check_grouping = FALSE through `...`: the grouping
+# guard already ran on the pooled method's input, so the per-fit calls skip
+# the (redundant) re-check.
 pooled_point <- function(object, newdata, adder, value_col, ...) {
   preds <- lapply(object[["fits"]], function(f) {
     adder(newdata, f, ci = FALSE, ...)[[value_col]]
@@ -31,29 +34,6 @@ pooled_point <- function(object, newdata, adder, value_col, ...) {
   rowMeans(do.call(cbind, preds))
 }
 
-# Pooled cumulative quantities (cumulative hazard, survival, CIF) are
-# accumulated *within* each group of `newdata`. If a group still contains
-# repeated time points, distinct covariate (or cause) combinations are being
-# pooled together -- almost always a forgotten group_by() -- which silently
-# produces wrong results. Warn so the problem is visible (the default add_*
-# methods share the same group_by() requirement).
-warn_if_pooled_undergrouped <- function(newdata, time_var) {
-  if (!time_var %in% names(newdata)) {
-    return(invisible(NULL))
-  }
-  by_group <- split(newdata[[time_var]], group_indices(newdata))
-  if (any(vapply(by_group, anyDuplicated, integer(1)) > 0L)) {
-    warning(
-      "`newdata` has repeated '", time_var, "' values within a group, so ",
-      "distinct covariate/cause combinations are pooled together when ",
-      "accumulating cumulative quantities. Separate them with group_by() to ",
-      "avoid incorrect results.",
-      call. = FALSE
-    )
-  }
-  invisible(NULL)
-}
-
 ic_prediction_grid <- function(object, newdata, time_var, interval_length) {
   fit1 <- object[["fits"]][[1]]
   tv <- resolve_time_var(time_var, fit1, newdata)
@@ -202,10 +182,13 @@ add_cumu_hazard.pamm_ic <- function(
   nsim = 500L,
   time_var = NULL,
   interval_length = "intlen",
+  check_grouping = TRUE,
   ...
 ) {
   time_var <- resolve_time_var(time_var, object[["fits"]][[1]], newdata)
-  warn_if_pooled_undergrouped(newdata, time_var)
+  if (check_grouping) {
+    stop_if_undergrouped_for_cumulation(newdata, time_var, "add_cumu_hazard")
+  }
   newdata <- drop_cumulative_boundary(newdata, time_var)
   grid <- ic_prediction_grid(object, newdata, time_var, interval_length)
   joindata <- grid[["data"]]
@@ -218,6 +201,7 @@ add_cumu_hazard.pamm_ic <- function(
     time_var = time_var,
     interval_length = interval_length,
     boundary = FALSE,
+    check_grouping = FALSE,
     ...
   )
   if (ci) {
@@ -256,10 +240,13 @@ add_surv_prob.pamm_ic <- function(
   nsim = 500L,
   time_var = NULL,
   interval_length = "intlen",
+  check_grouping = TRUE,
   ...
 ) {
   time_var <- resolve_time_var(time_var, object[["fits"]][[1]], newdata)
-  warn_if_pooled_undergrouped(newdata, time_var)
+  if (check_grouping) {
+    stop_if_undergrouped_for_cumulation(newdata, time_var, "add_surv_prob")
+  }
   newdata <- drop_cumulative_boundary(newdata, time_var)
   grid <- ic_prediction_grid(object, newdata, time_var, interval_length)
   joindata <- grid[["data"]]
@@ -272,6 +259,7 @@ add_surv_prob.pamm_ic <- function(
     time_var = time_var,
     interval_length = interval_length,
     boundary = FALSE,
+    check_grouping = FALSE,
     ...
   )
   if (ci) {
@@ -412,6 +400,7 @@ add_cif.pamm_ic <- function(
   cause_var = "cause",
   time_var = NULL,
   interval_length = "intlen",
+  check_grouping = TRUE,
   ...
 ) {
   if (!identical(object[["type"]], "cr")) {
@@ -424,7 +413,9 @@ add_cif.pamm_ic <- function(
   m <- length(object[["fits"]])
   per <- ceiling(nsim / m)
   time_var <- resolve_time_var(time_var, fit1, newdata)
-  warn_if_pooled_undergrouped(newdata, time_var)
+  if (check_grouping) {
+    stop_if_undergrouped_for_cumulation(newdata, time_var, "add_cif")
+  }
   newdata <- drop_cumulative_boundary(newdata, time_var)
   joindata <- reconstruct_cutpoints(newdata, fit1, time_var, interval_length)
 
diff --git a/man/add_cif.Rd b/man/add_cif.Rd
index 260148ef..d8a92db6 100644
--- a/man/add_cif.Rd
+++ b/man/add_cif.Rd
@@ -18,6 +18,7 @@ add_cif(newdata, object, ...)
   cause_var = "cause",
   time_var = NULL,
   interval_length = "intlen",
+  check_grouping = TRUE,
   ...
 )
 
@@ -30,6 +31,7 @@ add_cif(newdata, object, ...)
   cause_var = "cause",
   time_var = NULL,
   interval_length = "intlen",
+  check_grouping = TRUE,
   ...
 )
 }
@@ -64,16 +66,26 @@ to \code{"tend"}.}
 
 \item{interval_length}{\code{Character}, defaults to \code{"intlen"}.
 contains the interval length in \code{newdata}.}
+
+\item{check_grouping}{Logical. If \code{TRUE} (default), \code{stop} if
+\code{newdata} is not grouped so that the time variable is unique within
+each group, guarding against silently accumulating the cumulative
+incidence across distinct covariate profiles or causes. Note that
+\code{check_grouping = FALSE} only skips this profile-level safeguard;
+the independent check that \code{newdata} is grouped by \code{cause}
+(inside \code{get_cif()}) still applies.}
 }
 \description{
 Add cumulative incidence function to data
 }
 \details{
 When computing cumulative incidence for multiple groups, the input data must
-be grouped via \code{group_by()} before calling this function. Omitting
-\code{group_by()} will not produce an error or warning but will return
-silently incorrect results, as the cumulative incidence will be accumulated
-over the entire dataset rather than within each group.
+be grouped via \code{group_by()} (by cause and any covariates) before calling
+this function. If \code{newdata} still contains several profiles per group
+(repeated \code{time_var} values within a group, typically a forgotten
+\code{group_by()}), the function now \strong{stops with an error} rather than
+returning silently incorrect results, as the cumulative incidence would
+otherwise be accumulated across profiles rather than within each group.
 
 The returned data contains one boundary row per group at \code{time_var = 0}
 for plotting cumulative incidence from the time origin. On this row,
diff --git a/man/add_hazard.Rd b/man/add_hazard.Rd
index 7da081f3..9a37b9e3 100644
--- a/man/add_hazard.Rd
+++ b/man/add_hazard.Rd
@@ -37,6 +37,7 @@ add_cumu_hazard(newdata, object, ...)
   time_var = NULL,
   interval_length = "intlen",
   boundary = TRUE,
+  check_grouping = TRUE,
   ...
 )
 
@@ -58,6 +59,7 @@ add_cumu_hazard(newdata, object, ...)
   nsim = 500L,
   time_var = NULL,
   interval_length = "intlen",
+  check_grouping = TRUE,
   ...
 )
 }
@@ -123,6 +125,13 @@ Can be either bare unquoted variable name or character. Defaults to \code{"intle
 so that cumulative hazards start at the natural origin (consistent with
 \code{\link{add_surv_prob}}, \code{\link{add_cif}} and
 \code{\link{add_trans_prob}}).}
+
+\item{check_grouping}{Logical. If \code{TRUE} (default), the function checks
+that \code{newdata} is grouped so that the time variable is unique within
+each group and \code{\link[base]{stop}}s otherwise, guarding against
+silently accumulating cumulative quantities across distinct covariate
+profiles (a forgotten \code{group_by()}). Set to \code{FALSE} to skip the
+check (e.g.\ for internal calls on already-validated data).}
 }
 \description{
 Add (cumulative) hazard based on the provided data set and model.
@@ -143,10 +152,17 @@ the standard errors reported by \code{scam} itself.
 \details{
 When computing cumulative hazards or survival probabilities across groups,
 the input data must be grouped via \code{group_by()} prior to calling
-\code{add_cumu_hazard()} or \code{add_surv_prob()}. Omitting
-\code{group_by()} will not produce an error or warning but will return
-silently incorrect results, as the cumulative hazard will be accumulated
-over the entire dataset rather than within each group.
+\code{add_cumu_hazard()} or \code{add_surv_prob()}, so that the cumulative
+quantity is accumulated within each covariate profile rather than across the
+whole dataset. If \code{newdata} still contains several profiles per group
+(i.e.\ repeated \code{time_var} values within a group, typically a forgotten
+\code{group_by()}), the functions now \strong{stop with an error} rather than
+returning silently incorrect results. Set \code{check_grouping = FALSE} to
+skip this safeguard. Note the check detects the \emph{common} mis-grouping
+cases -- in particular any grid built with \code{\link{make_newdata}}, where
+profiles share time values -- but cannot catch hand-built grids that stack
+profiles with disjoint time grids, as these are indistinguishable from a
+single profile with time-varying covariates.
 See the \href{https://adibender.github.io/pammtools/articles/convenience.html#cumulative-hazard}{workflow vignette}
 for a worked example.
 }
diff --git a/man/add_surv_prob.Rd b/man/add_surv_prob.Rd
index e9e737e3..a4a1b4d7 100644
--- a/man/add_surv_prob.Rd
+++ b/man/add_surv_prob.Rd
@@ -17,6 +17,7 @@ add_surv_prob(newdata, object, ...)
   time_var = NULL,
   interval_length = "intlen",
   boundary = TRUE,
+  check_grouping = TRUE,
   ...
 )
 
@@ -28,6 +29,7 @@ add_surv_prob(newdata, object, ...)
   nsim = 500L,
   time_var = NULL,
   interval_length = "intlen",
+  check_grouping = TRUE,
   ...
 )
 }
@@ -66,6 +68,13 @@ so that cumulative hazards start at the natural origin (consistent with
 \code{\link{add_surv_prob}}, \code{\link{add_cif}} and
 \code{\link{add_trans_prob}}).}
 
+\item{check_grouping}{Logical. If \code{TRUE} (default), the function checks
+that \code{newdata} is grouped so that the time variable is unique within
+each group and \code{\link[base]{stop}}s otherwise, guarding against
+silently accumulating cumulative quantities across distinct covariate
+profiles (a forgotten \code{group_by()}). Set to \code{FALSE} to skip the
+check (e.g.\ for internal calls on already-validated data).}
+
 \item{alpha}{Significance level for pooled confidence intervals.}
 
 \item{nsim}{Total number of pooled posterior draws used for the interval.}
@@ -79,10 +88,17 @@ for the specified covariate and follow-up information (and CIs
 \details{
 When computing cumulative hazards or survival probabilities across groups,
 the input data must be grouped via \code{group_by()} prior to calling
-\code{add_cumu_hazard()} or \code{add_surv_prob()}. Omitting
-\code{group_by()} will not produce an error or warning but will return
-silently incorrect results, as the cumulative hazard will be accumulated
-over the entire dataset rather than within each group.
+\code{add_cumu_hazard()} or \code{add_surv_prob()}, so that the cumulative
+quantity is accumulated within each covariate profile rather than across the
+whole dataset. If \code{newdata} still contains several profiles per group
+(i.e.\ repeated \code{time_var} values within a group, typically a forgotten
+\code{group_by()}), the functions now \strong{stop with an error} rather than
+returning silently incorrect results. Set \code{check_grouping = FALSE} to
+skip this safeguard. Note the check detects the \emph{common} mis-grouping
+cases -- in particular any grid built with \code{\link{make_newdata}}, where
+profiles share time values -- but cannot catch hand-built grids that stack
+profiles with disjoint time grids, as these are indistinguishable from a
+single profile with time-varying covariates.
 See the \href{https://adibender.github.io/pammtools/articles/convenience.html#cumulative-hazard}{workflow vignette}
 for a worked example.
 
diff --git a/man/add_trans_prob.Rd b/man/add_trans_prob.Rd
index 68604584..1e11367a 100644
--- a/man/add_trans_prob.Rd
+++ b/man/add_trans_prob.Rd
@@ -14,6 +14,7 @@ add_trans_prob(
   time_var = "tend",
   interval_length = "intlen",
   transition = "transition",
+  check_grouping = TRUE,
   ...
 )
 }
@@ -51,6 +52,16 @@ contains the interval length in \code{newdata}.}
 \item{transition}{\code{Character}, defaults to \code{"transition"}.
 contains the transition labels in \code{newdata}.}
 
+\item{check_grouping}{Logical. If \code{TRUE} (default), the function checks
+that \code{newdata} is grouped so that the time variable is unique within
+each group once \code{transition} is part of the grouping, and
+\code{\link[base]{stop}}s otherwise, guarding against silently accumulating
+transition probabilities across distinct covariate profiles (a forgotten
+\code{group_by()}). Set to \code{FALSE} to skip the check. As for the
+other cumulative \code{add_*} functions, hand-built grids stacking
+profiles with disjoint time grids cannot be detected (see
+\code{\link{add_cumu_hazard}}).}
+
 \item{...}{Further arguments passed to underlying methods.}
 }
 \description{
@@ -60,10 +71,13 @@ The function builds on cumulative hazards \code{cumu_hazard} and \code{mgcv::gam
 }
 \details{
 When computing transition probabilities for multiple groups, the input data must
-be grouped via \code{group_by()} before calling this function. Omitting
-\code{group_by()} will not produce an error or warning but will return
-silently incorrect results, as the transition probability will be accumulated
-over the entire dataset rather than within each group.
+be grouped via \code{group_by()} before calling this function. If \code{newdata}
+still contains several covariate profiles per (group, transition) -- i.e.\
+repeated \code{time_var} values within a group once \code{transition} is added
+to the grouping, typically a forgotten \code{group_by()} -- the function now
+\strong{stops with an error} rather than returning silently incorrect results,
+as the transition probability would otherwise be accumulated across profiles
+rather than within each group.
 
 The returned data contains one boundary row per group and transition at
 \code{time_var = 0} for plotting transition probabilities from the time
diff --git a/man/stop_if_undergrouped_for_cumulation.Rd b/man/stop_if_undergrouped_for_cumulation.Rd
new file mode 100644
index 00000000..5679bf41
--- /dev/null
+++ b/man/stop_if_undergrouped_for_cumulation.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/add-functions.R
+\name{stop_if_undergrouped_for_cumulation}
+\alias{stop_if_undergrouped_for_cumulation}
+\title{Guard against silently-wrong cumulative results from mis-grouped \code{newdata}}
+\usage{
+stop_if_undergrouped_for_cumulation(newdata, time_var, fun = "add_*")
+}
+\description{
+The cumulative post-processing functions (\code{add_cumu_hazard},
+\code{add_surv_prob}, \code{add_cif}, \code{add_trans_prob}) accumulate hazards over the
+rows of each group of \code{newdata}. If several covariate profiles (or causes /
+transitions) share a group -- typically a forgotten \code{group_by()} -- the
+accumulation silently runs across them and returns wrong curves. A correctly
+grouped prediction grid has a unique \code{time_var} value per group, so a
+duplicate within a group signals this situation. (The converse does not
+hold: profiles stacked with \emph{disjoint} time grids have unique \code{time_var}
+values and pass undetected -- such input is indistinguishable from a
+legitimate single profile with time-varying covariates, so it cannot be
+guarded against. Grids built with \code{make_newdata()} always share time values
+across profiles and are therefore always caught.) This generalises
+the \code{stopifnot()} guard used for the RMST helper to the whole \verb{add_*}
+family, turning a silent numerical error into an explicit one.
+}
+\keyword{internal}
diff --git a/tests/testthat/test-grouping-guard.R b/tests/testthat/test-grouping-guard.R
new file mode 100644
index 00000000..80dc74dd
--- /dev/null
+++ b/tests/testthat/test-grouping-guard.R
@@ -0,0 +1,119 @@
+# Regression tests for the safeguard that turns silently-wrong cumulative
+# results (from a forgotten group_by()) into an explicit error. See
+# stop_if_undergrouped_for_cumulation(). Uses only synthetic data.
+
+make_strata_ped <- function(seed = 1L, n = 400L) {
+  set.seed(seed)
+  df <- data.frame(
+    time = rexp(n, 0.1),
+    status = rbinom(n, 1, 0.7),
+    grp = factor(sample(c("a", "b"), n, replace = TRUE))
+  )
+  as_ped(df, Surv(time, status) ~ grp)
+}
+
+make_cr_ped <- function(seed = 2L, n = 300L) {
+  set.seed(seed)
+  df <- data.frame(time = rexp(n, 0.2), status = sample(0:2, n, replace = TRUE))
+  as_ped(df, Surv(time, status) ~ ., id = "id") %>%
+    dplyr::mutate(cause = as.factor(cause))
+}
+
+test_that("add_surv_prob / add_cumu_hazard error on ungrouped multi-profile newdata", {
+  ped <- make_strata_ped()
+  pam <- pamm(ped_status ~ s(tend) + grp, data = ped)
+  nd <- ped %>% make_newdata(tend = unique(tend), grp = unique(grp))
+
+  expect_error(add_surv_prob(nd, pam), "repeated 'tend' values within a group")
+  expect_error(
+    add_cumu_hazard(nd, pam),
+    "repeated 'tend' values within a group"
+  )
+
+  # correctly grouped -> no error
+  expect_error(nd %>% dplyr::group_by(grp) %>% add_surv_prob(pam), NA)
+  expect_error(nd %>% dplyr::group_by(grp) %>% add_cumu_hazard(pam), NA)
+})
+
+test_that("single-profile newdata needs no grouping (no false positive)", {
+  ped <- make_strata_ped()
+  pam0 <- pamm(ped_status ~ s(tend), data = ped)
+  nd0 <- ped %>% make_newdata(tend = unique(tend))
+  expect_error(add_surv_prob(nd0, pam0), NA)
+  expect_error(add_cumu_hazard(nd0, pam0), NA)
+})
+
+test_that("check_grouping = FALSE opts out of the safeguard", {
+  ped <- make_strata_ped()
+  pam <- pamm(ped_status ~ s(tend) + grp, data = ped)
+  nd <- ped %>% make_newdata(tend = unique(tend), grp = unique(grp))
+  # opting out reproduces the old (unsafe) behaviour without erroring
+  expect_error(add_surv_prob(nd, pam, check_grouping = FALSE), NA)
+  expect_error(add_cumu_hazard(nd, pam, check_grouping = FALSE), NA)
+})
+
+test_that("add_cif errors on ungrouped newdata and works when grouped by cause", {
+  pedc <- make_cr_ped()
+  pamc <- pamm(ped_status ~ s(tend, by = cause), data = pedc)
+  ndc <- pedc %>% make_newdata(tend = unique(tend), cause = unique(cause))
+
+  expect_error(add_cif(ndc, pamc), "repeated 'tend' values within a group")
+  expect_error(ndc %>% dplyr::group_by(cause) %>% add_cif(pamc, ci = FALSE), NA)
+  # check_grouping = FALSE skips the profile guard, but get_cif()'s own
+  # (pre-existing) cause check still catches ungrouped multi-cause input
+  expect_error(
+    add_cif(ndc, pamc, ci = FALSE, check_grouping = FALSE),
+    "group by cause"
+  )
+})
+
+test_that("add_trans_prob errors when covariate profiles are not separated", {
+  skip_if_not_installed("mstate")
+  data("prothr", package = "mstate")
+  prothr <- prothr %>%
+    dplyr::mutate(
+      transition = as.factor(paste0(from, "->", to)),
+      treat = as.factor(treat)
+    ) %>%
+    dplyr::filter(Tstart != Tstop, id <= 80) %>%
+    dplyr::select(-trans)
+  ped <- as_ped(
+    prothr,
+    Surv(Tstart, Tstop, status) ~ .,
+    transition = "transition",
+    id = "id",
+    timescale = "calendar"
+  )
+  pam <- mgcv::bam(
+    ped_status ~ s(tend, by = transition) + transition * treat,
+    data = ped,
+    family = poisson(),
+    offset = offset,
+    method = "fREML",
+    discrete = TRUE
+  )
+  ndf <- make_newdata(
+    ped,
+    tend = unique(tend),
+    treat = unique(treat),
+    transition = unique(transition)
+  )
+
+  # grouped only by transition, but two `treat` profiles remain -> error
+  expect_error(
+    ndf %>% dplyr::group_by(transition) %>% add_trans_prob(pam),
+    "repeated 'tend' values within a group"
+  )
+  # grouped by both profile columns -> no error
+  expect_error(
+    ndf %>% dplyr::group_by(treat, transition) %>% add_trans_prob(pam),
+    NA
+  )
+  # check_grouping = FALSE opts out (old unsafe behaviour, no error)
+  expect_error(
+    ndf %>%
+      dplyr::group_by(transition) %>%
+      add_trans_prob(pam, check_grouping = FALSE),
+    NA
+  )
+})
diff --git a/tests/testthat/test-interval-censored.R b/tests/testthat/test-interval-censored.R
index 68a1eb69..09b57506 100644
--- a/tests/testthat/test-interval-censored.R
+++ b/tests/testthat/test-interval-censored.R
@@ -329,7 +329,7 @@ test_that("competing-risks IC pipeline yields valid pooled CIFs", {
   cif2 <- add_cif(nd, fcr, ci = FALSE, nsim = 120)
   expect_equal(cif1$cif, cif2$cif, tolerance = 1e-12)
   expect_false(any(c("cif_lower", "cif_upper") %in% names(cif1)))
-  expect_error(add_cif(dplyr::ungroup(nd), fcr, ci = FALSE), "group by cause")
+  expect_error(add_cif(dplyr::ungroup(nd), fcr, ci = FALSE), "group_by")
 })
 
 test_that("competing-risks imputation samples exact unknown causes", {
@@ -554,26 +554,47 @@ test_that("add_inspections right-censoring is non-informative", {
   expect_error(add_inspections(sdf_na, rate = 1), "missing")
 })
 
-test_that("pooled pamm_ic adders warn on undergrouped newdata", {
+test_that("pooled pamm_ic adders error on undergrouped newdata", {
   set.seed(11)
   df <- data.frame(z = c(rep(0, 200), rep(1, 200)))
   sdf <- sim_pexp(~ -2 + 1.0 * z, df, cut = seq(0, 8, by = 0.25))
   icd <- add_inspections(sdf, rate = 1, max_time = 8)
   cut <- seq(0, 8, by = 0.5)
-  fit <- pamm_ic(Surv(L, R, type = "interval2") ~ s(tend) + z, icd, cut = cut, m = 2)
+  fit <- pamm_ic(
+    Surv(L, R, type = "interval2") ~ s(tend) + z,
+    icd,
+    cut = cut,
+    m = 2
+  )
   ped <- as_ped(icd, Surv(L, R, type = "interval2") ~ z, cut = cut)
 
-  # multiple covariate groups but no group_by() -> distinct curves get pooled
+  # multiple covariate groups but no group_by() -> distinct curves would get
+  # pooled; the safeguard now errors instead of silently returning wrong curves
   nd_multi <- make_newdata(ped, tend = unique(tend), z = c(0, 1))
-  expect_warning(add_surv_prob(nd_multi, fit, ci = FALSE), "group_by")
-  expect_warning(add_cumu_hazard(nd_multi, fit, ci = FALSE), "group_by")
+  expect_error(add_surv_prob(nd_multi, fit, ci = FALSE), "group_by")
+  expect_error(add_cumu_hazard(nd_multi, fit, ci = FALSE), "group_by")
 
-  # group_by() silences the warning and yields correct per-group curves
+  # check_grouping = FALSE opts out for the pooled methods as well
+  expect_error(
+    add_surv_prob(nd_multi, fit, ci = FALSE, check_grouping = FALSE),
+    NA
+  )
+  expect_error(
+    add_cumu_hazard(nd_multi, fit, ci = FALSE, check_grouping = FALSE),
+    NA
+  )
+
+  # group_by() silences the error and yields correct per-group curves
   nd_g <- dplyr::group_by(nd_multi, z)
-  expect_warning(s <- add_surv_prob(nd_g, fit, ci = TRUE, nsim = 40), regexp = NA)
+  expect_warning(
+    s <- add_surv_prob(nd_g, fit, ci = TRUE, nsim = 40),
+    regexp = NA
+  )
   s <- dplyr::ungroup(s)
   expect_true(all(tapply(s$surv_prob, s$z, function(v) all(diff(v) <= 1e-8))))
-  expect_true(all(s$surv_lower <= s$surv_prob + 1e-8 & s$surv_prob <= s$surv_upper + 1e-8))
+  expect_true(all(
+    s$surv_lower <= s$surv_prob + 1e-8 & s$surv_prob <= s$surv_upper + 1e-8
+  ))
 
   # a single-curve prediction grid must not warn
   nd1 <- make_newdata(ped, tend = unique(tend))