RでBayesLiNGAM

RでBayesLiNGAM
===
> 07May'20: Written


## はじめに
現在、因果探索アプローチに用いられるLiNGAMはRおよびPythonで実装されている。
> いろいろなLiNGAMモデルのR and/or Python
> https://sites.google.com/site/sshimizu06/lingam

今回はその中でもBayesLiNGAMを用いたモデルをRで使用してみる。
原理は以下の文章に記載されている。

>https://pdfs.semanticscholar.org/7122/03d0696feb5f03b1d046e6db59f7e4087cff.pdf
>https://www.jstage.jst.go.jp/article/pjsai/JSAI2014/0/JSAI2014_2G13/_pdf


## 導入環境
OS: Mojave (version; 10.14.6)
R version: 3.6.3 (2020-02-29)
R Studio version: 1.2.5033

## 手順
1. パッケージのダウンロード・解凍
2. ソースコードの書き換え
3. 実行に必要な他のパッケージのインストール
4. 動作確認

## パッケージのダウンロード・解凍
以下のサイトの「Code Package」から「Version 1.1」をダウンロード(07May'20)。
https://www.cs.helsinki.fi/group/neuroinf/lingam/bayeslingam/

「bayeslingam.tar.gz」がダウンロードされるので、これを解凍。
→ 「bayeslingam」

## ソースコードの書き換え
このままRで起動しようとしてもErrorになるので、一部ソースコードを書き換える。
書き換えるものは2つ。
・loud.R
・loadbayeslingam.R

### loud.R
bayeslingam/main/loud.Rの一部を書き換え。
元は以下。

```
cauzality_path<<-"./../.."
```
修正後は以下。

```
cauzality_path<<-"/Users/（作業ディレクトリ）/bayeslingam"
```
作業ディレクトリには任意で。

### loadbayeslingam.R
bayeslingam/bayeslingam/loadbayeslingam.Rの一部を書き換え。
元は以下。

```
common_dir<-sprintf("%s/trunk/common",cauzality_path)

bayeslingam_dir<-sprintf("%s/trunk/bayeslingam",cauzality_path)    

rbayeslingam_dir<-sprintf("%s/trunk/bayeslingam/rbayeslingam",cauzality_path)

rlingam_dir<-sprintf("%s/trunk/rlingam",cauzality_path)

rdeal_dir<-sprintf("%s/trunk/rdeal",cauzality_path)

rpc_dir<-sprintf("%s/trunk/rpc",cauzality_path)

etudag_dir<-sprintf("%s/trunk/ETUDAG",cauzality_path)
```
trunk/を消去する。
修正後は以下。

```
common_dir<-sprintf("%s/common",cauzality_path)

bayeslingam_dir<-sprintf("%s/bayeslingam",cauzality_path)    

rbayeslingam_dir<-sprintf("%s/bayeslingam/rbayeslingam",cauzality_path)

rlingam_dir<-sprintf("%s/rlingam",cauzality_path)

rdeal_dir<-sprintf("%s/rdeal",cauzality_path)

rpc_dir<-sprintf("%s/rpc",cauzality_path)

etudag_dir<-sprintf("%s/ETUDAG",cauzality_path)
```

## 実行に必要な他のパッケージのインストール
以下のパッケージをインストールする。
今回は'fastICA'のみでも可能。

|Package |Feature
|--------|:-----------------------------
|deal    | GH-algorithm
|mclust  | MoG Bayeslingam
|combinat| PC-algorithm
|ggm     | PC-algorithm
|fastICA | Standard LiNGAM
|gtools  | Some method's in ETUDAG use this package (not used by BL)

```
install.packages('PACKAGE')
library(PACKAGE)
```

## 動作確認
RStudioでbayeslingamを読み込む。

```
setwd(~/作業ディレクトリ)
library('fastICA')
source('bayeslingam/main/loud.R')
loud()
```

試しに実行してみる。
今回は以下の因果グラフを推定してみる。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/497950/80dafc36-ac2e-6fdc-7ba6-992684979f8b.png)


まずはデータ生成。


```
x1 <- runif(5000, 0, 1)
x2 <- x1 + runif(5000, 0, 1)
x3 <- x1 + x2 + runif(5000, 0, 1)
x4 <- x1 + x3 + runif(5000, 0, 1)
```

データフレームに格納し、BayesLiNGAMを実施。

```
d <- data.frame(x1=x1, x2=x2, x3=x3, x4=x4)
result <- bayeslingam(d)
```

DAGsをみてみる。

```
> result$DAGs
       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
  [1,]    1    2    3    4    0    0    0    0    0     0
  [2,]    1    2    3    4    0    0    0    0    0     1
  [3,]    1    2    3    4    0    0    0    0    1     0
  [4,]    1    2    3    4    0    0    0    0    1     1
  [5,]    1    2    3    4    0    0    0    1    0     0
```
1列目から４列目: 各変数の並び順。
e.g.) 　[2,]はx1,x2,x3,x4を表す。

５列目から１０列目: 有向線の有無を0,1で表す。1列目〜４列目で並べた変数から左の変数への有向線の有無を表す。
e.g.) 　[5,]では[,5]はx1からx2への有向線なし。[,6]はx1からx3への有向線なし。[,7]はx1からx4への有向線なし。[,8]はx2からx3への有向線あり。[,9]はx2からx4への有向線なし。[,10]はx3からx4への有向線なし。

また、result$probはその行が表す因果グラフである推定確率を表す。


最後に、因果グラフを作る。

```
library(igraph)
par(mfrow=c(2,3))
prob <- round(result$prob,digits=4) * 100
n.node <- max(result$components$node)
case <- nrow(result$DAGs)
node <- result$DAGs[,1:n.node]
res <- as.matrix(result$DAGs[,-c(1:n.node)],nrow=case)
name <-paste("X",1:n.node,sep="")

for(i in order(prob,decreasing=T)[1:6]){
    amat <- matrix(0,n.node,n.node)
    index <- node[i,]
    amat[lower.tri(amat,diag=FALSE)] <- res[i,]
    amat <- t(amat[index,index])
    g <- graph.adjacency(amat,mode="directed")
    E(g)$label <- NA
    pos <- layout.circle(g)
    rot <- matrix(c(0,1,-1,0),nrow=2,byrow=T)
    la <- pos %*% rot
    if(n.node == 2)la <- matrix(c(-1,0.5,1,0.5),,nrow=2,byrow=T)
    plot(g,layout=la,vertex.size=70,vertex.label=name)
    mtext(paste(prob[i],"%"),line=0)
}
```

この結果が以下。


![](https://i.imgur.com/9gN9KDo.png)

98%で因果グラフを特定することができた。

## おわりに
基本的には、@m__k さんの「LiNGAMのベイズ的アプローチであるBayesLiNGAMを動かして見た」(https://qiita.com/m__k/items/288910248efafa106863)に沿って実行したらうまくいった。
あと、因子が5つ以上の場合に「greedybayeslingam()」を使用するらしいがうまくいかなかった。

## 参照
・LiNGAMのベイズ的アプローチであるBayesLiNGAMを動かして見た (@m__k さん)
https://qiita.com/m__k/items/288910248efafa106863
・NagoyaR19_因果分析　LiNGAMを試してみた (吉野睦 さん)
https://drive.google.com/file/d/0B5obiqfxcbRcNWZCaTloQlhoai1FN1JJRWJFUk9wSmRTdVFZ/view
・A→Bなのか、B→Aなのかをデータから見抜くことはできるだろうか？（LiNGAMのシミュレーションをしてみた）
https://blog.albert2005.co.jp/2015/02/24/a%E2%86%92b%E3%81%AA%E3%81%AE%E3%81%8B%E3%80%81b%E2%86%92a%E3%81%AA%E3%81%AE%E3%81%8B%E3%82%92%E3%83%87%E3%83%BC%E3%82%BF%E3%81%8B%E3%82%89%E8%A6%8B%E6%8A%9C%E3%81%8F%E3%81%93%E3%81%A8%E3%81%AF/
・BayesLiNGAM (IMOTO Yusuke さん)
https://sites.google.com/view/y-imoto/%E5%82%99%E5%BF%98%E9%8C%B2/bayeslingam