使用SAS 9.2 繪製ROC曲線和比較

在2008年4月的時候，曾經寫過一篇繪製ROC曲線和比較兩條ROC曲線下面積的小文章，這篇是以MACRO來畫圖和比較，而在新版的SAS 9.2中已經可以於PROC LOGISTICS直接繪製高品質的ROC曲線(Receiver Operating Characteristic Curves, ROC Curves)。

在SAS 9.2中ODS GRAPHICS支援的程序增加許多，其中也包括PROC LOGISTICS。

先產生一筆虛擬的資料 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 52  0 41 66
0 43 52  0 36 54  0 37 46  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 68  1 64 53  1 62 41  1 39 57
1 61 51  1 61 55  1 57 49  1 79 67
;
run;

接下來就直接套用繪製的ROC曲線的語法 ods graphics on;
proc logistic data=age_data plots(only)=(roc) rocoptions(id=prob) ;
model disease (event='1')=age;
run;
ods graphics off;

ods graphics on：為宣告要繪製ods graphics的圖，而ods graphics off宣告結束，幾乎ods系列的語法都是類似的語法。
plots(only)=(roc)：在proc logistic程序中，新增了幾個選項，plots為宣告要繪製統計圖，括號中的only為要求SAS只要畫出所需的圖就好，而在後面的【 (roc) 】和SAS說我們需要的圖為ROC曲線。
rocoptions(id=prob)：表示在圖形上的點要放哪些資訊，目前有兩種可以選擇，【prob】為predicted probability ，而【obs】為observation number。
這幾個簡單的指令便可畫出高品質的ROC曲線。 而且還順便秀出AUC（Area Under Curve）的值。

 

如果在model後面多放一個變項（BW），則畫出此模型（age和bw的模型）的ROC曲線。 ods graphics on;
proc logistic data=age_data plots(only)=(roc) rocoptions( id=prob) ;
model disease (event='1')=age bw;
run;
ods graphics off;

此模型的ROC曲線為下圖：

不過也許大家比較想看到的是，單圖age或bw的圖，或是一張圖裡面同時有age和bw「個別」的ROC曲線，而不是模型的結果。

ods graphics on;
proc logistic data=age_data plots(only)=(roc) rocoptions( id=prob) ;
model disease (event='1')=age bw;
roc '年齡' age;
roc '體重' bw;
run;
ods graphics off;

roc：個別宣告要繪製哪些變項的圖形，這裡我想個別看age及bw的ROC曲線。有單引號包住的文字，為宣告此圖形中變項的名稱，分別為「年齡」及「體重」。不過如果只使用這個語法的話，SAS同時還會多給一張model的圖，和age+bw+model的圖。

上面這張圖，雖然多了一個model的曲線，但這不一定是每次都想見到的，我的經驗，大多臨床工作者，都不太瞭解model這條曲線的意義，即使解釋了很久都還是不會用，因此單純的age和bw的比較也許就很足夠了。這個時候將程式修改成下面的寫法即可。

ods graphics on;
proc logistic data=age_data plots(only)=(roc) rocoptions( id=prob) ;
model disease (event='1')=age bw/nofit;
roc '年齡' age;
roc '體重' bw;
run;
ods graphics off;

其實就只是多了一【nofit】，這會讓SAS不要做整個model估計，只分別對age和bw做兩個logisitc regression，此時Log視窗也會提醒，「NOTE: The NOFIT option is specified.  No parameters are estimated and all other MODEL statement options (except FIRTH ,LINK= ,NOINT ,OFFSET= ,ROC ,and TECHNIQUE=) are ignored.」。結果如下圖：

最後就是比較age和bw的AUC有無統計上的差異了。

ods graphics on;
proc logistic data=age_data plots(only)=(roc) rocoptions( id=prob  ) ;
model disease (event='1')=age bw/nofit;
roc '年齡' age;
roc '體重' bw;
roccontrast  reference('體重')/ estimate e;
run;
ods graphics off;

roccontrast：比較ROC曲線。結果如下：

SAS 9.2此項新功能，我自己覺得使用上比MACRO直覺許多！

Firefox 3.5 中文版 (100% 純天然有機軟體)

咦…這個版本就突然上市了，這個很號稱很快的版本(Firefox 3.5測試：速度是Firefox 3的兩倍快)，我今天才用到，好吧，反正也沒有太大的感覺，不過終究是花很多時間研發，不斷的延期又延期的Firefox 3.5，還是要用力推一下啦！

Firefox 3.5中文版

奔牛節

某一天看到報紙「奔牛節的街頭藝術品被破壞」，才知道台中市有這個項活動，這應該是屬於街頭裝置的藝術品吧，感覺還滿有趣的。

到了國美館，在綠色的草地上就出現了數頭顏色鮮豔的大牛！

分別給其中隻三來個特寫。

這幾隻牛上的小東西真的是很有趣，牛那麼大，在小時候的印象總是需要一些小動物的～

其他室外的牛也是很有趣的，甚至躺著的牛 XD

除了室外，室內的展區也是有很多奇妙的牛。這一隻應該也是牛吧，有一堆小牛圍在旁邊。
（照片中的為國美館的服務人員，還很辛苦的在數小牛數目有無正確。）

這裡也是有一堆小牛，由一個大圓石走出來。

當然室內也是有大牛，有閃亮花俏的大牛，有小人控制的大牛…

還有「天下第一味」的…牛肉麵？

有興趣想看這群牛的，時間不多了，展出到2009/7/5。

在SAS中使用已經整理好的統計量

實務上常遇到研究者只提供已算好的數字(統計量)，例如列聯表中的個數和百分比或兩樣本的人數、平均數和標準差等，如果要利用這些統計量做某些檢定，例如卡方檢定或雙樣本T檢定等，可使用PROC FREQ和PROC TTEST處理僅有統計量的問題。

PROC FREQ

有兩群病患分別使用A藥和B藥，再觀察兩組實驗組症狀有無改善，結果如下表，如研究者想問A藥和B藥病患症狀改善的比率有無不同，可採用卡方檢定(Chi-square test)

	A藥	B藥
症狀無改善	9	12
症狀改善	21	18

當然第一步要先把資料輸入至SAS中，程式可以這樣寫：

DATA TRIAL_1;
INPUT DRUG $ EFFECT COUNT;
DATALINES;
A 0 9
A 1 21
B 0 12
B 1 18
;
RUN;
PROC SQL;
SELECT *
FROM TRIAL_1;
QUIT;

接下來是卡方檢定的部分

PROC FREQ DATA=TRIAL_1;
WEIGHT COUNT;
TABLE EFFECT*DRUG/CHISQ;
RUN;

PROC TTEST

同樣兩群病患使用A藥和B藥，使用後測量某個臨床指標（連續型變項），研究者欲使用T檢定，探討兩組病患的臨床指標平均值有無差異。資料如下表

	A藥(N=30)	B藥(N=30)
臨床指標	89(12)	95(16)

mean(std)

因為這裡我們只有拿到兩組病患的人數、平均值和標準差，所以在輸入資料時要稍微複雜一點，必須符合SAS程式的規定，一定要有一個變項名稱為『_TYPE_』或『_STAT_』，在這個變項中，必須要放入‘N’、‘MEAN’、 ‘STD’這三筆資料，資料輸入程式如下

DATA TRIAL_2;
INPUT DRUG $ _TYPE_ $ EFFECT;
DATALINES;
A N 30
A MEAN 89
A STD 12
B N 30
B MEAN 95
B STD 18
;
RUN;

接下來執行PROC TTEST程序就簡單了，和一般寫法沒有太大差異

PROC TTEST DATA=TRIAL_2;
CLASS DRUG;
VAR EFFECT;
RUN;

以上程式可於SAS 9.1.3和SAS 9.2執行。

熬夜看書

數統最後一學期的期末考終於到了，
每次考試幾乎都要熬夜看書，可是又不敢不睡，
隔天還有許多待續工作要完成。

胡言亂語中，做個紀念吧！
明天（今天？）也要考試的伙伴加油吧！

Stay Hungry, Stay Foolish

真是一場發人深思的演說。

這個片段的起頭就說到了死亡，中間也提及了Your time is limited..，死亡是人類重要，但又深深隱身在人群之中，舉個例，某個工安事件發生了，新聞|社會最關注的就是有無人員死亡，誰要對這個死者負責，家屬萬分哀痛等。每個人都知道死亡一定會公平的發生在每個人身上，有可能就是下一秒，但平常生活中，我們卻常刻意忽略這個重要的議題，假裝死亡不存在，只有在她真的發生時才覺得重要，但對個人來說太晚了。

Steve Jobs曾經得到死亡率超高的胰臟癌，但他活下來了，可能因此體會到死亡的重要，重新去排序生活中的重要的事情，他也提到了不要被傳統的教條束縛，要多聽聽自己的聲音，在死亡前多完成「自己想做的事情」，也許只是件傻事。（不過我相信Steve Jobs在發現胰臟癌之前就懂得尊重自己）

Who is Steve Jobs，請自行查閱Wiki的資訊吧（http://zh.wikipedia.org/w/index.php?title=Steve_Jobs&variant=zh-tw）

完整的演說-1（只有英文字幕）：

完整的演說-2（只有英文字幕）：

遇見未來的

這4個研究所的小朋友，是我於2007年就期待看到的一群，當時還在猜測他們是誰，會哪裡考來。

現在他們已經入學超過了半年，也認識快9個月了，我無法否認，對這群小朋友有不一樣的感覺，因此也花了比較多時間在他們身上，把一些我覺得重要的東西讓他們瞭解，但當下一屆出現時，也許我這個熱情就會修正了。終究，他們是我曾經期待的一群。〔將自己放在未來，總有減壓的和跳離現實痛苦的感覺，享受過就好了，不要沈迷在未來與過去之間。〕

我沒讀過國立大學，我不知道那邊的資源，但我的過往經驗告訴我，他們真的很幸福，老師數比學生多，每次諮詢報告老師全到的陣容，還真是完整，除了老師的資源，我同事總是能給這群小朋豐富又實務的資源，這也許就是某種「產學合作」吧。如果當時有這樣的資源，也許我就不用花那麼多時間自修了，想當初自己嘗試「認識」MACRO或是說入門，就花了好多好多時間了。

兩隻笨貓終於可以重逢了

橘子貓今天終於可以在遙遠的那一方見到另外一隻笨貓了，雖然這一切真的是很突然。
永遠祝福那可愛的兩隻笨貓。

想念的影片

這個已經不知道在那個世界的貓，這是唯一能夠想起她身影的方式了，
一位很愛她的女生曾經說過「她是全世界最可愛的貓」。

當我想要一個人的時候，我都會去她休息的地方坐坐，曬曬太陽，
如果你們看到一個人，默默在那邊舔者傷口，請不要打擾他，
也許他身旁正有一隻全世界最可愛的貓。

SAS Macro ( 巨集 ) 簡易教學

SAS的巨集雖然屬於P3等級的內容，不過對於實際分析資料的使用者來說，應該是相當實用的語法，且必要學習的。巨集沒有想像中的難，只是多了一些平常不會使用的符號，例如『%』或『&』等，這一篇教學將會舉一個簡單的例子，說明巨集的使用方式。


先來假設一個狀況，有4個資料檔{a1,a2,a3,a4}，每個資料檔中均包含4個變項{y1,x1,x2,x3}，而我想以PROC REG來分析4個資料檔中「y1和x1」、「y1和x2」及「y1和x3」的相關性。

第一步：先建立4個資料檔並包含3個變數。
先以do loop建立4個資料檔，而y1、x1、x2、x3的值將隨機由uniform分佈中抽取。程式如下：

data a1 a2 a3 a4 ;
do i=1 to 4;
do j=1 to 100;
y1=uniform(123);
x1=uniform(1234); x2=uniform(12345); x3=uniform(12345);
drop i j;
if i=1 then output a1;
else if i=2 then output a2;
else if i=3 then output a3;
else if i=4 then output a4;
end;
end;
run;


第二步：用傳統方式寫PROC REG程式。
當然可以傳統語法，逐筆將資料檔名稱和變項寫入程式中，以資料a1為範例：

proc reg data=a1;
model y1=x1;
run;
proc reg data=a1;
model y1=x2;
run;
proc reg data=a1;
model y1=x3;
run;

然後再將這一小段程式複製3次，將data分別改成a2、a3、a4即可。

第三步：改用巨集的語法。
其實可以很清楚的找到一個邏輯性；PROC REG的程式每次只需要修改下面藍色的部分就好，
proc reg data=a1;
model y1=x3;
run;

因此抓到這個邏輯以後，就可以寫出一段小小的巨集。

%macro My_reg (mydata, myvar);
proc reg data=&mydata;
model y1=&myvar;
run;
%mend My_reg;

這裡出現了奇怪的符號「%」和「&」。『%mcaro』為告訴SAS要開始寫一段巨集程式，有開始就會有結束，如同do loop會有do...end一樣，在巨集語法中結束即為『%mend』。

在『%macro』後面的My_reg為宣告這段巨集程式的名稱，為何要宣告，因為我們可能會同時使用上很多個巨集程式，如沒有宣告，SAS程式也搞不清楚現在要使用那個巨集程式。

接下來的『(mydata, myvar)』，就是在巨集有幾個「變數」會換（前面PROC REG藍色的2部分），並給這些變數一個名稱。

在被巨集程式的包裹住的語法，絕大部分都和之前寫的程式一樣，只有多了『&mydata』和『&myvar』這兩個，這個由『&』開頭的兩個部分稱為『巨集變項』，這為等一下要被取代的部分。

第四步：使用My_reg巨集。
已經寫好的巨集叫『My_reg』使用時，在巨集名稱前多加『%』，寫法如下：

%My_reg (a1,x1);
%My_reg (a1,x2);
%My_reg (a1,x3);

%My_reg (a2,x1);
%My_reg (a2,x2);
%My_reg (a2,x3);

%My_reg (a3,x1);
%My_reg (a3,x2);
%My_reg (a3,x3);

%My_reg (a4,x1);
%My_reg (a4,x2);
%My_reg (a4,x3);

在巨集名稱後面（）要放入取代巨集變項的「文字」，記得要以『,』區隔。執行後SAS便會將所寫入的文字逐步取代到『My_reg』巨集中，然後執行。

第五步：多一點的寫法。
這裡主要希望SAS能自己重複執行這個結果，不需要如第四步一樣手動改文字，此程式先不在這裡說明了，有興趣的請自行研究了。

%macro My_reg1 ;
%do i=1 %to 4;
%do j=1 %to 3;
proc reg data=a&i;
model y1=x&j;
run;
%end;
%end;
%mend My_reg1;
%My_reg1;