---
title: "Classificazione I (notebook 4)"
author: "Dario Trevisan"
date: "15/10/2025"
output:
  html_notebook:
    toc: true
    toc_depth: 3
    toc_float: true
    theme: readable
    df_print: paged
    download_handler: true
  html_document:
    toc: true
    toc_depth: '3'
    df_print: paged
subtitle: "Statistica II - 750AA"
---




# Introduzione

In questo notebook R iniziamo a studiare la classificazione statistica su dataset generati e un dataset reale (pre-caricato). Il notebook si concentra sull'utilizzo dell'algoritmo $k$-NN, ma allo stesso tempo discute concetti generali come errori di training/test, cross-validation, matrici di confusione e curve ROC, che saranno poi utili anche per altri metodi di classificazione.

Iniziamo caricando le librerie necessarie -- usare ``install.packages("nome pacchetto")`` per installare quelle non già presenti.

```{r}
# ggplot2 per grafiche migliorate
library(ggplot2)

# libreria con comandi per classificazione di base
library(class)

# libreria con comandi avanzati per classificazione e regressione
library(caret)


# libreria per il dataset medicale
library(MedDataSets)
```

## Dataset Generati

In questa sezione creiamo un dataset fittizio, usando comandi di campionamento come ``rnorm`` ma anche ``rbinom``, ``rpois`` per avere caratteristiche discrete.

```{r}

# Creazione di un dataset fittizio 
# impostare il seed per riproducibilità
#set.seed(42)

# numero di osservazioni
n_neg <- 200 # numero di negativi
n_pos <- 100 # numero di positivi

# n è la taglia del campione
n <- n_neg+n_pos

# generiamo le features
x1 <- rnorm(n_neg, mean=-1, sd=1)
x2 <- rnorm(n_neg, mean=-1, sd=1)
x3 <- rnorm(n_pos, mean=1, sd=1)
x4 <- rnorm(n_pos, mean=1, sd=1)

# generiamo le classi
classe <- factor(c(rep(-1, n_neg), rep(1, n_pos) ))

data_generated <- data.frame(c(x1, x3), c(x2, x4), classe)

colnames(data_generated) <- c("feature_1", "feature_2", "label")

# Visualizzazione del dataset
ggplot(data_generated, aes(x =feature_1, y = feature_2, color = label)) +
  geom_point() +
  labs(title = "Dataset Fittizio", x = "feature 1", y = "feature 2")
```

Come discusso a lezione, procediamo inizialmente dividendo il dataset in training e test sets. Scegliamo una percentuale significativa per il training ($80\%$ ad esempio, ma comunque anche da valutare in base a quanti dati disponiamo). 

Usiamo la funzione ``sample()`` per campionare in modo casuale, per evitare di prendere ad esempio dati da una stessa classe se già ordinati per classe (come nel nostro caso). 


```{r}
# Suddivisione del dataset in training e test set
#set.seed(42)

train_index <- sample(1:n, size = 0.8 * n)


train_generated <- data_generated[train_index, ]
test_generated <- data_generated[-train_index, ]

head(train_generated)
```

Visualizziamo i dati distinguendo train e test usando una geometria diversa.

```{r}

# creiamo un vettore che indica se il dato è test o train

is_train <- rep(FALSE, n)
is_train[train_index] <- TRUE

# Visualizzazione del dataset distinguendo tra train e test
ggplot(test_generated, aes(x =feature_1, y = feature_2, color = label) )+ #shape=is_train )) +
  geom_point() +
  labs(title = "Dataset Fittizio", x = "feature 1", y = "feature 2") # shape="è training point?")

```

Ora *alleniamo* il modello $k$-NN sui dati del training set e lo valutiamo sui dati del test set generato. Usiamo la funzione ``knn()`` dal pacchetto ``class``. In alternativa si può usare ``knn3()`` dal pacchetto ``caret``.



```{r}
# Allenamento del modello KNN

# Scegliamo il valore di k che preferiamo, meglio dispari per classificazione binaria altrimenti deve risolvere le parità a caso
k <- 1

punto_test <- c(0,0)
# salviamo le previsioni
knn_pred <- knn(train_generated[, 1:2], test_generated[, 1:2], train_generated$label, k = k)

table(knn_pred)
```

Visualizziamo le etichette attribuite al test set (usiamo delle forme) confrontandole con quelle originali (usiamo i colori originali).

```{r}

test_generated$knn <- knn_pred

ggplot(test_generated, aes(x =feature_1, y = feature_2, color = label, shape=knn )) +
  geom_point() +
  labs(title = "Dataset Fittizio", x = "feature 1", y = "feature 2", shape="knn")
```

### Indicatori di performance

Una metrica generale è l'accuratezza sul test set, calcoliamola confrontando le previsioni con le vere etichette.

```{r}
# accuratezza

accuracy <- mean(test_generated$knn== test_generated$label )

round(accuracy, 6)

```

Il *test error* è quindi `r round(1-accuracy,2)`.

Trattandosi di una classificazione binaria valutiamo la performance in una *matrice di confusione* (questo si può fare anche per classificazioni non binarie, ma in tal caso la taglia della matrice aumenta ed è meno riassuntiva).


```{r}
# Matrice di confusione
conf_matrix <- table("Reale" = test_generated$label, "Classificato" = test_generated$knn)
print(conf_matrix)
```

Convenientemente abbiamo chiamato le classi $-1$ per i Negativi, $1$ per i Positivi, ma possiamo anche cambiare i nomi alle righe e colonne della tabella.

```{r}
colnames(conf_matrix) <- c("CN", "CP")
row.names(conf_matrix) <-  c("N", "P")

print(conf_matrix)

```

Estraiamo *a mano* le componenti della matrice e calcoliamo le metriche principali

```{r}

reali <- test_generated$label
previsti <- test_generated$knn

# positivi reali
P <- sum( reali == 1)
# negativi reali
N <- sum( reali == -1)

# prevalenza
prevalence <- P/(P+N) 

# veri negativi
TN <- sum( reali == -1 & previsti == -1)
# veri positivi
TP <-sum( reali == 1 & previsti == 1)
# falsi negativi, mancato allarme
FN <- sum( reali == 1 &  previsti == -1)
# falsi positivi, falso allarme
FP <- sum( reali == -1 & previsti == 1)

# precisione (probabilità che sia vero positivo sapendo che il test indica positivo)
precision <- TP/(TP+FP)

# TPR, recall/richiamo, sensibilità (probabilità che il test indichi positivo sapendo che il soggetto è realmente positivo)
TPR <- TP/P
recall <- TPR

# TNR, specificità  (probabilità che il test indichi negativo sapendo che è realmente negativo)

TNR <- TN/N

# FPR, errore di tipo I (probabilità che il test indichi positivo sapendo che è realmente negativo)

FPR <- FP/N

# FNR, errore di tipo II (probabilità che il test indichi negativo sapendo che è realmente positivo)

FNR <- FN/P

```


Possiamo rappresentare il risultato del classificatore sul test set come un punto nel grafico FPR/TPR. Più il punto è vicino all'angolo in alto a sinistra $(0,1)$ migliore è la performance del classificatore (sul test set).

```{r}

ROC <- data.frame( "FPR" = FPR, "TPR" = TPR, "k" = k)

ggplot(data = ROC, aes(x=FPR, y=TPR, colour = k))+geom_point()+xlim(0,1)+ylim(0,1)+geom_abline(slope=1, intercept = 0, col="red")

# togliamo il commento per aggiungere la bisettrice per indicare la performance dei classificatori puramente casuali

```

Ci sono altri indicatori che cercano di *riassumere* la performance di un classificatore binario (o meglio, di un test) con un singolo numero. Uno è anche l'accuratezza, ma potrebbe essere influenzata dalla bassa prevalenza. Indicatori migliori in questo senso sono la *misura $F_1$* ($F_1$ score) che media  tra precisione e richiamo. Idealmente vorremmo sia recisione che recall vicine ad $1$, quindi mediandole avremo un indicatore che più si avvicina ad $1$ meglio è. La media che viene utilizzata è la media armonica:
$$ F_1 =\left( \frac{precision^{-1} + recall^{-1}}{2} \right)^{-1} $$

```{r}
F1 <- 2/(1/precision+1/recall)
```

La misura $F_1$ in questo caso vale `r F1`. Più in generale è possibile introdurre un *peso* relativo tra precisione e richiamo con la misura $F_\beta$, dove $\beta>0$ indica un fattore che pesa il richiamo $\beta$ volte più importante della precisione.

$$ F_\beta =\left( \frac{\beta^{-1} \cdot precision^{-1} + \beta \cdot recall^{-1}}{\beta^{-1}+\beta} \right)^{-1} $$
Ad esempio, si usa  $F_2$ che pone più peso al richiamo, oppure il caso $F_{1/2}$ che pone più peso alla precisione.


```{r}
F_beta <- function(precision, recall, beta=1){
  1/((1/(beta *precision) + beta/recall)/(1/beta+beta))
}

F1 <- round(F_beta(precision, recall), 2)
F2 <- round(F_beta(precision, recall, beta=2), 2)
F0.5 <- round(F_beta(precision, recall, beta=0.5), 2)

```

Troviamo in questo caso $F_1=$ `r F1`, $F_2=$ `r F2`, $F_{1/2}=$ `r F0.5` (in questo caso precision = `r round(precision, 2)`, recall=`r round(recall, 2)`).

Un ulteriore indicatore è il coefficiente $\phi$, noto anche come *Matthews correlation coefficient (MCC)* in machine learning. Per definizione, si tratta di calcolare il coefficiente di correlazione di Pearson tra le variabili aleatorie (a valori in $-1$, $1$) con densità discreta di probabilità congiunta data dalle entrate della matrice di confusione divise per la taglia del test set. Si trova la formula
$$ \varphi = \frac{ TP \cdot TN - FN \cdot FN}{\sqrt{ P \cdot N \cdot (TP+FP) \cdot (TN+FN)}}. $$

Trattandosi di un coefficiente di correlazione, valori estremi sono vicini a $1$ o $-1$ (e indicano una forte associazone tra previsione ed etichette reali), mentre valori vicini a $0$ indicano una debole associazione. Tipicamente si osservano valori positivi comunque (un classificatore con valore $\phi$ negativo indica che sbaglia sistematicamente).


```{r}

phi = (TP*TN - FP*FN)/(sqrt( P*N* (TP+FP)*(TN+FN) ) )

```

In questo caso il coefficiente vale $\varphi =$ `r round(phi, 2)`.

Non è necessario scrivere il codice per calcolarsi tutti questi indicatori ogni volta *ex novo*: possiamo usare la funzione ``confusionMatrix()`` dalla libreria ``caret``.

```{r}
#carichiamo la libreria caret se non già fatto sopra
library(caret)

conf_matrix <- confusionMatrix(test_generated$knn, test_generated$label,positive="1")

# visualizziamo la matrice di confusione
print(conf_matrix$table)

# gli indicatori principali sono pure calcolati (e accessibili singolarmente)

print(conf_matrix)

```

Nel pacchetto ``caret`` sono disponiblii anche ulteriori funzioni, come ``recall()``, ``precision()`` e ``F_meas()`` dal nome autoesplicativo (leggerne l'help per informazioni più precise).


```{r}

precision(conf_matrix$table, relevant="1")

recall(conf_matrix$table, relevant="1")

F_meas(conf_matrix$table, relevant="1", beta=1)
```



### Cross-Validation

La cross-validation divide il training set in più sottogruppi, utilizzando alcuni di essi per addestrare il modello e altri per testarlo. Questo processo viene ripetuto più volte, cambiando ad ogni iterazione i dati utilizzati per l'addestramento e il test, in modo da avere uno stimatore dell'errore di generalizzazione, il rischio atteso.

Per la funzione ``knn()`` del pacchetto ``class`` è già implementata la *leave-one-out* CV, con la funzione ``knn.cv()``.

```{r}

k <- 3
generated_knn_cv <- knn.cv(train_generated[, 1:2], train_generated$label, k)

# aggiungiamo il risultato della cv come colonna

train_generated$knn_cv <- generated_knn_cv

# calcoliamo l'errore di cross validation

error_cv_generated = mean(train_generated$knn_cv != train_generated$label)

# confrontiamolo con l'errore di training e di test

error_train_generated = mean( knn(train_generated[, 1:2], train_generated[, 1:2], train_generated$label, k) != train_generated$label)

error_test_generated = mean( knn(train_generated[, 1:2], test_generated[, 1:2], train_generated$label, k) != test_generated$label)

print(paste("errore di cv (LoO):", round(error_cv_generated, 2)))
print(paste("errore di training:", round(error_train_generated, 2)))
print(paste("errore di test:", round(error_test_generated, 2)))


```


Un approccio comune è la $k$-fold cross-validation, dove il dataset viene suddiviso in $k$ parti, e il modello viene addestrato $k$ volte, ogni volta utilizzando un diverso fold come set di test.


```{r}

# impostiamo il k della k-fold CV.

k_cv <- 10

n_train <- nrow(train_generated)

# prendiamo una permutazione casuale degli indici

permutation_train_index <- sample(1:n_train)

folds <- round(seq(1, n_train, length.out=k_cv), 0)


#carichiamo gli errori in un vettore

error_cv_kfold <- c()

for( i in 1:(length(folds)-1)){
  # seleziona gli indici della fold i
  fold_index <- permutation_train_index[folds[i]:(folds[i+1]-1)]
  # definisce il test/train set per la fold corrente
  test_fold_generated <- train_generated[fold_index, ]
  train_fold_generated <- train_generated[-fold_index, ]
  
  #calcola l'errore
  error_fold <-  mean( knn(train_fold_generated[, 1:2], test_fold_generated[, 1:2], train_fold_generated$label, k) != test_fold_generated$label)
  
  # aggiunge l'errore alla sequenza
  error_cv_kfold <- c(error_cv_kfold, error_fold)
}

# calcola la media degli errori
mean(error_cv_kfold)

```


Anche per la CV, tuttavia, è possibile e *fortemente consigliato* utilizzare delle funzioni già predisposte in R. Nel pacchetto ``caret`` possiamo usare la funzione ``train()`` e ``trainControl()`` specificando il metodo ``cv`` (altri metodi sono pure già implementati)


```{r}
# Esecuzione della cross-validation
#set.seed(42)

#specifichiamo i valori per cui calcolare l'errore di CV

control <- trainControl(method = "cv", number = 10)

# usare method = "LOOCV" per leave-one-out CV

knn_cv <- train(train_generated[,1:2], train_generated$label, method = "knn", trControl = control, tuneLength = 20)

# usare l'opzione tuneGrid = data.frame("k"=1:10) per specificare i valori di k da esplorare

# Visualizzazione dei risultati della cross-validation
print(knn_cv)
plot(knn_cv)
# oppure ggplot(knn_cv)
```


## Dataset Precaricato

Applichiamo il clustering con il metodo $k$-NN al dataset ``Melanoma_df`` del pacchetto ``MedDataSets``.

```{r}
#installare il pacchetto install.packages("MedDataSets")

mel_data <- MedDataSets::Melanoma_df

head(mel_data)
```
Un obiettivo dell'analisi di classificazione potrebbe essere di determinare se il paziente sopravvive oltre $3$ anni dalla diagnosi (che possiamo ottenere dalla prima colonna) in base alle caratteristiche come il sesso, l'età, lo spessore e la presenza di ulcerazioni.


```{r}
#introduciamo una colonna per indicare se il paziente è sopravvive a 3 anni dalla diagnosi

mel_data$label <- mel_data$time > (3*365)

mel_data_tidy <- mel_data[c(3, 4, 6, 7, 8)]

head(mel_data_tidy)

# rappresentiamo con un boxplot l'associazione tra spessore e label
ggplot(mel_data_tidy, aes(y = thickness, group= label, fill=label))+geom_boxplot()

# rappresentiamo con una tabella di contingenza l'associazone tra sesso e label

table("sesso" = mel_data_tidy$sex, "sopravvivenza" = mel_data_tidy$label)
```

Passiamo alla classificazione con $k$-NN. Poiché il metodo si basa sulle distanze (euclidee) tra le caratteristiche, è buona pratica *standardizzare* i dati, quantomeno le variabili numeriche. **Per esercizio:** ripetere l'analisi senza standardizzare i dati, oppure standardizzando tutte le variabili.

```{r}

mel_data_tidy[2:3] <- scale(mel_data_tidy[2:3])

head(mel_data_tidy)
```
Dividiamo in training e test set.

```{r}

n <- nrow(mel_data_tidy)

train_index <- sample(1:n, size = 0.8 * n)


train_mel <- mel_data_tidy[train_index, ]
test_mel <- mel_data_tidy[-train_index, ]

```

Alleniamo intanto il modello con un $k$ scelto da noi.


```{r}
# k scelto da noi
k <- 5

# salviamo le previsioni
knn_pred <- knn(train_mel[, -5], test_mel[, -5], train_mel$label, k = k)

```

Valutiamo la performance con la matrice di confusione. In questo caso impostiamo come *positivo* un paziente che *non* sopravvive a $5$ anni dalla diagnosi (vogliamo che il nostro classificatore fornisca una valutazione del rischio di recivida/complicazioni ecc.).

```{r}


conf_matrix <- confusionMatrix(as.factor(knn_pred), as.factor(test_mel$label),positive="FALSE")

# visualizziamo la matrice di confusione
print(conf_matrix$table)

# e gli indicatori principali sono pure calcolati (e accessibili singolarmente)

print(conf_matrix)
```

Vediamo l'$F_1$ score.
```{r}
 F1 <- F_meas(conf_matrix$table, relevant="FALSE", beta=1)
F1
```

Plottiamo il risultato nel piano FPR/TPR.

```{r}

# TPR è la sensitivity
TPR = conf_matrix$byClass[1]
# FPR è 1-specificity
FPR = 1-conf_matrix$byClass[2]

ROC <- data.frame( "FPR" = FPR, "TPR" = TPR, "k" = k)

ggplot(data = ROC, aes(x=FPR, y=TPR, colour = k))+geom_point()+xlim(0,1)+ylim(0,1) +geom_abline(slope=1, intercept = 0, col="red")


```
Per il valore di $k$ scelto il test performa male per la TPR (sensibilità). Esploriamo altri valori di $k$ usando però la cross validation.

```{r}

k_max = 10

performance <- data.frame( "FPR" = c(), "TPR" = c(), "accuracy"=c(), "F1" = c(), "k" = c())

for( k in 1:k_max){
  knn_cv <- knn.cv(train_mel[, 1:4], train_mel$label, k = k)
  
  accuracy <- mean( knn_cv == train_mel$label)
  
  conf_matrix <- confusionMatrix(as.factor(knn_cv), as.factor(train_mel$label),positive="FALSE")
  # TPR è la sensitivity
  TPR <- conf_matrix$byClass[1]
  # FPR è 1-specificity
  FPR <- 1-conf_matrix$byClass[2]
  
  F1 <- F_meas(conf_matrix$table, relevant="FALSE", beta=1)


  performance = rbind(performance, data.frame("FPR" = FPR, "TPR" = TPR, "accuracy"=accuracy, "F1"=F1, "k" = k) )
}

ggplot(performance, aes(x=k, y=accuracy))+geom_line()+ggtitle("validation error (LoOCV) al variare di k, dataset melanoma_df")

```
Eseguiamo una cross validation con $k$-fold, $k=10$.

```{r}

control <- trainControl(method = "cv", number = 10)

# usare method = "LOOCV" per leave-one-out CV

knn_cv <- train(train_mel[, -5], factor(train_mel$label), method = "knn", trControl = control,  tuneGrid = data.frame("k"=1:10))

# usare l'opzione tuneGrid = data.frame("k"=1:10) per specificare i valori di k da esplorare

# Visualizzazione dei risultati della cross-validation
plot(knn_cv)

```


Sembra che $k=3$ dia una accuratezza leggermente preferibile, stando ad entrambe le CV. Plottiamo gli $F$-score al variare di $k$.

```{r}
ggplot(performance, aes(x=k, y=F1))+geom_line()+ggtitle("F1 score al variare di k, dataset melanoma_df")
```

Plottiamo infine i punti sul piano FPR/TPR al variare di k.

```{r}

ggplot(data = performance, aes(x=FPR, y=TPR, colour = k, label=k))+geom_text(hjust=0, vjust=0)+xlim(0,1)+ylim(0,1) +geom_abline(slope=1, intercept = 0, col="red")


```

La cross-validation e gli altri indicatori ci suggeriscono di passare a $k=2$. Scegliamo quindi questo modello e valutiamo infine il test error.


```{r}
# k individuato tramite CV

k <- 2

# salviamo le previsioni
knn_pred <- knn(train_mel[, -5], test_mel[, -5], train_mel$label, k = k)


conf_matrix <- confusionMatrix(as.factor(knn_pred), as.factor(test_mel$label),positive="FALSE")

# visualizziamo la matrice di confusione
print(conf_matrix$table)

# e gli indicatori principali sono pure calcolati (e accessibili singolarmente)

print(conf_matrix)

```


## Esercizi


1.Genera un dataset fittizio con due classi utilizzando distribuzioni normali con deviazioni standard diverse (e medie diverse). Ogni classe deve avere una dimensione di 100 campioni. Usa ggplot2 per visualizzare il dataset. Calcola l'accuratezza del modello KNN e seleziona un k ottimale basandoti sugli indicatori.

2. Utilizza il dataset Iris per applicare il KNN (label è la colonna delle specie). Suddividi il dataset in un set di training (70%) e un set di test (30%). Calcola la matrice di confusione e valuta l'accuratezza del modello. 


3. Utilizza il dataset mtcars per applicare il KNN e ottimizzare il valore di k attraverso la k-fold cross-validation (ad esempio 10-fold). Traccia il grafico dell'accuratezza media in funzione dei diversi valori di $k$ e identifica il valore ottimale.

4. Usa il dataset ``wine`` (disponibile online <https://www.kaggle.com/datasets/yasserh/wine-quality-dataset> o nel pacchetto ``rattle``) per applicare il KNN. Suddividi il dataset in un set di training e uno di test. Calcola la precisione, il richiamo e il punteggio F1 per ciascuna classe. Visualizza il risultato con una heatmap.

5.  Utilizza il dataset ``mtcars`` e confronta le prestazioni del modello KNN su dati non standardizzati e standardizzati. Calcola e confronta l'accuratezza del modello su entrambi i set di dati. Discutere come la standardizzazione influisce sulle prestazioni del KNN.


6.  Utilizza il dataset ``iris`` per applicare la PCA e ridurre la dimensionalità a due componenti principali. Successivamente, utilizza KNN per classificare le specie. Confronta la performance con il risultato senza PCA.

7. Usa il dataset ``mtcars`` per applicare un algoritmo di clustering (ad esempio K-means) per identificare gruppi di automobili. Successivamente, utilizza un modello di classificazione per prevedere il gruppo a cui appartiene ogni automobile.