TIGRAMITEによる因果探索　(条件独立性テスト)

　TIGRAMITEは時系列分析Pythonモジュールです。PCMCIフレームワークを基礎にした離散、連続時系列からグラフィカルモデル(条件独立グラフ)を再構築し、結果の高品質の図を生成します。このチュートリアルは、変数の異なる種別に対する条件独立テストを説明しています。

　このテーブルは、X⊥Y | Zのためのテスト、それらの関連した仮定の概要を与えます。

　概要のチュートリアルで議論した密度の図示は、どのテストを使うか選択するのに役立ちます。

　弱い仮定を作るのは、強い仮定でテストするよりも力の検出が典型的に低くなることに注意するのは重要です。例えば、ParCorrは、実際の線形の独立性では、CMIKnnよりもずっと高品質に機能します。さらに、ノンパラメトリックテストは、さらにいっそう計算上高価になります。

# Imports
import numpy as np
import matplotlib
from matplotlib import pyplot as plt
%matplotlib inline     
## use `%matplotlib notebook` for interactive figures
# plt.style.use('ggplot')
import sklearn

import tigramite
from tigramite import data_processing as pp
from tigramite.toymodels import structural_causal_processes as toys

from tigramite import plotting as tp
from tigramite.pcmci import PCMCI

from tigramite.independence_tests.parcorr import ParCorr
from tigramite.independence_tests.robust_parcorr import RobustParCorr
from tigramite.independence_tests.parcorr_wls import ParCorrWLS
from tigramite.independence_tests.gpdc import GPDC
from tigramite.independence_tests.cmiknn import CMIknn
from tigramite.independence_tests.cmisymb import CMIsymb
from tigramite.independence_tests.gsquared import Gsquared
from tigramite.independence_tests.regressionCI import RegressionCI

1.連続値の変数

1.1 ParCorr

　部分相関(ParCorrクラス)は概要チュートリアルで説明されました。

1.2 ロバストParCorr

　もし、まだ線形の従属性があれば、しかしノンガウス分布の、RobustParCorrテストを使うことを助言できます。それは、部分相関テストのために、データを事前に正規分布に変換します。

　ここで、私たちは、新しい関数toys.geneate_structual_causal_processを扱います。それは、パラメータの数とともに構造的因果過程モデルをランダムに生成することができます。以下のノイズ項は、極端なWeibullと一様分布からです。さらに、私たちは、いくつかのポリノミアル変換を適用します。それは、データをもっと極端に生成します。

# np.random.seed(1)
T = 500
links, noises = toys.generate_structural_causal_process(
                N=3, 
                L=4, 
                dependency_funcs=['linear'], 
                dependency_coeffs=[-0.5, 0.5], 
                auto_coeffs=[0.3, 0.5], 
                contemp_fraction=0.,
                max_lag=2, 
                noise_dists=['weibull', 'uniform'],
                noise_means=[0.],
                noise_sigmas=[1., 1.5],
                noise_seed=5,
                seed=9)
for j in links:
    print(j, [(link[0], link[1], link[2].__name__) for link in links[j]])
    print(noises[j].__name__)
data, nonstat = toys.structural_causal_process(links,
             T=T, noises=noises)
data[:,0] = data[:,0] + 0.3*data[:,0]**3
data[:,1] = data[:,1] + 0.3*data[:,1]**5
data[:,2] = np.sign(data[:,2])*data[:,2]**2
var_names = [r'$X^0$', r'$X^1$', r'$X^2$']

dataframe = pp.DataFrame(data, var_names=var_names)
tp.plot_timeseries(dataframe)

0 [((0, -1), 0.5, 'linear'), ((2, -1), -0.5, 'linear'), ((1, -2), 0.5, 'linear')]
weibull
1 [((1, -1), 0.3, 'linear'), ((0, -2), -0.5, 'linear')]
weibull
2 [((2, -1), 0.3, 'linear'), ((0, -2), 0.5, 'linear')]
uniform

# Init
pcmci_robust_parcorr = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=RobustParCorr())

　私たちは、密度図の特徴で結合と周辺密度の歪みを説明することができます。ここで、私たちは、強力に非ガウス分布として表示するために、歪んだデータを標準化します。

# tp.plot_densityplots(dataframe=dataframe, add_densityplot_args={'matrix_lags':matrix_lags}); plt.show()
# correlations = pcmci_robust_parcorr.get_lagged_dependencies(tau_max=20, val_only=True)['val_matrix']
# matrix_lags = np.argmax(np.abs(correlations), axis=2)
matrix_lags = None

matrix = tp.setup_density_matrix(N=dataframe.N, 
        var_names=dataframe.var_names)

# Standardize to better compare skewness with gaussianized data
dataframe.values[0] -= dataframe.values[0].mean(axis=0)
dataframe.values[0] /= dataframe.values[0].std(axis=0)

matrix.add_densityplot(dataframe=dataframe, 
    matrix_lags=matrix_lags, label_color='red', label="standardized data",
    snskdeplot_args = {'cmap':'Reds', 'alpha':1., 'levels':4})

# Transform data to normal marginals
data_normal = pp.trafo2normal(data)
dataframe_normal = pp.DataFrame(data_normal, var_names=var_names)

matrix.add_densityplot(dataframe=dataframe_normal, 
    matrix_lags=matrix_lags, label_color='black', label="gaussianized data",
    snskdeplot_args = {'cmap':'Greys', 'alpha':1., 'levels':4})
matrix.adjustfig()
# plt.show()

pcmci_parcorr = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=ParCorr())
results = pcmci_parcorr.run_pcmci(tau_max=2)

pcmci_robust_parcorr = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=RobustParCorr())
results_robust = pcmci_robust_parcorr.run_pcmci(tau_max=2)

fig, axes = plt.subplots(nrows=1, ncols=3)
axes[0].set_title("True graph")
tp.plot_graph(
    graph=pcmci_robust_parcorr.get_graph_from_dict(links, tau_max=None),
    var_names=var_names,
    fig_ax=(fig, axes[0]),
    show_colorbar=False,
    )


axes[1].set_title("ParCorr")
tp.plot_graph(
    val_matrix=results['val_matrix'],
    graph=results['graph'],
    var_names=var_names,
    fig_ax=(fig, axes[1]),
    show_colorbar=False,
    )

axes[2].set_title("RobustParCorr")
tp.plot_graph(
    val_matrix=results_robust['val_matrix'],
    graph=results_robust['graph'],
    var_names=var_names,
    fig_ax=(fig, axes[2]),
    show_colorbar=False,
    )
plt.show()

　RobustParCorr(右）は、周辺分布が前もって正規分布に変換されるので、標準のParCorr(左)よりもっと信頼性のある結果に分岐します。

　チュートリアルtigmaite_tutorial_heteroskedastic_ParCorrWLSは、私たちは他のバージョンの部分相関をカバーします。それは、ParCorrWLSと呼ばれ、不均一なデータに適用できます。

非線形テスト

　もし、非線形従属性が与えられると、ノンパラメトリックテストを使うことが助言されます。以下のモデルを考えましょう。

random_state = np.random.default_rng(seed=42)
data = random_state.standard_normal((500, 3))
for t in range(1, 500):
    data[t, 0] += 0.4*data[t-1, 1]**2
    data[t, 2] += 0.3*data[t-2, 1]**2
var_names = [r'$X^0$', r'$X^1$', r'$X^2$']

dataframe = pp.DataFrame(data, var_names=var_names)
tp.plot_timeseries(dataframe); plt.show()

pcmci_parcorr = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=ParCorr(),
    verbosity=0)
results = pcmci_parcorr.run_pcmci(tau_max=2, pc_alpha=0.2, alpha_level = 0.01)
tp.plot_graph(
    val_matrix=results['val_matrix'],
    graph=results['graph'],
    var_names=var_names,
    show_colorbar=False,
    )

　ここで、ParCorrは、二つのやり方で失敗します。(1)二つの非線形リンクを検出できない。(2)リンクX_t-1⁰ -> X_t² を間違って検出します。なぜなら、非線形従属性を制約外にできません。

1.3 GPDC

　Trigmaiteは、剰余の距離の相関とガウス過程回帰を基礎にしたテストの非線形加法従属性をカバーします。GPDCにとって距離相関(DC)の非分析、null-分布が利用できます。重要なテストは、パラメータsignificance ='analytic'に設定したTigramiteは、各サンプルサイズの分布を事前に計算します。そこでは、各条件独立テストのための高価な計算である並べ替えのテストを避けます。GP回帰は、skleanでデフォルトのパラメータで実行されます。カーネルを除いて、それは、ここで、放射状の関数＋ホワイトカーネル(両方のハイパーパラメータは内部で最適化されます。)そしてノイズレベルalphaを仮定します。それは、私たちがホワイトカーネルを追加しているのでゼロに設定されます。これらとその他のパラメータは、gp_params辞書を経由して、設定することができます。より詳しい議論はskleanのドキュメントを参照していください。また、GPUで高速に計算するためにgpytorchを利用するモジュール(gpdc_torch.py)があります。

gpdc = GPDC(significance='analytic', gp_params=None)
pcmci_gpdc = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=gpdc,
    verbosity=0)

　ParCorrと対照的に、非線形リンクは、GPDCとともに正しく検出されます。

results = pcmci_gpdc.run_pcmci(tau_max=2, pc_alpha=0.1, alpha_level = 0.01)
tp.plot_graph(
    val_matrix=results['val_matrix'],
    graph=results['graph'],
    var_names=var_names,
    show_colorbar=False,
    )

　少し脱線して、私たちは、スキャッタープロットを見ることによって、いかにGPDCが働くかわかります。

array, dymmy, dummy, dummy = gpdc._get_array(X=[(0, -1)], Y=[(2, 0)], Z=[(1, -2)], tau_max=2)
x, meanx = gpdc._get_single_residuals(array, target_var=0, return_means=True)
y, meany = gpdc._get_single_residuals(array, target_var=1, return_means=True)

fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(8,3))
axes[0].scatter(array[2], array[0], color='grey', alpha=0.3)
axes[0].scatter(array[2], meanx, color='black')
axes[0].set_title("GP of %s on %s" % (var_names[0], var_names[1]) )
axes[0].set_xlabel(var_names[1]); axes[0].set_ylabel(var_names[0])
axes[1].scatter(array[2], array[1], color='grey', alpha=0.3)
axes[1].scatter(array[2], meany, color='black')
axes[1].set_title("GP of %s on %s" % (var_names[2], var_names[1]) )
axes[1].set_xlabel(var_names[1]); axes[1].set_ylabel(var_names[2])
axes[2].scatter(x, y, color='red', alpha=0.3)
axes[2].set_title("DC of residuals:" "\n val=%.3f / p-val=%.3f" % (gpdc.run_test(
            X=[(0, -1)], Y=[(2, 0)], Z=[(1, -2)], tau_max=2)) )
axes[2].set_xlabel("resid. "+var_names[0]); axes[2].set_ylabel("resid. "+var_names[2])
plt.tight_layout()

ノイズを乗数したモデルのもっと非線形な従属性をいくつか見てみましょう。

random_state = np.random.default_rng(seed=11)
data = random_state.standard_normal((600, 3))
for t in range(1, 600):
    data[t, 0] *= 0.2*data[t-1, 1]
    data[t, 2] *= 0.3*data[t-2, 1]
dataframe = pp.DataFrame(data, var_names=var_names)
tp.plot_timeseries(dataframe); plt.show()

　乗数したノイズは、GPDC下で加法従属性の仮定に違反するので、制約外にできないために、不確かなリンクX_t-1⁰ -> X_t²は、間違って検出します。これと対照的にParCorrは、しかし、DCが何らかの従属性を検出するために、二つの真のリンクが検出されます。(skleanのガウス過程関数の元から警告があります)

pcmci_gpdc = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=gpdc)
results = pcmci_gpdc.run_pcmci(tau_max=2, pc_alpha=0.1, alpha_level = 0.01)
tp.plot_graph(
    val_matrix=results['val_matrix'],
    graph=results['graph'],
    var_names=var_names,
    show_colorbar=False,
    )

　ここで、私たちは、スキャッタープロットを見ることができます。ガウス仮定は、従属性と剰余を適合できないので、このように独立でなくなります。

array, dymmy, dummy, dummy = gpdc._get_array(X=[(0, -1)], Y=[(2, 0)], Z=[(1, -2)], tau_max=2)
x, meanx = gpdc._get_single_residuals(array, target_var=0, return_means=True)
y, meany = gpdc._get_single_residuals(array, target_var=1, return_means=True)

fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(8,3))
axes[0].scatter(array[2], array[0], color='grey')
axes[0].scatter(array[2], meanx, color='black')
axes[0].set_title("GP of %s on %s" % (var_names[0], var_names[1]) )
axes[0].set_xlabel(var_names[1]); axes[0].set_ylabel(var_names[0])
axes[1].scatter(array[2], array[1], color='grey')
axes[1].scatter(array[2], meany, color='black')
axes[1].set_title("GP of %s on %s" % (var_names[2], var_names[1]) )
axes[1].set_xlabel(var_names[1]); axes[1].set_ylabel(var_names[2])
axes[2].scatter(x, y, color='red', alpha=0.3)
axes[2].set_title("DC of residuals:" "\n val=%.3f / p-val=%.3f" % (gpdc.run_test(
            X=[(0, -1)], Y=[(2, 0)], Z=[(1, -2)], tau_max=2)) )
axes[2].set_xlabel("resid. "+var_names[0]); axes[2].set_ylabel("resid. "+var_names[2])
plt.tight_layout()

1.4 CMIknn

　Tigramiteのもっとも一般的な条件独立テストの実装は、k-最近傍推定で推定された条件相互情報を基礎にしたCMIknnです。このテストは以下の論文で記述されます。

　Runge, Jakob 2018"Conditional Independence Testing Based on a Nearest-Neighbor Estimator of Conditional Mutual Information.” In Proceedings of the 21st International Conference on Artificial Intelligence and Statistics.

　CMIknnは、従属性に関する仮定を含みません。パラメータknnは、ハイパーキューブ、例えば、(データ適合)局所的長さスケールで決定されます。ここで、私たちは、CMIknnは、GPDCのように剰余を基礎にしておらず、null-分布はもっと多くの因子に依存するため、null分布で事前の計算すらできません。私たちは、個別のテストに置いて、それを生成するために、sigicance='shuffle_test'を使います。I(X;Y|Z)=0をシャッフルテストするために、Xの値を局所的にシャッフルします。各サンプルポイントi番目のx-値は、サブ空間Zのその近傍(shuffle_neighborsパラメータ)の一つへランダムにマップされます。そのほかの自由パラメータはtransformです。それは、データがCMI推定の前に、変換されるかどうかを特徴づけます。新しいデフォルト値、transform=ranksです。それは、古いtransform=standardizeよりもよく働きます。以下のセルは実行に数分かかります。

cmi_knn = CMIknn(significance='shuffle_test', knn=0.1, shuffle_neighbors=5, transform='ranks', sig_samples=200)
pcmci_cmi_knn = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=cmi_knn,
    verbosity=0)
results = pcmci_cmi_knn.run_pcmci(tau_max=2, pc_alpha=0.05, alpha_level = 0.01)

tp.plot_graph(
    val_matrix=results['val_matrix'],
    graph=results['graph'],
    var_names=var_names,
    link_colorbar_label='cross-MCI',
    node_colorbar_label='auto-MCI',
    vmin_edges=0.,
    vmax_edges = 0.3,
    edge_ticks=0.05,
    cmap_edges='OrRd',
    vmin_nodes=0,
    vmax_nodes=.5,
    node_ticks=.1,
    cmap_nodes='OrRd',
    ); plt.show()

　ここで、CMIknnは正確に真のリンクを検出します。そしてまた、不確かなリンクを明らかにします。CMIknnは、最高の独立テストの選択として見えるので、一般的に、従属性が実際に任意のパラメトリック様式に従うケースでは、コストにおいてもっと低いパワーから生まれることに注意しなければなりません。そのため、ParCorrまたはGPDCは、もっと強力な測定です。もちろん、ParCorrは線形のリンクをGPDCよりよく検出します。

　その他の注記：CMIknnは、シャッフルテストを実行するために、莫大な計算コストがかかります。あなたは、sig_samples(p-値の推定で大きなエラーのコスト)によってシャッフルサンプル数を減少させることができます。

　一方で、あなたは、任意の条件独立テストにおいて、オプションsignificance ='fixed_thres'を使い、閾値fixed_thres Iの選択するでしょう。その後、条件独立で帰無仮説テストが伝達されます。静的なI(X;Y|Z)のテストで条件独立の基準は、

  　　　　　　　　　　　　　 I(X;Y|Z) < I*

　I^＊は、そこで、ハイパーパラメータとみなされます。有限のサンプルの推定バイアスのせいで、I(X;Y|Z)の値が変数の次元に依存するため、CMIKnnのためのIの選択は、巧妙（微妙）です。良い点は、CMIknnはかなり高速です。

　このオプションだけが、条件独立テストが並べ替えのテストに関連することがわかります。CMIknnとCMIsymbがあります。

# Pick fixed_thres
fixed_thres = 0.05

cmi_knn = CMIknn(significance='fixed_thres', fixed_thres=fixed_thres, transform='ranks', knn=0.1)
pcmci_cmi_knn = PCMCI(dataframe=dataframe, cond_ind_test=cmi_knn, verbosity=0)
results = pcmci_cmi_knn.run_pcmciplus(tau_max=2)  # if fixed_thres is used, pc_alpha (and alpha_level) is ignored.

tp.plot_graph(
    val_matrix=results['val_matrix'],
    graph=results['graph'],
    var_names=var_names,
    link_colorbar_label='cross-MCI',
    node_colorbar_label='auto-MCI',
    vmin_edges=0.,
    vmax_edges = 0.3,
    edge_ticks=0.05,
    cmap_edges='OrRd',
    vmin_nodes=0,
    vmax_nodes=.5,
    node_ticks=.1,
    cmap_nodes='OrRd',
    ); plt.show()

2a.分類別/象徴的時系列

　分類別のまたは、象徴的(または離散)データは、値の間で順序と測度を持ちません。りんごと西洋梨を考えてください。そのような時系列を適応させるには、TigramiteはGsuaredとCMIsymb条件独立テストをを含んでいます。両方、離散値のヒストグラムから直接推定されます。null分布を得るために。Gsquaredは、X²を(ゼロエントリーに調整された自由度)を使います。それは漸近(データに依存)にだけ正しくなります。より小さなサンプルサイズはCMIsymbクラスを使います。それは、条件相互情報を基礎にした局所的な並べ替えテストを含みます。CMYsymbは、その時もっと多くの計算時間がかかります。

　下の過程、X⁰,X¹は、両方、二つのカテゴリー、X²は三つのカテゴリーを持ちます。X⁰,X²の確率は、交絡として振る舞うX¹のそれらに依存します。

T = 1000
def get_data(T, seed=1):
    random_state = np.random.default_rng(seed)

    data = random_state.binomial(n=1, p=0.4, size=(T, 3))
    for t in range(T):
        prob = 0.4+data[t-1, 1].squeeze()*0.2
        data[t, 0] = random_state.choice([0, 1], p=[prob, 1.-prob])
        prob = 0.4+data[t-2, 1].squeeze()*0.2
        data[t, 2] = random_state.choice([0, 1, 2], p=[prob, (1.-prob)/2., (1.-prob)/2.])
    return data
dataframe = pp.DataFrame(get_data(T), var_names=var_names)
tp.plot_timeseries(dataframe, figsize=(10,4)); plt.show()

　ここで、私たちは、Gsquaredに焦点を当てます。G2-静的値はサンプルサイズでスケールされることに注意してください。しかし、nがサンプルサイズのとき、G=2nCMIによって、条件相互情報に関連づけられます。このように私たちは、val_matrixをここで解釈できる(そして図示できる)量を得るために、変換します。

gsquared = Gsquared(significance='analytic')
pcmci_cmi_symb = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=gsquared)
results = pcmci_cmi_symb.run_pcmci(tau_min = 1, tau_max=2, pc_alpha=0.2, alpha_level = 0.01)

val_matrix = results['val_matrix']
val_matrix /= (2.*T)

tp.plot_graph(
    val_matrix=val_matrix,
    graph=results['graph'],
    var_names=var_names,
    link_colorbar_label='cross-MCI',
    node_colorbar_label='auto-MCI',
    vmin_edges=0.,
    vmax_edges = 0.3,
    edge_ticks=0.05,
    cmap_edges='OrRd',
    vmin_nodes=0,
    vmax_nodes=.5,
    node_ticks=.1,
    cmap_nodes='OrRd',
    ); plt.show()

　分けたチュートリアル、CMIsymbとまたCMIknnMixedを参照してくだい。それは、連続のカテゴリカル変数を処理することを議論しています。そして、変数は両方の種別を含んでいます。

　CMIsymbはここで、同じ結果を得るのにもっと長い時間かかります。しかし、より少ないサンプルサイズと高次の条件セットのために、Gsquaredのnull分布への漸近は、もはや、正しくありません。その時　CMIsymbが好まれます。

cmi_symb = CMIsymb(significance='shuffle_test')
pcmci_cmi_symb = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=cmi_symb)
results = pcmci_cmi_symb.run_pcmci(tau_min = 1, tau_max=2, pc_alpha=0.2, alpha_level = 0.01)

tp.plot_graph(
    val_matrix=val_matrix,
    graph=results['graph'],
    var_names=var_names,
    link_colorbar_label='cross-MCI',
    node_colorbar_label='auto-MCI',
    vmin_edges=0.,
    vmax_edges = 0.3,
    edge_ticks=0.05,
    cmap_edges='OrRd',
    vmin_nodes=0,
    vmax_nodes=.5,
    node_ticks=.1,
    cmap_nodes='OrRd',
    ); plt.show()

2b.混合分類別/連続時系列

　しばしば、いくつかの変数が分類別で、いくつかの変数が連続であるデータセットの状況に扱うことがあるでしょう。このケースは、RegressionCI独立テストで扱うことができます。その時、各変数の種別は、data_type引数によってデータフレームで設定されなければなりません。これは、ここでもっと一般的に、同じシェイプのバイナリデータ配列として実装されます。それは、変数の個別のサンプルが、連続か離散であるかどうか、連続であれば0、離散変数であれば、１。ここで、すべての変数のサンプルは、同じ種別です。

# Generate some mixed-type data with a binary variable (can also be multinomial) causing two continuous ones.
random_state = np.random.default_rng(42)
T = 1000
data = np.zeros((T, 3))
data[:, 1] = random_state.binomial(n=1, p=0.5, size=T)
for t in range(2, T):
    data[t, 0] = 0.5 * data[t-1, 0] + random_state.normal(0.2 + data[t-1, 1] * 0.6, 1)
    data[t, 2] = 0.4 * data[t-1, 2] + random_state.normal(0.2 + data[t-2, 1] * 0.6, 1)

data_type = np.zeros(data.shape, dtype='int')
# X0 is continuous, encoded as 0 in data_type
data_type[:,0] = 0
# X1 is discrete, encoded as 1 in data_type
data_type[:,1] = 1
# X2 is continuous, encoded as 0 in data_type
data_type[:,2] = 0

dataframe = pp.DataFrame(data,
                         data_type=data_type,
                         var_names=var_names)
tp.plot_timeseries(dataframe, figsize=(10,4)); plt.show()

regressionCI = RegressionCI(significance='analytic')

pcmci = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=regressionCI)
results = pcmci.run_pcmci(tau_min = 1, tau_max=2, pc_alpha=0.2, alpha_level = 0.01)

# Also here the test statistic values are not very useful
val_matrix = results['val_matrix']
val_matrix /= (2.*T)

tp.plot_graph(
    val_matrix=val_matrix,
    graph=results['graph'],
    var_names=var_names,
    link_colorbar_label='cross-MCI',
    node_colorbar_label='auto-MCI',
    vmin_edges=0.,
    vmax_edges = 0.1,
    edge_ticks=0.05,
    cmap_edges='OrRd',
    vmin_nodes=0,
    vmax_nodes=.5,
    node_ticks=.1,
    cmap_nodes='OrRd',
    ); plt.show()