【Mostly Harmless Ch.3.3.2】傾向スコアを使った共変量のコントロール

はじめに

この記事では傾向スコアマッチング（Propensity Score matching）ついて扱います。内容はJoshua D. Angrist & Jorn-steffen Pischke (2008)『Mostly Harmless Econometrics』Ch.3.3.2"Control for Covariates Using the Propensity Score"を参考にしています。

はじめに
傾向スコアとは何か
- 傾向スコア定理（The Propensity Score Theorem）
傾向スコアとIPW推定量
IPW推定量の性質
傾向スコア推定量と回帰係数
終わりに・感想
参考資料

傾向スコアとは何か

回帰の理論の中で最も重要なものの1つが欠落変数バイアスです。回帰モデルに含まれる説明変数 $x_{1i}$ の係数 $β_{1i}$ は、 $x_{1i}$ と回帰モデルから除外された（欠落した）説明変数 $x_{2i}$ との相関がない限りは、何の影響も受けません。Rosenbaum and Rubin(1983)*1による傾向スコア定理は、因果関係を示す説明変数が処置変数ダミー（treatment dummy ）である回帰分析の考え方を、マッチングへと拡張したものです。

傾向スコア定理は、潜在結果が多次元共変量ベクトル $X_i$ を条件付けた上で介入変数と独立であるとき、潜在結果が多次元共変量ベクトルのスカラー関数を条件付けた上で介入変数と独立であることを意味します。共変量ベクトルによるコントロールの機能が、共変量ベクトルの情報を含むスカラーによっても満たされます。このスカラーが傾向スコア（Propensity Score）です。

傾向スコアは共変量を $X_i$ ,介入変数を $D_i$ とするとき、 $p(X_i)=E[D_i|X_i]=P[D_i=1|X_i]$ と定義され、介入ダミーの条件付き期待値＝介入 $X_i$ が想起する条件付き確率を意味します。

傾向スコア定理（The Propensity Score Theorem）

条件付き独立の仮定（CIA）が潜在結果 $Y_{0i},Y_{1i}$ に対して成立するとき、
$Y_{0i},Y_{1i}$ は傾向スコア $p(X_i)$ で条件付けた介入変数 $D_i$ と独立となる。

$\{Y_{0i},Y_{1i}\} \bot D_i|\textit{p}(X_i)\tag{1}$

傾向スコア定理は、介入が想起する可能性 $P[D_i=1|Y_{ji},p(X_i)]$ が潜在結果 $Y_{ji},j=\{0,1\}$ に依存しないことを示せば証明できます。

${ \begin{eqnarray} P[D_i=1|Y_{ji},\textit{p}(X_i)] &=& E[D_i|Y_{ji},\textit{p}(X_i)] \\ &=& E\{E[D_i|Y_{ji},\textit{p}(X_i),X_i]|Y_{ji},\textit{p}(X_i)\} \\ &=& E\{E[D_i|Y_{ji},X_i]|Y_{ji},\textit{p}(X_i)\} \\ &=& E\{E[D_i|X_i]|Y_{ji},\textit{p}(X_i) \} \\ &=& E\{\textit{p}(X_i)|Y_{ji},\textit{p}(X_i)]\} \\ &=& p(X_i) \tag{2} \end{eqnarray} }$

5番目の等号では $E[D_i|X_i]=p(X_i)$ の関係を用いています。したがって傾向スコア $p(X_i)$ を所与として、潜在結果 $Y_{0i},Y_{1i}$ の期待値の差をとると、平均処置効果（ATT）を求めることができるのです。

傾向スコアとIPW推定量

欠落変数バイアスと同様、傾向スコア定理の意味するのは、介入変数の想起する確率（probability of treatment）に影響を与える共変量のみコントロールすれば介入変数と潜在結果の条件付き独立は成立する、ということです。加えて傾向スコア定理は、本当にコントロールが必要な唯一の共変量は介入の確率そのものだ、ということも意味します。
傾向スコアを推定に用いる場合、二段階推定が必要です。

ロジットやプロビットモデルによる傾向スコア $p(X_i)$ の推定。
マッチングまたは重み付けによる1.の $p(X_i)$ の推定値の調整。*2

マッチングには実際に複数の手法があります。共変量マッチングは近しい共変量の値を持つ処置群と対照群のサンプルをペアにするという方法であり、傾向スコアマッチングの場合は $p(X_i)$ の値の近さをもってペアを作ります。

傾向スコア定理と条件付き独立の仮定（CIA）により、処置群における平均処置効果（ATT）は、

${ \begin{eqnarray} E[Y_{1i}-Y_{0i}|D_i=1]&=& E\{E[Y_i|\textit{p}(X_i),D_i=1]-E[Y_i|\textit{p}(X_i),D_i=0]|D_i=1\} \tag{3} \end{eqnarray} }$

と表せます。

ATTの推定は傾向スコアの大小によってサンプルをいくつかの層に分け（stratifying）、各層におけるATTを(3)の内側の期待値に代入するという方法や、傾向スコアの同じまたは類似する処置群と介入群のサンプル同士をペアにして、その差を値をATTと見なす、という方法があります。*3もう1つは、 $E[Y_{i}|p(X_i),D_i]$ の推定値をこれらの条件付き期待値の代わりに使用し、(3)の外側の $E[ ・ ]$ を $\sum$ に置き換えるものです。*4

我々はCIAの示唆する $E \left[ \frac{Y_i D_i}{p(X_i)} \right]=E[Y_{1i}],E \left[ \frac{Y_i (1-D_i)}{(1-p(X_i))} \right]=E[Y_{0i}]$ という関係を用いることで、面倒なマッチングの作業をスキップすることができます。 $p(X_i)$ の推定方法を所与として、平均処置効果（ATE）の推定値を、標本対応（sample analoge）つまり観察データから構築することができるのです。平均処置効果（ATE）は、

${ \begin{eqnarray} E[Y_{1i}-Y_{0i}]&=& E\left[\frac{Y_{i}D_{i}}{\textit{p}(X_i)} - \frac{Y_i(1-D_i)}{1-\textit{p}(X_i)}\right] \\ &=& E\left[\frac{(D_i-\textit{p}(X_i))Y_{i}}{(1-\textit{p}(X_i))}\right] \tag{4} \end{eqnarray} }$

と表すことができます。(4)はNewey(1990) *5とRobins, Mark and Newey(1992) *6で示されたものです。同様に処置群への平均処置効果（ATT）を計算すると、

${ \begin{eqnarray} E[Y_{1i}-Y_{0i}|D_i=1]&=& E\left[\frac{(D_i-\textit{p}(X_i))Y_{i}}{(1-\textit{p}(X_i))(P(D_i=1))}\right] \tag{5} \end{eqnarray} }$

(4),(5)では、処置群の潜在結果には傾向スコアの逆数 $\frac{1}{p(X_i)}$ 、介入群の潜在結果には $\frac{1}{1-p(X_i)}$ の重みが付けられています。これらの潜在結果は実際に観察することができないため、観察できる標本を使って推定することになります。傾向スコアの逆数で重みを付けた潜在結果の期待値(6)(7)をIPW（Inverse Probability Weighting）推定量と呼びます。

${ \begin{eqnarray} E[\hat{Y_{1i}}]&=& \frac{\sum_{i=1}^{N} \frac{Y_iD_i}{p(X_i)}}{\sum_{i=1}^{N} \frac{D_i}{p(X_i)}} \tag{6} \\ E[\hat{Y_{0i}}]&=& \frac{\sum_{i=1}^{N} \frac{Y_i (1-D_i)}{1-p(X_i)}}{\sum_{i=1}^{N} \frac{(1-D_i)}{1-p(X_i)}} \tag{7} \end{eqnarray} }$

IPW推定量の性質

$E \left[ \frac{Y_i D_i}{p(X_i)} \right]=E[Y_{1i}],E \left[ \frac{Y_i (1-D_i)}{(1-p(X_i))} \right]=E[Y_{0i}]$ はIPW推定量の不偏性（unbiasedness）という性質を用いています。これを確認するため、(6)の期待値をとってみましょう。

${ \begin{eqnarray} E \left[\frac{\sum_{i=1}^{N} \frac{Y_iD_i}{p(X_i)}}{\sum_{i=1}^{N} \frac{D_i}{p(X_i)}} \right]&=& E \left[\frac{\frac{1}{N}\sum_{i=1}^{N} \frac{Y_iD_i}{p(X_i)}}{\frac{1}{N}\sum_{i=1}^{N} \frac{D_i}{p(X_i)}} \right]\\ &=& E \left[ \frac{1}{N}\sum_{i=1}^{N} \frac{Y_iD_i}{p(X_i)} \right] \\ &=& E\left[E \left[ \frac{Y_i D_i}{p(X_i)} |X_i \right]\right] \tag{8}\\ &=& E\left[ \frac{Y_i D_i}{p(X_i)}\right] \tag{9} \end{eqnarray} }$

(8)の内側の期待値の中身を計算すると、

${ \begin{eqnarray} E \left[ \frac{Y_i D_i}{p(X_i)} |X_i \right] &=& E[Y_i|X_i]E\left[\frac{D_i}{p(X_i)}|X_i\right] \\　 &=& E[Y_i|D_i=1,X_i] E\left[\frac{D_i}{p(X_i)}|X_i\right] \\　 &=& \frac{E[Y_i|D_i=1,X_i]p(X_i)}{p(X_i)} \\ &=& E[Y_i|D_i=1,X_i] \\ &=& E[Y_{1i}|X_i] \tag{10} \end{eqnarray} }$

1つ目の等号は $Y_i$ と $D_i$ は $X_i$ の条件付き独立であるためです。(10)を(8)に代入すると、繰り返し期待値の法則により

${ \begin{eqnarray} E \left[ \frac{Y_i D_i}{p(X_i)} \right] &=& E\left\{ E[Y_{1i}|X_i]\right\} \\ &=& E[Y_{1i}] \tag{11} \end{eqnarray} }$

(4),(5)で用いた関係を導出できました。(7)の場合も同様です。

IPW推定量のもう一つの性質は一致性（consistency）です。再び $Y_{1i}$ の場合を示します。

${ \begin{eqnarray} \frac{\sum_{i=1}^{N} \frac{Y_iD_i}{p(X_i)}}{\sum_{i=1}^{N} \frac{D_i}{p(X_i)}} = \frac{\frac{1}{N}\sum_{i=1}^{N} \frac{Y_iD_i}{p(X_i)}}{\frac{1}{N}\sum_{i=1}^{N} \frac{D_i}{p(X_i)}} \stackrel{p}{\longrightarrow }　\frac{E\left[\frac{Y_i D_i}{p(X_i)}\right]}{E\left[\frac{D_i}{p(X_i)}\right]} =E[Y_{1i}] \tag{12} \end{eqnarray} }$

サンプルサイズが十分に大きなとき、IPW推定量は潜在結果の期待値に確率収束します。(7)の $Y_{0i}$ の場合も同様に示すことができます。

傾向スコア推定量と回帰係数

傾向スコアの逆数で重み付けするという考えは、Horvitz and Thompson(1952)*7によるものです。

Horvitz-Thompsonによるマッチング方法は、面倒なマッチングの手順が必要ありません。また傾向スコアマッチングと回帰の間のつながりも示唆しています。回帰の推定対象は以下のように表せます。

${ \begin{eqnarray} δ_R &=& \frac{E[(D_i-\textit{p}(X_i))Y_{i}]}{E[\textit{p}(X_i)(1-\textit{p}(X_i))]} \tag{13} \end{eqnarray} }$

Horvitz-Thompsonマッチングの推定対象と回帰の推定対象はどちらも、Hirano,Imbens and Ridder（2003）*8によって検討された加重平均推定対象に分類されます。

${ \begin{eqnarray} E\left\{ g(X_i)\left[\frac{Y_{i}D_{i}}{\textit{p}(X_i)}-\frac{Y_i(1-D_i)}{(1-\textit{p}(X_i))} \right]\right\} \tag{14} \end{eqnarray} }$

$g(X_i)$ は既知の重み付け関数です。(14)は推定対象（estimand）を表しており、実際に推定量を算出する段階では $p(X_i)$ を一致推定量 $\hat{p(X_i)}$ に、期待値 $E[・]$ を合計 $\sum_i$ にそれぞれ置き換えて考えます。平均処置効果（ATE）を推定するには $g(X_i)=1$ 、処置群への平均処置効果（ATT）を推定するには $g(X_i)=\frac{p(X_i)}{P(D_i=1)}$ 、そして回帰の場合には以下のように仮定します。

${ \begin{eqnarray} g(X_i) &=& \frac{\textit{p}(X_i)(1-\textit{p}(X_i))}{E[\textit{p}(X_i)(1-\textit{p}(X_i))]} \tag{15} \end{eqnarray} }$

前回の記事で、回帰係数は説明変数の分散で重み付けられたものだと説明しましたが、まさに(15)の分子と分母には二項変数 $X_i$ の分散が登場しています。傾向スコアを明示的にモデルに組み込まない限りは、回帰とマッチングは別物ではないということが示唆されます。

終わりに・感想

大変読みにくい章です。IPW推定量の部分は合ってるのかどうか自信がありません。読んでいただいてありがとうございます。

参考資料

Mostly Harmless Econometrics: An Empiricist's Companion

作者:Angrist, Joahua D.,Pischke, Jorn-steffen
Princeton University Press

Amazon

日本語で読める傾向スコアマッチングの教科書としてとても有名ですが、以下の文献を参考にさせていただきました。

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

作者:星野崇宏
岩波書店

Amazon

以下の論文も傾向スコアを用いた調整に関するモチベーションが整理されており、非常に勉強になりました。

星野崇宏，繁桝算男（2004）. 傾向スコア解析法による因果効果の推定と調査データの調整について，行動計量学, 31, 43–61.

*1:PAUL R. ROSENBAUM, DONALD B. RUBIN, The central role of the propensity score in observational studies for causal effects, Biometrika, Volume 70, Issue 1, April 1983, Pages 41–55, https://doi.org/10.1093/biomet/70.1.41

*2:Imbens G. Nonparametric Estimation of Average Treatment Effects under Exogeneity: A Review. Review of Economics and Statistics. 2004.

*3:Dehejia and Wahba(1999)((Rajeev H. Dehejia & Sadek Wahba (1999) Causal Effects in Nonexperimental Studies: Reevaluating the Evaluation of Training Programs, Journal of the American Statistical Association, 94:448, 1053-1062, DOI: 10.1080/01621459.1999.10473858

*4:Heckman,Ichimura,Todd(1998)((Heckman, J. J., Ichimura, H., & Todd, P. (1998). Matching As An Econometric Evaluation Estimator. Review of Economic Studies, 65(2), 261-294. https://doi.org/10.1111/1467-937X.00044

*5:Newey, Whitney K, 1990. "Semiparametric Efficiency Bounds," Journal of Applied Econometrics, John Wiley & Sons, Ltd., vol. 5(2), pages 99-135, April-Jun.

*6:Robins, J., Mark, S., & Newey, W. (1992). Estimating Exposure Effects by Modelling the Expectation of Exposure Conditional on Confounders. Biometrics, 48(2), 479-495. doi:10.2307/2532304

*7:D. G. Horvitz & D. J. Thompson (1952) A Generalization of Sampling Without Replacement from a Finite Universe, Journal of the American Statistical Association, 47:260, 663-685, DOI: 10.1080/01621459.1952.10483446

*8:Keisuke Hirano & Guido W. Imbens & Geert Ridder, 2003. "Efficient Estimation of Average Treatment Effects Using the Estimated Propensity Score," Econometrica, Econometric Society, vol. 71(4), pages 1161-1189, July.

Goodな生活

INTP型コンサルタントが好奇心の受け皿を探す