Regularized Joint Estimation of Related Vector Autoregressive Models

A. Skripnikov; G. Michailidis

doi:10.1016/j.csda.2019.05.007

计算统计数据分析。作者手稿；PMC 2020年11月1日提供。

以最终编辑形式发布为：

计算统计数据分析。2019年11月；139: 164–177.

在线发布2019年5月22日。数字对象标识：2016年10月10日/j.csda.2019.05.007

PMCID公司：项目管理委员会7079674

尼姆斯：美国国家卫生研究院1530078

PMID：32189818

相关向量自回归模型的正则联合估计

A.斯科里普尼科夫^a、，^1,^*和G.米查利迪斯^一

作者信息版权和许可证信息 PMC免责声明

关联数据

补充资料: 1
NIHMS1530078-补充-1.pdf（219000）
GUID:772D09D9-D637-4EE7-AD67-E18946F4882C
2
NIHMS1530078-增补-2.R（1.9万）
GUID:A6143D1A-BDB2-4112-932F-A41790C30E44
三。
NIHMS1530078-补充-3.R（1.7万）
GUID:B02E145F-AE00-41FA-8DDA-E2B3516490DB
4
NIHMS1530078-增补-4.R（1300万）
GUID:4B04F6BE-69EE-4B27-97ED-BB8D92BF787E
5
美国国立卫生研究院1530078增补件5.R（27K）
GUID:D1565AA7-3314-4E88-B7E4-57CA1704D883
6
NIHMS1530078-增补-6.R（16000）
GUID:7DB950D2-721D-4652-A16A-B713037E13DE

摘要

在许多应用程序中，可以访问多个相关主题的高维时间序列数据。一个令人鼓舞的应用领域来自神经成像领域，例如从不同组患有特定神经疾病的受试者（病例/对照）获得的大脑fMRI时间序列数据。讨论了多相关向量自回归（VAR）模型的正则化联合估计问题，利用群拉索惩罚和正则拉索惩罚，通过在模型之间借用强度来提高估计的统计效率。开发了一个建模框架，该框架允许相关受试者的群体级和特定受试者效应，使用群体套索惩罚来估计前者。介绍了一种估计过程，其性能在合成数据上进行了说明，并与其他最先进的方法进行了比较。此外，该方法还用于静息状态fMRI数据的分析。特别是，利用ADHD-200全球竞争数据库中的数据，对注意力缺陷多动障碍（ADHD）患者与对照组相比的大脑区域间时间效应进行了群体级描述性分析。

关键词：注意缺陷多动障碍，组套索，正规估计，静止状态fMRI，稳定性选择，向量自回归

1.简介

随着技术的最新进步和可用数据的不断增加（例如点击生成的网络浏览数据、社交网络、图像和视频数据），人们对高维时间序列数据的建模和分析产生了浓厚的兴趣。应用领域包括基因调控网络推断[29]，大脑fMRI数据[37]宏观经济时间序列预测与结构分析[1,24]举几个例子。它们的共同特点是相对于可用时间点分析大量变量关系，从而导致高维统计估计问题。在许多情况下，自回归模型很好地捕捉到了所考虑数据的时间动态，因此使用向量自回归模型（VAR）可以对变量之间的时间相关性进行建模。然而，在存在大量要估计的参数且只有几个时间点的情况下，需要将适当的稀疏性假设纳入VAR建模框架。为了加强稀疏性，通常使用ℓ₁套索惩罚[2,三]，其理论性质已在前文中确立。

如前所述，在许多应用中，除了典型的高维设置之外，还必须对中等数量到大量相关主题的时间序列进行估计。作为一个鼓舞人心的例子，医学研究领域带来了许多实验设置，随着时间的推移对多个受试者进行监测，例如收集一组患者的fMRI时间序列数据。观察患有特定疾病/障碍（如阿尔茨海默病）的患者，预计大脑区域之间的连接呈现出共同的结构模式。然而，有充分的证据表明，尽管有相同的疾病，患者仍然表现出个体模式[44,4]这导致了学科的异质性。

fMRI时间序列数据的标准分析管道通常包括分别估计每个受试者的网络，随后累积估计值以进行进一步的群体级分析[30]. 这种方法已被应用于研究阿尔茨海默病患者的大脑活动[22]、自闭症[31]、帕金森氏病[25]等等。虽然该方法提供了群体层面的分析工具，但它并没有将同一群体中不同受试者之间的相似性假设（例如案例或对照）纳入评估程序。

这项工作的一个关键贡献是开发了一个联合建模框架，该框架将在评估过程中实施相似性假设，同时也能够评估特定主题的网络效应。该方法通过借用相关时间序列的强度来增加公共结构估计的有效样本量。此外，该框架允许检测每个受试者（如果存在）最显著的个人效应。联合估计问题最近在文献中受到了关注，主要集中在多个图形模型的估计上。这些方法利用各种惩罚，鼓励稀疏性和跨模型参数的联合估计；请参阅中使用的分级惩罚[19]或融合套索惩罚[14,27]. 在这项工作中，我们采用了一种群体套索惩罚，因为它能够清楚地识别多个受试者之间的共同结构，同时让影响的大小发生变化。检测到共同结构后，将应用标准套索程序，以获得对特定主题效果的稀疏估计。

虽然引入的联合估计程序可用于其他时间序列设置（具有共同制造业和工业特征的城市或州的经济数据、匹配患者的基因表达数据、类似商店的销售数据），主要的激励应用是静止状态fMRI时间序列数据，用于研究各种认知障碍的自发脑时间动力学。先前引用的关于认知障碍群体级推理的文献侧重于评估功能性大脑连通性（研究大脑区域信号之间的相关性），而不是推断大脑区域之间的前导滞后关系。静止状态数据表示脑内波动的单调性，VAR模型适用于捕捉此类数据中的时间动态。此外，我们假设每个患者的VAR模型是对所考虑的群体（无论是患有疾病的受试者还是健康对照者）的一些常见的潜在VAR模型的扰动，其结构将通过我们基于群体套索惩罚的联合程序来估计。注意，文献中对神经影像数据中VAR模型的使用提出了某些反对意见[12,32]. 然而，我们试图在我们的ADHD研究中解决其中一些问题第4节例如，在fMRI数据应用中使用VAR模型的六个问题之一[32]关注来自多个研究参与者大脑区域的不同信号强度，尽管他们都共享一个共同的抽象处理结构。在我们的框架中，这个问题是通过一个群体套索惩罚直接解决的，该惩罚在实施公共结构的同时，允许公共效应的大小发生变化。值得一提的是，通过正则化惩罚的联合估计方法以前曾用于脑fMRI时间序列数据[5,11]，但仅用于功能连接性估计，而我们将其应用于推断前导滞后关系。

要引入关节建模框架，请考虑第页-可变平稳时间序列 ${X（X）}_{t吨}^{(k个)} = {({x个}_{1, t吨}^{(k个)}, \dots, {x个}_{第页, t吨}^{(k个)})}^{⊤}$ ,t吨= 1, …，T，k= 1, …,K、，对于K（K）相关科目。具有滞后次序的VAR模型D类，或V AR（车辆配置总成）(D类)，由给出

{X（X）}_{t吨}^{(k个)} = {A类}_{1}^{(k个)} {X（X）}_{t吨 - 1}^{(k个)} + \dots + {A类}_{D类}^{(k个)} {X（X）}_{t吨 - D类}^{(k个)} + ϵ_{t吨}^{(k个)}, ϵ_{t吨}^{(k个)} ~ N个 (0, σ_{(k个)}^{2} 我_{第页}), t吨 = D类, \dots, T型, k个 = 1, \dots, K（K）,

(1)

哪里 ${A类}_{d日}^{(k个)}$ 是一个第页×第页捕捉订单时间效应的转移矩阵d日在第页主题变量k、 d日= 1,…,D类,k个= 1, …,K（K）。我们进一步假设对角误差协方差矩阵 $Σ_{k个} = σ_{(k个)}^{2} 我_{第页}$ ，它允许我们将问题（1）分解为第页可以并行解决的更简单的子问题。在这项工作中，我们重点讨论了具有滞后阶的VAR模型的情况(D类=1），以便强调研究联合估计过程的性质，而不是滞后阶选择的方面。联合估计方法首先假设每个VAR模型的公共和单独成分： ${A类}_{d日}^{(k个)} = {A类}_{d日, C类}^{(k个)} + {A类}_{d日, 我}^{(k个)}$ ,d日= 1, …,D、 k个= 1, …,K（K）之后，提出了一种迭代两阶段估计算法，该算法由一组套索优化程序组成，用于联合估计公共分量 ${{A类}_{d日, C类}^{(k个)}}$ 的K（K）受试者在第一阶段，然后使用稀疏套索优化程序来估计各个成分 ${{A类}_{d日, 我}^{(k个)}}$ 在第二阶段。群套索惩罚有效地将所有转移矩阵的各个元素分组K（K）并在模型中保留或排除整个组，从而保证生成的公共成分估计的共享结构。同时，公共分量信号的残差被用作数据，通过标准套索优化来估计单个结构，代表特定主题的效果。

本文的其余部分组织如下：第2节描述了联合建模框架并介绍了两阶段估计过程，第3节演示了各种设置下联合估计程序的模拟研究结果，并将其性能与其他最先进的方法进行了比较，第4节为ADHD研究的静止状态fMRI数据提供了所引入方法的大量实证应用，而第5节包含结束语和对未来工作的讨论。

2.问题制定

我们首先编写假定的VAR模型(1)以标准回归形式。首先，我们放弃k个从符号来看，k个= 1,…,K、，并显示单个VAR（D）模型所需的代数变换序列：

{X（X）}_{t吨} = {A类}_{1} {X（X）}_{t吨 - 1} + \dots + {A类}_{D类} {X（X）}_{t吨 - D类} + ϵ_{t吨}, ϵ_{t吨} ~ N个 (0, σ^{2} 我_{第页}), t吨 = D类, \dots, T型 .

(2)

误差协方差为对角cov的假设(ϵ_t吨) =σ²我_第页对于ϵ_t吨= (ϵ_1，吨, …,ϵ_{p、 t吨})^⊤，吨=D类, …,T型，允许我们表示每个第页变量为以下方程组：

{x个}_{j个, t吨} = \sum_{我 = 1}^{第页} ({A类}_{1} [j个, 我] {x个}_{我, t吨 - 1} + \dots + {A类}_{D类} [j个, 我] {x个}_{我, t吨 - D类}) + ϵ_{j个, t吨}, ϵ_{j个, t吨} ~ N个 (0, σ^{2}), t吨 = D类, \dots, T型, j个 = 1, \dots, 第页,

(3)

哪里A类_d日[j、我]是订单-d日时间效应我^第个变量打开j个^第个,l、 j个= 1, …,第页。如果我们让A类_d日[j个, .] = (A类_d日[j个, 1], …,A类_d日[j、第页])^⊤，天= 1, …,D类，然后全部T型−D类+1个方程式来自(三)可以用紧凑的矩阵形式表示每个变量j个,j个=1,…,第页，分别为：

\underset{(T型 - D类 + 1) \times 1}{\underset{︸}{{\tilde{X（X）}}_{j个}}} = \underset{(T型 - D类 + 1) \times (D类 第页)}{\underset{︸}{Z}} \underset{(D类 第页) \times 1}{\underset{︸}{A类 [j个, .]}} + \underset{(T型 - D类 + 1) \times 1}{\underset{︸}{{\tilde{ϵ}}_{j个},}} {\tilde{ϵ}}_{j个} ~ N个 (0, σ^{2} 我_{T型 - D类 + 1}) .

(4)

哪里

\underset{(T型 - D类 + 1) \times 1}{\underset{︸}{{\tilde{X（X）}}_{j个}}} = (\begin{matrix} {x个}_{j个, T型} \\ \dots \\ {x个}_{j个, D类} \end{matrix}), A类 [j个, .] = (\begin{matrix} {A类}_{1} [j个, .] \\ \dots \\ {A类}_{D类} [j个, .] \end{matrix}), {\tilde{ϵ}}_{j个} = (\begin{matrix} ϵ_{j个, T型} \\ \dots \\ ϵ_{j个, D类} \end{matrix}), Z = (\begin{matrix} {x个}_{1, T型 - 1} & \dots & {x个}_{第页, T型 - 1} \\ \dots & \dots & \dots \\ {x个}_{1, D类 - 1} & \dots & {x个}_{第页, D类 - 1} \end{matrix} | \begin{matrix} \dots \\ \dots \\ \dots \end{matrix} | \begin{matrix} {x个}_{1, T型 - D类} & \dots & {x个}_{第页, T型 - D类} \\ \dots & \dots & \dots \\ {x个}_{1, 0} & \dots & {x个}_{第页, 0} \end{matrix}) .

接下来，重新引入k、 k个= 1, …,K（K），回到符号并使用标准回归表示法(4)为所有人K（K）在考虑VAR模型时，我们获得（j个= 1, …,第页):

{\begin{array}{l} {\tilde{X（X）}}_{j个}^{(1)} = Z^{(1)} {A类}^{(1)} [j个, .] + {\tilde{ϵ}}_{j个}^{(1)}, & {\tilde{ϵ}}_{j个}^{(1)} ~ N个 (0, σ_{(1)}^{2} 我_{T型 - D类 + 1}), \\ \dots \\ {\tilde{X（X）}}_{j个}^{(K（K）)} = Z^{(K（K）)} {A类}^{(K（K）)} [j个, .] + {\tilde{ϵ}}_{j个}^{(K（K）)}, & {\tilde{ϵ}}_{j个}^{(K（K）)} ~ N个 (0, σ_{(K（K）)}^{2} 我_{T型 - D类 + 1}), \end{array}

(5)

哪里A类^(k个)[j个, .] 与时间效应相对应D类订单）每个第页变量在j个^第个变量用于k个^第个主题，k个= 1, …,K（K）.

2.1. 分解为通用和不同步组件

一个关键的建模假设是跨所有K（K）学科。在模型中，这通过跨K（K）转换矩阵。此外，由于受试者之间的自然变异性，我们必须以特定受试者特定效应的形式解释异质性的存在。通过将每个主题的转换矩阵分解为两部分，可以捕获这两个方面：

{A类}_{d日}^{(k个)} = {A类}_{d日, C类}^{(k个)} + {A类}_{d日, 我}^{(k个)}, d日 = 1, \dots, D类, k个 = 1, \dots, K（K）,

(6)

哪里 ${A类}_{d日, C类}^{(k个)}$ 是订单的常见组成部分-d日时间效应k个^第个主题，而 ${A类}_{d日, 我}^{(k个)}$ 是特殊成分。将此表示应用于中的方程式(5)，我们得到（j个= 1, …,第页):

{\begin{array}{l} {\tilde{X（X）}}_{j个}^{(1)} = Z^{(1)} & ({A类}_{C类}^{(1)} [j个, .] + {A类}_{我}^{(1)} [j个, .]) + {\tilde{ϵ}}_{j个}^{(1)}, & {\tilde{ϵ}}_{j个}^{(1)} ~ N个 (0, σ_{(1)}^{2} 我_{T型 - D类 + 1}), \\ \dots \\ {\tilde{X（X）}}_{j个}^{(K（K）)} = Z^{(K（K）)} & ({A类}_{C类}^{(K（K）)} [j个, .] + {A类}_{我}^{(K（K）)} [j个, .]) + {\tilde{ϵ}}_{j个}^{(K（K）)}, & {\tilde{ϵ}}_{j个}^{(K（K）)} ~ N个 (0, σ_{(K（K）)}^{2} 我_{T型 - D类 + 1}) . \end{array}

(7)

假设每个主题的VAR模型是公共基础VAR模型的扰动，我们继续在 ${A类}_{C类}^{(1)} [j个, .], \dots, {A类}_{C类}^{(K（K）)} [j个, .]$ ，表示为 ${A类}_{C类}^{(1)} [j个, .] \approx {A类}_{C类}^{(2)} [j个, .] \approx \dots \approx {A类}_{C类}^{(K（K）)} [j个, .]$ ，j= 1, …,第页更正式地定义向量的支持 $β = {(β_{1}, \dots, β_{米})}^{⊤} \in ℝ^{米}$ 作为支持(β) = {我:β_我≠0}这个假设等价于 $秒 u个第页第页哦第页 t吨 ({A类}_{C类}^{(1)} [j个, .]) \equiv 秒 u个第页第页哦第页 t吨 ({A类}_{C类}^{(2)} [j个, .]) \equiv \dots \equiv 秒 u个第页第页哦第页 t吨 ({A类}_{C类}^{(K（K）)} [j个, .])$ 此外，稀疏性假设被施加在各个组件上 ${A类}_{我}^{(1)} [j个, .], \dots, {A类}_{我}^{(K（K）)} [j个, .]$ 恢复最重要的特定主题效果。最后，对于参数可识别性，我们假设 ${A类}_{C类}^{(k个)} [j个, .] ⊥ {A类}_{我}^{(k个)} [j个, .] ({A类}_{C类}^{(k个)} [j个, .]$ “垂直”于 ${A类}_{我}^{(k个)} [j个, .])$ 这意味着 ${A类}_{C类}^{(k个)} [j个, .]$ 和 ${A类}_{我}^{(k个)} [j个, .]$ 是空的，或者更正式地说，支持 $({A类}_{C类}^{(k个)} [j个, .]) \cap$ 支持 $({A类}_{我}^{(k个)} [j个, .]) \equiv 0$ ,k个= 1, …,K（K）,j个= 1, …,第页系统的完整约束(7)概述如下：

{\begin{array}{l} {A类}_{C类}^{(1)} [j个, .] \approx {A类}_{C类}^{(2)} [j个, .] \approx \dots & \approx {A类}_{C类}^{(K（K）)} [j个, .], j个 = 1, \dots, 第页, \\ {A类}_{我}^{(k个)} [j个, .] - 稀疏 & k个 = 1, \dots, K（K）, j个 = 1, \dots, 第页, \\ {A类}_{C类}^{(k个)} [j个, .] ⊥ {A类}_{我}^{(k个)} [j个, .], & k个 = 1, \dots, K（K）, j个 = 1, \dots, 第页 . \end{array}

(8)

2.2. 两阶段估算程序

估计参数 ${A类}_{C类}^{(k个)} [j个, .]$ , ${A类}_{我}^{(k个)} [j个, .]$ ，千= 1, …,K、 j个= 1, …,第页，在中显示(7)，同时执行上述一组约束(8)，我们将其表示为一个双凸优化任务，该任务可以通过迭代两阶段方法求解，该方法始终保证收敛到局部极小值[18,17]. 首先，我们要处理令人讨厌的参数{σ_(k个)，千= 1, …,K（K）}通过插入最大似然估计量 ${{\hat{σ}}_{(k个)}, k个 = 1, \dots, K（K）}$ 计算方法如下[26]对于系统方程式（1）。接下来，我们设置 $β_{j个, C类} = {({A类}_{C类}^{(1)} [j个, .], \dots, {A类}_{C类}^{(K（K）)} [j个, .])}^{⊤}$ , $β_{j个, 我} = {({A类}_{我}^{(1)} [j个, .], \dots, {A类}_{我}^{(K（K）)} [j个, .])}^{⊤}$ , ${\tilde{X（X）}}_{j个} = {({\tilde{X（X）}}_{j个}^{(1)}, \dots, {\tilde{X（X）}}_{j个}^{(K（K）)})}^{⊤}$ ,j个= 1, …,第页, ${\hat{D类}}_{{\hat{σ}}^{2}} = 诊断 (\underset{T型 - D类 + 1}{\underset{︸}{{\hat{σ}}_{(1)}^{2}, \dots, {\hat{σ}}_{(1)}^{2}}}, \underset{T型 - D类 + 1}{\underset{︸}{{\hat{σ}}_{(2)}^{2}, \dots, {\hat{σ}}_{(2)}^{2}}}, \dots, \underset{T型 - D类 + 1}{\underset{︸}{{\hat{σ}}_{(K（K）)}^{2}, \dots, {\hat{σ}}_{(K（K）)}^{2}}})$ 。请注意 $\tilde{Z} \in ℝ^{K（K） (T型 - D类 + 1) \times K（K） (D类第页)}$ 是一个块-对角矩阵k个^第个块等于 $Z^{(k个)} \in ℝ^{(T型 - D类 + 1) \times (D类第页)}$ 从方程式（5）,k个= 1, …,英国。此外，对于任意向量 $β = (β_{1}, \dots, β_{米}) \in ℝ^{米}$ 让我们表示 $‖ β ‖_{2}^{2} = \sum_{我 = 1}^{米} β_{我}^{2}, ‖ β ‖_{2} = \sqrt{\sum_{我 = 1}^{米} β_{我}^{2}}, ‖ β ‖_{1} = \sum_{我 = 1}^{米} | β_{我} |$ .与求解相对应的大型优化问题方程式（7）带有约束(8)如下所示

\underset{β_{j个, C类}, β_{j个, 我}}{最小值} （f） (β_{j个, C类}, β_{j个, 我}) = \underset{β_{j个, C类}, β_{j个, 我}}{最小值} {‖ {\hat{D类}}_{{\hat{σ}}^{2}}^{- \frac{1}{2}} ({\tilde{X（X）}}_{j个} - \tilde{Z} [β_{j个, C类} + β_{j个, 我}]) ‖}_{2}^{2} + λ_{j个}^{G公司} \sum_{我 = 1}^{第页} {‖ {(β_{j个, C类}^{(1)} [我], \dots, β_{j个, C类}^{(K（K）)} [我])}^{⊤} ‖}_{2} + + λ_{j个}^{秒} {‖ β_{j个, 我} ‖}_{1} + λ^{\infty} \sum_{我 = 1}^{第页} \sum_{k个 = 1}^{K（K）} | β_{j个, C类}^{(k个)} [我] \cdot β_{j个, 我}^{(k个)} [我] |,

(9)

哪里 $β_{j个, C类}^{(k个)} [我] = {A类}_{C类}^{(k个)} [j个, 我]$ , $β_{j个, 我}^{(k个)} [我] = {A类}_{我}^{(k个)} [j个, 我]$ ,k个= 1, …,K（K）,我= 1, …,第页。这里，来自(8)通过相应的惩罚条款来解决。集体套索惩罚 $λ_{j个}^{G公司} \sum_{我 = 1}^{第页} {‖ {(β_{j个, C类}^{(1)} [我], \dots, β_{j个, C类}^{(K（K）)} [我])}^{T型} ‖}_{2}$ ，引入于[45]，要么缩小我^第个元素全部为零K（K）向量 $β_{j个, C类}^{(1)}, \dots, β_{j个, C类}^{(K（K）)}$ ，其中 $β_{j个, C类}^{(k个)} = {A类}_{C类}^{(k个)} [j个, .]$ ，或估计它对所有人来说都是非零的K（K）矢量。这保证了所有K（K）共同成分估计， ${A类}_{C类}^{(1)} [j个, .] \approx {A类}_{C类}^{(2)} [j个, .] \approx \dots \approx {A类}_{C类}^{(K（K）)} [j个, .]$ 稀疏套索惩罚 $λ_{j个}^{秒} | | β_{j个, 我} ‖_{1}$ ，首次引入于[38]，指向每个 ${A类}_{我}^{(k个)} [j个, .]$ ,k个= 1, …,K（K）,j个= 1, …,第页，只有几个非零元素，因此很稀疏。罚款期限 $λ^{\infty} \sum_{我 = 1}^{第页} \sum_{k个 = 1}^{K（K）} | β_{j个, C类}^{(k个)} [我] \cdot β_{j个, 我}^{(k个)} [我] |$ 具有调谐参数值λ^∞设置足够高，以便哈达玛产品β_{j、 C类}和β_{j、我}等于0。因此，它导致 $β_{j个, C类}^{(k个)} [我] \neq 0$ 暗示 $β_{j个, 我}^{(k个)} [我] = 0$ ，反之亦然， $β_{j个, 我}^{(k个)} [我] \neq 0$ 暗示 $β_{j个, C类}^{(k个)} [我] = 0$ .这一保证支持 $(β_{j个, C类}^{(k个)}) \cap$ 支持 $(β_{j个, 我}^{(k个)}) \equiv 0$ ，导致 ${A类}_{C类}^{(k个)} [j个, .] ⊥ {A类}_{我}^{(k个)} [j个, .] \forall j个, k个$ ,j个= 1, …,第页,k个= 1, …,K（K）.

对于关节功能（f）(β_{j、 C类},β_{j、我})正在进行优化(9)，使用凸性的定义和保凸运算的知识[6]，我们可以显示：修复β_{j、我}以某种价值 ${\hat{β}}_{j个, 我}, （f） (β_{j个, C类}, {\hat{β}}_{j个, 我})$ 是的凸函数β_{j、 C类}; 固定β_{j、 C}以某种价值 ${\hat{β}}_{j个, C类}$ , $（f） ({\hat{β}}_{j个, C类}, β_{j个, 我})$ 是的凸函数β_{j、我}根据定义[18]，它导致（f）(β_{j、 C类},β_{j、我})构成其自变量的双凸函数β_{j、 C类}和β_{j、我}。下面我们给出了一个两阶段算法，该算法执行另一种凸搜索方法[18,43]，其中，对于一般双凸函数（f）(x个,年)，一个替代更新x个和年按以下方式：修复年在初始化值时 $\hat{年}, 年 \equiv \hat{年}$ 求解凸优化问题x个, $\hat{x个} = \underset{x个}{argmin（最小值）} （f） (x个, \hat{年})$ ; 修理x个以优化器值 $\hat{x个}$ , $x个 \equiv \hat{x个}$ ; 求解凸优化问题年, $\hat{年} = \underset{年}{argmin（最小值）} （f） (\hat{x个}, 年))$ ; 等等。

两阶段估计算法(用于任意 j、 j个= 1, …,第页)

0.初始化 ${\hat{β}}_{j个, 我}$ 带有零矢量， ${\hat{β}}_{j个, 我} \equiv 0$ .

1第一阶段：继续解决 $\underset{β_{j个, C类}}{米我 n个} （f） (β_{j个, C类}, {\hat{β}}_{j个, 我})$ ，相当于求解以下凸群lasso优化准则

\underset{β_{j个, C类}}{最小值} {‖ {\hat{D类}}_{{\hat{σ}}^{2}}^{- \frac{1}{2}} ([{\tilde{X（X）}}_{j个} - \tilde{Z} {\hat{β}}_{j个, 我}] - {\tilde{Z}}^{(C类)} β_{j个, C类}) ‖}_{2}^{2} + λ_{j个}^{G公司} \sum_{我 = 1}^{第页} {‖ {(β_{j个, C类}^{(1)} [我], \dots, β_{j个, C类}^{(K（K）)} [我])}^{⊤} ‖}_{2},

(10)

其中矩阵 ${\tilde{Z}}^{(C类)}$ 是这样的 ${\tilde{Z}}^{(C类)} [., 我] \equiv \tilde{Z} [, 我]$ 如果 ${\hat{β}}_{j个, 我} [我] = 0$ 、和 ${\tilde{Z}}^{(C类)} [., 我] \equiv 0$ 如果 ${\hat{β}}_{j个, 我} [我] \neq 0$ .使用此类矩阵 ${\tilde{Z}}^{(C类)}$ 保证优化器 ${\hat{β}}_{j个, C类} \equiv \underset{β_{j个, C类}}{argmin（最小值）} （f） (β_{j个, C类}, {\hat{β}}_{j个, 我})$ 将具有非重叠支持 ${\hat{β}}_{j个, 我}, {\hat{β}}_{j个, C类} ⊥ {\hat{β}}_{j个, 我}$ 。我们这样做，而不是明确包括罚款条款 $λ^{\infty} \sum_{我 = 1}^{第页} \sum_{k个 = 1}^{k个} | β_{j个, C类}^{(k个)} [我] \cdot β_{j个, 我}^{(k个)} [我] |$ 来自(9)简化了优化过程。

2第二阶段：让 ${\hat{β}}_{j个, C类}$ 表示对β_{j、 C类}从第一阶段开始。继续解决最小值 $\underset{β_{j个, 我}}{米我 n个} （f） ({\hat{β}}_{j个, C类}, β_{j个, 我})$ ，它等价于以下凸套索问题

\underset{β_{j个, 我}}{最小值} {‖ {\hat{D类}}_{{\hat{σ}}^{2}}^{- \frac{1}{2}} ([{\tilde{X（X）}}_{j个} - {\tilde{Z}}^{(我)} {\hat{β}}_{j个, C类}] - \tilde{Z} β_{j个, 我}) ‖}_{2}^{2} + λ_{j个}^{秒} {‖ β_{j个, 我} ‖}_{1},

(11)

其中矩阵 ${\tilde{Z}}^{(我)}$ 是这样的 ${\tilde{Z}}^{(我)} [., 我] \equiv \tilde{Z} [, 我]$ 如果 ${\hat{β}}_{j个, C类} [我] = 0$ 、和 ${\tilde{Z}}^{(我)} [., 我] \equiv 0$ 如果 ${\hat{β}}_{j个, C类} [我] \neq 0$ 类似于 ${\tilde{Z}}^{(C类)}$ 在第一阶段，使用这样的矩阵 ${\tilde{Z}}^{(我)}$ 保证 ${\hat{β}}_{j个, C类} ⊥ {\hat{β}}_{j个, 我}$ .

3.如果是第二次迭代（或更高）：表示 ${\hat{β}}_{j个} = {\hat{β}}_{j个, C类} + {\hat{β}}_{j个, 我}$ 作为当前迭代的完整估计，以及 ${\hat{β}}_{j个}^{第页第页}$ -前一次迭代的完全估计；如果出现以下情况，则停止算法 ${‖ {\hat{β}}_{j个} - {\hat{β}}_{j个}^{第页第页} ‖}_{2}^{2} < 10^{- 4}$ 。否则，返回第一阶段，使用 ${\hat{β}}_{j个, 我}$ 在第二阶段进行计算。

所述的两阶段迭代方法执行交替的凸搜索方法，这将导致双凸函数的局部优化器（f）(β_{j、 C类},β_{j、我})调整参数的固定值 $λ_{j个}^{G公司}$ 和 $λ_{j个}^{秒}$ 。的特定值 $λ_{j个}^{G公司}$ 和 $λ_{j个}^{秒}$ 反过来，可以通过启发式方法进行选择，即将上述两阶段算法运行几次迭代，但向两个阶段中的每个阶段添加一个调整参数选择步骤，如下所示：

1第一阶段（续）：计算的求解路径 $λ_{j个}^{G公司}$ 优化任务的值(10)如中所示[45]. 选择调谐参数值 ${\hat{λ}}_{j个}^{G公司}$ 最小化贝叶斯信息标准（BIC）

B类 我 C类 (λ_{j个}) = n个 日志 (‖ {\hat{D类}}_{{\hat{σ}}^{2}}^{- \frac{1}{2}} ([{\tilde{X（X）}}_{j个} - \tilde{Z} {\hat{β}}_{j个, 我}] - \tilde{Z} {\hat{β}}_{j个, C类} (λ_{j个}^{G公司}))] ‖_{2}^{2} + 日志 (n个) {数据流}_{λ_{j个}^{G公司}},

(12)

哪里n个=K（K）(T型−D类+ 1), ${\hat{β}}_{j个, C类} (λ_{j个}^{G公司})$ -的估计β_{j、 C类}对应于值 $λ_{j个}^{G公司}$ 解决方案路径的， ${数据流}_{λ_{j个}^{G公司}}$ -估计自由度 ${\hat{β}}_{j个, C类} (λ_{j个}^{G公司})$ ，计算如下[7].

2第二阶段（续）：获取的解决方案路径 $λ_{j个}^{秒}$ 优化任务的值(11)如中所示[39]与第一阶段类似，选择调整参数值 ${\hat{λ}}_{j个}^{G公司}$ 使BIC最小化。自由度 ${数据流}_{λ_{j个}^{秒}}$ 在这种情况下，计算为估计中非零元素的数量 ${\hat{β}}_{j个, 我} (λ_{j个}^{秒})$ .

我们在最初的几次迭代中实现了这种两阶段方法的扩展版本（在所介绍的数值计算中，我们在前五次迭代中使用了它），以便获得良好的数据驱动的调整参数值选择 ${\hat{λ}}_{j个}^{G公司}$ 和 ${\hat{λ}}_{j个}^{秒}$ 然后，设置 $λ_{j个}^{G公司} \equiv {\hat{λ}}_{j个}^{G公司}$ 和 $λ_{j个}^{秒} \equiv {\hat{λ}}_{j个}^{秒}$ 用于优化任务(9) – (11)，我们只需执行原始版本的两阶段算法，该算法保证收敛到（f）(β_{j、 C类},β_{j、我}). 为了提高在全局最小值中着陆的可能性，我们可能需要查看多个不同的初始化值进行估计 ${\hat{β}}_{j个, 我}$ 在我们两阶段方法的第0步，这可能是未来工作的一个有趣的主题（请参阅中的相关工作[23]).

3.绩效评估

我们的评估研究重点是V AR（车辆配置总成）订单模型D类=1，以评估拟议联合估算程序的性能，而不是侧重于模型订单选择：

{X（X）}_{t吨}^{(k个)} = {A类}^{(k个)} {X（X）}_{t吨 - 1}^{(k个)} + ϵ_{t吨}^{(k个)}, ϵ_{t吨}^{(k个)} ~ N个 (0, σ_{(k个)}^{2} 我_{第页}), t吨 = 1, \dots, T型, k个 = 1, \dots, K（K） .

(13)

采用的绩效指标如下所示，其中 $\hat{A类} = {({\hat{一}}_{我, j个})}_{第页 \times 第页}$ 表示转移矩阵估计，以及 $A类 = {(一_{我, j个})}_{第页 \times 第页}$ -真值矩阵：

假阳性（FP）和真阴性（TN）率：
$F类 P（P） = \frac{\sum_{1 \leq 我, j个 \leq 第页} 我 (一_{我, j个} = 0, {\hat{一}}_{我, j个} \neq 0)}{\sum_{1 \leq 我, j个 \leq 第页} 我 (一_{我, j个} = 0)}, T型 N个 = 1 - F类 P（P） .$
假阴性（FN）和真阳性（TP）率：
$F类 N个 = \frac{\sum_{1 \leq 我, j个 \leq 第页} 我 (一_{我, j个} \neq 0, {\hat{一}}_{我, j个} = 0)}{\sum_{1 \leq 我, j个 \leq 第页} 我 (一_{我, j个} \neq 0)}, T型 P（P） = 1 - F类 N个 .$
马修斯相关系数（MC，FP和FN的几何平均值）：
$M（M） C类 = \frac{T型 P（P） \times T型 N个 - F类 P（P） \times F类 N个}{\sqrt{(T型 P（P） + F类 P（P）) (T型 P（P） + F类 N个) (T型 N个 + F类 P（P）) (T型 N个 + F类 N个)}} .$
归一化Frobenius差（NFD）：
$N个 F类 D类 = \frac{‖ \hat{A类} - A类 ‖_{F类}}{‖ A类 ‖_{F类}} = \frac{\sqrt{\sum_{1 \leq 我, j个 \leq 第页} {({\hat{一}}_{我, j个} - 一_{我, j个})}^{2}}}{\sqrt{\sum_{1 \leq 我, j个 \leq 第页} 一_{我, j个}^{2}}} .$

前三个度量值（FP、FN、MC）描述了矩阵结构（非零元素的定位）的估计情况，这是我们在本工作中的优先考虑事项。FP、F N的低值（接近0）和MC的高值（接近1）表示性能良好。同时，NFD评估矩阵元素大小的估计程度。NFD值越低（接近0），性能越好。当仅评估中引入的两阶段估计算法的性能时第2节，所有四个度量值都是针对公共分量估计值计算的 ${\hat{A类}}_{C类}^{(k个)}$ 和单个组件估算 ${\hat{A类}}_{我}^{(k个)}$ ,k个= 1, …,K（K）同时，当将此两阶段方法与其他最先进的方法进行比较时，由于这些其他方法没有实现对通用和单个组件的分解，因此我们仅计算了完整矩阵估计值的四个度量值 ${\hat{A类}}^{(k个)} = {\hat{A类}}_{C类}^{(k个)} + {\hat{A类}}_{我}^{(k个)}$ ,k个= 1, …,K（K）.

模拟研究的设计如下：矩阵 ${{A类}^{(k个)} = {A类}_{C类}^{(k个)} + {A类}_{我}^{(k个)}, k个 = 1, \dots, K（K）}$ 以0.4的光谱半径生成，非零效应的幅度至少为0.2。这将保证生成的时间序列的平稳性，从信号中分离噪声，并在很大程度上复制fMRI数据的特征第4节（没有估计的影响大于0.4）。通过添加独立的N个（0，1）错误，并将信噪比设置为1（如ADHD研究中第4节). 通过改变时间序列的数量来考虑多种设置第页每个受试者和受试者数量K（K）在一个组中。生成的所有对角线都是非零的，而公共分量中元素的边缘密度被设置为（o对角线元素的）5%第页=20，2%第页= 30, 40. 单个组件的边缘密度设置为矩阵元素总数的2–3%，这意味着存在中等程度的异质性。将0.02的阈值应用于结果估计的元素，以消除噪声条目。最后，通过平均50次以上的重复获得指标。

鉴于我们的联合评估程序是以跨多个学科借力为目标而开发的，当逐渐增加数量时，其效率最佳K（K）相关科目。关于时间序列长度增加情况下的模拟设置T型，见附录中的结果。在表1下面，我们演示了当数字K（K）联合评估的受试者从10人增加到50人。随着每组受试者人数的增加，我们发现对共同成分结构的估计有了明显的改进（FP减少，MC增加），这也导致了对个体效应的更好估计（FN增加，MC减少）。这方面源于对公共和单个组件施加的可识别性约束：公共组件中的每个误报可能对应于单个组件中错误放置的真阳性，从而直接导致其中的误报。因此，在估计公共成分时，误报率越低，导致单个成分估计的误报率就越低。尽管如此，单个组件的FN数并不是非常接近于零。这是强制执行以下假设的直接结果(8)单个组件非常稀疏。放松这一假设，尽管FN降低，但会导致单个组件的FP增加，根据所采用的可识别性假设，这会导致缺失真正的组级效应（见补充确认这一说法的额外模拟结果）。鉴于我们优先考虑正确检测群体层面的影响，并且只想说明最关键的特定主题影响，因此最好在单个组件中使用FN，而不是在公共组件中使用。此外，它可能源于在两阶段方法中使用零矢量初始化单个组件，从而导致联合函数的局部最小值(9). 另一方面，当考虑总体估算的好坏时 $\hat{A类} = {\hat{A类}}_{C类} + {\hat{A类}}_{我}$ 与真矩阵有关A类，您可能会注意到，对于某些设置，FP和FN都接近于0（请参见表2其中，我们将两阶段方法与三种竞争方法进行了比较）。由于它属于算法收敛特性，在完成调谐参数选择迭代后，算法对每个设置进行收敛，在50个重复中平均需要2-3次迭代。

表1：

增加受试者群体规模的两阶段程序表现K（K）.

K（K）	FP（通信）	FN（通信）	MC（通信）	NFD（通信）	FP（指示）	FN（指示）	MC（印度）	NFD（印度）
				第页=20	T型=80
10	0.06(0.02)	0.05(0.08)	0.89(0.06)	0.5(0.18)	0(0)	0.92(0.07)	0.19(0.08)	0.98(0.02)
20	0.03(0.02)	0.06(0.08)	0.92(0.06)	0.51(0.19)	0(0)	0.65(0.06)	0.46(0.05)	0.89(0.02)
50	0.01(0.01)	0.06(0.1)	0.94(0.08)	0.51(0.21)	0(0)	0.43(0.11)	0.63(0.08)	0.8(0.05)
				第页=30	T型=120
10	0.06(0.01)	0(0)	0.94(0)	0.35(0.04)	0(0)	0.95(0.03)	0.15(0.04)	0.99(0.01)
20	0.02(0)	0(0)	0.98(0)	0.3(0.01)	0(0)	0.51(0.04)	0.56(0.03)	0.83(0.02)
50	0.01(0)	0(0)	0.99(0)	0.3(0.07)	0(0)	0.18(0.02)	0.83(0.02)	0.66(0.01)
				第页=40	T型=150
10	0.06(0.01)	0(0)	0.94(0.01)	0.39(0.1)	0(0)	0.96(0.04)	0.12(0.07)	0.99(0.01)
20	0.02(0.01)	0(0)	0.97(0)	0.34(0.11)	0(0)	0.5(0.09)	0.57(0.07)	0.83(0.04)
50	0.01(0)	0(0)	0.99(0)	0.29(0.1)	0(0)	0.19(0.04)	0.82(0.04)	0.64(0.01)

在单独的窗口中打开

注：行的三元组对应于相同的设置第页和T型例如，第3-5行对应于以下情况第页= 20，T型=80，第7-9行对应设置第页= 30，T型=120等第3节，每个设置有50个重复。

表2：

四种方法的比较模拟研究：规则套索（lasso）、稀疏群套索（SGL）、群桥（GB）和我们的两阶段方法（2ST）。

设置	方法	FP公司	FN公司	国会议员	非金融机构
	拉索	0.02(0.02)	0.42(0.42)	0.61(0.33)	0.71(0.21)
K（K）= 20	SGL公司	0.3(0.05)	0.05(0.08)	0.68(0.04)	0.55(0.11)
第页= 20,T型= 80.	GB（英国）	0.03(0.01)	0.13(0.09)	0.85(0.07)	0.5(0.16)
	第2阶段	0.02(0.01)	0.14(0.08)	0.85(0.06)	0.62(0.16)
	拉索	0.01(0)	0.11(0.18)	0.9(0.15)	0.49(0.14)
K（K）= 20,	SGL公司	0.2(0.03)	0.02(0.04)	0.79(0.02)	0.45(0.08)
第页= 40,T型= 150.	GB（英国）	0.02(0)	0.03(0.02)	0.95(0.02)	0.32(0.08)
	第二	0.02(0.01)	0.03(0.01)	0.95(0.01)	0.49(0.08)
	拉索	0.02(0.02)	0.43(0.41)	0.61(0.32)	0.71(0.21)
K（K）= 50,	SGL公司	0.36(0.05)	0.04(0.05)	0.64(0.04)	0.55(0.11)
第页= 20,T型= 80.	GB（英国）	0.03(0.01)	0.17(0.13)	0.81(0.1)	0.54(0.18)
	第2阶段	0.01(0.01)	0.13(0.1)	0.87(0.08)	0.61(0.17)
	拉索	0.01(0)	0.08(0.17)	0.92(0.14)	0.47(0.13)
K（K）= 50,	SGL公司	0.23(0.03)	0.02(0.02)	0.76(0.01)	0.45(0.05)
第页= 40,T型= 150.	GB（英国）	0.03(0)	0.05(0.02)	0.92(0.02)	0.35(0.08)
	第2阶段	0.01(0)	0.02(0)	0.97(0)	0.42(0.07)

在单独的窗口中打开

注：每四行对应于应用于相同设置的四种方法。例如，第2-5行对应于设置第页= 20，T型= 80，K=20，第6-9行对应设置第页= 40，T型= 150，K=20等第3节，每个设置50个重复。

为了评估拟议联合建模方法与其他既定正则化估计方法相比的性能，我们进行了比较仿真研究，包括：

常规拉索方法[38]仅使用稀疏惩罚标准‖β‖₁= ∑我|β_我|以分别获得每个受试者的稀疏估计，而不考虑分组方面。
稀疏群套索[36]，它在常规组套索的顶部添加了一个单变量稀疏惩罚项[45]从而在群体和个人层面上都考虑到稀疏性。虽然听起来与我们提出的方法类似，但稀疏组套索并没有将其分解为通用和单个组件。相反，它将组级和个体级稀疏性惩罚直接应用于原始VAR转换矩阵。
团体桥梁[21]，它实现了桥接型组级惩罚规范 $‖ β ‖_{1}^{γ} = {(\sum_{我} | β_{我} |)}^{γ}$ , 0<γ<1，而不是经典的欧几里德范数 $‖ β ‖_{2} = \sqrt{\sum_{我} β_{我}^{2}}$ 在原来的组套索中使用[45]. 类似地，稀疏组套索[36]，它不会导致VAR转换矩阵分解为通用和单个组件。

常规套索和群桥的调整参数值是通过BIC选择的，正如我们在两阶段方法中所做的那样。同时，对于稀疏组套索，由于自由度公式不可用，我们使用了交叉验证。表2下面演示了所有四种方法（常规套索、稀疏组套索、组桥、两阶段方法）在估计完整VAR转移矩阵时的性能A类.

就结构估计（FP、FN、MC）而言，对于具有K（K）=20名受试者，组桥法的表现与我们的两阶段方法相当，而其他两种方法则远远落后。同时，当将受试者数量增加到K（K）=50，我们发现所提议的两阶段方法明显优于所有竞争对手，因为它在许多类似学科中都具有优势。对于K（K）= 50，第页= 20，T型=80和K（K）= 50，第页= 40，T型分别=150个设置。另一方面，就单元幅度估计（通过NFD测量）而言，群桥算法比所有其他方法都具有优势。然而，在K（K）= 50，第页= 40，T型=150我们的两阶段方法在NFD中优于常规套索和稀疏组套索，并且在所有方法中已经提供了最佳的结构估计质量。

4.大脑fMRI数据在ADHD研究中的应用

我们分析了两组患者的静止状态fMRI数据：一组是健康对照组，另一组是表现出某种认知障碍的患者。预计两组受试者的大脑区域之间的前导-滞后关系模式相似，同时也表现出特定的受试者效应。共享结构的思想以前曾用于fMRI数据（参见[11], [5])但重点是大脑区域之间的同期依赖性（功能连通性）。我们的估算程序试图描述这些区域的时间动态（有效连通性）。

注意，脑内跨区域关系的VAR建模最近受到了文献的关注。尽管已经有工作应用VAR模型来评估大脑的有效连通性([15], [16])，已经指出了这种方法的各种缺点。例如，一篇关于静止状态fMRI分析的进展和缺陷的综述性论文[12]强调各区域血流动力学延迟的变化，这可能会在任何估算时间效应的尝试中引入偏差。与此同时[32]在从fMRI数据推断时间因果关系时，继续列出了六个谨慎的详细理由，其中之一是不同实验地点之间的潜在异质性，以及由于神经活动过程发生得比fMRI测量的采样率更快而无法捕获因果关系。然而，我们的建模方法缓解了上面讨论的一些问题。首先，套索组直接解决了[32]通过捕捉常见的抽象处理结构，例如大脑的哪些区域影响其他区域，同时考虑到这些影响在患者中的不同强度。作为提醒，它选择了一种对所有患者都常见的特定关系，但确切的影响程度可能会因患者而异。我们还继续选择功能磁共振成像测量重复时间相同的研究，否则可能会联合估计不同延迟的时间效应。此外，个别成分还说明了特定主题的影响，例如年龄或惯用手。最后，为了避免实验产生不同受试者的不同感兴趣区域（ROI），我们确保在所有患者中使用标准的统一脑图谱进行区域分配。

4.1. ADHD研究设置

我们考虑了20名ADHD患者和20名对照组的静止状态脑fMRI时间序列数据，这些数据来自ADHD-200全球竞赛，可以使用Python模块检索尼罗河受试者来自五个实验点，其中三个（纽约大学儿童研究中心、北京大学和Radboud大学神经影像学研究）的TR为2.0秒，另外两个（俄勒冈州健康与科学大学、肯尼迪克里格研究所）TR更长，为2.5秒。由于我们的模型旨在推断大脑中的时间效应，因此我们将最后两项研究排除在考虑范围之外，以便在所有受试者中具有一致的测量重复时间。这就给我们留下了12名ADHD受试者（所有男性；年龄平均值±SD=13.85±3.83）和12名对照者（所有雄性；年龄平均数±SD=13.72±3.72）。所有站点的信噪比均为1。数据预处理步骤包括校正切片采集和运动的延迟、滤波以去除高频、数据标准化和去趋势（有关更多详细信息，请参阅[41]). 如上所述，我们根据多学科词典学习图谱（MSDL）将大脑分成39个区域，然后按照[41]. 我们通过体素时间序列的平均值总结这些区域上的信号，并通过灰度概率分割进行加权。区域的解剖位置和由此提取的时间序列都可以在图1和和2，2，详细介绍了包含在表3属于附录6.1此外，图3描述了对应于为其中一个区域提取的时间序列的自相关图。PACF（右下角图）第一个时间滞后时的明显峰值是一个强大的AR（1）信号，它存在于绝大多数考虑中的大脑区域。这使我们相信VAR（1）模型就足够了，这也激发了在第3节.

保存图片、插图等的外部文件。对象名为nihms-1530078-f0001.jpg

在单独的窗口中打开

图1：

根据MSDL图谱的大脑区域（参见表3用于区域名称）。

保存图片、插图等的外部文件。对象名为nihms-1530078-f0002.jpg

提取单个受试者39个脑区的时间序列。

单区域时间序列的自相关图：顶部时间序列图、左下角自相关图（ACF）、右下角部分自相关性图（PACF）。

表3：

MSDL图谱脑区名称。

地区#	地区名称	地区#	地区名称
1	L听觉皮层	21	初级视皮层
2	R听觉皮层	22	R Lat Occ复合体
三	终纹	23	蚂蚁背皮层
4	L默认模式网络	24	腹侧Ant Cing皮层
5	Med Prefr皮层	25	R蚂蚁岛皮层
6	Fr默认模式网络	26	L Sup-Temp沟
7	R默认模式网络	27	R Sup温度沟
8	Occ叶	28	L Temp-Par接头
9	运动皮层	29	前额叶Broca区
10	R Dorso-Lat Prefr皮层	30	Sup-Fr公司
11	右极前叶	31	R Temp-Par接头
12	R Par叶	32	Opercularis部分
13	R Inf温度皮层	33	小脑
14	基底神经节	34	背侧后皮层
15	L Par叶	35	L岛叶皮层
16	L多尔索·拉特	36	Cing Cortex公司
17	L Polar Fr叶	37	R岛叶皮层
18	L车内沟	38	L Ant手掌内侧
19	R车内沟	39	R Ant车内沟
20	L Lat Occ综合体

在单独的窗口中打开

注：常用缩写为“L”-左、“R”-右、“Lat”-侧、“Ant”-前、“Fr”-前，“Med”-中、“Post”-后、“Sup”-上、“Inf”-下、“Temp”-颞、“Par”-顶骨、“Occ”-枕、“Cing”-扣带。

考虑到缺乏关于分组套索程序估计值的分布特征和渐近性质的文献，我们采用稳定性选择来评估大脑区域之间估计时间效应的一致性。稳定性选择的概念在[28]其主要思想是使用bootstrap，然后将结果累加到所有bootstrapped样本上，以供进一步分析。在有限样本量的多重测试中，子采样导致控制族I型错误率，这对于高维问题来说比观测数量趋于无穷大的渐近语句重要得多。在处理时间序列数据时，必须生成子样本，同时保留观测值之间的相关性结构。时间序列块自举技术([20], [8])负责这项任务，为我们提供B类=沿途每名受试者100次评估。它允许评估其稳定性，并导致对大脑区域之间的时间关系进行更可靠的描述性分析。

4.2. ADHD研究结果

我们从第2.2节为了联合估计同一组患者跨大脑区域的时间依赖性，并将我们的注意力集中在共同成分估计上。图4下面以有向图的形式描述了大脑区域之间产生的时间关系：左侧39个节点对应于每个大脑区域的血氧水平依赖信号（BOLD）t吨，右侧39个节点-时间t吨+ 1. 从左到右节点的每一条有向边都表示在患者组（ADHD或对照组）的所有引导样本的公共分量估计中一致出现的时间效应。特别是，我们计算了每个时间效应在公共分量估计中出现的时间比例（在自举估计的总数中），并将结果阈值设置为0.75，仅显示出最一致的效应。

保存图片、插图等的外部文件。对象名为nihms-1530078-f0004.jpg

在单独的窗口中打开

图4：

ADHD患者（上）和对照组（下）的时间效应有向图。当时左脑区域的节点t吨，正确-在时间t吨+ 1. 蓝边区域的自相关效应；红色-两组患者均存在区域间效应；绿色群体的区域间效应。

图4描述了ADHD组和对照组对每个大脑区域（蓝色边缘）都有很强的自相关效应，这是意料之中的。由于它与区域间的时间效应有关，我们有相当多的时间效应存在于两组（红边）和特定组（绿边）中。检测到的时间效应的大小为（平均值±SD）0.23±0.1（自相关效应），0.07±0.05（区域间效应）。从上到下，我们从听觉皮层网络（上两个区域）开始，观察左右听觉皮层对对照组的相互影响，而对ADHD患者则只朝一个方向工作。如中所述[34]，许多研究表明ADHD确诊患者存在听觉问题，这可能反映在他们组的听觉皮层网络缺乏沟通图4研究表明，ADHD患者的默认模式网络（DMN）存在不规则性，这会导致认知能力中断，导致注意力缺失（该病的主要症状之一）。特别是[13]发现了过度激活的模式，而[42]和[33]指出了其失活的缺陷，这降低了对刺激的敏感性。我们的研究结果部分反映了这些想法，显示了多动症组DMN的更多活动：除了自相关效应外，它还包含右侧DMN对内侧前额叶皮层的跨区域时间效应。

两组的腹侧注意网络，右（RVAN，从右背外侧前额叶到右颞下皮质）和左（LVAN，从左顶叶到左极额叶）都显示出大量不同的活动模式。如中所述[13]在儿童ADHD研究中，他们发现VAN中存在激活不足和激活过度的区域。低活性区域在检测和调整环境异常方面表现出ADHD相关缺陷，而高活性区域则强调注意力分散性，这是ADHD最重要的症状之一。同时，在背侧注意网络（DAN，左侧和右侧顶叶内沟）中，我们可以看到更好的控制沟通，这加强了[9]和[40]两者都表明DAN的异常和缺乏交互作用是ADHD患者的特征之一。此外[35]强调与ADHD患者相比，对照组的顶沟内激活增强。至于视觉网络（VN、初级视皮层、左右枕复合体），在比较ADHD组和对照组的[46]，使用[10]揭示ADHD受试者视觉和枕叶皮质的低连通性模式。在我们的研究中，我们发现对照组的活动要多得多，左右枕复合体的时间交叉效应，视觉皮层也会影响这两个区域。

对照组的突起网络（SN、背/腹前扣带回皮层和前岛叶皮层）显示出较高的内源性活动（腹前扣带皮层影响其他两个区域），同时也受到扣带回岛叶网络（CIN、扣带回、右和左岛叶皮质）的严重影响。根据[13]这一网络在决策和处理来自外部因素的信息等执行过程中发挥着至关重要的作用，而这些不足正是ADHD的特点。此外[42]注意缺失与前扣带区刺激前活动减少之间的联系，这证明了我们研究中ADHD组SN的这些区域之间缺乏沟通（同时，对于对照组，也显示出更多的沟通，并且受到CIN网络成员的影响）。对照组的左右颞上沟区域显示出它们之间的时间交叉效应，而ADHD只显示出从左到右的效应，这再次证实了[13]他们声称多动症患者颞区活动减少。最后一个表现出相当活跃的网络是上述扣带回岛网络（CIN）。在这里，ADHD患者表现出扣带皮层区域对右侧岛叶皮层的明显影响，这与ADHD受试者扣带皮层过度激活的假设相一致，如[13]. 同时，控制在右侧和左侧岛叶皮层之间具有时间交叉效应，这两个区域都影响显著性网络的成员。

关于特殊成分估计，这是由我们在第2.2节，这两个研究参与者中的任何一个都没有一致的受试者特定影响（在超过20%的自举样本中没有选择任何一个）。这可能部分归因于我们的数据集缺乏异质性（所有受试者都是十几岁的男孩），以及第3节对于类似的设置(第页= 40，T型= 150，K=10）显示个别成分估计中的假阴性倾向。

5.结论和未来工作

在这项工作中，我们提出了一个联合估计程序，其中多个相关VAR模型是单个基础公共VAR模型的扰动。一个强大的激励性应用是分析大脑功能磁共振成像时间序列数据，用于精神障碍研究，患者共享相同的精神状态（疾病或健康控制）。最终估计值由两阶段算法提供，该算法将时间信号分解为公共分量和单个分量，使用所有受试者通过组套索联合估计公共分量，并通过稀疏套索减去公共VAR信号来估计单个分量。模拟数据的性能表明，在大多数设置中，公共分量的性能是一致的，而单独分量的估计随着联合估计受试者数量的增加而逐渐变得更好。与其他最先进的方法相比，提出的方法在矩阵结构估计方面也显示出优越的性能。此外，由于明确定义了共同成分的概念，我们的两阶段估计方法使得在ADHD患者和对照患者的脑功能磁共振成像数据研究中更容易捕捉和解释群体级效应。

最重要的扩展是为联合估算程序开发一个假设测试框架第2节关于由群套索程序得出的估计值的分布特性，目前还没有足够的文献，这对于测试单个群中时间关系的重要性以及比较不同群之间时间关系的强度都是必要的。关于正则化估计的群级推理问题的论文之一是[31]，其中它被称为群体选择后推断。它概述了一个群体级效应的测试程序，该程序解释了正则化和主体间变异性引入的不确定性。不幸的是，他们没有使用联合估计方法和组套索惩罚，用常规套索分别估计大脑网络，因此不能直接应用于我们的程序。其他需要进一步研究的方面包括：探索单个分量估计的不同初始化，以增加所提出的两阶段迭代过程收敛到全局最小值的机会；评估联合估计程序处理各种滞后阶VAR模型的能力，以及开发滞后阶选择方法。

补充材料

致谢

第一作者的工作得到了研究生奖学金和国家科学基金会的部分支持[授予DMS-1632730]。第二位作者的工作得到了国家科学基金会（CCF 1540093和IIS 1632730）和国家卫生研究院（5 R01 GM114029–03）的部分支持。

6.附录

6.1.MSDL Atlas注释

在表3下面是多学科词典学习（MSDL）地图集中大脑区域的完整学名。枚举对应于区域的显示顺序图1.

脚注

出版商免责声明：这是一份未经编辑的手稿的PDF文件，已被接受出版。作为对客户的服务，我们正在提供这份早期版本的手稿。手稿在以最终可引用的形式出版之前，将经过编辑、排版和校对结果证明。请注意，在制作过程中可能会发现可能影响内容的错误，适用于该期刊的所有法律免责声明均适用。

补充资料

文件补充.pdf包含优化任务（9）的二凸性和额外模拟结果的详细信息。

工具书类

[1]Ban bura M、Giannone D和Reichlin L（2010年）。大型贝叶斯向量自回归.应用计量经济学杂志 25(1), 71–92.[谷歌学者]

[2]Basu S、Michailidis G等（2015）。稀疏高维时间序列模型的正则化估计.统计年鉴 43(4), 1535–1567.[谷歌学者]

[3]Basu S、Shojaie A和Michailidis G（2015年）。具有内在分组结构的网络格兰杰因果关系.机器学习研究杂志 16(1), 417–453.[PMC免费文章][公共医学][谷歌学者]

[4]Beckmann CF、Jenkinson M和Smith SM（2003年）。fmri群体分析的通用多级线性模型.神经成像 20(2), 1052–1063. [公共医学][谷歌学者]

[5]Belilovsky E、Varoqueux G和Blaschko MB（2016年）。高斯图形模型的差异测试：大脑连通性的应用神经信息处理系统研究进展第595-603页。[谷歌学者]

[6]Boyd S和Vandenberghe L（2004年）。凸优化剑桥大学出版社。[谷歌学者]

[7]Breheny P和Huang J（2009）。二级变量选择的惩罚方法.统计及其界面 2(三), 369.[PMC免费文章][公共医学][谷歌学者]

[8]Bühlmann P和Künsch HR（1999）。时间序列引导中的块长度选择.计算统计与数据分析 31(三), 295–310.[谷歌学者]

[9]Castellanos FX、Margulies DS、Kelly C、Uddin LQ、Ghaffari M、Kirsch A、Shaw D、Shehzad Z、Di Martino A、Biswal B等人（2008年）。扣带回相互作用：成人注意缺陷/多动障碍的一个新的功能障碍位点.生物精神病学 63(三), 332–337.[PMC免费文章][公共医学][谷歌学者]

[10]Castellanos FX和Proal E（2012年）。多动症的大尺度脑系统：超越前额叶-纹状体模型.认知科学趋势 16(1), 17–26.[PMC免费文章][公共医学][谷歌学者]

[11]Chu S-H、Lenglet C和Parhi KK（2015）。基于扩散和功能磁共振成像数据的联合脑连通性估计SPIE医学成像第941321-941321页。国际光学与光子学学会。[谷歌学者]

[12]Cole DM、Smith SM和Beckmann CF（2010年）。静止状态fmri数据分析和解释的进展和缺陷.系统神经科学前沿 4.[PMC免费文章][公共医学][谷歌学者]

[13]Cortese S、Kelly C、Chabernaud C、Proal E、Di Martino A、Milham MP和Castellanos FX（2012年）。adhd的系统神经科学：55项fmri研究的荟萃分析.美国精神病学杂志 169(10), 1038–1055.[PMC免费文章][公共医学][谷歌学者]

[14]Danaher P、Wang P和Witten DM（2014年）。多类协方差逆估计的联合图形套索.英国皇家统计学会期刊：B辑（统计方法论） 76(2), 373–397.[PMC免费文章][公共医学][谷歌学者]

[15]Friston KJ、Harrison L和Penny W（2003年）。动态因果建模.神经成像 19(4), 1273–1302. [公共医学][谷歌学者]

[16]Goebel R、Roebroeck A、Kim D-S和Formisano E（2003年）。利用向量自回归模型和格兰杰因果映射研究时间分辨fmri数据中的定向皮层相互作用.磁共振成像 21(10), 1251–1261. [公共医学][谷歌学者]

[17]Goh K、Turan L、Safonov M、Papavassilopulos G和Ly J（1994年）。Biaffine矩阵不等式性质及计算方法《1994年美国控制会议记录-ACC'94》，第1卷，第850-855页。电气与电子工程师协会。[谷歌学者]

[18]Gorski J、Pfeuffer F和Klamroth K（2007年）。双凸集与双凸函数优化：综述与推广.运筹学的数学方法 66(三), 373–407.[谷歌学者]

[19]郭J、列维娜E、米查利迪斯G和朱J（2011）。多图形模型的联合估计.生物特征，asq060。[PMC免费文章][公共医学][谷歌学者]

[20]Härdle W、Horowitz J和Kreiss J-P（2003）。时间序列的Bootstrap方法.国际统计评论 71(2), 435–459.[谷歌学者]

[21]黄J、马S、谢H和张C-H（2009）。一种用于变量选择的群桥方法.生物特征 96(2), 339–355.[PMC免费文章][公共医学][谷歌学者]

[22]Huang S，Li J，Sun L，Ye J，Fleisher A，Wu T，Chen K，Reiman E，Initiative ADN等（2010）。稀疏逆协方差估计学习阿尔茨海默病的脑连接性.神经影像 50(三), 935–949.[PMC免费文章][公共医学][谷歌学者]

[23]Lin J、Basu S、Banerjee M和Michailidis G（2016）。多层高斯图形模型的惩罚最大似然估计.机器学习研究杂志 17(1), 5097–5147.[谷歌学者]

[24]Lin J和Michailidis G（2017）。高维多块向量自回归模型的正则化估计与检验.机器学习研究杂志 18(1), 4188–4236.[谷歌学者]

[25]Liu A、Chen X、Wang ZJ和McKeown MJ（2014）。基于fmri信号的时变脑连接性建模《声学、语音和信号处理（ICASSP）》，2014年IEEE国际会议，第2089-2093页。电气与电子工程师协会。[谷歌学者]

[26]Lütkepohl H（2005）。多时间序列分析新介绍Springer科学与商业媒体。[谷歌学者]

[27]Ma J和Michailidis G（2016）。多图形模型的联合结构估计.机器学习研究杂志 17(1), 5777–5824.[谷歌学者]

[28]Meinshausen N和Bühlmann P（2010）。稳定性选择.英国皇家统计学会期刊：B辑（统计方法论） 72(4), 417–473.[谷歌学者]

[29]Michailidis G和d'AlchéBuc F（2013）。基因调控网络推断的自回归模型：稀疏性、稳定性和因果关系问题.数学生物科学 246(2), 326–334. [公共医学][谷歌学者]

[30]Narayan M和Allen GI（2016）。重采样网络统计的混合效应模型提高了统计能力，以发现多主题功能连接的差异.神经科学前沿 10, 108.[PMC免费文章][公共医学][谷歌学者]

[31]Narayan M、Allen GI和Tomson S（2015年）。图形模型总体的两样本推理及其在函数连通性中的应用.arXiv预打印arXiv:1502.03853.[谷歌学者]

[32]Ramsey JD、Hanson SJ、Hanson C、Halchenko YO、Poldrack RA和Glymour C（2010年）。fmri因果推理的六个问题.神经影像 49(2), 1545–1558. [公共医学][谷歌学者]

[33]Sato JR、Hoexter MQ、Castellanos XF和Rohde LA（2012年）。成人多动症患者的异常脑连接模式：一项一致性研究.公共科学图书馆 7(9)，e45671。[PMC免费文章][公共医学][谷歌学者]

[34]Serrallach B、GroßC、Bernhofs V、Engelmann D、Benner J、Gündert N、Blatow M、Wengenroth M、Seitz A、Brunner M等（2016）。儿童阅读障碍、注意力缺陷多动症和听觉皮层中的神经生物标记物.神经科学前沿 10.[PMC免费文章][公共医学][谷歌学者]

[35]Sigi Hale T、Bookheimer S、McGough JJ、Phillips JM和McCracken JT（2007年）。成人adhd在简单和复杂水平处理过程中的非典型脑激活：一项fmri研究.注意力障碍杂志 11(2), 125–139. [公共医学][谷歌学者]

[36]Simon N、Friedman J、Hastie T和Tibshirani R（2013）。稀疏群套索.计算与图形统计杂志 22(2), 231–245.[谷歌学者]

[37]Song S、Zhan Z、Long Z、Zhang J和Yao L（2011）。基于体素选择的支持向量机方法在fmri数据目标分类中的比较研究.公共科学图书馆 6(2)，e17191。[PMC免费文章][公共医学][谷歌学者]

[38]Tibshirani R（1996）。通过套索回归收缩和选择.英国皇家统计学会杂志。B系列（方法学）, 267–288.[谷歌学者]

[39]Tibshirani RJ、Taylor JE、Candes EJ和Hastie T（2011）。广义套索的解路径斯坦福大学。[谷歌学者]

[40]Tomasi D和Volkow ND（2012年）。注意缺陷/多动障碍儿童的异常功能连接.生物精神病学 71(5), 443–450.[PMC免费文章][公共医学][谷歌学者]

[41]Varoquaux G和Craddock RC（2013年）。学习并比较不同科目的功能连接体.神经影像 80, 405–415. [公共医学][谷歌学者]

[42]Weissman DH、Roberts K、Visscher K和Woldorff M（2006年）。瞬时注意力缺失的神经基础.自然神经科学 9(7), 971. [公共医学][谷歌学者]

[43]Wendell RE和Hurter AP Jr（1976年）。不相交约束下不可分离目标函数的最小化.运筹学 24(4), 643–657.[谷歌学者]

[44]Woolrich MW、Behrens TE、Beckmann CF、Jenkinson M和Smith SM（2004年）。基于贝叶斯推理的fmri群组分析多级线性模型.神经成像 21(4), 1732–1747. [公共医学][谷歌学者]

[45]袁明、林毅（2006）。分组变量回归中的模型选择与估计.英国皇家统计学会期刊：B辑（统计方法论） 68(1), 49–67.[谷歌学者]

[46]朱C-Z、臧Y-F、曹Q-J、闫C-G、何Y、蒋T-Z、隋M-Q和王Y-F（2008）。注意缺陷/多动障碍静止状态脑功能的Fisher判别分析.神经成像 40(1), 110–120. [公共医学][谷歌学者]