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

# Datasets

Lavoreremo con un dataset creato appositamente, uno invece già caricato in R (ad esempio ``ChickWeights``) e uno di dati reali (censimento di persone senza tetto per regioni italiane dall'ISTAT).

### Dataset pseudo-casuale

Per generare il dataset, usiamo il comand ``rnorm()``.

```{r}
# Possiamo importare un seed per la riproducibilità, altrimenti ognuno avrà un dataset leggermente diverso
# set.seed(42)

# Generazione di dati fittizi con 200 osservazoni, 
n <- 100
x1 <- rnorm(n, mean = 5, sd = 1)
y1 <- rnorm(n, mean = 5, sd = 1)
x2 <- rnorm(n, mean = 8, sd = 1)
y2 <- rnorm(n, mean = 8, sd = 1)


# Creiamo un dataframe, includiamo anche il gruppo (non useremo per il clustering)
data_random <- data.frame(
  x = c(x1, x2),
  y = c(y1, y2),
  group = factor(rep(1:2, each = n))
)
```

Visualizziamo i dati generati con il comando ``plot()`` di base oppure usando ggplot2.

```{r}
# plot di base
plot(data_random$x, data_random$y, col=data_random$group)

# usando ggplot

library("ggplot2")

ggplot(data_random, aes(x, y, color = group)) +
  geom_point() +
  labs(title = "Dataset generato", x = "X", y = "Y")
```

### Dataset ChickWeight

Questo dataset ``ChickWeight`` (cercare sull'help) contiene informazioni relative alla crescita di polli con diverse diete. Il clustering permette di esplorare i dati ad esempio individuando due (o più) categorie di diete in base all'effetto sulla crescita.


```{r}
head(ChickWeight)
```
Notiamo però che è presente una colonna _Time_ in cui è registrato il tempo dalla crescita dell'individuo colonna _Chick_. Per applicare il clustering  voremmo invece avere solo una riga per individuo contenente i pesi nei vari giorni (0, 2, 4, 6, ecc.). Questo accade spesso con dati reali: non sono strutturati come vorremmo (perché raccolti da altri, ad esempio con altri scopi). Per _ripulire_ i dati, possiamo usare il pacchetto ``tidyr`` della suite ``tidyverse`` (<https://www.tidyverse.org/>). Questo foglio descrive brevemente alcuni comandi principali: <https://leadousset.github.io/intro-to-R/cheatsheet_tidy.pdf>


```{r}
# dare il comando install.packages("tidyverse") se non è già installato
library("tidyr")
```

Nel nostro caso vogliamo _allargare_ il data frame in modo che abbia più colonne (una per ciascuna età). Usiamo quindi il comando ``pivot_wider()``.


```{r}
chick_tidy <- pivot_wider(data=ChickWeight, names_from = Time, values_from = weight)
```

Visualizziamo la tabella _ripulita_.

```{r}
chick_tidy
```

Notiamo che ci sono dei valori mancanti ``NA``. Questo accade spesso con dati reali e bisogna tenerne conto. L'approccio più semplice è di rimuovere completamente le righe per cui almeno una osservazione è mancante. In ``tidyr`` esiste il comando ``drop_na()`` che fa appunto questo.

```{r}
chick_tidy <- drop_na(chick_tidy)

chick_tidy

```
Le righe con almeno un dato mancante sono state rimosse. Attenzione però al rischio di cadere in errori logici come il [Survivorship bias](https://en.wikipedia.org/wiki/Survivorship_bias) oppure il [cherry picking](https://en.wikipedia.org/wiki/Cherry_picking): rimuovendo una porzione di dati rilevanti al problema, si potrebbero trarre conclusioni del tutto errate. Ad esempio, i polli che crescono di più grazie ad una particolare dieta   vengono uccisi una settimana prima e quindi non si registra il peso dell'ultima settimana. Se fosse così, rimuovendo i valori ``NA`` si rimuove proprio il segnale che si sta cercando! 
 
### Dataset senzatetto

Carichiamo infine un dataset reale, proveniente dall'Istituto nazionale di statistica (ISTAT), relativo al numero di persone senza tetto e senza fissa dimora per regione, anno 2021. I dati sono scaricati dalla pagina <https://esploradati.istat.it/databrowser/#/it/censpop/categories/DCSS_SENZA_TETTO_TV/IT1,DF_DCSS_SENZA_TETTO_TV_1,1.0> ma si trovano anche sul team o il sito del corso. 

I dati possono essere salvati in vari formati, il più comune e standard è il formato ``.csv``(*comma-separated values*) in cui (eccetto al più qualche riga di commento iniziare), i dati sono scritti in testo semplice, ciascuna riga relativa ad una riga del data frame e le colonne separate dalla virgola (*comma* in inglese).  Il comando di base per leggere i file in questo formato è ``read.csv()``, mettendo come argomento una stringa con il nome del file (eventualmente nelle sottocartelle del progetto).

```{r}
# dopo aver scaricato il file senza_tetto_italia.csv e averlo copiato in una sottocartella del progetto chiamata datasets, carichiamo il file

senza_tetto_caricato <- read.csv("datasets/senza_tetto_italia.csv")

head(senza_tetto_caricato)

```

Anche in questo caso vediamo che i dati vanno _ripuliti_, selezionando solo le colonne di interesse ed eventualmente allargando le righe per le classi di età. 


#### Digressione su input/output di dati


Notiamo intanto che l'uso della virgola potrebbe essere un problema con i decimali: in inglese si usa invece il punto $\pi= 3.1415..$, mentre in italiano e [altre lingue](https://en.wikipedia.org/wiki/Decimal_separator#/media/File:DecimalSeparator.svg) potrebbe creare letture sbagliate. Per questo ci sono anche formati alternativi, come il ``.tsv`` (*tab-separated values*) in cui si usa una spaziatura ``tab`` per separare i valori. Il comando diventa allora ``read.tsv()``

Altra cosa invece sono i formati di Excel che sono proprietari contengono molte più informazioni, anche circa la storia delle operazioni che sono state effettuate sui dati (ad esempio, un caso di dati _manomessi_ è stato appunto sollevato proprio studiando le operazioni effettuate sull'Excel fornito da un gruppo di autori <https://datacolada.org/109>). Il pacchetto standard dedicato all'input di dati da Excel in R è ``readxl`. Se preferite usare l'interfaccia grafica di RStudio, nel tab _Environment_ trovate il bottone _Import Dataset_ che permette di importare da csv, Excel e anche fare delle prime operazioni sul data frame che verrà assegnato.

Per gestire in modo più semplice l'input di dati da vari formati, personalmente consiglio il pacchetto ``rio`` che ha qualche funzionalità automatizzata per l'input (comando ``import()``) e l'output (comando ``export``).

```{r}
# per installare rio usare il comando install.packages("rio")
library("rio")

# esportiamo il dataset dei polli ripulito in formato excel semplicemente dando l'estensione .xlsx al nome del file

export(chick_tidy, "datasets/polli_puliti.xlsx")

# se navigate nella sottocartella datasets troviamo il file salvato. Possiamo salvare anche in altri formati, digitando l'estensione corretta.

export(chick_tidy, "datasets/polli_puliti.csv")

# con il comando import carichiamo invece da molteplici formati, incluso Excel.

polli_excel <- import("datasets/polli_puliti.xlsx")

head(polli_excel)
```
#### Pulizia dataset senza tetto

Torniamo al dateset senza tetto che abbiamo caricato.

```{r}
head(senza_tetto_caricato)
```

Selezioniamo solo le colonne corrispondenti alla regione _Territorio_ al codice di età e alla frequenza rilevata _Osservazione_.

```{r}

senza_tetto_selezione <- data.frame( "regione" = senza_tetto_caricato$Territorio, "age" = senza_tetto_caricato$AGE_CLASS, "frequenza" = senza_tetto_caricato$Osservazione)

head(senza_tetto_selezione)
```
A questo punto usiamo di nuovo ``pivot_wider`` per creare delle colonne relative alle varie età.

```{r}
senza_tetto_tidy <- pivot_wider(senza_tetto_selezione, names_from = age, values_from = frequenza)

head(senza_tetto_tidy)
```

Possiamo anche salvare il dataset _ripulito_ nel caso ci servisse più avanti (o dovessimo caricarlo, ad esempio per il progetto della proma di esame).

```{r}
export(senza_tetto_tidy, "datasets/senza_tetto_tidy.csv")
```



# Metodi non gerarchici 

Discutiamo prima i metodi di clustering non gerarchici, in particolare K-means ``kmeans()`` e Partitioning Around Medoids (``pam()`` dalla libreria ``cluster``). 

```{r setup, include=FALSE}
# Carichiamo la libreria cluster
library(cluster)

# mostreremo sia i plot di base che i plot ottenuti con ggplot2
library(ggplot2)
```


## K-means clustering



Applichiamo K-means al dataset generato. Togliamo la terza colonna (quella che contiene già il gruppo, per come l'abbiamo generato).

```{r}
head(data_random)
```

```{r}
kmeans_data_random <- kmeans(data_random[-3], centers = 2)

# Il risultato è una lista contenente varie informazioni

kmeans_data_random
```

Ricordiamo che erano 100 punti per gruppo e le medie (che ci aspettiamo essere i centri) erano rispettivamente $(5,5)$, $(8,8)$.  Aggiungiamo la colonna del cluster al data frame originale.

```{r}

data_random <- cbind( data_random, "kmeans" =  factor(kmeans_data_random$cluster ))
```

Possiamo confrontare quindi con una tabella di contingenza ``table()`` i due cluster (originale e quello di kmeans).

```{r}
table( "originale" = data_random$group, "kmeans"=data_random$kmeans)
```
Ovviamente non è detto che il cluster etichettato con $1$ da k-means corrisponda con la classe $1$ originale!

Possiamo fare un plot usando come colore il cluster trovato con k-means.

```{r}

plot(data_random$x, data_random$y, col=data_random$kmeans)
```

Con ``ggplot`` possiamo assegnare la forma del gruppo originale e il colore del cluster trovato con k-means. Riuscite a vedere i punti non _classificati correttamente_?

```{r}
ggplot(data = data_random, aes(x=x, y=y, colour = kmeans, shape = group)) + geom_point()
```

Come determinare $k$? l'output di k-means fornisce il WCSS per cluster e il WCSS totale.
```{r}

kmeans_data_random$withinss

kmeans_data_random$tot.withinss

```

Plottiamo il WCSS per vari valori di $k$ (cercando di applicare il cosiddetto _metodo elbow_).

```{r}

WCSS <- c()

for(k in 2:10){
  WCSS <- c(WCSS, kmeans(data_random[1:2], centers=k)$tot.withinss)
}

plot(2:10,  WCSS, type='l', xlab="numero di clusters", ylab="WCSS")

```

Provate a riapplicare il blocco di codice sopra. Cosa notate? L'algoritmo di k-means (e pure pam) è in realtà molto sensibile ai centri iniziali. Per questo è meglio applicare il metodo diverse istanze, indicare il valore medio e un intervallo di confidenza (in questo caso bilatero al livello $95\%$).

```{r}

WCSS <- data.frame("k"=numeric(), "mean"=numeric(), "sd" =numeric() )

numero_runs <- 10

for(k in 2:10){
  wcss_runs <- c()
  for (i in 1:numero_runs){
    wcss_runs <- c(wcss_runs,  kmeans(data_random[1:2], centers=k)$tot.withinss)
  }
  WCSS <- rbind( WCSS, data.frame(k, mean(wcss_runs), sd(wcss_runs)))
}

plot(2:10,  WCSS$mean, type='l', xlab="numero di clusters", ylab="WCSS", col="red")
lines(2:10, WCSS$mean+qt(0.975, df=numero_runs-1) * WCSS$sd/sqrt(numero_runs), col="grey")
lines(2:10, WCSS$mean-qt(0.975, df=numero_runs-1) * WCSS$sd/sqrt(numero_runs), col="grey")
```

Per **esercizio**, creare un dataset in cui il numero di *clusters* sia diverso da 2 e verificare l'andamento del WCSS riconoscendo se possibile il punto di *gomito*.


## Partitioning Around Medoids (PAM)

Applichiamo ora PAM e confrontiamo il risultato con k-means al dataset ``chick_tidy``. Il comando è ``pam()`` dal pacchetto ``cluster``. Come kmeans, bisogna specificare $k$. Togliamo le prime due colonne (numero identificativo individuo e tipo di dieta).

```{r}

pam_chick <- pam(chick_tidy[-(1:2)], k = 2)

pam_chick
```

Possiamo anche in questo case aggiungere il vettore di clustering trovato al data frame.

```{r}

chick_tidy$pam_cluster <- pam_chick$clustering

```

In questo caso non abbiamo delle classi con cui naturalmente confrontare il risultato, quindi possiamo ad esempio confrontare con $k$-means.

```{r}

# rimuoviamo anche l'ultima colonna (che contiene il clustering di pam)

kmeans_chick <- kmeans(chick_tidy[-c(1, 2, 15)],  2)

kmeans_chick
```

Aggiungiamo anche questi risultati al data frame e confrontiamo con una tabella di contingenza.

```{r}

chick_tidy$kmeans_cluster <- kmeans_chick$cluster

table("pam"=chick_tidy$pam_cluster, "kmeans"=chick_tidy$kmeans_cluster)
```

Vediamo che le classi trovate sono molto simili. La scelta dei medoid può influenzare i risultati, ma PAM è generalmente più stabile in presenza di rumore.


Un vantaggio notevole di ``pam`` rispetto a k-means è la possibilità di usare metriche diverse. È possibile dare come input invece del data frame una matrice di dissimilarità, oppure specificare ``metric = "manhattan"`` per usare la distanza $\ell_1$ (taxicab). 

```{r}

pam_chick <- pam(chick_tidy[-c(1,2, 15, 16)], k = 2, metric="manhattan")

chick_tidy$pam_manhattan_cluster <- pam_chick$cluster

```

Confrontiamo i risultati di pam con le due metriche.

```{r}
table("pam"=chick_tidy$pam_cluster, "pam_manhattan"=chick_tidy$pam_manhattan_cluster)

```

Possiamo infine usare la silhouette per valutare i cluster trovati (anche per lo stesso $k$). Il comando è ``silhouette()`` dal pacchetto ``cluster``, a cui bisogna dare come input il vettore del cluster e una matrice di dissimilarità (ad esempio calcolata con la funzione ``dist()``)

```{r}
sil_pam <- silhouette(chick_tidy$pam_cluster, dist(chick_tidy[-c(1,2,15, 16, 17)]))

sil_pam_manhattan <- silhouette(chick_tidy$pam_manhattan_cluster, dist(chick_tidy[-c(1,2,15, 16, 17)]))
                                
sil_kmeans <- silhouette(chick_tidy$kmeans_cluster, dist(chick_tidy[-c(1,2,15, 16, 17)]))
```

Possiamo plottare la silhouette per ciascun metodo o per ciascun cluster, oppure limitarci alla silhouette media.

```{r}

plot(sil_pam)

boxplot(sil_pam[, 3])

# con la funzione summary otteniamo un riassunto delle varie silhouette in cui possiamo visualizzare la silhouette media

summary(sil_pam)

```

Confrontiamo le tre silhouette medie.

```{r}

print(paste("silhouette media per PAM:", round(mean(sil_pam[, 3]), 4)))

print(paste("silhouette media per PAM Manhattan:", round( mean(sil_pam_manhattan[,3]), 4)))

print(paste("silhouette media per k-means:", round(mean(sil_kmeans[,3]), 4)))


```
Per **esercizio**: completate aggiungendo la deviazione standard delle silhouette calcolate.



# Metodi gerarchici


## Agnes (agglomerativo)

Il comando ``hclust()`` permette di utilizzare diversi metodi. In alternativa, possiamo essere più specifici e usare ``agnes()`` del pacchetto ``cluster``.


Consideriamo il dataset relativo alle persone senza tetto e usiamo  ``agnes()``, che di base usa il metodo di _average linkage_ con distanza Euclidea (togliamo la prima colonna che contiene i nomi).

```{r}

agnes_senza_tetto <- agnes(senza_tetto_tidy[-1])
```

Possiamo visualizzare il risultato con il dendrogramma.

```{r}
plot(agnes_senza_tetto)
```
Per visualizzare il plot possiamo usare prima colonna come nome delle righe del data frame e usare ``ggdendro``(estensione di ``ggplot2``).

```{r}
senza_tetto_labels <- as.data.frame(senza_tetto_tidy[-1])
rownames(senza_tetto_labels) <- senza_tetto_tidy$regione

#ggdendro permette di visualizzare meglio i dendrogrammi usando la grammatica di ggplot2

library("ggdendro")

dg <- dendro_data(agnes(senza_tetto_labels))

ggdendrogram(dg)

```

Una osservazione _importante_ sui dati: stiamo confrontando le frequenze assolute dei senza tetto, quindi è naturale che le regioni più popolose avranno più persone senza fissa dimora (stiamo quindi implicitamente facendo clustering in base alla popolazione totale nella regione). 

**Esercizio**: recuperare dal sito ISTAT il numero di abitanti per regione ed eseguire un clustering usando la frequenza relativa dei senza tetto sulla popolazione totale. Confrontare i dendrogrammi ottenuti.

Per ovviare a questo problema, effuttiamo clustering soltanto sulla frequenza _relativa_ della popolazione nelle varie classi di età.

```{r}
senza_tetto_relative <- data.frame( senza_tetto_labels[1:4]/senza_tetto_labels[,5])

agnes_senza_tetto_relative <- agnes(senza_tetto_relative)

ggdendrogram(dendro_data(agnes_senza_tetto_relative))

```


Tornando al problema, il comando ``cutree()`` permette di ricavare il vero e proprio clustering _tagliando_ il dendrogramma: si può specificare l'altezza o il numero di cluster desiderati. Consideriamo ad esempio $k=5$ clusters.

```{r}

senza_tetto_relative$agnes_k_5 <- cutree(agnes_senza_tetto_relative, k=5)

```

Studiamo la silhouette: notiamo un valore medio non molto alto.

```{r}
sil_agnes_senza_tetto <- silhouette(senza_tetto_relative$agnes_k_5, dist(senza_tetto_relative[1:4]))

summary(sil_agnes_senza_tetto)
```

## Diana (divisivo) 

Con il comando ``diana()`` applichiamo il metodo divisivo. Vediamo che non ci sono grandi differenze.

```{r}

diana_senza_tetto_relative <- diana(senza_tetto_relative)

ggdendrogram(dendro_data(diana_senza_tetto_relative))

```

Calcoliamo al solito la silhouette, per $k=5$.

```{r}

senza_tetto_relative$diana_k_5 = cutree(diana_senza_tetto_relative, k=5)

sil_diana_senza_tetto <- silhouette(senza_tetto_relative$diana_k_5, dist(senza_tetto_relative[1:4]))

summary(sil_diana_senza_tetto)

```

Confrontiamo con un metodo non gerarchico, $k$-means.


```{r}
senza_tetto_relative$kmeans <- kmeans(senza_tetto_relative, 5)$cluster

summary(silhouette(senza_tetto_relative$kmeans,dist(senza_tetto_relative[1:4]) ))

```

Confrontiamo i plot trovati con una tabella di contingenza.

```{r}
table(data.frame("kmeans"= factor(senza_tetto_relative$kmeans), "agnes"=factor(senza_tetto_relative$agnes_k_5)))

```

# Esercizi

1. Generare una tabella di 3 colonne e 120 righe, in modo tale che la terza colonna indichi l’appartenenza ad un cluster, e sia pari a 1 per le prime 50 righe e pari a 2 per le ultime 70 righe. Implementare il calcolo diretto della silhouette dell’individuo corrispondente alla prima riga, usando come distanza tra individui la distanza euclidea tra i punti le cui coordinate sono i fattori degli individui.

2. Generare una tabella in modo che la silhouette ottenuta come risultato di una analisi di clustering mostri la pessima attribuzione di un individuo.

3. Generare una tabella in modo che la silhouette ottenuta come risultato di una analisi di clustering mostri il pessimo punteggio di un cluster.

4. Generare un campione i cui individui siano caratterizzati da 6
diverse caratteristiche, e tali che in una analisi di clustering tipo pam la scelta di un numero di cluster inferiore a 4 non risulti buona. Implementare anche l’analisi della bontà del metodo.

5. Svolgere una analisi di clustering sul dataset ``USArrests`` utilizzando il metodo partition around medoids con distanza manhattan.

6. Analizzare il problema del clustering sul dataset ``USArrests`` utilizzando metodi gerarchici.

7. Analizzare il problema del clustering sul dataset ``iris`` utilizzando metodi gerarchici, valutando i differenti casi ottenuti al variare delle possibili distanze tra punti e tra cluster.

8. Analizzare il problema del clustering sul dataset ``votes.repub`` utilizzando metodi gerarchici, valutando i differenti casi ottenuti al variare delle possibili distanze tra punti e tra cluster.

9. Analizzare il problema del clustering per il dataset ``agriculture`` (presente nel pacchetto cluster), relativo a dati su PIL e percentuale di impiegati nell’agricoltura nei paesi UE nel 1993.

10. Analizzare il problema del clustering per il dataset ``flower`` (presente nel pacchetto cluster), relativo a otto caratteristiche di alcuni fiori.

11. Analizzare il problema del clustering per il dataset ``UScereals``.
