はじめに
Pearl, Raubinときたので、LiNGAMについて知っておこうと思い以下の書籍を読みました。
- 作者:清水 昌平
- 発売日: 2017/05/25
- メディア: 単行本(ソフトカバー)
いつも通り簡単にまとめて、Rで実行実行してみます。
統計的因果探索
ある2つの変数の因果関係を考えた時に、ある共変量となるの存在が考えられる。
次の三つの構造的因果モデルが考えられる。
この候補となる因果モデルの中で、得られているデータを生成したと考えられるモデルを特定したい。
このように、得られているデータを用いて、候補となる構造的因果モデルの中から適当なモデルを推定することを統計的因果探索と呼ぶ。
因果探索のの三つのアプローチ
ノンパラメトリック
線形性や誤差項に分布を仮定しないアプローチ
パラメトリック
線形性や誤差項にガウス分布を仮定するアプローチ 二変数の因果構造において、分布は同じものになってしまう。 そのため、因果構造を推定することができない。
セミパラメトリック
関数形には分布を仮定し、誤差項に分布を仮定しないアプローチ。 特にLiNGAMは関数形には分布を仮定し、誤差項に非ガウス分布を仮定する。 非ガウス分布を仮定することにより因果構造の推定を行うことができる。
LiNGAM
LiNGAMはセミパラメトリックアプローチの因果グラフの識別が可能なモデルである。
p個の観測変数に対するLiNGAMモデルは次のようになる。
誤差変数が独立で、非ガウス連続分布に従い、未観測共通原因がないことを意味する。
なお、行列表現は次にようになる。
この時、因果グラフは有向非巡回路グラフであると仮定する。
LiNGAMの識別可能性
関数系には非線形性、誤差変数には非ガウス分布の仮定を置いている。
ここ、次のような二つのモデルの識別することを考える。
- モデル1
- モデル2
もし、誤差変数が二つともガウス分布に従っているとする。
この時の観測変数の分布は次のように同じものとなってしまう。
一方で非ガウス分布(ここでは一様分布)を仮定すると、観測変数の分布は次のようになる。
異なる分布をしていることから、モデルの識別を行うことができる。
非正規分布に従っていると仮定することがLiNGAMの重要なポイントの一つ。
LiNGAMモデルの推定
LiNGAMモデルの推定には次のような手順をとる。
変数の因果的順序k(i)を推定する。 因果的順序に従って変数を並び替えると、係数行列は下三角行列となる。
係数行列の下三角部分の成分を推定する。 この時独立成分分析の方法を用いて推定する。
因果順序を推定
ここで係数ベクトルは厳密な下三角行列と仮定します。
を左辺に移行すると次のようになる。
行列の逆行列を両辺にかける。
上記の式は、独立成分分析モデルと解釈できる。行列は独立成分分析モデルの混合行列にあたる。
行列の逆行列であり、係数行列とはという関係がある。
ICAで推定することできる復元行列の限界は、次のような である。
は行の順序を表す置換行列であり、は行の尺度を表す対角行列である。
もし、とが既知であれば、得たい復元行列を推定することができる。
係数行列は、から求めることができる。
Pの算出
推定された、は行の順序が一意に推定されていない。
そこで、LiNGAMモデルの非巡回性の仮定を利用して、推定された復元行列の対角成分に0が来ないように並び替える。
具体的には対角成分の値の絶対値ができるだけ大きくな流ように、の行を並び替える。
式でかくと、次のような置換行列を探すことにあたる。
ここで、分母は行列の第成分を表す。
Dの算出
行列 の対角成分が全て 1 であることを用いると、となる。
これを用いて、次のようにを推定することができる。
三角行列の成分の推定
因果順序の推定ができたら、その順序を利用して回帰分析を行い係数行列を推定します。
次にような構造方程式に従い係数を推定する。
因果的順序が < となるような変数候補を説明変数にとって偏回帰係数を推定する。
係数の実際の推定値が0だったとしても、通常の推定方法では0になることは少ない。
そのため、adaptive Lassoのようなスパース回帰分析を用い係数が0になるように推定する。
adaptive Lassoでは、次のような目的関数を用いる。
ここで、とは調整パラメータであり、はの一致推定量の値である。
これにより、因果構造における不要な有効辺を削除できる。これを枝刈り(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)
因果構造の推定をしてみます。
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 . . .
因果構造を推定できていますね。
次に、誤差項が一様分布以外の非ガウス分布に従う3変数の因果構造を推定します。
# パターン2 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
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 . . .
因果構造を推定できていますね。
少し真の構造を変えてみます。
# パターン3 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)
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の仮定を満たせていないためだと考えられる。
サンプルサイズと因果構造の検出率
次に、サンプル数を変化させた時に因果方向の検出率がどう変わるのかをみてみました。
パターン3の因果構造を真の構造とし、サンプルサイズを50から1000まで10刻みで変化させ、各サンプルサイズでリサンプリング回数を100としてブートストラップ法により検出率を計算した。
その結果を示す。
変数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 . . . . . . . .
わかりにくいが、いくつかの辺の検出ができていない。
おわりに
非ガウス性の仮定と非巡回性を仮定することで独立成分分析を適用することで、うまく因果構造を推定することができるのは素晴らしい考え方です。
適用先は色々ありそうで、うまく利用できれば効果を発揮しそうですね。