GroupDirectLiNGAM

Import and settings

In this example, we need to import numpy, pandas, and graphviz in addition to lingam.

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

import warnings
warnings.filterwarnings('ignore')

print([np.__version__, pd.__version__, graphviz.__version__, lingam.__version__])

np.set_printoptions(precision=3, suppress=True)

['1.26.4', '2.2.3', '0.20.3', '1.11.0']

Test data

First, we generate a causal structure with 6 variables. Then we create a dataset with 6 variables from x0 to x5. These variables are grouped as follows:

  • Group 1: x0, x1

  • Group 2: x2

  • Group 3: x3, x4

  • Group 4: x5

np.random.seed(0)

n_samples = 1000
x0 = np.random.uniform(size=n_samples)
x1 = 2.0 * x0 + np.random.uniform(size=n_samples)
x2 = 0.5 * x1 + np.random.uniform(-1, 1, size=n_samples)
x3 = 0.3 * x1 + 0.7 * x2 + np.random.uniform(-2, 2, size=n_samples)
x4 = 1.5 * x0 + 0.8 * x3 + np.random.uniform(-2, 2, size=n_samples)
x5 = -0.6 * x3 - 0.5 * x4 + np.random.uniform(-3, 3, size=n_samples)

X = pd.DataFrame(np.array([x0, x1, x2, x3, x4, x5]).T ,columns=['x0', 'x1', 'x2', 'x3', 'x4', 'x5'])
groups = [[0, 1], [2], [3, 4], [5]]

X.head()
x0 x1 x2 x3 x4 x5
0 0.548814 1.690507 1.468291 1.190806 0.946433 -1.987018
1 0.715189 1.440442 0.672389 1.421278 2.475880 -3.304962
2 0.602763 1.681353 0.886988 2.239635 1.245511 -4.554939
3 0.544883 1.798537 0.400310 2.226009 1.996981 -3.218930
4 0.423655 0.891285 0.655729 1.992046 0.441985 -3.023044 

m = np.array([
    [  0,   0,   0,   0,   0, 0],
    [2.0,   0,   0,   0,   0, 0],
    [  0, 0.5,   0,   0,   0, 0],
    [  0, 0.3, 0.7,   0,   0, 0],
    [1.5,   0,   0, 0.8,   0, 0],
    [  0,   0,   0,-0.6,-0.5, 0]])

dot = make_dot(m, labels=['x0', 'x1', 'x2', 'x3', 'x4', 'x5'])

# Save pdf
dot.render('dag')

# Save png
dot.format = 'png'
dot.render('dag')

dot
../_images/group_lingam.svg

Causal Discovery

To run causal discovery, we create a GroupDirectLiNGAM object and call the fit method.

model = lingam.GroupDirectLiNGAM()
model.fit(X, groups)

<lingam.group_direct_lingam.GroupDirectLiNGAM at 0x1d17f5af890>

Using the causal_order_ properties, we can see the causal order of the groups as a result of the causal discovery.

model.causal_order_

[0, 1, 2, 3]

The causal order of the variables is as follows:

[groups[group_idx] for group_idx in model.causal_order_]

[[0, 1], [2], [3, 4], [5]]

Also, using the adjacency_matrix_ properties, we can see the adjacency matrix as a result of the causal discovery.

model.adjacency_matrix_

array([[ 0.   ,  0.   ,  0.   ,  0.   ,  0.   ,  0.   ],
       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.   ,  0.   ],
       [ 0.   ,  0.482,  0.   ,  0.   ,  0.   ,  0.   ],
       [ 0.   ,  0.194,  0.792,  0.   ,  0.   ,  0.   ],
       [ 1.882,  0.   ,  0.572,  0.   ,  0.   ,  0.   ],
       [ 0.   ,  0.   ,  0.   , -0.622, -0.487,  0.   ]])

make_dot(model.adjacency_matrix_)
../_images/group_lingam2.svg

Bootstrapping

We call bootstrap() method instead of fit(). Here, the third argument specifies the number of bootstrap sampling.

model = lingam.GroupDirectLiNGAM()
result = model.bootstrap(X, groups, 100)

Causal Directions

Since BootstrapResult object is returned, we can get the ranking of the causal directions extracted by get_causal_direction_counts() method. In the following sample code, n_directions option is limited to the causal directions of the top 8 rankings, and min_causal_effect option is limited to causal directions with a coefficient of 0.01 or more.

cdc = result.get_causal_direction_counts(n_directions=8, min_causal_effect=0.01, split_by_causal_effect_sign=True)

We can check the result by utility function.

print_causal_directions(cdc, 100)

x2 <--- x1 (b>0) (100.0%)
x3 <--- x2 (b>0) (100.0%)
x4 <--- x2 (b>0) (100.0%)
x5 <--- x3 (b<0) (100.0%)
x5 <--- x4 (b<0) (100.0%)
x4 <--- x0 (b>0) (99.0%)
x3 <--- x1 (b>0) (49.0%)
x3 <--- x0 (b>0) (23.0%)

Directed Acyclic Graphs

Also, using the get_directed_acyclic_graph_counts() method, we can get the ranking of the DAGs extracted. In the following sample code, n_dags option is limited to the dags of the top 3 rankings, and min_causal_effect option is limited to causal directions with a coefficient of 0.01 or more.

dagc = result.get_directed_acyclic_graph_counts(n_dags=3, min_causal_effect=0.01, split_by_causal_effect_sign=True)

We can check the result by utility function.

print_dagc(dagc, 100)

DAG[0]: 35.0%
    x2 <--- x1 (b>0)
    x3 <--- x1 (b>0)
    x3 <--- x2 (b>0)
    x4 <--- x0 (b>0)
    x4 <--- x2 (b>0)
    x5 <--- x3 (b<0)
    x5 <--- x4 (b<0)
DAG[1]: 27.0%
    x2 <--- x1 (b>0)
    x3 <--- x2 (b>0)
    x4 <--- x0 (b>0)
    x4 <--- x2 (b>0)
    x5 <--- x3 (b<0)
    x5 <--- x4 (b<0)
DAG[2]: 19.0%
    x2 <--- x1 (b>0)
    x3 <--- x0 (b>0)
    x3 <--- x2 (b>0)
    x4 <--- x0 (b>0)
    x4 <--- x2 (b>0)
    x5 <--- x3 (b<0)
    x5 <--- x4 (b<0)

Probability

Using the get_probabilities() method, we can get the probability of bootstrapping.

prob = result.get_probabilities(min_causal_effect=0.01)
print(prob)

[[0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.  ]
 [0.08 1.   0.   0.   0.   0.  ]
 [0.23 0.49 1.   0.   0.   0.  ]
 [0.99 0.07 1.   0.   0.   0.  ]
 [0.01 0.03 0.   1.   1.   0.  ]]

Total Causal Effects

Using the get_total_causal_effects() method, we can get the list of total causal effect. The total causal effects we can get are dictionary type variable. We can display the list nicely by assigning it to pandas.DataFrame. Also, we have replaced the variable index with a label below.

causal_effects = result.get_total_causal_effects(min_causal_effect=0.01)

# Assign to pandas.DataFrame for pretty display
df = pd.DataFrame(causal_effects)
labels = [f'x{i}' for i in range(X.shape[1])]
df['from'] = df['from'].apply(lambda x : labels[x])
df['to'] = df['to'].apply(lambda x : labels[x])
df
from to effect probability
0 x1 x2 0.483013 1.00
1 x1 x3 0.504680 1.00
2 x2 x3 0.813637 1.00
3 x1 x4 0.278019 1.00
4 x2 x4 0.555515 1.00
5 x1 x5 -0.456973 1.00
6 x2 x5 -0.786690 1.00
7 x3 x5 -0.629688 1.00
8 x4 x5 -0.480164 1.00
9 x0 x4 1.855977 0.99
10 x0 x5 -0.913577 0.99
11 x0 x3 0.385821 0.29
12 x0 x2 -0.395791 0.08

We can easily perform sorting operations with pandas.DataFrame.

df.sort_values('effect', ascending=False).head()
from to effect probability
9 x0 x4 1.855977 0.99
2 x2 x3 0.813637 1.00
4 x2 x4 0.555515 1.00
1 x1 x3 0.504680 1.00
0 x1 x2 0.483013 1.00 

df.sort_values('probability', ascending=True).head()
from to effect probability
12 x0 x2 -0.395791 0.08
11 x0 x3 0.385821 0.29
9 x0 x4 1.855977 0.99
10 x0 x5 -0.913577 0.99
0 x1 x2 0.483013 1.00

Because it holds the raw data of the total causal effect (the original data for calculating the median), it is possible to draw a histogram of the values of the causal effect, as shown below.

import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
%matplotlib inline

from_index = 0 # index of x0
to_index = 5 # index of x5
plt.hist(result.total_effects_[:, to_index, from_index])

(array([ 1.,  0.,  0.,  0.,  8., 25., 48., 13.,  4.,  1.]),
 array([-2.527, -2.274, -2.021, -1.769, -1.516, -1.263, -1.011, -0.758,
        -0.505, -0.253,  0.   ]),
 <BarContainer object of 10 artists>)
../_images/group_lingam3.svg

Bootstrap Probability of Path

Using the get_paths() method, we can explore all paths from any variable to any variable and calculate the bootstrap probability for each path. The path will be output as an array of variable indices. For example, the array [3, 0, 1] shows the path from variable X3 through variable X0 to variable X1.

from_index = 0 # index of x0
to_index = 5 # index of x5

pd.DataFrame(result.get_paths(from_index, to_index))
path effect probability
0 [0, 4, 5] -0.899122 0.99
1 [0, 3, 5] -0.280465 0.23
2 [0, 2, 3, 5] 0.181984 0.08
3 [0, 2, 4, 5] 0.118028 0.08
4 [0, 5] -1.663298 0.01