dream#

voom_with_dream_weights() gives voom’s precision weights from mixed models of every gene, fitted by maximum likelihood, as variancePartition’s voomWithDreamWeights does. moderated_fixed_effect_tests() moderates the Satterthwaite or Kenward-Roger t-tests of fixed_effect_tests() by empirical Bayes, as variancePartition’s eBayes does for dream fits, adds topTable’s z-statistic, and tests the contrasts at once by the moderated F-statistic. dream_voom() runs these steps on every cell type of a cube, with dream’s choice of the degrees of freedom.

dream: voom weights from mixed models, and moderated tests of their fixed effects.

dream (Hoffman and Roussos 2021) is limma-voom with a linear mixed model in place of the linear model, for designs with repeated measures: several samples of a donor, tissues of one patient, time courses. The donor (or another grouping factor) gets a random intercept, whose variance is estimated for every gene, and contrasts of the fixed effects are tested with Kenward and Roger’s (1997) degrees of freedom in small studies and Satterthwaite’s (1946) in larger ones. dreamlet (Hoffman et al. 2026) runs it on the pseudobulk profiles of every cell type. The steps, as variancePartition runs them (voomWithDreamWeights, dream, eBayes and topTable):

Precision weights, voom_with_dream_weights(). voom’s (Law et al. 2014): the log2 counts per million are fitted, the square roots of the residual standard deviations are smoothed by LOWESS against the average log2 counts, and each observation gets the inverse fourth power of the trend at its fitted value (cytepool.stats.voom). The fit is the mixed model of every feature, by maximum likelihood, and the fitted values include the predicted random effects.
Fits and tests: fit_mixed_model() by REML with those weights, and fixed_effect_tests(), with the degrees of freedom of dream’s ddf: by default ("adaptive") Kenward and Roger’s for fewer than 20 samples, with their adjusted standard errors, and Satterthwaite’s otherwise.
Moderation, moderated_fixed_effect_tests(). The residual variance of every feature, that of an observation of its average weight, is shrunk towards a prior by limma’s empirical Bayes (Smyth 2004) on the approximate residual degrees of freedom of the mixed model, tr((I - H)’ (I - H)) for its hat matrix H (MixedModelFit.residual_df). The moderated t-statistic divides an estimate by its standard error with the posterior variance in place of the residual variance, and has the degrees of freedom of the test plus the prior’s. Several contrasts are tested at once by the moderated F-statistic t’ R^-1 t / r, limma’s, for their moderated t-statistics and the correlation matrix R of their estimates, with r the number of linearly independent contrasts.

Conventions that the documentation leaves open were found by running variancePartition 1.36.2, with limma 3.62.1 and lme4 1.1-38, as a black box:

voomWithDreamWeights fits the mixed models by maximum likelihood (its weights agree with those of maximum likelihood fits to lme4’s tolerance, and differ from those of REML fits by up to 12% in an example), and reads the trend at lme4’s fitted values, which include the predicted random effects. Without random effects it is limma’s voom.
It divides each feature’s prior weights by their mean before the fits, for sample weights and observation weights alike. By default it returns the trend’s weights alone (rescaleWeightsAfter = FALSE), as limma’s voom does; with rescaleWeightsAfter = TRUE, the trend’s weights times the scaled prior weights, which are the weights a fit needs and the ones returned here. They do not change when the prior weights are multiplied by a constant.
dream also fits each feature with its weights divided by their mean, so its sigma is the residual standard deviation of an observation of the feature’s average weight: residual_variance / mean_weights of a fit to the weights as given. With voom’s weights this is about the variance that the trend predicts at the feature’s expression, so the variances that eBayes moderates keep the mean-variance trend that the weights removed; trend=True lets the prior follow it. Its sigma therefore differs from limma’s with the same weights; with every variance component at zero, the fit is limma’s weighted fit, and the moderated tests are limma’s for the weights of each feature divided by their mean.
eBayes moderates those variances on the approximate residual degrees of freedom rdf, with limma’s squeezeVar and its default estimator: the likelihood estimator (Chen et al. 2025) when the degrees of freedom differ between features, as they usually do, and the method of moments otherwise (legacy = NULL). The moderated t-statistics and the B-statistics are those of limma’s eBayes for these variances on these degrees of freedom; the p-values use Satterthwaite’s degrees of freedom plus the prior’s, which, unlike limma’s total degrees of freedom, are not capped at the sum of the residual degrees of freedom.
topTable adds z.std, sign(t) times the standard normal quantile of half the two-sided p-value: the z-statistic with the p-value of the t-statistic, so that tests on different degrees of freedom can be compared. It computes it from the p-value floored at 1e-300.
dream(ddf = "adaptive"), the default, uses Kenward and Roger’s method for 19 samples and Satterthwaite’s for 20; the help page says “<= 20 samples”. With ddf = "Kenward-Roger" it tests with lmerTest’s contest, which calls pbkrtest: Kenward and Roger’s adjustment is computed with every observation at the average weight (see cytepool.stats.mixed_models), the standard errors are the adjusted ones, and eBayes moderates the tests as it does Satterthwaite’s. A gene whose adjusted covariance matrix of the coefficients is not positive definite or has a condition number above 1e6, by R’s estimate kappa, is tested by Satterthwaite’s method instead (variancePartition’s NEWS for 1.31.16 states the rule; the bound was found by probing). The condition number depends on the units of continuous covariates.
topTable with several coefficients gives the coefficients, AveExpr, F, P.Value, adj.P.Val and F.std: F is limma’s moderated F-statistic, from the moderated t-statistics and each gene’s covariance matrix of the coefficients, Kenward and Roger’s adjusted one under their method, and F.std the F-statistic on infinite denominator degrees of freedom with the same p-value.

Deliberate differences from variancePartition, reported by validation/dream.py and validation/kenward_roger_and_f_tests.py:

The LOWESS span is chosen from the number of features by default, as voom() and limma’s voom (since limma 3.66) choose it. voomWithDreamWeights uses 0.5, which span=0.5 reproduces.
Beyond a p-value of 1e-300, the z-statistic is computed from the logarithm of the tail probability, so it stays exact where topTable stops at 37.07.
A feature whose mixed model has no minimum, which the fixed and random effects reproduce, is not fitted (see MixedModelFit): it is left out of the trend, its weights are read at its values, which the effects reproduce, and its tests are NaN. lme4 fits it with a residual standard deviation near zero, which enters the trend, and dream tests it.
topTable computes the p-values of an F-test of every gene on the total degrees of freedom of the first gene for the first coefficient listed (found by probing), so that they change with the order of the genes and of the coefficients. Here each gene has its own: F is the mean of the squares of r uncorrelated moderated t-statistics, those of the eigenvectors of R, whose degrees of freedom are those of their contrasts plus the prior’s, and the denominator degrees of freedom combine them as lmerTest combines those of its F-tests (2 E / (E - r), with E = sum nu_m / (nu_m - 2); see cytepool.stats.mixed_models). With one contrast, the F-test is the moderated t-test. F.std uses the number of coefficients where it differs from their rank; here the rank, the F-test’s numerator degrees of freedom.
dream_voom decides which genes fall back to Satterthwaite’s method from the exact condition number, where dream uses R’s estimate: they agree except within about 2% of 1e6.

Random slopes are not implemented yet; variancePartition’s variance partitioning is cytepool.variance_partitioning.

References

Hoffman GE, Roussos P (2021). Dream: powerful differential expression analysis for repeated measures designs. Bioinformatics 37, 192-201. doi:10.1093/bioinformatics/btaa687

Hoffman GE, Lee D, Bendl J, Prashant NM, Hong A, Casey C, Alvia M, Shao Z, Argyriou S, Therrien K, Venkatesh S, Voloudakis G, Haroutunian V, Fullard JF, Roussos P (2026). Efficient differential expression analysis of large-scale single-cell transcriptomics data using Dreamlet. Nature Communications 17, 9597. doi:10.1038/s41467-026-75680-8

Law CW, Chen Y, Shi W, Smyth GK (2014). voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology 15, R29. doi:10.1186/gb-2014-15-2-r29

Kenward MG, Roger JH (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics 53, 983-997. doi:10.2307/2533558

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

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

Satterthwaite FE (1946). An approximate distribution of estimates of variance components. Biometrics Bulletin 2, 110-114. doi:10.2307/3002019

class cytepool.stats.dream.ModeratedFixedEffectTests(tests, residual_variance_at_mean_weight, residual_df, prior_variance, prior_df, posterior_variance, t, df_total, p_value, lods, z, coefficient_prior_variance, proportion, prior_estimator, f_statistic, f_df1, f_df2, f_p_value, standardized_f_statistic)#

dream’s moderated t-tests of contrasts of the fixed effects of mixed models.

Rows of the per-feature arrays are features and columns are contrasts, as in tests. The arrays are read-only.

Variables:

tests – The tests moderated, with their fit, contrasts and Satterthwaite’s degrees of freedom.
residual_variance_at_mean_weight – (n_features,) the variances moderated: the residual variance of an observation of each feature’s average weight, fit.residual_variance / fit.mean_weights (dream’s sigma^2).
residual_df – (n_features,) their degrees of freedom, the fit’s approximate residual degrees of freedom (dream’s rdf).
prior_variance – The prior variance (s2.prior): a float, or one value per feature with a trend.
prior_df – The prior degrees of freedom (df.prior): a float, possibly infinite, or one value per feature from a robust fit that finds outliers.
posterior_variance – (n_features,) posterior variances (s2.post).
t – (n_features, n_contrasts) moderated t-statistics: the estimates over their standard errors times (residual_variance_at_mean_weight / posterior_variance)^(1/2).
df_total – (n_features, n_contrasts) their degrees of freedom: Satterthwaite’s plus the prior’s (df.total).
p_value – Two-sided p-values of t on df_total degrees of freedom.
lods – Log-odds that each contrast is non-zero (the B-statistic), as limma’s eBayes gives them for these variances on residual_df.
z – Signed z-statistics with the same p-values, sign(t) times the standard normal quantile of half the p-value (topTable’s z.std).
coefficient_prior_variance – (n_contrasts,) unscaled prior variance of a non-zero contrast (var.prior).
proportion – The assumed proportion of non-zero contrasts behind lods.
prior_estimator – The estimator of the prior: "moments" or "likelihood".
f_statistic – (n_features,) moderated F-statistics testing every contrast of tests at once: t’ R^-1 t / f_df1, for the moderated t-statistics t and the correlation matrix R of their estimates (topTable’s F for these contrasts).
f_df1 – The numerator degrees of freedom of the F-tests, the number of linearly independent contrasts.
f_df2 – (n_features,) their denominator degrees of freedom; see the module docstring.
f_p_value – (n_features,) p-values of f_statistic on f_df1 and f_df2 degrees of freedom.
standardized_f_statistic – (n_features,) the F-statistic on f_df1 and infinite degrees of freedom with the same p-value (topTable’s F.std), so that tests on different degrees of freedom compare.

Explained in: Mixed models, REML and dream

cytepool.stats.dream.moderated_fixed_effect_tests(tests, *, proportion=0.01, stdev_coef_lim=(0.1, 4.0), trend=False, robust=False, winsor_tail_p=(0.05, 0.1), prior_estimator=None)#

Empirical Bayes moderated t-tests of mixed model contrasts, as dream’s eBayes.

variancePartition’s eBayes for dream fits, with topTable’s z-statistic: the residual variance of an observation of each feature’s average weight is shrunk towards a prior on the approximate residual degrees of freedom by limma’s empirical Bayes (moderated_statistics(), with its estimators and options), each contrast’s t-statistic uses the posterior variance, and its degrees of freedom are those of the tests (Satterthwaite’s or Kenward and Roger’s) plus the prior’s. The moderated F-statistic tests every contrast of tests at once, as topTable does for several coefficients. See the module docstring.

Parameters:

tests – Tests from fixed_effect_tests(), of a REML fit as dream makes them by default. To F-test some of the contrasts, test those alone.
proportion – Assumed proportion of non-zero contrasts, in (0, 1).
stdev_coef_lim – Limits for the standard deviation of a non-zero contrast, as for moderated_statistics.
trend – Whether the prior variance follows the mean expression of the features (eBayes(trend = TRUE)), or an (n_features,) covariate.
robust – Whether to estimate the prior robustly (robust = TRUE).
winsor_tail_p – Winsorized tail proportions of the robust moment estimator.
prior_estimator – "moments", "likelihood" or None for limma’s default (legacy = NULL), as for moderated_statistics.

Returns:

The moderated tests.

Raises:

TypeError – If tests is not a FixedEffectTests, or an option has the wrong type or is a numpy masked array.
ValueError – If tests has no contrasts, an option is invalid, or no feature has a fit.

Explained in: Mixed models, REML and dream

cytepool.stats.dream.voom_with_dream_weights(counts, design, *, library_sizes=None, norm_factors=None, span=None, weights=None, method='ml', n_threads=None)#

voom’s precision weights from the mean-variance trend of mixed model fits.

variancePartition’s voomWithDreamWeights: voom (see voom()), with the linear model replaced by the mixed model of design, fitted to every feature by maximum likelihood, and the trend read at the fitted values with the predicted random effects; the module docstring gives the conventions. Continue with fit_mixed_model(result.log_cpm, design, weights=result.weights), fixed_effect_tests and moderated_fixed_effect_tests().

Parameters:

counts – (n_samples, n_features) counts, as for voom: finite and non-negative, at least two features.
design – A Design with random intercepts, as fit_mixed_model() takes; for a design without random effects use voom.
library_sizes – Library size of each sample; defaults to the total counts.
norm_factors – Normalization factor of each sample, such as TMM factors; defaults to ones.
span – The span of the LOWESS trend, in (0, 1]; by default chosen from the number of features, as voom chooses it. voomWithDreamWeights uses 0.5.
weights – Prior weights: (n_samples,) sample weights, such as each pseudobulk sample’s number of cells, or (n_samples, n_features) observation weights; finite and positive. Each feature’s weights are divided by their mean, then weight the fits behind the trend, and multiply the weights returned.
method – "ml", maximum likelihood, as voomWithDreamWeights fits, or "reml".
n_threads – The most threads for the fits, as for fit_mixed_model.

Returns:

The log2 counts per million, their weights (the trend’s times the scaled prior weights), the effective library sizes, and the trend with its span.

Raises:

TypeError – If an argument has the wrong type, or is or holds a numpy masked array.
ValueError – If counts is not a 2-D array of finite non-negative values with at least two features, the design has no random effects or cannot be fitted (as for fit_mixed_model), its samples are not the rows of counts, a library size or normalization factor is not positive, span is not in (0, 1], a prior weight is not finite and positive, fewer than two features with counts can be fitted, or the trend is not positive where weights are read from it.

Explained in: Mixed models, REML and dream

dream

Contents

dream#