Cours2

II-‐ Autocorréla-on Spa-ale Généralités •  Lorsqu’on étudie une variable géoréférencée {Y (s), s ∈ S } , ce qui se
passe dans une localisation particulière dépend de ce qui se passe dans
d’autres localisations. Ces interactions sont d’autant plus fortes que les
localisations concernées sont plus proches.
•  Autocorrélation spatiale= corrélation d’une variable géoréférencée avec
elle-même.
Permet de mesurer le degré de ressemblance entre observa-ons voisines •  L’ Autocorrélation spatiale peut être:
–  Positive: les valeurs de la variable dans des lieux proches se
ressemblent et dans des lieux éloignés sont dissemblables
–  Négative dans le cas contraire (phénomène d’évitement)
Généralités Généralités • 
Dans la suite du cours, on s’intéresse à des données laticielles : S est discret
fixé à n points (par exemple, S est un ensemble de zones géographiques)
• 
On note (y1,…yn) les données dont on dispose sur le phénomène Y
sur les n points de S. On considère qu’il s’agit d’une réalisation du vecteur
aléatoire {Y (s), s ∈ S } .
• 
On étudie la loi de {Y (s), s ∈ S }
Généralités •  Exemple : données sids du package spdep (mort subite du nourrisson en caroline du nord) ü  Données de type "Spa-alPolygonsDataFrame" . Ce type d’objet peut contenir à la fois des données de localisa-ons géographiques et des données classiques sous la forme d’un tableau de données. Les localisa-ons sont repérées par le couple longitude-‐la-tude et u-lisent ici la référence NAD27. nc_ﬁle=system.ﬁle("etc/shapes/sids.shp", package="spdep")[1] llCRS=CRS("+proj=longlat +datum=NAD27") nc.sids = readShapeSpa-al(nc_ﬁle, ID="FIPSNO", proj4string=llCRS) ü  Le format .shp présuppose qu’on a 3 ﬁchiers : *.shp, *.shx, and *.dbf, Le premier con-ent les données géométriques (cartographie), le deuxième les index des enregistrements du ﬁchier .shp, le troisième les a_ributs. Les centres des régions sont accessible avec la fonc-on coordinates(), et peuvent être u-lisées pour placer des labels sur la carte. Généralités Contenu de nc.sids : SP_ID : Spa-alPolygons ID CNTY_ID : county ID east : eas-ngs, county seat, miles, local projec-on north : northings, county seat, miles, local projec-on L_id : Cressie and Read (1985) L index M_id : Cressie and Read (1985) M index names : County names AREA : County polygon areas in degree units PERIMETER : County polygon perimeters in degree units CNTY_ : Internal county ID NAME : County names FIPS : County ID FIPSNO : County ID CRESS_ID : Cressie papers ID BIR74 : births, 1974-‐78 SID74 : SID deaths, 1974-‐78 NWBIR74 : non-‐white births, 1974-‐78 BIR79 : births, 1979-‐84 SID79 : SID deaths, 1979-‐84 NWBIR79 : non-‐white births, 1979-‐84 names(nc.sids) [1] "SP_ID" "CNTY_ID" "east" "north" "L_id" "M_id" "names" "AREA" "PERIMETER" "CNTY_" "NAME" "FIPS" "FIPSNO" "CRESS_ID" "BIR74" [16] "SID74" "NWBIR74" "BIR79" "SID79" "NWBIR79" Généralités Contours des unités géographiques (comtés) plot(nc.sids, axes=TRUE) text(coordinates(nc.sids), label=nc.sids$FIPSNO, cex=0.5) -tle("Fron-ères des comtés, Caroline du nord") 37005
37073 37029
37009
37131
37171 37169 37157 3703337145
37053
37185
37181
37091
37139
37077
37083
37189
37143
37193 37197
37041
37067 37081
37069
37011
37015
37135
37001
37063
37121
37127 37065
3702737003
37059
37115 37199
37117
3718737177
37097
37057
37183
37055
37023
37151
37195
37037
37111
37035
37159
37147
37013
37087 37021
37095
37101
37079
37173
37109
37105
37161
37025
37085
37191
37075
37089
37045
37123
37167
37125
37071
37099
37149
37119
37107 37049
37175
37113
37137
37039
37043
37103
37093 37051 37163
37179 3700737153
37061
37165
37031
37133
37155
37017
37141
37047
34°N
35°N
36°N
37°N
38°N
Frontières des comtés, Caroline du nord
37129
32°N
33°N
37019
84°W
82°W
80°W
78°W
76°W
Généralités •  Les diﬀérentes étapes d’une étude de l’autocorréla@on spa@ale ü  Choix d’un critère de voisinage : Déterminer les aires liées entre elles (dépend de la connaissance du problème étudié) ü  Choix d’ une matrice de poids : Assigner des poids aux aires qui sont liées ü  Mise en oeuvre de tests sta-s-ques en u-lisant ce_e matrice de poids, pour examiner l’autocorréla-on spa-ale. No-on de voisinage • 
Généralités ü  Pour mesurer la ressemblance entre lieux proches, il faut déﬁnir une no-on de voisinage entre localisa-ons ́de S. ü  On peut déﬁnir diﬀérentes no-ons de voisinage, selon l’idée que l’on se fait du problème étudié, qui impacteront le résultat des analyses de l’autocorréla-on. ü  Un type donné de voisinage pourra être représenté par un graphe, appelé graphe de voisinage : à chaque observa-on spa-ale (zone géographique par ex ), on associe un nœud du graphe (centre de la zone). Les nœuds correspondant à deux zones voisines sont reliés par un arc. No-on de voisinage Exemple : Ajout du système de voisinage aux données nc.sids (voisinage proposé par Cressie et Chan (1989), stocké dans le ﬁchier ncCC89.gal) •  Créa-on de l’objet ncCC89 de classe neigborhood gal_ﬁle=system.ﬁle("etc/weights/ncCC89.gal", package="spdep")[1] ncCC89=read.gal(gal_ﬁle, region.id=nc.sids$FIPSNO) Les systèmes de voisinage sont généralements stockés dans des ﬁchiers .gal lus sous R avec la fonc-on read.gal(). L’objet ainsi crée est un objet de classe « nb » (neigborhood). Il s’agit d’une liste de listes, une pour chaque observa-on, contenant une liste de voisins et d’autres informa-ons, comme l’iden-ﬁant de la région. Contenu résumé de l’objet ncCC89 est un objet de classe “nb” summary(ncCC89) Neighbour list object: Number of regions: 100 Number of nonzero links: 394 Percentage nonzero weights: 3.94 Average number of links: 3.94 2 regions with no links: 37055 37095 Link number distribu-on: 0 1 2 3 4 5 6 7 8 2 3 14 18 28 21 7 6 1 3 least connected regions: 37019 37031 37133 with 1 link 1 most connected region: 37041 with 8 links No-on de voisinage Créa@on du graphe de voisinage plot(nc.sids, border="grey") plot(ncCC89, coordinates(nc.sids), add=TRUE, col="blue") No-on de voisinage • 
Les diﬀérentes no@ons de voisinage v  Con@güité ü  No-on fréquemment u-lisée lorsque les données portent sur des zones géographiques. ü  Distance de con-güité entre deux zones = nombre minimal de fron-ères qu’il faut franchir pour aller de l’intérieur de l’une à l’intérieur de l’autre. ü  La no-on de con-güité repose sur celle de fron-ère partagée q  Deux zones sont dites con-gües à l’ordre k (lag k) si leur distance de con-güité est égale à k (analogue spa-al de la déﬁni-on des retards (lag) en séries temporelles) q  CP: deux zones sont dites con-gües (à l’ordre 1) si elles ont une fron-ère commune. 2/15/20
No-on de voisinage Spatial weights matrices
 Neighborhoods
can be defined in a number of ways
ü  Déﬁni-on d’une fron-ère partagée 
Contiguity (common boundary)


Distance (distance band, K-nearest neighbors)


What is a “shared” boundary?
How many “neighbors” to include, what distance do we use?
General weights (social distance, distance decay)
La zone centrale La zone centrale La zone centrale a 4 voisins dans a 4 voisins dans a 8 voisins dans les direc-ons les direc-ons les direc-ons N-‐S et E-‐O N-‐E/S-‐O et N-‐O/S-‐E N-‐S, E-‐O N-‐E/S-‐O, N-‐O/S-‐E Common weights measures

Most common is using binary connectivity based on
No-on de voisinage ü  Exemple : voisins de 5? 1
4
7
• 
• 
• 
2 3
5 6
8 9
Voisinage de type « rook » : 2,4,6,8 Voisinage de type « bishop »: 1,3,7,9 Voisinage de type « queen » 1,2,3,4,6,7,8,9. No-on de voisinage ü  Exemple Créa+on d’un voisinage de type queen sids_nbq=poly2nb(nc.sids) ; sids_nbq Neighbour list object: Number of regions: 100 Number of nonzero links: 490 Percentage nonzero weights: 4.9 Average number of links: 4.9 Créa+on d’un voisinage de type rook sids_nbr=poly2nb(nc.sids, queen=FALSE) Sids_nbr Neighbour list object: Number of regions: 100 Number of nonzero links: 462 Percentage nonzero weights: 4.62 Average number of links: 4.62 Graphes de voisinage coords=coordinates(nc.sids) par(mfrow=c(3,1)) plot(nc.sids) plot(sids_nbq, coords, add=T, col=2) legend("bo_om","queen", text.col=2) plot(nc.sids) plot(sids_nbr, coords, add=T, col=3) legend("bo_om","rook", text.col=3) plot(nc.sids) plot(sids_nbq, coords, add=T, col=2) plot(sids_nbr, coords, add=T, col=3) legend("bo_om",c("queen","rook"), text.col=c(2,3)) queen
rook
queen
rook
No-on de voisinage Répar++on des liens de voisinage qq=table(card(sids_nbq)); qq 2 3 4 5 6 7 8 9 8 15 17 23 19 14 2 2 barplot(qq, col=2,main="répar--on du nombre de liens par région") 5
10
15
20
répartition du nombre de liens par région
0
Contenu du ﬁchier “nb” créé summary(sids_nbq) Neighbour list object: Number of regions: 100 Number of nonzero links: 490 Percentage nonzero weights: 4.9 Average number of links: 4.9 Link number distribu-on: 2 3 4 5 6 7 8 9 8 15 17 23 19 14 2 2 8 least connected regions: 37041 37043 37053 37055 37129 37137 37149 37177 with 2 links 2 most connected regions: 37097 37125 with 9 links 2
3
4
5
6
7
8
9
No-on de voisinage Distance based neighbors
k nearest neighbors
 dCan
also choose the k nearest points as neighbors
v  Voisinage basé sur une istance > coords<-coordinates(sids_SP)
> u
IDs<-row.names(as(sids_SP,
"data.frame"))
Nécessite que l’on ait choisi ne distance dij entre lieux i et j (euclidienne, ou autre). > sids_kn1<-knn2nb(knearneigh(coords, k=1), row.names=IDs)
sids_kn2<-knn2nb(knearneigh(coords,
k=2),
ü  Méthode des k voisins >les plus proches (k-‐nearest neighbors) : On row.names=IDs)
choisit comme voisinage les k voisins l>es sids_kn4<-knn2nb(knearneigh(coords,
plus proches de l’observa-on i au sk=4),
ens drow.names=IDs)
e la distance dij. > plot(sids_SP)
k=2
> plot(sids_kn2, coords, add=T)
k=3
k=1
ü  Méthode des voisins à distance radiale : On chosit comme voisins de i les points j situés à moins d’une distance d: 0 ≤ d ij ≤ d
ü  Voisinage global : On choisit l’ensemble de tous les autres points de S comme voisins de i. 1
No-on de voisinage ü  Exemple Créa+on d’un voisinage basé sur les k (=1,2,4) plus proches voisins IDs=row.names(as(nc.sids, "data.frame")) library(RANN) sids_kn1=knn2nb(knearneigh(coords, k=1), row.names=IDs) Neighbour list object: Number of regions: 100 Number of nonzero links: 100 Percentage nonzero weights: 1 Average number of links: 1 Non-‐symmetric neighbours list sids_kn2=knn2nb(knearneigh(coords, k=2), row.names=IDs) sids_kn4=knn2nb(knearneigh(coords, k=4), row.names=IDs) par(mfrow=c(3,1)) plot(nc.sids, main="k=1") plot(sids_kn1, coords, add=T) plot(nc.sids, main="k=2") plot(sids_kn2, coords, add=T) plot(nc.sids, main="k=4") plot(sids_kn4, coords, add=T) k=1
k=2
k=4
Créa+on d’un voisinage à distance radiale Fonc-on dnearneigh (): Iden-ﬁe les voisins d’un point par la distance euclidienne entre une borne inférieure et une borne supérieure •  Fonc-on nbdists (): Etant donnée la liste des voisins de chaque point (donnée dans un ﬁchier de voisinage de classe nb, cf sec-on voisinage) retourne les distances euclidiennes entre le point et chaque voisins via une liste de même forme. nbdists(sids_kn1, coords) sont voisins les points à une distance < à max(dist)
[[1]] [1] 0.2766561 [[2]] [1] 0.2607982 [[3]] [1] 0.2914539 …………... dist=unlist(nbdists(sids_kn1, coords)) summary(dist) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.1197 0.2582 0.2964 0.2910 0.3239 0.4226 sids_kd2=dnearneigh(coords, d1=0, d2=1*max(dist), row.names=IDs) • 
No-on de voisinage • 
Choix per@nent d’un voisinage : • 
Lorsque les zones sont très hétérogènes, les k plus proches voisins ne conviennent pas car les distances entre points peuvent être très dissemblables • 
La con-güité est inadaptée à des données portant sur des points : contrairement au cas des zones, la détermina-on de la con-guïté entre deux points n’a guère de sens. • 
• 
Il faut choisir une méthode ne laissant pas trop de points isolés En réalité, on peut généralement supposer que toutes les observa-ons (et pas seulement celles qui sont les plus proches les unes des autres) interagissent les unes avec les autres, même si leur ac-on s’amoindrit avec la distance. Considérer que le voisinage est S en-er est donc per-nent lorsque S n’est pas trop grand. Matrice de poids • 
Généralités ü  Une matrice de poids permet de quan-ﬁer les rela-ons de voisinage entre points i.e. d’aﬀecter un poids à chaque voisin. Sa structure dépend de la no-on de voisinage choisie ü  Diﬀérents types de matrice de poids selon le voisinage choisi: v  Matrices binaires: Matrice de con-guité, Matrice de distance radiale, Matrice basée sur les k plus proches voisins,… v  Matrice de poids générales: lorsque le voisinage est S tout en-er. Elles sont généralement basées sur une distance. Matrice de poids • 
Forme générale d’une matrice de poids n= nombre de localisa-ons. !
#
#
#
#
#
W =#
#
#
#
##
"
0
... w1i
... 0
wi1 ...
... w1 j
...
0
...
...
...
wij
... ... ...
w j1 ... w ji
0
...
...
0
... ... ...
wn1 ... wni
... ...
... wnj
... w1n $&
... ... &
&
... win &
&
... ... &
... w jn &
&
0 ... &
... 0 &&
%
ü  wij est d’autant plus élevé que j a une inﬂuence importante sur i. n
ü  Généralement la matrice est standardisée en ligne : ∑ wij = 1
j=1
Matrice de poids v  Propriétés : ü  Les matrices de poids sont souvent standardisées en ligne, en par-culier lorsqu’elles sont binaires. -‐ Permet de créer des poids propor-onnels dans le cas où les aires n’ont pas le même nombre de voisins : Une matrice de poids non standardisée fait varier l’inﬂuence des observa-ons. Les zones ayant beaucoup de voisins sont sur-‐ pondérées en comparaison avec celles qui en ont peu. -‐ Facilite les calculs des sta-s-ques de tests et permet de comparer les résultats obtenus avec diﬀérentes matrices de poids ü  Les matrices de poids peuvent être ou non symétriques (la non symétrie est approprié lorsque les rela-ons étudiées sont des rela-ons de diﬀusion ou des rela-ons de type centre -‐ périphérie). Matrice de poids • 
Les diﬀérents types de matrices de poids v  Matrices binaires ²  Matrice de con-guité : !# 0 si i et j ne sont pas contigus
W = (wij ), wij = "
#$ 1 si i et j sont contigus
Con-guité au sens d’une déﬁni-on donnée de la fron-ère commune (rook, queen, bishop). n
ü  Le nombre de régions con-gues à i est alors Li = ∑ wij
j=1
ü  Une matrice de con-guité est généralement symétrique Matrice de poids ²  Matrice de distance radiale " 1 si 0 ≤ d ≤ d
$
ij
W = (wij ), wij = #
$% 0 si d ij > d
ü  généralement symétrique ²  Matrice basée sur les k-‐plus proches voisins : Soit N k (i)
L’ensemble des k voisins les plus proches de i. "$ 1 si j ∈ N (i)
k
W = (wij ), wij = #
$% 0 sinon
ü  généralement pas symétrique Matrice de poids v  Matrices de distances générales ²  Distance puissance ²  Distance exponen-elle W = (wij ), wij = dij− β , β > 0
W = (wij ), wij = exp(−β dij ), β > 0
²  Autres distances : Ex: Les poids reﬂètent l’appartenance de deux individus à la même catégorie sociale. W = (wij ), wij =
1
yi − y j
où yi et yj sont des observa-ons sur des caractéris-ques socio-‐économiques per-nentes (revenu par tête , % de la popula-on dans un groupe racial ou ethnique,…). B 100 10000 490 980 10696
Matrice de poids v  Exemple de calcul d’une matrice de poids Binary vs. row-standardized

A binary weights matrix looks like:





Observation 1 has
Observation 2 has
Observation 3 has
Observation 4 has
neighbor 2
neighbors 3 and 4
neighbors 1 and 2
neighbor 2, 3 and 4
0
1
0
0
0
0
1
1
1
1
0
0
0
1
1
1
A row-standardized matrix it looks like:
0
1
0
0
0
0
.5
.5
.5
.5
0
0
0
.33
.33
.33
Matrice de poids v  Créa@on d’une matrice de poids avec R: nb2listw(neighbours, glist=NULL, style="W", zero.policy=FALSE) ü  Arguments: •  neighbours: un objet de classe ’nb’ (voisinage) •  glist: liste de poids généraux correspondant aux voisins •  style: peut prendre les valeurs W, B, C, U, and S (W si l’on veut une matrice standardisée, B si l’on veut une matrice binaire, le défaut est W) •  zero.policy: si FALSE, arrête avec un message d’erreur dès qu’il y a un point sans voisins. Si TRUE permetde prendre en compte les points isolés (sans voisin). ü 
• 
• 
• 
Renvoie un objet de classe listw, comprenant: style : le style neighbours: la liste des voisins de chaque localisa-on weights : la liste des poids des voisins de chaque localisa-on Matrice de poids v  Exemples : ü  Créa+on d’une matrice de con+guité standardisée (type queen) sids_nbq_w=nb2listw(sids_nbq) Characteris-cs of weights list: Neighbour list object: Number of regions: 100 Number of nonzero links: 490 Percentage nonzero weights: 4.9 Average number of links: 4.9 Weights style: W Weights constants summary: N nn S0 S1 S2 W 100 10000 100 44.65023 410.4746 nn, S0, S1 and S2 (cf précédemment) sont plusieurs indicateurs calculés sur la matrice des poids et u-lisés dans le calcul des sta-s-ques d’autocorréla-ons spa-ales (cf après) Matrice de poids ü  Créa+on d’une matrice de con+guité non standardisée (type queen) sids_nbq_wb= nb2listw(sids_nbq, style="B") sids_nbq_wb$weights Créa+on d’une matrice des deux plus proches voisins non standardisée sids_kn2_w= nb2listw(sids_kn2, style="B") sids_kn2_w$neig [[1]] [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 Neighbour list object: 0.1666667 Number of regions: 100 [[2]] Number of nonzero links: 200 [1] 0.25 0.25 0.25 0.25 [[3]] Percentage nonzero weights: 2 [1] 0.3333333 0.3333333 0.3333333 Average number of links: 2 [[4]] Non-‐symmetric neighbours list [1] 0.25 0.25 0.25 0.25 [[5]] [1] 0.3333333 0.3333333 0.3333333 ……… Chaque sous-‐liste correspond à une localisa-on et ses poids en ligne Matrice de poids ü  Créa+on d’une matrice de distance radiale sids_kd2_w=nb2listw(sids_kd2, style="B") sids_kd2_w$weight [[1]] [1] 1 1 1 1 [[2]] [1] 1 1 1 1 [[3]] [1] 1 1 [[4]] [1] 1 1 1 Matrice de poids • 
Variable décalée (lag) : ü  La i° variable décalée mesure de l’eﬀet global sur la i° observa-on des autres points de l’espace ; Si W est une matrice de poids et Y = (Y1, ..., Yn)’ est le vecteur colonne des n observa-ons de la variable spa-alisée Y(s), la i° variable décalée n
vaut: zi = (WY )i = ∑ wijY j
j=1
On note dans la suite n
S0 = ∑
n
∑w
ij
i=1 i=1,i≠ j
S1 =
n
wi. = ∑ wij
1 n n
(wij + w ji ) 2
∑
∑
2 i=1 i= j
j=1
n
(
S2 = ∑ wi. + w.i
i=1
)
ü  Lorsque la matrice de poids W est une matrice standardisée, la i°variable spa-ale décalée con-ent la moyenne pondérée des valeurs prises par Y sur les régions voisines à la région i. De plus, n n n n et z=
1
1
1
zi = ∑∑ wijY j = ∑Y j = y
∑
n i=1
n j=1 i=1
n j=1
S0 = n
Matrice de poids Le modèle CAR stationnaire
Exemple de construction de W
• 
jj On définit pour tout i un facteur de pondération tel que wij est
grand lorsque i est proche de j .
Exemple : Voisinage "tour" d’ordre 1 :
1
4
7
2
5
8
3
6
9
w51 = w53 = w57 = w59 = 0
w52 = w54 = w56 = w58 = 1/4.
Donc
z5 =
1
(Y2 + Y4 + Y6 + Y8 )
4
1
z1 = (Y2 + Y4 )
2
Matrice de poids ü  Créa-on de la variable SID décalée pour chaque lieu avec matrice de con-guité queen standardisée wy=lag.listw(sids_nbq_w, nc.sids$SID79) [1] 11.000000 10.500000 4.333333 7.750000 3.666667 6.400000 5.166667 5.200000 21.400000 9.666667 5.142857 10.142857 12.800000 7.666667 2.666667 14.333333 6.500000 9.833333 [19] 12.875000 1.333333 1.000000 2.000000 15.400000 10.000000 6.000000 9.333333 1.500000 1.000000 13.428571 8.000000 Tests d’autocorréla-on spa-ale •  Les diﬀérents tests mesurant l’autocorréla@on spa@ale : ü  Tests pour variables qualita-ves: test de couleur des cartes ü  Tests globaux d’autocorréla-on spa-ale : les valeurs sont-‐elles globalement (sur l’ensemble du site) posi-vement (néga-vement) corrélées? Deux indices principaux : I de Moran et C de Geary. I préféré à C, en raison d’une stabilité générale plus grande. ü  Tests locaux d’associa-on spa-ale (LISA): La valeur observée en i est-‐elle associée posi-vement aux valeurs des localisa-ons voisines (ressemblances) ou néga-vement (dissemblance)? Indices LISA : Décomposent l’indice global de façon à iden-ﬁer la contribu-on individuelle de chaque lieu. Perme_ent de détecter les poches locales d’autocorréla-on spa-ale Tests d’autocorréla-on spa-ale • 
Tests globaux d’autocorréla@on spa@ale : On veut tester : H0 : independance entre les valeurs de Y prises sur les diﬀérentes localisa-ons (spa-al randomness) contre H1 : associa-on ü  associa+on posi+ve : les valeurs de Y proches ont tendance à être regroupées dans l’espace (agréga-on) ü  associa+on néga+ve : les voisins tendent à être dissemblables (structure en damier) RQ : Les résultats de ces tests dépendent du choix de la matrice W 

variance
fornormal
each region
Regional x/y mean
valuesand
all come
from
distributions w/same
Regional x/y values all come from normal distributions w/same
mean and variance
for each
regionthe data on map and compute
Randomly
rearrange
I many
mean and
variance for each region
times, would
haveon
a normal
Randomly rearrange
the data
map anddistribution
compute I many
 Randomly rearrange the data on map and compute I many
times, would have
a normal
Why?
Because distribution
we use the normal distribution totimes,
calculate
the p-value
would
have a normal distribution

Why? Because we use the normal distribution to calculate the p-value
 Why? Because we use the normal distribution to calculate the p-value
Tests d’autocorréla-on spa-ale v  I de Moran (1950) Moran’sla Iplus classique d’associa-on spa-ale, analogue au Durbin-‐Watson de ü 
Sta-s-que Moran’s I
I the mean
Product Moran’s
of the deviation from
séries temporelles for all pairs of adjacent regions (w =1)
n
Product of the deviation from the mean
nfor all pairs of adjacent regions (wij=1)
ij
Product of the deviation from the mean
for all pairs of adjacent regions (wij=1)
n
n
n wij ( yi  y )( y j  y )
n
n
w
(
y
−
y)(
y
− y)
n wij ( yi i 1i yj )( y j ny )
n
ij
i
j wij ( yi  y )( y j  y )
I

i 1 i  j
i 1 i  j


I
 n
2  i=1 j=1
I

n
I
=


(
y

y
)
w


 ij

 i


 n
2
  ( yi  y)i21  wij i  j n

2 ( yi  y )   wij 


 i 1
 i  j 
S0 ( yi − y)
 i 1
 i  j 
n
n
∑∑
∑
Essentially a measure of
Sum of the weights (count
variance
regions
Essentially a measure
of across theSum
of all(count
adjacenti=1
pairs)
of the weights
Essentially a measure of
variance across the regions
of all adjacent pairs)
variance across the regions


Sum of the weights (count
n = number of regions
of all adjacent pairs)
• n = number
I a des ofvregions
aleurs comprises entre -‐1  ent =1number
of regions
 w = measure of spatial proximity between region i and
j
• wij = I measure
est dij’autant lus grand que d
es valeurs s
emblables tendent à se regrouper dans des of spatialpproximity
between
region
iwand
j
ij = measure of spatial proximity between region i and j

• 
• 
lieux proches (a_rac-on) I est d’autant plus pe-t que des valeurs dissemblables tendent à se regrouper dans des lieux proches (répulsion) Lorsque I est proche de 0, la répar--on des valeurs de Y se fait au hasard : il y a absence de corréla-on spa-ale. 18
18
AVen@on : le I de Moran dépend du choix de la matrice W et peut être aﬀecté par le niveau d’aggréga-on (eﬀet d’échelle) ainsi que par la forme des unités spa-ales 18
Tests d’autocorréla-on spa-ale ü  Loi de la sta-s-que de test sous H0 sous l’hypothèse de normalité (N) q  (N): on suppose que les observa-ons (y1,..yn) sont des réalisa-ons d’un vecteur gaussien à composantes i.i.d sous H0. avec Z=
I − E N (I )
→ N (0,1) sous H0
VN (I )
1
n −1
n 2 S1 − nS2 + 3S02
2
VN (I ) =
−
E
N (I )
2
2
S0 (n −1)
E N (I ) = −
Tests d’autocorréla-on spa-ale ü  Loi de la sta-s-que de test sous H0 sous Hypothèse de Randomisa9on (R) q  (R): Lorsqu’on répar-t plusieurs fois les données au hasard sur la carte (permuta-on des lieux) et qu’on calcule I à chaque fois, la répar--on des valeurs I obtenues est à peu près gaussienne. Z=
I − E R (I )
→ N (0,1) sous H0
VR (I )
1
n −1
n (n 2 − 3n + 3)S1 − nS2 + 3S02 − K (n 2 − n)S1 − 2nS2 + 6S02
E R (I ) = −
VN (I ) =
(
)
2
0
S (n −1)(n − 2)(n − 3)
n
K=
(
n∑ ( yi − y ) 4
i=1
# n
&
%% ∑ ( yi − y ) 2 ((
$ i=1
'
2
)−E
2
R
(I )
Tests d’autocorréla-on spa-ale ü  Loi de la sta-s-que de test sous H0 en l’absence d’hypothèse sur la distribu@on de I q  Si on a un doute sur N ou R, on simule la loi de I par Bootstrap pour es-mer la p-‐
value. Algorithme : 1)  On calcule I. 2)  On exécute B fois la tâche suivante : On permute aléatoirement les observa-ons et on recalcule I : on ob-ent un ensemble de B valeurs de l’indice de Moran (I1 ,..., I B )
3)  La p-‐value du test est le pourcentage de valeurs générées supérieures à I : Nombre de b / I b > I
B
Tests d’autocorréla-on spa-ale ü  Exemples Test de Moran avec hypothèse N : moran.test(nc.sids$SID79, listw=sids_nbq_w, alterna-ve="two.sided", randomisa-on=FALSE) Moran's I test under normality Pour un test à 5%, on reje_e H0 dès que p-‐value<5% H0: a
bsence data: nc.sids$SID79 d’autocorréla-on weights: sids_nbq_w H1: corréla-on Moran I sta-s-c standard deviate = 2.5927, p-‐value = 0.009524 alterna-ve hypothesis: two.sided sample es-mates: Ecart-‐type de la sta-s-que de Valeur de I Moran I sta-s-c Expecta-on Variance VarN(I) moran sous H0 lorsque (N) vraie 0.158977893 -‐0.010101010 0.004252954 EN(I) Conclusion : on reje_e H0 : il existe une autocorréla-on spa-ale globale, qui est posi-ve (signe de l’indice) Tests d’autocorréla-on spa-ale Test de Moran avec hypothèse R : moran.test(nc.sids$SID79, listw=sids_nbq_w, alterna-ve="two.sided") Moran's I test under randomisa-on data: nc.sids$SID79 weights: sids_nbq_w Moran I sta-s-c standard deviate = 2.6926, p-‐value = 0.00709 alterna-ve hypothesis: two.sided sample es-mates: Moran I sta-s-c Expecta-on Variance 0.158977893 -‐0.010101010 0.003943149 Conclusion : on reje_e H0 : il existe une autocorréla-on spa-ale globale, qui est posi-ve (signe de l’indice) Tests d’autocorréla-on spa-ale Test de Moran par bootstrap: bperm=moran.mc(nc.sids$SID79,listw=sids_nbq_w, nsim=999) Monte-‐Carlo simula-on of Moran's I data: nc.sids$SID79 B Rang de la valeur de I parmi weights: sids_nbq_w les valeurs générées number of simula-ons + 1: 1000 sta-s-c = 0.159, observed rank = 990, p-‐value = 0.01 alterna-ve hypothesis: greater H0 : absence d’autocorréla-on H1 : autocorréla-on posi-ve Tests d’autocorréla-on spa-ale Monte Carlo Moran’s I
• 
> mean(bperm$res[1:999])
[1] -0.0130007
> var(bperm$res[1:999])
[1] 0.004201103
> summary(bperm$res[1:999])
Min, Max, Median, Quartiles
> hist(bperm$res, freq=TRUE,
breaks=20, xlab="Simulated
Moran's I")
> abline(v=0, col=“red”)
Tests d’autocorréla-on spa-ale v  Indice de Geary C (1954) n
C=
n
(n −1)∑∑ wij ( yi − y j ) 2 )
i=1 i=1
n
2S0 ∑ ( yi − y) 2
i=1
ü 
ü 
ü 
ü 
C est compris entre O et 3 Des valeurs <1 indique une autocorrélation spatiale positive
Des valeurs >1 indique une autocorrélation spatiale négative
Des valeurs proches de 1 indiquent une absence d’autocorrélation positive
Tests d’autocorréla-on spa-ale ü  Loi de la sta-s-que sous H0 q  Sous les hypothèses (N) ou (R), q  Sinon, bootstrap Z=
C − E(C)
→ N (0,1) sous H0
V (C)
Tests d’autocorréla-on spa-ale Test de Geary C sous l’hypothèse R geary.test(nc.sids$SID79, listw=sids_nbq_w, alterna-ve="two.sided") Geary's C test under randomisa-on data: nc.sids$SID79 weights: sids_nbq_w Geary C sta-s-c standard deviate = 2.1501, p-‐value = 0.03155 alterna-ve hypothesis: two.sided sample es-mates: Geary C sta-s-c Expecta-on Variance 0.829864635 1.000000000 0.006261663 Conclusion : Au seuil 5% on reje_e H0: il existe une autocorréla-on spa-ale (posi-ve car C<1) Tests d’autocorréla-on spa-ale v  Choix de la sta@s@que de test ü  l’indice de Geary est plus sensible aux points aberrants ü  l’approxima-on gaussienne est meilleure pour I que pour G Tests d’autocorréla-on spa-ale v  Autocorrélogramme Une autocorréla-on spa-ale posi-ve signiﬁe qu’à une valeur élevée de yi correspondent des valeurs élevées dans les zones con-guës ou courte distance. Un moyen u-le de représenter les valeurs successives de l’indice d’autocorréla-on spa-ale (généralement I) en fonc-on de la distance ou de diﬀérents ordres de con-guïté est le corrélogramme. ü  L’Autocorrélogramme de Moran permet de visualiser la façon dont l’autocorréla-on change avec l’éloignement des lieux. La distance entre deux lieux est mesurée via la matrice de poids. • 
Pour chaque valeur h de distance (qui dépend du type de voisinage choisi), on détermine un W(h) qui vaut 1 si la distance vaut h, 0 sinon • 
on calcule un indice de moran I(h) appelé corrélogramme de Moran • 
On regarde alors comment évolue l'autocorréla-on avec l'ordre du voisinage. Cela permet d'appréhender des structures plus complexes, opérant a diﬀérentes échelles. Tests d’autocorréla-on spa-ale ü  CP : Distance de con-guité : Lorsqu'on mesure l'autocorréla-on, on mesure le lien entre la valeur observée en un point et les valeurs observées aux points voisins. Dans le graphe de con-guité, les voisins sont déﬁnis par deux points reliés par une arête. On parle de voisinage d'ordre 1. On peut aussi déﬁnir un voisinage d'ordre 2 qui consiste à prendre comme voisin d'un point, les voisins d'ordre 1 de ses voisins d'ordre 1. Ainsi de suite. On a alors distance h=con-guité d’ordre h Tests d’autocorréla-on spa-ale 0.2
0.1
0.0
-0.2
-0.1
Moran's I
• 
0.3
nc.sids$SID79
1
2
3
4
5
6
7
8
lags
Le proﬁl (ou la courbure) de l’autocorrélogramme est en général monotone décroissant. La courbure des indices correspondant à la vitesse à laquelle décroît l’indice. • 
La valeur maximale (à distance minimale), ici de 0,157, est le I de Moran (moran. Test) du premier voisinage. • 
La portée (ici voisinage à l’ordre 5) est le niveau au delà duquel l’indice s’annule ou change de signe . Tests d’autocorréla-on spa-ale ü  Rq : Graphe des voisinages de con-guité d’ordres 1 à h nb.lag=nblag(sids_nbq,2) (h=2) [[1]] Neighbour list object: Number of regions: 100 Number of nonzero links: 490 Percentage nonzero weights: 4.9 Average number of links: 4.9 [[2]] Neighbour list object: Number of regions: 100 Number of nonzero links: 868 Percentage nonzero weights: 8.68 Average number of links: 8.68 par(mfrow=c(2,1)) plot(nc.sids) ; plot(nb.lag[[1]], coords, add=T, col=2, main="h=1") plot(nc.sids) ;plot(nb.lag[[2]], coords, add=T, col=2, main="h=2") Tests d’autocorréla-on spa-ale cor=sp.correlogram(sids_nbq, nc.sids$SID79, order=8, method="I", style="W") Spa-al correlogram for nc.sids$SID79 method: Moran's I es-mate expecta-on variance standard deviate Pr(I) two sided 1 0.1589779 -‐0.0101010 0.0039431 2.6926 0.00709 ** 2 0.1048517 -‐0.0101010 0.0023014 2.3962 0.01657 * 3 0.0596812 -‐0.0101010 0.0017727 1.6574 0.09744 . 4 0.0127441 -‐0.0101010 0.0016063 0.5700 0.56868 5 0.0312701 -‐0.0101010 0.0016498 1.0185 0.30842 6 -‐0.0377728 -‐0.0101010 0.0018493 -‐0.6435 0.51991 7 -‐0.0759001 -‐0.0101010 0.0021352 -‐1.4240 0.15446 8 -‐0.1188161 -‐0.0101010 0.0024119 -‐2.2136 0.02685 * -‐-‐-‐ Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Tests d’autocorréla-on spa-ale nc.sids$SID79
0.0
-0.1
-0.2
Moran's I
0.1
0.2
0.3
plot(cor) •  La portée de l’autocorréla-on posi-ve est de 5 (voisinage d’ordre 5). On observe une légère auto-‐corréla-on néga-ve après lag 6. 1
2
3
4
5
lags
6
7
8
2
Tests d’autocorréla-on spa-ale • 
2/
Tests locaux d’autocorréla@on spa@ale (LISA) : Perme_ent de porter un regard plus précis sur les varia-ons sous-‐régionales, indépendamment de la structure spa-ale d’ensemble. Local Moran’s I
Local
II
Local Moran’s
Moran’s
 Used to determine
if local autocorrelation
exists around
ü  Permet, pour chaque région i, dregion
e déterminer si il existe ifulocal
ne autocorréla-on locale  Used
determine
autocorrelationexists
exists
around
each
I to

Used
to
determine
if
local
autocorrelation
around
Local
Moran’s
I
dans un voisinage de i. region
 Returns aneach
I value
for IIeach region
each
region
to sdetermine
if local
autocorrelation
exists around
ü  Souvent u-lisé après un test global pour voir i an II value
Used
an
each region
region
 Returns
Returns
value for
for each
Often used after global Moran’s I to see if:

region
I
–  le phénomène Y est homogène sur S (dans each
ce cas, les sta-s-ques locales ont des valeurs similaires sur  Often
used
after
global
Moran’s
tosimilar
seeif:if:across regions)
Often
used
after
global
Moran’s
II to
see
Returns
an I value for (local
each
region
 The study area
is homogeneous
statistics
chaque région), ou non.  The study
area
isisglobal
homogeneous
(local
statistics
similaracross
regions)
studyafter
area
homogeneous
(local
statistics
similar
regions)
 Often
used
to
see if:
There are
local
that
contribute
to
a Isignificant
globalacross
statistic
–  il existe des abéra-ons locales responsables de lThe
a outliers
signiﬁca-vité de la Moran’s
sta-s-que globale.  There
are
local
that contribute
contribute
toaasignificant
significant
globalstatistic
statistic
are area
localisoutliers
outliers
that
to
global
The
study
homogeneous
(local statistics
similarNeighborhood
across
regions)
–  Rq : peuvent être pra-qués même n
en dThere
’autocorréla-on globale ( yl’absence 
y
)
i
 There
that contribute to a significant Specified
global
statistic
Neighborhood
Neighborhood
nn((yyiiare
wyijlocal
)( youtliers
by the
Ii 
j  y)
2
Specified
I

w
y
)
Neighborhood
Specified
bybythethe
I

(
y

y
)
weights matrix
( yi i i y ) n( yj i  y ) 22
ijij
jj
v  Indices locaux de Moran: n
1
I = ∑ Ii
S0 i=1
ü  Rq : Ii =

Ii 
j
n


j wij ( y j  y)
jj((yyii  y)2 
 ( yThe ithy )point for
j
i
j
n( yi − y )∑ wij ( y j − y )
The ith point
pointfor
for
point
for
The ith
which
we
calculate
I
i
i=1
which we
calculateIIiIi i
we calculate
calculate
which
n
∑( y
j=1
j
in the
− y ) 2 The j pointsThe
jj points
points
The
points ininthe
the
local neighborhood
local
neighborhood
local neighborhood
neighborhood
Specified
by
the
weightsmatrix
matrix
weights
weights matrix
2/
v  Tests locaux de Moran ü  Pour chaque localisa-on i, on teste : H(i)0 : independance entre la valeur de Y dans la localisa-on i et celle de ses voisins contre H(i)1 : associa-on locale entre la valeur de Y prise par i et celles de ses voisins ü  La sta-s-que de test u-lisée est Ii. Elle est posi-ve (resp. néga-f) dans les régions ayant des valeurs similaires (resp. diﬀérentes) ü  La distribu-on de Ii sous H0 n’est pas Gaussienne et la randomisa-on ne marche pas. On u-lise l’approxima-on par « points de selles » (cf. Tiefelsdorf (2000; 2002)) pour approcher ce_e loi. Tests d’autocorréla-on spa-ale ﬁps=order(nc.sids$FIPSNO ) nclocI=localmoran(nc.sids$SID79, sids_nbq_w) pp=printCoefmat(data.frame(nclocI[ﬁps,], row.names=nc.sids$NAME[ﬁps]), check.names=FALSE) Ii
E.Ii
Var.Ii
Z.Ii
Pr.z...0.!
Alamance
0.07913669 -0.01010101 0.14471804 0.23457787
0.4073!
Alexander
-0.15454000 -0.01010101 0.22126296 -0.30706484
0.6206!
Alleghany
0.24506455 -0.01010101 0.29780788 0.46757790
0.3200!
Anson
0.03019857 -0.01010101 0.22126296 0.08567343
0.4659!
Ashe
0.44551026 -0.01010101 0.29780788 0.83488444
0.2019!
Avery
0.18605116 -0.01010101 0.17533601 0.46844387
0.3197!
Beaufort
0.15808868 -0.01010101 0.14471804 0.44211784
0.3292!
Bertie
0.12055810 -0.01010101 0.17533601 0.31203559
0.3775!
Bladen
-0.49749291 -0.01010101 0.17533601 -1.16397260
0.8778!
Brunswick
-0.03501441 -0.01010101 0.29780788 -0.04565254
0.5182!
Buncombe
-0.35214166 -0.01010101 0.12284806 -0.97587361
0.8354!
…………………………………….. Donne une indica-on sur La contribu-on de chaque région I à l’autocorréla-on globale I, et le sens de celle-‐ci. Une p-‐valeur <5% met en évidence une autocorréla-on locale signiﬁca-ve (répulsion si néga-ve, a_rac-on si posi-ve) RQ : sum(pp[,1])/100 0.1589779 (I de Moran) Tests d’autocorréla-on spa-ale v  Diagramme de Moran (Moran sca_erplot) : Permet à la fois une analyse globale et une analyse locale par la répar--on des points au sein de quatre quadrants. ü  Relie les valeurs de Y (axe des x) à celles de la variable décalée WY, représentant la valeur moyenne de Y sur les voisins (axe des y) Les trois dimensions du processus d’intégration régionale en Europe…
Figure 1 : Le diagramme de Moran
Décalage spatial standardisé
ü  Les quatre quadrants du diagramme correspondent aux 4 types d’associa-on locale pouvant exister entre les régions et leurs voisins LH
HH
0
Variable standardisée
LL
HL
²  Associa-ons posi-ves (Ii>0) : aggloméra-on locale de valeurs semblables -‐ HH (high-‐high): région à valeur de Y élevée entourée de régions à valeurs de Y similaires et élevées (ces régions sont appelées hot spots) -‐ LL (low-‐low) : région à valeur de Y faible entourée de régions à valeurs de Y similaires et faibles (cold spots) Les quadrants HH5 et LL sont associés à une autocorrélation spatiale positive car ils
indiquent un regroupement spatial de valeurs similaires. En revanche, les quadrants LH et HL
représentent une autocorrélation spatiale négative car ils indiquent un regroupement spatial de
valeurs dissemblables. Par conséquent, le diagramme de Moran peut être utilisé pour
visualiser les localisations atypiques, c'est-à-dire les régions qui se trouvent dans les quadrants
LH et HL.
Une des limites de cet outil est qu’il ne fournit pas d'information sur la significativité
des regroupements spatiaux. Celle-ci est obtenue à l'aide d'indicateurs locaux d’association
spatiale.
1.3 Les indicateurs locaux d’association spatiale (LISA)
Pour toutes nos analyses, nous utilisons des matrices standardisées en lignes, dès lors, la
moyenne des statistiques locales de Moran est égale à la statistique globale du I de Moran6.
Anselin [1995] définit un indicateur local d’association spatiale ou « LISA » comme toute
statistique satisfaisant deux critères. Premièrement, le LISA donne une indication sur le
regroupement spatial significatif de valeurs similaires autour de chaque observation.
Deuxièmement, la somme des LISA associés à toutes les observations est proportionnelle à un
indicateur global d’association spatiale.
²  Associa-ons néga-ves (Ii<0) : aggloméra-on locale de valeurs dissemblables -‐ HL (high-‐low) : région à valeur de Y élevée entourée de régions à valeurs de Y similaires et faibles (atypiques locaux) -‐ LH (low-‐high) : : région à valeur de Y faible entourée de régions à valeurs de Y similaires et élevées(atypiques locaux) Anselin fournit une version locale de la statistique
région i, elle s'écrit de la façon suivante :
Ii
5
6
xi
x
m0
wij ( x j
j
x ) avec m0
xi
x
2
, du I de Moran. Pour chaque
/N
i
A noter que « élevé » (resp. « faible ») signifie au-dessus (resp. en dessous) de la moyenne européenne.
La démonstration de cette propriété est établie par Anselin [1995].
-8-
[1]
Tests d’autocorréla-on spa-ale ü  Propriétés du diagramme de Moran ²  La pente de la droite de régression du nuage de points (regression de WY sur Y) correspond au I de Moran ²  La valeur moyenne de (Y, WY) se situe à l’intersec-on des axes (x,y) ü  U-lisa-on ²  Permet de détecter des points aberrants ²  La densité des quadrants permet voir le processus spa-al dominant ²  La non linéarité permet d’appréhender le degré d’homogénéité du processus spa-al : Un ajustement faible du nuage (R2 de la régression de WY sur Y faible) indique une inhomogénéité importante du phénomène Y (poches locales) et de sa dépendance spa-ale (non-‐sta-onarité) • 
Remarque : il est intéressant de normaliser Y avant de faire le graphique pour pouvoir ainsi comparer plusieurs diagrammes de Moran entre eux. Le graphique est alors centré en (0,0) Tests d’autocorréla-on spa-ale 10 15 20
Mecklenburg
Wake
0
Guilford
Cumberland
10
20
30
40
50
nc.sids$SID79
1.5
Commentaires : les Outliers contribuent à aﬀaiblir l’associa-on spa-ale posi-ve -‐ Hoke, Bladen, Lincon : faible taux de morts subites entourées de voisins ayant un fort taux -‐ Wake, Guilford, Cumberland : fort taux de morts subites entourées de voisins ayant un faible taux Hoke
Bladen
Lincoln
0.5
Gaston
Mecklenburg
Wake
Guilford
Cumberland
-0.5
spatially lagged as.vector(nc)
Gaston
5
nci=moran.plot (nc.sids$SID79, sids_nbq_w, labels=as.character (nc.sids$NAME)) nc=scale(nc.sids$SID79) ncis=moran.plot (as.vector(nc), sids_nbq_w, labels=as.character (nc.sids$NAME)) Hoke
Bladen
Lincoln
0
valeurs sta+s+quement signiﬁca+ves (p<0.05) spatially lagged nc.sids$SID79
ü  Exemple : ²  Tracer du diagramme de Moran avec surlignage des -1
0
1
2
as.vector(nc)
3
4
5
Tests d’autocorréla-on spa-al ²  Cartographie des valeurs du I de Morans locI=cut(nclocI[,1], breaks=c(min(nclocI[,1]), 0, max(nclocI[,1])), labels=c("nega-f", "posi-f"), include.lowest=TRUE) plot(nc.sids, col=c(3,5)[locI]) legend("topright", legend=c("I Nega-f", "I Posi-f"), ﬁll=c(3,5), bty="n", y.intersp=0.8) I Negatif
I Positif
Tests d’autocorréla-on spa-al ²  Cartographie des voisinages inﬂuants (points atypiques de la régression de wy sur y) inﬂ=apply(nci$is.inf, 1, any) (permet de repérer les points atypiques de la régression) y =nc.sids$SID79 lhy =cut(y, breaks=c(min(y), mean(y), max(y)), labels=c("L", "H"), include.lowest=TRUE) wy= lag(sids_nbq_w, nc.sids$SID79) lhwy=cut(wy, breaks=c(min(wy), mean(wy), max(wy)), labels=c("L", "H"), include.lowest=TRUE) lhlh = interac-on(lhy, lhwy, inﬂ, drop=TRUE) cols =rep(1, length(lhlh)) cols[lhlh == "H.L.TRUE"] = 2 cols[lhlh == "L.H.TRUE"] =3 cols[lhlh == "H.H.TRUE"] =4 plot(nc.sids, col=c(5,3,8,6)[cols]) legend("topright", legend=c("None", "HL", "LH", "HH"), ﬁll=c(5,3,8,6), bty="n", y.intersp=0.8) None
HL
LH
HH

Download Report