機械学習と情報技術

因果推論の感度分析 — 未観測交絡の影響を定量的に評価する

因果推論

観察研究で「新薬が血圧を平均10mmHg下げる」という推定結果を得たとします。しかし、この推定は「すべての交絡因子が適切に制御されている」という仮定に依存しています。もし、測定していない重要な交絡因子（例えば、遺伝的素因や生活習慣）が存在していたらどうなるでしょうか。推定結果は信頼できるのでしょうか。

因果推論における最大の脅威は未観測の交絡変数（unmeasured confounders）の存在です。傾向スコア法や回帰調整などの手法は、観測された交絡因子を制御することはできますが、観測されていない交絡因子に対しては無力です。そして、未観測交絡の有無はデータから検証することができません — これは原理的な制約です。

感度分析（sensitivity analysis）は、この根本的な不確実性に対処するための方法論です。「未観測交絡がどの程度存在すれば推定結果が覆るか」を定量的に評価し、推定結果のロバスト性（頑健性）を判断するためのツールです。

感度分析は以下の場面で不可欠です。

  • 医学研究: 観察研究で得られた治療効果の信頼性を評価する。規制当局（FDAなど）への申請で求められることがある
  • 疫学: 暴露とアウトカムの関連が交絡によって説明されうるかを評価する
  • 経済学: 政策評価において、操作変数法やDIDの仮定が破れた場合の影響を定量化する
  • 社会科学: 観察研究の結果の信頼性を示すために、ジャーナル論文で報告が推奨されている

本記事の内容

  • 感度分析の基本的なアイデアと必要性
  • ローゼンバウムの感度分析 — $\Gamma$ パラメータと感度パラメータの解釈
  • E-value — 未観測交絡の強さの直感的な尺度
  • オミッテッド変数バイアス（OVB）の定量化
  • 部分 $R^2$ に基づく感度分析（Cinelli and Hazlett, 2020）
  • 感度分析の実践的なガイドライン
  • Pythonによるシミュレーションと可視化

前提知識

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

感度分析の基本アイデア

未観測交絡はなぜ問題か

因果推論入門で学んだように、観察データから因果効果を推定するためには無交絡仮定（unconfoundedness）が必要です。

$$ (Y(0), Y(1)) \perp\!\!\!\perp T \mid \bm{X} $$

この仮定は「観測された共変量 $\bm{X}$ で条件付ければ、処置の割り当てはポテンシャルアウトカムと独立になる」と言っています。しかし、もし $\bm{X}$ に含まれていない交絡変数 $U$ が存在すれば、この仮定は破れ、推定結果にバイアスが生じます。

重要なのは、この仮定はデータから検証不可能であるという点です。どれだけ多くの共変量を観測して制御しても、「まだ制御していない交絡変数が存在するかもしれない」という可能性を排除できません。

感度分析のアプローチ

感度分析は、この検証不可能な仮定に対して次のように対処します。

  1. 未観測交絡の強さをパラメータ化する: 未観測交絡の強さを1つまたは複数のパラメータで表現する
  2. パラメータの値を変化させて推定結果の変化を追跡する: 未観測交絡が弱い場合（パラメータが小さい）から強い場合（パラメータが大きい）まで、推定結果がどう変わるかを系統的に調べる
  3. 臨界値を特定する: 推定結果が覆る（効果がゼロになる、p値が有意水準を超える、信頼区間がゼロを含む）ために必要な未観測交絡の強さを特定する
  4. 臨界値の妥当性を議論する: 臨界値が「現実的に起こりうる」強さなのかどうかを、ドメイン知識や観測された交絡因子の強さと比較して判断する

大雑把に言えば、感度分析は「結論を覆すために必要な未観測交絡のハードル」を示すものです。このハードルが高いほど（強い未観測交絡がないと結論が覆らないほど）、推定結果はロバストだと言えます。

感度分析の重要性は、近年の方法論的進展とともに急速に認知されています。『New England Journal of Medicine』や『Journal of the American Statistical Association』などの一流ジャーナルでは、観察研究の報告において感度分析の結果を含めることが事実上の標準となりつつあります。感度分析がなければ、推定結果がどの程度の未観測交絡に耐えうるかが不明であり、結論の信頼性に関する建設的な議論が困難になるためです。

では、具体的な感度分析の手法を見ていきましょう。最初に紹介するのは、最も古典的かつ広く使われているローゼンバウムの枠組みです。

ローゼンバウムの感度分析

基本的な設定

ポール・ローゼンバウム（Paul Rosenbaum）が開発した感度分析の枠組みは、マッチングに基づく因果推論に対する感度分析として最も広く使われています。

傾向スコアマッチングを行った後の設定を考えます。マッチングにより、観測された共変量 $\bm{X}$ が類似した処置群と対照群のペアが構成されています。無交絡仮定が成り立つなら、各ペア内で処置を受ける確率は等しいはずです。

しかし、未観測の交絡変数 $U$ が存在する場合、$\bm{X}$ が同じでも処置確率が異なります。ローゼンバウムはこの乖離を $\Gamma$ パラメータで定量化します。

$\Gamma$ パラメータの定義

観測された共変量 $\bm{X}$ が同じ2つの個体 $i$ と $j$ について、処置のオッズ比の限界を $\Gamma$ で制約します。

$$ \begin{equation} \frac{1}{\Gamma} \leq \frac{\pi_i / (1 – \pi_i)}{\pi_j / (1 – \pi_j)} \leq \Gamma \end{equation} $$

ここで $\pi_i = P(T_i = 1 \mid \bm{X}_i, U_i)$ は、未観測交絡を含めた場合の個体 $i$ の真の処置確率です。

$\Gamma$ の解釈は直感的です。

  • $\Gamma = 1$: 未観測交絡なし。$\bm{X}$ が同じなら処置確率も同じ
  • $\Gamma = 2$: 未観測交絡により、$\bm{X}$ が同じ個体間で処置オッズが最大2倍異なりうる
  • $\Gamma = 3$: 処置オッズが最大3倍異なりうる
  • $\Gamma = 5$: 処置オッズが最大5倍異なりうる

感度分析の手順

ローゼンバウムの感度分析は以下のステップで行います。

ステップ1: $\Gamma = 1$（未観測交絡なし）で推定結果（p値、信頼区間）を計算する。これがベースラインの結果です。

ステップ2: $\Gamma$ を1から徐々に増加させ、各 $\Gamma$ のもとでp値の上限（worst-case p-value）を計算する。$\Gamma$ が大きくなるほど、未観測交絡の可能性をより多く許容するため、p値の上限は大きくなります。

ステップ3: p値の上限が有意水準 $\alpha$（通常0.05）を超える最小の $\Gamma$ を見つける。これを感度パラメータの臨界値 $\Gamma^*$ と呼びます。

$$ \Gamma^* = \min\{\Gamma : p_{\text{upper}}(\Gamma) > \alpha\} $$

ステップ4: $\Gamma^*$ の値を、観測された交絡因子の強さや先行研究の知見と比較して、推定結果のロバスト性を判断する。

臨界値の解釈

$\Gamma^*$ が大きいほど、推定結果はロバストです。

  • $\Gamma^* < 1.5$: 推定結果は未観測交絡に対して脆弱。わずかな未観測交絡で結果が覆る
  • $1.5 \leq \Gamma^* < 2.5$: 中程度のロバスト性。ある程度の未観測交絡があっても結果は維持される
  • $\Gamma^* \geq 2.5$: 推定結果はかなりロバスト。かなり強い未観測交絡がないと結果は覆らない

ただし、これらの目安は絶対的なものではありません。重要なのは、$\Gamma^*$ の値が「現実的に起こりうる未観測交絡の強さ」と比較してどうかを議論することです。例えば、観測された最も強い交絡因子のオッズ比が1.5程度であるならば、$\Gamma^* = 3$ はかなりロバストと判断できます。

ローゼンバウムの感度分析の数学的背景

ローゼンバウムの感度分析がp値の上限を計算する仕組みをもう少し詳しく見てみましょう。

マッチされたペア $i$ について、処置群の個体が処置を受ける確率を $\pi_i$、対照群の個体が処置を受ける確率を $\pi_{i’}$ とします。$\Gamma$ パラメータのもとで、処置のオッズ比は次の範囲に制約されます。

$$ \frac{1}{\Gamma} \leq \frac{\pi_i(1 – \pi_{i’})}{\pi_{i’}(1 – \pi_i)} \leq \Gamma $$

帰無仮説（処置効果ゼロ）のもとで、各マッチドペアにおける検定統計量の分布は $\pi_i / (\pi_i + \pi_{i’})$ に依存します。$\Gamma = 1$ のとき、この確率は0.5であり、標準的なマッチドペア検定（符号検定やウィルコクソンの符号順位検定）が得られます。

$\Gamma > 1$ のとき、各ペア内の処置確率は等しくなくなり、検定統計量の分布が変化します。p値の上限は、帰無仮説と $\Gamma$ の制約のもとで最悪のケース（p値が最大になるケース）を計算することで得られます。この計算は組合せ最適化問題であり、効率的なアルゴリズムが開発されています。

ローゼンバウムの感度分析はマッチングデータに特化した方法ですが、より一般的な設定で使える指標として、次に紹介するE-valueがあります。

E-value — 直感的な感度分析指標

E-valueの定義

E-value（VanderWeele and Ding, 2017）は、推定された因果効果をゼロに減少させるために必要な未観測交絡の最小の強さを表す指標です。ローゼンバウムの $\Gamma$ パラメータと比べて、より一般的な設定で適用でき、解釈も直感的です。

リスク比（Risk Ratio, RR）で表された推定効果に対して、E-valueは次のように定義されます。

$$ \begin{equation} \text{E-value} = \text{RR}_{\text{obs}} + \sqrt{\text{RR}_{\text{obs}} \times (\text{RR}_{\text{obs}} – 1)} \end{equation} $$

ここで $\text{RR}_{\text{obs}}$ は観測された（調整後の）リスク比です。$\text{RR}_{\text{obs}} < 1$ の場合は $1/\text{RR}_{\text{obs}}$ に置き換えてから計算します。

E-valueの解釈

E-valueの解釈は明快です。

「E-value = $k$」は、未観測の交絡変数 $U$ が、処置とアウトカムのそれぞれに対してリスク比で $k$ 以上の関連を持たなければ、観測された効果を完全に説明できないことを意味します。

例えば、E-value = 3であれば、「未観測交絡が処置のリスクを3倍以上にし、かつアウトカムのリスクも3倍以上にしなければ、観測された効果は交絡では説明できない」ということです。

逆に、E-value = 1.2であれば、処置とアウトカムにそれぞれリスク比1.2程度の関連を持つ未観測交絡が存在するだけで効果が説明されてしまうため、推定結果は脆弱です。

信頼区間のE-value

推定効果の点推定だけでなく、信頼区間の端点に対してもE-valueを計算することが推奨されます。

例えば、リスク比の推定値が2.5で、95%信頼区間の下限が1.8であれば:

  • 点推定のE-value: $2.5 + \sqrt{2.5 \times 1.5} \approx 4.44$
  • 信頼区間下限のE-value: $1.8 + \sqrt{1.8 \times 0.8} \approx 3.00$

信頼区間下限のE-valueは「推定結果を統計的に非有意にするために必要な未観測交絡の強さ」を表しており、点推定のE-valueよりも保守的な指標です。

連続アウトカムの場合

連続アウトカムで効果量がコーエンの $d$ として報告される場合、近似的なリスク比への変換が必要です。

$$ \begin{equation} \text{RR}_{\text{approx}} = \exp(0.91 \times d) \end{equation} $$

この変換は、正規分布の仮定のもとでの近似です。変換後のRRをE-valueの公式に代入します。

他にも、オッズ比、ハザード比、平均差などの効果量に対するE-valueの計算法が提案されています。

E-valueの実践的な活用

E-valueの実践的な強みは、異なる研究間でロバスト性を比較できることです。例えば、ある研究でE-value = 5.2、別の研究でE-value = 1.8であれば、前者の方が未観測交絡に対してはるかにロバストであると判断できます。

E-valueのもう一つの利点は、ドメイン知識との橋渡しが容易なことです。研究者は「この分野で知られている最も強い交絡因子のリスク比はどの程度か」を考え、それとE-valueを比較することで、推定結果のロバスト性を実質的に評価できます。

例えば、喫煙と肺がんの関係に関する研究で、観測された最も強い交絡因子（例えばアスベスト暴露）のリスク比が約2であったとします。このとき、E-value = 4であれば「アスベスト暴露の2倍の強さの未観測交絡がないと結果は覆らない」と解釈でき、推定結果は比較的ロバストと判断できます。

ここまで紹介した2つの手法は、いずれも「未観測交絡の強さ」を1つのパラメータで要約するアプローチでした。次に、より精密なアプローチとして、オミッテッド変数バイアスの定量化を紹介します。

オミッテッド変数バイアス（OVB）の定量化

OVBの数学的構造

回帰分析において、重要な変数 $U$ を省略した場合に生じるバイアスをオミッテッド変数バイアス（Omitted Variable Bias, OVB）と呼びます。

真のモデルが

$$ Y = \beta_0 + \beta_1 T + \beta_2 U + \epsilon $$

であるのに、$U$ を省略して

$$ Y = \tilde{\beta}_0 + \tilde{\beta}_1 T + \tilde{\epsilon} $$

を推定した場合、OLSの性質から $\tilde{\beta}_1$ のバイアスは次のように表されます。

$$ \begin{equation} \text{Bias} = E[\tilde{\beta}_1] – \beta_1 = \beta_2 \cdot \delta_{TU} \end{equation} $$

ここで $\delta_{TU}$ は $T$ を $U$ に回帰した際の係数、すなわち $T$ と $U$ の（他の共変量で条件付けた）関連の強さです。

OVBの構造は直感的に理解できます。バイアスは2つの要素の積です。

  1. $\beta_2$: 省略された変数 $U$ がアウトカム $Y$ に与える影響の強さ
  2. $\delta_{TU}$: 省略された変数 $U$ と処置 $T$ の関連の強さ

両方が大きい場合にのみ、大きなバイアスが生じます。$U$ が $Y$ に強く影響しても $T$ と無関係なら（$\delta_{TU} = 0$）、バイアスはゼロです。

部分 $R^2$ に基づく感度分析（Cinelli and Hazlett, 2020）

Cinelli and Hazlett（2020）は、OVBの大きさを部分 $R^2$（partial $R^2$）を用いて表現する感度分析のフレームワークを提案しました。このアプローチは、未観測交絡の「強さ」を、観測された交絡因子の強さとの比較で評価できるという利点があります。

未観測交絡変数 $U$ の影響は、以下の2つの部分 $R^2$ で特徴づけられます。

$$ R^2_{Y \sim U | T, \bm{X}}: \quad U \text{がアウトカム}Y\text{の分散をどれだけ説明するか（}T\text{と}\bm{X}\text{で条件付け）} $$ $$ R^2_{T \sim U | \bm{X}}: \quad U \text{が処置}T\text{の分散をどれだけ説明するか（}\bm{X}\text{で条件付け）} $$

処置効果の推定値が $\hat{\beta}_1$ であるとき、未観測交絡 $U$ によるバイアスの大きさは、これら2つの部分 $R^2$ の関数として上界が求まります。

感度コンターの可視化

部分 $R^2$ に基づく感度分析の強みは、2つの部分 $R^2$ の空間上にコンタープロット（等値線図）を描けることです。横軸に $R^2_{T \sim U | \bm{X}}$（$U$ と $T$ の関連の強さ）、縦軸に $R^2_{Y \sim U | T, \bm{X}}$（$U$ と $Y$ の関連の強さ）を取り、推定された処置効果がゼロになる曲線（バイアス = 推定値となる曲線）を描きます。

この曲線の外側（右上側）は推定結果が覆る領域、内側（左下側）は推定結果が維持される領域を表します。さらに、観測された各共変量の部分 $R^2$ をこのプロット上にプロットすることで、「未観測交絡が観測された交絡因子と比較してどの程度強ければ結果が覆るか」を視覚的に評価できます。

例えば、最も強い観測された交絡因子の部分 $R^2$ が $(0.05, 0.08)$ であり、推定結果が覆る曲線が $(0.15, 0.20)$ の領域を通っているなら、「最も強い観測された交絡因子の約3倍の強さを持つ未観測交絡がないと結果は覆らない」と結論できます。

バイアスの公式の導出

Cinelli and Hazlett（2020）のフレームワークにおけるバイアスの式をより具体的に導出しましょう。回帰モデル $Y = \beta_1 T + \bm{\gamma}^\top \bm{X} + \delta U + \epsilon$ で、$U$ が省略されている場合、$\beta_1$ の推定量のバイアスは以下のように表されます。

$$ |\text{Bias}| \leq \hat{\text{SE}} \cdot \sqrt{\frac{R^2_{Y \sim U | T, \bm{X}} \cdot R^2_{T \sim U | \bm{X}}}{1 – R^2_{T \sim U | \bm{X}}}} \cdot \frac{1}{\hat{\text{SE}}} $$

ここで $\hat{\text{SE}}$ は処置係数の推定標準誤差です。この式は、バイアスが2つの部分 $R^2$ の幾何平均に比例することを示しています。一方の部分 $R^2$ がゼロであれば（未観測交絡が処置かアウトカムのどちらか一方に全く影響しない場合）、バイアスはゼロになります。

この公式の重要な帰結は、「処置の説明力が高い未観測交絡」と「アウトカムの説明力が高い未観測交絡」の両方が同時に存在する場合にのみ大きなバイアスが生じるということです。どちらか一方だけでは不十分であり、この点は感度分析の解釈において重要です。

ここまでの理論を、Pythonシミュレーションで実装してみましょう。

Pythonシミュレーション

シミュレーション1: E-valueの計算と解釈

import numpy as np
import matplotlib.pyplot as plt

# --- E-valueの計算関数 ---
def e_value(RR):
    """リスク比からE-valueを計算する"""
    if RR < 1:
        RR = 1 / RR
    return RR + np.sqrt(RR * (RR - 1))

def e_value_from_d(d):
    """コーエンのdからE-valueを計算（近似変換）"""
    RR = np.exp(0.91 * abs(d))
    return e_value(RR)

# E-valueの計算例
risk_ratios = [1.2, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0]
print("=== リスク比別 E-value ===")
print(f"{'RR':>6}  {'E-value':>10}")
print("-" * 20)
for rr in risk_ratios:
    ev = e_value(rr)
    print(f"{rr:>6.1f}  {ev:>10.3f}")

print("\n=== コーエンのd別 E-value ===")
print(f"{'d':>6}  {'RR_approx':>12}  {'E-value':>10}")
print("-" * 34)
for d in [0.2, 0.5, 0.8, 1.0, 1.5]:
    rr = np.exp(0.91 * d)
    ev = e_value_from_d(d)
    print(f"{d:>6.1f}  {rr:>12.3f}  {ev:>10.3f}")

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

# (a) E-value vs RR
ax = axes[0]
rr_range = np.linspace(1.01, 6, 200)
ev_range = [e_value(rr) for rr in rr_range]
ax.plot(rr_range, ev_range, 'b-', linewidth=2)
ax.fill_between(rr_range, rr_range, ev_range, alpha=0.15, color='blue',
                label='Additional confounding needed')
ax.plot(rr_range, rr_range, 'gray', linestyle='--', linewidth=1,
        label='E-value = RR (identity)')
ax.set_xlabel('Observed Risk Ratio (RR)', fontsize=12)
ax.set_ylabel('E-value', fontsize=12)
ax.set_title('E-value as a Function of Risk Ratio', fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(1, 6)
ax.set_ylim(1, 12)

# (b) E-value vs Cohen's d
ax = axes[1]
d_range = np.linspace(0.1, 2.0, 200)
ev_d_range = [e_value_from_d(d) for d in d_range]
ax.plot(d_range, ev_d_range, 'r-', linewidth=2)

# 参考線
for d_ref, label, color in [(0.2, 'Small', 'green'), (0.5, 'Medium', 'orange'),
                              (0.8, 'Large', 'red')]:
    ev_ref = e_value_from_d(d_ref)
    ax.plot(d_ref, ev_ref, 'o', color=color, markersize=8)
    ax.annotate(f'{label}\n(d={d_ref}, E={ev_ref:.2f})',
                xy=(d_ref, ev_ref), xytext=(d_ref + 0.15, ev_ref + 0.5),
                fontsize=9, arrowprops=dict(arrowstyle='->', color=color))

ax.set_xlabel("Cohen's d", fontsize=12)
ax.set_ylabel('E-value', fontsize=12)
ax.set_title("E-value as a Function of Cohen's d", fontsize=12)
ax.grid(True, alpha=0.3)

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

左のグラフはリスク比とE-valueの関係を示しています。E-valueは常にRRより大きく（灰色の破線が恒等線）、RRが大きいほどE-valueも急速に増加します。青い塗りつぶし領域は「RRを超えてさらに必要な交絡の強さ」を表しています。RR=2.0の場合、E-value=3.41となり、未観測交絡は処置とアウトカムのそれぞれに対してRR=3.41以上の関連を持つ必要があります。

右のグラフはコーエンのdとE-valueの関係を示しています。小さな効果量（d=0.2）ではE-value=1.47と低く、中程度以下の未観測交絡で結果が覆ります。大きな効果量（d=0.8）ではE-value=3.85と高く、かなり強い未観測交絡が必要です。

シミュレーション2: ローゼンバウムの感度分析の実装

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

np.random.seed(42)
n = 500

# --- データ生成 ---
X = np.random.normal(0, 1, n)  # 観測された共変量
U = np.random.normal(0, 1, n)  # 未観測の交絡変数

# 処置確率（XとUに依存）
gamma_U = 0.8  # 未観測交絡の強さ
logit_e = 0.3 + 0.5 * X + gamma_U * U
prob_T = 1 / (1 + np.exp(-logit_e))
T = np.random.binomial(1, prob_T, n)

# アウトカム（真の因果効果 = 2.0）
true_effect = 2.0
delta_U = 1.0  # U -> Y
Y = 5 + true_effect * T + 1.0 * X + delta_U * U + np.random.normal(0, 1, n)

# --- 傾向スコアマッチング（簡易版）---
from sklearn.linear_model import LogisticRegression
ps_model = LogisticRegression().fit(X.reshape(-1, 1), T)
ps = ps_model.predict_proba(X.reshape(-1, 1))[:, 1]

# ナイーブなATT推定（Xのみで制御）
from sklearn.linear_model import LinearRegression
reg = LinearRegression().fit(np.column_stack([T, X]), Y)
naive_effect = reg.coef_[0]

# 正しい推定（Uも制御）
reg_correct = LinearRegression().fit(np.column_stack([T, X, U]), Y)
correct_effect = reg_correct.coef_[0]

print(f"真の効果: {true_effect:.4f}")
print(f"ナイーブ推定（Xのみ）: {naive_effect:.4f}")
print(f"正しい推定（X+U）: {correct_effect:.4f}")
print(f"バイアス: {naive_effect - true_effect:.4f}")

# --- ローゼンバウム型の感度分析（シミュレーション） ---
gammas = np.linspace(1, 5, 50)
p_upper_values = []
effect_ranges = []

# 各Gammaでのworst-case推定
n_sims = 500
for Gamma in gammas:
    effects_worst = []
    for _ in range(n_sims):
        # Gammaに対応する未観測交絡を生成
        U_sim = np.random.normal(0, 1, n)
        logit_sim = 0.3 + 0.5 * X + np.log(Gamma) * U_sim
        prob_sim = 1 / (1 + np.exp(-logit_sim))
        T_sim = np.random.binomial(1, prob_sim, n)

        Y_sim = 5 + true_effect * T_sim + 1.0 * X + np.random.normal(0, 1, n)
        reg_sim = LinearRegression().fit(np.column_stack([T_sim, X]), Y_sim)
        effects_worst.append(reg_sim.coef_[0])

    effects_worst = np.array(effects_worst)
    # p値の上限（推定値がゼロを含む確率）
    z = np.mean(effects_worst) / np.std(effects_worst)
    p_upper = 2 * (1 - stats.norm.cdf(abs(z)))
    p_upper_values.append(p_upper)
    effect_ranges.append((np.percentile(effects_worst, 5),
                          np.percentile(effects_worst, 95)))

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

# (a) 推定効果の範囲 vs Gamma
ax = axes[0]
lower = [r[0] for r in effect_ranges]
upper = [r[1] for r in effect_ranges]
ax.fill_between(gammas, lower, upper, alpha=0.3, color='steelblue',
                label='90% range of estimates')
ax.plot(gammas, [(l+u)/2 for l, u in effect_ranges], 'b-', linewidth=2,
        label='Mean estimate')
ax.axhline(true_effect, color='green', linestyle='--', linewidth=1.5,
           label=f'True effect = {true_effect}')
ax.axhline(0, color='red', linestyle=':', linewidth=1.5, label='Zero effect')
ax.set_xlabel('$\\Gamma$ (sensitivity parameter)', fontsize=12)
ax.set_ylabel('Estimated treatment effect', fontsize=12)
ax.set_title('Sensitivity of Treatment Effect to $\\Gamma$', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# (b) p値 vs Gamma
ax = axes[1]
ax.plot(gammas, p_upper_values, 'b-', linewidth=2)
ax.axhline(0.05, color='red', linestyle='--', linewidth=1.5, label='$\\alpha = 0.05$')
ax.set_xlabel('$\\Gamma$ (sensitivity parameter)', fontsize=12)
ax.set_ylabel('p-value (upper bound)', fontsize=12)
ax.set_title('Upper Bound on p-value', fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, max(0.2, max(p_upper_values) + 0.02))

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

左のグラフは、感度パラメータ $\Gamma$ に対する推定効果の範囲を示しています。$\Gamma = 1$（未観測交絡なし）では推定値が真の効果付近にありますが、$\Gamma$ の増加に伴って推定値の範囲が広がります。ゼロ効果のラインとの交差点が臨界値です。

右のグラフは、$\Gamma$ に対するp値の上限を示しています。$\Gamma$ の増加とともにp値が上昇し、$\alpha = 0.05$ のライン（赤の破線）を超える点が臨界値 $\Gamma^*$ に対応します。

シミュレーション3: OVBとコンターの可視化

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

np.random.seed(42)
n = 2000

# 共変量
X1 = np.random.normal(0, 1, n)
X2 = np.random.normal(0, 1, n)
X3 = np.random.normal(0, 1, n)

# 処置（X1, X2, X3に依存）
T_logit = 0.3 * X1 + 0.5 * X2 + 0.2 * X3
T = (np.random.uniform(0, 1, n) < 1 / (1 + np.exp(-T_logit))).astype(float)

# アウトカム（真の効果 = 1.5）
true_beta = 1.5
Y = true_beta * T + 0.8 * X1 + 1.2 * X2 + 0.5 * X3 + np.random.normal(0, 1, n)

# --- 各共変量の部分R^2を計算 ---
def partial_r2(Y, X_full, X_reduced):
    """部分R^2を計算"""
    reg_full = LinearRegression().fit(X_full, Y)
    reg_reduced = LinearRegression().fit(X_reduced, Y)
    SS_full = np.sum((Y - reg_full.predict(X_full))**2)
    SS_reduced = np.sum((Y - reg_reduced.predict(X_reduced))**2)
    return 1 - SS_full / SS_reduced

# X1, X2, X3それぞれの部分R^2
X_all = np.column_stack([T, X1, X2, X3])

# T~Xi|残りのX（処置との関連）
pr2_T_X1 = partial_r2(T, np.column_stack([X1, X2, X3]),
                       np.column_stack([X2, X3]))
pr2_T_X2 = partial_r2(T, np.column_stack([X1, X2, X3]),
                       np.column_stack([X1, X3]))
pr2_T_X3 = partial_r2(T, np.column_stack([X1, X2, X3]),
                       np.column_stack([X1, X2]))

# Y~Xi|T,残りのX（アウトカムとの関連）
pr2_Y_X1 = partial_r2(Y, X_all, np.column_stack([T, X2, X3]))
pr2_Y_X2 = partial_r2(Y, X_all, np.column_stack([T, X1, X3]))
pr2_Y_X3 = partial_r2(Y, X_all, np.column_stack([T, X1, X2]))

print("=== 各共変量の部分R^2 ===")
print(f"{'変数':>4}  {'R^2(T~Xi|X-i)':>14}  {'R^2(Y~Xi|T,X-i)':>17}")
print("-" * 40)
for name, r2t, r2y in [('X1', pr2_T_X1, pr2_Y_X1),
                         ('X2', pr2_T_X2, pr2_Y_X2),
                         ('X3', pr2_T_X3, pr2_Y_X3)]:
    print(f"{name:>4}  {r2t:>14.4f}  {r2y:>17.4f}")

# --- コンタープロット ---
fig, ax = plt.subplots(figsize=(8, 7))

# 推定結果
reg_full = LinearRegression().fit(X_all, Y)
beta_hat = reg_full.coef_[0]
residuals = Y - reg_full.predict(X_all)
se = np.sqrt(np.sum(residuals**2) / (n - 4) * np.linalg.inv(X_all.T @ X_all)[0, 0])

# コンター: 推定効果がゼロになる（近似的な）曲線
r2_t_range = np.linspace(0.001, 0.3, 200)
r2_y_range = np.linspace(0.001, 0.5, 200)
R2T, R2Y = np.meshgrid(r2_t_range, r2_y_range)

# バイアスの近似: bias ≈ beta_hat * sqrt(R2T * R2Y / (1-R2T))
# 推定値がゼロになる条件: bias = beta_hat
bias_surface = beta_hat * np.sqrt(R2T * R2Y / (1 - R2T))
zero_contour = beta_hat  # この値のコンター

cs = ax.contour(R2T, R2Y, bias_surface,
                levels=[beta_hat * 0.25, beta_hat * 0.5, beta_hat * 0.75, beta_hat],
                colors=['lightblue', 'blue', 'darkblue', 'red'],
                linewidths=[1, 1, 1.5, 2.5])
ax.clabel(cs, fmt={beta_hat * 0.25: f'Bias={beta_hat*0.25:.2f}',
                    beta_hat * 0.5: f'Bias={beta_hat*0.5:.2f}',
                    beta_hat * 0.75: f'Bias={beta_hat*0.75:.2f}',
                    beta_hat: f'Effect=0'},
          fontsize=9)

# 観測された共変量をプロット
for name, r2t, r2y, color in [('X1', pr2_T_X1, pr2_Y_X1, 'green'),
                                ('X2', pr2_T_X2, pr2_Y_X2, 'orange'),
                                ('X3', pr2_T_X3, pr2_Y_X3, 'purple')]:
    ax.plot(r2t, r2y, 'o', color=color, markersize=10, label=name)

ax.set_xlabel('Partial $R^2$ of $U$ with Treatment ($R^2_{T \\sim U | X}$)', fontsize=12)
ax.set_ylabel('Partial $R^2$ of $U$ with Outcome ($R^2_{Y \\sim U | T, X}$)', fontsize=12)
ax.set_title('Sensitivity Contour Plot', fontsize=13)
ax.legend(fontsize=10, loc='upper left')
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 0.3)
ax.set_ylim(0, 0.5)

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

コンタープロットから以下のことが読み取れます。

  1. 赤い実線（Effect=0のコンター）が「推定結果が覆る境界」を示しています。未観測交絡の部分 $R^2$ がこの曲線の右上にある場合、推定された効果はゼロまたは負になります。
  2. 観測された共変量（X1, X2, X3の点）はコンターの左下に位置している: これは、観測された共変量と同程度の強さの未観測交絡では推定結果が覆らないことを示唆しています。
  3. 推定結果を覆すには、最も強い観測された交絡因子（X2）よりもかなり大きな部分 $R^2$ を持つ未観測交絡が必要: これは推定結果が比較的ロバストであることを示唆しています。

このようなベンチマーキング（観測された交絡因子との比較）は、感度分析の結果を解釈する上で極めて有用です。

感度分析の実践ガイドライン

報告すべき事項

観察研究の報告において、感度分析の結果を適切に報告するために、以下の事項を含めることが推奨されます。

  1. 使用した手法: どの感度分析手法を使ったか（ローゼンバウム、E-value、OVBなど）
  2. 臨界値: 推定結果が覆るために必要な未観測交絡の強さ（$\Gamma^*$、E-value等）
  3. ベンチマーキング: 臨界値を観測された交絡因子の強さと比較した議論
  4. 信頼区間の感度分析: 点推定だけでなく信頼区間に対する感度分析の結果

複数の手法の併用

異なる感度分析手法は異なる仮定に基づいており、異なる側面からロバスト性を評価します。可能であれば、複数の手法を併用することで、より包括的なロバスト性の評価が可能になります。

感度分析の限界

感度分析自体にも限界があることを認識しておく必要があります。

  • 未観測交絡の「形」を仮定している: 多くの手法は、未観測交絡が1つの変数であることを暗黙に仮定している。複数の弱い未観測交絡の合成効果は想定されていない場合がある
  • パラメトリックな仮定への依存: 一部の手法は線形関係や特定の分布を仮定している
  • ロバスト性の判断は主観的: 臨界値が「十分に大きい」かどうかの判断は、ドメイン知識に依存する

因果推定手法ごとの感度分析

ここまで主に回帰調整やマッチングに基づく因果推定の感度分析を議論してきましたが、他の因果推定手法にも固有の感度分析が存在します。

操作変数法の感度分析: 操作変数法の排除制約が「近似的に」しか成り立たない場合、IV推定量のバイアスの大きさを評価するための感度分析が開発されています。Conley, Hansen, and Rossi（2012）は、排除制約からの乖離をパラメータ化し（$Y = \beta T + \gamma Z + \epsilon$ で $\gamma \neq 0$）、$\gamma$ の値を変化させて推定結果の変化を追跡する方法を提案しました。

差の差法（DID）の感度分析: 差の差法の平行トレンド仮定が破れた場合の影響を評価する感度分析として、Rambachan and Roth（2023）は、トレンドの差が最大でどの程度であるかを制約するアプローチを提案しました。処置前のトレンドの差を基準にして、処置後のトレンドの差が同程度以内であれば推定結果がどう変わるかを体系的に評価できます。

回帰不連続デザイン（RDD）の感度分析: RDDでは、帯域幅の選択が推定結果に影響するため、帯域幅を変化させたときの推定結果の安定性を確認することが最も基本的な感度分析です。さらに、ランニング変数の操作（manipulation）に対するマッキャリー密度検定や、共変量のバランス検定も広義の感度分析として機能します。

各手法に適した感度分析を報告することで、推定結果の信頼性について多角的な評価が可能になります。

まとめ

本記事では、因果推論における感度分析の理論と手法を解説しました。

  • 感度分析の必要性: 無交絡仮定はデータから検証できない。未観測交絡がどの程度存在すれば推定結果が覆るかを定量的に評価する必要がある
  • ローゼンバウムの感度分析: $\Gamma$ パラメータで処置オッズの偏りを制約し、$\Gamma$ を増加させてp値の上限を追跡する。臨界値 $\Gamma^*$ が推定結果のロバスト性の指標
  • E-value: 推定効果をゼロにするために必要な未観測交絡の最小の強さを、リスク比で表現する直感的な指標
  • オミッテッド変数バイアス: バイアスは省略変数の「$Y$ への影響 $\times$ $T$ との関連」で決まる。部分 $R^2$ を用いたコンタープロットで視覚的に評価可能
  • ベンチマーキング: 臨界値を観測された交絡因子の強さと比較することで、ロバスト性の実質的な評価が可能

感度分析は、因果推定の結果に「どの程度の信頼を置けるか」を定量的に示すツールです。「未観測交絡がないことを証明する」ことはできませんが、「推定結果を覆すためにはどの程度の未観測交絡が必要か」を示すことで、推定結果の信頼性について建設的な議論を可能にします。

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

補足: 感度分析の歴史的背景

感度分析のアイデアは古くからあります。Cornfield et al.（1959）は、喫煙と肺がんの関係に関する議論の中で、「観測された喫煙と肺がんの関連を完全に説明するためには、未観測交絡が喫煙と肺がんのそれぞれに対して少なくとも9倍のリスク比を持つ必要がある」と論じました。この議論は現代のE-valueの先駆けであり、感度分析の基本的な考え方を含んでいます。

その後、Rosenbaum and Rubin（1983）が傾向スコア法の枠組みで感度分析を形式化し、Rosenbaum（2002）がマッチドペアに対する体系的な感度分析を発展させました。VanderWeele and Ding（2017）のE-valueは、これらの先行研究を統合し、簡便で広く適用可能な指標として提案されたものです。

感度分析の発展は、因果推論が「仮定に依存する」という本質的な限界を認めつつ、その限界を定量的に評価する方法論の成熟を反映しています。

関連記事

前の記事

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

次の記事

選択バイアスとヘックマン補正 — サンプル選択問題の数理と解法

関連タグ: E-value ローゼンバウム限界 ロバスト性 因果推論 感度分析 未観測交絡