偏相関係数の概要と求め方（python）

偏相関係数とは

ある変数（$x_1,x_2,\dots, x_n$）がある時、他の変数全てを固定した時の $x_i$ と $x_j$ の相関係数を $x_i$ と $x_j$ の 偏相関係数という。

なぜ偏相関係数を計算する必要があるのか？

相関分析において着目した2つの変数間に相関がある（相関係数の絶対値が大きい）場合、実際2つの変数間に相関がある場合もあるが、「第3の因子」の影響によって相関係数が高くなっているだけの場合も考えられる。

例えば、各都道府県ごとの単位面積あたりのスーパーの数とコンビニの数の相関係数が高くなったする。しかし、この2つの数には直接的には関係あるとは考えにくい。実はこの場合、人口密度という因子が「第3の因子」として、スーパーおよびコンビニの数に影響しているため、見かけ上、スーパーの数とコンビニの数の間に相関があるように見えてしまう。
（このような相関を疑似相関といい、第3の因子を交絡因子という）

このような交絡因子の影響を除いた上での相関係数が偏相関係数あり、二つの変数間の相関関係を正しく見積もることができる。

偏相関係数の計算方法

交絡因子含めた3つの変数から求める

3つの因子をそれぞれ、$x$, $y$, $z$ とおき、それぞれの因子間の相関係数を以下のようにおく。

$\rho_{xy}, \rho_{yz}, \rho_{zx}$

この時、$z$ 因子の影響を除いた $x$ と $y$ の偏相関係数 $\rho_{xy\cdot z}$ は、以下の式から求められる。

$$ \rho_{xy\cdot z} = \frac{\rho_{xy} - \rho_{yz}\rho_{zx}}{\sqrt{1-\rho_{zx}^2}\sqrt{1-\rho_{yz}^2}} $$

愚直に求める

例えば（$x_1,x_2,\dots, x_n$）が多変量正規分布に従う時。

$x_1$ と $x_2$ の偏相関係数を求めるには、$x_1$と$x_2$以外の成分を使って、$x_1$と$x_2$を目的変数とした回帰分析を行った際の、残差の相関係数から求められる。

回帰分析の残差は、他の成分で説明できない（影響がない）成分であるため。
（$x_3, \dots, x_n$成分の影響を除いた$x_1, x_2$独自の成分と言える）

$$x_1 = a_{13}x_3 + \cdots + a_{1n}x_n + b_1 + \varepsilon_1$$

$$x_2 = a_{23}x_3 + \cdots + a_{2n}x_n + b_2 + \varepsilon_2$$

$$ \rho_{1,2,rest} = cor(\varepsilon_1, \varepsilon_2)$$

より簡単な方法

$\rho_{i,j,rest}$ は、$x_1,x_2,\dots, x_n$の共分散行列の逆行列から簡単に計算できる。
（AIcia Solid曰く、この証明がすごく面白いらしい)
（下記コードの「愚直に計算（モジュール不使用）」はこの方法）

偏相関の検定

標本数を$n$、固定する変数を$q$、変数$a,b$の偏相関係を$\rho_{i,j,rest}$とする時、
以下の検定料は、自由度が$n-q-2$ の$t$分布に従う。
棄却される時偏相関係数は0でないとみなす

$$ t_0 = \frac{|\rho_{i,j,rest}|\sqrt{n-q-2}}{\sqrt{1-\rho_{i,j,rest}^2}}$$

コード

参考動画内でも使用している、AlciaSolidProjectチャンネルにおける動画の情報（動画時間、視聴回数など）を使用した。（githubから取得）

import os
import numpy as np
import pandas as pd
from scipy.stats import chi2

pd.options.display.float_format = '{:.4f}'.format

df = pd.read_csv('AIcia_videos.csv')
df['動画時間_s'] = pd.to_timedelta(df['動画時間']).apply(lambda x: x.seconds)
df = df[['動画時間_s', '視聴回数', 'コメント', '高評価件数', '低評価件数']]
df.head()

愚直に計算（モジュール不使用）

# 相関係数を計算
corr_matrix = df.corr()

# 逆行列を計算
corr_matrix_inv = np.linalg.inv(corr_matrix)

# 標準化するための対角行列を計算
diag = np.diag(1/np.sqrt(np.diag(corr_matrix_inv)))

# 標準化してpandas dataframe化
partial_corr_array = diag.dot(corr_matrix_inv).dot(diag)
partial_corr_matrix = pd.DataFrame(partial_cor_array, index=column_names, columns=column_names)
partial_corr_matrix = partial_corr_matrix.applymap(lambda x: -x if x < 1-1e-4 else 1)
partial_corr_matrix

penguin を使う方法

pandas やscikit-learnには偏相関係数を計算するメソッドがないため、統計解析パッケージPingouinを使って計算する。
（偏相関係数以外にも色々な統計解析ができるらしい）

インストール

conda install -c conda-forge pingouin

import pingouin as pg

#偏相関行列の計算
partial_corr_matrix = pg.pcorr(df)
partial_corr_matrix

参考

参考HP

1. グラフィカルモデリングで変数の相関関係を把握する(AlciaSolidProject|youtube)
2. データから因果を推定する
3. 偏相関係数(統計WEB)
4. Pythonで相関行列・偏相関行列（２）

関連書籍

Udemy

【ゼロからおさらい】統計学の基礎

いちばん理解できる統計学ベーシック講座その１
【確率分布・推定・検定】

いちばん理解できる統計学ベーシック講座その２
【相関分析・回帰分析】

豊富な演習でゼロから学ぶ統計学入門コース
検定3級レベルの知識をしっかり習得しよう