この記事はEnigmo Advent Calendar 2018の22日目です。
はじめに
https://atarimae.biz/archives/7374
交番と犯罪件数が正の相関があるからといって、交番を減らして犯罪件数は減らないですよね。
さて、データ分析を行う上では、相関関係と因果関係を切り分けることが重要になることがあります。
例えば、KPIとある数値xが相関しているとします。
x → KPI という因果関係であれば、xの操作でKPI向上の施策を検討することができます。
逆に、KPI → x という因果関係であれば、xを操作してもKPIは変化しません。
y → x なのか、x → y なのか、xとyの相関関係の有無だけでは、因果関係は分かりません。
この記事では、機械学習ブロフェッショナルシリーズ、統計的因果探索を参考にしています。
なお、、z → x , z → y という未観測の共通原因(交絡因子とも呼ばれる)zが存在する場合についても、書籍では扱われていますが、今回の記事では未観測の共通原因については扱わないことにさせていただきます。
LiNGAM
$x_{1}, x_{2}$の２変数の関係bを以下のように構造方程式とう呼ばれるモデルであると仮定します。
\begin{align}
x_{1} &= b_{12}x_{2} + e_{1} \\
x_{2} &= b_{21}x_{1} + e_{2} \\
\end{align}
$x_{1}, x_{2}$の２変数で現れ、$e_{1}, e_{2}$は外生変数がノイズ項です。
なお、LiNGAMでは非ガウス分布として扱います。
因果関係の推定の結果、$x_{1}→ x_{2}$という因果関係がある場合は、$b_{12}=0$となり下記のようになります。
\begin{align}
x_{1} &= e_{1} \\
x_{2} &= b_{21}x_{1} + e_{2} \\
\end{align}
ここで、LiNGAMは以下の方法があります
- (1)独立成分によるによるアプローチ
- (2)回帰分析と独立性によるアプローチ
独立成分によるによる推定
既にPythonで実装している方がいらっしゃるので、このアプローチに関してはそちらを参考にさせていただきます。
モジュールとともに、lingam.pyをimportします。
import pandas as pd
import numpy as np
from lingam import LiNGAM
scikit-learn付属のボストン不動産データセットを使います。
価格をyにするのが自然ですが、敢えて
価格:x
部屋数:y
とします。
boston = load_boston()
df = pd.DataFrame(boston.data, columns=boston.feature_names)
x = boston.target
y = df.RM.values
推定を実行します。
data = pd.DataFrame(np.asarray([x, y]).T, columns=['target', 'rooms'])
lg = LiNGAM()
lg.fit(data)
推定結果
rooms ---|9.102|---> target
部屋数が上がると、価格が上がる。という正しそうな因果関係を推定することができました。
回帰分析と独立性によるアプローチ
本をPythonで実装しました。ここでは、因果の向きのみの推定を行います。
なお、エントロピーは近似であることが本では述べられており、近似式の導出については、参考文献を読む必要があります。
def calc_r(x, y):
return ((x - np.mean(x * y) - np.mean(x)*np.mean(y)) / np.var(y) * y,
(y - np.mean(x * y) - np.mean(x)*np.mean(y)) / np.var(x) * x)
def normalize(x):
return (x - np.mean(x)) / np.std(x)
def entropy(u):
Hv = (1 + np.log(2 * np.pi)) / 2
k1 = 79.047
k2 = 7.4129
gamma = 0.37457
return Hv - k1*(np.mean(np.log(np.cosh(u))) - gamma)**2 - k2 * (np.mean(u * np.exp(-1 * u**2 /2)))**2
def lingam_reg(x, y, columns):
xr, yr = calc_r(x, y)
m = entropy(normalize(x)) + entropy(normalize(xr) / np.std(xr)) \
- entropy(normalize(y)) - entropy(normalize(yr) / np.std(yr))
if m >= 0:
print('{0} ---> {1}'.format(columns[0], columns[1]))
else:
print('{0} ---> {1}'.format(columns[1], columns[0]))
推定を実行します。
data = pd.DataFrame(np.asarray([x, y]).T, columns=['target', 'rooms'])
lg = LiNGAM()
lg.fit(data)
推定結果
rooms ---> target
こちらも正しそうな推定結果を偉ました。
部屋数以外の変数でもやってみる
- INDUS: 小売業以外の商業が占める面積の割合
y = df.INDUS
X = pd.DataFrame(np.asarray([x, y]).T, columns=['target', 'INDUS'])
lg = LiNGAM()
lg.fit(X)
lingam_reg(x, y, columns=['target', 'INDUS'])
推定結果
INDUS ---|-0.648|---> target
INDUS ---> target
小売業以外の商業、つまりオフィスの面積の割合が増えると、不動産価値は上がる。正しい気がする。
- 税率
y = df.TAX
X = pd.DataFrame(np.asarray([x, y]).T, columns=['target', 'TAX'])
lg = LiNGAM()
lg.fit(X)
lingam_reg(x, y, columns=['target', 'TAX'])
推定結果
target ---|-8.586|---> TAX
target ---> TAX
不動産価値が高いと、税率が上がる。今までは逆の因果関係だが、この場合は正しそう。
最後に
ごく簡単な例ですが、統計的因果探索であるLiNGAMを体験することができました。
前述の構造方程式の通り、変数間の関係が線形である等、仮定があることは注意しなければいけませんが、あくまで因果関係の向きを推定する上では、その仮定の上でも議論が進められると思います。
また、こちらの取り組みについても興味深いです。
https://qiita.com/MorinibuTakeshi/items/402cb905e70655724d35
未観測の共通要因については今回割愛してしまいましたが、LiNGAMは因果推論という枠組みの一つの取り組みであり、傾向スコアやグレンジャー因果性といった取り組みについても今後学習を進めていきたいと思います。