Unpenalized Covariates and Penalized Genomic Features

Why Separate Covariates From Genomic Features?

A genomic prediction model may need to account for variables that are specified in advance, such as age, sex, study site, or technical batch. These covariates serve a different role from the large set of genomic predictors: the analysis may require their effects to remain in the model regardless of whether they improve prediction enough to survive a feature-selection penalty. The estimation problem is therefore not simply “penalize every coefficient.” It is to estimate covariate effects without shrinkage while regularizing the genomic coefficients.

This distinction is especially important when the number of genomic features $p$ is much larger than the sample size $n$. Penalization can stabilize genomic prediction in this setting, but it should not silently change the treatment of required covariates. The model below makes that distinction explicit. It is a prediction model; including a covariate does not, by itself, give its coefficient a causal interpretation.

Definition: Let $y\in\mathbb{R}^n$ be a continuous outcome, $Z\in\mathbb{R}^{n\times q}$ the matrix of prespecified covariates, and $X\in\mathbb{R}^{n\times p}$ the genomic-feature matrix. The partially penalized linear model is $$ y=Z\gamma+X\beta+\varepsilon, $$ where $\gamma\in\mathbb{R}^q$ contains unpenalized covariate coefficients and $\beta\in\mathbb{R}^p$ contains penalized genomic coefficients. Include a column of ones in $Z$ when the model has an intercept. A common estimation objective is $$ (\hat{\gamma},\hat{\beta})= \arg\min_{\gamma,\beta} \left\{ \frac{1}{2n}\|y-Z\gamma-X\beta\|_2^2 +\lambda\left[\rho\|\beta\|_1+\frac{1-\rho}{2}\|\beta\|_2^2\right] \right\}, $$ where $\lambda\geq 0$ controls the penalty and $0\leq\rho\leq1$ sets its form. Only $\beta$ is penalized. When $\rho=1$, this is the lasso penalty; when $0<\rho<1$, it is an elastic-net penalty; when $\rho=0$, it is a ridge penalty.

The factor $1/(2n)$ is a scaling convention for the squared-error loss. It makes the penalty parameter’s numerical scale depend on that convention, so a value of $\lambda$ should always be understood together with the stated objective. For a fixed $\lambda$, the procedure estimates the covariate coefficients jointly with the genomic coefficients; “unpenalized” does not mean that the covariates are fitted first and then ignored.

Assumptions and Profiling Out the Covariates

For the squared-error objective to be a suitable model, the outcome should be reasonably represented by an additive linear prediction $Z\gamma+X\beta$, with the outcome measured on a scale for which squared-error loss is appropriate. For the usual interpretation of rows as independent observations, each row of $Z$, $X$, and $y$ must refer to the same individual, and feature processing must not use information from held-out outcomes. The optimization itself does not require normally distributed errors. Normality and other assumptions would be needed for particular inferential procedures, not merely to compute this penalized prediction fit.

Suppose $Z$ has full column rank. For any fixed $\beta$, the best-fitting unpenalized covariate coefficient is ordinary least squares applied to $y-X\beta$:

\hat{\gamma}(\beta)=(Z^\top Z)^{-1}Z^\top(y-X\beta).

Define the projection onto the columns of $Z$ as $P_Z=Z(Z^\top Z)^{-1}Z^\top$, and define the residual-maker matrix as $M_Z=I-P_Z$. Substituting $\hat{\gamma}(\beta)$ into the residual gives $M_Z(y-X\beta)$. Thus, for a fixed $\lambda$, estimating the penalized coefficients is equivalent to minimizing the penalized loss on the outcome and genomic features after projecting both off the covariate space:

\hat{\beta}=\arg\min_{\beta} \left\{ \frac{1}{2n}\|M_Zy-M_ZX\beta\|_2^2+ \lambda\left[\rho\|\beta\|_1+\frac{1-\rho}{2}\|\beta\|_2^2\right] \right\}.

This is a computational and conceptual result, not a recommendation to omit the covariates from the original model. The residualization is exact because $\gamma$ is unrestricted and the penalty depends only on $\beta$. After estimating $\beta$, recover $\gamma$ using $\hat{\gamma}(\hat{\beta})$. If $Z$ is rank deficient, use the projection onto its column space (for example, from a QR decomposition) rather than the inverse formula. The fitted values may still be identifiable even when individual coefficients in a redundant $Z$ are not.

Justification: For fixed $\beta$, the loss as a function of $\gamma$ is $\|y-X\beta-Z\gamma\|_2^2$. Its least-squares fitted component is the projection $P_Z(y-X\beta)$, leaving residual $M_Z(y-X\beta)$. Since the penalty is unchanged by $\gamma$, minimizing over $\gamma$ first and then over $\beta$ gives the same minimum as joint minimization. The equivalence would not hold in this form if the covariate coefficients were also penalized.

Worked Fits

Worked Example: A Covariate and Two Orthogonal Genomic Features

Consider four hypothetical individuals. Let $Z$ contain an intercept and a covariate $c=(-1,-1,1,1)^\top$. Let the two genomic columns be $x_1=(1,-1,1,-1)^\top$ and $x_2=(1,-1,-1,1)^\top$. Each vector is centered, and the three non-intercept columns are mutually orthogonal, with squared length 4. Suppose the outcome is $y=10\mathbf{1}+2c+1.5x_1+0.4x_2$. Its four values are $9.9, 6.1, 13.1, 10.9$. Fit a lasso model with $\lambda=0.5$.

Because each genomic column is orthogonal to $Z$ and to the other genomic column, its contribution to the loss separates. For a genomic column with $x_j^\top x_j/n=1$, its unpenalized score is $x_j^\top M_Zy/n$. Here the scores are 1.5 and 0.4. The one-coordinate lasso estimate is the soft-thresholded score: $\hat{\beta}_j=\operatorname{sign}(q_j)\max(|q_j|-\lambda,0)$, where $q_j=x_j^\top M_Zy/n$. Therefore,

\hat{\beta}_1=\max(1.5-0.5,0)=1.0,\qquad \hat{\beta}_2=\max(0.4-0.5,0)=0.

Recovering the unpenalized coefficients gives $\hat{\gamma}_0=10$ and $\hat{\gamma}_c=2$, because the covariate is orthogonal to both genomic columns. The fitted values are $10\mathbf{1}+2c+1.0x_1$. The residual is $0.5x_1+0.4x_2$; orthogonality gives its squared length divided by $n$ as $0.5^2+0.4^2=0.41$. The objective value is consequently $0.41/2+0.5(1.0)=0.705$. The fit retains the covariate effect without shrinkage, while the smaller genomic coefficient is set to zero by the lasso penalty.

Worked Example: Adjustment Changes the Genomic Estimate

Now let $c=(-1,-1,1,1)^\top$, $u=(1,-1,1,-1)^\top$, and $x=c+u=(0,-2,2,0)^\top$. Use $Z$ with an intercept and $c$, and suppose the hypothetical outcome is $y=5\mathbf{1}+2c+u=(4,4,8,8)^\top$. Fit a lasso model with one genomic feature $x$ and $\lambda=0.25$. The vectors $c$ and $u$ are centered and orthogonal, each with squared length 4.

The genomic feature is correlated with the covariate: its component along $c$ is removed by adjustment. Specifically, $M_Zx=u$ and $M_Zy=u$. Thus $x^\top M_Zx/n=1$ and $x^\top M_Zy/n=1$, so the residualized one-feature lasso estimate is

\hat{\beta}=\max(1-0.25,0)=0.75.

The covariate coefficient must then be recomputed conditional on this shrunken genomic estimate. Since $y=5\mathbf{1}+2c+u$ and $x=c+u$, the fitted genomic component is $0.75c+0.75u$. The remaining coefficient of $c$ is $2-0.75=1.25$, so $\hat{\gamma}_0=5$ and $\hat{\gamma}_c=1.25$. The residual is $0.25u$, giving loss $\frac12(0.25^2)=0.03125$; the penalty is $0.25(0.75)=0.1875$. The objective value is $0.21875$. Penalizing the genomic coefficient does not mean the covariate-adjusted estimate stays at its unpenalized value: $\gamma$ is unpenalized, but its fitted value is jointly adjusted as $\beta$ changes.

Worked Example: Standardizing a Genomic Feature Before Penalization

Suppose $n=4$, $Z$ contains an intercept and $c=(-1,-1,1,1)^\top$, and the raw genomic feature is $x_{\text{raw}}=2u$, where $u=(1,-1,1,-1)^\top$. Let the outcome be $y=3\mathbf{1}+2c+1.2u$. The covariate and genomic feature are orthogonal. The raw feature has $x_{\text{raw}}^\top x_{\text{raw}}/n=4$; standardize it to $x=x_{\text{raw}}/2=u$, which has $x^\top x/n=1$. Use lasso with $\lambda=0.4$, penalizing the coefficient of the standardized feature.

After removing $Z$, the outcome residual is $1.2u$, so $x^\top M_Zy/n=1.2$. The standardized-coordinate lasso estimate is $\hat{\beta}=1.2-0.4=0.8$. Its coefficient is in outcome units per one standardized feature unit. Since $x=x_{\text{raw}}/2$, the equivalent coefficient per raw feature unit is $\hat{\beta}_{\text{raw}}=0.8/2=0.4$. The fitted covariate coefficients remain $\hat{\gamma}_0=3$ and $\hat{\gamma}_c=2$. The residual is $0.4u$, so the loss is $\frac12(0.4^2)=0.08$; the penalty is $0.4(0.8)=0.32$, and the objective value is $0.40$.

This calculation shows why the scale used for genomic predictors must be documented. Penalization acts on the coefficient in the chosen feature scale. In contrast, changing the units of an unpenalized covariate changes its coefficient’s numerical value but does not impose shrinkage on its fitted contribution.

Practical Estimation, Diagnostics, and Alternatives

In practice, the genomic features are often centered and scaled using quantities calculated in the training data. The intercept is normally retained as an unpenalized column of $Z$, rather than standardized as though it were a genomic measurement. Prespecified covariates should be coded consistently with the scientific or operational specification. If the penalty parameter is selected by cross-validation, all centering, scaling, and other data-dependent preprocessing must be learned inside each training fold and then applied to that fold’s validation observations. Otherwise, information from validation data can leak into model fitting.

A useful diagnostic is to verify that the fitted residual is orthogonal to the columns of $Z$: $Z^\top(y-Z\hat{\gamma}-X\hat{\beta})$ should be approximately zero for an unweighted least-squares fit with unrestricted $\gamma$. Also check the feature scales, the coding and rank of $Z$, and whether the reported genomic coefficients correspond to standardized or raw predictors. If a covariate is accidentally included in the penalized matrix rather than $Z$, its coefficient can be shrunk or set to zero. If the genomic predictors are scaled inconsistently, the penalty applies unevenly for reasons unrelated to their predictive evidence.

The main alternatives answer different modeling needs. Penalizing all coefficients can be appropriate when no covariate is required to remain in the model, but it does not implement the present specification. Fitting the covariates first and then fitting genomic predictors to residuals is equivalent only when residualization is done consistently for both the outcome and genomic columns, the covariate effects remain unrestricted, and the same training data and objective are used. Omitting the covariates after residualizing only the outcome generally changes the problem. Ordinary least squares without penalization is another alternative, but when $p$ exceeds $n$, the genomic coefficients are not generally uniquely determined by the data.

This model does not guarantee that every prespecified covariate is scientifically appropriate, that a linear form is adequate, or that the selected genomic features are stable. Strong correlation among genomic features can make individual coefficients or selected features sensitive to the data, even when predictions are similar. Penalized coefficients are not ordinary least-squares estimates, and their inclusion in a predictive model does not automatically support unadjusted coefficient tests or causal claims. Those questions require separate methods and assumptions.

Common Pitfall

Common mistake: “The intercept and age are in the model, so they are adjusted for; I can apply the lasso to every column without changing that.” This is wrong: if the intercept and age columns are included in the penalized matrix, their coefficients are subject to shrinkage. Put them in $Z$, keep the penalty on $\beta$ only, and verify that any residualization or software penalty-factor settings preserve this distinction.

Check Your Understanding

Answer each question, then compare your reasoning with the brief solutions.

In the stated objective, which coefficients are penalized when $Z$ contains an intercept and age, and $X$ contains genomic features?
Why must both $y$ and the columns of $X$ be residualized against $Z$ when using the profiled objective?
Suppose a standardized, residualized genomic column has $x^\top x/n=1$ and $x^\top M_Zy/n=0.3$. What is its lasso estimate at $\lambda=0.4$?
For a cross-validation analysis, where should genomic-feature means and scales be calculated?

Solutions. (1) Only the genomic coefficients $\beta$ are penalized; the intercept and age coefficients in $\gamma$ are unpenalized. (2) Profiling out $\gamma$ leaves $M_Z(y-X\beta)=M_Zy-M_ZX\beta$; omitting either projection does not give the same profiled residual. (3) Soft-thresholding gives $\max(0.3-0.4,0)=0$, so the estimated coefficient is zero. (4) Calculate means and scales using each training fold only, then apply those training-fold quantities to its validation data.