Causal Embeddingsの解説と追試

はじめに

前回は, ログデータの観測確率が一様ではない場合に傾向スコアで補正した損失関数を用いるPropensity Matrix Factorizationを紹介しました. しかし, 2018年のRecsysにてBest Paper Awardを受賞したCausal Embeddings for Recommendation [Bonner+ 2018]で, Propensity MFを実験的に上回る手法が提案されました. 本記事では人工データを用いて, その提案手法の追試を行います.

目次

定式化のおさらい

推薦アルゴリズムの学習の定式化をおさらいします. しかし, ここでは[Bonner+ 2018]のモデル化に合わせるため, 前回記事とは少し異なる部分があることに注意してください.

各ユーザーを u \in \{1, ..., U \} = \mathcal{U}, 各アイテムを i \in \{ 1, ...., I \} = \mathcal{I}とします. また, P_{u,i}をユーザー uがアイテム iに対して有する真のPreferenceとします. P_{u,i}はユーザーとアイテムの組み合わせにのみ依存して決定的な値を取るとします. 推薦アルゴリズムの学習によって達成したいのは, 次のように定義される真の損失関数 \mathcal{L}を最小化するような予測値集合 \hat{P} = \{ \hat{P}_{u, i} \}_{(u,i) \in \mathcal{U} \times \mathcal{I}} を得ることです.

\begin{aligned} \mathcal{L} \left( \hat{P} \right) = \frac{1}{U \cdot I} \sum_{u,i} \delta \left( P_{u,i}, \hat{P}_{u,i} \right) \end{aligned}


ここで, \delta (\cdot, \cdot): \mathbb{R} \times \mathbb{R} \rightarrow \mathbb{R}は真のRatingとそれに対する予測値を入力とする関数でした.

さらに, 自分たちが観測できるPreferenceと観測できないそれを区別するために, Recommendationの頭文字をとって新たな確率変数 Rを導入します. この Rは2値確率変数であり各ユーザーとアイテムのペアに対応して独立に存在します. R_{u,i} = 1 ならば, 過去にユーザー uにアイテム iが推薦されたことを意味し, そのペアについてのPreferenceが観測されることを意味します. [Bonner+ 2018]では, 観測されるPreferenceを Y_{u,i}^{obs} とした時, 次のようなモデルを想定しています.

\begin{aligned} Y_{u,i}^{obs} = P_{u,i} \cdot R_{u,i} \end{aligned}


このモデル化では, 過去に推薦が発生したペア（ R_{u,i} = 1）については, Y^{obs}_{u,i} = P_{u,i} となり真のPreferenceが観測されますが, それ以外のペアについては Y_{u,i}^{obs} = 0となりPreferenceが観測されないということになります.（ユーザーとアイテムのMatrixを考えた時に, Preferenceが観測されず欠損となってしまっている部分を0で表しているというイメージです.）

以降, 推薦有無を表す確率変数 R_{u,i}は議論の中で重要な役割を担うので, その期待値にnotationを用意しておきましょう.

\begin{aligned} \mathbb{E} \left [ R_{u,i} \right ] = \pi_c (u,i) \end{aligned}


\pi_c (u,i)はログデータが集められた際に走っていた推薦方策を表し, logging policyと呼びます. これを用いると私たちが学習の際に用いることができる（真の損失関数の推定に用いることができる）学習データセット \mathcal{S}を次のように表すことができます.

\begin{aligned} \mathcal{S} = \left\{ \left( u, i, Y_{u,i}^{obs} \right): R_{u,i} = 1 \right\} \end{aligned}

Radomized Policy下でのNaive Loss

前章で定義した真の損失関数 \mathcal{L}を知ることができるならばそれ最適化するパラメータを得れば良いのですが, それは不可能です. よって手元にあるログデータ \mathcal{S}から損失関数を推定して, それを信じた上で, 最適化を実行する必要があります.

まずNaiveな損失関数 \hat{\mathcal{L}}_{naive} \left( \hat{P} \right)は次のように定義されます. これは, 過去に推薦が発生したデータについての損失を単純に平均するという非常に単純な推定量です.

\begin{aligned} \hat{\mathcal{L}}_{naive} \left( \hat{P} \right) = \frac{1}{ | \mathcal{S} | } \sum_{(u,i) \in \mathcal{S}} \delta \left( Y^{obs}_{u,i}, \hat{P}_{u,i} \right) \end{aligned}


[Bonner+ 2018]は, このNaive推定量 \hat{\mathcal{L}}_{naive} がある条件下で望ましい性質を満たすことに注目します. まず, ある条件とは次のようなものです.

\begin{aligned} \pi_c (u,i) = \pi_{rand} = \frac{C \cdot | \mathcal{S} |}{U \cdot I} \quad (constant) \end{aligned}


つまり, logging policyがユーザーとアイテムに非依存で, 一様な定数であるという要求です. このようなlogging policyを特にrandomized policyと呼ぶことにし, \pi_{rand}で表します. 定数部分は実はなんでも良いのですが, ここでは後の式変形が綺麗な形になる定数を置いています. さて, このlogging policyがrandomized policyであるという仮定が成り立っているときに, Naive推定量 \hat{\mathcal{L}}_{naive}の期待値をとってみましょう.

\begin{aligned} \mathbb{E}_R \left[ \hat{\mathcal{L}}_{naive} \left( \hat{P} \right) \right ] & = \mathbb{E}_R \left[ \frac{1}{ | \mathcal{S} | } \sum_{(u,i) \in \mathcal{S}} \delta \left( Y^{obs}_{u,i}, \hat{P}_{u,i} \right) \right] \\ & = \mathbb{E}_R \left[ \frac{1}{ | \mathcal{S} | } \sum_{(u,i) \in \mathcal{U} \times \mathcal{I} } R_{u,i} \cdot \delta \left( P_{u,i}, \hat{P}_{u,i} \right) \right] \\ & = \frac{1}{ | \mathcal{S} | } \sum_{(u,i) \in \mathcal{U} \times \mathcal{I}} \mathbb{E}_{R_{u,i}} \left[ R_{u,i} \right] \cdot \delta \left( P_{u,i}, \hat{P}_{u,i} \right)\\ & = \frac{1}{ | \mathcal{S} | } \sum_{(u,i) \in \mathcal{U} \times \mathcal{I}} \pi_{rand} \cdot \delta \left( P_{u,i}, \hat{P}_{u,i} \right) \\ & = \frac{1}{ | \mathcal{S} | } \sum_{(u,i) \in \mathcal{U} \times \mathcal{I}} \left( \frac{C \cdot | \mathcal{S} |}{U \cdot I} \right) \cdot \delta \left( P_{u,i}, \hat{P}_{u,i} \right) \\ & = C \left( \frac{1}{U \cdot I} \sum_{u,i} \delta \left( P_{u,i}, \hat{P}_{u,i} \right) \right) \\ & = C \cdot \mathcal{L} \left( \hat{P} \right) \\ & \propto \mathcal{L} \left( \hat{P} \right) \end{aligned}


これにより, randomized policyで集められたログデータを用いて計算されるNaiveな損失関数は真の損失関数に比例します. よって, \hat{\mathcal{L}}_{naive}を最適化するような学習は妥当であると言えそうです.

Causal Embeddings (CausE)

[Bonner+ 2018]の提案手法であるCausEは, 前章で紹介したNaive推定量の性質を活用します. CausEは学習データが次のように表せることを要求します.

\begin{aligned} \mathcal{S} = \mathcal{S}_t \cup \mathcal{S}_c \end{aligned}


ここで, \mathcal{S}_t はrandomized policyによって集められた学習データで, \mathcal{S}_cは uと iに依存するlogging policyによって集められた学習データを表します. ただし, \mathcal{S}_tは \mathcal{S}_cよりも要素数の少ない集合で良いです. また, それぞれ異なるembeddingsによって構成される2つの予測値集合 \hat{P}_t = \{ \theta_{u, t}^{\top} \beta_{i, t} \}_{ (u,i) \in \mathcal{U} \times \mathcal{I} } , \hat{P}_c = \{ \theta_{u, c}^{\top} \beta_{i, c} \}_{ (u,i) \in \mathcal{U} \times \mathcal{I} } を用いた次の損失関数をCausEは最適化します.

\begin{aligned} \hat{\mathcal{L}}_{CausE} & = \underbrace{ \frac{1}{| \mathcal{S}_t |} \sum_{ (u, i) \in S_t } \delta \left( P_{u,i}, \theta_{u, t}^{\top} \beta_{i, t} \right) }_{(1)} + \underbrace{ \frac{1}{| \mathcal{S}_c |} \sum_{ (u, i) \in S_c } \delta \left( P_{u,i}, \theta_{u, c}^{\top} \beta_{i, c} \right) }_{(2)} \\ & + \underbrace{ \lambda \left( \sum_{u \in \mathcal{U}} \| \theta_{u, t} - \theta_{u, c} \|_p + \sum_{i \in \mathcal{I}} \| \beta_{i, t} - \beta_{i, c} \|_p \right)}_{(1)} \end{aligned}


(1)はrandomized policyによって構成された \mathcal{S}_tに対するNaiveな損失（Biasなし）, (2)はlogging policyによって構成された \mathcal{S}_cに対するNaiveな損失（Biasあり）, (3)は(2)によって得られるembeddings ( \theta_{u, c}, \beta_{i, c} ) が望ましい損失である(1)から得られるembeddings ( \theta_{u, t}, \beta_{i, t} ) と乖離することに対してかける正則化で, \lambdaはその強さを調節するハイパーパラメータです. （ \theta_{u, c}, \beta_{i, c} を更新するときのみ, (3)の正則化を考慮に入れるアルゴリズムになっています.）

確かにこの損失関数からは, Biasがあるが考慮できるデータが多い損失(2)から得られたembeddingsが, Biasのない損失(1)から得られたembeddingsと離れないようにしようというお気持ちが読み取れます. しかし, この正則化について理論的な考察が論文でなされているわけではなく, どこかもどかしい気もします.

さて, 最終的な予測値集合は次のように作ります.（これが実験的に最も良かったらしいです.）

\begin{aligned} \hat{P}_{CausE} = \{ \theta_{u, c}^{\top} \beta_{i, c} \}_{ (u,i) \in \mathcal{U} \times \mathcal{I} } \end{aligned}

簡易実験

CausEとPropensity Matrix Factorizationの性能を同様の人工データを用いて比較します.

実験設定

設定に関しては前回の実験と同じです.

しかし, 前回記事ではあえて書いていなかったのですが学習データの約3%はrandomized policy ( \pi_{rand}) によって生成されています. CausEは, この少量の \mathcal{S}_tを明示的に活用することで精度向上を目指します.

また, 各手法について, Optimizerはエポック数1024のMini-batch Momentum, 学習率には初期値を0.1としたexponential decayを施しています. CausEの \mathcal{S}_tと \mathcal{S}_cで得られたembeddingsの乖離に対してかける正則化はL1 normを用いました（L2 normよりも若干良かった）.

実験結果

評価方法は前回と同様で, 各手法の25,000エポックの内, validation dataの評価が最も良かったエポックのtest dataに対するMSEの相対改善度を性能評価として用います. ただし, test dataのratingの平均で全埋めしたときのMSEを相対改善度のベースラインとしています. （[Bonner+ 2018]と同じような評価方法です.）

validation時の評価指標は, NaiveなMatix FactorizationはNaiveな評価指標を, Propensity Matrix FactorizatioはUnbiasedな評価指標をそれぞれ対応する指標として用いました. CausEに関しては, validationの際も少量の S_tを用いています.

結果は次の通りです. Best on Valはvalidation dataの評価指標が最も良かったときのエポックのtest dataに対する性能で, Best on Testは, test dataを対する真の性能が最も良かったときのエポックのtest dataに対する性能です. この2つの差が大きい場合, validationに使っていた評価指標が信頼できるものではなかった, ということです.

Best on Val Best on Test
Naive MF - 66.05 % 18.71 %
Propensity MF 16.22 % 19.31 %
Causal Embeddings (CausE) 23.02 % 23.11 %

Naive MFとPropensity MFは, 前回記事と同様の結果を掲載しています. CausEは, 少量の \mathcal{S}_tを明示的に活用することで精度を向上できていることがわかります. CausEの学習の様子は次のようになりました. 青い線は \mathcal{S}_tに対するNaiveな損失ですが, 学習データの3%ほどから計算しているのでエポック間でのブレが大きくなっています.

f:id:usaito:20190416073151p:plain
CausEのlearning curve

さいごに

これまではexplicit feedback recommendationのBias除去の話題について触れてきました. これからはimplicit feedback recommendationに関する同様の話題についても触れていきたいです.

参考

[Bonner+ 2019] Stephen Bonner and Flavian Vasile. Causal embeddings for recommendation: An Extended Abstract. arXiv:1904.05165. 2019.
[Bonner+ 2018] Stephen Bonner and Flavian Vasile. Causal embeddings for recommendation. In Proceedings of the 12th ACM Conference on Recommender Systems, pages 104– 112. ACM, 2018.
[Schnabel+ 2016] Tobias Schnabel, Adith Swaminathan, Ashudeep Singh, Navin Chandak, and Thorsten Joachims. Recommendations as treatments: Debiasing learning and evaluation. In Proceedings of the 33rd International Conference on International Conference on Machine Learning - Volume 48, ICML’16, pages 1670–1679, 2016.

因果推論で推薦システムを問い直す（学習アルゴリズム編）

はじめに

以前, こちらに本記事の評価指標編を書きました. 今回は, 同様の問題が推薦アルゴリズムの学習時にも発生し得ることを指摘し, その解決方法について議論します. 評価指標編を読んでいただいている方は, 重複する内容も多いのですんなり読んでいただけると思います.

目次

Toy Example

問題のイメージを持っていただくために, 評価指標編と同様の例を載せておきます.

今, Horror Lovers・Romance Loversというユーザー属性とHorror・Rommance・Dramaという3つのジャンルのみが存在するシンプルな映画推薦の系を考えます. さらに各ユーザーがそれぞれの映画に対して付すRatingがユーザー属性と映画ジャンルの組み合わせのみによって決まる（つまりRatingの付き方は高々6通り）とします. この時, 予測1 ( \hat{Y}_1) と予測2 ( \hat{Y}_2)（例えば予測1は, Horror LoversがRomance映画に対してRating 1という予測を出力するという見方です）の性能を自分たちが所持している推薦ログデータからMAEを用いて評価し, その評価（損失）に基づいて予測を最適化することを考えます.

f:id:usaito:20190416063747p:plain
[Schnabel+ 2016]のFigure 1より

ここで重要なのが評価に用いるログデータの生成過程です. この例では, ログデータが右上の過去の推薦方策と名付けられている行列のように各ユーザー属性と映画ジャンルの組み合わせごとに異なる推薦確率を持った方策によって集められたログデータであるとします. これにより, 例えば3200個のRatingが含まれるログデータについては各ユーザーと映画ジャンルの観測回数は図のようになり, 要素ごとの観測回数がかなり異なることがわかります. この観測回数の不均衡はひとえにこのデータが各ユーザー属性と映画ジャンルの組み合わせごとに異なる推薦確率を持った方策によって集められたことに起因します.

f:id:usaito:20190416063816p:plain
[Schnabel+ 2016]のFigure 1より

このような各要素ごとに観測回数（生成確率）が一様ではないデータを推薦システムの評価に使ってしまうと図のように, 真の性能は予測1の方が良い（MAE = 0.67）にも関わらず, ログデータ上で計算されたMAEは予測2の方が良くなってしまうという非常に厄介な問題が生じます. つまり, naiveに損失関数を設計してしまうと, 誤った方向に学習が進んでしまう可能性があるのです.

推薦アルゴリズムの定式化

前章の例で見た問題の原因を探るために, 推薦アルゴリズムの学習を定式化します. 今各ユーザーを u \in \{1, ..., U \} = \mathcal{U}, 各アイテムを i \in \{ 1, ...., I \} = \mathcal{I}とします. また, Y_{u,i}をユーザー uがアイテム iに対して有する真のRatingとします. 本記事では, Y_{u,i}がユーザーとアイテムの組み合わせによって決定的な値を取るとします（仮にRatingを確率変数として扱った場合, 以降の議論が同様に成り立つためには, 追加的な仮定が必要になります. ）.

推薦アルゴリズムの学習によって達成したいのは, 次のように定義される真の損失関数 \mathcal{L}を最小化するような予測値集合 \hat{Y} = \{ \hat{Y}_{u, i} \}_{(u,i) \in \mathcal{U} \times \mathcal{I}} を得ることです.

\begin{aligned} \mathcal{L} \left( \hat{Y} \right) = \frac{1}{U \cdot I} \sum_{u,i} \delta \left( Y_{u,i}, \hat{Y}_{u,i} \right) \end{aligned}


ここで, \delta (\cdot, \cdot): \mathbb{R} \times \mathbb{R} \rightarrow \mathbb{R}は真のRatingとそれに対する予測値を入力とする関数で, 例えばMSEを評価指標とするなら \delta(x, y) = (x - y)^{2} とすれば良いです.

もちろん, 真の損失関数 \mathcal{L}を知ることができるならばそれを用いてアルゴリズムのパラメータを最適化すれば良くなんら困難は生じません. しかし, 私たちは全てのユーザーとアイテムのインタラクションを観測できるわけではなく, 手元にあるログデータから損失関数を推定して, それを信じた上で, 最適化を実行するしかありません.

ここで自分たちが観測しているRatingと観測していないそれを区別するために新たな変数 O を導入します. この Oは2値確率変数であり各ユーザーとアイテムのペアに対応して独立に存在します. O_{u,i} = 1 ならば, ユーザー uとアイテム iのペアについてのRatingを観測していることを意味します. また, O_{u,i}の期待値は, uが iに対して有する真のRating Y_{u,i}が観測される確率を意味します.

以降, どのようにして真の損失関数の代理として妥当な, そして手持ちのログデータから計算可能な損失関数を構築すべきかについて議論します.

Naive Lossに存在するBias

まず最初に紹介するのが, Naiveな損失関数 \hat{\mathcal{L}}_{naive} \left( \hat{Y} \right)です. これは, 観測されているデータについての損失を単純に平均するという非常に単純な推定量と言えますが, 多くの論文ではこのNaiveな損失関数の推定量をいかにして最適化するかに力点が置かれているように思います.

\begin{aligned} \hat{\mathcal{L}}_{naive} \left( \hat{Y} \right) = \frac{1}{ \sum_{u,i} O_{u,i} } \sum_{(u,i): O_{u,i} = 1} \delta \left( Y_{u,i}, \hat{Y}_{u,i} \right) \end{aligned}


このNaiveな推定量は O_{u,i} = 1つまり, ログデータとして観測されているRatingに対する損失を単純に平均しています. ごくごく自然な損失関数の推定方法だと思いますが, 実はログデータの Y_{u,i}が観測される確率が一様ではない場合, このNaiveな損失関数の推定量 \hat{\mathcal{L}}_{naive}は真の損失関数 \mathcal{L} に対してBiasを持つことが次のように示されます.

今分母の \sum_{u, i} O_{u, i} を任意の正の定数 Nで置き換えた上で, \hat{\mathcal{L}}_{naive}の確率変数 Oについての期待値をとります.

\begin{aligned} \mathbb{E}_O \left[ \hat{\mathcal{L}}_{naive} \left( \hat{Y} \right) \right] & = \mathbb{E}_O \left[ \frac{1}{ N } \sum_{(u,i): O_{u,i} = 1} \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right) \right] \\ & = \mathbb{E}_O \left[ \frac{1}{ N } \sum_{u,i} O_{u,i} \cdot \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right) \right] \\ & = \sum_{u,i} \frac{\mathbb{E}_{O_{u,i}} \left[ O_{u,i} \right]} {N} \cdot \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right) \end{aligned}


ここで, \hat{\mathcal{L}}_{naive}\left( \hat{Y} \right) \propto \mathcal{L} \left( \hat{Y} \right) が成り立つためには, 任意の uと iついて

\begin{aligned} \frac{\mathbb{E}_{O_{u,i}} \left[ O_{u,i} \right]} {N} = \frac{C}{U \cdot I} \end{aligned}


が成り立つ必要があります ( Cは何かしらの正の定数.). しかし観測確率が一様ではない場合, この条件を満たす定数 Nは存在しません. よって, \hat{\mathcal{L}}_{naive} \left( \hat{Y} \right) \neq \mathcal{L} \left( \hat{Y} \right) であり, naive推定量は真の評価値に対してBiasを持つ（比例関係にない）ことがわかりました. また, そのBiasは観測確率が一様ではないことに起因することも示唆されました.

UnbiasedなLossの構築

前章では, naiveな推定量がBiasを持ってしまうという悲しい事実を紹介しましたが, 評価指標の構築時と同様にIPS (Inverse Propensity Score)と呼ばれる方法でBiasの除去が可能です. ここで, Propensity Score (傾向スコア) とは uが iに対して有する真のRating Y_{u,i}が観測される確率のことです. 表記の簡略化のため, Propensity Scoreを

\begin{aligned} P_{u,i} = \mathbb{E} \left[ O_{u,i} \right] \end{aligned}


と置くと, IPSを用いた推定量 \hat{\mathcal{L}}_{IPS} \left( \hat{Y} \right)は, 次のように定義されます.

\begin{aligned} \hat{\mathcal{L}}_{IPS} \left( \hat{Y} \right) = \frac{1}{ U \cdot I } \sum_{(u,i): O_{u,i} = 1} \frac{ \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right)} {P_{u,i}} \end{aligned}


summationが O_{u,i} = 1 で条件付けられたデータに対するものであることから, この \hat{\mathcal{L}}_{IPS} \left( \hat{Y} \right) は観測可能なデータから推定可能な統計量です. 少し天下り的に与えてしまった \hat{\mathcal{L}}_{IPS} \left( \hat{Y} \right) ですが, その期待値が真の評価値 \mathcal{L} \left( \hat{Y} \right) に一致することを次のように示すことができます.

\begin{aligned} \mathbb{E}_O \left[ \hat{\mathcal{L}}_{IPS} \left( \hat{Y} \right) \right] & = \mathbb{E}_O \left[ \frac{1}{ U \cdot I } \sum_{(u,i): O_{u,i} = 1} \frac{ \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right)} {P_{u,i}} \right] \\ & = \mathbb{E}_O \left[ \frac{1}{ U \cdot I } \sum_{u,i} O_{u,i} \cdot \frac{ \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right)} {P_{u,i}} \right] \\ & = \frac{1}{ U \cdot I } \sum_{u,i} \frac{\mathbb{E}_{O_{u,i}} \left[ O_{u,i} \right]} {P_{u,i}} \cdot \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right) \\ & = \frac{1}{ U \cdot I } \sum_{u,i} \delta \left( Y_{u,i}, \hat{Y}_{u, i} \right) \\ & = \mathcal{L} \left( \hat{Y} \right) \end{aligned}


よって Biasの観点で議論すると, 真の損失に対して比例関係を満たさないnaive推定量 \hat{\mathcal{L}}_{naive}ではなくて, 不偏性を満たすIPS推定量を最適化することでアルゴリズムを学習するのが良さそうと言えます.

Propensity Matrix Factorization

[Schnabel+ 2016]は, 真の損失関数に対する不偏性を満たす \hat{\mathcal{L}}_{IPS} を用いてユーザーとアイテムのembeddings （ \{ \theta_u \}_{u \in \mathcal{U}} , \{ \beta_i \}_{i \in \mathcal{I}}）を得るというPropensity Matrix Factorizationを提案しました. これは, 次のような式でembeddingsを逐次更新することにより, Unbiasedな勾配に基づいたembeddingsを獲得することができます. （実際には正則化項も入れます.）

\begin{aligned} \theta_u & \leftarrow \theta_u - \nabla_{\theta} \hat{\mathcal{L}}_{IPS} \left( \hat{Y} \right) \\ \beta_i & \leftarrow \beta_i - \nabla_{\beta} \hat{\mathcal{L}}_{IPS} \left( \hat{Y} \right) \end{aligned}

簡易実験

Propensity Matrix FactorizationとNaiveなMatrix Factorizationの性能を人工データを用いて比較します.

実験設定

まず, 各ユーザーとアイテムに対応する20次元のベクトルを独立にサンプリングします. \lambdaは適当に0.5としました.

\begin{aligned} \theta_u \sim \mathcal{N} \left( {\bf 0}, \lambda {\bf I} \right), \; \beta_i \sim \mathcal{N} \left( {\bf 0}, \lambda {\bf I} \right) \end{aligned}


次にユーザーバイアス b_uとアイテムバイアス b_iを独立にサンプリングします. こちらも \lambdaは0.5に設定しました.

\begin{aligned} b_u \sim \mathcal{N} \left( 0, \lambda \right), \; b_i \sim \mathcal{N} \left( 0, \lambda \right) \end{aligned}


これらのサンプリングされたベクトルやスカラー値を用いてRatingを次のように生成したあとで, 整数値に直し, Ratingの最小値が1, 最大値が5になるような処理を加えました. 末尾のrating_biasは, ratingの全体平均を調整するために入れています. 今回の実験では3に設定しています.

\begin{aligned} Y_{u,i} = \theta_u^{\top} \beta_i + b_u + b_i + rating\_bias \end{aligned}


コードで書くとこのような感じになります. 

# theta, beta, user_bias, item_biasからratingを生成.
preference = theta @ beta.T + user_bias + item_bias.T + rating_bias
rate = np.clip(np.round(preference), a_min=1., a_max=5.)

最後に傾向スコアを次のように生成します. ここでは, 真のRatingが高いほどデータが観測されやすいという考えに基づき次の生成過程を用いました. 実験では, k=0.275, \alpha=0.2に設定しています. これでsparcityがおおよそ4%となりました.

\begin{aligned} P_{u,i} = k \cdot \alpha^{ \left( 5 - Y_{u,i} \right) } \end{aligned}


上記のような過程で生成されたRatingの分布は次の通りです. 学習データのRatingの平均は, である一方で, テストデータのRatingの平均は, です. 大きなズレが生じていることが一目瞭然です.

f:id:usaito:20190416061820p:plain
TrainとTestのRating分布

実験結果

Naiveな損失を使ったMatrix FactorizationとIPSによるUnbiasedな損失を用いたPropensity Matrix Factorizationを比較します.

両者25,000エポックの内, validation dataの評価が最も良かったエポックのtest dataに対するMSEの相対改善度を性能評価として用います. ただし, test dataのratingの平均で全埋めしたときのMSEを相対改善度のベースラインとしています. （[Bonner+ 2018]と同じような評価方法です.）

validation時の評価指標は, NaiveなMatix FactorizationはNaiveな評価指標を, Propensity Matrix FactorizatioはUnbiasedなMatrix Factorizationをそれぞれ対応する指標として用いています（評価指標について詳しくはこちら）.

結果は次の通りです. Best on Valはvalidation dataの評価指標が最も良かったときのエポックのtest dataに対する性能で, Best on Testは, test dataを対する真の性能が最も良かったときのエポックのtest dataに対する性能です. つまり, この2つの差が大きい場合, validationに使っていた評価指標が信頼できるものではなかった, ということになります.

Best on Val Best on Test
Naive MF - 66.05 % 18.71 %
Propensity MF 16.22 % 19.31 %

validationを対応する評価指標で行なった場合, Propensity MFが大きく上回っています. 一方で, test dataに対して最適だったときの性能はさほど大きな差はありません. 損失関数のBiasを取り除くことが性能向上に寄与することに加えて, 評価指標がいかに重要であるかを物語る結果となりました. さらにNaive MFをNaiveな評価指標で評価するという組み合わせは, baselineよりも悪い結果となってしまう可能性があることも示唆されました.

それぞれの学習の様子は次の通りです. Naive MFは途中でtest dataに対する性能が上昇し始めてしまっていることが見えます. 一方で, Propensity Matrix Factorizationの方はほぼ単調にtest dataに対する損失が減少しています. ちなみに, Propensity Matrix Factorizationの損失は観測確率の逆数で重み付けすることもあって分散が大きくなる（≒各エポックで損失の推定がぶれる）ので, 今回はバッチサイズが1024のミニバッチ学習を用いています.

f:id:usaito:20190416061647p:plain
Naive Matrix Factorizationのlearning curve

f:id:usaito:20190416061708p:plain
Propensity Matrix Factorizationのlearning curve

さいごに

本記事では, 簡単な例から入り理論事実と簡易実験により推薦システムの学習をnaiveな損失関数に基づいて行ってしまうと学習過程でBiasが生じること, 傾向スコアを用いたIPS推定量によりそのBiasを除去し予測性能を向上できることを示してきました. 実は, 2018年のRecsysでPropensity Matrix Factorizationを実験的に上回った手法が提案されています. 次回は, その手法の追試を行っていくつもりです.

参考

[Bonner+ 2018] Stephen Bonner and Flavian Vasile. Causal embeddings for recommendation. In Proceedings of the 12th ACM Conference on Recommender Systems, pages 104– 112. ACM, 2018.
[Schnabel+ 2016] Tobias Schnabel, Adith Swaminathan, Ashudeep Singh, Navin Chandak, and Thorsten Joachims. Recommendations as treatments: Debiasing learning and evaluation. In Proceedings of the 33rd International Conference on International Conference on Machine Learning - Volume 48, ICML’16, pages 1670–1679, 2016.

Domain Adversarial Neural Networksの解説

はじめに

最近自分の研究分野との親和性が高いこともあり, Unsupervised Domain Adaptation (教師なしドメイン適応)の理論を勉強しています. その理論を応用した手法に, Domain Adversarial Neural Networks (DANN) というものがあり自分でも動かしてみました.

目次

  • Unsupervised Domain Adaptationとは
  • H-divergenceを用いた汎化誤差上界
  • Domain Adversarial Neural Networks
  • 簡易実験
  • さいごに

Unsupervised Domain Adaptationとは

まずUnsupervised Domain Adaptation (UDA) を定式化します. 入力空間を \mathcal{X}, 出力空間を \mathcal{Y} = \{0, 1\}とします. ここで, あるdomain Dとは, 入力の分布 P_Dとlabeling function f_D: \mathcal{X} \rightarrow \mathcal{Y}のpair \left( P_D, f_D \right)のことを指します. UDAは, labelが得られているsource domainからのサンプル \mathcal{S} = \{ x_i, y_i \}_{i=1}^{n_S} とlabelが得られていないtarget domainからのサンプル \mathcal{T} = \{ x_j \}_{j=1}^{n_T} から, 次のtarget domainにおける期待判別誤差をできるだけ小さくするような仮説 hを得ることを目指します.

\begin{aligned} R_T^l (h, f_T) = \mathbb{E}_{X \sim P_T} \left[ l \left( h(X), f(X) \right) \right] \end{aligned}


以降は, 損失関数 lとして0-1損失のみを考えます. この時,

\begin{aligned} R_T (h, f_T) = \mathbb{E}_{X \sim P_T} \left[ \left| h(X) - f(X) \right| \right] \end{aligned}


です. このように表されるtarget domainにおける期待判別誤差を（一様に）boundしたいというのが理論的なモチベーションになります. 学習データとテストデータのDomainが同一であるような通常の教師あり機械学習の場合, 次のような形で予測判別誤差をboundするのが一般的です.

任意の h \in \mathcal{H}について, 少なくとも 1 - \delta ( \forall \delta \in (0, 1)) の確率で次の不等式を満たす.

\begin{aligned} R_D (h, f_D) \leqq \hat{R}_D (h, f_D) + complexity \left( \mathcal{H} \right) + confidence\_level (\delta) \end{aligned}


しかし, UDAではテストデータのlabelがサンプルとして得られていないので, target domainの経験判別誤差を用いることができません. よって, 次のようなboundを得ることを目指すこととします.

任意の h \in \mathcal{H}について, 少なくとも 1 - \delta ( \forall \delta \in (0, 1)) の確率で次の不等式を満たす.

\begin{aligned} R_T (h, f_D) \leqq & \underbrace{\hat{R}_S (h, f_S)}_{(1)} + \underbrace{complexity \left( \mathcal{H} \right)}_{(2)} + \underbrace{confidence\_level (\delta)}_{(3)} \\ & + \underbrace{descrepancy (P_S, P_T)}_{(4)} + \underbrace{difference (f_S, f_D)}_{(5)} \end{aligned}


(1)は, Source Domainにおける経験判別誤差. (2)は, 仮説集合の複雑さ. (3)は \deltaに依存する項. (4)はSourceとTargetの入力分布の乖離度. (5)はSourceとTargetのlabeling functionの乖離度です. 次節では, H-divergenceと呼ばれるdiscrepancyを用いたboundについて説明します.

H-divergenceを用いた汎化誤差上界

H-divergenceは次のように定義されるdiscrepancyの一種です.

\begin{aligned} d_{\mathcal{H}} \left( P_S, P_T \right) & = 2 \sup_{h \in \mathcal{H}} \left| R_S (h, 1 ) - R_T (h, 1 ) \right| \\ & = 2 \sup_{h \in \mathcal{H}} \left| \left( 1 - R_S (h, 0 ) \right) - R_T (h, 1 ) \right|\\ & = 2 \sup_{h \in \mathcal{H}} \left| 1 - \left( R_S (h, 0 ) + R_T (h, 1 ) \right) \right| \end{aligned}

仮説集合が対称であるとき, empiricalには

\begin{aligned} d_{\mathcal{H}} \left( \hat{P}_S, \hat{P}_T \right) = 2 \left( 1 - \min_{h \in \mathcal{H}} \left( \frac{1}{n_S} \sum_{i = 1}^{n_S}\mathbb{I} \left[ h(x_i) = 0 \right] + \frac{1}{n_T} \sum_{j = 1}^{n_T} \mathbb{I} \left[ h(x_j) = 1 \right] \right) \right) \end{aligned}


です. よって, H-divergenceは仮説集合 \mathcal{H}がsource domainとtarget domainのデータを入力から判別する性能に依存することがわかります. 入力からどちらのDomainからのサンプルかが判別できるほど, 2つのDomainのH-Divergenceが大きくなるというイメージです.

このH-divergenceを用いると次のようなboundが得られます. ([Ganin+ 2015]のTheorem 2, [Ben David+ 2010]のTheorem 2を参考にした)

任意の h \in \mathcal{H}について, 少なくとも 1 - \delta ( \forall \delta \in (0, 1)) の確率で次の不等式を満たす.

\begin{aligned} R_T (h, f_D) \leqq & R_S (h, f_S) + \frac{1}{2} d_{\mathcal{H}} \left( P_S, P_T \right) + \beta \end{aligned}


ただし, \beta \geqq \inf_{h \in \mathcal{H}} \left( R_S(h, f_S) + R_T (h, f_T) \right)です. つまり, \betaは R_S(h, f_S)と R_T(h, f_T)の和の下限の上界です.


さらに, 既存の統計的機械学習における結果を用いてempricalに推定可能なboundを次の通りに得ます.

\mathcal{H}を有限のVC次元 dを持つ仮説集合とする. 任意の h \in \mathcal{H}について, 少なくとも 1 - \delta ( \forall \delta \in (0, 1)) の確率で次の不等式を満たす. ただし, n_S = n_T = nとした.

\begin{aligned} R_T (h, f_D) \leqq & \hat{R}_S (h, f_S) + \sqrt{ \frac{4}{n} \left( d \log \left( \frac{2en}{d} \right) + \log \left( \frac{4}{\delta} \right) \right) } \\ & + \frac{1}{2} d_{\mathcal{H}} \left( \hat{P}_S, \hat{P}_T \right) + 4 \sqrt{ \frac{1}{n} \left( 2d \log \left( 2n \right) + \log \left( \frac{4}{\delta} \right) \right) } + \beta \end{aligned}


VC次元やそれを用いた予測判別誤差と経験判別誤差の差の一様boundについては, MLPシリーズ『統計的機械学習』のChapter 2に説明があります.

Domain Adversarial Neural Networks (DANN)

ようやくDANNの説明に入ります. 前節で得た R_T (h, f_T)の上界のうち, 私たちがどうにかできるのはsource domainにおける経験判別誤差 \hat{R}_S (h, f_S)と P_S, P_Tの経験H-Divergence d_{\mathcal{H}} \left( \hat{P}_S, \hat{P}_T \right) です. これらの和を小さくするために3つのlayer G_f, G_y, G_dを考えます. G_f (\cdot ; \theta_f)はrepresentation layerと呼び, 入力空間 \mathcal{X}をある望ましい空間 \mathcal{R}に写像する役割を持ちます. G_y (\cdot ; \theta_y)はprediction layerと呼び, G_fで得た特徴表現 \mathcal{R}からlabelを予測します. 最後に, G_d (\cdot ; \theta_d)はdomain layerと呼び, 特徴表現 \mathcal{R}からSource Domainから得られたサンプルなのかTarget Domainで得られたサンプルなのかを判別します.

これらを用いてDANNの損失関数は次のように定義されます.

\begin{aligned} E \left( \theta_f, \theta_y, \theta_d \right) = \underbrace{ \frac{1}{n_S} \sum_{i = 1}^{n_S} \mathcal{L}^i_y \left( \theta_f, \theta_y \right)}_{(1)} - \underbrace{ \lambda \left( \frac{1}{n_S} \sum_{i = 1}^{n_S} \mathcal{L}^i_d \left( \theta_f, \theta_d \right) + \frac{1}{n_T} \sum_{j = 1}^{n_T} \mathcal{L}^j_d \left( \theta_f, \theta_d \right) \right) }_{(2)} \end{aligned}


ここで,

\begin{aligned} & \mathcal{L}^i_y \left( \theta_f, \theta_y \right) = \mathcal{L}_y \left( G_y \left( G_f (x_i, \theta_f), \theta_y \right), y_i \right) \\ & \mathcal{L}^i_d \left( \theta_f, \theta_d \right) = \mathcal{L}_d \left( G_d \left( G_f (x_i, \theta_f), \theta_d \right), d_i \right) \end{aligned}


はそれぞれサンプル iに対するprediction lossとdomain lossです. 損失関数 E \left( \theta_f, \theta_y, \theta_d \right) のうち, (1)はsource domainにおける経験判別誤差を表しており(2)は, G_fによって生成される表現 \mathcal{R}上での経験H-Divergenceと読めます. \lambdaはそのどちらをどれだけ重視するかを司るハイパーパラメータです. 要は E \left( \theta_f, \theta_y, \theta_d \right) は, R_T(h, f_T)の上界のうち私たちがどうにかできる項と言えます. 3つのパラメータ \theta_f, \theta_y, \theta_gはそれぞれ次のように更新します．

\begin{aligned} \theta_f & \leftarrow \theta_f - \mu \cdot \left( \frac{\partial \mathcal{L}_y^i }{\partial \theta_f} - \lambda \frac{\partial \mathcal{L}_d^i }{\partial \theta_f} \right) \\ \theta_y & \leftarrow \theta_y - \mu \cdot \frac{\partial \mathcal{L}_y^i }{\partial \theta_y} \\ \theta_d & \leftarrow \theta_d - \mu \cdot \frac{\partial \mathcal{L}_d^i }{\partial \theta_d} \end{aligned}


\muは学習率です. 3つの更新式の中で最も重要なのは, \theta_fの更新式でしょう. \theta_fは, prediction lossを小さくするような勾配とdomain lossを大きくするような勾配によって更新されていることがわかります. これにより, G_fはLabelの予測には役立つ ( \hat{R}_S (h, f_S)を小さくする) がDomainの予測には役立たない ( d_{\mathcal{H}}を大きくする) ような入力表現 \mathcal{R}を得るための写像に近づいていくことが期待されます.

DANNのarchitectureは次の通りです. Domain Adversarialの名は, prediction layer G_yとdomain layer G_dが敵対的な関係にあることに由来すると思われます.

f:id:usaito:20190413045050p:plain
[Ganin+ 2015]のFigure 1

簡易実験

DANNのイメージをより鮮明に持つため, 人工データを用いた簡易実験を行ってみます. 本節は大いにこちらのrepositoryを参考にしました.

まず, scikit learnのmake_blobsを用いて人工データを生成します. sはsource, tはtargetを表しています.

Xs, ys = make_blobs(500, centers=[[0, 0], [0, 1]], cluster_std=0.2)
Xt, yt = make_blobs(500, centers=[[1, -1], [1, 0]], cluster_std=0.2)

描画すると次のような感じです.

f:id:usaito:20190413060233p:plain
入力の初期分布

このうち学習時にsource domainの入力とラベル, target domainの入力のみを用いて, テストデータにおけるラベルを精度よく予測したいというのがUDAの目標でした. いよいよDANNを学習します. 学習とテストは8:2で分け, OptimizerはMomentum (learning_rate=0.01, momentum=0.6), batchサイズは32, epoch数は5,000としました.

結果は次の通りです. 表の結果は, テストデータにおける最終epochの結果です. ちゃんとvalidationを用意して検証すればもう少し良い結果が出ると思いますが, target domainのラベルを学習時に全く用いていないのにも関わらず, 90%以上の精度を達成しています.

Source Target Domain
Cross Entropy 0.03359 0.20099 0.68040
Accuracy (%) 99.019 90.815 55.371

f:id:usaito:20190413060745p:plain
学習の様子

一方で, domainの判別はうまくいっていないことから, representation layerでdomainの判別が付かないような（ d_{\mathcal{H}}が小さいような）入力表現を得ることができていそうです. 実際, representation layerでの表現を抜き出してPCAで2次元に圧縮して描画してみると次のようになりました.

f:id:usaito:20190413060320p:plain
representation layerにおける特徴表現

これを見ると, source domain (赤, 薄赤) と target domain (青, 緑) が上と下に分かれていそうですが, 初期表現と比べるとかなり判別しにくくなっていることがわかります. 一方で, class 0とclass 1は綺麗に左右に分かれており, source domain target domainに関わらず, labelの判別はうまくいきそうなことがわかります. もちろんかなりシンプルな人工データを使ったからうまくいっているのですが, DANNのイメージが湧きやすい結果が出たのではないでしょうか.

さいごに

今回は, [Ganin+ 2015]で提案されたDomain Adversarial Neural Networksのarchitectureを理論背景も含めて整理し, 人工データを用いた追試を行ってみました. 個人的なモチベーションとしてはDANNそのものではなく, その別の分野への応用です. その話題についても今後触れようと思います.

参考

[Ben David+ 2007] Ben-David, S.; Blitzer, J.; Crammer, K.; and Pereira, F. 2007. Analysis of representations for domain adaptation. In NIPS, 137–144.
[Ben David+ 2010] Ben-David, S.; Blitzer, J.; Crammer, K.; Kulesza, A.; Pereira, F.; and Vaughan, J. W. 2010. A theory of learning from different domains. Machine Learning 79(1-2):151– 175.
[Ganin+ 2015] Ganin, Y.; Ustinova, E.; Ajakan, H.; Germain, P.; Larochelle, H.; Laviolette, F.; Marchand, M.; and Lempitsky, V. 2016. Domain-adversarial training of neural networks. Journal of Machine Learning Research 17(1):2096–2030.
[Kota Matsui 2019] Recent Advances on Transfer Learning and Related Topics. (https://www.slideshare.net/KotaMatsui/recent-advances-on-transfer-learning-and-related-topics)

EconMLパッケージの紹介 (meta-learners編)

はじめに

近年計量経済学と機械学習の融合分野の研究が盛り上がりを見せています. 例えば, KDD2018やNeurIPS2018で関連のTutorialが開催されるなどしています. その流れの一つとしてMicrosoft ResearchがEconMLというパッケージを公開していて非常に有用だと思ったので簡単に紹介します.

目次

  • Conditional Average Treatment Effects Estimation
  • EconMLとは
  • Meta-Learners
  • 用法と簡易実験
  • さいごに
  • 参考

Conditional Average Treatment Effects Estimation

ある特徴量で条件付けた際の介入の因果効果の期待値を Conditional Average Treatment Effects Estimation (CATE)と呼び, 次のように表されます.

\begin{aligned} \tau(X) = \mathbb{E} \left[ Y^{(1)} - Y^{(0)} \, | \, X \right] \end{aligned}

ここで, Y^{(1)}, Y^{(0)} はpotential outcomesです. 怪しい方は, 自分が以前書いた記事等をご参照いただけたらと思います.

CATEを推定することができれば, 嬉しいことがたくさんあります. 例えば, 因果効果がプラスであるような特徴量を持つ人だけに広告を打つことで商品の購入確率を最大化したり, 投薬計画を最適化することで生存率を改善できるかもしれません.

似たような目的を持つ分野にUplift Modelingと呼ばれるものがあります（参考1, 参考2）が, Uplift ModelingはA/Bテスト (RCT)によって収集された学習データがあることを前提とします. しかし多くの場合, A/Bテストを走らせて学習データを集めるようなことはコストの面から望ましくなく, 容易に実適用可能な技術とは言えないでしょう.

今回紹介するEconMLパッケージはより安価に手に入るobservational data (過去の介入方策が観測されている特徴量に依存しているようなデータ)から CATEを推定するための手法が豊富に実装されており, 低コストに個別化された介入施策を導くための非常に有用なツールになりうると思います.

EconMLとは

このパッケージについてはTwitterでも共有させていただきました.

github repositoryはこちら, documentはこちら にあります. documentにはCATE推定手法の概要も掲載されており, 分野を概観するのにも役立つでしょう. 今後blogでも何回か扱っていけたらと思っています.

Meta-Learners

今回はEconMLに収録されているCATE推定方法の中でも最もシンプルなMeta-Learnearsモジュールを使ってみます. このモジュールは, 既存の機械学習アルゴリズムを内部で用いるMeta的なCATEの推定方法を提供します. 収録されているMeta-Learnersは5つです. 以下簡単に説明を付しますが, documentを読んだ方が早い可能性があります. また, 機械学習アルゴリズム Mで Yを Xに回帰することで構築された予測器を M \left( Y \sim X \right)と記します.

T-Learner

Potential Outcomes ( Y^{(1)}, Y^{(0)})をそれぞれ個別にモデリングします.

\begin{aligned} \hat{\tau}(x) & = \hat{\mu}_1 (x) - \hat{\mu}_0(x) \end{aligned}

ここで, \hat{\mu_0} = M_1 \left( Y^{(0)} \sim X^{(0)} \right), \hat{\mu}_1 = M_2 \left( Y^{(1)} \sim X^{(1)} \right) です. また, X^{(0)}, X^{(1)}はそれぞれ, Xを T=0, T=1で条件付けた時の分布と一致する確率変数とし, 傾向スコアを e(X) = \mathbb{P} \left( T = 1 \, | \, X \right)としました.

S-Learner

特徴量に介入有無を表す変数 Tを含めます.

\begin{aligned} \hat{\tau}(x) & = \hat{\mu} (x, 1) - \hat{\mu}(x, 0) \end{aligned}

ここで, \hat{\mu} = M \left( Y^{obs} \sim \left( X, T \right) \right)です.

X-Learner

X-Learnerは少し複雑な手順を踏みます. [Kunzel+ 2017]で提案されているようです (まだ, ちゃんと読んでません).

\begin{aligned} \hat{D}^{(0)} & = \hat{ \mu }_1 (X^{(0)}) - Y^{(0)}, \quad \hat{D}^{(1)} = Y^{(1)} - \hat{ \mu }_0 (X^{(1)}) \\ \hat{\tau}_0 & = M_3 \left( \hat{D}^{(0)} \sim X^{(0)} \right) \\ \hat{\tau}_1 & = M_4 \left( \hat{D}^{(1)} \sim X^{(1)} \right) \\ \hat{\tau}(x) & = e(x) \hat{\tau}_0 (x) + (1 - e(x)) \hat{\tau}_1(x) \end{aligned}

T-Learnerの時と同様に, \hat{\mu_0} = M_1 \left( Y^{(0)} \sim X^{(0)} \right), \hat{\mu}_1 = M_2 \left( Y^{(1)} \sim X^{(1)} \right) です.

DA-Learner (Domain Adaptation Learner)

DA-Learnerは, X-Learnerにおける \hat{\mu_0}, \hat{\mu}_1の学習に共変量シフトを用いた手法です. 具体的には,

\begin{aligned} \hat{\mu_0} & = M_1 \left( Y^{(0)} \sim X^{(0)}, weights = \frac{e \left(X^{(0)} \right)}{1 - e \left( X^{(0)} \right) } \right) \\ \hat{\mu}_1 & = M_2 \left( Y^{(1)} \sim X^{(1)} , weights = \frac{1 - e \left(X^{(1)} \right)}{e \left( X^{(1)} \right) } \right) \\ \hat{D}^{(0)} & = \hat{ \mu }_1 (X^{(0)}) - Y^{(0)}, \quad \hat{D}^{(1)} = Y^{(1)} - \hat{ \mu }_0 (X^{(1)}) \\ \\ \hat{\tau} & = M_3 \left( \left( \hat{D}^{(0)}, \hat{D}^{(1)} \right) \sim \left( \hat{X}^{(0)}, \hat{X}^{(1)} \right) \right) \end{aligned}

\hat{\mu_0}, \hat{\mu}_1の学習の際に傾向スコアに依存する重みを考慮に入れた損失を用いています. この重みは, importance weightと呼ばれることもあるやつです.

DR-Learner (Doubly Robust Learner)

DR-Learnerは, Doubly Robustを用いてCATEを代替するようなsurrogate outcomeを作り, それを Xに回帰する方法です. Y_{i,t}^{DR}を次のように定義します.

\begin{aligned} Y_{i,t}^{DR} & = \hat{E} \left[ Y \, | \, T = t, x_i \right] - \mathbb{I} \{ T_i = t \} \cdot \frac{ Y_i^{obs} - \hat{E} \left[ Y \, | \, T = t, x_i \right] }{e (x) } \\ \hat{\tau} & = M \left( \left( Y_{*, 1}^{DR} - Y_{*,0}^{DR} \right) \sim \left( \hat{X}^{(0)}, \hat{X}^{(1)} \right) \right) \end{aligned}

hatが付いている部分は, 自分たちでモデリングする必要があります.

用法と簡易実験

EconMLパッケージは, scikit-learnと同じようなインターフェースで実装されており, とても扱いやすいです. 例えば,

# T-learnerを初期化.
controls_model = GradientBoostingRegressor(n_estimators=100, max_depth=6, min_samples_leaf=15)
treated_model = GradientBoostingRegressor(n_estimators=100, max_depth=6, min_samples_leaf=15)
T_learner = TLearner(controls_model, treated_model)
# 学習 (fit methodの引数に介入有無を表す変数Tがあるのがsklearnとの違い)
T_learner.fit(y_train, t_train, X_train)
# CATEを推定.
cate_pred = T_learner.effect(X_test)

EconMLに収録されている5つのMeta-Learnersの予測性能を, [Powers+ 2017]で用いられている8つの人工データセットを用いて評価してみます． ただし今回はこのパッケージを触ってみたという意味合いが強いので, ハイパーパラメータチューニングはしておらず, 正確な性能評価ではないのでご注意ください. 非常に簡単に使えるので, 実際に触ってみるのが良いと思います.

結果は次の通りです. 人工データを使ったため, 真のCATEに対するRMSEで評価しました. 学習データ, テストデータは数はともに1,000で, 同様の実験を15回繰り返し箱ひげ図をplotしました.

f:id:usaito:20190407205402p:plain

f:id:usaito:20190407210453p:plain

X-Learnerが良さそうな雰囲気が漂っているので, 論文 ([Kunzel+ 2017])を調査してみる必要がありそうです. また, 今回はoutcomeのモデリングにはGradient Boosting Regressorを, 傾向スコアのモデリングにはRandom Forest Regressorを用いました. scenarioごとにこれらのモデリング方法を変えた方が良いと思われるので, 実際はちゃんとチューニングしましょう.

さいごに

EconMLパッケージを簡単に触ってみました. まだ, 実装されている手法やアルゴリズムが提案された論文を読みきれていないのでちゃんと追っていこうと思います. 今回はMeta-Learnersを扱いましたが, 今後より複雑なアルゴリズムを理解しながら触っていくつもりです.

参考