网络干预分析（NIA）——把"治疗分组"作为网络中的一个节点

# 模拟一项随机对照试验（模拟数据，仅用于演示方法）：治疗只直接改善"睡眠问题"，
# 睡眠改善再通过"疲乏"间接改善"情绪低落"。治疗后把分组变量（0/1）与症状放在同一网络中估计偏相关。
# 输出解读：分组与三个症状的简单相关都不为 0（总效应），但在偏相关网络中，
# 只有"分组–睡眠问题"有明显的边——条件于其他症状后，治疗对疲乏与情绪的关联消失，
# 提示其效应可能经由睡眠传递。这是 Blanken 等 (2019) NIA 的基本逻辑；
# 横截面的单个时间点不能确定先后，原研究用逐周重复测量来看效应如何在网络中展开。
import numpy as np

rng = np.random.default_rng(5)
n = 600
treat = np.repeat([0, 1], n // 2)
sleep = 1.0 - 0.8 * treat + rng.standard_normal(n)
fatigue = 0.6 * sleep + rng.standard_normal(n)
mood = 0.5 * fatigue + 0.2 * sleep + rng.standard_normal(n)
worry = 0.4 * mood + rng.standard_normal(n)
names = ["分组(治疗=1)", "睡眠问题", "疲乏", "情绪低落", "担忧"]
X = np.column_stack([treat, sleep, fatigue, mood, worry])

R = np.corrcoef(X, rowvar=False)
K = np.linalg.inv(np.cov(X, rowvar=False))
d = np.sqrt(np.diag(K))
P = -K / np.outer(d, d)
print(f"{'症状':<8}{'与分组的相关':>10}{'与分组的偏相关':>12}")
for j in range(1, 5):
    print(f"{names[j]:<8}{R[0, j]:>12.2f}{P[0, j]:>14.2f}")

import networkx as nx
import matplotlib.pyplot as plt

def draw_signed(ax, M, labels, thr, title, node_color="#DDDDDD", node_size=180, pos=None):
    """边宽 ∝ |权重|，蓝色为正边、红色为负边。"""
    n = len(labels)
    pos = nx.circular_layout(range(n)) if pos is None else pos
    G = nx.Graph(); G.add_nodes_from(range(n))
    G.add_weighted_edges_from((i, j, M[i, j]) for i in range(n) for j in range(i + 1, n) if abs(M[i, j]) > thr)
    ws = [d["weight"] for _, _, d in G.edges(data=True)]
    nx.draw_networkx_edges(G, pos, ax=ax, width=[max(1.0, 8 * abs(w)) for w in ws],
                           edge_color=["#4477AA" if w > 0 else "#CC6677" for w in ws])
    nx.draw_networkx_nodes(G, pos, ax=ax, node_size=node_size, node_color=node_color, edgecolors="k")
    lab_pos = {k: (1.22 * np.asarray(v)) for k, v in pos.items()}
    nx.draw_networkx_labels(G, lab_pos, labels=dict(enumerate(labels)), ax=ax, font_size=9)
    ax.plot([], [], color="#4477AA", lw=2, label="positive")
    ax.plot([], [], color="#CC6677", lw=2, label="negative")
    ax.legend(loc="lower right", fontsize=7, frameon=False)
    ax.set_title(title, fontsize=10); ax.set_aspect("equal"); ax.axis("off")
    ax.set_xlim(-1.45, 1.45); ax.set_ylim(-1.45, 1.45)

abbr = ["Treat", "Slp", "Fat", "Sad", "Wry"]
fig, axes = plt.subplots(1, 2, figsize=(10, 4.2))
x = np.arange(1, 5)
axes[0].bar(x - 0.2, R[0, 1:], 0.4, color="#BBBBBB", label="correlation with treatment")
axes[0].bar(x + 0.2, P[0, 1:], 0.4, color="#4477AA", label="partial correlation")
axes[0].axhline(0, color="k", lw=0.6)
axes[0].set_xticks(x, abbr[1:])
axes[0].set(ylabel="Association with treatment", title="Total vs direct association")
axes[0].legend(fontsize=8)
P0 = P.copy(); np.fill_diagonal(P0, 0)
draw_signed(axes[1], P0, abbr, 0.1, "Partial correlation network (|pcor| > 0.1)",
            node_color=["#EE8866"] + ["#DDDDDD"] * 4)
plt.tight_layout(); plt.show()