1 Why smooth terms?

Continuous factors in linear models are tricky, becausegenerally you need to make assumptions about their linearity of higher polinomial trend. A spline replaces a single linear column with a small, flexible basis that lets the curve bend, while a penalty keeps the bend conservative

sccomp adopts the same syntax as brms / mgcv:

~ type + s(pseudotime, k = 5)

Under the hood every smooth is split into two pieces:

2 What is supported

Feature Notes
s(x, k = ...) in formula_composition
s(x, z, bs = "fs", k = ...) (factor smooth)
Multiple smooth terms in the same formula
Mix smooths with parametric terms (~ type + s(x))
Mix smooths with explicit (... | g) random effects
s(x, by = factor) (independent smooths per level) Max 2 smooth terms allowed 
library(dplyr)
library(tidyr)
library(ggplot2)
library(sccomp)
data("counts_obj")

3 Continuous covariate

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()

3.1 1. Fit ~ type + s(pseudotime, k = 5)

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
  )

The output will include one new linear parameter

s(pseudotime, k = 5)__lin1

and three penalised basis

s(pseudotime, k = 5)___basis01
s(pseudotime, k = 5)___basis02
s(pseudotime, k = 5)___basis03
fit |>
  sccomp_test() |>
  select(cell_group, parameter, c_lower, c_effect, c_upper, c_FDR)

These parameters cannot be interpreted directly, on the contrary of other parameter in the linear model.

3.2 2. Visualise the curve

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):

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.

We show the six cell groups with the largest pseudotime range in the posterior mean.

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) +
  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.

3.3 3. Test differences along the curve

The individual s(pseudotime, k = 5)__* parameters from sccomp_test() cannot be interpreted one-by-one. The scientific question is usually “does composition differ between two points (or intervals) of pseudotime?”. sccomp_test_smooth() answers exactly that: it predicts the mean composition (logit scale) at the two positions, forms the per-draw contrast mu(to) - mu(from), and runs it through the same probability-of-null / FDR machinery as sccomp_test() — so c_effect, c_pH0 and c_FDR mean the same thing here as everywhere else in sccomp.

Because the linear predictor is additive, any covariate held constant between the two positions (e.g. type) cancels out of the contrast, so the test isolates the pseudotime effect.

fit |>
  sccomp_test_smooth(
    smooth = "pseudotime",
    from   = 1,
    to     = 5
  ) |>
  arrange(c_FDR) |>
  select(cell_group, from, to, c_effect, c_lower, c_upper, c_FDR)

A positive c_effect means the cell group is (on the logit scale) more abundant at pseudotime 5 than at pseudotime 1; c_FDR controls the false-discovery rate across cell groups.

Because the curve can wiggle, comparing two single points can be sensitive to the exact positions. Use comparison = "average" to contrast the average of the curve over two intervals instead (each sampled at resolution points):

fit |>
  sccomp_test_smooth(
    smooth     = "pseudotime",
    from       = c(0, 2),   # early interval
    to         = c(4, 6),   # late interval
    comparison = "average"
  ) |>
  arrange(c_FDR) |>
  select(cell_group, from, to, c_effect, c_lower, c_upper, c_FDR)

For factor smooths (next section) the grouping factor selects which curve is tested — pass it via at, e.g. at = list(tissue = "tumor").

3.4 4. Let’s use 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.

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
  )
pred_k3 <- fit_k3 |>
  sccomp_predict(new_data = grid, number_of_draws = 200)
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) +
  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.

3.5 5. Under the hood: 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.

smooth_meta    <- sccomp:::get_smooth_results(fit)
smooth_specs   <- smooth_meta$smooth_specs
smooth_re_objs <- smooth_meta$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 linear trend, 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.

4 Hierarchical 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 shorthand for this (also used in mgcv) is s(x, factor, bs = "fs"). It builds one common wiggly curve and a per-level deviation, governed jointly by a small set of smoothing parameters.

4.1 1. Fit ~ s(pseudotime, tissue, bs = “fs”, k = 5)

We extend counts_obj with a synthetic tissue factor (3 levels) on top of the existing pseudotime:

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_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.

4.2 2. Predict and plot per-tissue curves

sccomp_predict() takes a grid that crosses pseudotime with tissue:

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:

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.

4.3 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. 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.

5 Multiple smooths and mixing with random effects

sccomp can model up to two splines together with random effects.

~ type +
  s(pseudotime, k = 5) +
  s(age,        k = 4) +
  (1 | donor)

6 Session info

#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB              LC_COLLATE=C              
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] tidyr_1.3.2       forcats_1.0.1     ggplot2_4.0.3     sccomp_2.5.3     
#> [5] instantiate_0.2.4 dplyr_1.2.1      
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.1            farver_2.1.2               
#>  [3] S7_0.2.2                    fastmap_1.2.0              
#>  [5] SingleCellExperiment_1.35.2 tensorA_0.36.2.1           
#>  [7] dotCall64_1.2               digest_0.6.39              
#>  [9] lifecycle_1.0.5             SeuratObject_5.4.0         
#> [11] processx_3.9.0              magrittr_2.0.5             
#> [13] posterior_1.7.0             compiler_4.6.1             
#> [15] rlang_1.3.0                 sass_0.4.10                
#> [17] tools_4.6.1                 yaml_2.3.12                
#> [19] knitr_1.52                  S4Arrays_1.13.2            
#> [21] sp_2.2-4                    DelayedArray_0.39.8        
#> [23] RColorBrewer_1.1-3          abind_1.4-8                
#> [25] withr_3.0.3                 purrr_1.2.2                
#> [27] BiocGenerics_0.59.13        grid_4.6.1                 
#> [29] stats4_4.6.1                future_1.82.0              
#> [31] progressr_1.1.0             globals_0.20.0             
#> [33] scales_1.4.0                dichromat_2.0-1            
#> [35] SummarizedExperiment_1.43.0 cli_3.6.6                  
#> [37] rmarkdown_2.32              crayon_1.5.3               
#> [39] generics_0.1.4              otel_0.3.0                 
#> [41] future.apply_1.21.0         tzdb_0.5.0                 
#> [43] cachem_1.1.0                stringr_1.6.0              
#> [45] parallel_4.6.1              XVector_0.53.0             
#> [47] matrixStats_1.5.0           vctrs_0.7.3                
#> [49] Matrix_1.7-6                jsonlite_2.0.0             
#> [51] callr_3.8.0                 IRanges_2.47.5             
#> [53] hms_1.1.4                   patchwork_1.3.2            
#> [55] S4Vectors_0.51.10           ggrepel_0.9.8              
#> [57] listenv_1.1.0               jquerylib_0.1.4            
#> [59] spam_2.11-4                 parallelly_1.50.0          
#> [61] ggside_0.4.1                glue_1.8.1                 
#> [63] codetools_0.2-20            ps_1.9.3                   
#> [65] distributional_0.9.0        stringi_1.8.9              
#> [67] gtable_0.3.6                GenomicRanges_1.65.4       
#> [69] tibble_3.3.1                pillar_1.11.1              
#> [71] htmltools_0.5.9             Seqinfo_1.3.2              
#> [73] R6_2.6.1                    evaluate_1.0.5             
#> [75] lattice_0.23-1              Biobase_2.73.2             
#> [77] readr_2.2.0                 backports_1.5.1            
#> [79] bslib_0.12.0                Rcpp_1.1.2                 
#> [81] SparseArray_1.13.4          checkmate_2.3.4            
#> [83] xfun_0.62                   fs_2.1.0                   
#> [85] MatrixGenerics_1.25.0       prettydoc_0.4.1            
#> [87] pkgconfig_2.0.3