<jupyter_start><jupyter_text>VARLiNGAM  Import and settingsIn this example, we need to import `numpy`, `pandas`, and `graphviz` in addition to `lingam`.<jupyter_code>import numpy as np
import pandas as pd
import graphviz
import lingam
from lingam.utils import make_dot, print_causal_directions, print_dagc

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

np.set_printoptions(precision=3, suppress=True)
np.random.seed(0)<jupyter_output>['1.16.2', '0.24.2', '0.11.1', '1.5.2']<jupyter_text>Test dataWe create test data consisting of 5 variables.<jupyter_code>B0 = [
    [0,-0.12,0,0,0],
    [0,0,0,0,0],
    [-0.41,0.01,0,-0.02,0],
    [0.04,-0.22,0,0,0],
    [0.15,0,-0.03,0,0],
]
B1 = [
    [-0.32,0,0.12,0.32,0],
    [0,-0.35,-0.1,-0.46,0.4],
    [0,0,0.37,0,0.46],
    [-0.38,-0.1,-0.24,0,-0.13],
    [0,0,0,0,0],
]
causal_order = [1, 0, 3, 2, 4]

# data generated from B0 and B1
X = pd.read_csv('data/sample_data_var_lingam.csv')<jupyter_output><empty_output><jupyter_text>Causal DiscoveryTo run causal discovery, we create a `VARLiNGAM` object and call the `fit` method.<jupyter_code>model = lingam.VARLiNGAM()
model.fit(X)<jupyter_output><empty_output><jupyter_text>Using the `causal_order_` properties, we can see the causal ordering as a result of the causal discovery.<jupyter_code>model.causal_order_<jupyter_output><empty_output><jupyter_text>Also, using the `adjacency_matrices_` properties, we can see the adjacency matrix as a result of the causal discovery.<jupyter_code># B0
model.adjacency_matrices_[0]
# B1
model.adjacency_matrices_[1]
model.residuals_<jupyter_output><empty_output><jupyter_text>Using `DirectLiNGAM` for the `residuals_` properties, we can calculate B0 matrix.<jupyter_code>dlingam = lingam.DirectLiNGAM()
dlingam.fit(model.residuals_)
dlingam.adjacency_matrix_<jupyter_output><empty_output><jupyter_text>We can draw a causal graph by utility funciton.<jupyter_code>labels = ['x0(t)', 'x1(t)', 'x2(t)', 'x3(t)', 'x4(t)', 'x0(t-1)', 'x1(t-1)', 'x2(t-1)', 'x3(t-1)', 'x4(t-1)']
make_dot(np.hstack(model.adjacency_matrices_), ignore_shape=True, lower_limit=0.05, labels=labels)<jupyter_output><empty_output><jupyter_text>Independence between error variablesTo check if the LiNGAM assumption is broken, we can get p-values of independence between error variables. The value in the i-th row and j-th column of the obtained matrix shows the p-value of the independence of the error variables $e_i$ and $e_j$.<jupyter_code>p_values = model.get_error_independence_p_values()
print(p_values)<jupyter_output>[[0.    0.065 0.068 0.038 0.249]
 [0.065 0.    0.13  0.88  0.57 ]
 [0.068 0.13  0.    0.321 0.231]
 [0.038 0.88  0.321 0.    0.839]
 [0.249 0.57  0.231 0.839 0.   ]]<jupyter_text>Bootstrap  BootstrappingWe call `bootstrap()` method instead of `fit()`. Here, the second argument specifies the number of bootstrap sampling.<jupyter_code>model = lingam.VARLiNGAM()
result = model.bootstrap(X, n_sampling=100)<jupyter_output><empty_output><jupyter_text>Causal DirectionsSince `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.3 or more.<jupyter_code>cdc = result.get_causal_direction_counts(n_directions=8, min_causal_effect=0.3, split_by_causal_effect_sign=True)<jupyter_output><empty_output><jupyter_text>We can check the result by utility function.<jupyter_code>print_causal_directions(cdc, 100, labels=labels)<jupyter_output>x0(t) <--- x0(t-1) (b<0) (100.0%)
x1(t) <--- x1(t-1) (b<0) (100.0%)
x1(t) <--- x3(t-1) (b<0) (100.0%)
x1(t) <--- x4(t-1) (b>0) (100.0%)
x2(t) <--- x2(t-1) (b>0) (100.0%)
x2(t) <--- x4(t-1) (b>0) (100.0%)
x3(t) <--- x0(t-1) (b<0) (100.0%)
x2(t) <--- x0(t) (b<0) (99.0%)<jupyter_text>Directed Acyclic GraphsAlso, 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.2 or more.<jupyter_code>dagc = result.get_directed_acyclic_graph_counts(n_dags=3, min_causal_effect=0.2, split_by_causal_effect_sign=True)<jupyter_output><empty_output><jupyter_text>We can check the result by utility function.<jupyter_code>print_dagc(dagc, 100, labels=labels)<jupyter_output>DAG[0]: 57.0%
	x0(t) <--- x0(t-1) (b<0)
	x0(t) <--- x3(t-1) (b>0)
	x1(t) <--- x1(t-1) (b<0)
	x1(t) <--- x3(t-1) (b<0)
	x1(t) <--- x4(t-1) (b>0)
	x2(t) <--- x0(t) (b<0)
	x2(t) <--- x2(t-1) (b>0)
	x2(t) <--- x4(t-1) (b>0)
	x3(t) <--- x1(t) (b<0)
	x3(t) <--- x0(t-1) (b<0)
	x3(t) <--- x2(t-1) (b<0)
DAG[1]: 42.0%
	x0(t) <--- x0(t-1) (b<0)
	x0(t) <--- x3(t-1) (b>0)
	x1(t) <--- x1(t-1) (b<0)
	x1(t) <--- x3(t-1) (b<0)
	x1(t) <--- x4(t-1) (b>0)
	x2(t) <--- x0(t) (b<0)
	x2(t) <--- x2(t-1) (b>0)
	x2(t) <--- x4(t-1) (b>0)
	x3(t) <--- x0(t-1) (b<0)
	x3(t) <--- x2(t-1) (b<0)
DAG[2]: 1.0%
	x0(t) <--- x0(t-1) (b<0)
	x0(t) <--- x3(t-1) (b>0)
	x1(t) <--- x1(t-1) (b<0)
	x1(t) <--- x3(t-1) (b<0)
	x1(t) <--- x4(t-1) (b>0)
	x2(t) <--- x0(t) (b<0)
	x2(t) <--- x2(t-1) (b>0)
	x2(t) <--- x4(t-1) (b>0)
	x3(t) <--- x1(t) (b<0)
	x3(t) <--- x0(t-1) (b<0)
	x3(t) <--- x2(t-1) (b<0)
	x4(t) <--- x0(t) (b>0)<jupyter_text>ProbabilityUsing the `get_probabilities()` method, we can get the probability of bootstrapping.<jupyter_code>prob = result.get_probabilities(min_causal_effect=0.1)
print('Probability of B0:\n', prob[0])
print('Probability of B1:\n', prob[1])<jupyter_output>Probability of B0:
 [[0.   0.98 0.   0.02 0.  ]
 [0.   0.   0.   0.   0.  ]
 [1.   0.   0.   0.   0.01]
 [0.1  1.   0.   0.   0.  ]
 [0.51 0.   0.02 0.08 0.  ]]
Probability of B1:
 [[1.   0.   0.02 1.   0.  ]
 [0.   1.   1.   1.   1.  ]
 [0.03 0.   1.   0.05 1.  ]
 [1.   0.16 1.   0.   1.  ]
 [0.   0.   0.   0.25 0.  ]]<jupyter_text>Total Causal EffectsUsing the `get_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.<jupyter_code>causal_effects = result.get_total_causal_effects(min_causal_effect=0.01)
df = pd.DataFrame(causal_effects)

df['from'] = df['from'].apply(lambda x : labels[x])
df['to'] = df['to'].apply(lambda x : labels[x])
df<jupyter_output><empty_output><jupyter_text>We can easily perform sorting operations with pandas.DataFrame.<jupyter_code>df.sort_values('effect', ascending=False).head()<jupyter_output><empty_output><jupyter_text>And with pandas.DataFrame, we can easily filter by keywords. The following code extracts the causal direction towards x1(t).<jupyter_code>df[df['to']=='x1(t)'].head()<jupyter_output><empty_output><jupyter_text>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.<jupyter_code>import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
%matplotlib inline

from_index = 7 # index of x2(t-1). (index:2)+(n_features:5)*(lag:1) = 7
to_index = 2 # index of x2(t). (index:2)+(n_features:5)*(lag:0) = 2
plt.hist(result.total_effects_[:, to_index, from_index])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17985898629772515, 'fraction_numerical': 0.0714380736996142, 'mean_word_length': 4.475600873998543, 'pattern_counts': {'":': 0, '<': 134, '<?xml version=': 0, '>': 84, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '414821', 'license': 'MIT', 'n_tokens_mistral': 3509, 'n_tokens_neox': 3279, 'n_words': 736, 'path': 'examples/VARLiNGAM.ipynb', 'repo_name': 'Koji-Kurihara/lingam', 'size': '7517'}