目前日期文章:200707 (2)

瀏覽方式: 標題列表 簡短摘要
今天很久沒有聯絡的學妹突然傳 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) 人氣()

找更多相關文章與討論