<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 MultiGroupDirectLiNGAM
Import and settings
In this example, we need to import numpy, pandas, and graphviz in addition to lingam.
Step1:  Test data
We generate two datasets consisting of 6 variables.
Step2:  We create a list variable that contains two datasets.
Step3:  Causal Discovery
To run causal discovery for multiple datasets, we create a MultiGroupDirectLiNGAM object and call the fit method.
Step4:  Using the causal_order_ properties, we can see the causal ordering as a result of the causal discovery.
Step5:  Also, using the adjacency_matrix_ properties, we can see the adjacency matrix as a result of the causal discovery. As you can see from the following, DAG in each dataset is correctly estimated.
Step6:  To compare, we run DirectLiNGAM with single dataset concatenating two datasets.
Step7:  You can see that the causal structure cannot be estimated correctly for a single dataset.
Step8:  Independence between error variables
To 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$.
Step9:  Bootstrapping
In MultiGroupDirectLiNGAM, bootstrap can be executed in the same way as normal DirectLiNGAM.
Step10:  Causal Directions
The bootstrap method returns a list of multiple BootstrapResult, so we can get the result of bootstrapping from the list. We can get the same number of results as the number of datasets, so we specify an index when we access the results. We can get the ranking of the causal directions extracted by get_causal_direction_counts().
Step11:  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.
Step12:  Probability
Using the get_probabilities() method, we can get the probability of bootstrapping.
Step13:  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.
Step14:  We can easily perform sorting operations with pandas.DataFrame.
Step15:  And with pandas.DataFrame, we can easily filter by keywords. The following code extracts the causal direction towards x1.
Step16:  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.
Step17:  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.
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import pandas as pd
import graphviz
import lingam
from lingam.utils import print_causal_directions, print_dagc, make_dot
print([np.__version__, pd.__version__, graphviz.__version__, lingam.__version__])
np.set_printoptions(precision=3, suppress=True)
np.random.seed(0)
Explanation: MultiGroupDirectLiNGAM
Import and settings
In this example, we need to import numpy, pandas, and graphviz in addition to lingam.
End of explanation
x3 = np.random.uniform(size=1000)
x0 = 3.0*x3 + np.random.uniform(size=1000)
x2 = 6.0*x3 + np.random.uniform(size=1000)
x1 = 3.0*x0 + 2.0*x2 + np.random.uniform(size=1000)
x5 = 4.0*x0 + np.random.uniform(size=1000)
x4 = 8.0*x0 - 1.0*x2 + np.random.uniform(size=1000)
X1 = pd.DataFrame(np.array([x0, x1, x2, x3, x4, x5]).T ,columns=['x0', 'x1', 'x2', 'x3', 'x4', 'x5'])
X1.head()
m = np.array([[0.0, 0.0, 0.0, 3.0, 0.0, 0.0],
              [3.0, 0.0, 2.0, 0.0, 0.0, 0.0],
              [0.0, 0.0, 0.0, 6.0, 0.0, 0.0],
              [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
              [8.0, 0.0,-1.0, 0.0, 0.0, 0.0],
              [4.0, 0.0, 0.0, 0.0, 0.0, 0.0]])
make_dot(m)
x3 = np.random.uniform(size=1000)
x0 = 3.5*x3 + np.random.uniform(size=1000)
x2 = 6.5*x3 + np.random.uniform(size=1000)
x1 = 3.5*x0 + 2.5*x2 + np.random.uniform(size=1000)
x5 = 4.5*x0 + np.random.uniform(size=1000)
x4 = 8.5*x0 - 1.5*x2 + np.random.uniform(size=1000)
X2 = pd.DataFrame(np.array([x0, x1, x2, x3, x4, x5]).T ,columns=['x0', 'x1', 'x2', 'x3', 'x4', 'x5'])
X2.head()
m = np.array([[0.0, 0.0, 0.0, 3.5, 0.0, 0.0],
              [3.5, 0.0, 2.5, 0.0, 0.0, 0.0],
              [0.0, 0.0, 0.0, 6.5, 0.0, 0.0],
              [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
              [8.5, 0.0,-1.5, 0.0, 0.0, 0.0],
              [4.5, 0.0, 0.0, 0.0, 0.0, 0.0]])
make_dot(m)
Explanation: Test data
We generate two datasets consisting of 6 variables.
End of explanation
X_list = [X1, X2]
Explanation: We create a list variable that contains two datasets.
End of explanation
model = lingam.MultiGroupDirectLiNGAM()
model.fit(X_list)
Explanation: Causal Discovery
To run causal discovery for multiple datasets, we create a MultiGroupDirectLiNGAM object and call the fit method.
End of explanation
model.causal_order_
Explanation: Using the causal_order_ properties, we can see the causal ordering as a result of the causal discovery.
End of explanation
print(model.adjacency_matrices_[0])
make_dot(model.adjacency_matrices_[0])
print(model.adjacency_matrices_[1])
make_dot(model.adjacency_matrices_[1])
Explanation: Also, using the adjacency_matrix_ properties, we can see the adjacency matrix as a result of the causal discovery. As you can see from the following, DAG in each dataset is correctly estimated.
End of explanation
X_all = pd.concat([X1, X2])
print(X_all.shape)
model_all = lingam.DirectLiNGAM()
model_all.fit(X_all)
model_all.causal_order_
Explanation: To compare, we run DirectLiNGAM with single dataset concatenating two datasets.
End of explanation
make_dot(model_all.adjacency_matrix_)
Explanation: You can see that the causal structure cannot be estimated correctly for a single dataset.
End of explanation
p_values = model.get_error_independence_p_values(X_list)
print(p_values[0])
print(p_values[1])
Explanation: Independence between error variables
To 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$.
End of explanation
results = model.bootstrap(X_list, n_sampling=100)
Explanation: Bootstrapping
In MultiGroupDirectLiNGAM, bootstrap can be executed in the same way as normal DirectLiNGAM.
End of explanation
cdc = results[0].get_causal_direction_counts(n_directions=8, min_causal_effect=0.01)
print_causal_directions(cdc, 100)
cdc = results[1].get_causal_direction_counts(n_directions=8, min_causal_effect=0.01)
print_causal_directions(cdc, 100)
Explanation: Causal Directions
The bootstrap method returns a list of multiple BootstrapResult, so we can get the result of bootstrapping from the list. We can get the same number of results as the number of datasets, so we specify an index when we access the results. We can get the ranking of the causal directions extracted by get_causal_direction_counts().
End of explanation
dagc = results[0].get_directed_acyclic_graph_counts(n_dags=3, min_causal_effect=0.01)
print_dagc(dagc, 100)
dagc = results[1].get_directed_acyclic_graph_counts(n_dags=3, min_causal_effect=0.01)
print_dagc(dagc, 100)
Explanation: 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.
End of explanation
prob = results[0].get_probabilities(min_causal_effect=0.01)
print(prob)
Explanation: Probability
Using the get_probabilities() method, we can get the probability of bootstrapping.
End of explanation
causal_effects = results[0].get_total_causal_effects(min_causal_effect=0.01)
df = pd.DataFrame(causal_effects)
labels = [f'x{i}' for i in range(X1.shape[1])]
df['from'] = df['from'].apply(lambda x : labels[x])
df['to'] = df['to'].apply(lambda x : labels[x])
df
Explanation: 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.
End of explanation
df.sort_values('effect', ascending=False).head()
Explanation: We can easily perform sorting operations with pandas.DataFrame.
End of explanation
df[df['to']=='x1'].head()
Explanation: And with pandas.DataFrame, we can easily filter by keywords. The following code extracts the causal direction towards x1.
End of explanation
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
%matplotlib inline
from_index = 3
to_index = 0
plt.hist(results[0].total_effects_[:, to_index, from_index])
Explanation: 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.
End of explanation
from_index = 3 # index of x3
to_index = 1 # index of x0
pd.DataFrame(results[0].get_paths(from_index, to_index))
Explanation: 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.
End of explanation
<END_TASK>