{"word_count": 505, "character_count": 6264}
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 commencerlibrary(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 [latex]Y[/latex] qui va nous intéresser. Et je suppose qu'elle dépend de [latex]X[/latex] 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$comn=1000set.seed(123)id=sample(1:nrow(base),size=n)simbase=data.frame(insee=base$insee[id])X=(46-base$lat[id])/2simbase$Y=X+rnorm(n)Ici, on a la base suivante (mais on va supposer que je n'observe pas [latex]X[/latex] - qui n'a servi qu'à simuler des données)&#62; head(simbase)  insee          Y1 61499 -1.85733632 19181 -0.70591913 55307 -0.59236494 74030  0.60987735 30328 -0.45847956 10050 -1.2361240Classiquement, [latex]Y[/latex] peut être un résidu (normalisé) de régression, par exemple. Pour visualiser, il faut un fond de cartedownload.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_2d[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 quatrebk=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 inseesimbase=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 FranceP1=FR0@polygons[[1]]@Polygons[[355]]@coordsP2=FR0@polygons[[1]]@Polygons[[27]]@coordsplot(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 rectanglegrille&#60;-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))/100paslat=(max(simbase$lat)-min(simbase$lat))/100On retient alors juste les points qui sont dans les polygônesf=function(i){ (point.in.polygon (grille[i, 1]+paslong/2 , grille[i, 2]+paslat/2 , P1[,1],P1[,2])&#62;0)+(point.in.polygon (grille[i, 1]+paslong/2 , grille[i, 2]+paslat/2 , P2[,1],P2[,2])&#62;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ésultatMaintenant, 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&#60;=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 &#60;- 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 suivantemais là encore, on peut retenir deux niveauxbk=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 suivantesbk=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 niveauxA 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éepred=function(z){  d=distHaversine(z,grille[,1:2], r=6378.137)  grille[which.min(d),"cuty"]}Facile, non ?