相関係数の計算と可視化

データ解析・ケモインフォマティクスでは、ある化合物の物性など(目的変数)に対して、実験条件や記述子など(説明変数)が関係しているか調査するために、相関係数を計算することが必須の作業になる。
本記事では相関係数の種類や計算方法についてまとめた。

相関係数について

  • 2つの変数の間にある線形な関係の強弱を測る指標(-1 ~ 1)
  • 主にピアソンの(積率)相関係数, スピアマンの(順位)相関係数の2種類がよく使われる
  • ピアソンの相関係数は、間隔尺度、比尺度に対して利用だが、名義尺度、順序尺度に対しては利用できない(数字の値には意味がないため)
  • スピアマンの相関係数はピアソンの相関係数において相関係数を計算する前にデータを順位に変換した場合に当たる。
  • スピアマンの相関係数は名義尺度に対しても利用可能

利用上の注意

  • 相関係数はあくまでの(確率)変数の線形な関係の尺度であり、因果関係を説明するものではない
  • 相関係数は順序尺度であり比尺度ではないので、「相関係数0.2と0.4だと後者は2倍の相関がある」という主張は間違い。
  • 外れ値があると、その外れ値に引っ張られて相関がある様に見えてしまう(相関係数の絶対値が大きくなる)ので、注意が必要。(この問題を避けるためにも一度ペアプロット等を作成し、データの分布をチェックすることが重要)

参考)変数の尺度について

  • 名義尺度:カテゴリーなど大小に意味のないもの(電話番号など)
  • 順序尺度:順序や大小の意味はあるが、間隔には意味がないもの(競争の順位など)
  • 間隔尺度:値が等間隔になっており、その間隔に意味があるもの(温度など)
  • 比例尺度:0が原点であり、間隔と比率に意味があるもの(速度など)

pythonコード

データの準備

scikit-learn のBostonの家賃データを利用する。

import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
import seaborn as sns

# データ読み込み
boston = load_boston()
df_X = pd.DataFrame(boston.data, columns=boston.feature_names)

#可視化するには変数が多いので少し減らす
df_X = df_X.iloc[:,::4]
df_y = pd.DataFrame(boston.target,columns=['RENT'])
df = pd.concat([df_X,df_y],axis=1)

全変数の分布と相関をざっくり確認する

seabornのpairplotを使用する

sns.set_context('paper', font_scale=2.0)
pp = sns.pairplot(df)
pp.savefig('sns_pairplot.png')


2変数の相関をみる場合にはseabornのjoinplotを利用する

#pairplot で相関のありそうだった2変数でjointplotを作成してみる
jp = sns.jointplot(df.loc[:,'LSTAT'],df.loc[:,'RENT'],kind="reg")
jp.savefig('sns_joinplot.png')


相関係数の計算と可視化

ピアソン、スピアマンの相関係数ともに、pandas のcorrメソッドで一括計算計算可能

#全ての変数ペアの相関係数を求める。
df_pearson = df.corr(method='pearson')
df_spearman = df.corr(method='spearman')

#結果をヒートマップで可視化する
fig, ax = plt.subplots(figsize=(10, 8))
sns.heatmap(df_pearson, square=True, annot=True, fmt='.2f',annot_kws={'size':10},vmax=1, vmin=-1, center=0)
plt.title('Correlation coefficient (pearson)',fontsize=18)
plt.ylim(df_pearson.shape[1],0)
plt.savefig('Correlation_coefficient_pearson.png')
plt.show()

fig, ax = plt.subplots(figsize=(10, 8))
sns.heatmap(df_spearman, square=True, annot=True, fmt='.2f',annot_kws={'size':10},vmax=1, vmin=-1, center=0)
plt.title('Correlation coefficient (spearman)',fontsize=18)
plt.ylim(df_spearman.shape[1],0)
plt.savefig('Correlation_coefficient_spearman.png')
plt.show()

ピアソン・スピアマンの相関係数のヒートマップはそれぞれ以下の通り。




詳細(計算方法)

ピアソンの相関係数

2組の数値からなるデータ列 $ (x_1,y_1),(x_2,y_2), \cdots , (x_n, y_n) $ が与えられたとき、

標本共分散を $ s_{xy} $ 標本標準偏差を $s_x, s_y$ とおくと、相関係数 $r$ は下式で求められる。

$$ \begin{align} r & = \frac{s_{xy}}{s_xs_y}\\ & = \frac{\sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{\sum_{j=1}^{n}(x_j-\bar{x})^2} \sqrt{\sum_{k=1}^{n}(y_k-\bar{y})^2}} \end{align} $$

import numpy as np
import pandas as pd
from sklearn.datasets import load_boston

# データ読み込み
boston = load_boston()
df_X = pd.DataFrame(boston.data, columns=boston.feature_names)
df_y = pd.DataFrame(boston.target,columns=['RENT'])
df = pd.concat([df_X,df_y],axis=1)

#CRIM とRENT の相関係数を計算する
x = df.loc[:,'CRIM']
y = df.loc[:,'RENT']
sx = x.std()
sy = y.std()
sxy = np.cov(x,y)[0][1]
r = sxy / (sx * sy)
print('相関係数r:',round(r,2))
# >>> 相関係数r: -0.39

スピアマンの相関係数

  1. 生のスコアを順位に変換する
  2. 各観察(各ペア)における2つの変数の順位の差 D を計算する。
  3. ペアの数 N D を用いて以下の式から相関係数 $\rho$ を求める

$$ \rho = 1 - \frac{6\sum D^2}{N^3-N} $$

同順位がある場合

X、Yにおける同順位の個数をそれぞれ $n_x$、$n_y$、それらの順位を$t_i$、$t_j$ ($i=1, 2, \cdots, n_x; \quad j = 1, 2, \cdots, n_y$)として、以下の式を用いる。

ただし、同順位が少なければそれらを無視して最初の式を用いても影響は小さいとのこと。

$$ \rho = \frac{T_x + T_y - \sum D^2}{2\sqrt{T_x T_y}} $$

$$ T_x = \frac{N^3-N-\sum(t_i^3 - t_i)}{12} $$

$$T_y = \frac{N^3-N-\sum(t_j^3 - t_j)}{12} $$

import numpy as np
import pandas as pd
from sklearn.datasets import load_boston

boston = load_boston()
df_X = pd.DataFrame(boston.data, columns=boston.feature_names)
df_y = pd.DataFrame(boston.target,columns=['RENT'])
df = pd.concat([df_X,df_y],axis=1)

#CRIM とRENT の相関係数を計算する
x = df.loc[:,'CRIM'].values
y = df.loc[:,'RENT'].values

#それぞれの配列で順位づけを行う。(関数定義して計算する)
def calc_rank(num):
    num_sort = np.sort(num)
    num_rank = np.empty(len(num))

    i = 0
    while i < len(num_sort):
        val = num_sort[i]
        pos = np.where(num==val)[0]
        num_rank[pos] = i+1 / (len(pos))
        i += len(pos)

    return np.array(num_rank)

x_rank = calc_rank(x)
y_rank = calc_rank(y)
print(x_rank[:10])
print(y_rank[:10])
"""ourput
[  1.  25.  24.  34. 112.  29. 138. 199. 235. 213.]
[349.5        263.5        455.         450.5        464.5
 409.33333333 308.25       394.5        118.5        169.25]
"""

#各ペアにおける二つの変数の順位の差を計算する
D = x_rank - y_rank
D2 = D**2
N = len(D)
rho = 1 - (6 * D2.sum()) / (N**3-N)
print('スピアマンの相関係数ρ:',round(rho,2))
# >>> スピアマンの相関係数ρ: -0.56

参考

関連書籍