未観測交絡への感度分析

🎓 レベル：発展　|　重要度：A（必須）

📎 前提：識別の仮定　|　二重頑健推定AIPW

要点（BLUF）

観察研究の因果主張は「未観測交絡がない（条件付き交換可能性 $Y(a)\perp X\mid C$ ）」という検証できない仮定に立つ。感度分析はこの仮定が崩れたとき、結論がどれだけ揺らぐかを定量化する。
E値（E-value, VanderWeele & Ding 2017） ＝観測された関連を未観測交絡だけで完全に説明し去るために、その交絡が処置とも結果とも持っていなければならない最小の関連強度（リスク比スケール）。リスク比 $\mathrm{RR}\ge 1$ なら $E = \mathrm{RR} + \sqrt{\mathrm{RR}(\mathrm{RR}-1)}$ 。
真の効果が厳密にゼロのデータで見かけの $\mathrm{RR}$ が出る例を作り、そこから E値を計算する。**観測済み共変量の交絡強度と比べて「その強さの未観測交絡はもっともらしいか」**を判断するのが感度分析の使い方だ。

1. 概念：識別は検証できない――だから感度分析

二重頑健推定AIPW までの調整法は、すべて条件付き交換可能性（観測した共変量 $C$ で層別すれば処置がランダムと見なせる）に依存していた。だがこの仮定はデータから検証できない。未観測の交絡 $U$ が処置 $X$ と結果 $Y$ の両方を動かしていれば、どんなに丁寧に $C$ を調整しても推定値は偏る。

flowchart LR
    U["未観測交絡 U"] --> X["処置 X"]
    U --> Y["結果 Y"]
    C["観測共変量 C（調整済）"] --> X
    C --> Y
    X --> Y

感度分析は問いを反転させる。「未観測交絡はあるか/ないか」（答えられない）ではなく、**「観測された関連を未観測交絡だけで消し去るには、その交絡はどれほど強くなければならないか」**を計算する。求めた強さが「ありえないほど強い」なら結論は頑健、「いかにもありそう」なら結論は脆い。これを一目で伝える指標が E値だ。

2. 識別：E値の導出（バイアス係数 → 対称な最悪ケース）

二値の未観測交絡 $U$ を考える。Ding と VanderWeele（2016）は、交絡 $U$ が作れる最大のバイアス係数を次の2つの関連強度で抑えた。

$\mathrm{RR}_{EU}$ ：処置群と対照群のあいだで $U$ がどれだけ偏るか（処置 → $U$ のリスク比）
$\mathrm{RR}_{UD}$ ： $U$ の有無で結果のリスクがどれだけ変わるか（ $U$ → 結果のリスク比）

このとき、観測リスク比 $\mathrm{RR}_{\text{obs}}$ と真のリスク比 $\mathrm{RR}_{\text{true}}$ の比（＝交絡が水増しした倍率）は、高々

B = \frac{\mathrm{RR}_{EU}\,\mathrm{RR}_{UD}}{\mathrm{RR}_{EU} + \mathrm{RR}_{UD} - 1}

までしか大きくならない（これがシャープな上界である点が VanderWeele の貢献）。観測された関連を「実は効果ゼロ（ $\mathrm{RR}_{\text{true}}=1$ ）」へ完全に説明し去るには、 $B \ge \mathrm{RR}_{\text{obs}}$ が必要になる。

ここで「最小の強さ」を求めるため、 $\mathrm{RR}_{EU}$ と $\mathrm{RR}_{UD}$ の大きいほうを最小化する。対称な点 $\mathrm{RR}_{EU}=\mathrm{RR}_{UD}=E$ で最小になり、 $B=\mathrm{RR}_{\text{obs}}$ と置くと

\frac{E^2}{2E-1} = \mathrm{RR}_{\text{obs}} \quad\Longrightarrow\quad E^2 - 2\,\mathrm{RR}_{\text{obs}}\,E + \mathrm{RR}_{\text{obs}} = 0.

正の根を取って

E = \mathrm{RR}_{\text{obs}} + \sqrt{\mathrm{RR}_{\text{obs}}\,(\mathrm{RR}_{\text{obs}}-1)}.

これが E値だ。 $\mathrm{RR}_{\text{obs}}<1$ （予防的な効果）のときは $1/\mathrm{RR}_{\text{obs}}$ に直してから同じ式に入れる（E値は $1$ をはさんで対称）。

❓ 確認：E値が大きいほど、結論は未観測交絡に対して頑健か脆いか？ ―― 頑健。大きい E値は「説明し去るのにそれだけ強い交絡が要る＝そう簡単には覆らない」を意味する。

3. 実証（1）：E値は観測効果とともに増える

まず E値の関数を定義し、いくつかの観測 $\mathrm{RR}$ に対する E値を表にする。

import numpy as np

# === 観測リスク比から E 値（E-value）を計算する ===
def e_value(rr):
    # RR<1 のときは逆数にして RR>=1 に揃える（E値は1未満/超で対称）
    if rr < 1.0:
        rr = 1.0 / rr
    return rr + np.sqrt(rr * (rr - 1.0))

for rr in [1.0, 1.25, 1.5, 2.0, 3.0, 5.0]:
    print(f"観測RR = {rr:5.2f}  ->  E値 = {e_value(rr):.3f}")

出力は次の通り。

観測RR =  1.00  ->  E値 = 1.000
観測RR =  1.25  ->  E値 = 1.809
観測RR =  1.50  ->  E値 = 2.366
観測RR =  2.00  ->  E値 = 3.414
観測RR =  3.00  ->  E値 = 5.449
観測RR =  5.00  ->  E値 = 9.472

出力の意味――観測 $\mathrm{RR}=1$ （関連なし）なら E値も $1$ （交絡ゼロでも説明できる）。観測効果が大きくなるほど E値は急に増える。 $\mathrm{RR}=2$ を消すには両方向に $\mathrm{RR}\approx 3.4$ の交絡が、 $\mathrm{RR}=5$ を消すには $\mathrm{RR}\approx 9.5$ もの交絡が要る。強い関連ほど、それを偶然の交絡で片付けるのは難しいことが数値で見える。

4. 実証（2）：真の効果ゼロのデータで E値を読む

ここが核心だ。処置 $X$ は結果 $Y$ に一切影響しない（真の効果ゼロ）が、未観測交絡 $U$ が両方を動かすデータを作る。素朴な（ $U$ を無視した）リスク比から E値を計算し、実際に仕込んだ交絡の強さと突き合わせる。

import numpy as np

# === 真の処置効果ゼロ。未観測交絡 U だけで見かけのリスク比(RR>1)が出る ===
rng = np.random.default_rng(7)
n = 200_000

U = rng.binomial(1, 0.5, size=n)                 # 未観測交絡（二値）
p_treat = np.where(U == 1, 0.8, 0.2)             # U=1 で処置されやすい
X = rng.binomial(1, p_treat)                     # 処置
p_out = np.where(U == 1, 0.40, 0.10)             # 結果は U に依存（X は不在）
Y = rng.binomial(1, p_out)                       # 真の処置効果は厳密にゼロ

def e_value(rr):
    if rr < 1.0:
        rr = 1.0 / rr
    return rr + np.sqrt(rr * (rr - 1.0))

# --- 素朴な(交絡を無視した)リスク比とその信頼区間 ---
a = Y[X == 1].sum(); n1 = (X == 1).sum()
c = Y[X == 0].sum(); n0 = (X == 0).sum()
p1, p0 = a / n1, c / n0
RR_obs = p1 / p0
se_log = np.sqrt(1/a - 1/n1 + 1/c - 1/n0)        # log RR の標準誤差（Katz）
lo = np.exp(np.log(RR_obs) - 1.96 * se_log)
hi = np.exp(np.log(RR_obs) + 1.96 * se_log)
print("真のRR = 1.000（処置は結果に無関係）")
print(f"観測RR = {RR_obs:.3f}  95%CI = [{lo:.3f}, {hi:.3f}]")
print(f"点推定の E値 = {e_value(RR_obs):.3f}")
print(f"CI下限の E値 = {e_value(lo):.3f}")

# --- 実際に仕込んだ交絡 U の強さ（リスク比スケール）---
RR_UD = Y[U == 1].mean() / Y[U == 0].mean()      # U→結果のリスク比
RR_EU = U[X == 1].mean() / U[X == 0].mean()      # 処置→U のリスク比
print(f"\n仕込んだ交絡の強さ: U→結果 RR = {RR_UD:.2f}, 処置→U RR = {RR_EU:.2f}")
B = (RR_EU * RR_UD) / (RR_EU + RR_UD - 1.0)
print(f"この交絡が作れる最大バイアス係数 B = {B:.3f} >= 観測RR {RR_obs:.3f}")

# --- もし U を観測できれば層別調整 → RR は 1 に戻る ---
num = den = 0.0
for u in [0, 1]:
    m = (U == u); w = m.mean()
    num += w * Y[m & (X == 1)].mean()
    den += w * Y[m & (X == 0)].mean()
print(f"U で層別調整した RR = {num / den:.3f}（真値 1.0 に戻る）")

出力は次の通り。

真のRR = 1.000（処置は結果に無関係）
観測RR = 2.127  95%CI = [2.092, 2.163]
点推定の E値 = 3.675
CI下限の E値 = 3.603

仕込んだ交絡の強さ: U→結果 RR = 3.99, 処置→U RR = 3.98
この交絡が作れる最大バイアス係数 B = 2.280 >= 観測RR 2.127
U で層別調整した RR = 1.003（真値 1.0 に戻る）

出力の意味――いくつもの数字が一本につながる。

真の効果はゼロなのに、 $U$ を無視すると観測 $\mathrm{RR}=2.13$ （信頼区間も $1$ をまたがず「有意」）。素朴な分析なら「処置に効果あり」と誤って結論する。
そこから計算した E値 $=3.68$ 。これは「処置とも結果とも $\mathrm{RR}\approx 3.7$ の関連を持つ未観測交絡があれば、この見かけの効果は説明し去れる」という意味。
実際に仕込んだ交絡 $U$ は 両方向とも $\mathrm{RR}\approx 4.0$ 。E値 $3.68$ を上回っている。だからこの交絡は見かけの効果を作るのに十分で、事実 $B=2.28\ge \mathrm{RR}_{\text{obs}}=2.13$ を満たす。
そして $U$ を観測できて層別調整すれば、RR は $1.00$ に戻る――関連は最初から交絡の産物だった。

E値は「未観測交絡があるか」には答えない。だが「もしあるなら、これだけ強くなければ説明できない」という具体的な基準を与える。あとは分野知識で「 $\mathrm{RR}=3.7$ ほど強い未測定要因はありそうか」を問えばよい。本例では現にそれが存在したので、結論は脆かった。

5. 図：E値曲線

E値が観測効果とともにどう増えるかを描く。

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

# === E値曲線: 観測効果が大きいほど、説明し去るのに強い交絡が要る ===
def e_value(rr):
    return rr + np.sqrt(rr * (rr - 1.0))

rr_grid = np.linspace(1.0, 5.0, 200)
ev = e_value(rr_grid)

fig, ax = plt.subplots(figsize=(7, 4.5))
ax.plot(rr_grid, ev, lw=2.2, color="#1f77b4", label="E値")
ax.plot(rr_grid, rr_grid, ls="--", color="gray", label="参照線 y=x")
for r in [1.5, 2.0, 3.0]:
    ax.plot([r], [e_value(r)], "o", color="crimson")
    ax.annotate(f"RR={r}\nE値={e_value(r):.2f}", (r, e_value(r)),
                textcoords="offset points", xytext=(8, -2), fontsize=9)
ax.set_xlabel("観測リスク比 RR")
ax.set_ylabel("E値（必要な未観測交絡の強さ）")
ax.set_title("観測効果が大きいほど、説明し去るのに強い交絡が要る")
ax.legend()
plt.tight_layout()
plt.show()

E値曲線は常に参照線 $y=x$ より上にある（ $\sqrt{\mathrm{RR}(\mathrm{RR}-1)}>0$ ）。つまり必要な交絡の強さは、観測されたリスク比そのものより必ず大きい。これが「弱い関連は交絡で消えやすく、強い関連は消えにくい」ことの幾何学的な表れだ。

6. 使いどころ

報告のしかた：点推定の E値と、信頼区間の（ $1$ に近い側の）端の E値の両方を添える。本例なら「E値 $3.68$ 、CI下限の E値 $3.60$ 」。CI下限の E値は、標本の偶然も加味した上でどれだけの交絡が要るかを示す。
比較対象：E値は観測済み共変量の交絡強度と比べて読む。測定した交絡がどれも $\mathrm{RR}\approx 2$ 程度なら、E値 $3.7$ の未観測交絡は「測ったどれよりも強い未知の要因」を要求するので、結論は比較的頑健。逆に E値が $1.3$ なら、ありふれた弱い交絡で覆る。
どんな推定値でも使える：オッズ比・ハザード比は近似的にリスク比へ変換してから（VanderWeele は $\mathrm{HR},\mathrm{OR}$ の変換式も与えている）。連続アウトカムは標準化効果量から近似 $\mathrm{RR}$ を作る。

⚠️ よくある誤解・落とし穴

「E値が大きい＝交絡はない」ではない。E値は交絡の有無を判定しない。「もし覆すならこれだけ強い交絡が要る」という条件付きの強さを述べるだけ。最終判断は分野知識による。
E値は未観測”交絡”に対する指標。選択バイアス（衝突点バイアスと選択バイアス）・測定誤差・モデル誤特定には別の感度分析が要る。万能の頑健性スコアではない。
シャープだが最悪ケース。E値は「交絡が最も効率よくバイアスを作る」配置を仮定した境界。現実の交絡はそこまで効率的でないことも多く、E値は保守的（必要強度を低めに見積もる）になりうる。
大きな観測効果に引きずられない。E値が大きいのは結論が頑健だからとは限らず、単に観測効果が大きいだけのこともある。必ず測定済み交絡の強さと相対比較する。