はしくれエンジニアもどきのメモ

情報系技術・哲学・デザインなどの勉強メモ・備忘録です。

トップ > 確率・統計 > pythonで偏相関係数行列(pcor)を計算

pythonで偏相関係数行列(pcor)を計算

確率・統計 math 数式

pythonで偏相関係数行列(pcor)を計算

以前，3変数(X,Y,Z)の場合の偏相関係数の式を導出した．

cartman0.hatenablog.com

今回は３変数以上の多変数の場合の偏相関係数行列を求めてみる． ただscipy.stats, pandas, statsmodels,scikit-learnなどで関数がないようなので， 実装して計算してみる．

Rではcor2pcorという相関係数行列を偏相関係数行列に変換してくれる関数があるぽいので， これに近いものを実装してみる．

偏相関係数と似たL1ノルムを利用したスパースな精度行列なるものであればscikit-learnのほうにGraphicalLassoが実装されている．

自分用ノート： https://github.com/Cartman0/MultivariateAnalysis/blob/master/PartialCorrelationCoefficient_%E5%81%8F%E7%9B%B8%E9%96%A2%E4%BF%82%E6%95%B0%E3%81%AE%E5%B0%8E%E5%87%BA.ipynb

３変数の場合の偏相関係数（おさらい）

詳細はこちら．

cartman0.hatenablog.com

$X,Y,Z$の３変数でZの成分を取り除いたX,Yの偏相関係数とは， Zを説明変数，Xを目的変数とする単回帰モデルの予測値$\tilde{X}$はZ成分で説明されるXとなる． 偏相関係数では，Z成分を取り除いたXをその残差$e_{X} = X-\tilde{X}$としている． $Y$も同様で$e_{Y} = Y-\tilde{Y}$．これらの残差はZ成分が取り除かれているのでこれらの残差の相関を偏相関係数としている．

\begin{eqnarray} \rho_{XY \cdot Z} &=& \frac{ \rho_{XY} - \rho_{XZ}\rho_{YZ} }{\sqrt{1 - \rho_{XZ}^{2}} \sqrt{1 - \rho_{YZ}^{2} } }\\\\ \end{eqnarray}

偏相関係数行列

先程の例は，Z成分を除いた偏相関係数1つであったが， 今回は，X成分を除いた，Y成分を除いた偏相関係数も考えた行列を求めてみる．

また，X,Y,Z,Wなどの4次元以上の場合でも考え方は同じ． Z,WでXを表す重回帰モデルから残差を考えればZ,W成分を除いたXとなりYも同様で， 残差の相関はZ,W成分を除いた$\rho_{XY \cdot ZW}$となる．

相関係数行列の逆行列

結論から言うと，相関係数行列の逆行列の各要素を対応する対角成分で正規化すると求まる．（イメージは共分散行列から相関係数行列を求めるときと同じ） 3x3で計算すれば上の式も出てくる．

偏相関係数行列は，相関係数行列の逆行列を利用するので精度行列という言い方のほうが正しいのかもしれない．

$ m \times m $の相関係数行列の逆行列の各要素を以下の記号で定めておく．

\mathbf{\Pi}^{-1} = \begin{pmatrix} \rho^{11} & \rho^{12} & \cdots & \rho^{1m} \\ \rho^{21} & \rho^{22} & \cdots & \rho^{2m} \\ \vdots & \vdots & \ddots & \vdots \\ \rho^{m1} & \rho^{m2} & \cdots & \rho^{mm} \\ \end{pmatrix}

３次元の場合で求めてみる．

3次正方行列Aとする．

\mathbf{A} = \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \\ \end{pmatrix}

その逆行列は，以下で求まる．

\mathbf{A}^{-1} = \frac{1}{|\mathbf{A}|} \tilde{\mathbf{A}}

より，

\begin{eqnarray} |\mathbf{A}| = a_{11}a_{22}a_{33} + a_{12}a_{23}a_{31} + a_{13}a_{21}a_{32} - a_{31}a_{22}a_{13} - a_{21}a_{12}a_{33} - a_{11}a_{32}a_{23} \end{eqnarray}
\mathbf{\tilde{A}}_{ij} = (-1)^{i+j} |M_{ji}|
\begin{eqnarray} \mathbf{\tilde{A}} = \begin{pmatrix} a_{22}a_{33} - a_{23}a_{32} &-(a_{12}a_{33}-a_{13}a_{32})& a_{12}a_{23}-a_{13}a_{22}\\ -(a_{21}a_{33} - a_{23}a_{31})&a_{11}a_{33}-a_{13}a_{31} &-(a_{11}a_{23}-a_{13}a_{21})\\ a_{21}a_{32}-a_{22}a_{31} &-(a_{11}a_{32}-a_{12}a_{31})& a_{11}a_{22}-a_{12}a_{21}\\ \end{pmatrix} \end{eqnarray}

以上より，相関係数行列の逆行列を求める．

相関係数行列を次とする．

\begin{eqnarray*} \mathbf{\Pi}_{3} &=& \begin{pmatrix} 1 & \rho_{XY} & \rho_{XZ} \\ \rho_{YX}& 1 & \rho_{YZ} \\ \rho_{ZX}& \rho_{ZY} & 1 \\ \end{pmatrix}\\\\ \end{eqnarray*}
\begin{eqnarray*} |\mathbf{\Pi}_{3}| &=& 1 + \rho_{XY}\rho_{YZ}\rho_{ZX} + \rho_{XZ}\rho_{YX}\rho_{ZY} -\rho_{XY}^{2} -\rho_{XZ}^{2} - \rho_{YZ}^{2}\\ &=& 1 + 2\rho_{XY}\rho_{YZ}\rho_{ZX} -\rho_{XY}^{2} -\rho_{XZ}^{2} - \rho_{YZ}^{2} \end{eqnarray*}

逆行列は

\begin{eqnarray} \mathbf{\Pi}_{3}^{-1} &=& \begin{pmatrix} 1 & \rho_{XY} & \rho_{XZ} \\ \rho_{YX}& 1 & \rho_{YZ} \\ \rho_{ZX}& \rho_{ZY} & 1 \\ \end{pmatrix}^{-1}\\\\ &=& \frac{1}{|\mathbf{\Pi}_{3}|} \begin{pmatrix} 1-\rho_{YZ}^{2} & \rho_{XZ}\rho_{YZ}-\rho_{XY} & \rho_{XY}\rho_{YZ}- \rho_{XZ}\\ \rho_{YZ}\rho_{XZ}-\rho_{XY} & 1 - \rho_{XZ}^{2} & \rho_{XZ}\rho_{XY} - \rho_{YZ} \\ \rho_{XY}\rho_{YZ}- \rho_{XZ} & \rho_{XZ}\rho_{XY} - \rho_{YZ} & 1 - \rho_{XY}^{2} \\ \end{pmatrix}\\ \end{eqnarray}

例

次の相関係数行列の逆行列をpythonで求めてみる．

\begin{eqnarray} \mathbf{R} = \begin{pmatrix} 1. & 0.6 & 0.6 & 0.6 \\ 0.6 & 1. & 0.36 & 0.36 \\ 0.6 & 0.36 & 1. & 0.36 \\ 0.6 & 0.36 & 0.36 & 1. \end{pmatrix} \end{eqnarray} 
import scipy as sp
from scipy import linalg

R = sp.matrix([[1, 0.6, 0.6, 0.6],
              [0.6, 1, 0.36, 0.36],
              [0.6, 0.36, 1, 0.36],
              [0.6, 0.36, 0.36, 1]])
R

> matrix([[1.  , 0.6 , 0.6 , 0.6 ],
        [0.6 , 1.  , 0.36, 0.36],
        [0.6 , 0.36, 1.  , 0.36],
        [0.6 , 0.36, 0.36, 1.  ]])

行列を可視化する．

import matplotlib.pyplot as plt

%matplotlib inline

plt.matshow(R, vmin=-1, vmax=1, cmap="jet")
plt.colorbar()
plt.show()

相関係数行列

linalg.inv(R)

array([[ 2.6875, -0.9375, -0.9375, -0.9375],
       [-0.9375,  1.5625,  0.    , -0.    ],
       [-0.9375,  0.    ,  1.5625, -0.    ],
       [-0.9375,  0.    ,  0.    ,  1.5625]])

偏相関係数行列を求める

偏相関係数行列の各要素は以下になる．

\rho_{ij \cdot rest} = - \frac{ \rho^{ij} }{ \sqrt{\rho^{ii}\rho^{jj}} }

pythonで計算してみる．

以下を参考にした．

https://digitalcommons.wpi.edu/cgi/viewcontent.cgi?article=6311&context=iqp-all

import scipy as sp
from scipy import linalg

def cor2pcor(R):
    inv_cor = linalg.inv(R)
    rows = inv_cor.shape[0]
    regu_1 = 1 / sp.sqrt(sp.diag(inv_cor))
    regu_2 = sp.repeat(regu_1, rows).reshape(rows, rows)
    pcor = (-inv_cor) * regu_1 * regu_2
    sp.fill_diagonal(pcor, 1)
    return pcor

pcor = cor2pcor(R)
pcor

> array([[ 1.        ,  0.45749571,  0.45749571,  0.45749571],
       [ 0.45749571,  1.        , -0.        ,  0.        ],
       [ 0.45749571, -0.        ,  1.        ,  0.        ],
       [ 0.45749571, -0.        , -0.        ,  1.        ]])

可視化

import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

sns.heatmap(pcor, annot=True, square=True, vmin=-1, vmax=1, fmt=".2f", cmap="jet")
plt.savefig("pcor.png")
plt.show()

偏相関係数行列

正規分布を仮定しているなら相関が0と独立は同値より他の変数が存在するときの条件付き独立として扱える．

  • 変数1,4があるなら変数2と3は独立，
  • 変数1,3があるなら変数2と4は独立
  • 変数1,2があるなら変数3と4は独立

この例だと変数1はすべての変数と関係があるが，他の変数2,3,4同士は独立といえる．