結局ROC AUCはなにものかPython実装して定義への理解を深める

今回はROC AUC (Receiver Operating Characteristic curve Area Under the Curve)についてです。

2値分類の機械学習モデルの予測スコア精度指標としてめちゃめちゃ使われる印象があるROC AUCですが、結構説明は難しいですよね。

説明しようとすると、ROC AUCとはこのカーブの下側の面積です、このカーブがROCカーブといって、TPr、FPrというものをプロットして、、閾値を変えて、、、閾値というのは、、、といったように一筋縄では説明できないな、と思ったことがあります。

ただ、実はROC AUCは、「正例(Pos)のサンプルの予測確率が、負例(Neg)のサンプルの予測確率に比べて高くなっているか」、つまり正しい順序で各サンプルの予測確率が並んでいるかの割合、というとても理解しやすい定義で表すことができます。加えてその定義をみると、ROC  AUCはクラス比率の逆数によって、それぞれのクラスの精度を重みづけている指標ともとれることがわかります。 

また、その定義にのっとって自分でROC AUCをPython実装してsklearn実装と比較してみました。

ROC AUCの図の見方

混合行列の整理

ROC AUCについて見る前に、TP, FP, FN, TNについて整理します。下図に表しました。

2値分類において、予測クラスと実際のクラスが合っている場合は、Trueで表され、間違っている場合はFalse。そして予測クラスがPositiveクラス（正例, 1）である場合にPositive、Negativeクラス（負例, 1）である場合Negativeとすると、上のようにTP, FP, FN, TNが定義されます。 

ROC Curve

それでは、ROC AUCを見ていきます。ROC AUCの概要を下図で表現しました。$x$軸はFPrで、$y$軸はTPrです。このFPr, TPrの点を結んだ曲線がROC Curveであり、その下側の面積がROC AUCとなります。またROC AUCの重要な性質として必ず0から1の範囲になるという性質があります。

上図の指標の整理を見ると、$x$軸であるFPrは、実際のクラスがNegのうち間違えてPosと予測されてしまった割合。$y$軸であるTPrは、実際のクラスがPosのうち正しくPosと予測された割合 ( = Recall )で表されます。

 

FPr, TPrの定義を確認すると、予測値がすべてNeg(0)であった場合には分子が0になるため、$(0,0)$となることがわかります。予測値がすべてPos(1)である場合には、TN, FNが0となり、分母分子が一致するので$(1,1)$を通ることになります。また、完全にランダム予測の場合は、ROC Curveが$y = x$線になりROC AUCは0.5となります。

ROC AUCの定義

それではROC AUCをよりシンプルな形で表現をしたいと思います。

準備として、実際のクラスがPos(1)であるサンプルの集合を$S_1$とします。そのサンプル数は$|S_1| =TP + FN$となります。また、実際のクラスがNeg(0)であるサンプルの集合を$S_0$とします。そのサンプル数は$|S_0| = FP + TN$となります。学習した分類モデルを$f$として、サンプル$s$に対しての予測スコアを$f(s)$で表します。

これらを用いるとROC AUCは、

$$\text{ROCAUC} = \dfrac{ \sum_{s_0\in S_0,s_1 \in S_1}{\mathbf{1}[f(s_0) < f(s_1)]} + 0.5\times\sum_{s_0\in S_0,s_1 \in S_1}{\mathbf{1}[f(s_0) = f(s_1)]}}{|S_0| \times |S_1|}$$ 

と表すことができます。*1

分母は、実際のクラスがNeg(0)のサンプル数$|S_0|$、実際のクラスがPos(1)のサンプル数$|S_1|$の全組み合わせ数$|S_0| \times |S_1|$です。

分子の第1項目は$S_0$と$S_1$の全サンプルの組み合わせの予測スコア比較で、$S_1$のサンプルの方が予測スコアが高い数 ( = 正しい順序で判別できる数 )、第2項は全サンプルの組み合わせのうち予測スコアが同じになっている数に$0.5$をかけたものです。

つまり、ROC AUCは、Pos(1), Neg(0)それぞれのクラスに属するサンプルの全組み合わせの中で、どれだけ予測スコアの順序が正しく保たれているかの指標、と見なせます。

また、この定義をよく見てみると、$|S_0|, |S_1|$に精度がとても依存していることがわかります。もしクラス比率が$(|S_0| :  |S_1|) = (8, 2)$の場合には、$S_1$の1サンプルに対応する$S_0$のサンプルは$8/2 = 4$倍あることになります。そこで分子の第一項をみると、Pos(1)を1サンプルを完全に分類できるようになったときの精度の上昇幅は、Neg(0)を1サンプルを完全に分類できるようになった場合の上昇幅の4倍となります。*2

以上のことからROC AUCは各サンプルに対して、クラス比率の逆数分だけ重みづけている指標だと考えられます。

Pythonで実装してみる

自作関数の定義

定義に従って関数を作ってみます。各クラスの全サンプルの組み合わせはitertoolsのproductを用いて実装します。そして愚直に全サンプルの予測スコアを比較して最終的にサンプル数で割る(平均をとる)実装にします。コードはこちらです。

# -*- coding: utf-8 -*-
import numpy as np
from itertools import product

def my_roc_auc(true_label: np.array, pred_proba : np.array) -> float:
    """
    ROCAUC = Pr (true_labelが0のサンプルの予測スコア < true_labelが1のサンプルの予測スコア) --- (a)
            + 0.5 * Pr (true_labelが0のサンプルの予測スコア = true_labelが1のサンプルの予測スコア) --- (b)
    from sklearn.metrcis import roc_auc_scoreの結果と一致することを確認済み
    """
    # サンプルを抽出
    pred_proba_true_label_0 = pred_proba[true_label == 0]
    pred_proba_true_label_1 = pred_proba[true_label == 1]

    # 全組み合わせを直積から作成
    pred_proba_product = product(pred_proba_true_label_0, pred_proba_true_label_1)
    pred_proba_product = np.array(list(pred_proba_product))

    # 正しい順序で分類できているフラグ (a)
    success_classified_flg = pred_proba_product[:,0] < pred_proba_product[:,1]

    # true_labelが0, 1のサンプル間で予測確率が一致しているフラグ (b)
    equal_score_flg = pred_proba_product[:,0] == pred_proba_product[:,1]

    # どちらもintにして, 足す
    roc_auc_array = success_classified_flg.astype(int) + equal_score_flg.astype(int) * 0.5
    roc_auc = np.mean(roc_auc_array)
    return roc_auc

sklearn実装との比較

自作の関数とskelarn.metricsのroc_auc_scoreを比較します。比較のためにassert_almost_equalというnumpy.testingの関数を使って比較します。assert_almost_equalは小数点第何位までの一致を判定できる関数で、一致していない場合はAssertionErrorがでます。

from numpy.testing import assert_almost_equal # 比較用の関数
from sklearn.metrics import roc_auc_score # sklearn ROC AUC
from myfunc.my_roc_auc import my_roc_auc # 自作 ROC AUC

n_samples = 1000
np.random.seed(n_samples)

# generate sample data
true_label = np.random.choice([0,1], size=n_samples)
pred_proba = np.random.choice(np.linspace(0, 1, 100), size=n_samples)

# calc ROCAUC
roc_auc_sklearn = roc_auc_score(true_label, pred_proba)
roc_auc_mine = my_roc_auc(true_label, pred_proba)

print("ROCAUC\n\tsklearn:\t%.15f\n\tmine:\t%.15f"%(roc_auc_sklearn, roc_auc_mine))
assert_almost_equal(roc_auc_sklearn, roc_auc_mine, decimal=10)

# 出力(小数点以下15桁までの一致を確認)
ROCAUC
	sklearn:	0.523924202292845
	mine:	　  0.523924202292845

結果を見てみると、小数点15位まで一致していることがわかります。

これを一応以下のコードで1000回ほどやってもAssertionErrorがでなかったので、一致していると見なします。

# 実装がsklearnのAUCと一致するかテスト, エラーがでなければ一致している

n_trials = 1000
n_samples = 200

for n in tqdm(range(n_trials)):
    np.random.seed(n)
    # generate sample data
    true_label = np.random.choice([0,1], size=n_samples)
    pred_proba = np.random.choice(np.linspace(0, 1, 100), size=n_samples)
    
    # calc ROCAUC
    roc_auc_sklearn = roc_auc_score(true_label, pred_proba)
    roc_auc_mine = my_roc_auc(true_label, pred_proba)

    # check almost equal
    assert_almost_equal(roc_auc_sklearn, roc_auc_mine, decimal=10)

1サンプル分類成功するごとの上昇幅の確認

次にクラス比率とROC AUCの関係を確認するために10サンプルのデータで検証します。クラス比率はPos(1) : Neg(0) = 2:8とします、この場合に、Posを1サンプル完全に分類できることによってNegを1サンプル完全に分類できるよりも4倍精度が向上しやすいことを確認します。

すべて0.5と予測した場合のROC AUC

# example 1
y_true_label = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0])
y_pred_proba = np.array([0.5] * 10)

# calc ROCAUC
roc_auc_sklearn = roc_auc_score(y_true_label, y_pred_proba)
roc_auc_mine = my_roc_auc(y_true_label, y_pred_proba)

print("ROCAUC\n\tsklearn:\t%.5f\n\tmine:\t%.5f"%(roc_auc_sklearn, roc_auc_mine))

全て同じ値の予測値なので、全組み合わせが分子の第2項に該当するためROC AUCが0.5となります。

# 出力
ROCAUC
	sklearn:	0.50000
	mine:	    0.50000

少数クラスであるPosを1サンプル完全に分類できるようになった場合のROC AUC

実際のラベルがPos(1)であるサンプルの予測値を1にしてROC AUCを計算してみます。Pos(1)は全体が10サンプルある中の2サンプルのみなのでクラス比率が低いクラスとなっています。

# example 2 (Posを1サンプル完全に分類できるようになった場合の精度)
y_true_label =   np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0])
y_pred_proba = np.array([0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1.0])

# calc ROCAUC
roc_auc_sklearn = roc_auc_score(y_true_label, y_pred_proba)
roc_auc_mine = my_roc_auc(y_true_label, y_pred_proba)

print("ROCAUC\n\tsklearn:\t%.5f\n\tmine:\t%.5f"%(roc_auc_sklearn, roc_auc_mine))

結果は以下のように0.75となり、0.25上昇しました。

# 出力 (Posを1サンプル完全に分類できるようになった場合の精度)
ROCAUC
	sklearn:	0.75000
	mine:	    0.75000

多数クラスであるNegを1サンプル完全に分類できるようになった場合のROC AUC

実際のラベルがNeg(0)であるサンプルの予測値を0にしてROC AUCを計算してみます。Neg(0)は全体が10サンプルある中の8サンプルなのでクラス比率が高いクラスとなっています

# example 3 (Negを1サンプル完全に分類できるようになった場合の精度)
y_true_label =   np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0])
y_pred_proba = np.array([0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.0, 0.5, 0.5])

# calc ROCAUC
roc_auc_sklearn = roc_auc_score(y_true_label, y_pred_proba)
roc_auc_mine = my_roc_auc(y_true_label, y_pred_proba)

print("ROCAUC\n\tsklearn:\t%.5f\n\tmine:\t%.5f"%(roc_auc_sklearn, roc_auc_mine))

結果は以下のように0.5625となり、0.0625の上昇となりました。Posの時の上昇幅である0.25と比べるとだいぶ小さいですね。

# 出力 (Negを1サンプル完全に分類できるようになった場合の精度)
ROCAUC
	sklearn:	0.56250
	mine:	    0.56250

1サンプル完全に分類できるようになった場合のROC AUC上昇幅の比較

上記の結果を比較してしてみます。Posの上昇幅とNegの上昇幅は、

# 1サンプルを分類成功することによるROC AUC上昇幅の比率を確認 
# (少数クラスPosの上昇幅) / (多数クラスNegの上昇幅)

(0.75000 - 0.5) / (0.56250 - 0.5)

> 4

となり、4となります。またこれは、下記のようにクラス比率の逆数と一致することがわかりました。

# 不均衡データの比率を計算する
n_major = np.sum(y_true_label == 0)
n_minor = np.sum(y_true_label == 1)
minor_ratio = n_minor / n_major
print("少数クラスの比率 :\t\t\t %.3f"%minor_ratio)
print("少数クラスの比率の逆数 :\t %.3f"%(1/minor_ratio))

少数クラスの比率 :		0.250
少数クラスの比率の逆数 :	 4.000

上記から1サンプルを完全に分類できるようになった場合のROC AUCの上昇幅は、クラス比率の逆数と一致しました。上記の実験コードは全てこちらです。

まとめ

ROC AUCが各クラスのサンプルの全組み合わせで予測スコアの順序が保たれている度合いで表せるという直感的に理解しやすい定義を紹介しました。またROC AUCがクラス比率の逆数によって重み付けされているのは、驚きでした。超極端な不均衡データ(1000:1とか) に対して ROC AUCは、その逆数分だけ重みづけてしまうので、その少数クラスのスコアの変動によってROC AUCが高くも低くもブレやすいのは納得です。ただ逆数で補正している分、各クラスを公平にみているとも取れますね。

ちなみに、自作実装のROC AUCとsklearnは結果は一致したものの、おそらく計算の仕方が異なっていて、sklearnの方がめちゃめちゃ早いです。

*1:Pepe, M. S., Longton, G., & Janes, H. (2009). Estimation and comparison of receiver operating characteristic curves. The Stata Journal, 9(1), 1-16.を参考に筆者作成

*2:ここで完全に分類できるようになるとは、そのサンプルが所属するクラスの予測スコアのうち最も極端な値をとるようになる (ROC AUCの分子の第1項の不等号が完全に成り立つようになる) ことを表しています。つまり、Posである確率が全サンプルのPosである予測スコアの最大値より大きくなる、Negの場合、Posである予測確率の最小値以下となることを指しています。

Pythonによる因果グラフ推定 -causalnexの紹介 その2-

Image by Pixabay

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

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

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

準備

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

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によって、正しく因果グラフの矢印の向きが推定されていることがわかります。正則化がうまく働いていますね。

# 出力 (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))

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

# 出力 (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)

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

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は各クラスの所属確率になっています。

また、因果グラフは正しく推定できていることがわかります。

今回は設定しませんでしたが、DAGRegressor, DAGClassifierもpytorch版のfrom_pandas同様にhidden_layer_unitsをlistで指定すると中でニューラルネットを用いて、NOTEARS推定することができます。

結構時間がかかる印象ですし、係数の意味があんまりわからなくなるので、個人的にはあまり使い所は難しいかなとは思いましたが、DAG構造が推定しつつ、とにかく精度を上げたい場合などは良いかもしれないですね。

まとめ

今回はより前回紹介しきれなかった、causalnexの機能を紹介していきました。sklearnライクに使えるのは使い勝手的に良いですね。今回使用したコードはこちらに公開しています。

*1:リンク先のNotearsMLPというクラスがpytorch版のfrom_pandas, from_numpyで使われていて、そこでMLPが設定されています。

causalnex/core.py at develop · quantumblacklabs/causalnex · GitHub

PythonでNOTEARS・ベイジアンネットによる因果グラフ推定 -causalnexの紹介-

Image by Pixabay

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

因果探索は、データを与えることで、そのデータの変数間に潜む因果構造を推定しようという手法です。メジャーな手法としては、離散変数に対してベイジアンネットワークや、非ガウス連続変数に対してのLiNGAMなどの手法があります。

その中でも機械学習系のトップカンファレンスNeurIPSの2018年論文「DAGs with NO TEARS: Continuous Optimization from Structure Learning」で提案されたNOTEARSという手法について紹介したいです！NOTEARSは、実質損失関数を加えるのみで、今までの機械学習ライクなシンプルなアプローチで因果グラフを推定できる手法なので紹介してみます。

また、Python実装がNOTEARS実装がcausalnexでQuantumBlackから公開されているので使用方法も含め紹介していきます。

DAGと因果グラフ

DAGとは

ベイジアンネットやLiNGAMといった因果探索手法では、推定する因果グラフにDAG構造を想定することが多いです。NOTEARSも例外ではなく、推定される因果グラフにDAGを想定しています。

DAGとはDirect Acyclic Graph、日本語では有向非巡回グラフというもので、つまりは矢印がついていて巡回していないグラフを指します。

例を出してみます。下の左図は、有向ではあるものの、巡回している（どこか任意のノードから見て自分に返ってくるパスがある）ため、DAGではないです。右図に関しては、有向かつ巡回していないため、DAGの例となります。

 因果グラフの意味合い

誤解を恐れずにいうと、因果グラフは、NOTEARSやLiNGAMの場合は、データ生成過程とみなすとわかりやすいかなと思います。例えば、上の右図の因果グラフで、$x_1$に関しては、$x_1 = f_{31}(x_3) + e_1$という風にデータが生成され、$x_2$に関しては、$x_2 = f_{12}(x_1) + f_{32} (x_3) + e_2$というように$x_1, x_3$を基にデータが生成されるとみなします。

では、因果グラフにDAGを仮定するのはどういうことかというと、もちろん制約として入れることで問題を解きやすくするというのもあるのでしょうが、データが巡回して生成されないということなので、自然な仮定なのかなと思われます。

NOTEARSとは

以下ではNOTEARSの基本的な考え方を紹介します。NO TEARS はNon-combinatorial Optimization via Trace Exponential and Augmented lagRangian for Structure learning の略です。

手法名めちゃめちゃ長いですが、NOTEARSの重要なポイントとしては、

1. DAGであることを定量的に評価。
2. それを最適化問題の制約とすることで、組み合わせ爆発を避けてDAGを探索することが可能になる。

の2点かなと思います。

そもそもDAGを、データの生成過程と捉えると、データ$X$は、下図のように

$X = WX$のように表すことが可能になります。この図では、$W$について一旦線形を仮定してます。また、$W$の対角成分はDAGの制約（巡回なし）から数値が入らないです。

NOTEARSではこの$W$についてDAGであることを定量的に表すことで、最適化問題の制約として、それを加えることで、因果グラフ(DAG)を推定するということを考えます。

結論をいうと、DAGであることの定量評価(DAG-ness)としては、

$$h(W) = \text{tr}\left( e^{W\circ W} \right) -d = 0$$ 

が提案されています。ここで$W$は$d\times d$の重み行列で、$\circ$はアダマール積（各要素を対応する要素で積をとるもの)になります。$W$がDAGである場合には、$h(W) = 0$となります。つまり、$h(W)$が0になるような$W$を推定できればそれはDAGであることになります。*1

最終的には、これを制約として加えて、$W$を以下の問題を推定することを考えます。

 \begin{align}& \min_{W\in\mathbb{R}^{d\times d}} 　&   \dfrac{1}{2n}\| X -  WX \|_{F}^2 + \lambda \| W\|_1\\\ & \text{subject to }   &  \text{tr}\left( e^{W\circ W}  \right)-d = 0\end{align}

この式は、単純に$X$の推定式である$WX$と$W$の2乗誤差（と$\ell_1-$penalty）となっており、制約として$W$がDAGであることが加わっています。単純にそれぞれの変数に対して他の変数からのみ（自身を含めずに）回帰する、かつ制約としてその前回の関係式はDAGであるという最適化問題になります。

これによって推定される$W$の意味合いとしては、そもそも$W$はデータ生成仮定を表す行列である上、DAGであることを制約として加えているため、この$W$が因果グラフとなります。

式自体は線形の式を紹介しましたが、もちろんこれは非線形にも$WX$を変えるだけで容易に対応できます。causalnexでもニューラルネットを用いて非線形にしたバージョンなどが実装されています。

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

ここからcausalnexを使ってNOTEARSで因果グラフを、その推定結果を元にベイジアンネットを構築し、予測まで実行していきます。コードはこちらで公開しています。 

サンプルデータ作成

実験用のサンプルデータを生成します。

サンプルデータとしては、以下の因果グラフを従う5変数($y, x_1, x_2, c_1, c_2$)を作成して、それだけを入力とした時に正しく因果グラフが推定できるかどうかを確認していきます。

def generate_sample_data(n_samples=1000):
    """
    サンプルデータ生成
    | c1 --> x1
    | c2 --> x1
    | c1 --> y
    | c2 --> y
    | x1 --> y
    | x2 は独立
    """
    np.random.seed(0)
    c1 = np.random.normal(0, 0.5, size=n_samples)
    c2 = np.random.choice(2, size=n_samples)

    np.random.seed(1)
    x1 = np.random.normal(1, 2, size=n_samples) + c1 + c2/10

    np.random.seed(2)
    y = np.random.uniform(-1,2, size=n_samples) + x1 + c1 + c2 /10

    np.random.seed(3)
    x2 = np.random.uniform(-2, 2, size=n_samples)
    raw_data = pd.DataFrame({"x1": x1, "y":y, "x2": x2, "c1": c1, "c2":c2})
    return raw_data

# データ作成
n_samples = 2000
struct_data = generate_sample_data(n_samples)

上の因果グラフに従う2000サンプルのデータを作成しました。

結果を表示するための関数を準備

causalnexでは、plotツールが準備されているのですが、pygraphvizなど環境をつくるのが面倒なインストールする必要があります。それを避けるために自作のプロット関数を作っておきます。

causalnexで推定した因果グラフはnetworkxを元に作成されているため、networkxのグラフから因果グラフをplotするutilを作りました。こちらも同様にこちらのリンクに入っています (util/networkx_util.py)。 

# -*- coding: utf-8 -*-
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

class NetworkxUtil:
    def plot_structure_model(self, structure_model, layout_method="spring", layout_seed=1, figsize=(5,5), 
                              node_shape="o", node_size=1500, node_color="#1abc9c", font_size=18, 
                             edge_color="#34495e", min_edge_width=2, max_edge_width =3.5, arrowsize=25, 
                              plot_edge_weights=True, edge_weights_color="#e74c3c", edge_weights_fontsize=12, alpha=0.8):
      # plot settings
      fig, ax = plt.subplots(figsize=figsize)
      pos = self._load_layout(structure_model=structure_model, layout_method=layout_method, layout_seed=layout_seed)

      # edgeの重み係数をplotするか否かのフラグ
      if plot_edge_weights:
        nx.draw_networkx_edge_labels(
            structure_model, pos,
            edge_labels={(u, v): round(d["weight"], 4) for (u,v,d) in structure_model.edges(data=True)},
            font_color=edge_weights_color,
            font_size = edge_weights_fontsize
        )

      # edgeの重みに応じて太さを変更、最低でもmin_edge_widthに設定
      edge_width = [
                    np.min([np.max([d["weight"], min_edge_width]), max_edge_width] )
                    for (u, v, d) in structure_model.edges(data=True)
      ]
      # networkxでの推定したモデルのplot
      nx.draw_networkx(
          structure_model,
          ax = ax,
          pos = pos,
          node_shape = node_shape,
          node_color = node_color,
          node_size = node_size,
          edge_color = edge_color,
          width = edge_width,
          alpha = alpha, 
          font_size = font_size,
          arrowsize = arrowsize,
          with_labels = True,
          # connectionstyle="arc3, rad=0.1" # if curve
      )
      plt.axis("off")
      fig.show()

    def _load_layout(self, structure_model, layout_method, layout_seed=0):
      """layout_methodを指定して、pos: position keyed by nodeを取得"""
      if layout_method == "circular":
        pos = nx.circular_layout(structure_model)
      elif layout_method == "spring":
        pos = nx.spring_layout(structure_model, seed=layout_seed)
      elif layout_method == "planar":
        pos = nx.planar_layout(structure_model)
      elif layout_method == "shell":
        pos = nx.shell_layout(structure_model)
      elif layout_method == "random":
        pos = nx.random_layout(structure_model, seed=layout_seed)
      else:
        raise ValueError(f"Method {layout_method} is not expected.")
      return pos

    @staticmethod
    def print_weights(structure_model, cutoff=0):
      """エッジとその重みをprint"""
      for (u, v, d) in structure_model.edges(data=True):
        if np.abs(d["weight"]) >= cutoff:
          print(f"[ {u} ] ---> [ {v} ]\t\tWeight\t{d['weight']:.5f}")
  

NOTEARSで因果グラフ推定

では早速、NOTEARSを実行していきます。 NOTEARSはfrom_pandasというmethodでデータを入力することで簡単に実行できます。

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

# NOTEARSを実行, from_pandas_lassoでL1 penaltyをつけて推定することが可能
sm = from_pandas( 
    struct_data,
    tabu_edges = [],
    tabu_parent_nodes = None,
    tabu_child_nodes = None,
)
sm.threshold_till_dag() # DAGになるように閾値をあげる
sm.remove_edges_below_threshold(0.2) # 係数の閾値を0.2にする
nx_util.plot_structure_model(structure_model=sm, layout_method="circular") # この条件で正しく推定できていることがわかる。

from_pandasの結果のみからでは誤差の影響で完全なDAGにならないケースがあります。そのためthresholdをDAGになるまで上げる(threshold_till_dag)や、小さすぎる重みを持つエッジを削除する(remove_edges_below_threshold)などある程度操作する必要があります。

そうして出たプロット結果をみてみると、ただ正しく因果グラフが推定されていることがわかります。

また今回の$x2$のように他のノードとの関連がないノードが含まれてしまっている場合など、因果グラフをきれいにするため最大のサブグラフに限定して抽出することもできます。

# largest subgraphを抽出
sm = sm.get_largest_subgraph()
nx_util.plot_structure_model(sm, layout_method="random", figsize=(9,6))

# 係数をprint
nx_util.print_weights(sm)
 

# 出力
[ x1 ] ---> [ y ]		Weight	1.05582
[ c1 ] ---> [ x1 ]		Weight	0.93608
[ c1 ] ---> [ y ]		Weight	0.91724
[ c2 ] ---> [ x1 ]		Weight	1.05973
[ c2 ] ---> [ y ]		Weight	0.49489

以上のようにNOTEARSによって正しく因果グラフ（ここでは、データ生成過程）を推定できることがわかりました。

NOTEARSでTabuを指定して因果グラフを推定

NOTEARSでは、あらかじめ「この矢印はないだろう」といった因果グラフに関しての人間の知見をtabu_edgesという形で取り込むことができます。

他にも、このノードは親ノードにはなり得ない、子ノードにはなり得ないといった知見を加えることができるので、実際にやってみて因果グラフを推定してみます。

# struct data
sm_with_tabu = from_pandas(
    struct_data,
    tabu_edges = [("c1", "x1")], # c1 --> x1 に線が引かれなくなる。
    tabu_parent_nodes = ["c2"], # 親ノードにならなくなる。矢印元にならなくなる。
    tabu_child_nodes = ["y"], # 子ノードにならなくなる。矢印が刺されなくなる。
)
sm_with_tabu.threshold_till_dag()
sm_with_tabu.remove_edges_below_threshold(0.2)

nx_util.plot_structure_model(structure_model=sm_with_tabu, layout_method="circular") # 
nx_util.print_weights(sm_with_tabu) # 係数を print

プロット結果を見ると、指定した$c_1$から$y$までのパスはtabu_edgesを指定するまではあったのですが、消えていることがわかります。また他にもtabu_child_nodesで指定した$y$が矢印がささっておらず子ノードになっていないこと、tabu_parents_nodesで指定した$c2$から矢印が出ておらず親ノードになっていないことが確認できます。

以下ではtabuを指定せず、正しく因果グラフを推定できたグラフを元にベイジアンネットを構築します。

推定した因果グラフ構造を元にベイジアンネットを推定

causalnexではNOTEARSで推定した因果グラフを元にベイジアンネットワーク(pgmpy実装を内部で利用している)を推定することができます。

ただし、ベイジアンネットワークは離散化した変数に対して適用する必要があるため、NOTEARSを推定したデータそのものではなく、適宜離散化する必要がありますcausalnexではDiscretiserという離散化のための便利ツールが用意されているので、今回はこれを使っていきます。

from causalnex.discretiser import Discretiser

# ベイジアンネットを推定するために離散化データを作成する。
discretised_data = struct_data.copy().drop(["x2"], axis=1)

# quantileでnum_buckets数に分割。
discretised_data["y"] = Discretiser(method="quantile", num_buckets=5).fit_transform(discretised_data["y"].values) 
# ほぼ同数num_buckets数に分割。
discretised_data["x1"] = Discretiser(method="uniform", num_buckets=5).fit_transform(discretised_data["x1"].values)
# percentile_split_pointsを指定して分割。それぞれ<=が含まれる集合になる
discretised_data["c1"] = Discretiser(method="percentiles", percentile_split_points=[0.25, 0.5, 0.75]).fit_transform(discretised_data["c1"].values)
# split_pointsで分割。それぞれ<=が含まれる集合になる
discretised_data["c2"] = Discretiser(method="fixed", numeric_split_points=[0,1]).transform(discretised_data["c2"].values) 

NOTEARSで推定した因果グラフと、離散化したデータを用いてベイジアンネットワークを推定します。 推定したfit_cpdsによって条件付き確率を計算することができます。計算した後は、対象の子ノードを指定してbn.cpds["y"]といった形で条件付き確率表を出力することができます。

from causalnex.network import BayesianNetwork
# load estimated causal graph
bn = BayesianNetwork(sm)

bn = bn.fit_node_states(discretised_data) # 取りうるすべての状態を指定、 本来testデータに新出のカテゴリがある場合には判定し、変換する必要がある
bn = bn.fit_cpds(discretised_data, method="BayesianEstimator", bayes_prior="K2")

ベイジアンネットでの推論

 推定したベイジアンネットワークでは、ラベルの予測(predict), そのラベルへの予測所属確率(predict_probability)を推定することができます。今回予測対象を$y$にしますが、推定されたデータ生成過程では$x2$などは$y$には不要なので、必要なカラムのみをスライスして予測してみます。

pred_target_col = "y"
parents_columns = bn.cpds[pred_target_col].columns.names # 予測対象カラムの親ノードのみ必要であるためsliceする、sliceせずとも実行可能
y_pred = bn.predict(discretised_data[parents_columns], pred_target_col) # ラベル予測
y_predict_proba = bn.predict_probability(discretised_data[parents_columns], pred_target_col) # 予測確率

# 予測結果
print("True Label :\t", discretised_data["y"].values[:5])
print("Pred Label :\t", np.ravel(y_pred)[:5])
print("Pred Proba [y=4] :\t", y_predict_proba["y_4"].values[:5])

# 出力結果
True Label :	 [4 0 2 2 4]
Pred Label :	 [4 1 3 1 4]
Pred Proba [y=4] :	 [0.88235294 0.01298701 0.30075188 0.01639344 0.85074627]

予測結果がでてきますね。こうして、NOTEARSで推定した因果グラフを元に、ベイジアンネットワークを推定し、予測まで行うことができました。

フルのコードはこちらにおいています。 

まとめ

ベイジアンネットやLiNGAMなどの因果探索手法は理論面が難しいなーと思っていたのですがNOTEARSのように実質損失関数を追加するだけで、因果グラフを候補を出すことができるのは、説明もしやすいしとても便利ですね。また、tabu_edgeなどで人の知識を取り入れることができるので、より納得感のあるモデルを作りやすそうですね。

人の知識をいれられるようなモデルは説明しやすい面白いので、今後も積極的にキャッチアップできればと思います。causalnex内でもまだまだ色んな機能があるので、次回紹介できればと思います。

今回は詳細を端折ってしまいましたが、そもそも因果・相関って何が違うのといった話や因果推論のイメージについてはこちらの本がわかりやすいので超おすすめです。

「原因と結果」の経済学―――データから真実を見抜く思考法

他にも同じ因果探索手法であるLiNGAMについて、LiNGAMの発明者である清水先生の本がとても行間が少なくわかりやすい本なのでオススメです。

統計的因果探索 (機械学習プロフェッショナルシリーズ)

*1:詳細のアルゴリズムや証明は元論文をご参照ください。

https://papers.nips.cc/paper/2018/file/e347c51419ffb23ca3fd5050202f9c3d-Paper.pdf

機械学習の公平性を担保する手法紹介 -Pythonパッケージfairlearnをつかった実装-

Image by Pixabay

今回の記事は、機械学習の公平性 (Fairness)についてです。

機械学習の公平性は、ざっくりいうと、AIによって差別を生まない・助長しないためにAIを公平に作れないかといったことです。例えば、予測結果がセンシティブな変数（状況やタスクによるものだとは思いますが、人種や性別など）に左右されないようにしたいなど。

ただそれらはセンシティブな変数に関するデータを集めない or モデルに入力しないといったことだけでは実現できない場合が多いです。

そもそもの公平性の定義やそれを実現するアルゴリズムとして様々なものがあるのですが、今回はPythonパッケージであるfairlearnに実装されているアルゴリズムを中心に紹介していきます。

そもそも公平性って？

定義

そもそも公平性って何だって話ですが、いろんな定義があるようです。サーベイ論文とかみると10個くらいあるんですが*1*2、（集団）公平性の定義メジャーどころとしてDemographic Parity、Equalized Oddsがあります。

2値分類を例に挙げて、$Y$が真のラベル、$S$がセンシティブな変数、$\hat{Y}$を予測ラベルで、それぞれ$\left\{0,1\right\}$だとすると、

1. Demographic parity (Statistical Parity)
   $$ \text{Pr}\left\{  \hat{Y} = 1 | S = 0 \right\} = \text{Pr}\left\{  \hat{Y} = 1 | S = 1 \right\}$$
   
2. Equalized odds
   $$ \text{Pr}\left\{  \hat{Y} = 0 | S = 0, Y=y \right\} = \text{Pr}\left\{  \hat{Y} = 1 | S = 1, Y=y \right\}, y \in \left\{0,1\right\}$$

と表されます。

ちなみに、機械学習では$y=1$の場合に、採用、合格、昇進など"advanced"なアウトカムを持ってくることが多いことから、Equalized Oddsで$y=1$のみのケースに限定したものは、Equal Opportunityと言われるそうです *3。 

定義の違いと解釈

定義をそのままみると、Demographic Parityは、センシティブな変数で条件づけられた予測ラベルの分布が一致する。Equalized oddsは、センシティブな変数に加え、真のラベルで条件づけられた予測ラベルの分布が一致するという定義になっています。

また、よくみてみると、Equalized oddsは、$y=1$の場合は$S$ごとのTrue Positive Rateの差が0、$y=0$の場合はFalse Positive Rateの差が0、と読むことができ、センシティブ変数によって精度（TPR、FPR）が変わらないことを要求していることがわかります。

これらの公平性の解釈を自分なりにしてみると、

Demographic Parityは、センシティブ変数による予測確率に差がないものの、本来の能力や適性といった$\text{Pr} (Y = 1)$や$\text{Pr} (Y = 0)$がセンシティブ変数によって本来偏っているべきものの場合、逆差別を生んでしまう可能性がある。

ex. 採用条件としては本来女性の方が向いているのに、男性を採用するようにした。

Equalized oddsは、センシティブ変数によってTPR、TPRに差がなく同精度を要求しているものの、本来の$\text{Pr} (Y = 1)$や$\text{Pr} (Y = 0)$が差別を生んでしまっている場合にそれを許容してしまいます。

ex. 採用条件としては本来女性の方が向いているので、女性を採用したいがそれは差別を生んでしまっている。

といったことが言えるのではないかなと。

どうやって公平性を実現するか

ではこのDemographic ParityやEqualized Oddsといった公平性をどうやって実現するのかの話です。単純にセンシティブな変数をモデルに食わせなきゃいいじゃん！と思いますが、その場合センシティブな変数と、他のモデルに入れる変数に関連があると、間接的にセンシティブな変数の情報を含んでしまうので公平な予測結果を実現できません。

ではどうやるのかというと、これには大きく分けて3つの手法の枠組みがあります。

1. Pre-Processing
   モデルにfitする前に前処理としてセンシティブ変数に関連するような情報を除去しようとするアプローチ。モデルは既存のものを想定。例えば、センシティブ変数でそれ以外を回帰して、その結果を引いたデータをつくるなど。
   
   
2. In-Processing
   フィットするモデルを修正して公平性をアルゴリズム内で担保しようとするアプローチ。例えば、モデルの最適化の際に制約を加えるなど。
   
   
3. Post-Processing
   モデルは既存のものを使ってでてきた予測分布を補正しようとするアプローチ。例えば、センシティブ変数でグループ分けして、公平性の基準で閾値をグループごとに決定し、予測ラベルを決めるなど。

fairlearnで実装されているメソッド紹介

 fairlearnで実装されているものの中で以下の3つのメソッドを紹介してきます。

1. preprocess.CorrelationRemover（Pre-Processing）
   

   このメソッドでは単純にセンシティブ変数で、それ以外の変数を推定する線形回帰を行い、元のデータセットから予測結果引くことで、センシティブ変数とはそれぞれ無相関なデータセットを得る方法です。
   
   元論文が見つけられなかったのですが、fairlearnのサイトと実装をみると、センシティブ変数群$S$を中心化$S_\text{center} = S - \bar{S}$して、それを元にセンシティブ変数以外の変数群$X_\text{orig}$を回帰し、係数$W$の推定値
   $$\hat{W} = \text{arg}\min_W \| X_\text{orig} - WS_\text{center} \| $$を得ます。そうして得た$\hat{W}$を用いて、
   $$X_\text{filtered} =X_\text{orig}  - \hat{W}S_\text{center} $$を推定。最終的に、どれだけリークさせるかの調整パラメータ$\alpha$を用いて、
   $$X_\text{tfm}  = \alpha X_\text{filtered} - (1 - \alpha) X_\text{orig}$$が返ってきます。
   
   この時、センシティブ変数は落ちて返ってくるので変数数は減って返ってくることに注意です。また、公式ドキュメントにも書いてありますが、無相関は独立とは別なので、公平性の文脈では（一般化）線形回帰などの線形関係を扱うモデルの前処理として使うべき手法だと思われます。
   
   他のメソッドと違い、明示的にDemographic Parity, Equalized Oddsを実現するわけではないですが、$\alpha = 1$の場合でかつ線形回帰の前処理として用いるとしたら、完全にセンシティブ変数の情報が入力データからは消えるので、Demographic Parityは保てるのではないかなと思いました。(間違っていたらごめんなさい）

2. reduction.DemographicParity, EqualizedOdds, TruePositiveRateParity（In-Processing）
   元論文はこちらで提案手法はReduction Approachといわれています。ざっくり手法の方針を書くと、通常の分類器が解く最適化問題にそれぞれDemographicParity, EqualizedOdds, EqualOpportunity(fairlearnではTruePositiveRateParity)を制約として入れることを考えます。
   
   この定式化の際に、ラグランジュの未定乗数法を使うため、「通常のloss」 + $\lambda$「公平性のloss」を最小化する問題となります。この$\lambda$はハイパーパラメータなので、指定する必要があるのですが、fairlearn.reduction.GridSearchにて$\lambda$を探索することができます。
   
   緩和した問題を解くため、それぞれの性質を保証できるわけではないことに注意が必要です。LightGBMなどポピュラーな分類器を元のestimatorとして指定して使用することができるものの、sample_weight（とy）を用いて制約を表現するため、fit時に指定できるモデルでなければ現状は使えない点に注意が必要です。
   
   この手法の強みとして、学習時にはセンシティブ変数が必要だが、予測時（test-time access）には必要ないことが挙げられます。今まではセンシティブなデータを集めていたけども、今後はセンシティブデータを集めないようにしたいなど場面で活躍しそうですね。
   
   
3. postprocess.ThresholdOptimizer（Post-Processing）
   元論文はこちらです。
   
   ざっくり手法の方針を書くと、何らかの方法で学習済みの2値分類モデルによる予測ラベルを変換して、Demographic Parity、Equalized OddsやEqual Opportunityなど指定した公平性を満たすような"良い"閾値を各群で見つけることになります。既にある学習済みのモデルの予測ラベルを変換することで、それを達成しようとします。
   

   学習済みの2値分類モデルによる予測ラベルを$\hat{Y}$として、それをpost-processingにより変換した予測$\tilde{Y}$、センシティブ変数$S\in \left\{0,1\right\}$だとします。この予測$\tilde{Y}$を求めるのがpost-processingの目的になります。その可動域から、それぞれの求めたい制約に一致する点を設定し、解を求めることになります。

   

   元論文(https://arxiv.org/abs/1610.02413), fig1

   
   元論文の図がわかりやすいのですが、TPr, FPrに対して元の予測$\hat{Y}$を変化させることで達成しうる可動域がそれぞれのセンシティブ変数の値で緑と青とで分かれています。このときEqualized Oddsを達成するのは、$x, y$が一致する点（TPr, FPrが一致する点がEqualOddsなので）となります。$x, y$が一致するのはOverlap全体になるわけなのですが、その中でもロスが最小化される点（左図の赤点など）に対して、各群でその点を達成する点を算出するのが流れになります。
   
   ただし、これらを最小化問題として解くため、こちらのメソッドに関しても、それぞれの性質を完璧に保証できるわけではないことに注意が必要です。
   
   こちらもLightGBMなどポピュラーな分類器の学習結果を使用することができます。記事作成の2021/5月時点でセンシティブ変数が1つ、2値分類のみに対応しています。

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

前置きが長くなりましたが、早速2値分類の問題に対してfairlearnを使ってみます。コードはすべてこちらに置いています。

公平性に関するMetricsの設定 

fairlearnでは公平性に関する様々なmetricsが用意されているようです。

今回の記事ではDemographicParityとEqualizedOddsに注目したためそれに関連する、以下の4つのmetricsを使用するのですが、理解を深める意味で自分でも実装してみようとおもいます！

• demographic_parity_difference
• true_positive_parity_difference
• false_positive_parity_difference
• equalized_odds_difference

 まずはfairlearnではfairlearn.metrics配下に各metricsが用意されているのでそれをimportとします。

# metrics import
from fairlearn.metrics import (
    demographic_parity_difference,
    true_positive_rate_difference, 
    false_positive_rate_difference,
    equalized_odds_difference)

# テストデータを準備
y_true = np.array([0,0,1,1,1,0,0])
y_pred = np.array([0,0,0,1,1,1,1])
s = np.array([1,0,1,0,1,0,0])
    
# fairlearnでのmetrics実装, 自作funcとの[abs_diff]との一致が確認できる
print("Demographic Parity Diff :\t\t %.3f"% demographic_parity_difference(y_true=y_true, y_pred=y_pred, sensitive_features=s))
print("TruePositiveRate Parity Diff :\t %.3f"% true_positive_rate_difference(y_true=y_true, y_pred=y_pred, sensitive_features=s))
print("FalsePositiveRate Parity Diff :\t %.3f"% false_positive_rate_difference(y_true=y_true, y_pred=y_pred, sensitive_features=s))

# TruePositiveRate DiffとFalsePositiveRate Diffのうち大きい方の値が入る
print("EqualizedOdds Parity Diff :\t\t %.3f"% equalized_odds_difference(y_true=y_true, y_pred=y_pred, sensitive_features=s))

# print結果
Demographic Parity Diff :		 0.417
TruePositiveRate Parity Diff :	 0.500
FalsePositiveRate Parity Diff :	 0.667
EqualizedOdds Parity Diff :		 0.667

fairlearnでの実行結果が確認できたので、今度は自分で実装してみます、それぞれの公平性の定義を元に各センシティブ変数の群で差をとったものを実装します。

# 公平性の各種metricsをprint
def print_result_summary(y_true, y_pred, s):
    # Demographic Parityで用いる差
    dp_ave =  y_pred.mean()
    dp_s1 =  y_pred[s==1].mean()
    dp_s0 = y_pred[s==0].mean()
    dp_diff = np.abs(dp_s1 - dp_s0)
    
    # TruePositiveRate Parityで用いる差 (Equal Opportunity)
    tpr_ave = y_pred[y_true==1].mean()
    tpr_s1 = y_pred[np.all([y_true==1, s==1],axis=0)].mean()
    tpr_s0 = y_pred[np.all([y_true==1, s==0],axis=0)].mean()
    tpr_diff = np.abs(tpr_s1 - tpr_s0)
    
    # FalsePositiveRate Parityで用いる差
    fpr_ave = y_pred[y_true==0].mean()
    fpr_s1 = y_pred[np.all([y_true==0, s==1],axis=0)].mean()
    fpr_s0 = y_pred[np.all([y_true==0, s==0],axis=0)].mean()
    fpr_diff = np.abs(fpr_s1 - fpr_s0)
    
    # print result
    dp_text = f"Demographic Parity:\t\t[mean] {dp_ave:.3f},\t[s=1] {dp_s1:.3f},\t[s=0] {dp_s0:.3f},\t[abs_diff] {dp_diff:.3f}"
    tpr_text = f"TruePositiveRate Parity:\t[mean] {tpr_ave:.3f},\t[s=1] {tpr_s1:.3f},\t[s=0] {tpr_s0:.3f},\t[abs_diff] {tpr_diff:.3f}"
    fpr_text = f"FalsePositiveRate Parity:\t[mean] {fpr_ave:.3f},\t[s=1] {fpr_s1:.3f},\t[s=0] {fpr_s0:.3f},\t[abs_diff] {fpr_diff:.3f}"
    print(dp_text)
    print(tpr_text)
    print(fpr_text)

実行してみます。

# 自作metricsのprint
print_result_summary(y_true, y_pred, s)

# 結果
Demographic Parity:		[mean] 0.571,	[s=1] 0.333,	[s=0] 0.750,	[abs_diff] 0.417
TruePositiveRate Parity:	[mean] 0.667,	[s=1] 0.500,	[s=0] 1.000,	[abs_diff] 0.500
FalsePositiveRate Parity:	[mean] 0.500,	[s=1] 0.000,	[s=0] 0.667,	[abs_diff] 0.667

 結果を見るとabs_diffの結果がfairlearn実装のmetricsと一致していることがわかります。equalized_odds_differenceはtpr, fpr_differenceのうち大きい方をとっているため少々トリッキーですが、それ以外は単純に定義に従って差を取っているだけなので、どれだけその公平性の定義に従っているかの定量評価指標になっていることがわかります。

サンプルデータ作成

これから上で紹介した3つのメソッドを実行するのですが、そのための実験用の擬似データを作成しておきます。作成するのは3変数で$x, y, s$を作成します。$s$はセンシティブ変数として、$x$は$s$と関連がある特徴量です。そして$y$はラベルで$x, s$の両方とも関連があるようにします。

import numpy as np
def logistic(x, a=1):
    return 1/(1 + np.exp(-a*x))

# サンプルデータ生成
def make_sample_data(N=10000, p_s=0.1):
    # p_s: sensitive変数の平均値
    s = np.r_[np.zeros(int(N*p_s)) , np.ones(int(N*(1-p_s)))] # sensitive
    np.random.seed(0)
    x = np.random.normal(0, 1, size=N) - 0.2*s  # correlated with s

    np.random.seed(1)
    y = 0.3*x + 0.5*s + np.random.normal(0,1,size=N)  # outcome y is correlated with x and s
    y -= np.mean(y)
    y = np.array(logistic(y) > 0.5, dtype=int) # flg

    np.random.seed(2)
    train_idx = np.random.choice(N, size=N//2, replace=False) # random split
    test_idx = np.array(list(set(np.arange(N)) - set(train_idx)))

    return x[train_idx], y[train_idx], s[train_idx], x[test_idx], y[test_idx], s[test_idx]
    
# sample data
x_train, y_train, s_train, x_test, y_test, s_test = make_sample_data(N=10000, p_s=0.1)
X_train = np.c_[x_train, s_train]
X_test = np.c_[x_test, s_test]

 上記のコードでサンプル数10000のデータを作成しました。以下ではこれを用いてfitしていきます。 

ベースラインの実行

fairlearnで公平性を実現する前に、単純に何も処理せずモデルを学習した時の各metricsの値を見てみます。LGBMClassifierを用いてfitしてみます。

# model fit
baseline_model = LGBMClassifier(random_state=1)
baseline_model.fit(X_train, y_train)

# predict
y_train_pred_baseline = baseline_model.predict(X_train)
y_test_pred_baseline = baseline_model.predict(X_test)

# print result
print("Train")
print_result_summary(y_true=y_train, y_pred=y_train_pred_baseline, s=s_train)
print("\nTest")
print_result_summary(y_true=y_test, y_pred=y_test_pred_baseline, s=s_test)

# 結果 （ベースラインモデル）
Train
Demographic Parity:		[mean] 0.548,	[s=1] 0.575,	[s=0] 0.304,	[abs_diff] 0.271
TruePositiveRate Parity:	[mean] 0.681,	[s=1] 0.689,	[s=0] 0.580,	[abs_diff] 0.109
FalsePositiveRate Parity:	[mean] 0.411,	[s=1] 0.451,	[s=0] 0.130,	[abs_diff] 0.321

Test
Demographic Parity:		[mean] 0.529,	[s=1] 0.558,	[s=0] 0.270,	[abs_diff] 0.288
TruePositiveRate Parity:	[mean] 0.607,	[s=1] 0.624,	[s=0] 0.372,	[abs_diff] 0.252
FalsePositiveRate Parity:	[mean] 0.451,	[s=1] 0.486,	[s=0] 0.216,	[abs_diff] 0.270

Demographic Parity, TPR, FPR Parityの全ての各群で大きな差が出てしまっていることがわかります。次に、これをセンシティブ変数sをdropした状態でfitした結果を見てみると、

# 結果 （ベースラインモデル、センシティブ変数drop）
Train
Demographic Parity:		[mean] 0.563,	[s=1] 0.558,	[s=0] 0.608,	[abs_diff] 0.050
TruePositiveRate Parity:	[mean] 0.682,	[s=1] 0.671,	[s=0] 0.824,	[abs_diff] 0.153
FalsePositiveRate Parity:	[mean] 0.440,	[s=1] 0.435,	[s=0] 0.472,	[abs_diff] 0.037

Test
Demographic Parity:		[mean] 0.543,	[s=1] 0.539,	[s=0] 0.578,	[abs_diff] 0.039
TruePositiveRate Parity:	[mean] 0.608,	[s=1] 0.604,	[s=0] 0.669,	[abs_diff] 0.065
FalsePositiveRate Parity:	[mean] 0.476,	[s=1] 0.468,	[s=0] 0.530,	[abs_diff] 0.062

となり、センシティブ変数を入れていた時よりもだいぶ差は緩和されるものの、まだ最小でも5%の差が群間であることがわかります。

各アルゴリズムの実行 (CorrelationRemover)

さて、それではfairlearnを用いていきます。最初にCorrelationRemoverを用いていきます。使い方は簡単で、CorrelationRemoverの中のsensitive_feature_idsにセンシティブ変数のインデックスを指定して、あとはfit、transformを行うだけです。 

from lightgbm import LGBMClassifier
from fairlearn.preprocessing import CorrelationRemover
# 相関を除去
corr_remover = CorrelationRemover(sensitive_feature_ids=[1]) # X_trainのうち2列目がsに該当
X_train_rmcorr = corr_remover.fit_transform(X_train)
X_test_rmcorr = corr_remover.transform(X_test)

# model fit
clf = LGBMClassifier()
clf.fit(X_train_rmcorr, y_train)

# predict
y_train_pred_rmcorr = clf.predict(X_train_rmcorr)
y_test_pred_rmcorr = clf.predict(X_test_rmcorr)

# print result
print("Train")
print_result_summary(y_true=y_train, y_pred=y_train_pred_rmcorr, s=s_train)
print("\nTest")
print_result_summary(y_true=y_test, y_pred=y_test_pred_rmcorr, s=s_test)

予測値の結果を見てみると、

# 結果 （CorrelationRemover）
Train
Demographic Parity:		[mean] 0.539,	[s=1] 0.539,	[s=0] 0.540,	[abs_diff] 0.001
TruePositiveRate Parity:	[mean] 0.661,	[s=1] 0.655,	[s=0] 0.736,	[abs_diff] 0.080
FalsePositiveRate Parity:	[mean] 0.413,	[s=1] 0.412,	[s=0] 0.417,	[abs_diff] 0.005

Test
Demographic Parity:		[mean] 0.526,	[s=1] 0.524,	[s=0] 0.544,	[abs_diff] 0.020
TruePositiveRate Parity:	[mean] 0.594,	[s=1] 0.588,	[s=0] 0.674,	[abs_diff] 0.087
FalsePositiveRate Parity:	[mean] 0.458,	[s=1] 0.456,	[s=0] 0.476,	[abs_diff] 0.020

となり、どれも大幅に改善されていることがわかります。

特にDemographic Parity Differenceに関しては、trainで0.1%、testでも2%となっていることがわかります。今回はシンプルな線形データを用いたからだとは思いますが、このようにわかりやすい手法で良い結果が出るのはいいですね。

各アルゴリズムの実行 (Reduction Approach)

ReductionMethodに関しては、GridSearchを用いて実行していきます。肝となるのはGridSearchのconstraintsという引数になります。これに対してreduction.DemographicParityなどを呼び出して指定することで、指定された公平性をモデルに加味することができます。

他にもEqualizedOdds, TruePositiverateParityなど様々なものが用意されているので、それを切り替えるだけで加味する公平性を変えることができます。その他の引数として、estimatorを指定して（sample_weightをfit時に指定できるモデルに限定）、グリッドサーチのパラメータを設定することができます。

実行はsklearnライクにfit, predictで行うことができます。

# model fit
from fairlearn.reductions import DemographicParity # EqualizedOdds, TruePositiveRateParity
from fairlearn.reductions import GridSearch
sweep = GridSearch(
                   estimator=LGBMClassifier(random_state=1), 
                   constraints=DemographicParity(),  # EqualizedOdds, TruePositiveRateParityなども使用可
                   grid_limit = 1, # lambdaを-grid_limitからgrid_limitまでに設定
                   grid_size = 50  # 何分割して実行するか

)
sweep.fit(X_train, y_train, sensitive_features=s_train)

# predict
y_train_pred_indp = sweep.predict(X_train) # 最もDemographicParityが小さいものが選ばれる
y_test_pred_indp = sweep.predict(X_test)

# print result
print("Train")
print_result_summary(y_true=y_train, y_pred=y_train_pred_indp, s=s_train)
print("\nTest")
print_result_summary(y_true=y_test, y_pred=y_test_pred_indp, s=s_test)

結果として、

# 結果 （Reduction with DemographicParity）
Train
Demographic Parity:		[mean] 0.544,	[s=1] 0.544,	[s=0] 0.536,	[abs_diff] 0.008
TruePositiveRate Parity:	[mean] 0.669,	[s=1] 0.660,	[s=0] 0.788,	[abs_diff] 0.128
FalsePositiveRate Parity:	[mean] 0.414,	[s=1] 0.420,	[s=0] 0.378,	[abs_diff] 0.042

Test
Demographic Parity:		[mean] 0.524,	[s=1] 0.524,	[s=0] 0.520,	[abs_diff] 0.004
TruePositiveRate Parity:	[mean] 0.591,	[s=1] 0.589,	[s=0] 0.616,	[abs_diff] 0.027
FalsePositiveRate Parity:	[mean] 0.456,	[s=1] 0.454,	[s=0] 0.470,	[abs_diff] 0.016

となり、指定したDemographic Parityが高い精度で保たれていることがわかります。

ちなみに上のコードのpredictで用いられるモデルはGridSearchで探索した範囲で、指定したMetricsを最小にするモデルになります。一応確かめるために各predictorを取得して、demographic_parity_differenceを計算して、最終的にそれをソートしてみます。

sweep_preds = [predictor.predict(X_train) for predictor in sweep.predictors_] # 各predictorの予測値を取得
dp_diff_list = [
    demographic_parity_difference(y_train, preds, sensitive_features=s_train)
    for preds in sweep_preds
]
print(np.sort(dp_diff_list)[:5]) #top-5のdemographic_parity_diffを確認, sweep.predictによる結果と同一であることを確認できる

> array([0.00844444, 0.01511111, 0.038     , 0.05511111, 0.05911111])

となり、最小のモデルがTrainのDemographic Parityのabs_diffと一致しているため、指定したmetricsでの最良のモデルが選択されていることがわかります。

各アルゴリズムの実行 (ThresholdOptimizer)

最後にThresholdOptimizerを使用してみます。こちらも使い方はシンプルで、constraintsの引数に加味したい公平性を"demographic_parity"や"equalized_odds"などを指定します。

他の引数としては、すでに学習しているモデルをestimatorとして使用する場合にはprefit=Trueとして、学習済みでない場合はインスタンンスをestimatorで設定しprefit=Falseに設定します。 そしてfit, predictを行えば実行できます。

1点注意すべき点として、ThresholdOptimizerに関しては、predictが乱数依存になっているので、再現性を確保するためにはrandom_stateを指定する必要があります。

from fairlearn.postprocessing import ThresholdOptimizer
# set 
optimizer = ThresholdOptimizer(
    estimator=baseline_model, 
    constraints="demographic_parity", # 他に’{false,true}_{positive,negative}_rate_parity’’equalized_odds’ が使用可能
    prefit=True # 学習済みモデルを渡している場合 True, prefit=Falseでestimator=LGBMClassifier(random_state=1)でも同じ結果
)

# fit optimizer
optimizer.fit(X=X_train,y=y_train.reshape(-1,1), sensitive_features=s_train)
y_test_pred_postdp = optimizer.predict(X_test, sensitive_features=s_test, random_state=20) #乱数固定
y_train_pred_postdp = optimizer.predict(X_train, sensitive_features=s_train, random_state=100) # 乱数固定

# print result
print("Train")
print_result_summary(y_true=y_train, y_pred=y_train_pred_postdp, s=s_train)
print("\nTest")
print_result_summary(y_true=y_test, y_pred=y_test_pred_postdp, s=s_test)

結果は、以下のようになり、特に指定したDemographic Parityについて、GridSearchほどでもないですが、保たれていることがわかります。これはpredictの乱数によって変わるのであまりどっちがよいなどというものではなさそうですね。

# 結果 （ThresholdOptimizer with DemographicParity）
Train
Demographic Parity:		[mean] 0.578,	[s=1] 0.575,	[s=0] 0.600,	[abs_diff] 0.025
TruePositiveRate Parity:	[mean] 0.693,	[s=1] 0.689,	[s=0] 0.741,	[abs_diff] 0.051
FalsePositiveRate Parity:	[mean] 0.459,	[s=1] 0.451,	[s=0] 0.511,	[abs_diff] 0.060

Test
Demographic Parity:		[mean] 0.557,	[s=1] 0.558,	[s=0] 0.554,	[abs_diff] 0.004
TruePositiveRate Parity:	[mean] 0.625,	[s=1] 0.624,	[s=0] 0.640,	[abs_diff] 0.015
FalsePositiveRate Parity:	[mean] 0.489,	[s=1] 0.486,	[s=0] 0.509,	[abs_diff] 0.023

以上3つのメソッドをみてきましたが、とても簡単に使えるようになっていますね。記事中のコードや、記事内では紹介できなかったEqualizedOddsを指定したバージョンなどのコードはこちらに置いています。

まとめ

fairlearnはまだまだ発展性がありそうなパッケージですし、現状でもある程度の問題には対処できそうです。公平性は社会的ニーズが強そうな話題なので今後も研究が盛んにされていきそうです。

個人的には、センシティブなデータの影響を除くためにセンシティブなデータを使うっていうのは、確かにそうなんだろうけど、センシティブなデータがなければ学習できない、かつ手法によっては引き続き収集しなきゃいけないってことがちょっと微妙っぽいなと思いました。その中でも、predict時にセンシティブなデータを使わなくて済むReduction Approachに関しては比較的使いやすいなと思いました。ただモデル更新するときどうするのかはどうしようって感じですが。。。

また、今まで人間にもできていないケースがあったであろう公平性を機械学習モデルなどに求めるために、人間の理想とする公平性が定量化されて機械学習モデルに使えるようになっていくといったプロセスは非常に興味深いです。

公平性(fairness)に関して、この記事では集団公平性というものにしか触れられていませんが、個人公平性など様々な研究がなされているようです。日本語ですばらしいまとめ*4*5があるので興味がある方はぜひ。 

*1:レビュー論文1

[1908.09635] A Survey on Bias and Fairness in Machine Learning

*2:レビュー論文2

[2010.04053] Fairness in Machine Learning: A Survey

*3: fairlearnのThresholdOptimizerの元論文

[1610.02413] Equality of Opportunity in Supervised Learning

*4:公平性に配慮した学習と理論的課題 (2018) 機械学習の公平性の定義も踏まえ技術的な話題に詳しいです。

https://ibisml.org/ibis2018/files/2018/11/fukuchi.pdf

*5:機械学習と公平性(2020), 公平性の背景などに詳しいです。

http://ai-elsi.org/wp-content/uploads/2020/01/20200109-fairness_sympo.pdf