Dans le cours d’actuariat de l’assurance non-vie, on avait évoqué rapidement l’idée de faire un zonier. Autrement dit, on souhaite créer une variable polytomique (avec 4 ou 5 classes) représentant le critère spatial du risque. On va tenter de segmenter le territoire, en regroupant l’ensemble des territoires ‘proches’ (en terme de risque) étudiées en un nombre de classes prédéfini, permettant de tenir compte d’un (potentiel) critère spatial dans la tarification. Pour les aspects pratiques, comme on va manipuler des données spatiales, chargeons quelques packages pour commencer
library(maptools) library(rgeos) library(rgdal) library(ggplot2) library(plyr) library(maptools) library(cartography)
On va aussi avoir besoin de données spatialisées (disons en France métropolitaire). Pour faire simple, je vais tirer des codes insee de communes au hasard, un millier, et une variable Y qui va nous intéresser. Et je suppose qu’elle dépend de X que je n’observe pas, mais qui est liée à des caractéristiques spatiales (en gros la latitude, i.e. un positionnement nord/sud)
download.file("http://freakonometrics.free.fr/popfr19752010.csv","popfr.csv") base = read.csv("popfr.csv",header=TRUE) base$insee = base$dep*1000+base$com n=1000 set.seed(123) id=sample(1:nrow(base),size=n) simbase=data.frame(insee=base$insee[id]) X=(46-base$lat[id])/2 simbase$Y=X+rnorm(n)
Ici, on a la base suivante (mais on va supposer que je n’observe pas X – qui n’a servi qu’à simuler des données)
> head(simbase) insee Y 1 61499 -1.8573363 2 19181 -0.7059191 3 55307 -0.5923649 4 74030 0.6098773 5 30328 -0.4584795 6 10050 -1.2361240
Classiquement, Y peut être un résidu (normalisé) de régression, par exemple. Pour visualiser, il faut un fond de carte
download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds","FRA_adm0.rds") FR0=readRDS("FRA_adm0.rds")
Ensuite, on met les points (en fusionnant nos données avec la base insee donnant latitude et longitude des codes insee de communes)
plot(FR0) simbase = merge(simbase,base[,c("insee","long","lat")]) cols = rev(carto.pal(pal1 = "red.pal",n1 = 10,pal2="green.pal",n2=10)) bk = seq(-5,4.5,length=21) cuty = cut(simbase$Y,breaks=bk,labels=1:20) points(simbase$long,simbase$lat, col=cols[cuty],pch=19,cex=.5)
On va voir eux techniques pour faire un zonier. Le premier est de travailler par zone prédéfinie, comme le département,
simbase$dpt=trunc(simbase$insee/1000) A=aggregate(x = simbase$Y,by=list(simbase$dpt),mean) names(A)=c("dpt","y") download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm2.rds","FRA_adm2.rds") FR=readRDS("FRA_adm2.rds") donnees_carte=data.frame(FR@data) d=donnees_carte$CCA_2 d[d=="2A"]="201" d[d=="2B"]="202" donnees_carte$dpt=as.numeric(as.character(d)) donnees_carte=merge(donnees_carte,A,all.x=TRUE) donnees_carte=donnees_carte[order(donnees_carte$OBJECTID),] bk=seq(-2.75,2.75,length=21) donnees_carte$cuty=cut(donnees_carte$y,breaks=bk,labels=1:20) cols = rev(carto.pal(pal1 = "red.pal",n1 = 10,pal2="green.pal",n2=10)) plot(FR, col=cols[donnees_carte$cuty],xlim=c(-5.2,12))
autrement dit, on fait une carte choroplèthe
(il faut s’assurer que les couleurs sont mises au bon endroit). Ensuite, on peut définir deux zones,
bk=seq(-2.75,2.75,length=3) donnees_carte$cuty=cut(donnees_carte$y,breaks=bk,labels=1:2) plot(FR, col=cols[c(4,16)][donnees_carte$cuty],xlim=c(-5.2,12))
ou quatre
bk=seq(-2.75,2.75,length=5) donnees_carte$cuty=cut(donnees_carte$y,breaks=bk,labels=1:4) plot(FR, col=cols[c(3,8,12,17)][donnees_carte$cuty],xlim=c(-5.2,12))
On va créer uen variable à 4 modalités, à partir des départements.
Une autre solution consiste à utiliser les données spatiales, obtenues par fusion avec la base insee
simbase=merge(simbase,base[,c("insee","long","lat")])
On va ensuite se donner une grille, en France. C’est un peu pénible. Il nous faut le polynôme de contour de la France
P1=FR0@polygons[[1]]@Polygons[[355]]@coords P2=FR0@polygons[[1]]@Polygons[[27]]@coords plot(FR0,border=NA) polygon(P1) polygon(P2)
(ici on a juste la France métropolitaine, et la Corse – deux polygônes donc)
et on part d’une grille sur un rectangle
grille<-expand.grid(seq(min(simbase$long),max(simbase$long),length=101),seq(min(simbase$lat),max(simbase$lat),length=101)) paslong=(max(simbase$long)-min(simbase$long))/100 paslat=(max(simbase$lat)-min(simbase$lat))/100
On retient alors juste les points qui sont dans les polygônes
f=function(i){ (point.in.polygon (grille[i, 1]+paslong/2 , grille[i, 2]+paslat/2 , P1[,1],P1[,2])>0)+(point.in.polygon (grille[i, 1]+paslong/2 , grille[i, 2]+paslat/2 , P2[,1],P2[,2])>0) } indic=unlist(lapply(1:nrow(grille),f)) grille=grille[which(indic==1),] points(grille[,1]+paslong/2,grille[,2]+paslat/2,cex=.4,pch=19,col="blue")
voilà le résultat
Maintenant, on peut utiliser du krigeage mais j’ai plutôt voulu tenter des plus proches voisins. Pour chaque point de la grille, on prend la moyenne des plus proches voisins (à vol d’oiseau, i.e. avec une norme Euclidienne – sur la sphère)
library(geosphere) knn=function(i,k=20){ d=distHaversine(grille[i,1:2],simbase[,c("long","lat")], r=6378.137) r=rank(d) ind=which(r<=k) mean(simbase[ind,"Y"]) } grille$y=Vectorize(knn)(1:nrow(grille)) bk=seq(-2.75,2.75,length=21) grille$cuty=cut(grille$y,breaks=bk,labels=1:20) cols <- rev(carto.pal(pal1 = "red.pal",n1 = 10,pal2="green.pal",n2=10)) points(grille[,1]+paslong/2,grille[,2]+paslat/2,col=cols[grille$cuty],pch=19)
Ici, ca donne la carte suivante
mais là encore, on peut retenir deux niveaux
bk=seq(-2.75,2.75,length=3) grille$cuty=cut(grille$y,breaks=bk,labels=1:2) plot(FR0,border=NA) polygon(P1) polygon(P2) points(grille[,1]+paslong/2,grille[,2]+paslat/2,col=cols[c(4,16)][grille$cuty],pch=19)
ce qui donne les zones suivantes
bk=seq(-2.75,2.75,length=5) grille$cuty=cut(grille$y,breaks=bk,labels=1:4) plot(FR0,border=NA) polygon(P1) polygon(P2) points(grille[,1]+paslong/2,grille[,2]+paslat/2,col=cols[c(3,8,12,17)][grille$cuty],pch=19)
ou bien, si on retient quatre niveaux
A partir de là, pour créer une variable de zone, le plus simple est à un point donné de chercher le point de la grille le plus proche, et de lui donner la valeur de la couleur associée
pred=function(z){ d=distHaversine(z,grille[,1:2], r=6378.137) grille[which.min(d),"cuty"] }
Facile, non ?