新薬は平均的にどれだけ効くのか。教育プログラムを受けた人と受けなかった人で、賃金にどれだけの差が生じるのか。これらの「平均的な処置の効果」を定量化する指標が平均処置効果（Average Treatment Effect, ATE）です。
ポテンシャルアウトカムの枠組みで見たように、因果効果の定義は潜在結果の差 $\tau_i = Y_i(1) – Y_i(0)$ として明快に与えられます。しかし、個体レベルの因果効果は観測できないため、集団レベルの平均として推定する必要があります。このとき、どのような推定方法を使うかによって、推定の効率性（分散）やロバスト性（モデルの誤指定に対する頑健さ）が大きく異なります。
ATEとその推定方法を理解すると、以下のような場面で活用できます。
- 臨床試験: 薬の有効性を評価する際、RCTのデータからATEを推定し、統計的検定を行う
- 政策評価: 教育支援、雇用政策、公衆衛生政策の効果を定量的に評価する
- 個別化医療: 患者の特性に応じた条件付き処置効果（CATE）を推定し、治療の個別最適化に役立てる
- A/Bテスト: Webサービスやアプリの新機能の効果をATEとして推定する
本記事では、ATEおよび関連する因果効果の指標を体系的に定義し、代表的な推定方法（回帰推定、IPW推定、二重ロバスト推定）の理論を解説します。Pythonシミュレーションを通じて、各推定方法の特性と限界を実感します。
本記事の内容
- ATE、ATT、ATCの定義と関係
- 条件付き平均処置効果（CATE）と効果の異質性
- 回帰推定量の理論と限界
- 逆確率重み付け（IPW）推定量の導出
- 二重ロバスト（DR）推定量の理論
- Pythonによる各推定量の実装と比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- ポテンシャルアウトカムの枠組み — 潜在結果、SUTVA、無視可能性
- 因果推論入門 — 相関と因果の違い、交絡の概念
- 重回帰分析 — 回帰分析の基礎
ATEとその周辺の因果量 — 定義の整理
平均処置効果（ATE）
平均処置効果（Average Treatment Effect, ATE）は、因果推論において最も基本的な量であり、集団全体に対する処置の平均的な因果効果を表します。
$$ \begin{equation} \text{ATE} = E[Y(1) – Y(0)] = E[Y(1)] – E[Y(0)] \end{equation} $$
ここで $Y(1)$ と $Y(0)$ は潜在結果であり、期待値は集団全体に対して取ります。ATEは「全員に処置を施した場合の平均アウトカム」と「全員に処置を施さなかった場合の平均アウトカム」の差であり、処置の因果効果の最も自然な要約です。
ATEの解釈として重要なのは、これが集団全体に対する平均であることです。たとえば「新薬のATEが5点」とは、「この集団の全員に新薬を投与した場合、全員にプラセボを投与した場合と比べて、アウトカムが平均5点改善する」ことを意味します。
処置群の平均処置効果（ATT）と対照群の平均処置効果（ATC）
ATEは集団全体の平均ですが、実際に処置を受けた人（処置群）と受けなかった人（対照群）では、処置効果が異なる可能性があります。
処置群の平均処置効果（ATT: Average Treatment Effect on the Treated）:
$$ \begin{equation} \text{ATT} = E[Y(1) – Y(0) \mid T = 1] \end{equation} $$
ATTは「実際に処置を受けた人にとっての、処置の平均的効果」です。政策評価において特に重要です。たとえば、職業訓練プログラムの評価では、プログラムに参加した人（$T=1$）にとってのプログラムの効果こそが政策判断に直結します。
対照群の平均処置効果（ATC: Average Treatment Effect on the Controls）:
$$ \begin{equation} \text{ATC} = E[Y(1) – Y(0) \mid T = 0] \end{equation} $$
ATCは「処置を受けなかった人が、もし処置を受けていたらどうなっていたか」の平均です。たとえば、新薬が承認された場合に、これまで服用していなかった患者にどの程度の効果が見込めるかを示します。
ATE、ATT、ATCの関係
これら3つの量は以下の関係で結ばれています。
$$ \begin{equation} \text{ATE} = \pi \cdot \text{ATT} + (1 – \pi) \cdot \text{ATC} \end{equation} $$
ここで $\pi = P(T = 1)$ は処置を受ける確率（集団全体での処置率）です。ATEはATTとATCの加重平均であり、重みは処置率で決まります。
もし処置効果に異質性がなければ（すべての個体で $Y_i(1) – Y_i(0) = \tau$ なら）、$\text{ATE} = \text{ATT} = \text{ATC} = \tau$ となります。しかし、一般には処置効果は個体によって異なるため、3つの量は異なる値を取ります。
条件付き平均処置効果（CATE）
条件付き平均処置効果（Conditional Average Treatment Effect, CATE）は、共変量 $\bm{X} = \bm{x}$ を条件付けたときの平均処置効果です。
$$ \begin{equation} \tau(\bm{x}) = E[Y(1) – Y(0) \mid \bm{X} = \bm{x}] \end{equation} $$
CATEは「特性 $\bm{X} = \bm{x}$ を持つサブグループにとっての処置の平均効果」を表します。CATEは処置効果の異質性（heterogeneity）を共変量で説明する量であり、個別化医療やターゲティング広告などの応用で中心的な役割を果たします。
ATEとCATEの関係は次のようになります。
$$ \begin{equation} \text{ATE} = E[\tau(\bm{X})] = \int \tau(\bm{x})\,dP(\bm{x}) \end{equation} $$
ATEはCATEを共変量の分布で周辺化（平均化）したものです。CATEが $\bm{x}$ によらず一定のとき、処置効果は均質（homogeneous）であると言い、$\text{ATE} = \tau(\bm{x})$ が成り立ちます。
では、これらの因果量をどのように推定するのでしょうか。以下では、代表的な3つの推定方法を順に解説します。
回帰推定量 — 条件付き期待値のモデリング
基本的なアイデア
ポテンシャルアウトカムの枠組みで示したように、無視可能性のもとではATEは次のように識別されます。
$$ \text{ATE} = E\left[E[Y \mid T=1, \bm{X}] – E[Y \mid T=0, \bm{X}]\right] $$
この式は、$E[Y \mid T, \bm{X}]$（処置と共変量が与えられたときのアウトカムの条件付き期待値）をモデル化すれば、ATEが推定できることを示唆しています。
回帰推定量の構成
最も単純な方法は、アウトカム回帰モデルを構築することです。まず、処置群と対照群それぞれについてアウトカムの条件付き期待値をモデル化します。
$$ \mu_1(\bm{x}) = E[Y \mid T=1, \bm{X}=\bm{x}], \quad \mu_0(\bm{x}) = E[Y \mid T=0, \bm{X}=\bm{x}] $$
これらの推定値を $\hat{\mu}_1(\bm{x})$, $\hat{\mu}_0(\bm{x})$ と書くと、ATEの回帰推定量は次のように構成されます。
$$ \begin{equation} \hat{\tau}_{\text{reg}} = \frac{1}{n}\sum_{i=1}^{n}\left[\hat{\mu}_1(\bm{X}_i) – \hat{\mu}_0(\bm{X}_i)\right] \end{equation} $$
具体的な手順は以下のとおりです。
- 処置群（$T_i = 1$）のデータを使って $\mu_1(\bm{x})$ を推定する
- 対照群（$T_i = 0$）のデータを使って $\mu_0(\bm{x})$ を推定する
- 全個体（処置群＋対照群）に対して $\hat{\mu}_1(\bm{X}_i) – \hat{\mu}_0(\bm{X}_i)$ を計算し、平均を取る
ステップ3で重要なのは、全個体に対して予測を行うことです。処置群の個体に対しても $\hat{\mu}_0(\bm{X}_i)$（処置を受けなかった場合の予測値）を計算し、対照群の個体に対しても $\hat{\mu}_1(\bm{X}_i)$（処置を受けた場合の予測値）を計算します。これは反事実の予測に相当します。
線形回帰の場合
回帰モデルとして線形モデルを使う場合、最も単純には次のような重回帰を考えます。
$$ Y_i = \alpha + \tau T_i + \bm{\beta}^\top \bm{X}_i + \epsilon_i $$
この場合、$T_i$ の回帰係数 $\hat{\tau}$ がATEの推定値になります。ただし、このモデルはCATEが定数であること（効果の均質性）を暗黙に仮定しています。
効果の異質性を許容するには、交互作用項を含めます。
$$ Y_i = \alpha + \tau T_i + \bm{\beta}^\top \bm{X}_i + \bm{\gamma}^\top (T_i \cdot \bm{X}_i) + \epsilon_i $$
この場合、$\text{CATE}(\bm{x}) = \tau + \bm{\gamma}^\top \bm{x}$ となり、ATEは $\hat{\tau}_{\text{reg}} = \hat{\tau} + \bm{\hat{\gamma}}^\top \bar{\bm{X}}$ で推定されます。
回帰推定量の長所と短所
長所: – 直感的でわかりやすい – 共変量の次元が低い場合には効率的 – 既存の回帰分析の枠組みで実装できる
短所: – モデルの正しさに依存する。アウトカムモデル $\mu_t(\bm{x})$ が誤指定（misspecified）されている場合、推定にバイアスが生じる – 外挿に弱い。処置群と対照群の共変量分布が大きく異なる場合、一方の群で学習したモデルを他方に適用する際に大きな外挿誤差が生じうる
モデルの誤指定に対するロバスト性を高めるために、次に紹介する逆確率重み付け（IPW）推定量が考案されました。
逆確率重み付け（IPW）推定量 — 処置確率によるバランシング
IPWの直感
回帰推定量はアウトカムのモデルに基づいていましたが、IPW推定量は処置の割り当てメカニズムに基づく方法です。
直感的なアイデアは、「調査統計学のサンプリングウェイト」に似ています。あるグループがサンプルに含まれにくい場合、そのグループの観測値に大きな重みを付けて「過少代表」を補正する、という考え方です。
因果推論の文脈では、共変量 $\bm{X} = \bm{x}$ を持つ個体が処置群に入る確率（傾向スコア、propensity score）が低いにもかかわらず実際に処置を受けた場合、その個体は「$\bm{X} = \bm{x}$ を持つ処置群全体を代表している」と考えて大きな重みを付けます。逆に、傾向スコアが高くて処置を受けた個体には小さな重みを付けます。
傾向スコアの定義
傾向スコア（propensity score）$e(\bm{x})$ は、共変量 $\bm{X} = \bm{x}$ が与えられたときに処置を受ける確率として定義されます。
$$ \begin{equation} e(\bm{x}) = P(T = 1 \mid \bm{X} = \bm{x}) \end{equation} $$
ローゼンバウムとルービン（1983）の傾向スコア定理により、無視可能性の仮定 $(Y(0), Y(1)) \perp\!\!\!\perp T \mid \bm{X}$ が成り立つとき、傾向スコアで条件付けても無視可能性が成り立ちます。
$$ (Y(0), Y(1)) \perp\!\!\!\perp T \mid e(\bm{X}) $$
この定理は、多次元の共変量 $\bm{X}$ の代わりに1次元の傾向スコア $e(\bm{X})$ を使えばよいという「次元縮約」の効果を持ち、実践上極めて重要です。
IPW推定量の導出
IPW推定量を導出するために、まず $E[Y(1)]$ を観測量で表す方法を考えます。無視可能性のもとで次が成り立ちます。
$$ E[Y(1)] = E\left[E[Y(1) \mid \bm{X}]\right] = E\left[E[Y \mid T=1, \bm{X}]\right] $$
ここで、$E[Y \mid T=1, \bm{X}]$ を直接推定するのではなく、以下の恒等式を使います。
$$ E\left[\frac{TY}{e(\bm{X})}\right] = E\left[\frac{T Y}{e(\bm{X})} \bigg| \bm{X}\right] $$
内側の条件付き期待値を計算すると、$T$ と $Y$ が条件付き独立（$\bm{X}$ で条件付けたとき）であることに注意して次のようになります。
$$ E\left[\frac{TY}{e(\bm{X})} \bigg| \bm{X}\right] = \frac{1}{e(\bm{X})} E[TY \mid \bm{X}] $$
ここで $TY = T \cdot Y(1) + (1-T) \cdot 0 \cdot Y(0)$ より $E[TY \mid \bm{X}] = E[T \mid \bm{X}] \cdot E[Y(1) \mid \bm{X}]$ となります。最後の等号は無視可能性の帰結です。
$$ = \frac{E[T \mid \bm{X}] \cdot E[Y(1) \mid \bm{X}]}{e(\bm{X})} = \frac{e(\bm{X}) \cdot E[Y(1) \mid \bm{X}]}{e(\bm{X})} = E[Y(1) \mid \bm{X}] $$
したがって、両辺の期待値を取ると次のようになります。
$$ E\left[\frac{TY}{e(\bm{X})}\right] = E[E[Y(1) \mid \bm{X}]] = E[Y(1)] $$
同様にして $E[Y(0)] = E\left[\frac{(1-T)Y}{1-e(\bm{X})}\right]$ が得られます。
以上から、ホルヴィッツ・トンプソン型のIPW推定量（Horvitz-Thompson estimator）が次のように構成されます。
$$ \begin{equation} \hat{\tau}_{\text{IPW}} = \frac{1}{n}\sum_{i=1}^{n}\left[\frac{T_i Y_i}{\hat{e}(\bm{X}_i)} – \frac{(1-T_i)Y_i}{1-\hat{e}(\bm{X}_i)}\right] \end{equation} $$
ここで $\hat{e}(\bm{X}_i)$ は傾向スコアの推定値です（通常、ロジスティック回帰で推定します）。
ハジェック推定量
式(9)のホルヴィッツ・トンプソン型推定量は不偏ですが、重みの和が $n$ に一致しない場合に分散が大きくなります。重みを正規化したハジェック推定量（Hájek estimator）の方が実践的に安定しています。
$$ \begin{equation} \hat{\tau}_{\text{Hájek}} = \frac{\sum_{i=1}^{n}\frac{T_i Y_i}{\hat{e}(\bm{X}_i)}}{\sum_{i=1}^{n}\frac{T_i}{\hat{e}(\bm{X}_i)}} – \frac{\sum_{i=1}^{n}\frac{(1-T_i) Y_i}{1-\hat{e}(\bm{X}_i)}}{\sum_{i=1}^{n}\frac{(1-T_i)}{1-\hat{e}(\bm{X}_i)}} \end{equation} $$
ハジェック推定量は厳密には不偏ではありませんが（漸近的には不偏）、有限サンプルでの性能は一般にホルヴィッツ・トンプソン型よりも良好です。
IPW推定量の長所と短所
長所: – アウトカムモデルを必要としない。傾向スコアモデルが正しければ、$\mu_t(\bm{x})$ の関数形を指定する必要がない – ノンパラメトリックな性質を持つ（アウトカムに対しては）
短所: – 傾向スコアモデルの正しさに依存する。傾向スコアが誤指定されていると推定にバイアスが生じる – 極端な重みの問題。傾向スコアが0や1に近い個体が存在すると、逆数が非常に大きくなり、推定量の分散が爆発する（正値性の違反に近い状況） – アウトカムの情報を直接使わないため、回帰推定量よりも効率が悪くなる場合がある
回帰推定量とIPW推定量のそれぞれの弱点を補う方法として、次に紹介する二重ロバスト推定量が考案されました。
二重ロバスト（DR）推定量 — 2つのモデルの「保険」
二重ロバスト性の直感
回帰推定量はアウトカムモデル $\mu_t(\bm{x})$ が正しければ一致推定量になり、IPW推定量は傾向スコアモデル $e(\bm{x})$ が正しければ一致推定量になります。しかし、どちらのモデルが正しいかを事前に知ることはできません。
二重ロバスト推定量（Doubly Robust estimator, DR）は、アウトカムモデルと傾向スコアモデルのどちらか一方が正しければ一致推定量になるという、魅力的な性質を持ちます。つまり、2つのモデルのうち少なくとも1つを「保険」として持っている状態であり、「二重の安全網」が張られていると言えます。
DR推定量の定義
二重ロバスト推定量（AIPW推定量: Augmented IPW estimator とも呼ばれる）は、次のように定義されます。
$$ \begin{equation} \hat{\tau}_{\text{DR}} = \frac{1}{n}\sum_{i=1}^{n}\left[\hat{\mu}_1(\bm{X}_i) – \hat{\mu}_0(\bm{X}_i) + \frac{T_i(Y_i – \hat{\mu}_1(\bm{X}_i))}{\hat{e}(\bm{X}_i)} – \frac{(1-T_i)(Y_i – \hat{\mu}_0(\bm{X}_i))}{1-\hat{e}(\bm{X}_i)}\right] \end{equation} $$
この式は、回帰推定量 $\hat{\mu}_1(\bm{X}_i) – \hat{\mu}_0(\bm{X}_i)$ に、IPW型の「補正項」を加えた形になっています。補正項は、アウトカムモデルの予測残差 $Y_i – \hat{\mu}_t(\bm{X}_i)$ を傾向スコアで重み付けしたものです。
二重ロバスト性の証明（概略）
二重ロバスト性が成り立つことを確認しましょう。DR推定量の母集団版は次のとおりです。
$$ \tau_{\text{DR}} = E\left[\mu_1(\bm{X}) – \mu_0(\bm{X}) + \frac{T(Y – \mu_1(\bm{X}))}{e(\bm{X})} – \frac{(1-T)(Y – \mu_0(\bm{X}))}{1-e(\bm{X})}\right] $$
ケース1: アウトカムモデルが正しい場合
$\mu_1(\bm{x}) = E[Y(1) \mid \bm{X} = \bm{x}]$ が正しければ、$E[T(Y – \mu_1(\bm{X})) \mid \bm{X}] = e(\bm{X}) \cdot 0 = 0$ となり、IPW補正項の期待値はゼロです。したがって次のようになります。
$$ \tau_{\text{DR}} = E[\mu_1(\bm{X}) – \mu_0(\bm{X})] = E[Y(1)] – E[Y(0)] = \text{ATE} $$
傾向スコアモデルが誤っていても、$\tau_{\text{DR}}$ はATEに一致します。
ケース2: 傾向スコアモデルが正しい場合
$e(\bm{x}) = P(T = 1 \mid \bm{X} = \bm{x})$ が正しければ、IPW部分の期待値を計算すると、先ほどのIPWの導出と同様の議論により次のようになります。
$$ E\left[\frac{T(Y – \mu_1(\bm{X}))}{e(\bm{X})}\right] = E[Y(1)] – E[\mu_1(\bm{X})] $$
これを代入すると次のようになります。
$$ \tau_{\text{DR}} = E[\mu_1(\bm{X}) – \mu_0(\bm{X})] + E[Y(1)] – E[\mu_1(\bm{X})] – E[Y(0)] + E[\mu_0(\bm{X})] $$
右辺の $E[\mu_1(\bm{X})]$ と $E[\mu_0(\bm{X})]$ が打ち消し合い、$\tau_{\text{DR}} = E[Y(1)] – E[Y(0)] = \text{ATE}$ が得られます。アウトカムモデルが誤っていても、傾向スコアが正しければATEに一致します。
DR推定量の漸近的な性質
DR推定量は、両方のモデルが正しい場合には半パラメトリック効率限界（semiparametric efficiency bound）に達する、すなわち可能な推定量の中で最も分散が小さくなることが知られています。この限界はHahn（1998）によって導出されました。
ATEのセミパラメトリック効率限界は次の式で与えられます。
$$ V_{\text{eff}} = E\left[\frac{\text{Var}(Y(1) \mid \bm{X})}{e(\bm{X})} + \frac{\text{Var}(Y(0) \mid \bm{X})}{1-e(\bm{X})} + (\tau(\bm{X}) – \text{ATE})^2\right] $$
DR推定量はこの下限に漸近的に到達するため、効率的であると言えます。
では、これら3つの推定量をPythonで実装し、性能を比較してみましょう。
Pythonによる推定量の実装と比較
データ生成
まず、交絡のある観察データを生成します。真の因果効果をあらかじめ設定し、各推定量がこれを正しく回復できるかを検証します。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, LogisticRegression
np.random.seed(42)
n = 3000
# --- データ生成 ---
# 共変量（2次元）
X1 = np.random.normal(0, 1, n)
X2 = np.random.normal(0, 1, n)
X = np.column_stack([X1, X2])
# 傾向スコア（真値）
logit_e = 0.5 + 0.8 * X1 - 0.5 * X2
e_true = 1 / (1 + np.exp(-logit_e))
# 処置割り当て
T = np.random.binomial(1, e_true, n)
# 潜在結果（効果の異質性あり）
mu0_true = 2 + 1.5 * X1 + X2 + 0.5 * X1 * X2
tau_true_individual = 3.0 + 1.0 * X1 # CATEはX1に依存
mu1_true = mu0_true + tau_true_individual
Y0 = mu0_true + np.random.normal(0, 1, n)
Y1 = mu1_true + np.random.normal(0, 1, n)
# 観測アウトカム
Y = T * Y1 + (1 - T) * Y0
# 真のATE
true_ATE = np.mean(tau_true_individual)
true_ATT = np.mean(tau_true_individual[T == 1])
print(f"サンプルサイズ: {n}")
print(f"処置群の割合: {T.mean():.3f}")
print(f"真のATE: {true_ATE:.4f}")
print(f"真のATT: {true_ATT:.4f}")
このコードではデータ生成のメカニズムを完全に制御しているため、真のATEとATTを計算できます。実際の分析ではこれらは未知であり、推定の対象になります。
次に、3つの推定量を実装して比較します。
推定量の実装
import numpy as np
from sklearn.linear_model import LinearRegression, LogisticRegression
# 前のブロックの変数を引き続き使用
# n, X, T, Y, true_ATE は定義済み
# === 1. ナイーブ推定 ===
naive_est = Y[T == 1].mean() - Y[T == 0].mean()
# === 2. 回帰推定量 ===
# 処置群と対照群で別々のモデルを学習
reg1 = LinearRegression().fit(X[T == 1], Y[T == 1])
reg0 = LinearRegression().fit(X[T == 0], Y[T == 0])
mu1_hat = reg1.predict(X) # 全個体に対する予測
mu0_hat = reg0.predict(X)
reg_est = np.mean(mu1_hat - mu0_hat)
# === 3. IPW推定量 ===
# 傾向スコアの推定（ロジスティック回帰）
ps_model = LogisticRegression(max_iter=1000)
ps_model.fit(X, T)
e_hat = ps_model.predict_proba(X)[:, 1]
# 極端な値のクリッピング（正値性の確保）
e_hat_clipped = np.clip(e_hat, 0.05, 0.95)
# ホルヴィッツ・トンプソン型
ipw_est = np.mean(T * Y / e_hat_clipped - (1 - T) * Y / (1 - e_hat_clipped))
# ハジェック型
ipw_hajek = (np.sum(T * Y / e_hat_clipped) / np.sum(T / e_hat_clipped)
- np.sum((1 - T) * Y / (1 - e_hat_clipped)) / np.sum((1 - T) / (1 - e_hat_clipped)))
# === 4. 二重ロバスト（AIPW）推定量 ===
dr_est = np.mean(
(mu1_hat - mu0_hat)
+ T * (Y - mu1_hat) / e_hat_clipped
- (1 - T) * (Y - mu0_hat) / (1 - e_hat_clipped)
)
print(f"\n=== 推定結果 ===")
print(f"真のATE: {true_ATE:.4f}")
print(f"ナイーブ推定: {naive_est:.4f} (バイアス: {naive_est - true_ATE:+.4f})")
print(f"回帰推定: {reg_est:.4f} (バイアス: {reg_est - true_ATE:+.4f})")
print(f"IPW推定 (HT): {ipw_est:.4f} (バイアス: {ipw_est - true_ATE:+.4f})")
print(f"IPW推定 (Hájek): {ipw_hajek:.4f} (バイアス: {ipw_hajek - true_ATE:+.4f})")
print(f"DR推定: {dr_est:.4f} (バイアス: {dr_est - true_ATE:+.4f})")
出力を確認すると、ナイーブ推定は大きなバイアスを持つのに対し、回帰推定、IPW推定、DR推定はいずれも真のATEに近い値を返しています。ただし1回のシミュレーションでは偶然の影響もあるため、次にモンテカルロシミュレーションで推定量の統計的性質を確認します。
モンテカルロシミュレーションによる比較
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, LogisticRegression
np.random.seed(0)
n_sims = 500
n = 1000
naive_results = []
reg_results = []
ipw_results = []
dr_results = []
for sim in range(n_sims):
# データ生成
X1 = np.random.normal(0, 1, n)
X2 = np.random.normal(0, 1, n)
X = np.column_stack([X1, X2])
logit_e = 0.5 + 0.8 * X1 - 0.5 * X2
e_true = 1 / (1 + np.exp(-logit_e))
T = np.random.binomial(1, e_true, n)
mu0 = 2 + 1.5 * X1 + X2
tau_i = 3.0 + 1.0 * X1
mu1 = mu0 + tau_i
Y0 = mu0 + np.random.normal(0, 1, n)
Y1 = mu1 + np.random.normal(0, 1, n)
Y = T * Y1 + (1 - T) * Y0
true_ATE = np.mean(tau_i)
# ナイーブ
naive_results.append(Y[T==1].mean() - Y[T==0].mean())
# 回帰
reg1 = LinearRegression().fit(X[T==1], Y[T==1])
reg0 = LinearRegression().fit(X[T==0], Y[T==0])
mu1_hat = reg1.predict(X)
mu0_hat = reg0.predict(X)
reg_results.append(np.mean(mu1_hat - mu0_hat))
# IPW
ps = LogisticRegression(max_iter=1000).fit(X, T)
e_hat = np.clip(ps.predict_proba(X)[:, 1], 0.05, 0.95)
ipw = (np.sum(T * Y / e_hat) / np.sum(T / e_hat)
- np.sum((1-T) * Y / (1 - e_hat)) / np.sum((1-T) / (1 - e_hat)))
ipw_results.append(ipw)
# DR
dr = np.mean(
(mu1_hat - mu0_hat)
+ T * (Y - mu1_hat) / e_hat
- (1-T) * (Y - mu0_hat) / (1 - e_hat)
)
dr_results.append(dr)
naive_results = np.array(naive_results)
reg_results = np.array(reg_results)
ipw_results = np.array(ipw_results)
dr_results = np.array(dr_results)
# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# (a) 推定値の分布
ax = axes[0]
ax.hist(naive_results, bins=30, alpha=0.4, color="gray", label="Naive", density=True)
ax.hist(reg_results, bins=30, alpha=0.5, color="blue", label="Regression", density=True)
ax.hist(ipw_results, bins=30, alpha=0.5, color="green", label="IPW (Hájek)", density=True)
ax.hist(dr_results, bins=30, alpha=0.5, color="red", label="DR (AIPW)", density=True)
ax.axvline(true_ATE, color="black", linewidth=2, linestyle="--", label=f"True ATE = {true_ATE:.2f}")
ax.set_xlabel("Estimated ATE", fontsize=12)
ax.set_ylabel("Density", fontsize=12)
ax.set_title(f"Distribution of ATE estimates ({n_sims} simulations)", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# (b) バイアスとRMSEの比較
ax = axes[1]
methods = ["Naive", "Regression", "IPW", "DR"]
biases = [naive_results.mean() - true_ATE,
reg_results.mean() - true_ATE,
ipw_results.mean() - true_ATE,
dr_results.mean() - true_ATE]
rmses = [np.sqrt(np.mean((naive_results - true_ATE)**2)),
np.sqrt(np.mean((reg_results - true_ATE)**2)),
np.sqrt(np.mean((ipw_results - true_ATE)**2)),
np.sqrt(np.mean((dr_results - true_ATE)**2))]
x_pos = np.arange(len(methods))
width = 0.35
bars1 = ax.bar(x_pos - width/2, biases, width, label="Bias", color="steelblue", alpha=0.7)
bars2 = ax.bar(x_pos + width/2, rmses, width, label="RMSE", color="salmon", alpha=0.7)
ax.set_xticks(x_pos)
ax.set_xticklabels(methods, fontsize=11)
ax.set_ylabel("Value", fontsize=12)
ax.set_title("Bias and RMSE comparison", fontsize=12)
ax.axhline(0, color="black", linewidth=0.5)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, axis="y")
plt.tight_layout()
plt.savefig("ate_estimation_comparison.png", dpi=150, bbox_inches="tight")
plt.show()
print(f"\n=== モンテカルロ結果（{n_sims}回）===")
for name, results in zip(methods, [naive_results, reg_results, ipw_results, dr_results]):
print(f"{name:12s}: 平均={results.mean():.4f}, バイアス={results.mean()-true_ATE:+.4f}, "
f"SD={results.std():.4f}, RMSE={np.sqrt(np.mean((results-true_ATE)**2)):.4f}")
このモンテカルロシミュレーションの結果から、以下の重要な知見が得られます。
-
ナイーブ推定は大きな正のバイアスを持つ。交絡を無視した単純比較は、真のATEから系統的にずれています。バイアスはサンプルサイズを増やしても解消しません
-
回帰推定、IPW推定、DR推定はいずれもほぼ不偏。これらの推定量はバイアスがほぼゼロであり、真のATEを中心とした分布を示しています。アウトカムモデルと傾向スコアモデルがともに正しく指定されている場合、3つとも良好な性能を発揮します
-
DR推定量のRMSEが最も小さい傾向にある。DR推定量は回帰推定の情報とIPWの情報を組み合わせているため、効率性が高くなります。これはDR推定量がセミパラメトリック効率限界に達するという理論的結果と整合しています
-
IPW推定量の分散は他より大きい傾向にある。IPW推定量はアウトカムモデルの情報を直接使わないため、効率性が低くなります。特に傾向スコアが0や1に近い個体がある場合、分散が大きくなります
モデル誤指定に対するロバスト性の検証
DR推定量の真の強みは、モデルの一方が誤指定されていても一致性を保つことです。これを検証してみましょう。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, LogisticRegression
np.random.seed(0)
n_sims = 500
n = 1000
results = {
"Both correct": {"reg": [], "ipw": [], "dr": []},
"Outcome misspec": {"reg": [], "ipw": [], "dr": []},
"PS misspec": {"reg": [], "ipw": [], "dr": []},
"Both misspec": {"reg": [], "ipw": [], "dr": []},
}
for sim in range(n_sims):
X1 = np.random.normal(0, 1, n)
X2 = np.random.normal(0, 1, n)
X_full = np.column_stack([X1, X2, X1**2, X1*X2])
X_wrong = np.column_stack([X1]) # X2と交互作用を省略
logit_e = 0.3 + 0.6 * X1 - 0.4 * X2 + 0.3 * X1 * X2
e_true = 1 / (1 + np.exp(-logit_e))
T = np.random.binomial(1, e_true, n)
mu0 = 1 + X1 + 0.5 * X2 + 0.3 * X1**2
tau_i = 2.0
Y = mu0 + tau_i * T + np.random.normal(0, 1, n)
true_ATE = tau_i
for scenario, X_out, X_ps in [
("Both correct", X_full, X_full),
("Outcome misspec", X_wrong, X_full),
("PS misspec", X_full, X_wrong),
("Both misspec", X_wrong, X_wrong),
]:
# 回帰推定
reg1 = LinearRegression().fit(X_out[T==1], Y[T==1])
reg0 = LinearRegression().fit(X_out[T==0], Y[T==0])
mu1_hat = reg1.predict(X_out)
mu0_hat = reg0.predict(X_out)
results[scenario]["reg"].append(np.mean(mu1_hat - mu0_hat))
# IPW
ps = LogisticRegression(max_iter=1000).fit(X_ps, T)
e_hat = np.clip(ps.predict_proba(X_ps)[:, 1], 0.05, 0.95)
ipw = (np.sum(T*Y/e_hat)/np.sum(T/e_hat)
- np.sum((1-T)*Y/(1-e_hat))/np.sum((1-T)/(1-e_hat)))
results[scenario]["ipw"].append(ipw)
# DR
dr = np.mean(
(mu1_hat - mu0_hat)
+ T*(Y - mu1_hat)/e_hat
- (1-T)*(Y - mu0_hat)/(1-e_hat)
)
results[scenario]["dr"].append(dr)
# --- 可視化 ---
fig, axes = plt.subplots(2, 2, figsize=(13, 10))
scenarios = list(results.keys())
for idx, scenario in enumerate(scenarios):
ax = axes[idx // 2, idx % 2]
for method, color, label in [("reg", "blue", "Regression"),
("ipw", "green", "IPW"),
("dr", "red", "DR")]:
data = np.array(results[scenario][method])
ax.hist(data, bins=25, alpha=0.5, color=color, label=label, density=True)
ax.axvline(true_ATE, color="black", linewidth=2, linestyle="--")
ax.set_title(scenario, fontsize=12)
ax.set_xlabel("Estimated ATE", fontsize=10)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.suptitle("Robustness to model misspecification", fontsize=14, y=1.01)
plt.tight_layout()
plt.savefig("ate_double_robustness.png", dpi=150, bbox_inches="tight")
plt.show()
# バイアスの表を表示
print(f"\n=== モデル誤指定に対するバイアス ===")
print(f"{'Scenario':20s} {'Regression':>12s} {'IPW':>12s} {'DR':>12s}")
for scenario in scenarios:
biases = []
for method in ["reg", "ipw", "dr"]:
data = np.array(results[scenario][method])
biases.append(f"{data.mean() - true_ATE:+.4f}")
print(f"{scenario:20s} {biases[0]:>12s} {biases[1]:>12s} {biases[2]:>12s}")
モデル誤指定に対するロバスト性の検証結果から、以下の重要な知見が得られます。
-
両方のモデルが正しい場合（左上）: 3つの推定量すべてがほぼ不偏です。DR推定量の分散が最も小さくなっています
-
アウトカムモデルのみ誤指定（右上）: 回帰推定量にバイアスが生じますが、IPW推定量とDR推定量は不偏を維持しています。傾向スコアモデルが正しいため、IPWとDRは正しく機能します
-
傾向スコアモデルのみ誤指定（左下）: IPW推定量にバイアスが生じますが、回帰推定量とDR推定量は不偏を維持しています。アウトカムモデルが正しいため、回帰とDRは正しく機能します
-
両方のモデルが誤指定（右下）: 3つの推定量すべてにバイアスが生じます。DR推定量の「二重のロバスト性」はどちらか一方が正しい場合にのみ発揮され、両方が誤っていれば万能ではありません
この結果は、DR推定量の「二重ロバスト性」を如実に示しています。2つのモデルのうちどちらが正しいかを事前に知ることはできないため、DR推定量は実践上の安全策として極めて有用です。
まとめ
本記事では、平均処置効果（ATE）をはじめとする因果効果の指標の定義と推定方法を体系的に解説しました。
- ATEは集団全体の平均的な因果効果であり、ATT（処置群のATE）とATC（対照群のATE）の加重平均として分解できる
- CATEは共変量に条件付けた処置効果であり、効果の異質性を捉える
- 回帰推定量はアウトカムモデルに基づく推定であり、モデルが正しければ一致推定量になる
- IPW推定量は傾向スコアに基づく重み付けにより推定を行い、傾向スコアモデルが正しければ一致推定量になる
- DR推定量（AIPW推定量）はアウトカムモデルと傾向スコアモデルのどちらか一方が正しければ一致推定量になるという二重ロバスト性を持つ
- モンテカルロシミュレーションにより、DR推定量のロバスト性と効率性を確認した
これらの推定量は、因果推論の実践における基本的なツールキットを構成しています。次の記事では、因果効果を推定するための最も直接的なデザインであるランダム化比較試験（RCT）を詳しく解説します。
次のステップとして、以下の記事も参考にしてください。
- ランダム化比較試験（RCT） — ATEを推定する「ゴールドスタンダード」
- 傾向スコア法 — IPW推定量の詳細とマッチング法
- 因果フォレスト — CATEの非線形な推定手法