-
Notifications
You must be signed in to change notification settings - Fork 0
/
Food_production_analysis.r
185 lines (138 loc) · 5.11 KB
/
Food_production_analysis.r
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
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
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
179
180
181
182
183
184
185
# Initialisation ----
# Chargemeent des données
donnees <- fread("Food_Production.csv")
donnees <- as.data.frame(donnees)
# Nomme les colonnes par le nom du produit alimentaire correspondant
names <- donnees$`Food product`
rownames(donnees) <- names
# Elimination des produits qui ont une valeur négative pour la variable Lande use change
donnees <- donnees %>% filter(`Land use change`>=0)
donnees <- donnees[,c("Land use change", "Animal Feed","Farm","Processing", "Transport","Packging","Retail","Total_emissions")]
# Satistiques descriptives ----
set.caption('Statistiques descriptives') #Légende du tableau
panderOptions('table.split.table', 300) # Pour éviter de spliter les tableaux
pander(apply(donnees,2,stat.desc), digits = 2, cex = 0.7)
boxplot(donnees[,1:7], ylab = "Kg d'équivalent CO2 par kg")
pairs(donnees)
# Corrélations entre les variables
set.caption('Matrice des corrélations')
pander(cor(donnees), digits = 2, cex = 0.1)
panderOptions('table.split.table', 300)
Data_ACP <- subset(donnees, select = -c(Total_emissions))
# ACP ----
# Ajout d'une variable catégorielle
x=c()
for (i in 1:29) {
x[i] <- "No"
}
y= c()
for (i in 1:10) {
y[i] <- "Yes"
}
Product_type <- c(x,y)
Data_ACP <- cbind(Data_ACP,Product_type)
res <- PCA(Data_ACP[,-8], graph=FALSE)
# Choix composantes principales
set.caption('Valeurs propres de la matrice de corrélation')
eig <- res$eig
pander(eig)
## Variables ----
# Coordonnées factorielles des variables
coord <- round(res$var$coord[,1:3],3)
set.caption('Coordonnées factorielles des variables')
pander(coord)
# Qualité de représentation des variables
cos2 <- round(res$var$cos2[,1:3],3)
set.caption('Qualité de représentation des variables (Cos^2)')
pander(cos2)
# Contribution des variables
contrib <- round(res$var$contrib[,1:3],3)
set.caption("Contribution des variables à l'élaboration des composantes principales")
pander(contrib)
# Cercles des corrélations
options(ggrepel.max.overlaps = Inf)
fviz_pca_var(res,axes = c(1,2), col.var = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE )
options(ggrepel.max.overlaps = Inf)
fviz_pca_var(res, axes = c(1,3), col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE )
options(ggrepel.max.overlaps = Inf)
fviz_pca_var(res, axes = c(2,3), col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE )
## Individus ----
# Coordonnées factorielles des individus
coord_ind <- round(res$ind$coord[,1:3],3)
set.caption('Coordonnées factorielles des individus')
pander(coord_ind)
# Qualité de représentation des observations
coord_ind <- round(res$ind$cos2[,1:3],3)
set.caption('Qualité de représentation des observations')
pander(coord_ind)
# Contributions des observations
contrib_ind <- round(res$ind$contrib[,1:3],3)
set.caption('Contribution des observations')
pander(contrib_ind)
# Contributions des observations
dist_ind <- round(res$ind$dist,3)
set.caption("Distance à l'origine des observations")
pander(dist_ind)
# Distance à l'origine
dist_ind <- round(res$ind$dist,3)
set.caption("Distance à l'origine des observations")
pander(dist_ind)
# Cartes des individus
options(ggrepel.max.overlaps = Inf)
fviz_pca_ind(res, axes = c(1,2), col.ind = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Évite le chevauchement de texte
)
options(ggrepel.max.overlaps = Inf)
fviz_pca_ind(res, axes = c(1,3), col.ind = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Évite le chevauchement de texte
)
options(ggrepel.max.overlaps = Inf)
fviz_pca_ind(res, axes = c(2,3), col.ind = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Évite le chevauchement de texte
)
options(ggrepel.max.overlaps = Inf)
fviz_pca_ind(res,
geom.ind = "point",
col.ind = Data_ACP$Product_type, # colorer by groups
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
addEllipses = TRUE,
mean.point = TRUE,
legend.title = "Product type"
)
# Clustering ----
# Standardisation des données
donnees <- scale(donnees)
donnees <- donnees[,1:7]
# Distances euclidiennes
dist_donnees <- dist(donnees, method = "euclidean")
# Classification hiérarchique algorithme de Ward
res_clust <- hclust(dist_donnees, method = "ward.D")
plot(res_clust,
xlab = "",
ylab = "Niveau d'agrégation",
cex = 0.7)
rect.hclust(res_clust, k=6, border="red")
rect.hclust(res_clust, k=3, border="green")
rect.hclust(res_clust, k=2, border="blue")
# Barplot
barplot(res_clust$height,
xlab = "Nombre de classes",
names.arg = (nrow(donnees)-1):1,
ylab = "Niveau d'agrégation")
#qualité partition - 6 classes
BSS_W <- sum(tail(res_clust$height,n=(6-1)))
TSS_W <- sum(res_clust$height)
BSS_W/TSS_W*100
#qualité partition - 3 classes
BSS_W <- sum(tail(res_clust$height,n=(3-1)))
TSS_W <- sum(res_clust$height)
BSS_W/TSS_W*100
#qualité partition - 2 classes
BSS_W <- sum(tail(res_clust$height,n=(2-1)))
TSS_W <- sum(res_clust$height)
BSS_W/TSS_W*100