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.

Valeurs manquantes

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.

Estimation des NA

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.

ACP

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 :

  1. Les pays à forte croissance démographique ont une consommation calorique plus faible. différence Nord/Sud.
  2. Les plus gros consommateurs de volaille sont probablement des pays exportateurs : forte corrélation entre les deux variables
  3. Les pays importateurs sont moins représentatifs et sûrement en nombre plus restreint.

Clustering

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)

Cluster 4 : Pays importateurs

# 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.

Test statistique

# 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