統計的因果探索の勉強を始めたので、その基本的な手法であるLiNGAMを理論に沿って整理します。数学的な議論は厳密なものではありません。詳しく知りたい方は以下の参考文献を参照してください。
間違っている箇所や気づいたことがあればご指摘いただけると幸いです。　

参考文献

LiNGAMモデルの元論文
https://www.jmlr.org/papers/volume7/shimizu06a/shimizu06a.pdf
清水昌平「統計的因果探索」, 講談社
www.kspub.co.jp
上論文のFirstAutherである清水昌平先生によって書かれた著書です。LiNGAMモデルを含む統計的因果探索の基礎を極めて丁寧に解説されています。
適応型Lasso回帰の元論文
http://users.stat.umn.edu/~zouxx019/Papers/adalasso.pdf
LiNGAMのライブラリ
https://github.com/cdt15/lingam

LiNGAMの仮定

LiNGAMモデルでは以下の構造方程式について考えます。

$\displaystyle{ x_i = \sum_{j\not=i} b_{ij}x_j + e_{i} \quad (i = 1,\ldots , p) }$

$x$ は観測変数、 $p$ は観測変数 $x$ の個数、 $b_{ij}$ は観測変数 $x_j$ が $x_i$ に及ぼす直接的な影響の大きさを表す係数、 $e_i$ は誤差変数を表しています。

これは行列表現で

$x = Bx + e\tag{0}$

と書けます。

LiNGAMは $x$ の分布 $p(x)$ を入力として $B$ の値を一意に推定できます。

また、LiNGAMは次のような仮定を置いています。

各変数は連続変数である。
$e$ は互いに独立である(＝未観測共通原因が存在しない)かつ非ガウス分布に従う。
観測変数が他の観測変数と誤差変数の線形和で表現される。
観測変数は非巡回である。

データ生成過程

説明に用いるデータ生成過程の例は以下のものとします。

$\displaystyle x_2 = e_2 \\\ x_1 = 2.0x_{2} + e_1 \\\ x_0 = 5.0x_{1} + 4x_{2} + e_0 \\\ x_5 = 1.5x_{0} + e_5 \\\ x_4 = 3.0x_{0} + e_4 \\\ x_3 = 1.0x_{4} + 2.0x_{2} + e_3$

DAGは以下のようになります。

独立成分分析の導入

$B$ を求めるために、まずは行列式を $x$ について解きます。

$\displaystyle x = Bx + e \\\ \rightarrow(I-B)x = e \\\ \rightarrow x = (I-B)^{-1}e$

$(I-B)^{-1}$ を $A$ と置いて

$\displaystyle x = Ae$

さらに $A$ の逆行列を $W$ $(= I - B)$ として

$\displaystyle e = Wx \tag{1}$

となり、 $W$ がわかれば $B$ をもとめることができそうです。

これを解くために独立成分分析(Independent component analysis, ICA)という手法を援用します。 ICAは主に信号源分離に用いられている計算手法です。
観測変数ベクトル $x$ が、非ガウス分布に従う独立成分ベクトル $s$ と係数 $a$ によって生成されるというモデルを以下のように定義します。

$\displaystyle{ x_i = \sum_{j=1}^q a_{ij}s_j \quad (i = 1,\ldots , q) }$

係数行列 $A$ を用いて行列式に書き直すと次のようになります。

$\displaystyle{ x = As }\tag{2}$

ICAはこの $x$ から

$\displaystyle{ x = A_{ica} s _ {ica} }$

の形に復元することができます。式 $(2)$ 中の $s$ が非ガウス分布に従う独立ベクトルであることに注目して、LiNGAMの式 $(1)$ に適応できます。ここで注意しないといけないのは、ICAによって復元された $A _ {ica}$ は、本来の混合行列 $A$ と列の順序と尺度が異なる可能性があるということです。つまり、ICAで得られた $A _ {ica}$ の逆行列 $W_ {ica}$ は、本来の $B$ と単位行列 $I$ から $W=I-B$ として得られた真の $W$ と比べ、列の順序と尺度が異なる可能性があります。したがって、 $W$ と $W _ {ica}$ の関係は、置換行列 $P$ と対角行列 $D$ を用いて次のように表せます。

$\displaystyle{ W_{ica} = PDW }\tag{3}$

ここから $P$ , $D$ を推定して $W$ の値を求めたいです。独立成分分析において $P$ , $D$ は識別不可能です。しかしLiNGAMモデルにおいては $B$ の性質を考慮することによって $P$ , $D$ を求めることができます。

因果的順序と非巡回性の行列表現

LiNGAMモデルでは 4.観測変数は非巡回という仮定を置いたのでした。これについて考えてみましょう。
初めに観測変数の因果的順序について説明します。因果的順序とは観測関数のデータが生成される順番のことです。
今回のデータ生成過程においては、まずはじめに非ガウス分布に従って $x2$ が生成されます。これはほかの観測変数に依存しない独立な変数で、このことは係数行列 $B$ の $x2$ 行の成分がすべて $0$ であることからも確認できます。次に、 $x1$ は生成された $x2$ に係数 $2.0$ をかけたものと、非ガウス分布によって生成された誤差変数の線形和として生成されます。以降、同様にして、 $x0$ , $x5$ , $x4$ , $x3$ と生成されます。
したがって、このデータの因果的順序の一つは $x2$ , $x1$ , $x0$ , $x5$ , $x4$ , $x3$ となります。この因果的順序の順番にデータが生成される過程はDAGを見ても直感的に理解出来るかと思います。以降は観測変数 $x_i$ の因果的順序を $k(i)　 (i=1,...p)$ とあらわすことにします。
非巡回とは、データの生成過程で因果的順序が下流の変数が上流の変数と自分自身の値に影響を及ぼさない状態のことです。
さて、データ生成の行列表現は式 $(0)$ です。

$x = Bx + e\tag{0}$

式 $(0)$ で観測変数 $x$ ,誤差変数 $e$ を因果的順序 $x _ {ord}$ , $e _ {ord}$ に並び替え、それに対応して係数ベクトル $B$ を $B _ {ord}$ に変換します。

$x_{ord} = B_{ord}x_{ord} + e_{ord}\tag{4}$

$B _ {ord}$ とデータ生成過程を確認します。非巡回は自分自身の値に影響しませんから、それぞれの行において観測変数自身にかかる係数成分は0です。したがって $B _ {ord}$ の対角成分は0となります。では行列式 $(0)$ に生成したデータを代入します。

式 $(5)$ を因果的順序に並び変えてみます。

観測変数は因果的順序が自身よりの下流の変数の影響を受けない(=係数成分が0)であるためそれぞれの行で対角成分をなす成分より右の成分は0となることがわかります。したがって、係数行列 $B_{ord}$ は厳密な下三角を成します。
続いて $W (= I - B)$ を考えます。 $B$ の対角成分は0ですから $W$ の対角成分はすべて1となります。具体例で確認します。

この条件を用いて置換行列 $P$ と対角行列 $D$ を推定します。
置換行列は行の順序を変更するのでした。行の順序が正しくない場合、変換後の行列の対角成分のいずれかが0になることが知られています。なので、 $P$ によって変更された順序を、すべての対角成分が0ではないような順序へ戻すような置換行列 $\tilde{P}$ を考えます。
$\tilde{P}$ を式 $(3)$ の右からかけます。 $\tilde{P}P = I$ であることを利用して

$\tilde{P} W_{ica} = \tilde{P}PDW= DW\tag{7}$

となり、行の順序が $W$ と同じになるように復元できます。(復元する順序は $W$ の順序であって因果的順序に並び替えた $W_{ord}$ の順序ではありません。)
また、実際の対角成分の値が0であっても計算によって生じる推定誤差により0とならない場合があるため、置換行列 $\tilde{P}$ は対角成分ができるだけ0から離れた値になるような推定値 $\hat{\tilde{P}}$ を以下のように求めます。

$\displaystyle \hat{\tilde{P}}={\underset{\tilde{P}}{\textrm{argmin}}} \sum_{i=1}^{p} \dfrac{1}{|( \tilde{P} \hat{W_{ica}})_{ii}|) }$

次に $D$ を復元します。これは復元行列 $W$ の対角成分は1なので、 $\tilde{P} W _ {ica}$ のそれぞれの行の対角成分は行の尺度の値とわかります。したがって式 $(7)$ ]の左から対角行列 $D$ の逆行列 $D^{-1}$ をかけて

$D^{-1} \tilde{P} W_{ica} = D^{-1} DW = W\tag{7}$

となります。あとは $W = I - B$ から $B$ が推定できます。以降ではこの推定値を $\hat{B}$ と表記します。
以上より独立成分分析によって得られた式

$\displaystyle{ W_{ica} = PDW }$

のみでは識別不可能だった $P$ , $D$ が、LiNGAMモデルにおける因果的順序の非巡回性の仮定を条件に識別可能となり、結果として $W _ {ica}$ から $W$ を求め、 $B$ の推定値 $\hat{B}$ を推定することができました。
この時点で $B$ は推定値 $\hat{B}$ として復元されていますが、 $\hat{B}$ の成分は本当の値が0であったとしても推定誤差のために0ではないため、このままだと巡回な関係となります。なので、前述した因果的順序を推定し、並び替えた $B _ {ord}$ の上三角部分の成分を0とすることで厳密な下三角とし、本来の非巡回の関係を推定します。
具体的な手順は次の通りです。 $B _ {ord}$ の上三角成分の値がすべて0ということは、 $B$ の $p^{2}$ 個の成分のうち、すくなくとも対角成分の個数 $p$ 個と、残りの成分の半分にあたる $\dfrac{p^{2}-p}{2}$ 個の合計 $\dfrac{p(p+1)}{2}$ 個の成分の値が0である必要があります。なので $\hat{B}$ の成分のうち絶対値が小さい順に $\dfrac{p(p+1)}{2}$ 個の成分を0で置き換えます。そして $\hat{B}$ が厳密な下三角になるような並び替えが可能かを調べます。もし可能なら並び替え後の行列を $B _ {ord}$ とし、不可能であれば成分のうち次に絶対値が小さい値一つを0としてもう一度並び替えを調べます。これを厳密な下三角が出来るまで繰り返し、因果的順序および $B_ {ord}$ を推定します。

適応型Lasso回帰による冗長性の排除

上の手順で推定した $B_ {ord}$ の下三角にはまだ、本来の値は0であるが推定誤差によって0になってない値があるかもしれません。観測変数の線形回帰によって正しい値を求めたいところですが、最小二乗法といった単純な線形回帰分析を行っても残されてしまいます。
適応型Lasso回帰の式を成り立ちを最小二乗法による線形回帰から軽く見ていきます。なお、適応型Lasso回帰について詳しく知りたい方は参考文献の元論文を参考にしてください。
まずは最小二乗法の式からです。

$\displaystyle \hat{\beta}={\underset{b}{\textrm{argmin}}\|{x_ {i}-{\sum_{k(j) < k(i)}}{b_{ij}x_{j}}}\|^{2} }$

この式の時点で行っていることは係数の推定のみで、推定誤差による冗長な因果関係を排除できていません。 Lasso回帰では、右辺の目的関数にL1正則化項を足し合わせたものになります。

$\displaystyle \hat{\beta}={\underset{b}{\textrm{argmin}}\|{x_ {i}-{\sum_{k(j) < k(i)}}{b_{ij}x_{j}}}\|^{2} + \lambda {\sum_{k(j) < k(i)}}|b_{ij}| }$

L1正則化項が最小となるのはすべての係数が0となるときですから、最小化をしていくと係数が0に近づいていくことが直感的にわかります。実際いくつかの係数は0となるので、本来の値が0である係数を0とすることができます。イメージとしては本当に重要な係数だけが選択され0ではない値として生き残る感じです。つまり、Lasso回帰は係数の推定と同時に、変数の選択(=因果関係の選択)を行っていることになります。
Lasso回帰では係数の選択に一貫性が無い場合があるためさらに発展させた適応型Lasso回帰を導入します。適応型Lasso回帰では、L1正則化項に含まれていた係数ベクトル $b _ {ij}$ それぞれに重み $w_{ij}$ を掛け合わせます。

$\displaystyle \hat{\beta}={\underset{b}{\textrm{argmin}}\|{x_ {i}-{\sum_{k(j) < k(i)}}{b_{ij}x_{j}}}\|^{2} + \lambda {\sum_{k(j) < k(i)}}{\hat{w_{ij}}}|b_{ij}| }$

一般的に、この係数 $\hat{w _ {ij}}$ は推定値 $\hat{b _ {ij}}$ とパラメータ $\gamma$ を用いて次のように定義します。

$\displaystyle \hat{w _ {ij}} = \dfrac{1}{|\hat{b _ {ij}}|^{\gamma}}$

書き直して最終的な適応型Lasso回帰の式を得ます。

$\displaystyle \hat{\beta}={\underset{b}{\textrm{argmin}}\|{x_ {i}-{\sum_{k(j) < k(i)}}{b _ {ij}x_{j}}}\|^{2} + \lambda {\sum_{k(j) < k(i)}} \dfrac{|b_{ij}|}{|\hat{b_{ij}}|^{\gamma}}}$

やってみる！

以上を踏まえて実際に推定をしてみましょう。サンプルサイズは10000としました。

import numpy as np
import pandas as pd
import graphviz
import lingam
from lingam.utils import print_causal_directions, print_dagc, make_dot

x2 = np.random.uniform(size=10000)
x1 = 2.0*x2 + np.random.uniform(10000)
x0 = 4.0*x2 + 5.0*x1 +np.random.uniform(size=10000)
x5 = 1.5*x0 + np.random.uniform(size=10000)
x4 = 3.0*x0  + np.random.uniform(size=10000)
x3 = 2.0*x2 + 1.0*x4 + np.random.uniform(size=10000)

X = pd.DataFrame(np.array([x0,x1,x2,x3,x4,x5]).T, columns=['x0','x1','x2','x3','x4','x5'])
X.head()

model = lingam.ICALiNGAM()
model.fit(X)
model.causal_order_

B_ord =  model.adjacency_matrix_
B_ord

make_dot(model.adjacency_matrix_)

結果のDAGを見ると、想定していたDAGが正しい因果的順序で高い精度で復元できていることがわかります。
異なっている箇所は $x_2$ から $x_5$ への有向辺ですが、推定係数は0.02と小さい値でおさまっていると思います。それから $x_0$ に対する $x_2$ の有向辺の表記が異なりますが、想定しているDAGでは $x_2$ の単位量の変化が $x_0$ に及ぼす影響は $x2$ からの直接の係数と $x_1$ を経由して及ぼす間接的な係数の合計(=4.00+2.005.00=14.00)で、これは復元したDAGでの $x_2$ から $tex_1$ を経由して $x_0$ にかかる係数(2.007.00=14.00)一致しています。