Statistical Methods

This page describes the statistical framework underlying SPLISOSM’s spatial variability (SV) and differential isoform usage (DU) tests. For a full derivation and theoretical justification, please refer to the Supplementary Notes of the SPLISOSM paper [SQS+26], and of [SFW+26].

Overview

SPLISOSM treats every gene as a multivariate object: a gene with \(p\) isoforms at \(n\) spatial locations is represented by an \(n \times p\) matrix. The objective of the SV test is to determine whether the multivariate \(p\)-dimensional vector varies across the \(n\) spatial locations, while the DU test is to determine whether variation in this vector is associated with other spatial covariates (e.g., spatial domains, RBP expression), potentially after conditioning on spatial autocorrelation.

The main statistical tool for both tests is a kernel-based measure of association called the Hilbert-Schmidt Independence Criterion (HSIC) [GBSScholkopf05]. Intuitively,

• The spatial kernel \(K \in \mathbb{R}^{n \times n}\) encodes the spatial structure of the tissue.

• The response kernel \(L \in \mathbb{R}^{n \times n}\) encodes similarity between isoform profiles of different spots/cells.

• HSIC measures whether spots that are spatially close (large \(K_{ii'}\)) also tend to have similar isoform profiles (large \(L_{ii'}\)).

Spatial Kernel: CAR Model

SPLISOSM uses a Conditional Autoregressive (CAR) model [SRF+23] to define spatial covariance. Given a k-mutual-nearest-neighbor adjacency matrix \(W\) built from spot coordinates, the CAR precision matrix is

\[M = I - \rho D^{-1/2} W D^{-1/2},\]

where \(D = \mathrm{diag}(\sum_{i'} W_{ii'})\) is the degree matrix and \(\rho \in (0, 1)\) is the spatial autocorrelation coefficient. The spatial covariance (kernel) matrix is \(K = M^{-1}\), standardised to unit marginal variance.

We choose the CAR kernel for the following properties:

• Sparse precision: only the k-NN graph is stored explicitly; \(K\) is never formed for large datasets (implicit LU-solve mode when \(n > 5000\)).

• Irregular geometries: the k-NN graph is translation-invariant and naturally handles non-grid layouts (Visium spots, single cells, segmented tissue).

• Graph Fourier basis: eigenvectors of \(K\) form a graph Fourier basis ordered by spatial frequency.

• Polynomial spectrum decay: the eigenvalues of \(K\) decay polynomially, which yields higher power than exponential-decay kernels for mix-frequency patterns.

For general theory and the equivalence of spatial variability testing methods, see our recent work on the consistent and scalable detection of spatial patterns [SFW+26].

Spatial Variability (SV) Tests

Three SV tests are available, differing only in how the response matrix \(Y \in \mathbb{R}^{n \times p}\) is constructed.

Test types

Let \(c_{ij}\) denote the raw count at spot \(i \in \{1, \ldots, n\}\) of isoform \(j \in \{1, \ldots, p\}\), and \(r_{ij} = c_{ij} / \sum_{j'} c_{ij'}\) the corresponding isoform usage ratio.

Test	Response \(Y\)	Null hypothesis	Typical use case
HSIC-IR	Centred isoform usage ratios \(r_{ij}\) (optionally log-ratio transformed)	Isoform usage is spatially uniform	Identify SVP genes (spatially variable RNA processing)
HSIC-GC	Centred total gene count \(\sum_j c_{ij}\) (single column)	Gene expression is spatially uniform	Identify SVE genes; drop-in for SPARK-X
HSIC-IC	Centred raw isoform counts \(c_{ij}\)	Isoform counts are spatially uniform	Reflects joint changes in expression and processing

HSIC-IR is the recommended test for discovering spatially variable RNA processing (SVP genes). HSIC-GC is equivalent to a gene-level spatial variability test and can serve as a drop-in replacement for SPARK-X with improved statistical power. HSIC-IC tests the joint null; significance can arise from either differential expression or differential processing. HSIC-IC does not test each isoform individually and instead yields a single test statistic per gene. In practice, the results of HSIC-IC and HSIC-GC are often similar, as changes in overall gene expression are the main driver of isoform expression changes for many genes.

Ratio transformations

For HSIC-IR, the usage ratios are optionally transformed before computing the HSIC statistic. The ratio_transformation argument controls this:

• 'none' (default): raw proportions \(r_{ij}\), mean-centred per isoform.

• 'clr': centred log-ratio \(\log(r_{ij}) - \frac{1}{p}\sum_{j'} \log(r_{ij'})\).

• 'ilr': isometric log-ratio (orthonormal Helmert contrast in log-simplex).

• 'alr': additive log-ratio relative to the last isoform.

• 'radial': radial transformation \(r_{ij} / \|r_i\|\) [PYPA22], where \(r_i = (r_{i1}, \ldots, r_{ip})\) is the isoform-ratio vector at spot \(i\). Empirically, it is not calibrated and thus not recommended.

Due to excessive sparsity, log-ratio-based transformations require pseudo-counts to avoid zero ratios, which may lead to unwanted artefacts. Our empirical results suggest that the untransformed ratios ('none') are well-calibrated, robust and often more powerful for SVP detection.

Test statistic

SPLISOSM uses a quadratic form V-statistic, which computes HSIC with a linear kernel on \(Y\):

\[\widehat{\mathrm{HSIC}} = \frac{1}{(n-1)^2} \mathrm{tr}(Y^\top K Y),\]

where \(K\) is the (double-centred, optionally standardised) CAR kernel. This equals the sum of squared spatial autocorrelations of each linear combination of isoform profiles. Since we have

\[\mathrm{tr}(Y^\top K Y) = \mathrm{tr}(K \cdot Y Y^\top),\]

the test statistic can also be viewed as a weighted sum of spatial autocorrelations of the principal components of \(Y\).

Null distribution

To compute p-values, we need to compute the distribution of \(\widehat{\mathrm{HSIC}}\) under the null hypothesis of no spatial variability. Three methods are available via the null_method argument to test_spatial_variability().

1. Liu’s approximation of chi-square mixture (default): null_method='liu'

Let \(Q = \mathrm{tr}(Y^\top K Y) = (n-1)^2 \widehat{\mathrm{HSIC}}\) denote the unnormalised HSIC V-statistic, write \(m=n-1\) for the dimension of the centred spot space, with \(\lambda_1^K \geq \cdots \geq \lambda_n^K\) the eigenvalues of \(K\) and \(\lambda_1^Y \geq \cdots \geq \lambda_p^Y\) those of \(Y^\top Y\). Under the null, \(Q\) asymptotically follows a weighted sum of independent \(\chi^2_{1}\) variables,

\[Q \;\overset{d}{\approx}\; \frac{1}{n-1} \sum_{i=1}^{n} \sum_{j=1}^{p} \lambda_i^K \, \lambda_j^Y \; Z_{ij}, \qquad Z_{ij} \overset{\text{iid}}{\sim} \chi^2_{1},\]

where the double sum runs over all pairs of spatial and response eigenvalues. SPLISOSM evaluates Liu et al. [LTZ09] from the first four cumulants \(c_r = \mathrm{tr}(K^r)\mathrm{tr}((Y^\top Y)^r)/(n-1)^r\), so it does not need to materialize the full pairwise eigenvalue product. See splisosm.utils.hsic.liu_sf_from_cumulants() for implementation details.

Note

• The previous name null_method='eig' is retained as a deprecated alias and will emit a DeprecationWarning.

• Full eigen-decomposition of \(K\) is \(O(n^3)\) and is not feasible for large datasets. For implicit CAR kernels with no realised dense covariance, SPLISOSM estimates \(\mathrm{tr}(K^r)\) with Hutchinson Rademacher probes by default. Use null_configs={"n_probes": m} to control that budget.

• Since v1.2.0, large SplisosmNP SV tests use full-rank cumulants by default rather than the v1.1.x low-rank spatial approximation. This preserves sensitivity to both global and local spatial patterns. If the analysis should intentionally emphasize global smooth patterns, prefer increasing the CAR smoothness parameter (for example rho=0.999) instead of truncating the spatial rank. See the Visium HD SV hyperparameter tutorial for a practical comparison against an FFT regular-grid reference.

• When nan_filling='mean' (default), the spatial cumulants are cached once and reused for all subsequent genes.

• When nan_filling='none' for HSIC-IR, SPLISOSM drops zero-coverage spots per gene and applies a masked implicit spatial kernel. This avoids dense parent-kernel realization, but per-gene masked cumulants must still be estimated and the path is slower than mean filling.

2. Welch-Satterthwaite scaled chi-squared approximation: null_method='welch'

Alternatively, we may use only the first two moments of the null distribution, which requires \(\mathrm{tr}(K)\) and \(\mathrm{tr}(K^2)\) but not higher cumulants or a full eigendecomposition. Taking moments of the \(\chi^2_1\) mixture above gives

\[\begin{split}\mu_0 &= \mathbb{E}[Q] = \frac{1}{m}\,\mathrm{tr}(K)\,\mathrm{tr}(Y^\top Y), \\[4pt] \sigma_0^2 &= \mathrm{Var}(Q) = \frac{2}{m^2}\,\mathrm{tr}(K^2)\,\mathrm{tr}\!\bigl((Y^\top Y)^2\bigr).\end{split}\]

With all positive eigenvalues, the chi-squared mixture null can be approximated by the Welch-Satterthwaite method, using one scaled chi-squared variable \(g \, \chi^2_h\) with scale parameter \(g\) and degrees of freedom \(h\). The parameters are chosen to match the first two moments of the null \((\mu_0, \sigma_0^2)\).

\[g = \frac{\sigma_0^2}{2\mu_0}, \qquad h = \frac{2\mu_0^2}{\sigma_0^2},\]

and the p-value is \(\mathbb{P}\!\left(\chi^2_h \geq Q/g\right)\).

Note

• The retired null_method='clt' and null_method='trace' names are retained as deprecated aliases and automatically use this Welch approximation.

• Using only the first two cumulants \(\mathrm{tr}(K)\) and \(\mathrm{tr}(K^2)\), Welch is typically less accurate in the right tail (small p-values) than the Liu approximation, which uses four cumulants. In practice the difference is often small.

3. Batched permutation test: null_method='perm'

Generates a null distribution by repeatedly permuting the rows of \(Y\) (breaking the spatial structure) and recomputing \(\mathrm{tr}(Y_\pi^\top K Y_\pi)\) for each permutation \(\pi\). P-values use the finite-permutation correction \((1 + \#\{Q_\pi \ge Q_{\mathrm{obs}}\}) / (B + 1)\).

To avoid \(O(n^2)\) memory, SPLISOSM batches \(B\) permutations into a single \((n, Bp)\) matrix and calls xtKx() once per batch. Per-permutation traces are recovered as diagonal blocks of the \((Bp, Bp)\) result matrix without materialising the full kernel.

Configuration via null_configs:

• n_perms_per_gene (default 1000): total permutations.

• perm_batch_size (default 50): permutations per xtKx() call.

FFT acceleration for regular grids

SplisosmFFT exploits the translation-invariance of regular grids (Visium HD, Xenium binned data) to accelerate the SV test via the 2-D Fast Fourier Transform. Specifically, it reduces the kernel eigen-decomposition required for the Liu’s method from \(O(n^3)\) to \(O(n \log n)\), and the quadratic form computation from \(O(n^2 p)\) to \(O(n p \log n)\).

On a regular \(H \times W\) grid with \(n = HW\) spots, uniform degree \(d\) and periodic boundaries, the CAR precision \(M = I - (\rho/d)\,W\) is block-circulant. Its eigenvalues (and hence those of \(K = M^{-1}\)) are the 2-D DFT of its first row, computable in \(O(n \log n)\) time:

\[\lambda_{(h,w)}^K = \frac{1}{1 - \rho \, \hat{W}_{hw} / d}, \qquad h = 0, \ldots, H-1,\; w = 0, \ldots, W-1,\]

where \(\hat{W}_{hw} = (\mathcal{F}\,w)_{hw}\) is the unnormalised 2-D DFT of the first row of the adjacency matrix \(W\).

The quadratic form \(\mathrm{tr}(Y^\top K Y)\) then reduces to a pointwise product in Fourier space. Reshape each isoform image \(y_j \in \mathbb{R}^n\) to \(H \times W\) and let \(\hat{Y}_{hw} = \bigl((\mathcal{F} y_1)_{hw},\, \ldots,\, (\mathcal{F} y_p)_{hw}\bigr) \in \mathbb{C}^{p}\) denote the vector of 2-D DFT coefficients at frequency \((h, w)\) stacked across all \(p\) isoforms. Then

\[\mathrm{tr}(Y^\top K Y) \;=\; \frac{1}{n} \sum_{h=0}^{H-1} \sum_{w=0}^{W-1} \lambda_{(h,w)}^K \; \|\hat{Y}_{hw}\|_2^2,\]

where the leading \(1/n\) comes from the unnormalised DFT convention (scipy.fft.fft2) that satisfies \(\mathcal{F}^{-1} = (1/n)\,\mathcal{F}^{*}\). This reduces the quadratic-form computation from \(O(n^2 p)\) to \(O(n p \log n)\).

Furthermore, the spatial eigenvalues \(\{\lambda_{(h,w)}^K\}\) are shared across all genes (computed once), so the cost of the eigenvalue null is also \(O(1)\) per gene in terms of kernel operations.

Note

For irregular coordinates, use SplisosmNP for SV tests. The FINUFFT backend described below is currently used for GP residualization in conditional DU tests, not for the SV spatial kernel.

Differential Isoform Usage (DU) Tests

DU tests ask whether isoform usage is associated with a covariate (e.g., spatial domain label, RBP expression), potentially conditioned on spatial autocorrelation.

Unconditional test (method='hsic')

The unconditional test directly applies HSIC between the centred covariate vector \(Z \in \mathbb{R}^n\) and the isoform ratio matrix \(Y\):

\[\widehat{\mathrm{HSIC}}_{\text{uncond}} = \frac{1}{(n-1)^2} \mathrm{tr}(Y^\top K_Z Y),\]

where \(K_Z = ZZ^\top\) is a rank-1 linear kernel on the covariate, so that \(\mathrm{tr}(Y^\top K_Z Y) = \|Z^\top Y\|^2\). Up to a scalar, the statistic is equivalent to the multivariate RV coefficient between \(Z\) and \(Y\). When \(Z\) is binary, this is also equivalent to a two-step procedure where two-sample T-test is first applied to each isoform separately, and the resulting p-values are combined (e.g., via Fisher’s method method='t-fisher').

Limitation: If both isoform usage and the covariate are spatially autocorrelated (which is common), the unconditional test can be anti-conservative, meaning that it will report false positive associations that merely reflect shared spatial structure.

Conditional test (method='hsic-gp')

Given spatial coordinates \(X\), the conditional test assesses the association between \(Z | X\) and \(Y | X\). This is a difficult problem because in most spatial data, each spatial location is observed only once (i.e., the above conditionals are point mass). To condition on spatial autocorrelation, we adopt a residualisation approach [ZPJScholkopf12].

Specifically, SPLISOSM first residualises the covariate against a Gaussian Process (GP) spatial model, then tests the independence of the residuals against isoform usage:

1. Fit a GP to the covariate \(z\) using spatial coordinates \(x\):

   \[z = f(x) + \varepsilon, \quad f \sim \mathcal{GP}(0, k_\theta),\]

   where \(k_\theta\) is a Constant x RBF + WhiteNoise kernel. SPLISOSM provides dense sklearn, GPyTorch, FFT, and FINUFFT-backed NUFFT GP backends for different data geometries and scales.

2. Compute covariate residuals \(\tilde{Z} = Z - \hat{f}(X)\), capturing the part of covariate variation not explained by spatial position.

3. Test \(\widehat{\mathrm{HSIC}}(\tilde{Z},\, Y)\) using the rank-1 linear covariate kernel \(K_{\tilde{Z}} = \tilde{Z}\tilde{Z}^\top\) and a similar linear kernel \(K_Y = Y Y^\top\) for the response. Since both kernels are low-rank, the Liu null can be computed from the nonzero covariate and response eigenvalues.

Note

Theoretically, a fully equivalent test of conditional independence \(Z \perp\!\!\!\perp Y \mid X\) would residualise both \(Z \mid X\) and \(Y \mid X\) against the spatial GP. However, \(\tilde{Z} \perp\!\!\!\perp Y\) implies \(Z \perp\!\!\!\perp Y \mid X\) (but not vice versa), so any significant dependency found by residualising only the covariate is also a true positive for the conditional association.

In practice, we found that the covariate-only residualisation delivers massive computational savings and remains well-calibrated. For flexibility, we provide optional control via the residualize argument:

• residualize='cov_only' (default): only \(Z\) the spatial covariate is residualised; \(Y\) the isoform ratio matrix is used as-is.

• residualize='both': both \(Z\) and \(Y\) are GP-residualised before testing.

Note

The GP fitting step is the dominant computational cost of the DU test. Three backend families are supported:

• gpr_backend='sklearn' (default): dense sklearn GaussianProcessRegressor with optional subset-of-data hyperparameter fitting.

• gpr_backend='gpytorch': exact or FITC sparse GP with n_inducing inducing points and optional GPU support.

• gpr_backend='nufft' / 'finufft': FINUFFT-backed implicit RBF grid-kernel for irregular 2-D coordinates, recommended for large-scale spatial data.

For sklearn and gpytorch, pass gpr_configs={"covariate": {"n_inducing": 1000}} to control the subset or inducing-point budget. For the NUFFT backend, max_auto_modes caps the automatically inferred full effective grid, and lml_approx_rank controls the irregular-grid likelihood approximation used during hyperparameter fitting. These options are passed through gpr_configs; see splisosm.SplisosmNP.test_differential_usage() for the user-facing configuration table and GP Regression Backends for backend class details.

NUFFT backend for irregular-coordinate GP residualization

The gpr_backend='nufft' / 'finufft' path targets irregular 2-D coordinates where dense GP matrices are too expensive. Coordinates are affinely mapped into a periodic box \([-\pi,\pi)^2\), and the stationary RBF covariance is represented on a Fourier mode set \(\Omega\) as

\[k_\theta(x_i, x_j) \approx \sigma_f^2 \sum_{\omega \in \Omega} a_\theta(\omega) \exp\{i\omega^\top(t_i - t_j)\} \;+\; \sigma_\varepsilon^2 \mathbf{1}_{i=j},\]

where \(t_i\) are the mapped coordinates and \(a_\theta(\omega)\) are non-negative spectral weights determined by the RBF length scale and Fourier grid spacing. With n_modes=None, SPLISOSM uses the full effective Fourier grid inferred from the point count and coordinate aspect ratio; FINUFFT’s internal oversampling is controlled separately via nufft_opts.

With this representation, we can compute the matrix-vector product in \(O(n \log n)\) time by

\[K_s v \;\approx\; F_X^* \left(a_\theta \odot F_X v\right),\]

where \(F_X v = \sum_i v_i \exp(-i\omega^\top t_i)\) is a type-1 NUFFT and \(F_X^*\) is the corresponding type-2 NUFFT. GP residualization solves

\[(K_s + \sigma_\varepsilon^2 I)\alpha = z, \qquad \tilde z = \sigma_\varepsilon^2 \alpha,\]

using conjugate gradients, because \(K_s\) is only available through NUFFT matvecs.

On compatible regular grids, the same spectral representation reduces to the FFT GP path because the Fourier modes are exact eigenvectors. On irregular coordinates the nonuniform Fourier features are not orthogonal eigenvectors, so hyperparameter fitting approximates the log marginal likelihood with leading NUFFT eigensummaries plus trace and trace-square tail corrections. The lml_approx_rank parameter controls this eigensummary rank; it affects hyperparameter fitting accuracy and memory, not the Fourier grid used for each matvec. Memory scales as \(O(nr)\) for rank r = lml_approx_rank (about 512 MB for 1M spots and r=64 in float64), while the NUFFT matvecs avoid forming the \(n \times n\) dense GP matrix.

FFT-accelerated conditional DU test

SplisosmFFT accelerates the conditional DU test (method='hsic-gp') by replacing the dense GP with an FFT-based Gaussian process.

On a regular grid with periodic boundaries the GP covariance kernel \(k_\theta(x, x')\) is stationary, so the Gram matrix \(K_\theta\) is block-circulant with eigenvalues \(\hat{k}_{\theta,(h,w)} = (\mathcal{F}\,k_\theta)_{hw}\) (2-D DFT). All operations required for GP fitting and prediction — matrix-vector products, log-determinants, and gradient computations — can be performed in the spectral domain:

\[K_\theta \cdot v \;=\; \tfrac{1}{n}\, \mathcal{F}^{*}\!\bigl(\hat{k}_\theta \odot \mathcal{F}(v)\bigr), \qquad \log |K_\theta| \;=\; \sum_{h,w} \log \hat{k}_{\theta,(h,w)},\]

where the \(1/n\) in the matrix-vector product reflects the unnormalised DFT convention \(\mathcal{F}^{-1} = (1/n)\,\mathcal{F}^{*}\), and the log-determinant identity follows from \(\hat{k}_{\theta,(h,w)} > 0\) for a positive-definite stationary kernel. This reduces the per-step GP cost from \(O(n^3)\) (dense Cholesky) or \(O(nM^2)\) (inducing-point) to \(O(n \log n)\) per L-BFGS iteration, with no approximation error for stationary kernels on the grid.

The GP hyperparameters (signal variance, length scale, noise variance) are optimised by maximising the marginal log-likelihood using L-BFGS with gradients computed in Fourier space. Covariate residuals \(\tilde{Z} = Z - \hat{f}(X)\) are then obtained via the FFT-based posterior mean.

After residualisation, the HSIC test itself proceeds identically to the SplisosmNP path: linear-kernel HSIC via \(\mathrm{tr}(Y^\top K_{\tilde{Z}} Y)=\|\tilde{Z}^\top Y\|^2\), with p-values from Liu’s chi-squared mixture null using the eigenvalues of \(K_{\tilde{Z}} = \tilde{Z}\tilde{Z}^\top\) (rank-1, giving a single nonzero eigenvalue \(\|\tilde{Z}\|^2\)) and \(K_Y = Y Y^\top\).

Parametric test: SplisosmGLMM

SplisosmGLMM provides a parametric alternative based on a multinomial GLMM with a Gaussian random field (GRF) random effect. Using a reference-category multinomial-logit link with \(q = p - 1\) free dimensions per spot, the model is

\[\begin{split}Y_i \mid \eta_i &\sim \mathrm{Multinomial}(N_i,\, \mathrm{softmax}([\eta_i, 0])), \quad i = 1, \ldots, n, \\ \eta_i &= X_i \beta + b + u_i \;\in\; \mathbb{R}^{q}, \\ \mathrm{vec}(U) &\sim \mathcal{N}\!\bigl(0,\; \sigma^2 \, \Sigma \otimes I_q \bigr), \quad \Sigma \;=\; \theta\, K \;+\; (1-\theta)\, I_n,\end{split}\]

where \(X_i\) is the row of covariates for spot \(i\), \(\beta \in \mathbb{R}^{d \times q}\) are fixed effects, \(b \in \mathbb{R}^{q}\) is an intercept, and \(U = (u_1, \ldots, u_n)^\top \in \mathbb{R}^{n \times q}\) stacks the per-spot random effects. The covariance mixes the CAR spatial kernel \(K\) with an i.i.d. component via \(\theta \in [0, 1]\), and assumes independence across the \(q\) logit dimensions (\(\otimes I_q\)).

Note

The SV test (H₀: \(\theta = 0\)) is implemented as a likelihood ratio test (LRT) in test_spatial_variability(). However, it is not well-calibrated due to technical challenges in model fitting. The equivalent score test version also takes a quadratic form similar to the HSIC test statistic but with spot-specific adjustments. See [SFW+26] for detailed analysis.

For DU testing, we use a score test comparing coefficient gradients at the null model (no fixed-effect covariates, \(\beta = 0\)), which avoids fitting the full model for each covariate. However, it still requires estimating nuisance parameters (intercept \(b\), total variance \(\sigma^2\), and spatial variance proportion \(\theta\)). To compute the maximum likelihood estimates, we approximate the marginal likelihood via Laplace’s method at the mode of the random effects. See splisosm.glmm.MultinomGLMM for implementation details.

Compared to SplisosmNP, the GLMM approach:

• Makes stronger distributional assumptions (multinomial counts), which may or may not be well-satisfied in practice. See [SFW+26] for a theoretical analysis of the effect of link function on test power.

• Can be more interpretable (when fitting with effect sizes \(\hat\beta\) and standard errors).

• Is more computationally intensive (per-gene numerical optimisation).

Multiple testing correction

All SV and DU tests apply Benjamini-Hochberg (BH) FDR correction [BH95] across genes (for SV) or gene-covariate pairs (for DU, each covariate factor is tested and adjusted separately). The adjusted p-values are stored in the pvalue_adj column of the results DataFrame returned by get_formatted_test_results().

Summary

Test	Class	Response \(Y\)	Conditioning	Null method
HSIC-IR (SVP)	SplisosmNP / FFT	Isoform ratios	None	liu [1] / welch / perm
HSIC-GC (SVE)	SplisosmNP / FFT	Gene counts	None	liu [1] / welch / perm
HSIC-IC	SplisosmNP / FFT	Isoform counts	None	liu [1] / welch / perm
DU (unconditional)	SplisosmNP / FFT	Isoform ratios	None	liu
DU (GP-conditional)	SplisosmNP / FFT	Isoform ratios	GP spatial residualisation	liu
DU (GLM score)	SplisosmGLMM	Isoform counts	None	Chi-squared (score)
DU (GLMM score)	SplisosmGLMM	Isoform counts	GRF random effect	Chi-squared (score)

[1] (1,2,3)

SplisosmFFT does not take a null_method parameter; 'liu' is the default and only option for FFT-accelerated tests, since the full kernel spectrum is efficiently computed via FFT.