nb_profile_psi – SpaNorm

nb_profile_psi(
    Y,
    W,
    alpha,
    psi,
    psi_range=(0.001, 1000.0),
    block_size=None,
    n_jobs=1,
    offset=None,
)

Profile-ML dispersion at fixed coefficients, blocked over genes (R nbProfilePsi).

For each gene, the NB dispersion that maximises the likelihood at the mean exp(W alpha + offset) with the coefficients held fixed: one profile search per gene and no Newton, by the same bisection on the score that the batched engine uses. A gene whose optimum sits on a bound of psi_range (an under-dispersed or near-empty gene), or whose search is not finite, keeps its input dispersion. A caller that re-polishes the mean at a held dispersion (polish_nb(warm=True)) uses this to put the dispersion back at the reported mean.

Parameters

Name	Type	Description	Default
Y		`genes x cells` counts (dense, sparse or tensor).	required
W		`cells x p` design or an `NBBlockDesign`.	required
alpha		`genes x p` coefficients.	required
psi		per-gene dispersions (length `nrow(Y)`, or one value), kept for a gene whose optimum is on a bound.	required
psi_range		search interval for the dispersion.	`(0.001, 1000.0)`
block_size	int \| None	gene blocking and dispatch, as in :func:`polish_nb`.	`None`
n_jobs	int \| None	gene blocking and dispatch, as in :func:`polish_nb`.	`None`
offset		`None`, a per-cell log-scale offset or a genes x cells matrix. The dispersion is profiled at the mean it gives, so it must be the offset the coefficients were fitted with.	`None`

Returns

Name	Type	Description
	numpy.ndarray	one dispersion per gene.

Examples

import numpy as np
import spanormpy as sn
rng = np.random.default_rng(1)
Y = rng.poisson(5, size=(20, 50))  # counts, genes x cells
W = np.column_stack([np.ones(50), np.linspace(-1, 1, 50)])  # intercept + one covariate
fit = sn.fit_nb(Y, W, verbose=False, backend="cpu")
sn.nb_profile_psi(Y, W, fit["alpha"], fit["psi"])[:3]  # profile-ML dispersions

array([0.00372586, 0.00743151, 0.00184694])