Worked Ridge Fits for High-Dimensional Expression Data

From Expression Measurements to a Ridge Prediction

In expression prediction, a data set may contain measurements on thousands of genes but only a modest number of samples. A ridge fit gives a prediction rule by balancing agreement with the observed responses against the size of the coefficient vector. Its coefficients are generally all retained, but their magnitudes are pulled toward zero. A worked calculation makes clear what that shrinkage does to predictions and what it does not tell us about the biological role of an individual gene.

Let \(X\) be the \(n\times p\) design matrix of gene measurements and let \(y\) be the \(n\)-vector of observed responses. For this tutorial, use the squared-error ridge objective

\hat{\beta}_{\lambda} = \arg\min_{\beta\in\mathbb{R}^{p}} \left\{\|y-X\beta\|_2^2+\lambda\|\beta\|_2^2\right\}, \qquad \lambda>0.

Here \(\|\cdot\|_2^2\) is the sum of squared entries, \(\lambda\) is the penalty strength, and \(\hat{\beta}_{\lambda}\) is the fitted coefficient vector. The objective as written has no intercept; if an intercept is used in an analysis, it is ordinarily left unpenalized. The examples below use measurements on a common scale and omit an intercept to keep the calculations transparent.

Definition: A ridge prediction for a new sample with gene-measurement vector \(x_{\mathrm{new}}\) is \(x_{\mathrm{new}}^{T}\hat{\beta}_{\lambda}\) when the model has no intercept. Ridge shrinkage means that the penalty discourages large coefficient vectors; it does not generally set individual coefficients exactly to zero.

The optimization has a unique coefficient solution for every \(\lambda>0\), including when \(p>n\). Indeed, its Hessian with respect to \(\beta\) is \(2(X^{T}X+\lambda I_p)\). For any nonzero vector \(v\),

v^{T}(X^{T}X+\lambda I_p)v = \|Xv\|_2^2+\lambda\|v\|_2^2 > 0.

Thus \(X^{T}X+\lambda I_p\) is positive definite. Setting the gradient to zero gives the normal equations \((X^{T}X+\lambda I_p)\hat{\beta}_{\lambda}=X^{T}y\). When \(p\) is much larger than \(n\), it is often more convenient to calculate via the dual form.

Result: For \(\lambda>0\), the ridge solution can be written

\hat{\beta}_{\lambda} = X^{T}(XX^{T}+\lambda I_n)^{-1}y.

This form requires solving an \(n\times n\) system rather than a \(p\times p\) system. It is valid even if the columns of \(X\) are linearly dependent: \(XX^{T}+\lambda I_n\) is positive definite because \(u^{T}(XX^{T}+\lambda I_n)u=\|X^{T}u\|_2^2+\lambda\|u\|_2^2>0\) for every nonzero \(u\).

Worked Fit With More Genes Than Samples

Worked Example: Two Samples, Three Expression Predictors

Suppose a hypothetical expression data set has two samples and three gene predictors. The response is a standardized assay score, and the entries below are illustrative measurements on a common scale. The model has no intercept.

Sample	Gene A	Gene B	Gene C	Response
1	1	0	1	2
2	0	1	1	1

Thus \(X\) is \(2\times3\), so \(p>n\), and \(y=(2,1)^T\). First calculate the sample-space matrix:

XX^{T} = \begin{pmatrix} 2&1\\ 1&2 \end{pmatrix}.

At \(\lambda=1\), the dual system is \(\begin{pmatrix}3&1\\1&3\end{pmatrix}q=y\). Its determinant is \(3(3)-1=8\), so

q = \begin{pmatrix} 3&1\\ 1&3 \end{pmatrix}^{-1} \begin{pmatrix}2\\1\end{pmatrix} = \frac{1}{8} \begin{pmatrix} 3&-1\\ -1&3 \end{pmatrix} \begin{pmatrix}2\\1\end{pmatrix} = \begin{pmatrix}5/8\\1/8\end{pmatrix}.

Then \(\hat{\beta}_1=X^Tq=(5/8,1/8,6/8)^T=(0.625,0.125,0.750)^T\). At \(\lambda=4\), the system is \(\begin{pmatrix}6&1\\1&6\end{pmatrix}q=y\), whose determinant is \(36-1=35\). Therefore

q = \frac{1}{35} \begin{pmatrix} 6&-1\\ -1&6 \end{pmatrix} \begin{pmatrix}2\\1\end{pmatrix} = \begin{pmatrix}11/35\\4/35\end{pmatrix}, \qquad \hat{\beta}_4 = \begin{pmatrix}11/35\\4/35\\15/35\end{pmatrix} \approx \begin{pmatrix}0.3143\\0.1143\\0.4286\end{pmatrix}.

For a new hypothetical sample with expression vector \(x_{\mathrm{new}}=(1,1,0)^T\), the prediction at \(\lambda=1\) is \(1(0.625)+1(0.125)+0(0.750)=0.750\) response units. At \(\lambda=4\), it is \(11/35+4/35=15/35\approx0.4286\) response units. The stronger penalty reduces this prediction for this particular sample.

The coefficient magnitudes have decreased from \(\lambda=1\) to \(\lambda=4\) in this example, but shrinkage does not guarantee that every individual prediction moves toward zero. A prediction is a signed sum of the coefficients weighted by a sample's gene values; cancellations and the signs of those values matter.

Why the Coefficients Shrink

A useful way to understand the fit is to examine the directions of variation in \(X\). If \(X=UDV^T\) is a singular value decomposition, with positive singular values \(d_j\), the ridge coefficient along right singular vector \(v_j\) is multiplied relative to the least-squares component by \(d_j^2/(d_j^2+\lambda)\). This factor is between zero and one. Directions with smaller \(d_j^2\) are shrunk more strongly at the same \(\lambda\), while a direction with no data support receives no fitted component.

The same idea appears in the fitted responses. The matrix taking \(y\) to \(X\hat{\beta}_{\lambda}\) is \(XX^T(XX^T+\lambda I_n)^{-1}\). Along an eigenvector of \(XX^T\) with eigenvalue \(d^2\), its multiplier is \(d^2/(d^2+\lambda)\). This explains why increasing \(\lambda\) suppresses weakly supported directions and can stabilize predictions when genes are highly correlated.

Worked Example: Shrinkage Along Two Sample-Space Directions

For the data above, \(XX^T=\begin{pmatrix}2&1\\1&2\end{pmatrix}\) has eigenvectors \(u_+=(1,1)^T/\sqrt{2}\) and \(u_-=(1,-1)^T/\sqrt{2}\), with eigenvalues 3 and 1. The response decomposes as \(y=(3/2,3/2)^T+(1/2,-1/2)^T\), its components in those two directions.

At \(\lambda=1\), the fitted-response multipliers are \(3/(3+1)=3/4\) and \(1/(1+1)=1/2\). Hence the fitted response is \((3/4)(3/2,3/2)^T+(1/2)(1/2,-1/2)^T=(1.375,0.875)^T\). At \(\lambda=4\), the multipliers are \(3/7\) and \(1/5\), giving \((3/7)(3/2,3/2)^T+(1/5)(1/2,-1/2)^T\approx(0.7429,0.5429)^T\). These are also obtained directly as \(X\hat{\beta}_1\) and \(X\hat{\beta}_4\).

The eigenvalue-1 direction is more strongly attenuated than the eigenvalue-3 direction at either penalty. This is a statement about the fitted-data directions, not a claim that one particular gene is less important: directions generally combine several genes.

This calculation is a direct justification of the shrinkage factors. If \(XX^Tu=d^2u\), then

XX^T(XX^T+\lambda I_n)^{-1}u = \frac{d^2}{d^2+\lambda}u.

The denominator exceeds the numerator for \(\lambda>0\), so the multiplier is less than one. As \(\lambda\) increases, the multiplier decreases; as \(\lambda\) approaches zero for a positive eigenvalue, it approaches one. This does not imply that the ridge coefficient vector is a list of gene-specific importance measures. Correlated predictors can share or trade off coefficient magnitude, and the penalty makes no causal interpretation.

Units, Assumptions, and Practical Checks

For a statistical prediction interpretation, a common working model is \(y=X\beta^*+\varepsilon\), with errors having conditional mean zero given \(X\) and finite variance. Independent errors with common variance are a familiar additional assumption for simple risk calculations. The ridge minimizer itself is a deterministic calculation from \(X\), \(y\), and \(\lambda\); it does not require normally distributed errors. A model assumption should not be confused with evidence that the genes cause the response.

Gene predictors should be placed on an appropriate common scale before penalization, with centering and scaling parameters estimated from training data alone. Otherwise, a change in a predictor's measurement scale changes how strongly its coefficient is penalized. If cross-validation is used to choose \(\lambda\), each training fold must perform its own preprocessing; using held-out samples to determine scaling leaks information into validation. Inspect the design dimensions, response and predictor units, penalty convention, and the fitted predictions before interpreting a result.

Worked Example: Changing the Response Units

Suppose the response in the first example is reported in tenths of its original unit, so \(y'=10y=(20,10)^T\). At \(\lambda=1\), the dual solution is \(q'=10q=(50/8,10/8)^T\), and the coefficient vector is \(\hat{\beta}'_1=10\hat{\beta}_1=(6.25,1.25,7.50)^T\). The new-sample prediction is \(7.50\) in the new response units, equal to ten times the old prediction.

For any candidate coefficient vector \(\beta\), the corresponding rescaled coefficient vector is \(10\beta\). Its new objective is

\|10y-X(10\beta)\|_2^2+\lambda\|10\beta\|_2^2 = 100\left(\|y-X\beta\|_2^2+\lambda\|\beta\|_2^2\right).

Every candidate objective is multiplied by the same positive constant, so the minimizing numerical value of \(\lambda\) is unchanged under this squared-error formulation. A squared-error validation score is also multiplied by 100 at each candidate \(\lambda\), so its minimizer is unchanged when the same folds and candidate penalties are used. By contrast, rescaling predictors changes the penalty's effect on their coefficients and can change a useful numerical \(\lambda\).

In a real analysis, the penalty is commonly selected based on prediction performance using data not used to fit that candidate model. Ridge is often compared with ordinary least squares, which does not penalize coefficients and can be unstable or nonunique when \(p>n\), and with methods using a different penalty that can produce exact zero coefficients. Ridge's strength is stable, distributed shrinkage; its limitation is that it does not by itself provide a sparse list of genes. The choice depends on whether the goal is prediction, compact feature representation, or another scientific objective.

Common Mistakes and Failure Modes

Interpreting shrinkage as gene selection. A smaller ridge coefficient is not evidence that a gene is irrelevant, and a nonzero coefficient is not evidence of a gene's causal effect.
Comparing coefficient sizes across mismatched scales. A one-unit change in a gene measured on a different scale means something different. Scale predictors consistently and report the scale used.
Assuming every prediction must decrease as \(\lambda\) increases. The direction of a prediction depends on its gene values and coefficient signs. Compare predictions for the specified sample rather than inferring them from coefficient shrinkage alone.
Using a numerical penalty without its objective convention. An objective with a factor of \(1/n\) on the squared-error term uses a different numerical penalty scale. State the objective before reporting \(\lambda\).
Trusting an implausible fit without checking the calculation. Confirm that the dual system is \(XX^T+\lambda I_n\), that the penalty is positive, and that computed coefficients reproduce the stated predictions.
Preprocessing before splitting validation data. Estimate scaling parameters using each training portion only; otherwise the validation assessment can be optimistic.

Key takeaway: Ridge fits can be computed efficiently in the sample space when \(p>n\). Increasing \(\lambda\) attenuates supported directions of the fit, while the resulting prediction depends on how a new sample projects onto those directions. Interpret the model with its predictor scales, response units, and validation procedure in view.

Check Your Understanding

Try each question, then check the solutions below.

In a problem with \(n=40\) samples and \(p=2{,}000\) genes, which dimension of the dual linear system must be solved?
For an eigenvalue \(d^2=2\) and \(\lambda=2\), what is the fitted-response shrinkage factor along that direction?
If every response is multiplied by 5, how do the ridge coefficients and squared-error objective change at fixed \(\lambda\)? Does the minimizing \(\lambda\) change under the stated objective?
Why can a ridge coefficient's magnitude not be read as a causal effect of that gene on the response?

Solutions. (1) Solve a \(40\times40\) system, since the dual matrix is \(XX^T+\lambda I_{40}\). (2) The factor is \(d^2/(d^2+\lambda)=2/(2+2)=1/2\). (3) The coefficients multiply by 5, and both the squared-error and penalty terms multiply by 25 under the corresponding coefficient transformation; the minimizing numerical \(\lambda\) is unchanged. (4) Ridge describes a penalized predictive association conditional on the fitted design and model, not an intervention or a causal comparison; correlated genes can also share coefficient magnitude.