Two-stage estimation: benchmarked against the published simulation study

What is being compared

condition is assigned per patient, so patients are the experimental units. [fitSpiDE()] fits one cell-level GLM and corrects the pseudo-replication with patient random effects (random = "none" | "intercept" | "slope"); [twoStageSpiDE()] instead estimates each patient’s niche slope and contrasts those slopes at the patient level — on this study’s largest design point, 30 numbers per (gene, index, niche) triplet rather than 9,000 cells. The model and its equations are documented in the companion vignette The spiDE model; this vignette measures the estimator.

Methods

Simulator and design points. The two-stage arm is run on the published simulation study’s own generator (research/R/sim_model.R): 10,000 genes (500 markers per cell type, 1,500 housekeeping, background remainder), 900 planted niche-responsive genes, niche bandwidth $\sigma = 50$ , and the same design grid — a null family (beta = 0) over samples $S \in \{4, 10, 16, 24, 30\}$ × {gradient, clustered}, and a power family sweeping samples, cells per sample, effect size and spatial layout one at a time around the fixed point ( $S = 24$ , 300 cells, $\beta = 2$ , gradient). Every design point is repeated over 40 replicates, and discoveries are re-gated at four FDR thresholds (0.01–0.2) from one fit per dataset.

Pairing. The two-stage harness (research/plasmode/benchmark_twostage_hpc.R) seeds each replicate from the published study’s generator, so the two-stage rows are computed on the same simulated datasets as the published fixed/intercept/slope rows — the comparison is paired, not two independent studies placed side by side. The published arms’ numbers are read from the published tables unchanged, never re-run.

Arms. fixed is random = "none"; intercept and slope are the mixed modes at their preferred Satterthwaite reference df. The two-stage arm shown throughout is stage1 = "ols" on raw counts — log-CPM regressed jointly on the centred niche columns per (sample, index) subset, unit weights — with patient pooling and the weighted moderated limma contrast exactly as twoStageSpiDE() ships them.

Measured quantities. Type-I error is the fraction of genes whose p-value for the single CellTypeA:conditionResponder:B triplet falls below 0.05 under the null. TPR (recall) and FDP score discoveries restricted to that same A:B triplet against the 900 planted genes. One methodological point matters when reading TPR/FDP: each arm is scored by the multiplicity machinery it actually ships. The GLM arms gate discoveries through spiDE’s three-level hierarchical BH cascade; the two-stage arm applies a plain BH across triplets (fdr.niche), because that is what twoStageSpiDE() reports. Raw power — the fraction of planted genes with uncorrected $p < 0.05$ — involves no multiplicity rule at all and is therefore the cleanest cross-arm comparator of whether the estimator sees the effect; it is read first below.

Why the main grid is the ols arm. stage1 = "spanorm" — the package default and the pipeline a real user runs — cannot be measured on the published simulator: it draws library sizes i.i.d. (rlnorm(n_cells, sdlog = 0.2), no spatial or biological structure). Across 30 null simulations, log library size is explained essentially not at all by structure — mean $R^2$ 0.0007 from cell type, 0.0062 from a smooth spatial surface and 0.0050 from sample (maxima 0.005, 0.022 and 0.010) — while at $\beta = 2$ the planted effect induces a correlation of $+0.687$ (range $+0.655$ to $+0.724$ over 8 simulations) between log library size and the spatial coordinate among affected cells, against $\approx 0$ in both control groups. The planted effect is therefore the only spatially structured contributor to sequencing depth, and SpaNorm’s library-size term removes it along with the depth it is designed to remove — a limitation of the simulator, not a property of the estimator, which leaves the GLM arms untouched only because they never normalise. Two follow-up runs close the gap and are reported in Stage-1 pipelines under the microscope below: the nb reference arm on the published simulator (samples sweep), and a structured library-size re-run of all three stage-1 arms, which restores what the i.i.d. simulator starves. Note the input difference when comparing across runs: in the structured-LS design the ols and nb arms analyse SpaNorm-adjusted expression (the realistic post-normalisation workflow), whereas on the published simulator they ran on raw counts — only the spanorm arm (raw counts plus the stored fit in both) is directly comparable across the two.

Adversarial scenarios. Three additional scenarios plant failure modes the published simulator cannot express, each on the gradient layout at $S = 10$ , 300 cells per sample, scored on 50 planted index-A marker genes:

maineffect — a niche-B slope shared by every patient plus a flat 2× baseline shift in responders. No condition × niche interaction exists, so any call on a planted gene is a false positive (the “removable interaction” leak that identity-scale estimators suffer).
depth — every gene’s counts scaled by a per-cell depth factor that tracks niche-B density, twice as strongly in responders; expression itself untouched. Tests library-size confounding.
attribution — a decoy cell type whose density is positively correlated with B’s, with the true effect running through B only. Tests whether the joint stage-1 fit attributes the effect to the right niche.

Null calibration

Type-I error at a nominal 0.05 under the null, versus the number of patient samples, for two layouts. Above the dashed line is anti-conservative.

Type-I error, nominal 0.05, 40 replicates per point.
layout	S	fixed	intercept	slope	twostage
clustered	4	0.0065	0.0489	0.0457	0.0491
clustered	10	0.0085	0.0628	0.0728	0.0474
clustered	16	0.0085	0.0633	0.0760	0.0436
clustered	24	0.0084	0.0639	0.0733	0.0399
clustered	30	0.0090	0.0634	0.0744	0.0371
gradient	4	0.0045	0.0360	0.0360	0.0457
gradient	10	0.0046	0.0518	0.0597	0.0438
gradient	16	0.0044	0.0532	0.0569	0.0384
gradient	24	0.0042	0.0542	0.0621	0.0319
gradient	30	0.0044	0.0531	0.0653	0.0283

The two-stage arm is the best-calibrated arm in this simulation study (on the real cohort the ordering reverses — the niche-shuffle companion report measures both estimators against their own nulls there). On the clustered layout it reads 0.049 at $S = 4$ , essentially on the nominal line, and drifts mildly conservative as patients accrue (to 0.037 at $S = 30$ ; the gradient layout is slightly more conservative still, ending at 0.028). The mixed GLM arms drift the other way — anti-conservative as $S$ grows, up to 0.064 (intercept) and 0.076 (slope) — and fixed is drastically conservative throughout (0.004–0.009). The two-stage drift is in the safe direction at every measured point: the arm never crosses nominal on this grid.

Raw power: does the estimator see the effect?

Fraction of the 900 planted genes with uncorrected p < 0.05, samples sweep, gradient layout. No FDR rule is involved, so the arms are directly comparable.

The two-stage estimator sees the planted effect at every sample size — raw power 0.13 at $S = 4$ rising to 0.48 at $S = 30$ — sitting above fixed (0.32 at $S = 30$ ) from $S = 10$ onward, but below the mixed GLM arms throughout (intercept reaches 0.66). Part of that gap is calibration rather than information: the null figure above shows intercept/slope buying detections with an inflated type-I error, while the two-stage test sits at or under nominal everywhere. The remainder is real — the GLM arms use every cell in one fit, while the two-stage estimator compresses each patient into one slope before testing.

Delivered discoveries: TPR and FDP at 5 % FDR

Recall (TPR) and realized false-discovery proportion at a nominal 5% FDR, samples sweep, gradient layout. Each arm is scored by the FDR machinery it ships.

Samples sweep, gradient layout, FDR 0.05.
S	method	tpr	fdp	n_disc
4	fixed	0.009	0.240	10.575
4	intercept	0.040	0.346	53.500
4	slope	0.041	0.354	54.475
4	twostage	0.000	0.000	0.000
10	fixed	0.029	0.104	29.175
10	intercept	0.068	0.072	65.575
10	slope	0.064	0.070	61.525
10	twostage	0.001	0.008	1.150
16	fixed	0.074	0.054	70.425
16	intercept	0.193	0.041	181.350
16	slope	0.198	0.044	186.000
16	twostage	0.020	0.035	18.825
24	fixed	0.146	0.040	136.750
24	intercept	0.315	0.034	293.075
24	slope	0.320	0.033	297.650
24	twostage	0.084	0.061	81.200
30	fixed	0.184	0.033	170.900
30	intercept	0.361	0.033	336.175
30	slope	0.370	0.038	346.175
30	twostage	0.136	0.075	133.075

Under each pipeline’s own delivered FDR control the ordering does not change: at $S = 30$ the two-stage arm recovers 0.14 of the planted genes against intercept’s 0.36, and below $S = 16$ it delivers almost nothing (TPR 0.001 at $S = 10$ ). The FDP panel shows the two failure modes living at opposite ends of the sweep. The GLM arms are worst where they are proudest: at $S = 4$ roughly a third of intercept’s discoveries are false (FDP 0.35), settling under nominal only from $S = 16$ . The two-stage arm makes no calls there at all, and once it does discover its FDP creeps modestly above nominal (0.061 at $S = 24$ , 0.075 at $S = 30$ , against intercept’s 0.033) — the flat BH over triplets is not immune at the margin. The practical reading: on this simulator the two-stage estimator trades delivered power for error control whose worst case is mild; it does not repeat its predecessor’s apparent dominance at large $S$ , which was measured on a different estimator with unmatched seeds and is superseded by this paired study.

Cells per sample: the estimator is cells-hungry

Recall at 5% FDR versus cells per sample (S = 24, gradient). Stage 1 estimates one slope per patient, so per-patient precision – cells – is what it consumes.

The samples sweep above makes the two-stage arm look patient-limited, but the cells sweep shows its real appetite: recall climbs from 0 at 100 cells per sample to 0.39 at 1,000 (intercept: 0.039 to 0.56). Every cell sharpens the one slope its patient contributes, so the estimator is cells-per-patient hungry, not only patient-hungry. At the smallest published grid point (25 cells per sample) the two-stage arm has no measurement at all — with min.cells = 30, no (sample, index) subset is estimable, which is the designed behaviour rather than a failure: subsets too small to support a slope are excluded instead of contributing an unstable one.

Spatial layouts and the negative control

Layout sweep at FDR 0.20 (S = 24, 300 cells).
layout	method	tpr	fdp	n_disc
clustered	fixed	0.208	0.071	201.750
clustered	intercept	0.424	0.160	457.150
clustered	slope	0.444	0.217	518.575
clustered	twostage	0.257	0.180	282.650
gradient	fixed	0.193	0.065	186.075
gradient	intercept	0.400	0.111	406.300
gradient	slope	0.420	0.158	454.150
gradient	twostage	0.184	0.184	203.800
multiniche	fixed	0.205	0.089	202.425
multiniche	intercept	0.399	0.120	411.575
multiniche	slope	0.421	0.179	472.775
multiniche	twostage	0.206	0.134	214.750
random	fixed	0.002	0.710	5.450
random	intercept	0.001	0.646	3.150
random	slope	0.001	0.669	3.800
random	twostage	0.000	0.075	0.075

On the structured layouts realized FDP is comparable across arms (0.07–0.22 at a nominal 0.20). The random layout — the no-true-structure negative control, where recall is correctly ~0 for every arm — separates them: the GLM arms’ FDP is 0.71 (fixed), 0.65 (intercept) and 0.67 (slope), against 0.075 for two-stage — a substantial specificity advantage. Read it with its denominator, though: mean discoveries per replicate there are 0.07 for two-stage (and single digits for the GLM arms), so the FDP ratio is formed over a handful of calls and is unstable replicate to replicate.

A niche-independent condition response

The ctresp scenario plants a condition main effect on the index cell type (beta_ct up to 2) with no niche dependence, so any niche call on a planted gene is spurious — the leak an identity-scale estimator would suffer, because a flat fold-change times a niche-correlated baseline looks like a slope.

At FDR 0.05 the two-stage arm’s spurious-call rate is 8.3^{-5} at worst across $\beta_{ct} \in \{0.5, 1, 2\}$ — effectively identical to intercept (0.0011), i.e. no leak. The log-scale stage-1 response is what buys this: a flat fold-change moves a gene’s intercept, not its slope.

Adversarial scenarios

Adversarial null scenarios (nominal 0.05), stage1 = “ols” on raw counts.
scenario	type1	type1_sd	type1_n
depth	0.096	0.042	40
maineffect	0.080	0.038	40

Attribution: raw p < 0.05 rate on genes whose effect runs through the true niche (B) versus the correlated decoy.
scenario	disc_true	disc_decoy	disc_true_n
attribution	0.126	0.054	40

Two of the three scenarios expose mild anti-conservatism in this (unnormalised) arm, and both readings are informative:

depth type-I is 0.096, roughly twice nominal. The scenario couples sequencing depth to niche density, condition-specifically — and depth confounding is exactly what normalisation exists to remove. The arm measured here regresses log-CPM of raw counts, whose global library-size scaling cannot absorb a spatially structured depth gradient. This is the empirical case for the stage1 = "spanorm" default on real data: its working response subtracts the fitted library-size field before the slope is estimated.
maineffect type-I is 0.08 against a nominal 0.05. The removable-interaction leak is largely — not completely — closed by the log-scale response; a residual excess of roughly half the nominal rate remains, and pinning down its mechanism is on the list for the normalised (spanorm) re-run.
Attribution works: planted genes are called through the true niche at 0.126 versus 0.054 through the correlated decoy — a 2.3-fold separation. The joint stage-1 fit forces the correlated columns to compete for the effect rather than each absorbing it fully, which is why the slopes are estimated jointly rather than one niche at a time.

The three stage-1 pipelines

The main grid above fixes stage1 = "ols". The structured-library-size simulator gives sequencing depth genuine spatial and biological structure, so normalisation has real work to do, and all three stage-1 arms were run on it with 40 replicates per point on one frozen package snapshot, sharing one per-sample SpaNorm fit per replicate: ols and nb analyse the SpaNorm-adjusted expression a real user would, spanorm takes the raw counts plus the stored fit and passes its technical linear predictor to the subset GLM as a fixed offset. (An earlier grid on the published i.i.d.-depth simulator is superseded: there the planted effect was the only structured contributor to depth, which starved the offset path and flattered ols; its rows remain in the tables under ls.model == "iid".)

Structured library sizes, samples sweep, 40 replicates per point. Left: raw power on the 900 planted genes by stage-1 pipeline. Right: the raw null rate on non-planted genes against the nominal 0.05 (dashed) – the calibration axis that separates the two fallbacks.

Three readings, in order of importance:

ols is the strongest stage-1 arm. Its null rate is flat in sample size (0.056 at $S = 4$ , 0.059 at $S = 30$ , nominal 0.05) at the highest raw power (0.45 at $S = 30$ against nb’s 0.29 and spanorm’s 0.27). A paired ablation on identical datasets also settled pool.psi = TRUE (type-I 0.067 pooled against 0.073 unpooled over 30 pairs, $p < 10^{-4}$ ).
The spanorm default is the weakest and most liberal arm. Its raw power scales with sample size (0.106 at $S = 4$ to 0.27 at $S = 30$ ) but trails both fallbacks, and its null rate 0.065–0.099 is the highest of the three (nb: 0.045–0.05). Prefer ols.
Little survives multiplicity at these sample sizes. Mean discoveries per replicate at $S = 30$ , FDR 0.05: 36.9 (ols), 37 (spanorm), 0.4 (nb), with TPR at or below 0.04. The raw p-values carry the signal; the estimator needs a pre-specified, restricted hypothesis.

On the real cohort

The simulation’s ordering — two-stage the best-calibrated arm — reverses on real data, and the reason is instructive. Under the niche-shuffle null of The real cohort report (the niche rows permuted within (sample, cell type), so that every call is false), the two-stage estimator’s statistic spread is far above 1 across all index types and inside its own null: the inflation is manufactured by thin (patient, index) subsets, about a dozen niche columns on a few dozen cells, and is present whether or not any niche association exists.

Statistic spread, sd(t), on the real cohort at bandwidth 30 (diamond and dashed line) against niche-shuffle nulls (points), for the two-stage estimator and the intercept GLM as it was then shipped, in all index types and in the two well-populated ones. The two-stage estimator sits inside its own null; its inflation across all types disappears in the populous ones.

Across all index types the two-stage spread is 1.36 on the real data against 1.35–1.38 shuffled; restricted to Tumor and Fibroblast it is 1.08, the same as the GLM’s. The calibration of the estimator on a cohort is governed by cells per (patient, index) subset: per index type its spread tracks the median subset size at $r = -0.84$ (Tumor, 388 cells per subset, at 1.05; Mast, 38 cells, at 1.65), against $r = -0.66$ for the patient count, and neither the bandwidth (1.363 at $\sigma = 30$ against 1.364 at 70) nor restricting to patients with more cores moves it. The lever is min.cells or the niche resolution, and a flat min.cells that reaches calibration is expensive: on this cohort 30 keeps 11 index types, 75 keeps 7, 100 keeps 5, and 100 deletes the entire T cell, dendritic, monocyte and mast compartment.

sd(t) against gene expression decile, Tumor + Fibroblast. Grey lines are shuffles; the coloured line is the real run. The rise with expression is reproduced by the null in both estimators: it is variance under-estimation scaling with precision, not biology.

Separately, @diagnostics$inclusion fires on this cohort: which patients contribute cells of a type is associated with their outcome in B cell (Fisher $p = 0.027$ ), dendritic cells (0.006) and monocytes (0.016). No threshold or variance correction fixes informative missingness; only Tumor and Fibroblast have complete inclusion and calibrated variance. The GLM, which fits all cells jointly so that thin types borrow strength, is the better-calibrated estimator across all index types on this cohort — and, since it is within-sample by construction, the two-stage estimator never carried the composition confound that the GLM needed a nested intercept to remove.

Reading this for a real study

What the estimator buys. On the simulation the ols and nb stage-1 arms are the best-calibrated methods in the study — at or under nominal at every null design point — and on the spatially random negative control the estimator makes essentially no false calls where the GLM arms accumulate them. The error structure is honest by construction: the contrast is formed over patients, which are the units the condition is assigned to.

What it costs. It is uniformly less powerful than the mixed GLM arms — at delivered 5% FDR it needs $S \gtrsim 16$ patients before it reports anything, and at $S = 30$ recovers about 0.38 of intercept’s recall — and it is cells-hungry, so its calibration degrades in exactly the thin subsets where the GLM borrows strength.

Multiplicity is untouched. Two-stage fixes the error structure, not the size of the hypothesis space. A full-panel space of about 1.8 million triplets buries real signal, and combining a gene’s mostly-null triplets dilutes a single true effect to no better than Bonferroni. Restrict index, niche and the gene set to a pre-specified hypothesis — genes rather than index types, since restricting index types raises BH’s threshold while the pathological genes stay — and prefer stage1 = "ols" with pool.psi = TRUE, on index types with complete inclusion and enough cells per patient.

And, as for every result in this suite: these are simulations sharing the package’s own generative model, so they cannot expose a mismatch common to simulator and method; the real-cohort null above is the check that can.

#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_4.0.3    spiDE_0.99.17    BiocStyle_2.40.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] SummarizedExperiment_1.42.0 gtable_0.3.6               
#>  [3] rjson_0.2.23                xfun_0.60                  
#>  [5] bslib_0.12.0                Biobase_2.72.0             
#>  [7] lattice_0.22-9              vctrs_0.7.3                
#>  [9] tools_4.6.1                 generics_0.1.4             
#> [11] stats4_4.6.1                parallel_4.6.1             
#> [13] tibble_3.3.1                pkgconfig_2.0.3            
#> [15] Matrix_1.7-5                RColorBrewer_1.1-3         
#> [17] S7_0.2.2                    desc_1.4.3                 
#> [19] S4Vectors_0.50.2            lifecycle_1.0.5            
#> [21] compiler_4.6.1              farver_2.1.2               
#> [23] textshaping_1.0.5           statmod_1.5.2              
#> [25] Seqinfo_1.2.0               codetools_0.2-20           
#> [27] htmltools_0.5.9             sass_0.4.10                
#> [29] yaml_2.3.12                 pkgdown_2.2.1              
#> [31] pillar_1.11.1               jquerylib_0.1.4            
#> [33] BiocParallel_1.46.0         SingleCellExperiment_1.34.0
#> [35] DelayedArray_0.38.2         cachem_1.1.0               
#> [37] limma_3.68.5                magick_2.9.1               
#> [39] abind_1.4-8                 tidyselect_1.2.1           
#> [41] digest_0.6.39               dplyr_1.2.1                
#> [43] bookdown_0.48               labeling_0.4.3             
#> [45] fastmap_1.2.0               grid_4.6.1                 
#> [47] cli_3.6.6                   SparseArray_1.12.2         
#> [49] magrittr_2.0.5              S4Arrays_1.12.0            
#> [51] withr_3.0.3                 scales_1.4.0               
#> [53] rmarkdown_2.32              XVector_0.52.0             
#> [55] matrixStats_1.5.0           otel_0.2.0                 
#> [57] ragg_1.5.2                  SpatialExperiment_1.22.0   
#> [59] evaluate_1.0.5              knitr_1.52                 
#> [61] GenomicRanges_1.64.0        IRanges_2.46.0             
#> [63] rlang_1.3.0                 Rcpp_1.1.2                 
#> [65] glue_1.8.1                  BiocManager_1.30.27        
#> [67] BiocGenerics_0.58.1         jsonlite_2.0.0             
#> [69] R6_2.6.1                    MatrixGenerics_1.24.0      
#> [71] systemfonts_1.3.2           fs_2.1.0

Dharmesh D. Bhuva