Solving the chinese postman problem

Some pre-Halloween post today. It started actually while I was in Barcelona : kids wanted to go back to some store we’ve seen the first day, in the gothic part, and I could not remember where it was. And I said to myself that would be quite long to do all the street of the neighborhood. And I discovered that it was actually an old problem. In 1962, Meigu Guan was interested in a postman delivering mail to a number of streets such that the total distance walked by the postman was as short as possible. How could the postman ensure that the distance walked was a minimum?

A very close notion is the concept of traversable graph, which is one that can be drawn without taking a pen from the paper and without retracing the same edge. In such a case the graph is said to have an Eulerian trail (yes, from Euler’s bridges problem). An Eulerian trail uses all the edges of a graph. For a graph to be Eulerian all the vertices must be of even order.

An algorithm for finding an optimal Chinese postman route is:

  1. List all odd vertices.
  2. List all possible pairings of odd vertices.
  3. For each pairing find the edges that connect the vertices with the minimum weight.
  4. Find the pairings such that the sum of the weights is minimised.
  5. On the original graph add the edges that have been found in Step 4.
  6. The length of an optimal Chinese postman route is the sum of all the edges added to the total found in Step 4.
  7. A route corresponding to this minimum weight can then be easily found.

For the first steps, we can use the codes from Hurley & Oldford’s Eulerian tour algorithms for data visualization and the PairViz package. First, we have to load some R packages

require(igraph)
require(graph)
require(eulerian)
require(GA)

Then use the following function from stackoverflow,

make_eulerian = function(graph){
  info = c("broken" = FALSE, "Added" = 0, "Successfull" = TRUE)
  is.even = function(x){ x %% 2 == 0 }
  search.for.even.neighbor = !is.even(sum(!is.even(degree(graph))))
  for(i in V(graph)){
    set.j = NULL
    uneven.neighbors = !is.even(degree(graph, neighbors(graph,i))) 
if(!is.even(degree(graph,i))){ 
if(sum(uneven.neighbors) == 0){ 
if(sum(!is.even(degree(graph))) > 0){
          info["Broken"] = TRUE
          uneven.candidates <- !is.even(degree(graph, V(graph)))
          if(sum(uneven.candidates) != 0){
            set.j <- V(graph)[uneven.candidates][[1]]
          }else{
            info["Successfull"] <- FALSE
          }
        }       
      }else{
        set.j <- neighbors(graph, i)[uneven.neighbors][[1]]
      }
    }else if(search.for.even.neighbor == TRUE & is.null(set.j)){
      info["Added"] <- info["Added"] + 1     
      set.j <- neighbors(graph, i)[ !uneven.neighbors ][[1]]
      if(!is.null(set.j)){search.for.even.neighbor <- FALSE}
    }
    if(!is.null(set.j)){
      if(i != set.j){
        graph <- add_edges(graph, edges=c(i, set.j))
        info["Added"] <- info["Added"] + 1
      }
    }
  }
  (list("graph" = graph, "info" = info))}

Then, consider some network, with 12 nodes

g1 = graph(c(1,2, 1,3, 2,4, 2,5, 1,5, 3,5, 
4,7, 5,7, 5,8, 3,6, 6,8, 6,9, 9,11, 8,11, 
8,10, 8,12, 7,10, 10,12, 11,12), directed = FALSE)

To plot that network, use

V(g1)$name=LETTERS[1:12]
V(g1)$color=rgb(0,0,1,.4)
ly=layout.kamada.kawai(g1)
plot(g1,vertex.color=V(newg)$color,layout=ly)

Then we convert it to some traversable graph by adding 5 vertices

eulerian = make_eulerian(g1)
eulerian$info
     broken       Added Successfull 
          0           5           1 
g = eulerian$graph

as shown below

ly=layout.kamada.kawai(g)
plot(g,vertex.color=V(newg)$color,layout=ly)

We cut those 5 vertices in two part, and therefore, we add 5 artificial nodes

A=as.matrix(as_adj(g))
A1=as.matrix(as_adj(g1))
newA=lower.tri(A, diag = FALSE)*A1+upper.tri(A, diag = FALSE)*A
for(i in 1:sum(newA==2)) newA = cbind(newA,0)
for(i in 1:sum(newA==2)) newA = rbind(newA,0)
s=nrow(A)
for(i in 1:nrow(A)){
  Aj=which(newA[i,]==2)
  if(!is.null(Aj)){
      for(j in Aj){
        newA[i,s+1]=newA[s+1,i]=1
        newA[j,s+1]=newA[s+1,j]=1
        newA[i,j]=1
        s=s+1
      }}}

We get the following graph, where all nodes have an even number of vertices !

newg=graph_from_adjacency_matrix(newA)
newg=as.undirected(newg)
V(newg)$name=LETTERS[1:17]
V(newg)$color=c(rep(rgb(0,0,1,.4),12),rep(rgb(1,0,0,.4),5))
ly2=ly
transl=cbind(c(0,0,0,.2,0),c(.2,-.2,-.2,0,-.2))
for(i in 13:17){
  j=which(newA[i,]>0)
  lc=ly[j,]
  ly2=rbind(ly2,apply(lc,2,mean)+transl[i-12,])
}
plot(newg,layout=ly2)

Our network is now the following (new nodes are small because actually, they don’t really matter, it’s just for computational reasons)

plot(newg,vertex.color=V(newg)$color,layout=ly2,
     vertex.size=c(rep(20,12),rep(0,5)),
     vertex.label.cex=c(rep(1,12),rep(.1,5)))

Now we can get the optimal path

n <- LETTERS[1:nrow(newA)]
g_2 <- new("graphNEL",nodes=n) for(i in 1:nrow(newA)){ for(j in which(newA[i,]>0)){
    g_2 <- addEdge(n[i],n[j],g_2,1) 
  }}
etour(g_2,weighted=FALSE)
 [1] "A" "B" "D" "G" "E" "A" "C" "E" "H" "F" "I" "K" "H" "J" "G" "P" "J" "L" "K" "Q" "L" "H" "O" "F" "C"
[26] "N" "E" "B" "M" "A"

or

edg=attr(E(newg), "vnames")
ET=etour(g_2,weighted=FALSE)
parcours=trajet=rep(NA,length(ET)-1)
for(i in 1:length(parcours)){
  u=c(ET[i],ET[i+1])
  ou=order(u)
  parcours[i]=paste(u[ou[1]],u[ou[2]],sep="|")
  trajet[i]=which(edg==parcours[i])
}
parcours
 [1] "A|B" "B|D" "D|G" "E|G" "A|E" "A|C" "C|E" "E|H" "F|H" "F|I" "I|K" "H|K" "H|J" "G|J" "G|P" "J|P"
[17] "J|L" "K|L" "K|Q" "L|Q" "H|L" "H|O" "F|O" "C|F" "C|N" "E|N" "B|E" "B|M" "A|M"
trajet
 [1]  1  3  8  9  4  2  6 10 11 12 16 15 14 13 26 27 18 19 28 29 17 25 24  7 22 23  5 21 20

Let us try now on a real network of streets. Like Missoula, Montana.

I will not try to get the shapefile of the city, I will just try to replicate the photography above.

If you look carefully, you will see some problem : 10 and 93 have an odd number of vertices (3 here), so one strategy is to connect them (which explains the grey line).

But actually, to be more realistic, we start in 93, and we end in 10. Here is the optimal (shortest) path which goes through all vertices.

Now, we are ready for Halloween, to go through all streets in the neighborhood !

Acheter un billet de train (pas trop cher)

Hier, je suis tombé sur article qui discutait des prix des billets de train, en France (et du prix très élevé, a certaines dates, genre pendant les vacances d’hiver). Tous ceux qui ont l’habitude de prendre le train savent que le prix que l’on paye dépend du moment ou on achète le billet (et de la souplesse que l’on peut avoir sur l’heure (voire la date) du trajet). Cet été, dans le cadre du projet de la formation Data Science pour l’Actuariat, pour mon cours, Pierre proposait de moissonner le site https://www.oui.sncf/ pour suivre un peu l’évolution du prix des billets.

Le principal soucis est que le site https://www.oui.sncf/ s’appuie sur du javascript pour les formulaires de saisie et pour l’affichage des résultats, ce qui empêche l’utilisation classique du package rvest, par exemple. J’avais évoqué dans un autre billet l’utilisation de wdman, pour scraper le site des incendies de forets. Ici, Pierre proposait de passer par casperjs, et je vais reprendre un peu ici sa stratégie:

  • on va utiliser casperjs, un émulateur de navigateur écrit en javascript. Il permet d’émuler un véritable navigateur (même moteur que google chrome) et de résoudre le javascript intégrés à la page
  • on va utiliser un petit code en bash pour lancer le code

Pour le code en bash, c’est juste que je suis sous mac et linux. Sous mac, on peut faire pas mal de choses en bash… Tout passe par des variables, que l’on peut définir, et afficher, par exemple

qui nous donne l’heure. Je peux aussi définir une variable, et l’incrémenter (pratique pour faire des boucles)

Plus intéressant, on peut planifier des taches. Pour ça, on tape

qui va ouvrir un éditeur,

Je demande ici de lancer un code R, tous les heures, a 13:50, 14:50, 15:50, etc. Pour demander tous les jours a 13:50, je tape

on va ensuite sauver l’instruction

On voit que la commande sera lancée, tous les jours : elle est dans la liste des taches a faire

notons qu’on peut lancer un code tout en passant des arguments : ici je lui dis quel objet manipuler (le premier argument) et le second sert pour créer un fichier (et pour le nommer).

Bref, c’est assez facile de lancer automatiquement des codes, pour scraper. Sous windows, on passe par le planificateur de taches. Un premier script permet d’extraire les trains ouverts dans la journée,

et le second, les trains ouverts depuis plus de 24 heures,

Ensuite, on va utiliser http://docs.casperjs.org/en/latest/ pour coder notre émulateur de navigateur internet (le code est ici en ligne).

On va ainsi créer plein de fichiers, contenant les informations que l’on veut ! Je vais passer un peu le retraitement, et juste présenter les informations qu’on peut en tirer. En particulier, Pierre avait stocke les données entre mars et juin dernier.

Ici, on est juste sur quelques trajets de grande ligne,

library(readr)
library(rgdal)
nomFichier = tempfile(fileext = ".zip")
  download.file("https://freakonometrics.free.fr/CarteFrance.zip", destfile = nomFichier, mode = "wb")
  unzip(zipfile = nomFichier, exdir = getwd())
  download.file("https://freakonometrics.free.fr/LgTroncons.csv", destfile = "LgTroncons.csv", mode = "wb")
  download.file("https://freakonometrics.free.fr/CoordVilles.csv", destfile = "CoordVilles.csv", mode = "wb")
  fra0 = readOGR(dsn = paste(getwd(), "/CarteFrance", sep = ""), layer = "gadm36_FRA_0", verbose = F)
  LgTroncons = read_delim("LgTroncons.csv",";", escape_double = FALSE,locale = locale(decimal_mark = ","), trim_ws = TRUE)
CoordVilles = read_delim("CoordVilles.csv",";", escape_double = FALSE,locale = locale(decimal_mark = ","), trim_ws = TRUE)
NomsVilles = CoordVilles[CoordVilles$NOM_A_AFFICHER==1,]
library(ggplot2)
fr_df = fortify(fra0)
ggp = ggplot() + geom_polygon(data=fr_df, aes(long, lat,group = group), fill = "#3A8EBA") 
  ggp = ggp + geom_path(data = LgTroncons, aes(x = LONG, y = LAT, group = ID_TRONCON), colour= "#CC5500", lineend = "round", size=3) + geom_path(data = LgTroncons, aes(x=LONG, y=LAT, group = ID_TRONCON), colour="white", lineend = "round",  size=1.75)
  ggp = ggp + geom_point(data = NomsVilles, aes(x=LONG, y=LAT), colour = "blue", fill = "white", shape=21, size = NomsVilles$PT_SIZE, stroke = NomsVilles$PT_STROKE) + theme_void()
  ggp = ggp + geom_text(data = NomsVilles, aes(x=LONG, y=LAT, label=NOM),hjust = NomsVilles$H_AJUST,
vjust=NomsVilles$V_AJUST, colour = "white", fontface = "bold", size =3.25)+coord_fixed(1.47)
  ggp = ggp + ggtitle("Représentation des trajets étudiées") + theme(plot.title = element_text(hjust = 0.5, face="bold"))
  print(ggp)

On va travailler sur les trajets suivants,

Comparons ici les billets pour un trajet Pars-Rennes, a partir des informations moissonnées pendant 3 mois, le vendredi soir, en particulier 2 vendredi de juin 2018 (les 15 et 22 juin). Pour ces deux jours, il y a eu 6 trains, entre 17 et 20 heures. Pour le 15 juin, les trois premiers ont commencé avec un prix de 45 €. Pour le 22 juin, le premier a commence a 45 €, mais les deux suivant on été lancés a 33 €. Assez rapidement, les prix sont monte a 45 €.

On peut regarder l’evolution du prix

Si on regarde plusieurs destination, on observe des comportement très différentes,

  • pour Le Mans, les prix montent très vite, commençant a 15 €, montant a 18 € 10 heures après le lancement, 21 € le lendemdemain, 27 € au bout d’une semaine. En un mois, les prix ont presque doublé.
  • pour Rennes, on observe une evolution similaire, passant de 20 € a 25 € en quelques heures, et 40 € deux semaines après !
  • pour Toulouse en revanche, le prix initial est plus haut, 43 €, monte de 3 € en 10 heures, 6 € en 16 heures, pour rester a 49 €

Mais pour toutes les destinations, les prix sont croissants.

soit graphiquement

On peut aussi faire une carte. Si on regarde les prix a l’ouverture, Lille, Le Mans et Rennes sont peu chères.

Et les plus fortes variations sur 10 heures sont observées sur Nantes et Bordeaux.

Amusant non ?

Monte Carlo techniques to create counterfactuals

In the previous STT5100 course, last week, we’ve seen how to use monte carlo simulations. The idea is that we do observe in statistics a sample \{y_1,\cdots,y_n\}, and more generally, in econometrics \{(y_1,\mathbf{x}_1),\cdots,(y_n,\mathbf{x}_n)\}. But let’s get back to statistics (without covariates) to illustrate. We assume that observations y_i are realizations of an underlying random variable Y_i. We assume that Y_i are i.id. random variables, with (unkown) distribution F_{\theta}. Consider here some estimator \widehat{\theta} – which is just a function of our sample \widehat{\theta}=h(y_1,\cdots,y_n). So \widehat{\theta} is a real-valued number like . Then, in mathematical statistics, in order to derive properties of the estimator \widehat{\theta}, like a confidence interval, we must define \widehat{\theta}=h(Y_1,\cdots,Y_n), so that now, \widehat{\theta} is a real-valued random variable. What is puzzling for students, is that we use the same notation, and I have to agree, that’s not very clever. So now, \widehat{\theta} is .

There are two strategies here. In classical statistics, we use probability theorem, to derive properties of \widehat{\theta} (the random variable) : at least the first two moments, but if possible the distribution. An alternative is to go for computational statistics. We have only one sample, \{y_1,\cdots,y_n\}, and that’s a pity. But maybe we can create another one \{y_1^{(1)},\cdots,y_n^{(1)}\}, as realizations of F_{\theta}, and another one \{y_1^{(2)},\cdots,y_n^{(2)}\}, anoter one \{y_1^{(3)},\cdots,y_n^{(3)}\}, etc. From those counterfactuals, we can now get a collection of estimators, \widehat{\theta}^{(1)},\widehat{\theta}^{(2)}, \widehat{\theta}^{(3)}, etc. Instead of using mathematical tricks to calculate \mathbb{E}(\widehat{\theta}), compute \frac{1}{k}\sum_{s=1}^k\widehat{\theta}^{(s)}That’s what we’ve seen last friday.

I did also mention briefly that looking at densities is lovely, but not very useful to assess goodness of fit, to test for normality, for instance. In this post, I just wanted to illustrate this point. And actually, creating counterfactuals can we a good way to see it. Consider here the height of male students,

Davis=read.table(
  "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
Davis[12,c(2,3)]=Davis[12,c(3,2)]
X=Davis$height[Davis$sex=="M"]

We can visualize its distribution (density and cumulative distribution)

u=seq(155,205,by=.5)
par(mfrow=c(1,2))
hist(X,col=rgb(0,0,1,.3))
lines(density(X),col="blue",lwd=2)
lines(u,dnorm(u,178,6.5),col="black")
Xs=sort(X)
n=length(X)
p=(1:n)/(n+1)
plot(Xs,p,type="s",col="blue")
lines(u,pnorm(u,178,6.5),col="black")

Since it looks like a normal distribution, we can add the density a Gaussian distribution on the left, and the cdf on the right. Why not test it properly. To be a little bit more specific, I do not want to test if it’s a Gaussian distribution, but if it’s a \mathcal{N}(178,6.5^2). In order to see if this distribution is relevant, one can use monte carlo simulations to create conterfactuals

hist(X,col=rgb(0,0,1,.3))
lines(density(X),col="blue",lwd=2)
  Y=rnorm(n,178,6.5)
  hist(Y,col=rgb(1,0,0,.3))
  lines(density(Y),col="red",lwd=2)
Ys=sort(Y)
plot(Xs,p,type="s",col="white",lwd=2,axes=FALSE,xlab="",ylab="",xlim=c(155,205))
polygon(c(Xs,rev(Ys)),c(p,rev(p)),col="yellow",border=NA)
lines(Xs,p,type="s",col="blue",lwd=2)
lines(Ys,p,type="s",col="red",lwd=2)

We can see on the left that it is hard to assess normality from the density (histogram and also kernel based density estimator). One can hardly think of a valid distance, between two densities. But if we look at graph on the right, we can compare the empirical distribution cumulative distribution \widehat{F} obtained from \{y_1,\cdots,y_n\} (the blue curve), and some conterfactual, \widehat{F}^{(s)} obtained from \{y_1^{(s)},\cdots,y_n^{(s)}\} generated from F_{\theta_0} – where \theta_0 is the value we want to test. As suggested above, we can compute the yellow area, as suggest in Cramer-von Mises test, or the Kolmogorov-Smirnov distance.

d=rep(NA,1e5)
for(s in 1:1e5){
d[s]=ks.test(rnorm(n,178,6.5),"pnorm",178,6.5)$statistic
}
ds=density(d)
plot(ds,xlab="",ylab="")
dks=ks.test(X,"pnorm",178,6.5)$statistic
id=which(ds$x>dks)
polygon(c(ds$x[id],rev(ds$x[id])),c(ds$y[id],rep(0,length(id))),col=rgb(1,0,0,.4),border=NA)
abline(v=dks,col="red")

If we draw 10,000 counterfactual samples, we can visualize the distribution (here the density) of the distance used a test statistic \widehat{d}^{(1)}, \widehat{d}^{(2)}, etc, and compare it with the one observe on our sample \widehat{d}. The proportion of samples where the test-statistics exceeded the one observed

mean(d>dks)
[1] 0.78248

is the computational version of the p-value

ks.test(X,"pnorm",178,6.5)
 
	One-sample Kolmogorov-Smirnov test
 
data:  X
D = 0.068182, p-value = 0.8079
alternative hypothesis: two-sided

I thought about all that a couple of days ago, since I got invited for a panel discussion on “coding”, and why “coding” helped me as professor. And this is precisely why I like coding : in statistics, either manipulate abstract objects, like random variables, or you actually use some lines of code to create counterfactuals, and generate fake samples, to quantify uncertainty. The later is interesting, because it helps to visualize complex quantifies. I do not claim that maths is useless, but coding is really nice, as a starting point, to understand what we talk about (which can be very usefull when there is a lot of confusion on notations).

October, grant proposal season

In 2012, Danielle Herbert, Adrian Barnett, Philip Clarke and Nicholas Graves published an article entitled “on the time spent preparing grant proposals: an observational study of Australian researchers“, whose conclusions had been included in Nature under a more explicit title, “Australia’s grant system wastes time” ! In this study, they included 3700 grant applications sent to the National Health and Medical Research Council, and showed that each application represented 37 working days: “Extrapolating this to all 3,727 submitted proposals gives an estimated 550 working years of researchers’ time (95% confidence interval, 513-589)“. But in these times when I have to write my funding application, I find that losing 37 days of work is huge. Because it’s become the norm! And somehow, it’s sad.

Forget about the crazy idea that I would rather, in fact, spend more time doing my research. In fact, the thought I had this morning was that it is rather sad that in the Faculty of Science, mathematicians are asked to spend a considerable amount of time, comparable to that required of physicists or chemists, for often smaller amounts of funding… And I thought it could be easily verified. We start by retrieving the discipline codes

url="http://www.nserc-crsng.gc.ca/NSERC-CRSNG/FundingDecisions-DecisionsFinancement/ResearchGrants-SubventionsDeRecherche/ResultsGSC-ResultatsCSS_eng.asp"
download.file(url,destfile = "GSC.html")
library(XML)
tables=readHTMLTable("GSC.html")
GSC=tables[[1]]$V1
GSC=as.character(GSC[-(1:2)])
namesGSC=tables[[1]]$V2
namesGSC=as.character(namesGSC[-(1:2)])

We’re going to need a small function, to remove the $ and other symbols that pollute the data (and prevent them from being treated as numbers)

library(stringr)
Correction = function(x) as.numeric(gsub('[$,]', '', x))

We will now read the 12 pages, and harvest (we will just take the 2017 data, but we could go back a few years before)

grants= function(gsc){
     url=paste("http://www.nserc-crsng.gc.ca/NSERC-CRSNG/FundingDecisions-DecisionsFinancement/ResearchGrants-SubventionsDeRecherche/ResultsGSCDetail-ResultatsCSSDetails_eng.asp?Year=2017&GSC=",gsc,sep="")
    download.file(url,destfile = "GSC.html")
    library(XML)
    tables=readHTMLTable("GSC.html")
    X=as.character(tables[[1]]$"Awarded Amount")
    A=as.numeric(Vectorize(Correction)(X))
return(c(median(A),mean(A),as.numeric(quantile(A,(1:99)/100))))
}
M=Vectorize(grants)(GSC[1:12])

The average amounts of individual grants can be compared,

barplot(M[2,])

In mathematics, the average grant amount is $24400. If we normalize by this quantity, we obtain

barplot(M[2,]/M[2,8])

In other words, the average amount of a (individual) grant in chemistry (to pay for students, conferences, etc.) is twice that in mathematics, 60% higher in physics than in maths…

We can also look at the median values (rather than the averages)

barplot(M[1,])

Here again, it is in mathematics that it is the weakest….

barplot(M[1,]/M[1,8])

in comparable proportions. If we think that the time spent writing should be proportional to the amount allocated, we should spend half as much time in math as in chemistry.

Cumulative functions can also be ploted,

plot(M[3:101,8],(1:99)/100,type="s",xlim=range(M))
lines(M[3:101,5],(1:99)/100,type="s",col="red")
lines(M[3:101,4],(1:99)/100,type="s",col="blue")

with math in black, physics in red, and chemistry in blue. What is surprising is the bottom part: a “bad” researcher in chemistry or physics will earn more than the median researcher in mathematics…

Now that my intuition is confirmed, I have to go back, writing my proposal… and explain to my coauthors that I have to postpone some research projects because, well, you know…

Combining automatically factor levels in R

Each time we face real applications in an applied econometrics course, we have to deal with categorial variables. And the same question arise, from students : how can we combine automatically factor levels ? Is there a simple R function ?

I did upload a few blog posts, over the pas years. But so far, nothing satistfying. Let me write down a few lines about what could be done. And if some wants to write a nice R function, that would be awesome. To illustrate the idea, consider the following (simulated dataset)

n=200
set.seed(1)
x1=runif(n)
x2=runif(n)
y=1+2*x1-x2+rnorm(n,0,.2)
LB=sample(LETTERS[1:10])
b=data.frame(y=y,x1=x1,
             x2=cut(x2,breaks=
             c(-1,.05,.1,.2,.35,.4,.55,.65,.8,.9,2),
             labels=LB))
str(b)
'data.frame':	200 obs. of  3 variables:
 $ y : num  1.345 1.863 1.946 2.481 0.765 ...
 $ x1: num  0.266 0.372 0.573 0.908 0.202 ...
 $ x2: Factor w/ 10 levels "I","A","H","F",..: 4 4 6 4 3 6 7 3 4 8 ...
table(b$x2)[LETTERS[1:10]]
 
 A  B  C  D  E  F  G  H  I  J 
11 12 23 34 23 36 12 32  3 14

There is one (continuous) dependent variable y, one continuous covariable x_1 and one categorical variable x_2, with here ten levels. We can plot the data using

plot(b$x1,y,col="white",xlim=c(0,1.1))
text(b$x1,y,as.character(b$x2),cex=.5)

The output of a linear regression yield the following predictions

for(i in 1:10){
p=function(x) predict(lm(y~x1+x2,data=b),newdata=data.frame(x1=x,x2=LETTERS[i]))
u=seq(-1,1.065,by=.01)
v=Vectorize(p)(u)
lines(u,v)}

the slope for x_1 is the same, we simply add a different constant for each level. As we can see, some levels are very very close, so it seems legitimate to combine them into one single category. Here is the output of the linear regression,

summary(lm(y~x1+x2,data=b))
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.843802   0.119655   7.052 3.23e-11 ***
x1           1.992878   0.053838  37.016  < 2e-16 ***
x2A          0.055500   0.131173   0.423   0.6727    
x2H          0.009293   0.121626   0.076   0.9392    
x2F         -0.177002   0.121020  -1.463   0.1452    
x2B         -0.218152   0.130192  -1.676   0.0955 .  
x2D         -0.206970   0.121294  -1.706   0.0896 .  
x2G         -0.407417   0.129999  -3.134   0.0020 ** 
x2C         -0.526708   0.123690  -4.258 3.24e-05 ***
x2J         -0.664281   0.128126  -5.185 5.54e-07 ***
x2E         -0.816454   0.123625  -6.604 3.94e-10 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.2014 on 189 degrees of freedom
Multiple R-squared:  0.8995,	Adjusted R-squared:  0.8942 
F-statistic: 169.1 on 10 and 189 DF,  p-value: < 2.2e-16
AIC(lm(y~x1+x2,data=b))
[1] -60.74443
BIC(lm(y~x1+x2,data=b))
[1] -21.16463

Here the reference category is “I”. And it looks like we could actually combine that category with several others. One strategy here would be to select all categories that seem to be not significantly different, and to run a (multiple) test

library(car)
linearHypothesis(lm(y~x1+x2,data=b), c("x2A = 0", "x2H = 0", "x2F = 0"))
 
Hypothesis:
x2A = 0
x2H = 0
x2F = 0
 
Model 1: restricted model
Model 2: y ~ x1 + x2
 
  Res.Df    RSS Df Sum of Sq      F Pr(>F)    
1    192 8.4651                               
2    189 7.6654  3   0.79971 6.5726  3e-04 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

It seems that we can combine those four categories together.

Here, we can see what’s going on when we change the reference category (actually, loop on all categories)

P=matrix(NA,nlevels(b$x2),nlevels(b$x2))
colnames(P)=rownames(P)=LETTERS[1:10]
plot(1:nlevels(b$x2),1:nlevels(b$x2),col="white",xlab="",ylab="",axes=F,xlim=c(0,10.5),
     ylim=c(0,10.5))
text(1:10,0,LETTERS[1:10])
text(0,1:10,LETTERS[1:10])
for(i in 1:nlevels(b$x2)){
#levels(b$x2)=LETTERS[1:10]
b$x2=relevel(b$x2,LETTERS[i])
p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4]
names(p)=substr(names(p),3,3)
P[LETTERS[i],names(p)]=p
p=P[LETTERS[i],]
idx=which(p>.05)
points(((1:10))[idx],rep(i,length(idx)),pch=1,cex=2)
idx=which(p>.1)
points(((1:10))[idx],rep(i,length(idx)),pch=19,cex=2)}

We are glad to see that it is symmetric : if “H” should be combined with “I”, “I” should also be combined with “H”.

Here black points are related with the 10% p-value, and white points the 5% p-value. This graph is actually hard to read… And actually, this reminds us of  Bertin (1967).

Here, we can predefine manually some ordering (we will see below how it might be automatised)

LETTERSord=c("I","A","H","F","B","D","G","C","J","E")
P=matrix(NA,nlevels(b$x2),nlevels(b$x2))
colnames(P)=rownames(P)=LETTERSord
plot(1:nlevels(b$x2),1:nlevels(b$x2),col="white",xlab="",ylab="",axes=F,xlim=c(0,10.5),
     ylim=c(0,10.5))
ct=c(3,3,2,1,1)
abline(v=.5+c(0,cumsum(ct)),lty=2)
abline(h=.5+c(0,cumsum(ct)),lty=2)
text(1:10,0,LETTERSord)
text(0,1:10,LETTERSord)
for(i in 1:nlevels(b$x2)){
  #levels(b$x2)=LETTERS[1:10]
  b$x2=relevel(b$x2,LETTERSord[i])
  p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4]
  names(p)=substr(names(p),3,3)
  P[LETTERSord[i],names(p)]=p
  p=P[LETTERSord[i],]
  idx=which(p>.05)
  points(((1:10))[idx],rep(i,length(idx)),pch=1,cex=2)
  idx=which(p>.1)
  points(((1:10))[idx],rep(i,length(idx)),pch=19,cex=2)
}

Here we get the following

It looks like we have our combined categories…

Actually, it is possible to use another strategy. We start from some level, say “A”. Then, we merge it with all non-significantly different levels. If “B” is not one of them, we use it as the new reference. Etc.

for(i in 1:nlevels(b$x2)){
  if(LETTERS[i]%in%levels(b$x2)){
  b$x2=relevel(b$x2,LETTERS[i])
  p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4]
  names(p)=substr(names(p),3,nchar(p))
  idx=which(p>.05)
  mix=c(LETTERS[i],names(p)[idx])
  b$x2=recode(b$x2, paste("c('",paste(mix,collapse = "','"),"')='",paste(mix,collapse = "+"),"'",sep=""))
}}

The final categories are

table(b$x2)
 
A+I+H B+D+F   C+G     E     J 
   46    82    35    23    14

with the following regression output

summary(lm(y~x1+x2,data=b))
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.86407    0.03950  21.877  < 2e-16 ***
x1           1.99180    0.05323  37.417  < 2e-16 ***
x2B+D+F     -0.21517    0.03699  -5.817 2.44e-08 ***
x2C+G       -0.50545    0.04528 -11.164  < 2e-16 ***
x2E         -0.83617    0.05128 -16.305  < 2e-16 ***
x2J         -0.68398    0.06131 -11.156  < 2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.2008 on 194 degrees of freedom
Multiple R-squared:  0.8975,	Adjusted R-squared:  0.8948 
F-statistic: 339.6 on 5 and 194 DF,  p-value: < 2.2e-16
AIC(lm(y~x1+x2,data=b))
[1] -66.76939
BIC(lm(y~x1+x2,data=b))
[1] -43.68117

Which is consistent with the group we got before. But actually, if we change the order, we can get different combinations. For instance, if we go from “J” to “A”, instead of “A” to “J”, we obtain

for(i in nlevels(b$x2):1){
  #levels(b$x2)=LETTERS[1:10]
  if(LETTERS[i]%in%levels(b$x2)){
  b$x2=relevel(b$x2,LETTERS[i])
  p=summary(lm(y~x1+x2,data=b))$coefficients[-(1:2),4]
  names(p)=substr(names(p),3,nchar(p))
  idx=which(p>.05)
  mix=c(LETTERS[i],names(p)[idx])
  b$x2=recode(b$x2, paste("c('",paste(mix,collapse = "','"),"')='",paste(mix,collapse = "+"),"'",sep=""))
}}
table(b$x2)
 
          E         G+C I+A+B+D+F+H           J 
         23          35         128          14

with different information criteria here

AIC(lm(y~x1+x2,data=b))
[1] -36.61665
BIC(lm(y~x1+x2,data=b))
[1] -16.82675

I guess it would be necessary to run randomly the order we go through the levels. Last, but not least, one can use regression trees (even if it not per se in the syllabus of the course). The problem is that there is another explanatory variable that might interphere. So I would suggest (1) to fit a linear model y=\beta_0+\beta_1x_1+u_i, to calculate the residuals, \widehat{u}_i (2) to run a regression tree, to explain \widehat{u}_i with categorical variable x_2 (I did explain how trees are build when the explanatory variable is a categorical one in a previous post)

library(rpart)
library(rpart.plot)
b$e=residuals(lm(y~x1,data=b))
arbre=rpart(e~x2,data=b)
prp(arbre,type=2,extra=1)

Observe that the leaves have the same groups as the one we got.

arbre
n= 200 
 
node), split, n, deviance, yval
      * denotes terminal node
 
1) root 200 22.563500  7.771561e-18  
  2) x2=G,C,J,E 72  4.441495 -3.232525e-01  
    4) x2=J,E 37  1.553520 -4.578492e-01 *
    5) x2=G,C 35  1.509068 -1.809646e-01 *
  3) x2=I,A,H,F,B,D 128  6.366628  1.818295e-01  
    6) x2=F,B,D 82  2.983381  1.048246e-01 *
    7) x2=I,A,H 46  2.030229  3.190993e-01 *

I guess that it should be possible to put all that in an R function, to suggest combinations of level that might improve the regression.

Traitement des valeurs manquantes : remplacer les NA par une constante ?

Un rapide billet pour répondre à une question posée à la fin du cours de ce matin (ST5100), par Jean-Pierre Liégeois, jeune lecteur du Var (pour préserver un peu d’anonymat)

Dans mon stage, quand on avait des valeurs manquantes, on me disait de remplacer par -1, puis de rajouter une indicatrice comme quoi la variable vaut -1. Ça permet de ne supprimer ni variable, ni observations. On peut faire ça ?

Si je formalise un peu, on va simuler ici des données, disons x_1 et x_2, on génère ensuite des données suivant un modèle, de la forme y=\beta_0+\beta_1x_1+\beta_2x_2+\varepsilon. Une proportion \alpha de x_1 sera transformée en NA.  Ce que suggérais Jean-Pierre, c’est de remplacer les valeurs manquantes par -1, puis d’ajuster un modèle y=\beta_0+\beta_1x_1+\beta_{-1}\mathbf{1}(x_1=-1)+\beta_2x_2+\varepsilon. Côté code, c’est assez simple. Par défaut, la stratégie de R est de supprimer les valeurs manquantes. Si 50% des données de x_1 sont manquantes, la moitié des lignes sont supprimées

n=1000
x1=runif(n)
x2=runif(n)
e=rnorm(n,.2)
y=1+2*x1-x2+e
alpha=.05
indice=sample(1:n,size=round(n*alpha))
base=data.frame(y=y,x1=x1)
base$x1[indice]=NA
reg=lm(y~x1+x2,data=base)

Au lieu de générer un unique échantillon, on va en simuler 10,000, et regarder la distribution de \widehat{\beta}_1,

m=10000
B=rep(NA,m)
for(s in 1:m){
  x1=runif(n)
  x2=runif(n)
  e=rnorm(n,.2)
  y=1+2*x1-x2+e
  alpha=.5
  indice=sample(1:n,size=round(n*alpha))
  base=data.frame(y=y,x1=x1,x2=x2)
  base$x1[indice]=NA
  reg=lm(y~x1+x2,data=base)
  B[s]=coefficients(reg)[2]
}
hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white",xlab="missing values = 50%")
lines(density(B),lwd=2,col="blue")
abline(v=2,lty=2,col="red")

Bien sur, avec un taux de valeurs manquantes plus faibles – disons \alpha=5\% – on perd moins d’observations, et donc l’estimateur a une variance plus faible.

Tentons maintenant la stratégie consistant à remplacer les valeurs manquantes par des valeurs numériques fixes, et de rajouter une indicatrice,

B=rep(NA,m)
for(s in 1:m){
  x1=runif(n)
  x2=runif(n)
  e=rnorm(n,.2)
  y=1+2*x1-x2+e
  alpha=.5
  indice=sample(1:n,size=round(n*alpha))
  base=data.frame(y=y,x1=x1,x2=x2)
  base$x1[indice]=-1
  reg=lm(y~x1+x2+I(x1==(-1)),data=base)
  B[s]=coefficients(reg)[2]
}
hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white")
lines(density(B),lwd=2,col="blue")
abline(v=2,lty=2,col="red")

Ce qui ne change pas grand chose, on en conviendra…  y compris si le taux de valeurs manquantes passe à 5%,

On peut se demander ce qui se passe si le shift n’est plus de 1 mais de 10 (a priori, c’est arbitraire, la variable x_1  pouvant être plus ou moins dispersée… -1 pour une variable entre 0 et 1, ou entre 0 et 1000, ça n’est pas tout à fait pareil). Mais non, par exemple avec toujours 5% de valeurs manquantes, on a

Si on regarde notre échantillon, en particulier le nuage de points  (x_1,y), on observe

ici, les valeurs manquantes sont choisies au hasard, de manière totalement indépendante,

x1=runif(n)
x2=runif(n)
e=rnorm(n,.2)
y=1+2*x1-x2+e
alpha=.3333333
indice=sample(1:n,size=round(n*alpha))
clr=rep("black",n)
clr[indice]="red"
plot(x1,y,col=clr)

(ici avec 1/3 de valeurs manquantes, en rouge). Mais on pourrait supposer que les valeurs manquantes sont les plus grandes valeurs de x_1, par exemple,

x1=runif(n)
x2=runif(n)
e=rnorm(n,.2)
y=1+2*x1-x2+e
alpha=.3333333
indice=sample(1:n,size=round(n*alpha),prob = x1^3)
clr=rep("black",n)
clr[indice]="red"
plot(x1,y,col=clr)

On peut se demander ce que ça donnerait sur l’estimateur \widehat{\beta}_1

Ça ne change pas grand chose, mais on a plus de variance, si on regarde bien. Dernier essai : que se passe-t-il si les variables x_1 et x_2 sont maintenant corrélées,

B=rep(NA,m)
library(mnormt)
r=.8
S = matrix(c(1,r,r,1),2,2)
for(s in 1:m){
  x=rmnorm(n,varcov = S)
  x1=pnorm(x[,1])
  x2=pnorm(x[,2])
  e=rnorm(n,.2)
  y=1+2*x1-x2+e
  alpha=.2
  indice=sample(1:n,size=round(n*alpha),prob = x1^3)
  base=data.frame(y=y,x1=x1,x2=x2)
  base$x1[indice]=-1
  reg=lm(y~x1+x2+I(x1==(-1)),data=base)
  B[s]=coefficients(reg)[2]
}
hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white")
lines(density(B),lwd=2,col="blue")
abline(v=2,lty=2,col="red")

Cette fois, on a un estimateur biaisé (de l’ordre de 10% sur cet exemple numérique). Manifestement, cette technique n’est pas très concluante…

Je pourrais ajouter que cette méthode ne revient pas à la première, même si la distribution des estimateurs est proche

set.seed(1)
x=rmnorm(n,varcov = S)
x1=pnorm(x[,1])
x2=pnorm(x[,2])
e=rnorm(n,.2)
y=1+2*x1-x2+e
alpha=.2
indice=sample(1:n,size=round(n*alpha),prob = x1^3)
base=data.frame(y=y,x1=x1,x2=x2)
base$x1[indice]=-1
reg1=lm(y~x1+x2+I(x1==(-1)),data=base)
coefficients(reg1)
      (Intercept)                x1                x2 I(x1 == (-1))TRUE 
        1.0988005         1.7454385        -0.5149477         3.1000668 
base$x1[indice]=NA
reg2=lm(y~x1+x2,data=base)
coefficients(reg2)
(Intercept)          x1          x2 
  1.1123953   1.8612882  -0.6548206

Comme je le disais (lors de la discussion qui a suivi le cours) une méthode plus prometteuse est l’imputation. L’idée est de prédire une valeur pour les x_1 manquants. On pourrait être tenté de prendre  \overline{x}_1, la moyenne des x_1 observés. Mais on sait que les valeurs manquantes sont justement pour les grandes valeurs de x_1, ici, donc on doit pouvoir faire mieux ! On sait aussi que x_1  et x_2 sont corrélés ici. Corrélés positivement, en plus. Autrement dit, si x_2 est grand, on sait que le x_1 (non observé) devait être grand. Le plus simple est de faire un modèle linéaire, x_1=\alpha_0+\alpha_2x_2+\eta_i, calibré sur les valeurs non-manquantes. Puis on utilise \widehat{x}_1=\widehat{\alpha}_0+\widehat{\alpha}_2x_2 pour les valeurs manquantes. C’est simpliste, mais pourquoi pas ? On estime alors le modèle sur cette nouvelle base.

for(s in 1:m){
  x=rmnorm(n,varcov = S)
  x1=pnorm(x[,1])
  x2=pnorm(x[,2])
  e=rnorm(n,.2)
  y=1+2*x1-x2+e
  alpha=.2
  indice=sample(1:n,size=round(n*alpha),prob = x1^3)
  base=data.frame(y=y,x1=x1,x2=x2)
    base$x1[indice]=NA
    reg0=lm(x1~x2,data=base[-indice,])
    base$x1[indice]=predict(reg0,newdata=base[indice,])
  reg=lm(y~x1+x2,data=base)
  B[s]=coefficients(reg)[2]
}
hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white")
lines(density(B),lwd=2,col="blue")
abline(v=2,lty=2,col="red")

sur l’exemple numérique, on obtient

base$x1[indice]=NA
reg0=lm(x1~x2,data=base[-indice,])
base$x1[indice]=predict(reg0,newdata=base[indice,])
reg3=lm(y~x1+x2,data=base)
coefficients(reg3)
(Intercept)          x1          x2 
  1.1593298   1.8612882  -0.6320339

Cette méthode a au moins pu corriger du biais…

Après, si on regarde attentivement, on a exactement la même valeur qu’avec la première méthode qui consiste à supprimer les lignes avec des valeurs manquantes !

summary(reg3)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.15933    0.06649  17.435  < 2e-16 ***
x1           1.86129    0.21967   8.473  < 2e-16 ***
x2          -0.63203    0.20148  -3.137  0.00176 ** 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.051 on 997 degrees of freedom
Multiple R-squared:  0.1094,	Adjusted R-squared:  0.1076 
F-statistic: 61.23 on 2 and 997 DF,  p-value: < 2.2e-16 
 
summary(reg2) 
 
Coefficients: Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.11240    0.06878  16.173  < 2e-16 ***
x1           1.86129    0.21666   8.591  < 2e-16 ***
x2          -0.65482    0.20820  -3.145  0.00172 ** 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.037 on 797 degrees of freedom
  (200 observations deleted due to missingness)
Multiple R-squared:  0.1223,	Adjusted R-squared:   0.12 
F-statistic:  55.5 on 2 and 797 DF,  p-value: < 2.2e-16

Au lieu de faire une régression linéaire, on peut utiliser une autre méthode d’imputation. Par exemple prendre la moyenne des k valeurs de x_1 (observées) pour les k individus ayant des x_2 les plus proches du x_2 de l’individu ayant x_1 manquant :

kpp=function(i,basena,k=5){
  x2=basena$x2[i]
  sb=basena[!is.na(basena$x1),]
  idx=rank(abs(sb$x2-x2))
  mean(sb[which(idx<=k),"x1"])
}

Sur notre base simulée on obtient

base$x1[indice]=NA
base0=base
for(j in indice) base0$x1[j]=kpp(j,base0,k=5)
reg4=lm(y~x1+x2,data=base)
coefficients(reg4)
(Intercept)          x1          x2 
   1.197944    1.804220   -0.806766

Si on regarde ce que ça donne sur nos 10,000 simulations, on a (c’est un peu long, car j’ai codé ça très rapidement, et pas du tout de manière optimale)

for(s in 1:m){
  x=rmnorm(n,varcov = S)
  x1=pnorm(x[,1])
  x2=pnorm(x[,2])
  e=rnorm(n,.2)
  y=1+2*x1-x2+e
  alpha=.2
  indice=sample(1:n,size=round(n*alpha),prob = x1^3)
  base=data.frame(y=y,x1=x1,x2=x2)
  base$x1[indice]=NA
  base0=base
    for(j in indice) base0$x1[j]=kpp(j,base0,k=5)
  reg=lm(y~x1+x2,data=base0)
  B[s]=coefficients(reg)[2]
}
hist(B,probability=TRUE,col=rgb(0,0,1,.4),border="white")
lines(density(B),lwd=2,col="blue")
abline(v=2,lty=2,col="red")

Le biais semble ici plus faible que sans imputation… autrement dit, les méthodes d’imputation me semblent plus robustes que cette stratégie visant à remplacer des NA par une valeur arbitraire, et de rajouter une indicatrice dans la régression.

« Dans toute statistique, l’inexactitude du nombre est compensée par la précision des décimales »

Le statisticien et économiste Alfred Sauvy est resté dans les mémoires pour avoir inventé en 1952 le terme “tiers-monde”. Mais on lui a aussi attribué la paternité de la phrase suivante « dans toute statistique, l’inexactitude du nombre est compensée par la précision des décimales ». J’ai du l’entendre alors que j’étais étudiant, et depuis, elle me suit partout.

J’y repensais l’autre jour, quand Mathieu Gallard mentionnait sur Twitter le graphique suivant (correspondant a la popularité du président de la république, en France, dans les 18 mois qui suivent l’élection).

Le tweet disait (entre autres) “ est à ce stade de son quinquennat légèrement moins populaire que “. Pour rappel, ces courbes “de popularité” sont construites par un sondage, avec a chaque fois environ 1000 personnes interrogées (“sur la taille de l’échantillon on est toujours entre 950 et 970 interviews” me disait Mathieu). Bon, tous ceux qui ont des souvenirs de cours de stats se souviennent qu’avec 1000 personnes interrogées, 2 points de différence, c’est rarement significatif. Mais plus globalement, compte tenu de la marge d’erreur, je me suis demande pourquoi les courbes n’étaient pas lissées ? Ça éviterait les discussions stériles pour une variation de 2 points par exemple..

Si on reprend les données brutes (merci Mathieu), on a ici

rate=read.csv2("http://freakonometrics.free.fr/satisfaction.csv")
plot(rate[,2],type="b",ylim=c(.2,.75),xlim=c(0,20),pch=19)
lines(rate[,3],type="b",col="red")
lines(rate[,4],type="b",col="blue")
lines(rate[,5],type="b",col="dark green")
text(18.15,rate[17,2],"EM")
text(18.15,rate[17,3],"FH",col="red")
text(18.15,rate[17,4],"NS",col="blue")
text(18.15,rate[17,5],"JC",col="dark green")

ce qui donne le même que dans le tweet (même si je n’ose pas interpoler linéairement les valeurs manquantes – il y en a deux dans mon fichier)

Prenons la courbe la plus récente, celle d’Emmanuel Macron (avec en plus les valeurs manquantes pour pimenter un peu) et rajoutons les intervalles de confiance ponctuels.

plot(rate$EM,type="b",ylim=c(.2,.5))
p=rate$EM
n=length(p)
arrows(1:n,p-2/sqrt(1000)*sqrt(p*(1-p)),1:n,p+2/sqrt(1000)*sqrt(p*(1-p)),code=3,angle=90,length=.1,col="blue")

Personnellement, j’aurais bien voulu (1) lisser tout ça, (2) rajouter quelque chose qui s’apparente a des bandes de confiance. Mais avec des erreurs de mesure (c’est comme ça qu’on peut interpréter le fait que les points viennent d’un sondage), je ne sais pas trop quoi faire. J’ai tenté la méthode suivante : le tire au hasard des points dans l’intervalle de confiance, puis je lisse sur ce nouveau jeu de points. Et je répète mille fois

library(mgcv)
Y=matrix(NA,73,1000)
for(s in 1:1000){
  x=(1:n)[!is.na(p)]
  pna=p[!is.na(p)]
  ps=rnorm(length(x),pna,1/sqrt(1000)*sqrt(pna*(1-pna)))
  b=data.frame(x=x,y=ps)
  reg=gam(y~s(x),data=b)
  yp=predict(reg,newdata=data.frame(x=seq(0,18,by=.25)))
  Y[,s]=yp
  if(s<100) lines(seq(0,18,by=.25),yp,col="light blue")
} 
lines(seq(0,18,by=.25),apply(Y,1,mean),col="red",lwd=2)
lines(seq(0,18,by=.25),apply(Y,1,function(x) quantile(x,.95)),col="red",lty=2)
lines(seq(0,18,by=.25),apply(Y,1,function(x) quantile(x,.05)),col="red",lty=2)

On voit que notre courbe lissée est réaliste, voire même les pseudo-bandes de confiance autour. Pour obtenir ces trois courbes, on peut utiliser la fonction suivante

courbe=function(j=1){
p=rate[,1+j]
n=length(p)
Y=matrix(NA,73,1000)
for(s in 1:1000){
  x=(1:n)[!is.na(p)]
  pna=p[!is.na(p)]
  ps=rnorm(length(x),pna,1/sqrt(1000)*sqrt(pna*(1-pna)))
  b=data.frame(x=x,y=ps)
  reg=gam(y~s(x),data=b)
  yp=predict(reg,newdata=data.frame(x=seq(0,18,by=.25)))
  Y[,s]=yp
} 
data.frame(
x=seq(0,18,by=.25),
pred=apply(Y,1,mean),
upr=apply(Y,1,function(x) quantile(x,.975)),
lwr=apply(Y,1,function(x) quantile(x,.025)))
}

Sur les quatre colonnes de notre tableau, ça donne

plot(rate[,2],type="b",ylim=c(0,.8),xlim=c(0,20),col="white")
Y=courbe(4)
polygon(c(Y$x,rev(Y$x)),c(Y$upr,rev(Y$lwr)),col=rgb(0,1,0,.4),border=NA)
lines(Y$x,Y$pred,col="dark green",lwd=2)
text(18.65,Y$pred[73],"JC",col="dark green")
Y=courbe(3)
polygon(c(Y$x,rev(Y$x)),c(Y$upr,rev(Y$lwr)),col=rgb(0,0,1,.4),border=NA)
lines(Y$x,Y$pred,col="blue",lwd=2)
text(18.65,Y$pred[73],"NS",col="blue")
Y=courbe(2)
polygon(c(Y$x,rev(Y$x)),c(Y$upr,rev(Y$lwr)),col=rgb(1,0,0,.4),border=NA)
lines(Y$x,Y$pred,col="red",lwd=2)
text(18.65,Y$pred[73],"FH",col="red")
Y=courbe(1)
polygon(c(Y$x,rev(Y$x)),c(Y$upr,rev(Y$lwr)),col="grey",border=NA)
lines(Y$x,Y$pred,col="black",lwd=2)
text(18.65,Y$pred[73],"EM",col="black")

Pourquoi les instituts de sondages, qui produisent les courbes de popularité, ne montrent pas ce genre de courbes ? Elles sont – a mon avis – aussi justes que celles qu’ils fournissent, au centième près, jouant sur une précision que l’incertitude ne devrait pas autoriser…

Régression sur une variable qualitative et ANOVA

Ce matin, pour le cours STT5100, on évoquait la régression sur une variable catégorielle. En particulier, on avait commencé par regarder ce que donnerait la régression sans la constante, et son interprétation. On s’était appuyé sur la base des poids et des tailles des élèves, et de la variable de genre.

Davis=read.table(
  "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
Davis[12,c(2,3)]=Davis[12,c(3,2)]
Davis=data.frame(Y=Davis$weight * 2.204622,
                 X1=Davis$sex)

On voulait estimer le modèle y_i =\beta_F\boldsymbol{1}_F(x_i)+\beta_H\boldsymbol{1}_H(x_i)+\varepsilon_iOn avait vu que l’on pouvait passer par l’écriture matricielle

 X=cbind(Davis$X1=='F',Davis$X1=='M') 
 Y=Davis$Y

car la matrice \mathbf{X}^T\mathbf{X} est inversible (une fois que l’on enlève la constante)

 solve(t(X)%*%X)
            [,1]       [,2]
[1,] 0.008928571 0.00000000
[2,] 0.000000000 0.01136364

et donc l’estimateur par moindres carrés est (classiquement)\widehat{\mathbf{\beta}} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}

 solve(t(X)%*%X) %*% (t(X)%*%Y)
         [,1]
[1,] 125.4272
[2,] 167.3258

ce qui correspond effectivement à la sortie de R,

 reg=lm(Y~0+X1,data=Davis)
 summary(reg)
 
Coefficients:
    Estimate Std. Error t value Pr(>|t|)    
X1F  125.427      1.960   64.00   <2e-16 ***
X1M  167.326      2.211   75.68   <2e-16 ***

Considérons maintenant les deux sous-populations, avec le poids des femmes, et le poids des hommes

x=Y[X[,1]==1]
y=Y[X[,2]==1]
nx=length(x)
ny=length(y)

On avait vu en cours que les \widehat{\mathbf{\beta}} avaient une interprétation très simple, puisque\widehat{{\beta}}_M = \frac{1}{n_M}\sum_{i:x_i=M} y_iautrement dit \widehat{{\beta}}_M   est le poids moyen des hommes. Et en effet

 mean(y)
[1] 167.3258

C’est finalement très naturel, ou intuitif.

On peut maintenant s’interroger sur l’écart-type de l’estimateur de \widehat{{\beta}}_M . Intuitivement, on aurait envie d’avoir la variance de l’estimateur de la moyenne, soit ici

 sqrt(var(y)/ny)
[1] 2.794391
 sqrt(1/(ny-1)*sum( (y-mean(y))^2 )/ny)
[1] 2.794391

car pour rappel\text{Var}[\overline{y}]=\frac{\text{Var}(y)}{n}Comme on l’a vue dans le modèle de régression multiple, la variance de l’estimateur de \mathbf{\beta} est proportionnel à \sigma^2 , la variance globale des résidus (c’est l’hypothèse d’homoscédasticité ! les deux groupes doivent avoir la même variance). On va donc calculer l’estimateur naturel de \sigma^2

 s2=1/(nx+ny-2)*(sum( (x-mean(x))^2 )+sum( (y-mean(y))^2))
 sqrt(s2/ny)
[1] 2.210863

et en effet, on retombe sur la valeur donnée dans le tableau de régresion

 sqrt(s2/nx)
[1] 1.959721

(pareil pour l’autre coefficient).

On avait ensuite regardé la régression telle qu’elle faite classiquement, sous R : on garde la constante, et on enlève une des variables indicatrices (qui devient alors la “modalité de référence”).

 X=cbind(1,Davis$X1=='M')

Là encore, le modèle devient identifiable, et on obtient ici

 solve(t(X)%*%X) %*% (t(X)%*%Y)
          [,1]
[1,] 125.42724
[2,]  41.89855

On avait noté qu’il y avait un interprétation de cette seconde valeur, comme un différentiel par rapport à la modalité de référence

mean(y)-mean(x)
[1] 41.89855

La sortie de régression devient ici

 reg2=lm(Y~X1,data=Davis)
 summary(reg2)
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  125.427      1.960   64.00   <2e-16 ***
X1M           41.899      2.954   14.18   <2e-16 ***

Et comme je l’avais dit, le test de Student correspond ici à un test d’égalité entre la taille moyenne des hommes et celle des femmes. Et en effet, si on fait le test, on voit que la différence est significative, comme attendu (pour la même raison qu’au dessus, on suppose la même variance dans les deux groupes)

 t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE)
 
	Two Sample t-test
 
data:  Y[X[, 1] == 1] and Y[X[, 2] == 1]
t = -6.4475, df = 286, p-value = 4.826e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -30.62603 -16.30035
sample estimates:
mean of x mean of y 
 143.8626  167.3258

Je suis par contre un peu surpris que les p-values soient différente. Mon interprétation est que les p-values sont (de toutes façons) très faibles, et donc ça a peu d’importance. En fait, si on rend les deux variables indépendantes (par exemple en mélangeant la variable \mathbf{y} ), ça marche ! Posons

 Davis$Y=sample(Davis$Y)

ce qui revient à permuter toutes les observations de la variable dépendante (mais pas les autres !). La régression donne ici

 reg2=lm(Y~X1,data=Davis)
 summary(reg2)
 
Call:
lm(formula = Y ~ X1, data = Davis)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-57.458 -22.184  -5.512  17.809 118.912 
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 143.4382     2.7820   51.56   <2e-16 ***
X1M           0.9645     4.1940    0.23    0.818    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 29.44 on 198 degrees of freedom
Multiple R-squared:  0.000267,	Adjusted R-squared:  -0.004782 
F-statistic: 0.05289 on 1 and 198 DF,  p-value: 0.8183

autrement dit, le genre n’est plus significatif, avec une p-value de 81.8%. Ce qui est bien au dessus de 5%. Si on fait maintenant le test de comparaison de moyenne, sur les deux sous-groupes, on obtient

 Y=Davis$Y
 t.test(Y[X[,1]==1],Y[X[,2]==1],var.equal=TRUE)
 
	Two Sample t-test
 
data:  Y[X[, 1] == 1] and Y[X[, 2] == 1]
t = -0.22998, df = 198, p-value = 0.8183
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -9.235209  7.306165
sample estimates:
mean of x mean of y 
 143.4382  144.4027

et le test a ici également une p-value de 81.8%. Les deux tests sont donc rigoureusement équivalents.

Les transports publics parisiens

Histoire de continuer la série de billets sur la visualisation, et la manipulation de données ouvertes, je vais reprendre de codes de Tony, de la formation Data Science pour l’Actuariat, pour visualiser le transport dans Paris (et la région parisienne). Si j’ai le temps, dans les jours a venir, je ferais une analyse du réseau de métro, compare aux autres grandes villes européennes. Pour commencer, on va récupérer les données, fournies par le site d’open data du stif, le syndicat des transports d’Ile de France (https://opendata.stif.info). Les données sont découpées par semestre, ce qui rend le code un peu lourd… mais bon, ça n’est pas plus complique pour autant.

library(dplyr)
library(stringr)
library(ggplot2)
library(xlsx)
library(ggmap)

On commence par lire tous les fichiers en ligne

nbvalid = list()
download.file("https://opendata.stif.info/explore/dataset/emplacement-des-gares-idf-data-generalisee/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true","Gares.csv")
gares = read.csv("Gares.csv", sep=";", header = TRUE)
distr_pers = list()
download.file("https://opendata.stif.info/explore/dataset/validations-sur-le-reseau-ferre-profils-horaires-par-jour-type-1er-sem/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true","horaires1.csv")
distr_pers$S1 = read.csv("horaires1.csv", sep=";", header = TRUE)
download.file("https://opendata.stif.info/explore/dataset/validations-sur-le-reseau-ferre-profils-horaires-par-jour-type-2e-sem/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true","horaires2.csv")
distr_pers$S2 = read.csv("horaires2.csv", sep=";", header = TRUE)
download.file("https://opendata.stif.info/explore/dataset/validations-sur-le-reseau-ferre-nombre-de-validations-par-jour-1er-sem/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true","validations1.csv")
nbvalid$S1 = read.csv("validations1.csv", sep=";", header = TRUE)
download.file("https://opendata.stif.info/explore/dataset/validations-sur-le-reseau-ferre-nombre-de-validations-par-jour-2e-sem/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true","validations2.csv")
nbvalid$S2 = read.csv("validations2.csv", sep=";", header = TRUE)
download.file("https://freakonometrics.free.fr/Correspondance_NOM.csv","Correspondance_NOM.csv")
Cooresp = read.csv("Correspondance_NOM.csv", sep=";", header = TRUE)

On a ensuite besoin de définir les dates des vacances, pour 2017

Vacances = list()
Vacances$Noel = append(seq(from = as.Date("01/01/2017", format="%d/%m/%Y"), to=as.Date("02/01/2017", format="%d/%m/%Y"), by=1),seq(from = as.Date("24/12/2017", format="%d/%m/%Y"), to=as.Date("31/12/2017", format="%d/%m/%Y"), by=1))
Vacances$Ski = seq(from = as.Date("04/02/2017", format="%d/%m/%Y"), to=as.Date("19/02/2017", format="%d/%m/%Y"), by=1)
Vacances$Printemps = seq(from = as.Date("02/04/2017", format="%d/%m/%Y"), to=as.Date("17/04/2017", format="%d/%m/%Y"), by=1)
Vacances$Ete = seq(from = as.Date("08/07/2017", format="%d/%m/%Y"), to=as.Date("03/09/2017", format="%d/%m/%Y"), by=1)
Vacances$Toussaint = seq(from = as.Date("21/10/2017", format="%d/%m/%Y"), to=as.Date("05/11/2017", format="%d/%m/%Y"), by=1)
Vacances$All=Reduce(append,Vacances)

Après, un peu de nettoyage est nécessaire, avec des gares en double (par exemple quand passe a la fois le RER et le métro), et pour recuperer leur localisation spatiale (latitude et longitude)

gares$NOM_LONG = as.character(gares$NOM_LONG)
DD = (gares$NOM_LONG[duplicated(gares$NOM_LONG)])
i = (gares$NOM_LONG %in% DD) & gares$MODE_=="Metro"
gares$NOM_LONG[i] = paste(gares$NOM_LONG[i],"M", sep="-")
i = (gares$NOM_LONG %in% DD) & gares$MODE_=="RER"
gares$NOM_LONG[i] = paste(gares$NOM_LONG[i],"R", sep="-")
gares$NOM_LONG=factor(gares$NOM_LONG)
 
a=as.character(gares$Geo.Point)
gares$Y=as.numeric(str_extract_all(a,"^[0-9]+.[0-9]+"))
gares$X=as.numeric(str_extract_all(a,"[0-9]+.[0-9]+$"))

On compte ensuite les nombres de validation de tickets, par gare

Manip_nbvalid = function(Data,DD,gares) {
  i=grep("^[a-zA-Z]+",as.character(Data$NB_VALD))
  Data$NB_VALD[i]=as.integer(5)
  i=is.na(Data$NB_VALD)
  Data$NB_VALD[i]=as.integer(5)
  Data$LIBELLE_ARRET=as.character(Data$LIBELLE_ARRET)
  i=(Data$LIBELLE_ARRET %in% DD) & Data$CODE_STIF_TRNS=="100"
  Data$LIBELLE_ARRET[i]=paste(Data$LIBELLE_ARRET[i],"M", sep="-")
  i=(Data$LIBELLE_ARRET %in% DD) & Data$CODE_STIF_TRNS=="800"
  Data$LIBELLE_ARRET[i]=paste(Data$LIBELLE_ARRET[i],"R", sep="-")
 
  for (i in seq(1,nrow(Cooresp))) { Data$LIBELLE_ARRET=gsub(as.character(Cooresp$nbval[i]),as.character(Cooresp$gares[i]),Data$LIBELLE_ARRET)
  }
gares$NOM_LONG=as.character(gares$NOM_LONG)
Data=dplyr::left_join(Data,gares[,c("NOM_LONG","X","Y")],by=c("LIBELLE_ARRET"="NOM_LONG"))
  Data=Data[is.na(Data$CODE_STIF_ARRET)==FALSE,]
  Data=Data[Data$CODE_STIF_ARRET!="ND",]
  Data$NB_VALD=as.integer(as.character(Data$NB_VALD))
  Data$JOUR=as.Date(Data$JOUR)
  Data$CODE_STIF_TRNS=factor(Data$CODE_STIF_TRNS)
  Data$CODE_STIF_RES=factor(Data$CODE_STIF_RES)
  Data$CODE_STIF_ARRET=factor(Data$CODE_STIF_ARRET)
  Data$LIBELLE_ARRET=factor(Data$LIBELLE_ARRET)
  Data$ID_REFA_LDA=factor(Data$ID_REFA_LDA)
  Data$CATEGORIE_TITRE=factor(Data$CATEGORIE_TITRE)
  Data$JOURSEM=weekdays(Data$JOUR)  
  return(Data)
}
nbvalid=lapply(nbvalid, Manip_nbvalid,DD=DD,gares=gares)

On a ainsi tous les comptages, pour toutes les gares. On fait ensuite un découpage par tranche horaire

Manip_dist_pers = function(DataFrame) {
  DataFrame=DataFrame[(DataFrame$TRNC_HORR_60)!="ND",]
  DataFrame$TRNC_HORR_60=factor(DataFrame$TRNC_HORR_60, levels = c("0H-1H", "1H-2H", "2H-3H", "3H-4H", "4H-5H", "5H-6H", "6H-7H", "7H-8H", "8H-9H", "9H-10H", "10H-11H", "11H-12H", "12H-13H", "13H-14H", "14H-15H", "15H-16H", "16H-17H", "17H-18H", "18H-19H", "19H-20H", "20H-21H", "21H-22H", "22H-23H", "23H-0H")) 
  DataFrame=DataFrame[(DataFrame$CODE_STIF_ARRET)!="ND",]
  DataFrame$CODE_STIF_ARRET=factor(DataFrame$CODE_STIF_ARRET)
DataFrame$TRANCHE=str_extract(as.character(DataFrame$TRNC_HORR_60),"^[0-9]{1,2}")
  return(DataFrame)
}
distr_pers=lapply(distr_pers, Manip_dist_pers)

On peut ensuite recuperer la distribution des validation, par jour

distr_JOURV=list()
distr_JOURV$S1 = nbvalid$S1 %>% group_by(JOUR, JOURSEM,CATEGORIE_TITRE) %>% summarise(NB_VALD=sum(NB_VALD))
distr_JOURV$S2 = nbvalid$S2 %>% group_by(JOUR, JOURSEM,CATEGORIE_TITRE) %>% summarise(NB_VALD=sum(NB_VALD))
distr_JOURV$Y=rbind(distr_JOURV$S1,distr_JOURV$S2)
distr_JOUR=list()
distr_JOUR$S1 = nbvalid$S1 %>% group_by(JOUR, JOURSEM) %>% summarise(NB_VALD=sum(NB_VALD))
distr_JOUR$S2 = nbvalid$S2 %>% group_by(JOUR, JOURSEM) %>% summarise(NB_VALD=sum(NB_VALD))
distr_JOUR$Y=rbind(distr_JOUR$S1,distr_JOUR$S2)
distr_JOUR_Station=list()
distr_JOUR_Station$S1 = nbvalid$S1 %>% group_by(JOUR, JOURSEM,CODE_STIF_ARRET,LIBELLE_ARRET) %>% summarise(NB_VALD=sum(NB_VALD), X=max(X), Y=max(Y))
distr_JOUR_Station$S2 = nbvalid$S2 %>% group_by(JOUR, JOURSEM,CODE_STIF_ARRET,LIBELLE_ARRET) %>% summarise(NB_VALD=sum(NB_VALD), X=max(X), Y=max(Y))
Manip_dist_Jour = function(DataFrame) {
  DataFrame$JOURSEM=factor(DataFrame$JOURSEM,levels = c("lundi","mardi","mercredi","jeudi","vendredi","samedi","dimanche"))
  DataFrame$TypeJ=character(nrow(DataFrame))
  DataFrame$TypeJ[DataFrame$JOUR %in% Vacances$Ete]="Ete"
  DataFrame$TypeJ[DataFrame$JOUR %in% Vacances$Noel]="Noel"
  DataFrame$TypeJ[DataFrame$JOUR %in% Vacances$Ski]="Ski"
  DataFrame$TypeJ[DataFrame$JOUR %in% Vacances$Printemps]="Printemps"
  DataFrame$TypeJ[DataFrame$JOUR %in% Vacances$Toussaint]="Toussaint"
  DataFrame$TypeJ[DataFrame$JOUR %in% Vacances$All == FALSE]="HorsVacances"
  DataFrame$CAT_JOUR=character(nrow(DataFrame))
  DFr=list()
  ii=(DataFrame$JOURSEM!="samedi" & DataFrame$JOURSEM!="dimanche") & DataFrame$TypeJ!="HorsVacances"
  DataFrame$CAT_JOUR[ii]="JOVS"
  DFr$JOVS$Data = DataFrame[ii,]
  DFr$JOVS$Nom="Jours ouvrés Vacances Scolaires"
  ii=(DataFrame$JOURSEM!="samedi" & DataFrame$JOURSEM!="dimanche") & DataFrame$TypeJ=="HorsVacances"
  DataFrame$CAT_JOUR[ii]="JOHV"
  DFr$JOHV$Data = DataFrame[ii,]
  DFr$JOHV$Nom="Jours ouvés Hors Vacances Scolaires"
  ii=DataFrame$JOURSEM=="samedi" & DataFrame$TypeJ!="HorsVacances"
  DataFrame$CAT_JOUR[ii]="SAVS"
  DFr$SAVS$Data = DataFrame[ii,]
  DFr$SAVS$Nom="Samedi VS"
  ii=DataFrame$JOURSEM=="samedi" & DataFrame$TypeJ=="HorsVacances"
  DataFrame$CAT_JOUR[ii]="SAHV"
  DFr$SAHV$Data = DataFrame[ii,]
  DFr$SAHV$Nom="Samedi HV"
  ii=DataFrame$JOURSEM=="dimanche"
  DataFrame$CAT_JOUR[ii]="DIJFP"
  DFr$DIJFP$Data = DataFrame[ii,]
  DFr$DIJFP$Nom="Dimanche"
  return(list("TypeJ"=DFr,"Distr"=DataFrame))
}
res=Manip_dist_Jour(distr_JOUR$Y)
distr_TypeJ=res$TypeJ
distr_JOUR$Y=res$Distr
res=Manip_dist_Jour(distr_JOURV$Y)
distr_TypeJV=res$TypeJ
distr_TypeJ_Station=list()
res=Manip_dist_Jour(distr_JOUR_Station$S1)
distr_TypeJ_Station$S1=res$TypeJ
distr_JOUR_Station$S1=res$Distr
res=Manip_dist_Jour(distr_JOUR_Station$S2)
distr_TypeJ_Station$S2=res$TypeJ
distr_JOUR_Station$S2=res$Distr
rm(res)

On peut alors tracer toutes sortes de graphiques, par exemple le nombre de validations, par jour, entre le 1er janvier et le 31 décembre, en fonction du jour de la semaine.

g0 = ggplot(distr_JOUR$Y, aes(x=JOUR, y=NB_VALD, color = JOURSEM)) + geom_point()
g0 = g0 + labs(title="Nombres de validations chaque jours de 2017", x="Date", y="Nombre de validations")
g0

On peut voir la très forte baisse les jours de semaine pendant les vacances d’été. Au lieu de regarder sur l’année, on peut regarder sur la journée

Fct_FqH = function(DataFrame,distr_pers) {
DataFrame=dplyr::full_join(DataFrame,distr_pers[,c("CAT_JOUR","CODE_STIF_ARRET","pourc_validations","TRANCHE","TRNC_HORR_60")],by=c("CODE_STIF_ARRET"="CODE_STIF_ARRET","CAT_JOUR"="CAT_JOUR"))
  DataFrame$NB_VALD=DataFrame$NB_VALD*DataFrame$pourc_validations
  return(DataFrame)
}
distr_JOUR_Station$S1=Fct_FqH(distr_JOUR_Station$S1, distr_pers$S1)
distr_JOUR_Station$S2=Fct_FqH(distr_JOUR_Station$S2, distr_pers$S2)
distr_JOUR_Station$Y=rbind(distr_JOUR_Station$S1,distr_JOUR_Station$S2)
distr_JOUR_Station$Y=distr_JOUR_Station$Y[is.na(distr_JOUR_Station$Y$NB_VALD)==FALSE,]

On peut alors faire un graphique, en fonction de la tranche horaire, pour certaines périodes, par exemple en dehors de vacances scolaires, en semaine (par heure, on a ici un boxplot)

Graphique_HOR = function(DataFrame,TypeJ,NomJ) {
  # Graphique de la distribution de l'affluence par tranche horaire et type de jours
  g1 = ggplot(DataFrame[DataFrame$CAT_JOUR==TypeJ,], aes(x=TRNC_HORR_60, y=pourc_validations, color = TRNC_HORR_60,las=2)) + geom_boxplot() + ylim(c(0,100))
  g1 = g1 + labs(title=paste(c("Distribution des validations par tranche horaire ",NomJ), sep="", collapse = ""), x="Jours", y="Nombre de validations") +
  theme(axis.text.x= element_text(size = 8, angle = 45))
  g1
}
Graphique_HOR(distr_JOUR_Station$Y,"JOHV","Jours ouvrés Hors Vacances Scolaires")

ou bien le samedi

Graphique_HOR(distr_JOUR_Station$Y,"SAHV","Samedi Hors Vacances Scolaires")

On peut tenter un peu de cartographie. Comme nombre de métros/bus, dans le monde, on a souvent uniquement accès aux nœuds d’entrée dans le réseau (et pas aux nœuds de sortie). Mais ça reste intéressant, et très informatif

get_Paris1 = get_map(c(2.3448688,48.8613029), zoom = 11)
Paris1 = ggmap(get_Paris1)

Par gare, et par heure, on peut regarder le nombre de validations de tickets

Median_Valid = distr_JOUR_Station$Y %>% group_by(CAT_JOUR, LIBELLE_ARRET, X, Y) %>% summarise(NB_VALD=median(NB_VALD))
Median_Valid_Station = distr_JOUR_Station$Y %>% group_by(CAT_JOUR, TRNC_HORR_60,LIBELLE_ARRET, X, Y) %>% summarise(NB_VALD=median(NB_VALD))
 
Carte_Densite = function(Nom,Carte,TypeJ,HOR,DataFrame) {
if (HOR=="") {
    ii=DataFrame$CAT_JOUR==TypeJ
    NomSave=paste("Densité des validations",Nom,TypeJ)
  }
  else {
    ii=DataFrame$CAT_JOUR==TypeJ & DataFrame$TRNC_HORR_60==HOR
    NomSave=paste("Densité des validations",Nom,TypeJ,HOR)
  }
  U=DataFrame[ii,]
  n=round(log10(median(U$NB_VALD)))-1
  n=max(1,10^n)
  Nb_Repete_Stations=ceiling(U$NB_VALD/n)
  U$Size_Stations=U$NB_VALD/max(U$NB_VALD)
  Z=U[rep(1:nrow(U),Nb_Repete_Stations),]
  Carte_A= Carte + geom_point(aes(x=X,y=Y),data=Z,col="coral", size=10*Z$Size_Stations) +
    geom_density2d(data = Z, aes(x=X,y=Y), size = 0.5) + 
    stat_density2d(data = Z, aes(x=X,y=Y,fill = ..level.., alpha = ..level..),size = 0.01, bins = 16, geom = "polygon") +
    scale_fill_gradient(low = "chartreuse", high = "red",guide = FALSE) + 
    scale_alpha(range = c(0, 0.3), guide = FALSE) + ggtitle(NomSave) +
    theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x= element_blank(), axis.text.y = element_blank())
 
  suppressWarnings(print(Carte_A))
}

Par exemple, si on regarde les points de validations de tickets entre 5 et 6 heures du matin, on obtient

L=levels(Median_Valid_Station$TRNC_HORR_60)
Carte_Densite("dans la petite ceinture",Paris1,"JOHV",L[6],Median_Valid_Station)

avec beaucoup de ville dans la banlieue proche. Plus tard en journée, entre 11 heures et midi, les gares de validation sont davantage dans le cœur de Paris, avec la Défense a gauche et Saint-Denis au nord

Carte_Densite("dans la petite ceinture",Paris1,"JOHV",L[12],Median_Valid_Station)

En fin de journée, c’est Paris et surtout la Défense qui ressortent

Carte_Densite("dans la petite ceinture",Paris1,"JOHV",L[19],Median_Valid_Station)

Amusant, non ?

Modélisation “appliquée”

Depuis la rentrée, je donne un cours de modélisation appliquée, ou pour être plus précis, intitulé “modèles linéaires appliqués“. Ça a été l’occasion de me rendre compte que le mot “appliqué” pouvait prêter à confusion. Dans ma tête, on allait appliquer une théorie (la théorie des modèles linéaires) à cas concrets. Généralement, je présente les grandes lignes de la théorie : les hypothèses, les principaux résultats, et je fais des démonstrations “avec les mains” pour expliquer de manière intuitive certains passages (genre la perte des degrés de liberté dans le test de Student, le signe de la covariance des estimateurs dans le modèle simple, etc) en renvoyant à des livres de référence pour les preuves complètes. Et j’essaye, autant que faire se peut, d’illustrer sur des données, avec du code R (par exemple vendredi dernier, c’était sur la notion de “modalité de référence” quand on régresse sur une variable catégorielle, montrant que tous les modèles sont équivalent et que ce n’est vraiment qu’une histoire de conventions). J’ai aussi donné un gros devoir, avec des vraies données, 2400 observations, 34 variables (quantitatives mais aussi qualitatives). Bref, un vrai cas d’étude. Et j’ai même proposé de revenir sur leurs modélisations en fin de cours…

Et j’ai été surpris quand un groupe d’élèves est venu me voir à la pause, pour me demander s’il y aurait des “chiffres” : si au lieu d’écrire {\displaystyle {\hat {\beta }}_{1}={\frac {\sum x_{i}\sum y_{i}-n\sum x_{i}y_{i}}{\left(\sum x_{i}\right)^{2}-n\sum x_{i}^{2}}}={\frac {\sum (x_{i}-{\bar {x}})(y_{i}-{\bar {y}})}{\sum (x_{i}-{\bar {x}})^{2}}}}ou encore {\displaystyle {\hat {\beta }}_{0}={\frac {\sum y_{i}-{\hat {\beta }}_{1}\sum x_{i}}{n}}={\bar {y}}-{\hat {\beta }}_{1}{\bar {x}}}il y aurait des calculs, à la calculatrice, avec des chiffres au lieu des x et des y. J’ai été très supris, parce qu’on est en 2018, que tout le monde parle de “big data” et revenir à la calculatrice pour calculer des sommes de 10 nombres, c’est possible. En tous cas, certains profs le font encore. J’ai dit que si ça pouvait les rassurer, je pouvais chercher des bases dans des vieux livres, mais que c’était pas ma manière de comprendre le terme “appliqué” dans le nom du cours. Je me suis souvenu du livre de Robert Pindyck et Daniel Rubinfeld, Econometric Models and Economic Forecasts, qui était le livre de référence de la SOA il y a 20 ans, qui avait des tableaux de ce genre. Mais je ne l’ai pas trouvé en ligne. Par contre, je suis tombé sur la Table 4.3. dans le livre de John Neter, Applied Linear Regression Models, paru en 1974 (repris dans la Table 3.4. de Applied Linear Statistical Models)

On a ici juste 11 observations, un modèle linéaire simple. Il est possible de faire tous les calculs à la main. Mais je me suis dit que ce n’était pas un cadeau : on ne peut pas faire grand chose sans faire un minimum de maths… Faisons juste une régression linéaire sur ces données, pour “voir”

base = data.frame(x=c(125,100,200,75,150,175,75,175,125,200,100),
                  y=c(160,112,124,28,152,156,42,124,150,104,136))
plot(base,xlab="Size of Minumum Deposit",
     ylab="Number of New Accounts",pch=19,ylim=c(30,180))
abline(lm(y~x,data=base),col="red")

La droite de régression (qui minimise la somme des carrés des erreurs) est la courbe (ou ici “droite”, on est sur un modèle linéaire) rouge

summary(lm(y~x,data=base))
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  50.7225    39.3979   1.287     0.23
x             0.4867     0.2747   1.772     0.11

On ne voit rien, littéralement ! Surtout que sur cet exemple, la question qui nous brule les lèvres, c’est “où est le maximum ?” car on voit bien qu’on a une courbe qui semble concave, croissante au début, décroissante sur la fin… On peut faire du lissage non-paramétrique pour s’en convaincre

with(base, 
     scatter.smooth(x, y,
        lpars = list(col = "red"),
        xlab="Size of Minumum Deposit",
       ylab="Number of New Accounts",
       ylim=c(30,180),pch=19))

Peut-aller plus loin, pour répondre à la question “où est le maximum ?”. En particulier, peut-on suggérer une valeur, voir trouver un intervalle de confiance ?

On peut continuer avec un modèle paramétrique – l’essentiel de ce qu’on voit dans les cours de statistique tourne autour des modèles paramétrique. Mais si on est plus pragmatique, on peut considérer un modèle quadratique, autrement dit, notre prévision sera une parabole. Ça pourrait être cohérent avec ce qu’on voit quand on trace les points, et puis souvent, les paraboles ont un maximum

with(base, 
     scatter.smooth(x, y,
        lpars = list(col = "red",lty=2),
        xlab="Size of Minumum Deposit",
        ylab="Number of New Accounts",
        ylim=c(30,180),pch=19))
reg=lm(y~x+I(x^2),data=base)
vx=seq(min(base$x),max(base$x),length=251)
vy=predict(reg,newdata=data.frame(x=vx))
lines(vx,vy,col="red")

On peut voir que le modèle pourrait sembler réaliste, non seulement visuellement, mais en plus si on regarde la sortie de la régression (même si on pourrait discuter ici l’interprétation de la significativité…)

summary(reg)
 
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.255e+02  5.589e+01  -5.824 0.000394 ***
x            6.569e+00  8.744e-01   7.513 6.84e-05 ***
I(x^2)      -2.203e-02  3.143e-03  -7.011 0.000111 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Or on peut montrer que pour une parabole de la forme y=\beta_2x^2+\beta_1x+\beta_0 l’optimum (qui sera un maximum si \beta_2<0 – la courbe part en -\infty à gauche et à droite) sera obtenu en x^\star=\theta=-\beta_1/2\beta_2. Un estimateur naturel de \theta est alors \widehat{\theta}=-\widehat{\beta}_1/2\widehat{\beta}_2 où les estimateurs sont obtenus ici par minimisation de la somme des carrés des erreurs. Mais comment avoir la variance de cet estimateur ? On peut naturellement tenter la \Delta-méthode (décrite sur wikipedia), en considérant \displaystyle{\theta=g(\beta_1,\beta_2)=\frac{-\beta_1}{2\beta_2}}Comme \displaystyle{\frac{\partial \theta}{\partial \beta_1}=\frac{-1}{2\beta_2}}et \displaystyle{\frac{\partial \theta}{\partial \beta_2}=\frac{\beta_1}{2\beta_2^2}}la variance (asymptotique) vaut ici \left[ \frac{-1}{2\beta_2}~~\frac{\beta_1}{2\beta_2^2} \right] \left[ \begin{array}{cc} \sigma_{1}^2 & \sigma_{12} \\ \sigma_{12} & \sigma_{2}^2 \end{array} \right] \left[\frac{-1}{2\beta_2}~~\frac{\beta_1}{2\beta_2^2} \right]^{\text{ T}}soit \frac{\sigma_1^2\beta_2^2-2\beta_1\beta_2\sigma_{12}+\beta_1^2\sigma_2^2}{4\beta_2^2}On obtient alors un estimateur de cette variance asymptotique en remplaçant les valeurs inconnues par leurs estimateurs. Ici, on aurait, numériquement

beta=coefficients(reg)
n=nrow(base)
sigma=n*vcov(reg)[2:3,2:3]
dg = c(-.5/beta[3],.5*beta[2]/beta[3]^2)
theta=-beta[2]/(2*beta[3])
theta
[1]  149.0676 
s2=t(dg) %*% sigma %*% dg
s2
         [,1]
[1,] 94.59911
sqrt(s2)
         [,1]
[1,] 9.726207

Autrement dit, si on suppose la normalité de notre estimateur (ce que nous dit d’ailleurs la \Delta-méthode, tout du moins asymptotiquement), on a l’intervalle de confiance suivant

theta+qt(c(.025,.975),n-3)*sqrt(s2)
[1] 126.6389 171.4963
vx=theta
vy=predict(reg,newdata=data.frame(x=vx))
points(vx,vy,pch=19,cex=1.5,col="red")
arrows(vx-qt(.975,n-3)*sqrt(s2),vy,
       vx+qt(.975,n-3)*sqrt(s2),vy,code=3,angle=90,length=.1,
       col="red",lwd=1.2)

On peut aussi tenter une autre stratégie, basée sur l’utilisation de la vraisemblance profilée. Je renvoie à mon précédent billet pour expliquer l’idée de la méthode, malheureusement pas assez enseignée dans les cours de statistique… Pour notre modèle quadratique, on maximisait la log-vraisemblance, classiquement, sous une hypothèse de résidus Gaussiens et homoscédastiques

logL = function(pm){
  b0=pm[1]
  b1=pm[2]
  b2=pm[3]
  b3=pm[4]
  -sum(log(dnorm( (base$y-(b0+b1*base$x+b2*base$x^2)) ,0,b3)))
}

Ici, la première astuce est d’introduire \theta (là où se trouve le maximum de la parabole) comme un des paramètres du modèle – par exemple à la place de \beta_2

logL = function(pm){
  b0=pm[1]
  b1=pm[2]
  theta=pm[3]
  b3=pm[4]
  -sum(log(dnorm( (base$y-(b0+b1*base$x-.5*b1/theta*base$x^2)) ,0,b3)))
}

Si on cherche juste le maximum de la vraisemblance, on obtient

 optim(par=c(-323,6.5,140,5),logL)
$par
[1] -325.521244    6.569328  149.070597   13.693488
 
$value
[1] 44.39629

ce qui est cohérent avec notre calcul précédant (et c’est normal, c’est la propriété du cours de statistique qui dit qu’on peut reparamétrer le modèle, ça ne change rien). La seconde astuce est celle de la vraisemblance profilée : on dit que dans notre paramètre multivarié, un est plus important que les autres. Et les autres sont juste maximisés, à \theta donné. Techniquement, la log-vraisemblance profilée du paramètre \theta est

logL=function(theta){
  f=function(pm){
  b0=pm[1]
  b1=pm[2]
  b3=pm[3]
  -sum(log(dnorm( (base$y-(b0+b1*base$x-.5*b1/theta*base$x^2)) ,0,b3)))
  }
  optim(par=c(-323,6.5,5),f)$value
}

On peut d’ailleurs la tracer

v1 = seq(130,190)
v2 = Vectorize(logL)(v1)
plot(v1,-v2,type="l",xlim=range(base$x),xlab="",ylab="log vraisemblance")

Le maximum est atteint, ici en

opt
$minimum
[1] 149.068

Ce qui est cohérent avec nos calculs. On peut ensuite un résultat sur le test du rapport de vraisemblance, qui nous dit que la différence entre la log-vraisemblance au maximum et en la vraie valeur suit la moitié du quantile d’une loi du chi-deux à 1 degré de liberté

points(opt$minimum,-opt$objective,pch=19,cex=1.5,col="red")
ref=-opt$objective-.5*qchisq(.95,1)
abline(h=ref,lty=2,col="red")

On peut alors alors tracer l’intervalle de confiance sur le graphique initial

with(base, 
     scatter.smooth(x, y,
                    lpars = list(col = "red",lty=2),
                    xlab="Size of Minumum Deposit",
                    ylab="Number of New Accounts",
                    ylim=c(30,180),pch=19))
reg=lm(y~x+I(x^2),data=base)
vx=seq(min(base$x),max(base$x),length=251)
vy=predict(reg,newdata=data.frame(x=vx))
lines(vx,vy,col="red")
vx=opt$minimum
vy=predict(reg,newdata=data.frame(x=vx))
points(vx,vy,pch=19,cex=1.5,col="red")
arrows(min(v1[id]),vy,max(v1[id]),vy,code=3,angle=90,length=.1,
       col="red",lwd=1.2)

Voilà au moins deux méthodes, basées sur deux résultats du cours de statistique. Après, comme la majorité des techniques de statistique, ce sont des résultats asymptotiques, et rien ne garantie leur validité avec seulement 11 observations (“rien” est un peu fort… on pourrait faire des développements sur ce point).

Une autre solution est d’utiliser des simulations. On suppose que les observations, c’est un modèle, et du bruit. On peut prendre comme modèle le modèle non-paramétrique (lissage local) et supposer le bruit Gaussien. Pour générer d’autres échantillons, on va garder nos observations en x, par contre, pour y, on va utiliser \widehat{y}+\varepsilon\varepsilon sera tiré suivant une loi normale

reg=loess.smooth(x = base$x, y= base$y, evaluation = 501)
pred=function(x) reg$y[which.min(abs(x-reg$x))]
ypred=Vectorize(pred)(base$x)
epsilon = base$y-ypred
simu = function(graph=FALSE){
newbase = data.frame(x= base$x,
                     y= ypred + rnorm(n,mean(epsilon),sd(epsilon)))
reg=loess.smooth(x = newbase$x, y= newbase$y, evaluation = 501)
lines(reg$x,reg$y,col=rgb(1,0,0,.2))
reg$x[which.max(reg$y)]}
for(i in 1:20) simu(TRUE)
lines(loess.smooth(x = base$x, y= base$y, evaluation = 501),lwd=1.2,col="red")

On récupère ainsi plein de modèles, estimés sur ces nouvelles données. Compte tenu de l’asymétrie des données, on va oublier le modèle quadratique, et utiliser un modèle non-paramétrique, là encore. Et on calcule numériquement le maximum. Et on répète 10,000 fois.

V=Vectorize(simu)(1:1e5)
hist(V,probability = TRUE, col="light blue", border="white",ylim=c(0,.07))
lines(density(V),col="blue")
quantile(V,c(.025,.975))

On a ici la distribution empirique du maximum, observée sur nos 10,000 échantillons simulés. On peut même obtenir un intervalle de confiance, en prenant les quantiles empiriques

vy=pred(mean(V))
arrows(quantile(V,.025),vy,quantile(V,.975),vy,code=3,angle=90,length=.1,
       col="red",lwd=1.2)

Bref, j’ai l’impression que pour faire des applications, il vaut mieux bien connaître la théorie… non ?

Analyse des résultats au baccalauréat des séries générales

Pour continuer sur les manipulation de données publiques, je voulais m’inspirer du projet de Cédric, de la formation Data Science pour l’Actuariat sur les résultats au baccalauréat. Les données nécessaires à cette étude sont disponibles sur plusieurs sites,

Il ne s’agit aucunement d’une analyse poussée des résultats, juste un peu de visualisation, sans aucune autre prétention ! Ah oui, même si on ne va pas faire de carte (je les trouve peu lisibles) on va quand même utiliser les données spatiales : les établissements scolaires sont géolocalises, et on peut obtenir des informations locales, sur le taux de chômage, ou le revenu médian. Et faire des graphiques.

Ce préambule passé, on peut commencer.

library(dplyr)
library(readxl)
library(sp)
library(ggmap)
library(raster)
library(leaflet)
library(DT)
library(cowplot)
library(gstat)
library(tmap)

On va commencer par récupérer par établissement, les résultats au bac.

url_resultat_etab = "https://data.education.gouv.fr/explore/dataset/fr-en-indicateurs-de-resultat-des-lycees-denseignement-general-et-technologique/download/?format=csv&amp;timezone=Europe/Berlin&amp;use_labels_for_header=true"
download.file(url_resultat_etab,destfile = paste0(librairie,"import_resultat_etab.csv"), method="curl")
df_resultat_etab = read.csv("import_resultat_etab.csv",header=TRUE, sep= ";", encoding="UTF-8")

Comme bien souvent avec les données des administrations françaises, on a souvent des soucis de typographie. Pour simplifier, on va supprimer les accents, et uniformiser un peu les noms

MiseEnForme_Colonnes = function(text) {
  text &lt;- gsub("è", "e", text)  
  text &lt;- gsub("é", "e", text)         
  text &lt;- gsub("_", ".", text)
  text &lt;- gsub("serie.", "", text)
  text &lt;- gsub("Effectif.Presents.", "Effectif.", text)
  text &lt;- gsub("Taux.","Tx.",text)
  text &lt;- gsub("Brut.de.Reussite.", "Admis.Etab.", text)
  text &lt;- gsub("Reussite.Attendu.", "Admis.", text)
  text &lt;- gsub("brut", "Etab", text)
  text &lt;- gsub("attendu", "Academie", text)
  text &lt;- gsub("toutes.", "TOTAL", text)
  text &lt;- gsub("Total.", "TOTAL", text)
  text &lt;- gsub("..Etablissement", ".Etab", text)
  text &lt;- gsub("Pourcentage", "Tx", text)
  return(text)
}
for(i in 1:ncol(df_resultat_etab)){
  colnames(df_resultat_etab)[i] &lt;- MiseEnForme_Colonnes(names(df_resultat_etab)[i])
}

On va ensuite supprimer les départements et régions d’outre-mer,

df_resultat_etab = df_resultat_etab[-which(toupper(df_resultat_etab$Departement) %in% c("GUADELOUPE","MAYOTTE","MARTINIQUE","REUNION","GUYANE")),]

récupérer les noms des colonnes

Colonnes = colnames(df_resultat_etab)

et comme on s’intéresse aux premières variables

Colonnes_Generiques = Colonnes[1:8]

on les recupere, pour construire quelques statistiques pour colonnes relatives aux séries L, ES et S

Colonnes_Series = Colonnes[grepl("([a-zA-Z]*?.)*\\.S$|([a-zA-Z]*?.)*\\.ES$|([a-zA-Z]*?.)*\\.L$|([a-zA-Z]*?.)*\\.TOTAL$",Colonnes)]

Et on finit avec les autres

Colonnes_Autres = Colonnes[grepl("(Tx.Bacheliers.*)|(Tx.acces.*)|(Effectif.de.*)|(libelle.region)|(code.region)|(element)",Colonnes)] 
df_resultat_etab = cbind(df_resultat_etab[Colonnes_Generiques],df_resultat_etab[Colonnes_Series],df_resultat_etab[Colonnes_Autres])

On peut aussi localiser les établissements

url_carto_etab &lt;- "https://www.data.gouv.fr/s/resources/adresse-et-geolocalisation-des-etablissements-denseignement-du-premier-et-second-degres/20160526-143453/DEPP-etab-1D2D.csv"
 
download.file(url_carto_etab,destfile=paste0(librairie,"import_carto_etab.csv"))
df_carto_etab = read.csv2("import_carto_etab.csv",header=TRUE

On récupère ici la géolocalisation de 66556 établissements nationaux ! On peut croiser avec des données socio-économiques des communes

nom_base_emploi = "base-cc-emploi-pop-act-2014"
url_baseemploi_popactive = paste0("https://www.insee.fr/fr/statistiques/fichier/2862207/",nom_base_emploi,".zip")
download.file(url_baseemploi_popactive,destfile=paste0(librairie,nom_base_emploi,".zip"))
unzip(paste0(nom_base_emploi,".zip"),overwrite = TRUE) 
df_base_emploi_source = read_excel(paste0(nom_base_emploi,".xls"),sheet="COM_2014",skip=5)

On va exclure les territoires d’outre-mer ici

df_base_emploi_source &lt;- df_base_emploi_source[-which(df_base_emploi_source$DEP %in% c("971","972","973","974","975")),]
df_base_emploi_colonnes = c("CODGEO","P14_POP1564","P14_H1564","P14_F1564","P14_ACT1564","P14_ACTOCC1564","P14_CHOM1564","P14_INACT1564", "P14_ETUD1564", "P14_RETR1564", "P14_AINACT1564", "P14_HCHOM1524", "P14_FCHOM1524", "C14_ACT1564","C14_ACT1564_CS1","C14_ACT1564_CS2","C14_ACT1564_CS3","C14_ACT1564_CS4","C14_ACT1564_CS5","C14_ACT1564_CS6","P14_POP15P")
df_base_emploi = df_base_emploi_source[,names(df_base_emploi_source) %in% df_base_emploi_colonnes]

et corriger les soucis classiques de la Corse,

MiseEnForme_CodeGeo = function(text) {
  text &lt;- gsub("2A", "20", text)  
  text &lt;- gsub("2B", "20", text)  
}
df_base_emploi$CODGEO = MiseEnForme_CodeGeo(df_base_emploi$CODGEO)

On peut aussi utiliser des données de revenus, par communes

nom_base_revenus = "indic-struct-distrib-revenu-2014-COMMUNES"
url_baserevenus = paste0("https://www.insee.fr/fr/statistiques/fichier/3126151/",nom_base_revenus,".zip")
download.file(url_baserevenus,destfile=paste0(librairie,nom_base_revenus,".zip"))
unzip(paste0(nom_base_revenus,".zip"),overwrite = TRUE)
df_base_revenus = read_excel("FILO_DISP_COM.xls",sheet="ENSEMBLE",skip=5)[,c(1,4,7)]
df_base_revenus$CODGEO = MiseEnForme_CodeGeo(df_base_revenus$CODGEO)

On recupere des donnees spatiales relatives aux communes

url_geoloc_communes = "http://www.nosdonnees.fr/wiki/images/b/b5/EUCircos_Regions_departements_circonscriptions_communes_gps.csv.gz"
download.file(url_geoloc_communes,destfile=paste0(librairie,"geoloc_communes.csv.gz"))
df_geoloc_communes = read.csv2(gzfile("geoloc_communes.csv.gz"),header=TRUE, stringsAsFactors = FALSE,encoding="UTF-8")

et comme toujours, un peu de corrections s’imposent

df_geoloc_communes = df_geoloc_communes[-which(df_geoloc_communes$numéro_département %in% c("971","972","973","974","975")),]
df_geoloc_communes = df_geoloc_communes[,names(df_geoloc_communes) %in% c("code_insee","latitude","longitude","codes_postaux")]
df_geoloc_communes_nb &lt;- nrow(df_geoloc_communes)

On va ensuite creer une fonction de remplacement des valeurs manquantes, et de correction des séparateurs décimaux

MiseEnForme_CoordonneesGeo = function(valeur){
pretraitement = ifelse(as.character(valeur)=="-","0",as.character(valeur))
traitement = as.numeric(ifelse(pretraitement==".","0",gsub(pattern=",",replacement=".",pretraitement)))
  return(traitement)
}
df_geoloc_communes$latitude = MiseEnForme_CoordonneesGeo(df_geoloc_communes$latitude)
df_geoloc_communes$longitude = MiseEnForme_CoordonneesGeo(df_geoloc_communes$longitude)

On passe ensuite a l’élimination des lignes en double

df_geoloc_communes = unique(df_geoloc_communes)

On va ensuite changer les noms des colonnes pour harmoniser avec les autres bases

names(df_geoloc_communes) = c("Codes_Postaux","CODGEO","coordonnee_y","coordonnee_x")

On peut ensuite rechercher les lignes en double sur les codes insee

liste_CODGEO2 = aggregate(x=df_geoloc_communes$Codes_Postaux,by=list(df_geoloc_communes$CODGEO),FUN="length")
list_geoloc_communes_CODGEO2 = liste_CODGEO2[liste_CODGEO2$x&gt;1,1]
df_geoloc_communes_CODGEO2 = df_geoloc_communes[df_geoloc_communes$CODGEO %in% list_geoloc_communes_CODGEO2,1:2]

Ici, un correction manuelle s’impose pour 4 configurations : les données propres à Lyon, Paris et Marseille ne sont pas géolocalisées, les données propres à la ville de Laguépie sont géolocalisées en doubles

df_geoloc_communes_propre = df_geoloc_communes[!df_geoloc_communes$CODGEO %in% list_geoloc_communes_CODGEO2,]
df_geoloc_communes_corrige = data.frame(Codes_Postaux=c("13001","69001","75001","82250"),                        CODGEO=c("13055","69123","75056","82088"),
coordonnee_y=c(43.3,45.75,48.85,44.15),
coordonnee_x=c(5.4,4.85,2.31,1.97))
df_geoloc_communes = rbind(df_geoloc_communes_propre,df_geoloc_communes_corrige)

On peut enfin fusionner les bases

df_etab = merge(df_resultat_etab,df_carto_etab,by="Cod.Etab")

Certains établissements ne peuvent être géolocalisés pour certaines années

df_etab_total_nongeolocalises &lt;- df_resultat_etab[!df_resultat_etab$Cod.Etab %in% df_carto_etab$Cod.Etab,]

Comme l’étude ne porte que sur les seuls lycées d’enseignement général et technologique, le dataframe est réduit aux observations relatives d’une part aux lycées, d’autre part aux établissements d’enseignement polyvalent, général ou général et technologique.

df_etab = df_etab[grep("LYCÉE",toupper(df_etab$nature_uai_libe)),]
df_etab = df_etab[grep("GÉNÉRAL|POLYVALENT",toupper(df_etab$nature_uai_libe)),]
df_etab_nongeolocalises = df_etab[df_etab$Cod.Etab %in% df_etab_total_nongeolocalises$Cod.Etab,]
df_etab_geolocalise = df_etab[!is.na(df_etab$coordonnee_x),]
df_etab_geolocalise = df_etab_geolocalise[!is.na(df_etab_geolocalise$coordonnee_y),]

Enfin, on va convertir les code géographiques (ici reconnus comme facteur) en chaines de caractères (de 5 caractères) pour pouvoir fusionner les tables

ConvertCODGEO = function(code) {
  if(is.character(code)) {
    code_character = ifelse(nchar(code)&lt;5, paste0("0",code), code)
    return(code_character)
  }
  else if(is.factor(code)){
    code_character = ifelse(code&lt;10000, paste0("0",as.numeric(as.character(code))), as.numeric(as.character(code)))
    return(code_character)
  }
  else if(is.numeric(code)){
    code_character = ifelse(code&lt;10000, paste0("0",code), as.character(code))
    return(code_character)
  } 
}
df_etab_geolocalise$Code.commune = ConvertCODGEO(df_etab_geolocalise$Code.commune) 
df_etab_geolocalise$Secteur.Public.Prive = sapply(df_etab_geolocalise$Secteur.Public.Prive,function(nature) {ifelse(nature=="PU","Lycées Publics","Lycées Privés")})

On conservation alors les établissements dont la commune n’est pas manquante

df_etab_geolocalise = df_etab_geolocalise[!is.na(df_etab_geolocalise$Code.commune),]

Pour finir, on va creer une base, pour ensuite faire une graphique

tbl_etab_nature_res_source = df_etab_geolocalise[,c(3,8,9,10,11,13,14,15)]
for(i in c(6,7,8)){
  temp = tbl_etab_nature_res_source[!is.na(tbl_etab_nature_res_source[i]),c(1,2,i-3,i)]
  temp$Serie = ifelse(i==6,"L",ifelse(i==7,"ES","S"))
  names(temp)[2:4] = c("Nature","Effectif","Tx.Admis")
  if(i==6){
    tbl_etab_nature_result = temp
  }
  else{
    tbl_etab_nature_result = rbind(tbl_etab_nature_result,temp)
  }
}
graph = ggplot(tbl_etab_nature_result,aes(x=Effectif,y=Tx.Admis,colour=factor(Annee))) 
graph = graph + geom_point(alpha=0.45)
graph = graph + facet_grid(Serie~Nature)
graph = graph + xlab("Effectifs de l'établissement en terminale (par série)") + ylab("Taux d'admission (%)") 
graph = graph + scale_color_discrete(name="Année des\nrésultats")
graph = graph + theme(legend.title = element_text(size=9,face="bold"),
       legend.text = element_text(size=9),
       strip.background = element_rect(colour="black", fill="gray95"),
       panel.border = element_rect(linetype = "solid"),
       panel.grid.major = element_line(colour = "gray75",linetype = "dashed"),
       panel.grid.minor = element_line(colour = "gray95",linetype = "dashed"),
       axis.title.x = element_text(size=9, face="bold"),
       axis.text.x  = element_text(size=8),
       axis.title.y = element_text(size=9, face="bold"),
       axis.text.y  = element_text(size=8))
graph

On a ici l’evolution des resultats en fonction de la taille des etablissements.

df_communes_CorrNaN = df_communes_Corr[which(!df_communes_Corr$TxChomage == "NaN" &amp; !df_communes_Corr$TxCadres == "NaN" &amp; !df_communes_Corr$TxOuvriers == "NaN" &amp; !df_communes_Corr$NbPopulation == "NaN" &amp; !df_communes_Corr$TxSenior == "NaN" &amp; !df_communes_Corr$RevenusMedians == "NaN"),]
df_communes_sp = SpatialPointsDataFrame(coords = df_communes_CorrNaN[, c("coordonnee_x", "coordonnee_y")], data = df_communes_CorrNaN) 
Grille              = as.data.frame(makegrid(df_communes_sp, nsig=2, cellsize = 0.1))
names(Grille)       = c("X", "Y")
coordinates(Grille) = c("X", "Y")
gridded(Grille)     = TRUE  
fullgrid(Grille)    = TRUE  
proj4string(Grille) = proj4string(df_communes_sp)

On peut ensuite faire du krigeage, histoire de lisser un peu nos donnees de chomage et de revenu

df_communes_sp.TxChomage = krige(TxChomage ~ 1, df_communes_sp, Grille, nmax=1)
df_communes_sp.RevenusMedians = krige(RevenusMedians ~ 1, df_communes_sp, Grille, nmax=1)
sp_lycee_WGS84@data$TxChomage      = extract(R.TxChomage,sp_lycee_WGS84)
sp_lycee_WGS84@data$TxCadres       = extract(R.TxCadres,sp_lycee_WGS84)
sp_lycee_WGS84@data$TxOuvriers     = extract(R.TxOuvriers,sp_lycee_WGS84)
sp_lycee_WGS84@data$NbPopulation   = extract(R.NbPopulation,sp_lycee_WGS84)
sp_lycee_WGS84@data$TxSenior       = extract(R.TxSenior,sp_lycee_WGS84)
sp_lycee_WGS84@data$RevenusMedians = extract(R.RevenusMedians,sp_lycee_WGS84)

On peut enfin conclure, en faisant une fonction generique de visualisation

Creation_Graphique = function(df, AnneeObs_Ouv, AnneeObs_Clo, Effectifs, Abscisses, Ordonnees, TitreAbs, TitreOrd, CouleurGraph, CouleurLiss, Serie) {
  df_temp = df[which(df$Annee&gt;=AnneeObs_Ouv &amp; df$Annee&lt;=AnneeObs_Clo),]
  df_temp = df_temp[which(!is.na(df_temp[,Effectifs]) &amp; !is.na(df_temp[,Abscisses]) &amp; !is.na(df_temp[,Ordonnees])),]
  df_temp = df_temp[!df_temp[,Effectifs]==0,]
  df_temp = df_temp[,c(Abscisses,Ordonnees)]
  graphique = ggplot(df_temp,aes(x = df_temp[,Abscisses],y = df_temp[,Ordonnees])) 
  graphique = graphique + geom_point(data = df_temp, aes(x = df_temp[,Abscisses],y = df_temp[,Ordonnees]),size=1, color=CouleurGraph,alpha=0.25) 
  graphique = graphique + geom_density2d(aes(colour=..level..),show.legend=F) + scale_colour_gradient(low="gray55",high="gray25") 
  graphique = graphique + scale_y_continuous(breaks= seq(80,100,by=2), limits = c(80,100))
  graphique = graphique + xlab(TitreAbs) + ylab(TitreOrd) 
  graphique = graphique + ggtitle(Serie) 
  graphique = graphique + theme(plot.title   = element_text(size=13,color=CouleurLiss, face="bold", hjust=0),
       axis.title.x = element_text(size=8, face="bold"),
       axis.text.x  = element_text(size=8),
       axis.title.y = element_text(size=8, face="bold"),
       axis.text.y  = element_text(size=8),
       panel.border = element_rect(linetype = "solid"),
       panel.grid.major = element_line(colour = "gray55",linetype = "dashed"),
       panel.grid.minor = element_line(colour = "gray75",linetype = "dashed")) 
graphique = graphique + stat_smooth(method = "loess",fill=CouleurLiss,color=CouleurLiss)
  return(graphique)
}
Production_Graphique_VI_1 = function(df, Titre_General, Axe_Abscisses, Titre_Abscisses, Annee_Observee_Ouv, Annee_Observee_Clo){
  Graph_S = Creation_Graphique(df, Annee_Observee_Ouv, Annee_Observee_Clo, "Effectif.S", Axe_Abscisses, "Tx.Admis.Etab.S", Titre_Abscisses, "Taux d'admission (%)", "dodgerblue3","dodgerblue4","Série S")
  Graph_ES = Creation_Graphique(df, Annee_Observee_Ouv, Annee_Observee_Clo, "Effectif.ES", Axe_Abscisses, "Tx.Admis.Etab.ES", Titre_Abscisses, "Taux d'admission (%)","darkorange2","darkorange3","Série ES")
  Graph_L = Creation_Graphique(df, Annee_Observee_Ouv, Annee_Observee_Clo, "Effectif.L", Axe_Abscisses, "Tx.Admis.Etab.L", Titre_Abscisses, "Taux d'admission (%)", "chartreuse4","darkgreen","Série L")
  Graph_TS = Creation_Graphique(df, Annee_Observee_Ouv, Annee_Observee_Clo, "Effectif.Etab", Axe_Abscisses, "Tx.Admis.Etab", Titre_Abscisses, "Taux d'admission (%)", "indianred1","red4","Toutes séries")
  p = plot_grid(Graph_S, Graph_ES, Graph_L, Graph_TS, ncol = 2, nrow = 2,align = 'hv',
  scale = c(0.95, 0.95, 0.95, 0.95),vjust = 0.9, hjust=-0.5)
  titre &lt;- ggdraw() + draw_label(Titre_General,fontface="bold", size=10)
  plot_grid(titre, p, ncol = 1, rel_heights=c(.25,5))
}

On note ici

df_lycee &lt;- sp_lycee_WGS84@data

et on peut faire un premier graphique, avec le taux de chomage

Production_Graphique_VI_1(df              = df_lycee,
                          Titre_General   = "Taux d'admission par série en fonction du taux de chômage \n dans la population active - Tous lycées confondus",
                          Axe_Abscisses   = "TxChomage",
                          Titre_Abscisses = "Taux de chômage dans la population active (%)",
                          Annee_Observee_Ouv     = "2013",
                          Annee_Observee_Clo     = "2015")

On ne va pas enfoncer les portes ouvertes de l’inference ecologique, en affirmant des choses aussi stupides que “on a moins de chances d’avoir le bac quand on est au chômage”. Mais on peut noter que dans les zones avec un fort taux de chômage, les résultats au bac sont moins bons.

On peut ensuite regarder en fonction du revenu de la commune du lycée

Production_Graphique_VI_1(df                     = df_lycee,
                          Titre_General          = "Taux d'admission par série en fonction du niveau des revenus disponibles médians - Tous lycées confondus",
                          Axe_Abscisses          = "RevenusMedians",
                          Titre_Abscisses        = "Quantile du niveau des revenus disponibles médians (%)",
                          Annee_Observee_Ouv     = "2013",
                          Annee_Observee_Clo     = "2015")

Fascinant, non ? Mais c’est clairement juste une première approche… il faudrait aller plus loin ensuite !

Scraper, ou pas ?

Ce matin, je mettais en ligne un billet scraper la base d’incendies de forêts expliquant comment scraper la base Promethee, en remplaçant automatique une formulaire. Dans la soirée, sur Twitter me faisait remarquer que c’était intéressant, mais peut être inutile (sur cet exemple particulier en tous cas).

Sur Firefox, il existe un “Moniteur Réseau“, qui s’ouvre en cliquant sur le clavier CtrlMaj + E ( CommandOption + E quand on est est sur un Mac). Ça ouvre un espace en bas de la fenêtre

Si maintenant je lance a la main ma requête, on voit ce qui est fait

Une methode GET apparait. En cliquant dessus, on a toute l’information nécessaire a droite

Dans l’onglet ‘Reponse’, on voit l’information sur le fichier json créé,

on a plus d’information bien entendu, si on agrandit

et surtout si on va dans l’onglet ‘En-tetes’

on a l’URL complet de la requête ! On peut ensuite copier le lien, et l’ouvrir

On a ainsi a accès a presque toute la table. Quand on regarde l’URL, seuls 20 incendies sont renvoyés, a cause du &nbrLigne=20& mais on peut tenter de le changer. propose de mettre &nbrLigne=10000& et il affirmait que ça marchait. J’avoue avoir essayé, sans succès.

Mais je retiens le tuyau en tous cas. Parfois, on peut faire simple !

Scraper la base d’incendies de forêts

Aller, deux billets pour parler de scraping, ou harvesting serait peut être plus juste ? Disons qu’on va faire un code automatique pour récupérer des informations en ligne. Plus particulièrement, on va s’intéresser au cas ou les données sont obtenues par une requête manuelle, et le site n’a pas d’API pour générer un joli fichier csv. Pierre-Marie, de la formation Data Science pour l’Actuariat a propose un exemple dont je vais m’inspirer (avec une autre methode pour scraper car j’ai eu des soucis pour installer le logiciel).

En 1973 a été lancée en France Prométhée, une base de données sur les incendies de forêts de la région méditerranéenne (couvrant 15 départements du Sud-Est). Oui, Προμηθεύς est ce titan qui a voulu jouer avec le feu… On peut consulter la base sur http://promethee.com/default/incendies,

En remplissant les champs, on peut faire une extraction suivant tel ou tel critère. Par exemple, en 2016, en Corse, il y a eu 6 incendies de nature “accidentelle”,

On le voit, le lien de la page ne change pas… il faut donc trouver une methode pour remplir manuellement les champs, et récupérer les informations. En particulier, on voit qu’on peut avoir accès aux informations relatives aux différents sinistres en cliquant sur la loupe

Idéalement, on voudrait un code R, automatique, pour faire ça ! Une première solutions, proposée par Pierre-Marie, est de passer par RSelenium. Le package R (non maintenu sur le CRAN depuis quelques semaines) est base sur Selenium, “Selenium automates browsers” comme dit la publicité. Le code pour scraper une page comme par

require("devtools")
devtools::install_github("ropensci/RSelenium")
require(RSelenium)
remDr = remoteDriver(remoteServerAddr = "localhost" 
                      , port = 4445L
                      , browserName = "firefox")

Le soucis est que ce code ne tourne pas sur mon mac. Il y a de nombreuses pages qui parlent des spécificités de mac (en particulier le soucis du port par défaut 4445, qu’il faut changer). J’ai tout essayé, sans succès.

Une autre piste est de passer par un autre package

require("devtools")
devtools::install_github("ropensci/wdman")

wdman, pour “Webdriver Manager“, propose aussi de le faire, en s’appuyant sur diverses technologies, selenium, chromedriver, phantomJS binary ou encore geckodriver. Moyennant l’installation de chrome, la seconde option est intéressante

library(wdman)
cDrv = chrome()
eCaps = list(chromeOptions = list(
  args = c('--headless', '--disable-gpu', '--window-size=1280,800')
))
remDr = remoteDriver(browserName = "chrome", port = 4567, 
                     extraCapabilities = eCaps)
remDr$open()

qui marche ! On peut lui demander d’aller sur la base promethee, et de faire une copie d’écran

remDr$navigate("http://www.promethee.com/default/incendies")
remDr$screenshot(display = TRUE)

On va ensuite remplir les champs, pour notre requête, sur la page. Par exemple, pour la date de début, on peut remonter au début des années 80

webElem = remDr$findElement(using = 'name', 'dtAlerteDeb')
webElem$clearElement()
webElem$sendKeysToElement(list('01/01/1981'))

et terminer fin 2017

webElem = remDr$findElement(using = 'name', 'dtAlerteFin')
webElem$clearElement()
webElem$sendKeysToElement(list('31/12/2017'))

On se focalise seulement sur la Corse

webElem = remDr$findElement(using = 'name', 'codeg')
webElem$sendKeysToElement(list('CORSE'))

et on veut 50 observations par page

webElem = remDr$findElement(using = 'name', 'nbrLigne')
webElem$sendKeysToElement(list('50'))

Ah oui, et on va temporiser entre chaque requête

webElem = remDr$findElement(using = 'name', 'btnSubmit')
webElem$submitElement()
Sys.sleep(15)
remDr$screenshot(TRUE)

(c’est une règle assez classique quand on scrap: on évite d’aller trop vite, pour pas se faire repérer…). On va ensuite commencer a regarder ce que l’on demande

text_base = read_html(remDr$getPageSource()[[1]])
list_tr = html_nodes(text_base, "tr")
a = html_nodes(list_tr[length(list_tr)], "a")
page_max = as.numeric(html_text(a[10]))

On a environ 620 pages ! avec au moins 10 secondes pour changer de pages, et en rajoutant 15 secondes (on en parlait avant), ca fait quelques heures pour scraper.

Que veut-on récupérer en grattant la page ? On peut récupérer les coordonnées de l’incendie (latitude/longitude), les caractéristiques de l’incendie (dates, lieu, surface), et retourner le tout dans une table

scrapping = function(feux, page) {
  list_a = html_nodes(read_html(feux), "a")
  if (length(list_a) &gt;= 2) {
    att = html_attrs(list_a[2])
    longitude = as.numeric(str_extract(att, "[0-9]\\.([0-9])*"))
    latitude = as.numeric(str_extract(att, "([0-9]){2}\\..([0-9])*"))
  } else {
    longitude = 0
    latitude = 0
  }
  list_span = html_nodes(read_html(feux), "span")
  date = html_text(list_span[2])
  commune = html_text(list_span[4])
  surface = as.numeric(html_text(list_span[5]))  
  list_div = html_nodes(read_html(feux), "div")
  if (length(list_a) &gt;= 2) {
    infos1 = html_text(list_div[3])
    infos2 = html_text(list_div[5])
  } else {
    infos1 = html_text(list_div[2])
    infos2 = html_text(list_div[4])
  }
  alerte = gsub("First alert : |First intervention", "", str_extract(infos1, "First alert : (.)*First intervention"))
  deb_inter = gsub("First intervention : |End of intervention", "", str_extract(infos1, "First intervention : (.)*End of intervention"))
  fin_inter = gsub("End of intervention : |Origins of the alert", "", str_extract(infos1, "End of intervention : (.)*Origins of the alert"))
  code_INSEE = gsub("INSEE code : |Locality", "", str_extract(infos2, "INSEE code : (.)*Locality"))
  lieudit = str_trim(gsub("Locality : ", "", str_extract(infos2, "Locality : (.)*")))
  liste_feux = data.frame(commune, date, code_INSEE, longitude, latitude, lieudit, surface, alerte, deb_inter, fin_inter, page)
  return(liste_feux)
}

Reste a faire tourner sur les 620 pages. Pour chaque ligne, on récupère les informations !

df = data.frame()
for (page in 1:page_max) {
  text_base = read_html(remDr$getPageSource()[[1]])
  list_tr = html_nodes(text_base, "tr")
  if (length(html_nodes(read_html(as.character(list_tr[14])), "a")) == 8) {
    deb = 15
  } else {
    deb = 14
  }
  for (i in deb:(deb+49)) {
    feux = as.character(list_tr[i])
    df = rbind(df, scrapping(feux, page))
  }
  if (page == 1) {
    webElem = remDr$findElement(using = 'xpath', '//*[contains(concat( " ", @class, " " ), concat( " ", "page_img", " " ))]')
  } else {
    webElem = remDr$findElement(using = 'xpath', '//*[contains(concat( " ", @class, " " ), concat( " ", "page_img", " " ))][3]')
  }
  webElem$clickElement()
  heure = (floor((page_max-page)*15/60) - floor((page_max-page)*15/60) %% 60) / 60
  minute = floor((page_max-page)*15/60) %% 60
  pb = progress_bar$new(total = 100)
  for (i in 1:100) {
    if (i==1) {cat("page", page, "sur", page_max," -  (reste", heure, "h", minute, "min)",  "\n")}
    pb$tick()
    Sys.sleep(15 / 100)
  }
  remDr$screenshot(TRUE)
}
feux = df

Ça y est, on a une base, avec presque 31 000 lignes ! on va ensuite ranger un peu, et clarifier

colnames(feux)[colnames(feux)=="longitude"] = "x"
colnames(feux)[colnames(feux)=="latitude"] = "y"
feux[, date := as.Date(date, "%d/%m/%Y")]
feux[, date_alerte := as.Date(substr(alerte, 1, 10), "%d/%m/%Y")]
feux$heure_alerte[nchar(feux$alerte) &gt; 10] = substr(feux$alerte[nchar(feux$alerte) &gt; 10], 12, 16)
feux$heure_alerte[!(nchar(feux$alerte) &gt; 10)] = ""
feux[, date_deb_inter := substr(deb_inter, 1, 10)]
feux[, date_deb_inter := as.Date(date_deb_inter, "%d/%m/%Y")]
feux$heure_deb_inter[nchar(feux$deb_inter) &gt; 10] = substr(feux$deb_inter[nchar(feux$deb_inter) &gt; 10], 12, 16)
feux$heure_deb_inter[!(nchar(feux$deb_inter) &gt; 10)] = ""
feux$heure_deb_inter[feux$heure_deb_inter == "00:00"] = ""
feux[, date_fin_inter := substr(fin_inter, 1, 10)]
feux[, date_fin_inter := as.Date(date_fin_inter, "%d/%m/%Y")]
feux$heure_fin_inter[nchar(feux$fin_inter) &gt; 10] = substr(feux$fin_inter[nchar(feux$fin_inter) &gt; 10], 12, 16)
feux$heure_fin_inter[!(nchar(feux$fin_inter) &gt; 10)] = ""
feux$heure_fin_inter[feux$heure_fin_inter == "00:00"] = ""
feux$annee = year(feux$date_alerte)
feux$month = month(feux$date_alerte)

Pour voir un peu nos données, rien de tel que des statistiques descriptives. On peut commencer avec du barplot du nombre de feux de forêts par année

agg_annee = aggregate(alerte~annee, FUN=length, data=feux)
colnames(agg_annee)[colnames(agg_annee)=="alerte"] = "Nombre"
nb = ggplot(agg_annee) +
  geom_bar(aes(x=annee, y=Nombre, fill="Feux de forêts"), stat="identity", width = 0.8) +
  scale_colour_manual(name=NULL, values=c("orange")) +
  ggtitle("Nombre de feux par an") +
  theme(plot.title = element_text(hjust = 0.5, size=14),
        legend.position = "none",
        axis.text.x = element_text(angle=50, vjust=0.5, size=10),
        axis.text.y = element_text(size=12),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  scale_fill_discrete("") +
  scale_x_discrete(limit=seq(1981, 2017, by = 2))

ou un barplot de la surface de feux de forêts par anee

agg_surface = aggregate(surface~annee, FUN=sum, data=feux)
surf = ggplot(agg_surface) +
  geom_bar(aes(x=annee, y=surface, fill="Feux de forêts"), stat="identity", width = 0.8) +
  scale_colour_manual(name=NULL, values=c("orange")) +
  ggtitle("Surface totale incendiée par an (en ha)") +
  theme(plot.title = element_text(hjust = 0.5, size=14),
        legend.position = "none",
        axis.text.x = element_text(angle=50, vjust=0.5, size=10),
        axis.text.y = element_text(size=12),
        axis.title.x = element_blank(),
        axis.title.y = element_blank()) +
  scale_fill_discrete("") +
  scale_x_discrete(limit=seq(1981, 2017, by = 2)) +
  scale_y_discrete(limit=seq(0, 25000, by = 5000))

Et on fait le graphique

require(grid)
grid.newpage()
pushViewport(viewport(layout = grid.layout(1, 2)))
define_region = function(row, col){
  viewport(layout.pos.row = row, layout.pos.col = col)
} 
print(nb, vp = define_region(1, 1))
print(surf, vp = define_region(1, 2))

Passons aux choses serieuses, avec des cartes (la partie la plus serieuse, on en conviendra, c’etait la phase initiale de scraping). On commence par des informations relatives a la Corse

url="https://raw.githubusercontent.com/data-PM/Promethee/master/communes_corse-shp.zip"
temp="corse.zip"
download.file(url, "corse.zip")
unzip("corse.zip")
communes_corse = sp::spTransform(shapefile("/Users/arthur/communes_corse-shp/communes_corse.shp"), CRS("+proj=longlat +datum=WGS84"))

On se limite aux feux de plus de 1000m2

feux_min1000m2 = feux[feux$surface &gt; 0.1, ]

Maintenant, on peut faire la carte

feux_INSEE = aggregate(x = feux_min1000m2$code_INSEE, by = list(unique.values = feux_min1000m2$code_INSEE), FUN = length)
colnames(feux_INSEE)[colnames(feux_INSEE)=="unique.values"] = "insee"
communes_corse = merge(x = communes_corse, y = feux_INSEE, by = "insee", all.x = TRUE)
communes_corse$x[is.na(communes_corse$x)] = 0
max = max(communes_corse$x)
communes_corse$x2 = communes_corse$x/max * 100
corse2 = gSimplify(communes_corse, tol = 0.001, topologyPreserve = TRUE)
corse2 = SpatialPolygonsDataFrame(corse2, communes_corse@data, match.ID = FALSE)
opar = par(mar = c(0,0,1.2,0))
plot(corse2, col = "#DAE3E6", border = "#8A0641", lwd = 0.7, bg = "#B5D0D0")

Le fond de carte (avec les communes) ressemble a ca

On peut maintenant passer a la carte choroplethe

choroLayer(spdf = corse2, df = corse2@data, var = "x", border = NA, 
           col = carto.pal("wine.pal", 6), legend.pos = "topleft", 
           add = TRUE, method = "fisher-jenks", nclass = 6, legend.title.txt = NULL)
top_villes = corse2[corse2$nom %in% c("Ajaccio", "Bastia", "Corte", "Porto-Vecchio", "Propriano", "Bonifacio", "Calvi"), ]
graphics::points(coordinates(top_villes), pch = 20, cex = 0.5)
labelLayer(top_villes, top_villes@data, txt = "nom", cex = 1, pos = 2, font = 4, offset = 0.5)
layoutLayer(title = "Nombre de feux de forêts depuis 1981 (minimum 1000 m2)",
            source = "Source : 1973-2018 Prométhée ©",
            author = NULL, scale = 10, frame = FALSE, col = "#cdd2d4", coltitle = "#8A5543")

Comme souvent avec R, plusieurs formats sont envisageables. On peut passer a du ggplot2, mais on commence par aller chercher un shapefile global de la Corse

corse = shp_zip("https://raw.githubusercontent.com/data-PM/Promethee/master/corse-shp.zip", "corse")
corse = spTransform(corse, CRS("+proj=longlat +datum=WGS84"))
corse = fortify(corse)
p=data %&gt;%
  ggplot() +
  geom_polygon(data=corse, aes(x=long, y = lat, group = group), fill="grey", alpha=0.3) +
  geom_point(data=feux_min100ha, aes(x=x, y=y, size=surface, color=surface, text=mytext, alpha=surface) ) +
  labs(title = "\nPrincipaux feux de forêts depuis 1981 \n (plus de 100 Ha)") +
  xlim(8.5, 9.6) + ylim(41.3,43.1) +
  scale_size_continuous(range=c(1,5), ) +
  scale_color_viridis(option="inferno", trans="log", name="Surface", limits=c(99, 5700), breaks=c(100, 200, 500, 1000, 2500, 5000)) +
  scale_alpha_continuous(trans="log") +
  theme_void() +
  theme(legend.title = element_text(size=10, face="bold"))
t=800
div(ggplotly(tooltip="text", width=t/1.7, height=t), align = "center")

Joli, non ?

Carte météo dynamique

Pour aller un peu plus loin par rapport au précédent billet, on peut faire une carte dynamique, pour visualiser une tempête, sur une journée. Je vais partir du fait qu’on a encore les objets du précédant billet en mémoire.

On commence par garder uniquement les données relatives a la journée qui nous intéresse,

date_sel = "2016-10-13"
date_sel_OK = paste0(substr(date_sel,9,10),"/",substr(date_sel,6,7),"/",substr(date_sel,1,4))
donnees_carte_0 = subset(data, data$date==as.POSIXct(date_sel))
horaires = sort(unique(donnees_carte_0$heure))

(mais on pourrait prendre une semaine). On recupere ensuite le contour de la France,

download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds","FRA_adm0.rds")
FR0=readRDS("FRA_adm0.rds")
P1=FR0@polygons[[1]]@Polygons[[355]]@coords
P2=FR0@polygons[[1]]@Polygons[[27]]@coords

car seule la France métropolitaine nous intéressera (incluant la Corse, soit 2 gros polygones) dans l’esprit du précédant billet sur le zonier. On va ainsi definir un maillage sur lesquel on va lisser la pluviometrie,

grille = expand.grid(seq(min(donnees_carte_0$longitude),max(donnees_carte_0$longitude),length=101),seq(min(donnees_carte_0$latitude),max(donnees_carte_0$latitude),length=101))
paslong=(max(donnees_carte_0$longitude)-min(donnees_carte_0$longitude))/100
paslat=(max(donnees_carte_0$latitude)-min(donnees_carte_0$latitude))/100
f=function(i){ (point.in.polygon (grille[i, 1]+paslong/2 , grille[i, 2]+paslat/2 , P1[,1],P1[,2])&gt;0)+(point.in.polygon (grille[i, 1]+paslong/2 , grille[i, 2]+paslat/2 , P2[,1],P2[,2])&gt;0) }
indic=unlist(lapply(1:nrow(grille),f))
grille=grille[which(indic==1),]

Pour lisser, on utilise les k-plus proches voisins

knn=function(i,k=10){
  d=distHaversine(grille[i,1:2],donnees_carte[,
    c("longitude","latitude")], r=6378.137) 
  r=rank(d)
  ind=which(r&lt;=k)
  weighted.mean(donnees_carte[ind,"rr3"],(1/d[ind])/sum(1/d[ind]))}

On a ensuite la fonction suivante, pour faire une carte

carto_prec&lt;-function(){
  grille2&lt;-grille
  grille2$rr=Vectorize(knn)(i=1:nrow(grille2))
  bk=seq(0,50,length=21)
  grille2$cuty=cut(grille2$rr,breaks=bk,labels=1:20)
  cols = rev(carto.pal(pal1 = "blue.pal", n1=20, pal2 = "white.pal", n2=1))
  plot(FR0,border=NA)
  polygon(P1)
  polygon(P2)
  points(grille2[,1]+paslong/2,grille2[,2]+paslat/2,col=cols[grille2$cuty],pch=19)
  points(donnees_carte$longitude,donnees_carte$latitude, col="black",pch=19,cex=.5)
  title(main = paste0(date_sel_OK," à ",heure_sel,"H"),
        sub = "Précipitations des 3 dernières heures",
        cex.main = 1.5,   font.main= 4,
        cex.sub = 1, font.sub = 3)
  legend(8.2, 50, legend=seq(0,50,length=5), title='en mm',
         fill=cols[seq(1,20,length=5)], cex=0.8)
}

On va ensuite utiliser 8 fichiers – car on regarde toutes les 3 heures

for (hh in horaires) {
  heure_sel&lt;-hh
  donnees_carte &lt;- subset(donnees_carte_0, donnees_carte_0$heure==heure_sel)
  png(paste0("CartePrec",date_sel,"_",hh,".png"))
  carto_prec()
  dev.off()
}

(mais ça peut se raffiner, avec un lissage spatio-temporel). Et on va ensuite créer notre animation en concaténant les images

frames &lt;- image_morph(
  c(image_scale(image_read("CartePrec2016-10-13_1.png")),
image_scale(image_read("CartePrec2016-10-13_3.png")),
image_scale(image_read("CartePrec2016-10-13_6.png")),
image_scale(image_read("CartePrec2016-10-13_9.png")),
image_scale(image_read("CartePrec2016-10-13_12.png")),
image_scale(image_read("CartePrec2016-10-13_15.png")),
image_scale(image_read("CartePrec2016-10-13_18.png")),
image_scale(image_read("CartePrec2016-10-13_21.png"))), frames = 16)
frames.anim &lt;- image_animate(frames)
image_write(frames.anim, paste0("CartePrecDay",date_sel,".gif"))

Il faut avouer que ça a de gueule… On voit clairement la tempête arriver…

Faire ses cartes météo

Allez, un peu de scraping aujourd’hui (promis, il y aura d’autres billets d’ici la fin de la semaine sur le sujet). Je vais m’inspirer des codes de Romain de la formation Data Science pour l’Actuariat. Le but est de récupérer des données de températures, de précipitations, de vitesse de vent, journalières, en France. On commence par charger quelques librairies,

library(plyr)
library(stringr)
library(OpenStreetMap)
library(leaflet)
library(shiny)
library(rsconnect)
library(mapview)
library(png)
library(magick)
library(yaml)

On va récupérer les deux dernières années, mais comme toujours, le code peut facilement s’adapter.

annee = c(2016:2017)
mois = c('01','02','03','04','05','06','07','08','09','10','11','12')
aaaamm = sort(sub(pattern=' ',replacement = '', x=outer(annee,mois,paste)))
myCols = c(NA,NA,"NULL","NULL","NULL",NA,NA,"factor","NULL","NULL","NULL",NA,"NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL",NA,NA,NA,NA,NA,"NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL","NULL")
for (mois in aaaamm) {
  fichier = paste0("synop_",mois,".csv.gz")
  if (!file.exists(fichier)) {
    page = paste0('https://donneespubliques.meteofrance.fr/donnees_libres/Txt/Synop/Archive/synop.',mois,'.csv.gz')
    download.file(page, fichier, quiet = TRUE, cacheOK = FALSE)
  }
}  
for (mois in aaaamm) {
  fichier = paste0("synop_",mois,".csv.gz")
  if (mois=='201601') {
    data = read.table(gzfile(fichier),header=TRUE, sep=";",quote='"',fileEncoding="UTF-8",colClasses = myCols)
  } else {
    data = rbind(data,read.table(gzfile(fichier),header=TRUE, sep=";",quote='"',fileEncoding="UTF-8",colClasses = myCols))
  }
}

On supprime une des colonnes qui ne sert a rien

data[,"X"] = NULL

et on bricole un peu (sinon on a des soucis avec minuit, qui est décalé de 24 heures)

data$temp = as.numeric(substr(data$date,1,14))
decalage_minuit = function(dateheure){
  r = dateheure
  if (dateheure-trunc(dateheure/1000000)*1000000==0){
    r = dateheure + 10000
  }
  return(r)
}
data$temp = sapply(data$temp,decalage_minuit)
data$jourheure_obs = strptime(data$temp, format = '%Y%m%d%H%M%S', 'UTC')
data$date = as.POSIXct((strptime(data$jourheure_obs,format = "%Y-%m-%d")))
data$heure = as.numeric(substr(data$jourheure_obs, start=12, stop=13))
data[,"temp"] = NULL
data$t = as.numeric(as.character(data$t))-273.15

oui, on doit convertir en degrés Celsius. On peut aussi travailler sur les vitesses de vent

data$dd = as.numeric(as.character(data$dd))
data$ff = as.numeric(as.character(data$ff))/1000*3600
data$ff = round(data$ff/5)*5

que l’on va convertir en km/h, et enfin les précipitations

data$rr1 = pmax(0,as.numeric(as.character(data$rr1)))
data$rr3 = pmax(0,as.numeric(as.character(data$rr3)))
data$rr6 = pmax(0,as.numeric(as.character(data$rr6)))
data$rr12 = pmax(0,as.numeric(as.character(data$rr12)))
data$rr24 = pmax(0,as.numeric(as.character(data$rr24)))

On peut d’ailleurs commencer par une simple visualisation de séries de températures

data_temperature = subset(data, !is.na(data$t))
temp_min = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=min)
names(temp_min)[2] = "temp_min"
temp_max = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=max)
names(temp_max)[2] = "temp_max"
temp_moy = aggregate(x=data_temperature$t,by=list(date=data_temperature$date),FUN=mean)
names(temp_moy)[2] = "temp_moy"
temp_stat = cbind(temp_min,temp_max[],temp_moy)
rm(temp_min,temp_max,temp_moy)
temp_stat[,3] = NULL
temp_stat[,4] = NULL

avec le maximum observé, par jour, sur toutes les stations, le minimum, et la moyenne des stations

plot(x=temp_stat$date,y=temp_stat$temp_moy,'l',ylim=c(-15,40),xlab='Date',ylab='Température (°C)')
title('Températures françaises métropolitaines de 01/2016 à 12/2017')
points(x=temp_stat$date,y=temp_stat$temp_min,col="blue",'l')
points(x=temp_stat$date,y=temp_stat$temp_max,col="red",'l')

Pour faire des jolies cartes, il nous faut des informations sur les stations météo.

stations = read.table('https://donneespubliques.meteofrance.fr/donnees_libres/Txt/Synop/postesSynop.csv',header=TRUE, sep=";",quote='"',fileEncoding="UTF-8")
stations = rename(stations, replace=c("ID"="numer_sta"))
stations = rename(stations, replace=c("Nom"="nom"))
stations = rename(stations, replace=c("Latitude"="latitude"))
stations = rename(stations, replace=c("Longitude"="longitude"))
stations = rename(stations, replace=c("Altitude"="altitude"))
data = merge(x=data,y=stations)

Pour la carte, on va se restreindre aux données en France métropolitaine

data = subset(data, data$latitude &gt;= 39)
data = subset(data, data$latitude &lt;= 53)
data = subset(data, data$longitude &gt;= -5.3)
data = subset(data, data$longitude &lt;= 9.8)
coordonnees_stations = unique(data[, c("nom","latitude","longitude")])

Commencons par visualiser le maximum observee le 19 juillet 2016

date_sel = "2016-07-19"
date_sel_OK = paste0(substr(date_sel,9,10),"/",substr(date_sel,6,7),"/",substr(date_sel,1,4))
donnees_carte_0 = subset(data, data$date==as.POSIXct(date_sel))
temp = aggregate(x=donnees_carte_0$t,by=list(numer_sta=donnees_carte_0$numer_sta),FUN="max")
donnees_carte_0 = merge(x=donnees_carte_0,y=temp)
donnees_carte_0 = rename(donnees_carte_0, replace=c("x"="TEMP"))
donnees_carte_0$TEMP = round(donnees_carte_0$TEMP)
donnees_carte = unique(donnees_carte_0[,c('numer_sta','nom','longitude','latitude','TEMP')])
get_france = get_map(c(lon=2.25,lat=46), zoom=5, col='bw')
colfunc = colorRampPalette(c("darkblue","royalblue","cyan","lightblue","orange","red","darkred"))
liste_couleurs = data.frame(TEMP=as.numeric(seq(-49,50,1)),couleur=as.character(colfunc(100)))
donnees_carte = merge(x=donnees_carte,y=liste_couleurs)
france_tmax = ggmap(get_france) + 
  ggtitle(paste0('Températures maximales du ', date_sel_OK)) +
  scale_x_continuous(limits = c(-5, 10), expand = c(0, 0)) +
  scale_y_continuous(limits = c(41.3, 51.1), expand = c(0, 0))
france_tmax = france_tmax + geom_label(data=donnees_carte,aes(x=longitude,y=latitude,label=TEMP),col=as.character(donnees_carte$couleur))
france_tmax

On peut tenter aussi une visualisation des précipitations. On se place le 13 octobre 2016 (oui, il y a eu des inondations dans la region de Montpelier ce jour la)

date_sel = "2016-10-13"
colfunc = colorRampPalette(c('white','darkblue'))
data$precip_en_mm = data$rr3
data_precipitation = subset(data, !is.na(data$precip_en_mm))
precipitations = aggregate(x=data_precipitation$precip_en_mm,by=list(date=data_precipitation$date,station=data_precipitation$nom),FUN=sum)
names(precipitations)[3] = "pp"
precipitations$pp = pmax(precipitations$pp,0)
donnees_carte = subset(precipitations, precipitations$date==date_sel)
qt_prec = quantile(round(donnees_carte$pp),seq(0, 1, 0.05))
donnees_carte$Couleur_prec = colfunc(21)[(findInterval(round(donnees_carte$pp), qt_prec, all.inside=TRUE))]
donnees_carte = rename(donnees_carte, replace=c("station"="nom"))
donnees_carte = merge(x=donnees_carte,y=coordonnees_stations)

On peut tenter un leaflet ici (ça sera plus joli)

france = leaflet() %&gt;% addTiles() %&gt;% fitBounds(lng1=-4.412, lat1=41.92, lng2=9.485, lat2=50.57)
france = france %&gt;% addCircles(lng=donnees_carte$longitude, lat=donnees_carte$latitude, color=donnees_carte$Couleur_prec, opacity = 1, fillColor=donnees_carte$Couleur_prec, fillOpacity = 1, radius=donnees_carte$pp*500)
france

Comme souvent, j’ai des soucis pour integrer du leaflet en wordpress, alors on se contentera dans le billet d’une copie d’ecran (mais un copier/coller du code suffit).

Et pour finir, la meme chose avec les vitesses de vent. On se place le 6 mars 2017 (la encore, la date n’est pas choisie au hasard, c’était la tempete Zeus)

date_sel = "2017-03-06"
date_sel_OK = paste0(substr(date_sel,9,10),"/",substr(date_sel,6,7),"/",substr(date_sel,1,4))
heure_sel = 18
donnees_carte_0 = subset(data, data$date==as.POSIXct(date_sel))
donnees_carte_0 = subset(donnees_carte_0, donnees_carte_0$heure==heure_sel)
donnees_carte = donnees_carte_0

Pour les precipitations, des cercles plus ou moins grands suffisaient. La, on va mettre des fleches

logo = image_trim(image_read("https://image.freepik.com/icones-gratuites/fleche-noire-vers-le-haut_318-30934.jpg"))
if (file.exists('flecheblanche90.png')==FALSE){
  for (angle in seq(10,360, by=10)) {
    temp = image_rotate(logo,angle)
    image_write(temp,paste0('flecheblanche',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='green')
    image_write(temp,paste0('flecheverte',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='yellow')
    image_write(temp,paste0('flechejaune',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='orange')
    image_write(temp,paste0('flecheorange',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='red')
    image_write(temp,paste0('flecherouge',angle,'.png'))
    temp = image_colorize(image=temp,opacity=40,color='black')
    image_write(temp,paste0('flechenoire',angle,'.png'))
  }
}
arrowIcons = icons(
  iconUrl = ifelse(donnees_carte$dd == 0,
            ifelse(donnees_carte$ff&lt;39,
                    'flecheverte360.png',
            ifelse(donnees_carte$ff&lt;79,
                    'flechejaune360.png',
            ifelse(donnees_carte$ff&lt;119,
                    'flecheorange360.png',
            ifelse(donnees_carte$ff&lt;159,
                    'flecherouge360.png',
            ifelse(donnees_carte$ff&gt;159,
                    'flechenoire360.png',
                    'flecheblanche360.png'))))),
            ifelse(donnees_carte$ff&lt;39,
                   paste0('flecheverte',donnees_carte$dd,'.png'),
            ifelse(donnees_carte$ff&lt;79,
                    paste0('flechejaune',donnees_carte$dd,'.png'),
            ifelse(donnees_carte$ff&lt;119,
                    paste0('flecheorange',donnees_carte$dd,'.png'),
            ifelse(donnees_carte$ff&lt;159,
                    paste0('flecherouge',donnees_carte$dd,'.png'),
            ifelse(donnees_carte$ff&gt;159,
                    paste0('flechenoire',donnees_carte$dd,'.png'),
           paste0('flecheblanche',donnees_carte$dd,'.png'))))))),
  iconWidth = 30, iconHeight = 30)

Cette fois, on est bon ! On peut faire la carte

france = leaflet(options = leafletOptions()) %&gt;% addTiles() %&gt;% fitBounds(lng1=-4.412, lat1=41.92, lng2=9.485, lat2=50.57)
france = addMarkers(map = france,lng=donnees_carte$longitude,lat=donnees_carte$latitude,icon=arrowIcons,popup=as.character(donnees_carte$ff))
france

An Open Lab-Notebook Experiment