中断時系列デザインの理論 — 政策評価のための準実験デザイン

ある都市が飲酒運転の罰則を強化した後、交通事故件数はどう変わったか。ある国が予防接種プログラムを開始した後、感染症の発生率はどう変わったか。国が最低賃金を引き上げた後、雇用にどのような影響があったか。これらの問いに共通するのは、特定の時点で行われた介入（政策変更）の因果効果を、時系列データから推定したいということです。

理想的にはランダム化比較試験（RCT）を行いたいところですが、政策介入を都市や国にランダムに割り当てることは通常不可能です。また、差の差法（DID）を使いたくても、適切な対照群が見つからないこともあります。

このような状況で威力を発揮するのが中断時系列デザイン（Interrupted Time Series, ITS）です。ITSは、介入前の時系列トレンドを「反事実」として利用し、介入後のデータとの乖離から因果効果を推定する準実験デザイン（quasi-experimental design）です。対照群が不要であるという大きな利点を持ちます。

ITSは以下の分野で広く活用されています。

公衆衛生: 禁煙法施行後の心筋梗塞入院率の変化、予防接種プログラムの感染症減少効果
交通安全: 速度制限変更やシートベルト法の交通事故死亡率への影響
医療政策: 診療ガイドライン変更後の処方パターンの変化
経済政策: 税制変更や最低賃金引き上げの経済指標への影響
環境政策: 排出規制導入後の大気質の変化

本記事の内容

中断時系列デザインの基本的なアイデアと直感的理解
セグメンテッド回帰分析の数学的定式化
レベル変化とトレンド変化の解釈
介入効果の反事実的推定
ITSの前提条件と識別仮定
自己相関への対処法
対照群付きITS（Controlled ITS）
非線形トレンドへの拡張
多重介入への拡張
ITSの限界とDIDとの比較
Pythonによるシミュレーションと推定

前提知識

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

差の差法（DID）の理論 — 政策評価手法の比較対象
因果推論入門 — 因果推論の基本的な枠組み
重回帰分析 — 回帰分析の基礎

中断時系列デザインの基本アイデア

時系列トレンドを「対照群」として使う

差の差法では、介入を受けなかった対照群のトレンドを反事実として利用しました。ITSではこれを一歩進めて、介入前の同じ集団のトレンドを反事実として利用します。

考え方はシンプルです。もし介入がなければ、介入前のトレンド（上昇傾向や減少傾向）がそのまま続いていたはずだと仮定します。介入前のトレンドを外挿（extrapolation）して「介入がなかった世界」のアウトカムを予測し、実際の介入後データとの差を因果効果として推定するのです。

日常的なアナロジーで考えると、体重を記録している人が、ある日からダイエットを始めた場面を想像してください。ダイエット前は毎月0.5kgずつ体重が増えていたとします。ダイエット開始後、体重の増加が止まったり減り始めたりしたら、「ダイエットがなければ続いていたはずの増加傾向」との差がダイエットの効果です。

2つの効果: レベル変化とトレンド変化

ITSで推定する因果効果は、大きく2つの成分に分けられます。

レベル変化（level change）: 介入の即時効果。介入時点でアウトカムが「ジャンプ」するかどうか。例えば、禁煙法が施行された月に心筋梗塞入院率がガクンと下がるような場合です。

トレンド変化（trend change）: 介入の持続的効果。介入後にアウトカムの変化速度（傾き）が変わるかどうか。例えば、予防接種プログラム開始後、感染症の減少速度が以前より速くなるような場合です。

多くの介入は両方の効果を持ちます。法律の施行は即座に行動を変える（レベル変化）と同時に、長期的な傾向も変える（トレンド変化）可能性があります。一方、効果が徐々に現れる介入ではレベル変化が小さくトレンド変化が主な効果となることもあります。

この2つの効果を同時に推定する方法が、次に説明するセグメンテッド回帰分析です。

セグメンテッド回帰分析の定式化

基本モデル

ITSの最も標準的な分析方法はセグメンテッド回帰（segmented regression）です。時系列データを介入前と介入後の2つのセグメントに分割し、それぞれに異なる線形トレンドを当てはめます。

$$ \begin{equation} Y_t = \beta_0 + \beta_1 t + \beta_2 D_t + \beta_3 P_t + \epsilon_t \end{equation} $$

各変数と係数の意味を詳しく説明します。

$Y_t$: 時点 $t$ でのアウトカム（例: 月ごとの事故件数）
$t$: 時点インデックス（$t = 1, 2, \ldots, T$）。介入前からの経過時間を表す
$T_0$: 介入が行われた時点
$D_t$: 介入後ダミー変数。$D_t = \mathbf{1}(t \geq T_0)$（介入時点 $T_0$ 以降で1、それ以前で0）
$P_t$: 介入後の経過時間。$P_t = (t – T_0) \cdot D_t$（介入前は0、介入後は $t – T_0$）
$\epsilon_t$: 誤差項

各回帰係数の意味は次の通りです。

$\beta_0$: 介入前のベースラインレベル（$t = 0$ での切片）
$\beta_1$: 介入前のトレンド（1時点あたりのアウトカムの変化量）
$\beta_2$: レベル変化 — 介入直後のアウトカムの即座の変化量。介入の即時効果
$\beta_3$: トレンド変化 — 介入後のトレンドの変化量。介入前のトレンド $\beta_1$ に $\beta_3$ が加わり、介入後のトレンドは $\beta_1 + \beta_3$

介入前後のトレンドの明示

モデルを介入前と介入後に分けて書くと、構造がより明確になります。

介入前（$t < T_0$）:

$$ E[Y_t] = \beta_0 + \beta_1 t $$

アウトカムは切片 $\beta_0$ から時間 $t$ に比例して変化します。

介入後（$t \geq T_0$）:

$$ \begin{align} E[Y_t] &= \beta_0 + \beta_1 t + \beta_2 + \beta_3(t – T_0) \\ &= (\beta_0 + \beta_2 – \beta_3 T_0) + (\beta_1 + \beta_3) t \end{align} $$

介入後の切片は $\beta_0 + \beta_2 – \beta_3 T_0$、傾きは $\beta_1 + \beta_3$ です。

$\beta_2$ は介入時点 $t = T_0$ での回帰直線のジャンプの大きさを表し、$\beta_3$ は傾きの変化を表しています。

反事実の構成と因果効果の推定

介入がなかった場合の反事実（counterfactual）は、介入前のトレンドを介入後に外挿して構成します。

$$ \begin{equation} \hat{Y}_t^{\text{cf}} = \hat{\beta}_0 + \hat{\beta}_1 t \quad (t \geq T_0) \end{equation} $$

介入の因果効果は、実際の予測値と反事実の差として推定されます。

$$ \begin{equation} \hat{\tau}_t = \hat{E}[Y_t] – \hat{Y}_t^{\text{cf}} = \hat{\beta}_2 + \hat{\beta}_3(t – T_0) \quad (t \geq T_0) \end{equation} $$

この式は重要な含意を持ちます。

$t = T_0$ での因果効果: $\hat{\tau}_{T_0} = \hat{\beta}_2$（レベル変化そのもの）
$t > T_0$ での因果効果: 時間の経過とともに $\hat{\beta}_3$ の分だけ効果が蓄積する
$\hat{\beta}_3 < 0$（アウトカムが減少する方向の効果）の場合: 時間が経つほど因果効果（減少の大きさ）が増大する

ここで重要なのは、反事実の構成が介入前トレンドの外挿に完全に依存しているということです。介入前トレンドが正しく推定されていなければ、因果効果の推定も歪みます。この点は後ほど仮定と限界のセクションで詳しく議論します。

それでは、ITSが成り立つための条件を数学的に整理しましょう。

ITSの識別仮定

仮定1: 介入前トレンドの継続性（反事実の妥当性）

ITSの最も核心的な仮定は、介入がなければ介入前のトレンドがそのまま継続していたはずというものです。

$$ \begin{equation} E[Y_t(0)] = \beta_0 + \beta_1 t \quad \text{for all } t \end{equation} $$

ここで $Y_t(0)$ は介入がなかった場合のポテンシャルアウトカムです。この仮定が破れる典型的な状況は以下の通りです。

非線形トレンド: 実際のトレンドが線形ではなく、介入前の線形近似がそのまま外挿できない場合。例えば、感染症の流行が自然にピークを迎えて減少に転じるタイミングに介入が重なった場合
季節変動: 季節的なパターンを考慮しないと、介入効果と季節変動を混同してしまう
平均回帰: アウトカムが偶然に極端な値を取ったタイミングで介入が行われた場合、平均回帰を介入効果と誤認するリスクがある

仮定2: 同時イベントの不在

介入時点 $T_0$ の前後で、介入以外の大きなイベントが発生していないことを仮定します。もし介入と同時期に別の重要な変化（他の政策変更、経済危機、自然災害など）があれば、観測される変化が介入の効果なのか他のイベントの効果なのか区別できません。

これは因果推論の一般的な問題ですが、ITSでは対照群がないため特に深刻です。DIDでは対照群も同じ外部ショックを受けるため、その影響が差し引かれますが、ITSではそのようなメカニズムがありません。

仮定3: 介入の予測と先行効果の不在

介入が事前に予測されておらず、介入前に効果が先行して現れていないことを仮定します。

もし人々が介入を予測して行動を変えれば（例: 罰則強化が公表された時点で飲酒運転を控え始める）、介入前のデータが既に介入の影響を受けており、介入前トレンドの推定が歪みます。

仮定4: 構成変化の不在

時系列データの対象集団が介入前後で変わらないことを仮定します。例えば、ある病院のデータを使う場合、介入前後で患者の構成（年齢、重症度など）が大きく変わっていないことが必要です。

これらの仮定は直接検証できないものもありますが、いくつかの診断方法が存在します。その点は後ほど議論します。

次に、ITSの実践上の重要な課題である自己相関への対処について説明します。

自己相関への対処

時系列データの自己相関

OLS（通常の最小二乗法）によるセグメンテッド回帰は、誤差項 $\epsilon_t$ が互いに独立であることを仮定しています。しかし、時系列データでは通常、自己相関（autocorrelation）が存在します。つまり、$\text{Corr}(\epsilon_t, \epsilon_{t-k}) \neq 0$（$k \neq 0$）です。

例えば、月ごとの交通事故件数は、ある月の値が前月の値に影響される傾向があります。天候や季節要因の持続、事故報告の遅延なども自己相関を生みます。

自己相関がある場合、OLS推定量は依然として不偏ですが、標準誤差が過小推定されます。その結果、検定統計量が過大になり、介入効果が統計的に有意であるという誤った結論（偽陽性）を導きやすくなります。

対処法1: ニューイ・ウェスト標準誤差

ニューイ・ウェスト標準誤差（Newey-West standard errors）は、HAC（Heteroskedasticity and Autocorrelation Consistent）推定量とも呼ばれ、自己相関と不均一分散に対してロバストな標準誤差を提供します。

OLSの共分散行列推定量を次のように修正します。

$$ \begin{equation} \hat{V}_{\text{NW}} = (\bm{X}^T\bm{X})^{-1} \hat{\bm{\Omega}} (\bm{X}^T\bm{X})^{-1} \end{equation} $$

ここで $\hat{\bm{\Omega}}$ は自己相関を考慮したカーネル重み付き推定量です。

$$ \hat{\bm{\Omega}} = \sum_{j=-L}^{L} w(j/L) \hat{\bm{\Gamma}}_j $$

$w(\cdot)$ はバートレットカーネル等の重み関数、$L$ はラグの打ち切り次数（bandwidth）、$\hat{\bm{\Gamma}}_j$ は $j$ 次のサンプル自己共分散行列です。

対処法2: AR(p)誤差モデル

誤差項にAR(p)構造を仮定してモデルを拡張する方法です。

$$ \begin{align} Y_t &= \beta_0 + \beta_1 t + \beta_2 D_t + \beta_3 P_t + u_t \\ u_t &= \phi_1 u_{t-1} + \phi_2 u_{t-2} + \cdots + \phi_p u_{t-p} + \epsilon_t \end{align} $$

ここで $\epsilon_t \sim N(0, \sigma^2)$ はホワイトノイズです。この方法は、コクラン・オーカット法やプレイス・ウィンステン法で推定できます。

対処法3: 一般化最小二乗法（GLS）

自己相関構造を明示的にモデル化し、GLSで推定する方法です。誤差項の共分散行列 $\bm{\Sigma}$ を推定し、変換したモデルにOLSを適用します。

$$ \bm{\Sigma}^{-1/2} \bm{Y} = \bm{\Sigma}^{-1/2} \bm{X} \bm{\beta} + \bm{\Sigma}^{-1/2} \bm{\epsilon} $$

実用上は、AR(1)構造を仮定して自己相関係数 $\rho$ を推定するアプローチがよく使われます。

これらの対処法の中で、ニューイ・ウェスト標準誤差が最も広く使われています。OLS推定量自体は変更せず、標準誤差のみを修正するため、実装が容易です。

次に、ITSの頑健性を高めるための重要な拡張として、対照群を組み合わせたデザインを紹介します。

対照群付きITS（Controlled ITS）

単独ITSの限界

前述の通り、単独のITS（single-group ITS）は「介入以外に何も変わっていない」という強い仮定に依存しています。この仮定を弱めるために、対照群を導入したITSデザインが開発されました。

CITS（Controlled Interrupted Time Series）の定式化

対照群付きITS（CITS）は、介入を受けたグループと受けなかったグループの両方について時系列データを収集し、両者のトレンド変化の差を推定します。

$$ \begin{equation} Y_{gt} = \beta_0 + \beta_1 t + \beta_2 D_t + \beta_3 P_t + \beta_4 G_g + \beta_5 G_g \cdot t + \beta_6 G_g \cdot D_t + \beta_7 G_g \cdot P_t + \epsilon_{gt} \end{equation} $$

ここで: – $G_g$: 介入群ダミー（介入群 = 1、対照群 = 0） – $\beta_4$: 介入群と対照群のベースラインの差 – $\beta_5$: 介入群と対照群の介入前トレンドの差 – $\beta_6$: 介入群に固有のレベル変化（因果効果のレベル成分） – $\beta_7$: 介入群に固有のトレンド変化（因果効果のトレンド成分）

$\beta_6$ と $\beta_7$ がCITSで推定したい因果効果です。対照群も同じ外部環境の影響を受けるため、$\beta_2$（両群共通のレベル変化）と $\beta_3$（両群共通のトレンド変化）を通じてその影響が制御されます。

CITSとDIDの関係

CITSは、DIDとITSを統合したデザインと見なせます。

デザイン	時点	グループ	推定される効果
ITS	複数時点	1群のみ	レベル変化 + トレンド変化
DID	2時点（前後）	2群（介入・対照）	レベル変化
CITS	複数時点	2群（介入・対照）	レベル変化 + トレンド変化

CITSは時間軸の情報（複数時点のトレンド）とグループ間比較の両方を活用するため、ITSよりも頑健であり、DIDよりも柔軟です。ただし、適切な対照群の存在が前提条件となります。

ここまでの議論を拡張して、基本的な線形トレンドでは対処しきれない場合の非線形トレンドへの対応を考えてみましょう。

非線形トレンドと季節性への拡張

非線形トレンドモデル

介入前のトレンドが線形でない場合、線形セグメンテッド回帰の外挿が不適切になります。例えば、感染症の流行曲線は指数関数的な成長と減衰を示し、線形モデルでは捉えきれません。

非線形トレンドに対処する方法として、以下が挙げられます。

多項式トレンド:

$$ Y_t = \beta_0 + \beta_1 t + \beta_{1′} t^2 + \beta_2 D_t + \beta_3 P_t + \beta_{3′} P_t^2 + \epsilon_t $$

二次項を加えることで、曲線的なトレンドを捉えます。ただし、過剰適合のリスクがあるため、次数の選択には注意が必要です。

スプライン関数:

より柔軟な非線形トレンドには、制限付き3次スプライン（restricted cubic spline）が有効です。ノット（節点）の位置を適切に設定することで、複雑なトレンドを滑らかに近似できます。

季節性の制御

月次や週次の時系列データには季節変動が含まれることが多いです。季節性を制御しないと、季節的な上昇・下降を介入効果と混同するリスクがあります。

季節性はダミー変数やフーリエ級数で制御します。

$$ Y_t = \beta_0 + \beta_1 t + \beta_2 D_t + \beta_3 P_t + \sum_{s=1}^{S-1} \gamma_s M_{st} + \epsilon_t $$

ここで $M_{st}$ は月ダミー変数（月次データの場合 $S = 12$）です。フーリエ級数を使う場合は、

$$ Y_t = \beta_0 + \beta_1 t + \beta_2 D_t + \beta_3 P_t + \sum_{k=1}^{K} \left[\alpha_k \sin\left(\frac{2\pi k t}{S}\right) + \delta_k \cos\left(\frac{2\pi k t}{S}\right)\right] + \epsilon_t $$

フーリエ級数の利点は、少ないパラメータで季節変動を近似できることです。$K = 1$ で基本的な年周期パターンを、$K = 2$ でより複雑なパターンを捉えられます。

次に、1つの時系列上に複数の介入がある場合の拡張を考えます。

多重介入への拡張

複数の介入がある場合

現実には、時系列上に複数の介入が存在することが珍しくありません。例えば、交通安全対策として、最初にシートベルト法が導入され（時点 $T_1$）、その後飲酒運転の罰則が強化される（時点 $T_2$）場合です。

多重介入モデルは、各介入に対応するダミー変数と経過時間変数を追加して構成します。

$$ \begin{equation} Y_t = \beta_0 + \beta_1 t + \sum_{k=1}^{K} [\beta_{2k} D_{kt} + \beta_{3k} P_{kt}] + \epsilon_t \end{equation} $$

ここで $D_{kt} = \mathbf{1}(t \geq T_k)$、$P_{kt} = (t – T_k) \cdot D_{kt}$ は第 $k$ 介入に対応する変数です。

$\beta_{2k}$ は第 $k$ 介入のレベル変化、$\beta_{3k}$ はトレンド変化を表します。各介入の効果を分離して推定できます。

ただし、介入間の間隔が短いと、各介入の効果を安定的に推定するのが難しくなります。一般的に、各セグメント（介入と介入の間の期間）に十分なデータポイントが必要です。

ITSの診断と妥当性の評価

介入前トレンドの検証

ITSの信頼性を高めるために、以下の診断を行うことが推奨されます。

視覚的検査: まず時系列データをプロットし、介入前のトレンドが線形で安定しているかを目視で確認します。明らかな非線形パターンや外れ値がないかを確認します。

介入前トレンドの安定性検定: 介入前データのみを使って、トレンドの安定性を統計的に検定します。例えば、介入前期間をさらに2つに分割し、前半と後半でトレンドが変わっていないかを検定します。

プラセボ検定

プラセボ検定（placebo test）は、介入がなかった時点に仮想的な介入点を設定し、ITSを適用する方法です。もし実際の介入以外の時点でも「効果」が検出されれば、推定結果の信頼性は低下します。

具体的には、介入前の複数の時点を仮想介入点として設定し、各時点でのレベル変化とトレンド変化を推定します。これらが統計的に有意でなければ、介入前トレンドの安定性が支持されます。

ダービン・ワトソン検定

残差の自己相関を検出するために、ダービン・ワトソン検定（Durbin-Watson test）を実施します。

$$ DW = \frac{\sum_{t=2}^{T} (e_t – e_{t-1})^2}{\sum_{t=1}^{T} e_t^2} $$

$DW \approx 2$ であれば自己相関なし、$DW < 2$ は正の自己相関、$DW > 2$ は負の自己相関を示唆します。有意な自己相関が検出された場合は、前述の対処法（ニューイ・ウェスト標準誤差やGLS）を適用する必要があります。

ここまでの理論を、Pythonシミュレーションで具体的に確認しましょう。

Pythonシミュレーション

シミュレーション1: 基本的なITSの実装

まず、基本的なセグメンテッド回帰によるITSを実装し、レベル変化とトレンド変化の推定を行います。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

np.random.seed(42)

# --- データ生成 ---
T_total = 100  # 全時点数
T0 = 50        # 介入時点

t = np.arange(1, T_total + 1)
D = (t >= T0).astype(float)          # 介入後ダミー
P = np.maximum(t - T0, 0) * D        # 介入後の経過時間

# 真のパラメータ
beta_0 = 50.0    # ベースライン
beta_1 = -0.3    # 介入前トレンド（緩やかな減少）
beta_2 = -8.0    # レベル変化（即時効果: 8ポイント減少）
beta_3 = -0.5    # トレンド変化（減少傾向が加速）

# データ生成
Y = beta_0 + beta_1 * t + beta_2 * D + beta_3 * P + np.random.normal(0, 2, T_total)

# --- セグメンテッド回帰 ---
X_design = np.column_stack([t, D, P])
reg = LinearRegression().fit(X_design, Y)

beta_0_hat = reg.intercept_
beta_1_hat, beta_2_hat, beta_3_hat = reg.coef_

print("=== 推定結果 ===")
print(f"beta_0 (ベースライン):    {beta_0_hat:.4f} (真値: {beta_0})")
print(f"beta_1 (介入前トレンド):  {beta_1_hat:.4f} (真値: {beta_1})")
print(f"beta_2 (レベル変化):      {beta_2_hat:.4f} (真値: {beta_2})")
print(f"beta_3 (トレンド変化):    {beta_3_hat:.4f} (真値: {beta_3})")

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

# (a) データとフィット
ax = axes[0]
ax.scatter(t, Y, s=10, alpha=0.5, color="gray", label="Observed data")

# 回帰線（介入前）
t_pre = t[t < T0]
y_pre_fit = beta_0_hat + beta_1_hat * t_pre
ax.plot(t_pre, y_pre_fit, "b-", linewidth=2, label="Pre-intervention trend")

# 回帰線（介入後）
t_post = t[t >= T0]
y_post_fit = (beta_0_hat + beta_1_hat * t_post
              + beta_2_hat + beta_3_hat * (t_post - T0))
ax.plot(t_post, y_post_fit, "r-", linewidth=2, label="Post-intervention trend")

# 反事実（介入前トレンドの外挿）
y_cf = beta_0_hat + beta_1_hat * t_post
ax.plot(t_post, y_cf, "b--", linewidth=1.5, alpha=0.7, label="Counterfactual")

# 介入線
ax.axvline(T0, color="gray", linestyle=":", linewidth=1.5, label="Intervention")

# レベル変化の矢印
ax.annotate("", xy=(T0, y_post_fit[0]), xytext=(T0, y_cf[0]),
            arrowprops=dict(arrowstyle="<->", color="green", lw=2))
ax.text(T0 + 2, (y_post_fit[0] + y_cf[0]) / 2,
        f"Level change\n= {beta_2_hat:.1f}", fontsize=9, color="green")

ax.set_xlabel("Time", fontsize=12)
ax.set_ylabel("Outcome Y", fontsize=12)
ax.set_title("Interrupted Time Series Analysis", fontsize=12)
ax.legend(fontsize=8, loc="lower left")
ax.grid(True, alpha=0.3)

# (b) 推定された因果効果の時間推移
ax = axes[1]
effect = beta_2_hat + beta_3_hat * (t_post - T0)
ax.plot(t_post, effect, "g-", linewidth=2)
ax.fill_between(t_post, 0, effect, alpha=0.2, color="green")
ax.axhline(0, color="black", linewidth=0.5)
ax.axhline(beta_2_hat, color="red", linestyle="--", linewidth=1,
           label=f"Immediate effect = {beta_2_hat:.1f}")
ax.set_xlabel("Time", fontsize=12)
ax.set_ylabel("Estimated causal effect", fontsize=12)
ax.set_title("Estimated intervention effect over time", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

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

左のグラフでは、灰色の点が観測データ、青の実線が介入前のフィット、赤の実線が介入後のフィット、青の破線が反事実（介入前トレンドの外挿）を表しています。緑の矢印がレベル変化 $\hat{\beta}_2$ に対応しています。

右のグラフは、推定された因果効果の時間推移を示しています。介入直後のレベル変化（約-8ポイント）に加えて、トレンド変化（約-0.5ポイント/時点）により効果が時間とともに蓄積していく様子がわかります。

推定値が真の値に近いことから、セグメンテッド回帰が正しくパラメータを回復していることが確認できます。

シミュレーション2: 自己相関の影響と対処

自己相関がある場合にOLSの標準誤差がどれだけ過小推定されるかを確認し、ニューイ・ウェスト標準誤差との比較を行います。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(0)

T_total = 100
T0 = 50
n_simulations = 2000

t = np.arange(1, T_total + 1)
D = (t >= T0).astype(float)
P = np.maximum(t - T0, 0) * D
X = np.column_stack([np.ones(T_total), t, D, P])

# 真のパラメータ（介入効果なし: beta_2 = 0, beta_3 = 0）
beta_true = np.array([50.0, -0.3, 0.0, 0.0])

rho_values = [0.0, 0.3, 0.6, 0.9]  # 自己相関係数

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

for idx, rho in enumerate(rho_values):
    ax = axes[idx // 2][idx % 2]
    ols_pvalues_level = []

    for _ in range(n_simulations):
        # AR(1)誤差を生成
        eps = np.zeros(T_total)
        eps[0] = np.random.normal(0, 2)
        for tt in range(1, T_total):
            eps[tt] = rho * eps[tt-1] + np.random.normal(0, 2 * np.sqrt(1 - rho**2))

        Y = X @ beta_true + eps

        # OLS推定
        XtX_inv = np.linalg.inv(X.T @ X)
        beta_hat = XtX_inv @ X.T @ Y
        residuals = Y - X @ beta_hat
        sigma2_hat = np.sum(residuals**2) / (T_total - 4)
        se_ols = np.sqrt(sigma2_hat * np.diag(XtX_inv))

        # beta_2のt検定（OLS標準誤差）
        t_stat = beta_hat[2] / se_ols[2]
        p_val = 2 * (1 - stats.t.cdf(np.abs(t_stat), T_total - 4))
        ols_pvalues_level.append(p_val)

    # 名目5%水準での棄却率
    rejection_rate = np.mean(np.array(ols_pvalues_level) < 0.05)

    ax.hist(ols_pvalues_level, bins=50, density=True, alpha=0.6,
            color='steelblue', edgecolor='white')
    ax.axhline(1.0, color='red', linestyle='--', linewidth=1.5,
               label='Uniform (no autocorrelation)')
    ax.set_xlabel('p-value', fontsize=11)
    ax.set_ylabel('Density', fontsize=11)
    ax.set_title(f'$\\rho = {rho}$, Rejection rate = {rejection_rate:.3f}', fontsize=12)
    ax.legend(fontsize=9)
    ax.set_xlim(0, 1)
    ax.grid(True, alpha=0.3)

plt.suptitle('Distribution of p-values under H0 (no intervention effect)\nUsing OLS standard errors',
             fontsize=13, y=1.02)
plt.tight_layout()
plt.savefig('its_autocorrelation.png', dpi=150, bbox_inches='tight')
plt.show()

このシミュレーションでは、真の介入効果がゼロ（$\beta_2 = 0$, $\beta_3 = 0$）の状況で、異なる自己相関係数 $\rho$ のもとでOLSのp値の分布を確認しています。

結果から以下のことが読み取れます。

$\rho = 0$（自己相関なし）の場合: p値は一様分布に従い（赤の破線と一致）、名目5%水準での棄却率も約5%です。OLSの標準誤差は正しく機能しています。
$\rho = 0.3$ の場合: p値の分布がやや左に偏り、棄却率が5%を上回り始めます。軽度の自己相関でも影響が出ています。
$\rho = 0.6$ の場合: 棄却率が大幅に上昇し、10%以上になることがあります。「効果がないのに効果あり」と結論してしまう偽陽性率が2倍以上になっています。
$\rho = 0.9$ の場合: 棄却率が非常に高くなり、20%を超えることもあります。強い自己相関のもとでOLSの標準誤差をそのまま使うことは極めて危険です。

シミュレーション3: 対照群付きITS（CITS）

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

np.random.seed(123)

T_total = 80
T0 = 40

t = np.arange(1, T_total + 1)
D = (t >= T0).astype(float)
P = np.maximum(t - T0, 0) * D

# 真のパラメータ
# 共通ショック: 介入時点で外部要因により両群のアウトカムが+3増加
common_shock_level = 3.0
common_shock_trend = 0.0

# 介入効果: 介入群のみに適用
true_effect_level = -5.0
true_effect_trend = -0.3

# 介入群
Y_treat = (40 + 0.2 * t
           + (common_shock_level + true_effect_level) * D
           + (common_shock_trend + true_effect_trend) * P
           + np.random.normal(0, 1.5, T_total))

# 対照群（共通ショックは受けるが介入効果は受けない）
Y_control = (35 + 0.15 * t
             + common_shock_level * D
             + common_shock_trend * P
             + np.random.normal(0, 1.5, T_total))

# --- 単独ITS（介入群のみ）---
X_single = np.column_stack([t, D, P])
reg_single = LinearRegression().fit(X_single, Y_treat)
single_level = reg_single.coef_[1]
single_trend = reg_single.coef_[2]

# --- CITS（対照群付き）---
# データを結合
Y_all = np.concatenate([Y_treat, Y_control])
G = np.concatenate([np.ones(T_total), np.zeros(T_total)])
t_all = np.concatenate([t, t])
D_all = np.concatenate([D, D])
P_all = np.concatenate([P, P])

X_cits = np.column_stack([t_all, D_all, P_all, G, G * t_all, G * D_all, G * P_all])
reg_cits = LinearRegression().fit(X_cits, Y_all)
cits_level = reg_cits.coef_[5]   # G * D の係数
cits_trend = reg_cits.coef_[6]   # G * P の係数

print("=== 推定結果の比較 ===")
print(f"真の介入効果: レベル = {true_effect_level}, トレンド = {true_effect_trend}")
print(f"単独ITS:      レベル = {single_level:.3f}, トレンド = {single_trend:.3f}")
print(f"CITS:         レベル = {cits_level:.3f}, トレンド = {cits_trend:.3f}")

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

# (a) 単独ITS
ax = axes[0]
ax.scatter(t, Y_treat, s=8, alpha=0.4, color="steelblue", label="Treatment group")
t_pre = t[t < T0]
t_post = t[t >= T0]
y_pre = reg_single.intercept_ + reg_single.coef_[0] * t_pre
y_post = (reg_single.intercept_ + reg_single.coef_[0] * t_post
          + reg_single.coef_[1] + reg_single.coef_[2] * (t_post - T0))
y_cf = reg_single.intercept_ + reg_single.coef_[0] * t_post
ax.plot(t_pre, y_pre, 'b-', linewidth=2)
ax.plot(t_post, y_post, 'r-', linewidth=2)
ax.plot(t_post, y_cf, 'b--', linewidth=1.5, alpha=0.7, label='Counterfactual')
ax.axvline(T0, color='gray', linestyle=':', linewidth=1.5)
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('Outcome', fontsize=12)
ax.set_title(f'Single ITS: Level = {single_level:.2f}, Trend = {single_trend:.3f}',
             fontsize=11)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# (b) CITS
ax = axes[1]
ax.scatter(t, Y_treat, s=8, alpha=0.4, color="steelblue", label="Treatment")
ax.scatter(t, Y_control, s=8, alpha=0.4, color="darkorange", label="Control")
ax.axvline(T0, color='gray', linestyle=':', linewidth=1.5)
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('Outcome', fontsize=12)
ax.set_title(f'CITS: Level = {cits_level:.2f}, Trend = {cits_trend:.3f}',
             fontsize=11)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

plt.suptitle(f'True effect: Level = {true_effect_level}, Trend = {true_effect_trend}',
             fontsize=13, y=1.02)
plt.tight_layout()
plt.savefig('cits_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

このシミュレーションでは、介入時点で共通ショック（外部要因による+3の増加）が両群に影響する状況を設定しています。真の介入効果はレベル変化-5、トレンド変化-0.3です。

結果から重要な点が読み取れます。

単独ITSのレベル変化の推定値は共通ショックの影響を含んでいる: 真のレベル変化は-5ですが、単独ITSでは共通ショック（+3）が含まれるため、推定値は約-2（= -5 + 3）になります。これは因果効果の過小推定です。
CITSは共通ショックを除去できる: 対照群も共通ショックの影響を受けるため、介入群と対照群の差を取ることで共通ショックが相殺されます。CITSのレベル変化の推定値は真の値-5に近くなっています。

この例は、対照群を組み合わせることでITSの頑健性が大幅に向上することを示しています。

ITSとDIDの比較

ITSとDIDは政策評価で最もよく使われる準実験デザインです。両者の関係を整理しておきましょう。

特徴	ITS	DID
対照群の要否	不要（単独ITS）/ 任意（CITS）	必要
時点数	多数の時点が必要	2時点（前後）でも可
トレンドの仮定	介入前トレンドの継続性	平行トレンド仮定
推定できる効果	レベル変化 + トレンド変化	レベル変化
外部ショックへの頑健性	低い（単独ITS）	高い
非線形トレンドへの対処	柔軟に対応可能	2時点では対処困難

どちらを使うべきかは、利用可能なデータと研究状況に依存します。対照群が利用可能で時点数が限られている場合はDID、対照群がなく時系列データが豊富な場合はITSが適しています。もちろん、両者を組み合わせたCITSが最も頑健です。

まとめ

本記事では、中断時系列デザイン（ITS）について以下を解説しました。

基本アイデア: 介入前のトレンドを外挿して反事実を構成し、介入後データとの乖離から因果効果を推定する
セグメンテッド回帰: $Y_t = \beta_0 + \beta_1 t + \beta_2 D_t + \beta_3 P_t + \epsilon_t$ により、レベル変化 $\beta_2$ とトレンド変化 $\beta_3$ を同時に推定する
識別仮定: 介入前トレンドの継続性、同時イベントの不在、先行効果の不在、構成変化の不在
自己相関への対処: ニューイ・ウェスト標準誤差、AR誤差モデル、GLS。自己相関を無視するとOLS標準誤差が過小推定され、偽陽性率が大幅に上昇する
対照群付きITS（CITS）: 対照群を導入して外部ショックを制御し、推定の頑健性を高める
拡張: 非線形トレンド（多項式、スプライン）、季節性の制御、多重介入への対応
ITSとDIDの比較: 対照群の有無、時点数、推定対象の違い

ITSは、特に公衆衛生や政策評価の分野で極めて実用的なデザインです。対照群を必要としないという大きな利点がある一方、介入前トレンドの継続性仮定に強く依存するため、その妥当性の検証（プラセボ検定など）が不可欠です。

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

差の差法（DID）の理論 — 対照群を用いた政策評価手法
合成コントロール法 — 対照群を人工的に構築する手法
実験計画法 — 実験デザインの基礎
回帰不連続デザイン — 閾値に基づく因果推論

中断時系列デザインの理論 — 政策評価のための準実験デザイン

前提知識

中断時系列デザインの基本アイデア

時系列トレンドを「対照群」として使う

2つの効果: レベル変化とトレンド変化

セグメンテッド回帰分析の定式化

基本モデル

介入前後のトレンドの明示

反事実の構成と因果効果の推定

ITSの識別仮定

仮定1: 介入前トレンドの継続性（反事実の妥当性）

仮定2: 同時イベントの不在

仮定3: 介入の予測と先行効果の不在

仮定4: 構成変化の不在

自己相関への対処

時系列データの自己相関

対処法1: ニューイ・ウェスト標準誤差

対処法2: AR(p)誤差モデル

対処法3: 一般化最小二乗法（GLS）

対照群付きITS（Controlled ITS）

単独ITSの限界

CITS（Controlled Interrupted Time Series）の定式化

CITSとDIDの関係

非線形トレンドと季節性への拡張

非線形トレンドモデル

季節性の制御

多重介入への拡張

複数の介入がある場合

ITSの診断と妥当性の評価

介入前トレンドの検証

プラセボ検定

ダービン・ワトソン検定

Pythonシミュレーション

シミュレーション1: 基本的なITSの実装

シミュレーション2: 自己相関の影響と対処

シミュレーション3: 対照群付きITS（CITS）

ITSとDIDの比較

まとめ

実験計画法 — 完全無作為化・乱塊法・ラテン方格

ネイマン・ピアソンの基本補題をわかりやすく導出する

中断時系列デザインの理論 — 政策評価のための準実験デザイン

前提知識

中断時系列デザインの基本アイデア

時系列トレンドを「対照群」として使う

2つの効果: レベル変化とトレンド変化

セグメンテッド回帰分析の定式化

基本モデル

介入前後のトレンドの明示

反事実の構成と因果効果の推定

ITSの識別仮定

仮定1: 介入前トレンドの継続性（反事実の妥当性）

仮定2: 同時イベントの不在

仮定3: 介入の予測と先行効果の不在

仮定4: 構成変化の不在

自己相関への対処

時系列データの自己相関

対処法1: ニューイ・ウェスト標準誤差

対処法2: AR(p)誤差モデル

対処法3: 一般化最小二乗法（GLS）

対照群付きITS（Controlled ITS）

単独ITSの限界

CITS（Controlled Interrupted Time Series）の定式化

CITSとDIDの関係

非線形トレンドと季節性への拡張

非線形トレンドモデル

季節性の制御

多重介入への拡張

複数の介入がある場合

ITSの診断と妥当性の評価

介入前トレンドの検証

プラセボ検定

ダービン・ワトソン検定

Pythonシミュレーション

シミュレーション1: 基本的なITSの実装

シミュレーション2: 自己相関の影響と対処

シミュレーション3: 対照群付きITS（CITS）

ITSとDIDの比較

まとめ

関連記事

do計算とバックドア基準 — パールの因果推論

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

媒介分析 — 直接効果と間接効果の因果的分解

実験計画法 — 完全無作為化・乱塊法・ラテン方格

ネイマン・ピアソンの基本補題をわかりやすく導出する