合成コントロール法の理論 — 反事実を人工的に構築する因果推論

2000年、カリフォルニア州はアメリカで最も包括的なタバコ規制法「Proposition 99」を導入しました。この法律がタバコの売上にどのような影響を与えたかを評価したいとします。しかし、カリフォルニアにこの政策を適用した一方で「もし適用しなかった場合」のカリフォルニアを観察することはできません。他の州を対照群にすることを考えますが、カリフォルニアは人口構成、所得水準、都市化率などの点で他のどの州とも大きく異なります。単一の州を対照群に選ぶのは不適切です。

合成コントロール法(Synthetic Control Method, SCM)は、この問題に対するエレガントな解決策を提供します。複数のコントロールユニット(他の州)の重み付き平均として、処置ユニット(カリフォルニア)の「もし処置を受けなかったら」という反事実を人工的に構築するのです。例えば、「ネバダの20% + コネチカットの30% + ユタの15% + …」のように組み合わせることで、処置前のカリフォルニアのアウトカムをよく再現する「合成カリフォルニア」を作り、処置後のギャップから因果効果を推定します。

Abadie and Gardeazabal(2003)およびAbadie, Diamond, and Hainmueller(2010)によって体系化されたこの手法は、政策評価の分野で急速に普及し、現在では以下の幅広い分野で活用されています。

  • 公衆衛生: タバコ規制法、銃規制法、大麻合法化の効果評価
  • 経済政策: 経済特区、最低賃金引き上げ、貿易協定の経済的影響
  • 環境政策: 排出規制の環境・経済効果
  • 社会政策: 移民政策、教育改革の影響
  • 企業分析: 企業合併・買収や戦略変更の影響評価

本記事の内容

  • 合成コントロール法のアイデアと直感的理解
  • パネルデータにおける反事実の構築原理
  • 重みの推定: 制約付き最適化問題の定式化
  • 予測変数を用いた重み推定の一般化
  • 合成コントロール法の前提条件と識別
  • 推論方法: プラセボテストとRMSPE比
  • DIDとの関係と比較
  • 最近の発展: ペナルティ付きSCM、一般化SCM
  • Pythonシミュレーションによる実装と可視化

前提知識

この記事を読む前に、以下の知識があると理解が深まります。

合成コントロール法のアイデア

問題の設定

$J + 1$ 個のユニット(国、地域、企業など)を考えます。ユニット1が時点 $T_0$ で処置(政策介入)を受け、ユニット $2, \ldots, J+1$ は処置を受けないコントロールユニットです。これらのコントロールユニットの集合をドナープール(donor pool)と呼びます。

各ユニット $j$ の時点 $t$ でのアウトカムを $Y_{jt}$ とします。時点は $t = 1, 2, \ldots, T$ で、処置前期間は $t = 1, \ldots, T_0$、処置後期間は $t = T_0 + 1, \ldots, T$ です。

ポテンシャルアウトカムの枠組みで表現すると、$Y_{1t}(0)$ はユニット1が処置を受けなかった場合のアウトカム(反事実)、$Y_{1t}(1)$ は処置を受けた場合のアウトカムです。処置後の期間($t > T_0$)では $Y_{1t} = Y_{1t}(1)$ が観測され、$Y_{1t}(0)$ は観測されません。

推定したい処置効果は次の通りです。

$$ \begin{equation} \tau_{1t} = Y_{1t}(1) – Y_{1t}(0) \quad (t > T_0) \end{equation} $$

$Y_{1t}(1)$ は直接観測されるので、$Y_{1t}(0)$(反事実)を推定することが課題です。

反事実の構築: 重み付き平均

合成コントロール法の核心は、コントロールユニットの重み付き平均として反事実 $Y_{1t}(0)$ を推定することです。

$$ \begin{equation} \hat{Y}_{1t}(0) = \sum_{j=2}^{J+1} w_j Y_{jt} \end{equation} $$

ここで $\bm{w} = (w_2, \ldots, w_{J+1})^T$ は以下の制約を満たす重みベクトルです。

$$ w_j \geq 0 \quad \forall j, \quad \sum_{j=2}^{J+1} w_j = 1 $$

非負性制約($w_j \geq 0$)と和が1の制約は、合成コントロールがドナープールの凸結合(convex combination)であることを保証します。この制約は重要な意味を持ちます。凸結合は外挿(extrapolation)を防ぎ、合成コントロールがドナープールの「内側」に位置することを保証します。外挿はモデルの仮定に強く依存するため、凸結合の制約は推定の頑健性を高めます。

なぜ単一の対照群ではなく「合成的な」対照群か

単一のコントロールユニットを選ぶ場合、処置ユニットとの類似性が限られます。例えば、カリフォルニアに最も似た州がネバダだとしても、ネバダとカリフォルニアは人口、産業構造、文化など多くの点で異なります。

合成コントロール法は、複数のユニットを適切に組み合わせることで、どの単一のユニットよりも処置ユニットに類似した「合成的な」対照群を構築できます。直感的には、カリフォルニアの人口規模にはニューヨーク、西海岸の特性にはオレゴン、産業構造にはイリノイが似ているかもしれません。これらの特徴を重み付きで組み合わせることで、カリフォルニアの多面的な特性をより正確に再現できるのです。

処置効果の推定値は次のようになります。

$$ \begin{equation} \hat{\tau}_{1t} = Y_{1t} – \hat{Y}_{1t}(0) = Y_{1t} – \sum_{j=2}^{J+1} w_j^* Y_{jt} \quad (t > T_0) \end{equation} $$

では、重みをどのように決定するかという核心的な問題に進みましょう。

重みの推定

処置前フィットに基づく推定

最も基本的なアプローチは、処置前の期間でアウトカムの再現誤差を最小化するように重みを選ぶことです。

$$ \begin{equation} \bm{w}^* = \arg\min_{\bm{w}} \sum_{t=1}^{T_0}\left(Y_{1t} – \sum_{j=2}^{J+1} w_j Y_{jt}\right)^2 \end{equation} $$ $$ \text{subject to} \quad w_j \geq 0, \quad \sum_{j=2}^{J+1} w_j = 1 $$

行列表記では、処置前のアウトカムを $\bm{Y}_1 = (Y_{11}, \ldots, Y_{1,T_0})^T$ と $\bm{Y}_0 = [Y_{jt}]_{j=2,\ldots,J+1;\ t=1,\ldots,T_0}$ で表すと、

$$ \bm{w}^* = \arg\min_{\bm{w}} \|\bm{Y}_1 – \bm{Y}_0^T \bm{w}\|_2^2 \quad \text{s.t.} \quad \bm{w} \geq \bm{0}, \quad \bm{1}^T\bm{w} = 1 $$

これは制約付き二次計画問題(constrained quadratic program)であり、凸最適化問題なので大域的最適解が効率的に求まります。

予測変数を用いた一般化

Abadie et al.(2010)は、アウトカムだけでなく予測変数(predictors)の再現も考慮した一般化を提案しました。

各ユニット $j$ について、$K$ 個の予測変数の値を含むベクトル $\bm{X}_j \in \mathbb{R}^K$ を定義します。予測変数には、処置前の期間のアウトカムの他に、人口、GDP、教育水準などの共変量を含めることができます。

重みは次の距離を最小化するように選ばれます。

$$ \begin{equation} \bm{w}^* = \arg\min_{\bm{w}} (\bm{X}_1 – \bm{X}_0 \bm{w})^T \bm{V} (\bm{X}_1 – \bm{X}_0 \bm{w}) \end{equation} $$ $$ \text{s.t.} \quad w_j \geq 0, \quad \sum w_j = 1 $$

ここで $\bm{V}$ は $K \times K$ の正定値対称行列で、各予測変数の重要度を表します。$\bm{V}$ が対角行列であれば、各予測変数に異なる重みを割り当てることに対応します。

$\bm{V}$ の選択方法としては、以下のアプローチがあります。

  1. クロスバリデーション: 処置前期間の一部を検証用に残し、予測精度を最大化する $\bm{V}$ を選ぶ
  2. 回帰ベース: アウトカムに対する各予測変数の説明力に基づいて $\bm{V}$ を決定する
  3. 恒等行列: $\bm{V} = \bm{I}$(全予測変数を等しく扱う)

重みのスパース性

凸結合の制約と処置前フィットの最適化の結果、最適重み $\bm{w}^*$ は多くの場合スパース(多くの成分がゼロ)になります。つまり、少数のコントロールユニットだけが正の重みを持ち、残りのユニットの重みはゼロです。

このスパース性は解釈の観点からも望ましい性質です。「合成カリフォルニアはネバダ30%、コネチカット25%、ユタ20%、コロラド15%、モンタナ10%で構成される」のような、透明性の高い記述が可能になります。

ここまでで合成コントロールの構築方法を理解しました。次に、合成コントロール法が正当な因果推定を与えるための条件を整理しましょう。

前提条件と識別

処置前フィットの良さ

合成コントロール法が有効に機能するための最も重要な条件は、処置前の期間で合成コントロールが処置ユニットのアウトカムをよく再現していることです。

処置前のフィットが良い(処置前のギャップが小さい)ことは、合成コントロールが処置ユニットの「反事実」として妥当であることの必要条件です。処置前のフィットが悪い場合、合成コントロール法の推定結果は信頼できません。

フィットの良さは処置前二乗平均平方根誤差(pre-treatment RMSPE)で定量化されます。

$$ \begin{equation} \text{RMSPE}_{\text{pre}} = \sqrt{\frac{1}{T_0}\sum_{t=1}^{T_0}(Y_{1t} – \hat{Y}_{1t}(0))^2} \end{equation} $$

SUTVA(安定ユニット処置値仮定)

ポテンシャルアウトカムの枠組みと同様に、合成コントロール法はSUTVAを必要とします。特に重要なのは干渉なし(no interference)の仮定です。処置ユニットへの介入がコントロールユニットのアウトカムに影響しないことを仮定します。

例えば、カリフォルニアのタバコ規制法がネバダのタバコ売上に影響(カリフォルニアの住民がネバダでタバコを購入する「スピルオーバー効果」)を与えれば、SUTVAは破れます。この場合、ドナープールから影響を受けるユニットを除外する必要があります。

コントロールユニットが処置を受けない

ドナープール内のユニットが、処置期間中に類似の処置を受けていないことを仮定します。もしコントロールユニットの一部が同様の政策を導入していれば、合成コントロールの反事実としての妥当性が損なわれます。

凸包の条件

処置ユニットの特性ベクトル $\bm{X}_1$ が、コントロールユニットの特性ベクトル $\{\bm{X}_j\}_{j=2}^{J+1}$ の凸包(convex hull)の内部にあることが望ましいです。処置ユニットがコントロールユニットの凸包の外にある場合、良好な処置前フィットが得られない可能性があります。

これらの条件が満たされたとして、推定結果の統計的有意性をどのように評価するのでしょうか。次にプラセボテストを紹介します。

推論方法: プラセボテスト

標準的な推論の困難さ

合成コントロール法では、処置ユニットが1つだけであることが普通です。処置ユニットの数が1では、通常の統計的検定($t$ 検定など)に必要なサンプルサイズの要件を満たしません。

プラセボテスト(In-space placebo test)

プラセボテストはパーミュテーション(置換)の考え方に基づく推論手法です。各コントロールユニットに対して「もしそのユニットが処置を受けていたら」という仮想的な分析を行い、推定された効果の分布を構成します。

手順:

  1. 処置ユニット(ユニット1)に対して合成コントロールを構築し、処置効果 $\hat{\tau}_{1t}$ を推定する
  2. 各コントロールユニット $j = 2, \ldots, J+1$ を「仮想的な処置ユニット」と見なし、残りのユニットから合成コントロールを構築する
  3. 各ユニット $j$ のプラセボ効果 $\hat{\tau}_{jt}$ を計算する
  4. 処置ユニットの推定効果が、プラセボ効果の分布の中でどの程度極端かを評価する

RMSPE比に基づくp値

効果の大きさを処置前フィットの良さで正規化したRMSPE比を用いることが推奨されます。

$$ \begin{equation} \text{ratio}_j = \frac{\text{RMSPE}_{\text{post}, j}}{\text{RMSPE}_{\text{pre}, j}} \end{equation} $$

ここで $\text{RMSPE}_{\text{pre}, j}$ と $\text{RMSPE}_{\text{post}, j}$ はユニット $j$ の処置前・処置後のRMSPEです。

処置前のフィットが悪いユニットは大きなプラセボ効果を持つ可能性がありますが、それは処置効果ではなくフィットの悪さに起因します。RMSPE比はこの問題を正規化によって緩和します。

$p$ 値は次のように計算されます。

$$ \begin{equation} p = \frac{|\{j : \text{ratio}_j \geq \text{ratio}_1\}|}{J + 1} \end{equation} $$

$J = 20$ のドナープールの場合、$p$ 値の最小値は $1/21 \approx 0.048$ です。

In-time プラセボテスト

空間方向のプラセボテストに加えて、時間方向のプラセボテスト(in-time placebo test)も重要な診断ツールです。処置前の任意の時点を仮想的な処置時点として設定し、合成コントロール法を適用します。もし処置前の仮想処置時点で「効果」が検出されれば、合成コントロール法のモデルが不適切であることを示唆します。

ここまでの理論を整理したところで、合成コントロール法とDIDの関係を明確にしておきましょう。

DIDとの関係と比較

DIDは合成コントロール法の特殊ケース

差の差法は、全てのコントロールユニットに均等な重み $w_j = 1/J$ を暗黙に割り当てていると解釈できます。合成コントロール法は、この重みをデータから最適に推定する一般化です。

特徴 DID SCM
対照群の構成 均等重み or 固定のグループ データから最適重みを推定
処置前フィットの要件 平行トレンド仮定 処置前の良好なフィット
処置ユニット数 複数可 通常1つ
時点数の要件 2時点(前後)でも可 多数の処置前時点が必要
推論方法 標準的な統計推論 プラセボテスト
透明性 対照群の選択が主観的 重みが明示的

合成コントロール法の利点

  1. 処置前フィットの可視化: 合成コントロールが処置ユニットをどの程度再現しているかを視覚的に確認でき、推定の妥当性を直感的に評価できる
  2. 重みの透明性: どのコントロールユニットがどの程度貢献しているかが明示される
  3. 平行トレンド仮定の緩和: DIDの平行トレンド仮定(全てのコントロールユニットが同じトレンドを持つ)よりも柔軟な仮定のもとで因果推定が可能

合成コントロール法の限界

  1. 処置ユニットが1つの場合に限定: 複数のユニットが同時に処置を受ける場合への拡張は活発な研究テーマ
  2. ドナープールのサイズに依存: コントロールユニットが少ないと、良好なフィットが得られない
  3. 負の重みを許さない: 凸結合の制約により、処置ユニットがドナープールの凸包の外にある場合に対応できない

最近の発展

ペナルティ付きSCM(Penalized SCM)

Abadie and L’Hour(2021)は、重みの推定にペナルティ項を追加する方法を提案しました。通常のSCMではスパースな重みが得られますが、処置前フィットが同程度に良い複数の重みベクトルが存在する場合(重みの非一意性)、推定結果が不安定になることがあります。

ペナルティ付きSCMは、処置前フィットの最小化に加えて、重みの分散(スパース性)にペナルティを加えることで、この非一意性の問題を緩和します。

一般化SCM(Augmented SCM)

Ben-Michael, Feller, and Rothstein(2021)は、拡張合成コントロール法(Augmented Synthetic Control Method)を提案しました。これは、合成コントロールの残差にバイアス補正項を追加することで、処置前フィットが完全でない場合でもバイアスを削減する手法です。

SCMとDIDの統合

Arkhangelsky et al.(2021)は、合成差の差法(Synthetic Difference-in-Differences, SDID)を提案しました。SDIDはSCMとDIDの利点を組み合わせ、単位方向と時間方向の両方に重みを最適化する手法です。

これらの発展を踏まえつつ、基本的なSCMをPythonで実装してみましょう。

Pythonシミュレーション

シミュレーション1: 基本的な合成コントロール法

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

np.random.seed(42)

# --- パラメータ ---
J = 15       # コントロールユニット数
T = 40       # 全期間数
T0 = 25      # 処置開始時点
true_effect = 5.0

# --- データ生成 ---
# 各ユニットの潜在ファクター
factors = np.random.normal(0, 1, (3, T))  # 3つの共通ファクター
loadings = np.random.normal(0, 1, (J + 1, 3))  # ファクターローディング

Y = np.zeros((J + 1, T))
for j in range(J + 1):
    trend = 0.2 * np.arange(T) + np.random.uniform(5, 15)
    Y[j] = trend + loadings[j] @ factors + np.random.normal(0, 0.5, T)

# 処置効果を追加(ユニット0のみ、処置後期間)
for t in range(T0, T):
    Y[0, t] += true_effect

# --- 合成コントロールの推定 ---
Y_pre_treated = Y[0, :T0]
Y_pre_control = Y[1:, :T0]

def objective(w):
    synthetic = Y_pre_control.T @ w
    return np.sum((Y_pre_treated - synthetic)**2)

constraints = [{"type": "eq", "fun": lambda w: np.sum(w) - 1}]
bounds = [(0, 1)] * J
w0 = np.ones(J) / J

result = minimize(objective, w0, method="SLSQP",
                  bounds=bounds, constraints=constraints)
w_opt = result.x

# 合成コントロール
Y_synthetic = Y[1:, :].T @ w_opt

# 処置効果の推定
gaps = Y[0] - Y_synthetic
estimated_effect = gaps[T0:]

# RMSPE
rmspe_pre = np.sqrt(np.mean(gaps[:T0]**2))
rmspe_post = np.sqrt(np.mean(gaps[T0:]**2))

print("=== 合成コントロール法の結果 ===")
print(f"処置前RMSPE: {rmspe_pre:.4f}")
print(f"処置後RMSPE: {rmspe_post:.4f}")
print(f"推定効果の平均: {estimated_effect.mean():.4f} (真値: {true_effect})")
print(f"\n重み(正の重みのみ):")
for idx in np.argsort(w_opt)[::-1]:
    if w_opt[idx] > 0.01:
        print(f"  ユニット {idx+2}: w = {w_opt[idx]:.4f}")

# --- 可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# (a) アウトカムの推移
ax = axes[0]
for j in range(1, J + 1):
    ax.plot(range(T), Y[j], color='gray', alpha=0.2, linewidth=0.5)
ax.plot(range(T), Y[0], "r-", linewidth=2, label="Treated unit")
ax.plot(range(T), Y_synthetic, "b--", linewidth=2, label="Synthetic control")
ax.axvline(T0, color="black", linestyle=":", linewidth=1.5, label="Intervention")
ax.set_xlabel("Time", fontsize=12)
ax.set_ylabel("Outcome", fontsize=12)
ax.set_title("Synthetic Control Method", fontsize=12)
ax.legend(fontsize=9, loc='upper left')
ax.grid(True, alpha=0.3)

# (b) ギャップ(処置効果)
ax = axes[1]
ax.plot(range(T), gaps, "k-", linewidth=2)
ax.axhline(0, color="gray", linewidth=0.5)
ax.axvline(T0, color="black", linestyle=":", linewidth=1.5, label="Intervention")
ax.axhline(true_effect, color="red", linestyle="--", linewidth=1.5,
           label=f"True effect = {true_effect}")
ax.fill_between(range(T0, T), 0, gaps[T0:], alpha=0.2, color="green")
ax.set_xlabel("Time", fontsize=12)
ax.set_ylabel("Gap (Treated - Synthetic)", fontsize=12)
ax.set_title("Estimated treatment effect", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# (c) 重みの分布
ax = axes[2]
nonzero_idx = np.where(w_opt > 0.01)[0]
labels = [f"Unit {i+2}" for i in nonzero_idx]
weights = w_opt[nonzero_idx]
ax.barh(labels, weights, color='steelblue', alpha=0.7, edgecolor='black')
ax.set_xlabel("Weight", fontsize=12)
ax.set_title("Synthetic Control Weights", fontsize=12)
ax.grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.savefig("synthetic_control_basic.png", dpi=150, bbox_inches="tight")
plt.show()

左のグラフでは、灰色の線がコントロールユニット、赤い実線が処置ユニット、青い破線が合成コントロールです。処置前の期間で合成コントロールが処置ユニットをよく追従しており、処置後に明確なギャップが生じていることがわかります。

中央のグラフは処置効果の時系列です。処置前のギャップはほぼゼロ(良好なフィット)であり、処置後に約5のギャップが安定的に維持されています。これは真の効果5.0と整合的です。

右のグラフは合成コントロールの重みを示しています。少数のユニットが大きな重みを持ち、残りのユニットの重みはゼロに近いというスパースなパターンが確認できます。

シミュレーション2: プラセボテスト

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

np.random.seed(42)

J = 15
T = 40
T0 = 25
true_effect = 5.0

# データ生成(前のシミュレーションと同じ)
factors = np.random.normal(0, 1, (3, T))
loadings = np.random.normal(0, 1, (J + 1, 3))
Y = np.zeros((J + 1, T))
for j in range(J + 1):
    trend = 0.2 * np.arange(T) + np.random.uniform(5, 15)
    Y[j] = trend + loadings[j] @ factors + np.random.normal(0, 0.5, T)
for t in range(T0, T):
    Y[0, t] += true_effect

def fit_synthetic(treated_idx, Y, T0, J_total):
    """指定されたユニットを処置ユニットとして合成コントロールを推定"""
    control_idx = [j for j in range(J_total + 1) if j != treated_idx]
    J_c = len(control_idx)

    Y_pre_t = Y[treated_idx, :T0]
    Y_pre_c = Y[control_idx, :T0]

    def obj(w):
        return np.sum((Y_pre_t - Y_pre_c.T @ w)**2)

    constraints = [{"type": "eq", "fun": lambda w: np.sum(w) - 1}]
    bounds = [(0, 1)] * J_c
    w0 = np.ones(J_c) / J_c

    result = minimize(obj, w0, method="SLSQP", bounds=bounds,
                      constraints=constraints, options={'maxiter': 1000})
    w = result.x

    Y_synth = Y[control_idx, :].T @ w
    gaps = Y[treated_idx] - Y_synth
    return gaps, w

# 全ユニットに対してプラセボテスト
all_gaps = {}
rmspe_ratios = {}

for j in range(J + 1):
    gaps, _ = fit_synthetic(j, Y, T0, J)
    all_gaps[j] = gaps
    rmspe_pre = np.sqrt(np.mean(gaps[:T0]**2))
    rmspe_post = np.sqrt(np.mean(gaps[T0:]**2))
    if rmspe_pre > 0:
        rmspe_ratios[j] = rmspe_post / rmspe_pre
    else:
        rmspe_ratios[j] = 0

# p値
ratio_treated = rmspe_ratios[0]
p_value = sum(1 for j in range(J + 1) if rmspe_ratios[j] >= ratio_treated) / (J + 1)

print(f"処置ユニットのRMSPE比: {ratio_treated:.4f}")
print(f"p値: {p_value:.4f}")

# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))

# (a) プラセボギャップ
ax = axes[0]
for j in range(1, J + 1):
    # 処置前RMSPEが大きすぎるユニットは除外
    rmspe_pre_j = np.sqrt(np.mean(all_gaps[j][:T0]**2))
    if rmspe_pre_j < 3 * np.sqrt(np.mean(all_gaps[0][:T0]**2)) + 1:
        ax.plot(range(T), all_gaps[j], color='gray', alpha=0.3, linewidth=0.8)
ax.plot(range(T), all_gaps[0], 'r-', linewidth=2.5, label='Treated unit')
ax.axvline(T0, color='black', linestyle=':', linewidth=1.5)
ax.axhline(0, color='black', linewidth=0.3)
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('Gap', fontsize=12)
ax.set_title(f'Placebo Test (p = {p_value:.3f})', fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# (b) RMSPE比の分布
ax = axes[1]
ratios = [rmspe_ratios[j] for j in range(J + 1)]
colors = ['red' if j == 0 else 'steelblue' for j in range(J + 1)]
sorted_idx = np.argsort(ratios)[::-1]
labels = [f'Unit {j+1}' if j > 0 else 'Treated' for j in sorted_idx]
vals = [ratios[j] for j in sorted_idx]
bar_colors = ['red' if j == 0 else 'steelblue' for j in sorted_idx]

ax.barh(range(len(vals)), vals, color=bar_colors, alpha=0.7, edgecolor='black')
ax.set_yticks(range(len(vals)))
ax.set_yticklabels(labels, fontsize=8)
ax.set_xlabel('Post/Pre RMSPE Ratio', fontsize=12)
ax.set_title('RMSPE Ratio Distribution', fontsize=12)
ax.grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.savefig('synthetic_control_placebo.png', dpi=150, bbox_inches='tight')
plt.show()

左のグラフは、処置ユニット(赤い太線)とコントロールユニットのプラセボギャップ(灰色の細線)を重ねたものです。処置ユニットのギャップが処置後に明確に他のユニットから乖離していることがわかります。プラセボテストのアイデアは、「処置を受けていないユニットでも同程度のギャップが生じるなら、処置ユニットのギャップは処置効果ではなく偶然かもしれない」ということです。

右のグラフはRMSPE比のランキングです。処置ユニット(赤)が最大のRMSPE比を持てば、p値は $1/(J+1)$ と最小になります。処置ユニットのRMSPE比が十分に大きく、プラセボユニットとの差が明確であれば、統計的に有意な処置効果が存在すると結論できます。

シミュレーション3: 処置前フィットの重要性

処置前フィットが良い場合と悪い場合で、推定結果がどう変わるかを比較します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

np.random.seed(0)

T = 40
T0 = 25
true_effect = 4.0

fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))

for idx, (J, title) in enumerate([(20, 'Good Fit (J=20)'),
                                    (3, 'Poor Fit (J=3)')]):
    factors = np.random.normal(0, 1, (3, T))
    loadings = np.random.normal(0, 1, (J + 1, 3))
    Y = np.zeros((J + 1, T))
    for j in range(J + 1):
        trend = 0.2 * np.arange(T) + np.random.uniform(5, 15)
        Y[j] = trend + loadings[j] @ factors + np.random.normal(0, 0.5, T)
    for t in range(T0, T):
        Y[0, t] += true_effect

    # 合成コントロール
    Y_pre_t = Y[0, :T0]
    Y_pre_c = Y[1:, :T0]

    def obj(w):
        return np.sum((Y_pre_t - Y_pre_c.T @ w)**2)

    constraints = [{"type": "eq", "fun": lambda w: np.sum(w) - 1}]
    bounds = [(0, 1)] * J
    w0 = np.ones(J) / J
    result = minimize(obj, w0, method="SLSQP", bounds=bounds,
                      constraints=constraints)
    w_opt = result.x
    Y_synth = Y[1:, :].T @ w_opt
    gaps = Y[0] - Y_synth

    rmspe_pre = np.sqrt(np.mean(gaps[:T0]**2))
    est_effect = gaps[T0:].mean()

    ax = axes[idx]
    ax.plot(range(T), Y[0], 'r-', linewidth=2, label='Treated')
    ax.plot(range(T), Y_synth, 'b--', linewidth=2, label='Synthetic')
    ax.axvline(T0, color='black', linestyle=':', linewidth=1.5)
    ax.set_xlabel('Time', fontsize=12)
    ax.set_ylabel('Outcome', fontsize=12)
    ax.set_title(f'{title}\nPre-RMSPE={rmspe_pre:.2f}, Est. effect={est_effect:.2f}',
                 fontsize=11)
    ax.legend(fontsize=10)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('scm_fit_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

2つのパネルから、ドナープールのサイズと処置前フィットの関係が明確に読み取れます。

  1. $J = 20$(左): 十分なコントロールユニットがあり、処置前のフィットが良好です。合成コントロールが処置ユニットのトレンドを忠実に再現しており、処置後のギャップが因果効果として信頼できます。
  2. $J = 3$(右): コントロールユニットが少ないため、処置前のフィットが悪化しています。処置前からギャップが存在し、処置後のギャップのうちどれだけが処置効果でどれだけがフィットの悪さに起因するかが不明確です。

この比較は、合成コントロール法において十分な数のコントロールユニットが処置前フィットの良さに不可欠であることを示しています。

まとめ

本記事では、合成コントロール法の理論を以下の流れで解説しました。

  • 基本アイデア: コントロールユニットの重み付き平均として処置ユニットの反事実を構築する。重みは処置前フィットを最適化するように制約付き凸最適化で推定
  • 凸結合の制約($w_j \geq 0$, $\sum w_j = 1$)は外挿を防ぎ、推定の頑健性を高める。結果として重みはスパースになりやすい
  • 前提条件: 処置前の良好なフィット、SUTVA、コントロールユニットが処置を受けないこと
  • プラセボテスト: 各コントロールユニットを仮想的な処置ユニットとして分析し、RMSPE比の分布から統計的有意性を評価する
  • DIDとの関係: SCMはDIDの一般化であり、重みをデータから最適化することで平行トレンド仮定を緩和する
  • 最近の発展: ペナルティ付きSCM、拡張SCM、合成差の差法(SDID)

合成コントロール法は、処置ユニットが少数で適切な対照群が見つからない政策評価の場面で極めて有用な手法です。推定の透明性(重みの明示性)と処置前フィットの可視化が、結果の信頼性評価を直感的に可能にする点が大きな強みです。

次のステップとして、以下の記事も参考にしてください。