Faire une ACP avec R

Analyse factorielle · R Lab — Antoine Géré

Exercice 1 (guidé) — ACP de la fertilité des sols de quinze sites du Cameroun

Un projet d'appui à la filière maïs veut savoir quels types de sols on rencontre entre la zone forestière du Sud et les savanes de l'Adamaoua, et si ces types de sols expliquent les différences de rendement. Sur 15 sites, on a analysé un échantillon composite de sol (horizon 0–20 cm). Le fichier data_stat-0109.xlsx contient, dans la feuille donnees, une ligne par site :

• huit variables d'analyse de sol : pH (pH eau), MO (matière organique, %), N (azote total, g/kg), P (phosphore assimilable, mg/kg), K (potassium échangeable, cmol⁺/kg), CEC (capacité d'échange cationique, cmol⁺/kg), Argile et Sable (%) ;
• Rdt_mais : rendement moyen du maïs sur le site (t/ha) ;
• Zone : zone agro-écologique (Forêt, Hauts plateaux, Transition, Savane).

La feuille description rappelle le sens et l'unité de chaque colonne (les données sont fictives).

1. Préparer la séance.
   1. Dans votre dossier de travail OneDrive, placer le fichier data_stat-0109.xlsx (encadré « Récupérer les données » ci-dessus) dans data_raw. Lancer RStudio à partir du fichier _stat_Istom.Rproj.
   2. Installer les packages readxl et FactoMineR avec RStudio, puis les charger (à chaque séance) :

      library(readxl)       # pour lire les fichiers Excel
      library(FactoMineR)   # pour l'ACP

   3. Importer les données :

      sols <- as.data.frame(read_excel("data_raw/data_stat-0109.xlsx", sheet = "donnees"))
      rownames(sols) <- sols$Site       # les sites deviennent les noms de lignes
      sols$Site <- NULL                 # on supprime la colonne devenue inutile
      sols$Zone <- as.factor(sols$Zone) # la zone est une variable qualitative (facteur)
      str(sols)                         # structure du tableau
      head(sols)                        # les 6 premieres lignes

      ✓

      Contrôlestr affiche 15 obs. of 10 variables ; les neuf premières colonnes sont num (numériques) et Zone est un Factor w/ 4 levels.
      Pourquoi faut-il mettre les noms des sites en noms de lignes, et non les laisser dans une colonne ? Quelle est maintenant la position (le numéro de colonne) de Rdt_mais ? de Zone ?
2. Regarder les données avant l'ACP.
   1. Afficher un résumé de chaque variable et les écarts-types :

      summary(sols)
      round(apply(sols[, 1:9], 2, sd), 2)   # ecart-type de chacune des colonnes 1 a 9

      Lecture

      apply(tableau, 2, fonction) applique la fonction à chaque colonne (le 2 désigne les colonnes) ; round(..., 2) arrondit à deux décimales. Attention : sd divise par \(I-1\), alors que l'ACP utilise des variances « population » (division par \(I\)) ; la différence est ici faible.

      Les variables ont-elles la même unité ? le même ordre de grandeur de dispersion ? En déduire qu'il faut faire une ACP normée (sur données centrées réduites). Que se passerait-il sinon ?
   2. Afficher la matrice des corrélations des huit variables de sol :

      round(cor(sols[, 1:8]), 2)

      Repérer les quatre couples de variables les plus corrélés (en valeur absolue). Ces liaisons vous semblent-elles logiques pour un agronome ?
   3. Regarder le rendement selon la zone :

      boxplot(Rdt_mais ~ Zone, data = sols, ylab = "Rendement (t/ha)")
      tapply(sols$Rdt_mais, sols$Zone, mean)   # rendement moyen par zone

      Quelle zone a les meilleurs rendements ? Gardez ce résultat en tête : on cherchera à l'expliquer par les sols à la question 7.
3. Lancer l'ACP.
   1. On veut construire les axes avec les huit variables de sol seulement, et utiliser le rendement et la zone pour interpréter. Taper :

      res <- PCA(sols, scale.unit = TRUE, quanti.sup = 9, quali.sup = 10,
                 graph = FALSE)

      Expliquer le rôle de chacun des arguments scale.unit = TRUE, quanti.sup = 9, quali.sup = 10 et graph = FALSE.
   2. Pourquoi est-il judicieux de mettre le rendement en variable supplémentaire plutôt qu'en variable active ?
   3. Que contient res ? Taper names(res). L'objet res est une liste qui range tous les résultats de l'ACP ; on y accède avec le signe $. Les notations sont celles du cours et des exercices faits à la main :

      Dans R	Ce que c'est
      res$eig	valeurs propres \(\lambda_s\), % d'inertie, % cumulé
      res$var$coord (ou $cor)	coordonnées des variables : \(r(X_k,F_s)\)
      res$var$cos2, res$var$contrib	\(\cos^2_s(X_k)=r^2(X_k,F_s)\) et contributions (%)
      res$ind$coord	coordonnées des individus : \(F_s(i)\)
      res$ind$cos2, res$ind$contrib	qualité de représentation et contributions (%)
      res$ind$dist	distance au centre de gravité \(d(i,G)\)
      res$quanti.sup, res$quali.sup	résultats des éléments supplémentaires

      Dans ces tableaux, chaque ligne est un élément (variable ou site) et chaque colonne un axe (Dim.1, Dim.2, …). Ainsi res$ind$coord[, 1:2] donne les coordonnées de tous les sites sur les axes 1 et 2, et res$ind$coord["Meiganga", 1] la coordonnée de Meiganga sur l'axe 1. Que donne summary(res) ?
4. Combien d'axes faut-il regarder ?
   Rappels

   En ACP normée, chacune des \(K\) variables actives a une variance égale à \(1\) : l'inertie totale vaut \(K\). Le critère de Kaiser retient les axes tels que \(\lambda_s\ge 1\), c'est-à-dire ceux qui résument au moins autant qu'une variable de départ (ici, avec \(K=8\), au moins \(\tfrac18 = 12{,}5\,\%\) de l'inertie). Le critère du coude retient les axes situés avant la cassure de l'éboulis, là où les barres deviennent petites et à peu près de même hauteur.

   1. Afficher les valeurs propres et l'éboulis :

      round(res$eig, 2)
      barplot(res$eig[, 2], names.arg = 1:nrow(res$eig),
              ylab = "% d'inertie", xlab = "axe")

      Lecture

      La première colonne de res$eig contient les valeurs propres, la deuxième le pourcentage d'inertie, la troisième le pourcentage cumulé.

      ✓

      Contrôle\(\lambda_1 \approx 4{,}61\).
      Que vaut la somme des valeurs propres (sum(res$eig[, 1])) ? Pourquoi ?
   2. Quel pourcentage d'inertie le plan \((1,2)\) restitue-t-il ?
   3. Combien d'axes retenez-vous selon le critère de Kaiser ? selon le critère du coude ?
5. Le cercle des corrélations : interpréter les axes.
   1. Tracer le cercle et afficher les résultats des variables :

      plot(res, choix = "var")          # cercle des correlations
      round(res$var$cor[, 1:2], 2)      # coordonnees = correlations avec les axes
      round(res$var$cos2[, 1:2], 2)     # qualite de representation
      round(res$var$contrib[, 1:2], 1)  # contributions (en %)

      À noter

      Sur le cercle, la flèche en pointillés est celle du rendement, variable supplémentaire ; on l'utilisera à la question 7.

      Tracer ensuite le cercle en coloriant les flèches selon leur qualité de représentation, puis selon leur contribution :

      plot(res, choix = "var", habillage = "cos2")     # couleur = cos2 sur le plan
      plot(res, choix = "var", habillage = "contrib")  # couleur = contribution au plan

      Lecture

      La couleur va du bleu (valeur faible) au rouge (valeur forte), selon l'échelle affichée à droite. Avec "cos2", c'est la qualité de représentation sur le plan, \(\cos^2_1+\cos^2_2\). Avec "contrib", c'est la contribution au plan \((1,2)\) : la moyenne des contributions aux axes 1 et 2, pondérée par \(\lambda_1\) et \(\lambda_2\).

   2. Axe 1. Quelles variables lui sont fortement corrélées (en valeur absolue, disons au-dessus de \(0{,}8\)) ? Dans quel sens ? Lesquelles y contribuent le plus (comparer à la contribution moyenne \(100/8 = 12{,}5\,\%\)) ? Proposer un nom pour cet axe.
   3. Axe 2. Mêmes questions.
   4. Sur le cercle colorié selon le \(\cos^2\), quelle variable est la moins bien représentée sur le plan \((1,2)\) ? Calculer sa qualité de représentation \(\cos^2_1+\cos^2_2\), ou l'obtenir directement avec round(rowSums(res$var$cos2[, 1:2]), 2). Sur quel axe se trouve le reste de son information ? (Regarder round(res$var$cos2[, 3], 2).)
   5. Sur le cercle colorié selon la contribution, lire l'échelle de couleurs : les contributions au plan sont-elles très différentes d'une variable à l'autre ? Comparer avec les contributions axe par axe affichées en a). Pourquoi vaut-il mieux lire les contributions axe par axe pour nommer les axes ?
   6. Faire le lien avec le calcul à la main. Pour la variable MO, vérifier avec les valeurs affichées que :
      • \(\cos^2_1(\text{MO}) = r^2(\text{MO},F_1)\) ;
      • \(\operatorname{contrib}_1(\text{MO}) = 100\times r^2(\text{MO},F_1)/\lambda_1\) ;
      • \(r(\text{MO},F_1) = \sqrt{\lambda_1}\,u_{1k}\), où \(u_{1k}\) est la coordonnée de MO dans le premier vecteur propre \(U_1\), donnée par res$svd$V[, 1] (MO est la 2ᵉ variable).
   7. Sur le cercle, les flèches de P et K sont presque confondues. Calculer le produit scalaire de leurs coordonnées sur le plan (\(r(\text{P},F_1)\,r(\text{K},F_1)+r(\text{P},F_2)\,r(\text{K},F_2)\)) et le comparer à leur corrélation \(r(\text{P},\text{K})\), calculée à la question 2b. Pourquoi les deux valeurs sont-elles proches sans être égales ?
6. La carte des sites.
   1. Tracer la carte des sites coloriée selon la zone, et afficher les résultats des sites :

      plot(res, choix = "ind", habillage = 10)   # couleur = zone (10e colonne)
      round(res$ind$coord[, 1:2], 2)             # coordonnees F1, F2
      round(res$ind$cos2[, 1:2], 2)              # qualite de representation
      round(res$ind$contrib[, 1:2], 1)           # contributions (en %)

      À noter

      Sur la carte, les petits carrés dont le nom est écrit en italique sont les barycentres des zones (la position moyenne des sites de chaque zone).

   2. Décrire la position des quatre zones sur la carte. En vous aidant de la question 5, caractériser les sols de chaque zone en une phrase.
   3. Quels sont les deux sites qui contribuent le plus à l'axe 1 ? à l'axe 2 ? Comparer à la contribution moyenne \(100/15 \approx 6{,}7\,\%\). Un site domine-t-il à lui seul un axe ? Tracer ensuite la carte coloriée selon la contribution au plan :

      plot(res, choix = "ind", habillage = "contrib")  # couleur = contribution au plan

      Les sites en rouge sont-ils bien ceux que vous avez repérés ? Où se trouvent-ils sur la carte ?
   4. Afficher la qualité de représentation de chaque site sur le plan, tracer la carte coloriée selon le \(\cos^2\), puis ne tracer que les sites bien représentés :

      round(rowSums(res$ind$cos2[, 1:2]), 2)            # cos2 sur le plan (1,2)
      plot(res, choix = "ind", habillage = "cos2")      # couleur = cos2 sur le plan
      plot(res, choix = "ind", select = "cos2 0.6")     # sites avec cos2 > 0.6

      Quels sites disparaissent de la carte ? Quel est le plus mal représenté ? Peut-on interpréter sa position près du centre ?
   5. Faire le lien avec le calcul à la main. Pour Meiganga, vérifier que \[ \operatorname{contrib}_1 = 100\times\frac{F_1(i)^2}{I\,\lambda_1} \qquad\text{et}\qquad \cos^2_1 = \frac{F_1(i)^2}{d^2(i,G)}, \] avec res$ind$coord["Meiganga", 1], res$eig[1, 1] et res$ind$dist["Meiganga"].
   6. Retour aux données brutes. Afficher les lignes de trois sites situés à trois coins différents de la carte (par exemple Mbouda, Meiganga, Ebolowa) avec

      sols[c("Mbouda", "Meiganga", "Ebolowa"), ]

      et vérifier que leurs valeurs sont cohérentes avec les noms donnés aux axes.
7. Les éléments supplémentaires : et le rendement ?
   1. Afficher la corrélation du rendement avec les axes :

      round(res$quanti.sup$cor[, 1:2], 2)
      round(res$quanti.sup$cos2[, 1:2], 2)

      Où se trouve la flèche du rendement sur le cercle ? Est-elle bien représentée ? Quels sols semblent les plus favorables au maïs ? Est-ce cohérent avec la boîte à moustaches de la question 2 ?
   2. Afficher les coordonnées des barycentres des zones et leurs valeurs-test :

      round(res$quali.sup$coord[, 1:2], 2)
      round(res$quali.sup$v.test[, 1:2], 2)

      Lecture

      Une valeur-test supérieure à \(2\) en valeur absolue indique que les sites de la zone ont, sur cet axe, une position moyenne significativement différente de \(0\) (au seuil d'environ \(5\,\%\)).

      Quelles zones sont significativement d'un côté ou de l'autre de chaque axe ?
   3. Obtenir la description automatique des axes :

      dimdesc(res, axes = 1:2)

      Lecture

      Pour chaque axe, $quanti liste les variables significativement corrélées à l'axe (de la plus positive à la plus négative) ; $quali donne le \(R^2\) de la zone, c'est-à-dire la part de la variabilité des coordonnées sur l'axe expliquée par la zone, avec sa \(p\)-valeur ; $category liste les zones significatives.

      Le rendement est-il significativement lié à un axe ? Quelle part de la variabilité de l'axe 1 la zone explique-t-elle ? Que peut-on en conclure ?
8. Synthèse. Rédiger, pour le chef de projet, une synthèse de l'ACP : les deux grands gradients qui structurent les sols, les zones qui s'opposent, le lien avec le rendement et ce que cela suggère pour la fertilisation du maïs. Préciser les limites de l'étude.

Exercice 2 (en autonomie) — typologie d'exploitations familiales au Burkina Faso

Dans l'Ouest du Burkina Faso, une fédération de coopératives de producteurs veut adapter les services qu'elle propose à ses membres (conseil agricole, crédit de campagne, approvisionnement en intrants, commercialisation). Ses conseillers agricoles ont conduit une enquête auprès de 32 exploitations familiales de trois villages : Koumbia (zone cotonnière), Bondoukuy et Dano. Avant de discuter des services avec les producteurs, ils veulent savoir quels grands types d'exploitations coexistent dans la zone et ce qui les distingue.

Le fichier data_stat-0110.xlsx (feuille donnees) contient une ligne par exploitation et les colonnes suivantes ; la feuille description précise les unités. Les données sont fictives, construites pour l'exercice.

1. Préparer les données. Importer le fichier, mettre les identifiants en noms de lignes, déclarer correctement le type de chaque variable et faire un premier résumé du tableau.
2. Choisir les variables. Quelles variables mettez-vous en actives (celles qui définissent la ressemblance entre exploitations) ? Lesquelles en supplémentaires ? Justifier en particulier le statut du revenu. Faut-il réduire les variables ?
3. Réaliser l'ACP et choisir le nombre d'axes à interpréter, en croisant au moins deux critères.
4. Interpréter les axes retenus à l'aide du cercle des corrélations, des qualités de représentation, des contributions des variables et de la description automatique des axes. Donner un nom à chaque axe. Une variable est-elle liée à plusieurs axes à la fois ?
5. Décrire la carte des exploitations. Situer les trois villages, l'accès au crédit et l'équipement en traction animale : lesquelles de ces variables sont liées à quel axe ? Pour chaque grand type d'exploitation que vous identifiez, citer une ou deux exploitations représentatives et vérifier leur profil sur les données brutes.
6. Le revenu. Quels types d'exploitations ont les revenus les plus élevés ? les plus faibles ? La structure de l'exploitation suffit-elle à expliquer le revenu ?
7. Contrôler la fiabilité de la lecture. Repérer les exploitations mal représentées sur le plan et celles qui pèsent le plus dans la construction des axes. Ces dernières posent-elles problème ? Le vérifier.
8. Note pour la fédération. En une demi-page, à destination du conseil d'administration de la fédération, présenter la typologie obtenue et proposer, pour chaque type d'exploitation, des services adaptés à discuter avec les producteurs concernés. Indiquer les limites de l'analyse.

Indications

Exercice 1

1. Si read_excel ne trouve pas le fichier, vérifier qu'il est bien dans le dossier data_raw, que RStudio a été lancé depuis le fichier .Rproj (getwd() affiche le dossier de travail) et que le chemin s'écrit avec une barre oblique /, jamais avec \. Les colonnes sont numérotées après la suppression de Site.
2. Comparer les écarts-types de K et de Sable. Les corrélations les plus fortes dépassent \(0{,}8\) en valeur absolue.
3. Une variable supplémentaire est projetée après le calcul des axes : elle ne les influence pas.
4. L'inertie totale d'une ACP normée vaut le nombre de variables actives.
5. Pour nommer un axe, partir des variables de corrélation supérieure à \(0{,}8\) en valeur absolue, en regardant leur signe. La qualité de représentation d'une variable sur le plan est la somme de ses \(\cos^2\) sur les axes 1 et 2.
6. La contribution moyenne d'un site vaut \(100/15\). res$ind$dist donne la distance \(d(i,G)\), pas son carré.
7. Une variable qualitative supplémentaire se lit sur la carte des sites (barycentres des zones), pas sur le cercle.

Exercice 2

1. Les trois variables oui/non et le village doivent être des facteurs.
2. Une variable que l'on cherche à expliquer (un résultat économique, par exemple) se met plutôt en supplémentaire : on construit la typologie sans elle, puis on regarde où elle se place. Comparer les unités et les écarts-types des variables.
3. Comparer les valeurs propres à \(1\) et regarder l'éboulis.
4. Une flèche « en diagonale » est liée aux deux axes à la fois. Un axe ne classe pas forcément les exploitations de la meilleure à la moins bonne : il peut opposer deux façons de produire.
5. Pour lier une variable qualitative à un axe, lire la rubrique $quali de dimdesc, ou les valeurs-test (v.test) des modalités.
6. La fonction tapply calcule une moyenne par groupe, par exemple le revenu moyen selon l'accès au crédit.
7. Avec 32 exploitations, la contribution moyenne vaut \(100/32 \approx 3{,}1\,\%\). L'argument ind.sup de PCA permet de refaire l'analyse sans quelques exploitations.

Mise en forme par Antoine Géré.