R Notebook - PROJET 5 - PEYRONNE

—————————————————————————————-

Rappel du projet

L’objectif de ce projet est de réaliser une étude de marché à l’aide d’informations que l’on peut choisir librement pour trouver des pays dans lesquels il serait possible d’exporter des poulets depuis la France. Ainsi, idéalement, nous recherchons des pays riches et sûrs, voisins de la France dans lesquesl la consommation de protéines animale est élevée mais la production y est faible.

—————————————————————————————-

—————————

1 Importation des données

—————————

1.1 Sources des données

Les données sont extraites du site de la FAO qui fournit de nombreuses informations au niveau mondial sur des thèmes tels que la population, l’agriculture, la sécurité etc.

Il y a des variables imposées :

• différence de population entre une année antérieure (au choix) et l’année courante, exprimée en pourcentage
• proportion de protéines d’origine animale par rapport à la quantité totale de protéines dans la disponibilité alimentaire du pays
• disponibilité alimentaire en protéines par habitant
• disponibilité alimentaire en calories par habitant

Et des variables au choix : - PIB par habitant ( pour avoir une idée de la possibilité des vendre quelquechose ) - Distance depuis la France ( pour évaluer le coût de transport et l’impact écologique ) - Nom des pays au format ISO trois caractères ( pour combiner le tableau des distances ) - Production de poulets ( pour voir la concurrence potentielle ) - Indice de sécurité ( pour évaluer le risque de s’implanter dans ce pays )

1.2 Installation des packages

-> L’installation des packages ne se fait qu’une seule fois, je les ai commentées pour éviter de réisntaller à chaque chargement du script

#install.packages("FactoMineR", dependendcies=T)
#install.packages("factoextra",dependencies=T)
#install.packages(c("Factoshiny"))
#install.packages("ggpubr")

1.3 Installation des libraries

-> L’activation des librairies se fait à chaque démarrage du script
-> corrplot doit être appelé après ggplot2 pour éviter un warning

library(FactoMineR)

library(factoextra)

library(ggpubr)
library(dplyr)

library(ggplot2)
library(corrplot) 

options(ggrepel.max.overlaps = Inf)

1.4 Configuration des variables pour chaque fichier

-> définition d’une variable répertoire pour faciliter le changement d’OS

repertoire <- "/home/user/Documents/formations/oc/Formation Analyst/P5/sources/"
fic_pop <- paste (repertoire,"population.csv",sep = "")
fic_pop2017 <- paste (repertoire,"pop_2017.csv",sep = "")
fic_cal <- paste(repertoire,"kcal_personne.csv",sep="")
fic_prot_volaille <- paste(repertoire,"prot_volaille.csv",sep="")
fic_poulet_prod <- paste(repertoire,"poulets_production_2019.csv",sep="")
fic_prot_ani <- paste(repertoire,"prot_ani.csv",sep="")
fic_prot_veg <- paste(repertoire,"prot_veg.csv",sep="")
fic_pib <- paste (repertoire,"pib_gdp.csv",sep = "")
fic_dis <- paste (repertoire,"distances.csv",sep = "")
fic_iso <- paste (repertoire,"sql-pays.csv",sep = "")
fic_prod <- paste (repertoire,"production_poulet.csv",sep = "")
fic_secu <- paste (repertoire,"security.csv",sep = "")
fic_full <- paste (repertoire,"full.csv", sep = "")
fic_export_pays_cluster <- paste (repertoire,"../P5_pays_cluster.csv", sep = "")
fic_export_centroid <- paste (repertoire,"../P5_centroid.csv", sep = "")

1.5 Création des Data Frame par lecture des fichiers

-> les critères d’importation, séparateurs, virgules décimales, entête etc, ont étés définis après observation des fichiers sources à l’aide d’une éditeur de texte ou d’un tableur

df_pop <- read.csv2(fic_pop,header = TRUE, sep = ",", dec = ".")
df_pop2017 <- read.csv2(fic_pop2017,header = TRUE, sep = ",", dec = ".")
df_cal <- read.csv2(fic_cal,header = TRUE, sep = ",", dec = ".")
df_prot_volaille <- read.csv2(fic_prot_volaille,header = TRUE, sep = ",", dec = ".")
df_prot_ani <- read.csv2(fic_prot_ani,header = TRUE, sep = ",", dec = ".")
df_prot_veg <- read.csv2(fic_prot_veg,header = TRUE, sep = ",", dec = ".")
df_pib <- read.csv2(fic_pib,header = TRUE, sep = ",", dec = ".")
df_dis <- read.csv2(fic_dis,header = TRUE, sep = ",", dec = ".")
df_iso <- read.csv2(fic_iso,header = FALSE, sep = ",", dec = ".")
df_prod <- read.csv2(fic_prod,header = TRUE, sep = ",", dec = ".")
df_secu <- read.csv2(fic_secu,header = TRUE, sep = ",", dec = ".")

-> Vérification d’un DF créé par importation

head(df_prot_ani)

1.6 Nettoyage des DataFrames en supprimant et renommant des colonnes

-> Je supprime les colonnes superflues et normalise les noms des variables ou les rend plus explicites

df_pop <- subset(df_pop, select = c(Code.zone, Area, Value))
names(df_pop)[names(df_pop)=="Area"]<-"Zone"
names(df_pop)[names(df_pop)=="Value"]<-"Hab1000"
df_pop2017 <- subset(df_pop2017, select = c(Code.zone, Valeur))
names(df_pop2017)[names(df_pop2017)=="Valeur"]<-"Hab2017x1000"
df_cal <- subset(df_cal, select = c(Code.zone, Valeur))
names(df_cal)[names(df_cal)=="Valeur"]<-"kcal_pers_jour"
df_prot_veg <- subset(df_prot_veg, select = c(Code.zone, Valeur))
names(df_prot_veg)[names(df_prot_veg)=="Valeur"]<-"prot_veg"
df_prot_ani <- subset(df_prot_ani, select = c(Code.zone, Valeur))
names(df_prot_ani)[names(df_prot_ani)=="Valeur"]<-"prot_ani"
df_prot_volaille <- subset(df_prot_volaille, select = c(Code.zone, Valeur))
names(df_prot_volaille)[names(df_prot_volaille)=="Valeur"]<-"prot_volaille"
df_pib <- subset(df_pib, select = c(Code.zone, Valeur))
names(df_pib)[names(df_pib)=="Valeur"]<-"PIB_usd_hab"
df_dis <- subset(df_dis, select = c(X, FRA))
names(df_dis)[names(df_dis)=="X"]<-"ISO"
names(df_dis)[names(df_dis)=="FRA"]<-"distance"
df_iso <- subset(df_iso, select = c(V6, V4))
names(df_iso)[names(df_iso)=="V6"]<-"Zone"
names(df_iso)[names(df_iso)=="V4"]<-"ISO"
df_prod <- subset(df_prod, select = c(Code.zone, Valeur))
names(df_prod)[names(df_prod)=="Valeur"]<-"Poulets1000"
df_secu <- subset(df_secu, select = c(Area.Code, Value))
names(df_secu)[names(df_secu)=="Value"]<-"securite"
names(df_secu)[names(df_secu)=="Area.Code"]<-"Code.zone"

-> Vérification d’une modification de DF

head(df_prot_ani)

-> Suppression des provinces de la Chine pour éviter les doublons car on a déjà vu que la chine était représentée deux fois, une fois le pays et une fois les provinces mais aussi la France car ce n’est pas un pays d’exportation

subset(df_pop,  Zone == "France")

df_pop <- df_pop[!df_pop$Code.zone %in% c("96","128","41","214","68"),]

1.7 Compilation des df pour former la DF globale nommée df_full

-> J’aurais été tenté de l’appeler data, mais n’ayant pas encore novice, je préfère lui donner un nom plus personnel en commençant par df ( data frame ) puis full pour plein et après j’y rajouterais les différentes formes ou compléments
-> Je suis parti des pays de la population et j’ai complété avec les autres données

df_full <- merge(df_pop,df_pop2017,by="Code.zone")
df_full <- merge(df_full,df_iso,by="Zone")
df_full <- merge(df_full,df_cal,by="Code.zone")
df_full <- merge(df_full,df_prot_veg,by="Code.zone")
df_full <- merge(df_full,df_prot_ani,by="Code.zone")
df_full <- merge(df_full,df_prot_volaille,by="Code.zone")
df_full <- merge(df_full,df_pib,by="Code.zone")
df_full <- merge(df_full,df_prod,by="Code.zone")
df_full <- merge(df_full,df_dis,by="ISO")
df_full <- merge(df_full,df_secu,by="Code.zone")

-> Vérification de la compilation des DF

head(df_full)

1.8 Colonnes calculées ( évolution population et rapport protéines )

df_full$pop_evo <-100-(100*df_full$Hab2017x1000/df_full$Hab1000)
df_full$prot_rapport <-df_full$prot_ani/(df_full$prot_veg+df_full$prot_ani)

-> suppression des colonnes qui ont servi pour les calculs

df_full <- subset(df_full, select = -c(Code.zone, ISO, Hab2017x1000,prot_veg,prot_ani, prot_volaille))

Nous vérifions les variables qui auraient des manques.

print(paste0('Il y a ', dim(subset(df_full, is.na(Hab1000)))[1]," données manquantes pour le nombre d'habitants."))

[1] "Il y a 0 données manquantes pour le nombre d'habitants."

print(paste0('Il y a ', dim(subset(df_full, is.na(kcal_pers_jour)))[1]," données manquantes pour le kcal_pers_jour."))

[1] "Il y a 0 données manquantes pour le kcal_pers_jour."

print(paste0('Il y a ', dim(subset(df_full, is.na(PIB_usd_hab)))[1]," données manquantes pour le PIB_usd_hab."))

[1] "Il y a 0 données manquantes pour le PIB_usd_hab."

print(paste0('Il y a ', dim(subset(df_full, is.na(Poulets1000)))[1]," données manquantes pour le Poulets1000."))

[1] "Il y a 0 données manquantes pour le Poulets1000."

print(paste0('Il y a ', dim(subset(df_full, is.na(distance)))[1]," données manquantes pour la distance."))

[1] "Il y a 0 données manquantes pour la distance."

print(paste0('Il y a ', dim(subset(df_full, is.na(securite)))[1]," données manquantes pour la securite."))

[1] "Il y a 0 données manquantes pour la securite."

print(paste0('Il y a ', dim(subset(df_full, is.na(pop_evo)))[1]," données manquantes pour le pop_evo."))

[1] "Il y a 0 données manquantes pour le pop_evo."

print(paste0('Il y a ', dim(subset(df_full, is.na(prot_rapport)))[1]," données manquantes pour le prot_rapport."))

[1] "Il y a 0 données manquantes pour le prot_rapport."

# on aurait pu ne garder que les pays ayant toutes les variables définies avec la commande suivante
df_full <- df_full[complete.cases(df_full),]

Vérification

head(df_full)

Ecriture d’un fichier CSV avec cette DF complête pour un contrôle extérieur si nécessaire

write.csv2(df_full, file = fic_full)

DF trié par nom de pays et nom du pays comme index

df_full_order <- df_full[order(df_full[,"Zone"],decreasing=F),]
df_full_individus<- data.frame(df_full_order[,-1], row.names=df_full_order[,1])
# vérification
head(df_full_individus) 

—————————

2 - CAH Classification Assendante Hiérarchique

—————————

Objetcif : Identifier des groupes similaires dans un jeu de données

data frame centré réduit avec scale - réduit le poids des variables à forte variance

df_full_centre_reduit <- scale(df_full_individus) 

HCPC par rapport à l’ACP sur 5 clusters

-> je fais dabord un premier calcul de l’ACP pour utiliser le résultat avec la fonction HCPC
-> j’exploite le res.pca issu de la fonction PCA sur la data frame centré réduit

ACP sur la table centré réduit

res.pca <- PCA(df_full_centre_reduit, scale.unit = FALSE, graph = FALSE)

HCPC

res.hcpc <- HCPC(res.pca,nb.clust = 5, graph = FALSE)

dendograme en 5 clusters à partir de l’HCPC

img_dendogramme <- fviz_dend(res.hcpc,
          labelsize = 2,
          cex = 0.7,                     # Taille du text
          palette = "jco",               # Palette de couleur ?ggpubr::ggpar
          rect = TRUE, rect_fill = TRUE, # Rectangle autour des groupes
          rect_border = "jco",           # Couleur du rectangle
          labels_track_height = 0.8      # Augment l'espace pour le texte
          ) + ggsave("P5_img_dendogramme.png", width = 11, height = 8)
plot(img_dendogramme)

Observation des clusters

fviz_cluster(res.hcpc,
             labelsize = 3,
             repel = TRUE,            # Evite le chevauchement des textes
             show.clust.cent = TRUE, # Montre le centre des clusters
             palette = "jco",         # Palette de couleurs, voir ?ggpubr::ggpar
             ggtheme = theme_minimal(),
             main = "Factor map"
)

rajout de la variable cluster au df

df_full_centre_reduit_cluster <-res.hcpc$data.clust
pays_cluster <- subset(df_full_centre_reduit_cluster, select=c(clust))
write.csv2(pays_cluster, file = fic_export_pays_cluster)
head(pays_cluster)

Centroides - grouby sur les cluster et moyennes de chaque varibles

df_full_cluster <- cbind(df_full_individus,pays_cluster$clust)
colnames(df_full_cluster)[colnames(df_full_cluster)=="pays_cluster$clust"] <- "clust"
df_full_centroides <- df_full_cluster %>%
  group_by(clust) %>%
  summarize(moy_hab = mean(Hab1000, na.rm = TRUE), moy_dis = mean(distance,na.rm=TRUE), moy_PIB=mean(PIB_usd_hab,na.rm=TRUE),moy_poulets=mean(Poulets1000,na.rm=TRUE),
            moy_secu=mean(securite,na.rm=TRUE),moy_kcal=mean(kcal_pers_jour,na.rm=TRUE), moy_pop_evo=mean(pop_evo,na.rm=TRUE),moy_prot=mean(prot_rapport,na.rm=TRUE) )
write.csv2(df_full_centroides, file = fic_export_centroid)
head(df_full_centroides)

df_centroid <- df_full_centroides[,2:9]
df_centroid_scaled <- scale(df_centroid)
Matrice <- data.matrix(df_centroid_scaled)
print(Matrice)

        moy_hab     moy_dis    moy_PIB moy_poulets    moy_secu   moy_kcal moy_pop_evo
[1,] -0.4009869 -0.09033865 -0.6414561  -0.4676152 -0.82938665 -1.2559588  1.57151929
[2,]  1.7880118  0.86341248 -0.5549279   1.7883461 -1.03896752 -0.3772466  0.06496343
[3,] -0.4537773  1.12973859 -0.3495361  -0.4324813  0.50847179 -0.3622087 -0.13029278
[4,] -0.4462014 -1.23398449 -0.2178675  -0.4156580 -0.04444845  0.7164432 -1.18600755
[5,] -0.4870461 -0.66882794  1.7637878  -0.4725916  1.40433082  1.2789709 -0.32018239
       moy_prot
[1,] -1.4253574
[2,] -0.4060094
[3,]  0.2949569
[4,]  0.2504042
[5,]  1.2860057
attr(,"scaled:center")
      moy_hab       moy_dis       moy_PIB   moy_poulets      moy_secu      moy_kcal 
 1.178189e+05  5.947875e+03  1.958858e+04  4.148921e+05 -1.261467e-01  2.943993e+03 
  moy_pop_evo      moy_prot 
 1.194850e+00  4.364698e-01 
attr(,"scaled:scale")
     moy_hab      moy_dis      moy_PIB  moy_poulets     moy_secu     moy_kcal  moy_pop_evo 
2.196203e+05 2.991056e+03 2.798917e+04 7.913869e+05 8.699533e-01 3.790431e+02 7.637357e-01 
    moy_prot 
1.405852e-01 

Description de la répartition des clusters

res.hcpc$desc.var

Link between the cluster variable and the quantitative variables
================================================================
                    Eta2      P-value
PIB_usd_hab    0.7467896 3.244785e-37
Poulets1000    0.6596635 4.762308e-29
prot_rapport   0.6233566 2.959460e-26
kcal_pers_jour 0.5553040 1.093472e-21
pop_evo        0.5416622 7.381464e-21
distance       0.5352245 1.781485e-20
Hab1000        0.4606421 2.109411e-16
securite       0.4350143 3.894074e-15

Description of each cluster by quantitative variables
=====================================================
$`1`
                  v.test Mean in category  Overall mean sd in category Overall sd
pop_evo         7.801112        0.9786282 -2.714581e-17     0.54097450  0.9962335
PIB_usd_hab    -4.315366       -0.5413509 -3.201413e-17     0.08106279  0.9962335
securite       -6.068048       -0.7612200  1.330389e-18     0.84216224  0.9962335
kcal_pers_jour -6.729919       -0.8442500  1.357519e-16     0.75417992  0.9962335
prot_rapport   -8.588433       -1.0773955  7.492571e-17     0.48658880  0.9962335
                    p.value
pop_evo        6.136397e-15
PIB_usd_hab    1.593388e-05
securite       1.294742e-09
kcal_pers_jour 1.697575e-11
prot_rapport   8.816324e-18

$`2`
              v.test Mean in category  Overall mean sd in category Overall sd      p.value
Poulets1000 9.317068         4.587946 -1.383214e-17       3.002534  0.9962335 1.196001e-20
Hab1000     7.778593         3.830364  1.943412e-18       3.933976  0.9962335 7.333563e-15

$`3`
                v.test Mean in category  Overall mean sd in category Overall sd
distance      7.196330        0.8039690  5.154607e-17      0.8196829  0.9962335
securite      4.480013        0.5005039  1.330389e-18      0.6781611  0.9962335
prot_rapport  4.169296        0.4657908  7.492571e-17      0.6649568  0.9962335
pop_evo      -2.098623       -0.2344567 -2.714581e-17      0.6774155  0.9962335
                  p.value
distance     6.185483e-13
securite     7.463847e-06
prot_rapport 3.055425e-05
pop_evo      3.585011e-02

$`4`
                  v.test Mean in category  Overall mean sd in category Overall sd
kcal_pers_jour  5.459027        0.9379086  1.357519e-16      0.6784478  0.9962335
prot_rapport    2.478485        0.4258254  7.492571e-17      0.6974020  0.9962335
pop_evo        -5.744707       -0.9869909 -2.714581e-17      0.9007136  0.9962335
distance       -6.337045       -1.0887597  5.154607e-17      0.2718960  0.9962335
                    p.value
kcal_pers_jour 4.787501e-08
prot_rapport   1.319416e-02
pop_evo        9.208006e-09
distance       2.342134e-10

$`5`
                  v.test Mean in category  Overall mean sd in category Overall sd
PIB_usd_hab     9.521384         3.064530 -3.201413e-17      1.2229449  0.9962335
kcal_pers_jour  4.493226         1.446179  1.357519e-16      0.2544259  0.9962335
prot_rapport    4.209311         1.354799  7.492571e-17      0.4681281  0.9962335
securite        4.180053         1.345382  1.330389e-18      0.2388849  0.9962335
distance       -1.976700        -0.636216  5.154607e-17      1.3324673  0.9962335
                    p.value
PIB_usd_hab    1.708871e-21
kcal_pers_jour 7.015228e-06
prot_rapport   2.561507e-05
securite       2.914417e-05
distance       4.807553e-02

Heatmap des variables dans chaque cluster pour voir leur importance

corrplot(Matrice, is.corr = FALSE, method = "circle")

Interprétation : -> plus le point est bleu, plus la variable caractérise le cluster
-> plus le point est grand, plus la variable est forte
-> plus le point est rouge moins la variable caractérise le cluster
-> plus le point est petit plus la variable est faible

CLUSTER 1
- évolution de la population est forte
- kcal et protéines sont très faibles
CLUSTER 2 - nombre d’habitants et de poulets élevés
- distance un peu favorable
- la sécurité est en revanche défavorable
CLUSTER 3
- seul la distance est l’atout de cet ensemble
CLUSTER 4
- distance et évolution de la population sont deux variables particulièrement fortes et négatives pour cet ensemble
CLUSTER 5
- des pays riches et sûrs
- une alimentation riche en protéines et en kcal

-> Le groupe 5 présente les plus intéressants paramètres pour notre projet car on y trouve les meilleurs résultats en matière de richesse et de sécurité, ce qui permet d’envisager un commerce plus stable et moins risqué. De plus l’alimentation est plutôt riche.

—————————

3 ACP - Analyse de Composantes Principales

—————————

res.desc <- dimdesc(res.pca, axes = c(1,2), proba = 0.05)
res.desc$Dim.1

$quanti
               correlation      p.value
prot_rapport     0.9031001 6.236248e-50
kcal_pers_jour   0.7653118 7.768570e-27
securite         0.7587525 3.734427e-26
PIB_usd_hab      0.7404124 2.337017e-24
Hab1000         -0.2167121 1.222769e-02
pop_evo         -0.7111666 8.718108e-22

attr(,"class")
[1] "condes" "list"  

Eboulis des valeurs propres : pourcentage de variance expliquée par chaque dimension

eig.val <- get_eigenvalue(res.pca) 
eig.val

      eigenvalue variance.percent cumulative.variance.percent
Dim.1  3.0903096        38.921512                    38.92151
Dim.2  1.4739074        18.563418                    57.48493
Dim.3  1.2049839        15.176407                    72.66134
Dim.4  0.7415876         9.340071                    82.00141
Dim.5  0.5860269         7.380831                    89.38224
Dim.6  0.3829458         4.823087                    94.20533
Dim.7  0.2666754         3.358696                    97.56402
Dim.8  0.1934130         2.435978                   100.00000

fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50))

-> On observe un coude au niveau de la seconde dimension. Après quoi la diminution est plus régulière.
Mais il sera peut-être nécessaire d’observer les 3 premières composantes principales pour totaliser 70% des informations.

Comparaison de la qualité de la représentation des variables en fonction des dimensions

var <- get_pca_var(res.pca)
corrplot(var$cos2, is.corr=FALSE,addCoef.col = "black", title="Qualité de représentation des variables", mar=c(0,0,1,0))

-> dans la première composante, ce sont kcal, pib, sécurité et rapport en protéines qui sont le mieux représentées
-> dans la seconde, c’est le nombre de poulets et d’habitants.
-> dans la troisième, c’est la distance qui est le mieux représentée

Graphique des corrélations pour les dimensions 1 et 2

fviz_pca_var (res.pca, col.var = "cos2",
              gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
              repel = TRUE # Évite le chevauchement de texte
)

fviz_pca_var(res.pca, col.var = "cos2",gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),repel = TRUE)+ labs(title = "Cercle des corrélations") + theme(plot.title = element_text(color="#3876C2", size=20, face="bold",hjust = 0.5),axis.title.x = element_text(color="black", size=14, face="bold",hjust = 0.5),axis.title.y = element_text(color="black", size=14, face="bold",vjust = 0.5))

Le graphique de corrélation des variables montre les relations entre toutes les variables et les axes factoriels. Il peut être interprété comme suit:

Les variables positivement corrélées sont regroupées :
    disp_ali_prot_g_pers_jour, disp_ali_kcal_pers_jour
    prop_ani_disp_ali, pib_hab
Les variables négativement corrélées sont positionnées sur les côtés opposés de l'origine du graphique (quadrants opposés) :
    indice_risque_pays et evolution_pop_2015_2016 par rapport aux autres
La distance entre les variables et l'origine mesure la qualité de représentation des variables. Les variables qui sont loin de l'origine sont bien représentées par l'ACP (couleur rouge / orange). Toutes les variables (un peu moins pour evolution_pop_2015_2016 : couleur bleue) sont bien expliquées par les deux premières composantes principales (Dim.1 & Dim.2) car elles sont positionnées près du le cercle de corrélation.

Soit le repère orthogonal porté par les axes (distance_fr_fm et un axe normal à ce dernier).

La dimension 1 symbolise la qualité de vie dans le pays : les pays les plus à droite sont les pays où la qualité de vie est la meilleur (fort PIB par habitant, forte disponibilité alimentaire en kcal / protéines, forte consommation de produits d'origine animale, bonne stabilité politique).
La dimension 2 symbolise la distance du pays à la France. : plus un pays est en bas, plus il est proche de la France.
Les pays qui nous intéressent sont donc les pays les plus à droite (prioritairement car la distance reste secondaire au vue de la qualité de vie dans le pays) et les plus en bas.

De manière générale, le cercle des corrélations n’est pas fait pour interpréter les corrélations entre les variables initiales. Il est fait pour interpréter les axes d’inertie F1, F2, etc. Pour vérifier la corrélation entre les variables, les coefficients de corrélations entre les variables sont calculés et affichés dans la matrice des corrélations.

qualité de représentation ( cos2 ) barplot

Variables

fviz_cos2(res.pca, choice = "var")

—————————

4 Clusters

—————————

4.1 Analyse des différents groupes

cluster 1 - évolution de la population sinificativement supérieur à la moyenne

cluster 2 - nombre de poulets sinificativement supérieur à la moyenne

cluster 3 - distances sinificativement supérieur à la moyenne

cluster 4 - évolution des protéines sinificativement supérieur à la moyenne

cluster 5 - PIB sinificativement supérieur à la moyenne

Représentation visuelle des Clusters

df_full_centre_reduit_cluster_trans <- transform(df_full_centre_reduit_cluster, clust = as.factor(clust)) 
fviz_pca_biplot(res.pca, pointsize = 1, labelsize = 2, axes = c(1,2), select.ind = list(cos2 = 0.50),addEllipses = TRUE, ellipse.level=0.60,col.ind = df_full_centre_reduit_cluster$clust, color_palette=c("#00AFBB", "#E7B800", "#FC4E07","#FC4E07","#FC4E07")) +  theme(plot.title = element_text(color="#3876C2", size=20, face="bold",hjust = 0.5),axis.title.x = element_text(color="black", size=14, face="bold",hjust = 0.5),axis.title.y = element_text(color="black", size=14, face="bold",vjust = 0.5)) + ggsave("P5_img_biplot_pca.png", width = 11, height = 8)

-> Cette représentation confirme le choix du cluster 5 car c’est le cluster qui s’approche le plus des consommations souhaitées, mais aussi pour la richesse et la sécurité

4.2 Choix du cluster 5

df_full_individus

df_full_individus$Area <- rownames(df_full_individus)
df_clust5 <- subset(df_full_centre_reduit_cluster, clust == 5)
#df_clust5_full <- df_clust5[,0]
#df_clust5_full$Area <- rownames(df_clust5_full)
#df_clust5_full <- merge(df_clust5_full,df_full_individus, by = 'Area')
#df_clust5 <- subset(df_full_centre_reduit_cluster, clust == 5)
print(df_clust5[,0])

—————————

5 Cluster numéro 5 - Application de nouvelles analyses

—————————

5.1 Suppression de variables

clust car on garde tous les pays
PIB_usb_hab car tous ces pays sont riches
securite car tous ces pays sont sûrs
Poulets1000 car ils n’ont pas une grande quantité de volailles
pop_evo car l’évolution est semblable dans tous ces pays ( plutôt faible )

df_clust5 <- subset(df_clust5, select = -c(clust,PIB_usd_hab,Poulets1000, pop_evo, prot_rapport,securite))

5.2 Ajout d’une nouvelle variable pour essayer de distinguer les individus de ce groupe

# ajout de la production de poulets
fic_poulet_prod <- paste(repertoire,"poulets_production_2019.csv",sep="")
df_poulet_prod <- read.csv2(fic_poulet_prod,header = TRUE, sep = ",", dec = ".")
df_poulet_prod <- subset(df_poulet_prod,select = c(Area,Value))
names(df_pop)[names(df_pop)=="Zone"]<-"Area"
df_poulet_prod <- merge(df_poulet_prod,df_pop,by="Area")
df_poulet_prod$prod_poulet_hab <-df_poulet_prod$Value/df_poulet_prod$Hab1000
df_poulet_prod <- subset(df_poulet_prod,select = c(Area,prod_poulet_hab))

df_clust5$Area <- rownames(df_clust5)
#df_clust5_prod <- merge(df_clust5,df_poulet_prod,by="Area")
#names(df_clust5_prod)[names(df_clust5_prod)=="Value"]<-"Prod_PouletX1000"

#df_clust5_prod<- data.frame(df_clust5_prod[,-1], row.names=df_clust5_prod[,1])
#df_full_individus
df_full_individus$Area <- rownames(df_full_individus)

df_clust5_full <- df_clust5[,0]
df_clust5_full$Area <- rownames(df_clust5_full)
df_clust5_full <- merge(df_clust5_full,df_full_individus, by = 'Area')
df_clust5_full <- merge(df_clust5_full,df_poulet_prod, by = 'Area')
df_clust5_full <- subset(df_clust5_full, select = -c(PIB_usd_hab,Poulets1000, pop_evo, prot_rapport,securite))
df_clust5_full<- data.frame(df_clust5_full[,-1], row.names=df_clust5_full[,1])

df_clust5 <- subset(df_clust5, select = -c(Area))

#df_clust5 <- subset(df_full_centre_reduit_cluster, clust == 5)
#print(df_clust5[,0])

Affichage des tables pour contrôle

df_pop

df_clust5

df_poulet_prod

#df_clust5_prod
df_clust5_full

6.3 ACP sur la table centré réduit du cluster 5

res5.pca <- PCA(df_clust5, scale.unit = FALSE, graph = FALSE)

6.4 HCPC par rapport à l’ACP sans préciser le nombre de clusters

res5.hcpc <- HCPC(res5.pca, graph = FALSE) 

6.5 dendograme à partir de l’HCPC

fviz_dend(res5.hcpc, 
          cex = 0.7,                     # Taille du text
          palette = "jco",               # Palette de couleur ?ggpubr::ggpar
          rect = TRUE, rect_fill = TRUE, # Rectangle autour des groupes
          rect_border = "jco",           # Couleur du rectangle
          labels_track_height = 0.8      # Augment l'espace pour le texte
)

-> Le dendagramme a constitué 3 clusters ce qui semble un nombre correct de groupes compte tenu du nombre de pays dans chacun d’eux.

Justifier le choix des variables pib et sécurité corrélés mais tous ces pays sont riches et sécurisés

Observation des clusters

fviz_cluster(res5.hcpc,
             labelsize = 8 ,
             show.clust.cent = TRUE, # Montre le centre des clusters
             palette = "jco",         # Palette de couleurs, voir ?ggpubr::ggpar
             ggtheme = theme_minimal(),
             main = "Factor map"
)

fviz_pca_var (res5.pca, col.var = "cos2",
              gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
              repel = TRUE # Évite le chevauchement de texte
)

df_clust5_scaled <- as.data.frame(df_clust5)
res5.pca <- PCA(df_clust5_scaled, graph = FALSE)

éboulis de valeurs propres

eig5.val <- get_eigenvalue(res5.pca)

perte d’inertie en fonction des dimensions

fviz_eig(res5.pca, addlabels = TRUE, ylim = c(0, 50))

var5 <- get_pca_var(res5.pca)
corrplot(var5$cos2, is.corr=FALSE, title="Qualité des variables", mar=c(0,0,1,0))

# Colorer en fonction du cos2: qualité de représentation
fviz_pca_var(res5.pca, col.var = "cos2",gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),repel = TRUE) + labs(title = "Figure 10 - Cercle des corrélations") + theme(plot.title = element_text(color="#3876C2", size=20, face="bold",hjust = 0.5),axis.title.x = element_text(color="black", size=14, face="bold",hjust = 0.5),axis.title.y = element_text(color="black", size=14, face="bold",vjust = 0.5)) 

NA
NA

fviz_pca_biplot(res5.pca, col.var = "#2E9FDF", pointsize = 1, labelsize = 3, axes = c(1,2)) +  theme(plot.title = element_text(color="#3876C2", size=20, face="bold",hjust = 0.5),axis.title.x = element_text(color="black", size=14, face="bold",hjust = 0.5),axis.title.y = element_text(color="black", size=14, face="bold",vjust = 0.5))

df_clust_scaled_cluster <-res5.hcpc$data.clust
df_clust_scaled_cluster

pays_cluster <- subset(df_clust_scaled_cluster, select=c(clust))
# vérification
# summary(df_full_cluster)
pays_cluster

df_clust_scaled_cluster_trans <- transform(df_clust_scaled_cluster, clust = as.factor(clust)) 
fviz_pca_biplot(res5.pca, pointsize = 1, labelsize = 2, axes = c(1,2), select.ind = list(cos2 = 0.30),addEllipses = TRUE, ellipse.level=0.60,col.ind = df_clust_scaled_cluster$clust, color_palette=c("#00AFBB", "#E7B800", "#FC4E07","#FC4E07","#FC4E07")) +  theme(plot.title = element_text(color="#3876C2", size=20, face="bold",hjust = 0.5),axis.title.x = element_text(color="black", size=14, face="bold",hjust = 0.5),axis.title.y = element_text(color="black", size=14, face="bold",vjust = 0.5))

df_clust_scaled_cluster_trans <- transform(df_clust_scaled_cluster, clust = as.factor(clust)) 
fviz_pca_biplot(res5.pca, pointsize = 1, labelsize = 2, axes = c(1,3), select.ind = list(cos2 = 0.30),addEllipses = TRUE, ellipse.level=0.60,col.ind = df_clust_scaled_cluster$clust, color_palette=c("#00AFBB", "#E7B800", "#FC4E07","#FC4E07","#FC4E07")) +  theme(plot.title = element_text(color="#3876C2", size=20, face="bold",hjust = 0.5),axis.title.x = element_text(color="black", size=14, face="bold",hjust = 0.5),axis.title.y = element_text(color="black", size=14, face="bold",vjust = 0.5))

NA
NA

En faisant un ACP sur l’axe 1.2 et 1.3, on observe un groupe se distinguer des autres, le cluster 3 Islande, Luxembourg, Suisse, Finlande et Norvège.

si l’on souhaite réduire encore l’échantillon, je propose de prendre en compte la distance, ce qui réduira la sélection au Luxembourg et la Suisse, lesquels sont frontaliers, riches et qui plus est francophones.

trouver les informations concernant la production, l’importation et éventuellement la consommation

test avec df_clust5_full # ################################################### ### ACP sur la table centré réduit du cluster 5

df_clust5_full

df_clust5_full_centre_reduit

                        Hab1000 kcal_pers_jour    distance prod_poulet_hab
Australia             0.5396770     -0.5428360  2.87116921      0.21115057
Canada                1.1777444      1.1636845  0.51793967      0.24849363
Denmark              -0.4636119     -0.4453206 -0.55487032      0.03870601
Finland              -0.4756407     -1.0109102 -0.36466099     -0.59947471
France                2.6405942      0.5395855 -0.77616541     -0.03294302
Iceland              -0.7473900      2.0218205 -0.29529018     -0.90547345
Israel               -0.3258267      0.7931257 -0.05861119      0.60724951
Kuwait               -0.5482312      0.2372876  0.17909646      2.75647840
Luxembourg           -0.7333708      0.1787783 -0.71422707     -1.03864342
Norway               -0.4853149     -0.7378670 -0.48708771     -0.23466830
Switzerland          -0.3182759     -0.9036432 -0.67070839     -0.69702870
United Arab Emirates -0.2603535     -1.2937051  0.35341593     -0.35384653
attr(,"scaled:center")
        Hab1000  kcal_pers_jour        distance prod_poulet_hab 
   14599.353750     3446.666667     3608.653838        3.772712 
attr(,"scaled:scale")
        Hab1000  kcal_pers_jour        distance prod_poulet_hab 
   19083.264400      102.547845     4649.336044        3.422019 

df_clust5_full_centre_reduit <- scale(df_clust5_full)
res5.pca <- PCA(df_clust5_full_centre_reduit, scale.unit = FALSE, graph = FALSE)

HCPC par rapport à l’ACP et je laisse libre le nombre de clusters

res5.hcpc <- HCPC(res5.pca, graph = FALSE) # nb.clust = 4,

dendograme à partir de l’HCPC

fviz_dend(res5.hcpc, 
          cex = 0.7,                     # Taille du text
          palette = "jco",               # Palette de couleur ?ggpubr::ggpar
          rect = TRUE, rect_fill = TRUE, # Rectangle autour des groupes
          rect_border = "jco",           # Couleur du rectangle
          labels_track_height = 0.8      # Augment l'espace pour le texte
)

Observation des clusters

fviz_cluster(res5.hcpc,
             labelsize = 8 ,
             show.clust.cent = TRUE, # Montre le centre des clusters
             palette = "jco",         # Palette de couleurs, voir ?ggpubr::ggpar
             ggtheme = theme_minimal(),
             main = "Factor map"
)

fviz_pca_var (res5.pca, col.var = "cos2",
              gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
              repel = TRUE # Évite le chevauchement de texte
)

res5.pca <- PCA(df_clust5_full_centre_reduit, graph = FALSE)

éboulis de valeurs propres

eig5.val <- get_eigenvalue(res5.pca)

perte d’inertie en fonction des dimensions

fviz_eig(res5.pca, addlabels = TRUE, ylim = c(0, 50))

var5 <- get_pca_var(res5.pca)
corrplot(var5$cos2, is.corr=FALSE, title="Qualité des variables", mar=c(0,0,1,0))

# Colorer en fonction du cos2: qualité de représentation
fviz_pca_var(res5.pca, col.var = "cos2",gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),repel = TRUE) + labs(title = "Figure 10 - Cercle des corrélations") + theme(plot.title = element_text(color="#3876C2", size=20, face="bold",hjust = 0.5),axis.title.x = element_text(color="black", size=14, face="bold",hjust = 0.5),axis.title.y = element_text(color="black", size=14, face="bold",vjust = 0.5)) 

NA
NA

fviz_pca_biplot(res5.pca, col.var = "#2E9FDF", pointsize = 1, labelsize = 3, axes = c(1,2)) +  theme(plot.title = element_text(color="#3876C2", size=20, face="bold",hjust = 0.5),axis.title.x = element_text(color="black", size=14, face="bold",hjust = 0.5),axis.title.y = element_text(color="black", size=14, face="bold",vjust = 0.5))

fviz_pca_biplot(res5.pca, pointsize = 1, labelsize = 3, axes = c(1,2), select.ind = list(cos2 = 0.30),addEllipses = TRUE, ellipse.level=0.60,col.ind = df_clust_scaled_cluster$clust, color_palette=c("#00AFBB", "#E7B800", "#FC4E07","#FC4E07","#FC4E07")) +  theme(plot.title = element_text(color="#3876C2", size=20, face="bold",hjust = 0.5),axis.title.x = element_text(color="black", size=14, face="bold",hjust = 0.5),axis.title.y = element_text(color="black", size=14, face="bold",vjust = 0.5))

df_clust_scaled_cluster_trans <- transform(df_clust_scaled_cluster, clust = as.factor(clust)) 
fviz_pca_biplot(res5.pca, pointsize = 1, labelsize = 2, axes = c(1,3), select.ind = list(cos2 = 0.30),addEllipses = TRUE, ellipse.level=0.60,col.ind = df_clust_scaled_cluster$clust, color_palette=c("#00AFBB", "#E7B800", "#FC4E07","#FC4E07","#FC4E07")) +  theme(plot.title = element_text(color="#3876C2", size=20, face="bold",hjust = 0.5),axis.title.x = element_text(color="black", size=14, face="bold",hjust = 0.5),axis.title.y = element_text(color="black", size=14, face="bold",vjust = 0.5))

En faisant un ACP sur l’axe 1.2 et 1.3, on observe un groupe se distinguer des autres, le cluster 3 Islande, Luxembourg, Suisse, Finlande et Norvège.

si l’on souhaite réduire encore l’échantillon, je propose de prendre en compte la distance, ce qui réduira la sélection au Luxembourg et la Suisse, lesquels sont frontaliers, riches et qui plus est francophones.

trouver les informations concernant la production, l’importation et éventuellement la consommation

—————————

7 - TESTS STATISTIQUES

—————————

Affichage des histogrames de plusieurs variables pour voir si l’on peut y retrouver la loi normale

# multiplot : Pour afficher plusieurs graphiques 

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

p1 <- ggplot(df_full_individus, aes(x=kcal_pers_jour)) + geom_histogram(aes(y=..density..), colour="black", fill="white",bins=30)+
 geom_density(alpha=.2, fill="#FF6666") +ggsave("img/figure12_histo.png", width = 15, height = 8)
p2 <- ggplot(df_full_individus, aes(x=prot_rapport)) + geom_histogram(aes(y=..density..), colour="black", fill="white",bins=30)+
 geom_density(alpha=.3, fill="#FF6666") 
p3 <- ggplot(df_full_individus, aes(x=Poulets1000))+ geom_histogram(aes(y=..density..), colour="black", fill="white",bins=30)+
 geom_density(alpha=.3, fill="#FF6666") 
p4 <- ggplot(df_full_individus, aes(x=pop_evo))+ geom_histogram(aes(y=..density..), colour="black", fill="white",bins=20)+
 geom_density(alpha=.3, fill="#FF6666") 
p5 <- ggplot(df_full_individus, aes(x=PIB_usd_hab)) + geom_histogram(aes(y=..density..), colour="black", fill="white",bins=6)+
 geom_density(alpha=.3, fill="#FF6666") 
p6 <- ggplot(df_full_individus, aes(x=distance)) + geom_histogram(aes(y=..density..), colour="black", fill="white",bins=50)+
 geom_density(alpha=.3, fill="#FF6666") 
p7 <- ggplot(df_full_individus, aes(x=securite)) + geom_histogram(aes(y=..density..), colour="black", fill="white",bins=8)+
 geom_density(alpha=.3, fill="#FF6666") 
multiplot(p1, p2, p3, p4,p5,p6,p7, cols=2)

La loi normale se représente par une courbe semblable à une cloche d’église avec un point culminant au milieu

Deux histogrammes s’apparentent à une loi normale, la sécurité et l’évolution de la population

Nous allons donc réaliser des test statistiques sur ces deux variables pour voir si elles respectent la loi normale

—————————

7.1 - TEST D’ADEQUATION

—————————

Pour la sécurité

df_secu_test <- transform(df_secu, securite = as.numeric(securite))
ks.test(df_secu$securite,"pnorm",mean=mean(df_secu$securite))

aucun ex-aequo ne devrait être présent pour le test de Kolmogorov-Smirnov

    One-sample Kolmogorov-Smirnov test

data:  df_secu$securite
D = 0.064163, p-value = 0.3888
alternative hypothesis: two-sided

La p-valeur est supérieure à 5%. On ne peut donc pas rejeter l’hypothèse de normalité au niveau du test à 5%. L’indice de la sécurité suit une loi normale.

Pour l’évolution de la population

df_pop_test <- transform(df_pop, Hab1000 = as.numeric(Hab1000))
ks.test(df_pop$Hab1000,"pnorm",mean=mean(df_pop$Hab1000))

    One-sample Kolmogorov-Smirnov test

data:  df_pop$Hab1000
D = 0.82759, p-value < 2.2e-16
alternative hypothesis: two-sided

Quand à la population, avec une p-valeur de 2.2e-16 on peut dire que l’indice de la population ne suit pas une loi normale.

loi normale sur cluster 4 et 5 pour sécurité

clust4_secu <- filter(df_full_centre_reduit_cluster, clust == '4')$securite
clust5_secu <- filter(df_full_centre_reduit_cluster, clust == '5')$securite

loi normale sur la sécurité du cluster 4

ks.test(clust4_secu,"pnorm",mean=mean(clust4_secu))

aucun ex-aequo ne devrait être présent pour le test de Kolmogorov-Smirnov

    One-sample Kolmogorov-Smirnov test

data:  clust4_secu
D = 0.12197, p-value = 0.8166
alternative hypothesis: two-sided

La p-valeur est supérieure à 5%. On ne peut donc pas rejeter l’hypothèse de normalité au niveau du test à 5%. L’indice de la sécurité suit une loi normale.

loi normale sur la sécurité du cluster 5

ks.test(clust5_secu,"pnorm",mean=mean(clust5_secu))

    One-sample Kolmogorov-Smirnov test

data:  clust5_secu
D = 0.3666, p-value = 0.1366
alternative hypothesis: two-sided

Pour le cluster 5. La p-valeur est inférieure à 5%. On peut donc rejeter l’hypothèse de normalité au niveau du test à 5%. L’indice de la sécurité ne suit pas une loi normale.

—————————

7.2 - TEST DE COMPARAISON

—————————

D’après le graphe on peut voir que les cluster 4 et 5 sont proches, nous allons faire une comparaison des variances sur la variable de sécurité pour ces deux clusters

var.test(clust4_secu,clust5_secu)

    F test to compare two variances

data:  clust4_secu and clust5_secu
F = 10.747, num df = 26, denom df = 8, p-value = 0.00167
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
  2.736687 29.330586
sample estimates:
ratio of variances 
          10.74663 

La p-valeur valant 0.2833, on ne rejette donc pas l’égalité des variances au niveau de test 5%.

TEST D’ÉGALITÉ DES MOYENNES

t.test(clust5_secu,clust4_secu,var.equal=TRUE,alternative="greater")

    Two Sample t-test

data:  clust5_secu and clust4_secu
t = 4.8187, df = 34, p-value = 1.474e-05
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
 0.8868746       Inf
sample estimates:
  mean of x   mean of y 
 1.34538169 -0.02095088 

La p-valeur étant proche de 0 on rejette l’hypothèse d’égalité des moyennes. Les moyennes μ4 et μ5 ne sont donc pas égales.

En conclusion :

—————–

La variable sécurité a été soumise à des tests statistiques dans le but de savoir si elle suivait les mêmes influences sur deux clusters différents mais proches.

L’aspect gaussien des variables sécurité a été vérifié par un test d’adéquation via une loi normale par le test de Kolmogorov-Smirnov. Pour les clusters 4 et 5, l’hypothèse d’égalité des variances a été vérifiée.

contrairement à l’hypothèse d’égalité des moyennes de la variable sécurité pour les clusters 4 et 5 qui ne suivent pas la même distribution sur les deux clusters et sont donc distincts.