『ビジネスデータサイエンスの教科書』6章の「不均一な処置効果」解説

Matt Taddy著『ビジネスデータサイエンスの教科書』6章の「不均一な処置効果」の節が少々分かりづらかったため，解説的な記事を書く．

目次

  • Notation
  • この節の概要
  • 順を追った要約

Notation

基本的には教科書のnotationに従うが，私の趣味で確率変数は大文字を使って表す． また，DMLのモデルは教科書よりも若干generalな形で， Y_i = D_\i \tau + m(X_i) + u_i, \quad D_i = g(X_i) + v_i というモデルで表現する．

不均一な処置効果（異質な処置効果）の概論

日経新聞に竹村『現代数理統計学』の広告を出すか出さないかを意思決定するためには，日経新聞購読者全体に対しての平均介入効果を調査すれば良い*1．つまり，日経新聞購読者を母集団として， E[Y_i(\text{広告})-Y_i(\text{出さない})]を知れば良い．しかし，web広告ではどうだろう．web広告では新聞広告と異なり，ユーザー毎に異なる広告戦略を採用することができる．例えば，Springerへのアクセス履歴がある人には広告を見せて，TikTokにしかアクセス履歴がない人には広告を見せないという戦略をとることも可能である（偏見）．ここで我々が知りたいのは， E[Y_i(\text{広告})-Y_i(\text{出さない})|\text{アクセス履歴}]である．

不均一な処置効果（Heterogeneous Treatment Effect）という言葉は，人や性質によって処置効果が異なるということを明示的に表現した言葉である*2．ここまでの文脈では平均介入効果のみにフォーカスしていたため，それと対比するために用いられている．この章はその人の介入効果 Y_i(1) - Y_i(0)をなるべく正確に予測しようという試みが行われる．よりフォーマルには， E[Y_i(1)-Y_i(0)|X_i=x]を推定することがこの章の目標である．

これまでの流れを改めて抑え直そう．完全な実験のとき，つまり Y_i(1), Y_i(0) \perp D_iのとき，ATEを求めるためには Y_i = \tau D_i + u_iというモデルの回帰式を回せばよいだけであった（このOLSの解は介入群と対照群でのアウトカムの標本平均の差と一致する）．次に，不完全な実験のとき，つまり Y_i(1), Y_i(0) \perp D_i | X_iのとき，ATEを求めるためには Y_i = \tau D_i + X_i^{\top} \beta + u_iというモデルの回帰式を回すのであった．もしくは，DMLを用いれば，はじめに D_i = g(X_i) + v_i, Y_i = m(X_i) + w_iという形で機械学習による回帰をまわして \hat{D}_i = \hat{g}(X_i), \hat{Y}_i = \hat{m}(X_i)を計算してから， Y_i - \hat{m}(X_i) = (D_i-\hat{g}(X_i)) \tau + u_iを回帰すれば良いのであった．

では，完全な実験のときにHTEを求めるためにはどうすれば良いのか．完全な実験のときにせよ，不完全な実験のときにせよ，以下の交差項を含めた回帰式を回すのが自然に思える．

\displaystyle Y_i = \tau D_i + X_i^{\top} \beta + X_i^{\top} D_i \gamma + u_i

つまり， X_i=xの人の介入効果を \tau + x^{\top} \gammaという風に推定しており，右辺第一項は平均的な介入効果，右辺第二項はその人のベースライン（私が広告を見なかったときに本を買っていたか否か），第三項はその人特有の（異質な）介入効果を意味する*3．

具体的な係数の求め方としては，（高次元な X_iのときは）これまでの流れの通りLassoかDMLを用いるというのが著者の主張である．通常のLassoの実装がp.169の最後で，DMLによる実装がp.176くらいに書かれている．DMLによる実装を簡単に書くと，はじめに D_iを X_iで予測してやって，

\displaystyle Y_i - \hat{m}(X_i) = \tau (D_i - \hat{g}(X_i)) + X_i^{\top} (D_i - \hat{g}(X_i)) \gamma + u_i

を回帰する．

「不均一な処置効果」の節の要約

p.167-169中盤まで: 欠損値の話．正直kaggle本の内容で十分な気がする．
  • カテゴリ変数 -欠損している値はNAというカテゴリとして扱う．
  • 数値変数
    • 0がNA以外の値のなかで半分以上を占めていた場合は0埋め．疎行列の性質を保てることがハッピー．
    • そうでない場合は平均埋め
    • 他の変数から予測して埋めるのも興味深い方法だぜ
p.169-p.170: 交差項モデルを通常のLassoを用いて回帰する．

\displaystyle Y_i = \tau D_i + X_i^{\top} \beta + X_i^{\top} D_i \gamma + u_i

あまりにも見づらいが，dが介入 D_iに対応して，xhteが高次元特徴量 X_i，dxhteは交差項 D_i X_iに対応している．gamという変数はdとdxhteの係数のみを見ていて，round(sort(gam)[1,6], 4)は係数の小さい変数を6つ取り出している．例えば，d.race_asian_12mYesの係数が-0.04であることは，アジア人に介入を受けさせても平均的な介入効果（次のコードで出てくるdの係数）より4\%ほど介入効果が小さいことを示唆している．

p.170-171: 価格弾力性の話．

式6.16がどれくらい妥当なのかがよく分からない．

p.171- スーパーマーケット『ドミニクス』でのビールの売上量の話
  • データの特性
    • 各ビールの製品コードのについて各週の単位売上量を63店舗分集めたデータ
    • 店舗コード, 製品コード, 週番号, 価格, 単位売上量のデータ
    • 各製品コードに関して，内容量，短い説明書きのデータ
    • 観測値としては160万以上あるが，シミュレーションのためにそこから5000行だけ抽出する．残りの160万のデータに対しては，機械学習の手法とか使わずに，OLSを用いて弾力性を求めて，あたかも答えであるかのように扱う．暴論だと思う．
  • まず，平均的な価格弾力性を求める p.172最後-p.175
    • 最初に \log(q_{i}) = \log(p_{i}) \tau + u_{i}で回帰するが，弾力性が高く，信憑性がないと著者は主張
    • 価格と売上の両方に影響を与えていそうな，店舗コード・製品コード・週番号・商品の説明文に含まれている単語を統制変数 xとみなして \log(q_{i}) = \log(p_{i}) \tau + X_i^{\top} \beta + u_{i}でLasso回帰する．この結果にも著者は納得しない．
    • DMLを用いて，はじめに Xから \log q_{i}と \log p_{i}を予測して， \log(q_{i}) - \hat{m}(X_i) = (\log(p_{i})-\hat{g}(X_i))\tau + X_i^{\top} \beta + u_{i}でLasso回帰する．この結果には，ご満悦のようである．
  • 次に，各製品の異質な価格弾力性を推定する（ここが本題） p.175-p.178
    • ここからは統制変数 xは短い説明書きのデータのみを考え，店舗コード, 製品コード, 週番号は変数から落とす．
    • 例にもれず，交差項で考える． \log(q_{i}) = \log(p_{i}) \tau + X_i^{\top} \beta + X_i^{\top} D_i \gamma +u_{i}
    • とりまDMLを使う． \log(q_{i}) - \hat{m}(X_i) = (\log(p_{i})-\hat{g}(X_i))\tau + X_i^{\top} \log(p_{i}) \beta + X_i^{\top} (\log(p_{i})-\hat{g}(X_i)) \gamma+ u_{i}
    • p.177以降で，二個上の交差項モデルをLassoで回帰するパターンとOLS（最尤法）で回帰するパターンも実施．160万行のデータを用いて「答え」とされている値と比較してDMLこそ神だと主張．
    • ちなみに．図6.3-6.5はそれぞれのビールの価格弾力性の推定値のヒストグラムで，図6.6はそれぞれのビールの価格弾力性の推定値とそれぞれのビールの価格弾力性の答え（160万行からの推定値）をプロットしたもので， y=xの直線上に載っているのが望ましい．なんで R^2で評価しようとしたんや．

*1:厳密には日経新聞購読者の人数などの情報も必要だろうが．

*2:個人的な肌感覚としては，性別や年齢といった比較的低次元なグループレベルで条件付けたときの介入効果はGroup Average Treatment Effect(GATE)と表現され，webデータの情報や遺伝子情報など高次元の情報で条件付けるときはConditional Average Treatment Effect(CATE)やIndividual Treatment Effect(ITE)という用語が用いられるように思われる．HTEはGATEやCATEを内包した表現だと思っている．

*3:ここで X_i^{\top} \betaが介入効果に含まれないのは因果推論で大事なことで，仮に私が広告を見なくても本を買うのであれば（換言すれば， x^{\top} \betaが大きな値であれば），私に広告を見せる必要は特にないかもしれない．

Doubly Robust の直感的解釈（を書こうとしたけど無理だった）

因果推論において、Doubly Robust法とは、傾向スコアまたはresponse functionのどちらかがmisspecifiedされていても、CATE(またはATEなど何かしらの介入効果)の不偏推定量を求めることができる手法です。

Notation, Definition, Assumption

$$ \begin{align} &\text{共変量（covariate）}:X \in \mathbb{R}^p \\ &\text{介入（treatment）}: T \in {0,1}\\ &\text{potential outcome}: Y(t) \in \mathbb{R}\\ &\text{被説明変数（outcome）}: Y = Y(1)T - Y(0)(1-T) \in \mathbb{R}\\ &\text{傾向スコア（propensity score）}:p(x) = E[T=1|X=x]\\ &\text{反応関数(?)（response function）}: \mu_t(x) = E[Y|T=t,X=x] = E[Y(t)|T=t,X=x]\\ &\text{CATE（Conditional Average Treatment Effect）}: E[Y(1)|X=x] - E[Y(0)|X=x] = \mu_1(x) - \mu_0(x) \end{align} $$

なお、以下では、Conditional unconfoundedness: $Y(1),Y(0) \perp T | X$およびoverlap condition: $\forall x :0<p(x)<1$は満たされていると仮定します。この仮定のもとでは、$E[Y(1)T|X=x] = E[Y(1)|X=x]E[T|X=x]=E[Y(1)]$

Doubly Robust推定量

Doubly Robust推定量とは、先述の通り「傾向スコアまたはresponse functionのどちらかがmisspecifiedされているときでも成立する、CATEの不偏推定量」です。その性質さえ満たしていればDoubly Robustなので、Doubly Robustな推定量自体はたくさん開発されていますが、その根幹になる式は以下のものだと理解しています。 $$ \begin{align} E[Y(1)|X=x] &= E\left[ \frac{TY}{p(X)} | X= x\right] - E \left [\frac{T-p(X)}{p(X)} \mu_1(X) | X=x \right] \quad (1) \end{align} $$

(1)のsample analogue（サンプルでの平均）をとってやると、いわゆるAugumented Inverse Probability Weightingに帰着します。この式が、$p(x),\mu_1(x)$のどちらかが間違っていたとしても成立することを今から説明していきます。

\eqref{1}の右辺を書き換えてみましょう。 $$ \begin{align} &E\left[ \frac{TY}{p(X)} | X= x\right] - E \left [\frac{T-p(X)}{p(X)} \mu_1(X) | X=x \right] \\ = &\frac{ E\left[ TY | X= x\right]}{p(x)} - \frac{E \left [T-p(x) | X=x \right] }{p(x)} \mu_1(x) \quad (2) \\ = &\mu_1(x) - \frac{E[T(Y-\mu_1(x))|X=x]}{p(x)} \quad (3) \end{align} $$

この変形は$p(x),\mu_1(x)$が特定化されていなくても成立することに注意して下さい。つまり、たとえ$\mu_1(x)$や$p(x)$がでたらめな関数でも成立する変形です。

では、$p(x)$は正しく特定化されており、$\mu(x)$はでたらめな関数だったとしましょう。でたらめであることを明示的にするために、$\tilde{\mu}_1(x)$と書きます。(2)の第一項は

$$ \begin{align} \frac{E[TY|X=x]}{p(x)} &= \frac{E[T(TY(1) - (1-T)Y(0))|X=x]}{p(x)} \\ &= \frac{E[TY(1)|X=x]}{p(x)} \\ &= \frac{E[T|X=x]}{p(x)}E[Y(1)|X=x] \\ &= E[Y(1)|X=x] \end{align} $$

となり、第2項が0になれば\eqref{1}が成立しそうです。では、第二項を見てみると $$ \begin{align} \frac{ E\left[ T-p(x)| X= x \right]}{p(x)} \tilde{\mu}1(x) = \frac{p(x) - p(x)}{p(x)} \tilde{\mu}1(x) = 0 \end{align} $$

となり、確かに成立します！すごい！

つぎに、$p(x)$はでたらめな$\tilde{p}(x)$で、$\mu(x)$は正しく特定化されているとしましょう。でたらめとは言いましたが、一応、$\tilde{p}(x)>0$は満たされているとします。すると、今度は(3)に注目してやると、第一項は$E[Y(1)|X=x]$なので、第二項が0に帰着してほしくなります。すると、実際に $$ \begin{align} \frac{E[T(Y-\mu_1(x))|X=x]}{\tilde{p}(x)} &= \frac{E[T|X=x]E[Y-\mu_1(x)|X=x]}{\tilde{p}(x)}\\ &= \frac{p(x) \cdot 0}{\tilde{p}(x)} = 0 \end{align} $$

しゅごい！

阪急：大阪梅田駅から京阪：淀屋橋駅まで歩く

阪急梅田駅から京阪淀屋橋まで歩いたので、真の最短ルートではないかもしれませんが、備忘録の役割も含めて共有します。 間違いなく、以下の動画を見るのが最も分かりやすいので、家を出る前にこの動画を見ておいて、あとは私が以下に書くメモを見ながら実際に歩けば良いのではないでしょうか。


www.youtube.com


www.youtube.com

大まかな流れとしては、大阪駅前第3ビルに行く→地上に上がって南下するです。

大阪駅前第3ビルまで

少々古いですが、写真付きの こちら もご参照下さい。

阪急梅田駅３階改札（電車のホームと同じ階の改札！）

→エスカレーターで下る

→ムービングウォーク

→南北コンコースを直進（阪急梅田百貨店を右手に歩く）

→エスカレーターで地下へと下る

→ホワイティに入り、すぐに左折→すぐに右折して地下鉄東梅田駅

→すぐ目に入る北東改札ではなく、ずっと直進して南改札手前の右手階段（8番出入口）を上がる

→大阪駅前第4ビルに入る→そのまま3ビル奥までいく

大阪駅前第3ビルから淀屋橋駅まで

（3ビルの奥まで歩いた続き）

→吉野家らへんで外にでる。さらに地下への階段（梅田新道横断地下道）があるので、そこを下って、直進して、地上に上がる。大和証券のでっかいビルがあればok。

→大江橋駅が途中にあるが無視して直進→地下鉄・京阪淀屋橋駅がすぐに見える

ただし、目的地によっては、大江橋で乗って京橋乗り換えという作戦もあります。