観察研究で因果効果を推定したいが、重要な交絡変数が観測できない——このような状況は現実にしばしば起こります。傾向スコア法や回帰調整は、すべての交絡変数が観測されているという仮定(無視可能性)に依存しています。しかし、個人の動機、遺伝的素因、測定されていない環境要因などの未観測交絡が存在する場合、これらの方法ではバイアスが残ります。
操作変数法(Instrumental Variables method, IV)は、未観測の交絡が存在する場合でも因果効果を推定できる強力な手法です。このアイデアは、処置と直接アウトカムには影響せず、処置の割り当てのみに影響する「てこ(instrument)」のような変数を利用するものであり、計量経済学の分野で広く発展してきました。
操作変数法は以下の場面で活用されます。
- 経済学: 教育の収益率の推定(カードの兵役くじの研究)、制度の効果分析
- 疫学: メンデルのランダム化(遺伝子変異を操作変数として使用)
- マーケティング: 広告の因果効果の推定(広告枠の価格変動を操作変数として使用)
- 公衆衛生: 近接性を操作変数とした医療サービスの効果分析
本記事では、操作変数の定義と3つの条件、ワルド推定量の導出、2段階最小二乗法(2SLS)の理論、そしてLATE(局所的平均処置効果)の解釈を解説します。
本記事の内容
- 操作変数の定義と3条件(関連性、排除制約、独立性)
- ワルド推定量の導出と直感
- 2段階最小二乗法(2SLS)の理論と導出
- LATE解釈と遵守者の概念
- 弱い操作変数の問題
- Pythonによる2SLSの実装とシミュレーション
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 交絡バイアスとその制御法 — 未観測交絡の問題
- 平均処置効果(ATE) — 因果効果の定義
- 重回帰分析 — OLS推定量の性質
内生性と操作変数の必要性
内生性問題
まず、操作変数法が必要になる場面を確認しましょう。次の線形モデルを考えます。
$$ \begin{equation} Y_i = \beta_0 + \beta_1 T_i + \epsilon_i \end{equation} $$
ここで $T_i$ は処置変数、$Y_i$ はアウトカム、$\beta_1$ は推定したい因果効果です。もし $T_i$ と $\epsilon_i$ が相関していれば($\text{Cov}(T, \epsilon) \neq 0$)、OLS推定量 $\hat{\beta}_1$ は $\beta_1$ の一致推定量になりません。この状態を内生性(endogeneity)と呼びます。
内生性が生じる主な原因は以下の3つです。
- 未観測の交絡: 省略された変数が $T$ と $Y$ の両方に影響する
- 逆因果: $Y$ が $T$ に影響する
- 測定誤差: $T$ が誤差を含んで観測される
操作変数法は、これらの内生性問題に対処するための手法です。
操作変数のアイデア
操作変数法の直感的なアイデアを理解するために、具体例を考えましょう。
問い: 教育年数 $T$ が賃金 $Y$ に与える因果効果 $\beta_1$ を推定したい。
問題: 能力 $U$(未観測)が教育年数にも賃金にも影響するため、$T$ と $\epsilon$ が相関している。能力が高い人ほど長く教育を受け、かつ能力自体が高賃金に直結する。
解決策: ベトナム戦争時の徴兵くじ(draft lottery)を操作変数 $Z$ として使う(Angristの研究)。くじに当たった若者は兵役に就くため教育が中断され、教育年数が減少する。くじの番号はランダムに決まるため、能力とは無関係。
この例での操作変数 $Z$(くじ番号)の役割は以下のとおりです。
- $Z$ は $T$(教育年数)に影響する(くじに当たると教育が中断される)
- $Z$ は $Y$(賃金)に直接は影響しない(くじ番号自体は賃金に無関係。影響するとすれば教育年数を通じてのみ)
- $Z$ は未観測の交絡 $U$(能力)と無関係(くじはランダム)
このような性質を持つ変数を使えば、$T$ の変動のうち $Z$ に起因する「外生的な」部分のみを取り出し、その部分と $Y$ の関係から因果効果を推定できるのです。
操作変数の直感: 自然実験
操作変数法の魅力は、未観測の交絡が存在する場合でも因果効果を推定できることですが、その代価として「良い操作変数」を見つける必要があります。良い操作変数は、しばしば自然実験(natural experiment)として解釈できます。
自然実験とは、研究者が意図的にランダム化を行ったわけではないが、制度的・自然的な要因により処置の割り当てが「ほぼランダム」になっている状況のことです。徴兵くじはまさにその例であり、くじの番号は個人の特性とは無関係にランダムに割り当てられます。
他の有名な操作変数の例を挙げましょう。Angrist and Krueger(1991)は、生まれた四半期を操作変数として教育年数の効果を推定しました。義務教育制度のもとでは、誕生日の四半期によって退学可能な年齢に達する時期が異なるため、早生まれの生徒は遅生まれの生徒よりも若干短い教育を受ける傾向があります。生まれた四半期はランダムに決まるため、能力とは無関係であると考えられます。
このように、操作変数法は「自然が提供してくれた準ランダム化」を利用する手法であり、適切な自然実験が存在する場面で最も説得力があります。
では、操作変数の条件を正式に定義しましょう。
操作変数の定義と3条件
操作変数の3条件
変数 $Z$ が処置 $T$ のアウトカム $Y$ に対する因果効果を推定するための操作変数(instrumental variable)であるとは、以下の3つの条件を同時に満たすことです。
条件1: 関連性(relevance)
$$ \begin{equation} \text{Cov}(Z, T) \neq 0 \end{equation} $$
操作変数 $Z$ は処置 $T$ と相関がなければなりません。$Z$ が $T$ に影響しなければ、$T$ の外生的な変動を取り出すことができません。
条件2: 排除制約(exclusion restriction)
$$ \begin{equation} Z \text{ は } T \text{ を通じてのみ } Y \text{ に影響する} \end{equation} $$
操作変数 $Z$ はアウトカム $Y$ に直接的な効果を持ってはいけません。$Z$ の $Y$ への影響は、すべて $T$ を経由するものでなければなりません。構造方程式で書けば、$Y = \beta_0 + \beta_1 T + \epsilon$ のモデルに $Z$ が直接現れないことを意味します。
条件3: 独立性(independence / exogeneity)
$$ \begin{equation} \text{Cov}(Z, \epsilon) = 0 \end{equation} $$
操作変数 $Z$ はモデルの誤差項 $\epsilon$(未観測の交絡変数を含む)と無相関でなければなりません。
DAGによる表現
3つの条件をDAGで表すと次のようになります。
$$ Z \to T \to Y, \quad U \to T, \quad U \to Y $$
ここで $U$ は未観測の交絡変数です。重要な点は以下のとおりです。
- $Z \to T$ のパスがある(条件1: 関連性)
- $Z \to Y$ の直接パスがない(条件2: 排除制約)
- $Z$ と $U$ の間にパスがない(条件3: 独立性)
条件の検証可能性
3つの条件のうち、条件1(関連性)のみがデータから検証可能です。$Z$ と $T$ の相関はデータから直接計算できます。
条件2(排除制約)と条件3(独立性)はデータからは検証できません。これらはドメイン知識と理論的な議論に基づいて正当化する必要があります。操作変数法の分析では、これらの仮定の妥当性を丁寧に議論することが不可欠です。
排除制約の議論の重要性
排除制約はIV分析における最も議論の余地がある仮定です。排除制約が成り立たない例をいくつか考えてみましょう。
教育年数の効果を推定するために徴兵くじを操作変数として使う場合、くじに当たって兵役に就いた経験そのものが賃金に影響する可能性があります(退役軍人優遇制度、軍隊でのスキル獲得、PTSDによる労働能力の低下など)。これらの経路が存在すれば、くじ番号は教育年数を通じてだけでなく、兵役経験を通じても直接的に賃金に影響するため、排除制約が破れます。
このように、排除制約の妥当性は常に議論の余地があり、研究者はこの仮定に対する脅威を丁寧に列挙し、それぞれについて反論を提示する必要があります。完全に反論できない脅威がある場合は、感度分析によって排除制約からの乖離がどの程度まで結果に影響するかを評価することが推奨されます。
操作変数の具体例
実証研究で用いられてきた代表的な操作変数を紹介しましょう。
地理的近接性を操作変数として使う研究は多くあります。Card(1993)は、大学までの距離を操作変数として教育の収益率を推定しました。大学に近い地域に住む若者は大学に進学しやすいが、大学の近さ自体は賃金に直接影響しないという仮定です。
気象条件も操作変数として頻繁に使われます。Miguel, Satyanath, and Sergenti(2004)は、降水量を操作変数としてアフリカにおける経済成長と内戦の関係を分析しました。降水量は農業生産を通じて経済成長に影響するが、内戦を直接引き起こすのは経済成長の低下を通じてのみであるという仮定です。
メンデリアンランダム化(Mendelian randomization)は、遺伝子変異を操作変数として用いる疫学的手法です。遺伝子変異は受精時にランダムに割り当てられるため、環境要因とは独立であり(条件3を満たしやすい)、特定の表現型(例えば飲酒量)に影響を与えるが、アウトカム(例えば心疾患)に直接は影響しないという仮定のもとで用いられます。
条件が定義できたところで、最もシンプルなケースであるワルド推定量を導出しましょう。
ワルド推定量 — 最もシンプルなIV推定
二値の操作変数の場合
操作変数 $Z$ が二値($Z \in \{0, 1\}$)で、処置 $T$ も二値($T \in \{0, 1\}$)の場合、因果効果の推定量としてワルド推定量(Wald estimator)が使えます。
$$ \begin{equation} \hat{\beta}_{\text{Wald}} = \frac{E[Y \mid Z=1] – E[Y \mid Z=0]}{E[T \mid Z=1] – E[T \mid Z=0]} = \frac{\text{$Z$の$Y$への誘導効果}}{\text{$Z$の$T$への効果}} \end{equation} $$
ワルド推定量の直感
ワルド推定量の分子 $E[Y \mid Z=1] – E[Y \mid Z=0]$ は、$Z$ の値が変わったときのアウトカム $Y$ の変化量です。排除制約(条件2)により、$Z$ は $T$ を通じてのみ $Y$ に影響するため、この変化は「$Z$ が $T$ を変え、その変わった $T$ が $Y$ を変えた」ことによるものです。
分母 $E[T \mid Z=1] – E[T \mid Z=0]$ は、$Z$ の値が変わったときの処置 $T$ の変化量です。
したがって、ワルド推定量は「$Z$ による $Y$ の変化を、$Z$ による $T$ の変化で割る」ことで、$T$ が1単位変わったときの $Y$ の変化、すなわち因果効果 $\beta_1$ を推定しているのです。
ワルド推定量の導出
式(1)のモデル $Y = \beta_0 + \beta_1 T + \epsilon$ のもとで、条件3($\text{Cov}(Z, \epsilon) = 0$)を使ってワルド推定量が $\beta_1$ を推定することを確認します。
$Z$ が二値のとき、$\text{Cov}(Z, Y)$ を計算すると次のようになります。
$$ \text{Cov}(Z, Y) = \text{Cov}(Z, \beta_0 + \beta_1 T + \epsilon) = \beta_1 \text{Cov}(Z, T) + \text{Cov}(Z, \epsilon) $$
条件3より $\text{Cov}(Z, \epsilon) = 0$ なので次のようになります。
$$ \text{Cov}(Z, Y) = \beta_1 \text{Cov}(Z, T) $$
したがって次の式が得られます。
$$ \beta_1 = \frac{\text{Cov}(Z, Y)}{\text{Cov}(Z, T)} $$
二値の $Z$ に対しては、$\text{Cov}(Z, Y) = P(Z=1)(1-P(Z=1))(E[Y|Z=1] – E[Y|Z=0])$ であり、$\text{Cov}(Z, T)$ も同様の形なので、$P(Z=1)(1-P(Z=1))$ が打ち消し合い、式(5)のワルド推定量が得られます。
連続的な操作変数の場合に対応するのが2段階最小二乗法です。
2段階最小二乗法(2SLS)
2SLSの手順
2段階最小二乗法(Two-Stage Least Squares, 2SLS)は、操作変数法の最も標準的な推定方法です。名前のとおり、2つの段階から構成されます。
第1段階(first stage): 処置 $T$ を操作変数 $Z$(と外生的な共変量 $\bm{X}$)に回帰し、$T$ の予測値 $\hat{T}$ を求める。
$$ \begin{equation} T_i = \pi_0 + \pi_1 Z_i + \bm{\pi}_X^\top \bm{X}_i + \nu_i \end{equation} $$
$$ \hat{T}_i = \hat{\pi}_0 + \hat{\pi}_1 Z_i + \hat{\bm{\pi}}_X^\top \bm{X}_i $$
第2段階(second stage): アウトカム $Y$ を $\hat{T}$(と外生的な共変量 $\bm{X}$)に回帰し、$\hat{T}$ の係数を因果効果の推定値とする。
$$ \begin{equation} Y_i = \beta_0 + \beta_1 \hat{T}_i + \bm{\beta}_X^\top \bm{X}_i + \text{error}_i \end{equation} $$
なぜ2SLSが機能するのか
2SLSが因果効果を正しく推定できる理由を直感的に理解しましょう。
第1段階で $T$ を $Z$ に回帰すると、$T$ は「$Z$ で説明される部分」$\hat{T}$ と「$Z$ で説明されない部分」$\hat{\nu}$ に分解されます。
$$ T = \hat{T} + \hat{\nu} $$
$\hat{T}$ は $Z$ の線形結合であり、条件3($\text{Cov}(Z, \epsilon) = 0$)により $\hat{T}$ は $\epsilon$ と無相関です。つまり、$\hat{T}$ は $T$ の変動のうち「外生的な部分」のみを取り出したものです。
$\hat{\nu}$ は $Z$ では説明できない $T$ の変動であり、ここに未観測交絡の影響が含まれています。
第2段階で $Y$ を $\hat{T}$ に回帰することで、$T$ の外生的な変動のみを使って $Y$ との関係を推定しているため、内生性によるバイアスが除去されるのです。
この直感をさらに噛み砕きましょう。元の $T$ は「$Z$ による外生的な変動」と「$U$ による内生的な変動」の混合物です。OLSが $T$ をそのまま使うと、$T$ と $\epsilon$ の相関($U$ を通じた)がバイアスを引き起こします。2SLSの第1段階は、この混合物から $Z$ に由来する「きれいな」部分のみを抽出するフィルターの役割を果たしています。第2段階はこの「きれいな」$\hat{T}$ と $Y$ の関係を推定するだけであり、$U$ の汚染が入り込む余地はありません。
2SLSの標準誤差
2SLSの標準誤差の計算には注意が必要です。第2段階の回帰で通常の標準誤差を計算すると、$\hat{T}$ が推定値であることの不確実性が無視されるため、正しい標準誤差が得られません。
正しい標準誤差は、以下のように計算されます。第2段階の残差は $\hat{e}_i = Y_i – \hat{\beta}_0 – \hat{\beta}_1 T_i$($\hat{T}_i$ ではなく元の $T_i$ を使う)で計算し、これを用いて分散共分散行列を構成します。
$$ \widehat{\text{Var}}(\hat{\beta}_{2SLS}) = \hat{\sigma}^2 (\hat{\bm{T}}^\top \hat{\bm{T}})^{-1} $$
ここで $\hat{\sigma}^2 = \frac{1}{n-k}\sum_i \hat{e}_i^2$ です。多くの統計ソフトウェア(Stata、R、Pythonのlinearmodels等)は2SLSの正しい標準誤差を自動的に計算しますが、手動で2段階に分けて推定する場合はこの修正が必要であることに注意してください。
2SLSの一致性の証明
2SLS推定量が $\beta_1$ の一致推定量であることを代数的に示しましょう。
まず、IV推定量は次のように書けます。$\bm{X}$ を省略して簡略化します。
$$ \hat{\beta}_{1,\text{IV}} = \frac{\sum_i (Z_i – \bar{Z})(Y_i – \bar{Y})}{\sum_i (Z_i – \bar{Z})(T_i – \bar{T})} $$
これを確率極限として評価するために、$Y_i = \beta_0 + \beta_1 T_i + \epsilon_i$ を代入すると次のようになります。
$$ \hat{\beta}_{1,\text{IV}} = \frac{\sum_i (Z_i – \bar{Z})(\beta_1 T_i + \epsilon_i – \beta_1 \bar{T} – \bar{\epsilon})}{\sum_i (Z_i – \bar{Z})(T_i – \bar{T})} $$
定数項の差は分子では打ち消し合うので、次のように整理できます。
$$ \hat{\beta}_{1,\text{IV}} = \beta_1 + \frac{\frac{1}{n}\sum_i (Z_i – \bar{Z})(\epsilon_i – \bar{\epsilon})}{\frac{1}{n}\sum_i (Z_i – \bar{Z})(T_i – \bar{T})} $$
大数の法則により、分子 $\to \text{Cov}(Z, \epsilon) = 0$(条件3)、分母 $\to \text{Cov}(Z, T) \neq 0$(条件1)であるため、次の結果が得られます。
$$ \hat{\beta}_{1,\text{IV}} \xrightarrow{p} \beta_1 + \frac{0}{\text{Cov}(Z, T)} = \beta_1 $$
すなわち、IV推定量は $\beta_1$ の一致推定量です。
弱い操作変数の問題
条件1(関連性)は $\text{Cov}(Z, T) \neq 0$ を要求しますが、この相関が弱い場合(弱い操作変数、weak instrument)、深刻な問題が生じます。
IV推定量の有限サンプルのバイアスは次のように近似されます。
$$ \text{Bias}(\hat{\beta}_{\text{IV}}) \approx \frac{\text{Cov}(Z, \epsilon)}{\text{Cov}(Z, T)} \cdot \frac{1}{F} $$
ここで $F$ は第1段階の $F$ 統計量です。条件3が厳密に成り立てば($\text{Cov}(Z, \epsilon) = 0$)バイアスはゼロですが、実際には$Z$ と $\epsilon$ のわずかな相関(条件3の近似的な成立)がある場合、$F$ が小さいとバイアスが拡大します。
経験的なルールとして、第1段階の $F$ 統計量が10未満であれば弱い操作変数の問題が深刻であると見なされます(Staiger and Stock, 1997の経験則)。
弱い操作変数の問題は、2SLSの推定量の分布が正規分布から大きくずれること、信頼区間のカバレッジが名目水準を大きく下回ること、推定量がOLSのバイアスの方向に引き寄せられることなど、多くの統計的問題を引き起こします。
弱い操作変数への対処法
弱い操作変数が疑われる場合の対処法としては、Anderson-Rubin(AR)検定と条件付き尤度比(CLR)検定が推奨されます。これらは操作変数の強さに関わらず正しいサイズを維持する検定であり、通常のワルド型の信頼区間が弱い操作変数のもとで崩壊する問題を回避できます。
また、LIML推定量(Limited Information Maximum Likelihood)は2SLSと比べて弱い操作変数に対するバイアスが小さいことが知られています。LIMLは、2SLSと同じく$Z$による$T$の外生的変動を利用しますが、尤度に基づく推定であるためメジアンバイアスが小さく、弱い操作変数の問題に対する頑健性が高い推定量です。
LATE — 操作変数法が推定するもの
遵守者の概念
操作変数法の因果的な解釈を深めるために、Imbens and Angrist(1994)が導入したLATE(Local Average Treatment Effect、局所的平均処置効果)の枠組みを紹介します。
操作変数 $Z$ と処置 $T$ がともに二値の場合、各個体は $Z$ の値に対する反応パターンにより4つのタイプに分類できます。
| タイプ | $T(Z=0)$ | $T(Z=1)$ | 説明 |
|---|---|---|---|
| 遵守者(complier) | 0 | 1 | $Z$に従って処置を受ける/受けない |
| 拒否者(defier) | 1 | 0 | $Z$と逆の行動を取る |
| 常に処置(always-taker) | 1 | 1 | $Z$に関わらず常に処置を受ける |
| 常に未処置(never-taker) | 0 | 0 | $Z$に関わらず常に処置を受けない |
単調性の仮定
LATE解釈には、操作変数の3条件に加えて単調性(monotonicity)の仮定が必要です。
$$ T_i(Z=1) \geq T_i(Z=0) \quad \text{for all } i $$
これは「操作変数 $Z$ が1になると、処置を受ける確率が増えることはあっても減ることはない」という仮定であり、拒否者(defier)が存在しないことを意味します。
LATEの定義
単調性のもとで、ワルド推定量は遵守者(complier)に対する平均処置効果を推定します。
$$ \begin{equation} \beta_{\text{Wald}} = \frac{E[Y|Z=1] – E[Y|Z=0]}{E[T|Z=1] – E[T|Z=0]} = E[Y(1) – Y(0) \mid \text{complier}] = \text{LATE} \end{equation} $$
LATEはATEとは一般に異なります。LATEは「操作変数の変化によって処置の受否が変わる個体(遵守者)にとっての」平均処置効果です。
LATEがATEと一致するのは、処置効果が均質な場合(すべての個体で同じ効果を持つ場合)のみです。効果に異質性がある場合、LATEとATEは異なる値を取りうるため、推定結果の解釈には注意が必要です。
LATEの外的妥当性
LATEの解釈における重要な限界は、推定された効果が遵守者のみに適用されることです。遵守者は操作変数の値によって処置の受否が変わるサブグループであり、常受容者や常拒否者とは異なる特性を持つ可能性があります。
例えば、徴兵くじの研究では、くじに当たることで教育が中断される遵守者は「くじの結果がなければ教育を続けたが、兵役義務によって中断せざるを得なかった」人々です。教育への投資意欲が非常に高い人(くじに関係なく教育を続ける常受容者)や、教育への関心が低い人(くじに関係なく教育を中断する常拒否者)とは異なる特性を持つでしょう。したがって、教育年数の賃金への効果は遵守者とそれ以外のグループで異なりうるのです。
LATEの外的妥当性を高めるためには、遵守者の特性を可能な限り記述し、遵守者のサブグループが政策対象として関心のある集団とどの程度重なるかを議論することが重要です。
過剰識別制約の検定
複数の操作変数が利用可能な場合($q > p$、操作変数の数が内生変数の数より多い場合)、過剰識別制約(overidentification)の検定が可能になります。
サーガン検定(Sargan test)/ ハンセンJ検定(Hansen J test)は、すべての操作変数が排除制約を満たしているかを検定します。帰無仮説は「すべての操作変数が妥当」であり、棄却は少なくとも1つの操作変数が排除制約を違反していることを示唆します。
$$ J = n \cdot \hat{\epsilon}_{2SLS}^\top \bm{Z}(\bm{Z}^\top \bm{Z})^{-1}\bm{Z}^\top \hat{\epsilon}_{2SLS} \xrightarrow{d} \chi^2(q – p) $$
ここで $\hat{\epsilon}_{2SLS}$ は2SLSの残差です。$J$ 統計量が有意であれば、操作変数の一部が妥当でない可能性があり、モデルの再検討が必要です。ただし、この検定は少なくとも1つの操作変数が妥当であることを前提としており、すべての操作変数が妥当でない場合にはこの検定自体が信頼できない点に注意が必要です。
Pythonシミュレーションで2SLSの実装と性質を確認してみましょう。
Pythonシミュレーション
シミュレーション1: 2SLSの実装
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
np.random.seed(42)
n = 5000
# --- データ生成 ---
# 未観測の交絡変数
U = np.random.normal(0, 1, n)
# 操作変数(Uと独立)
Z = np.random.normal(0, 1, n)
# 処置(ZとUの両方に依存)
T = 0.5 + 0.8 * Z + 1.2 * U + np.random.normal(0, 0.5, n)
# アウトカム(Tの因果効果 + Uの直接効果)
beta_true = 2.0
Y = 1.0 + beta_true * T + 1.5 * U + np.random.normal(0, 1, n)
# --- OLS推定(バイアスあり)---
ols = LinearRegression().fit(T.reshape(-1, 1), Y)
beta_ols = ols.coef_[0]
# --- 2SLS推定 ---
# 第1段階
first_stage = LinearRegression().fit(Z.reshape(-1, 1), T)
T_hat = first_stage.predict(Z.reshape(-1, 1))
pi_1 = first_stage.coef_[0]
# 第1段階のF統計量
T_residual = T - T_hat
SSR_first = np.sum(T_residual**2)
SSM_first = np.sum((T_hat - T.mean())**2)
F_stat = SSM_first / (SSR_first / (n - 2))
# 第2段階
second_stage = LinearRegression().fit(T_hat.reshape(-1, 1), Y)
beta_2sls = second_stage.coef_[0]
# ワルド推定量(連続版 = IV推定量)
beta_iv = np.cov(Z, Y)[0, 1] / np.cov(Z, T)[0, 1]
print(f"真の因果効果: {beta_true:.4f}")
print(f"OLS推定値: {beta_ols:.4f} (バイアス: {beta_ols - beta_true:+.4f})")
print(f"2SLS推定値: {beta_2sls:.4f} (バイアス: {beta_2sls - beta_true:+.4f})")
print(f"IV推定値: {beta_iv:.4f} (バイアス: {beta_iv - beta_true:+.4f})")
print(f"第1段階の係数: {pi_1:.4f}")
print(f"第1段階のF統計量: {F_stat:.1f}")
# --- 可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a) 第1段階: ZとTの関係
ax = axes[0]
ax.scatter(Z, T, alpha=0.1, s=5, color="steelblue")
z_sort = np.sort(Z)
ax.plot(z_sort, first_stage.predict(z_sort.reshape(-1, 1)),
"r-", linewidth=2, label=rf"$\hat{{\pi}}_1 = {pi_1:.3f}$")
ax.set_xlabel("Z (instrument)", fontsize=12)
ax.set_ylabel("T (treatment)", fontsize=12)
ax.set_title(f"(a) First stage (F = {F_stat:.1f})", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# (b) OLS vs IV
ax = axes[1]
ax.scatter(T, Y, alpha=0.05, s=5, color="gray")
t_sort = np.sort(T)
ax.plot(t_sort, ols.predict(t_sort.reshape(-1, 1)),
"b-", linewidth=2, label=f"OLS: {beta_ols:.3f}")
ax.plot(t_sort, second_stage.predict(
first_stage.predict(
np.polyfit(Z, t_sort, 0).reshape(-1, 1) if False else
np.full((len(t_sort), 1), 0)
)),
"r--", linewidth=2)
ax.plot(t_sort, beta_true * t_sort + (Y.mean() - beta_true * T.mean()),
"g--", linewidth=2, label=f"True: {beta_true:.3f}")
ax.set_xlabel("T (treatment)", fontsize=12)
ax.set_ylabel("Y (outcome)", fontsize=12)
ax.set_title("(b) OLS vs True causal effect", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# (c) 推定結果の比較
ax = axes[2]
methods = ["True", "OLS", "2SLS"]
values = [beta_true, beta_ols, beta_2sls]
colors = ["green", "red", "blue"]
bars = ax.bar(methods, values, color=colors, alpha=0.7, edgecolor="black")
ax.axhline(beta_true, color="green", linestyle="--", linewidth=1.5)
for bar, val in zip(bars, values):
ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.03,
f"{val:.3f}", ha="center", fontsize=11)
ax.set_ylabel("Estimated coefficient", fontsize=12)
ax.set_title("(c) Comparison", fontsize=12)
ax.grid(True, alpha=0.3, axis="y")
plt.tight_layout()
plt.savefig("iv_2sls_basic.png", dpi=150, bbox_inches="tight")
plt.show()
2SLSの結果から、以下の重要な点が確認できます。
-
OLS推定値は真の因果効果を過大に推定しています。未観測の交絡 $U$ が $T$ と $Y$ の両方に正の影響を与えるため、$T$ と $\epsilon$ の正の相関が生じ、上方バイアスが発生しています
-
2SLS推定値は真の因果効果にほぼ一致しています。操作変数 $Z$ が $U$ と独立であるため、$T$ の変動のうち $Z$ に起因する外生的な部分のみを使って因果効果を正しく推定できています
-
第1段階の $F$ 統計量が大きい。$F \gg 10$ であり、弱い操作変数の問題はありません。操作変数 $Z$ が処置 $T$ と十分に強い関連を持っています
シミュレーション2: 弱い操作変数の影響
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
np.random.seed(0)
n_sims = 1000
n = 500
beta_true = 2.0
# 操作変数の強さを変化させる
pi_values = [0.05, 0.1, 0.2, 0.5, 0.8, 1.5]
fig, axes = plt.subplots(2, 3, figsize=(15, 9))
for idx, pi_1 in enumerate(pi_values):
ols_list = []
iv_list = []
F_list = []
for _ in range(n_sims):
U = np.random.normal(0, 1, n)
Z = np.random.normal(0, 1, n)
T = 0.5 + pi_1 * Z + 1.0 * U + np.random.normal(0, 1, n)
Y = 1.0 + beta_true * T + 1.5 * U + np.random.normal(0, 1, n)
# OLS
ols = LinearRegression().fit(T.reshape(-1, 1), Y)
ols_list.append(ols.coef_[0])
# IV
cov_zy = np.cov(Z, Y)[0, 1]
cov_zt = np.cov(Z, T)[0, 1]
if abs(cov_zt) > 1e-10:
iv_list.append(cov_zy / cov_zt)
# F統計量
fs = LinearRegression().fit(Z.reshape(-1, 1), T)
T_hat = fs.predict(Z.reshape(-1, 1))
SSR = np.sum((T - T_hat)**2)
SSM = np.sum((T_hat - T.mean())**2)
F_list.append(SSM / (SSR / (n-2)))
iv_arr = np.array(iv_list)
mean_F = np.mean(F_list)
ax = axes[idx // 3, idx % 3]
ax.hist(ols_list, bins=40, alpha=0.4, color="red", density=True, label="OLS")
iv_trimmed = iv_arr[(iv_arr > -5) & (iv_arr < 10)]
ax.hist(iv_trimmed, bins=40, alpha=0.4, color="blue", density=True, label="IV")
ax.axvline(beta_true, color="black", linewidth=2, linestyle="--")
ax.set_title(rf"$\pi_1 = {pi_1}$, mean F = {mean_F:.1f}", fontsize=11)
ax.set_xlabel("Estimate", fontsize=10)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_xlim(-3, 8)
plt.suptitle("Weak vs Strong instruments", fontsize=14, y=1.01)
plt.tight_layout()
plt.savefig("iv_weak_instruments.png", dpi=150, bbox_inches="tight")
plt.show()
弱い操作変数のシミュレーション結果から、以下の重要な知見が得られます。
-
$\pi_1$ が小さいとき(弱い操作変数)、IV推定量の分布が極めて広がっている。$\pi_1 = 0.05$ や $0.1$ の場合、IV推定量の分布は真の値の周辺に集中しておらず、非常に不安定です。これは分母 $\text{Cov}(Z, T)$ がゼロに近いため、わずかな誤差が推定値を大きく変動させることの反映です
-
弱い操作変数の場合、IV推定量のバイアスがOLSのバイアスの方向に引き寄せられる。$\pi_1$ が小さいほど、IV推定量の中央値がOLSの値(真の値よりも大きい)に近づく傾向があります
-
$\pi_1$ が大きくなるにつれて、IV推定量の分布が急速に改善。$\pi_1 = 0.5$ 以上では分布が真の値を中心に狭く集中し、OLSとの差が明確になります。$F$ 統計量も10を大きく超えています
-
OLS推定量は一貫して上方バイアスを持ち、$\pi_1$ の値に依存しません。OLSのバイアスは操作変数とは無関係であり、未観測交絡の構造にのみ依存します
操作変数法と他の因果推定手法の比較
操作変数法の位置づけを明確にするために、他の因果推定手法と比較してみましょう。
回帰調整や傾向スコア法は無視可能性(すべての交絡が観測されている)を仮定します。この仮定が成り立てば効率的な推定が可能ですが、未観測交絡が存在する場合にはバイアスが残ります。
操作変数法は無視可能性を仮定しません。代わりに、「良い操作変数が存在する」という異なるタイプの仮定を置きます。この仮定のもとでは未観測交絡が存在しても一致推定が可能ですが、推定されるのはLATEであり、ATEとは一般に異なります。
差の差法(DID)は平行トレンド仮定を置き、時間不変の未観測交絡を除去します。操作変数法とDIDは、どちらも未観測交絡に対処できる手法ですが、必要な仮定と推定対象が異なります。
回帰不連続デザイン(RDD)はカットオフにおける連続性仮定を置き、カットオフ近傍での局所的な因果効果を推定します。ファジーRDDは操作変数法の特殊ケースとして解釈でき、カットオフの上下を操作変数として用いる2SLS推定に帰着します。
実務的には、利用可能なデータと妥当な仮定に基づいて手法を選択し、複数の手法の結果を比較することでロバスト性を確認することが推奨されます。
操作変数法の実践的チェックリスト
操作変数法を実際に適用する際に確認すべき項目を整理します。
デザイン段階
まず、操作変数の候補がなぜ3条件を満たすと考えられるかを理論的に議論します。特に排除制約と独立性はデータから検証できないため、ドメイン知識に基づいた丁寧な議論が不可欠です。排除制約が破れうる経路をすべて列挙し、それぞれについて反論を提示することが望ましいです。
推定段階
第1段階の回帰結果を報告し、操作変数と処置の関連を確認します。第1段階のF統計量が10以上であることを確認します。F統計量が10未満の場合は弱い操作変数の問題が疑われるため、LIML推定量やAnderson-Rubin検定の使用を検討します。
複数の操作変数がある場合は、過剰識別制約のサーガン検定を実施して、操作変数の妥当性を間接的に検証します。ただし、この検定は少なくとも一つの操作変数が妥当であることを前提としている点に注意が必要です。
報告段階
OLS推定量と2SLS推定量を並べて報告し、推定量の差が内生性の存在を示唆するかどうかを議論します。ハウスマン検定(Hausman test)はOLSとIVの推定量の差が統計的に有意かを検定する方法であり、内生性の有無の判断材料になります。
推定された効果がLATEであることを明示し、遵守者がどのようなサブグループであるかを可能な限り記述します。LATEをATEとして解釈するための追加的な仮定が必要かどうかを議論します。
まとめ
本記事では、操作変数法の理論を体系的に解説しました。
- 操作変数は、関連性、排除制約、独立性の3条件を満たす変数であり、未観測の交絡が存在する場合でも因果効果の推定を可能にする
- ワルド推定量は二値の操作変数に対する最もシンプルなIV推定量であり、誘導効果を処置の変化量で割った値
- 2段階最小二乗法(2SLS) は一般的なIV推定の標準的手法。第1段階で処置の外生的変動を取り出し、第2段階でその変動とアウトカムの関係を推定する
- LATEは操作変数法が推定する因果効果の正確な解釈であり、遵守者に対する平均処置効果
- 弱い操作変数は推定の不安定性と有限サンプルバイアスを引き起こすため、第1段階の $F$ 統計量による診断が重要
次のステップとして、以下の記事も参考にしてください。