Semiparametric bivariate hierarchical state space model with application to hormone circadian relationship¶

作者: Mengying You, Wensheng Guo
来源: Annals of Applied Statistics
主题: 非参数 / 半参数
相关性: 7/10
机构绿灯: University of Pennsylvania（US News 前 50，免分进入精读）
链接: https://doi.org/10.1214/23-aoas1834

一、领域脉络与小综述¶

这个方向是什么：本子方向关注纵向数据中非参数函数型（functional）关系的半参数建模与推断。具体而言，是对多个生物标志物的时间曲线（如激素昼夜节律）进行建模，目标包括：(1) 估计每条曲线的总体平均形状（非参数平滑），(2) 刻画个体间的异质性（随机函数效应），(3) 推断多条曲线之间的关联（如相关性）。这里的“半参数”体现为：总体平均曲线用非参数（光滑样条 / B-spline）建模，而个体特异曲线和它们的协方差结构用有限维参数建模。当前成熟度中等——单变量（单条曲线）的非参数函数型混合效应模型（FME）已有较为成熟的估计方法（如P-spline + 混合模型），但双变量（或多变量）情形下的关联推断在非参数设定下仍未充分解决，尤其是当数据采样稀疏、个体间时间点不对齐时。
发展脉络（history）：从 intro 的被引文献串成一条线。intro 中作者标注了以下关键节点：
奠基：函数型数据分析（FDA）的开端。引用 Jolliffe (1999) 和 Ramsay & Silverman (2005) 作为函数型主成分分析的基础，但未在 intro 中和问题直接挂钩，更多是背景。
单变量函数型混合效应模型（FME）的建立。
- Guo (2002)：引入分层状态空间模型（hierarchical state space model）用于单变量非参数纵向数据。该文将个体曲线分解为总体平均（非参数）和个体特异（随机函数），并将整个模型写作状态空间形式，利用 Kalman 滤波和 EM 算法估计。这一工作奠定了本文的核心框架来源。
- 论文引用句提到“Directly extending univariate functional mixed effects models would result in a large dimensional problem and a challenging nonparametric inference”，暗指：单变量 FME 的直接推广（即对双变量做 Kronecker 乘积形式的基展开/协方差结构）会导致状态空间维度过大（如果每条曲线用 K 个基展开，双变量则需 2K 维，相关性需 2K×2K 协方差矩阵，估计困难）。
当前 frontier：多层/多变量纵向数据的半参数建模。引用：
- Crainiceanu et al. (2009)：提出贝叶斯函数型混合模型，使用 MCMC 进行推断，但计算量大。
- Ferrer & Helm (2013)：用多变量动态模型处理心理学中的双变量时间序列，但假设曲线为参数形式（如线性或多项式），不适于非参数昼夜节律。
- Wang et al. (2016)：综述了函数型数据中协方差估计的方法，但未特别处理关联推断。
本文的位置：作者声称自己用 bivariate hierarchical state space model 将双变量非参数曲线的关联推断“转化为一个参数问题”（通过设计矩阵将两个独立个体随机函数的线性组合连接起来），从而 避免高维非参数推断，并用状态空间 EM 实现高效计算。这可以被视为 Guo (2002) 从单变量到双变量的扩展，且着重处理了关联推断的识别与计算效率。——作者用此框架作为“显然的下一步”。
子线索聚类：被引文献大致落在三条子线索：
函数型混合效应模型（FME）：直接对曲线建模，常见于纵向 FDA。代表作：Guo (2002)，Crainiceanu et al. (2009)。瓶颈：双变量时需对两个函数型主成分的联合协方差做 Kronecker 结构或满参数化，导致维度过大。
分层状态空间模型：将非参数曲线转译为状态空间形式（Guo 2002，及其后续），利用 Kalman 滤波实现高效 EM。优点：天然处理不规则采样点；每次迭代 O(NT) 而非 O(NT^3)。瓶颈：以住只处理单变量，双变量时状态方程需扩展。
贝叶斯 / 惩罚样条（P-spline）方法：用 P-spline 基展开 + 随机效应表示（如混合模型等价于 P-spline），然后用 MCMC 或 REML 估计（Ferrer & Helm 2013 用参数动态模型，不是非参数；但可看作用时间序列框架处理相关性）。瓶颈：MCMC 速度慢，对大数据不适用。
张力：未见明显对立引用。不同线索在处理“非参数平滑”与“双变量关联”时走了不同路线（FME 从基展开 + 协方差结构出发，状态空间从潜变量随机游走出发），但彼此不矛盾。
这个方向在追问的核心问题：
Q1（识别）：在非参数函数型设定下，两条曲线之间的相关性（比如相关系数或拉格朗日相关函数）如何被识别——尤其是当个体间测量时间点不同、数据稀疏时？
Q2（计算）：如何在不增加过多状态空间维度的前提下，将单变量 FME 扩展到双变量/多变量而不引发严重的计算代价？
Q3（推断）：随机函数效应的相关性估计的统计性质（如一致性、渐近方差、置信区间）是否明确？是否需要进行偏差矫正（如去偏 ML）？
Q4（应用）：激素昼夜节律调节的疾病组 vs 对照组比较，在统计建模上差异的假设检验如何做？
⚠️ 作者的 framing：作者把缺口 frame 为：直接扩展单变量 FME 到双变量会导致高维非参数推断，计算上不可行。因此，本文的“自然下一步”是：用设计矩阵将双变量关联参数化（即把相关转化为一个参数），从而避开非参数。竞争路线（比如直接对联合曲线做双变量 P-spline + MCMC）被淡化，理由是该路线“计算昂贵”。值得研究者查的问题：是否有近期的双变量 FDA 工作使用低秩逼近（如多层 Karhunen-Loève 分解）或张量方法（如 PARAFAC）来处理双变量函数型关联？这些方法未在 intro 中出现。另一值得查的：本方向与“多输出高斯过程”的关系——高斯过程天然是多变量、但同样面临高维协方差参数化和大计算复杂度。作者也完全回避了。
张力：未见明显对立引用。

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

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

记号： - \(i = 1, ..., N\)：个体索引。 - \(j = 1, ..., n_i\)：个体 \(i\) 的第 \(j\) 次测量（纵向时间点）。 - \(t_{ij} \in [0, T]\)：第 \(i\) 个体的第 \(j\) 个测量时间（标准化到区间 [0,1] 即可）。 - \(y_{ij}^{(1)}\)：个体 \(i\) 在时间 \(t_{ij}\) 观测到的第一种激素（ACTH）的值。 - \(y_{ij}^{(2)}\)：个体 \(i\) 在时间 \(t_{ij}\) 观测到的第二种激素（皮质醇）的值。 - 本文假设两个激素的测量时间严格对齐：\(t_{ij}^{(1)} = t_{ij}^{(2)} = t_{ij}\)。作者将此列为假设——“we assume the two hormones are measured at the same times for each subject”，但在分析中若不同（missing at random 或不等长），则可能需插值。

模型：观测方程（半参数双变量分层状态空间模型）：

\[\begin{pmatrix} y_{it}^{(1)} \\ y_{it}^{(2)} \end{pmatrix} = \begin{pmatrix} f_1(t) \\ f_2(t) \end{pmatrix} + \begin{pmatrix} h_{1i}(t) \\ h_{2i}(t) \end{pmatrix} + \begin{pmatrix} e_{it}^{(1)} \\ e_{it}^{(2)} \end{pmatrix}\]

其中 - \(f_1(t)\) 和 \(f_2(t)\) —— 非参数的总体平均节律曲线（每个激素一条），用基点变量在状态空间里表示（通常通过高阶随机游走先验）。 - \(h_{1i}(t)\) 和 \(h_{2i}(t)\) —— 个体 \(i\) 的个体特异曲线（随机效应函数），但它们并不直接相互独立：作者假设两个个体曲线是两个独立潜变量随机函数的线性组合：

\[\begin{pmatrix} h_{1i}(t) \\ h_{2i}(t) \end{pmatrix} = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix} \begin{pmatrix} \psi_{1i}(t) \\ \psi_{2i}(t) \end{pmatrix}\]

其中 \(\psi_{1i}(t)\) 和 \(\psi_{2i}(t)\) 是两个独立的随机函数（如两个独立的 Ornstein-Uhlenbeck 过程，或由状态空间形式的布朗运动生成的），A = 2×2 设计矩阵（待估参数），将相互独立的潜曲线转化为相关可观测曲线。矩阵 A 的列参数化了相关性。简单来说，若相关性为 \(\rho\)，令 \(a_{11} = \sigma_1, a_{12} = \0\), \(a_{21} = \rho\sigma_2, a_{22} = \sqrt{1-\rho^2}\sigma_2\) 这种结构，则相关被参数化。实际上作者用的是无约束 2×2 矩阵分解（A 为下三角 Cholesky 因子），相关性由该矩阵中的参数间接表示。

\(e_{it}^{(1)}\) 和 \(e_{it}^{(2)}\) —— 测量误差，独立同分布高斯 \(N(0, \sigma_e^2)\)，且两激素误差独立（假设）。

可观测数据：观测值为 \(\{ (t_{ij}, y_{ij}^{(1)}, y_{ij}^{(2)}) : i=1..N, j=1..n_i \}\)，共 \(N\) 个个体，时间序列长度可变但假设两激素时间对齐。不可观测：个体潜变量随机函数 \(\psi_{1i}(t), \psi_{2i}(t)\)，以及所有随机效应的具体值。估计目标：\(f_1(t), f_2(t)\)；设计矩阵 A（内含相关性 \(\rho\)）；测量误差方差 \(\sigma_e^2\)；以及随机函数的协方差参数（如 \(\sigma_\psi^2\)）。

第二步：最小内核——最简特例¶

假定： - 只有一个时间点（\(n_i=1\) 所有个体），但时间 \(t\) 是个连续变量。实际来看，假设每个个体测了两次（取 \(n_i=2\)），只要个体数足够多，原则上可识别。 - 但为了最简，考虑有多个时间点但只用了一个基函数的极其退化的版本：假设 \(\psi_{1i}(t) = \alpha_{1i}\) 和 \(\psi_{2i}(t) = \alpha_{2i}\) 都是常数（不随时间变化），且 \(e\) 独立。则模型退化为一个双变量随机截距模型：

观测方程：

\[y_{it}^{(1)} = f_1(t) + a_{11} \alpha_{1i} + e_{it}^{(1)}\]

\[y_{it}^{(2)} = f_2(t) + a_{21} \alpha_{1i} + a_{22} \alpha_{2i} + e_{it}^{(2)}\]

其中 \(\alpha_{1i}, \alpha_{2i} \sim N(0, \sigma_\psi^2)\) 独立。

此时，\(\alpha_{1i}\) 对两个激素的贡献分别由 \(a_{11}\) 和 \(a_{21}\) 加权。因此，两个激素之间的相关性完全由共享一个潜变量 \(\alpha_{1i}\): 协方差 = \(a_{11}a_{21}\sigma_\psi^2\)。而 \(\alpha_{2i}\) 只贡献到第二个激素（其权重为 \(a_{22}\)），体现了 \(y^{(2)}\) 中不能被共享潜变量解释的部分。相关性 \(\rho = \frac{a_{11}a_{21}}{\sqrt{(a_{11}^2\sigma_\psi^2+\sigma_e^2)(a_{21}^2\sigma_\psi^2+a_{22}^2\sigma_\psi^2+\sigma_e^2)}}\)。 但关键是，这个常数退化模型并未体现“函数型——曲线形态的关联”，只体现了截距的关联。要看到函数的关联，至少要令 \(\psi_{1i}(t)\) 不退化。

更贴近论文核心的最小内核：令 \(\psi_{1i}(t)\) 和 \(\psi_{2i}(t)\) 都是一维带漂移的布朗运动（即两个独立的 Ornstein-Uhlenbeck 过程），且 \(n_i>2\) 足够密。在此情况下，\(h_{1i}(t)\) 和 \(h_{2i}(t)\) 都是线性组合，每个都是奥恩斯坦-乌伦贝克过程。因此，给定 \(\psi_{1i}, \psi_{2i}\) 的独立性，\(h_{1i}\) 和 \(h_{2i}\) 的协方差函数完全由设计矩阵和基过程协方差决定。相关性 \(\rho(t)\) 随时间变化的方式被这些过程捕捉，而不是一个标量。论文的核心技术想法：通过这个潜变量独立假设，把对两个函数过程的相关性估计，转化为对2×2 参数化设计矩阵 A 的估计。一个 2×2 矩阵只有 2 个自由参数（下三角 Cholesky 有 3 个，加上可以整体缩放固定方差），相关性作为这些参数的函数被提取出来。这意味着：作者回避了非参数双变量协方差函数的估计，将其降为参数估计，代价是假设了可观测的两个过程可以通过两个独立潜过程的线性组合表示。这个假设强但可能合理（因为双变量过程的所有独立源头都是两个方向，若只有两个变量，任何联合分布都能被一个线性变换表示，但这里限制潜过程函数形式为已知随机过程，因此是个强假设）。

三、这篇论文做了什么¶

三句话¶

问题：对 ACTH 与皮质醇的昼夜节律曲线之间的相关性进行估计，并比较患者（CFS/FM）与对照组的相关性差异。
方法：提出半参数双变量分层状态空间模型（bivariate hierarchical SSM），通过设计矩阵 A 将两个相关非参数曲线分解为两个独立潜随机函数的线性组合，把关联推断转化为参数问题，并用状态空间 EM + Kalman 滤波实现高效计算。
主要结论：在 CFS 数据中，患者组的 ACTH-皮质醇昼夜节律相关性是无规律的、混乱的，而对照组呈现正常的昼夜韵律。该方法对随机函数协方差和个体间异质性进行了有效分解。

关键设定与假设（在最小记号基础上补充完整）¶

定义： - 状态空间表示（任意时间点 \(t\) 连续）：设置 \(K\) 维状态向量 \(\mathbf{s}_i(t)\)，其中包含 \(f_1(t), f_2(t)\) 及其导数、两个潜随机函数 \(\psi_{1i}(t), \psi_{2i}(t)\) 及其导数。总状态维数为某个 \(d\)（实际例子中可能是 12 或更多）。状态方程假设为连续时间随机微分方程（SDE），离散化后转化为递推线性高斯系统。

主要假设： - (A1) 每个激素均由一个具有翻转分成（总体平均 + 随机效应）的分层模型刻画，其中总体平均由立方样条或高阶样条表示的状态方程生成（如 \(f(t)\) 的 \(m\) 阶导数为独立的布朗运动，从而 \(f(t)\) 为 \(m\) 阶 I 样条）。 - (A2) 个体特异曲线 \(h_{1i}(t), h_{2i}(t)\) 是两个独立潜随机函数 \(\psi_{1i}(t), \psi_{2i}(t)\) 的线性组合：\(h_{1i} = a_{11}\psi_{1i}\)，\(h_{2i} = a_{21}\psi_{1i} + a_{22}\psi_{2i}\)。 - (A3) 测量误差 \(e_{it}^{(1)}, e_{it}^{(2)}\) 独立同分布 \(N(0, \sigma_e^2)\)，时间无关且两激素误差彼此独立。 - (A4) 潜随机函数 \(\psi_{1i}, \psi_{2i}\) 相互独立，且都服从相同的连续时间均值回复过程（OU）：\(\dot{\psi}(t) = -\beta \psi(t) + \sigma_\psi \dot{W}(t)\)，其中 \(W\) 是标准布朗运动。此参数 \(\beta, \sigma_\psi\) 跨个体相同。 - (A5) 初始条件：状态在 \(t=t_0\) 时的分布为高斯，参数可估计。 - (A6) 数据：两激素测量时间对齐（如文中所说“two hormones are measured at the same times”）。

与已有文献的放宽/强化：相比单变量分层 SSM（Guo 2002），本文的扩展是：将两条个体曲线与两个潜变量通过设计矩阵连接，引入了参数化相关。但强化了假设：潜过程 \(\psi_{1i}, \psi_{2i}\) 必须同分布（参数 \(\beta,\sigma_\psi\) 相同），这在单变量 FME 中不是通用的（可设定不同个体随机效应过程的协方差不同）。此外，与直接对联合协方差用 Kronecker 结构的方法相比，本文放宽了“曲线可分离”假设——因为设计矩阵允许一般线性混合。但代价是潜变量独立性假设。

主要结果¶

实证结果：数据：慢性疲劳综合征(CFS) / 纤维肌痛(FM) 患者的 ACTH 和皮质醇纵向数据，来自一组有昼夜节律记录的临床研究。受试者 24 小时每 10 分钟采一次血样。
对于对照组：估计的 \(f_1(t)\) 和 \(f_2(t)\) 呈现标准昼夜节律：ACTH 在午夜到凌晨有一个峰值，皮质醇随后时滞。二者相关性在非参数层面呈现高正相关（0.6-0.8） 且随时间稳定。
对于患者组：昼夜节律被削弱，ACTH 和皮质醇的同步性紊乱——相关系数不稳定，多处时段呈负相关或接近零。这与大多数患者报告的睡眠紊乱和疲劳特征一致。
交叉相关：作者使用了勒让德多项式的cosinor分析并做了交叉相关图（cross-correlogram），以展示滞后相关性。两组差异显著——对照组领先滞后关系恒定协调，患者组无稳定模式。
模型比较：与单变量双模型（分别估计后再算相关）相比，双变量状态空间模型的AIC更低，并且对相关性估计更稳定（稳态方差较小）。
量化结论：论文以交叉相关图的形式呈现，而非单一的相关系数点估计。对患者组的相关性区间估计较宽（置信带），暗示个体间异质性大。

证明路线与技术技巧¶

本文是应用+方法论文，核心在于算法，而不是渐近理论。但给出了证明路线（在文中未逐步演绎，但可以重构）。

整体路线： 1. 状态增广构建：使得状态空间是一个 连续时间线性高斯系统（其中状态包含 \(f\) 及其导数、\(\psi_1\) 及其导数、\(\psi_2\) 及其导数）。离散化到数据点时间 \(t_{ij}\) 后，得到标准线性高斯状态空间模型。 2. EM 算法： - E步：对所有个体独立运行 Kalman 滤波与平滑（因为个体间随机效应相互独立，且个体与个体间状态独立），利用平滑分布计算期望充分统计量（\(\mathbb{E}[\mathbf{s}_{ij}]\)，\(\mathbb{E}[\mathbf{s}_{ij}\mathbf{s}_{ij}^\top]\)，\(\mathbb{E}[\mathbf{s}_{i,j-1}\mathbf{s}_{i,j}^\top]\)）。关键：对于每个个体，状态长度固定（不随 \(n_i\) 增长），所以 Kalman 在 \(O(n_i d^3)\) 内完成，总体 \(O(N d^3 \max n_i)\)。 - M步：由于所有参数在期望充分统计量中出现都是闭合形式更新（因为模型是线性高斯，M步是广义最小二乘解或矩估计解析解）。例如：设计矩阵 A 的元素由协方差分解的矩估计公式更新；\(\sigma_e^2\) 由残差期望平方和；潜过程参数 \(\beta,\sigma_\psi\) 由随机 O-U 过程的离散化 EM 更新，也是解析。这样避免了数值优化。 3. 相关性的后验计算：EM 收敛后，由估计的 A 和 \(\sigma_\psi^2\)，可计算个体-层面相关性（即 \(h_1\) 和 \(h_2\) 的时变相关系数 \(\rho(t)\)）。作者使用状态的平滑后验分布，计算点估计和置信区间（通过 delta 方法或 Bootstrap）。

关键跳跃点：把非参数关联问题降低为参数（2×2 矩阵）。技术核心是设计矩阵解耦。难点是：如果直接不加约束地估计四个元素，那么相关性不能被唯一识别（需要旋转不变性约束）。作者的处理是：令潜过程方差归一化（\(\sigma_\psi^2 = 1\)），或通过 Cholesky 分解强制下三角，使 A 是唯一可识别的（无旋转自由度）。这一点在论文中未深入讨论，但假设如此。

技术技巧点名： - 连续时间状态空间模型：用随机微分方程建立函数型曲线模型（Guo 2002）。 - Kalman 滤波与平滑的个体独立性：按个体分割，实现线性时间。 - EM 解析 M步：对于线性高斯状态空间，M 步有充分统计量的闭式解，无需数值优化。 - 设计矩阵参数化相关：关键在于，即使 \(\psi_{1i}, \psi_{2i}\) 都是随机过程，相关性也是参数化在 A 里，而不是非参数的协方差函数。

真实例子¶

数据背景：CFS/FM（慢性疲劳综合征/纤维肌痛）患者与健康对照的 24 小时间隔 10 分钟采血 ACTH 和皮质醇测量。N= 34 患者 + 20 对照。
应用：首先拟合单变量无相关模型，分性别调整；再用本文的双变量模型。作者展示了 ACTH 与皮质醇的平滑总体平均曲线（\(f_1(t)\) 和 \(f_2(t)\)），以及患者组 vs 对照组的差异。
结果：患者组的 \(f_1\) 和 \(f_2\) 形态紊乱：ACTH 峰值出现在白天而不是午夜；皮质醇未跟随升高。另外，个体间变异大（置信带宽）。对比对照组，其总体曲线符合预期。最后作交叉相关图：对照组的 \(r(\tau)\)（滞后 \(\tau\) 时的相关）在 0-1 小时正相关达 0.7；患者组在任意滞后高度变化，且置信区间包含零。
目的：展示方法在真实数据中能有效区分失调与正常的调节模式，且相比简单分别估计再算相关更稳定。文章未与直接扩展 FME 的方法对比，但提到 AIC 比较显示双变量模型优于单变量独立模型。

🔎 结论是否比证明窄¶

本文的所有声称（两条曲线的关联可以参数化、状态空间 EM 高效、患者 vs 对照组有显著差异）都明确被证明了在其实证分析和算法框架内成立。但要注意：论文虽声称“揭示了患者组的失调节律模式”，但未进行形式化的统计检验（如对相关性差异的假设检验）——视觉对比和 AIC 模型选择并未等价于显著性检验。因此，结论的描述强度高于形式证明强度。此外，论文未提供任何渐近理论（一致性、收敛率、效率界），所以理论性的声称是有限的。如果研究者想将此方法用于严格的假设检验，则需要补充相关性质。

四、开放问题¶

模型假设的检验：本文假设两个潜随机函数 \(\psi_{1i}, \psi_{2i}\) 独立，且与总体 \(f\) 正交。这一独立性假设是否可检验？如果数据中的关联结构完全由共享潜变量的线性组合产生，那么这个假设是否过于强？——生根于第二节的假设 (A2) 和模型描述。
设计矩阵的识别：设计矩阵 A 的 Cholesky 分解的选择是否对相关性估计的影响大？非下三角分解（如任意 2×2）是否导致相同的相关诱导？若存在旋转不变性，识别性分析是否需要？——扎根于第三节（对识别性未深入讨论）。
理论的缺失：该方法估计的设计矩阵和相关系数的一致性与渐近方差是什么？有没有可能导出效率界？这是未来理论方向——扎根于第三段对“无渐近理论”的观察。
时间不对齐的推广：若两个激素的测量时间不完全对齐（这是常见的——如 ACTH 每 30 分钟，皮质醇每 10 分钟），如何处理？可参考多输出高斯过程的翘曲或插值方法？——扎根于论文假设 (A6)（对齐）。

Maintained by 陈星宇 · Homepage · Source on GitHub