動画概要

因果推論の第4弾までは、EconMLを使うための外堀を埋めるための基礎的事項を取り上げてきました。第5弾となる本動画では、EconMLの実際の活用を念頭に、その推論時のメソッドであるeffectメソッドとmarginal_effectメソッドについて解説します。また動画後半では、2020年4月に公開されたPyCaretを機械学習ベースの因果推論業務のどの時点でどのような温度感で使えそうかについて考察しました。

参照

・EconML Problem Setup and API Design
https://econml.azurewebsites.net/spec/api.html

・WhyNot Opioid RCT data generator
https://github.com/zykls/whynot/blob/master/examples/causal_inference/heterogeneous_treatment_effects.ipynb

・PyCaret
https://pycaret.org/

利用コード

利用ライブラリのインポート。whynotはデータ生成用途としてのみここでは使います。

import pandas as pd
import numpy as np

import econml
import whynot as wn

from IPython.display import Image, display
import warnings
warnings.filterwarnings('ignore')

from pycaret.regression import *

numpy配列をpandasのデータフレームにしてリターンするだけの関数です。

def get_ncols(dataset):
    try:
        return dataset.shape[1]
    except:
        return 1

def get_df(dataset, Cap):
    return pd.DataFrame(dataset, columns=[Cap+str(i) for i in range(get_ncols(dataset))])

オピオイドのRCTデータを生成しています。各行処置変数(T)、交絡因子(X)、アウトカム変数(Y)、シミューレション上の実効果(true effect：処置なしから処置ありに変化させたときの効果）が生成されます。

n_sample = 5000
experiment = wn.opioid.RCT
data = experiment.run(num_samples=n_sample)

Y = data.outcomes     # continuous
T = data.treatments   # binary
X = data.covariates   # continuous
E = data.true_effects # individual treatment effect (true effect)

Y_, T_, X_ = get_df(Y,'Y'), get_df(T,'T'), get_df(X,'X')
Dataset = Y_.join(T_).join(X_)

print(Y_.shape, T_.shape, X_.shape)
display(Dataset.head())

EconMLのDML(Double Machine Learning)ベースのCATE推定器をセットします。簡単化のためDMLに限定し、処置変数とアウトカム変数の推定もツリー系とリニア系の２種類にしています。

from econml.dml import BaseCateEstimator, DMLCateEstimator
from sklearn.linear_model import Lasso, Ridge, ElasticNet, LogisticRegression
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
from sklearn.ensemble import GradientBoostingRegressor, GradientBoostingClassifier

estimators = {
    'dml_rf':
        DMLCateEstimator(
            model_y=RandomForestRegressor(),
            model_t=RandomForestClassifier(),
            model_final=RandomForestRegressor(),
            discrete_treatment=True), # default: False

    'dml_ridge':
        DMLCateEstimator(
            model_y=Ridge(),
            model_t=LogisticRegression(),
            model_final=Ridge(),
            discrete_treatment=True),

    'dml_lasso':
        DMLCateEstimator(
            model_y=Lasso(),
            model_t=LogisticRegression(),
            model_final=Lasso(),
            discrete_treatment=True),

    'dml_elnet':
        DMLCateEstimator(
            model_y=ElasticNet(),
            model_t=LogisticRegression(),
            model_final=ElasticNet(),
            discrete_treatment=True),
}

DMLによるATE(Average Treatment Effect)推定値の評価のため、一部データをホールドアウトし、評価用データにおけるATEの推定値を取得します。この作業を30回繰り返しました。

# Repeated holdout for dml-estimators' evaluation
ate, train_, test_ = {}, [], []
for i in range(30):
    # holdout
    rnd_idx = np.array(range(n_sample))
    np.random.shuffle(rnd_idx)
    train_idx, test_idx = rnd_idx[:int(n_sample*0.7)], rnd_idx[int(n_sample*0.7):]
    Y_train, T_train, X_train = Y[train_idx], T[train_idx], X[train_idx]
    Y_test, T_test, X_test = Y[test_idx], T[test_idx], X[test_idx]

    # prepare T0/T1 to estimator ate from zero to one
    rows_train = T_train.shape[0]
    rows_test = T_test.shape[0]
    T0_train, T1_train = np.zeros((rows_train, get_ncols(T_train))), np.ones((rows_train, get_ncols(T_train)))
    T0_test, T1_test  = np.zeros((rows_test, get_ncols(T_test))), np.ones((rows_test, get_ncols(T_test)))

    # evaluation by ate(averate treatment effect) comparison
    for est_name, est in estimators.items():
        # print(est_name)
        est.fit(Y_train, T_train, X_train)
        
        if i==0:
            ate[(est_name,'train')] = []
            ate[(est_name,'test')] = []

        ate_train = np.mean(est.effect(X_train, T0=T0_train, T1=T1_train))
        ate_test = np.mean(est.effect(X_test, T0=T0_test, T1=T1_test))
        ate[(est_name,'train')].append(ate_train)
        ate[(est_name,'test')].append(ate_test)

最後のターンだけですが、訓練用と評価用データセットの実際のATEを出力します。この部分は拡張してもう少し丁寧に見ることができますが、傾向としてどちらの数値もほぼ同一です。

print('true ate(train): {:.3f}'.format(np.mean(E[train_idx])))
print('true ate(test) : {:.3f}'.format((np.mean(E[test_idx]))))

繰り返しホールドアウトの推定ATEの平均値と標準偏差を取得します。

mean_ = pd.Series(ate).unstack().applymap(np.mean).reindex(['train','test'], axis=1).sort_values('train')
std_ = pd.Series(ate).unstack().applymap(np.std).reindex(['train','test'], axis=1).sort_values('train')
display(mean_.join(std_, lsuffix='_avg', rsuffix='_std'))

本記事は以上です。