1. k-nearest neighbors

Si tratta di un metodo di stima non-parametrico: non si fanno ipotesi esplicite sulla forma funzionale di \(f()\) ma si lasciano parlare i dati, in modo più flessibile per catturare forme anche strane della relazione funzionale. Il lato negativo, invece, è che occorrono molte più osservazioni per stimare \(f()\) in quanto la stima non dipende più da pochi parametri (sufficienti ad individuare la funzione particolare all’interno della famiglia generale).

L’algoritmo richiede all’utente di specificare un parametro di tuning \(k\) che determina la flessibilità del modello (più è piccolo, più è flessibile). La previsione che viene fatta in un particolare punto \(x^*_1\) è la media dei \(k\)-valori di y delle \(k\)-osservazioni più vicine a \(x^*_1\).

La nozione di vicinanza è data dalla minore distanza Euclidea: prendo le \(k\)-osservazioni con valori \(x_i\) più prossimi a \(x^*_1\) in Norma Euclidea.

\[\| x_i - x^*_1 \|=\sqrt{(x_i - x^*_1)^\mathsf{T}(x_i - x^*_1)}\]

\[\hat{f}(x^*_1) = \frac{1}{k}\sum_{i \in N_k(x^*_1)} y_i\]

dove il simbolo \(N_k(x^*_1)\) indica il set delle \(k\) osservazioni più prossime nello spazio generato da \(X\) al punto in cui avviene la previsione \(x^*_1\).

Poiché trattiamo di distanze, di solito si scalano (normalizzo) e centrano (sottraggo la media) i dati prima di applicare l’algoritmo.

Inoltre per applicare il metodo bisogna specificare una funzione kernel, che nel caso standard non pesato è quella rettangolare. Altre scelte possibili sono triangular, epanechnikov, biweight (or beta(3,3)), triweigth (or beta(4,4)), cos, inv, gaussian, rank, optimal kernels.

Nel Fixed-X setting si dimostra che l’ottimismo di previsione \(OptF\) diventa:
\[\mathrm{OptF} = \mathbb{E}(\mathrm{MSE}_{\mathrm{Te}}) - \mathbb{E}(\mathrm{MSE}_{\mathrm{Tr}}) = \frac{2}{n}\sum_{i=1}^{n}\mathbb{C}\mathrm{ov}(y_i,\hat{f}(x_i))= \frac{2\sigma^2}{k}\] dove l’ultima uguaglianza è giustificata dal fatto che ognuno degli \(n\) termini di covarianza è uguale a \(\frac{\sigma^2}{k}\).

Si stimano per diversi valori del parametro di tuning \(k\) i modelli non parametrici kNNs e si trova il k che rende minimo l’errore di previsione:

\[ \widehat{\mathrm{Err}} = \mathrm{MSE}_{\mathrm{Tr}} + \frac{2\sigma^2}{k}\]

oppure quello che minimizza l’errore riducibile, \([Bias(\hat{f})]^2 + Varianza(\hat{f})\), dove la varianza della funzione stimata è pari alla media delle \(n\) varianze puntuali, ciascuna delle quali è pari a sua volta al rapporto \(\frac{\sigma^2}{k}\).

Nell’esercizio seguente si stimano modelli kNN per diversi valori del parametro di tuning e si sceglie il \(k\) ottimo per cui l’errore riducibile e/o l’errore di previsione sono minimizzati.

set.seed(123)
y <- ftrue + rnorm(n, 0, sigma)
ystar <- ftrue + rnorm(n, 0, sigma)

train <- data.frame(x=x, y=y)
test <- data.frame(x=x, y=ystar)

library(kknn)

# stima per k=21
fit <- kknn(y ~ x, train, test, kernel = "rectangular", k = 21)

# MSEs 
MSEtr <- mean((train$y - fitted(fit))^2)
MSEtest <- mean((test$y - fitted(fit))^2)

# stima per k=1,...,40 e trova k* che rende minimo l'Errore di Previsione
K <- 1:40

fits <- lapply(K, function(k) 
  kknn(y ~ x, train, test, kernel = "rectangular", k = k)
  )

MSEs <- unlist(lapply(fits, function(fit)
  mean((train$y - fitted(fit))^2)
  ))

MSEs_test <- unlist(lapply(fits, function(fit)
  mean((test$y - fitted(fit))^2)
  ))

par(mfrow=c(1,2))
plot(x=1:40, y=MSEs, type="b", col="blue")
plot(x=1:40, y=MSEs_test, type="b", col="red")

which.min(MSEs_test)
[1] 9
ErrF <- MSEs + (2*sigma^2)/(1:40)

which.min(ErrF)
[1] 8
variance <- (sigma^2)/(1:40)

Biass <- unlist(lapply(K, function(k)
  mean((ftrue - fitted(kknn(ftrue ~ x, train, test, kernel = "rectangular", k = k)))^2)
  ))

err_reducible <- Biass + variance

which.min(err_reducible)
[1] 9
barplot(err_reducible, names.arg=1:40)
barplot(variance, add=T, col="lightblue")
legend("topright", c("Bias","Var"), col=c("gray","lightblue"), pch=c(19,19))


# best model according to reducible error : k=9
plot(x=x, y=y, col="blue")
lines(x=x, fitted(kknn(y ~ x, train, test, kernel = "rectangular", k = 9)), lwd=4, col="red")

2. CARET Package for Classification & Regression Training

Si tratta del pacchetto forse più famoso per creare modelli cross-validati, in quanto permette di splittare i dati in train e test, di fare in modo automatizzato il tuning dei parametri e di selezionare - dopo diversi resampling - il miglior modello secondo la misura di perdita specificata.

La funzione di base del pacchetto per le fasi di cui sopra è train, esistono però altre funzioni che useremo come specifica di argomenti opzionali, quali resamples, bwplot e trainControl.

Il pacchetto include anche la funzione featurePlot per la parte di data visualization; si tratta di una funzione wrapper delle principali funzionalità del pacchetto base lattice, con cui si possono realizzare scatterplots personalizzati, confronti tra le variabili continue con faceting, box-plots e curve di densità.

Al link http://topepo.github.io/caret/ si trova la documentazione ufficiale del pacchetto e di tutte le sue funzioni, divise per ambiti applicativi (visualization, pre-processing, data splitting, model training and tuning, feature selection).

# data import 
set.seed(123)

y <- ftrue + rnorm(n, 0, sigma)
ystar <- ftrue + rnorm(n, 0, sigma)
train <- data.frame(x=x, y=y)
test <- data.frame(x=x, y=ystar)

# data visualization 
featurePlot(x=x, y=y, plot="scatter", type=c("p","smooth"), degree=2, span=0.6)

# kNN model training (with Control) & tuning for k
fit.knn <- train(
   y ~ ., train, method = "knn", 
   trControl = trainControl(method = "repeatedcv", number = 5, repeats = 2),
   tuneLength = 10, metric = "RMSE"
)

plot(fit.knn)

fit.knn
k-Nearest Neighbors 

250 samples
  1 predictor

No pre-processing
Resampling: Cross-Validated (5 fold, repeated 2 times) 
Summary of sample sizes: 200, 199, 200, 200, 201, 199, ... 
Resampling results across tuning parameters:

  k   RMSE        Rsquared   MAE        
   5  0.01078883  0.7439864  0.008583497
   7  0.01069480  0.7507754  0.008479432
   9  0.01070968  0.7515992  0.008409578
  11  0.01072294  0.7528524  0.008402358
  13  0.01086012  0.7498961  0.008439059
  15  0.01108839  0.7417333  0.008533102
  17  0.01125582  0.7362683  0.008550564
  19  0.01154106  0.7240076  0.008730682
  21  0.01177913  0.7143308  0.008842280
  23  0.01210702  0.6993483  0.008961112

RMSE was used to select the optimal model using the smallest value.
The final value used for the model was k = 7.
# Polynomial model training (with Control) 
fit.poly <- train(
   y ~ poly(x, degree=5), train, method = "lm", 
   trControl = trainControl(method = "repeatedcv", number = 5, repeats = 2),
  metric = "RMSE"
)

# comparison 
models <- list(
   kNN = fit.knn,
   Polynomial = fit.poly
)

resampling <- resamples(models)

bwplot(resampling, metric="RMSE")  

# performance evaluation on Test Set
yhats <- predict(models, newdata=test)

lapply(yhats, function(yhat) 
   sqrt( mean( (test$y - yhat)^2 ) ) 
   )
$kNN
[1] 0.01091321

$Polynomial
[1] 0.01012645

Commento dei risultati : grazie al pacchetto caret abbiamo addestrato 2 modelli sul training set con metodologia repeated cross validation e abbiamo usato un campione delle loro predictions per vedere la metrica di performance (Root Mean Squared Error) nel test set grazie alle funzioni resampling e bwplot.
Dopodiché abbiamo visto come performano entrambi i modelli nel test set, creando i vettori di previsioni con la funzione predict (che accetta più di un modello in input) e applicando questi vettori alla funzione che definisce l’RMSE.
Anche nel test set vediamo che il modello di regressione polinomiale di 5° grado performa un po’ meglio del modello antagonista di k-Nearest Neighbors.

Si può generalizzare questo esempio a real-world-applications in cui a priori non si sa che tipo di modello potrebbe performare meglio: si addestrano diversi modelli chiamando più volte la funzione train(); questi vengono poi inseriti in una lista che contiene tutti i modelli da cui campioniamo per vedere chi performa meglio nel training e che passiamo come argomento a predict() per generare le previsioni dei vari modelli. In base alla funzione di perdita utilizzata (che può essere custom) vediamo quale modello performa meglio nell’insieme di validazione: sarà quello il best model per l’applicazione in esame.