The spiDE model

Overview

This vignette states the model behind spiDE once, from the top down: the quantity it estimates, the design that isolates that quantity, why cells are not replicates and what the random effects do about it, how each gene’s fit is converged, and how a Wald statistic becomes a row of results(). It is the reference for the other vignettes. The quickstart (vignette("spiDE")) shows the calls; the calibration vignette (vignette("spiDE-calibration")) holds every benchmark number quoted in this package; the two-stage estimator, a different route to the same question, has its own vignette (vignette("spiDE-twostage")).

One idea organises the whole design and is used throughout: the Frisch–Waugh–Lovell theorem. Every nuisance block in the design removes one component from the tested slope, and the model is best read as a sequence of such removals. It is stated in full in the section What each block does to the tested slope.

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

The question, and the quantity estimated

Within an index cell type, does a gene’s response to an experimental condition change with the local density of a niche cell type? That is one hypothesis per (gene, index type, niche type) triplet, and the quantity behind it is a slope: how much the condition effect on log expression changes per unit of niche density, among cells of the index type.

Three quantities have to be kept apart, and the design gives each its own column so that the tested slope carries only the third:

A flat cell-type response. The gene is up in every index cell of every Responder, regardless of neighbourhood. This is a between-patient contrast on cells of one type; results(type = "celltype") reports it.
A condition-independent niche slope. Expression in the index type tracks the niche’s density in every patient, Responder or not. This is a nuisance in the conditioned analysis and the tested quantity in niche-only mode (condition = NULL).
The niche-dependent response. The difference in that slope between conditions — the three-way cell type : condition : niche term.

Two more things are deliberately not part of the third quantity, and the model’s random effects exist to keep them out. Patients differ in baseline expression of every gene, and patients differ in composition: a patient whose index cells sit, on average, in a denser niche of some type may also express the gene differently in that type for reasons that have nothing to do with the neighbourhood of any single cell. Both are patient-level associations with one observation per patient. The second is a legitimate question of its own and compositionTest() asks it at the patient level; it is not neighbourhood-dependent expression and the niche slope must not carry it.

The effective niche

Every quantity in this section is computed within one sample at a time, entirely independently of every other sample. Let $\mathcal C_s$ be the cells of sample $s$ . For each cell $c \in \mathcal C_s$ the local intensity of cell type $t'$ is a kernel-smoothed count of that type’s cells around $c$ , with a Gaussian kernel $K_\sigma$ of bandwidth $\sigma$ :

$\Lambda_{c,t',s}(\sigma) \;=\; e_s(c) \sum_{c' \in \mathcal C_s} K_\sigma\!\left(d_{c,c'}\right)\,\mathbb{1}(t_{c'} = t'),$

where $d_{c,c'}$ is the planar distance between two cells of the same sample and $e_s(c)$ is Diggle’s edge correction for that sample’s own spatial window. buildNiches() never forms a distance matrix across samples: it splits the SpatialExperiment per sample and builds a fresh spatstat point pattern, window and density surface for each subset, so two cells from different samples that happen to share coordinates never interact. The expected log-count of type $t'$ in a disc of radius $\sigma$ is then

$N_{c,t',s}(\sigma) \;=\; \log\!\left(\pi\sigma^2\,\Lambda_{c,t',s}(\sigma) + 1\right),$

and stacking over types gives the effective niche of the cell, $N_{c,s}(\sigma)$ , a vector with one entry per cell type. Per-sample results are row-bound into one cells × types matrix per bandwidth, stored as reducedDim(spe, "Niche<sigma>"). Several bandwidths can be built at once; inference combines them later.

The model

Notation

$Y_{c,g}$ is the count of gene $g$ in cell $c$ ; $t_c$ its cell type and $s_c$ its sample. The condition $r_s \in \{0, 1\}$ is a property of the sample — constant across all cells of a sample, which drives everything below. $k$ indexes the (sample, cell type) group a cell belongs to. $S$ is the number of samples, $C$ the number of cells, $G$ the number of genes.

The negative binomial GLM, stated once

For one bandwidth $\sigma$ , gene $g$ has a negative binomial distribution with dispersion $\psi_g$ , $\mathrm{Var}(Y_{c,g}) = \mu_{c,g} + \psi_g\mu_{c,g}^2$ , and log-mean

$\begin{aligned} \log \mu_{c,g} \;=\;& \underbrace{\gamma_g^\top x_c \;+\; \zeta_{g,t_c} \;+\; \sum_{t'} N_{c,t'}\,\alpha_{g,t_c,t'}}_{\text{within-sample, fixed}} \;+\; \underbrace{\eta_{g,t_c}\, r_{s_c} \;+\; r_{s_c}\sum_{t'} N_{c,t'}\,\beta_{g,t_c,t'}}_{\text{response effects, fixed, tested}} \\[4pt] &+\; \underbrace{u_{g,s_c} \;+\; v_{g,s_c t_c} \;+\; \sum_{t'} N_{c,t'}\, w_{g,s_c,t_c,t'}}_{\text{patient-level random effects}} , \end{aligned}$

with

$u_{g,s} \sim \mathcal N(0, \tau^2_{\text{sample}}),\qquad v_{g,sk} \sim \mathcal N(0, \tau^2_{\text{sample}\times\text{type}}),\qquad w_{g,s,t,t'} \sim \mathcal N(0, \tau^2_{\text{slope}}).$

term	meaning	design tag
$\gamma_g^\top x_c$	cell-level covariates (`covariates =`), e.g. log library size	`Other`
$\zeta_{g,t}$	baseline expression of the gene in cell type $t$	`CellType`
$\alpha_{g,t,t'}$	condition-independent slope of expression in type $t$ on the density of type $t'$	`Niche`
$\eta_{g,t}$	the flat response of type $t$ : quantity 1 above	`ResponseCellType`
$\beta_{g,t,t'}$	the niche-dependent response: quantity 3, the tested slope	`ResponseNiche`
$u_{g,s}$	random intercept per sample	`Random` (`SampleInt`)
$v_{g,sk}$	random intercept per (sample, cell type) — `re.celltype = TRUE`, the default	`Random` (`SampleCellTypeInt`)
$w_{g,s,t,t'}$	random niche slopes per sample — `random = "slope"` only	`Random` (`SampleSlope`)

The variance components $\tau^2$ carry no gene index: one value of each is shared by every gene, estimated from all of them. The realised deviations $u_{g,s}$ , $v_{g,sk}$ vary by gene. random = "intercept" (the default) adds $u$ and, with re.celltype = TRUE, $v$ ; random = "slope" adds $w$ as well; random = "none" adds nothing and should not be used for inference (see Why cells are not replicates).

We test, per triplet,

$H_0:\ \beta_{g,t,t'} = 0 \qquad\text{against}\qquad H_1:\ \beta_{g,t,t'} \ne 0 .$

Cell-means coding

nicheDesign() assembles the fixed part as

~ 0 + <covariates> + CellType + CellType:condition
    + CellType:(niches) + CellType:condition:(niches) + niches

There is no bare condition main effect: the response contrast is carried by one CellType:condition column per cell type, each readable directly as that type’s flat response, and the three-way columns are one per (index type, niche type) with the index type’s own density dropped (a type’s density of itself is not a niche). The random-effect columns are ordinary design columns — an indicator per sample, an indicator per (sample, cell type), and, for slopes, sample indicators times the CellType:niche bases — tagged Random and identified by the ridge penalty described under Random effects as ridge. Every column carries a tag (fits(res)[[1]]$covtype), and the tested tag is ResponseNiche (or Niche in niche-only mode).

What each block does to the tested slope: Frisch–Waugh–Lovell

Consider a weighted least-squares fit of $y$ on two blocks of columns, $[X_1 \mid X_2]$ , with weights $D$ . The Frisch–Waugh–Lovell theorem says the coefficient on $X_2$ equals the coefficient from a regression of $M_1 y$ on $M_1 X_2$ , where

$M_1 \;=\; I - X_1\big(X_1^\top D X_1\big)^{-1} X_1^\top D$

removes from any vector its weighted projection on $X_1$ . In words: the coefficient on a column is estimated only from the part of that column, and of the response, that the other columns cannot explain. Every nuisance block in the design above is a choice of $X_1$ , and the theorem says exactly what it removes from $\beta$ .

The (sample × cell type) indicators. When $X_1$ is a set of group indicators, $M_1$ subtracts each group’s weighted mean: $(M_1 x)_c = x_c - \bar x_{k(c)}$ . With the nested block in the design, the tested slope is therefore estimated from the within-(sample, cell type) deviations of niche density and of expression, and the between-group part — patients whose index cells sit in denser niches also expressing the gene differently in that type — is removed exactly and assigned to the intercepts, where it belongs. Without the block, only a per-sample intercept is present, shared across cell types; nothing removes the between-sample covariance within an index type, so it loads onto $\beta$ and is reported with a standard error computed from $C$ cells rather than $S$ patients. A permutation of the niche rows within (sample, cell type) cannot detect this, because it preserves every group’s mean niche density and so preserves the confound. That is what the nested block was added for, and it is why random = "slope" is not a substitute: per-sample slopes on the niche bases leave the group means untouched.

The CellType:condition column. Now let $X_1$ be the flat response of the index type. $M_1$ removes the condition’s mean effect in that type from both the response and the niche column, so $\beta$ is the gradient of the response along the niche and not its level. Without the column the flat response has nowhere to go but the uncentred slope, and the slope’s standard error is that of a regression through the origin — smaller, because it is bought by assuming there is no level to estimate. A design without this column reports more niche calls; it is estimating a broader quantity under the niche label, not testing the same one better. The measurement is in the calibration vignette.

The IRLS fit. The theorem is about least squares, and the negative binomial fit is not least squares — but each step of the iteratively reweighted least-squares (and, later, Newton) iteration is a weighted least-squares problem in the working response, so the theorem holds step by step and hence at the converged fit, with $D$ the working weights $w_{c,g} = \mu_{c,g}/(1+\psi_g\mu_{c,g})$ .

One caveat, because it is honest. The nested block is not free: it is ridge-penalised (next section), and a penalised block does not centre completely. With penalty $\lambda$ on the group indicators, group $k$ ’s mean is removed by the factor $n_k \bar w_k / (n_k \bar w_k + \lambda)$ , where $n_k \bar w_k$ is the group’s total working weight — exact centring only as $\lambda = 1/\tau^2 \to 0$ (large estimated between-group variance) or as groups become large. On the real cohort that this was built for, the shrunken centring reproduced the exact-centring null within measurement; the numbers are in the calibration vignette.

Why cells are not replicates

The design has one row per cell. That is fine for the within-sample coefficients ( $\zeta$ , $\alpha$ ), but $\eta$ and $\beta$ are contrasts on the patient-level label $r_s$ : there are $S$ independent patients, yet a Wald information matrix that sums over $C$ cells makes their standard errors shrink like $1/\sqrt C$ instead of $1/\sqrt S$ . With $C \gg S$ the statistics are badly inflated. This is what random = "none" does, and the reason it is no longer the default.

The demonstration below simulates a null: every cell of a patient shares a per-patient shift, and there is no response effect at all, so every response test should be null.

sim_clustered <- function(n_samples = 10, n_per = 90, n_genes = 300,
                          field = 500, sd_patient = 0.7, seed = 1) {
  set.seed(seed)
  cts <- c("A", "B", "C"); gn <- sprintf("G%d", seq_len(n_genes))
  sids <- sprintf("S%d", seq_len(n_samples))
  cond <- rep(c("Responder", "Non-responder"), length.out = n_samples)
  names(cond) <- sids
  cells <- lapply(sids, function(sid) {
    x <- runif(n_per, 0, field); y <- runif(n_per, 0, field)
    ct <- ifelse(x > 0.6 * field & runif(n_per) < 0.7, "B",
                 sample(c("A", "C"), n_per, replace = TRUE))
    data.frame(sample_id = sid, condition = cond[[sid]],
               x = x, y = y, cell_type = ct)
  })
  cd <- do.call(rbind, cells); n <- nrow(cd)
  cd$cell_id <- sprintf("c%d", seq_len(n))
  base <- matrix(rnorm(n_genes * 3, 1.5, 0.6), n_genes, dimnames = list(gn, cts))
  lmu <- base[, cd$cell_type, drop = FALSE]
  # per-(gene, patient) random intercept: real clustering, NO response effect
  u <- matrix(rnorm(n_genes * n_samples, 0, sd_patient), n_genes,
              dimnames = list(gn, sids))
  lmu <- lmu + u[, cd$sample_id, drop = FALSE]
  counts <- matrix(rnbinom(length(lmu), mu = as.vector(exp(lmu)), size = 5),
                   n_genes, dimnames = list(gn, cd$cell_id))
  list(spe = SpatialExperiment(assays = list(counts = counts),
                               colData = S4Vectors::DataFrame(cd),
                               spatialCoords = as.matrix(cd[, c("x", "y")])),
       u = u)
}

sim_null <- sim_clustered()
spe <- buildNiches(sim_null$spe, sigma = 30)

300 genes rather than a handful, because $\tau^2$ is a single value shared across every gene: estimating a variance from a few genes is noisy in the same way as estimating it from a few observations.

A helper that fits one random mode and returns the response p-values and the fitted object:

fit_mode <- function(random) {
  # df.method = "between" (the scalar S - 2) is pinned here: this section is
  # about the correction itself, and the between reference keeps the demo
  # cheap. The default "satterthwaite" per-coefficient df has its own section.
  f <- fitSpiDE(spe, "condition", sigma = 30, random = random,
                df.method = "between", verbose = FALSE)
  res <- testSpiDE(f, spe = spe)          # populates the Wald inference slots
  ff <- fits(res)[[1]]
  # the response contrast is carried by the CellType:condition columns; match
  # on column NAMES -- `covtype` spans every design column, `t_stat` only the
  # tested ones
  ct <- as.character(ff$covtype)
  rc <- intersect(colnames(ff$t_stat),
                  ff$coefmap$covariate[ct == "ResponseCellType"])
  tt <- ff$t_stat[, rc, drop = FALSE]
  df <- if (is.null(ff$df)) Inf else if (length(ff$df) == 1L) ff$df else ff$df[rc]
  dfm <- matrix(rep(df, each = nrow(tt)), nrow = nrow(tt))
  list(fit = ff, p = 2 * pt(-abs(tt), dfm))
}
fixed <- fit_mode("none")
mean(fixed$p < 0.05)  # should be ~0.05 under the null
#> [1] 0.6233333

The fixed-effects fit rejects a large fraction of null tests at $\alpha = 0.05$ .

Why not add sample as a fixed effect?

The instinct from fixed-effects modelling is to add sample_id as a factor, the way one adds batch. It cannot work here, and the reason motivates the random effect:

set.seed(42)
S <- 12; n_per <- 12
toy_sample <- factor(rep(seq_len(S), each = n_per))
toy_resp   <- rep(rbinom(S, 1, 0.5), each = n_per)   # constant within a sample
true_tau2  <- 0.49
u_true     <- rnorm(S, 0, sqrt(true_tau2))
eta_true   <- 0.3 * toy_resp + u_true[as.integer(toy_sample)]
toy <- data.frame(y = rpois(length(eta_true), exp(eta_true)),
                  resp = toy_resp, sample_id = toy_sample)
coef(glm(y ~ sample_id + resp, family = poisson, data = toy))["resp"]
#> resp 
#>   NA

resp is NA because every cell of a sample shares the response label: resp is an exact linear combination of the sample dummies, and a design with a free column per sample has nothing left to estimate the contrast from. spiDE enforces the mirror image — checkSample() errors if a sample-constant covariate is passed with random != "none", since the random intercept makes it exactly as inestimable. A free per-sample parameter would also fit each sample perfectly and say nothing about how much patients differ in general, which is precisely the quantity the corrected test needs. A random effect sits between the two extremes: each patient gets its own deviation $u_s$ , but the deviations are draws from a shared $\mathcal N(0, \tau^2)$ , with $\tau^2$ estimated from the data and each $\hat u_s$ shrunk toward zero by an amount that reflects how reliable that patient’s own data is.

Random effects as ridge

The identity

One IRLS step for a GLM is a weighted least-squares regression of the working response $z_{c,g} = \eta_{c,g} + (Y_{c,g} - \mu_{c,g})/\mu_{c,g}$ on the design, with weights $D_g = \mathrm{diag}(w_{c,g})$ . Write the expanded design as $\tilde W = [\,W \mid Z\,]$ with $Z$ the random-effect columns. On this working scale, a Gaussian random effect on a block of columns is exactly an $L_2$ (ridge) penalty on those columns — an identity, not an approximation. The penalised objective

$\lVert z - \tilde W\theta \rVert^2_{D_g} + \theta^\top \Lambda\,\theta, \qquad \Lambda = \mathrm{blockdiag}\big(0_{\text{fixed}},\ \lambda_{\text{sample}} I,\ \lambda_{\text{sample}\times\text{type}} I,\ \lambda_{\text{slope}} I\big),$

has the closed-form solution $\hat\theta = (\tilde W^\top D_g \tilde W + \Lambda)^{-1}\tilde W^\top D_g z$ with

$\boxed{\ \lambda_k = 1/\tau^2_k\ } .$

Read it either way: a ridge penalty is the log-density of a Gaussian prior $b \sim \mathcal N(0, 1/\lambda)$ , so the penalised fit is the posterior mode under the random-effect distribution; or, the penalty pulls every $\hat u_s$ toward zero, which is what borrowing strength across patients means numerically. The penalised estimates of the random columns are the model’s BLUPs. This is the penalised quasi-likelihood construction (Breslow & Clayton, 1993; Schall, 1991) that lme4, mgcv and MASS::glmmPQL() build on, and it needs no change to the fitting engine: SpaNorm::fitNB accepts a per-column penalty lambda.a, set to $0$ on fixed columns and $1/\tau^2_k$ on random ones.

Estimating the variance components

The one thing the identity does not supply is $\tau^2$ . .fitNBmixed() gets it with two nested loops. The inner loop is fitNB at a fixed penalty: the penalised estimate of every gene at once, BLUPs included. The outer loop re-estimates each component by a method-of-moments update and refits:

$\hat\tau^2_k = \frac{\sum_g \lVert \hat b_{g,k}\rVert^2} {G\,\big(q_k - \lambda_k\,\mathrm{tr}\,C_{kk}\big)}, \qquad C = \big(\tilde W^\top \bar D\,\tilde W + \Lambda\big)^{-1},$

where $\hat b_{g,k}$ are gene $g$ ’s BLUPs in block $k$ , $q_k$ the block’s number of columns, $\bar D$ a gene-averaged working-weight matrix (one representative information matrix, since $\tau^2$ is shared), $C_{kk}$ the block’s penalised covariance, and $q_k - \lambda_k\,\mathrm{tr}\,C_{kk}$ its effective degrees of freedom — less than $q_k$ because shrunken BLUPs use up less than a parameter each. That denominator is what makes $\hat\tau^2$ an estimate of the true between-patient variance rather than the variance of the shrunken BLUPs. The loop stops when $\log\hat\tau^2$ moves by less than re.tol, or after re.maxit iterations (default 2: for the intercept model one iteration is indistinguishable from ten, and a small cap avoids a two-cycle the estimate can fall into; pass re.maxit = 10 for random = "slope", whose slope component converges slowly). re.prop < 1 runs the inner fits on a stratified cell subsample; the coefficients inference uses always come from a final fit on all cells, but the shared $\tau^2$ is biased downward by subsampling, so the default is re.prop = 1.

Watching it work

The same construction, on one gene of the toy data, with an independent implementation — MASS::glmmPQL():

library(MASS)
fit_pql <- glmmPQL(y ~ resp, random = ~ 1 | sample_id, family = poisson,
                   data = toy, verbose = FALSE)
summary(fit_pql)$tTable
#>                  Value Std.Error  DF    t-value   p-value
#> (Intercept) 0.33464279 0.6069414 132 0.55135935 0.5823195
#> resp        0.03611324 0.7018176  10 0.05145673 0.9599749

resp is now estimable, and its reference DF is $S - 2 = 10$ , the between-patient degrees of freedom. glmmPQL() exposes its outer-iteration cap as niter, so refitting with niter = 1, 2, … replays the trajectory .fitNBmixed()’s loop takes:

tau2_trace <- vapply(1:8, function(k) {
  f <- suppressWarnings(glmmPQL(y ~ resp, random = ~ 1 | sample_id,
                                family = poisson, data = toy, niter = k,
                                verbose = FALSE))
  as.numeric(nlme::VarCorr(f)[1, 1])
}, numeric(1))
ggplot(data.frame(iteration = seq_along(tau2_trace), tau2 = tau2_trace),
      aes(iteration, tau2)) +
  geom_hline(yintercept = true_tau2, linetype = 2, colour = spide_pal$ref) +
  geom_line(colour = spide_pal$accent) + geom_point(colour = spide_pal$accent) +
  labs(x = "outer iteration", y = expression(hat(tau)^2),
       title = "One gene's variance-component estimate")

The outer loop converging on the patient-level variance component for one gene. The dashed line is the simulated value.

It converges in a handful of iterations and does not land on the true value: with $S = 12$ patients, one gene’s $\hat\tau^2$ has real sampling noise. That is why spiDE pools the numerator across every gene; on the 300-gene null from earlier the pooled estimate sits close to the simulated $0.49$ :

mixed <- fit_mode("intercept")
mixed$fit$tau2[["SampleInt"]]
#> [1] 4.644419

Shrinkage is the other signature of a random effect. Compare each sample’s BLUP with the fully unpooled per-sample estimate from the saturated model y ~ sample_id - 1:

naive_u <- coef(glm(y ~ sample_id - 1, family = poisson, data = toy))
naive_u <- naive_u - mean(naive_u)
blup_u <- nlme::ranef(fit_pql)[[1]]
ord <- order(naive_u)
shrink_df <- data.frame(sample = factor(seq_len(S)[ord], levels = seq_len(S)[ord]),
                        naive = naive_u[ord], blup = blup_u[ord])
ggplot(shrink_df) +
  geom_segment(aes(x = sample, xend = sample, y = naive, yend = blup),
              arrow = arrow(length = unit(0.12, "cm")), colour = "grey55") +
  geom_hline(yintercept = 0, linetype = 2, colour = spide_pal$ref) +
  geom_point(aes(sample, naive, colour = "unpooled"), size = 2.2) +
  geom_point(aes(sample, blup, colour = "BLUP"), size = 2.2) +
  scale_colour_manual(
    values = c("unpooled" = unname(spide_pal$method[["fixed"]]),
              "BLUP" = unname(spide_pal$method[["intercept"]])), name = NULL) +
  labs(x = "sample (sorted by unpooled deviation)",
       y = "deviation from population mean",
       title = sprintf("var(unpooled) = %.2f, var(BLUP) = %.2f",
                       var(naive_u), var(blup_u))) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

Each sample’s unpooled deviation and its BLUP, connected by an arrow. Every arrow points toward zero, and the least reliable estimate is pulled in the most.

Shrinkage is strongest where a sample’s own data is least reliable relative to $\tau^2$ — the toy uses few cells and a low mean so that it is visible. A well-measured gene gets little help from the random effect and only a modest standard-error correction; a sparse gene, whose per-sample counts are noisy, is where the correction does most of its work.

Converging each gene

Why the shared fit is not at every gene’s optimum

SpaNorm::fitNB fits every gene in one IRLS loop. To make a whole-transcriptome fit affordable it shares one gene-averaged cell-weight vector across genes, judges step-halving and convergence on the aggregate log-likelihood, and clamps coefficients across genes. For most genes that is indistinguishable from the gene’s own optimum. For a bright, cell-type-restricted gene, whose working weights look nothing like the average, the aggregate criterion is met before that gene’s own score reaches zero, and the dispersion estimated at that point is off. converge = TRUE (the default) finishes the job per gene, and polishSpiDE() applies the same stage to an existing fit — made with converge = FALSE, or saved by an older version — without refitting. The two routes share one implementation.

Damped Newton on the gene’s own likelihood

For gene $g$ , with the design $\tilde W$ and penalty $\Lambda$ fixed, the objective is its own penalised negative binomial log-likelihood,

$\ell_g(\theta;\psi_g) \;=\; \sum_c \log \mathrm{NB}\!\left(Y_{c,g};\ \mu_{c,g} = e^{\tilde w_c^\top\theta},\ \psi_g\right) \;-\; \tfrac12\,\theta^\top\Lambda\,\theta ,$

whose score and (expected) information are

$s(\theta) = \tilde W^\top \frac{Y_g - \mu_g}{1 + \psi_g\mu_g} - \Lambda\theta, \qquad I(\theta) = \tilde W^\top \mathrm{diag}\!\left(\frac{\mu_g}{1+\psi_g\mu_g}\right)\tilde W + \Lambda .$

Each Newton step solves $I(\theta)\,\delta = s(\theta)$ and halves the step until $\ell_g$ increases; the iteration stops when the gain falls below tol relative to $|\ell_g|$ . A gene whose starting fit is degenerate — a fitted mean below $e^{-10}$ somewhere — restarts from the per-cell-type log means, because Newton from a degenerate point diverges. Each polished gene depends only on its own counts, so the stage is blocked over genes and parallelised (block.size, BPPARAM), exactly as inference is.

Profile maximum likelihood for the dispersion

fitNB estimates $\psi_g$ with edgeR’s empirical-Bayes moderation, which shrinks each gene’s dispersion toward a cross-gene trend. The polish replaces it with the gene’s own profile maximum-likelihood value. “Profile” means: hold the coefficients at their converged value $\hat\theta$ , so that $\hat\mu_{c,g}$ is fixed, and view the log-likelihood as a function of the dispersion alone,

$\ell_g(\psi \mid \hat\mu) \;=\; \sum_c \Big[ \log\Gamma\!\big(Y_{c,g} + \tfrac1\psi\big) - \log\Gamma\!\big(\tfrac1\psi\big) - \log\Gamma(Y_{c,g}+1) + Y_{c,g}\log\frac{\psi\hat\mu_{c,g}}{1+\psi\hat\mu_{c,g}} - \frac1\psi\log\big(1+\psi\hat\mu_{c,g}\big)\Big] ,$

then maximise it over $\log\psi$ on a bounded interval (psi.range, by default $[10^{-3}, 10^{3}]$ ) with a one-dimensional search (stats::optimize). The coefficients are then re-polished at the new $\psi_g$ , and the pair of steps is repeated once, so the returned $(\hat\theta, \hat\psi_g)$ is a joint optimum to the stated tolerance.

Three consequences follow from giving up moderation.

A gene that carries no information about overdispersion — all zero, or effectively Poisson — has a profile likelihood that is monotone in $\psi$ , so the search returns a bound, not an estimate. Such a gene keeps fitNB’s moderated value and is flagged psi_bound in @polish.
An unmoderated $\psi_g$ is no longer a safe scale for the standard error, so inference uses the Pearson working dispersion (next section) on every path.
With few samples the loss of shrinkage is measurable as a mild anti-conservatism that fades as samples accrue (calibration vignette). A moderated dispersion evaluated at the converged mean is the indicated refinement.

The polish also converges on the unclamped likelihood: fitNB’s winsorisation is a robustness device for a fit that shares one weight vector across genes; a per-gene converged fit does not need it, and inference evaluates the same unclamped mean the polish converged on. With converge = TRUE the winsor argument therefore only sets the Newton starting point. The negative binomial likelihood is undefined off the integers, so the stage requires integer counts and refuses a back-transformed assay rather than silently collapsing every dispersion to a bound.

Absorbing the nested block: the Schur complement

The nested block adds one column per non-empty (sample, cell type) group — several hundred on a cohort — and a naive Newton step would factor an information matrix of that size for every gene at every iteration. It does not have to. Partition the design into the dense columns $X$ and the group indicators $Z$ , so that the information is

$I(\theta) \;=\; \begin{pmatrix} A & B \\ B^\top & C \end{pmatrix}, \qquad A = X^\top D X + \Lambda_X,\quad B = X^\top D Z,\quad C = Z^\top D Z + \Lambda_Z ,$

with $D$ the working weights. Because every cell belongs to exactly one group, $Z^\top D Z$ is diagonal — entry $k$ is the total working weight of group $k$ — and so is $C$ once the penalty is added. The Newton step for the dense columns is then the $p \times p$ system

$S\,\delta_X = s_X - B\,C^{-1} s_Z, \qquad S = A - B\,C^{-1}B^\top ,$

and the group steps follow by back-substitution, $\delta_Z = C^{-1}\big(s_Z - B^\top\delta_X\big)$ . $S$ is the Schur complement of $C$ . Forming it costs one dense gram $X^\top D X$ plus group sums of the weighted columns, whatever the number of groups, which is what makes the per-gene polish affordable on a cohort.

The Schur complement is the Frisch–Waugh–Lovell theorem in matrix form: $S$ is precisely the information of the dense columns after the group indicators have been partialled out — the centring described earlier, done by algebra rather than by subtracting means. .newtonSolver() verifies the partition (each row of $Z$ sums to one) before relying on it, and .blockedInference() uses the same absorption to form each gene’s covariance.

What the stage reports

@polish is a per-gene table: iterations, psi_fitnb (the dispersion before polishing), restarted, capped, singular, psi_bound and polished. polished = FALSE marks a gene the Newton could not run on at all, which keeps the multi-gene fit’s estimate rather than substituting a worse one. Both flags are expected in small numbers; a large count is worth investigating before trusting that bandwidth. On bright genes the converged standard errors are larger than the unconverged ones — the unconverged dispersion, estimated off the optimum, was hiding inflation behind deflation — and the larger number is the honest one.

Inference

Covariance and dispersion

With the random effects and the converged fit in place, the Wald statistic for a tested coefficient uses three things.

The full penalised covariance, $\widehat{\mathrm{Var}}(\hat\theta) = \varphi_g\,[(\tilde W^\top D_g \tilde W + \Lambda)^{-1}]$ , with the tested block extracted. Because $Z$ absorbs the within-patient correlation, this block is on a between-patient scale — larger than the fixed-effects standard error, which is the point.
The Pearson working dispersion, $\varphi_g = \frac{1}{C - p}\sum_c \frac{(Y_{c,g} - \mu_{c,g})^2}{\mu_{c,g} + \psi_g\mu_{c,g}^2},$ the same quasi-likelihood device summary(glm(..., family = quasipoisson)) reports: trust the variance function, rescale it by the residual scatter this fit shows.
A $t$ reference distribution with a per-coefficient degrees of freedom, next.

Every per-gene matrix inversion is guarded: a gene whose information matrix is singular drops out as NA rather than aborting the run.

Reference degrees of freedom

The flat response $\eta_{g,t}$ is a pure between-patient contrast, so its reference df is close to $S - 2$ , the bookkeeping of a two-sample $t$ -test on patients. The niche slopes $\beta_{g,t,t'}$ also draw on within-patient variation in the niche covariate, so they carry a larger effective sample size, and one $S - 2$ cannot be right for both. df.method = "satterthwaite" (the default) gives each tested coefficient its own df,

$\nu_j \;=\; \frac{2\,\big(\widehat{\mathrm{Var}}\,\hat\theta_j\big)^2} {\mathrm{Var}\!\big(\widehat{\mathrm{Var}}\,\hat\theta_j\big)},$

with the denominator propagating the uncertainty of the variance parameters $(\varphi, \tau^2_{\text{sample}}, \tau^2_{\text{sample}\times\text{type}}, \tau^2_{\text{slope}})$ through the standard error, from the reduced-form REML information of the working model — a small matrix built from the $p \times p$ penalised inverse, never an $n \times n$ one. The ratio is invariant to $\varphi_g$ , so one df vector per bandwidth serves every gene. df.method = "between" uses the scalar $S - 2$ for every coefficient and is severely conservative for the slopes when patients are few.

fs <- fitSpiDE(spe, "condition", sigma = 30, random = "intercept",
               df.method = "satterthwaite", verbose = FALSE)
ff <- fits(fs)[[1]]
ct <- as.character(ff$covtype)
resp  <- ff$coefmap$covariate[ct == "ResponseCellType"]
niche <- ff$coefmap$covariate[ct == "ResponseNiche"]
c(ResponseCellType.median = round(median(ff$df[resp]), 1),
  ResponseNiche.median = round(median(ff$df[niche]), 1))
#> ResponseCellType.median    ResponseNiche.median 
#>                    11.2                   857.0

The flat-response df sits near $S - 2 = 8$ and the slope df is much larger. @df is therefore a named vector aligned to the columns of @t_stat; index the column you are testing rather than passing the vector to pt().

The null calibration of the corrected fit, against the fixed-effects one from earlier:

c(fixed = mean(fixed$p < 0.05), mixed = mean(mixed$p < 0.05))
#>     fixed     mixed 
#> 0.6233333 0.0000000

qq_points <- function(p, mode) {
  data.frame(mode = mode, expected = sort(-log10(ppoints(length(p)))),
             observed = sort(-log10(p)))
}
qq_df <- rbind(qq_points(fixed$p, "fixed"), qq_points(mixed$p, "mixed"))
ggplot(qq_df, aes(expected, observed, colour = mode)) +
  geom_abline(slope = 1, intercept = 0, linetype = 2, colour = spide_pal$ref) +
  geom_point(size = 1.5) +
  scale_colour_manual(values = unname(spide_pal$method[c("fixed", "intercept")]),
                      name = NULL) +
  labs(x = expression(-log[10]("expected p")), y = expression(-log[10]("observed p")))

Null response-effect p-values: the fixed-effects fit against the mixed fit.

From coefficients to results

A gene has one tested slope per (index type, niche type), and a gene’s niche columns are correlated because niche densities are spatially smooth. Three steps turn the per-coefficient statistics into the results() table.

Within a gene, across niches. The two-sided p-values of a gene’s slopes are combined by the Cauchy combination (combine = "cauchy", the default): $T = \sum_j w_j \tan\{(\tfrac12 - p_j)\pi\}$ , whose null distribution is Cauchy whatever the correlation between the $p_j$ , so no correlation matrix has to be estimated. It is computed for the gene as a whole and per index cell type. combine = "brown" is the alternative: Brown’s extension of Fisher’s method, which consumes the coefficient correlation matrix and combines one-sided p-values, separately for the up and down directions. The sidedness differs by necessity — the tangent transform diverges as a one-sided $p \to 1$ , so under Cauchy a gene up in one niche and down in another would cancel, whereas Brown’s $-2\log p$ is bounded and safe one-sided. Both populate the same slots, so everything downstream is combiner-agnostic.

Across bandwidths. The gene-level and per-index combined p-values from each bandwidth are combined by a second Cauchy combination, weighted by the gene’s relative log-likelihood across bandwidths, $w_\sigma \propto \exp(\ell_{g,\sigma} - \max_\sigma \ell_{g,\sigma})$ with small weights thresholded, so the bandwidth that fits a gene best carries the most weight.

The hierarchical FDR. A three-level nested Benjamini–Hochberg cascade: gene level; then index cell type within gene; then (index, niche) within those — the outer levels gated at half the nominal rate per direction and merged, the innermost at the nominal rate. What survives is a row keyed by (gene, index type, niche type) with the best bandwidth, its direction and its three FDR values.

Niche-only mode

With condition = NULL the design drops the condition terms and the bare niche main effects,

~ 0 + <covariates> + CellType + CellType:(niches)

and the CellType:niche slopes — nuisance in the conditioned model — become the tested quantity, cell-means coded. The main effects must go: a niche’s main effect equals the sum of its per-type interactions, so keeping it would alias one interaction per niche away and leave that cell type untestable. Everything downstream is shared. The random effects apply as before; note that the tested slope is now a within-sample contrast on a spatially autocorrelated covariate, which spiDE does not model, so this mode is mildly anti-conservative and borderline calls deserve corresponding caution.

Reading a fit

A fit with random != "none" carries these slots on each SpiDEFit (fits(res)[[1]], or $):

slot	what it holds	how to read it
`@tau2`	the variance components, named `SampleInt`, `SampleCellTypeInt` (default `re.celltype = TRUE`) and `SampleSlope` (`random = "slope"`)	on the log-mean scale; index by name, never by position
`@df`	the reference df: a named per-coefficient vector (`satterthwaite`) or a scalar $S-2$ (`between`)	index the column you test; small values mean few patients and a wide-tailed $t$ , which is expected
`@penalty`	the per-column ridge penalty used, $0$ on fixed columns and $1/\hat\tau^2_k$ on random ones	a very small penalty on a block means that source of between-patient variation converged large
`@polish`	the per-gene convergence diagnostics	see What the stage reports

f <- fits(fs)[[1]]
f$tau2
#>         SampleInt SampleCellTypeInt 
#>          4.638875          1.322273
colSums(f$polish[, c("restarted", "capped", "singular", "psi_bound")])
#> restarted    capped  singular psi_bound 
#>         0         0         0         0

Three caveats govern use.

Sample-level covariates (Age, Sex, Stage) are collinear with the per-sample random intercept and are rejected under random != "none"; the intercept subsumes them. Adjust for cell-level covariates only — and do adjust for per-cell library size, which the cell-type intercepts cannot carry.
Few samples. The variance components and the reference df are only as good as the between-sample information; a warning is issued when samples are very few, and the trade-offs by sample size are in the calibration vignette.
The mixed fit is not “fixed with extras”. Even when $\hat\tau^2 \to 0$ it does not reproduce random = "none": it also replaces the sub-design covariance and the NB dispersion scaling with the full penalised covariance and the Pearson dispersion, so its standard errors are the correctly scaled ones and "none"’s apparent extra power is partly spurious.

The patient-level question

What the nested intercept removes from the niche slope is a real association in the data. “Do patients whose tumour compartment is fibroblast-rich express this gene differently in tumour cells?” is a legitimate question with one observation per patient, and compositionTest() asks it there: per (sample, index type), a pseudobulk log-CPM profile of the index cells and the mean niche density around them; then, per (index, niche) pair, a limma moderated $t$ across samples for the pooled association ("niche") and for its difference between conditions ("condition:niche", the patient-level counterpart of $\beta$ ), with patient-level covariates allowed — exactly the ones fitSpiDE() rejects. On a cohort it typically reports that neighbour types’ own markers rise in the index cells’ pseudobulk with the neighbours’ density, uniformly across conditions: composition, plausibly with a segmentation-spillover component, and the reason the niche test keeps it out. Report the two results side by side, and never one as the other.

Choosing the estimator

fitSpiDE(random = "intercept") with the defaults is the estimator for discovery: it fits all cells jointly, so thin cell types borrow strength, and its hierarchical FDR is built for the full hypothesis space. random = "none" is kept for reproducing older results only. twoStageSpiDE() estimates each patient’s niche slopes first and contrasts them across patients, removing pseudo-replication by construction; it is the estimator for a pre-specified hypothesis on well-populated index types, and it is documented in vignette("spiDE-twostage"). The measured operating characteristics of every choice named in this vignette are collected in vignette("spiDE-calibration").

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] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] MASS_7.3-65                 ggplot2_4.0.3              
#>  [3] SpatialExperiment_1.22.0    SingleCellExperiment_1.34.0
#>  [5] SummarizedExperiment_1.42.0 Biobase_2.72.0             
#>  [7] GenomicRanges_1.64.0        Seqinfo_1.2.0              
#>  [9] IRanges_2.46.0              S4Vectors_0.50.2           
#> [11] BiocGenerics_0.58.1         generics_0.1.4             
#> [13] MatrixGenerics_1.24.0       matrixStats_1.5.0          
#> [15] spiDE_0.99.17               BiocStyle_2.40.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] deldir_2.0-4           rlang_1.3.0            magrittr_2.0.5        
#>  [4] otel_0.2.0             compiler_4.6.1         spatstat.geom_3.8-2   
#>  [7] systemfonts_1.3.2      callr_3.8.0            vctrs_0.7.3           
#> [10] pkgconfig_2.0.3        fastmap_1.2.0          magick_2.9.1          
#> [13] XVector_0.52.0         labeling_0.4.3         rmarkdown_2.32        
#> [16] ps_1.9.3               ragg_1.5.2             torch_0.17.0          
#> [19] bit_4.6.0              xfun_0.60              cachem_1.1.0          
#> [22] beachmat_2.28.0        jsonlite_2.0.0         goftest_1.2-3         
#> [25] DelayedArray_0.38.2    spatstat.utils_3.2-4   BiocParallel_1.46.0   
#> [28] irlba_2.3.7            parallel_4.6.1         R6_2.6.1              
#> [31] bslib_0.12.0           RColorBrewer_1.1-3     spatstat.data_3.1-9   
#> [34] limma_3.68.5           spatstat.univar_3.2-0  jquerylib_0.1.4       
#> [37] Rcpp_1.1.2             bookdown_0.48          knitr_1.52            
#> [40] tensor_1.5.1           Matrix_1.7-5           splines_4.6.1         
#> [43] tidyselect_1.2.1       abind_1.4-8            yaml_2.3.12           
#> [46] codetools_0.2-20       spatstat.random_3.5-1  spatstat.explore_3.8-2
#> [49] processx_3.9.0         lattice_0.22-9         tibble_3.3.1          
#> [52] withr_3.0.3            S7_0.2.2               evaluate_1.0.5        
#> [55] desc_1.4.3             SpaNorm_1.7.9          polyclip_1.10-7       
#> [58] pillar_1.11.1          BiocManager_1.30.27    scales_1.4.0          
#> [61] coro_1.1.0             glue_1.8.1             tools_4.6.1           
#> [64] ScaledMatrix_1.20.0    locfit_1.5-9.12        fs_2.1.0              
#> [67] grid_4.6.1             edgeR_4.10.4           nlme_3.1-169          
#> [70] BiocSingular_1.28.0    cli_3.6.6              rsvd_1.0.5            
#> [73] spatstat.sparse_3.2-0  textshaping_1.0.5      S4Arrays_1.12.0       
#> [76] dplyr_1.2.1            gtable_0.3.6           sass_0.4.10           
#> [79] digest_0.6.39          SparseArray_1.12.2     rjson_0.2.23          
#> [82] farver_2.1.2           htmltools_0.5.9        pkgdown_2.2.1         
#> [85] lifecycle_1.0.5        statmod_1.5.2          bit64_4.8.6

Dharmesh D. Bhuva