MACRO MORANTEST : Test d'autocorrélation spatiale de Moran et de Geary sous SAS

**Haache** · 28/07/2016, 21h10

Bonjour,

Je vous propose un nouvel élément à utiliser : MACRO MORANTEST : Test d'autocorrélation spatiale de Moran et de Geary sous SAS

Bonjour chers programmeurs.

Je vous propose dans ce blog une macro que je nomme morantest. Cette macro effectue les deux tests d'autocorrélation spatiale les plus connus, notamment le test de Moran (1950) et celui de Geary (1954). Une théorie sur ces tests est disponible ici ou ici.

La macro prend trois paramètres : Data où on renseigne la table des données, Var où on indique le variable dont l'autocorrélation spatiale sera testée et Poids où on indique la matrice de poids (une matrice carrée qui mesure la proximité entre les points, ça peut être une matrice de contiguïté ou de distance). Voici le code de la macro :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
 
 
%macro MoranTest(Data= ,Var= , Poids=)/STORE Des="Tests pour détecter une autocorrélation spatiale";options nonotes;
 
ods exclude all;
 
proc iml;
 
Title "Tests d'autocorrélation spatiale de MORAN et de GEARY";
 
/******************************************* Lecture des données dans une matrice******************************************************************/
 
use &Data; read all var {&Var} into Y;
 
/*****************************************Matrice de poids dans le même ordre que les Id**********************************************************/
 
use &Poids; read all into W;
 
/*Standardiser les lignes*/
 
W_row=W[,+];
 
do j=1 to nrow(W_row);
 
if W_row[j,1]=0 then W_row[j,1]=1;
 
end;
 
W=W/W_row;
 
 
 
 
 
N=nrow(Y);
 
Wcol=W[+,];
 
Wrow=W[,+];
 
K=(sum((Y-mean(Y))##4)/N)/((sum((Y-mean(Y))##2))/N)**2;
 
/******************************************* Calcul de la statistique*******************************************************************/
 
MORAN=(t(Y-mean(Y)) * W * (Y-mean(Y)))/(t((Y-mean(Y))) * (Y-mean(Y))); 
 
EMORAN=-1/(N-1);   /*Espérance de moran*/
 
S0=sum(W);
 
*S1=0.5*sum((W+t(W))##2);
 
S1=sum(W#W+W#t(W));
 
S2=sum((Wcol+t(Wrow))##2);
 
/*VMORAN=(N**2*S1-N*S2+3*S0**2)/((N-1)*(N+1)*S0**2)-(1/(N-1))**2; /*Variance de moran*/
 
VMORAN=(((N*((N**2-3*N+3)*S1-3*S2+3*S0**2))-(K*((N**2-3)*S1-2*N*S2+6*S0**2))))/((N-1)*(N-2)*(N-3)*S0**2)-EMORAN**2;
 
SMORAN=(VMORAN)**0.5;
 
RESULTAT=MORAN||EMORAN||VMORAN;
 
NAME={"Statistique","Espérance","Variance"};
 
OUT= nrow(Y);
 
OUTCN={"Nombre d'observations"};
 
OUT1={&VAR};
 
options notes;
 
ods select all;
 
print OUT1[label=" " colname={"Variable"}] OUT[label=" " colname=OUTCN] ;
 
print RESULTAT[label="Indice de Moran" colname=NAME  rowname=" "];
 
options nonotes;
 
ods exclude all;
 
MORANST=(MORAN-EMORAN)/SMORAN;
 
PROB1=CDF("normal",MORANST,0,1);   /*unilatérale à gauche*/
 
PROB2=1-CDF("normal",MORANST,0,1);   /*unilatérale à droite*/
 
PROB3=2-2*CDF("normal",abs(MORANST),0,1);   /*bilatérale*/
 
RESULTAT=MORANST || PROB1 || PROB2 ||PROB3;
 
NAME={"Statistique standardisée", "Prob unilatérale à gauche H1: I0", "Prob bilatéral: Alternative générale"};
 
options notes;
 
ods select all;
 
print RESULTAT[label="Test de Moran" colname=NAME rowname=" " format=10.5];
 
options nonotes;
 
ods exclude all;
 
/********************************************************************GEARY INDEX*****************************************************/
 
Y2=Y##2;
 
Wcol=W[+,];
 
GEARY=(N-1)*(sum(Y2)-2*t(Y)*W*Y+Wcol*Y2)/(2*S0*sum((Y-mean(Y))##2));
 
EGEARY=1;
 
/*VGEARY=( (2*S1+S2) * (N-1) -4*S2 )/( 2*(N+1)*S0 );*/
 
/*VGEARY=S0/(2*S1+S2);*/
 
VGEARY=(((N-1)*S1*(N**2-3*N+3-(N-1)*K))-(0.25*((N-1)*S2*(N**2+3*N-6-(N**2-N+2)*K)))+(S0**2*(N**2-3-(N-1)**2*K)))/(N*(N-2)*(N-3)*S0**2);
 
SGEARY=VGEARY**0.5;
 
RESULTAT=GEARY||EGEARY||VGEARY;
 
NAME={"Statistique","Espérance","Variance"};
 
options notes;
 
ods select all;
 
print RESULTAT[label="Indice de Geary" colname=NAME  rowname=""];
 
options nonotes;
 
ods exclude all;
 
GEARYST=(GEARY-EGEARY)/SGEARY;
 
PROB1=CDF("normal",GEARYST,0,1);   /*unilatérale à gauche*/
 
PROB2=1-CDF("normal",GEARYST,0,1);   /*unilatérale à droite*/
 
PROB3=2-2*CDF("normal",abs(GEARYST),0,1);   /*bilatérale*/
 
RESULTAT=GEARYST || PROB1 || PROB2 ||PROB3;
 
NAME={"Statistique standardisée", "Prob unilatérale à gauche H1: G1", "Prob bilatéral: Alternative générale"};
 
options notes;
 
ods select all;
 
print RESULTAT[label="Test de Geary" colname=NAME rowname=" " format=10.5];
 
quit;
 
title " ";
 
%put ***************************************************************************************;
 
%put *              MERCI D AVOIR UTILISE LES MACROS "REGRESSIONS SPATIALES"               *;                                                      
 
%put *                          DE Elysée Aristide HOUNDETOUNGAN                           *;  
 
%put *                          Courriel :  ariel92and@gmail.com                           *; 
 
%put *                  N hésitez pas à envoyer un mail en cas de besoin                   *; 
 
%put ***************************************************************************************;
 
 
 
 
 
%mend;

Comme exercice, supposons que nous disposons de la base de données ci-contre : Ici

Si vous décompressez ce fichier, vous y trouverez deux bases SAS nommées Base et matrice. Base contient les données sur les pays de l'Afrique décrits par 5 variables : une variable pays qui indique le libellé du pays, l'identifiant des pays, le code des pays, la température moyenne et les taxes. Les données sur la température et sur les taxes ne sont pas des données réelles.

Les sont simulées de sorte que la température soit spatialement autocorrélée et les taxes non autocorrélées. Nous disposons également du shapefile des pays d'Afrique dans la librairie Maps de SAS (ceci vient par défaut avec toutes les versions de SAS. La table s'appelle Maps.Africa). Décidons d'abord de représenter les variables sur des cartes

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
pattern6 v=s c=green;pattern5 v=s c=cx55aa11;
 
pattern4 v=s c=cx229922;
 
pattern3 v=s c=yellow;
 
pattern2 v=s c=cxffaa00;
 
pattern1 v=s c=cxff5555;legend1 shape=bar(.15in,.15in) frame cshadow=graylabel=(position=middle height=3.5pct font="arial/bold")value=(height=3.5pct font="arial");proc gmap data=Base map=Maps.Africa;format Precipitation Taxe comma10.1;id ID;choro Temperature / coutline=cx888888 legend=legend1; choro Taxe/ coutline=cx888888 legend=legend1; run;

On soupçonne déjà une autocorrélation spatiale avec la température. Les valeurs les plus élevées sont au sud de l'Afrique. En Afrique de l'Ouest les valeurs sont plus faibles. Par ailleurs, les taxes ne révèlent pas une telle caractéristiques (mélange de valeurs; il n'y a pas a priori un ordre sur la carte). Pour en être sûr, nous allons utiliser notre fameuse macro pour tester l'autocorrélation spatiale. A cette étape nous aurons besoin de la seconde base (Matrice qui est la matrice de contiguïté) des pays.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
 
 
%Morantest(Data=Base, Var=Temperature, Poids=Matrice);%Morantest(Data=Base, Var=Taxe, Poids=Matrice);

Pour la variable Taxe, ni les tests (unilatéral à gauche, à droite et bilatéral) de Moran et de Geary ne sont significatifs. Mais pour la variable température, le test unilatéral à droite de Moran est significatif. Autrement dit, la température est à une autocorrélation spatiale positive (I>0). Ce résultat est également confirmé par le test de Geary (G

Qu'en pensez-vous ?