本記事では、 g-estimation でパラメータ推定を行う際に推定方程式を用いる方法について説明します。what if book の14章で紹介されている話の一部を深く掘り下げていきます。

記事のゴール

以下の疑問に答えることをゴールにします（what if 勉強会でも疑問として紹介されていました）:

what if book の本文で紹介されているロジットモデルの係数検定の話と、Technical Point 14.2 で紹介されている推定方程式の話が等価になる理由がわからない
なぜ推定方程式の形で表現したいのかわからない

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

「 $H(\psi)$ と $A$ の相関がないような制約をかけたい」という気持ちを共分散がゼロになるという式で表現して経験近似するとそうなるから
(a) 「処置モデルの係数の検定」という制約に縛られる必要がなくなり、処置モデルについてより柔軟なケースを扱えるから (b) 推定方程式の一般論を用いて、推定量の漸近的な性質について語れるから

新規性

g-estimation と推定方程式の関係性を日本語で一つ一つ解説している記事は見つけられませんでしたので、それを新規性とさせていただきます（正確性は保証できないです）。

前提知識

what if book の14章に軽く目を通していること（rank preservation, conditional exchangeability, consistency などは記事の中で説明しないです）
共分散や検定論といった基本的な統計の概念がわかっていること

what if book の g-estimation のおさらい

what if book 14章では、 g-estimation の定義はふわっとしか与えられていない認識です。

章の構成としては、ざっくり以下のようになっています:

structural nested mean model (SMM) について導入する
rank preservation 仮定を紹介する*1
rank preservation と conditional exchangeability と consistency が成立していたと仮定すると、処置モデルの係数の検定結果によって因果効果を推定できることを示す
こうした推定のための手続きを g-estimation と呼んでいる

ひとまずは、SMM のパラメータ推定を行うためのトリッキーな手法なんだ、くらいに思っていただけるとよいと思います。

以下、記事のゴールに対して必要最低限と思われる範囲で、 g-estimation の説明を行います。

SMMの導入

SMMは、potential outcome の差分の期待値をモデリングする手法のことを指します。

例えば、以下のようなモデルを立てて、パラメータ $\beta_1$ を平均因果効果として推定することが目標になります:

\begin{eqnarray} \mathbb{E}[Y^{a} - Y^{a=0} \mid L] = \beta_1 a \end{eqnarray}

rank preservation と conditional exchangeability の仮定のもとで、potential outcome の差分が以下のような effect modification のない線形モデルで表現できると仮定します:

\begin{eqnarray} Y^{a} = Y^{a=0} + \psi_1 a \tag{1} \end{eqnarray}

なお、rank preservation の仮定によって $\beta_1 = \psi_1$ となりますので、以降はパラメータ $\psi_1$ を求めることにフォーカスして考えます。

g-estimation

g-estimation の一歩目として、consistency の仮定から、観測データを用いて式(1)を書き換えます:

\begin{eqnarray} Y^{a=0} = Y - \psi_1 A \tag{2} \end{eqnarray}

次に、この式をもう少し抽象化して考えます。

まず、 $H(\psi) = Y - \psi A$ とすると、真の $\psi_1$ に対して $H(\psi_1) = Y^{a=0}$ になることは定義より自明です。

また、 conditional exchangeability の仮定により、 $Y^{a=0}$ と $A$ は $L$ を与えたときに独立になります。

よって、真の $\psi_1$ に対して $H(\psi_1)$ と $A$ は $L$ を与えたときに独立になります。

逆に言えば、 $\psi$ を色々入れ替えながら「 $H(\psi)$ と $A$ が $L$ に対して条件つき独立になるような $\psi$ を見つけられれば、それが真の $\psi_1$ の候補になる（必要条件）」ことがわかります。

さらに、処置モデルがロジットモデルで特定できていると仮定すると、この条件つき独立性が成り立つとき、 $H(\psi)$ についての係数 $\alpha_1$ がゼロになります。処置モデルの例は以下です:

\begin{eqnarray} logit(\Pr[A=1 \mid H(\psi), L]) = \alpha_0 + \alpha_1H(\psi) + \alpha_2L \tag{3} \end{eqnarray}

すなわち、以下の手順で因果効果の候補を推定します:

$\psi$ として色々な値を用意する
それぞれの値に対して式(3)のモデルを学習させ、係数 $\alpha_1$ が0であるという帰無仮説を検定する
検定結果が有意にならなかった $\psi$ を解候補として出力する

実装時には計算機のパワーでグリッドサーチすることが多いです。

Technical Point 14.2: 推定方程式の紹介

簡単のため、censoringは無視して考えます。

Technical Point 14.2 では、以下の式(4)の推定方程式の解を求めることと、14章でやってきた g-estimation の手続きは等価だと述べられています:

\begin{eqnarray} \sum_{i=1}^{n} H_i(\psi) (A_i - \mathbb{E}[A \mid L_i]) \tag{4} \end{eqnarray}

ここで、 $H_i(\psi)$ はi番目のサンプルの観測値によって計算された $H(\psi)$ の値です。

期待値の部分は処置モデルを推定して予測すればよいようです。

なぜ推定方程式が出てくるのか

ここからが本記事の本題です。

式の導出

とは言っても考え方は非常にシンプルで、処置のロジットモデルのことは忘れて、「 $H(\psi)$ と $A$ が独立になるような $\psi$ を求めたい」という原点に立ち帰ればよいです。

この独立性の条件を必要条件に緩めると「 $H(\psi)$ と $A$ が無相関になるような $\psi$ を求めたい」になり、その条件を満たす $\psi$ の集合が解候補になると言えます。

無相関ということは相関係数がゼロということであり、すなわち共分散がゼロということになります。

$H(\psi)$ と $A$ の共分散がゼロであることは、以下の式(5)で表現できます。

\begin{eqnarray} \mathbb{E} \biggl[ \bigl( H(\psi) - \mathbb{E}[H(\psi)\mid L]\bigr) \bigl( A - \mathbb{E}[A \mid L] \bigr) \biggr] &=& \mathbb{E} \biggl[ H(\psi) \bigl( A - \mathbb{E}[A \mid L] \bigr) \biggr] \\ &=& 0 \tag{5} \end{eqnarray}

よって、結論としては、冒頭に書いた通り、「 $H(\psi)$ と $A$ の相関がないような制約をかけたい」という気持ちを共分散がゼロになるという式で表現して経験近似したら推定方程式(4)が出てきます。

解析解

式(2)の形でモデルを立てていた場合、推定方程式(5)にHの式を代入して移項と四則演算をするだけで、解 $\psi$ を陽に得られることがわかります。

Hajek estimator と似たような形の解になりますね:

\begin{eqnarray} \hat{\psi} = \frac{\sum_{i=1}^{n} Y_i (A_i - \mathbb{E}[A \mid L_i])} {\sum_{i=1}^{n} A_i (A_i - \mathbb{E}[A \mid L_i])} \tag{6} \end{eqnarray}

推定方程式として扱うことの嬉しさ

推定方程式では、ロジットモデルの係数というものを考えていないので、処置モデルの形状については明らかにagnosticです。

よって、推定方程式の形で議論することにより、現実にあったより柔軟なモデルを扱うことが可能、というメリットがあると考えます。

また、Technical Point 14.2 のように、censoringを考慮した定式化にも容易に拡張できたり、推定方程式の $H(\psi)$ の部分を $H(\psi) - \mathbb{E}[H(\psi)\mid L$ ] によって置き換えることで doubly robust な推定結果を得られたりするようです（証明は読んでいません）。