Consideriamo tre tipi di norme (anche se in realtà la prima di queste non è una norma) \(l_0\), \(l_1\), \(l_2\) definite rispettivamente come:
\[\|\boldsymbol{\beta} \|_0 = \sum_{j=1}^{p}1\{\beta_j\neq 0\}, \quad \|\boldsymbol{\beta} \|_1 = \sum_{j=1}^{p} |\beta_j |, \quad \|\boldsymbol{\beta} \|_2 =\Big( \sum_{j=1}^{p} \beta_j^2 \Big)^{1/2}\]
L’utilizzo di ciascuna di queste norme si riflette in un problema di penalizzazione (se aggiunte come funzione di penalty, premoltiplicate per \(0 \leq \lambda\)) o in un problema di minimo vincolato secondo il seguente schema:
\(l_0\) -> Best Subset Selection
Penalized form
\[ \min_{\boldsymbol{\beta} \in \mathbb{R}^p} \| \mathbf{y} - \mathbf{X}\boldsymbol{\beta} \|^2_2 + \lambda\| \boldsymbol{\beta}\|_0 \]
Constrained form
\[ \min_{\boldsymbol{\beta} \in \mathbb{R}^p} \| \mathbf{y} - \mathbf{X}\boldsymbol{\beta} \|^2_2 \mathrm{\,\,subject\,\,to\,\,}\|\boldsymbol{\beta}\|_0 \leq k \]
Solo in questo caso le due forme non sono equivalenti, i.e. la soluzione di un problema in forma vincolata non è necessariamente uguale alla soluzione del problema in forma penalizzata. Il problema di best subset selection è NP-Hard, la tipologia di problemi non convessi più difficile da risolvere.
Per le norme di tipo 1 e di tipo 2, invece, le due soluzioni sono equivalenti: si tratta di problemi convessi. In particolare il problema ridge regression è strettamente convesso ed ammette forma chiusa, mentre il lasso non è strettamente convesso e per questo è non lineare in y e non ha una forma chiusa.
\(k\) è un parametro numerico intero, mentre per i problemi convessi \(t\) è uno scalare.
\(l_1\) -> LASSO
Penalized form
\[ \min_{\boldsymbol{\beta} \in \mathbb{R}^p} \| \mathbf{y} - \mathbf{X}\boldsymbol{\beta} \|^2_2 + \lambda\| \boldsymbol{\beta}\|_1 \]
Constrained form
\[ \min_{\boldsymbol{\beta} \in \mathbb{R}^p} \| \mathbf{y} - \mathbf{X}\boldsymbol{\beta} \|^2_2 \mathrm{\,\,subject\,\,to\,\,}\|\boldsymbol{\beta}\|_1 \leq t \]
\(l_2\) -> RIDGE
Penalized form
\[ \min_{\boldsymbol{\beta} \in \mathbb{R}^p} \| \mathbf{y} - \mathbf{X}\boldsymbol{\beta} \|^2_2 + \lambda\| \boldsymbol{\beta}\|_2 \]
Constrained form
\[ \min_{\boldsymbol{\beta} \in \mathbb{R}^p} \| \mathbf{y} - \mathbf{X}\boldsymbol{\beta} \|^2_2 \mathrm{\,\,subject\,\,to\,\,}\|\boldsymbol{\beta}\|_2 \leq k \]
Il concetto di sparsità consiste nel fatto che la maggior parte delle variabili esplicative presenti in un dataframe sono in realtà rumore, quindi vorremmo che i loro coefficienti \(\beta_j\) siano nulli.
Lo stimatore Ridge non è un metodo di stima sparso, mentre gli altri due stimatori trattati in questa lezione (Lasso e BSS) sì. L’intuizione grafica stà nel tipo di norma utilizzata: con la norma \(l_1\) lo spazio bidimensionale di una coppia di coefficienti vincolati ad una distanza \(t\) è un rombo, quindi è più probabile avere soluzioni d’angolo in cui uno dei due coefficienti vale 0.
Per scegliere il parametro di tuning \(\lambda\) si ricorre di solito alla Cross-Validazione: si divide il training set in K folds e si usa a turno un fold diverso come insieme di validazione del modello addestrato sugli altri (K-1) folds; si sceglie quindi il parametro \(\lambda\) che minimizza l’errore di Cross-Validazione.
Un metodo alternativo e più parsimonioso è quello del one-standard-error rule che consiste nel scegliere il modello più semplice il cui errore di Cross-Validazione sia superiore al massimo di una quantità pari ad 1 standard error.
# data
library(ISLR)## Warning: package 'ISLR' was built under R version 3.6.3data(Hitters)
Hitters = Hitters[complete.cases(Hitters),]
Hitters[,"League"]=(Hitters[,"League"]=="A")*1
Hitters[,"Division"]=(Hitters[,"Division"]=="E")*1
Hitters[,"NewLeague"]=(Hitters[,"NewLeague"]=="A")*1
set.seed(1)
train.id = sample(nrow(Hitters),nrow(Hitters)/2)
train = Hitters[train.id,]
names(train)[19] = "y"
# train
X = as.matrix(train[,-19])
y = train$y
n = nrow(X)
p = ncol(X)
# test
test = Hitters[-train.id,]
names(test)[19] = "y"
X.star = as.matrix(test[,-19])
y.star = test$y
m = nrow(X.star)
# package
require(glmnet)## Loading required package: glmnet## Warning: package 'glmnet' was built under R version 3.6.3## Loading required package: Matrix## Loaded glmnet 3.0-2K = 10
set.seed(123)
# RIDGE
ridge.cv <- cv.glmnet(X,y,alpha=0, nfolds = K, grouped=FALSE)
hatlambda <- ridge.cv$lambda.min #min = minimizza error CV
yhat.ridge <- predict(ridge.cv, s=hatlambda, newx=X.star, exact=TRUE)
# LASSO
cv.lasso <- cv.glmnet(X, y, alpha=1, nfolds = K)
hatlambda <- cv.lasso$lambda.1se #1se = 1 standard error rule
yhat.lasso <- predict(cv.lasso, s=hatlambda, newx=X.star, exact=T)
# BEST SUBSET SELECTION
require(leaps)## Loading required package: leaps## Warning: package 'leaps' was built under R version 3.6.3fit.bests <- regsubsets(y~.,train, nvmax=p)
summary.bests <- summary(fit.bests)
predict.regsubsets <- function(object, newdata, id, ...){
form=as.formula(object$call[[2]])
mat=model.matrix(form, newdata)
coefi =coef(object, id=id)
xvars =names(coefi)
mat[,xvars]%*%coefi
}
yhat.bestBIC <- predict.regsubsets(fit.bests, newdata=test, id=which.min(summary.bests$bic))