「効果検証入門」を読む　～第2章介入効果を測るための回帰分析 ③回帰分析の導入の分析例～

# (1) 必要なパッケージのインポート
import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
from urllib.request import urlopen
from io import StringIO

# (3) データの読み込み
url = "http://www.minethatdata.com/Kevin_Hillstrom_MineThatData_E-MailAnalytics_DataMiningChallenge_2008.03.20.csv"
email_data = pd.read_csv(url)

# (4) 女性向けメールが配信されたデータを除外し、treatment列を追加
male_df = email_data[email_data["segment"] != "Womens E-Mail"].copy()
male_df["treatment"] = np.where(male_df["segment"] == "Mens E-Mail", 1, 0)

# (5) セレクションバイアスのあるデータを作成
np.random.seed(1)
obs_rate_c = 0.5
obs_rate_t = 0.5

# 条件に応じて観測率を変更
male_df["obs_rate_c"] = np.where(
    (male_df["history"] > 300) | (male_df["recency"] < 6) | (male_df["channel"] == "Multichannel"),
    obs_rate_c, 1)
male_df["obs_rate_t"] = np.where(
    (male_df["history"] > 300) | (male_df["recency"] < 6) | (male_df["channel"] == "Multichannel"),
    1, obs_rate_t)
male_df["random_number"] = np.random.uniform(size=male_df.shape[0])

biased_data = male_df[
    ((male_df["treatment"] == 0) & (male_df["random_number"] < male_df["obs_rate_c"])) |
    ((male_df["treatment"] == 1) & (male_df["random_number"] < male_df["obs_rate_t"]))
].copy()

回帰式

本節で対象としているモデルは，以下のようなモデルである。

$\begin{align} Spend_i = \beta_0 + \beta_{treatment} treatment_i + \beta_{history} history_i \\ \end{align}$

ただし，

$Spend$ : ユーザの購入額(目的変数)
$treatment$ : メールの有無(介入変数)
$history$ : 過去の購入額(共変量)

である。

バイアスのあるデータでの回帰分析

バイアスのあるデータ(biased_data)で回帰分析を行ない，介入の効果 $\beta_{treatment}$ を確認する。モデルには，介入変数に加えて共変量も加えている。

サンプルコード

# (6) バイアスのあるデータでの回帰分析
biased_reg = smf.ols("spend ~ treatment + history", data=biased_data).fit()
print(f'{biased_reg.params["treatment"]:.4f}')

出力

0.7293

2.2 回帰分析におけるバイアス

2.2.1 共変量の追加による効果への作用

バイアスのないデータでの回帰分析

共変量の追加による効果への作用を確認するために，まず理想的なデータである，バイアスのないデータmale_dfによる回帰分析の結果を確認する。
今回扱うモデルは，共変量のないモデル

$\begin{align} Spend_i = \beta_0 + \beta_{treatment} treatment_i \\ \end{align}$

である。

サンプルコード

# (7) RCTデータでの回帰分析とバイアスのあるデータでの回帰分析の比較
rct_reg = smf.ols("spend ~ treatment", data=male_df).fit()
print(f'{rct_reg.params["treatment"]:.4f}')

出力

0.7698

バイアスのあるデータでの回帰分析

次に，先ほどと同様に共変量のないモデルを，バイアスのあるデータbiased_dataを用いて回帰分析を行なう。

サンプルコード

nonrct_reg = smf.ols("spend ~ treatment", data=biased_data).fit()
print(f'{nonrct_reg.params["treatment"]:.4f}')

出力

0.8296

このように，バイアスのあるデータでは，介入の効果が大きく見積もられることが確認できた。

共変量を加えた回帰分析

そして，バイアスのあるデータを対象に，共変量を加えたモデルによる回帰分析を行なう。

サンプルコード

nonrct_mreg = smf.ols("spend ~ treatment + recency + channel + history", data=biased_data).fit()
print(f'{nonrct_mreg.params["treatment"]:.4f}')

出力

0.6272

共変量を加えることにより，バイアスのないデータによる結果に近づくことが確認できた。

2.2.3 OVBの確認

本節では，脱落変数バイアス(Omitted Variable Bias : OVB)の効果について確認する。以下のような2つの回帰式を考える。

$\begin{align} &Spend_i = \alpha_0 + \alpha_1 treatment_i + \alpha_2 recency_i + \alpha_3 channel_i + e_i \text{(モデルA)} \\ \\ &Spend_i = \beta_0 + \beta_1 treatment_i + \beta_2 recency_i + \beta_3 channel_i + \beta_4 history_i + u_i \text{(モデルB)} \\ \end{align}$

モデルAでは，モデルBに含まれている変数 $history$ が無いため，OVBが発生している。このOVBの式は，

$\begin{align} \beta_1 - \alpha_1 = \gamma_1 \beta_4 \\ \end{align}$

で表される。ただし $\gamma_1$ は，以下の回帰式(モデルC)によって推定された値である。

$\begin{align} history_i = \gamma_0 + \gamma_1 treatment_i + \gamma_2 recency_i + \gamma_3 channel_i + \varepsilon_i \text{(モデルC)} \\ \end{align}$

計算によって， $\gamma_1 \beta_4$ と $\beta_1 - \alpha_1$ を比較する。

サンプルコード

# (8) OVBの確認
# モデルA
short_model = smf.ols("spend ~ treatment + recency + channel", data=biased_data).fit()
alpha_1 = short_model.params["treatment"]

# モデルB
long_model = smf.ols("spend ~ treatment + recency + channel + history", data=biased_data).fit()
beta_1 = long_model.params["treatment"]
beta_4 = long_model.params["history"]

# モデルC
omitted_model = smf.ols("history ~ treatment + channel + recency", data=biased_data).fit()
gamma_1 = omitted_model.params["treatment"]

OVB = beta_4 * gamma_1
coef_gap = alpha_1 - beta_1
print(f'OVB(gamma_1 * beta_4):{OVB:.4f}')
print(f'alpha_1 - beta_1):{coef_gap:.4f}')

出力

OVB(gamma_1 * beta_4):0.0365
alpha_1 - beta_1):0.0365

上記のように， $\gamma_1 \beta_4$ と $\beta_1 - \alpha_1$ の値が同じであることが確認できた。

まとめと感想

今回は，「第2章介入効果を測るための回帰分析」における，回帰分析の導入の分析例についてまとめた。

共変量を加えることによって，データに含まれるバイアスの影響が低減できることが確認できた。

また実装においては，Pythonのstatsmodels.formula.apiを用いることで，Rにおける記法とかなり似た書き方で線形回帰が実装できることが分かった。RのコードをPythonで置き換えるときには活用していきたい。

本記事を最後まで読んでくださり，どうもありがとうございました。

参考サイト

ランキング参加中

統計

jiku log

データサイエンスの核心を掴む : 学びと発見の記録

「効果検証入門」を読む　～第2章介入効果を測るための回帰分析 ③回帰分析の導入の分析例～

はじめに

本記事を読むことで得られること

目次

2.1 回帰分析の導入

2.1.5 メールマーケティングデータの分析(回帰編)

バイアスのあるデータの作成

回帰式

バイアスのあるデータでの回帰分析

2.2 回帰分析におけるバイアス

2.2.1 共変量の追加による効果への作用

バイアスのないデータでの回帰分析

バイアスのあるデータでの回帰分析

共変量を加えた回帰分析

2.2.3 OVBの確認

まとめと感想

参考サイト