🙊

偏相関係数の概要と求め方(python)

2022/07/26に公開

偏相関係数とは

ある変数(x_1,x_2,\dots, x_n)がある時、他の変数全てを固定した時の x_ix_j の相関係数を x_ix_j偏相関係数という。

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

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

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

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

偏相関係数の計算方法

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

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

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

この時、z 因子の影響を除いた xy の偏相関係数 \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_1x_2 の偏相関係数を求めるには、x_1x_2以外の成分を使って、x_1x_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-2t分布に従う。 棄却される時偏相関係数は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()

.dataframe tbody tr th:only-of-type { vertical-align: middle; } <div></div> .dataframe tbody tr th { vertical-align: top; } <div></div> .dataframe thead th { text-align: right; }

動画時間_s 視聴回数 コメント 高評価件数 低評価件数
0 1082 714 18 69 0
1 1025 1533 14 98 0
2 1489 1219 23 111 0
3 530 2118 28 164 2
4 1283 1462 18 112 1

愚直に計算(モジュール不使用)

# 相関係数を計算
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

.dataframe tbody tr th:only-of-type { vertical-align: middle; } <div></div> .dataframe tbody tr th { vertical-align: top; } <div></div> .dataframe thead th { text-align: right; }

動画時間_s 視聴回数 コメント 高評価件数 低評価件数
動画時間_s 1.0000 -0.0577 0.2602 0.1509 -0.1315
視聴回数 -0.0577 1.0000 0.4250 0.7019 0.7329
コメント 0.2602 0.4250 1.0000 0.1184 -0.1821
高評価件数 0.1509 0.7019 0.1184 1.0000 -0.3367
低評価件数 -0.1315 0.7329 -0.1821 -0.3367 1.0000

penguin を使う方法

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

インストール

conda install -c conda-forge pingouin
import pingouin as pg

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

.dataframe tbody tr th:only-of-type { vertical-align: middle; } <div></div> .dataframe tbody tr th { vertical-align: top; } <div></div> .dataframe thead th { text-align: right; }

動画時間_s 視聴回数 コメント 高評価件数 低評価件数
動画時間_s 1.0000 -0.0577 0.2602 0.1509 -0.1315
視聴回数 -0.0577 1.0000 0.4250 0.7019 0.7329
コメント 0.2602 0.4250 1.0000 0.1184 -0.1821
高評価件数 0.1509 0.7019 0.1184 1.0000 -0.3367
低評価件数 -0.1315 0.7329 -0.1821 -0.3367 1.0000

参考

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

関連書籍

以上

Discussion