Double/Debiased Machine Learning

Double/Debiased Machine Learning#

Paper: Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., & Robins, J. (2018). Double/debiased machine learning for treatment and structural parameters.
Python Package: DoubleML

概要#

motivation#

世の中の多くの現象は非線形な関係性が想定される。回帰分析は線形モデルであるため、モデルの定式化の誤りに起因するバイアスが生じかねない。

実際に関心のあるパラメータは少なく、交絡のコントロールのために入れている局外母数（nuisance parameters）は高次元になりがち。

局外母数を非線形モデルで表し、関心のあるパラメータは線形モデルで表現する部分線形モデル

\[ Y = \theta D + g(X) + e \]

を作り、非線形モデル$g(X)$を機械学習で構築したい

先行研究#

ノンパラメトリックモデルの文脈で研究があるが、Donsker条件を満たすような、複雑性の低い関数クラスでないとうまくいかないという理論解析がある

機械学習のようなDonsker条件を満たさない複雑な関数でも使えるようにしたい

課題#

$g(X)$を機械学習で作って線形回帰するだけだと推定量は$\sqrt{n}$の収束レートにならない

2つのバイアスがある

正則化バイアス（Reguralization bias）
過学習によるバイアス（Bias induced by overfifting）

それぞれの対策として、

正則化バイアス → 残差回帰（直交化）で対応
過学習 → Cross-Fittingで対応

を行う

前提知識#

知ってると理解を深めやすくなる前提知識まとめ

（参考）残差回帰#

DMLでは、FWL定理を利用した残差回帰を一般化する。

残差回帰#

通常の残差回帰は線形回帰モデル

\[ Y = X_1 \beta_1 + X_2 \beta_2 + e \]

を用いて

\[ Y - X_2 \hat{\gamma} = (X_1 - X_2 \hat{\delta}) \beta_1 \]

のようにして関心のあるパラメータ$\beta_1$を推定する。

線形回帰モデルは回帰関数$E[Y|X]$を近似するため、上記の残差回帰は期待値を用いて

\[ Y - E[Y|X_2] = (X_1 - E[X_1|X_2]) \beta_1 \]

と表すことができる。

部分線形モデルへの拡張#

部分線形モデル

\[ Y = \theta D + g(X) + e \]

においても同様に

\[ Y - E[Y|X] = (D - E[D|X]) \beta_1 \]

とするアプローチをとる

FWL定理

目的変数のベクトル$Y \in \mathbb{R}^{n\times 1}$、説明変数の行列$X \in \mathbb{R}^{n\times d}$と誤差項$e \in \mathbb{R}^{n\times 1}$による線形回帰モデル

\[ Y = X \beta + e \]

があるとする。

説明変数を$X = (X_1 | X_2)$と2つのグループに分割し、回帰係数ベクトル$\beta$も合わせて$\beta = (\beta_1 | \beta_2)^T$と2つに分割して

\[ Y = X_1 \beta_1 + X_2 \beta_2 + e \]

と表す。

この回帰モデルの$\beta_1$は、残差回帰（residual regression）と呼ばれる以下の手順に従うことでも得ることができる。

$X_1$を$X_2$に回帰して残差$\tilde{X}_1$を得る：$\tilde{X}_1 = X_1 - X_2 \hat{\delta}$
$Y$を$X_2$に回帰して残差$\tilde{Y}$を得る：$\tilde{Y} = Y - X_2 \hat{\gamma}$
$\tilde{Y}$を$\tilde{X}_1$に回帰させる：$\tilde{Y} = \tilde{X}_1 \beta_1$

（参考）モーメント法#

確率変数$X$が$\theta$というパラメータをもつ分布に従うとする。

\[ E[\psi(X; \theta)] = 0 \]

という条件（直交条件）を満たすスコア関数$\psi(X; \theta)$があるとき、標本$X_1, \dots, X_n$を使った直交条件

\[ \frac{1}{n} \sum^n_{i=1} \psi(X_i; \theta) = 0 \]

を解いて$\theta$を推定する方法を モーメント法 （method of moments）という。

ナイーブな推定量と正則化バイアス#

話を簡潔にするために、サンプルを2つにランダム分割するとしよう：メインパートは$n$サンプルで$i\in I$のインデックスで表し、補助パートは$N-n$として$i \in I^c$とする。単純のため$n=N/2$とする。補助パートのサンプルで$\hat{g}_0$を獲得し、$\theta_0$の推定にメインパートのサンプルを使うことにする

\[ \hat{\theta}_0 = \left( \frac{1}{n} \sum_{i\in I} D_i^2 \right)^{-1} \frac{1}{n} \sum_{i\in I} D_i (Y_i - \hat{g}_0(X_i)) \tag{1.3} \]

推定量$\hat{\theta}_0$は一般に$1/\sqrt{n}$より遅い収束レート、つまり

\[ |\sqrt{n} (\hat{\theta}_0 - \theta_0)| \overset{p}{\to} \infty \]

となる。

この「劣った」振る舞いの背後には$g_0$の学習におけるバイアスがある。

ヒューリスティックにこの$\hat{g}_0$の学習のバイアスのインパクトを説明すると、スケールされた$\hat{\theta}_0$の推定誤差は

\[ \sqrt{n} (\hat{\theta}_0 - \theta_0) = \underbrace{ \left( \frac{1}{n} \sum_{i\in I} D_i^2 \right)^{-1} \frac{1}{\sqrt{n}} \sum_{i\in I} D_i U_i }_{:=a} + \underbrace{ \left( \frac{1}{n} \sum_{i\in I} D_i^2 \right)^{-1} \frac{1}{\sqrt{n}} \sum_{i\in I} D_i (g_0(X_i) - \hat{g}_0(X_i)) }_{:=b} \]

である。

第1項$a$は$a \rightsquigarrow N(0, \bar{\Sigma})$となるので問題ない。第2項の$b$の項は正則化バイアス項で、一般に中心にならず発散する。first orderで以下を得る

\[ b = (E[D_i^2])^{-1} \frac{1}{\sqrt{n}} \sum_{i\in I} m_0(X_i)(g_0(X_i) - \hat{g}_0(X_i)) + o_P(1) \]

ヒューリスティックには、$b$は平均がゼロにならない$m_0(X_i)(g_0(X_i) - \hat{g}_0(X_i))$の$n$個の総和で、$\sqrt{n}$で割られる。これらの項は非ゼロの平均になる。なぜなら高次元あるいは複雑な問題設定のために我々はlasso, ridge, boostingあるいはpenalized neural netsなどの正則化推定量を採用するためである。

それらの推定量において正則化は推定量の分散が爆発しないようにするものの相当なバイアスを引き起こす。とりわけ、$g_0$への$\hat{g}_0$のバイアスの収束レートは、RMSEにおいて$n^{-\phi_g}$（$\phi_g < 1/2$）である。ゆえに$b$は$D_i$が$m_0(X_i)\neq 0$で中心化されるとき$\sqrt{n} n^{-\phi_g} \to \infty$の確率的オーダーになることが期待され、よって$|\sqrt{n} (\hat{\theta}_0 - \theta_0)| \overset{p}{\to} \infty$となる

直交化による正則化バイアスの打破#

直交化されたregressor $V=D-m_0(X)$を得るためにDからXの効果をpartialling outすることにより得られる直交化された式（orthogonalized formulation）を使う。具体的には $$ \hat{V} = D - \hat{m}_0(X) $$ ここで$\hat{m}_0$は補助的サンプル（auxiliary sample）からを用いたML推定量。

DからXの効果をpartialling out したあとは、main sampleを使って$\theta_0$のdebiased ML (DML)を構築する $$ \check{\theta}_0 = \left( \frac{1}{n} \sum_{i\in I} \hat{V}_i D_i \right)^{-1} \frac{1}{n} \sum_{i\in I} \hat{V}_i (Y_i - \hat{g}_0(X_i)) $$ 近似的にDをXについて直交化し、$g_0$の推定値を引くことで近似的に交絡の直接効果を除去することで、$\check{\theta}_0$は(1.3)の正則化バイアスを除去している。

$\check{\theta}_0$の性質#

スケールされた推定誤差は３つの要素に分解できる

\[ \sqrt{n}(\check{\theta}_0 - \theta_0) = a^* + b^* + c^* \]

$a^*$はmild conditionsのもとで以下を満たす

\[ a^* = (E[V^2])^{-1} \frac{1}{\sqrt{n}} \sum_{i\in I} V_i U_i \rightsquigarrow N(0, \Sigma) \]

$b^*$は$g_0, m_0$の推定における正則化バイアスの影響を捉える。具体的には

\[ b^* = (E[V^2])^{-1} \frac{1}{\sqrt{n}} \sum_{i\in I} (\hat{m}_0(X_i) - m_0(X_i)) (\hat{g}_0(X_i) - g_0(X_i)) \]

で、これは$\hat{m}_0$と$\hat{g}_0$の推定誤差の積に依存する。そのため、幅広い範囲のデータ生成過程のもとで消失させることが可能である。

実際、この項は$\sqrt{n}n^{-(\phi_m+\phi_g)}$で上界になり、ここで$n^{-\phi_m}, n^{-\phi_g}$はそれぞれ$\hat{m}_0, \hat{g}_0$の$m_0, g_0$への収束レートである。これは両者が比較的遅い収束レートで推定されたとしても、消失しうる。

$c^*$は

\[ c^* = o_P(1) \]

となる。これが弱い条件のもとで成り立つことを保証するにあたって、sample splittingが重要な役割を果たす。

Cross Fittingによる過学習のバイアスの除去#

sample splittingは$c^*$が確率的に消失するために使われる。

$c^*$は $$ \frac{1}{\sqrt{n}} \sum_{i \in I} V_i (\hat{g}_0(X_i) - g_0(X_i)) $$ などの項を含む。この項は推定誤差と構造的未観測要因の積の和を$1/\sqrt{n}$-normalizedしたものである。

$E[V_i|X_i]=0$であることを思い出すと、この項は平均ゼロで分散は $$ \frac{1}{n} \sum_{i \in I} (\hat{g}_0(X_i) - g_0(X_i))^2 \overset{p}{\to} 0 $$ であることがわかり、チェビシェフの不等式を使って確率的には消失することがわかる。

sample splittingの欠点は推定に使用するサンプル数が減ることによる効率性の低下である。しかし、mainとauxiliaryの2つでそれぞれ推定を行い、両者の平均を取ればfull efficiencyを取り戻す。この手続き（mainとauxiliaryの役割を取り替えて複数の推定値を取得しそれらの平均をとる）を「cross-fitting」と呼ぶことにする。一般にk-foldにすることもできる。Section 3で詳細を述べる。

Neyman Orthogonality#

partialling outしたモーメント条件による推定量とナイーブな推定量は何が違うのか？

→ ネイマン直交性（Neyman orthogonality）がカギになる。

ネイマン直交性

サンプル$W$、関心のあるパラメータ$\theta_0$、局外母数$\eta_0$についてのスコア関数$\psi(W; \theta_0, \eta_0)$のベクトル$\psi = (\psi_1, \dots, \psi_d)^T$があるとする。このスコア関数の直交条件

\[ E[\psi(W; \theta_0, \eta_0)] = 0 \]

について、ガトー微分が存在し、

ネイマン直交性（Neyman orthogonality）

We show that if the population risk satisfies a condition called Neyman orthogonality, the impact of the nuisance estimation error on the excess risk bound achieved by the meta-algorithm is of second order.

Foster, D. J., & Syrgkanis, V. (2023). Orthogonal statistical learning. The Annals of Statistics, 51(3), 879-908.

モーメント条件#

モーメント条件と非線形モデル#

GMM推定量は非線形モデルに対しても用いることができる（末石2015, p.78）

部分線形モデル#

結果変数$Y$、説明変数$X, Z$についての次のようなモデルを**部分線形モデル（Partially Linear Model）**という。

\[\begin{split} \begin{align} Y &= X^T \beta + g(Z) + u\\ E[u|X,Z] &= 0\\ E[u^2|X,Z] &= \sigma^2(X, Z)\\ \end{align} \end{split}\]

ここで、$g(\cdot)$は任意の関数（非線形の関数でもよい）である。$\sigma^2(\cdot, \cdot)$も未知の関数で、不均一分散を許容する。

関心のあるパラメータ$\beta$を解釈性の高い線形モデルで推定（パラメトリック推定）しつつ、影響を統制するためだけにモデルに投じている$Z$は関数形を特定せず推定（ノンパラメトリック推定）することができるため、部分線形モデルはセミパラメトリックモデルと呼ばれる。

Robinson (1983)の推定量#

$Z$で両辺の条件付き期待値をとった

\[ E[Y|Z] = E[X|Z]^T \beta + g(Z) \]

を差し引くと

\[ Y - E[Y|Z] = (X - E[X|Z])^T \beta + u \]

$\tilde{Y}_i = Y_i - E[Y_i|Z]$、$\tilde{X}_i = X_i - E[X_i|Z]$とおけば、OLS推定量の形になる

\[ \hat{\beta}_{\mathrm{inf}}=\left(\sum_{i=1}^n \tilde{X}_i \tilde{X}_i^{\prime}\right)^{-1} \sum_{i=1}^n \tilde{X}_i \tilde{Y}_i \]

が$E[Y|Z], E[X|Z]$が未知のため実行不可能である。$E[Y|Z], E[X|Z]$をそれぞれのノンパラメトリック推定量で置換して推定を行うことは可能であり、その推定量はroot-N consistentで漸近正規性を持つ

参考#

末石直也 - セミ・ノンパラメトリック計量分析 (京都大学) 第5回部分線形モデルとセミパラメトリック推定量の性質
西山・人見（2023）『ノン・セミパラメトリック統計解析』、共立出版

（疑問）Robinsonの推定量とDMLの違い、そうなった原因は？#

Robinsonの推定量とDMLの違い#

期待値を差し引くというFWL-style or Neiman-orthogonalityを使うのは一緒
モーメント条件が議論の根幹なのでOLS推定量の形には限らない？？
差し引く期待値の推定をMLにしたのがDML

新規性は#

ML推定ゆえoverfitting biasが入る → Cross Fitting
- なんでMLだとbiasがはいる？高次元や複雑なモデルはDonsker conditionを満たさないから？

4. Inference in Partially Linear Models

PLR model

\[\begin{split} \begin{aligned} Y&=D \theta_0+g_0(X)+U, & E_P[U \mid X, D]=0 \\ D&=m_0(X)+V, & E_P[V \mid X]=0 \end{aligned} \end{split}\]

の$\theta_0$の推定方法は2つある

1つめは

\[ \psi(W ; \theta, \eta):=\{Y-D \theta-g(X)\}(D-m(X)), \quad \eta=(g, m) \]

というスコア関数のもとでのDML

2つ目はRobinsonの”partialling-out”のスコア関数

\[ \psi(W ; \theta, \eta):=\{Y-\ell(X)-\theta(D-m(X))\}(D-m(X)), \quad \eta=(\ell, m) \]

を使うもの

Donsker条件#

機械学習モデルの収束レートが遅い（Donsker条件）→Cross Fitting

XユーザーのMasahiro Katoさん: 「Double/debiased machine learningはどういう手法かというと，機械学習はDonsker条件を満たさないので，サンプルを分割することによって，異なる部分集合間で独立だと思えるような局外母数の推定量を構築することで，収束率だけで漸近正規性を出す手法です．」 / X

考察#

機械学習を使うと常にバイアスが入る？？#

X-learnerとかどうなる？

参考#

DMLによるDID#

Double/debiased machine learning for difference-in-differences models | The Econometrics Journal | Oxford Academic

講義動画（Youtube）#

Double Machine Learning for Causal and Treatment Effects - YouTube

MLでのcausal parametersの推定は良いとは限らない
double or orthogonalized MLとsample splittingによって、causal parametersの高精度な推定が可能

Partially Linear Modelを使う

\[ Y = D \theta_0 + g_0(Z) + U , \hspace{1em} E[U|Z, D] = 0 \]

MLをそのまま使うと一致推定量にならない（例えば$Y - D$で$g_0(Z)$をRandom Forestで学習しても、予測精度は良いがバイアスがある）
FWL定理を用いて、残差の回帰にするとよい

\[\begin{split} \hat{W} = Y - \hat{E[Y|Z]}\\ \hat{V} = D - \hat{E[D|Z]} \end{split}\]

モーメント条件

Regression adjustment: $E[(Y - D \theta_0 - g_0(Z) ) D] = 0$
“propensity score adjustment”: $E[(Y - D \theta_0) (D - E[D|Z])] = 0$
Neyman-orthogonal (semi-parametrically efficient under homoscedasticity): $E[(\hat{W} - \hat{V}\theta_0) \hat{V}] = E[\{(Y - E[Y|Z]) - (D - E[D|Z])\theta_0\} (D - E[D|Z])] = 0$

3は不偏

Sample Splitting

Splittingによるefficiencyの低下問題

2個に分けて2回やって平均とればfull efficiency → k個に分けての分析をk回やって平均とってもfull efficiency

応用研究#

[2002.08536] Debiased Off-Policy Evaluation for Recommendation Systems

PR: AI Lab、推薦システム分野におけるトップカンファレンス「RecSys2021」にて共著論文採択ー高次元なデータを使った意思決定評価の信頼性を改善ー | 株式会社サイバーエージェント

部分線形モデルに対するモーメント条件#

異なる推定量でどれだけバイアスが入るか確認したい

import numpy as np

def g(X):
    # linear
    return 2 * X

n = 1000
np.random.seed(0)
X = np.random.uniform(size=n)
theta = 3
D = X + np.random.normal(size=n)
Y = D * theta + g(X) + np.random.normal(size=n)

import pandas as pd
df = pd.DataFrame({"Y": Y, "D": D, "X": X})
X_mat = X.reshape(-1, 1)
D_mat = D.reshape(-1, 1)

ナイーブな推定量

\[ E[UD] = E[(Y - D\theta - g(X)) D] = 0 \]

\[ \hat{\theta}_0=\left(\frac{1}{n} \sum_{i \in I} D_i^2\right)^{-1} \frac{1}{n} \sum_{i \in I} D_i\left(Y_i-\hat{g}_0\left(X_i\right)\right) \]

かりに、$\hat{g}(X_i)$は$Y_i$を$X$に回帰する（予測値$\hat{Y}_i = \hat{g}(X_i)$を作る）とするなら

from lightgbm import LGBMRegressor

# Y ~ X
g = LGBMRegressor(max_depth=4, verbose=-1)
g.fit(X_mat, Y)

# Y_res := Y - g(X)
Y_res = Y - g.predict(X_mat)
theta_hat = np.mean(D**2) ** (-1) * np.mean(D * Y_res)
print(f"θ={theta_hat:.3f}")

# Y_res ~ D のOLS
import statsmodels.formula.api as smf
final_model = smf.ols(
    formula='Y_res ~ -1 + D',
    data=df.assign(
        Y_res = Y - g.predict(X_mat)
    )
).fit()
final_model.summary().tables[1]

θ=1.981

	coef	std err	t	P>\|t\|	[0.025	0.975]
D	1.9811	0.049	40.446	0.000	1.885	2.077

from sklearn.linear_model import LinearRegression

# gをDonsker条件を満たしそうな線形モデルにしたらどうなるか
# これは普通に残渣回帰でもないしだめか
# Y ~ X
g = LinearRegression()
g.fit(X_mat, Y)

# Y_res := Y - g(X)
Y_res = Y - g.predict(X_mat)
theta_hat = np.mean(D**2) ** (-1) * np.mean(D * Y_res)
print(f"θ={theta_hat:.3f}")

θ=2.177

g.coef_

array([4.65807719])

# TODO: モーメント条件の方向微分をplotできないものか？
# 横軸はthetaとr

疑問：$\hat{g}(X)$はどう得る？

$Y = g(X)$で学習させる？でもそうなると$E[Y|X] = \ell(X)$を使うRobinsonとの違いは…？

DML score function#

ネイマン直交性を満たす

\[ E[UV] = E[ \{Y-D \theta-g(X)\}(D-m(X)) ]=0 \]

$\hat{V}= D-\hat{m}_0(X)$

\[ \check{\theta}_0=\left(\frac{1}{n} \sum_{i \in I} \widehat{V}_i D_i\right)^{-1} \frac{1}{n} \sum_{i \in I} \widehat{V}_i\left(Y_i-\hat{g}_0\left(X_i\right)\right) \]

from lightgbm import LGBMRegressor
from sklearn.model_selection import cross_val_predict

# Y ~ X
g = LGBMRegressor(max_depth=4, verbose=-1)
g.fit(X_mat, Y)

# D ~ X
m = LGBMRegressor(max_depth=4, verbose=-1)
m.fit(X_mat, D)

# y_res := Y - g(X)
# V_hat := D - m(X)
V_hat = D - m.predict(X_mat)
Y_res = Y - g.predict(X_mat)

theta_hat = np.mean(V_hat * D) ** (-1) * np.mean(V_hat * Y_res)
print(f"θ={theta_hat:.3f}")

import statsmodels.formula.api as smf
final_model = smf.ols(
    formula='Y_res ~ -1 + V_hat',
    data=df.assign(
        Y_res = Y_res,
        V_hat = V_hat,
    )
).fit()
final_model.summary().tables[1]

θ=2.872

	coef	std err	t	P>\|t\|	[0.025	0.975]
V_hat	2.9672	0.032	93.561	0.000	2.905	3.029

Robinson-style “partialling-out” score function#

Robinson (1988)の推定量

部分線形モデル

\[ Y = D \theta_0 + g(X) + U , \quad E[U|D,X]=0 \]

の両辺を$X$で条件づけて期待値をとると

\[ E[Y|X] = E[D|X] \theta_0 + g(X) \]

これをモデルから差し引くと

\[ Y - E[Y|X] = \theta_0 (D - E[D|X]) + U \]

という線形回帰の形になる。

ただし、$E[Y|X], E[D|X]$は未知なのでそれぞれノンパラメトリック推定量$\hat{\ell}(X), \hat{m}(X)$で置き換える。

こちらもネイマン直交性を満たす。

\[ \mathrm{E}\left[\left\{(Y-E[Y \mid X])-(D-E[D \mid X]) \theta_0\right\}(D-E[D \mid X])\right]=0 \]

$V=D - E[D|X] = D-m_0(X)$として、$\ell_0(X) := E[Y|X]$とすると

\[ \mathrm{E}\left[\left((Y-\ell_0(X))-V \theta_0\right)V\right]=0 \]

cross fitting#

# cross-fittingあり
from lightgbm import LGBMRegressor
from sklearn.model_selection import cross_val_predict

m = LGBMRegressor(max_depth=6, verbose=-1)
df["V_hat"] = df["D"] - cross_val_predict(m, df[["X"]], df["D"], cv=5)

g = LGBMRegressor(max_depth=6, verbose=-1)
df["Y_res"] = df["Y"] - cross_val_predict(g, df[["X"]], df["Y"], cv=5)

import statsmodels.formula.api as smf
final_model = smf.ols(formula='Y_res ~ -1 + V_hat', data=df).fit()
final_model.summary().tables[1]

	coef	std err	t	P>\|t\|	[0.025	0.975]
V_hat	2.9471	0.032	92.429	0.000	2.885	3.010