Double/Debiased Machine Learning による因果効果推定（２）

はじめに

東北大学/株式会社Nospareの石原です．前回の記事に引き続き，機械学習を用いた平均処置効果 (average treatment effect; ATE) の推定方法を紹介します．前回の記事では，条件付き期待値関数と傾向スコアを機械学習手法で推定することによってどのような問題が起きるかを説明しました．本記事では，Chernozhukov et al. (2017) がこの問題点をどのように解決しているかを説明します．

既存の漸近理論

前回の記事と同様に，次のパラメータ $\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$ は未知のベクトル値関数です．このとき，$\eta_0$ の推定量を $\hat{\eta}$ として，

\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))

という推定量を考えます．前回の記事で説明した通り，

\begin{align}
\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} & (1)\\
v_n(\hat{\eta}) - v_n(\eta_0) = o_p(1) \hspace{2in} & (2)
\end{align}

という2つの条件が成り立てば，$\sqrt{n} (\hat{\theta} - \theta_0)$ は漸近的に正規分布に従います．ここで，$\hat{\Psi}_n(\eta) = \frac{1}{n} \sum_{i=1}^n \psi(W_i;\eta)$，$\Psi(\eta) = E[\psi(W_i;\eta)]$，$v_n(\eta) = \sqrt{n} \left\{ \hat{\Psi}_n(\eta) - \Psi(\eta) \right\}$であり，$\phi$ は $E[\phi(W_i)] = 0$ を満たす適当な関数です．

Double/Debiased Machine Learning

前回の記事で紹介した通り，$\hat{\eta}$ の推定に機械学習の手法を用いる場合，(1)，(2) の条件は成り立たない可能性があります．一般的に機械学習手法を用いた推定量は非常に複雑なので，(1) の後者の条件を満たすかどうかを簡単に確かめることはできません．また，例えば $X_i$ が高次元の場合，$\hat{\eta}$ が入る関数のクラスの複雑度が高くなり，(2) の条件が成り立つとは限りません．そこで，Chernozhukov et al. (2017) は，Neyman 直交性とクロスフィッティング (cross-fitting) を用いてこの問題に対処しています．このような手法は double/debiased machine learning とよばれています．

今回の設定では，$\psi$ が次の条件を満たしているとき，Neyman 直交性が成り立つといいます：

\frac{\partial}{\partial r} \Psi(\eta_0 + r \tilde{\eta}) \big|_{r=0} \ = \ 0, \ \ \ \forall \tilde{\eta} \hspace{1in} (3)

この条件は $\eta$ が真の関数 $\eta_0$ から少しずれていても $\Psi(\eta)$ の値は $\Psi(\eta_0)$ と大きく変わらないということを意味しています．つまり，このような $\psi$ を用いれば，最終的な $\theta_0$ の推定量は $\eta_0$ の推定の影響を受けにくくなります．後で確認するように，ターゲットパラメータ $\theta_0$ の識別・推定方法が複数存在する場合があります．Chernozhukov et al. (2017) は，そのような場合に (3) の条件を満たすような $\psi$ を用いた推定方法を選択することを提案しています．実際に，doubly robust 推定量は (3) を満たしており，他の推定量は (3) を満たしていないことを後で確認します．

前回の記事で解説したように，$\eta_0$ を1次元の関数とすると，$\sqrt{n} \left\{ \Psi(\hat{\eta}) - \Psi(\eta_0) \right\}$ は次のように近似することができます：

\begin{align}
\sqrt{n} \left\{ \Psi(\hat{\eta}) - \Psi(\eta_0) \right\} &\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階微分です．このとき，Neyman 直交性は $E\left[ \psi'(Z_i,\eta_0(X_i)) | X_i \right] = 0$ を意味するので，

\sqrt{n} \left\{ \Psi(\hat{\eta}) - \Psi(\eta_0) \right\} \ \approx \ \frac{\sqrt{n}}{2} \int \Big\{\psi''(z,\eta_0(x)) \cdot (\hat{\eta}(x)-\eta_0(x))^2 \Big\} dF_{Z,X}(z,x)

となります．したがって，$\hat{\eta}$ の収束レートが $o_p(n^{-1/4})$ なら，Neyman 直交性から (1) の前者の条件を導けることが分かります．

Chernozhukov et al. (2017) は，(3) の条件を満たすような $\psi$ を用いて，クロスフィッティングを用いた以下のようなアルゴリズムで $\theta_0$ を推定することを提案しています．

1. $n$ 個のデータ $\{1, \ldots , n\}$ をランダムに $K$ 個のグループ $I_1, \ldots, I_K$ に分割する．ここで，簡単のため，各グループのサンプルサイズは $n / K$ であるとし，$I_k^c = \{1, \ldots , n\} \setminus I_k$ と定義する．
2. 各 $k = 1, \ldots , K$ に対して，$I_k^c$ のデータのみを用いて機械学習の手法で $\eta_0$ の推定量 $\hat{\eta}_k$ を得る．この $\hat{\eta}_k$ を用いて

   \check{\theta} (I_k, I_k^c) \ = \ \frac{K}{n} \sum_{i \in I_k} \psi(W_i, \hat{\eta}_k)
   

   という推定量を構成する．
3. $K$ 個の推定量を平均して最終的な推定量を得る：

   \tilde{\theta} \ = \ \frac{1}{K} \sum_{k=1}^K \check{\theta} (I_k, I_k^c)
   

このアルゴリズムでは，2つ目のステップで $\check{\theta} (I_k, I_k^c)$ を構成するときに，$W_i$ とは異なる観測値を用いて $\eta_0$ の推定量 $\hat{\eta}_k$ を構成しています．$\eta_0$ の推定に用いるサンプルと $\psi$ の値を評価するサンプルを分けることで，オーバーフィッティングによるバイアスを取り除いています．このような標本分割の方法をクロスフィッティングとよびます．

このような標本分割によって (2) がどのように満たされるか説明するために，(2) の $\hat{\eta}$ はデータ $\mathcal{D} = \{W_1, \ldots, W_n\}$ とは独立のデータ $\tilde{\mathcal{D}} = \{\tilde{W}_1, \ldots, \tilde{W}_m\}$ を用いて構成されている状況を考えます．このとき，

\begin{align}
Var \left( v_n(\hat{\eta}) - v_n(\eta_0) \Big| \tilde{\mathcal{D}} \right) &= Var \left( \frac{1}{\sqrt{n}} \sum_{i=1}^n \{ \psi(W_i;\hat{\eta}) - \psi(W_i;\eta_0) \} \Big| \tilde{\mathcal{D}}\right) \\
&\leq \int \{ \psi(w;\hat{\eta}) - \psi(w;\eta_0) \}^2 dF_{W}(w) \hspace{0.9in} (4)
\end{align}

となります．例えば，$\{ \psi(w;\hat{\eta}) - \psi(w;\eta_0) \}^2 \leq C \| \hat{\eta}(x) - \eta_0(x) \|^2$ が成り立っているとすると，

\int \| \hat{\eta}(x) - \eta_0(x) \|^2 dF_{X}(x) \ = \ o_p(1)

なら，(4) の右辺は $0$ に確率収束することが分かります．したがって，$\hat{\eta}$ が $\eta_0$ の一致推定量であれば，チェビシェフの不等式から $v_n(\hat{\eta}) - v_n(\eta_0) = o_p(1)$ となることが分かります．クロスフィッティングの場合はもう少し複雑になりますが，$\eta_0$ の推定に用いるサンプルと $\psi$ の値を評価するサンプルが異なっているという事実を用いることで (2) を示すことができます．

以上から，推定量 $\tilde{\theta}$ は適当な仮定の下で漸近正規性をもつことが分かります．また，その漸近分散は

\hat{\sigma}^2 \ = \ \frac{1}{n} \sum_{i=1}^n \left( \hat{\psi}_i - \tilde{\theta} \right)^2

で推定することができます．ただし，$k(i) = \{ k \in \{1, \ldots, K \} : i \in I_k \}$ であり，$\hat{\psi}_i = \psi(W_i, \hat{\eta}_{k(i)})$ です．したがって，Neyman 直交性とクロスフィッティングを用いることで，簡単に95%信頼区間 $[\tilde{\theta} \pm 1.96 \hat{\sigma} / \sqrt{n}]$ を構成することができます．

ATE の推定

以降では，ATE の推定問題での Neyman 直交性と $\hat{\eta}$ の満たすべき条件について考えます．前回の記事で紹介したように，次の $\psi$ を考えます：

\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) \\
&= m_1(X_i) - m_0(X_i) + \frac{D_i(Y_i-m_1(X_i))}{p(X_i)} - \frac{(1-D_i)(Y_i - m_0(X_i))}{1-p(X_i)}
\end{align}

ここで，$Z_i = (Y_i,D_i)$，$\eta_0 = (m_0,m_1,p)$ であり，$m_0,m_1,p$ は $m_d(x) = E[Y_i|D_i=d,X_i=x]$，$p(x) = P(D_i=1|X_i=x)$ を表しています．

まず，この $\psi$ が Neyman 直交性 (3) を満たしていることを確認します．簡単のため，$m_0, m_1$ は既知であり，傾向スコア $p$ のみが未知であるとします．つまり，$\eta_0 = p$ という状況を考えます．このとき，

\begin{align}
\Psi(p+r \tilde{p}) &= E[m_1(X_i)-m_0(X_i)] + E\left[ E\left[ \frac{D_i(Y_i-m_1(X_i))}{p(X_i) + r \tilde{p}(X_i)} \Big| X_i \right] \right] \\
& \hspace{2in} - E\left[ E\left[ \frac{(1-D_i)(Y_i-m_0(X_i))}{1 - p(X_i) - r \tilde{p}(X_i)} \Big| X_i \right] \right] \\
&= E[m_1(X_i)-m_0(X_i)] + E\left[ E\left[ \frac{p(X_i)(m_1(X_i)-m_1(X_i))}{p(X_i) + r \tilde{p}(X_i)} \Big| X_i \right] \right] \\
& \hspace{2in} - E\left[ E\left[ \frac{(1-p(X_i))(m_0(X_i)-m_0(X_i))}{1 - p(X_i) - r \tilde{p}(X_i)} \Big| X_i \right] \right] \\
&= E[m_1(X_i)-m_0(X_i)]
\end{align}

となるので，(3) 式が成り立つこと分かります．これは，以前の記事で紹介した doubly robust 推定量の二重頑健性という性質と関連しています．今回は簡単のために $m_0, m_1$ を既知としましたが，$m_0, m_1$ が未知であるとしても同様に Neyman 直交性が成り立ちます．

以前の記事で紹介したように，doubly robust 推定量以外にも ATE の推定方法は存在します．例えば，ATE の推定方法として inverse probability weighting (IPW) 法が有名ですが，IPW 推定量は Neyman 直交性を満たしません．実際に，

\psi(W_i;\eta_0) \ = \ \frac{D_iY_i}{p(X_i)} + \frac{(1-D_i)Y_i}{1-p(X_i)}

とすると，

\begin{align}
\Psi(p+r \tilde{p}) &= E\left[ E\left[ \frac{D_iY_i}{p(X_i) + r \tilde{p}(X_i)} \Big| X_i \right] \right] - E\left[ E\left[ \frac{(1-D_i)Y_i}{1 - p(X_i) - r \tilde{p}(X_i)} \Big| X_i \right] \right] \\
&= E\left[ \frac{p(X_i)m_1(X_i)}{p(X_i) + r \tilde{p}(X_i)}  \right] - E\left[ \frac{(1-p(X_i))m_0(X_i)}{1- p(X_i) - r \tilde{p}(X_i)}  \right]
\end{align}

となります．よって，微分と期待値の順序交換ができるとすると，

\begin{align}
\frac{\partial}{\partial r} \Psi(p+r \tilde{p}) \Big|_{r=0} &= - E\left[ \frac{\tilde{p}(X_i)p(X_i)m_1(X_i)}{p(X_i)^2} \right] - E\left[ \frac{\tilde{p}(X_i)(1-p(X_i))m_0(X_i)}{(1- p(X_i))^2}  \right] \\
&= - E\left[ \frac{\tilde{p}(X_i)m_1(X_i)}{p(X_i)} \right] - E\left[ \frac{\tilde{p}(X_i)m_0(X_i)}{(1- p(X_i))^2}  \right] 
\end{align}

となるので，一般的には Neyman 直交性が成り立ちません．つまり，傾向スコアを機械学習の手法で推定した場合，IPW 推定量では上手く統計的推測を行えない可能性があります．したがって，条件付き期待値関数と傾向スコアの推定に機械学習の手法を用いる場合は，IPW 推定量ではなく doubly robust 推定量で ATE を推定する必要があります．このように，Neyman 直交性は $\eta_0$ の誤特定に対してより頑健な推定方法を用いる方がよいということを意味しています．

上で確認したように，Neyman 直交性が満たされている場合，$\hat{\eta}$ の収束レートが $o_p(n^{-1/4})$ なら (1) の前者の条件が成り立ちます．つまり，$m_0,m_1,p$ の推定量が全て $o_p(n^{-1/4})$ であれば，クロスフィッティングを用いた doubly robust 推定量は漸近正規性をもちます．しかし，実際には $m_0,m_1,p$ の推定量の全てが $o_p(n^{-1/4})$ でなくても，条件付き期待値関数 $m_d$ と傾向スコア $p$ の収束レートの積が $o_p(n^{-1/2})$ であればよいということが論文で指摘されています．つまり，もし傾向スコア $p$ の収束が十分に速ければ，条件付き期待値関数 $m_d$ の収束が $o_p(n^{-1/4})$ より遅くても条件が満たされます．特に，傾向スコアが既知である場合や傾向スコア $p(x)$ が $x$ に依存せず一定である（$n^{-1/2}$ のレートで推定できる）場合は，条件付き期待値関数 $m_d$ が一致性をもてば収束レートがどれだけ遅くても，クロスフィッティングを用いた doubly robust 推定量は漸近正規性をもちます．

最後に

今回の記事では，Chernozhukov et al. (2017) で提案された機械学習を用いた平均処置効果の推定方法を紹介しました．今回は平均処置効果の推定問題についてのみ解説しましたが，Chernozhukov et al. (2018) では機械学習を用いた一般的なセミパラメトリックモデルの推定方法を提案しており，その推定方法を double/debiased machine learning と名付けています．この手法については末石 (2024) でより詳しく解説されています．また，最近では，因果推論に機械学習の手法を取り入れた研究も多くなっているので，興味のある方は是非勉強してみてください．

株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospareまでお問い合わせください．

参考文献

1. 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.
2. Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., & Robins, J. (2018). Double/debiased machine learning for treatment and structural parameters. Econometrics Journal, 21, C1–C68.
3. 末石直也 (2024). 『データ駆動型回帰分析』日本評論社．