問いと設定
観測していない共通原因があると、通常のbackdoor調整では因果効果を識別できないことがあります。操作変数（IV）は、別の仮定を使う設計です。「未観測交絡をソフトウェアで消す」方法ではありません。
Zが介入Tに影響するrelevance、Zと結果側の未観測要因の独立性、Tを経由しないZ→Yの効果がないexclusionを検討します。二値Z・Tで効果が異質な場合、monotonicityなどの条件の下で識別するLATEはcomplierの効果であり、一般の集団ATEではありません。
今回の生成式は効果を全員で1.5に固定します。有効なIV、exclusion違反、弱いIVを作り、同じ推定器がどう変わるかを比較します。
実行できるコード
Python 3.12.14で実行しました。解析はすべて合成データを使う研究・学習用の例です。
専用の仮想環境で、次の版をインストールして実行します。数値は2026-10-02の実行結果です。
python -m pip install numpy==2.4.6 pandas==3.0.6 scipy==1.15.3 scikit-learn==1.6.1 statsmodels==0.15.0 dowhy==0.14
次のコードは上から一括実行できます。seedを固定していますが、依存ライブラリや実装の変更で細かな数値は変わり得ます。
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy.special import expit
from dowhy import CausalModel
rng = np.random.default_rng(20261002)
n = 12000
x, u = rng.normal(size=(2, n))
z = rng.binomial(1, .5, n)
noise = rng.normal(size=n)
graph = '''digraph {X -> T; X -> Y; U -> T; U -> Y; Z -> T; T -> Y;}'''
rows = []
for name, strength, direct in [('valid',1.6,0), ('invalid exclusion',1.6,.8), ('weak',.08,0)]:
t = rng.binomial(1, expit(-.2+strength*z+.8*x+u))
y = 1.5*t+.7*x+1.2*u+direct*z+noise
df = pd.DataFrame(dict(X=x, Z=z, T=t, Y=y))
first = sm.OLS(t, sm.add_constant(df[['Z','X']])).fit()
fitted = first.fittedvalues
# 2SLS point estimate; second-stage OLS SE is NOT a valid IV SE.
second = sm.OLS(y, np.column_stack([np.ones(n), fitted, x])).fit()
ols = sm.OLS(y, sm.add_constant(df[['T','X']])).fit()
model = CausalModel(data=df, treatment='T', outcome='Y', graph=graph)
estimand = model.identify_effect(proceed_when_unidentifiable=True)
iv = model.estimate_effect(estimand, method_name='iv.instrumental_variable', method_params={'iv_instrument_name':'Z'})
rows.append(dict(scenario=name, true_effect=1.5, ols=ols.params['T'],
two_stage_point=second.params[1], dowhy_wald=iv.value,
first_stage_partial_F=float(first.tvalues['Z']**2)))
print(pd.DataFrame(rows).round(4).to_string(index=False))
# Sharp RD: score directly affects Y. Do NOT treat score as an excluded IV.
score = rng.uniform(-1, 1, n)
t = (score >= 0).astype(int)
y = 2*t+.8*score+.5*score**2+rng.normal(size=n)
for h in [.05, .1, .2]:
keep = abs(score) <= h
design = np.column_stack([np.ones(keep.sum()),t[keep],score[keep],t[keep]*score[keep]])
fit = sm.OLS(y[keep], design).fit(cov_type='HC1')
print('RD bandwidth, n, jump, robust SE:', h, keep.sum(), round(fit.params[1],4), round(fit.bse[1],4))
実行結果
| 条件 | OLS | 共変量調整2SLSの点推定 | DoWhyのWald推定 | 第一段階partial F |
|---|---|---|---|---|
| 有効なIV | 2.4816 | 1.5332 | 1.5674 | 1134.28 |
| exclusion違反 | 2.6917 | 4.4192 | 4.4348 | 1138.88 |
| 弱いIV | 2.5413 | 1.8838 | 2.2236 | 7.55 |
真の効果はすべて1.5です。ZはX・Uから独立に生成しています。DoWhyの二値IVのWald計算と、Xを調整した2SLSは有限標本で同じ計算ではありません。
強い第一段階でも無効なIVは救えない
exclusion違反の条件では、第一段階Fが大きくても推定値は約4.4でした。ZからYへの直接効果0.8を加えたためです。コードでは解析者が誤ってZ→Yを省いたDAGを与えています。識別結果が返ることは、その省略が正しいことの証明になりません。
弱いIVでは標本変動が増えます。F<10は古典的な目安の一つですが、10を超えれば常に安全という基準ではありません。この1標本の値からバイアスの向きや安定性は断定できません。反復シミュレーションや弱IVに頑健な推論を別途検討します。コード中の第二段階OLSは点推定の説明用で、その標準誤差をIVの標準誤差として使ってはいけません。
回帰不連続デザインは別の仮定を使う
後半は、score≥0でTが決まるsharp RDです。score自体がYに滑らかに影響するので、scoreをexclusionを満たすIVとみなすグラフは置きません。閾値両側で傾きを分けた局所線形回帰を使います。
| 帯域幅 | 閾値近傍のn | 推定ジャンプ | HC1標準誤差 |
|---|---|---|---|
| 0.05 | 663 | 2.2344 | 0.1649 |
| 0.10 | 1275 | 2.0355 | 0.1150 |
| 0.20 | 2434 | 1.9698 | 0.0815 |
真の閾値での効果は2.0です。狭い帯域幅では標本数が減り、推定は揺れます。このHC1標準誤差は帯域幅選択や局所近似バイアスを補正したRD推論ではありません。
実際のRDでは、潜在結果の閾値近傍での連続性、scoreの操作、他の制度も同じ閾値で変わらないか、帯域幅・関数形への感度を確認します。自然実験という名称だけで設計の妥当性が得られるわけではありません。