Everything in Its Place: Efficient Geostatistical Analysis with SAS/STAT ® Spatial Procedures

原文載點：http://support.sas.com/resources/papers/proceedings10/337-2010.pdf

SAS 在處理空間統計和空間模型上一向相當薄弱，根據這幾年全球最大的統計會議 Joint Statistical Meeting 所顯示，越來越多人在做空間統計分析（我自己也是）。所以特別利用這一篇文章來介紹一份 2010 年發表的 SAS 技術文件，來介紹在 SAS 裡面到底開怎麼樣進行空間統計分析。

SAS 軟體裡面提供了三個程序，可拿來處理空間統計分析：
(1) PROC VARIOGRAM: 計算資料在空間上的相關程度
(2) PROC KRIG2D: 空間資料差補
(3) PROC SIM2D: 空間資料模擬

除了上述這些程序外，SAS 的 ODS/GRAPH 也提供一些已經最佳化的圖形讓使用這能夠輕鬆輸出。因此，本篇將使用一個有關全球太陽輻射度的資料來示範該如何使用這些程序：

data daGSIdata;
input Easting Northing daGSI @@;
label daGSI='Daily Average Global Solar Irradiance';
datalines;
13.6 33.4 17.4906 116.6 40.9 17.9709 86.0 215.3 17.0232
248.1 266.9 18.8413 40.2 149.3 17.6422 249.9 219.5 18.4011
298.6 19.8 18.3800 161.6 262.3 18.6537 165.9 204.9 18.8376
236.3 249.9 18.2338 141.5 108.1 18.3977 253.0 156.4 18.8378
218.6 81.2 18.6493 46.5 77.8 17.7590 292.8 47.7 18.5304
83.0 164.0 17.7159 110.8 225.6 17.9139 202.9 189.1 18.7072
199.5 17.5 18.8718 103.7 158.7 17.9594 184.7 252.5 18.6570
76.0 200.3 16.9481 51.5 35.6 17.7746 158.1 220.3 18.4775
230.9 55.0 18.9772 158.5 226.4 18.3128 252.8 131.5 18.6633
61.2 20.6 17.8255 279.6 212.6 19.0338 294.1 42.3 18.6408
80.9 82.3 17.6355 124.8 253.5 18.3031 18.0 286.4 17.7586
229.7 154.7 18.5950 71.8 47.0 17.9481 280.9 210.9 19.1626
92.9 125.8 17.6257 28.6 275.5 17.5359 186.4 111.7 18.5684
16.1 167.7 17.8476 114.9 130.2 17.7519 140.7 272.4 18.1968
122.5 155.3 17.7771 63.8 259.7 17.4883 262.1 12.3 18.5385
257.6 262.2 18.9482 128.8 229.1 17.9638 269.8 15.9 18.4247
288.6 236.9 18.6137 124.5 130.2 17.9114 7.3 117.2 17.3567
153.0 299.9 18.0288 75.8 85.2 17.5313 67.9 247.3 17.6824
81.7 199.0 17.2191 275.3 260.2 18.9168 180.6 39.6 18.6395
28.8 10.3 17.6430 217.2 76.1 18.7481 282.4 95.2 18.0833
143.0 196.7 18.5825 42.0 245.1 17.2660 179.1 187.1 18.6990
41.3 134.6 17.5906 217.0 130.8 18.6883 166.1 152.1 18.5032
186.4 59.2 18.9922 264.9 294.8 18.3716 228.6 42.1 18.7740
214.1 38.5 19.0145 123.5 66.2 18.5217 26.3 20.4 17.4753
39.7 78.8 17.6835 113.8 247.1 18.0232 225.0 63.3 18.8826
242.4 29.8 19.0110 260.5 158.1 18.7582 297.8 270.0 19.0336
111.8 90.1 18.2008 79.6 43.7 17.7467 102.8 90.7 18.0879
51.1 184.0 17.4158 216.3 85.0 18.7803 11.7 261.0 17.5346
291.6 265.7 19.2310 119.8 125.1 17.7093 127.4 176.8 18.1153
152.7 74.3 18.5789 202.1 70.3 18.7343 67.7 218.1 17.0247
138.5 38.3 17.8701 72.9 163.9 17.6970 189.7 144.7 18.4584
144.3 121.6 18.3561 120.7 246.9 18.0552 216.6 22.4 19.3410
;


這筆資料總共上百個觀測值，每個觀測值有三個變數，分別是座標（Easting & Northing）和太陽輻射度（daGSI）而這些觀測值大約散佈在一個300公里X300公里的平面上，所以首先我們想要了解如何在一張地圖上呈現出這些觀測值。此外，任意兩個點之間的距離也是我們想要知道的數據。因此，這兩個要求，可以用 PROC VARIOGRAM 程序輕鬆完成。程式如下所示：
ods graphics on;
proc variogram data=daGSIdata plots=equate;
    compute novariogram nhclass=30;
    coord xc=Easting yc=Northing;
    var daGSI;
run;
ods graphics off;
程式解說：
(1) plots=equate 表示做出來的地圖要有等距的X軸和Y軸。
(2) compute novariogram nhclass=30 表示去做兩兩觀測值距離的直方圖，而這些距離將會被分成30等份。而分好的距離被定義為延遲距離(lag distance)，從 0 到 30。數字越大距離越長。
(3) coord xc=Easting yc=Northing 表示設定 X 軸是 Easting 變數，Y 軸是 Northing 變數。
(4) var daGSI 表示觀測值為 daGSI 變數。

因此，利用 ODS GRAPHICS 的威能，SAS 能立刻生出下面兩個圖形：

其中上圖就是觀測值在平面上的散布圖，每個點的深淺顏色所代表的太陽輻射度可以參考右邊的漸層棒。下圖則是把所有任意兩點之間的距離算完收集好後製作成的直方圖，從右上角的圖例可知兩個軸的最大值以及任意兩點的最大距離。特別注意的是，X 軸已經改成 30 等份的距離，所以每一個等份相當於 13.69 公里。因此，資料本身在地圖上的敘述統計量大概就這樣完成了。

接下來所要關心的，就是任意兩點間的相關程度。最簡單的方法，就是去計算 empirical semivariogram。計算公式如下：

其中 Si, Sj 是兩點位置，而 h 就是兩點的距離。Z(S) 則是定義在一個二維平面裡面的隨機域(random field)，通常假設為二階平穩(second-order stationary)的隨機域。程式如下：
ods graphics on;
proc variogram data=daGSIdata plots(only)=semivar;
    compute lagdistance=15 maxlags=12 cl;
    coord xc=Easting yc=Northing;
    var daGSI;
run;
ods graphics off;
程式解說：
(1) plots(only)=semivar 表示只要畫 semivariogram 的圖。
(2) compute lagdistance=15 maxlags=12 cl 表示定義每一個延遲距離為 15 公里，但圖形最多呈現到第 12 個延遲距離。最後的 cl 表示要計算信賴區間。

圖形如下：

注意一點就是，雖然程式裡面設定只秀出第十二個延遲距離，但圖上卻有十三條線，那是因為最左邊的那個是 lag = 0。剛剛有提到 empirical semivariogram 是由二階平穩的隨機域所計算出來的，而若能滿足二階平穩隨機域的假設，則 semivariogram 的圖形會有兩個特徵，一是圖形會隨著距離的增加單調遞增到某一點，而通常從原點到這一個點的距離就稱之為 range。二是經過這個點之後，semivariogram 的走勢會趨緩成類似一個平行線上下震盪走動，而這趨緩的程度則稱之為 scale 或 sill。如果資料裡面含有某種空間趨勢的話，則 sill 就不會發生，因此 semivariogram 就會如同上圖一樣一直往上飄去。

所以，當要用這筆資料配飾任何統計模型時，我們都要設法移除掉空間上的相關性，這個道理就跟在做長期趨勢模型的時候要移除掉時間上的相關性是一樣的。比方說，以 PROC GLM 程序為例。程式如下：
proc glm data=daGSIdata plots=none;
     store out=trendStore / label='Trend Analysis Information';
     model daGSI = Easting Easting*Easting Northing Northing*Northing;
     output out=resdaGSIdata predicted=pred residual=resdaGSI;
run;
此例所用來移除空間上的相關性的方法是把兩個地理變數各自加上一個二次項，如紅字所示。當然這不是最好的方法，不過就當成一個範例來看。配飾好後把殘差項輸出成 resdaGSI 這個新資料檔。為了檢驗這個殘差項已經沒有空間相關性的干擾，最直接的方法就是把這個檔案丟回去剛剛那個計算 semivariogram 的程式去看：
ods graphics on;
proc variogram data=resdaGSIdata;
     compute lagdistance=15 maxlags=12 cl;
     coord xc=Easting yc=Northing;
     var resdaGSI;
run;
ods graphics off;
圖形如下所示：

明顯得可以看到大概在第六條線(lag=5)之後，semivariogram 就變成下降的趨勢了，然後就開始上下擺盪，因此可以說用地理變數的二次項可以移除掉一點資料的空間相關性。

除此之外，殘差項還有兩個現象也需要觀察，首先是 semivariogram 的 sill 會不會因為方向的不同而改變，或者是在同一個角度但不同方向的 range 會不會改變。若 sill 不會變但 range 會變，這種現象稱為 anisotropy，估狗翻譯為"各向異性"。程式如下：
ods graphics on;
proc variogram data=resdaGSIdata plots(only)=semivar;
     compute lagdistance=15 maxlags=12 cl ndirections=4;
     coord xc=Easting yc=Northing;
     var resdaGSI;
run;
ods graphics off;
此程式碼跟上一個程式碼差別只有 compute statement 後面多了個 ndirections=4。這個指令會把 180 度切成四等份，因此四個角度便依序是 0 度、45 度、90 度、135 度。圖形如下所示：

上述的 semivariogram 圖形都只能算是描述性的統計量，裡面並沒有包含任何參數和模型。如何用統計模型去配飾 semivariogram 並且能產生平滑的曲線，則需要在 PROC VARIOGRAM 程序裡面額外呼叫 model statement 來進行。該成續提供了多種不同的半變異模型(semivariation model)，使用者可以同時在一個 PROC VARIOGRAM 程序裡面進行模型估計，SAS 還能幫忙把配飾出來的 semivariogram 曲線畫成一張圖來方便比較。程式如下：
ods graphics on;
proc variogram data=resdaGSIdata plots(only)=(fit);
     compute lagdistance=15 maxlags=12 cl;
     coord xc=Easting yc=Northing;
     var resdaGSI;
     model form=auto(mlist=(gau,exp,sph,she,mat) nest=1 to 2) cl;
     store out=daGSIsemiv / label='Global Solar Irradiance Semivariances';
run;
ods graphics off;
從紅字部分可以發現，在 PROC VARIOGRAM 程序底下的 model 並不用特地把變數名稱寫上，因為變數名稱已經定義在 var statement 裡面了，所以 model statement 只需要把想要配飾的半變異模型的名稱寫上就可以了。在這個範例裡面，我們一次配飾五個模型型態，分別是高斯(gau)、指數(exp)、球面(sph)、正弦孔效應(she，全名sin-hole effect)以及Matérn(mat，這是一個人名)。更威的地方是，ODS GRAPHICS會自動將這五個模型做出來的 semivariogram 線畫在同一張圖上，如下所示：

回顧本文最一開始的那張散佈圖，我們可以隱隱約約發現，全球太陽輻射度在那個平面上有呈現某種趨勢，但因為只有這麼幾個點有觀測值，所以直觀上給人感覺也不是非常強烈。如果能找出一個方法來準確的預測空白處的全球太陽輻射度的話，整張平面就會出現不同的色彩，也提供了一個比較容易判斷空間分布的結果。要完成這種目標，最簡單的方法就是用 kriging（克利金）法來做空間差補。在 SAS 裡面，有提供  PROC KRIGE2D 程序來完成空間的差補與作圖。程式如下：
proc krige2d data=resdaGSIdata outest=resdaGSIpred plots=equate;
     restore in=daGSIsemiv / info(det);
     coordinates xc=Easting yc=Northing;
     predict var=resdaGSI;
     model storeselect;
     grid x=0 to 300 by 7.5 y=0 to 300 by 7.5;
run;
特別注意的是，PROC KRIGE2D 所使用的資料不是原始的 daGSIdata，而是從之前那個 PROC GLM 另存的殘差項 resdaGSIdata。然後要把上一個程式另存出來的 daGSIsemiv 叫回來，以提供之後要用哪種模型內差的依據。coordinates 設定 X軸和 Y軸變數，而要預測的變數 resdaGSI 就要宣告在 predict 語法裡面。model storeselect 則是可以選擇要用什麼模型去預測。我們在上一例已經用了五種不同的模型，相關數據也都存在 daGISsemiv，如果沒有特別指定，SAS 就會去把最好的模型拿來使用。最後把差補範圍和差補距離宣告在 grid 語法就可以了。送出程式後 SAS 就會自動生出下面這張圖：

以本例來看，SAS 使取 SHE-SPH 模型去做克利金。如果想要特別指定某種模型，就直接在 model storeselect寫上模型名稱即可，如下所示：
proc krige2d data=resdaGSIdata outest=resdaGSIpredAlt plots=equate;
     restore in=daGSIsemiv;
     coordinates xc=Easting yc=Northing;
     predict var=resdaGSI;
     model storeselect(model=gau);
     grid x=0 to 300 by 7.5 y=0 to 300 by 7.5;
run;
圖形：

上述都是在講殘差像的地圖作法，但有時候我們真正關心的則是原本DAGSI的預測值，所以也可以用之前的方法來繪圖。原文內還有提到另一種方法，但因為結果很類似，所以不多加累述。DAGSI的圖如下所示：

CONTACT INFORMATION
Your comments and questions are valued and encouraged. Contact the author:
Alexander Kolovos
SAS Institute Inc.
SAS Campus Drive
Cary, NC 27513
919-531-2165
alexander.kolovos@sas.com
(PS) 此人是我朋友的朋友，現已從SAS離職，目前在加州聖地牙哥自己開公司中。可用下面這個email和他聯繫：alexander.kolovos@spacetimeworks.com