diff --git a/DESCRIPTION b/DESCRIPTION
index da8a8cb7..59410f74 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,14 +1,13 @@
 Package: sccomp
 Type: Package
 Title: Differential Composition and Variability Analysis for Single-Cell Data
-Version: 2.1.32
+Version: 2.1.33
 Date: 2026-05-11
 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
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.3
 Depends: 
     R (>= 4.3.0),
     instantiate (>= 0.2.3)
@@ -39,6 +38,7 @@ Imports:
     SingleCellExperiment,
     posterior
 Suggests:
+    mgcv,
     knitr,
     rmarkdown,
     BiocStyle,
@@ -63,3 +63,4 @@ Config/testthat/edition: 3
 Config/testthat/parallel: false
 Config/testthat/snapshot/parallel: false
 Config/testthat/snapshot/parallel/workers: 2
+Config/roxygen2/version: 8.0.0
diff --git a/R/sccomp_remove_outliers.R b/R/sccomp_remove_outliers.R
index b9133be6..9b101747 100644
--- a/R/sccomp_remove_outliers.R
+++ b/R/sccomp_remove_outliers.R
@@ -196,8 +196,7 @@ sccomp_remove_outliers.sccomp_tbl = function(.estimate,
     
     # This is for the new data generation with selected factors to do adjustment
     data = 
-      .estimate |>
-      attr("model_input") |> 
+      data_for_model |> 
       c(
         list(
           # Add subset of coefficients
diff --git a/R/sccomp_replicate.R b/R/sccomp_replicate.R
index fbc87b66..941e9b3f 100644
--- a/R/sccomp_replicate.R
+++ b/R/sccomp_replicate.R
@@ -135,6 +135,16 @@ sccomp_replicate.sccomp_tbl = function(fit,
 #' @param formula_variability Formula for the variability model
 #' @param new_data New data to generate predictions for. If NULL, uses the original data
 #' @param original_count_data Original count data from the model
+#' @param smooth_results Optional named list of smooth-term metadata captured
+#'   at fit time. In the standard call path this is read from
+#'   `get_smooth_results(.data)` (an attribute on `model_input`, not Stan data).
+#'   Contains
+#'   `smooth_specs`, the `mgcv::smoothCon()` objects used to evaluate bases at
+#'   `new_data`; `smooth_re_objs`, the parallel `mgcv::smooth2random()` results
+#'   used to recover the fitted `Xf` / `Xr` parameterisation; and
+#'   `smooth_labels`, the original term labels used to name generated smooth
+#'   columns so they match the fit-time design. `NULL` when the composition
+#'   formula has no smooths.
 #' 
 #' @return A list containing:
 #' - model_input: The prepared model input data
@@ -156,7 +166,8 @@ prepare_replicate_data = function(X,
                                   original_formula_composition,
                                   formula_variability,
                                   new_data = NULL,
-                                  original_count_data) {
+                                  original_count_data,
+                                  smooth_results = NULL) {
   
   
   .sample = enquo(.sample)
@@ -226,13 +237,17 @@ prepare_replicate_data = function(X,
   
   new_data =  old_data |> bind_rows( new_data ) 
   
+  # Smooth columns are evaluated separately via PredictMat below; keep the
+  # random-effect clauses so the later RE parsing still sees the same formula.
+  formula_composition = strip_smooth_terms(formula_composition)
+  
   new_X =
     new_data |>
     get_design_matrix(
       # Drop random intercept
       formula_composition |>
+        strip_random_effect_terms() |>
         as.character() |>
-        str_remove_all("\\+ ?\\(.+\\|.+\\)") |>
         paste(collapse="") |>
         as.formula(),
       !!.sample, 
@@ -242,6 +257,17 @@ prepare_replicate_data = function(X,
     # Remove columns that are not in the original design matrix
     .[,colnames(.) %in% colnames(X), drop=FALSE]
   
+  # Evaluate any smooth bases on the replicate rows and merge the resulting
+  # design pieces (unpenalised columns appended to `new_X`, one RE slot per
+  # penalised block) using the fit-time `smoothCon` / `smooth2random` objects.
+  smooth_design = build_smooth_replicate_design(
+    parametric_X   = new_X,
+    new_data_tail  = new_data |> tail(nrow_new_data),
+    smooth_results = smooth_results
+  )
+  new_X                  = smooth_design$new_X
+  smooth_replicate_slots = smooth_design$smooth_replicate_slots
+  
   # Check that all effect combination were present when the model was fitted
   check_missing_parameters(
     new_X |> colnames(), 
@@ -263,8 +289,8 @@ prepare_replicate_data = function(X,
     get_design_matrix(
       # Drop random intercept
       formula_variability |>
+        strip_random_effect_terms() |>
         as.character() |>
-        str_remove_all("\\+ ?\\(.+\\|.+\\)") |>
         paste(collapse="") |>
         as.formula(),
       !!.sample, 
@@ -363,6 +389,25 @@ prepare_replicate_data = function(X,
   
   replicate_slots = map(seq_len(4L), build_replicate_slot)
   
+  # Append smooth-derived replicate slots (one per smooth term in the
+  # composition formula). They occupy whichever slots come after the
+  # explicit RE clauses, mirroring the slot ordering used at fit time.
+  if (length(smooth_replicate_slots) > 0) {
+    n_explicit_re = length(original_grouping_names)
+    n_smooth      = length(smooth_replicate_slots)
+    n_used        = n_explicit_re + n_smooth
+    if (n_used > 4L) {
+      stop(sprintf(
+        "sccomp says: the replicate model needs %d RE slot(s) but only 4 are available.",
+        n_used
+      ))
+    }
+    # Replace placeholder slots `[n_explicit_re + 1 .. n_used]` with smooths.
+    for (k in seq_len(n_smooth)) {
+      replicate_slots[[n_explicit_re + k]] = smooth_replicate_slots[[k]]
+    }
+  }
+  
   # setup default unknown_grouping variable for generated quantities
   unknown_grouping = rep(0L, 4L)
   
@@ -464,7 +509,8 @@ replicate_data = function(.data,
     original_count_data =
       .data |>
       attr("count_data") |>
-      .subset(!!.sample) 
+      .subset(!!.sample),
+    smooth_results = get_smooth_results(.data)
   )
   
   # Original input
diff --git a/R/smooths.R b/R/smooths.R
new file mode 100644
index 00000000..75de99a9
--- /dev/null
+++ b/R/smooths.R
@@ -0,0 +1,520 @@
+# Smooth-term utilities for sccomp formulas (brms-style splines)
+#
+# This file implements the R-side machinery that turns smooth terms
+# (`s()`, `t2()`) in a sccomp formula into design-matrix columns that the
+# existing Stan model can consume without modification.
+#
+# The decomposition is the standard mgcv / brms one:
+#   * `Xf` — unpenalised null-space columns of the smooth (linear part).
+#     These are appended to the fixed-effect design matrix `X` and
+#     estimated as ordinary `beta` coefficients.
+#   * `Xr` — penalised "wiggly" basis columns. These occupy a sccomp
+#     random-effect slot with `n_factors = 1` and `n_groups = ncol(Xr)`,
+#     so the slot's per-cell-group SD plays the role of brms's `sds_*`.
+#
+# Built via:
+#   1. `mgcv::smoothCon(s(...), data, absorb.cons = TRUE,
+#                       diagonal.penalty = TRUE)[[1]]` — basis + penalty.
+#   2. `mgcv::smooth2random(sm, vnames = "", type = 2)` — re-parameterise
+#      into `Xf` + standardised `rand$Xr` (so the wiggly coefficients have
+#      a `N(0, sds * I)` prior).
+#
+# At prediction time, `predict_smooth_at_newdata()` evaluates the basis at
+# new covariate values via `mgcv::PredictMat()` and re-applies the same
+# `trans.D` (diagonal rescaling) and `trans.U` (orthonormal rotation)
+# captured at fit time, so the new columns line up with the columns the
+# model was trained on.
+
+#' Detect smooth specials in a formula
+#'
+#' @param fm A one-sided formula.
+#' @return TRUE if `fm` contains any `s(...)` or `t2(...)` term.
+#' @keywords internal
+#' @noRd
+has_smooth_terms <- function(fm) {
+  if (is.null(fm)) return(FALSE)
+  trm <- stats::terms(fm, specials = c("s", "t2"))
+  length(unlist(attr(trm, "specials"))) > 0
+}
+
+
+#' Strip smooth specials from a formula
+#'
+#' Returns a new formula with `s()` / `t2()` term labels removed. All other
+#' terms, including random-effect clauses (`(... | g)`), are preserved. The
+#' intercept setting (`+ 0` / `- 1`) is preserved. If nothing remains, returns
+#' `~ 1` (or `~ 0`).
+#'
+#' @param fm A one-sided formula.
+#' @return A one-sided formula without smooth term labels.
+#' @keywords internal
+#' @noRd
+strip_smooth_terms <- function(fm) {
+  if (is.null(fm)) return(fm)
+  trm <- stats::terms(fm, specials = c("s", "t2"))
+  smooth_idx <- unlist(attr(trm, "specials"))
+  
+  vars <- attr(trm, "variables")
+  smooth_labels <- if (length(smooth_idx) > 0) {
+    vapply(
+      smooth_idx + 1L,
+      function(i) deparse(vars[[i]], width.cutoff = 500L),
+      character(1)
+    )
+  } else {
+    character(0)
+  }
+  
+  all_labels <- attr(trm, "term.labels")
+  keep_labels <- setdiff(all_labels, smooth_labels)
+  keep_labels <- ifelse(grepl("\\|", keep_labels), paste0("(", keep_labels, ")"), keep_labels)
+  
+  has_intercept <- attr(trm, "intercept") == 1L
+  
+  if (length(keep_labels) == 0) {
+    return(if (has_intercept) ~ 1 else ~ 0)
+  }
+  
+  rhs <- paste(keep_labels, collapse = " + ")
+  if (!has_intercept) rhs <- paste0(rhs, " + 0")
+  stats::as.formula(paste("~", rhs), env = environment(fm))
+}
+
+
+#' Strip random-effect clauses from a formula
+#'
+#' Returns a formula with `(... | g)` term labels removed, preserving all
+#' ordinary fixed-effect terms and the intercept setting. This is used just
+#' before calling `model.matrix()`, which cannot evaluate sccomp/brms-style
+#' random-effect clauses.
+#'
+#' @param fm A one-sided formula.
+#' @return A one-sided formula without random-effect term labels.
+#' @keywords internal
+#' @noRd
+strip_random_effect_terms <- function(fm) {
+  if (is.null(fm)) return(fm)
+  trm <- stats::terms(fm)
+  
+  keep_labels <- attr(trm, "term.labels")
+  keep_labels <- keep_labels[!grepl("\\|", keep_labels)]
+  
+  has_intercept <- attr(trm, "intercept") == 1L
+  
+  if (length(keep_labels) == 0) {
+    return(if (has_intercept) ~ 1 else ~ 0)
+  }
+  
+  rhs <- paste(keep_labels, collapse = " + ")
+  if (!has_intercept) rhs <- paste0(rhs, " + 0")
+  stats::as.formula(paste("~", rhs), env = environment(fm))
+}
+
+
+#' Parse smooth terms from a sccomp formula
+#'
+#' For each `s(...)` / `t2(...)` term, builds an mgcv smooth via
+#' `smoothCon(..., absorb.cons = TRUE, diagonal.penalty = TRUE)` and the
+#' brms-style re-parameterisation via `smooth2random(..., type = 2)`. Single-
+#' penalty smooths (e.g. plain `s(x)`) yield one penalised block; multi-
+#' penalty smooths (e.g. `t2(x, y)`, `s(x, z, bs = "fs")`) yield 2+ blocks,
+#' one per penalty. Each block is mapped to its own sccomp random-effect
+#' slot, with its own per-cell-group smoothing SD — matching the way brms
+#' generates one `sds_*` parameter per penalty.
+#'
+#' @param fm A one-sided formula, possibly containing `s()` / `t2()` terms.
+#' @param data A data frame containing every variable referenced inside the
+#'   smooth terms (in the same row order as the design matrices that will be
+#'   built downstream).
+#' @return A list with components:
+#'   * `parametric_formula` — `fm` with smooth specials stripped.
+#'   * `smooth_labels`      — character vector of original smooth labels
+#'     (one per `s()` / `t2()` term).
+#'   * `Xf_list`            — list of `N × Ks_k` matrices (unpenalised cols),
+#'     one per smooth term.
+#'   * `Xr_list`            — **flat** list of `N × k_b` matrices, one per
+#'     penalty block. A single-penalty smooth contributes 1 entry; a multi-
+#'     penalty smooth contributes 2+. Each consumes one RE slot.
+#'   * `Xr_to_smooth`       — integer vector parallel to `Xr_list` mapping
+#'     each block to its parent smooth index.
+#'   * `Xr_slot_labels`     — character vector parallel to `Xr_list` giving
+#'     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).
+#'   * `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
+#'   `parametric_formula` is `fm` unchanged.
+#'
+#' @keywords internal
+#' @noRd
+parse_formula_smooths <- function(fm, data) {
+  empty <- list(
+    parametric_formula = fm,
+    smooth_labels      = character(0),
+    Xf_list            = list(),
+    Xr_list            = list(),
+    Xr_to_smooth       = integer(0),
+    Xr_slot_labels     = character(0),
+    smooth_specs       = list(),
+    smooth_re_objs     = list()
+  )
+  if (!has_smooth_terms(fm)) return(empty)
+  
+  check_and_install_packages("mgcv")
+  
+  trm <- stats::terms(fm, specials = c("s", "t2"))
+  vars <- attr(trm, "variables")
+  # +1 to skip the `list` call head when indexing into `vars`.
+  smooth_idx <- unlist(attr(trm, "specials")) + 1L
+  
+  smooth_calls  <- lapply(smooth_idx, function(i) vars[[i]])
+  smooth_labels <- vapply(
+    smooth_calls,
+    function(e) deparse(e, width.cutoff = 500L),
+    character(1)
+  )
+  smooth_labels <- unname(smooth_labels)
+  
+  Xf_list        <- vector("list", length(smooth_calls))
+  smooth_specs   <- vector("list", length(smooth_calls))
+  smooth_re_objs <- vector("list", length(smooth_calls))
+  Xr_list        <- list()   # flat across smooths × blocks
+  Xr_to_smooth   <- integer(0)
+  Xr_slot_labels <- character(0)
+  
+  for (k in seq_along(smooth_calls)) {
+    label <- smooth_labels[k]
+    
+    # Evaluate the s()/t2() call against the data frame. `s` and `t2` are
+    # bound in a small lookup env so the call resolves even when the
+    # caller's formula environment doesn't have mgcv attached.
+    fm_env  <- environment(fm)
+    if (is.null(fm_env)) fm_env <- baseenv()
+    eval_env <- new.env(parent = fm_env)
+    eval_env$s  <- mgcv::s
+    eval_env$t2 <- mgcv::t2
+    sm_spec <- eval(smooth_calls[[k]], envir = as.list(data), enclos = eval_env)
+    
+    # absorb.cons = TRUE removes the constant function from the basis so the
+    # smooth doesn't fight the intercept (matches brms).
+    # diagonal.penalty = TRUE re-parameterises so the penalty matrix is
+    # diagonal, which makes the random-effect prior N(0, sds^2 I) clean.
+    sm <- tryCatch(
+      mgcv::smoothCon(
+        sm_spec,
+        data            = as.data.frame(data),
+        absorb.cons     = TRUE,
+        diagonal.penalty = TRUE
+      )[[1]],
+      error = function(e) stop(sprintf(
+        "sccomp says: failed to build smooth basis for `%s`: %s",
+        label, conditionMessage(e)
+      ))
+    )
+    
+    re <- mgcv::smooth2random(sm, vnames = "", type = 2)
+    
+    Xf <- if (is.null(re$Xf)) matrix(0, nrow = nrow(data), ncol = 0) else re$Xf
+    attr(Xf, "s.label") <- NULL
+    if (ncol(Xf) > 0) {
+      colnames(Xf) <- sprintf("%s__lin%d", label, seq_len(ncol(Xf)))
+    }
+    
+    # Pull all penalty blocks (`re$rand$Xr`, `re$rand$Xr.0`, `re$rand$Xr.1`,
+    # ...). Each one becomes its own sccomp RE slot.
+    if (is.null(re$rand) || length(re$rand) == 0L) {
+      stop(sprintf(
+        "sccomp says: smooth `%s` produced no penalised basis (unsupported basis type?).",
+        label
+      ))
+    }
+    block_names    <- names(re$rand)
+    n_blocks_smooth <- length(block_names)
+    
+    for (b in seq_along(block_names)) {
+      Xr_b <- re$rand[[block_names[b]]]
+      attr(Xr_b, "s.label") <- NULL
+      # Single-block smooths keep the historical `<label>` slot label (and the
+      # matching `<label>___basis<jj>` column naming); multi-block smooths
+      # disambiguate with `<label>__b<b>`.
+      slot_label <- if (n_blocks_smooth == 1L) label
+                    else sprintf("%s__b%d", label, b)
+      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)
+    }
+    
+    Xf_list[[k]]        <- Xf
+    smooth_specs[[k]]   <- sm
+    smooth_re_objs[[k]] <- re
+  }
+  
+  list(
+    parametric_formula = strip_smooth_terms(fm),
+    smooth_labels      = smooth_labels,
+    Xf_list            = Xf_list,
+    Xr_list            = Xr_list,
+    Xr_to_smooth       = Xr_to_smooth,
+    Xr_slot_labels     = Xr_slot_labels,
+    smooth_specs       = smooth_specs,
+    smooth_re_objs     = smooth_re_objs
+  )
+}
+
+
+#' Evaluate the per-level base basis (`Xb`) of a factor-smooth at new data
+#'
+#' `mgcv::smooth2random.fs.interaction()` ignores `object$X` and rebuilds the
+#' basis from `object$Xb` (the per-level base basis) whenever `Xb` is non-NULL.
+#' To predict at new data we must therefore recompute `Xb` ourselves, which is
+#' the same logic `mgcv:::Predict.matrix.fs.interaction()` runs internally
+#' before its level-by-level expansion. We replicate just that initial step.
+#'
+#' @keywords internal
+#' @noRd
+compute_fs_Xb_at_newdata <- function(sm, newdata) {
+  newdata <- as.data.frame(newdata)
+  fac_dropped <- newdata
+  fac_dropped[[sm$fterm]] <- NULL
+  
+  obj <- sm
+  class(obj)         <- sm$base$bs
+  obj$rank           <- sm$base$rank
+  obj$null.space.dim <- sm$base$null.space.dim
+  obj$bs.dim         <- sm$base$bs.dim
+  obj$term           <- sm$base$term
+  mgcv::Predict.matrix(obj, fac_dropped) %*% sm$P
+}
+
+
+#' Evaluate a smooth basis at new covariate values
+#'
+#' Reproduces the (Xf, Xr_blocks) decomposition obtained at fit time, but on
+#' `newdata`. The approach is *clone-and-rerun*: we feed `mgcv::smooth2random`
+#' a shallow clone of the smooth object whose basis fields point at the new
+#' data, then re-invoke the same `smooth2random(..., type = 2)` decomposition
+#' that ran at fit time. Because `smooth2random` is purely a function of the
+#' smooth's basis matrices and penalty / level structure, this reproduces the
+#' fit-time `re$rand$Xr*` and `re$Xf` column-for-column when fed the original
+#' training rows.
+#'
+#' Supported smooth families:
+#'   * `tp.smooth`, `cr.smooth`, single-penalty smooths: clone-and-rerun on
+#'     `$X = PredictMat(sm, newdata)`.
+#'   * `fs.interaction` (`bs = "fs"`): factor smooths store a per-level base
+#'     basis in `$Xb` that `smooth2random` reads in preference to `$X`. We
+#'     recompute `Xb` at new rows via `compute_fs_Xb_at_newdata()` and also
+#'     refresh `$fac` so the level masks `fac == flev[j]` line up.
+#'
+#' Tensor-product `t2.smooth` is not yet supported at prediction time because
+#' `mgcv::PredictMat` applies the t2 constraint differently from `smoothCon`,
+#' producing a basis that does not match the fit-time decomposition.
+#'
+#' @param sm      An mgcv smooth object (the result of `smoothCon()[[1]]`).
+#' @param re_obj  The matching `smooth2random()` result (unused, kept for
+#'   API parity with prior versions).
+#' @param newdata A data frame containing every variable the smooth uses.
+#' @return A list with `Xf` (`N_new × Ks`) and `Xr_blocks`, a list of
+#'   `N_new × k_b` matrices (one per penalty block, in `re$rand` order).
+#'
+#' @keywords internal
+#' @noRd
+predict_smooth_at_newdata <- function(sm, re_obj, newdata) {
+  check_and_install_packages("mgcv")
+  if (inherits(sm, "t2.smooth")) {
+    stop(
+      "sccomp says: prediction with t2() smooths is not yet supported ",
+      "(constraint absorption mismatch between PredictMat() and smoothCon()). ",
+      "Use s(...) or s(..., bs = \"fs\") for now."
+    )
+  }
+  
+  newdata <- as.data.frame(newdata)
+  
+  sm_new <- sm
+  if (inherits(sm, "fs.interaction")) {
+    # fs smooths: smooth2random uses `$Xb` + `$base$S` rather than `$X` + `$S`,
+    # so we refresh `$Xb` (per-level base basis at new rows) and `$fac`.
+    fac_new <- newdata[[sm$fterm]]
+    if (is.null(fac_new)) {
+      stop(sprintf(
+        "sccomp says: factor `%s` is missing from new data (required by smooth `%s`).",
+        sm$fterm, sm$label
+      ))
+    }
+    sm_new$Xb  <- compute_fs_Xb_at_newdata(sm, newdata)
+    sm_new$fac <- factor(fac_new, levels = levels(sm$fac))
+  } else {
+    Xp <- mgcv::PredictMat(sm, data = newdata)
+    attr(Xp, "s.label") <- NULL
+    sm_new$X <- Xp
+  }
+  
+  re_new <- mgcv::smooth2random(sm_new, vnames = "", type = 2)
+  
+  Xf_new <- if (is.null(re_new$Xf)) {
+    matrix(0, nrow = nrow(newdata), ncol = 0)
+  } else {
+    as.matrix(re_new$Xf)
+  }
+  attr(Xf_new, "s.label") <- NULL
+  
+  Xr_blocks <- if (is.null(re_new$rand)) {
+    list()
+  } else {
+    lapply(re_new$rand, function(m) {
+      m <- as.matrix(m)
+      attr(m, "s.label") <- NULL
+      m
+    })
+  }
+  
+  list(Xf = Xf_new, Xr_blocks = Xr_blocks)
+}
+
+
+#' Wrap an Xr matrix into a sccomp random-effect slot
+#'
+#' Each smooth occupies exactly one slot with `n_factors = 1` and
+#' `n_groups = ncol(Xr)`: the slot's per-cell-group SD then plays the role of
+#' brms's smooth-wise `sds_*`. The slot structure (X / X_unseen / ncol / gfi /
+#' n_groups / n_factors) mirrors what `prepare_re_slot()` returns for an
+#' ordinary random-effect clause.
+#'
+#' @param Xr     A numeric matrix (N rows × K wiggly basis columns).
+#' @param label  Smooth label, used as the slot's single "factor" name.
+#' @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) {
+  stopifnot(is.matrix(Xr) || is.data.frame(Xr))
+  Xr <- as.matrix(Xr)
+  K  <- ncol(Xr)
+  
+  if (K == 0L) {
+    return(list(
+      X         = Xr,
+      X_unseen  = Xr[, integer(0), drop = FALSE],
+      ncol      = 0L,
+      gfi       = matrix(integer(0), nrow = 0L, ncol = 0L),
+      n_groups  = 0L,
+      n_factors = 0L
+    ))
+  }
+  
+  # One factor (= the smooth), K "groups" (= the basis columns).
+  # `gfi[1, j] = j` means basis-column j is at position j in X_random_effect_k.
+  gfi <- matrix(
+    seq_len(K),
+    nrow = 1L, ncol = K,
+    dimnames = list(label, colnames(Xr))
+  )
+  
+  list(
+    X         = Xr,
+    X_unseen  = Xr[, integer(0), drop = FALSE],
+    ncol      = K,
+    gfi       = gfi,
+    n_groups  = K,
+    n_factors = 1L
+  )
+}
+
+
+#' Build replicate-time smooth design pieces from fit-time metadata
+#'
+#' For every `s()` / `t2()` term recorded in `smooth_results`, evaluates the
+#' basis at `new_data_tail` via [predict_smooth_at_newdata()] and packages the
+#' two design contributions that `prepare_replicate_data()` needs:
+#'
+#' * the unpenalised (null-space) columns are renamed `<label>__lin*` and
+#'   appended to `parametric_X` (so they line up positionally with the
+#'   fit-time `beta` coefficients);
+#' * each penalised block becomes one random-effect slot, using the same
+#'   `<label>___basis*` (single-penalty) or `<label>__b<b>___basis*`
+#'   (multi-penalty) column naming as fit time, so matching against the
+#'   original `X_random_effect_*` matrices still succeeds.
+#'
+#' When `smooth_results` is empty (no smooth terms in the original formula)
+#' the function returns `parametric_X` unchanged and an empty slot list.
+#'
+#' @param parametric_X Numeric matrix; the parametric (non-smooth) design
+#'   matrix evaluated on the replicate rows. New unpenalised smooth columns
+#'   are appended on the right.
+#' @param new_data_tail Data frame of the replicate rows (must contain every
+#'   variable referenced by any smooth term).
+#' @param smooth_results A list with `smooth_specs`, `smooth_re_objs`, and
+#'   `smooth_labels` (parallel lists captured at fit time), or `NULL` /
+#'   length-0 to signal "no smooths".
+#' @return A list with:
+#'   * `new_X` — `parametric_X` with the unpenalised smooth columns appended.
+#'   * `smooth_replicate_slots` — flat list of slot lists (one per penalised
+#'     block, in fit-time order), each with `X`, `X_unseen`, `which`.
+#'
+#' @keywords internal
+#' @noRd
+build_smooth_replicate_design <- function(parametric_X,
+                                          new_data_tail,
+                                          smooth_results) {
+  smooth_specs   <- smooth_results$smooth_specs
+  smooth_re_objs <- smooth_results$smooth_re_objs
+  smooth_labels  <- smooth_results$smooth_labels
+  
+  if (length(smooth_specs) == 0L) {
+    return(list(new_X = parametric_X, smooth_replicate_slots = list()))
+  }
+  
+  # One basis evaluation per smooth term, carrying the label alongside the
+  # (Xf, Xr_blocks) pair so the downstream maps don't need parallel vectors.
+  smooth_evals <- purrr::pmap(
+    list(spec = smooth_specs, re = smooth_re_objs, label = smooth_labels),
+    function(spec, re, label) {
+      pred <- predict_smooth_at_newdata(spec, re, new_data_tail)
+      pred$label <- label
+      pred
+    }
+  )
+  
+  # Append unpenalised (Xf) columns in fit-time order.
+  new_X <- purrr::reduce(smooth_evals, .init = parametric_X, function(X, e) {
+    if (ncol(e$Xf) == 0L) return(X)
+    colnames(e$Xf) <- sprintf("%s__lin%d", e$label, seq_len(ncol(e$Xf)))
+    cbind(X, e$Xf)
+  })
+  
+  # One slot per penalised block, flattened across smooths.
+  smooth_replicate_slots <- smooth_evals |>
+    purrr::map(function(e) {
+      n_blocks <- length(e$Xr_blocks)
+      purrr::map2(e$Xr_blocks, seq_along(e$Xr_blocks), function(Xr_b, b) {
+        colnames(Xr_b) <- if (n_blocks == 1L) {
+          sprintf("%s___basis%02d", e$label, seq_len(ncol(Xr_b)))
+        } else {
+          sprintf("%s__b%d___basis%02d", e$label, b, seq_len(ncol(Xr_b)))
+        }
+        list(
+          X        = Xr_b,
+          X_unseen = Xr_b[, integer(0), drop = FALSE],
+          which    = seq_len(ncol(Xr_b)) |> as.array()
+        )
+      })
+    }) |>
+    unlist(recursive = FALSE)
+  
+  list(new_X = new_X, smooth_replicate_slots = smooth_replicate_slots)
+}
+
+#' Smooth metadata attached at fit time (R-only, not Stan data)
+#'
+#' @param x A fitted `sccomp_tbl` or its `model_input` list.
+#' @keywords internal
+#' @noRd
+get_smooth_results <- function(x) {
+  mi <- if (inherits(x, "sccomp_tbl")) attr(x, "model_input") else x
+  if (is.null(mi)) return(NULL)
+  attr(mi, "smooth_results")
+}
\ No newline at end of file
diff --git a/R/utilities.R b/R/utilities.R
index 1c0b633f..dbb03465 100755
--- a/R/utilities.R
+++ b/R/utilities.R
@@ -48,6 +48,25 @@ add_attr = function(var, attribute, name) {
   var
 }
 
+#' Ensure optional CRAN packages are installed (lazy dependency pattern)
+#'
+#' Used for features that need extra packages (e.g. `mgcv` for smooth terms)
+#' without listing them in `Imports:`. Mirrors the tidybulk pattern:
+#' `rlang::check_installed()` with a non-interactive install action.
+#'
+#' @param packages Character vector of package names.
+#' @keywords internal
+#' @noRd
+check_and_install_packages <- function(packages) {
+  rlang::check_installed(
+    pkg = packages,
+    action = function(...) {
+      install.packages(..., repos = getOption("repos"))
+    }
+  )
+}
+
+
 #' Subset results by model factor
 #'
 #' @param .data A sccomp results tibble
@@ -153,10 +172,16 @@ incorporate_parameters_into_sccomp_object = function(obj, parameters_to_load = N
 
 #' Formula parser
 #'
+#' Returns the names of data columns referenced by a formula. Random-effect
+#' grouping terms (anything containing `|`) are excluded. Smooth specials
+#' produced by mgcv (`s()` and `t2()`) are expanded to their underlying
+#' variable names so that, e.g., `~ type + s(pseudotime, k = 5)` returns
+#' `c("type", "pseudotime")` rather than the raw term label.
+#'
 #' @param fm A formula
 #'
 #' @importFrom stringr str_subset
-#' @importFrom magrittr extract2
+#' @importFrom purrr map
 #' @importFrom stats terms
 #'
 #' @return A character vector
@@ -165,14 +190,20 @@ incorporate_parameters_into_sccomp_object = function(obj, parameters_to_load = N
 #' @noRd
 parse_formula <- function(fm) {
   stopifnot("The formula must be of the kind \"~ factors\" " = attr(terms(fm), "response") == 0)
-  
-  as.character(attr(terms(fm), "variables")) |>
-    str_subset("\\|", negate = TRUE) %>%
-    
-    # Does not work the following
-    # |>
-    # extract2(-1)
-    .[-1]
+
+  # `variables` attribute lists each top-level term (interactions like `a:b`
+  # are not listed, which is the behaviour callers rely on). We drop the
+  # `list` call-head and any `(... | g)` random-effect clauses, then expand
+  # smooth wrappers (`s()`, `t2()`) to their underlying column names so
+  # downstream `select(all_of(...))` resolves against real columns.
+  attr(terms(fm), "variables") |>
+    as.character() |>
+    _[-1] |>
+    str_subset("\\|", negate = TRUE) |>
+    map(~ if (grepl("^(s|t2)\\(", .x)) all.vars(parse(text = .x)) else .x) |>
+    unlist() |>
+    as.character() |>
+    unique()
 }
 
 
@@ -991,28 +1022,55 @@ data_spread_to_model_input =
     
     if (length(.grouping_for_random_effect)==0 ) .grouping_for_random_effect = "random_effect"
     
+    # ----------------------------------------------------------------------
+    # Smooth terms: parse `s()` / `t2()` out of the composition formula so
+    # the parametric design `X` can be built from a smooth-free formula and
+    # the smooth basis pieces (`Xf`, `Xr_k`) are wired in explicitly below.
+    # Smooths in the variability formula are not supported yet — see TODO
+    # in `dev/splines_exploration/PLAN_sccomp_splines.md`.
+    # ----------------------------------------------------------------------
+    if (has_smooth_terms(formula_variability)) {
+      stop(
+        "sccomp says: smooth terms s()/t2() in `formula_variability` are not yet supported."
+      )
+    }
+    smooth_pieces = parse_formula_smooths(formula, .data_spread)
+    formula_parametric = smooth_pieces$parametric_formula
+    
     
     X  =
       
       .data_spread |>
       get_design_matrix(
         # Drop random intercept
-        formula |>
+        formula_parametric |>
+          strip_random_effect_terms() |>
           as.character() |>
-          str_remove_all("\\+ ?\\(.+\\|.+\\)") |>
           paste(collapse="") |>
           as.formula(),
         !!.sample,
         accept_NA_as_average_effect = accept_NA_as_average_effect
       )
     
+    # Append unpenalised (null-space) spline columns to the parametric X.
+    # These are estimated as ordinary `beta` rows by the existing Stan model
+    # — sccomp's `prior_mean_coefficients` already applies to non-intercept
+    # columns, which is what we want.
+    if (length(smooth_pieces$Xf_list) > 0) {
+      Xf_all = do.call(cbind, smooth_pieces$Xf_list)
+      if (ncol(Xf_all) > 0) {
+        rownames(Xf_all) = rownames(X)
+        X = cbind(X, Xf_all)
+      }
+    }
+    
     Xa  =
       .data_spread |>
       get_design_matrix(
         # Drop random intercept
         formula_variability |>
+          strip_random_effect_terms() |>
           as.character() |>
-          str_remove_all("\\+ ?\\(.+\\|.+\\)") |>
           paste(collapse="") |>
           as.formula() ,
         !!.sample,
@@ -1093,13 +1151,6 @@ data_spread_to_model_input =
           ~ get_random_effect_design3(.data_spread, .x, .y, !!.sample)
         ))
       
-      if (nrow(random_effect_grouping) > N_RE_SLOTS) {
-        stop(sprintf(
-          "sccomp says: at the moment sccomp supports up to %d random-effect groupings; the formula has %d. Combine related clauses or split the model.",
-          N_RE_SLOTS, nrow(random_effect_grouping)
-        ))
-      }
-      
       is_random_effect = 1
       n_random_eff = nrow(random_effect_grouping)
       
@@ -1121,6 +1172,38 @@ data_spread_to_model_input =
       re_slots         = list()
     }
     
+    # Each smooth term occupies one additional RE slot. `n_factors = 1` and
+    # `n_groups = ncol(Xr_k)` so the slot's per-cell-group SD plays the role
+    # of brms's `sds_*`. No Stan-side change is needed — the smooth slot is
+    # structurally identical to any other RE block from the Stan model's
+    # 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).
+    if (length(smooth_pieces$Xr_list) > 0) {
+      smooth_slots = map2(
+        smooth_pieces$Xr_list,
+        smooth_pieces$Xr_slot_labels,
+        build_smooth_slot
+      )
+      re_slots         = c(re_slots, smooth_slots)
+      is_random_effect = 1L
+      n_random_eff     = n_random_eff + length(smooth_slots)
+    }
+    
+    if (length(re_slots) > N_RE_SLOTS) {
+      stop(sprintf(
+        "sccomp says: the formula needs %d random-effect slot(s) (%d explicit clause(s) + %d smooth block(s) across %d smooth term(s); multi-penalty smooths such as t2() / bs=\"fs\" produce several blocks). The Stan model has %d. Reduce the number of smooths or merge RE clauses.",
+        length(re_slots),
+        length(re_slots) - length(smooth_pieces$Xr_list),
+        length(smooth_pieces$Xr_list),
+        length(smooth_pieces$smooth_labels),
+        N_RE_SLOTS
+      ))
+    }
+    
     # Pad to exactly N_RE_SLOTS so the Stan data list is uniformly shaped
     while (length(re_slots) < N_RE_SLOTS) re_slots = c(re_slots, list(empty_re_slot))
     
@@ -1286,6 +1369,17 @@ data_spread_to_model_input =
     # Default all grouping known (four RE slots; see glm_multi_beta_binomial_generate_data.stan)
     data_for_model$unknown_grouping = rep(0L, 4L)
     
+    # Smooth-term metadata is R-only (mgcv `smoothCon` / `smooth2random`
+    # objects, used by prediction / replicate helpers). It is NOT shipped to
+    # Stan, so we store it as an attribute on `model_input` (not a list element).
+    # Downstream code reads it via `get_smooth_results()`.
+    if (length(smooth_pieces$smooth_specs) > 0) {
+      attr(data_for_model, "smooth_results") <- list(
+        smooth_specs   = smooth_pieces$smooth_specs,
+        smooth_re_objs = smooth_pieces$smooth_re_objs,
+        smooth_labels  = smooth_pieces$smooth_labels
+      )
+    }
     
     # Return
     data_for_model
diff --git a/man/sccomp-package.Rd b/man/sccomp-package.Rd
index 0d214a20..7404392e 100644
--- a/man/sccomp-package.Rd
+++ b/man/sccomp-package.Rd
@@ -65,6 +65,7 @@ Useful links:
 
 Authors:
 \itemize{
+  \item Stefano Mangiola \email{stefano.mangiola@unimelb.edu.au}
   \item Alexandra J. Roth-Schulze
   \item Marie Trussart
   \item Enrique Zozaya-Valdés
diff --git a/man/sccomp_stan_models_cache_dir.Rd b/man/sccomp_stan_models_cache_dir.Rd
index 4738baeb..8a465a78 100644
--- a/man/sccomp_stan_models_cache_dir.Rd
+++ b/man/sccomp_stan_models_cache_dir.Rd
@@ -1,12 +1,8 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/utilities.R
-\docType{data}
 \name{sccomp_stan_models_cache_dir}
 \alias{sccomp_stan_models_cache_dir}
 \title{Default cache directory for Stan models}
-\format{
-An object of class \code{character} of length 1.
-}
 \usage{
 sccomp_stan_models_cache_dir
 }
@@ -36,4 +32,3 @@ sccomp_stan_models_cache_dir
 \seealso{
 \code{\link{sccomp_estimate}}, \code{\link{sccomp_replicate}}, \code{\link{clear_stan_model_cache}}
 }
-\keyword{datasets}
diff --git a/tests/testthat/test-smooths.R b/tests/testthat/test-smooths.R
new file mode 100644
index 00000000..f08c31ee
--- /dev/null
+++ b/tests/testthat/test-smooths.R
@@ -0,0 +1,418 @@
+## Tests for brms-style smooth terms (s(), t2()) in sccomp formulas.
+##
+## The unit tests at the top do not need cmdstan: they only exercise the
+## R-side machinery in `R/smooths.R` and the patched `parse_formula()`.
+## The end-to-end smoke test does need cmdstan and is gated by
+## `skip_cmdstan()`.
+
+library(dplyr)
+
+# -----------------------------------------------------------------------------
+# Pure-R unit tests
+# -----------------------------------------------------------------------------
+
+test_that("has_smooth_terms detects s() / t2() in a formula", {
+  expect_false(sccomp:::has_smooth_terms(~ 1))
+  expect_false(sccomp:::has_smooth_terms(~ type + age))
+  expect_true(sccomp:::has_smooth_terms(~ type + s(age)))
+  expect_true(sccomp:::has_smooth_terms(~ s(x, k = 5) + s(y, k = 3)))
+  expect_true(sccomp:::has_smooth_terms(~ t2(x, z)))
+})
+
+
+test_that("strip_smooth_terms removes s()/t2() and preserves other terms", {
+  strip   <- sccomp:::strip_smooth_terms
+  labels  <- function(fm) attr(stats::terms(fm), "term.labels")
+  has_int <- function(fm) attr(stats::terms(fm), "intercept") == 1L
+  
+  # ---- No-op cases --------------------------------------------------------
+  
+  # NULL passes through unchanged.
+  expect_null(strip(NULL))
+  
+  # Smooth-free formulas are unchanged in term content and intercept setting.
+  fm_plain <- ~ type + age + (1 | sample)
+  expect_equal(labels(strip(fm_plain)), labels(fm_plain))
+  expect_true(has_int(strip(fm_plain)))
+  
+  fm_noint <- ~ 0 + type + age
+  expect_equal(labels(strip(fm_noint)), labels(fm_noint))
+  expect_false(has_int(strip(fm_noint)))
+  
+  # `~ 1` and `~ 0` round-trip.
+  expect_true (has_int(strip(~ 1)))
+  expect_false(has_int(strip(~ 0)))
+  
+  # ---- Removes s() / t2(), keeps the rest ---------------------------------
+  
+  fm <- ~ type + s(pseudotime, k = 5)
+  expect_false(sccomp:::has_smooth_terms(strip(fm)))
+  expect_equal(labels(strip(fm)), "type")
+  expect_true(has_int(strip(fm)))
+  
+  # `+ 0` is preserved when smooths are stripped.
+  fm_noint_s <- ~ 0 + type + s(pseudotime)
+  expect_equal(labels(strip(fm_noint_s)), "type")
+  expect_false(has_int(strip(fm_noint_s)))
+  
+  # Smooth-only formula collapses to `~ 1` / `~ 0`.
+  expect_equal(labels(strip(~ s(x))), character(0))
+  expect_true (has_int(strip(~ s(x))))
+  expect_equal(labels(strip(~ 0 + s(x))), character(0))
+  expect_false(has_int(strip(~ 0 + s(x))))
+  
+  # `t2()` is removed just like `s()`.
+  expect_equal(labels(strip(~ type + t2(x, z, k = c(3, 4)))), "type")
+  
+  # Multiple smooths and a `by =` argument: all smooth calls vanish, parametric
+  # terms remain. Variables that only appeared inside `s(..., by = ...)` are
+  # gone because they were not parametric to begin with.
+  fm_many <- ~ type + age + s(x, k = 4) + s(y, by = condition) + t2(u, v)
+  expect_equal(labels(strip(fm_many)), c("type", "age"))
+  
+  # ---- Preserves random-effect clauses verbatim ---------------------------
+  
+  # Single bar term sits alongside the parametric terms.
+  fm_re <- ~ type + s(pseudotime) + (1 | sample)
+  expect_equal(labels(strip(fm_re)), c("type", "1 | sample"))
+  
+  # Multiple bar terms, including multi-LHS bars, are all kept.
+  fm_re_multi <- ~ 0 + type + s(x) + (1 | sample) + (age | dataset)
+  expect_setequal(
+    labels(strip(fm_re_multi)),
+    c("type", "1 | sample", "age | dataset")
+  )
+  expect_false(has_int(strip(fm_re_multi)))
+  
+  # ---- Preserves transforms and interactions ------------------------------
+  
+  # Function transforms and interactions survive intact.
+  fm_tx <- ~ type * age + I(age^2) + log(exposure) + s(pseudotime)
+  out_tx <- strip(fm_tx)
+  # `type * age` expands to type + age + type:age in term.labels; check the
+  # main effects and the interaction are all retained, and no smooth label
+  # remains.
+  expect_true (all(c("type", "age", "type:age", "I(age^2)", "log(exposure)")
+                   %in% labels(out_tx)))
+  expect_false(any(grepl("^s\\(",  labels(out_tx))))
+  expect_false(any(grepl("^t2\\(", labels(out_tx))))
+  
+  # Explicit interaction term with `:` is preserved.
+  expect_setequal(
+    labels(strip(~ a:b + s(x))),
+    "a:b"
+  )
+})
+
+
+test_that("parse_formula expands smooth specials to underlying variables", {
+  expect_equal(
+    sccomp:::parse_formula(~ type + s(pseudotime, k = 5)),
+    c("type", "pseudotime")
+  )
+  # `s(x, by = condition)` should expose both x and condition
+  expect_setequal(
+    sccomp:::parse_formula(~ s(x, by = condition)),
+    c("x", "condition")
+  )
+  # Random-effect bars are still excluded
+  expect_equal(
+    sccomp:::parse_formula(~ type + s(age, k = 4) + (1 | dataset)),
+    c("type", "age")
+  )
+})
+
+
+test_that("parse_formula_smooths produces Xf + Xr with brms shapes (tp, 1 block)", {
+  skip_if_not_installed("mgcv")
+  set.seed(1)
+  dat <- data.frame(x = sort(stats::runif(50, 0, 10)))
+  res <- sccomp:::parse_formula_smooths(~ s(x, bs = "tp", k = 8), dat)
+  
+  expect_length(res$Xf_list, 1L)
+  expect_length(res$Xr_list, 1L)              # 1 penalty -> 1 block
+  expect_equal(res$Xr_to_smooth, 1L)
+  expect_length(res$smooth_specs, 1L)
+  expect_length(res$smooth_re_objs, 1L)
+  expect_equal(res$smooth_labels, "s(x, bs = \"tp\", k = 8)")
+  
+  Xf <- res$Xf_list[[1]]
+  Xr <- res$Xr_list[[1]]
+  expect_equal(nrow(Xf), nrow(dat))
+  expect_equal(nrow(Xr), nrow(dat))
+  # tp + absorb.cons: 1 unpenalised column + (k-2) wiggly columns
+  expect_equal(ncol(Xf), 1L)
+  expect_equal(ncol(Xr), 8L - 2L)
+  expect_true(all(grepl("___basis", colnames(Xr))))
+})
+
+
+test_that("parse_formula_smooths flattens t2 / fs into multiple Xr blocks", {
+  skip_if_not_installed("mgcv")
+  set.seed(3)
+  
+  # t2: 3 penalties -> 3 random blocks
+  dat_t <- data.frame(
+    x = sort(stats::runif(40, 0, 10)),
+    y = sort(stats::runif(40, 0, 5))
+  )
+  res_t <- sccomp:::parse_formula_smooths(~ t2(x, y, k = c(4, 3)), dat_t)
+  expect_length(res_t$Xr_list, 3L)
+  expect_equal(res_t$Xr_to_smooth, c(1L, 1L, 1L))
+  # All blocks share the same parent smooth label
+  expect_length(res_t$smooth_labels, 1L)
+  # Multi-block: column names are disambiguated with __b<b>___basis
+  expect_true(any(grepl("__b1___basis", colnames(res_t$Xr_list[[1]]))))
+  expect_true(any(grepl("__b2___basis", colnames(res_t$Xr_list[[2]]))))
+  expect_true(any(grepl("__b3___basis", colnames(res_t$Xr_list[[3]]))))
+  
+  # fs: factor smooth with 3 levels and k=5 -> 3 penalties -> 3 blocks
+  dat_f <- data.frame(
+    x = seq(0, 6, length.out = 24),
+    z = factor(rep(c("a", "b", "c"), length.out = 24))
+  )
+  res_f <- sccomp:::parse_formula_smooths(~ s(x, z, bs = "fs", k = 5), dat_f)
+  expect_length(res_f$Xr_list, 3L)
+  expect_equal(res_f$Xr_to_smooth, c(1L, 1L, 1L))
+  # Sum of block widths matches the raw fs basis dimension (k * n_levels = 15)
+  expect_equal(sum(vapply(res_f$Xr_list, ncol, integer(1))), 5L * 3L)
+  # fs has no null-space block (everything is penalised)
+  expect_equal(ncol(res_f$Xf_list[[1]]), 0L)
+})
+
+
+test_that("parse_formula_smooths returns identity for smooth-free formulas", {
+  dat <- data.frame(type = letters[1:5], age = 1:5)
+  res <- sccomp:::parse_formula_smooths(~ type + age, dat)
+  expect_length(res$Xf_list, 0L)
+  expect_length(res$Xr_list, 0L)
+  expect_length(res$smooth_specs, 0L)
+  # Parametric formula identical to input
+  expect_equal(deparse(res$parametric_formula), deparse(~ type + age))
+})
+
+
+test_that("predict_smooth_at_newdata reproduces Xr at training data (tp)", {
+  skip_if_not_installed("mgcv")
+  set.seed(2)
+  dat <- data.frame(x = sort(stats::runif(40, 0, 10)))
+  res <- sccomp:::parse_formula_smooths(~ s(x, bs = "tp", k = 7), dat)
+  
+  sm     <- res$smooth_specs[[1]]
+  re_obj <- res$smooth_re_objs[[1]]
+  Xr_fit <- res$Xr_list[[1]]
+  Xf_fit <- res$Xf_list[[1]]
+  
+  pred <- sccomp:::predict_smooth_at_newdata(sm, re_obj, dat)
+  
+  expect_length(pred$Xr_blocks, 1L)
+  expect_equal(dim(pred$Xr_blocks[[1]]), dim(Xr_fit))
+  expect_equal(dim(pred$Xf), dim(Xf_fit))
+  
+  # Same columns, same rows: round-trip should be numerically identical
+  expect_equal(unname(pred$Xr_blocks[[1]]), unname(Xr_fit), tolerance = 1e-10)
+  expect_equal(unname(pred$Xf),            unname(Xf_fit), tolerance = 1e-10)
+})
+
+
+test_that("predict_smooth_at_newdata reproduces all Xr blocks at training data (fs)", {
+  skip_if_not_installed("mgcv")
+  set.seed(4)
+  
+  dat <- data.frame(
+    x = seq(0, 6, length.out = 24),
+    z = factor(rep(c("a", "b", "c"), length.out = 24))
+  )
+  res <- sccomp:::parse_formula_smooths(~ s(x, z, bs = "fs", k = 5), dat)
+  sm     <- res$smooth_specs[[1]]
+  re_obj <- res$smooth_re_objs[[1]]
+  
+  pred <- sccomp:::predict_smooth_at_newdata(sm, re_obj, dat)
+  expect_length(pred$Xr_blocks, length(res$Xr_list))
+  
+  for (b in seq_along(pred$Xr_blocks)) {
+    expect_equal(
+      unname(pred$Xr_blocks[[b]]),
+      unname(res$Xr_list[[b]]),
+      tolerance = 1e-10
+    )
+  }
+  expect_equal(unname(pred$Xf), unname(res$Xf_list[[1]]), tolerance = 1e-10)
+})
+
+
+test_that("predict_smooth_at_newdata errors cleanly on t2()", {
+  skip_if_not_installed("mgcv")
+  dat <- data.frame(
+    x = sort(stats::runif(30, 0, 10)),
+    y = sort(stats::runif(30, 0, 5))
+  )
+  res <- sccomp:::parse_formula_smooths(~ t2(x, y, k = c(4, 3)), dat)
+  expect_error(
+    sccomp:::predict_smooth_at_newdata(
+      res$smooth_specs[[1]], res$smooth_re_objs[[1]], dat
+    ),
+    regexp = "t2\\(\\) smooths is not yet supported"
+  )
+})
+
+
+test_that("build_smooth_slot produces a sccomp RE-slot structure", {
+  Xr <- matrix(stats::rnorm(20), nrow = 5, ncol = 4,
+               dimnames = list(NULL, paste0("c", 1:4)))
+  slot <- sccomp:::build_smooth_slot(Xr, "s(x)")
+  expect_equal(slot$ncol, 4L)
+  expect_equal(slot$n_factors, 1L)
+  expect_equal(slot$n_groups, 4L)
+  expect_equal(dim(slot$gfi), c(1L, 4L))
+  expect_equal(rownames(slot$gfi), "s(x)")
+  expect_equal(as.vector(slot$gfi), 1:4)
+})
+
+
+# -----------------------------------------------------------------------------
+# End-to-end smoke test (requires cmdstan)
+# -----------------------------------------------------------------------------
+
+test_that("sccomp_estimate works with a single s() smooth on a real dataset", {
+  skip_if_not_installed("mgcv")
+  skip_cmdstan()
+  data("seurat_obj")
+  
+  fit <- sccomp_estimate(
+    seurat_obj,
+    formula_composition = ~ s(continuous_covariate, k = 4),
+    formula_variability = ~ 1,
+    sample              = "sample",
+    cell_group          = "cell_group",
+    cores               = 1,
+    inference_method    = "pathfinder",
+    max_sampling_iterations = 200,
+    verbose = FALSE
+  )
+  
+  mi <- attr(fit, "model_input")
+  sr <- sccomp:::get_smooth_results(fit)
+  
+  # Smooth metadata lives on `model_input` (R-only, not a Stan list element).
+  expect_false(is.null(sr))
+  expect_equal(length(sr$smooth_specs), 1L)
+  expect_equal(length(sr$smooth_re_objs), 1L)
+  expect_equal(length(sr$smooth_labels), 1L)
+  # Not a Stan list element, and not duplicated on the fitted tibble.
+  expect_null(mi$smooth_results)
+  expect_null(attr(fit, "smooth_results"))
+  
+  # Model input slot 1 (no explicit RE clauses) carries the wiggly basis.
+  expect_gt(mi$ncol_X_random_eff[1], 0L)
+  expect_equal(mi$how_many_factors_in_random_design[1], 1L)
+  
+  # Fixed-effect design has the unpenalised null-space column too.
+  expect_true(any(grepl("__lin", colnames(mi$X))))
+})
+
+
+test_that("fs factor smooth fits and predicts end-to-end (multi-block)", {
+  skip_if_not_installed("mgcv")
+  skip_cmdstan()
+  data("counts_obj")
+  
+  # Attach a per-sample pseudotime and a 3-level tissue factor (both sample-
+  # level). counts_obj has 20 samples; we just assign tissues cyclically.
+  set.seed(7)
+  samples <- levels(counts_obj$sample)
+  covs <- tibble::tibble(
+    sample     = samples,
+    pseudotime = seq(0, 6, length.out = length(samples)) +
+                   rnorm(length(samples), 0, 0.05),
+    tissue     = factor(rep(c("blood", "lymph", "tumor"),
+                            length.out = length(samples)))
+  )
+  counts_pt <- counts_obj |> dplyr::left_join(covs, by = "sample")
+  
+  fit <- sccomp_estimate(
+    counts_pt,
+    formula_composition = ~ s(pseudotime, tissue, bs = "fs", k = 5),
+    formula_variability = ~ 1,
+    sample              = "sample",
+    cell_group          = "cell_group",
+    abundance           = "count",
+    cores               = 1,
+    inference_method    = "pathfinder",
+    max_sampling_iterations = 200,
+    verbose = FALSE
+  )
+  
+  mi <- attr(fit, "model_input")
+  sr <- sccomp:::get_smooth_results(fit)
+  
+  # One mgcv smooth -> 3 penalty blocks -> 3 RE slots occupied
+  expect_equal(length(sr$smooth_specs), 1L)
+  expect_equal(sum(mi$ncol_X_random_eff > 0L), 3L)
+  
+  # All wiggly columns landed in the RE slots; no `__lin` cols (fs has empty
+  # null-space)
+  expect_false(any(grepl("__lin", colnames(mi$X))))
+  
+  # Predict on a grid that spans pseudotime x tissue. Must succeed (this is
+  # the failure mode the multi-block refactor fixes).
+  grid <- expand.grid(
+    pseudotime = seq(min(covs$pseudotime), max(covs$pseudotime),
+                     length.out = 20),
+    tissue     = levels(covs$tissue)
+  ) |>
+    tibble::as_tibble() |>
+    dplyr::mutate(sample = sprintf("grid_%03d", dplyr::row_number()))
+  
+  pred <- sccomp_predict(fit, new_data = grid, number_of_draws = 100)
+  expect_true(all(c("pseudotime", "tissue", "proportion_mean") %in% names(pred)))
+  expect_equal(nrow(pred), nrow(grid) * length(unique(counts_obj$cell_group)))
+  # No NAs in the posterior summary
+  expect_true(all(is.finite(pred$proportion_mean)))
+})
+
+
+test_that("smooth slot count beyond N_RE_SLOTS errors cleanly", {
+  skip_if_not_installed("mgcv")
+  skip_cmdstan()
+  data("seurat_obj")
+  
+  # 4 smooths + 1 RE clause needs 5 slots; budget is 4. Use the same
+  # continuous_covariate four times under different `k` (parser doesn't
+  # de-duplicate intentionally — each call is a fresh basis).
+  expect_error(
+    sccomp_estimate(
+      seurat_obj,
+      formula_composition =
+        ~ s(continuous_covariate, k = 3) +
+          s(continuous_covariate, k = 4) +
+          s(continuous_covariate, k = 5) +
+          s(continuous_covariate, k = 6) +
+          (1 | group__),
+      formula_variability = ~ 1,
+      sample = "sample", cell_group = "cell_group",
+      cores = 1, inference_method = "pathfinder",
+      max_sampling_iterations = 50, verbose = FALSE
+    ),
+    regexp = "random-effect slot",
+    fixed  = FALSE
+  )
+})
+
+
+test_that("smooths in formula_variability are rejected", {
+  skip_cmdstan()
+  data("seurat_obj")
+  
+  expect_error(
+    sccomp_estimate(
+      seurat_obj,
+      formula_composition = ~ 1,
+      formula_variability = ~ s(continuous_covariate),
+      sample = "sample", cell_group = "cell_group",
+      cores = 1, inference_method = "pathfinder",
+      max_sampling_iterations = 50, verbose = FALSE
+    ),
+    regexp = "smooth terms .* `formula_variability`"
+  )
+})
diff --git a/vignettes/splines.Rmd b/vignettes/splines.Rmd
new file mode 100644
index 00000000..e99c5cc8
--- /dev/null
+++ b/vignettes/splines.Rmd
@@ -0,0 +1,499 @@
+---
+title: "Smooth (spline) terms in sccomp"
+author: "sccomp authors"
+date: "`r Sys.Date()`"
+package: sccomp
+output:
+  prettydoc::html_pretty:
+    theme: cayman
+    toc: yes
+    toc_depth: 2
+    number_sections: yes
+    fig_caption: yes
+    df_print: paged
+abstract: >
+  This vignette introduces `brms`-style smooth terms (`s()`, `t2()`) in
+  `sccomp`. Smooths let you model non-linear effects of a continuous
+  covariate — pseudotime, age, dose, time-since-treatment — on cell-type
+  proportions, while keeping the rest of the `sccomp` modelling pipeline
+  (compositional likelihood, variability modelling, outlier detection,
+  random effects) unchanged. We start from the bundled `counts_obj`
+  dataset, add a synthetic continuous covariate, fit a smooth, predict on
+  a dense grid, and visualise both the fitted curves and the underlying
+  basis decomposition.
+vignette: >
+  %\VignetteIndexEntry{Smooth (spline) terms in sccomp}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment  = "#>",
+  fig.width  = 7,
+  fig.height = 4.5,
+  dpi = 110
+)
+# Smooth vignette chunks need mgcv (Suggests) and CmdStan.
+run_smooth_vignette <-
+  requireNamespace("mgcv", quietly = TRUE) &&
+  instantiate::stan_cmdstan_exists()
+```
+
+# Why smooth terms?
+
+`sccomp` fits a (logit-)linear model on the unconstrained scale and then maps
+it to the simplex via softmax. If you only have categorical covariates this is
+all you need. But with a **continuous covariate** — pseudotime, age, dose,
+days since treatment — the relationship between covariate and proportion is
+rarely linear. A spline (`s()` from `mgcv`) replaces a single linear column
+with a small, flexible basis that lets the curve bend, while a penalty keeps
+the bend honest.
+
+`sccomp` adopts the same syntax as `brms` / `mgcv`:
+
+```{r show-formulae, eval=FALSE}
+~ type + s(pseudotime, k = 5)        # smooth of one variable
+~ type + t2(pseudotime, dose)        # tensor product of two variables
+```
+
+Under the hood every smooth is split into two pieces:
+
+* an unpenalised **null-space** block `Xf` (the linear "trend" part), and
+* a penalised **wiggly** block `Xr` (the bend-y part).
+
+`Xf` joins the design matrix as ordinary fixed-effect columns; `Xr` joins one
+of `sccomp`'s random-effect slots and is governed by a single smoothing
+parameter `sds`. This is exactly how `brms` handles smooths, and it means the
+Stan model does **not** change — only the R-side data assembly.
+
+# What is supported
+
+| Feature                                            | Status                                                |
+|----------------------------------------------------|-------------------------------------------------------|
+| `s(x, k = ...)` in `formula_composition`           | yes (1 slot)                                          |
+| `s(x, z, bs = "fs", k = ...)` (factor smooth)      | yes (3 slots: 1 wiggly + 2 per-level deviation blocks)|
+| Multiple smooth terms in the same formula          | yes (each penalty block consumes one RE slot)         |
+| Mix smooths with parametric terms (`~ type + s(x)`)| yes                                                   |
+| Mix smooths with explicit `(... | g)` random effects | yes (total ≤ 4 slots)                              |
+| Smooths in `formula_variability`                   | not yet                                               |
+| `t2(x, y, ...)` (tensor-product smooth)            | fits, but `sccomp_predict()` not yet supported        |
+| `s(x, by = factor)` (independent smooths per level) | works but consumes one slot **per level**            |
+
+The 4-slot budget comes from the existing `sccomp` Stan model. If your
+formula needs more than 4 random-effect slots in total (explicit `( | g)`
+clauses + smooth penalty blocks), `sccomp_estimate()` will stop with an
+informative error. Note that multi-penalty smooth types like `bs = "fs"`
+emit several blocks; in mgcv / brms these correspond to several `sds_*`
+parameters.
+
+# Worked example
+
+```{r load-libs, message=FALSE, warning=FALSE}
+library(dplyr)
+library(tidyr)
+library(ggplot2)
+library(sccomp)
+data("counts_obj")
+```
+
+## 1. Synthetic continuous covariate
+
+`counts_obj` has no continuous covariate, so we attach one — a *pseudotime*
+between 0 and 6 — that varies across samples but is constant within a
+sample (a covariate must always be sample-level in `sccomp`).
+
+```{r make-pseudotime}
+set.seed(1)
+samples <- levels(counts_obj$sample)
+pt_map <- tibble(
+  sample     = samples,
+  pseudotime = seq(0, 6, length.out = length(samples)) +
+                 rnorm(length(samples), 0, 0.05)
+)
+
+counts_pt <- counts_obj |> left_join(pt_map, by = "sample")
+counts_pt |> distinct(sample, type, pseudotime) |> head()
+```
+
+## 2. Fit `~ type + s(pseudotime, k = 5)`
+
+```{r fit-smooth, message=FALSE, warning=FALSE, eval = run_smooth_vignette}
+fit <- counts_pt |>
+  sccomp_estimate(
+    formula_composition = ~ type + s(pseudotime, k = 5),
+    formula_variability = ~ 1,
+    sample              = "sample",
+    cell_group          = "cell_group",
+    abundance           = "count",
+    cores               = 1,
+    verbose             = FALSE,
+    max_sampling_iterations = 2000
+  )
+```
+
+Look at the design-matrix column names that `sccomp_estimate()` prints. You
+should see one new fixed column
+
+```
+s(pseudotime, k = 5)__lin1
+```
+
+and three penalised basis columns living in an RE slot
+
+```
+s(pseudotime, k = 5)___basis01
+s(pseudotime, k = 5)___basis02
+s(pseudotime, k = 5)___basis03
+```
+
+The `__lin*` columns are the null-space block (`Xf`); the `___basis*` columns
+are the wiggly block (`Xr`). The single smoothing parameter `sds` for this
+smooth controls how much the `___basis*` coefficients are pulled toward
+zero — exactly the role of `sds_*` in `brms`.
+
+## 3. Tabular output
+
+`sccomp_test()` works unchanged. It returns one row per (cell_group,
+parameter), and smooth basis columns appear like any other coefficient:
+
+```{r test-table, eval = run_smooth_vignette}
+fit |>
+  sccomp_test() |>
+  select(cell_group, parameter, c_lower, c_effect, c_upper, c_FDR) |>
+  head(10)
+```
+
+`c_effect` on the `___basis*` rows is rarely interpreted directly — the
+basis is internal — but the linear `__lin1` row is a clean "is there an
+overall monotone trend in pseudotime for this cell group?" test, and the
+basis-level `c_FDR` values let you spot non-linear effects.
+
+## 4. Predict on a dense grid
+
+To draw the fitted curve, predict at a regular grid of `pseudotime` values
+(remember to also supply any other covariates the formula needs — here,
+`type`):
+
+```{r predict-grid, message=FALSE, warning=FALSE, eval = run_smooth_vignette}
+grid <- expand.grid(
+  pseudotime = seq(min(pt_map$pseudotime), max(pt_map$pseudotime), length.out = 60),
+  type       = levels(counts_obj$type)
+) |>
+  as_tibble() |>
+  mutate(sample = sprintf("grid_%03d", row_number()))
+
+pred <- fit |>
+  sccomp_predict(new_data = grid, number_of_draws = 200)
+
+head(pred)
+```
+
+`sccomp_predict()` already returns the `new_data` columns (`pseudotime`,
+`type`, `sample`) alongside `proportion_mean` and the 95% credible interval
+`proportion_lower` / `proportion_upper`.
+
+## 5. Plot: fitted curves over the raw proportions
+
+We show the six cell groups with the largest pseudotime range in the
+posterior mean.
+
+```{r plot-curves, fig.width=8, fig.height=5, eval = run_smooth_vignette}
+top_cgs <- pred |>
+  group_by(cell_group) |>
+  summarise(rng = max(proportion_mean) - min(proportion_mean)) |>
+  arrange(desc(rng)) |>
+  slice_head(n = 6) |>
+  pull(cell_group)
+
+raw <- counts_pt |>
+  group_by(sample) |>
+  mutate(prop = count / sum(count)) |>
+  ungroup() |>
+  filter(cell_group %in% top_cgs)
+
+pred |>
+  filter(cell_group %in% top_cgs) |>
+  ggplot(aes(pseudotime, proportion_mean, colour = type, fill = type)) +
+  geom_ribbon(aes(ymin = proportion_lower, ymax = proportion_upper),
+              alpha = 0.15, linewidth = 0) +
+  geom_line(linewidth = 0.7) +
+  geom_point(data = raw, aes(y = prop), alpha = 0.5, size = 1.2) +
+  facet_wrap(~ cell_group, scales = "free_y", ncol = 3) +
+  labs(x = "pseudotime", y = "proportion",
+       title = "Fitted compositional curves vs pseudotime") +
+  theme_minimal(base_size = 10)
+```
+
+The blue/red curves are posterior means, ribbons are 95% credible intervals,
+points are observed sample proportions. Note the wiggle — e.g. CD8 1 peaks
+around pseudotime ≈ 4 — which a purely linear `~ type + pseudotime` model
+could not capture.
+
+## 6. Same model with a smaller basis (`k = 3`)
+
+The same parametric + smooth layout works with a lower `k`. Fewer basis
+functions mean a smoother curve and less risk of over-fitting when sample
+size is modest — here `k = 3` yields one unpenalised linear column and a
+single wiggly basis column (instead of three with `k = 5`).
+
+```{r fit-smooth-k3, message=FALSE, warning=FALSE, eval = run_smooth_vignette}
+fit_k3 <- counts_pt |>
+  sccomp_estimate(
+    formula_composition = ~ type + s(pseudotime, k = 3),
+    formula_variability = ~ 1,
+    sample              = "sample",
+    cell_group          = "cell_group",
+    abundance           = "count",
+    cores               = 1,
+    verbose             = FALSE,
+    max_sampling_iterations = 2000
+  )
+```
+
+The printed design matrix should list `(Intercept)`, `typecancer`, one
+`__lin` column, and a single `___basis01` column for the smooth.
+
+```{r predict-grid-k3, message=FALSE, warning=FALSE, eval = run_smooth_vignette}
+pred_k3 <- fit_k3 |>
+  sccomp_predict(new_data = grid, number_of_draws = 200)
+```
+
+```{r plot-curves-k3, fig.width=8, fig.height=5, eval = run_smooth_vignette}
+pred_k3 |>
+  filter(cell_group %in% top_cgs) |>
+  ggplot(aes(pseudotime, proportion_mean, colour = type, fill = type)) +
+  geom_ribbon(aes(ymin = proportion_lower, ymax = proportion_upper),
+              alpha = 0.15, linewidth = 0) +
+  geom_line(linewidth = 0.7) +
+  geom_point(data = raw, aes(y = prop, colour = type), alpha = 0.5, size = 1.2) +
+  facet_wrap(~ cell_group, scales = "free_y", ncol = 3) +
+  labs(x = "pseudotime", y = "proportion",
+       title = "Fitted curves: ~ type + s(pseudotime, k = 3)",
+       subtitle = "One wiggly basis column — smoother than k = 5") +
+  theme_minimal(base_size = 10)
+```
+
+Compared with the `k = 5` figure above, the credible bands are often a
+little tighter and the lines less wiggly: the model has less capacity to
+bend, which is appropriate when you want a gentle pseudotime trend on top
+of the `type` effect.
+
+# What `s()` actually decomposes into
+
+The smooth metadata used to fit the model is stored on the result, so we
+can re-evaluate the basis at any new pseudotime value. This is exactly what
+`sccomp_predict()` does internally; we can also do it by hand to *see* the
+basis.
+
+```{r basis-decomp, fig.width=8, fig.height=6, eval = run_smooth_vignette}
+smooth_specs   <- attr(fit, "smooth_specs")
+smooth_re_objs <- attr(fit, "smooth_re_objs")
+
+grid_basis <- tibble(
+  pseudotime = seq(min(pt_map$pseudotime), max(pt_map$pseudotime),
+                   length.out = 200),
+  type = levels(counts_obj$type)[1]  # any constant works; smooth ignores it
+)
+
+basis <- sccomp:::predict_smooth_at_newdata(
+  smooth_specs[[1]], smooth_re_objs[[1]], grid_basis
+)
+# Plain s() has one penalised block — for multi-penalty smooths
+# (e.g. bs = "fs", t2(...)) `Xr_blocks` would contain several matrices.
+Xr_single <- basis$Xr_blocks[[1]]
+
+basis_long <- bind_rows(
+  as.data.frame(basis$Xf) |>
+    setNames(sprintf("lin%d", seq_len(ncol(basis$Xf)))) |>
+    mutate(pseudotime = grid_basis$pseudotime) |>
+    pivot_longer(-pseudotime, names_to = "column", values_to = "value") |>
+    mutate(block = "Xf  (unpenalised, joins fixed effects)"),
+  as.data.frame(Xr_single) |>
+    setNames(sprintf("basis%d", seq_len(ncol(Xr_single)))) |>
+    mutate(pseudotime = grid_basis$pseudotime) |>
+    pivot_longer(-pseudotime, names_to = "column", values_to = "value") |>
+    mutate(block = "Xr  (penalised, joins one RE slot, shared sds)")
+)
+
+ggplot(basis_long, aes(pseudotime, value, colour = column)) +
+  geom_line(linewidth = 0.7) +
+  facet_wrap(~ block, scales = "free_y", ncol = 1) +
+  labs(x = "pseudotime", y = "basis value",
+       title = "Decomposition of s(pseudotime, k = 5)",
+       subtitle = "Linear null-space (Xf) is a fixed effect; wiggly basis (Xr) is one RE slot") +
+  theme_minimal(base_size = 10)
+```
+
+The top panel is the **null-space**: for a default thin-plate `s()` this is
+just the linear term in `pseudotime`. The bottom panel is the **wiggly
+basis**: each of these columns gets one coefficient per cell group, and a
+*single* smoothing standard deviation `sds` per cell group controls how
+much wiggle is allowed. Increasing `k` would add more wiggly basis columns;
+the penalty would then keep the extra capacity from over-fitting.
+
+# Factor smooths: per-group curves with shared structure
+
+A common ask is: "fit a smooth per tissue (or per condition, per donor, ...)
+but partial-pool across groups". The `mgcv` shorthand for this is
+`s(x, factor, bs = "fs")` — a *factor smooth*. It builds one common
+wiggly curve and a per-level deviation, governed jointly by a small set of
+smoothing parameters (three, in mgcv's standard reparameterisation: the
+common wiggly part and two per-level deviation blocks). In `sccomp` each
+of these blocks occupies its own RE slot.
+
+## A worked example
+
+We extend `counts_obj` with a synthetic `tissue` factor (3 levels) on top of
+the existing `pseudotime`:
+
+```{r make-tissue, message=FALSE}
+set.seed(11)
+samples <- levels(counts_obj$sample)
+covs <- tibble::tibble(
+  sample     = samples,
+  pseudotime = seq(0, 6, length.out = length(samples)) +
+                 rnorm(length(samples), 0, 0.05),
+  tissue     = factor(rep(c("blood", "lymph", "tumor"),
+                          length.out = length(samples)))
+)
+counts_tissue <- counts_obj |> dplyr::left_join(covs, by = "sample")
+counts_tissue |> dplyr::distinct(sample, type, pseudotime, tissue) |> head()
+```
+
+Fit `~ s(pseudotime, tissue, bs = "fs", k = 5)`:
+
+```{r fs-fit, message=FALSE, warning=FALSE, eval = run_smooth_vignette}
+fit_fs <- counts_tissue |>
+  sccomp_estimate(
+    formula_composition = ~ s(pseudotime, tissue, bs = "fs", k = 5),
+    formula_variability = ~ 1,
+    sample              = "sample",
+    cell_group          = "cell_group",
+    abundance           = "count",
+    cores               = 1,
+    verbose             = FALSE,
+    max_sampling_iterations = 2000
+  )
+
+# This one smooth occupies 3 RE slots (one per mgcv penalty block);
+# you can confirm via:
+attr(fit_fs, "model_input")$ncol_X_random_eff
+```
+
+The three non-zero entries correspond to the wiggly main block (9 columns
+for `k = 5` × 3 tissues) and two smaller per-level deviation blocks.
+
+## Predict and plot per-tissue curves
+
+`sccomp_predict()` takes a grid that crosses `pseudotime` with `tissue`:
+
+```{r fs-predict, message=FALSE, warning=FALSE, eval = run_smooth_vignette}
+grid_fs <- expand.grid(
+  pseudotime = seq(min(covs$pseudotime), max(covs$pseudotime),
+                   length.out = 60),
+  tissue     = levels(covs$tissue)
+) |>
+  as_tibble() |>
+  mutate(sample = sprintf("fs_grid_%03d", row_number()))
+
+pred_fs <- fit_fs |>
+  sccomp_predict(new_data = grid_fs, number_of_draws = 200)
+
+head(pred_fs)
+```
+
+Pick the most pseudotime-responsive cell groups and plot one curve per tissue:
+
+```{r fs-plot, fig.width=8.5, fig.height=5.5, eval = run_smooth_vignette}
+top_cgs_fs <- pred_fs |>
+  group_by(cell_group) |>
+  summarise(rng = max(proportion_mean) - min(proportion_mean)) |>
+  arrange(desc(rng)) |>
+  slice_head(n = 6) |>
+  pull(cell_group)
+
+raw_fs <- counts_tissue |>
+  group_by(sample) |>
+  mutate(prop = count / sum(count)) |>
+  ungroup() |>
+  filter(cell_group %in% top_cgs_fs)
+
+pred_fs |>
+  filter(cell_group %in% top_cgs_fs) |>
+  ggplot(aes(pseudotime, proportion_mean,
+             colour = tissue, fill = tissue)) +
+  geom_ribbon(aes(ymin = proportion_lower, ymax = proportion_upper),
+              alpha = 0.15, linewidth = 0) +
+  geom_line(linewidth = 0.7) +
+  geom_point(data = raw_fs, aes(y = prop, colour = tissue),
+             alpha = 0.5, size = 1.2) +
+  facet_wrap(~ cell_group, scales = "free_y", ncol = 3) +
+  labs(x = "pseudotime", y = "proportion",
+       title = "Factor smooth s(pseudotime, tissue, bs = 'fs', k = 5)",
+       subtitle = "Per-tissue curves sharing one factor-smooth structure") +
+  theme_minimal(base_size = 10)
+```
+
+The three coloured curves are the same smooth applied to different tissues,
+allowed to deviate at each pseudotime point. Where the deviations are small
+(e.g. CD8 1) the curves track each other; where they are large
+(e.g. CD8 2) one tissue diverges. Whether the deviations are "real" is
+expressed by the credible-interval ribbons — overlap means the data
+don't support a difference.
+
+## fs vs `s(x, by = tissue)`
+
+`s(x, by = tissue)` is an alternative spelling that *also* gives a per-tissue
+curve, but it treats each tissue's smooth as independent: each level gets
+its own smoothing parameter, and the curves can have totally unrelated
+shape. With L tissue levels this consumes L sccomp slots (one per level).
+Use `bs = "fs"` when you want a single shared smoothness scale and partial
+pooling across levels; use `by =` when the tissues are clearly heterogeneous
+and a shared smoothness assumption would be misleading.
+
+
+# Multiple smooths and mixing with random effects
+
+`sccomp_estimate()` is happy with several smooths, and you can also keep
+explicit `(... | g)` random-effect clauses:
+
+```{r multi-smooth-formula, eval=FALSE}
+~ type +
+  s(pseudotime, k = 5) +
+  s(age,        k = 4) +
+  (1 | donor)
+```
+
+Each `s()` / `t2()` term consumes one of the four available RE slots.
+The remaining slots are shared with any explicit `( | g)` clauses, so the
+total must stay ≤ 4. If you exceed it, `sccomp_estimate()` will stop with
+a message naming the smooth and the slot budget.
+
+# Caveats and current limitations
+
+* **No smooths in `formula_variability`** yet. `sccomp_estimate()` errors
+  out if you try.
+* **`t2()` tensor smooths** fit fine but `sccomp_predict()` is not yet
+  supported for them: `mgcv::PredictMat()` applies the `t2` constraint
+  differently from `smoothCon()`, producing a basis that doesn't line up
+  with the fitted coefficients. The predict call errors out cleanly. Use
+  `s(x, z, bs = "fs")` (shown above) when you want a smooth indexed by a
+  factor, or stack univariate `s(x) + s(y)` smooths for additive
+  multivariate effects.
+* **`s(x, by = factor)`** decomposes into one independent smooth per factor
+  level, each costing a separate RE slot. With many levels you will run
+  out of the 4-slot budget. For partial-pooling across factor levels use
+  `bs = "fs"` (shown above) — it fits in 3 slots regardless of the number
+  of levels.
+* The smooth columns are part of `sccomp_test()`'s output. The basis-level
+  `c_effect` / `c_FDR` are diagnostics, not biological effect sizes;
+  interpret the *curve* (predicted proportions) rather than individual
+  basis coefficients.
+
+# Session info
+
+```{r session-info, echo=FALSE}
+sessionInfo()
+```