偏相関係数の意味
先日、ある先生が偏相関係数について説明していたけど、適切とはお思えない方法で(少なくとも非常にわかりにくい方法で)、その数式の意味を教えていたので、偏相関係数について改めて考える。
偏相関係数とは、変数xと変数y間の関係において、変数z影響を取り除いた、変数xと変数yの相関係数のことである。
その式は以下の通りになる。
なぜ、この式で変数zの影響を取り除いた変数xと変数yの相関係数が導きだせるのかを説明する上でもっともわかりやすいのは、以下の説明であると思う。
すなわち、
「xを従属変数、zを独立変数とした時の回帰分析において、zでは説明できない部分、すなわち誤差(以下、e1)と、yを従属変数、zを独立変数とした時の回帰分析において、zでは説明できない部分、すなわち誤差(以下、e2)との相関が、変数z影響を取り除いた変数xと変数yの偏相関係数である」という説明である。
直感的にもわかりやすいと思う。
実際に、Rで試してみる。
結果を先取りして言うと、全く同じにはならなかった。
おそらく、誤差e1や誤差e2の平均が0になるとか回帰分析におけるいくつかの仮定を完全に満たしていない事や、そもそも計算の過程における四捨五入などが影響している可能性がある。
そこで、同じ工程を1000回繰り返すことで、どれくらい異なるか調べる。
その工程を簡単に述べると以下の通りである。
①変数zを正規分布からランダムに導出(n=1000)
②変数zと0.3の相関関係にある変数xを導出(n=1000)
③変数zと0.3の相関関係にある変数yを導出(n=1000)
④zをxに回帰させた時の誤差e1を取り出す
⑤zをyに回帰させた時の誤差e2を取り出す
⑥誤差e1と誤差e2の相関係数を求める(=res1)
⑦上記の定義式から、変数zの影響を取り除いた変数xと変数yの偏相関係数を求める(res=2)
まず、res1とres2がどのような関係にあるのかプロットしてみた。
その結果は以下の通り。
当然だが、綺麗な相関関係にある。
要するに、res1とres2は一定の誤差(res1とres2の差)を持って、ほとんど同じ動きをする。
すなわち、一定の誤差をもって「xを従属変数、zを独立変数とした時の回帰分析において、zでは説明できない部分、すなわち誤差(以下、e1)と、yを従属変数、zを独立変数とした時の回帰分析において、zでは説明できない部分、すなわち誤差(以下、e2)との相関が、変数z影響を取り除いた変数xと変数yの偏相関係数」という説明が今回の条件のもとではみられる。
それでは、次にどのくらいそこに誤差がみられたのかを確認しよう。
その誤差(res1とres2の差)のヒストグラムが以下の通り。
これをみると、その誤差の多くは-0.01から0.01の間に収まっている。そこまで、大きな誤差ではないと言える(もちろん分野によって異なるが)。
以上から、数学的には「xを従属変数、zを独立変数とした時の回帰分析において、zでは説明できない部分、すなわち誤差(以下、e1)と、yを従属変数、zを独立変数とした時の回帰分析において、zでは説明できない部分、すなわち誤差(以下、e2)との相関が、変数z影響を取り除いた変数xと変数yの偏相関係数」という説明は全く正しい一方で、実際にRで実験したところ、少なくとも今回の条件下では0.01程度以下の誤差を含む(多くの場合)ことがわかった。
以上のコードは以下の通り。
N<-1000 #試行回数
n<-1000 #サンプル数
res1<-numeric(N) #収納
res2<-numeric(N) #収納
for(i in 1:N){
z<-rnorm(n) #変数z
r<-0.3
k<-rnorm(n)
x<-z*r+sqrt(1-r^2)*k #変数x
u<-rnorm(n)
y<-z*r+sqrt(1-r^2)*u #変数y
lm1<-lm(x~z)
beta1<-lm1$coefficient[2]
e1<-x-(beta1*z) #zをxに回帰させた時の誤差
lm2<-lm(y~z)
beta2<-lm2$coefficient[2]
e2<-y-(beta2*z) #zをyに回帰させた時の誤差
r_xy<-cor(x,y)
r_xz<-cor(x,z)
r_yz<-cor(y,z)
res1[i]<-(r_xy-(r_xz*r_yz))/(sqrt(1-r_xz^2)*sqrt(1-r_yz*2)) #偏回帰係数の定義式から
res2[i]<-cor(e1,e2) #e1とe2の相関
}
#図示
library(ggplot2)
g<-ggplot()
g<-g+geom_point(mapping=aes(x=res1,y=res2))
print(g)
res<-res1-res2
g<-ggplot()
g<-g+geom_histogram(mapping=aes(x=res),binwidth=0.001)
print(g)