ToTo 奇妙の冒險

去年年底，俺好友「Dr. 史考特。哈蘇」(見上圖) 從密爾瓦基打電話給我，問我有沒有「狗屎」...啊不對，是「Gauss」啦～

我跟他講這套軟體基本上從俺離開清大統研所後就再也沒有見到有人在用了。

餓死抬頭，這些年來很多人一直找我做過許多power analysis來求power或樣本數，而SAS內建的PROC POWER和PROC GLMPOWER這兩個程序已經能夠把最基本的各種假設檢定、ANOVA和GLM下的power analysis給搞定，順便還可以畫出相當精美的圖形。

其實也不算是密技啦！只是前面掛這兩個字，感覺這篇文章好像比較威似的。

這是繼 [密技] 如何在 Vista 上安裝 SAS 之後另一篇關於統套軟體安裝的教學。

將近一年前分享了這篇文章：

[密技] 如何在 Vista 安裝SAS

成為本部落格最熱門的文章。這必須要"洩洩" Vista 的高度不相容性以及 SAS 公司遲遲沒有釋出更新檔。不過這一年來，還是陸續看到一些朋友們按照文內的步驟安裝但仍舊安裝失敗的情況。於是我又回原來發佈這個方法的部落格瀏覽一下，發現有人在回文內分享了另種解決方案。

昨天在幫一個博士班學生處理 mixed model 的 model diagnostics 問題時，使用 PROC MIXED 程序並配合 ODS GRAPHICS 來產生圖形時，發現 log 裡面出現這一行錯誤訊息：

Java class generated an exception

由於兩年前我在考資格考時還有用過這道指令，記得當時可以順利生出 model diagnostics 的圖形，不過這兩年一直都沒再使用這個功能，所以完全搞不清楚這個錯誤訊息是怎麼跑出來的，因此只好去 SAS 官網上面找看看。果真看到 SAS support 網頁下有說明這個錯誤訊息的源由：

暑假買了一台筆電，裡面安裝的作業系統是 Vista Home Edition。最近想要提昇模擬程式的效率，所以決定讓這台筆電也加入模擬程式的陣容。不過從學校那邊申請的 SAS 只能用在 Vista 的 Business 和 Ultimate Premium 版本，因此一直在苦思該怎樣把 SAS 安裝上去。最初想到的方法是把筆電的作業系統換回 Win XP，可是光想到要搞定筆電的驅動程式就讓人傷透腦筋。於是我把主意轉向Linux系統，尤其是最近很紅的 Ubuntu 推出新的 7.10 版本，所以想說用 SAS Linux 版本去灌看看，結果因為不熟悉 Linux 的操作方式，到現在還沒有安裝成功。這幾天我一直用 google 搜尋 SAS Linux Version 相關的安裝訊息，但一直沒有找到明確的解決方案，可是卻無意間在一個國外的 blog 上面看到有人發表如何把 SAS 安裝在 Vista Home Edition & Home Premium 的版本上。因此，我僅將文字的安裝步驟列出來，給需要在 Vista 上安裝 SAS 的人一個參考：

步驟一：到 Control Panel（控制臺）點選 Programs。

步驟二：接著選擇「Use an older program with this version of Windows」。

步驟三：此時會啟動 Program Compatibility Wizard 視窗，有三個選項，請選擇「I want to select the program manually」，然後按「Next」。

步驟四：把 SAS 的 setup disk 放入光碟機，接著在下一個對話視窗上面把 SAS setup.exe 的路徑寫上去（或點選「Browse」按鈕來選擇）。選好後點「Next」。

步驟五：接著在選擇作業系統的選項介面中選擇「Microsoft Windows XP (Service Pack 2)」，然後點「Next」。

步驟六：下一個視窗是選擇解析度，不要選直接按「Next」跳到下一個視窗。

步驟七：勾選「Run this program as an administrator」，然後按「Next」。

步驟八：接著就會開始進入 SAS 安裝介面。SAS 可能會要求你重開機，請重開機後重新重複上述步驟來安裝即可。

圖文版可至這個網址查看：

http://lazytechie.wordpress.com/2007/03/19/how-to-install-sas-under-windows-vista/

今天很久沒有聯絡的學妹突然傳 MSN 給我，當然我知道「無事不登三寶殿」的道理，所以早就有心理準備。果不其然有 SAS 的問題，所以來我這討救兵。

這問題還滿有點挑戰性的，所以當作自我鍛鍊就答應幫她了。問題的大意如下：

有一群病人被診斷出七種症狀，但不是每個人都會出現所有症狀。考慮C7取4、C7取5、C7取6和C7取7的不同組合，共計有35、21、7和1種組合。要如何看這些組合下共計有多少病人。

由於我手邊沒有資料，所以就先用電腦來模擬一百個病人，每個病人有七種症狀（symp1~symp7），此為二分法的數值變數，1表示有症狀、0表示無症狀。另外，為了要識別不同的排列組合，又另外定義了七個文字變數（sympc1~sympc7）。當symp=1時，sympc=Y，反之則sympc=N。

接著，把symp1到symp7加總起來(Total），則此變數就可以表示每個病人出現多少個症狀。再把sympc1到sympc7用連結符號（||）給黏起來，這樣就可以替每個病人出現症狀的前後次序命名。如「NYYNYYN」表示該病人出現第二、三、五、六號症狀。

再來就是定義四個新變數（combine4~combine7）來表示病人是否出現N個以上的症狀。如果病人出現四個以上症狀（即 Total >=4，則 combine4=1，反之為零）。

把上面的變數都建立好後，就可以用 PROC FREQ 去分別對 combine4=1、combine5=1、combine6=1和combine7=1 去做次數分配表，這樣一來就可以算出每種不同組合下有多少病人了。

程式碼如下：

data test;
    array symp{7} symp1-symp7;
    array sympc{7} $1 sympc1-sympc7;
    array c{4} combine4-combine7;
    seed=1234;
    do i=1 to 100;
        do j=1 to 7;
            symp(j)=round(ranuni(seed));
            if symp(j)=1 then sympc(j)='Y';
                else sympc(j)='N';
        end;
        total=symp1+symp2+symp3+symp4+symp5+symp6+symp7;
        combine=sympc1||sympc2||sympc3||sympc4||sympc5||sympc6||sympc7;
        do k=1 to 4;
            if total>=k+3 then c(k)=1;
                else c(k)=0;
        end;   
        output;
    end;
    drop seed j k;
run;

%macro combine(indata,var,outdata);
proc freq data=&indata;
    where &var=1;
    tables combine;
    ods output OneWayFreqs=&outdata;
run;
%mend;

%combine(test,combine4,combine4freq);
%combine(test,combine5,combine5freq);
%combine(test,combine6,combine6freq);
%combine(test,combine7,combine7freq);

之前在幫長庚護理系的葉昭幸教授作一些資料分析時，曾經寫了一個 macro 以便於將 PROC FREQ 的結果拿來製造 82 個資料集。因此這個 macro 就跑了 82 次。今天我把這個程式拿來看一看，覺得有點 stupid。心想如果有更方便的方法來讓 macro 可以自動重複的把 PROC procedure 跑完的話，那就更完美了。



因此，上網路稍微搜尋一下，馬上找到答案。原始的 stupid 程式如下：



%macro variable(var,outdata);

proc freq data=symptom.msas noprint;

tables &var /out=&outdata;

run;

data &outdata;

set &outdata;

drop &var percent;

rename count=&var;

run;

%mend;

/*Distress*/

%variable(distress1, symptom.distress1);

...

...

...（省略）

%variable(distress30, symptom.distress30);



/*Severity*/

%variable(severity1, symptom.severity1);

...

...

...（省略）

%variable(severity30, symptom.severity30);



/*Freq*/

%variable(freq1, symptom.freq1);

...

...

...（省略）

%variable(freq22, symptom.freq22);



後來改成這樣：



%macro variable(var,num);

%do i = 1 %to #

proc freq data=symptom.msas noprint;

tables &var&i /out=&var&i;

run;

data &var&i;

set &var&i;

drop &var&i percent;

rename count=&var&i;

run;

%end;

%mend;

%variable(distress,30);

%variable(severity,30);

%variable(freq,22);



主要修改的地方有兩個：

一、把原本要存入symptom這個library的檔案全部改存到work裡面去，這樣這些暫存的檔案在經過下一步驟的合併之後並離開SAS後就會消失，不會佔用硬碟空間。

二、在 PROC FREQ 和 data step 外面加上一個 do loop，然後把要執行的圈數設成一個新的 macro 參數叫做 num，這樣就能拿來控制要重複執行的次數。

我在二月八日發表了一篇關於如何大量重新命名新變數的文章（How to concatenate two variables with array）。當時寫完覺得自己很屌，但今天看到一篇文章後頓時覺得自己很遜。原先的設計是要將上千個變數切成兩半，前半段包含舊變數的第一個字母，並且重新命名為「舊變數a」，後半段包含舊變數的最後一個字母，並且重新命名為「舊變數b」。



原始的程式主要是將變數名稱另存一個新檔，針對那個檔來做重新命名的工作，然後再把新變數名稱存出來貼到另一個程式去做切割舊變數的動作。但有個神人寫了一個可以自動重新命名的 macro，讓上述動作在數行程式碼之間迅速完成！！



這篇文章可至 SUGI 的官方網站下載，但我已經大部分的文件弄成一個簡單的中文說明放在我的 SUGI CLUB 裡面。有興趣的人可以去看一看。



舊的程式碼我就不列了，可以到之前的文章看。新的程式碼如下：



%include <把add_string macro的路徑放入>;

%let old_vars=<裡面把所有舊變數名稱全部放入>;

%let N=<放入變數總數>;



data newset;

set oldset;

array old[*] &old_vars;

array recode1[*] %add_string(&old_vars, a);

array recode2[*] %add_string(&old_vars, b);

do i = 1 to &N;

recode1[i]=substr(old[i],1,1);

recode2[i]=substr(old[i],3,1);

end;

run;



裡面的 %add_string 就可以將所有 old_vars 裡面的舊變數後面多加一個 a 或 b。然後這新的變數放進兩個不同的陣列，再用 substr 指令去割舊變數就完成了。

-----

我偶爾會去逛一個叫做「統計生活館」的網站。這個大概是我覺得國內最完整的一個在介紹統計學的網站。裡面有一個討論區，經常有人會上去問些問題。如果我知道答案，也會在那邊幫點小忙。



前幾日有位仁兄問了一個問題。他用亂數生成的方法模擬了一個資料庫包含1068筆資料，然後用 SAS 的 PROC CLUSTER 做 Cluster analysis，最後再把樹狀圖畫出來。這個流程其實是很簡單的，但是他要重複上述的動作 10000 次。他問程式該怎樣寫。



我本來想說用個 Macro 來搞定，但是卡在怎樣在亂數模擬變數那一階段請 SAS 生出 10000 個資料庫，因為 Do loop 不能掛在 Data procedure 的外面。後來想到一個簡單的方法，不用 Macro 就可以搞定。



原始程式如下：

data simulation;

do no=0001 to 1068;

array Q[*] Q01-Q25;

do j=1 to 25;

k=rand('uniform');

Q(j)=ranbin(1,1,k );

end;

output ;

end;

proc distance data=simulation method=djaccard absent=0 out=distjacc;

var anominal(Q01-Q25);

run;

proc cluster data=distjacc method=ward rsq rmsstd pseudo outtree=tree;

var _numeric_;

run;

proc tree data=tree out=out;

run;

修改後的程式如下：

data simulation;

do index=1 to 10000;

do no=0001 to 1068;

array Q[*] Q01-Q25;

do j=1 to 25;

k=rand('uniform');

Q(j)=ranbin(1,1,k );

end;

output ;

end;

end;

run;



proc distance data=simulation method=djaccard absent=0 out=distjacc;

by index;

var anominal(Q01-Q25);

run;



proc cluster data=distjacc method=ward rsq rmsstd pseudo outtree=tree;

by index;

var Dist1-Dist1068;

run;



proc tree data=tree;

by index;

run;

差別主要在生成資料那一步驟，只需要在外面多掛一個迴圈，讓裡面生亂數的迴圈再跑 10000 次，然後資料裡面會多了一個叫做 index 的變數。這個變數從 1 到 10000，就等於是我把 10000 個資料庫全部放在一起。之後的任何 Procedure 就可以用 by statement 來搞定。稍微需要注意的是，在 PROC CLUSTER 的 Var statement 中不能再用通用變數 _numeric_，而要改成Dist1-Dist1068。原因是若用 _numeric_ 通用變數的話，會把 index 也放進去分析，這樣會導致距離矩陣不是 n x n，而是 (n+1) x n。目前唯一還不算解決的是最後的 Proc tree 程序，不能再使用 out option，我想這是可以解決的，只是懶得再想而已。:P













-----