sccomp_remove_outliers() is not reproducible despite a fixed mcmc_seed

[shihsama]

Hi!

I am observing non-reproducible results from sccomp_remove_outliers() in sccomp 2.1.30, even when:

• the input object is exactly the same;
• mcmc_seed is explicitly fixed;
• cores = 1;
• the same fitted sccomp_estimate() object is passed to both calls.

We first tested whether the non-reproducibility originated from sccomp_estimate() with the Pathfinder inference method. The results from two independent sccomp_estimate() runs were numerically identical when the seed and number of cores were fixed.

We then isolated sccomp_remove_outliers() by calling it twice on the exact same fitted object. The two calls identified different sample × cell_group observations as outliers and consequently produced different fitted parameter estimates.

This suggests that the issue originates inside sccomp_remove_outliers(), rather than in the initial sccomp_estimate() fit.

Step 1: sccomp_estimate() is reproducible

fit1 <- sce |>
  sccomp_estimate(
    formula_composition = ~ 0 + condition,
    sample = "sample_name",
    cell_group = "L1_celltype",
    inference_method = "pathfinder",
    cores = 1,
    mcmc_seed = 12345,
    verbose = FALSE
  ) |>
  sccomp_test(
    contrasts = c(
      "conditionB - conditionA",
      "conditionC - conditionA",
      "conditionC - conditionB"
    )
  )

fit2 <- sce |>
  sccomp_estimate(
    formula_composition = ~ 0 + condition,
    sample = "sample_name",
    cell_group = "L1_celltype",
    inference_method = "pathfinder",
    cores = 1,
    mcmc_seed = 12345,
    verbose = FALSE
  ) |>
  sccomp_test(
    contrasts = c(
      "conditionB - conditionA",
      "conditionC - conditionA",
      "conditionC - conditionB"
    )
  )

The statistical result columns were identical:

result_columns <- c( "L1_celltype", "parameter", "c_lower", "c_effect", "c_upper", "c_pH0", "c_FDR" )
vapply( result_columns, \(x) identical(fit1[[x]], fit2[[x]]), logical(1) )

All columns returned TRUE.

The only differences reported by comparing the full objects were CmdStanR metadata such as timestamps, temporary output-file paths, and elapsed time.

Therefore, in this setting, sccomp_estimate() using Pathfinder is reproducible when mcmc_seed and cores are fixed.

Step 2: isolate sccomp_remove_outliers()

A single estimate object was fitted once:

fit0 <- sce |>
  sccomp_estimate(
    formula_composition = ~ 0 + condition,
    sample = "sample_name",
    cell_group = "L1_celltype",
    inference_method = "pathfinder",
    cores = 1,
    mcmc_seed = 12345,
    verbose = FALSE
  )

The exact same object was then passed twice to sccomp_remove_outliers() with an explicitly fixed seed:

out1 <- fit0 |>
  sccomp_remove_outliers(
    cores = 1,
    mcmc_seed = 67890,
    verbose = FALSE
  )

out2 <- fit0 |>
  sccomp_remove_outliers(
    cores = 1,
    mcmc_seed = 67890,
    verbose = FALSE
  )

However, the truncation indices were different:

identical(
  attr(out1, "model_input")$truncation_not_idx,
  attr(out2, "model_input")$truncation_not_idx
)

Output:

[1] FALSE

The outlier classifications also differed:

dplyr::full_join(
  attr(out1, "outliers") |>
    dplyr::select(
      sample_name,
      L1_celltype,
      outlier_1 = outlier
    ),
  attr(out2, "outliers") |>
    dplyr::select(
      sample_name,
      L1_celltype,
      outlier_2 = outlier
    ),
  by = c("sample_name", "L1_celltype")
) |>
  dplyr::filter(outlier_1 != outlier_2)

Output:

# A tibble: 2 × 4
  sample_name L1_celltype outlier_1 outlier_2
  <chr>       <chr>       <lgl>     <lgl>
1 Sample11        NK          FALSE     TRUE
2 Sample3      Plasma      TRUE      FALSE

Because different observations were excluded, the fitted composition parameters from out1 and out2 were also slightly different.

Expected behaviour

Given:

• the exact same input sccomp_estimate() object;
• the same mcmc_seed;
• cores = 1;
• the same software and hardware environment;

I would expect repeated calls to sccomp_remove_outliers() to produce the same:

• posterior predictive quantities;
• outlier classifications;
• truncation_not_idx;
• final fitted parameters.
  Actual behaviour

Repeated calls to sccomp_remove_outliers() produce different outlier classifications and different truncation indices, despite an explicitly fixed mcmc_seed.

Possible source of the issue

From inspecting the implementation, mcmc_seed appears to be passed to the model-fitting steps inside sccomp_remove_outliers().

However, the posterior predictive generated-quantities steps used to identify outliers appear to involve stochastic RNG operations, and the public mcmc_seed does not appear to be explicitly forwarded to all of those generated-quantities calls.

If those calls use an independently generated seed, the posterior predictive counts and their quantiles can differ between runs. Observations near the outlier threshold may then switch between TRUE and FALSE, which changes truncation_not_idx and the final model fit.

Could you please confirm whether the mcmc_seed argument is propagated to every stochastic step inside sccomp_remove_outliers(), including the generated-quantities calls used for posterior predictive outlier detection?

Relevance to issue #135

I noticed issue #135 and the related change allowing custom posterior output samples for variational Bayes.

However, this example uses:

inference_method = "pathfinder"

with sccomp 2.1.30, and the initial Pathfinder fit is reproducible. The non-reproducibility appears specifically after calling sccomp_remove_outliers().

Environment

packageVersion("sccomp")
[1] ‘2.1.30’
packageVersion("cmdstanr")
[1] ‘0.9.0’
cmdstanr::cmdstan_version()
[1] "2.39.0"

sessionInfo()
R version 4.5.2 (2025-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)

Matrix products: default
LAPACK version 3.12.1

locale:
[1] LC_COLLATE=Chinese (Simplified)_China.utf8 LC_CTYPE=Chinese (Simplified)_China.utf8 LC_MONETARY=Chinese (Simplified)_China.utf8
[4] LC_NUMERIC=C LC_TIME=Chinese (Simplified)_China.utf8

time zone: Asia/Shanghai
tzcode source: internal

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] sccomp_2.1.30 instantiate_0.2.3 anndataR_1.0.2 reticulate_1.44.1 RColorBrewer_1.1-3 Nebulosa_1.0.1
[7] data.table_1.18.2.1 patchwork_1.3.2 qs_0.27.3 lubridate_1.9.5 forcats_1.0.1 stringr_1.6.0
[13] dplyr_1.2.0 purrr_1.2.1 readr_2.1.6 tidyr_1.3.2 tibble_3.3.1 ggplot2_4.0.2
[19] tidyverse_2.0.0 Seurat_5.5.0 SeuratObject_5.4.0 sp_2.2-0

loaded via a namespace (and not attached):
[1] RcppAnnoy_0.0.23 splines_4.5.2 later_1.4.8 polyclip_1.10-7
[5] fastDummies_1.7.5 lifecycle_1.0.5 vroom_1.7.0 globals_0.19.1
[9] processx_3.8.6 lattice_0.22-7 MASS_7.3-65 backports_1.5.0
[13] magrittr_2.0.4 plotly_4.12.0 httpuv_1.6.16 otel_0.2.0
[17] sctransform_0.4.3 spam_2.11-3 spatstat.sparse_3.1-0 cowplot_1.2.0
[21] pbapply_1.7-4 abind_1.4-8 Rtsne_0.17 GenomicRanges_1.62.1
[25] BiocGenerics_0.56.0 pracma_2.4.6 tensorA_0.36.2.1 IRanges_2.44.0
[29] S4Vectors_0.48.0 ggrepel_0.9.6 irlba_2.3.7 listenv_0.10.1
[33] spatstat.utils_3.2-1 goftest_1.2-3 RSpectra_0.16-2 spatstat.random_3.4-4
[37] fitdistrplus_1.2-6 parallelly_1.46.1 codetools_0.2-20 DelayedArray_0.36.0
[41] RApiSerialize_0.1.4 tidyselect_1.2.1 farver_2.1.2 matrixStats_1.5.0
[45] stats4_4.5.2 spatstat.explore_3.7-0 Seqinfo_1.0.0 jsonlite_2.0.0
[49] ks_1.15.1 progressr_0.18.0 ggridges_0.5.7 survival_3.8-3
[53] tools_4.5.2 pak_0.9.2 ica_1.0-3 Rcpp_1.1.1
[57] glue_1.8.0 gridExtra_2.3 SparseArray_1.10.8 xfun_0.57
[61] cmdstanr_0.9.0 MatrixGenerics_1.22.0 distributional_0.6.0 withr_3.0.2
[65] fastmap_1.2.0 callr_3.7.6 digest_0.6.39 timechange_0.4.0
[69] R6_2.6.1 mime_0.13 scattermore_1.2 tensor_1.5.1
[73] dichromat_2.0-0.1 spatstat.data_3.1-9 utf8_1.2.6 generics_0.1.4
[77] httr_1.4.8 htmlwidgets_1.6.4 S4Arrays_1.10.1 uwot_0.2.4
[81] pkgconfig_2.0.3 gtable_0.3.6 lmtest_0.9-40 S7_0.2.1
[85] SingleCellExperiment_1.32.0 XVector_0.48.0 htmltools_0.5.9 dotCall64_1.2
[89] scales_1.4.0 Biobase_2.68.0 png_0.1-8 posterior_1.6.1
[93] spatstat.univar_3.1-6 knitr_1.51 rstudioapi_0.18.0 tzdb_0.5.0
[97] reshape2_1.4.5 checkmate_2.3.4 nlme_3.1-168 zoo_1.8-15
[101] KernSmooth_2.23-26 parallel_4.5.2 miniUI_0.1.2 vipor_0.4.7
[105] ggrastr_1.0.2 pillar_1.11.1 grid_4.5.2 vctrs_0.7.2
[109] RANN_2.6.2 promises_1.5.0 stringfish_0.18.0 xtable_1.8-4
[113] cluster_2.1.8.1 beeswarm_0.4.0 evaluate_1.0.5 mvtnorm_1.3-3
[117] cli_3.6.5 compiler_4.5.2 crayon_1.5.3 rlang_1.1.7
[121] future.apply_1.20.2 labeling_0.4.3 mclust_6.1.2 ps_1.9.1
[125] fs_1.6.6 plyr_1.8.9 ggbeeswarm_0.7.3 stringi_1.8.7
[129] viridisLite_0.4.3 deldir_2.0-4 lazyeval_0.2.2 spatstat.geom_3.7-0
[133] Matrix_1.7-4 RcppHNSW_0.6.0 hms_1.1.4 bit64_4.6.0-1
[137] future_1.70.0 shiny_1.12.1 SummarizedExperiment_1.40.0 ROCR_1.0-12
[141] igraph_2.2.2 RcppParallel_5.1.11-1 bit_4.6.0

[shihsama]

Hi,

Additionally, I tested a local patch in R/sccomp_remove_outliers.R.

Specifically, I added deterministic seeds to the two generate_quantities() calls used for posterior predictive outlier detection.

For the first generated-quantities call (rng), I added:

seed = as.integer(mcmc_seed + 1L)

For the second generated-quantities call (rng2), I added:

seed = as.integer(mcmc_seed + 2L)

Conceptually, the modified calls are:

rng <- mod_rng |>
  sample_safe(
    generate_quantities_fx,
    attr(.estimate, "fit")$draws(format = "matrix"),
    data = ...,
    parallel_chains = ...,
    threads_per_chain = cores,
    sig_figs = sig_figs,
    seed = as.integer(mcmc_seed + 1L)
  )

and:

rng2 <- mod_rng |>
  sample_safe(
    generate_quantities_fx,
    fit2$draws(format = "matrix"),
    data = ...,
    parallel_chains = ...,
    threads_per_chain = cores,
    sig_figs = sig_figs,
    seed = as.integer(mcmc_seed + 2L)
  )

After reinstalling the locally patched package, I repeated the previous reproducibility test using:

• the exact same sccomp_estimate() object;
• cores = 1;
• the same mcmc_seed in both sccomp_remove_outliers() calls.

The following comparisons now return TRUE:

identical(
  attr(out1, "model_input")$truncation_not_idx,
  attr(out2, "model_input")$truncation_not_idx
)

identical(
  attr(out1, "outliers"),
  attr(out2, "outliers")
)

The final parameter estimates were also identical between the two runs.

This local test suggests that explicitly propagating deterministic seeds to both generated-quantities steps resolves the reproducibility issue.

I used two deterministic offsets (mcmc_seed + 1L and mcmc_seed + 2L) so that the two posterior predictive stages use different but reproducible random-number streams.

[stemangiola]

Thanks @shihsama for the well done issue. I'm currently travelling. Do you feel like opening a PR?

[stemangiola]

@shihsama, could you please run your tests for the PR above? It should have solved the issue.

[shihsama]

@shihsama, could you please run your tests for the PR above? It should have solved the issue.

Hi!
Apologies for the delayed reply. I have tested your modified version (#280) and it solved perfectly.
Thanks for your efforts! Great trip!