Residual Variance and Residual Degrees of Freedom

Why Residual Variance Depends on the Design

A gene-wise linear model separates systematic expression differences, described by the design, from remaining sample-to-sample variation. The residual variance estimate summarizes that remaining variation for a particular gene. It also depends on how much information the fit used: estimating more independent model directions leaves fewer residual degrees of freedom.

This distinction matters in expression analysis because each gene has its own observed residual variation, while the sample design is shared across genes. The fitted values and residual sum of squares vary by gene; if all genes use the same design and observations, the residual degrees of freedom are the same. The earlier material on contrasts describes questions about model coefficients; here the focus is how the fit supplies a residual variance estimate and how its degrees of freedom are determined.

The entries of \(y_g\), \(\hat{y}_g\), and \(e_g\) are expression values, often on a transformed scale such as log base 2. Accordingly, \(s_g^2\) is in squared expression units. Its square root \(s_g\) is in expression units and describes the estimated residual standard deviation. The residual variance is not the variance of the fitted values or the total variation in the observations: it is the variation left after fitting the specified design.

Model Conditions and the Least-Squares Residuals

For a gene-wise model with a common residual variance, assume the errors have mean zero and covariance matrix \(\sigma_g^2 I\). This means the observations have equal error variance and uncorrelated errors under the model. For exact normal-theory results, additionally assume \(\epsilon_g\sim N(0,\sigma_g^2 I)\). The normality assumption is not needed to define the least-squares fit or calculate its residual sum of squares.

When \(X\) has full column rank, the least-squares coefficient estimate is \(\hat{\beta}_g=(X'X)^{-1}X'y_g\), and the fitted values are \(\hat{y}_g=X\hat{\beta}_g\). If columns of \(X\) are redundant, this coefficient formula using an ordinary inverse does not apply. Least-squares fitted values are still defined as the projection of \(y_g\) onto the column space of \(X\), and the rank \(r\) counts the independent directions in that space.

The residuals are orthogonal to the fitted model space: \(X'e_g=0\). Thus, when \(X\) is full rank, fitting its \(q\) independent columns imposes \(q\) independent constraints on the \(n\) residuals. More generally there are \(r\) independent constraints, even if the design contains more than \(r\) columns. The residuals therefore lie in a space of dimension \(n-r\).

Deriving the Residual Variance Estimate

Let \(P_X\) denote the projection matrix onto the column space of \(X\), and let \(M_X=I-P_X\) denote the projection onto the residual space. These matrices are symmetric and idempotent, and \(P_X\) has rank \(r\). The residual vector is \(e_g=M_Xy_g\). Since \(M_XX=0\), the model component vanishes from the residuals:

Under the equal-variance, uncorrelated-error assumption, the expected residual sum of squares is obtained by taking the trace of the residual projection:

The final equality follows because \(M_X\) has \(n-r\) eigenvalues equal to 1 and \(r\) equal to 0, so its trace is \(n-r\). Dividing by this expectation factor gives \(E(s_g^2)=\sigma_g^2\). Thus \(s_g^2=\mathrm{RSS}_g/(n-r)\) is unbiased under the stated covariance assumptions. Using \(n\) or \(n-q\) in the denominator instead would generally not produce this unbiased estimate.

Under the additional normality assumption, rotate the errors into an orthonormal basis whose first \(r\) directions span the model space and whose remaining \(n-r\) directions span the residual space. In those residual directions, the standardized errors are independent standard normal random variables, and the residual sum of squares is their squared sum multiplied by \(\sigma_g^2\). Therefore:

This exact distribution relies on normal, independent errors with a common variance. Without normality, the formula for \(s_g^2\) still defines the least-squares residual variance estimate, and it remains unbiased under the stated covariance assumption, but the chi-squared distribution need not hold.

How Design Rank Determines Residual Information

The rank of \(X\) is the number of linearly independent columns, or equivalently the dimension of the model space. With a full-rank design, this equals the number of columns \(q\). With redundant columns, rank is smaller than \(q\); the fitted values remain identifiable as a projection, but some individual coefficients cannot be uniquely determined. In either case, the residual degrees of freedom are \(n-r\).

For a group-means design with \(k\) represented groups and an intercept-style group parameterization, the model space has rank \(k\). If there are \(n\) observations in total, residual degrees of freedom are \(n-k\): within-group deviations remain after estimating the group means. This count depends on the observations and the independent model directions, not on the number of genes.

When There Is Too Little Residual Information

If \(n=r\), the design is saturated: its model space fills the observation space. Least squares can then fit every possible response vector exactly, leaving \(\mathrm{RSS}_g=0\) and \(df_{\mathrm{res}}=0\). The quotient \(\mathrm{RSS}_g/df_{\mathrm{res}}\) is undefined, so the model provides no residual variance estimate. A zero residual sum of squares in this setting is not evidence that the measurement process has no variability; it is a consequence of having no residual directions.

Even when \(df_{\mathrm{res}}>0\), a very small value means the variance estimate is based on little residual information and can be unstable. Under the normal model, the chi-squared distribution above makes this variability explicit: with few degrees of freedom, the residual sum of squares can fluctuate substantially around its expectation. The number of genes does not add residual degrees of freedom to a gene-specific fit; each gene’s \(s_g^2\) is calculated from that gene’s sample residuals.

Rank also reflects the actual design and observations used. If an intended group has no samples, its effect cannot be estimated from that dataset as a distinct direction. If samples are removed, \(n\) may decrease; if the retained design loses independent columns, \(r\) may also decrease. Recalculate the rank and the residual degrees of freedom for the data actually fitted rather than assuming that the planned column count still applies.

Diagnostics, Alternatives, and Common Mistakes

Useful checks begin with the fit itself. Verify that the reported residual degrees of freedom equal \(n-\operatorname{rank}(X)\), inspect the design for duplicate or dependent columns, and confirm that the residuals are orthogonal to the model columns up to numerical precision. Residual plots can reveal structure left unexplained by the model, such as a trend across a predictor or markedly different spread across groups. Such patterns call into question whether a common-variance linear error model is adequate. A residual plot cannot establish normality or independence by itself, and technical or biological dependence must be considered from the study design.

The ordinary gene-wise estimate uses only that gene’s residual sum of squares and divides by its residual degrees of freedom. This is directly interpretable but can be unstable when the residual degrees of freedom are small. Other approaches can model unequal variances, dependence, or share information across genes, but they require additional assumptions and are not the same calculation as the ordinary residual variance estimate. In particular, borrowing information across genes does not change the rank-based residual degrees of freedom of this least-squares fit.

• Using the number of columns instead of rank. The wrong reasoning is that every written column costs one degree of freedom. The correction is \(df_{\mathrm{res}}=n-r\); redundant columns do not add independent fitted directions.
• Dividing by \(n\). This treats the fitted model as if it used no information. Under the stated assumptions, divide by \(n-r\) to obtain the unbiased residual variance estimate.
• Calling a saturated fit’s zero RSS a zero variance estimate. With \(n=r\), the residual variance formula has a zero denominator and is undefined, regardless of the fitted RSS.
• Interpreting residual variance as total expression variability. It measures variation around the fitted model, in squared expression units; variation explained by the model is not part of the residual sum of squares.