公告

[公告]
2014/01/17
由於已經是faculty的關係,不太有足夠時間寫部落格。因此更新的速度會相當緩慢。再加上近幾年來SAS GLOBAL FORUM沒有出現讓我覺得驚艷的技術文件,所以能分享的文章相對也減少許多。若有人推薦值得分享的SAS技術文件,請利用『問題討論區』告知。

2013/07/19
臉書留言板的功能因為有不明原因故障,因此特此移除。而intensedebate的留言板因管理不易,也一併移除。目前已經開啟內建的 G+ 留言系統,所以請有需要留言的朋友,可直接至『問題討論區』裡面留言。


2012年3月31日 星期六

Computation of Correlation Coefficient and Its Confidence Interval in SAS®

原文載點:http://www2.sas.com/proceedings/sugi31/170-31.pdf

在一些特殊的情況下,相關係數需要搭配信賴區間。PROC CORR程序不像其他程序語法一樣在後面的選項加上一個 CL 就可以自動生出信賴區間,而是需要用別的選項語法才能產生。而且,這個功能是只有 SAS v9.1.3 之後的版本才有,舊版本則需要靠使用者自己寫程式去算。David Shen 和 Zaizai Lu 在 SUGI31 發表了一篇技術文件來解說如何在新舊版本下算出相關係數的信賴區間。
。理論
由於相關係數並不是常態分配,所以要經過三個步驟,先去做變數變換,算出變換後的信賴區間,然後再轉回去。步驟如下:
(1) Fisher 變數變換:第一步驟先把相關係數 r 用 Fisher 變數變換轉成一個逼近常態分配的變數 z:
(2) 計算 z 的信賴區間:
(3) 轉回成 r 的信賴區間:

。新版本(v9.1.3之後)
程式碼:
proc corr data=a  fisher (biasadj=no);
  var  x y;
run;


如同前述,要計算相關係數的信賴區間,必須要先經過 Fisher 變數變換,所以語法就是直接在 PROC CORR 後面加上 fisher 選項。如果想要調整一些誤差,可以再在後面打括弧來指定調整方法。本例是不做任何調整,所以寫上 biasadj = no。報表如下:
                 Calculation and Test of Correlations, 95% CI
            With
Obs   Var   Var    NObs     Corr     ZVal        Lcl        Ucl   pValue
 1     x     y       30  0.87982  1.37496   0.760653   0.941620   <.0001


所以本例的相關係數是 0.87982,而其 95% 信賴區間是 (0.790653, 0.941620)。如果想要把信賴區間畫出來,則 ODS 可以輕鬆辦到:
ods html;
ods graphics on;
proc corr data=a noprint plots=scatter(alpha=.2 0.3 );
 var x y;
run;
ods graphics off;
ods html close;


這個程式碼重點在 plots=scatter 選項可以指定信賴區間的長度,如果都不寫的話則預設值就是 95% 信賴區間。而本例則是畫了 80% 和 70% 的信賴區間,所以 alpha 值設定為  0.2 以及 0.3。結果如下所示:

。舊版本(v8 之前)
舊版本的寫法是先用 PROC CORR 算好相關係數矩陣後另存新檔。然後再用這個新檔去跑前述那三個步驟:
proc corr data=a outp=corr; *outp with pearson correlation coefficient;
  var x y;
run;

data corr_ci;
   set corr (rename=(x=corr) drop=y _name_);
   retain n;
   if _type_='N' then n=corr;
   if _type_='CORR'and corr ^= 1;
   fishersz=0.5*(log(1+corr)-log(1-corr)); *Fisher Z transformation;
   sigmaz=1/sqrt(n-3);                     *variance;
   l95=fishersz-1.96*sigmaz;               *α=0.05, i.e. at 95% level;
   u95=fishersz+1.96*sigmaz;
                                                                   
   l95=(exp(2*l95)-1)/(exp(2*l95)+1);      *inverse of Fisher Z ;
   u95=(exp(2*u95)-1)/(exp(2*u95)+1);      *transformation to get CI;
run;


報表如下:
Obs    _TYPE_      corr      n    fishersz     sigmaz      l95        u95
 1      CORR     0.87982    30     1.37496    0.19245    0.76065    0.94162


不過舊版本如何繪圖,這篇技術文件就沒有額外說明了。

CONTACT INFORMATION 
Zaizai Lu
zz_lu@hotmail.com
AstraZeneca Pharmaceuticals
Wilmington, Delaware

沒有留言:

張貼留言

要問問題的人請在文章下方的intensedebate欄位留言,請勿使用blogger預設的意見表單。今後用blogger意見表單留言的人我就不回應了。

CODE { display: block; /* fixes a strange ie margin bug */ font-family: Courier New; font-size: 8pt; overflow:auto; background: #f0f0f0 url(http://klcintw.images.googlepages.com/Code_BG.gif) left top repeat-y; border: 1px solid #ccc; padding: 10px 10px 10px 21px; max-height:200px; height:200px; // for IE6 line-height: 1.2em; }