Wine Clustering

In questo caso noi stiamo lavorando con dati multivariati… Vuoldire che noi abbiamo una matrice di dati dove sulle righe ci sono le unità che sono le componenti chimiche e sulle colonne ci sono le variabil, Noi vogliamo trovare dei gruppi, delle classi di individui simili fra loro, quindi vogliamo vedere se nei nostri dati sono presentidei gruppi. La cluster analysis viene chiamata anche analisi dei gruppi Stiamo lavorando con variabili quantitative numeriche,fatta eccezione per la variabile “type”(variabile qualitativa). Abbiamo un dataset analitico di 178 vini piemontesi appartenenti a tre diversi vitigni (Barbera, Grignolino, Barolo) questo dataset contiene 28 variabili di cui 27 numeriche che rappresentano la composizione chimica e una variabile categorica “type” che rappresenta i diversi tipi di vino (che sono 3) 1 = Barolo 2 = Grignolino 3 = Barbera Lo scopo principale del nostro lavoro è stato dimostrare che queste 27 variabili numeriche chimiche,sono in grado di distinguere diverse “tipologie” o marche di vino.

#Importiamo le librerie

library(gridExtra)
library(tidytext)
library(tidyverse)
library(psych)
library(corrr)
library(tidymodels)
library(GGally)  
library(janitor)             
library(kohonen)
library(factoextra)
library(cluster)
library(FactoMineR)
library(Factoshiny)
library(readr)
library(stringr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(skimr)
library(gt)
library(VIM)
library(corrplot)
library(NbClust)

#lettura e visualizzazione dati

View(wine)
skim(wine) # Skim è una funziona che si trova all'interno del pacchetto Skimr, molto simile al summary

── Data Summary ────────────────────────
                           Values
Name                       wine  
Number of rows             178   
Number of columns          28    
_______________________          
Column type frequency:           
  numeric                  28    
________________________         
Group variables            None  

creiaamo nuovo dataset e togliamo la prima colonna (type)

wine1 <- wine[2:28] 

#Check dei valori mancanti

Grafico che ci mostra che effettivamente non ci sono valori mancanti (abbiamo usato una funsione che si chiama aggr che si trova nel pacchetto VIM)

aggr(wine1)

#Costruiamo una tabella di frequenza Questa tabella di frequenza è stata presa in considerazione dalla variabile Type del dataset, ci servirà in un secondo momento per confrontare i valori dei cluster (abbiamo usato la funzione tabyl che si trova nel pacchetto janitor e gt che si trova nel pacchetto gt)

tab_freq <- tabyl(wine, Type) %>% 
  adorn_totals() %>% 
  adorn_pct_formatting() %>% 
  mutate(Vini = recode(Type,
                               '1'= 'Barolo',
                               '2' = 'Grignolino',
                               '3'= 'Barbera')) %>% 
  relocate(Vini, .after = Type) %>% 
  gt()
tab_freq

Type	Vini	n	percent
1	Barolo	59	33.1%
2	Grignolino	71	39.9%
3	Barbera	48	27.0%
Total	Total	178	100.0%

#Check outlier e visualizzazione dataset

Quando facciamo l’analisi dei gruppi le variabili devono essere misurate sulla stessa scala, quindi dobbiamo standardizzare le variabili se sono molto diverse dalle altre…

notiamo subito che che ci sono dei valori che si discostano molto dagli altri

wine1 %>%
  gather(Attributes, values) %>%
  ggplot(aes(x=reorder(Attributes, values, FUN=median), y=values, fill=Attributes)) +
  geom_boxplot(show.legend=FALSE) +
  labs(title="Wines Attributes - Boxplots") +
  theme_bw() +
  theme(axis.title.y=element_blank(),
        axis.title.x=element_blank()) +
  ylim(0, 800) +
  coord_flip()

Avvertimento: Removed 239 rows containing non-finite values (`stat_boxplot()`).

Con la funzione scale() per ciascun dato viene calcolata la deviata normale standardizzata. In pratica questa funzione prima calcola per i dati di ciascuna colonna/variabile la media e la deviazione standard, poi calcola per ciascuno dato x la corrispondente deviata normale standardizzata

z = (x – media) / deviazione standard

Quando usiamo scale, nel caso in cui volessimo un dataframe dobbiamo specificarlo, altrimenti ritorna una matrice di dati

wine2_s <- as.data.frame(scale(as.data.frame(wine1)))

BoxPlot con dataset standardizzato

wine2_s %>%
  gather(Attributes, values) %>%
  ggplot(aes(x=reorder(Attributes, values, FUN=median), y=values, fill=Attributes)) +
  geom_boxplot(show.legend=FALSE) +
  labs(title="Wines Attributes - Boxplots") +
  theme_bw() +
  theme(axis.title.y=element_blank(),
        axis.title.x=element_blank()) +
  ylim(0, 10) +
  coord_flip()

Avvertimento: Removed 2529 rows containing non-finite values (`stat_boxplot()`).

Vogliamo cercare di fare una selezione di variabili attraverso la visualizzazzione dei grafici

Partendo dal presupposto che abbiamo 27 variabili, facendo un plot di tutte le variabili sarebbe dificile riuscire ad individuare i cluster per le variabili, quindi abbiamo pensato di fare 3 subplot con

** Impostiamo numeri random e utiliazziamo la funzione kmeans** L’algoritmo K-means è un algoritmo di analisi dei gruppi partizionale che permette di suddividere un insieme di oggetti in k gruppi sulla base dei loro attributi

set.seed(123)
wine_dens <- kmeans(wine2_s, centers=3)
summary(wine_dens)

             Length Class  Mode   
cluster      178    -none- numeric
centers       81    -none- numeric
totss          1    -none- numeric
withinss       3    -none- numeric
tot.withinss   1    -none- numeric
betweenss      1    -none- numeric
size           3    -none- numeric
iter           1    -none- numeric
ifault         1    -none- numeric

**Stiamo cercando di fare una selezione delle features attraverso i grafici, e abbiamo provato a farlo anche attraverso la correlazione, siamo certi esistono metodi statistici per la selezione delle features, (abbiamo visto un pacchetto che si chiama caret) ma in fine abbiamo deciso di tenere in considerazione tutte le variabili eccetto type

ggcorr(wine2_s, low = "navy", high = "darkred")

**Con questo istogramma vediamo quante volte si ripetono valori delle variabili

wine2_s %>% 
  gather(key = value_groups, value = value) %>% 
  ggplot(aes(x=value)) +
  geom_histogram() +
  facet_wrap(.~value_groups, scales = "free")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

** Con un qq plot vediamo le distribuzioni delle variabili chimiche e notiamo che sono poche le variabili che sono distribuite come una normale**

wine2_s %>% 
  gather(key = value_groups, value = value) %>% 
  ggplot(aes(sample=value, colour=value_groups)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(.~value_groups, scales = "free_y")

dobbiamo definire una matrice di dissimilarità .. andiamo a vedere quanto sono simili o diversi partendo dalla nostra matrice di dati. Per ciascuna coppia di individui, andiamo a vedere per ciascuna variabile quanto questi due individui sono diversi o simili, a seconda del valore R ( il peso che noi diamo alle distanze) otteniamo delle metriche diverse, la più famosa è la norm euclidea dove R=2

Defianiamo le distanze (euclidea,manhattan,minkowski)

mat <- hclust(dist(wine2_s, method = "euclidean"), method="ward.D") # matrice delle distanze euclidee
mat1 <- hclust(dist(wine2_s, method = "manhattan"), method="ward.D")# distanza di manhattan
mat2 <- hclust(dist(wine2_s, method = "minkowski"), method="ward.D")# distanza di minkowski

**Queste funzioni ci consentono di capire la quantità ottimale dei cluster che possiamo fare con il dataset

fviz_nbclust(wine2_s, FUNcluster = kmeans, nstart=50, "silhouette") # Il coefficente di Silhouette misura quanto e' simile un oggetto x con gli altri oggetti nel proprio cluster contro quelli negli altri cluster. I valori piu' vicini a 1 indicano un cluster migliore

fviz_nbclust(wine2_s, FUNcluster = kmeans, nstart=80, "wss") +
  geom_vline(xintercept = 3, linetype=2)

fviz_nbclust(wine2_s, FUNcluster = kmeans, nstart=30, "gap_stat", nboot = 50)

Clustering k = 1,2,..., K.max (= 10): .. done
Bootstrapping, b = 1,2,..., B (= 50)  [one "." per sample]:
.................................................. 50 

** Questa funzione è un test che genera in automatico il numero di cluster che sarebbe meglio tirar fuoriin questo caso il numero delle classsi migliori è 3**

x <- NbClust(wine2_s, distance = "euclidean", min.nc = 2, max.nc = 10, method = "ward.D2")

*** : The Hubert index is a graphical method of determining the number of clusters.
                In the plot of Hubert index, we seek a significant knee that corresponds to a 
                significant increase of the value of the measure i.e the significant peak in Hubert
                index second differences plot. 
 

*** : The D index is a graphical method of determining the number of clusters. 
                In the plot of D index, we seek a significant knee (the significant peak in Dindex
                second differences plot) that corresponds to a significant increase of the value of
                the measure. 
 
******************************************************************* 
* Among all indices:                                                
* 4 proposed 2 as the best number of clusters 
* 16 proposed 3 as the best number of clusters 
* 1 proposed 8 as the best number of clusters 
* 1 proposed 10 as the best number of clusters 

                   ***** Conclusion *****                            
 
* According to the majority rule, the best number of clusters is  3 
 
 
******************************************************************* 

** Visualizziamo tutte le proposte **

x$All.index

       KL      CH Hartigan    CCC     Scott      Marriot   TrCovW   TraceW Friedman  Rubin Cindex     DB Silhouette   Duda Pseudot2  Beale Ratkowsky      Ball Ptbiserial    Frey McClain   Dunn Hubert SDindex Dindex   SDbw
2  1.3943 43.3719  32.9472 0.2449  382.9864 1.970498e+52 39291.48 3834.146  31.5640 1.2464 0.3400 1.7959     0.2039 0.7883  33.3041 5.0571    0.2592 1917.0732     0.4622  0.5073  0.5711 0.2152 0.0008  0.8833 4.5062 0.8437
3  5.6787 41.9793   8.6883 2.9309  696.9339 7.599268e+51 23170.59 3229.570  44.8292 1.4798 0.2916 1.8784     0.1984 0.8841   7.9972 2.4483    0.3079 1076.5235     0.5233  0.5038  1.5205 0.2102 0.0008  0.8438 4.1069 0.7603
4  0.8909 32.0873   8.5378 1.6012  859.8732 5.408798e+51 20509.53 3076.814  49.6488 1.5532 0.3589 2.2321     0.1773 0.8529   8.6209 3.2085    0.2840  769.2036     0.5175  0.4614  1.9870 0.2543 0.0009  0.9683 4.0169 0.7608
5  1.5750 27.2234   5.9738 1.3038 1034.9995 3.159641e+51 18163.48 2932.903  54.1784 1.6294 0.3443 2.3906     0.1486 0.8951   7.1500 2.1889    0.2661  586.5807     0.5085  0.7322  2.4358 0.2543 0.0010  0.9575 3.9167 0.7305
6  0.9464 23.5884   5.8566 0.7741 1144.9602 2.453081e+51 16677.85 2835.009  56.3149 1.6857 0.3336 2.4890     0.1099 0.8650   5.6198 2.8829    0.2520  472.5016     0.4816 -0.0319  3.0509 0.2543 0.0011  1.0231 3.8527 0.8097
7  0.9142 21.1786   6.0199 0.5208 1251.0749 1.839504e+51 15771.56 2741.656  59.3104 1.7431 0.3642 2.3391     0.1116 0.8032   5.6347 4.4563    0.2403  391.6651     0.4956 -0.0124  3.1347 0.2790 0.0011  0.9806 3.8049 0.7880
8  1.0787 19.5372   5.5609 0.5966 1389.9304 1.101278e+51 14749.02 2648.421  61.3904 1.8045 0.3707 2.1743     0.1137 0.8478   5.9227 3.3064    0.2308  331.0526     0.5034  0.2329  3.1914 0.2868 0.0012  0.9765 3.7456 0.7288
9  1.0894 18.2414   5.1187 0.7401 1492.9426 7.813888e+50 13587.34 2564.531  64.0273 1.8635 0.3614 2.1318     0.1118 0.8810   6.6175 2.5121    0.2222  284.9479     0.5024  1.2192  3.4373 0.2868 0.0012  0.9872 3.6895 0.7170
10 0.9953 17.1722   5.0641 0.8378 1584.5417 5.766250e+50 12510.82 2489.140  66.2517 1.9199 0.3558 2.1039     0.1048 0.7946   5.1698 4.6727    0.2145  248.9140     0.4698  0.2131  4.1187 0.2868 0.0012  0.9795 3.6321 0.6917

Sulla base delle distanze che abbiamo andiamo a considerare ● Single linkage —> Nearest Neighbor (il vicino più vicino) cerca di mettere ● Complete linkage —> il vicino più lontano cerca di mettere insieme gruppi tra di loro che sono vicini comparando le distanze fra di loro ● Average linkage —> il linkage medio fa una media dei due ● Ward’s minimum variance linkage —> il word ion genere è quello che funziona meglio…word minimizza la varianza within…perché il metodo di word non si basa sulla distanza ma ha l’idea di minimizzare l’inerzia ovvero la varianza….quindi lui tende a far si che la varianza whitin non sia troppo grande

Definiamo le variabili per la visualizzazione dei dendogrammi

d_eucli <- dist(wine2_s, method = "euclidean") # distanza euclide
wine_hc <- hclust(d_eucli, method = "ward.D2")
fviz_dend(wine_hc, k = 3, show_labels = F, rect = T) #dendogramma migliore

{
vini_complete <- hclust(d_eucli, method = "complete")   ## cluster with "complete" linkage
vini_single <- hclust(d_eucli, method = "single")     ## cluster with "single" linkage
vini_average <- hclust(d_eucli, method = "average")   ## cluster with "average" linkage
vini_ward <- hclust(d_eucli, method = "ward.D2")      ## cluster with "ward" linkage
}

Facciamo una comparazione dei dendogrammi per vedere quale ci sembra il migliore

par(mfrow = c(2, 2))
plot(vini_complete, cex = 0.8, hang = -1, main = "Hierarchical Clustering with Complete Linkage")
plot(vini_single, cex = 0.8, hang = -1, main = "Hierarchical Clustering with Single Linkage")

plot(vini_average, cex = 0.8, hang = -1, main = "Hierarchical Clustering with Average Linkage")
plot(vini_ward, cex = 0.8, hang = -1, main = "Hierarchical Clustering with Ward Linkage")

Abbiamo optato per il dendogramma migliore,abbiamo evidenziato il dendogramma con colori diversi e abbiamo tolto i numeri per una visualizzazione migliore

Il centroide è un punto appartenente allo spazio delle features che media le distanze tra tutti i dati appartenenti al cluster ad esso associato. Rappresenta quindi una sorta di baricentro del cluster ed in generale, proprio per le sue caratteristiche, non è uno dei punti del dataset. Non conoscendo le classi presenti nel dataset di ingresso, la prima cosa da fare è decidere il numero di classi, Si scelgono in modo casuale K centroidi appartenenti allo spazio delle features. Si calcola la distanza di ogni punto del dataset rispetto ad ogni centroide Ogni punto del dataset viene associato al cluster collegato al centroide più vicino Si ricalcola la posizione di ogni centroide facendo la media delle posizioni di tutti i punti del cluster associato Si itera dal punto 3 fino a quando non ci sarà più alcun ingresso che cambia di cluster.

set.seed(123)
vini_kmeans <- kmeans(wine2_s, centers = 3)
vini_kmeans

K-means clustering with 3 clusters of sizes 69, 60, 49

Cluster means:
     Alcohol Sugar-free Extract Fixed Acidity Tartaric Acid Malic Acid Uronic Acids           pH        Ash Alcalinity of Ash  Potassium    Calcium   Magnesium  Phosphate   Chloride Total Phenols  Flavanoids Non-flavanoid Phenols Proanthocyanins
1 -0.8856865         -0.4979609    -0.2282955    -0.1499780 -0.3579000   -0.3954178 -0.006752522 -0.5161014         0.1225921 -0.2533765  0.5580202 -0.41519242 -0.4450871  0.1693383    -0.0804181  0.01146806           -0.01373209      0.05540199
2  0.8666244          0.6638956    -0.5146858    -0.5791356 -0.3253961   -0.3561667  0.240038184  0.3905674        -0.6163311  0.2068074 -0.4749177  0.51870720  0.7143103  0.1655889     0.8975314  0.99353094           -0.56798207      0.54674519
3  0.1860184         -0.1117231     0.9517048     0.9203392  0.9024258    0.9929352 -0.284415653  0.2485092         0.5820616  0.1035619 -0.2042516 -0.05049296 -0.2479104 -0.4412179    -0.9857762 -1.23271740            0.71482528     -0.74749896
  Color Intensity        Hue OD280/OD315 of Diluted Wines OD280/OD315 of Flavanoids    Glycerol 2-3-Butanediol Total Nitrogen    Proline    Methanol
1      -0.8730101  0.4230425                    0.2317587                 0.3918922 -0.59945726     -0.5310122     -0.1339752 -0.7333939 -0.18931901
2       0.1989588  0.4836584                    0.7934116                 0.5529190  0.66943876      0.1927984      0.5168134  1.1528997  0.03215541
3       0.9857177 -1.1879477                   -1.2978785                -1.2288919  0.02441276      0.5116722     -0.4441737 -0.3789756  0.22721810

Clustering vector:
  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 
  2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   1   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   1   1   1   1   1 
 65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 
  1   1   1   1   1   1   1   1   1   2   1   1   1   1   1   1   1   1   1   3   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   2   1   1   1   1   1   1 
129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 
  1   1   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3 

Within cluster sum of squares by cluster:
[1] 1373.5013  919.6663  893.0953
 (between_SS / total_SS =  33.3 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"    "size"         "iter"         "ifault"      

table(wine$Type,vini_kmeans$cluster)

   
     1  2  3
  1  1 58  0
  2 68  2  1
  3  0  0 48

fviz_cluster(vini_kmeans, data = wine2_s,
             palette = c("#2E9FDF", "#00AFBB", "#E7B800"), 
             geom = "point",
             ellipse.type = "convex", 
             ggtheme = theme_bw()
)

print(vini_kmeans)

K-means clustering with 3 clusters of sizes 69, 60, 49

Cluster means:
     Alcohol Sugar-free Extract Fixed Acidity Tartaric Acid Malic Acid Uronic Acids           pH        Ash Alcalinity of Ash  Potassium    Calcium   Magnesium  Phosphate   Chloride Total Phenols  Flavanoids Non-flavanoid Phenols Proanthocyanins
1 -0.8856865         -0.4979609    -0.2282955    -0.1499780 -0.3579000   -0.3954178 -0.006752522 -0.5161014         0.1225921 -0.2533765  0.5580202 -0.41519242 -0.4450871  0.1693383    -0.0804181  0.01146806           -0.01373209      0.05540199
2  0.8666244          0.6638956    -0.5146858    -0.5791356 -0.3253961   -0.3561667  0.240038184  0.3905674        -0.6163311  0.2068074 -0.4749177  0.51870720  0.7143103  0.1655889     0.8975314  0.99353094           -0.56798207      0.54674519
3  0.1860184         -0.1117231     0.9517048     0.9203392  0.9024258    0.9929352 -0.284415653  0.2485092         0.5820616  0.1035619 -0.2042516 -0.05049296 -0.2479104 -0.4412179    -0.9857762 -1.23271740            0.71482528     -0.74749896
  Color Intensity        Hue OD280/OD315 of Diluted Wines OD280/OD315 of Flavanoids    Glycerol 2-3-Butanediol Total Nitrogen    Proline    Methanol
1      -0.8730101  0.4230425                    0.2317587                 0.3918922 -0.59945726     -0.5310122     -0.1339752 -0.7333939 -0.18931901
2       0.1989588  0.4836584                    0.7934116                 0.5529190  0.66943876      0.1927984      0.5168134  1.1528997  0.03215541
3       0.9857177 -1.1879477                   -1.2978785                -1.2288919  0.02441276      0.5116722     -0.4441737 -0.3789756  0.22721810

Clustering vector:
  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 
  2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   1   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   1   1   1   1   1 
 65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 
  1   1   1   1   1   1   1   1   1   2   1   1   1   1   1   1   1   1   1   3   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   2   1   1   1   1   1   1 
129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 
  1   1   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3 

Within cluster sum of squares by cluster:
[1] 1373.5013  919.6663  893.0953
 (between_SS / total_SS =  33.3 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"    "size"         "iter"         "ifault"      

vini_kmeans

K-means clustering with 3 clusters of sizes 69, 60, 49

Cluster means:
     Alcohol Sugar-free Extract Fixed Acidity Tartaric Acid Malic Acid Uronic Acids           pH        Ash Alcalinity of Ash  Potassium    Calcium   Magnesium  Phosphate   Chloride Total Phenols  Flavanoids Non-flavanoid Phenols Proanthocyanins
1 -0.8856865         -0.4979609    -0.2282955    -0.1499780 -0.3579000   -0.3954178 -0.006752522 -0.5161014         0.1225921 -0.2533765  0.5580202 -0.41519242 -0.4450871  0.1693383    -0.0804181  0.01146806           -0.01373209      0.05540199
2  0.8666244          0.6638956    -0.5146858    -0.5791356 -0.3253961   -0.3561667  0.240038184  0.3905674        -0.6163311  0.2068074 -0.4749177  0.51870720  0.7143103  0.1655889     0.8975314  0.99353094           -0.56798207      0.54674519
3  0.1860184         -0.1117231     0.9517048     0.9203392  0.9024258    0.9929352 -0.284415653  0.2485092         0.5820616  0.1035619 -0.2042516 -0.05049296 -0.2479104 -0.4412179    -0.9857762 -1.23271740            0.71482528     -0.74749896
  Color Intensity        Hue OD280/OD315 of Diluted Wines OD280/OD315 of Flavanoids    Glycerol 2-3-Butanediol Total Nitrogen    Proline    Methanol
1      -0.8730101  0.4230425                    0.2317587                 0.3918922 -0.59945726     -0.5310122     -0.1339752 -0.7333939 -0.18931901
2       0.1989588  0.4836584                    0.7934116                 0.5529190  0.66943876      0.1927984      0.5168134  1.1528997  0.03215541
3       0.9857177 -1.1879477                   -1.2978785                -1.2288919  0.02441276      0.5116722     -0.4441737 -0.3789756  0.22721810

Clustering vector:
  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 
  2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   1   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   1   1   1   1   1 
 65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 
  1   1   1   1   1   1   1   1   1   2   1   1   1   1   1   1   1   1   1   3   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   2   1   1   1   1   1   1 
129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 
  1   1   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3 

Within cluster sum of squares by cluster:
[1] 1373.5013  919.6663  893.0953
 (between_SS / total_SS =  33.3 %)

Available components:

[1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss" "betweenss"    "size"         "iter"         "ifault"      

bar <- subset(wine, Type == "1")
bar2 <- subset(wine, Type == "2")
bar3 <- subset(wine, Type == "3")
mean(bar$Alcohol)

[1] 13.74475

mean(bar2$Alcohol)

[1] 12.27873

mean(bar3$Alcohol)

[1] 13.15375

bar <- subset(wine, Type == "1")
bar2 <- subset(wine, Type == "2")
bar3 <- subset(wine, Type == "3")
mean(bar$Proline)

[1] 1115.712

mean(bar2$Proline)

[1] 519.507

mean(bar3$Proline)

[1] 629.8958

cluster basato su modelli identifichiamo la classe, ovvero il tipo dei vini:

class <- factor(wine$Type, levels = 1:3, labels = c("Barolo", "Grignolino", "Barbera"))
table(class)

class
    Barolo Grignolino    Barbera 
        59         71         48 

Selezione delle variabili con clustvarsel

library(clustvarsel)#La metodologia consente di trovare il sottoinsieme (localmente) ottimale di variabili in un set di dati che hanno informazioni di gruppo/cluster.
out <- clustvarsel(wine1)      # utilizzo la funzione clustvarsel per selezionare le varibili ottimali

iter 1
+ adding step

iter 2
+ adding step

iter 3 
+ adding step
- removing step

iter 4 
+ adding step
- removing step

iter 5 
+ adding step
- removing step

iter 6 
+ adding step
- removing step

iter 7 
+ adding step
- removing step

iter 8 
+ adding step
- removing step

iter 9 
+ adding step
- removing step

final iter
* fitting model on selected subset

out                              # ai fini del clustering, in modo da utilizzarle nei grafici

------------------------------------------------------ 
Variable selection for Gaussian model-based clustering
Stepwise (forward/backward) greedy search
------------------------------------------------------ 

Selected subset: Chloride, Malic Acid, Flavanoids, Color Intensity, Proline, Uronic Acids, OD280/OD315 of Flavanoids

out$subset

                 Chloride                Malic Acid                Flavanoids           Color Intensity                   Proline              Uronic Acids OD280/OD315 of Flavanoids 
                       14                         5                        16                        19                        26                         6                        22 

selected.var <- wine1[out$subset]
selected.var

# Le variabili selezionate sono: Chloride, Malic Acid, Flavanoids, Color Intensity, Proline, Uronic Acids,
#OD280/OD315 of Flavanoids
#la prolina indica la secchezza del vino, l'acido malico è l'acido più diffuso nel regno vegetale
#il gruppo dei flavonoidi ha una concentrazione più elevata nei vini rossi, e svolgono un azione protettiva contro le radiazioni ultraviolette

Analisi esplorativa delle variabili selezionate

selezione delle variabili con la libreria “BORUTA”,Boruta è un algoritmo di classificazione e selezione delle caratteristiche basato su un algoritmo di foreste casuali. vantaggio di Boruta è che decide chiaramente se una variabile è importante o meno e aiuta a selezionare variabili statisticamente significative.

selezione delle variabili con il metodo SelvarClustLasso, questa funzione implementa la selezione delle variabili nel clustering basato su modelli

library(sparcl) 
library(vscc) 
library(SelvarMix) 
require(Rmixmod)
require(glasso)
lasso <- SelvarClustLasso(x=wine1, nbcluster=1:3,criterio="ICL")

variable  ranking
SRUW selection with ICL criterion
model selection  with ICL criterion

summary(lasso)

Criterion: ICL 
Criterion value: -10843.71 
Number of clusters: 3 
Gaussian mixture model: Gaussian_pk_Lk_Ck 
Regression covariance model: LC 
Independent covariance model:  
The SRUW model:
 S: 1 6 9 10 12 14 16 18 19 20 21 22 23 25 26 27 
 R: 9 10 12 16 19 20 23 
 U: 2 3 4 5 7 8 11 13 15 17 24 
 W: 

lasso$S #variabili accettate

 [1]  1  6  9 10 12 14 16 18 19 20 21 22 23 25 26 27

lasso$R #variabili per regressione

[1]  9 10 12 16 19 20 23

lasso$U #variabili ridondanti

 [1]  2  3  4  5  7  8 11 13 15 17 24

lasso$nbcluster

[1] 3

table(lasso$partition,wine[,1])

   
     1  2  3
  1 59 71 48

** MODEL - BASED CLUSTERING applicato a tutto il dataset wine **

Il miglior modello è EVI, con tre cluster, EVI ha distrubuzione diagonale, volume uguale e forma variabile (aggiungere foto)

rappresentazione dei risultati

rappresentazione dei valori Bic utilizzati per la scelta del numero dei cluster

plot(mod, what = "BIC", ylim = range(mod$BIC[,-(1:2)], na.rm = TRUE), legendArgs = list(x = "bottomleft"))

** rappresentazione dei risultati con la libreria factoextraClassificazione: grafico che mostra il raggruppamento**

library(factoextra)
fviz_mclust(mod, "classification", geom = "point", 
            pointsize = 1.5, palette = "jco")

Classificazione dell’incertezza

fviz_mclust(mod, "uncertainty", palette = "jco")

** Mclust con le variabili selezionate con il metodo di clustvarsel**

fviz_mclust(mclust.selected, "classification", geom = "point", 
            pointsize = 1.5, palette = "set2")

fviz_mclust(mclust.selected, "uncertainty", palette = "set2")

** mclust con le variabile selezionate ci dice cge il miglior modello è EVI, ma con 7 cluster**

ora rifacciamo la stessa funzione, ma impostando il numero di cluster a 3

plot(mclust.selected1, what = "classification")

mclustDA Analisi discriminante basata sulla modellazione di miscele finite gaussiane

wineMclustDA <- MclustDA(selected.var, class, modelType = "EDDA", modelNames = "EEE")

fitting ...

  |                                                                                                                                                                                                                                                            
  |                                                                                                                                                                                                                                                      |   0%
  |                                                                                                                                                                                                                                                            
  |======================================================================================================================================================================================================================================================| 100%

summary(wineMclustDA, parameters = TRUE)

------------------------------------------------ 
Gaussian finite mixture model for classification 
------------------------------------------------ 

EDDA model summary: 

            
Classes       n     % Model G
  Barolo     59 33.15   EEE 1
  Grignolino 71 39.89   EEE 1
  Barbera    48 26.97   EEE 1

Class prior probabilities:
    Barolo Grignolino    Barbera 
 0.3314607  0.3988764  0.2696629 

Class = Barolo

Means:
                                  [,1]
Chloride                    71.7118644
Malic Acid                   2.0106780
Flavanoids                   2.9823729
Color Intensity              5.5283051
Proline                   1115.7118644
Uronic Acids                 0.8110169
OD280/OD315 of Flavanoids    3.4484746

Variances:
[,,1]
                              Chloride    Malic Acid   Flavanoids Color Intensity      Proline Uronic Acids OD280/OD315 of Flavanoids
Chloride                  2309.3076627  -5.292799727  2.892527934      4.13800916   420.801184 -0.868100227              -3.584535341
Malic Acid                  -5.2927997   0.872588143 -0.009250663     -0.25433997   -32.501326  0.039696376               0.128537935
Flavanoids                   2.8925279  -0.009250663  0.270077612      0.28184664     3.367962  0.005806350               0.074876906
Color Intensity              4.1380092  -0.254339969  0.281846643      2.24641309    66.946034 -0.010459694              -0.143612129
Proline                    420.8011842 -32.501326194  3.367961843     66.94603423 29206.990603  3.068294871             -18.334531286
Uronic Acids                -0.8681002   0.039696376  0.005806350     -0.01045969     3.068295  0.038401609               0.008848205
OD280/OD315 of Flavanoids   -3.5845353   0.128537935  0.074876906     -0.14361213   -18.334531  0.008848205               0.338611842

Class = Grignolino

Means:
                                 [,1]
Chloride                   71.6901408
Malic Acid                  1.9326761
Flavanoids                  2.0808451
Color Intensity             3.0866197
Proline                   519.5070423
Uronic Acids                0.8253521
OD280/OD315 of Flavanoids   3.3316901

Variances:
[,,1]
                              Chloride    Malic Acid   Flavanoids Color Intensity      Proline Uronic Acids OD280/OD315 of Flavanoids
Chloride                  2309.3076627  -5.292799727  2.892527934      4.13800916   420.801184 -0.868100227              -3.584535341
Malic Acid                  -5.2927997   0.872588143 -0.009250663     -0.25433997   -32.501326  0.039696376               0.128537935
Flavanoids                   2.8925279  -0.009250663  0.270077612      0.28184664     3.367962  0.005806350               0.074876906
Color Intensity              4.1380092  -0.254339969  0.281846643      2.24641309    66.946034 -0.010459694              -0.143612129
Proline                    420.8011842 -32.501326194  3.367961843     66.94603423 29206.990603  3.068294871             -18.334531286
Uronic Acids                -0.8681002   0.039696376  0.005806350     -0.01045969     3.068295  0.038401609               0.008848205
OD280/OD315 of Flavanoids   -3.5845353   0.128537935  0.074876906     -0.14361213   -18.334531  0.008848205               0.338611842

Class = Barbera

Means:
                                 [,1]
Chloride                   41.2708333
Malic Acid                  3.3337500
Flavanoids                  0.7814583
Color Intensity             7.3962500
Proline                   629.8958333
Uronic Acids                1.1743750
OD280/OD315 of Flavanoids   1.8860417

Variances:
[,,1]
                              Chloride    Malic Acid   Flavanoids Color Intensity      Proline Uronic Acids OD280/OD315 of Flavanoids
Chloride                  2309.3076627  -5.292799727  2.892527934      4.13800916   420.801184 -0.868100227              -3.584535341
Malic Acid                  -5.2927997   0.872588143 -0.009250663     -0.25433997   -32.501326  0.039696376               0.128537935
Flavanoids                   2.8925279  -0.009250663  0.270077612      0.28184664     3.367962  0.005806350               0.074876906
Color Intensity              4.1380092  -0.254339969  0.281846643      2.24641309    66.946034 -0.010459694              -0.143612129
Proline                    420.8011842 -32.501326194  3.367961843     66.94603423 29206.990603  3.068294871             -18.334531286
Uronic Acids                -0.8681002   0.039696376  0.005806350     -0.01045969     3.068295  0.038401609               0.008848205
OD280/OD315 of Flavanoids   -3.5845353   0.128537935  0.074876906     -0.14361213   -18.334531  0.008848205               0.338611842

Training confusion matrix:
            Predicted
Class        Barolo Grignolino Barbera
  Barolo         56          3       0
  Grignolino      1         69       1
  Barbera         0          1      47
Classification error = 0.0337 
Brier score          = 0.0279 

summary(wineMclustDA, newdata = selected.var, newclass = class)

------------------------------------------------ 
Gaussian finite mixture model for classification 
------------------------------------------------ 

EDDA model summary: 

            
Classes       n     % Model G
  Barolo     59 33.15   EEE 1
  Grignolino 71 39.89   EEE 1
  Barbera    48 26.97   EEE 1

Training confusion matrix:
            Predicted
Class        Barolo Grignolino Barbera
  Barolo         56          3       0
  Grignolino      1         69       1
  Barbera         0          1      47
Classification error = 0.0337 
Brier score          = 0.0279 

Test confusion matrix:
            Predicted
Class        Barolo Grignolino Barbera
  Barolo         56          3       0
  Grignolino      1         69       1
  Barbera         0          1      47
Classification error = 0.0337 
Brier score          = 0.0279 

plot(wineMclustDA)

Model-based discriminant analysis plots: 

1: scatterplot
2: classification
3: train&test
4: error

1

Model-based discriminant analysis plots: 

1: scatterplot
2: classification
3: train&test
4: error

2

Model-based discriminant analysis plots: 

1: scatterplot
2: classification
3: train&test
4: error

3

Model-based discriminant analysis plots: 

1: scatterplot
2: classification
3: train&test
4: error

4

Model-based discriminant analysis plots: 

1: scatterplot
2: classification
3: train&test
4: error

0

NA
NA

visualizziamo i cluster su un biplot PCA

wine_pca <- prcomp(X, scale. = TRUE) #Esegue un'analisi dei componenti principali sulla matrice di dati data e restituisce i risultati come oggetto di classe prcomp.

library(ggbiplot)
mc_clusters <- factor(mod$classification)
ggbiplot(wine_pca, groups = mc_clusters) +
  labs(color = "mclust clusters") +
  scale_color_brewer(palette = "Set2") +
  theme(aspect.ratio = 0.8, legend.position = "top")

** analisi delle componenti principaòi con il dataset ridotto**

wine_pca1 <- prcomp(selected.var, scale. = TRUE)
mc_clusters1 <- factor(mclust.selected1$classification)
ggbiplot(wine_pca1, groups = mc_clusters1) +
  labs(color = "mclust clusters") +
  scale_color_brewer(palette = "Set2") +
  theme(aspect.ratio = 0.8, legend.position = "top")