数学推导与算法

本节从多元正态分布出发推导“偏相关为 0 ⇔ 条件独立”，写出图形 LASSO 与 EBIC 的完整形式，给出 Ising 伪似然估计、病例删除 bootstrap 与 CS 系数、NCT 置换检验和 GVAR 三种网络的算法。符号：$p$ 个症状（节点），$n$ 名被试，$X\in\mathbb R^{n\times p}$，样本协方差（或相关）矩阵 $S$，总体协方差 $\Sigma$，精度矩阵 $K=\Sigma^{-1}$。偏相关等于精度矩阵归一化负值的回归推导（不需要正态假设）见脑连接组的数学推导。

A. 高斯图模型：为什么 $K_{ij}=0$ 意味着条件独立

设 $x\sim\mathcal N(0,\Sigma)$，密度为

$$f(x)\propto\exp\Big(-\tfrac12 x^\top Kx\Big)=\exp\Big(-\tfrac12\sum_{a}K_{aa}x_a^2-\sum_{a<b}K_{ab}x_ax_b\Big).$$

指数中含 $x_i$ 与 $x_j$ 的交叉项只有 $-K_{ij}x_ix_j$。若 $K_{ij}=0$，密度可以分解成“不含 $x_i$ 的因子”乘“不含 $x_j$ 的因子”（其余变量固定），由因子分解定理即得 $x_i\perp x_j\mid x_{\text{rest}}$。反之亦然。这就是 GGM 中“无边 = 条件独立”的来源，它依赖多元正态假设；对非正态数据，$K_{ij}=0$ 只保证“偏相关为 0”，不保证条件独立。

再看回归形式：把 $x_i$ 对其余变量回归，$x_i=\sum_{j\ne i}\beta_{ij}x_j+\varepsilon_i$，可以证明

$$\beta_{ij}=-\frac{K_{ij}}{K_{ii}},\qquad \operatorname{Var}(\varepsilon_i)=\frac{1}{K_{ii}},\qquad \rho_{ij\cdot\text{rest}}=\operatorname{sign}(\beta_{ij})\sqrt{\beta_{ij}\beta_{ji}}.$$

所以 GGM 也可以看成 $p$ 个相互约束的回归；“邻域选择”（逐节点 LASSO 回归）正是利用这一点的另一种估计方法。由第二式，对标准化变量（$\operatorname{Var}(x_i)=1$）有可预测性 $R^2_i=1-\operatorname{Var}(\varepsilon_i)=1-1/K_{ii}$。

B. 图形 LASSO 与 EBIC

高斯似然（忽略常数）为 $\ell(K)=\tfrac n2\big[\log\det K-\operatorname{tr}(SK)\big]$。图形 LASSO 求解

$$\hat K_\lambda=\arg\max_{K\succ0}\ \log\det K-\operatorname{tr}(SK)-\lambda\sum_{i\ne j}|K_{ij}|.$$

两个有用的事实：（1）当 $\lambda\ge\max_{i\ne j}|S_{ij}|$ 时解是对角阵，即空网络，这给出 $\lambda$ 路径的上端 $\lambda_{\max}$；（2）问题是凸的，常用块坐标下降求解（Friedman 等, 2008），每次更新一行/列相当于解一个 LASSO 回归。

块坐标下降的一步：记 $W=\hat K^{-1}$ 为当前协方差估计。把第 $j$ 行/列分出来，$W_{11}$ 为去掉 $j$ 后的子矩阵，$s_{12}$ 为 $S$ 的第 $j$ 列（去掉第 $j$ 个元素）。求解 LASSO 问题

$$\hat\beta=\arg\min_\beta\ \tfrac12\beta^\top W_{11}\beta-s_{12}^\top\beta+\lambda\|\beta\|_1 ,$$

用坐标下降：$\beta_k\leftarrow \operatorname{soft}\big(s_{12,k}-\sum_{l\ne k}(W_{11})_{kl}\beta_l,\ \lambda\big)/(W_{11})_{kk}$，其中 $\operatorname{soft}(u,\lambda)=\operatorname{sign}(u)\max(|u|-\lambda,0)$。更新 $w_{12}=W_{11}\hat\beta$。所有列收敛后，由分块求逆公式恢复精度矩阵：$\hat K_{jj}=1/(W_{jj}-w_{12}^\top\hat\beta)$，$\hat K_{-j,j}=-\hat\beta\,\hat K_{jj}$。收敛解满足 KKT 条件：$\hat K_{ij}\ne0$ 时 $W_{ij}-S_{ij}=\lambda\operatorname{sign}(\hat K_{ij})$，$\hat K_{ij}=0$ 时 $|W_{ij}-S_{ij}|\le\lambda$——可以用来检验实现是否正确。

EBIC（Foygel & Drton, 2010）在 BIC 上加了对“候选模型数量”的惩罚：

$$\text{EBIC}_\gamma(\lambda)=-2\,\ell(\hat K_\lambda)+E_\lambda\log n+4\gamma\,E_\lambda\log p,$$

$E_\lambda$ 是 $\hat K_\lambda$ 中非零边数。$\gamma=0$ 时退化为普通 BIC；$\gamma$ 越大越偏向稀疏网络。$\gamma=0.5$ 是常用默认，研究目标偏向“少报假边”时取大，偏向探索（如样本量较小但希望保留弱边）时可取 0.25 甚至 0，但必须事先说明。

算法 1  EBICglasso
输入：相关矩阵 S（p×p，可为多分格相关）；样本量 n；γ；路径长度 L（如 100）；比例 r（如 0.01）
输出：选定的偏相关网络 W
1. λ_max ← max_{i≠j} |S_ij|；λ_min ← r·λ_max
2. 在 [λ_min, λ_max] 上取 L 个对数等距的 λ
3. 对每个 λ：用 glasso 解出 K_λ（可用上一个 λ 的解作热启动）
4.     E_λ ← K_λ 非对角非零元素对数
5.     EBIC(λ) ← −n[log det K_λ − tr(S K_λ)] + E_λ log n + 4γ E_λ log p
6. λ* ← argmin EBIC(λ)
7. W_ij ← −K*_ij / sqrt(K*_ii K*_jj)，i ≠ j
复杂度：每个 λ 的 glasso 约 O(p³)（每轮块坐标下降），总计约 O(L·p³)；p≈40 时几乎瞬时完成

C. Ising 模型的伪似然估计

Ising 模型的对数似然含 $\log Z(\tau,\omega)$，$Z$ 是 $2^p$ 项之和。伪似然用全部条件分布的乘积代替联合分布：

$$\text{PL}(\tau,\omega)=\sum_{k=1}^n\sum_{i=1}^p\log P\big(x_{ki}\mid \mathbf x_{k,-i}\big),\qquad P(x_i=1\mid\mathbf x_{-i})=\sigma\Big(\tau_i+\sum_{j\ne i}\omega_{ij}x_j\Big).$$

由于 $\omega_{ij}$ 同时出现在节点 $i$ 和 $j$ 的条件分布中，最大化伪似然可以分解为 $p$ 个独立的逻辑回归（各自估计 $\beta_{ij}$），再用 AND/OR 规则把 $\beta_{ij}$ 与 $\beta_{ji}$ 合并。每个回归加 $\ell_1$ 惩罚，以节点级 EBIC 选 $\lambda$：

$$\text{EBIC}_i(\lambda)=-2\,\ell_i(\hat\beta_\lambda)+k_\lambda\log n+2\gamma\,k_\lambda\log(p-1),$$

$k_\lambda$ 为节点 $i$ 回归中的非零系数个数。Gibbs 抽样可以从给定参数的 Ising 模型生成数据：依次对每个节点按上面的条件概率重新抽取 0/1，多轮后收敛到联合分布。

D. 强度、预期影响与桥中心性

由 $W$ 定义：强度 $s_i=\sum_j|W_{ij}|$；一步预期影响 $\text{EI}^{(1)}_i=\sum_jW_{ij}$；Robinaugh 等（2016）还定义了两步预期影响 $\text{EI}^{(2)}_i=\text{EI}^{(1)}_i+\sum_jW_{ij}\,\text{EI}^{(1)}_j$，把邻居的影响也计入。若网络全为正边，$s_i=\text{EI}^{(1)}_i$；负边越多，两者差别越大。给定社团划分，把求和限制在其他社团的节点上即得桥强度与桥 EI（见桥中心性与共病）。偏相关边与因果作用的区别（包括对撞偏倚）见横断面统计依赖与因果作用。

E. 病例删除 bootstrap 与 CS 系数

算法 2  中心性稳定性（case-dropping）与 CS 系数
输入：数据 X（n×p）；中心性函数 c(·)（如强度）；删除比例集合 P = {0.05, 0.1, …, 0.75}；
      每个比例的重抽样次数 B（如 1000）；相关阈值 r0 = 0.7；置信水平 q = 0.95
输出：CS 系数
1. c0 ← c(EBICglasso(X))                       // 全样本中心性向量，长度 p
2. 对每个 π ∈ P：
3.     对 b = 1..B：
4.         无放回抽取 ⌊(1−π)n⌋ 名被试，得子样本 X_b
5.         c_b ← c(EBICglasso(X_b))
6.         ρ_b ← cor(c_b, c0)
7.     prop(π) ← #{ρ_b ≥ r0} / B
8. CS ← max{ π ∈ P : prop(π) ≥ q }（若无满足者则 CS = 0）
9. 报告：CS ≥ 0.25 为最低可接受，≥ 0.5 较好
复杂度：|P|·B 次 glasso 估计，每次约 O(L·p³)；可完全并行

注意三点：（1）这里用无放回子抽样，不是常规 bootstrap；（2）子样本中每次都要重新选 $\lambda$，否则低估不稳定性；（3）CS 只衡量“排序稳不稳”，不衡量“排序对不对”，有偏样本中排序也可以很稳。边权准确性的非参数 bootstrap 则是有放回抽取 $n$ 人，对每条边取 2.5% 和 97.5% 分位数；由于 LASSO 会把边压向 0，这些区间并不是严格意义上的置信区间，只宜用来判断边的大小是否可区分。

F. 网络比较检验的置换算法

算法 3  NCT（独立两组）
输入：X1（n1×p），X2（n2×p）；网络估计函数 est(·)；置换次数 R（如 1000–5000）
输出：M、S 与各边的置换 p 值
1. W1 ← est(X1)，W2 ← est(X2)
2. M0 ← max_{i<j} |W1_ij − W2_ij|；S0 ← | Σ_{i<j}|W1_ij| − Σ_{i<j}|W2_ij| |；D0_ij ← |W1_ij − W2_ij|
3. 合并 X ← [X1; X2]
4. 对 r = 1..R：
5.     随机打乱 X 的行，前 n1 行为 X1*，其余为 X2*
6.     W1* ← est(X1*)，W2* ← est(X2*)；计算 M_r、S_r、D_r
7. p_M ← (1 + #{M_r ≥ M0}) / (R + 1)；p_S 同理；p_ij ← (1 + #{D_r,ij ≥ D0_ij}) / (R + 1)
8. 对 {p_ij} 做多重比较校正（如 Holm）
复杂度：2R 次网络估计

置换检验的有效性依赖“零假设下组标签可交换”：两组除网络结构外，变量的边际分布（均值以外的方差等）若明显不同，也可能导致显著结果，解读时需要同时检查。

G. GVAR：三种网络的估计

对第 $k$ 个人、第 $t$ 次测量的 $\mathbf y_{kt}$：

$$\mathbf y_{kt}=\boldsymbol\mu_k+B(\mathbf y_{k,t-1}-\boldsymbol\mu_k)+\boldsymbol\varepsilon_{kt},\qquad \boldsymbol\varepsilon_{kt}\sim\mathcal N(\mathbf 0,\Sigma_\varepsilon),\qquad \boldsymbol\mu_k\sim\mathcal N(\boldsymbol\mu,\Sigma_\mu).$$

• 时间网络：$B$（有向；也可标准化为偏有向相关）；
• 同期网络：$\Theta=\Sigma_\varepsilon^{-1}$ 标准化取负，$\rho^{\text{contemp}}_{ij}=-\Theta_{ij}/\sqrt{\Theta_{ii}\Theta_{jj}}$；
• 个体间网络：$\Sigma_\mu^{-1}$ 标准化取负。

mlVAR 的两步近似：（1）用每人的样本均值 $\bar{\mathbf y}_k$ 代替 $\boldsymbol\mu_k$，对滞后预测变量做个体中心化，逐变量拟合多层回归，得到 $B$ 的固定效应与随机效应；（2）由残差估计 $\Sigma_\varepsilon$（同期网络），由各人均值之间的回归估计 $\Sigma_\mu^{-1}$ 的结构（个体间网络）。直接把全部观测合并计算协方差，得到的是 $\Sigma_\mu$ 与个体内协方差的混合，这就是代码示例中合并相关接近 0 的原因。