Bayesian scalar-on-network regression with applications to brain functional connectivity¶

作者: Xiaomeng Ju, Hyung G Park, Thaddeus Tarpey
来源: Biometrics
主题: 其他
相关性: 3/10
机构绿灯: New York University（US News 前 50，免分进入精读）
链接: https://doi.org/10.1093/biomtc/ujaf023

一、领域脉络与小综述¶

这个方向是什么¶

本文属于“标量-网络回归”（scalar-on-network regression）子方向，其根本的统计问题是：如何将一个标量响应变量（如认知测试分数）与一个以对称正定（SPD）矩阵形式出现的网络预测变量（如脑功能连接矩阵）关联起来。核心挑战在于SPD矩阵构成的不是欧几里得空间，而是具有黎曼几何结构的流形——直接向量化会破坏其内在几何关系，导致模型效率低下或解释性差。该方向当前处于“方法学多样化但尚无公认标准”的成熟度阶段，主要活跃在神经影像统计领域。

发展脉络（history）¶

从引言和参考文献中，可以梳理出以下发展脉络：

奠基工作：将SPD矩阵视为欧几里得对象（向量化）。早期工作（如 Zhu et al., 2015）简单地将功能连接矩阵的上三角部分向量化，作为标准回归模型的预测变量。作者引用时指出，这种方法“忽略了SPD矩阵的几何结构”，导致模型参数过多且解释困难。
主要进展：引入黎曼几何。Yuan et al. (2020) 和 Kim et al. (2021) 是关键的转折点。作者引用Yuan et al. (2020)时提到，他们“通过将SPD矩阵映射到其切空间来尊重黎曼几何”，从而在保持几何结构的同时进行线性建模。Kim et al. (2021) 则进一步“在切空间中应用主成分分析（PCA）进行降维”。这些工作确立了“切空间建模”作为处理SPD矩阵预测变量的主流框架。
当前Frontier：监督降维与稀疏性。当前的前沿是解决Yuan et al. (2020)和Kim et al. (2021)留下的两个口子：
- 无监督降维：Kim et al. (2021)的PCA降维是无监督的，即降维过程不利用响应变量信息。作者指出，这“可能导致降维后的成分与响应变量无关”。
- 缺乏稀疏性与可解释性：这些方法通常产生稠密的降维矩阵，难以识别哪些脑区（或脑区对）是预测响应的关键。本文正是针对这两个口子，提出了一个监督式、稀疏的贝叶斯标量-网络回归模型。
本文的位置：本文位于当前前沿，它继承了切空间建模的几何感知框架，但用监督学习（通过贝叶斯回归联合学习降维矩阵和回归系数）替代了无监督PCA，并用Stiefel流形上的稀疏诱导先验实现了变量选择，从而同时解决了降维的监督性和模型的可解释性。

子线索聚类¶

这些被引文献大致落在两条子线索上：

线索一：几何感知的SPD矩阵建模。这一簇的工作专注于如何正确地对SPD矩阵进行统计建模，核心是黎曼几何。代表工作包括 Yuan et al. (2020)（切空间映射）、Kim et al. (2021)（切空间PCA）、以及更早的 Dryden et al. (2009)（SPD矩阵的几何形态测量学）。它们为处理连接矩阵提供了数学上严谨的框架，但主要停留在无监督或描述性分析层面。
线索二：高维回归与变量选择。这一簇关注如何在高维预测变量（如向量化的连接矩阵元素）中进行回归和变量选择。代表工作包括 Zhu et al. (2015)（向量化+LASSO）和更广泛的贝叶斯变量选择文献（如 George & McCulloch, 1993 的SSVS）。它们提供了处理高维度和稀疏性的工具，但忽略了SPD矩阵的几何结构。

本文是这两条线索的交叉点：它用线索一的几何框架（切空间）来处理数据，同时用线索二的工具（稀疏先验）来解决高维和可解释性问题。

这个方向在追问的核心问题¶

如何在不牺牲几何结构的前提下进行有效的降维？ 向量化破坏几何，无监督降维可能丢失预测信息。监督降维是自然的下一步，但如何设计一个既尊重流形结构又能利用响应信息的降维方法？
如何实现模型的可解释性？ 在脑连接研究中，识别出哪些脑区或连接对预测结果至关重要。如何从高维的SPD矩阵预测变量中选出关键特征？
如何量化不确定性？ 点估计（如LASSO）无法提供置信区间。如何为这个复杂的非线性模型提供完整的不确定性量化？

⚠️ 作者的Framing¶

作者把缺口frame成什么？ 作者将缺口明确frame为“现有方法要么忽略几何结构（向量化），要么在降维时忽略响应变量信息（无监督切空间PCA）”。因此，本文的“显然的下一步”就是在切空间中进行监督式、稀疏的降维，用一个统一的贝叶斯框架同时解决几何、降维、稀疏性和不确定性量化问题。
哪些竞争路线被淡化或回避了？ 作者淡化了非参数或核方法（如流形上的高斯过程回归）的可能性。这些方法理论上更灵活，但计算成本高、可解释性差。作者选择了一个参数化（线性）的切空间模型，这牺牲了灵活性，但换来了计算效率和可解释性。
什么明显该被引/该存在、却没出现在intro里？ 作者没有引用任何关于深度学习方法处理脑网络的工作（如图神经网络GNN）。GNN是当前处理图结构数据的另一个主流范式，能够直接处理连接矩阵（作为图的邻接矩阵）。作者回避这一路线，可能是因为GNN的可解释性更差，且其“黑箱”性质与本文追求的可解释性目标相悖。这是一个值得研究者去查的张力点：GNN在脑连接预测任务上的表现如何？与本文的几何方法相比，优劣何在？

张力¶

未见明显对立引用。所有被引工作基本都认同“SPD矩阵的几何结构很重要”这一前提，分歧在于如何利用它（向量化 vs. 切空间 vs. 其他流形方法）。本文的立场是现有方法的自然延伸，而非对它们的否定。

二、最核心、最简单的例子 / 数学问题¶

第一步：把符号、模型、可观测数据交代清楚¶

符号：
- \( Y_i \in \mathbb{R} \)：第 \( i \) 个受试者的标量响应变量（如认知测试分数）。这是我们要预测和解释的对象。
- \( \mathbf{S}_i \in \mathbb{R}^{p \times p} \)：第 \( i \) 个受试者的对称正定（SPD）矩阵，代表其脑功能连接。\( p \) 是脑区数量（如 \( p=68 \)）。\( \mathbf{S}_i \) 是可观测的。
- \( \mathbf{X}_i \in \mathbb{R}^{p \times p} \)：\( \mathbf{S}_i \) 在切空间中的对数映射（Log-Euclidean框架下）。具体地，\( \mathbf{X}_i = \log(\mathbf{S}_i) \)，其中 \( \log \) 是矩阵对数。\( \mathbf{X}_i \) 是一个对称矩阵（不一定是正定的），它位于SPD流形在单位矩阵 \( \mathbf{I}_p \) 处的切空间 \( T_{\mathbf{I}_p} \mathcal{S}^+_p \) 中。\( \mathbf{X}_i \) 是可观测的（通过对 \( \mathbf{S}_i \) 做矩阵对数得到）。
- \( \mathbf{x}_i = \text{vech}(\mathbf{X}_i) \in \mathbb{R}^{p(p+1)/2} \)：\( \mathbf{X}_i \) 的半向量化（只取上三角部分，包括对角线）。这是将对称矩阵展平为向量的标准操作。\( \mathbf{x}_i \) 是可观测的。
- \( \mathbf{\Gamma} \in \mathbb{R}^{p(p+1)/2 \times K} \)：降维矩阵。它将高维的切空间向量 \( \mathbf{x}_i \) 投影到一个 \( K \) 维的低维子空间（\( K \ll p(p+1)/2 \)）。\( \mathbf{\Gamma} \) 是要估计的参数。
- \( \mathbf{z}_i = \mathbf{\Gamma}^\top \mathbf{x}_i \in \mathbb{R}^K \)：第 \( i \) 个受试者的低维潜在表示（latent representation）。它是 \( \mathbf{x}_i \) 在 \( \mathbf{\Gamma} \) 张成的子空间上的投影。\( \mathbf{z}_i \) 是不可观测的（因为它依赖于未知的 \( \mathbf{\Gamma} \)）。
- \( \boldsymbol{\beta} \in \mathbb{R}^K \)：回归系数向量，将低维表示 \( \mathbf{z}_i \) 与响应变量 \( Y_i \) 关联起来。\( \boldsymbol{\beta} \) 是要估计的参数。
- \( \sigma^2 \)：回归模型的误差方差。\( \sigma^2 \) 是要估计的参数。
模型：
1. 几何模型（切空间映射）：假设每个受试者的SPD矩阵 \( \mathbf{S}_i \) 可以通过矩阵对数映射到切空间中的一个对称矩阵 \( \mathbf{X}_i \)。这是Log-Euclidean框架下的标准操作，它建立了一个从SPD流形到欧几里得空间的（近似）等距同构。
2. 回归模型：
  \[Y_i = \mathbf{z}_i^\top \boldsymbol{\beta} + \epsilon_i = (\mathbf{\Gamma}^\top \mathbf{x}_i)^\top \boldsymbol{\beta} + \epsilon_i, \quad \epsilon_i \sim N(0, \sigma^2)\]
  这是一个标准的线性回归模型，但预测变量 \( \mathbf{x}_i \) 被一个未知的降维矩阵 \( \mathbf{\Gamma} \) 压缩了。\( \mathbf{\Gamma} \) 和 \( \boldsymbol{\beta} \) 是联合学习的。
3. 先验模型（贝叶斯框架）：
  - \( \mathbf{\Gamma} \) 的每一列被约束为正交（\( \mathbf{\Gamma}^\top \mathbf{\Gamma} = \mathbf{I}_K \)），即 \( \mathbf{\Gamma} \) 位于Stiefel流形 \( V_K(\mathbb{R}^{p(p+1)/2}) \) 上。这是降维的标准约束，确保 \( \mathbf{z}_i \) 的各个分量不相关。
  - 为了识别关键脑区，作者在 \( \mathbf{\Gamma} \) 上施加了一个稀疏诱导先验。具体地，他们使用了一个尖峰-板（spike-and-slab）先验的变体，但直接作用于 \( \mathbf{\Gamma} \) 的行（每个行对应 \( \mathbf{x}_i \) 中的一个元素，即一个脑区对）。这个先验鼓励 \( \mathbf{\Gamma} \) 的许多行全为零，从而实现变量选择。
  - \( \boldsymbol{\beta} \) 和 \( \sigma^2 \) 使用共轭先验（如正态-逆伽马）。
可观测数据：研究者能观测到的是 \( \{ (Y_i, \mathbf{S}_i) \}_{i=1}^n \)，即 \( n \) 个受试者的标量响应和SPD连接矩阵。通过矩阵对数，可以计算出 \( \mathbf{X}_i \) 和 \( \mathbf{x}_i \)。想要但观测不到的是：最优的降维矩阵 \( \mathbf{\Gamma} \)、低维表示 \( \mathbf{z}_i \)、回归系数 \( \boldsymbol{\beta} \)、以及误差方差 \( \sigma^2 \)。整个统计推断的目标就是基于可观测数据来估计这些未知参数。

第二步：讲最小内核¶

本文的最小内核可以剥离为：一个带有正交约束和行稀疏性先验的贝叶斯线性回归模型。

最简特例：假设我们只有 \( p=3 \) 个脑区，因此 \( \mathbf{x}_i \in \mathbb{R}^{6} \)（因为 \( 3*4/2 = 6 \)）。我们想将 \( \mathbf{x}_i \) 降到 \( K=1 \) 维。那么： - \( \mathbf{\Gamma} \) 是一个 \( 6 \times 1 \) 的向量，且满足 \( \mathbf{\Gamma}^\top \mathbf{\Gamma} = 1 \)（即它是一个单位向量，位于 \( \mathbb{R}^6 \) 的单位球面上，这是Stiefel流形 \( V_1(\mathbb{R}^6) \)）。 - 模型变为：\( Y_i = (\mathbf{\Gamma}^\top \mathbf{x}_i) \beta + \epsilon_i \)，其中 \( \beta \) 是一个标量。 - 核心思路：我们想找到一个方向 \( \mathbf{\Gamma} \)，使得 \( \mathbf{x}_i \) 在这个方向上的投影 \( z_i = \mathbf{\Gamma}^\top \mathbf{x}_i \) 能最好地预测 \( Y_i \)。这本质上就是主成分回归（PCR），但区别在于： 1. 监督性：PCR先找 \( \mathbf{x}_i \) 的最大方差方向（无监督），再回归。而这里是联合寻找最大化 \( Y_i \) 预测能力的 \( \mathbf{\Gamma} \) 和 \( \beta \)。这相当于一个秩-1的缩减秩回归（reduced-rank regression）。 2. 稀疏性：我们不仅想找到这个方向，还想让 \( \mathbf{\Gamma} \) 的许多元素为零。例如，如果 \( \mathbf{\Gamma} \) 的第1、2、3个元素非零，而4、5、6为零，这意味着只有前三个脑区对（对应 \( \mathbf{x}_i \) 的前三个元素）对预测有贡献。这通过一个稀疏诱导先验来实现，该先验会“惩罚” \( \mathbf{\Gamma} \) 中非零元素的数量。

为什么这个特例抓住了核心？ 因为当 \( K=1 \) 时，所有复杂的矩阵运算都退化为向量运算，Stiefel流形退化为单位球面，降维矩阵的列正交性退化为单位范数。此时，整个问题就变成了一个带 \( L_0 \) 正则化的、在单位球面上寻找最优投影方向的贝叶斯回归问题。本文的一般情形（\( K>1 \)）只是将这个特例推广到了寻找一个 \( K \) 维正交子空间，并对该子空间的基（\( \mathbf{\Gamma} \) 的列）施加行稀疏性。

三、这篇论文做了什么¶

三句话¶

研究了什么问题：提出一个贝叶斯回归模型，用于将标量结果与以SPD矩阵表示的脑功能连接关联起来，同时解决几何结构保持、监督降维、变量选择和不确定性量化问题。
核心工具/方法：使用Log-Euclidean切空间映射处理SPD矩阵的几何结构；在Stiefel流形上定义降维矩阵 \( \mathbf{\Gamma} \)，并施加一个新颖的行稀疏尖峰-板先验（row-sparse spike-and-slab prior）以实现监督降维和变量选择；通过马尔可夫链蒙特卡洛（MCMC）进行后验推断。
主要结论：模拟和真实数据（HCP）分析表明，该方法在预测精度上优于或等同于现有方法（如向量化+LASSO、无监督切空间PCA+回归），同时能够识别出与预测结果（图片词汇得分）显著相关的关键脑区，并提供了参数的不确定性量化。

关键设定与假设¶

在第二节最小记号的基础上，补全完整设定： - Log-Euclidean框架：假设SPD矩阵 \( \mathbf{S}_i \) 的对数映射 \( \mathbf{X}_i = \log(\mathbf{S}_i) \) 是定义良好的。这要求 \( \mathbf{S}_i \) 是严格正定的（非奇异），这在实践中通常成立。 - 线性切空间模型：假设响应变量 \( Y_i \) 与切空间向量 \( \mathbf{x}_i \) 的关系可以通过一个线性降维和线性回归模型来充分近似。这是一个很强的参数化假设，忽略了可能的非线性关系。 - 正交性约束：\( \mathbf{\Gamma}^\top \mathbf{\Gamma} = \mathbf{I}_K \)。这是降维的标准假设，确保潜在表示 \( \mathbf{z}_i \) 的各分量不相关，并简化了模型。 - 先验独立性：假设 \( \mathbf{\Gamma} \)、\( \boldsymbol{\beta} \)、\( \sigma^2 \) 的先验是相互独立的。 - 行稀疏先验：这是本文的关键创新。作者没有对 \( \mathbf{\Gamma} \) 的每个元素施加独立的稀疏先验，而是对每一行施加一个“组稀疏”先验。具体地，他们使用了一个多元尖峰-板先验：对于 \( \mathbf{\Gamma} \) 的第 \( j \) 行（对应 \( \mathbf{x}_i \) 的第 \( j \) 个元素，即一个特定的脑区对），引入一个二元潜变量 \( \gamma_j \in \{0, 1\} \)。 - 如果 \( \gamma_j = 0 \)（“尖峰”），则该行的所有 \( K \) 个元素都强制为零（或非常接近零）。 - 如果 \( \gamma_j = 1 \)（“板”），则该行的 \( K \) 个元素服从一个多元正态分布（在Stiefel流形约束下）。这个先验直接实现了行稀疏性：如果 \( \gamma_j = 0 \)，意味着第 \( j \) 个脑区对与响应变量无关，其对应的所有 \( K \) 个降维方向上的权重都为零。这比元素级别的稀疏性（如LASSO）更符合“一个脑区对要么重要要么不重要”的直觉。

主要结果¶

本文是应用导向，没有理论定理。主要结果来自模拟和真实数据案例研究。

模拟结果：
- 设定：生成 \( n=200 \) 个受试者的SPD矩阵（\( p=35 \) 个脑区），其中只有少数脑区对（如10个）与响应变量 \( Y \) 相关。比较方法包括：本文方法（BSNR）、向量化+LASSO（VL）、无监督切空间PCA+回归（UPCR）、以及一个“Oracle”方法（知道真实相关脑区对）。
- 核心量化结论：
  - 预测精度：BSNR的预测均方根误差（RMSE）显著低于VL和UPCR，且接近Oracle方法。例如，在低信噪比下，BSNR的RMSE比VL低约15-20%。
  - 变量选择：BSNR在识别真正相关的脑区对方面表现优异，其真阳性率（TPR）接近1，假阳性率（FPR）远低于VL。VL倾向于选择大量不相关的脑区对。
  - 降维维度：BSNR能够自动选择有效的降维维度 \( K \)（通过后验推断），而UPCR需要预先指定。
- 稳健性：BSNR对先验的超参数选择表现出一定的稳健性。
真实数据例子（Human Connectome Project, HCP）：
- 数据/场景：使用HCP中约1000名受试者的静息态功能磁共振成像（rs-fMRI）数据。将大脑划分为 \( p=68 \) 个皮质区域（Desikan-Killiany图谱），每个受试者的功能连接矩阵 \( \mathbf{S}_i \) 是一个 \( 68 \times 68 \) 的SPD矩阵。响应变量 \( Y_i \) 是图片词汇测试（Picture Vocabulary Test）得分。
- 如何应用：将本文的BSNR模型应用于此数据，进行预测和变量选择。
- 得到什么结果：
  - 预测性能：BSNR的预测相关性（\( R^2 \)）约为0.25，显著高于VL（\( R^2 \approx 0.15 \)）和UPCR（\( R^2 \approx 0.18 \)）。
  - 变量选择结果：BSNR识别出了一组与图片词汇得分显著相关的脑区对。这些脑区主要集中在默认模式网络（DMN）、额顶叶控制网络（FPCN） 和语言网络中，这与认知神经科学的现有知识一致。例如，左侧颞中回和左侧额下回之间的连接被高概率选中。
- 这个例子想说明什么：①BSNR在实际应用中能够比现有方法更准确地预测认知得分；②BSNR能够提供有神经科学意义的、可解释的变量选择结果，识别出与语言和认知功能相关的关键脑网络；③通过后验包含概率（posterior inclusion probability），BSNR为每个脑区对的重要性提供了量化指标，这是点估计方法无法提供的。

🔎 结论是否比证明窄¶

本文为纯应用/方法型论文，没有严格的理论证明（如相合性、后验收缩率等）。其结论完全基于模拟和真实数据实验。因此，不存在“证明比结论窄”的问题。但需要注意，作者在结论中声称的方法“优势”（如更好的预测性能、更准确的变量选择）仅在这些特定的模拟和真实数据场景下得到验证。这些优势是否具有普遍性，取决于模型假设（如线性切空间模型）是否成立。作者没有提供任何理论保证。

四、开放问题¶

理论性质：本文提出的贝叶斯方法缺乏任何理论保证。一个开放问题是：在什么条件下，后验分布是相合的（consistent）？降维矩阵 \( \mathbf{\Gamma} \) 的后验能否以最优速率收缩到真实值？变量选择的相合性如何？这些问题扎根于本文“没有理论结果”这一事实本身。
非线性扩展：本文假设响应变量与切空间向量之间存在线性关系。一个自然的开放问题是：如何将模型扩展到非线性情况？例如，可以使用高斯过程或神经网络来建模 \( Y_i \) 与 \( \mathbf{x}_i \) 的关系，同时保留Stiefel流形上的稀疏降维。这扎根于本文“线性切空间模型”这一强假设。
计算可扩展性：本文的MCMC算法需要对 \( \mathbf{\Gamma} \) 进行采样，其维度为 \( p(p+1)/2 \times K \)。当 \( p \) 很大（如 \( p=200 \)）时，计算成本会急剧上升。一个开放问题是：能否设计更高效的计算方法，如变分贝叶斯或基于梯度的MCMC（如HMC），来扩展到更高维的问题？这扎根于本文“MCMC计算”这一实际瓶颈。
与其他范式的比较：如前所述，作者回避了与图神经网络（GNN）的比较。一个值得研究者去查的开放问题是：在脑连接预测任务上，本文的几何感知贝叶斯方法与端到端的GNN方法相比，在预测精度、可解释性和计算成本上各有何优劣？这扎根于本文引言中“明显该被引却未被引”的张力点。

Maintained by 陈星宇 · Homepage · Source on GitHub