Empirical Bayes moderation

Empirical Bayes moderation of feature-wise variances, as in limma’s eBayes.

A linear model fit gives each feature g a residual variance s_g^2 on d_g degrees of freedom. The hierarchical model of Smyth (2004) assumes

s_g^2 | sigma_g^2 ~ sigma_g^2 chi^2(d_g) / d_g, 1 / sigma_g^2 ~ chi^2(d0) / (d0 s0^2),

so that s_g^2 ~ s0^2 F(d_g, d0) marginally. The prior (s0^2, d0) is estimated from all features, and each variance is replaced by its posterior value s~_g^2 = (d0 s0^2 + d_g s_g^2) / (d0 + d_g), the inverse of the posterior mean of 1 / sigma_g^2 (Smyth 2004). A moderated t-statistic, which uses s~_g in place of s_g, follows a t distribution on d0 + d_g degrees of freedom under the null hypothesis; an infinite d0 gives every feature the prior variance.

Estimation of the prior, by one of two estimators (prior_estimator):

• The method of moments ("moments", limma’s legacy=TRUE): fit_f_distribution() matches the mean and variance of e_g = log s_g^2 - psi(d_g / 2) + log(d_g / 2), whose variance is psi’(d0 / 2) plus the mean of psi’(d_g / 2) (Smyth 2004, section 6.2). When the variance of the e_g is no larger than chi-squared sampling explains, d0 is infinite and s0^2 is the mean variance. With a covariate (limma-trend uses the mean log-expression) the moments are taken about a natural cubic spline fit of e_g, and s0^2 becomes a trend (Phipson et al. 2016, section 10.1). fit_f_distribution_robust() matches the moments of the log variances Winsorized at their 5% and 90% quantiles to those of a Winsorized log F variable, computed by Gauss-Legendre quadrature, so that hypervariable features do not lower d0. Each feature then gets its own prior df, smaller for outliers (Phipson et al. 2016). With a covariate, the log variances are first detrended by LOWESS.

• Weighted moments and profile likelihood ("likelihood", limma’s fitFDistUnequalDF1): fit_f_distribution_likelihood() estimates log s0^2 by a weighted mean of the e_g, which gives less weight to variances on fewer degrees of freedom, and d0 by maximising the likelihood of the variances under s0^2 F(d_g, d0) with that scale (Chen et al. 2025). It suits degrees of freedom that differ between features, or are small or fractional. With a covariate the weighted mean becomes a weighted LOWESS trend; its robust variant downweights outliers by prior weights instead of Winsorizing.

By default (prior_estimator=None) the likelihood estimator is used when the positive degrees of freedom differ between features, as with missing values or zero observation weights, and the method of moments otherwise: limma’s default, legacy=NULL, since limma 3.61.8.

Statistics: moderated_statistics() computes moderated t and F statistics, their p-values and the log-odds of differential expression (the B-statistic of Loennstedt and Speed 2002, in the form of Smyth 2004, section 5); treat() tests whether coefficients exceed a threshold in absolute value (McCarthy and Smyth 2009).

Results match limma 3.62.1 (fitFDist, fitFDistRobustly, fitFDistUnequalDF1, squeezeVar, eBayes and treat, with legacy NULL, TRUE or FALSE). Conventions that limma’s documentation leaves open were settled by running it as a black box; comments mark them “found by probing”. The robust likelihood estimator deliberately differs from limma’s where some features have no usable variance or prior weights are given (see fit_f_distribution_likelihood()).

References

Smyth GK (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology 3, article 3. doi:10.2202/1544-6115.1027

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

McCarthy DJ, Smyth GK (2009). Testing significance relative to a fold-change threshold is a TREAT. Bioinformatics 25, 765-771. doi:10.1093/bioinformatics/btp053

Loennstedt I, Speed TP (2002). Replicated microarray data. Statistica Sinica 12, 31-46.

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

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

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

Olver FWJ, Olde Daalhuis AB, Lozier DW, Schneider BI, Boisvert RF, Clark CW, Miller BR, Saunders BV, Cohl HS, McClain MA, eds (2026). NIST Digital Library of Mathematical Functions, release 1.2.8 of 2026-09-15. https://dlmf.nist.gov/

class cytepool.stats.empirical_bayes.FDistributionFit(scale, df2)

A scaled F distribution fitted to variances: s_g^2 ~ scale F(df_g, df2).

Variables:

• scale – The prior variance s0^2: a float, or one value per feature when the fit has a covariate.

• df2 – The prior degrees of freedom d0, possibly infinite.

Explained in: Empirical Bayes moderation and TREAT

class cytepool.stats.empirical_bayes.LikelihoodFDistributionFit(scale, df2, df2_shrunk, df2_outlier)

A scaled F distribution fitted by weighted moments and profile likelihood.

Variables:

• scale – The prior variance s0^2: a float, or one value per feature when it follows a covariate (NaN where the covariate is NaN).

• df2 – The prior degrees of freedom d0, between 2 and 9998.

• df2_shrunk – (n,) prior degrees of freedom of each feature: df2, except for the outliers of a robust fit, which get smaller values.

• df2_outlier – Prior degrees of freedom of the most extreme outlier of a robust fit; NaN when the fit is not robust or finds no outlier.

Explained in: Empirical Bayes moderation and TREAT

class cytepool.stats.empirical_bayes.ModeratedFit(fit, prior_variance, prior_df, posterior_variance, df_total, t, p_value, lods, coefficient_prior_variance, f_statistic, f_p_value, f_df1, proportion, prior_estimator)

Moderated t and F statistics and log-odds of a linear model fit.

Rows of the per-feature arrays are features and columns are the coefficients (or contrasts) of fit. The arrays are read-only.

Variables:

• fit – The linear model fit.

• prior_variance – The prior variance s0^2 (limma’s s2.prior): a float, or one value per feature with a trend or with at most two features.

• prior_df – The prior degrees of freedom d0 (df.prior): a float, possibly infinite, or one value per feature from a robust fit (from the likelihood estimator only when it finds outliers, as in limma).

• posterior_variance – (n_features,) posterior variances (s2.post).

• df_total – (n_features,) degrees of freedom of the moderated t-statistics, df_residual + prior_df capped at the residual degrees of freedom summed over all features (df.total).

• t – (n_features, n_coefficients) moderated t-statistics.

• p_value – Two-sided p-values of t on df_total degrees of freedom.

• lods – Log-odds that each coefficient is non-zero (the B-statistic), on the natural log scale.

• coefficient_prior_variance – (n_coefficients,) unscaled prior variance v0 of a non-zero coefficient, which is N(0, v0 sigma_g^2) (var.prior).

• f_statistic – (n_features,) moderated F-statistics testing all coefficients equal to zero.

• f_p_value – P-values of f_statistic on f_df1 and df_residual + prior_df degrees of freedom; limma does not cap the denominator degrees of freedom here (found by probing).

• f_df1 – Numerator degrees of freedom of the F-tests: the rank of the correlation matrix of the coefficients.

• proportion – The assumed proportion of non-zero coefficients behind lods.

• prior_estimator – The estimator of the prior: "moments" or "likelihood" (see squeeze_variances()); with at most two features neither is used.

Explained in: Empirical Bayes moderation and TREAT

class cytepool.stats.empirical_bayes.RobustFDistributionFit(scale, df2, df2_shrunk, df2_outlier, tail_p_value, prob_outlier)

A scaled F distribution fitted robustly, with prior df for each feature.

Variables:

• scale – The prior variance s0^2: a float, or one value per feature when the fit has a covariate.

• df2 – Prior degrees of freedom of the features that are not outliers.

• df2_shrunk – (n,) prior degrees of freedom of each feature: df2, except for outliers, which get smaller values.

• df2_outlier – Prior degrees of freedom of the most extreme outlier; NaN when df2 is infinite or the moment estimate is returned.

• tail_p_value – (n,) P(F > s_g^2 / s0^2), the probability of a variance as large; NaN for unusable features, and when the moment estimate is returned.

• prob_outlier – (n,) estimated probability that each feature is an outlier, before it is made monotonic in tail_p_value; NaN where tail_p_value is.

Explained in: Empirical Bayes moderation and TREAT

class cytepool.stats.empirical_bayes.SqueezedVariances(posterior_variance, prior_variance, prior_df)

Posterior variances and the prior they were shrunk towards.

Variables:

• posterior_variance – (n,) posterior values of the variances: the inverses of the posterior means of their reciprocals.

• prior_variance – The prior variance: one value per feature when it follows a covariate or with at most two variances, a float otherwise.

• prior_df – The prior degrees of freedom: a float, or one value per feature from a robust fit (from the likelihood estimator only when it finds outliers, as in limma).

Explained in: Empirical Bayes moderation and TREAT

class cytepool.stats.empirical_bayes.TreatFit(fit, lfc, prior_variance, prior_df, posterior_variance, df_total, t, p_value, prior_estimator)

Moderated t-tests relative to a threshold on the absolute coefficients.

Variables:

• fit – The linear model fit.

• lfc – The threshold tau of the null hypothesis |beta| <= tau, on the scale of the coefficients (log2 fold changes for log2 expression).

• prior_variance – The prior variance, as in ModeratedFit.

• prior_df – The prior degrees of freedom, as in ModeratedFit.

• posterior_variance – (n_features,) posterior variances.

• df_total – (n_features,) degrees of freedom of the tests.

• t – (n_features, n_coefficients) moderated t-statistics shrunk towards zero by the threshold: sign(t) max(|t| - tau / se, 0), with se the moderated standard error.

• p_value – P-values of the tests of |beta| <= tau.

• prior_estimator – The estimator of the prior, as in ModeratedFit.

Explained in: Empirical Bayes moderation and TREAT

cytepool.stats.empirical_bayes.fit_f_distribution(variances, df, *, covariate=None)

Fit a scaled F distribution to variances by the method of moments.

Each variance is modelled as scale * F(df_g, df2). Matches limma’s fitFDist. The moments are those of the logs of the variances, each bounded below by 1e-5 times their median (limma does the same). With a covariate, the scale follows a natural cubic spline in it: with n usable features, the spline has 4 degrees of freedom including the intercept for n >= 30, 3 for n from 6 to 29, 2 for n from 3 to 5 and 1 (a constant) below, and never more than the covariate has distinct values; its interior knots are at quantiles of the covariate.

A variance is usable when it is finite and its df is positive. Unusable features do not enter the estimate; with a covariate, their scale is the spline evaluated at their covariate (linear beyond the range of the usable features). A single usable variance gives df2 = 0 and itself as the scale.

Parameters:

• variances – (n,) non-negative variances; NaN marks a missing value.

• df – Degrees of freedom of the variances: one value, or (n,). Values below 1e-15 count as zero.

• covariate – (n,) covariate for a trended scale, such as the mean log-expression. Infinite values are moved one unit beyond the finite range, as limma does; NaN is allowed only for unusable features.

Returns:

The scale and df2, which is infinite when the variances vary no more than chi-squared sampling explains; the scale is then their mean, or the exponentiated trend.

Raises:

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

• ValueError – If shapes disagree, a variance is negative or infinite, a df is negative or not finite, a usable feature has a NaN covariate, or no variance is usable.

Explained in: Empirical Bayes moderation and TREAT

cytepool.stats.empirical_bayes.fit_f_distribution_likelihood(variances, df, *, covariate=None, prior_weights=None, robust=False)

Fit a scaled F distribution by weighted moments and profile likelihood.

limma’s fitFDistUnequalDF1 (limma 3.61.8 and later), which squeezeVar, eBayes and treat use by default when the degrees of freedom of the variances differ. Each variance is modelled as scale * F(df_g, df2). The mean of e_g = log s_g^2 - psi(d_g / 2) + log(d_g / 2) is log s0^2 - psi(d0 / 2) + log(d0 / 2) whatever d_g, so for each d0 the scale is estimated from the weighted mean of the e_g with weights pw_g / psi’(d_g / 2): psi’(d_g / 2) is the sampling variance of log s_g^2, so variances on fewer degrees of freedom count less (pw_g are the prior weights, 1 by default). d0 then maximises the likelihood of the variances, each term weighted by pw_g, under the scaled F distribution with that scale: a profile likelihood (Chen et al. 2025). With a covariate, the weighted mean becomes a weighted trend of the e_g on the covariate.

A variance is usable when it is finite, its df is at least 0.01 and its prior weight is positive; at least two usable variances must be above zero. The conventions of limma 3.62.1 that its documentation leaves open were found by probing it:

• Variances at or below zero are raised to 1e-12 times the median of the positive usable variances.

• d0 minimises minus twice the log-likelihood by Brent’s method, as R’s optimize with its default tolerance, over v = d0 / (d0 + 2) in [1/2, 0.9998]. d0 therefore lies between 2 and 9998: variances that vary no more than chi-squared sampling explains give d0 near 8,135, not infinity, and the scale is that of this d0.

• The trend is limma’s loessFit(e, covariate, weights, span = chooseLowessSpan(n, small.n = 500), iterations = 1) with the weights divided by their upper quartile and kept in [1e-8, 100]. Every feature with a covariate enters it, an unusable one with the smallest weight, 1e-8, its variance (zero when missing) raised as above and a df of 1 when below 0.01. Equal weights give R’s lowess. Unequal weights give a weighted least squares line through at most four features, and through more limma’s weighted LOWESS, whose windows hold a share of the total weight.

• With at most two usable variances above zero the scale is one value whatever the covariate, and a robust fit is the first fit. When some variances are not usable as well, limma returns NaN; here they are left out.

• A feature whose covariate is NaN (allowed only when its variance is not usable) is left out of the trend and gets a NaN scale, where limma stops; infinite covariates are moved one unit beyond the finite range, as for fit_f_distribution(), where limma returns NaN.

Robust fits (robust=True) downweight outliers by prior weights instead of Winsorizing, as limma does:

1. Each usable variance gets the two-sided tail probability of s_g^2 / s0^2 under F(d_g, d0) from the fit above, adjusted by Benjamini and Hochberg’s (1995) method; those below 0.3 are outliers (large or small variances).

2. The fit is repeated with each outlier’s prior weight multiplied by its adjusted probability.

3. pi_g = min(1, p_g / ((r_g - 0.5) / n)), with p_g the upper tail probability of the first fit and r_g the rank of s_g^2 / s0^2 from the largest among the n usable variances, is made non-decreasing in p_g as in fit_f_distribution_robust() and taken as the probability of not being an outlier. df2_outlier takes two steps of d <- d log(1/2) / log P(F > r) (Phipson et al. 2016, section 10.5) for the variance with the smallest p_g, from d = df2 of the second fit with its first p_g, then with its ratio r to the first fit’s scale. df2_shrunk is pi_g df2 + (1 - pi_g) df2_outlier.

Without outliers the first fit is returned, and when every pi_g is 1 the second one, with df2_shrunk equal to df2. This deliberately differs from limma 3.62.1 in two respects. Only the usable variances enter steps 1 and 3, where limma also counts the others, a missing variance as a zero one: it becomes the most extreme outlier and lowers the adjusted probabilities of the others, so missing values change which features are outliers. In a simulation where 40 of 2,000 genes had no residual df, limma’s robust prior df was 11.5% larger with a constant prior, where it also gave 45 genes a smaller prior df as outliers and this function gives every gene the same, and 6% to 10% larger with a trend. Without a covariate the result is limma’s fit of the usable variances alone. And the second fit keeps the prior weights, which limma’s drops, bringing back the features they leave out. The two agree when every variance is usable and there are no prior weights.

With df of 1e7 and more, limma’s likelihood loses precision: its d0 differed from this function’s by up to 1e-5 relative at 1e7 and 5e-2 at 1e12, where this function’s agrees with the likelihood computed to 60 digits.

Parameters:

• variances – (n,) non-negative variances; NaN marks a missing value.

• df – Degrees of freedom of the variances: one value, or (n,). They may be fractional or below 1.

• covariate – (n,) covariate for a trended scale, such as the mean log-expression; NaN is allowed only where the variance is not usable.

• prior_weights – (n,) non-negative weights of the features; 0 leaves a feature out.

• robust – Whether to downweight outliers.

Returns:

The scale, df2, and with robust the prior df of each feature and of the most extreme outlier.

Raises:

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

• ValueError – As fit_f_distribution(), if prior_weights is invalid, or if fewer than two usable variances are above zero (limma returns NA), in the fit or, with robust, in its second fit.

Explained in: Empirical Bayes moderation and TREAT

cytepool.stats.empirical_bayes.fit_f_distribution_robust(variances, df, *, covariate=None, winsor_tail_p=(0.05, 0.1))

Fit a scaled F distribution robustly, with prior df for each feature.

The robust empirical Bayes estimator of Phipson et al. (2016), matching limma’s fitFDistRobustly:

1. Variances are bounded below by 1e-12 times their median. When their df differ, each is mapped to the largest df d by matching tail probabilities under the moment estimate (s0^2, d0): s0^2 Q(P(s^2 / s0^2; d_g, d0); d, d0), with Q the F quantile (section 10.2).

2. With a covariate, the log variances are detrended by LOWESS (span 0.4, three robustness iterations), and the scale is the exponentiated trend times the estimated scale of the detrended values.

3. The log variances are Winsorized at their winsor_tail_p quantiles (R’s type 7), and df2 solves var(Winsorized log variances) = var(log win(F)), with the Winsorized moments of log F computed by 128-point Gauss-Legendre quadrature (sections 3 and 10.3). It is never below the moment estimate (a rule found by probing). log s0^2 is the Winsorized mean minus the mean of the Winsorized log F.

4. Each feature gets a tail probability p_g = P(F > s_g^2 / s0^2), and pi_g = min(1, p_g / ((r_g - 0.5) / n)), with r_g the rank of s_g^2 / s0^2 from the largest, is its probability of not being an outlier; prob_outlier is 1 - pi_g. The pi_g are made non-decreasing in p_g: in increasing order of p_g, the leading values up to the first minimum of their cumulative mean are set to that minimum, and a cumulative maximum follows. df2_shrunk is pi_g df2 + (1 - pi_g) d_outlier, where s0^2 F(d, d_outlier) has the largest ratio as median (two fixed-point iterations, as in limma). When df2 is infinite, the cumulative maximum alone makes the pi_g non-decreasing, and df2_shrunk is infinite where pi_g = 1 and pi_g n d elsewhere, the df of all n variances on the common df d (both found by probing).

Unusable features (NaN variance or zero df) are left out of the fit and get df2; with a covariate, their scale is interpolated linearly in the covariate (constant beyond the range), as limma does. The moment estimate itself is returned, with df2_shrunk equal to its df2, when the number of usable features times the larger tail proportion is below one (fewer than 10 features with the default proportions), and when the moment estimate of df2 is infinite but the Winsorized one is not (both rules found by probing). With more than two usable features, more than half of which have a zero variance, the estimator is undefined and raises an error, as limma does.

limma solves the equation of step 3 only to a tolerance, so its df2 agrees with this exact solution to about 1e-7 relative. Two limitations are shared with limma, whose results agree with these. The quadrature of step 3 is accurate to better than 1e-12 when both df are at least 1 and the tail proportions are the defaults, but not when a df is well below 1 (at df = 0.3 the Winsorized variance of log F is off by 14% to 32%) or, near df = 1, with much smaller tail proportions (2e-4 at 0.01). And with about 20 or fewer features and a covariate, LOWESS passes exactly through some points, whose detrended log variances then tie up to rounding, so that their ranks, and with them prob_outlier and df2_shrunk, can depend on the last bits of the input.

Parameters:

• variances – (n,) non-negative variances; NaN marks a missing value.

• df – Degrees of freedom of the variances: one value, or (n,).

• covariate – (n,) covariate for a trended scale; see fit_f_distribution().

• winsor_tail_p – Proportions of the variances Winsorized in the lower and upper tails, each strictly between 0 and 1/2; one value applies to both.

Returns:

The robust fit.

Raises:

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

• ValueError – As fit_f_distribution(), if winsor_tail_p is invalid, or if more than half of three or more usable variances are zero.

Explained in: Empirical Bayes moderation and TREAT

cytepool.stats.empirical_bayes.moderated_statistics(fit, *, proportion=0.01, stdev_coef_lim=(0.1, 4.0), trend=False, robust=False, winsor_tail_p=(0.05, 0.1), f_statistic=None, prior_estimator=None)

Empirical Bayes moderated t and F statistics and log-odds of a fit.

Matches limma’s eBayes(fit), except that by default the F-statistics of fits with missing values or observation weights are exact rather than limma’s approximation (see below and f_statistic). The residual variances of the fit are squeezed by squeeze_variances(); the moderated t-statistics divide the coefficients by stdev_unscaled times the posterior standard deviations and have df_total degrees of freedom. By default the prior is estimated as limma’s eBayes estimates it (legacy=NULL): by the likelihood estimator when the residual degrees of freedom differ between the features that have some, as with missing values or zero observation weights, and by the method of moments otherwise.

The log-odds (B-statistic) of coefficient j of feature g is (Smyth 2004)

log(p / (1 - p)) + log(r) / 2 + (1 + d) / 2 log((t^2 + d) / (t^2 r + d)),

with p = proportion, d = df_total, r = v / (v + v0), v the squared unscaled standard error and v0 the coefficient’s prior variance; with an infinite prior df the last term is its limit t^2 (1 - r) / 2 (found by probing). v0 is estimated from the largest |t| (Smyth 2004, section 6.3): with G the number of finite t-statistics, each of the ceiling(p G / 2) largest gives an estimate, clipped to [lim1^2, lim2^2] / median(prior variance) from stdev_coef_lim, and v0 is their mean. The proportion used there is at least the share of the largest values, and t-statistics on unequal df are first mapped to the largest df by matching tail probabilities (found by probing).

The moderated F-statistic tests all coefficients of a feature equal to zero: t’ R^+ t / r, with R the correlation matrix of the coefficients and r its rank, counting eigenvalues above 1e-8 times the largest. f_statistic="limma" uses the correlation of fit.cov_coefficients for every feature, as limma does; "exact" uses each feature’s own covariance matrix when the fit holds them (missing values or observation weights), which is the correct F-statistic (Smyth 2004, section 7). The two agree when one covariance matrix serves every feature. A feature with a NaN coefficient has a NaN F-statistic.

Parameters:

• fit – A fit from fit_linear_model() or contrasts_fit().

• proportion – Assumed proportion of non-zero coefficients, in (0, 1).

• stdev_coef_lim – Lower and upper limits for the standard deviation of a non-zero coefficient, used to bound v0; 0 <= lower <= upper, and the upper limit may be infinite.

• trend – Whether the prior variance follows the mean expression of the features (limma-trend), or an (n_features,) covariate to use instead. Features without observed values may have a NaN covariate.

• robust – Whether to estimate the prior robustly.

• winsor_tail_p – Winsorized tail proportions of the robust estimator.

• f_statistic – “exact” or “limma”; by default “exact” unless the fit’s standard errors are limma’s approximation from contrasts_fit(method="limma"), whose F-statistics can only be limma’s.

• prior_estimator – "moments", "likelihood" or None, as in squeeze_variances(); the result records the estimator used.

Returns:

The moderated statistics.

Raises:

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

• ValueError – If an option is invalid, no feature has residual degrees of freedom, or “exact” F-statistics are asked of a fit whose standard errors are limma’s approximation.

Used in: Quick start

Explained in: Choosing a method, Empirical Bayes moderation and TREAT, Precision weights

cytepool.stats.empirical_bayes.squeeze_variances(variances, df, *, covariate=None, robust=False, winsor_tail_p=(0.05, 0.1), prior_estimator=None)

Shrink variances towards a prior by empirical Bayes.

Matches limma’s squeezeVar. The prior is estimated by the method of moments, fit_f_distribution() or, when robust, fit_f_distribution_robust(); or by weighted moments and profile likelihood, fit_f_distribution_likelihood(), robustly when robust. Each variance is replaced by its posterior value (d0 s0^2 + d_g s_g^2) / (d0 + d_g), the inverse of the posterior mean of its reciprocal (Smyth 2004). A feature with zero df, or an infinite prior df, gets the prior variance; a NaN variance with positive df stays NaN. With at most two variances there is no moderation: the prior df is 0 and each variance is its own prior, as in limma 3.61.12 and later. A robust fit gives each feature its own prior df: always with the method of moments, and with the likelihood estimator when it finds outliers, as in limma. The likelihood estimator counts a variance on zero df as missing, as limma’s squeezeVar does.

Parameters:

• variances – (n,) non-negative variances; NaN marks a missing value.

• df – Degrees of freedom of the variances: one value, or (n,).

• covariate – (n,) covariate for a trended prior variance.

• robust – Whether to estimate the prior robustly.

• winsor_tail_p – Winsorized tail proportions of the robust moment estimator; the likelihood estimator does not use them, as in limma.

• prior_estimator – "moments" (limma’s legacy=TRUE), "likelihood" (legacy=FALSE), or None (the default) for limma’s default legacy=NULL: the likelihood estimator when the positive degrees of freedom are not all equal, and the method of moments otherwise.

Returns:

The posterior variances, the prior variance and the prior df.

Raises:

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

• ValueError – As fit_f_distribution() or fit_f_distribution_likelihood(), or if an option is invalid.

Explained in: Empirical Bayes moderation and TREAT

cytepool.stats.empirical_bayes.treat(fit, *, lfc=None, fc=None, trend=False, robust=False, winsor_tail_p=(0.05, 0.1), prior_estimator=None)

Moderated t-tests that coefficients exceed a threshold in absolute value.

Tests H0: |beta| <= tau against |beta| > tau (McCarthy and Smyth 2009), matching limma’s treat(fit). With the moderated t-statistic t, its standard error se and delta = tau / se, the p-value is the worst case over the null interval, P(T > |t| - delta) + P(T > |t| + delta) with T on df_total degrees of freedom; it is at least the ordinary moderated p-value and equals it for tau = 0. The reported t-statistics are sign(t) max(|t| - delta, 0). The variances are moderated as by moderated_statistics(). Unlike limma, which reports t = 0 for a coefficient that is NaN, and for an all-zero contrast tested with a zero threshold, the t-statistic is NaN there.

Give the threshold as a log2 fold change lfc or as a fold change fc, with lfc = log2(fc); the default is fc = 1.2, as in limma. The threshold is |lfc|, so fc = 0.8 and fc = 1.25 are the same test. Typical values are fc of 1.1, 1.2 or 1.5: estimated fold changes must be clearly above the threshold to be significant.

Parameters:

• fit – A fit from fit_linear_model() or contrasts_fit() of log2 values.

• lfc – Threshold as a log2 fold change, finite.

• fc – Threshold as a fold change, positive and finite.

• trend – As in moderated_statistics().

• robust – Whether to estimate the prior robustly.

• winsor_tail_p – Winsorized tail proportions of the robust estimator.

• prior_estimator – "moments", "likelihood" or None, as in squeeze_variances(); the result records the estimator used.

Returns:

The tests relative to the threshold.

Raises:

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

• ValueError – If both lfc and fc are given, either is invalid, or no feature has residual degrees of freedom.

Explained in: Choosing a method, Empirical Bayes moderation and TREAT