図:CausalImpactの推定結果例 (Kay H. Brodersen, et. al., 2015)[5]
実データを用いたPearlとRubinの因果推論
上記で概要を一通り述べましたので、実際のデータを利用して簡単なデモをしてみようと思います。
データ作成
まずは今回のデモで利用するライブラリをインポートします。見かけないライブラリpgmpy
は、Pearlの因果推論フレームワークにて、データからDAGを推定する際に利用します。
import numpy as np
import pandas as pd
from scipy.special import expit
import itertools
import warnings
warnings.simplefilter('ignore', FutureWarning)
from pgmpy.estimators import ConstraintBasedEstimator
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import auc, roc_curve
次に、データの生成です。以下のコードでは、x1, x2がzを生成し、x1, zがyを生成するようなデータを生成しています。また、データ生成のストーリーをDAGでも表現してみました。
class SampleDataset:
@classmethod
def generate_df(cls, num_data=10000, seed=0):
np.random.seed(seed=seed)
x1 = np.random.choice([0, 1], num_data, p=[0.3, 0.7])
x2 = np.random.choice([0, 1], num_data, p=[0.6, 0.4])
z = cls._get_logistic(num_data, [0.5, 1, -2], [1, x1, x2])
y = cls._get_logistic(num_data, [0.2, -1.5, 1], [1, x1, z], prob=False)
y0 = cls._get_logistic(num_data, [0.2, -1.5, 1], [1, x1, 0], prob=True)
y1 = cls._get_logistic(num_data, [0.2, -1.5, 1], [1, x1, 1], prob=True)
df = pd.DataFrame({
'x1': x1,
'x2': x2,
'z': z,
'y': y,
'y0': y0,
'y1': y1,
})
return df
def _get_logistic(num_data, b_list, x_list, prob=False):
e_z = np.random.randn(num_data)
res = 0
for b, x in zip(b_list, x_list):
res += b * x
res += e_z
z_prob = expit(res)
func = lambda x: np.random.choice([0, 1], 1, p=[1-x, x])
if prob:
return z_prob
else:
return np.vectorize(func)(z_prob)
df = SampleDataset().generate_df()
print(df.shape)
display(df.head())
図:生成したデータ例
図:デモで利用したデータ生成に対応するDAG
Pearlの因果推論
PearlはSCMに基づいた因果推論のフレームワークでした。ここでは、データのみを用いて、DAGを推定することに挑戦します。
class DAGEstimator:
def __init__(self, df):
self.est = ConstraintBasedEstimator(df)
self.cols = df.columns
self.num_variable = len(self.cols)
self.chk_pair = 2
self.test_0 = {c: [] for c in self.cols}
def test_0_dependency(self):
print('='*10 + '\n0次の独立性検定\n' + '='*10)
for pair in list(itertools.combinations(range(self.num_variable), self.chk_pair)):
c1 = self.cols[pair[0]]
c2 = self.cols[pair[1]]
cond = []
result = self.est.test_conditional_independence(
c1, c2, cond,
method='chi_square',
tol=0.05
)
if not result:
self.test_0[c1].append(c2)
self.test_0[c2].append(c1)
self._print_result(c1, c2, cond, result)
def test_1_dependency(self):
print('='*10 + '\n1次の独立性検定\n' + '='*10)
for key, value in self.test_0.items():
if len(value)<=1:
continue
c1 = key
for i, v in enumerate(value):
c2 = v
for j in range(len(value)):
if i==j:
continue
else:
cond = [value[j]]
result = self.est.test_conditional_independence(
c1, c2, cond,
method='chi_square',
tol=0.05
)
self._print_result(c1, c2, cond, result)
def test_n_dependency(self, c1, c2, cond):
print('='*10 + '\n任意の独立性検定\n' + '='*10)
result = self.est.test_conditional_independence(
c1, c2, cond,
method='chi_square',
tol=0.05
)
self._print_result(c1, c2, cond, result)
def _print_result(self, c1, c2, cond, result):
print("{}\tvs\t{}\t| {}\t: {}".format(
c1, c2, ', '.join(cond),
'独立' if result else '関連'
))
dag = DAGEstimator(df[['x1', 'x2', 'z', 'y']])
簡単に解析用のクラスを用意して、まずは全変数間の組み合わせに対して、独立性の検定を行ってみます。
結果としては、x1-x2の組み合わせのみが独立と判定できたため、グラフ構造のうち点線のエッジが切り離せることがわかります。
dag.test_0_dependency()
図:変数間の0次の独立性検定
dag.test_0_dependency()
図:検定結果から推定したグラフ
続いて、上記で残ったエッジに対して、もう一つ変数の条件付き確率を計算してみます。
すると、zが条件に付いたときx2-yが独立となるので、x2-y間のエッジは切り離せることがわかります。
dag.test_1_dependency()
図:変数間の1次の独立性検定
図:検定結果から推定したグラフ
最後に、エッジの向きを推定してみます。
グラフにV字構造を持つとき、条件付き確率を計算することでエッジの向きを計算することができます。今回はx2-z-x1にV字構造が存在しており、zを条件につけてx2-x1間の独立性検定を行ってみましょう。
すると、zを条件につけると、実際にはエッジが存在しないx2-x1間に関連があるという結果になりました。これはグラフ理論で登場する、合流点の問題となり、x2, x1からzへの有向エッジを書くことができるようになります。
dag.test_n_dependency('x1', 'x2', ['z'])
図:特定の変数間の1次の独立性検定
図:検定結果から推定したグラフ
データからDAGを推定できるのはここまでです。残りはデータ生成のストーリーに沿ってエッジの向きを考えていくことになります。
今回のデモの状況を考えると、「少なくとも効果yを知りたい状況なので、zからyにエッジの向きがあるだろう」「もしyからx1にエッジの向きがあると、巡回グラフとなってしまう」という2点を考慮して、最終的にはデータから正しいDAGを導くことになります。
この先の分析ストーリーとしては、先に登場した調整化公式を利用して介入効果(P(Y=y|do(Z=z)))を求めるわけなのですが、上述の通り、調整化公式自体がRubinの枠組みに登場するIPW推定と同じように定式化されるため、そちらに譲ることとします。
また、今回のDAG作成方法は、ベイジアンネットワークの手法を流用しています。ベイジアンネットワーク自体にも未観測のデータを推論する手法が存在します。
Rubinの因果推論
SCMには必ずしも基づくわけではないのがRubinの因果推論フレームワークでした。ここまでに生成したDAGは一旦忘れて、作成したデータセットからIPW推定を行ってみましょう。
class POEstimator:
def __init__(self, df):
self.df = df
self.df['z_0'] = 0
self.df['z_1'] = 1
def generate_prob_col(self, x: list, z: str, col: str, model=LogisticRegression()):
X = self.df[x]
y = self.df[z]
model.fit(X, y)
self.df[col] = model.predict_proba(X)[:, 1]
fpr, tpr, _ = roc_curve(y_true=y, y_score=self.df[col])
print('AUC ({})\t: {:.5f}'.format(col, auc(fpr, tpr)))
def get_true(self):
res = (self.df['y1'] - self.df['y0']).mean()
print('True causal effect\t: {:.5f}'.format(res))
def get_ipw(self):
res = self.df.apply(lambda x: self._calc_ipw(x), axis=1).mean()
print('IPW estimation\t: {:.5f}'.format(res))
def get_dr(self):
res = self.df.apply(lambda x: self._calc_dr(x), axis=1).mean()
print('DR estimation\t: {:.5f}'.format(res))
def _calc_ipw(self, x):
def _estimator(y, p_z, z):
return y / p_z * z
y1_hat = _estimator(x['y'], x['prob_z'], x['z'])
y0_hat = _estimator(x['y'], (1 - x['prob_z']), (1 - x['z']))
return y1_hat - y0_hat
def _calc_dr(self, x):
def _estimator(y, p_z, z, y_hat):
return (y / p_z) * z + (1 - z / p_z) * y_hat
y1_hat = _estimator(x['y'], x['prob_z'], x['z'], x['y1_hat'])
y0_hat = _estimator(x['y'], (1 - x['prob_z']), (1 - x['z']), x['y0_hat'])
return y1_hat - y0_hat
po = POEstimator(df)
先程作成したデータフレームを与えることで各種推定が実施できるようなクラスを準備して、いざ、IPW推定です。
po.generate_prob_col(['x1', 'x2'], 'z', 'prob_z')
po.get_true()
po.get_ipw()
図:IPW推定の結果
結果を見てみると、真の効果量とおおよそ同じ値が推定できていることがわかります。
さて、次にDR推定量を求めてみましょう。Rubinのフレームワークでは特段SCMを意識しないため、以下のようにyのモデリングにx2を追加してしまうかもしれません。
po.generate_prob_col(['x1', 'x2', 'z_0'], 'y', 'y0_hat')
po.generate_prob_col(['x1', 'x2', 'z_1'], 'y', 'y1_hat')
po.get_true()
po.get_dr()
図:DR推定の結果
モデリングには不要な変数x2を追加してしまっていますが、推定結果には特段影響がなさそうにみえます。これが中間変数のような変数だった場合、効果を過小見積もりしてしまうケースがありますが、今回はそのような結果は確認できませんでした。
Rubinの特徴として、SCMに基づくポテンシャルアウトカムに特段とらわれていないというように述べました。ここが、調整化公式だけでなく、種々の機械学習手法を導入できるポイントにもなります。
例えば、yの推定にランダムフォレストを用いることもできます。
po.generate_prob_col(['x1', 'z_0'], 'y', 'y0_hat', RandomForestClassifier())
po.generate_prob_col(['x1', 'z_1'], 'y', 'y1_hat', RandomForestClassifier())
po.get_true()
po.get_dr()
図:ランダムフォレストを用いたDR推定の結果
先に紹介したMeta-Learnerもこの実装の延長になりますので、ぜひ手を動かしてみて結果を確認いただければと思います。
最後に
統計的因果推論の体系を知りたいというモチベーションから調査を開始しましたが、調べれば調べるほど「先人たちよ、もう少し実務者にわかりやすく手法を一覧化しておくれよ」という気持ちが溢れてきました。学問体系として複数の流派が存在することは素晴らしいことですが、お互いに交通整理をして、統計的因果推論がよりビジネスに浸透していくことを期待しています。
次世代システム研究室では、ビッグデータ解析プラットホームの設計・開発を行うアーキテクトとデータサイエンティストを募集しています。興味を持って頂ける方がいらっしゃいましたら、ぜひ 募集職種一覧 からご応募をお願いします。
一緒に勉強しながら楽しく働きたい方のご応募をお待ちしております。
参考資料
[1] A Survey on Causal Inference (L. Yao, et. al., 2020)
[2] 入門 統計的因果推論(https://www.asakura.co.jp/detail.php?book_code=12241 )
[3] 岩波データサイエンス Vol. 3(https://www.iwanami.co.jp/book/b243764.html )
[4] 調査観察データの統計科学(https://www.iwanami.co.jp/book/b257892.html )
[5] Inferring causal impact using Bayesian structural time-series models (Kay H. Brodersen, et. al., 2015) (https://google.github.io/CausalImpact/CausalImpact.html )