Fit the spiDE negative binomial model

Builds the neighbourhood-interaction design for each niche bandwidth and fits a per-gene negative binomial GLM over all genes using the SpaNorm fitNB engine. The heavy fit uses fitNB's compute backend (CPU or GPU); the counts Y may be dense, sparse, or a DelayedArray and is not densified up front. The returned object holds one SpiDEFit per bandwidth. Neighbourhood effects are tested separately with testSpiDE().

Usage

fitSpiDE(spe, condition = NULL, ...)

# S4 method for class 'ANY'
fitSpiDE(
  spe,
  condition = NULL,
  index = NULL,
  niche = NULL,
  covariates = character(),
  sigma = NULL,
  assay = "counts",
  cell_type = "cell_type",
  sample_id = "sample_id",
  random = c("intercept", "none", "slope"),
  winsor = 4,
  lambda.a = 0,
  backend = c("auto", "cpu", "gpu"),
  name = "Niche",
  re.maxit = 2L,
  re.tol = 0.001,
  tau2.init = 1,
  re.prop = 1,
  re.maxit.psi = 1L,
  re.min.cells = 100L,
  df.method = c("satterthwaite", "between"),
  re.celltype = TRUE,
  converge = TRUE,
  converge.maxit = 50L,
  converge.tol = 1e-08,
  block.size = NULL,
  BPPARAM = BiocParallel::SerialParam(),
  verbose = TRUE,
  ...
)

Arguments

spe

a SpatialExperiment with niche reducedDims (see buildNiches()).

condition

a character, the colData column of the tested condition (must have exactly two levels), or NULL for a condition-free (niche-only) analysis. With NULL the condition terms are dropped from the design and the two-way CellType:niche interactions become the tested effects: within index cell type c, how expression changes with the local density of niche cell type n. The results(type = "celltype") and results(type = "patient") tables are empty in that mode, there being no condition to contrast.

...

further arguments forwarded to fitNB.

index, niche

character vectors restricting the index / niche cell types considered (NULL = all).

covariates

a character vector of nuisance colData columns to adjust for (e.g. library size, age, sex).

sigma

a numeric vector of bandwidths to fit; NULL (default) uses every Niche<sigma> reducedDim present.

assay

a character, the counts assay to model.

cell_type

a character, the colData column of cell type labels.

sample_id

a character, the colData column identifying samples (patients); used only when random != "none".

random

one of "intercept" (the default), "none" or "slope". Adds patient-level random effects (implemented as ridge-penalised design columns) to correct anti-conservative inference caused by cell-level pseudo-replication.

The default changed from "none" to "intercept". "none" forms Wald standard errors from cell-level information, which treats every cell as an independent replicate of a patient-level contrast. On a null with per-sample intercepts it rejects at \(\approx 0.71\) against a nominal 0.05, worsening as cell counts become imbalanced, where the random intercept holds \(\approx 0.04\) against a calibrated pseudobulk reference of \(\approx 0.04\). "none" is retained for back-compatibility and for reproducing pre-correction results, but it should not be used for inference. "intercept" adds a per-sample random intercept; "slope" additionally adds per-sample random slopes on the niche covariates, which protects the response x niche tests when the niche-slope varies between samples. Note that "slope" estimates an extra variance component (tau2 for SampleSlope) that is collinear with the tested fixed effect and, with few samples, is less stable than the intercept-only variance; prefer "intercept" at small S, and "slope" when between-sample niche-slope variation is expected. The fit is also stochastic (fitNB subsamples cells for the dispersion estimate), so set a seed for reproducible variance components. See the mixed-effects and simulation vignettes.

In a condition-free analysis (condition = NULL) the random-slope block sits on exactly the CellType:niche columns being tested, so "slope" is the natural correction there when between-sample variation in niche slopes is plausible. Be aware of a limitation specific to that mode: the tested slope is a within-sample contrast on a spatially autocorrelated covariate, and spiDE does not model spatial autocorrelation. Neighbouring cells are therefore not independent replicates of the slope, and neither random-effect structure can recover that – on a null fixture with a planted per-sample intercept, the fixed-effects fit made 37 calls, "intercept" 5 and "slope" 5 (see longtests/testthat/test-nicheOnly-mixed.R). Random effects remove most of the inflation but niche mode remains mildly anti-conservative; treat borderline calls with corresponding caution.

winsor, lambda.a

fitting parameters forwarded to fitNB (coefficient winsorisation and the base ridge penalty on the fixed columns).

backend

a character, the fitNB compute backend ("auto", "cpu", or "gpu").

name

a character, the niche reducedDim prefix.

re.maxit, re.tol

iteration cap and relative tolerance for the random-effect variance-component (Schall/PQL) loop (used when random != "none"). The tolerance is a relative change on log(tau2).

The default was lowered from 10 to 2 on measurement: for random = "intercept" one iteration is indistinguishable from ten on null type-I error (to three decimal places) and on tau2 (to two), because the loop typically converges in a couple of steps. It also largely dissolves a hazard of the larger cap – tau2 can enter a 2-cycle, making the result depend on the parity of re.maxit.

This evidence covers the intercept model only. Under random = "slope" the slope variance component decays monotonically across all ten iterations without meeting re.tol, so a cap of 2 leaves it far from where 10 leaves it. Pass re.maxit = 10 for slope fits.

tau2.init

initial random-effect variance component.

re.prop

the cell-subsampling proportion used to speed up the variance-component (PQL) loop, sampled per cell type within each sample (random != "none"). For a stratum of n cells, min(n, max(ceil(re.prop * n), re.min.cells)) are used. 1 (the default) disables subsampling (all cells, fully reproducible). The final fit that feeds inference always uses all cells; only the shared tau2 estimate is affected. No seed is set internally — set one externally for reproducibility of re.prop < 1 runs. Lowering re.prop trades accuracy for speed, and the trade is worse than it looks: a replicate study on real data (vignettes/spiDE-mixed- benchmark.Rmd) found that subsampling noise in tau2 does not shrink as re.prop rises from 0.2 to 0.8 (it stays comparable to or larger than genuine between-patient variation), and — more importantly — tau2 is systematically biased downward at every re.prop < 1 tested, an attenuation that averaging replicates cannot fix, only shrinking as re.prop approaches 1. Using re.prop < 1 is therefore rarely advised; only do so when the variance component's absolute scale doesn't matter (e.g. a quick feasibility check) and treat its tau2 as a lower bound, not a point estimate.

re.maxit.psi

dispersion iterations for the inner PQL loop fits (the final all-cell fit always uses full dispersion). 1 (default) skips the redundant re-estimation of the barely-moving dispersion each iteration.

re.min.cells

the per-stratum floor for re.prop subsampling.

df.method

one of "satterthwaite" (default) or "between"; only used when random != "none". "satterthwaite" derives a separate df per tested column from the shared variance-component fit, distinguishing between-sample contrasts (Response: small df, close to "between") from within-sample contrasts (ResponseNiche: larger df, more power) rather than applying S - 2 to both; @df is then a named per-column vector. "between" tests every Response/ResponseNiche coefficient against the same scalar between-sample reference df (S - 2), the original back-compatible behaviour, and @df is a scalar.

The default changed to "satterthwaite" after the benchmark study (research/) measured both arms on identically seeded data: "between" is severely over-conservative when samples are few (null type-I \(\approx 0.001\) at \(S = 4\) against a nominal 0.05, with correspondingly near-zero power), while "satterthwaite" holds type-I in \(0.042\)-\(0.065\) across the whole sampled range and gains \(\approx 0.10\) mean TPR. The trade is a mild liberal drift at larger \(S\) (worst measured \(\approx 0.065\)); use "between" when strict conservatism matters more than power, or for back-compatibility. Ignored when random == "none".

In a condition-free analysis (condition = NULL) the "between" reference df changes, because the tested CellType:niche slope is a within-sample contrast rather than a between-condition one: it is ncells - p_fixed under random = "intercept" (cells are the replicates) and S - 1 under random = "slope" (the per-sample random slopes sit on the tested columns, moving the contrast into the between-sample stratum). The name "between" is therefore a misnomer in the intercept case; it is retained for back-compatibility. "satterthwaite" computes this distinction from the fitted variance components and is preferred.

re.celltype

logical; when random != "none", add a nested (sample x cell type) random intercept alongside the per-sample one. Default TRUE. Without it the tested niche slopes are estimated from the total covariance of niche density and expression within an index cell type, so they also carry the between-sample composition effect: samples whose index cells sit in denser niche surroundings also differ in mean expression there. That is a patient-level association with S units, not neighbourhood-dependent differential expression, and a shuffle null that permutes within (sample, cell type) preserves it – which is why real data and such a null were indistinguishable on the YTMA cohort. With the block present every niche slope is a within-group slope and the shuffle null is calibrated in every expression band. Set FALSE to reproduce pre-correction fits. Ignored when random = "none".

converge

logical; after fitNB returns, converge each gene to its own penalised negative-binomial optimum and re-estimate its dispersion there. Default TRUE. fitNB fits every gene in one IRLS loop with a single gene-averaged cell weight vector and an aggregate convergence criterion, which is what makes a whole-transcriptome fit affordable; for bright, cell-type-restricted genes it stops 1-4 standard errors short of that gene's own optimum, with a dispersion about 1.6 times too large. The stage is per-gene, so it is blocked and parallelised over block.size / BPPARAM exactly as inference is. Note that it replaces edgeR's cross-gene moderated dispersion with a per-gene profile maximum-likelihood dispersion at the converged mean: with many cells per gene that is well determined, but it is a deliberate departure from fitNB's moderation. Set FALSE to reproduce pre-correction fits.

converge.maxit, converge.tol

iteration cap and relative log-likelihood tolerance for the per-gene convergence stage.

block.size

genes per block in the per-gene convergence stage (see testSpiDE() for the same argument at inference time). NULL fits every gene in one block.

BPPARAM

a BiocParallelParam. Used by the per-gene convergence stage (converge = TRUE), which is blocked over genes and dispatched with it; each worker densifies its own gene block, so peak memory scales with the number of workers. Absent an explicit block.size, one block per worker is used.

verbose

a logical, whether to print fitting progress.

Value

a SpiDEResults object (inference not yet computed).

Examples

data(toySpiDE)
spe <- toySpiDE
spe <- buildNiches(spe, sigma = 20)
fit <- fitSpiDE(spe, condition = "condition", sigma = 20, verbose = FALSE)
fit
#> SpiDEResults
#> Bandwidths (sigma): 20
#> Genes: 20
#> Condition: condition
#> Index cell types: A, B, C
#> Niche cell types: A, B, C
#> Tested: FALSE (0 rows in results table)

fit0 <- fitSpiDE(spe, condition = NULL, sigma = 20, random = "none",
                 verbose = FALSE)
fit0
#> SpiDEResults
#> Bandwidths (sigma): 20
#> Genes: 20
#> Mode: niche-only (no condition)
#> Index cell types: A, B, C
#> Niche cell types: A, B, C
#> Tested: FALSE (0 rows in results table)