Hernan and Robins によって書かれた因果推論の本（本記事では "what if book" と呼びます。著者のサイトにPDFが無料で公開されています）を読んでいく中で、IPW (Inverse Probability Weighting; 逆確率重みづけ) 周りの細かいところを整理するのが難しいなと感じたので、自習のついでに記事としてまとめようと思います。現時点では十分な数の関連論文を読めているわけではないので、調査に抜け漏れがある可能性についてはご容赦いただけると幸いです。

アイキャッチとして使った紀伊国屋のリンクも貼っておきます:

目次はこちらです

記事のゴール
新規性
前提知識
Notation
IPW による Potential outcome の推定の基本: Horvitz-Thompson estimator
Hajek (modified Horvitz-Thompson) estimator
- Hajek estimator とは
- Hajek estimatorの性質
  - 漸近分散以外の部分の議論の整理 (bias of Hajek estimator)
  - 漸近分散の解析
終わりに

記事のゴール

本記事は、IPW に関する連載のスタートの記事です。

本記事では、 single time point exposure の条件下での平均因果効果の推定にあたって、 IPW に関連する以下の4つの主張のうち、1と2について説明を行うことをゴールにします。今後の連載の中で、3と4の説明、ならびに人工データ実験による検証を行う予定です。

IPW によって処置 $a' \in \mathcal{A}$ の Potential outcome $\mathbb{E}[Y^{(a')}$ ] を推定するときの基本は、 $\mathbb{E} [ \frac{I(A, a')Y}{f_{A \mid L}(a' \mid L)}$ ] を estimand とした Horvitz-Thompson estimator の推定である
Hajek (modified Horvitz-Thompson) estimator は不偏推定量ではないが、いくつかの仮定*1のもとで Horvitz-Thompson estimator よりも漸近分散が小さい
Marginal Structural Model のパラメータ推定を行う際に、逆確率で重み付けた最小二乗法を用いることが多いが、これは Hajek estimator の推定に相当する
Stabilized weights は Hajek estimator を推定するときに有用かもしれないが、他の推定量を推定するときには注意が必要 (特に、Doubly robust では適切な使用方法がわからない)

要するに、私が以前 Twitterで投稿していた内容を整理します。Notation はすぐ後で与えます。推定手法については適宜説明します。

<今の理解>
1. single time pointのIPWの基本はE[I(A=a)Y/f(A|L)]をestimandとしたHorvitz-Thompson estimator (HTe) による推定
2. modified HTeの方がHTeより安定
3. IP weighted LS = modified HTeの推定
4. stabilized weightはmodified HTeを推定するときに利用可 (→Doubly robustでは不可)
— fullflu (@fullflu) April 29, 2020

新規性

新規性が気になる方のために、既存の資料と比較しておきます。

以下の点で新規性があると考えています:

因果推論についての有名な資料間での、推定量についての記述の相違を整理した
「Hajek estimator は Potential outcome の期待値についての不偏推定量である」と主張する what if book や星野先生の著書の証明に不備がある可能性を指摘した
Horvitz-Thompson estimator の漸近分散を解析し、 Hajek estimator の方が漸近分散が小さいことを保証できる仮定（Potential outcome の符号一致）を調べた

前提知識

本記事の理解のために最低限必要な前提知識は2つあると考えます。

1つ目は、Potential outcome を用いた因果推論の枠組みを勉強した経験があることです。因果推論そのものの説明は本記事では軽く導入する程度です。

2つ目は、学部教養レベルの統計学を理解していることです。例えば、期待値や不偏性といった概念の説明は本記事では行わないです。

本記事の最後に登場する漸近分散の解析においてのみ、M推定量やモーメント法についての知識を仮定しますが、興味のない方はスルーしていただいても概要の理解には差し支えないかと思います。

Notation

必要なNotationを与えます。

$Y \in \mathbb{R}$ : Outcomeの確率変数です。実数の連続値をとるとします。
$A \in \mathcal{A}$ : 処置の確率変数です。連続か離散か二値かは都度説明します。
$L \in \mathbb{R}^d$ : 共変量の確率変数です。 $d \in \{1, 2, \cdots \}$ 次元のベクトルとします。
$V \subset L$ : Effect modifierの確率変数です。簡単のため、共変量のsubsetである場合を考えます。
$Y^{(a)}$ : 処置 $a \in \mathcal{A}$ についての Potential outcome を意味する確率変数です*2。
$I(\cdot, \cdot)$ : 指示関数です。一つ目の引数と二つ目の引数が等しいときに1を、そうでないときに0を返します。
$f_{A \mid L}(a \mid L)$ : 処置が離散のとき、共変量の値が $L$ であったときの処置 $a$ の確率分布です。処理が連続のとき、 $L$ のもとでの処置 $a$ の確率密度関数です。
$f_A(a)$ : $f_{A \mid L}(a \mid L)$ から共変量についての条件を除いたものです。
$f_{Y \mid A, L}(a, l)$ : 共変量Lと処置Aを引数として、Outcomeを出力する関数です。
single time point exposure: 1時点の介入のみ考慮する実験設定のことです。
$\pi$ : 処置が二値のときの処置確率を簡易表記するためのものです。

その他細かいルールとして、実際のデータを扱う際には、サンプルサイズを $n$ とし、確率変数に対して添え字 $i \in [n$ ] をつけることによって実現値を表現します。例えば、 $i$ 行目のレコードの処置の実現値は $A_i$ として表現します。

また、関数 $f$ の推定結果を $\hat{f}$ で表現します。

IPW による Potential outcome の推定の基本: Horvitz-Thompson estimator

ここでは、冒頭に述べた主張の一つ目の話、Horvitz-Thompson estimator について説明します。

必要な仮定

本題に入る前に、因果効果を適切に推定するために必要な仮定について説明します。

ここからの話の中では、以下の6つの条件が成立することを仮定します。詳細は本記事では説明しないので、気になる方は what if book などを確認してみてください。

SUTVA: 「ある個人への処置が、別の個人のoutcomeに影響しない」ことと「処置にバリエーションがない」ことの仮定です。
Positivity: 共変量 $L$ で条件づけたときに、すべての処置 $a \in \mathcal{A}$ に対して、処置を受ける確率がゼロより大きいことを仮定します（特に断りがない限り、本記事では処置が二値のときを考えます）。
Consistency: 実際に観測されたoutcomeは、実際に行われた処置に対応するoutcomeと一致することを仮定します。
Conditional exchangeability: 共変量で条件づけたときに、Potential outcomesと処置変数は確率変数の意味で独立であるという仮定です。
No selection bias: 選択バイアスが存在しないことを仮定します。
No measurement error: 測定誤差が存在しないことを仮定します。

ちなみに、Conditional exchangeability を式で書くとこのようになります:

$Y^{(a)} \perp\!\!\!\perp A \mid L \tag{1}$

IPWを用いた因果効果の推定

本記事では、因果効果として、 Potential outcomes の期待値の差分による定義（Average Treatment Effect; ATE）を採用します。特に、処置が二値の場合は、以下のように定義されます:

\begin{eqnarray} ATE := \mathbb{E}[Y^{a=1}] - \mathbb{E}[Y^{a=0}] \tag{2} \end{eqnarray}

定義の詳細が気になる方は what if book の1章や2章あたりをご覧ください。

6つの条件の仮定のもとで、Potential outcome についての以下の等式が成立することが知られています。

\begin{eqnarray} \mathbb{E}[Y^{(a')}] = \mathbb{E} \Bigl[ \frac{I(A, a')Y}{f_{A \mid L}(a' \mid L)} \Bigr] \tag{3} \end{eqnarray}

等式の左側の期待値は Potential outcome についてとっていて、右側の期待値は $A, L, Y$ についてとっています。

証明は what if book の Technical Point 2.3 に書かれているので省略します。本記事とは少し表記が異なるので適宜読み替えてください。

とにかく、等式の右辺を Estimand*3 として何かしらの方法で推定することで、因果効果を推定できそうだということがわかりました。

Horvitz-Thompson estimator とは

what if book の Technical Point 12.1 に書かれているように、式 (2) の右辺を観測データを用いて以下のように推定（経験近似）したものは Horvitz-Thompson estimator と呼ばれています:

\begin{eqnarray} \hat{\mathbb{E}} \Bigl[ \frac{I(A, a')Y}{f_{A \mid L}(a' \mid L)} \Bigr] \tag{4} \end{eqnarray}

これは Potential outcome の期待値の不偏推定量かつ一致推定量になっています。

実際に計算するときには、傾向スコアのモデル $\hat{f}$ を推定した後に、以下のように $\hat{f}$ をプラグインして計算します。

\begin{eqnarray} \frac{1}{n}\sum_{i=1}^{n}\frac{I(A_i, a')Y_i}{\hat{f}_{A \mid L}[a' \mid L_i]} \tag{5} \end{eqnarray}

最後に、処置群と統制群で式 (4) の計算結果の差をとります。

ここまでが基本的な IPW による因果効果の推定になります。

以上を踏まえつつ、ここから先は、推定の改良を試みる手法について説明していきます。

Hajek (modified Horvitz-Thompson) estimator

ここでは、冒頭に述べた主張の二つ目の話、Hajek (modified Horvitz-Thompson) estimator について説明します。

Hajek estimator とは

what if book の Technical Point 12.1 を見ると、 Horvitz-Thompson estimator のすぐ下で、 modified Horvitz-Thompson estimator という以下のような推定量が紹介されているかと思います:

\begin{eqnarray} \frac{\hat{\mathbb{E}} \Bigl[ \frac{I(A, a')Y}{f_{A \mid L}(a' \mid L)} \Bigr]} {\hat{\mathbb{E}} \Bigl[ \frac{I(A, a')}{f_{A \mid L}(a' \mid L)} \Bigr]} \tag{6} \end{eqnarray}

どうやら、Horvitz-Thompson estimator の期待値の中から $Y$ を削除したものを分母につけているようです。

IPW の本質は Data augmentation だと思っているのですが、「観測データの個数 $n$ ではなく、実際にどのくらい augmentation したのかに相当する項で割った方がよさそう」という気持ちを感じます。

what if book では推定量の名前は紹介されていませんでしたが、この推定量には Hajek estimator という名前がついているらしいので、以降 Hajek estimator と呼ぶこととします。

Hajek estimatorの性質

Hajek estimator の性質についての本記事での立場は以下となります:

Hajek estimator は Potential outcome の期待値の不偏推定量とは限らないが、一致推定量ではある
Potential outcome の符号が等しいとき、 Horvitz-Thompson estimator より Hajek estimator の方が漸近分散が小さい
Hajek estimator は $Y$ が二値である場合の推定結果が 0 と 1 の間に収まる（直感的に理解しやすい）

まず、漸近分散以外の部分について、参考文献とともに議論を整理します。

漸近分散以外の部分の議論の整理 (bias of Hajek estimator)

以前の記事で紹介した星野先生の著書（調査観察データの統計科学）の P70 では、関連内容についての記述が存在します。

調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

作者:星野崇宏
発売日: 2009/07/29
メディア: 単行本

該当部分の記述を本記事の用語に翻訳すると、以下のようにまとめられます:

Horvitz-Thompson estimator は Potential outcome の期待値の不偏推定量であり、一致推定量でもある
Hajek estimator は Potential outcome の期待値の不偏推定量であり*4、一致推定量でもある
Horvitz-Thompson estimator より、 Hajek estimator の方が推定精度が高い

また、 what if book の Technical Point 12.1 では以下のように述べられています:

Positivity の仮定が成立するとき、 Hajek estimator は Potential outcome の期待値の不偏推定量である
Hajek estimator は $Y$ が二値である場合の推定結果が 0 と 1 の間に収まるので好ましい

一方、世界的に有名な今井先生の資料 (https://imai.fas.harvard.edu/teaching/files/weighting.pdf) では、以下のように述べられています:

Weights are normalized but no longer unbiased （Hajek estimator で重みは正規化できるけど、不偏ではなくなる）

すなわち、不偏性についての記述に矛盾が起きているように感じます。

詳しく見ていくと、what if book では不偏性の証明が与えられていません。

また、星野先生の著書には証明が与えられていますが、 Hajek estimator の不偏性を示すべき文脈で、 Horvitz-Thompson estimator の不偏性を示す式 (3.9) が書かれてしまっています:

\begin{eqnarray} \mathbb{E}\Bigl[\frac{AY}{\pi} \Bigr] = \mathbb{E}[Y^{a=1}] \tag{3.9} \end{eqnarray}

つまり、Hajek estimator の不偏性については、 what if book ならびに星野先生の著書では証明が与えられていないと考えられます。

証明を与える必要性を感じたため、似たような推定量の形を調べてみると、どうやら Ratio estimator というクラスとの関係性が深そうだとわかりました。

実際にRatio estimator の資料 (https://jkim.public.iastate.edu/teaching/book8.pdf) を読んでみると、以下のことが示せます:

Ratio estimator には一般に、式 (8.5) の形でバイアスが存在する
Hajek estimator は Ratio estimator で $x_i = 1$ とおいた場合に相当する

よって、本記事の中では「what if book や星野先生の著書の記述は誤りであり、 Hajek estimator には不偏性があるとは言えないのではないか」という立場をとります*5。

バイアスの具体的な形については、英語のブログ: Bias of the Hajek estimator (potential errors in Technical Point 12.1 of Causal inference: what if) - fullflu-english の後半に自分の計算結果を記載しておきました。

漸近分散の解析

星野先生の著書の中で「Horvitz-Thompson estimator より、 Hajek estimator の方が推定精度が高い」という主張がされていましたが、この部分について説明します。

著書の中ではこの主張の証明は行われていないため、Horvitz-Thompson estimator と Hajek estimator の漸近分散を比較することで、この主張に妥当性がありそうだということを示します。

ここから先はM推定量の理論解析を行うため、詳細に興味のある方のみご覧ください。

なお、傾向スコアについては真のモデルが given である場合を想定して計算を簡略化します*6。

Hajek estimator の漸近分散の整理

星野先生の著書の中では、処置が二値の場合の Hajek estimator の漸近分散の解析結果が与えられています。この議論の流れにのっとります。

漸近分散の解析においてM推定量の理論を使うのですが、こちらのスライドの27枚目以降の補足が非常にわかりやすいので、本記事でも参照させていただきます:

星野「調査観察データの統計科学」第3章 from Shuyo Nakatani

さて、まず星野先生の著書で与えられている Hajek estimator の漸近分散の主張（上記スライドの27枚目のスライドを復唱する形になり恐縮ですが）を本記事の Notation を使って整理します。

$\theta=(\theta_1, \theta_0)^{t}$ を真値 $\theta^{*}$ が $(\mathbb{E}[Y^{(a=1)}], \mathbb{E}[Y^{(a=0)}])$ である母数とし、以下の推定関数を考えます:

\begin{eqnarray} m(Y, \theta) := \Bigl(\frac{A}{\pi}(Y-\theta_1), \frac{1-A}{1-\pi}(Y-\theta_0) \Bigr)^{t} \tag{7} \end{eqnarray}

このとき、推定量 $\hat{\theta}$ の漸近分散（分散共分散行列） $Var(\theta^{*})$ は $J^{-1}K(J^{-1})^{t}$ になるという事実がM推定量についての一般論として知られています。ただし、

\begin{eqnarray} J(\theta^{*}) &:=& \mathbb{E}\Bigl[-\frac{\partial m(Y, \theta)}{\partial \theta^{t}} \Bigr]_{\theta=\theta^{*}}\\ K(\theta^{*}) &:=& \mathbb{E}[m(Y, \theta^{*}) m(Y, \theta^{*})^{t}] \tag{8} \end{eqnarray}

と定義しました。

これを計算すると、行列の要素は以下のようになります（ $J$ は単位行列になるので、 $K$ だけを考えることになります）:

\begin{eqnarray} Var(\theta^{*})_{1,1} &=& \mathbb{E}\Bigl[ \frac{(Y^{a=1}-\theta_1)^2}{\pi} \Bigr] \\ Var(\theta^{*})_{2,2} &=& \mathbb{E}\Bigl[ \frac{(Y^{a=0}-\theta_0)^2}{1-\pi} \Bigr] \\ Var(\theta^{*})_{1,2} &=& Var(\theta^{*})_{2,1} = 0 \tag{9} \end{eqnarray}

よって、平均因果効果の漸近分散は以下のように求まります:

\begin{eqnarray} Var_{Hajek}(\hat{\theta}_{1} - \hat{\theta}_{0}) &=& Var(\theta^{*})_{1,1} + Var(\theta^{*})_{2,2} - 2Var(\theta^{*})_{1,2}\\ &=& \mathbb{E}\Bigl[ \frac{(Y^{a=1}-\theta_1)^2}{\pi} + \frac{(Y^{a=0}-\theta_0)^2}{1-\pi} \Bigr] \tag{3.10改} \end{eqnarray}

これで復唱が終了したので、同様の方法で Horvitz-Thompson estimator の漸近分散の解析を行います。

Horvitz-Thompson estimator の漸近分散 (asymptotic variance of Horvitz-Thompson estimator)

Horvitz-Thompson estimator の漸近分散について、公開されている資料でわかりやすい解説がすぐには見つからなかったので、自分で解析してみます*7。

まず、このような推定関数を考えます:

\begin{eqnarray} m(Y, \theta) := \Bigl(\frac{AY}{\pi} -\theta_1, \frac{(1-A)Y}{1-\pi}-\theta_0 \Bigr)^{t} \tag{10} \end{eqnarray}

M推定量の枠組みで考えると、Hajek estimator の場合と同様に $J$ は単位行列になるので、 $K$ だけを考えればよくなります。

以下のように行列の要素を計算できます:

\begin{eqnarray} Var(\theta^{*})_{1,1} &=& \mathbb{E}\Bigl[ \frac{A^2{Y^{a=1}}^{2}}{\pi^2} + \theta_1^2 - 2 \frac{AY^{a=1}}{\pi}\theta_1 \Bigr] \\ &=& \mathbb{E}\Bigl[ \frac{{Y^{a=1}}^{2}}{\pi} + \theta_1^2 - 2Y^{a=1}\theta_1 \Bigr]\\ Var(\theta^{*})_{2,2} &=& \mathbb{E}\Bigl[ \frac{(1-A)^2{Y^{a=0}}^{2}}{(1-\pi)^2} + \theta_0^2 - 2 \frac{(1-A)Y^{a=0}}{1-\pi}\theta_0 \Bigr] \\ &=& \mathbb{E}\Bigl[ \frac{{Y^{a=0}}^{2}}{1-\pi} + \theta_0^2 - 2Y^{a=0}\theta_0 \Bigr]\\ Var(\theta^{*})_{1,2} &=& Var(\theta^{*})_{2,1}\\ &=& \mathbb{E}[\theta_1\theta_0 - Y^{a=1}\theta_0 - Y^{a=0}\theta_1] \tag{11} \end{eqnarray}

よって、 $\mathbb{E}[Y^{a=1}] = \theta_1$ であることを使うと、平均因果効果の漸近分散は以下のようにまとめられます:

\begin{eqnarray} Var_{HT}(\hat{\theta}_{1} - \hat{\theta}_{0}) &=& Var(\theta^{*})_{1,1} + Var(\theta^{*})_{2,2} - 2Var(\theta^{*})_{1,2}\\ &=& \mathbb{E} \Bigl[ \frac{{Y^{a=1}}^{2}}{\pi} + \frac{{Y^{a=0}}^{2}}{1-\pi} + (\theta_1 - \theta_0)^2 + 2Y^{a=1}\theta_0 + 2Y^{a=0}\theta_1 - 2Y^{a=1}\theta_1 - 2Y^{a=0}\theta_0 \Bigr] \\ &=& \mathbb{E} \Bigl[ \frac{{Y^{a=1}}^{2}}{\pi} + \frac{{Y^{a=0}}^{2}}{1-\pi} - (\theta_1 - \theta_0)^2 \Bigr] \tag{12} \end{eqnarray}

Hajek estimator と Horvitz-Thompson estimator の漸近分散の比較

最後に、式 (3.10改) と式 (12) を比較してみます。

\begin{eqnarray} Var_{HT}(\hat{\theta}_{1} - \hat{\theta}_{0}) - Var_{Hajek}(\hat{\theta}_{1} - \hat{\theta}_{0}) &=& \mathbb{E} \Bigl[ \frac{{Y^{a=1}}^{2}}{\pi} + \frac{{Y^{a=0}}^{2}}{1-\pi} - (\theta_1 - \theta_0)^2 \Bigr] - \Bigl(\mathbb{E}\Bigl[ \frac{(Y^{a=1}-\theta_1)^2}{\pi} + \frac{(Y^{a=0}-\theta_0)^2}{1-\pi} \Bigr] \Bigr)\\ &=& \mathbb{E}\Bigl[ \frac{{Y^{a=1}}^{2}}{\pi} - \frac{(Y^{a=1}-\theta_1)^2}{\pi} \Bigr] + \mathbb{E}\Bigl[ \frac{{Y^{a=0}}^{2}}{1-\pi} - \frac{(Y^{a=0}-\theta_0)^2}{1-\pi} \Bigr] - (\theta_1 - \theta_0)^2 \tag{13} \end{eqnarray}

のように計算できます。

2次のモーメント ≧ 分散の不等式を使うことにより、期待値が残っている部分については非負であることが示せるものの、最後の定数項がゼロ以下であるので、何かしらの不等式を示すにはもう少し展開する必要がありそうです。

ここで、さらに2次モーメントと分散の関係を明示的に利用することで、以下のように式を展開できると考えられます:

\begin{eqnarray} (13) &=& \mathbb{E}\Bigl[ \frac{\theta_1^{2}}{\pi} \Bigr] + \mathbb{E}\Bigl[ \frac{\theta_0^{2}}{1-\pi} \Bigr] - (\theta_1 - \theta_0)^2 \ \ \ (\because \text{Relationship between variance and moment}) \\ &>& \theta_1^2 + \theta_0^2 - (\theta_1 - \theta_0)^2 \ \ \ (\because 0 < \pi < 1) \\ &=& 2\theta_1\theta_0 \tag{14} \end{eqnarray}

よって、 Potential outcome の符号が等しければ右辺が0以上になり、「Horvitz-Thompson estimator より Hajek estimator の方が漸近分散が小さい」という主張が成立すると考えられます。

なお、ここで行った解析は専門家のレビューを通過しているわけではないため、以下の可能性が残っているということを示唆して終わりとします:

漸近分散の解析を改良すればより少ない仮定で不等式が示せる
漸近分散以外の観点で見ればより意味のある不等式が示せる
私の解析が誤っていて、このような不等式は示せない

終わりに

本記事では、因果推論で重要となる IPW について、いくつかの代表的な推定量の性質を解析しました。

一通り読んでいただいた方には、因果推論の著名な本の中でさらっと述べられていて盲目的に信じてしまいがちなところに意外な落とし穴があるかもしれない、ということをお伝えできたかと思います。

今後は、他の推定量の性質や、実験による主張の検証についての記事を書く予定をしています。

記事を読んでいただいた皆さんの豊かな因果推論ライフに貢献できれば幸いです。

*1:本記事では、6つの仮定 + Potential outcome の符号一致仮定のもとでの主張を考えます

*2:個人的には $Y^{(a=a')}$ を処置 $a'$ の Potential outcomeと書いた方が他の表記との整合性があると思うのですが、ここらへんは慣習の問題もありそうなのでよしなに読み替えてください

*3:Estimand の定義については what if book の10章をご覧ください

*4:「周辺平均の不偏推定量になる」と書かれています

*5:私の理解に誤りがあればご指摘いただけると幸いです

*6:多くの場合、傾向スコアモデルのパラメータ推定を含んだ形での解析を行う必要があるのですが、複雑になるので割愛しました。ご容赦ください

*7:漸近分散の解析に慣れていないので誤りがあればご指摘いただけると幸いです

果物をいっぱい食べたい

統計、機械学習周りの勉強や提案やOSS開発の記録

色々なIPWの整理〜 Causal inference: what if の12章を中心に〜 (1)