Tag Archives: ACP

Analyse des Correspondances et ACP

Considérons les données suivantes, correspondant à des professions au sein d’un couple. On dispose du tableau de contingence suivant

> base=read.table(
"http://freakonometrics.free.fr/epoux.csv",
sep=";",header=TRUE)
> rownames(base)=base$Nom

Classiquement, avec ce genre de données, on est habitué à utiliser la distance du chi-deux, et les contributions au chi-deux pour voir les modalités qui correspondent fortement

> M=base[1:9,2:10]
> CT=chisq.test(M)

On peut le visualiser de la manière suivante,

> mosaicplot(t(M),main="",las=2,shade=TRUE)

avec les maris en ligne, et les épouses en colonnes. Les contribution importantes sont en bleu ou rouge foncé, les couleurs correspondant respectivement à un lien ‘positif’ (probabilité jointe plus importante que le cas indépendant) ou un lien ‘négatif’ (probabilité jointe moins importante que le cas indépendant)

ou dans l’autre sens

> mosaicplot(M,main="",las=2,shade=TRUE)

mais avec la même conclusion que précédemment: il y a de fortes valeurs bleues sur la diagonale.

Autrement dit, les couples sont relativement homogènes en terme de professions.

En analyse des correspondance, on va regarder notre tableau de contingence, en ligne ou en colonne. Par exemple, on peut définir les profils lignes, qui sont des vecteurs de probabilité

> N=M 
> L=N/apply(N,1,sum)

En notant , on peut écrire

Le barycentre de nos vecteurs de lignes est ici

Là encore, en notant , on peut écrire ce barycentre de manière matricielle, .

> Lbar=apply(N,2,sum)/sum(N)
> matL0=t(t(L)-Lbar)

A chaque point, on a associé la fréquence (relative) comme poids, , ce qui correspond à utiliser la matrice . Pour mesurer la distance entre deux points, on va pondérer la distance euclidienne par l’inverse des des probabilités, , ce qui revient à utiliser la matrice . La distance entre deux lignes est

On va alors faire une analyse en composantes principales avec ces différents poids. D’un point de vue matriciel, on va alors l’analyse spectrale de

En particulier, on note  les vecteurs propres, et on défini les composantes principales

La projection sur les deux premières composantes des lignes donne ici

> library(FactoMineR)
> acpL=PCA(matL0,scale.unit=FALSE,
+          row.w=(apply(N,1,sum)),
+          col.w=1/(apply(N,2,sum)))

L’idée est de visualiser le nuage des individus correspondants aux lignes. Dans un second temps, on va faire exactement pareil, mais sur les colonnes

> C=t(t(N)/apply(N,2,sum))

Le barycentre est ici

> Cbar=apply(N,1,sum)/sum(N)
> matC0=C-Cbar

et on peut alors faire une analyse en composante principale

> acpC=PCA(t(matC0),scale.unit=FALSE,
+          row.w=(apply(N,2,sum)),
+          col.w=1/(apply(N,1,sum)))

ce qui donne, en regardant le nuage des individus,

La magie de l’analyse des correspondances est qu’on ‘peut’ représenter les deux projections des individus sur le même plan,

> matC=acpC$ind$coord
> matL=acpL$ind$coord
> plot(matC[,1:2],col="red",xlim=c(-.002,.015))
> text(matC[,1:2],rownames(matC),col="red",pos=4)
> points(matL[,1:2],col="blue")
> text(matL[,1:2],rownames(matL),col="blue")

C’est exactement ce que fait la fonction

> afc=CA(N)

Dualité Individus / Variables en ACP

Dans le cours d’analyse des données, on va commencer l’analyse en composantes principales. Comme on l’a déjà évoqué à plusieurs reprises, il y a une dualité entre les individus dans l’espace des variables, et les variables dans l’espace des individus. Histoire de faire des dessins, on va se contenter de 3 individus, avec 3 variables,

> BD=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/securite.txt",
header=TRUE)
>  M=as.matrix(BD[,2:10])
>  colnames(M)=names(BD)[2:10]
>  rownames(M)=BD$dep
>  M=M[1:3,1:3]
>  n=nrow(M)  
>  M_s=scale(M,scale=TRUE)

 

On peut commencer par visualiser les individus

Les directions principales sont ici

> U=eigen(t(M_s)%*%(M_s))$vectors

 

Les deux premières nous donne le plan principal

La projection sur le plan donne ici

Mais on peut aussi visualiser les données dans l’autre dimension.

Rappelons que les variables ont été centrées et réduites. La encore, l’analyse spectrale nous donne les vecteurs principaux (qui vont maximiser l’inertie en projection),

>  V=eigen((M_s)%*%t(M_s))$vectors

 

Là encore, on peut s’intéresser au plan principal

Là encore on peut regarder la projection dans cette dimension,

Comme on a centré et réduit nos variables, on voit apparaître le “cercle de corrélations“.

Le principe de dualité dit que l’on peut superposer les projections

 

ce qui permettra de faire une réelle analyse des données

Manipulation matricielle de données

Pour travailler un peu la manipulation de matrices, dans le cadre du cours d’analyse de données, quelques lignes de code.

La matrice de données est ici

> BD=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/securite.txt",header=TRUE)
> M=as.matrix(BD[,2:10])
> colnames(M)=names(BD)[2:10]
> rownames(M)=BD$dep
> head(M[,1:5])
    infra  vols  eco crim  vma
D1  44.11 29.51 4.14 2.79 0.08
D2  45.97 26.29 3.97 4.54 0.07
D3  38.83 21.90 4.61 3.33 0.03
D4  49.68 30.11 5.28 4.54 0.07
D5  47.67 27.19 4.94 3.91 0.05
D6 109.21 77.71 6.87 5.25 0.24

La décomposition spectrale de

> round(t(M)%*%M)
       infra   vols   eco  crim vma  vvi  camb  
infra 308599 191433 29179 21290 642 7042 30951 
vols  191433 119961 17859 12967 404 4427 19335 
eco    29179  17859  2946  2042  59  648  2942  
crim   21290  12967  2042  1564  42  448  2129  
vma      642    404    59    42   2   19    62   
vvi     7042   4427   648   448  19  237   665  
camb   30951  19335  2942  2129  62  665  3200  
roul   55363  34745  5089  3792 113 1217  5612 
auto   24841  15747  2214  1676  53  545  2520

(j’ai viré les deux dernières colonnes pour des raisons de présentation) donne ici

> eigen(t(M)%*%M)
$values
[1] 4.471766e+05 1.239465e+03 3.462936e+02 1.764038e+02 1.303316e+02 6.623592e+01 2.624739e+01 1.761801e+01 1.564292e-01
 
$vectors

(etc). Pour rappel, les liens entre la diagonalisation et l’analyse spectrale se voit avec l’écriture

> P=eigen(t(M)%*%M)$vectors
> P%*%diag(eigen(t(M)%*%M)$values)%*%t(P)

Maintenant, deux remarques. Tout d’abord le lien entre la décomposition spectrale de cette matrice, et la décomposition en valeurs singulières

> sqrt(eigen(t(M)%*%M)$values)
[1] 668.7126681  35.2060415  18.6089656  13.2817097  11.4162872   8.1385454   5.1232210   4.1973809   0.3955113
> svd(M)
$d
[1] 668.7126681  35.2060415  18.6089656  13.2817097  11.4162872   8.1385454   5.1232210   4.1973809   0.3955113
 
$u
             [,1]          [,2]          [,3]          [,4]         [,5]          [,6]         [,7]         [,8]          [,9]
 [1,] -0.08156446 -0.0632514738 -2.596265e-02  0.0960587227 -0.064759395 -0.0430947134 -0.024719699 -0.074347281  0.0550356487
 [2,] -0.08087892  0.0541859480 -5.777492e-02 -0.0628249751 -0.099084701 -0.1311683858  0.083484434  0.011872601 -0.0195034904

(etc) et la décomspotion spectrale de l’autre produit matriciel

> sqrt(eigen(M%*%t(M))$values)
 [1] 6.687127e+02 3.520604e+01 1.860897e+01 1.328171e+01 1.141629e+01 8.138545e+00 5.123221e+00 4.197381e+00 3.955113e-01 7.546505e-06 6.656212e-06 5.285131e-06 4.270035e-06

(etc).

Maintenant, pour comprendre un peu mieux la recherche de composantes pricinpales, considérons les deux variables

> sM=M[,c(1,3)]
> plot(sM)

 

On a tout intérêt à centrer et réduire les variables (ou changer la métrique, comme vu en cours)

> sMcr=sM
> for(j in 1:2) sMcr[,j]=(sMcr[,j]-mean(sMcr[,j]))/sd(sMcr[,j])
> plot(sMcr)

Avant de faire des projections sur un axe, introduisons deux fonctions

> prod_scal=function(x,u){
+   ps=rep(NA,nrow(x))
+   for(i in 1:nrow(x)) ps[i]=sum(x[i,]*u)
+   return(ps)
+ }

> proj=function(x,u){
+   px=x
+   for(j in 1:length(u)){
+     px[,j]=prod_scal(x,u)/sqrt(sum(u^2))*u[j]  
+   }
+   return(px)
+ }

Par exemple, si on projette sur l’axe des abscisses,

> points(proj(sMcr,c(1,0)),col="blue")

On peut ensuite chercher la direction de l’axe qui nous donnera le nuage de points qui aura la plus grande inertie

> inertie=function(x) sum(x^2)
> Theta=seq(0,3.141592657,length=101)
> V=unlist(lapply(Theta,function(theta) inertie(proj(sMcr,c(cos(theta),sin(theta))))))
> plot(Theta,V,type='l')

> (angle=optim(0,function(theta) -inertie(proj(sMcr,c(cos(theta),
sin(theta)))))$par)
[1] 0.7853516

Visuellement, on obtient

> plot(sMcr)
> psMcr=proj(sMcr,c(cos(angle),sin(angle)))
> points(psMcr,col="blue")

https://freakonometrics.hypotheses.org/files/2016/01/proj1.png

On notera que l’axe qui donne l’inertie maximale est lié aux vecteurs propres de la décomposition spectrale (celui associé à la plus grande valeur propre)

> c(cos(angle),sin(angle))
[1] 0.7071397 0.7070738
> eigen(t(sMcr)%*%sMcr)
$values
[1] 159.64663  26.35337
 
$vectors
          [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068

On continuera demain, à manipuler des matrices de données, de faire des projections, avant de commencer l’analyse en composantes principales, à proprement parler, la semaine prochaine

Analyse des Correspondances, suite

Hier soir, on voyait comment faire une analyse des correspondances à partir du tableau de contingence. Mais on peut avoir le problème autrement, à partir des individus. Certes, ces derniers ne sont pas observés (vraiment), mais peu importe.

> data(HairEyeColor)
> N = HairEyeColor[,,"Male"] + HairEyeColor[,,"Female"]
> Hair = rep(rep(rownames(N),ncol(N)),
+ as.vector(N))
> Eye =  rep(colnames(N),apply(N,2,sum))
> df = data.frame(Ind=1:sum(N),Hair,Eye)
> tail(df)
    Ind  Hair   Eye
587 587 Blond Green
588 588 Blond Green
589 589 Blond Green
590 590 Blond Green
591 591 Blond Green
592 592 Blond Green

Continue reading Analyse des Correspondances, suite

Le “cercle des corrélations” en ACP

Tout au long de l’analyse en composantes principales, nous avons parlé du “cercle des corrélations“. Cette notion est fondamentale en ACP, et je voulais insister un peu sur ce point pour le démystifier,

A chaque point-variable, on associe un point dont la coordonnée sur un axe factoriel est une mesure de la corrélation entre cette variable et le facteur. Dans l’espace de dimension p la distance des points-variables à l’origine est égale à 1. Donc par projection sur un plan factoriel les points-variables s’inscrivent dans un cercle de rayon 1 – le cercle des corrélations – et sont d’autant plus proche du bord du cercle que le pointvariable est bien représenté par le plan factoriel, c’est-à-dire que la variable est bien corrélée avec les deux facteurs constituant ce plan.

L’angle entre 2 variables, mesuré par son cosinus, est égal au coefficient de corrélation linéaire entre les 2 variables: cos(angle) = r(X1,X2)

  • si les points sont très proches (angle proche de 0) : cos(angle) = r(X1,X2) = 1 donc X1 et X2 sont très fortement corrélés positivement
  • si a est égal à 90° , cos(angle) = r(X1,X2) = 0 alors pas de corrélation linéaire entre X1 et X2
  • si les points sont opposés, a vaut 180°, cos(angle) = r(X1,X2) = -1 : X1 et X2 sont très fortement corrélés négativement

Attention, on ne peut interpréter que les variables situées au bord du cercle. En effet, seuls ces points là sont “bien” projetés, comme le rappelle le dession ci-dessous,

ACP avec XLstat, un peu de terminologie

Pour compléter (et clarifier) le premier TD, je voulais revenir sur le type d’ACP proposé sous XLstat,

Pour l’analyse de la covariance, on travaille sur les données brutes. L’analyse de Pearson propose de travailler sur les données centrées-réduites. L’analyse de Spearman est basée sur les rangs. Enfin, l’analyse polychorique propose de travailler sur les rangs normalisés. Les quatre transformations sont proposées ci-dessous,

Les rangs normalisés correspondent à l’étude d’un facteur Gaussien sous-jacent: on retrouve le quantile normal (centré-réduit) pour des seuils de 20%, 40%, 60% et 80%, c’est à dire le rang divisés par la taille de l’échantillon+1. La “robustesse” de ces deux dernières méthodes a été vu en particulier en introduisant un point aberrant dans l’échantillon.

Analyse en Composantes Principales

Les slides de la première partie sont en ligne sur mon blog. Parmi les cours en ligne, à l’agro (Grignon), de Stéphane Guidon, ou plus formalisé, des notes de cours de Joseph Gergaud.

Sinon, pour accéder aux bases de données, les codes R sont les suivants,
elections<-read.table("https://perso.univ-rennes1.fr/arthur.charpentier/election95.csv",header=TRUE,sep=";")
insecurite<-read.table("https://perso.univ-rennes1.fr/arthur.charpentier/ADD-ex-villes.txt",header=TRUE)

La base de données qui seront utilisées en TD est la base departement.xls. On pourra aussi reprendre quelques exemple mentionnés en cours, en particulier l’étude de l’express sur l’insécurité, villes.xls, ou encore la base de décathlon mentionnée en cours, decathlon.xls.

Ce premier cours sera relativement technique puisque seront présentées les principales notions sur les distances, les projections, et utilisera plusieurs résultats “classiques” d’algèbre linéaire et matricielle. En particulier, sur la décomposition en valeurs singulières, je renvoie les élèves intéressés à l’article de wikipedia (svd). Je mets aussi un lien vers des slides sur la décomposition factoriel, en liaison avec la recherche de valeurs propres, et de la décomposition en valeurs singulières.

Sinon pour quelques groupes, les TD se feront sous XLStat, qui comme on le verra, permet de générer un certain nombre de graphiques (que l’on s’attachera à commenter)