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