强度中心性与预期影响

# 在偏相关网络上计算强度 s_i = Σ_j |pcor_ij| 与预期影响 EI_i = Σ_j pcor_ij，并与相关网络上的强度比较。
# 输出解读：偏相关网络中疲乏、情绪低落的强度最高，与生成数据的真实结构一致；
# 本例真实网络没有负边，但未正则化的样本偏相关含有一些小的负估计值（抽样噪声），使部分症状的 EI 低于强度；
# 相关网络的强度排序前几位相同，中后段与偏相关网络不同。
# 中心性描述的是条件关联的多少，不说明干预该症状能改变其他症状。
import numpy as np
import matplotlib.pyplot as plt

names = ["兴趣减退", "情绪低落", "睡眠问题", "疲乏", "食欲改变",
         "自我评价低", "注意力差", "精神运动改变", "自杀意念"]
TRUE_EDGES = {(0, 1): .40, (1, 5): .30, (1, 8): .25, (5, 8): .25, (2, 3): .40,
              (3, 6): .30, (3, 7): .20, (3, 4): .20, (0, 3): .20, (6, 7): .20}

def simulate_symptoms(n, seed=0):
    K = np.eye(9)                        # 模拟数据，仅用于演示方法
    for (i, j), w in TRUE_EDGES.items():
        K[i, j] = K[j, i] = -w
    rng = np.random.default_rng(seed)
    return rng.multivariate_normal(np.zeros(9), np.linalg.inv(K), size=n)

def pcor(X):
    K = np.linalg.inv(np.cov(X, rowvar=False))
    d = np.sqrt(np.diag(K))
    P = -K / np.outer(d, d)
    np.fill_diagonal(P, 0)
    return P

X = simulate_symptoms(500)
P = pcor(X)
R = np.corrcoef(X, rowvar=False)
np.fill_diagonal(R, 0)

strength = np.abs(P).sum(axis=1)
ei = P.sum(axis=1)
s_corr = np.abs(R).sum(axis=1)
print("偏相关强度排序：", [names[i] for i in np.argsort(strength)[::-1]])
print("相关网络强度排序：", [names[i] for i in np.argsort(s_corr)[::-1]])
print("强度与预期影响的最大差值 =", round(float(np.max(strength - ei)), 3))

x = np.arange(9)
plt.bar(x - 0.2, strength, 0.4, label="strength (pcor)")
plt.bar(x + 0.2, s_corr / s_corr.max() * strength.max(), 0.4, label="strength (corr, rescaled)")
plt.xticks(x, [f"PHQ{i + 1}" for i in range(9)], rotation=45)
plt.legend()
plt.title("Strength centrality (simulated data)")
plt.show()