Gene-Wise Linear Models for Microarray Expression

Why Model Expression One Gene at a Time?

A microarray experiment measures expression for many genes across a set of biological samples. After appropriate preprocessing, the measured values are commonly analyzed on a logarithmic scale. The central modeling question is how expression for a given gene varies with sample characteristics recorded in the experiment, such as treatment, batch, or a clinical covariate. A gene-wise linear model gives a separate set of coefficients for each gene while using the same sample-level information to describe them.

The term gene-wise is important: the response values and coefficients belong to one gene at a time, but the design matrix records the same samples and predictors for all genes. This separates the expression measurements from the experimental structure. The model is a framework; the design matrix specifies which sample characteristics are represented, while the coefficient vector describes their relationship with a particular gene’s expected expression.

What the Model Says

For sample \(i\) and gene \(g\), the model can be written as

Here \(x_{ij}\) is the value of design-matrix column \(j\) for sample \(i\), and \(\beta_{jg}\) is the coefficient for that column and gene \(g\). The expected expression, conditional on the design, is the linear predictor \(\sum_j x_{ij}\beta_{jg}\), provided the errors have conditional mean zero. Each gene has its own coefficient vector; the model does not require the coefficients to be equal across genes. In contrast, \(X\) is shared because the samples and their recorded predictors are the same when each gene is analyzed.

For example, if expression is recorded as \(\log_2\) intensity, a coefficient of 1 on a predictor that changes by one unit corresponds to a one-unit difference in expected log-base-2 expression. On the original intensity scale, a one-unit log-base-2 difference corresponds to a factor of \(2\). This describes the model’s scale relationship; it does not establish that a predictor causes a change in expression.

Error Assumptions and Coefficient Estimation

The model’s error assumptions describe variation around the expected expression, conditional on the design matrix. The most basic mean assumption is \(E(\boldsymbol{\varepsilon}_g\mid X)=\mathbf{0}\). It says that the model’s linear predictor is the conditional mean, so systematic shifts left out of the predictor are not being placed in the errors. A common classical assumption for a gene is equal error variance across samples and zero covariance between different samples:

This covariance condition means that each sample error has variance \(\sigma_g^2\), and that errors for distinct samples have zero covariance. It does not, by itself, imply independence. Independence is a separate assumption when it is required. If errors are jointly normal conditional on \(X\), then zero covariance between components does imply their independence; without joint normality, zero covariance alone is not enough. In the standard normal-error version of the model, \(\boldsymbol{\varepsilon}_g\mid X\sim N_n(\mathbf{0},\sigma_g^2 I_n)\). Normality is an additional distributional assumption, not a consequence of the covariance equation.

The variance \(\sigma_g^2\) is gene-specific: different genes may have different error variances. The model does not require errors from different genes to be independent. Their expression measurements may be correlated across genes, even when the model is fit gene by gene. Treating each gene’s coefficients separately is therefore distinct from asserting that genes are biologically or statistically unrelated.

For the ordinary least-squares coefficient estimate to be uniquely defined, \(X\) must have full column rank: its \(p\) columns cannot be exact linear combinations of one another. Under that condition, the estimate for gene \(g\) is

This formula follows by minimizing the residual sum of squares. Write \(Q(\boldsymbol{\beta}_g)=(\mathbf{y}_g-X\boldsymbol{\beta}_g)^\mathsf{T}(\mathbf{y}_g-X\boldsymbol{\beta}_g)\). Differentiating with respect to the coefficient vector and setting the derivative to zero gives the normal equations:

When \(X\) has full column rank, \(X^\mathsf{T}X\) is invertible, so solving the normal equations yields the stated estimate. If the columns are not full rank, the data do not uniquely determine all \(p\) coefficients; a different parameterization or additional constraints are needed. Full rank does not guarantee that every coefficient is scientifically meaningful, but it is the algebraic condition for this unique least-squares solution.

Checking the Model’s Fit and Scope

The error assumptions are about the pattern of unexplained variation, not just whether the coefficient formula can be computed. Residuals are the observed values minus their fitted values. For a given gene, residuals plotted against fitted values or against sample order can reveal patterns such as changing spread or a trend across samples. Such patterns can indicate that the assumed equal variance or lack of sample structure is not credible. They are diagnostics rather than definitive tests: a small sample may not reveal a real problem, and a pattern may have several possible causes.

Sample dependence deserves particular attention. For example, repeated measurements from the same person may have correlated errors. If the model is fit as though those errors had covariance \(\sigma_g^2 I_n\), the coefficient estimate may still be calculable, but conventional uncertainty calculations based on the independent, equal-variance model can be misleading. A covariance matrix with nonzero off-diagonal entries represents a different error structure and calls for a model that accounts for that dependence.

The gene-wise model is related to, but more general in setup than, fitting a separate two-group test for every gene. A two-group comparison can be expressed as a linear model with a suitable design matrix; other predictors or multiple modeled terms can also be represented in that framework. This common representation helps keep the expression response, sample structure, and gene-specific coefficients explicit. It does not mean that every possible design is appropriate or that a linear model automatically handles repeated measures, unequal variances, or nonlinear relationships.

Common Pitfalls and Limitations

• Inferring independence from the covariance condition. A diagonal covariance matrix establishes zero covariance, not general independence. State independence separately, or invoke joint normality when using the result that uncorrelated normal components are independent.
• Confusing the dimensions of \(Y\) and \(X\). Rows refer to samples in both matrices. Columns of \(Y\) refer to genes, whereas columns of \(X\) refer to modeled sample-level predictors.
• Assuming the errors must be independent across genes. The gene-wise model specifies an error structure for each gene. It does not require the errors for different genes to be unrelated.
• Using a design matrix with redundant columns. If \(X\) is not full column rank, the formula using \((X^\mathsf{T}X)^{-1}\) is not available and the full coefficient vector is not uniquely estimable.
• Reading coefficients as causal effects. A coefficient describes the modeled association on the chosen expression scale, conditional on the predictors in the model. Causal interpretation requires additional design and identification assumptions.

The model is also limited by its chosen scale and linear predictor. A log transformation can make multiplicative differences additive, but it does not guarantee that the conditional mean is correctly represented or that the error variance is constant. Model assumptions must be considered for the particular samples and genes being analyzed, rather than treated as automatic consequences of using a linear model.