Example05 / 08
网络边检测的统计功效——Fisher z 近似与蒙特卡洛模拟
第一部分：用 Fisher z 变换近似计算“检出一条真实（偏）相关边”所需样本量， n ≈ ((z_{1-α/2} + z_{1-β}) / artanh(ρ))² + 3 + k，k 为偏相关中被控制的变量数。 第二部分：模拟一个 15 节点、稀疏的真实偏相关网络，对不同样本量重复抽样，用偏相关的 Fisher z 检验 （Bonferroni 校正）判定边，统计灵敏度（检出真实边的比例）和假阳性边数。 解读：边越弱、节点越多（多重比较越严、控制变量越多），所需样本量越大；把功效分析写进分析计划。
Output输出图
Output文字输出
ρ = 0.1: 单条边 α=0.05 需 n ≈ 783；105 条边 Bonferroni 且控制 13 个变量需 n ≈ 1884 ρ = 0.2: 单条边 α=0.05 需 n ≈ 194；105 条边 Bonferroni 且控制 13 个变量需 n ≈ 474 ρ = 0.3: 单条边 α=0.05 需 n ≈ 85；105 条边 Bonferroni 且控制 13 个变量需 n ≈ 213 真实边数 17，真实偏相关范围 0.20–0.25 n = 100: 灵敏度 0.12，平均假阳性边 0.04 n = 200: 灵敏度 0.48，平均假阳性边 0.10 n = 400: 灵敏度 0.91，平均假阳性边 0.04 n = 800: 灵敏度 1.00，平均假阳性边 0.08 n = 1600: 灵敏度 1.00，平均假阳性边 0.02
Code代码
# 第一部分：用 Fisher z 变换近似计算“检出一条真实（偏）相关边”所需样本量，
# n ≈ ((z_{1-α/2} + z_{1-β}) / artanh(ρ))² + 3 + k，k 为偏相关中被控制的变量数。
# 第二部分：模拟一个 15 节点、稀疏的真实偏相关网络，对不同样本量重复抽样，用偏相关的 Fisher z 检验
# （Bonferroni 校正）判定边，统计灵敏度（检出真实边的比例）和假阳性边数。
# 解读：边越弱、节点越多（多重比较越严、控制变量越多），所需样本量越大；把功效分析写进分析计划。
import numpy as np
from scipy.stats import norm
def n_required(rho, alpha=0.05, power=0.80, k=0):
za, zb = norm.ppf(1 - alpha / 2), norm.ppf(power)
return int(np.ceil(((za + zb) / np.arctanh(rho)) ** 2 + 3 + k))
p = 15; m = p * (p - 1) // 2
for rho in [0.1, 0.2, 0.3]:
print(f"ρ = {rho}: 单条边 α=0.05 需 n ≈ {n_required(rho)}；"
f"{m} 条边 Bonferroni 且控制 {p-2} 个变量需 n ≈ {n_required(rho, alpha=0.05/m, k=p-2)}")
rng = np.random.default_rng(7)
# 构造稀疏精度矩阵：环状结构，偏相关约 0.2–0.3
K = np.eye(p)
for i in range(p):
j = (i + 1) % p; K[i, j] = K[j, i] = -0.25
for i, j in [(0, 7), (3, 11)]:
K[i, j] = K[j, i] = -0.2
Sigma = np.linalg.inv(K)
d = np.sqrt(np.diag(K)); true_pcor = -K / np.outer(d, d)
iu = np.triu_indices(p, 1); truth = np.abs(true_pcor[iu]) > 1e-8
print(f"真实边数 {truth.sum()}，真实偏相关范围 {np.abs(true_pcor[iu][truth]).min():.2f}–{np.abs(true_pcor[iu][truth]).max():.2f}")
crit = norm.ppf(1 - 0.05 / m / 2)
power_curve = [] # 保存模拟结果用于作图
for n in [100, 200, 400, 800, 1600]:
sens, fp = [], []
for _ in range(50):
X = rng.multivariate_normal(np.zeros(p), Sigma, size=n)
P = np.linalg.inv(np.cov(X.T)); dd = np.sqrt(np.diag(P))
r = (-P / np.outer(dd, dd))[iu]
z = np.arctanh(r) * np.sqrt(n - 3 - (p - 2))
detected = np.abs(z) > crit
sens.append(detected[truth].mean()); fp.append(detected[~truth].sum())
print(f"n = {n:5d}: 灵敏度 {np.mean(sens):.2f}，平均假阳性边 {np.mean(fp):.2f}")
power_curve.append((n, np.mean(sens), np.mean(fp)))
# 作图：解析的所需样本量曲线；模拟得到的灵敏度与假阳性边数随 n 的变化
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
rhos = np.linspace(0.08, 0.5, 40)
ax[0].plot(rhos, [n_required(r_) for r_ in rhos], "-", label="single edge, alpha = 0.05")
ax[0].plot(rhos, [n_required(r_, alpha=0.05 / m, k=p - 2) for r_ in rhos], "--",
label=f"Bonferroni over {m} edges, partial corr.")
ax[0].set_yscale("log"); ax[0].set_xlabel("true (partial) correlation rho"); ax[0].set_ylabel("n for 80% power")
ax[0].set_title("Required sample size (Fisher z)", fontsize=10); ax[0].legend(fontsize=8)
pc = np.array(power_curve)
ax[1].plot(pc[:, 0], pc[:, 1], "o-", label="sensitivity")
ax[1].set_xscale("log"); ax[1].set_xlabel("sample size n"); ax[1].set_ylabel("sensitivity"); ax[1].set_ylim(0, 1.05)
ax2 = ax[1].twinx(); ax2.plot(pc[:, 0], pc[:, 2], "s--", color="C1"); ax2.set_ylabel("false-positive edges (orange)")
ax[1].set_title("Simulated edge detection, 15-node network", fontsize=10)
plt.tight_layout(); plt.show()Run在浏览器中运行
首次加载 Python 环境约需 20–40 秒；可以修改参数后重新运行。