Sample quality weights

Sample quality weights from a heteroscedastic linear model, as in limma’s arrayWeights.

Samples differ in quality: pooled from fewer cells, sequenced less deeply, or taken from an outlying donor, some follow the linear model less closely than others. Ritchie et al. (2006) estimate a weight for each sample from the residuals of every gene, and Liu et al. (2015) extend the model to groups of samples and combine it with voom.

For gene g and sample j, with design rows x_j and prior weights w_gj (one when none are given), the model is

E(y_gj) = x_j’ b_g, var(y_gj) = sigma_g^2 exp(u_j) / w_gj,

with the log variance factors u = Z gamma shared by every gene. Z is an orthonormal basis of the columns of the variance design projected off the constant vector, so u sums to zero and only the span of the variance design matters (found by probing: limma gives the same weights for any parametrisation of it, with or without an intercept). The default variance design gives each sample a variance of its own, and variance_groups one per group. The sample quality weights are exp(-u_j), whose product is one.

gamma is estimated by residual maximum likelihood (REML) with the gene variances sigma_g^2 profiled out. For a weighted least squares fit of gene g at weights w_gj exp(-u_j) over its samples, with standardised residuals e_gj, leverages h_gj, residual degrees of freedom d_g and s_g^2 = sum_j e_gj^2 / d_g, the score (twice the gradient of the REML log likelihood) and the information are

U_g = Z’ z_g, z_gj = e_gj^2 / s_g^2 - (1 - h_gj), I_g = Z’ (V_g - (1 - h_g)(1 - h_g)’ / d_g) Z,

where V_g has elements (delta_jk - h_g,jk)^2, h_g,jk being those of the hat matrix (Smyth 2002). The weights are squeezed toward equality by a prior equivalent to prior_n genes at equal weights: it adds prior_n Z’(exp(-u) - 1) to the score and prior_n Z’Z, the identity, to the information (found by probing). Two algorithms:

• "reml": Fisher scoring from equal weights, gamma <- gamma + d with the step d = I^-1 U, U and I summed over the genes and the prior, until the criterion d’U / ((G + prior_n) q) is below tolerance, for q columns of Z and G genes; where the criterion rises instead, the iterations stop before that step (found by probing).

• "gene_by_gene": one pass through the genes in order (Ritchie et al. 2006). Each gene adds its information, with V_g replaced by the diagonal matrix of its row sums, diag(1 - h_g), to an information matrix that starts at prior_n Z’Z, and moves gamma by that matrix inverted times its score, evaluated at the current gamma. The result depends on the order of the genes.

• "auto": REML without prior weights or missing values, and gene by gene otherwise, as limma chooses.

Results match limma 3.62.1. Conventions that its documentation leaves open, found by probing (running limma as a black box):

• Gene by gene, a gene is skipped where it has fewer than two residual degrees of freedom, or a residual variance below 1e-15 at the current weights. Its samples are those with an observed value and a positive weight.

• REML uses only the genes with every value observed and every prior weight positive. Without prior weights, it leaves out the genes whose residual variance at equal weights is below 1e-15. With prior weights, every gene adds its information, a gene whose residual variance at the current weights is below 1e-15 adds no score, and the iterations stop when d’U / ((G + prior_n)^2 q) is below tolerance, going on where the criterion rises.

• Every weight is one when there are fewer than two genes, when the design leaves fewer than two residual degrees of freedom, or, for REML without prior weights, when fewer than two genes have a positive residual variance; CytePool warns.

REML stops where limma stops, so that the weights are limma’s. With the default tolerance they are then within about half a percent of the REML estimate without prior weights, and up to about 15% from it with prior weights, where the criterion falls below the tolerance after one or two iterations (as with the voom weights of voom_with_quality_weights()), or where the criterion rises first. A smaller tolerance converges further where the iterations converge. A step that would give weights that are not finite is not taken: REML stops before it, with a warning.

Deliberate deviations, where limma returns values without meaning or fails:

• variance_design and variance_groups together are an error; limma uses var.group and ignores var.design.

• prior_n must be at least zero; limma accepts a negative prior.

• A design that is not of full column rank is an error, as in fit_linear_model(); limma drops aliased columns.

• REML with fewer than two complete genes gives weights of one, with a warning; limma fails when one is left.

• An information matrix that is singular is an error: some sample carries no information about its variance, such as the only sample of a group of the design, and prior_n is zero or too small to make up for it. limma’s REML returns weights that rounding error decides.

References

Ritchie ME, Diyagama D, Neilson J, van Laar R, Dobrovic A, Holloway A, Smyth GK (2006). Empirical array quality weights in the analysis of microarray data. BMC Bioinformatics 7, 261. doi:10.1186/1471-2105-7-261

Liu R, Holik AZ, Su S, Jansz N, Chen K, Leong HS, Blewitt ME, Asselin-Labat ML, Smyth GK, Ritchie ME (2015). Why weight? Modelling sample and observational level variability improves power in RNA-seq analyses. Nucleic Acids Research 43, e97. doi:10.1093/nar/gkv412

Smyth GK (2002). An efficient algorithm for REML in heteroscedastic regression. Journal of Computational and Graphical Statistics 11, 836-847. doi:10.1198/106186002871

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

cytepool.stats.sample_weights.sample_quality_weights(y, design, *, weights=None, variance_design=None, variance_groups=None, prior_n=10.0, method='auto', max_iterations=50, tolerance=1e-5)

Estimate a quality weight for each sample from how well it fits the model.

limma’s arrayWeights (Ritchie et al. 2006; Liu et al. 2015); the module docstring gives the model, the algorithms and the conventions. Samples whose values have larger residuals across genes get lower weights. Multiply them into the precision weights of a fit, as voom_with_quality_weights() does for voom.

Parameters:

• y – (n_samples, n_features) log-expression values, such as the log2 counts per million of voom(). NaN marks a missing value; infinite values are rejected.

• design – (n_samples, n_coefficients) design matrix of full column rank, or a Design without random effects. A column of ones treats the samples as replicates.

• weights – Prior weights: (n_samples,) sample weights, (n_samples, n_features) observation weights (such as voom’s), or (1, n_features) feature weights. Finite and non-negative where y is observed; a zero weight leaves the value out. A vector is always sample weights. limma’s arrayWeights takes a matrix only.

• variance_design – (n_samples, n_columns) design of the variance model, a vector for one column, or a Design. Only its span matters, and its constant part is ignored. Defaults to a variance for each sample.

• variance_groups – One label per sample: samples with the same label share a weight. Another way to give a variance design; give at most one of the two.

• prior_n – The number of genes at equal weights that the prior is worth; larger values squeeze the weights more toward one. At least 0.

• method – "reml", "gene_by_gene", or "auto", which chooses REML when there are no prior weights and no missing values.

• max_iterations – The most REML iterations, at least 1.

• tolerance – REML stops when the criterion of the module docstring falls below it, as limma’s does. Positive; lower it to converge further than limma.

Returns:

(n_samples,) positive weights whose product is one.

Warns:

RuntimeWarning – If every weight is one for lack of genes or of residual degrees of freedom, if REML leaves out genes with missing values or zero weights, or if REML stops without converging: at max_iterations, or where its criterion rises or its step is not finite.

Raises:

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

• ValueError – If shapes disagree, y has infinite values, the design is not of full column rank or has random effects, a weight is negative or not finite where y is observed, the variance design is not finite, does not vary between samples, is given twice or is a Design with random effects, variance_groups has fewer than two groups, an option is out of range, or the information matrix is singular (some sample carries no information about its variance, and prior_n is too small).

Explained in: Precision weights