はじめに
傾向スコアというものがあると聞いたので、マーケティングとかのバイアス補正に使えるんじゃないかと思って調べてみました。
まずは、理論的な面から説明してみて、その後 python 使って例題やってみます。
理論的なこと
忍び寄るバイアス
まずは、簡単な例を考えてみましょう。
ある広告を打って、ユーザーの購買行動へ効果があったか調査するとします。広告を打たれたユーザは $N_1$ 人でその購買額は $y_1$、広告を打たれなかったユーザは $N_0$ 人で購買額は $y_0$ だったとします。 広告を打たれたユーザの平均購買額は
\frac{1}{N_1} \sum_{i=1}^{N_1}y_{1i}
広告を打たれなかったユーザの平均購買額は
\frac{1}{N_0} \sum_{i=1}^{N_0}y_{0i}
となるので広告の効果は
\frac{1}{N_1} \sum_{i=1}^{N_1}y_{1i} - \frac{1}{N_0} \sum_{i=1}^{N_0}y_{0i} \tag{1.1}
となりそうな気がします。。。 本当でしょうか？
もちろんこれで十分な場合もあります。しかし、(1.1) では、広告以外の影響が紛れ込んでいることが多々あります。
例えば、普段からよく買い物をしてくれる優良ユーザにだけ広告が打たれていたとしたらどうでしょう。(1.1) 式の1項目は優良ユーザの購買額に大きく引きずられます。一方で、2項目は普段からあまり買い物をしてくれないユーザが集まり値が小さくなります。引き算をしているので、これも
(1.1) の値を底上げする方向に働きます。以上の結果として、広告効果を過大評価してしまうことになるでしょう。
このようにデータには本来は見るつもりのなかった偏りが忍び込んでくることが多々あります。この偏りをバイアスと呼び、広告効果を正しく評価するためには、バイアスをいかにして除去するかがカギとなります。
並行世界からの救済「反実仮想」
さて問題は、「バイアスを除去するにはどうしたらいいか」ですが、Rubin(1974) などによって導入された「反実仮想」というフレームワークで捉えなおしてみると見通しが良くなります。
では、(1.1) に戻って考えてみましょう。先に説明した優良ユーザバイアスの本質は何でしょうか。それは、(1.1) の比較が、優良ユーザーvs非優良ユーザの差分になってしまっていることです。そうであるならば、バイアスを除去するためには傾向の似たユーザ同士の比較に持ち込めればよさそうです。もっと言えば、同一のユーザに対して広告に接触した状態の購買額と、接触していない状態の購買額の差分を評価することができれば純粋な広告効果そのものを取り出せそうです。。。
広告に接触した世界と、接触していない世界・・・現実にはどちらか一方のみ実現するので、もう一方は現実に反した想定です。・・・はい、反実仮想いただきました。
平行世界に迷い込んでいる感はありますが、もうちょっと考えてみましょう。先の例では、広告に接触したユーザの購買額を $y_1$、接触していないユーザの購買額を $y_0$ としました。ここまで来たのなら思い切って、各ユーザが両方の値を持っているとしてしまいましょう。つまり、各ユーザについて広告に接触した場合の購買額が $y_1$ 、接触しなかった場合の購買額が $y_0$ であるとします。そして、全ユーザ数を $N$ として
\frac{1}{N} \sum_{i=1}^{N}(y_{1i} - y_{0i}) \tag{2.1}
を考えてみてください。(1.2) の第1項、2項ともに全ユーザを含んでいるため以前のような優良ユーザvs非優良ユーザによるバイアスは存在しません。
もっとも、このままでは単一世界に住む私たちの感覚に合わないので、各ユーザが広告接触者として選択された場合には $y_1$ のみが、非接触者として選択された場合は $y_0$ のみが観測されるのだ、と考えます。同時に、広告接触/非接触者のどちらとして選択されたかを表す変数 $z$ を導入しておくと便利です。便宜上、$z=1$ が接触、$z=0$ が非接触に対応するものとします。すると、実際に計測することのできる購買額 $y$ は $y_1$, $y_0$, $z$ の関数として
y(z, y_1, y_0) = z y_1 + (1-z) y_0 \tag{2.2}
と表せます。
さて、(2.1) を評価できることが理想だとしても、計測不可能な値を含んでいるため直接求めることはできません。しかし、直接は求められなくとも、何らかの「妥当な」仮定のもと近似的にでも求められれば嬉しいはずです。
数式で整理してみよう
では、どうすれば(2.1)を近似的にでも求められるか。先ほど導入した記号を使って整理してみましょう。
ぼくらのみかた「大数の法則」
その前に、以下の式変形で核となる大数の法則について触れておきます。
私たちは普段、平均として「サンプル平均」を用います。つまり、ある確率変数 $x$ があったとして $f(x)$ の平均を
\frac{1}{N}\sum_{i=1}^{N}f(x_i) \tag{3.1}
として計算します。しかし、上式の値 $x_1,x_2, ... ,x_N$ は $x$ の従う確率分布 $p(x)$ からのサンプリング結果として得られたものにすぎず、真の平均は
E(f(x)) = \sum_x p(x) f(x) \tag{3.2}
となるはずです。この2つの関係は何でしょうか。そう、大数の法則です。すなわち、サンプル数無限大の極限で
\frac{1}{N}\sum_{i=1}^{N}f(x_i) \to \sum_x p(x) f(x) = E(f(x)) \quad (N \to \infty) \tag{3.3}
となるのです。よって、真の平均値 (3.2) の近似値としての期待をこめて、(3.1) を用いているわけです。ちなみに、(3.3) の関係が成り立つことを、(3.1) は (3.2) の普遍推定量であると言います。
(3.1) のようにサンプリングされた世界の代わりに、(3.2) のような確率分布の世界で話をすると理論的な話がすっきりします。
普遍推定量を探せ
さて、ここまでで出てきた数式でサンプル数無限大の極限を取ってみましょう。先に導入した記号(2.2) を使うと便利です。
まず、(1.1) の1項目に注目します
\begin{align}
&\frac{1}{N_1} \sum_{i=1}^{N_1}y_{1i} = \frac{\sum_{i=1}^{N}z_iy_i}{\sum_{i=1}^{N}z_i} = \frac{\frac{1}{N}\sum_{i=1}^{N}z_iy_i}{\frac{1}{N}\sum_{i=1}^{N}z_i} \\
\to & \frac{E(zy)}{E(z)} = \frac{\sum_{y_1}y_1p(y_1, z=1)}{p(z=1)}
=E(y_1 | z=1) \tag{3.4}
\end{align}
同様に(1.1)の2項目は
\begin{align}
&\frac{1}{N_0} \sum_{i=1}^{N_0}y_{0i} = \frac{\frac{1}{N}\sum_{i=1}^{N}(1-z_i)y_i}{\frac{1}{N}\sum_{i=1}^{N}(1-z_i)} \\
\to & \frac{E((1-z)y)}{E(1-z)} = \frac{\sum_{y_0}y_0p(y_0, z=0)}{p(z=0)}
=E(y_0 | z=0) \tag{3.5}
\end{align}
となることがわかるので、(1.1) は
E(y_1 | z=1) - E(y_0 | z=0) \tag{3.6}
の普遍推定量だったのだと分かります。広告接触者（$z=1$）の購買額 $y_1$ ひくことの非接触者（$z=0$）の購買額 $y_0$ なので直感的にも納得だと思います。
次に、(2.1) をみてみましょう。こちらは簡単で、
\begin{align}
\frac{1}{N} \sum_{i=1}^{N}y_{1i} \to E(y_1), \quad
\frac{1}{N} \sum_{i=1}^{N}y_{0i} \to E(y_0) \tag{3.7}
\end{align}
なので、
E(y_1) - E(y_0) \tag{3.8}
の普遍推定量なのでした。
一般に、ある真の平均値一つに対して、普遍推定量は一意に決まるわけではありません。ということは、バイアスを除いた広告効果を評価するためには、実際に計測可能な値のみをで構成された (3.8) の普遍推定量がわかればよいわけです。少しずつ道が開けてきましたね。
IPW 推定量
(3.8) の普遍推定量を探すにあたって、まずは、なぜバイアスなどというものが生じるのか考えてみましょう。
共変量と、強く無視できる割り当て
もし、$z=1$ と $z=0$ の割り当て方が何にも依存しないならば、すなわち完全にランダムであるならば、バイアスなど存在しようがないはずです。逆に言えば、バイアスが存在するからには、それを特徴づける何かが背後に隠れているはずです。これを共変量 $x$ と名付けます。
$x$ がバイアスを特徴づけることをどう数式で表現したらよいでしょうか。バイアスが存在しないならば、(3.6) と (3.8) は一致するはずです。$E(y_1|z=1)$ と $E(y_1)$、$E(y_0|z=0)$ と $E(y_0)$ がそれぞれ一致すると期待してもよいでしょう。そしてこれらは、$p(y_1, y_0, z) = p(y_1, y_0)p(z)$ すなわち、$z$ と {$y_1,y_0$} が確率変数として独立
\left\{ y_1, y_0 \right\} \perp z
であれば成立します。では、共変量 $x$ が存在している場合はどうでしょうか。$z=0,1$ の割り当てが $x$ に支配されているならば、$x$ を固定してしまえばバイアスが生じる自由度は消えるはずです。つまり、$x$ を条件づけた時に $z$ と {$y_1,y_0$} が確率変数として独立
\left\{ y_1, y_0 \right\} \perp z | x \tag{4.1}
であることがカギになっていそうです。(4.1) は「強く無視できる割り当て」と呼ばれ、この条件を満たす共変量 $x$ を用いることで (3.8) の普遍推定量を求めることが可能となります。
確かに、共変量 $x$ が必ずしも手元にあるとは限りません。ですので、前に述べた「妥当な条件」のひとつが、「共変量が利用可能なこと」となります。利用できない場合は、残念ですが他の手法を探るしかありません。
傾向スコア
さて、(3.8) の普遍推定量を求めていくわけですが、何か指針がほしいところです。そこで、計測可能な値のみで構成する意味でも (1.1) を出発点にしてみましょう。(1.1) の問題点を思い出してみると、核心部分は「1項目に含まれる優良ユーザと、2項目に含まれる非優良ユーザが期待値を押し上げる」ことでした。ということは、これらの優良/非優良ユーザの weight を下げてしまえばよいのではないでしょうか。でも、どうやって重みを計算しましょう。
幸い、私たちの手もとには共変量 $x$ があります。$z=1$ は、ほぼ優良ユーザですが、彼らは $x$ によって特徴づけられています。よって、 $x$ を説明変数とした優良ユーザーのスコアリングを行い、weight の調整を行うのはどうでしょう。
これらを踏まえ、傾向スコアと呼ばれる以下の関数を導入します。
e(x) := p(z=1|x) = E_{z}(z|x) \tag{4.2}
傾向スコア(4.2)は共変量 $x$ を 与えた時の $z=1$ に割り当てられる確率となります。もちろん、傾向スコアの真の値などは知りようがありませんが、通常はロジスティック回帰などを使って見積もることが多いかと思います。
IPW 推定量
優良ユーザーの weight を落とすため (1.1) の1項目を
\hat{E}_1 := \frac{\sum_{i=1}^N\frac{z_i}{e_i}{y_i}}{\sum_{i=1}^N\frac{z_i}{e_i}} \tag{4.3}
と修正してみましょう。同様に、非優良ユーザーの weight を落とすため (1.1) の2項目を
\hat{E}_0 := \frac{\sum_{i=1}^N\frac{1-z_i}{1-e_i}{y_i}}{\sum_{i=1}^N\frac{1-z_i}{1-e_i}} \tag{4.4}
としてみます。全て計測可能な値で書かれていることに注意してください。実は、これらが欲しかった普遍推定量になっています。いわゆる、IPW推定量と呼ばれるものです。ですので、普遍推定量であることを確かめれば目的達成です。
まず、$e(x)$ の定義 (4.2) に注意して、以下の関係式を導きます。
E\Big(\frac{z}{e}\Big) = \sum_x p(x)\frac{1}{e(x)}\sum_z zp(z|x)
= \sum_x p(x)\frac{1}{e(x)} p(z=1|x) = 1 \tag{4.5}
また、$z=0,1$ であることより、$z^2=z$ となるので
zy = z(zy_1 + (1 - z)y_0) = xy_1
であることがわかります。これらを用いて
\begin{align}
E\Big(\frac{zy}{e}\Big) = E\Big(\frac{zy_1}{e}\Big)
&= E_x\Big(E_{zy_1}\Big(\frac{zy_1}{e}\Big|x\Big)\Big) \\
&= E_x\bigg(\frac{1}{e(x)}E_{zy_1}\Big(zy_1\Big|x\Big)\bigg) \\
&= E_x\bigg(\frac{1}{e(x)}E_z(z|x)E_{y_1}(y_1|x)\bigg) \\
&= E_x\Big(E_{y_1}(y_1|x)\Big) = E(y_1) \tag{4.6}
\end{align}
と計算しておきます。1行目から2行目の変形では $x$ が固定されていることを、2行目から3行目の変形では条件 (4.1) を用いています。これで道具がそろいました。
さっそく、(4.3) で $N \to \infty$ の極限をとってみると
\hat{E}_1 = \frac{\frac{1}{N}\sum_{i=1}^{i=N} \frac{z_i}{e_i}y_i}{\frac{1}{N}\sum_{i=1}^{i=N} \frac{z_i}{e_i}} \to \frac{E\big(\frac{zy}{e}\big)}{E\big(\frac{z}{e}\big)} = E\Big(\frac{zy}{e}\Big) = E(y_1) \tag{4.7}
が得られます。$\hat{E}_0$ でも同様にして $\hat{E}_0(y) \to E(y_0)$ が示されます。
以上により、
\hat{E}_0 - \hat{E}_1 \tag{4.8}
が (3.8) の普遍推定量であることが示されました。
長い道のり、本当にお疲れ様でした。
ATE と愉快な仲間たち
ちなみに、(3.8) はATEと呼ばれていまして、(4.8) が普遍推定量
\text{ATE} := E(y_1) - E(y_0) \sim \frac{\sum_{i=1}^N\frac{z_i}{e_i}{y_i}}{\sum_{i=1}^N\frac{z_i}{e_i}} - \frac{\sum_{i=1}^N\frac{1-z_i}{1-e_i}{y_i}}{\sum_{i=1}^N\frac{1-z_i}{1-e_i}} \tag{4.9}
なのはここまで見てきた通りなのですが、他にも
\begin{align}
\text{ATT} := E(y_1|z=1) - E(y_0|z=1) \tag{4.10} \\
\text{ATU} := E(y_1|z=0) - E(y_0|z=0) \tag{4.11}
\end{align}
というのがあって、ATT は広告を打ったユーザーに対しての効果。ATU は今回広告を打たなかったユーザに打ってみた場合に期待できる効果という意味になります。普遍推定量は (3.4) と (3.5) に加えて、ベイズ定理から
\begin{align}
E(y_1|z=0) = \frac{E(y_1) - p(z=1)E(y_1|z=1)}{p(z=0)} \tag{4.12}\\
E(y_0|z=1) = \frac{E(y_0) - p(z=0)E(y_0|z=0)}{p(z=1)} \tag{4.13}\\
\end{align}
と求めたものを使えばよいかと思います。
以上で、理論編終了です。
pythonで例題
数式だけではイメージがつきにくいかと思いますので、簡単な例題を python など使いつつ試してみましょう。自作のデータを使った方が、何が起きているのか把握しやすい気がするので、ダミーデータの生成から始めましょう。
データ作成
たっぷりとバイアスの乗ったデータを作ります。パッケージは定番の pandas と numpy があれば十分です。
import pandas as pd
import numpy as np
理論編での「広告打って購買額への影響調べました」の例を念頭に置いてデータを作ります。作るデータはユーザリストです。ユーザーが広告に接触したかどうかのフラグ、そして、その結果としての購買額か何かの連続値が入っているものを生成します。それから、共変量としてのデモグラ情報もあれば例題として十分でしょう。デモグラは、以下のように年齢・性別とします。
もちろん、フラグは広告接触と解釈する必然性はありませんし、連続値も購買額である必要はありません。なんらかの「値」です。
age = ["10代", "20代", "30代", "40代", "50代", "60代", "70代"]
gender = ["女性", "男性"]
以下、データ生成用の関数です。上から順に、デモグラ付きユーザリストを作成、（広告接触）フラグの付与、購買額か何かの連続値の付与と続きます。フラグは 'x_flag'、連続値は 'x_rate' と言う名前で 'X_rate' の値は 0 ~ 1 の間に入るように規格化しています。
ユーザ生成時にデモグラ情報は一様分布に従って振ります。そのため、年代・性別間でのユーザー数の偏りはありません。一方で、'x_flag' と 'x_rate' にはバイアスをのせます。関数の targets 変数に設定したデモグラ値に対して、バイアスがか乗るようになっています。
'X_flag' については、フラグ '1' の付与確率は通常ユーザーでは 10% なのに対して、targets で設定したデモグラ値を持つユーザーには 90% となります。また、'x_rate' の値については、フラグ '0' のユーザーの場合は平均 0.2 標準偏差 0.1 の正規分布に従いますが、フラグ '1' のユーザーでは 0.2 ポイント分だけかさ上げされます。この 0.2 ポイントが広告効果です。さらに、targets 変数で設定したデモグラ値を持つユーザに対しては、フラグの値によらず、平均 0.8
標準偏差 0.2 の正規分布に従う値をセットします。超優良ユーザーに広告なんて関係ないというイメージです。1 以上 0 以下の値は丸めます。
def generate_users(n_users, **features):
"""make user list where columns are 'user_id' and feature columns ('age', 'gender',... etc)
values of the feature columns are randomly chosen from the given list
"""
user_id = [i + 1 for i in range(n_users)]
df = pd.DataFrame({"user_id":user_id})
np.random.seed(100)
for col in features:
feature = features[col]
values = np.random.choice(feature, n_users)
df[col] = values
return df
def add_flag(df, **targets):
"""add flag column named as 'x_flag'
bias is added to the users specified in the variables '**targets'
"""
n_users = df.shape[0]
# specify users to whom bias is added
is_target = [True for _ in range(n_users)]
for col in targets:
target_values = targets[col]
is_value = [False for _ in range(n_users)]
for value in target_values:
tmp_bool = (df[col] == value)
is_value = np.logical_or(is_value, tmp_bool)
is_target = np.logical_and(is_target, is_value)
#assign flag
flag = [0, 1]
#for ordinary users, 10% of them are assigned to 'flag=1'
weight = [0.9, 0.1]
np.random.seed(200)
df["x_flag"] = np.random.choice(flag, n_users, p=weight)
#for users specified above, 90% of them are assigned to 'flag=1'
weight = [0.1, 0.9]
np.random.seed(300)
df.ix[is_target, "x_flag"] = np.random.choice(flag, n_users, p=weight)
return df
def add_rate(df, **targets):
"""add real valued column named as 'x_rate', where the values are normalized between 0 and 1
for the users with 'flag=1', values in 'x_rate' are increased
for the users specified by '**targets', values in 'x_rate' are much higher, regardless of 'flag=1 or 0'
"""
n_users = df.shape[0]
#specify users to whom high 'x_rate' values are set
is_target = [True for _ in range(n_users)]
for col in targets:
target_values = targets[col]
is_value = [False for _ in range(n_users)]
for value in target_values:
tmp_bool = (df[col] == value)
is_value = np.logical_or(is_value, tmp_bool)
is_target = np.logical_and(is_target, is_value)
#for users with 'flag=0', 'x_rate' values are sampled from Normal(0.2, 0.1)
np.random.seed(200)
df["x_rate"] = np.random.normal(0.2, 0.1, n_users)
#for users with 'flag=1', 'x_rate' values are increased by 0.2 point
is_flag = (df["x_flag"] == 1)
df.ix[is_flag, "x_rate"] += 0.2
#for users specified by '**targets', 'x_rate' values are sampled from Normal(0.8, 0.2)
np.random.seed(300)
df.ix[is_target, "x_rate"] = np.random.normal(0.8, 0.2, n_users)
#for convenience normalize to [0 ,1]
df.ix[df["x_rate"] > 1, "x_rate"] = 1
df.ix[df["x_rate"] < 0, "x_rate"] = 0
return df
それでは、実際にユーザを生成してみます。ユーザ数は 10万人とします。バイアスは、20~30代の男性に乗せます。つまり、この層は広告接触率が異常に高いにも関わらず、すでに超優良ユーザのため広告が効かないという設定です。
n_users = 100000
df = generate_users(n_users, age=age, gender=gender)
df = add_flag(df, age=["20代", "30代"], gender=["男性"])
df = add_rate(df, age=["20代", "30代"], gender=["男性"])
作ったデータを確認してみよう
さて、生成したデータを確認してみましょう。
df.head(10)
きちんと出来てそうなので集計してみましょうか。まずは、デモグラ×フラグ別ユーザ数です。
df_count = df.groupby(["x_flag", "gender", "age"])["user_id"].count().unstack(level=[1,2]).sort_index(ascending=False)
df_count.ix["ratio of 1 (%)"] = 100 * df_count.ix[1, :] / (df_count.ix[1, :] + df_count.ix[0, :])
df_count.astype(int)
ratio of 1 はフラグ '1' ユーザの割合です。確かに20~30代男性でフラグ '1' のユーザ割合が異様に高くなっていることが見て取れます。
次に、'x_rate' の平均値を比較してみましょう。
df_rate = df.groupby(["x_flag", "gender", "age"])["x_rate"].mean().unstack(level=[1,2]).sort_index(ascending=False)
df_rate.ix['diff (1-0)', :] = df_rate.ix[1, :] - df_rate.ix[0, :]
df_rate.round(3)
diff の行はフラグ '1' と '0' で平均値の差分を取ったものです。基本的にフラグ '1' が立っている方が 0.2 ポイントほど値が高くなっています。ただし、20~30代男性では、フラグの値によらず大きな値となっています。また、フラグ '1' と '0' の間に差は見受けられません。
この表を見ても、広告効果は 0.2 ポイントと考えるのが妥当かなという感じがします。
素朴に広告効果を評価してみると
では、理論編で見たように広告効果を評価してみましょう。まずは素朴に (1.1) で計算してみます。
E1 = np.mean(df.ix[df["x_flag"] == 1, "x_rate"])
E0 = np.mean(df.ix[df["x_flag"] == 0, "x_rate"])
E = E1 - E0
print("Expected Effect = {0} (x_flag=1 : {1}, x_flag=0 : {2}) ".format(round(E, 3), round(E1, 3), round(E0, 4)))
0.418 ポイントです。予想通り、過大評価されているように思えます。明らかに原因は20~30代男性の強力なバイアスです。いよいよ本題です。果たして、傾向スコアを用いると機械的にバイアスの補正ができるのでしょうか？
傾向スコアでバイアス除去してみよう
傾向スコアを求めた上で、IPW 推定量を求めます。傾向スコアの算出には、ロジスティック回帰を用います。
from sklearn.linear_model import LogisticRegression
以下、傾向スコアの計算から ATE, ATT, ATU の算出までを実行するクラスです。先に作成したユーザリストのデータフレームで初期化し、self.fit で傾向スコア算出、self.estimate で ATE などを求めます。fit の引数には傾向スコア算出時に使う共変量を設定してください。
class PropensityScore:
def __init__(selef, df):
selef.df_data = df.copy()
def fit(self, X_vars):
#calculate propensity score 'e' by Logistic Regression
self.model = LogisticRegression()
self.X = pd.get_dummies(self.df_data[X_vars]).values
self.y = self.df_data["x_flag"].values
self.model.fit(self.X, self.y)
self.df_data["e"] = self.model.predict_proba(self.X)[:, np.where(self.model.classes_ == 1)].flatten()
def estimate(self):
#calculate 'IPWE'(Inverse Probability Weighting Estimator)
w1 = self.df_data["x_flag"] / self.df_data["e"]
w0 = (1 - self.df_data["x_flag"]) / (1 - self.df_data["e"])
self.E1 = (np.sum(self.df_data["x_rate"] * w1)) / np.sum(w1) # E(y1)
self.E0 = (np.sum(self.df_data["x_rate"] * w0)) / np.sum(w0) # E(y0)
#calculate other expectation values
w1 = self.df_data["x_flag"]
w0 = (1 - self.df_data["x_flag"])
self.E11 = (np.sum(self.df_data["x_rate"] * w1)) / np.sum(w1) # E(y1|z=1)
self.E00 = (np.sum(self.df_data["x_rate"] * w0)) / np.sum(w0) # E(y0|z=0)
self.p1 = np.sum(w1) / np.sum(w0 + w1) # p(z=1)
self.p0 = np.sum(w0) / np.sum(w0 + w1) # p(z=0)
self.E10 = (self.E1 - self.E11 * self.p1) / self.p0 if self.p0 != 0 else None # E(y1|z=0)
self.E01 = (self.E0 - self.E00 * self.p0) / self.p1 if self.p1 != 0 else None # E(y0|z=1)
self.ATE = self.E1 - self.E0 # E(y1 - y0)
self.ATT = self.E11 - self.E01 # E(y1 - y0|z=1)
self.ATU = self.E10 - self.E00 # E(y1 - y0|z=0)
では、計算してみましょう。
ps = PropensityScore(df)
ps.fit(["age", "gender"])
ps.estimate()
結果を表示します。本題のバイアス補正済み広告効果 (2.1) は (4.9) に従って計算されて
print("ATE = E(y1 - y0) = {0}".format(round(ps.ATE, 3)))
となりました。すばらしい、0.2 ポイントきましたね。
ついでに ATT, ATU も表示してみます。
print("ATT = E(y1 - y0|z=1) = {0}".format(round(ps.ATT, 3)))
print("ATU = E(y1 - y0|z=0) = {0}".format(round(ps.ATU, 3)))
う～む。正直、この結果はどう捉えるべきかよく分からないですね。
では、このへんで。