Every default spiDE chose on evidence from a
real dataset was chosen here: a 55-patient CosMx lung cohort,
77,454 cells, 13,348 genes, condition Responder against Non-responder,
index cell types from Tumor (tens of thousands of cells) to Mast (a few
hundred). The simulation study on this site measures what the estimator
does when the truth is known; this report is about what it did when the
truth was not, and how a null built from the real data itself found a
defect that no simulation contained. It is written by argument, not by
date. The notebook that records the order of discovery is
fdr-ordering/FINDINGS.md in the research repository.
Six things are established, each with the figure that shows it:
Permuting the patient label — the classical null for a
between-patient contrast — leaves every patient’s real niche slopes
intact and randomises only who is a Responder. It tests the second stage
of the inference and can never detect a spurious slope. The
niche-shuffle null permutes the rows of the niche
matrix within (sample, cell type), so the true niche slope is zero by
construction and every call is false. It comes in two forms:
free (random rows) and block (a toroidal shift
of each sample’s niche field, which preserves local spatial smoothness
so that any architecture or field-of-view confound survives). Expression
is never shuffled across cells, because that would break each cell’s
pairing with its own library size and test the normalisation at the same
time.
The spread of the test statistic, sd(t), for the real cohort (diamond and dashed line) against its own shuffled nulls (points), for the two estimators, in all index types and in the two well-populated ones. A calibrated estimator sits at 1.0 (dotted). Under the design shipped before 0.99.17 the real data sit inside the null range: nothing distinguishes them.
The intercept GLM’s real-data spread was 1.053 against a shuffle maximum of 1.041; the two-stage estimator’s 1.36 against 1.35–1.38. The two forms of shuffle agree, which rules out unmodelled spatial structure as the driver. On this cohort, at this resolution, neither estimator detected niche-dependent differential expression — and, as the rest of this report shows, that verdict was an artefact of the design, not a fact about the tissue. A second fact is visible in the same figure and stands: across all index types the GLM is far better calibrated than the two-stage estimator, because it fits all cells jointly and thin cell types borrow strength; in the two populous types the two are equally, mildly liberal (the Two-stage estimator report takes this up).
Left: how far the null p-value distribution exceeds uniform at each threshold, on raw-count shuffles, with no standardisation, standardised per column (index x niche), per gene, or both. Standardising each gene by its own null spread removes almost the whole excess; standardising per column removes almost none of it. Right: the sorted per-gene and per-column null spreads — genes reach 2.14, columns stop at 1.21.
The heavy tail that defeats every FDR procedure is scale heterogeneity between genes, not between columns: standardising each gene’s statistic by its own null spread takes the excess at \(p < 10^{-6}\) from × to ×, where a per-column rescale reaches only ×. Most of the exceedances come from a small fraction of genes, and those genes are the brightest ones.
Median null spread by expression decile, on the SpaNorm-adjusted assay and on raw counts, with and without a library-size term. Every configuration inflates with expression; raw counts with a single library-size slope is the only arm centred on 1, and the best-calibrated substrate measured in this investigation.
Two facts about the substrate came out of this figure and were
adopted. The counts assay of the cohort object was a back-transform of
SpaNorm’s adjusted values, not counts: it reached totals well above any
cell’s true library size and destroyed sparsity. The raw counts still
existed and were restored. And per-cell sequencing depth was not in the
design at all — the cell-type intercepts cannot carry a per-cell
quantity — so a single global library-size slope was added through
covariates; the slope the data want on raw counts is 0.98,
a conventional offset. Cell-type-specific size factors buy nothing:
The per-cell-type library-size slope on the adjusted assay and on raw counts. The apparent cell-type structure of the adjusted assay’s slopes is an artefact of the adjustment; on raw counts every type sits at 0.94–0.98 and modelling them separately moves the null gradient from 1.149 to 1.144.
Twelve candidate causes on the variance side were measured and refuted, and the refutations are recorded so that none is re-run:
| candidate | why it was refuted |
|---|---|
the calculateMu() winsorisation clamp |
removing it does not move the tail |
| a patient-clustered sandwich variance | predicts a spread of 0.92–1.0 where 1.25–1.42 is observed |
| an edgeR-v4 quasi-likelihood dispersion | built and oracle-tested: it steepens the expression gradient |
| the ordering of the FDR cascade | every ordering fails on the same tail |
| within-gene variance misspecification | it makes the coefficient conservative, with the wrong sign for the excess |
| residual spatial autocorrelation | free and block shuffles agree |
| the shared-weight IRLS inefficiency | the converged fit lifts the tail rather than lowering it (below) |
| the extreme fits themselves | 35 genes; removing them leaves the band inflated |
| a cell-type-level dispersion | per-cell-type Pearson scale predicts 0.92–1.0 |
| a cell-level HC0 sandwich at the converged fit | the same 0.92–1.0 |
| a µ-tercile Pearson ratio | no marginal-variance quantity predicts the observed inflation |
| dropping the brightest genes | a real but 1.5× effect, and circular when defined on the grid it is scored on |
Observed against negative-binomial-predicted variance by fitted mean, split by the gene’s own abundance tercile. The variance model is wrong at low fitted means — but the cells with catastrophically low fitted means belong to the HIGH-expression genes, whose fitted means in the wrong cell types collapse toward zero. The variance function is wrong in every gene and with the wrong sign to explain the excess; the defect is in the mean, not the variance.
The multi-gene fit shares one gene-averaged weight vector, judges convergence on the aggregate log-likelihood and clamps coefficients across genes. For the brightest genes that is not their optimum:
Left: how far each gene’s production coefficients sit from its own penalised optimum, in units of its production standard error, against expression; the top-5% band sits 1–4 SE away, the 35 extreme fits further. Middle and right: neither that distance nor the inflated dispersion at the production point ranks the null spread within the band — the fit defect is real, and it is not the cause.
Converging each gene to its own optimum. Left: the production fit’s null spread collapses for bright genes (an inflated dispersion hiding inflation behind deflation); the converged fit lifts them above 1. Middle: no variant of dispersion or standard-error scaling restores a spread of 1 in the top 5%. Right: the extreme tail grows. Converging is necessary — the honest null variance of bright genes is higher — and it is not the cure.
Converging is now the default (converge = TRUE) because
the honest fit is the one to test, and because it sharpens real signal;
but on its own it makes the bright-gene null worse, which is
the observation that turned the investigation from the variance to the
estimand.
Hold a gene’s converged fit completely fixed and permute the niche rows within (sample, cell type). The permutation distribution of each tested \(t\) has exactly the spread a heteroscedasticity-robust sandwich predicts — and a per-column mean that runs to \(\pm 3\) and equals the real-data \(t\), column by column:
Left: for eight genes and 132 tested columns, the real-data t against the mean of t over 200 free shuffles with the fit held fixed (r = 0.61); the tested columns of Tumor, B cell and Fibroblast in red. Middle: the mean squared null t split into its variance and its squared mean, per index type — the excess over the sandwich prediction (crosses) is bias, and it is confined to the populous types. Right: centring the niche columns within (sample, cell type) removes it in every index type.
The mechanism is the design. The tested
cell type : condition : niche slope is estimated from the
covariance of niche density and expression among the cells of one index
type, and that covariance has two parts: within a sample, cells differ
in how dense their surroundings are (the effect spiDE measures), and
between samples, patients differ in how dense the surroundings
of their index cells are on average and also in how much of the gene
those cells express. The second is a patient-level association with 55
units. The shipped design had one ridge-penalised intercept per sample,
shared across cell types, and nothing per (sample × cell type); so the
between-patient part loaded onto the slope and was reported with a
cell-level standard error. The free shuffle permutes within
each (sample, cell type) group and so preserves every group’s mean niche
density — it preserves the association exactly, which is why the real
data and their null were indistinguishable: both carried the same
confound. The bias scales with the number of cells in the index type and
is largest for bright genes, exactly the pattern of the tail.
The Frisch–Waugh–Lovell theorem says what removes it. The coefficient
on a column in a weighted least-squares fit is estimated only from the
part of that column, and of the response, that the other columns cannot
explain; when the other columns include an indicator per (sample, cell
type), that residual is the within-group deviation, and the
between-group covariance is gone by construction. Each IRLS step is such
a fit, so it holds at the converged solution. The model vignette
(vignette("spiDE-model", package = "spiDE")) carries the
full statement.
Centring the niche-dependent columns within (sample, cell type) — the exact form of the theorem — was validated first, on all twelve raw grids:
Left: per-gene null spread against expression for the production fit, the converged fit, and the converged fit with the niche columns centred within (sample, cell type). Middle: by expression band, under block and free shuffles — the centred null is flat under free shuffles in every band; block shuffles keep a residual gradient in the brightest band, the spatially smooth component. Right: the extreme tail, 1,091 exceedances for the converged fit against 136 centred.
The package implements it as a ridge-penalised intercept per
non-empty (sample, cell type),
fitSpiDE(re.celltype = TRUE), carrying its own variance
component so the Satterthwaite degrees of freedom account for it. A
penalised block centres only partially — group means are shrunk toward
the pooled mean by the estimated variance component — so the packaged
fix was re-measured against the exact centring on the same grids:
The packaged fix on the real cohort at bandwidth 30 (spiDE 0.99.17, frozen snapshot). Top left: the median per-gene null spread by expression band under free and block shuffles (bars: the range over null grids), with the real data as diamonds. Top right: sd(t) per index type, real against the block-null range. Bottom: the number of |t| > 4.89 exceedances on each null grid, with the real data as the red line.
Under free shuffles the null is flat: the median
per-gene spread runs from 0.96 to 0.98 across the five expression bands,
with 0 exceedances of 101,508 in total over the free grids. Under
block shuffles the spatially smooth component survives in
the brightest genes (0.98 in the lowest band to 1.19 in the highest), so
calibration is read against block. The real data give
97 exceedances against a block-null range of 5–12, and
their spread sits above the null range in 8 of 12 index types.
That is the reversed verdict. The old confound was present in the real data and in every shuffle, so both inflated equally and matched; removing it from both leaves a difference. One hundred random control genes from below the top expression band sit inside the null range, so the difference is not a residual scale artefact.
Calibrated per gene — each gene’s real statistic divided by its own null spread over the block grids, then a threshold set so that the expected number of null exceedances is at most five percent of the calls — the base configuration gives 84 calls at bandwidth 30. Three artefact arms then asked whether those calls are stable:
Top: calls at per-gene empirical FDP <= 0.05 by arm and index type — the base configuration, five imaging covariates (cell area, nuclear area, DAPI intensity, aspect ratio, negative-probe background; all confounded with the niche covariate by construction), Plasma split out of the B cell compartment, and the base configuration at bandwidths 10, 50 and 70. Bottom: how many of each arm’s calls are also called in the base configuration.
The imaging covariates keep 78 of 84 calls; the Plasma split keeps 70 and removes every immunoglobulin call, so the earlier reading of plasma-cell activity was an artefact of the merge; bandwidths 10, 50 and 70 give 75, 63 and 61 calls, of which 1 triplet survives every bandwidth.
The bandwidth-robust triplets:
| gene | ct_index | ct_niche |
|---|---|---|
| NDRG1 | Tumor | Fibroblast |
Of these, 1 is also called under the imaging covariates and the Plasma split: NDRG1 in Tumor against a Fibroblast niche.
The estimator is calibrated in every arm — the null spread stays at 0.99–1.03 across bandwidths, compartment definitions and covariate sets — and the real grid exceeds its null in every arm. What is not stable is which genes are called. Segmentation spillover, confounded with the niche covariate by construction because both scale with neighbour density, was tested directly and is not supported at the cell level: genes called up are depleted, not enriched, in the neighbouring type. A raw discovery count from this cohort should not be quoted; what survives every perturbation can be.
What the nested intercept removes from the niche slope is a real
association in the data. compositionTest() asks the
question on 55 units: per (sample, index type) a pseudobulk profile of
the index cells and the mean niche density around them, then a
limma moderated \(t\)
across patients for the pooled association and for its difference
between conditions.
The leading composition associations on the cohort, one point per patient: pseudobulk log2-CPM of the gene in the index cells against the mean niche density around them, coloured by response. Every leading association is a neighbour type’s own marker appearing in the index cells’ pseudobulk in proportion to that neighbour’s density, with responders and non-responders on the same line.
The pooled association is strong — 1,462 calls at FDR 0.05 within (index, niche) blocks, 71 at a global FDR 0.05 — and its difference by condition essentially absent (6 and 2). The immunoglobulin-in-tumour signal the old design reported as niche-dependent DE is between-patient composition, plausibly with a segmentation-spillover component at the patient scale, uniform across response. It is reported beside the niche result, never as it.
| gene | ct_index | ct_niche | coef | t | fdr.global | n_samples |
|---|---|---|---|---|---|---|
| COL1A1 | Macrophage | Fibroblast | 2.447 | 7.33 | 0.00013 | 54 |
| MZB1 | Fibroblast | Plasma | 1.306 | 6.87 | 0.00228 | 55 |
| SPARC | Tumor | Fibroblast | 1.699 | 6.94 | 0.00280 | 55 |
| SLC6A8 | Tumor | Smooth.muscle.cell | 4.087 | 6.56 | 0.00776 | 55 |
| LYZ | NK | Macrophage | 2.775 | 5.80 | 0.00904 | 29 |
| LYZ | Mast | Macrophage | 3.079 | 5.57 | 0.02080 | 37 |
| COL1A1 | Monocyte | Fibroblast | 2.045 | 5.63 | 0.02080 | 51 |
| B3GALT2 | Smooth.muscle.cell | B.cell | -3.166 | -5.55 | 0.02080 | 37 |
The simulation study on this site places niche cells by the same potential in every sample, so it plants no composition effect; on it the nested block and the convergence step can only cost, and they do — a mild liberal shift on the null, decomposed by switch and by sample size in that report’s calibration section. Both facts hold: the defaults were chosen for a confound that real cohorts carry, because patients differ in composition.
The pre-specified final configuration is the 0.99.17 defaults on raw counts with the library-size slope, the five imaging covariates, Plasma split out of the B-cell compartment, bandwidth 30, calibrated per gene against five of its own block nulls. Two things about the calibration matter for reading its count. The threshold is set where the null exceedances, averaged over grids, fall to five percent of the real ones, and with only a handful of null triplets above \(|z| \approx 5\) that cut is discrete: a change of one or two null exceedances can move it by more than a unit of \(z\) and the count by an order of magnitude. Real and null exceedances at a fixed \(|z|\) are the stable comparison, so they are shown beside the calls.
Top: triplets above a fixed |z| at bandwidth 30, the real grid against the mean per block null, for the base configuration, the two single perturbations and their combination. The real grid exceeds its null by an order of magnitude in every arm at every threshold; the combination has fewer real exceedances than either perturbation alone. Bottom: the final configuration across bandwidths, calls at empirical FDP 0.05 and 0.10, with the real and null exceedances above |z| > 4.89 in the subtitle.
At bandwidth 30 the final configuration gives 9 calls at FDP 0.05 (|z| > 5.57) and 84 at FDP 0.10 (|z| > 4.25), against 84 and 146 for the base configuration. The drop is only partly the threshold: above a fixed |z| > 4.89 the final configuration has 38 real triplets to 2.4 per null grid, where the base has 63 to 1.6, so combining the imaging covariates with the Plasma split removes about 40% of the real exceedances (67 and 52 for either perturbation alone) while the null is unchanged. The real grid still exceeds its null 16-fold at that threshold. Across bandwidths 10 / 30 / 50 / 70 the final configuration calls 43 / 9 / 0 / 0 triplets at FDP 0.05 and 66 / 84 / 27 / 0 at FDP 0.10 (two null grids each beyond bandwidth 30, so those thresholds are the noisier ones).
At FDP 0.05 no triplet is called at more than 2 of the four bandwidths; those called at 2 are H2BC4 in B.cell against a Macrophage niche; H3C2 in B.cell against a Macrophage niche; NDRG1 in Tumor against a Fibroblast niche.
At FDP 0.10, 3 triplets are called at three or more bandwidths:
| gene | ct_index | ct_niche | n_bandwidths | bw10 | bw30 | bw50 | bw70 |
|---|---|---|---|---|---|---|---|
| H2BC12L | B.cell | Tumor | 3 | TRUE | TRUE | TRUE | FALSE |
| NDRG1 | Tumor | Fibroblast | 3 | TRUE | TRUE | TRUE | FALSE |
| SFTPA1 | Tumor | T.cell | 3 | TRUE | TRUE | TRUE | FALSE |
The all-gene run of the same configuration (13,348 genes, one real grid and five block nulls at bandwidth 30) is scored the same way when it completes, and is the run whose count would be quotable; the 769-gene panel here is enriched for the pathological tail.
ql-arms, fdr-ordering/FINDINGS.md,
2026-09-08): a moderated dispersion at the converged mean changes the
null by nothing, whereas a quasi-likelihood scale in the standard error
holds the nominal level at every sample size including \(S = 4\), at the no-switch design’s power,
and on this cohort calls the same triplets. Both the moderated
dispersion and the QL scale are the defaults from spiDE 0.99.18, with
the machinery in SpaNorm 1.7.10 on both backends. The nested block’s
reference degrees of freedom remain unmeasured.#> 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] S4Arrays_1.10.0 tools_4.5.2 parallel_4.5.2
#> [28] BiocParallel_1.44.0 dplyr_1.1.4 vctrs_0.6.5
#> [31] R6_2.6.1 lifecycle_1.0.4 magick_2.9.0
#> [34] fs_1.6.6 pkgconfig_2.0.3 pillar_1.11.1
#> [37] bslib_0.9.0 gtable_0.3.6 glue_1.8.0
#> [40] Rcpp_1.1.2 statmod_1.5.1 tidyselect_1.2.1
#> [43] tibble_3.3.0 xfun_0.54 knitr_1.50
#> [46] farver_2.1.2 rjson_0.2.23 htmltools_0.5.8.1
#> [49] labeling_0.4.3 rmarkdown_2.30 compiler_4.5.2
#> [52] S7_0.2.1