Marginal Structural Model のパラメータ推定に重み付き最小二乗法を用いることの正当性

本記事では、 Marginal Structural Model (MSM) のパラメータ推定に重み付き最小二乗法 (IP-Weighted Least Squared method; IP-WLS method) を用いることの正当性についての解釈を試みます。what if book の12章で紹介されている話の一部を深く掘り下げていきます。

記事のゴール

以下の3つの問いに答えることを目指します（自分が12章を読んだときに感じた疑問でした）:

1. 処置が2値の場合には、IP-WLS の推定量が Hajek estimator と等しいという主張が what if book P152 で書かれているが、これは本当か？
2. 処置が多値の場合にも、IP-WLS の推定量は Hajek estimator と等しいのか？
3. 処置が連続の場合にも、何かしら似たような解釈はできるのか？

執筆終了時点での私の回答としては以下になります:

1. 本当であるということを証明します。
2. 星野先生の本の議論を参考にしつつ、等しいことを証明します。
3. モーメント法・M推定量の文脈から、自分なりに正当化を試みます。

なお、本記事のタイトルの「正当性」という単語には、

• 特定（処置が離散）の場合に、既知の推定量と同じ結果が得られるという意味で手法の正当化ができる
• 処置が連続の場合も、離散の場合の拡張として捉えると、一定の妥当性があるように思える

といった気持ちを込めています。

新規性

モーメント法・M推定量の文脈で処置が連続の場合のお気持ちを書いてある日本語記事は見つけられませんでしたので、それを新規性とさせていただきます（正確性は保証できないです）。

前提知識

• 過去の自分の記事 における、 Hajek estimator のモーメント法による定式化の雰囲気がわかっていること

what if book のMSMのおさらい

what if book 12章については過去の記事でも紹介しましたが、MSMやそのパラメータ推定方法については述べていなかったので、軽くおさらいします。

fullflu.hatenablog.com

モデリング

MSMとは、potential outcome の期待値を（周辺化して）直接予測するようなモデルを立てることを意味します。

過去の記事でも紹介した exchangeability などの仮定のもとで、そのモデルのパラメータの推定結果を用いて因果効果を良い感じに*1推定することが可能です。

処置が2値のとき、以下のような MSM のモデルが 12.4 節で紹介されています:

\begin{eqnarray} \mathbb{E}[Y^{a}] = \beta_0 + \beta_1 a \end{eqnarray}

個人的には を使って以下のように書いた方が綺麗だと思うのですが、必要に迫られるまでは慣習にしたがって を使わずに書きます。

\begin{eqnarray} \mathbb{E}[Y^{a=a'}] = \beta_0 + \beta_1 a' \end{eqnarray}

このとき、以下のように因果効果を計算できます:

\begin{eqnarray} \mathbb{E}[Y^{a=1}] &=& \beta_0 + \beta_1\\ \mathbb{E}[Y^{a=0}] &=& \beta_0\\ ATE &=& \beta_1 \end{eqnarray}

また、処理が連続のとき、以下のようなモデルが紹介されています:

\begin{eqnarray} \mathbb{E}[Y^{a}] = \beta_0 + \beta_1 a + \beta_2 a^2 \end{eqnarray}

あくまでこれは2次式までを考慮した一例であり、実際にはより複雑なモデルを立てることもありえるということに注意してください。

パラメータ推定

MSM のパラメータ推定にあたって、 IPW によって重み付けした最小二乗法、通称 IP-WLS を用いる手法が 12.4 節で紹介されています。

実装は、Python の場合は statsmodels、R の場合は lm や gee を使うことになるでしょう。

パラメータの分散を推定する際には、厳密に解析するか、 bootstrap を使うか、サンドウィッチ分散を使うかの三択になるとも紹介されています。

また、P152 の左側には、 IP-WLS の結果が Hajek estimator と一致すると書かれています。

処置が2値の場合、本当に IP-WLS は Hajek estimator と一致するのか？

まず、一つ目の疑問に答えていきます。愚直に解析します。

個人 が処置 を受ける確率 ] を と表記すると、 IP-WLS は

\begin{eqnarray} L = \sum_{i=1}^{n} \Biggl[ \frac{A_i}{\pi_i} \bigl(Y_i - \beta_0 - \beta_1 A_i \bigr)^2 + \frac{1-A_i}{1-\pi_i} \bigl(Y_i - \beta_0 - \beta_1 A_i \bigr)^2 \Biggr] \end{eqnarray}

を最小化する問題を考えればよいことになります。偏微分を行った式

\begin{eqnarray} \frac{\partial{L}}{\partial{\beta_1}} &=& 2\sum_{i=1}^{n} \Biggl[ \frac{A_i}{\pi_i}\bigl(Y_i - \beta_0 - \beta_1 \bigr)(-A_i) + \frac{1-A_i}{1-\pi_i} \bigl(Y_i - \beta_0 - \beta_1 A_i \bigr)(-A_i) \Biggr] \\ &=& 2\sum_{i=1}^{n} \Biggl[ \frac{-A_i}{\pi_i}\bigl(Y_i - \beta_0 - \beta_1 \bigr) \Biggr] \\ &=& 0 \end{eqnarray}

を整理すると、] の推定量は

\begin{eqnarray} \sum_{i=1}^{n} \Biggl[ \frac{A_i}{\pi_i}\bigl(Y_i - \hat{\theta_1} \bigr) \Biggr] &=& 0 \end{eqnarray}

の解となります。

これを踏まえてさらに、

\begin{eqnarray} \frac{\partial{L}}{\partial{\beta_0}} &=& 2\sum_{i=1}^{n} \Biggl[ \frac{-A_i}{\pi_i}\bigl(Y_i - \beta_0 - \beta_1 \bigr) - \frac{1-A_i}{1-\pi_i} \bigl(Y_i - \beta_0 - \beta_1 A_i \bigr) \Biggr] \\ &=& -2\sum_{i=1}^{n} \Biggl[ \frac{1-A_i}{1-\pi_i} \bigl(Y_i - \beta_0 \bigr) \Biggr] \\ &=& 0 \end{eqnarray}

を整理すると、] の推定量は

\begin{eqnarray} \sum_{i=1}^{n} \Biggl[ \frac{1-A_i}{1-\pi_i} \bigl(Y_i - \hat{\theta_0} \bigr) \Biggr] &=& 0 \end{eqnarray}

の解になります。

これらの式を、前提知識の章のリンクで説明した Hajek estimator のモーメント法による定式化と比較すると、全く同じ形になっているということが確認できると思います。

したがって、 処置が2値のとき、 IPW で重み付けをした最小二乗法の結果は、確かに Hajek estimator と等しいことが示せました。

処置が多値のときはどうなるか？

what if book には処置が多値の場合については説明がされていないので、2値の場合をよしなに拡張して定式化を行います。

ある処置のレベル に対して、周辺構造モデルを以下のように表現します:

\begin{eqnarray} \mathbb{E}[Y^{a=j}] = \beta_j \end{eqnarray}

個人 が処置 を受ける確率 ] を と表記すると、 IP-WLS では

\begin{eqnarray} L = \sum_{i=1}^{n} \sum_{j=1}^{J} \frac{\delta(A_i-j)}{\pi_{ij}} \bigl(Y_i - \beta_j \bigr)^2 \tag{1} \end{eqnarray}

を最小化する問題を考えればよいことになります。 はディラックのデルタ関数です。

\begin{eqnarray} \frac{\partial{L}}{\partial{\beta_j}} &=& -2\sum_{i=1}^{n} \frac{\delta(A_i-j)}{\pi_{ij}} \bigl(Y_i - \beta_j \bigr)\\ &=& 0 \end{eqnarray}

を整理すると、] の推定量は

\begin{eqnarray} \sum_{i=1}^{n} \frac{\delta(A_i-j)}{\pi_{ij}} \bigl(Y_i - \hat{\theta_j} \bigr) &=& 0 \end{eqnarray}

の解になります。

多値の場合の Hajek estimator の定義は what if book には載っていなかったはずですが、この推定量は、2値の場合の Hajek estimator を素朴に拡張した式になっていることが確認できるかと思います。

M推定量の文脈による正当化（処置が離散の場合）

これまで、Hajek estimator との対応によって MSM の正当化を試みてきましたが、ここではM推定量の立場から改めて整理することを試みます。

以下の星野本の P81-82 を参考にします。

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

まず、天下り式に、処置 に対する potential outcome とパラメータを引数にとった非負のスコア関数 というものを導入します。

二乗誤差の最小化を考える場合、以下のようなスコア関数を考えることになるでしょう:

\begin{eqnarray} m_j(Y^{a=j}, \theta_j) &:=& -(Y^{a=j} - \theta_j)^2 \end{eqnarray}

さらに、以下のようなスコア関数の和の期待値を最大化したい、と考えたとします:

\begin{eqnarray} \mathbb{E}_{\tilde{Y}} \Biggl[ \sum_{j=1}^{J} m_j(Y^{a=j}, \theta_j) \Biggr] \tag{2} \end{eqnarray}

ここでは、期待値は potential outcomes についてとっています。

この式からは、すべての処置の水準の重要度を等しく扱った上で、処置ごとの potential outcome をできる限り正しく推定したい という気持ちを感じることができますし、この気持ち自体をある意味で「自然」なものだとして正当化することは可能だと考えられます（私個人の感覚ですが）。

しかし残念なことに、式 (2) には counterfactual outcome が含まれているので、そのまま経験近似を行うことが不可能です。

そこで、式 (2) と期待値の意味で等しくて（不偏で）、かつ経験近似が容易に行える ような別の式を考えたくなります。

そのような性質を持つ式の一つが以下のような形のもの（期待値の中身）です。不偏性は conditional exchangeability の仮定によって示せます。

\begin{eqnarray} \mathbb{E}_{\tilde{Y}, L, A} \Biggl[ \sum_{j=1}^{J} \frac{\delta(A-j)}{\pi_{j}} m_j(Y^{a=j}, \theta_j) \Biggr] \tag{3} \end{eqnarray}

そして、consistency の仮定を用いてこの式 (3) を経験近似したものが、式 (1) のIP-WLS であると解釈できます。

また、式 (3) をパラメータで微分してから、微分と積分が交換できることを仮定して経験近似をすると、M推定量を考えていることになるような気持ちになれます。

改めてまとめます。

• 「すべての処置の水準の重要度を等しく扱った上で、処置ごとの potential outcome をできる限り正しく推定したい」という気持ちを式 (2) として表現することを認めた上で
• 「式 (2) と期待値の意味で等しくて、かつ経験近似が容易に行える」ような便利な式を考えると
• 式 (3) の形は自然と導出され、推定の正しさを二乗誤差で測定することに同意すれば、その経験近似により IP-WLS の式が得られる（& M推定量として解釈できる）

というのがここでの主張になります。

一見突拍子もないように感じた IP-WLS が自然な導出により得られる、ということが共有できたかと思います。

M推定量の文脈による正当化（処置が連続の場合）

最後に、処置が連続の場合についても直感的な解釈を与えていきます。

ここから先の話はどの本にも載っていなかったので、自分の感覚を信じて論理を組み立てていきます。

まず、先ほどと同様に式 (2) のようなものを考えてみると、処置が連続の場合、処置について和ではなく積分を考える必要が出てくることがわかります。

また、すべての処置の水準の重要度を等しく扱おうとすると、処置の区間が有限でない場合積分値が無限に飛んでいってしまいます。

そこで、処置の密度関数 によって重み付けを行った上で、処置ごとの potential outcome をできる限り正しく推定したい という気持ちで式 (2) を修正してみます:

\begin{eqnarray} \mathbb{E}_{\tilde{Y}} \Biggl[ \int_{a' \in \mathcal{A}} f(a') m(Y^{a=a'}, \theta, a') da' \Biggr] \tag{4} \end{eqnarray}

なお、ここから先は potential outcome を扱うときに を使った表記に切り替えます。

また、スコア関数 の引数に処置の情報を含むようにします。

式 (4) を見てみると、当然 potential outcomes は観測できないので、この最大化をそのまま考えるのは困難であり、代理の式を探したくなります。

処置が離散の場合を参考にして、 をかけて期待値をとってみると、

\begin{eqnarray} \mathbb{E}_{\tilde{Y}, L, A} \Biggl[ \int_{a' \in \mathcal{A}} f(a') m(Y^{a=a'}, \theta, a') \frac{\delta(A - a')}{f(a' \mid L)} da' \Biggr] &=& \int_{a' \in \mathcal{A}} \mathbb{E}_{L}\Biggl[ \mathbb{E}_{A \mid L} \bigl[\frac{\delta(A - a')}{f(a' \mid L)} \bigr] \mathbb{E}_{Y^{a=a'}}[m(Y^{a=a'}, \theta, a')] \Biggr] f(a')da' \\ &=& \int_{a' \in \mathcal{A}} \Biggl[\mathbb{E}_{Y^{a=a'}}[m(Y^{a=a'}, \theta, a')] \Biggr] f(a')da' \tag{5} \end{eqnarray}

となり、式 (4) に対する不偏性を持つことが確認できます*2。

したがって、これを経験近似した以下の式 (6) を最大化するという操作には、処置が離散の場合と同様にM推定量の文脈からその正当性を与えられそうです。

\begin{eqnarray} \sum_{i=1}^{n} \Biggl[ \int_{a' \in \mathcal{A}} \biggl( \frac{\delta(A_i - a') f(a')} {f(a' \mid L_i)} m(Y_{i}^{a=a'}, \theta, a') \biggr) da' \Biggr] &=& \sum_{i=1}^{n} \Biggl[ \frac{f(A_i)} {f(A_i \mid L_i)} m(Y_{i}^{a=A_i}, \theta, A_i) \Biggr]\\ &=& \sum_{i=1}^{n} \Biggl[ \frac{f(A_i)} {f(A_i \mid L_i)} m(Y_{i}, \theta, A_i) \Biggr] \tag{6} \end{eqnarray}

個人的には、不偏性を示すときと経験近似をするときとで、デルタ関数の処理が異なってくるというのがエモポイントです。

そして、式 (6) のスコア関数に二乗誤差を用いると stabilized weights を用いた IP-WLS そのものになります。

以上の議論により、処置が連続の場合の IP-WLS についても、M推定量の文脈での正当化を行うことができ、その意味で「自然な」推定手法だと解釈できると考えられます。

おわりに

本記事では、 IP-WLS で MSM のパラメータ推定を行うことの正当性について、Hajek estimator や M推定量との対応をとりながら紹介を行いました。

少し記述が回りくどくなってしまいましたが、モヤモヤを持たれている方の助けになれば幸いです。

誤りを見つけられた方は、そっと教えていただけると幸いです。

*1:一致性を持ち、漸近分散が小さいという意味で良い感じです

*2:密度関数の notation が崩れてますが、この辺は what if book を見て察していただけると幸いです