生物医学网络科学Biomedical Network Science
Example07 / 11

从时间序列推断有向边——Granger 因果、条件 GC 与转移熵（高斯情形）

所属专题
网络构建与统计验证
下载
Notebook · .py

模拟 3 变量一阶向量自回归 VAR(1)：Z 同时驱动 X 与 Y（Z 先于 X 一步、先于 Y 两步），X 与 Y 之间没有直接作用。 用最小二乘拟合受限/完整自回归模型，计算 GC = ln(受限残差方差 / 完整残差方差) 与 F 检验。 (1) 两两 GC 会给出虚假的 X→Y（来自共同驱动 Z 的时滞差）；条件 GC（把 Z 的过去也放入两个模型）消除它； (2) 对高斯变量，转移熵 TE = GC / 2（Barnett 等, 2009），这里用高斯熵公式直接计算 TE 来核对； (3) 时间降采样（每 k 步取一个点）削弱真实方向的 GC；只给“源”Z 加测量噪声，还会产生高度显著的反向 X→Z。 这正是 fMRI 中 GC 的主要难点（Seth 等, 2013）。p 值显示为 0 表示小于浮点下限。

Output输出图

从时间序列推断有向边——Granger 因果、条件 GC 与转移熵（高斯情形）：输出图 1

Output文字输出

两两 GC（不条件于第三个变量）：
  X→Y: GC=0.2385  F=  403.0  p=1.1e-155
  Y→X: GC=0.0046  F=    6.9  p=1.1e-03
  Z→X: GC=0.5678  F= 1143.8  p=0.0e+00
  Z→Y: GC=0.6058  F= 1246.1  p=0.0e+00
条件 GC（条件于另一变量的过去）：
  X→Y | Z: GC=0.0001  F=    0.1  p=0.9
  Y→X | Z: GC=0.0007  F=    1.1  p=0.34
  Z→X | Y: GC=0.5639  F= 1132.9  p=0
  Z→Y | X: GC=0.3674  F=  663.9  p=2.4e-239
Z→Y|X：GC=0.3674，高斯 TE=0.1837，GC/2=0.1837（二者相等）
真实方向 Z→X；比较条件 GC(Z→X|Y) 与 GC(X→Z|Y)：
  每 1 步取样、Z 的测量噪声 sd=0.0: Z→X GC=0.5832 (p=0)   X→Z GC=0.0001 (p=0.4)
  每 4 步取样、Z 的测量噪声 sd=0.0: Z→X GC=0.0941 (p=1e-102)   X→Z GC=0.0005 (p=0.3)
  每 1 步取样、Z 的测量噪声 sd=3.0: Z→X GC=0.0588 (p=9e-256)   X→Z GC=0.0053 (p=9e-24)
  每 4 步取样、Z 的测量噪声 sd=3.0: Z→X GC=0.0105 (p=4e-12)   X→Z GC=0.0002 (p=0.6)

Code代码

# 模拟 3 变量一阶向量自回归 VAR(1)：Z 同时驱动 X 与 Y（Z 先于 X 一步、先于 Y 两步），X 与 Y 之间没有直接作用。
# 用最小二乘拟合受限/完整自回归模型，计算 GC = ln(受限残差方差 / 完整残差方差) 与 F 检验。
# (1) 两两 GC 会给出虚假的 X→Y（来自共同驱动 Z 的时滞差）；条件 GC（把 Z 的过去也放入两个模型）消除它；
# (2) 对高斯变量，转移熵 TE = GC / 2（Barnett 等, 2009），这里用高斯熵公式直接计算 TE 来核对；
# (3) 时间降采样（每 k 步取一个点）削弱真实方向的 GC；只给“源”Z 加测量噪声，还会产生高度显著的反向 X→Z。
#     这正是 fMRI 中 GC 的主要难点（Seth 等, 2013）。p 值显示为 0 表示小于浮点下限。
import numpy as np
from scipy import stats

rng = np.random.default_rng(11)

def simulate(T, rng):
    Z = np.zeros(T); X = np.zeros(T); Y = np.zeros(T)
    for t in range(2, T):
        Z[t] = 0.6 * Z[t - 1] + rng.standard_normal()
        X[t] = 0.4 * X[t - 1] + 0.8 * Z[t - 1] + rng.standard_normal()
        Y[t] = 0.4 * Y[t - 1] + 0.8 * Z[t - 2] + rng.standard_normal()
    return np.c_[X, Y, Z]

def lagmat(D, p):
    """把 T×k 序列展开为滞后 1..p 的设计矩阵（加截距）。"""
    T = D.shape[0]
    cols = [D[p - l:T - l] for l in range(1, p + 1)]
    return np.column_stack([np.ones(T - p)] + cols)

def gc(data, src, tgt, cond=(), p=2):
    """src→tgt 的 Granger 因果（可条件于 cond 中的变量）；返回 GC 值、F 统计量与 p 值。"""
    y = data[p:, tgt]
    restricted = [tgt] + list(cond)
    full = restricted + [src]
    Xr, Xf = lagmat(data[:, restricted], p), lagmat(data[:, full], p)
    rss_r = np.sum((y - Xr @ np.linalg.lstsq(Xr, y, rcond=None)[0]) ** 2)
    rss_f = np.sum((y - Xf @ np.linalg.lstsq(Xf, y, rcond=None)[0]) ** 2)
    df1, df2 = p, len(y) - Xf.shape[1]
    F = ((rss_r - rss_f) / df1) / (rss_f / df2)
    return np.log(rss_r / rss_f), F, stats.f.sf(F, df1, df2)

def gaussian_te(data, src, tgt, cond=(), p=2):
    """高斯转移熵：TE = H(Y_t | Y_past, Z_past) - H(Y_t | Y_past, Z_past, X_past)，用条件方差计算。"""
    def cond_var(target_cols):
        y = data[p:, tgt]
        Xd = lagmat(data[:, target_cols], p)
        r = y - Xd @ np.linalg.lstsq(Xd, y, rcond=None)[0]
        return r.var()
    v_r = cond_var([tgt] + list(cond))
    v_f = cond_var([tgt] + list(cond) + [src])
    return 0.5 * np.log(v_r / v_f)        # 高斯熵 H = 0.5 ln(2πe σ²)，常数项相消

D = simulate(3000, rng)
names = ["X", "Y", "Z"]
print("两两 GC（不条件于第三个变量）：")
gc_pair, gc_cond, ds_res = [], [], []          # 保存结果用于作图
for s, t in [(0, 1), (1, 0), (2, 0), (2, 1)]:
    g, F, pv = gc(D, s, t)
    gc_pair.append(g)
    print(f"  {names[s]}→{names[t]}: GC={g:.4f}  F={F:7.1f}  p={pv:.1e}")
print("条件 GC（条件于另一变量的过去）：")
for s, t, c in [(0, 1, (2,)), (1, 0, (2,)), (2, 0, (1,)), (2, 1, (0,))]:
    g, F, pv = gc(D, s, t, cond=c)
    gc_cond.append(g)
    print(f"  {names[s]}→{names[t]} | {names[c[0]]}: GC={g:.4f}  F={F:7.1f}  p={pv:.2g}")

g = gc(D, 2, 1, cond=(0,))[0]
te = gaussian_te(D, 2, 1, cond=(0,))
print(f"Z→Y|X：GC={g:.4f}，高斯 TE={te:.4f}，GC/2={g / 2:.4f}（二者相等）")

# 降采样与测量噪声：真实方向 Z→X（X 不影响 Z），比较 GC(Z→X|Y) 与 GC(X→Z|Y) 及其 F 检验 p 值
print("真实方向 Z→X；比较条件 GC(Z→X|Y) 与 GC(X→Z|Y)：")
Dlong = simulate(20000, rng)
for k, noise_z in ((1, 0.0), (4, 0.0), (1, 3.0), (4, 3.0)):
    Dk = Dlong[::k].copy()
    Dk[:, 2] += noise_z * rng.standard_normal(len(Dk))      # 只给“源”Z 加测量噪声
    f, _, pf = gc(Dk, 2, 0, cond=(1,))
    b, _, pb = gc(Dk, 0, 2, cond=(1,))
    ds_res.append((k, noise_z, f, b))
    print(f"  每 {k} 步取样、Z 的测量噪声 sd={noise_z}: Z→X GC={f:.4f} (p={pf:.1g})   "
          f"X→Z GC={b:.4f} (p={pb:.1g})")

# 作图：两两 GC 与条件 GC；降采样与测量噪声下真实方向与反向 GC
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
x = np.arange(4)
ax[0].bar(x - 0.2, gc_pair, 0.4, label="pairwise GC")
ax[0].bar(x + 0.2, gc_cond, 0.4, label="conditional GC")
ax[0].set_xticks(x); ax[0].set_xticklabels(["X->Y", "Y->X", "Z->X", "Z->Y"])
ax[0].set_ylabel("Granger causality (log variance ratio)"); ax[0].legend(fontsize=8)
ax[0].set_title("Conditioning on Z removes the spurious X->Y", fontsize=10)
ax[1].bar(x - 0.2, [r[2] for r in ds_res], 0.4, label="Z->X | Y (true)")
ax[1].bar(x + 0.2, [r[3] for r in ds_res], 0.4, label="X->Z | Y (reverse)")
ax[1].set_xticks(x); ax[1].set_xticklabels([f"step {k}\nnoise sd {nz:g}" for k, nz, _, _ in ds_res])
ax[1].set_yscale("log"); ax[1].set_ylabel("conditional GC"); ax[1].legend(fontsize=8)
ax[1].set_title("Downsampling and source noise distort direction", fontsize=10)
plt.tight_layout(); plt.show()

Run在浏览器中运行

首次加载 Python 环境约需 20–40 秒；可以修改参数后重新运行。

从时间序列推断有向边——Granger 因果、条件 GC 与转移熵（高斯情形）