偏相関係数の話

見せかけの相関に注意しましょう。
データを弄ると必ず聞く言葉です。
偏相関係数を見ましょう。
時々聞く言葉です。というわけで偏相関係数の説明をします。
偏相関係数は、共通の特徴量を持つと思われる量の相関を測るものです。
相関係数が高いのに偏相関係数が低い状態の事を、見せかけの相関がある状態という訳です。
相関係数の復習から始めて、偏相関係数の定義を学んで、時系列分析での応用を説明します。
最後に、pythonで計算させてみます。
参考文献は2冊です。

https://amzn.to/2vhViH3

amzn.to

Amazon.co.jp

amzn.to

記事で使っているソースコードはgithub に置いてあります。
https://github.com/msamunetogetoge

相関係数

相関係数の定義を復習します。
期待値を取る操作を$E[-]$で表すとき¹ 、$x,y$の相関係数$r_{x,y}$は,色々な量を定義して、以下のように計算されます。
$$\begin{array}{rcl}\mathrm{C}\mathrm{o}\mathrm{v}(x,y)& =& E[(x-E[x])(y-E[y])]\\ \mathrm{V}\mathrm{a}\mathrm{r}(x)& =& \mathrm{C}\mathrm{o}\mathrm{v}(x,x)=E[x^2]–{(E[x])}^2\\ r_{x,y}& =& \frac{\mathrm{C}\mathrm{o}\mathrm{v}(x,y)}{\sqrt{\mathrm{V}\mathrm{a}\mathrm{r}(x)\mathrm{V}\mathrm{a}\mathrm{r}(y)}}\end{array}$$
$x\cdot y=\mathrm{C}\mathrm{o}\mathrm{v}(x,y)$を内積だと思うと、相関係数は$x$と$y$のなす角度を測っていると解釈できます。そうすると、$\Vert r_{x,y}\Vert$が1に近いという事は、$x,y$が似ていると解釈できます。
耳にタコが出来るほど聞く話ですが、相関係数が高いからと言って、$x,y$の間に因果関係があるとは限らない事に注意が必要です。
例えば、体重と年収が相関があったとしましょう。体重が重い程、年収が高いようなデータがあったとします。意味不明だなと思って調べてみると、何故かデータ収集の対象年齢が12歳からとなっていることが分かりました。12歳の子供は体重が軽く、年収がほぼ0なので、そんなデータになっていたと分かってめでたしめでたしです。
未知の事象に挑戦している時に、同じような事が起きたときに冷静に行動するための量として、偏相関係数があります。

偏相関係数

偏相関係数は、$x,y,z,z_1,\cdots ,z_n$というデータがあるとき、$x,y$から、$z$たちの影響を除いた相関係数です。この意味を説明します。
$x,y$から、$z$たちの影響を除いた 偏相関係数$r_{x,y,z}$は以下のように定義されます。
$x,y$ を$z$たちから回帰した量を $\hat{x},\hat{y}$とします。
$$\begin{array}{rcl}{y}^{\prime }& =& y-\hat{y}\\ z^{\prime }& =& z-\hat{z}\\ r_{x,y,z}& =& \frac{\mathrm{C}\mathrm{o}\mathrm{v}(x^{\prime },y^{\prime })}{\sqrt{\mathrm{V}\mathrm{a}\mathrm{r}(x^{\prime })\mathrm{V}\mathrm{a}\mathrm{r}(y^{\prime })}}\end{array}$$
変数が大量にある場合でも同じ定義です。$\hat{y}-y$を行う事で、$y$から$x$ の影響を排除していると考えます。
実際、$y^{\prime },x$は無相関である事が示せます。単回帰の場合に計算してみます。
単回帰の時、回帰式は以下のようになっていました。
²
$$\begin{array}{rcl}\hat{y}& =& \widehat{a_y}x+\widehat{b_y}\\ \widehat{a_y}& =& \frac{\mathrm{C}\mathrm{o}\mathrm{v}(x,y)}{\mathrm{V}\mathrm{a}\mathrm{r}(x)}=r_{x,y}\\ \widehat{b_y}& =& E[y]–\widehat{a_y}E[x]\end{array}$$
これらの式を元に共分散を計算します。
$$\begin{array}{rcl}\mathrm{C}\mathrm{o}\mathrm{v}(x,y^{\prime })& =& E[(x-E[x])(y^{\prime }-0)]\\ & =& E[(x-E[x])\{(y-E[y])-r_{x,y}(x-E[x])\}]\\ & =& \mathrm{C}\mathrm{o}\mathrm{v}(x,y)–r_{x,y}E[(x-E[x]{)}^2]\\ & =& \mathrm{C}\mathrm{o}\mathrm{v}(x,y)–\mathrm{C}\mathrm{o}\mathrm{v}(x,y)=0\end{array}$$
$r_{x,y,z}$を計算してみます。
まずは$x,y$が$z$だけに依っていると仮定した場合³です。
$$\begin{array}{rcl}{r}_{x,y,z}& =& \frac{\mathrm{C}\mathrm{o}\mathrm{v}(x^{\prime },y^{\prime })}{\sqrt{\mathrm{V}\mathrm{a}\mathrm{r}(x^{\prime })\mathrm{V}\mathrm{a}\mathrm{r}(y^{\prime })}}\\ & =& \frac{E[x^{\prime }{y}^{\prime }]}{\sqrt{E[x^{\prime 2}]E[y^{\prime 2}]}}\end{array}$$
$E[x^{\prime }{y}^{\prime }],E[x^{\prime 2}]$を計算すると以下のようになります。
$$\begin{array}{rcl}E[x^{\prime }{y}^{\prime }]& =& (r_{x,y}-r_{x,z}{r}_{y,z})\mathrm{V}\mathrm{a}\mathrm{r}(x)\mathrm{V}\mathrm{a}\mathrm{r}(y)\\ E[x^{\prime }]& =& \mathrm{V}\mathrm{a}\mathrm{r}(x{)}^2(1-r_{x,z}^2)\end{array}$$
偏相関係数の式に代入すると、
$$\begin{array}{r}{r}_{x,y,z}=\frac{r_{x,y}-r_{x,z}{r}_{y,z}}{\sqrt{1-r_{x,z}^2}\sqrt{1-r_{y,z}^2}}\end{array}$$
です。
$x,y$が$z_1,\cdots ,z_n$に依る場合も少しだけ計算しておきます⁴ 。回帰式を$$\begin{array}{r}\hat{y}={\beta }_0+\sum _{i=1}^n{\beta }_i{z}_i\end{array}$$
と書く時、$\beta$たちは、
$$\begin{array}{r}\beta =(z^{T}z{)}^{-1}(z^{T}y)\end{array}$$
と書けます。ただし、$z$は$(1,z_1,\cdots ,z_n)$と$1$と$z_i$たちを横に並べて得られる (データ数)$\times (n+1)$ 行列です。また、以下のようにも書けます。
$z$の分散共分散行列$S$と、$y$と $z_i$たちの共分散$S_y$は、以下のように計算出来ます。
$$\begin{array}{rcl}{S}_{i,j}& =& \mathrm{C}\mathrm{o}\mathrm{v}(z_i,z_j)\\ {S_y}_i& =& \mathrm{C}\mathrm{o}\mathrm{v}(y,z_i)\end{array}$$
この行列と、 ${\beta }_0=-E[y]+\sum {\beta }_{i}E[z_i]$で、 ${\beta }_1,\cdots ,{\beta }_n$は、
$$\begin{array}{r}\beta =S^{-1}{S}_y\end{array}$$
となります。行列とベクトルで書くと、回帰式は
$$\begin{array}{r}\hat{y}=z{S}^{-1}{S}_y\end{array}$$
となります。偏相関係数は、この回帰式を使って
$$\begin{array}{r}{r}_{x,y,z_1,\cdots ,z_n}=\frac{\mathrm{C}\mathrm{o}\mathrm{v}(x^{\prime },y^{\prime })}{\sqrt{\mathrm{V}\mathrm{a}\mathrm{r}(x^{\prime })\mathrm{V}\mathrm{a}\mathrm{r}(y^{\prime })}}\end{array}$$
と書かれます。手計算するのは大変ですが、データがある場合にパソコンに計算させるのは簡単です。

時系列分析での応用

偏相関係数は時系列分析で重要な役割を果たします。具体的には、定常なAR過程か否かを見極めるのに使う事が出来ます。つまり、定常な$AR(p)$過程 なら、 ある次数から$p+1$から偏相関係数が0 となります。⁵
簡単の為に、$AR(1)$過程の2次の偏相関係数を計算してみます。$\{y_t\}$を定常な $AR(1)$過程とします。
$$\begin{array}{r}{y}_t=\varphi y_{t-1}+{ϵ}_t\end{array}$$
ただし、$\{{ϵ}_t\}$は分散0のホワイトノイズとします。
2次の偏相関係数を計算するために、$y_{t-1}$による$y_t,y_{t-2}$ の回帰式を考えます。⁶
$$\begin{array}{r}{\hat{y}}_{t-2}=\beta y_{t-1}\end{array}$$
最小二乗法から、$\beta =\varphi$が分かります。$y_t$についても同様です。このことから、
$$\begin{array}{rcl}{r}_2& =& \mathrm{C}\mathrm{o}\mathrm{v}(y_{t-2}–{\hat{y}}_{t-2},y_t–{\hat{y}}_t)\\ & =& \mathrm{C}\mathrm{o}\mathrm{v}(y_{t-2}–\varphi y_{t-1},\varphi y_{t-1}+{ϵ}_t–\varphi y_{t-1})\\ & =& \mathrm{C}\mathrm{o}\mathrm{v}(y_{t-2}–\varphi y_{t-1},{ϵ}_t)=0\end{array}$$
同じような計算で、3次以上でも0になることが分かります。

pythonでの計算

初めに、irisデータで偏相関係数を計算してみます。
patal length とpatal width の相関が非常に高いので、他の変数の影響を取り除いて偏相関係数を計算してみます。

from sklearn.datasets import load_iris
import pandas as pd
import numpy as np
iris = load_iris()

df = pd.DataFrame(iris.data, columns=iris.feature_names)

X=df.iloc[:,:2]
X["beta_0"]=1
X=X[["beta_0","sepal length (cm)","sepal width (cm)" ]]

#patal length とpatal width の相関が0.96なので、偏相関係数を計算
y=df[["petal length (cm)","petal width (cm)"]]
beta = np.dot(np.dot(np.linalg.inv( np.dot(X.T , X) ) , X.T), y)
#重回帰による予測値
y_hat = np.dot(X, beta)
#patal length とpatal width の相関が0.96
Y = y - y_hat
S_Y =np.zeros((2,2))
for j in [0,1]:
    Y_j =Y.iloc[:,j] - np.mean(Y.iloc[:,j])
    for i in [0,1] :
        Y_i =Y.iloc[:,i] - np.mean(Y.iloc[:,i])  
        S_Y[j,i]=np.dot(Y_i, Y_j)

pac = S_Y[0,1]/np.sqrt(S_Y[0,0]*S_Y[1,1])
print(pac)
#0.870769773836145

偏相関係数は0.87 と大きな値になりました。
次に、適当に$AR(1)$過程を生成して相関係数と偏相関係数を計算してみます。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

#相関係数計算の関数
def autcor(X, k):
    X_bar = np.mean(X)
    X_bar_k = X_bar * np.ones(k)
    X_k = np.append(X_bar_k ,X)
    k_X = np.append(X, X_bar_k )
    r= np.dot(X_k -X_bar , k_X - X_bar )
    r/= np.linalg.norm(X-X_bar)**2
    
    return r
#偏相関係数計算の関数
def Pac(Y,k):
    Y = pd.DataFrame(Y, columns=["y_t"])
    for i in range(k):
        if i ==0:
            Y["y_t-1"] = Y["y_t"].shift().fillna(0)
        else:    
            Y["y_t-"+str(i+1)]=Y["y_t-"+str(i)].shift().fillna(0)
    X =Y.iloc[:,1:k-1]
    y=Y.iloc[:,[0,k]]
    beta = np.dot(np.dot(np.linalg.inv( np.dot(X.T , X) ) , X.T), y)
    y_hat = np.dot(X, beta)
    y_til = y - y_hat
    S_Y =np.zeros((2,2))
    for j in [0,1]:
        Y_j =y_til.iloc[:,j] - np.mean(y_til.iloc[:,j])
        for i in [0,1] :
            Y_i =y_til.iloc[:,i] - np.mean(y_til.iloc[:,i])  
            S_Y[j,i]=np.dot(Y_i, Y_j)
    pac = S_Y[0,1]/np.sqrt(S_Y[0,0]*S_Y[1,1])
    return pac

 #AR(1)の生成
phi = 0.6
epsilon = np.random.randn(200000)
T=1000
Y=np.zeros(T)
for t in range(T):
    Y[t]=phi*Y[t-1] +epsilon[t]   

#10次まで計算
n=10
pac=[]
ac=[]
for k in range(n):
    pac=np.append(pac,Pac(Y,k+1))
    ac = np.append(ac,autcor(Y,k+1))

AR過程では、指数関数的に相関係数が小さくなっていく様子が見えます。また、偏相関係数は、2次以降でほぼ0となっています。MA過程では次のようなグラフになります。

上の2つの図から分かるように、偏相関係数は、AR過程か否かを判定する材料には出来ますが、MA過程にはあまり効果がありません。

まとめ

• 相関係数について復習した
• 偏相関係数の定義を確認した
• 時系列分析における偏相関係数の役割を解説した
• python で計算してみた

1. $x$が長さ$n$の、 データからなるベクトルだったら、$E[x]=\sum _{i=1}^n{x}_i/n$です。$p(x)$に従う確率変数なら$E[x]=\int xp(x)dx$です
2. 計算は単回帰分析の記事を読んでください 。
3. 単回帰の場合です。
4. 省略されている計算は重回帰分析の記事に書いてます。
5. 定常な$MA(q)$過程 なら、 相関係数が次数$q+1$から0でした。MA過程に関する記事を読んでみて下さい。
6. 一般にk次の偏相関係数を計算するための回帰係数について、${\beta }_i={\beta }_{k-i}$が成り立ちます。