1. Validation Set Approach

Usare gli stessi dati per addestrare e testare il modello porta a dei risultati troppo ottimisti, che difficilmente resterebbero tali di fronte un nuovo insieme di dati. Per questo motivo seguiamo un approccio di validazione in cui splittiamo l’universo dei dati in 2 insiemi disgiunti: uno di questi verrà usato come training set, mentre l’altro verrà usato come insieme di test per misurare la performance del modello addestrato nella fase precedente. Le regole di split più utilizzate sono 1/2 - 1/2 oppure 2/3 - 1/3, dove il primo valore si riferisce alla proporzione di osservazioni utilizzate nel training set; per le restanti osservazioni si prevede la risposta utilizzando il modello stimato e si calcola la funzione di perdita confrontando le previsioni con i veri valori assunti dalla variabile target \(Y\).

Indichiamo con \(T\) l’insieme delle osservazioni di training e con \(V\) l’insieme delle osservazioni di validazione. \(\hat{f}^{-V}\) è il modello stimato con le sole osservazioni di training e, se usiamo una perdita quadratica, l’errore di previsione sarà:

\[\widehat{\mathrm{Err}} = \frac{1}{\# V} \sum_{i \in V} (y_i - \hat{f}^{-V}(x_i))^2 \]

dove \(\# V\) è la cardinalità dell’insieme di validazione. Quando la totalità dei dati \(n\) è molto elevata, è consigliabile ricorrere alla Cross-Validazione semplice. Se non abbiamo abbastanza dati, la perdita d’informazione potrebbe essere troppa… in tal caso ricorriamo alla K-Fold CV o addirittura al Leave One Out CV (caso estremo di K-fold CV con \(K=n\)).

# import dei dati
# fixed X setting 
library(readr)
df <- read_table2("http://azzalini.stat.unipd.it/Book-DM/yesterday.dat", col_types = cols())[-31,]
train <- data.frame(x=df$x, y=df$y.yesterday)
n <- nrow(train)
d <- 3

# splitting rule: 1/2 - 1/2
index <- sample.int(n=n, size=n/2, replace=F)
test <- train[-index,]
train <- train[index,]

# fit sul training set 
mod <- lm(y ~ poly(x, degree=d), data=train) 

# performance sul validation set
yhat <- fitted(mod)
err <- sum((test$y - yhat)^2)/(n/2)

2. K-Fold CV

Generalizzazione del caso precedente e di quello successivo. L’approccio di cross-validazione con K-Folds è quello di usare diverse porzioni di dati come insieme di validazione, una volta ciascuno, e tutte le restanti \(K-1\) come insieme di training. Ogni porzione di dati (fold) contiene lo stesso numero di osservazioni.
N.B. se K=2 allora torniamo al caso precedente.

Ad ogni round uso il k-esimo fold come insieme di validazione, nel senso che prevedo la variabile di risposta per quelle osservazioni usando il modello addestrato con i restanti dati e calcolo il valore della funzione di perdita. La media dei K-valori così calcolati è la stima dell’errore di previsione.

\[\widehat{\mathrm{Err}} = \frac{1}{K} \sum_{k=1}^{K}\left[ \frac{1}{\# V_k} \sum_{i \in V_k} (y_i - \hat{f}^{-V_k}(x_i))^2 \right]\]

Scelte comuni sono K=5 o K=10, ma non esiste una regola sempre valida per ottenere migliori performance previsionali.

rm(list=ls())
library(readr)
df <- read_table2("http://azzalini.stat.unipd.it/Book-DM/yesterday.dat")[-31,]
train <- data.frame(x=df$x, y=df$y.yesterday)
n <- nrow(train)
# K-fold CV
d = 3
K = 5
set.seed(123)
# create folds
folds <- sample( rep(1:K,length=n) )
# initialize vector
KCV <- vector()
# loop
for (k in 1:K){
  fit <- lm(y~poly(x,degree=d), train, subset=which(folds!=k))
  x.out <- train$x[which(folds==k)]
  yhat <- predict(fit, newdata=list(x=x.out))
  y.out <- train$y[which(folds==k)]
  KCV[k]<- mean( ( y.out - yhat )^2 )
}
# KCV estimate 
mean(KCV)
[1] 0.0003160667
# PLOT 
ds <- 1:12 
ps <- ds + 1
library(boot)
set.seed(123)

KCV = sapply(ds, function(d) 
     cv.glm(train, glm(y~poly(x,degree=d), train, family = gaussian), K=K )$delta[1] )

plot(ps, KCV, type="b", log="y")

3. Leave One Out CV

Caso particolare di K-Fold CV con \(K=n\). Si usa 1 osservazione (\(i\)) come insieme di validazione e le restanti \(n-1\) osservazioni come insieme di training per stimare la funzione \(f\) ed ottenere \(\hat{f}^{-i}\); calcolo la funzione di perdita con l’\(i\)-esima osservazione \((y_i - \hat{f}^{-i}(x_i))^2\) e ripeto per ogni \(\quad i=1,\ldots,n\).
Infine prendo la media degli \(n\) valori di perdita ottenuti negli \(n\) insiemi singleton di validazione.
\[\widehat{\mathrm{Err}} = \frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{f}^{-i}(x_i))^2\]

Confronto bias-variance trade-off per K-Fold CV e LOOC :
Le stime con K-Fold CV tendono ad essere più distorte rispetto a quelle con LOOCV perché usano meno informazioni per la fase di train (MSE in media più alto), ma anche con una varianza più bassa in quanto costruite su porzioni di dati non troppo correlati tra loro (mentre con LOOC la varianza è alta perché il training set differisce di 1 sola osservazione ad ogni round).

Short-cut nei modelli lineari per LOOCV :
Per calcolare il LOOCV nel caso di modello lineare esiste una short-cut basata sulla matrice di proiezione \(\underset{n\times n}{\mathbf{H}} = \mathbf{X}(\mathbf{X}^\mathsf{T}\mathbf{X})^{-1}\mathbf{X}^\mathsf{T}\).
Anziché costruire \(n\)-insiemi di validazione si stima una sola volta il modello utilizzando tutti i dati e si aggiusta la funzione di perdita dividendo ogni valore per \((1-h_{ii})^2\) , dove \(h_{ii}\) è l’\(i\)-esimo elemento nella diagonale della matrice di proiezione \({\mathbf{H}}\).

\[\frac{1}{n} \sum_{i=1}^{n}\left( \frac{y_i - \hat{f}(x_i)}{1-h_{ii}} \right)^2\]

rm(list=ls())
library(readr)
df <- read_table2("http://azzalini.stat.unipd.it/Book-DM/yesterday.dat", col_types = cols())[-31,]
train <- data.frame(x=df$x, y=df$y.yesterday)
n <- nrow(train)
d <- 3


# LOOCV
oneout <- vector()
# loop
for (i in 1:n){
fit_i <- lm( y~poly(x,degree=d), data=train[-i,])
yhat_i <- predict(fit_i, newdata=data.frame(x=train$x[i]) )
oneout[i] <- ( train$y[i] -  yhat_i )^2
}
# LOOCV estimate
mean(oneout)
[1] 0.0003439458
# short-cut for linear models 
fit <-  lm(y~poly(x,d), train)
# design matrix
X <- model.matrix(fit)
# projection matrix 
H <- X%*%solve(t(X)%*%X)%*%t(X)
# LOOCV estimate
mean(
  ( (train$y - predict(fit)) / (1-diag(H))  )^2
)
[1] 0.0003439458
# PLOT 
ds <- 1:12 
ps <- ds + 1
library(boot)
set.seed(123)

LOOCV = sapply(ds, function(d) 
     cv.glm(train, glm(y~poly(x,degree=d), 
     train, family = gaussian) )$delta[1] )
plot(ps, LOOCV, type="b")

4. Generalized CV

Nella cross-validazione generalizzata usiamo tutti i dati del train set per calcolare l’\(\mathrm{MSE}_{\mathrm{Tr}}\) e correggiamo questa misura di errore dividendola per un’approssimazione della quantità vista nella short-cut precedente, dove ogni \(h_{ii}\) è approssimato dalla sua media \(\frac{p}{n}\)

Si usa per selezionare il miglior modello per diversi valori del parametro; la formula dell’Errore diventa:

\(\widehat{\mathrm{Err}} =\frac{ \mathrm{MSE}_{\mathrm{Tr}}}{\left(1 - \frac{p}{n}\right)^2}\)

rm(list=ls())
library(readr)
df <- read_table2("http://azzalini.stat.unipd.it/Book-DM/yesterday.dat", col_types = cols())[-31,]
train <- data.frame(x=df$x, y=df$y.yesterday)
n <- nrow(train)
d <- 3
p <- d+1 
fit <- lm(y~poly(x,d), train)
MSE <- deviance(fit)/n
# GCV 
MSE/((1 - p/n))^2
[1] 0.0002776357
# PLOT 
ds <- 1:12 
ps <- ds + 1
library(boot)
set.seed(123)

fun <- function(d) if (d==0) lm(y~1, train) else lm(y~poly(x,degree=d), train)
fits <- lapply(ds, fun)

MSEs.tr <- unlist( lapply(fits, deviance) )/n
GCV = MSEs.tr/(1-(ps)/n )^2

plot(ps, GCV, type="b")