こんにちは、1月よりR&Dユニットのユニットリーダーになりました今井です。
変数A,Bの間に(非線形も含む)相関関係がある場合、(1) A→Bの因果がある、(2) B→Aの因果がある、(3) A→BとB→Aの双方向の因果がある、(4) AとBの間に因果関係は無いという4パターンが考えられます。相関関係があるが因果関係が無いという現象についてなじみが無い方は、以下のサイトが参考になると思います。
因果関係がないのに相関関係があらわれる４つのケースをまとめてみたよ（質問テンプレート付き） ―Take a Risk：林岳彦の研究メモ
今回は話を簡単にするため、「双方向の因果」及び「因果は無いが相関がある」というケースを除いて考えましょう。つまり、A→Bの因果がある、B→Aの因果がある、またはAとBは独立であるという3パターンのどれかであるとします。
ここで質問です。この単純化された世界において、データだけからA→Bの因果があるのかB→Aの因果があるのか（または関係がないか）を識別できるでしょうか？
時系列データであればGranger因果のように時間の不可逆性を用いて向きを決めるということが考えられます。では時系列データ以外でも識別できるでしょうか？
この識別問題は長年難しいと考えられてきましたが、ある条件化では識別可能であるという研究が近年盛んに行われております。（Shimizu et al., 2006 など） こちらに関する概説は、清水先生の2012年行動計量学会のチュートリアルの以下の資料をご覧になってください。
構造方程式モデルによる因果推論: 因果構造探索に関する最近の発展
上記資料にあるとおり様々な条件化での識別性が研究されていますが、今回は単純にLiNGAM (Linear Non-Gaussian Acyclic Model： 線形非ガウシアン非巡回モデル)で、どれだけ識別がうまくいくのかを見てみましょう。
今回識別に用いる手法はBayesian LiNGAMと呼ばれるもので、正規分布とラプラス分布を併せた分布又は混合正規分布を用いて候補の因果ダイアグラムの中でどれが妥当であるかの事後確率を計算します（以下、事前確率は各候補モデルで等しくなるよう設定しています）。こちらの手法の元論文及びRのコードは以下のサイトからダウンロードできます。
Bayeslingam – Bayesian discovery of linear acyclic causal models
x1→x2に線形的な因果関係があり、ノイズが一様分布であるデータを以下のように100個生成させて、相関係数及び散布図を見てみましょう。
set.seed(123) x1 <- runif(100,0,1) x2 <- x1 + runif(100,0,1) cor(x1,x2)
[1] 0.7060377
線形的な関係があることは散布図から分かりますが、これを見てもx1→x2なのかx2→x1なのかはパッと見判断が付きません。そこで、上記のRの関数を用いてBayesian LiNGAMを使ってみましょう。（デフォルトでは混合正規分布を用いたものになります。オプション指定で正規分布とラプラス分布を併せた分布も使えます）
#事前にbayeslingamの関数を読み込んでおいてください d <- data.frame(x1=x1, x2=x2) result <- bayeslingam(d)
ここから得られる因果ダイアグラムの候補の事後確率を可視化したものが以下になります。
x1→x2である事後確率が95.75%と正しい向きが捉えられています！
なぜ向きを捉えることができるのかを少し考察してみましょう。次の2つ例を考えてみます。
#example 1: x1→x2 x1 <- runif(10000,-sqrt(3),sqrt(3)) x2 <- sqrt(2/3)*x1 + runif(10000,-1,1) #example 2: x2→x1 x2 <- runif(10000,-sqrt(3),sqrt(3)) x1 <- sqrt(2/3)*x2 + runif(10000,-1,1)
#example 1 > mean(x1) [1] -0.001519604 > mean(x2) [1] -0.00810296 > var(x1) [1] 0.9945374 > var(x2) [1] 0.9952946 #example 2 > mean(x1) [1] -0.002352132 > mean(x2) [1] -0.001133402 > var(x1) [1] 1.021738 > var(x2) [1] 1.011572
どちらの例もx1,x2の平均はほぼ0、分散はほぼ1に等しいので、周辺分布の2次程度のモーメントを見ても区別が付きそうにありません。しかしこれを可視化してみると次のようになります。
このとおり同時分布の形状は明らかに異なっていることが分かります。この非対称性をモデルにうまく組み込むことができれば、その識別が付くという仕組みです。
では次に、最初に用いたデータ生成のサンプルサイズを変化させて、どれだけ正しく関係が捉えられるかを、サンプルサイズ20,40,60,80,100においてそれぞれ1000回ずつシミュレーションを行い、正解率の推移を見てみましょう。
ここで横軸がサンプルサイズ、縦軸が正解率となっています。サンプルサイズが100の時には96.5%正解しています。このシミュレーション結果を見ると、サンプルサイズが十分大きければ正しい関係を捉えられそうですね。
同様に、関係が無い場合にも正しく関係が捉えられるかも見てみましょう。
この場合もうまくいっていることが分かります。
今回は単純なケースで因果構造探索の一端を見てみましたが、いかがでしたでしょうか。データだけからA→B / B→Aの識別ができるなんてすごいと思いませんか？
もしこのような話に興味があるようでしたら、2015年3月19日(木)-20日(金)に電気通信大学にて開催予定の数学協働プログラム「確率的グラフィカルモデル」（要事前登録）に参加すると面白い講演が聞けると思います。私も20日(金)の午後に登壇いたします。
最後になりましたが、弊社は先日東京証券取引所マザーズ市場へ新規上場をしまして、更なる事業拡大に向けてデータ分析者を積極採用しています。詳細は求人ページを見ていただければと思います。