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













-----

 昨日在幫 Ann 完成 GEE model diagnostic 之後，她雖很心滿意足地回去繼續她的論文，但留了一個洞給我繼續挖。任何的 model diagnostic 中，每個觀測值都會產生一些統計值，如 Cook's D 或 Leverage，都需要畫一些圖表來顯示每個統計值的高低，讓讀者可以很明白地看出哪些值是 influential data。因此，除了要描點之外，還要特別把幾個 influential data 的編號標示出來。這種圖示我就不知道該怎樣去做了。因此，趁著寫感謝信給 Dr. John Pressier 之餘，順便在信中問他那種圖該怎樣畫（因為他在 paper 裡面秀了很多張這類的圖表）。但如我一開始的直覺假設，他那些圖果然是用 S-PLUS 弄出來的。這時要我去研究怎樣用 S-PLUS 去做圖是不可能的，因為我從上個世紀（也就是我念研究所的時代）就不喜歡用 S-PLUS，即便他後來的功能改的很強。因此，正當我在翻閱 SAS/GRAPH 手冊尋求任何配套程式寫法時，John 的另一位 co-author，Dr. Bradley Hammill，寄了一封信給我，教我如何用 SAS 來畫那種圖。原來只要用 POINTLABEL 這個指令就可以標示想要的數字出來（但 SAS/GRAPH 原文手冊裡面居然沒有列出這個指令！）。在此分享一下 Brad 提供的範例：

* Generate some random data for plotting;

data a;

    do idnum = 1 to 100;

        sres = rannor(23456);

        output;

    end;

run;

* Label all points in index plot with ID number;

symbol v=dot c=black pointlabel=(h=2pct "#idnum");

proc gplot data=a;

    plot sres * idnum;

run;

* Setup a new var that only has the ID number if abs(SRES) > 2;

data a;

    set a;

    if abs(sres) > 2 then idhigh = idnum;

    else idhigh = .;

run;

* Label only high-SRES points in index plot with ID number;

symbol v=dot c=black pointlabel=(h=2pct "#idhigh");

* W/out this option, missing labels in plot will appear as .;

options missing=" ";

proc gplot data=a;

    plot sres * idnum;

run;

在問卷調查中，為了避免輸入錯誤，比較有經費的計畫大都會一次雇用兩個人手做double
entry的程序，就是說一份問卷由不同的人各輸入一次，然後比較看這兩個人有沒有輸入錯誤。如果兩人的輸入完全一致，就比較有高度的可靠度來說這筆資料是正確地輸入的。當然也可能會發生兩人同時出錯在同一個地方，但這機率很小。




問題來了，輸入好的資料如何交由電腦來做自動比對。SAS很貼心地有一個PROC
COMPARE的程序可以進行交叉比對。這個工具很有威力的理由是，比對數值變數其實還是小case，更強的地方是可以比對文字變數，而且可以精確到計算出裡面有多少個typo。但其實我們不需要知道太詳盡，只要兩兩比對率沒有到100%，就可以肯定一定其中某個人有錯誤輸入。但SAS沒有辦法自動校正，畢竟連我們也不知道哪個人輸入的才是正確的版本，所以只能將有問題的資料重新調出來用人工修正。




可以在下面這個連結找到完整語法和幾個實用範例：



http://www.sussex.ac.uk/its/help/guides/sas/proc/z0057814.htm



如果是SPSS資料格式，可以用下面程式把SPSS檔案叫進SAS：



proc import datafile=xxx.sav out=xxx dbms=SAV replace;



run;

-----

Yesterday midnight, when I enjoyed watching X-files in the living room, a call disturbed me. That is my friend who is in the trouble of doing some data management in her dissertation. She has some character variables which looks like:



A/A

C/C

G/G

T/T



She needs to create two new variables. The first new variable contains the first word in old variable, and the second variable contains the last word in old variable. Both of them should also be converted into numeric variables. The criteria is A for 1, C for 2, G for 3, and T for 4.



Her first request is very easy, we can use a special function named "SUBSTR" to extract some word in specific location. The names of two new variables should be names as, for example:



old variable = rs1234567

1st new variable = rs1234567a

2nd new variable = rs1234567b



That's not difficult as well because we can just define the new names in data procedure. However, she said, "I have 3000 variables."



That's is a very very big dataset. If she define new variables by herself in SAS, she needs to write 6000 new variables. She thought a SAS macro could probably save her time, but it is still not efficient. Finally, I found out a good way to finish her job.



Step 1: Export her dataset into Excel, and the first line of the Excel file lists all variable names.

Step 2: Cut down those 3000 variables, and create a new dataset with 3000 temporary variables containing those 3000 variables.

Step 3: Put all 3000 temporary variables in an array named "old", and create 6000 new variables in another two array named "newa" and "newb".

Step 4: Use "Do...End" to concatenate 3000 old variables with "a" and "b", then put them into array newa and newb, respectively.

Step 5: Export this dataset into Excel, so all 6000 new variables are created.



The following program is a simplified example with only 5 old variables:



data test;

input (old1-old5) ($) ;

array old[5] $ old1-old5;

array newa[5] $12 newa1-newa5; /*adjust the lenght of new variables longer*/

array newb[5] $12 newb1-newb5;

do i=1 to 5;

newa[i]=right(old[i])||"a"; /*keep old variable right to prevent from a gap between old variable and a*/

newb[i]=right(old[i])||"b";

end;

keep newa1-newa5 newb1-newb5;

cards;

rs12 rs34 rs56 rs78 rs90

run;

proc print;run;

proc export data=test outfile="C:\Temp\newtest.xls"; run;

quit;

This is a new category in my blog which includes some practical issues in using SAS. Even though I use SAS for many years, I still can't say I totally understand SAS. The depth of SAS programming is far away from what we can imagine. The powerful of SAS is absolutely, but the complication of it also troubles many users. Many Taiwanese students in UNC-CH often ask some SAS problems. I'd be glade to help them because they also give me the chance to face some challenges that I didn't encounter before. I am sincere to share my solutions in this blog. It is a sort of record, but I rather regard it as a cumulative experience of using SAS in my life.