Pythonによる因果グラフ推定 -causalnexの紹介その2- - データサイエンティスト（仮）の暫定解更新ブログ

f:id:Yuminaga:20210523151713j:plain — Image by Pixabay

今回の記事は前回記事に引き続き因果探索についてです。

前回紹介できなかったcausalnexの$\ell_1-$制約付きのNOTEARS推定、pytorchによるfrom_pandas実装、sklearnライクにNOTEARSを使うことができるDAGRegressor, DAGClassifierについて紹介していきます。

Causalnexを用いてPythonでサンプルデータ実験
まとめ

Causalnexを用いてPythonでサンプルデータ実験

準備

サンプルデータや、前回と同様に下図の因果グラフを想定したデータを使います。またコードはこちらに公開しています。

f:id:Yuminaga:20210523021252p:plain


from util.networkx_util import NetworkxUtil
nx_util = NetworkxUtil() # networkxの図示, printのためのutil

# サンプル生成
n_samples = 2000
struct_data = generate_sample_data(n_samples)

# sklearnっぽくX, yに分割
X = struct_data.drop(["y"], axis=1)
y = struct_data["y"]

今回はsklearnっぽくデータを特徴量$X$, 予測対象の$y$という分割したデータも作っておきます。

Lassoを用いたNOTEARS

まず$\ell_1-$penaltyを加えてNOTEARSを推定してみます。通常のNOTEARSはfrom_pandasで推定できました。$\ell_1-$penalty付きはfrom_pandas_lassoで、from_pandasと同じように推定することができます。($\ell_1-$penaltyの強さのパラメータbetaが加わります。)

実際に推定するコードは以下です。


from causalnex.structure.notears import from_pandas_lasso
# NOTEARSを実行, from_pandas_lassoでL1 penaltyをつけて推定することが可能
sm = from_pandas_lasso( 
    struct_data,
    beta = 0.15, # L1 penalty の強さ
    tabu_edges = [],
    tabu_parent_nodes = None,
    tabu_child_nodes = None,
)
sm.threshold_till_dag() # DAGになるように閾値をあげる

# weights print
nx_util.print_weights(sm)

# plot sturecture model
nx_util.plot_structure_model(sm, layout_method="shell", layout_seed=1, figsize=(5,5))

結果を見ると、今回はremove_edges_below_thresholdを用いずとも、$\ell_1-$penaltyによって、正しく因果グラフの矢印の向きが推定されていることがわかります。正則化がうまく働いていますね。

f:id:Yuminaga:20210523170443p:plain


# 出力 (lasso)
[ x1 ] ---> [ y ]		Weight	1.08730
[ c1 ] ---> [ x1 ]		Weight	0.29201
[ c1 ] ---> [ y ]		Weight	0.24506
[ c2 ] ---> [ x1 ]		Weight	0.73727
[ c2 ] ---> [ y ]		Weight	0.14375

Pytorchを用いたNOTEARS

次にpytorchを用いたNOTEARS推定を試してみます。import先が異なるだけでfrom_pandasというmethod名は同じです。今回は区別するためにfrom_pandas_pytorchという名前で呼び出します。

pytorch版のfrom_pandasはデータの種類(binary or continuous)や$\ell_1, \ell_2-$penaltyなど、より細かい項目を設定することができるようです。実際につかってみます。


from causalnex.structure.pytorch import from_pandas as from_pandas_pytorch

dist_type_schema = {
    "x1" : "cont",
    "x2" : "cont",
    "c1" : "cont",
    "c2" : "bin",
    "y" : "cont"
}

# data typeを加味してpytorchを用いて推定
sm_pytorch = from_pandas_pytorch(
    struct_data, 
    dist_type_schema = dist_type_schema, 
    hidden_layer_units = None,
    lasso_beta = 0,  # L1
    ridge_beta = 0 , # L2
    w_threshold = 0.,  # edge_weight閾値
    use_bias = False
)
sm_pytorch.threshold_till_dag()
sm_pytorch.remove_edges_below_threshold(0.1) # 

# weights print
nx_util.print_weights(sm_pytorch)

# plot sturecture model
nx_util.plot_structure_model(sm_pytorch, layout_method="shell", layout_seed=1, figsize=(5,5))

結果を見てみると、こちらもいい感じに正しく推定できました。

f:id:Yuminaga:20210523163456p:plain


# 出力 (pytorch)
[ x1 ] ---> [ y ]		Weight	1.02769
[ c1 ] ---> [ x1 ]		Weight	0.82515
[ c1 ] ---> [ y ]		Weight	0.96580
[ c2 ] ---> [ x1 ]		Weight	0.96923
[ c2 ] ---> [ y ]		Weight	0.55306

今回は設定しませんでしたが、from_pandas_pytorchでhidden_layer_unitsを[10]のように指定した場合には、活性化関数がsigmoid関数で、10ノードもつ中間層となるニューラルネット（MLP）が設定され、推定に使われるようです。*1

Scikit-learnライクなNOTEARS推定インターフェース

sklearnライクにNOTEARSを使えるインターフェースとして、DAGRegressorとDAGClassifierが用意されています。どちらともインスタンスを作ってfit, predictで使うことができる上、pytorch版のfrom_pandasのように細かい設定まで設定することができるので使い勝手がよさそうです。

DAGRegressor

以下にDAGRegressorを使ったコードを示します。設定はfrom_pandas系とは名前が変わっているものの、基本的には同じです。 alpha, betaで$\ell_1, \ell_2-$penaltyを設定できる上、enforce_dag = Trueにより推定結果がDAGになるようにすることができます。またdependent_target = Trueを用いることで$y$が子ノードにならないような制約を加えることもできます。


from causalnex.structure import DAGClassifier, DAGRegressor

# sklearn
reg = DAGRegressor(
    threshold= 0.1,  # edge weightの閾値
    #dist_type_schema = dist_type_schema,  dist_typeを指定することも可能
    alpha = 0., # L1
    beta = 0., # L2 
    hidden_layer_units = None,
    standardize=False,
    enforce_dag=True,
    tabu_edges=None,
    tabu_parent_nodes=["y"],
    tabu_child_nodes=None,
    dependent_target = True # yは親ノードにならない制約
)
# sklearnライクでfit, predict
reg.fit(X, y)
y_pred = reg.predict(X)
print(y_pred)

という形で予測まで行うことができます。


# y_pred 出力
[ 6.62431383  0.81902099  1.38599503 ... -0.01556239  1.82589388
  2.93284011]

因果グラフに関してはreg.graph_を参照すれば見れるため、同様にplotしてみます。


# plot structure model
sm = reg.graph_
sm = sm.get_largest_subgraph()
nx_util.plot_structure_model(sm, layout_method="random", figsize=(6,6), layout_seed=1)

出力するととなり、真の因果グラフの矢印を推定できている上で、予測値を算出できていることがわかります。

f:id:Yuminaga:20210523165255p:plain

DAGClassifier

DAGClassifierも同様に使うことができます。2値分類タスクにするためにデータを変換しておきます。


# 2値分類タスクに変換
y = y > np.median(y)
y = y.astype(int)

clf = DAGClassifier(
    threshold= 0.1, # edge weightの閾値
    #dist_type_schema = dist_type_schema,
    alpha = 0.0, # L1
    beta = 0.0,
    hidden_layer_units =None,
    standardize=False,
    enforce_dag=True,
    tabu_edges=None,
    tabu_parent_nodes=["y"],
    tabu_child_nodes=None,
    dependent_target = True # yは親ノードにならない制約
)

clf.fit(X, y)
y_pred = clf.predict(X)
y_pred_proba = clf.predict_proba(X)

print("Pred Label" ,y_pred)
print("Pred Proba ",y_pred_proba)

# plot structure model
sm = clf.graph_
sm = sm.get_largest_subgraph()
nx_util.plot_structure_model(sm, layout_method="random", figsize=(6,6))

推定結果は以下になります。


# 出力    
Pred Label:  [1 0 0 1 1]
Pred Proba : [[1.72674656e-04 9.99827325e-01]
 [8.47911775e-01 1.52088225e-01]
 [6.12419337e-01 3.87580663e-01]
 [4.54858541e-01 5.45141459e-01]
 [2.28357315e-03 9.97716427e-01]]

通常のsklearn内のclassifierと同様の使い勝手ですね。predictでの結果はラベルの予測値で、predict_probaは各クラスの所属確率になっています。

f:id:Yuminaga:20210523165312p:plain