library(mlbench)
data(Glass)
      # funkcja obliczająca błąd
err <- function(y.true, y.pred) { sum(y.pred!=y.true)/length(y.true) }

      # pakiet cluster zawiera implementację kilku algorytmów grupowania
if (! "cluster" %in% row.names(installed.packages()))
  install.packages("cluster")
library(cluster)

      # będziemy grupować na podstawie atrybutów bez uwzględniania kategorii zawartych w danych

      # pam -- algorytm k środków (środki reprezentowane przez medoidy)
      # k=2, bez standaryzacji atrybutów
g.pam2 <- pam(Glass[,-10], 2, stand=F)
      # wykres grup zrzutowanych do dwóch wymiarów (składowych głównych)
clusplot(g.pam2)
      # wykres sylwetkowy (średnia szerokość sylwetki może być miarą jakości grupowania)
plot(silhouette(g.pam2))
      # k=2, ze standaryzacją atrybutów
g.pam2s <- pam(Glass[,-10], 2, stand=T)
clusplot(g.pam2s)
plot(silhouette(g.pam2s))
      # podobnie dla k=3, 5, 7
g.pam3 <- pam(Glass[,-10], 3, stand=F)
clusplot(g.pam3)
plot(silhouette(g.pam3))
g.pam3s <- pam(Glass[,-10], 3, stand=T)
clusplot(g.pam3s)
plot(silhouette(g.pam3s))
g.pam5 <- pam(Glass[,-10], 5, stand=F)
clusplot(g.pam5)
plot(silhouette(g.pam5))
g.pam5s <- pam(Glass[,-10], 5, stand=T)
clusplot(g.pam5s)
plot(silhouette(g.pam5s))
g.pam7 <- pam(Glass[,-10], 7, stand=F)
clusplot(g.pam7)
plot(silhouette(g.pam7))
g.pam7s <- pam(Glass[,-10], 7, stand=T)
clusplot(g.pam7s)
plot(silhouette(g.pam7s))

      # jak podział na grupy ma się do (nieużywanej) kategorii?
table(g.pam2s$clustering, Glass$Type)

      # agnes -- hierarchiczne grupowanie wstępujące
      # bez standaryzacji atrybutów, pojedyncze łączenie
g.as <- agnes(Glass[,-10], stand=F, method="single")
pltree(g.as)
      # ze standaryzacją atrybutów, pojedyncze łączenie
g.ass <- agnes(Glass[,-10], stand=T, method="single")
pltree(g.ass)
      # podobnie dla innych metod łączenia, np.:
g.ac <- agnes(Glass[,-10], stand=F, method="complete")
pltree(g.ac)
g.acs <- agnes(Glass[-10],, stand=T, method="complete")
pltree(g.acs)
g.aw <- agnes(Glass[,-10], stand=F, method="ward")
pltree(g.aw)
g.aws <- agnes(Glass[,-10], stand=T, method="ward")
pltree(g.aws)

      # wytnijmy z drzewa grupowania (np. g.acs) 7 grup
g.aws7 <- cutree(g.aws, 7)
      # jak ten podział ma się do (nieużywanych) kategorii
table(g.aws7, Glass$Type)

      # dla wycinków grupowania hierarchicznego można też sporządzić wykresy sylwetkowe
      # macierz niepodobieństwa
gd <- daisy(Glass[,-10], stand=F)
gds <- daisy(Glass[,-10], stand=T)
plot(silhouette(cutree(g.aw, 2), gd))
plot(silhouette(cutree(g.aw, 3), gd))
plot(silhouette(cutree(g.aw, 5), gd))
plot(silhouette(cutree(g.aw, 7), gd))

plot(silhouette(cutree(g.aws, 2), gds))
plot(silhouette(cutree(g.aws, 3), gds))
plot(silhouette(cutree(g.aws, 5), gds))
plot(silhouette(cutree(g.aws, 7), gds))