factorisation de ca en fonctions réutilisables

This commit is contained in:
Pierre-Edouard Portier 2023-01-08 16:54:44 +01:00
parent 59772d55d0
commit a509b95c5a
2 changed files with 197 additions and 165 deletions

View File

@ -7,10 +7,10 @@ source('05_c_svd_ca_code.R')
## Profils lignes et profils colonnes
Soit une matrice de données $\mathbf{N} \in \mathbb{R}^{I \times J}$. Elle peut représenter $I$ individus décrits par $J$ variables ($n_{i,j}$ est alors la valeur (qui doit être positive pour que les analyses à venir fassent sens) de la variable $j$ pour l'individu $i$) ou le croisement d'une première variable avec $I$ modalités et d'une seconde variable avec $J$ modalités (par exemple, couleurs des yeux et couleurs des cheveux d'un groupe de personnes, alors $n_{i,j}$ est le décompte du nombre de personnes ayant des yeux de la $i$-ème couleur et des cheveux de la $j$-ème couleur). Dans ce dernier cas, nous parlons de \emph{matrice de contingence}. Ce chapitre se concentre sur l'analyse des correspondances d'une matrice de contingence. Dans le chapitre suivant, nous montrerons comment analyser un jeu de données tabulaire multivarié (i.e. croisant individus et variables). Notons :
Soit une matrice de données $\mathbf{N} \in \mathbb{R}^{I \times J}$. Elle peut représenter $I$ individus décrits par $J$ variables. $n_{i,j}$ est alors la valeur, qui doit être positive pour que les analyses à venir soient justes, de la variable $j$ pour l'individu $i$. Elle peut aussi représenter le croisement d'une première variable catégorielle avec $I$ modalités et d'une seconde variable catégorielle avec $J$ modalités. Par exemple, si la première variable décrit la couleur des yeux et la seconde variable la couleur des cheveux d'un groupe de personnes, $n_{i,j}$ est alors le décompte du nombre de personnes ayant des yeux de la $i$-ème couleur et des cheveux de la $j$-ème couleur. Dans ce dernier cas, nous parlons de \emph{matrice de contingence}. Ce chapitre se concentre sur l'analyse des correspondances d'une matrice de contingence. Dans le chapitre suivant, nous montrerons comment analyser un jeu de données tabulaire multivarié (i.e. croisant individus et variables). Notons :
$$n \triangleq \sum_{i=1}^{I} \sum_{j=1}^{J} n_{i,j}$$
Pour fixer les idées, nous considérons un exemple développé par Jean-Paul Benzécri [@benzecri1985phosphates] dans un article qui introduit à l'analyse des correspondances sur la base de données du commerce mondial des phosphates (voir Table \@ref(tab:mat-N)). A l'intersection de la ligne $i$ et de la colonne $j$ du tableau de données, nous lisons, en milliers de tonnes de $P_2O_5$, la quantité importée par le pays $i$ depuis le pays $j$ durant 8 années, de 1973 à 1980. Les labels des pays importateurs et exportateurs sont donnés dans les Tables \@ref(tab:lblI) et \@ref(tab:lblJ).
Pour fixer les idées, nous considérons un exemple développé par Jean-Paul Benzécri [@benzecri1985phosphates] dans un article qui introduit à l'analyse des correspondances à partir de données du commerce mondial des phosphates (voir Table \@ref(tab:mat-N)). A l'intersection de la ligne $i$ et de la colonne $j$ du tableau de données, nous lisons, en milliers de tonnes de $P_2O_5$, la quantité importée par le pays $i$ depuis le pays $j$ durant 8 années, de 1973 à 1980. Les labels des pays importateurs et exportateurs sont donnés dans les Tables \@ref(tab:lblI) et \@ref(tab:lblJ).
```{r, include=FALSE}
#N <- matrix( c(1,2,1,3,2,5,1,7,3,6,2,8),
@ -195,7 +195,7 @@ Ainsi, l'inertie du nuage des profils lignes par rapport à son centre de masse
## GSVD des profils lignes et des profils colonnes
Nous souhaitons représenter les profils lignes dans un repère dont les axes orthonormaux sont orientés selon les direction d'inertie maximale. Nous nommerons ces axes les \emph{axes principaux} ou encore les \emph{facteurs}. Nous devons ainsi calculer le SVD :
Nous souhaitons représenter les profils lignes dans un repère dont les axes orthonormaux sont orientés selon les direction d'inertie maximale. Nous nommerons ces axes les \emph{axes principaux} ou encore les \emph{axes factoriels}. Nous devons ainsi calculer le SVD :
$$\mathbf{P} - \mathbf{r}\mathbf{c}^T = \mathbf{A}\boldsymbol{\Delta}\mathbf{B}^T \; \text{avec } \; \mathbf{A}^T\mathbf{D_r^{-1}}\mathbf{A} = \mathbf{B}^T\mathbf{D_c^{-1}}\mathbf{B} = \mathbf{I}$$
C'est une forme généralisée de la décomposition en valeurs singulières (GSVD) puisque $\mathbf{A}$ et $\mathbf{B}$ sont des matrices orthogonales dans des espaces euclidiens associés aux métriques respectives $\mathbf{D_r^{-1}}$ (pour prendre en compte les masses associées à chaque profil ligne) et $\mathbf{D_c^{-1}}$ (pour prendre en compte les échelles différentes de chaque colonne).
@ -209,7 +209,7 @@ Nous retrouvons la même inertie des écarts à l'indépendance que celle décou
## Calcul du GSVD
Par ailleurs, remarquons que l'inertie des écarts à l'indépendance, avec l'espace des colonnes $\mathbb{R}^J$ associé à la métrique $\mathbf{D_c^{-1}}$ et l'espace des lignes $\mathbb{R}^I$ associé à la métrique $\mathbf{D_r^{-1}}$, est aussi l'inertie des écarts à l'indépendance standardisés, avec les espaces des colonnes et des lignes associés à la métrique identité usuelle :
Par ailleurs, remarquons que l'inertie des écarts à l'indépendance, avec l'espace des colonnes $\mathbb{R}^J$ associé à la métrique $\mathbf{D_c^{-1}}$ et l'espace des lignes $\mathbb{R}^I$ associé à la métrique $\mathbf{D_r^{-1}}$, est aussi l'inertie des écarts à l'indépendance standardisés, avec les espaces des colonnes et des lignes associés à la métrique identité :
$$\sum_i \sum_j \frac{\left(P_{i,j} - r_ic_j\right)^2}{r_i c_j} \; = \; \sum_i \sum_j \left(\frac{P_{i,j} - r_ic_j}{\sqrt{r_i}\sqrt{c_j}}\right)^2$$
Donc, nous pouvons découvrir un sous-espace de projection optimal (au sens de la préservation de l'inertie), pour les profils lignes et pour les profils colonnes, en calculant le SVD (classique) :
@ -241,7 +241,7 @@ De même pour montrer $\mathbf{B}^T\mathbf{D_c^{-1}}\mathbf{B} = \mathbf{I}$.
## Décomposition de l'inertie des profils
Nous avons montré plus haut que pour représenter les profils lignes dans un repère dont les axes orthonormaux sont orientés selon les direction d'inertie maximale (nous appelons ces axes les \emph{facteurs}), il faut calculer le SVD :
Nous avons montré plus haut que pour représenter les profils lignes dans un repère dont les axes orthonormaux sont orientés selon les direction d'inertie maximale (nous appelons ces axes les \emph{axes factoriels}), il faut calculer le SVD :
$$\mathbf{P} - \mathbf{r}\mathbf{c}^T = \mathbf{A}\boldsymbol{\Delta}\mathbf{B}^T \; \text{avec } \; \mathbf{A}^T\mathbf{D_r^{-1}}\mathbf{A} = \mathbf{B}^T\mathbf{D_c^{-1}}\mathbf{B} = \mathbf{I}$$
Montrons que cette décomposition correspond bien à celle des profils lignes centrés (rappel : $\mathbf{c}$ est le centre de masse des profils lignes).
@ -251,11 +251,11 @@ Montrons que cette décomposition correspond bien à celle des profils lignes ce
& \mathbf{P} - \mathbf{r}\mathbf{c}^T = \left(\mathbf{D_r}\boldsymbol\alpha\right) \boldsymbol\Delta \boldsymbol\beta^T \; \text{avec} \; \left(\mathbf{D_r}\boldsymbol\alpha\right)^T\mathbf{D_r^{-1}}\left(\mathbf{D_r}\boldsymbol\alpha\right) = \boldsymbol\beta^T\mathbf{D_c^{-1}}\boldsymbol\beta = \mathbf{I} \\
\end{align*}
## Facteurs et coordonnées principales des profils
## Axes factoriels et facteurs
Ainsi, les colonnes de $\mathbf{B}$ sont les facteurs (ou \emph{axes principaux}), $\mathbf{D_c^{-1}}$-orthogonaux, pour le nuage centré des profils lignes. Nous montrerions de même que les colonnes de $\mathbf{A}$ sont les facteurs, $\mathbf{D_r^{-1}}$-orthogonaux, pour le nuage centré des profils colonnes.
Ainsi, les colonnes de $\mathbf{B}$ sont les \emph{axes factoriels}, $\mathbf{D_c^{-1}}$-orthogonaux, pour le nuage centré des profils lignes. Nous montrerions de même que les colonnes de $\mathbf{A}$ sont les axes factoriels, $\mathbf{D_r^{-1}}$-orthogonaux, pour le nuage centré des profils colonnes.
Notons $\mathbf{F}$ (respectivement, $\mathbf{G}$) une matrice dont les lignes sont les coordonnées principales (i.e., les coordonnées selon les facteurs) des profils lignes centrés (respectivement, des profils colonnes centrés).
Notons $\mathbf{F}$ (respectivement, $\mathbf{G}$) une matrice dont les lignes sont les facteurs (i.e., les coordonnées selon les axes factoriels) des profils lignes centrés (respectivement, des profils colonnes centrés).
\begin{align*}
& \mathbf{F} \\
= \{& \text{Par définition de $\mathbf{F}$} \} \\
@ -285,7 +285,7 @@ Notons $\mathbf{F}$ (respectivement, $\mathbf{G}$) une matrice dont les lignes s
= \{& \} \\
& \mathbf{P}^T = \mathbf{B}\boldsymbol\Delta\mathbf{A}^T + \mathbf{c}\mathbf{r}^T \\
\Rightarrow \{& \text{Multiplication à droite par $\mathbf{F}$ ; $\mathbf{F} = \mathbf{D_r^{-1}}\mathbf{A}\boldsymbol\Delta$ ; $\mathbf{A}^T\mathbf{D_r^{-1}}\mathbf{A} = \mathbf{I}$ ;} \\
\phantom{\Rightarrow}\phantom{\{}& \text{Les lignes de $\mathbf{F}$ sont les coordonnées principales des profils lignes centrés. Donc $\mathbf{r}^T \mathbf{F} = \mathbf{0}$.} \} \\
\phantom{\Rightarrow}\phantom{\{}& \text{Les lignes de $\mathbf{F}$ sont les facteurs des profils lignes \emph{centrés}. Donc $\mathbf{r}^T \mathbf{F} = \mathbf{0}$.} \} \\
& \mathbf{P}^T\mathbf{F} = \mathbf{B}\boldsymbol\Delta^2 \\
= \{& \} \\
& \mathbf{P}^T\mathbf{F}\boldsymbol\Delta^{-1} = \mathbf{B}\boldsymbol\Delta \\
@ -298,7 +298,7 @@ Nous montrerions de même :
& \mathbf{P}\mathbf{F}\boldsymbol\Delta^{-1} = \mathbf{A}\boldsymbol\Delta \\
\end{align*}
Nous découvrons ainsi des formules de transition entre les coordonnées principales des profils lignes et celles des profils colonnes.
Nous découvrons ainsi des formules de transition entre les facteurs des profils lignes et ceux des profils colonnes.
\begin{align*}
& \mathbf{G} \\
= \{& \text{Voir précédente dérivation} \} \\
@ -320,19 +320,19 @@ Nous découvrons ainsi des formules de transition entre les coordonnées princip
& \mathbf{R}\mathbf{G}\boldsymbol\Delta^{-1}
\end{align*}
## Inertie des profils selon les facteurs
## Inertie des profils selon les axes factoriels
Vérifions maintenant que la décomposition de l'inertie du nuage des profils lignes (ou du nuage des profils colonnes) selon les facteurs se représente par une matrice diagonale.
Vérifions maintenant que la décomposition de l'inertie du nuage des profils lignes (ou du nuage des profils colonnes) selon les axes factoriels se représente par une matrice diagonale.
$$\mathbf{F}^T\mathbf{D_r}\mathbf{F} = \left(\mathbf{D_r^{-1}} \mathbf{A}\boldsymbol\Delta\right)^T\mathbf{D_r}\mathbf{D_r^{-1}} \mathbf{A}\boldsymbol\Delta = \boldsymbol\Delta^2$$
Ainsi, l'inertie du nuage centré des profils lignes est : $\sum_i r_i \sum_k f_{ik}^2 = \sum_k \delta_k^2$ (où $\delta_k \triangleq \boldsymbol\Delta_{k,k}$).
Les lignes de $\boldsymbol\Phi \triangleq \mathbf{F}\boldsymbol\Delta^{-1}$ sont appelées les coordonnées standard des profils lignes. Selon ces dernières, les profils lignes centrés ont une inertie unité selon les facteurs ($\boldsymbol\Phi^T\mathbf{D_r}\boldsymbol\Phi = \mathbf{I}$).
Les lignes de $\boldsymbol\Phi \triangleq \mathbf{F}\boldsymbol\Delta^{-1}$ sont appelées les coordonnées standard des profils lignes. Selon ces dernières, les profils lignes centrés ont une inertie unité selon les axes factoriels ($\boldsymbol\Phi^T\mathbf{D_r}\boldsymbol\Phi = \mathbf{I}$).
$$\mathbf{G}^T\mathbf{D_c}\mathbf{G} = \left(\mathbf{D_c^{-1}} \mathbf{B}\boldsymbol\Delta\right)^T\mathbf{D_c}\mathbf{D_c^{-1}} \mathbf{B}\boldsymbol\Delta = \boldsymbol\Delta^2$$
Ainsi, l'inertie du nuage centré des profils colonnes est : $\sum_j c_j \sum_k g_{jk}^2 = \sum_k \delta_k^2$.
Les lignes de $\boldsymbol\Gamma \triangleq \mathbf{G}\boldsymbol\Delta^{-1}$ sont appelées les coordonnées standard des profils colonnes. Selon ces dernières, les profils colonnes centrés ont une inertie unité selon les facteurs ($\boldsymbol\Gamma^T\mathbf{D_c}\boldsymbol\Gamma = \mathbf{I}$).
Les lignes de $\boldsymbol\Gamma \triangleq \mathbf{G}\boldsymbol\Delta^{-1}$ sont appelées les coordonnées standard des profils colonnes. Selon ces dernières, les profils colonnes centrés ont une inertie unité selon les axes factoriels ($\boldsymbol\Gamma^T\mathbf{D_c}\boldsymbol\Gamma = \mathbf{I}$).
## Relations barycentriques
En observant l'égalité $\mathbf{G} = \mathbf{C}\mathbf{F}\boldsymbol\Delta^{-1}$, nous remarquons que la $j$-ème ligne de $\mathbf{G}$ (i.e., les coordonnées principales du $j$-ème profil colonne) est une combinaison linéaire (viz., un barycentre) des lignes de $\mathbf{F}$ (i.e., c'est à dire de l'ensemble des profils lignes exprimés dans leurs coordonnées principales) suivi par un changement d'échelle $1/\delta_k$ selon chacun des $K$ facteurs (d'où l'expression de quasi-barycentre). Nous ferions une observation symétrique à partir de l'égalité $\mathbf{F} = \mathbf{R}\mathbf{G}\boldsymbol\Delta^{-1}$.
En observant l'égalité $\mathbf{G} = \mathbf{C}\mathbf{F}\boldsymbol\Delta^{-1}$, nous remarquons que la $j$-ème ligne de $\mathbf{G}$ (i.e., les facteurs ou coordonnées principales du $j$-ème profil colonne) est une combinaison linéaire (viz., un barycentre) des lignes de $\mathbf{F}$ (i.e., c'est à dire de l'ensemble des profils lignes exprimés dans leurs coordonnées principales) suivi par un changement d'échelle $1/\delta_k$ selon chacun des $K$ facteurs (d'où l'expression parfois employée de \emph{quasi-barycentre}). Nous ferions une observation symétrique à partir de l'égalité $\mathbf{F} = \mathbf{R}\mathbf{G}\boldsymbol\Delta^{-1}$.
Comme $\mathbf{G} = \mathbf{C}\boldsymbol\Phi$, et que la somme des éléments d'une ligne de $\mathbf{C}$ (viz. un profil colonne), est égale à 1, un profil colonne exprimé dans ses coordonnées principales (viz. une ligne de $\mathbf{G}$) est un barycentre des profils lignes exprimés dans leurs coordonnées standard (viz. les lignes de $\boldsymbol\Phi$). Symétriquement, comme $\mathbf{F} = \mathbf{R}\boldsymbol\Gamma$, un profil ligne exprimé dans ses coordonnées principales est un barycentre des profils colonnes exprimés dans leurs coordonnées standard.
@ -348,199 +348,72 @@ Récapitulons les résultats des dérivations ci-dessus :
\item $\mathbf{F} = \mathbf{D_r^{-1/2}}\tilde{\mathbf{A}}\boldsymbol\Delta = \boldsymbol\Phi\boldsymbol\Delta$ ; $\mathbf{G} = \mathbf{D_c^{-1/2}} \tilde{\mathbf{B}}\boldsymbol\Delta = \boldsymbol\Gamma\boldsymbol\Delta$
\end{itemize}
Représentons sur une carte les coordonnées principales des dimensions n°1 et n°2 des profils lignes et des profils colonnes (voir Figure \@ref(fig:05-c-map-1-2)).
```{r 05-c-map-1-2, fig.width = 6, fig.cap = "Carte selon les facteurs 1 (x) et 2 (y)"}
scaleR <- 1/sqrt(r)
scaleC <- 1/sqrt(c)
S <- sweep(P - r %*% t(c), 1, scaleR, FUN = "*")
S <- sweep(S, 2, scaleC, FUN = "*")
dec <- svd(S)
Phi <- sweep(dec$u, 1, scaleR, FUN="*")
F <- sweep(Phi, 2, dec$d, FUN="*")
Gam <- sweep(dec$v, 1, scaleC, FUN="*")
G <- sweep(Gam, 2, dec$d, FUN="*")
# par(pty="s") # square plotting region, useful in interactive mode
plot(c(F[,1], G[,1]), c(F[,2], G[,2]), type = "n",
xlab="", ylab="", asp = 1, xaxt = "n", yaxt = "n")
text(c(F[,1], G[,1]), c(F[,2], G[,2]), c(rownames(P), colnames(P)),
adj = 0, cex = 0.6)
points(0, 0, pch = 3)
La fonction `ca` (\emph{analyse des correspondance}), voir code source en annexe de ce chapitre, effectue l'ensemble de ces calculs.
```{r}
cam <- ca(N) # cam pour 'correspondence analysis model'
```
Calculons la proportion d'inertie cumulée exprimée par les facteurs (voir Table \@ref(tab:05-c-vec-delta2)).
Calculons la proportion d'inertie cumulée exprimée par les axes factoriels (voir Table \@ref(tab:05-c-vec-delta2)).
```{r 05-c-vec-delta2}
tblTotalInr <- cumsum(round_preserve_sum(1000 * (dec$d^2 / sum(dec$d^2))))
tblTotalInr <- matrix(tblTotalInr, nrow=1)
colnames(tblTotalInr) <- 1:length(dec$d)
kbl(tblTotalInr, caption = "Proportion de l'inertie cumulée exprimée \
par les facteurs successifs (en millièmes)", booktabs = TRUE)
tblTotalInr <- tblTotalInr.ca(cam)
printTbl(tblTotalInr)
```
La carte de la Figure \@ref(fig:05-c-map-1-2) représente donc `r tblTotalInr[2]`‰ de l'inertie totale.
Représentons sur une carte les coordonnées principales des dimensions n°1 et n°2 des profils lignes et des profils colonnes (voir Figure \@ref(fig:05-c-map-1-2)). Nous avons calculé les parts de l'inertie du plan factoriel 1-2 expliquées par chaque profil ligne et chaque profil colonne (voir code source en annexe de ce chapitre). Nous utilisons ces valeurs pour fixer la taille de la police des labels des profils. Nous empruntons cette technique à [@greenacre2010biplots] (pp.181-182).
```{r 05-c-map-1-2, fig.width = 6, fig.cap = "Carte des axes factoriels 1 (x) et 2 (y)"}
plot(cam)
```
D'après la Table \@ref(tab:05-c-vec-delta2), la carte de la Figure \@ref(fig:05-c-map-1-2) représente `r tblTotalInr$tbl[2]`‰ de l'inertie totale.
## Importance des axes et des profils
Nous pouvons mesurer la contribution d'un profil ligne $i$ (resp. d'un profil colonne $j$) à l'inertie exprimée par l'axe principal $k$ :
Nous pouvons mesurer la contribution d'un profil ligne $i$ (resp. d'un profil colonne $j$) à l'inertie exprimée par l'axe factoriel $k$ :
\[CTR_{i,k} = \frac{r_i f_{i,k}^2}{\sum_i r_i f_{i,k}^2} = \frac{r_i f_{i,k}^2}{\delta_k^2}\]
\[CTR_{j,k} = \frac{c_j g_{j,k}^2}{\sum_j c_j g_{j,k}^2} = \frac{c_j g_{j,k}^2}{\delta_k^2}\]
Les contributions ont des valeurs entre $0$ et $1$. La somme des contributions des profils (ligne ou colonne) à un axe est égale à $1$.
Une heuristique utile : la contribution d'un profil ligne (resp. colonne) peut être considéré important quand sa valeur dépasse $1/I$ (resp. $1/J$).
```{r}
temp <- sweep(F^2, 1, r, FUN = "*")
CTRI <- sweep(temp, 2, dec$d^2, FUN = "/")
temp <- sweep(G^2, 1, c, FUN = "*")
CTRJ <- sweep(temp, 2, dec$d^2, FUN = "/")
```
Nous pouvons également mesurer la corrélation des axes principaux à un profil donné, c'est-à-dire la proportion de variance d'un profil qui peut être attribuée à un axe principal :
Nous pouvons également mesurer la corrélation des axes factoriels à un profil donné, c'est-à-dire la proportion de variance d'un profil qui peut être attribuée à un axe factoriel :
\[COR_{i,k} = \frac{f_{i,k}^2}{\sum_k f_{i,k}^2}\]
\[COR_{j,k} = \frac{g_{j,k}^2}{\sum_k g_{j,k}^2}\]
Les corrélations ont des valeurs entre $0$ et $1$. La somme des corrélations des axes à un profil (ligne ou colonne) est égale à $1$.
```{r}
temp <- temp <- F^2
CORI <- sweep(temp, 1, rowSums(temp), FUN = "/")
temp <- temp <- G^2
CORJ <- sweep(temp, 1, rowSums(temp), FUN = "/")
```
## Tableaux synthétiques des profils
Nous calculons la part d'inertie expliquée par chaque profil ligne :
\[ INR(i) = \sum_j r_i F_{i,j}^2 \]
Nous choisissons d'exprimer cette part d'inertie en proportion de l'inertie totale.
```{r}
INRI <- rowSums(sweep(F^2, 1, r, FUN="*"))
INRI <- INRI / sum(INRI)
INRJ <- rowSums(sweep(G^2, 1, c, FUN="*"))
INRJ <- INRJ / sum(INRJ)
```
Nous rappelons le poids d'un profil ligne : $PDS(i) = r_i$.
Nous synthétisons ces données dans un tableau pour les profils lignes sur 4 premiers facteurs (voir Table \@ref(tab:05-c-synth-row)).
```{r}
nbFact <- 4 # print data for the first nbFact factors
tempFactData <- c()
tempFactNames <- c()
for (k in 1:nbFact) {
tempFactData <- c(tempFactData,
round_preserve_sum(1000*F[,k]),
round_preserve_sum(1000*CORI[,k]),
round_preserve_sum(1000*CTRI[,k]))
tempFactNames <- c(tempFactNames, paste('F#', k, sep=""), 'COR', 'CTR')
}
tblI <- matrix( c(round_preserve_sum(1000*r),
round_preserve_sum(1000*INRI),
tempFactData),
ncol = 2 + nbFact * 3)
rownames(tblI) <- rownames(N)
colnames(tblI) <- c('PDS', 'INR', tempFactNames)
```
Nous synthétisons ces données dans un tableau pour les profils lignes sur les 4 premiers axes factoriels (voir Table \@ref(tab:05-c-synth-row)).
```{r 05-c-synth-row}
kbl(tblI, caption = "Synthèse des profils lignes \
sur les 4 premiers facteurs", booktabs = TRUE) %>%
add_header_above(c(" " = 3, "F#1" = 3, "F#2" = 3, "F#3" = 3, "F#4" = 3)) %>%
kable_styling(latex_options = c("striped", "scale_down"))
tblI <- tblRowProfiles.ca(cam, 4)
printTblWithHeaders(tblI)
```
Nous donnons aussi un tableau synthétique pour les profils colonnes sur les 4 premiers facteurs (voir Table \@ref(tab:05-c-synth-col)).
```{r}
nbFact <- 4 # print data for the first nbFact factors
tempFactData <- c()
tempFactNames <- c()
for (k in 1:nbFact) {
tempFactData <- c(tempFactData,
round_preserve_sum(1000*G[,k]),
round_preserve_sum(1000*CORJ[,k]),
round_preserve_sum(1000*CTRJ[,k]))
tempFactNames <- c(tempFactNames, paste('F#', k, sep=""), 'COR', 'CTR')
}
tblJ <- matrix( c(round_preserve_sum(1000*c),
round_preserve_sum(1000*INRJ),
tempFactData),
ncol = 2 + nbFact * 3)
rownames(tblJ) <- colnames(N)
colnames(tblJ) <- c('PDS', 'INR', tempFactNames)
```
Nous donnons aussi un tableau synthétique pour les profils colonnes sur les 4 premiers axes factoriels (voir Table \@ref(tab:05-c-synth-col)).
```{r 05-c-synth-col}
kbl(tblJ, caption = "Synthèse des profils colonnes \
sur les 4 premiers facteurs", booktabs = TRUE) %>%
add_header_above(c(" " = 3, "F#1" = 3, "F#2" = 3, "F#3" = 3, "F#4" = 3)) %>%
kable_styling(latex_options = c("striped", "scale_down"))
tblJ <- tblColProfiles.ca(cam, 4)
printTblWithHeaders(tblJ)
```
## Éléments d'interprétation des cartes des profils améliorées
Nous calculons les parts de l'inertie du plan factoriel 1-2 expliquées par chaque profil ligne et chaque profil colonne. Nous utilisons ces valeurs pour fixer la taille de la police des labels des profils. Nous empruntons cette technique à [@greenacre2010biplots] (pp.181-182).
```{r}
# Parts de l'inertie du plan factoriel 1-2 expliquée par chaque profil ligne
rowcon <- r * (F[,1]^2 + F[,2]^2) / sum(dec$d[1:2]^2)
rnames <- rownames(P)
rnames[rowcon < 0.01] <- "."
rsize <- log(1 + exp(1) * rowcon^0.3)
rsize[rowcon < 0.01] <- 1
# Parts de l'inertie du plan factoriel 1-2 expliquée par chaque profil colonne
colcon <- c * (G[,1]^2 + G[,2]^2) / sum(dec$d[1:2]^2)
cnames <- colnames(P)
cnames[colcon < 0.01] <- "."
csize <- log(1 + exp(1) * colcon^0.3)
csize[colcon < 0.01] <- 1
```
Nous générons la Figure \@ref(fig:05-c-map-1-2-better) qui est une carte améliorée du plan factoriel 1-2.
```{r 05-c-map-1-2-better, fig.width = 6, fig.cap = "Carte améliorée selon les facteurs 1 (x) et 2 (y)"}
plot(c(F[,1], G[,1]), c(F[,2], G[,2]), type = "n",
xlab="", ylab="", asp = 1, xaxt = "n", yaxt = "n")
text(c(F[,1], G[,1]), c(F[,2], G[,2]), c(rnames, cnames),
adj = 0, cex = c(rsize, csize))
points(0, 0, pch = 3)
```
* Les plus gros pays exportateurs (resp. Maroc, URSS et USA) forment trois pôles d'attraction distincts.
* $eCC$ (URSS en tant qu'exportateur) et $iEE$ (autres pays d'Europe de l'Est en tant qu'importateurs) sont isolés à une extrémité du premier facteur et contribuent de façon majeure à l'inertie de ce dernier.
* $eCC$ (URSS en tant qu'exportateur) et $iEE$ (autres pays d'Europe de l'Est en tant qu'importateurs) sont isolés à une extrémité du premier axe factoriel et contribuent de façon majeure à l'inertie de ce dernier.
* Les USA et l'URSS s'opposent sur le premier axe factoriel.
* Les USA et le Maroc s'opposent sur le second axe factoriel.
* Sur le plan factoriel 1-2, autour des USA gravitent les importateurs : d'abord le Canada (presque entièrement expliqué par le plan factoriel 1-2), mais aussi le Japon, le Brésil et l'Inde (le Japon et le Brésil sont très bien expliqués par le plan factoriel 1-2, trandis que l'inertie de l'Inde est fortement expliquée par l'axe factoriel n°4 sur lequel elle s'oppose en particulier à l'Espagne, voir Table \@ref(tab:05-c-synth-row)).
On poursuivrait l'analyse en explorant par exemple la carte du plan factoriel 2-3 (voir Figure \@ref(fig:05-c-map-2-3)).
```{r}
# Parts de l'inertie du plan factoriel 2-3 expliquée par chaque profil ligne
rowcon <- r * (F[,2]^2 + F[,3]^2) / sum(dec$d[2:3]^2)
rnames <- rownames(P)
rnames[rowcon < 0.01] <- "."
rsize <- log(1 + exp(1) * rowcon^0.3)
rsize[rowcon < 0.01] <- 1
# Parts de l'inertie du plan factoriel 2-3 expliquée par chaque profil colonne
colcon <- c * (G[,2]^2 + G[,3]^2) / sum(dec$d[2:3]^2)
cnames <- colnames(P)
cnames[colcon < 0.01] <- "."
csize <- log(1 + exp(1) * colcon^0.3)
csize[colcon < 0.01] <- 1
```
```{r 05-c-map-2-3, fig.width = 6, fig.cap = "Carte des facteurs 2 (x) et 3 (y)"}
plot(c(F[,2], G[,2]), c(F[,3], G[,3]), type = "n",
xlab="", ylab="", asp = 1, xaxt = "n", yaxt = "n")
text(c(F[,2], G[,2]), c(F[,3], G[,3]), c(rnames, cnames),
adj = 0, cex = c(rsize, csize))
points(0, 0, pch = 3)
```{r 05-c-map-2-3, fig.width = 6, fig.cap = "Carte des axes factoriels 2 (x) et 3 (y)"}
plot(cam, d1=2, d2=3)
```
## Annexe code source

View File

@ -14,3 +14,162 @@ round_preserve_sum <- function(x, digits = 0) {
y[indices] <- y[indices] + 1
y / up
}
ca <-
function(N)
{
n <- sum(sum(N))
P <- (1/n) * N
r <- rowSums(P)
c <- colSums(P)
R <- sweep(P, MARGIN = 1, 1/r, '*')
C <- sweep(t(P), MARGIN = 1, 1/c, '*')
scaleR <- 1/sqrt(r)
scaleC <- 1/sqrt(c)
S <- sweep(P - r %*% t(c), 1, scaleR, FUN = "*")
S <- sweep(S, 2, scaleC, FUN = "*")
dec <- svd(S)
Phi <- sweep(dec$u, 1, scaleR, FUN="*")
F <- sweep(Phi, 2, dec$d, FUN="*")
Gam <- sweep(dec$v, 1, scaleC, FUN="*")
G <- sweep(Gam, 2, dec$d, FUN="*")
temp <- sweep(F^2, 1, r, FUN = "*")
CTRI <- sweep(temp, 2, dec$d^2, FUN = "/")
temp <- sweep(G^2, 1, c, FUN = "*")
CTRJ <- sweep(temp, 2, dec$d^2, FUN = "/")
temp <- temp <- F^2
CORI <- sweep(temp, 1, rowSums(temp), FUN = "/")
temp <- temp <- G^2
CORJ <- sweep(temp, 1, rowSums(temp), FUN = "/")
INRI <- rowSums(sweep(F^2, 1, r, FUN="*"))
INRI <- INRI / sum(INRI)
INRJ <- rowSums(sweep(G^2, 1, c, FUN="*"))
INRJ <- INRJ / sum(INRJ)
r <- list(P=P, r=r, c=c, F=F, G=G, Phi=Phi, Gam=Gam, sv=dec$d,
CTRI=CTRI, CTRJ=CTRJ, CORI=CORI, CORJ=CORJ,
INRI=INRI, INRJ=INRJ)
class(r) <- "ca"
return(r)
}
plot.ca <-
function(o, d1 = NULL, d2 = NULL)
{
if(class(o) != "ca") {
warning("Object is not of class 'ca'")
UseMethod("plot")
return(invisible(NULL))
}
if(is.null(d1)) d1<-1
if(is.null(d2)) d2<-2
# Part de l'inertie du plan factoriel d1-d2 expliquée par chaque profil ligne
rowcon <- o$r * (o$F[,d1]^2 + o$F[,d2]^2) / (o$sv[d1]^2 + o$sv[d2]^2)
rnames <- rownames(o$P)
rnames[rowcon < 0.01] <- "."
rsize <- log(1 + exp(1) * rowcon^0.3)
rsize[rowcon < 0.01] <- 1
# Part de l'inertie du plan factoriel d1-d2 expliquée par chaque profil colonne
colcon <- o$c * (o$G[,d1]^2 + o$G[,d2]^2) / (o$sv[d1]^2 + o$sv[d2]^2)
cnames <- colnames(o$P)
cnames[colcon < 0.01] <- "."
csize <- log(1 + exp(1) * colcon^0.3)
csize[colcon < 0.01] <- 1
plot(c(o$F[,d1], o$G[,d1]), c(o$F[,d2], o$G[,d2]), type = "n",
xlab="", ylab="", asp = 1, xaxt = "n", yaxt = "n")
text(c(o$F[,d1], o$G[,d1]), c(o$F[,d2], o$G[,d2]), c(rnames, cnames),
adj = 0, cex = c(rsize, csize))
points(0, 0, pch = 3)
}
tblTotalInr.ca <-
function(o)
{
if(class(o) != "ca") {
warning("Object is not of class 'ca'")
return(invisible(NULL))
}
tblTotalInr <- cumsum(round_preserve_sum(1000 * (o$sv^2 / sum(o$sv^2))))
tblTotalInr <- matrix(tblTotalInr, nrow=1)
colnames(tblTotalInr) <- 1:length(o$sv)
r <- list(tbl=tblTotalInr, title="Proportion de l'inertie cumulée exprimée \
par les axes factoriels successifs (en millièmes)")
return(r)
}
printTbl <-
function(o)
{
kbl(o$tbl, caption = o$title, booktabs = TRUE)
}
printTblWithHeaders <-
function(o)
{
kbl(o$tbl, caption = o$title, booktabs = TRUE) %>%
add_header_above(o$headers) %>%
kable_styling(latex_options = c("striped", "scale_down"))
}
tblRowProfiles.ca <-
function(o, nbFact)
{
if(class(o) != "ca") {
warning("Object is not of class 'ca'")
return(invisible(NULL))
}
tempFactData <- c()
tempFactNames <- c()
subHeaders <- c(" " = 3) # first empty subheader for row names, PDS and INR
for (k in 1:nbFact) {
tempFactData <- c(tempFactData,
round_preserve_sum(1000*o$F[,k]),
round_preserve_sum(1000*o$CORI[,k]),
round_preserve_sum(1000*o$CTRI[,k]))
factIdStr <- paste('F#', k, sep="")
tempFactNames <- c(tempFactNames, factIdStr, 'COR', 'CTR')
subHeaders[[factIdStr]] <- 3
}
tblI <- matrix( c(round_preserve_sum(1000*o$r),
round_preserve_sum(1000*o$INRI),
tempFactData),
ncol = 2 + nbFact * 3)
rownames(tblI) <- rownames(o$P)
colnames(tblI) <- c('PDS', 'INR', tempFactNames)
r <- list(tbl=tblI, title="Synthèse des profils lignes", headers=subHeaders)
}
tblColProfiles.ca <-
function(o, nbFact)
{
if(class(o) != "ca") {
warning("Object is not of class 'ca'")
return(invisible(NULL))
}
tempFactData <- c()
tempFactNames <- c()
subHeaders <- c(" " = 3) # first empty subheader for row names, PDS and INR
for (k in 1:nbFact) {
tempFactData <- c(tempFactData,
round_preserve_sum(1000*o$G[,k]),
round_preserve_sum(1000*o$CORJ[,k]),
round_preserve_sum(1000*o$CTRJ[,k]))
factIdStr <- paste('F#', k, sep="")
tempFactNames <- c(tempFactNames, factIdStr, 'COR', 'CTR')
subHeaders[[factIdStr]] <- 3
}
tblJ <- matrix( c(round_preserve_sum(1000*o$c),
round_preserve_sum(1000*o$INRJ),
tempFactData),
ncol = 2 + nbFact * 3)
rownames(tblJ) <- colnames(o$P)
colnames(tblJ) <- c('PDS', 'INR', tempFactNames)
r <- list(tbl=tblJ, title="Synthèse des profils colonnes", headers=subHeaders)
}