はじめに

合成統制法（Abadie and Gardeazabal 2003）やCausalImpact（Brodersen et al. 2015）などの複数の観測対象の長い時系列データを利用した因果推論は政治学、経済学などのアカデミアの世界だけでなく、ビジネスでも広告などの施策の効果検証の際によく登場する。複数の観測対象の長い時系列のことを政治学の世界では「時系列横断面データ」（time-series cross-sectional data）と呼ぶ（Xu 2023）。また、時系列の長さに関わらず、複数の観測対象の時系列データのことをパネルデータと呼ぶことも多い。

本記事は計量政治学における時系列横断面データ因果推論の最先端の手法に焦点を当てたいため、基礎的な合成統制法やCausalImpactの説明や実装は割愛する。詳細は合成統制法やCausalImpactに関連する書籍と記事を参考してください。

時系列横断面データのメリット

時系列横断面データを利用した因果推論にどんなメリットがあるのかというと、まず分析の解釈可能性と透明性が挙げられる。処置群に処置がなかった場合の結果変数の推移を統制群のみを使って予測を行うというアイデアで、処置効果がどうやって推定されたのかが説明しやすい。また、処置群の状況を統制群でうまく再現できるかも処置開始前のデータを使えば可視化・説明しやすい。

伝統的な合成統制法・CausalImpactの限界

伝統的な合成統制法・CausalImpactの最大の弱みは、単一の処置対象しか容認しないところだと思われる。

ビジネスの現場でよく出てくるTVCMの効果検証の話で説明すると、現実の世界では複数の都道府県で同時のTVCMを出向するのが一般的だが、CausalImpact・合成統制法は単一の都道府県に出稿する状況にしか対応しない。

もちろん、Abadie（2021, p4）の脚注3にも書いているように、出稿している都道府県ごとにモデルを作るという対策もあるが、これは若干多重検定の問題が生じる恐れがありそうなので、個人的にはあまりお勧めしたくない。

そこで、次は複数の処置対象を容認する計量政治学の手法と、それらの手法の仮定について説明する。

2種類の仮説と手法の紹介

因果推論は仮説なしで成立しない。ここではXu（2023）を参考に、厳密な外生性（strict exogeneity）と逐次無視可能性（sequential ignorability）という「時系列横断面データ」を使った因果推論でよく置かれる二つの仮説について説明して、上記の二つの仮説を用いた分析手法も一つずつ紹介する。

ただし、説明のポイントは、厳密な外生性と逐次無視可能性はそれぞれどういう状況で成立しそうなのかを日本語で説明するところになるため、数式を使った導出は手法を提案した論文を参考にしてください。

また、本記事の目的は厳密な外生性と逐次無視可能性とそれぞれの手法の比較なので、手法の詳細に関してはまた別記事で取り上げる予定である。

厳密な外生性（strict exogeneity）と一般化合成統制法（generalized synthetic control）

概念

Xu（2023）の整理によると、厳密な外生性（strict exogeneity）とは

（１）時間と変動する交絡変数がない
（２）過去の結果変数の値は現在の結果変数の値に影響を与えない
（３）過去の結果変数の値と共変量の値は現在の処置変数に影響を与えない

を意味する。この3点の中で、1と2は検証が難しい場合もあるが、3はドメイン知識で検証できる。

具体的には、TVCMやOOH広告の効果検証の場合、どの地域にいつ出稿するかはかなり早い段階で決める必要があり、かつ事後的に変更するのは極めて困難であることを、業界経験で知った。

この場合処置変数は契約形式によって定数になるため、過去の結果変数の値と共変量の値からの影響は一切受けない。よって（３）の仮定を問題なく満たしている。

だが、（３）を満たさない状況もビジネスの現場でよく見られる。

下の事例は私の経験から少し改変したものになるので現実とは大きくかけ離れた記述も入っているかもしれない。

去年A社のデータサイエンスチームと営業が協力して、顧客に新しい商品を無料で試していただくことで効果を測定しようとした。営業側の都合でA社と関係性の良い顧客企業を対象にすることになったが、二つの問題がある。

まず関係性は時間と共に変動するものなので、過去の共変量の値（極端な話をすると前日の商談結果）が処置変数（新しい商品の無料テストの対象にするか）に大きく影響する。

また、もし新しい商品に顧客企業が不満を感じたら適用を直ちに中止することになったので、過去の結果変数の値（KPI）も処置変数（新しい商品の無料テストの対象にするか）に大きく影響する。

このように、アジャイル的に施策効果を検証する際に、（３）が満たされないことが多いため、厳密な外生性を元にする手法を利用する際に注意する必要がある。

手法

厳密な外生性の仮定の元で利用できる手法として、一般化合成統制法（generalized synthetic control）を紹介したい（Xu 2017）。

一般化合成統制法が想定する状況は下記の通りである。

まず観測できないトレンドが複数あり、観測値の結果変数は処置効果、共変量及びトレンドから生成されたと仮定する。具体的には、トレンドが二つある場合、式は下のようになる。

結果変数_{観測値i, 時間t} = 処置効果_{観測値i, 時間t} + 処置_{観測値i, 時間t} + 共変量_{観測値i, 時間t} * \beta + トレンド1 + トレンド1のウエイト_{観測値i} + トレンド2 + トレンド2のウエイト_{観測値i}

モデルは、統制群と処置前の処置群のデータを利用して、トレンドと各観測値が持つ各トレンドのウエイトを推定する。最後に、モデルは処置群の処置がない場合の結果変数の反実仮想予測値を推定して、実際の処置群の結果変数と比較して効果を推定する。

また、トレンドの数は交差検証で決める。

詳細はまた後で説明するが、各処置対象の個別処置効果もわかるので、より緻密な結果分析が可能になる。

逐次無視可能性（sequential ignorability）とパネルマッチ（PanelMatch）

概念

ここも同じくXu（2023）の整理からの引用になるが、逐次無視可能性（sequential ignorability）の仮定の元で、厳密な外生性で許容しない現象が許容される。たとえば過去の結果変数と共変量が現在の処置変数に影響しても良いし、過去の結果変数が現在の結果変数に影響する状況も許容する。さらに、過去の処置変数の値が現在の共変量と結果変数に影響する状況も問題なく対応できるので、利用できるシーンが厳密な外生性より多いことがわかる。

上で紹介したデータサイエンティストと営業の事例で説明すると、「関係性の良い顧客企業を対象」（過去の共変量が処置変数に影響）と「不満が出たら即中止」（過去の結果変数が処置変数に影響）の状況にも対応するするため、逐次無視可能性は柔軟性がかなり高いといえる。

だが、逐次無視可能性も万能なものではなく、たとえばImai, Kim and Wang（2021）が指摘するように、未観測な交絡変数がある際に成立しない。また、下で紹介するパネルマッチでいうと、分析者が指定したラグ数より前の処置変数・結果変数・共変量が現在の処置と現在の将来の結果変数の交絡因子になる場合も成立しない。

よって、逐次無視可能性は厳密な外生性より弱い仮説だが、利用する際にはやはり注意が必要である。

手法

逐次無視可能性の仮定を利用した手法として、Imai, Kim and Wang（2021）が提案したパネルマッチ（PanelMatch）を紹介する。

大前提として、この手法は処置と統制の任意な変動（今日は施策あり、明日は施策なし、明後日は施策あり）を許容するので、下記の説明では観測対象自体ではなく、観測対象 x 時間を一つの観測値とする。たとえば処置対象というのは、ある時間において処置の対象になっている観測値のことで、統制対象というのはある時間において処置の対象になっていない観測値のことである。また、たとえば処置対象を削除するというのは、その地点のその観測値のデータを消すことだけであって、その観測値の過去と未来のデータも消すことを意味しないことに留意していただきたい。

パネルマッチは横断面データの因果推論でよく使われるマッチングの一種で、何かしらの基準で処置対象と似ている統制対象を選ばないといけない。パネルマッチでは、分析者が指定したラグ数までの処置の経歴が処置対象と同一の統制対象を選び、マッチ集合（matched set）を作成する。ただし、処置の経歴が同一の統制対象のない処置対象は分析から除外される。詳細はImai, Kim and Wang（2021）のFigure 2を参考にした方がわかりやすいかもしれない。

マッチ集合が作成された後は、傾向スコアマッチングやマハラノビス距離などを使って、他の共変量を調整する。

最後に、調整ずみのマッチ集合を使ってDiDで因果効果を推定する。

だが、Xu（2023）が指摘するように、パネルマッチは多くのデータを捨てることになるだけでなく、パネルデータの時系列構造を無視するなどの問題点がある

データ分析

ここではXu（2017）が利用したデータで、選挙日登録制度（Election Day Registration。以下ではEDRと略す）の導入は投票率の向上につながるかを分析する。

まずはデータの読み込みと処置の状況の可視化から始める。

library(tidyverse)
library(gsynth)
library(panelView)
library(PanelMatch)

data(gsynth)

panelview(turnout ~ policy_edr, data = turnout, 
          index = c("abb","year"), pre.post = TRUE, 
          by.timing = TRUE)

図の凡例を説明すると

Controls：処置の対象になったことのない観測値
Treated（Pre）：処置状態になっていないがどこかのタイミングで処置を受ける観測値
Treated（Post）：処置の状態になっている観測値

この図を確認すればわかるように、処置の対象になる州が複数あるため、CausalImpactで分析しようとすると、処置の対象になった9つの州についてそれぞれモデルを構築しないといけない。

次に一般化合成統制法とパネルマッチを使って分析して、結果を比較する。

一般化合成統制法

ここでは一般化合成統制法にパッケージgsynthの利用方法を説明するサイトに従ってモデル推定を行う。

m_gsynth <- gsynth(turnout ~ policy_edr + policy_mail_in + policy_motor, 
                   data = turnout,  index = c("abb","year"), 
                   se = TRUE, inference = "parametric", 
                   r = c(0, 5), CV = TRUE, force = "two-way", 
                   nboots = 1000, seed = 12345)
Parallel computing ...
Cross-validating ... 
 r = 0; sigma2 = 77.99366; IC = 4.82744; PC = 72.60594; MSPE = 22.13889
 r = 1; sigma2 = 13.65239; IC = 3.52566; PC = 21.67581; MSPE = 12.03686
 r = 2; sigma2 = 8.56312; IC = 3.48518; PC = 19.23841; MSPE = 10.31254*
 r = 3; sigma2 = 6.40268; IC = 3.60547; PC = 18.61783; MSPE = 11.48390
 r = 4; sigma2 = 4.74273; IC = 3.70145; PC = 16.93707; MSPE = 16.28613
 r = 5; sigma2 = 4.02488; IC = 3.91847; PC = 17.05226; MSPE = 15.78683

 r* = 2

Bootstrapping ... ...

交差検証の結果、トレンド数は2になった。

推定の結果は下記の通りで、処置群における平均処置効果の平均（ATT.avg）は4.896である。

> m_gsynth
Call:
gsynth.formula(formula = turnout ~ policy_edr + policy_mail_in + 
    policy_motor, data = turnout, index = c("abb", "year"), force = "two-way", 
    r = c(0, 5), CV = TRUE, se = TRUE, nboots = 1000, inference = "parametric", 
    seed = 12345)

Average Treatment Effect on the Treated:
        Estimate S.E. CI.lower CI.upper p.value
ATT.avg    4.896 2.21   0.5649    9.227 0.02672

   ~ by Period (including Pre-treatment Periods):
         ATT   S.E. CI.lower CI.upper p.value n.Treated
-13 -0.66711 0.9867 -2.60109    1.267 0.49899         0
-12 -0.11117 1.1043 -2.27558    2.053 0.91981         0
-11 -0.11003 1.0365 -2.14156    1.921 0.91546         0
-10  0.46223 0.9863 -1.47090    2.395 0.63932         0
-9  -0.49238 0.9741 -2.40159    1.417 0.61323         0
-8   0.73385 0.8296 -0.89216    2.360 0.37639         0
-7   0.30492 0.9975 -1.65016    2.260 0.75984         0
-6  -0.55177 0.9084 -2.33211    1.229 0.54356         0
-5  -0.40395 0.8189 -2.00903    1.201 0.62183         0
-4  -1.19117 0.7486 -2.65839    0.276 0.11156         0
-3   1.39420 0.7419 -0.05983    2.848 0.06020         0
-2  -0.26077 0.6881 -1.60950    1.088 0.70473         0
-1   0.03661 0.8159 -1.56250    1.636 0.96421         0
0    0.37851 0.9941 -1.56984    2.327 0.70338         0
1    2.59173 1.2310  0.17902    5.004 0.03526         9
2    3.36681 1.6694  0.09483    6.639 0.04372         8
3    3.66693 2.0089 -0.27047    7.604 0.06795         6
4    3.22794 2.4412 -1.55672    8.013 0.18608         6
5    4.95432 2.7287 -0.39377   10.302 0.06942         6
6    5.31104 3.8039 -2.14456   12.767 0.16266         3
7    8.75843 4.1147  0.69375   16.823 0.03329         3
8   10.30966 4.3940  1.69755   18.922 0.01896         3
9    7.40992 5.0847 -2.55582   17.376 0.14503         3
10   9.35570 5.4324 -1.29161   20.003 0.08503         3

Coefficients for the Covariates:
                  beta   S.E. CI.lower CI.upper p.value
policy_mail_in  0.1545 0.8157   -1.444   1.7534  0.8498
policy_motor   -1.0515 0.7768   -2.574   0.4709  0.1758

図で可視化すると下記のようになる。

plot(m_gsynth)

この図は処置群の実際の結果変数の値とモデルが予測した処置がない場合の処置群の結果変数の反実仮想推定量の差で、処置前に差があまりなく、処置後に差が出たことが視覚できに確認できる。CausalImpactを利用しことのある方にとっては馴染みのある図かもしれない。

また、全体だけでなく、個々の州の状況も可視化できる。ここではコネチカット州（abb = CT）とウィスコンシン州（abb = WI）で説明する。

plot(m_gsynth, type = "gap", id = "CT")

plot(m_gsynth, type = "gap", id = "WI")

まず、コネチカット州の場合、処置後でも実際の結果変数と反実仮想推定量に差がなく、EDRの効果があまりないといえる。ウィスコンシン州の場合は処置後に両者に差が出たので、効果があるといえる。

ここでもわかるように、CausalImpactで同じようなことをしようと思ったら、全ての処置対象の州についてモデルを作らないといけないが、一般化合成統制法の場合は一気に全ての処置対象の州の効果を推定できる。

また、モデルが推定したトレンドも可視化できる。

plot(m_gsynth, type = "factors")

Xu（2017）はアメリカの南北の差でトレンドを説明する一方、推定されたトレンドが必ずしも解釈可能にわけではないと指摘する。ただ、私の経験からすると、ドメイン知識を使ってトレンドの意味を説明できた方が、ステイクホルダーの分析結果の対する納得感が上がる。

パネルマッチ

次に、パネルマッチで同様の分析をパネルマッチのパッケージのGithubのサイト

を参考に行う。

まず、パネルマッチのパッケージを利用する際は州の略をIDに変換しないといけないので、まず州と年のマスターテーブルを用意する。

abb_master <- tibble(
  abb = sort(unique(turnout$abb)),
  abb_id = 1:length(unique(turnout$abb))
)

year_master <- tibble(
  year = sort(unique(turnout$year)),
  year_id = 1:length(unique(turnout$year))
)

次に上記の二つのマスター表をデータにJOINするが、パネルマッチを利用する際はデータフレイムに推定で使わないカラムを入れない方がいいので、必要なカラムのSELECTも行う。

turnout_for_pm <- turnout %>%
  left_join(abb_master, by = "abb") %>%
  left_join(year_master, by = "year") %>%
  select(year_id, abb_id, turnout, policy_edr, policy_mail_in, policy_motor)

次に、マッチ集合を作成する。

m_panelmatch <- PanelMatch(lag = 5, time.id = "year_id", unit.id = "abb_id", 
                        treatment = "policy_edr", refinement.method = "mahalanobis", 
                        data = turnout_for_pm, match.missing = TRUE, 
                        covs.formula = ~ I(lag(policy_mail_in, 0:4)) + I(lag(policy_motor, 0:4)) + I(lag(turnout, 1:4)), 
                        size.match = 5, qoi = "att" ,outcome.var = "turnout",
                        lead = 0:5, forbid.treatment.reversal = FALSE)

ここで一般化合成統制法とは異なり、共変量の過去の値と結果変数の過去の値も共変量としてモデルに入れている。

ではマッチ集合で効果量の推定を実際に行う。

m_panelmatch_estimate <- PanelEstimate(sets = m_panelmatch, data = turnout_for_pm)

結果は次のとおりである。

> summary(m_panelmatch_estimate)
Weighted Difference-in-Differences with Mahalanobis Distance
Matches created with 5 lags

Standard errors computed with 1000 Weighted bootstrap samples

Estimate of Average Treatment Effect on the Treated (ATT) by Period:
$summary
    estimate std.error       2.5%    97.5%
t+0 5.468665  1.673681  4.0221957 10.37221
t+1 7.049888  3.136944  4.4837569 16.60847
t+2 5.499993  3.409878  2.8840032 16.35548
t+3 5.589277  5.716803  1.4566003 23.85922
t+4 8.649934  3.367160  3.6111715 17.29016
t+5 6.329012  8.727903 -0.1090847 33.51898

$lag
[1] 5

$iterations
[1] 1000

$qoi
[1] "att"

処置群における平均処置効果（ATT）は処置開始からの時間によって異なるが、だいたい5から8までの範囲で推移する。また、処置開始から5年が経ったら効果は統計的に有意でなくなる。

分析結果を図で可視化すると下記のようになる。

plot(m_panelmatch_estimate)

比較

上記の結果からわかるように、一般化合成統制法でもパネルマッチでも、EDRは投票率の向上に統計的に有意に寄与する結論になるため、今回の事例に関しては手法・仮定によって結論が大きく異なることは見られない。

その中で、一般化合成統制法は個々の観測値の個別処置効果も出せるので、ビジネス上ネクストアクションを検討する際に有用な材料になりうる一方、厳密な外生性という厳しい仮定を満たさないといけないため、利用が正当化されるシーンが限られる。ただ、前述したように、TVCMやOOH広告の分析など高額かつ重要な施策の効果検証で使えるので、ビジネスの場で活用できる場面が多いと思われる。

パネルマッチの場合、より弱い逐次無視可能性を満たせば利用できるので、アジャイル的に施策を検証する際に利用でき、ビジネスの場面で活用できる場面は一般化合成統制法より多いが、一般化合成統制法のように個別処置効果を推定するのは原理的に難しいので、手法自体で得られる示唆が限られる可能性がある。

参考文献

Abadie, Alberto. "Using synthetic controls: Feasibility, data requirements, and methodological aspects." Journal of Economic Literature 59.2 (2021): 391-425.

Abadie, Alberto, and Javier Gardeazabal. "The economic costs of conflict: A case study of the Basque Country." American economic review 93.1 (2003): 113-132.

Brodersen, Kay H., Fabian Gallusser, Jim Koehler, Nicolas Remy and Steven L. Scott
. "Inferring causal impact using Bayesian structural time-series models." The Annals of Applied Statistics (2015): 247-274.

Imai, Kosuke, In Song Kim, and Erik H. Wang. "Matching methods for causal inference with time‐series cross‐sectional data." American Journal of Political Science (2021).

Xu, Yiqing. "Generalized synthetic control method: Causal inference with interactive fixed effects models." Political Analysis 25.1 (2017): 57-76.

Xu, Yiqing, Causal Inference with Time-Series Cross-Sectional Data: A Reflection (April 2, 2023). Available at SSRN: https://ssrn.com/abstract=3979613 or http://dx.doi.org/10.2139/ssrn.3979613