データから因果関係を導く！〜ルービン？パール？全部のせ〜

Pearlの因果推論

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']])

dag.test_0_dependency()

dag.test_0_dependency()

続いて、上記で残ったエッジに対して、もう一つ変数の条件付き確率を計算してみます。
すると、zが条件に付いたときx2-yが独立となるので、x2-y間のエッジは切り離せることがわかります。

dag.test_1_dependency()

最後に、エッジの向きを推定してみます。
グラフにV字構造を持つとき、条件付き確率を計算することでエッジの向きを計算することができます。今回はx2-z-x1にV字構造が存在しており、zを条件につけてx2-x1間の独立性検定を行ってみましょう。
すると、zを条件につけると、実際にはエッジが存在しないx2-x1間に関連があるという結果になりました。これはグラフ理論で登場する、合流点の問題となり、x2, x1からzへの有向エッジを書くことができるようになります。

dag.test_n_dependency('x1', 'x2', ['z'])

データから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()

 
結果を見てみると、真の効果量とおおよそ同じ値が推定できていることがわかります。

さて、次に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()

 
モデリングには不要な変数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()

先に紹介した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)