Python3
統計学
statsmodels

クリスマス(季節性)に負けない統計的因果推論

Retty Advent Calendar2018の24日目の記事です。
先日は @shindo-taichi さんのRettyマネタイズを支えるプログラマティック広告運用 という記事でした。

24日という美味しい枠を貰いました、分析チーム19卒内定者の@wtnVengaと申します!
「Rettyの分析チームの内情を知りたい」 という方は、是非コチラの記事をご覧ください:santa_tone1:

はじめに

Rettyは生活の基本となる「衣・食・住」のうち「食」を扱うため、ほぼ常に季節性が加わった観測データを扱う事になります。特に クリスマス忘年会のある12月は飲食店需要の高まりから、目算で施策などの良し悪しを判断することは不可能になります。

しかし「季節性のせいで効果検証は行えません」などの言葉を発するのは、分析者として髀肉之嘆の限り。 「季節性に負けた、すなわちクリスマスに負けた」 と言っても過言ではないでしょう。1 季節性に負けない検証で、意思決定をサポートすること が分析者には求められます。

そこで今回はRettyの分析チームが重宝している、季節性を考慮できて便利な DID分析についてご紹介します。

TL; DR

  • 季節性を考慮した検証において、簡単便利な DID分析 をご紹介。
  • DID分析は 平行トレンド仮定という仮定を満たせば心強い分析手法。
  • マサカリプレゼント をお待ちしております:santa_tone2:

DID分析とは?

"DID(Difference in Difference)分析" とは、その名の通り「差分の差分法分析」になります。
以下の図のように、施策群と比較群それぞれの前後差から施策効果を検証する分析手法です。
スクリーンショット 2018-12-20 23.17.22.png

また一見すると扱い易く見えますが、 以下のような強い仮定のもとで成立する分析手法です。

  • 平行トレンド仮定…施策実施がないとき、両群ともに観測データは同傾向のまま
  • 共通ショック仮定…期間前後において、両群ともに施策以外に同様の処置のみ受けている
  • スピルオーバー効果がない…両群間に施策影響の波及・漏洩が起きていない

一応は平行トレンド仮定を満たしていれば、季節性などもはや敵ではありませんが、この仮定を満たす証明が困難です。しかし有難い事に非平行トレンド時の検証法も数多く提案されているため2、まずは基礎のDID分析を知ることで、より季節性に強い分析者としての端緒になります。

DID分析の実践

DID分析は基本的には回帰式で表現され、ダミー変数の群種($D_g$)・期間($D_t$)が説明変数になります。3

 y = β_0 + β_1D_g + β_2D_t + β_3D_gD_t + ε

モデリング時には、前提におく仮定と効果差分の種類によって手法が変わります。今回は上2種を実践します。

  • 平行トレンド仮定・数量の効果差分 → OLS回帰
  • 平行トレンド仮定・比率の効果差分 → GLMロジスティック回帰(二項分布)45
  • 非平行トレンド仮定 → 共変量を説明変数に追加、ロジット関数により推計した処置率の差を加重し推計6、 etc...

(以降の検証は全てコチラのcolaboratory資料 から実践できます。)

数量を目的変数とする場合

「Rettyアプリで新規フォローされる週間UU数を増やしたい。施策として4種のユーザー区分において各50%のユーザーにキャンペーン通知を表示した。結果としてキャンペーン通知は効果的であったか」という例で効果検証します。

#サンプルデータの準備
import pandas as pd
import numpy as np
np.random.seed(1224)

outcome = []
for i in range(4): 
  outcome.extend([ 
      #before
      np.random.normal(5000,100),
      np.random.normal(5000,100),

      #after :両群ともに倍増し対象群の方がさらに10ptの施策効果があった
      np.random.normal(10000,100),
      np.random.normal(10000,100) *1.1
  ])

params = pd.DataFrame({
    'Treated' : np.array([0,1,0,1]*4), # D_t 項: [比較群, 対象群] で [0,1,0,1...] 
    'Period' : np.array([0,0,1,1]*4), # D_p項 : [前期間, 前期間, 後期間, 後期間] で [0,0,1,1,...]
    'DID' : np.array([0,1,0,1]*4) * np.array([0,0,1,1]*4) # D_t*D_g項
})
#施策効果(平均効果差分)の計算
join_df = pd.concat([pd.Series(outcome), params], axis=1)
avg_df = join_df.groupby(['Treated', 'Period'], as_index=False).mean()

G1 = avg_df.query('Treated == 1 and Period == 1').iloc[:, 2].values - avg_df.query('Treated == 1 and Period == 0').iloc[:, 2].values
G0 = avg_df.query('Treated == 0 and Period == 1').iloc[:, 2].values - avg_df.query('Treated == 0 and Period == 0').iloc[:, 2].values

print(G1 - G0)

> [974.36110163]
#OLS回帰による検証
import statsmodels.api as sm
lm = sm.OLS(outcome,  sm.add_constant(params)).fit()
print(lm.summary())

> OLS Regression Results
...
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const       4957.9688     68.284     72.608      0.000    4809.190    5106.748
DID          974.3611    136.569      7.135      0.000     676.803    1271.919
Period      5050.1552     96.569     52.296      0.000    4839.750    5260.560
Treated       26.0136     96.569      0.269      0.792    -184.392     236.419
==============================================================================

以上から、結果として施策効果(平均効果差分)は約974UUの差がつき、P>|t|列の値からDID変数が(もし有意水準0.05と設定するならば)統計的に有意であり、この施策は結果に関係がある可能性が見えてきました。
一方で群種を意味するTreated変数が有意でないため、ここでは予定通りPeriodとDIDの2つの変数が着眼点になります。7

比率を目的変数とする場合

「『〇〇が美味しいお店20選』というページの直帰率を下げたい。施策として4種のページ上部に50%の確率でランダムに出し分けるキャンペーンバナーを追加した。結果としてキャンペーンは効果的であったか」という例で効果検証します。

import pandas as pd
import numpy as np
np.random.seed(1224)

#二項分布を仮定した一般化線形モデルを用いるため、[success, fail]の目的変数を作成。
outcome_n = [10000,10000,15000,15000]*4

outcome_s = []
for i in range(4): 
  outcome_s.extend([ 
      #before
      np.random.binomial(10000, 0.4,1)[0],
      np.random.binomial(10000, 0.4,1)[0],      #施策前、直帰率40%で開始

      #after :両群ともに1.5倍増し、対象群は2pt直帰率が減少した
      np.random.binomial(15000, 0.4,1)[0],
      np.random.binomial(15000, 0.38,1)[0],
  ])

outcome =[]
for n, s in zip(outcome_n, outcome_s):
    outcome.append([s, n-s]) #[success, all - success]

params = pd.DataFrame({
    'Treated' : np.array([0,1,0,1]*4), 
    'Period' : np.array([0,0,1,1]*4), 
    'DID' : np.array([0,1,0,1]*4) * np.array([0,0,1,1]*4)
})
# 施策効果(平均効果差分)の計算
join_df = pd.concat([pd.Series(outcome), params], axis=1)
rate_li = []
for i in join_df.iloc[:,0]: rate_li.append(i[0] / (i[0] + i[1]))
join_df_cal = pd.concat([pd.Series(rate_li), join_df], axis=1)

avg_df = join_df_cal.groupby(['Treated', 'Period'], as_index=False).mean()
G1 = avg_df.query('Treated == 1 and Period == 1').iloc[:, 2].values - avg_df.query('Treated == 1 and Period == 0').iloc[:, 2].values
G0 = avg_df.query('Treated == 0 and Period == 1').iloc[:, 2].values - avg_df.query('Treated == 0 and Period == 0').iloc[:, 2].values

print(G1- G0)
> [-0.012475]
#GLMロジスティック回帰による検証
import statsmodels.api as sm

#statmodelsのGLMには切片が必要なので、 sm.ad_constant()で追加。
glm = sm.GLM(outcome, sm.add_constant(params), family=sm.families.Binomial()).fit()
print(glm.summary())

> Generalized Linear Model Regression Results 
...
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.3929      0.010    -38.542      0.000      -0.413      -0.373
DID           -0.0526      0.019     -2.817      0.005      -0.089      -0.016
Period        -0.0236      0.013     -1.795      0.073      -0.049       0.002
Treated       -0.0105      0.014     -0.728      0.466      -0.039       0.018
==============================================================================

以上から、結果として施策効果(平均効果差分)は約1.2%の直帰率の差が見え、この場合でもP>|z|列の値から(もし有意水準0.05と設定するならば)DID変数が統計的に有意であり、施策が結果に影響があると考えられます。
またPeriodとTreated変数は有意でない点から、施策効果の可能性がより見えてきます。

おわりに

ざっくり説明ですが、簡単便利かつ奥深いDID分析での検証方法について実践しました。
前述の通りこの検証方法には強い仮定があり、またそれ以前にデータが適切に設計・取得されている大前提があります。

「アナリスト・エンジニア・プランナーの3者で密なコミュニケーションが取れる結果、ビジネスの場でより簡潔かつ正確にデータ分析ができる」というのがRetty分析チームの文化の1つなので、ご興味のある方は是非以下もお目通しください

それでは、皆さまからのマサカリプレゼントをお待ちしております。:santa_tone2:


主な参考元


  1. 過言です。 

  2. 例えば以下。
    Mora, R., & Reggio, I., "Treatment Effect Identification Using Alternative Parallel Assumptions", https://e-archivo.uc3m.es/bitstream/handle/10016/16065/we1233.pdf?sequence=1, 2012
    Mora, R., & Reggio, I., "Alternative diff-in-diffs estimators with several pretreatment periods", https://www.tandfonline.com/doi/abs/10.1080/07474938.2017.1348683, 2017 

  3. 期間・施策の項はともに連続値を取ることも可能です。 

  4. Donna Spiegelman Ellen Hertzmark, "Easy SAS Calculations for Risk or Prevalence Ratios and Differences", https://academic.oup.com/aje/article/162/3/199/171116, 2004 

  5. Guangyong Zou, "A Modified Poisson Regression Approach to Prospective Studies with Binary Data", https://academic.oup.com/aje/article/159/7/702/71883, 2004 

  6. Alberto Abadie, "Semiparametric Difference-in-Differences Estimators", https://economics.mit.edu/files/11869, 2005 

  7. 勿論その他の変数も統計的に有意であることも解釈時に考慮すべきです。期間差はそもそもどの程度なのか、群間差を連続値で表した変数の追加、ダミー変数の多重共線性の問題などなど…