はじめに

機械学習モデル等を用いた因果推論手法であるDouble Machine Learningについて、表面的には理解していましたが、詳しい部分まで知らなかったためいくつかのドキュメントを読みました。

簡単にまとめたいと思います。

ベースはMSのEconMLのページを参考にしています。

econml.azurewebsites.net

Double Machine Learning

特徴

二段階のモデル推定を行う
推定に用いる手法やモデルは柔軟に設定することができる
ネイマン直交条件を満たすスコア関数の零点を求めることで処置結果やその信頼区間を推定できる
モデルが誤推定されていなければ推定結果は√n-consistent になる

モデル

Double Machine Learningでは、部分線形モデルを考えることになる。部分線形モデルは線形モデルより表現力が高い。

$\displaystyle{ \begin{split} Y =~& \theta \cdot T + g(X, W) + \epsilon \\ \end{split} }$

ここでそれぞれの変数は、アウトカム $Y$ ,
※ 線形モデル

$\displaystyle{ Y = \theta T + b^T X + c^T W + \epsilon }$

Double Machine Learningでは、データ生成プロセスに関して次の構造方程式の仮定しています。

$\displaystyle{ \begin{split} Y =~& \theta(X) \cdot T + g(X, W) + \epsilon \\ T =~& f(X, W) + \eta \\ \end{split} }$

$\displaystyle{ \begin{split} E [ \epsilon | X, W ] = 0 \\ E [ \eta | X, W ] = 0 \\ E [ \eta \cdot \epsilon | X, W ] = 0 \end{split} }$

このモデルを持ちいてmarginal CATE : $\theta (X)$ を推定することを目指している。

１段階目の推定では条件付き期待関数 $g, f$ を推定する。このモデルには任意のノンパラメトリックモデル・機械学習モデルを用いることができる。

$\displaystyle{ \begin{split} g(X, W) = E [ Y | X, W ] \\ f(X, W) = E [ T | X, W ] \end{split} }$

2段階目では、 $\theta (X)$ を推定する。構造方程式を $T$ と $Y$ の残差同士の回帰の形に書き直す。これにより式から $g(X, W)$ を削除することができる。

$\displaystyle{ Y - E [ Y | X, W ] = \theta (X) \cdot (T - E [ T | X, W ] ) + \epsilon }$

それぞれの残差は次のように計算できる。

$\displaystyle{ \tilde {Y} = Y - q(X, W) }$

$\displaystyle{ \tilde {T} = T - f(X, W) = \eta }$

残差同士の回帰は以下のようになる。

$\displaystyle{ \tilde {Y} = \theta (X) \cdot \tilde {T} + \epsilon }$

最終的には次の目的関数を満たす問題に帰着することができる。

$\displaystyle{ \hat { \theta } = \arg \min_{ \theta \in \Theta} E _n \left [ ( \tilde {Y} - \theta (X) \cdot \tilde {T})^2 \right ] }$

EconMLでは、上記の目的関数に正則化項を加えて推定を行うこともできる。

部分線形モデルについて

部分線形モデルにおけるnaibeな推定量では、√n-consistentではなくなるため、推定精度が良くない。

要因としては $g,f$ の推定精度が関係している。

機械学習モデルの一般的な収束レートは、n^1/4-consistent()である。この時の $\theta$ の収束レートは()となってしまい、うまく推定することができない。

このバイアスをRegularization Biasと呼ばれている。

対策として操作変数法の考え方を用いたモーメント条件を考えると、推定量の収束レートは $g,f$ の推定精度の掛け算項が支配的になっている。

MLモデルの一般的な収束レートはn^1/4-consistent()であるため、その掛け算であるため収束レートは√n-consistent()となる。

これを達成する方法が上記に記したYとTの残差同士の回帰であり、Robinsonの手法と呼ばれている。

EconMLでの実行

CATEを推定するためのPythonライブラリ：EconMLをMicrosoftが提供しているです。

様々な因果推論+機械学習のモデルが簡単に実行することができる。

先に用いるサンプルデータを生成する。基本的にはEconMLのサンプルコードをベースにさせてもらっている。

github.com

基本的な使い方

提供されているモデルのタイプは大まかに3タイプある模様。いくつかのパターンの挙動を確認する。

LinearDML 線形の因果効果を仮定する任意のscikit-learn線形推定器等を最終段階として定義することができる方法
CausalForestDML Causal Forest を最終モデルとして使用する方法
NonParamDML 最終的な因果効果の推定にモデルを仮定しない方法

import numpy as np
from econml.dml import DML, LinearDML, CausalForestDML
from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor

それぞれのモデルの場合の処理を行う。

はじめにLinearDML。

結果変数を推定するモデルと処置を推定するモデルを指定する必要がある。今回はRandomForestRegressorを利用。因果効果は線形 $\theta$ であると考えることになる。

est1 = LinearDML(model_y = RandomForestRegressor(),
                　　　　model_t = RandomForestRegressor())
est1.fit(Y, T, X = X, W = W)  # モデルの学習
te_pred1 = est1.effect(X_test) # CATEの推定
te_lower1, te_upper1 = est1.effect_interval(X_test, alpha=0.01)

非線形の因果効果 $\theta (X)$ を推定する場合は、以下のように2段階目の推定に多項式モデルを仮定する。

# CausalForestDML
est3 = CausalForestDML(model_y=RandomForestRegressor(), 
                       model_t=RandomForestClassifier(min_samples_leaf=10), 
                       discrete_treatment=True, 
                       n_estimators=200, 
                       min_impurity_decrease=0.001, 
                       verbose=0, 
                       cv=6)
est3.tune(Y, T, X=X, W=W)
est3.fit(Y, T, X=X, W=W)
te_pred3 = est3.effect(X_test)
te_lower3, te_upper3 = est3.effect_interval(X_test, alpha=0.01) # 1-α %の信頼区間を出す

またCausalForestDMLでも非線形の因果効果 $\theta (X)$ を推定することができる。

# CausalForestDML
est3 = CausalForestDML(model_y=RandomForestRegressor(), 
                       model_t=RandomForestClassifier(min_samples_leaf=10), 
                       discrete_treatment=True, 
                       n_estimators=200, 
                       min_impurity_decrease=0.001, 
                       verbose=0, 
                       cv=6)
est3.tune(Y, T, X=X, W=W)
est3.fit(Y, T, X=X, W=W)
te_pred3 = est3.effect(X_test)
te_lower3, te_upper3 = est3.effect_interval(X_test, alpha=0.01)

NonParamDMLでは最終的な因果効果の推定にモデルを仮定しておらず、自由に推定モデルを設定することができる。

model_final引数にsklearnのモデルを設定する。

est4 = NonParamDML(model_y = GradientBoostingRegressor(),
                   model_t = GradientBoostingRegressor(),
                   model_final = GradientBoostingRegressor())
est4.fit(Y, T, X=X, W=W)
te_pred4 = est4.effect(X_test)

CausalForestDMLと同様にCATEを推定することができているが、分散がかなり大きいように見える。

推定モデルのサマリを確認する

est1.summary()

est3.summary()

2段階目のモテリングにおける交差検証で算出された誤差も確認することができる。

print("Linear Model score:",round(est1.score_,4))
print("Poly Model score:",round(est2.score_,4))
print("DML Model score:",round(est3.score_,4))
print("Non Model score:",round(est4.score_,4))

Linear Model score: 1.3845
Poly Model score: 1.0754
DML Model score: 0.936
Non Model score: 0.9862

そのほかの設定

SHAPによる寄与の算出

因果効果を推定する2段階目のモデルに対してSHAPを適用することで、因果効果の推定における各特徴量の寄与を推定することができる。

import shap
shap_values = est3.shap_values(X)
shap.plots.beeswarm(shap_values['Y0']['T0_1'])

単一な決定木による説明の獲得

CATEの説明をよりわかりやすくするために、獲得した因果効果を単一な決定木により表現する。これにより複雑な因果効果の場合分けを比較的人間が理解できる結果に変換してくれる。

from econml.cate_interpreter import SingleTreeCateInterpreter
plt.figure(figsize=(5, 8))

intrp = SingleTreeCateInterpreter(include_model_uncertainty=True, max_depth=2, min_samples_leaf=10)
intrp.interpret(est3, X_test)
plt.figure(figsize=(25, 5))
intrp.plot(feature_names=["X1","X2","X3","X4"], fontsize=12)

パイプラインの利用

モデルの設定で特徴量変換を設定することができる。定義した変換処理をfeaturizerに渡すだけで良い。

from sklearn.ensemble import RandomForestRegressor
class LogFeatures(object):
    def fit(self, X, y=None):
        return self
    def transform(self, X):
        return np.concatenate((X, np.log(1+X)), axis=1)
    def fit_transform(self, X, y=None):
        return self.fit(X).transform(X)

est = LinearDML(model_y=RandomForestRegressor(),
                model_t=RandomForestRegressor(),
                featurizer=LogFeatures())
est.fit(y, T, X=X, W=W)

featurizerにはパイプラインも設定することもできる。

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
est = LinearDML(model_y=RandomForestRegressor(),
                model_t=RandomForestRegressor(),
                featurizer=Pipeline([('log', LogFeatures()),
                                     ('poly', PolynomialFeatures(degree=2))]))

固定効果の導入

from sklearn.preprocessing import OneHotEncoder
X_oh = OneHotEncoder(sparse_output=False).fit_transform(X)[:, 1:]

est = LinearDML(model_y=RandomForestRegressor(),
                             model_t=RandomForestRegressor())
est.fit(y, T, X=X_oh, W=W)

複数の処置がある場合

est = LinearDML()
est.fit(y, np.concatenate((T0, T1), axis=1), X=X, W=W)

なお、複数の処理と複数の結果変数がある場合も対応することができる。

適用例

こちらのノートブックでより具体的なサンプルデータにおける分析の流れを説明している。 (サンプルデータの生成についてはここでは割愛) github.com

ここでは所得に応じた需要の価格弾力性を推定することを目指す。

用いるのは以下のようなモデル。

ここで、 $log(Y)$ は需要の自然対数、 $\theta (X)$ 所得 $X$ に基づいた需要の価格弾力性、 $log(T)$ は価格の自然対数、 $f(X,W)$ は、所得 $X$ とユーザー特性 $W$ に基づいた需要に対する他の要因の関数、 $g(X,W)$ は $X, W$ に基づいた価格決定に対する要因の関数。

多項式回帰を仮定したパラメトリック推定を行う。実行するには以下のとおり。

est = LinearDML(model_y = GradientBoostingRegressor(),
                model_t = GradientBoostingRegressor(),
                featurizer=PolynomialFeatures(degree=5, include_bias=False))
est.fit(Y, T, X = X, W = W) 
te_pred = est.effect(X_test) 
te_pred_interval = est.effect_interval(X_test, alpha=0.01) # 1-α %の信頼区間を出す

推定されたモデルは以下のようになっている。

Coefficient Results
point_estimate  stderr  zstat   pvalue  ci_lower    ci_upper
income  -50.646    7.043  -7.191 0.0    -64.45 -36.842
income^2   83.662 9.41   8.891  0.0    65.219 102.106
income^3   -43.952    5.3    -8.293 0.0    -54.34 -33.565
income^4   9.301  1.28   7.265  0.0    6.792  11.81
income^5   -0.683 0.109  -6.289 0.0    -0.896 -0.47
CATE Intercept Results
point_estimate  stderr  zstat   pvalue  ci_lower    ci_upper
cate_intercept  -8.698 1.751  -4.967 0.0    -12.131    -5.266

次にCausal Forestを用いたノンパラメトリックな推定を行う。実行するには以下のとおり。

est = CausalForestDML(
    model_y=GradientBoostingRegressor(), model_t=GradientBoostingRegressor()
)
est.fit(log_Y, log_T, X=X, W=W, inference="blb")

te_pred = est.effect(X_test)
te_pred_interval = est.effect_interval(X_test)

次に収益を最大化する施策割り当て方策(ポリシー)を考えたい。

EconMLライブラリでは、シンプルな説明可能なポリシーを生成するためにSingleTreePolicyInterpreterが提供されている。以下のように実行すると介入費用と介入効果の関係から、ターゲットにするべき対象を判別するルールを生成できる。

intrp = SingleTreePolicyInterpreter(risk_level=0.05, max_depth=2, min_samples_leaf=1, min_impurity_decrease=0.001)
intrp.interpret(est, X_test, sample_treatment_costs=-1)
plt.figure(figsize=(25, 5))
intrp.plot(feature_names=X.columns, treatment_names=["Discount", "No-Discount"], fontsize=12)