What was tried and rejected
Dharmesh D. Bhuva
Why this page exists
Every default in spiDE was chosen against alternatives,
and several of the alternatives were built, measured and removed. A
method that only shows what it kept invites the same alternatives to be
proposed again. This page records each one with the question it was
meant to answer, the measurement that decided it, and the number — so
that the record can be found without reading a notebook. The reports
that own the surviving choices are linked from each section.
A between-sample error stratum for the condition contrast
The idea. The condition is a patient-level label,
and a cell-level Wald standard error treats thousands of cells as
independent replicates of one patient-level contrast (it rejects 71% of
the time on a null with per-sample intercepts, against a nominal 5%).
Textbook split-plot logic says: code the samples as
contr.sum contrasts nested within condition, form the
between-sample mean square from them, and test the contrast against that
stratum instead of the residual.
What was right about it. Two things, both verified independently: the stratum must be cell-type specific (pooling it to sample level over-inflates the variance by the number of cell types), and it should be estimated from residual aggregates rather than extra design columns, which avoids a \(k(S-2)\) column blow-up and a double count that makes the design-based form fourteen times over-conservative.
Why it was rejected. Across an 18-cell factorial with a calibrated pseudobulk oracle and a known-broken negative control at every point, the construction is calibrated at exactly one design point — the one it was first tested on — and collapses to zero rejections as cells per sample grow. The decomposition against the true sampling distribution shows why:
The inflation the condition contrast requires (the true sampling SD of the coefficient over the cell-level SE) against the inflation the between-sample mean square applies, by cells per sample, at S = 12 with 30 replicates and fixed geometry. The required inflation is constant; the applied one grows as the square root of the cells per stratum. The scale is the wrong quantity, not a mis-scaled version of the right one.
| cells per sample | pseudobulk oracle | between-sample stratum | cell-level SE (control) |
|---|---|---|---|
| 60 | 0.0511 | 0.0198 | 0.6966 |
| 160 | 0.0472 | 0.0004 | 0.7091 |
| 400 | 0.0498 | 0.0000 | 0.7225 |
The random-effects model inflates through the penalised covariance
rather than by substituting a scale, and is calibrated throughout the
same factorial. The stratum was removed. What generalises: a negative
control at every design point, more than one design point before
declaring success, and measuring the pieces, not only the outcome. The
full account is fdr-ordering’s predecessor study in the
research repository (reports/data/,
calibration/).
Sample as a nested fixed effect
The same construction’s cousin: keep the nested
contr.sum block as ordinary fixed columns and test at the
cell-level standard error. Nesting keeps the condition contrast
estimable — a free indicator per sample would alias it away — but
nesting alone is not the correction: tested at the cell SE it is as
anti-conservative as ignoring samples entirely (type-I 0.76 balanced,
0.83 at ten-fold cell imbalance). The stratum above was the attempt to
supply the rest, and it failed for the reason above. The random
intercept supplies the between-sample stratum exactly, and the
Satterthwaite reference df gives the niche interactions the larger df
they earn; that is what ships.
Subsampling cells for the variance-component loop
The idea. The Schall/PQL loop that estimates the
shared variance components refits every gene per iteration, so fit the
inner iterations on a stratified cell subsample
(re.prop < 1), as fitNB already subsamples
cells for its dispersion step, and let a final all-cell fit supply the
coefficients.
Left: fit time against the subsampled fraction of cells (S = 24, 1,000 cells per sample, 40 replicates). Right: the correlation of the response-niche t-statistics with the all-cell fit. The saving is modest and the statistics move substantially.
Why it was rejected. On the real 55-patient cohort,
five refits at each of re.prop = 0.2, 0.5 and 0.8 gave a
mean variance component that rose monotonically with the fraction —
0.034, 0.052, 0.059 against 0.074 at re.prop = 1 — a
systematic downward attenuation, not scatter, so averaging replicates
cannot remove it; and the replicate noise itself (sd 0.015–0.017)
matched or exceeded the genuine between-patient signal the component is
meant to capture. re.prop = 1 is the default; the knob is
kept for feasibility checks, as a lower bound on the variance component.
The spiDE model vignette describes the loop.
Random slopes as the fix for the composition confound
The idea. If the tested niche slope carries a
between-patient component, give each patient its own random niche slope
(random = "slope") and test the fixed slope against
between-patient slope variation.
Why it does not do it. The composition confound is a
between-group mean effect: a patient’s index cells sit in a
denser niche on average and express the gene differently on average.
Per-sample slopes on the niche bases leave the group means untouched —
by Frisch–Waugh–Lovell, only a column that spans the group means
partials them out — so the confound stays on the fixed slope. Measured
on the cohort, the slope variance component converges to its floor
(there is no between-patient slope heterogeneity to absorb), and on the
simulation the random slope costs calibration at small \(S\) because a genuine effect seen in a few
responders is indistinguishable from slope noise. The nested (sample ×
cell type) intercept is the fix (The real cohort
report), and random = "slope" is kept for cohorts that do
carry slope heterogeneity, with re.maxit = 10 because its
component converges slowly.
A quasi-likelihood dispersion
The idea. The heavy null tail in the real cohort might be a dispersion estimated with too few degrees of freedom; an edgeR-v4 quasi-likelihood dispersion, with its own moderation, would be the principled replacement.
Why it was rejected. It was built and oracle-tested,
and it steepens the expression gradient of the null spread
rather than flattening it; the tail was not a dispersion problem at all
(The real cohort). The spec is retained in the package’s
design/specs/ with status refuted.
Filtering the brightest genes, and a circular six-fold
The idea. The null tail lives in the brightest genes; drop the genes whose null spread exceeds 1.3 and the false-call rate should fall.
False calls on a complete null at nominal BH 0.05, without a filter, with the hot-gene list defined on the same shuffle grid it is scored on, and with the list defined on other grids. The same-grid filter gives a six-fold reduction that is circular; the honest reduction is about 1.5-fold.
What survived. The filter is worth applying — it cuts complete-null false calls about 1.5-fold at \(\alpha = 0.05\) and 1.7-fold at 0.01 when the list is built on grids other than the one scored, and it costs no power — but it treats a symptom, and it is not the 6× that a single grid suggests: about two thirds of any one grid’s hot list is that grid’s own noise. It does not remove the need to fix the estimand, which The real cohort report does.
Restricting the index types to help the FDR
The idea. If only two index types are well calibrated, test only those, and the multiplicity burden falls with the hypothesis space.
| hypothesis space | tests | false calls, BH 0.01 | false calls, BH 0.05 |
|---|---|---|---|
| all 12 index x 13,348 genes | 1,815,328 | 116.5 | 223.5 |
| Tumor+Fibroblast only | 293,656 | 127.5 | 294.0 |
| drop sd(t)>1.3 genes | 1,798,736 | 6.5 | 36.0 |
| drop top-10% expressed genes | 1,633,768 | 29.5 | 88.5 |
| 500-gene panel x all index | 68,000 | 1.5 | 2.5 |
| 500-gene panel x Tumor+Fibroblast | 11,000 | 2.5 | 5.0 |
| panel, T+F, drop over-dispersed | 10,890 | 1.0 | 3.5 |
Why it makes the null worse. BH’s threshold is \(\alpha R / m\): shrinking \(m\) while the pathological genes stay raises the bar those genes are held to and increases the false calls. Restricting index types is a calibration-of-the-estimator argument, not a multiplicity one; restrict genes to a pre-specified hypothesis instead.
Cell-type-specific library-size factors
On the SpaNorm-adjusted assay, the per-cell-type library-size slopes spread from 0.46 to 0.72 in a biologically plausible order, which invited a size factor per cell type. On raw counts every type sits at 0.94–0.98, modelling them separately moves the null gradient from 1.149 to 1.144, and the two orderings correlate at Spearman 0.10: the adjusted-assay structure is an artefact of how the adjustment interacts with cell type. A single global library-size slope on raw counts is what ships (The real cohort, the substrate figures).
The design without CellType:condition
The design shipped before 0.99.7 reports about twice the niche recall of the current one on the simulation study. It is not a better test: without the column the flat cell-type response loads onto the niche slope and the slope’s standard error shrinks to that of a regression through the origin. The decomposition is in the Simulation study report’s calibration section, and the estimand argument in The spiDE model vignette.
Identity-scale niche densities, and quantile winsorisation
The idea. The effective niche enters the model as \(\log(1 + \text{count})\). Perhaps the raw density is the better covariate, and if outliers are the worry, winsorise them.
Paired comparison of the planted effect’s t-statistic under log1p and identity niche densities, on a clean fixture and on two fixtures with a discordant dense hotspot, by bandwidth (30 paired replicates). On clean data the transform barely matters; under the hotspot the identity scale loses more of the effect and can flip its sign.
What was decided. log1p stays: it costs
nothing detectable on well-behaved densities and protects the
sign of the estimate under a dense nest. Quantile winsorisation
is a no-op against a dense region by construction — the fence must come
from the middle of the distribution, not its rank order — and neither
transform nor cap recovers a hotspot: detection falls to 13–20% either
way, and the honest position is that those results are lost. The study
is niche-transform/ in the research repository, with the
package-side patch it was measured on.
Standardising the null per column rather than per gene
When the heavy tail was found, the natural repair was a scale per tested column (index × niche). It removes almost none of the excess — from about 60× at \(p < 10^{-6}\) to 37× — where a per-gene scale takes it to under 3×. The tail is a gene property; calibration is per gene (The real cohort).
Session info
#> R version 4.5.2 (2025-10-31)
#> Platform: x86_64-conda-linux-gnu
#> Running under: Red Hat Enterprise Linux 9.8 (Plow)
#>
#> Matrix products: default
#> BLAS/LAPACK: /home/uqdbhuva/miniconda3/envs/latest-r/lib/libopenblasp-r0.3.30.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Australia/Brisbane
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] patchwork_1.3.2 ggplot2_4.0.1
#> [3] SpatialExperiment_1.20.0 SingleCellExperiment_1.32.0
#> [5] SummarizedExperiment_1.40.0 Biobase_2.70.0
#> [7] GenomicRanges_1.62.0 Seqinfo_1.0.0
#> [9] IRanges_2.44.0 S4Vectors_0.48.0
#> [11] BiocGenerics_0.56.0 generics_0.1.4
#> [13] MatrixGenerics_1.22.0 matrixStats_1.5.0
#> [15] spiDE_0.99.18
#>
#> loaded via a namespace (and not attached):
#> [1] sass_0.4.10 SparseArray_1.10.1 lattice_0.22-9
#> [4] digest_0.6.38 magrittr_2.0.4 evaluate_1.0.5
#> [7] grid_4.5.2 RColorBrewer_1.1-3 fastmap_1.2.0
#> [10] jsonlite_2.0.0 Matrix_1.7-4 limma_3.66.0
#> [13] viridisLite_0.4.2 scales_1.4.0 codetools_0.2-20
#> [16] jquerylib_0.1.4 abind_1.4-8 cli_3.6.5
#> [19] rlang_1.1.6 XVector_0.50.0 withr_3.0.2
#> [22] cachem_1.1.0 DelayedArray_0.36.0 yaml_2.3.10
#> [25] poolr_1.2-0 S4Arrays_1.10.0 tools_4.5.2
#> [28] parallel_4.5.2 BiocParallel_1.44.0 dplyr_1.1.4
#> [31] mathjaxr_1.8-0 vctrs_0.6.5 R6_2.6.1
#> [34] lifecycle_1.0.4 magick_2.9.0 fs_1.6.6
#> [37] pkgconfig_2.0.3 pillar_1.11.1 bslib_0.9.0
#> [40] gtable_0.3.6 glue_1.8.0 Rcpp_1.1.2
#> [43] statmod_1.5.1 tidyselect_1.2.1 tibble_3.3.0
#> [46] xfun_0.54 knitr_1.50 farver_2.1.2
#> [49] rjson_0.2.23 htmltools_0.5.8.1 labeling_0.4.3
#> [52] rmarkdown_2.30 compiler_4.5.2 S7_0.2.1