diff --git a/DESCRIPTION b/DESCRIPTION
index 5fa885cb..798fc786 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: sccomp
 Type: Package
 Title: Differential Composition and Variability Analysis for Single-Cell Data
-Version: 2.5.0
-Date: 2026-05-11
+Version: 2.5.1
+Date: 2026-08-07
 Authors@R: c(person("Stefano", "Mangiola", email = "stefano.mangiola@unimelb.edu.au", role = c("aut", "cre")), person("Alexandra J.", "Roth-Schulze", role = "aut"), person("Marie", "Trussart", role = "aut"), person("Enrique", "Zozaya-Valdés", role = "aut"), person("Mengyao", "Ma", role = "aut"), person("Zijie", "Gao", role = "aut"), person("Alan F.", "Rubin", role = "aut"), person("Terence P.", "Speed", role = "aut"), person("Heejung", "Shim", role = "aut"), person("Anthony T.", "Papenfuss", role = "aut"))
 Description: Comprehensive R package for differential composition and variability analysis in single-cell RNA sequencing, CyTOF, and microbiome data. Provides robust Bayesian modeling with outlier detection, random effects, and advanced statistical methods for cell type proportion analysis. Features include probabilistic outlier identification, mixed-effect modeling, differential variability testing, and comprehensive visualization tools. Perfect for cancer research, immunology, developmental biology, and single-cell genomics applications.
 License: GPL-3
diff --git a/NAMESPACE b/NAMESPACE
index c5866835..5b81b228 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -123,8 +123,10 @@ importFrom(posterior,summarise_draws)
 importFrom(purrr,as_mapper)
 importFrom(purrr,map)
 importFrom(purrr,map2)
+importFrom(purrr,map2_chr)
 importFrom(purrr,map2_dfc)
 importFrom(purrr,map2_lgl)
+importFrom(purrr,map_chr)
 importFrom(purrr,map_dbl)
 importFrom(purrr,map_dfr)
 importFrom(purrr,map_int)
diff --git a/R/sccomp_estimate.R b/R/sccomp_estimate.R
index 1a326688..9d62c962 100644
--- a/R/sccomp_estimate.R
+++ b/R/sccomp_estimate.R
@@ -1036,9 +1036,8 @@ sccomp_glm_data_frame_counts = function(.data,
       random_effect_elements
     )
   
-  # Print design matrix
-  message(sprintf("sccomp says: the composition design matrix has columns: %s", data_for_model$X %>% colnames %>% paste(collapse=", ")))
-  message(sprintf("sccomp says: the variability design matrix has columns: %s", data_for_model$Xa %>% colnames %>% paste(collapse=", ")))
+  # Preview fixed- and random-effect design matrices.
+  message_design_matrices(data_for_model)
   
   # Force outliers, Get the truncation index
   data_for_model$user_forced_truncation_not_idx = 
diff --git a/R/sccomp_remove_outliers.R b/R/sccomp_remove_outliers.R
index 9b101747..9b2f2776 100644
--- a/R/sccomp_remove_outliers.R
+++ b/R/sccomp_remove_outliers.R
@@ -457,9 +457,8 @@ sccomp_remove_outliers.sccomp_tbl = function(.estimate,
   
   message("sccomp says: outlier-free model fitting - step 2/2")
   
-  # Print design matrix
-  message(sprintf("sccomp says: the composition design matrix has columns: %s", data_for_model$X |> colnames() |> paste(collapse=", ")))
-  message(sprintf("sccomp says: the variability design matrix has columns: %s", data_for_model$Xa |> colnames() |> paste(collapse=", ")))
+  # Preview fixed- and random-effect design matrices.
+  message_design_matrices(data_for_model)
   
   fit3 =
     data_for_model |>
diff --git a/R/smooths.R b/R/smooths.R
index 75de99a9..038b1cb6 100644
--- a/R/smooths.R
+++ b/R/smooths.R
@@ -141,6 +141,8 @@ strip_random_effect_terms <- function(fm) {
 #'     each block its sccomp RE-slot label: just the smooth's label for
 #'     single-penalty smooths, or `<label>__b<b>` for multi-penalty smooths
 #'     (so per-block `gfi` rownames stay unique).
+#'   * `Xr_slot_terms`      — character vector parallel to `Xr_list` describing
+#'     in plain words what each block models, for user-facing messaging.
 #'   * `smooth_specs`       — list of mgcv `smoothCon` objects (for predict).
 #'   * `smooth_re_objs`     — list of mgcv `smooth2random` objects (rotations).
 #'   When `fm` has no smooth specials, lists are length-0 and
@@ -156,6 +158,7 @@ parse_formula_smooths <- function(fm, data) {
     Xr_list            = list(),
     Xr_to_smooth       = integer(0),
     Xr_slot_labels     = character(0),
+    Xr_slot_terms      = character(0),
     smooth_specs       = list(),
     smooth_re_objs     = list()
   )
@@ -182,6 +185,7 @@ parse_formula_smooths <- function(fm, data) {
   Xr_list        <- list()   # flat across smooths × blocks
   Xr_to_smooth   <- integer(0)
   Xr_slot_labels <- character(0)
+  Xr_slot_terms  <- character(0)
   
   for (k in seq_along(smooth_calls)) {
     label <- smooth_labels[k]
@@ -232,6 +236,11 @@ parse_formula_smooths <- function(fm, data) {
     block_names    <- names(re$rand)
     n_blocks_smooth <- length(block_names)
     
+    # `bs = "fs"` and `by = <factor>` both fit one curve per level of a
+    # grouping factor, which is worth spelling out to the user.
+    is_grouped <- inherits(sm, "fs.interaction") ||
+      !identical(as.character(sm$by), "NA")
+    
     for (b in seq_along(block_names)) {
       Xr_b <- re$rand[[block_names[b]]]
       attr(Xr_b, "s.label") <- NULL
@@ -240,10 +249,15 @@ parse_formula_smooths <- function(fm, data) {
       # disambiguate with `<label>__b<b>`.
       slot_label <- if (n_blocks_smooth == 1L) label
                     else sprintf("%s__b%d", label, b)
+      slot_term <- sprintf("%s %s", if (is_grouped) "grouped smooth" else "smooth", label)
+      if (n_blocks_smooth > 1L) {
+        slot_term <- sprintf("%s, penalty block %d/%d", slot_term, b, n_blocks_smooth)
+      }
       colnames(Xr_b) <- sprintf("%s___basis%02d", slot_label, seq_len(ncol(Xr_b)))
       Xr_list[[length(Xr_list) + 1L]] <- Xr_b
       Xr_to_smooth                    <- c(Xr_to_smooth, k)
       Xr_slot_labels                  <- c(Xr_slot_labels, slot_label)
+      Xr_slot_terms                   <- c(Xr_slot_terms, slot_term)
     }
     
     Xf_list[[k]]        <- Xf
@@ -258,6 +272,7 @@ parse_formula_smooths <- function(fm, data) {
     Xr_list            = Xr_list,
     Xr_to_smooth       = Xr_to_smooth,
     Xr_slot_labels     = Xr_slot_labels,
+    Xr_slot_terms      = Xr_slot_terms,
     smooth_specs       = smooth_specs,
     smooth_re_objs     = smooth_re_objs
   )
@@ -385,12 +400,14 @@ predict_smooth_at_newdata <- function(sm, re_obj, newdata) {
 #'
 #' @param Xr     A numeric matrix (N rows × K wiggly basis columns).
 #' @param label  Smooth label, used as the slot's single "factor" name.
+#' @param term   Plain-words description of what the slot models, used when
+#'   reporting the design to the user.
 #' @return A list with the same fields as `empty_re_slot` in
 #'   `data_spread_to_model_input()`.
 #'
 #' @keywords internal
 #' @noRd
-build_smooth_slot <- function(Xr, label) {
+build_smooth_slot <- function(Xr, label, term = label) {
   stopifnot(is.matrix(Xr) || is.data.frame(Xr))
   Xr <- as.matrix(Xr)
   K  <- ncol(Xr)
@@ -402,7 +419,8 @@ build_smooth_slot <- function(Xr, label) {
       ncol      = 0L,
       gfi       = matrix(integer(0), nrow = 0L, ncol = 0L),
       n_groups  = 0L,
-      n_factors = 0L
+      n_factors = 0L,
+      term      = term
     ))
   }
   
@@ -420,7 +438,8 @@ build_smooth_slot <- function(Xr, label) {
     ncol      = K,
     gfi       = gfi,
     n_groups  = K,
-    n_factors = 1L
+    n_factors = 1L,
+    term      = term
   )
 }
 
diff --git a/R/utilities.R b/R/utilities.R
index dbb03465..7a8bd330 100755
--- a/R/utilities.R
+++ b/R/utilities.R
@@ -258,6 +258,83 @@ formula_to_random_effect_formulae <- function(fm) {
   
 }
 
+#' Format a compact preview of design-matrix parameters
+#'
+#' @param x A matrix-like object with column names.
+#' @param max_parameters Maximum number of parameter names to display.
+#'
+#' @return A single character string.
+#' @keywords internal
+#' @noRd
+format_design_matrix_preview <- function(x, max_parameters = 10L) {
+  parameters <- colnames(x)
+  n_parameters <- length(parameters)
+  shown <- utils::head(parameters, max_parameters)
+  suffix <- if (n_parameters > max_parameters) {
+    sprintf(", ... (+%d more)", n_parameters - max_parameters)
+  } else {
+    ""
+  }
+
+  sprintf(
+    "%d parameter%s: %s%s",
+    n_parameters,
+    if (n_parameters == 1L) "" else "s",
+    paste(shown, collapse = ", "),
+    suffix
+  )
+}
+
+
+#' Report the fixed- and random-effect design matrices
+#'
+#' @param model_input The list returned by `data_spread_to_model_input()`.
+#' @param max_parameters Maximum number of parameter names shown per matrix.
+#'
+#' @return `model_input`, invisibly.
+#' @keywords internal
+#' @noRd
+message_design_matrices <- function(model_input, max_parameters = 10L) {
+  lines <- c(
+    "sccomp says: model design matrices",
+    sprintf(
+      "  composition X - %s",
+      format_design_matrix_preview(model_input$X, max_parameters)
+    ),
+    sprintf(
+      "  variability Xa - %s",
+      format_design_matrix_preview(model_input$Xa, max_parameters)
+    )
+  )
+  
+  active_slots <- which(model_input$ncol_X_random_eff > 0L)
+  
+  if (length(active_slots) == 0L) {
+    lines <- c(lines, "  random effects - none")
+  } else {
+    # Fits produced before the design terms were recorded still report their
+    # matrices, just without the originating formula clause.
+    recorded <- model_input$random_effect_design_terms
+    slot_terms <- rep(NA_character_, length(model_input$ncol_X_random_eff))
+    if (!is.null(recorded)) slot_terms[recorded$slot] <- recorded$term
+    slot_terms[is.na(slot_terms)] <- sprintf("slot %d", which(is.na(slot_terms)))
+    
+    lines <- c(lines, sprintf(
+      "  random effect X%d [%s] - %s",
+      active_slots,
+      slot_terms[active_slots],
+      map_chr(
+        sprintf("X_random_effect_%d", active_slots),
+        ~ format_design_matrix_preview(model_input[[.x]], max_parameters)
+      )
+    ))
+  }
+  
+  message(paste(lines, collapse = "\n"))
+  invisible(model_input)
+}
+
+
 #' Formula parser
 #'
 #' @param fm A formula
@@ -932,6 +1009,9 @@ calculate_na_fraction_contribution = function(my_design_matrix, na_cols, design_
 #' @importFrom stringr str_remove_all
 #' @importFrom purrr reduce
 #' @importFrom purrr map_int
+#' @importFrom purrr map_chr
+#' @importFrom purrr map2_chr
+#' @importFrom purrr pmap
 #' @importFrom stats as.formula
 #'
 #' Match variability to composition design columns
@@ -1111,10 +1191,11 @@ data_spread_to_model_input =
       ncol     = 0L,
       gfi      = matrix(integer(0), nrow = 0, ncol = 0),
       n_groups = 0L,
-      n_factors = 0L
+      n_factors = 0L,
+      term     = NA_character_
     )
     
-    prepare_re_slot = function(design_matrix_tbl, sample_name) {
+    prepare_re_slot = function(design_matrix_tbl, sample_name, term) {
       X = design_matrix_tbl |> column_to_rownames(sample_name)
       
       is_NA_col = str_detect(colnames(X), "___NA$")
@@ -1138,7 +1219,8 @@ data_spread_to_model_input =
         ncol      = ncol(X),
         gfi       = gfi,
         n_groups  = ncol(gfi),
-        n_factors = nrow(gfi)
+        n_factors = nrow(gfi),
+        term      = term
       )
     }
     
@@ -1149,6 +1231,13 @@ data_spread_to_model_input =
         mutate(design = map2(
           formula, grouping,
           ~ get_random_effect_design3(.data_spread, .x, .y, !!.sample)
+        )) |>
+        
+        # Keep the clause the user wrote, e.g. `(1 + age | donor)`, next to the
+        # slot it produces, so reporting never re-parses the formula.
+        mutate(term = map2_chr(
+          formula, grouping,
+          ~ sprintf("(%s | %s)", deparse1(.x[[2]]), .y)
         ))
       
       is_random_effect = 1
@@ -1163,8 +1252,11 @@ data_spread_to_model_input =
             mutate(across(everything(), ~ .x |> replace_na(0)))
         ))
       
-      re_slots = random_effect_grouping$design_matrix |>
-        map(prepare_re_slot, sample_name = quo_name(.sample))
+      re_slots = map2(
+        random_effect_grouping$design_matrix,
+        random_effect_grouping$term,
+        ~ prepare_re_slot(.x, quo_name(.sample), .y)
+      )
       
     } else {
       is_random_effect = 0
@@ -1179,13 +1271,17 @@ data_spread_to_model_input =
     # point of view.
     # Each penalised basis block becomes one RE slot. `parse_formula_smooths()`
     # already computed the per-block slot label (`<label>` for single-penalty
-    # smooths, `<label>__b<b>` for multi-penalty ones) so we just zip the two
-    # parallel lists. Smooth terms also flip `is_random_effect = 1` because
-    # they reuse the Stan RE machinery even when the user wrote no (... | g).
+    # smooths, `<label>__b<b>` for multi-penalty ones) and its plain-words
+    # description, so we just zip the parallel lists. Smooth terms also flip
+    # `is_random_effect = 1` because they reuse the Stan RE machinery even
+    # when the user wrote no (... | g).
     if (length(smooth_pieces$Xr_list) > 0) {
-      smooth_slots = map2(
-        smooth_pieces$Xr_list,
-        smooth_pieces$Xr_slot_labels,
+      smooth_slots = pmap(
+        list(
+          smooth_pieces$Xr_list,
+          smooth_pieces$Xr_slot_labels,
+          smooth_pieces$Xr_slot_terms
+        ),
         build_smooth_slot
       )
       re_slots         = c(re_slots, smooth_slots)
@@ -1227,6 +1323,12 @@ data_spread_to_model_input =
     n_groups                          = map_int(re_slots, "n_groups")
     how_many_factors_in_random_design = map_int(re_slots, "n_factors")
     
+    # What each occupied slot models, in the user's own words. Empty padding
+    # slots model nothing, so they are dropped.
+    random_effect_design_terms =
+      tibble(slot = seq_along(re_slots), term = map_chr(re_slots, "term")) |>
+      filter(!is.na(term))
+    
     
     y = .data_spread %>% select(-any_of(factor_names), -exposure, -!!.grouping_for_random_effect) %>% column_to_rownames(quo_name(.sample)) %>% as.matrix()
     
@@ -1275,6 +1377,7 @@ data_spread_to_model_input =
         group_factor_indexes_for_covariance_2 = group_factor_indexes_for_covariance_2,
         group_factor_indexes_for_covariance_3 = group_factor_indexes_for_covariance_3,
         group_factor_indexes_for_covariance_4 = group_factor_indexes_for_covariance_4,
+        random_effect_design_terms = random_effect_design_terms,
         
         # For parallel chains
         grainsize = 1,
diff --git a/inst/NEWS.rd b/inst/NEWS.rd
index cf79f2a3..3dba0af8 100644
--- a/inst/NEWS.rd
+++ b/inst/NEWS.rd
@@ -1,6 +1,13 @@
 \name{NEWS}
 \title{News for Package \pkg{sccomp}}
 
+\section{News in version 2.5.1}{
+\itemize{
+    \item \strong{Clearer design-matrix reporting.} \code{sccomp_estimate()} and \code{sccomp_remove_outliers()} now report the composition and variability design matrices with their parameter count followed by a preview of the first 10 parameters, so wide designs no longer flood the console.
+    \item Random-effect design matrices are now reported as well, one per line, each labelled with what it models: the clause the user wrote (e.g. \code{(1 + type | donor)}), a smooth term, or a single penalty block of a grouped smooth such as \code{s(x, group, bs = "fs")}. This makes it visible how many random-effect slots a formula consumes.
+    \item \code{model_input} gained a \code{random_effect_design_terms} table mapping each occupied random-effect slot to the term it models.
+}}
+
 \section{News in version 2.1.34}{
 \itemize{
     \item Exported \code{sccomp_scatterplot()} for visualising cell-group proportions against a continuous covariate, complementing \code{sccomp_boxplot()} for discrete factors. The function accepts \code{.data}, \code{factor}, \code{significance_threshold}, and \code{remove_unwanted_effects}, and is used by the \code{plot()} method for numeric covariates.
diff --git a/tests/testthat/test-design-matrix-messages.R b/tests/testthat/test-design-matrix-messages.R
new file mode 100644
index 00000000..b5225999
--- /dev/null
+++ b/tests/testthat/test-design-matrix-messages.R
@@ -0,0 +1,110 @@
+library(testthat)
+library(sccomp)
+
+
+test_that("design matrix parameter previews are compact", {
+  x <- matrix(0, nrow = 2, ncol = 12)
+  colnames(x) <- paste0("parameter_", seq_len(ncol(x)))
+  
+  expect_equal(
+    sccomp:::format_design_matrix_preview(x),
+    paste0(
+      "12 parameters: ",
+      paste0("parameter_", seq_len(10), collapse = ", "),
+      ", ... (+2 more)"
+    )
+  )
+})
+
+
+test_that("design matrix messages identify every modelled matrix", {
+  model_input <- list(
+    X = matrix(
+      0,
+      nrow = 2,
+      ncol = 2,
+      dimnames = list(NULL, c("(Intercept)", "conditiontreated"))
+    ),
+    Xa = matrix(
+      0,
+      nrow = 2,
+      ncol = 1,
+      dimnames = list(NULL, "(Intercept)")
+    ),
+    ncol_X_random_eff = c(4L, 3L, 0L, 0L),
+    X_random_effect_1 = matrix(
+      0,
+      nrow = 2,
+      ncol = 4,
+      dimnames = list(
+        NULL,
+        c(
+          "(Intercept)___A", "(Intercept)___B",
+          "age___A", "age___B"
+        )
+      )
+    ),
+    X_random_effect_2 = matrix(
+      0,
+      nrow = 2,
+      ncol = 3,
+      dimnames = list(
+        NULL,
+        c(
+          "s(age, k = 5)___basis01",
+          "s(age, k = 5)___basis02",
+          "s(age, k = 5)___basis03"
+        )
+      )
+    )
+  )
+  model_input$random_effect_design_terms <- tibble::tibble(
+    slot = c(1L, 2L),
+    term = c("(1 + age | donor)", "smooth s(age, k = 5)")
+  )
+  
+  messages <- character(0)
+  withCallingHandlers(
+    sccomp:::message_design_matrices(model_input),
+    message = function(m) {
+      messages <<- c(messages, conditionMessage(m))
+      invokeRestart("muffleMessage")
+    }
+  )
+  
+  expect_length(messages, 1L)
+  expect_match(messages, "composition X - 2 parameters", fixed = TRUE)
+  expect_match(messages, "variability Xa - 1 parameter", fixed = TRUE)
+  expect_match(
+    messages,
+    "random effect X1 [(1 + age | donor)] - 4 parameters",
+    fixed = TRUE
+  )
+  expect_match(
+    messages,
+    "random effect X2 [smooth s(age, k = 5)] - 3 parameters",
+    fixed = TRUE
+  )
+  expect_false(grepl("random effect X3", messages, fixed = TRUE))
+})
+
+
+test_that("design matrix messages survive fits without recorded terms", {
+  model_input <- list(
+    X = matrix(0, nrow = 2, ncol = 1, dimnames = list(NULL, "(Intercept)")),
+    Xa = matrix(0, nrow = 2, ncol = 1, dimnames = list(NULL, "(Intercept)")),
+    ncol_X_random_eff = c(1L, 0L, 0L, 0L),
+    X_random_effect_1 = matrix(
+      0,
+      nrow = 2,
+      ncol = 1,
+      dimnames = list(NULL, "(Intercept)___A")
+    )
+  )
+  
+  expect_message(
+    sccomp:::message_design_matrices(model_input),
+    "random effect X1 [slot 1] - 1 parameter",
+    fixed = TRUE
+  )
+})