Shellfish Spatial Autocorrelation
Moran's I and Geary's C Statistics

dm "output;clear;log;clear";
**************************************************************;
* ShellfishAutocorr.sas -- Compute Moran's I and Geary's C   *;
* statistics for the pipis shellfish data reported in Manly  *;
*                                                            *;
* Data: Manly, B.F.J. 2001. Statistics for Environmental     *;
* Science and Management. Chapman & Hall, Boca Raton, p.224. *;
*                                                            *;
* Barry Moser, Dept. Experimental Statistics, Louisiana      *;
* State University, Baton Rouge, LA 70803-5606.              *;
* Date: 02/20/2001                                           *;
**************************************************************;
Options PS=55 LS=80 PageNo=1 NoDate FullSTimer
        FORMCHAR='|----|+|---+=|-/\<>*';
GOptions Device=gif Transparency NoBorder NoPrompt
         VSize=6 in HSize=6 in
         HText=1 FText=Swiss HTitle=1 FTitle=Swiss;


Title1 "Shellfish Spatial Autocorrelation";

Data Pipis;
 Input LowWaterDistance @;
 Do BeachDistance=0 To 200 By 20;
  Input Count @;
  Quadrat+1;
  Output;
 End;
 Label LowWaterDistance="Distance From Low Water (m)"
       BeachDistance="Distance Along Beach (m)"
       Count="Counts of Pipis (Paphies australis)"
       Quadrat="Sampling Quadrat ID";
Datalines;
 0   1   0   4   0   0   0   3   0   2   0   0
10   0   0   0   0 104   0   0   0   1   0   0
20   7  24   0   0 240   0   0 103   1   0   0
30  20   0   0   0   0   0   3 250   7   0   0
40  20   0   2   4   0 222   0 174   4   0  58
50   0   0  11   0   0 126   0  62   7   6  29
60   0   0   7   0   0   0   0   0  23   7  29
70   0   0   0   0  89   0   0   7   8   0  30
;

Proc Print Data=Pipis;
 Id Quadrat;
Run;

Title2 "Map of the Data";
Proc Plot Data=Pipis;
 Plot LowWaterDistance*BeachDistance=" " $ Count / Box
      VRef=5 15 25 35 45 55 65
      HRef=10 30 50 70 90 110 130 150 170 190;
Run;
Quit;

Title2 "Summary Statistics";
Proc Means Data=Pipis N Mean Var Min Max;
 Var Count;
Run;

Title2 "Negative Binomial Model of Counts";
Proc Genmod Data=Pipis;
 Model Count = / Dist=NegBin Link=Log;
 Estimate "Population Mean mu" Intercept 1 / Exp;
Run;

/*
 * Generate the distance information and output as 
 * Quadrat i, Quadrat j, and Distance ij. There will
 * be n(n-1)/2 distances.
 *
 * This data set can then be used to assign weights
 * for the autocorrelation program as seen below.
 */

Title2 "Inter-quadrat Distances";
Proc IML;
 Start Dist;
 Use Pipis;
 Read All Var{LowWaterDistance BeachDistance} Into X;
 Read All Var{Quadrat} Into Q;
 Close Pipis;
 n=NRow(X);
 D=J(1,3,0);
 Names={"Qi" "Qj" "Dij"};
 Create Distance From D[Colname=Names];
 Do i=1 To n-1;
  D[1]=Q[i]; 
  Do j=i+1 To n;
   D[2]=Q[j];
   Diff=X[i,]-X[j,];
   D[3]=Sqrt(Diff*Diff`); /* Compute Euclidean Distance From Centers */
   Append From D;
  End;
 End;
 Close Distance;
 Finish Dist;
 Run Dist;
Quit;

Title3 "First 40 Observations";
Proc Print Data=Distance(Obs=40);
Run;

Title3 "Distribution of Distances";
Proc Freq Data=Distance;
 Table Dij;
Run;

/*
 * Construct weight information. Weights could be
 * based on joins or other function of distance.
 * You might try (Dij=20) to give weight=1 to units
 * that are exactly 20m apart, or try (1/(Dij/10)**2)
 * to base weights on inverse of squared distance.
 * Examples:
 * Wt=(Dij=20);
 * Wt=(Dij<=20);
 * Wt=1/((Dij/20)**2);
 */

Data Weights;
 Set Distance;
 i=Qi; j=Qj;
 Wt=1/((Dij/20)**2);
Run;

Title2 "Moran's I and Geary's C Statistics";
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: If the amount of data (number of plots) is large,    *;
* computational time can be long.                            *;
**************************************************************;

START LOADDATA;
 USE Pipis VAR {Count};
 READ ALL Into X;
 CLOSE Pipis;
 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               */

       Plot of LowWaterDistance*BeachDistance$Count.  Symbol used is ' '.       
                                                                                
     --+------+------+------+------+------+------+------+------+------+------+--
  70 + 0   |  0   |  0   |  0   |  89  |  0   |  0   |  7   |  8   |  0   |  30+
     |     |      |      |      |      |      |      |      |      |      |    |
     |     |      |      |      |      |      |      |      |      |      |    |
     |-----+------+------+------+------+------+------+------+------+------+----|
     |     |      |      |      |      |      |      |      |      |      |    |
     |     |      |      |      |      |      |      |      |      |      |    |
  60 + 0   |  0   |  7   |  0   |  0   |  0   |  0   |  0   |  23  |  7   |  29+
     |     |      |      |      |      |      |      |      |      |      |    |
D    |     |      |      |      |      |      |      |      |      |      |    |
i    |-----+------+------+------+------+------+------+------+------+------+----|
s    |     |      |      |      |      |      |      |      |      |      |    |
t    |     |      |      |      |      |      |      |      |      |      |    |
a 50 + 0   |  0   |  11  |  0   |  0   | 126  |  0   |  62  |  7   |  6   |  29+
n    |     |      |      |      |      |      |      |      |      |      |    |
c    |     |      |      |      |      |      |      |      |      |      |    |
e    |-----+------+------+------+------+------+------+------+------+------+----|
     |     |      |      |      |      |      |      |      |      |      |    |
F    |     |      |      |      |      |      |      |      |      |      |    |
r 40 + 20  |  0   |  2   |  4   |  0   | 222  |  0   | 174  |  4   |  0   |  58+
o    |     |      |      |      |      |      |      |      |      |      |    |
m    |     |      |      |      |      |      |      |      |      |      |    |
     |-----+------+------+------+------+------+------+------+------+------+----|
L    |     |      |      |      |      |      |      |      |      |      |    |
o    |     |      |      |      |      |      |      |      |      |      |    |
w 30 + 20  |  0   |  0   |  0   |  0   |  0   |  3   | 250  |  7   |  0   |  0 +
     |     |      |      |      |      |      |      |      |      |      |    |
W    |     |      |      |      |      |      |      |      |      |      |    |
a    |-----+------+------+------+------+------+------+------+------+------+----|
t    |     |      |      |      |      |      |      |      |      |      |    |
e    |     |      |      |      |      |      |      |      |      |      |    |
r 20 + 7   |  24  |  0   |  0   | 240  |  0   |  0   | 103  |  1   |  0   |  0 +
     |     |      |      |      |      |      |      |      |      |      |    |
(    |     |      |      |      |      |      |      |      |      |      |    |
m    |-----+------+------+------+------+------+------+------+------+------+----|
)    |     |      |      |      |      |      |      |      |      |      |    |
     |     |      |      |      |      |      |      |      |      |      |    |
  10 + 0   |  0   |  0   |  0   | 104  |  0   |  0   |  0   |  1   |  0   |  0 +
     |     |      |      |      |      |      |      |      |      |      |    |
     |     |      |      |      |      |      |      |      |      |      |    |
     |-----+------+------+------+------+------+------+------+------+------+----|
     |     |      |      |      |      |      |      |      |      |      |    |
     |     |      |      |      |      |      |      |      |      |      |    |
   0 + 1   |  0   |  4   |  0   |  0   |  0   |  3   |  0   |  2   |  0   |  0 +
     --+------+------+------+------+------+------+------+------+------+------+--
       0     20     40     60     80     100    120    140    160    180    200 
                                                                                
                              Distance Along Beach (m)