Spatial Autocorrelation Example

**************************************************************;
* AUTOCORR.SAS - Spatial Autocorrelation Analysis Program.   *;
* Compute Moran's I and Geary's C spatial autocorrelation    *;
* measures. Generate asymptotic tests as well as Monte Carlo *;
* tests for significance. Requires as input, both a Weight   *;
* data set and a Response Data Set.                          *;
* Barry Moser, Dept. Experimental Statistics, Louisiana      *;
* State University, Baton Rouge, LA 70803-5606.              *;
* VOICE: 504-388-8303   E-Mail: barry@insight.agadm.lsu.edu  *;
**************************************************************;

TITLE1 'SPATIAL AUTOCORRELATION';

/*
 * Here we have 16 items with a response X measured on each.
 * These items are linked together according to the weight
 * matrix in data WEIGHTS. Use Moran's I and Geary's C statistics
 * to test for spatial autocorrelation among the items.
 */

DATA RESPONSE;
**************************************************************;
* I is the location index and X is the response value.       *;
**************************************************************;
 INPUT I X;  /* X'S ARE ASSUMED IN THE CORRECT ORDER */
 LIST;
CARDS;
 1  4
 2  4
 3  4
 4  4
 5  4
 6 -2
 7 -1
 8  5
 9  4
10 -1
11 -2
12  5
13  5
14  5
15  5
16  5
;

DATA WEIGHTS;
**************************************************************;
* I and J are the location indices for the pair and WT is    *;
* the weight assigned to the pair. Specify all pairwise      *;
* links between units considered to be directly correlated.  *;
**************************************************************;
 INPUT I J WT;
 LIST;
CARDS;
 1  6 1
 2  6 1
 2  7 1
 3  6 1
 3  7 1
 4  7 1
 5  6 1
 5 10 1
 6  9 1
 7  8 1
 7 12 1
 8 11 1
 9 10 1
10 13 1
10 14 1
10 15 1
11 14 1
11 15 1
11 16 1
11 12 1
;

PROC IML; /* SAS INTERACTIVE MATRIX LANGUAGE */
**************************************************************;
* Use IML to do the actual computations. This lets us easily *;
* do the Monte Carlo test as well. For the Monte Carlo test  *;
* you will be given the list of results for each simulation  *;
* as well as the rank of each result. The final value will   *;
* be the observed value and rank for the data. If the rank   *;
* is extreme then you will reject Ho that the data are       *;
* distributed as under the simulation (randomly).            *;
* NOTE: The printing format for the weights might need to be *;
* changed if you are using non-integer valued weights. I used*;
* 3.0 to keep the amount of the output down.                 *;
**************************************************************;

START LOADDATA;
 USE RESPONSE VAR {X};
 READ ALL;
 CLOSE RESPONSE;
 N=NROW(X);
 W=J(N,N,0); /* INITIALIZE WEIGHT MATRIX TO ZEROS */
 USE WEIGHTS VAR {I J WT};
 READ ALL;
 /* PLACE WEIGHTS INTO THE WEIGHT MATRIX */
 DO K=1 TO NROW(I);
  W[I[K],J[K]]=WT[K];
  W[J[K],I[K]]=WT[K];  /* ASSUMING SYMMETRY IN WEIGHTS */
 END;
 /* DROP THE I J AND WT MATRICES */
 FREE I J K WT;
 PRINT / 'INITIAL DATA ' ,, 'RESPONSE',X,, 'WEIGHTS', W[FORMAT=3.0];
   /* $$$    THE FORMAT MAY NEED TO BE CHANGED IF    $$$ */
   /* $$$     FRACTIONAL WEIGHTS ARE TO BE USED.     $$$ */
FINISH;

START STATS;
 XBAR=X[+]/N;
 Z=X-REPEAT(XBAR,N,1);
 ZSQ=Z`*Z;
 XVAR=ZSQ/(N-1);
 VMRATIO=XVAR/XBAR;
 Z2=Z#Z; Z4=Z2`*Z2;
 B2=N*Z4/ZSQ**2; /* KURTOSIS */
 FREE Z2 Z4;
 SUMW=0; S1=0;
 DO K=1 TO N;
  DO J=1 TO N;
   IF K <> J THEN
    DO;
      SUMW=SUMW + W[K,J];
      S1=S1 + (W[K,J]+W[J,K])**2;
    END;
  END;
 END;
 S1=S1/2;
 S2=W[,+]+W[+,]`;
 S2=S2`*S2;
 SUMW2=SUMW*SUMW;
 TEMP=N||XBAR||XVAR||VMRATIO||B2;
 RNAME={" "};
 CNAME={"   N" "   XBAR" "VARIANCE" "V/M RATIO" "KURTOSIS"};
 PRINT / 'SUMMARY STATISTICS ' TEMP[ROWNAME=RNAME COLNAME=CNAME];
 FREE XBAR XVAR VMRATIO K J RNAME CNAME TEMP;
FINISH;

START AUTOCORR;
 I=0; C=0;
 DO K=1 TO N;
  DO J=1 TO N;
    IF K <> J THEN
     DO;
       I = I + W[K,J] * Z[K]*Z[J];
       C = C + W[K,J] * (X[K]-X[J])**2;
     END;
  END;
 END;
 I=(N/SUMW)*(I/ZSQ);
 C=((N-1)/(2*SUMW))*(C/ZSQ);
 FREE K J;
FINISH;

START AUTOSTAT;
 RNAME={" "}; CNAME={"I   " "C   "};
 TEMP=I||C;
 PRINT ,, 'SPATIAL AUTOCORRELATION STATISTICS'
      TEMP[ROWNAME=RNAME COLNAME=CNAME];
 FREE RNAME CNAME TEMP;
FINISH; /* AUTOSTAT */

START TESTOFI;
 EXPECT=1/(1.0-N);
 RNAME={" "}; CNAME={"  E(I)  " "OBSERVED"};
 TEMP=EXPECT||I;
 PRINT "TEST OF I", "EXPECTED AND OBSERVED VALUES"
       TEMP[ROWNAME=RNAME COLNAME=CNAME];
 VAR=((N**2*S1-N*S2+3*SUMW2)/(SUMW2*N**2-1))-EXPECT**2;
 ZZ=(I-EXPECT)/SQRT(VAR);
 P=1-PROBNORM(ABS(ZZ));
 CNAME={"VARIANCE" "   Z   " "  P>|Z|  "};
 TEMP=VAR||ZZ||P;
 PRINT  'TEST BASED ON  NORMALITY:'
         TEMP[ROWNAME=RNAME COLNAME=CNAME];
 VAR=(N*((N**2-3*N+3)*S1-N*S2+3*SUMW2)-
       B2*((N**2-N)*S1-2*N*S2+6*SUMW2))/
       ((N-1)*(N-2)*(N-3)*SUMW2) - EXPECT**2;
 ZZ=(I-EXPECT)/SQRT(VAR);
 P=1-PROBNORM(ABS(ZZ));
 TEMP=VAR||ZZ||P;
 PRINT  'TEST BASED ON  RANDOMIZATION:'
         TEMP[ROWNAME=RNAME COLNAME=CNAME];
 FREE EXPECT VAR ZZ P RNAME CNAME TEMP;
FINISH;

START TESTOFC;
 EXPECT=1.0;
 RNAME={" "}; CNAME={"  E(C)  " "OBSERVED"};
 TEMP=EXPECT||C;
 PRINT "TEST OF C", "EXPECTED AND OBSERVED VALUES"
       TEMP[ROWNAME=RNAME COLNAME=CNAME];
 VAR=((2*S1+S2)*(N-1)-4*SUMW2)/(2*(N+1)*SUMW2);
 ZZ=(C-EXPECT)/SQRT(VAR);
 P=1-PROBNORM(ABS(ZZ));
 CNAME={"VARIANCE" "   Z   " "  P>|Z|  "};
 TEMP=VAR||ZZ||P;
 PRINT  'TEST BASED ON  NORMALITY:'
         TEMP[ROWNAME=RNAME COLNAME=CNAME];
 /* ZZ AND P BELOW ARE TEMPORARIES */
 ZZ=(N-1)*S1*(N**2-3*N+3-(N-1)*B2)-(N-1)/4*S2;
 P=(N**2+3*N-6-(N**2-N+2)*B2)+SUMW2*(N**2-3-(N-1)**2*B2);
 VARC=ZZ*P/(N*(N-2)*(N-3)*SUMW2);
 ZZ=(C-EXPECT)/SQRT(VAR);
 P=1-PROBNORM(ABS(ZZ));
 TEMP=VAR||ZZ||P;
 PRINT  'TEST BASED ON  RANDOMIZATION:'
         TEMP[ROWNAME=RNAME COLNAME=CNAME];
 FREE EXPECT VAR ZZ P RNAME CNAME TEMP;
FINISH;

START MONTE;
 ALLI=I; ALLC=C;
 DO REP=1 TO 99;
    R=UNIFORM(REPEAT(0,N,1));
    Y=X; X[RANK(R),]=Y;
    Y=Z; Z[RANK(R),]=Y;
    RUN AUTOCORR;
    ALLI=I//ALLI; ALLC=C//ALLC;
 END;
 PRINT / 'MONTE CARLO TEST - OBSERVED VALUE IS LAST VALUE';
 R=RANK(ALLI); ALLI=ALLI||R;
 CNAME={"  I " "RANK"};
 PRINT 'RANKS FOR I VALUES:' ALLI[COLNAME=CNAME];
 PRINT / 'MONTE CARLO TEST - OBSERVED VALUE IS LAST VALUE';
 R=RANK(ALLC); ALLC=ALLC||R;
 CNAME={"  C " "RANK"};
 PRINT 'RANKS FOR C VALUES:' ALLC[COLNAME=CNAME];
 FREE ALLI ALLC REP R Y ;
FINISH;

START SAVEW;
DO I=1 TO NROW(W);
 W[I,I]=2;
END;
CREATE WMATRIX FROM W;
APPEND FROM W;
CLOSE WMATRIX;
FINISH;

RESET NOLOG NONAME FW=10;
RUN LOADDATA;             /* READ THE DATA INTO THE MATRICES */
RUN STATS;                /* COMPUTE SUMMARY STATISTICS      */
RUN AUTOCORR;             /* COMPUTE AUTOCORRELATION STATS.  */
RUN AUTOSTAT;             /* PRINT COMPUTED STATS.           */
RUN TESTOFI;              /* COMPUTE TESTS OF I              */
RUN TESTOFC;              /* COMPUTE TESTS OF C              */
RUN MONTE;                /* PERFORM MONTE CARLO SIMULATIONS */

RUN SAVEW;                /* SAVE WEIGHT MATRIX FOR MDS      */

QUIT;                     /* DONE WITH SAS/IML               */


/*
 * Use the weights as similarities among the units. Use MDS
 * to display units as linked together by the weights.
 */

DATA WMatrix;
 Set WMatrix;
 Id+1;
RUN;

Proc MDS Data=WMatrix SIMILAR=2 LEVEL=Ordinal DIM=2 PCONFIG OUT=Config;
 ID Id;
 Var Col1-Col16;
Run;

%let plotitop = gopts = device = gif, Color = Black, Colors = Black;

title2 'Plot of Weight Matrix Configuration';
%Plotit(Data=Config(Where=(_type_='CONFIG')), Datatype=mds,
           LabelVar=Id, VtoH=1.75);


/*
 * Add lines in graph to joined units.
 */
Data New;
 Set Weights;
 PType+1;
 Id=I; Output;
 Id=J; Output;
 Keep PType Id;
Run;

Proc Sort Data=New;
 By Id;
Run;

Proc Sort Data=Config;
 By Id;
Run;

Data New;
 Merge New Config(Where=(_type_='CONFIG'));
 By Id;
Run;

Data New;
 Set Config(In=In1) New;
 If In1 Then PType=0;
Run;

GOptions Reset=Symbol Reset=Axis;

Title2 "Weight Configuration With Connections";
Proc GPlot Data=New;
 Plot Dim2*Dim1=PType / NoLegend VAxis=Axis1 HAxis=Axis2;
 Axis1 Order=(-2 To 2 By 1) Length=4in Label=(A=90);
 Axis2 Order=(-2 To 2 by 1) Length=4in;
 Symbol1 C=Black V=Dot I=None;
 Symbol2 C=Black V=None I=Join L=1 R=1000;
Run;
Quit;

 

SPATIAL AUTOCORRELATION

INITIAL DATA

 

RESPONSE

 

4
4
4
4
4
-2
-1
5
4
-1
-2
5
5
5
5
5

 

WEIGHTS

 

0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	1	1	0	0	0	0	0	0	0	0	0
0	0	0	0	0	1	1	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	1	0	0	0	0	0	0	0	0	0
0	0	0	0	0	1	0	0	0	1	0	0	0	0	0	0
1	1	1	0	1	0	0	0	1	0	0	0	0	0	0	0
0	1	1	1	0	0	0	1	0	0	0	1	0	0	0	0
0	0	0	0	0	0	1	0	0	0	1	0	0	0	0	0
0	0	0	0	0	1	0	0	0	1	0	0	0	0	0	0
0	0	0	0	1	0	0	0	1	0	0	0	1	1	1	0
0	0	0	0	0	0	0	1	0	0	0	1	0	1	1	1
0	0	0	0	0	0	1	0	0	0	1	0	0	0	0	0
0	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	1	1	0	0	0	0	0
0	0	0	0	0	0	0	0	0	1	1	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0

 

---

 

SPATIAL AUTOCORRELATION

		N	XBAR	VARIANCE	V/M RATIO	KURTOSIS
SUMMARY STATISTICS		16	3	7.46666667	2.48888889	2.37755102

 

		I	C
SPATIAL AUTOCORRELATION STATISTICS		-0.9642857	2.44419643

 

TEST OF I

 

		E(I)	OBSERVED
EXPECTED AND OBSERVED VALUES		-0.0666667	-0.9642857

 

		VARIANCE	Z	P>|Z|
TEST BASED ON NORMALITY:		0.0360244	-4.7292651	1.12667E-6

 

		VARIANCE	Z	P>|Z|
TEST BASED ON RANDOMIZATION:		0.0368859	-4.6737112	1.47903E-6

 

TEST OF C

 

		E(C)	OBSERVED
EXPECTED AND OBSERVED VALUES		1	2.44419643

 

		VARIANCE	Z	P>|Z|
TEST BASED ON NORMALITY:		0.07647059	5.22250725	8.82583E-8

 

		VARIANCE	Z	P>|Z|
TEST BASED ON RANDOMIZATION:		0.07647059	5.22250725	8.82583E-8

 

---

 

SPATIAL AUTOCORRELATION

MONTE CARLO TEST - OBSERVED VALUE IS LAST VALUE

 

	I	RANK
RANKS FOR I VALUES:	-0.2142857	26
	-0.3285714	13
	-0.2928571	16
	0.45	100
	-0.1714286	35
	-0.2	28
	-0.1428571	39
	0.16428571	89
	-0.2428571	21
	0.23571429	96
	-0.1071429	45
	-0.1071429	46
	-0.0071429	68
	-0.2	30
	0.32857143	98
	-0.0142857	67
	-0.1	48
	-0.3785714	8
	-0.1285714	43
	-0.3642857	9
	-0.0571429	57
	0.11428571	81
	-0.2214286	25
	-0.1071429	47
	0.15	84
	-0.5785714	3
	0.35	99
	-0.1857143	32
	0.29285714	97
	0.20714286	93
	0.15	85
	-0.0357143	61
	-0.15	37
	0.15714286	88
	-0.2928571	17
	-0.15	38
	0.10714286	80
	0.2	92
	-0.2071429	27
	-0.0214286	64
	-0.1	49
	-0.0214286	65
	-0.0714286	54
	-0.5	4
	-0.2642857	19
	0.02142857	76
	-0.3	15
	-0.2571429	20
	0.06428571	78
	0.20714286	95
	0.14285714	83
	-0.1428571	40
	0.00714286	74
	0.20714286	94
	-0.2357143	23
	-0.0071429	71
	0.11428571	82
	-0.05	60
	-0.2785714	18
	-0.3571429	11
	-0.0571429	59
	-0.2285714	24
	-0.0571429	58
	0.15714286	87
	0.15714286	86
	-0.7214286	2
	-0.2357143	22
	-0.1214286	44
	-0.0214286	63
	0.06428571	77
	-0.3142857	14
	-0.0714286	53
	-0.1357143	41
	-0.0071429	70
	0.00714286	73
	-0.0571429	56
	-0.0928571	50
	-0.0071429	69
	-0.0571429	55
	-0.3571429	10
	-0.4928571	5
	0.18571429	91
	-0.1285714	42
	-0.2	29
	-0.3857143	7
	0.18571429	90
	-0.0857143	52
	-0.3428571	12
	-0.0857143	51
	-0.4142857	6
	-0.0142857	66
	0	72
	-0.1857143	31
	0.09285714	79
	0.02142857	75
	-0.0214286	62
	-0.15	36
	-0.1785714	34
	-0.1785714	33
	-0.9642857	1

 

---

 

SPATIAL AUTOCORRELATION

MONTE CARLO TEST - OBSERVED VALUE IS LAST VALUE

 

	C	RANK
RANKS FOR C VALUES:	1.04799107	66
	1.27566964	84
	1.171875	75
	0.56584821	1
	1.31919643	85
	1.00446429	57
	0.890625	42
	0.88392857	37
	1.32589286	86
	0.63616071	3
	1.04799107	67
	0.84709821	28
	0.87388393	32
	1.20535714	77
	0.59933036	2
	0.6796875	8
	1.00111607	56
	1.68415179	96
	0.88727679	40
	1.60044643	95
	1.01116071	59
	0.85044643	30
	1.02455357	61
	0.84709821	29
	0.87723214	34
	1.72098214	97
	0.63950893	4
	0.99107143	54
	0.84375	26
	0.88392857	38
	0.796875	20
	1.04129464	65
	1.20870536	78
	0.6796875	9
	1.2421875	80
	1.12834821	73
	1.02790179	62
	0.76004464	14
	0.97098214	52
	0.67633929	7
	0.87053571	31
	0.87723214	35
	0.79352679	17
	1.51674107	94
	1.11495536	71
	0.87723214	36
	1.48995536	93
	1.24888393	82
	0.91741071	49
	0.68303571	11
	0.79352679	18
	0.91071429	44
	1.03125	63
	0.79352679	19
	1.24888393	83
	0.67299107	5
	0.890625	43
	1.15513393	74
	1.40959821	91
	1.2421875	81
	0.75	12
	1.05133929	68
	0.97098214	53
	0.6796875	10
	0.77008929	15
	1.90513393	99
	1.43973214	92
	1.01116071	60
	0.67633929	6
	0.75669643	13
	1.23214286	79
	0.84375	27
	0.88392857	39
	0.92410714	50
	1.11160714	70
	1.07142857	69
	0.83370536	23
	0.79352679	16
	0.9609375	51
	1.34263393	87
	1.80133929	98
	0.84375	25
	0.88727679	41
	1.36607143	88
	1.38950893	89
	0.83370536	22
	1.0078125	58
	1.18861607	76
	1.03794643	64
	1.39620536	90
	0.80022321	21
	1.12834821	72
	1.00111607	55
	0.84040179	24
	0.91741071	48
	0.87723214	33
	0.91741071	47
	0.9140625	46
	0.9140625	45
	2.44419643	100

 

---

 

SPATIAL AUTOCORRELATION

Multidimensional Scaling: Data=WORK.WMATRIX

Shape=TRIANGLE Condition=MATRIX Level=ORDINAL

Coef=IDENTITY Dimension=2 Formula=1 Fit=1

Mconverge=0.01 Gconverge=0.01 Maxiter=100 Over=2 Ridge=0.0001

Iteration	Type	Badness- of-Fit Criterion	Change in Criterion	Convergence Measures
				Monotone	Gradient
0	Initial	0.3892	.	.	.
1	Monotone	0.4052	-0.0160	0.0740	0.5697
2	Gau-New	0.3278	0.0774	.	.
3	Monotone	0.3380	-0.0102	0.0736	0.3556
4	Gau-New	0.3264	0.0116	.	.
5	Monotone	0.3322	-0.005751	0.0279	0.3490
6	Gau-New	0.3186	0.0136	.	.
7	Monotone	0.3231	-0.004572	0.003299	0.3475
8	Gau-New	0.3016	0.0216	.	0.0923
9	Gau-New	0.2997	0.001890	.	0.0436
10	Gau-New	0.2993	0.000418	.	0.0208
11	Gau-New	0.2992	0.0000975	.	0.0111
12	Gau-New	0.2992	0.0000296	.	0.007106

 

Convergence criteria are satisfied.

 

Configuration
	Dim1	Dim2
1	1.31	1.10
2	0.86	0.02
3	1.56	-0.05
4	1.26	-1.16
5	-0.29	1.38
6	0.84	0.87
7	0.83	-0.87
8	0.29	-1.38
9	0.33	1.40
10	-0.83	0.87
11	-0.84	-0.87
12	-0.33	-1.40
13	-1.26	1.16
14	-1.56	0.05
15	-0.86	-0.02
16	-1.31	-1.10

 

---

 

 

---