[SAS] 在SAS/Base實作R


看起來跟iPad2好像是沒什麼差別……


話說,在上一篇的SAS文章裡頭我才在訐譙SAS這個坑錢巨獸,很多功能都分開來賣,比方說可以輕鬆地叫用R的矩陣語言環境SAS/IML就是其中一例(順便推薦一下The Do Loop這個部落格,作者是IML的開發人員,不過裡面有許多具啟發性的統計學文章)

但是!

生命會找到她的出路。在SAS 9之後有個SAS Base裡東西叫做Unnamed Pipe是讓SAS可以在自己的程序中叫用執行其他程式的工具,而R——這個其他程式——當然也是可以被執行的。

這意味著理論上我們一定可以在SAS裡叫用R,問題是要如何達成滑順、無痛的互動,例如R如何接SAS的資料、SAS又如何接回R的運算結果,R的code又應該怎麼管理……等等。官方開發的SAS/IML為此提供了一個無縫整合,但是它要$$,所以我們當作沒這回事。(喂)

就在我寫完上一篇文章之後沒多久,我在Journal of Statistical Software上命運般地邂逅了Xin Wei (2012)這篇論文,作者寫了一個MACRO實作了SAS Base-R的整合,當然立馬抓下來測試看看。

測試的結果,感覺非常棒!

這個作者取名為%PROC_R的巨集,基本概念就是運用Unnamed Pipe,並且讓我們可以直接在SAS文本下寫R的語言,將這些R的文本以一行一行字串資料的形式存成表單,想辦法把表單丟給R,叫用R執行,並且透過表單轉檔csv的機制達成運算結果可以互丟的歡樂整合。

廢話不多說,立刻分享我的簡單測試。
首先是SAS底下的一個情境——

*--Generate synthetic data set in SAS;
data raw;
array origin[4] origin1-origin4;
do i = 1 to 100;
do j = 1 to 4;
origin[j] = rand("Gaussian", 0, 1);
end;
output;
end;
drop i j;
run;

*--Conduct PCA;
proc princomp data = raw out = princomp std;
var _NUMERIC_;
run;

*--Compute the Euclidean distances for the standardized principal components;
data princomp;
set princomp (keep = prin:);
maha2 = uss(of prin:);
run;

我從一個標準常態裡抽了四個變數,每個一百筆,形成一筆測試資料。
這是一個四維度的資料空間,假設我想在SAS裡頭計算這一百筆資料的馬氏距離(Mahalanobis Distance),該怎麼做?事實是,SAS裡頭似乎沒有任何一個程序可以直接計算馬氏距離,即使在多變量分析裡它是那麼地基本。

SAS的DATA STEP並不適合拿來做這樣的操作,即便馬氏距離的公式並不複雜,以下取自WIKIpedia:

Formally, the Mahalanobis distance of a multivariate vector from a group of values with mean and covariance matrix is defined as:

SAS裡看似最有希望的PROC DISTANCE也沒有提供馬氏距離的選項。為了計算每筆資料向量(距離其平均向量的)的馬氏距離,必須採取迂迴策略:我們先對原始資料作一個主成分分析(Principal Component Analysis, PCA),取得標準化之後的特徵向量,把資料做轉換之後,再對這個轉換後的資料計算其傳統的歐式距離(Euclidean Distance),這個歐式距離,從數學上可以看得出來也就是原始資料的馬氏距離。上面的SAS程序就是在做這件事情(我算的是馬氏距離的平方)。

現在輪到我們的主角了!
可以達成SAS-R無痛整合的%PROC_R,其使用方式非常簡單,假設我們仍然是用SAS製作了測試資料,那麼我只要寫下面這個東西就可以得到馬氏距離:

%include "C:\SASdata\Proc_R.sas";
%Proc_R (SAS2R = raw, R2SAS = maha2);
datalines4;
raw <- as.matrix(raw)
invcov <- solve(cov(raw)) # Get the inverse of covraiance matrix
X <- sweep(raw, 2, apply(raw, 2, mean), '-') # Demean the data
maha2 <- rowSums((X %*% invcov) * X)
;;;;
%quit;

這個巨集有兩個引數,從名字上很容易明白其意圖。在叫用R的時候如果你想將SAS data set丟給R來處理,就把表單的名字丟給SAS2R;同理,在R裡頭有什麼結果想回傳成SAS data set的形式,就把其物件名稱丟給R2SAS。

紅色標記的部分,就是完全的R語言!

這個巨集丟給R的資料會是dataFrame,為了做矩陣操作會先將之轉換成Matrix物件。
第二行計算共變矩陣的反矩陣。
第三行對資料矩陣做減除平均數的動作。
第四行,完成馬氏距離的計算!我們將結果回丟給SAS,同時在R裡頭產生的log也會存入SAS形成一個字串表單。

MAGIC!

看起來我只是把SAS可以做到的事情丟給R來做,然後就高潮了。不過那是因為這個測試只涉及一個簡單的運算工作。現在設想一個更複雜也更實務的情境。如果我面對一個大型表單,想要偵測高維度下的離群值,我必須計算馬氏距離,但我想計算的是Robust Estimate,也就是不受離群值影響馬氏距離,在SAS裡你很快就卡住了,因為沒有矩陣環境,諸如Minimum Covariance Determinant這類的演算法都難以實現。雖然PROC ROBUSTREG裡面有一個離群值檢測的選項,可以對你的自變數向量做Robust馬氏距離的估計,而且它採用的(根據官方文件)正是MCD演算法!但是在面對大型資料時的效率非常差,我到目前為止還沒成功跑完過,但我又無法自己調校這個演算法,因為它寫死在PROC ROBUSTREG程序裡了,遑論這是一個以回歸估計為訴求的程序。

但是現在有了R,這一切都變得更有彈性!
要使用這個驚人的%PROC_R巨集,只要去下載作者的原始碼,然後記得把裡面指向R的路徑修正為你電腦裡的情況,並且在每次叫用時以%INCLUDE巨集指向這份原始碼,一切就OK了!

基本上是如此。
一些小問題例如定義R的工作空間與互丟資料時的暫存檔命名、存放位置等,都可以自己再微調,但這都不是重點。重點就是,IT'S MAGIC!

0 comment(s):

Post a Comment

回應文章前請注意下列三勿原則:

1)勿拍照;(→會有靈異的照片從你的相機裡跑出來...
2)勿餵食;(→會有飢渴的猛獸從我的網誌裡跑出來...
3)勿告白。(→會有奇怪的東西從站長體內裡跑出來...

謝謝大家的配合。
( > ー <)b