TIGRAMITE 是一个时间序列数据分析的python包，它基于PCMCI框架，可以从离散或连续值的时间序列中重建图形模型（条件独立性图），并创建高质量的结果图。

本教程通过演示示例来解释主要功能。内容包括：

1. 基本用法
2. 绘图
3. 整合（专家）关于链接的假设
4. 基准测试和验证
5. 因果效应估计
6. 数据集挑战
7. 滑动窗口分析
   有关理论背景，请参阅以下论文：
   Runge, Jakob. 2018. “Causal Network Reconstruction from Time Series: From Theoretical Assumptions to Practical Estimation.” Chaos: An Interdisciplinary Journal of Nonlinear Science 28 (7): 075310.

最后，以下的《Nature Review Earth and Environment》论文提供了关于时间序列因果推断的概述：

基本用法

# 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 sklearnimport tigramite
from tigramite import data_processing as pp
from tigramite.toymodels import structural_causal_processes as toysfrom tigramite import plotting as tp
from tigramite.pcmci import PCMCI
from tigramite.lpcmci import LPCMCIfrom 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

Tigramite提供了来自PCMCI框架的几种因果发现方法，可以在不同的假设集下使用。一个应用始终由一个方法和选择的条件独立性测试组成，例如PCMCI与ParCorr一起。以下两个表格概述了涉及的假设:

方法	假设	输出
	(除了因果马尔可夫条件和忠实性之外)
PCMCI	因果平稳性，没有同时发生的因果链接，没有隐藏变量	有向滞后链接，无向同时发生链接（对于tau_min=0）
PCMCIplus	因果平稳性，没有隐藏变量	有向滞后链接，有向和无向同时发生链接（时间序列CPDAG）
LPCMCI	因果平稳性	时间序列PAG
RPCMCI	没有同时发生的因果链接，没有隐藏变量	每个区制变量的区域变量和因果图，带有有向滞后联系，无向同时联系（对于tau_min=0）
J-PCMCI+	多个数据集，因果平稳性，没有隐藏系统混淆，除非与上下文相关	有向滞后链接，有向和无向同时发生链接（联合时间序列CPDAG）

条件独立性测试	假设
ParCorr	单变量，具有线性依赖和高斯噪声的连续变量
RobustParCorr	具有线性依赖，针对不同边际分布鲁棒的单变量连续变量
ParCorrWLS	具有线性依赖，可以考虑异方差数据的单变量连续变量
GPDC / GPDCtorch	具有加法依赖的单变量连续变量
CMIknn	具有更一般依赖关系的多变量连续变量（基于排列的测试）
Gsquared	单变量离散/分类变量
CMIsymb	多变量离散/分类变量（基于排列的测试）
RegressionCI	包含单变量离散/分类和（线性）连续变量的混合数据集
PairwiseMultCI	元条件独立性测试，将上述每个单变量测试转化为包括可能有助于增加效果大小的多变量测试的方法

这些方法的参考资料在方法的文档字符串和相应的教程中。

条件独立性测试的参考资料在相应的教程conditional_independence_tests中。接下来的步骤将引导您完成典型的因果分析过程。

简单玩玩

假设时间序列来自下面的公式
$$\begin{array}{rl}{X}_t^0& =0.7X_{t-1}^0-0.8X_{t-1}^1+{\eta }_t^0\\ X_t^1& =0.8X_{t-1}^1+0.8X_{t-1}^3+{\eta }_t^1\\ X_t^2& =0.5X_{t-1}^2+0.5X_{t-2}^1+0.6X_{t-3}^3+{\eta }_t^2\\ X_t^3& =0.7X_{t-1}^3+{\eta }_t^3\end{array}$$
其中 $\eta$ 是独立同分布的零均值单位方差随机变量。我们的目标是重建每个变量的驱动因素。在 Tigramite 中，可以使用函数 toys.structural_causal_process 生成此类过程。

seed = 42
np.random.seed(seed)     # Fix random seed
def lin_f(x): return x
links_coeffs = {0: [((0, -1), 0.7, lin_f), ((1, -1), -0.8, lin_f)],1: [((1, -1), 0.8, lin_f), ((3, -1), 0.8, lin_f)],2: [((2, -1), 0.5, lin_f), ((1, -2), 0.5, lin_f), ((3, -3), 0.6, lin_f)],3: [((3, -1), 0.4, lin_f)],}
T = 1000     # time series length
data, _ = toys.structural_causal_process(links_coeffs, T=T, seed=seed)
T, N = data.shape# Initialize dataframe object, specify time axis and variable names
var_names = [r'$X^0$', r'$X^1$', r'$X^2$', r'$X^3$']
dataframe = pp.DataFrame(data, datatime = {0:np.arange(len(data))}, var_names=var_names)

万年不变的第一步:画出来

使用tp.plot_timeseries进行画图

tp.plot_timeseries(dataframe); plt.show()

它是平稳的（可以进行额外检查）且不包含缺失值（在教程missing_masking中讲解）。

现在我们做出因果稳定性、没有隐藏变量以及仅存在滞后依赖的假设。因此，我们选择PCMCI作为因果发现方法。接下来，我们需要选择条件独立测试和超参数，例如最大时间滞后tau_max。来吧，我们一起调查数据。

调查数据依赖性和滞后函数

为了调查依赖性的类型，我们使用 plot_scatterplots 和 plot_densityplots 函数来查看依赖性是否真的是线性的。通过将 matrix_lags 参数设置为一个 (N, N) 整数 NumPy 数组，您可以选择每一对变量使用哪个滞后。在这里，它被设置为 None，意味着滞后为零。

parcorr = ParCorr(significance='analytic')
pcmci = PCMCI(dataframe=dataframe, cond_ind_test=parcorr,verbosity=1)
correlations = pcmci.get_lagged_dependencies(tau_max=20, val_only=True)['val_matrix']

输出

##
## Estimating lagged dependencies 
##Parameters:independence test = par_corr
tau_min = 0
tau_max = 20

matrix_lags = None #np.argmax(np.abs(correlations), axis=2)
tp.plot_scatterplots(dataframe=dataframe, add_scatterplot_args={'matrix_lags':matrix_lags}); plt.show()

（这里的对角面板显示了变量与其自身的零滞后散点图。）
接下来，我们调查联合密度和边际（对角面板）密度的核密度估计。

tp.plot_densityplots(dataframe=dataframe, add_densityplot_args={'matrix_lags':matrix_lags})
plt.show()

由于在我们这个小模型中，依赖关系看起来相当线性，而且分布是高斯分布，我们使用实现线性偏相关条件的独立性测试 ParCorr。当 significance='analytic’时，假设零分布是学生$t$分布。

接下来，一个好主意是绘制滞后无条件依赖关系，例如，使用 ParCorr 类绘制滞后相关性。这可以帮助确定在因果算法中选择哪个最大时间滞后 tau_max。
为此，我们使用 dataframe 和 ParCorr 作为 cond_ind_test 初始化 PCMCI 方法。

parcorr = ParCorr(significance='analytic')
pcmci = PCMCI(dataframe=dataframe, cond_ind_test=parcorr,verbosity=1)

correlations = pcmci.get_lagged_dependencies(tau_max=20, val_only=True)['val_matrix']
lag_func_matrix = tp.plot_lagfuncs(val_matrix=correlations, setup_args={'var_names':var_names, 'x_base':5, 'y_base':.5}); plt.show()

除了调查滞后为零的散点图和密度 above，您可以选择依赖关系具有最大绝对值的滞后。当然，您可能想要使用非线性条件独立性测试来评估具有最大依赖性的滞后。也就是说，使用pcmci.get_lagged_dependencies初始化非线性度量（例如，CMIknn或GPDC，如条件独立性测试教程中介绍）的PCMCI。

PCMCI 因果发现

在这里，由于上述滞后函数图中依赖关系在最大滞后大约8之后衰减，我们选择tau_max=8用于PCMCI。另一个主要参数是pc_alpha，它设置了条件选择步骤中的显著性水平。在这里，我们让PCMCI选择最优值，通过将其设置为pc_alpha=None。然后，PCMCI将使用Akaike信息准则在合理默认的值列表中（例如，pc_alpha = [0.05, 0.1, 0.2, 0.3, 0.4, 0.5]）对ParCorr情况进行优化。参数alpha_level=0.01表明我们在这个显著性水平上对结果的p值矩阵进行阈值处理，以获得图形。

pcmci.verbosity = 1
results = pcmci.run_pcmci(tau_max=8, pc_alpha=None, alpha_level=0.01)

##
## Step 1: PC1 algorithm for selecting lagged conditions
##Parameters:
independence test = par_corr
tau_min = 1
tau_max = 8
pc_alpha = [0.05, 0.1, 0.2, 0.3, 0.4, 0.5]
max_conds_dim = None
max_combinations = 1## Resulting lagged parent (super)sets:Variable $X^0$ has 5 link(s):[pc_alpha = 0.1]($X^0$ -1): max_pval = 0.00000, |min_val| =  0.816($X^1$ -1): max_pval = 0.00000, |min_val| =  0.729($X^3$ -4): max_pval = 0.04439, |min_val| =  0.064($X^2$ -5): max_pval = 0.06669, |min_val| =  0.059($X^3$ -1): max_pval = 0.08886, |min_val| =  0.054Variable $X^1$ has 2 link(s):[pc_alpha = 0.05]($X^1$ -1): max_pval = 0.00000, |min_val| =  0.700($X^3$ -1): max_pval = 0.00000, |min_val| =  0.522Variable $X^2$ has 3 link(s):[pc_alpha = 0.1]($X^2$ -1): max_pval = 0.00000, |min_val| =  0.450($X^1$ -2): max_pval = 0.00000, |min_val| =  0.429($X^3$ -3): max_pval = 0.00000, |min_val| =  0.171Variable $X^3$ has 3 link(s):[pc_alpha = 0.2]($X^3$ -1): max_pval = 0.00000, |min_val| =  0.350($X^1$ -3): max_pval = 0.14772, |min_val| =  0.046($X^3$ -8): max_pval = 0.18318, |min_val| =  0.042##
## Step 2: MCI algorithm
##Parameters:independence test = par_corr
tau_min = 0
tau_max = 8
max_conds_py = None
max_conds_px = None## Significant links at alpha = 0.01:Variable $X^0$ has 2 link(s):($X^1$ -1): pval = 0.00000 | val = -0.619($X^0$ -1): pval = 0.00000 | val =  0.559Variable $X^1$ has 2 link(s):($X^3$ -1): pval = 0.00000 | val =  0.638($X^1$ -1): pval = 0.00000 | val =  0.628Variable $X^2$ has 3 link(s):($X^3$ -3): pval = 0.00000 | val =  0.452($X^1$ -2): pval = 0.00000 | val =  0.429($X^2$ -1): pval = 0.00000 | val =  0.422Variable $X^3$ has 1 link(s):($X^3$ -1): pval = 0.00000 | val =  0.350

从输出中，您可以看到PCMCI为每个变量选择了不同的pc_alpha。run_pcmci的结果是一个包含p值矩阵、测试统计值矩阵（这里是MCI部分相关性）、可选的置信区间（在初始化ParCorr时可以指定），以及graph矩阵的字典。p_matrix和val_matrix的形状为(N, N, tau_max+1)，其中元素$(i,j,\tau )$表示对链接$X_{t-\tau }^i\to X_t^j$的测试。对于$\tau =0$的MCI值并不排除其他同时效应，只以过去的变量为条件。具有相同形状的graph数组是通过在指定的alpha_level处阈值化p_matrix获得的。它是一个字符串数组，通过–>表示显著的滞后因果链接，并通过o-o表示同时链接（其中方向无法通过PCMCI确定）。PCMCIplus方法也可以对同时链接进行定向。

注意： 测试统计值（例如，部分相关性）可以提供依赖关系强度的定性直觉，但对于正确的因果效应分析，请参考CausalEffects类和教程。

print("p-values")
print (results['p_matrix'].round(3))
print("MCI partial correlations")
print (results['val_matrix'].round(2))

错误发现率控制

如果我们想控制在这里进行的 $N^2{\tau }_{max}$ 次测试，我们可以进一步校正p值，例如，通过错误发现率（FDR）控制来得到q_matrix。然后，可以使用 调整 p_matrix和不同的alpha_level更新图形，使用get_graph_from_pmatrix()。

q_matrix = pcmci.get_corrected_pvalues(p_matrix=results['p_matrix'], tau_max=8, fdr_method='fdr_bh')
pcmci.print_significant_links(p_matrix = q_matrix,val_matrix = results['val_matrix'],alpha_level = 0.01)
graph = pcmci.get_graph_from_pmatrix(p_matrix=q_matrix, alpha_level=0.01, tau_min=0, tau_max=8, link_assumptions=None)
results['graph'] = graph

## Significant links at alpha = 0.01:Variable $X^0$ has 2 link(s):($X^1$ -1): pval = 0.00000 | val = -0.619($X^0$ -1): pval = 0.00000 | val =  0.559Variable $X^1$ has 2 link(s):($X^3$ -1): pval = 0.00000 | val =  0.638($X^1$ -1): pval = 0.00000 | val =  0.628Variable $X^2$ has 3 link(s):($X^3$ -3): pval = 0.00000 | val =  0.452($X^1$ -2): pval = 0.00000 | val =  0.429($X^2$ -1): pval = 0.00000 | val =  0.422Variable $X^3$ has 1 link(s):($X^3$ -1): pval = 0.00000 | val =  0.350

进一步相关的方法学教程

教程assumptions解释了因果发现的基础假设以及这些假设的违反如何影响方法。
在这里，我们展示了run_pcmci，它假设没有同时因果关系。请查看相应教程pcmciplus中的run_pcmciplus，该方法还可以检测（并解释由）同时因果关系引起的混杂。

教程latent-pcmci解释了类LPCMCI，如果您想允许除了同时关系之外的隐藏混杂，可以使用该类。

教程tigramite_tutorial_jpcmciplus展示了如何利用多个数据集学习联合因果图，并克服一些隐藏混杂的问题。

教程conditional_independence_tests给出了tigramite中可用的所有条件独立性测试的概述，包括非线性和分类变量的测试。

教程tigramite_tutorial_bootstrap_aggregation展示了如何使用引导法用于更稳健的因果图和链接的置信度估计。

教程pcmci_fullci比较了PCMCI与矢量自回归模型估计器的替代方法。

画图

Tigramite 提供了一些绘图选项：滞后函数矩阵（如上所示）、时间序列图和过程图，后者聚合了时间序列图中的信息。这两个选项都需要graph数组作为参数，可选地还需要val_matrix以及其他链接属性。

在过程图中，节点的颜色表示自回归最小条件影响（auto-MCI）值，链接的颜色表示交叉最小条件影响（cross-MCI）值。如果两个变量之间在多个滞后出现链接，那么链接的颜色表示最强的那个，并且标签列表按其强度顺序列出了所有显著的滞后。此外，设置show_autodependency_lags=True将在相应节点标签下方显示显著的自依赖滞后。

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',show_autodependency_lags=False); plt.show()

# Plot time series graph    
tp.plot_time_series_graph(figsize=(6, 4),val_matrix=results['val_matrix'],graph=results['graph'],var_names=var_names,link_colorbar_label='MCI',); plt.show()

尽管过程图看起来更美观，但时间序列图更能代表时空依赖结构，从中可以解读出因果路径。您可以使用 node_size 和 node_aspect 参数调整节点的大小和长宽比，还可以修改许多其他属性，具体请查看 plot_graph 和 plot_time_series_graph 的参数。
另外，链接也可以导出到一个 csv 文件中：

tp.write_csv(val_matrix=results['val_matrix'],graph=results['graph'],var_names=var_names,save_name='test_graph.csv',digits=5,
)

通常，我们可能具有关于链接的存在或缺失及其方向的先验知识。这种专家知识可以通过 link_assumptions 参数进行整合。

整合专家对链条的假设

link_assumptions 是一个形式为 {j:{(i, -tau): link_type, ...}, ...} 的字典，用于指定链接的假设。内部处理这个参数时，会将图初始化为 graph[i,j,tau] = link_type，对于字典条目 link_assumptions[j][(i, -tau)] = link_type。允许以下指定：

• graph[i,j,0] = '-->'：在滞后0时，从i到j存在一个有向链接
• graph[i,j,0] = '-?>'：让方法测试是否存在邻接，如果存在，则其方向是从i到j在滞后0时
• graph[i,j,0] = 'o-o'：在滞后0时，i和j之间存在邻接，但让方法确定其方向
• graph[i,j,0] = 'o?o'：让方法测试邻接及其方向（默认值）

链接假设需要保持一致性，即 graph[i,j,0] = '-->' 要求 graph[j,i,0] = '<--' 并且必须满足无环性。如果字典中没有链接，则假定该链接不存在。也就是说，如果 link_assumptions 不为 None，则所有链接必须指定，否则假定链接不存在。

在某些情况下，你可能只有关于可能链接中一小部分的先验知识，并且构建完整的嵌套字典很繁琐，其中背景知识的缺失需要通过 link_assumptions[j][(i, tau)] = '-?>' 来指定 tau < 0 的情况，以及通过 link_assumptions[j][(i, 0)] = 'o?o' 来指定 tau = 0 的情况。在这种情况下，你可以在代码或笔记本中使用静态方便函数 pcmci.build_link_assumptions(...)，该函数允许通过 link_assumptions_absent_link_means_no_knowledge 构建实现的链条假设（假定不存在的条目代表不存在链接），从而构建链接假设。

考虑以下示例。

links_coeffs = {0: [((0, -1), 0.7, lin_f)],1: [((1, -1), 0.7, lin_f), ((0, 0), 0.2, lin_f), ((2, -2), 0.2, lin_f)],2: [((2, -1), 0.9, lin_f)],}
T = 100     # time series length
data, _ = toys.structural_causal_process(links_coeffs, T=T, seed=8)
T, N = data.shape# Initialize dataframe object
dataframe = pp.DataFrame(data)
pcmci = PCMCI(dataframe=dataframe, cond_ind_test=ParCorr(),verbosity=0)

我们首先运行PCMCIplus（允许检测同时期链接）而不考虑链接假设（默认设置）。

tau_max = 2
pc_alpha = 0.05
link_assumptions = None
results = pcmci.run_pcmciplus(tau_max=tau_max, pc_alpha=pc_alpha,link_assumptions=link_assumptions,)
tp.plot_graph(val_matrix=results['val_matrix'],graph=results['graph'],); plt.show()

在这里，由于样本量小且链接强度弱，我们漏了从0到1（滞后0）和从2到1（滞后2）的两个链接。

现在我们假设在滞后0时从0到1的链接必须存在，但其方向尚未确定，而从2到1在滞后2时的链接（其方向由时间顺序确定）。此外，我们还假设在0和2之间不存在任何链接（在任何滞后情况下）：

link_assumptions = {j:{(i, -tau):'o?o' for i in range(N) for tau in range(tau_max + 1) if (i, -tau) != (j, 0)} for j in range(N)}# Exclude all links between 0 and 2
link_assumptions[0] = {(i, -tau):'o?o' for i in range(N) for tau in range(tau_max + 1) if ((i, -tau) != (0, 0)and i != 2)} 
link_assumptions[2] = {(i, -tau):'o?o' for i in range(N) for tau in range(tau_max + 1) if ((i, -tau) != (2, 0)and i != 0)} # Set link 1 o-o 0 at lag 0
link_assumptions[1][(0, 0)] = 'o-o'
link_assumptions[0][(1, 0)] = 'o-o'    # Required for consistency of contemporaneous links, would be internally added if not present
# Set link 2 --> 1 at lag 2
link_assumptions[1][(2, -2)] = '-->'results = pcmci.run_pcmciplus(tau_max=tau_max, pc_alpha=pc_alpha,link_assumptions=link_assumptions,)
print(results['graph'][:,:,0])
tp.plot_graph(val_matrix=results['val_matrix'],graph=results['graph'],); plt.show()

[['' '-->' '']['<--' '' '']['' '' '']]

这给出了正确的图形。请注意，我们没有指定链接1 <-- 0的方向，只指定了邻接性1 o-o 0。方向是由PCMCIplus在制表规则中推断出来的。请参阅PCMCIplus相应的教程。

基准测试和验证

建议您自己创建一个玩具模型数据集，该数据集包含与您的真实数据相同的挑战，但其中已知真实情况，这样您可以评估哪种方法效果最佳，并选择超参数。

请参阅相应文件夹中的教程。

因果效应估计

上述图表将链接的 因果强度 表示为相应独立性测试的测试统计值。然而，这不能直接解释为因果效应的度量。这由 CausalEffect 和 LinearMediation 类覆盖。

请参阅相应文件夹中的教程。

数据集挑战

最后，您可能会有多个数据集以及/或者诸如缺失或遮蔽等挑战，这些挑战由tigramite的DataFrame功能覆盖。
请参阅相应文件夹中的教程。

滑动窗口分析

教程sliding_window_analysis解释了函数PCMCI.run_sliding_window_of，这是一个方便的功能，允许在多变量时间序列的滑动窗口上运行所有PCMCI因果发现方法。