はじめに
東北大学/株式会社Nospareの石原です．前回の記事に引き続き，共変量を用いた平均処置効果 (average treatment effect; ATE) の推定方法について紹介します．本記事では，inverse probability weighting 法 と doubly robust 法 という2つの推定方法を紹介します．
設定
前回の記事と同様に，処置を受けたときに実現するであろう潜在的な結果変数を $Y_{i}(1)$，処置を受けなかったときに実現するであろう潜在的な結果変数を $Y_{i}(0)$ と表記し，$D_i$ を処置を受けたかどうかを表すダミー変数とします．このとき，観測できる結果変数 $Y_i$ は
Y_i = D_i Y_i(1) + (1-D_i) Y_i(0)
と書くことができ，$\tau_{\text{ATE}} = E[Y_i(1)-Y_i(0)]$ が ATE となります．また，観測できる共変量を $X_i$ とし，条件付き独立の仮定
\newcommand{\indep}{\perp\hspace{-1em}\perp} D_i \ \indep \ (Y_i(0), Y_i(1)) \ | \ X_i \hspace{1in} (1)
が成り立つとします．さらに，条件付き独立の仮定に加えて，オーバーラップの仮定
0 < P(D_i=1|X_i) < 1 \hspace{1in} (2)
も成り立つとします．ここで，$p(X_i) = P(D_i = 1 | X_i)$ は傾向スコア (propensity score) と呼ばれています．
Inverse probability weighting 法
傾向スコアを用いることで，regression adjustment 法のときに紹介したものとは異なる方法で ATE を識別することができます．$D_i Y_i = D_i Y_i(1)$なので，(1) と (2) の仮定から
\begin{align}
E\left[ \frac{D_i Y_i}{p(X_i)} \right] &= E\left[ \frac{E[D_i Y_i|X_i]}{p(X_i)} \right] = E\left[ \frac{E[D_i Y_i(1)|X_i]}{p(X_i)} \right] \\
&= E\left[ \frac{E[D_i|X_i] \cdot E[Y_i(1)|X_i]}{p(X_i)} \right] = E\left[ \frac{p(X_i) \cdot E[Y_i(1)|X_i]}{p(X_i)} \right] \\
&= E[E[Y_i(1)|X_i]] = E[Y_i(1)]
\end{align}
が成り立ちます．同様に，$E\left[ \frac{(1-D_i) Y_i}{1-p(X_i)} \right] = E[Y_i(0)]$ も成り立ちます．したがって，
E\left[ \frac{D_i Y_i}{p(X_i)} \right] - E\left[ \frac{(1-D_i) Y_i}{1-p(X_i)} \right] = \tau_{\text{ATE}}\hspace{1in} (3)
となります．以上から，傾向スコアを用いることで ATE が識別できることが分かります．
もし $p(X_i)$ が分かれば，(3) 式の左辺は観測できる変数から計算することができます．そこで，まず傾向スコア $p(X_i)$ を推定し，$p(X_i)$ を推定値に置き換えて (3) 式の左辺の標本対応を計算すれば ATE を推定することができると考えられます．この推定方法は inverse probability weighting (IPW) 法 と呼ばれています．具体的には，次のような手順で推定量を求めます：
- パラメトリックまたはノンパラメトリックな手法で $p(x) = P(D_i=1|X_i=x)$ を推定し，その推定値 $\hat{p}(x)$ を構成する．
- 手順1で求めた $\hat{p}(x)$ を用いて，以下の推定量で ATE を推定する：
\hat{\tau}_{IPW} = \frac{1}{n} \sum_{i=1}^n \left\{ \frac{D_i Y_i}{\hat{p}(X_i)} \right\} - \frac{1}{n} \sum_{i=1}^n \left\{ \frac{(1-D_i) Y_i}{1-\hat{p}(X_i)} \right\}
このようにして得られる推定量が IPW 推定量です．多くの場合，手順1の傾向スコアの推定には，ロジットモデルやプロビットモデルが用いられています．
Doubly robust 法
regression adjustment 推定量と IPW 推定量を構成するには，それぞれ条件付き期待値関数 $m_d(x) = E[Y_i|D_i=d,X_i=x]$ と傾向スコア $p(x) = P(D_i=1|X_i=x)$ の推定値が必要となります．ノンパラメトリックな手法を使えばモデルの特定化の誤りはありませんが，$X_i$ の次元が大きい場合にはノンパラメトリックな手法は不安定になります．今から紹介する doubly robust 法 は，条件付き期待値関数と傾向スコアのどちらかの推定が失敗していても ATE を推定することができる推定方法です．
条件付き独立の仮定 (1) とオーバーラップの仮定 (2) の下で，次の式が成り立ちます：
\begin{align}
\tau_{\text{ATE}} &= E\left[ \frac{D_i Y_i}{p(X_i)} + \left( 1 - \frac{D_i}{p(X_i)} \right) m_1(X_i) \right] \nonumber \\
& \hspace{0.5in} - E\left[ \frac{(1-D_i) Y_i}{1-p(X_i)} + \left( 1 - \frac{1-D_i}{1-p(X_i)} \right) m_0(X_i) \right] \hspace{1in} (4)
\end{align}
したがって，$m_d(x)$ の推定値を $\hat{m}_d(x)$，$p(x)$ の推定値を $\hat{p}(x)$ とすると，
\begin{align}
\hat{\tau}_{DR} &= \frac{1}{n} \sum_{i=1}^n \left\{ \frac{D_i Y_i}{\hat{p}(X_i)} + \left( 1 - \frac{D_i}{\hat{p}(X_i)} \right) \hat{m}_1(X_i) \right\} \nonumber \\
& \hspace{0.5in} - \frac{1}{n} \sum_{i=1}^n \left\{ \frac{(1-D_i) Y_i}{1-\hat{p}(X_i)} + \left( 1 - \frac{1-D_i}{1-\hat{p}(X_i)} \right) \hat{m}_0(X_i) \right\}
\end{align}
で ATE を推定することができると考えられます．これが doubly robust 推定量です．
上で紹介したように，条件付き期待値関数 $m_d(x)$ の推定と傾向スコア $p(x)$ の推定のどちらかが失敗しても ATE を推定することができます．そのため，二重頑健 (doubly robust) 推定量と呼ばれています．まず，条件付き期待値関数 $m_d(x)$ の推定が失敗していてもよいことを示すために，(4) 式の右辺第1項の $m_1(x)$ を間違った関数 $\tilde{m}_1(x)$ と置き換えた場合を考えます．このとき，$p(X_i) = E[D_i|X_i]$ から
\begin{align}
E\left[ \frac{D_i Y_i}{p(X_i)} + \left( 1 - \frac{D_i}{p(X_i)} \right) \tilde{m}_1(X_i) \right] &= E\left[ \frac{D_i Y_i}{p(X_i)} \right] + E \left[ \left( 1 - \frac{D_i}{p(X_i)} \right) \tilde{m}_1(X_i) \right] \\
&= E[Y_i(1)] + E \left[ E \left[ \left( \frac{p(X_i)-D_i}{p(X_i)} \right) \tilde{m}_1(X_i) \Big| X_i \right] \right] \\
&= E[Y_i(1)] + E \left[ \left( \frac{p(X_i)-E[D_i|X_i]}{p(X_i)} \right) \tilde{m}_1(X_i) \right] \\
&= E[Y_i(1)]
\end{align}
となります．同様に，(4) 式の右辺第2項の $m_0(x)$ を間違った関数 $\tilde{m}_0(x)$ と置き換えても，
E\left[ \frac{(1-D_i) Y_i}{1-p(X_i)} + \left( 1 - \frac{1-D_i}{1-p(X_i)} \right) \tilde{m}_0(X_i) \right] = E[Y_i(0)]
となります．ここから，もし (4) 式の右辺の条件付き期待値関数 $m_d(x)$ に全く別の関数を入れても，傾向スコア $p(x)$ が正しければ (4) 式の右辺が ATE となることが分かります．
次に，傾向スコア $p(x)$ の推定が失敗していてもよいことを示すために，(4) 式の右辺第1項の $p(x)$ を間違った関数 $\tilde{p}_1(x)$ と置き換えた場合を考えます．このとき，$E[m_1(X_i)] = E[E[Y_i(1)|X_i]]=E[Y_i(1)]$ となるので，
\begin{align}
& E\left[ \frac{D_i Y_i}{\tilde{p}(X_i)} + \left( 1 - \frac{D_i}{\tilde{p}(X_i)} \right) m_1(X_i) \right] \\
= \ & E\left[ \frac{D_i \{Y_i - m_1(X_i)\}}{\tilde{p}(X_i)} \right] + E[m_1(X_i)] \\
= \ & E\left[ \frac{D_i \{ E[Y_i|D_i=1,X_i] - m_1(X_i)\}}{\tilde{p}(X_i)} \right] + E[Y_i(1)] \ = \ E[Y_i(1)]
\end{align}
となります．同様に，(4) 式の右辺第2項も，
E\left[ \frac{(1-D_i) Y_i}{1-\tilde{p}(X_i)} + \left( 1 - \frac{1-D_i}{1-\tilde{p}(X_i)} \right) m_0(X_i) \right] = E[Y_i(0)]
となります．以上から，条件付き期待値関数 $m_d(x)$ の場合と同様，傾向スコアを間違った関数 $\tilde{p}(x)$ と置き換えても (4) 式の右辺が ATE となることが分かります．以上から，doubly robust 推定量は，条件付き期待値関数 $m_d(x)$ や傾向スコア $p(x)$ の特定化の失敗による影響を受けにくい推定量であることが分かります．
最後に
前回に引き続き，条件付き独立の仮定の下での平均処置効果の推定方法を紹介しました．2回の記事で regression adjustment 法，inverse probability weighting 法，doubly robust 法という3つの推定方法を紹介しましたが，これ以外にもマッチング法などの様々な推定方法が提案されています．また，近年では，条件付き期待値関数や傾向スコアの推定に機械学習の手法を用いる doubly robust 推定量も提案されています．
株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospareまでお問い合わせください．
