Causal Impact -paper summary-

A summary of the paper: Brodersen KH, Gallusser F, Koehler J, Remy N, Scott SL. Inferring causal impact using Bayesian structural time-series models. Ann Appl Stat 2015;9. https://doi.org/10.1214/14-AOAS788.

Updated: Oct. 2022

Transcript

1. @Maxwell_110
   INFERRING CAUSAL IMPACT
   USING BAYESIAN STRUCTURAL
   TIME-SERIES MODELS
   Paper Summary
   Ver. 2022
   

   View Slide

2. 01
   Introduction
   02
   Bayesian
   structual
   time-series
   models
   03
   Application
   to
   simulated
   data
   04
   Application
   to
   empirical
   data
   

   View Slide

3. 01
    Causal Impact とは，Google による因果推論のフレームワーク
    A/B テストなどの RCT ができない環境下において，反実仮想的予測を行い介入効果（(A)TT: (Average)
   Treatment Effect on the Treated）を推定する
    時系列ベイズ状態空間モデル：
   Simulation Smoother + Gibbs sampling にて，状態 および パラメータ をサンプリング
    局所的トレンド・季節性・共変量などを組み込み可能
    Spike and slab 事前分布による共変量の選択機能
    Bayesian Model Averaging による過学習の抑制
   Introduction
   Resources:
   https://research.google/pubs/pub41854/
   https://youtu.be/GTgZfCltMm8
   http://google.github.io/CausalImpact/
   Summary
   

   View Slide

4. Causal Impact は
   「反実仮想（counterfactual）の時系列を
   扱えるように DID を一般化したもの」と謳っている
   Causal Effect
   in DID (traditional DID)
   Difference In Difference (DID) Causal Impact (CI)
   post-intervention
   pre-intervention post-intervention
   pre-intervention
   Causal Effect at any given time point
   Counterfactual time series
   Observed time series
   Treatment group
   Control group
   • Parallel trends assumption
   • Synthetic control (Abadie et al., 2010)
   𝑌𝑖
   = 𝛽0
   + 𝛽1
   𝑇𝑟𝑒𝑎𝑡𝑖
   + 𝛽2
   𝑃𝑜𝑠𝑡𝑖
   + 𝜷𝟑
   𝑇𝑟𝑒𝑎𝑡𝑖
   × 𝑃𝑜𝑠𝑡𝑖
   + 𝛽4
   𝐶𝑜𝑣𝑎𝑟𝑖𝑎𝑡𝑒𝑠𝑖
   + 𝑒𝑖
   実際には，DID の中でも，特に Synthetic Control に inspire されているといって良さそう
   論文中で使用される "control" という word は，普通の統制群の意味ではなく，この synthetic control の control の意に近いと思う
   

   View Slide

5. Synthetic Control
    Causal Impact と目的を同じくする手法で，Causal Impact は Synthetic Control に inspire されているといっても
   良さそう
    A/B test 等の場合と異なり非介入群が存在しない場合，介入をうけていない reference data から「任意の時系列の
   加重平均で求まる，介入群と似た性質をもつ時系列」を生成して比較に使用する手法 => Causal Impact と目的は同
   じで，統制群がなくても impact を推定したい時に使う
    例えば，Abadie and Gardeazabal, 2003 はスペインの Basque 州において発生したテロが経済に与えた影響を，
   synthetic control で推定した．その際に，2 つの別のスペインの地域の経済の時系列を組み合わせて，統制群を作
   成した．他にも US の California における禁煙キャンペーンの影響を synthetic control で推定した例もある
   (Abadie et al., 2010)
    擬似統制群の構築に使用可能な情報
   1. 介入前の目的変数そのもの（トレンド・季節性など）
   2. 他の時系列
   - 介入の影響が無い
   - 介入前の目的変数に対する予測力がある
   3. 事前分布に使用可能な知識（過去の研究結果など）
   これらは，Causal Impact でも同様に使用可能
   Using units from the
   rest of the U.S. to
   make a synthetic
   control.
   A synthetic control
   matches to the
   interved time-seriese.
   Abadie et al., 2010
   (Abadie et al., 2010; Abadie and Gardeazabal, 2003)
   

   View Slide

6. Limitations for Difference In Difference
    通常の DID は，data が i.i.d. であることを前提とした静的な回帰モデルを
   使用するため，推定した係数の標準誤差の過小評価に伴う Type I error が大
   きい
    大部分の DID は，介入前後の二時点のみの比較であり，介入効果がその後ど
   うなっていくかを分析しない（できない）
    DID の枠組み内で，統制群が存在しない場合に使用可能な "Synthetic
   Control" は，統制群の構成の仕方に数理的な制約がある．詳しくは以下参照．
   (Abadie et al., 2010; Abadie and Gardeazabal, 2003)
   

   View Slide

7. 行列表示ではなく，季節性をもつ場合の式を陽に表すと・・・
   𝑦𝑡
   = 𝜇𝑡
   + 𝛾𝑡
   + 𝛽𝑇 𝑥𝑡
   + 𝜀𝑡
   𝜇𝑡
   = 𝜇𝑡−1
   + 𝛿𝑡−1
   + 𝜂𝜇, 𝑡−1
   𝛿𝑡
   = 𝛿𝑡−1
   + 𝜂𝛿, 𝑡−1
   𝛾𝑡
   = − 𝑠=1
   𝑆−1 𝛾𝑡−𝑠
   + 𝜂𝛾, 𝑡−1
   状態空間モデルの一般式（行列表示）
   観測方程式： 𝑦𝑡
   = 𝑍𝑡
   𝑇 𝛼𝑡
   + 𝜀𝑡
   𝜀𝑡
   ~ 𝑁 0, 𝜎𝑡
   2
   状態方程式： 𝛼𝑡+1
   = 𝑇𝑡
   𝛼𝑡
   + 𝑅𝑡
   𝜂𝑡
   𝜂𝑡
   ~ 𝑁 0, 𝑄𝑡
   Causal Impact の実体は状態空間モデル
    AR や ARIMA などの時系列モデルより表現力が高い
   （例えば AR や ARIMA は状態空間モデルの特殊形）
    計算コストが低い
   bsts は simulation smoother (Durbin et al., 2002) + Gibbs sampling で推定
   状態方程式
   観測方程式
   線形ガウス状態空間モデル
   02
   Bayesian structural time-series models
   𝑦𝑡
   : observed data
   𝛼𝑡
   : state vector
   𝑍𝑡
   : output vector
   𝑇𝑡
   : transition matrix
   𝑅𝑡
   : control matrix
   𝜀𝑡
   : scalar observation error
   𝜎𝑡
   2: noise variance
   𝜂𝑡
   : system error
   𝑄𝑡
   : state-diffusion matrix
   

   View Slide

8. Basic structural model
   Modular
   ・ 線形モデルであるため，コンポーネント 2 や 3 を加え
   たり外したりしてモデル構築が可能
   Fully Bayesian Approach
   ・ 全コンポーネントをベイズの枠組みの中で計算可能
   𝛾𝑡
   = − 𝑠=1
   𝑆−1 𝛾𝑡−𝑠
   + 𝜂𝛾, 𝑡
   𝜂𝛾, 𝑡
   ~ 𝑁 0, 𝜎𝛾
   2
   周期 S の季節性
   春夏秋冬であれば S = 4 でいずれか一つの季節は reference（係数は 0）
   𝛾𝑡+1
   𝛾𝑡
   𝛾𝑡−1
   =
   −1 −1 −1
   1 0 0
   0 1 0
   𝛾𝑡
   𝛾𝑡−1
   𝛾𝑡−2
   +
   𝜂𝛾, 𝑡
   0
   0
   静的（static）・動的（dynamic）な係数をもつ共変量
   「(spill-over 効果なども含めた) 介入の影響を受けない」共変量 𝑋𝑡
   や係数 𝛽 が時間とともに変化する「動的係数 𝛽𝑡
   」の設定が可能
   𝛽𝑇𝑋𝑡
   → 𝛽𝑡
   𝑇𝑋𝑡
   =
   𝑗=1
   𝐽 𝛽𝑗, 𝑡
   𝑋𝑗, 𝑡
   𝛽𝑗, 𝑡
   = 𝛽𝑗, 𝑡−1
   + 𝜂𝛽, 𝑗, 𝑡−1
   𝜂𝛽, 𝑗, 𝑡
   ~ 𝑁 0, 𝜎𝛽𝑗
   2
   Random walk →
   dynamic
   ローカル線形トレンド（確率的トレンド）
   トレンド成分 𝛿𝑡
   が変化可能．例えば，最初は増加でその後，減少トレンドに転じるような場合も表現可能
   𝜇𝑡
   = 𝜇𝑡−1
   + 𝛿𝑡−1
   + 𝜂𝜇, 𝑡−1
   𝜂𝜇, 𝑡
   ~ 𝑁 0, 𝜎𝜇
   2
   𝛿𝑡
   = 𝛿𝑡−1
   + 𝜂𝛿, 𝑡−1
   𝜂𝛿, 𝑡
   ~ 𝑁 0, 𝜎𝛿
   2
   𝑦𝑡
   = 𝜇𝑡
   + 𝛾𝑡
   + 𝛽𝑇𝑥𝑡
   + 𝜀𝑡
   1 2 3
   1
   2
   3
   Random walk →
   

   View Slide

9. 1. 介入前において，共変量と介入群の関係性が安定 => 静的係数
   => Forward-filtering と backward-sampling の枠組みの中で，spike-and-slab 事前分布（後述）を効率的に実装
   できるため，変数選択による過学習の抑制が期待できる
   2. 時間の経過に伴い共変量と介入群の線形関係が変化すると考えられる => 動的係数
   => 共変量候補が多いと計算コストが大きくなってしまうが，いくつか代替案を述べている (Nakajima and West,
   2013)
   => しかし，bsts の CRAN manual の bsts::add.dynamic.regression をみる限りでは計算コストの低減対策はさ
   れておらず，他の状態変数と同様の扱いの模様
   介入効果を受けない共変量を多く用意すること自体が難しいので，計算コストはそこまで心配しなくて良い？
   共変量に対して静的係数と動的係数のどちらを選択すべきか？
   

   View Slide

10. Example: Local Linear Trend model with a static regression component
    𝑦𝑡
    = 𝜇𝑡
    + 𝛽𝜌
    𝑇𝑥𝑡
    + 𝜀𝑡
    𝜀𝑡
    ~ 𝑁 0, 𝜎𝑦
    2
    𝜇𝑡
    = 𝜇𝑡−1
    + 𝛿𝑡−1
    + 𝜂𝜇, 𝑡
    𝜂𝜇, 𝑡
    ~ 𝑁 0, 𝜎𝜇
    2
    𝛿𝑡
    = 𝛿𝑡−1
    + 𝜂𝛿, 𝑡
    𝜂𝛿, 𝑡
    ~ 𝑁 0, 𝜎𝛿
    2
    1
    𝜎2
    ~ 𝐺𝑎𝑚𝑚𝑎 𝜈
    2
    , 𝑠
    2
    共変量 𝛽 に関する項を除けば，MCMC でサンプリングするパラメータは観測・状態方程式内の各誤差項の 𝜎2
    𝜎2 の事前分布は，典型的なものとしては，(正規分布に対して共役な) 逆ガンマ分布を採用することが多い
    𝜈 や 𝑠 の値の設定の仕方は論文を参照 (基本的には 𝑦𝑡
    の不偏分散の値を利用して設定する)
    Brodersen et al., 2015
    Figure. 2
    b
    a
    a 状態方程式 (ローカル線形トレンド)
    b 観測方程式
    

    View Slide

11. 共変量の選択機能をもつ spike-and-slab 事前分布
    𝑝 𝛽, 𝜎𝜀
    2, 𝜚 = 𝑝 𝛽𝜚
    |𝜚, 𝜎𝜀
    2 𝑝 𝜎𝜀
    2|𝜚 𝑝 𝜚
    Spike: 𝛽𝑗
    が 0 になるかどうかを表す二項分布
    𝑝 𝜚 =
    𝑗=1
    𝐽 𝜋
    𝑗
    𝜚𝑗 1 − 𝜋𝑗
    1−𝜚𝑗
    Indicator vector: 𝜚 = 𝜚1
    , 𝜚2
    , … , 𝜚𝐽
    𝜚𝑗
    = 1 if 𝛽𝑗
    ≠ 0
    𝜚𝑗
    = 0 if 𝛽𝑗
    = 0
    Slab: 𝛽𝜚
    は 𝜌𝑗
    = 1 の 𝛽𝑗
    の集合
    𝛽𝜚
    | 𝜎𝜀
    2 ~ 𝑁 𝑏𝜚
    , 𝜎𝜀
    2 𝛴𝜚
    1
    𝜎𝜀
    2
    ~ 𝐺𝑎𝑚𝑚𝑎 𝜈𝜀
    2
    , 𝑠𝜀
    2
    観測方程式: 𝑦𝑡
    = 𝜇𝑡
    + 𝛾𝑡
    + 𝛽𝑇𝑥𝑡
    + 𝜀𝑡
    共変量： 𝛽 = 𝛽1
    , 𝛽2
    , … , 𝛽𝐽
    MCMC で 𝛽 の事後分布を sampling するにあたり，
    以下の事前分布を設定
    Spike-and-slab 事前分布
    spike
    slab
    (Gaussian with large variance)
    spike は確率質量関数なので
    実際の幅はゼロ
    𝛽𝑗
    spike
    slab
    𝜋𝑗
    や 𝛴𝜚
    などの値の設定の
    仕方は論文を参照
    

    View Slide

12. Inference
    定常事後分布が 𝑝 𝜃, 𝛼 𝑦1:𝑛
    ) で表される 𝜃, 𝛼 を以下の手続きで得る
    (a) と (b) を交互に繰り返し 𝜃, 𝛼 のサンプリングを繰り返す
    𝜃: モデルパラメータ (𝛽, 𝜎𝑦
    2, 𝜎𝜇
    2, 𝜎𝛿
    2, ...)
    𝛼: 状態ベクトル
    (a) Data augmentation
    𝑝 𝛼 𝑦1:𝑛
    , 𝜃) から 𝛼 を simulation smoother でサンプリング（Durbin et al., 2002）
    (b) Parameter simulation
    共変量 𝛽 と 𝜎𝑦
    2 以外のパラメータ 𝜃′ を 𝑝 𝜃′ 𝑦1:𝑛
    , 𝛼) からサンプリング
    (状態 𝛼 の分散 𝜎2 は事前分布の共役性から計算が容易)
    その後，各共変量 𝛽𝑗
    に関する 各 𝜚𝑗
    および 𝜎𝑦
    2 を Gibbs サンプリングで得る
    1. Posterior Simulation
    2. Posterior Predictive Simulation (Bayesian Model Averaging)
    介入後の反実仮想の観察量 𝑝 𝑦𝑛+1:𝑚
    | 𝑦1:𝑛
    , 𝑥1:𝑛
    を上記の 𝜃, 𝛼 のサンプルから計算し，それらの平均値をとる
    （介入後の全時点を同時分布としてサンプリングしている）
    異なった「共変量の組み合わせ」や「パラメータ」を平均化している
    => Bayesian Model Averaging による sparsity や uncertainty の推論への反映
    3. Evaluating Impact
    最後に各サンプリングと各時点における Causal Impact 𝜙
    𝑡
    (𝜏) を 𝜙
    𝑡
    (𝜏) = 𝑦𝑡
    − 𝑦
    𝑡
    (𝜏) で求める
    

    View Slide

13. "bsts" backend
    CausalImpact
    impact_analysis.R in "CausalImpact"
    RunWithData RunWithBstsModel
    https://github.com/cran/bsts/blob/master/src/bsts.cc
    with data arg?
    yes no
    ConstructModel
    Main components of "CausalImpact" package
    impact_model.R
    in "CausalImpact"
    CompilePosteriorInferences
    impact_inferece.R
    in "CausalImpact"
    Boom::SdPrior
    bsts::AddLocalLevel
    bsts::AddSeasonal
    bsts::bsts
    bsts::AddDynamicRegression
    Dynamic regression?
    yes
    no
    ※ The details of this flow are
    not shown here.
    ComputeResponseTrajectories
    GetPosteriorStateSamples
    ComputePointPredictions
    ComputeCumulativePredictions
    "bsts"
    functions
    

    View Slide

14. 03
    Application to simulated data
    𝑦𝑡
    = 𝛽𝑡, 1
    𝑧𝑡, 1
    + 𝛽𝑡, 2
    𝑧𝑡, 2
    + 𝜇𝑡
    + 𝜀𝑡
    𝛽𝑡
    ~ 𝑁 𝛽𝑡−1
    , 0.012
    𝛽0
    = 1
    𝑧𝑡, 1 or 2
    : 90 or 360 days の周期の sin 波
    𝜇𝑡
    ~ 𝑁 𝜇𝑡−1
    , 0.12
    𝜇0
    = 0
    𝜀𝑡
    ~ 𝑁 0, 0.12
    擬似データで Causal Impact の推定結果を検証
    上記の時系列を生成し，介入後 (2014 - 01 以降) は，
    介入効果として， 𝑦𝑡
    に 1 + 𝑒 をかけ， 𝑒 を Causal Impact で推定
     効果の大きさ 0%, 0.1%, 1%, 10%, 25%, 50%, 100% それぞれに対して，256 回の simulation を行った
     効果の検知の有無は，impact の事後分布の信用区間が zero を含むか含まないかで判定
     効果が 25% 以上だと，概ね 80% 以上の確率で検知可能であった (Fig 3 - b)
     介入期間の長さを変化させても，効果の 95% 信用区間に真の効果の大きさ (10%) が含まれる割合は殆ど変化しなかった (Fig 3 - c)
     但し，擬似データの生成モデルが介入後の途中で変化するような場合 (例えば，動的係数 𝛽 のノイズ成分の大きさが変化する)，推定精度は大きく悪化する (論
    文中の Fig. 4)
     良い感じに推定できているが，使用したモデルは生成モデルと同じ表現式を使用しているので，少し leaky な評価か？
    (要は，真のモデルを知った上で推定しているが，それが分からない時でも同じような推定ができるのか？)
    Brodersen et al., 2015. Fig 3. Adequacy of posterior uncertainty
    

    View Slide

15. 04
    Application to empirical data
    実データで Causal Impact の推定結果を検証
    US における広告効果のデータ (Vaver and Koehler, 2012) を使って検証
     広告は，190 地域のうち 95 地域をランダムに選択し，6 週間実施
     介入効果は広告の click 数で測定
     広告介入した 95 地域の売り上げを合計して，一つの介入地域を生成
     10,000 MCMC samples
     非介入 95 地域を静的係数とした local level model を採用
    (季節性は非介入地域の時系列に含まれており，それらを共変量として使用しているため，組み込まず)
    (a) 広告介入前をよく予測できており，介入後では，counterfactual な推
    定値と観測値との乖離がおおきくなっていき，明らかに広告効果があるこ
    とがみてとれる
    (b) 広告の効果を point-wise で表しており，3 週間付近でピークを向か
    えた後，段々と効果が薄れている
    (c) 広告の累積効果は 95% BCI (bayesian credible interval) でみても有
    意であり，treatment-control comparison による推定累積効果 (Vaver
    and Koehler, 2012) と比較してもほぼ差がない
    Causal Impact による推定累積効果: 22% [13%, 30%], 88,400 clicks
    Vaver and Koehler による推定累積効果: 21% [19%, 22%], 84,700 clicks
    Brodersen et al., 2015.
    ← Model fit Prediction →
    

    View Slide

16. 広告の対象となっていない地域 (controls) における効果の推定結果
     もし，広告の効果が介入された地域にしかないのであれば，controls における効果
    はないはず
     広告元の業界に関連するキーワード検索 (Google Trends) を共変量として使用
    有意な効果は一切みられず
    推定累積効果: 2% [-6%, 10%]
    

    View Slide

17. 1. 『効果検証入門 正しい比較のための因果推論／計量経済学の基礎』
    安井翔太 著，技術評論社，2020
    2. 『時系列分析と状態空間モデルの基礎』馬場真哉 著，プレアデス出版，2018
    3. 『経済・ファイナンスデータの計量時系列分析』沖本竜義 著，朝倉書店，2010
    4. George et al., 1997. Approaches for Bayesian variable selection. Statist.Sinica, 7, 339–374.
    5. Durbin et al., 2002. A simple and efficient simulation smoother for state space time series analysis.
    Biometrika. 89:603–16.
    6. Abadie et al., 2010. Synthetic control methods for comparative case studies: Estimating the effect of California’s tobacco control program.
    J. Amer. Statist. Assoc. 105 493–505.
    7. Scott et al., 2014. Predicting the present with Bayesian structural time series. International Journal of Mathematical Modeling and
    Optimization 5 4–23.
    8. Nakajima J, West M. Bayesian Analysis of Latent Threshold Dynamic Models. Journal of Business & Economic Statistics 2013;31:151–64.
    https://doi.org/10.1080/07350015.2012.747847.
    9. Vaver J, Koehler J. Periodic Measurement of Advertising Effectiveness Using Multiple-Test-Period Geo Experiments. Google Inc.; 2012.
    References
    

    View Slide