【論文】Treatment Effect Estimation with Data-Driven Variable Decomposition
0．Abstract
因果推論における根本的な問題の一つは、変数間に交絡がある場合の観察研究における治療効果の推定です。交絡の影響を制御する手法としては一般的に傾向スコアが用いられますが、この手法では全ての観察された変数を交絡因子として扱い、治療には影響を与えないがアウトカムを予測する調整変数を無視してしまいます。近年、調整変数は推定される治療効果の分散を低減する上で効果的であることが示されました。しかし、観察研究において、交絡因子と調整変数を自動的に分離する方法はまだ確立されていません。特に、ビッグデータ時代には高次元の変数を取り扱うことが一般的ですが、そのような場合の分離は特に困難です。本論文では、Data-Driven Variable Decomposition (D2VD) アルゴリズムを提案します。このアルゴリズムは、(1) データ駆動型のアプローチで交絡因子と調整変数を自動的に分離し、(2) 同時に高次元変数を持つ観察研究における治療効果を推定することができます。標準的な仮定の下で、我々は提案する D2VD アルゴリズムが合成データと実際のオンライン広告データセットの両方において、変数を自動的に正確に分離し、最先端の手法よりも高い精度で治療効果を推定し、より狭い信頼区間を得られることを実験的に示します。
1．Introduction
因果推論とは
因果推論とは、ある事象 (効果) が発生する条件に基づいて因果関係を結論付けるプロセスであり (Holland 1986)、説明的分析のための強力な統計モデリングツールです。因果推論における金字塔的手法は、ランダム化実験 (Lewis と Reiley 2009 年) です。A/B テストなどがこれにあたり、異なる処置がユニット (対象) にランダムに割り当てられます。しかし、完全なランダム化実験は通常非常にコストが高くつく (Kohavi と Longbotham 2011 年) 、場合によっては実施不可能 (Bottou et al. 2013 年) なことが多々あります。したがって、観察研究において治療効果を推論するための自動的な統計手法の開発が強く求められています。
観察研究における治療効果推定
これまでの研究では、(Rosenbaum と Rubin 1983 年) が、傾向スコア調整に基づく治療効果推定のための統計的枠組みを提案しました。この枠組みは、マッチング、層化、逆重み付け、傾向スコアによる回帰など、観察研究における因果関係の分析で広く用いられています (Austin 2011 年; Lunceford と Davidian 2004 年; Kuang et al. 2016 年)。逆確率重み付けは最も一般的に用いられる手法であり、周辺構造モデルと呼ばれる大規模な因果モデルの一種です (Hernán, Brumback, and Robins 2000 年; 2002 年)。(Bang と Robins 2005 年) は、逆傾向重み付けと回帰を組み合わせることで、二重に頑健な推定量を提案しました。これらの手法は、経済学 (Stuart 2010 年)、疫学 (Funk et al. 2011 年)、医療 (Reis et al. 2015 年)、社会科学 (Lechner 1999 年)、広告 (Sun et al. 2015 年) など様々な分野で広く用いられています。
これらの手法の核心は、交絡変数の影響を排除して、治療効果推定の精度を大幅に向上させることです。しかし、多くの研究では、傾向スコアを推定する際、観察されたすべての変数を交絡変数として扱っています。その結果、高次元変数のシナリオでは、一部の変数は交絡変数ではなく、図 1 の因果ダイアグラムで示したように アウトカムを予測する調整変数（図 1における Z）になります。調整変数を無視すると、推定される治療効果は不正確になり、分散が過大評価されてしまいます。
調整変数の重要性
近年、一部の研究者は調整変数の重要性を指摘しています。(Brookhart et al. 2006 年; VanderWeele と Shpitser 2011 年) は、調整変数を因果推論モデルに含めるべきだと主張しています。(Sauer et al. 2013 年) は、このような調整変数に対する条件付けはバイアスを取り除くために必須ではなく、治療効果推定の分散を低減することができると示唆しています。ランダム化実験の設定では、(Bloniarz et al. 2016 年) は、Lasso を用いて調整変数を調整することで、推定される治療効果の分散を低減できることを証明しています。
観察研究におけるこれらの手法はすべて、因果構造 (つまり、変数が治療の原因なのか、アウトカムの原因なのか) が事前にわかっていることを前提としています。しかし、因果構造は、特にビッグデータ時代の高次元変数のシナリオでは、ほとんどの場合、事前知識だけでは十分に定義できません。観察研究において、交絡変数と調整変数を自動的に分離する方法はまだ確立されていません。
新しいアルゴリズムの提案
この問題に対処するために、我々はデータ駆動型変数分解 (D2VD) アルゴリズムを提案し、交絡変数の分離と平均治療効果 (ATE) の推定を同時に最適化します。具体的には、結合された直交性とスパース性正規化項を持つ正規化統合回帰モデルを提案します。このモデルは、同時に 1) データ駆動型のアプローチで交絡変数と調整変数を分離し、2) 過学習を防ぐために交絡変数でも調整変数でもない無関係な変数を排除し、3) 観察研究における ATE を推定することを可能にします。治療効果を推定する際、分離された交絡変数はそれらの処置効果への交絡影響を効果的に排除することができ、一方、調整変数はアウトカムの調整を通じて推定された平均治療効果の分散を著しく低減することができます。これにより、ベースライン手法よりも真の平均治療効果をより正確に、より狭い信頼区間で推定することが可能になります。
本論文の主要な貢献は以下の通りです。
- 観察研究における平均治療効果推定の精度と信頼区間にとって重要である、交絡変数と調整変数の自動分離という新しい問題に取り組みます。
- 正規化統合回帰モデルを提示することにより、交絡変数の分離と平均治療効果の推定を同時に可能にする、新規のデータ駆動型変数分解アルゴリズムを提案します。
- D2VD アルゴリズムの利点は、合成データと実世界データの両方で実証されています。また、ソーシャルマーケティング、ヘルスケア、公共政策など、他の因果推論研究にも容易に適用できます。
2. 調整済み平均治療効果（ATE）の推定
1) 記号と仮定
図 1 の因果ダイアグラムで説明したように、処置を乱変数 T、特定の処置 T = t に対応する潜在的アウトカムを Y(t) と定義します。この論文では、バイナリ処置 (t ∈ {0, 1}) のみ考えます。処置を受けたユニット (T = 1) を処置群、処置を受けていないユニット (T = 0) を対照群と定義します。次に、i = 1, 2, ..., m でインデックス付けされた各ユニットについて、処置 Ti、観察されたアウトカム Yi^obs、変数ベクトル Ui を観測します。ユニット i の観察されたアウトカム Yi^obs は以下の式で表すことができます。
Yi^obs = Yi(Ti) = Ti · Yi(1) + (1 − Ti) · Yi(0) (式 1)
観察研究における平均治療効果 (ATE) 推定のための 3 つの標準的な仮定 (Rosenbaum と Rubin, 1983) があります。
仮定 1: 安定単位処置効果 (Stable Unit Treatment Value Assumption) あるユニットの潜在的アウトカムの分布は、観察された変数を与えた上で、別のユニットに割り当てられた特定の処置の影響を受けないものと仮定します。
仮定 2: 非交絡性 (Unconfoundedness Assumption) 観察された変数を与えられた場合、処置の分布は潜在的アウトカムに独立であると仮定します。式で表すと、T⊥(Y(0), Y(1))|U です (⊥は独立を表す記号)。
仮定 3: 重なり (Overlap Assumption) すべてのユニットは、観察された変数を与えられた場合、いずれかの処置状況を受ける確率がゼロではないと仮定します。式で表すと、0 < p(T = 1|U) < 1 です。
2) 調整済み平均治療効果推定量
観察研究における因果推論の重要な目的は、アウトカム Y に対する平均治療効果 (ATE) を評価することです。ATE は、処置群と対照群のユニットの潜在的アウトカム間の平均 (差) を表します。式で表すと、ATE は以下のように定義されます。
ATE = E[Y(1) - Y(0)] (式 2)
ここで、Y(1) と Y(0) はそれぞれ、処置を受けた状態 (T = 1) のユニットと、処置を受けていない状態 (T = 0) のユニットの潜在的アウトカムを表します。E(·) は期待値を表す記号です。
式 (2) は実現不可能です。なぜなら、各ユニットについて、処置を受けた状態か対照群の状態か、対応する 1 つの潜在的アウトカムしか観測できないからです。これを「反実仮想問題」(Chan et al. 2010) と呼びます。
この反実仮想問題は、観察されていない潜在的アウトカムを近似させることで対処できます。最も簡単な方法は、処置群と対照群の平均アウトカムを直接比較することです。しかし、観察研究においては、処置の割り当てがランダムでない場合、交絡の影響が考慮されていないため、2 つのサンプルを直接比較するとバイアスが生じる可能性があります (Chan et al. 2010)。
観察研究において、ATE を偏らずに評価するためには、交絡変数の影響を制御する必要があります。仮定 (1, 2, 3) のもとで、(Rosenbaum と Rubin, 1983) は、交絡変数を制御するために必要な情報をまとめる傾向スコアを導入しました。傾向スコアは e(U) で表され、すべての変数 U が与えられたときのユニットが処置を受ける確率 (T = 1) として定義されました。実際には、交絡変数 X だけが処置と関連しているため、
e(U) = p(T = 1|U) = p(T = 1|X) = e(X) (式 3)
となります。傾向スコアに基づいて、Rosenbaum (1987) は、逆確率重み付け (IPW) 推定量 ATEIPW を用いて、式 (2) の反実仮想問題に対処するための変換アウトカム Y' を提案しました (Athey と Imbens, 2016 参照)。変換アウトカム Y' は以下のように定義されます。
Y' = Y^obs · (T - e(U)) / (e(U) · (1 - e(U)))
= Y^obs · (T - e(X)) / (e(X) · (1 - e(X))) (式 4)
また、IPW 推定量は以下のように定義されます。
ATEIPW = E[Y'] = E[ Y^obs · (T - e(X)) / (e(X) · (1 - e(X))) ] (式 5)
しかし、従来の傾向スコアに基づく多くの手法では、傾向スコアを推定する際に、通常、すべての観察された変数を交絡変数として扱います。観察された変数の中には、交絡変数ではなく、アウトカムに直接影響を与える非交絡変数 (adjustment variable) が存在する場合もあり、このような変数をすべて交絡変数として扱ってしまうと、推定される治療効果は不正確になり、分散が過大評価されてしまいます。
したがって、図 1 の因果ダイアグラムに基づき、観察されたすべての変数 U を、交絡変数 X、調整変数 Z、および無関係変数 I (図 1 からは省略) の 3 つの集合に分割することを提案します。そして、以下の仮定の下で、推定される平均治療効果の分散を低減するために、調整変数を組み込んだ新しい調整済み推定量を提案します。
仮定 4: 分離可能性 (Separateness Assumption) 観察された変数 U は、交絡変数 X、調整変数 Z、および無関係変数 I の 3 つの集合 (U = (X, Z, I)) に分割できます。ここで、X は交絡変数、Z は調整変数、I は無関係変数です。
仮定 4 を用いて、以下の定義に基づいて Y' から導かれる新しい調整済み変換アウトカム Y+ を導入します。
Y+ = (Y^obs - φ(Z)) · (T - e(X)) / (e(X) · (1 - e(X))) (式 6)
ここで、φ(Z) は Z に関連する Y' の分散を低減するのに役立ちます。
次に、以下の式で表される調整済み推定量 ATEadj を提案します。
ATEadj = E[Y+] = E[ (Y^obs - φ(Z)) · (T - e(X)) / (e(X) · (1 - e(X))) ] (式 7)
この調整済み推定量は以下の特性を持っています。
- 第一に、仮定 1 から 4 までを満たす場合、ATEadj は不偏推定量です。つまり、
E(Y+|X) = E(Y(1) - Y(0)|X) (式 8)
です。この特性は、因果ダイアグラムにおいて、条件付け集合 X が処置 T とアウトカム Y を結ぶすべてのバックドアパスを遮断し、かつ T の子孫変数を含まないため、パールのバックドア基準 (Pearl 2009) から明らかです。
- 第二に、ATEadj の漸近分散は、IPW 推定量 ATEIPW の分散以下です。
式 (5) と (7) を比較すると、IPW 推定量はすべての変数 U が与えられた場合、交絡変数 X だけを考慮していることがわかります。一方、私たちの推定量は、分散を低減するためにアウトカム Y に対する調整を行うために調整変数 Z を利用しています。同様な調整は、ランダム化実験における ATE 推定の分散を低減することができることが (Bloniarz et al. 2016) によって証明されています。
3. 自動変数分解 D2VD アルゴリズム
式 (8) を用いて、E(Y+) = E(Y(1) - Y(0)) を得ることができ、以下の目的関数を最小化することで、調整済み変換アウトカム Y+ を変数 U に対して回帰させることによって、推定される平均治療効果 (ATE) を求めることができます。
目的関数 min ||Y+ - h(U)||^2 (式 9)
この式を最小化することで、h(U) の期待値を推定して ATE を求めることができます。
実際的な計算においては、φ(Z) と h(U) を係数ベクトル α と γ を用いた線形関数として指定します。
φ(Z) = Zα, h(U) = Uγ (式 10)
また、係数ベクトル β を用いた線形ロジスティック回帰を使用して、傾向スコア e(X) を推定します。
e(X) = p(T = 1|X) = 1 / (1 + exp(-Xβ)) (式 11)
式 (10) と (11) の設定では、変数の分解 U = (X, Z, I) を既知であると仮定しました。しかし、実際には正確な分解はわかりません。したがって、観察されたすべての変数 U を使用して X と Z を置き換え、交絡変数と調整変数を自動的に分離するデータ駆動型のアプローチを提案します。
式 (9) の目的関数を以下のように更新します。
目的関数 min ||Y^obs - Uα||^2 ⊗ W(β) - Uγ||^2_2 (式 12)
制約条件として、
Σ_(i=1)^m log(1 + exp((1 - 2Ti) * Uiβ)) < τ,
||α||_1 ≤ λ,
||β||_1 ≤ δ,
||γ||_1 ≤ η,
||α⊗β||^2_2 = 0
W(β) := T - e(U) / (e(U) * (1 - e(U))) と Σ_(i=1)^m log(1 + exp((1 - 2Ti) * Uiβ)) は、傾向スコアを推定する際の損失関数を表します。 ⊗ はアダマール積を表します。 ||α⊗β||^2_2 = 0 という式により、係数ベクトル α は調整変数を分離するように最適化され、β は交絡変数を分離するように最適化されます。特に、α と β に対して直交正規化項を採用することで、交絡変数と調整変数の分離を保証します。さらに、α、β、γ に L1 ペナルティを課すことで、無関係な変数 I を排除して分散をさらに低減し、変数の疎ら性問題に対処します。
用語解説
- L1 ペナルティ: パラメータベクトルの要素の絶対値の和をペナルティとして加える手法
- アダマール積: 行列同士で対応する要素同士をかける積
- 疎ら性問題: 多くの変数が 0 になるようなモデル選択問題
これらにより、以下の最適化問題が得られます。これは J(α, β, γ) を最小化することと同義です。
式 (13): 最適化問題 J(α, β, γ) = f(α, β, γ) + g(α, β, γ)
ここで、
f(α, β, γ) = || (Y^obs - Uα) ⊗ W(β) - Uγ||^2_2 + μ ||α ⊗ β||^2_2 + τ * Σ_(i=1)^m log(1 + exp((1 - 2Ti) * Uiβ))
g(α, β, γ) = λ ||α||_1 + δ ||β||_1 + η ||γ||_1
です。
f(α, β, γ) は目的関数、g(α, β, γ) はペナルティ関数を表します。
近接勾配アルゴリズム (Parikh and Boyd 2013) の分割可能演算子 (operator splitting property) を用いると、t ステップ目における最適化されたパラメータ (α(t+1) など) を、関数 g(·) の近接作用素 prox_κ(t)g とステップサイズ κ(t) を用いて求めることができます。
式 (14): パラメータ更新式
α(t+1) = prox_κ(t)g [ α(t) - κ(t) ∂f(·) / ∂α ]
ここで、∂f(·) / ∂α は変数 α に関する関数 f(·) の勾配を表します。
式 (15): 近接作用素
prox_κ(t)g(x)
= {
xi - κ(t) * λ, ( xi ≥ κ(t) * λ )
0, ( |xi| ≤ κ(t) * λ )
xi + κ(t) * λ, ( xi ≤ -κ(t) * λ)
}
式 (15) の λ は、関数 g(·) におけるパラメータ α の係数です。β を最適化するのであれば、係数は δ、γ を最適化するのであれば η を用います。
近接勾配アルゴリズムを用いて、式 (13) の目的関数を最小化することができます。つまり、α、β、γ をランダムに初期化し、互いに他の 2 つのパラメータを固定したまま交互に更新し、収束するまでステップを繰り返します。データ駆動型変数分解アルゴリズムは、アルゴリズム 1 に示されています。
アルゴリズム 1 の各反復処理では、アルゴリズム 2 で説明されているように、OPTIMIZATION によってパラメータ α、β、γ を更新します。ここで、関数 f^κ(·) は以下のように定義されます。
式 (16): f^κ(x, y) = f(y) + (x - y) * [∂f(·) / ∂x]^T + (1 / (2κ)) * ||x - y||^2
また、関数 f(α, β, γ) に関する変数 (α, β, γ) に対する勾配は以下の通りです。
勾配の計算式:
∂f(·) / ∂α = -2 * (W(β) * 1^T * U)^T * R + 2μ * α * ||β||^2 ∂f(·) / ∂β = 2 * [(Y - Uα) * 1^T * (∂W(β) / ∂β)^T * R + U^T (T - exp(Uβ)) + 2μ * α * ||β||^2] ∂f(·) / ∂γ = -2 * U^T * R
ここで、
R = (Y - Uα) * W(β) - Uγ
であり、
∂W(β) / ∂β = (2T - 1) * exp((1 - 2T) * Uβ) * (1 - 2T) * 1^T * U
です。
アルゴリズム 1 によって最適化されたパラメータ ˆα、ˆβ、ˆγ を用いて、交絡変数 X = {Ui : ˆβi = 0}、調整変数 Z = {Ui : ˆαi = 0} と分離し、平均治療効果 (ATE) を ˆATED2VD = E(Uˆγ) として推定することができます。
1) D2VDアルゴリズムの計算複雑度
D2VDアルゴリズムの計算複雑度を考える上で最も重要なのは、関数f(α, β, γ)に関して変数α、β、γの勾配を計算する手順です。∂f(·)/∂α、∂f(·)/∂β、∂f(·)/∂γの計算複雑度はいずれもO(mn)です。ここで、mはサンプル数（データポイント数）、nは観察される変数の総数です。
アルゴリズム内のforループとwhileループは、一定時間 (= 定数時間) で実行される操作のみで構成されているため、これらのループにかかる時間は全体的な計算複雑度に大きく影響しません。
したがって、支配的な要素である勾配計算と、一定時間のかかる操作のみを考慮すると、D2VDアルゴリズムの全体的な計算複雑度はO(mn)であると結論付けられます。つまり、データポイント数や変数の数が増えるにつれて、アルゴリズムに必要な時間も比例的に増加します。
2) 観察研究における平均治療効果推定のパラメータチューニング
観察研究における平均治療効果 (ATE) 推定のパラメータチューニングの主な課題は、実際には真の ATE についての基底真実は存在しないことです。
この課題に対処するために、(Athey and Imbens 2016) と同様に、マッチング法を使用して ATE を評価し、「近似的な基底真実」として設定しました。具体的には、各ユニット i について、反対の処置ステータスを持つユニットの中で最も近い一致を見つけます。
match(i) = arg min_{j:Tj=1−Ti} ||Ui - Uj||^2_2 (式 17)
マッチングを厳密なマッチングに近似させるために、match(i) が一定の閾値を超えた場合はユニット i を除外します。マッチングされた処置群と対照群の平均アウトカムを比較することで、マッチング推定量で ATE を推定し、これを「近似的な基礎的真値」として ATEmatching で表します。
この「近似的な基礎的真値」を用いて、交差検証でアルゴリズムのパラメータをチューニングすることができます。
5. 結論
本論文では、観察研究において、より正確に、より信頼区間の幅を狭めて平均治療効果を評価する方法に焦点を当てました。従来のプロペンシティスコアに基づく多くの因果推論手法は、通常すべての変数を交絡変数として扱うため、不十分であることを論じました。我々は因果ダイアグラムに基づき、観察されたすべての変数から交絡変数と調整変数を分離することを提案しました。そして、変数分解と平均治療効果推定を同時に最適化するためのデータ駆動型変数分解 (D2VD) アルゴリズムを提案しました。合成データと実世界データでの実験結果は、我々のモデルの有用性と、D2VD アルゴリズムが観察研究における平均治療効果推定に有効であることを実証しています。
コメント
コメントを投稿