ToTo 奇妙の冒險

暑假買了一台筆電，裡面安裝的作業系統是 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);

話說很久沒發正經一點的統計文章了。。。。

在迴歸分析中，最麻煩的事情莫過於模式鑑定。其中，又以模式是否違反變異數同質性是許多人感到最困難的。變異數同質性的道理基於殘差的變異數需呈現一個固定常數。但要達到這個需求，就必須要靠依變數的變異數 Var(Y) 不會隨著其期望值 E(Y) 的變動而變動。在一般大學的迴歸分析課程裡面，老師們通常會教學生畫一張 standardized residual (or  studentized residual) v.s. predicted value 的散佈圖。如果這張圖裡面的點是呈現水平帶狀散佈，就表示模式的變異數同質性沒有違反，如下所示：


如果是呈現某種特殊趨勢，如馬鞍狀或扇形，就表示模式違反了這個假設。圖形如下所示：

但是「看圖說故事」人人都會，有時候圖形明明出現一點趨勢，但還是可以硬凹講成水平帶狀。反正這是自由心證，也沒有科學數據來證實「多少程度水平帶狀才叫做真的水平帶狀」。此外，樣本太少也是導致誤判的因素之一。因為可能在樣本數少的時候呈現水平帶狀，但誰知道等到樣本數增加時會不會變成扇形。。。。

如果這個假設能夠使用一些真正可算出 p-value 的數據來幫助判別，那一定相當方便，也可避免硬凹的情況。不過，以前的老師沒教過，教科書上也鮮少提到。其實還是有人發明出來的，其實這個檢定早在 1979 年就由 Breusch  和 Pagan 發表出來了，就取名為 Breusch-Pagan Test。

Breusch, T. and Pagan, A. (1979), ``A Simple Test for Heteroscedasticity and Random Coefficient Variation," Econometrica, 47, 1287-1294.

接著 White 也發明了另一個檢定來進行變異數同質性檢定（取名為 White's Test）。

White, H. (1980), ``A Heteroscedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroscedasticity," Econometrica, 48, 817-838.

SAS 並沒有把這兩個檢定放進 PROC REG 裡面，而是放在 PROC MODEL 中。有興趣的人可以到下面的連結去 copy 程式碼。

A Simple Regression Model with Correction of Heteroscedasticity

好像很久沒有發 statistical consulting 的文章，因為這幾年在護理系什麼樣大大小小的統計問題都遇過了，後面的挑戰性就越來越低，所以也沒遇到什麼麻煩。

不過上星期一個 clinical staff 的人寫 email 給我請求統計協助。她說她找了兩個我們系上的學生問，可是他們都不會。我想說這是什麼樣的天大難題，居然讓兩個生統系博士班的學生束手無策，所以我就先答應了她，然後把她 email 給我的一個文件下載來看看。

結果，短短兩頁的說明文件，讓我看了一個多小時才看的懂。不是她寫的用字太難，而是整個敘述相當沒有條理。後來經過我額外整理出十點條目才把她的研究和問題釐清。

簡單來講，這是一個針對膝蓋手術後的療程試驗，有實驗組和對照組。實驗組用的是新療程，對照組當然是用傳統方法。但問題時，這六十幾個病人有些在進入這個臨床試驗時就已經用了額外的療法，有些病人則是在進入試驗後，自己額外去使用別的療法。總之這些病人的疼痛效應完全被跟實驗不相干的療程給干擾到了。

而基本上這種情況幾乎可以預想到統計分析的結果絕對會有偏差，而他們也不可能再用原始的組別（實驗組vs對照組）去作檢定或分析。但要他們去重新收資料已經是不可能了，畢竟金錢和人力都花下去了，現在就看有沒有辦法在 study design 上作任何改進。 

當時我第一個想到的就是重新分群，但每個人有各自不同複雜的情況，即便是這位 staff 寄給我資料我還是看不太懂。今天她跑來親自解釋給我聽這些病人的狀況，我再依照她的敘述把上星期整理的十個條目減為八個條目，然後設定八個 binary 變數來代表這八個條目。最後把這八個變數丟去跑 cluster analysis 後畫出像這樣子的樹狀圖：

然後就叫她把 output 帶回家自己去歸類並重新替新組別命名，而這部分我就幫不上忙了，因為只有她最瞭解每個病人的特徵。我請她去朝找尋同一組人有沒有共同治療特徵，再用有意義的名稱去定義新分群。而這也是目前我想到唯一能讓她繼續把接下來的分析弄完的方法。希望老天爺夠照顧她能讓她真的找到好的新分群。

之前在幫長庚護理系的葉昭幸教授作一些資料分析時，曾經寫了一個 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 ;

%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;

Ann Jessup，一個護理系即將畢業的博士班學生，其 committee member 之一正是我老闆 Dr. Bangdiwala。由於這一陣子我老闆跑去智利不知道去幹啥，因此 Ann 的博士論文中的統計部分就落到我的頭上了。



她的論文題目裡面用到所謂的 GEE model，這對我來說並不是太困難的事情。但之後卻發生了連我 RA 老闆 Mark 都不知道該怎樣解決的麻煩。第一個麻煩是，Ann 要從下面這三個 GEE model 中挑選最好的一個

：



Y=X1+X2+X3+W

Y=X1+X2+X3+V

Y=X1+X2+X3+Z



這三個模式除了最後一個獨立變數不一樣之外，其餘都一樣。簡單的作法就是分別做三個模式，如果 W,V,Z 其中一個顯著而另兩個不顯著，自然就可以認定含有那一個顯著變數的模式是比較好的。



但是，如果有兩個以上顯著呢？



當然，我沒有一開始就直接去跑這三個模式，所以不知道結果會是如何，但總是要假設一下各種可能性。如果是在 mixed model，這種 non-nested 的情況就直接用 AIC 或 BIC 這兩個數據來判定，越小的越好。那 GEE model 有沒有類似的準則可以判定呢？答案是有的！GEE model 也可以用 AIC，但不幸的是，目前 SAS 9.1.3 版的 PROC GENMOD 沒有內建這個功能。（聽說9.2版會內建了，但不知啥時才會 release 出來）。於是我就上網去看看有沒有其他人發明相關的方法，果真發現一個叫做 QIC 的數據，而且 SAS 公司已經把他的 macro 程式放在網上任意下載。所以我就拿 SAS 提供的程式搞定了這一部份。



第二個麻煩是，當選好了 GEE model，要如何做模式檢測來挑出可能的離群值或影響值呢？其他的統計模式都有些方法可以來做模式檢測，就唯獨 GEE model 沒看過什麼模式檢測。我翻了以前上課的講義，老師也沒講這一段。跑去問 Mark，他居然也不知道。此時，又得去拜一下 Google 大神了。沒想到不只有人把這段理論給弄了出來，連程式也寫好了。不過更幸運的是，發表這些 paper 的學者中，有一個共同名字，叫做 John Pressier。而這位仁兄，正巧就是我的 committee member 之一。



於是乎當然就很高興的把 paper 稍微讀一下，接著就馬上下載他寫好的 macro 來用。中間遇到一個問題，那就是他的程式裡面並沒有辦法設定某些變數是離散的。後來直接寫信去問他，他叫我把那些變數改成 dummy variable 就搞定了。



所以，就在這麼許多巧合和運氣之下，Ann 的兩大問題就這樣被我解決了。同時間自己也學到了很多，其實我應該要更感謝 Ann 才對，否則我永遠也不會知道 QIC 和 GEE model diagnostic 的存在。來護理系工作真是惠我良多啊。

在問卷調查中，為了避免輸入錯誤，比較有經費的計畫大都會一次雇用兩個人手做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;

-----