# 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

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

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

generate3=function(n,x,pb=c(-2,0)){ set.seed(x) X1=runif(n) X2=runif(n) X3=runif(n) s1=pb[1]+X1+X2 s2=pb[2]-X1+X2 P1=exp(s1)/(1+exp(s1)+exp(s2)) P2=exp(s2)/(1+exp(s1)+exp(s2)) Y=msample(0:2,n,cbind(1-P1-P2,P1,P2)) df=data.frame(Y=Y,X1=X1,X2=X2,X3=X3) return(df) }

Let us generate a training dataset and a validation one

 pb=c(.31,.42) DF1=generate3(1000,1,pb=pb) DF2=generate3(500,2,pb=pb)

With a multivariate logistic regression
$$\mathbb{P}[Y=A|\mathbf{x}]=\frac{\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]+\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}$$
$$\mathbb{P}[Y=B|\mathbf{x}]=\frac{\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]+\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}$$
$$\mathbb{P}[Y=B|\mathbf{x}]=\frac{1}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{\alpha}]+\exp[\mathbf{x}^{\text{T}}\mathbf{\beta}]}$$

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.  library(nnet) 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$. $$\mathbb{P}[Y_A=A|\mathbf{x}]=\frac{\exp[\mathbf{x}^{\text{T}}\mathbf{a}]}{1+\exp[\mathbf{x}^{\text{T}}\mathbf{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]$.  reg1=glm((Y==modalite[1])~.,data=DF1,family=binomial) reg2=glm((Y==modalite[2])~.,data=DF1[-which(DF1$Y==modalite[1]),],family=binomial) 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") mmp1=cbind(p11,(1-p11)*p12,(1-p11)*(1-p12)) mmp2=cbind(p21,(1-p21)*p22,(1-p21)*(1-p22)) colnames(mmp1)=colnames(mmp2)=modalite

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.

# Traveling Salesman

In the second part of the course on graphs and networks, we will focus on economic applications, and flows. The first series of slides are on the traveling salesman problem. Slides are available online.

# 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 = + FALSE)

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 = + FALSE) 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)$track-table 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

 BASE=NULL > 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

BASE=BASE[-duplicated(BASE),] BASE=BASE[-which(BASE$Year=="—"),] BASE$y=as.numeric(as.character(BASE$Year)) BASE$bpm=as.numeric(as.character(BASE$BPM)) BASE=BASE[BASE$y>=1940,]

and we end up with almost 50,000 tunes.

boxplot(BASE$bpm~as.factor(BASE$y), 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…)

library(mgcv) plot(BASE$y,BASE$bpm) reg=gam(bpm~s(y),data=BASE) B=data.frame(y=1950:2017) p=predict(reg,newdata=B) lines(B$y,p,lwd=3,col="red") which is confirmed with a (smoothed) regression p2=predict(reg,newdata=B,se.fit=TRUE) plot(B$y,p2$fit,lwd=3,col="red",type="l",ylim=c(90,140)) lines(B$y,p2$fit+p2$se.fit) lines(B$y,p2$fit-p2$se.fit) 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, require(xml2) require(downloader) require(stringr) require(classInt) require(plotrix) require(FactoMineR) require(sp) 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). path="http://data.assemblee-nationale.fr/static/openData/repository/AMO/deputes_senateurs_ministres_legislature/AMO20_dep_sen_min_tous_mandats_et_organes_XIV.csv.zip" dest="deputes.zip" download(path,destfile=dest,mode="wb") loc=paste(getwd(),"/",dest,sep="" ) unzip(loc) dest="acteurs.csv" loc=paste(getwd(),"/",dest,sep="" ) t2=read.csv(loc,sep=";") res=NULL for (i in 1:248){ test=str_c("mandats.1..mandat.",i,"..typeOrgane.1.") 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 t3=subset(t2,t2[,j]=="ASSEMBLEE") if(nrow(t3)!=0){ #Nom des différentes variables en concatenant les chaines circo=str_c("mandats.1..mandat.",i,"..election.1..lieu.1..numCirco.1.") dept=str_c("mandats.1..mandat.",i,"..election.1..lieu.1..numDepartement.1.") acteur=str_c("mandats.1..mandat.",i,"..acteurRef.1.") mandat=str_c("mandats.1..mandat.",i,"..uid.1.") # 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" t4=data.frame(as.numeric(as.character(t3[,k])),as.numeric(as.character(t3[,l])),str_c(t3[,m],t3[,n])) colnames(t4)=c("circo","dept","identifiant") #Et j'ajoute l'ensemble des députés res=rbind(res,t4) } } res$iden=str_c(res$dept,res$circo) res$nbVote=0 res$oui=0 res$non=0 res$abst=0 res$nomOrgane=0 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. path="http://data.assemblee-nationale.fr/static/openData/repository/LOI/scrutins/Scrutins_XIV.xml.zip" dest="Scrutins_XIV.xml.zip" download(path,destfile=dest,mode="wb") loc=paste(getwd(),"/",dest,sep="" ) unzip(loc) dest="Scrutins_XIV.xml" 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. t=read_xml(loc) liste=xml_children(t) 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){ if(xml_text(xml_children(a)[14])=="DecompteNominatif"){ numero=xml_text(xml_children(a)[2]) ventil=xml_children(a)[16] groupe=xml_children(xml_children(xml_children(ventil))) tempDataFrame=data.frame(res$identifiant,NA) colnames(tempDataFrame)=c("identifiant",str_c("scrutin",numero)) for (b in groupe){ intermediaire=xml_children(xml_children(b)) nomGroupe=xml_text(xml_children(b)[1]) for (i in 1:4){ df3=data.frame(xml_text(xml_children(xml_children(intermediaire[3])[i]))) if(nrow(df3)!=0){ for (j in df3[,]){ if (i==1){ j=strsplit(strsplit(strsplit(j,"MG")[[1]],"PSE")[[1]],"PAN")[[1]] } res[res$identifiant==j,]$nomOrgane=as.character(nomGroupe) if(i!=1) res[res$identifiant==j,]$nbVote=res[res$identifiant==j,]$nbVote+1 if (i==1) tempDataFrame[tempDataFrame$identifiant==j,2]="NV" if (i==2){ res[res$identifiant==j,]$oui=res[res$identifiant==j,]$oui+1 tempDataFrame[tempDataFrame$identifiant==j,2]="POUR" } if (i==3){ res[res$identifiant==j,]$non=res[res$identifiant==j,]$non+1 tempDataFrame[tempDataFrame$identifiant==j,2]="CTRE" } if (i==4){ res[res$identifiant==j,]$abst=res[res$identifiant==j,]$abst+1 tempDataFrame[tempDataFrame$identifiant==j,2]="ABST" }}}}} res=data.frame(res, tempDataFrame[,2]) }} res2=subset(res,res$dept<96) 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)) res2$iden=str_c(res2$dept,res2$circo) res2=subset(res2,!is.na(iden)) res2=res2[order(res2$iden),] 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 path="http://www.laspic.eu/data/documents/circosshp_v3.zip" dest="circosshp_v3.zip" download(path,destfile=dest,mode="wb") loc=paste(getwd(),"/",dest,sep="" ) unzip(loc) dest="circosSHP_v3.RData" loc=paste(getwd(),"/",dest,sep="" ) load(loc) 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)) parIdenScrutin$nomOrgane[i]=ifelse(parIdenScrutin$nomOrgane[i][[1]][1]!="0",as.character(parIdenScrutin$nomOrgane[i][[1]][1]),as.character(parIdenScrutin$nomOrgane[i][[1]][2])) for (i in 3:ncol(parIdenScrutin)) parIdenScrutin[,i]=as.factor(levels(parIdenScrutin[,i][[1]])[as.numeric(parIdenScrutin[,i])]) parIdenScrutin=replace(parIdenScrutin,is.na(parIdenScrutin),"NV") 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) head(acm$eig,5) 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 require(downloader) require(ggmap) require(stringr) options(encoding = "UTF-8") et on récupère les données, en ligne sur https://opendata.paris.fr path="https://opendata.paris.fr/explore/dataset/accidentologie/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true" dest="test" download(path,destfile=dest,mode="wb") loc=paste(getwd(),"/",dest,sep="" ) obj <- read.csv(file = loc, stringsAsFactor = FALSE,sep=";") file.remove(dest) On peut commencer par regarder quand on eu lieu les accidents, quelle jour et à quelle heure. On va alors agréger les données t=data.frame(obj,counter=1) 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=data.frame(); 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 names(vec_dat)[1]="Jour";names(vec_dat)[2]="Nbsinistre" On choisit éventuellement une palette de couleurs Palette<-colorRampPalette(c("lightblue","darkblue")) 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!= "") coord=data.frame(t(matrix(as.numeric(unlist(strsplit(t2$Coordonnées,"[,]"))),nrow=2))) t2$lat=coord[,1] t2$lon=coord[,2] 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") vis

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

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

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  library(maptools) library(maps) departements<-readShapeSpatial("DEPARTEMENT.SHP") plot(departements) points(departements@data$X_CENTROID,departements@data$Y_CENTROID,pch=19,col="red") Comme ça ne marche pas très bien, on va le refaire à la main, par exemple pour l’Ille-et-Vilaine, pos=which(departements@data[,"CODE_DEPT"]==35) Poly_35=departements[pos,] plot(departements) plot(Poly_35,col="yellow",add=TRUE) departements@data[pos,c("X_CENTROID","Y_CENTROID")] points(departements@data[pos,c("X_CENTROID","Y_CENTROID")],pch=19,col="red") library(rgeos) (ctd=gCentroid(Poly_35,byid=TRUE)) points(ctd,pch=19,col="blue") Comme ça marche mieux, on va utiliser ces centroïdes.  ctd=as.data.frame(gCentroid(departements,byid=TRUE)) plot(departements) points(ctd,pch=19,col="blue") 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" library(XML) voix=function(no=35){ testurl=FALSE i=1 vect_reg=c("084","027","053","024","094","044","032","028","075","076","052", "093","011","001","002","003","004") region=NA while((testurl==FALSE)&(i<=20)){ reg=vect_reg[i] nodpt=paste("0",no,sep="") # if(!is.na(as.numeric(no))){if(as.numeric(no)<10) nodpt=paste("00",no,sep="")} url=paste("http://elections.interieur.gouv.fr/presidentielle-2017/",reg,"/",nodpt,"/index.html",sep="") test=try(htmlParse(url),silent =TRUE) if(!inherits(test, "try-error")){testurl=TRUE region=reg} i=i+1 } tabs tab=tabs[[2]] nb=tab[tab[,1]==candidat,"Voix"] a<-unlist(strsplit(as.character(nb)," ")) as.numeric(paste(a,collapse=""))} On peut alors tester > voix(35) [1] 84648 Comme ça semble marcher, on le fait pour tous les départements  liste_dpt=departements@data$CODE_DEPT nbvoix=Vectorize(voix)(liste_dpt)

On peut alors visualiser sur une carte.

 plot(departements,border=NA) points(ctd,pch=19,col=rgb(1,0,0,.25),cex=nbvoix/50000)

Et on peut aussi tenter pour une autre candidate,

 candidat="Mme Marine LE PEN"

et on obtient la carte suivante

 plot(departements,border=NA) points(ctd,pch=19,col=rgb(0,0,1,.25),cex=nbvoix/50000)

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,

n=10000 idx=sample(1:4,size=n,replace=TRUE) pd=LETTERS[idx] lambda=1+(idx-1)/3 t=rexp(n,lambda) x=rexp(n) c=t>x y=pmin(t,x) df=data.frame(time=y,status=c,product=pd)

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

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

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)

library(ggplot2) library(survminer) 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

url="https://en.wikipedia.org/wiki/List_of_wars_involving_the_United_States" download.file(url,destfile = "warUS.html") url="https://en.wikipedia.org/wiki/List_of_wars_involving_France" download.file(url,destfile = "warFR.html") url="https://fr.wikipedia.org/wiki/Liste_des_guerres_de_la_France#Premi.C3.A8re_R.C3.A9publique" download.file(url,destfile = "guerre.html") url="https://en.wikipedia.org/wiki/List_of_wars_involving_Canada" 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
 library(stringr) ext_date=function(x){ dates12="[0-9]{4}" #grep(pattern = dates2, x = col1[1]) L=str_extract_all(as.character(x),dates12) return_L=list() if(length(L)>0){ 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]])))) } return(return_L)}

For the US, we get the following years

library(XML) tables=readHTMLTable("warUS.html") list_dates=list() for(i in 1:length(tables)){ if(!is.null(dim(tables[[i]]))){ if(ncol(tables[[i]])>1){ col1=tables[[i]][,1] list_dates[[i]]=lapply(col1,ext_date) } }} d=unique(unlist(list_dates))

(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.

tables=readHTMLTable("guerre.html") list_dates=list() for(i in 1:length(tables)){ if(!is.null(dim(tables[[i]]))){ if(ncol(tables[[i]])>1){ col1=tables[[i]][,1] col2=tables[[i]][,2] col12=paste(col1,col2) list_dates[[i]]=lapply(col12,ext_date) } }} d=unique(unlist(list_dates))

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,

tables=readHTMLTable("warCAN.html") list_dates=list() for(i in 1:length(tables)){ if(!is.null(dim(tables[[i]]))){ if(ncol(tables[[i]])>1){ col1=tables[[i]][,1] list_dates[[i]]=lapply(col1,ext_date) } }} d=unique(unlist(list_dates))

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

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éciﬁcité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…decnﬂwtülﬂws WWW… un»… M on m…… . … … .. m...… wma: .,… … V, …… … …;………yg…gn… …… pe- le…samemeuuœwpwv m… mum

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….

# (Advanced) R Crash Course, for Actuaries

In two weeks, the third year of the Data Science for Actuaries program will start. I will be there for the introduction to R. The slides are available online (created with slidify)

A markdown is also available. I have to thank Ewen for his help on slidify (especially for the online quizz, and the integration of leaflet maps or the rgl animated graph….)

# Evolution du nombre de demandeurs d’emploi

Allez, on continue avec les projets rendus dans le cadre du projet R de la formation en Data Science pour l’Actuariat, Meryem Schalck a proposé une jolie animation des statistiques des demandeurs d’emploi. Sur http://dares.travail-emploi.gouv.fr/, on récupère les séries d’évolution du nombre de demandeurs d’emploi, par département (entre autres), sous forme de fichier Excel. Histoire de simplifier, je mets en ligne une copie en csv,

ufc<- read.csv("http://freakonometrics.free.fr/DemandeV2.csv", header=TRUE, sep = ";")

Pour charger les packages, Meryem propose une jolie petite fonction

Chargepackages <- function(x){ x <- as.character(x) if(!require(x, character.only = TRUE)){ install.packages(x, dependencies = TRUE) require(x, character.only = TRUE) } } packages <- list("sp", "rgdal", "animation", "rgeos", "mapproj", "maptools", "dplyr", "ggplot2", "RColorBrewer", "classInt", "PBSmapping", "ggmap", "splancs", "osmar", "maps", "scales", "osmar", "geosphere", "RCurl", "bitops", "XML", "splancs", "PBSmapping", "htmltools" ) lapply(packages, Chargepackages)

On va retravailler un peu les données pour commencer

x<-ufc$nombre y<-ufc$demande dt<-as.Date(ufc$date,"%d/%m/%Y") n<-length(x) et ensuite on fait notre animation, colfunc <- colorRampPalette(c("green","yellow", "red")) YlOrBr <- data.frame(col=colfunc(length(unique(y))),num=unique(y)[order(unique(y))]) ufc %>% left_join(YlOrBr, by=c("demande"="num")) -> YlOrBr oopt = ani.options(interval = 0.5) library(scales) YlOrBr$Date<-as.Date(YlOrBr$date,"%d/%m/%Y") ani.record(reset = TRUE) # clear history before recording for (i in 1:length(unique(dt))) { YlOrBr[i>=x,] %>% ggplot() + geom_point(aes(x=Date,y=demande),colour=as.character(YlOrBr$col[i>=x]))+ scale_x_date(breaks=date_breaks("6 month"),limits = c(min(dt),max(dt))) + ylim(c(1000,max(y))) + ggtitle ("Nombre de demandeurs d'emploi par mois entre 01/2000 et 04/2016") + xlab("Date") + ylab("Nombre de demandeurs d'emploi cat A ensemble") + theme(axis.text.x  = element_text(angle=45, vjust=0.5, size=5)) -> p print(p) ani.record() } oopts = ani.options(interval = 0.2) ani.record() saveHTML(ani.replay(),img.name = "record_plot")

et voilà

# Cartographier le vote pour le Brexit

Je suis, ces temps-ci plongé dans les projets de R que j’avais donnés pour la formation Data Science pour l’Actuariat. Comme j’ai eu plein de choses intéressantes, je me suis dit que je pourrais tenter de reprendre dans des billets sur le blog (en plus, ça me permet de vérifier au passage que les codes tournent).

Pour commencer la série (oui, il y en aura d’autres), Flavien Thery a proposé une cartographie du vote pour le Brexit. La première étape était de récupérer un fond de carte (on va utiliser la seconde zone administrative ici)

library(sp) library(raster) download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/GBR_adm2.rds","GBR_adm2.rds") UK=readRDS("GBR_adm2.rds") UK@data[159,"HASC_2"]="GB.NR" plot(UK, xlim = c(-4,-2), ylim = c(50, 59), main="UK areas")

Cette carte est un peu surprenante… parce qu’on est sur le Royaume uni, ce qui inclu l’Angleterre, le Pays de Galles, l’Ecosse et l’Irlande du Nord. Et l’Irlande du Sud n’est pas sur la carte… Autant la rajouter (ça sera plus joli si on met du bleu pour l’eau)

download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/IRL_adm0.rds","IRL_adm0.rds") IRL=readRDS("IRL_adm0.rds") plot(IRL,add=TRUE)

Mais on commence à rajouter l’Irlande, pourquoi pas la France, qui est juste en bas à droite, et qu’on devrait voir un peu…

download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds","FRA_adm0.rds") FR=readRDS("FRA_adm0.rds") plot(FR,add=TRUE)

Bon, arrêtons là….

On peut ensuite récupérer les données du référendum sur le Brexit (que je stocke sur le blog histoire de gagner un peu de temps)

loc="https://f-origin.hypotheses.org/wp-content/blogs.dir/253/files/2016/12/EU-referendum-result-data.csv" referendum=read.csv(loc,header=TRUE,dec=".",sep=",",stringsAsFactors = FALSE) referendum=referendum[c(3,6,13,14)] library(plyr) referendum=ddply(referendum,.(Region,HASC_code),summarise,Remain=sum(Remain),Leave=sum(Leave))

On peut retrouver que le Brexit a gagné, avec 51,89% des votes exprimés (ce qui est conforme à ce que dit wikipedia)

> sum(referendum$Leave)/(sum(referendum$Leave)+sum(referendum$Remain)) [1] 0.5189184 On peut regarder, région par région, si c’est le leave ou le remain qui a gagné, en utilisant referendum=referendum[c(referendum$Region!="Northern Ireland"),] referendum=referendum[c(referendum$HASC_code!="Gibraltar"),] row.names(referendum)=seq(1,nrow(referendum),1) leave_or_remain=cbind(referendum,"Brexit?"=0) leave_or_remain[,"Brexit?"]=ifelse(leave_or_remain$Remain<leave_or_remain$Leave,rgb(1,0,0,.7),rgb(0,0,1,.7)) map_data=data.frame(UK@data) map_data=cbind(map_data,"Brexit"=0) for (i in 1:nrow(map_data)){ if(map_data[i,"NAME_1"]=="Northern Ireland"){ map_data[i,"Brexit"]="blue"}else{ map_data[i,"Brexit"]=as.character(leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],'Brexit?']) } } plot(UK, col = map_data$Brexit, border = "gray1", xlim = c(-4,-2), ylim = c(50, 59), main="How the UK has voted?", bg="#A6CAE0") plot(IRL, col = "lightgrey", border = "gray1",add=TRUE) plot(FR, col = "lightgrey", border = "gray1",add=TRUE)

legend(-1,59,c("Leave","Remain"),fill=c("red","blue"),bty="n")

(on rajoute une petite légende pour être plus clair). Mais on peut aller plus loin en représentant le pourcentage obtenu, par région, par le camp du leave. Pour ça, on peut utiliser le package cartography

library(cartography) cols <- carto.pal(pal1 = "red.pal",n1 = 5,pal2="green.pal",n2=5) map_data=data.frame(UK@data) map_data=cbind(map_data,"Percentage_Remain"=0) for (i in 1:nrow(map_data)){ if(map_data[i,"NAME_1"]=="Northern Ireland"){ map_data[i,"Percentage_Remain"]=55.78}else{         map_data[i,"Percentage_Remain"]=100*leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],"Remain"]/(leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],"Remain"]+leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],"Leave"]) } } plot(UK, col = "grey", border = "gray1", xlim = c(-4,-2), ylim = c(50, 59),bg="#A6CAE0") plot(IRL, col = "lightgrey", border = "gray1",add=TRUE) plot(FR, col = "lightgrey", border = "gray1",add=TRUE) choroLayer(spdf = UK, df = map_data, var = "Percentage_Remain", breaks = seq(0,100,10), col = cols, legend.pos = "topright", legend.title.txt = "", legend.values.rnd = 2, add = TRUE)

(là encore, une légende permet d’aider à lire la carte). Amusant, n’est ce pas ?