果物をいっぱい食べたい

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

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

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

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

20190411182013

目次はこちらです

記事のゴール

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

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

  1. 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 の推定である
  2. Hajek (modified Horvitz-Thompson) estimator は不偏推定量ではないが、いくつかの仮定*1のもとで Horvitz-Thompson estimator よりも漸近分散が小さい
  3. Marginal Structural Model のパラメータ推定を行う際に、逆確率で重み付けた最小二乗法を用いることが多いが、これは Hajek estimator の推定に相当する
  4. Stabilized weights は Hajek estimator を推定するときに有用かもしれないが、他の推定量を推定するときには注意が必要 (特に、Doubly robust では適切な使用方法がわからない)

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

新規性

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

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

  • 因果推論についての有名な資料間での、推定量についての記述の相違を整理した
  • 「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 では、関連内容についての記述が存在します。

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

  • 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枚目以降の補足が非常にわかりやすいので、本記事でも参照させていただきます:

さて、まず星野先生の著書で与えられている 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:漸近分散の解析に慣れていないので誤りがあればご指摘いただけると幸いです

調査観察データの統計科学 3.1の行間メモ

星野先生によって書かれた調査観察データの統計科学(以下本書と呼ぶ)を読んでいて、3.1節の式の導出に少し困ったので、メモを残します。 脚注に書いてある通り、厳密な証明は見つけられていないので、知見のある方はご指摘いただけると幸いです。

目次はこちらです:

3.1節の概要と本記事の目的

3.1節でやることはざっくり以下の三つです。

  1. バランシングスコアという概念を定義
  2. バランシングスコアを条件づけたときに、割り当て  z \in \{0, 1\}潜在的な結果変数  y_1, y_0 が独立になることを示す
  3. バランシングスコアの関数として傾向スコアを紹介

本記事では、バランシングスコアの必要条件についての記述、ならびにバランシングスコアを条件づけたときの独立性について補足します。

バランシングスコアの定義

共変量を  x とします。本書P60によると、バランシングスコア  b(x) とは、

それを条件付けすることにより、共変量と割り当てが独立になるような「共変量の関数」である

と定義されています。数式で書くと

 x \perp\!\!\!\perp z \mid b(x) \tag{3.1}

を満たす関数のことです。

バランシングスコアの必要条件

本書では、バランシングスコアの必要条件として、関数  g を使って  p(z=1 \mid x) = g(b(x)) と表現できることを挙げています。 その根拠として、以下の式の三つ目の等号成立条件において、上記のような表現が成立することが必要だから、というように書かれています。

\begin{eqnarray} p(z=1\mid b(x)) &=& \int p(z=1\mid x, b(x)) p(x\mid b(x)) dx \\ &=& \mathbb{E}_{x\mid b(x)} [p(z=1\mid x) ] \\ &=& \mathbb{E}_{x\mid b(x), \ p(z=1\mid x)} [p(z=1\mid x) ] \\ &=& p(z=1\mid x)\\ &=& p(z=1\mid x, b(x)) \tag{3.2} \end{eqnarray}

結論としては、この説明は必要条件ではなく十分条件の証明になっていると考えられます *1

本書の中では細かい式変形の行間が省略されているため、解説を試みます。

まず、一つ目の等式は、  x で条件つき確率を周辺化しています。こちらは周辺化の定義から自明かと思います。

次に、二つ目の等式は、  b(x) x の関数であることにより、  b(x) を条件づけから外せることを利用します *2

三つ目の等式は、もし  p(z=1 \mid x) = g(b(x)) と表現できたならば、二つ目の等式の条件づけを外したときと同じロジックで  p(z=1 \mid x) を条件づけから外せることを利用します。等式の左側が条件づけから外れた式で、右側が条件づけられた式になります。

四つ目の等式は、  p(z=1 \mid x) を条件づけた状態で  p(z=1 \mid x) の期待値を取っているので、期待値を外せる(であろう)ことを利用します*3

最後の五つ目の等式は、二つ目の等式を再び適用します。

三つ目の等式以外は条件つき確率の定義から計算できるであろうことがわかったので、三つ目の等式が十分条件になることがわかるかと思います。

バランシングスコアを条件づけたときの独立性

ここでは、バランシングスコアを条件づけたときに、割り当て  z潜在的な結果変数  y_1, y_0 が独立になること、すなわち

 (y_1, y_0) \perp\!\!\!\perp z \mid b(x) \\\tag{3.3}

を示します。

まず、先ほどと同様のやり方で、条件つき確率を展開します。

\begin{eqnarray} p(z=1\mid y_1, y_0, b(x)) &=& \int p(z=1\mid y_1, y_0, x, b(x)) p(x\mid y_1, y_0, b(x)) dx \\ &=& \mathbb{E}_{x\mid y_1, y_0, b(x)} [p(z=1\mid y_1, y_0, x)] \tag{a} \end{eqnarray}

一つ目の等式は周辺化、二つ目の等式は  b(x) x の関数であることを利用しています。

ここで、本書の2.5節で定義された"強く無視できる割り当て"条件が成立している、すなわち

 (y_1, y_0) \perp\!\!\!\perp z \mid x \\ \tag{2.15}

が成立していると仮定すると、  p(z=1\mid y_1, y_0, x) = p(z=1\mid x) なので、

\begin{eqnarray} \mathbb{E}_{x\mid y_1, y_0, b(x)} [p(z=1\mid y_1, y_0, x)] &=& \mathbb{E}_{x\mid y_1, y_0, b(x)} [p(z=1\mid x)] \tag{b} \end{eqnarray}

となります。

最後に、バランシングスコアの必要条件 *4 によって  p(z=1 \mid x) = g(b(x)) と表現できること、ならびに  b(x) で条件づけた確率で  g(b(x)) の期待値をとるとき期待値は外せる(であろう)ことを利用すると、

\begin{eqnarray} \mathbb{E}_{x\mid y_1, y_0, b(x)} [p(z=1\mid x)] &=& p(z=1\mid x) \tag{c} \end{eqnarray}

となり、式(3.2)と合わせると以下が言えます。

\begin{eqnarray} p(z=1\mid y_1, y_0, b(x)) &=& p(z=1\mid x) &=& p(z=1\mid b(x)) \tag{d} \end{eqnarray}

これにより、式(3.3) が成立します。

最後に

一つ一つ式を追っていくと、これ自明っぽいけど本当に自明なのかな?と迷うところがいくつか出てきて、鍛錬が足りないなぁと実感しました。 期待値の式変形を一つ一つ解説した日本語資料は簡単には見つけられなかったので、本記事が何かしらの助けになると幸いです。

*1:こちらのブログでも指摘がされています: 「調査観察データの統計科学」3.1章 傾向スコアの数式メモ(前半) - 木曜不足

*2:厳密に証明した資料は見つからなかったので正しさは保証できないですが、ここではこれが成り立つとします

*3:これも厳密な証明は見つけられていないです。適切な資料をご存知の方はご教授いただけると幸いです

*4:式(3.2)では十分条件しか示せなかったのですが、ここでは必要条件が成立することを認めるとします。気になる方は紹介されている論文を辿ってください