# Heuristics on Correspondance Analysis

This week, in the course on non-supervised techniques for data science, we’ve been using a dataset, with candidate for the presidential elections in 2002 (per row) and newpapers (per column). In order to visualize that dataset, consider three candidates, and three newspapers

> base=read.table(
> sb=base[,c(2,3,4)]
> sb=sb[c(4,12,7),]
> (N=sb)
LeFigaro Liberation LeMonde
Jospin        7         41      26
Chirac       35          9      18
Mamere        1         10       7

The first part is based on a description of rows. Consider here rows are conditional probabilities, in the set of newspapers,

> (L=N/apply(N,1,sum))
LeFigaro Liberation   LeMonde
Jospin 0.09459459  0.5540541 0.3513514
Chirac 0.56451613  0.1451613 0.2903226
Mamere 0.05555556  0.5555556 0.3888889

The “average row” is the marginal distribution of newspapers

> (Lbar=apply(N,2,sum)/sum(N))
LeFigaro Liberation    LeMonde
0.2792208  0.3896104  0.3311688

If we visualize those individuals, in the set of newpapers (in the simplexe in the newspapers space), we have Here it is, But actually, we will not stay in the simplexe. A PCA is considered, with weights on individuals, that take into account the importance of the different candidates, and weights for the scalar product (in order to have a distance related to the chi-square distance, and not a standard Euclidean distance)

> matL0=t(t(L)-Lbar)
> library(FactoMineR)
> acpL=PCA(matL0,scale.unit=FALSE,
+   row.w=(apply(N,1,sum)),
+   col.w=1/(apply(N,2,sum)))
> plot.PCA(acpL,choix="ind",ylim=c(-.02,.02)) The second part is based on a description of columns. Here Columns are conditional probabilities, in the set of candidates,

> (C=t(t(N)/apply(N,2,sum)))
LeFigaro Liberation   LeMonde
Jospin 0.16279070  0.6833333 0.5098039
Chirac 0.81395349  0.1500000 0.3529412
Mamere 0.02325581  0.1666667 0.1372549

Here again, we can compute the “average column”

> (Cbar=apply(N,1,sum)/sum(N))
Jospin    Chirac    Mamere
0.4805195 0.4025974 0.1168831

In the simplex, points are i.e. But here again, we won’t use that simplexe. We consider a PCA, with two vectors of weights, some to take into account the weights of the newspapers, and some to get a chi-square distance

> Cbar=apply(N,1,sum)/sum(N)
> matC0=C-Cbar
> acpC=PCA(t(matC0),scale.unit=FALSE,
+          row.w=(apply(N,2,sum)),
+          col.w=1/(apply(N,1,sum))) Now, we can almost overlap the two projections. Almost because we might, sometime switch right and left, top and bottom. Because if  is a (unit) eigenvector, so is . Here, for instance, we should swich them

> CA(N) # Séries Temporelles, Exercices

Quelques exercices de séries temporelles, histoire de travailler un peu le cours… En TP, on travaillera davantage les aspects computationnels.

# Inequalities and Poverty #3

The third part of the graduate course on inequalities and povertywill be based on the slides below

# Contraintes sur les autocorrélations d’une série temporelle stationnaire

En cours de séries temporelles, ce matin, on a fini avec les contraintes sur les auto-corrélations, compte tenu du fait que la matrice d’autocorrélations doit être positive. On l’a vu, la matrice d’autocorrélations et la matrice de Toeplitz de déterminant

Avec la règle qu’on apprend tout petit pour calculer les déterminants 3×3, on a la somme des produits des termes pourpres, bleus et verts

à laquelle on retranche la somme des produits des termes verts, bleus et oranges,

Bref, on obtient assez facilement que

soit

avec de plus les corrélations dans l’intervalle , La région admissible pour la paire de corrélations  est alors la région bleue ci-dessous # Niort, Data Day

Lundi et mardi, je serais à Niort où la MAIF organise un Data Day avec plusieurs conférences (tables rondes) autour de la donnée

L’arrivée du digital créée de nouvelles opportunités, de nouveaux besoins, et de nouveaux risques. Dans ce contexte fortement évolutif et très complexe comment les personnes que nous sommes vont-ils évoluer autour de l’information de son usage et de sa gouvernance ? Quelles tendances peut-on observer ? Quelles compétences à quels horizons et pour quelles attentes ? # 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",
> 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) # Les Espérances de Vie (à la naissance) Ce soir, nous avions une discussion passionnante avec @@ et @ au sujet d’un article paru dans Le Monde Comme souvent, seul un point très précis est retenu d’une longue étude sur la démographie française (ici l’espérance de vie à la naissance) publié par l’INSEE. Au delà du débat sur les chiffres, j’ai pu comprendre ce que l’INSEE mettait derrière le terme “espérance de vie à la naissance” à savoir “la durée de vie moyenne d’une génération fictive qui connaîtrait tout au long de son existence les conditions de mortalité par âge de l’année considérée“. En fait, je n’avais pas compris jusqu’à ce soir ce que cachait cette idée de “génération fictive soumise aux conditions de mortalité de l’année“. En formalisant un peu, ça veut dire que les calculs sont fait de manière transversale, et pas du tout prospective. Car pour faire simple, des données de mortalité, ce sont des données sous la forme suivante, A différentes dates, disons qu’on fait un recensement. On a alors le nombre de personnes en vie, et décédée à un âge donné. Par exemple, $p_{x,t}$ désigne ici la probabilité qu’une personne d’âge $x$ survive jusqu’à son prochain anniversaire, l’année $t$. Un peu de calculs permet de voire que $$e= \sum_{x=0}^{120} \prod_{n=0}^x p_{n,t}$$peut être interprété comme une espérance de vie à la naissance. C’est la version discrète de la formule $$\mathbb{E}[T]=\int_0^\infty \mathbb{P}(T>t)dt$$(quand on travaille sur une variable aléatoire positive) pour la somme, et le produit est juste lié au fait que $$\mathbb{P}[T>t]=\mathbb{P}[T>t\vert T>t-1]\cdot \mathbb{P}[T>t-1\vert T>t-2] \cdots$$ (avec des mots, la probabilité d’atteindre 60 ans, c’est la probabilité de survivre à son 1er anniversaire, puis à son second, …, puis à son 59ème). Bref, on a toutes ces probabilités de survie dans nos données. Des données, on peut en trouver par exemple sur mortality.org. > library(demography) > france <- hmd.mx("FRACNP", username, password, "France") La méthode décrite par l’INSEE ressemble à la suivante: On considère toutes les probabilités de survie brute, et on cumule, > kPx_obs=1-france$rate$total[,"2013"] > sum(cumprod(kPx_obs))  81.24053 (je ne fais pas ici de distinction Homme/Femme, les plus curieux pourront creuser par eux-mêmes). On note que si en 2013, il y a eu une grosse grippe, qui a touché tout le monde, l’espérance de vie en sera forcement impactée ! On n’a pas ici de correction ici, pour lisser les décès tout au long de l’année, mais c’est possible de le faire. Une autre stratégie peut être d’ajuster un modèle sur les données de mortalité, par exemple le modèle de Lee & Carter (1992), qui dit que la probabilité qu’une personne d’âge $x$ décède avant son prochain anniversaire, l’année $t$, est de la forme $$1-p_{x,t}=\exp[a_x+b_x\cdot \kappa_t]$$ On peut ajuster ce modèle et considérer alors $$\widehat{p}_{x,t}=1-\exp[\widehat{a}_x+\widehat{b}_x\cdot \widehat{\kappa}_t]$$ On peut utiliser une version transversale comme auparavant, > france.LC <- lca(france,adjust="e0",years=1950:2013) > france.fcast <- forecast(france.LC,jumpchoice="actual",h=150) > Ax=france.LC$ax
> Bx=france.LC$bx > Kt=c(france.LC$kt,france.fcast$kt$mean)
> years=c(1950:2013,2014:2163)

L’espérance de vie à la naissance (telle qu’on l’avait définie auparavant) est alors

> kPx_model=1-exp(Ax+Bx*Kt[years==2013])
> sum(cumprod(kPx_model))
 84.27459

J’avoue être ici surpris par l’écart avec le précédant. Mais pourquoi pas…

Enfin, autre méthode, basée sur un modèle prospectif. L’idée est de noter que le vieillissement d’un individu, dans une cohorte donnée, se fait en diagonale, et pas verticalement, dans le diagramme de Lexis On utilise ici le fait que quelqu’un né en 2015 autant 40 ans en 2055, et qu’il faut prévoir la probabilité de décéder l’année de ses 40 ans dans 40 ans. En intégrant les gains sur l’espérance de vie, on obtient cette fois l’espérance de vie à la naissance

> kPx_pred=1-exp(Ax+Bx*Kt[years%in 13:(2012+length(Ax))])
> sum(cumprod(kPx_pred))
 90.15392

pour quelqu’un né en 2013. Cette fois, on a un gain vraiment significatif, correspondant à des améliorations des conditions de vie avec le temps, et la baisse de la mortalité.

Je reste toutefois surpris que l’INSEE continue à utiliser une version transversale (et pas du tout prospective) dans ses calculs d’espérance de vie, même si, comme le notait @l’espérance de vie à la naissance est une statistique internationale publiée par tous les pays, et donc doit être calculée de façon simple (ou pas trop complexe). Et je renvoie vers un vieux billet pour ceux qui veulent aller plus sur ces calculs prospectifs, et le fameux “tous les ans, on gagne un trimestre d’espérance de vie“. Ou encore un article sur les tables de mortalité, paru dans Risques en 2007.

# Décomposition(s) de Séries Temporelles

La semaine dernière, pour le TP du cours de série temporelle, on avait parlé de décomposition de séries temporelles. On part d’une série temporelle , et on va considérer une décomposition de la forme

par exemple une tendance linéaire, et un cycle annuel,

Considérons la série suivante

> autoroute=read.table(
+ "http://freakonometrics.free.fr/autoroute.csv",
>  a7=autoroute$a007 > X=ts(a7,start = c(1989, 9), frequency = 12) correspondant à du trafic sur l’autoroute (A7) > plot(X,xlim=c(1989,2000)) Pour estimer la tendance, on va faire une régression linéaire, sur le temps, à partir de la base suivante > T=time(X) > base=data.frame(Y=a7,X1=as.numeric(T)) Le modèle est ici > reg1=lm(Y~X1,data=base) ce qui permet de faire une prévision sur les données passées, mais aussi futures, > tendance1=predict(reg1) > T=as.numeric(T) > lines(T,tendance1,col="red",lwd=3) > futur=seq(max(T)+1/12,length=36,by=1/12) > tendance2=predict(reg1, + newdata=data.frame(X1=futur)) > lines(futur,tendance2,col="red") Voilà pour la tendance. Pour le cycle, on va ajuste un modèle de la forme suivante sur la série des résidus > base$Y2=residuals(reg1)

Cette fois, on va définir le mois de la date

> base$mois=round(base$X1-t
>  polygon(P2$hull.loop,col="orange",border=NA) > segments((A+C)/2,(A+C)/2,B,B,col="grey",lty=2) > segments((A+B)/2,(A+B)/2,C,C,col="grey",lty=2) > segments((B+C)/2,(B+C)/2,A,A,col="grey",lty=2) > points(tripoint(b),pch=19,cex=2,col="red") # Inequalities and Poverty #2 The second part of the graduate course on inequalities and poverty will be based on the slides below # Multiattribute, Risk and Utility As discussed already, in June 2016, with Olivier L’Haridon, we will organize a (small) conference, in Rennes, on risk models in a multi-attribute framework. Related to that conference, we have a working group on related topics. At the end of September, I gave a brief survey on multivariate distributions. And yesterday, Olivier gave the first part of survey on multivariate decision making. The second part will be in two weeks, # Inequalities and Poverty #1 This Monday, we will start the graduate course on inequalities and poverty. Slides are below # Collectionner les Figurines Lego avec ses Ami(e)s Dans mon précédant billet, je regardais combien de figurines Lego je devais acheter pour avoir toute la série, avec ma fille. Pour la série de 16, il fallait, en moyenne, acheter 54 sachets, à 3 euros pièces. Dans les commentaires, @Guillaume et @Gaelle me faisaient noter que les échanges, ça pouvait aider. Bon, tout d’abord j’aime pas qu’on donne des conseils à mes filles pour avoir raison contre moi ! Parce que oui, c’est sûr qu’en s’y mettant à plusieurs, le coût individuel sera plus faible… Ensuite, parce que les maths me semblent plus compliqués…. Cela dit, en simulant des parties, ça tombe assez vite… Pour commencer, je fais une petite fonction (histoire d’alléger mes codes) pour incrémenter mes comptages de figurines, à chaque fois que j’en achète une > f=function(mat,i){ + mat[i]=mat[i]+1 + return(mat) + } Ensuite, on peut avoir une petite fonction, correspondant à une simulation, qui me renvoie le nombre de figurine de chaque type que j’ai du acheter pour avoir les 16 une fois. Ou plus. En fait, je veux avoir en tout k fois toutes les figurines (ce qui permet de partager parmi k ami(e)s) > tirage=function(k=1){ + test=FALSE + cpte=0 + M=rep(0,16) + while(!test){ + no=sample(1:16,size=1) + M=f(M,no) + cpte=cpte+1 + if(min(M)==k) test=TRUE + } + return(list(compte=M, + nbre=cpte)) + } Maintenant on peut faire des simulations. Par exemple, pour retrouver le fait qu’il faut, en moyenne, acheter 54 sachet pour que ma fille (seule) ait toutes les figurines > N=unlist(lapply(1:10000,function(x) tirage(1)$nbre))
> mean(N)
 54.1854

Si maintenant on le fait avec mes deux filles,

> N=unlist(lapply(1:10000,function(x) tirage(2)$nbre)) > mean(N)  82.7961 ce qui est (beaucoup) moins que le double. Si au lieux de faire ça entre mes deux filles, mais que ma fille collectionne les figurines avec une amie à elle (et que les papas partagent les frais), on a en moyenne 41 sachets à acheter, chacun, > mean(N)/2  41.39805 et pour trois, on passe à 36 sachets > N=unlist(lapply(1:10000,function(x) tirage(3)$nbre))
> mean(N)/3
 36.04363

etc

> V=unlist(lapply(1:10,function(k) mean(unlist(lapply(1:10000,function(x) tirage(k)$nbre))))) > plot(1:10,V/1:10,type="b",pch=19) On peut aller plus loin, en regardant les quantiles. Avec 4 ami(e)s, si on achète tous 45 sachets, on a 90% de chances que tout le monde ait les 16 figurines D’ailleurs, si on regarde l’évolution du quantile (à 90%) sur la moyenne, on voit qu’il baisse vite avec la taille…. Avec 2 ami(e)s, il faut 35% de sachets en plus pour s’assurer d’avoir toutes les figurines dans 90% de cas (par rapport au nombre moyen) alors qu’il en faut juste 20% de plus avec 8 ami(e)s. Maintenant, si quelqu’un a le courage (et le temps) de faire les maths, je suis preneur ! # TP, Séries Temporelles Deux bases pour illustrer le TP de séries temporelles. La première est du trafic sur autoroute, en France > autoroute=read.table( + "http://freakonometrics.free.fr/autoroute.csv", + header=TRUE,sep=";") > a7=autoroute$a007
> X=ts(a7,start = c(1989, 9), frequency = 12)
> T=time(X) et on utilisera aussi du traffic sur le réseau de la SNCF

> sncf=read.table(
+   "http://freakonometrics.free.fr/sncf.csv",
> train=as.vector(t(as.matrix(sncf[,2:13])))
> X=ts(train,start = c(1963, 1), frequency = 12) # 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",
>  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