library(tidyr)
library(dplyr)
library(broom)
library(MASS)
data(meatspec, package="faraway")

Ridge Regression

In contesti in cui l’informazione è ridondante, anche i metodi consolidati e semplici come la regressione lineare ordinaria (OLS) non sono efficaci. Cerchiamo di capire perché. Per la regressione lineare ordinaria, esiste una formula chiusa che restituisce esattamente i parametri del modello:

\[ \hat{\boldsymbol{\beta}}_{OLS}=(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}, \]

Dove \(\mathbf{X}\) è una matrice \(n \times p\), con \(n\) che rappresenta la dimensione del campione e \(p\) il numero di parametri, mentre \(\mathbf{y}\) è un vettore colonna di dimensione \(n \times 1\) contenente le realizzazioni della variabile dipendente.

La formula dello stimatore Ridge coinvolge l’inversione della matrice \(\mathbf{X}^T\mathbf{X}\), che non sempre esiste quando si lavora con dati ridondanti. Alcuni linguaggi di programmazione, come R, sono in grado di gestire tali problemi. Ciò è possibile poiché, invece di invertire la matrice, utilizzano un metodo di decomposizione chiamato decomposizione QR. Tuttavia, il problema persiste, poiché le stime così ottenute sono inaffidabili e associate ad un’elevata variabilità.

Per provare che ciò che abbiamo appena affermato è vero, cerchiamo di stimare un modello lineare che includa tutte le covariate presenti nel dataset. Dalle nostre analisi esplorative precedenti (output non è stato mostrato), abbiamo notato una forte correlazione tra le colonne della matrice dei dati.

#Fit del del modello di regressione
modlm <- lm(fat ~ ., meatspec)
my_summary <- summary(modlm)
my_df <- broom::tidy(my_summary)
my_df

Gli errori standard risultano essere troppo elevati per poter considerare le stime ottenute come affidabili. Pertanto, diventa necessario stabilizzare l’inversione della matrice e, conseguentemente, migliorare la precisione delle stime ottenute.

In questo contesto, è opportuno introdurre la regressione Ridge, che presuppone che le covariate \(X_1,\dots, X_p\) siano state centrate utilizzando la media e scalate tramite la deviazione standard, e che lo stesso procedimento sia stato applicato alla variabile risposta \(Y\). Il vettore dei coefficienti stimato \(\hat{\boldsymbol{\beta}}_{Ridge}\) ha la seguente forma

\[ \hat{\boldsymbol{\beta}}_{Ridge}=(\mathbf{X}^T\mathbf{X}+\lambda \mathbf{I})^{-1}\mathbf{X}^T\mathbf{y} \]

dove \(\lambda\) è un coefficiente reale che sfortunatamente non può essere stimato, metre \(\mathbf{I}\) è una matrice identità con \(1\) sulle posizioni diagonali e \(0\) altrimenti.

Esistono varie interpretazioni possibili della regressione Ridge, tra cui la più intuitiva consiste nell’aggiungere una costante positiva \(\lambda\) agli elementi diagonali della matrice \(\mathbf{X}^T\mathbf{X}\) prima di invertirla. Questa operazione permette di risolvere il problema di non invertibilità della matrice, che si presenta quando si lavora con informazioni ridondanti. Tale problema rappresenta la prima motivazione alla base della proposta di questo metodo (Hoerl e Kennard, 1970).

# fit della regressione ridge
rgmod <- lm.ridge(fat ~ ., meatspec, lambda = seq(0, 5e-8, len=21))
matplot(rgmod$lambda, 
        coef(rgmod), type="l", 
        xlab=expression(lambda) ,
        ylab=expression(hat(beta)),
        col=4)

Purtroppo, non esiste un metodo di stima diretto per \(\lambda\), quindi l’analista deve esercitare la pazienza e provare diversi valori possibili. Una rappresentazione grafica comunemente utilizzata è il ridge trace plot. Questo grafico rappresenta i coefficienti stimati in funzione dei diversi valori di \(\lambda\) (asse delle ascisse). Per selezionare il valore di \(\lambda\), un metodo è la generalized cross-validation (GCV).

Per fortuna, il pacchetto MASS che implementa la funzione lm.ridge e calcola già internamente il valore di \(\lambda\) che minimizza la generalized crossvalidation (GCV).

which.min(rgmod$GCV)
## 3.25e-08 
##       14

la quale restituisce il valore di \(\lambda\) per cui la GCV è minima.

{matplot(rgmod$lambda, 
         coef(rgmod), 
         type="l", 
         xlab=expression(lambda) ,
         ylab=expression(hat(beta)),
         col=4)
  abline(v=1.75e-08)}

ma quali sono le stima associate al valore di \(\lambda=1.75e-08\) ?

my_coef<- coef(rgmod)[which.min(rgmod$GCV),]
my_coef_dt<- broom::tidy(my_coef)
print(my_coef_dt)
## # A tibble: 101 × 2
##    names         x
##    <chr>     <dbl>
##  1 ""         7.14
##  2 "V1"    9693.  
##  3 "V2"  -12590.  
##  4 "V3"    2907.  
##  5 "V4"    5368.  
##  6 "V5"   -9038.  
##  7 "V6"     263.  
##  8 "V7"   -5984.  
##  9 "V8"    6133.  
## 10 "V9"      54.1 
## # … with 91 more rows

le stime così ottenute saranno distorte ma con una variabilità minore.

Lasso Regression

La regressione Lasso è apparentemente simile alla Ridge, ma ci sono alcune differenze concettuali che vale la pena discutere. Qui vogliamo scegliere \(\hat{\boldsymbol{\beta}}\) tale da minimizzare \[ (\mathbf{y}-\mathbf{X}\boldsymbol{\beta})^T(\mathbf{y}-\mathbf{X}\boldsymbol{\beta})^T,\,\,\,\sum_{j=1}^p|\beta_j|\le t \]

dove \(\mathbf{y}\) è il vettore delle variabili di risposta, \(\mathbf{X}\) è la matrice delle covariate, \(\boldsymbol{\beta}\) è il vettore dei coefficienti di regressione e \(t\) è un valore prefissato che rappresenta la somma degli scostamenti assoluti dei coefficienti.

La regressione Lasso si differenzia dalla Ridge per l’elemento penalizzante introdotto da Tibshirani nel 1996. Questo elemento, che impone la restrizione \(\sum_{j=1}^p|\beta_j|\le t\), non ha una soluzione chiusa per il vettore dei parametri. Notiamo inoltre che man mano che \(t\) aumenta, vengono aggiunte più variabili al modello e i loro coefficienti diventano più grandi. Per \(t\) sufficientemente grande, la restrizione di \(\sum_{j=1}^p|\beta_j|\le t\) diventa ridondante e viene restituita la soluzione dei minimi quadrati.

Questo metodo è particolarmente utile quando ci aspettiamo che gli effetti siano sparsi, ovvero quando solo pochi predittori sono in grado di spiegare la risposta, mentre gli altri non hanno alcun effetto significativo. Il segreto della regressione Lasso risiede nella capacità di eliminare i predittori che non sono rilevanti per la regressione. Quando \(\hat{\beta}_j\) è pari a zero, il corrispondente predittore \(x_j\) viene completamente eliminato dal modello. D’altra parte, se stai usando la regressione Ridge, non c’è eliminazione delle variabili; questo metodo riduce solo l’effetto dei coefficienti stimati \(\hat{\beta}_j\).

Divertente notare come termine “Lasso” non è un acronimo, ma fa riferimento alla corda utilizzata dai cowboy per catturare il bestiame, in analogia con il processo di selezione delle variabili dal modello da parte dello statistico.

Il dataset in questo esempio è state contenuto nel pacchetto datasets in R. Esso contiene informazioni su vari aspetti degli Stati Uniti, come la popolazione, il reddito pro capite, la spesa pubblica, l’istruzione e molti altri.

In particolare, il dataset include 50 osservazioni, una per ogni stato degli USA, e 25 variabili. Tra queste variabili, troviamo il nome dello stato, la sua abbreviazione, la sua capitale, la sua regione geografica, il numero di rappresentanti del Congresso, la popolazione totale, la densità di popolazione, la percentuale di popolazione di razza bianca, la percentuale di popolazione di razza nera, il reddito pro capite, la spesa pubblica per studente, il numero di scuole primarie e secondarie, e molte altre.

require(lars)
data(state)
statedata <- data.frame(state.x77,row.names=state.abb)
lmod <- lars(as.matrix(statedata[,-4]),statedata$Life)
plot(lmod)

il grafico mostrato è molto simile a quello che abbiamo già visto nel caso della Ridge regression. Tuttavia, ci sono alcune differenze notevoli. In particolare, si può notare che alcuni dei coefficienti, che dipendono da \(t\), sono costretti ad assumere il valore zero. Questo significa che le variabili corrispondenti vengono escluse dal modello, il che rende la regressione Lasso un metodo di selezione delle variabili molto potente.

Notiamo come l’asse orizzontale è stato scalato per il valore massimo del vettore dei coefficienti ottenuti tramite il metodo dei minimi quadrati. Per valori piccoli di \(t\) solo il predittore \(4\) (tasso di omicidio), risulta diverso da zero. Man mano che \(t\) aumenta, entra un secondo predittore in gico (tasso dei diplomati). Dunque, man mano che \(t\) aumenta, vediamo anche la dimensione dei coefficienti stimati aumentare.

Ok, ma come scegliamo \(t\) ? La selezione del parametro \(t\) (ricordiamo che non può esere stimato con metodi convenzionali) di solito viene effettuata utilizzando la Cross-Validation (CV). Nel pacchetto lars questo è implementato attraverso un 10-folds CV. In breve, i dati vengono suddivisi casualmente in 10 gruppi poi successivamente i restanti nove vengono impiegati per prevedere uno dei gruppi rimanenti. Questo per ciascuno dei 10 gruppi in modo iterativo, per poi calcolare la performance complessiva della previsione.

set.seed(123)
cvlmod <- cv.lars(as.matrix(statedata[,-4]),statedata$Life)

Il valore suggerito dalla CV risulta essere \(t_{CV}=0,65\). Il quale, in riferimento al grafico precedente si riferisce ad un modello contenente solo quattro predittori. A questo punto estraiamo i coefficienti

cvlmod$index[which.min(cvlmod$cv)]
## [1] 0.6464646
my_coef_lasso<- predict(lmod,s=0.65657,type="coef",mode="fraction")$coef
my_coef_lasso_dt<- broom::tidy(my_coef_lasso)
print(my_coef_lasso_dt)
## # A tibble: 7 × 2
##   names               x
##   <chr>           <dbl>
## 1 Population  0.0000235
## 2 Income      0        
## 3 Illiteracy  0        
## 4 Murder     -0.240    
## 5 HS.Grad     0.0353   
## 6 Frost      -0.00169  
## 7 Area        0

In conclusione, il Lasso e la Ridge sono due tecniche di regolarizzazione simili ma con obiettivi differenti. Il Lasso è più indicato quando si presume che solo un piccolo sottoinsieme di variabili sia rilevante per spiegare le \(Y_i\), mentre la Ridge è più adatta quando si presume che tutte le variabili siano rilevanti, ma alcune siano maggiormente rilevanti di altre. Inoltre, il Lasso è particolarmente utile quando si lavora con dataset ad alta dimensionalità, poiché aiuta a selezionare solo le variabili più importanti e a ridurre il rischio di overfitting. D’altra parte, la Ridge può essere più adatta quando le variabili sono fortemente correlate tra loro. In sintesi, la scelta tra Lasso e Ridge dipende dalle specifiche caratteristiche del dataset e dalle domande di ricerca dell’analista.

Riferimenti

  • Faraway, J. (2015). Linear Models with R Second Edition CHAPMAN & HALL/CRC Texts in Statistical Science.