【実践データ分析】LiNGAMを用いて因果探索　〜個人の属性が年収に与える影響について〜

概要
-
・LiNGAMとは因果探索するための手法の１つ。
・因果探索とは因果ダイアグラムを発見するためのもの。
・因果ダイアグラムは下に掲載しているようなグラフで、今回の実践分析の最終目標でもある。これをpythonで出力する。
![dag2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/33e506dc-3dd6-8320-2296-bdfb1af63326.png)
・LiNGAMのすごいところは、因果ダイアグラムの作成という人が今まで手作業でやってきたことを、客観的に数字で一意に特定できるようになったこと。
・他、因果ダイアグラム特定の手法として、ベイジアンネットワークや相関分析、CAM_UVなどがある。詳しくは、下URL参照。

https://lingam.readthedocs.io/en/latest/tutorial/pk_bottom_up_parce.html

コード：https://github.com/harukicode0/Causal_Discovery

LiNGAMについて
-
Shohei Shimizuが2006年に提唱した因果探索（Causal Discovery）のための手法。

https://www.jmlr.org/papers/volume7/shimizu06a/shimizu06a.pdf

因果探索とは因果ダイアグラムを発見すること。一方で因果推論とは因果効果の発見や効果検証を意味するので、２つは分離して考えること。

LiNGAMは4つ仮定を置くことで因果ダイアグラムを一意に特定することが可能。

LiNGAMが置いている仮定
-
以下の４つである。

①誤差項が非ガウス分布であること（Non-Gaussian continuous error variables (except at most one)）
②DAG（Directed Acyclic Graph：非巡回有向グラフ）である（Acyclicity）
③未観測共通潜在変数が存在しない（No hidden common causes）
④因果関係は線形和で表現可能（Linearity）

それぞれの仮定について説明する。

## ①誤差項が非ガウス分布であること
勉強時間と成績の関係について考える。今回想定するのは下図。
![スクリーンショット 2023-03-07 11.53.08.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/054a0a8e-4409-c8c6-e1f6-850cc16368dd.png)
この仮定を置くと、因果ダイアグラムの矢印の方向が識別可能になる。

### 誤差項が非ガウス（正規分布ではない）であるとなぜ因果の方向を識別可能なのか

x１とx２の関係について考える。

#### ①誤差が正規分布である場合
２つの因果ダイアグラムを考える。
x１が誤差で与えられ、x２がx１×0.8と誤差の和で表現できるもの、つまりx１→x２。
x２が誤差で与えられ、x１がx２×0.8と誤差の和で表現できるもの、つまりx２→x１。
それぞれの散布図を確認して、因果ダイアグラムが識別可能かどうか確認する。

```python:x1→x2
x1 = np.random.normal(0,1,1000)
x2 = x1 * 0.8 + np.random.normal(0,1,1000)
plt.xlabel('x1', fontsize = 16)
plt.ylabel('x2', fontsize = 16)
plt.scatter(x1,x2)
```

![スクリーンショット 2023-03-07 12.06.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/465df138-3e4d-7281-b8f8-56f1e51a09ca.png)

```python:x2→x1
x2 = np.random.normal(0,1,1000)
x1 = x2 * 0.8 + np.random.normal(0,1,1000)
plt.xlabel('x1', fontsize = 16)
plt.ylabel('x2', fontsize = 16)
plt.scatter(x1,x2)
```
![スクリーンショット 2023-03-07 12.08.07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/b146c57d-552a-4183-a058-d8975f6c500c.png)

#### ②誤差が一様分布である場合
x１とx２の散布図を確認すると、因果ダイアグラムの矢印が違うと散布図の形が違うことが確認できる。つまり、散布図から因果ダイアグラムを識別可能である。

x1→x2の因果ダイアグラムを想定した時の散布図について確認。
```python:x1→x2
x1 = np.random.uniform(0,1,1000)
x2 = x1 * 0.8 + np.random.uniform(0,1,1000)
plt.xlabel('x1', fontsize = 16)
plt.ylabel('x2', fontsize = 16)
plt.scatter(x1,x2)
```

![スクリーンショット 2023-03-07 12.11.39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/ff7edd3d-01a8-1445-4373-151dd674ea76.png)

```python:x2→x1の場合
x2 = np.random.uniform(0,1,1000)
x1 = x2 * 2 + np.random.uniform(0,1,1000)
plt.xlabel('x1', fontsize = 16)
plt.ylabel('x2', fontsize = 16)
plt.scatter(x1,x2)
```
![スクリーンショット 2023-03-07 12.12.27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/55876014-1d50-2050-8932-d0cacc585196.png)


## ②DAG（Directed Acyclic Graph：非巡回有向グラフ）である

DAGとは因果の矢印が一方通行であり、変数から流れた効果はその変数に２度と戻らない構造のこと。

下図はDAGではないことの例。勉強時間が増えると成績が上がり、成績が上がるとモチベーションが高まり、さらに勉強時間が増えるというループが発生している構造。このような構造は因果が逆流しているのでDAGとは言えない。

![スクリーンショット 2023-03-07 12.18.31.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/72b50539-d602-8284-5ff9-4be123ef2f3f.png)

## ③未観測共通潜在変数が存在しない
LiNGAMでは因果に関係する変数は全て観測されていることを前提としている。下図のような、点線で囲まれた未観測の影響する変数が存在しないこと。

![スクリーンショット 2023-03-07 12.20.04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/a347bd77-e3ec-d217-6ec4-a503a00f8de2.png)

## ④因果関係は線形和で表現可能
つまり、因果関係は１次式で表現可能であることを想定している。交互作用や２次関数などの非線形な項は存在しないことを仮定している。

![スクリーンショット 2023-03-07 12.22.39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/b1819a37-17eb-3de2-7fa2-a952a6bb6421.png)

LiNGAMの仕組み
-
大まかな流れを説明する。詳しくは他サイトを確認。
①独立主成分分析(ICA)を実施。
②①の結果出てきた行列の対角成分が０になるように調整。
③所定の数まで行列の中の０の数を増やす。この時、０にするのは絶対値の小さいもの。

さらに大雑把に説明すると、ICAで出てきた行列をごちゃごちゃいじって、因果ダイアグラムを作成する。以下①②ではICAのイメージと、行列で因果ダイアグラムを表現する方法について説明し、LiNGAMの頭と尻尾のイメージについて持ってもらう。

### ①ICAのイメージ

ICAとはデータに新たに軸を引き直す作業において、各軸（各独立主成分）同士の相関が０、かつ、独立になるように操作したもの。

２次元平面におけるPCA（主成分分析）とICAの幾何学的イメージが下図。ICAは軸が第１主成分に直交するのではなく、自由に決定できるというイメージ。

![スクリーンショット 2023-03-07 12.30.57.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/6273078d-10c1-cbb5-10d3-d726a81ec3e6.png)


http://meg.aalip.jp/ICA/

### ②行列で因果ダイアグラムを表現

因果ダイアグラムは行列で表現できる。
ここで注目すべき点は対角成分が０であること、０が多いこと。
以下の２つを理解するとLiNGAMのイメージがしやすい。
・０になると変数同士の関係を断ち切れること
・対角成分の下三角形のところ数字が０ではないと因果の下流に変数があること

![いんがすいろん-12.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/827ce99a-266a-2b6c-81dd-1a5dc11b59d1.jpeg)

シミュレーション
-
以下の２つシチュエーションで各メソッドのシミュレーションする。
①未観測共通変数が存在しない場合
②未観測共通変数が存在する場合

シミュレーションするメソッドは４つ。
・DirectLiNGAM
・BottomUpParceLiNGAM（未観測変数が存在しても良いLiNGAM）
・CAM-UV（因果の途中経路に未観測変数が存在する場合）
・ベイジアンネットワーク

以下はデータ数８００で、各シチュエーションで各メソッドでシミュレーションした時の結果を掲載している。

## ①未観測共通変数が存在しない場合

正解の因果ダイアグラム
x７は関係のないノイズを表現している。
結論はDirectLiNGAMが最も正解に近い構造をシミュレーションしている。

![correct_dag.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/ec24ed10-5a5f-c0f2-acfb-10ec182e077e.png)

・DirectLiNGAM
![direct_LiNGAM_dag.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/663adb5e-abb5-166a-41cb-8139085a0182.png)

・BottomUpParceLiNGAM（未観測変数が存在しても良いLiNGAM）
![BottomUpParceLiNGAM_dag.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/7be7ccd8-fafb-1a98-7c34-eef00e59babf.png)

・CAM-UV（因果の途中経路に未観測変数が存在する場合）
![CAM_UV_dag.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/2305f185-1a06-31a5-c625-fe686cf16524.png)

・ベイジアンネットワーク
![スクリーンショット 2023-03-07 13.08.04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/b2c1a665-dd94-00df-c793-4df6b62a8379.png)

## ②未観測共通変数が存在する場合
正解の因果ダイアグラム。
x５を未観測共通変数としている。
DirectLiNGAMが最も因果構造を正しく推定している。

・DirectLiNGAM
![direct_LiNGAM_dag.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/849cb6a1-3913-252c-c43c-e16bf3179b3e.png)

・BottomUpParceLiNGAM（未観測変数が存在しても良いLiNGAM）
![BottomUpParceLiNGAM_dag.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/cc545f9d-401e-1abf-0470-c4112cbf5f55.png)

・CAM-UV（因果の途中経路に未観測変数が存在する場合）
![CAM_UV_dag.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/1621c146-2966-e7e5-bbcc-dd228bfe762d.png)

・ベイジアンネットワーク
![スクリーンショット 2023-03-07 13.19.41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/3ed6c2f4-2a13-453b-e13b-0c910d6f2020.png)

以上よりDirectLiNGAMを利用して因果構造を推定する。

今回利用する実データ
-
American Community Survey 2012を利用する。

https://www.openintro.org/data/index.php?data=acs12

このデータはアメリカの個人の属性データと年収を調査したもの。

分析に使うデータは下表を確認。行った前処理としてはダミー変数化と欠損値のある行の削除、欠損値の多いカラムの削除。

![スクリーンショット 2023-03-07 13.32.20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/0fedb324-2bfa-3969-0e9b-bb69835b09dc.png)

結果と解釈
-
![dag2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/2c4cd881-6686-7fb0-ce53-8ee99eb5fce9.png)


概ね先行研究で示された結果が確認できた。
見方を少しだけ確認する

・年齢が１つ上がるごとに平均＋３７５ドルの影響がある。
・週の労働時間が１時間増えると年収が平均＋９３６ドルの影響がある
・男性の方が女性よりも平均＋１.３万ドル稼ぐ

といったことがわかる。

この結果に１点気になる点がある。それはアジア系であることが年収に対して平均＋１.７万ドルの影響があること。このことについて検証する。


アジア人であることが年収に影響を与えるのか〜交差項について検証
-
先行研究を調べた結果、アジア系であることが年収に正の影響があることに対して、様々な仮説を立てることができる。



仮説１：裕福なアジア系が世界１の国を目指して移住してくるから？
→検証するためには移民の資産や収入等について調査が必要。アジア系２世代の収入の確認も必要。

仮説２：アジア系は学業努力を惜しまない文化・コミュニティを持っており、アジア系は高学歴。Explaining Asian Americans’ academic advantage over whites（Amy Hsin and Yu Xie,2014）
→交互作用を検証すれば、何らかの知見を得られそう。手元データからも実施できる。検証するにはコミュニティや文化が学業努力に与える影響について調査すること。

仮説３：高収入（エンジニア、金融、医者等の専門職）についている傾向がある。
→検証するためには職業についての調査が必要。

仮説４：アジア系は物価の高いところに住みがち？（https://www.washingtonpost.com/news/wonk/wp/2016/12/29/the-asian-american-advantage-that-is-actually-an-illusion/）

結論としてはアジア系であることが直接高収入に繋がるというわけではなさそう。今回、仮説２に関しては手元データから検証できそうなので、実施する。

### 収入に対する学歴と人種の交互作用を検証
目的はアジア系であることの効果を調べること。
結論は手元のデータから確度の高いことは何も言えない。

目的を調査するために、アジア系で、かつ高学歴である場合は他の人種よりも稼ぎやすいかどうか交互作用検証し、アジア系であることの特異性を発見する。

検証した結果が下グラフ。下グラフより交互作用があることがわかる。

![ERII.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/0b896720-7279-4c74-ea34-3ae61fb3ff6a.png)

これより、高学歴なアジア系が年収を急激に上げていることがわかるが、念の為、高学歴におけるアジア系か否かの別で見た年収のヒストグラムを確認する。

![e=2rhist.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/4b989ce0-6cbc-a3b5-d263-44da88bdb377.png)

すると、アジア系である場合のデータ数が１１と少なく、外れ値に影響されやすいことがわかった。結論としては、このデータからでは妥当性のあることは何も言えないことがわかる。

LiNGAMの弱点
-
今回使用したLiNGAMの弱点を下記

・LiNGAMは未観測共通潜在変数について仮定しないので、今回年収の原因を考えるにあたって当たり前に考えなければならない共変量「能力」や「モチベーション」、「職業」などの交絡変数について考慮していないこと。ただし、BottomUpParceLiNGAMやCAM-UVは考慮している。

・線形和の因果関係を扱うので、交互作用について考慮していない。

・ループが発生するような因果関係についても考慮できない。

あくまで、データさえあれば因果ダイアグラムを数値で客観的に特定できる手法というだけで（未知なデータに対しても実行できるので十分すごい）、因果関係を確実に特定できるものではないことに注意したい。

また、LiNGAMで因果ダイアグラムを特定するには、因果関係があると考えられているデータを観測者が仮定した上で収集する必要がある。

参考文献
-
https://lingam.readthedocs.io/en/latest/tutorial/pk_bottom_up_parce.html

https://www.jmlr.org/papers/volume7/shimizu06a/shimizu06a.pdf

https://zenn.dev/megane_otoko/articles/072_casual_resoning_by_lignam

https://lingam.readthedocs.io/en/latest/tutorial/draw_graph.html

https://www.openintro.org/data/index.php?data=acs12

https://www.census.gov/programs-surveys/acs/technical-documentation/code-lists.2020.html#list-tab-HEIUEV8BP30QK3REE9

https://qiita.com/kirikei/items/3a47aecc317d6a580fbf

https://techblog.cccmk.co.jp/entry/2022/03/22/102151

https://vigne-cla.com/24-1/ 

https://www.slideshare.net/sshimizu2006/ss-60015497

http://meg.aalip.jp/ICA/

つくりながら学ぶ! Pythonによる因果分析 ~因果推論・因果探索の実践入門 (Compass Data Science)