🦉

Pythonで部分識別①

2022/01/29に公開

Python

統計

tech

はじめに

計量経済を勉強していた後輩から以下の本を勧められました。

部分識別という、計量経済学で急速に発展している研究手法がまとめられた本です。
発想が面白かったので簡単にご紹介しようと思います。技術系記事感を出すために、Pythonでちょこっとコードを書いてみました。

部分識別とは?

識別

計量経済学で使われる言葉ですが、そもそも「識別」とは何でしょうか。
「識別できる」とは、「データが無限個ある(=母集団全てが標本)時に、求めたいパラメータの真の値が決定できる」ということらしいです。識別はできる or できないで表されます。
一方で、部分識別は求めたいパラメータをバウンド(幅)で識別する方法だそうです。あまり自分もこの辺は腹落ちできておらず上手く説明できないので、上述の本を読んでいただければと思います。

データセットを用いた説明

話を分かりやすくするために、具体例を出して話を進めます。
今回は効果検証入門でも取り扱われていた MineThatData E-Mail Analytics And Data Mining Challenge datasetというデータを利用します。ECサイトのユーザーに対してRCTでメール送信したデータです。

今回はこのデータセットに対して部分識別の考え方を適用することで、「メール送付による売上への効果」のバウンドを算出します。

なお、ここでは以下のように記号を使用します。

$z$ : 処置変数。メールの送付有無であり、0と1の二値変数である。
$y(0)$ : メール送付されなかった場合の購買金額
$y(1)$ : メール送付された場合の購買金額

全てのユーザに対して $y(0)$ と $y(1)$ がありますが、実際に観測されるのは処置に応じてどちらか一方しか観測されません。

また、簡便のために女性向けメールが送信されたユーザーのデータを除去し、男性向けメールが送信されたユーザーは $z=1$ 、メールを受け取っていないユーザーは $z=0$ とします。

import pandas as pd


df = pd.read_csv('http://www.minethatdata.com/Kevin_Hillstrom_MineThatData_E-MailAnalytics_DataMiningChallenge_2008.03.20.csv')
# 女性向けメールが送信されたユーザーを除去
df = df.query('segment != "Womens E-Mail"')
# 処置変数を作成
df['treatment'] = df['segment'].apply(lambda x: 1 if x == 'Mens E-Mail' else 0)

以下数式がしばらく続きます。

算出したい効果

今回知りたいのは平均処置効果(ATE)です。具体的には「全てのユーザーが、メールを送付されることで平均的に購買金額がどれくらい増加するのか」ということです。これは以下の式で表せます。

E[y(1) - y(0)] = E[y(1)] - E[y(0)]

$E[y(1)]$ の識別

ここで $E[y(1)]$ の識別を考えます。 $E[y(1)]$ は以下のように分解できます。

E[y(1)] = E[y(1)|z=1]P(z=1) + E[y(1)|z=0]P(z=0)

各項を日本語で説明すると

$E[y(1)|z=1]$ …メール送付されたユーザーが、メール送付された場合の購買金額の期待値
$E[y(1)|z=0]$ …メール送付されていないユーザーが、メール送付された場合のユーザーの購買金額の期待値

となります。

上の説明で明らかですが、前者は観測することができますが、後者は観測することができません。これを計量経済学では「 $E[y(1)|z=0]$ は識別できない」と言います。

$E[y(0)]$ の識別

$E[y(1)]$ と同様に式変形をします。

E[y(0)] = E[y(0)|z=1]P(z=1) + E[y(0)|z=0]P(z=0)

各項を説明すると

$E[y(0)|z=1]$ …メール送付されたユーザーが、メール送付されていない場合の購買金額の期待値
$E[y(0)|z=0]$ …メール送付されていないユーザーが、メール送付されていない場合のユーザーの購買金額の期待値

となります。

こちらでは $E[y(0)|z=1]$ を観測することができません。そのため、 $E[y(0)|z=1]$ は識別できません。

伝統的な計量経済学での効果算出

伝統的な計量経済学では、 $y$ と $z$ が独立になるような条件が成り立つようにすることで、因果効果を推定することが多いと書かれています。
今回のデータセットはRCTを行っているデータのため、 $y$ と $z$ が独立になります。そのため以下の式が成り立ちます。

E[y(1)] = E[y(1)|z=1] = E[y(1)|z=0] \\ E[y(0)] = E[y(0)|z=1] = E[y(0)|z=0]

つまり観測できる値を使って

E[y(1)] - E[y(0)] = E[y(1)|z=1] - E[y(0)|z=0]

以上のようにメール送付の効果を算出することができます。

Pythonにて効果を算出してみます。

agged = df.groupby('treatment')['spend'].mean().reset_index()
agged.query('treatment == 1')['spend'].values[0] - agged.query('treatment == 0')['spend'].values[0]
### => return
0.7698271558945375

メール送付による売上増加は、0.77ドルということが分かります。t検定をすると有意であることが確認されます。

RCTを行っていないデータでは、傾向スコアを用いたマッチングや操作変数法などを用いて、観測できない値を観測できる値から推測することで、求めたい値を識別しようとします。

部分識別

これに対し、部分識別は全く違うアプローチで識別問題に立ち向かいます。
何も仮定しなければ、識別できない値は $y(\cdot)$ の値域に入るだろうと考えます。

ここで $y(\cdot)$ の値域を $[\underline{y}, \overline{y}]$ とします。すると、識別できなかった $E[y(1)|z=0]$ や $E[y(0)|z=1]$ の値域も $[\underline{y}, \overline{y}]$ となります。

これによって、識別できなかった $E[y(1)]$ のバウンドは以下のような不等式で表されます。

E[y(1)|z=1]P(z=1) + \color{red}\underline{y}\color{black}P(z=0) \leq E[y(1)] \leq E[y(1)|z=1]P(z=1) + \color{red}\overline{y}\color{black}P(z=0)

識別できなかった $E[y(1)|z=0]$ が $\color{red}\underline{y}, \overline{y}$ で置換されていることが分かります。

同様に、識別できなかった $E[y(0)]$ のバウンドは以下の不等式で表せます。

\color{red}\underline{y}\color{black}P(z=1) + E[y(0)|z=0]P(z=0) \leq E[y(0)] \leq \color{red}\overline{y}\color{black}P(z=1) + E[y(0)|z=0]P(z=0)

識別できなかった $E[y(0)|z=1]$ が $\color{red}\underline{y}, \overline{y}$ で置換されていることが分かります。

ATEのバウンド算出

以上によって $E[y(1)]$ と $E[y(0)]$ のバウンドが求まりました。
今回知りたい、メール送信による売上増加効果は $E[y(1)]-E[y(0)]$ でした。このバウンドは

E[y(1)]の下限 - E[y(0)]の上限 \leq E[y(1)] - E[y(0)] \leq E[y(1)]の上限 - E[y(0)]の下限

で算出することができます。これが何も仮定しない場合のバウンドとなります。

Pythonにて実装

これまで説明してきた、何も仮定を置かない場合の部分識別について、バウンドの算出をPythonにて実装しました。

def calc_bound(df, y, z):
    # yの値域を算出
    y_min = df[y].min()
    y_max = df[y].max()

    # P(z=1), P(z=0), E[y(1)|z=1], E[y(0)|z=0] を算出する
    agged = df.groupby(z)[y].agg(['mean', 'count']).reset_index()
    agged['prob'] = agged['count'] / agged['count'].sum()

    # z=0, 1でデータを分けておく
    a_z0 = agged[agged[z] == 0]
    a_z1 = agged[agged[z] == 1]

    # E[y(1)]のバウンドを算出
    e_y1_bound_min = a_z1['mean'].values[0] * a_z1['prob'].values[0] + y_min * a_z0['prob'].values[0]
    e_y1_bound_max = a_z1['mean'].values[0] * a_z1['prob'].values[0] + y_max * a_z0['prob'].values[0]

    # E[y(0)]のバウンドを算出
    e_y0_bound_min = a_z0['mean'].values[0] * a_z0['prob'].values[0] + y_min * a_z1['prob'].values[0]
    e_y0_bound_max = a_z0['mean'].values[0] * a_z0['prob'].values[0] + y_max * a_z1['prob'].values[0]

    # E[y(1)] - E[y(0)]のバウンドを算出
    bound_min = e_y1_bound_min - e_y0_bound_max
    bound_max = e_y1_bound_max - e_y0_bound_min

    return bound_min, bound_max
    
    
# メール送信による売上増加効果を算出
calc_bound(df, y='spend', z='treatment')
### => return
(-249.1209170910286, 249.87908290897144)

メール送付による売上増加効果は -249 ~ 249ドル という結果になります。何も仮定していない分、バウンドが広すぎるのであまり嬉しくない結果になってしまいました。

実装していて気になった点としては、 $y$ の上限下限をどのように設定するのが良いのかなぁなんて思いました。今回は最小値と最大値を利用しましたが、購買金額の分布を見るとデータのほとんどがゼロで最大値は499という分散の大きいデータでした。明らかに売上増加効果が499ドルとは考えにくいです。だからといって変に上限を小さくするのは、それはそれで恣意性が入るので、このあたりどうするのかなぁと気になりました。

おわりに

今回は部分識別の簡単な説明と、メールマーケティングデータに対してメール送付による売上増加効果のバウンドを求めるためのPythonコードを書きました。
部分識別の発想自体は非常にシンプルなものの、これまで勉強してきた因果推論の考え方とは全く違う考え方で、個人的には目から鱗でした。とても面白かったです。
今回の結果を見ても分かる通り、実用性の観点で言うと、何も仮定せずに算出するとバウンドの広さから分かる通り正直使い物にならないでしょう。

これに対して『部分識別入門』の2章以降は、皆が認めやすいような弱い仮定を課した場合に「どのようにバウンドを算出するのか」、「追加する仮定によってバウンドがどれくらい狭まるのか？」が書かれています。元気があればこちらもPythonで実装できたらなと思います。

以上です。