OLS variance

Se \(n \gg p\) allora lo stimatore OLS tenderà ad avere una varianza più bassa, cioé mi aspetto stime simili per i coefficienti delle variabili esplicative pur stimando il modello con dati diversi;
se \(n \approx p\) allora lo stimatore OLS tenderà ad andare in overfitting e restituire stime non affidabili una volta cambiato l’insieme di dati;
se \(n=p\) il modello overfitta completamente i dati, interpolandoli, e restituisce delle pessime previsioni;
infine se \(n<p\) non esiste una soluzione unica per le stime dei coefficienti \(\beta\) e non possiamo usare OLS.

Dobbiamo quindi in definitiva tenere sotto controllo la varianza dello stimatore OLS e, nell’ultimo caso, trovare una soluzione ricorrendo ad uno dei seguenti metodi:

1.Subset Selection
Identifichiamo le variabili che riteniamo essere in relazione con la variabile target ed addestriamo il modello solo sulle variabili così selezionate.

2.Dimension Reduction
Anziché selezionare variabili, usiamo tutte e \(p\) le variabili per proiettarle in un nuovo spazio di \(q < p\) variabili, dove ciascuna delle nuove variabili è una combinazione lineare di quelle originarie. Usiamo le \(q\) variabili per stimare il modello.

3. Regularization
Uso tutte e \(p\) le variabili, ma riformulo il problema OLS in uno di minimizzazione vincolata, che - relativamente al caso standard - porta a delle stime dei \(\beta\) shrinked verso lo 0.

RIDGE Regression

Caso di regressione con regolarizzatore, dove usiamo tutte e \(p\) le variabili ma imponiamo il vincolo che i \(\beta\) in norma non possono allontanarsi troppo dall’origine.
In termini formali il problema è:

\[\begin{aligned} & \min_{\boldsymbol{\beta} \in \mathbb{R}^p} \| \mathbf{y} - \mathbf{X}\boldsymbol{\beta} \|^2_2\\ & \mathrm{\,\,subject\,\,to\,\,}\|\boldsymbol{\beta}\|^2_2 \leq t \end{aligned}\]

che può riformulato in un problema di ottimizzazione libera semplicemente aggiungendo al vincolo il moltiplicatore di Kuhn-Tucker \(\lambda \in [0,\infty)\).

\[\min_{\boldsymbol{\beta} \in \mathbb{R}^p} \left\{ \sum_{i=1}^{n}(y_i - x_i^\mathsf{T}\boldsymbol{\beta} )^2 + \lambda\sum_{j=1}^{p}\beta_j^2 \right\} = \min_{\boldsymbol{\beta} \in \mathbb{R}^p} \left\{ \underbrace{ \| \mathbf{y} - \mathbf{X}\boldsymbol{\beta} \|^2_2 }_{\mathrm{RSS} }+ \underbrace{\lambda\| \boldsymbol{\beta}\|^2_2}_{\mathrm{penalty}} \right\}\]

Il minimo per la parte RSS si ottiene per \(\beta = \beta_{ols}\) mentre la funzione di penalty è minimizzata per \(\boldsymbol{\beta} = \boldsymbol{0_p}\), per cui la soluzione del problema sarà La soluzione del problema è la stima Ridge del vettore dei coefficienti: per calcolarla analiticamente basta considerare le derivate parziali rispetto a \(\beta\) e porle uguali a zero (First Order Conditions).

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

Questa soluzione sarà sempre unica per ogni design matrix \(X\), in quanto la matrice \((\mathbf{X}^\mathsf{T}\mathbf{X} + \lambda \mathbf{I}_p )\) è sempre invertibile per ogni \(\lambda > 0\).

All’aumentare del moltiplicatore, sarà più forte la penalità assegnata a \(\beta\) “grandi” in valore assoluto e quindi prevarrà l’effetto di shrinking verso lo 0.

Poiché il termine di penalty pesa allo stesso modo tutti i \(\beta\) bisogna prima standardizzare le variabili esplicative, i.e. bisogna centrarle (renderle a media 0) sottraendo a ciascuna la propria media campionaria e bisogna portare ad 1 la deviazione standard dividendo ciascuna per la propria deviazione standard campionaria.

Di solito anche la variabile risposta viene centrata e scalata e l’intercetta viene esclusa dalla funzione di penalty, in modo che la media dei valori di \(y\) coincide con la stima di \(\beta_0\).

Come scegliere il moltiplicatore \(\lambda\) ?
Nella pratica il parametro di penalty è ignoto e viene scelto con K fold CV o cross-validazione generalizzata.
Con il pacchetto glmnet si imposta il parametro \(\alpha=0\) per fare una regressione Ridge.
Con la funzione cv.glmnet() si esegue la cross-validazione e si specifica il numero di folds nell’argomento \(nfolds=K\). Si sceglie quindi il parametro ottimale chiamando la componente $lambda.min del modello addestrato.

Nel caso di LOOCV esiste una shortcut che permette di calcolare il valore ottimale di \(\lambda\) senza dover addestrare K-volte il modello, ma sfruttando la seguente relazione tra l’RSS che c’è nel K-esimo fold con l’RSS fatto sulla totalità dei dati:
\[ \frac{1}{n}\sum_{i=1}^{n}\Big[ y_i - x_i^\mathsf{T} \hat{\boldsymbol{\beta}}^{-i}(\lambda) \Big]^2 = \frac{1}{n}\sum_{i=1}^{n} \left[ \frac{y_i - x_i^\mathsf{T} \hat{\boldsymbol{\beta}}(\lambda)}{1-\{\mathbf{H}(\lambda)\}_{ii}} \right]^2\]

dove \(\mathbf{H}(\lambda)\) è la matrice di proiezione dei valori fittati, i.e. la matrice che post-moltiplicata per il vettore risposta fornisce il vettore delle previsioni.
Più formalmente è definita come \(\mathbf{H}(\lambda) = \mathbf{X}(\mathbf{X}^\mathsf{T}\mathbf{X} + \lambda \mathbf{I}_p )^{-1}\mathbf{X}^\mathsf{T}\)
Possiamo quindi calcolare la performance associata ad un particolare valore di \(\lambda\) a partire dalla matrice di proiezione e dal vettore \(y\), senza dover eseguire le \(K=n\) regressioni ridge (vantaggio computazionale notevole).

Bias-Variance Trade-off dello stimatore Ridge Il teorema n.2 di Theobald (1974) dimostra che, preso come criterio decisionale la somma \((Bias)^2 + Var\), esiste sempre un certo valore \(\lambda >0\) tale per cui la somma del Bias-Variance dello stimatore Ridge è inferiore a quella dello stimatore OLS.
In altre parole, esiste un valore del parametro per cui l’aumento di distorsione associato allo stimatore Ridge è inferiore alla riduzione di varianza, portando ad un livello complessivo di errore inferiore rispetto al caso standard OLS. Il problema è che tale valore del parametro dipende dalle quantità ignote \(\beta\) e \(\sigma\).

Passare dallo stimatore OLS allo stimatore Ridge

L’operatore lineare \(W = [\mathbf{I}_{p\times p} + \lambda(\mathbf{X}^\mathsf{T}\mathbf{X})^{-1}]^{-1}\) post-moltiplicato il vettore \(\beta_{ols}\) dà come risultato lo stimatore ridge \(\beta(\lambda)=H(\lambda)y\).

Dato valore di \(\lambda\) e noto lo stimatore OLS \(\beta_{OLS}\) si ricava quindi lo stimatore Ridge con la formula: \(\hat{\beta}(\lambda)=W \beta_{OLS}\).

library(ISLR)
data(Hitters)

# remove na
Hitters = na.omit(Hitters)
# transform into data frame the data matrix
Hitters = data.frame(model.matrix(lm(Salary ~ ., Hitters))[,-1], Salary=Hitters$Salary)
head(Hitters)
                  AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns CRBI CWalks LeagueN DivisionW PutOuts Assists Errors NewLeagueN Salary
-Alan Ashby         315   81     7   24  38    39    14   3449   835     69   321  414    375       1         1     632      43     10          1  475.0
-Alvin Davis        479  130    18   66  72    76     3   1624   457     63   224  266    263       0         1     880      82     14          0  480.0
-Andre Dawson       496  141    20   65  78    37    11   5628  1575    225   828  838    354       1         0     200      11      3          1  500.0
-Andres Galarraga   321   87    10   39  42    30     2    396   101     12    48   46     33       1         0     805      40      4          1   91.5
-Alfredo Griffin    594  169     4   74  51    35    11   4408  1133     19   501  336    194       0         1     282     421     25          0  750.0
-Al Newman          185   37     1   23   8    21     2    214    42      1    30    9     24       1         0      76     127      7          0   70.0
# length of test set
m = 100
# length of train set
n = nrow(Hitters)-m

set.seed(123)
test.id = sample(nrow(Hitters),m)
test = Hitters[test.id,-20]
ytest = Hitters[test.id,20]
train = Hitters[-test.id,]

lambda.vals = exp(seq(-10,10,by=0.1))


library(MASS)
fits = lm.ridge(Salary~.,data=train, lambda=lambda.vals)
preds = as.matrix(test) %*% t(coef(fits)[,-1]) + rep(1,m) %o% coef(fits)[,1]
resids = matrix(data=ytest, nrow=dim(preds)[1], ncol=dim(preds)[2], byrow=F)-preds
RMSEs = sqrt( apply(resids^2, 2, sum)/m )
plot(fits$lambda, RMSEs, type="l", xlab="lambda", ylab="RMSE", log="x")
points(fits$lambda[which.min(RMSEs)],min(RMSEs))

min = fits$lambda[which.min(RMSEs)]
coef(lm.ridge(Salary~.,data=train, lambda=min))
                      AtBat          Hits         HmRun          Runs           RBI         Walks         Years        CAtBat         CHits        CHmRun         CRuns          CRBI        CWalks       LeagueN 
 205.88399506   -1.89686063    4.93197606    3.43108687   -1.32839807    2.81862442    5.96348988   -7.31347964   -0.03074841   -0.02593659   -1.14173125    0.96682887    0.76012119   -0.80554313   51.35929558 
    DivisionW       PutOuts       Assists        Errors    NewLeagueN 
-151.70829282    0.40577805    0.43152667   -5.49071796  -23.96469332

Principal Component Regression

set.seed(123)

X = model.matrix(lm(Salary ~ ., Hitters))
n = nrow(X)
Z = scale(X[,-1])[,]*sqrt((n-1)/n)
# SVD
U = svd(Z)$u
D = diag(svd(Z)$d)
V = svd(Z)$v
# beta principal component
b = V%*%D%*%t(U)%*%Hitters$Salary

library(pls)
fit.pcr=pcr(Salary ~ ., data=Hitters, scale=TRUE, validation ="CV")
validationplot(fit.pcr, val.type="MSEP")

fit=pcr(Salary ~ ., data=Hitters, scale=TRUE, ncomp=6)
resids = predict(fit, newdata=test) - ytest

RMSE = sqrt(mean(resids^2))
RMSE >= min(RMSEs)
[1] TRUE
# true logical condition : the pcr perform worse than ridge with lambda selected by cv