我不清楚有多少人和我一樣，用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裡面寫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");範例一；表格
這個範例是教如何在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
沒有留言:
張貼留言
要問問題的人請在文章下方的intensedebate欄位留言，請勿使用blogger預設的意見表單。今後用blogger意見表單留言的人我就不回應了。