Calibration: reading lambda, and what the benchmarks measured

Dharmesh D. Bhuva

library(spiDE)
library(ggplot2)
source(system.file("spide_palette.R", package = "spiDE"))

Overview

This vignette has two halves. The first is a diagnostic: how to read the genomic inflation factor $\lambda$ from your own fit, per cell type, and what to do about it. The second is the evidence: every number this package quotes about its defaults — null calibration, power, and what was measured on a real cohort — collected in one place with pointers to the reports that carry the full tables (https://bhuvalab.github.io/spiDE-research/). The model those numbers are about is in vignette("spiDE-model").

Reading lambda

A spiDE fit tests thousands of genes against hundreds of celltype:condition:niche coefficients, and those coefficients are not equally well estimated. A tumour compartment with tens of thousands of cells spread over every patient supports a very different quality of inference from a dendritic-cell compartment with a few hundred cells concentrated in a handful of patients — yet both are reported in the same table, against the same FDR threshold. $\lambda$ detects that.

• $\lambda$ measures whether your test statistics follow the null distribution you assumed. It is a property of the test, not of the biology.
• In spiDE it varies within a single fit, across cell types, because cell types differ in how many patients contribute cells.
• $\lambda < 1$ costs power, not validity. A null result for a poorly calibrated cell type is weak evidence of absence.
• The right response is to restrict which coefficients you test, not to rescale the statistics and not to change the model.

What lambda is

For a one-degree-of-freedom test, the squared statistic $t^2$ is a $\chi^2_1$ variate under the null. The genomic inflation factor compares the observed median with the null’s:

$$\lambda = \frac{\mathrm{median}(t^2)}{\mathrm{median}(\chi^2_1)} = \frac{\mathrm{median}(t^2)}{0.4549}$$

The median is the important choice: a minority of true positives, however large, barely moves it, so $\lambda$ reads the bulk of the distribution, which is overwhelmingly nulls (Devlin and Roeder, 1999).

$\lambda$	reading
$= 1$	the observed null matches the theoretical one; p-values mean what they say
$> 1$	statistics too large: confounding or anti-conservatism, or genuinely widespread signal
$< 1$	statistics too small: the test is conservative

lambda <- function(t) stats::median(t^2, na.rm = TRUE) / stats::qchisq(0.5, 1)

Lambda below 1 does not mean “no signal”

Under a pure null with a calibrated test, $\lambda = 1$ exactly. Absence of biology gives 1, not less than 1:

z <- matrix(rnorm(200 * 50), 200, 50)
round(lambda(z), 3)
#> [1] 1.004

So $\lambda = 0.68$ is never a statement that a cell type is quiescent. It says the reference distribution is wrong for the statistic’s sampling behaviour — here, that the test is conservative. Wald tests are asymptotic: when few patients contribute cells, the variance is biased upward and $|t|$ comes out systematically too small. Because $|t|$ scales with $\sqrt\lambda$, a deflated test needs a proportionally larger effect to clear the same bar:

t_true <- stats::qnorm(1 - 0.001 / 2)          # |t| = 3.29 for p = 0.001
lam <- c(1.00, 0.90, 0.80, 0.70, 0.60, 0.50)
cost <- data.frame(lambda = lam,
                   t_observed = round(t_true * sqrt(lam), 2),
                   p_reported = signif(2 * stats::pnorm(-abs(t_true * sqrt(lam))), 2))
cost$times_worse <- round(cost$p_reported / 0.001, 1)
cost
#>   lambda t_observed p_reported times_worse
#> 1    1.0       3.29     0.0010         1.0
#> 2    0.9       3.12     0.0018         1.8
#> 3    0.8       2.94     0.0032         3.2
#> 4    0.7       2.75     0.0059         5.9
#> 5    0.6       2.55     0.0110        11.0
#> 6    0.5       2.33     0.0200        20.0

A cell type at $\lambda = 0.5$ reports a genuine $p = 0.001$ effect as about $p = 0.02$; after a genome-wide correction such a compartment is effectively unable to produce a discovery.

What drives it

The condition contrast is between patients, and the random intercept absorbs patient-level shifts, so what buys precision is the number of patients contributing usable cells of that type, not the raw cell count. $\lambda$ can depart from 1 in either direction, and the direction says which failure dominates:

• Many cells over few patients — pseudo-replication; statistics too large, $\lambda > 1$. This is what the mixed-effects correction fixes.
• Few cells and few patients — small-sample Wald conservatism; statistics too small, $\lambda < 1$.

Both resolve toward 1 as information accumulates. A third driver is not a variance at all: a between-patient composition association loading onto the niche slopes when the design has no (sample × cell type) intercept. It inflates $\lambda$ in exactly the pattern a real programme shift would — concentrated in the compartments with the most cells, largest for the brightest genes — and re.celltype = TRUE (the default) removes it; refit before interpreting an inflated compartment from an older fit. The measured form is under The real cohort below. A small measured cost of the defaults on a clean null also means a $\lambda$ a few percent above 1 is expected from the method and is not evidence of a data problem.

Diagnosing your own fit

Read $\lambda$per index cell type, never as one number for the fit — a global value is an abundance-weighted mixture that hides the heterogeneity.

lambda_by_celltype <- function(fit, type = "ResponseNiche") {
  cm <- as.data.frame(fit@coefmap)
  # `covtype` spans every design column, `t_stat` only the tested ones:
  # always match on column names
  cols <- intersect(colnames(fit@t_stat), cm$covariate[cm$type == type])
  idx  <- cm$index[match(cols, cm$covariate)]
  out  <- lapply(split(cols, idx), function(cc) {
    data.frame(n_col = length(cc), lambda = lambda(fit@t_stat[, cc, drop = FALSE]))
  })
  cbind(index = names(out), do.call(rbind, out))
}
# lambda_by_celltype(fits(res)[[1]])

Pair it with the number of patients contributing cells of each type, the quantity that predicts it:

tab <- table(spe$cell_type, spe$sample_id)
data.frame(index = rownames(tab), n_cells = rowSums(tab),
           n_patients = rowSums(tab > 0), n_pat_ge10 = rowSums(tab >= 10))

What to do about it

Do not rescale by lambda. Classical genomic control divides every statistic by $\lambda$; it assumes one inflation factor for all tests, whereas here $\lambda$ varies between cell types within a fit. A global correction deflates the calibrated compartments and leaves the sparse ones under-powered. The heterogeneity is the finding.

Do not reach for a bigger model. A random slope, a different penalty or a refit will not move $\lambda$ when the cause is small-sample conservatism. Two checks before blaming the model: a fixed-effects fit should be worse calibrated than the mixed fit, not better; and the dispersion multiplier at inference should have a median below 1, i.e. be shrinking standard errors. Conservatism at small $n$ is expected behaviour for a Wald test on a sparse negative binomial fit.

Restrict what you test. Stop testing coefficients that cannot support inference, filtering on patients contributing cells rather than on total cells. The under-powered tail stops consuming the multiple-testing budget, and the output stops presenting uninformative nulls as evidence. Report the threshold, and treat a null from an excluded compartment as weak evidence of absence.

When lambda exceeds 1. Because it is a median, a handful of strong hits cannot move it; a large fraction of genes must shift. Residual confounding affects null genes indiscriminately; a genuine programme-level shift concentrates in one compartment and one direction and reproduces in an independent contrast. Where one cell type shows $\lambda > 1$ and the rest of the fit sits at or below 1, a coordinated shift in that compartment is the more parsimonious reading — and it argues for interpreting at the level of the programme rather than gene by gene.

What the benchmarks measured

Each paragraph is one measured fact with the number that matters, in the order the defaults are introduced in the model vignette. The simulation study uses a realistic simulator (Gamma gene means, a mean–variance dispersion trend, library sizes of about 1,500–2,000 counts per cell, configurable spatial layouts), 40 replicates per design point, sample sizes $S$ from 4 to 30 and 300 cells per sample unless stated. Full tables: the Simulation study report on the benchmark site. All figures are for random = "intercept" with df.method = "satterthwaite" unless stated.

The random intercept

On a null with per-sample intercepts and no response effect, random = "none" rejects the response contrast at about 0.71 against a nominal 0.05, worse as cell counts become imbalanced across samples; the random intercept holds about 0.04, matching a calibrated pseudobulk reference. That is why "none" is not the default and should not be used for inference.

The reference degrees of freedom

df.method = "between" tests every coefficient against the scalar $S - 2$ and is severely over-conservative for the niche slopes when patients are few: null type-I about 0.001 at $S = 4$ against 0.05, with near-zero power. "satterthwaite" holds type-I in 0.042–0.065 over the whole sampled range and gains about 0.10 in recall, with a mild liberal drift at larger $S$ (worst 0.065). A Kenward–Roger correction of the covariance is the indicated refinement.

The small-sample regime

fitSpiDE(random = "intercept") is the most powerful calibrated estimator once samples are plentiful — recall 0.66 at $S = 30$ against 0.45 for the best two-stage arm — but it fails at small $S$ in a way the null table does not show: realised FDP 0.35 at $S = 4$ against a nominal 0.05, recovering to 0.04 by $S = 16$. Any sub-analysis that thins the patient count re-enters this regime.

The 0.99.17 defaults: nested intercept and per-gene convergence

The two defaults were adopted on real-cohort evidence (below). Measured as switch arms on the simulation study, which places niche cells by the same potential in every sample and therefore plants no composition effect, they can only cost there, and they do:

null type-I, mean over $S \ge 10$	at $S = 4$
design term alone: 0.058	0.036–0.049
convergence only: 0.066	0.086–0.095
nested block only: 0.068	0.037–0.043
both (the defaults): 0.071	0.085–0.096

(Per-cell standard error about 0.0016.) The two switches act at different sample sizes. Convergence’s cost is a small-$S$ effect — the per-gene dispersion gives up cross-gene moderation, which matters most when data are few — and fades onto the design-term arm by $S = 24$. The nested block is harmless at $S = 4$ and holds a constant $+0.010$ from $S = 16$, which points at its reference df or the estimation of its variance component. The pair gains recall (0.402 against 0.361 at $S = 30$) with slightly better FDP from $S = 10$ and a worse one at $S = 4$ (0.59 against 0.35). None of this is a reason to revert: the defaults were chosen for a confound real cohorts carry, because patients differ in composition. The two candidate refinements — a moderated dispersion at the converged mean, and the nested block’s reference df — are what this measurement will judge.

The CellType:condition term, and a design that looks more powerful

The design spiDE shipped before 0.99.7 had no CellType:condition column, and on the simulation study it reports far more niche calls: recall 0.36 against 0.17 at FDR 0.05 on the same dataset, concentrated in marker genes. That is not a better test. One power dataset ($S = 24$, $\beta = 2$) was fit under four package builds including one that differs from that design only by the added column; the build makes no difference (slopes correlate at 1.000), the column is the whole gap. Without it the niche slope exceeds the term-restored one by 0.25× the flat-response coefficient, gene by gene ($r = 0.88$ on planted genes, 0.92 on null genes) — the mean of the planted effect counted as niche dependence — which inflates the estimate by only about 10%; the $t$ is 1.5–1.6× larger because the standard error is 0.67–0.72×, a regression through the origin. On flat-shift genes with no niche effect the same channel produces fabricated niche calls; the recall saturates at about 0.5 from $\beta = 2$; and on the real cohort the term-restored design recovered 41% more triplets. Two-sided Cauchy combination gave more discoveries than one-sided Brown in every build, so sidedness is not involved. The recall of a design is comparable across designs only when they estimate the same quantity.

The combiner

Cauchy (two-sided) matches or beats Brown’s method on power while holding type-I error under correlation without estimating the correlation matrix; report: Combining niche p-values.

The two-stage estimator

On the simulation study its "ols" stage-1 path is the best-calibrated arm (null type-I 0.043–0.054, flat in $S$), with raw power 0.451 at $S = 30$ against 0.293 ("nb") and 0.274 ("spanorm", the most liberal at 0.055–0.071 and the least powerful); FDP 0.032 against 0.173 for "spanorm" at $S = 16$. Pooling the dispersion per sample (pool.psi = TRUE) is consistently better calibrated (0.067 against 0.073 paired over 30 datasets). On the real cohort the picture reverses for thin index types: its $\mathrm{sd}(t)$ tracks cells per (patient, index) subset at $r = -0.84$ — Tumor at 388 cells per subset is calibrated (1.05), Mast at 38 cells is not (1.65) — and the mixed GLM, which fits all cells jointly, is the better calibrated estimator across all index types (1.04 against 1.36). Dropout confounded with condition also fires there in B cell, DC and Monocyte. Report: Two-stage estimation.

The gene-set layer

spiGSEA(test = "competitive") is calibrated; the self-contained form is not (it called 20.6 of 208 null sets per replicate, all false) and is kept only to reproduce older scripts. Report: Gene-set inference.

The real cohort

The sharpest test in this package’s history, and the one to repeat before any claim of niche-dependent DE: a niche-shuffle null on a 55-patient CosMx cohort, permuting the rows of the niche matrix within (sample, cell type) so that the true slope is zero by construction, in free (random rows) and block (toroidal shift, preserving spatial smoothness) forms. Report: The real cohort.

• Under the design shipped before 0.99.17, the real data were indistinguishable from their null: $\mathrm{sd}(t)$ 1.100 against a shuffle maximum of 1.101. The heavy tail was per-gene — standardising each gene by its own null spread took the $p < 10^{-6}$ excess from 59.9× to 2.6× — and confined to the brightest genes in the most populous index types. Twelve variance-side explanations were refuted. The cause was the composition confound of the model vignette: with a gene’s converged fit held fixed, the permutation $t$ had the expected spread but a per-column mean up to $\pm 3$ that equalled the real-data $t$ column by column ($r = 0.66$); the shuffle preserves each group’s mean niche density and so preserves the confound, which is why real and null matched.
• With the nested intercept and convergence the free null is flat at 0.96–0.98 in every expression band with zero $|t| > 4.89$ exceedances of 101,508 where the old design produced 95; the block null keeps a 0.99 → 1.21 gradient in the ~100 brightest genes, the spatially smooth component, so calibrate against block. The real cohort now exceeds its null — 97 exceedances against a replicated null range of 5–11, with 100 random control genes inside the range — and gives about a hundred calls at a per-gene empirical FDP of 0.05.
• The count is not a set of findings. Splitting Plasma out of the B cell compartment keeps the count and loses every immunoglobulin call; five imaging covariates keep 84 of 104 and drop the keratin-in-immune calls; bandwidths 10, 50 and 70 share only four triplets. One triplet survives every perturbation. The patient-level composition test on the same cohort shows neighbour types’ own markers rising in the index cells’ pseudobulk with the neighbours’ density, uniformly across conditions — which is what the old design’s immunoglobulin-in-tumour calls were.
• Raw counts with a global library-size slope are the best-calibrated substrate measured: the expression gradient of the null went from 1.207 to 1.149 and the median null $\mathrm{sd}(t)$ from 1.055 to 1.008. A back-transformed assay is not counts, and cell-type-specific size factors buy nothing.

Session info

sessionInfo()
#> 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               fastmap_1.2.0              
#> [45] grid_4.6.1                  cli_3.6.6                  
#> [47] SparseArray_1.12.2          magrittr_2.0.5             
#> [49] S4Arrays_1.12.0             withr_3.0.3                
#> [51] scales_1.4.0                rmarkdown_2.32             
#> [53] XVector_0.52.0              matrixStats_1.5.0          
#> [55] otel_0.2.0                  ragg_1.5.2                 
#> [57] SpatialExperiment_1.22.0    evaluate_1.0.5             
#> [59] knitr_1.52                  GenomicRanges_1.64.0       
#> [61] IRanges_2.46.0              rlang_1.3.0                
#> [63] Rcpp_1.1.2                  glue_1.8.1                 
#> [65] BiocManager_1.30.27         BiocGenerics_0.58.1        
#> [67] jsonlite_2.0.0              R6_2.6.1                   
#> [69] MatrixGenerics_1.24.0       systemfonts_1.3.2          
#> [71] fs_2.1.0