VARMALiNGAM

Model

VARMALiNGAM [3] is an extension of the basic LiNGAM model [1] to time series cases. It combines the basic LiNGAM model with the classic vector autoregressive moving average models (VARMA). It enables analyzing both lagged and contemporaneous (instantaneous) causal relations, whereas the classic VARMA only analyzes lagged causal relations. This VARMALiNGAM model also is an extension of the VARLiNGAM model [2]. It uses VARMA to analyze lagged causal relations instead of VAR. This VARMALiNGAM makes the following assumptions similarly to the basic LiNGAM model [1]: #. Linearity #. Non-Gaussian continuous error variables (except at most one) #. Acyclicity of contemporaneous causal relations #. No hidden common causes between contempraneous error variables

References

[1](1, 2) S. Shimizu, P. O. Hoyer, A. Hyvärinen, and A. J. Kerminen. A linear non-gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7:2003-2030, 2006.
[2]A. Hyvärinen, K. Zhang, S. Shimizu, and P. O. Hoyer. Estimation of a structural vector autoregression model using non-Gaussianity. Journal of Machine Learning Research, 11: 1709-1731, 2010.
[3]Y. Kawahara, S. Shimizu and T. Washio. Analyzing relationships among ARMA processes based on non-Gaussianity of external influences. Neurocomputing, 74(12-13): 2212-2221, 2011. [PDF]

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 make_dot, print_causal_directions, print_dagc

import warnings
warnings.filterwarnings('ignore')

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

np.set_printoptions(precision=3, suppress=True)
np.random.seed(0)
['1.16.2', '0.24.2', '0.11.1', '1.5.2']

Test data

We create test data consisting of 5 variables.

psi0 = np.array([
    [ 0.  ,  0.  , -0.25,  0.  ,  0.  ],
    [-0.38,  0.  ,  0.14,  0.  ,  0.  ],
    [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
    [ 0.44, -0.2 , -0.09,  0.  ,  0.  ],
    [ 0.07, -0.06,  0.  ,  0.07,  0.  ]
])
phi1 = np.array([
    [-0.04, -0.29, -0.26,  0.14,  0.47],
    [-0.42,  0.2 ,  0.1 ,  0.24,  0.25],
    [-0.25,  0.18, -0.06,  0.15,  0.18],
    [ 0.22,  0.39,  0.08,  0.12, -0.37],
    [-0.43,  0.09, -0.23,  0.16,  0.25]
])
theta1 = np.array([
    [ 0.15, -0.02, -0.3 , -0.2 ,  0.21],
    [ 0.32,  0.12, -0.11,  0.03,  0.42],
    [-0.07, -0.5 ,  0.03, -0.27, -0.21],
    [-0.17,  0.35,  0.25,  0.24, -0.25],
    [ 0.09,  0.4 ,  0.41,  0.24, -0.31]
])
causal_order = [2, 0, 1, 3, 4]

# data generated from psi0 and phi1 and theta1, causal_order
X = np.loadtxt('data/sample_data_varma_lingam.csv', delimiter=',')

Causal Discovery

To run causal discovery, we create a VARMALiNGAM object and call the fit() method.

model = lingam.VARMALiNGAM(order=(1, 1), criterion=None)
model.fit(X)
<lingam.varma_lingam.VARMALiNGAM at 0x1acfc3fa6d8>

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

model.causal_order_
[2, 0, 1, 3, 4]

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

# psi0
model.adjacency_matrices_[0][0]
array([[ 0.   ,  0.   , -0.238,  0.   ,  0.   ],
       [-0.392,  0.   ,  0.182,  0.   ,  0.   ],
       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.   ],
       [ 0.523, -0.149,  0.   ,  0.   ,  0.   ],
       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.   ]])
# psi1
model.adjacency_matrices_[0][1]
array([[-0.145, -0.288, -0.418,  0.041,  0.592],
       [-0.324,  0.027,  0.024,  0.231,  0.379],
       [-0.249,  0.191, -0.01 ,  0.136,  0.261],
       [ 0.182,  0.698,  0.21 ,  0.197, -0.815],
       [-0.486,  0.063, -0.263,  0.112,  0.26 ]])
# omega0
model.adjacency_matrices_[1][0]
array([[ 0.247, -0.12 , -0.128, -0.124,  0.037],
       [ 0.378,  0.319, -0.12 , -0.023,  0.573],
       [-0.107, -0.624,  0.012, -0.303, -0.246],
       [-0.22 ,  0.26 ,  0.313,  0.227, -0.057],
       [ 0.255,  0.405,  0.41 ,  0.256, -0.286]])

Using DirectLiNGAM for the residuals_ properties, we can calculate psi0 matrix.

dlingam = lingam.DirectLiNGAM()
dlingam.fit(model.residuals_)
dlingam.adjacency_matrix_
array([[ 0.   ,  0.   , -0.238,  0.   ,  0.   ],
       [-0.392,  0.   ,  0.182,  0.   ,  0.   ],
       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.   ],
       [ 0.523, -0.149,  0.   ,  0.   ,  0.   ],
       [ 0.   ,  0.   ,  0.   ,  0.   ,  0.   ]])

We can draw a causal graph by utility funciton

labels = ['y0(t)', 'y1(t)', 'y2(t)', 'y3(t)', 'y4(t)', 'y0(t-1)', 'y1(t-1)', 'y2(t-1)', 'y3(t-1)', 'y4(t-1)']
make_dot(np.hstack(model.adjacency_matrices_[0]), lower_limit=0.3, ignore_shape=True, labels=labels)
../_images/varma_dag.svg

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\).

p_values = model.get_error_independence_p_values()
print(p_values)
[[0.    0.517 0.793 0.004 0.001]
 [0.517 0.    0.09  0.312 0.071]
 [0.793 0.09  0.    0.058 0.075]
 [0.004 0.312 0.058 0.    0.011]
 [0.001 0.071 0.075 0.011 0.   ]]

Bootstrap

Bootstrapping

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

model = lingam.VARMALiNGAM()
result = model.bootstrap(X, n_sampling=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.4 or more.

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

We can check the result by utility function.

labels = ['y0(t)', 'y1(t)', 'y2(t)', 'y3(t)', 'y4(t)', 'y0(t-1)', 'y1(t-1)', 'y2(t-1)', 'y3(t-1)', 'y4(t-1)', 'e0(t-1)', 'e1(t-1)', 'e2(t-1)', 'e3(t-1)', 'e4(t-1)']
print_causal_directions(cdc, 100, labels=labels)
y0(t) <--- y2(t-1) (b<0) (100.0%)
y0(t) <--- y4(t-1) (b>0) (100.0%)
y1(t) <--- e4(t-1) (b>0) (100.0%)
y2(t) <--- e1(t-1) (b<0) (100.0%)
y3(t) <--- y0(t) (b>0) (100.0%)
y3(t) <--- y1(t-1) (b>0) (100.0%)
y3(t) <--- y4(t-1) (b<0) (100.0%)
y4(t) <--- y0(t-1) (b<0) (100.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.3 or more.

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

We can check the result by utility function.

print_dagc(dagc, 100, labels=labels)
DAG[0]: 40.0%
    y0(t) <--- y2(t-1) (b<0)
    y0(t) <--- y4(t-1) (b>0)
    y1(t) <--- y0(t) (b<0)
    y1(t) <--- y0(t-1) (b<0)
    y1(t) <--- y4(t-1) (b>0)
    y1(t) <--- e0(t-1) (b>0)
    y1(t) <--- e1(t-1) (b>0)
    y1(t) <--- e4(t-1) (b>0)
    y2(t) <--- e1(t-1) (b<0)
    y2(t) <--- e3(t-1) (b<0)
    y3(t) <--- y0(t) (b>0)
    y3(t) <--- y1(t-1) (b>0)
    y3(t) <--- y4(t-1) (b<0)
    y3(t) <--- e2(t-1) (b>0)
    y4(t) <--- y0(t-1) (b<0)
    y4(t) <--- e1(t-1) (b>0)
    y4(t) <--- e2(t-1) (b>0)
DAG[1]: 19.0%
    y0(t) <--- y2(t-1) (b<0)
    y0(t) <--- y4(t-1) (b>0)
    y1(t) <--- y0(t) (b<0)
    y1(t) <--- y0(t-1) (b<0)
    y1(t) <--- y4(t-1) (b>0)
    y1(t) <--- e0(t-1) (b>0)
    y1(t) <--- e4(t-1) (b>0)
    y2(t) <--- e1(t-1) (b<0)
    y2(t) <--- e3(t-1) (b<0)
    y3(t) <--- y0(t) (b>0)
    y3(t) <--- y1(t-1) (b>0)
    y3(t) <--- y4(t-1) (b<0)
    y3(t) <--- e2(t-1) (b>0)
    y4(t) <--- y0(t-1) (b<0)
    y4(t) <--- e1(t-1) (b>0)
    y4(t) <--- e2(t-1) (b>0)
DAG[2]: 7.0%
    y0(t) <--- y2(t) (b<0)
    y0(t) <--- y2(t-1) (b<0)
    y0(t) <--- y4(t-1) (b>0)
    y1(t) <--- y0(t) (b<0)
    y1(t) <--- y0(t-1) (b<0)
    y1(t) <--- y4(t-1) (b>0)
    y1(t) <--- e0(t-1) (b>0)
    y1(t) <--- e1(t-1) (b>0)
    y1(t) <--- e4(t-1) (b>0)
    y2(t) <--- e1(t-1) (b<0)
    y2(t) <--- e3(t-1) (b<0)
    y3(t) <--- y0(t) (b>0)
    y3(t) <--- y1(t-1) (b>0)
    y3(t) <--- y4(t-1) (b<0)
    y3(t) <--- e2(t-1) (b>0)
    y4(t) <--- y0(t-1) (b<0)
    y4(t) <--- e1(t-1) (b>0)
    y4(t) <--- e2(t-1) (b>0)

Probability

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

prob = result.get_probabilities(min_causal_effect=0.1)
print('Probability of psi0:\n', prob[0])
print('Probability of psi1:\n', prob[1])
print('Probability of omega1:\n', prob[2])
Probability of psi0:
 [[0.   0.   1.   0.   0.  ]
 [1.   0.   0.95 0.   0.  ]
 [0.   0.   0.   0.   0.  ]
 [1.   0.96 0.24 0.   0.  ]
 [0.16 0.03 0.1  0.04 0.  ]]
Probability of psi1:
 [[1.   1.   1.   0.   1.  ]
 [1.   0.   0.   1.   1.  ]
 [1.   1.   0.   1.   1.  ]
 [1.   1.   1.   1.   1.  ]
 [1.   0.19 1.   0.96 1.  ]]
Probability of omega1:
 [[1.   0.77 1.   0.96 0.  ]
 [1.   1.   1.   0.   1.  ]
 [1.   1.   0.   1.   1.  ]
 [1.   1.   1.   1.   0.04]
 [1.   1.   1.   1.   1.  ]]

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)
df = pd.DataFrame(causal_effects)

df['from'] = df['from'].apply(lambda x : labels[x])
df['to'] = df['to'].apply(lambda x : labels[x])
df
from to effect probability
0 y4(t-1) y2(t) 0.377029 1.00
1 y2(t) y3(t) -0.238642 1.00
2 y1(t) y3(t) -0.213468 1.00
3 y0(t) y3(t) 0.563522 1.00
4 y3(t-1) y4(t) 0.343541 1.00
5 y0(t-1) y2(t) -0.254723 1.00
6 y4(t-1) y1(t) 0.438051 1.00
7 y3(t-1) y1(t) 0.266735 1.00
8 y1(t-1) y1(t) 0.312631 1.00
9 y0(t-1) y4(t) -0.531720 1.00
10 y1(t-1) y4(t) 0.226082 1.00
11 y2(t) y1(t) 0.231064 1.00
12 y0(t) y1(t) -0.310366 1.00
13 y4(t-1) y0(t) 0.210816 1.00
14 y3(t-1) y0(t) 0.375119 1.00
15 y2(t-1) y0(t) -0.377158 1.00
16 y2(t-1) y4(t) -0.368007 1.00
17 y0(t-1) y1(t) -0.419723 1.00
18 y1(t-1) y2(t) 0.329416 0.99
19 y0(t-1) y0(t) -0.188156 0.99
20 y1(t-1) y3(t) 0.120133 0.98
21 y0(t-1) y3(t) 0.217037 0.98
22 y4(t-1) y3(t) -0.186410 0.97
23 y3(t-1) y2(t) 0.184045 0.97
24 y4(t-1) y4(t) 0.287224 0.92
25 y2(t) y0(t) -0.147135 0.91
26 y3(t) y4(t) 0.056672 0.73
27 y3(t-1) y3(t) -0.139039 0.63
28 y0(t) y4(t) 0.086335 0.46
29 y2(t-1) y1(t) 0.081208 0.41
30 y1(t-1) y0(t) -0.040277 0.26
31 y2(t) y4(t) -0.088182 0.20
32 y2(t-1) y2(t) -0.052064 0.19
33 y1(t) y4(t) -0.056033 0.05
34 y4(t) y3(t) 0.057538 0.04
35 y2(t-1) y3(t) -0.261473 0.02
36 y4(t) y1(t) 0.013746 0.01

We can easily perform sorting operations with pandas.DataFrame.

df.sort_values('effect', ascending=False).head()
from to effect probability
3 y0(t) y3(t) 0.563522 1.0
6 y4(t-1) y1(t) 0.438051 1.0
0 y4(t-1) y2(t) 0.377029 1.0
14 y3(t-1) y0(t) 0.375119 1.0
4 y3(t-1) y4(t) 0.343541 1.0

And with pandas.DataFrame, we can easily filter by keywords. The following code extracts the causal direction towards y2(t).

df[df['to']=='y2(t)'].head()
from to effect probability
0 y4(t-1) y2(t) 0.377029 1.00
5 y0(t-1) y2(t) -0.254723 1.00
18 y1(t-1) y2(t) 0.329416 0.99
23 y3(t-1) y2(t) 0.184045 0.97
32 y2(t-1) y2(t) -0.052064 0.19

Because it holds the raw data of the 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 = 5 # index of y0(t-1). (index:0)+(n_features:5)*(lag:1) = 5
to_index = 2 # index of y2(t). (index:2)+(n_features:5)*(lag:0) = 2
plt.hist(result.total_effects_[:, to_index, from_index])
../_images/varma_hist.png