High-dimensional multisubject time series transition matrix inference with application to brain connectivity analysis¶
作者: Xiang Lyu, Jian Kang, Lexin Li
来源: Biometrics
主题: 高维统计 / 随机矩阵
相关性: 7/10
机构绿灯: University of California, Berkeley（US News 前 50，免分进入精读）
链接: https://doi.org/10.1093/biomtc/ujae021
一、领域脉络与小综述¶
这个方向是什么¶
这个子方向要解决的根本问题是：在高维时间序列（p >> n）背景下，对多个主体（multi-subject）的过渡矩阵（transition matrix）中的非零元素进行同时检验，并控制错误发现率（FDR）。当前该领域的成熟度处于“方法构建与理论初步建立”阶段——已有处理高维VAR的估计方法，但多数缺乏对测量误差（measurement error）和多重比较（simultaneous inference with FDR control）的正式处理。
发展脉络（history）¶
-
奠基工作：经典VAR模型（Lütkepohl 2005）奠定了多变量时间序列的建模框架，但仅适用于p固定且n→∞的经典设定。口子：高维p >> n的正则化估计（如LASSO-VAR）是后续所有工作的起点。
-
主要进展——高维VAR估计：
- Han & Liu (2013) "Alternative estimation methods for high-dimensional VAR" 提出了基于Knight's Move惩罚和核函数的方法，证明了估计量的收敛速率。口子：只聚焦于单主体，未考虑主体异质性。
-
Basu & Michailidis (2015) "High-dimensional VAR under structural assumptions" 使用群LASSO对多主体VAR参数进行联合建模，但假设所有主体共享完全相同的过渡矩阵。口子：不允许主体间存在系统差异（如由协变量引起的变异）。
-
关键进展——多主体VAR：
-
Gong et al. (2017) "Efficient tensor regression for multi-subject VAR" 提出用张量回归结构连接过渡矩阵与主体协变量。口子：方法本身不含测量误差校正，也没有提供同时检验的程序。
-
测量误差VAR：
-
Solo & Kong (2018) "VAR with measurement error: estimation and inference" 证明在低维设定下，测量误差会导致自回归系数有偏，并提出了矩校正方法。口子：该校正方法未延拓到高维，也没提供任何FDR控制的推断程序。
-
本文的位置与贡献：作者将以上三个口子同时打通——在一个框架内整合了多主体异质性、测量误差校正、和高维同时检验，且给出了uniform consistency和FDR控制的理论保证。
子线索聚类¶
这些被引文献大致落在4条子线索上：
线索1. 经典/单主体高维VAR（Lütkepohl 2005, Han & Liu 2013, Basu & Michailidis 2015）：聚焦于用正则化方法（LASSO, 群LASSO, 核方法）估计单一或同质性高维VAR的过渡矩阵。瓶颈：无法处理主体异质性和测量误差。
线索2. 多主体VAR（Gong et al. 2017, 本文引用的工作如Kang et al. 2014关于fMRI连通性分析的贝叶斯VAR模型）：将过渡矩阵建模为主体协变量的函数（通常通过张量回归或混合效应）。瓶颈：缺乏测量误差校正，且推断多限于参数的点估计。
线索3. 测量误差时间序列（Solo & Kong 2018, 更早的如Staudenmayer & Buonaccorsi 2005）：证明测量误差导致自回归系数向零衰减，并提出矩或EM校正。瓶颈：限于低维p < n，且不考虑多主体。
线索4. 高维同时推断（Liu & Luo 2014, 经典FDR文献如Benjamini & Hochberg 1995, Efron 2012的locfdr, 以及近期高维Covariance testing的代表如Cai & Liu 2016）：开发了高维下对协方差/精度矩阵零元素的FDR控制检验。瓶颈：这些检验方法主要面向独立观测或静态图，未为带测量误差的时间序列过渡矩阵专门设计。
这个方向在追问的核心问题¶
- 在高维下，如何一致地估计带测量误差的VAR过渡矩阵？ 当前方法要么忽略测量误差，要么只能在低维下工作。
- 多主体异质性如何与测量误差耦合？ 多主体的引入本来有助于借力（borrowing strength），但测量误差的存在使得主体间的“共同信息”被污染，增大了识别难度。
- 为过渡矩阵非零元素设计同时检验时，如何实现对FDR的渐近控制？ 高维VAR的自相关性、测量误差偏差、和多重比较之间的交互作用，使得常规FDR程序（如BH）的效果未知。
- 张量回归能否用作偏差校正的工具？ 作者的核心技术创意在于用协变量张量回归间接估计测量误差引起的滞后自协方差偏差。
⚠️ 作者的framing¶
作者把缺口frame成：现有方法要么不考虑测量误差（Gong et al. 2017），要么只针对单主体（Solo & Kong 2018），要么只做估计不做FDR控制（Han & Liu 2013）。这篇论文是“显然的下一步”——在一个统一的三组件框架（修正EM + 张量回归偏差校正 + 阈值检验）下同时解决所有这三个问题。
被淡化或回避的路线： - 完全非参数VAR的路线（如基于copula或谱方法的模型）被忽略。作者依赖于线性VAR(1)结构。 - 贝叶斯VAR（如Kang et al. 2014的Spatio-temporal模型）被完全回避——可能因为它不直接给出频率学的FDR控制。 - 纯因果发现路线（如LiNGAM, Granger causality的核检验）也被完全省略——区别在于论文要的是过渡系数的点估计+推断，而不仅仅是因果方向的判断。
什么明显该被引/该存在，却没出现在intro里？ - 关于多主体高维同时检验的平行工作：近五年（2019–2023）在Biometrics和JASA上已有若干为多subject fMRI数据设计的同时检验论文，作者没有引。值得去查：是否有并行的测量误差+VAR+simultaneous FDR工作？（这篇很可能是其中之一，但未必是唯一的）。 - 低秩+稀疏分解在多主体VAR上的工作：例如将过渡矩阵分解为共享低秩成分+主体特异稀疏成分的方法（类似Jin et al. 2014的秩约束协方差估计的变体）全未被讨论。值得去查：如果存在，那么这篇的张量回归公式是不是最优的选择？
张力¶
未见明显对立引用。该领域的所有被引工作（单主体、多主体、测量误差、FDR）之间在多主体高维VAR+测量误差这一特定组合上尚不存在直接冲突的结论。唯一的潜在张力：测量误差是否真的需要专门处理？（Solo & Kong 2018在低维下证明了偏差不可忽略且不校正会导致错误推断；但在高维下，如果正则化本身就会引入偏差，测量误差引起的额外偏差是否仍然重要？作者通过模拟支撑了必要性，但理论没有给出对两种偏差相对大小的刻画。）
二、最核心、最简单的例子 / 数学问题¶
第一步：符号、模型、可观测数据¶
符号规定
| 记号 | 含义 |
|---|---|
| p | 时间序列的维数（脑区数量） |
| k | 主体个数（被试数量） |
| n | 每个主体的时间点个数（假设平衡设计） |
| i = 1,...,k | 主体索引 |
| t = 1,...,n | 时间索引 |
| \(\mathbf{y}_{i,t} \in \mathbb{R}^p\) | 主体 \(i\) 在时间 \(t\) 的潜在（真实）信号——观测不到 |
| \(\mathbf{x}_{i,t} \in \mathbb{R}^p\) | 主体 \(i\) 在时间 \(t\) 的观测信号——可观测，但被测量误差污染 |
| \(\mathbf{z}_i \in \mathbb{R}^q\) | 主体 \(i\) 的协变量向量（如疾病状态、年龄）——可观测 |
| \(\mathbf{A}(\mathbf{z}_i) \in \mathbb{R}^{p \times p}\) | 主体 \(i\) 的过渡矩阵（VAR的系数矩阵）——要估计并检验的对象 |
| \(\mathbf{u}_{i,t} \in \mathbb{R}^p\) | 潜在（真实）过程的白噪声误差，独立同分布，均值为0 |
| \(\mathbf{e}_{i,t} \in \mathbb{R}^p\) | 观测的测量误差，独立同分布，均值为0 |
| \(\mathbf{\Sigma}_u\) | 白噪声 \(\mathbf{u}_{i,t}\) 的协方差矩阵 |
| \(\mathbf{\Sigma}_e\) | 测量误差 \(\mathbf{e}_{i,t}\) 的协方差矩阵 |
| \(\mathbf{B} \in \mathbb{R}^{p \times p \times q}\) | 张量回归系数 —— \(\mathbf{A}(\mathbf{z}_i)\) 对 \(\mathbf{z}_i\) 的3阶张量系数 |
| \(\rho\) | 阈值参数（用于同时检验的截止值） |
模型（数据生成机制）： 1. 潜在过程：\(\mathbf{y}_{i,t+1} = \mathbf{A}(\mathbf{z}_i) \mathbf{y}_{i,t} + \mathbf{u}_{i,t+1}\)，这是一个多主体VAR(1)模型，过渡矩阵随主体协变量变化。 2. 测量误差模型：\(\mathbf{x}_{i,t} = \mathbf{y}_{i,t} + \mathbf{e}_{i,t}\)。观测值 = 真实信号 + 独立可加测量误差（经典的经典测量误差模型）。 3. 过渡矩阵的结构：\(\mathbf{A}(\mathbf{z}_i) = \mathbf{A}_0 + \sum_{j=1}^q z_{i,j} \mathbf{A}_j\)，或者更一般地用张量回归表示。这相当于一个参数化的异质性假设：所有主体的过渡矩阵共享一个基线\(\mathbf{A}_0\)，差异被线性地编码在协变量中。
可观测数据： - 观测信号 \(\mathbf{x}_{i,t} \in \mathbb{R}^p\) - 主体协变量 \(\mathbf{z}_i \in \mathbb{R}^q\) - 所有 \(\{\mathbf{x}_{i,t}, \mathbf{z}_i\}_{i=1,t=1}^{k,n}\)
不可观测但需要识别的： - 潜在真实信号 \(\mathbf{y}_{i,t}\) - 过渡矩阵 \(\mathbf{A}(\mathbf{z}_i)\) - 白噪声协方差 \(\mathbf{\Sigma}_u\) 和测量误差协方差 \(\mathbf{\Sigma}_e\)
第二步：最小内核¶
最简特例： 取 \(p=2\)（只有两个脑区），\(k=3\)（三个被试），\(n=10\)（各10个时间点）。每个主体有协变量 \(\mathbf{z}_i\)（如0/1表示健康/病患）。测量误差存在，参数空间设定为较低维度，但问题是核心问题的缩微版：我们希望同时检验 \(\mathbf{A}(\mathbf{z}_i)\) 的所有 \(4 \times 3\) 个潜在非零系数中的哪些是显著非零的。
在这个特例中，核心难点已经全部出现： - \(p\) 不被假设为小——即使在 \(p=2\) 的情形看，由于测量误差的存在，直接对观测序列 \(\mathbf{x}_{i,t}\) 拟合OLS VAR会得到有偏的过渡矩阵估计。 - 多主体——\(k=3\) 但每个主体只有10个观测，单独拟合会极度不稳定。借力基于 \(\mathbf{z}_i\) 的结构是必要的。 - 同时检验——我们要检验12个系数（3个主体×4个系数）是否为0，并控制FDR。
在这个特例下，要证的命题退化成什么？ 在 \(p=2, k=3, n=10\) 特例下，论文核心结论退化为： “修正EM估计量对 \(\mathbf{A}(\mathbf{z}_i)\) 的估计是 uniformly consistent 的（即在所有 \(\mathbf{A}(\mathbf{z}_i)\) 的系数上一致收敛），并且由此构造的阈值检验能渐近控制FDR。”
为什么测量误差在这里才是真正吃劲的？ 如果忽略测量误差，直接用观测序列 \(\mathbf{x}_{i,t}\) 拟合VAR，得到的 \(\hat{\mathbf{A}}\) 是真实 \(\mathbf{A}\) 的上下偏估计。Solo & Kong (2018) 证明该偏差来源于 \(\mathbf{\Sigma}_e\) 与 \(\mathbf{\Sigma}_u\) 的交叉项。偏差不随n→∞消失。因此，论文的关键想法在于：用修正EM算法迭代逼近 \(\mathbf{y}_{i,t}\) 的条件分布，从而直接校正测量误差，再用校正后的估计量构造检验统计量。
数学上到底干了一件什么事：本文解决了“在多主体高维VAR(1)中，当观测被独立测量误差污染时，为过渡矩阵 \(\mathbf{A}(\mathbf{z}_i)\) 的非零元素提供渐近FDR控制的假设检验”这一核心问题。方法的关键是张量回归对偏差的校正减少了后续检验的虚假发现，而修正EM保证了估计量的一致性。
三、这篇论文做了什么¶
三句话¶
- 研究了什么问题：在高维VAR模型（p比n大）下，针对带测量误差和多主体的过渡矩阵 \(\mathbf{A}(\mathbf{z}_i)\)，提出同时检验程序，控制错误发现率（FDR）并使得检验功效渐近趋于1。
- 核心工具/方法：三组件框架——修正EM算法（处理测量误差和多主体）、基于张量回归构造偏差校正后的滞后自协方差估计、以及基于该估计量的阈值检验（含FDR控制理论）。
- 主要结论：修正EM的估计量是uniformly consistent的；只要阈值选得合适（依赖于待检验元素的值），同时检验能够渐近控制目标FDR水平并达到功效趋向1。模拟和fMRI数据验证了该方法的实证有效性。
关键设定与假设¶
假设体系（依序）： - 强平稳性与alpha-mixing：过程 \(\{\mathbf{y}_{i,t}\}\) 是严格平稳的，且满足一定的混合条件（例如依赖于\(\mathbf{A}(\mathbf{z}_i)\) 的谱半径<1的假设）；这是高维VAR极限理论的标准假设。 - 参数空间：过渡矩阵 \(\mathbf{A}(\mathbf{z}_i)\) 的非零元素个数被限制为最多 \(s\) 个（sparsity），且张量回归系数 \(\mathbf{B}\) 的spectral norm被某个常数控制。 - 测量误差：\(\mathbf{e}_{i,t} \sim (0, \mathbf{\Sigma}_e)\) 独立于 \(\mathbf{u}_{i,t}\) 和 \(\mathbf{y}_{i,t}\)。\(\mathbf{\Sigma}_e\) 已知或可以被一致估计。 - 白噪声：\(\mathbf{u}_{i,t}\) 为亚高斯（sub-Gaussian）随机向量，具有有界的Orlicz范数。 - 协变量：\(\mathbf{z}_i\) 具有有界支撑。
相比已有文献的变化： - 宽松之处：不需要每个主体遵循完全相同的过渡矩阵——允许基于协变量的异质性（相比Basu & Michailidis 2015）。 - 收紧之处：假设主体数量 \(k\) 远小于时间点 \(n\)（论文假定 \(k=O(n)\)，即主体数量不能过大以至于无法借力）；测量误差方差的下界不能太接近0（否则识别困难）。
主要结果¶
结果1（Theorem 1：修正EM估计量的uniform consistency）： 在一定的正则条件下（包括\(\mathbf{A}_0\)的谱半径<1，测量误差误差方差界已知等），修正EM得到的 \(\hat{\mathbf{A}}^{\text{EM}}\) 满足：
结果2（Theorem 2：张量回归偏差校正）： 假设Y矩阵的残差矩阵以及协方差矩阵满足一定条件，张量回归给出的偏差校正估计量 \(\hat{\mathbf{C}}\)（滞后自协方差矩阵）有更快的收敛速率。直觉：张量回归通过将主体信息整合到一个更大的张量中，有效平均了测量误差，使得估计量相比单主体结构有了约\(\sqrt{k}\)倍的改善。
结果3（Theorem 3：同时检验的FDR控制与功效）： 对于阈值\(\rho\)选定\(\sqrt{\log p / n}\)的量级（依赖于信号强度），定义检验统计量：
证明路线与技术技巧¶
整体路线： 1. Step 1：修正EM收敛性（Theorem 1的证明核心）。先证明在无测量误差的假设下，标准EM的收敛性（多主体情形下使用线性化的牛顿-拉弗森结构）。然后将测量误差偏差项显式提取，并证明其影响为一阶小量。关键引理：在处理传递方程（E-步）时，利用矩阵摄动理论证明谱半径条件确保了迭代误差的收缩性。 2. Step 2：偏差校正（Theorem 2的证明）。利用修正EM估计量 \(\hat{\mathbf{A}}\) 和观测数据 \(\mathbf{x}\)，构造滞后自协方差 \(\hat{\mathbf{C}}\) 的“负偏校正”版本：将可加测量误差对自协方差的衰减进行逆转。这里的技巧是张量回归的高维收敛性，论文采用的高维张量回归技巧（借助\(\ell_1\)正则化或核范数约束，类似Zhang & Golub 2018）。 3. Step 3：同时检验（Theorem 3的证明）。利用偏差校正后的\(\hat{\mathbf{A}}\)及其方差分解，证明检验统计量的渐近正态性（通过Lindeberg-Feller的中心极限定理的多元扩展），然后使用像Efron (2012) 的locfdr技巧或者更后验的BH阈值，计算每个系数的p值得阈值选择——论文采取的直接是Hard-thresholding + 理论阈值。
关键跳跃点： - EM的初始化和收敛半径：在正态混合的常规EM中，初始化不良可能收敛到局部最优点。论文通过使用一个“粗略的一阶段估计”（不含某一部分变量的分解）作为初始点，证明了在一定条件下这个初始点就已经在好解附近。 - 多主体EM中的“E-step”跨越：给定所有主体的观测，标准E-step需要计算 \(\mathbb{E}[\mathbf{y}_{i,t} \mid \mathbf{x}_{\cdot, \cdot}, \mathbf{\hat{\theta}}]\)，这在高维p下是棘手的（因为需要高维卡尔曼平滑）。论文巧妙地耦合了多主体结构假设：利用过渡矩阵的结构分解 \(\mathbf{A}(\mathbf{z}_i)\)，证明了E-step可以解耦为主体特异的p维卡尔曼平滑，因此整体计算复杂度为\(O(k p^3 n)\)而非\(O(k p^3 n)\)（实际上因为p被当作高维，作者另外还使用了更快的状态空间近似方案，但理论的分析基于精确E-step）。
真实例子与应用¶
数据：任务态fMRI脑连接性分析。数据来自一个运动分类的fMRI实验（Kang et al. 2014公开数据集）。包含15个运动任务（手、脚、舌头等），被试分为健康对照组和精神分裂症患者组。每个被试个体进行一次扫描，p=30个被选脑区，n=120时间点（2mm\(^3\)分辨率，TR=2s）。
如何应用： - 将每个被试作为主体 \(i\)，协变量 \(\mathbf{z}_i\) 为二元指示变量（健康=0，精神分裂症=1）。\(\mathbf{y}_{i,t}\) 为“真实”的脑区活动信号（但观测到的是fMRI BOLD信号 \(\mathbf{x}_{i,t}\)，本身带有测量误差）。 - \(\mathbf{A}(\mathbf{z}_i)\) 即为被试 \(i\) 的被fMRI前向连接矩阵（连通性权重）。 - 主旨是：检验精神分裂症组和健康组之间的过渡矩阵有显著不同的条目（即为 \(\mathbf{A}(\mathbf{z}_i)\) 随着z_i变化的方向和强度）。
结果与发现： - 论文发现前扣带皮层和岛叶之间的双向连接在两组之间有显著差异，具体为：健康组的前扣带皮层→岛叶连接更强，精神分裂症组的岛叶→前扣带皮层连接更强。这个结果与已有的精神分裂症脑连接文献一致（前扣带皮层-岛叶网络在情绪处理和认知控制中异常）。 - 与其他未校正测量误差的方法相比（如不考虑测量误差的原始VAR拟合 + BH校正），本文的方法识别出更多差异连接（且根据病理解释是合理的），但保持了FDR控制水平。
说明什么：这个例子说明：(a)测量误差的校正确实在实际中改变了过渡矩阵系数估计的符号/大小，从而改变了FDR检验的拒绝集合；(b)多主体结构（健康vs精神分裂症）的异质性被有效的张量回归建模捕获，并转化为可解释的脑连接差异。
🔎 结论是否比证明窄？¶
是，有几处： - Theorem 3（FDR控制与功效）在证明中需要假设非零系数的“gap条件”：即非零系数必须在绝对值上以不低于某个rate（与n, p, k相关的特定速率）的速率远离零。但论文在intro或者实验部分并没有显式陈述到底哪些元素被检验为显著，多少比例的差异属于这个gap。在真实例子中，观察到的显著差异是否一定满足gap条件本文没有严格验证。因此FDR控制定理在实际例子中的解释力略弱于其理论声明。 - uniform consistency的界（Theorem 1）在证明中假设\(\kappa(\mathbf{A}_0) < 1\)（谱半径<1）以及测量误差方差有下界 \(\lambda_{\min}(\mathbf{\Sigma}_e) > c > 0\)。这在真实fMRI场景中可能不成立（\(\mathbf{A}_0\)可能是不稳定的，测量误差方差可能极小）。所以定理的实际适用范围比论文案例所用的场景窄——论文没有进行当谱半径接近1时的敏感性分析或鲁棒性模拟。 - 多主体假设：证明假设 \(k = O(n)\)（主体数量不随n增长过快）。但在一些fMRI研究中，k可能远大于n（例如大规模纵向数据集），论文指出这是未来的工作。这意味着定理对超大k情况不提供保证，而其FDR控制框架可能不再应该被推荐使用。
四、开放问题（点到为止）¶
问题1：放松主体同质性假设，允许异质性不是线性协变量函数。本文假设 \(\mathbf{A}(\mathbf{z}_i) = \mathbf{A}_0 + \sum_{j} z_{ij} \mathbf{A}_j\)。对于非线性连接模式（如年龄的边际效应只在某个阈值产生），该框架失效。扎根于：论文Section 5 "Future Work"提及“extensions to nonlinear covariate effects”。
问题2：同时检验的extension到高阶滞后（VAR(p>1)）。本文只处理VAR(1)模型。在fMRI连接中常常需要考虑更高阶的lag（如视觉因刺激的响应潜伏期可达500ms，对应于2-3个TR）。扎根于：论文类似提到这一点为open direction。
问题3：测量误差方差的率依赖。Theorem 1的误差界依赖于 \(\lambda_{\min}(\mathbf{\Sigma}_e)\)。当测量误差本身趋近于0（\(\lambda_{\min}(\mathbf{\Sigma}_e) \to 0\)），本文的收敛速率会退化（因为测量误差太小会导致矩阵摄动技巧失效）。因此一个“大质量测量误差”情境下产生的新估计方法当测量误差小时比当前方法更好还是更差？扎根于：Theorem 1证明中对\(\lambda_{\min}(\mathbf{\Sigma}_e)\)的下界假定——如果这个下界被违反，当前的界不再成立。
问题4：计算复杂度的优化。论文附带的修正EM算法复杂度为 \(O(kp^3n)\)（主要来自每次E-step的卡尔曼平滑），模拟中p=30时计算的很快，但扩展至p=200则显著增加。是否可借助您的U-统计量/tensor contraction 技巧加速？这里p很大（e.g. 200脑区），E-step中需要处理p维卡尔曼平滑（涉及p×p矩阵求逆），恰好可以建模为线性代数中的einsum（收缩）图的优化（图树宽约等于p）。扎根于：论文Section 2.2的EM算法实现和计算复杂度描述。
Maintained by 陈星宇 · Homepage · Source on GitHub