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:

an unpenalised linear trend block Xf, and
a penalised wiggly block Xr (the non linear part).

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

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()
#> # A tibble: 6 × 3
#>   sample       type   pseudotime
#>   <chr>        <fct>       <dbl>
#> 1 10x_6K       benign    -0.0313
#> 2 10x_8K       benign     0.325 
#> 3 GSE115189    benign     0.590 
#> 4 SCP345_580   benign     1.03  
#> 5 SCP345_860   benign     1.28  
#> 6 SCP424_pbmc1 benign     1.54

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.

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

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

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.

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()
#> # A tibble: 6 × 4
#>   sample       type   pseudotime tissue
#>   <chr>        <fct>       <dbl> <fct> 
#> 1 10x_6K       benign    -0.0296 blood 
#> 2 10x_8K       benign     0.317  lymph 
#> 3 GSE115189    benign     0.556  tumor 
#> 4 SCP345_580   benign     0.879  blood 
#> 5 SCP345_860   benign     1.32   lymph 
#> 6 SCP424_pbmc1 benign     1.53   tumor

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.

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.

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.

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)

Session info

#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [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: Etc/UTC
#> 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.0     
#> [5] instantiate_0.2.3 dplyr_1.2.1       prettydoc_0.4.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                 utf8_1.2.6                 
#> [19] yaml_2.3.12                 knitr_1.52                 
#> [21] S4Arrays_1.13.2             sp_2.2-3                   
#> [23] DelayedArray_0.39.8         RColorBrewer_1.1-3         
#> [25] abind_1.4-8                 withr_3.0.3                
#> [27] purrr_1.2.2                 BiocGenerics_0.59.12       
#> [29] sys_3.4.3                   grid_4.6.1                 
#> [31] stats4_4.6.1                future_1.76.0              
#> [33] progressr_1.0.0             globals_0.19.1             
#> [35] scales_1.4.0                SummarizedExperiment_1.43.0
#> [37] cli_3.6.6                   rmarkdown_2.32             
#> [39] crayon_1.5.3                generics_0.1.4             
#> [41] otel_0.2.0                  future.apply_1.20.2        
#> [43] tzdb_0.5.0                  cachem_1.1.0               
#> [45] stringr_1.6.0               parallel_4.6.1             
#> [47] XVector_0.53.0              matrixStats_1.5.0          
#> [49] vctrs_0.7.3                 Matrix_1.7-6               
#> [51] jsonlite_2.0.0              callr_3.8.0                
#> [53] IRanges_2.47.5              hms_1.1.4                  
#> [55] patchwork_1.3.2             S4Vectors_0.51.10          
#> [57] ggrepel_0.9.8               listenv_1.1.0              
#> [59] maketools_1.3.2             jquerylib_0.1.4            
#> [61] spam_2.11-4                 parallelly_1.48.0          
#> [63] ggside_0.4.1                glue_1.8.1                 
#> [65] codetools_0.2-20            ps_1.9.3                   
#> [67] distributional_0.9.0        stringi_1.8.9              
#> [69] gtable_0.3.6                GenomicRanges_1.65.4       
#> [71] tibble_3.3.1                pillar_1.11.1              
#> [73] htmltools_0.5.9             Seqinfo_1.3.2              
#> [75] R6_2.6.1                    evaluate_1.0.5             
#> [77] lattice_0.23-1              Biobase_2.73.2             
#> [79] readr_2.2.0                 backports_1.5.1            
#> [81] bslib_0.12.0                Rcpp_1.1.2                 
#> [83] SparseArray_1.13.4          checkmate_2.3.4            
#> [85] xfun_0.61                   fs_2.1.0                   
#> [87] MatrixGenerics_1.25.0       buildtools_1.0.0           
#> [89] pkgconfig_2.0.3

Smooth (spline) terms in sccomp