Votre mission Pour identifier les pays propices à une insertion dans le marché du poulet, il vous a été demandé de cibler les pays. Il vous faudra également étudier les régimes alimentaires de chaque pays, notamment en termes de protéines d’origine animale et en termes de calories.
Construisez votre échantillon contenant l’ensemble des pays disponibles, chacun caractérisé par ces variables :
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. Construisez un dendrogramme contenant l’ensemble des pays étudiés, puis coupez-le afin d’obtenir 5 groupes.
Caractérisez chacun de ces groupes selon les variables cités précédemment, et facultativement selon d’autres variables que vous jugerez pertinentes (ex : le PIB par habitant). Vous pouvez le faire en calculant la position des centroïdes de chacun des groupes, puis en les commentant et en les critiquant au vu de vos objectifs.
Donnez une courte liste de pays à cibler, en présentant leurs caractéristiques. Un découpage plus précis qu’en 5 groupes peut si besoin être effectué pour cibler un nombre raisonnable de pays.
Visualisez vos partitions dans le premier plan factoriel obtenu par ACP.
Dans votre partition, vous avez obtenu des groupes distincts. Vérifiez donc qu’ils diffèrent réellement. Pour cela, réalisez les tests statistiques suivants :
un test d’adéquation : parmi les 4 variables, ou parmi d’autres variables que vous trouverez pertinentes, trouvez une variable dont la loi est normale ; un test de comparaison de deux populations (dans le cas gaussien) : choisissez 2 clusters parmi ceux que vous aurez déterminé. Sur ces 2 clusters, testez la variable gaussienne grâce à un test de comparaison.
VIM::aggr(df9,only.miss = T)
visdat:: vis_miss(df9)
## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## Please use `gather()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
On utilise le package MissMDA pour estimer les valeurs manquantes avant l’ACP
df9completed <- missMDA::imputePCA(select(df9, -code_pays, -dispo_poulet), ncp=2,quali.sup=c(1))$completeObs
head(df9completed)
## pays evolution_pop Dispo_kcal_p_j dispo_prot_p_j
## 1 Égypte 1.088460 3229 97.49
## 2 Émirats arabes unis 1.045233 3192 110.12
## 3 Équateur 1.070996 2509 65.94
## 4 États-Unis d'Amérique 1.026431 3657 113.58
## 5 Éthiopie 1.113464 2295 66.33
## 6 Afghanistan 1.113906 2020 55.50
## ratio_prot PIB_habitant inflation dispo_kcal_volaille Export Quantity
## 1 0.2179194 25.412600 13.2918710 2386334.0 2.0000
## 2 0.5260204 7.859682 3.5448870 3302500.2 58.0000
## 3 0.4470577 1.326688 -1.6191844 549640.0 0.0000
## 4 0.6340153 4.820792 0.4559698 26632711.0 3776.0000
## 5 0.3109497 1.747037 12.6079012 -596247.5 -125.6965
## 6 0.2648328 -5.691074 -1.0459523 -253312.1 0.0000
## Import Quantity
## 1 46
## 2 609
## 3 0
## 4 128
## 5 1
## 6 25
On ne modifie pas le nom des pays malgré les caractères spéciaux. Ces derniers sont retranscrits comme il se doit une fois exportés au format csv/excel.
acp9<-PCA(df9completed,quali.sup=c(1),graph=FALSE)
plot.PCA(acp9,choix='var',title="Graphe des variables de l'ACP")
plot.PCA(acp9,title="Graphe des individus de l'ACP")
fviz_eig(acp9, addlabels = TRUE, ylim = c(0, 50)) # 83 % de variance sur 5 prèmieres dimensions
83 % de variance sur 5 premières dimensions.
var <- get_pca_var(acp9)
head(var$cos2)
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## evolution_pop 0.431410639 0.01037855 0.2719627001 0.009464019 0.0168063517
## Dispo_kcal_p_j 0.725494200 0.01558782 0.0004319632 0.017246111 0.0069425156
## dispo_prot_p_j 0.750885473 0.05486327 0.0088962594 0.016899207 0.0140821552
## ratio_prot 0.596366948 0.08127893 0.0230402716 0.045164685 0.0134959480
## PIB_habitant 0.038683137 0.32053170 0.0136810437 0.307699512 0.2958007102
## inflation 0.005192124 0.30444070 0.5058944333 0.090736427 0.0003137293
head(var$coord)
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## evolution_pop -0.65681857 0.1018752 0.52150043 -0.09728319 -0.12963931
## Dispo_kcal_p_j 0.85175947 -0.1248512 0.02078373 -0.13132445 -0.08332176
## dispo_prot_p_j 0.86653648 -0.2342291 -0.09431998 -0.12999695 -0.11866826
## ratio_prot 0.77224798 -0.2850946 -0.15179022 -0.21251985 -0.11617206
## PIB_habitant 0.19668029 -0.5661552 0.11696599 0.55470669 0.54387564
## inflation -0.07205639 0.5517614 -0.71126256 0.30122488 -0.01771240
head(var$contrib)
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## evolution_pop 12.5865693 0.685178 25.23736802 1.048235 2.07133470
## Dispo_kcal_p_j 21.1665689 1.029087 0.04008496 1.910180 0.85564516
## dispo_prot_p_j 21.9073690 3.621999 0.82554767 1.871757 1.73558529
## ratio_prot 17.3992323 5.365925 2.13807192 5.002441 1.66333693
## PIB_habitant 1.1285952 21.161070 1.26956209 34.080804 36.45659019
## inflation 0.1514822 20.098763 46.94557007 10.049968 0.03866623
corrplot::corrplot(var$cos2, is.corr=FALSE)
ggcorrplot::ggcorrplot(var$cos2,lab = T)
Dimension 4 & 5 : Corrélées avec les quantités importées.
# Cos2 total des variables sur Dim.1 et Dim.2
fviz_cos2(acp9, choice = "var", axes = 1:2)
# Colorer en fonction du cos2: qualité de représentation
fviz_pca_var(acp9, col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Évite le chevauchement de texte
)
corrplot:: corrplot(var$contrib, is.corr=FALSE)
# Contributions des variables à PC1
fviz_contrib(acp9, choice = "var", axes = 1, top = 10)
# Contributions des variables à PC2
fviz_contrib(acp9, choice = "var", axes = 2, top = 10)
# dim 1 et 2
fviz_contrib(acp9, choice = "var", axes = 1:2, top = 10)
# Contribution
fviz_pca_var(acp9, col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")
)
lecture :
clus_acp9 <- HCPC (acp9, nb.clust = 5, graph = TRUE)
library(dendextend)
##
## ---------------------
## Welcome to dendextend version 1.15.2
## Type citation('dendextend') for how to cite the package.
##
## Type browseVignettes(package = 'dendextend') for the package vignette.
## The github page is: https://github.com/talgalili/dendextend/
##
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
## You may ask questions at stackoverflow, use the r and dendextend tags:
## https://stackoverflow.com/questions/tagged/dendextend
##
## To suppress this message use: suppressPackageStartupMessages(library(dendextend))
## ---------------------
##
## Attachement du package : 'dendextend'
## L'objet suivant est masqué depuis 'package:stats':
##
## cutree
# Dendogram
fviz_dend(clus_acp9, k = 5, show_labels = FALSE, rect = TRUE)
# plan factoriel
fviz_cluster(clus_acp9,
repel = TRUE, # Evite le chevauchement des textes
show.clust.cent = TRUE, # Montre le centre des clusters
palette = "jco",
ggtheme = theme_minimal(),
main = "Factor map"
)
# Nb de cluster
k4 <- HCPC (acp9, nb.clust = 4, graph = F)
k5 <- HCPC (acp9, nb.clust = 5, graph = F)
k6 <- HCPC (acp9, nb.clust = 6, graph = F)
k7 <- HCPC (acp9, nb.clust = 7, graph = F)
p4 <- fviz_cluster(k4,
repel = TRUE,
show.clust.cent = TRUE,
palette = "jco",
ggtheme = theme_minimal(),
main = "k4")
p5 <- fviz_cluster(k5,
repel = TRUE,
show.clust.cent = TRUE,
palette = "jco",
ggtheme = theme_minimal(),
main = "k5")
p6 <- fviz_cluster(k6,
repel = TRUE,
show.clust.cent = TRUE,
palette = "jco",
ggtheme = theme_minimal(),
main = "k6")
p7 <- fviz_cluster(k7,
repel = TRUE,
show.clust.cent = TRUE,
palette = "jco",
ggtheme = theme_minimal(),
main = "k7")
library(gridExtra)
##
## Attachement du package : 'gridExtra'
## L'objet suivant est masqué depuis 'package:dplyr':
##
## combine
grid.arrange(p4,p5,p6,p7, nrow = 2)
## Warning: ggrepel: 165 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 165 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 165 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 165 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
# On a bien les pays consommateurs et importateurs en 5 clusters,
#k=4 ne nous donne pas cette information , k=6 nous donne des informations non pertinentes
# Liste des pays par cluster
liste <- list(c(3,4,2,5,2,2,4,3,2,4,2,3,4,3,3,3,3,3,2,3,2,3,2,3,3,2,2,3,2,5,3,2,2,2,2,2,3,3,4,3,5,3,3,2,3,3,3,3,2,3,2,3,3,2,3,2,3,4,3,2,2,2,3,2,2,2,2,3,2,2,3,2,2
,2,2,3,3,3,3,3,4,2,3,2,2,2,3,2,3,2,2,3,3,3,2,3,2,2,2,3,3,3,2,4,2,3,2,3,2,2,2,2,2,3,2,3,2,2,2,2,2,3,2,2,4,2,3,3,3,2,2,2,3,3,3,2,3,4,2,2,3,3,3,3,2,3
,2,3,3,2,2,3,3,2,2,3,2,2,2,2,3,3,3,3,3,3,2,1,4,2,2,2))
liste <- as.data.frame(liste)
colnames(liste) <- c("cluster")
Z <- cbind(df9completed,liste) %>% select(pays,cluster)
head(Z)
## pays cluster
## 1 Égypte 3
## 2 Émirats arabes unis 4
## 3 Équateur 2
## 4 États-Unis d'Amérique 5
## 5 Éthiopie 2
## 6 Afghanistan 2
# On exporte la liste des pays par cluster au format csv
write.table(Z,"cluster_pays.csv",sep=";",row.names=FALSE)
# On visualise le cluster 4
cbind(df9completed,liste) %>% filter(cluster == 4) %>% view()
res.cluster4<-PCA(cbind(df9completed, liste) %>% filter(cluster == 4),quali.sup=c(1),graph=FALSE)
plot.PCA(res.cluster4,invisible=c('ind','ind.sup'),
title="Pays importateurs",cex=1.1,cex.main=1.1,cex.axis=1.1,col.quali='#703A6D',label =c('quali'))
res.HCPC4<-HCPC(res.cluster4,nb.clust=6,consol=FALSE,graph=FALSE)
plot.HCPC(res.HCPC4,choice='tree',title='Arbre hiérarchique')
fviz_cluster(res.HCPC4,
repel = TRUE,
show.clust.cent = TRUE,
palette = "jco",
ggtheme = theme_minimal(),
main = "")
Pays à sélectionner : France , Allemagne , Royaume Uni , Hollande , Émirats Arabes Unis.
Stratégie européenne : Allemagne, RU, et Hollande.
# TESTE STATISTIQUE -----------
plot_num(df9completed)
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
# La variable Dispo_kcal_p_j a la semble suivre la forme d'une Loi Normale
hist(df9completed$Dispo_kcal_p_j)
# Vérifions notre variable candidate à l'aide d'un test statistique
##A Kolmogorov-Smirnoc n >50-60 --------
# Dispo_kcal_p_j
ks.test(df9completed$Dispo_kcal_p_j,"pnorm",mean=mean(df9completed$Dispo_kcal_p_j),sd=sd(df9completed$Dispo_kcal_p_j))
## Warning in ks.test(df9completed$Dispo_kcal_p_j, "pnorm", mean =
## mean(df9completed$Dispo_kcal_p_j), : ties should not be present for the
## Kolmogorov-Smirnov test
##
## One-sample Kolmogorov-Smirnov test
##
## data: df9completed$Dispo_kcal_p_j
## D = 0.041145, p-value = 0.9329
## alternative hypothesis: two-sided
# p-value = 0.9329 soit 93% de chance de rejetter l'hypothèse H0 à tort
# On peut conclure que la variable Dispo_kcal_p_j suit une Loi Normale
## B Test de comparaison entre deux clusters -------------
# 2 échantillons Var = cluster 3 et 4
E1 <- cbind(df9completed,liste) %>% select(pays,cluster,Dispo_kcal_p_j) %>% filter(cluster == 4)
E2 <- cbind(df9completed,liste) %>% select(pays,cluster,Dispo_kcal_p_j) %>% filter(cluster == 3)
# Egalité des variance ?
var.test(E2$Dispo_kcal_p_j,E1$Dispo_kcal_p_j)
##
## F test to compare two variances
##
## data: E2$Dispo_kcal_p_j and E1$Dispo_kcal_p_j
## F = 1.6366, num df = 74, denom df = 10, p-value = 0.4006
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.5152145 3.6443216
## sample estimates:
## ratio of variances
## 1.636555
# P-value = 0.4006 On ne rejette donc pas Ho au seuil de 0.05 %
# Egalité des moyennes ?
t.test(E2$Dispo_kcal_p_j,E1$Dispo_kcal_p_j,var.equal = T)
##
## Two Sample t-test
##
## data: E2$Dispo_kcal_p_j and E1$Dispo_kcal_p_j
## t = -0.98195, df = 84, p-value = 0.3289
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -284.83417 96.52387
## sample estimates:
## mean of x mean of y
## 2965.027 3059.182
# p-value = 0.2562 On ne rejette pas H0 au seuil de 0.05 %
# On vérifie notre résulat à l'aide du teste de Wilcoxon-Mann-Whitney :
wilcox.test(E2$Dispo_kcal_p_j,E1$Dispo_kcal_p_j, paired = F)
##
## Wilcoxon rank sum test with continuity correction
##
## data: E2$Dispo_kcal_p_j and E1$Dispo_kcal_p_j
## W = 327, p-value = 0.2717
## alternative hypothesis: true location shift is not equal to 0
# p-value = 0.2717 On ne rejette pas H0
# conclusion : 2 échantillons suivent une même loi normale