Regularized regression and held-out prediction

Practical 20: Genomics, systems biology and multi-omics

Learning objectives

Implement training-only centering and scaling; compare a supplied ridge model with a baseline; explain tuning and deployment-population limitations.

Biological question and data

Do the hypothetical molecular measurements improve prediction for a held-out group of animals?

The response uses arbitrary units and depends on a shared latent factor by construction. Animals 1–32 train the model and 33–48 form a supplied holdout. This split also changes farm composition; its representativeness should be discussed. Lambda is supplied, not chosen using test outcomes.

No biological download is required. Numerical inputs are hypothetical, not observed animal data. The runnable examples use base R. Count and abundance units are specified where used.

1. Define the split and preprocess using training only

Combine available log2 protein and metabolite features, then estimate their means and scales in training.

helper_dir <- if (dir.exists("R")) "R" else "../R"
stopifnot(file.exists(file.path(helper_dir,"teaching_data.R")))
source(file.path(helper_dir,"teaching_data.R"))
d <- teaching_data()
X <- cbind(t(d$protein),t(d$metabolite));y <- d$metadata$response
train <- 1:32;test <- 33:48
stopifnot(length(intersect(train,test))==0,identical(rownames(X),d$metadata$animal))
center <- colMeans(X[train,,drop=FALSE]);spread <- apply(X[train,,drop=FALSE],2,sd)
stopifnot(all(spread>0))
Z <- scale(X,center,spread)
stopifnot(max(abs(colMeans(Z[train,,drop=FALSE])))<1e-12)
table(d$metadata$farm[train]);table(d$metadata$farm[test])

farm1 farm2 farm3 
   16    16     0 

farm1 farm2 farm3 
    0     0    16 

The test animals do not influence the fitted transformation. The farm composition identifies the transfer task rather than being hidden.

2. Fit ridge with a separate unpenalized intercept

Use the minimization convention sum of squared errors plus lambda times squared coefficient norm.

lambda <- 4
ybar <- mean(y[train])
beta <- solve(crossprod(Z[train,,drop=FALSE])+lambda*diag(ncol(Z)),
 crossprod(Z[train,,drop=FALSE],y[train]-ybar))
pred <- as.vector(ybar+Z[test,,drop=FALSE]%*%beta)
baseline <- rep(ybar,length(test))
normal_equation_error <- max(abs((crossprod(Z[train,,drop=FALSE])+lambda*diag(ncol(Z)))%*%beta-
 crossprod(Z[train,,drop=FALSE],y[train]-ybar)))
stopifnot(normal_equation_error<1e-9)

The penalty shrinks standardized feature coefficients. Its magnitude depends on the stated objective scaling; software conventions can differ.

3. Evaluate on the same untouched animals

Compare error and inspect predictions without re-tuning on these scores.

rmse <- function(observed,predicted) sqrt(mean((observed-predicted)^2))
data.frame(model=c("training mean","supplied ridge"),
 test_RMSE=c(rmse(y[test],baseline),rmse(y[test],pred)))

           model test_RMSE
1  training mean  1.348718
2 supplied ridge  1.359551

plot(pred,y[test],pch=19,xlab="Held-out predicted response",ylab="Observed response (arbitrary units)")
abline(0,1,col="#d95f02")

One split shows the evaluation procedure. The generating model makes useful molecular association plausible, but the test score does not establish practical value in a real livestock population.

4. Tune only within the training animals

For demonstration, partition the original training group into fitting and internal validation subsets. Every candidate model re-estimates preprocessing from its fitting subset. The final test set remains outside this selection.

source(file.path(helper_dir,'prediction_helpers.R'))
# Inspect the helper interface before replacing the transparent calculation above.
args(fit_teaching_ridge)

function (X, y, lambda = 4, weights = rep(1, ncol(X))) 
NULL

The complete internal-selection exercise below is written transparently to preserve the penalty convention.

inner_fit <- 1:24;inner_valid <- 25:32
candidate_lambda <- c(.5,4,20)
inner_center <- colMeans(X[inner_fit,,drop=FALSE])
inner_scale <- apply(X[inner_fit,,drop=FALSE],2,sd)
inner_Z <- scale(X,inner_center,inner_scale)
inner_mean <- mean(y[inner_fit])
validation_rmse <- sapply(candidate_lambda,function(penalty) {
 b <- solve(crossprod(inner_Z[inner_fit,,drop=FALSE])+penalty*diag(ncol(X)),
 crossprod(inner_Z[inner_fit,,drop=FALSE],y[inner_fit]-inner_mean))
 prediction <- as.vector(inner_mean+inner_Z[inner_valid,,drop=FALSE]%*%b)
 sqrt(mean((y[inner_valid]-prediction)^2))
})
data.frame(lambda=candidate_lambda,internal_RMSE=validation_rmse)

  lambda internal_RMSE
1    0.5      1.288845
2    4.0      1.182934
3   20.0      1.139183

5. Refit the selected procedure and evaluate

Choose the penalty from internal validation only, refit on all original training animals and apply the training transformation to the original test group.

selected_lambda <- candidate_lambda[which.min(validation_rmse)]
selected_model <- fit_teaching_ridge(X[train,,drop=FALSE],y[train],lambda=selected_lambda)
selected_prediction <- predict_teaching_ridge(selected_model,X[test,,drop=FALSE])
c(internally_selected_lambda=selected_lambda,
 selected_test_RMSE=rmse(y[test],selected_prediction))

internally_selected_lambda         selected_test_RMSE 
                 20.000000                   1.246502 

This score evaluates the supplied selection procedure on this realization. Do not use it to reconsider the chosen penalty.

6. Distinguish selection from final evaluation

After choosing a penalty internally, refit preprocessing and coefficients using all original training animals, then evaluate the selected procedure once on the test group. The small contiguous internal split is only a supplied illustration; its farm composition and uncertainty still matter. More candidate penalties create more selection opportunity, not more independent evaluation evidence.

Interpretation challenge

Can the original held-out test score be used to override the internally selected penalty?

Suggested reasoning

Doing so uses the test set for selection. Its score would no longer assess an untouched selected procedure. Use a new evaluation set or an appropriate nested design.

Interpretation and limits

Feature selection, scaling, imputation and lambda tuning must occur inside training-only validation. If lambda is optimized on these test animals, they become a tuning set. For real prediction, use a grouping/time split matching deployment and quantify uncertainty across independent units.

Exercises

1. Why is the intercept not penalized here?
2. What happens if test outcomes select lambda?
3. Could an accurate predictor be a noncausal proxy?

Suggested answers

1. The model centers the response separately, leaving the mean outside the coefficient penalty.
2. The reported score is no longer an untouched test evaluation.
3. Yes; predictive usefulness does not establish intervention or mediation.

Continue learning

B12 Statistics, B15 Integration and B18 Synthesis.

Continue with integrated analysis. Return to the practical index or the unified map.