Detecting Causality in Complex Ecosystems（检测复杂生态系统中的因果关系附件）

作者：George Sugihara, Robert May, Hao Ye,1Chih-hao Hsieh, Ethan Deyle, Michael Fogarty, Stephan Munch
翻译：Wendy

摘要

  此附件包含：材料和方法；图S1-S5；GC计算S1-S5；表S1-S26;方框S1；视频S1-S3的图例；参考文献。

材料和方法

  Box S1给出了一个简单的例子来说明不可分离性，并提供了GC和Takens定理之间冲突的例子。它为Takens定理提供了一个代数说明。

  S1、S2和S3是简短的动画，用来说明本文中使用的SSR和动力系统概念。
符号：

变量	描述
φ	动态过程(时间流)在系统中随着时间“进化”状态。φ出现在某个d维流形M上(例如一个吸引子)。
m(t)	M上的向量表示系统在t时刻的状态。如果m (t)是M上的一个点，φ为流量，则m (t+1) = φ(m (t))。
M	流形(吸引子)由系统的所有轨迹和可能的状态x(t)组成。M是嵌入在E维状态空间(d≤E)中的d维流形，状态空间包含流形及其动力学，由系统的原始E笛卡尔坐标(基本变量)组成。
X	系统的观测函数(X通常被认为是φ下的“变量”)。X可能对应于包含M的实际E维状态空间的笛卡尔坐标，但更一般地是一个函数(例如原始E笛卡尔坐标的旋转或线性组合)，它将M中的点映射到一个观察点(即X: M→R)。
{X}	变量X随时间变化的值所对应的时间序列。因此，{X}可以被认为是动态的连续投影，在流形M上发生在一个特定的坐标轴X上，记录沿该维度沿时间的位移。
M_X	利用时间滞后{X}重构的影子。
x(t)	M_X上的向量对应于系统在t时刻的状态。
φX	动态过程φ的表示，但发生在重构阴影流形M_X上的点。
Ŷ(t)	MX

动态系统中因果关系的背景定义

  在这里，我们对在本工作中提出的因果标准提出了一个正式的阐述。这些思想是由Crutchfield首先指出的时延嵌入理论建立起来的，随后由Takens证明，随后推广。

  考虑一个动态过程φ定义在E维状态空间中点的时间演化，其轨迹收敛于d维(d≤E)流形M，使得φ: M→M。也就是说，如果m (t)是M上的一个点，则m (t+1) = φ(m (t))。

  设X为φ的观测函数。X通常被认为是一个系统变量，但通常是一个函数(例如原始E坐标的旋转或线性组合)，它将M中的点映射为一个实值标量(例如X: M→R)。每个X对应一个时间序列，{X} ={X(1)，…， X(L)}，追踪映射到实数序列的M中的点的轨迹，即X(t+1) = X(m (t))。时间段的长度(库大小)为L。

  一个滞后坐标嵌入使用E个{X}的滞后值作为坐标轴或维度来重建一个阴影吸引器流形MX。这个流形中的点，用x(t)表示，由e维向量集合组成:x(t) = < x(t)， x(t -τ)， x(t -2τ)，…， X(t-(E-1)τ)>，其中时间滞后τ为正。一般来说，MX上的点x(t)与M上的点m(t)映射为1:1，因此MX是原始吸引子流形M的微分重构。我们注意到有一些特殊的情况(下面讨论)，流形之间的映射是不对称的，而不是1:1的。

  如果由历史时间序列值({X}的子集)构造的流形MX能够熟练地预测样本外补体中的点，则可以确认嵌入的有效性。这测试吸引子是否充分展开，使轨迹不交叉(没有奇点)。由X的滞后坐标向量构造的流形MX在样本外可预测性显著时捕获了X的动态(通过误差统计量测量:MAE, RMSE，或观测值和预测值之间的相关性，ρ)。系统完全确定性时间序列观察无噪声的就会往完美的可预测性的L→∞时(即MAE~ 0和ρ~ 1)。在实际应用中,MX是一个近似,将显示通过观察误差收敛到水平集,过程噪声和数据长度。

  有一些特殊情况不满足流形重建的要求。例如，在典型洛伦兹吸引子中，坐标Z不会产生有效的阴影流形。因为叶的吸引子是对称的Z和的两个不动点吸引子(两叶的中心)有相同的Z值,影子歧管MZ只包含一个不动点,从而无法复制完整洛伦兹系统的动力学。然而，即使做一个小的旋转(例如，如果在Z上增加0.01 X)，新的MZ将与M微分同构。

收敛交叉映射(CCM)

  滞后坐标嵌入的一个一般性质是MX上的点x(t) 会1:1映射到m地图上的点m(t)， MX地图上的局部邻域到m地图上的局部邻域。此前,对于动态耦合的两个变量X和Y,当地社区重建各自的滞后,MX和MY,将映射到彼此，由于X和Y本质上是选择观察常见的吸引子流形M(图S1A)。聚合的交叉映射决定了MX上的局部邻域与MY上的局部邻域的对应程度。为了做到这一点，一个流形MX由变量X的滞后值构造而成，并用来估计Y的同时期值。因为MX是M的微分形，Y的估计在L趋于无穷时收敛(图S2)。在实际应用中，MX是一种近似，它会收敛到观测误差和处理噪声所设定的水平。因此，CCM的估计精度(或相关性)随着L的上升而上升，并达到一个平台。

CCM的一种算法

  我们的收敛交叉映射算法是基于单纯形投影。单纯形投影是一种最近邻算法，涉及到重建流形上附近点的指数加权距离来做核密度估计。

  考虑两个时间序列的长度,{X} = {X(1),(2),…,X (L)}和{Y} = {Y Y(1),(2),…,Y (L)}。我们首先形成滞后坐标向量x(t) = < x(t)， x(t -τ)， x(t -2τ)，…，x(t -1)τ)> 其中 t = 1+(E-1)τ 到 t = l。为了生成Y(t)的交叉映射估计，用Ŷ(t) | MX表示，我们首先在MX, x(t)上定位同时代的滞后坐标向量，并找到它的E+1最近的邻居。注意，E+1是在E维空间中建立一个包围单形所需的最小点数。接下来，表示x(t)的E+1最近的邻居的时间指数(从最近到最远)为t1，…tE+1。这些时间指数对应于MX上x(t)的最近邻居，用于识别Y(假定的邻居)中的点(邻居)，从E+1 Y(ti)值的局部加权平均值估计Y(t)。
$$\hat{Y}(t)\mid {\mathbf{M}}_{\mathbf{X}}=\sum w_{i}Y\left(t_{\mathrm{i}}\right)i=1\dots E+1$$
其中wi是一个基于x(t)和它在MX上的第i个最近邻居之间的距离的权重，Y(ti)是Y的同时期值。权重由
$$w_i=u_i/\sum u_j;j=1\dots E+1$$
其中
u i = exp ⁡ { − d [ x ‾ ( t ) , x ‾ ( t i ) ] / d [ x ‾ ( t ) , x ‾ ( t 1 ) ] } u_{i}=\exp \left\{-d\left[\underline{x}(t), \underline{x}\left(t_{\mathrm{i}}\right)\right] / d\left[\underline{x}(t), \underline{x}\left(t_{1}\right)\right]\right\} ui​=exp{−d[x​(t),x​(ti​)]/d[x​(t),x​(t1​)]}
d[x(s)， x(t)]是两个向量之间的欧氏距离。从Y到X的交叉映射也有类似的定义。

  如果X和Y是动态耦合的，MX的最近邻居应该识别MY上相应的最近邻居的时间指数。随着L的增加，吸引子流形填充，E+1最近的邻居之间的距离缩小。因此，Ŷ(t) | MX收敛于Y(t)和Xˆ(t) | MY收敛于X(t)。这样，我们使用最近邻居的收敛来测试MX上的状态和MY上的状态是否存在对应关系。

定向耦合

  双向因果关系意味着CCM在两个方向上，即MX交叉映射Y和MY交叉映射X(图3和图S1A)。这是Takens所讨论的一般1:1情况。注意,随着X对Y动力学的因果影响的增加，更多关于X动力学的信息被编码在我从Y的固定数量的观察中构建的流形中可以得到。相反,考虑一个变量X的动力学对Y,但反过来是不正确的。这可能发生，例如，如果变量X对应于外部强迫(图S1B)。

  面板（A）显示了通用双向情况的交叉映射，其中重构歧管MX与原始歧管和MY是同形的，因为变量X包含有关Y动力学的信息。 在时间t（表示为黑色圆圈）上，MX（橙色三角形）中最接近的点映射到MY。 在MY上，这些点的质心是目标预测（绿色圆圈），它给出Ŷ（t）|MX的值 。 灰色圆圈的值的范围对应于预测中的不确定性。 由于这些点在流形MY中保持接近（MX是微晶的），因此预测Ŷ（t）|MX的不确定性很低。 和Xˆ（t）| MY。

  在面板（B）中，预测与面板（A）相同，但是对于特殊的非一般情况，即X对Y的状态不敏感（即，Y对X不起任何作用）。 由于在这种情况下MX与M不等价，因此MX（橙色三角形）上的最近邻居无法准确映射到MY上的最近邻居，并且相距很远。 更具体地说，在MY上，由于X包含有关Y动力学的不完整信息，因此相应的按时间索引的邻居无法指定状态（在MY上随机分布）。 MX具有很高的不确定性。

  在这种情况下，MY是一个有效的流形它是原始流形M的微分形，我们期望用MY对变量X进行收敛的交叉映射。然而，这种特殊情况下，重构流形MX和M之间并没有1:1的映射，因为强制变量X不包含关于Y的动力学信息。因此，当L趋于无穷时，使用MX的变量Y的交叉映射在极限上不会收敛。

  使用MX对Y可能有显著的可预测性，这取决于条件概率，Pr[Y(t) | X(t)， X(t-τ)，…(X (t) - e 1τ)]。然而，使用MX的这种可预测性并不收敛(图4B)。只在一个方向收敛的交叉映射是单向因果关系的标准。而且，由于Y使用X的交叉映射并不收敛，这说明Y到X的信息流是不完整的，所以Y与X不是因果关系。

  请注意，本工作的重点是典型的情况下，强制变量的内在动力学仍然重要，因此强制变量的动力学反映的不仅仅是驱动系统。在这里，完整的流形可以被认为包括驱动系统和强制变量，而驱动系统反过来可以被认为是在完整系统的一个低维子流形上运行。然而，当被强迫变量的内在动力学被驱动子系统“奴役”时，就有可能将完整的流形分解为子流形。这种情况会发生在病态的强强迫下，因此被强迫变量的动态对其自身状态的依赖性不再显著。这描述了众所周知的“同步”现象。在这种情况下，整个系统折叠为驱动系统子歧管，而强制变量有效地成为驱动系统的一个观测函数。这里，与上面的不对称耦合情况不同，Takens定理将适用于此，以便收敛的交叉映射在两个方向上运作，因为两个变量对整个系统的动力学都很敏感。这种病理特殊的病例与双向病例很难区分。

时间序列的长度讨论

  一般来说，状态空间重构(SSR)方法在系统是非线性且可以在少数维近似时工作得最好，特别是在观测噪声不太大时。虽然已经有很多关于SSR在生态学中的应用的报道，但CCM的应用范围有多广还有待观察。尽管如此，令人鼓舞的是，它与时间序列数据一起工作，就像那些来自渔业的数据一样稀疏、有噪声和潜在的复杂。正如Ruelle所指出的，对于用于估计吸引子维数的重建，所需的数据量取决于吸引子维数，吸引子的分辨率大致随其维数的反幂变化。(请注意，CCM更关心“嵌入维数”(E)，而不是吸引子维数(d)，但这些都与Whitney嵌入定理(E≤2d+1)有关。)

  尽管我们无法预测该方法的广泛适用性，但对于渔业（在今年春季（2012年春季，在斯克里普斯海洋研究所研究非线性时间序列方法的NOAA / NMFS NLTS研讨会上），SIO-WHOI CAMEO小组发现这些方法有效 在所研究的大多数渔业案例（〜72种）中，数据嘈杂且时间序列短（35-40点）。 Hsieh等发现了相关的状态空间重构方法，在对20世纪相似的有限长度的海洋生物时间序列的广泛调查中很有效。

模拟和数据集

  这里我们提供了主要文本中分析的数据集的细节。章节a)提供了模型示例的完整规范，b)提供了草履虫-二碘时间序列和分析的详细信息，c)提供了沙丁鱼-凤尾鱼-温度数据的详细信息。

a) 模型描述

1. 耦合动态系统中的短暂状态(图1)
     图1显示了具有常数系数的耦合2种非线性Logistic差分系统单次运行的三个样本与海市蜃楼相关的现象：方程（1）rx = 3.8，ry = 3.5，β(x,y )= 0.02，β(y,x) = 0.1，起始条件X（1）= 0.4，Y（1）= 0.2。
2. 一般双向耦合动力系统(图2B)
     交叉映射的相对收敛速度指示了交互强度。 图2B使用1000个参数化总结了方程（1）的结果，其中r和β分别在间隔[3.6，4.0]和[0，0.4]中均匀选择。 对于每个系统，使用L = 400和E = 2的吸引子重构，针对1000个点计算交叉图估计。
3. 非对称耦合动力系统(图2C-D)
     图2C-D演示了CCM在非对称耦合动力系统中确定因果关系影响方向的用法：方程（1）中的β(x,y )= 0（rx = 3.7，ry = 3.7，β(y,x) = 0.32） 条件X（1）= 0.2，Y（1）= 0.4。 使用L = 1000和E = 2的吸引子重构绘制1000个交叉图估计值。
4. 例1:非耦合变量的外部强迫(图3B)
     我们分析了具有共同环境强迫的两个非相互作用种群的标准渔业模型系统。 人群遵循扩展的Schaffer模型，其中包括生殖延迟和成年生存。 成人（N）的招聘（R）包括白噪声环境变量T。
   $$\begin{array}{l}{R}_x(t+1)=X(t)[3.1(1-X(t))]\exp (-0.3T)\\ R_{\mathrm{y}}(t+1)=Y(t)[2.9(1-Y(t))]\exp (-0.36T)\\ X(t+1)=0.4X(t)+\max \left(R_x(t-3),0\right)\\ Y(t+1)=0.35Y(t)+\max \left(R_y(t-3),0\right)\end{array}$$
     我们比较了人口X（t）和Y（t）在不同时间序列长度L上的互相关和交叉映射。对于L的每个值，交叉映射可预测性或互相关值在500个随机选择的长度段上取平均值 来自8000点时间序列的L。 然后，对所有交叉映射计算都使用留一法交叉验证来估计同一段。 从X和Y的单变量SSR可以找到E = 4的最佳嵌入尺寸。
5. 例2:5种模型(图3C)
     下面的方程式为模型提供了显式的方程式系统。 CCM结果显示在图4的面板（C）中。交叉图技能是L = 300的观测值与估计值之间的相关性。300点的不同时间序列段用于吸引子重建。
     子系统Y1，Y2和Y3是强制子系统，它与Y4和Y5单向交互。 Y4和Y5不互相影响，也不影响Y1，Y2和Y3。
   $$\begin{array}{l}{Y}_1(t+1)=Y_1(t)\left[4-4Y_1(t)-2Y_2(t)-0.4Y_3(t)\right]\\ Y_2(t+1)=Y_2(t)\left[3.1-0.31Y_1(t)-3.1Y_2(t)-0.93Y_3(t)\right]\\ Y_3(t+1)=Y_3(t)\left[2.12+0.636Y_1(t)+0.636Y_2(t)-2.12Y_3(t)\right]\\ Y_4(t+1)=Y_4(t)\left[3.8-0.111Y_1(t)-0.011Y_2(t)+0.131Y_3(t)-3.8Y_4(t)\right]\\ Y_5(t+1)=Y_5(t)\left[4.1-0.082Y_1(t)-0.111Y_2(t)-0.125Y_3(t)-4.1Y_5(t)\right]\end{array}$$

b)二草履虫系统因果关系的经验检验

  我们使用经典的草履虫-二叠纪原生动物捕食者-捕食者系统的实验时间序列演示了CCM程序。 Veilleux在对Gause 和Luckinbill 的早期工作进行扩展的基础上，确定了导致掠食者（Didinium nasutum）和猎物（Paramecium aurelia）持续振荡的条件。 该方案涉及在培养皿中维持在27°C，pH值为6.7-6.8的培养皿，培养皿中含有6mL叶绿素浓度为0.9 g / L的培养基和3.5 g / L的甲基纤维素。 每两天更换50％的培养基。 草履虫和Did虫的初始密度分别为15和5个人/毫升。 通过充分混合容器并计数8个0.1 mL样品中的个体数，每12小时进行一次丰度测量。 这8个样品的平均值用于确定密度，误差<3.5％。 分析的数据显示在Jost和Ellner中，可以在以下位置找到http://robjhyndman.com/tsdldata/data/veilleux.dat.

  删除了前10个数据点，以消除实验初期的瞬态行为。 然后在分析之前将时间序列归一化为单位均值和方差。 最佳嵌入维数为E = 3。 为了说明有限的时间序列长度，将留一法交叉验证用于CCM分析。 也就是说，要在MX中找到x（t）的最近邻居，必须先从MX中删除x（t）和所有包含x（t）坐标的向量，以使映射保持样本外。 所有交叉图相关性均被截断为0。

c）对加利福尼亚现行沙丁鱼-An鱼-SST系统因果关系的实证检验

  在岸站Scripps码头和Newport码头测量了海面温度（SST）（数据可从以下网站获得：http：//shorestation.ucsd.edu/active/index_active.html）。 继Lasker和MacCall之后，对月均值进行平均以形成年度时间序列，最后的时间序列由SST的3年运行平均值组成。

  太平洋沙丁鱼（Sardinops sagax）和北部an鱼（Engraulis mordax）在加利福尼亚的着陆数据来自以下两个来源：

1. (1928 – 2002) NOAA Southwest Fisheries Science Center (http://las.pfeg.noaa.gov:8080/las_fish1/servlets/dataset?catitem=2)
2. (2003 – 2006) California Department of Fish and Game
   (http://www.dfg.ca.gov/marine/landings05.asp)

  为沙丁鱼和an鱼构造了1928年至2006年的着陆时间序列。 继Hsieh等在分析之前先对数据进行一阶差分处理。 CCM测试了沙丁鱼与an鱼，沙丁鱼与Scripps Pier SST之间的关联，以及an鱼与Newport Pier SST之间的关联。 所有CCM估算均使用嵌入维E = 3。 所有交叉图相关性均被截断为0。为解决有限的时间序列长度，分析采用了留一法交叉验证。

补充文本

经沙丁鱼-cho鱼-SST CCM测试的替代数据

  正文的图4中给出的结果基于着陆数据。 作为检查，我们检查了我们从其他来源获得的太平洋沙丁鱼和北部an鱼丰度数据分析的稳健性； 即独立于渔业的CalCOFI鱼鳞浮游生物调查和每日产蛋数据。

  根据Hsieh等人（1951年至2008年）建立了每年的鱼类浮游生物记录（1951-2008年）。并已被证明是成年产卵生物量的良好指标。 因为它们独立于渔业，所以它们消除了着陆数据中可能出现的潜在混杂因素。 使用与先前分析相同的方法，使用CCM测试沙丁鱼与an鱼，沙丁鱼与Scripps Pier SST之间的关系，以及an鱼与Newport Pier SST之间的关系。 尽管鱼鳞浮游生物时间序列中存在抽样缺口（调查始于1950年代，并且仅在1966-1984年每三年进行一次），但未在着陆数据中出现，但CCM结果仍然相似（图S3）。

  每日产蛋量（DEP）数据是成人生物量的第三种代表指标。 我们使用CCM分析来检验DEP作为an鱼丰度的替代物与SST之间的因果关系。 从图S4中可以看出，CCM结果清楚地证明了SST驱动着cho鱼种群的大小，而an鱼DEP并不影响SST。 在加利福尼亚州现行系统中，沙丁鱼的每日鸡蛋产量数据不可用。

  因此，在对加利福尼亚沙丁鱼-An鱼-SST系统的所有可用替代数据上重复进行CCM分析之后，我们发现结果与图5中描绘的相同因果关系网络一致。在此基础上，我们认为结果是可靠的。 因此，海面温度是加利福尼亚沙丁鱼和an鱼的常见因果变量，但这些物种之间似乎没有负相关关系，尽管它们之间没有因果关系。

用于检测因果关系的备用状态空间重构(SSR)技术

  基于SSR的思想，已经有几种尝试来推断不可分离系统的因果关系。但是，这些都倾向于集中在同步性上，并且几乎所有的人都假设格兰杰的可分离性范例。此外，没有人考虑收敛的关键特性，也没有提供检测弱耦合到中等耦合（或存在观察噪声的地方）因果关系的有效标准。由Arnhold等人提供的方案和Le Van Quyen等在推断因果关系的方向时均不正确。同样，Marinazzo等使用基于核的方法（类似于Sugihara（53）S-map）将Granger准则应用到非线性系统中。与格兰杰一样，这些作者假设在一个动态系统中，Y受X的影响，反之亦然，而不受X的影响，时间序列{X}将包含有关Y的信息（耦合程度低于同步）。如整个工作所示（特别是方框S1），在一般动态系统中应该颠倒这一标准。

  涉及由多个时间序列上的反复试验和错误搜索构成的嵌入的SSR方法对于系统识别（查找关键变量）很有用。 但是，它们无法将相关的非因果变量与真正的动态因果关系区分开，并且类似地，也可能无法区分因果关系的方向-特别是如果它们包含的驱动变量较弱或什至是中等强度时。

  Schiff等人仅在SSR研究中对因果关系做出正确的方向推断； 但是他们没有考虑收敛性，并且基于矢量的方法缺乏具有弱到中等交互强度的技能（特别是与生态系统相关）。 他们在图3A中说明了此缺点，其中他们的方法无法检测到单向耦合的Henon映射中的因果关系。 但是，如下表S1所示，CCM在此示例中表现良好（使用与（41）相同的模型，并在其图3A中给出了确切的参数）。

  模型中耦合Henon映射(C = 0.1, B = 0.3) CCM (E = 3, L = 4000)的有限收敛结果。无意义的(或消极的)ρŶ(t) | MX表明{X}缺乏信息的动态Y .因此,X是不会受到Y的影响 。然而,CCM的有限收敛性结果(L = 4000)清楚地表明,{Y}包含X,信息显示,Y是受到X假定值获得通过使用功率谱与1000年对比代理生成时间序列保存技术。

表S1

Cross mapping	ρ	p-value
Ŷ(t) / MX	-0.113	NaN
Xˆ(t) /MY	0.808	< 0.01

  最后，所有以前的基于SSR的方法都未能认识到收敛对于识别流形的重要性。 需要收敛来证明因果关系，而未能说明这一点已导致文献中的矛盾结果。

  为了说明这一点，考虑Le-Van-Quyen等人讨论的耦合Rossler-Lorenz系统:
Rossler system:
$$\begin{array}{l}dU/dt=-\alpha (V+W)\\ dV/dt=\alpha (U+0.2V)\\ dW/dt=\alpha (0.2+W(U–5.7)\end{array}$$
Lorenz system:
$$\begin{array}{l}dX/dt=10(-X+Y)\\ dY/dt=28X-Y-XZ+C{V}^2\\ dZ/dt=XY-8/3Z\end{array}$$
其中C是两个吸引子之间单向耦合的强度，而α控制驱动系统的特征时标。

  在不考虑收敛的情况下，Le-Van-Quyen等人错误地发现了驱动的Rossler系统对驱动的Lorenz 系统的更强的交叉预测。 然后他们误解了这个不正确的结果，以正确地得出Rossler系统是因果关系的结论。（令人困惑，但我们在下面解释）。

  我们使用收敛单纯形投影（CCM的预测类似物）（α= 6，如（42）中所示）重新分析该耦合系统。为了获得对初始起点不敏感的可靠结果，使用在每个L处随机选择的100个起点来分析时间序列。图S5中C = 2.0的整体结果在质量上类似于以下的耦合强度（C）同步。取决于起点，相对交叉可预测性可以随着从非常低的L值（较短的时间序列长度）到较大的L的变化而反转。在L的低值处，误差线重叠，并且相对交叉可预测性足够模棱两可，可以得出相反的结果。分析中的交叉可预测性。有趣（且令人困惑）的Le Van Quyen等基于他们的因果关系方向的结论基于他们的预测方法获得的相对交叉可预测性的（不正确）方向，并得出了“明显”正确的结果（将两个二进制误差相乘以得到正确的结果）。正确的解释是，Lorenz系统对Rossler系统的更强的“收敛”交叉预测表明，Rossler系统有因果关系迫使Lorenz系统，而不是（42）中给出的相反解释。图S5显示，仅靠交叉预测本身而不考虑收敛性不足以识别因果关系。

格兰杰因果关系

  如格兰杰所述，当将X从所有可能的因果变量U中去除时，如果目标变量Y的可预测性下降(以更高的方差表示)，则变量X被称为“格兰杰原因”Y。也就是说，当
$${\sigma }^2\left\{(Y|\bar{U})\right\}<{\sigma }^2\left\{(Y|\overline{U-X})\right\}$$
  （横线表示仅使用因果变量的先前值进行预测。）Granger的关键要求是可分离性，即可以从影响中减去有关因果因素的信息（如等式右侧所示 ）。 在线性系统中满足了可分离性要求，GC已用于随机系统，具有稳定点或极限环的线性系统，以及用于检测非线性系统中强耦合变量之间的相互作用。

  但是，正如Granger在早期意识到的那样，这种方法在一般的非线性动力系统中可能会出现问题（特别是在弱耦合到中等耦合的情况下）。例如，可以证明对于动态耦合变量（如图1所示），该方法会产生含糊的结果（有关详细分析，请参见下文）。这是因为一般而言，格兰杰的可分离性条件（Y | U-X）通常无法在动态系统中实现。也就是说，如果X是导致Y的原因，则有关X的完整信息将被编码在目标变量Y本身中，因此X不能从U上正式删除-这实际上是Takens定理的必然结果。因此，GC仅限于可以满足可分离性的情况。正如格兰杰本人所怀疑的，这不包括一般动力系统。的确，当将其应用于不可分离的系统（违反了Granger的定义）时，GC计算不再与因果关系有明确的联系：遗漏了在此类系统中检测因果关系的问题（方框S1包含一个简单的示例，说明了这些想法）。

格兰杰因果关系方法应用于本工作所使用的时间序列数据

  在这里，我们总结了将三种主要Granger方法（向量自回归（VAR），频谱分解（SD）和条件互信息（CMI））应用于本文中分析的三个模型示例和两个经验示例的结果（表S2） ）。 严格来说，由于在此类示例中不满足Granger对可分离性的假设，因此此类计算与Granger定义的因果关系没有必然联系。 （它们属于Granger职权范围之外的一类系统）。 因此，由于可分离性问题与计算GC因果关系的任何特定方法（例如VAR，光谱或非线性方法）无关，因此不保证使用详尽无遗的特定方法样本。 尽管如此，出于启发式原因，我们仍将这些计算包括在内，以说明将Granger的方法应用于原本非预期的不可分离系统中的危险。

表S2

  本文讨论的五个时间序列数据集的格兰杰因果检验汇总表。L是时间序列的长度。对于向量自回归(VAR)检验，最大滞后是通过似然比检验确定的。对于条件互信息(CMI)检验，采用修正的Baek & Brock方法，ε = 1.5, m = 1, LX = LY = 1，直到最大VAR滞后。对于光谱方法，我们检查了100个傅里叶区间，小波基参数为6(后(56))和15(后(57))，光谱外推极限为0.02。在α < 0.05的水平上对500迭代振幅适应傅里叶变换(IAAFT)替代(58)进行显著性检验。当整个因果网络被正确地描述时，结果被标记为正确的，如果测试错误地识别因果交互为不存在，或非交互为假阳性，则结果不正确。结果被标记为不一致，当使用多个参数设置，只有一些参数设置组合出现给出正确的因果网络。


  分析表明，当违反可分离性假设时，Granger检验并不稳健。 在这里考虑的每种情况下，它们要么始终错误地识别因果网络，要么给出对测试参数选择（例如，滞后长度，傅立叶间隔数，小波基参数）含糊不清的结果，或者给出随时间序列长度变化的结果。 再次，该分析不是对G因果关系方法的批评，因为这些示例不在Granger职权范围之内。

  我们注意到，在这些应用中的某些应用中，可能会找到测试参数的选择，这些参数似乎可以识别正确的因果网络。 但是，这样做必须事先知道网络的外观，以便选择产生正确网络的参数（任意拟合）。 因此，即使可分离性不会使这些检验无效（违反G因果关系的定义），多个结果的歧义性也只会使它们仅在事后才有用，而不能作为在未知答案时先验确定因果关系的一种方法。

GC 计算 S1-S5

描述

  以前的大部分关于因果关系检测的工作集中在随机和可分离系统(包括在频域内容易线性化或分解的系统)。这些方法遵循格兰杰的框架——如果一个变量提高了可预测性，就将其确定为因果关系。应用程序通常使用向量自回归(VAR)，条件互信息(CMI)，或光谱方法。这些方法非常适合于基于格兰杰因果关系的测试适用的线性和随机系统，但它们并不总是适合于分析一般的动态系统(Takens定理适用的系统)，特别是具有弱到中度耦合的不同组件的非线性系统。这里我们将详细描述每一种方法，并将它们应用到正文中讨论的示例中。我们注意到，由于这些例子不满足可分性，这些检验结果与格兰杰定义的因果关系之间没有必然的联系。尽管如此，我们相信，出于启发式的原因，将它们包括在内是有益的。因此，本节不是对技能的比较回顾，而是作为一个示例，涉及各种方法(包括最常用的方法)，将格兰杰的思想应用到不可分离系统。

GC计算S1.2种群模型示例

向量自回归

  协整的向量自回归是最常用的格兰杰因果检验。我们将此检验应用于由正文中动态系统(1)生成的时间序列，该系统由两个耦合logistic差分方程(图1)组成，初始条件(X(1) = 0.4, Y(1) = 0.2):
$$\begin{array}{l}X(t+1)=X(t)(3.8–3.8X(t)–0.02Y(t))\\ Y(t+1)=Y(t)(3.5–3.5Y(t)–0.1X(t))\end{array}$$
  尽管300点在生态学中是一个不现实的长时间序列(40-60点(年)的时间序列被认为很长)，我们在这里的目的是给格兰杰方法一个保守有利的演示。首先，我们通过一系列似然比检验(表S3)确定最优滞后长度，并使用Sims提出的自由修正因子。我们选择最大的滞后(在本例中为9)，对此，似然比检验报告所有短滞后的可预测性显著降低(p < 0.05水平):

表S3

  X和Y的时间序列不同滞后长度的似然比检验结果。因为χ2统计量对滞后10和滞后9不显著，但对所有较小的滞后显著，所以我们使用滞后长度9来分析该系统。

Test	χ2 statistic	p-value
lag 12 vs. lag 11	2.608	0.626
lag 11 vs. lag 10	4.458	0.348
lag 10 vs. lag 9	2.457	0.652
lag 9 vs. lag 8	12.655	0.013
lag 8 vs. lag 7	25.851	<0.001
lag 7 vs. lag 6	20.379	<0.001
lag 6 vs. lag 5	35.137	<0.001
lag 5 vs. lag 4	61.885	<0.001
lag 4 vs. lag 3	56.837	<0.001
lag 3 vs. lag 2	109.258	<0.001
lag 2 vs. lag 1	166.963	<0.001

  接下来，我们使用Johansen 的方法来确定协整关系(表S4)。由于时间序列在0和1之间自然有界，我们只测试常数项而不是时间序列中的高阶多项式。

表S4

  确定协整关系的Johansen程序的结果。检验同时使用迹统计量和最大特征值统计量来确定协整关系数r的显著性。

Null	test	test statistic	90% crit. value	95% crit. value	99% crit. value
r ≤ 0	trace	92.021	13.429	15.494	19.935
r ≤ 1	trace	30.506	2.705	3.841	6.635
r ≤ 0	eigenvalue	61.515	12.297	14.264	18.520
r ≤ 1	eigenvalue	30.506	2.705	3.841	6.635

  由于有显著的协整证据，我们使用误差修正模型(60)和格兰杰因果检验使用f检验。如前所述，模型包括每个变量的9个滞后。虽然测试(表S5)似乎确定因果关系从X到Y,他们不从X Y检测耦合我们相信这是由于系统的参数化,它使用一个增长率对物种Y (r = 3.5),不在范围混沌动力学~ (r > 3.57)。本质上，围绕稳定极限环的4个点进行线性化可以很好地捕捉Y的行为，然后X的滞后可以用来估计扰动。

  在后面的所有示例中，我们使用相同的过程来使用VAR来测试格兰杰因果关系，但我们只报告结果的摘要，并包括显示ρ值的最终表。

表S5

  结果采用向量自回归模型对300个数据的两种模型样本进行格兰杰因果检验。每个单元格包含检验行列“Granger原因”列变量的p值(小值更显著)。不显著的条目(0.05水平)由NaN表示。结果错误地指示了从X到Y的单向因果关系(缺少从Y到X的耦合)。

Variable	X	Y
X	<0.001	<0.001
Y	NaN	<0.001

  我们重复了上面的分析，但使用的是3000个数据点。有了这么多的数据，VAR分析发现了两个方向的因果关系的重要证据;然而，我们注意到，该程序不合理地将滞后值(高达23)确定为显著(基于似然比检验)。因此,VAR模型包含49个系数(23为每个变量系数,两种纠错方面由于显著的协整,和一个常数项),表明该系统是非常高维的,即使它只由两个耦合差分方程与一个滞后。这一发现并不令人惊讶，因为即使实际系统可以在几个维度上表征，我们也试图将一个线性模型拟合到一个具有混沌动力学的非线性系统中。然而，当一个明显不合适的模型与数据相匹配时，这引起了人们对解释结果的关注，以及获得正确答案所需的数据要求。

基于频率的格兰杰因果检验

  许多作者也通过使用傅立叶和小波变换方法分解时间序列，将格兰杰因果检验扩展到频域。在这里，我们使用100个傅里叶区间，小波基参数为6(后(56))和15(后(57))，以及光谱外推极限0.02来测试所描述的方法。在500个迭代振幅适应傅里叶变换(IAAFT)替代值的p < 0.05水平上检验显著性。

  以300个数据点和小波基参数为6,格兰杰因果关系检验是无意义的光谱频率在Y⇒X方向,和重要的65的100的频率在X⇒Y方向(表S6),建议使用小波基不对称耦合从X到Y参数15产生相似的结果(表S6)。这个结果与上面的VAR结果相似，并且错误地只检测到这个系统的一个方向的耦合。

表S6

  对300个数据点的两种模型样本进行光谱格兰杰因果检验结果。表项给出了针对500个IAAFT代理在0.05水平上显著的频率数。在这里(以及在以后的表中)使用NA来指示一个缺失的值，因为该测试不能用来指示变量是否影响自身。使用6或15作为小波基参数，结果错误地指示了从X到Y的单向因果关系(缺少从Y到X的耦合)。
小波基参数= 6：

Variable	X	Y
X	NA	0
Y	65	NA

小波基参数= 15:

Variable	X	Y
X	NA	0
Y	63	NA

  将分析扩展到3000个数据点，小波基参数为6的光谱格兰杰因果检验对100个X宗方向的频率中有46个显著，100个X宗方向的频率中有81个显著(表S7)。这个测试表明存在双向因果关系，尽管我们注意到有效频率的数量不一定对应于相对耦合强度。使用小波基参数15产生类似的结果(表S5)。因此，随着样本量的增加，明显的因果网络发生了变化(这应该令人不安)。同样，由于不可分离性，这种计算与“因果关系”本身的解释是不清楚的。

表S7

  对3000点数据的两种模型样本进行光谱格兰杰因果检验结果。表项给出了针对500个IAAFT代理在0.05水平上显著的频率数。
小波基参数= 6：

Variable	X	Y
X	NA	46
Y	81	NA

小波基参数= 15:

Variable	X	Y
X	NA	28
Y	69	NA

条件互信息

  利用条件互信息(CMI)的概念和等效的传递熵的概念来检测时间序列中的因果关系已经有大量的工作。在过去的二十年里，这些材料大多出现在神经科学和物理学文献中。尽管CMI显然是独立开发(引用经济学早期工作中基本上没有CMI文学),这些统计方法在数学上等价于早些时候门敏&布鲁克方法(经济学)的使用相关积分和条件概率来确定因果关系在格兰杰的框架。

  与SSR预测方法类似，这些方法基于过去的行为可以作为时间序列未来行为的指标的思想。考虑一个由连续的时间滞后Y构造的E维向量:Y: y(t, E) = <Y(t), Y(t-1), …, Y(t-(E-1))>。如果我们有过去观察向量,类似于相应的向量为一个特定的时间,这些向量的前进轨迹可能包含有用的信息对未来值Y: Pr[|Y(t+1) – Y(s+1)| < ε | ||y(t, LY) – y(s, LY)|| < ε] > Pr[|Y(t+1) – Y(s+1)| < ε] 对于某个ε > 0和LY≥1。这里，||x||是最大范数，尽管任何距离度量都可以使用。从本质上讲，过去轨迹接近(||y(t, LY) - y(s, LY)|| < ε)的事实，提高了未来轨迹接近(| y(t +1) - y(s +1)| < ε)的可能性。

  将这一观点扩展到格兰杰因果关系，Baek & Brock检验，与单独条件作用于Y相比，条件作用于另一个时间序列的过去值(如X)是否能够提高Y未来值的可预测性。因此，如果X的历史值提高了对Y的预测，那么X格兰杰导致Y。对于给定的LX和LY≥1，和ε > 0, X格兰杰导致Y如果
$$Pr[|Y(t+1)–Y(s+1)|<\epsilon |||x(t,LX)–x(s,LX)||<\epsilon ,||y(t,LY)–y(s,LY)||<\epsilon ]>Pr[|Y(t+1)–Y(s+1)|<\epsilon |||y(t,LY)–y(s,LY)||<\epsilon ].$$

  对于上述2种系统，此方法得出的结果不明确。 表S8给出了修改后的Baek＆Brock程序的结果，该程序系统地应用到上面测试的相同的300个点，是从正文中的系统生成的（图1）。

表S8

基于条件互信息的因果关系检验结果。

Direction	Lx	Ly	ε	test-statistic	p-value
X ⇒ Y	1	1	1.5	0.290	0.081
X ⇒ Y	2	2	1.5	0.432	0.024
X ⇒ Y	3	3	1.5	0.383	0.018
X ⇒ Y	4	4	1.5	0.445	0.007
X ⇒ Y	5	5	1.5	0.516	0.004
X ⇒ Y	6	6	1.5	0.531	0.003
X ⇒ Y	7	7	1.5	0.419	0.017
X ⇒ Y	8	8	1.5	0.484	0.012
X ⇒ Y	9	9	1.5	0.344	0.040
Y ⇒ X	1	1	1.5	0.086	0.106
Y ⇒ X	2	2	1.5	0.063	0.054
Y ⇒ X	3	3	1.5	0.012	0.107
Y ⇒ X	4	4	1.5	0.002	0.250
Y ⇒ X	5	5	1.5	0	NaN
Y ⇒ X	6	6	1.5	0	NaN
Y ⇒ X	7	7	1.5	0	NaN
Y ⇒ X	8	8	1.5	0	NaN
Y ⇒ X	9	9	1.5	0	NaN

  这种计算对参数LX、LY和ε的选择很敏感，并且没有关于如何先验选择这些参数的指南(除了事后选择最接近已知正确答案的参数)。在这种情况下,即使事后过程失败:一个系统的搜索在LX = LY ={1,2,…,9}和ε= 1.5建议Hiemstra &琼斯(55)发现没有证据表明在任何滞后Y到X的格兰杰因果关系,以及因果关系从X到Y仅落后的一些测试。实际上，在图1的主要文本中，人们并不期望这个测试能够成功地识别系统的因果关系。更重要的是，由于不可分离性，任何发现的关系都与格兰杰所定义的因果关系存在不确定关系。

  将分析扩展到3000个数据点，非线性格兰杰因果关系的检验在不同的滞后检验中仍然不一致。对于大多数滞后(LX = LY = 3…23)，从X到Y的因果关系显著，而只有少数滞后(LX = LY = 2…4)从Y到X的因果关系显著。换句话说，在测试的23个滞后中，只有2个在两个方向上正确地识别了因果关系。因此，当得到相互冲突的结果时，必须提前知道答案，以便选择(事后)给出正确答案的延迟时间。这说明了在不可分离系统中应用格兰杰方法的危险性，而这并不是它的本意。

GC计算S2.例1(由一个共同的驱动程序外部强制两个变量)

向量自回归

  至于2种系统,我们测试格兰杰因果在示例1(外部强迫两个变量的一个共同的司机——渔业模型例子讨论的主要文本)应用VAR。再一次,这代表了一类系统以外的分离系统,格兰杰能适当的地址。我们使用与CCM分析相同的协议(确定最大滞后，确定协整关系的数量，在适当的情况下使用修正错误的VAR，等等)。下表S9给出的结果(N=300时)表明VAR可以产生与真实因果交互网络相一致的结果。然而，当数据扩展到3000个数据点时(表S10)， VAR错误地发现Y2导致Y1，这是不正确的——增加时间序列长度导致了假阳性。我们注意到，我们只能在事后确定VAR结果的正确性，因为我们知道生成数据的模型。尽管如此，由于不稳定性，这些计算关于因果关系的意义最终是不确定的。

表S9

   用向量自回归模型对实例1的300个数据点进行格兰杰因果检验的结果。每个单元格表示检验行变量“Granger导致”列变量的p值(小值更显著)。结果正确地表明，变量Y1和Y2都受到白噪声(WN)信号的影响，彼此之间不相互作用。(注意，这里的条目和其他地方一样，都是影响列变量的行变量，NaN表示不重要)。

Variable	Y1	Y2	WN
Y1	<0.001	NaN	NaN
Y2	NaN	<0.001	NaN
WN	<0.001	<0.001	NaN

表S10

  使用向量自回归模型对3000点数据的例1进行格兰杰因果检验的结果。每个单元格表示检验行变量“Granger导致”列变量的p值(小值更显著)。结果错误地指出非交互变量Y1和Y2相互作用。

Variable	Y1	Y2	WN
Y1	<0.001	NaN	NaN
Y2	0.04	<0.001	NaN
WN	<0.001	<0.001	NaN

格兰杰因果关系的光谱检验

  我们使用Detto等人对Granger因果关系进行光谱检验重复上述分析。在300个点和小波基参数为6的情况下，光谱格兰杰因果检验对Y1引起Y2的100个频率中的23个，以及白噪声强迫对Y1和Y2的影响的大多数频率都是显著的(表S11)。这些结果表明Y1对Y2有一定的影响，这是不正确的。当使用小波基参数为15时，光谱格兰杰因果检验在区分白噪声信号对Y1和Y2的影响以及Y1对Y2的早期假阳性影响方面表现得更好(表S11)。

表S11

  实例1： 300个数据点的格兰杰光谱因果检验结果。表项给出了针对500个IAAFT代理在0.05水平上显著的频率数。在小波基参数为6时，结果错误地暗示了Y1到Y2之间的单向因果关系，但在其他方面是正确的。在小波基参数为15时，结果正确地识别出白噪声信号是系统中唯一的因果影响源。
小波基参数= 6：

Variable	Y1	Y2	WN
Y1	NA	23	1
Y2	0	NA	0
WN	64	89	NA

小波基参数= 15：

Variable	Y1	Y2	WN
Y1	NA	2	0
Y2	2	NA	0
WN	49	71	NA

  将分析扩展到3000个数据点，小波基参数为6的光谱格兰杰因果检验似乎增加了信号强度，同时也引入了假阳性(表S12)。特别地，分析似乎表明当Y1和Y2之间不存在耦合时。使用小波基参数15产生类似的假结果(表S12)。

表S12

小波基参数= 6：

Variable	Y1	Y2	WN
Y1	NA	25	5
Y2	38	NA	6
WN	96	96	NA

小波基参数= 15：

Variable	Y1	Y2	WN
Y1	NA	18	8
Y2	18	NA	0
WN	94	95	NA

  在这里，我们没有按照(57)的建议进行条件谱格兰杰检验，因为在实践中，条件变量不能预先知道(需要知道答案才能进行检验)。相比之下，CCM不需要这些信息。

格兰杰因果关系的非线性检验

  在此，我们使用修正的Baek & Brock(61)非线性Granger因果检验重复上述分析。(注意Baek & Brock是一种CMI测试)。对于每一对变量，我们测试了LX和LY的多个参数值(ε = 1.5, LX = LY = 1到5，因为5是VAR分析中似有比检验所支持的最大滞后)。表S13包含每个可能变量对的5种不同测试的最小ρ值，使用300个数据点的时间序列。当最小的p值(考虑到所有的参数组合)不正确地表明所有的变量是相互作用的，一些参数值导致了不显著的结果。表中用*标记的条目表示该交互既有显著结果，也有不显著结果。唯一一致的结果是Y2不影响Y1。然而，没有一个单一的滞后值可以正确地确定两个生物时间序列正受到白噪声的影响，同时也不包含提示Y1和Y2相互作用的假阳性，或者外部信号正受到Y1或Y2的影响。同样，这个测试的结果是模棱两可的。
   所有测试的滞后值(最多5个，与上面相同的方法确定)，表明Y1影响其外部驱动因素，外部驱动因素不影响Y1或Y2。所有其他的交互作用都不一致。

表S13

  例1(一个共同驱动因素对两个变量的外部强迫)中非线性Granger因果关系的Baek & Brock检验的结果(因果关系检验的最小p值超过了LX和LY参数的5个值)有300分。带有*的条目表示pvalue既低于0.05也高于0.05，这取决于参数设置。这一分析错误地识别了因果网络(它表明所有变量都是相互作用的)。

Variable	Y1	Y2	WN
Y1	NA	0.023*	0.002*
Y2	0.282	NA	<0.001*
WN	0.009*	0.022*	NA

表S14

  例1(由共同驱动因素对两个变量施加的外部压力)的非线性格兰杰因果关系的Baek & Brock检验结果(因果关系检验的最小p值超过了LX和LY参数的5个值)，有3000点。带有*的条目表示根据参数设置，p值均低于或高于0.05。这一分析错误地识别了因果网络(它表明所有变量都是相互作用的)。

Variable	Y1	Y2	WN
Y1	NA	0.012*	<0.001
Y2	<0.001*	NA	<0.001*
WN	0.148	0.078	NA

GC计算S3.例2(5种复杂模型)

向量自回归

  在例2(5种模型例)中，我们对CCM分析所用的相同长度的数据(L = 300)应用VAR检验格兰杰因果关系。下表S15的结果表明VAR没有正确识别因果交互网络。特别是，它表明Y2和Y3都没有内部动力学(因为它们不影响自己)，而y4的动力学独立于任何其他变量。

表S15

  因果关系检验结果(p值为格兰杰因果关系)基于向量自回归模型，使用300个数据点。结果表明，变量Y1、Y2、Y3是耦合的，但有一些环节缺失，而Y4不受其他变量的影响。测试还错误地指出，影响Y5的唯一其他变量是Y2。

Variable	Y1	Y2	Y3	Y4	Y5
Y1	<0.001	<0.001	<0.001	NaN	NaN
Y2	<0.001	NaN	NaN	NaN	<0.001
Y3	<0.001	<0.001	NaN	NaN	NaN
Y4	NaN	NaN	NaN	<0.001	NaN
Y5	NaN	NaN	NaN	NaN	<0.001

  我们试图将这种分析扩展到3000个点的时间序列，就像我们在前面的示例中所做的那样。但是，通过似然比检验，我们发现每个变量需要包含44个滞后因子，这导致了VAR模型的拟合出现了误差。使用更少的滞后(最大滞后= 10)使模型收敛，结果见表S16。然而，我们注意到10个滞后表明系统是非常高维的，因此不适合在这种情况下检测因果关系。结果与真实系统的匹配要好于只使用300个点(因为Y1、Y2和Y3是完全耦合的)，但是测试现在检测Y4和Y5之间的因果关系，这是不正确的。

表S16

  因果关系检验结果(p值为格兰杰因果关系)基于向量自回归模型，使用3000点。结果表明，变量Y1、Y2、Y3是完全耦合的，Y2对Y5是唯一的外部影响，Y5对Y4是唯一的外部影响。

Variable	Y1	Y2	Y3	Y4	Y5
Y1	<0.001	<0.001	<0.001	NaN	NaN
Y2	<0.001	0.02	0.04	NaN	0.00
Y3	<0.001	<0.001	<0.001	NaN	NaN
Y4	NaN	NaN	NaN	<0.001	NaN
Y5	NaN	NaN	NaN	0.05	<0.001

格兰杰因果关系的光谱检验

  我们使用Detto等人对Granger因果关系进行光谱检验，重复上述分析。在300个点，小波基参数为6的情况下，光谱格兰杰因果检验发现，在合理的频率比例下，Y1、Y2、Y3之间的大部分变量交互作用显著(表S17)。但是，检验也表明Y4对Y1和Y3有很强的因果影响，这是不正确的，说明没有变量影响Y4和Y5。当使用小波基参数15时，结果是相似的(表S17)。

表S17

  实例2 :300点数据的格兰杰光谱因果检验结果。表项给出了针对500个IAAFT代理在0.05水平上显著的频率数。用小波基参数6或15得到的结果是相似的，错误地认为Y1、Y2或Y3对Y4和Y5没有强迫作用，Y4影响Y1和Y3。
小波基参数= 6:

variable	Y1	Y2	Y3	Y4	Y5
Y1	NA	52	25	8	5
Y2	30	NA	30	0	0
Y3	8	30	NA	3	4
Y4	44	0	41	NA	0
Y5	0	0	0	5	NA

小波基参数= 15:

variable	Y1	Y2	Y3	Y4	Y5
Y1	NA	42	35	0	2
Y2	26	NA	41	1	2
Y3	26	36	NA	5	2
Y4	36	3	15	NA	2
Y5	6	3	4	11	NA

  将分析扩展到3000个数据点，小波基参数为6的光谱格兰杰因果检验(表S18)表明，只有很少的相互作用，主要是Y1、Y2和Y3之间的相互作用。虽然测试确实表明Y1对Y4有微弱的影响，但这与Y5对Y4和Y5对Y1的影响是相同的(不正确)。使用小波基参数15产生类似的结果(表S18)。光谱检验在提示Y2导致Y5方面做得更好，但是仍然忽略了Y1、Y2、Y3对Y4和Y5的大部分影响，除了(错误地)提示Y4对Y2、Y5对Y4可能有微弱的影响。

表S18

  例2 ：3000点数据的格兰杰光谱因果检验结果。表项给出了针对500个IAAFT代理在0.05水平上显著的频率数。在小波基参数为6的情况下，只有少数相互作用是显著的，很难区分Y2对Y3的影响是真实的还是Y5对Y4的影响是假阳性的。在小波基参数为15的情况下，Y1、Y2、Y3之间的强迫更强，但假阳性信号(Y4影响Y2或Y5影响Y4)也更强。
小波基参数= 6:

variable	Y1	Y2	Y3	Y4	Y5
Y1	NA	1	17	12	3
Y2	28	NA	10	1	2
Y3	1	21	NA	1	3
Y4	0	2	0	NA	0
Y5	9	0	0	11	NA

小波基参数= 15:

variable	Y1	Y2	Y3	Y4	Y5
Y1	NA	15	30	3	7
Y2	64	NA	36	6	29
Y3	46	24	NA	10	4
Y4	6	16	5	NA	2
Y5	4 1 4 15	NA

格兰杰因果关系的非线性检验

  在此，我们使用修正的Baek & Brock(61)非线性Granger因果检验重复上述分析。对于每一对变量，我们测试了LX和LY的多个参数值(ε = 1.5, LX = LY = 1到10，因为10是VAR分析中似有比检验所支持的最大滞后)。表S19包含每个可能变量对的10个不同测试的最小p值。而最小的p值表明了许多相互作用，一些参数值导致了不显著的结果。表中用a *标记的条目表明该变量交互作用的结果有显著性和非显著性。在10个不同的滞后长度中，每一对相互作用在某一点上都是显著的。此外，Y1、Y2、Y3之间的部分交互作用在各参数设置时均显著。因此，该测试在解析Y1、Y2和Y3之间相互作用的网络方面是部分成功的，但对于系统中较弱的耦合(Y1、Y2和Y3对Y4和Y5的影响)，不能充分区分假阳性和真相互作用。

  我们通过对来自同一系统的3000个点时间序列进行非线性格兰杰因果检验来扩展分析(表S19)。使用高达44的滞后(由相同数据的VAR分析的似然比检验确定)，所有的交互作用在某些参数设置上是显著的。在所有参数设置中唯一显著的交互作用是Y1、Y2和Y3之间的6个交互作用中的3个。因此，在系统中的12个真实的交互作用中，有3个被一致地识别出来，而在系统中的8个非交互作用中，所有8个都被错误地识别为一些滞后的显著交互作用。

表S19

  5种模型系统的Baek & Brock非线性Granger因果检验结果(因果检验的最小p值超过10个LX和LY参数值)，有300个点。与上述向量自回归分析的结果一样，该分析错误地识别了因果网络(唯一一致识别的因果关系是Y1、Y2和Y3之间，但可能的6个因果关系中只有4个是一致显著的)。

Variable	Y1	Y2	Y3	Y4	Y5
Y1	NA	<0.001	<0.001*	0.023*	0.008*
Y2	<0.001*	NA	<0.001	0.046*	0.025*
Y3	<0.001	0.000	NA	0.026	0.042*
Y4	0.038*	0.041*	0.040*	NA	<0.001*
Y5	<0.001*	<0.001*	<0.001*	<0.001*	NA

表S20

  5种模型系统的Baek & Brock非线性格兰杰因果检验结果(因果检验的最小p值超过44个LX和LY参数)，有3000点。所有的交互作用在某个滞后值是显著的，但12个实际连接中只有3个被一致识别。

Variable	Y1	Y2	Y3	Y4	Y5
Y1	NA	<0.001	<0.001*	<0.001*	<0.001*
Y2	<0.001*	NA	<0.001*	<0.001*	<0.001*
Y3	<0.001	<0.001	NA	<0.001*	<0.001*
Y4	<0.001*	<0.001*	<0.001*	NA	<0.001*
Y5	<0.001*	<0.001*	<0.001*	<0.001*	NA

GC计算S4.草履虫系统（Didinium-Paramecium system ）

向量自回归

  我们将VAR应用到用于CCM分析的相同实验数据集(17)上，在Didinium-Paramecium系统中检验格兰杰因果关系。表S21的结果表明VAR没有正确识别因果交互网络。事实上，它根本找不到任何相互作用!

表S21

  因果关系检验结果(p值为格兰杰因果关系)基于向量自回归模型。结果错误地表明草履虫与二碘之间不存在相互作用。

Variable	Paramecium	Didinium
Paramecium	NaN	NaN
Didinium	NaN	NaN

格兰杰因果关系的光谱检验

  我们使用Detto等人(57)对Granger因果关系进行光谱检验，重复上述分析。以300分6和小波基础参数,试验确定了影响草履虫的栉毛虫属(表S22),但发现意义只有6 100频率测试栉毛虫属对草履虫的影响时,大约在non-interactions上面所讨论的范围。使用小波基参数15，测试表明草履虫和二碘之间根本没有相互作用(表S22)。

表S22

  二氮草履虫捕食-食饵系统的光谱格兰杰因果检验结果。表项给出了针对500个IAAFT代理在0.05水平上显著的频率数。使用小波基参数6，结果表明该系统完全由自下而上的强迫控制，仅检测草履虫到二氮虫的单向因果关系。在小波基参数为15的情况下，草履虫与二碘之间不存在相互作用。
小波基参数= 6：

variable	Paramecium	Didinium
Paramecium	NA	22
Didinium	6	NA

小波基参数= 15：

variable	Paramecium	Didinium
Paramecium	NA	6
Didinium	1	NA

格兰杰因果关系的非线性检验

  在此，我们重复上述分析，使用修正的Baek & Brock(55)检验非线性格兰杰因果关系。对于每一对变量，我们测试了LX和LY的值，直到VAR分析中的似然比检验显示的最大滞后(在这种情况下，只使用了1个滞后)。其他参数与上述分析一致(ε = 1.5)。表S23包含了交互的pvalue。我们注意到，在这个例子中测试是成功的，并确定了草履虫和二碘之间的双向耦合。然而，与CCM分析不同的是，这种显著性度量不能量化耦合的不对称性。

表S23

  非线性Granger因果关系修正的Baek & Brock检验结果(行变量对列变量影响因果关系检验的p值)。该测试在识别两个变量之间的因果关系方面是成功的，但没有提供耦合强度中可能不对称的度量。

variable	Paramecium	Didinium
Paramecium	NaN	0.006
Didinium	0.004	NaN

GC计算S5.Sardine-anchovy系统

向量自回归

  我们将VAR应用于用于CCM分析的相同数据集(太平洋沙丁鱼和北方鳀鱼的着陆数据，以及在斯克里普斯码头和纽波特码头测量的海面温度)，在沙丁鱼-鳀鱼系统中测试格兰杰因果关系。下表S24给出的结果表明VAR没有正确识别因果交互网络。特别是，它未能确定生物变量中的任何动态，但错误地发现沙丁鱼影响了这两个海洋表面温度变量。

表S24

  因果关系检验结果(p值为格兰杰因果关系)基于向量自回归模型。结果奇怪地表明，沙丁鱼是海洋表面温度的驱动因素，而在这个系统中没有其他的动力。

Variable	Anchovy	Sardine	NP SST	SIO SST
Anchovy	NaN	NaN	NaN	NaN
Sardine	NaN	NaN	0.02	0.04
NP	SST	NaN	NaN	NaN
SIO	SST	NaN	NaN	NaN

格兰杰因果关系的光谱检验

  我们使用Detto等人对Granger因果关系进行光谱检验，重复上述分析。当小波基参数为6时，测试发现在沙丁鱼-凤尾鱼-海洋表面温度系统中没有因果相互作用(表S25)。使用小波基参数15，结果是相似的(表S25)，尽管测试表明了一些证据表明海表温度对沙丁鱼的影响。

表S6

  沙丁鱼-凤尾鱼-海洋表面温度系统的格兰杰光谱因果检验结果。表项给出了针对500个IAAFT代理在0.05水平上显著的频率数。使用小波基参数6，结果表明在这个系统中变量之间不存在因果相互作用。小波基参数为15时，结果表明，在这个系统中，海面温度可能对沙丁鱼有一定的影响，但变量之间没有其他的因果相互作用。
wavelet base parameter = 6 :

variable	Anchovy	Sardine	NP SST	SIO SST
Anchovy	NA	1	2	1
Sardine	1	NA	1	0
NP SST	1	1	NA	1
SIO SST	1	1	2	NA

wavelet base parameter = 15 :

variable	Anchovy	Sardine	NP SST	SIO SST
Anchovy	NA	4	2	1
Sardine	3	NA	1	1
NP SST	2	10	NA	7
SIO SST	1	8	4	NA

格兰杰因果关系的非线性检验

  在此，我们重复上述分析，使用修正的Baek & Brock(55)检验非线性格兰杰因果关系。对于每一对变量，我们测试了LX和LY的多个参数值(LX = LY = 1到2 ε = 1.5)。表S26包含了每个可能的变量对在这两种不同的测试上的最小p值。虽然Newport Pier SST对SIO Pier SST的影响是略微显著的(p < 0.1)，但在任何测试的滞后长度下都没有发现交互作用是显著的。

表S26

  非线性Granger因果关系的Baek & Brock检验结果(因果关系的最小p值大于LX和LY参数的2个值)。没有发现显著的相互作用。

variable	Anchovy	Sardine	NP SST	SIO SST
Anchovy	NA	0.956	0.903	0.857
Sardine	0.985	NA	0.664	0.721
NP SST	0.306	0.124	NA	0.088
SIO SST	0.357	0.259	0.139	NA

Box1.通过一个简单的例子构建不可分离性的直觉

  这项工作的一个基本想法是单方面的因果关系时,X⇒Y (X驱动Y以防ii),那么它可以估计X与Y,但没有从X ,Y这违背直觉(格兰杰因果关系),并建议如果天气使鱼类数量,我们可以用鱼来预测天气但不是亦然。请注意，CCM本身并不涉及预测，而是预测(估计)导致变量的当前或过去的状态。因此，如果鱼时间序列包含可以估计过去天气状态的历史信息(滞后)，那么如果明确地将天气添加到预测fish的模型中，那么这一信息(与fish相关的天气信息)将完全是多余的。因此，在格兰杰的模型中，天气(不正确)不会被视为因果关系，因为它可以从模型中添加或删除，而不会产生任何影响。不可分离性产生于冗余的因果信息已经完全包含在受影响的变量(Takens定理的一个结果)。

  为了进一步阐明这是如何工作的，考虑以下两种逻辑模型的特殊情况:
$$\begin{array}{l}X(t+1)=3.9X(t)[1–X(t)–\beta Y(t)](B1)\\ Y(t+1)=3.7Y(t)[1–Y(t)–0.2X(t)]\end{array}$$
参数β控制X到Y的映射，当β的变化的敏感性大,X和Y强烈降低增长，当β= 0,X是Y的独立(X的动力学不依赖于Y)。在这个简单的系统,我们可以确定lagged-coordinate代数模型,我们可以用Y的值(t + 1)和Y (t)的影响中恢复过来,因此价值,X (t)(反之亦然):
$$\begin{array}{l}\beta Y(t)=1–X(t)–X(t+1)/3.9X(t)(B2)\\ 0.2X(t)=1–Y(t)–Y(t+1)/3.7Y(t)\end{array}$$
  为了恢复动态交叉映射(请注意，仅用公式(B2)是不够的，因为它需要未来的知识!)，我们可以简单地将公式(B2)代入(B1)，得到X(t)用Y(t)和Y(t-1)表示的方程，反之亦然:

  注意，当β接近0时，X的交叉映射估计仍然表现良好。但对于Y有一个奇点，当β = 0时，Y的交叉映射模型爆炸。这可以直接从方程(B2)中得到，很明显，当β = 0时，没有办法用X表示Y。非正式地说，当β = 0时，X在状态空间中移动，而不考虑Y的当前位置。因此，X的历史与决定Y无关。

  此外，在双向因果(β≠0)的情况下，可以将(B3)代回(B1)，得到X(t)和X(t-1)的方程。因此,尽管X和Y是耦合的,我们可以写一个精确的模型来预测单X的基础上其历史(即不包括任何的观察Y)。然而,在GC框架中,因为Y可以从假设的因果变量的设置删除没有减少X的可预测性,我们将得出结论(错误),Y不会引起X。

  这个简单的模型还为为什么交叉映射估计的收敛性是检测因果关系的必要条件以及为什么收敛率提供了交互强度的指标提供了另一个论据。要看到这一点，请注意，对于一阶，(B3)表示Y中的方差与Var(X)/β成比例2。因为对Y的十字图估计是基于X的最近邻居，所以当X的邻居很大时，Y的不确定性就会很大。随着时间序列长度的增加，吸引子填充，最近邻域变小，Y中的方差减小。因此，Y中的方差是相互作用强度的函数;为了获得相同的精度，在β较小时比β较大时需要更多的数据。

参考文献

  该部分请看原文献