tigramite教程（五）使用TIGRAMITE 进行自助聚合和链接置信度量化

使用TIGRAMITE 进行自助聚合和链接置信度量化

自助聚合（Bagging）和置信度估计
例子
- 数据生成模型
- 基本的PCMCI+
- Bagged-PCMCI+
- - 使用优化后的`pc_alpha`进行自举聚合
  - 使用优化的pc_alpha进行CMIknn的自举聚合

TIGRAMITE是一个用于时间序列分析的Python模块。它基于PCMCI框架，允许从离散或连续值时间序列中重建因果图模型，并创建结果的高质量图表。

这篇《自然-地球与环境》评论论文总结了一般时间序列因果推断的概况

本教程解释了时间序列因果发现的自助聚合（Bagging），该方法在PCMCIbase.run_bootstrap_of函数中实现。自助聚合是一种通用的元算法，可与TIGRAMITE的大多数因果发现方法结合使用，例如run_pcmci、run_pcalg_non_timeseries_data、run_pcmciplus、run_lpcmci等，包括整个条件独立性检验范围。您可以参考以下预印本获取更多信息。

import numpy as np

from matplotlib import pyplot as plt
%matplotlib inline     

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.lpcmci import LPCMCI
from tigramite.pcmci import PCMCI
from tigramite.independence_tests.parcorr import ParCorr
from tigramite.independence_tests.cmiknn import CMIknn
from tigramite.pcmci_base import PCMCIbase

自助聚合（Bagging）和置信度估计

在基于自助法的自助聚合（bagging）中，训练集中的随机样本是有放回地选择的，这意味着在每个自助样本中，每个数据点可以被多次抽取。以这种方式生成多个数据样本以产生一组复制品（也称为重新采样）。机器学习模型然后分别在每个复制品上进行训练，最终为预测任务对输出进行平均或为分类任务进行聚合（例如通过多数投票）。在我们的情况下，训练集是输入时间序列，机器学习模型是因果发现方法，输出是因果图。

自助聚合在时间序列因果发现中的主要兴趣在于改善输出图的稳健性，并为图中的连接提供置信度估计。由于时间序列因果发现对时间依赖性敏感，保持重采样过程中的时间依赖性是至关重要的。然而，标准重采样不可避免地会破坏（至少部分地）时间滞后依赖。为解决这个问题，我们采用以下所示的重采样策略：

（图片展示了重采样过程）

最终，我们的自助聚合方法结合时间序列因果发现方法（这里是PCMCI+）可以总结如下：原始时间序列被重采样 $B$ （自助法复制次数）次。在每个重采样中，PCMCI+被独立运行，产生一个输出图。然后通过相对多数投票对每对边的类型进行聚合，得到 $B$ 个输出图。对于每个边缘的类型，通过大多数投票对所有这些图进行聚合，得到最终输出图，解决连接冲突的优先级顺序（无连接，×−×和◦−◦）；如果→和←之间存在冲突的连接，则返回冲突连接×−×。

Bagged-PCMCI+的返回结果包括：

通过将PCMCI+应用于通过保留时间依赖性进行重采样获得的 $B$ 个数据集得到的 $B$ 个因果图集成，
将所有这些图通过多数投票（在每个单独边缘的级别）聚合到最终输出图中，
为最终聚合图提供连接频率，以提供连接的置信度量。

通过输出图中连接的宽度来表示此置信度量：时间滞后2时刻的 $X_2 \rightarrow X_3$ 的置信度量（与箭头宽度成比例）大于时间滞后1时刻的 $X_4 \rightarrow X_1$ 的置信度量。

在这里插入图片描述
与所选因果发现算法相比，被聚合的版本有一个进一步的参数：自助法复制次数 $B$ 。在实践中，我们建议尽可能使用较大的 $B$ 。该实施使用joblib进行并行化。在我们的数值实验中， $B = 25$ 已经获得了良好的结果，但我们建议使用 $B\geq 100$ 。

例子

这一部分演示和解释了Bagged-PCMCI+在合成数据上的应用。

数据生成模型

# Set seed for reproducibility
seed = 1111
# Choose the time series length
T = 100

# Specify the model (note that here, unlike in the typed equations, variables
# are indexed starting from 0)
def lin(x): return x

links = {0: [((0, -1), 0.3, lin), ((2, 0), 0.5, lin), ((3, -1), -0.5, lin)],   # X1
         1: [((1, -1), 0.3, lin)],                                             # X2
         2: [((2, -1), 0.3, lin), ((1, -2), 0.4, lin)],                        # X3
         3: [((3, -1), 0.3, lin)]                                              # X4                            
        }

var_names = [r'$X^{%d}$' % j for j in range(1, len(links)+1)]

# Show ground truth causal graph
tp.plot_graph(
    graph = PCMCI.get_graph_from_dict(links),
    var_names=var_names,
    )

在这里插入图片描述
现在往里面灌数据

# Generate data according to the full structural causal process
data, nonstationarity_indicator = toys.structural_causal_process(
    links=links, T=T, seed=seed)
assert not nonstationarity_indicator

# Number of variables
N = data.shape[1]

# Initialize dataframe object, specify variable names
dataframe = pp.DataFrame(data, var_names=var_names)

基本的PCMCI+

先展示以下基本的PCMCI+算法

tau_max = 2
pc_alpha = 0.01

pcmci = PCMCI(dataframe=dataframe,
            cond_ind_test=ParCorr(),
            verbosity=0,
            )
results_pcmciplus = pcmci.run_pcmciplus(tau_max=tau_max, pc_alpha=pc_alpha)
tp.plot_graph(
    graph = results_pcmciplus['graph'],
    val_matrix= results_pcmciplus['val_matrix'],
    var_names=dataframe.var_names,
    ); plt.show()

在这里插入图片描述
在这里，PCMCI+漏了连接 $X^2\to X^3$ ，并且无法确定连接 $X^3\to X^1$ 的方向。

Bagged-PCMCI+

从我们的数值实验中，我们发现Bagged-PCMCI+改善了邻接精度和方向识别率。然而，不能直接将Bagged-PCMCI+与相同的pc_alpha下的PCMCI+结果进行比较（请参见论文中的精度-召回曲线）。在这里，我们选择了一个更高的pc_alpha用于Bagged-PCMCI+。在下文中，我们将说明如何使用模型选择来选择pc_alpha。

Bagged-PCMCI+的另一个参数是自助法的块长度（默认为1）。可以选择性地用它来更好地处理自相关性，但其效果尚未评估。

pc_alpha_bootstrap = 0.1
boot_samples = 200

# The block-length of the bootstrap can optionally be used to better deal with autocorrelation, 
# but its effect was not yet evaluated.
boot_blocklength = 1

## Create PCMCI object to call run_bootstrap_of
pcmci = PCMCI(dataframe=dataframe,
        cond_ind_test=ParCorr(),
        verbosity=0,
        )

# Call bootstrap for the chosen method (here 'run_pcmciplus') and pass method arguments  
results = pcmci.run_bootstrap_of(
        method='run_pcmciplus', 
        method_args={'tau_max':tau_max, 'pc_alpha':pc_alpha_bootstrap}, 
        boot_samples=boot_samples,
        boot_blocklength=boot_blocklength,
        seed=123)

# Output graph, link frequencies (confidence measure), and mean test statistic values (val_mat)
boot_linkfreq = results['summary_results']['link_frequency']
boot_graph = results['summary_results']['most_frequent_links']
val_mat = results['summary_results']['val_matrix_mean']

# Plot
tp.plot_graph(
    graph = boot_graph,
    val_matrix= val_mat,
    link_width = boot_linkfreq,
    var_names=dataframe.var_names,
    ); plt.show()

在这里插入图片描述
在这个例子中，我们可以看到Bagged-PCMCI+获得了正确的图（在 $B$ 个重新采样中进行链接多数投票），并且还通过箭头的宽度提供了对链接的置信度的有用估计。

使用优化后的`pc_alpha`进行自举聚合

对于一系列算法和条件独立性测试，还可以将pc_alpha设置为指定的列表或为None（然后将使用列表[0.001, 0.005, 0.01, 0.025, 0.05]）。pc_alpha将根据在cond_ind_test.get_model_selection_criterion()中计算的得分在指定列表中的值上进行优化。这对于所有条件独立性测试都是不可能的。

这种方法也可以与自举聚合相结合。pc_alpha会在每次 $B$ 个自举重新采样中进行内部和个别优化。

pc_alpha_bootstrap = [0.001, 0.01, 0.05, 0.1, 0.2] # This can be adapted 
boot_samples = 200

# The block-length of the bootstrap can optionally be used to better deal with autocorrelation, 
# but its effect was not yet evaluated.
boot_blocklength = 1

## Create PCMCI object to call run_bootstrap_of
pcmci = PCMCI(dataframe=dataframe,
        cond_ind_test=ParCorr(),
        verbosity=0,
        )

# Call bootstrap for the chosen method (here 'run_pcmciplus') and pass method arguments  
results = pcmci.run_bootstrap_of(
        method='run_pcmciplus', 
        method_args={'tau_max':tau_max, 'pc_alpha':pc_alpha_bootstrap}, 
        boot_samples=boot_samples,
        boot_blocklength=boot_blocklength,
        seed=123)

# Output graph, link frequencies (confidence measure), and mean test statistic values (val_mat)
boot_linkfreq = results['summary_results']['link_frequency']
boot_graph = results['summary_results']['most_frequent_links']
val_mat = results['summary_results']['val_matrix_mean']

# Plot
tp.plot_graph(
    graph = boot_graph,
    val_matrix= val_mat,
    link_width = boot_linkfreq,
    var_names=dataframe.var_names,
    ); plt.show()

在这里插入图片描述

使用优化的pc_alpha进行CMIknn的自举聚合

最后，我们展示了通过优化pc_alpha和非线性条件独立性测试CMIknn来实现Bagged-PCMCI+方法。请参阅相应的教程。使用标准初始化CMIknn(significance='shuffle_test')将导致每个测试执行计算密集型的置换检验方案以获得零分布。另一种具有较少统计严谨性的替代方法是仅基于条件互信息上的固定阈值做出测试决策，使用CMIknn(significance='fixed_thres')，请参阅教程中的解释。这在这里进行了说明。

我们按以下方式创建非线性数据：

# Choose the time series length
T = 500

# Specify the model (note that here, unlike in the typed equations, variables
# are indexed starting from 0)
def lin(x): return x
def nonlin(x): return .2 * (x + 5. * x**2 * np.exp(-x**2 / 20.))

links = {0: [((0, -1), 0.3, lin), ((2, 0), 0.5, lin), ((3, -1), -0.7, nonlin)],   # X1
         1: [((1, -1), 0.3, lin)],                                             # X2
         2: [((2, -1), 0.3, lin), ((1, -2), 0.4, nonlin)],                        # X3
         3: [((3, -1), 0.3, lin)]                                              # X4                            
        }

# Generate data according to the full structural causal process
data, nonstationarity_indicator = toys.structural_causal_process(
    links=links, T=500, seed=seed)
# Initialize dataframe object, specify variable names
dataframe = pp.DataFrame(data, var_names=var_names)

# Use a range of fixed thresholds, these are used as pc_alpha (with a slight abuse of the parameter name)
# This can be adapted, higher thresholds lead to stricter link decisions and, hence, sparser graphs
fixed_thresholds = [0.01, 0.025, 0.05, 0.1]
boot_samples = 100

# The block-length of the bootstrap can optionally be used to better deal with autocorrelation, 
# but its effect was not yet evaluated.
boot_blocklength = 1

## Create PCMCI object to call run_bootstrap_of
pcmci = PCMCI(dataframe=dataframe,
        cond_ind_test=CMIknn(significance='fixed_thres', model_selection_folds=5),
        verbosity=0,
        )

# Call bootstrap for the chosen method (here 'run_pcmciplus') and pass method arguments  
results = pcmci.run_bootstrap_of(
        method='run_pcmciplus', 
        method_args={'tau_max':tau_max, 'pc_alpha':fixed_thresholds}, 
        boot_samples=boot_samples,
        boot_blocklength=boot_blocklength,
        seed=123)

# Output graph, link frequencies (confidence measure), and mean test statistic values (val_mat)
boot_linkfreq = results['summary_results']['link_frequency']
boot_graph = results['summary_results']['most_frequent_links']
val_mat = results['summary_results']['val_matrix_mean']

# Plot
tp.plot_graph(
    graph = boot_graph,
    val_matrix= val_mat,
    link_width = boot_linkfreq,
    var_names=dataframe.var_names,
    vmin_edges=0.,
    vmax_edges = 0.2,
    edge_ticks=0.05,
    cmap_edges='OrRd',
    vmin_nodes=0,
    vmax_nodes=.2,
    node_ticks=.1,
    cmap_nodes='OrRd',
    ); plt.show()