Analyse Spatiale

Jean Gaudart
Aix-Marseille Université
UMR 912, SESSTIM (AMU, INSERM, IRD)

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 1

 plan
1. Introduction
2. Classes spatiales
3. Visualisation
4. Interpolation et Géostatistique
5. Analyses de co-facteurs
6. Autocorrélation

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 2

 I. Introduction
Cartographie
⇒Informations sur un phénomène géographique
⇒Mais pas suffisant pour extraire ses caractéristiques et comprendre le processus
sous-jacent
Inférence statistique indispensable:
⇒« Pattern » spatial de l’incidence d’une maladie => existe-t-il une agrégation
(cluster)?
⇒En cas de clusters, sont-ils expliqués par des co-facteurs (âge, pauvreté,
pollution, environnement…)?
⇒En analysant des échantillons de sol, peut-on déterminer quelle zone est
polluée?
⇒En analysant des prélèvements d’air (pour mesurer sa qualité), comment
estimer le nombre de personnes exposées à de haut niveau de particules fines? Et
où vivent-elles?
⇒Est-ce que les gouvernements ont tendance à appliquer des politiques
semblables à leurs voisins ou bien sont-ils indépendants?
⇒…

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 3

Utiliser R pour l’analyse spatiale:

⇒Nombre de packages +++ adaptés aux différentes hypothèses sur l’organisation

spatiale des données
⇒Package « sp » (Pebesma & Bivand, 2005)
⇒Classes et méthodes spécifiques
⇒Interfaces avec des Systèmes d’Informations Géographiques (SIG)

http://cran.r-project.org/web/views/Spatial.html

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 4

Système d’informations Géographique (SIG/GIS):
⇒Ensemble d’outils pour collecter, conserver, extraire à volonté, transformer, croiser,
manipuler, analyser des données spatialisées issues du monde réel, afin de répondre
un ensemble de questions
⇒Base de données particulière: géographie
⇒Ex: ArcGIS, QGIS, GRASS…
⇒Nombreuses données accessibles gratuitement sur internet

•Christman N. (2002) Exploring Geographic Information Systems, Wiley.

•Heywood I, Cornelius S, Carver S. (2006) An introduction to Geographic Information
Systems. Pearson Education.
•Longley PA, Goodchild MF, Maguire DJ, Rhind DW. (2005) Geographic Information
Systems and Science. Wiley.

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 5

Types de données spatiales
Données spatiales => références spatiales: Coordonnées dans un système de référence.
Exemple: localisation des volcans dans le monde. Coordonnées:
⇒longitude degrés décimaux / méridien 0 (Greenwich)
⇒Latitude degrés décimaux / 0 équateur

Système: World Geodetic System WGS84

Cartographier les volcans actifs 1980-2000:
⇒Points (long,lat)
⇒Carte aplatie => problème de projection:
transformer un système sphérique en un système plan
Par ex. projection de Mollweide

Données supplémentaires associées, non spatiales : Attributs

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 6

•Points: localisation ponctuelle, GPS, adresse géocodées Vectoriel:
représentation la plus
•Lignes: ensemble de points ordonnés, connectés par des précise et fiable possible
segments de droites

•Polygones: zones définies par des lignes fermées (+- trous)

•Grilles: ensemble de points ou de cellules rectangulaires Raster:

organisés en un treillis régulier Représentation de surfaces
continues par une tessalisation
régulière

Attributs d’une cellule: moyenne, valeur centrale, %, …

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 7

 II. Classe pour données spatiales -R
chooseCRANmirror(ind=6)
install.packages("sp")
library(sp)

II.1 Objet spatial

La classe spatial de base comprend 2 slots (parties):
-Un cadre délimité (bounding box bbox): matrice de coordonnées numérique, colonnes
(min, max), lignes (Est=x, Nord=Y)
-Système de projection (proj4string) : CRS-class coordinate reference system
getClass("Spatial")
Class "Spatial" [package "sp" ]
Slots:
Name: bbox proj4string
Class: matrix CRS
Known Subclasses:
Class "SpatialPoints", directly
Class "SpatialGrid", directly
Class "SpatialLines", directly
Class "SpatialPolygons", directly
… (c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 8
 Construction d’un objet Spatial simple:
m1<-matrix(c(0,0,1,1),nc=2,nr=2,dimnames=list(NULL,c("min","max")))
crs1<-CRS(projargs=as.character(NA))
SpOb1<-Spatial(bbox=m1, proj4string=crs1)
SpOb1 An object of class "Spatial"
Slot "bbox":
min max
[1,] 0 1
[2,] 0 1
Slot "proj4string":
CRS arguments: NA
m2<-matrix(c(350,85,360,90),nc=2,nr=2,dimnames=list(NULL,c("min","max")))
crs2<-CRS("+proj=longlat +datum=WGS84")
SpOb2<-Spatial(bbox=m2, proj4string=crs2)
SpOb2
An object of class "Spatial"
Slot "bbox":
min max
[1,] 350 360
[2,] 85 90
Slot "proj4string":
CRS arguments:
+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 9
II.2 SpatialPoints Sous-classe de la classe Spatial

Un point est précisément localisé sur le globe par ses coordonnées (x,y), mesurées en
degrés
y=90° à -90° Nord->Sud
X=0° à 360° ou 180 à -180 Est->Ouest

Ex. localisation de CRANmirrors

CRANdf<-read.table("CRAN051001a.txt",header=TRUE)
CRAN_mat<-cbind(CRANdf$long,CRANdf$lat)
row.names(CRAN_mat)<-1:nrow(CRAN_mat)
crs3<-CRS("+proj=longlat +ellps=WGS84")
CRAN_sp<-SpatialPoints(CRAN_mat,proj4string=crs3)
summary(CRAN_sp)
Object of class SpatialPoints
Coordinates:
min max
coords.x1 -122.95000 153.0333
coords.x2 -37.81667 57.0500
Is projected: FALSE
proj4string : [+proj=longlat +ellps=WGS84]
Number of points: 54
(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 10
Pour retrouver les caractéristique d’un objet Spatial:
bbox(CRAN_sp) min max
proj4string(CRAN_sp) coords.x1 -122.95000 153.0333
coords.x2 -37.81667 57.0500
plot(CRAN_sp) [1] "+proj=longlat +ellps=WGS84"

wrld<-map("world")
plot(CRAN_sp,add=T,col="blue")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 11

 Extraire les coordonnées d’un point particulier
bresil<-which(CRANdf$loc=="Brazil")
bresil [1] 4 5 6 7 8
coordinates(CRAN_sp)[bresil,]
coords.x1 coords.x2
4 -49.26667 -25.41667
5 -42.86667 -20.75000
6 -43.20000 -22.90000
7 -47.63333 -22.71667
summary(CRAN_sp[bresil,]) 8 -46.63333 -23.53333
Object of class SpatialPoints
Coordinates:
min max
coords.x1 -49.26667 -42.86667
coords.x2 -25.41667 -20.75000
Is projected: FALSE
proj4string : [+proj=longlat +ellps=WGS84]
Number of points: 5

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 12

II.3 SpatialLines Class "Line" [package "sp"]
Slots:
getClass("Line") Name: coords
Class: matrix
Known Subclasses: "Polygon"
getClass("SpatialLines")
Class "SpatialLines" [package "sp"]
Slots:
Name: lines bbox proj4string
Class: list matrix CRS
Extends: "Spatial", "SpatialLinesNULL"
Known Subclasses: "SpatialLinesDataFrame"
Exemple Japon

library(maps)
Japan<-map("world","japan",plot=FALSE)
proj<-CRS("+proj=longlat +ellps=WGS84")
Sp_japan<-map2SpatialLines(japan,proj4string=proj)

plot(Sp_japan)

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 13

II.4 SpatialPolygons
getClass("Polygon")
Class "Polygon" [package "sp"]
Slots:
Name: labpt area hole ringDir coords
Class: numeric numeric logical integer matrix
Extends: "Line"
getClass("SpatialPolygons")

Class "SpatialPolygons" [package "sp"]

Slots:
Name: polygons plotOrder bbox proj4string
Class: list integer matrix CRS
Extends: "Spatial", "SpatialPolygonsNULL"
Known Subclasses: "SpatialPolygonsDataFrame"

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 14

II.5 SpatialGrid
getClass("GridTopology")
Class "GridTopology" [package "sp"]
Slots:
Name: cellcentre.offset cellsize cells.dim
Class: numeric numeric integer

getClass("SpatialGrid")

Class "SpatialGrid" [package "sp"]

Slots:
Name: grid bbox proj4string
Class: GridTopology matrix CRS
Extends: "Spatial"
Known Subclasses: "SpatialGridDataFrame"

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 15

 III. Visualiser des données spatiales
Commandes plot et image
Exemple: bassin de la Meuse (pollutions par les métaux lourds)
Points
library(sp)
data(meuse)
coordinates(meuse)<-c("x","y")
plot(meuse)
title("Points")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 16

 Lignes
coord.meuse<-coordinates(meuse)
list.l.meuse<-list(Lines(list(Line(coord.meuse)),"1"))
lines.meuse<-SpatialLines(list.l.meuse)
plot(lines.meuse)
title("Lines")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 17

 Polygones
data(meuse.riv)
list.p.meuse<-list(Polygons(list(Polygon(meuse.riv)),"meuse.riv"))
pol.meuse<-SpatialPolygons(list.p.meuse)
plot(pol.meuse,col="blue")
title("Polygones")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 18

 Grilles

data(meuse.grid)
coordinates(meuse.grid)<-c("x","y")
meuse.grid<-as(meuse.grid,"SpatialPixels")
image(meuse.grid,col="grey")
title("Grille")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 19

 Combinaison des éléments

image(meuse.grid,col="grey")
plot(pol.meuse,col="blue",add=TRUE)
plot(meuse,add=TRUE)

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 20

 Combinaison des éléments, avec un attribut

Ex. contamination par le zinc

image(meuse.grid,col="grey")
plot(pol.meuse,col="blue",add=TRUE)
plot(meuse, pch=1,cex=sqrt(meuse$zinc)/20,add=TRUE)

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 21

IV. Interpolation et Géostatistique
Observations en un nombre limité de localisations=> échantillons spatiaux
Ex: particules de l’air, teneur aurifère …

Objectif: connaitre la distribution spatiale de la variable, à partir d’échantillons spatiaux

Géostatistique: analyse de champs aléatoires Z(s), avec Z une VA aléatoire sur une
localisation s non aléatoire.
Estimer et prédire Z(s) sur une région définie, en tenant compte de l’autocorrélation
spatiale

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 22

Exemple contamination par le zinc, Meuse
spplot(meuse,"zinc",do.log=TRUE,add=T)

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 23

Contamination fonction de la distance à la Meuse
xyplot(log(zinc)~sqrt(dist), as.data.frame(meuse))

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 24

 Contamination fonction de la distance à la Meuse : modèle linéaire
zn.lm <- lm(log(zinc)~sqrt(dist), meuse)
meuse$estim<-predict(zn.lm, meuse)
meuse$res <- zn.lm$res
spplot(meuse, "estim",main="Estimations")
spplot(meuse, "res",main="Résidus")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 25

 Interpolation:

1/création de la surface : grille d’interpolation

data(meuse.grid)
coordinates(meuse.grid) <- c("x", "y")
meuse.grid <- as(meuse.grid, "SpatialPixelsDataFrame")

2/Choix de la méthode d’interpolation:

Méthodes non-géostatistiques
-Inverse Distance Weighted Interpolation
-Régression linéaire

Méthodes géostatistiques:
-Variogramme
-Modèle additif

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 26

 IV.1 Inverse Distance Weighted Interpolation:
Prédiction par une moyenne pondérée de la distance
n

∑ w (s )Z(s ) i i
Ẑ(s 0 ) = i =1
n

∑ w (s )
i =1
i

Poids:
1
w (s i ) =
d(s i − s 0 )
p

distance
p: paramètre de puissance.
Le choix de p, fonction
-degré de lissage désiré pour l'interpolation,
-densité et distribution des échantillons interpolés,
-distance maximum d’influence des points environnants.
grandes valeurs => influence++ valeurs proches du point interpolé.
Pour p > 1 le pic plus pointu au-dessus du point interpolé,
Pour 0 < p < 1 => pics lissés au-dessus du point interpolé.
(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 27
idw.meuse<-idw(zinc~1,meuse,meuse.grid,idp=2.5)
spplot(idw.meuse,"var1.pred",main="idw pred, idp=2.5")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 28

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 29
 IV.2 Interpolation par Régression linéaire:
Ẑ(s 0 ) = α + β d(s 0 − s i ) + ε

zn.lm <- lm(log(zinc)~sqrt(dist), meuse)

meuse.grid$pred <- exp(predict(zn.lm, meuse.grid))
spplot(meuse.grid,"pred",main="reg lin")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 30

  − d(s 0 − s i ) 
Ẑ(s 0 ) = α + β exp +ε
 δ 
zn.lm2 <- lm(log(zinc)~exp(-dist/0.2), meuse)
meuse.grid$pred2 <- exp(predict(zn.lm2, meuse.grid))
spplot(meuse.grid,"pred2",main="reg lin, exp,0.2")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 31

 Ẑ(s 0 ) = α + β1 x + β 2 y + β 3 x ² + β 4 y ² + β 5 x × y + ε
zn.lm3 <- lm(log(zinc)~x+y+I(x^2)+I(y^2)+I(x*y), meuse)
meuse.grid$pred3 <- exp(predict(zn.lm3, meuse.grid))
spplot(meuse.grid,"pred3",main="reg lin, polynome2")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 32

 IV.3 Variogramme
⇒Corrélation entre 2 observations à une distance donnée

Rappel: série temporelle, plusieurs observations espacées d’un lag k

Ex: t=0 ↔ t=k; t=10↔t=10+k …

Mais données spatiales, ce n’est pas toujours le cas: une distance de h peut être observée
une seule fois.
=> Hypothèse de stationnarité obligatoire
Z(s ) = m + ε(s ) E[Z(s )] = m = cst ∀ s
Variogramme (ou semi-variogramme)

1
γ (h ) = E[Z(s ) − Z(s + h )]
2

=> Variance constante et variogramme indépendant de s

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 33

Si on suppose que le variogramme est indépendant de la direction (isotropie),
Et si on regroupe des lags similaires => Nh paires

Estimation du variogramme

( )
Nh
~ 1 ~
γˆ h j = ∑ [Z(s ) − Z(s + h )]2
∀h ∈ h j
2N h i =1

Le modèle linéaire devient Stationnaire,

Variogramme à modéliser
Z(s ) = α + βX + ε(s )

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 34

Exploration du variogramme
hscat(log(zinc)~1,meuse,(0:9)*100)

?interprétez?

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 35

plot(variogram(log(zinc)~1,meuse,cloud=TRUE)

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 36

plot(variogram(log(zinc)~1,meuse,cloud=TRUE)

Quels sont les

paires
correspondantes?

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 37

sele<-
plot(variogram(log(zinc)~1,meuse,cloud=TRUE),digitize=TRUE)
plot(sele,meuse)

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 38

Forme du variogramme
plot(variogram(log(zinc)~1,meuse,cloud=FALSE)

Attention:
Nombre de paires

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 39

Forme du variogramme dans 4 directions (anisotropie)

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 40

Modélisation du variogramme
show.vgms()

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 41

show.vgms(model="Sph",nugget=0.2,range=60,sill=3,max=100)

Sill=plateau

Nugget=pépite
Range=portée

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 42

Variogramme « pépite pure » (~bruit blanc)
show.vgms(model="Nug",nugget=0.2,sill=1,max=100)

C, si h > 0
γ (h ) = 
 0 si h = 0

C=pépite

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 43

Variogramme Exponentiel
show.vgms(model="Exp",nugget = 0.2,sill=3,range=30,max=100)

 − 
h

γ (h ) = C1 − e a 
 

C=pépite
a=portée

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 44

Variogramme Gaussien
show.vgms(model="Gau",nugget=0.2,sill=3,range=30,max=100)

 
2
h
− 
γ (h ) = C1 − e  a  
 
 

C=pépite
a=portée

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 45

 Variogramme Sphérique
show.vgms(model="Sph",nugget=0.2,sill=3,range=30,max=100)

 C, si h ≥ a

γ (h ) =   3 h 1  h  
3

C − si 0 ≤ h ≤ a
  2 a 2  a  
  

C=pépite
a=portée

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 46

Variogramme Linéaire
show.vgms(model="Lin", nugget=0.2,sill=3,range=30,max=100)

h
γ (h ) = C
a

C=pépite
a=portée

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 47

Variogramme Puissance
show.vgms(model="Pow", nugget=0.2,sill=3,range=2,max=100)

b
h
γ (h ) = C 
a

C=pépite
a=portée

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 48

Variogramme Matérn
show.vgms(model="Mat",nugget=0.2,sill=3,range=30,max=100)

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 49

 Modèles de variogramme les plus utilisés:
Sphérique, Gaussien, Exponentiel, Matérn, Puissance

1. Choisir un modèle, en fonction du variogramme observé

2. Choisir des valeurs initiales: portée, pépite, plateau
3. Modéliser les observations par le modèle choisie

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 50

 Modèle sphérique
svgm <- variogram(log(zinc) ~ 1, meuse)
v.fit<-fit.variogram(svgm, vgm(1, "Sph", 800, 1))
plot(svgm, v.fit, pch = 3)

attr(v.fit, "SSErr")
[1] 9.011194e-06

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 51

 Modèle Gaussien
svgm <- variogram(log(zinc) ~ 1, meuse)
v.fit2<-fit.variogram(svgm, vgm(1, "Gau", 800, 1))
plot(svgm, v.fit2, pch = 3)

attr(v.fit2, "SSErr")
[1] 1.915069e-05

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 52

 Interpolation: Krigeage Ordinaire:
meuse.zinc.ok<-krige(log(zinc)~1,meuse,meuse.grid,v.fit)
summary(meuse.zinc.ok)

Object of class SpatialPixelsDataFrame

Coordinates:
min max
x 178460 181540
y 329620 333740
Is projected: NA
proj4string : [NA]
Number of points: 3103
Grid attributes:
cellcentre.offset cellsize cells.dim
x 178460 40 78
y 329620 40 104
Data attributes:
var1.pred var1.var
Min. :4.777 Min. :0.08549
1st Qu.:5.238 1st Qu.:0.13728
Median :5.573 Median :0.16218
Mean :5.707 Mean :0.18533
3rd Qu.:6.172 3rd Qu.:0.21161
Max. :7.440(c)2013,Max. :0.50028
Jean Gaudart Aix-Marseille Univ, UMR912 53
spplot(meuse.zinc.ok,"var1.pred",main="Ord. Krig. Sph")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 54

Krigeage Universel:
meuse.zinc.uk<-krige(log(zinc)~sqrt(dist),meuse,meuse.grid,v.fit)
spplot(meuse.zinc.uk,"var1.pred",main="Univ. Krig. Sph")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 55

Krigeage Universel (suite):
meuse.uk2<-krige(log(zinc)~exp(-dist/25),meuse,meuse.grid,v.fit)
spplot(meuse.uk2,"var1.pred",main="Univ. Krig. Sph")

?Que remarquez vous?

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 56

Diagnostique de la modélisation: validation croisée
meuse.cv<-krige.cv(log(zinc)~1,meuse,v.fit,nfold=5)
bubble(meuse.cv, "residual", main = "5-CV résidus")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 57

meuse.cv2<-
krige.cv(log(zinc)~sqrt(dist),meuse,v.fit,nfold=5)
bubble(meuse.cv2, "residual", main = "5-CV résidus")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 58

 V. Analyse de Co-facteur
Contrairement à l’interpolation, l’objectif ici est d’étudier le lien entre des facteurs de
risque X et une variable à expliquer Y, en tenant compte de la spatialisation

Approche par les modèles GAM

• GLM: µ i = E (Yi ) avec g (µ i ) = X i β

Loi de la famille exponentielle (Gaussian, Poisson, Binomial …)
• Modèle Additif Généralisé (GAM)
g (µ i ) = X i β + f1 (x1,i ) + f 2 (x 2,i ) + f 3 (x1,i , x 2,i )

où les f(.) sont des fonctions de lissage (splines)

=> modélise des relations non linéaires
• Si x1 et x2 sont les coordonnées de chaque localisation
g (µ i ) = X i β + f (x1,i , x 2,i )

=> modélise la variabilité spatiale des Yi

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 59

• « Thin plate splines » : Splines à plaque mince
– Estimer une fonction de plusieurs prédicteurs
– Sans besoin de spécifier les nœuds ni la base
– Recherche une surface lissée la moins déformée possible:
=>minimiser le degré de courbure tout en passant à proximité des observations

2 2 2
2 ∂ f
2
  ∂ f 2
 ∂ f 2

y− f + λ ∫∫  2  + 2  +  2  dx1 dx 2
 ∂x1   ∂x1 ∂x 2   ∂x 2 

• Logiciel R, package mgcv

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 60

 install.packages("mgcv")
library(mgcv)
summary(meuse)
mod.lm1<-lm(log(zinc)~dist+elev+soil, data=meuse)
summary(mod.lm1)

Call: lm(formula = log(zinc) ~ dist + elev + soil, data = meuse)

Residuals:
Min 1Q Median 3Q Max
-1.05511 -0.26034 0.01597 0.25479 1.08302
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.49161 0.30166 28.149 < 2e-16 ***
dist -1.96409 0.24341 -8.069 2.11e-13 ***
elev -0.25911 0.03951 -6.557 8.31e-10 ***
soil2 -0.09565 0.09167 -1.043 0.298
soil3 0.12545 0.16560 0.758 0.450
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4273 on 150 degrees of freedom
Multiple R-squared: 0.6587, Adjusted R-squared: 0.6496
F-statistic: 72.38 on 4 and 150 DF, p-value: < 2.2e-16
(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 61
mod.lm2<lm(log(zinc)~
dist+elev+soil+x+y+I(x^2)+I(y^2)+I(x*y),data=meuse)
summary(mod.lm2)

Residuals:
Min 1Q Median 3Q Max
-0.90671 -0.23200 0.02764 0.22892 1.08983
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.894e+03 5.404e+03 0.721 0.47234
dist -1.005e+00 3.523e-01 -2.854 0.00496 **
elev -2.793e-01 3.924e-02 -7.118 4.68e-11 ***
soil2 -1.033e-01 1.246e-01 -0.829 0.40846
soil3 -3.598e-02 1.854e-01 -0.194 0.84643
x 1.224e-01 4.560e-02 2.684 0.00812 **
y -9.005e-02 4.975e-02 -1.810 0.07232 .
I(x^2) 3.536e-07 1.928e-07 1.834 0.06870 .
I(y^2) 3.412e-07 1.320e-07 2.585 0.01072 *
I(x * y) -7.545e-07 2.650e-07 -2.847 0.00505 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4103 on 145 degrees of freedom
Multiple R-squared: 0.6958, Adjusted R-squared: 0.6769
F-statistic: 36.84 on (c)2013,
9 and 145 DF, p-value: < 2.2e-16
Jean Gaudart Aix-Marseille Univ, UMR912 62
mod.gam1<-gam(log(zinc)~
dist+elev+soil, family="gaussian", data=meuse)
summary(mod.gam1)

Family: gaussian
Link function: identity

Formula:
log(zinc) ~ dist + elev + soil

Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.49161 0.30166 28.149 < 2e-16 ***
dist -1.96409 0.24341 -8.069 2.11e-13 ***
elev -0.25911 0.03951 -6.557 8.31e-10 ***
soil2 -0.09565 0.09167 -1.043 0.298
soil3 0.12545 0.16560 0.758 0.450
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) = 0.65 Deviance explained = 65.9%

GCV score = 0.18867 Scale est. = 0.18258 n = 155
(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 63
mod.gam1p<-gam(zinc~
dist+elev+soil, family="poisson", data=meuse)
summary(mod.gam1p)

Family: poisson
Link function: log

Formula:
zinc ~ dist + elev + soil

Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.494003 0.029794 285.090 <2e-16 ***
dist -2.513710 0.031208 -80.546 <2e-16 ***
elev -0.232647 0.004036 -57.639 <2e-16 ***
soil2 -0.146465 0.011446 -12.796 <2e-16 ***
soil3 -0.034154 0.023993 -1.424 0.155
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) = 0.612 Deviance explained = 67.4%

UBRE score = 78.902 Scale est. = 1 n = 155
(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 64
mod.gam2p<-gam(zinc~dist+elev+soil
+x+y+I(x^2)+I(y^2)+I(x*y), family="poisson", data=meuse)
summary(mod.gam2p)

Family: poisson Link function: log

Formula:
zinc ~ dist + elev + soil + x + y + I(x^2) + I(y^2) + I(x * y)
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.000e+00 0.000e+00 NA NA
dist -1.584e+00 4.249e-02 -37.285 < 2e-16 ***
elev -2.666e-01 4.316e-03 -61.779 < 2e-16 ***
soil2 -1.862e-01 1.401e-02 -13.297 < 2e-16 ***
soil3 -3.000e-01 3.002e-02 -9.993 < 2e-16 ***
x 1.329e-01 4.684e-03 28.368 < 2e-16 ***
y -7.221e-02 2.544e-03 -28.385 < 2e-16 ***
I(x^2) 1.501e-07 2.133e-08 7.036 1.99e-12 ***
I(y^2) 2.625e-07 9.306e-09 28.205 < 2e-16 ***
I(x * y) -5.644e-07 2.559e-08 -22.060 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.632 Deviance explained = 70.6%
UBRE score = 71.18 Scale est. = 1 n = 155
(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 65
mod.gam3p<-gam(zinc~dist+elev+soil
+s(x,y), family="poisson", data=meuse)
summary(mod.gam3p)

Family: poisson Link function: log

Formula: zinc ~ dist + elev + soil + s(x, y)
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.629247 0.077135 111.872 < 2e-16 ***
dist -2.029306 0.240797 -8.427 < 2e-16 ***
elev -0.265391 0.005592 -47.458 < 2e-16 ***
soil2 -0.126823 0.020675 -6.134 8.57e-10 ***
soil3 -0.293970 0.041853 -7.024 2.16e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:

edf Ref.df Chi.sq p-value
s(x,y) 28.82 29 8240 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.848 Deviance explained = 88.7%
UBRE score = 27.05 Scale est. = 1 n = 155
(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 66
SIR associés aux cofacteurs:

exp(mod.gam3p$coeff[2:5])

dist elev soil2 soil3

0.1314267 0.7669060 0.8808892 0.7452989

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 67

 VI. Recherche de Cluster
Méthodes
– autour d’un facteur de risque identifié, géo-localisé
ex. centrales nucléaires / leucémies

– Méthodes générales
recherche d’agrégats sans facteur de risque géo-localisé
• Globales – « clustering »
• Locales

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 68

• Détection globale
– Le pattern observé globalement est-il concordant avec
l'hypothèse nulle ? 1 test
Moran
• corrélations
• comparaisons de distributions Tango

• Détection locale
– Y a-t-il localement un excès de cas ? k tests
• corrélations (voisinage)
• comparaisons de groupes de proximités LISA

Satscan

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 69

Hypothèse
3 questions possibles :
– Les distributions des cas de chaque zone sont-elles
indépendantes?

– Quelles sont ces distributions ?

– Le risque est-il constant sur l’ensemble de la zone

géographique?

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 70

Hypothèses nulles possibles:
– Distribution des cas Uniforme sur le plan
Complete Spatial Randomness
– Distribution de Poisson hétérogène
Constant Risk Hypothesis

Ei = λni
nombre de Effectif région i
cas attendus
dans la région i
risque constant
k
nombre total de cas observés : O+ = ∑ Oi
~ O+
λ= k
i =1

n+ effectif total : n+ = ∑ ni
i =1
(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 71
 VI.1 Méthodes Globales
Coefficient de Moran (1917-1988)

• Rappel Coefficient de corrélation de Pearson

n

∑ (X i − X )(Yi − Y )
r= i
n n

∑ (X −X) ∑ (Y − Y )
2 2
i i
i =1 i =1

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 72

• Coefficient de corrélation
• Pondéré par les distances
• Similarités entre régions:
écart à la moyenne de la région i
⇔ écart à la moyenne de la région j

K
K × ∑ wij (Yi − Y )(Y j − Y )
Plus les zones i et j sont
i, j
I= K
éloignées, moins le poids

w+ × ∑ (Yi − Y )
2 est important

i =1

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 73

 Préparation pour l’analyse
library(DCluster)

charger les
SPAT<-read.csv2("D:/Mon_dossier/DATA1.csv",header=TRUE)
attach(SPAT) données

TAB<-data.frame(Observed=cas) construire le tableau

TAB<-cbind(TAB, Expected=n*sum(cas)/sum(n)) de données adapté
SPAT.xy = SPAT[c("x", "y")]
coordinates(SPAT.xy) <- c("x","y") construire le tableau
class(SPAT.xy) des coordonnées

coords<-coordinates(SPAT.xy) construire le tableau

des distances
dlist<-dnearneigh(coords, 0, Inf)
et la matrice des poids
dlist<-include.self(dlist)
dlist.d<-nbdists(dlist, coords)
col.W<-nb2listw(
dlist,glist=lapply(dlist.d,function(x){exp(-x)}),style="C")

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 74

 Estimation du coefficient de Moran
I<-
moranI.stat(TAB, listw=col.W, n=length(dlist), S0=Szero(col.W) )

[1] 0.356332

k
S 0 = w+ = ∑ wi
i =1

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 75

Test:
H0: I=0; malades spatialement indépendants
H1: I>0
I − E (I )
– Sous H0: z= ∼>N(0;1)
var(I )
– Condition : distribution des Y normale
• insoutenable
• unité statistique: région, avec k petit

– Loi de I inconnue ⇒ Monte-Carlo ou Bootstrap

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 76

 test du coefficient de Moran

IT<-moranI.test(Observed~offset(log(Expected)), TAB,
model="poisson", R=999, listw=col.W, n=length(dlist),
S0=Szero(col.W))

Moran's I test of spatial autocorrelation

Type of boots.: parametric

Model used when sampling: Poisson
Number of simulations: 999
Statistic: 0.356332
p-value : 0.324

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 77

plot(IT)
I observé

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 78

Interprétation :
 I >0: régions voisines: mêmes écarts à la moyenne = pattern
sous forme de clusters
 I <0: régions voisines: ≠ écarts à la moyenne, = pattern régulier.
 I =0: aucune corrélation spatiale.

 Mesure de l'écart à la moyenne générale:

pas d'interprétation locale possible.
(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 79
 Distances et Poids

• Distances entre localisations (ou centres des régions)

2 2
d ij = ( x j − xi ) + ( y j − yi )
• Différents poids = différents résultats !!

1, si d ij < δ  d ij 
wij =   − 
 0, si non wij = e  τ 

wij = d ij nv

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 80

 distance

45,00
40,00
35,00
30,00
25,00
20,00
15,00
10,00

Exemple: 5,00
0,00
Pour le point le 1 19 37 55 73 91 109 127 145 163 181 199 217 235 253 271 289

plus à l’ouest 1
n/5d
exp (-d/5)
exp(-d/25)
0,8 seuil

0,6

0,4

0,2

0
1 19 37 55 73 91 109 127 145 163 181 199 217 235 253 271 289

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 81

 Méthodes Locales
LISA de Anselin
Local Indicator of Spatial Autocorrelation
Application locale du coefficient de Moran
A proximité d'un cas observé, les cas sont-ils regroupés?

Global Local: pour la région i

K
K × ∑ wij (Yi − Y )(Y j − Y ) K
I=
i, j
K
I i = (Yi − Y )× ∑ wij (Y j − Y )
w+ × ∑ (Yi − Y )
2
j =1
i =1

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 82

Test: Pour chaque région
H0: Ii=0 indépendance des régions/ voisins
H1: Ii>0 I i − E (I i )
z=
var(I i )
– Sous H0 et condition N : ∼>N(0;1)

Condition insoutenable:
=>Loi de Ii inconnue => Monte-Carlo ou Bootstrap

α
– Tests multiples et corrélés αi =
=> correction de Bonferroni nv

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 83

 plot(IT)

dlist2<-dnearneigh(coords, 0, Inf)
dlist2.d<-nbdists(dlist2, coords)

LM<-localmoran(TAB$Observed, nb2listw(dlist2),
p.adjust.method="bonferroni")

Ii E.Ii Var.Ii Z.Ii Pr(z>0)

1 -4.328844e-03 -0.003344482 1.720577e-05 -0.2373111 1
2 -5.323609e-03 -0.003344482 1.720577e-05 -0.4771299 1
3 -1.009687e-02 -0.003344482 1.720577e-05 -1.6278736 1
4 -1.284904e-05 -0.003344482 1.720577e-05 0.8031933 1

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 84

Interprétation
– Ii >0: régions voisines similaires à région i =
cluster local
– Ii <0: régions voisines: ≠ à région i , =
région i particulière.
– Ii =0: aucune corrélation spatiale entre la région i
et ses voisins.

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 85

 Sites
•RS. Bivand, EJ. Pebesma, V. Gomez-Rubio, Applied Spatial Data Analysis with R,
eds. Springer, Use R!

•L. Waller, C. Gotway, Applied Spatial Statistics for Public Health Data, eds. Wiley

•Open Source Geospatial Fondation www.osgeo.org

•Geospatial Data Abstraction Library www.gdal.org
•Cartographic Projections library trac.osgeo.org/proj/
•http://www.asdar-book.org/
•http://trac.osgeo.org/

(c)2013, Jean Gaudart Aix-Marseille Univ, UMR912 86