相関係数の標準誤差

母集団から抽出された標本より統計量を得るとき、統計量の標準偏差を標準誤差という。

相関係数の標準誤差としては、相関係数の有意性の検定の文脈で、帰無仮説$\rho=0$のもとでの以下がよく紹介される。

SE_{\rho} = \sqrt{ 1-\rho^2 \over n-2 }　　　　　　(1)

これを確認するためにシミュレーションを行い標準誤差を計算すると、$\rho$が大きくなるに従い上式にはあわなくなる。上式は帰無仮説のもとでの標準誤差なので当然である。なおここでいうシミュレーションとは、乱数で複数セットの標本を用意しそれぞれについて相関係数を計算し標準誤差を計算することをいう。

そこで調べたところ一般には$\rho$の標準誤差は次式になるようである。

SE_{\rho} = {1-\rho^2 \over \sqrt{ n-2 }}　　　　　　(2)

これは以下の二つのwebサイトとその引用文献による。
[1] https://teramonagi.hatenablog.com/entry/20170928/1506611938
[2] https://stats.stackexchange.com/questions/226380/derivation-of-the-standard-error-for-pearsons-correlation-coefficient

シミュレーションすると不偏性はない。とくに標本サイズが小さいときにズレが大きい。このズレを小さくするために分母を$\sqrt{ n-3 }$にしたほうが、よいのではないかという議論もあった。一致性はあるようである。

matlabによる検証

(2)式、シミュレーション、ブートストラップ法により標準誤差をmatlabで計算した例を下記に示す。標本サイズを100とし、シミュレーションでは100セットの標本を用いた。ブートストラップは標本サイズ100のデータを用い1000回のリサンプリングを行った(このあたりの数字は特に調整を行なっておらず意味がない)。

(2)式とシミュレーションの結果はよく一致するが、$\rho$が小さいときにズレが大きい傾向がある。この場合(1)式の方がシミュレーションによく一致する。

ブートストラップ法はただ100個のデータのみを用いてある程度標準誤差を評価できる。乱数生成したデータによりブートストラップの結果はゆらぐ。

母集団の相関係数 rho = 0.300000 0.500000 0.700000 0.900000 
(1-rho^2)/sqrt(n-2) 0.091924 0.075761 0.051518 0.019193 
シミュレーション　     0.099243 0.077898 0.051749 0.020465 
ブートストラップ　　　  0.105727 0.067170 0.042296 0.017960 

N = 100;
L = 100;
std1 = [];
std2 = [];
std3 = [];
for rho = [ 0.3 0.5 0.7 0.9]
  std1(end+1) = (1-rho^2)/sqrt(N-2);
  std2(end+1) = sim( N, L, rho );
  std3(end+1) = boot(N, L, rho );
end

fprintf( '母集団の相関係数 rho = %f %f %f %f \n', 0.3, 0.5, 0.7, 0.9 );
fprintf( '(1-rho^2)/sqrt(n-2) %f %f %f %f \n', std1 );
fprintf( 'シミュレーション　     %f %f %f %f \n', std2 );
fprintf( 'ブートストラップ　　　  %f %f %f %f \n', std3 ); fprintf( '\n' );


function ret = sim( N, L, rho )
cors = zeros( L, 1 );
for i = 1:L
    x = mvnrnd( [0;0], [1 rho; rho 1], N);
    tmp = corr( x );
    cors(i,1) = tmp( 1, 2 );
end
ret = std( cors );
end

function ret = boot( N, L, rho )
   x = mvnrnd( [0;0], [1 rho; rho 1], N);
   [ci, st] = bootci( L*10, @corrx, x ); 
   ret = std( st );
end

function ret = corrx( x )
  tmp = corr( x );
  ret = tmp( 1, 2 );
end