Tag Archives: Paris

Douce France

Hier, sur Twitter, je mentionnais une carte parue sur Times Higher Education à propos des “100 meilleurs universités” (dans le monde entier). Je ne sais plus par quel biais j’ai vu des cartes, par pays, permettant de voir où étaient les 100 meilleures universités (d’ordinaire, on a un classement, mais là, pour la première fois, on avait des cartes),

Le premier tweet, en anglais, qui montrait la carte est passé assez inaperçu. Ce n’est pas le cas du second, en français, qui a été beaucoup repris, dans lequel je posais la question de manière plus directe, “sauras-tu remarquer une particularité géographique de l’enseignement supérieur français?“. Car je trouve la carte française assez choquante. disait de manière assez élégante que la “ville lumière” faisait de l’ombre aux autres villes françaises. Ce n’est pas peu dire ! Dans tous les autres pays comparables, avec 4, 5 ou 6 universités dans ce classement, on a une relative dispersion sur le territoire. Mais en France, on a Paris. Et c’est tout.

C’est amusant car cette semaine, j’ai commencé le livre Unknown Quantity de John Derbyshire qui raconte, de manière très vulgarisée (et donc parfaite pour moi) l’histoire de l’algèbre. Et dans un des premiers chapitres (on est encore en Mésopotamie) il y a ce paragraphe qui m’a beaucoup rappelé cette carte

Effectivement, il faudrait que je creuse un peu, mais j’ai du mal à voir ce qui va ressortir de cette méga-concentration géographique typiquement française de l’enseignement supérieur. Rien de bon, je crains…

Crowded Cities, Paris, Hong Kong and Montréal

Over the past years, I’ve been living in different cities, all of them being completely different, compared with the others. I have been living in Paris, which is a big city in Europe, with a large suburban neighborhood, too (la banlieue).

Then, I’ve been living in Hong Kong, which is a larger city, in Asia.

It was crowded. I mean, it was the feeling I had, while I was living there. And more recently, I’ve been living in Montréal, in North America. Montreal is a large city. Or to be more specific, an island,

The three cities are quite different. Paris, 2.211 million unhabitants, and 105,4 km² (density 21,057 unhabitants per km²). Montréal, 1.621 million unhabitants, and three times wider 365.1 km² (density 4,441 unhabitants per km²). Hong Kong, 7.234 million unhabitants, and again three times wider 1,104 km² (density 6,553 unhabitants per km²). In Hong Kong, there are several hill where it is not possible to build anything: on a large part of the island, the density is null.

Continue reading Crowded Cities, Paris, Hong Kong and Montréal

Short (but Buzy) Trip in Europe

I will flight tomorrow evening to go to Belgium. I will attend a PhD defense this Monday, at Louvain-la-Neuve. On Tuesday, I will give a crash course on risk measures, for PhD Students, in Louvain-la-Neuve. Then, Wednesday morning, I have a meeting with some friends in an insurance company, in Brussels, and I will be around noon in Paris, to give a talk on Big Data, for the Big Data Task Force of the (French) Institute of Actuaries. I will then have some meetings in the afternoon. I should flight back on Thursday, and give my lecture on Predictive Modeling in the evening…

I guess that Saturday will be a mix between two very different activities, enjoying some time with the kids (and watch the game), and trying to finalize my slides…

Le Vote par Procuration en France

La Vie des Idées a mis en ligne, ce matin, un court texte, écrit par Baptiste Coulmont (a.k.a. @coulmont) et Joël Gombin (a.k.a. @joelgombin), auquel j’ai très modestement contribué, intitulé “Un homme, deux voix. Le vote par procuration“.

Alors que sur son blog, Baptiste a rajouté pas mal d’information sur le vote par procuration en France (et le contexte général, en particulier pourquoi autant de partis courtisent certaines personnes en les incitant à voter par procuration), et sur les bases de données, je voulais en profiter pour mettre en ligne quelques codes utilisés dans l’article, et en particulier, mentionner des graphiques non-utilisés car plus difficile à interpréter, mais à mon avis plus juste en terme de modèle (comme les conclusions étaient les mêmes, on a retenu des graphiques plus classiques). Rappelons tout d’abord qu’on analyse non pas le vote à partir de données individuelles (ceci ne peut s’obtenir, le vote étant encore secret en France), mais à partir des résultats des différents bureaux de vote (c’est la notion de corrélation écologique évoquée dans le texte, à cause du problème potentiel d’ecological fallacy). Moyennant toutes ces précautions d’usage, on a essayé d’analyser les données à notre disposition.

  • Modèle de régression, et recherche de variables explications

Pour faire une régression, et expliquer le taux de procurations dans un bureau de vote, ma première idée était de dire que  où  est le nombre de procurations dans le bureau de vote , et où  est (au choix) le nombre d’électeurs inscrits ou le nombre d’électeurs ayant voté. On suppose ici que , la proportion d’électeurs qui a voté par procuration peut être fonction de divers variables explicatives. Les variables, elles sont dans la base suivante (je renvoie vers le blog de Baptiste pour les bases que l’on utilise, en particulier à partir des données de insee.fr et d’opendata.paris.fr)

> bt1=read.table("paris2007-pres-t1.csv",header=TRUE,sep=";")
> bt2=read.table("paris2007-pres-t2.csv",header=TRUE,sep=";")
> bv=read.table("paris-bv-insee-07.csv",header=TRUE,sep=";")
> bv$BV=bv$BVCOM
> baset1=merge(bt1,bv,by="BV")
> baset2=merge(bt2,bv,by="BV")
> baset1$LOGEMENT=baset1$PROPRIO+baset1$LOCNONHLM+baset1$LOCHLM+baset1$GRATUIT
> baset2$LOGEMENT=baset2$PROPRIO+baset2$LOCNONHLM+baset2$LOCHLM+baset2$GRATUIT

Si on suppose que  est fonction  le taux de logements occupés par leur propriétaire, dans le quartier (associé à un bureau de vote),

> variable="PROPRIO"
> reference="LOGEMENT"
> baset1$taux=baset1[,variable]/baset1[,reference]
> baset2$taux=baset2[,variable]/baset2[,reference]

il est légitime de tenter une régression logistique,

voire un lissage par splines, si on pense que le lien peut ne pas être linéaire,

Ceci se fait à l’aide du code suivant, pour la version lissée (par splines cubiques)

> b=hist(baset1$taux,plot=FALSE)
> library(splines)
> regt1=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset1)
> regt2=glm(PROCURATIONS/INSCRITS~bs(taux,6),family=binomial,weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL[4])

et, pour la régression logistique standard (linéaire)

> lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2)
> regt1l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset1)
> regt2l=glm(PROCURATIONS/INSCRITS~taux,family=binomial,weights=INSCRITS,data=baset2)
> ND=data.frame(taux=ug)
> pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response")
> pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response")
> lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2)
> legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)],
+ lwd=2,lty=1,border=NA)

(en rajoutant une petite légende, avec une visualisation pour les deux tours de l’élection présidentielle, avec un intervalle de confiance sur la prévision de mon taux de procuration).

On peut faire la même chose sur le taux de logement HLM, dans le quartier,

Les dessins sont parlant, mais dans la sortie du modèle de régression, l’interprétation de  laisse à désirer (toute suggestion est la bienvenue !).

> summary(regt1l)

Call:
glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial, 
    data = baset1, weights = INSCRITS)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-12.9549   -1.5722    0.0319    1.6292   13.1303  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.70811    0.01516  -244.6   <2e-16 ***
taux         1.49666    0.04012    37.3   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 12507  on 836  degrees of freedom
Residual deviance: 11065  on 835  degrees of freedom
AIC: 15699

Number of Fisher Scoring iterations: 4

> summary(regt2l)

Call:
glm(formula = PROCURATIONS/INSCRITS ~ taux, family = binomial, 
    data = baset2, weights = INSCRITS)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-15.4872   -1.7817   -0.1615    1.6035   12.5596  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -3.24272    0.01230 -263.61   <2e-16 ***
taux         1.45816    0.03266   44.65   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9424.7  on 836  degrees of freedom
Residual deviance: 7362.3  on 835  degrees of freedom
AIC: 12531

Number of Fisher Scoring iterations: 4

 

On a alors voulu comparer avec un modèle qui me semble moins juste, mais qui est plus simple à interpréter, où on suppose que le taux de procuration (par bureau) est expliqué par un modèle linéaire

(que l’on peut aussi lisser, pour vérifier que le lien est effectivement linéaire). Le code est ici

> regt1=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset1)
> regt2=lm(PROCURATIONS/INSCRITS~bs(taux,6),weights=INSCRITS,data=baset2)
> u=seq(min(baset1$taux)+.015,max(baset1$taux)-.015,by=.001)
> ND=data.frame(taux=u)
> ug=seq(0,max(baset1$taux)+.05,by=.001)
> pt1=predict(regt1,newdata=ND,se=TRUE,type="response")
> pt2=predict(regt2,newdata=ND,se=TRUE,type="response")
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> plot(ug,ug*1,col="white",xlab=nom,ylab="Taux de procuration",
+ ylim=c(0,.1))
> for(i in 1:(length(b$breaks)-1)){
+ polygon(b$breaks[i+c(0,0,1,1)],c(0,b$counts[i],b$counts[i],0)
+ /max(b$counts)*.05,col="light yellow",border=NA)}
> polygon(c(u,rev(u)),c(pt1$fit+2*pt1$se.fit,rev(pt1$fit-2*pt1$se.fit)),
+ border=NA,density=30,col=CL[4])
> lines(u,pt1$fit,col=CL[6],lwd=2)
> polygon(c(u,rev(u)),c(pt2$fit+2*pt2$se.fit,rev(pt2$fit-2*pt2$se.fit)),
+ border=NA,density=30,col=CL[3])
> lines(u,pt2$fit,col=CL[1],lwd=2)
> regt1l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset1)
> regt2l=lm(PROCURATIONS/INSCRITS~taux,weights=INSCRITS,data=baset2)
> ND=data.frame(taux=ug)
> pt1l=predict(regt1l,newdata=ND,se=TRUE,type="response")
> pt2l=predict(regt2l,newdata=ND,se=TRUE,type="response")
> lines(ug,pt1l$fit,col=CL[5],lty=2)
> lines(ug,pt2l$fit,col=CL[2],lty=2)
> legend(0,.1,c("Second Tour","Premier Tour"),col=CL[c(1,6)],
+ lwd=2,lty=1,border=NA)

(j’ai tout mis d’un coup cette fois, les modèles lissés et linéaires, l’un à la suite de l’autre)

Cette fois, on a une sortie de régression plus classique,

> summary(regt1l)

Call:
lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset1, weights = INSCRITS)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-1.9994 -0.2926  0.0011  0.3173  3.2072 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.021268   0.001739   12.23   <2e-16 ***
taux        0.054371   0.004812   11.30   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.646 on 835 degrees of freedom
Multiple R-squared:  0.1326,	Adjusted R-squared:  0.1316 
F-statistic: 127.7 on 1 and 835 DF,  p-value: < 2.2e-16

> summary(regt2l)

Call:
lm(formula = PROCURATIONS/INSCRITS ~ taux, data = baset2, weights = INSCRITS)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-2.9029 -0.4148 -0.0338  0.4029  3.4907 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.033909   0.001866   18.17   <2e-16 ***
taux        0.079749   0.005165   15.44   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.6934 on 835 degrees of freedom
Multiple R-squared:  0.2221,	Adjusted R-squared:  0.2212 
F-statistic: 238.4 on 1 and 835 DF,  p-value: < 2.2e-16

On note plusieurs choses de ces graphiques. (i) les deux types de régression donnent des modèles pour les taux de procuration très très proches. Donc autant prendre le plus simple à interpréter. (ii) le lissage n’apporte rien, et le modèle linéaire semble pertinent. On a ainsi regardé plusieurs variables, et on en a retenu un certain nombre, pour tenter un modèle multiple. Avant de parler des résidus de notre modèle, je devrais peut être prendre quelques lignes pour parler un peu de cartographie (une nouvelle fois, Baptiste m’a fait découvrir de belles fonctions).

  • Visualiser les bureaux de vote à Paris

Pour récupérer le fond de carte, avec les bureaux de vote, on utilise (là encore, je renvoie au blog de Baptiste, qui explique l’utilisation de données de cartelec.net)

> library(maptools)
> library(rgdal)
> library(classInt)
> paris=readShapeSpatial("paris-cartelec.shp")

Si on veut visualiser, par exemple, le taux de procuration (disons la moyenne entre les deux tours), on utilise les données suivantes

> elec=data.frame()
> elec=cbind(bt1$BV,(bt1$PROCURATIONS+bt2$PROCURATIONS),(bt1$EXPRIMES+bt2$EXPRIMES))
> colnames(elec)=c("BV","PROCURATIONS","EXPRIMES")
> elec=as.data.frame(elec)
> elec$BV=bt1$BV

Ensuite, viennent les fonctions graphiques, où on va passer d’un taux à une classe et d’une classe à une couleur,

> m=match(paris$BUREAU,elec$BV)
> plotvar=100*elec$PROCURATIONS/elec$EXPRIMES
> nclr=7
> plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] 
> class=classIntervals(plotvar[m], nclr, style="fisher",dataPrecision=1)
> colcode=findColours(class, plotclr)

Reste à conclure, en faisant une visualisation graphique de nos données

> par(mar=c(1,1,1,1))
> plot(paris,col=colcode,border=colcode)
> legend(656274.9, 6867308,legend=names(attr(colcode,"table")), 
+ fill=attr(colcode, "palette"), cex=1, bty="n",
+ title="Frequence procurations (%)")

Histoire de conclure, on peut regarder un peu nos résidus, obtenus sur un modèle linéaire. Considérons un modèle sur seulement trois variables explicatives,

> regt1=lm(PROCURATIONS/INSCRITS~I(POP65P/POP)+
+ I(PROPRIO/LOGEMENT)+I(CS3/POP1564),weights=INSCRITS,data=baset1)

Dans ce cas, la visualisation des résidus donne

> m=match(paris$BUREAU,elec$BV)
> plotvar=100*residuals(regt1)
> nclr=7
> plotclr=brewer.pal(nclr,"RdYlBu")[nclr:1] 
> class=classIntervals(plotvar[m], nclr, style="fisher",dataPrecision=1)
> colcode=findColours(class, plotclr)
> par(mar=c(1,1,1,1))
> plot(paris,col=colcode,border=colcode)
> legend(656274.9, 6867308,legend=names(attr(colcode,"table")), 
+ fill=attr(colcode, "palette"), cex=1, bty="n",title="Residus")

Idéalement, il faudrait avoir un beau bruit (spatial), c’est à dire avoir des couleurs réparties aléatoirement sur Paris. Il reste encore pas mal de régions dont les voisins sont de la même couleurs, et on repère quelques quartiers atypiques, avec soit des résidus importants négativement, ou positivement. Comme toujours, en modélisation, on pourrait passer des heures pour essayer de capturer tous les effets, mais aucun modèle avec les variables à notre disposition nous a permis de faire réellement mieux.

Playing with robots

My son would be extremely proud if I tell him I can spend hours building robots. Well, my robots are not as fancy as Dr Tenma’s, but they usually do what I ask them to do. For instance, it is extremely simple to build a robot with R, to extract data from websites. I have mentioned it here (one tennis matches), but it failed there (on NY Marathon). To illustrate the use of robots, assume that one wants to build his own dataset to study prices of airline tickets. First, we have to choose a departure city (e.g. Paris) and an arrival city (e.g. Montreal). Then, one wants to look at all possible dates from April first (I ran it last month) till the end of December (so we create a vector with all leaving dates, namely a vector for the day, one for the month, and one for the year). Then, we choose a return date (say 3 days after).

DEP="Paris"
ARR="Montreal"
DATE1D=rep(c(1:30,1:31,1:30,1:31,1:31,1:30,1:31,1:30,
1:31,1:31,1:29),3)
DATE1M=rep(c(rep(4,30),rep(5,31),rep(6,30),rep(7,31),
rep(8,31),rep(9,30),rep(10,31),rep(11,30),rep(12,31),
rep(1,31),rep(2,29)),3)
DATE1Y=rep(c(rep(2011,30+31+30+31+31+30+31+
30+31+31+28),rep(2012,31+29)),3)
k=3
DATE3D=c((1+k):30,1:31,1:30,1:31,1:31,1:30,1:31,
1:30,1:31,1:31,1:29,1:k)
DATE3M=c(rep(4,30-k),rep(5,31),rep(6,30),rep(7,31),rep(8,31),
rep(9,30),rep(10,31),rep(11,30),rep(12,31),rep(1,31),rep(2,29),
rep(3,k))
DATE3Y=c(rep(2011,30+31+30+31+31+30+31+30+31+
31+28-k),re
p(2012,31+29+k))

It is also possible (for a nice robot), to skip all prior dates

skip=max(as.numeric(Sys.Date()-as.Date("2011-04-01")),1)

Then, we need a website where requests can be written nicely (with cities and dates appearing explicitly). Here, I cannot not mention the website that I used since it is stated on the website that it is strictly forbidden to run automatic requests… Anyway, consider a loop create a url address (actually I chose the value of the date randomly, since I had been told that those websites had memory: if you ask too many times for the same thing during a short period of time, prices would go up),

URL=paste("http://www.♦♦♦♦/dest.dll?qscr=fx&flag=q&city1=",
DEP,"&citd1=",ARR,"&",
"date1=",DATE1D[s],"/",DATE1M[s],"/",DATE1Y[s],
"&date2=",DATE3D[s],"/",DATE3M[s],"/",DATE3Y[s],
"&cADULT=1",sep="")

then, we just have to scan the webpage, looking for ticket prices (just looking for some specific names)

page=as.character(scan(URL,what="character"))
I=which(page%in%c("Price0","Price1","Price2"))
if(length(I)>0){
PRIX=substr(page[I+1],2,nchar(page[I+1]))
if(PRIX[1]=="1"){PRIX=paste(PRIX,page[I+2],sep="")}
if(PRIX[1]=="2"){PRIX=paste(PRIX,page[I+2],sep="")}

Here, we have to be a bit cautious, if prices exceed 1000. Then, it is possible to start a statistical study. For instance, if we compare to destination (from Paris), e.g. Montréal and New York, we obtain the following patterns (with high prices during holidays),

It is also possible to run the code twice (here it was run last month, and a couple of days ago), for the same destination (from Paris to Montréal),

Of course, it would be great if I could run that code say every week, to build up a nice dataset, and to study the dynamic of prices…

The problem is that it is forbidden to do this. In fact, on the website, it is mentioned that if we want to extract data (for an academic purpose), it is possible to ask for an extraction. But if we do tell that we study specific prices, data might be biased. So the good idea would be to use several servers, to make several requests, randomly, and to collect them (changing dates and destination). But here, my computing skills – unfortunately – reach a limit….

Warming in Paris: minimas versus maximas ?

Recently, I received comments (here and on Twitter) about my previous graphs on the temperature in Paris. I mentioned in a comment (there) that studying extremas (and more generally quantiles or interquantile evolution) is not the same as studying the variance. Since I am not a big fan of the variance, let us talk a little bit about extrema behaviour.

In order to study the average temperature it is natural to look at the linear (assuming that it is linear, but I proved that it could reasonably be assumed as linear in the paper) regression, i.e. least square regression, which gives the expected value. But if we care about extremes, or almost extremes, it is natural to look at quantile regression.

For instance, below, the green line is the least square regression, the red one is 97.5% quantile, and the blue on the 2.5% quantile regression.

It looks like the slope is the same, i.e. extremas are increasing as fast as the average…

tmaxparis=read.table("temperature/TG_SOUID100845.txt",
skip=20,sep=",",header=TRUE)
head(tmaxparis)
Dparis=as.Date(as.character(tmaxparis$DATE),"%Y%m%d")
Tparis=as.numeric(tmaxparis$TG)/10
Tparis[Tparis==-999.9]=NA
I=sample(1:length(Tparis),size=5000,replace=FALSE)
plot(Dparis[I],Tparis[I],col="grey")
abline(lm(Tparis~Dparis),col="green")
library(quantreg)
abline(rq(Tparis~Dparis,tau=.025),col="blue")
abline(rq(Tparis~Dparis,tau=.975),col="red")

(here I plot randomly some points to avoid a too heavy figure, since I have too many observations, but I keep all the observations in the regression !).

Now, if we look at the slope for different quantile level (Fig 6 in the paper, here, but on minimum daily temperature, here I look at average daily temperature), the interpretation is different.

s=0
COEF=SD=rep(NA,199)
for(i in seq(.005,.995,by=.005)){
s=s+1
REG=rq(Tparis~Dparis,tau=i)
COEF[s]=REG$coefficients[2]
SD[s]=summary(REG)$coefficients[2,2]
}

with the following graph below,

s=0
plot(seq(.005,.995,by=.005),COEF,type="l",ylim=c(0.00002,.00008))
for(i in seq(.005,.995,by=.005)){
s=s+1
segments(i,COEF[s]-2*SD[s],i,COEF[s]+2*SD[s],col="grey")
}
REG=lm(Tparis~Dparis)
COEFlm=REG$coefficients[2]
SDlm=summary(REG)$coefficients[2,2]
abline(h=COEFlm,col="red")
abline(h=COEFlm-2*SDlm,lty=2,lw=.6,col="red")
abline(h=COEFlm+2*SDlm,lty=2,lw=.6,col="red")

Here, for minimas (quantiles associated to low probabilities, on the left), the trend has a higher slope than the average, so in some sense, warming of minimas is stronger than average temperature, and on other hand, for maximas (high probabilities on the right), the slope is smaller – but positive – so summer are warmer, but not as much as winters.
Note also that the story is different for minimal temperature (mentioned in the paper) compared with that study, made here on average daily temperature (see comments)… This is not a major breakthrough in climate research, but this is all I got…

More climate extremes, or simply global warming ?

In the paper on the heat wave in Paris (mentioned here) I discussed changes in the distribution of temperature (and autocorrelation of the time series).

During the workshop on Statistical Methods for Meteorology and Climate Change today (here) I observed that it was still an important question: is climate change affecting only averages, or does it have an impact on extremes ? And since I’ve seen nice slides to illustrate that question, I decided to play again with my dataset to see what could be said about temperature in Paris.
Recall that data can be downloaded here (daily temperature of the XXth century).

tmaxparis=read.table("/temperature/TX_SOUID100124.txt",
skip=20,sep=",",header=TRUE)
Dmaxparis=as.Date(as.character(tmaxparis$DATE),"%Y%m%d")
Tmaxparis=as.numeric(tmaxparis$TX)/10
tminparis=read.table("/temperature/TN_SOUID100123.txt",
skip=20,sep=",",header=TRUE)
Dminparis=as.Date(as.character(tminparis$DATE),"%Y%m%d")
Tminparis=as.numeric(tminparis$TN)/10
Tminparis[Tminparis==-999.9]=NA
Tmaxparis[Tmaxparis==-999.9]=NA
annee=trunc(tminparis$DATE/10000)
MIN=tapply(Tminparis,annee,min)
plot(unique(annee),MIN,col="blue",ylim=c(-15,40),xlim=c(1900,2000))
abline(lm(MIN~unique(annee)),col="blue")
abline(lm(Tminparis~unique(Dminparis)),col="blue",lty=2)
annee=trunc(tmaxparis$DATE/10000)
MAX=tapply(Tmaxparis,annee,max)
points(unique(annee),MAX,col="red")
abline(lm(MAX~unique(annee)),col="red")
abline(lm(Tmaxparis~unique(Dmaxparis)),col="red",lty=2)

On the plot below, the dots in red are the annual maximum temperatures, while the dots in blue are the annual minimum temperature. The plain line is the regression line (based on the annual max/min), and the dotted lines represent the average maximum/minimum daily temperature (to illustrate the global tendency),

It is also possible to look at annual boxplot, and to focus either on minimas, or on maximas.

annee=trunc(tminparis$DATE/10000)
boxplot(Tminparis~as.factor(annee),ylim=c(-15,10),
xlab="Year",ylab="Temperature",col="blue")
x=boxplot(Tminparis~as.factor(annee),plot=FALSE)
xx=1:length(unique(annee))
points(xx,x$stats[1,],pch=19,col="blue")
abline(lm(x$stats[1,]~xx),col="blue")
annee=trunc(tmaxparis$DATE/10000)
boxplot(Tmaxparis~as.factor(annee),ylim=c(15,40),
xlab="Year",ylab="Temperature",col="red")
x=boxplot(Tmaxparis~as.factor(annee),plot=FALSE)
xx=1:length(unique(annee))
points(xx,x$stats[5,],pch=19,col="red")
abline(lm(x$stats[5,]~xx),col="red")

Plain dots are average temperature below the 5% quantile for minima, or over the 95% quantile for maxima (again with the regression line),

We can observe an increasing trend on the minimas, but not on the maximas !
Finally, an alternative is to remember that we focus on annual maximas and minimas. Thus, Fisher and Tippett theory (mentioned here) can be used. Here, we fit a GEV distribution on a blog of 10 consecutive years. Recall that the GEV distribution is

http://freakonometrics.blog.free.fr/public/perso/gev1.png
install.packages("evir")
library(evir)
Pmin=Dmin=Pmax=Dmax=matrix(NA,10,3)
for(s in 1:10){
X=MIN[1:10+(s-1)*10]
FIT=gev(-X)
Pmin[s,]=FIT$par.ests
Dmin[s,]=FIT$par.ses
X=MAX[1:10+(s-1)*10]
FIT=gev(X)
Pmax[s,]=FIT$par.ests
Dmax[s,]=FIT$par.ses
}

The location parameter http://freakonometrics.blog.free.fr/public/perso/gev4.png is the following, with on the left the minimas and on the right the maximas,

while the scale parameter http://freakonometrics.blog.free.fr/public/perso/gev3.png is

and finally the shape parameter http://freakonometrics.blog.free.fr/public/perso/gev2.png is

On those graphs, it is very difficult to say anything regarding changes in temperature extremes… And I guess this is a reason why there is still active research on that area…

Des étés pluvieux en Bretagne ? une réalité statistique…

Pour compléter le précédant billet (ici) on peut se demander en quoi la Bretagne est différente des autres régions françaises… Nous avions vu ici le niveau de précipitation moyen, jour après jours pendant les mois d’été, en Bretagne. A Rennes. En revanche, à Paris on obtient la moyenne suivante,

que l’on peut comparer à Marseille,

ou encore à Strasbourg,

Sur la figure ci-dessous, on voit que la probabilité d’avoir de la pluie à Paris (au sens au moins 0.1 mm d’eau dans la journée, en trait gras bleu, au moins 2 mm d’eau dans la journée, en trait bleu) est supérieure à la probabilité d’avoir de la pluie à Rennes (respectivement en bleu clair gras, et en bleu clair fin)

On est certes très au dessus de Marseille,

mais très en dessous de Strasbourg,

Mais au delà des lois marginales, ces villes sont différentes de la Bretagne si l’on regarde les matrices de transition.

  • Transition d’un jour sur l’autre

Pour Rennes, si on regarde jour après jour, on obtient

jour 
beau
temps
pluie
jour 
beau temps
1955 612 2567
pluie
606 723 1329

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
76,15 % 23,85 %
pluie
45,60 % 54,40 %

Pour Paris, la probabilité de transition jour après jour a la forme suivante

jour 
beau
temps
pluie
jour 
beau temps
2689 959 3648
pluie
946 1466 2412

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
73,71 % 26,29 %
pluie
39,22 % 60,78 %

Pour Marseille, la probabilité de transition jour après jour a la forme suivante

jour 
beau
temps
pluie
jour 
beau temps
2527 375 2902
pluie
362 216 578

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
87,08 % 12,92 %
pluie
62,63 % 37,37 %

Pour Strasbourg, la probabilité de transition jour après jour a la forme suivante

jour 
beau
temps
pluie
jour 
beau temps
31 128 159
pluie
132 1464 1596

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
19,50 % 80,50 %
pluie
8,27 % 91,73 %
  • Transition d’une semaine sur l’autre

Si en revanche on regarde les matrices de transition semaine par semaine, on a des résultats assez différents. Une bonne semaine signifie aucun jour avec plus de 2 mm de pluie.
Pour Rennes, si on regarde semaine après semaine

jour 
beau
temps
pluie
jour 
beau temps
379 25 404
pluie
26 7 33

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
93,81 % 6,19 %
pluie
78,79 % 21,21 %

Pour Paris, la probabilité de transition semaine après semaine a la forme suivante

jour 
beau
temps
pluie
jour 
beau temps
576 46 622
pluie
53 4 57

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
92,60 % 7,40 %
pluie
92,98 % 7,02 %

Pour Marseille, la probabilité de transition semaine après semaine a la forme suivante

jour 
beau
temps
pluie
jour 
beau temps
274 59 333
pluie
47 9 56

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
82,28 % 17,72 %
pluie
83,93 % 16,07 %

Pour Strasbourg, la probabilité de transition semaine après semaine a la forme suivante

jour 
beau
temps
pluie
jour 
beau temps
1494 614 2018
pluie
613 939 1552

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
70,87 % 29,13 %
pluie
39,50 % 60,50 %

Les tests du chi deux, d’indépendance d’une semaine sur l’autre donnent

  • à Rennes, une statistique du chi-deux de 8,054, soit une p-value de 0,45%
  • à Paris, une statistique du chi-deux de 0,025, soit une p-value de 87,26%
  • à Marseille, une statistique du chi-deux de 0,012, soit une p-value de 91,24%
  • à Strasbourg, une statistique du chi-deux de 0,7649, soit une p-value de 38,18%

autrement dit l’hypothèse d’indépendance est acceptée partout, sauf à Rennes….

  • Moralité ?

De manière assez paradoxale, on prétend que la Bretagne a un temps changeant, et pour reprendre le titre du précédant billet, effectivement, en Bretagne, il peut faire beau plusieurs fois par jour. Mais sur le long terme, d’une semaine sur l’autre, le temps est au contraire très corrélé, contrairement aux autres régions. A Paris, Marseille ou Strasbourg, qu’il ait fait beau ou qu’il ait plus la semaine précédente, cela n’apporte aucune information sur la probabilité d’avoir de la pluie la semaine où l’on vient en vacances…. Mais pas en Bretagne: manifestement, il existe donc des étés pourris, où il pourra pleuvoir toutes les semaines, et des étés superbes où il ne pleut jamais….