お久しぶりです。冬休みなので、以前Twitterでとったアンケートで一番人気だった内容について書きます。
今回は統計“モデル”を使うことの意味について因果推論の視点からまとめてみようと思います。普段なんとなく回帰分析を使っている人は一読をおススメします。なぜモデルを使うのか、（多くの場合無自覚に）置かれている仮定は何かなどを中心に書いていきます。
なお、本記事の内容はハーバード公衆衛生大学院でMiguel Hernan氏が教える中級因果推論の授業で教えられる内容を基にしています。来学期はこの授業のTAをしますので質問等あればご気軽に。
統計学をきちんと勉強した方であれば基本的な内容になりますので悪しからず、、、
因果推論における「推定(Estimation)」とは
統計的にみることができる因果関係、"平均因果効果"については過去のブログ記事で説明しました。今回はこの定義を知っていることを前提に話を進めていくので、必要あれば以下の記事を読み返してください。
端的にいえば、
反事実モデルのもと定義される平均因果効果（因果関係）はある一定の仮定のもとで条件付き期待値の差（もしくは比）と一致する
という内容でした。
ポイントは反事実モデルで定義される「因果関係」は現実世界で観測不可能なのに対して、「条件付き期待値」はデータから実際に計算可能ということです。
この、「ある一定の仮定」（疫学の用語でExchangeabilityという仮定などを含みます。総じて”Identifiability”と呼ばれることもあります)というのが成立しているかどうかが統計的因果推論の肝なわけですが、今回はとりあえずこれらの仮定が成立しているとしましょう。
すると、ランダム化比較試験のようにMarginalな（みている集団全体に対して）Exchangeabilityが期待できる場合は
となります。つまり、「A=1のひとたちにおけるYの期待値とA=0の人たちにおけるYの期待値の差が、Aを０から１に変えることがYに与える平均因果効果と一致する」ということです。
なお、平均因果効果を比で定義した場合は条件付き期待値の比を計算することになります。今回はA=1または０となる二値のtreatmentを考えていますが、「A=a0 or a1とすることでAがa1-a0変化したときの因果効果」と一般化することもできます。
観察研究では、あらゆる交絡因子C（DAGを使って考える）を条件づけることでConditional Exchangeabilityを目指します。その仮定のもとで、
となります。要は条件付き期待値を計算すればいいという点は同じです。
左側の因果効果も条件付きになっているのですが、この点については後ほど。
ここまでのプロセスをIdentificationといいますが、ではこの条件付き期待値を実際にどうやって求めるのか。ここからが"推定"（Estimation）と呼ばれるプロセスです。
条件付き期待値の推定にはいろいろな方法があるのですが、統計モデルを使った回帰分析はそのうちの一つです。*1
ノンパラメトリックvsパラメトリック
超シンプルな回帰”モデル”
なんとなく使われることの多い回帰分析ですが、条件付き期待値をモデルしているものと捉えることができます。「モデル」とは数式を使って共変量がある値をとったときのアウトカムYの期待値がどのように決まるかを記述したものです。例えば、次のシンプルな回帰”モデル”を考えます。
左側が条件付き期待値になっていることに注目してください。Aが二値変数だとすると、この回帰“モデル”のβ0と β1を（OLSや最尤法などを使って）推定することで
E[Y|A=0] = β0 + β1*0 = β0
E[Y|A=1] = β0 + β1*1 = β0 + β1
といった具合に条件付き期待値の推定値を得ることができるのです。
（Marginal）Exchangeabilityの仮定のもとで
E[Y|A=1] - E[Y|A=0] = β0 + β1 - β0 = β1
と、AがYに与える平均因果効果が一致するわけです。今回は単純な線形回帰分析（linear regression)でしたが、その他の回帰分析（例：ロジスティック、ポアソンなど）でも左側がlogit pやlog Yになったりするだけで、本質的な違いはありません。
さて、このシンプルな例ですがわざわざ回帰”モデル”を使わずとも条件付き期待値を推定することができます。
層別分析で条件付き期待値を計算
そもそも条件付き期待値の意味するところを考えれば非常に簡単です。
例えば、E[Y|A=0]は「A=0であったときのYの期待値」「A=0の人たちにおけるYの期待値」という意味なので、
A=0の人たち、A=1の人たちそれぞれで期待値を推定してあげればいいだけです。このように各Aの値をサブグループごとに期待値を推定することを層別分析といいます。何のひねりもありません。
実は今回のケースだと、層別分析と回帰分析、それぞれから得られる条件付き期待値の推定値は完全に一致します。
実際にRで簡単なシミュレーションをしてみました。
gist50530c2663f95f614b37282447c320b1
A=0のときのYの期待値（E[Y|A=0] = 20)とA=1のときのYの期待値(E[Y|A=1] = 15)を層別分析・回帰分析を使って推定した結果が次です。
回帰分析の結果を使うと、
E[Y|A=0] = β0 + β1*0 = β0 = 20.953
E[Y|A=1] = β0 + β1*1 = β0 + β1 = 20.953 - 5.985 = 14.968
となり、層別分析の結果と一致しています。
結果を図示すると以下のようなイメージです。赤点が層別分析から得られた各グループにおけるYの期待値の推定値。
ノンパラメトリックな“モデル”
なぜ層別分析と回帰分析の結果が一致するのか。それは今回用いた回帰モデルがSaturated Modelと呼ばれるものだからです。
Saturated modelとは、共変量（モデルの右側に含まれている変数）の値の組み合わせのパターンの数と回帰モデルのなかに含まれているパラメーター（β）の数が一致しているモデルのことです。
先ほどの例では、モデルの右側には二値変数Aしか含まれていません。つまり、考えうる共変量の組み合わせはA=0かA=1の2パターン。モデルにはβ0とβ1の二つパラメーターが入っているのでこれはSaturated modelです。Saturated Modelでは数学的な制約なく、注目しているすべての条件付き期待値をユニークに推定することができるので層別分析と結果が一致します。
別の例を見てみましょう。二値変数Aのほかに、新たな二値変数Cを調整した条件付き期待値、E[Y|A, C]を考えます。モデルの右側には二値変数が二個あるので、考えうる共変量の値の組み合わせは
① A=0, C=0, ② A=0, C=1, ③ A=1, C=0, ④ A=1, C=1
の4パターンです。よって、対応するSaturated Modelは４つのパラメーターをもつ
になります。この“モデル”から推定される条件付き期待値は①～④のサブグループそれぞれで期待値を推定したものと同じ結果になります。
層別分析のようにモデル（アウトカムYと共変量の間の関係に一定の制約を想定する）を使わないアプローチをノンパラメトリックといいます。それに対してモデルを使ったアプローチの多くはパラメトリックと呼ばれます。
Saturated Modelは“モデル”と名前にあるだけに、パラメトリックなアプローチのように思えます。しかし、Yの条件付き期待値と共変量の関係性になんの制約も置いていない（全ての共変量のパターンに対して、Yの条件付き期待値は自由な値をとれる）ので本質的にはノンパラメトリックといえます。
なぜパラメトリックモデルを使うのか
Saturated Modelよりもパラメータの数が少なくなったモデルは、Unsaturated Modelと呼ばれます。Unsaturated Modelは条件付き期待値が取りうる値に対して数学的な制約を置くのでパラメトリックなアプローチです。例えば、次のようなモデル。
このとき、共変量の組み合わせ全パターンに対応する条件付き期待値は
E[Y|A=0, C=0] = β0
E[Y|A=0, C=1] = β0 + β2
E[Y|A=1, C=0] = β0 + β1
E[Y|A=1, C=1] = β0 + β1 + β2
となります。どのような制約があるかわかりますか？
実は、最初の三つの条件付き期待値が決まった時点で β0, β1, β2が確定するので、自動的に最後の期待値（E[Y|A=1, C=1])の値も確定してしまいます。自由な値がとれないわけです。層別分析をイメージすると、A=1,C=1の人たちにおけるYの期待値が他のサブグループにおける期待値によって勝手に決まっていることになります。
このような制約がある推定方法はできることならば避けたいですよね。ではなぜ層別分析やSaturated Modelのようなノンパラメトリックな方法ではなく、制約の多いパラメトリックなアプローチをとるのか。パラメトリックモデルが必要になる状況は大きくわけて二つあります。どちらの場合も共通しているのはCurse of Dimentionalityという問題です。
パラメトリックモデルが必要な状況①：共変量が多い
観察データからの因果推論を行う際、操作変数や自然実験など経済学的なアプローチが使える場合を除いて、あらゆる交絡因子を調整することでConditional Exchangeabilityを目指していくというのは前述のとおりです。
よくある臨床研究の論文などで調整されている交絡因子の数は、２つや３つではありません。例えば、二値変数の交絡因子を10個回帰分析で調整したとしましょう。その場合、考えうる共変量の組み合わせパターンは2^10=1024通りです。
例えば層別分析をしようとすると、サンプルサイズが十分に大きくないと、1024通りのうちの何パターンはサブグループ内に一人もデータが存在しないなんてことがあり得ます。回帰モデルを使うにしても1024個のβを推定する必要がでてきます。
さらに共変量の数が増えていったり、二値変数だけでなく三種類以上の値を取りうるカテゴリ変数が含まれてくると組み合わせパターンはもっと増えていきます。このような状況は統計学でCurse of Dimentionality（次元の呪い）と呼ばれています。
このような理由から、調整している変数が複数ある場合ではノンパラメトリックなアプローチが厳しく、条件付き期待値のふるまいに対して仮定を置くパラメトリックモデルに頼らざるをえなくなることがあります。
パラメトリックモデルが必要な状況②：連続変数を扱っている
たとえ調整している変数が一つしかなかったとしても、その変数が連続変数であった場合はパラメトリックモデルが必要になります。
例えば、E[Y|A]のAが連続変数である場合。Aがとりうるすべての値に対して層別分析は不可能だというのは簡単に想像できると思います。A=0.001の人たち、A=0.002の人たち、、、といった具合にサブグループを作っていくととんでもない数になってしまいます。同じようにSaturated Modelを作ろうとすると、とてつもない数のパラメータβを推定する必要がでてきます。
因果推論を目的としたパラメトリックモデルを使った回帰分析の仮定
回帰分析をつかって、「〇〇の効果を推定しました」という論文は多いですがどのような仮定に基づいているのでしょうか。
仮定１：Conditional Exchangeabilityの成立
重要なのでもう一度いいます。Conditional Exchangeabilityが成立していないのであれば、条件付き期待値をつかって見られる関連は因果関係と一致しません。
言い換えると、モデルに含まれている共変量Cを条件づけることでDAG上のすべてのBackdoor pathを閉じることができるという仮定です。
DAGについては過去記事を参照ください。
もちろん、このような強い仮定がperfectに成立することは現実的に難しそうです。そのような場合は、未測定の交絡因子Cがどの程度結果に影響を与えそうなのか検討をつける、（交絡バイアスのための）感度分析といった考え方があります。
注意していただきたいのは、この仮定はIdentifiabilityに関係するものなので、推定の段階でノンパラメトリックであろうがパラメトリックであろうが必要な条件であるということです。
この条件下では
となります。仮定のもとでデータから見ることができる平均因果効果がCによって条件づけられていることがポイントです。つまり、「Cが同じ値をとっている集団におけるAの平均因果効果」を推定しています。これをConditional Average Treatment Effectと呼びます。当然ながら、集団全体に介入を行ったときの因果効果(Marginal Average Treatment Effect)とは必ずしも一致しません。
仮定２：掛け算項の有無
たとえば二値変数Cだけ調整すれば仮定１が成立するとします（現実にはそんなことはありえませんが）。この場合、ノンパラメトリックなモデルは
パラメトリックなモデルを使うと、
のような回帰分析を行うことになります。パラメトリックモデルではA*Cの掛け算項がありません。これは以前も紹介した（AとCの間の）「交互作用（厳密には効果修飾）」を表す項です。Aの効果がCのレベルによってどの程度異なるのかを表す項です。
掛け算項がないパラメトリックモデルでは、Aの効果はCの値によらず一定であるおいう仮定がおかれています。
ちなみにこのモデルが正しいのであれば、推定されているConditional Effect (E[Y_a=1|C] - E[Y_a=0|C])とMarginal Effectは一致するはずです。Cの値と効果の大きさが無関係なのですから。
しかし、効果の異質性(heterogeneity)が存在すると考えられる場合は掛け算項を含めたモデルを使わなければ、モデルの誤設定（misspecification)が生じていることになります。
さらに分析の目的がConditional EffectではなくMarginal Effectの推定にある場合は、別途そのための方法を使う必要があります。疫学の世界ではparametric g-formulaやMarginal Structural Modelといったたぐいのものが使われます。また近々解説したいと思います。
用いている手法がConditional Effectを推定しているのか、Marginal Effectを推定しているのか、Constant Effectの仮定のもとでConditional Effectを推定することでMarginal Effectも見れていると主張しているのか。これらを区別し、自分の分析が何を推定しているのかについて自覚的であることは個人的にとても重要だと思っています。
仮定３：連続変数とアウトカムの関係性
モデルに連続変数が含まれていたとします。さきほど書いたように連続変数に対してノンパラメトリックなアプローチはほとんど無意味です。次の２つのパラメトリックモデルを考えましょう。
いずれも連続変数AとYの条件付き期待値の関係性に一定の制約（functional form)を置いています。
例えば一つ目のモデルでは、E[Y|A]とAの間にlinearity（直線関係）が存在すると仮定されています。言い換えると、Aが一単位増加したときの効果が一定という仮定です。A=0からA=1になるのと、A=100からA=101になるのも同じ効果ということです。
二つ目のモデルもアウトカムとAの関係が二次関数的(quadratic)に記述されるという仮定を置いています。
（仮定４：同一カテゴリ内ではアウトカムと一定の関係）
カテゴリ変数を使ったときの仮定です。例えば「０＝非喫煙者、１＝喫煙者」といったカテゴリ変数や、年収について「1＝〇〇万円から△△万円、２＝◇◇万円から・・・」といった具合に連続変数をカテゴリにする場合があります。
このモデルでは、カテゴリ変数Cで測定されているものとアウトカムの関係が、同じCの値をとる人たちの間で一定という仮定をおいています。
例えば「喫煙」とアウトカムの関係は「喫煙者」の間では一定ということです。喫煙者であれば一日に何箱吸ってようがアウトカムとの関係性が同じという仮定です。年収の例も同様です。同じカテゴリ（例：年収500-800万円)内であったとしても、500万円の人と800万円の人ではアウトカムとの関係が違うと考えられるかもしれませんが、このカテゴリ内において年収とアウトカムの関係は一定と仮定されています。
このように粗い(Crude)カテゴリを使うと、強い仮定をおくことになります。Cを交絡因子として調整している場合は、バイアスを完全に除去できていない（Residual Confounding)の可能性がでてきます。
なお、これはモデルの仮定というよりは仮定１のConditional Exchangeabilityに関連するもの、すなわちIdentificationのための仮定になります。
仮定５：アウトカムと共変量の関係が加法的
モデルにおける各共変量のアウトカムへの貢献具合が加法的に決まるという仮定です。例えば、次のモデル
では条件付き期待値がAとCそれぞれのContribution (i.e. β1とβ2)の足し算で表現されています。ですが、もしかすると真のモデルは次のようなものかもしれません。
この場合E[Y|A]に対するAとCの貢献が加法的ではありません。実は両辺logをとると再び(logE[Y|A]に対する)additiveモデルになるのですが、それもまた仮定にすぎません。
まとめ
回帰分析をつかって因果推論をおこなう場合、おなじみの因果推論のための仮定に加えて、さらにモデルに対する仮定を多くおくことになります。
これらモデルの仮定を理解することは、よりadvancedなモデルをつかった因果推論（Marginal Structural Model、doubly robustなど）の意義を考えるうえでも重要です。
また強調しておきたいのは、ExchangeabilityのようなIdetificationに関わる仮定は人間の頭をつかって必要な条件を考えていく必要があるのに対して、Modeling Assumptionsについては機械学習のようなアプローチが有効な可能性があるということです。
近々ブログに書きたいと思っているのですが、UC Berkleyの研究者がごり押ししているSuperLearnerというアプローチもモデリングに機械学習のテクニックを取り入れることで単純なパラメトリック推定よりも良い推定ができるというアイデアです。
因果推論と機械学習の融合は、こういったモデルの仮定に対する突破口になるのかもしれないと個人的に期待しています。
*1:回帰モデルの係数を推定する方法にもいろいろ（OLS, 最尤法など）ありますが、今回はその違いには触れません
