本記事は効果検証入門3章の備忘録となります。
概要
傾向スコアとは
介入グループと非介入グループのデータの性質を近くすることで、介入の効果を明らかにする方法です。介入が割り振られる傾向スコアが同一のユーザの中ではランダムに決めてるに等しいと言う考えロジスティック回帰を用いて介入変数を目的変数として計算する
活用用途について
- 介入がランダムに行なわれる場合。
傾向スコアは介入を確率で求め類似する性質のものを比較する方法です。しかし、実際の分析時にはメール配信の意思決定にランダム性がないことが考えられます。例えば予測値が一定以上の場合メールが配信されるが、一定以下の場合配信されないなど。この場合、傾向スコアは常に1か０になるためここで紹介した方法での分析はできません。このような場合５章で説明するような回帰不連続デザインと言う方法を利用します。
www.case-k.jp
- 学習に必要な過去のデータが十分にある場合
介入確率を算出する上で十分な過去データが必要になります。なので一定量のデータがない場合適用するのは難しいです。
傾向スコアマッチング
得られた傾向スコアを利用してサンプル同士をマッチングさせる方法。セレクションバイアスがなくなるようなサンプルのペアそれぞれの効果を算出して平均する方法です。
IPW
傾向スコアをサンプルの重みとして利用する逆確率重み付き推定。傾向スコアをサンプルの重みとして利用して、与えられたデータ全体での介入を受けた場合の結果の期待値と介入を受けなかった場合の期待値を推定します。IPWはグループ間の傾向スコアの偏りを重みで補正します。傾向スコアが高いデータ(P(X)=0.7)は必然的に介入が行なわれたデータ(Z=1)が多く、傾向スコアが低いデータ*1は加入が行なわれなかったデータ(Z=0)が多く含まれます。この時仮に傾向スコア(P(X))と目的変数(Y)に正の相関があった場合、目的変数の値が小さいデータほど介入データは含まれないことになります。それぞれのグループ内で目的変数の期待値を計算すると、介入が行なわれたグループで算出した平均の期待値はより大きくなり、介入が行なわれなかったグループは平均の期待値が小さくなります。IPWはこのような問題に対して重み付きの平均を求めることで対処しています。バイアスのある平均に対して、サンプルサイズ(n)を傾向スコアの逆数(1/傾向スコア)倍することで、バイアスを取り除いたヒストグラムを実現します。
実践
事前準備
事前にデータセットを準備します。
# install package library("tidyverse") # read data email_data <- read_csv("http://www.minethatdata.com/Kevin_Hillstrom_MineThatData_E-MailAnalytics_DataMiningChallenge_2008.03.20.csv")
データ項目を確認してみます。
- recency：最後の購入からの経過月数
- history_segment：昨年の購入額の階層
- history：昨年の購入額
- mens：昨年に男物の商品を購入しているか
- womens：昨年に女物の商品を購入しているか
- zip_code ：地区分類コード
- newbie：過去12ヶ月に新しくユーザになったか
- channel：昨年に置いてどのチャンネルから購入したか
- segment：どのメールが配信されたか
- visit：メールが配信されてから２週間いないにサイトに訪問したか
- conversion：メールが配信されてから２週間いないに購入したか
- spend：購入した際の購入額
head(email_data,5) # A tibble: 5 x 12 recency history_segment history mens womens zip_code newbie channel segment visit conversion spend <dbl> <chr> <dbl> <dbl> <dbl> <chr> <dbl> <chr> <chr> <dbl> <dbl> <dbl> 1 10 2) $100 - $200 142. 1 0 Surburban 0 Phone Womens E-Mail 0 0 0 2 6 3) $200 - $350 329. 1 1 Rural 1 Web No E-Mail 0 0 0 3 7 2) $100 - $200 181. 0 1 Surburban 1 Web Womens E-Mail 0 0 0 4 9 5) $500 - $750 676. 1 0 Rural 1 Web Mens E-Mail 0 0 0 5 2 1) $0 - $100 45.3 1 0 Urban 0 Web Womens E-Mail 0 0 0
介入は男性向けメールと女性向けメールがあるので、まず女性向けメールのレコードを取り除きます。
男性向けのメールが配信されたサンプルと配信されなかったサンプルに絞り込みます。(後から色々と検証できるのでログは残しておいた方が良さそうです)
# remove delivery message for women and add columns of delivery message or not male_df <- email_data %>% filter(segment != "Womens E-Mail") %>% mutate(treatment= if_else(segment == "Mens E-Mail", 1, 0))
女性に対して配信されたメールが削除されて、施策の介入有無を意味するカラムが追加されていることを確認できます。
head(male_df, n=5) # A tibble: 5 x 13 recency history_segment history mens womens zip_code newbie channel segment visit conversion spend treatment <dbl> <chr> <dbl> <dbl> <dbl> <chr> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> 1 6 3) $200 - $350 329. 1 1 Rural 1 Web No E-Mail 0 0 0 0 2 9 5) $500 - $750 676. 1 0 Rural 1 Web Mens E-Mail 0 0 0 1 3 9 5) $500 - $750 675. 1 1 Rural 1 Phone Mens E-Mail 0 0 0 1 4 2 2) $100 - $200 102. 0 1 Urban 0 Web Mens E-Mail 1 0 0 1 5 4 3) $200 - $350 241. 0 1 Rural 1 Multichannel No E-Mail 0 0 0 0
バイアスのあるデータを準備します。
biased_data <- male_df %>% mutate(obs_rate_c=if_else( (history > 300) | (recency < 6) | (channel == "Multichannel"), obs_rate_c,1), obs_rate_t=if_else( (history > 300) | (recency < 6) | (channel == "Multichannel"), 1, obs_rate_t), random_number=runif(n=NROW(male_df))) %>% filter((treatment==0 & random_number < obs_rate_c)| (treatment==1 & random_number < obs_rate_t))
中身を確認しておく
head(biased_data,5) # A tibble: 5 x 16 recency history_segment history mens womens zip_code newbie channel segment visit conversion spend treatment obs_rate_c obs_rate_t <dbl> <chr> <dbl> <dbl> <dbl> <chr> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> 1 6 3) $200 - $350 329. 1 1 Rural 1 Web No E-M… 0 0 0 0 0.5 1 2 9 5) $500 - $750 676. 1 0 Rural 1 Web Mens E… 0 0 0 1 0.5 1 3 9 5) $500 - $750 675. 1 1 Rural 1 Phone Mens E… 0 0 0 1 0.5 1 4 2 2) $100 - $200 102. 0 1 Urban 0 Web Mens E… 1 0 0 1 0.5 1 5 4 3) $200 - $350 241. 0 1 Rural 1 Multic… No E-M… 0 0 0 0 0.5 1
ps_model Call: glm(formula = treatment ~ recency + history + channel, family = binomial, data = biased_data) Coefficients: (Intercept) recency history channelPhone channelWeb 0.9854765 -0.1248497 0.0006195 -0.3093588 -0.2862948 Degrees of Freedom: 31896 Total (i.e. Null); 31892 Residual Null Deviance: 44050 Residual Deviance: 41990 AIC: 42000
傾向スコアマッチング
ライブラリの読み込み
install.packages("MatchIt") library("MatchIt")
傾向スコアを用いたマッチング
m_near <- matchit(formula=treatment ~ recency + history + channel, data=biased_data, method="nearest", replace=TRUE )
結果を確認
m_near Call: matchit(formula = treatment ~ recency + history + channel, data = biased_data, method = "nearest", replace = TRUE) Sample sizes: Control Treated All 14791 17106 Matched 8223 17106 Unmatched 6568 0 Discarded 0 0
マッチング後のデータ作成
matched_data <- match.data(m_near)
中身を確認
head(matched_data,5) recency history_segment history mens womens zip_code newbie channel segment visit conversion spend treatment obs_rate_c 1 6 3) $200 - $350 329.08 1 1 Rural 1 Web No E-Mail 0 0 0 0 0.5 2 9 5) $500 - $750 675.83 1 0 Rural 1 Web Mens E-Mail 0 0 0 1 0.5 3 9 5) $500 - $750 675.07 1 1 Rural 1 Phone Mens E-Mail 0 0 0 1 0.5 4 2 2) $100 - $200 101.64 0 1 Urban 0 Web Mens E-Mail 1 0 0 1 0.5 5 4 3) $200 - $350 241.42 0 1 Rural 1 Multichannel No E-Mail 0 0 0 0 0.5 obs_rate_t random_number distance weights 1 1 0.1087145 0.5384078 0.4807085 2 1 0.2857189 0.4985455 1.0000000 3 1 0.4806688 0.4926623 1.0000000 4 1 0.3423975 0.6253792 1.0000000 5 1 0.4290973 0.6537643 1.4421256
マッチング後のデータで効果の推定
matchit()では基本的にATTの推定が行なわれます。
- ATT (Average Treatment effect on Treated)
ATTは介入を受けたサンプルにおける効果となります。そのため平均的な効果を推定した値と結果が異なる可能性が考えられます。
- ATE ( Average Treatment effect)
これはATTのように介入グループに対してのみマッチングが行なわれるのではなく、非介入グループに対しても同様にマッチングを行います。マッチングは計算時間が長くなります。また、説明変数が多くなると完全に傾向スコアが一致するケースは少ないため類似してるものでペアを作る必要があります。
P値の値は小さく、介入の効果を確認してみると、0.888であることが確認できます。
matched_data <- match.data(m_near) head(matched_data,5) PSM_result<-matched_data %>% lm(spend~treatment,data=.)%>% tidy() PSM_result # A tibble: 2 x 5 term estimate std.error statistic p.value <chr> <dbl> <dbl> <dbl> <dbl> 1 (Intercept) 0.710 0.179 3.97 0.0000736 2 treatment 0.795 0.218 3.65 0.000262
IPW
IPWで効果を調べてみます。
ライブラリをインストール
install.packages("WeightIt") library("WeightIt")
重みの推定
weighting <- weightit(formula=treatment ~ recency + history + channel, data=biased_data, method="ps", estimand="ATE")
確認
weighting A weightit object - method: "ps" (propensity score weighting) - number of obs.: 31897 - sampling weights: none - treatment: 2-category - estimand: ATE - covariates: recency, history, channel
重みつきデータでの効果の推定
IPW_result <- lm(data=biased_data, formula=spend~treatment, weights=weighting$weights) %>% tidy()
メール介入で0.569程度の売上向上。(数値が小さすぎるのと書籍と数値が一致してないので間違ってる？)
IPW_result # A tibble: 2 x 5 term estimate std.error statistic p.value <chr> <dbl> <dbl> <dbl> <dbl> 1 (Intercept) 0.825 0.122 6.78 1.21e-11 2 treatment 0.569 0.173 3.30 9.77e- 4
傾向スコアと回帰分析はどちらを使うべきか
Yの値がどのような仕組みで決まるか知識がある場合は、回帰分析を使い、Yに関するドメイン知識が少ない場合には傾向スコアを用いるのが望ましい。
回帰分析のメリットとデメリットは次の通りです。
メリット
- 手軽で取り組みやすい
デメリット
- 入念なモデリングが必要
傾向スコアのメリットとデメリットは次の通りです。
メリット
デメリット
- 計算時間
set.seed(1) train_flag <- sample(NROW(male_df), NROW(male_df)/2, replace = FALSE) male_df_train<-male_df[train_flag,] %>% filter(treatment==0) mael_df_test<-male_df[-train_flag,] predict_model <-glm( data=male_df_train, formula=conversion~recency+history_segment+channel+zip_code, family=binomial ) pred_cv<-predict(predict_model, newdata=male_df_test, type="response") pred_cv_rank<-percent_rank(pred_cv) pred_cv_rank
所感
傾向スコアに置いては、モデルの推定確率を元に分類するためモデルの推定精度は重要になるように思いました。
実務で使える機会はあまりなさそうな気がしましたが、介入がランダムに割り当てられる事象があった場合の効果を見るのに使えそうです。
*1:P(X)=0.3
