Tag Archives: computer

Optimal Portfolios #1

This afternoon, I will start a crash course on financial portfolio optimization, with application in R. This week, we start with simple things, with the theoretical setup, without and with a risk free asset. We will discuss then the problem of estimating parameters, in a robust way. Then we introduce the idea of consider a more general criteria to quantify risk than the variance (but it means more general distributions… this point will be discussed further next time). The slides are available here, and R codes from there (in a Markdown)

Traffic Flow of Kota Kinabalu (with R)

This morning, we had our first practicals on network flows, using  an example mentioned in some papers published by Noraini Abdullah and Ting Kien Hua, max flow min cut theorem to minimize traffic congestion in Kota Kinabalu and application of the Shortest Path and Maximum Flow with Bottleneck in Traffic Flow of Kota Kinabalu. From the roads mentioned in the articles, I did try my best to locate the nodes on a map,

m=matrix(c(0,5.995910, 116.105520,
1,5.992737, 116.093718,
2,5.992066, 116.109883,
3,5.976947, 116.095760,
4,5.985766, 116.091580,
5,5.988940, 116.080112,
6,5.968318, 116.080764,
7,5.977454, 116.075460,
8,5.974226, 116.073604,
9,5.969651, 116.073753,
10,5.972341, 116.069270,
11,5.978818, 116.072880),3,12)

we can be visualized below

map = openmap(c(lat= 6.000, lon= 116.06),
c(lat= 5.960, lon= 116.12))
points(t(m[3:2,]),col="black", pch=19, cex=3 )

If the source is realistic (up north), I do not feel very confortable with the location of the sink (on the west). But let’s pretend it’s find (to do the maths, at least).

To extract information about edge capacity, on that network use the following code that will extract the three tables from the paper

location <- 'http://www.jistm.com/PDF/JISTM-2017-04-06-02.pdf'
out <- extract_tables(location)

with Windows, it seems to be necessary to download another package first

location <- 'http://www.jistm.com/PDF/JISTM-2017-04-06-02.pdf'
out <- extract_tables(location)

Now we can get out data frame with capacities


We can add those edges on our map (without the arrows to indicate the direction, it would be to heavy to read)

points(t(m[3:2,]),col="black", pch=19, cex=3 )
for(i in 1:nrow(E)){

To get the graph with capacities, an alternative is to use


but it does not respect geographical locations of nodes. It can actually be done using

plot(g, layout=as.matrix(B[,c("x","y")]))

To get a better understanding of the capacities of the road, use

plot(g, layout=as.matrix(B[,c("x","y")]),

From that network with capacities, the goal is to determine maximum flow on that network, from the source to the sink. This can be done with R using

> (m=max_flow(graph=g, source="s", target="t"))
[1] 2571

[1] 1191 1380 1422 1380 231 0 231 0 1149 1422 1149 0 0 1149 1422
[16] 1149

Our maximum flow is here 2571, which is different from was is actually claimed both in the two papers  max flow min cut theorem to… and application of the Shortest Path… (“the maximum flow for the capacitated network with 12 nodes and 16 edges of the selected scope in this study was 2598 vehicles per hour“) where there are clearly typos since values in the table and on the graph are different. Here I did use the ones from the tables.

plot(g, layout=as.matrix(B[,c("x","y")]),

That is nice, but rather odd. Actually, a much simpler flow can be considered, but the same global value

plot(g, layout=as.matrix(B[,c("x","y")]),

Nice, isn’t it. It is actually possible to do exactly the same on another paper they have, on the same city, traffic congestion problem of road networks in Kota Kinabalu.

location <- 'http://www.worldresearchlibrary.org/up_proc/pdf/999-150486366625-30.pdf'
out <- extract_tables(location)

Here the value of the maximal flow is 4017, just as they found in the original paper

Multinomial Logit as an Iterated Logit Regression

For the second section of the course at ENSAE, yesterday, we’ve seen how to run a multinomial logistic regression model. It is simply an extension of the binomial logistic regression. But actually, it is also possible to consider iterative binomial regressions.

Consider here a response variable Y with a multinomial distribution (3 factors to have something more general than the binomial), taking values \{A,B,C\}, with respective probabilities \mathbf{p}=(p_A,p_B,p_C). Here is a code to generate some multinomial variables

for(i in 1:B){Y[i]=sample(A,size=1,prob=C[i,])}

and here is a code to generate a dataset with n rows,


Let us generate a training dataset and a validation one


With a multivariate logistic regression

For convenience, consider the most popular factor in our training dataset

modalite=names(sort(table(DF1$Y),decreasing = TRUE))

Consider a regression model on the simulated dataset (with several covariates), let us estimate it, and let us get predictions.

reg=multinom(as.factor(Y) ~ ., data = DF1)
mp1=predict (reg, DF1, "probs")
mp2=predict (reg, DF2, "probs")

An alternative can be the following.
consider a first regression model on the Bernoulli variable Y_A=\mathbf{1}(Y=A). Actually, we will consider the most important factor, but for convenience, assume that it is A.
On our dataset, estimate that model, and get predictions. In the case where Y\neq A, define another Bernoulli variable Y_B=\mathbf{1}(Y=B|Y\neq A). We can estimate that model and derive two probabilities, \mathbb{P}(Y=B|Y\neq A) and \mathbb{P}(Y=C|Y\neq A) (the sum of the two being equal to 1). Based on those two models, it is possible to compute the three probabilities we are looking for. \mathbb{P}[Y=A] is obtained from the first model, and we can derive the other two from \mathbb{P}[Y=B|Y\neq A]\cdot\mathbb{P}[Y\neq A] and \mathbb{P}[Y=C|Y\neq A]\cdot\mathbb{P}[Y\neq A].

p11=predict (reg1, newdata=DF1, type="response")
p12=predict (reg2, newdata=DF1, type="response")
p21=predict (reg1, newdata=DF2, type="response")
p22=predict (reg2, newdata=DF2, type="response")

Let us compare the predicted probabilites, on the same dataset (here the training dataset)

> mmp1[1:9,c("0","1","2")]
0 1 2
1 0.19728737 0.4991805 0.3035321
2 0.17244580 0.5648537 0.2627005
3 0.19291753 0.5971058 0.2099767
4 0.09087176 0.7787304 0.1303978
5 0.23400225 0.4083022 0.3576955
6 0.18063647 0.6637352 0.1556283
7 0.13188881 0.7402710 0.1278401
8 0.13776970 0.6524959 0.2097344
9 0.12325864 0.6790336 0.1977078
> mp1[1:9,c("0","1","2")]
0 1 2
1 0.19691036 0.5022692 0.3008205
2 0.17123189 0.5680647 0.2607034
3 0.19293066 0.5984402 0.2086291
4 0.08821851 0.7813318 0.1304497
5 0.23470739 0.4109990 0.3542936
6 0.18249687 0.6602168 0.1572863
7 0.13128711 0.7400898 0.1286231
8 0.13525341 0.6553618 0.2093848
9 0.12090016 0.6815915 0.1975084

The two are very close. So yes, it is possible to see the multinomial regression as some sequential binomial regressions.

A la recherche des homonymes…

Il y a quelques mois, Baptiste Coulmont m’avait contacté avec une question passionnante (comme chaque fois qu’il me contacte).

Baptiste voulait travailler sur la proportion de personnes qui ont un homonymes dans une population de taille donnée, ou sur la probabilité de ne pas avoir d’homonymes dans un bureau de vote, par exemple. Or ce dernier problème n’est pas sans rappeler le “paradoxe des anniversaires“.

Dans un groupe de 23 personnes, il y a plus d’une chance sur deux pour que deux personnes de ce groupe aient leur anniversaire le même jour. Ou dit autrement, il est plus surprenant de ne pas avoir deux personnes qui ont leur anniversaire le même jour que d’avoir deux personnes qui ont leur anniversaire le même jour.

On parle de “paradoxe” parce que le nombre de jours possibles est relativement grand (k=365), mais que la taille du groupe est elle relativement petite (n=23). Ici c’est un peu la même chose, sauf qu’on est sur des nombres beaucoup plus grands: on a des milliers de noms possibles (en fait, on s’intéresse au couple prénom+nom, ce qui augmente très largement le nombre de possibilités, que l’on notera k ici), on considère des groupes relativement grands (la taille du groupe sera elle notée n), et surtout, l’hypothèse (implicite) d’uniformité dans la distribution des dates ne marche pas sur les noms de famille, et les prénoms. On peut malgré tout tenter de formaliser le problème.

Considérons une collection de k éléments, notée \mathcal{X}=\lbrace x_1,\cdots,x_k\rbrace, et faisons n tirages avec remise. Les tirages obtenus sont notés ici X_1,\cdots,X_n et correspondent à des variables aléatoires iid, suivant des lois multinomiales sur \mathcal{X} avec comme vecteur de probabilités \mathbf{p}=(p_1,\cdots,p_k). Dans le cas des anniversaires, \mathcal{X} correspond aux dates, avec k=365, et \mathbf{p} est uniforme sur \mathcal{X} – ou encore \mathbf{p}=\mathbf{1}/k . Le problème classique est de calculer la probabilité que dans un groupe de taille donnée tout le monde ait des dates d’anniversaire différentes

\mathbb{P}[\forall j\neq i, X_i\neq X_j]

mais on peut aussi calculer (comme le voulait Baptiste) la proportion de personne ayant un alter-ego (une personne dans le groupe qui a son anniversaire le même jour)
\frac{1}{n}\sum_{i=1}^n \mathbf{1}(\exists j\neq i,X_i=X_j)

On a donc tenté d’approcher le second, quand on a un groupe de taille n. Notons Z_i le nombre de personnes (parmi les n-1 restantes) qui partagent la même date d’anniversaire que i,
Z_i=\sum_{j\in\{1,\cdots,n\}\backslash\{i\}} \mathbf{1}(X_j \neq X_i)
Aussi, la personne i a un alter-ego si Z_i\geq 1. Et donc, la proportion de personnes avec un alter-ego dans le groupe s’écrit
P_{n,k}=\frac{1}{n}\sum_{i=1}^n \mathbf{1}(Z_i\geq 1)
Si \mathbf{p} correspond à la loi uniforme sur \mathcal{X}, alors
1-P_{n,k}=\binom{k}{n}\frac{1}{k^n}\sim \exp\left[\binom{n+1}{2}\frac{1}{k}\right]\sim e^{-n^2/2k}
comme le montre Chatterjee et al. (2004). Je pourrais aussi renvoyer aux articles de Inoue & Aki (2008) et plus récemment de  Cortina Borja (2013) pour des revues de littérature sur les aspects computationnels.

Mais la plupart des résultats sont obtenu dans le cas où \mathbf{p} correspond à une distribution uniforme. Les propriétés plus générales sont beaucoup plus rares(on pourra lire Munford (1977), Das Gupta (2004),  Inoue & Aki (2008) ou Nunnikhoven (1992) pour des tentatives. Utiliser des méthodes de Monte Carlo est probablement le plus naturel (si on n’est pas trop pressé) en se donnant un ensemble \mathcal{X} (qui sont les dates dans le problème des anniversaires et les noms dans le cas des homonymes) de taille k, et une distribution \mathbf{p} sur \mathcal{X}. On a ainsi pu approcher P_{n,k}, la proportion de personnes ayant un alter-ego dans une population de taille n: sur les figures ci-dessous, on a cette proportion en fonction de  n (en abscisse), avec des courbes différentes pour différents k. Histoire d’être un peu réaliste, on va chercher le taux de personnes ayant un homonyme dans une popoluation de plusieurs centaine de milliers de personnes, et quand aussi des milliers, des dizaines de milliers, de possibilités.

Dans le cas d’une distribution uniforme, la proportion évolue comme sur la figure ci-dessous

On a tenté de quitter l’indépendance en utilisant une loi paramétrique, plus précisément une loi de type Pareto, ou Zipf, pour \mathbf{p}. En fait, c’est le type de distribution qu’avait observé Li (2012) sur la distribution des prénoms, avec p_i \propto i^{-\alpha} (ou plus récemment, Höhle (2017)). On a alors l’évolution suivante

On est ici en échelle logarithmique pour la taille n parce que c’est naturel. Mais pour l’instant, on était resté sur P_{n,k}. Baptiste insistait pour passer aussi à une échelle logarithmique sur l’axe des ordonnées, mais ca n’etait pas concluant. On a alors tenté deux autres représentations, plus classiques quand on travaille sur des probabilités (ou des proportions) : représenter \text{logit}(P_{n,k}) en fonction du logarithme de n (où \text{logit}(u)=\log[u/(1-u)]) et une version probit, à savoir \Phi^{-1}(P_{n,k}) (où ici \Phi est la fonctionc cumulative de la loi normale centrée et réduite).

En grande dimension, on a l’impression d’avoir une relation linéaire entre la transformée probit de la proportion de personnes ayant un homonyme et le logarithme de la taille de la population considérée

\Phi^{-1}(P_{n,k}) \sim a_{k}+b_{k} \log[n]

C’est joli, mais on se demandait comment utiliser ca sur nos données, pour travailler vraiment sur les homonymes. Le soucis, c’est qu’on n’a pas de fichier national, avec les noms et les prénoms. Alors Baptiste a sorti de son sac les listes électorales parisiennes. Le soucis est que la base est (relativement) petite. Le nombre de couple est très limité. Au lieu de faire des tirages avec remise dans les couples (“prénom”, “nom”) – encore une fois, la base est trop petite – on a voulu tenter de tirer indépendemment le “prénom” et le “nom”. On y gagne en ayant à disposition une base plus grande, mais on introduit un biais colossal, les “prénom” et “nom” étant assez corrélés. Dans notre cas, on a k_1 “prénom” et k_2 “prénom”. Il y a quelques années Chalabi & Flowers (2014)} (repris plus en détails sur ce blog) était revenu sur la non indépendance des noms et des prénoms, en s’interrogeant sur le nom le plus populaire aux Etats-Unis. Selon les auteurs “Michael Smith” n’était probablement pas le nom le plus populaire, même si “Michael” et “Smith” étaient respectivement les “prénom” et “nom” les plus populaires. En France “Jean” et “Marie” sont les prénoms les plus populaires pour les garcons et les filles, respectivement, et “Martin” est le nom le plus populaire (au moins sur les listes électorales parisiennes), mais en fait, le nom le plus populaire est “Thi Nguyen”, avec 350 personnes à Paris, alors que seulement 84 “Jean Martin” sont observés, le troisième plus haut score. Ah oui, je dois préciser qu’on se limite ici au premier prénom, en cas de prénoms composés. Tout comme Chalabi & Flowers (2014), on a fait un test du chi-deux pour tester l’indépendance, sur les 20 “noms” et “prénoms” les plus populaires. Sur la figure ci-dessous, on a les résidus normalisés du test. Par exemple “Michel Michel”est en rouge, car c’est beaucoup moins donné que ce qu’on aurait si  les deux composantes étaient indépendantes

Bref, on n’a pas trop le droit de tirer le “prénom” et le “nom” indépendament. Mais on va quand même le faire, pour voir. Ici \mathcal{X} est l’ensemble \mathcal{X}_1\times\mathcal{X}_2 avec k_1 choix pour le “prénom” et k_2 pour le “nom”.

Soit \mathbf{p}=(p_{i_1,i_2}) la probabilité empirique sur \mathcal{X}. On imagine que \mathbf{p} contient énormément de 0, car beaucoup de couples ne sont pas observés. Pour les simulations, on va noter \mathbf{p}^{\perp} la pobabilité jointe sous hypothèse d’indépendance
p^{\perp}_{i_1,i_2}=\left(\frac{1}{k_2}\sum_{j=1}^{k_2} p_{i_1,j}\right)\cdot\left(\frac{1}{k_1}\sum_{j=1}^{k_1} p_{j,i_2}\right)=p_{i_1,\cdot}p_{\cdot,i_2}

Lors de nos simulations, on va considérer deux cas : on dire des \mathbf{x}=(x_1,x_2) soit suivant la distribution \mathbf{p}, soit suivant \mathbf{p}^{\perp}. Bon, j’avais formalisé pour clarifier un peu dans ma tête et aussi pour impressionner Baptiste, mais techniquement, dans un cas on tire une paire observée – genre “Baptiste Coulmont” – et dans le second cas, on tire d’un côté un prénom – genre “Baptiste” – et de l’autre un nom de famille – genre “Charpentier” (même s’il n’existe aucun “Baptiste Charpentier”, on dit juste qu’il serait possible qu’il y en ait un).

Bref, on tirant des échantillons (plusieurs milliers de taille donnée, puis en moyennant), on obtient la figure ci-dessous (avec une échelle logarithme pour l’axe des ordonnées)

La courbe noir sur-estime sûrement à cause du faible nombre de choix possibles (on va tirer jusqu’à 50 millions de noms). En revanche, l’indépendance n’étant pas réaliste, la courbe rouge sous-estime probablement la proportion. Mais le plus intéressant est le graphique ci-dessous, obtenu en représentant la version probit de la probabilité

La courbe rouge est incroyablement linéaire ! Bien plus linéaire que ce qu’on espérait avoir suite à nos simulations dans le cas de la loi de Pareto.

En fait, on note que si on se limite au cas où n n’est pas trop grand, la courbe noire semble aussi linéaire. Ensuite on perd la linéarité parce qu’on tire avec remise dans une base trop petite, mais au début (lorsque n n’est pas trop grand), on pourrait accepter l’hypothèse d’une relation linéaire entre la transformée probit, et le logarithme de n. Bref, on estime une relation linéaire, et on extrapole pour les valeurs de n plus grandes

Si on revient sur la représentation plus classique, de la proportion de personnes ayant un homonyme en fonction de n, on obtient

Autrement dit, en France, on peut estimer que 80% des francais ont un homonyme.

Deux points pour poursuivre, sur la robustesse de l’approche. Tout d’abord, si on change un peu la taille de l’intervalle des n pour lesquels la courbe noire semble linéaire, la pente change peu

et on a un ordre de grandeur assez proche. Bref, ce 80% semble assez robuste,

Ensuite, comme cette relation linéaire me surprenait, on a aussi tenté sur des données américaines… et on a exactement le même genre de résultats…

Pour aller plus loin, on a tout mis sur github (avec des graphs en ggplot2 pour faire plus joli…) et une version (en anglais) est en ligne sur arxiv.

Networks with R

In order to practice with network data with R, we have been playing with the Padgett (1994) Florentine’s wedding dataset (discussed in the lecture). The dataset is available from

> library(network)
> data(flo)
> nflo=network(flo,directed=FALSE)
> plot(nflo, displaylabels = TRUE,
+ boxed.labels =

The next step was to move from the network package to igraph. Since we have the adjacency matrix, we can use it

> library(igraph)
> iflo=graph_from_adjacency_matrix(flo,
+ mode = "undirected")
> plot(iflo)

The good thing is that a lot of functions are available, for instance we can get shortest paths, between two specific nodes. And we can give appropriate colors to the nodes that we’ll cross

> AP=all_shortest_paths(iflo,
+ from="Peruzzi",
+ to="Ginori")
> L=AP$res[[1]]
> V(iflo)$color="yellow"
> V(iflo)$color[L[2:4]]="light blue"
> V(iflo)$color[L[c(1,5)]]="blue"
> plot(iflo)

We can also visualize edges, but I found it slightly more complicated (to extract edges from the output)

> liens=c(paste(as.character(L)[1:4],
+ "--",
+ as.character(L)[2:5],sep=""),
+ paste(as.character(L)[2:5],
+ "--",
+ as.character(L)[1:4],sep=""))
> df=as.data.frame(ends(iflo,E(iflo)))
> names(df)=c("src","target")
> lstn=sort(unique(c(as.character(df[,1]),as.character(df[,2]),"Pucci")))
> Eliens=paste(as.numeric(factor(df[,1],levels=lstn)),"--",
+ as.numeric(factor(df[,2],levels=lstn)),sep="")
> EU=unlist(lapply(Eliens,function(x) x%in%liens))
> E(iflo)$color=c("grey","black")[1+EU]
> plot(iflo)

But it works. It is also possible to use some D3js visualization

> library( networkD3 )
> simpleNetwork (df)

Then the next question was to add a vertice to the network. The most simple way to do it is probability through the adjacency matrix

> flo2=flo
> flo2["Pucci","Bischeri"]=1
> flo2["Bischeri","Pucci"]=1
> nflo2=network(flo2,directed=FALSE)
> plot(nflo2, displaylabels = TRUE,
+ boxed.labels =

Then, we’ve been playing with centrality measures.

> plot(iflo,vertex.size=betweenness(iflo))

The goal was to see how related they were. Here, for all of them, “Medici” is the central node. But what about the others?

> B=betweenness(iflo)
> C=closeness(iflo)
> D=degree(iflo)
> E=eigen_centrality(iflo)$vector
> base=data.frame(betw=B,close=C,deg=D,eig=E)
> cor(base)
betw close deg eig
betw 1.0000000 0.5763487 0.8333763 0.6737162
close 0.5763487 1.0000000 0.7572778 0.7989789
deg 0.8333763 0.7572778 1.0000000 0.9404647
eig 0.6737162 0.7989789 0.9404647 1.0000000

Those measures are quite correlated. It is also possible to use a hierarchical graph to visualize how close those centrality measures can be

> H=hclust(dist(t(base)),
+ method="ward")
> plot(H)

Instead of looking at values of centrality measures, it is possible to looks are ranks

> rbase=base
> for(i in 1:4) rbase[,i]=rank(base[,i])
> H=hclust(dist(t(rbase)),
+ method="ward")
> plot(H)

Here the eigenvector measure is very close to the degree of vertices.

Finally, it is possible to seek clusters (in the context of coalition here, in case a war should start between those families)

> kc <- fastgreedy.community ( iflo )

Here we have 3 classes (+1 for the node that is disconnected from the other families)

> V(iflo)$color=c("yellow","orange",
+ "light blue")[membership ( kc )]
> plot(iflo)

> plot(kc,iflo)

I Got The Feelin’

Last week, I’ve been going through my CD collection, trying to find records I haven’t been listing for a while. And I got the feeling that music I listen to nowadays is slower than the one I was listening to in my 20’s. I was wondering if that was an age issue, or it was simply the fact that music in the 90s was “faster” than the one released in 2015. So I tried to scrap the BPM database to get a more appropriate answer about this “feeling” I have. I did extract two information: the BPM (beat per minute) and the year (of release).

Here is a function to extract information from the website,

> library(XML)
> extractbpm = function(VBP,P){
+ url=paste("https://www.bpmdatabase.com/music/search/?artist=&title=&bpm=",VBP,"&genre=&page=",P,sep="")
+ download.file(url,destfile = "page.html")
+ tables=readHTMLTable("page.html")
+ return(tables)}

For instance

> extractbpm(115,13)
Artist Title
1 Eros Ramazzotti y Claudio Guidetti Dimelo A Mi
2 Everclear Volvo Driving Soccer Mom
3 Evils Toy Dear God
4 Expose In Walked Love
5 Fabolous ft. 2 Chainz When I Feel Like It
6 Fabolous ft. 2 Chainz When I Feel Like It
7 Fabolous ft. 2 Chainz When I Feel Like It
8 Fanny Lu Fanfarron
9 Featurecast Ain't My Style
10 Fem 2 Fem Obsession
11 Fernando Villalona Mi Delito
12 Fever Ray Triangle Walks
13 Firstlove Freaky
14 Fito Blanko Pegadito Suavecito
15 Flechazo Del Norte Mariposa Traicionera
16 Fluke Switch/Twitch
17 Flyleaf Something Better
18 FM Static The Next Big Thing
19 Fonseca Eres Mi Sueno
20 Fonseca ft. Maffio & Nayer Eres Mi Sueno
21 Francesca Battistelli Have Yourself A Merry Little Christmas
22 Frankie Ballard Young & Crazy
23 Frankie J. More Than Words
24 Frank Sinatra The Hucklebuck
25 Franz Ferdinand The Dark Of The Matinée
Mix BPM Genre Label Year
1 — 115 — Sony 2009
2 — 115 — Capitol Records 2003
3 — 115 — — —
4 — 115 — Arista Records 1994
5 Explicit 115 Urban Def Jam/Island Def Jam 2013
6 — 115 Urban Def Jam/Island Def Jam 2013
7 Radio Edit 115 Urban Def Jam/Island Def Jam 2013
8 — 115 Latin Pop Universal Latino 2011
9 Psychemagik Dub 115 — Jalapeno 2012
10 — 115 — Critique Records 1993
11 — 115 — Mt&vi Records/caminante Records 2001
12 Rex The Dog Remix 115 — Little Idiot/Mute 2012
13 — 115 — Jwp Music 2000
14 — 115 Merengue Mambo Crown Loyalty 2012
15 — 115 — Hacienda 2010
16 Album Version 115 — One Little Indian Records 2004
17 — 115 Alternative A&M/Octone 2013
18 — 115 — Tooth & Nail Records 2007
19 — 115 Merengue Mambo 10 2012
20 Urban Version 115 — 10 2012
21 — 115 — Word/Fervent/Warner Bros 2009
22 — 115 Country Warner Bros 2015
23 Mynt Rocks Radio Edit 115 — Columbia 2005
24 — 115 Jazz Columbia 1950
25 — 115 New Wave — 2004

We have here one of the few old songs, a 1950 tune by Frank Sinatra. If we scrap the website, with a simple loop (where the bpm is from 40 to 200). Start with

> vbp=40
> p=1

and then, a loop based on

> while(vbp<=2017){
+ F=extractbmp(vbp,p)
+ if(length(F)==1){
+ BASE=rbind(BASE,F[[1]][,c("Artist","Title","BPM","Year")])
+ p=p+1}
+ if(length(F)==0){
+ vbp=vbp+1
+ p=1}}

Then we should clean the dataset


and we end up with almost 50,000 tunes.

col="light blue")

Over the past 20 years, it looks like speed of tunes has declined (let us forget tunes of 2017, clearly we have a problem here…)


which is confirmed with a (smoothed) regression


even when incorporating the confidence band. Bumps are probably related to smoothing parameters, but indeed, it looks like the average speed of music tune has decreased, from 110-115 in the 90’s to less than 100 nowadays. Now to be honest, I would love to have access to personal information from itunes, deezer or spotify, to get a better understanding (eg when in the week, in the day, do we like to listen to faster music for instance). But so far, I could not have access to such data. Too bad…

Les votes à l’assemblée nationale

Un dernier petit billet basé sur les projets de R que j’avais donnés pour la formation Data Science pour l’Actuariat. Aujourd’hui, je reviens sur des codes tirés du projet de Raphaël qui a scrapé les données de l’Assemblée Nationale. On commence par charger les librairies qui nous seront utiles,


La première partie de l’importation des données vise à importer l’ensemble des députés, c’est à dire l’ensemble des personnes  (compris entre 1 et 248) dont le mandat est à l’ASSEMBLEE. Pour chacun de ces députés, on récupére le numéro de circonscription et de département, afin de pouvoir construire le code IDEN, identifiant unique de chaque circonscription (DDCC, deux chiffres pour le département, deux pour la circonscription), la référence acteur (acteurRef) et la référence mandat (uid) afin de pouvoir construire un identifiant pour chaque député du type (PA???XXXPMXXXXXX, chaque X representant un chiffre, les ? representant un chiffre ou l’absence de caractères – les PA n’ont pas tous le même nombre de caractères). On va ensuite créer des variables qui serviront par la suite : le nombre de vote total, le nombre de vote POUR, CONTRE, ABSTENTION et le nom de l’organe (groupe auquel appartient le deputé à l’assemblée).

loc=paste(getwd(),"/",dest,sep="" )
loc=paste(getwd(),"/",dest,sep="" )
for (i in 1:248){
j=which(colnames(t2) == test)
#Je ne récupère que les personnes dont le i ème mandat est un mandat de député i.e siegeant à l'Assemblée Nationale
#Nom des différentes variables en concatenant les chaines
# Récupération du numero de la colonne du data frame dont le nom correspond a chacune des variable voulues
k=which(colnames(t3) == circo)
l=which(colnames(t3) == dept)
m=which(colnames(t3) == acteur)
n=which(colnames(t3) == mandat)
#Je cree un data frame correspondant "au numero de mandat"
#Et j'ajoute l'ensemble des députés

Ensuite, pour récupérer les votes de chacun des députés à chaque scrutin, on récupére le fichier XML du site data.assemblee-nationale.fr.

loc=paste(getwd(),"/",dest,sep="" )
loc=paste(getwd(),"/",dest,sep="" )

En parcourant la hiérarchie du fichier XML (à l’aide de la fonction xml_children) on obtient les différents fils d’un noeud.


Le premier niveau de noeud correspond à l’ensemble des scrutins. En bouclant sur chacun d’eux, on récupéree dans les fils le résultat du vote de chaque scrutin pour chaque député. On va se limiter aux scrutins dont le mode de publication scrutin (14ème fils) est “Decompte nominatif” : en effet, pour les scrutins par “Decompte Dissident”, il n’y avait que les votes des députés qui n’avaient pas voté dans le sens la majorité, ce qui introduirait un biais ensuite.

Le 16ème fils correspond à la variable ventilation du vote. Ce noeud recense toutes les données concernant le résultat du vote.
Pour chacun des groupes représentés à l’Assemblée présent dans le fils ventilation vote, on récupère le vote de chaque député du groupe dont les modalités sont : “Non Votant”, “Pour”, “Contre” et “Abstention”.
En bouclant sur l’ensemble des scrutins on comptabilise également le nombre de vote réalisé pour chaque député et le nom de son groupe (organeRef).

for (a in liste){
for (b in groupe){
for (i in 1:4){
for (j in df3[,]){
if (i==1){
if (i==1)
if (i==2){
if (i==3){
if (i==4){
res=data.frame(res, tempDataFrame[,2])
res2$circo=ifelse(nchar(res2$circo)==1, str_c("0",res2$circo) , str_c("",res2$circo))
res2$dept=ifelse(nchar(res2$dept)==1, str_c("0",res2$dept) , str_c("",res2$dept))
parIden=aggregate(res2$nbVote, by=list(res2$iden), sum)

Par soucis de visualisation, on va se limiter ici à la France métropolitaine. Il a fallu uniformiser le code IDEN de chaque département pour qu’il corresponde au format utilisé par la carte de France utilisé (DDCC, deux chiffres pour le département, deux pour la circonscription). Le fond de carte est le suivant

loc=paste(getwd(),"/",dest,sep="" )
loc=paste(getwd(),"/",dest,sep="" )

On agrége enfin par par député afin d’obtenir un nombre de vote par circonscription. Notons que dans certaines circonscriptions, le député démissionait ou mourrait et était donc remplacé par un autre député. J’ai donc choisi de sommer le vote de l’ensemble des députés ayant représenté une circonscription pour le représenter graphiquement : on raisonne par circonscription, et par par représentant.

nuancier <- findColours(classIntervals(parIden$x, 6, style = "quantile"), smoothColors("white",98,"#0C3269"))
plot(fdc, col=nuancier)
leg <- findColours(classIntervals( round(parIden$x), 6, style="quantile"), smoothColors("white",98,"#0C3269"), under="moins de", over="plus de", between="–",cutlabels=FALSE)
legend("bottomleft",fill=attr(leg, "palette"), legend=names(attr(leg,"table")),title = "Nombre de Votes",bty="n")
title( main="Nombre de votes par circonscription",cex.main=1.5)

parIdenScrutin=aggregate(res2[,9:ncol(res2)], by=list(res2$iden),na.omit)
for (i in 1:nrow(parIdenScrutin))
for (i in 3:ncol(parIdenScrutin))
for (i in 1:nrow(parIdenScrutin)) {
if (parIdenScrutin$nomOrgane[i]=="PO656014" || parIdenScrutin$nomOrgane[i]=="PO713077" || parIdenScrutin$nomOrgane[i]=="PO656002")
parIdenScrutin$nomOrgane[i] <- "SER"
if (parIdenScrutin$nomOrgane[i] == "PO656006" || parIdenScrutin$nomOrgane[i] =="PO707869")
parIdenScrutin$nomOrgane[i] <- "LR"
if (parIdenScrutin$nomOrgane[i] == "PO656022")
parIdenScrutin$nomOrgane[i] <- "RRDP"
if (parIdenScrutin$nomOrgane[i] == "PO656010")
parIdenScrutin$nomOrgane[i] <- "UDI"
if (parIdenScrutin$nomOrgane[i] == "PO656018")
parIdenScrutin$nomOrgane[i] <- "GDR"
if (parIdenScrutin$nomOrgane[i] == "PO645633")
parIdenScrutin$nomOrgane[i] <- "NI"
parIdenScrutin$Group.1 <- as.factor(parIdenScrutin$Group.1)
parIdenScrutin$nomOrgane <- as.factor(as.character(parIdenScrutin$nomOrgane))

On peut finir par une petite analyse des correspondances : les députés ont voté a 644 scrutins (en colonne). Pour chaque scrutin, un député peut avoir voté pour (POUR), contre (CTRE), s’être abstenu (ABST) ou ne pas s’être présenté. Si on ne garde que les deux premiers axes pricinpaux, on obtient

acm <- MCA(parIdenScrutin, quali.sup=2, graph=FALSE)
plot(acm$ind$coord[, 1:2], type= "n", xlab=paste0("Axe 1 (" , round(acm$eig[1,2], 1), " %)"), ylab=paste0("Axe 2 (", round(acm$eig[2,2], 1), " %) "), main= "Nuage des individus selon les partis", cex.main=1, cex.axis=1, cex.lab=1, font.lab=3)
abline(h=0, v=0, col= "grey", lty=3, lwd=1)
points(acm$ind$coord[,1:2], col=as.numeric(parIdenScrutin$nomOrgane), pch=19, cex=0.5)
legend("topleft", legend=levels(parIdenScrutin$nomOrgane), bty= "o", text.col=1:10, col=1:10, pch=18, cex=0.8)
text(acm$ind$coord[,1:2], labels=rownames(acm$ind$coord), col=as.numeric(parIdenScrutin$nomOrgane), cex=0.7, pos=4)

On a ici le nuage des individus. Ces derniers sont placés en fonction de leur tendance de votes. Les individus les plus proches sont ceux qui votent sensiblement de la même façon tandis que les individus éloignés votent différemment à chaque scrutin. Les couleurs représente le parti auquel appartient le député. Elles sont présentent à titre indicatif pour valider ou non l’hypothèse de “vote en groupes“. Sur ce graphique, parmi les groupes représentés, deux sont très nets à gauche et à droite. Un 3ème se remarque dans la partie centre haute du graphique et enfin un 4e groupe dans la supérieure droite du graphique.

Si on se limite aux couleurs (les partis), les deux groupes très nets qui s’opposent sont Les Républicains (LR – Droite, à gauche sur la projection) et le groupe Socialiste, Ecologiste et Républicain (SER – Gauche, à droite sur la projection). Le 3ème groupe qui se distingue représente le groupe d’extrème Gauche GDR (Gauche démocrate et républicaine). Ces trois groupes sont disctinct car ils votent de façon opposé. Ensuite en regardant toujours les partis, on remarque que l’Union des Démocrates et Indépendants (UDI – Droite) ont des votes proches de ceux du groupe LR. Le groupe Radical, Républicain, démocrate et progressiste ont une façon de voter similaire à SER. Les votes des députés non inscrits varient énormément entre les scrutins, ils ne se rapporchent d’aucun parti.

Accidents de voiture et de deux roues à Paris

Il y a quelques mois, j’avais commencé à reprendre des codes que j’avais vu passer dans les projets de R que j’avais donnés pour la formation Data Science pour l’Actuariat. J’avais commencé avec un billet sur sur cartographier le vote pour le Brexit (qui avait été repris – et amélioré – sur le site des voisins de rgeomatic). Aujourd’hui, en rangeant un peu mon ordinateur, je suis tombé sur un code de Quentin sur les acccidents de la route à Paris.

Comme toujours on charge quelques packages

options(encoding = "UTF-8")

et on récupère les données, en ligne sur https://opendata.paris.fr

loc=paste(getwd(),"/",dest,sep="" )
obj <- read.csv(file = loc, stringsAsFactor = FALSE,sep=";")

On peut commencer par regarder quand on eu lieu les accidents, quelle jour et à quelle heure. On va alors agréger les données

parJour=aggregate(t$counter, by=list(format(as.Date(t$Date),"%A")), sum)

et créer une base qui contient toutes les informations

vec_dat[ (1:3),1]=parJour$Group.1[3:5]
vec_dat[ (1:3),2]=parJour$x[3:5]
vec_dat[ 4,1]=parJour$Group.1[2]
vec_dat[ 4,2]=parJour$x[2]
vec_dat[ 5,1]=parJour$Group.1[7]
vec_dat[ 5,2]=parJour$x[7]
vec_dat[ 6,1]=parJour$Group.1[6]
vec_dat[ 6,2]=parJour$x[6]
vec_dat[ 7,1]=parJour$Group.1[1]
vec_dat[ 7,2]=parJour$x[1]

on va rajouter un nom, pour que ce soit plus clair


On choisit éventuellement une palette de couleurs


et on fait un barplot, pour compter le nombre d’accidents, par jour de la semaine, sur Paris

barplot(vec_dat$Nbsinistre,names.arg=vec_dat$Jour,main="Répartition du nombre d'accidents par jour de la semaine ", col=Palette(10), xlab="jour de la semaine", ylab="nombre d'accidents")

On peut ensuite regarder par heure de la journée (en données brutes, sans lissage)

parHeure=aggregate(t$counter, by=list(format(strptime(t$Heure,format="%H:%M"),"%H")), sum)
plot(parHeure,type="l", col=c("blue"), main=" Répartition du nombre d'accidents par heure de la journée", xlab="heure de la journee", ylab="nombre d'accidents")

Passons maintenant à de la géolocalisation

t2 = subset(t, Coordonnées!= "")
t2 = subset(t2, (
(48.81 <= lat) & (lat<=48.90) & (2.25<=lon) &(lon<=2.43) ) )

On peut regarder sur une carte où ont eu lieu les accidents

map<-get_map(location = "Paris", zoom=12, maptype="roadmap", color="bw")
vis <- ggmap(map) + geom_point(data = t2, aes(x=lon, y=lat, color=1)) + theme( legend.position="none") + labs(title="Cartographie des accidents à Paris")

On peut tenter une densité (brute, non corrigée d’une densité de routes, ou de population)

overlay<-stat_density2d(data=t2, aes(x=lon, y=lat, fill=..level.., alpha=..level..), contour=T, n=100, geom="polygon")
densi<- ggmap(map) + overlay + scale_fill_gradient("Accident Density") +scale_alpha(range=c(0.4,0.75),guide=FALSE) + guides(fill=guide_colorbar(barwidth=1.5, barheight = 10))
densi = densi + labs(title="Densite des accidents à Paris")

et non-corrigée ici d’effets de bords (ce qui donne l’illusion qu’il y a plus d’accidents intra-muros que sur le boulevard périphérique par exemple). On peut également regarder en fonction de variables explicatives, comme le type de véhicules impliqués, ou la gravité

densi_byGravUs<- densi + facet_wrap(~ Usager.1.Grav, nrow=2)
densi_byGravUs = densi_byGravUs + labs(title="Densité des accidents à Paris par gravité d'accident")

Bref, on peut faire facilement plein de chose à partir de ces données.

Visualiser les résultats du premier tour

Juste quelques lignes de code, pour visualiser les résultats du premier tour en France. L’idée est de faire une carte assez minimaliste, avec des cercles centrés sur les centroïdes des départements. On commence par récupérer les données pour le fond de carte, un fichier 7z sur le site de l’ign.

download.file("https://wxs-telechargement.ign.fr/oikr5jryiph0iwhw36053ptm/telechargement/inspire/GEOFLA_THEME-DEPARTEMENTS_2016$GEOFLA_2-2_DEPARTEMENT_SHP_LAMB93_FXX_2016-06-28/file/GEOFLA_2-2_DEPARTEMENT_SHP_LAMB93_FXX_2016-06-28.7z",destfile = "dpt.7z")

On a dans ce fichier des informations sur les centroïdes


Comme ça ne marche pas très bien, on va le refaire à la main, par exemple pour l’Ille-et-Vilaine,


Comme ça marche mieux, on va utiliser ces centroïdes.


Maintenant, il nous faut les résultats des élections, par département. On peut aller scraper le site du ministère d’ intérieur. On a une page par département, donc c’est facile de parcourir. Par contre, dans l’adresse url, il faut le code région. Comme je suis un peu fainéant, au lieu de faire une table de correspondance, on teste tout, jusqu’à ce que ça marche. L’idée est d’aller cherche le nombre de voix obtenues par un des candidats.

candidat="M. Emmanuel MACRON"
# if(!is.na(as.numeric(no))){if(as.numeric(no)<10) nodpt=paste("00",no,sep="")}
test=try(htmlParse(url),silent =TRUE)
if(!inherits(test, "try-error")){testurl=TRUE
tabs tab=tabs[[2]]
a<-unlist(strsplit(as.character(nb)," "))

On peut alors tester

> voix(35)
[1] 84648

Comme ça semble marcher, on le fait pour tous les départements


On peut alors visualiser sur une carte.


Et on peut aussi tenter pour une autre candidate,

candidat="Mme Marine LE PEN"

et on obtient la carte suivante


Visualizing (censored) lifetime distributions

There are now more than 10,000 R packages available from CRAN, much more if you include those available only on github. So, to be honest, it become difficult to know all of them. But sometimes, you discover a nice function in one of them, and that is really awesome. Consider for instance some (standard) censored lifetime data,


(yes, I will generate them here). Consider Kaplan-Meier estimator of the survival function,

km.base = survfit( Surv(time,status) ~ 1  , data = df )

This week end, Anat (currently finishing the Data Science for Actuaries program) made me discover a nice R function, to add information to that graph (well, not that graph, since it will be a ggplot version, but the same survival distribution plot)

ggsurvplot(km.base, main = "", color = "blue" , censor = FALSE, xlim = c(0,3), risk.table = TRUE ,
risk.table.col = "blue" , risk.table.height = 0.2, risk.table.title = "" , legend.labs = "All" , legend.title = "" , break.time.by = 1, xlab = "" , ylab = "")

This is more interesting when we have different lifetimes

km.prod = survfit( Surv(time,status) ~ product  , data = df )
ggsurvplot(km.prod, main = "", censor = FALSE, xlim = c(0,3), risk.table = TRUE , risk.table.col = "strata" , risk.table.height = 0.3, risk.table.title = "" , legend.labs = LETTERS[1:4] , legend.title = "" , break.time.by = 1, xlab = "" , ylab = "")

or, with a different time granularity

ggsurvplot(km.prod, main = "", censor = FALSE, xlim = c(0,3), risk.table = TRUE , risk.table.col = "strata" , risk.table.height = 0.3, risk.table.title = "" , legend.labs = LETTERS[1:4] , legend.title = "" , break.time.by = .5, xlab = "" , ylab = "")

Nice, isn’t it?

The U.S. Has Been At War 222 Out of 239 Years

This morning, I discovered an interesting statistic, America Has Been At War 93% of the Time – 222 Out of 239 Years – Since 1776,  i.e. the U.S. has only been at peace for less than 20 years total since its birth. I wanted to check, get a better understanding and look at other countries in the world.

As always, we can try to extract information from wikipedia, since there are pages dedicated to that information

download.file(url,destfile = "warUS.html")
download.file(url,destfile = "warFR.html")
download.file(url,destfile = "guerre.html")
download.file(url,destfile = "warCAN.html")

If we look at the US page, there are tables, so it should be easy to extract it. For instance,

Even if the war did last 1 day, we will say that the US were at war in 1811. The information we want to confirm can be “there were 21 full years – from Jan 1st till Dec 31st – where the US were not at war, once, during those years“. From the row above, we can claim that the US were at war in 1811. Most of the time, we have

I.e. there is a beginning (here 1775) and an end (1783). So here, the US are said to be at war in 1775, 1776, 1777, 1778, 1779, 1780, 1781, 1782, 1783. To extract the information, we look for regular expressions in the first column, with number, on 4 digits.


Well, sometimes it can be a bit tricky, since we have 3 dates, 1941, 1945 and (in the legend) 1944. But if we consider the minimal and the maximal dates, we have our range of dates.

Now that we we how to extract information, let’s do it. The code will be

#grep(pattern = dates2, x = col1[1])
for(j in 1:length(L))
if(length(L[[j]])==1) return_L[[j]]=as.numeric(L[[j]])
if(length(L[[j]])>=2) return_L[[j]]=seq(min(as.numeric(L[[j]])),max((as.numeric(L[[j]]))))

For the US, we get the following years

for(i in 1:length(tables)){

(red means at war, while green means no-war) and indeed,

> length(d)
[1] 222

there were 222 years with war.  Now, what about another country. Like France. Here I use the French wiki page, since information is not in tables in the English one.

for(i in 1:length(tables)){

On the same period of time (starting in 1775), France was also on war most of the time.

Less than the US, but still: 185 years with war,

> length(d[d>=1775])
[1] 185

And on a longer period of time? Why not start, say, around the Hundred Years’s War,

meaning that since 1337, there were (only) 174 years without a single war where France was involved.

Let’s try another one. Like Canada,

for(i in 1:length(tables)){

Guess what… there’s a lot of green on that graph. Surprised?

Reading text automatically

It is now very easy to read (automatically) some text that can be found in a pdf file. For instance, consider the program of the conference we had yesterday – and today – in Rennes

> library(pdftools)
> scan_pdf <- pdf_text("http://crem.univ-rennes1.fr/Documents/Docs_sem_divers/2017_03_10-11_JJD/JDD_prog.pdf")
> cat(scan_pdf)
Journées Jeunes Docteurs
Programme du jeudi 9 mars 2017
Faculty of Economics - Rennes - Amphi Henri Krier
9h- 9h30 - Accueil
9h30-10h15 :      Présentation du CREM, de la faculté et des activités de recherche liées du ou laboratoire
10h15-10h50 :     Emmanuel LORENZON (Université de Bordeaux, GREThA)
Collusion with a rent seeking agency in sponsored search auctions
10h50-11h25 :     Julien BERTHOUMIEU (Université de Bordeaux, GREThA)
The Impact of “At-the-Border” and “Behind-the-Border” Policies on Cost-Reducing Research
and Development
Co-écrit avec Antoine Bouët

(etc). As you can see, it is working well, even in French, where we have those weird letters (with accents). Here, it is working well because the pdf is vectorized, i.e. it was generated properly, by open office.

But sometimes, we can have only a scanned version of a letter

or just a picture with some typed text. I will not mention hand-writing because it is much more complex.

The other day, my friend Fleur did show me a picture, and some very simple lines of code,

> library('tesseract')
> pic1="https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2017/03/pic1.png"
> text_fr <- ocr(pic1, engine = tesseract("fra"))
> cat(text_fr)
Près de 14.400 décès

Si [épidémie de grippe est un phénomène récurrent. celle de
2016—2017 présente plusieurs spécificités. outre sa virulence :
une survenue plus précoce que d‘habitude. une activité
modérée en médecine ambulatoire. mais un impact fort en
milieu hospitalier.

It looks like we’ve be able to extract typed text from a picture ! I want to check. I have to admit, first of all that installation on a linux machine is tricky: one has to install first leptonica, and then follow some guidelines to install tesseract (see also Artem‘s advices). It took me some time, but I’ve been able to install the package.

The first important step, it to train the algorithm with some texts in French (because it is in French in my picture)

> library('tesseract')
> tesseract_download("fra")

Then, I did try with the picture that Fleur did send me (the picture was inserted in the core of the message)

> pic2="https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2017/03/pic2.png"
> text_fr <- ocr(pic2, engine = tesseract("fra"))
> cat(text_fr)
Près de 14400 decès

s. mm….agw«… ………«…m ……
a……u…u Dhs—ur; ;pmum…. ;: u……
… W»: »… w…q…na… … ……
…… ………u…_ …… mm……
…… nwm/u

… … mm…—mg…»— sa…… su a.…….…
: :mmræwesdæ ; ; m…decnflwtülflws WWW…
un»… M on m…… . … … .. m...… wma:
.,… … V, …… … …;………yg…gn…
…… pe- le…samemeuuœwpwv m…

Clearly, something went wrong here. When I got that output, I thought that I did not train properly the function. But it was not the answer. As described in that post (in French) it is necessary to have a clean picture, to read it properly

And actually, if we zoom in our picture – the first one, used by Fleur, to show me that package – we have

while for the second one – with a lower resolution – we have

It is necessary to have a scan of a typed text with high resolution… And you have to admit that it is awesome….

The good thing is that I have to work with a judge, in France, to assess quality of experts. And since most of the reports are typed, and then scanned, I am glad to have such a function. I just have to make sure that the resolution is high enough…

Juxtaposer des cartes et ‘confounding’

Pour aller un peu plus loin que mon précédant billet, élections, pièges à com(mentaires à la con) je voulais juste montrer quelques cartes, obtenues par simulations. Mais avant, revenons deux minutes sur la notion de ‘confounding‘. Pour reprendre le dessin de wikipedia, l’idée est que, conditionnellement à Z (notre ‘confounding factor’), les variables X et Y sont indépendantes. Mais cette indépendance conditionnelle peut (comme bien souvent) impliquer une forte corrélation entre X et Y. On peut ainsi imaginer un modèle de la forme X=Z+\varepsilon \varepsilon désigne une composante idiosyncratique, et Y=Z+\eta (où là aussi \eta est une composante idiosyncratique).

Par exemple, on peut supposer que Z est une variable liés à la lattitude du département, par exemple la température moyenne (histoire d’avoir une vague histoire à raconter, et d’être moins rébarbatif)

> n=115
> X=matrix(rnorm(3*n),n,3)
> X[,1]=sort(X[,1])
> library(maps)
> map("france",fill =TRUE,
col= rgb(1,0,0,pnorm(X[,1])))

Maintenant, supposons qu’on ait une variable X très corrélée à Z, par exemple la vente de crème solaire (en volume, toujours pour raconter une histoire qui illustrerait un peu mon propos)

> r=.8
> X[,2]=r*X[,1]+sqrt(1-r^2)*X[,2]
> map("france",fill =TRUE,
col= rgb(0,1,0,pnorm(X[,1])))

ce qui donne le nuage de points suivant (avec une forte corrélation avec le ‘confounding factor’)

On suppose maintenant qu’on a une autre variable Y , elle aussi très corrélée à Z, par exemple le prix d’une plaquette de beurre (qui pourrait être très cher dans le sud, où l’on cuisine plutôt à l’huile)

> r=.8
> X[,3]=r*X[,1]+sqrt(1-r^2)*X[,3]
> map("france",fill =TRUE,
col= rgb(0,1,0,pnorm(X[,1])))

avec une corrélation du même ordre, avec le ‘confounding factor’

Si on juxtapose les deux cartes, on a des variables corrélées,

mais on ne peut pas en conclure, juste en regardant ces cartes, que les ventes de crème solaire influencent le prix d’une plaquette de beurre ! Et on peut même accentuer encore l’effet, en jouant un peu sur les couleurs. Par exemple, avec

> map("france",fill =TRUE,col = rgb(1,0,0,qbeta(pnorm(X[,1]),.2,.5)))

on accentue un peu l’écart entre le nord et le sud (je joue ici sur l’échelle des couleurs, pas sur les valeurs !)

et sur nos deux variables, on a encore une juxtaposition qui donne envie de raconter un jolie histoire, non ?

Ah oui, ma corrélation entre les variables X et Y est relativement forte, de l’ordre de 0.6 (qui est supérieure à ce qu’on avait observé hier sur le taux de chômage et le taux du FN aux élections régionales)

> cor(X)
[,1]      [,2]      [,3]
[1,] 1.0000000 0.7863838 0.7727523
[2,] 0.7863838 1.0000000 0.6312878
[3,] 0.7727523 0.6312878 1.0000000

mais la corrélation est très sensible aux lois marginales, ça ne veut pas dire grand chose. Même si, visuellement, le nuage de points peut faire penser à celui qu’on avait hier

Bref, juxtaposer des cartes, c’est joli. De là à en déduire quoi que ce soit….