### Cluster Analysis: K-means & PAM

## 30/04/26
## K-Means

#install.packages("cluster")
library(cluster)
dat<-read.csv("food_spending.csv", header=TRUE)
head(dat)
plot(dat$FreqAway, dat$Home)

# funzione kmeans()
?kmeans
kmeans(dat, centers=2)
#apply(dat, 2, mean)
#apply(dat, 2, var)

# standardizzazione delle variabili
x.scaled<-scale(dat)

# k-means con k=2
set.seed(4)
km.out<-kmeans(x.scaled, centers=2)
km.out$cluster   # vettore del clustering
km.out$centers
km.out$withinss
plot(dat, col=(km.out$cluster+1), pch=19)

# k-means con k=3
set.seed(6)
km.out<-kmeans(x.scaled, 3)  
km.out$tot.withinss
plot(dat, col=(km.out$cluster), pch=19)
# nstart>1 --> diverse inizializzazioni per l'algoritmo
km.out2<-kmeans(x.scaled, 3, nstart = 20)  
plot(dat, col=(km.out2$cluster), pch=19)
km.out2$tot.withinss

### Esercizio con i dati "sanita.txt"
sanita<-read.table("sanita.txt", header =TRUE, row.names=1)

# 1. standardizzare le variabili ed ottenere il risultato di k-means con K=3 
# 2. estrarre le regioni appartenenti a ciascun gruppo 

set.seed(4)
sanita_st<-scale(sanita)
K<-3   # numero di clusters
res<-kmeans(sanita_st, K, nstart=20)
which(res$cluster==1)

rownames(sanita_st)[which(res$cluster==1)] #regioni in Cluster 1
rownames(sanita_st)[which(res$cluster==2)] #regioni in Cluster 2
rownames(sanita_st)[which(res$cluster==3)] #regioni in Cluster 3

## 06/05/26
sanita<-read.table("sanita.txt", header = TRUE, row.names = 1)
sanita_st<-scale(sanita)

# selezione del numero di gruppi per K-means --> indici di validità 
## Esempio: CH index

#install.packages("fpc")
library(fpc)
library(ggplot2)
rangeK<- 2:8
dim(sanita_st)

res<-sapply(rangeK, function(nc, data)
calinhara(data, kmeans(data, nc, nstart=20)$cluster), 
data=sanita_st)

ggplot(data=data.frame(x=2:(length(res)+1), y=res), 
mapping=aes(x=x, y=y)) +
  geom_point() +
  geom_line() +
  xlab("k") +
  ylab("CH index")

# in alternativa, si può utilizzare un ciclo for()
ch.v<-c()
for(i in 1:length(rangeK)){
 km<-kmeans(sanita_st, rangeK[i], 
           nstart=20)$cluster
 ch.v[i]<-calinhara(sanita_st, km) 
}
plot(2:8, ch.v, type="o")