Variance partitioning#

variance_fractions() partitions the variance of every gene among the variables of a Design and the residuals, as variancePartition’s fitExtractVarPartModel does: each gene is fitted with a linear mixed model, by maximum likelihood, whose random intercepts are the categorical variables (donor, condition, batch) and whose fixed effects are the continuous ones (age), with precision weights such as voom’s, and the fraction of a variable is its variance component, or the variance of its fitted effect, over their sum. variance_partitioning() does it for every cell type of a cube, after the steps of dream_voom() up to the precision weights, and records the samples each fit used. plot_variance_fractions() and plot_gene_variance_fractions() draw the fractions as variancePartition’s plotVarPart and plotPercentBars do.

Variance partitioning: the share of each gene’s variance that each variable explains.

Which variables drive the differences between samples, and which must be adjusted for? For every gene, variance partitioning (Hoffman and Schadt 2016) fits one linear mixed model with every variable of the design, so that each is assessed together with all the others, and reports the fraction of the gene’s variance that each variable explains, and the residual fraction. It is what variancePartition’s fitExtractVarPartModel computes.

variance_fractions() partitions the variance of every feature of a samples x features matrix of log-expression values, and variance_partitioning() of every group (cell type) of a cube, after the steps of dream_voom() up to voom’s precision weights.

The model:

Every feature is fitted with the model of cytepool.stats.mixed_models,

y = X beta + Z_1 b_1 + … + Z_K b_K + e, b_k ~ N(0, sigma_k^2 I), e ~ N(0, sigma^2 W^-1),

where categorical variables, such as the donor, the condition or the batch, are random intercepts, (1 | donor), and continuous variables, such as age, are fixed effects (~ age). Without random effects it is a linear model, fitted by weighted least squares. The variances are estimated by maximum likelihood, as fitExtractVarPartModel does by default (REML = FALSE): its documentation explains that REML, which estimates the random effects after removing the fixed effects, biases the fractions of variance, while maximum likelihood estimates the fixed and random effects together. method="reml" fits by REML instead.

The fractions:

Each variable contributes a variance, and its fraction is that variance over their sum:

a random effect k, its variance component sigma_k^2;
the residual, sigma^2 / mean(w): the residual variance of an observation of the feature’s average weight (sigma^2 itself without weights);
a fixed-effect term j, the variance of its fitted contribution X_j beta_j over the samples, the post hoc variance of Hoffman and Schadt (2016). The variances of fitted values are weighted, sum_i w_i (v_i - v_w)^2 / sum_i w_i about the weighted mean v_w, with the denominator of a maximum likelihood variance (n without weights). With random effects, the terms share the variance of all fitted fixed effects, X beta, in proportion to the variances of their own contributions; without random effects, each term has the variance of its own contribution. The two agree when the terms are uncorrelated over the samples.

The fractions of a feature sum to one, and a random effect whose variance is estimated at zero (a singular fit in lme4) has a fraction of exactly zero. They do not change when the weights of a feature are multiplied by a constant, when the levels of a factor are relabelled, or when the values of a feature are multiplied by a constant or shifted.

Conventions found by probing variancePartition 1.36.2, with lme4 1.1-38, as a black box, where its documentation leaves them open:

The variances of the fixed effects are those above: weighted, with the denominator of a maximum likelihood variance, and shared among the terms in a mixed model but not in a linear model (calcVarPart of lmer and lm fits).
fitExtractVarPartModel divides each gene’s precision weights by their mean, as dream does, so the residual variance is that of an observation of average weight and the fractions do not depend on the scale of the weights. Without random effects the residual variance is the weighted mean of the squared residuals, the maximum likelihood estimate.
A design with random effects must give every categorical variable as a random effect: a factor among its fixed effects is refused, as variancePartition refuses it (“Must model either _all_ or _no_ categorical variables as random effects”), and so is a formula without an intercept. Without random effects, factors are fixed effects of the linear model.
variancePartition lists the random effects in alphabetical order, then the fixed effects; here they are in the order of the design.

Deliberate differences from variancePartition, reported by validation/variance_partitioning.py:

A feature that the fixed and random effects reproduce, such as one that is constant within every donor, has no maximum likelihood estimate (see MixedModelFit) and its fractions are NaN. lme4 fits it with a residual variance near zero, which gives it a fraction near one for the donor.
method="reml" without random effects takes the REML estimate of the residual variance, the residual sum of squares over n - p. variancePartition fits such a design with lm and takes the residual sum of squares over n whatever REML says.

Limits of this version:: Random slopes, which variancePartition uses to partition the variance within subsets of the samples, such as (tissue + 0 | donor), generalized linear mixed models, and bootstrap confidence intervals of the fractions (varPartConfInf) are not implemented yet.

References

Hoffman GE, Schadt EE (2016). variancePartition: interpreting drivers of variation in complex gene expression studies. BMC Bioinformatics 17, 483. doi:10.1186/s12859-016-1323-z

Bates D, Maechler M, Bolker B, Walker S (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software 67(1), 1-48. doi:10.18637/jss.v067.i01

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

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

Patterson HD, Thompson R (1971). Recovery of inter-block information when block sizes are unequal. Biometrika 58, 545-554. doi:10.1093/biomet/58.3.545

class cytepool.variance_partitioning.VarianceFractions(fractions, variances, variables, random_effects, fixed_effects, method, converged, fit)#

The fraction of every feature’s variance that each variable explains.

Rows of the arrays are features, in the order of the columns of y, and columns are variables, in the order of variables. The arrays are read-only. A feature without a fit has NaN fractions: one that the fixed and random effects reproduce, as MixedModelFit describes.

Variables:

fractions – (n_features, n_variables) fraction of each feature’s variance that each variable explains; each row sums to one.
variances – (n_features, n_variables) the variances the fractions are shares of (calcVarPart(fit, returnFractions = FALSE)): the variance component of each random effect, the variance of the fitted contribution of each fixed-effect term, and the residual variance of an observation of the feature’s average weight, on the scale of the values.
variables – The names of the variables: the random effects, as lme4 names their grouping factors (such as "donor" or "sample:donor"), in the order of the design; the fixed-effect terms, as the formula labels them (such as "age"); and "residuals", last.
random_effects – The names of the random effects among variables.
fixed_effects – The names of the fixed-effect terms among variables.
method – "ml" or "reml", how the variances were estimated.
converged – (n_features,) False where a feature was not fitted, or its mixed model fit did not converge (its fractions are those where the fit stopped).
fit – The fits: a MixedModelFit, or for a design without random effects a LinearModelFit.

Explained in: Mixed models, REML and dream

class cytepool.variance_partitioning.VariancePartitioningResult(variables, summary, samples, tables, provenance)#

The fractions of variance of every group of a cube, from variance_partitioning().

Tables hold, for every group that was fitted, one row per feature that passed the expression filter, in the cube’s feature order: feature_id and, when the cube’s feature table has one, symbol; then the fraction of the feature’s variance that each variable explains, one column per variable, named as in variables (NaN for a feature whose mixed model has no maximum, which the fixed and random effects reproduce).

The summary has one row per group: group_id; status, "fitted" or why the group was not fitted ("no_samples", "no_residual_df", "fixed_effects_not_estimable", "random_effects_not_estimable", "no_features", "too_few_features", "empty_library" or "failed"); reason, the explanation (empty when fitted); n_samples_used and n_samples_left_out; and n_features, the number of features fitted (0 where the group was not fitted).

The sample table has one row per group and design sample, as for DifferentialExpressionResult: group_id, sample_id, n_cells (when the cube records cell numbers), used, reason, and the library_size and norm_factor that entered the fit; with prior_weights, the prior_weight of each sample.

Built by variance_partitioning(). Immutable.

Raises:

TypeError – If variables is not a sequence of names, summary or samples not a Table, tables not a mapping of group id to Table, or provenance not a mapping that JSON can write.
ValueError – If provenance holds a NaN or an infinity.

Explained in: Mixed models, REML and dream

property variables#: The variables, in the order of the columns of the tables, residuals last.

property groups#: The ids of the groups analysed, in order.

property summary#: The status, samples and features of every group.

property samples#: Every design sample of every group, used in its fit or left out, and why.

property tables#: The table of fractions of every group that was fitted, by group id.

property settings#: The settings of the analysis, as a new dict (also in the provenance).

property provenance#

A JSON-serialisable record of what produced every number, as a new dict.

It holds the method; the CytePool, numpy, numba and Python versions; the cube, by its fingerprint, digests and provenance; the layer read; the design (formula, columns, levels, random effects, and the values of the variables of both for every sample); the settings; and, for every group, its status, the samples that entered its fit in the order they were fitted, those left out with the reason, the number of features fitted, the span of its trend (voom_span), the number of features whose mixed model has no maximum (features_not_fitted) and of those whose fit stopped without converging (features_not_converged), and with prior_weights the prior weights of its samples. The cost of the cube’s record is that of limma_voom()’s (“Provenance”).

table(group)#

Return the table of fractions of one group.

Parameters:: group – A group id.
Raises:: KeyError – If the group was not analysed or not fitted; the message says why.

Explained in: Mixed models, REML and dream

cytepool.variance_partitioning.variance_fractions(y, design, *, weights=None, method='ml', n_threads=None)#

The fraction of every feature’s variance that each variable of a design explains.

variancePartition’s fitExtractVarPartModel: every feature (column of y) is fitted with the design’s fixed effects and random intercepts, by maximum likelihood as lme4’s lmer(REML = FALSE) fits it (see fit_mixed_model()), and its variance is partitioned among the random effects, the fixed-effect terms and the residual, as the module documentation describes. Categorical variables are given as random effects and continuous ones as fixed effects, as in Design(samples, "~ age", random="(1 \| donor) + (1 \| condition)"). A design without random effects is fitted by weighted least squares, and its factors are fixed effects.

Parameters:

y – (n_samples, n_features) expression values on a log scale, such as the log-CPM values of voom. Every value must be finite.
design – A Design with an intercept, whose random effects are random intercepts, such as random="(1 | donor) + (1 | batch)", and whose fixed effects are then continuous covariates; or a design without random effects. Its fixed effects must have full column rank, and its random effects must pass design.check().
weights – Precision weights, (n_samples,) for sample weights or (n_samples, n_features) for observation weights, such as those of voom_with_dream_weights(): the residual variance of an observation is sigma^2 / weight. They must be finite and positive; only their ratios within a feature matter.
method – "ml" for maximum likelihood, variancePartition’s default, or "reml" for restricted maximum likelihood.
n_threads – The most threads for fitting features, as for fit_mixed_model; the result does not depend on it.

Returns:

The fractions and variances of every feature, with the fits.

Raises:

TypeError – If design is not a Design, an input is not real-valued, method is not a string, n_threads is not an integer, or an input is or holds a numpy masked array.
ValueError – If shapes disagree, y is not finite, a weight is not finite and positive, method is unknown, the formula has no intercept, a design with random effects has a factor among its fixed effects or a random effect that is not a random intercept, the fixed effects are not of full column rank, a random effect fails the design checks, the fixed and random effects together span every sample, or a variable is named "residuals" or two variables have the same name.

Example

>>> import numpy as np
>>> from cytepool.design import Design
>>> samples = {"donor": ["D1", "D1", "D2", "D2", "D3", "D3", "D4", "D4"],
...            "age": [31.0, 31.0, 45.0, 45.0, 52.0, 52.0, 60.0, 60.0],
...            "condition": ["ctrl", "stim"] * 4}
>>> design = Design(samples, "~ age",
...                 random="(1 | donor) + (1 | condition)")
>>> y = np.random.default_rng(1).normal(5.0, 1.0, (8, 3))
>>> result = variance_fractions(y, design)
>>> result.variables
('donor', 'condition', 'age', 'residuals')

Explained in: Mixed models, REML and dream

cytepool.variance_partitioning.variance_partitioning(cube, design, *, groups=None, min_cells=10, min_counts=0, expression_filter=True, normalization='tmm', span=None, prior_weights=None, method='ml', n_threads=None)#

Partition the variance of every feature in every group of a cube.

What drives the differences between samples in each cell type: for every group of the cube, the fraction of each feature’s variance that each variable of the design explains, as dreamlet runs variancePartition’s fitExtractVarPartModel on the voomWithDreamWeights output of every cell type. For each group, the samples are left out, the features filtered and the counts normalized as for dream_voom(); voom’s precision weights come from the mean-variance trend of the design’s mixed model, fitted by maximum likelihood (voom_with_dream_weights()), or from limma’s voom for a design without random effects; and the fractions of every feature are those of variance_fractions() of its log2 counts per million with these weights.

Every group’s mixed model uses the random effects of the samples left in it, as in dream_voom. A group is not fitted, and its status says why, when no sample is left, its samples cannot estimate a random effect (one donor left, or one sample per donor), they cannot separate the fixed effects (a covariate constant over them), fewer than two features with counts pass the filter, or the trend cannot give weights. A run never stops for one group.

Parameters:

cube – A cube with a counts layer.
design – A Design with an intercept, matched to the cube by sample_id as for limma_voom(), whose categorical variables are random intercepts and continuous variables fixed effects, such as Design(samples, "~ age", random="(1 \| donor) + (1 \| condition)"); or a design without random effects.
groups – A group id or a list of them; all groups by default.
min_cells – The fewest cells a sample needs in a group to enter its fit; ignored for cubes without cell numbers.
min_counts – The fewest counts a sample needs in a group, over every feature of the cube, to enter its fit; 0 by default.
expression_filter – True for edgeR’s filterByExpr defaults with the group’s fixed effects, a mapping of changed settings, or False.
normalization – "tmm" for TMM normalization factors, or None.
span – The LOWESS span of voom’s trend, in (0, 1], for every group; None chooses it in each group from the number of features, as voom does. variancePartition’s voomWithDreamWeights uses 0.5.
prior_weights – "n_cells" to give each sample its number of cells in the group as a prior weight, which the precision weights multiply, as dreamlet weights samples; or None. Needs a cube that records cell numbers.
method – "ml" (maximum likelihood, variancePartition’s default) or "reml", for the fits behind the fractions; the fits behind the precision weights are by maximum likelihood, as voomWithDreamWeights fits them.
n_threads – The most threads for the fits of the features, at least 1; None keeps numba’s setting. The results do not depend on it.

Returns:

The tables of fractions of every fitted group, a summary of every group, the samples used and left out, the settings and a provenance record.

Raises:

TypeError – If cube is not a Cube, design not a Design, an option has the wrong type, or groups, min_cells, expression_filter or n_threads is or holds a numpy masked array.
ValueError – If the design cannot be partitioned (as for variance_fractions()), has no sample_id column, names samples the cube does not have or names a sample twice, a variable is named "feature_id" or "symbol", an option is invalid, or prior_weights is "n_cells" and the cube does not record cell numbers.
KeyError – If a group is not in the cube.

Explained in: Mixed models, REML and dream

Variance partitioning

Contents

Variance partitioning#