楕円状の確率分布、上から見るか、斜めから見るか -- 回帰直線とp値と相関係数のお話し --

2次元に散布図を描いで、回帰直線とp値と相関係数を求めて関係性を調べる解析は、しばしば行われます。 結果を見て、「p値は 10-3 のように有意な値が出たのに、相関係数があまり良くない」なんてことが良くありませんでしょうか。そこで、実験的に、楕人口的な傾きを持った楕円状の分布を使って、サンプル数、p値、相関係数 R の関係を調べてみることにします。

分野にもよるでしょうが、観測値は正規分布に近いけれど離れた観測値も観察されることが、現実には多いと考えます。こんな時は、t分布の方が、現実にはあっているかもしれません。以前、当ブログで、正規分布とt分布の比較もしていますので、御参照下さい。今回は、t分布に従う分布を乱数として使うことにします。

まずは、確率密度分布に従う傾いた楕円状の分布の観測点を得ます。方法は、

  1. 200個など数の (x, y) を t分布に従う乱数から作ります
  2. y軸方向を何通りかにグシャっと圧縮して、X軸上に近づけます
  3. 原点を中心に回転させます
以上です。

この方法だと、回帰直線の傾きは、回転させた角度 (正確には傾きは、回転角の正接 [タンジェント, tangent])になります。 そこで、回帰直線の値は、あえて求めたり、議論はしないことにします。回帰直線の切片は複雑になるので、0にします。

例によって、julia で計算していきます。プログラムに興味がない方は、コードは読み飛ばしても大丈夫です。


using Statistics, StatsBase, Distributions, Plots, DataFrames, GLM, ColorTypes

# 回転行列
function rotation_m(θ::Float64, A0::Array)
  A1 = zeros(size(A0));
  for i in 1:size(A0,1)
    A1[i,:] = [cos(θ) -sin(θ); sin(θ) cos(θ)]*A0[i,:]
  end
  return A1
end

# x, y ごとに自由度を設定したt分布の(x、y) のサンプル N 個を得る
function get_A(xdf::Int64, ydf::Int64, N::Int64, ratio::Float64)
  x0 = rand(TDist(xdf),N,1);
  y0 = rand(TDist(ydf),N,1)*ratio;
  A0 = hcat(x0,y0)
  return A0
end

xdf = 10 ; #x,yの自由度を設定する
ydf = 10 ; #x,yの自由度を設定する
ratio = 0.5 # y軸方向に潰す割り合い
A0 = get_A(xdf, ydf, 200, ratio); # 潰す前の観測点

A1 = rotation_m(π/6 ,A0); #  回転行列を掛けます

これで、A1 に 200 個の観測点が格納されます。

参考までに、200 個の観測点を求めた場合、回転前にY軸方向に圧縮する Flat rate 別に標本の分布は、図1 のよううになります。乱数ですから、当然、再現性は限定的です。

図1. 30度回転前にY軸方向に圧縮する Flat rate 別に描画した標本の分布 (n = 200)

主観で「相関あり」と感じるのは、Flat rate = 0.5 つまり、Y軸方向を回転前に半分に圧縮した辺り位の人が、多いのではないでしょうか。

次に、p値と相関係数を求めます。この場合、 


data = DataFrame(x=A1[:,1], y=A1[:,2]) # DataFrameを経由します
ols = lm(@formula(y ~ x), data) # 傾きと切片の情報が得られます
pvalue_i = coeftable(ols).cols[4][2] # p値は、これで、取り出せます

cor_xy =cor(A1[:,1], A1[:,2]) # 相関係数は Statistics.cor を使います

こんな感じで求められます。

では、観測点を200個、回転前のY軸方向の圧縮率 Flat rate を変化させ、回転すなわち傾きを 30度 (π/6) に固定した時の、相関係数を求めてみましょう。それぞれの Flat rate で20回繰替えして、相関係数を求めます。


function get_cor()
  A_rep = zeros(200, 2) # 値を格納する空のアレイを作っておく (vcat() で継ぐと遅い))
  count_n =0 # 値を格納する際の行番号を示す
    for ratio_i in 0.1:0.1:1  # ratio は、for ループで変更
      for j in 1:20 # 20回、くりかえし観測
        count_n +=1 # 行番号を返す
        xdf,ydf = 10,10; # 自由度が10なので、かなりバラつく
        ratio = ratio_i
        A0 = get_A(xdf, ydf, 200, ratio); # ratio は、for ループで変更
        A1 = rotation_m(π/6 ,A0); # 回転は π/6 で固定
        cor_xy =cor(A1[:,1], A1[:,2])
        A_rep[count_n,:] = [ratio, cor_xy]
        # println(count_n) # モニター
      end
    end
  return A_rep
end

A_cor = get_cor(); # A_cor に返り値のアレイを格納する

scatter(A_cor[:,1].+(randn(size(A_cor,1)).*3 .*10^-3) ,A_cor[:,2], msize=3, legend=false,
ylabel = "R", xlabel="Flat rate"
) # Plots.jl でグラフを描く. 1回目は恐しく遅い

これで、図1のような分布で、どの位の相関係数が出るかの目安が得られます。(図2)

図2. 30度回転前のY軸方向の圧縮率 Flat rate と相関係数 R の関係 (n = 200)

相関係数 R が 0.8 を越えると、かなり強い相関といわれます。かなり Y軸方向を圧縮しないと R = 0.8 などという値にはなりません。 Flat rate 0.5 のあたりで、R の平均が 0.5 位なのは、多分、偶然です。再現性には検証が必要です。30度の傾きを前提としていますが、Flat rate が 0.8 位だと、かなり厳しい R の値が帰ってきます。

次に、サンプル数によって、どのくらい p値や相関係数 R が変化するか調べてみます。


function get_pval(N::Int64) # N はサンプル数
  A_rep0 = zeros(480,3) # 角度、p値、相関係数を順に格納するアレイ
  countn = 0
  for i in 1:24
    for j in 1:20
      countn +=1
      A1 = get_A_rotated(10, 10, N , 0.4 ,π/48 * i) #自由度は10で、Y軸方向は 0.4  倍に圧縮
      data = DataFrame(x=A1[:,1], y=A1[:,2]);
      ols = lm(@formula(y ~ x), data); # model, GLM.lm ie. linear model
      pvalue_i = coeftable(ols).cols[4][2] # p value
      A_rep0[countn,:] = [π/24 * i, pvalue_i, cor(A1[:,1], A1[:,2])]
    end
  end
  return A_rep0
end

回転前のY軸の圧縮は、0.5 と 0.4 の2つで調べてみました。

図3. 4つの標本数での回転量とp値と相関係数の関係 (回転前のY 軸圧縮率は、0.4[左] と 0.5[右])

相関係数 R は、色で値を示しています。

いずれの場合も、角度が45度付近で、p値は低く、相関係数 R は高くなります。P値は特に、サンプル数の影響が大きくでるようです。

以上、人口的な傾きを持った楕円状の分布を使って、サンプル数、p値、相関係数 R の関係を調べてみました。

直線回帰の分析の場合、自分の持論は、サンプル数、p値, 相関係数に加えて、回帰直線の傾きの値のすべてに加えて、散布図まで見て、慎重に解釈をするという、とてつもなく、面倒かつ慎重な意見です。分析に意味が、どのくらいあるのかは、あくまで、自分で考えるべきだと、考えています。

いくつか、計算をしてみましたが、直線回帰分析は、思っている以上に解釈が難しい解析だと思いました。

直線回帰の場合、相関係数 R に比べると、傾きの p = 0.05 という有意水準は、少し甘めなのかもしれません。 p値の帰無仮説は、「傾きが 0 である」というものですから、0.00001 でも 0 から傾くなら、検出しないといけません。 この意味では、役割を果していると言うべきだろうと考えます。 また、p値はいつも標本数と密接な関係にありますが、図3を見ると、直線回帰の際は、かなり標本数が支配的に感じられました。 やはり、p値だけで判断せずに、どのくらい回帰直線が傾いているかは確認する方が良いでしょう

今回用いた t 分布では、傾きは 45度付近が最も統計学的な指標が強く出ました。 傾き 45度で統計指標が良いかは、t分布と異なる分布や、平均の異なる分布の混合などでも確認するべきかもしれません。 しかし、かなり一般性が高い可能性があります。 解析の前処置として、正規化などで、x軸とy軸のスケールをしっかり揃えてから解析しないと、有意な変化を検出し損ねる可能性があるようです。 多変量解析でも、単純に p値 だけを気にする場合は、各因子のスケーリングは重要な可能性があります。

おしまい

 B! LINE