Correction — Faire une ACP avec R

Exercice 1 (guidé) · fertilité des sols, Cameroun — Antoine Géré

1. Préparer la séance.

   c) str(sols) affiche 'data.frame': 15 obs. of 10 variables : neuf colonnes num (de pH à Rdt_mais) et Zone : Factor w/ 4 levels.

   • Les noms des sites sont des étiquettes, pas une variable mesurée : en noms de lignes, ils servent à nommer les points de la carte, et PCA ne les traite pas comme une donnée (une colonne de texte ferait d'ailleurs échouer la fonction).
   • Après suppression de Site, les colonnes 1 à 8 (pH, MO, N, P, K, CEC, Argile, Sable) sont les variables actives, Rdt_mais est la colonne 9 (quanti.sup) et Zone la colonne 10 (quali.sup) : ce sont ces numéros qu'on utilisera dans PCA.
2. Regarder les données avant l'ACP.
   1. Unités et dispersions. Les unités sont toutes différentes (pH sans unité, %, g/kg, mg/kg, cmol⁺/kg) et les dispersions très inégales :

      	pH	MO	N	P	K	CEC	Argile	Sable	Rdt
      Moyenne	5,67	3,25	1,75	11,75	0,42	11,93	30,6	47,0	2,55
      Écart-type	0,47	1,14	0,66	4,41	0,16	4,67	11,10	12,01	0,68
      Variance	0,22	1,29	0,43	19,43	0,03	21,79	123,26	144,14	—

      Sans réduction, l'inertie totale serait la somme des variances (\(\approx 310\)), et chaque variable y pèserait selon sa variance : l'argile et le sable feraient à eux seuls 86 % de l'inertie, le pH ou le potassium presque rien — simplement à cause de leur unité. L'ACP normée donne à chaque variable une variance égale à 1, donc le même poids (\(100/8 = 12{,}5\,\%\)) : elle est indispensable ici.

      Part de l'inertie sans réduction	pH	MO	N	P	K	CEC	Argile	Sable
      %	0,07	0,42	0,14	6,26	0,01	7,02	39,69	46,41

   2. Corrélations. La matrice se lit mieux en couleurs (rouge : liaison positive, bleu : négative ; plus foncé, plus fort). Les quatre couples les plus corrélés sont surlignés :

      	pH	MO	N	P	K	CEC	Argile	Sable
      pH	1	−0,59	−0,49	0,43	0,33	−0,23	−0,57	0,65
      MO	−0,59	1	0,95	0,28	0,30	0,76	0,72	−0,88
      N	−0,49	0,95	1	0,36	0,37	0,78	0,74	−0,87
      P	0,43	0,28	0,36	1	0,82	0,69	0,12	−0,07
      K	0,33	0,30	0,37	0,82	1	0,63	0,04	−0,01
      CEC	−0,23	0,76	0,78	0,69	0,63	1	0,57	−0,63
      Argile	−0,57	0,72	0,74	0,12	0,04	0,57	1	−0,70
      Sable	0,65	−0,88	−0,87	−0,07	−0,01	−0,63	−0,70	1

      Couple	\(r\)	Explication agronomique
      MO – N	0,95	l'azote du sol est surtout organique : MO et N vont ensemble
      MO – Sable	−0,88	les sols sableux retiennent mal la matière organique
      N – Sable	−0,87	conséquence des deux lignes précédentes
      P – K	0,82	P et K sont tous deux abondants dans les sols volcaniques

      Ces liaisons annoncent déjà deux « familles » de variables (matière organique et texture d'un côté, P et K de l'autre), que l'ACP va retrouver.
   3. Rendement selon la zone. Rendements moyens : Hauts plateaux 3,50 t/ha, Forêt 2,34, Transition 2,17, Savane 2,00. Les Hauts plateaux se détachent nettement (tous leurs sites au-dessus de 3 t/ha) ; les trois autres zones se chevauchent largement autour de 2 à 2,5 t/ha. On cherchera à l'expliquer par les sols.
3. Lancer l'ACP.
   1. Rôle des arguments.

      Argument	Rôle
      scale.unit = TRUE	ACP normée : chaque variable est centrée et réduite
      quanti.sup = 9	la colonne 9 (Rdt_mais) est une variable quantitative supplémentaire
      quali.sup = 10	la colonne 10 (Zone) est une variable qualitative supplémentaire
      graph = FALSE	pas de graphique automatique : on les tracera nous-mêmes, un par un

   2. Pourquoi le rendement en supplémentaire ? La question est « les types de sols expliquent-ils le rendement ? ». Il faut donc construire la typologie des sols sans le rendement, puis regarder où il se place. En variable active, le rendement participerait à la définition des axes et le raisonnement deviendrait circulaire ; il dépend en outre d'autres facteurs que le sol (pluie, variété, conduite). Les éléments supplémentaires ne modifient pas les axes : ils sont placés ensuite, pour aider à les interpréter.
   3. Que contient res ? names(res) affiche : "eig" "var" "ind" "svd" "quanti.sup" "quali.sup" "call". summary(res) résume toute l'ACP : valeurs propres, puis coordonnées (Dim.1…), contributions (ctr) et \(\cos^2\) des premiers individus et variables, et résultats des éléments supplémentaires.
4. Combien d'axes ?
   1. Valeurs propres.

      Axe	1	2	3	4	5	6	7	8
      \(\lambda_s\)	4,61	2,36	0,38	0,29	0,19	0,08	0,06	0,03
      % inertie	57,7	29,5	4,7	3,6	2,4	1,0	0,7	0,4
      % cumulé	57,7	87,2	91,9	95,6	97,9	98,9	99,6	100

      ✓

      Contrôle\(\lambda_1 \approx 4{,}61\) ; le pourcentage est \(100\times\lambda_s/8\) (ex. \(100\times 4{,}613/8 = 57{,}7\,\%\)).
      sum(res$eig[, 1]) vaut 8 : en ACP normée, l'inertie totale est égale au nombre de variables actives.
   2. Le plan \((1,2)\) restitue \(57{,}7+29{,}5 = \mathbf{87{,}2\,\%}\) de l'inertie : la carte sera fidèle pour la plupart des sites.
   3. Nombre d'axes.
      • Kaiser : \(\lambda_1 = 4{,}61\) et \(\lambda_2 = 2{,}36\) dépassent 1, \(\lambda_3 = 0{,}38\) non : deux axes.
      • Coude : la chute est forte jusqu'à l'axe 2, puis les barres deviennent petites et presque égales : coude très net après le deuxième axe, même conclusion.
      On interprète donc le plan \((1,2)\).
5. Le cercle des corrélations.
   1. Résultats des variables (en gras : contributions supérieures à la moyenne \(12{,}5\,\%\)) :

      Variable	\(r(X_k,F_1)\)	\(r(X_k,F_2)\)	\(\cos^2_1\)	\(\cos^2_2\)	plan	contrib 1	contrib 2
      pH	−0,54	0,73	0,30	0,54	0,83	6,4	22,7
      MO	0,96	−0,10	0,92	0,01	0,93	19,9	0,4
      N	0,97	−0,01	0,93	0,00	0,93	20,2	0,0
      P	0,40	0,87	0,16	0,77	0,93	3,5	32,4
      K	0,38	0,85	0,15	0,72	0,86	3,2	30,4
      CEC	0,86	0,39	0,75	0,15	0,90	16,2	6,4
      Argile	0,80	−0,27	0,64	0,07	0,71	13,8	3,1
      Sable	−0,88	0,33	0,77	0,11	0,88	16,8	4,6
      Rdt (sup.)	0,59	0,47	0,35	0,22	0,57	—	—

      ✓

      Contrôlechaque colonne de \(\cos^2\) totalise la valeur propre de l'axe (\(4{,}61\) ; \(2{,}36\)), et chaque colonne de contributions totalise 100 %. Le rendement, supplémentaire, n'a pas de contribution.

   2. Axe 1 (57,7 %). Corrélations fortes et positives avec N (0,97), MO (0,96), CEC (0,86) et Argile (0,80) ; forte et négative avec Sable (−0,88). Ce sont aussi les cinq variables qui y contribuent plus que la moyenne (87 % de l'axe à elles cinq). Il oppose les sols argileux, riches en matière organique et à forte CEC (à droite) aux sols sableux et pauvres (à gauche) : c'est un axe « réserve organique et texture ».
   3. Axe 2 (29,5 %). Corrélations fortes et positives avec P (0,87), K (0,85) et, dans une moindre mesure, pH (0,73) ; ces trois variables font 85,5 % de l'axe. Il oppose, à niveau organique donné, des sols peu acides et riches en P et K (en haut) à des sols acides et pauvres (en bas) : c'est un axe « acidité et richesse minérale », indépendant du premier.
   4. Qualité de représentation. L'argile est la moins bien représentée sur le plan : \(\cos^2 = 0{,}64+0{,}07 = \mathbf{0{,}71}\) ; MO, N et P sont les mieux représentées (0,93). Le reste de l'information de l'argile est surtout sur l'axe 3 (\(\cos^2_3 = 0{,}26\)).
   5. Le piège des couleurs « contribution ». La contribution au plan est la moyenne des contributions aux deux axes, pondérée par les valeurs propres ; pour MO : \[ \frac{4{,}613\times 19{,}93 + 2{,}363\times 0{,}42}{4{,}613+2{,}363} \approx 13{,}3\,\%. \] Au plan, toutes les variables sont proches de 12,5 % (de 10,2 % pour l'argile à 13,4 % pour N). Mais axe par axe, c'est tranché : MO et N font 40 % de l'axe 1 et presque rien sur l'axe 2 ; P et K font 63 % de l'axe 2 et presque rien sur l'axe 1. Pour nommer un axe, on lit donc les contributions axe par axe.
   6. Lien avec le calcul à la main (MO). Avec \(r(\text{MO},F_1)=0{,}959\), \(\lambda_1=4{,}613\) et \(u_{12}=0{,}446\) : \[ \cos^2_1 = 0{,}959^2 = 0{,}919;\quad \operatorname{contrib}_1 = 100\times\tfrac{0{,}919}{4{,}613} = 19{,}9\,\%;\quad r = \sqrt{4{,}613}\times 0{,}446 = 0{,}958. \] Le logiciel applique exactement les formules vues à la main.
   7. P et K. Produit scalaire des deux flèches sur le plan : \(0{,}400\times0{,}382 + 0{,}875\times0{,}847 = \mathbf{0{,}89}\), contre \(r(\text{P},\text{K}) = \mathbf{0{,}82}\). La corrélation est la somme de ces produits sur tous les axes ; les axes suivants apportent de petites corrections négatives. Les deux valeurs sont proches (flèches longues) mais pas égales (une petite part de l'information de P et K est hors du plan).
6. La carte des sites.
   1. Carte des sites, couleur = zone ; les losanges sont les barycentres des zones.

   2. Les quatre zones (moyennes brutes par zone, en gras les écarts marqués à la moyenne) :

      Zone	barycentre	pH	MO	P	K	Sable	Rdt
      Forêt	(1,28 ; −1,34)	5,18	4,06	9,48	0,36	38,8	2,34
      Hauts plateaux	(1,97 ; 1,66)	5,80	4,03	17,60	0,60	40,5	3,50
      Transition	(−1,56 ; −0,69)	5,87	2,60	8,57	0,30	51,7	2,17
      Savane	(−3,20 ; 0,71)	6,10	1,53	10,93	0,42	64,7	2,00
      Moyenne	(0 ; 0)	5,67	3,25	11,75	0,42	47,0	2,55

      • Forêt (bas-droite) : sols argileux et organiques, mais acides et pauvres en P et K (ferrallitiques désaturés).
      • Hauts plateaux (haut-droite) : organiques et riches en P et K, peu acides (sols volcaniques).
      • Savane (gauche) : sableux et pauvres en matière organique, mais peu acides.
      • Transition (gauche, un peu en bas) : position intermédiaire, plutôt sableux et pauvres.
   3. Contributions des sites (moyenne \(100/15 \approx 6{,}7\,\%\)). Axe 1 : Meiganga (21,3 %) et Ngaoundéré (13,8 %) — les sols les plus sableux et pauvres. Axe 2 : Mbouda (22,1 %) et Bafoussam (15,7 %), suivis d'Ebolowa (15,0 %). Aucun site ne domine à lui seul un axe : les deux gradients sont portés par plusieurs sites, tous loin du centre. La contribution dépend du carré de la distance au centre.
   4. Qualité de représentation. Avec le seuil \(\cos^2 > 0{,}6\) de select = "cos2 0.6", trois sites disparaissent : Foumban (0,21), Akonolinga (0,39) et Sangmélima (0,56), tous près du centre. Foumban est le plus mal représenté : sa position centrale ne doit pas s'interpréter comme « site moyen » (79 % du carré de sa distance est sur les axes suivants, dont \(\cos^2_3 = 0{,}22\)).
   5. Lien avec le calcul à la main (Meiganga). Avec \(F_1=-3{,}840\), \(\lambda_1=4{,}613\), \(d(i,G)=4{,}006\) : \[ \operatorname{contrib}_1 = 100\times\frac{(-3{,}840)^2}{15\times 4{,}613} \approx 21{,}3\,\%;\qquad \cos^2_1 = \frac{14{,}75}{16{,}05} \approx 0{,}92. \] On retrouve les valeurs de R.
   6. Retour aux données brutes (bleu : nettement au-dessus de la moyenne, rouge : nettement au-dessous) :

      Site	position	pH	MO	N	P	K	CEC	Argile	Sable
      Mbouda	haut-droite	6,2	4,2	2,57	21,1	0,73	17,3	44	40
      Meiganga	gauche	6,4	0,8	0,73	10,8	0,29	5,8	21	65
      Ebolowa	bas-droite	4,9	4,5	2,50	7,3	0,30	14,3	46	33
      Moyenne		5,67	3,25	1,75	11,75	0,42	11,93	30,6	47,0

      Mbouda est organique et très riche en P et K ; Meiganga très sableuse et pauvre ; Ebolowa organique et argileuse mais acide et pauvre en P. Les données brutes confirment les noms des axes.
7. Les éléments supplémentaires.
   1. Le rendement. \(r(\text{Rdt},F_1)=0{,}59\) et \(r(\text{Rdt},F_2)=0{,}47\) : la flèche pointe vers le haut-droite, du côté des Hauts plateaux. Les meilleurs rendements s'observent sur les sols à la fois organiques et riches en P et K, peu acides (3,2 à 3,6 t/ha à Bafoussam, Dschang, Mbouda, Foumban). Les sols de forêt, pourtant organiques, plafonnent (acidité, carence en P) ; savane et transition sont pauvres. C'est cohérent avec la boîte à moustaches. La flèche reste loin du cercle (\(\cos^2 = 0{,}57\)) : le sol n'explique qu'une partie du rendement.
   2. Barycentres et valeurs-test (|v.test| > 2 : significatif) :

      Zone	Dim 1	Dim 2	v.test 1	v.test 2
      Forêt	1,28	−1,34	1,58	−2,31
      Hauts plateaux	1,97	1,66	2,06	2,43
      Savane	−3,20	0,71	−2,79	0,86
      Transition	−1,56	−0,69	−1,36	−0,83

      Axe 1 : les Hauts plateaux sont significativement à droite, la Savane à gauche. Axe 2 : les Hauts plateaux significativement en haut, la Forêt en bas. La Transition n'est significative sur aucun axe (trois sites, position intermédiaire).
   3. Description automatique des axes (dimdesc) :

      	Axe 1	Axe 2
      $quanti	N (0,97), MO (0,96), CEC (0,86), Argile (0,80), Rdt (0,59, p = 0,02) ; pH (−0,54), Sable (−0,88)	P (0,87), K (0,85), pH (0,73)
      $quali (Zone)	\(R^2 = 0{,}89\), \(p \approx 10^{-5}\)	\(R^2 = 0{,}65\), \(p \approx 0{,}008\)
      $category	Hauts plateaux (+), Savane (−)	Hauts plateaux (+), Forêt (−)

      Le rendement n'est significativement lié qu'à l'axe 1 (sa corrélation avec l'axe 2, \(r = 0{,}47\), n'est pas significative avec 15 sites, \(p \approx 0{,}08\)). La zone explique 89 % de la variabilité des coordonnées sur l'axe 1 : la typologie des sols colle très étroitement à la zonation agro-écologique.
8. Synthèse.
   Exemple de synthèse

   Les sols des quinze sites s'organisent selon deux gradients qui résument 87 % de l'information. Le premier oppose les sols argileux et riches en matière organique de la forêt et des hauts plateaux aux sols sableux et pauvres des savanes de l'Adamaoua. Le second oppose les sols volcaniques des hauts plateaux, peu acides et riches en phosphore et potassium, aux sols acides et carencés en phosphore de la zone forestière. Le rendement du maïs suit surtout le premier gradient et il est maximal sur les hauts plateaux, qui cumulent les deux atouts. En zone forestière, la priorité serait la correction de l'acidité (chaulage) et l'apport de phosphore ; en savane, l'apport de matière organique et une fertilisation complète.

   Limites : quinze sites seulement, un seul prélèvement par site, des données fictives, et un rendement qui dépend aussi du climat et des pratiques ; l'ACP décrit des associations, elle ne prouve pas de causalité.

Mise en forme par Antoine Géré.