LiNGAMまで超特急：因果推論入門のはずだった

はじめに

因果推論についてまとめようと思ったけど良い教科書がたくさん出てきたので不要になった記事。

おすすめ教科書

• つくりながら学ぶ! Pythonによる因果分析
• 統計的因果探索
• 入門 統計的因果推論

処置効果

TVCMの効果を測定したい。CMを見た人と見ていない人で商品の購入金額$Y$にどれくらい差があるか考える。
$Z$は割り付けを表す確率変数。TVのCMの効果を検証する例では、$Z=1$はCMを"見た"人ではなく、統計的な抽出対象としてCMを見た集団に"割り付けられた"人のことを指す。$Z=1$に属する人達のことを処置群と呼ぶ。

もし、割り付けを実験者が操作可能であれば、RCTのような割り付けを実施することでこのような複雑なことを考える必要はない。しかし、CMの効果を測定する場合、CMを見たかどうかは操作できるものではなく、観察によって得られるものなので、割り付けが操作できない。このような場合、割り付けの変数$Z$と他の説明変数（性別や年齢）が相関してしまい、それが目的変数$Y$と相関して影響を及ぼしてしまう。これらの影響を考慮して処置効果を計算しなければならない。

CMを見た／見ていない世界での購入金額を表す確率変数には各世界での添字をつけて$Y_i,i\in\{0,1\}$ で表すことにする。
$E[Y_1|Z=1]$はCMを見た集団（が、CMを見た場合）の購入金額を表す。$E[Y_0|Z=1]$は、CMを見た集団が、もしCMを見ていなかった場合の購入金額（観測不可能！）を表す。表でまとめると次のようになる。

ID	割付（$Z$）	購入金額（$Y$）	CMを見た場合の購入金額（$Y_1$）	CMを見なかった場合の購入金額（$Y_0$）	性別（$X_1$）	年齢（$X_2$）
1	1	$Y_{1}=900$	$Y_{1,1}=900$	$Y_{1,0}=$ 観測不可能	F	53
2	1	$Y_{2}=600$	$Y_{2,1}=600$	$Y_{2,0}=$ 観測不可能	F	32
3	1	$Y_{3}=700$	$Y_{3,1}=700$	$Y_{3,0}=$ 観測不可能	M	55
4	1	$Y_{4}=1000$	$Y_{4,1}=1000$	$Y_{4,0}=$ 観測不可能	F	52
5	1	$Y_{5}=800$	$Y_{5,1}=800$	$Y_{5,0}=$ 観測不可能	F	48
6	1	$Y_{6}=700$	$Y_{6,1}=700$	$Y_{6,0}=$ 観測不可能	F	47
7	0	$Y_{7}=300$	$Y_{7,1}=$ 観測不可能	$Y_{7,0}=300$	M	23
8	0	$Y_{8}=0$	$Y_{8,1}=$ 観測不可能	$Y_{8,0}=0$	M	31
9	0	$Y_{9}=0$	$Y_{9,1}=$ 観測不可能	$Y_{9,0}=0$	M	21
10	0	$Y_{10}=200$	$Y_{10,1}=$ 観測不可能	$Y_{10,0}=200$	F	25
11	0	$Y_{11}=400$	$Y_{11,1}=$ 観測不可能	$Y_{11,0}=400$	M	43
12	0	$Y_{12}=100$	$Y_{12,1}=$ 観測不可能	$Y_{12,0}=100$	M	27

処置効果はいくつか種類があり、具体的には次のようになる。

項目	式
ITE（Individual Treatment Effect; 個別処置効果）	$Y_{i, 1} - Y_{i, 0}$
ATE（Average Treatment Effect; 平均処置処置効果）	$E[Y_1 - Y_0] = \frac{1}{12}\sum_{i=1}^{12}(Y_{i, 1} - Y_{i, 0})$
ATT（Average Treatment effect on the Treated; 処置群における平均処置効果）	$E[Y_1 - Y_0 | Z=1] = \frac{1}{6}\sum_{i=1}^{6}(Y_{i, 1} - Y_{i, 0})$
ATU（Average Treatment effect on the Untreated; 対照群における平均処置効果）	$E[Y_1 - Y_0 | Z=0] = \frac{1}{6}\sum_{i=7}^{12}(Y_{i, 1} - Y_{i, 0})$

全データに対して$Z$をランダムに割り付けた場合、各変数の期待値は$Z$によらない。したがって次の関係が成り立つ。

\begin{aligned}
E[\cdot|Z=1] =& E[\cdot],\\
E[\cdot|Z=0] =& E[\cdot]
\end{aligned}

介入

集団全員に対して、ある確率変数の値を決めてしまうことを介入と呼ぶ。$E[\cdot | \mathrm{do}(Z=1)]$は集団全員が$Z=1$であったと（仮想的に）考えたときの期待値を表す。

\begin{aligned}
E[Y_1|\mathrm{do}(Z=1)] =& \frac{1}{12}\sum_{i=1}^{12}Y_{i, 1}, \\
E[Y_0|\mathrm{do}(Z=0)] =& \frac{1}{12}\sum_{i=1}^{12}Y_{i, 0}
\end{aligned}

doオペレータで条件付けられた確率には調整化公式が存在する。

\begin{aligned}
P(Y=y|\mathrm{do}(Z=z)) =& \sum_x P(Y=y|Z=z, X=x)P(X=x) \\
=& \sum_x \frac{P(Y=y, Z=z, X=x)}{P(Z=z|X=x)}
\end{aligned}

介入を利用して、ATEを計算する事ができる。

\begin{aligned}
\mathrm{ATE} =& E[Y_1] - E[Y_0] \\
=& E[Y_1|\mathrm{do}(Z=1)] - E[Y_0|\mathrm{do}(Z=0)] \\
=& \sum_x E(Y_1|Z=1, X=x)P(X=x) - \sum_x E(Y_0|Z=0, X=x)P(X=x) \\
=& \sum_x \frac{E(Y_1 | Z=1, X=x)P(Z=1, X=x)}{P(Z=1|X=x)} - \sum_x \frac{E(Y_0 | Z=0, X=x)P(Z=0, X=x)}{P(Z=0|X=x)} \\
=& \frac{1}{N}\sum_{i=1}^N \frac{y_i Z_i}{P(Z=1|X=x_i)} - \frac{1}{N}\sum_{i=1}^N \frac{y_i (1-Z_i)}{P(Z=0|X=x_i)}
\end{aligned}

独立成分分析

観測された変数 $x_i$ が未観測変数 $s_j$ の結合であると考える。

\begin{aligned}
x_i = \sum_{j=1}^q a_{ij}s_j = a_{i1}s_1 + a_{i2}s_2 + ... + a_{iq}s_q
\end{aligned}

行列で表示する。

\begin{aligned}
\boldsymbol{x} = \boldsymbol{A}\boldsymbol{s}
\end{aligned}

$\boldsymbol{A}$を混合行列、$\boldsymbol{x}$ を観測変数ベクトル、$\boldsymbol{s}$ を未観測変数ベクトルと呼ぶ。混合行列は $p\times q$ の行列で、各列は互いに直交する。

$\boldsymbol{s}$ は未観測変数であるため、$\boldsymbol{s}$ の成分同士の入れ替えおよび成分のスケールを自由に決定する自由度がある。したがって、混合行列にも自由度があって、一意に決定できない。

\begin{aligned}
\boldsymbol{A}' &= \boldsymbol{A}\boldsymbol{D}\boldsymbol{P}, \\
\boldsymbol{s}' &= \boldsymbol{P}^{\top}\boldsymbol{D}^{-1}\boldsymbol{s}.
\end{aligned}

$\boldsymbol{D}$ は $q\times q$ の対角行列、$\boldsymbol{P}$は $q\times q$ の置換行列である。置換行列は直交行列であるので、$\boldsymbol{P}^{-1} = \boldsymbol{P}^\top$である。

未観測変数ベクトル $\boldsymbol{s}$ が多次元ガウス分布に従う場合、$\boldsymbol{s}$ の各成分の独立は無相関を意味するので、$\boldsymbol{s}$ の分散共分散行列は対角行列になる。したがって、未観測変数ベクトルの基底の直交変換の自由度があるため、この自由度も混合行列に存在する。

\begin{aligned}
\boldsymbol{A}' &= \boldsymbol{A}\boldsymbol{D}\boldsymbol{P}\boldsymbol{Q}, \\
\boldsymbol{s}' &= \boldsymbol{Q}^\top\boldsymbol{P}^{\top}\boldsymbol{D}^{-1}\boldsymbol{s}, \\
\exp\left[ \boldsymbol{s}'^{\top}\boldsymbol{s}' \right] &= \exp\left[ \boldsymbol{s}^{\top}\boldsymbol{s} \right]
\end{aligned}

混合行列$\boldsymbol{A}$を推定する。行列$\boldsymbol{W}$を$\boldsymbol{x}$にかけて出来たベクトル$\boldsymbol{s}$が最大限独立になるように$\boldsymbol{W}$を決定する。

\begin{aligned}
\boldsymbol{s} = \boldsymbol{W}\boldsymbol{x}
\end{aligned}

$\boldsymbol{W}$の逆行列が$\boldsymbol{A}$である。

$\boldsymbol{s}$の独立度合いを測る指標として相互情報量を使用する。

\begin{aligned}
I(\boldsymbol{s}) &=  \sum_{j=1}^q H(s_j) - H(\boldsymbol{s}), \\
H(\boldsymbol{s}) &= E\left[-\log p(\boldsymbol{s}) \right], \\
H(s_j) &= E\left[-\log p(s_j)\right].
\end{aligned}

推定量は次のように書ける。

\begin{aligned}
I(\boldsymbol{s}) &= \sum_{j=1}^q H(s_j) - H(\boldsymbol{s}), \\
&=  \sum_{j=1}^q \frac{1}{N}\sum_{n=1}^N \left\{ -\log p\left(\boldsymbol{w}_j^\top \boldsymbol{x}^{(n)}\right) \right\} - \frac{1}{N}\sum_{n=1}^N \left\{ -\log p\left(\boldsymbol{W} \boldsymbol{x}^{(n)}\right) \right\} 
\end{aligned}

LiNGAM

各特徴量が他の特徴量の結合で書けるとする。つまり対角成分が0の行列$B$を用いて次のようにかける。

\begin{aligned}
& \boldsymbol{x} = \boldsymbol{B}\boldsymbol{x} + \boldsymbol{e} \\
\Leftrightarrow & (\boldsymbol{I} - \boldsymbol{B})\boldsymbol{x} = \boldsymbol{e} \\
\Leftrightarrow & \boldsymbol{x} = (\boldsymbol{I} - \boldsymbol{B})^{-1}\boldsymbol{e} \\ 
\end{aligned}

これは$\boldsymbol{W}=\boldsymbol{A}^{-1}=\boldsymbol{I}-\boldsymbol{B}$とした独立成分分析になる。

独立成分分析で$\boldsymbol{W}$を1つ推定できたとする。ただし$\boldsymbol{W}$には次の自由度が存在するのであった。

\begin{aligned}
\boldsymbol{W}'^{-1} &= \boldsymbol{W}^{-1}\boldsymbol{D}\boldsymbol{P}, \\
\Leftrightarrow \boldsymbol{W} &= \boldsymbol{D} \boldsymbol{P} \boldsymbol{W}', \\
\end{aligned}

したがって自由度を残した解は次のように書ける。

\begin{aligned}
\boldsymbol{B} &= \boldsymbol{I} - \boldsymbol{W} \\
&= \boldsymbol{I} - \boldsymbol{D} \boldsymbol{P} \boldsymbol{W}'
\end{aligned}

$\boldsymbol{B}$を下三角になるようにしたい。$\boldsymbol{W}'$の対角成分の絶対値が大きくなるように転置行列$\boldsymbol{P}$を選ぶ。

\boldsymbol{P} = \mathrm{argmin}_\boldsymbol{P} \sum_{i=1}^p \frac{1}{\left|\left(\boldsymbol{P}\boldsymbol{W}'\right)_{ii}\right|}

次に対角成分が$1$になるように対角行列を決定する。

\boldsymbol{D} = \mathrm{diag}\left( \boldsymbol{P}\boldsymbol{W}' \right)^{-1}

これでとりあえず$\boldsymbol{B}$が推定された。推定された$\boldsymbol{B}$は対角成分は0であるが、下三角行列になっていない場合が多い。これは非対角成分で0になるべき成分が数値計算の誤差で0にならないからである。そこで非対角成分のうち絶対値が小さい順に$p(p-1)/2$個の成分を0にする。

次に、元の方程式に左から置換行列$\boldsymbol{\tilde{P}}$をかけて行の順番を変えることができる。

\boldsymbol{\tilde{P}}\boldsymbol{x} = \boldsymbol{\tilde{P}}\boldsymbol{B}\boldsymbol{x} + \boldsymbol{\tilde{P}}\boldsymbol{e}

したがって、推定した$\boldsymbol{B}$の行の順番を下三角行列になるように並び替える。下三角行列にならない場合は、絶対値が小さい成分を0にし、並び替えを試す。これを下三角行列になるまで繰り返す。以上で非巡回型の因果関係が求められる。