|
|
楼主
楼主 |
发表于 2012-3-31 16:09:18
|
只看该作者
Computation of Correlation Coefficient and Its Confidence In
From LCChien's blog on blogspot
原文載點：<a href="http://www2.sas.com/proceedings/sugi31/170-31.pdf"><!-- m --><a class="postlink" href="http://www2.sas.com/proceedings/sugi31/170-31.pdf">http://www2.sas.com/proceedings/sugi31/170-31.pdf</a><!-- m --></a><br /><br />在一些特殊的情況下，相關係數需要搭配信賴區間。PROC CORR程序不像其他程序語法一樣在後面的選項加上一個 CL 就可以自動生出信賴區間，而是需要用別的選項語法才能產生。而且，這個功能是只有 SAS v9.1.3 之後的版本才有，舊版本則需要靠使用者自己寫程式去算。David Shen 和 Zaizai Lu 在 SUGI31 發表了一篇技術文件來解說如何在新舊版本下算出相關係數的信賴區間。<br /><a name='more'></a><b><span style="font-size: large;">。理論</span></b><br />由於相關係數並不是常態分配，所以要經過三個步驟，先去做變數變換，算出變換後的信賴區間，然後再轉回去。步驟如下：<br />(1) Fisher 變數變換：第一步驟先把相關係數 r 用 Fisher 變數變換轉成一個逼近常態分配的變數 z：<br /><div style="text-align: center;"><img src="http://latex.codecogs.com/gif.latex?z=0.5\times log_{e}(\frac{1+r}{1-r})" title="z=0.5\times log_{e}(\frac{1+r}{1-r})" /></div>(2) 計算 z 的信賴區間：<br /><div style="text-align: center;"><img src="http://latex.codecogs.com/gif.latex?L_{z}=z-z_{(1-\alpha /2)}\sqrt{\frac{1}{n-3}}" title="L_{z}=z-z_{(1-\alpha /2)}\sqrt{\frac{1}{n-3}}" /></div><div style="text-align: center;"><img src="http://latex.codecogs.com/gif.latex?U_{z}=z+z_{(1-\alpha /2)}\sqrt{\frac{1}{n-3}}" title="U_{z}=z+z_{(1-\alpha /2)}\sqrt{\frac{1}{n-3}}" /></div>(3) 轉回成 r 的信賴區間：<br /><div style="text-align: center;"><img src="http://latex.codecogs.com/gif.latex?L_{r}=tanh(L_{z})=\frac{exp(2L_{z})-1}{exp(2L_{z})+1}" title="L_{r}=tanh(L_{z})=\frac{exp(2L_{z})-1}{exp(2L_{z})+1}" /></div><div style="text-align: center;"><img src="http://latex.codecogs.com/gif.latex?U_{r}=tanh(U_{z})=\frac{exp(2U_{z})-1}{exp(2U_{z})+1}" title="U_{r}=tanh(U_{z})=\frac{exp(2U_{z})-1}{exp(2U_{z})+1}" /></div><br /><b><span style="font-size: large;">。新版本（v9.1.3之後）</span></b><br />程式碼：<br /><code>proc corr data=a <span style="color: red;">fisher (biasadj=no)</span>;<br /> var x y;<br />run;</code><br /><br />如同前述，要計算相關係數的信賴區間，必須要先經過 Fisher 變數變換，所以語法就是直接在 PROC CORR 後面加上 fisher 選項。如果想要調整一些誤差，可以再在後面打括弧來指定調整方法。本例是不做任何調整，所以寫上 biasadj = no。報表如下：<br /><code> Calculation and Test of Correlations, 95% CI<br /> With<br />Obs Var Var NObs Corr ZVal Lcl Ucl pValue<br /> 1 x y 30 0.87982 1.37496 0.760653 0.941620 <.0001</code><br /><br />所以本例的相關係數是 0.87982，而其 95% 信賴區間是 (0.790653, 0.941620)。如果想要把信賴區間畫出來，則 ODS 可以輕鬆辦到：<br /><code>ods html;<br />ods graphics on;<br />proc corr data=a noprint <span style="color: red;">plots=scatter(alpha=.2 0.3 )</span>;<br /> var x y;<br />run;<br />ods graphics off;<br />ods html close;</code><br /><br />這個程式碼重點在 plots=scatter 選項可以指定信賴區間的長度，如果都不寫的話則預設值就是 95% 信賴區間。而本例則是畫了 80% 和 70% 的信賴區間，所以 alpha 值設定為 0.2 以及 0.3。結果如下所示：<br /><div class="separator" style="clear: both; text-align: center;"><a href="http://i.imgur.com/w53vN.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="484" src="http://i.imgur.com/w53vN.png" width="640" /></a></div><br /><b><span style="font-size: large;">。舊版本（v8 之前）</span></b><br />舊版本的寫法是先用 PROC CORR 算好相關係數矩陣後另存新檔。然後再用這個新檔去跑前述那三個步驟：<br /><code>proc corr data=a outp=corr; *outp with pearson correlation coefficient;<br /> var x y;<br />run;<br /><br />data corr_ci;<br /> set corr (rename=(x=corr) drop=y _name_);<br /> retain n;<br /> if _type_='N' then n=corr;<br /> if _type_='CORR'and corr ^= 1;<br /> fishersz=0.5*(log(1+corr)-log(1-corr)); *Fisher Z transformation;<br /> sigmaz=1/sqrt(n-3); *variance;<br /> l95=fishersz-1.96*sigmaz; *α=0.05, i.e. at 95% level; <br /> u95=fishersz+1.96*sigmaz;<br /> <br /> l95=(exp(2*l95)-1)/(exp(2*l95)+1); *inverse of Fisher Z ; <br /> u95=(exp(2*u95)-1)/(exp(2*u95)+1); *transformation to get CI;<br />run;</code><br /><br />報表如下：<br /><code>Obs _TYPE_ corr n fishersz sigmaz l95 u95<br /> 1 CORR 0.87982 30 1.37496 0.19245 0.76065 0.94162</code><br /><br />不過舊版本如何繪圖，這篇技術文件就沒有額外說明了。<br /><br /><b>CONTACT INFORMATION </b><br />Zaizai Lu<br /><!-- e --><a href="mailto:zz_lu@hotmail.com">zz_lu@hotmail.com</a><!-- e --><br />AstraZeneca Pharmaceuticals<br />Wilmington, Delaware<div class="blogger-post-footer"><img width='1' height='1' src='https://blogger.googleusercontent.com/tracker/6268919072942670865-6117377632566650273?l=sugiclub.blogspot.com' alt='' /></div> |
|
