Linear models

Linear models#

Linear models fitted to many features at once.

Every feature (gene) of a samples x features matrix is fitted with the same design by ordinary, weighted or generalised least squares, and the fitted coefficients can be re-expressed as contrasts. The results match limma’s lmFit and contrasts.fit. For feature g with observed responses y_g,

E(y_g) = X a_g, var(y_g) = sigma_g^2 W_g^(-1/2) R W_g^(-1/2),

where W_g is the diagonal matrix of the feature’s precision weights and R is the identity, or the correlation matrix of a known correlation between samples of the same block (Smyth 2004; Smyth, Michaud and Scott 2005). A feature with missing values or weights of its own is fitted to its observed samples by a Householder QR decomposition of its own (Golub and Van Loan 2013).

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

Smyth GK, Michaud J, Scott HS (2005). Use of within-array replicate spots for assessing differential expression in microarray experiments. Bioinformatics 21, 2067-2075. doi:10.1093/bioinformatics/bti270

Golub GH, Van Loan CF (2013). Matrix Computations, 4th edition. Johns Hopkins University Press, Baltimore.

class cytepool.stats.linear_models.LinearModelFit(coefficients, stdev_unscaled, sigma, df_residual, cov_coefficients, mean_expression, coefficient_names, design, cov_coefficients_per_feature=None, contrasts=None, exact_standard_errors=True)#

Least squares fits of one design to many features.

Rows of the per-feature arrays are features, in the order of the columns of y; columns are coefficients, or contrasts after contrasts_fit(). The arrays are read-only views; the constructor does not copy them.

Variables:

  • coefficients – (n_features, n_coefficients) estimates. NaN marks a coefficient that the observed values of a feature cannot estimate.

  • stdev_unscaled – (n_features, n_coefficients) unscaled standard errors: the standard errors are stdev_unscaled * sigma[:, None].

  • sigma – (n_features,) residual standard deviations, NaN when df_residual is zero.

  • df_residual – (n_features,) residual degrees of freedom: the number of samples in a feature’s fit minus the rank of the design over them. Stored as floats.

  • cov_coefficients – (n_coefficients, n_coefficients) unscaled covariance matrix of the coefficients. It is the covariance matrix of every feature when cov_coefficients_per_feature is None and exact_standard_errors is True. Otherwise it is the covariance matrix of the design without weights, (X' R^-1 X)^-1 with R the block correlation (or the identity), transformed by any contrasts: what limma reports, and where contrasts_fit(method="limma") takes its correlations from.

  • mean_expression – (n_features,) unweighted mean of the observed values of each feature (limma’s Amean), including values whose weight is zero. The sums are accumulated with error compensation and divided with one rounding, which gives R’s rowMeans to the last bit on log-expression data.

  • coefficient_names – Names of the coefficients or contrasts.

  • design – (n_samples, n_design_coefficients) design matrix.

  • cov_coefficients_per_feature – (n_features, n_coefficients, n_coefficients) unscaled covariance matrix of each feature, or None. Rows and columns of coefficients that a feature cannot estimate are NaN.

  • contrasts – (n_design_coefficients, n_coefficients) matrix mapping the design’s coefficients to coefficients, or None when the fit holds the design’s own coefficients.

  • exact_standard_errors – False when stdev_unscaled are limma’s approximation, from contrasts_fit(method="limma") of a fit with per-feature covariance matrices; such a fit has none.

Raises:

  • TypeError – If coefficient_names is a string, exact_standard_errors is not a bool, or a field is a numpy masked array.

  • ValueError – If shapes disagree, names are not unique, or a fit marked exact has neither per-feature covariance matrices nor standard errors that cov_coefficients implies.

Explained in: Designs, linear models and contrasts

cytepool.stats.linear_models.contrasts_fit(fit, contrasts, *, names=None, method='exact')#

Re-express a fit in terms of contrasts of its coefficients.

For a contrast matrix C whose column c_j defines contrast j, the new coefficients of feature g are C’ a_g. With method="exact" their unscaled variances are c_j’ V_g c_j, with V_g the unscaled covariance matrix of the feature. method="limma" reproduces limma’s contrasts.fit, which does not use V_g: it combines the feature’s unscaled standard errors u_g with the correlation matrix P of cov_coefficients, as sum_kl C_kj C_lj u_gk u_gl P_kl.

The two methods agree when the fit factorised one design for every feature (no missing values; no weights, or positive sample weights), and for contrasts of uncorrelated coefficients, such as group means of a one-way design without blocks. Otherwise limma’s standard errors are approximate, because P then comes from the whole design without weights rather than from each feature’s own samples and weights. On a design whose coefficients are correlated, such as group comparisons adjusted for a covariate, they can be off by a factor of two or more; and with sample weights, a single missing value anywhere, or a single zero weight, makes them approximate for every feature.

A contrast with a non-zero multiplier on a coefficient that is NaN for a feature is NaN for that feature, as in limma. This includes contrasts that the feature’s data do determine although the coefficients they combine are not estimable, such as group C minus group B when the reference group is missing; fit the group means without an intercept to obtain those.

Contrasts that name the coefficients they weigh are matched to fit.coefficient_names by name, whatever their order: a Contrasts (from make_contrasts or Design.contrasts) by its columns, and a pandas DataFrame or Series by its index. They must name every coefficient of the fit exactly once. Weights without such names (numpy arrays, lists, the vectors that indexing a Contrasts returns, and pandas objects with the default index 0, 1, …, n - 1) are taken in the order of fit.coefficient_names; a pandas index of other labels that are not strings is an error. limma’s contrasts.fit always takes the rows in order, and only warns when their names differ from the coefficient names (found by probing).

Parameters:

  • fit – A fit from fit_linear_model() or from this function.

  • contrasts – The contrasts, one per column of a (n_coefficients, n_contrasts) matrix: a Contrasts, or any object with matrix, names and columns attributes laid out the same way; a pandas DataFrame with one contrast per column; a mapping of contrast name to vector or pandas Series; a matrix; or a single vector or Series.

  • names – Contrast names, when contrasts is not a mapping or a Contrasts. Defaults to the column names of a DataFrame or the name of a Series, else contrast1, contrast2, .... Names rename the DataFrame’s columns in order, so none may be the name of another column.

  • method – "exact" or "limma".

Returns:

A fit whose coefficients are the contrasts. sigma, df_residual, mean_expression and design are unchanged, and contrasts maps the design’s coefficients to the new ones. The exact method carries per-feature covariance matrices over. The limma method drops them, as they would contradict its standard errors, which it marks as approximate unless one covariance matrix served every feature.

Raises:

  • TypeError – If fit is not a LinearModelFit, contrasts is not real-valued, it is a design rather than contrasts, or contrasts or names is or holds a numpy masked array.

  • ValueError – If contrasts is not finite or has the wrong number of rows, names coefficients that are not exactly the fit’s, the names are invalid, method is unknown, or the exact method is asked of a fit whose standard errors are approximate.

Used in: Quick start

Explained in: Choosing a method, Designs, linear models and contrasts

cytepool.stats.linear_models.fit_linear_model(y, design, *, weights=None, block=None, correlation=None, coefficient_names=None, n_threads=None)#

Fit the same linear model to every feature by least squares.

Each column of y is regressed on design: by ordinary least squares, by weighted least squares with precision weights, and by generalised least squares when block and correlation set a known correlation between samples of the same block. Results match limma’s lmFit.

NaN in y marks a missing value of one feature, and a zero weight leaves an observed value out of that feature’s fit (it still counts in mean_expression). Each feature is fitted to its remaining samples; its df_residual is their number minus the rank of the design over them, and coefficients they cannot estimate are NaN. A coefficient is estimable when its design column, over those samples, is not a linear combination of the other columns. limma instead drops the later of collinear columns and reports the rest, so it can return a value for a coefficient that the data do not determine, such as the intercept when every sample of the reference group is missing; such values are NaN here, and all estimable values agree with limma.

Parameters:

  • y – (n_samples, n_features) log-scale expression values. NaN marks a missing value; infinite values are rejected.

  • design – (n_samples, n_coefficients) design matrix of full column rank, an object with matrix and columns attributes such as a Design, or a table with columns; column names name the coefficients. A Design with random effects is refused, as they would be dropped, unless its only random effect is a random intercept, such as (1 | donor), block groups the samples as it does and correlation is at least zero; the message says what can be fitted instead.

  • weights – Precision weights, (n_samples,) for sample weights or (n_samples, n_features) for observation weights (as from voom). They must be finite and non-negative wherever y is observed.

  • block – One label per sample; samples with the same label are correlated. Give together with correlation. With a correlation of at least zero this is a random intercept for each block, with the correlation known.

  • correlation – Correlation between samples of the same block, in (-1, 1). With blocks of up to m samples it must exceed -1 / (m - 1), so that the correlation matrix is positive definite.

  • coefficient_names – Names of the design columns, in order. Defaults to the design’s column names when it has them, else x1, x2, .... They may rename the design’s columns, but not give a column the name of another: contrasts are matched to coefficients by name.

  • n_threads – The most threads for fitting features one at a time. None uses numba’s setting for the calling thread (see numba.set_num_threads). The result does not depend on it, and with one thread no parallel region runs.

Returns:

The fit. Without missing values, and with no weights or positive sample weights, one factorisation serves every feature. Otherwise each feature is fitted on its own, in parallel on up to n_threads threads, and the fit stores the covariance matrix of every feature: n_features * n_coefficients**2 floats.

Raises:

  • TypeError – If an input is not real-valued, a label is not comparable, or an input is or holds a numpy masked array (mark missing values of y with NaN).

  • ValueError – If shapes disagree, y has infinite values, design has random effects that the fit would drop, is not finite or is not of full column rank (the message names the aliased columns), a weight is negative or not finite where y is observed, block and correlation are invalid, or coefficient_names move a design column’s name to another column.

Used in: Quick start

Explained in: Choosing a method, Designs, linear models and contrasts, Precision weights