From a805110563bb5644c016fccf2806d1b855abc893 Mon Sep 17 00:00:00 2001 From: A Wokaty Date: Wed, 29 Oct 2025 11:13:54 -0400 Subject: [PATCH 01/11] bump x.y.z version to even y prior to creation of RELEASE_3_22 branch --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 78729895..8d9f1571 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: sccomp Type: Package Title: Differential Composition and Variability Analysis for Single-Cell Data -Version: 2.1.18 +Version: 2.2.0 Date: 2024-01-15 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. From ea505684b96ffc0a1909b3e08b6be4213a164a45 Mon Sep 17 00:00:00 2001 From: A Wokaty Date: Wed, 29 Oct 2025 11:13:54 -0400 Subject: [PATCH 02/11] bump x.y.z version to odd y following creation of RELEASE_3_22 branch --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 8d9f1571..8613958b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: sccomp Type: Package Title: Differential Composition and Variability Analysis for Single-Cell Data -Version: 2.2.0 +Version: 2.3.0 Date: 2024-01-15 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. From e0b24e07d53c5c30c2057f8600955117c8c4cde3 Mon Sep 17 00:00:00 2001 From: A Wokaty Date: Tue, 28 Apr 2026 08:57:29 -0400 Subject: [PATCH 03/11] bump x.y.z version to even y prior to creation of RELEASE_3_23 branch --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 8613958b..f2313dd7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: sccomp Type: Package Title: Differential Composition and Variability Analysis for Single-Cell Data -Version: 2.3.0 +Version: 2.4.0 Date: 2024-01-15 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. From 67110133c81134c0eea020eab7b916279ce23187 Mon Sep 17 00:00:00 2001 From: A Wokaty Date: Tue, 28 Apr 2026 08:57:30 -0400 Subject: [PATCH 04/11] bump x.y.z version to odd y following creation of RELEASE_3_23 branch --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index f2313dd7..f06a9890 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: sccomp Type: Package Title: Differential Composition and Variability Analysis for Single-Cell Data -Version: 2.4.0 +Version: 2.5.0 Date: 2024-01-15 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. From 90cd05060670cf78b02e820b79e4e53143e09ac9 Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Thu, 23 Jul 2026 10:00:59 +0930 Subject: [PATCH 05/11] Add sccomp_test_smooth function for hypothesis testing along smooth terms --- NAMESPACE | 5 + R/sccomp_test_smooth.R | 314 ++++++++++++++++++++++++++++++ man/sccomp_test_smooth.Rd | 124 ++++++++++++ tests/testthat/test-smooth-test.R | 152 +++++++++++++++ vignettes/splines.Rmd | 54 ++++- 5 files changed, 647 insertions(+), 2 deletions(-) create mode 100644 R/sccomp_test_smooth.R create mode 100644 man/sccomp_test_smooth.Rd create mode 100644 tests/testthat/test-smooth-test.R diff --git a/NAMESPACE b/NAMESPACE index c5866835..bd119351 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -13,6 +13,7 @@ S3method(sccomp_remove_outliers,sccomp_tbl) S3method(sccomp_remove_unwanted_effects,sccomp_tbl) S3method(sccomp_replicate,sccomp_tbl) S3method(sccomp_test,sccomp_tbl) +S3method(sccomp_test_smooth,sccomp_tbl) S3method(simulate_data,tbl) export(clear_draw_files) export(clear_stan_model_cache) @@ -32,6 +33,7 @@ export(sccomp_replicate) export(sccomp_scatterplot) export(sccomp_stan_models_cache_dir) export(sccomp_test) +export(sccomp_test_smooth) export(sccomp_theme) export(simulate_data) import(dplyr) @@ -68,6 +70,8 @@ importFrom(dplyr,rowwise) importFrom(dplyr,select) importFrom(dplyr,slice) importFrom(dplyr,summarise) +importFrom(dplyr,tibble) +importFrom(dplyr,ungroup) importFrom(dplyr,where) importFrom(dplyr,with_groups) importFrom(forcats,fct_inorder) @@ -148,6 +152,7 @@ importFrom(rlang,quo_is_symbolic) importFrom(rlang,quo_name) importFrom(rlang,quo_squash) importFrom(rlang,set_names) +importFrom(rlang,sym) importFrom(scales,trans_new) importFrom(stats,C) importFrom(stats,as.formula) diff --git a/R/sccomp_test_smooth.R b/R/sccomp_test_smooth.R new file mode 100644 index 00000000..8948d09c --- /dev/null +++ b/R/sccomp_test_smooth.R @@ -0,0 +1,314 @@ +# Hypothesis testing for smooth (spline) terms. +# +# `sccomp_test_smooth()` tests whether the mean composition differs between two +# positions (or two intervals) of a continuous covariate entered as a smooth +# `s()` term. It is a thin wrapper that: +# 1. predicts the *mean* linear predictor (`mu_unconstrained`, the logit-scale +# mean — no beta-binomial overdispersion) at the requested covariate +# positions via the existing replicate-data / generated-quantities pass, +# 2. forms the per-draw contrast `mu(to) - mu(from)` for each cell group, and +# 3. feeds those contrast draws into the same `draws_to_statistics()` engine +# used by `sccomp_test()`, so the returned `c_pH0` / `c_FDR` mean exactly +# what they mean elsewhere in sccomp. +# +# The contrast lives on the logit (unconstrained) scale, which is the scale +# `test_composition_above_logit_fold_change` is calibrated for. + +#' Test differences along a smooth term +#' +#' @description +#' Tests whether the mean cell-group composition differs between two positions +#' (or two intervals) of a continuous covariate modelled with a smooth `s()` +#' term. The test is a contrast of the model's mean prediction on the logit +#' scale, evaluated at the two positions and summarised with the same +#' probability-of-null (`pH0`) and false-discovery-rate (`FDR`) machinery as +#' [sccomp_test()]. +#' +#' @details +#' Because the smooth curve can be non-linear, "difference across an interval" +#' is ambiguous. Two comparison modes are provided: +#' \itemize{ +#' \item \code{"endpoints"} (default): \code{from} and \code{to} are single +#' values; the contrast is \eqn{\mu(\code{to}) - \mu(\code{from})}. +#' \item \code{"average"}: \code{from} and \code{to} are length-2 intervals +#' \code{c(lo, hi)}; the contrast is the average of the curve over the +#' \code{to} interval minus the average over the \code{from} interval +#' (each interval sampled at \code{resolution} points). +#' } +#' +#' Any covariate not being varied cancels out of the contrast because the +#' logit-scale predictor is additive; it therefore does not matter what value +#' those covariates take, as long as it is held constant. For factor smooths +#' (\code{s(x, g, bs = "fs")}) or by-factor smooths (\code{s(x, by = g)}) the +#' grouping factor selects \emph{which} curve is tested; set it via \code{at}, +#' e.g. \code{at = list(tissue = "tumor")}. +#' +#' The test is about the \strong{mean} composition: it uses the expected +#' linear predictor (\code{mu_unconstrained}), not the overdispersed +#' beta-binomial realisation. It therefore answers "does the expected +#' composition differ between these covariate positions?", not "would an +#' individual sample differ?". +#' +#' @param fit The result of [sccomp_estimate()] with a smooth term in +#' \code{formula_composition}. +#' @param smooth Character. The continuous covariate inside the smooth to test +#' (e.g. \code{"pseudotime"}). Can be \code{NULL} when the model has exactly +#' one smooth term over one continuous covariate, in which case it is +#' auto-detected. +#' @param from Numeric. The reference position. A single value for +#' \code{comparison = "endpoints"}, or a length-2 interval \code{c(lo, hi)} +#' for \code{comparison = "average"}. +#' @param to Numeric. The comparison position, same shape as \code{from}. +#' @param comparison One of \code{"endpoints"} or \code{"average"}. See details. +#' @param at Optional named list fixing the value of other covariates (e.g. +#' the grouping factor of a factor smooth). Covariates not supplied are held +#' at their first observed value (which cancels for non-grouping covariates). +#' @param resolution Integer. Number of grid points used to approximate each +#' interval average when \code{comparison = "average"}. +#' @param test_composition_above_logit_fold_change Positive numeric. Effect +#' threshold for the hypothesis test, on the logit scale — identical meaning +#' to the argument of [sccomp_test()]. +#' @param percent_false_positive Numeric in (0, 100). Used for the credible +#' interval width, as in [sccomp_test()]. +#' @param number_of_draws Integer. Number of posterior draws used for the +#' prediction pass. +#' @param mcmc_seed Integer. Seed for the generated-quantities pass. +#' @param robust Logical. Currently unused placeholder for API parity with +#' [sccomp_predict()]; the effect is always the posterior mean of the +#' contrast. +#' +#' @return A tibble with one row per cell group: +#' \itemize{ +#' \item \code{cell_group} — the cell group tested. +#' \item \code{smooth} — the continuous covariate tested. +#' \item \code{from}, \code{to} — the compared positions (as labels). +#' \item \code{c_lower}, \code{c_effect}, \code{c_upper} — 95% CI and posterior +#' mean of the logit-scale contrast \eqn{\mu(\code{to}) - \mu(\code{from})}. +#' \item \code{c_pH0} — probability the effect is within the null region. +#' \item \code{c_FDR} — false-discovery rate across cell groups. +#' } +#' +#' @examples +#' \donttest{ +#' if (instantiate::stan_cmdstan_exists()) { +#' data("counts_obj") +#' # add a continuous covariate +#' counts_obj$pseudotime <- as.numeric(factor(counts_obj$sample)) +#' +#' fit <- sccomp_estimate( +#' counts_obj, +#' ~ s(pseudotime, k = 4), ~ 1, "sample", "cell_group", "count", +#' cores = 1 +#' ) +#' +#' fit |> sccomp_test_smooth(from = 2, to = 8) +#' } +#' } +#' +#' @export +sccomp_test_smooth <- function(fit, + smooth = NULL, + from, + to, + comparison = c("endpoints", "average"), + at = NULL, + resolution = 20L, + test_composition_above_logit_fold_change = 0.1, + percent_false_positive = 5, + number_of_draws = 500, + mcmc_seed = sample_seed(), + robust = FALSE) { + check_and_install_cmdstanr() + UseMethod("sccomp_test_smooth", fit) +} + +#' @importFrom dplyr tibble bind_rows left_join group_by summarise ungroup arrange mutate select distinct pull +#' @importFrom tidyr pivot_wider nest unnest +#' @importFrom rlang enquo quo_name sym +#' @export +sccomp_test_smooth.sccomp_tbl <- function(fit, + smooth = NULL, + from, + to, + comparison = c("endpoints", "average"), + at = NULL, + resolution = 20L, + test_composition_above_logit_fold_change = 0.1, + percent_false_positive = 5, + number_of_draws = 500, + mcmc_seed = sample_seed(), + robust = FALSE) { + + # Define the variables as NULL to avoid CRAN NOTES + M <- NULL + .chain <- NULL + .iteration <- NULL + .draw <- NULL + .value <- NULL + .smooth_group <- NULL + parameter <- NULL + + comparison <- match.arg(comparison) + .sample <- attr(fit, ".sample") + .cell_group <- attr(fit, ".cell_group") + model_input <- attr(fit, "model_input") + smooth_specs <- get_smooth_results(fit)$smooth_specs + + if (is.null(smooth_specs) || length(smooth_specs) == 0) + stop("sccomp says: this fit has no smooth (s()/t2()) terms; use sccomp_test() for parametric contrasts.") + + # ---- Resolve which continuous covariate to test ------------------------- + # For each smooth, the continuous variable(s) are `term` minus the grouping + # factor `fterm` (present only for factor smooths). + smooth_continuous <- lapply(smooth_specs, function(sm) setdiff(sm$term, sm$fterm)) + all_continuous <- unique(unlist(smooth_continuous)) + + if (is.null(smooth)) { + if (length(all_continuous) != 1) + stop(sprintf( + "sccomp says: the model has several smooth covariates (%s); specify one via `smooth = `.", + paste(all_continuous, collapse = ", ") + )) + smooth <- all_continuous + } else if (!smooth %in% all_continuous) { + stop(sprintf( + "sccomp says: `%s` is not a continuous smooth covariate in this model (available: %s).", + smooth, paste(all_continuous, collapse = ", ") + )) + } + + # ---- Validate from/to shapes against comparison mode -------------------- + if (comparison == "endpoints") { + if (length(from) != 1L || length(to) != 1L) + stop("sccomp says: for comparison = 'endpoints', `from` and `to` must be single values.") + from_grid <- from + to_grid <- to + from_label <- as.character(from) + to_label <- as.character(to) + } else { + if (length(from) != 2L || length(to) != 2L) + stop("sccomp says: for comparison = 'average', `from` and `to` must be length-2 intervals c(lo, hi).") + from_grid <- seq(from[1], from[2], length.out = resolution) + to_grid <- seq(to[1], to[2], length.out = resolution) + from_label <- sprintf("[%s, %s]", from[1], from[2]) + to_label <- sprintf("[%s, %s]", to[1], to[2]) + } + + # ---- Build new_data at the requested positions -------------------------- + new_data <- build_smooth_test_newdata( + fit, smooth, from_grid, to_grid, at, .sample + ) + + # ---- Predict the mean linear predictor (logit scale) -------------------- + rng <- replicate_data( + fit, + formula_composition = NULL, # use the fitted model + formula_variability = ~ 1, + new_data = new_data, + number_of_draws = number_of_draws, + mcmc_seed = mcmc_seed + ) + + sample_col <- quo_name(.sample) + sample_names <- new_data |> pull(!!.sample) + + # Per-draw logit-scale mean (mu_unconstrained), keeping chain/iteration so + # the diagnostics in draws_to_statistics remain meaningful. + draws <- + rng |> + draws_to_tibble_x_y("mu_unconstrained", "M", "N") |> + # Map N -> sample name -> group ("from"/"to") + left_join( + tibble( + N = seq_along(sample_names), + !!sample_col := sample_names + ), + by = "N" + ) |> + left_join( + new_data |> select(!!.sample, .smooth_group), + by = sample_col + ) |> + # Map M -> cell group name + left_join( + tibble( + M = seq_len(ncol(model_input$y)), + !!quo_name(.cell_group) := colnames(model_input$y) + ), + by = "M" + ) + + # ---- Form the contrast per (cell_group, draw) --------------------------- + # Average over grid points within each group (a no-op for endpoints, where + # each group has one point), then difference to - from. + contrast_draws <- + draws |> + group_by(!!.cell_group, M, .chain, .iteration, .draw, .smooth_group) |> + summarise(.value = mean(.value), .groups = "drop") |> + pivot_wider(names_from = .smooth_group, values_from = .value) |> + mutate(.value = to - from) |> + mutate(parameter = sprintf("%s: %s -> %s", smooth, from_label, to_label)) |> + select(!!.cell_group, M, parameter, .chain, .iteration, .draw, .value) + + # ---- Reuse the standard pH0 / FDR engine -------------------------------- + result <- + contrast_draws |> + draws_to_statistics( + percent_false_positive / 100, + test_composition_above_logit_fold_change, + !!.cell_group, + "c_" + ) |> + select(-M, -parameter) |> + mutate(smooth = smooth, from = from_label, to = to_label) |> + select(!!.cell_group, smooth, from, to, everything()) + + result +} + + +#' Build the sample-wise new_data for a smooth test +#' +#' Two groups of rows ("from" and "to"), each covering the requested grid of +#' the continuous covariate. All other covariates in the composition formula +#' are held at a constant value (from `at`, else the first observed value), so +#' they cancel in the contrast; a grouping factor supplied via `at` selects +#' which curve of a factor / by-factor smooth is tested. +#' +#' @keywords internal +#' @noRd +#' @importFrom dplyr tibble bind_rows distinct mutate row_number all_of +build_smooth_test_newdata <- function(fit, smooth, from_grid, to_grid, at, .sample) { + + sample_col <- quo_name(.sample) + count_data <- attr(fit, "count_data") + formula_composition <- attr(fit, "formula_composition") + + covariates <- parse_formula(formula_composition) + other_covs <- setdiff(covariates, smooth) + + # Defaults for non-varied covariates: `at` override, else first observed. + ref_row <- count_data |> distinct(!!.sample, .keep_all = TRUE) |> head(1) + cov_value <- function(v) { + if (!is.null(at) && v %in% names(at)) return(at[[v]]) + if (v %in% colnames(count_data)) return(ref_row[[v]]) + stop(sprintf("sccomp says: covariate `%s` needed by the model is missing from the data and not supplied via `at`.", v)) + } + other_values <- lapply(other_covs, cov_value) + names(other_values) <- other_covs + + make_group <- function(x_grid, group_label) { + df <- tibble(!!smooth := x_grid, .smooth_group = group_label) + for (v in other_covs) df[[v]] <- other_values[[v]] + df + } + + new_data <- + bind_rows( + make_group(from_grid, "from"), + make_group(to_grid, "to") + ) |> + mutate(!!sample_col := sprintf("smoothtest_%03d", row_number())) + + new_data +} diff --git a/man/sccomp_test_smooth.Rd b/man/sccomp_test_smooth.Rd new file mode 100644 index 00000000..9ee96be9 --- /dev/null +++ b/man/sccomp_test_smooth.Rd @@ -0,0 +1,124 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/sccomp_test_smooth.R +\name{sccomp_test_smooth} +\alias{sccomp_test_smooth} +\title{Test differences along a smooth term} +\usage{ +sccomp_test_smooth( + fit, + smooth = NULL, + from, + to, + comparison = c("endpoints", "average"), + at = NULL, + resolution = 20L, + test_composition_above_logit_fold_change = 0.1, + percent_false_positive = 5, + number_of_draws = 500, + mcmc_seed = sample_seed(), + robust = FALSE +) +} +\arguments{ +\item{fit}{The result of \code{\link[=sccomp_estimate]{sccomp_estimate()}} with a smooth term in +\code{formula_composition}.} + +\item{smooth}{Character. The continuous covariate inside the smooth to test +(e.g. \code{"pseudotime"}). Can be \code{NULL} when the model has exactly +one smooth term over one continuous covariate, in which case it is +auto-detected.} + +\item{from}{Numeric. The reference position. A single value for +\code{comparison = "endpoints"}, or a length-2 interval \code{c(lo, hi)} +for \code{comparison = "average"}.} + +\item{to}{Numeric. The comparison position, same shape as \code{from}.} + +\item{comparison}{One of \code{"endpoints"} or \code{"average"}. See details.} + +\item{at}{Optional named list fixing the value of other covariates (e.g. +the grouping factor of a factor smooth). Covariates not supplied are held +at their first observed value (which cancels for non-grouping covariates).} + +\item{resolution}{Integer. Number of grid points used to approximate each +interval average when \code{comparison = "average"}.} + +\item{test_composition_above_logit_fold_change}{Positive numeric. Effect +threshold for the hypothesis test, on the logit scale — identical meaning +to the argument of \code{\link[=sccomp_test]{sccomp_test()}}.} + +\item{percent_false_positive}{Numeric in (0, 100). Used for the credible +interval width, as in \code{\link[=sccomp_test]{sccomp_test()}}.} + +\item{number_of_draws}{Integer. Number of posterior draws used for the +prediction pass.} + +\item{mcmc_seed}{Integer. Seed for the generated-quantities pass.} + +\item{robust}{Logical. Currently unused placeholder for API parity with +\code{\link[=sccomp_predict]{sccomp_predict()}}; the effect is always the posterior mean of the +contrast.} +} +\value{ +A tibble with one row per cell group: +\itemize{ +\item \code{cell_group} — the cell group tested. +\item \code{smooth} — the continuous covariate tested. +\item \code{from}, \code{to} — the compared positions (as labels). +\item \code{c_lower}, \code{c_effect}, \code{c_upper} — 95\% CI and posterior +mean of the logit-scale contrast \eqn{\mu(\code{to}) - \mu(\code{from})}. +\item \code{c_pH0} — probability the effect is within the null region. +\item \code{c_FDR} — false-discovery rate across cell groups. +} +} +\description{ +Tests whether the mean cell-group composition differs between two positions +(or two intervals) of a continuous covariate modelled with a smooth \code{s()} +term. The test is a contrast of the model's mean prediction on the logit +scale, evaluated at the two positions and summarised with the same +probability-of-null (\code{pH0}) and false-discovery-rate (\code{FDR}) machinery as +\code{\link[=sccomp_test]{sccomp_test()}}. +} +\details{ +Because the smooth curve can be non-linear, "difference across an interval" +is ambiguous. Two comparison modes are provided: +\itemize{ +\item \code{"endpoints"} (default): \code{from} and \code{to} are single +values; the contrast is \eqn{\mu(\code{to}) - \mu(\code{from})}. +\item \code{"average"}: \code{from} and \code{to} are length-2 intervals +\code{c(lo, hi)}; the contrast is the average of the curve over the +\code{to} interval minus the average over the \code{from} interval +(each interval sampled at \code{resolution} points). +} + +Any covariate not being varied cancels out of the contrast because the +logit-scale predictor is additive; it therefore does not matter what value +those covariates take, as long as it is held constant. For factor smooths +(\code{s(x, g, bs = "fs")}) or by-factor smooths (\code{s(x, by = g)}) the +grouping factor selects \emph{which} curve is tested; set it via \code{at}, +e.g. \code{at = list(tissue = "tumor")}. + +The test is about the \strong{mean} composition: it uses the expected +linear predictor (\code{mu_unconstrained}), not the overdispersed +beta-binomial realisation. It therefore answers "does the expected +composition differ between these covariate positions?", not "would an +individual sample differ?". +} +\examples{ +\donttest{ + if (instantiate::stan_cmdstan_exists()) { + data("counts_obj") + # add a continuous covariate + counts_obj$pseudotime <- as.numeric(factor(counts_obj$sample)) + + fit <- sccomp_estimate( + counts_obj, + ~ s(pseudotime, k = 4), ~ 1, "sample", "cell_group", "count", + cores = 1 + ) + + fit |> sccomp_test_smooth(from = 2, to = 8) + } +} + +} diff --git a/tests/testthat/test-smooth-test.R b/tests/testthat/test-smooth-test.R new file mode 100644 index 00000000..24e2b661 --- /dev/null +++ b/tests/testthat/test-smooth-test.R @@ -0,0 +1,152 @@ +test_that("sccomp_test_smooth() errors on a fit without smooth terms", { + skip_if_not_installed("mgcv") + skip_cmdstan() + data("counts_obj") + + fit <- sccomp_estimate( + counts_obj, + formula_composition = ~ type, + formula_variability = ~ 1, + sample = "sample", + cell_group = "cell_group", + abundance = "count", + cores = 1, + inference_method = "pathfinder", + max_sampling_iterations = 200, + verbose = FALSE + ) + + expect_error( + sccomp_test_smooth(fit, from = 1, to = 2), + "no smooth" + ) +}) + +test_that("sccomp_test_smooth() endpoints contrast returns the standard columns", { + skip_if_not_installed("mgcv") + skip_cmdstan() + data("counts_obj") + + samples <- levels(counts_obj$sample) + covs <- tibble::tibble( + sample = samples, + pseudotime = seq(0, 6, length.out = length(samples)) + ) + counts_pt <- counts_obj |> dplyr::left_join(covs, by = "sample") + + fit <- sccomp_estimate( + counts_pt, + formula_composition = ~ s(pseudotime, k = 4), + formula_variability = ~ 1, + sample = "sample", + cell_group = "cell_group", + abundance = "count", + cores = 1, + inference_method = "pathfinder", + max_sampling_iterations = 200, + verbose = FALSE + ) + + res <- fit |> + sccomp_test_smooth(from = 1, to = 5, number_of_draws = 100) + + n_groups <- counts_pt |> dplyr::distinct(cell_group) |> nrow() + expect_equal(nrow(res), n_groups) + + expect_true(all( + c("cell_group", "smooth", "from", "to", + "c_lower", "c_effect", "c_upper", "c_pH0", "c_FDR") %in% colnames(res) + )) + + # auto-detected smooth variable, correct labels + expect_equal(unique(res$smooth), "pseudotime") + expect_equal(unique(res$from), "1") + expect_equal(unique(res$to), "5") + + # probabilities and FDR are in [0, 1] + expect_true(all(res$c_pH0 >= 0 & res$c_pH0 <= 1)) + expect_true(all(res$c_FDR >= 0 & res$c_FDR <= 1)) + # CI brackets the point effect + expect_true(all(res$c_lower <= res$c_effect & res$c_effect <= res$c_upper)) +}) + +test_that("sccomp_test_smooth() average mode requires length-2 intervals", { + skip_if_not_installed("mgcv") + skip_cmdstan() + data("counts_obj") + + samples <- levels(counts_obj$sample) + covs <- tibble::tibble( + sample = samples, + pseudotime = seq(0, 6, length.out = length(samples)) + ) + counts_pt <- counts_obj |> dplyr::left_join(covs, by = "sample") + + fit <- sccomp_estimate( + counts_pt, + formula_composition = ~ s(pseudotime, k = 4), + formula_variability = ~ 1, + sample = "sample", + cell_group = "cell_group", + abundance = "count", + cores = 1, + inference_method = "pathfinder", + max_sampling_iterations = 200, + verbose = FALSE + ) + + # wrong shape for average + expect_error( + sccomp_test_smooth(fit, from = 1, to = 5, comparison = "average"), + "length-2" + ) + + res <- fit |> + sccomp_test_smooth( + from = c(0, 2), to = c(4, 6), + comparison = "average", resolution = 6, number_of_draws = 100 + ) + + n_groups <- counts_pt |> dplyr::distinct(cell_group) |> nrow() + expect_equal(nrow(res), n_groups) + expect_equal(unique(res$from), "[0, 2]") + expect_equal(unique(res$to), "[4, 6]") +}) + +test_that("sccomp_test_smooth() selects a curve via `at` for factor smooths", { + skip_if_not_installed("mgcv") + skip_cmdstan() + data("counts_obj") + + samples <- levels(counts_obj$sample) + covs <- tibble::tibble( + sample = samples, + pseudotime = seq(0, 6, length.out = length(samples)), + 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 + ) + + res <- fit |> + sccomp_test_smooth( + smooth = "pseudotime", from = 1, to = 5, + at = list(tissue = "tumor"), number_of_draws = 100 + ) + + n_groups <- counts_pt |> dplyr::distinct(cell_group) |> nrow() + expect_equal(nrow(res), n_groups) + expect_equal(unique(res$smooth), "pseudotime") +}) diff --git a/vignettes/splines.Rmd b/vignettes/splines.Rmd index 80dd6371..92a86434 100644 --- a/vignettes/splines.Rmd +++ b/vignettes/splines.Rmd @@ -180,7 +180,57 @@ 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`) +## 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. + +```{r test-smooth-endpoints, eval = run_smooth_vignette} +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): + +```{r test-smooth-average, eval = run_smooth_vignette} +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")`. + +## 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 @@ -223,7 +273,7 @@ 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 +## 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 From 7ebe11e3d004544c201a1b707208303663f5222d Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Wed, 5 Aug 2026 12:01:34 +0930 Subject: [PATCH 06/11] Update NEWS to reflect changes in significance colouring defaults for plotting functions. The default now uses pH0 instead of FDR, with guidance on retaining FDR-based colouring. Adjusted messaging for Bayesian FDR to display only when applicable. --- inst/NEWS.rd | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/inst/NEWS.rd b/inst/NEWS.rd index 9f4f15bd..cf79f2a3 100644 --- a/inst/NEWS.rd +++ b/inst/NEWS.rd @@ -57,8 +57,9 @@ \section{News in version 2.1.25}{ \itemize{ - \item Improved plotting significance colouring controls. Added \code{significance_statistic} to \code{sccomp_boxplot()} with default \code{c("pH0", "FDR")}, so colouring now defaults to posterior probability while still supporting FDR-based colouring. - \item Bayesian FDR messaging is now shown only when FDR is selected, avoiding FDR-specific text when probability-based significance is used. + \item \strong{Breaking plotting-default change: significance colouring now uses pH0 instead of FDR by default.} This applies consistently to \code{plot()}, \code{sccomp_plot_intervals_1D()}, \code{sccomp_plot_intervals_2D()}, and \code{sccomp_boxplot()}. In each plot, an effect is highlighted when its posterior null-hypothesis probability satisfies \code{pH0 < significance_threshold}. + \item To retain the previous FDR-based colouring, explicitly set \code{significance_statistic = "FDR"}. The plot will then highlight effects satisfying \code{FDR < significance_threshold}. Set \code{significance_statistic = "pH0"} explicitly when reproducible plotting behaviour across package versions is required. + \item Bayesian FDR explanatory text is shown only when \code{significance_statistic = "FDR"} (and \code{show_fdr_message = TRUE}); it is not shown for the default pH0-based colouring. \item Fixed the package version dot-numbering (\url{https://github.com/MangiolaLaboratory/sccomp/issues/256}). }} From 77ace249ee0df5cf5855914104c95a95c3d317f3 Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Thu, 23 Jul 2026 10:00:59 +0930 Subject: [PATCH 07/11] Add sccomp_test_smooth function for hypothesis testing along smooth terms --- NAMESPACE | 5 + R/sccomp_test_smooth.R | 314 ++++++++++++++++++++++++++++++ man/sccomp_test_smooth.Rd | 124 ++++++++++++ tests/testthat/test-smooth-test.R | 152 +++++++++++++++ vignettes/splines.Rmd | 54 ++++- 5 files changed, 647 insertions(+), 2 deletions(-) create mode 100644 R/sccomp_test_smooth.R create mode 100644 man/sccomp_test_smooth.Rd create mode 100644 tests/testthat/test-smooth-test.R diff --git a/NAMESPACE b/NAMESPACE index c5866835..bd119351 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -13,6 +13,7 @@ S3method(sccomp_remove_outliers,sccomp_tbl) S3method(sccomp_remove_unwanted_effects,sccomp_tbl) S3method(sccomp_replicate,sccomp_tbl) S3method(sccomp_test,sccomp_tbl) +S3method(sccomp_test_smooth,sccomp_tbl) S3method(simulate_data,tbl) export(clear_draw_files) export(clear_stan_model_cache) @@ -32,6 +33,7 @@ export(sccomp_replicate) export(sccomp_scatterplot) export(sccomp_stan_models_cache_dir) export(sccomp_test) +export(sccomp_test_smooth) export(sccomp_theme) export(simulate_data) import(dplyr) @@ -68,6 +70,8 @@ importFrom(dplyr,rowwise) importFrom(dplyr,select) importFrom(dplyr,slice) importFrom(dplyr,summarise) +importFrom(dplyr,tibble) +importFrom(dplyr,ungroup) importFrom(dplyr,where) importFrom(dplyr,with_groups) importFrom(forcats,fct_inorder) @@ -148,6 +152,7 @@ importFrom(rlang,quo_is_symbolic) importFrom(rlang,quo_name) importFrom(rlang,quo_squash) importFrom(rlang,set_names) +importFrom(rlang,sym) importFrom(scales,trans_new) importFrom(stats,C) importFrom(stats,as.formula) diff --git a/R/sccomp_test_smooth.R b/R/sccomp_test_smooth.R new file mode 100644 index 00000000..8948d09c --- /dev/null +++ b/R/sccomp_test_smooth.R @@ -0,0 +1,314 @@ +# Hypothesis testing for smooth (spline) terms. +# +# `sccomp_test_smooth()` tests whether the mean composition differs between two +# positions (or two intervals) of a continuous covariate entered as a smooth +# `s()` term. It is a thin wrapper that: +# 1. predicts the *mean* linear predictor (`mu_unconstrained`, the logit-scale +# mean — no beta-binomial overdispersion) at the requested covariate +# positions via the existing replicate-data / generated-quantities pass, +# 2. forms the per-draw contrast `mu(to) - mu(from)` for each cell group, and +# 3. feeds those contrast draws into the same `draws_to_statistics()` engine +# used by `sccomp_test()`, so the returned `c_pH0` / `c_FDR` mean exactly +# what they mean elsewhere in sccomp. +# +# The contrast lives on the logit (unconstrained) scale, which is the scale +# `test_composition_above_logit_fold_change` is calibrated for. + +#' Test differences along a smooth term +#' +#' @description +#' Tests whether the mean cell-group composition differs between two positions +#' (or two intervals) of a continuous covariate modelled with a smooth `s()` +#' term. The test is a contrast of the model's mean prediction on the logit +#' scale, evaluated at the two positions and summarised with the same +#' probability-of-null (`pH0`) and false-discovery-rate (`FDR`) machinery as +#' [sccomp_test()]. +#' +#' @details +#' Because the smooth curve can be non-linear, "difference across an interval" +#' is ambiguous. Two comparison modes are provided: +#' \itemize{ +#' \item \code{"endpoints"} (default): \code{from} and \code{to} are single +#' values; the contrast is \eqn{\mu(\code{to}) - \mu(\code{from})}. +#' \item \code{"average"}: \code{from} and \code{to} are length-2 intervals +#' \code{c(lo, hi)}; the contrast is the average of the curve over the +#' \code{to} interval minus the average over the \code{from} interval +#' (each interval sampled at \code{resolution} points). +#' } +#' +#' Any covariate not being varied cancels out of the contrast because the +#' logit-scale predictor is additive; it therefore does not matter what value +#' those covariates take, as long as it is held constant. For factor smooths +#' (\code{s(x, g, bs = "fs")}) or by-factor smooths (\code{s(x, by = g)}) the +#' grouping factor selects \emph{which} curve is tested; set it via \code{at}, +#' e.g. \code{at = list(tissue = "tumor")}. +#' +#' The test is about the \strong{mean} composition: it uses the expected +#' linear predictor (\code{mu_unconstrained}), not the overdispersed +#' beta-binomial realisation. It therefore answers "does the expected +#' composition differ between these covariate positions?", not "would an +#' individual sample differ?". +#' +#' @param fit The result of [sccomp_estimate()] with a smooth term in +#' \code{formula_composition}. +#' @param smooth Character. The continuous covariate inside the smooth to test +#' (e.g. \code{"pseudotime"}). Can be \code{NULL} when the model has exactly +#' one smooth term over one continuous covariate, in which case it is +#' auto-detected. +#' @param from Numeric. The reference position. A single value for +#' \code{comparison = "endpoints"}, or a length-2 interval \code{c(lo, hi)} +#' for \code{comparison = "average"}. +#' @param to Numeric. The comparison position, same shape as \code{from}. +#' @param comparison One of \code{"endpoints"} or \code{"average"}. See details. +#' @param at Optional named list fixing the value of other covariates (e.g. +#' the grouping factor of a factor smooth). Covariates not supplied are held +#' at their first observed value (which cancels for non-grouping covariates). +#' @param resolution Integer. Number of grid points used to approximate each +#' interval average when \code{comparison = "average"}. +#' @param test_composition_above_logit_fold_change Positive numeric. Effect +#' threshold for the hypothesis test, on the logit scale — identical meaning +#' to the argument of [sccomp_test()]. +#' @param percent_false_positive Numeric in (0, 100). Used for the credible +#' interval width, as in [sccomp_test()]. +#' @param number_of_draws Integer. Number of posterior draws used for the +#' prediction pass. +#' @param mcmc_seed Integer. Seed for the generated-quantities pass. +#' @param robust Logical. Currently unused placeholder for API parity with +#' [sccomp_predict()]; the effect is always the posterior mean of the +#' contrast. +#' +#' @return A tibble with one row per cell group: +#' \itemize{ +#' \item \code{cell_group} — the cell group tested. +#' \item \code{smooth} — the continuous covariate tested. +#' \item \code{from}, \code{to} — the compared positions (as labels). +#' \item \code{c_lower}, \code{c_effect}, \code{c_upper} — 95% CI and posterior +#' mean of the logit-scale contrast \eqn{\mu(\code{to}) - \mu(\code{from})}. +#' \item \code{c_pH0} — probability the effect is within the null region. +#' \item \code{c_FDR} — false-discovery rate across cell groups. +#' } +#' +#' @examples +#' \donttest{ +#' if (instantiate::stan_cmdstan_exists()) { +#' data("counts_obj") +#' # add a continuous covariate +#' counts_obj$pseudotime <- as.numeric(factor(counts_obj$sample)) +#' +#' fit <- sccomp_estimate( +#' counts_obj, +#' ~ s(pseudotime, k = 4), ~ 1, "sample", "cell_group", "count", +#' cores = 1 +#' ) +#' +#' fit |> sccomp_test_smooth(from = 2, to = 8) +#' } +#' } +#' +#' @export +sccomp_test_smooth <- function(fit, + smooth = NULL, + from, + to, + comparison = c("endpoints", "average"), + at = NULL, + resolution = 20L, + test_composition_above_logit_fold_change = 0.1, + percent_false_positive = 5, + number_of_draws = 500, + mcmc_seed = sample_seed(), + robust = FALSE) { + check_and_install_cmdstanr() + UseMethod("sccomp_test_smooth", fit) +} + +#' @importFrom dplyr tibble bind_rows left_join group_by summarise ungroup arrange mutate select distinct pull +#' @importFrom tidyr pivot_wider nest unnest +#' @importFrom rlang enquo quo_name sym +#' @export +sccomp_test_smooth.sccomp_tbl <- function(fit, + smooth = NULL, + from, + to, + comparison = c("endpoints", "average"), + at = NULL, + resolution = 20L, + test_composition_above_logit_fold_change = 0.1, + percent_false_positive = 5, + number_of_draws = 500, + mcmc_seed = sample_seed(), + robust = FALSE) { + + # Define the variables as NULL to avoid CRAN NOTES + M <- NULL + .chain <- NULL + .iteration <- NULL + .draw <- NULL + .value <- NULL + .smooth_group <- NULL + parameter <- NULL + + comparison <- match.arg(comparison) + .sample <- attr(fit, ".sample") + .cell_group <- attr(fit, ".cell_group") + model_input <- attr(fit, "model_input") + smooth_specs <- get_smooth_results(fit)$smooth_specs + + if (is.null(smooth_specs) || length(smooth_specs) == 0) + stop("sccomp says: this fit has no smooth (s()/t2()) terms; use sccomp_test() for parametric contrasts.") + + # ---- Resolve which continuous covariate to test ------------------------- + # For each smooth, the continuous variable(s) are `term` minus the grouping + # factor `fterm` (present only for factor smooths). + smooth_continuous <- lapply(smooth_specs, function(sm) setdiff(sm$term, sm$fterm)) + all_continuous <- unique(unlist(smooth_continuous)) + + if (is.null(smooth)) { + if (length(all_continuous) != 1) + stop(sprintf( + "sccomp says: the model has several smooth covariates (%s); specify one via `smooth = `.", + paste(all_continuous, collapse = ", ") + )) + smooth <- all_continuous + } else if (!smooth %in% all_continuous) { + stop(sprintf( + "sccomp says: `%s` is not a continuous smooth covariate in this model (available: %s).", + smooth, paste(all_continuous, collapse = ", ") + )) + } + + # ---- Validate from/to shapes against comparison mode -------------------- + if (comparison == "endpoints") { + if (length(from) != 1L || length(to) != 1L) + stop("sccomp says: for comparison = 'endpoints', `from` and `to` must be single values.") + from_grid <- from + to_grid <- to + from_label <- as.character(from) + to_label <- as.character(to) + } else { + if (length(from) != 2L || length(to) != 2L) + stop("sccomp says: for comparison = 'average', `from` and `to` must be length-2 intervals c(lo, hi).") + from_grid <- seq(from[1], from[2], length.out = resolution) + to_grid <- seq(to[1], to[2], length.out = resolution) + from_label <- sprintf("[%s, %s]", from[1], from[2]) + to_label <- sprintf("[%s, %s]", to[1], to[2]) + } + + # ---- Build new_data at the requested positions -------------------------- + new_data <- build_smooth_test_newdata( + fit, smooth, from_grid, to_grid, at, .sample + ) + + # ---- Predict the mean linear predictor (logit scale) -------------------- + rng <- replicate_data( + fit, + formula_composition = NULL, # use the fitted model + formula_variability = ~ 1, + new_data = new_data, + number_of_draws = number_of_draws, + mcmc_seed = mcmc_seed + ) + + sample_col <- quo_name(.sample) + sample_names <- new_data |> pull(!!.sample) + + # Per-draw logit-scale mean (mu_unconstrained), keeping chain/iteration so + # the diagnostics in draws_to_statistics remain meaningful. + draws <- + rng |> + draws_to_tibble_x_y("mu_unconstrained", "M", "N") |> + # Map N -> sample name -> group ("from"/"to") + left_join( + tibble( + N = seq_along(sample_names), + !!sample_col := sample_names + ), + by = "N" + ) |> + left_join( + new_data |> select(!!.sample, .smooth_group), + by = sample_col + ) |> + # Map M -> cell group name + left_join( + tibble( + M = seq_len(ncol(model_input$y)), + !!quo_name(.cell_group) := colnames(model_input$y) + ), + by = "M" + ) + + # ---- Form the contrast per (cell_group, draw) --------------------------- + # Average over grid points within each group (a no-op for endpoints, where + # each group has one point), then difference to - from. + contrast_draws <- + draws |> + group_by(!!.cell_group, M, .chain, .iteration, .draw, .smooth_group) |> + summarise(.value = mean(.value), .groups = "drop") |> + pivot_wider(names_from = .smooth_group, values_from = .value) |> + mutate(.value = to - from) |> + mutate(parameter = sprintf("%s: %s -> %s", smooth, from_label, to_label)) |> + select(!!.cell_group, M, parameter, .chain, .iteration, .draw, .value) + + # ---- Reuse the standard pH0 / FDR engine -------------------------------- + result <- + contrast_draws |> + draws_to_statistics( + percent_false_positive / 100, + test_composition_above_logit_fold_change, + !!.cell_group, + "c_" + ) |> + select(-M, -parameter) |> + mutate(smooth = smooth, from = from_label, to = to_label) |> + select(!!.cell_group, smooth, from, to, everything()) + + result +} + + +#' Build the sample-wise new_data for a smooth test +#' +#' Two groups of rows ("from" and "to"), each covering the requested grid of +#' the continuous covariate. All other covariates in the composition formula +#' are held at a constant value (from `at`, else the first observed value), so +#' they cancel in the contrast; a grouping factor supplied via `at` selects +#' which curve of a factor / by-factor smooth is tested. +#' +#' @keywords internal +#' @noRd +#' @importFrom dplyr tibble bind_rows distinct mutate row_number all_of +build_smooth_test_newdata <- function(fit, smooth, from_grid, to_grid, at, .sample) { + + sample_col <- quo_name(.sample) + count_data <- attr(fit, "count_data") + formula_composition <- attr(fit, "formula_composition") + + covariates <- parse_formula(formula_composition) + other_covs <- setdiff(covariates, smooth) + + # Defaults for non-varied covariates: `at` override, else first observed. + ref_row <- count_data |> distinct(!!.sample, .keep_all = TRUE) |> head(1) + cov_value <- function(v) { + if (!is.null(at) && v %in% names(at)) return(at[[v]]) + if (v %in% colnames(count_data)) return(ref_row[[v]]) + stop(sprintf("sccomp says: covariate `%s` needed by the model is missing from the data and not supplied via `at`.", v)) + } + other_values <- lapply(other_covs, cov_value) + names(other_values) <- other_covs + + make_group <- function(x_grid, group_label) { + df <- tibble(!!smooth := x_grid, .smooth_group = group_label) + for (v in other_covs) df[[v]] <- other_values[[v]] + df + } + + new_data <- + bind_rows( + make_group(from_grid, "from"), + make_group(to_grid, "to") + ) |> + mutate(!!sample_col := sprintf("smoothtest_%03d", row_number())) + + new_data +} diff --git a/man/sccomp_test_smooth.Rd b/man/sccomp_test_smooth.Rd new file mode 100644 index 00000000..9ee96be9 --- /dev/null +++ b/man/sccomp_test_smooth.Rd @@ -0,0 +1,124 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/sccomp_test_smooth.R +\name{sccomp_test_smooth} +\alias{sccomp_test_smooth} +\title{Test differences along a smooth term} +\usage{ +sccomp_test_smooth( + fit, + smooth = NULL, + from, + to, + comparison = c("endpoints", "average"), + at = NULL, + resolution = 20L, + test_composition_above_logit_fold_change = 0.1, + percent_false_positive = 5, + number_of_draws = 500, + mcmc_seed = sample_seed(), + robust = FALSE +) +} +\arguments{ +\item{fit}{The result of \code{\link[=sccomp_estimate]{sccomp_estimate()}} with a smooth term in +\code{formula_composition}.} + +\item{smooth}{Character. The continuous covariate inside the smooth to test +(e.g. \code{"pseudotime"}). Can be \code{NULL} when the model has exactly +one smooth term over one continuous covariate, in which case it is +auto-detected.} + +\item{from}{Numeric. The reference position. A single value for +\code{comparison = "endpoints"}, or a length-2 interval \code{c(lo, hi)} +for \code{comparison = "average"}.} + +\item{to}{Numeric. The comparison position, same shape as \code{from}.} + +\item{comparison}{One of \code{"endpoints"} or \code{"average"}. See details.} + +\item{at}{Optional named list fixing the value of other covariates (e.g. +the grouping factor of a factor smooth). Covariates not supplied are held +at their first observed value (which cancels for non-grouping covariates).} + +\item{resolution}{Integer. Number of grid points used to approximate each +interval average when \code{comparison = "average"}.} + +\item{test_composition_above_logit_fold_change}{Positive numeric. Effect +threshold for the hypothesis test, on the logit scale — identical meaning +to the argument of \code{\link[=sccomp_test]{sccomp_test()}}.} + +\item{percent_false_positive}{Numeric in (0, 100). Used for the credible +interval width, as in \code{\link[=sccomp_test]{sccomp_test()}}.} + +\item{number_of_draws}{Integer. Number of posterior draws used for the +prediction pass.} + +\item{mcmc_seed}{Integer. Seed for the generated-quantities pass.} + +\item{robust}{Logical. Currently unused placeholder for API parity with +\code{\link[=sccomp_predict]{sccomp_predict()}}; the effect is always the posterior mean of the +contrast.} +} +\value{ +A tibble with one row per cell group: +\itemize{ +\item \code{cell_group} — the cell group tested. +\item \code{smooth} — the continuous covariate tested. +\item \code{from}, \code{to} — the compared positions (as labels). +\item \code{c_lower}, \code{c_effect}, \code{c_upper} — 95\% CI and posterior +mean of the logit-scale contrast \eqn{\mu(\code{to}) - \mu(\code{from})}. +\item \code{c_pH0} — probability the effect is within the null region. +\item \code{c_FDR} — false-discovery rate across cell groups. +} +} +\description{ +Tests whether the mean cell-group composition differs between two positions +(or two intervals) of a continuous covariate modelled with a smooth \code{s()} +term. The test is a contrast of the model's mean prediction on the logit +scale, evaluated at the two positions and summarised with the same +probability-of-null (\code{pH0}) and false-discovery-rate (\code{FDR}) machinery as +\code{\link[=sccomp_test]{sccomp_test()}}. +} +\details{ +Because the smooth curve can be non-linear, "difference across an interval" +is ambiguous. Two comparison modes are provided: +\itemize{ +\item \code{"endpoints"} (default): \code{from} and \code{to} are single +values; the contrast is \eqn{\mu(\code{to}) - \mu(\code{from})}. +\item \code{"average"}: \code{from} and \code{to} are length-2 intervals +\code{c(lo, hi)}; the contrast is the average of the curve over the +\code{to} interval minus the average over the \code{from} interval +(each interval sampled at \code{resolution} points). +} + +Any covariate not being varied cancels out of the contrast because the +logit-scale predictor is additive; it therefore does not matter what value +those covariates take, as long as it is held constant. For factor smooths +(\code{s(x, g, bs = "fs")}) or by-factor smooths (\code{s(x, by = g)}) the +grouping factor selects \emph{which} curve is tested; set it via \code{at}, +e.g. \code{at = list(tissue = "tumor")}. + +The test is about the \strong{mean} composition: it uses the expected +linear predictor (\code{mu_unconstrained}), not the overdispersed +beta-binomial realisation. It therefore answers "does the expected +composition differ between these covariate positions?", not "would an +individual sample differ?". +} +\examples{ +\donttest{ + if (instantiate::stan_cmdstan_exists()) { + data("counts_obj") + # add a continuous covariate + counts_obj$pseudotime <- as.numeric(factor(counts_obj$sample)) + + fit <- sccomp_estimate( + counts_obj, + ~ s(pseudotime, k = 4), ~ 1, "sample", "cell_group", "count", + cores = 1 + ) + + fit |> sccomp_test_smooth(from = 2, to = 8) + } +} + +} diff --git a/tests/testthat/test-smooth-test.R b/tests/testthat/test-smooth-test.R new file mode 100644 index 00000000..24e2b661 --- /dev/null +++ b/tests/testthat/test-smooth-test.R @@ -0,0 +1,152 @@ +test_that("sccomp_test_smooth() errors on a fit without smooth terms", { + skip_if_not_installed("mgcv") + skip_cmdstan() + data("counts_obj") + + fit <- sccomp_estimate( + counts_obj, + formula_composition = ~ type, + formula_variability = ~ 1, + sample = "sample", + cell_group = "cell_group", + abundance = "count", + cores = 1, + inference_method = "pathfinder", + max_sampling_iterations = 200, + verbose = FALSE + ) + + expect_error( + sccomp_test_smooth(fit, from = 1, to = 2), + "no smooth" + ) +}) + +test_that("sccomp_test_smooth() endpoints contrast returns the standard columns", { + skip_if_not_installed("mgcv") + skip_cmdstan() + data("counts_obj") + + samples <- levels(counts_obj$sample) + covs <- tibble::tibble( + sample = samples, + pseudotime = seq(0, 6, length.out = length(samples)) + ) + counts_pt <- counts_obj |> dplyr::left_join(covs, by = "sample") + + fit <- sccomp_estimate( + counts_pt, + formula_composition = ~ s(pseudotime, k = 4), + formula_variability = ~ 1, + sample = "sample", + cell_group = "cell_group", + abundance = "count", + cores = 1, + inference_method = "pathfinder", + max_sampling_iterations = 200, + verbose = FALSE + ) + + res <- fit |> + sccomp_test_smooth(from = 1, to = 5, number_of_draws = 100) + + n_groups <- counts_pt |> dplyr::distinct(cell_group) |> nrow() + expect_equal(nrow(res), n_groups) + + expect_true(all( + c("cell_group", "smooth", "from", "to", + "c_lower", "c_effect", "c_upper", "c_pH0", "c_FDR") %in% colnames(res) + )) + + # auto-detected smooth variable, correct labels + expect_equal(unique(res$smooth), "pseudotime") + expect_equal(unique(res$from), "1") + expect_equal(unique(res$to), "5") + + # probabilities and FDR are in [0, 1] + expect_true(all(res$c_pH0 >= 0 & res$c_pH0 <= 1)) + expect_true(all(res$c_FDR >= 0 & res$c_FDR <= 1)) + # CI brackets the point effect + expect_true(all(res$c_lower <= res$c_effect & res$c_effect <= res$c_upper)) +}) + +test_that("sccomp_test_smooth() average mode requires length-2 intervals", { + skip_if_not_installed("mgcv") + skip_cmdstan() + data("counts_obj") + + samples <- levels(counts_obj$sample) + covs <- tibble::tibble( + sample = samples, + pseudotime = seq(0, 6, length.out = length(samples)) + ) + counts_pt <- counts_obj |> dplyr::left_join(covs, by = "sample") + + fit <- sccomp_estimate( + counts_pt, + formula_composition = ~ s(pseudotime, k = 4), + formula_variability = ~ 1, + sample = "sample", + cell_group = "cell_group", + abundance = "count", + cores = 1, + inference_method = "pathfinder", + max_sampling_iterations = 200, + verbose = FALSE + ) + + # wrong shape for average + expect_error( + sccomp_test_smooth(fit, from = 1, to = 5, comparison = "average"), + "length-2" + ) + + res <- fit |> + sccomp_test_smooth( + from = c(0, 2), to = c(4, 6), + comparison = "average", resolution = 6, number_of_draws = 100 + ) + + n_groups <- counts_pt |> dplyr::distinct(cell_group) |> nrow() + expect_equal(nrow(res), n_groups) + expect_equal(unique(res$from), "[0, 2]") + expect_equal(unique(res$to), "[4, 6]") +}) + +test_that("sccomp_test_smooth() selects a curve via `at` for factor smooths", { + skip_if_not_installed("mgcv") + skip_cmdstan() + data("counts_obj") + + samples <- levels(counts_obj$sample) + covs <- tibble::tibble( + sample = samples, + pseudotime = seq(0, 6, length.out = length(samples)), + 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 + ) + + res <- fit |> + sccomp_test_smooth( + smooth = "pseudotime", from = 1, to = 5, + at = list(tissue = "tumor"), number_of_draws = 100 + ) + + n_groups <- counts_pt |> dplyr::distinct(cell_group) |> nrow() + expect_equal(nrow(res), n_groups) + expect_equal(unique(res$smooth), "pseudotime") +}) diff --git a/vignettes/splines.Rmd b/vignettes/splines.Rmd index 80dd6371..92a86434 100644 --- a/vignettes/splines.Rmd +++ b/vignettes/splines.Rmd @@ -180,7 +180,57 @@ 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`) +## 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. + +```{r test-smooth-endpoints, eval = run_smooth_vignette} +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): + +```{r test-smooth-average, eval = run_smooth_vignette} +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")`. + +## 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 @@ -223,7 +273,7 @@ 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 +## 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 From e1b4b0550e8e5bf7d560df12fb6d77611dc10410 Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Fri, 7 Aug 2026 12:56:40 +1000 Subject: [PATCH 08/11] Increase random-effect slot budget from 4 to 5, allowing for coexistence of multiple smooth terms with explicit random effects. Adjusted relevant functions and documentation accordingly. --- DESCRIPTION | 4 +- R/sccomp_estimate.R | 4 +- R/sccomp_remove_outliers.R | 24 +++---- R/sccomp_replicate.R | 27 ++++---- R/sccomp_test.R | 8 +-- R/utilities.R | 18 +++-- inst/NEWS.rd | 5 ++ inst/stan/glm_multi_beta_binomial.stan | 59 ++++++++++++---- ...glm_multi_beta_binomial_generate_data.stan | 47 ++++++++++--- tests/testthat/test-replicate_data.R | 4 +- tests/testthat/test-smooths.R | 68 ++++++++++++++++++- 11 files changed, 201 insertions(+), 67 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5fa885cb..798fc786 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: sccomp Type: Package Title: Differential Composition and Variability Analysis for Single-Cell Data -Version: 2.5.0 -Date: 2026-05-11 +Version: 2.5.1 +Date: 2026-08-07 Authors@R: c(person("Stefano", "Mangiola", email = "stefano.mangiola@unimelb.edu.au", role = c("aut", "cre")), person("Alexandra J.", "Roth-Schulze", role = "aut"), person("Marie", "Trussart", role = "aut"), person("Enrique", "Zozaya-Valdés", role = "aut"), person("Mengyao", "Ma", role = "aut"), person("Zijie", "Gao", role = "aut"), person("Alan F.", "Rubin", role = "aut"), person("Terence P.", "Speed", role = "aut"), person("Heejung", "Shim", role = "aut"), person("Anthony T.", "Papenfuss", role = "aut")) Description: Comprehensive R package for differential composition and variability analysis in single-cell RNA sequencing, CyTOF, and microbiome data. Provides robust Bayesian modeling with outlier detection, random effects, and advanced statistical methods for cell type proportion analysis. Features include probabilistic outlier identification, mixed-effect modeling, differential variability testing, and comprehensive visualization tools. Perfect for cancer research, immunology, developmental biology, and single-cell genomics applications. License: GPL-3 diff --git a/R/sccomp_estimate.R b/R/sccomp_estimate.R index 1a326688..fcb01ec4 100644 --- a/R/sccomp_estimate.R +++ b/R/sccomp_estimate.R @@ -1085,8 +1085,8 @@ sccomp_glm_data_frame_counts = function(.data, "beta", "alpha", "prec_intercept_1", "prec_slope_1", "prec_intercept_2", "prec_slope_2", "prec_sd", - # Random effect outputs - one per slot (1..4) - "random_effect_1", "random_effect_2", "random_effect_3", "random_effect_4", + # Random effect outputs - one per slot (1..5) + "random_effect_1", "random_effect_2", "random_effect_3", "random_effect_4", "random_effect_5", "log_lik" ), sig_figs = sig_figs, diff --git a/R/sccomp_remove_outliers.R b/R/sccomp_remove_outliers.R index 9b101747..4ffc87f2 100644 --- a/R/sccomp_remove_outliers.R +++ b/R/sccomp_remove_outliers.R @@ -212,21 +212,21 @@ sccomp_remove_outliers.sccomp_tbl = function(.estimate, # the slot's columns; unseen matrices are empty). ncol_X_random_eff_new = data_for_model$ncol_X_random_eff, length_X_random_effect_which = data_for_model$ncol_X_random_eff, - ncol_X_random_eff_unseen = rep(0L, 4L), + ncol_X_random_eff_unseen = rep(0L, 5L), create_intercept = FALSE ), # Identity which-indices per slot, generated programmatically setNames( - lapply(seq_len(4L), function(k) + lapply(seq_len(5L), function(k) seq_len(data_for_model$ncol_X_random_eff[k]) |> as.array()), - paste0("X_random_effect_which_", seq_len(4L)) + paste0("X_random_effect_which_", seq_len(5L)) ), # Empty unseen design matrices per slot setNames( - replicate(4L, matrix(0, nrow = nrow(data_for_model$X), ncol = 0), + replicate(5L, matrix(0, nrow = nrow(data_for_model$X), ncol = 0), simplify = FALSE), - paste0("X_random_effect_", seq_len(4L), "_unseen") + paste0("X_random_effect_", seq_len(5L), "_unseen") ) ), @@ -339,7 +339,7 @@ sccomp_remove_outliers.sccomp_tbl = function(.estimate, pars = c( "beta", "alpha", "prec_intercept_1", "prec_slope_1", "prec_intercept_2", "prec_slope_2", "prec_sd", - "random_effect_1", "random_effect_2", "random_effect_3", "random_effect_4" + "random_effect_1", "random_effect_2", "random_effect_3", "random_effect_4", "random_effect_5" ), sig_figs = sig_figs, cache_stan_model = cache_stan_model, @@ -365,19 +365,19 @@ sccomp_remove_outliers.sccomp_tbl = function(.estimate, # Per-slot random-effect pass-throughs (see notes in the first call site) ncol_X_random_eff_new = data_for_model$ncol_X_random_eff, length_X_random_effect_which = data_for_model$ncol_X_random_eff, - ncol_X_random_eff_unseen = rep(0L, 4L), + ncol_X_random_eff_unseen = rep(0L, 5L), create_intercept = FALSE ), setNames( - lapply(seq_len(4L), function(k) + lapply(seq_len(5L), function(k) seq_len(data_for_model$ncol_X_random_eff[k]) |> as.array()), - paste0("X_random_effect_which_", seq_len(4L)) + paste0("X_random_effect_which_", seq_len(5L)) ), setNames( - replicate(4L, matrix(0, nrow = nrow(data_for_model$X), ncol = 0), + replicate(5L, matrix(0, nrow = nrow(data_for_model$X), ncol = 0), simplify = FALSE), - paste0("X_random_effect_", seq_len(4L), "_unseen") + paste0("X_random_effect_", seq_len(5L), "_unseen") ) ), @@ -476,7 +476,7 @@ sccomp_remove_outliers.sccomp_tbl = function(.estimate, pars = c( "beta", "alpha", "prec_intercept_1", "prec_slope_1", "prec_intercept_2", "prec_slope_2", "prec_sd", - "random_effect_1", "random_effect_2", "random_effect_3", "random_effect_4", "log_lik" + "random_effect_1", "random_effect_2", "random_effect_3", "random_effect_4", "random_effect_5", "log_lik" ), cache_stan_model = cache_stan_model, ... diff --git a/R/sccomp_replicate.R b/R/sccomp_replicate.R index 941e9b3f..50327dc7 100644 --- a/R/sccomp_replicate.R +++ b/R/sccomp_replicate.R @@ -126,7 +126,7 @@ sccomp_replicate.sccomp_tbl = function(fit, #' @param Xa Original variability design matrix #' @param N Original number of samples #' @param intercept_in_design Whether intercept is in design -#' @param X_random_effect_slots Length-4 list of original random-effect design +#' @param X_random_effect_slots Length-5 list of original random-effect design #' matrices (one per slot). Empty slots are zero-column matrices. #' @param .sample Quosure for the sample identifier column #' @param .cell_group Quosure for the cell group column @@ -150,7 +150,7 @@ sccomp_replicate.sccomp_tbl = function(fit, #' - model_input: The prepared model input data #' - X_which: Indices for the composition design matrix #' - XA_which: Indices for the variability design matrix -#' - X_random_effect_which_1..4: per-slot indices into the original RE design matrix +#' - X_random_effect_which_1..5: per-slot indices into the original RE design matrix #' - create_intercept: Boolean indicating if intercept should be created #' #' @noRd @@ -334,7 +334,7 @@ prepare_replicate_data = function(X, # ---------------------------------------------------------------------- # Build the per-slot replicate design matrices. # - # For each of the 4 slots: if the slot was active in the original fit + # For each of the 5 slots: if the slot was active in the original fit # (original_grouping_names[k] exists) and the new formula references it, # build a new design matrix restricted to the columns the model saw, and # an index vector mapping new columns back to those of the original matrix. @@ -387,7 +387,7 @@ prepare_replicate_data = function(X, list(X = X_new, X_unseen = X_new_unseen, which = which_idx) } - replicate_slots = map(seq_len(4L), build_replicate_slot) + replicate_slots = map(seq_len(5L), build_replicate_slot) # Append smooth-derived replicate slots (one per smooth term in the # composition formula). They occupy whichever slots come after the @@ -396,9 +396,9 @@ prepare_replicate_data = function(X, n_explicit_re = length(original_grouping_names) n_smooth = length(smooth_replicate_slots) n_used = n_explicit_re + n_smooth - if (n_used > 4L) { + if (n_used > 5L) { stop(sprintf( - "sccomp says: the replicate model needs %d RE slot(s) but only 4 are available.", + "sccomp says: the replicate model needs %d RE slot(s) but only 5 are available.", n_used )) } @@ -409,7 +409,7 @@ prepare_replicate_data = function(X, } # setup default unknown_grouping variable for generated quantities - unknown_grouping = rep(0L, 4L) + unknown_grouping = rep(0L, 5L) list( X = new_X, @@ -422,16 +422,19 @@ prepare_replicate_data = function(X, X_random_effect_2 = replicate_slots[[2]]$X, X_random_effect_3 = replicate_slots[[3]]$X, X_random_effect_4 = replicate_slots[[4]]$X, + X_random_effect_5 = replicate_slots[[5]]$X, X_random_effect_1_unseen = replicate_slots[[1]]$X_unseen, X_random_effect_2_unseen = replicate_slots[[2]]$X_unseen, X_random_effect_3_unseen = replicate_slots[[3]]$X_unseen, X_random_effect_4_unseen = replicate_slots[[4]]$X_unseen, + X_random_effect_5_unseen = replicate_slots[[5]]$X_unseen, X_random_effect_which_1 = replicate_slots[[1]]$which, X_random_effect_which_2 = replicate_slots[[2]]$which, X_random_effect_which_3 = replicate_slots[[3]]$which, X_random_effect_which_4 = replicate_slots[[4]]$which, + X_random_effect_which_5 = replicate_slots[[5]]$which, ncol_X_random_eff_new = map_int(replicate_slots, ~ ncol(.x$X)), ncol_X_random_eff_unseen = map_int(replicate_slots, ~ ncol(.x$X_unseen)), @@ -497,7 +500,7 @@ replicate_data = function(.data, Xa = model_input$Xa, N = model_input$N, intercept_in_design = model_input$intercept_in_design, - X_random_effect_slots = lapply(seq_len(4L), function(k) + X_random_effect_slots = lapply(seq_len(5L), function(k) model_input[[paste0("X_random_effect_", k)]]), .sample = !!.sample, .cell_group = !!.cell_group, @@ -523,8 +526,8 @@ replicate_data = function(.data, model_input$N = prepared_data$N model_input$exposure = prepared_data$exposure - # Per-slot RE design + unseen + which-indices (4 slots) - for (k in seq_len(4L)) { + # Per-slot RE design + unseen + which-indices (5 slots) + for (k in seq_len(5L)) { model_input[[paste0("X_random_effect_", k)]] = prepared_data[[paste0("X_random_effect_", k)]] model_input[[paste0("X_random_effect_", k, "_unseen")]] = prepared_data[[paste0("X_random_effect_", k, "_unseen")]] model_input[[paste0("X_random_effect_which_", k)]] = prepared_data[[paste0("X_random_effect_which_", k)]] @@ -539,9 +542,9 @@ replicate_data = function(.data, model_input$X_which = prepared_data$X_which model_input$XA_which = prepared_data$XA_which - # Length-4 vector of which-index lengths for the random-effect slots + # Length-5 vector of which-index lengths for the random-effect slots model_input$length_X_random_effect_which = - map_int(seq_len(4L), ~ length(prepared_data[[paste0("X_random_effect_which_", .x)]])) + map_int(seq_len(5L), ~ length(prepared_data[[paste0("X_random_effect_which_", .x)]])) # Should I create an intercept for generate quantities? model_input$create_intercept = prepared_data$create_intercept diff --git a/R/sccomp_test.R b/R/sccomp_test.R index f9208b84..808a02de 100644 --- a/R/sccomp_test.R +++ b/R/sccomp_test.R @@ -352,8 +352,8 @@ sccomp_summarise_posterior_for_estimate <- function( prefix = "c_" ) ) - # Random effect blocks: append a summary for each non-empty slot (1..4). - for (k in seq_len(4L)) { + # Random effect blocks: append a summary for each non-empty slot (1..5). + for (k in seq_len(5L)) { if (model_input$ncol_X_random_eff[k] == 0) next X_slot <- model_input[[paste0("X_random_effect_", k)]] abundance_parts <- c( @@ -495,7 +495,7 @@ build_stan_parameter_subset <- function(contrasts, design_columns, stan_paramete } # ---------------------------------------------------------------------- -# Random effect draws: extract one slot at a time (1..4) and left-join into +# Random effect draws: extract one slot at a time (1..5) and left-join into # `draws`. Per-slot logic is identical, so we loop over a helper instead of # duplicating the block once per slot. # ---------------------------------------------------------------------- @@ -611,7 +611,7 @@ get_abundance_contrast_draws = function(.data, contrasts = NULL){ random_effect_covariates_all = character(0) - for (k in seq_len(4L)) { + for (k in seq_len(5L)) { if (model_input$ncol_X_random_eff[k] == 0) next res <- add_random_effect_draws(draws, contrasts, model_input, k, attr(.data, "fit")) draws <- res$draws diff --git a/R/utilities.R b/R/utilities.R index dbb03465..1cf816d4 100755 --- a/R/utilities.R +++ b/R/utilities.R @@ -1102,7 +1102,7 @@ data_spread_to_model_input = # so the Stan-side data list is uniformly shaped regardless of how many # clauses the user wrote. # ---------------------------------------------------------------------- - N_RE_SLOTS = 4L + N_RE_SLOTS = 5L n_rows_design = nrow(.data_spread) empty_re_slot = list( @@ -1213,6 +1213,7 @@ data_spread_to_model_input = X_random_effect_2 = re_slots[[2]]$X X_random_effect_3 = re_slots[[3]]$X X_random_effect_4 = re_slots[[4]]$X + X_random_effect_5 = re_slots[[5]]$X # NOTE: per-slot $X_unseen matrices are available in `re_slots[[k]]$X_unseen` # but are not shipped via data_for_model (downstream replicate / outlier @@ -1222,6 +1223,7 @@ data_spread_to_model_input = group_factor_indexes_for_covariance_2 = re_slots[[2]]$gfi group_factor_indexes_for_covariance_3 = re_slots[[3]]$gfi group_factor_indexes_for_covariance_4 = re_slots[[4]]$gfi + group_factor_indexes_for_covariance_5 = re_slots[[5]]$gfi ncol_X_random_eff = map_int(re_slots, "ncol") n_groups = map_int(re_slots, "n_groups") @@ -1261,20 +1263,22 @@ data_spread_to_model_input = bimodal_mean_variability_association = bimodal_mean_variability_association, use_data = use_data, - # Random intercept - 4 uniform slots (see Stan glm_multi_beta_binomial.stan) + # Random intercept - 5 uniform slots (see Stan glm_multi_beta_binomial.stan) is_random_effect = is_random_effect, n_random_eff = n_random_eff, - ncol_X_random_eff = ncol_X_random_eff, # length 4 - n_groups = n_groups, # length 4 - how_many_factors_in_random_design = how_many_factors_in_random_design, # length 4 + ncol_X_random_eff = ncol_X_random_eff, # length 5 + n_groups = n_groups, # length 5 + how_many_factors_in_random_design = how_many_factors_in_random_design, # length 5 X_random_effect_1 = X_random_effect_1, X_random_effect_2 = X_random_effect_2, X_random_effect_3 = X_random_effect_3, X_random_effect_4 = X_random_effect_4, + X_random_effect_5 = X_random_effect_5, group_factor_indexes_for_covariance_1 = group_factor_indexes_for_covariance_1, group_factor_indexes_for_covariance_2 = group_factor_indexes_for_covariance_2, group_factor_indexes_for_covariance_3 = group_factor_indexes_for_covariance_3, group_factor_indexes_for_covariance_4 = group_factor_indexes_for_covariance_4, + group_factor_indexes_for_covariance_5 = group_factor_indexes_for_covariance_5, # For parallel chains grainsize = 1, @@ -1366,8 +1370,8 @@ data_spread_to_model_input = nrow() } - # Default all grouping known (four RE slots; see glm_multi_beta_binomial_generate_data.stan) - data_for_model$unknown_grouping = rep(0L, 4L) + # Default all grouping known (five RE slots; see glm_multi_beta_binomial_generate_data.stan) + data_for_model$unknown_grouping = rep(0L, 5L) # Smooth-term metadata is R-only (mgcv `smoothCon` / `smooth2random` # objects, used by prediction / replicate helpers). It is NOT shipped to diff --git a/inst/NEWS.rd b/inst/NEWS.rd index cf79f2a3..865e7512 100644 --- a/inst/NEWS.rd +++ b/inst/NEWS.rd @@ -1,6 +1,11 @@ \name{NEWS} \title{News for Package \pkg{sccomp}} +\section{News in version 2.5.1}{ +\itemize{ + \item Increased the random-effect slot budget from 4 to 5 so multiple smooth terms (e.g. a global \code{s()} plus a factor-smooth \code{bs = "fs"}) can coexist with an explicit RE clause. +}} + \section{News in version 2.1.34}{ \itemize{ \item Exported \code{sccomp_scatterplot()} for visualising cell-group proportions against a continuous covariate, complementing \code{sccomp_boxplot()} for discrete factors. The function accepts \code{.data}, \code{factor}, \code{significance_threshold}, and \code{remove_unwanted_effects}, and is used by the \code{plot()} method for numeric covariates. diff --git a/inst/stan/glm_multi_beta_binomial.stan b/inst/stan/glm_multi_beta_binomial.stan index e74600de..c60d9498 100755 --- a/inst/stan/glm_multi_beta_binomial.stan +++ b/inst/stan/glm_multi_beta_binomial.stan @@ -88,16 +88,18 @@ functions{ matrix beta, int M, - // Random effects (up to 4 uniform slots) + // Random effects (up to 5 uniform slots) array[] int ncol_X_random_eff, matrix X_random_effect_1, // Sliced matrix X_random_effect_2, // Sliced matrix X_random_effect_3, // Sliced matrix X_random_effect_4, // Sliced + matrix X_random_effect_5, // Sliced matrix random_effect_1, matrix random_effect_2, matrix random_effect_3, matrix random_effect_4, + matrix random_effect_5, // truncation array[,] int truncation_not_idx_minimal @@ -112,6 +114,7 @@ functions{ if(ncol_X_random_eff[2]>0) mu = mu + (X_random_effect_2[idx_y,] * random_effect_2)'; if(ncol_X_random_eff[3]>0) mu = mu + (X_random_effect_3[idx_y,] * random_effect_3)'; if(ncol_X_random_eff[4]>0) mu = mu + (X_random_effect_4[idx_y,] * random_effect_4)'; + if(ncol_X_random_eff[5]>0) mu = mu + (X_random_effect_5[idx_y,] * random_effect_5)'; for(n in 1:N) mu[,n] = softmax(mu[,n]); @@ -334,32 +337,33 @@ data{ int intercept_in_design; // ---------------------------------------------------------------------- - // Random effect blocks: up to 4 uniform "slots", one block per slot. + // Random effect blocks: up to 5 uniform "slots", one block per slot. // Each slot has its own n_factors (K) so that no padding is needed across // slots. A slot with ncol_X_random_eff[k] == 0 is unused (its arrays have // length 0 and no parameters get sampled). // - // slot 1 -> *_1 slot 2 -> *_2 slot 3 -> *_3 slot 4 -> *_4 + // slot 1 -> *_1 ... slot 5 -> *_5 // - // To go past 4 slots, paste another "slot 4" block in this file and bump - // the array length below from 4 to 5. There is no other architectural - // limit. + // To go past 5 slots, paste another "slot 5" block in this file and bump + // the array length below. There is no other architectural limit. // ---------------------------------------------------------------------- int is_random_effect; - array[4] int ncol_X_random_eff; + array[5] int ncol_X_random_eff; matrix[N, ncol_X_random_eff[1]] X_random_effect_1; matrix[N, ncol_X_random_eff[2]] X_random_effect_2; matrix[N, ncol_X_random_eff[3]] X_random_effect_3; matrix[N, ncol_X_random_eff[4]] X_random_effect_4; + matrix[N, ncol_X_random_eff[5]] X_random_effect_5; // Covariance setup (per slot) - array[4] int n_groups; - array[4] int how_many_factors_in_random_design; + array[5] int n_groups; + array[5] int how_many_factors_in_random_design; array[how_many_factors_in_random_design[1], n_groups[1]] int group_factor_indexes_for_covariance_1; array[how_many_factors_in_random_design[2], n_groups[2]] int group_factor_indexes_for_covariance_2; array[how_many_factors_in_random_design[3], n_groups[3]] int group_factor_indexes_for_covariance_3; array[how_many_factors_in_random_design[4], n_groups[4]] int group_factor_indexes_for_covariance_4; + array[how_many_factors_in_random_design[5], n_groups[5]] int group_factor_indexes_for_covariance_5; // LOO int enable_loo; @@ -374,6 +378,7 @@ transformed data{ int ncol_X_random_eff_safe_2 = max(ncol_X_random_eff[2], 1); int ncol_X_random_eff_safe_3 = max(ncol_X_random_eff[3], 1); int ncol_X_random_eff_safe_4 = max(ncol_X_random_eff[4], 1); + int ncol_X_random_eff_safe_5 = max(ncol_X_random_eff[5], 1); // For parallelisation array[N] int array_N; @@ -397,7 +402,7 @@ parameters{ real mix_p; // ---------------------------------------------------------------------- - // Random effect parameters - 4 uniform slots. + // Random effect parameters - 5 uniform slots. // // For each slot k: // * random_effect_raw_k : sum_to_zero_vector[M] per design column @@ -405,7 +410,7 @@ parameters{ // * sigma_correlation_factor_k: per-category Cholesky of correlation matrix // (n_factors[k] x n_factors[k]; 1x1 = no LKJ work) // - // Hyperprior scalars sigma_mu / sigma_sigma are shared in length-4 arrays. + // Hyperprior scalars sigma_mu / sigma_sigma are shared in length-5 arrays. // ---------------------------------------------------------------------- // Slot 1 @@ -428,9 +433,14 @@ parameters{ array[M * (ncol_X_random_eff[4]>0)] vector[how_many_factors_in_random_design[4]] random_effect_sigma_raw_4; array[M * (ncol_X_random_eff[4]>0)] cholesky_factor_corr[how_many_factors_in_random_design[4] * (ncol_X_random_eff[4]>0)] sigma_correlation_factor_4; + // Slot 5 + array[ncol_X_random_eff[5] * (ncol_X_random_eff[5]>0)] sum_to_zero_vector[M] random_effect_raw_5; + array[M * (ncol_X_random_eff[5]>0)] vector[how_many_factors_in_random_design[5]] random_effect_sigma_raw_5; + array[M * (ncol_X_random_eff[5]>0)] cholesky_factor_corr[how_many_factors_in_random_design[5] * (ncol_X_random_eff[5]>0)] sigma_correlation_factor_5; + // Shared hyperprior scalars (one mu, one sigma per slot) - array[4 * (is_random_effect>0)] real random_effect_sigma_mu; - array[4 * (is_random_effect>0)] real random_effect_sigma_sigma; + array[5 * (is_random_effect>0)] real random_effect_sigma_mu; + array[5 * (is_random_effect>0)] real random_effect_sigma_sigma; // For models with a single group (kept from the original design) array[is_random_effect>0] real zero_random_effect; @@ -472,6 +482,7 @@ transformed parameters{ matrix[ncol_X_random_eff_safe_2 * (is_random_effect>0), M] random_effect_2; matrix[ncol_X_random_eff_safe_3 * (is_random_effect>0), M] random_effect_3; matrix[ncol_X_random_eff_safe_4 * (is_random_effect>0), M] random_effect_4; + matrix[ncol_X_random_eff_safe_5 * (is_random_effect>0), M] random_effect_5; if (ncol_X_random_eff[1] > 0) { array[ncol_X_random_eff[1]] vector[M] raw_vec; @@ -516,6 +527,17 @@ transformed parameters{ random_effect_sigma_raw_4, sigma_correlation_factor_4 ); } + + if (ncol_X_random_eff[5] > 0) { + array[ncol_X_random_eff[5]] vector[M] raw_vec; + for (i in 1:ncol_X_random_eff[5]) raw_vec[i] = to_vector(random_effect_raw_5[i]); + random_effect_5 = build_re_block( + M, n_groups[5], how_many_factors_in_random_design[5], ncol_X_random_eff[5], + group_factor_indexes_for_covariance_5, raw_vec, + random_effect_sigma_mu[5], random_effect_sigma_sigma[5], + random_effect_sigma_raw_5, sigma_correlation_factor_5 + ); + } } model{ // Fit main distribution @@ -541,16 +563,18 @@ model{ beta, M, - // Random effects (4 uniform slots) + // Random effects (5 uniform slots) ncol_X_random_eff, X_random_effect_1, X_random_effect_2, X_random_effect_3, X_random_effect_4, + X_random_effect_5, random_effect_1, random_effect_2, random_effect_3, random_effect_4, + random_effect_5, //truncation truncation_not_idx_minimal @@ -643,6 +667,12 @@ model{ for (m in 1:M) random_effect_sigma_raw_4[m] ~ std_normal(); for (m in 1:M) sigma_correlation_factor_4[m] ~ lkj_corr_cholesky(2); } + + if (ncol_X_random_eff[5] > 0) { + for (m in 1:M) random_effect_raw_5[,m] ~ normal(0, inv(sqrt(1 - inv(M)))); + for (m in 1:M) random_effect_sigma_raw_5[m] ~ std_normal(); + for (m in 1:M) sigma_correlation_factor_5[m] ~ lkj_corr_cholesky(2); + } } generated quantities { // LOO @@ -663,6 +693,7 @@ generated quantities { if(ncol_X_random_eff[2]>0) mu = mu + (X_random_effect_2 * random_effect_2)'; if(ncol_X_random_eff[3]>0) mu = mu + (X_random_effect_3 * random_effect_3)'; if(ncol_X_random_eff[4]>0) mu = mu + (X_random_effect_4 * random_effect_4)'; + if(ncol_X_random_eff[5]>0) mu = mu + (X_random_effect_5 * random_effect_5)'; // Calculate proportions diff --git a/inst/stan/glm_multi_beta_binomial_generate_data.stan b/inst/stan/glm_multi_beta_binomial_generate_data.stan index 4427d86c..501b78cc 100755 --- a/inst/stan/glm_multi_beta_binomial_generate_data.stan +++ b/inst/stan/glm_multi_beta_binomial_generate_data.stan @@ -45,40 +45,44 @@ data { int is_random_effect; - // Four uniform random-effect slots (same layout as glm_multi_beta_binomial.stan) - array[4] int ncol_X_random_eff; - array[4] int ncol_X_random_eff_new; + // Five uniform random-effect slots (same layout as glm_multi_beta_binomial.stan) + array[5] int ncol_X_random_eff; + array[5] int ncol_X_random_eff_new; matrix[N, ncol_X_random_eff_new[1]] X_random_effect_1; matrix[N, ncol_X_random_eff_new[2]] X_random_effect_2; matrix[N, ncol_X_random_eff_new[3]] X_random_effect_3; matrix[N, ncol_X_random_eff_new[4]] X_random_effect_4; + matrix[N, ncol_X_random_eff_new[5]] X_random_effect_5; - array[4] int n_groups; - array[4] int how_many_factors_in_random_design; + array[5] int n_groups; + array[5] int how_many_factors_in_random_design; array[how_many_factors_in_random_design[1], n_groups[1]] int group_factor_indexes_for_covariance_1; array[how_many_factors_in_random_design[2], n_groups[2]] int group_factor_indexes_for_covariance_2; array[how_many_factors_in_random_design[3], n_groups[3]] int group_factor_indexes_for_covariance_3; array[how_many_factors_in_random_design[4], n_groups[4]] int group_factor_indexes_for_covariance_4; + array[how_many_factors_in_random_design[5], n_groups[5]] int group_factor_indexes_for_covariance_5; // Per-slot column counts for X_random_effect_which_* arrays (from R: ncol per slot) - array[4] int length_X_random_effect_which; + array[5] int length_X_random_effect_which; array[length_X_random_effect_which[1]] int X_random_effect_which_1; array[length_X_random_effect_which[2]] int X_random_effect_which_2; array[length_X_random_effect_which[3]] int X_random_effect_which_3; array[length_X_random_effect_which[4]] int X_random_effect_which_4; + array[length_X_random_effect_which[5]] int X_random_effect_which_5; int create_intercept; int A_intercept_columns; - array[4] int unknown_grouping; + array[5] int unknown_grouping; - array[4] int ncol_X_random_eff_unseen; + array[5] int ncol_X_random_eff_unseen; matrix[N, ncol_X_random_eff_unseen[1]] X_random_effect_1_unseen; matrix[N, ncol_X_random_eff_unseen[2]] X_random_effect_2_unseen; matrix[N, ncol_X_random_eff_unseen[3]] X_random_effect_3_unseen; matrix[N, ncol_X_random_eff_unseen[4]] X_random_effect_4_unseen; + matrix[N, ncol_X_random_eff_unseen[5]] X_random_effect_5_unseen; } transformed data { matrix[N,1] X_intercept; @@ -88,6 +92,7 @@ transformed data { int ncol_X_random_eff_safe_2 = max(ncol_X_random_eff[2], 1); int ncol_X_random_eff_safe_3 = max(ncol_X_random_eff[3], 1); int ncol_X_random_eff_safe_4 = max(ncol_X_random_eff[4], 1); + int ncol_X_random_eff_safe_5 = max(ncol_X_random_eff[5], 1); } parameters { // Unconstrained vectors so posterior CSVs (rounded sig_figs) still validate; @@ -116,8 +121,12 @@ parameters { array[M * (ncol_X_random_eff[4]>0)] vector[how_many_factors_in_random_design[4]] random_effect_sigma_raw_4; array[M * (ncol_X_random_eff[4]>0)] cholesky_factor_corr[how_many_factors_in_random_design[4] * (ncol_X_random_eff[4]>0)] sigma_correlation_factor_4; - array[4 * (is_random_effect>0)] real random_effect_sigma_mu; - array[4 * (is_random_effect>0)] real random_effect_sigma_sigma; + array[ncol_X_random_eff[5] * (ncol_X_random_eff[5]>0)] vector[M] random_effect_raw_5; + array[M * (ncol_X_random_eff[5]>0)] vector[how_many_factors_in_random_design[5]] random_effect_sigma_raw_5; + array[M * (ncol_X_random_eff[5]>0)] cholesky_factor_corr[how_many_factors_in_random_design[5] * (ncol_X_random_eff[5]>0)] sigma_correlation_factor_5; + + array[5 * (is_random_effect>0)] real random_effect_sigma_mu; + array[5 * (is_random_effect>0)] real random_effect_sigma_sigma; array[is_random_effect>0] real zero_random_effect; } transformed parameters { @@ -142,6 +151,7 @@ transformed parameters { matrix[ncol_X_random_eff_safe_2 * (is_random_effect>0), M] random_effect_2; matrix[ncol_X_random_eff_safe_3 * (is_random_effect>0), M] random_effect_3; matrix[ncol_X_random_eff_safe_4 * (is_random_effect>0), M] random_effect_4; + matrix[ncol_X_random_eff_safe_5 * (is_random_effect>0), M] random_effect_5; if (ncol_X_random_eff[1] > 0) { array[ncol_X_random_eff[1]] vector[M] raw_vec; @@ -183,6 +193,16 @@ transformed parameters { random_effect_sigma_raw_4, sigma_correlation_factor_4 ); } + if (ncol_X_random_eff[5] > 0) { + array[ncol_X_random_eff[5]] vector[M] raw_vec; + for (i in 1:ncol_X_random_eff[5]) raw_vec[i] = normalize_sum_to_zero(random_effect_raw_5[i]); + random_effect_5 = build_re_block( + M, n_groups[5], how_many_factors_in_random_design[5], ncol_X_random_eff[5], + group_factor_indexes_for_covariance_5, raw_vec, + random_effect_sigma_mu[5], random_effect_sigma_sigma[5], + random_effect_sigma_raw_5, sigma_correlation_factor_5 + ); + } } model { } @@ -253,6 +273,13 @@ generated quantities { mu = add_unseen_random_effect_contribution_rng(mu, X_random_effect_4_unseen, ncol_X_random_eff_unseen[4], M); } + if (ncol_X_random_eff[5] > 0) { + mu = add_seen_random_effect_contribution(mu, X_random_effect_5, random_effect_5, + X_random_effect_which_5); + if (ncol_X_random_eff_unseen[5] > 0) + mu = add_unseen_random_effect_contribution_rng(mu, X_random_effect_5_unseen, + ncol_X_random_eff_unseen[5], M); + } matrix[M, N] mu_unconstrained = mu; diff --git a/tests/testthat/test-replicate_data.R b/tests/testthat/test-replicate_data.R index 34e12e52..dff55220 100644 --- a/tests/testthat/test-replicate_data.R +++ b/tests/testthat/test-replicate_data.R @@ -3,9 +3,9 @@ library(dplyr) library(tidyr) library(sccomp) -# Four-slot RE design (see prepare_replicate_data(X_random_effect_slots = )) +# Five-slot RE design (see prepare_replicate_data(X_random_effect_slots = )) re_slots_from_mi <- function(mi) { - lapply(seq_len(4L), function(k) mi[[paste0("X_random_effect_", k)]]) + lapply(seq_len(5L), function(k) mi[[paste0("X_random_effect_", k)]]) } test_that("replicate_data works correctly", { diff --git a/tests/testthat/test-smooths.R b/tests/testthat/test-smooths.R index f08c31ee..bf437c53 100644 --- a/tests/testthat/test-smooths.R +++ b/tests/testthat/test-smooths.R @@ -372,13 +372,76 @@ test_that("fs factor smooth fits and predicts end-to-end (multi-block)", { }) +test_that("global smooth + fs smooth + RE fills all 5 slots end-to-end", { + # Motivating case for the 5-slot budget: one global age smooth, one + # factor-smooth (3 penalty blocks), and one explicit RE clause. + skip_if_not_installed("mgcv") + skip_cmdstan() + data("counts_obj") + + set.seed(42) + samples <- levels(counts_obj$sample) + covs <- tibble::tibble( + sample = samples, + age_days_scaled = as.numeric(scale(seq_along(samples))), + tissue_groups = factor(rep(c("blood", "lymph", "tumor"), + length.out = length(samples))), + dataset_id = factor(rep(c("ds1", "ds2"), + length.out = length(samples))) + ) + counts_age <- counts_obj |> dplyr::left_join(covs, by = "sample") + + fit <- sccomp_estimate( + counts_age, + formula_composition = + ~ s(age_days_scaled, k = 5) + + s(age_days_scaled, tissue_groups, bs = "fs", k = 5) + + (1 | dataset_id), + formula_variability = ~ 1, + sample = "sample", + cell_group = "cell_group", + abundance = "count", + cores = 1, + inference_method = "pathfinder", + max_sampling_iterations = 200, + mcmc_seed = 42, + verbose = FALSE + ) + + mi <- attr(fit, "model_input") + sr <- sccomp:::get_smooth_results(fit) + + # 1 global smooth + 3 fs blocks + 1 RE = 5 occupied slots + expect_equal(length(sr$smooth_labels), 2L) + expect_equal(sum(mi$ncol_X_random_eff > 0L), 5L) + expect_equal(mi$n_random_eff, 5L) + expect_true(all(mi$ncol_X_random_eff > 0L)) + + # Predict across age x tissue; dataset_id held at a seen level. + grid <- expand.grid( + age_days_scaled = seq(min(covs$age_days_scaled), + max(covs$age_days_scaled), length.out = 15), + tissue_groups = levels(covs$tissue_groups), + dataset_id = covs$dataset_id[1] + ) |> + tibble::as_tibble() |> + dplyr::mutate(sample = sprintf("grid_%03d", dplyr::row_number())) + + pred <- sccomp_predict(fit, new_data = grid, number_of_draws = 50) + expect_true(all(c("age_days_scaled", "tissue_groups", "proportion_mean") + %in% names(pred))) + expect_equal(nrow(pred), nrow(grid) * length(unique(counts_obj$cell_group))) + 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 + # 5 smooths + 1 RE clause needs 6 slots; budget is 5. Use the same + # continuous_covariate five times under different `k` (parser doesn't # de-duplicate intentionally — each call is a fresh basis). expect_error( sccomp_estimate( @@ -388,6 +451,7 @@ test_that("smooth slot count beyond N_RE_SLOTS errors cleanly", { s(continuous_covariate, k = 4) + s(continuous_covariate, k = 5) + s(continuous_covariate, k = 6) + + s(continuous_covariate, k = 7) + (1 | group__), formula_variability = ~ 1, sample = "sample", cell_group = "cell_group", From 0e1a4504c9321015ef5bbc379b561df3cd2324a8 Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Fri, 7 Aug 2026 13:00:56 +1000 Subject: [PATCH 09/11] Update random-effect initialization to accommodate 5 slots, enhancing support for multiple smooth terms with explicit random effects. Adjusted related code accordingly. --- R/model_fitting.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/model_fitting.R b/R/model_fitting.R index e463af97..debe4f68 100644 --- a/R/model_fitting.R +++ b/R/model_fitting.R @@ -74,7 +74,7 @@ fit_model = function( if (data_for_model$n_random_eff > 0) { init_list$zero_random_effect = rep(0, size = 1) |> as.array() - for (k in seq_len(4L)) { + for (k in seq_len(5L)) { if (data_for_model$ncol_X_random_eff[k] == 0) next K = data_for_model$how_many_factors_in_random_design[k] From 2be536ddfdc39f3f89b15d36f96b40c0da0c19b9 Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Fri, 7 Aug 2026 13:04:20 +1000 Subject: [PATCH 10/11] update comment --- R/model_fitting.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/model_fitting.R b/R/model_fitting.R index debe4f68..12a88da7 100644 --- a/R/model_fitting.R +++ b/R/model_fitting.R @@ -69,7 +69,7 @@ fit_model = function( init_list$prec_slope_2 = rep(0, data_for_model$A) } - # Random effect inits - 4 uniform slots (one per non-empty random-effect block). + # Random effect inits - 5 uniform slots (one per non-empty random-effect block). # Each slot gets zero-initialised raws + an identity-like correlation Cholesky. if (data_for_model$n_random_eff > 0) { init_list$zero_random_effect = rep(0, size = 1) |> as.array() From 0a29093788b5901ab46a089d92ecb8c240343c6f Mon Sep 17 00:00:00 2001 From: Stefano Mangiola Date: Tue, 11 Aug 2026 11:35:17 +1000 Subject: [PATCH 11/11] Enhance continuous covariate scaling in design matrix functions. Introduced a new `scale_numeric_covariates` function to ensure z-scoring is based on fitted samples, preventing prediction shifts due to varying grid densities. Updated relevant functions and added tests to verify correct behavior. --- NAMESPACE | 1 + R/sccomp_replicate.R | 13 ++- R/utilities.R | 90 ++++++++++++++++++- inst/NEWS.rd | 1 + .../test-continuous-covariate-scaling.R | 70 +++++++++++++++ 5 files changed, 168 insertions(+), 7 deletions(-) create mode 100644 tests/testthat/test-continuous-covariate-scaling.R diff --git a/NAMESPACE b/NAMESPACE index bd119351..a08ad3c0 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -158,6 +158,7 @@ importFrom(stats,C) importFrom(stats,as.formula) importFrom(stats,model.matrix) importFrom(stats,quantile) +importFrom(stats,sd) importFrom(stats,terms) importFrom(stringr,str_detect) importFrom(stringr,str_remove) diff --git a/R/sccomp_replicate.R b/R/sccomp_replicate.R index 50327dc7..9d75fab4 100644 --- a/R/sccomp_replicate.R +++ b/R/sccomp_replicate.R @@ -251,7 +251,12 @@ prepare_replicate_data = function(X, paste(collapse="") |> as.formula(), !!.sample, - accept_NA_as_average_effect = TRUE + accept_NA_as_average_effect = TRUE, + # Continuous covariates are z-scored here. The centre and the scale have + # to come from the fitted samples alone: taking them from `new_data`, + # which carries the old rows too, would make the prediction at a given + # covariate value depend on the range and density of the requested grid. + scaling_reference = old_data ) |> tail(nrow_new_data) %>% # Remove columns that are not in the original design matrix @@ -294,7 +299,8 @@ prepare_replicate_data = function(X, paste(collapse="") |> as.formula(), !!.sample, - accept_NA_as_average_effect = TRUE + accept_NA_as_average_effect = TRUE, + scaling_reference = old_data ) |> tail(nrow_new_data) %>% # Remove columns that are not in the original design matrix @@ -328,7 +334,8 @@ prepare_replicate_data = function(X, mutate(design = map2( formula, grouping, ~ get_random_effect_design3(new_data, .x, .y, !!.sample, - accept_NA_as_average_effect = TRUE) + accept_NA_as_average_effect = TRUE, + scaling_reference = old_data) )) # ---------------------------------------------------------------------- diff --git a/R/utilities.R b/R/utilities.R index 1cf816d4..4ad10571 100755 --- a/R/utilities.R +++ b/R/utilities.R @@ -670,13 +670,18 @@ summary_to_tibble = function(fit, par, x, y = NULL, probs = c(0.025, 0.25, 0.50, #' @noRd get_random_effect_design3 = function( .data_, formula, grouping, .sample, - accept_NA_as_average_effect = FALSE + accept_NA_as_average_effect = FALSE, + scaling_reference = NULL ){ # Define the variables as NULL to avoid CRAN NOTES .sample = enquo(.sample) - mydesign = .data_ |> get_design_matrix(formula, !!.sample, accept_NA_as_average_effect = accept_NA_as_average_effect) + mydesign = .data_ |> get_design_matrix( + formula, !!.sample, + accept_NA_as_average_effect = accept_NA_as_average_effect, + scaling_reference = scaling_reference + ) # Create a matrix of group assignments group_matrix = .data_ |> @@ -731,6 +736,83 @@ get_random_effect_design3 = function( result_long } +#' Z-score the continuous covariates of a design +#' +#' Continuous covariates are z-scored before the design matrix is built. When +#' the design matrix is built for new data, the centre and the scale have to be +#' those of the fitted samples, otherwise the columns end up on a different +#' scale from the one the coefficients were estimated against. `reference` is +#' the data frame those two statistics are taken from, and defaults to the data +#' being scaled, which is the right choice at fit time. +#' +#' @details +#' Only the numeric columns named in `variables` are touched; factors, the +#' sample column and anything absent from `reference` are returned untouched. A +#' covariate that is constant in `reference` becomes zeros rather than `NaN`. +#' +#' ``` +#' training = tibble(sample = c("S1", "S2", "S3"), age = c(30, 40, 50)) +#' +#' scale_numeric_covariates(training, "age") +#' # # A tibble: 3 x 2 +#' # sample age +#' # +#' # 1 S1 -1 +#' # 2 S2 0 +#' # 3 S3 1 +#' +#' # The same age is given the same column value whatever else is in the grid +#' scale_numeric_covariates( +#' tibble(sample = "g1", age = c(40)), "age", reference = training +#' ) +#' # # A tibble: 1 x 2 +#' # sample age +#' # +#' # 1 g1 0 +#' ``` +#' +#' @param .data_spread A data frame with one row per sample. +#' @param variables Character vector of covariate names to consider. +#' @param reference Data frame the centre and scale are computed from. When +#' `NULL` they are computed from `.data_spread` itself. +#' +#' @return `.data_spread`, with its numeric covariates z-scored. +#' +#' @importFrom dplyr select +#' @importFrom dplyr where +#' @importFrom dplyr any_of +#' @importFrom stats sd +#' @noRd +scale_numeric_covariates = function(.data_spread, variables, reference = NULL){ + + if(is.null(reference)) reference = .data_spread + + columns = + .data_spread |> + select(any_of(variables)) |> + select(where(is.numeric)) |> + colnames() |> + intersect(colnames(reference)) + + for(column in columns){ + centre = mean(reference[[column]], na.rm = TRUE) + spread = sd(reference[[column]], na.rm = TRUE) + + # A covariate constant across samples carries no information; keep it + # finite rather than dividing by zero + if(is.na(spread) || spread == 0) spread = 1 + + .data_spread[[column]] = (.data_spread[[column]] - centre) / spread + } + + .data_spread +} + +#' @param scaling_reference Data frame whose continuous covariates give the +#' centre and the scale, passed to [scale_numeric_covariates()]. When `NULL` +#' they come from `.data_spread` itself; supply the fitted samples when +#' building a design matrix for new data. +#' #' @importFrom glue glue #' @importFrom dplyr select #' @importFrom dplyr mutate @@ -739,14 +821,14 @@ get_random_effect_design3 = function( #' @importFrom dplyr where #' @importFrom rlang enquo #' @noRd -get_design_matrix = function(.data_spread, formula, .sample, accept_NA_as_average_effect = FALSE){ +get_design_matrix = function(.data_spread, formula, .sample, accept_NA_as_average_effect = FALSE, scaling_reference = NULL){ .sample = enquo(.sample) .data_spread = .data_spread %>% select(!!.sample, parse_formula(formula)) |> - mutate(across(where(is.numeric), scale)) + scale_numeric_covariates(parse_formula(formula), reference = scaling_reference) # Check for NAs in the data has_na = any(is.na(.data_spread |> select(parse_formula(formula)))) diff --git a/inst/NEWS.rd b/inst/NEWS.rd index 865e7512..5f1f0811 100644 --- a/inst/NEWS.rd +++ b/inst/NEWS.rd @@ -4,6 +4,7 @@ \section{News in version 2.5.1}{ \itemize{ \item Increased the random-effect slot budget from 4 to 5 so multiple smooth terms (e.g. a global \code{s()} plus a factor-smooth \code{bs = "fs"}) can coexist with an explicit RE clause. + \item Fixed the scaling of continuous covariates in \code{sccomp_predict()} and \code{sccomp_replicate()}. Continuous covariates were being z-scored using the fitted samples and the new data pooled together, so the prediction at a given covariate value shifted with the range and the density of the requested grid. The centre and scale of the fitted samples are now used, making predictions a function of the model and the covariate value alone. Smooth terms were never affected, as they are evaluated against the fit-time basis. }} \section{News in version 2.1.34}{ diff --git a/tests/testthat/test-continuous-covariate-scaling.R b/tests/testthat/test-continuous-covariate-scaling.R new file mode 100644 index 00000000..144367f5 --- /dev/null +++ b/tests/testthat/test-continuous-covariate-scaling.R @@ -0,0 +1,70 @@ +library(testthat) +library(sccomp) + +test_that("covariates are z-scored against the reference samples", { + training <- tibble::tibble( + sample = paste0("S", 1:5), + age = c(10, 20, 30, 40, 50), + condition = c("a", "b", "a", "b", "a") + ) + grid <- tibble::tibble(sample = "g1", age = 50, condition = "a") + + scaled <- grid |> + sccomp:::scale_numeric_covariates( + c("age", "condition"), + reference = training + ) + + expect_equal(scaled$age, (50 - 30) / sd(c(10, 20, 30, 40, 50))) + expect_equal(scaled$condition, "a") +}) + +test_that("a covariate constant across samples is not divided by zero", { + training <- tibble::tibble(sample = paste0("S", 1:3), age = c(40, 40, 40)) + + scaled <- training |> sccomp:::scale_numeric_covariates("age") + + expect_equal(scaled$age, c(0, 0, 0)) +}) + +test_that("design matrix columns do not depend on the prediction grid", { + training <- tibble::tibble( + sample = paste0("S", 1:5), + age = c(10, 20, 30, 40, 50) + ) + + # Two grids that share the age of interest but differ in range and density + narrow_grid <- tibble::tibble(sample = c("N1", "N2"), age = c(35, 45)) + wide_grid <- tibble::tibble( + sample = paste0("W", 1:4), + age = c(0, 35, 45, 200) + ) + + column_at_35 <- function(grid) { + design <- training |> + dplyr::bind_rows(grid) |> + sccomp:::get_design_matrix(~ age, sample, scaling_reference = training) + design[grid$sample[grid$age == 35], "age"] + } + + expect_equal(column_at_35(narrow_grid), column_at_35(wide_grid)) + expect_equal( + unname(column_at_35(narrow_grid)), + (35 - 30) / sd(c(10, 20, 30, 40, 50)) + ) +}) + +test_that("design matrix scaling is unchanged when no reference is supplied", { + training <- tibble::tibble( + sample = paste0("S", 1:5), + age = c(10, 20, 30, 40, 50) + ) + + design <- training |> + sccomp:::get_design_matrix(~ age, sample) + + expect_equal( + unname(design[, "age"]), + as.vector(scale(c(10, 20, 30, 40, 50))) + ) +})