Tag Archives: 3wen

La Ville, Conférence CartoStats

Vendredi, je participerai à la journée CartoStats organisée à Paris Diderot à l’occasion des 60 ans du Comité Français de Cartographie par Françoise Bahoken, Laurent Jégou et Hélène Mathian (entre autres). On m’a proposé une intervention de 20 minutes, et ça sera l’occasion de revenir sur des travaux récents avec Lucas Vernet (doctorant à Sciences Po), Alfred Galichon (professeur à Columbia) et Ewen Gallic (qui vient de défendre son doctorat, la semaine dernière, et qui commencera son post-doc avec moi à la rentrée sur la démographie participative). Et histoire d’avoir une cohérence dans mon exposé, j’ai pris un fil rouge qui sera la ville en économie. Les slides sont en ligne (ils ont \’et\’e mis \`a jour par rapport au billet, avec des r\’ef\’erences compl\’ementaires), mais ils sont purement visuels, et je me suis dit qu’un billet permettrait de les compléter.

Plus que la ville, c’est la notion d’agglomération autour d’un centre qui va m’intéresser – ce que les géographes appellent CBD (Central Business District) ou centre des affaires – et la localisation des vendeurs dans l’espace (souvent dans des “villes”) compte tenu d’une répartition d’acheteurs. On peut remonter aux travaux de von Thünen ou Kohl qui, les premiers, justifient l’existence de villes monocentriques, autour d’un centre.

Mais pour tous les économistes, quand on parle de “ville”, on pense instinctivement à la “ville linéaire” introduite par Hotelling en 1929

avec cette ville incroyablement simpliste, on peut justifier la concentration des firmes en un “centre”. On suppose  ici la population répartie (uniformément) dans cette ville, et les firmes cherchent où se placer, de manière optimale. Et on peut montrer assez facilement qu’avec deux firmes, un équilibre stable est que les villes s’installent toutes les deux au centre.

Mais on l’a dit, cette ville linéaire n’est pas très réaliste, et pas très représentative de la plupart des villes. Cela n’empêche pas certains auteurs de changer les hypothèses pour arriver à d’autres configurations d’équilibres, comme le montre Fujita & Thisse (1997), en reprenant Fujita & Ogawa (1982)

La première tentative pour passer en bi-dimensionnel a été la ville circulaire de Salop (1979), mais on le voit, c’est une fausse extension du cas précédant. On reste assez fondamentalement en une dimension,

Une autre extension pourrait être la ville “étoilée”, avec une population répartie en étoile, le long de routes. Mais là encore, si ce modèle est encore assez populaire, on se demande pourquoi ne pas envisager la création de route périphérique, autour du centre.

Autrement dit, l’économie spatiale et l’économie urbaine ne peuvent se penser sans faire intervenir la dimension des transports (comme le dit le théorème d’impossibilité de Starrett (1978)). L’extension logique est alors de voir la ville sous forme d’un réseau, comme l’avait fait Christaller (1933), afin de commencer une analyse microéconomique de la répartition de l’offre et de la demande.

Mais cette dimension rend le problème plus complexe, car il existe une infinité de réseaux, et tous ont leurs particularité. Et les réseaux ont une particularité : l’intuition que l’on a en terme de distance et de voisins tend à disparaître. Par exemple sur le graphique ci-dessous, 8 points sont répartis sur la région parisienne, et on y accède par métro. La figure du haut correspond au diagramme de Voronoï classique (c’est à dire aux zones d’affluence de magasins si les prix sont identiques) avec une norme euclidienne (à vol d’oiseau). Le graphique du dessous correspond au diagramme de Voronoï si on utilise la distance en métro. On voit que les zones d’affluence ne sont plus les mêmes.

Le problème économique que l’on peut envisager est alors celui d’un matching bipartite, où les “consommateurs” vont se rendre chez un “producteur” (et un seul). On peut imaginer des contraintes de capacité des producteurs, voire des problèmes de congestion sur le réseau, etc.

On a récemment étudié ce genre de modèles avec Alfred et Lucas, tout d’abord avec la localisation des maternités (le bien est ici un accouchement : on sait – par communes – le nombre de femme ayant accouché, et on a des informations sur les maternités) et les stations essence (le bien est ici un plein d’essence : on a la localisation des stations en Bretagne, et la localisation d’utilisateurs d’une api qui souhaitent connaître le prix de l’essence proche de leur localisation). Compte tenu de l’offre et de la demande, il est possible d’en déduire le “prix” de ces biens, et donc leur valeur. Ceci permet, par exemple, de voir où il serait nécessaire d’implanter une nouvelle maternité, par exemple.

Mais l’étude de réseaux ne peut se faire en statistique. L’important est d’ajouter une dimension temporelle. Dans l’exemple de la ville linéaire, on suppose les agents uniformément répartis dans la ville, puis les firmes se localisent, mais comment ne pas imaginer ensuite que cette localisation influencera a son tour la répartition des agents? Beckman (1976) note ainsi que les individus se répartiront au sein de la ville avec “a density distribution which is symmetric and decreasing from the center

C’est la notions d’externalités résultats des interactions entre firmes et/ou ménages.

Cette dimension est assez fondamentale, en témoigne la visualisation de l’épidémie de choléra à New York en 1832.

En effet, la prise en compte de la dimension spatiale permet de voir où est née l’épidémie de choléra, et de montrer l’importance des points d’eau, comme l’avait noté John Snow lors de l’épidémie de Londres en 1854 (raconté dans un article paru la semaine passée).

Les exposés de vendredi après midi étant autour des thèmes “nouvelles méthodes, nouveaux outils” et “nouveaux usages“, ça sera l’occasion de présenter des travaux récents avec Ewen sur la création de cartes évoluant dans le temps.

A titre indicatif, nous allons travailler cet été sur la création de telles cartes pour étudier le déplacement de descendants d’un couple, qui s’est marié en 1800 dans un village du Finistère, du bocage normand, ou du sud de la Bourgogne. A suivre donc….

Automatic Detection of the Language of a Tweet

Two days ago, in my post to extract automatically my own tweets, and to generate some html list, I mentioned that it would be great if there were a function that could be used to distinguish tweets in English, and tweets in French (usually, I tweet in one of those two languages). And one more time, @3wen came to rescue me! In my previous post, I spent some time to extract mentions to other twitter accounts, and urls. This time, I have to remove those two informations, that are not relevant to get the language of a tweet.

Continue reading Automatic Detection of the Language of a Tweet

Names in the U.S., from James Smith to Jose Rodriguez

Two weeks ago, @mona published an interesting post on her blog, about a difficult question, What’s The Most Common Name In America? There were stats about first names, in the U.S., and last names, too. Those informations are – somehow – easy to get. But usually, it is more complicated to get the first and the last name together. For confidentiality issues ! Datasets – the ones I deal with – are supposed to be anonymized, so I never see the first and the last names.  In a previous post, a few years ago, I did mention the so-called Social Security Death Master File. In that file, we have Social Security numbers, with the date of birth, the date of death as well as the first and the last name. So I did use those files to get stats about the first and the last names of American citizens. Of course, it is very restrictive. I have only U.S. citizens that have a Social Security number (which is not compulsary in the U.S. as far as I understood) and who passed away (as mentioned in the name of the dataset: the death master file). Another great thing about that dataset is that I have the date of birth, so I can look at some cohort effect (see opendata.stackexchange for an interesting discussion on that dataset).

Continue reading Names in the U.S., from James Smith to Jose Rodriguez

Shapefiles from Isodensity Curves

Recently, with @3wen, we wanted to play with isodensity curves. The problem is that it is difficult to get – numerically – the equation of the contour (even if we can easily plot it). Consider the following surface (just for fun, in order to illustrate the idea)

> f=function(x,y) x*y+(1-x)*(1-y)
> u=v=seq(0,1,length=21)
> v=seq(0,1,length=11)
> f=outer(u,v,f)
> persp(u,v,f,theta=angle,phi=10,box=TRUE,
+ shade=TRUE,ticktype="detailed",xlab="",
+ ylab="",zlab="",col="yellow")

For instance, assume that we want to locate areas where the density exceed 0.7 (here in the lower left corner, SW, and the upper right corner, NE)

> image(u,v,f)
> contour(u,v,f,add=TRUE,levels=.7)

Continue reading Shapefiles from Isodensity Curves

Extracting datasets from excel files in a zipped folder

The title of the post is a bit long, but that’s the problem I was facing this morning: importing dataset from files, online. I mean, it was not a “problem” (since I can always download, and extract manually the files), more a challenge (I should be able to do it in R, directly). The files are located on ressources-actuarielles.net, in a zip file. Those are mortality tables used in French speaking African countries, and I guess that one problem came from special characters, such as “é” or “è”… When you open the zip file, you see a folder

and in that folder, several files that I would like to import

Continue reading Extracting datasets from excel files in a zipped folder

Kernel Density Estimation with Ripley’s Circumferential Correction

The revised version of the paper Kernel Density Estimation with Ripley’s Circumferential Correction is now online, on hal.archives-ouvertes.fr/.

In this paper, we investigate (and extend) Ripley’s circumference method to correct bias of density estimation of edges (or frontiers) of regions. The idea of the method was theoretical and difficult to implement. We provide a simple technique — based of properties of Gaussian kernels — to efficiently compute weights to correct border bias on frontiers of the region of interest, with an automatic selection of an optimal radius for the method. We illustrate the use of that technique to visualize hot spots of car accidents and campsite locations, as well as location of bike thefts.

There are new applications, and new graphs, too

Most of the codes can be found on https://github.com/ripleyCorr/Kernel_density_ripley (as well as datasets).

We are winter

Non, pas de billet pour annoncer la reprise de Game of Thrones dans les jours à venir. Mais un billet rapide pour rappeler qu’au Canada, l’hiver, c’est un sujet important, et sérieux.

C’est ce thème qu’avait repris le Comité Olympique Canadien avant et pendant les Jeux.

Et lors de mon excursion au Mexique la semaine passée (où nous avons retrouvé des collègues chiliens d’ailleurs), je me suis rendu compte que les Jeux Olympiques d’hiver n’intéressaient que peu de monde, en dehors des Canadiens, et des nations qui subissent l’hiver, comme les Russes, ou les Norvégiens (entre autres). Un autre point comment de toutes ces nations est d’avoir un intérêt très fort pour un endroit important d’un point de vue géopolitique, le pôle nord. Et on m’a demandé d’intervenir dans une conférence de Thierry Garcin, et Philippe Reka. J’ai toujours un plaisir fou à discuter avec Philippe, donc c’est avec plaisir que je viendrais, mais j’avoue (une nouvelle fois) avoir du mal à trouver ma place (oui, je suis un peu le sparadrap du capitaine Haddock)

Il y a tout juste un an, Philippe publiait dans le Monde Diplo un article inspirant sur la Cartographie Radicale,

Depuis des décennies, la cartographie traditionnelle revendique le statut de science exacte s’appuyant sur des données fiables. Elle se targue de fournir une image neutre et fidèle de la réalité. Mais une telle approche fait l’impasse sur l’utilisation politique et sociale de la carte, et sur son rôle tant de propagande que de contestation. Depuis le début des années 2000, et dans le désordre, émerge une pratique cartographique qui se dit « radicale » (on parle aussi de « cartographie critique » ou de « contre-cartographie »), riche combinaison revendiquée d’art, de sciences, de géographie, de politique et de militantisme social. Les cartographes orthodoxes ont vu naître avec beaucoup de méfiance les premiers projets exploratoires qui bousculaient, voire cassaient, les codes. « On peut “crobarder” des schémas réducteurs donnant l’illusion de la synthèse, et en faire des figures de style prétentieuses, écrivait un cartographe professionnel en écho à l’une des toutes premières expositions de cette mouvance. La carte géographique, ce fantastique outil de savoir et d’appréhension de l’espace et du monde, se voit ravalée au rang de simple outil de prestidigitation. »

C’est peu dire que le symposium « Art and Cartography — Cartography and Art », organisé par l’Académie des beaux-arts et l’Université technique de Vienne (Autriche) en février 2008, fut le théâtre d’intenses discussions, où s’exprimèrent des positions inconciliables. La cartographie engagée a trouvé la voie d’une forme de contestation au service de la justice sociale, dénonçant des pratiques politiques et économiques douteuses. Elle est un exercice libre de déconstruction de l’espace et des phénomènes sociaux, pour lequel les protagonistes se permettent de pervertir les conventions les plus classiques. Les premiers à s’en emparer furent des artistes, mais aussi des architectes, des urbanistes et des militants, rejoints plus tard par des géographes. Les projets, individuels ou collectifs, se sont surtout multipliés en Argentine, aux Etats-Unis, en Suisse, en Allemagne et en France, couvrant un champ thématique assez large (finance, surveillance et sécurité, marketing, environnement, organisation des territoires, pour ne donner que quelques exemples), et profitant largement de la puissance des outils de la cartographie participative et des réseaux sociaux.

Déchiffrer pour dénoncer, mettre en images des processus peu visibles qui concourent à confisquer l’espace public (voire des biens publics), à compromettre les libertés individuelles, à détourner des lois, tels sont les principaux objectifs des promoteurs de cette cartographie. Dans cet esprit, informer n’est qu’un début. L’étape suivante de la démarche reste l’action pour le changement sur le terrain. Largement informelle, cette initiative se reconnaît quelques affinités avec les mouvements de résistance tels que celui des « indignés » ou Occupy. Ainsi, le collectif d’artistes argentins Grupo de Arte Callejero ne s’est pas contenté de dénoncer les bourreaux de l’ancienne dictature en cartographiant les lieux de leur paisible retraite : il a aussi organisé des itinéraires « touristiques » qui passaient par leurs domiciles et mené campagne sur des panneaux publicitaires… Le projet « Duty Free Shop », en ce qu’il se propose de déconstruire, rendre lisible et donc fragiliser un processus d’aliénation des citoyens dans un espace public, s’inscrit de plain-pied dans ce mouvement pétillant.

En tant qu’amateur (plus qu’expert) de visualisation de données, je rejoins complètement Philippe dans sa lecture critique (je ne vois les cartes que comme un support graphique parmi beaucoup d’autres). Philippe a ainsi proposé une très jolie visualisation de l’arctique.

Et comme souvent, Philippe a raison. Pour comprendre l’arctique, il faut avoir une visualisation claire. Car comme beaucoup de monde, j’ai été éduqué avec des cartes qui ressemblaient à ça

Autrement dit, pour moi, l’arctique et le pôle nord ont longtemps été des concepts flous, un peu comme l’horizon… c’est au nord, tout là haut. Par contre, j’ai bien conscience que l’antarctique existe, par exemple, mais les cartes usuelles ne nous permettent de prendre conscience de sa superficie (un peu moins de deux fois la taille de l’Australie). D’ailleurs, si on va sur https://www.google.ca/maps/place/North+Pole/ on peut voir (en butant contre la partie supérieure de la carte) qu’on ne voit rien, si on cherche à visualiser ce pôle nord (même s’il est sous la glace)

Comme le notait ma fille, le pôle nord est une enjeu stratégique important car c’est là que se trouve le village du père Noël (c’est manifestement la croyance répandue en Amérique). Et si le père Noël n’existe pas (c’est ce que pense ma fille), peut être que le pôle nord n’existe pas non plus ? Le premier défi est donc d’avoir une visualisation aussi pertinente que possible du pôle. Je vais donc revenir un peu sur les problématique de projection. Car les Soviétiques ne se trompaient pas lorsqu’ils représentèrent l’arctique avec un point de vue qui donnait l’impression que l’arctique était russe,

Pour appréhender les projecteurs, deux ouvrages fabuleux expliquent cela très clairement. Flattening the Earth – Two Thousand Years of Map Projections de John Parr Snyder, mais aussi (et surtout) How to Lie with Maps de Mark Monmonier. Mark Monmonier a également publié un livre passionnant sur le sujet (en particulier le chapitre 6), intitulé Air Apparent – How Meteorologists Learned to Map, Predict, and Dramatize Weather

Pour visualiser le pôle nord, on peut utiliser quelques fonctions simples (basées sur un travail en cours avec @3wen), pour visualiser différentes projections. On va ici mettre l’accent sur trois pays importantes quand on parle du pôle nord, le Canada, la Russie, et le Groenland (qui est techniquement un pays consultatif), et sur trois deux régions, une située à moins de 1500km du pôle nord, et l’autre à 3000km du pôle (soit un peu en dessous du cercle arctique)

Si on veut changer de point de vue, on peut visualiser l’animation suivante

D’ailleurs, si on compare la proportion du disque occupé par chacun des pays (j’ai pris ici chacune des huit pays de l’arctique), on obtient, avec en abscisse la distance au pôle, et en ordonnée, la proportion du disque occupée par chacun des pays, avec le Canada (rouge), la Russie (bleu) et le Groenland (vert)

Techniquement, le pôle pose de réels soucis car la plupart des données sont représentées par leur latitude, et leur longitude. Les pôles sont les seuls points – dans cette représentation – à avoir une infinité de coordonnées !

File:Latitude and Longitude of the Earth.svg

La latitude est de 90°, mais toutes les longitudes sont valides. On a ainsi des soucis y compris sur les cartes construites à partir de données satellites, qui pourtant pourraient balayer uniformément la sphère.

Le pôle nord peut ainsi compter beaucoup de discontinuités, uniquement à cause de la représentation des données par la latitude et la longitude : en quelques kilomètres, on passe d’une zone rouge à une zone blanche, ce qui pose des soucis de lissage de la carte (les données n’étant manifestement pas considérées ici comme circulaire pour le lissage).

Maintenant, plus qu’un point sur la carte, on peut voir le pôle nord comme une région, et non pas comme un point sur la sphère, défini très arbitrairement. Bon, j’exagère un peu ici, car le pôle nord est défini par l’axe de rotation de la terre,

qui va induire une espère de symétrie, par exemple sur le climat. Mais le second problème, selon moi, qui va alors apparaître est que le pôle nord vu comme un territoire (un territoire de glace) est alors un territoire qui va changer avec les saisons.

Par exemple, en septembre dernier, la glace occupait la région suivante, ce qui montre que peu de pays peuvent effectivement rejoindre le pôle nord ‘à pieds’ quel que soit le moment de l’année.

On continue de travailler avec @3wen sur le sujet, mais le sujet et vaste, ambitieux… et complexe ! à suivre donc…

Visualizing densities of spatial processes

We recently uploaded on http://hal.archives-ouvertes.fr/hal-00725090 a revised version of our work, with Ewen Gallic (a.k.a. @3wen) on Visualizing spatial processes using Ripley’s correction: an application to bodily-injury car accident location

In this paper, we investigate (and extend) Ripley’s circumference method to correct bias of density estimation of edges (or frontiers) of regions. The idea of the method was theoretical and di#cult to implement. We provide a simple technique – based of properties of Gaussian kernels – to compute e#efficiently weights to correct border bias on frontiers of the region of interest, with an automatic selection of an optimal radius for the method. An illustration on location of bodily-injury car accident (and hot spots) in the western part of France is discussed, where a lot of accident occur close to large cities, next to the sea.

Sketches of the R code can be found in the paper, to produce maps, an to describe the impact of our boundary correction. For instance, in Finistère, the distribution of car accident is the following (with a standard kernel on the left, and with correction on the right), with 186 claims (involving bodily injury)

and in Morbihan with 180 claims, observed in a specific year (2008 as far as I remember),

The code is the same as the one mentioned last year, except perhaps plotting functions. First, one needs to defi
ne a color scale and associated breaks

breaks <- seq( min( result $ZNA , na.rm = TRUE ) * 0.95 , max ( result$ZNA , na.rm = TRUE ) * 1.05 , length = 21)
col <- rev( heat . colors (20) )

to
finally plot the estimation

image . plot ( result $X, result $Y, result $ZNA , xlim = range (pol[,
1]) , ylim = range (pol[, 2]) , breaks = breaks , col = col ,
xlab = "", ylab = "", xaxt = "n", yaxt = "n", bty = "n",
zlim = range ( breaks ), horizontal = TRUE )

It is possible to add a contour, the observations, and the border of the polygon

contour ( result $X, result $Y, result $ZNA , add = TRUE , col = "grey ")
points (X[, 1], X[, 2], pch = 19, cex = 0.2 , col = " dodger blue")
polygon (pol , lwd = 2)

Now, if one wants to improve the aesthetics of the map, by adding a Google Maps base map, the
first thing to do – after loading ggmap package – is to get the base map

theMap <- get_map( location =c( left =min (pol [ ,1]) , bottom =min (pol[ ,2]) , right =max (pol [ ,1]) , 
top =max (pol [ ,2])), source =" google ", messaging =F, color ="bw")

Of course, data need to be put in the right format

getMelt <- function ( smoothed ){
res <- melt ( smoothed $ZNA)
res [ ,1] <- smoothed $X[res [ ,1]]
res [ ,2] <- smoothed $Y[res [ ,2]]
names (res) <- list ("X","Y","ZNA")
return (res )
}
smCont <- getMelt ( result )

Breaks and labels should be prepared

theLabels <- round (breaks ,2)
indLabels <- floor (seq (1, length ( theLabels ),length .out =5)) 
indLabels [ length ( indLabels )] <- length ( theLabels ) 
theLabels <- as. character ( theLabels [ indLabels ])
theLabels [ theLabels =="0"] <- " 0.00 "

Now, the map can be built

P <- ggmap ( theMap )
P <- P + geom _ point (aes(x=X, y=Y, col=ZNA), alpha =.3 , data =
smCont [!is.na( smCont $ZNA ) ,], na.rm= TRUE )

It is possible to add a contour

P <- P + geom _ contour ( data = smCont [!is.na( smCont $ZNA) ,] ,aes(x=
X, y=Y, z=ZNA ), alpha =0.5 , colour =" white ")

and colors need to be updated

P <- P + scale _ colour _ gradient ( name ="", low=" yellow ", high ="
red", breaks = breaks [ indLabels ], limits = range ( breaks ),
labels = theLabels )

To remove the axis legends and labels, the theme should be updated

P <- P + theme ( panel . grid . minor = element _ line ( colour =NA), panel
. grid . minor = element _ line ( colour =NA), panel . background =
element _ rect ( fill =NA , colour =NA), axis . text .x= element _ blank() ,
axis . text .y= element _ blank () , axis . ticks .x= element _ blank() ,
axis . ticks .y= element _ blank () , axis . title = element _ blank() , rect = element _ blank ())

The
final step, in order to draw the border of the polygon

polDF <- data . frame (pol)
colnames ( polDF ) <- list ("lon","lat")
(P <- P + geom _ polygon ( data =polDF , mapping =( aes(x=lon , y=lat)), colour =" black ", fill =NA))

Then, we’ve applied that methodology to estimate the road network density in those two regions, in order to understand if high intensity means that it is a dangerous area, or if it simply because there is a lot of traffic (more traffic, more accident),

We have been using the dataset obtained from the Geofabrik website which provides
Open-StreetMap data. Each observation is a section of a road, and contains a few points identifi
ed by their geographical coordinates that allow to draw lines. We have use those points to estimate a proxy of road intensity, with weight going from 10 (highways) to 1 (service roads).

splitroad <- function ( listroad , h = 0.0025) {
pts = NULL
weights <- types . weights [ match ( unique ( listroad $ type ), types .
weights $ type ), " weight "]
for (i in 1:( length ( listroad ) - 1)) {
d = diag (as. matrix ( dist ( listroad [[i]]))[, 2: nrow ( listroad
[[i ]]) ]))
}}
return (pts )
}

See Ewen’s blog for more details on the code, http://editerna.free.fr/blog/…. Note that Ewen did publish a poster of the paper (in French), for the http://r2013-lyon.sciencesconf.org/ conference, that will be organized in Lyon on June 27th-28th, see

All comments are welcome…

Border bias and weighted kernels

With Ewen (aka @3wen), not only we have been playing on Twitter this month, we have also been working on kernel estimation for densities of spatial processes. Actually, it is only a part of what he was working on, but that part on kernel estimation has been the opportunity to write a short paper, that can now be downloaded on hal.

The problem with kernels is that kernel density estimators suffer a strong bias on borders. And with geographic data, it is not uncommon to have observations very close to the border (frontier, or ocean). With standard kernels, some weight is allocated outside the area: the density does not sum to one. And we should not look for a global correction, but for a local one. So we should use weighted kernel estimators (see on hal for more details). The problem that weights can be difficult to derive, when the shape of the support is a strange polygon. The idea is to use a property of product Gaussian kernels (with identical bandwidth) i.e. with the interpretation of having noisy observation, we can use the property of circular isodensity curve. And this can be related to Ripley (1977) circumferential correction. And the good point is that, with R, it is extremely simple to get the area of the intersection of two polygons. But we need to upload some R packages first,

require(maps)
require(sp)
require(snow)
require(ellipse)
require(ks)
require(gpclib)
require(rgeos)
require(fields)

To be more clear, let us illustrate that technique on a nice example. For instance, consider some bodiliy injury car accidents in France, in 2008 (that I cannot upload but I can upload a random sample),

base_cara=read.table(
"http://freakonometrics.blog.free.fr/public/base_fin_morb.txt",
sep=";",header=TRUE)

The border of the support of our distribution of car accidents will be the contour of the Finistère departement, that can be found in standard packages

geoloc=read.csv(
"http://freakonometrics.free.fr/popfr19752010.csv",
header=TRUE,sep=",",comment.char="",check.names=FALSE,
colClasses=c(rep("character",5),rep("numeric",38)))
geoloc=geoloc[,c("dep","com","com_nom",
"long","lat","pop_2008")]
geoloc$id=paste(sprintf("%02s",geoloc$dep),
sprintf("%03s",geoloc$com),sep="")
geoloc=geoloc[,c("com_nom","long","lat","pop_2008")]
head(geoloc)
france=map('france',namesonly=TRUE,
plot=FALSE)
francemap=map('france', fill=TRUE, col="transparent",
plot=FALSE)
detpartement_bzh=france[which(france%in%
c("Finistere","Morbihan","Ille-et-Vilaine",
"Cotes-Darmor"))]
bretagne=map('france',regions=detpartement_bzh,
fill=TRUE, col="transparent", plot=FALSE,exact=TRUE)
finistere=cbind(bretagne$x[321:678],bretagne$y[321:678])
FINISTERE=map('france',regions="Finistere", fill=TRUE,
col="transparent", plot=FALSE,exact=TRUE)
monFINISTERE=cbind(FINISTERE$x[c(8:414)],FINISTERE$y[c(8:414)])

Now we need simple functions,

cercle=function(n=200,centre=c(0,0),rayon)
{theta=seq(0,2*pi,length=100)
m=cbind(cos(theta),sin(theta))*rayon
m[,1]=m[,1]+centre[1]
m[,2]=m[,2]+centre[2]
names(m)=c("x","y")
return(m)}
poids=function(x,h,POL)
{leCercle=cercle(centre=x,rayon=5/pi*h)
POLcercle=as(leCercle, "gpc.poly")
return(area.poly(intersect(POL,POLcercle))/
area.poly(POLcercle))}
lissage = function(U,polygone,optimal=TRUE,h=.1)
{n=nrow(U)
IND=which(is.na(U[,1])==FALSE)
U=U[IND,]
if(optimal==TRUE) {H=Hpi(U,binned=FALSE);
H=matrix(c(sqrt(H[1,1]*H[2,2]),0,0,
sqrt(H[1,1]*H[2,2])),2,2)}
if(optimal==FALSE){H= matrix(c(h,0,0,h),2,2)

before defining our weights.

poidsU=function(i,U,h,POL)
{x=U[i,]
poids(x,h,POL)}
OMEGA=parLapply(cl,1:n,poidsU,U=U,h=sqrt(H[1,1]),
POL=as(polygone, "gpc.poly"))
OMEGA=do.call("c",OMEGA)
stopCluster(cl)
}else
{OMEGA=lapply(1:n,poidsU,U=U,h=sqrt(H[1,1]),
POL=as(polygone, "gpc.poly"))
OMEGA=do.call("c",OMEGA)}

Note that it is possible to parallelize if there are a lot of observations,

if(n>=500)
{cl <- makeCluster(4,type="SOCK")
worker.init <- function(packages)
{for(p in packages){library(p, character.only=T)}
NULL}
clusterCall(cl, worker.init, c("gpclib","sp"))
clusterExport(cl,c("cercle","poids"))

Then, we can use standard bivariate kernel smoothing functions, but with the weights we just calculated, using a simple technique that can be related to one suggested in Ripley (1977),

fhat=kde(U,H,w=1/OMEGA,xmin=c(min(polygone[,1]),
min(polygone[,2])),xmax=c(max(polygone[,1]),
max(polygone[,2])))
fhat$estimate=fhat$estimate*sum(1/OMEGA)/n
vx=unlist(fhat$eval.points[1])
vy=unlist(fhat$eval.points[2])
VX = cbind(rep(vx,each=length(vy)))
VY = cbind(rep(vy,length(vx)))
VXY=cbind(VX,VY)
Ind=matrix(point.in.polygon(VX,VY, polygone[,1],
polygone[,2]),length(vy),length(vx))
f0=fhat
f0$estimate[t(Ind)==0]=NA
return(list(
X=fhat$eval.points[[1]],
Y=fhat$eval.points[[2]],
Z=fhat$estimate,
ZNA=f0$estimate,
H=fhat$H,
W=fhat$W))}
lissage_without_c = function(U,polygone,optimal=TRUE,h=.1)
{n=nrow(U)
IND=which(is.na(U[,1])==FALSE)
U=U[IND,]
if(optimal==TRUE) {H=Hpi(U,binned=FALSE);
H=matrix(c(sqrt(H[1,1]*H[2,2]),0,0,sqrt(H[1,1]*H[2,2])),2,2)}
if(optimal==FALSE){H= matrix(c(h,0,0,h),2,2)}
fhat=kde(U,H,xmin=c(min(polygone[,1]),
min(polygone[,2])),xmax=c(max(polygone[,1]),
max(polygone[,2])))
vx=unlist(fhat$eval.points[1])
vy=unlist(fhat$eval.points[2])
VX = cbind(rep(vx,each=length(vy)))
VY = cbind(rep(vy,length(vx)))
VXY=cbind(VX,VY)
Ind=matrix(point.in.polygon(VX,VY, polygone[,1],
polygone[,2]),length(vy),length(vx))
f0=fhat
f0$estimate[t(Ind)==0]=NA
return(list(
X=fhat$eval.points[[1]],
Y=fhat$eval.points[[2]],
Z=fhat$estimate,
ZNA=f0$estimate,
H=fhat$H,
W=fhat$W))}

So, now we can play with those functions,

base_cara_FINISTERE=base_cara[which(point.in.polygon(
base_cara$long,base_cara$lat,monFINISTERE[,1],
monFINISTERE[,2])==1),]
coord=cbind(as.numeric(base_cara_FINISTERE$long),
as.numeric(base_cara_FINISTERE$lat))
nrow(coord)
map(francemap)
lissage_FIN_withoutc=lissage_without_c(coord,
monFINISTERE,optimal=TRUE)
lissage_FIN=lissage(coord,monFINISTERE,
optimal=TRUE)
lesBreaks_sans_pop=range(c(
range(lissage_FIN_withoutc$Z),
range(lissage_FIN$Z)))
lesBreaks_sans_pop=seq(min(lesBreaks_sans_pop)*.95,
max(lesBreaks_sans_pop)*1.05,length=21)

plot_article=function(lissage,breaks,
polygone,coord){
par(mar=c(3,1,3,1))
image.plot(lissage$X,lissage$Y,(lissage$ZNA),
xlim=range(polygone[,1]),ylim=range(polygone[,2]),
breaks=breaks, col=rev(heat.colors(20)),xlab="",
ylab="",xaxt="n",yaxt="n",bty="n",zlim=range(breaks),
horizontal=TRUE)
contour(lissage$X,lissage$Y,lissage$ZNA,add=TRUE,
col="grey")
points(coord[,1],coord[,2],pch=19,cex=.1,
col="dodger blue")
polygon(polygone,lwd=2,)}

plot_article(lissage_FIN_withoutc,breaks=
lesBreaks_sans_pop,polygone=monFINISTERE,
coord=coord)

plot_article(lissage_FIN,breaks=
lesBreaks_sans_pop,polygone=monFINISTERE,
coord=coord)

If we look at the graphs, we have the following densities of car accident, with a standard kernel on the left, and our proposal on the right (with local weight adjustment when the estimation is done next to the border of the region of interest),

Similarly, in Morbihan,

With those modified kernels, hot spots appear much more clearly. For more details, the paper is online on hal.

Qui se ressemble se suit (sur Twitter au moins)

Un nouveau billet pour reprendre une analyse marrante faite par @3wen  (Ewen Gallic, qui travaille à Montréal alors que je profite de la Suisse). Suite à l’analyse amusante des trolls de Twitter, j’avais lancé l’idée qu’il pourrait être amusant de regarder parmi les députés français (que j’avais un peu suivis l’autre jour), qui tweete avec qui. Ewen a eu la bonne idée de regarder sur  http://lelab.europe1.fr/ ce qui lui a permis de récupérer la liste des comptes Twitter des députés, en France. L’idée est simple: parmi les députés français, on regarde qui suit qui. Quelqu’un de très suivi sera au centre du nuage, alors que quelqu’un qui se contente de suivre sera sur le bord (intensément connecté au reste du nuage). Pour les personnes peu familières, Twitter n’est pas Facebook: on n’a pas des “amis“, il y a des gens que l’on suit parce qu’ils racontent des choses qui peuvent nous intéresser.

En utilisant gephi Ewen a ensuite pu visualiser le graphique des interconnexions entre les députés.

Sans grande surprise, les députés de gauche suivent essentiellement les députés de gauche, et inversement. Quelques gros comptes font la passerelle entre les deux groupes parlementaires. En fait, si on regarde en détail (voire sur le fichier image complet), on peut observer un peu mieux ce qui se passe.

Bon, la grande difficulté est de lire ces interactions correctement. En particulier, on ne peut pas conclure (à la vue seule du graphique) que Cécile Duflot est de gauche ! Ce que cela nous dit est que ce qu’elle raconte intéresse les gens de gauche (ou en tous cas les députés du Parti Socialiste), et pas du tout les gens de droite (les députés de l’UMP) ! On notera aussi que le centre n’existe pas. Sur Twitter en tous cas. Et si on regarde tout en bas, à droite, on retrouve le Front National, et on a la confirmation que ce que raconte le Front National n’intéresse personne !

Do you still have time to sleep ?

Last week, @3wen (Ewen) helped me to write nice R functions to extract tweets in R and build datasets containing a lot of information. I’ve tried a couple of time on my own. Once on tweet contents, but it was not convincing and once on the activity on Twitter following an event (e.g. the death of someone famous). I have to admit that I am not a big fan of databases that can be generated using standard function to study tweets. For instance, we can only extract tweets, notre-tweets (which is also an important indicator of tweet-activity). @3wen suggested to use

require("RJSONIO")

The first step is to extract some information from a tweet, and store it in a dataset (details can be found on https://dev.twitter.com/)

obtenir_ligne <- function(unTweet){
date_courante=unTweet$created_at
id_courant=unTweet$id_str
text=unTweet$text
nb_followers=unTweet$user$followers_count
nb_amis=unTweet$user$friends_count
utc_offset=unTweet$user$utc_offset
listeMentions=unTweet$entities$user_mentions
return(c(list(c(id_courant,date_courante,text,
nb_followers,nb_amis,utc_offset)),
list(do.call("rbind",lapply(listeMentions,
function(x,id_courant) c(id_courant,
x$screen_name),unTweet$id_str)))))
}

Now that we  have the code to extract information from one tweet, let us find several tweets, from one user, say my account,

nom="Freakonometrics"

The (small) problem here, is that we have a limitation: we can only get 100 tweets per call of the function

n=100
tweets_courants=scan(paste(
"http://api.twitter.com/1/statuses/user_timeline.json?
include_entities=true&include_rts=true&screen_name=
",nom,"&count=",n,sep=""),what = "character",
encoding="latin1")
tweets_courants=paste(tweets_courants[
1:length(tweets_courants)],collapse=" ")
tweets_courants=fromJSON(tweets_courants,
method = "C")

Then, we use our function to build a database with 100 lines,

extracTweets <- lapply(tweets_courants,
obtenir_ligne)
mentions=do.call("rbind",lapply(extracTweets,
function(x) x[[2]]))
colnames(mentions)=list("id","screen_name")
res=t(sapply(extracTweets,function(x) x[[1]]))
colnames(res) <- list("id","date","text",
"nb_followers","nb_amis","utc_offset")

The idea then is simply to use a loop, based on the latest id observed

dernier_id=tweets_courants[[length(
tweets_courants)]]$id_str

So, here we go,

compteurLimite=100

while(compteurLimite<4100){
tweets_courants=scan(paste(
"http://api.twitter.com/1/statuses/user_timeline.json?
include_entities=true&include_rts=true&screen_name=
",nom,"&count=",n,"&max_id=",dernier_id,sep=""),
what = "character", encoding="latin1")
tweets_courants=paste(tweets_courants[
1:length(tweets_courants)],collapse=" ")
tweets_courants=fromJSON(tweets_courants,
method = "C")

extracTweets <- lapply(tweets_courants[
2:length(tweets_courants)],obtenir_ligne)
mentions=rbind(mentions,do.call("rbind",
lapply(extracTweets,function(x) x[[2]])))
res=rbind(res,t(sapply(extracTweets,function(x) x[[1]])))
t(sapply(extracTweets,function(x) x[[1]]))
dernier_id=tweets_courants[[length(
tweets_courants)]]$id_str
compteurLimite=compteurLimite+100
}

resFreakonometrics=res=
data.frame(res,stringsAsFactors=FALSE)

All the information about my own tweets (and re-tweets) are stored in a nice dataset. Actually, we have even more, since we have extracted also names of people mentioned in tweets,

mentionsFreakonometrics=
data.frame(mentions)

We can look at people I mention in my tweets

gazouillis=sapply(split(mentionsFreakonometrics,
mentions$screen_name),nrow)
gazouillis=gazouillis[order(gazouillis,
decreasing=TRUE)]

plot(gazouillis)
plot(gazouillis,log="xy")
> gazouillis[1:20]
tomroud freakonometrics       adelaigue       dmonniaux
155              84              77              56
J_P_Boucher         embruns      SkyZeLimit        coulmont
42              39              35              31
Fabrice_BM            3wen          obouba          msotod
31              30              29              27
StatFr     nholzschuch        renaudjf        squintar
26              25              23              23
Vicnent        pareto35        romainqc        valatini
23              22              22              22

If we plot those frequencies, we can clearly observe a standard Pareto distribution,

Now, let us spend some time with dates and time of tweets (it was the initial goal of this post)… One more time, there is a (small) technical problem that we have to deal with: language. We need a function to convert date in English (on Twitter) to dates in French (since I have a French version of R),

changer_date_anglais <- function(date_courante){
mois <- c("Jan","Fév", "Mar", "Avr", "Mai",
"Jui", "Jul", "Aoû", "Sep", "Oct", "Nov", "Déc")
months <- c("Jan", "Feb", "Mar", "Apr", "May",
"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
jours <- c("Lun","Mar","Mer","Jeu",
"Ven","Sam","Dim")
days <- c("Mon","Tue","Wed","Thu",
"Fri","Sat","Sun")
leJour <- substr(date_courante,1,3)
leMois <- substr(date_courante,5,7)
return(paste(jours[match(leJour,days)]," ",
mois[match(leMois,months)],substr(
date_courante,8,nchar(date_courante)),sep=""))
}

So now, it is possible to plot the times where I am online, tweeting,

DATE=Vectorize(changer_date_anglais)(res$date)
DATE=sapply(resSkyZeLimit$date,
changer_date_anglais,simplify=TRUE)

DATE2=strptime(as.character(DATE),
"%a %b %d %H:%M:%S %z %Y")
lt= as.POSIXlt(DATE2, origin="1970-01-01")
heure=lt$hour+lt$min/60
plot(DATE2,heure)

On this graph, we can see that I am clearly not online almost 6 hours a day (or at least not on Twitter). It is possible to visualize more precisely the period of the day where I might be on Twitter,

hist(heure,breaks=0:24,col="light green",proba=TRUE)
X=c(heure-24,heure,heure+24)
d=density(X,n = 512, from=0, to=24,bw=1)
lines(d$x,d$y*3,lwd=3,col="red")

or, if we want to illustrate with some kind of heat plot,

Note that we did it for my Twitter account, but we can also run the code on (almost) anyone on Twitter. Consider e.g. @adelaigue. Since Alexandre is tweeting in France, we have to play with time-zones,

res=extractR("adelaigue")
DATE=Vectorize(changer_date_anglais)(res$date)
DATE2=strptime(as.character(DATE),
"%a %b %d %H:%M:%S %z %Y",tz = "GMT")+2*60*60

or I can also look at @skythelimit who’s usually twitting from Singapore (I am in Montréal). I can seen clearly when we might have overlaps,

res=extractR("skythelimit")

Nice isn’t it. But it is possible to do much better… for instance, for those who do not ask specifically not to be Geo-located, we can see where they do tweet during the day, and during the night… I am quite sure a dozen posts with those functions can be written…

French dataset: population and GPS coordinates

A short post today based on recent work by @3wen (Ewen Gallic, graduate Student in Rennes, spending a year in Montreal). Since we were working on a detailed French dataset (per commune), we needed a dataset containing a list allcommunes, with population and location. GPS coordinates were extracted from Google, using the following php file, inspired by http://www.andrew-kirkpatrick.com/ on Google geocoding api with php webpage. Population was interpolated from INSEE’s datasets, i.e. http://www.insee.fr/ (since data are over a 35 year period, from 1975 to 2010, changes have been taken into account as carefully are possible – e.g. merges and splits of cities – based on thatdescription). A spline model has been used for all cities (with three degrees of freedom, and null and negative interpolation became one, since we’ll be using loglinear models afterwards). Names are from that dataset, still on INSEE’s website, http://www.insee.fr/.

A zipped file can be downloaded here popfr19752010.zip, but it is also possible to use the code below (it is a 24Mo dataset). Since it was hard to find such a dataset online (different files can be found, but we found none with population and location), we have decided to upload that dataset. Please let us know if there are problems with those data…

> base=read.csv(
+ "http://freakonometrics.free.fr/popfr19752010.csv",
+ header=TRUE)

Using that code, it is possible to locate all the communes in France (metropolitan), for instance

> library(maps)
> map("france")
> points(base$long,base$lat,cex=.1,col="red",pch=19)
> points(base$long,base$lat,cex=2*base$pop_2010/
+ max(base$pop_2010),col="blue",pch=19)

Several additional lines of code on that dataset (and also others) will be uploaded, soon.

Cette oeuvre est mise à disposition sous licence Paternité – Partage à l’Identique 3.0 non transposé. Pour voir une copie de cette licence, visitez http://creativecommons.org/. Date : 24 mai 2012, par Ewen GALLIC. Sources : INSEE, API Google Maps v3 et GeoHack (coordonnées GPS), propres calculs (estimation de population à partir des données INSEE).

  • reg : code region INSEE (character)
  • dep : code departement INSEE (character, corse 201 et 202 au lieu de 2A et 2B)
  • com : code commune INSEE (character)
  • article : article du nom de la commune (character)
  • com_nom : nom de la commune (character)
  • long : longitude (numeric)
  • lat : latitude (numeric)
  • pop_i : estimation de la population à la date i (ramenée à 1 si <=0), i=1975,…,2010 (numeric)