【書評/書籍紹介】施策デザインのための機械学習入門 〜データ分析技術のビジネス活用における正しい考え方

本記事では、2021年8月4日に技術評論社さまから発売された施策デザインのための機械学習入門〜データ分析技術のビジネス活用における正しい考え方（齋藤優太、安井翔太 著、株式会社ホクソエム 監修; 以下本書と呼びます）の紹介を行います。

gihyo.jp

なお、本書の詳しい内容や執筆背景については著者がすでにブログを書かれているので、まずは著者の記事をご覧になった上で、一読者の感想として私の書評記事を読んでいただければと思います。

usaito.hatenablog.com

紹介の経緯

以前からお世話になっていた著者の齋藤(@usait0)さんからありがたいことに電子版の献本をしていただき、お盆休みの数日を使い読ませていただきました。

献本への感謝の気持ちに加え、本書に内容を十分に理解した上で一緒にデータ分析の施策デザインの質を高めていく仲間が増えると私自身もハッピーになりそうな予感を強く感じたため、一人でもそのような仲間を増やせるとよいなという気持ちで記事を書くことにしました。本書が指南書・啓蒙書という形をとっていることにあやかり、本記事もそれに乗っかった形をとろうと思います。

自分で本記事を読み返してみると、「全体的にベタ褒めしすぎでは？」「記事の書き方が冗長すぎでは？」という印象を持たれそうな気がしましたが、心の赴くままに書いた結果なのでご容赦いただければと思います。

本書の目次

目次は以下のようになっています:

1章 機械学習実践のためのフレームワーク
2章 機械学習実践のための基礎技術
3章 Explicit Feedbackを用いた推薦システム構築の実践
4章 Implicit Feedbackを用いた推薦システムの構築
5章 因果効果を考慮したランキングシステムの構築
付録A 演習問題

詳細はHPで確認していただければと思いますが、個人的には、1章と2章が良い意味で期待を裏切られる内容でした。

章の粒度でタイトルだけを見ると「1-2章では 機械学習とは〜〜 のようなふわっとした話が展開されるのかな？」という気がしなくもないですが、非常にしっかりとした内容で、むしろこの2つの章が本書の根幹になる部分だと感じました（IPWLearnerでOff Policy Learningを行うという、知名度は低いが本書の根幹となる重要な手法もここで紹介されています）。

それに伴い、本書の読み方としては、1章から順に読んでいき、後半の具体例を眺めることで、自分の普段向き合っているビジネス課題との対応を考えながら進めていくというやり方が適切そうだと感じました。

対象読者

はじめに の章によると、本書のターゲットは

機械学習を活用したビジネス施策の実践に取り組んでいるすべての機械学習エンジニアやデータサイエンティスト

だと記載されています。ここで最も重要なのは すべての という部分で、特定の機械学習のタスク・技術にフォーカスした内容ではなく、ほとんどすべての機械学習施策で共通して役に立つ内容を扱っている点だと考えます。読者が解いているタスクがユーザのデモグラ予測だろうがアイテムの推薦や並び替えであろうが、またそのために使う手法が勾配ブースティングだろうが深層学習であろうが（なんならルールベースであろうが）、本書のフレームワークや思考方法を活用することで施策の質を高めてくれる可能性があるでしょう。

その理由は、本書が「機械学習を機能させるために、どのようなステップが必要であるか」を言語化した上で、現実に起こりうる様々なビジネス課題に対して、「どのように解くか」ではなく「何を解くか」を上流部分から定式化していくスタンスをとっていることにあると考えます。これについての詳細は後で説明します。

私個人の感想としては、すべての という表現は誇張ではないものの、読者の背景に合わせて訴求ポイントは変わりそうだと感じました。具体的には以下の3つのパターンを想定し、それぞれのパターンに対する訴求ポイントを、本書の内容について触れながら説明していこうと思います:

1. 「目次に出てくる 機械学習実践に潜む落とし穴 や バイアス や 反実仮想機械学習 に興味があり、ちゃんと勉強してみたい」と思っている人
2. 「データ分析において 何を解くか を考えるのが元々得意だし、反実仮想機械学習 のことも知ってる。でも結局この話はケースバイケースなので、わざわざ本書を買う必要はなさそう」と感じている人
3. 「普段機械学習をやっているが、詳細目次に出てくる 反実仮想機械学習 といった専門用語がピンとこないので買うか迷っている」人

訴求ポイント:

1. ご興味にぴったりの本だと思うので、マストバイだと思います。
2. 私の知る限りですが、 何を解くか を考えて定式化するという文脈において、本書以上に整理された資料や議論は存在しないです。ある程度自信のある人であっても、本書を活用することでより思考の質を高められると考えます（私自身、著者が公開されているスライドや論文の大部分には元々目を通していましたが、それでもなお多くの学びがありました）。
3. 本書の本質は 反実仮想機械学習 という方法論ではなく、 何を解くか を定式化する部分にあると思うので、その部分に興味がある場合はおすすめできます。また、 反実仮想機械学習 自体がここ5年くらいのトレンドなので、分野として追い始めるきっかけに使うという手もありそうです。

ちなみに、2021年8月13日時点でAmazonの情報学・情報科学全般関連書籍カテゴリーの中でベストセラー1位だったので、すでにかなり多くの方が購入されていると推測され、この記事の存在価値は書き始める前からすでに危ういです（書評記事に仲間集めの方向性を強めに置いたのはそれも理由の一つです）。

対象でない読者

個人の感想としては、「データ分析に関わるビジネス課題の 定式化 に興味がない人」は対象外かなと感じました。

特定の機械学習手法の研究者でそのクオリティを高めていきたい方や、あるいは分析プロジェクトに企画職として関わる方の中には、これに該当する方がいらっしゃるかと思います。

判断が難しいのは、本書の概要だけを理解して 何を解くか について専門家と議論をしていきたい企画職の方でしょうか。本書では、想定知識としては統計学や機械学習の基礎知識を仮定されています。機械学習の 損失関数の最小化 を中心としたアプローチに慣れていない方が本書に対してどのような印象を持つかはわからないですが、そこは我々専門家がコミュニケーションしながら個別に期待値調整をする必要があるかな、というのが私の感想です。

本の中身

本書を読んだ中で個人的に重要だと感じた点を列挙します:

1. ビジネス課題に対して「何を解くか」を、上流部分から定式化しているところ
2. ビジネス課題を解く際に有用な思考フレームワークを定義し、そのフレームワークに沿った議論が一貫してなされているところ
3. 一つの定式化や手法を正解として与えるわけではなく、試行錯誤を通してより良い定式化に辿り着くための実務家向け指南書としての側面が強調されているところ
4. 反実仮想機械学習（CFML）の手法を丁寧に解説しているところ
5. データやコードがすべて公開されており、実際に使ってみるハードルが低いところ

4と5 は議論の余地が少なそうなので、それ以外について以下で詳細に補足します。ただし、ここから先は未読者への紹介というよりは、既読者とのコミュニケーションと自己満足の側面が強くなっていきます。

定式化

対象読者のところでもお話しましたが、本書は「機械学習を機能させるために、どのようなステップが必要であるか」を言語化した上で、現実に起こりうる様々なビジネス課題に対して、「どのように解くか」ではなく「何を解くか」を上流部分から定式化していくという点で、他の本にない独自の価値を持っています。

ここで注目すべきは 定式化 という部分です。

世の中には「モデルの性能を追い求める前に、どういう問題を解くかを設計することが重要だ」ということを語る方は多く存在していますが、「具体的にどのような思考のフレームを使えばそのクオリティを高めていけるのか」という問いに対して明確な回答を持っている方を私はこれまで観測したことがありませんでした。私自身もそうだったのですが、このような議論は結局ケースバイケースになってしまい、「具体と抽象の間をうまく行き来しながらちょうど良い解像度で何を解くかを語ることは非常に難しいよね」というある種の諦めを持っている方は多いのではないでしょうか。

本書はこれに対し、実例をベースに「悪い問いの立て方」と「良い問いの立て方」の違いを数式の世界で表現しています。1章を読んだ時点で、定式化 についてのこれまでの認識を一歩前に進める世界観を提示しているような感想を抱きました。

具体例を挙げすぎるとネタバレになってしまうので難しいのですが、Xという属性の人の値Yを関数fで予測して損失関数lで評価したくなるであろうケースにおいて、

]

を解くべき問題としていいんだっけ？それって観測データで素朴に経験近似していいんだっけ？のような話を突き詰めていく感覚です。

フレームワーク

本書では、ビジネス課題を解くにあたって独自の思考フレームワークを導入し、そのフレームワークに基づいて一貫して議論が展開されています。以下の目次からなる 1.2 機械学習実践のためのフレームワーク で提案されている6段階からなるフレームワークです:

1.2.1 KPIを設定する
1.2.2 データの観測構造をモデル化する
1.2.3 解くべき問題を特定する
1.2.4 観測データを用いて解くべき問題を近似する
1.2.5 機械学習モデルを学習する
1.2.6 施策を導入する

データ分析のフレームワークというと、多くの方が第一に思い浮かべるのは CRISP-DM だと思います。私自身CRISP-DMのフレームワークに従って問題を考えることは多いですし、CRISP-DM自体の有用さは間違いないと思います。しかし、CRISP-DMだけでは、より良い定式化をするためのヒントを得ることは難しいと感じていました。参考のため、本書とCRISP-DMの対応を考えてみます。

1.2.1 KPIを設定する は、Business Understandingに対応しそうです。

1.2.2 データの観測構造をモデル化する は、Data Understanding に対応しそうです。ただし、CRISP-DMではこのステップでは可視化やEDAがメインになることが多いですが、本書ではデータ生成過程を式で表現するところがメインになります。CRISP-DMの Data Understanding では「生成過程を式で書くこと」を明言されていなかったという認識をしており、実際の業務の場でも、生成過程をしっかりと式で書き切っているケースばかりではないでしょう。本書はその部分を明示的に書いている、という点が重要な違いの一つだと感じました。

1.2.3 解くべき問題を特定する・ 1.2.4 観測データを用いて解くべき問題を近似する・1.2.5 機械学習モデルを学習する は Modeling・Evaluation に対応します。この部分の具体性の違いも重要な部分だと感じました。CRISP-DMでは、各ステップについてそこまで踏み込んだ議論はされていない認識です。一方本書では、「解くべき問題を間違えるとどうなるのか」について例示し、このステップの解像度を高めています。また、「解くべき問題の特定と観測データを使った近似の間にギャップが存在しうる」という点に注目し、その二つのステップを明確に分離しています。モデルの学習や評価においても、学習したモデルの評価を事前に行うOff Policy Evaluationという仕組みを詳しく説明しているという点で、CRISP-DMよりは厳密な評価を想定していると言えるでしょう。

最後の 1.2.6 施策を導入する は Deployment に対応しています。本書のユニークな特徴として、「施策の導入によって将来のデータが生成される」という点に着目し、MLOpsの文脈では語られることがほとんどない「将来のモデル改善のためのデータ生成装置としての機械学習モデルの導入」を実践しているところが挙げられるでしょう。

一方、CRISP-DMの Data Preparation に当たる内容は本書では明確に述べられていません。 1.2.5 機械学習モデルを学習する の中にデータを準備する工程も含まれていそうですが、前処理や特徴量設計については詳しい書籍が既に多く存在しているため、本書はその詳細については触れないというスタンスをとっているようでした。本書でも繰り返し述べられていますが、特定のフレームワークを絶対視することなく、目的に応じて適宜チューニングしていく姿勢が重要になります。いずれにせよ、思考を整理してくれる良質なフレームワークが増えるに越したことはないので、分析者にとって助けになる内容だと感じました。

試行錯誤・実務家向け指南書

本書の大きな特徴として、特定のフレームワークや手法を絶対視せずに、目的に応じてチューニングしていくことを何度も強調されている点が挙げられます。その上で、実務家が自分で応用できる力を養成すべく、「素朴に定式化したらどうなるか」「それでどういう問題が発生するか」「改善するためにはどうすればいいか」という流れで説明がなされています。

冒頭で述べた通り、何を解くか というテーマの議論は非常に難しくケースバイケースであるので、この点を強調されているのは非常にありがたいです。実際私も、本書を読んだことをきっかけに普段の定式化行動を省みてみましたが、反省すべき点や改善すべき点が多々存在していました。また、本書の内容をさらに進めて「こういうことも考えられそうだな」という案もいくつか出てきました。

全部の案は書ききれないので、参考までに私の思考の一部を以下で掲載しておきます。

2.2のcolumn

ログデータを収集した際に稼働していた意思決定モデルが行動選択に用いた特徴量を分析者がすべて把握していること という仮定が紹介されており、ビジネス現場ではその仮定は満たされることが多いと述べられていました。私の解釈では、この仮定はoutcome regression modelの交絡を防ぐためだと考えています（IPS推定量については処置確率だけわかっていればよい）。そして、outcome regression modelを作る場合には、過去施策の特徴量を把握しているだけではなく、過去施策の特徴量をすべて使って交絡を調整することまでが必要になるように思えました。これが何を意味しているかというと、「outcome regression modelを使ってモデル改善をする前提では、モデルがupdateされるたびに考慮するべき特徴量が増えていく」ということであり、「無闇矢鱈にモデルに使う特徴量を増やすと後々のモデル改善が大変になりかねない」というメッセージに繋がりそうだなと感じました。

Off Policy Evaluationとモデル選択

本書ではOff Policy Evaluation (OPE) によってモデルの性能を評価し、より良いモデルを選択するという手続きが説明されています。これについて、本書の思考に則ると、「より良いモデルを選択する ことがやりたいのであれば、 モデルの性能の評価 は解きたい問題と違うタスクを解いていることになるのではないか」という疑問を抱くのが自然かと思います。

この点について本書の著者に質問してみたところ、 Confident Off-Policy Evaluation and Selection through Self-Normalized Importance Weighting という論文がAISTATS 2021で発表されており、そのような疑問を解決していく方向性も存在しているようでした。自分がゼロから理論を組み立てるのはなかなか大変そうですが、参考にできるところを取り入れながら業務に応用していきたいところです。

regression-EM周り

regression-EMについて、著者が以前以下の疑問を書かれていました:

regression-EMでやっているrelevanceパラメータをGBDTで予測するパートだけど, それがそもそもUnbiased LTRでやりたいことなので, regression-EMのアルゴリズムで精度よくpropensity推定できるならそもそもpropensity推定のこと考える必要ないのでは？という本末転倒な感じがした...

本書にはこの疑問点についての言及は存在していなかったのですが、私自身regression-EMの論文を読んだときに同じような感想を持っていました。しかし、本書の思考をトレースしてみると、この2段階のアプローチをとっている理由は、最終的に解くべき問題に対して改めてチューニングするためなのかなという仮説に至りました。多少粗くてもいいいのでポジションバイアスを推定してくれれば、あとは後段の学習でなんとかするという気持ちがあるのかもしれないです。これについては、本書で紹介されているポジションバイアス推定の後続研究を読んでいくと何かしらのヒントが得られるかもしれないです。

おわりに

以上で紹介は終わりです。

改めてではありますが、献本していただいた齋藤(@usait0)さん、ならびに執筆に関わられたみなさま、ありがとうございました。

同じような問題に興味を持たれた方と、本書をもとに議論やお仕事をする機会があれば嬉しいなと思っております。

Marginal Structural Model のパラメータ推定に重み付き最小二乗法を用いることの正当性

本記事では、 Marginal Structural Model (MSM) のパラメータ推定に重み付き最小二乗法 (IP-Weighted Least Squared method; IP-WLS method) を用いることの正当性についての解釈を試みます。what if book の12章で紹介されている話の一部を深く掘り下げていきます。

記事のゴール

以下の3つの問いに答えることを目指します（自分が12章を読んだときに感じた疑問でした）:

1. 処置が2値の場合には、IP-WLS の推定量が Hajek estimator と等しいという主張が what if book P152 で書かれているが、これは本当か？
2. 処置が多値の場合にも、IP-WLS の推定量は Hajek estimator と等しいのか？
3. 処置が連続の場合にも、何かしら似たような解釈はできるのか？

執筆終了時点での私の回答としては以下になります:

1. 本当であるということを証明します。
2. 星野先生の本の議論を参考にしつつ、等しいことを証明します。
3. モーメント法・M推定量の文脈から、自分なりに正当化を試みます。

なお、本記事のタイトルの「正当性」という単語には、

• 特定（処置が離散）の場合に、既知の推定量と同じ結果が得られるという意味で手法の正当化ができる
• 処置が連続の場合も、離散の場合の拡張として捉えると、一定の妥当性があるように思える

といった気持ちを込めています。

新規性

モーメント法・M推定量の文脈で処置が連続の場合のお気持ちを書いてある日本語記事は見つけられませんでしたので、それを新規性とさせていただきます（正確性は保証できないです）。

前提知識

• 過去の自分の記事 における、 Hajek estimator のモーメント法による定式化の雰囲気がわかっていること

what if book のMSMのおさらい

what if book 12章については過去の記事でも紹介しましたが、MSMやそのパラメータ推定方法については述べていなかったので、軽くおさらいします。

fullflu.hatenablog.com

モデリング

MSMとは、potential outcome の期待値を（周辺化して）直接予測するようなモデルを立てることを意味します。

過去の記事でも紹介した exchangeability などの仮定のもとで、そのモデルのパラメータの推定結果を用いて因果効果を良い感じに*1推定することが可能です。

処置が2値のとき、以下のような MSM のモデルが 12.4 節で紹介されています:

\begin{eqnarray} \mathbb{E}[Y^{a}] = \beta_0 + \beta_1 a \end{eqnarray}

個人的には を使って以下のように書いた方が綺麗だと思うのですが、必要に迫られるまでは慣習にしたがって を使わずに書きます。

\begin{eqnarray} \mathbb{E}[Y^{a=a'}] = \beta_0 + \beta_1 a' \end{eqnarray}

このとき、以下のように因果効果を計算できます:

\begin{eqnarray} \mathbb{E}[Y^{a=1}] &=& \beta_0 + \beta_1\\ \mathbb{E}[Y^{a=0}] &=& \beta_0\\ ATE &=& \beta_1 \end{eqnarray}

また、処理が連続のとき、以下のようなモデルが紹介されています:

\begin{eqnarray} \mathbb{E}[Y^{a}] = \beta_0 + \beta_1 a + \beta_2 a^2 \end{eqnarray}

あくまでこれは2次式までを考慮した一例であり、実際にはより複雑なモデルを立てることもありえるということに注意してください。

パラメータ推定

MSM のパラメータ推定にあたって、 IPW によって重み付けした最小二乗法、通称 IP-WLS を用いる手法が 12.4 節で紹介されています。

実装は、Python の場合は statsmodels、R の場合は lm や gee を使うことになるでしょう。

パラメータの分散を推定する際には、厳密に解析するか、 bootstrap を使うか、サンドウィッチ分散を使うかの三択になるとも紹介されています。

また、P152 の左側には、 IP-WLS の結果が Hajek estimator と一致すると書かれています。

処置が2値の場合、本当に IP-WLS は Hajek estimator と一致するのか？

まず、一つ目の疑問に答えていきます。愚直に解析します。

個人 が処置 を受ける確率 ] を と表記すると、 IP-WLS は

\begin{eqnarray} L = \sum_{i=1}^{n} \Biggl[ \frac{A_i}{\pi_i} \bigl(Y_i - \beta_0 - \beta_1 A_i \bigr)^2 + \frac{1-A_i}{1-\pi_i} \bigl(Y_i - \beta_0 - \beta_1 A_i \bigr)^2 \Biggr] \end{eqnarray}

を最小化する問題を考えればよいことになります。偏微分を行った式

\begin{eqnarray} \frac{\partial{L}}{\partial{\beta_1}} &=& 2\sum_{i=1}^{n} \Biggl[ \frac{A_i}{\pi_i}\bigl(Y_i - \beta_0 - \beta_1 \bigr)(-A_i) + \frac{1-A_i}{1-\pi_i} \bigl(Y_i - \beta_0 - \beta_1 A_i \bigr)(-A_i) \Biggr] \\ &=& 2\sum_{i=1}^{n} \Biggl[ \frac{-A_i}{\pi_i}\bigl(Y_i - \beta_0 - \beta_1 \bigr) \Biggr] \\ &=& 0 \end{eqnarray}

を整理すると、] の推定量は

\begin{eqnarray} \sum_{i=1}^{n} \Biggl[ \frac{A_i}{\pi_i}\bigl(Y_i - \hat{\theta_1} \bigr) \Biggr] &=& 0 \end{eqnarray}

の解となります。

これを踏まえてさらに、

\begin{eqnarray} \frac{\partial{L}}{\partial{\beta_0}} &=& 2\sum_{i=1}^{n} \Biggl[ \frac{-A_i}{\pi_i}\bigl(Y_i - \beta_0 - \beta_1 \bigr) - \frac{1-A_i}{1-\pi_i} \bigl(Y_i - \beta_0 - \beta_1 A_i \bigr) \Biggr] \\ &=& -2\sum_{i=1}^{n} \Biggl[ \frac{1-A_i}{1-\pi_i} \bigl(Y_i - \beta_0 \bigr) \Biggr] \\ &=& 0 \end{eqnarray}

を整理すると、] の推定量は

\begin{eqnarray} \sum_{i=1}^{n} \Biggl[ \frac{1-A_i}{1-\pi_i} \bigl(Y_i - \hat{\theta_0} \bigr) \Biggr] &=& 0 \end{eqnarray}

の解になります。

これらの式を、前提知識の章のリンクで説明した Hajek estimator のモーメント法による定式化と比較すると、全く同じ形になっているということが確認できると思います。

したがって、 処置が2値のとき、 IPW で重み付けをした最小二乗法の結果は、確かに Hajek estimator と等しいことが示せました。

処置が多値のときはどうなるか？

what if book には処置が多値の場合については説明がされていないので、2値の場合をよしなに拡張して定式化を行います。

ある処置のレベル に対して、周辺構造モデルを以下のように表現します:

\begin{eqnarray} \mathbb{E}[Y^{a=j}] = \beta_j \end{eqnarray}

個人 が処置 を受ける確率 ] を と表記すると、 IP-WLS では

\begin{eqnarray} L = \sum_{i=1}^{n} \sum_{j=1}^{J} \frac{\delta(A_i-j)}{\pi_{ij}} \bigl(Y_i - \beta_j \bigr)^2 \tag{1} \end{eqnarray}

を最小化する問題を考えればよいことになります。 はディラックのデルタ関数です。

\begin{eqnarray} \frac{\partial{L}}{\partial{\beta_j}} &=& -2\sum_{i=1}^{n} \frac{\delta(A_i-j)}{\pi_{ij}} \bigl(Y_i - \beta_j \bigr)\\ &=& 0 \end{eqnarray}

を整理すると、] の推定量は

\begin{eqnarray} \sum_{i=1}^{n} \frac{\delta(A_i-j)}{\pi_{ij}} \bigl(Y_i - \hat{\theta_j} \bigr) &=& 0 \end{eqnarray}

の解になります。

多値の場合の Hajek estimator の定義は what if book には載っていなかったはずですが、この推定量は、2値の場合の Hajek estimator を素朴に拡張した式になっていることが確認できるかと思います。

M推定量の文脈による正当化（処置が離散の場合）

これまで、Hajek estimator との対応によって MSM の正当化を試みてきましたが、ここではM推定量の立場から改めて整理することを試みます。

以下の星野本の P81-82 を参考にします。

調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

まず、天下り式に、処置 に対する potential outcome とパラメータを引数にとった非負のスコア関数 というものを導入します。

二乗誤差の最小化を考える場合、以下のようなスコア関数を考えることになるでしょう:

\begin{eqnarray} m_j(Y^{a=j}, \theta_j) &:=& -(Y^{a=j} - \theta_j)^2 \end{eqnarray}

さらに、以下のようなスコア関数の和の期待値を最大化したい、と考えたとします:

\begin{eqnarray} \mathbb{E}_{\tilde{Y}} \Biggl[ \sum_{j=1}^{J} m_j(Y^{a=j}, \theta_j) \Biggr] \tag{2} \end{eqnarray}

ここでは、期待値は potential outcomes についてとっています。

この式からは、すべての処置の水準の重要度を等しく扱った上で、処置ごとの potential outcome をできる限り正しく推定したい という気持ちを感じることができますし、この気持ち自体をある意味で「自然」なものだとして正当化することは可能だと考えられます（私個人の感覚ですが）。

しかし残念なことに、式 (2) には counterfactual outcome が含まれているので、そのまま経験近似を行うことが不可能です。

そこで、式 (2) と期待値の意味で等しくて（不偏で）、かつ経験近似が容易に行える ような別の式を考えたくなります。

そのような性質を持つ式の一つが以下のような形のもの（期待値の中身）です。不偏性は conditional exchangeability の仮定によって示せます。

\begin{eqnarray} \mathbb{E}_{\tilde{Y}, L, A} \Biggl[ \sum_{j=1}^{J} \frac{\delta(A-j)}{\pi_{j}} m_j(Y^{a=j}, \theta_j) \Biggr] \tag{3} \end{eqnarray}

そして、consistency の仮定を用いてこの式 (3) を経験近似したものが、式 (1) のIP-WLS であると解釈できます。

また、式 (3) をパラメータで微分してから、微分と積分が交換できることを仮定して経験近似をすると、M推定量を考えていることになるような気持ちになれます。

改めてまとめます。

• 「すべての処置の水準の重要度を等しく扱った上で、処置ごとの potential outcome をできる限り正しく推定したい」という気持ちを式 (2) として表現することを認めた上で
• 「式 (2) と期待値の意味で等しくて、かつ経験近似が容易に行える」ような便利な式を考えると
• 式 (3) の形は自然と導出され、推定の正しさを二乗誤差で測定することに同意すれば、その経験近似により IP-WLS の式が得られる（& M推定量として解釈できる）

というのがここでの主張になります。

一見突拍子もないように感じた IP-WLS が自然な導出により得られる、ということが共有できたかと思います。

M推定量の文脈による正当化（処置が連続の場合）

最後に、処置が連続の場合についても直感的な解釈を与えていきます。

ここから先の話はどの本にも載っていなかったので、自分の感覚を信じて論理を組み立てていきます。

まず、先ほどと同様に式 (2) のようなものを考えてみると、処置が連続の場合、処置について和ではなく積分を考える必要が出てくることがわかります。

また、すべての処置の水準の重要度を等しく扱おうとすると、処置の区間が有限でない場合積分値が無限に飛んでいってしまいます。

そこで、処置の密度関数 によって重み付けを行った上で、処置ごとの potential outcome をできる限り正しく推定したい という気持ちで式 (2) を修正してみます:

\begin{eqnarray} \mathbb{E}_{\tilde{Y}} \Biggl[ \int_{a' \in \mathcal{A}} f(a') m(Y^{a=a'}, \theta, a') da' \Biggr] \tag{4} \end{eqnarray}

なお、ここから先は potential outcome を扱うときに を使った表記に切り替えます。

また、スコア関数 の引数に処置の情報を含むようにします。

式 (4) を見てみると、当然 potential outcomes は観測できないので、この最大化をそのまま考えるのは困難であり、代理の式を探したくなります。

処置が離散の場合を参考にして、 をかけて期待値をとってみると、

\begin{eqnarray} \mathbb{E}_{\tilde{Y}, L, A} \Biggl[ \int_{a' \in \mathcal{A}} f(a') m(Y^{a=a'}, \theta, a') \frac{\delta(A - a')}{f(a' \mid L)} da' \Biggr] &=& \int_{a' \in \mathcal{A}} \mathbb{E}_{L}\Biggl[ \mathbb{E}_{A \mid L} \bigl[\frac{\delta(A - a')}{f(a' \mid L)} \bigr] \mathbb{E}_{Y^{a=a'}}[m(Y^{a=a'}, \theta, a')] \Biggr] f(a')da' \\ &=& \int_{a' \in \mathcal{A}} \Biggl[\mathbb{E}_{Y^{a=a'}}[m(Y^{a=a'}, \theta, a')] \Biggr] f(a')da' \tag{5} \end{eqnarray}

となり、式 (4) に対する不偏性を持つことが確認できます*2。

したがって、これを経験近似した以下の式 (6) を最大化するという操作には、処置が離散の場合と同様にM推定量の文脈からその正当性を与えられそうです。

\begin{eqnarray} \sum_{i=1}^{n} \Biggl[ \int_{a' \in \mathcal{A}} \biggl( \frac{\delta(A_i - a') f(a')} {f(a' \mid L_i)} m(Y_{i}^{a=a'}, \theta, a') \biggr) da' \Biggr] &=& \sum_{i=1}^{n} \Biggl[ \frac{f(A_i)} {f(A_i \mid L_i)} m(Y_{i}^{a=A_i}, \theta, A_i) \Biggr]\\ &=& \sum_{i=1}^{n} \Biggl[ \frac{f(A_i)} {f(A_i \mid L_i)} m(Y_{i}, \theta, A_i) \Biggr] \tag{6} \end{eqnarray}

個人的には、不偏性を示すときと経験近似をするときとで、デルタ関数の処理が異なってくるというのがエモポイントです。

そして、式 (6) のスコア関数に二乗誤差を用いると stabilized weights を用いた IP-WLS そのものになります。

以上の議論により、処置が連続の場合の IP-WLS についても、M推定量の文脈での正当化を行うことができ、その意味で「自然な」推定手法だと解釈できると考えられます。

おわりに

本記事では、 IP-WLS で MSM のパラメータ推定を行うことの正当性について、Hajek estimator や M推定量との対応をとりながら紹介を行いました。

少し記述が回りくどくなってしまいましたが、モヤモヤを持たれている方の助けになれば幸いです。

誤りを見つけられた方は、そっと教えていただけると幸いです。

*1:一致性を持ち、漸近分散が小さいという意味で良い感じです

*2:密度関数の notation が崩れてますが、この辺は what if book を見て察していただけると幸いです

IPW 推定量の漸近分散についての補足

過去の記事の中で、 IPW の漸近分散について、以下の主張を行っていました。

よって、 Potential outcome の符号が等しければ右辺が0以上になり、「Horvitz-Thompson estimator より Hajek estimator の方が漸近分散が小さい」という主張が成立すると考えられます。

しかし、 @BluesNoNo さんとやりとりをしていく中で、 の値に着目することで別の条件を考えられる可能性を教えていただきました。

ありがとうございます．自信は全くありませんが，単純計算で，
w1=E[1/π], w2 = E[1/(1-π)]としたら，
(13)式 = w1 θ1^2 + w0 θ0^2 - (θ1 - θ0)^2
= (w1 - 1) θ1^2 - 2 θ0 θ1 + (w2 - 1) θ2^2
となり，これが0以上になるのは，1 < (w1 - 1)(w2 - 1)
となる気もします．間違えていたら，すみません．

— Tarotan (@BluesNoNo) June 2, 2020

その後、さらに式展開を行っていく中で、特別な条件なしで漸近分散の比較ができそうだということがわかってきたので、本記事の中で補足を行いたいと思います。

記事のゴール

本記事では、以下の主張を行うことをゴールにします:

• Horvitz-Thompson estimator より Hajek estimator の方が漸近分散が大きくなることはない（傾向スコアモデルgivenのとき）

新規性

以下の2点で新規性があると考えています。

• 上記の主張を（少なくとも日本語で簡単に）示したこと
• その過程で共分散についての簡単な補題を示したこと

補題については車輪の再発明も含んでいそうですが、証明を探すのに苦労した（結局自分で示しました）ので、日本語で必要な情報をサクッとまとめているという点で、一応新規性としておきます*1。

これ以降は証明をするだけなので、証明に興味がある方のみご覧いただければと思います。

過去の記事のおさらい

過去の記事では、傾向スコアモデルが与えられたもとで、以下のような式展開により、 「Potential outcome の符号が等しければ、 Horvitz-Thompson estimator より Hajek estimator の方が漸近分散が小さい」という主張に至っていました。

\begin{eqnarray} Var_{HT}(\hat{\theta}_{1} - \hat{\theta}_{0}) - Var_{Hajek}(\hat{\theta}_{1} - \hat{\theta}_{0}) &=& \mathbb{E}\Bigl[ \frac{\theta_1^{2}}{\pi} \Bigr] + \mathbb{E}\Bigl[ \frac{\theta_0^{2}}{1-\pi} \Bigr] - (\theta_1 - \theta_0)^2 \\ &>& \theta_1^2 + \theta_0^2 - (\theta_1 - \theta_0)^2 \ \ \ (\because 0 < \pi < 1) \\ &=& 2\theta_1\theta_0 \tag{1} \end{eqnarray}

新たな主張

記事の主張を定理にすると、以下になります:

[Theorem 1]

\begin{eqnarray} Var_{HT}(\hat{\theta}_{1} - \hat{\theta}_{0}) - Var_{Hajek}(\hat{\theta}_{1} - \hat{\theta}_{0}) &\ge& 0 \tag{2} \end{eqnarray}

以下の流れで証明を行っていきます:

• 式1の中でいきなり不等式を作るのではなく、 を使ったまま式を展開する
• 二次式が出てくるので最小値を調べてみる
• 共分散の符号を調べる補題に落とし込む

二次式を作るところまで

まず、改めて式1を展開をしていきます。

\begin{eqnarray} w_1 :&=& \mathbb{E}\bigl[\frac{1}{\pi}\bigr] \\ w_0: &=& \mathbb{E}\bigl[\frac{1}{1-\pi}\bigr] \\ \mathbb{E}\Bigl[ \frac{\theta_1^{2}}{\pi} \Bigr] + \mathbb{E}\Bigl[ \frac{\theta_0^{2}}{1-\pi} \Bigr] - (\theta_1 - \theta_0)^2 &=& \bigl(w_1- 1\bigr) \theta_1^2 + \bigl(w_0 - 1\bigr) \theta_0^2 + 2\theta_1\theta_0 \tag{3} \end{eqnarray}

ここで、 なので、両辺を で割っても符号は変わりません。実際に割ってみると、

\begin{eqnarray} \biggl(\theta_1 + \frac{\theta_0}{w_1-1}\biggr)^2 + \frac{w_0 - 1}{w_1 - 1}\theta_0^2 - \frac{1}{(w_1-1)^2}\theta_0^2 &\ge& \frac{\theta_0^2}{(w_1-1)^2}\bigl( (w_1-1)(w_0-1)-1\bigr) \tag{4} \end{eqnarray}

のように計算でき、Theorem 1 が成立するには以下の条件が成立すればよいことがわかります:

\begin{eqnarray} \bigl(w_1-1)(w_0-1) \ge 1 \tag{5} \end{eqnarray}

共分散の符号を調べる補題への落とし込み

式5に をかけた上で再度 を使って書いてみると、以下のように計算できます:

\begin{eqnarray} -\mathbb{E}\bigl[\frac{1}{\pi}\bigr] \mathbb{E}\bigl[\frac{1}{1-\pi}\bigr] + \biggl(\mathbb{E}\bigl[\frac{1}{\pi}\bigr] + \mathbb{E}\bigl[\frac{1}{1-\pi}\bigr] \biggr) &=& -\mathbb{E}\bigl[\frac{1}{\pi}\bigr] \mathbb{E}\bigl[\frac{1}{1-\pi}\bigr] + \mathbb{E}\bigl[\frac{1}{\pi}\frac{1}{1-\pi}\bigr]\\ &\le& 0 \tag{6} \end{eqnarray}

ここで、 はランダムな共変量によって決まる確率変数であるので、式6は と の共分散がゼロ以下であるという意味になります。

よって、共分散についての補題2を示すことにより、Theorem 1 が示せます*2。

補題 1

Lemma 1

\begin{eqnarray} \text{Let} \ C \ \text{be a covariance of two random variables} \\ \text{and} \ X \ \text{be a random variable s.t.} \ 0 < X.\ \text{Then,}\\ C\biggl(X, \frac{1}{X}\biggr) &\le& 0 \tag{7} \end{eqnarray}

Proof of Lemma 1

\begin{eqnarray} C\biggl(X, \frac{1}{X} \biggr) &=& \mathbb{E}\bigl[X \cdot \frac{1}{X}\bigr] - \mathbb{E}[X]\mathbb{E}\bigl[\frac{1}{X}\bigr] \\ &=& 1 - \mathbb{E}[X]\mathbb{E}\bigl[\frac{1}{X}\bigr] \\ &\le& 0 \ \ \ (\because \text{Jensen's inequality and} \ X > 0) \tag{8} \end{eqnarray}

補題 2

Lemma 2

\begin{eqnarray} \text{Let} \ C \ \text{be a covariance of two random variables} \\ \text{and} \ X \ \text{be a random variable s.t.} \ 0 < X < 1.\ \text{Then,}\\ C\biggl(\frac{1}{X}, \frac{1}{1-X}\biggr) &\le& 0 \tag{7} \end{eqnarray}

Proof of Lemma 2

\begin{eqnarray} t_1 &:=& \frac{1}{X}\\ t_0 &:=& \frac{1}{1-X} \tag{8} \end{eqnarray}

とおくと、以下の関係式を得られます:

\begin{eqnarray} t_1 &>& 1\\ t_0 &>& 1\\ \frac{1}{t_1} &=& 1 - \frac{1}{t_0}\\ t_0 &=& \frac{t_1}{t_1 - 1}\\ &=& 1 + \frac{1}{t_1 - 1} \tag{9} \end{eqnarray}

これらを用いると、期待値の線形性や共分散の性質から、

\begin{eqnarray} C\biggl(\frac{1}{X}, \frac{1}{1-X}\biggr) &=& \mathbb{E}\biggl[\frac{1}{X} \cdot \frac{1}{1-X}\biggr] - \mathbb{E}\biggl[\frac{1}{X}\biggr]\mathbb{E}\biggl[\frac{1}{1-X}\biggr] \\ &=& \mathbb{E}\biggl[ t_1 \cdot \biggl(1 + \frac{1}{t_1 - 1}\biggr)\biggr] - \mathbb{E}\biggl[ t_1 \biggr] \mathbb{E}\biggl[1 + \frac{1}{t_1 - 1}\biggr] \\ &=& \mathbb{E}\biggl[\frac{t_1}{t_1 - 1}\biggr] - \mathbb{E}\biggl[ t_1 \biggr] \mathbb{E}\biggl[\frac{1}{t_1 - 1}\biggr] \\ &=& C\biggl(t_1, \frac{1}{t_1 - 1}\biggr) \\ &=& C\biggl(t_1-1, \frac{1}{t_1 - 1}\biggr) \\ &\le& 0 \ \ \ (\because \text{Lemma 1}) \tag{10} \end{eqnarray}

となるので、題意が示されました*3。

まとめ

以上のように、推定量の漸近分散の比較を詳しく行い、傾向スコアモデルgivenのときには「Horvitz-Thompson estimator より Hajek estimator の方が漸近分散が大きくなることはない」ことを示しました。

共分散の符号の補題については、よく使いそうなものであるにも関わらずピンポイントな主張をすぐに見つけられず、自分で示した方が早いだろうと判断して証明してみました。

証明の誤りを見つけられた方、あるいは補題の証明について言及されているよりよい資料を見つけられた方がいらっしゃいましたら、そっとご指摘いただけると幸いです。

このような記事を書くきっかけを作ってくださった @BluesNoNo さん、ありがとうございます。

*1:厳密に新規性についてサーベイをできているわけではないのでご容赦ください

*2:Twitter上のやりとりでは「直感的に成立しそう」というところまでしか言えていませんでした

*3:2行目から3行目の展開で ] がキャンセルされます

因果推論における standardization と parametric g-formula の関係性

過去の記事では、因果推論の中で重要な手法の一つである IPW の概要や疑問点について整理してきました。

今回は、 IPW と双璧をなす重要な手法である Standardization（標準化） について簡単に紹介します。

これまで同様、what if book を中心に議論していきます。

目次はこちらです。以前の記事よりはボリュームを抑えています。

what if book の主張を紹介した後で、数式部分を少しだけ補足します。

記事のゴール

本記事では、single time point exposure の条件下での平均因果効果の推定にあたって、以下の主張を行うことをゴールにします:

• what if book 2章の標準化と13章の parametric g-formula は、 outcome model について共変量で期待値をとっているという意味で統一的に見ることができるが、共変量の分布の扱い方によっては推定結果に違いが生じうる

Notation や仮定は冒頭に紹介した過去の記事と同じものを使います。

新規性

記事の主張について明確に述べている資料は見つからなかったので、一応それを新規性とします（ほんの少し式を補足しただけですが）。

what if book 2章で紹介された標準化のおさらい

標準化は、 what if book では2章と13章で登場します。

ここでは、2章の主張をおさらいします。

2章では、共変量が離散の場合に、以下のように potential outcome を推定できることが紹介されています:

\begin{eqnarray} \mathbb{E}[Y^{(a)}] &=& \sum_{l} \mathbb{E}[Y^{(a)} \mid L=l] \ \Pr[L=l] \ \ \ (\because \text{definition of marginalization}) \\ &=& \sum_{l} \mathbb{E}[Y^{(a)} \mid L=l, A=a] \ \Pr[L=l] \ \ \ (\because \text{conditional exchangeability}) \\ &=& \sum_{l} \mathbb{E}[Y \mid L=l, A=a] \ \Pr[L=l] \ \ \ (\because \text{consistency}) \tag{1} \end{eqnarray}

最後の行の summation の中にある期待値と確率分布をそれぞれ観測データから推定すればよいということになり、以下の式で最終的な potential outcome の推定を行います:

\begin{eqnarray} \sum_{l} \hat{\mathbb{E}}[Y \mid L=l, A=a] \ \hat{\Pr}[L=l] \tag{2} \end{eqnarray}

what if book 2章の標準化と13章の parametric g-formula の対応

次に、what if book 2章と13章の対応を考え、それらを統一的に見た上で、違いが生じうる部分について簡単に考察します。

13章のおさらい

13章では以下の2段階で議論が展開されています:

• 標準化は、outcome model について共変量で期待値をとる手法として考えられる
• outcome regression の予測結果について sample mean をとることでその期待値をいい感じに推定できる

打ち切り変数 については無視した上で、議論を追っていきます。

共変量で期待値をとっているという部分は以下の式で書けます:

\begin{eqnarray} \mathbb{E}[Y^{(a)}] &=& \mathbb{E}_{l} [\mathbb{E}[Y^{(a)} \mid L=l]] \ \ \ (\because \text{definition of marginalization}) \\ &=& \mathbb{E}_{l} [\mathbb{E}[Y^{(a)} \mid L=l, A=a]] \ \ \ (\because \text{conditional exchangeability}) \\ &=& \mathbb{E}_{l} [\mathbb{E}[Y \mid L=l, A=a]] \ \ \ (\because \text{consistency}) \tag{3} \end{eqnarray}

によって共変量で期待値をとっていることを表現しました。

2章の標準化では共変量が離散の場合の表現しか与えられていませんでしたが、期待値の形で書くことで、共変量が連続の場合も表現できるようになりました。

その上で、 の期待値については outcome regression によって推定を行い、共変量で期待値をとる部分については観測データによる経験近似を行う以下の式が紹介されています:

\begin{eqnarray} \frac{1}{n} \sum_{i=1}^{n} \hat{\mathbb{E}}[Y \mid L=L_i, A=a] \tag{4} \end{eqnarray}

これが parametric g-formula です。

outcome のモデルが特定されていて、かつ観測データが i.i.d. で得られていることを仮定すれば、求めたい期待値をいい感じに推定できそうです。

共変量については観測データ の値を使い、処置については興味のある処置 の値を固定している（処置については観測データではない）というところが肝になります。

2章との対応

what if book を読まれた方の中には、2章の標準化と13.3節の parametric g-formula の数式（式2と式4）間の対応がピンと来ていない方がいるかもしれないので、簡単に説明します。

結論としては、共変量の分布をどのように扱うかによって、同値になったりならなかったりする という話になります。

共変量が離散の場合、共変量の分布は以下のように推定されることが多いでしょう（ただの集計です）:

\begin{eqnarray} n_l &:=& \sum_{i=1}^{n} I(L_i, l) \\ \hat{\Pr}[L=l] &=& \frac{n_l}{n} \tag{5} \end{eqnarray}

このとき、式2は以下のように展開できます:

\begin{eqnarray} \sum_{l} \hat{\mathbb{E}}[Y \mid L=l, A=a] \ \hat{\Pr}[L=l] &=& \sum_{l} \hat{\mathbb{E}}[Y \mid L=l, A=a] \ \frac{n_l}{n}\\ &=& \sum_{l} \frac{\sum_{i=1}^n I(L_i, l)\hat{\mathbb{E}}[Y \mid L=L_i, A=a]}{n_l} \ \frac{n_l}{n}\\ &=& \sum_{l} \frac{\sum_{i=1}^n I(L_i, l)\hat{\mathbb{E}}[Y \mid L=L_i, A=a]}{n}\\ &=& \frac{1}{n} \sum_l \sum_{i=1}^n I(L_i, l)\hat{\mathbb{E}}[Y \mid L=L_i, A=a] \\ &=& \frac{1}{n} \sum_{i=1}^n \hat{\mathbb{E}}[Y \mid L=L_i, A=a] \tag{6} \end{eqnarray}

最後の行は式4と全く同じ形になるので、式5のように共変量の分布を推定した場合、2章の標準化と13章の parametric g-formula は同じものとして対応がとれることがわかります。

逆に、2章の標準化にあたって共変量の分布の推定方法を変えた場合（例: 適当な prior を置いてベイズ推定する）、parametric g-formula による推定結果とは異なる結果が得られうることがわかります。

同様に、共変量が連続であるときに、「共変量の分布を推定した上でその分布を明示的に使って期待値をとるという手法」と 「parametric g-formula のように共変量の分布を推定しない手法」を比較すると、やはり推定結果が異なりうると言えるでしょう。

おまけ: conditional effect と標準化

標準化は marginal effect を推定するための手法であるため、基本的には conditional effect の推定には使えません*1。

一般に、 conditional effect を推定する際には、outcome regression や Marginal Structural Model や g-estimation を利用することになると思われます。

まとめ

what if book の2章と13章の数式を繋ぐことを目的に記事を書いてみましたが、意外と非自明なことが少なく、かなり内容の薄い記事になってしまいました。

もし同じような引っかかりを感じた方の理解の助けになれば幸いです。

*1:marginal effect を推定するために、 conditional effect を推定する手法の一つである outcome regression に修正を加えたものが標準化という認識です

色々なIPWの整理 〜 Stabilized weights の使い所 〜 (2)

前回の記事に引き続き、 IPW についての連載第2弾です。

今回も what if book を中心に議論していきます。

目次はこちらです。

IPW や Stabilized weights についての what if book の主張を紹介した後で、Stabilized weights の使い所についていくつかの角度から説明していきます。

記事のゴール

前回の記事では、single time point exposure の条件下での平均因果効果の推定にあたって、以下の2つの主張について説明しました:

1. IPW によって処置 の Potential outcome ] を推定するときの基本は、 ] を estimand とした Horvitz-Thompson estimator の推定である
2. Hajek (modified Horvitz-Thompson) estimator は不偏推定量ではないが、いくつかの仮定のもとで Horvitz-Thompson estimator よりも漸近分散が小さい

本記事では、前回の記事で4つ目の主張として紹介していた以下の内容について説明することをゴールにします:

• Stabilized weights は Hajek estimator を推定するときに有用かもしれないが、他の推定量を推定するときには注意が必要 (特に、Doubly robust では適切な使用方法がわからない)

Notationは前回の記事を踏襲します。

新規性

新規性が気になる方のために、既存の資料と比較しておきます。

以下の点で新規性があると考えています:

• IPW において、 stabilized weights を使うことが妥当と思われる条件を調べ、 stabilized weights の安易な利用に警鐘をならした
• what if book における IPW と exchangeability と stabilized weights の関係性の不明点を列挙した（列挙しただけ）
• Doubly robust 推定量において stabilized weights を使うべきではないということを理論的な背景から説明した

前提知識

前回の記事 (1) の Notation と、 Horvitz-Thompson estimator と Hajek estimator の定義を理解していることを前提知識にします。

IPWの目的

Stabilized weights を理解するには、「逆確率で重みづけすることの意味や目的」について整理する必要があると考えているので、その部分の説明から始めます。

what if book の2章の紹介

前回の記事でも軽く触れたのですが、 IPW の目的は「データを augmentation することによって、exchangeability が成立するような擬似データセットを作ること」と言えます。

この件について、以前勉強会で私が発表を行ったときに、「IPW はピコ太郎だ」という旨のスライドを作成しました:

IPWがデータを augmentation していることの直感的な説明

スライドに掲載している図はすべて what if book の2章から拾ってきている*1のですが、「介入が行われた世界線と、行われなかった世界線を考えて、それらを（ピコ太郎のように）合体させた pseudo population（擬似データセット） を作ること」を IPW の目的だと主張するものです。

逆確率を用いてサンプルを重みづけすれば、そのような擬似データセットを作成することができます。

擬似データセットの特徴としては、 what if book の Figure 2.3 の手前で以下のように述べられています:

Under conditional exchangeability in the original population, the treated and the untreated are (unconditionally) exchangeable in the pseudo-population because the L is independent of A

つまり、共変量と処置が独立になり、（条件づきではない）exchangeability が成立するということです。

そして、以下の文が続きます:

That is, the associational risk ratio in the pseudo-population is equal to the causal risk ratio in both the pseudo-population and the original population

exchangeability が成立することで、擬似データにおける関連のリスク比と、擬似データにおける因果のリスク比と、元のデータにおける因果のリスク比がすべて一致すると述べられています。

なお、処置が二値の場合、因果のリスク比は以下のように定義されます:

\begin{eqnarray} \frac{\mathbb{E}[Y^{a=1}] }{ \mathbb{E}[Y^{a=0}] } \tag{1} \end{eqnarray}

さらに、 Technical Point 2.3 では、逆確率を用いて擬似データセットを作ることで、リスク比が一致するだけでなく、 Potential outcome の期待値についての以下の等式も成立することが示されています:

\begin{eqnarray} \mathbb{E}[Y^{(a')}] = \mathbb{E} \Bigl[ \frac{I(A, a')Y}{f_{A \mid L}(a' \mid L)} \Bigr] \tag{2} \end{eqnarray}

この等式が成立することから、因果のリスク比が関連のリスク比と一致することも自明ですし、因果のリスク差（ATE）が関連のリスク差と一致することも自明なので、 IPW の有用性を感じられるかと思います。

2章の主張への疑問点

先ほど「exchangeability が成立する」と述べましたが、 what if book の中で厳密な証明を発見することができませんでした*2。

ただし、これまでの流れから察するに、共変量と処置の独立性からpseudo-population の世界において確率変数の意味での exchangeability が成立するということは言えるかもしれません。

また、式 (2) によって 「pseudo-populationでの期待値と元のデータのpotential outcomeの期待値が一致する」 という（2つのデータ間のmean exchangeability のような）ことは言えたと思います。

しかし、pseudo-population の世界において確率変数の意味での exchangeability が成立することが何に使えるのか ということは明記されておらず、なぜこのようなことを考えたいのかが不明瞭であるように感じました（Technical Point 2.3が pseudo-population の世界における確率変数の意味での exchangeability とどう結びつくのがわからないです）。

今後この記事の中で何度か exchangeability についての言及を行いますが、このあたりが曖昧なままになっているということを前提にお読みいただけると幸いです。

what if book 12.3 節の紹介

Stabilized weights について最も丁寧に述べている資料は what if book 12.3 節だと思うので、そこで述べられている内容を紹介します。

12.3 節の主張をまとめると「IPW の目的を踏まえると、必ずしも逆確率で重みづけを行わなくてもよい」というものになります。

色々な擬似データセットの生成

12.3節の冒頭では、 IPW の目的が共変量と処置の関連がないような（exchangeabilityが成立するような）擬似データセットの作成だということが再度強調されています。

The goal of IP weighting is to create a pseudo-population in which there is no association between the covariates L and treatment A

そして、逆確率による重みづけは、 そのような特徴を持つデータセットのクラスの中で、元のデータの2倍の大きさの擬似データセットを作成していることになると説明されています:

\begin{eqnarray} \mathbb{E}_{A}\Bigl[\frac{1} { f_{A \mid L}(A \mid L=l) } \Bigr] = \frac{\Pr[A=1 \mid L=l] }{ f_{A \mid L}(1 \mid L=l)} + \frac{ \Pr[A=0 \mid L=l] }{ f_{A \mid L}(0 \mid L=l)} = 2 \tag{3} \end{eqnarray}

しかし、ここで「共変量と処置が独立になるような擬似データセットを作る方法は他にもある」という主張がされています:

However, there are other ways to create a pseudo-population in which L and A are independent.

その直後に直感的な説明として、「逆確率に一律 をかけたもので重みづけをしても、 exchangeability への影響はなくて、かつ擬似データセットのサイズは元のデータと同じになる」という主張がされています。

詳しくは、擬似データセットについての定理により正当化されています。

擬似データセットについての定理

Technical Point 12.2 では、任意の関数 （正の値をとる必要がある）を使って重みづけを行うことに妥当性があるということを、Estimand が Potential outcome の期待値と等しくなる式によって説明しています。

すなわち、以下が成立することをもって、どんな関数を使っても妥当だという話になっています。

\begin{eqnarray} \frac{\mathbb{E} \Bigl[ \frac{I(A, a')Y}{f_{A \mid L}(a' \mid L)} g(A) \Bigr]} {\mathbb{E} \Bigl[ \frac{I(A, a')}{f_{A \mid L}(a' \mid L)} g(A) \Bigr]} &=& \frac{\mathbb{E}[Y^{(a')}] g(a')} {g(a')} \\ &=& \mathbb{E}[Y^{(a')}] \tag{4} \end{eqnarray}

どの関数を使うべきか

関数 として、処置の確率 ないしはその推定結果を使うべきだという主張がされています。

このとき、実際の推定は以下のように行われます:

\begin{eqnarray} \frac{\sum_{i=1}^{n} \Bigl[ \frac{\hat{f}(A_i)}{\hat{f}_{A \mid L}(A_i \mid L_i)}I(A_i, a')Y_i \Bigr]} {\sum_{i=1}^{n} \Bigl[ \frac{\hat{f}(A_i)}{\hat{f}_{A \mid L}(A_i \mid L_i)} I(A_i, a')\Bigr]} \tag{5} \end{eqnarray}

ここで使われた重みのことを stabilized weights と定義しています。

Stabilized weights を使うべき理由として、 Estimand が変わらないという条件を満たしつつ、一定の条件のもとで信頼区間が小さくなることを挙げられています。

また、effect modifier が存在するときは、effect modifier で条件づけた処置確率の推定結果を使うことが良しとされています。

what if book 12.3 節の補足と修正提案

以上をご覧になった方は、stabilized weights は素晴らしいもので、 IPW を考えるときにはなんでもかんでも stabilized weights を使えばいいような気持ちになるかもしれません。

しかし、ここから先では、 stabilized weights をそのまま使うと誤った解釈に至りうるケースについての指摘を行おうと思います。

ここでの結論は以下のようになります:

• Technical Point 12.2 は Hajek estimator を Estimand にした場合の話であり、それ以外の推定量について同様の等式が成り立つという保証はないのでは？
• 例えば、 Horvitz-Thompson estimator を Estimand にしたときは補正が必要になる（というかそもそも Horvitz-Thompson estimator で stabilized weights を使っても意味がない）のでは？

擬似データセットについての定理の補足

Technical Point 12.2 では、 Estimand が Potential outcome の期待値と等しくなる式を紹介しています。

しかし、前回の記事を読んでいただいた方はわかると思うのですが、これは Hajek estimator を Estimand とした場合の式展開になります。

Horvitz-Thompson estimator を Estimand にしたときは以下のように計算できます:

\begin{eqnarray} \mathbb{E} \Bigl[ \frac{I(A, a')Y}{f_{A \mid L}(a' \mid L)} g(A) \Bigr] &=& \mathbb{E}[Y^{(a')}] g(a') \\ \tag{6} \end{eqnarray}

つまり、stabilized weights を使った場合、 の逆数をかける必要が出てきます。

ただし、式を見るとわかるように、計算結果全体に対して を一律でかけただけになっているので、そもそもこの場合は stabilized weights を使う意味がありません。

このように、なんでもかんでも stabilized weights を使うべきとは言えない上に、誤った使い方をすると推定結果がおかしな値になってしまう可能性があります。

Technical Point 12.2 ではこうした前提条件が明示的に記述されていなかったため、誤った使い方をしてしまうリスクが一定存在すると考えられ、本記事で補足を行う価値があると判断しました*3。

Technical Point 12.2 以外の記述について

Stabilized weights を使うことの妥当性が Technical Point 12.2 に依存しているのであれば、 12.3 節の中で直感的な説明を行う部分についても、Technical Point 12.2 の前提条件（Hajek estimator の推定であること）を明示した方がよいと考えています*4。

逆にいうと、「他の推定量を扱う場合に、stabilized weights を使うことの妥当性がどのような直感で説明されるのか」という説明があるとより親切だと感じました*5。

また、 Hajek estimator を紹介するときに、 「重みづけするweight のスケーリングに対して、少なくとも Estimand の世界では計算結果が不変*6」といった特徴を記載してもよいのではないかと感じました。

Doubly robust 推定量と 様々な IPW 推定量

IPW を使った手法として、Doubly robust という便利なものが知られています。

記事の締めくくりとして、 Doubly robust 推定量 と 様々な IPW 推定量の関係性についての私見を述べます。

結論としては、以下になります:

• Horvitz-Thompson estimator 以外の IPW 推定量を Doubly robust で使うのは難しい（方法がわからない）
• 例えば、 Stabilized weights をむやみに Doubly robust で使うのは危険

Doubly robust 推定量

IPW では、傾向スコアモデルが正しく推定されていた場合、因果効果の推定結果にバイアスがのらないという保証はできません。

また、アウトカム回帰モデルを利用した Standardization*7 では、回帰モデルが誤って推定されていた場合、因果効果の推定結果にバイアスがのらないという保証はできません。

そこで、「傾向スコアモデルかアウトカム回帰モデルの少なくとも一方が正しく推定されていれば因果効果の推定結果にバイアスがのらないような推定量」として Doubly robust 推定量が提案されています*8。

処置が二値のとき、処置 の Potential outcome の期待値の Doubly robust 推定量は以下のように定義できます:

\begin{eqnarray} \hat{\mathbb{E}}_{DR}[Y^{(a')}] &:=& \frac{1}{n} \sum_{i=1}^{n}\Bigl[\frac{I(A_i, a')Y_i}{\hat{f}_{A \mid L}(a' \mid L_i)} + \Bigl(1 - \frac{I(A_i, a')}{\hat{f}_{A \mid L}(a' \mid L_i)}\Bigr) \hat{f}_{Y\mid A, L}(a', L_i)\Bigr] \tag{7} \end{eqnarray}

what if book では Technical Point 13.2 で紹介されており、 Horvitz-Thompson estimator を アウトカム回帰モデルを伴った関数で補正しているとみなせます:

which can also be viewed as a correction of the Horvitz-Thompson estimator by a function that involves the outcome regression model

式を見ると明らかですが、 summation の中身の一項目は Horvitz-Thompson estimator と同じ形をしています。

また、二項目において、 Standardization を傾向スコアにより補正しています。

要するに、IPW と Standardization の合わせ技が Doubly robust ということになります。

Horvitz-Thompson estimator 以外の IPW 推定量を Doubly robust で使えるか？

さて、IPW を学んできた立場としては、以下のことが気になるのは自然なことかと思います:

• Doubly robust で、 Stabilized weights は使えるのか？
• Doubly robust で、 Horvitz-Thompson estimator 以外の形の推定量も使えるのか？

これについて説明した資料はなかなか見つけられなかったのですが、私の意見は「自明ではなく、難しそう」です。

いくつかの具体例をもとに考察していきます。

例えば、stabilized weights を 式 (7) に導入する方法についていくつか考えてみます。

まず、逆確率の部分を直接処置確率で補正すると、このようになると思います:

\begin{eqnarray} \frac{1}{n} \sum_{i=1}^{n}\Bigl[\frac{I(A_i, a')\hat{f}_{A}(a')Y_i}{\hat{f}_{A \mid L}(a' \mid L_i)} + \Bigl(1 - \frac{I(A_i, a') \hat{f}_{A}(a') }{\hat{f}_{A \mid L}(a' \mid L_i)}\Bigr) \hat{f}_{Y\mid A, L}(a', L_i)\Bigr] \tag{8} \end{eqnarray}

しかしこの式は、「傾向スコアモデルが正しく推定できているがアウトカムモデルが誤って推定されている」ときにバイアスが生じてしまうので、 Doubly robust とは言えません。

また、式全体に対して処置確率で補正をするとこのようになります。

\begin{eqnarray} \frac{1}{n} \sum_{i=1}^{n}\Bigl[\frac{I(A_i, a')\hat{f}_{A}(a')Y_i}{\hat{f}_{A \mid L}(a' \mid L_i)} + \hat{f}_{A}(a') \Bigl(1 - \frac{I(A_i, a') }{\hat{f}_{A \mid L}(a' \mid L_i)}\Bigr) \hat{f}_{Y\mid A, L}(a', L_i)\Bigr] \tag{9} \end{eqnarray}

しかしこの式は、「傾向スコアモデルが誤って推定されているがアウトカムモデルが正しく推定されている」ときにバイアスが生じてしまう上、全体に処置確率をかけているだけなので、 Doubly robust でない上に stabilization の意味もありません。

また、 summation の中身の一項目（Horvitz-Thompson estimator）、あるいは 二項目（Standardization の補正）のどちらか一方にのみ補正を行う場合、色々なケースで明らかにバイアスが出てきてしまうので不適切です。

よって、 Doubly robust で Horvitz-Thompson estimator を扱う場合、 Stabilized weights を使うことは難しいと考えられます。

最後に、Doubly robust で Hajek estimator のようなものを扱う場合を考えてみます。

この記事の前半で述べたように、 Stabilized weights の有用性は Hajek estimator を扱うときに発揮されるので、 Doubly robust でも Hajek estimator のようなものを考えれば Stabilized weights を使えるかもしれないという発想に至るでしょう。

素朴に考えると、以下のように分母をつけたくなるかと思います。

\begin{eqnarray} \frac{ \sum_{i=1}^{n}\Bigl[\frac{I(A_i, a')Y_i}{\hat{f}_{A \mid L}(a' \mid L_i)} + \Bigl(1 - \frac{I(A_i, a')}{\hat{f}_{A \mid L}(a' \mid L_i)}\Bigr) \hat{f}_{Y\mid A, L}(a', L_i)\Bigr] } { \sum_{i=1}^{n} \Bigl[ \frac{I(A_i, a')}{\hat{f}_{A_i \mid L}(a' \mid L_i)} \Bigr] } \tag{10} \end{eqnarray}

しかし、前回の記事で指摘したように、これは ratio estimator の特殊ケースなので、傾向スコアモデルが正しく推定されていたとしても、不偏性は持たないというのが私の意見です。

またそもそも、傾向スコアモデルが誤って推定されていたとき、分母の期待値は とは大きく異なることが予想されるので、式 (10) は傾向スコアモデルに対してロバストではない（すなわち、 Doubly robust ではない）推定量と言えるかと思います。

よって、 Hajek estimator のような形の Doubly robust 推定量を作ることは難しく、式 (10) に対する stabilized weights を考えるまでもなくお蔵入りになってしまうでしょう。

まとめ

本記事では、 IPW の文脈でよく紹介される stabilized weights という重みづけ方法の特徴と、その注意点をまとめました。

一通り読んでいただいた方には、「なんでもかんでも stabilized weights を使えばいい」というわけではないことがお伝えできたかと思います。

今後は、より一般的な連続値の介入を想定した手法や、実験による主張の検証についての記事を書く予定をしています。

記事を読んでいただいた皆さんの豊かな因果推論ライフに貢献できれば幸いです。

*1:もちろんピコ太郎氏の写真は what if book には載っていません。Avex さんの HP (https://avex.jp/pikotaro/) から借りてきました

*2:ご存知の方は教えていただけると幸いです

*3:Technical Point 12.1 で「この章では Hajek esitmator を推定する」とは書かれていますが、これが stabilized weights を使うにあたっての本質的な前提条件であることに気付いている読者は多くはないと思います・・・

*4:これは先ほどの「pseudo-population の世界において確率変数の意味での exchangeability が成立することが何に使えるのかわからない」という疑問とも関連していて、pseudo-population の世界において確率変数の意味での exchangeability が成立するような stabilized weights を使ったとしても、結局は別の基準（不偏性や漸近分散）をもとに議論するのであれば、pseudo-population における exchangeability とはいったい何なんだろう・・・という気持ちになりました

*5:私も適切な説明はまだ思いついていないです

*6:「不偏」の間違いではなく invariant の意味です

*7:g-formula として what if book の 13章で紹介されています

*8:この記事では議論を簡略化するために「真の関数が得られている」という意味での「正しく推定されている」という条件下での不遍性を考えました。しかし、Doubly robust 推定量の性質を語るときには、モデルが特定されている」という条件下での漸近的な不偏性（これを一致性と呼んでいる文献もある気がします）を持つものとしての議論をすることが多いと思われます。適宜読み替えていただけると幸いです

色々なIPWの整理 〜 Causal inference: what if の12章を中心に 〜 (1)

Hernan and Robins によって書かれた因果推論の本（本記事では "what if book" と呼びます。著者のサイト にPDFが無料で公開されています）を読んでいく中で、IPW (Inverse Probability Weighting; 逆確率重みづけ) 周りの細かいところを整理するのが難しいなと感じたので、自習のついでに記事としてまとめようと思います。現時点では十分な数の関連論文を読めているわけではないので、調査に抜け漏れがある可能性についてはご容赦いただけると幸いです。

アイキャッチとして使った紀伊国屋のリンクも貼っておきます:

目次はこちらです

記事のゴール

本記事は、IPW に関する連載のスタートの記事です。

本記事では、 single time point exposure の条件下での平均因果効果の推定にあたって、 IPW に関連する以下の4つの主張のうち、1と2について説明を行うことをゴールにします。今後の連載の中で、3と4の説明、ならびに人工データ実験による検証を行う予定です。

1. IPW によって処置 の Potential outcome ] を推定するときの基本は、 ] を estimand とした Horvitz-Thompson estimator の推定である
2. Hajek (modified Horvitz-Thompson) estimator は不偏推定量ではないが、いくつかの仮定*1のもとで Horvitz-Thompson estimator よりも漸近分散が小さい
3. Marginal Structural Model のパラメータ推定を行う際に、逆確率で重み付けた最小二乗法を用いることが多いが、これは Hajek estimator の推定に相当する
4. Stabilized weights は Hajek estimator を推定するときに有用かもしれないが、他の推定量を推定するときには注意が必要 (特に、Doubly robust では適切な使用方法がわからない)

要するに、私が以前 Twitterで 投稿していた内容を整理します。Notation はすぐ後で与えます。推定手法については適宜説明します。

<今の理解>
1. single time pointのIPWの基本はE[I(A=a)Y/f(A|L)]をestimandとしたHorvitz-Thompson estimator (HTe) による推定
2. modified HTeの方がHTeより安定
3. IP weighted LS = modified HTeの推定
4. stabilized weightはmodified HTeを推定するときに利用可 (→Doubly robustでは不可)

— fullflu (@fullflu) April 29, 2020

新規性

新規性が気になる方のために、既存の資料と比較しておきます。

以下の点で新規性があると考えています:

• 因果推論についての有名な資料間での、推定量についての記述の相違を整理した
• 「Hajek estimator は Potential outcome の期待値についての不偏推定量である」と主張する what if book や星野先生の著書の証明に不備がある可能性を指摘した
• Horvitz-Thompson estimator の漸近分散を解析し、 Hajek estimator の方が漸近分散が小さいことを保証できる仮定（Potential outcome の符号一致）を調べた

前提知識

本記事の理解のために最低限必要な前提知識は2つあると考えます。

1つ目は、Potential outcome を用いた因果推論の枠組みを勉強した経験があることです。因果推論そのものの説明は本記事では軽く導入する程度です。

2つ目は、学部教養レベルの統計学を理解していることです。例えば、期待値や不偏性といった概念の説明は本記事では行わないです。

本記事の最後に登場する漸近分散の解析においてのみ、M推定量やモーメント法についての知識を仮定しますが、興味のない方はスルーしていただいても概要の理解には差し支えないかと思います。

Notation

必要なNotationを与えます。

• : Outcomeの確率変数です。実数の連続値をとるとします。
• : 処置の確率変数です。連続か離散か二値かは都度説明します。
• : 共変量の確率変数です。 次元のベクトルとします。
• : Effect modifierの確率変数です。簡単のため、共変量のsubsetである場合を考えます。
• : 処置 についての Potential outcome を意味する確率変数です*2。
• : 指示関数です。一つ目の引数と二つ目の引数が等しいときに1を、そうでないときに0を返します。
• : 処置が離散のとき、共変量の値が であったときの処置 の確率分布です。処理が連続のとき、 のもとでの処置 の確率密度関数です。
• : から共変量についての条件を除いたものです。
• : 共変量Lと処置Aを引数として、Outcomeを出力する関数です。
• single time point exposure: 1時点の介入のみ考慮する実験設定のことです。
• : 処置が二値のときの処置確率を簡易表記するためのものです。

その他細かいルールとして、実際のデータを扱う際には、サンプルサイズを とし、確率変数に対して添え字 ] をつけることによって実現値を表現します。 例えば、 行目のレコードの処置の実現値は として表現します。

また、関数 の推定結果を で表現します。

IPW による Potential outcome の推定の基本: Horvitz-Thompson estimator

ここでは、冒頭に述べた主張の一つ目の話、Horvitz-Thompson estimator について説明します。

必要な仮定

本題に入る前に、因果効果を適切に推定するために必要な仮定について説明します。

ここからの話の中では、以下の6つの条件が成立することを仮定します。詳細は本記事では説明しないので、気になる方は what if book などを確認してみてください。

• SUTVA: 「ある個人への処置が、別の個人のoutcomeに影響しない」ことと「処置にバリエーションがない」ことの仮定です。
• Positivity: 共変量 で条件づけたときに、すべての処置 に対して、処置を受ける確率がゼロより大きいことを仮定します（特に断りがない限り、本記事では処置が二値のときを考えます）。
• Consistency: 実際に観測されたoutcomeは、実際に行われた処置に対応するoutcomeと一致することを仮定します。
• Conditional exchangeability: 共変量で条件づけたときに、Potential outcomesと処置変数は確率変数の意味で独立であるという仮定です。
• No selection bias: 選択バイアスが存在しないことを仮定します。
• No measurement error: 測定誤差が存在しないことを仮定します。

ちなみに、Conditional exchangeability を式で書くとこのようになります:

IPWを用いた因果効果の推定

本記事では、因果効果として、 Potential outcomes の期待値の差分による定義（Average Treatment Effect; ATE）を採用します。 特に、処置が二値の場合は、以下のように定義されます:

\begin{eqnarray} ATE := \mathbb{E}[Y^{a=1}] - \mathbb{E}[Y^{a=0}] \tag{2} \end{eqnarray}

定義の詳細が気になる方は what if book の1章や2章あたりをご覧ください。

6つの条件の仮定のもとで、Potential outcome についての以下の等式が成立することが知られています。

\begin{eqnarray} \mathbb{E}[Y^{(a')}] = \mathbb{E} \Bigl[ \frac{I(A, a')Y}{f_{A \mid L}(a' \mid L)} \Bigr] \tag{3} \end{eqnarray}

等式の左側の期待値は Potential outcome についてとっていて、右側の期待値は についてとっています。

証明は what if book の Technical Point 2.3 に書かれているので省略します。本記事とは少し表記が異なるので適宜読み替えてください。

とにかく、等式の右辺を Estimand*3 として何かしらの方法で推定することで、因果効果を推定できそうだということがわかりました。

Horvitz-Thompson estimator とは

what if book の Technical Point 12.1 に書かれているように、式 (2) の右辺を観測データを用いて以下のように推定（経験近似）したものは Horvitz-Thompson estimator と呼ばれています:

\begin{eqnarray} \hat{\mathbb{E}} \Bigl[ \frac{I(A, a')Y}{f_{A \mid L}(a' \mid L)} \Bigr] \tag{4} \end{eqnarray}

これは Potential outcome の期待値の不偏推定量かつ一致推定量になっています。

実際に計算するときには、傾向スコアのモデル を推定した後に、以下のように をプラグインして計算します。

\begin{eqnarray} \frac{1}{n}\sum_{i=1}^{n}\frac{I(A_i, a')Y_i}{\hat{f}_{A \mid L}[a' \mid L_i]} \tag{5} \end{eqnarray}

最後に、処置群と統制群で式 (4) の計算結果の差をとります。

ここまでが基本的な IPW による因果効果の推定になります。

以上を踏まえつつ、ここから先は、推定の改良を試みる手法について説明していきます。

Hajek (modified Horvitz-Thompson) estimator

ここでは、冒頭に述べた主張の二つ目の話、Hajek (modified Horvitz-Thompson) estimator について説明します。

Hajek estimator とは

what if book の Technical Point 12.1 を見ると、 Horvitz-Thompson estimator のすぐ下で、 modified Horvitz-Thompson estimator という以下のような推定量が紹介されているかと思います:

\begin{eqnarray} \frac{\hat{\mathbb{E}} \Bigl[ \frac{I(A, a')Y}{f_{A \mid L}(a' \mid L)} \Bigr]} {\hat{\mathbb{E}} \Bigl[ \frac{I(A, a')}{f_{A \mid L}(a' \mid L)} \Bigr]} \tag{6} \end{eqnarray}

どうやら、Horvitz-Thompson estimator の期待値の中から を削除したものを分母につけているようです。

IPW の本質は Data augmentation だと思っているのですが、「観測データの個数 ではなく、実際にどのくらい augmentation したのかに相当する項で割った方がよさそう」という気持ちを感じます。

what if book では推定量の名前は紹介されていませんでしたが、この推定量には Hajek estimator という名前がついているらしいので、以降 Hajek estimator と呼ぶこととします。

Hajek estimatorの性質

Hajek estimator の性質についての本記事での立場は以下となります:

• Hajek estimator は Potential outcome の期待値の不偏推定量とは限らないが、一致推定量ではある
• Potential outcome の符号が等しいとき、 Horvitz-Thompson estimator より Hajek estimator の方が漸近分散が小さい
• Hajek estimator は が二値である場合の推定結果が 0 と 1 の間に収まる（直感的に理解しやすい）

まず、漸近分散以外の部分について、参考文献とともに議論を整理します。

漸近分散以外の部分の議論の整理 (bias of Hajek estimator)

以前の記事で紹介した星野先生の著書（調査観察データの統計科学）の P70 では、関連内容についての記述が存在します。

調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

該当部分の記述を本記事の用語に翻訳すると、以下のようにまとめられます:

• Horvitz-Thompson estimator は Potential outcome の期待値の不偏推定量であり、一致推定量でもある
• Hajek estimator は Potential outcome の期待値の不偏推定量であり*4、一致推定量でもある
• Horvitz-Thompson estimator より、 Hajek estimator の方が推定精度が高い

また、 what if book の Technical Point 12.1 では以下のように述べられています:

• Positivity の仮定が成立するとき、 Hajek estimator は Potential outcome の期待値の不偏推定量である
• Hajek estimator は が二値である場合の推定結果が 0 と 1 の間に収まるので好ましい

一方、世界的に有名な今井先生の資料 (https://imai.fas.harvard.edu/teaching/files/weighting.pdf) では、以下のように述べられています:

• Weights are normalized but no longer unbiased （Hajek estimator で重みは正規化できるけど、不偏ではなくなる）

すなわち、不偏性についての記述に矛盾が起きているように感じます。

詳しく見ていくと、what if book では不偏性の証明が与えられていません。

また、星野先生の著書には証明が与えられていますが、 Hajek estimator の不偏性を示すべき文脈で、 Horvitz-Thompson estimator の不偏性を示す式 (3.9) が書かれてしまっています:

\begin{eqnarray} \mathbb{E}\Bigl[\frac{AY}{\pi} \Bigr] = \mathbb{E}[Y^{a=1}] \tag{3.9} \end{eqnarray}

つまり、Hajek estimator の不偏性については、 what if book ならびに 星野先生の著書では証明が与えられていないと考えられます。

証明を与える必要性を感じたため、似たような推定量の形を調べてみると、どうやら Ratio estimator というクラスとの関係性が深そうだとわかりました。

実際にRatio estimator の資料 (https://jkim.public.iastate.edu/teaching/book8.pdf) を読んでみると、 以下のことが示せます:

• Ratio estimator には一般に、式 (8.5) の形でバイアスが存在する
• Hajek estimator は Ratio estimator で とおいた場合に相当する

よって、本記事の中では「what if book や星野先生の著書の記述は誤りであり、 Hajek estimator には不偏性があるとは言えないのではないか」という立場をとります*5。

バイアスの具体的な形については、英語のブログ: Bias of the Hajek estimator (potential errors in Technical Point 12.1 of Causal inference: what if) - fullflu-english の後半に自分の計算結果を記載しておきました。

漸近分散の解析

星野先生の著書の中で「Horvitz-Thompson estimator より、 Hajek estimator の方が推定精度が高い」という主張がされていましたが、この部分について説明します。

著書の中ではこの主張の証明は行われていないため、Horvitz-Thompson estimator と Hajek estimator の漸近分散を比較することで、この主張に妥当性がありそうだということを示します。

ここから先はM推定量の理論解析を行うため、詳細に興味のある方のみご覧ください。

なお、傾向スコアについては真のモデルが given である場合を想定して計算を簡略化します*6。

Hajek estimator の漸近分散の整理

星野先生の著書の中では、処置が二値の場合の Hajek estimator の漸近分散の解析結果が与えられています。この議論の流れにのっとります。

漸近分散の解析においてM推定量の理論を使うのですが、こちらのスライドの27枚目以降の補足が非常にわかりやすいので、本記事でも参照させていただきます:

星野「調査観察データの統計科学」第3章 from Shuyo Nakatani

さて、まず星野先生の著書で与えられている Hajek estimator の漸近分散の主張（上記スライドの27枚目のスライドを復唱する形になり恐縮ですが）を本記事の Notation を使って整理します。

を真値 が である母数とし、以下の推定関数を考えます:

\begin{eqnarray} m(Y, \theta) := \Bigl(\frac{A}{\pi}(Y-\theta_1), \frac{1-A}{1-\pi}(Y-\theta_0) \Bigr)^{t} \tag{7} \end{eqnarray}

このとき、推定量 の 漸近分散（分散共分散行列） は になるという事実がM推定量についての一般論として知られています。ただし、

\begin{eqnarray} J(\theta^{*}) &:=& \mathbb{E}\Bigl[-\frac{\partial m(Y, \theta)}{\partial \theta^{t}} \Bigr]_{\theta=\theta^{*}}\\ K(\theta^{*}) &:=& \mathbb{E}[m(Y, \theta^{*}) m(Y, \theta^{*})^{t}] \tag{8} \end{eqnarray}

と定義しました。

これを計算すると、行列の要素は以下のようになります（ は単位行列になるので、 だけを考えることになります）:

\begin{eqnarray} Var(\theta^{*})_{1,1} &=& \mathbb{E}\Bigl[ \frac{(Y^{a=1}-\theta_1)^2}{\pi} \Bigr] \\ Var(\theta^{*})_{2,2} &=& \mathbb{E}\Bigl[ \frac{(Y^{a=0}-\theta_0)^2}{1-\pi} \Bigr] \\ Var(\theta^{*})_{1,2} &=& Var(\theta^{*})_{2,1} = 0 \tag{9} \end{eqnarray}

よって、平均因果効果の漸近分散は以下のように求まります:

\begin{eqnarray} Var_{Hajek}(\hat{\theta}_{1} - \hat{\theta}_{0}) &=& Var(\theta^{*})_{1,1} + Var(\theta^{*})_{2,2} - 2Var(\theta^{*})_{1,2}\\ &=& \mathbb{E}\Bigl[ \frac{(Y^{a=1}-\theta_1)^2}{\pi} + \frac{(Y^{a=0}-\theta_0)^2}{1-\pi} \Bigr] \tag{3.10改} \end{eqnarray}

これで復唱が終了したので、同様の方法で Horvitz-Thompson estimator の漸近分散の解析を行います。

Horvitz-Thompson estimator の漸近分散 (asymptotic variance of Horvitz-Thompson estimator)

Horvitz-Thompson estimator の漸近分散について、公開されている資料でわかりやすい解説がすぐには見つからなかったので、自分で解析してみます*7。

まず、このような推定関数を考えます:

\begin{eqnarray} m(Y, \theta) := \Bigl(\frac{AY}{\pi} -\theta_1, \frac{(1-A)Y}{1-\pi}-\theta_0 \Bigr)^{t} \tag{10} \end{eqnarray}

M推定量の枠組みで考えると、Hajek estimator の場合と同様に は単位行列になるので、 だけを考えればよくなります。

以下のように行列の要素を計算できます:

\begin{eqnarray} Var(\theta^{*})_{1,1} &=& \mathbb{E}\Bigl[ \frac{A^2{Y^{a=1}}^{2}}{\pi^2} + \theta_1^2 - 2 \frac{AY^{a=1}}{\pi}\theta_1 \Bigr] \\ &=& \mathbb{E}\Bigl[ \frac{{Y^{a=1}}^{2}}{\pi} + \theta_1^2 - 2Y^{a=1}\theta_1 \Bigr]\\ Var(\theta^{*})_{2,2} &=& \mathbb{E}\Bigl[ \frac{(1-A)^2{Y^{a=0}}^{2}}{(1-\pi)^2} + \theta_0^2 - 2 \frac{(1-A)Y^{a=0}}{1-\pi}\theta_0 \Bigr] \\ &=& \mathbb{E}\Bigl[ \frac{{Y^{a=0}}^{2}}{1-\pi} + \theta_0^2 - 2Y^{a=0}\theta_0 \Bigr]\\ Var(\theta^{*})_{1,2} &=& Var(\theta^{*})_{2,1}\\ &=& \mathbb{E}[\theta_1\theta_0 - Y^{a=1}\theta_0 - Y^{a=0}\theta_1] \tag{11} \end{eqnarray}

よって、 であることを使うと、平均因果効果の漸近分散は以下のようにまとめられます:

\begin{eqnarray} Var_{HT}(\hat{\theta}_{1} - \hat{\theta}_{0}) &=& Var(\theta^{*})_{1,1} + Var(\theta^{*})_{2,2} - 2Var(\theta^{*})_{1,2}\\ &=& \mathbb{E} \Bigl[ \frac{{Y^{a=1}}^{2}}{\pi} + \frac{{Y^{a=0}}^{2}}{1-\pi} + (\theta_1 - \theta_0)^2 + 2Y^{a=1}\theta_0 + 2Y^{a=0}\theta_1 - 2Y^{a=1}\theta_1 - 2Y^{a=0}\theta_0 \Bigr] \\ &=& \mathbb{E} \Bigl[ \frac{{Y^{a=1}}^{2}}{\pi} + \frac{{Y^{a=0}}^{2}}{1-\pi} - (\theta_1 - \theta_0)^2 \Bigr] \tag{12} \end{eqnarray}

Hajek estimator と Horvitz-Thompson estimator の漸近分散の比較

最後に、式 (3.10改) と式 (12) を比較してみます。

\begin{eqnarray} Var_{HT}(\hat{\theta}_{1} - \hat{\theta}_{0}) - Var_{Hajek}(\hat{\theta}_{1} - \hat{\theta}_{0}) &=& \mathbb{E} \Bigl[ \frac{{Y^{a=1}}^{2}}{\pi} + \frac{{Y^{a=0}}^{2}}{1-\pi} - (\theta_1 - \theta_0)^2 \Bigr] - \Bigl(\mathbb{E}\Bigl[ \frac{(Y^{a=1}-\theta_1)^2}{\pi} + \frac{(Y^{a=0}-\theta_0)^2}{1-\pi} \Bigr] \Bigr)\\ &=& \mathbb{E}\Bigl[ \frac{{Y^{a=1}}^{2}}{\pi} - \frac{(Y^{a=1}-\theta_1)^2}{\pi} \Bigr] + \mathbb{E}\Bigl[ \frac{{Y^{a=0}}^{2}}{1-\pi} - \frac{(Y^{a=0}-\theta_0)^2}{1-\pi} \Bigr] - (\theta_1 - \theta_0)^2 \tag{13} \end{eqnarray}

のように計算できます。

2次のモーメント ≧ 分散 の不等式を使うことにより、期待値が残っている部分については非負であることが示せるものの、最後の定数項がゼロ以下であるので、何かしらの不等式を示すにはもう少し展開する必要がありそうです。

ここで、さらに2次モーメントと分散の関係を明示的に利用することで、以下のように式を展開できると考えられます:

\begin{eqnarray} (13) &=& \mathbb{E}\Bigl[ \frac{\theta_1^{2}}{\pi} \Bigr] + \mathbb{E}\Bigl[ \frac{\theta_0^{2}}{1-\pi} \Bigr] - (\theta_1 - \theta_0)^2 \ \ \ (\because \text{Relationship between variance and moment}) \\ &>& \theta_1^2 + \theta_0^2 - (\theta_1 - \theta_0)^2 \ \ \ (\because 0 < \pi < 1) \\ &=& 2\theta_1\theta_0 \tag{14} \end{eqnarray}

よって、 Potential outcome の符号が等しければ右辺が0以上になり、「Horvitz-Thompson estimator より Hajek estimator の方が漸近分散が小さい」という主張が成立すると考えられます。

なお、ここで行った解析は専門家のレビューを通過しているわけではないため、以下の可能性が残っているということを示唆して終わりとします:

• 漸近分散の解析を改良すればより少ない仮定で不等式が示せる
• 漸近分散以外の観点で見ればより意味のある不等式が示せる
• 私の解析が誤っていて、このような不等式は示せない

終わりに

本記事では、因果推論で重要となる IPW について、いくつかの代表的な推定量の性質を解析しました。

一通り読んでいただいた方には、因果推論の著名な本の中でさらっと述べられていて盲目的に信じてしまいがちなところに意外な落とし穴があるかもしれない、ということをお伝えできたかと思います。

今後は、他の推定量の性質や、実験による主張の検証についての記事を書く予定をしています。

記事を読んでいただいた皆さんの豊かな因果推論ライフに貢献できれば幸いです。

*1:本記事では、6つの仮定 + Potential outcome の符号一致仮定のもとでの主張を考えます

*2:個人的には を処置 の Potential outcomeと書いた方が他の表記との整合性があると思うのですが、ここらへんは慣習の問題もありそうなのでよしなに読み替えてください

*3:Estimand の定義については what if book の10章をご覧ください

*4:「周辺平均の不偏推定量になる」と書かれています

*5:私の理解に誤りがあればご指摘いただけると幸いです

*6:多くの場合、傾向スコアモデルのパラメータ推定を含んだ形での解析を行う必要があるのですが、複雑になるので割愛しました。ご容赦ください

*7:漸近分散の解析に慣れていないので誤りがあればご指摘いただけると幸いです

調査観察データの統計科学 3.1の行間メモ

星野先生によって書かれた調査観察データの統計科学（以下本書と呼ぶ）を読んでいて、3.1節の式の導出に少し困ったので、メモを残します。 脚注に書いてある通り、厳密な証明は見つけられていないので、知見のある方はご指摘いただけると幸いです。

調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

目次はこちらです:

3.1節の概要と本記事の目的

3.1節でやることはざっくり以下の三つです。

1. バランシングスコアという概念を定義
2. バランシングスコアを条件づけたときに、割り当て と潜在的な結果変数 が独立になることを示す
3. バランシングスコアの関数として傾向スコアを紹介

本記事では、バランシングスコアの必要条件についての記述、ならびにバランシングスコアを条件づけたときの独立性について補足します。

バランシングスコアの定義

共変量を とします。本書P60によると、バランシングスコア とは、

それを条件付けすることにより、共変量と割り当てが独立になるような「共変量の関数」である

と定義されています。数式で書くと

を満たす関数のことです。

バランシングスコアの必要条件

本書では、バランシングスコアの必要条件として、関数 を使って と表現できることを挙げています。 その根拠として、以下の式の三つ目の等号成立条件において、上記のような表現が成立することが必要だから、というように書かれています。

\begin{eqnarray} p(z=1\mid b(x)) &=& \int p(z=1\mid x, b(x)) p(x\mid b(x)) dx \\ &=& \mathbb{E}_{x\mid b(x)} [p(z=1\mid x) ] \\ &=& \mathbb{E}_{x\mid b(x), \ p(z=1\mid x)} [p(z=1\mid x) ] \\ &=& p(z=1\mid x)\\ &=& p(z=1\mid x, b(x)) \tag{3.2} \end{eqnarray}

結論としては、この説明は必要条件ではなく十分条件の証明になっていると考えられます *1。

本書の中では細かい式変形の行間が省略されているため、解説を試みます。

まず、一つ目の等式は、 で条件つき確率を周辺化しています。こちらは周辺化の定義から自明かと思います。

次に、二つ目の等式は、 が の関数であることにより、 を条件づけから外せることを利用します *2。

三つ目の等式は、もし と表現できたならば、二つ目の等式の条件づけを外したときと同じロジックで を条件づけから外せることを利用します。等式の左側が条件づけから外れた式で、右側が条件づけられた式になります。

四つ目の等式は、 を条件づけた状態で の期待値を取っているので、期待値を外せる（であろう）ことを利用します*3。

最後の五つ目の等式は、二つ目の等式を再び適用します。

三つ目の等式以外は条件つき確率の定義から計算できるであろうことがわかったので、三つ目の等式が十分条件になることがわかるかと思います。

バランシングスコアを条件づけたときの独立性

ここでは、バランシングスコアを条件づけたときに、割り当て と潜在的な結果変数 が独立になること、すなわち

を示します。

まず、先ほどと同様のやり方で、条件つき確率を展開します。

\begin{eqnarray} p(z=1\mid y_1, y_0, b(x)) &=& \int p(z=1\mid y_1, y_0, x, b(x)) p(x\mid y_1, y_0, b(x)) dx \\ &=& \mathbb{E}_{x\mid y_1, y_0, b(x)} [p(z=1\mid y_1, y_0, x)] \tag{a} \end{eqnarray}

一つ目の等式は周辺化、二つ目の等式は が の関数であることを利用しています。

ここで、本書の2.5節で定義された"強く無視できる割り当て"条件が成立している、すなわち

が成立していると仮定すると、 なので、

\begin{eqnarray} \mathbb{E}_{x\mid y_1, y_0, b(x)} [p(z=1\mid y_1, y_0, x)] &=& \mathbb{E}_{x\mid y_1, y_0, b(x)} [p(z=1\mid x)] \tag{b} \end{eqnarray}

となります。

最後に、バランシングスコアの必要条件 *4 によって と表現できること、ならびに で条件づけた確率で の期待値をとるとき期待値は外せる（であろう）ことを利用すると、

\begin{eqnarray} \mathbb{E}_{x\mid y_1, y_0, b(x)} [p(z=1\mid x)] &=& p(z=1\mid x) \tag{c} \end{eqnarray}

となり、式(3.2)と合わせると以下が言えます。

\begin{eqnarray} p(z=1\mid y_1, y_0, b(x)) &=& p(z=1\mid x) &=& p(z=1\mid b(x)) \tag{d} \end{eqnarray}

これにより、式(3.3) が成立します。

最後に

一つ一つ式を追っていくと、これ自明っぽいけど本当に自明なのかな？と迷うところがいくつか出てきて、鍛錬が足りないなぁと実感しました。 期待値の式変形を一つ一つ解説した日本語資料は簡単には見つけられなかったので、本記事が何かしらの助けになると幸いです。

*1:こちらのブログでも指摘がされています: 「調査観察データの統計科学」3.1章 傾向スコアの数式メモ（前半） - 木曜不足

*2:厳密に証明した資料は見つからなかったので正しさは保証できないですが、ここではこれが成り立つとします

*3:これも厳密な証明は見つけられていないです。適切な資料をご存知の方はご教授いただけると幸いです

*4:式(3.2)では十分条件しか示せなかったのですが、ここでは必要条件が成立することを認めるとします。気になる方は紹介されている論文を辿ってください