Parameter Estimation and Variable Selection for Big Systems of Linear Ordinary Differential Equations: A Matrix-Based Approach

Leqin Wu; Xing Qiu; Ya-xiang Yuan; Hulin Wu

doi:10.1080/01621459.2017.1423074

美国统计协会。作者手稿；PMC 2021 8月11日提供。

以最终编辑形式发布为：

美国统计协会2019年；114(526): 657–667.

2018年7月11日在线发布。数字对象标识：10.1080/01621459.2017.1423074

预防性维修识别码：PMC8357247号

NIHMSID公司：美国国家卫生研究院1502423

PMID：34385718

线性常微分方程大系统的参数估计和变量选择：基于矩阵的方法

吴乐琴,^# 邢秋,^# 雅香园、和胡林武

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

关联数据

补充资料: 补充1。
NIHMS1502423-补充-Supp1.pdf（30万）
GUID:59BADD17-8048-4CAC-824A-555CE8EE2AB5
补充2。
NIHMS1502423-补充-Supp2.pdf（24K）
GUID:B99EA3B7-0A05-48F7-AFC9-B9E86A5A713E
补充3。
NIHMS1502423-补充-Supp3.pdf（1300万）
GUID:C75A3BB6-731C-4D14-959F-C010EEB52BAD
补充4。
NIHMS1502423-补充-Supp4.pdf（1300万）
GUID:7407DEC4-2E41-4B76-8B8D-718A5ED21315
补充5。
NIHMS1502423-补充-Supp5.pdf（14K）
GUID:92F2DF95-57CE-4A31-B7BD-00D4226C69A3
补充6。
NIHMS1502423-补充-Supp6.pdf（1.9万）
GUID:DE8BD671-D882-42D9-B331-62510E667C80
补充7。
NIHMS1502423-补充-Supp7.pdf（24K）
GUID:5E77477A-A2F5-4B74-88AB-B3D9CDE7B34A
补充8。
NIHMS1502423-补充-支持.pdf（18K）
GUID:D6606A91-9168-4267-888A-45A530A1B02B
补充9。
NIHMS1502423-补充-支持.pdf（1100）
GUID:B7D12C31-8490-4B8C-9A03-589584503D40
补充10。
NIHMS1502423-补充-支持10.pdf（64K）
GUID:B2C1A8C1-D0FD-433B-AE48-1C74D473D792
补充11。
NIHMS150223-补充-Supp11.pdf（9.3万）
GUID:CA791A9B-6146-434C-8B23-B6C16C64AAE5
补充12。
NIHMS1502423-补充-补充12.xlsx（1300万）
GUID:965FBFDF-DFDC-4BF5-94EA-B59B02BCAA73
补充13。
NIHMS150223-补充-Supp13.xlsx（450万）
GUID:1599E699-9ACC-41FA-87EA-70A55E8EDEDE

摘要

常微分方程（ODE）广泛用于模拟复杂系统的动态行为。具有线性常微分方程的“大系统”的参数估计和变量选择非常具有挑战性，因为需要在超高维参数空间中进行非线性优化。本文基于相似变换和可分离最小二乘（SLS）的思想，提出了一种参数估计和变量选择方法。仿真研究表明，所提出的基于矩阵的SLS方法可以更准确地估计系数矩阵，并对具有数千维和数百万个参数的线性常微分方程系统进行变量选择，这比直接最小二乘法（LS）要好得多方法和当前可用的基于向量的两阶段方法。我们将这种新方法应用于两个实际数据集：一个是30维930未知参数的酵母细胞周期基因表达数据集，另一个是1250维1563750未知参数的标准普尔1500指数股价数据集，以说明所提出的参数估计和变量选择方法在实际大系统中的实用性和数值性能。

关键词：复杂系统，常微分方程，基于矩阵的变量选择，高维，特征值更新算法，可分离最小二乘

1 简介

常微分方程（ODE）广泛用于模拟复杂系统的动态行为(屠夫，2014;Commumes等人，2011年;德容，2002;Hemker，1972年;Holter等人，2001年;Huang等人，2006年;Lavielle等人，2011年;Li等人，2011年;Lu等人，2011年;拉姆齐等人，2007年). 在许多实际应用中，表征系统的参数通常必须根据数据进行估计。常微分方程的参数估计，也称为反问题，已经用最小二乘法进行了研究(Li等人，2005年;薛等，2010)，可能性(Commumes等人，2011年;Lavielle等人，2011年)和贝叶斯(Putter等人，2002年;Huang和Wu，2006年;Huang等人，2006年,2010)方法。其他几种方法，如主微分分析和广义剖面法(Ramsay等人，2007年;Poyton等人，2006年;拉姆齐，1996年;拉姆齐和西尔弗曼，1998年)和两阶段方法(Hemker，1972年;瓦拉，1982年;陈和武，2008a,b条;梁和武，2008)也提出了建议。

例如，ODE是量化动态基因调控网络（DGRN）的流行模型之一(Bonneau等人，2006年;Li等人，2011年;德容，2002;坂本和伊巴，2001年;杨等人，2002年;Voit，2000年;Holter等人，2001年;Spieth等人，2006年). 基于微阵列的高维时程基因表达数据(Schena等人，1995年;洛克哈特等人，1996年)和RNA-seq(Wang等人，2009年;Garber等人，2011年). 然而，由于计算成本高和模型可识别性问题，上述大多数参数估计方法只适用于最多包含几十个变量的小规模系统(德容，2002;坂本和伊巴，2001年;杨等人，2002年;Voit，2000年;Holter等人，2001年;Spieth等人，2006年). 最近，Lu等人（2011）开发了一种基于线性齐次ODE系统重建DGRN的程序。在这种方法中，差异表达基因（DEG）首先被聚集成共表达模块(Luan和Li，2003年;Ma等人，2006年)基于其表达式的时间模式，以降低维数和缓解可识别性问题。一般来说，对于d日-维线性ODE系统，有第页=d日²+d日需要估计的参数。在这里d日²是ODE系数矩阵中未知参数的总数d日是也需要估计的状态变量的初始条件的数量。例如，即使在降维后，基于ODE的DGRN在Lu等人（2011）仍然包含d日=41个维度（共表达模块）和1722个未知参数，需要从离散的、有噪声的时间进程基因表达数据中估计。我们想指出的一个重要事实是，线性常微分方程组的解是矩阵指数函数(屠夫，2014)高度重视非线性。如果我们直接使用标准非线性最小二乘法（NLS）(薛等，2010)为了估计线性常微分方程系统中的参数，我们需要计算矩阵指数，以评估观测变量与其基于模型的相应预测（或估计）之间的差异。然而，已知矩阵指数在数值上不稳定，无法有效计算(Moler和Van Loan，2003年). 另一种方法是，我们可以使用Runge-Kutta算法等方法对ODE系统进行反复数值求解，以评估NLS目标函数。然而，这样一个高维NLS问题不仅从计算角度很难解决，而且容易陷入局部最优，这可能与真正的全局解相去甚远。

基于上述考虑，Lu等人（2011）应用了两阶段方法(陈和武，2008a，b；梁和武，2008)将ODE系数矩阵解耦为d日数量d日-维向量；然后他们应用了SCAD方法(范和李，2001)分别针对每个行向量（方程）同时进行参数估计和变量选择。这种方法实现简单，计算效率高。然而，这种基于向量的变量选择方法忽略了ODE系数矩阵中固有的丰富结构信息，严重依赖于对测量误差敏感的状态变量导数的良好估计。因此，它经常导致参数估计不准确和变量选择结果不佳(丁和吴，2014)。

在本文中，我们提出了一种新的基于矩阵的方法，以避免基于向量的两阶段方法的估计不良和NLS方法的计算问题(Xue等人，2010年). 该方法的核心是一种特殊形式的可分离最小二乘（SLS）方法(Ruhe和Wedin，1980年)基于系数矩阵的Jordan标准分解（JCD），它本质上将原来的非线性优化问题转化为一个等价问题，其中只有d日特征值的数量，而不是全部d日²+d日参数，需要通过非线性优化算法进行估计。其余的可以通过一个封闭式公式得到，计算量很小。我们进一步利用了JCD中使用的相似变换后的线性ODE系统解的解析形式，以避免在评估NLS目标函数时对原始ODE系统进行数值求解。此外，导出的目标函数解析形式具有解析梯度，可以稳定有效地计算。原始未知参数的估计是从系数矩阵特征值估计的封闭函数中恢复出来的。在仿真研究中，我们表明，新方法不仅速度快，但与其他方法相比，更频繁地到达全局最优值，并产生更准确和稳定的估计。最后，我们将提出的方法应用于两个实际应用，一个是DGRN建模d日=30个尺寸和第页=930个未知参数，另一个是股票市场系统建模d日=1250尺寸和第页=1563750未知参数，以证明使用该方法可以很好地恢复大型线性ODE系统。

2 模型和方法

2.1、。型号说明

我们考虑以下高维齐次线性常微分方程系统的参数估计问题

{\begin{array}{l} \frac{d日 x个 (t吨)}{d日 t吨} = 一个 x个 (t吨), \forall t吨 \in [{T型}_{1}, {T型}_{2}], \\ x个 ({T型}_{1}) = {x个}_{0}, \end{array}

(1)

哪里

x个 (t吨) = {({x个}_{1} (t吨), {x个}_{2} (t吨), \dots, {x个}_{d日} (t吨))}^{T型}

(2)

是一个d日-满足范围的维状态变量向量

0 \leq {T型}_{1} < {T型}_{2} < \infty .

系数矩阵 $一个 \in ℝ^{d日 \times d日}$ 和初始条件 ${x个}_{0} \in ℝ^{d日}$ 是系统中需要使用观测数据估计的未知参数。

在实际应用中，我们假设x个(t吨)在有限时间点用独立误差测量 $({t吨}_{1}, {t吨}_{2}, \dots, {t吨}_{n个})$ ，并且每个时间点的测量误差遵循非奇异协方差矩阵的高斯分布 $Σ_{ϵ}$ 即。，

年_{我} ({t吨}_{j个}) = {x个}_{我} ({t吨}_{j个}) + ϵ_{我 j个}, ϵ_{\cdot j个} \in ℝ^{d日} \sim N个 (0, Σ_{ϵ}), 我 = 1, \dots, d日; j个 = 1, \dots, n个 .

(3)

为了方便起见，我们表示d日×n个-维数据矩阵 ${年_{我} ({t吨}_{j个})}$ 统称为年.

基于最大似然原理，估计的逆问题一个和x个₀可以表示为以下非线性加权最小二乘（NWLS）最小化问题

\underset{一个, {x个}_{0}}{最小值} {‖ 年 (t吨) - x个 (一个, {x个}_{0}) (t吨) ‖}_{Σ_{ϵ}}^{2} .

（4）

根据上下文， $年_{我} (t吨)$ 可以表示观测曲线或离散数据t吨; $x个 (一个, {x个}_{0}) (t吨) : = {e（电子）}^{t吨一个} {x个}_{0}$ 是ODE系统的解曲线(1)带参数一个和x个₀.上述优化空间的维数为第页=d日²+d日.两种典型的规范选择 $∥ \cdot ∥_{Σ_{ϵ}}$ 在中使用方程式（4）离散观测的加权欧几里德度量和加权L（左）²-函数的度量，在补充文本(截面S2)。

2.2. 相似变换与可分离最小二乘

优化问题（4）原则上，可以通过为非线性最小二乘问题设计的任何合适的非线性优化算法直接进行数值求解，在本研究中称为“直接LS方法”。实际上，当参数空间较大时(例如,d日≥100），非线性优化问题的维数可能非常大( $第页 = {d日}^{2} + d日 ⩾ 10, 100$ 如果d日⩾100），这很难用数值求解，并且很可能陷入局部解。在本小节中，我们提出了一种基于相似变换（ST）和可分离最小二乘（SLS）的方法，旨在显著降低非线性优化维数，从而扩展我们处理大型ODE系统的能力。

对于ODE系统(1)，我们进一步假设系数矩阵一个谱中没有多重性（没有两个特征值完全相同）。这个假设不会导致太大的通用性损失，因为这样的矩阵只形成一个零测度集（w.r.t.或dμ，标准勒贝格测量 $ℝ^{d日 \times d日}$ 或任何绝对连续w.r.t.d的概率测度μ，如与实际随机矩阵相关的概率测度，如Ginibre系综、高斯正交系综、Wishart系综等）d日-由-d日矩阵(Ginibre，1965年;莱曼和索默斯，1991年;陶，2012)。

备注。

我们想指出，如果一个谱中确实存在多重性，那么就存在其他可以生成完全相同曲线（或数据点）的系数矩阵。换句话说，这个系统在理论上是不可识别的。如果没有其他结构信息一个已给出先验的，我们将无法恢复一个即使我们得到无限多没有噪声的观测结果。

在这个假设下一个是

一个 = 问 Λ 问^{- 1},

(5)

哪里Λ是一个只有两种类型块的块对角矩阵：一个包含一个实特征值的1×1块一个; 或2×2块 $[\begin{matrix} 一 & b条 \\ - b条 & 一 \end{matrix}]$ 对应于一对共轭复特征值一±双.我们可以随时选择适当的安排问使得对角块Λ组织如下

Λ = [\begin{matrix} 一_{1} & {b条}_{1} \\ - {b条}_{1} & 一_{1} \\ ⋱ \\ 一_{k个} & {b条}_{k个} \\ - {b条}_{k个} & 一_{k个} \\ {c（c）}_{1} \\ ⋱ \\ {c（c）}_{d日 - 2 k个} \end{matrix}],

(6)

哪里 $k个 \leq \frac{n个}{2}$ 是一个非负整数，表示Λ.

定理2.1。

优化问题（4）等于

\underset{Λ, 问}{最小值} {‖ 年 (t吨) - 问^{T型} x个 (Λ, e（电子）) (t吨) ‖}_{Σ_{ϵ}}^{2},

(7)

哪里e（电子）是由k对(0,1)^T型 的s和d− 21的k如下

e（电子） = (\overset{k个 属于 (0, 1) 秒}{\overset{︷}{0, 1, \dots, 0, 1}}, \overset{d日 - 2 k个 属于 1^{'} 秒}{\overset{︷}{1, \dots, 1}})^{T型} .

(8)

证明。请参阅补充文本(第S6.1节)。

备注。

值得注意的是定理2.1转换原始优化问题（4）具有d日²+d日优化问题中的未知参数(7)，它还有一个总数d日²+d日未知参数（即。，d日²中的未知参数问和d日特征值Λ). 我们当前的选择e（电子）在里面等式（8）确保新转换的ODE系统有一个简单的解决方案（见下文），这使我们能够为问（请参见等式（13）)何时Λ给出了。然后我们可以应用可分离LS方法来估计Λ，这只是一个非线性优化d日未知参数。之后Λ是估计的，问（带有d日²未知参数）可以通过闭合解计算等式（13）只需很少的计算成本。

众所周知，对于给定的块二对角矩阵ΛODE系统

{\begin{array}{l} \frac{d日 x个 (t吨)}{d日 t吨} = Λ x个 (t吨), \forall t吨 \in [{T型}_{1}, {T型}_{2}], \\ x个 ({T型}_{1}) = e（电子）, \end{array}

(9)

有以下解决方案

{\begin{array}{l} {x个}_{2 j个 - 1} (Λ, e（电子）) (t吨) = 经验 (一_{j个} t吨) 罪 {b条}_{j个} t吨, \\ {x个}_{2 j个} (Λ, e（电子）) (t吨) = 经验 (一_{j个} t吨) 余弦 {b条}_{j个} t吨, \end{array} j个 = 1, 2, \dots, k个

(10)

{x个}_{j个} (Λ, e（电子）) (t吨) = 经验 ({c（c）}_{j个 - 2 k个} t吨), j个 = 2 k个 + 1, 2 k个 + 2, \dots, d日 .

(11)

在这里 ${x个}_{我} (Λ, e（电子）) (t吨)$ 是我解向量的第个分量 $x个 (Λ, e（电子）) (t吨)$ 、和e（电子）是给定的常量向量(8).

请注意，对于固定Λ，优化中的目标函数(7)关于问可以归结为一个具有封闭解的线性回归问题，即对于给定的Λ

\underset{问}{最小值} {‖ 年 (t吨) - 问^{T型} x个 (Λ, e（电子）) (t吨) ‖}_{Σ_{ϵ}}^{2}

（12）

给出了一个封闭的解决方案（请参见补充文本第S6.2节用于扣除）。

问 (Λ) = {〈 x个 (Λ, e（电子）) (t吨), x个 (Λ, e（电子）) (t吨) 〉}^{- 1} \cdot 〈 x个 (Λ, e（电子）) (t吨), 年 (t吨) 〉 .

(13)

使用封闭式解决方案(10),(11)、和(13)，优化问题(7)可以转化为只涉及Λ根据可分离LS原理

\underset{Λ}{最小值} {‖ 年 (t吨) - 问^{T型} (Λ) x个 (Λ, e（电子）) (t吨) ‖}_{Σ_{ϵ}}^{2} .

(14)

尽管事实上Λ是一个d日×d日矩阵，它只包含d日由于其特殊的结构，需要估计的未知参数。此外，目标函数(14)相对于中的参数是连续可微的Λ因此，我们可以推导出参数梯度的解析公式，这将大大加快非线性优化过程（参见补充文本,第S6.4节)。

一旦我们获得 $\hat{Λ}$ ，估计Λ使其最小化(14)，我们可以立即获得 $\hat{问}$ ，矩阵的回归估计问，来自方程式（13）原始ODE系数矩阵和初始条件的估计值可以计算为

\hat{一个} = \hat{问} \hat{Λ} {\hat{问}}^{- 1}, {\hat{x个}}_{0} = \hat{问} e（电子）,

(15)

哪里e（电子）在中给出(8).

以下定理表明 $\hat{一个}$ , ${\hat{x个}}_{0}$ 确实是（4）.

定理2.2。

$\hat{Λ}$ 是的最小值(14)当且仅当 $(\hat{一个}, {\hat{x个}}_{0})$ 由生成方程式（15）是的最小值（4）.

证明。请参阅补充文本(第S6.3节)。

基于原最小二乘问题局部极小值之间的双射（4）以及重新制定的问题(14)再加上这两个目标函数在相应的局部极小值处具有相同的值，我们立即得到了以下推论。

推论2.3。

$\hat{Λ}$ 是全球最小值(14)当且仅当 $(\hat{一个}, {\hat{x个}}_{0})$ 由生成方程式（15）是全球最小值（4）.

的本质定理2.2和推论2.3就是原来的NLS优化问题（4），尺寸为d日²+d日，等价于特征值估计问题，这是一个维数的非线性优化问题d日这是非线性优化问题的显著降维。我们需要不同数据点的数量n>d为了避免可辨识性问题，虽然未知参数的总数第页=d日²+d日可以是大得多比n个我们在补充文本(S1段). 当我们选择Levenberg-Marquardt算法（LMA）来解决重新表述的优化问题时(14)在ST-SLS算法中，基于其灵活性，原则上可以用任何合适的优化算法来代替。

我们想指出，尽管上述方法是为齐次线性ODE模型开发的(1)，它可以应用于异构线性ODE模型，只需在状态变量中添加一个附加常数项的简单数学技术x个(t吨). 详细讨论见补充文本(截面S4)。

2.3. 渐近方差与推断

在本节中，我们提供了以下项的渐近方差-方差矩阵估计 $\hat{一个}$ ，表示参数估计的不确定性。本节中两个定理的证明在补充文本(第S6.5节)。

我们考虑以下内容一个作为vec(一个)，这是一个d日²-参数的维向量，以便 $向量 {(一个)}_{(我 - 1) d日 + k个} : = {一个}_{我 k个}$ 。我们定义D类(一个,t吨)如下所示d日×d日²-维数矩阵函数

{D类}_{我, (我 - 1) d日 + k个} (一个, t吨) = \frac{\partial {({e（电子）}^{t吨 一个} {x个}_{0})}_{我}}{\partial {一个}_{k个 我}}, 我 = 1, \dots, d日; 我 = 1, \dots, d日; k个 = 1, \dots, d日 .

(16)

显然地， ${D类}_{我, (我 - 1) d日 + k个} (一个, t吨)$ 是解曲线的雅可比矩阵， $x个 (t吨) : = {e（电子）}^{t吨一个} {x个}_{0}$ ，关于一个在时间进行评估t吨然后，我们可以表示一个作为的函数D类(一个,t吨)。

定理2.4。

鉴于 x个 ₀ 和 $Σ_{ϵ}$ （测量误差的协方差矩阵），总Fisher信息矩阵属于ODE系统的估计 (1) 是

我 (一个) = \sum_{j个 = 1}^{n个} {D类}^{T型} (一个, {t吨}_{j个}) Σ_{ϵ}^{- 1} D类 (一个, {t吨}_{j个}) .

(17)

显然， ${D类}^{T型} (一个, {t吨}_{j个}) Σ_{ϵ}^{- 1} D类 (一个, {t吨}_{j个})$ 对所有人来说都是正定的j个。因为我(一个)是正semi-definite矩阵的总和，它必须为正semi-definite。因此，我(一个)只要它是满秩的，它就是正定的。

基于定理2.4，我们有以下渐近结果 $\hat{一个}$ .

定理2.5。

假设

$\hat{一个}$ 是唯一的全球最小值(14).
我(一个),总Fisher信息矩阵是全秩的（因此是正定的）。
ODE系统(1),考虑到这一点一个是真正的系统矩阵x个₀是真实的初始条件，在以下意义上是可识别的。如果有矩阵B类使得 ${e（电子）}^{{t吨}_{j个} B类} {x个}_{0} = {e（电子）}^{{t吨}_{j个} 一个} {x个}_{0}$ 为所有人 $j个 = 1, 2, \dots, n个$ ,然后B类=一个.

根据这些假设，当 $n个 \to \infty, \hat{一个}$ 在分布上收敛到具有正确平均值的正态分布(一个)和协方差矩阵 我(一个)⁻¹

\hat{一个} \sim 一个 N个 (一个, 我 {(一个)}^{- 1}) .

(18)

根据定义 ${\hat{一个}}_{我 j个}$ 是的对角线元素我(一个)⁻¹更具体地说， $无功功率，无功功率 ({\hat{一个}}_{我 j个}) = 我 {(一个)}_{(j个 - 1) d日 + 我, (j个 - 1) d日 + 我}^{- 1}$ 渐进地。通过这些方差估计，我们可以检验零假设 ${H（H）}_{0, 我 j个} : {一个}_{我 j个} = 0$ 与相应的替代假设相反 ${H（H）}_{1, 我 j个} : {一个}_{我 j个} \neq 0$ 通过使用标准化网络强度， ${z（z）}_{我 j个} : = \frac{{一个}_{我 j个}}{\sqrt{无功功率，无功功率 ({\hat{一个}}_{我 j个})}}$ ，作为测试统计。这种统计遵循渐近标准正态分布 ${H（H）}_{0, 我 j个}$ 。因为我们需要测试大量数据(d日²)假设，一个合适的多重测试程序，如Holm-Bonferroni程序霍尔姆（1979），Šidák过程(什伊达克，1967年)Benjamini-Hochberg程序(本杰米尼和霍奇伯格，1995年)需要应用于控制总体类型I错误。参数估计的置信区间也可以基于渐近结果来构造。

2.4. 变量选择

对于以下问题先验的系数矩阵的信息一个是一个稀疏矩阵，在系数矩阵估计中添加一个正则项将稀疏性强加于系数矩阵估计是有利的。

\underset{一个, {x个}_{0}}{最小值} {‖ 年 (t吨) - x个 (一个, {x个}_{0}) (t吨) ‖}_{Σ_{ϵ}}^{2} + ρ (一个) .

(19)

刑罚期限的可能选择ρ(一个)包括LASSO(Tibshirani，2011年)、SCAD(范和李，2001)、MCP(张，2010)等。

按照我们在前面小节中所做的相似性转换，将导致

\underset{Λ}{最小值} {‖ 年 (t吨) - 问^{T型} (Λ) x个 (Λ, e（电子）) (t吨) ‖}_{Σ_{ϵ}}^{2} + ρ (问 (Λ) Λ 问^{- 1} (Λ)),

(20)

哪里问(Λ)由定义(13)分别是。我们可以使用相同的优化算法来解决上述问题。

请注意，将中的目标函数最小化(20)应导致矩阵的稀疏估计一个理论上。然而，我们仍然需要使用可分离LS估计(13)对于问，这是（12），而不是(20)。此近似可能不会缩小对一个精确到零。因此，我们需要确定一个数字阈值c（c）这样，如果 $| {\hat{一个}}_{我 j个} | < c（c）$ ，我们替换 ${\hat{一个}}_{我 j个}$ 乘以零。类似的想法也被用于消除由于以下数值误差导致的小非零值估计L（左）₁正则化回归算法，如LASSO(Yukawa等人，2012年;Combettes and Wajs，2005年). 确定阈值的一个简单方法是在(17)和(18)在里面第2.3节表示渐近z（z）-测试以检查是否 ${一个}_{我 j个} \neq 0$ 然而，此方法不适用于大型网络或系统，因为它需要估计d日²×d日²-维协方差矩阵，当d日另一种方法是根据标准分类方法，如K近邻（KNN）算法，选择阈值将估计系数分为两组：零组和非零组。

注意，目标函数参数的闭合梯度(20)不可用。因此，NEWUOA等无导数优化（DFO）算法(鲍威尔，2006年;Zhang等人，2010年)在这种情况下，需要使用，并且与ST-SLS算法相比，计算成本更高（请参阅补充文本,S3剖面). 带有变量选择的ST-SLS伪码（ST-SLS-VS）在补充文本(S1段)。

三。模拟研究

3.1. 模拟实验设计

在本节中，我们通过不同维度和不同噪声水平的仿真研究，将基于特征值估计框架的拟议ST-SLS和ST-SLS-VS方法与现有方法进行了比较。请注意，设计高维ODE模拟实验并非易事。为了生成合理的ODE仿真模型，必须特别注意避免仿真系统的共线，并缓解可识别性问题。特别是系数矩阵的特征值一个需要彼此分开。基于这些考虑，我们得出一个首先生成具有良好性质的特征值，然后随机生成其特征向量。更具体地说，我们首先通过其实部和虚部分别生成特征值。为了使系统稳定，每个特征值的实部应该是非正的。但它不能太负，否则ODE解作为时间的指数函数将非常迅速地衰减到零，这可能会导致ODE系统在数值上无法识别(Miao等人，2011年). 在我们的模拟实验中，特征值的实部由[-0.7,0]上的均匀分布生成。特征值的虚部只需要彼此远离就可以了。例如，在我们的模拟研究中，一个典型的选择是±2π,±4π,...,±dπ添加了一个小的高斯噪声。一次Λ生成，我们将其乘以随机生成的非奇异矩阵问创建系数矩阵，即。，

一个 = 问 Λ 问^{- 1} .

(21)

从技术上讲问可以是任何可逆的平方矩阵，但是对于变量选择实验，系数矩阵应该是稀疏的。因此我们使用矩阵问具有特殊的块-对角结构，保证了两者的稀疏性问因此，它的逆函数一个可以生成稀疏矩阵。

一旦生成ODE系数矩阵，就可以从ODE模型中生成观测数据(1)使用其分析溶液。观测的时间点均匀分布在区间[0,1]上。将随机噪声添加到ODE系统的模拟数据中，即身份证号码。高斯噪声分布 $N个 (0, {(α σ)}^{2})$ ，其中σ取原始数据的样本标准偏差，α控制噪声级，分别取0、0.1或0.3；其中0代表无噪情况。仿真模型的尺寸设置为d日分别为301003001000。所有结果均基于1000个随机模拟给出，但1000维情况除外，该情况基于100个模拟，因为计算成本较高。所有仿真都是在一台运行Xubuntu 14.04操作系统的笔记本电脑上进行的，该笔记本电脑具有2.5GHz的CPU和8G的RAM。

3.2. 参数估计比较

在本小节中，我们给出了所提出的ST-SLS方法和仿真研究中的直接LS方法之间的参数估计比较结果。为了进行公平比较，采用Levenberg-Marquardt算法作为ST-SLS方法的优化求解器（迭代更新Λ)以及直接LS方法（估计一个直接最小化LS目标函数）。请参阅补充文本(S3剖面)有关优化算法的更多详细信息。

参数估计比较的模拟结果报告于表1针对不同的ODE系统尺寸，我们比较了这两种方法在计算成本、图像质量和参数估计精度方面的差异(d日)和不同的噪音水平(α). 计算成本由用于运行算法的CPU时间（秒）量化。通过模型拟合的相对残差平方和（RRSS）来评估菲特的良好性，即最终解的目标函数值除以数据矩阵的平方Frobenius范数。通过相对估计误差（REE）来评估总体参数估计值，其定义为

稀土元素 (一个) = \frac{∥ \hat{一个} - 一个 ∥_{F类}}{{‖ {一个}_{0} ‖}_{F类}} \times 100 %, 稀土元素 ({x个}_{0}) = \frac{{‖ {\hat{x个}}_{0} - {x个}_{0} ‖}_{2}}{{‖ {x个}_{0} ‖}_{2}} \times 100 %,

(22)

其中(一个,x个₀)是真正的参数和 $(\hat{一个}, {\hat{x个}}_{0})$ 是他们的相应估计。

表1：

直接最小二乘法和我们新的ST-SLS方法之间的参数估计比较。计算成本（CPU时间，单位为秒）、相对RSS（RRSS）和REE以百分比形式给出，基于大多数情况下1000次模拟运行的平均值（括号中的标准偏差），但具有以下特征的Direct-LS方法除外d日=100，ST-SLS方法d日=1000，由于计算成本高，仅使用了100次模拟运行。

		直接-LS				ST-SLS型
尺寸(第页)	噪音(α)	时间（秒）	RRSS（%）	稀土元素(一个)(%)	稀土元素(x个₀)(%)	时间（秒）	RRSS（%）	稀土元素(一个)(%)	稀土元素(x个₀)(%)
30	0	852	0	0	0	0.005	0	0	0
30	0.1	908(632)	30(12)	610(350)	43(15)	0.013(0.0005)	0.019(0.0018)	0.21(0.050)	0.020(0.0062)
30	0.3	866（349）	41(15)	1220(700)	58(31)	0.014(0.0008)	0.17(0.016)	2.1(0.64)	0.21（0.0056）

100	0	14982	0	0	0	1	0	0	0
100	0.1	16010(3972)	35(13)	760(300)	42(15)	1.2(0.49)	0.017(0.0054)	0.97(0.091)	0.023(0.0034)
100	0.3	55293(41769)	45(31)	1111(812)	73(51)	1.2(0.50)	0.17(0.0098)	1.7(0.99)	0.021(0.0030)

300	0	-	-	-	-	93	0	0	0
300	0.1	-	-	-	-	99(5.0)	0.017（0.0016）	2.3(0.41)	0.023(0.0019)
300	0.3	-	-	-	-	101(2.8)	0.015(0.0014)	4.0(0.74)	0.021(0.0017)

1000	0	-	-	-	-	7613	0	0	0
1000	0.1	-	-	-	-	8437(13)	0.018（0.0020）	8.8(2.9)	0.030(0.0094)
1000	0.3	-	-	-	-	8434(22)	0.015(0.0041)	20.6(16.7)	0.020(0.0069)

在单独的窗口中打开

发件人表1，我们看到，对于无噪声的情况(α=0），这两种方法对所有维度都产生了良好的参数估计，且拟合良好。计算时间随系统维数的增加而增加(d日)如预期。这表明，在无测量误差的理想情况下，这两种方法都是好的。在数据中加入测量噪声时，直接最小二乘法的结果较差。直接最小二乘法的拟合相对误差（RRSS）可高达30–45%，估计误差（REE）可为系数矩阵估计值与真实值之间的6至12倍差异一个这表明直接-LS方法可能收敛到远离真实解的局部解。对于更高维度的情况(d日=300和1000），直接最小二乘法不收敛，无法获得估计值。相反，我们的新ST-SLS方法对所有模拟情况都产生了合理的结果。模型拟合的相对RSS很低(<对于所有仿真案例，该方法都比直接最小二乘法小得多，表明我们的新方法非常适合模型。系数矩阵的REE一个范围从<1%至20.6%，初始值估算的REE甚至更小。因此，ST-SLS算法对所有未知参数都有很好的估计。我们还观察到，对于这两种方法，初始状态的估计 ${\hat{x个}}_{0}$ 比系数矩阵的估计要好得多一个这是因为估计x个₀当时是否对合适的解决方案进行评估t吨=0，并且模型拟合总是非常好。此外，所提出的ST-SLS速度非常快，对于低维或中等维的情况，只需几秒钟就能产生结果(d日=10至100），这需要多个小时的CPU时间用于直接-LS方法。对于高维情况(d日=300和1000），对于直接-LS方法无法获得结果的情况，ST-SLS算法仍能在几分钟内获得良好的结果(d日=300）或几个小时(d日=1000），这证明了我们处理大型系统的新方法的可扩展性。

3.3. 变量选择比较

对于高维ODE变量选择，现有的唯一可行计算方法是两阶段方法(Lu等人，2011年). 在本小节中，我们将我们新的ST-SLS-VS算法（配备了三个不同的正则化项）与两阶段方法在大常微分方程系统的变量选择性能方面进行了比较。

我们用不同的噪音水平对d日分别=30、100和300；和100次模拟运行d日=1000，因为计算成本高。基于这些模拟运行平均值的结果报告于表2在这项模拟研究中，我们比较了两阶段方法之间的灵敏度（SEN）和特异性（SPE），这两种方法分别测量变量选择的真阳性率和假阳性率(Lu等人，2011年)以及提出的ST-SLS-VS方法。

表2：

两阶段方法和我们的新ST-SLS-VS算法之间的ODE变量选择比较。灵敏度（SEN）和特异性（SPE）都是基于1000次案例模拟运行得出的百分比d日=3010300和100次模拟运行d日=1000，括号内为标准偏差。

		2阶段方法		ST-SLS-VS（拉索）
尺寸(d日)	噪音(α)	发送百分比	特殊目的实体%	发送百分比	特殊目的实体%
30	0	91.7	98.8	100	100
30	0.1	93.3(1.1)	63.6(1.3)	100(0)	100(0)
30	0.3	85.6（3.8）	43.5(1.5)	100(0)	100(0)

100	0	94	99.5	100	100
100	0.1	94.7(1.3)	67.3(0.3)	100(0)	100(0)
100	0.3	89.8（1.3）	65.9(0.2)	100(0)	99(0.1)

300	0	96.3	99.7	100	100
300	0.1	93.9（0.1）	62.1(0.2)	100(0)	99(0.1)
300	0.3	88.3(0.5)	62.5(0.2)	91(0.9)	98(0.1)

1000	0	98	99.9	100	100
1000	0.1	90.8(2.7)	60.2(1.5)	85(3.9)	99(0.1)
1000	0.3	85.1(4.3)	53.3(6.2)	80(5.5)	97(0.1)

在单独的窗口中打开

正如我们在中指出的那样第2.4节，没有目标函数的闭式梯度公式(20)可用于实现ST-SLS-VS算法。我们必须使用无导数优化算法，如NEWUOA(鲍威尔，2006年;Zhang等人，2010年)用于优化（请参见补充文本,S3剖面)，这需要较高的计算成本。例如，ST-SLS-VS算法为维度生成了平均约5分钟的结果d日=300箱，尺寸为5-6小时d日=1000个高噪声级的案例，比ST-SLS算法慢。ST-SLS-VS算法较慢的原因有两个：不能使用闭合形式的梯度，并且目标函数的评估更复杂。我们使用Fortran实现了ST-SLS-VS算法，而两阶段方法是在R中实现的，这就是为什么我们没有比较ST-SLS-VS算法和两阶段方法的计算成本。但是，一般来说，两阶段方法要快得多，因为它将线性ODE参数估计转换为线性回归模型拟合。

发件人表2和和3，三我们可以看到，两阶段方法的灵敏度在85%到95%之间，对于大多数模拟情况来说通常是好的，但对于有噪声的数据，其特异性很低（在44%到67%之间）。相比之下，所提出的ST-SLS-VS方法的性能非常稳定，三种正则化项的选择产生了类似的结果。我们的方法不仅在所有无噪声情况下都能准确地识别出正确的结果，而且在其他情况下也具有很好的敏感性（大多高于80%）和特异性（97-100%）。总的来说，我们的新方法在变量选择方面优于现有的两阶段方法。

表3：

两阶段方法和我们的新ST-SLS-VS算法之间的ODE变量选择比较。灵敏度（SEN）和特异性（SPE）都是基于1000次案例模拟运行得出的百分比d日=3010300和100次模拟运行d日=1000，括号内为标准偏差。

		ST-SLS-VS（SCAD）		ST-SLS-VS（MCP）
尺寸(d日)	噪音(α)	发送百分比	特殊目的实体%	发送百分比	特殊目的实体%
30	0	100	100	100	100
30	0.1	100(0)	100(0)	100(0)	100(0)
30	0.3	100(0)	100(0)	100(0)	100（0）

100	0	100	100	100	100
100	0.1	100(0)	100(0)	100(0)	100(0)
100	0.3	100(0)	100(0)	100(0)	100(0)

300	0	100	100	100	100
300	0.1	100(0)	99(0.1)	100(0)	99(0.1)
300	0.3	92(0.9)	98(0.1)	92（0.9）	98(0.1)

1000	0	100	100	100	100
1000	0.1	87(4.1)	98(0.1)	86(4.3)	98(0.1)
1000	0.3	81(5.4)	93(0.8)	82(5.4)	93(0.9)

在单独的窗口中打开

4 实际数据分析

我们将所提出的ODE参数估计和变量选择方法应用于两个应用数据集，以说明它们对于大规模系统的实用性和可扩展性。第一个是从固定相酵母培养物中收集的一组时间进程微阵列数据(Aragon等人，2006年)具有30个维度和930个未知参数的中型系统。第二个是由标准普尔指数化的股票的10年历史日价值组成，该指数具有1250个维度和1563750个未知参数的大系统。

4.1. 时程酵母芯片数据分析

第一个应用示例是时间进程基因表达数据的子序列（基因表达总括数GSE3688标准)从暴露于氧化应激的固定期培养的酵母细胞中收集(Aragon等人，2006年). 这些数据每1分钟收集一次，持续35分钟，最后一个时间点为60分钟（共37个时间点），使用微阵列。我们应用了函数主成分分析方法(Wu和Wu，2013年)并确定了与周期调控相关的前30个重要基因(Spellman等人，1998年). 我们的目标是使用线性ODE模型研究这30个基因之间的调控关系。

我们应用了我们提出的方法和开发的ST-SLS/ST-SLS-VS算法，如第2节并恢复了前30个重要酵母细胞周期相关基因的动态网络。我们得到了估计的动力系统系数矩阵( $\hat{一个}$ )每个边缘的标准偏差使用方程式（17）和(18)在里面第2.3节。如中所讨论和建议第2.4节，我们使用了双面z（z）-使用Holm-Bonferroni多重测试程序进行测试(霍尔姆，1979年)并通过将家庭错误率控制在0.05来确定网络稀疏性。生成的网络具有95%的稀疏性，如所示图1注意，此图中未包括14个分离基因（FLC2、PET9、RDH54、BEM1、BUD3、NDC80、MMR1、CAR2、SPT21、GCV2、WHI3、ARG1、GNT和YKR012C）。重建的基因调控网络在补充表S1.

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

在单独的窗口中打开

图1：

从时间进程酵母微阵列数据重建基因调控网络。不同大小的节点表示网络度，网络度定义为重构网络中节点的相邻边数。正（负）规用绿色（红色）表示。

发件人图1，我们可以看到SST2、PUT1、ZSP1、DSN1和SPC34是具有最大相邻边数（网络度）的中央集线器节点。根据酿酒酵母基因组数据库（SGD）(Cherry等人，1998年)，SST2编码GPA1P的GTPase激活蛋白，这是阻止交配途径的受体依赖性信号传导所必需的。该基因的零突变导致细胞大小增加和生长速度降低。PUT1编码脯氨酸氧化酶，当脯氨酸是唯一的氮源时，该基因的突变导致酵母无法生长。ZSP1是一种功能未知的蛋白质，但已知与磷酸代谢途径的成员基因PHO88相互作用。DSN1是MIND动粒复合体的重要组成部分，已知在纺锤体微管与参与减数分裂姐妹染色单体分离的动粒结合中发挥重要作用。SPC34是一个纺锤极成分，是Dam1复合体（DASH复合体）的一个基本亚单位。DSN1和SPC34都是动粒的组成部分，它们之间的联系已经建立(田中等人，2005年;Pramila等人，2006年). SPC34和SST2之间的连接也已记录在案(Montpetit等人，2005年)。

通过我们的方法确定的其他网络连接是新颖的，可能有助于生成进一步调查的假设。例如，ZSP1是一个研究不足的基因，已知它只与PHO88相互作用。我们发现它与磷酸代谢途径中的另一个成员基因PHO89有很强的联系。这一观察结果表明，ZSP1可能在磷酸盐代谢中发挥比我们目前所知更重要的作用。PUT1和SST2之间的紧密联系有点令人惊讶和有趣，因为PUT1与SST2似乎实现了截然不同的生物功能。PUT1对于酿酒酵母消化脯氨酸，脯氨酸是葡萄中最丰富的氮源，也是野生酵母的自然环境(黄和布兰德里斯，2000). SST2以其调节交配反应的功能而闻名，这似乎与脯氨酸消化无关。然而，已知SST2也参与细胞增殖(Lopez等人，1997年)和生长，尤其是在营养有限的环境中(Lopez等人，2001年;Boer等人，2003年). 我们的研究结果表明，PUT1和SST2在氮代谢和细胞生长之间的相互作用中可能有密切关系。此外，我们发现网络中的三个基因（YLR297W、YMR253C、YPR174C）在文献中没有明确的生物学注释。其中，YLR297W是PUT1和SST2这两个连接最紧密的集线器节点的调节器。我们的结果可能为这些基因生物学功能的未来实验研究提供有用的见解。

4.2。标准普尔股市数据分析

传统的随机微分方程（SDE）模型，如Black-Scholes-Marton模型(布莱克和斯科尔斯，1973年;默顿，1973年;Øksendal，2013年)用于建模股市价格数据。众所周知，相应的ODE模型可以用来描述SDE的平均行为(艾哈迈德，1998年)（第2章中的定理1）。在这里，我们将线性ODE模型应用于标准普尔1500（也称为标准普尔复合1500指数）的股价数据，以研究标准普尔1500中公司股价变化的长期动态交互作用。本研究使用的数据涵盖了这些股票2004年至2014年（2668个交易日）的10年期每日收盘价。原始指数包含1501只股票，其中251只因缺失和其他数据问题而从分析中删除。基于剩余的1250只股票，我们重建了一个线性ODE系统d日=1250维，或第页=1563750未知参数。我们的变量选择算法产生了97.3%的稀疏网络。重建后的网络见补充表S2.

表4列出了此图中网络度最高的前十家公司（节点）。一个有趣的观察结果是，这些高度关联的公司大多不市值最大的公司，如苹果公司（Apple Inc.）或埃克森美孚（Exxon Mobile）。相反，其中四家公司提供基本的IT基础设施，如电话服务或网络硬件；其中两个与医疗服务有关；三家提供金融服务，也可以被视为现代经济的基础设施。总之，大多数关联公司不是标准普尔指数中最大或最著名的公司，而是为整个经济提供基础设施的公司。

表4：

根据网络度排名前十位最具影响力的公司，网络度定义为重构网络中节点的相邻边数。

	单位	类别	网络学位

TDS公司	电话和数据系统公司	IT基础设施	1162
紫外线辐射	环球保险控股公司	财务	696
AKRX公司	Akorn公司	能源	648
保密协议	纳斯达克OMX集团	财务	613
地理信息系统	通用磨坊	食物	589
ABAX公司	Abaxis公司	保健	552
CMTL公司	Comtech电信	IT基础设施	548
NTGR公司	Netgear公司	IT基础设施	545
NTCT公司	Netcout Systems公司	IT基础设施	530
SBRA公司	Sabra Health Care房地产投资信托基金	金融/医疗	515

在单独的窗口中打开

为了从更集中的角度更好地了解这些公司的相互作用，我们将股票划分为多个板块，并为每个板块重建子网络。更具体地说，我们下载了截至2015年10月12日由标准普尔索引的500家大型公司发行的股票清单，其中421个节点（公司）不是孤立节点。根据全球行业分类标准（GICS），这些公司进一步分为九个部门^性虐待在每个部门的子网络建设中，我们保留绝对强度大于所有边缘95%的边缘，以使结果在各个部门之间具有可比性。我们将枢纽定义为每个行业中连接最紧密（以网络度衡量）的前10%公司，这些公司在表5这些子网络在各个图中示出，并且作为一个压缩文件提供(补充文件S3)。

表5：

每个部门中最具影响力的公司由网络度定义，即其在重构网络中的相邻边数。

	单位	部门	学位
HAR公司	哈曼国际工业公司	消费者自由裁量权	20
AMZN公司	亚马逊网站股份有限公司	消费者自由裁量权	16
韦恩	永利度假村有限公司	消费者自由裁量权	13
地理信息系统	通用磨坊公司	消费必需品	8
ADM公司	Archer-Daniels-Midland公司	消费必需品	7

CHK（检查）	切萨皮克能源公司	能源	9
执行	钻石海上钻井	能源	7
RRC公司	Range Resources公司。	能源	7

保密协议	券交易所	金融	43
王牌	ACE有限公司	金融	19
FITB公司	五三银行	金融	16
TMK公司	Torchmark公司。	金融	12
BXP公司	波士顿地产	金融	10
PNC公司	PNC金融服务	金融	10

ENDP公司	Endo国际	医疗保健	13
JNJ公司	强生公司	医疗保健	10
SYK公司	史崔克公司。	医疗保健	9

FLS公司	福斯公司	工业	20
飞行情报室	FLIR系统	工业	16
ITW公司	伊利诺伊工具厂	工业	14

nflex公司	Netflix公司。	信息技术	11
支付（PAYX）	Paychex公司。	信息技术	11

BLL公司	Ball公司	材料	9

T型	美国电话电报公司	电信服务	三

AEE公司	Ameren公司	公用设施	5
NEE公司	NextEra能源	公用设施	5
SCG公司	SCANA公司	公用设施	5

在单独的窗口中打开

我们注意到一些有趣的结果来自表5例如，Harman International Industries和Amazon不出所料地是消费者自由裁量权领域联系最密切的两家公司，因为两家公司的产品种类繁多，可能会影响或受到其他行业领导者的影响。然而，令人惊讶的是，高端酒店和赌场的开发商和运营商永利度假村（Wynn Resorts）在这一类别的63家公司中排名第三。进一步调查显示，与永利度假村相关的所有13个连接都是向内的这意味着永利度假村高度依赖于该行业许多其他公司的业绩，但其股价对其他公司的影响不大。这一观察结果可能表明，我们可以将酒店和赌场的表现作为消费者支出整体适宜性的“试金石”。

5 讨论

本文提出了一种新的常微分方程参数估计和模型选择框架，该框架基于线性常微分方程系数矩阵的特征值估计，而不是直接估计其条目。这种新方法大大降低了相应非线性优化问题的维数第页=d日²+d日到d日以及其他d日²参数可以从不需要大量计算的封闭式公式中获得。因此，我们提出的算法比竞争过程更快、更稳定，并且可以轻松扩展以处理大型ODE系统。此外，我们对问题的重新计算提供了目标函数的闭合梯度，可以用于进一步加速和稳定计算。

在仿真研究中，我们证明了新的ST-SLS方法比竞争方法更稳定、更快地定位高维优化问题的全局解，从而提高了大型ODE系统的参数估计性能。我们新的ST-SLS估计方法和相应的变量选择算法的优越性能不仅是因为它具有显著的降维能力和目标函数的闭式梯度的可用性，还因为它有效地利用了耦合的ODE信息。

我们还将新算法应用于两个实际应用程序，以说明其在实践中的可用性；一个是30维的酵母细胞周期基因表达数据，另一个是1250维的标准普尔指数股价数据。我们的分析结果表明，新方法可以基于观测到的时间历程数据有效地恢复高维动态网络。

我们提出的方法适用于理论上可识别的一般高维线性ODE模型，但在实际实现中应注意一些问题。在实践中，如果系数矩阵的特征值是不同的，则线性常微分方程在理论上是可识别的；但ODE模型可能存在数值或统计可识别性问题(Miao等人，2011年)当多个特征值具有零或接近零的虚部时（例如，存在2个以上的实特征值），这是因为更多的实特征根表示ODE解中更多的指数项，并且指数项的幂很难在数值上区分和识别，这类似于线性回归中的多重共线性问题。还要注意，我们提出的方法要求不同数据点的数量大于ODE系统的维数，即。，n>d，尽管未知参数的总数第页=d日²+d日可以大于n个需要此要求以避免可识别性问题。通常，在应用我们的方法之前，必须处理可识别性问题。通常需要修改模型或组合一些变量来减少可识别性问题，但这超出了本文的范围。有动机的读者可以在中找到有关此主题的更多信息Miao等人（2011）.

在这个大数据时代，基于越来越多可负担的频繁时程数据，在一个大系统中的许多组件或元素之间建立动态关系，以便对复杂网络进行重构和分析，是一项常见的任务(Liu等人，2011年;Barabasi等人，2011年). 线性ODE系统是一个简单而强大的模型，可以用来描述大系统中元素之间的动态关系。本文中类似思想的未来扩展到高维非线性ODE系统(Wu等人，2014年)和/或具有部分观测变量的系统(Wu等人，2015年)尽管具有挑战性，但这是有保证的。我们认为，尽管高维复杂动态系统在实践中得到了广泛应用，但其识别和分析领域仍处于起步阶段。我们希望，我们的工作将推动这方面的更多研究。

补充材料

补充1

单击此处查看。^{（300K，pdf格式）}

补充2

单击此处查看。^{（24K，pdf）}

补充3

单击此处查看。^{（13K，pdf格式）}

补充4

单击此处查看。^{（13K，pdf格式）}

补充8

单击此处查看。^{（18K，pdf格式）}

补充9

单击此处查看。^{（11K，pdf）}

补充10

单击此处查看。^{（6.4K，pdf）}

补充11

单击此处查看。^{（9.3K，pdf）}

补充12

单击此处查看。^{（13K，xlsx）}

补充13

单击此处查看。^{（450万，xlsx）}

致谢

作者要感谢薛红旗博士对本文早期版本的深思熟虑的讨论和有益的评论，以及Michelle Carey博士对股市数据的建议。这项工作得到了NIH RO1 AI087135、呼吸病原体研究中心（NIAID合同号HHSN272201200005C）、罗切斯特大学CTSA奖号UL1 TR002001的部分支持，该奖由美国国立卫生研究院转化科学促进中心颁发，罗切斯特大学艾滋病研究中心（NIH 5 P30 AI078498-08）和国家自然科学基金资助11526096和11601185。

参与者信息

吴乐勤，暨南大学数学系，广州，中国。

邢秋，美国纽约罗切斯特罗切斯特大学生物统计与计算生物学系。

袁亚香，中国科学院数学与系统科学研究院，北京，中国。

胡林武，美国德克萨斯州休斯顿德克萨斯大学健康科学中心生物统计学系。

参考文献

艾哈迈德·N（1998），科学家和工程师的线性和非线性滤波，世界科学。[谷歌学者]
Aragon AD、Quiñones GA、Thomas EV、Roy S和Werner Washburne M.（2006），“在固定相酿酒酵母中释放抗萃取mRNA会导致转录物丰度大幅增加，以应对压力,”基因组生物学,7，R9。[PMC免费文章][公共医学][谷歌学者]
Barabasi A、Gulbahce N和Loscalzo J.（2011），“网络医学：基于网络的人类疾病治疗方法,”遗传学自然评论,12, 56–68.[PMC免费文章][公共医学][谷歌学者]
Benjamini Y.和Hochberg Y.（1995），“控制错误发现率：一种实用而有效的多重测试方法,”皇家统计学会期刊。B系列（方法学）,57, 289–300.[谷歌学者]
Black F.和Scholes M.（1973），“期权和公司负债的定价,”政治经济学杂志,81, 637–654.[谷歌学者]
Boer VM、de Winde JH、Pronk JT和Piper MD（2003），“在碳、氮、磷或硫限制的好氧恒化器培养基中，生长在葡萄糖上的酿酒酵母的全基因组转录反应,”生物化学杂志,278, 3265–3274. [公共医学][谷歌学者]
Bonneau R、Reiss DJ、Shannon P、Facciotti M、Hood L、Baliga NS和Thorsson V.（2006），“推断器：从系统-生物数据集从头学习节约型调控网络的算法,”基因组生物学,7，R36。[PMC免费文章][公共医学][谷歌学者]
Butcher J.（2014），“常微分方程”，in沃尔特·高茨基，施普林格，卷。三第7-8页。[谷歌学者]
Chen J.和Wu H.（2008a），“确定性动力学模型中时变系数的有效局部估计及其在HIV-1动力学中的应用,”美国统计协会杂志,103, 369–384.[谷歌学者]
陈杰（2008b），“确定性动态模型中时变参数的估计,”中国统计局,18, 987–1006.[谷歌学者]
Cherry JM、Adler C、Ball C、Chervitz SA、Dwight SS、Hester ET、Jia Y、Juvik G、Roe T、Schroeder M等人（1998年），“SGD：酵母菌基因组数据库,”核酸研究,26, 73–79.[PMC免费文章][公共医学][谷歌学者]
Combettes PL和Wajs VR（2005），“近端前向背向分裂信号恢复,”SIAM多尺度建模与仿真杂志,4, 1168–1200.[谷歌学者]
赞扬D、Jolly D、Drylewicz J、Putter H和Thiébaut R.（2011），”基于层次似然的HIV动力学模型推断,”计算统计与数据分析,55, 446–456.[谷歌学者]
De Jong H.（2002），“遗传调控系统的建模与仿真：文献综述,”计算生物学杂志,9, 67–103. [公共医学][谷歌学者]
丁安和吴华（2014），“用约束局部多项式回归估计常微分方程参数,”中国统计局,24，1613年至1631年。[PMC免费文章][公共医学][谷歌学者]
Fan J.和Li R.（2001），“基于非冲突惩罚似然的变量选择及其oracle性质,”美国统计协会杂志,96, 1348–1360.[谷歌学者]
Garber M、Grabherr MG、Guttman M和Trapnell C.（2011），“使用RNA-seq进行转录组注释和量化的计算方法,”自然方法,8, 469–477. [公共医学][谷歌学者]
Ginibre J.（1965），“复矩阵、四元数矩阵和实矩阵的统计系综,”数学物理杂志,6, 440.[谷歌学者]
Hemker PW（1972），“系统仿真和参数估计中微分方程的数值方法,”生化系统的分析与模拟,28, 59–80.[谷歌学者]
霍尔姆·S（1979），“一种简单的顺序拒绝多重测试程序,”斯堪的纳维亚统计杂志,6, 65–70.[谷歌学者]
Holter NS、Maritan A、Cieplak M、Fedoroff NV和Banavar JR（2001），“基因表达数据的动态建模,”美国国家科学院院刊,98, 1693–1698.[PMC免费文章][公共医学][谷歌学者]
Huang HL和Brandriss MC（2000），“酵母脯氨酸利用途径的调节器根据氮源的质量进行不同程度的磷酸化,”分子和细胞生物学,20, 892–899.[PMC免费文章][公共医学][谷歌学者]
Huang Y、Liu D和Wu H.（2006），“纵向HIV动态系统参数估计的层次贝叶斯方法,”生物计量学,62, 413–423.[PMC免费文章][公共医学][谷歌学者]
Huang Y.和Wu H.（2006），“HIV动态模型抗病毒疗效的贝叶斯估计方法,”应用统计学杂志,33, 155–174.[谷歌学者]
Huang Y、Wu H和Acosta EP（2010），“多治疗因素HIV动态微分方程模型的层次贝叶斯推断,”生物医学杂志,52, 470–486.[PMC免费文章][公共医学][谷歌学者]
Lavielle M、Samson A、Karina Fermin A和MentréF（2011），“HIV长期动态模型的最大似然估计与抗病毒反应,”生物计量学,67, 250–259.[PMC免费文章][公共医学][谷歌学者]
Lehmann N.和Sommers H-J（1991），“随机实矩阵的特征值统计,”物理审查信函,67，941–944页。[公共医学][谷歌学者]
Li Z、Li P、Krishnan A和Liu J.（2011），“结合微分方程模型和局部动态贝叶斯网络分析的大规模动态基因调控网络推理,”生物信息学,27, 2686–2691. [公共医学][谷歌学者]
Li Z、Osborne MR和Prvan T.（2005），“常微分方程的参数估计,”IMA数值分析杂志,25, 264–285.[谷歌学者]
Liang H.和Wu H.（2008），“基于回归模型测量误差框架的微分方程模型参数估计,”美国统计协会杂志,103, 1570–1583.[PMC免费文章][公共医学][谷歌学者]
Liu Y、Slotine JJ和Barabasi AL（2011），“复杂网络的可控性,”自然,473, 167–73. [公共医学][谷歌学者]
Lockhart DJ、Dong H、Byrne MC、Follettie MT、Gallo MV、Chee MS、Mittmann M、Wang C、Kobayashi M、Norton H等人（1996年），“高密度寡核苷酸阵列杂交监测表达,”自然生物技术,14, 1675–1680. [公共医学][谷歌学者]
Lopez F、Estève J-P、Buscail L、Delesque N、Saint-Laurent N、Théveniau M、Nahmias C、Vaysse N和Susini C（1997），“酪氨酸磷酸酶SHP-1与sst2生长抑素受体结合，是sst2介导的抑制性生长信号传导的重要组成部分,”生物化学杂志,272, 24448–24454. [公共医学][谷歌学者]
Lopez F、Ferjoux G、Cordelier P、Saint-Laurent N、Estève J-P、Vaysse N、Buscail L和Susini C（2001），“神经型一氧化氮合酶：SHP-1参与sst2生长抑素受体生长抑制信号传导的底物,”实验生物学联合会会刊,15, 2300–2302. [公共医学][谷歌学者]
Lu T、Liang H、Li H和Wu H（2011），“用于动态基因调控网络识别的高维ODE与混合效应建模技术,”美国统计协会杂志,106, 1242–1258.[PMC免费文章][公共医学][谷歌学者]
Luan Y.和Li H.（2003），“基于B样条混合效应模型的时程基因表达数据聚类,”生物信息学,19, 474–482. [公共医学][谷歌学者]
Ma P、Castillo-Davis CI、Zhong W和Liu JS（2006），“时间进程基因表达数据的数据驱动聚类方法,”核酸研究,34, 1261–1269.[PMC免费文章][公共医学][谷歌学者]
默顿RC（1973），“理性期权定价理论,”贝尔经济与管理科学杂志,4, 141–183.[谷歌学者]
Miao H、Xia X、Perelson AS和Wu H.（2011），“非线性常微分方程模型的可辨识性及其在病毒动力学中的应用,”SIAM审查,53, 3–39.[PMC免费文章][公共医学][谷歌学者]
Moler C.和Van Loan C.（2003年），“二十五年后，十九种计算矩阵指数的可疑方法,”SIAM审查,45, 3–49.[谷歌学者]
Montpetit B、Thorne K、Barrett I、Andrews K、Jadusing R、Hieter P和Measday V（2005），”全基因组合成致死筛选鉴定酿酒酵母核膜蛋白APQ12P和动粒之间的相互作用,”遗传学,171, 489–501.[PMC免费文章][公共医学][谷歌学者]
Øksendal B.（2013），随机微分方程：应用简介施普林格科学与商业媒体，第5版。[谷歌学者]
Powell M.（2006），“无导数无约束优化的NEWUOA软件,”大规模非线性优化,83, 255–297.[谷歌学者]
Poyton A、Varziri MS、McAuley KB、McLellan P和Ramsay J.（2006），“基于主微分分析的连续动态模型参数估计,”计算机与化学工程,30, 698–708.[谷歌学者]
Pramila T、Wu W、Miles S、Noble WS和Breeden LL（2006），“Forkhead转录因子Hcm1调节染色体分离基因并填补细胞周期转录回路中的S期缺口,”基因与发育,20, 2266–2278.[PMC免费文章][公共医学][谷歌学者]
推杆H、海斯特坎普S、兰格J和德沃尔夫F（2002），“HIV动力学模型参数估计的贝叶斯方法,”医学统计学,21, 2199–2214. [公共医学][谷歌学者]
Ramsay J.（1996），“主微分分析：微分算子的数据简化,”英国皇家统计学会杂志。B系列（方法学）,58, 495–508.[谷歌学者]
Ramsay J.和Silverman BW（1998），“功能数据分析,”统计与计算,8, 401–403.[谷歌学者]
Ramsay JO、Hooker G、Campbell D和Cao J.（2007），“微分方程的参数估计：广义平滑方法（附讨论）,”英国皇家统计学会杂志,69, 741–796.[谷歌学者]
Ruhe A.和Wedin PA（1980），“可分离非线性最小二乘问题的算法,”SIAM审查,22, 318–337.[谷歌学者]
Sakamoto E.和Iba H.（2001），“用遗传编程推导基因调控网络的微分方程组，“in2001年进化计算大会会议记录，体积。1第720-726页。[谷歌学者]
Schena M、Shalon D、Davis RW和Brown PO（1995），“用互补DNA微阵列定量监测基因表达模式,”科学类,270, 467–470. [公共医学][谷歌学者]
Šidák Z.（1967），“多元正态分布均值的矩形置信域,”美国统计协会杂志,62, 626–633.[谷歌学者]
Spellman PT、Sherlock G、Zhang MQ、Iyer VR、Anders K、Eisen MB、Brown PO、Botstein D和Futcher B（1998），“基因芯片杂交技术综合鉴定酿酒酵母细胞周期调控基因,”细胞分子生物学,9，3273–3297。[PMC免费文章][公共医学][谷歌学者]
Spieth C、Hassis N和Streichert F（2006），“网络推理问题的数学模型比较，《第八届遗传和进化计算年会论文集》，ACM，第279-286页。[谷歌学者]
Tanaka K、Mukae N、Dewar H、van Breugel M、James EK、Prescott AR、Antony C和Tanaka TU（2005），“纺锤体微管捕获动粒的分子机制,”自然,434, 987–994. [公共医学][谷歌学者]
Tao T.（2012），随机矩阵理论专题，卷。132，美国数学学会普罗维登斯，RI。[谷歌学者]
Tibshirani R.（2011），“通过套索进行回归收缩和选择,”英国皇家统计学会杂志,73, 267–288.[谷歌学者]
Varah JM（1982），“微分方程数值参数估计的样条最小二乘法,”SIAM科学与统计计算杂志,三, 28–46.[谷歌学者]
Voit EO（2000），生物化学系统的计算分析：生物化学家和分子生物学家的实用指南剑桥大学出版社。[谷歌学者]
Wang Z、Gerstein M和Snyder M（2009），“RNA-Seq：转录组学的革命性工具,”遗传学自然评论,10, 57–63.[PMC免费文章][公共医学][谷歌学者]
Wu H、Lu T、Xue H和Liang H（2014），“用于动态基因调控网络建模的稀疏可加ODE,”美国统计协会杂志,109, 700–716.[PMC免费文章][公共医学][谷歌学者]
Wu H、Miao H、Xue H、Topham D和Zand M.（2015），“基于部分可观测状态变量和时变参数的多元非线性ODE模型量化流感病毒感染的免疫反应,”生物科学统计学,7, 147–166.[PMC免费文章][公共医学][谷歌学者]
Wu S.和Wu H.（2013），“使用功能主成分分析方法对时间进程基因表达数据进行更强大的重要测试,”BMC生物信息学,14, 6.[PMC免费文章][公共医学][谷歌学者]
薛浩、苗浩和吴浩（2010），“同时考虑数值误差和测量误差的非线性常微分方程模型中常系数和时变系数的筛选估计,”统计年刊,38, 2351–2387.[PMC免费文章][公共医学][谷歌学者]
Yeung MS、Tegnér J和Collins JJ（2002），“基于奇异值分解和稳健回归的基因网络逆向工程,”美国国家科学院院刊,99，6163–6168。[PMC免费文章][公共医学][谷歌学者]
Yukawa M、Tawara Y、Yamagishi M和Yamada I（2012），“基于lp形式软阈值技术的稀疏自适应滤波器2012年IEEE国际电路与系统研讨会，IEEE，第2749–2752页。[谷歌学者]
张驰（2010），“MINIMAX凹形惩罚下的近似无偏变量选择,”统计年刊,38, 894–942.[谷歌学者]
Zhang H、Conn AR和Scheinberg K.（2010），“最小二乘最小化的无导数算法,”SIAM优化杂志,20, 3555–3576.[谷歌学者]

线性常微分方程大系统的参数估计和变量选择：基于矩阵的方法

吴乐琴

邢秋

雅香园

胡林武

关联数据

摘要

1 简介

2 模型和方法

2.1、。型号说明

2.2. 相似变换与可分离最小二乘

备注。

定理2.1。

备注。

定理2.2。

推论2.3。

2.3. 渐近方差与推断

定理2.4。

定理2.5。

2.4. 变量选择

三。 模拟研究

3.1. 模拟实验设计

3.2. 参数估计比较

表1：

3.3. 变量选择比较

表2：

表3：

4 实际数据分析

4.1. 时程酵母芯片数据分析

4.2。标准普尔股市数据分析

表4：

表5：

5 讨论

补充材料

补充1

补充2

补充3

补充4

补充5

补充6

补充7

补充8

补充9

补充10

补充11

补充12

补充13

致谢

参与者信息

参考文献

三。模拟研究