使用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執行。

熬夜看書

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

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