1 Cross validation
Il criterio è applicabile in tutti i problemi di regressione lineare e non lineare parametrica e non parametrica, dal momento che consiste nel cercare un’approssimazione non distorta dell’errore di previsione
E’ un’approssimazione dell’errore quadratico medio di previsione \[ MSEP=\mathrm{E}\left[\sum_{i=1}^{n}(\mu_i-\hat{\mu}_i)^2\right] \] Se in luogo \(\hat{\mu}_i\) utilizzassimo il classico valore stimato, \(\hat{y}_i\), non faremmo altro che riottenere la somma dei quadrati degli scarti \(Res(\hat{y})\)=\(\sum(y_i-\hat{y}_i)^2\), quantità minimizzata per stimare i parametri: Se per esempio aumentiamo il numero di parametri questa quantità diminuisce sicuramente, per cui non possiamo usarla per confrontare modelli.
*Questo è dovuto al fatto che la stima \(\hat{y}_i\) è stata ottenuta utilizzando \(y_i\): Insomma non pssiamo utilizzare la stessa quantità \(Res(\hat{y})\) sia per la stima dei parametri che per la valutazione della bontà del modello in confronto ad altri.
Si può confrontare ciascun valore \(y_i\) osservato con la stima \(\hat{y}_{(-i)}\) ottenuta dalle altre \(n-1\) osservazioni (cioè utilizzando la stima \(\hat{\boldsymbol{ \beta}}_{(-i)}\), ottenuto omettendo la \(i-\)esima osservazione \(y_i\) e la \(i-\)esima riga di variabili esplicative, \(\mathbf{x}^{\mathsf{T}}_i\) e calcolando poi \(\hat{y}_{(-i)}=\mathbf{x}^{\mathsf{T}}_i \hat{\boldsymbol{ \beta}}_{(-i)}\) ). \[ MSE(CV)=\frac{\sum_{i=1}^{n}(y_i-\hat{y}_{(-i)})^2}{n} \] Si ha una stima corretta dell’errore di previsione perchè nel calcolo di \(\hat{y}_{(-i)}\) non interviene \(y_i\).
Questa tecnica è molto utilizzata nell’approccio moderno alla costruzione di modelli con buone capacità previsive. Qualche volta viene indicata con la sigla LOOCV (leave one out cross validation), perchè le osservazioni vengono levate dal campione ad una ad una (e poi rimesse…); altre tecniche più complesse prevedono l’esclusione di una frazione fissata di valori, ma non me ne occuperò in questo corso.
1.1 Esempio base di errore di cross validation
data(antropometric)
n=20
ese=antropometric[1:n,c(7,8)]
plot(ese,main= "Esempio per tre valori di un campione di 20 osservazioni")
ind=c(8,9,10)
#ind=1:n
icol=2
lmtot=lm(PESOKG ~ALTEZZA,data=ese)
abline(lmtot,lty=3)
for (i in ind){
points(ese[i,],col=icol,pch=19)
yitot=predict(lmtot,newdata=ese[i,])
lmcv.i=lm(PESOKG ~ALTEZZA,data=ese[-i,])
ycv.i=predict(lmcv.i,newdata=ese[i,])
abline(lmcv.i,col=icol)
points(ese$ALTEZZA[i],ycv.i,col=icol)
segments(ese[i,1],ycv.i,ese[i,1],ese[i,2],lty=3,col=icol)
icol=icol+1
}
A titolo di esempio sugli stessi dati mostro il calcolo:
plot(ese,main= "Punti osservati, stime sul modello totale, stime di cross validation",col=1,pch=19)
abline(lmtot,lty=3)
ycv=array(0,n)
for (i in 1:n){
points(ese[i,],col=1,pch=19)
yitot=predict(lmtot,newdata=ese[i,])
lmcv.i=lm(PESOKG ~ALTEZZA,data=ese[-i,])
ycv[i]=predict(lmcv.i,newdata=ese[i,])
points(ese$ALTEZZA[i],yitot,col=1)
points(ese$ALTEZZA[i],ycv[i],col=2)
segments(ese[i,1],ycv[i],ese[i,1],ese[i,2],lty=3,col=2)
}
E’ interessante notare che, come è logico, tutti i residui di cross-validation sono maggiori di quelli classici calcolati attraverso la retta di regressione generale. Le somme dei quadrati dei residui sono:
Devianza residua classica: \(\sum_{i=1}^{n}(y_i-\hat{y}_{i})^2=615.04\)
Devianza residua di CV: \(\sum_{i=1}^{n}(y_i-\hat{y}_{(-i)})^2= 722.06\)
L’errore di Cross Validation introdotto può sembrare molto oneroso dal punto di vista computazionale, ma in realtà nel caso del modello lineare (e in molte altre situazioni ad esso riconducibili), non è necessario per calcolare gli scarti stimare \(n\) modelli lineari ciascuno con \(n-1\) dati.
1.2 Formula semplificata per il calcolo
Si può vedere, in questo e in altri problemi simili, che la somma dei quadrati di cross validation può essere espressa senza bisogno di fare i calcoli espliciti e ricorrendo a quantità disponibili dopo la stima del modello originario: \[ n MSE(CV)=\sum_{i=1}^{n}\left(\frac{y_i-\hat{y}_i}{1-h_{ii}}\right)^2 \] essendo \(h_{ii}\) l’\(i-\)esimo elemento diagonale della hat matrix \(\mathbf{H}\)
Invece di ricorrere alle proprietà delle matrici orlate ragioniamo in questo modo, ricordando le proprietà della hat matrix \(\mathbf{H}=\mathbf{X}\left(\mathbf{X}^{\mathsf{T}}\mathbf{X}\right)^{-1}\mathbf{X}^{\mathsf{T}}\):
ossia \(\hat{\mathbf{y}}=\mathbf{H}\mathbf{y}\) e
\(\hat{y}_i=\sum_{j=1}^{n}h_{ij} y_i=h_{ii} y_i+ \sum_{i \ne j} h_{ij} y_j\).
Ora consideriamo la stima dei parametri nel modello in cui manca la riga \(i-\)esima, ossia \(\hat{\boldsymbol{ \beta}}_{(-i)}\). Costruiamo ora un modello lineare in cui la riga \(i-\)esima (originariamente costituita da \((y_i,\mathbf{x}^{\mathsf{T}}_i )\)), viene rimpiazzata dalla riga \((\hat{y}_{(-i)},\mathbf{x}^{\mathsf{T}}_i\). In questo nuovo modello lo stima dei minimi quadrati è sempre \(\hat{\boldsymbol{ \beta}}_{(-i)}\) (perchè il punto aggiunto \((\hat{y}_{(-i)},\mathbf{x}^{\mathsf{T}}_i\) fornisce un residuo empirico nullo, in quanto \(\hat{y}_{(-i)}=\mathbf{x}^{\mathsf{T}}_i \hat{\boldsymbol{ \beta}}_{(-i)}\) ed è un punto che giace già sull’iperpiano adattato agli \(n-1\) dati).
Nell’esempio visto prima sarebbe il punto rosso che giace sulla retta rossa
Questo nuovo modello ha la stessa matrice del disegno\(\mathbf{X}\) del modello originale, e in particolare per la riga \(i-\)esima, possiamo riapplicare la proprietà degli elementi \(h_{ij}\) della hat matrix \(\mathbf{H}=\mathbf{X}\left(\mathbf{X}^{\mathsf{T}}\mathbf{X}\right)^{-1}\mathbf{X}^{\mathsf{T}}\), considerando che in questo modello artificiale l’ \(i-\)esimo valore stimato (\(\hat{y}_{(-i)}\)) è per costruzione uguale al valore osservato!
Quindi si ha nell’\(i-\)esimo modello costruito, riapplicando la formula vista per il modello completo:
\(\hat{y}_{(-i)}=h_{ii} \hat{y}_{(-i)}+ \sum_{i \ne j} h_{ij} y_j\), da cui:
\(\sum_{i \ne j} h_{ij} y_j=(1-h_{ii})\hat{y}_{(-i)}\) Ora sostituiamo questa espressione alla stessa sommatoria per il valore stimato nel modello completo:
\(\hat{y}_i=h_{ii} y_i+ \sum_{i \ne j} h_{ij} y_j =h_{ii} y_i +(1-h_{ii})\hat{y}_{(-i)}\)
da qui con pochi passaggi ricaviamo una nuova espressione del residuo empirico:
\[y_i-\hat{y}_i=y_i-h_{ii} y_i -(1-h_{ii})\hat{y}_{(-i)} = (1-h_{ii}) (y_i-\hat{y}_{(-i)})\]
e quindi il risultato da dimostrare:
\[(y_i-\hat{y}_{(-i)})=\frac{y_i-\hat{y}_i}{(1-h_{ii})} \]
Questa relazione, utilissima per il calcolo, è molto importante anche dal punto di vista statistico per far vedere che il residuo di cross-validation è sempre superiore al residuo ordinario, tanto più quanto maggiore \(h_{ii}\) ## esempi
1.3 Altri criteri di cross-validation
In alcuni contesti si usa anche un errore di cross-validation generalizzato: \[ GCV(\boldsymbol{ \alpha})=\frac{\sum_{i=1}^{n}\left(y_i-\hat{y}_i\right)^2} {\left(1-\frac{tr(\mathbf{H})}{n}\right)^2} \]
Per problemi non lineari o comunque più complessi di quelli base trattati in questo corso, si utilizzano tecniche più complesse dal punto di vista computazionale, ma molto semplici concettualmente.
Suddividiamo i dati in due insiemi:
- uno lo utilizziamo per stimare i parametri
- l’altro per valutare la bontà dell’adattamento Usualmente il procedimento viene ripetuto più volte (k-folds cross validation)
Eventuale esempio in laboratorio
Altro esempio completamente svolto su “MASTER2020crossval1.R”