第4章 差分の差分法（DID）とCausalImpact：ユーザ単位で割り当てられない施策の効果を測る

4.1.1 DIDが必要になる状況

回帰分析と傾向スコアは、介入グループと非介入グループの両方に、同じような特徴を持つサンプルが含まれている状況を利用して比較する方法であった。しかし、実際に得られるデータセットには、介入・非介入それぞれのグループに同質なサンプルが存在しないことがよくある。これは、 介入が広い範囲に一斉に影響する 場合によく発生する。特定の地域における政策や法律の変更のほか、特定の地域で広告を出稿する場合や、価格を変更する場合などが典型である。

たとえば、ある都市で条例を変更した場合、その効果は都市内のすべての人に及ぶ。ある地域のすべてのスーパーで商品を10%値引きしたとき、その地域の中に値引きしない店舗は存在しない。効果を測る観点では、介入を個人や店舗の単位でランダムに割り当てるRCTが最適である。第2章・第3章の回帰分析や傾向スコアも、介入が個人や店舗の単位で選ばれているデータを前提にしている。しかし、介入が地域単位でしか決められない場合、これらの分析方法ではうまく対応できない。

介入を行った地域の中に非介入グループがないなら、他の地域のデータを非介入グループとして扱うことを考える。しかし後述するように、この方法では地域ごとの違いによるセレクションバイアスが生じる。他の地域がだめなら過去との比較を考えるが、この場合も「過去と今の差が介入だけである」という保証がない。

DIDは、介入が行われた地域の介入前後のデータに加えて、 他の地域の介入前後のデータ を使うことで、これらの問題を乗り越える方法である。大まかには、介入前後の差分を介入グループと非介入グループでそれぞれ計算し、さらにその差分同士の差をとる。2回の差分をとる手続きが、名前の由来である。

興味のある介入が特定の共変量、とくに地理的な情報と完全に相関するケースは、実店舗・マス広告・行政のように、管理運営が地理的な区画ごとに行われる場合によく見られる。経済学では、ある都市の政策の効果や、テロなどの事象が経済に与える影響の評価に使われてきた。Card (1990) は1980年にボートでマイアミに到着した移民が労働市場に与えた影響を分析している。Blake et al. (2015) は、eBayにおいてGoogleでの検索連動型広告の出稿を一部地域で取りやめる疑似的な実験を行い、DIDで分析した。その結果、eBayにとって検索連動型広告の売上への効果は限定的であることを示している。

Huluの推薦の文脈に置き換えると、「地域」にあたるのはデバイス種別、アプリのバージョン、契約経路（キャリア提携など）、そして「作品」である。たとえば新しい推薦モデルをAndroid TVのアプリだけに先行導入した場合、介入はデバイスの単位で完全に決まっている。同じデバイスの中に非介入ユーザはいない。

4.1.2 集計による効果検証とその欠点

地域ごとに複数の時期のデータが手に入った場合、二つの単純な集計による効果検証がよく行われる。

一つ目は、介入を受けた地域と受けなかった地域を比較する方法 である。しかしこの場合、地域と介入が完全に相関してしまう。地域A、B、C、Dがあり、Aだけでスーパーの価格変更を行った状況を考える。AとB・C・Dの平均を比べると、その差には本来の介入効果のほかに、地域固有の効果も含まれる。介入する地域を実験で選んでいても、何らかの意思決定で選んでいても、この問題は残る。データが店舗ごとに手に入ったとしても同じである。

二つ目は、同じ地域の値下げ前と値下げ後を比較する方法 で、いわゆる 前後比較 である（図4.1）。前後比較は同一の地域で比べているので、地域によるセレクションバイアスはなくなる。一見正しそうに見える。しかし、セレクションバイアスはサンプル間だけでなく、 同一サンプルのタイミング間 でも発生する。価格変更のあとに売上が変化したとしても、価格変更がなくとも同じぶんだけ売上が変化する状況であれば、その変化は価格によるものではなく、他の要因で自然に変わっただけである。同一サンプルの前後を比較した場合、その差には本来の効果のほかに、時間を通した自然な変化（トレンド）が含まれてしまう。

図4.1の実線は実際の売上の推移で、点線は「介入が起きなかった場合」の推移である。前後比較は、介入前の水準と介入後の実線の差を効果とみなす。しかし本来の効果は実線と点線の差であり、点線が介入前の水準から上に動いているぶんだけ、前後比較は効果を過大に見積もる。この差が「前後比較でのバイアス」である。

4.1.3 DIDのアイデアを用いた集計分析

同じ時期の他地域との比較にも、同じ地域の前後比較にもバイアスが生まれるなら、どう分析すればよいのか。

イギリスの医師John Snowは、19世紀中ごろにコレラの感染源を探る中でこの課題に直面し、DIDの基礎となるアイデアで分析を行った。そのアイデアは、 介入を受けるグループと受けないグループでそれぞれ前後比較を行い、その結果の差を見る というものである。コレラの感染源が空気にあるのか水源にあるのかを検証していたSnowは、死者数の多い地域で水の供給源が変わったことに着目し、水源の変化を介入として分析を試みた。

調査対象の地域は、Southwark and Vauxhall社とLambeth社という二つの会社によって水が供給されていた。1849年にはどちらの会社もテムズ川を水源としていたが、1852年にLambeth社は汚染が少ないと考えられる上流へ水源を移した。したがって、Southwark and Vauxhall社の供給地域では1849年も1854年も同じ水源が使われている一方、Lambeth社が供給する地域では1854年には別の水源に変わっている。もし感染源が水源にあるなら、Lambeth社が供給する地域では死者数が減るはずである。表4.1はSnow (1855) のTable 12から一部を抜粋したものである。

まずSouthwark and Vauxhall社の地域では、1849年から1854年にかけて死者数が増えている。水源が変わらなかった地域では、コレラが依然として猛威をふるっていたことが分かる。Lambeth社のみの地域は数が少ないので死者数自体は少ないが、162人から37人へと減少している。もし水源が変わらず、Southwark and Vauxhall社の地域と同じように死者が増えていたはずだとすれば、Lambeth社の地域では水源の変化によって死者が減ったことになる。両社が供給していた地域では、一部だけがLambeth社の水なので効果は限定的と想定されるが、それでも3,905人から2,547人へ大幅に減っており、水源の変化に効果があったことを示唆している。

このように、 「非介入グループの変化」と「介入グループが介入を受けなかった場合の変化」が一致するだろう という仮定を 平行トレンド仮定（Common Trend Assumption） と呼ぶ。詳細は4.1.6で扱う。

John Snowのデータでの集計結果

教科書は、Snowの論文からRでデータを作り直して集計している。本ノートでは同じ値をPython（Polars）で集計した。データは会社別・年別の死者数で、両社が供給する地域を介入グループ（ LSV = 1 ）としている。

company	LSV	year_1849	year_1854	gap	gap_rate
Southwark and Vauxhall	0	2,261	2,458	197	0.0871
Lambeth & Southwark and Vauxhall	1	3,904	2,547	−1,357	−0.3476

gap は各地域の1849年と1854年の死者数の差、つまり前後比較の結果である。 gap_rate はその差を1849年の死者数で割った比率である。非介入グループであるSouthwark and Vauxhall社の地域では死者が197名（8.7%）増えている。これが時間の効果であり、介入がなかった場合にどの程度死者が増えたかを表す。介入グループである両社供給地域では、死者が1,357名（34.8%）減っている。DIDによる効果量は、これらの差である。

$$\tau=(-1357)-(197)=-1554$$

比率で計算すれば $-34.76\%-8.71\%=-43.5\%$ である。水源が変わった地域では死者数が大幅に減少したことが確認できる。

図4.2は、この集計をプロットしたものである。(1)は1849年から1854年にかけて非介入グループで起きた死者数の変化の大きさ（+197）である。平行トレンド仮定が成り立つなら、介入グループでも介入がなければ同じぶんだけ変化していたはずで、それが(2)の点線の矢印として描かれている。実際に観測された介入グループの変化は(3)の −1,357 である。よって効果は、(3)の実際の変化と(2)の「あったはずの変化」の差、すなわち $-1357-197=-1554$ となる。

4.1.4 回帰分析を利用したDID

前節までで、DIDは「4つの平均を使った2回の引き算」であることが分かった。引き算だけで効果が出るなら、なぜここで回帰分析が出てくるのか。教科書はいきなり回帰モデルの式から始めるが、その前に「集計のDIDには何が足りないのか」と「第2章で学んだ回帰と、ここでの回帰は同じものなのか」を確認しておく。ここを押さえると、回帰モデルの式が自然に読めるようになる。

uv run output/ch4_examples.py

import numpy as np
import polars as pl


def load_john_snow() -> pl.DataFrame:
    """Snow (1855) Table 12 の地区別死者数。LSV=1 が両社供給地区（介入グループ）。"""
    sv1849 = [283, 157, 192, 249, 259, 226, 352, 97, 111, 8, 235, 92]
    sv1854 = [371, 161, 148, 362, 244, 237, 282, 59, 171, 9, 240, 174]
    lsv1849 = [256, 267, 312, 257, 318, 446, 143, 193, 243, 215, 544, 187, 153, 81, 113, 176]
    lsv1854 = [113, 174, 270, 93, 210, 388, 92, 58, 117, 49, 193, 303, 142, 48, 165, 132]

    def block(deaths: list[int], year: int, prefix: str, lsv: int) -> pl.DataFrame:
        return pl.DataFrame({
            "area": [f"{prefix}_{i + 1}" for i in range(len(deaths))],
            "year": [year] * len(deaths),
            "death": deaths,
            "LSV": [lsv] * len(deaths),
        })

    return pl.concat([
        block(sv1849, 1849, "sv", 0), block(sv1854, 1854, "sv", 0),
        block(lsv1849, 1849, "lsv", 1), block(lsv1854, 1854, "lsv", 1),
    ]).with_columns(D1854=(pl.col("year") == 1854).cast(pl.Int64))


js_df = load_john_snow()

# 会社別（LSV別）・年別に集計し、前後比較の差をとる
js_sum = js_df.group_by("LSV", "year").agg(pl.col("death").sum()).sort("LSV", "year")
js_grp = (
    js_sum.pivot(on="year", index="LSV", values="death")
    .rename({"1849": "year_1849", "1854": "year_1854"})
    .with_columns(
        gap=pl.col("year_1854") - pl.col("year_1849"),
        gap_rate=pl.col("year_1854") / pl.col("year_1849") - 1,
    )
)

def ols(frame: pl.DataFrame, outcome: str, predictors: list[str],
        cluster: str | None = None) -> OLSResult:
    """最小2乗法。cluster を指定するとクラスター標準誤差を返す。"""
    x, names = design_matrix(frame, predictors)   # 文字列列はダミー展開、"a:b" は積
    y = frame.get_column(outcome).cast(pl.Float64).to_numpy()
    n, k = x.shape
    xtx_inv = np.linalg.pinv(x.T @ x)
    beta = xtx_inv @ x.T @ y
    resid = y - x @ beta
    if cluster is None:
        sigma2 = resid @ resid / (n - k)
        var = sigma2 * xtx_inv
    else:
        groups = frame.get_column(cluster).to_numpy()
        meat = np.zeros((k, k))
        for g in np.unique(groups):          # クラスターごとにスコアを合計してから外積
            idx = groups == g
            score = x[idx].T @ resid[idx]
            meat += np.outer(score, score)
        n_groups = len(np.unique(groups))
        correction = (n_groups / (n_groups - 1)) * ((n - 1) / (n - k))
        var = correction * xtx_inv @ meat @ xtx_inv
    return OLSResult(names, beta, np.sqrt(np.diag(var)), resid, n, k)

js_sum = js_sum.with_columns(D1854=(pl.col("year") == 1854).cast(pl.Int64))
print(ols(js_sum, "death", ["LSV", "D1854", "LSV:D1854"]).tidy)

show(ols(js_df, "death", ["LSV", "area", "D1854", "LSV:D1854"]), drop_prefix=("area",))

js_log = js_df.with_columns(log_death=pl.col("death").log())
show(ols(js_log, "log_death", ["LSV", "area", "D1854", "LSV:D1854"]), drop_prefix=("area",))

4.1.5 DIDにおける標準誤差

通常の回帰分析では、大まかには一つの観測対象から一つのデータが得られることを想定している。一方でDIDは、同一の対象からいくつかの期間に取得したデータを利用する。このような場合、 自己相関（auto-correlation, serial correlation） を持つデータになる可能性がある。自己相関とは、ある時点で取得された値が、その近辺の時点の同じ変数の値と相関する状態をいう。

たとえば、あるコンビニの売上を毎日観測した場合、店舗近くの人通りや品ぞろえは大きく変わらないので、ある日とその前後の日の売上は非常に似た値になる。この売上データは自己相関を持っている可能性が高い。

自己相関を持つデータで回帰分析を行うと、誤差項の値が同一店舗内で似通ったものになる。これは、 「独立な観測がたくさんある」という回帰の前提よりも、実質的な情報量が少ない ということである。回帰の標準誤差は観測が独立であることを前提に計算されるので、結果として標準誤差が過小に算出される。有意差検定の観点では、統計的に有意な結果が過剰に得られてしまう。回帰分析が想定していないデータ構造に対して、強引に回帰を使ったために起きる問題である。

4.1.6 平行トレンド仮定（Common Trend Assumption）と共変量

本章のこれまでの分析では、「もしLambeth社が水源を変更しなかったなら、両社供給地域でもSouthwark and Vauxhall社の地域と同じように死者が増えていたはずだ」という仮定を置いている。

介入グループと非介入グループの目的変数の時間を通じた変化、いわゆるトレンドが同一であるという仮定を 平行トレンド仮定（Common Trend Assumption） と呼ぶ。これはDIDの分析結果を正しいと考えるために必要な、もっとも重要な仮定である（図4.3）。この仮定が成り立たない場合、推定された効果は本来の効果にトレンドの乖離分を加えたものになる。これは、回帰分析で重要な共変量が脱落した場合に発生するOVB（脱落変数バイアス）と同じ問題である。

図4.3では、水源が変わらない地域の死者数（上の実線）と、水源が変わる地域の死者数（下の実線）が描かれている。介入開始後、介入グループには「介入が行われた場合の推移」（実線、下がっていく）と「介入が行われなかった場合の推移」（破線、実際には観測されない）の二つがある。平行トレンド仮定は、この破線が非介入グループの実線と 同じ動き をすると仮定するものである。効果は、破線と介入グループの実線の差として測られる。

4.2.1 DIDの欠点

集計済みのデータしか持てない分析者にとって、DIDは最後の頼みの綱とも言える手法だが、二つの欠点がある。

一つは、効果を調べたい変数が複数の場所や時期で得られている必要がある ことである。スーパーの売上への価格変更の効果を考える際、売上データは介入地域だけでなく、介入のなかった別の地域でも手に入れる必要がある。しかし介入対象のデータしかないことはよくある。たとえばすべての地域に広告を配信した場合、手に入るのは広告の介入を受けた売上データだけである。Huluで全ユーザに訴求を行った場合も同じで、「訴求を受けなかったユーザの視聴数」は存在しない。

もう一つは、どのデータを分析に用いるかが分析者の仮説に依存する ことである。DIDは平行トレンド仮定を必要とする。つまり分析者は、この仮定が満たされるように非介入グループを構成するか、共変量を調整する必要がある。しかし4.1.6で見たように、仮定が実際に満たされているかは分からないので、近隣地域を使ったり、知り得る情報から共変量を選んだりすることになる。

4.2.2 CausalImpactのアイデア

CausalImpactは、複数の期間のデータを必要とする一方で、上の二つの欠点を補える分析方法である。

DIDのアイデアを言い直すと、 介入が行われたサンプルの「介入が行われなかった場合」の結果を、非介入グループのデータで補う というものであった。この見方を一歩進めると、本来この反実仮想の結果を何らかの方法で予測できるなら、 どのようなデータから予測しても構わない ということになる。もしうまく予測できるなら、実際に観測された値と予測値の差が介入の効果である。

CausalImpactは、さまざまな種類の変数 $X$ を利用して目的変数 $Y$ を予測するモデルを、 介入が行われる前の期間だけで 作る（図4.4）。図4.4では、商品Aの売上を商品Bの売上から予測するモデルを介入前に学習し、介入後は商品Bのデータをそのモデルに入れて「商品Aに介入がなかった場合の売上」を予測している。実線（実測）と破線（予測）の差が介入の効果である。

たとえば、ある店舗で価格を変えた効果を推定する際には、 $Y$ として売上を用いる一方で、 $X$ として価格を変えた商品名の検索回数などを利用できる。検索回数は、その商品に興味のあるユーザ数や認知度に応じて変動するので、商品の潜在的な売上と強い相関を持つ。したがって、価格を変えなかった場合の売上とも強い相関を持つ。 $X$ の種類には限りがないので、複数のキーワードの検索回数やWebページのアクセス回数といった情報も使える。

さらにCausalImpactは、介入前のデータを利用して、どの変数が $Y$ の予測に役立つかを判別し、使うデータを自動的に決めながらモデルを学習してくれる（内部ではスパイク・スラブ事前分布という仕組みで、多数の候補から少数の変数を選ぶ）。学習されたモデルは、 $X$ を入力すると介入前の状態の $Y$ を予測する。よって介入後の $X$ を入力すれば、介入が行われなかった場合の $Y$ の予測が出力され、この予測値と実際の $Y$ の差が効果として得られる。

4.3.1 データの準備

教科書はRの Ecdat パッケージの Cigar データセットを使う。本ノートでは同じデータをRdatasetsの公開CSVから取得し、 data/Cigar.csv に保存して使う。分析に使うのは全米各州における人口1人あたりのタバコの売上（箱数）である。カリフォルニア州（ state == 5 ）の売上を $Y_t$ とし、それ以外の州のデータを $X_{it}$ とする。本来の分析は1970〜2000年のデータを使うが、入手できるデータの関係から期間は1970〜1992年である（教科書では1993年までと書かれているが、 Cigar データは1992年までしかない）。

また本来の分析と同じく、1988年以降にタバコへ一定水準以上の増税を行った州は分析から除外する。このような州を含めると、カリフォルニア州以外でも別の介入が行われている状態になり、「介入がなかった場合のカリフォルニア州」との平行トレンド仮定が満たされなくなる。

SKIP_STATE = [3, 9, 10, 22, 21, 23, 31, 33, 48]  # 1988年以降に大幅増税した州など
CA = 5

cigar = (
    pl.read_csv("data/Cigar.csv")
    .filter(~pl.col("state").is_in(SKIP_STATE), pl.col("year") >= 70)
    .with_columns(
        area=pl.when(pl.col("state") == CA).then(pl.lit("CA")).otherwise(pl.lit("Rest of US")),
        post=(pl.col("year") > 87).cast(pl.Int64),
        ca=(pl.col("state") == CA).cast(pl.Int64),
        year_dummy=pl.format("D_{}", pl.col("year")),
    )
)


def weighted_sales(frame: pl.DataFrame, by: list[str]) -> pl.DataFrame:
    """人口（16歳以上）で重み付けした1人あたり売上。R の sum(sales*pop16)/sum(pop16)。"""
    return (
        frame.group_by(by)
        .agg(sales=(pl.col("sales") * pl.col("pop16")).sum() / pl.col("pop16").sum())
        .sort(by)
    )

「Rest of US」の売上は、州ごとの1人あたり売上を単純平均するのではなく、16歳以上人口で重み付けして平均している。人口の多い州の売上が全体をより強く代表するようにするためである。

4.3.2 DIDの実装

まず前後比較をする。介入が始まった1988年より前か後かで、1人あたり売上を比較する。

before_after = weighted_sales(
    cigar.with_columns(period=pl.when(pl.col("post") == 1).then(pl.lit("after")).otherwise(pl.lit("before"))),
    ["area", "period"],
).pivot(on="period", index="area", values="sales").with_columns(change=pl.col("after") - pl.col("before"))

図4.5はこの表をプロットしたものである。カリフォルニア州では40箱程度売上が減少しているが、他の州でも20箱程度減少している。このことから、大まかには $-39.5-(-20.0)=-19.6$ 箱程度の減少がProposition 99の影響だと考えられる。

did_sum = weighted_sales(cigar, ["post", "year", "year_dummy", "ca"])
show(ols(did_sum, "sales", ["ca", "post", "ca:post", "year_dummy"]), drop_prefix=("year_dummy",))

did_sum_log = did_sum.with_columns(log_sales=pl.col("sales").log())
show(ols(did_sum_log, "log_sales", ["ca", "post", "ca:post", "year_dummy"]), drop_prefix=("year_dummy",))

did_state = weighted_sales(cigar, ["post", "ca", "year", "year_dummy", "state"]).with_columns(
    state=pl.col("state").cast(pl.String)
)
plain = ols(did_state, "sales", ["ca", "state", "post", "ca:post", "year_dummy"])
clustered = ols(did_state, "sales", ["ca", "state", "post", "ca:post", "year_dummy"], cluster="state")

平行トレンド仮定は満たされているか

しかし、カリフォルニア州以外のデータは平行トレンド仮定を満たしているだろうか。これを確認するために、年ごとの売上をプロットする（図4.6）。

介入前のトレンドは、カリフォルニア州とそれ以外で特に平行に動いているわけではない。1970年代前半は両者がほぼ同じ水準だったが、1970年代後半から他州の売上が伸びる一方でカリフォルニア州は横ばいから減少に転じ、差が開いていっている。数値で確かめると、両者の差（CA − Rest of US）は次のように推移している。

年	1970	1975	1980	1985	1987	1988	1990	1992
CA − Rest	+8.0	−1.5	−14.2	−19.7	−19.0	−23.3	−27.5	−34.9

介入前の18年間で差が +8.0 から −19.0 まで一方向に広がっている。つまり、 介入がなくても差は広がり続けていた可能性が高い 。介入後の差の拡大の一部は、この既存のトレンドの延長かもしれない。DIDの −20.5 という推定値は、その分だけ効果を過大に見積もっている疑いがある。この確認の作業（イベントスタディと呼ばれる）は、追加手法の節で改めて扱う。

4.3.3 CausalImpactの実装

CausalImpactを試す。CausalImpactには、目的変数 $Y$ 、予測に使う変数 $X$ 、介入期間を表す情報を入力する。このとき $X$ は、これまで登場したような共変量ではなく、平行トレンド仮定に相当する条件を満たす（介入の影響を受けず、 $Y$ をよく予測する）変数である。目的変数はカリフォルニア州の売上で、予測変数には他の州の売上を使う。

教科書ではRの CausalImpact パッケージを使う。Pythonでは、Googleの実装をTensorFlow Probabilityに移植した tfcausalimpact パッケージが同等の機能を持つ。ただしTensorFlowへの依存が重いため、本ノートでは同じ構造のモデル（ローカルレベル + 回帰）を statsmodels の状態空間モデルで書き、考え方を再現する。

from statsmodels.tsa.statespace.structural import UnobservedComponents

wide = cigar.select("state", "year", "sales").pivot(on="state", index="year", values="sales").sort("year")
years = wide.get_column("year").to_numpy()
y = wide.get_column(str(CA)).to_numpy().astype(float)           # カリフォルニア州の売上
others = [c for c in wide.columns if c not in ("year", str(CA))]
x_all = wide.select(others).to_numpy().astype(float)              # 他の州の売上（列 = 州）
pre, post = years < 88, years >= 88

# R 版のスパイク・スラブによる変数選択の代わりに、介入前の相関が高い 3 州を使う
corr = np.array([np.corrcoef(y[pre], x_all[pre, j])[0, 1] for j in range(x_all.shape[1])])
top = np.argsort(-corr)[:3]
x = x_all[:, top]

# 介入前だけで「ローカルレベル + 回帰」モデルを学習し、介入後の X を入れて反実仮想を予測する
fit = UnobservedComponents(y[pre], level="llevel", exog=x[pre]).fit(disp=False)
forecast = fit.get_forecast(steps=int(post.sum()), exog=x[post])
pred = forecast.predicted_mean
effect = y[post] - pred                     # pointwise effect
cumulative = np.cumsum(effect)              # cumulative effect

結果は次のようになった。

教科書のRによる CausalImpact の出力は、平均の効果が −14（標準偏差7.5、95%信用区間 [−36, −0.79]）、累積が −70、比率で −15% であった。本ノートの簡易版でも平均 −14.3、累積 −71、比率 −15.6% と、点推定はほぼ一致している。一方で区間はRの出力よりかなり狭い。これは、本ノートの実装が最尤推定でパラメータを一点に決めてしまい、パラメータの不確かさを無視しているためである。R版はベイズ推定でパラメータの不確かさも含めて事後サンプルから区間を出すので、正直に広い区間になる。 区間の妥当性まで含めて結果を使うなら、 tfcausalimpact のようなベイズ実装を使うべきである。

図4.7はRの plot(impact) の出力で、上から「original」「pointwise」「cumulative」の3つのグラフからなる。

• original は、カリフォルニア州のタバコ売上（実線）と、その予測値（破線）を横軸に時間をとって描いたものである。灰色の帯は予測の信用区間である。介入開始を示す縦の破線の前では実線と破線がほぼ重なっており、介入後に乖離が生まれていることが分かる。
• pointwise は、originalの実測値と予測値の差を各時点で描いたもので、 各時点における介入の効果 を示す。介入前は0付近をうろつき、介入後に徐々に負の方向へ大きくなっている。
• cumulative は、pointwiseの値を介入開始以降で積み上げたものである。Proposition 99は売上を減らすので、下に向かって増えていく。

pointwiseのグラフから分かる通り、CausalImpactの特徴は 各時点で効果を推定できる 点にある。DIDでは介入開始後の平均の効果しか分からない。段階的に強化されるような介入であっても、その効果の推移を見ることができる。

Rの出力では、Absolute effectが効果を実際の箱数で表しており、左側（Average）に平均的な効果、右側（Cumulative）に積み上げの効果が示される。それぞれのカッコ内は標準偏差である。95%の信用区間に0が含まれるか否かを考えることは、有意差検定に類似した情報を与える。教科書の出力では区間 [−36, −0.79] に0を含んでいないため、効果が0である可能性は低いと判断している。比率での効果はRelative effectとして報告され、解釈はAbsolute effectと同じである。

4.3.4 分析結果の比較

この分析の元となったAbadie et al. (2010) は、 合成コントロール（Synthetic Control） という手法を使っている。合成コントロールでは、まず各州の経済や人口の状況が介入前にカリフォルニア州と類似している州をアルゴリズムで抽出する。抽出された州は介入前にカリフォルニア州と似たトレンドを持っているので、平行トレンド仮定を満たしていると考えられる。そして抽出された州のデータを重み付き平均して、「合成されたカリフォルニア州」の売上を予測する。

図4.8はAbadie et al. (2010) からの引用で、CausalImpactのpointwiseに対応するグラフである。縦軸は実際のカリフォルニア州と合成コントロールの差（効果量）、横軸は年である。介入を開始した1988年から徐々に売上が低下し、1992年で15箱程度、その後1995年ごろまで効果が強くなり、およそ20箱程度の減少となっている。介入前の期間で差がほぼ0に張り付いていることが、合成コントロールが介入前のカリフォルニア州をよく再現している証拠である。

これに対してDIDの結果は、1988〜1992年の平均的な効果が20箱程度の減少となっており、やや過大に評価されている。これは4.3.2で確認した通り、平行トレンド仮定がうまく満たされていないデータを使っているためだと考えられる。一方でCausalImpactの結果は、徐々に効果が強くなり1992年で18箱前後の減少（本ノートの再現では20.8箱）となっており、推定された効果の傾向にはある程度の類似が認められる。

from scipy import optimize

def loss(w: np.ndarray) -> float:
    return float(np.sum((y[pre] - x_all[pre] @ w) ** 2))

n_donor = x_all.shape[1]
result = optimize.minimize(
    loss, x0=np.full(n_donor, 1.0 / n_donor), method="SLSQP",
    bounds=[(0.0, 1.0)] * n_donor,
    constraints=[{"type": "eq", "fun": lambda w: np.sum(w) - 1.0}],
)
synth = x_all @ result.x
gap = y - synth

Huluでの同じ構造：介入が「他の腕」に漏れる

推薦システムでも、ユーザ単位のA/Bテストが「本来比較したい状況」を作れないことがある。

• 人気度に基づく特徴量の共有 。「いま人気」トレイのランキングが全ユーザの再生ログから計算されている場合、半分のユーザに新しいロジックを配信すると、その視聴が人気度に流れ込み、旧ロジックのユーザのランキングも変わる。両腕が互いに影響し合うので、「新ロジックを全員に配信した世界」と「旧ロジックを全員に配信した世界」の比較にはならない。
• 在庫や供給側の制約 。TVODの割引施策では、半分のユーザにだけ割引を出すと、割引を受けなかった側の購入も（口コミなどを通じて）変わり得る。
• 世帯内の複数プロフィール 。同じテレビを家族で使う場合、プロフィール単位の割り当ては相互に漏れる。

こうした「干渉」があるときは、介入の単位を 時間 （ある週は新ロジック、次の週は旧ロジックというスイッチバック）や デバイス・地域 に切り替え、DIDやCausalImpactで分析する方が、比較したい世界に近づける。

4.4.1 DIDのアイデアを用いた分析が使えないとき

DIDやCausalImpactでは、効果を分析したい介入が 他の介入や施策と同時に導入される 場合、その効果を分析できない。ある広告の効果を知りたいときに、広告の開始と同時に商品の値下げを行っているような状態を考える。このときDIDやCausalImpactで正の効果が推定されても、その結果は広告の効果と値下げの効果が混ざったものになる。価格が頻繁に変動する場合には、価格を共変量としてモデルに入れることでこの問題を解決できる。実際、このようなデータが手に入ることは珍しくない。

何かの介入を行う際に、その開始のタイミングを揃えることがよくある。マーケティングでは、商品の入れ替え、Webページの刷新、広告が同一のタイミングで始まる。こうしたケースでは、どれか一つの効果を推定することはほぼ不可能である。このような広告と他の施策とのタイミングの重なりが生むバイアスは、マーケティングの分析者の間で アクティビティバイアス（Activity Bias） と呼ばれている。

Huluの推薦開発でこれがもっとも起きやすいのは、 新作の配信開始日 である。配信開始、プッシュ通知、ホームでの訴求、地上波での放送、SNSでの話題化がすべて同じ日に集中する。ここで新しい推薦モデルの投入も同じ日に重ねると、モデルの効果は永久に分離できなくなる。対策は分析手法ではなく 運用 にある。推薦モデルの投入日は配信イベントから意図的にずらす、あるいは施策ごとに開始日を1〜2週間ずらすだけで、DIDやCausalImpactが使えるデータになる。

1. イベントスタディと事前トレンドの確認

平行トレンド仮定を「見る」もっとも標準的な方法が イベントスタディ である。介入前後の各時点について、介入グループと非介入グループの差を、介入直前の時点を基準として推定する。回帰では、年ダミーと介入グループの交互作用項を（基準年を除いて）すべて入れる。

$t^*$ は介入の開始時点で、 $\delta_k$ は介入から $k$ 期後（負なら前）における差の変化である。介入前の $\delta_k$ （ $k<-1$ ）がすべて0付近なら、平行トレンドが介入前に成り立っていたことの証拠になる。介入後の $\delta_k$ は、時点ごとの効果になる。

Proposition 99のデータで、1987年を基準にしたCAと他州の差の変化を計算すると次のようになる（追加手法のためスクリプトに含めてある）。

1983年以降は0付近で安定しているが、1970年代は大きく正であり、介入前を通して平行だったとは言えない。これは、DIDの比較期間を1983年以降に絞る、あるいは合成コントロールを使うといった判断の根拠になる。

事前トレンドの検定には限界もある。Roth (2022) は、「事前トレンドが有意でない」ことは「平行トレンドが成り立つ」ことを意味せず、検定の検出力が低いだけであることが多いと指摘している。Rambachan and Roth (2023) の HonestDiD は、「介入後のトレンドの乖離は、介入前に観測された乖離の $M$ 倍以内」という形の仮定を明示的に置き、その下で効果の区間を出す。平行トレンドを「成り立つか否か」ではなく「どれだけ破れても結論が変わらないか」という感度分析として扱う考え方である。

2. 介入時期がばらばらのとき：Staggered DID

実務では、介入の開始時期がユニットごとに異なることが多い。たとえば新しい推薦モデルを、Fire TVに4月、Android TVに5月、iOSに6月と順次展開する場合である。この場合、すべての交互作用項を一つの $\beta_3$ で表す 二方向固定効果（TWFE）回帰 をそのまま使うと、Goodman-Bacon (2021) が示したように、「すでに介入されたユニット」が「あとから介入されるユニット」の比較相手として使われてしまう。効果が時間とともに変化する場合、この比較は誤った符号すら生む。

対策として、Callaway and Sant’Anna (2021) は、介入開始時期のコホートごと・時点ごとの効果 $ATT(g,t)$ を「まだ介入されていないユニット」だけを比較相手にして推定し、それらを集計する方法を提案した。Sun and Abraham (2021) も同様に、コホートごとの交互作用項を使う推定量を提案している。Pythonでは pyfixest パッケージや differences パッケージがこれらを実装している。 順次展開の効果をDIDで測るなら、単純なTWFE回帰ではなく、これらの推定量を使うべきである。

3. 合成コントロールの発展：Synthetic DID と Augmented SC

4.3.4で試した合成コントロールにはいくつかの弱点がある。介入前の再現が完全でないと効果にバイアスが残ること、統計的な推論が難しいこと、介入ユニットが候補ユニットの「凸包」の外にあると再現できないことである。

• Augmented Synthetic Control（Ben-Michael, Feller, and Rothstein, 2021）は、合成コントロールの介入前の再現誤差を、リッジ回帰などの結果モデルで補正する。介入前の再現が悪いときにバイアスを減らせる。
• Synthetic DID（Arkhangelsky et al., 2021）は、合成コントロールの「ユニット方向の重み」に加えて「時点方向の重み」を導入し、DIDのように水準の差を許容しながら重み付きの差分の差分を計算する。合成コントロールとDIDの両方を特殊ケースとして含み、どちらよりも頑健な推定ができることが示されている。Rの synthdid パッケージが標準実装で、Pythonでは pysyncon が合成コントロール系の手法を提供している。

4. 介入ユニットが1つのときの推論：プラセボ検定

4.3.2で見たように、介入されたクラスターが1つ（カリフォルニア州だけ）のとき、クラスター標準誤差は信頼できない。Abadie et al. (2010) が用いた代替手段が プラセボ検定（permutation test） である。

考え方は単純で、「介入されていない各州を、順番に介入されたと見なして同じ分析を行い、本物の効果がその分布の中でどれだけ極端かを見る」というものである。効果の大きさの指標には、介入後の再現誤差と介入前の再現誤差の比（RMSPE比）を使う。介入前は再現がうまくいき、介入後だけ乖離が大きい州ほど、この比が大きくなる。

スクリプトでは36の候補州すべてで合成コントロールを繰り返し、RMSPE比を並べた。

カリフォルニア州は37州中5位で、プラセボ p 値はおよそ $5/37=0.135$ である。Abadie et al. (2010) の本来の分析では、共変量も使ってより良い再現をし、かつ2000年までの長い介入後期間を使っているため、カリフォルニア州が最上位（p 値 $1/39$ ）になっている。本ノートの簡易版では介入後が5年しかなく、効果が最大になる1995年以降を含んでいないため、順位が下がっている。 結果が違うこと自体が重要な教訓である。 合成コントロールの結論は、候補ユニットの選び方、一致させる変数、介入後の期間の長さに敏感であり、これらを変えた結果を併記するのが誠実な報告である。

DIDでクラスターが少ない場合には、 wild cluster bootstrap（Cameron, Gelbach, and Miller, 2008）も使われる。クラスターごとに残差の符号をランダムに反転させて疑似データを作り、統計量の分布を得る方法で、クラスター数が5〜10程度でも比較的うまく働くことが知られている。Conley and Taber (2011) は、介入クラスターが少数で非介入クラスターが多数ある場合に、非介入クラスターの残差の分布から介入効果の区間を作る方法を提案している。

5. 地域・時間単位の実験を設計する：Geo実験とスイッチバック

第4章の手法は「すでに起きた地域単位の介入」を事後分析するものだが、 最初から地域や時間の単位で実験を設計する 方法も発展している。

• Geo実験（Vaver and Koehler, 2011）は、地域をランダムに介入群と非介入群に割り当て、DIDや合成コントロールで効果を測る。Metaの GeoLift はこの設計を支援するオープンソースパッケージで、介入前データから「どの地域を介入群にすれば検出力が最大になるか」を選び、Augmented Synthetic Controlで分析する。
• スイッチバック実験 は、時間をブロック（たとえば1時間や1日）に分け、ブロックごとにランダムに介入と非介入を切り替える。フードデリバリーや配車サービスのように、ユーザ単位では干渉が避けられないマーケットプレイスで標準的に使われている。4.4で述べた「人気度の共有」のような干渉があるHuluの推薦施策でも、デバイス単位・日単位のスイッチバックは有効な選択肢である。分析にはブロック内の自己相関を考慮した推定（Bojinov, Simchi-Levi, and Zhao, 2023）が必要になる。

6. Pythonのライブラリ

本ノートでは学習のために回帰と最適化を自分で書いたが、実務では次のライブラリが使える。

設計例1：Fire TVだけに先行導入した新ホーム構成（DID）

この設計例を、profile_id × 時間のログから実際にパネルを組んで分析する手順（群の固定、0埋め、イベントスタディ、固定効果回帰、プラセボ検定、報告の書き方）は、実践編ノート ch4_hulu_device_rollout_did.md に架空データつきでまとめている。ここでは要点だけを示す。

アウトカム 。日次・デバイス別の「ホームのトレイ経由で再生したプロフィールの割合」または「1プロフィールあたり視聴シリーズ数」。前者はGA4のクリックログと content_viewing の結合から、後者は personalize.daily_personalize_viewing_kpi の series_cnt をデバイスで切って作れる。 content_viewing の device_low_category で fire_tv と android_tv を区別する。

比較相手 。介入前12週間のアウトカムの推移がFire TVともっとも平行なデバイスを選ぶ。候補はAndroid TV、Apple TV、ゲーム機など、同じ「テレビ画面」のデバイスである。スマートフォンは利用文脈が違い、トレンドの平行性が疑わしい。

モデル 。日次・デバイス別のパネルに対して、デバイス固定効果と日付固定効果を入れたDID回帰を行う。式は4.3.2の年ダミー付きモデルと同じで、 ca をFire TVダミー、年ダミーを日付ダミーに置き換える。

点検事項 。

• 介入前のイベントスタディで、 $\delta_k$ が0付近か。
• 同時施策の有無。Fire TVアプリのバージョンアップ、Fire TV向けキャンペーンなどが同時期になかったか。
• 干渉。世帯内でFire TVとスマホを併用するプロフィールが多いなら、スマホの視聴も新ホームの影響を受け得る。比較相手はテレビ系デバイスに限る理由の一つである。
• 標準誤差。デバイスは数個しかないので、クラスター標準誤差は信頼できない。他のデバイスを順に「介入された」と見なすプラセボ検定を併用する。

設計例2：全ユーザへの新シーズン訴求（CausalImpact）

アウトカム 。日次の「対象作品を25%以上視聴したプロフィール数」（ unique_viewed_series の登録基準に合わせる）。

予測変数 。介入前90日以上の期間で、次の系列を候補にする。訴求の影響を受けない（あるいは受けにくい）ことが条件である。

• サービス全体の日次アクティブプロフィール数。
• 対象作品と視聴者層の重なりが小さい他ジャンル作品群の日次視聴プロフィール数。同ジャンルの作品はカニバリゼーションで介入の影響を受ける可能性がある。
• 曜日・祝日のダミー。

やってはいけないこと 。対象作品の詳細ページ閲覧数や検索数を予測変数に入れること。これらは訴求の直接の結果であり、入れると効果が消える。

出力の読み方 。pointwiseの効果で「訴求の効果が何日間持続したか」を確認する。訴求の翌日にピークがあり1週間で0に戻るのか、新シーズンの評判で効果が持続するのかは、次の施策の頻度設計に直結する。累積効果を「訴求によって増えた視聴プロフィール数」として報告する。

限界の明示 。訴求と同日に地上波放送や配信開始があれば、その効果と分離できない（アクティビティバイアス）。報告書には同日に起きたイベントを必ず列挙する。

問1：前後比較で効果があったように見えた。DIDを使うと何が変わるか

前後比較の差には、介入の効果と「介入がなくても起きたはずの変化」（時間の効果）が混ざっている。DIDは非介入グループの前後の変化を時間の効果の推定値として使い、それを差し引く。介入がなくても両グループが同じだけ変化する（平行トレンド）という仮定のもとで、残るのは介入の効果だけになる。

問2：DID回帰で post の係数が有意でなかった。介入の効果はなかったと言えるか

言えない。 post の係数は「非介入グループにおける介入後の変化」、すなわち時間の効果であり、介入の効果ではない。介入の効果は交互作用項 treat:post の係数である。年ダミーを入れたモデルでは、 post の値はそもそも基準年の選び方で変わり、単独では意味を持たない。

問3：Fire TVの日次データを100日分使ったのに、なぜ標準誤差を通常の方法で計算してはいけないのか

同じデバイスの連続する日のアウトカムは互いに似ており（自己相関）、100日分は100個の独立な情報ではない。通常の標準誤差は独立を仮定するので、実際よりも小さく計算され、有意な結果が出すぎる。デバイスをクラスターとした標準誤差を使うのが原則だが、デバイスが数個しかなく介入デバイスが1つなら、クラスター標準誤差も信頼できない。その場合はプラセボ検定を使う。

問4：CausalImpactの予測変数に「対象作品の詳細ページ閲覧数」を使うと何が起きるか

訴求は詳細ページの閲覧を増やすので、この変数は介入の影響を受ける。介入後の予測値も訴求の影響で高くなり、実測との差である効果は過小に、極端な場合は0に推定される。第2章のPost treatment biasと同じ構造であり、介入の影響を受けない変数だけを候補にする。

問5：合成コントロールとDIDの推定値が大きく違った。どちらを信じるべきか

まず介入前の再現を比べる。DIDの比較相手（他州の平均）が介入前に平行でなく、合成コントロールが介入前をよく再現しているなら、合成コントロールの方が平行トレンドの根拠が強い。ただし合成コントロールも、候補ユニットの選び方や介入後の期間の長さに敏感である。両方の結果と、イベントスタディによる介入前の確認、プラセボ検定の順位を併記して報告するのが望ましい。