公告

[公告]
2014/01/17
由於已經是faculty的關係,不太有足夠時間寫部落格。因此更新的速度會相當緩慢。再加上近幾年來SAS GLOBAL FORUM沒有出現讓我覺得驚艷的技術文件,所以能分享的文章相對也減少許多。若有人推薦值得分享的SAS技術文件,請利用『問題討論區』告知。

2013/07/19
臉書留言板的功能因為有不明原因故障,因此特此移除。而intensedebate的留言板因管理不易,也一併移除。目前已經開啟內建的 G+ 留言系統,所以請有需要留言的朋友,可直接至『問題討論區』裡面留言。


2011年4月9日 星期六

Improving Your Statistical Consulting Prowess by Adding R to Your SAS® Repertoire

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

我不清楚有多少人和我一樣,用SAS分析資料,但是用R來做圖。其實我一開始堅持資料分析和做圖都在SAS內一次完成,但後來我終於明白SAS/GRAPH要做出等同用R畫出來的圖形品質,需要花費更多時間,於是我終於將做圖的工作完全交付給R,也展開了我要一直切換兩個軟體視窗的日子。但我一直相信可以在SAS裡面呼叫R來完成畫圖的工作,終於讓我在SAS Global Forum 2010上面發現到一篇技術文件,在此跟大家分享。

首先要讓SAS執行R的程式所需的第一件事情,是要去SAS官網下載SAS/IML Studio 3.2,這是一套龐大但免費的軟體,總容量是整整1GB,下載前需要去SAS官網註冊帳號,註冊好就可以直接下載了。網址是 http://support.sas.com/rnd/app/da/workshop.html

原文內雖然是說用SAS/IML Workshop 2.2,但我實際去官網查詢的結果,發現只有 SAS/IML Workshop 2.1的連結,而且沒有辦法用在最新的SAS 9.2版上面。所以請直接去下載SAS/IML Studio 3.2版即可。經過我的測試是可以執行以下的範例的。

下載完畢並安裝好後,就可以在SAS裡面寫R的程式並呼叫R來執行,而結果會顯示在SAS的視窗裡面。不過根據本人的經驗,事情沒有那麼簡單,因為當我裝好並執行一段官網上的範例時,視窗一直彈跳出軟體找不到R的錯誤訊息,經過與SAS官網客服人員的聯繫,才知道SAS/IML Studio 3.2目前只能與R的2.9.1和2.11.1這兩個版本相容。如果你已經安裝了最新的R 2.12.2版,目前唯一的解決方法就是移除最新版然後降級安裝成2.11.1版。如果你有某些特殊的函式比需要在R 2.12.2版才能執行,則必須再次安裝最新版的R才行。

所有的R程式,必須要用"submit /R;"開頭,用"endsubmit;"結尾。而以往在SAS裡面常用的注解符號"*"也必須換成井字符號"#"。在進行資料分析前,得先將SAS資料轉成R可以用的資料格式,可利用下面這段程式碼進行轉檔;
run ExportDataSetToR("libname.dataname", "R_dataname");
其中libname是寫SAS library的名稱,而dataname自然就是SAS資料的名稱。R_dataname則是填入資料集在R裡面所用的名稱。作者建議儘量少用底線當做用在R裡面的資料名稱,因為有可能會出問題。此外,如果是用自定的library,則必須照舊先在一開始用libname語法宣告自定library的路徑。底下的範例因為都是用SAS內部已經建好的Sashelp,所以沒有必要在一開始用libname來宣告。在此需特別注意這一點。以下用幾個範例來說明:

範例一;表格
這個範例是教如何在SAS裡面使用一些R的指令來整理表格,雖然我覺得SAS在這方面就比R強太多了,但這就是一個範例,程式如下;

run ExportDataSetToR("Sashelp.cars", "cars");
submit /R;
library(Hmisc, T); 
sink("c/sugi2010/car_summaries.txt");
car.sum <- summary.formula(DriveTrain ~ Type + Origin + MPG_City + MPG_Highway + Weight + Wheelbase + Length, data = cars, method = "reverse", overall = FALSE, test = TRUE,digits = 3, prmsd=TRUE, exclude1 = FALSE, long = TRUE);
print(car.sum, test = TRUE,digits = 3, prmsd=TRUE, exclude1 = FALSE, long = TRUE);
sink();
endsubmit;

為了顧及一些人不懂R的語法,所以在此稍微解釋一下,之後的範例會省略對R的介紹,有需要的人請自行參照R的線上手冊。

第一行就是作資料轉檔的語法,剛剛已經介紹過了。本例是拿Sashelp裡面的一個cars的檔案來做表格。"submit /R;"之後開始就全部都是R的語法了,裡面不會有任何SAS程式,直到最後用"endsubmit;"才是回到SAS語法。library() 是R裡面呼叫某個package的語法,Hmisc就是專門拿來做summary table的package。sink("....txt") 是先定義接下來所有的報表是要輸出到某txt文件裡面,記得結束要補上"sink();"在最後面才會完成存檔動作。summary.formula()就是製作的函數,裡面要放很多參數,在此就不一一介紹了,總之做完的結果是放在car.sum這個物件裡面。而print()是把car.sum這個物件印進剛剛用sink函數建立的純文字檔。而那個純文字檔打開後如下所示;

R繪製出來的圖形相當漂亮,但表格就很陽春,非常像古早PE2打出來的效果。因此這個時候就是SAS展現威能的時候了。用一個data step把那個純文字檔叫進來;
data cars2;
  infile "("c/sugi2010/car_summaries.txt"" dsd dlm="|"   
         lrecl=1024 firstobs=5 truncover;
  format junk $20. label $34. all $34. front $35. rear $35. test $35. ;
  input junk $ label $ all $ front $ rear $ test $ ; 
  if substr(junk,1,1)="+" then delete;
  drop junk;
run;
其中:
  • 用infile把純文字檔呼叫進SAS,用dlm="|"將每一行分段,用lrecl來定義整行最大長度,並用firstobs指定SAS從第五行開始讀。
  • 用format把剛剛用"|"將每一行切成數段,所以每一段就變成一個變數,因此需要定義變數的格式。在此將所有變數定義成文字格式,因為接下來不用再做分析了,所以只要能列印出及可。
  • 用input設定各變數的名稱,這些名稱要和format裡面所寫的變數名稱一樣。
  • 用junk變數把不必要的隔線刪除掉。

把不要的內容清乾淨後,用ods和proc report把結果印出來;
ods rtf file="c:\sugi2010\sas_r\car_summaries2_09oct20b.rtf" ;
title "Summaries by treatment group";
proc report data=cars2 nowindows;
  columns label all front rear test;
run; ods rtf close;
結果如下;

範例二;繪圖
在SAS裡面用R畫圖是我比較關心的,但和上面這個範例比起來並沒有太多不同。程式如下;
run ExportDataSetToR("Sashelp.CARS", "cars");
submit /R;
library(Hmisc, T);  
## Create data frame – cars2 – dropping cars with 3, missing, and 5 cylinders; 
cars2 <- cars[cars$Cylinders>3 & !is.na(cars$Cylinders) & cars$Cylinders != 5,] 
## place cars2 in search path 
attach(cars2); 
## Create summary data sets for Highway and city gas mileage MPG by number of 
cylinders and origin.  ; 
## Bootstrap resampling used to calculate 95% confidence intervals of mean; 
mean.hw <- summarize(MPG_Highway, llist(Cylinders, Origin), smean.cl.boot); 
mean.city <- summarize(MPG_City, llist(Cylinders, Origin), smean.cl.boot); 
## rename MPG_Highway and MPG_City to common name:  mean.mpg.; 
mean.hw2 <- upData(mean.hw, rename=list(MPG_Highway="mean.mpg")); 
mean.city2 <- upData(mean.city, rename=list(MPG_City="mean.mpg")); 
##create new variable type to label type of mean.mpg; 
mean.hw2$type[mean.hw2$mean.mpg>0]<- "Highway" 
mean.city2$type[mean.city2$mean.mpg>0]<- "City" 
## Vertical combine of the 2 data frames – new dataframe is named all; 
all <- rbind(mean.hw2, mean.city2);
##  Use Hmisc function xYplot to plot mean gas mileage and 95% confidence intervals; 
## by cylinders with separate panels for each of 3 origins.  ; 
## groups – separate curves by type within each panel.; 
## col = c(1,2) sets colors for 2 groups (1=black, 2=red). Lty sets line types; 
## layout=c(3,1)  determines number of rows and columns – c(3,1) is 3 columns and 1 row; 
xYplot(Cbind(mean.mpg, Lower, Upper) ~ Cylinders | Origin, data=all, groups=type, 
col=c(1,2), lty=c(1,2), layout=c(3,1), type='o'); 
endsubmit;

基本上你只要把前兩行和最後一行先在SAS裡面寫好,然後在第三行前面貼上R的程式即可。不過上面這段用R繪圖的語法,其實在SAS裡面是可以直接做的,下面是兩個軟體畫出來的圖形的比較;

左圖是用R畫出來的圖,右邊是用SAS畫的。乍看之下好像右圖比較好,但那是因為新版的SAS剛好有PROC SGPANEL這個程序把所有繪圖參數都最佳化了,如果遇到PROC SGPANEL沒有支援的圖形,那就不要太指望SAS能給你甚麼好圖了。而根據本人實際測試的結果,其實還是R做的圖比較好看,所以就大膽地在SAS/IML Studio底下執行R的程式吧!

CONTACT INFORMATION  
Name: Jeff Gossett
University of Arkansas for Medical Sciences. Department of Pediatrics.
One Children's Way
Little Rock, AR 72202
Work Phone: 501-364-6631
Fax: 501-364-1431
E-mail: GossettJeffreyM@uams.edu
Web: http://www.arpediatrics.org/research/biostatistics
CODE { display: block; /* fixes a strange ie margin bug */ font-family: Courier New; font-size: 8pt; overflow:auto; background: #f0f0f0 url(http://klcintw.images.googlepages.com/Code_BG.gif) left top repeat-y; border: 1px solid #ccc; padding: 10px 10px 10px 21px; max-height:200px; height:200px; // for IE6 line-height: 1.2em; }