Quasi-likelihood F-tests

Quasi-likelihood F-tests#

fit_quasi_likelihood() and quasi_likelihood_f_test() are edgeR’s glmQLFit and glmQLFTest as edgeR 4 runs them by default: the negative binomial dispersion is the Cox-Reid common dispersion of the most highly expressed genes (common_dispersion()), each gene’s quasi-dispersion comes from deviances adjusted for small counts, and limma’s prior for unequal degrees of freedom moderates them, towards a trend in the average log CPM (average_log_cpm()). The tests give edgeR’s columns logFC, logCPM, F, PValue and FDR for coefficients or named contrasts of a design.

from cytepool.design import Design, make_contrasts
from cytepool.normalization import tmm_factors
from cytepool.stats.quasi_likelihood import (
    fit_quasi_likelihood,
    quasi_likelihood_f_test,
)

design = Design(samples, "~ 0 + group")
fit = fit_quasi_likelihood(counts, design, norm_factors=tmm_factors(counts))
contrasts = make_contrasts({"b_vs_a": "groupB - groupA"}, design)
test = quasi_likelihood_f_test(fit, contrasts=contrasts)
test.log_fold_change, test.p_value, test.adjusted_p_value

edgeR’s quasi-likelihood F-tests: glmQLFit and glmQLFTest as edgeR 4 runs them.

The quasi-likelihood model of Lund, Nettleton, McCarthy and Smyth (2012) gives each gene g a negative binomial model for its counts with a second, gene-wise dispersion:

var(y_gi) = sigma_g^2 (mu_gi + phi mu_gi^2),

with phi the negative binomial (NB) dispersion, shared by the genes, and sigma_g^2 the quasi-dispersion. The quasi-dispersions are estimated from the deviances of the NB fits, moderated by empirical Bayes as limma moderates residual variances, and genes are tested with F-tests whose denominator is the posterior quasi-dispersion. fit_quasi_likelihood() and quasi_likelihood_f_test() follow edgeR 4.4’s glmQLFit and glmQLFTest with their defaults (legacy = FALSE), the method of Chen, Chen, Lun, Baldoni and Smyth (2025), as documented and, where the documentation is silent, as found by running edgeR 4.4.0 as a black box (“found by probing”).

The steps:

  1. Abundance. average_log_cpm() is edgeR’s aveLogCPM: the log2 of the mean count per million of a gene, fitted as one group of samples by a NB model after a prior count of 2, scaled by library size, is added to every count (edgeR’s addPriorCount). glmQLFit computes it once, at the given NB dispersion, or at 0.05 when it estimates the dispersion (found by probing); it ranks the genes and is the covariate of both trends below. The tests report it at the NB dispersion.

  2. NB dispersion. Unless given, it is the Cox-Reid common dispersion (common_dispersion(), edgeR’s estimateGLMCommonDisp) of the most highly expressed genes: the ceiling(p G) of the G genes with the largest average log CPM, with p = chooseLowessSpan(G sqrt(n - k), small.n = 20, min.span = 0.02) for n samples and k coefficients, as edgeR’s help page says (the count of genes is a ceiling, found by probing).

  3. Fits. Every gene gets the NB fit of fit_negative_binomial_model() at the working dispersion phi / A, where A is the average quasi-dispersion below.

  4. Bias-adjusted deviances. The unit deviance d_i of an observation is approximated by a scaled chi-squared variable s_i chi^2(k_i) whose mean and variance are those of the unit deviance (Chen et al. 2025): s_i = V_i / (2 E_i) and k_i = 2 E_i^2 / V_i, with E_i and V_i from negative_binomial_unit_deviance_moments(). The adjusted unit deviance is d_i / s_i and its degrees of freedom are k_i (1 - h_i), with h_i the leverage of the observation. An observation of leverage 1, which the fit reproduces, has a residual deviance and degrees of freedom of zero, as in exact arithmetic, not the rounding of its fitted value. A gene’s adjusted residual deviance and degrees of freedom are their sums over its samples, and its quasi-dispersion is their ratio, on degrees of freedom that are fractional and differ between genes. The moments are those of the divided counts (Baldoni et al. 2024): with quasi-dispersion A, var(y) = A (mu + (phi / A) mu^2), and y / A follows a NB distribution with mean mu / A and dispersion phi. E_i and V_i are taken at the mean mu_i / A and the dispersion phi, and d_i at the working dispersion (found by probing).

  5. Average quasi-dispersion (found by probing). A starts at 1 and is updated twice: each time the genes are fitted at phi / A and their quasi-dispersions formed as in step 4; the fourth roots of the quasi-dispersions of the genes with residual degrees of freedom are smoothed against the abundance by R’s lowess with a span of 0.5 (three robustness iterations), and A becomes the fourth power of the 0.9 quantile of the smoothed values, or 1 if that is smaller. The fits of step 3 use the final A, which is edgeR’s average.ql.dispersion: its help page says that it “scales the NB dispersion”.

  6. Moderation. The quasi-dispersions are squeezed by limma’s squeezeVar with its estimator for unequal degrees of freedom (fitFDistUnequalDF1, legacy = FALSE whatever the degrees of freedom; found by probing), through squeeze_variances() with prior_estimator="likelihood": towards a trend in the abundance by default, and robustly when asked. winsor.tail.p plays no part, as in edgeR (found by probing).

  7. Tests. quasi_likelihood_f_test() fits the null model at the working dispersion, without the tested coefficients or with the tested contrasts set to zero, and divides the likelihood ratio statistic by its degrees of freedom r (the rank of the contrasts) and by the posterior quasi-dispersion. The F-statistic is referred to the F distribution on r and min(d0 + d_g, sum of n - k over the genes) degrees of freedom, with d0 the prior and d_g the adjusted residual degrees of freedom (found by probing). Log fold changes are the coefficients shrunk with a prior count of 0.125 (Phipson 2013), divided by log 2, and adjusted p-values are Benjamini and Hochberg’s (1995), as edgeR’s topTags gives them.

Prior weights (found by probing). edgeR uses them in the fits as prior weights, and in the Cox-Reid likelihood as the number of independent counts each count averages (edgeR’s adjustedProfileLik). In the adjusted deviances it treats an observation of weight w as w technical replicates: its moments are taken at the mean w mu / A, its unit deviance at the dispersion w phi / A, and its leverage is that of the weight mu / (1 + w phi mu / A) rather than of the weight of the fit; its adjusted unit deviance is per unit of weight, so that the gene’s adjusted deviance is the sum of the unit values times the weights, and its degrees of freedom the plain sum. This module follows edgeR in each. The weights are therefore not relative: multiplying all of them by a constant changes the quasi-dispersions and the F-statistics, in both.

Deliberate differences from edgeR 4.4.0:

  • The mean and variance of the unit deviance are exact, to about 1e-14. edgeR approximates them (found by probing): the scale implied by its adjusted unit deviances is off by up to about 1.2e-3 at dispersions up to 0.2 and 1.3e-2 at a dispersion of 1, and the degrees of freedom by up to about 2e-4 and 4e-3. The adjusted deviances, quasi-dispersions and F-statistics inherit differences of up to about 3e-3.

  • The Cox-Reid dispersion is the maximum of the adjusted profile likelihood, located to about 1e-7 relative; edgeR’s search stops within its tolerance of 1e-5. Where edgeR approximates the deviance, below a dispersion of 1e-4, its likelihood is noisy (found by probing), so that for counts close to Poisson the two estimates differ by more than that, both being near zero.

  • The fits converge to the maximum, where edgeR’s Levenberg fits stop short of it by up to about 3e-4 in the coefficients (see cytepool.stats.negative_binomial_models), so that the log fold changes of designs that are not one-way layouts differ by up to about that. Fitted values that tend to zero, as in a group without counts, are zero in both for one-way layouts; otherwise they stop near 1e-8 in edgeR’s fits and far lower here, and both give such observations the leverages of the limit. Either way they carry almost no degrees of freedom. A gene with no others, such as one counted in only one condition of a paired design, has the ratio of their almost-zero deviances and degrees of freedom as its quasi-dispersion, which therefore differs (0.053 in edgeR and 0.021 here, for such a gene in eight paired samples). It enters the average quasi-dispersion in both.

  • The robust prior tests only the genes with residual degrees of freedom for outliers, as fit_f_distribution_likelihood() does; limma counts a gene without them, such as one without counts, as a zero variance, which then lowers the outlier probabilities of the others.

  • With offsets per observation, edgeR’s aveLogCPM gives the first gene a value that does not depend on its offsets (found by probing: a copy of the gene in another row gets the value computed here), which then enters the ranking, the trends and the log CPM of that gene.

Not implemented yet: edgeR’s legacy = TRUE methods (the deviances of glmQLFit in edgeR 3, with the degrees of freedom adjusted for groups without counts, and the Poisson bound of glmQLFTest). This module follows edgeR’s method for a matrix of counts. edgeR’s method for a DGEList agrees with it, except when the DGEList holds a trended dispersion from estimateDisp: it then takes one dispersion from that trend in place of the Cox-Reid estimate (found by probing).

References

Lund SP, Nettleton D, McCarthy DJ, Smyth GK (2012). Detecting differential expression in RNA-sequence data using quasi-likelihood with shrunken dispersion estimates. Statistical Applications in Genetics and Molecular Biology 11(5), article 8. doi:10.1515/1544-6115.1826

Chen Y, Chen L, Lun ATL, Baldoni PL, Smyth GK (2025). edgeR v4: powerful differential analysis of sequencing data with expanded functionality and improved support for small counts and larger datasets. Nucleic Acids Research 53, gkaf018. doi:10.1093/nar/gkaf018

Baldoni PL, Chen L, Smyth GK (2024). Faster and more accurate assessment of differential transcript expression with Gibbs sampling and edgeR v4. NAR Genomics and Bioinformatics 6, lqae151. doi:10.1093/nargab/lqae151

McCarthy DJ, Chen Y, Smyth GK (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Research 40, 4288-4297. doi:10.1093/nar/gks042

Cox DR, Reid N (1987). Parameter orthogonality and approximate conditional inference. Journal of the Royal Statistical Society B 49, 1-18. doi:10.1111/j.2517-6161.1987.tb01422.x

Loader C (2000). Fast and accurate computation of binomial probabilities. Technical note, Lucent Technologies, Murray Hill; the version of 5 February 2002 is available at https://www.r-project.org/doc/reports/CLoader-dbinom-2002.pdf

Phipson B, Lee S, Majewski IJ, Alexander WS, Smyth GK (2016). Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression. Annals of Applied Statistics 10, 946-963. doi:10.1214/16-AOAS920

Robinson MD, Smyth GK (2008). Small-sample estimation of negative binomial dispersion, with applications to SAGE data. Biostatistics 9, 321-332. doi:10.1093/biostatistics/kxm030

Brent RP (1973). Algorithms for Minimization without Derivatives. Prentice-Hall, Englewood Cliffs.

Neumaier A (1974). Rundungsfehleranalyse einiger Verfahren zur Summation endlicher Summen. Zeitschrift fuer Angewandte Mathematik und Mechanik 54, 39-51. doi:10.1002/zamm.19740540106

Giraud L, Langou J, Rozloznik M (2005). The loss of orthogonality in the Gram-Schmidt orthogonalization process. Computers & Mathematics with Applications 50, 1069-1075. doi:10.1016/j.camwa.2005.08.009

Berrut JP, Trefethen LN (2004). Barycentric Lagrange interpolation. SIAM Review 46, 501-517. doi:10.1137/S0036144502417715

Phipson B (2013). Empirical Bayes modelling of expression profiles and their associations. PhD thesis, University of Melbourne. https://hdl.handle.net/11343/38162

Benjamini Y, Hochberg Y (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society B 57, 289-300. doi:10.1111/j.2517-6161.1995.tb02031.x

class cytepool.stats.quasi_likelihood.QuasiLikelihoodFit(fit, dispersion, top_proportion, average_quasi_dispersion, unit_deviance_adjusted, unit_df_adjusted, leverage, deviance_adjusted, df_residual_adjusted, quasi_dispersion, prior_df, prior_quasi_dispersion, posterior_quasi_dispersion, covariate, average_log_cpm, robust, counts, max_iterations)#

Negative binomial fits of every gene with moderated quasi-dispersions.

Rows of the per-gene arrays are genes, in the order of the columns of the counts; arrays with a sample axis are samples x genes. The arrays are read-only. The names in parentheses are edgeR’s.

Variables:

  • fit – The negative binomial fits at the working dispersion dispersion / average_quasi_dispersion: coefficients, shrunk coefficients (coefficients), fitted values, deviances (deviance) and the offsets, dispersions and weights of the fits.

  • dispersion – The negative binomial dispersion (dispersion), as given, or estimated from the top genes: a float, one value per gene, or one per observation.

  • top_proportion – The share of genes the dispersion was estimated from (top.proportion); None when it was given.

  • average_quasi_dispersion – The average quasi-dispersion A (average.ql.dispersion), at least 1.

  • unit_deviance_adjusted – (n_samples, n_features) adjusted unit deviances (unit.deviance.adj).

  • unit_df_adjusted – (n_samples, n_features) their degrees of freedom (unit.df.adj).

  • leverage – (n_samples, n_features) leverages of the fits (leverage); those of fitted values that tend to zero are the limits.

  • deviance_adjusted – (n_features,) adjusted residual deviances, the sums of the adjusted unit deviances.

  • df_residual_adjusted – (n_features,) adjusted residual degrees of freedom (df.residual.adj), fractional; zero for a gene without counts.

  • quasi_dispersion – (n_features,) quasi-dispersions, the ratio of the two; NaN where the degrees of freedom are zero.

  • prior_df – The prior degrees of freedom (df.prior): a float, or one value per gene from a robust fit that finds outliers.

  • prior_quasi_dispersion – The prior quasi-dispersion (s2.prior): one value per gene with a trend, a float otherwise.

  • posterior_quasi_dispersion – (n_features,) posterior quasi-dispersions (s2.post).

  • covariate – (n_features,) covariate of the trend of the prior: by default the average log CPM that also ranks the genes, at the given NB dispersion or at 0.05 when the dispersion is estimated (edgeR’s abundance); None without a trend.

  • average_log_cpm – (n_features,) average log2 counts per million at the negative binomial dispersion, the logCPM of the tests.

  • robust – Whether the prior was estimated robustly.

  • counts – (n_samples, n_features) the counts fitted (counts).

  • max_iterations – The most Newton steps of each fit, as given.

Explained in: Count models

property coefficient_names#

Names of the coefficients, the design’s columns.

class cytepool.stats.quasi_likelihood.QuasiLikelihoodTest(fit, names, contrasts, log_fold_change, average_log_cpm, likelihood_ratio, f_statistic, p_value, adjusted_p_value, df_test, df_total)#

Quasi-likelihood F-tests of coefficients or contrasts, one per gene.

The columns of edgeR’s topTags table are log_fold_change (logFC), average_log_cpm (logCPM), f_statistic (F), p_value (PValue) and adjusted_p_value (FDR). The arrays are read-only.

Variables:

  • fit – The quasi-likelihood fit tested.

  • names – Names of the coefficients or contrasts tested together.

  • contrasts – (n_coefficients, n_tested) weights of each tested contrast on the coefficients; the columns of the identity for coefficients.

  • log_fold_change – (n_features, n_tested) log2 fold changes, each contrast of the shrunk coefficients divided by log 2.

  • average_log_cpm – (n_features,) average log2 counts per million.

  • likelihood_ratio – (n_features,) likelihood ratio statistics: the deviance of the null fit less that of the full fit, both at the working dispersion; at least 0.

  • f_statistic – (n_features,) F-statistics, the likelihood ratios over df_test and the posterior quasi-dispersions.

  • p_value – (n_features,) P(F > f) on df_test and df_total degrees of freedom.

  • adjusted_p_value – (n_features,) Benjamini and Hochberg’s adjusted p-values over the genes.

  • df_test – Numerator degrees of freedom, the rank of contrasts.

  • df_total – (n_features,) denominator degrees of freedom (df.total): the prior plus the adjusted residual degrees of freedom, at most the residual degrees of freedom summed over the genes.

Explained in: Count models

class cytepool.stats.quasi_likelihood.UnitDevianceMoments(mean, variance)#

Mean and variance of the negative binomial unit deviance.

Variables:

  • mean – E[d(Y, mu, phi)] for Y ~ NB(mu, phi).

  • variance – var[d(Y, mu, phi)].

Explained in: Count models

cytepool.stats.quasi_likelihood.average_log_cpm(counts, *, library_sizes=None, norm_factors=None, offsets=None, weights=None, prior_count=2.0, dispersion=0.05)#

Average log2 counts per million of each gene, as edgeR’s aveLogCPM.

Each gene’s samples are fitted as one group by the negative binomial model of fit_negative_binomial_model(), after prior_count times its library size over the mean library size is added to each count and twice that to its library size (edgeR’s addPriorCount). The result is log2 of the fitted count per million: log2(1e6 exp(b)) with b the fitted log mean per unit of library size. It is similar to the log2 of the mean of the counts per million, with larger libraries given more weight, and equal to it for a large dispersion (edgeR’s help page). The prior count keeps genes without counts finite.

Parameters:

  • counts – (n_samples, n_features) counts, finite and non-negative.

  • library_sizes – Library size of each sample, positive; defaults to the total counts.

  • norm_factors – Normalization factor of each sample, positive; defaults to ones. The library sizes are multiplied by them.

  • offsets – Log library sizes, in place of library sizes and normalization factors: one value, one per sample, or one per observation.

  • weights – Prior weights, positive: one value, one per sample, or one per observation.

  • prior_count – Average count added to each observation, at least 0.

  • dispersion – Negative binomial dispersion of the fits: one value, one per gene, or one per observation; 0.05 by default, as in edgeR.

Returns:

(n_features,) average log2 counts per million.

Raises:

  • TypeError – If an argument has the wrong type, or is a numpy masked array.

  • ValueError – As fit_negative_binomial_model(), or if prior_count is negative.

Explained in: Count models

cytepool.stats.quasi_likelihood.common_dispersion(counts, design, *, library_sizes=None, norm_factors=None, offsets=None, weights=None, interval=_DISPERSION_INTERVAL, min_row_sum=_MIN_ROW_SUM, subset=_SUBSET, coefficient_names=None, max_iterations=200, n_threads=None)#

Common negative binomial dispersion of all genes by Cox-Reid adjusted profile likelihood, as edgeR’s estimateGLMCommonDisp.

The dispersion phi maximizes the sum over genes of the adjusted profile log-likelihoods (McCarthy, Chen and Smyth 2012; Cox and Reid 1987)

APL_g(phi) = l_g(phi, beta_g(phi)) - 1/2 log det(X’ W_g X),

with beta_g(phi) the maximum likelihood fit at phi of fit_negative_binomial_model() and W_g the expected information of the fitted values, w mu / (1 + phi mu) with prior weights w. With prior weights each count is the average of w independent negative binomial counts, as for edgeR’s adjustedProfileLik: its log-likelihood is that of w y under NB(w mu, phi / w). The adjustment, which a REML criterion parallels, removes most of the downward bias of the maximum likelihood estimate (Robinson and Smyth 2008). Observations whose fitted values tend to zero, below 1e-12 of the gene’s largest, carry no information and are left out of the determinant, whose remaining factor does not depend on phi.

As edgeR does (its help pages): genes with fewer than min_row_sum counts in all are left out; when more than subset genes remain, a systematic subset of them, evenly spaced in average log CPM, is used (systematicSubset: every s-th gene in increasing average log CPM, s = floor(G / subset), from gene floor(s / 2) + 1; found by probing), which keeps all G genes while G < 2 subset and from subset to about 1.5 subset of them otherwise; and the maximum is searched over interval. The search, Brent’s (1973) method as R’s optimize, locates the maximum to about 1e-7 relative, as closely as the rounding of the likelihood allows; edgeR’s stops within its tolerance of 1e-5. A maximum at an end of the interval is returned with a warning.

Parameters:

  • counts – (n_samples, n_features) counts, finite and non-negative.

  • design – (n_samples, n_coefficients) design matrix of full column rank, an object with matrix and columns such as a Design, or a table.

  • library_sizes – Library size of each sample, positive; defaults to the total counts.

  • norm_factors – Normalization factor of each sample, positive; defaults to ones.

  • offsets – Log effective library sizes, in place of library sizes and normalization factors: one value, one per sample, or one per observation.

  • weights – Prior weights, positive: one value, one per sample, or one per observation.

  • interval – The smallest and largest dispersion searched, 0 <= low < high.

  • min_row_sum – Genes with fewer counts in all are left out; at least 0.

  • subset – The number of genes beyond which a systematic subset of them is used, as above; at least 1.

  • coefficient_names – Names of the design columns.

  • max_iterations – The most Newton steps of each fit.

  • n_threads – The most threads of the fits; None uses numba’s setting.

Returns:

The common dispersion.

Raises:

  • TypeError – If an argument has the wrong type, or is a numpy masked array.

  • ValueError – As fit_negative_binomial_model(), if an option is out of range, or if no gene has min_row_sum counts.

Explained in: Count models

cytepool.stats.quasi_likelihood.fit_quasi_likelihood(counts, design, *, dispersion=None, library_sizes=None, norm_factors=None, offsets=None, weights=None, trend=True, robust=False, prior_count=0.125, top_proportion=None, coefficient_names=None, max_iterations=200, n_threads=None)#

Fit negative binomial models with empirical Bayes quasi-dispersions, as edgeR’s glmQLFit.

edgeR 4’s default method (legacy = FALSE), described in the module description: the negative binomial dispersion phi is given or estimated as the Cox-Reid common dispersion of the most highly expressed genes; every gene is fitted at the working dispersion phi / A, with A the average quasi-dispersion; its quasi-dispersion is its bias-adjusted deviance over its adjusted residual degrees of freedom; and the quasi-dispersions are squeezed towards a prior that follows the average log CPM, by limma’s estimator for unequal degrees of freedom. edgeR’s help page says that a dispersion that is given should come from a global model, common or trended, not from each gene on its own: the quasi-dispersions carry the variation between genes.

Parameters:

  • counts – (n_samples, n_features) counts, finite and non-negative, integer or not.

  • design – (n_samples, n_coefficients) design matrix of full column rank with fewer columns than samples, an object with matrix and columns such as a Design, or a table. A Design with random effects is refused.

  • dispersion – Negative binomial dispersion, finite, from 0 to 1e4: one value, one per gene (n_features,) or one per observation (n_samples, n_features). By default it is estimated as the common dispersion (common_dispersion()) of the top top_proportion of the genes by average log CPM.

  • library_sizes – Library size of each sample, positive; defaults to the total counts.

  • norm_factors – Normalization factor of each sample, such as TMM factors from tmm_factors(), positive; defaults to ones.

  • offsets – Log effective library sizes, in place of library sizes and normalization factors: one value, one per sample, or one per observation.

  • weights – Prior weights, positive and finite: one value, one per sample, or one per observation; see the module description for how edgeR uses them.

  • trend – Whether the prior quasi-dispersion follows the average log CPM (edgeR’s abundance.trend), or (n_features,) covariate values for it to follow instead (covariate.trend).

  • robust – Whether to estimate the prior robustly, so that genes whose quasi-dispersions are outliers count less (Phipson et al. 2016), as edgeR’s User’s Guide recommends.

  • prior_count – Average count added to each observation for the shrunk coefficients, whose differences are the log fold changes of the tests; at least 0.

  • top_proportion – Share of the genes, by average log CPM, that the dispersion is estimated from, in (0, 1]. By default chooseLowessSpan(G sqrt(n - k), small.n = 20, min.span = 0.02) for G genes, n samples and k coefficients, which falls as genes and residual degrees of freedom are added, as in edgeR. Only for an estimated dispersion.

  • coefficient_names – Names of the design columns.

  • max_iterations – The most Newton steps of each fit, at least 1.

  • n_threads – The most threads of the fits; None uses numba’s setting. The result does not depend on it.

Returns:

The fits, the dispersions and the quasi-dispersions.

Raises:

  • TypeError – If an argument has the wrong type, or is or holds a numpy masked array.

  • ValueError – As fit_negative_binomial_model(), if the design leaves no residual degrees of freedom, an option is out of range, top_proportion is given with a dispersion, no gene has adjusted residual degrees of freedom, the fitted values times the prior weights exceed 1e100, or fewer than two genes have a positive quasi-dispersion on at least 0.01 adjusted residual degrees of freedom, so that limma’s prior cannot be estimated.

Used in: Bulk RNA-seq

Explained in: Choosing a method, Count models

cytepool.stats.quasi_likelihood.negative_binomial_unit_deviance_moments(means, dispersion)#

Mean and variance of the unit deviance of a negative binomial count.

For Y ~ NB(mu, phi), with var(Y) = mu + phi mu^2 (Poisson at phi = 0), the unit deviance d(Y, mu, phi) of negative_binomial_unit_deviance() has mean E = sum_y P(Y = y) d(y, mu, phi) and variance V = sum_y P(Y = y) (d(y, mu, phi) - E)^2. edgeR 4 approximates each unit deviance by the scaled chi-squared variable with these two moments, s chi^2(k) with s = V / (2 E) and k = 2 E^2 / V (Chen et al. 2025): k is near 1 for large means, larger for means near 0.5, and tends to 0 with the mean.

The probabilities are formed in Loader’s (2000) deviance form, P(Y = y) = sqrt(r / (2 pi y (y + r))) exp(-d / 2 + delta(y + r) - delta(r) - delta(y)) with r = 1 / phi and delta Stirling’s error, which has no cancellation for large counts or large r. Counts are summed one by one while that takes at most 600 terms; otherwise the counts below 64 are summed and the rest of the sum is the integral of its continuation, with Gregory’s end correction from the forward differences at 64, by Gauss-Legendre panels on the log of the count. The moments agree to about 1e-14 relative with sums at 30 digits for means up to 300, with the gamma limit of large means and with the large-mean expansion of the Poisson moments. edgeR’s approximations of them differ by up to about 1e-2 (see the module description).

Parameters:

  • means – Means, finite and non-negative; at most 1e100.

  • dispersion – Dispersions, finite and non-negative; at most 1e4.

Returns:

The means and variances, as float64 arrays of the broadcast shape, or floats when both arguments are scalars. A mean of zero gives zero for both.

Raises:

  • TypeError – If an argument is not real-valued, or is a numpy masked array.

  • ValueError – If an argument is negative or not finite, a mean exceeds 1e100 or a dispersion 1e4, or the shapes do not broadcast.

Explained in: Count models

cytepool.stats.quasi_likelihood.quasi_likelihood_f_test(fit, *, coefficients=None, contrasts=None, names=None)#

Quasi-likelihood F-tests that coefficients or contrasts are zero, as edgeR’s glmQLFTest.

The null model is fitted at the working dispersion of fit, without the tested coefficients, or on the design X N whose coefficients are those that the tested contrasts C leave free (N spans the vectors b with C’ b = 0, so that every full-rank parametrisation of the same null model gives the same deviance). The likelihood ratio statistic, the null deviance less the full one, over the rank r of the contrasts and the posterior quasi-dispersion, is referred to the F distribution on r and df_total degrees of freedom (see the module description). Several coefficients or contrasts are tested together, as edgeR’s matrix of contrasts is; test each named contrast in its own call for a table per contrast. A gene without counts gets an F-statistic of 0 and a p-value of 1 (NaN when the fit has at most two genes, whose quasi-dispersions are not moderated).

Contrasts that name their coefficients are matched to them by name, as in contrasts_fit(): a Contrasts (make_contrasts, Design.contrasts or its select), a pandas DataFrame or Series, or a mapping of name to weights; other weights are taken in the order of the coefficients.

Parameters:

  • fit – A fit from fit_quasi_likelihood().

  • coefficients – The coefficients to test, by name or position; by default the last one, as in edgeR. Not with contrasts.

  • contrasts – Weights on the coefficients of the contrasts to test: a vector, a (n_coefficients, n_contrasts) matrix, or named contrasts as above.

  • names – Names of the contrasts when contrasts carries none; by default contrast1, contrast2, ....

Returns:

The tests, with the columns of edgeR’s topTags table.

Raises:

  • TypeError – If fit is not a QuasiLikelihoodFit, or an argument has the wrong type or is a numpy masked array.

  • ValueError – If both coefficients and contrasts are given, a coefficient is unknown or repeated, the contrasts are not finite, have the wrong number of rows or are all zero, the tests would leave the null model without coefficients (as edgeR’s), or names are invalid.

Used in: Bulk RNA-seq

Explained in: Choosing a method, Count models