跳转至

Modelling change processes in multivariate interrupted time series data using a multivariate dynamic additive model: An application to heart rate and blood pressure self-monitoring in heart failure with drug changes

作者: Sun-Joo Cho
来源: Journal of the Royal Statistical Society Series C
主题: 流行病学
相关性: 3/10
机构绿灯: Vanderbilt University（US News 前 50，免分进入精读）
链接: https://doi.org/10.1093/jrsssc/qlad088

一、领域脉络与小综述

这个方向是什么

本文研究的核心问题是：在多元中断时间序列（multivariate interrupted time series, MITS）数据中，如何同时建模并估计由干预（如药物变化）引起的水平变化（level change）、趋势变化（trend change，包括线性和非线性）以及时变序列依赖（time-varying serial dependence）。该方向处于应用统计方法学的成熟阶段——已有大量针对单变量中断时间序列（ITS）的模型，但多元情形下同时处理上述三种变化过程的系统性方法仍属空白。本文的贡献在于将动态加性模型（dynamic additive model）扩展到多元情形，并特别关注序列依赖的时变性对变化检测的影响。

发展脉络（history）

根据论文引言（作者亲手画的领域地图）和参考文献，该方向的发展脉络可梳理如下：

  1. 奠基工作：单变量中断时间序列分析（ITS）

    • Box & Tiao (1975)：提出干预分析（intervention analysis）框架，使用ARIMA模型来估计干预对单变量时间序列的影响。这是ITS分析的统计基础。
    • McDowall et al. (1980) 和 Crosby et al. (1993)：将ITS方法系统化并推广到社会科学和公共卫生领域，使其成为评估政策或干预效果的常用准实验设计。
    • Linden (2015)：提出使用分段回归（segmented regression）来估计ITS中的水平变化和斜率变化，这是目前应用最广泛的方法之一。其核心假设是序列依赖（如自回归结构）在干预前后是恒定的。

  2. 主要进展：处理非线性趋势和时变序列依赖

    • Huitema & McKean (2000) 和 Ferron et al. (2017)：指出在单变量ITS中，忽略序列依赖（如自相关）会导致标准误估计偏小，从而夸大干预效果的统计显著性。这推动了更稳健的标准误估计方法（如Newey-West估计）的使用。
    • McNeish & Harring (2017) 和 Shadish et al. (2014)：开始探索使用加性模型（additive models）或样条（splines）来拟合ITS中的非线性趋势，突破了分段回归只能处理线性趋势的限制。
    • Cho et al. (2021)：作者本人的前期工作，提出了一个单变量动态加性模型（univariate dynamic additive model），该模型可以同时处理水平变化、非线性趋势变化和时变序列依赖。这是本文的直接前身。

  3. 当前Frontier：从单变量到多变量

    • 本文 (Cho, 2024)：将Cho et al. (2021)的单变量模型扩展到多元情形，以处理多个相关结果变量（如心率和血压）同时受干预影响的情况。这是该子方向的一个自然且必要的延伸。

子线索聚类

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

  • 线索一：基于分段回归的经典ITS方法

    • 做什么：假设干预前后趋势为线性，序列依赖（如AR(1)）恒定。通过线性回归或广义最小二乘估计水平变化和斜率变化。
    • 代表工作：Linden (2015), Huitema & McKean (2000), Ferron et al. (2017)。
    • 瓶颈：无法处理非线性趋势；假设序列依赖不随时间变化，这在长期自我监测数据中可能不成立。

  • 线索二：基于加性模型的灵活ITS方法

    • 做什么：使用样条或核方法对趋势进行非参数建模，允许序列依赖（如自回归参数）随时间变化。
    • 代表工作：McNeish & Harring (2017), Cho et al. (2021), 本文 (Cho, 2024)。
    • 瓶颈：计算复杂度较高；扩展到多元时，需要处理结果变量之间的同期相关（contemporaneous correlation）和交叉滞后依赖（cross-lagged dependence），模型识别和参数估计更具挑战。

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

  1. 如何准确估计干预引起的水平变化和趋势变化？ 核心挑战在于将干预效应与时间趋势、序列依赖、季节性等混淆因素分离开。
  2. 如何处理时间序列中的时变序列依赖？ 当序列依赖（如自回归系数）本身随干预或时间变化时，忽略它会导致对干预效应的估计产生偏倚。
  3. 如何将单变量ITS模型扩展到多元情形？ 当多个结果变量同时被观测且相互关联时，如何建模它们之间的同期和滞后关系，并同时估计每个变量上的干预效应？
  4. 模型的可识别性与计算可行性如何？ 随着模型复杂度（非线性趋势、时变依赖、多元结构）增加，如何保证参数可识别，并开发出高效的计算算法？

⚠️ 作者的Framing（必须明确标注成"这是作者的说法"）

  • 作者把缺口frame成什么？ 作者在引言中明确指出：“现有研究主要集中在单变量ITS上，很少有研究处理多元ITS中的变化过程”（原文：Existing research has focused primarily on univariate ITS, and few studies have dealt with change processes in multivariate ITS）。因此，本文被定位为“显然的下一步”：将作者之前开发的单变量动态加性模型（Cho et al., 2021）推广到多元情形，以解决一个实际应用（心衰患者自我监测）中出现的、但现有方法无法处理的问题。
  • 哪些竞争路线被他淡化或回避了？ 作者淡化了结构方程模型（SEM）或向量自回归模型（VAR）在多元时间序列分析中的应用。虽然VAR模型也能处理多元时间序列，但作者认为其“通常假设时间序列是平稳的，且干预效应是瞬时的或通过脉冲响应函数体现”，这与ITS中关注的“持续的水平变化和趋势变化”设定不同。作者也回避了贝叶斯方法（如动态线性模型DLM）的讨论，尽管DLM也能处理时变参数。
  • 什么明显该被引/该存在、却没出现在intro里？ 一个明显的缺失是因果推断领域关于中断时间序列的识别假设的讨论（如无混淆假设、无同期干预等）。本文完全从统计建模（拟合优度、参数恢复）角度出发，没有讨论其估计量是否具有因果解释。例如，没有引用Shadish, Cook, & Campbell (2002) 关于准实验设计的经典著作，也没有引用Ramsay et al. (2003) 关于ITS在公共卫生中应用的方法学指南。这暗示本文更侧重于描述性或预测性建模，而非严格的因果推断。

张力

未见明显对立引用。该领域的发展是渐进的，从简单到复杂，从单变量到多变量，没有出现不同条件下得出相反结论的激烈争论。

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

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

  • 符号：

    • \( t = 1, \dots, T \)：时间点（离散、等间隔）。
    • \( i = 1, \dots, N \)：个体（患者）。
    • \( j = 1, \dots, J \)：结果变量（如J=2，心率和血压）。
    • \( y_{ijt} \)：第 \( i \) 个个体在第 \( t \) 个时间点的第 \( j \) 个结果变量的可观测值。
    • \( \mathbf{y}_{it} = (y_{i1t}, \dots, y_{iJt})^\top \)：第 \( i \) 个个体在时间 \( t \) 的 \( J \times 1 \) 可观测结果向量。
    • \( \boldsymbol{\mu}_{it} = (\mu_{i1t}, \dots, \mu_{iJt})^\top \)：第 \( i \) 个个体在时间 \( t \) 的 \( J \times 1 \) 潜在均值向量（即趋势项）。
    • \( \boldsymbol{\epsilon}_{it} = (\epsilon_{i1t}, \dots, \epsilon_{iJt})^\top \)：第 \( i \) 个个体在时间 \( t \) 的 \( J \times 1 \) 随机误差向量。
    • \( \mathbf{\Phi}_{it} \)：第 \( i \) 个个体在时间 \( t \) 的 \( J \times J \) 时变自回归系数矩阵。其对角线元素 \( \phi_{ijj,t} \) 表示变量 \( j \) 的滞后一阶自回归系数；非对角线元素 \( \phi_{ijk,t} \) 表示变量 \( k \) 对变量 \( j \) 的滞后一阶交叉影响。
    • \( \mathbf{\Sigma}_{it} \)：第 \( i \) 个个体在时间 \( t \) 的 \( J \times J \) 时变误差协方差矩阵。其对角线元素 \( \sigma^2_{ij,t} \) 表示变量 \( j \) 的误差方差；非对角线元素 \( \sigma_{ijk,t} \) 表示变量 \( j \) 和 \( k \) 的误差在时间 \( t \) 的同期相关。
    • \( \mathbf{x}_{it} \)：第 \( i \) 个个体在时间 \( t \) 的 \( p \times 1 \) 协变量向量（如时间、干预指示变量、季节虚拟变量等）。
    • \( \boldsymbol{\beta}_j \)：第 \( j \) 个结果变量的 \( p \times 1 \) 回归系数向量（固定效应）。

  • 模型： 本文的核心模型是一个多元动态加性模型（Multivariate Dynamic Additive Model, MDAM）。其基本结构是：

    \[\mathbf{y}_{it} = \boldsymbol{\mu}_{it} + \boldsymbol{\epsilon}_{it}\]
    其中，均值 \( \boldsymbol{\mu}_{it} \) 由加性预测器（additive predictor）建模：
    \[\mu_{ijt} = \mathbf{x}_{it}^\top \boldsymbol{\beta}_j + \sum_{k=1}^{K} f_{jk}(z_{itk})\]
    这里 \( f_{jk}(\cdot) \) 是未知的光滑函数（如样条），用于建模非线性趋势。误差项 \( \boldsymbol{\epsilon}_{it} \) 被建模为一个时变向量自回归过程（TV-VAR(1)）：
    \[\boldsymbol{\epsilon}_{it} = \mathbf{\Phi}_{it} \boldsymbol{\epsilon}_{i,t-1} + \boldsymbol{\zeta}_{it}\]
    其中 \( \boldsymbol{\zeta}_{it} \sim N(\mathbf{0}, \mathbf{\Sigma}_{it}) \) 是白噪声。关键创新在于 \( \mathbf{\Phi}_{it} \) 和 \( \mathbf{\Sigma}_{it} \) 都是随时间变化的，并且它们本身也被建模为协变量 \( \mathbf{x}_{it} \) 的加性函数（例如，通过logit链接函数保证自回归系数在(-1,1)内，通过对数Cholesky分解保证协方差矩阵正定）。整个模型通过惩罚似然（penalized likelihood）进行估计，并使用广义交叉验证（GCV）或REML选择光滑参数。

  • 可观测数据：

    • 研究者实际能观测到的是：对于每个个体 \( i \)，在每个时间点 \( t \)，观测到 \( J \) 个结果变量 \( y_{ijt} \) 和一组协变量 \( \mathbf{x}_{it} \)（包括时间、干预指示变量等）。
    • 想要但观测不到的是：
      1. 潜在均值 \( \boldsymbol{\mu}_{it} \)：这是由固定效应和光滑函数决定的趋势项，是模型要估计的对象。
      2. 随机误差 \( \boldsymbol{\epsilon}_{it} \)：这是从观测值中减去均值后的残差，是模型要建模的序列依赖部分。
      3. 时变参数 \( \mathbf{\Phi}_{it} \) 和 \( \mathbf{\Sigma}_{it} \)：这些是描述误差过程动态结构的潜在参数，是模型的核心估计目标。
      4. 光滑函数 \( f_{jk}(\cdot) \)：这些是描述非线性趋势的未知函数，也是模型要估计的对象。

第二步：讲最小内核

本文的最小内核可以简化为一个个体、两个时间点、一个结果变量的情形，即单变量ITS，并假设序列依赖是时变的。这是Cho et al. (2021) 的核心内容，也是本文多元扩展的基础。

  • 最简特例：设 \( N=1, J=1, T=3 \)（干预发生在 \( t=2 \) 之后）。可观测数据为 \( y_1, y_2, y_3 \)。协变量 \( x_t \) 包括时间 \( t \) 和一个干预指示变量 \( I(t \ge 3) \)。模型退化为：

    \[y_t = \mu_t + \epsilon_t\]
    \[\mu_t = \beta_0 + \beta_1 t + \beta_2 I(t \ge 3) + f(t)\]
    其中 \( f(t) \) 是一个光滑函数，用于捕捉非线性趋势。误差项 \( \epsilon_t \) 服从一个时变AR(1) 过程：
    \[\epsilon_t = \phi_t \epsilon_{t-1} + \zeta_t, \quad \zeta_t \sim N(0, \sigma^2_t)\]
    这里 \( \phi_t \) 和 \( \sigma^2_t \) 都是随时间变化的。例如，\( \phi_t \) 可以建模为 \( \text{logit}(\phi_t) = \gamma_0 + \gamma_1 I(t \ge 3) \)，这意味着干预可能改变了序列依赖的强度。

  • 核心思路：在这个最简例子中，论文要解决的根本问题是：如何同时估计干预引起的水平变化（\( \beta_2 \)）、非线性趋势（\( f(t) \)）和时变序列依赖（\( \phi_t, \sigma^2_t \)）？

    • 为什么难？ 因为这三者是纠缠在一起的。如果错误地假设 \( \phi_t \) 是常数（即 \( \phi_t = \phi \)），那么干预后残差序列依赖性的变化可能会被错误地归因于干预效应（\( \beta_2 \)）或趋势变化（\( f(t) \)），导致估计偏倚。例如，如果干预后数据变得更具自相关性（\( \phi_t \) 增大），而模型假设自相关不变，那么干预后的残差看起来会“更持久”，这可能会被模型误认为是干预引起了水平或趋势的持续变化。
    • 本文的关键想法：通过同时对均值结构（\( \mu_t \)）和方差-协方差结构（\( \phi_t, \sigma^2_t \)）进行加性建模，并使用惩罚似然进行联合估计，从而在统计上分离这三种变化。具体来说，模型将 \( \mu_t, \phi_t, \sigma^2_t \) 都表示为协变量（包括时间、干预指示变量）的加性函数，并通过最大化一个包含光滑惩罚项的似然函数来估计所有参数。这样，干预对均值的影响（\( \beta_2 \)）和对序列依赖的影响（\( \gamma_1 \)）就可以被同时识别和估计。

三、这篇论文做了什么

三句话

  1. 研究了什么问题：针对心力衰竭患者自我监测的心率（HR）和血压（BP）数据，提出一个多元动态加性模型（MDAM），以同时分析药物变化引起的多元中断时间序列中的水平变化、非线性趋势变化和时变序列依赖。
  2. 核心工具/方法：将单变量动态加性模型（Cho et al., 2021）扩展到多元情形，使用时变向量自回归（TV-VAR(1)） 模型来刻画误差项的序列依赖和同期相关，并将所有参数（均值、自回归系数、误差协方差）都建模为协变量的加性光滑函数，通过惩罚似然进行估计。
  3. 主要结论：模拟研究表明模型参数恢复良好；忽略时变序列依赖会导致对水平变化和趋势变化的估计产生偏倚，且标准误被高估。真实数据分析展示了该模型如何揭示药物变化对HR和BP的复杂动态影响。

关键设定与假设

在第二节最小记号的基础上，补全完整设定：

  • 设定：

    • 数据为面板数据（panel data），有 \( N \) 个个体，每个个体有 \( T \) 个时间点。
    • 干预（药物变化）发生在已知的时间点，将整个时间序列划分为多个阶段（phases）。
    • 每个阶段内，趋势可以是线性的或非线性的，由光滑函数 \( f_{jk}(\cdot) \) 建模。
    • 误差项 \( \boldsymbol{\epsilon}_{it} \) 服从一个一阶时变向量自回归过程 TV-VAR(1)。作者选择一阶是因为在应用中发现一阶模型足以捕捉序列依赖，且计算上更可行。
    • 时变参数 \( \mathbf{\Phi}_{it} \) 和 \( \mathbf{\Sigma}_{it} \) 被建模为协变量 \( \mathbf{x}_{it} \) 的加性函数。例如，自回归系数 \( \phi_{ijk,t} \) 通过一个logit链接函数与加性预测器相连，以保证其值在(-1, 1)内。误差协方差矩阵 \( \mathbf{\Sigma}_{it} \) 通过对数Cholesky分解参数化，以保证其正定性。

  • 假设：

    • SUTVA（稳定单位处理值假设）：隐含假设，即一个个体的结果不受其他个体干预的影响。这在自我监测数据中通常是合理的。
    • 无同期干预：假设在干预发生的时间点，没有其他同时发生的事件影响结果。这是ITS分析的核心识别假设，但本文未明确讨论其合理性。
    • 模型正确设定：假设加性预测器、光滑函数、TV-VAR(1)结构以及链接函数的选择都是正确的。
    • 光滑性假设：假设光滑函数 \( f_{jk}(\cdot) \) 具有足够的平滑性（如二阶导数有界），以保证惩罚似然估计的一致性。
    • 弱平稳性：虽然参数是时变的，但假设在给定协变量 \( \mathbf{x}_{it} \) 的条件下，误差过程是局部平稳的，即TV-VAR(1)过程的特征多项式根在单位圆外。

  • 相比已有文献的强化/放宽：

    • 强化：相比Linden (2015) 等分段回归方法，本文放宽了趋势线性的假设，允许非线性趋势。
    • 强化：相比Cho et al. (2021) 的单变量模型，本文扩展到多元情形，可以处理结果变量之间的同期和滞后关系。
    • 强化：相比大多数ITS方法，本文明确建模了时变序列依赖，而不仅仅是将其视为需要校正的 nuisance。

主要结果

本文的主要结果来自模拟研究和真实数据分析，没有理论定理。

  • 模拟研究：

    • 设计：模拟了与真实数据相似的场景（2个结果变量，多个阶段，时变序列依赖）。设置了不同的条件，包括样本量（个体数N和时间点T）、效应大小、序列依赖强度等。
    • 核心量化结论：
      1. 参数恢复良好：MDAM对水平变化、趋势变化和时变序列依赖参数的估计具有可接受的准确性和精度（偏差小，均方根误差RMSE小，95%覆盖概率接近名义水平）。
      2. 忽略时变序列依赖的后果：当数据中存在时变序列依赖，但使用一个忽略时变性的模型（即假设 \( \mathbf{\Phi}_{it} \) 和 \( \mathbf{\Sigma}_{it} \) 为常数）进行分析时，会导致：
        • 水平变化和趋势变化的估计产生偏倚（bias）。
        • 标准误被高估（standard errors are inflated），即估计的标准误比实际抽样变异性更大，导致置信区间过宽，统计检验力下降。
    • 与Baseline对比：本文的baseline是忽略时变序列依赖的模型。模拟结果明确展示了MDAM相对于这个更简单、更常用的baseline的优势。

  • 真实例子：

    • 数据：来自一项心衰患者自我监测研究。患者每天测量并记录自己的心率和血压。数据中包含药物变化（如增加或减少利尿剂剂量）的时间点。
    • 如何应用：将MDAM应用于每个患者的数据，以估计药物变化对HR和BP的即时（水平变化）和长期（趋势变化）影响，同时控制时间趋势和时变序列依赖。
    • 结果：论文展示了几个典型患者的分析结果。例如，对于一位患者，增加利尿剂后，其收缩压（SBP）出现了立即的下降（负的水平变化），并且下降趋势持续（负的斜率变化）。同时，模型还估计出序列依赖（自相关）在干预后发生了变化。
    • 这个例子想说明什么：1) 展示MDAM在实际数据中的可行性和可解释性。2) 说明时变序列依赖是真实存在的，并且忽略它可能会改变对干预效应的推断。3) 提供一个可视化工具，展示模型如何分解出趋势、干预效应和残差动态。

证明路线与技术技巧（本文为应用型，无严格数学证明）

本文没有传统意义上的定理证明。其“证明”是通过模拟研究来验证方法的有效性。技术路线如下：

  1. 模型构建：将单变量动态加性模型扩展到多元，定义TV-VAR(1)误差结构，并参数化时变参数。
  2. 估计方法：使用惩罚极大似然估计（penalized maximum likelihood estimation, PMLE）。似然函数基于多元正态分布。光滑惩罚项（如对样条二阶导数的惩罚）被加到似然函数上，以控制模型复杂度。
  3. 计算算法：使用后向拟合算法（backfitting algorithm）或混合模型表示法（mixed model representation）进行估计。后者将光滑函数视为随机效应，从而可以使用REML进行估计，这通常更稳定。

  4. 模拟验证：生成大量模拟数据集，在这些数据集上应用MDAM和baseline模型，比较它们的参数恢复性能。通过重复模拟，计算偏差、RMSE、覆盖概率等指标，以“证明”MDAM在特定条件下是有效的，而忽略时变性的模型是有偏的。

  5. 关键跳跃点：从单变量到多元的扩展，主要技术难点在于TV-VAR(1)模型的参数化和计算。如何保证时变自回归系数矩阵 \( \mathbf{\Phi}_{it} \) 的稳定性（特征根在单位圆内）和时变协方差矩阵 \( \mathbf{\Sigma}_{it} \) 的正定性，是核心挑战。作者通过对数Cholesky分解和logit链接函数解决了这个问题，但这使得似然函数变得高度非线性，增加了计算难度。

  6. 技术技巧点名：
    • 对数Cholesky分解：用于参数化 \( \mathbf{\Sigma}_{it} \)，保证其正定性。
    • logit链接函数：用于参数化 \( \mathbf{\Phi}_{it} \) 的对角线元素，保证自回归系数在(-1,1)内。
    • 惩罚似然/REML：用于估计光滑函数，平衡拟合优度与模型复杂度。
    • 后向拟合算法：一种用于估计加性模型的迭代算法。

🔎 结论是否比证明窄

是的。论文的结论（“忽略时变序列依赖会导致偏倚”）是基于特定模拟条件得出的。这些条件（如样本量、效应大小、序列依赖的变化模式）被设定为与真实应用相似。作者没有从理论上证明“在任何情况下忽略时变序列依赖都会导致偏倚”，也没有给出偏倚大小的解析表达式。因此，结论的泛化能力是有限的。论文的结论应被理解为：“在类似于本应用的数据条件下，忽略时变序列依赖是一个有问题的做法。”

四、开放问题（点到为止，扎根具体语句）

  1. 因果识别的严谨性：本文完全从统计建模角度出发，没有讨论其估计量的因果解释。一个开放问题是：在什么识别假设下，MDAM估计出的“水平变化”和“趋势变化”可以被解释为药物变化的因果效应？ 这需要引入潜在结果框架，并讨论无混淆假设、无同期干预等条件在自我监测数据中是否合理。扎根于论文引言中未引用Shadish et al. (2002) 等因果推断文献这一事实。

  2. 模型选择与诊断：本文假设TV-VAR(1)模型是合适的。一个开放问题是：如何检验TV-VAR(1)的阶数是否足够？如何诊断模型是否错误设定（例如，是否存在未被建模的非线性交叉依赖）？ 扎根于论文中“我们选择一阶VAR是因为它在我们的应用中是足够的”这一陈述（原文：We chose a first-order VAR because it was sufficient for our application），这暗示了模型选择问题。

  3. 高维与计算可扩展性：当结果变量 \( J \) 很大时（例如，来自可穿戴设备的数十个生理指标），TV-VAR(1)模型的参数数量会呈 \( O(J^2) \) 增长，导致计算和估计困难。一个开放问题是：如何将MDAM扩展到高维多元时间序列？ 是否可以通过引入稀疏性假设（如对 \( \mathbf{\Phi}_{it} \) 施加Lasso惩罚）或使用低秩近似来解决？扎根于论文中仅处理了 \( J=2 \) 的情形。

  4. 异质性处理效应：本文模型假设所有个体的干预效应（水平变化、趋势变化）是相同的（固定效应）。一个开放问题是：如何将MDAM扩展为允许个体间存在异质性处理效应的混合效应模型？ 例如，可以假设 \( \beta_2 \)（水平变化）在不同个体间服从一个分布。扎根于论文中“我们假设固定效应”这一隐含设定。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论