統計学：偏相関係数って何？式の導出から整理します。

変数\(x\)と変数\(y\)に因果関係があり（\(x\)が原因、\(y\)が結果）、変数\(x\)と変数\(z\)の間にも因果関係がある（\(x\)が原因、\(z\)が結果）場合、\(y\)と\(z\)は相関があるような結果になります。しかし、その相関関係の中には\(x\)を介した影響が含まれています。偏相関係数とは、\(x\)の影響を除いた\(y\)と\(z\)の相関関係を表すものです。

偏相関係数の計算式は、以下のように表されるのですが、教科書でも突如この式が出てきたりして、何を意味しているのかわかりませんよね。この式が何を意味するものなのか整理しました。
$$
\rho_{yz\cdot{x}}=\frac{\rho_{yz}-\rho_{xy}\rho_{xz}}{\sqrt{1-\rho_{xy}^2}\sqrt{1-\rho_{xz}^2}}
$$
結論から言うと、\(y,z\)に対し、それぞれ、\(x\)に関する線形回帰モデルを作り、その残差を\(y’,z’\)とする場合、\(y’,z’\)の相関係数が偏相関係数です。具体的に式を計算して確かめてみます。

目次

• 準備（線形回帰モデル）
• 偏相関係数の導出
• Rで試してみる
  • ケース1（y’,z’が独立な場合）
  • ケース2（y’,z’が独立でない場合）

準備（線形回帰モデル）

偏相関係数の計算式の導出の前に、まず、線形回帰モデルの結果について整理しておきます（線形回帰モデルについては、別途記事を作成します）。

データ\((x_1,y_1),(x_2,y_2),…\)に対し、\(x\)を説明変数とした線形回帰モデルを作成すると、説明変数\(x\)の回帰係数\(a\),切片\(b\)、残差\(y’\)は以下のように表せます。

$$
\begin{eqnarray}
a&=&\frac{\rho_{xy}\sigma_y}{\sigma_x}\\
b&=&\mu_y-a\mu_x\\
\\
y’&=&y-(ax+b)\\
&=&y-\mu_y-a(x-\mu_x)\\
\\
&\mu_x&,\mu_y：x,yの平均\\
&\sigma_x&,\sigma_y：x,yの標準偏差\\
&\rho_{xy}&：x,yの相関係数
\end{eqnarray}
$$

また、残差\(y’\)に対しては、\(E(y’)=0\)が成立します。

これらの結果を用いて、偏相関係数の導出を行なっていきます。

偏相関係数の導出

\(y,z\)に対して線形回帰モデルを作成した場合の残差をそれぞれ\(y’,z’\)とし、それぞれの回帰係数を\(a,c\)とすると、\(y’,z’,a,c\)は以下のように表せます。
$$
\begin{eqnarray}
y’&=&y-\mu_y-a(x-\mu_x)\\
z’&=&z-\mu_z-c(x-\mu_x)\\
\\
a&=&\frac{\rho_{xy}\sigma_y}{\sigma_x}\\
c&=&\frac{\rho_{xz}\sigma_z}{\sigma_x}\\
\end{eqnarray}
$$

このとき、\(y’,z’\)の相関係数\(\rho_{y’z’}\)は以下のように表せます。

$$
\begin{eqnarray}
\rho_{y’z’}&=&\frac{E{(y’-\bar{y’})(z’-\bar{z’})}}{\sqrt{E{(y’-\bar{y’})^2}}\sqrt{E{(z’-\bar{z’})^2}}}\\
&=&\frac{E(y’z’)}{\sqrt{E(y’^2)}\sqrt{E(z’^2)}}\quad(∵\bar{y’}=0,\bar{z’}=0)\\
\\
\end{eqnarray}
$$

\(E(y’z’),E(y’^2),E(z’^2)\)は\(\sigma_x,\sigma_y,\sigma_z,\rho_{xy},\rho_{xz},\rho_{yz}\)を用いて整理すると以下のように表せます。

$$
\begin{eqnarray}
E(y’z’)&=&E[\{y-\mu_y-a(x-\mu_x)\}\{z-\mu_z-c(x-\mu_x)\}]\\
&=&E[(y-\mu_y)(z-\mu_z)]-aE[(x-\mu_x)(z-\mu_z)]\\
&& -cE[(x-\mu_x)(y-\mu_y)]+acE[(x-\mu_x)^2]\\
&=&\rho_{yz}\sigma_y\sigma_z-a\rho_{xz}\sigma_x\sigma_z-c\rho_{xy}\sigma_x\sigma_y+ac\sigma_x^2\\
&=&\rho_{yz}\sigma_y\sigma_z-\frac{\rho_{xy}\sigma_y}{\sigma_x}\rho_{xz}\sigma_x\sigma_z-\frac{\rho_{xz}\sigma_z}{\sigma_x}\rho_{xy}\sigma_x\sigma_y+\frac{\rho_{xy}\sigma_y}{\sigma_x}\frac{\rho_{xz}\sigma_z}{\sigma_x}\sigma_x^2\\
&=&\sigma_y\sigma_z(\rho_{yz}-\rho_{xy}\rho_{xz})\\
\\
E(y’^2)&=&E[\{y-\mu_y-a(x-\mu_x)\}^2]\\
&=&E[(y-\mu_y)^2]-2aE[(y-\mu_y)(x-\mu_x)]+a^2E[(x-\mu_x)^2]\\
&=&\sigma_y^2-2a\rho_{xy}\sigma_x\sigma_y+a^2\sigma_x^2\\
&=&\sigma_y^2-2\frac{\rho_{xy}\sigma_y}{\sigma_x}\rho_{xy}\sigma_x\sigma_y+\frac{\rho_{xy}^2\sigma_y^2}{\sigma_x^2}\sigma_x^2\\
&=&\sigma_y^2(1-\rho_{xy}^2)\\
\\
E(z’^2)&=&\sigma_z^2(1-\rho_{xz}^2)\\
\end{eqnarray}
$$

したがって、\(y’,z’\)の相関係数は以下のようになります。これが偏相関係数の正体です。

$$
\begin{eqnarray}
\rho_{y’z’}&=&\frac{\sigma_y\sigma_z(\rho_{yz}-\rho_{xy}\rho_{xz})}{\sqrt{\sigma_y^2(1-\rho_{xy}^2)}\sqrt{\sigma_z^2(1-\rho_{xz}^2)}}\\
&=&\frac{\rho_{yz}-\rho_{xy}\rho_{xz}}{\sqrt{1-\rho_{xy}^2}\sqrt{1-\rho_{xz}^2}}\\
&\equiv&\rho_{yz\cdot{x}}
\end{eqnarray}
$$

Rで試してみる

ケース1（y’,z’が独立な場合）

$$
\begin{eqnarray}
y&=&ax+b+y’\\
z&=&cx+d+z’\\
a&=&1,\quad b=2,\quad c=5,\quad d=10\\
\\
&y’,z’&は独立 \\
y’&\sim& N(0,2^2)\\
z’&\sim& N(0,1^2)\\
\end{eqnarray}
$$

library(ggplot2)
library(gridExtra)

a <- 1
b <- 2
c <- 5
d <- 10

n <- 20  # データ数

set.seed(1)
x <- runif(n, min=0, max=10)
y_dash <- rnorm(n, 0, 2)  # N(0, 2^2)
z_dash <- rnorm(n, 0, 1)  # N(0, 1^2)

y<- a*x + b + y_dash
z<- c*x + d + z_dash

model_y <- lm(y~x)  #線形回帰
model_z <- lm(z~x)  #線形回帰

#回帰モデルの残差計算
y_dash_est <- y - (model_y$coefficients[2] * x + model_y$coefficients[1])
z_dash_est <- z - (model_z$coefficients[2] * x + model_z$coefficients[1])

残差\(y’\)と\(z’\)の相関を計算してみると以下のようになります。

#残差y',z'の相関係数
cor(y_dash_est, z_dash_est)
## [1] -0.06684436

一方、偏相関係数を計算してみると以下のようになります。

#偏相関係数の計算
r_xy <- cor(x,y)
r_xz <- cor(x,z)
r_yz <- cor(y,z)

(r_yz - r_xy*r_xz)/sqrt(1-r_xy^2)/sqrt(1-r_xz^2)  #偏相関係数
## [1] -0.06684436

このように、両者は一致していることが確認できます。

ケース1の\(x-y,x-z,y-z\)の関係は以下のようになっています。

\(y’,z’\)の関係は以下のようになっています（真の値、推定値）。この図からも、\(y’,z’\)の相関（\(y,z\)の偏相関係数）は小さいことがわかります。

ケース2（y’,z’が独立でない場合）

$$
\begin{eqnarray}
y&=&ax+b+y’\\
z&=&cx+d+z’\\
a&=&1,\quad b=2,\quad c=5,\quad d=10\\
\\
&y’,z’&は独立でない \\
z’&=&y’+\epsilon\\
y’&\sim& N(0,2^2)\\
\epsilon &\sim& N(0,1^2)\\
\end{eqnarray}
$$

library(ggplot2)
library(gridExtra)

a <- 1
b <- 2
c <- 5
d <- 10

n <- 20  # データ数

set.seed(1)
x <- runif(n, min=0, max=10)
y_dash <- rnorm(n, 0, 2)  # N(0, 2^2)
eps <- rnorm(n, 0, 1)     # N(0, 1^2)
z_dash <- y_dash + eps

y<- a*x + b + y_dash
z<- c*x + d + z_dash

model_y <- lm(y~x)  #線形回帰
model_z <- lm(z~x)  #線形回帰

#回帰モデルの残差計算
y_dash_est <- y - (model_y$coefficients[2] * x + model_y$coefficients[1])
z_dash_est <- z - (model_z$coefficients[2] * x + model_z$coefficients[1])

残差\(y’\)と\(z’\)の相関を計算してみると以下のようになります。

#残差y',z'の相関係数
cor(y_dash_est, z_dash_est)
## [1] 0.9303714

一方、偏相関係数を計算してみると以下のようになります。

#偏相関係数の計算
r_xy <- cor(x,y)
r_xz <- cor(x,z)
r_yz <- cor(y,z)

(r_yz - r_xy*r_xz)/sqrt(1-r_xy^2)/sqrt(1-r_xz^2)  #偏相関係数
## [1] 0.9303714

こちらの場合も、両者が一致していることが確認できます。

ケース2の\(x-y,x-z,y-z\)の関係は以下のようになっています。

\(y’,z’\)の関係は以下のようになっています（真の値、推定値）。この図からも、\(y’,z’\)の相関（\(y,z\)の偏相関係数）が大きいことがわかります。

R言語とは？Rでできること。現役Rユーザーが語る。

R言語は、データ分析、統計解析に特化した言語であるとよく言われます。しかし「データ分析」「統計解析」しかできないと考えてしまっては損をしています。C/C++のような汎用のプログラミング言語同様の制御構造を備えていますし、CRANで公開されているパッケージを利用すれば、はるかに簡単に望みの機能を実装することができる言語です。

walkintheforest.net

まとめ

教科書等でいきなり現れる偏相関係数について、式の導出からその意味を整理してみました。偏相関係数とは線形回帰したときのの残差通しの相関係数のことです。