⬜︎⬜︎⬜︎

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

LiNGAMで因果グラフを探索してみる with R

R データ分析 因果推論

はじめに

Pearl, Raubinときたので、LiNGAMについて知っておこうと思い以下の書籍を読みました。

統計的因果探索 (機械学習プロフェッショナルシリーズ)

統計的因果探索 (機械学習プロフェッショナルシリーズ)

  • 作者:清水 昌平
  • 発売日: 2017/05/25
  • メディア: 単行本（ソフトカバー）

いつも通り簡単にまとめて、Rで実行実行してみます。

統計的因果探索

ある2つの変数x,yの因果関係を考えた時に、ある共変量となるzの存在が考えられる。

次の三つの構造的因果モデルが考えられる。

この候補となる因果モデルの中で、得られているデータを生成したと考えられるモデルを特定したい。

このように、得られているデータを用いて、候補となる構造的因果モデルの中から適当なモデルを推定することを統計的因果探索と呼ぶ。

因果探索のの三つのアプローチ

ノンパラメトリック

線形性や誤差項に分布を仮定しないアプローチ

パラメトリック

線形性や誤差項にガウス分布を仮定するアプローチ 二変数の因果構造において、分布は同じものになってしまう。 そのため、因果構造を推定することができない。

セミパラメトリック

関数形には分布を仮定し、誤差項に分布を仮定しないアプローチ。 特にLiNGAMは関数形には分布を仮定し、誤差項に非ガウス分布を仮定する。 非ガウス分布を仮定することにより因果構造の推定を行うことができる。

LiNGAM

LiNGAMはセミパラメトリックアプローチの因果グラフの識別が可能なモデルである。

p個の観測変数x_1,x_2,\dots, x_pに対するLiNGAMモデルは次のようになる。

\displaystyle{ x_i = \sum_{j \neq i} b_{ij} x_j + e_i (i=1,\dots,p) }

誤差変数e_iが独立で、非ガウス連続分布に従い、未観測共通原因がないことを意味する。

なお、行列表現は次にようになる。

\displaystyle{ x = Bx+e }

この時、因果グラフは有向非巡回路グラフであると仮定する。

LiNGAMの識別可能性

関数系には非線形性、誤差変数には非ガウス分布の仮定を置いている。

ここ、次のような二つのモデルの識別することを考える。

  • モデル１
\displaystyle{ x_1 = e_1 x_2 = 0.8x_1 + e_2 }
  • モデル２
\displaystyle{ x_1 = 0.8x_2 + e_1 x_2 = e_2 }

もし、誤差変数が二つともガウス分布に従っているとする。

この時の観測変数x_1,x_2の分布は次のように同じものとなってしまう。

f:id:saltcooky:20200415015411p:plain

一方で非ガウス分布(ここでは一様分布)を仮定すると、観測変数x_1,x_2の分布は次のようになる。

f:id:saltcooky:20200415015429p:plain

異なる分布をしていることから、モデルの識別を行うことができる。

非正規分布に従っていると仮定することがLiNGAMの重要なポイントの一つ。

LiNGAMモデルの推定

LiNGAMモデルの推定には次のような手順をとる。

  1. 変数x_iの因果的順序k(i)を推定する。 因果的順序に従って変数を並び替えると、係数行列は下三角行列となる。

  2. 係数行列の下三角部分の成分を推定する。 この時独立成分分析の方法を用いて推定する。

因果順序を推定

\displaystyle{ x = Bx+e }

ここで係数ベクトルBは厳密な下三角行列と仮定します。

Bxを左辺に移行すると次のようになる。

\displaystyle{ (I-B)x=e }

行列I-Bの逆行列を両辺にかける。

\displaystyle{ \begin{aligned} (I-B)^{-1}(I-B)x &=(I-B)^{-1} e \\ x&=(I-B)^{-1} e \\ x&=A e \\ \end{aligned} }

上記の式は、独立成分分析モデルと解釈できる。行列Aは独立成分分析モデルの混合行列にあたる。

行列Aの逆行列W=A^{-1}であり、係数行列BとはW=I-Bという関係がある。

ICAで推定することできる復元行列Wの限界は、次のような W_{ICA} である。

\displaystyle{ \begin{aligned} W_{ICA} = PDW =PD(I-B) \end{aligned} }

Pは行の順序を表す置換行列であり、Dは行の尺度を表す対角行列である。

もし、PとDが既知であれば、得たい復元行列Wを推定することができる。

\displaystyle{ \begin{aligned} P^{-1} D^{-1} W_{ICA} &= P^{-1}D^{-1}PDW \\ &=W \end{aligned} }

係数行列Bは、B=I-Wから求めることができる。

Pの算出

推定された、\hat W_{ICA}は行の順序が一意に推定されていない。

そこで、LiNGAMモデルの非巡回性の仮定を利用して、推定された復元行列\hat W_{ICA}の対角成分に０が来ないように並び替える。

具体的には対角成分の値の絶対値ができるだけ大きくな流ように、\hat W_{ICA}の行を並び替える。

式でかくと、次のような置換行列\tilde P' を探すことにあたる。

\displaystyle{ \tilde P' = min_{\tilde P} \frac{1}{|(\tilde P, \tilde W_{ICA})_{ii}|} }

ここで、分母は行列\tilde P, \tilde W_{ICA}の第(i, i)成分を表す。

Dの算出

行列 Wの対角成分が全て 1 であることを用いると、diag(D)=diag(DW)となる。

これを用いて、次のようにDを推定することができる。

\displaystyle{ D=diag(D)=diag(DW)=diag(\tilde P W_{ICA} ) }

三角行列の成分の推定

因果順序の推定ができたら、その順序を利用して回帰分析を行い係数行列Bを推定します。

次にような構造方程式に従い係数b_{ij}を推定する。

\displaystyle{ x_i=\sum_{k(j) < k(i)} b_{ij} x_{j} + e_{j} }

因果的順序が k ( j ) < k ( i ) となるような変数候補を説明変数にとって偏回帰係数を推定する。

係数b_{ij}の実際の推定値が０だったとしても、通常の推定方法では0になることは少ない。

そのため、adaptive Lassoのようなスパース回帰分析を用い係数が０になるように推定する。

adaptive Lassoでは、次のような目的関数を用いる。

\displaystyle{ \left|\left| x_i - \sum_{k(j) < k(i)} b_{ij} x_{j} \right|\right|^2 + \lambda \sum_{k(j) < k(i)} \frac{|b_{ij}|}{| \hat b_{ij} |^ \gamma } }

ここで、 \lambdaと\gammaは調整パラメータであり、\hat b _ { ij } は b _ { ij } の一致推定量の値である。

これにより、因果構造における不要な有効辺を削除できる。これを枝刈り(pruning)とも呼ぶ。

その他のトピック

通常のLiNGAMモデルの他に次のような場合のモデルが提案されている

  • 未観測共通変数があるLiNGAM
  • 巡回モデルを仮定する場合
  • 時系列モデル：VAR LiNGAM
  • 非線形モデル：加法誤差モデル
  • Baysian LiNGAM

RでLiNGAM

RでLiNGAMモデルを利用するにはpcalgパッケージを用います。

次のような因果構造を推定する関数を作成しました。

# LiNGAMモデルの推定
est_lingam <- function(X){

  model <- 
    X %>%
    lingam(verbose = F)
  
  arrow_amt <- 
    model %>% 
    as("amat") %>%
    as.matrix()
  
  colnames(arrow_amt) <- names(X)
  rownames(arrow_amt) <- names(X)
  
  return(arrow_amt)
}

三変数の場合

最初はLiNGAMの話でよくある因果構造を推定してみます。

# パターン1
N <- 100
x1 <- 0.7*runif(N,-1,1)
x2 <- 1.6*x1 + runif(N,-1,1)
x3 <- 1.2*x1 - 0.5*x2 + runif(N,-1,1)

graph_df <- data.frame(x1,x2,x3)
GGally::ggpairs(graph_df)

f:id:saltcooky:20200516180337p:plain

因果構造の推定をしてみます。

arrow_lingam <- est_lingam(graph_df)

> arrow_lingam
Adjacency Matrix 'amat' (3 x 3) of type ‘pag’:
   x1 x2    x3   
x1 .   TRUE  TRUE
x2 .  .      TRUE
x3 .  .     .

因果構造を推定できていますね。

次に、誤差項が一様分布以外の非ガウス分布に従う３変数の因果構造を推定します。

# パターン２
x1 <- 0.7*runif(N,-1,1)
x2 <- 1.6*x1 + rpois(N,3)
x3 <- 1.2*x1 + 0.5*x2 + rgamma(N,2,2)

graph_df <- data.frame(x1,x2,x3)
GGally::ggpairs(graph_df)

arrow_lingam <- est_lingam(graph_df)
arrow_lingam

f:id:saltcooky:20200516180406p:plain

LiNGAMモデルで因果構造推定を推定します。

arrow_lingam <- est_lingam(graph_df)

> arrow_lingam
Adjacency Matrix 'amat' (3 x 3) of type ‘pag’:
   x1 x2    x3   
x1 .   TRUE  TRUE
x2 .  .      TRUE
x3 .  .     .

因果構造を推定できていますね。

少し真の構造を変えてみます。

# パターン３
x1 <- 1.7*runif(N,-1,1)
x2 <- 0.6*x1 + rpois(N,3)
x3 <- 1.2*x1 + rgamma(N,2,2)

graph_df <- data.frame(x1,x2,x3)
GGally::ggpairs(graph_df)

f:id:saltcooky:20200516180421p:plain

LiNGAMモデルで因果構造推定を推定します。

arrow_lingam <- est_lingam(graph_df)

> arrow_lingam
Adjacency Matrix 'amat' (3 x 3) of type ‘pag’:
   x1 x2    x3   
x1 .   TRUE  TRUE
x2 .  .     .    
x3 .   TRUE .

この場合は、因果構造推定を検出できていません。

サンプルサイズが小さい or ガンマ分布に従うデータが正規分布に近いことからLiNGAMの仮定を満たせていないためだと考えられる。

サンプルサイズと因果構造の検出率

次に、サンプル数を変化させた時に因果方向の検出率がどう変わるのかをみてみました。

パターン３の因果構造を真の構造とし、サンプルサイズを50から1000まで10刻みで変化させ、各サンプルサイズでリサンプリング回数を100としてブートストラップ法により検出率を計算した。

その結果を示す。

f:id:saltcooky:20200419164955p:plain

変数1 -> 2への因果構造推定は、偏回帰係数の絶対値が小さいためサンプルサイズが1000程度まで増えないと正確に検出することは難しいようだ。

多変量の場合

多変量の場合として次の因果構造を推定する。

N<-10000
x1 <- 0.7*runif(N,-1,1)
x3 <- 0.9*rchisq(N,3,0)
x2 <- 1.7*x1 + 1.2*x3 - 1.4*rgamma(N,1,2)
x4 <- 1.6*x1 + 1.3*x3 + 1.4*rnbinom(N,2,0.5)
x5 <- 1.4*x1 + 1.2*x2 - 1.4*x3 + 1.2*x7 + 0.4*rgamma(N,2,1)
x6 <- 0.8*runif(N,0,1)
x7 <- 1.5*x4 - 1.2*x6 + 0.6*rpois(N,3)
x8 <- 1.1*x5 + 1.2*x7 + 0.5*rchisq(N,3,0)

graph_df <- data.frame(x1,x2,x3,x4,x5,x6,x7,x8)

推定してみる

arrow_lingam <- est_lingam(graph_df)

結果は次のようになった。

> arrow_lingam
Adjacency Matrix 'amat' (8 x 8) of type ‘pag’:
   x1 x2    x3 x4    x5    x6 x7    x8   
x1 .   TRUE .   TRUE  TRUE .  .     .    
x2 .  .     .  .      TRUE .  .     .    
x3 .   TRUE .   TRUE  TRUE .  .     .    
x4 .  .     .  .     .     .   TRUE .    
x5 .  .     .  .     .     .  .      TRUE
x6 .  .     .  .     .     .   TRUE .    
x7 .  .     .  .     .     .  .      TRUE
x8 .  .     .  .     .     .  .     .

わかりにくいが、いくつかの辺の検出ができていない。

おわりに

非ガウス性の仮定と非巡回性を仮定することで独立成分分析を適用することで、うまく因果構造を推定することができるのは素晴らしい考え方です。

適用先は色々ありそうで、うまく利用できれば効果を発揮しそうですね。

参考