はじめに

Pearl派の因果推論について勉強しています。勉強した内容の自分用のメモになります。いつも通り、少しづつ加筆していくスタイルで行きます。

また、理論面は(雰囲気で)なんとなくわかってきたのですが、実際分析を行っていくイメージが持てなかったので、手を動かそうと思いました。構造的因果推定に簡単にまとめるとともに、Rでそのように行うのか試してみました。

Assessing Causality from Observational Data using Pearl's Structural Causal Models | R-bloggers

入門統計的因果推論

作者:Judea Pearl,Madelyn Glymour,Nicholas P. Jewell
発売日: 2019/08/28
メディア: 単行本（ソフトカバー）

構造的因果モデルの基礎

作者:学, 黒木
発売日: 2017/08/24
メディア: 単行本

統計的因果推論―回帰分析の新しい枠組み (シリーズ・予測と発見の科学)

作者:宮川雅巳
発売日: 2004/04/01
メディア: 単行本

はじめに、上記のPearl先生の訳本に非常に素晴らしい言葉が載っていましたので引用します。

「なぜ因果を学ぶのか」という問いに対する答えは、「なぜ統計学を学ぶのか」という問いに対する答えと同じである。

至極うなずけるお言葉だなと感じます。

グラフィカルモデル

有向非巡回グラフ(DAG)

有向グラフは頂点と有向辺（方向を示す矢印付きの辺）からから構成されるグラフのこと。

このグラフにおいて、ある頂点から出発し、辺をたどり元の頂点に戻ってこない、巡回路を持たないグラフのことを有向非巡回グラフを呼ぶ。

因果ダイアグラム

因果ダイアグラムはデータの生成過程を視覚的に表現したもの。非巡回的有効グラフGとその頂点に対応する確率変数の集合、確率変数の集合が与えられる。グラフGが確率変数の関数関係を規定される。

$\displaystyle{ X_i=g_i(pa(X_i),\epsilon_i), i=1,\dots,p }$

確率変数がこの関数関係に従って自律的かつ定常的に生成され得る時、 $G$ を因果ダイアグラムと呼ぶ。

この時の同時分布は、因果ダイアグラム $G$ に従って、次のように逐次的に因数分解された形で表現できる。

$\displaystyle{ pr(x_1,...,x_p)=\prod_{i=1}^p pr(x_i|pa(x_i)) }$

因果ダイアグラムの例を次の図1に示す。

f:id:saltcooky:20200412013435p:plain — 図1. 因果ダイヤグラムの例

ベイジアンネットワーク

双方向線のない非巡回的有向グラフの同時分布をグラフ $G$ に従う逐次的因数分解の形、次のように規程できる時、ベイジアンネットワークという。

$\displaystyle{ pr(x_1,\dots ,x_p)=\prod^p_{i=1} pr(x_i|pa(x_i)) }$

$pr(x_i|pa(x_i))$ は $X_i$ の親の集合 $pa(X_i)=pa(x_i)$ を与えたときの $X_i$ の条件付き確率分布であり、 $pa(X_i)$ が空集合の時には $X_i$ の周辺分布 $pr(X_i)$ を意味する。

ベイジアンネットワークと因果ダイアグラム

ベイジアンネットワークと因果ダイアグラムでは次にような違いがある。

矢印があるということは、ベイジアンネットでは変数間に統計的従属関係が存在する可能性があること、因果ダイアグラムは変数間に何らかの条件付き独立関係が成り立っていることを意味する。

矢印がないということは、ベイジアンネットはそれらに直接的に因果関係が存在する可能性があること、因果ダイアグラムは変数間に直接的な因果関係が存在しないことを強調することを意味する。

d分離

ノードXからYへの道が複数ある場合がある。この時、XとYとは背反なノードの集合Zにより、XからYへの道を遮断する時にd分離（有効分離）しているという。

ZがXとYをd分離しているかの基準は、次のいづれかを満たす時である。

XとYを結ぶ経路の合流点でその合流点とその子孫がZに含まれないものがある。
XとYを結ぶ道に非合流点で、Zに含まれるものがある。

次のようなDAGを考える。

f:id:saltcooky:20200412170402p:plain

XとYを有効分離する集合は、 ${Z_1 },{Z_2 },{Z_3 },{Z_1,Z_2 }, {Z_2,Z_3 }, {Z_1,Z_2,Z_3 }$ である。

一方、 ${ Z_1,Z_3 }$ はXとYを有効分離しない。

合流点である ${ Z_3 }$ の固定してしまうと、 ${ Z_2 }$ の値からYが得られるようになり、有効分離出来ない。

${ Z_1 }$ の固定も同様。

そのため、 ${ Z_1,Z_3 }$ を固定してしまうと、 ${ Z_2 }$ を通じてXの値からYが得られるようになり、有効分離出来ない。

因果効果と識別方法

構造方程式モデル

数学的には、複数の変数間の従属関係を誤差を加えた連立一次方程式によって表現した統計モデル。

図1で示す構造方程式モデルの一例は次にようになる。

$\displaystyle{ Y=g_y(X,Z_1,Z_2,\epsilon_y) \\ X=g_x(Z_2,\epsilon_x) \\ Z_2=g_{z_2}(Z_1,\epsilon_{z_2}) \\ Z_1=g_{z_1}(\epsilon_{z_1}) \\ }$

構造方程式モデル単体では、因果関係を記述することができない。

因果ダイヤグラムはデータ生成過程に対する定性的な視覚表現であり、即座に構造方程式モデルを意味するものである

do演算子

構造方程式モデルにおいて、興味のある状態にモデルを変更すること考える。

図1において $X=x$ に固定に肯定したすると、 $Z_2$ から $X$ への矢印がなくなる。

$Z_2$ と $X$ は独立した状態に変化する。この操作を外的操作と呼ぶ。

この操作をを行った時に得られる確率変数の同時分布は次のようになる。

$\displaystyle{ pr(y,z_1,z_2|do(X=x))= pr(y|X=x,z_1,z_2) pr(z_2|z_1)pr(z_1) }$

ここで、 $do(X=x)$ は外的操作によりXを $X=x$ に固定する操作を意味しており、do演算子と呼ばれている。

構造方程式モデルでは $do(X=x)$ はXを右辺にもつ構造方程式の関数を $X=x$ という関数に置き換えることを意味する。

$X=x$ が観測されている状態と、外的操作により $X=x$ に固定する状態の違いは次のようになる。

$X=x$ が得られたときに $Z_1$ や $Z_2$ の確率分布はそれぞれ $pr(Z_1|x)$ と $pr(Z_2|x)$ となり、 $X=x$ 時に $Z_1$ や $Z_2$ の状態が推定できることを意味する。

一方で、外的操作により $X=x$ に固定する場合、 $Z_1$ や $Z_2$ の関する分布は次にようなる。

$X=x$ の外的操作により $Z_1$ や $Z_2$ の確率分布が変化しないことを意味する。

また、 $X=x_1$ や $X=x_0$ をランダムに割り当てていくような単純ランダム化に対応する。

因果効果

因果ダイヤグラムにおける頂点集合 $V={X,Y} Z$ とする時の、 $X=x$ から $Y=y$ への因果効果は次のように定義される。

$\displaystyle{ pr(y|do(X=x))=\sum_{z} \frac{pr(x,y,z)}{ pr(x|pa(x)))} }$

図1において、 $Z_2$ と $X$ が $Z_1$ によって条件付き独立になっている場合に、次のように変形できる。

$\displaystyle{ pr(y|do(X=x))=\sum_{z_1,z_2} pr(y|X=x,z_1,z_2)pr(z_1,z_2) }$

この式を調整化式と呼ぶ。

$Z_2$ を条件付けることにより $Z_1$ と $X$ が独立になることを用いると、次のように変形できる。

$\begin{aligned} pr(y|do(X=x)) &=\sum_{z_1,z_2} pr(y|X=x,z_1,z_2)pr(z_1,z_2) \\ &=\sum_{z_1,z_2} pr(y|X=x,z_1,z_2)pr(z_1|X=x, z_2) pr(z_2)\\ &=\sum_{z_2} pr(y|X=x,z_2)pr(z_2) \\ \end{aligned}$

因果効果は $Z_1$ の値によらないことことがわかる。

その他の因果効果の指標

そのほか関連する指標は以下の通り。

$X=x$ からYへの平均に対する因果効果

$\displaystyle{ E(Y|do(X=x))=\sum_{y} y pr(y|do(X=x)) }$

$X=x$ と $X=x'$ を比較した時の $X$ から $Y$ への平均的因果効果

$\displaystyle{ E(Y|do(X=x))-E(Y|do(X=x')) }$

$X=x$ と $X=x'$ を比較した時のリスク比

$\displaystyle{ \frac{E(y|do(X=x))}{E(y|do(X=x'))} }$

識別可能

特定の頂点 $X$ から $Y$ への因果効果を共変量調整などにより一意に算出できることを因果効果が識別可能であると呼ぶ。

推論規則

因果ダイアグラムで成立している独立性・条件付き独立性を利用した式変形を行うことで、因果効果についての基本的性質を得ることができる規則。

規則１

因果ダイアグラムGによりXへ向かう矢線を全て除いたグラフにおいて、XがYとZを有効分離するならば、次が成り立つ。

$\displaystyle{ pr(y|do(X=x),z)=pr(y|do(X=x)) }$

規則２

因果ダイアグラムGによりXから出る矢線を全て除いたグラフにおいて、ZがXとYを有効分離するならば、次が成り立つ。

$\displaystyle{ pr(y|do(X=x),z)=pr(y|x,z) }$

また、Zが空の時次が成り立つ。

$\displaystyle{ pr(y|do(X=x))=pr(y|x) }$

規則３

因果ダイアグラムGによりXへ向かう矢線を全て除いたグラフにおいて、ZがXとYを有効分離するならば、次が成り立つ。

$\displaystyle{ pr(y|do(X=x),z)=pr(y|z) }$

また、Zが空の時次が成り立つ。

$\displaystyle{ pr(y|do(X=x))=pr(y) }$

バックドア基準

XからYへの因果効果を考える際には、交絡が存在するかを考え、調整する必要がある。

この時にどの共変量に対して調整を行うかを判断する基準としてバックドア基準がある。

定義は以下の通り。

<定義> バックドア基準

非巡回的有効グラフ $G$ において $X$ は $Y$ の非子孫であるとする。この時、次の条件を満たす頂点集合が、バックドア基準を満たす。

$X$ から $Z$ の任意の要素への有効道がない
$G$ より $X$ から出る全ての矢印を除いたグラフにおいて、 $Z$ が $X$ と $Y$ を有効分離する

具体的に、図1においるXからYへの因果効果を識別することを考える。

XとYを有効分離する集合は{ $Z_2$ },{ $Z_1,Z_2$ }である。

また、いずれもXの非子孫であるためバックドア基準を満たす。

そのため、{ $Z_2$ }または{ $Z_1,Z_2$ }を調整することにより、XからYへの因果効果を推定することができる。

フロントドア基準

次の図2のようなXからYの有向道上に中間変量Sが存在する際、XからYへの因果効果を識別するのに必要な基準をフロントドア基準と呼ぶ。

f:id:saltcooky:20200326015312p:plain

定義は以下の通り。

<定義> フロントドア基準

非巡回的有効グラフ $G$ において $X$ は $Y$ の非子孫であるとする。この時、次の条件を満たす頂点集合が、フロントドア基準を満たす。

$X$ から $Y$ への任意の有効道上に $S$ の要素が存在する
$G$ より $X$ から出る全ての矢印を除いたグラフにおいて、空集合は $X$ が $S$ を有効分離する
$G$ より $S$ の任意の要素から出る全ての矢印を除いたグラフにおいて、 $X$ は $Y$ が $S$ を有効分離する

なお、頂点集合 $S$ に該当する頂点が一つのみである場合、上記の条件2と条件3は次のように、言い換えることができる。

　2'. $G$ において、空集合は $(X,S)$ についてバックドア基準を満たす

　3'. $G$ において、 $S$ は $(S,Y)$ についてバックドア基準を満たす

これらの条件からフロントドア基準は、意味合いとしてバックドア基準の二段階適用と言うことがわかる。

この時の因果効果は次にようにかける。

$\displaystyle{ pr(y|do(X=x))=\sum_{x',s}pr(y|x',s)pr(s|X=x)pr(x') }$

ここで、 $x'$ は $X$ の取りうる値である。

次のようにも変形することができる。

$\displaystyle{ pr(y|do(X=x))=\sum_{s} pr(s|X=x) \left\{ \sum_{x'} pr(y|x',S=s)pr(x') \right\} }$

ここで、条件2'により $pr(s|X=x)$ について次のことが得られる。

$\displaystyle{ pr(s|do(X=x))=pr(s|X=x) }$

また、条件3'から次のことが言える

$\displaystyle{ pr(y|do(S=s))= \sum_{x'} pr(y|x',S=s)pr(x') }$

ここからも、バックドア基準の二段階適用であることがわかる。

条件付き因果効果

ここまでは外的操作を $X=x$ としていたが、外的操作をどのようなものにするかを共変量 $Z$ に依存する場合が考えられる。図1の因果グラフが考えられる場合の条件付き因果効果をみる。

f:id:saltcooky:20200326014907p:plain

この時、 $X$ は $Z_1$ により決定される外的操作を次のように示す。

$\displaystyle{ X=h(Z) }$

この時の因果効果は次にように表現できる。

$\displaystyle{ pr(y|do(X=x))=\sum_{s} pr(s|X=x) \left\{ \sum_{x'} pr(y|x',S=s)pr(x') \right\} }$

これは次のような因果グラフへ変化させる操作であるといえる。

f:id:saltcooky:20200326014925p:plain

同時介入効果

変数ベクトル $\boldsymbol{X}$ を定数ベクトル $\boldsymbol{x}$ に固定した時の $Y$ の分布が同時介入効果という。

$\displaystyle{ pr(y|do(X=x))=\sum_{z} \frac{pr(\boldsymbol{x}, y, z)}{ \prod_{x \in \boldsymbol{x}} pr(x|pa(x)) } }$

この時の同時因果効果の識別可能性の判断には認容性基準に従って行う。

構造方程式モデルによる因果効果の識別

図1の因果ダイアグラムの線形構造方程式モデルを次のように考える。

$\displaystyle{ Y=\alpha_y+\alpha_{yx}X+\alpha_{yz_2}Z_1+\alpha_{yz_2}Z_2+\epsilon_y \\ X=\alpha_x + \alpha_{xz_2}Z_2+\epsilon_x \\ Z_2=\alpha_{z_2}+\alpha_{z_2z_1}Z_1+\epsilon_{z_2} \\ Z_1=\alpha_{z_1}+\epsilon_{z_1} \\ }$

ここで、各 $\epsilon$ は互いに独立である。

<あとで>

構造因果モデルから潜在反応モデルへ

<あとで>

Rで線形回帰モデルを用いた因果推論

バックドアパスの識別

因果ダイヤグラムの生成とバックドアパスの識別にはdagittyパッケージを用いる。

今回考える因果ダイヤグラムは次のようなものである。

library(tidyverse)
library(dagitty)
library(lavaan)

graph <- dagitty('dag {
             S [pos="0,0"]
             T [pos="1,2"]
             U [pos="0,4"]
             V [pos="2,1"]
             W [pos="2,3"]
             X [pos="3,0"]
             Y [pos="3,4"]
             Z [pos="4,2"]

  S -> X -> Z
  S -> T
  V  -> W
  T -> V -> X -> Z
  T -> V -> Z
  T -> W -> Z 
  U -> Y -> Z
  U -> T
  Y -> Z 
}')


plot(graph)

f:id:saltcooky:20200317005319p:plain

Xからyへのパスは、paths関数で列挙できる。

paths(graph, "X", "Y")

調整すべき頂点集合を探す。

type引数で"all"とすると調整すべきパスの集合の全てが表示される。

"minimal"とすると最小限の候補が表示される。

> adjustmentSets(graph, "X", "Y", type = "minimal")
 { V, W, Z }
 { T, V, Z }
 { U, W, Z }
 { T, U, Z }
 { S, Z }

ここで共変量による調整を行わない回帰分析を行ってみる。

> lm(Y ~ X, data = g_tbl) %>%
+   summary()

Call:
lm(formula = Y ~ X, data = g_tbl)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.0245 -1.2404  0.0074  1.2236  5.8945 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.01985    0.05830   -0.34    0.734    
X            1.42181    0.04632   30.69   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.842 on 998 degrees of freedom
Multiple R-squared:  0.4856, Adjusted R-squared:  0.4851 
F-statistic: 942.1 on 1 and 998 DF,  p-value: < 2.2e-16

推定された、効果は1.42と設定した真値より高い値で推定されてしまっている。

次にバックドアパスを閉じるための頂点集合であるZとSを共変量として、回帰分析を行う。

> lm(Y ~ X + Z + S, data = g_tbl) %>%
+   summary()

Call:
lm(formula = Y ~ X + Z + S, data = g_tbl)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.6457 -0.9319  0.0141  0.9189  3.9391 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.01762    0.04475  -0.394    0.694    
X            0.78253    0.04610  16.975   <2e-16 ***
Z            1.18934    0.04501  26.422   <2e-16 ***
S            0.04100    0.04917   0.834    0.405    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.413 on 996 degrees of freedom
Multiple R-squared:  0.6977, Adjusted R-squared:  0.6968 
F-statistic: 766.2 on 3 and 996 DF,  p-value: < 2.2e-16

XからYの因果効果は0.78と真の値に近づいた。

なお、バックドアを閉じるのに不十分な集合であるW, Vのみを共変量とした時の回帰分析を行うと次にようになる。

> lm(Y ~ X +  W + V, data = g_tbl) %>% 
+   summary()

Call:
lm(formula = Y ~ X + W + V, data = g_tbl)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.5458 -0.7703 -0.0353  0.7145  4.0233 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.004678   0.035884    0.13    0.896    
X           0.992874   0.031208   31.81   <2e-16 ***
W           0.845427   0.026312   32.13   <2e-16 ***
V           0.633488   0.032948   19.23   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.132 on 996 degrees of freedom
Multiple R-squared:  0.806,  Adjusted R-squared:  0.8054 
F-statistic:  1379 on 3 and 996 DF,  p-value: < 2.2e-16

銀行の顧客マーケティングデータに適用

次によくみる銀行の顧客のマーケティングデータに、適用してみたいと思います。

このデータを用いて婚姻が収入($50K以上か以下か)に与える因果効果を推定します。

私の感覚から次のような因果構造を仮定しました。

graph <- dagitty('dag {
   married -> income 
   sex -> income
   jobclass -> income
   sex -> jobclass-> married
   jobclass -> hoursweek
   married -> hoursweek
   sex -> job -> married
   job -> income
   edunum -> job
   edunum -> married
   edunum -> jobclass
   edunum -> income
   age -> income
   age -> jobclass
   age -> married
}')

plot(graphLayout(graph))

f:id:saltcooky:20200413010816p:plain

この因果構造を仮定したときの婚姻から収入へのバックドアパスを満たす変数集合は次にようになった。

> adjustmentSets(graph, "married", "income", type = "all")
 { age, edunum, job, jobclass }
 { age, edunum, hoursweek, job, jobclass }
 { age, edunum, job, jobclass, sex }
 { age, edunum, hoursweek, job, jobclass, sex }

こちらでも共変量による調整を行わないロジスティック回帰分析を行ってみる。

log_simple_model <- glm(income ~ married,
                 data = bank_df,
                 family = binomial(link = log))

summary(log_simple_model)

ブートストラップ法による偏回帰係数の分布の推定を行ってみます。

beta_simple.boot <- Boot(log_simple_model, R = 2000, method = "case")

> confint(beta_simple.boot, level=.95, type="bca")
Bootstrap percent confidence intervals

                2.5 %    97.5 %
(Intercept) -3.135666 -2.897781
married      2.083186  2.329583

次にバックドア基準を満たす共変量であるage, edunum, job, jobclassを加えたモデルを作成します。

# バックドア基準を満たす共変量を追加
log_bd_model <- glm(income ~ married + age + edunum + job + jobclass,
                    data = bank_df,
                    family = binomial(link = log))

summary(log_bd_model)

こちらもブートストラップ法による偏回帰係数の分布の推定を行ってみます。

beta.boot <- Boot(log_bd_model, R = 2000, method = "case")

> confint(beta.boot, level=.95, type="bca")
Bootstrap percent confidence intervals

                                2.5 %       97.5 %
(Intercept)               -8.68912813  -7.75489475
married                    2.41076550   2.71471287
age                        0.02357196   0.03255189
edunum                     0.27030066   0.32156546
# - 略 -

共変量を調整する前よりも、婚姻が収入に与える影響が大きいことがわかりました。

最後に

因果推論の考え方が少しわかってきた気がします。

今後も少しづつですが、理解していきます。

名前はまだない

データ分析とかの備忘録か, 趣味の話か, はたまた

構造因果モデルを考え描いて推定する時のメモ(未完)