2008年4月20日

繪製ROC曲線和比較兩條ROC曲線下面積

SAS實用巨集簡介:繪製ROC曲線和比較兩條ROC曲線下面積

目的為介紹SAS官方網頁中關於ROC曲線的巨集(Macro)程式,程式分別為「繪製ROC曲線」和「以無母數方式比較兩條ROC曲線下面積」。本文中所提及之程式,僅為原作者的部分語法,詳細語法請見原文說明,兩篇原文出處連結如下:

Sample 25017: Nonparametric comparison of areas under correlated ROC curves (http://support.sas.com/kb/25/017.html)
 roc.sas程式:http://support.sas.com/kb/25/addl/fusion_25017_6_roc.sas.txt

Sample 25018: Plot ROC curve with labelled points for a binary-response model (http://support.sas.com/kb/25/018.html)
 rocplot.sas程式:http://support.sas.com/kb/25/addl/fusion25018_4_rocplot.sas.txt


請先將程式下載至本機硬碟中,建議將Sample 25017的程式取名為roc.sas,而Sample 25018的程式取名為rocplot.sas。必須要注意的是,如果是想比較兩條ROC曲線下面積,也就是roc.sas這個巨集程式,必須安裝SAS/IML。

繪製ROC曲線

在SAS說明檔中其實已經有介紹如何繪製出簡易的ROC曲線,但在我們諮詢的個案中比較常被問到的問題為,圖上的某個點其「原始數值為何」,利用rocplot.sas 便可輕鬆完成這個工作。
假設有一組資料,研究者想以年齡和某一個BW指標來判斷罹患某疾病的機會。
問題:繪製以年齡作為判斷此是否罹患此疾病的ROC曲線。

第一步:先將資料讀入SAS中
data Age_data;
input disease age bw@@;
datalines;
0 50 65 0 39 61 0 21 70 0 61 67
0 30 55 0 35 63 0 25 72 0 41 66
0 43 52 0 36 54 0 37 76 0 25 61
0 41 53 0 62 55 0 28 70 0 33 68
1 52 45 1 49 61 1 47 42 1 62 31
1 55 67 1 70 61 1 75 55 1 77 52
1 81 58 1 64 53 1 62 41 1 39 57
1 61 51 1 61 55 1 57 49 1 79 47
;
run;

第二步:將rocplot.sas程式預先讀入SAS中
%include "c:\temp\ROC\rocplot.sas";
/*
%include “程式所放的位置”;
*/
透過%include這個巨集指令,可將rocplot.sas呼叫至SAS程式中。

第三步:執行proc logistic程式,並將部分結果輸出,以配合rocplot.sas程式
proc logistic data=age_data;
model disease (event='1')=age / outroc=roc1 roceps=0;
output out=outp p=phat;
run;
這個程式中,我們將disease變項中的「1」視為事件(罹患疾病),並將繪製ROC曲線所需的資料輸出為「roc1」,和每個原值(年齡)的預測機率資料輸出為「outp」並將預測機率的變項命名為「phat」。基本上,使用時只需要改動依變項(disease)和自變項(age)這兩個部分即可。

第四步:使用%rocplot這個巨集程式。
%rocplot (outroc = roc1, out = outp, p = phat, id = age);
這個巨集中,outroc= 資料檔,配合proc logistic的outroc輸出。
out= 資料檔,也同樣配合proc logistic 的out輸出。
p= 預測機率值變項名稱,在輸出預測機率時的命名。
Id= 變項名稱,要依據那個變項繪製ROC曲線上所標記的原始數值,通常應與自變項相同。

第五步:圖形輸出結果



比較兩條ROC曲線下面積

問題:比較以年齡或BW指標對於判斷是否罹患疾病,此兩指標ROC曲線下面積。

第一步:先將roc.sas程式呼叫入SAS中
%include "c:\temp\ROC\roc.sas";


第二步:執行自變項為age的logistic 模型
proc logistic data=age_data;
model disease(event='1')=age / outroc=roc1_age roceps=0;
output out=outp_age p=phat_age;
run;
為了區別age和bw這兩個結果,這裡我們將age的output輸出的檔案命名為「outp_age」而預測機率命名為「phat_age」。

第三步:執行自變項為BW的logistic 模型
proc logistic data=age_data;
model disease(event='1')=bw / outroc=roc1_bw roceps=0;
output out=outp_bw p=phat_bw;
run;


第四步:執行%roc 巨集
%roc (data = outp_age outp_bw,
var = phat_age phat_bw,
response = disease);
在data=資料檔,這後面主要放前面步驟透過output所輸出的檔案,因此目前所放的為「outp_age」和「outp_bw」兩個檔案。
var=變項名稱,在output檔中,預測機率的變項名稱。
response= 依變項名稱,宣告依變項名稱為何。

第五步:結果
於「ROC Curve Areas and 95% Confidence Intervals」的表格中,age的AUC(area under curve, 曲線下面積)為0.912,而bw的AUC為0.6836。age和bw這兩個自變項AUC比較,可於「Contrast Test Results」中看到p-value為0.0378,因此這兩條ROC曲線下面積為顯著的不同。

2008年2月18日

Yahoo 的 e-mail?

我小玩了Yahoo的知識網一段時間,心裡想說有一些建議想給「知識網」,就遇到我的第一個問題,找了老半天為何沒有聯絡方式,好歹有建議表可以讓我送出吧?

沒想到...

我真的找不到,只好回到Yahoo首頁,很興奮的看到了「關於Yahoo!奇摩」裡面有一個「聯絡我們」,沒想到卻是Yahoo公司台北的電話和地址....整個傻眼!

回頭努力找,還是找不到,心理一橫,我也來找找Google的聯絡方式。好吧!我也是找到Google的地址,但很快的,我在Gmail中就看到意見表了。

回頭找找Yahoo mail,在一個很奇怪的角落「Yahoo!奇摩隱私權保護政策」最下面「看法與建議」終於看到了意見表!很快的打完我的幾意見,當然包括找不到和Yyahoo互動的意見,興高采烈的按下「送出」,咦~~沒反應耶,這時我心理的OS就是....『不會吧,Yahoo不支援Firefox 2...』,只好很無奈的打開IE 7,把剛剛寫好的建議貼過去,才順利寄出。

整個過程我大約花了2個小時以上,真的沒有辦法想像!

難怪最近Yahoo公司要裁員,個人的感想就是,一個沒有與網友互動空間的「網路公司」,還能在Web 2.0的世界下生存嗎,看看Google到處都是和公司直接討論的空間,例如說「Google網上論壇」。
好吧,也許我對Yahoo瞭解的實在太少,難怪不知道如何和他們互動吧。

呼~不管如何,還是把想法寄給Yahoo了。

2008年1月27日

把多張圖放在一頁,同時轉成PDF檔

如何把多張圖片放在一頁裡面,這算困擾我很久的問題,當然我知道Greplay可以用,大致上我也會用,但我總覺得預設的template不是很好用,終於在SUGI發現了Kevin P. Delaney寫的文章,不只簡單教了如何把多張圖片放在同一頁上,還教如何將圖片轉成PDF檔。

主要參考文章:
Multiple Graphs on One Page: The Easy Way (PDF) and the Hard Way (RTF)
ODS LAYOUT: Arranging ODS Output as You See Fit

先提如何將圖片放在一個PDF檔案中,其實語法相當的簡單。可以分為三步驟,1:打開一個PDF檔案。2:想要畫的SAS圖形。3:完成繪圖把該PDF檔關閉。

1:ODS PDF file="你想放的路徑(例如:c:\temp\graph.pdf)";
2:SAS/Grasph 語法;
3:ODS PDF close;

來做一個簡單的範例,先產生一堆在不同自由度下T分佈的機率吧。

data t (drop=i j);
   do j=1 to 30 by 2;

      do i=-3 to 3 by 0.1;
         df=j;
         x=i;
         y=probt(x,df);
         output;
      end;
   end;
run;


接者用Gplot把圖畫出來然後放一個PDF檔中。

options nodate nonumber;
goption reset=all device=pdfc;
ods pdf file="C:\temp\CHart\chart_test_1.pdf" STARTPAGE=never;

   proc gplot data=t;
      title2 "Chart1"; 
     where df le 15;
      plot y*x=df;
      symbol1 i=l;
   run;quit;
ods pdf close;


我這裡只把DF小於15的資料畫出來。其實這樣就簡單完成把圖形放置PDF檔中的語法了,比較特別需要說明的為device=pdfc,這指令是和SAS說在PDF檔中圖片的格式,作者覺得pdfc的格式在PDF檔案中會有比較好的效果,相對於GIF或其他的圖片格式。



過來就是重頭戲了,如何把多張圖片放在一頁,在文中作者教了至少4中方法,我這裡只說明ODS LAYOUT這個指令,Layout是在SAS V9或V9.1才有的指令,我覺得還滿簡單的,相對於要去修改Greplay的版模。

在畫圖之前先簡單說明LAYOUT我會的語法。

ods layout Start width=20cm height=28cm;
ods region x=0% y=5% width=45% height=45%;
/*SAS/Graphs 語法1*/

ods region x=51% y=5% width=45% height=45%;
/*SAS/
Graphs 語法2*/

ods region x=0% y=51% width=45% height=45%;
/*SAS/
Graphs 語法3*/

ods region x=51% y=51% width=45% height=45%;
/*SAS/
Graphs
語法4*/

ods layout end;


ODS LAYOUT:Start就是宣告我們要開始使用layout的指令,請想像我們手上拿一張白,width就是設定紙的寬度,我設定了20公分,height設定紙的高度,我設定28公分。

ODS REGION:宣告我下面的圖片要放在那個位置,
『x=51% y=5%』是說我將圖片要放在由左上角算起,整張紙X軸的第51%和Y軸上第5%的交叉點,當然這裡也可以用實際的長度單位來設定,例如x=11cm y=1cm。
『width=45% height=45%』設定我這一張圖片的大小,我是用相對長度,為整張紙45%寬度和紙長45%的長度,同樣也可以設為實際長度,如width=8cm height=10cm。
ods layout end:這也是必要語法,宣告上面的ODS LAYOUT設定結束。

依序設計這四張圖片想要放的位置就完成了,我隨意將剛剛製造的資料檔分成四張圖繪製,然後配合上面的語法就可完成我想要的工作了,完整的程式如下:


options nodate nonumber;
goption reset=all device=pdfc ;
ods pdf file="C:\temp\CHart\chart_test_2.pdf" STARTPAGE=never;

   ods layout Start width=20cm height=28cm;
   symbol1 i=l;

   ods region x=0% y=5% width=45% height=45%;
   proc gplot data=t;
      title2 'Chart1';
      where df le 5;
      plot y*x=df/noframe;
   run;

   ods region x=51% y=5% width=45% height=45%;
   proc gplot data=t;
      title2 'Chart2';
      where df ge 5 and df le 15;
      plot y*x=df/noframe;
   run;

   ods region x=0% y=51% width=45% height=45%;
   proc gplot data=t;
      title2 'Chart3';
      where df ge 15 and df le 25;
      plot y*x=df/noframe;
   run;

   ods region x=51% y=51% width=45% height=45%;
   proc gplot data=t;
      title2 'Chart4';
      where df ge 25 and df le 30;
      plot y*x=df/noframe;
   run;quit;

   ods layout end;
ods pdf close;


結果如下圖:




利用PDF輸出和ODS LAYOUT的指令理論上就在一頁可繪製出多張高品值的圖形了,當然ODS LAYOUT應用不只用於圖形,想要結合表格和圖形,也是很OK的啦。

2008年1月10日

有質感的code+ & SAS如何匯入Excel檔案

很多BLOG在說明程式的時候都會用漂亮的框框把程式『框』起來,看來我偷偷學會了。

那就順便簡單說明一下的SAS如何匯入Excel檔案吧。

如果已經很熟libname的語法,那就不用特別去背proc import的語法了,所以那下面的語法對很多人來說就是不難了

libname myxls "c:\myexcel.xls";
data temp1;
set myxls."Sheet1$"n;
run;
libname myxls clear;

這裡比較需要注意的地方在於Sheet1$後面那個「$」和"Sheet1$"n後面的「n」,其他就很簡單了。
當然這個語法不是只有那麼陽春,功能也是滿強大的,有機會在介紹了喔。

code的那個框框,我是參考[筆記]在文章裡顯示優質的程式碼區