はじめに
東北大学/株式会社Nospareの石原です．本記事では，機械学習を用いた平均処置効果 (average treatment effect; ATE) の推定方法を紹介します．以前の記事で紹介したように，条件付き独立の仮定の下では条件付き期待値関数と傾向スコアを推定することができれば ATE を推定することができます．Chernozhukov et al. (2017) は，条件付き期待値関数と傾向スコアを機械学習手法で推定することで ATE を推定するという方法を提案しています．彼らの提案した方法は "Double/Debiased Machine Learning" と呼ばれています．本記事では，通常のノンパラメトリック推定方法を用いた場合の既存の理論結果を紹介し，機械学習手法を用いる場合の理論的な問題点を解説します．そして，次回の記事で，Chernozhukov et al. (2017) がこの問題点をどのように解決しているかを説明します．
ATE の推定
処置を受けたときに実現するであろう潜在的な結果変数を $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
の下で次のように ATE を識別することができます：
\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} (1)
\end{align}
ここで，$m_d(x) = E[Y_i|D_i=d,X_i=x]$，$p(x) = P(D_i=1|X_i=x)$ です．(1) の条件付き期待値関数 $m_d(\cdot)$ と傾向スコア $p(\cdot)$ を推定量 $\hat{m}_d(\cdot), \ \hat{p}(\cdot)$ に置き換えることで
\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\} \hspace{1in} (2)
\end{align}
という推定量を得ることができます．この推定量は doubly robust 推定量と呼ばれています．
観測される共変量の次元があまり高くない場合は $m_d(\cdot)$ と $p(\cdot)$ を通常のノンパラメトリック手法で推定することが可能ですが，次元が高い場合は通常のノンパラメトリック手法の精度はあまり良くないことが知られています．そこで，Chernozhukov et al. (2017) では，機械学習の手法を用いて $m_d(\cdot)$ と $p(\cdot)$ を推定する方法を提案しています．しかし，機械学習の手法で得られた $\hat{m}_d(\cdot)$ と $\hat{p}(\cdot)$ を直接代入した (2) のような推定量では，仮説検定や信頼区間などの統計的推測を正確に行うことができません．以降では，通常のノンパラメトリック手法を用いた場合の既存の理論結果を紹介し，機械学習の手法を使用した場合にどのような仮定が満たされないのか解説します．
既存の漸近理論
次のパラメータ $\theta_0 \in \mathbb{R}$ を推定したいという問題を考えます：
\theta_0 = E[\psi(W_i;\eta_0)] = E[\psi(Z_i,\eta_0(X_i))]
ここで，$W_i = (Z_i,X_i)$ であり，$\eta_0$ は未知のベクトル値関数です．上の設定では，$Z_i = (Y_i,D_i)$，$\eta_0 = (m_0,m_1,p)$，
\begin{align}
\psi(W_i;\eta_0) &= \frac{D_i Y_i}{p(X_i)} + \left( 1 - \frac{D_i}{p(X_i)} \right) m_1(X_i) \\
& \hspace{0.4in} - \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)
\end{align}
と考えることができます．このとき，$\eta_0$ の推定量 $\hat{\eta}$ があれば，パラメータ $\theta_0$ の自然な推定量は
\hat{\theta} = \frac{1}{n} \sum_{i=1}^n \psi(W_i;\hat{\eta}) = \frac{1}{n} \sum_{i=1}^n \psi(Z_i,\hat{\eta}(X_i))
となります．実際に，上で紹介した doubly robust 推定量 (2) は，この推定量と同じ構造をもっています．
統計的推測を行うためには，$\hat{\theta}$ の漸近分布を求める必要があります．以下では，$\hat{\theta}$ の漸近分布を求めるにはどのような仮定が必要になるかを解説します．関数 $\eta$ に対して $\hat{\Psi}_n(\eta) = \frac{1}{n} \sum_{i=1}^n \psi(W_i;\eta)$，$\Psi(\eta) = E[\psi(W_i;\eta)]$ と定義すると，$\sqrt{n} (\hat{\theta} - \theta_0)$ は次のように書き直すことができます：
\begin{align}
\sqrt{n} (\hat{\theta} - \theta_0) &= \sqrt{n} \left\{ \hat{\Psi}_n(\hat{\eta}) - \Psi(\eta_0) \right\} \\
&= \sqrt{n} \left\{ \hat{\Psi}_n(\eta_0) - \Psi(\eta_0) \right\} + \sqrt{n} \left\{ \hat{\Psi}_n(\hat{\eta}) - \Psi(\hat{\eta}) \right\} \\
& \hspace{0.5in} - \sqrt{n} \left\{ \hat{\Psi}_n(\eta_0) - \Psi(\eta_0) \right\} + \sqrt{n} \left\{ \Psi(\hat{\eta}) - \Psi(\eta_0) \right\} \\
&= \sqrt{n} \left\{ \hat{\Psi}_n(\eta_0) - \Psi(\eta_0) \right\} + \left\{ v_n(\hat{\eta}) - v_n(\eta_0) \right\} + \sqrt{n} \left\{ \Psi(\hat{\eta}) - \Psi(\eta_0) \right\}
\end{align}
ここで，$v_n(\eta) = \sqrt{n} \left\{ \hat{\Psi}_n(\eta) - \Psi(\eta) \right\}$ です．以上から，次の2つの条件があれば $\sqrt{n} (\hat{\theta} - \theta_0)$ が漸近的に正規分布に従うことが分かります：
\begin{align}
v_n(\hat{\eta}) - v_n(\eta_0) = o_p(1) \hspace{2in} & (3)\\
\sqrt{n} \left\{ \Psi(\hat{\eta}) - \Psi(\eta_0) \right\} = o_p(1) \ \ \text{or} \ \ \frac{1}{\sqrt{n}} \sum_{i=1}^n \phi(W_i) + o_p(1) \hspace{0.3in} & (4)
\end{align}
ここで，$\phi$ は $E[\phi(W_i)] = 0$ を満たす適当な関数です．実際に，(3)，(4) の仮定が成り立つなら，$\tilde{\psi}(W_i) = \psi(W_i;\eta_0) - E[\psi(W_i;\eta_0)]$ と表記すると，
\begin{align}
\sqrt{n} (\hat{\theta} - \theta_0) &= \frac{1}{\sqrt{n}} \sum_{i=1}^n \tilde{\psi}(W_i) + o_p(1)
\end{align}
または
\begin{align}
\sqrt{n} (\hat{\theta} - \theta_0) &= \frac{1}{\sqrt{n}} \sum_{i=1}^n \left\{ \tilde{\psi}(W_i) + \phi(W_i) \right\} + o_p(1)
\end{align}
が成り立ちます．以上から，(3) と (4) の仮定の下で，$\hat{\theta}$ が漸近正規性をもつことが分かります．特に，(4) の仮定が $\sqrt{n} \left\{ \Psi(\hat{\eta}) - \Psi(\eta_0) \right\} = o_p(1)$ である場合には，漸近的に $\eta_0$ が既知の場合と同じ漸近分散を得ることができます．
以下では，どのような条件の下で (3) と (4) の仮定が成り立つのかについて議論します．$\hat{\eta}$ が適当な関数のクラス $\mathcal{H}$ に入っており，その関数のクラス $\mathcal{H}$ の複雑度があまり高くないとき，(3) の仮定は成り立つことが知られています．例えば，$\mathcal{H}$ が十分に滑らかな関数のクラスであり，$P(\hat{\eta} \in \mathcal{H}) \to 1$ が成り立つなら，(3) の仮定は成り立ちます．次に，(4) の仮定がどのような条件の下で成り立つか議論します．簡単のため，$\eta_0$ を1次元の関数とし，$\psi(z, \hat{\eta}(x))$ を $\eta_0(x)$ の周りでテイラー展開すると，
\begin{align}
\sqrt{n} \left\{ \Psi(\hat{\eta}) - \Psi(\eta_0) \right\} &= \sqrt{n} E_{Z,X} \left[ \psi(Z, \hat{\eta}(X)) - \psi(Z, \eta_0(X)) \right] \\
&= \sqrt{n} \int \psi(z, \hat{\eta}(x)) - \psi(z, \eta_0(x)) dF_{Z,X}(z,x) \\
&\approx \sqrt{n} \int \Big\{ \psi'(z,\eta_0(x)) \cdot (\hat{\eta}(x)-\eta_0(x)) \\
& \hspace{0.8in} + \frac{1}{2} \psi''(z,\eta_0(x)) \cdot (\hat{\eta}(x)-\eta_0(x))^2 \Big\} dF_{Z,X}(z,x)
\end{align}
と近似できます．$\psi'(z,h)$ と $\psi''(z,h)$ はそれぞれ $\psi(z,h)$ の $h$ についての1階微分と2階微分です．このとき，$\hat{\eta}$ の収束レートが $o_p(n^{-1/4})$ なら，2階微分の項は $o_p(1)$ となり無視することができます．一方で，1階微分の項は
\begin{align}
& \sqrt{n} \int \Big\{ \psi'(z,\eta_0(x)) \cdot (\hat{\eta}(x)-\eta_0(x)) \Big\} dF_{Z,X}(z,x) \\
=& \sqrt{n} \int E[\psi'(Z,\eta_0(x))|X=x] \cdot f_X(x) \cdot (\hat{\eta}(x)-\eta_0(x)) dx \\
=& \sqrt{n} \int \delta(x) f_X(x) (\hat{\eta}(x)-\eta_0(x)) dx
\end{align}
と書くことができます．ここで，$\delta(x) = E[\psi'(Z,\eta_0(x))|X=x]$ です．したがって，$\delta(x) = 0$ なら，1階微分の項は 0 となるので $\sqrt{n} \left\{ \Psi(\hat{\eta}) - \Psi(\eta_0) \right\} = o_p(1)$ が成り立ちます．一方で，$\delta(x) \neq 0$ なら，
\begin{align}
\sqrt{n} \int \delta(x) f_X(x) (\hat{\eta}(x)-\eta_0(x)) dx = \frac{1}{\sqrt{n}} \sum_{i=1}^n \phi(W_i) + o_p(1) \hspace{0.5in} (5)
\end{align}
となる必要があります．(5) が成り立つかを調べることは大変ですが，通常のノンパラメトリック推定量のように推定量が陽に求められる場合は，(5) の条件を直接確かめることができます．
機械学習手法を用いることの問題点
上で説明したように，$\hat{\theta}$ の漸近正規性を示すには，仮定 (3)，(4) を示す必要があります．しかし，$\hat{\eta}$ の推定に機械学習の手法を用いる場合，(3)，(4) の仮定は成り立たない可能性があります．例えば，$X_i$ が高次元の場合，$\hat{\eta}$ が入る関数のクラスの複雑度が高くなり，(3) の仮定は成り立たないかもしれません．また，一般的に機械学習手法を用いた推定量は非常に複雑なので，直接 (5) の条件が成り立つかを確認することは難しいです．したがって，$\delta(x) \neq 0$ の場合には，機械学習手法を用いると (4) の仮定が成り立っていない可能性があります．Chernozhukov et al. (2017) では，クロスフィッティング (cross-fitting) と Neyman 直交性を用いて，この問題を解決しています．詳しくは，次回の記事で紹介します．
最後に
今回の記事では，通常のノンパラメトリック推定方法を用いた場合の既存の理論結果を紹介し，機械学習手法を用いる場合の理論的な問題点を解説しました．今回紹介した既存の漸近理論は Andrews (1994) の結果を参照しています．Chen et al. (2003) では，Andrews (1994) の結果を微分できない関数に拡張しており，セミパラメトリック推定量のより一般的な理論結果を求めています．また，今回紹介した内容は西山・人見 (2023) や末石 (2024) などの日本語のテキストでもより詳しく解説されているので，興味のある人はテキストを読んでみてください．
株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospareまでお問い合わせください．
参考文献
- Andrews, D. W. (1994). Asymptotics for semiparametric econometric models via stochastic equicontinuity. Econometrica, 43-72.
- Chen, X., Linton, O., & Van Keilegom, I. (2003). Estimation of semiparametric models when the criterion function is not smooth. Econometrica, 71(5), 1591-1608.
- Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., & Newey, W. (2017). Double/debiased/neyman machine learning of treatment effects. American Economic Review, 107(5), 261-265.
- 末石直也 (2024). 『データ駆動型回帰分析』日本評論社．
- 西山慶彦・人見光太郎 (2023). 『ノン・セミパラメトリック統計解析』 共立出版．