Tag Archives: AFC

ACT6100, analyse non-supervisée

On avance tranquillement dans le cours ACT6100 d’analyse des données en actuariat. Les supports de cours sont en ligne sur https://github.com/freakonometrics/ACT6100. Cette session (mais aussi celle d’hiver) étant en distanciel, le cours est asynchrone, et je poste régulièrement des capsules vidéos. Les capsules présentant les principales méthodes d’analyse non-supervisée sont en ligne

  1. ACP (1) (PCAvideo pdf (48:42)
  2. ACP (2) video pdf (31:10)
  3. ACP (3) video pdf (47:51)
  4. ACP (4) video pdf (39:59)
  5. ACP (5) video pdf (28:53)
  6. CA (1) video pdf (38:33)
  7. CA (2) video pdf (46:51)
  8. MCA (1) video pdf (28:03)
  9. Clusters (k-means) video pdf (48:20)
  10. Clusters (CAH) video pdf (37:38)
  11. NA & Imputations (k-nn) video pdf (17:47)
  12. NA & Imputations (ACP) video pdf (15:28)

Si les liens des vidéos ne marchent pas, je renvoie vers l’ensemble des capsules du cours, ici.

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é \boldsymbol{L}_i=\left(\frac{p_{i1}}{p_{i\cdot}},\cdots,\frac{p_{iJ}}{p_{i\cdot}}\right)\in\mathbb{R}^J

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

En notant \boldsymbol{D}_I=\text{diag}(p_{1\cdot},\cdots,p_{I\cdot}), on peut écrire \boldsymbol{L}=\boldsymbol{D}_I^{-1}\boldsymbol{N}Le barycentre de nos vecteurs de lignes est ici\overline{\boldsymbol{L}}=\left(p_{\cdot 1},\cdots,p_{\cdot J}\right)\in\mathbb{R}^J

Là encore, en notant \boldsymbol{D}_J=\text{diag}(p_{\cdot 1},\cdots,p_{\cdot J}), on peut écrire ce barycentre de manière matricielle, \overline{\boldsymbol{L}}=\mathbb{I}_{IJ}\boldsymbol{D}_J.\boldsymbol{L}_{0,i}=\boldsymbol{L}_i-\overline{\boldsymbol{L}}

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

A chaque point, on a associé la fréquence (relative) comme poids, p_{i\cdot}, ce qui correspond à utiliser la matrice \boldsymbol{D}_I. Pour mesurer la distance entre deux points, on va pondérer la distance euclidienne par l’inverse des des probabilités, p_{\cdot j}, ce qui revient à utiliser la matrice \boldsymbol{D}_J^{-1}. La distance entre deux lignes estd(\boldsymbol{L}_{i_1},\boldsymbol{L}_{i_2})=\sum_{j=1}^J\frac{1}{p_{\cdot j}}\left(\frac{p_{i_1j}}{p_{i_1\cdot}}-\frac{p_{i_2j}}{p_{i_2\cdot}}\right)^2

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\boldsymbol{H}=\boldsymbol{D}_J^{-1/2}\boldsymbol{L}_0^{\text{T}}\boldsymbol{D}_I\boldsymbol{L}_0\boldsymbol{D}_J^{-1/2}

En particulier, on note \boldsymbol{U}_1,\cdots,\boldsymbol{U}_k les vecteurs propres, et on défini les composantes principales\boldsymbol{C}=\boldsymbol{L}_0\boldsymbol{D}_J^{-1/2}\boldsymbol{U}

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\boldsymbol{C}_j=\left(\frac{p_{1j}}{p_{\cdot j}},\cdots,\frac{p_{jI}}{p_{\cdot j}}\right)\in\mathbb{R}^I

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

Le barycentre est ici\overline{\boldsymbol{C}}=\left(p_{1\cdot},\cdots,p_{I\cdot}\right)\in\mathbb{R}^I

> 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)

Note: en réalité, si on regarde les coordonnées des lignes par ACP on obtient

acpL$ind$coord
                    Dim.1     Dim.2     Dim.3     Dim.4     Dim.5
Agriculteur (M) -0.000056 -0.000414  0.002940  0.000150  0.000011
Artisan (M)     -0.000014  0.000152 -0.000057  0.000323  0.001279
Cadres (M)      -0.000262  0.001734 -0.000015  0.000359 -0.000219
Prof_Int (M)    -0.000194  0.000450 -0.000003 -0.000504  0.000041
Employe (M)     -0.000104 -0.000317 -0.000083 -0.000188  0.000011
Ouvrier (M)     -0.000052 -0.000764 -0.000080  0.000133 -0.000141
Retraite (M)     0.007771  0.000502 -0.000001 -0.000089 -0.000065
Inactif (M)     -0.000018 -0.000150 -0.000163  0.000751 -0.000007
Inconnu(M)       0.000261 -0.000492 -0.000554  0.000398  0.000648

alors que

CA(M, graph = FALSE)$row$coord
                 Dim 1  Dim 2  Dim 3  Dim 4  Dim 5
Agriculteur (M) -0.028 -0.209  1.481  0.076  0.005
Artisan (M)     -0.007  0.077 -0.029  0.163  0.644
Cadres (M)      -0.132  0.873 -0.007  0.181 -0.110
Prof_Int (M)    -0.098  0.227 -0.001 -0.254  0.020
Employe (M)     -0.053 -0.160 -0.042 -0.095  0.006
Ouvrier (M)     -0.026 -0.385 -0.040  0.067 -0.071
Retraite (M)     3.914  0.253 -0.001 -0.045 -0.033
Inactif (M)     -0.009 -0.075 -0.082  0.378 -0.004
Inconnu(M)       0.132 -0.248 -0.279  0.201  0.327

ce qui donne un ratio constant,

acpL$ind$coord/CA(M, graph = FALSE)$row$coord
                      Dim.1       Dim.2       Dim.3       Dim.4       Dim.5
Agriculteur (M) 0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Artisan (M)     0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Cadres (M)      0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Prof_Int (M)    0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Employe (M)     0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Ouvrier (M)     0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Retraite (M)    0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Inactif (M)     0.001985182 0.001985182 0.001985182 0.001985182 0.001985182
Inconnu(M)      0.001985182 0.001985182 0.001985182 0.001985182 0.001985182

qui correspond à 1/\sqrt{n}

1/sqrt(sum(M))
[1] 0.001985182

Comprendre les tableaux croisés, 1

Je voulais faire quelques billets rapides sur les tableaux croisés, suite à des demandes consécutives à mon billet sur le sexe des bébés (ici). Pour commencer, je vais m’intéresser au cas où les variables supposées indépendantes (d’autres billets seront l’occasion de tester d’autres types d’hypothèses). Histoire de trouver des données un peu originales, on peut traîner sur internet. Etant tombé (presque) par hasard sur le fichier pdf ici, j’ai voulu m’intéresser à l’origine des enseignants-chercheurs recrutés lors de la campagne 2008 (comme le dit le titre du document en ligne sur le site du ministère). Prenons deux variables catégorielles géographiques, par exemple

  • https://perso.univ-rennes1.fr/arthur.charpentier/latex/cont-01.png, ville (au sens large) où l’on a soutenu sa thèse
  • https://perso.univ-rennes1.fr/arthur.charpentier/latex/cont-02.png, ville (au sens large) où l’on a obtenu un poste de maître de conférence

Compte tenu des règles de non-recrutement local, on devrait pouvoir supposer que ces variables sont indépendantes. En fait, je ne raisonne pas ici par université, mais par ville, ce qui correspond à un recrutement local au sens large (passer de Paris 6 à Paris 7 sera ici considéré comme “local“).
L’étude a été faite à partir des données disponibles ici et . Petite précaution d’usage: ce sont les données de la campagne 2008, avant que les comités de sélection ne soient mis en place. Il conviendrait de refaire l’étude une fois les nouvelles données publiques. On a deux bases, une en droit-économie-gestion, et une en sciences. J’insiste sur le fait que l’on a des données sur une année seulement, donc avec peu d’observations… Mais comme l’indique le titre, je vais parler des outils visualisation des tableaux croisés, pas faire un billet sur le recrutement local. Voilà pour les précautions d’usage…
Pour visualiser ici, j’ai fait des regroupements, car avec 25 villes, ça ne tient pas dans mon billet (mais par la suite, l’analyse sera plus fine). Les données de base sont les suivantes, les lignes étant l’endroit où un doctorant a soutenu sa thèse (en économie ou droit), et en colonnes, l’endroit où a été recruté un maître de conférence,

Rennes Province
(autre)
Région
parisienne
total
Rennes 5 8 0 13
Province (autre) 3 124 23 150
Région parisienne 1 13 49 63
total 9 145 72 226

Classiquement, on peut retraduire ces chiffres (qui sont des comptages) sous la forme de probabilités,

Rennes Province
(autre)
Région
parisienne
total
Rennes 2,21% 3,54% 0,00% 5,75%
Province (autre) 1,32% 54,86% 10,18% 66,37%
Région parisienne 0.44% 5,75% 21,68% 27,88%
total 3,98% 64,15% 31,86% 100,00%

Les probabilités qui apparaissent sur les totaux sont les probabilités marginales. Aussi, au bout de la première ligne on trouve que 5,75% des doctorants qui ont eu un poste en 2008 avaient soutenu leur thèse à Rennes (ou dans la région, i.e. les universités de Rennes 1, de Rennes 2 et de Brest). De manière similaire, en bas de la première colonne, 3,98% des maîtres de conférences recrutés en 2008 ont été recrutés à Rennes.
Deux autres tableaux classiques sont alors les profils ligne et profils colonnes, c’est à dire les probabilités conditionnelles,

Rennes Province
(autre)
Région
parisienne
total
Rennes 38,46% 61,54% 0,00% 100,00%
Province (autre) 2,00% 82,67% 15,33% 100,00%
Région parisienne 1,58% 20,64% 77,78% 100,00%
total 3,98% 64,15% 31,86% 100,00%

par ligne, autrement dit, pour un étudiant qui a soutenu sa thèse à Rennes (et qui a eu un poste de Maître de Conférences), il avait  38,46% de chance de rester à Rennes (ou plutôt en Bretagne), et 61,54% de chance d’aller dans une autre ville de province. Côté colonnes,

Rennes Province
(autre)
Région
parisienne
total
Rennes 55,56% 5,52% 0,00% 5,75%
Province (autre) 33,33% 85,52% 31,94% 66,37%
Région parisienne 11,11% 8,96% 68,06% 27,88%
total 100,00% 100,00% 100,00% 100,00%

Ici, on lit que pour un Maître de Conférences recruté sur la région parisienne, il y avait 68,06% de chances pour qu’il ait fait en région parisienne.
Bref, à partir de ces tableaux, on peut essayer de tester l’hypothèse d’indépendance, et si cette dernière n’est pas vérifiée, on va essayer de comprendre pourquoi.
Mais avant, rappelons que l’hypothèse d’indépendance se traduit par le fait que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/cont-03.png

La valeur de gauche est la probabilité jointe (que l’on retrouve au coeur du second tableau), alors que les deux autres sont les probabilités marginales, qui apparaissent en bas, et à droite (dans les totaux). Bref, on pourrait construire le tableau des probabilités que l’on devrait avoir sous hypothèse d’indépendance,

Rennes Province
(autre)
Région
parisienne
total
Rennes 0,22% 3,69% 1,83% 5,75%
Province (autre) 2,64% 42,58% 21,14% 66,37%
Région parisienne 1,11% 17,88% 8,88% 27,88%
total 3,98% 64,15% 31,86% 100,00%

que l’on peut retraduire en terme d’effectif espéré, toujours sous l’hypothèse d’indépendance. Pour ceux qui connaissent les copules, on change juste la structure de dépendance ici, pas les lois marginales (ou les effectifs marginaux) qui restent inchangés.

Rennes Province
(autre)
Région
parisienne
total
Rennes 0,5 8,3 4,2 13
Province (autre) 6,0 96,2 47,8 150
Région parisienne 2,5 40,5 20 63
total 9 145 72 226

L’outils de base pour comparer ces données est la distance du chi-deux.
La statistique du chi-deux est

https://perso.univ-rennes1.fr/arthur.charpentier/latex/cont-04.png

ou encore

https://perso.univ-rennes1.fr/arthur.charpentier/latex/cont-05.png

Karl Pearson avait été un des premiers à étudier ces grandeurs, sous des hypothèses de normalité (des sommes de lois normales indépendantes au carré se ramenant à des lois du chi-deux).
On notera qu’une approche alternative est de supposer que les nombres (ou les effectifs) suivent une loi de Poisson. On peut alors se demander si

https://perso.univ-rennes1.fr/arthur.charpentier/latex/cont-06.png

Si on suppose avoir une loi de Poisson, on peut s’intéresser aux résidus de Pearson,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/cont-07.png

autrement dit, on retrouve des grandeurs très proches de la contribution au chi-deux (en prenant le carré). On aussi le signe qui nous indique s’il y a trop peu de monde, ou au contraire par assez….
Si je travaille sur la base totale, soit environ 25 “villes“, on obtient

> cont1=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/contingence-universite_1.csv",
+ header=TRUE,sep=";")
> cont2=cont1
> for(i in 1:ncol(cont2)){
+ for(j in 1:nrow(cont2)){
+ if(is.na(cont2[j,i]==TRUE)){cont2[j,i]=0} } }
> X=as.matrix(cont2[1:nrow(cont2),2:ncol(cont2)])
> rownames(X)=as.character(cont2[,1])
> colnames(X)=names(cont2[2:ncol(cont2)])
> X=X[,-2]
> X=X[-2,]
> X=X[,-9]
> X=X[-9,]
> Z=X[-nrow(X),]
> Z=Z[,-ncol(Z)]
> Z=Z[-nrow(Z),]
> chisq.test(Z)
        Pearson's Chi-squared test
data:  Z 
X-squared = 1342.546, df = 576, p-value < 2.2e-16

La base est assez sale, j’ai fait un copier/coller rapide du fichier pdf (et je vire les Antilles et la Réunion car après je voudrais faire un peu de visualisation sur une carte). Autrement dit, l’hypothèse d’indépendance ne semble pas validée par le test du chi-deux. Essayons maintenant de comprendre un peu mieux les interactions dans notre carré…
Si on regarde les problèmes de recrutement local, il faut regarder ce qui se passe sur la diagonale de la matrice des résidus,

> CT=chisq.test(Z)
Message d'avis :
In chisq.test(Z) : l'approximation du Chi-2 est peut-être incorrecte
> Dia=diag(CT$residuals)
> names(Dia)=rownames(Z)
> Dia
 AIX-MARSEILLE        BESANCON        BORDEAUX            CAEN        CLERMONT  
     6.43901593      3.26116568      8.37261202     10.53607373      5.12692873 
         CORSE           DIJON        GRENOBLE           LILLE         LIMOGES  
     8.56426340      4.10931024      9.18869705      7.59456681      8.56426340 
          LYON     MONTPELLIER           NANCY          NANTES            NICE  
     3.16119625      4.18366071      8.84765951      2.67087742      7.86973591 
       ORLEANS        POITIERS           REIMS          RENNES           ROUEN  
     1.79480842      5.81144228     -0.06651901      6.22963060      6.57435350 
    STRASBOURG        TOULOUSE           PARIS         CRETEIL      VERSAILLES  
     5.64850464      4.12289099      5.46092001     -0.22061860      4.01351350

Visuellement, sur une carte de France, on obtient la distribution ci-dessous, pour les postes en économie et droit,

et ci-dessous pour les postes en sciences et pharmacie,

Si on regarde le tableau des résidus on obtient les valeurs suivantes sur la ligne correspondant à Rennes, i.e. les endroits où vont les étudiants ayant soutenu leur thèse à Rennes,

avec en bleu les valeurs positives, et en rouge les valeurs négatives. A droite sont représentées les lignes, i.e. les endroits d’où viennent les Maîtres de Conférences recrutés à Rennes.
Au lieu de regarder ville par ville, on peut aussi faire une petite animation, soit sur les endroits où sont partis les doctorants,

ou les endroits d’où viennent les Maîtres de Conférence recrutés,

On peut également s’il existe des interactions, ou des échanges d’une ville à l’autre. Dans ce cas, les résidus croisées https://perso.univ-rennes1.fr/arthur.charpentier/latex/cont-08.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/cont-09.png doivent être positifs.

Bref, au delà du “recrutement local” on voit également apparaître de fortes interactions régionales (Rennes-Nantes, Toulouse-Bordeaux, Orléans-Paris).

Avant de conclure, il existe une autre technique permettant de mieux visualiser les interactions, c’est l’analyse factorielle des correspondances (je renvois icipour des slides). Si on regarde les 2 premiers axes principaux, en sciences (à gauche) et en droit et économie (à droite), on obtient les projections suivantes des modalités, avec en bleu les villes où les étudiants ont soutenu leur thèse, enrouge les villes où les postes de Maître de Conférences ont été obtenus,

et sur les les 2 axes suivants

Bref, comme souvent en statistique, afin de mieux comprendre un ajustement de modèle (ou comme ici lors la réalisation d’un test d’hypothèse) la visualisation des résidus est la technique de base (c’était l’idée d’ailleurs du commentaire que j’avais fait ).
Voilà pour la première partie, sur les modèles où l’on suppose l’indépendance. La prochaine fois, on verra d’autres formes d’hypothèses.

Analyse Factorielle des Correspondances

Les slides de la seconde partie sont en ligne sur mon blog.

Parmi les cours en ligne, on trouvera des notes de cours au CNAM,  ou de Jean Michel Loubes. Sinon, pour accéder aux bases de données, les codes R sont les suivants,

base1=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/Diplomes-emplois-1.csv",header=TRUE,sep=";")
base2=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/Diplomes-emplois-2.csv",header=TRUE,sep=";")
base=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/election2002.txt",header=TRUE)
base=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/csp-vac.csv",header=TRUE,sep=";")

Enfin, les bases de données qui seront utilisées en TD sont les suivantes,

base=read.table("https://perso.univ-rennes1.fr/arthur.charpentier/election.csv",header=TRUE,sep=";")