Regularized Linear Models (R) — glmnet¶
Ridge, Lasso and Elastic Net from the package that defined them¶
glmnet (Friedman, Hastie & Tibshirani, 2010) is the reference implementation of the coordinate-descent algorithm the Python notebook built from scratch — the package cross-check. Its one knob alpha selects the penalty: alpha=0 is Ridge (L2), alpha=1 is Lasso (L1), and in between is the Elastic Net. It fits the entire regularization path at once, standardizes predictors internally, and chooses the penalty by built-in cross-validation (cv.glmnet). We reproduce the Python notebook's arc: the high-dimensional sparse problem where regularization is essential, then California housing and credit default, closing on the Bayesian reading (Ridge = Gaussian prior, Lasso = Laplace prior — the frequentist face of the Variable-Selection arc). ROC-AUC is defined in the CART notebook (0.5 = chance, 1 = perfect).
options(repr.plot.width=13, repr.plot.height=4.6, warn=-1)
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages({library(glmnet); library(MASS)})
auc <- function(y,p){ r<-rank(p); n1<-sum(y==1); n0<-sum(y==0); (sum(r[y==1])-n1*(n1+1)/2)/(n1*n0) }
rmse <- function(a,b) sqrt(mean((a-b)^2))
cat("glmnet", as.character(packageVersion("glmnet")), "\n")
glmnet 5.0
1. When OLS fails — a high-dimensional sparse problem¶
The same setup as the Python notebook: $n=120$ observations, $p=200$ predictors, only 8 truly nonzero — the sparse, high-dimensional regime of the horseshoe notebook in the Variable-Selection arc. With $p>n$ the OLS fit is not unique (we use the minimum-norm pseudo-inverse) and predicts poorly; the question is whether a penalty recovers the 8 signals.
set.seed(0); n<-120; p<-200; k<-8
beta<-numeric(p); signals<-sample(p,k); beta[signals]<-sample(c(-2.5,-1.5,1.5,2.5),k,replace=TRUE)
X<-matrix(rnorm(n*p),n,p); y<-as.numeric(X%*%beta+rnorm(n,0,1.5))
Xte<-matrix(rnorm(600*p),600,p); yte<-as.numeric(Xte%*%beta+rnorm(600,0,1.5))
b_ols<-ginv(X)%*%(y-mean(y)); ols_rmse<-rmse(yte, Xte%*%b_ols+mean(y))
cat(sprintf("n=%d train, p=%d predictors, %d truly nonzero -> p>n, OLS not unique\n", n,p,k))
cat(sprintf("OLS (min-norm pseudo-inverse) test RMSE %.3f (noise sd 1.5 -- badly overfit)\n", ols_rmse))
n=120 train, p=200 predictors, 8 truly nonzero -> p>n, OLS not unique
OLS (min-norm pseudo-inverse) test RMSE 4.390 (noise sd 1.5 -- badly overfit)
2. The glmnet path and cross-validation¶
glmnet(X, y, alpha=1) fits the Lasso for a whole grid of penalties $\lambda$ in one call. The coefficient path (left) shows the 8 true signals leaving zero as $\lambda$ relaxes while the 192 noise predictors stay pinned near zero — Lasso doing variable selection. cv.glmnet (right) picks $\lambda$ by cross-validation: lambda.min (lowest CV error) and the sparser lambda.1se (the one-standard-error rule, the more parsimonious default).
fit<-glmnet(X,y,alpha=1)
cv<-cv.glmnet(X,y,alpha=1,nfolds=5)
par(mfrow=c(1,2), mar=c(4,4,3,1))
B<-as.matrix(fit$beta); ll<-log(fit$lambda)
matplot(ll, t(B), type="l", lty=1, col="grey80", xlab="log(lambda)", ylab="coefficient",
main="Lasso path (glmnet): 8 true signals in red")
matlines(ll, t(B[signals,]), lty=1, lwd=2, col="#c53030")
abline(v=log(cv$lambda.min), lty=2, col="#2f855a")
plot(cv); title("Cross-validation picks lambda", line=2.5)
par(mfrow=c(1,1))
bcv<-as.numeric(coef(cv,s="lambda.min"))[-1]; nnz<-sum(abs(bcv)>1e-8); rec<-sum(abs(bcv[signals])>1e-8)
cat(sprintf("CV-lasso (lambda.min) keeps %d/%d predictors, recovers %d/%d true signals; test RMSE %.3f (vs OLS %.3f)\n",
nnz,p,rec,k, rmse(yte, predict(cv,Xte,s="lambda.min")), ols_rmse))
CV-lasso (lambda.min) keeps 41/200 predictors, recovers 8/8 true signals; test RMSE 1.677 (vs OLS 4.390)
3. Ridge vs Lasso vs Elastic Net — out of sample¶
The three penalties on the sparse problem, each cross-validated. Ridge (alpha=0) shrinks every coefficient but keeps all 200 — no selection, so it cannot isolate the 8 signals and predicts poorly. Lasso (alpha=1) selects a sparse set and slashes the error. Elastic Net (alpha=0.5) sits between. When the truth is sparse, L1 selection beats L2 shrinkage.
fit_oos<-function(a){ m<-cv.glmnet(X,y,alpha=a,nfolds=5)
list(rmse=rmse(yte,predict(m,Xte,s="lambda.min")), nz=sum(abs(as.numeric(coef(m,s="lambda.min"))[-1])>1e-8)) }
set.seed(1); R<-list(Ridge=fit_oos(0), ElasticNet=fit_oos(0.5), Lasso=fit_oos(1))
nm<-c("OLS",names(R)); rm_<-c(ols_rmse,sapply(R,function(z)z$rmse)); nz<-c(p,sapply(R,function(z)z$nz))
par(mar=c(4,4,3,1)); bp<-barplot(rm_, names.arg=nm, col=c("#a0aec0","#2b6cb0","#2f855a","#c53030"),
ylab="test RMSE (lower=better)", main="Out-of-sample on the p>n sparse problem", ylim=c(0,max(rm_)*1.15))
text(bp, rm_+0.1, sprintf("%.2f\n%d vars", rm_, nz), cex=0.9)
cat(sprintf("Ridge stays dense (%d vars); Lasso selects %d and more than halves the OLS error.\n", nz[2], nz[4]))
Ridge stays dense (200 vars); Lasso selects 41 and more than halves the OLS error.
4. Real data I — California housing (regression)¶
California housing (20,640 block groups, predicting median value in $100k) has 8 predictors and $n\gg p$ — the opposite regime, where OLS is already stable and regularization can help only marginally. That is the lesson: regularization pays off in high dimensions, not when data is abundant relative to predictors. All the linear models also trail the previous subsection's tree ensembles (~0.49 RMSE) because a linear fit cannot capture the nonlinear income–price relationship.
h<-read.csv("cali_housing.csv"); hf<-setdiff(names(h),"MedHouseVal")
set.seed(0); i<-sample(nrow(h),0.7*nrow(h)); Xh<-as.matrix(h[i,hf]); yh<-h$MedHouseVal[i]
Xhte<-as.matrix(h[-i,hf]); yhte<-h$MedHouseVal[-i]
ols<-lm(yh~Xh); ph<-cbind(1,Xhte)%*%coef(ols)
cal<-c(OLS=rmse(yhte,ph)); pe<-NULL
for(nm_ in c("Ridge","Lasso","ElasticNet")){ a<-c(Ridge=0,Lasso=1,ElasticNet=0.5)[nm_]
m<-cv.glmnet(Xh,yh,alpha=a,nfolds=5); pr<-as.numeric(predict(m,Xhte,s="lambda.min"))
cal[nm_]<-rmse(yhte,pr); if(nm_=="ElasticNet") pe<-pr }
cat("California test RMSE ($100k):\n"); for(nm_ in names(cal)) cat(sprintf(" %-11s %.4f\n", nm_, cal[nm_]))
cat(sprintf("The penalties barely separate (spread %.3f) -- with n>>p regularization is near-neutral; the gap to trees\n", max(cal)-min(cal)))
cat("(~0.49) is NONLINEARITY. The predicted-vs-actual graph (regression 'proportions vs predictions') shows where it misses.\n")
options(repr.plot.width=6.4, repr.plot.height=5)
plot(pe, yhte, pch=".", col="#a0aec0", xlab="predicted value ($100k)", ylab="actual value ($100k)",
main=sprintf("California -- predicted vs actual (ElasticNet, RMSE %.3f)", cal["ElasticNet"]))
abline(0,1,lty=2)
qb<-quantile(pe,seq(0,1,length=11)); bb<-cut(pe,unique(qb),include.lowest=TRUE)
points(tapply(pe,bb,mean), tapply(yhte,bb,mean), type="b", pch=19, col="#c53030", lwd=2)
legend("topleft", c("perfect","decile means"), col=c("black","#c53030"), lty=c(2,1), pch=c(NA,19), lwd=2, bty="n")
California test RMSE ($100k):
OLS 0.7342 Ridge 0.7504 Lasso 0.7336 ElasticNet 0.7338
The penalties barely separate (spread 0.017) -- with n>>p regularization is near-neutral; the gap to trees
(~0.49) is NONLINEARITY. The predicted-vs-actual graph (regression 'proportions vs predictions') shows where it misses.
5. Real data II — regularized logistic on credit default¶
glmnet(..., family="binomial") fits penalized logistic regression. On the Taiwan credit-default data (30,000 clients, 23 predictors, ~22% default; goal: classify defaulters) this is again $n\gg p$, so at the cross-validated penalty L1 prunes nothing. Its value is the sparsity-vs-accuracy tradeoff along the path: the whole Lasso path is fit at once, and we read test AUC against the number of features each $\lambda$ keeps (fit$df) to find the compact scorecard that matches the full model.
d<-read.csv("credit_default.csv"); feat<-setdiff(names(d),"default")
set.seed(0); j<-sample(nrow(d),0.7*nrow(d)); Xc<-as.matrix(d[j,feat]); yc<-d$default[j]
Xct<-as.matrix(d[-j,feat]); yct<-d$default[-j]
plain<-glm(default~., data=d[j,], family=binomial); au_plain<-auc(yct, predict(plain,d[-j,],type="response"))
fitL<-glmnet(Xc,yc,alpha=1,family="binomial")
P<-predict(fitL,Xct,type="response"); aucs<-apply(P,2,function(pp)auc(yct,pp)); df<-fitL$df
ok<-which(aucs>=max(aucs)-0.003); spi<-ok[which.min(df[ok])]
rel<-function(y,p){ q<-quantile(p,seq(0,1,length=11)); b<-cut(p,unique(q),include.lowest=TRUE); list(pp=tapply(p,b,mean), ot=tapply(y,b,mean)) }
brier<-function(y,p) mean((p-y)^2)
options(repr.plot.width=16, repr.plot.height=4.6); par(mfrow=c(1,3), mar=c(4,4,3,1))
plot(df, aucs, type="b", pch=19, col="#c53030", xlab="number of features kept by L1", ylab="test AUC",
main="Lasso-logistic: sparsity vs accuracy")
abline(h=au_plain, lty=2, col="grey50"); points(df[spi], aucs[spi], pch=19, col="#2f855a", cex=1.8)
legend("bottomright", c(sprintf("unpenalized logistic (%.3f)",au_plain), sprintf("sparse pick: %d feats, AUC %.3f",df[spi],aucs[spi])),
col=c("grey50","#2f855a"), pch=c(NA,19), lty=c(2,NA), bty="n", cex=0.8)
# reliability: observed vs predicted default proportion (proportions vs predictions)
pf<-predict(plain,d[-j,],type="response"); ps<-P[,spi]; rp<-rel(yct,pf); rs<-rel(yct,ps); mm<-max(pf,ps)
plot(c(0,mm),c(0,mm),type="l",lty=2,xlab="predicted P(default) [decile bins]",ylab="observed default proportion",
main="Credit -- reliability (proportions vs predictions)")
points(rp$pp,rp$ot,type="b",pch=19,col="grey50",lwd=2); points(rs$pp,rs$ot,type="b",pch=19,col="#c53030",lwd=2)
legend("topleft", c(sprintf("full logistic (Brier %.3f)",brier(yct,pf)), sprintf("%d-feat L1 (Brier %.3f)",df[spi],brier(yct,ps))),
col=c("grey50","#c53030"), pch=19, lwd=2, bty="n", cex=0.8)
bsp<-as.numeric(fitL$beta[,spi]); ord<-order(abs(bsp),decreasing=TRUE); ord<-ord[abs(bsp[ord])>1e-8][1:min(12,df[spi])]
par(mar=c(4,7,3,1)); barplot(rev(bsp[ord]), names.arg=rev(feat[ord]), horiz=TRUE, las=1,
col=ifelse(rev(bsp[ord])>0,"#c53030","#2b6cb0"), main=sprintf("Sparse scorecard: %d features",df[spi]), xlab="coefficient")
par(mfrow=c(1,1))
cat(sprintf("unpenalized logistic AUC %.4f (tree-scoreboard baseline). L1 path: a %d-feature scorecard reaches AUC %.3f,\n",au_plain,df[spi],aucs[spi]))
cat("within 0.003 of the full model -- recent repayment status (PAY_*) carries almost all the signal.\n")
unpenalized logistic AUC 0.7120 (tree-scoreboard baseline). L1 path: a 16-feature scorecard reaches AUC 0.709,
within 0.003 of the full model -- recent repayment status (PAY_*) carries almost all the signal.
6. The Bayesian bridge, and summary¶
glmnet reproduces the from-scratch coordinate descent exactly: Ridge, Lasso and Elastic Net as one algorithm with a mixing knob alpha. On the $p>n$ sparse problem Lasso recovered the 8 signals and halved the OLS error while Ridge stayed uselessly dense; on abundant low-dimensional California and credit data ($n\gg p$) the penalties were near-neutral, and the L1 path still yielded a compact credit scorecard.
The penalties are priors in disguise: Ridge is the MAP under a Gaussian prior on the coefficients, Lasso the MAP under a Laplace prior — the penalty is exactly $-\log(\text{prior})$. The fully-Bayesian version of this shrinkage is the Variable-Selection arc: SSVS (spike-and-slab), the horseshoe for $p\gg n$ (the Bayesian relative of §1–3's sparse recovery, cross-checked here), and Bayesian model averaging — a posterior over which coefficients are nonzero, not a single penalized point estimate. The next notebook, SVM & kernel methods, adds nonlinearity and the second bridge: kernel ridge = Gaussian-process posterior mean (the BNP Gaussian-Process notebooks). This notebook is the package mirror of reglm_python.ipynb.