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

cchien 發表在 痞客邦 PIXNET 留言(13) 人氣()

今天很久沒有聯絡的學妹突然傳 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);

cchien 發表在 痞客邦 PIXNET 留言(2) 人氣()

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

在迴歸分析中,最麻煩的事情莫過於模式鑑定。其中,又以模式是否違反變異數同質性是許多人感到最困難的。變異數同質性的道理基於殘差的變異數需呈現一個固定常數。但要達到這個需求,就必須要靠依變數的變異數 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

cchien 發表在 痞客邦 PIXNET 留言(1) 人氣()

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

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

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

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

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

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

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

cchien 發表在 痞客邦 PIXNET 留言(0) 人氣()

之前在幫長庚護理系的葉昭幸教授作一些資料分析時,曾經寫了一個 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,這樣就能拿來控制要重複執行的次數。












cchien 發表在 痞客邦 PIXNET 留言(4) 人氣()

我在二月八日發表了一篇關於如何大量重新命名新變數的文章(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 指令去割舊變數就完成了。














-----

cchien 發表在 痞客邦 PIXNET 留言(0) 人氣()

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



前幾日有位仁兄問了一個問題。他用亂數生成的方法模擬了一個資料庫包含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













-----

cchien 發表在 痞客邦 PIXNET 留言(0) 人氣()

 昨日在幫 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;














cchien 發表在 痞客邦 PIXNET 留言(3) 人氣()

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 的存在。來護理系工作真是惠我良多啊。












cchien 發表在 痞客邦 PIXNET 留言(11) 人氣()

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















-----

cchien 發表在 痞客邦 PIXNET 留言(0) 人氣()