Tag Archives: computer

Les votes à l’assemblée nationale

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

require(xml2)
require(downloader)
require(stringr)
require(classInt)
require(plotrix)
require(FactoMineR)
require(sp)

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

path="http://data.assemblee-nationale.fr/static/openData/repository/AMO/deputes_senateurs_ministres_legislature/AMO20_dep_sen_min_tous_mandats_et_organes_XIV.csv.zip"
dest="deputes.zip"
download(path,destfile=dest,mode="wb")
loc=paste(getwd(),"/",dest,sep="" )
unzip(loc)
dest="acteurs.csv"
loc=paste(getwd(),"/",dest,sep="" )
t2=read.csv(loc,sep=";")
res=NULL
for (i in 1:248){
test=str_c("mandats.1..mandat.",i,"..typeOrgane.1.")
j=which(colnames(t2) == test)
#Je ne récupère que les personnes dont le i ème mandat est un mandat de député i.e siegeant à l'Assemblée Nationale
t3=subset(t2,t2[,j]=="ASSEMBLEE")
if(nrow(t3)!=0){
#Nom des différentes variables en concatenant les chaines
circo=str_c("mandats.1..mandat.",i,"..election.1..lieu.1..numCirco.1.")
dept=str_c("mandats.1..mandat.",i,"..election.1..lieu.1..numDepartement.1.")
acteur=str_c("mandats.1..mandat.",i,"..acteurRef.1.")
mandat=str_c("mandats.1..mandat.",i,"..uid.1.")
# Récupération du numero de la colonne du data frame dont le nom correspond a chacune des variable voulues
k=which(colnames(t3) == circo)
l=which(colnames(t3) == dept)
m=which(colnames(t3) == acteur)
n=which(colnames(t3) == mandat)
#Je cree un data frame correspondant "au numero de mandat"
t4=data.frame(as.numeric(as.character(t3[,k])),as.numeric(as.character(t3[,l])),str_c(t3[,m],t3[,n]))
colnames(t4)=c("circo","dept","identifiant")
#Et j'ajoute l'ensemble des députés
res=rbind(res,t4)
}
}
res$iden=str_c(res$dept,res$circo)
res$nbVote=0
res$oui=0
res$non=0
res$abst=0
res$nomOrgane=0

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

path="http://data.assemblee-nationale.fr/static/openData/repository/LOI/scrutins/Scrutins_XIV.xml.zip"
dest="Scrutins_XIV.xml.zip"
download(path,destfile=dest,mode="wb")
loc=paste(getwd(),"/",dest,sep="" )
unzip(loc)
dest="Scrutins_XIV.xml"
loc=paste(getwd(),"/",dest,sep="" )

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

t=read_xml(loc)
liste=xml_children(t)

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

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

for (a in liste){
if(xml_text(xml_children(a)[14])=="DecompteNominatif"){
numero=xml_text(xml_children(a)[2])
ventil=xml_children(a)[16]
groupe=xml_children(xml_children(xml_children(ventil)))
tempDataFrame=data.frame(res$identifiant,NA)
colnames(tempDataFrame)=c("identifiant",str_c("scrutin",numero))
for (b in groupe){
intermediaire=xml_children(xml_children(b))
nomGroupe=xml_text(xml_children(b)[1])
for (i in 1:4){
df3=data.frame(xml_text(xml_children(xml_children(intermediaire[3])[i])))
if(nrow(df3)!=0){
for (j in df3[,]){
if (i==1){
j=strsplit(strsplit(strsplit(j,"MG")[[1]],"PSE")[[1]],"PAN")[[1]]
}
res[res$identifiant==j,]$nomOrgane=as.character(nomGroupe)
if(i!=1)
res[res$identifiant==j,]$nbVote=res[res$identifiant==j,]$nbVote+1
if (i==1)
tempDataFrame[tempDataFrame$identifiant==j,2]="NV"
if (i==2){
res[res$identifiant==j,]$oui=res[res$identifiant==j,]$oui+1
tempDataFrame[tempDataFrame$identifiant==j,2]="POUR"
}
if (i==3){
res[res$identifiant==j,]$non=res[res$identifiant==j,]$non+1
tempDataFrame[tempDataFrame$identifiant==j,2]="CTRE"
}
if (i==4){
res[res$identifiant==j,]$abst=res[res$identifiant==j,]$abst+1
tempDataFrame[tempDataFrame$identifiant==j,2]="ABST"
}}}}}
res=data.frame(res, tempDataFrame[,2])
}}
res2=subset(res,res$dept<96)
res2$circo=ifelse(nchar(res2$circo)==1, str_c("0",res2$circo) , str_c("",res2$circo))
res2$dept=ifelse(nchar(res2$dept)==1, str_c("0",res2$dept) , str_c("",res2$dept))
res2$iden=str_c(res2$dept,res2$circo)
res2=subset(res2,!is.na(iden))
res2=res2[order(res2$iden),]
parIden=aggregate(res2$nbVote, by=list(res2$iden), sum)

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

path="http://www.laspic.eu/data/documents/circosshp_v3.zip"
dest="circosshp_v3.zip"
download(path,destfile=dest,mode="wb")
loc=paste(getwd(),"/",dest,sep="" )
unzip(loc)
dest="circosSHP_v3.RData"
loc=paste(getwd(),"/",dest,sep="" )
load(loc)

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

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

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

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

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

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

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

Accidents de voiture et de deux roues à Paris

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

Comme toujours on charge quelques packages

require(downloader)
require(ggmap)
require(stringr)
options(encoding = "UTF-8")

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

path="https://opendata.paris.fr/explore/dataset/accidentologie/download/?format=csv&timezone=Europe/Berlin&use_labels_for_header=true"
dest="test"
download(path,destfile=dest,mode="wb")
loc=paste(getwd(),"/",dest,sep="" )
obj <- read.csv(file = loc, stringsAsFactor = FALSE,sep=";")
file.remove(dest)

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

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

et créer une base qui contient toutes les informations

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

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

names(vec_dat)[1]="Jour";names(vec_dat)[2]="Nbsinistre"

On choisit éventuellement une palette de couleurs

Palette<-colorRampPalette(c("lightblue","darkblue"))

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

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

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

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

Passons maintenant à de la géolocalisation

t2 = subset(t, Coordonnées!= "")
coord=data.frame(t(matrix(as.numeric(unlist(strsplit(t2$Coordonnées,"[,]"))),nrow=2)))
t2$lat=coord[,1]
t2$lon=coord[,2]
t2 = subset(t2, (
(48.81 <= lat) & (lat<=48.90) & (2.25<=lon) &(lon<=2.43) ) )

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

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

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

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

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

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

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

Visualiser les résultats du premier tour

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

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

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

library(maptools)
library(maps)
departements<-readShapeSpatial("DEPARTEMENT.SHP")
plot(departements)
points(departements@data$X_CENTROID,departements@data$Y_CENTROID,pch=19,col="red")

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

pos=which(departements@data[,"CODE_DEPT"]==35)
Poly_35=departements[pos,]
plot(departements)
plot(Poly_35,col="yellow",add=TRUE)
departements@data[pos,c("X_CENTROID","Y_CENTROID")]
points(departements@data[pos,c("X_CENTROID","Y_CENTROID")],pch=19,col="red")
library(rgeos)
(ctd=gCentroid(Poly_35,byid=TRUE))
points(ctd,pch=19,col="blue")

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

ctd=as.data.frame(gCentroid(departements,byid=TRUE))
plot(departements)
points(ctd,pch=19,col="blue")

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

candidat="M. Emmanuel MACRON"
library(XML)
voix=function(no=35){
testurl=FALSE
i=1
vect_reg=c("084","027","053","024","094","044","032","028","075","076","052",
"093","011","001","002","003","004")
region=NA
while((testurl==FALSE)&(i<=20)){
reg=vect_reg[i]
nodpt=paste("0",no,sep="")
# if(!is.na(as.numeric(no))){if(as.numeric(no)<10) nodpt=paste("00",no,sep="")}
url=paste("http://elections.interieur.gouv.fr/presidentielle-2017/",reg,"/",nodpt,"/index.html",sep="")
test=try(htmlParse(url),silent =TRUE)
if(!inherits(test, "try-error")){testurl=TRUE
region=reg}
i=i+1
}
tabs tab=tabs[[2]]
nb=tab[tab[,1]==candidat,"Voix"]
a<-unlist(strsplit(as.character(nb)," "))
as.numeric(paste(a,collapse=""))}

On peut alors tester

> voix(35)
[1] 84648

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

liste_dpt=departements@data$CODE_DEPT
nbvoix=Vectorize(voix)(liste_dpt)

On peut alors visualiser sur une carte.

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

Et on peut aussi tenter pour une autre candidate,

candidat="Mme Marine LE PEN"

et on obtient la carte suivante

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

Visualizing (censored) lifetime distributions

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

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

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

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

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

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

This is more interesting when we have different lifetimes

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

or, with a different time granularity

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

Nice, isn’t it?

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

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

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

url="https://en.wikipedia.org/wiki/List_of_wars_involving_the_United_States"
download.file(url,destfile = "warUS.html")
url="https://en.wikipedia.org/wiki/List_of_wars_involving_France"
download.file(url,destfile = "warFR.html")
url="https://fr.wikipedia.org/wiki/Liste_des_guerres_de_la_France#Premi.C3.A8re_R.C3.A9publique"
download.file(url,destfile = "guerre.html")
url="https://en.wikipedia.org/wiki/List_of_wars_involving_Canada"
download.file(url,destfile = "warCAN.html")

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

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

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

https://freakonometrics.hypotheses.org/files/2017/03/guerre-us1.png

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

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

library(stringr)
ext_date=function(x){
dates12="[0-9]{4}"
#grep(pattern = dates2, x = col1[1])
L=str_extract_all(as.character(x),dates12)
return_L=list()
if(length(L)>0){
for(j in 1:length(L))
if(length(L[[j]])==1) return_L[[j]]=as.numeric(L[[j]])
if(length(L[[j]])>=2) return_L[[j]]=seq(min(as.numeric(L[[j]])),max((as.numeric(L[[j]]))))
}
return(return_L)}

For the US, we get the following years

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

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

> length(d)
[1] 222

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

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

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

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

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

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

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

Let’s try another one. Like Canada,

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

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

Reading text automatically

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Juxtaposer des cartes et ‘confounding’

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Evolution du nombre de demandeurs d’emploi

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

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

Pour charger les packages, Meryem propose une jolie petite fonction

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

On va retravailler un peu les données pour commencer

x<-ufc$nombre
y<-ufc$demande
dt<-as.Date(ufc$date,"%d/%m/%Y")
n<-length(x)

et ensuite on fait notre animation,

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

et voilà

Cartographier le vote pour le Brexit

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

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

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

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

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

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

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

Bon, arrêtons là….

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

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

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

> sum(referendum$Leave)/(sum(referendum$Leave)+sum(referendum$Remain))
[1] 0.5189184

On peut regarder, région par région, si c’est le leave ou le remain qui a gagné, en utilisant

referendum=referendum[c(referendum$Region!="Northern Ireland"),]
referendum=referendum[c(referendum$HASC_code!="Gibraltar"),]
row.names(referendum)=seq(1,nrow(referendum),1)
leave_or_remain=cbind(referendum,"Brexit?"=0)
leave_or_remain[,"Brexit?"]=ifelse(leave_or_remain$Remain<leave_or_remain$Leave,rgb(1,0,0,.7),rgb(0,0,1,.7))
map_data=data.frame(UK@data)
map_data=cbind(map_data,"Brexit"=0)
for (i in 1:nrow(map_data)){
if(map_data[i,"NAME_1"]=="Northern Ireland"){
map_data[i,"Brexit"]="blue"}else{
map_data[i,"Brexit"]=as.character(leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],'Brexit?'])
}
}
plot(UK, col = map_data$Brexit, border = "gray1", xlim = c(-4,-2), ylim = c(50, 59), main="How the UK has voted?", bg="#A6CAE0")
plot(IRL, col = "lightgrey", border = "gray1",add=TRUE)
plot(FR, col = "lightgrey", border = "gray1",add=TRUE)

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

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

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

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

Tests non-paramétriques et simulations

Lors du dernier cours de statistique, nous avons présenter les tests d’ajustment de lois. Nous avions illustré ces tests à partir de la taille d’individus (déjà utilisé pour présenter l’ajustement de lois, et l’estimation de densité), correspondant à https://latex.codecogs.com/gif.latex?\boldsymbol{x}=\{x_1,\cdots,x_n\}.

> 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
> n=length(X)

On notera https://latex.codecogs.com/gif.latex?%20(x_{i:n}) la statistique d’ordre, au sens où

https://latex.codecogs.com/gif.latex?%20x_{1:n}\leq%20x_{2:n}\leq\cdots\leq%20x_{n-1:n}\leq%20x_{n:n}

Parmi les outils graphiques, nous avons vu le PP plot (graphique probabilité-probabilité) et le QQ plot (graphique quantile). Le code pour créer un PP plot peut être le suivant

> PP=function(Y,F=pnorm){
+   n=length(Y)
+   x=F(sort(Y))
+   y=seq(1/n/2,1-1/n/2,by=1/n)
+   return(data.frame(x=x,y=y))
+ }

qui représente (à un détail près) le nuage de points

https://latex.codecogs.com/gif.latex?\left\{F_0(x_{i:n}),\frac{i}{n}\right\}

et celui pour le QQ plot serait

> QQ=function(Y,Q=qnorm){
+   n=length(Y)
+   x=Q(seq(1/n/2,1-1/n/2,by=1/n))
+   y=sort(Y)
+   return(data.frame(x=x,y=y))
+ }

qui représente (toujours à un détail près) le nuage de points

https://latex.codecogs.com/gif.latex?%20\left\{F_0^{-1}\left(\frac{i}{n}\right),x_{i:n}\right\}

où https://latex.codecogs.com/gif.latex?%20F_0 est la loi que l’on cherche à tester, au sens où on a https://latex.codecogs.com/gif.latex?%20H_0:F=F_0 avec comme hypothèse alternative https://latex.codecogs.com/gif.latex?%20H_1:F\neq%20F_0.

Continue reading Tests non-paramétriques et simulations

Interpréter un test de Wald

Jeudi, dans le cadre du cours de statistique mathématiques, nous avions passé un peu de temps à mettre en oeuvre le test de Wald, et discuter la construction et l’interprétation d’un test. On reste ici sur un test paramètrique. La vraie valeur du paramètre, on va la fixer, pour générer des données

> p=.5

En l’occurence, on va considérer un modèle de Bernoulli, i.e. une succession de tirages de pièces (pile ou face)

> n=20
> set.seed(1)
> echantillon=sample(c("pile","face"),
+    size=n,replace=TRUE,
+    prob=c(p,1-p))
> echantillon
[1] "face" "face" "pile" "pile" "face" "pile" "pile" "pile" "pile" "face" "face" "face" "pile" "face" "pile" "face" "pile" "pile" "face" "pile"

Classiquement, on va recoder en https://latex.codecogs.com/gif.latex?\{0,1\}

> X=(echantillon=="pile")*1
> X
[1] 0 0 1 1 0 1 1 1 1 0 0 0 1 0 1 0 1 1 0 1

La vraisemblance est ici donnée par

> L=function(p) prod(dbinom(X,1,p))
> logL=function(p) sum(log(dbinom(X,1,p)))
> xp=seq(0,1,by=.01)
> yL=Vectorize(L)(xp)
> plot(xp,yL,type="l")

mais on va plutôt travailler sur la log-vraisemblance

> ylogL=Vectorize(logL)(xp)
> plot(xp,ylogL,type="l")

Continue reading Interpréter un test de Wald

Être – en même temps – moyen et extraordinaire

Juste avant 1850, Adolphe Quételet publiait Sur l’appréciation des documents statistiques, et en particulier sur l’application des moyennes, et introdusait alors le concept d’homme moyen (l’ouvrage est cité dans La théorie des moyennes et son emploi dans les sciences d’observation, de Joseph Lottin, paru en 1909 – je n’ai pas trouvé en ligne cet ouvrage de Quételet). Ce concept sera détaillé encore davantage dans Anthropométrie, ou mesure des différentes facultés de l’homme, concept qui sera abondament illustré

Mais si la moyenne est un concept fondamental en statistique, il convient de garder en mémoire que la moyenne n’existe peut être pas, au sein de la population. C’est ce que rappelle la fameuse ‘blague’ “un mec qui à la tête dans le four et les pieds dans le frigo, en moyenne, il est bien“. La moyenne peut être quelquechose d’extraordinaire (au sens premier du terme).

Continue reading Être – en même temps – moyen et extraordinaire

Où ont été observées les fortes hausses du prix de l’essence ?

Vendredi, dans mon billet “le prix de l’essence en cas de pénurie“, j’essayais de visualiser l’explosion du prix de l’essence dans certaines stations services, ces dernières semaines. Mais plus généralement, on peut se demander où ont été observées les fortes hausses.

A partir de la base téléchargeable sur http://www.prix-carburants.economie.gouv.fr/, il est possible de récupérer les informations les prix, par jour, dans toutes les stations essence de France. Pour ça, on utilise la fonction suivante,

> spatial=function(dt){
+  base=NULL
+  for(no in 1:length(l)){  
+    prix=list()
+    date=list()
+    j=0
+    for(i in 1:length(l[[no]])){
+    v=names(l[[no]])
+    if(!is.null(v[i])){
+    if(v[i]=="prix"){
+    j=j+1
+    date[[j]]=as.character(l[[no]][[i]]["maj"])
+    }}
+  }
+ n=j
+ D=as.Date(substr(unlist(date),1,10),"%Y-%m-%d")
+ k=which(D==D[which.max(D[D<=dt])])
+ if(length(k)>0){
+ B=Vectorize(function(i) l[[no]][[k[i]]])
+   (1:length(k))
+ if("nom" %in%  rownames(B)){  
+   k=which(B["nom",]=="Gazole")
+   prix=as.numeric(B["valeur",k])/1000
+   if(length(prix)==0) prix=NA
+   base1=data.frame(indice=no,
+ lat=as.numeric(l[[no]]$.attrs["latitude"])
+    /100000,
+ lon=as.numeric(l[[no]]$.attrs["longitude"])
+    /100000,
+ gaz=prix,
+ cp=l[[no]]$.attrs["cp"])
+ base=rbind(base,base1)
+ }}}
+ return(base)}

Par exemple, on peut récupérer les prix fin mai

> B1=spatial(as.Date("2016-05-31"))

et début mai,

> B2=spatial(as.Date("2016-05-01"))

En fusionnant ces deux bases, on peut récupérer la variation du prix, sur un mois,

> names(B1)=c("indice","lat","long","fin","cp")
> names(B2)=c("indice","lat","long","debut","cp")
> B=merge(B1,B2)
> B$var=(B$fin-B$debut)/B$debut*100

Plus particulièrement, on va s’intéresser aux stations dans lesquelles le prix a augmenté de plus de 10%,

> idx=which((B$lon>(-10))&(B$lon<20)&
+             (B$lat>35)&(B$lat<55))
> B=B[idx,]
> B=B[!is.na(B$var),]
> B$Y=(B$var>10)*1

Notons déjà que dans plus de 90% des stations, le prix a augmenté entre début mai, et fin mai

> mean(B$var > 0)
[1] 0.9288522

mais les hausses de plus de 10% restent rares (moins de 5% des stations essence)

> mean(B$var > 10)
[1] 0.03305452

La distribution de la variation des prix (en %) est la suivante

> plot(density(B$var),xlim=c(-5,20))

Si on regarde maintenant où ont été observées ces fortes hausses, on observe qu’elles peuvent avoir lieu n’importe où

> library(maps)
> map("france")
> points(B$lon,B$lat,pch=19,col=c(rgb(0,0,1,.25),
+ rgb(1,0,0,.95))[1+(B$var>10*1)],cex=.5)

et c’est encore plus flagrant quand on regarde les hausses de plus de 5%

> points(B$lon,B$lat,pch=19,col=c(rgb(0,0,1,.35),
+ rgb(1,0,0,.95))[1+(B$var>5*1)],cex=.5)

Et si on regarde les endroits où la proportion des stations ayant fortement augmenté les prix est la plus forte, on obtient la carte suivante,

Autrement dit la hausse des prix a pu être observée un peu partout en France. On peut alors comparer avec la carte des pénuries, mise en ligne sur https://mon-essence.fr/

Si certaines régions avec des fortes hausses de prix correspondent à des pénuries d’essence, ce n’est pas le cas partout. Bref, comprendre ce qui s’est passé sur les prix ce mois de mai est loin d’être simple… à suivre donc….

Le prix de l’essence en cas de pénurie

Cet après midi, sur le site de 60 millions de consommateurs, Lionel Maugain revenait sur les stations services qui avaient de très fortes variations de prix pendant les jours de pénurie.

En utilisant la base des prix des carburants à la pompe, téléchargeable sur http://www.prix-carburants.economie.gouv.fr/, il est possible de visualiser ces fortes variations de prix. Pour lire la base, rien de plus simple

> rm(list=ls())
> annee=year=2016
> fichier="PrixCarburants_annuel_2016.xml"
> library(plyr)
> library(XML)
> library(lubridate)
> l=xmlToList(fichier)

Ensuite, on transforme la grosse liste en une série temporelle, à l’aide de la fonction suivante

> time_series=function(no,type_gas="Gazole"){
+    prix=list()
+    date=list()
+    nom=list()
+    j=0
+    for(i in 1:length(l[[no]])){
+      v=names(l[[no]])
+      if(!is.null(v[i])){
+      if(v[i]=="prix"){
+      j=j+1
+      date[[j]]=as.character(l[[no]][[i]]["maj"])
+      prix[[j]]=as.character(l[[no]][[i]]["valeur"])
+      nom[[j]]=as.character(l[[no]][[i]]["nom"])
+     }}
+    }
+    id=which(unlist(nom)==type_gas)
+    n=length(id)
+    jour=function(j) as.Date(substr(date[[id[j]]],1,10),"%Y-%m-%d")
+    jour_heure=function(j) as.POSIXct(substr(date[[id[j]]],1,19), format = "%Y-%m-%d %H:%M:%S", tz = "UTC")
+    ext_y=function(j) substr(date[[id[j]]],1,4)
+    ext_m=function(j) substr(date[[id[j]]],6,7)
+    ext_d=function(j) substr(date[[id[j]]],9,10)
+    ext_h=function(j) substr(date[[id[j]]],12,13)
+    ext_mn=function(j) substr(date[[id[j]]],15,16)
+    prix_essence=function(i) as.numeric(prix[[id[i]]])/1000
+    base1=data.frame(indice=no,
+   id=l[[no]]$.attrs["id"],
+   adresse=l[[no]]$adresse,
+   ville=l[[no]]$ville,
+   lat=as.numeric(l[[no]]$.attrs["latitude"])
    /100000,
+   lon=as.numeric(l[[no]]$.attrs["longitude"])
    /100000,
+   cp=l[[no]]$.attrs["cp"],
+   saufjour=l[[no]]$ouverture["saufjour"], 
+   Y=unlist(lapply(1:n,ext_y)),
+   M=unlist(lapply(1:n,ext_m)),
+   D=unlist(lapply(1:n,ext_d)),
+   H=unlist(lapply(1:n,ext_h)),
+   MN=unlist(lapply(1:n,ext_mn)),
+   prix=unlist(lapply(1:n,prix_essence)))
+   base1=base1[!is.na(base1$prix),]
+   date_d=paste(year,"-01-01 12:00:00",sep="")
+   date_f=paste(year,"-06-01 12:00:00",sep="")
+   vecteur_date=seq(as.POSIXct(date_d, format =
+   "%Y-%m-%d %H:%M:%S"),
+   as.POSIXct(date_f, format = 
+   "%Y-%m-%d %H:%M:%S"),by="days")
+   date=paste(base1$Y,"-",base1$M,"-",base1$D,
+   " ",base1$H,":",base1$MN,":00",sep="")
+   date_base=as.POSIXct(date, format = 
+   "%Y-%m-%d %H:%M:%S", tz = "UTC")
+   idx=function(t) sum(vecteur_date[t]>=date_base)
+      vect_idx=Vectorize(idx)(1:length(vecteur_date))
+      P=c(NA,base1$prix)
+      vp=P[1+vect_idx]
+      vp=vp*100/(vp[!is.na(vp)][1])
+      base2=ts(vp,start=year,frequency=365)
+      list(base=base1,ts=base2)
+    }

Dans l’article, la plus forte variation a été observée à Sevran. Pour récupérer la liste des stations services, on utilise

> ville=list()
> for(no in 1:length(l)) 
      ville[[no]]=l[[no]]$ville
> ville=unlist(ville)
> idx=which(ville=="Sevran")

Pour la première station, l’évolution du prix est relativement modérée

> plot(time_series(idx[1])$ts,col="red",type="l")
> abline(h=100,lty=2)

Même si le prix semble avoir augmenté de 5% pendant les dernières semaines. Mais ce n’est effectivement rien en comparaison de l’autre station,

> plot(time_series(idx[2])$ts,col="red",type="l")
> abline(h=100,lty=2)

avec un prix qui a augmenté de 50% en 2 mois, et presque 30% en quelques jours.

A Villeneuve-d’Asq, on retrouve un comportement similaire, mais moins important

> idx=which(ville=="Villeneuve-d'Ascq")
> plot(time_series(idx[3])$ts,col="red",type="l")
> abline(h=100,lty=2)