はじめに
Double/Debiased Machine Learning（以下，DML）を解説している書籍や記事を読んでいると，どうしても「機械学習手法を使った統計的因果推論の（比較的）新しい手法」として紹介されているような印象を受ける．
勿論，DMLではランダムフォレストや勾配ブースティング，Lasso，ニューラルネットワークのような機械学習手法を用いるため，機械学習との関係が深いことは言を俟たない．しかしながら，DMLの本質は「機械学習手法そのものを用いること」ではないように思う．寧ろDMLは，関心のあるパラメータと，関心はないが推定には必要な関数を分けて考え，後者の推定誤差が前者の推定に悪影響を及ぼさないようにするための方法として解釈できる．
その意味においては，DMLは回帰分析の古典的な考え方の延長線上にあると考えて良いだろう．特に重回帰分析におけるFrisch-Waugh-Lovell定理（以下，FWL定理）を理解していると，DMLの直観は非常に掴みやすくなる．FWL定理は，誤解を恐れずに言えば，重回帰分析において，関心のある変数とコントロール変数を区別して考え，コントロール変数で説明できる部分を取り除いた残差同士を回帰しても，元の重回帰分析と同じ係数が得られる，という定理である．
DMLの部分線形回帰モデルにおける基本的な考え方も，これに非常に近い．相違点は，コントロール変数の影響を線形関数によって取り除くのではなく，より柔軟な（非線形）関数で取り除く点にある．
この記事では，DMLをいきなり条件付き平均処置効果（Conditional Average Treatment Effect：CATE）や因果機械学習の文脈から説明するのではなく，まずは古典的な重回帰分析，FWL定理，部分線形回帰モデルの延長として説明する．その上で，シミュレーションデータを用いて，Pythonで（EconMLは使わずに）DMLを実装する．
重回帰分析の基本的理解
まず，通常の重回帰分析を考える．各個体$i = 1, 2, \cdots, n$について，アウトカムを$Y_{i}$，関心のある変数を$D_{i}$，$k$次元のコントロール変数ベクトルを$\boldsymbol{X}_{i}$とする．ここでは，$D_{i}$を処置変数のように読んでも良いし，価格，政策変数，説明変数のように読んでも良い．
重回帰モデルは次のように書ける．
$$
Y_{i} = \theta_{0} D_{i} + \boldsymbol{X}_{i}^{\top} \boldsymbol{\beta}_{0} + \varepsilon_{i}
$$
ここで，関心があるのは$D_{i}$の係数$\theta_{0}$である．一方で$\boldsymbol{\beta}_{0}$は，$\boldsymbol{X}_{i}$をコントロールするために必要なパラメータではあるものの，それ自体が主たる関心対象であるとは限らない．
このとき，最小二乗法（Ordinary Least Squares：OLS）では$Y_{i}$を$D_{i}$と$\boldsymbol{X}_{i}$に同時に回帰する．これ自体は非常に馴染み深い操作だが，この操作の背後では「$D_{i}$のうち$\boldsymbol{X}_{i}$で説明できる部分」と「$Y_{i}$のうち$\boldsymbol{X}_{i}$で説明できる部分」を取り除いたうえで，残った変動同士の関係を見ている，と解釈できる．このことを定式化したのが，後述するFWL定理である．
FWL定理の簡単な解説
以下，行列表記でFWL定理を説明する（ベクトル表記だと何かと面倒なので）．アウトカムを並べたベクトルを，
$$
\boldsymbol{Y} =
\begin{pmatrix}
Y_{1} & Y_{2} & \cdots & Y_{n}
\end{pmatrix}^{\top}
$$
とする．同様に，関心のある変数を並べたベクトルを
$$
\boldsymbol{D} =
\begin{pmatrix}
D_{1} & D_{2} & \cdots & D_{n}
\end{pmatrix}^{\top}
$$
とする．また，コントロール変数を並べた行列を
\boldsymbol{X} =
\begin{pmatrix}
\boldsymbol{X}_{1}^{\top} \\
\boldsymbol{X}_{2}^{\top} \\
\vdots \\
\boldsymbol{X}_{n}^{\top} \\
\end{pmatrix}
とする．このとき，重回帰モデルは
$$
\boldsymbol{Y} = \boldsymbol{D} \theta_{0} + \boldsymbol{X} \boldsymbol{\beta_{0}} + \boldsymbol{\varepsilon}
$$
と書ける．
ここで，$\boldsymbol{X}$への射影行列（Projection Matrix）を
$$
\boldsymbol{P_{X}} = \boldsymbol{X} \left( \boldsymbol{X}^{\top} \boldsymbol{X} \right)^{-1} \boldsymbol{X}^{\top}
$$
とする（当然，$\boldsymbol{P_{X}}$は$n$次正方行列である）．また，$\boldsymbol{X}$で説明できる部分を取り除く行列を
$$
\boldsymbol{M_{X}} = \boldsymbol{I} - \boldsymbol{P_{X}}
$$
とする（そして，これもまた射影行列である）．$\boldsymbol{P_{X}}$は，ある変数を$\boldsymbol{X}$で回帰したときの予測値を作る行列であり，$\boldsymbol{M_{X}}$は，ある変数を$\boldsymbol{X}$で回帰したときの残差を作る行列である．
したがって，
$$
\tilde{\boldsymbol{Y}} = \boldsymbol{M_{X} Y}
$$
は，$\boldsymbol{Y}$を$\boldsymbol{X}$に回帰した残差に他ならない．同様に，
$$
\tilde{\boldsymbol{D}} = \boldsymbol{M_{X} D}
$$
は，$\boldsymbol{D}$を$\boldsymbol{X}$に回帰した残差に他ならない．この辺りの議論をもう少し詳しく知りたい方はHansen (2022)を参照されたい．
閑話休題，FWL定理は，$\boldsymbol{Y}$を$\boldsymbol{D}$と$\boldsymbol{X}$に回帰して得られる$\boldsymbol{D}$の係数が，$\tilde{\boldsymbol{Y}}$を$\tilde{\boldsymbol{D}}$
に回帰して得られる係数と一致することを主張する．すなわち，
$$
\hat{\theta}_{OLS} = \left( \tilde{\boldsymbol{D}}^{\top} \tilde{\boldsymbol{D}} \right)^{-1} \tilde{\boldsymbol{D}}^{\top} \tilde{\boldsymbol{Y}}
$$
である．
つまり，重回帰分析における$D_{i}$の係数は，$Y_{i}$と$D_{i}$から$\boldsymbol{X}_{i}$で説明できる部分をそれぞれ取り除き，その残差同士を単回帰した係数として理解できる．
FWL定理の意味
FWL定理の意味を少し丁寧に見ておく．
重回帰分析で$D_{i}$の係数$\theta_{0}$を推定したいとき，問題になるのは$D_{i}$と$\boldsymbol{X}_{i}$が相関していることである．例えば，教育年数$D_{i}$が賃金$Y_{i}$に与える影響を考えるとき，家庭環境や能力のような変数$\boldsymbol{X}_{i}$が教育年数とも賃金とも関係しているかも知れない．このとき，$Y_{i}$を$D_{i}$だけに回帰すると，$D_{i}$の係数には$\boldsymbol{X}_{i}$の影響も混ざってしまう．
重回帰分析では，$\boldsymbol{X}_{i}$を同時に入れることでこの問題に対処する．しかし，FWL定理の見方をすると，これは次の2段階の操作として理解できる．
第1に，$Y_{i}$から$\boldsymbol{X}_{i}$で説明できる部分を取り除く．
$$
\tilde{Y}_{i} = Y_{i} - \hat{\mathbb{E}} [Y_{i} | \boldsymbol{X_{i}}]
$$
第2に，$D_{i}$から$\boldsymbol{X}_{i}$で説明できる部分を取り除く．
$$
\tilde{D}_{i} = D_{i} - \hat{\mathbb{E}} [D_{i} | \boldsymbol{X_{i}}]
$$
そのうえで，$\tilde{Y}_{i}$を$\tilde{D}_{i}$に回帰する．
$$
\tilde{Y}_{i} = \theta \tilde{D}_{i} + \nu_{i}
$$
ここで重要なのは，$Y_{i}$だけでなく，$D_{i}$からも$\boldsymbol{X}_{i}$の影響を取り除いている点である．$Y_{i}$だけを補正しても不十分である．$D_{i}$のうち，$\boldsymbol{X}_i$によって説明される部分を使って$\theta_{0}$を推定してしまうと，結局，$\boldsymbol{X}_{i}$の影響を通じた比較をしていることになるからである．
この「$Y_{i}$と$D_{i}$の両方を残差化する」という考え方が，DMLの部分線形回帰モデルにおいても中心的な役割を果たす．
線形回帰モデルから部分線形回帰モデルへ
ここまでは，$\boldsymbol{X}_{i}$の影響を線形で表していた．
$$
Y_{i} = \theta_{0} D_{i} + \boldsymbol{X}_{i}^{\top} \boldsymbol{\beta}_{0} + \varepsilon_{i}
$$
しかし，実際には$\boldsymbol{X}_{i}$の影響が線形であるとは限らない．例えば，年齢の効果は非線形かも知れない．所得の効果も逓減的かも知れない．変数同士の交互作用もあるかも知れない．高次元のテキスト特徴量や画像特徴量を使う場合には，そもそも線形回帰モデルを仮定すること自体がかなり強い仮定である．
そこで，次の部分線形回帰モデルを考える．
$$
Y_{i} = \theta_{0} D_{i} + g_{0} (\boldsymbol{X}_{i}) + \varepsilon_{i}
$$
ここで，$g_{0} (\boldsymbol{X}_{i})$は，$\boldsymbol{X}_{i}$が$Y_{i}$に与える影響を表す未知関数である．線形回帰モデルは，
$$
g_{0} (\boldsymbol{X}_{i}) = \boldsymbol{X}_{i}^{\top} \boldsymbol{\beta}_{0}
$$
と置いた特殊ケースである．
つまり，部分線形回帰モデルは，$D_{i}$については線形係数$\theta_{0}$を保ちつつ，$\boldsymbol{X}_{i}$の影響については柔軟な関数$g_{0}$で表すモデルである．関心のあるパラメータはあくまでも$\theta_{0}$であり，$g_{0}$そのものではない．
このような，関心のあるパラメータではないが，それを推定するために必要になる対象は，局外母数（Nuisance Parameter）と呼ばれる．DMLでは，この局外母数の推定が重要になる．
局外母数としての条件付き期待値
部分線形回帰モデルは，$g_{0} (\boldsymbol{X}_{i})$を推定すればよいように見える．しかし，DMLのPartialling-out型の定式化では，次の2つの条件付き期待値を考える．
$$
\ell_{0} (\boldsymbol{X}_{i}) = \mathbb{E} [Y_{i} | \boldsymbol{X}_{i}]
$$
$$
r_{0} (\boldsymbol{X}_{i}) = \mathbb{E} [D_{i} | \boldsymbol{X}_{i}]
$$
$\ell_{0} (\boldsymbol{X}_i)$は，$\boldsymbol{X}_{i}$で説明できる$Y_{i}$の成分である．一方，$r_{0} (\boldsymbol{X}_{i})$は，$\boldsymbol{X}_{i}$で説明できる$D_{i}$の成分である．
部分線形回帰モデル
$$
Y_{i} = \theta_{0} D_{i} + g_{0} (\boldsymbol{X}_{i}) + \varepsilon_{i}
$$
において，$\mathbb{E}[\varepsilon_i | \boldsymbol{X}_i] = 0$を仮定すると，
$$
\mathbb{E} [Y_{i} | \boldsymbol{X}_{i}] = \theta_{0} \mathbb{E} [D_{i} | \boldsymbol{X}_{i}] + g_{0} (\boldsymbol{X}_{i})
$$
である．したがって，
$$
\ell_{0} (\boldsymbol{X}_{i}) = \theta_{0} r_{0} (\boldsymbol{X}_{i}) + g_{0} (\boldsymbol{X}_{i})
$$
となる．
ここで，$Y_{i}$と$D_{i}$をそれぞれ条件付き期待値で残差化する．
$$
\tilde{Y}_{i} = Y_{i} - \ell_{0} (\boldsymbol{X}_{i})
$$
$$
\tilde{D}_{i} = D_{i} - r_{0} (\boldsymbol{X}_{i})
$$
すると，部分線形回帰モデルは，残差同士の関係として
$$
\tilde{Y}_{i} = \theta_{0} \tilde{D}_{i} + \varepsilon_{i}
$$
と書ける．したがって，$\theta_0$は
$$
\theta_{0} = \dfrac{\mathbb{E} [\tilde{D}_{i} \tilde{Y}_{i}]}{\mathbb{E} [\tilde{D}_{i}^{2}]}
$$
として識別できる．
これは，FWL定理のノンパラメトリック版のように読める．線形回帰モデルでは，$\boldsymbol{X}_{i}$で説明できる部分を線形射影で取り除いた．一方，DMLでは，$\boldsymbol{X}_{i}$で説明できる部分を条件付き期待値として取り除く．
ナイーブなPlug-in推定の問題
ここまで読むと，次のように考えたくなる．
まず，何らかの機械学習手法で$\ell_{0} (\boldsymbol{X}_{i})$と$r_{0} (\boldsymbol{X}_{i})$を推定する．
$$
\hat{\ell} (\boldsymbol{X}_{i}) \approx \ell_{0} (\boldsymbol{X}_{i})
$$
$$
\hat{r} (\boldsymbol{X}_{i}) \approx r_{0} (\boldsymbol{X}_{i})
$$
次に，残差を作る．
$$
\hat{\tilde{Y}}_{i} = Y_{i} - \hat{\ell} (\boldsymbol{X}_{i})
$$
$$
\hat{\tilde{D}}_{i} = D_{i} - \hat{r} (\boldsymbol{X}_{i})
$$
最後に，$\hat{\tilde{Y}}_{i}$を$\hat{\tilde{D}}_{i}$に回帰する．
$$
\hat{\theta} = \dfrac{\sum_{i=1}^{n} \hat{\tilde{D}}_{i} \hat{\tilde{Y}}_{i}}{\sum_{i=1}^{n} \hat{\tilde{D}}_{i}^{2}}
$$
この発想はかなり自然である．実際，DMLの実装も大まかにはこの形をしている．しかし，ここには注意点がある．$\hat{\ell}$や$\hat{r}$を同じデータで学習し，同じデータで残差化に使うと，過剰適合の影響が$\hat{\theta}$に入り込む可能性がある．また，局外母数の推定誤差が$\hat{\theta}$に直接影響すると，通常の標準誤差や信頼区間が信用できなくなる．
DMLは，この問題に対して，主に2つの工夫を用いる．
1つ目は，ネイマン直交スコア（Neyman Orthogonal Score）を用いることである．これは局外母数の推定誤差が，関心のあるパラメータの推定に一次の影響を与えにくいような推定方程式を使う，という考え方である．
2つ目は，交差適合（Cross Fitting）を用いることである．これは，局外母数を学習するデータと，関心のあるパラメータを推定するために使うデータを分ける，という考え方である．
ネイマン直交スコア
部分線形回帰モデルでは，次のスコアを考える．
$$
\psi (W_{i}; \theta, \eta) = [Y_{i} - \ell (\boldsymbol{X}_{i}) - \theta \lbrace D_{i} - r(\boldsymbol{X}_{i}) \rbrace] \lbrace D_{i} - r(\boldsymbol{X}_{i}) \rbrace
$$
ここで，
$$
W_{i} = (Y_{i}, D_{i}, \boldsymbol{X}_{i})
$$
であり，
$$
\eta = (\ell , r)
$$
である．真の局外母数は
$$
\eta_{0} = (\ell_{0}, r_{0})
$$
である．
このスコアは，真のパラメータ$\theta_{0}$と真の局外母数$\eta_{0}$のもとで
$$
\mathbb{E} [\psi (W_{i}; \theta_{0}, \eta_{0})] = 0
$$
というモーメント条件を満たす．したがって，サンプル版の推定方程式
$$
\dfrac{1}{n} \sum_{i=1}^{n} \psi (W_{i}; \theta, \hat{\eta}) = 0
$$
を解くことで，$\theta_{0}$を推定できる．
ここで重要なのは，このスコアがネイマン直交性（Neyman Orthogonality）を満たすことである．直観的には，$\ell_{0}$や$r_{0}$を少し間違えても，その誤差が$\theta_{0}$の推定に一次の影響を与えにくい，という性質である．
もう少し砕いて書くと，DMLは局外母数を完璧に推定することを要求しているわけではない．勿論，ある程度よく推定されている必要はある．しかし，推定方程式の作り方を工夫することで，局外母数の推定誤差に対して鈍感な形にしている．
この点が，単に機械学習で$Y_{i}$を予測して，残差を作り，$D_{i}$に回帰するというナイーブな方法との違いである．
交差適合
もう1つの重要な工夫が交差適合である．これは，サンプルをいくつかのfoldに分け，各foldについて，そのfold以外のデータで局外母数を推定し，残されたfold上で予測値を作る方法である．
例えば，サンプルを$K$個に分ける．
$$
\lbrace 1, 2, \cdots , n \rbrace = I_{1} \cup I_{2} \cup \cdots I_{K}
$$
$k$番目のfoldである$I_{k}$に属する個体について予測値を作るときは，$I_{k}$を除いたデータを使って$\hat{\ell}_{-k}$と$\hat{r}_{-k}$を学習する．そして，$i \in I_{k}$について$\hat{\ell}_{-k} (\boldsymbol{X}_{i})$と$\hat{r}_{-k} (\boldsymbol{X}_{i})$を作る．
このようにすると，各個体$i$について，その個体自身を使わずに学習されたモデルから予測値が作られる．これにより，局外母数の推定量と，推定方程式に入る観測値との過度な依存を避けることができる．
この操作は，機械学習における単なる予測精度評価のための交差検証（Cross Validation）とは少し目的が異なる．DMLにおける交差適合は，予測性能を評価するためというより，過学習によるバイアスを避け，推論を正当化しやすくするために用いられる．
DML推定量
部分線形回帰モデルにおけるDML推定量は，実装上は極めてシンプルである．まず交差適合によって，すべての$i$について$\hat{\ell}_{-k(i)} (\boldsymbol{X}_{i})$と$\hat{r}_{-k(i)} (\boldsymbol{X}_{i})$を作る．ここで，$k(i)$は個体$i$が属するfoldを表す．
次に残差を作る．
$$
\hat{\tilde{Y}}_{i} = Y_{i} - \hat{\ell}_{-k(i)} (\boldsymbol{X}_{i})
$$
$$
\hat{\tilde{D}}_{i} = D_{i} - \hat{\ell}_{-k(i)} (\boldsymbol{X}_{i})
$$
最後に残差同士を回帰する．
$$
\hat{\theta}_{DML} = \dfrac{\sum_{i=1}^{n} \hat{\tilde{D}}_{i} \hat{\tilde{Y}}_{i}}{\sum_{i=1}^{n} \hat{\tilde{D}}_{i}^{2}}
$$
これは見た目は非常にシンプルだが，この単純な式の中にFWL定理，部分線形回帰モデル，ネイマン直交性，交差適合という重要な考え方が詰まっている．
Pythonによるシミュレーション
ここからはシミュレーションデータを作成し，DMLでPythonを実装する．EconMLのようなパッケージは使わずに，scikit-learnとstatsmodelsだけで実装する．
まず，次のようなデータ生成過程を考える．
$$
D_{i} = r_{0} (\boldsymbol{X}_{i}) + v_{i}
$$
$$
Y_{i} = \theta_{0} D_{i} + g_{0} (\boldsymbol{X}_{i}) + \varepsilon_{i}
$$
ここで，真の処置効果を$\theta_{0} = 2$とする．$g_{0}$と$r_{0}$は非線形関数にしておく（i.e., 線形回帰では正しくモデル化できない状況を作る）．
import numpy as np
import pandas as pd
from sklearn.base import clone
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import KFold
import statsmodels.api as sm
def make_data(n=3000, p=10, theta=2.0, seed=123):
rng = np.random.default_rng(seed)
X = rng.normal(size=(n, p))
g = (
np.sin(X[:, 0])
+ 0.5 * X[:, 1] ** 2
- 0.75 * X[:, 2] * X[:, 3]
+ 0.3 * np.exp(X[:, 4] / 2)
)
r = (
0.7 * np.sin(X[:, 0])
- 0.4 * X[:, 1]
+ 0.3 * X[:, 2] ** 2
+ 0.2 * X[:, 3] * X[:, 4]
)
v = rng.normal(size=n)
eps = rng.normal(size=n)
D = r + v
Y = theta * D + g + eps
ell = theta * r + g
return Y, D, X, g, r, ell
ここで，$g$は$g_{0} (\boldsymbol{X}_{i})$，$r$は$r_{0} (\boldsymbol{X}_{i})$，ellは$\ell_{0} (\boldsymbol{X}_{i}) = \mathbb{E}[Y_{i} | \boldsymbol{X}_{i}]$である．通常，$g$，$r$，ellは観測できないが，シミュレーションでは真の値を知っているので，比較用に用いる．
比較対象①：最小二乗法
まずは，通常の線形回帰モデルに対して最小二乗法を適用する．
$$
Y_{i} = \theta D_{i} + \boldsymbol{X}_{i}^{\top} \boldsymbol{\beta} + \varepsilon_{i}
$$
ただし，実際の$g_{0}$と$r_{0}$は非線形関数なので，この線形回帰モデルはミススペックされている．
def ols_linear(Y, D, X):
Z = np.column_stack([D, X])
Z = sm.add_constant(Z)
res = sm.OLS(Y, Z).fit(cov_type="HC1")
theta_hat = res.params[1]
se = res.bse[1]
return theta_hat, se
比較対象②：Partialling-out
次に，真の$\ell_0 (\boldsymbol{X}_{i})$と$r_{0} (\boldsymbol{X}_{i})$が既知の場合を考える．これは実行不可能な推定量であるが，DMLが目指している理想的な残差化のベンチマークになる．
def partial_out_theta(Y, D, y_hat, d_hat):
y_tilde = Y - y_hat
d_tilde = D - d_hat
theta_hat = np.mean(d_tilde * y_tilde) / np.mean(d_tilde ** 2)
psi = d_tilde * (y_tilde - theta_hat * d_tilde)
se = np.sqrt(
np.mean(psi ** 2)
/ (np.mean(d_tilde ** 2) ** 2 * len(Y))
)
return theta_hat, se
DMLの実装
次に交差適合を用いてDMLを実装する．model_yは$\ell_{0} (\boldsymbol{X}_{i}) = \mathbb{E} [Y_{i} | \boldsymbol{X}_{i}]$を推定するモデルであり，model_dは$r_{0} (\boldsymbol{X}_{i}) = \mathbb{E} [D_{i} | \boldsymbol{X}_{i}]$を推定するモデルである．
def dml_plr(Y, D, X, model_y, model_d, n_splits=5, seed=123):
n = len(Y)
y_hat = np.zeros(n)
d_hat = np.zeros(n)
kf = KFold(
n_splits=n_splits,
shuffle=True,
random_state=seed
)
for train_idx, test_idx in kf.split(X):
m_y = clone(model_y)
m_d = clone(model_d)
m_y.fit(X[train_idx], Y[train_idx])
m_d.fit(X[train_idx], D[train_idx])
y_hat[test_idx] = m_y.predict(X[test_idx])
d_hat[test_idx] = m_d.predict(X[test_idx])
theta_hat, se = partial_out_theta(Y, D, y_hat, d_hat)
return theta_hat, se, y_hat, d_hat
やっていることは単純である．各foldについて，そのfold以外のデータで$Y_{i}$と$D_{i}$をそれぞれ予測するモデルを学習する．そして，残されたfold上で予測値を作る．最後に，$Y_{i}$と$D_{i}$をそれぞれ残差化し，残差同士を回帰する．
交差適合なしの推定
比較のために，交差適合を使わず，同じデータで局外母数を学習し，同じデータで残差化する推定も行う．
def no_split_partialling_out(Y, D, X, model_y, model_d):
m_y = clone(model_y)
m_d = clone(model_d)
m_y.fit(X, Y)
m_d.fit(X, D)
y_hat = m_y.predict(X)
d_hat = m_d.predict(X)
theta_hat, se = partial_out_theta(Y, D, y_hat, d_hat)
return theta_hat, se
この方法は，見た目にはDMLにかなり近い．しかしながら，局外母数の学習と関心のあるパラメータの推定に同じデータを使っているため，過学習の影響を受けやすい．
実行
Y, D, X, g, r, ell = make_data()
rf = RandomForestRegressor(
n_estimators=300,
min_samples_leaf=10,
random_state=123,
n_jobs=-1
)
theta_ols, se_ols = ols_linear(Y, D, X)
theta_oracle, se_oracle = partial_out_theta(
Y=Y,
D=D,
y_hat=ell,
d_hat=r
)
theta_no_split, se_no_split = no_split_partialling_out(
Y=Y,
D=D,
X=X,
model_y=rf,
model_d=rf
)
theta_dml, se_dml, y_hat_dml, d_hat_dml = dml_plr(
Y=Y,
D=D,
X=X,
model_y=rf,
model_d=rf,
n_splits=5,
seed=123
)
results = pd.DataFrame(
{
"method": [
"Linear OLS",
"Oracle partialling out",
"No-split ML partialling out",
"DML with cross-fitting",
],
"theta_hat": [
theta_ols,
theta_oracle,
theta_no_split,
theta_dml,
],
"standard_error": [
se_ols,
se_oracle,
se_no_split,
se_dml,
],
}
)
results
手元で実行すると以下のような結果になる．
method theta_hat standard_error
0 Linear OLS 2.022551 0.026131
1 Oracle partialling out 1.993059 0.018463
2 No-split ML partialling out 1.954083 0.023850
3 DML with cross-fitting 1.964409 0.022979
真の値は$\theta_0 = 2$である．乱数や機械学習モデルの設定によって数値は変わるが，ここで見たいのは，DMLが「$Y_{i}$と$D_{i}$を柔軟に残差化し，残差同士を回帰する」という形で実装できる点である．
また，交差適合なしの推定とDMLは似ているが，同じものではない．特に，より複雑な学習器を使ったり，サンプルサイズが小さかったり，高次元の特徴量を使ったりすると，交差適合の有無が結果に大きく効く場合がある．DMLでは，この交差適合が過学習バイアスを抑えるための重要な役割を果たす．
なぜ結果変数だけではなく処置変数も残差化するのか
DMLで躓きやすい点の1つは，なぜ$Y_i$だけでなく$D_i$も残差化するのか，という点であるように思う．もし目的が単なる予測であれば，$Y_{i}$を$\boldsymbol{X}_{i}$でよく予測できればそれでよいかも知れない．しかし，ここでの目的は$Y_{i}$の予測ではなく，$D_{i}$の係数$\theta_{0}$の推定である．$D_{i}$が$\boldsymbol{X}_{i}$と関係している場合，$D_{i}$の変動には，$\boldsymbol{X}_{i}$によって説明される成分と，$\boldsymbol{X}_{i}$では説明されない成分が混ざっている．$\theta_{0}$を推定するために使いたいのは，$\boldsymbol{X}_{i}$では説明されない$D_{i}$の変動である．
したがって，$D_i$についても
$$
\tilde{D}_{i} = D_{i} - \mathbb{E} [D_{i} | \boldsymbol{X}_{i}]
$$
を作る必要がある．
これはFWL定理とまったく同じ発想である．重回帰分析では，$D_{i}$から$\boldsymbol{X}_{i}$で線形に説明できる成分を取り除く．DMLでは，$D_{i}$から$\boldsymbol{X}_{i}$で非線形に説明できる成分を取り除く．この意味で，DMLは重回帰分析のコントロールを機械学習で柔軟にしたものと見ることができる．ただし，単に機械学習で残差化するだけではなく，ネイマン直交スコアと交差適合によって，推論が壊れにくい形にしている点が重要である．
DMLは何をしていて何をしていないのか
ここまでの議論をまとめると，DMLは次のことをしている．
第1に，関心のあるパラメータ$\theta_{0}$と，関心はないが推定に必要な局外母数$\eta_{0}$を分ける．
第2に，局外母数を機械学習などで柔軟に推定する．
第3に，ネイマン直交スコアを使うことで，局外母数の推定誤差が$\theta_{0}$の推定に一次の影響を与えにくくする．
第4に，交差適合を使うことで，過学習によるバイアスを抑える．
一方で，DMLは次のことを自動的に解決してくれるわけではない．
まず，DMLは識別仮定を与えてくれるわけではない．因果効果として$\theta_{0}$を解釈したいなら，何を条件付ければ交絡が十分にコントロールされるのか，そもそも$D_{i}$の変動をどのようなものとして解釈するのかを考える必要がある．
また，DMLは局外母数をどんなに雑に推定してもよい，という話でもない．$\hat{\ell}$や$\hat{r}$の推定精度が極端に悪ければ，当然，$\hat{\theta}_{DML}$も信用できない．
さらに，DMLは機械学習モデルの選択から自由にしてくれるわけでもない．ランダムフォレストを使うのか，勾配ブースティングを使うのか，RidgeやLassoを使うのか，ニューラルネットワークを使うのかによって，結果が変わることはあり得る．したがって，実務では複数の学習器で結果を比較することが望ましい．
おわりに
この記事では，DMLを重回帰分析とFWL定理の延長として説明した．
通常の重回帰分析では，$Y_{i}$と$D_{i}$から$\boldsymbol{X}_{i}$で線形に説明できる部分を取り除き，残差同士の関係から$\theta_{0}$を推定する．DMLの部分線形回帰モデルでは，この操作を非線形・高次元の設定に拡張する．すなわち，$\ell_{0} (\boldsymbol{X}_{i}) = \mathbb{E} [Y_{i} | \boldsymbol{X}_{i}]$と$r_{0} (\boldsymbol{X}_{i}) = \mathbb{E} [D_{i} | \boldsymbol{X}_{i}]$を機械学習で推定し，$Y_{i}$と$D_{i}$をそれぞれ残差化したうえで，残差同士を回帰する．
ただし，単に機械学習で残差化すればよいわけではない．DMLでは，ネイマン直交スコアによって局外母数の推定誤差の影響を抑え，交差適合によって過学習によるバイアスを抑える．この2つがあるからこそ，柔軟な機械学習手法を使いながら，関心のある低次元パラメータについて通常の推論を行うことができる．
したがって，DMLは「機械学習を使えば因果推論ができる」という話ではない．むしろ，DMLは，古典的な回帰分析やセミパラメトリック推論の考え方を，現代的な機械学習手法と接続するための方法である．少なくとも部分線形回帰モデルにおいては，DMLはFWL定理の自然な拡張として理解するのが最も見通しがよいと思う．
次回は，この記事で直観的に説明したネイマン直交性と交差適合をもう少し丁寧に見ていく．特に，なぜナイーブなPlug-in推定ではダメなのか，なぜネイマン直交性があると局外母数の推定誤差に強くなるのか，標準誤差や信頼区間をどのように計算するのかを整理する．
参考文献
- Hansen, B. 2022. Econometrics, Princeton University Press.