检测心理过程动态中突然变化：一种基于状态空间建模和时间可变参数的方法

《British Journal of Mathematical and Statistical Psychology》：Detecting Critical Change in Dynamics Through Outlier Detection with Time-Varying Parameters

本文综述了一种基于状态空间建模（SSM）框架的创新方法，用于检测心理过程动态中的突然变化或“转折点”。作者将时间可变参数（TVP）建模与异常值检测程序相结合，提出了一种计算高效、可形式化检验的解决方案，适用于识别个体特异性或群体水平的动态转变。该方法通过扩展卡尔曼滤波器和固定区间平滑器进行模型拟合和参数估计，并利用广义最小二乘（GLS）导出的t统计量来筛查创新性异常值（innovative outliers），从而定位动态参数的突变点。文章通过模拟研究和真实面部肌电（EMG）数据应用，验证了该方法在单主体、多指标以及群体模型中的有效性和灵活性，为研究非平稳心理过程（如情绪诱导期间的变化）提供了强有力的分析工具。

　　

引言：动态转折点的关键重要性

在心理学和行为科学中，许多核心构念，如情绪、认知和动机，本质上是动态的，其潜在的过程参数（例如，自回归强度或长期设定点）可能会随时间发生变化。识别这些动态过程中的突然转变或“转折点”具有至关重要的意义，因为它有助于阐明关键变化时期（如干预前后）的差异、调查影响因素并实施针对性干预。然而，仅凭视觉检查是不充分的，需要能够提供正式检验的分析技术。

先前的研究，如Albers和Bringmann（2020）的工作，应用了带有通过广义可加模型（GAM）框架中的惩罚样条函数建模的时变系数的自回归（AR）模型，并使用Akaike信息准则（AIC）和贝叶斯信息准则（BIC）来筛选个体时变系数出现异常偏差的转折点。然而，随着时间序列长度的增加或过程维度的增加，这种逐点识别可能转折点的方法在计算上效率低下。

一种更有效的识别这些潜在转折点的方法是在所选模型的背景下进行诊断性筛查，寻找异常值。沿着这一思路，de Jong和Penzer（1998）提出了一种用于检测和区分加性异常值（additive outliers）和创新性异常值（innovative outliers）的方法。加性异常值是指在某一个时间点显示出异常值或模式，且仅影响该时间点，可以理解为极大的测量误差。相比之下，创新性异常值是指标志着一种非典型变化模式开始的数据点，这种模式的影响会持续超过一个时间点。

本研究提出利用这种异常值检测框架，通过单次模型拟合来有效检测TVP中的转折点或其他异常偏差。该方法将TVP参数化地作为额外的潜变量纳入动态模型，并在状态空间建模（SSM）框架下使用滤波-平滑方法估计其轨迹。先前的研究已经开发了利用这种滤波-平滑过程的副产品得出的统计量，可用于检测过程和数据中异常值的存在。然而，之前的实现和评估仅限于涉及潜变量线性函数的动态模型。将TVP作为潜变量纳入，会导致系统潜变量与这些TVP之间存在乘法关系的任何模型都是非线性的。因此，需要在异常值检测过程中进行一些调整以适应非线性，如估计部分所述。

在本工作中，我们提出了将You等人（2020）评估的异常值检测方法进行扩展，使其适用于对潜变量可微的线性或非线性动态函数。我们考虑了一个具有表现出离散跳跃的TVP的自回归模型作为动机示例。通过模拟研究，我们检验了该方法在三种情景下是否能有效检测TVP的突然变化：（1）潜在的个体过程；（2）具有个体特异性转折点的潜在群体过程；（3）具有多指标和因子载荷结构的潜在个体过程。

动机模型：时变自回归模型

滞后1的自回归模型（AR(1)模型）及其多元扩展已广泛应用于心理学研究。因此，在本文中，我们以其作为TVP模型的基础。传统的AR(1)模型通常写作：x_t = ω + β x__t-1+ ξ_t，其中x代表感兴趣的过程，x_t是其在时间t的值，ω与过程的整体均值（即截距）相关，ξ_t是随机噪声项。直观上，AR(1)模型表明一个变量的值随时间变化，是其前一个值、某个整体均值和一个随机扰动的函数。

为了更容易解释参数，我们采用了一种略微不同的参数化方法，明确映射到具有吸引子的自回归过程：x_t = μ + β (x__t-1- μ) + ξ_t, ξ_t ~ N(0, σ_ξ²)。其中μ代表过程围绕波动的长期设定点或吸引子，β是自回归强度或惯性（即倾向于偏离设定点的趋势），ξ_t是服从正态分布的随机过程扰动。μ可以取任何实数值，而β通常只在(-1, 1)范围内考虑。该范围外的值会导致过程在足够时间后爆炸性地趋向正负无穷，而在边界值1或-1则会导致在实数轴上的随机游走或围绕设定点的嘈杂振荡。过程噪声ξ_t假设服从均值为零、方差为σ_ξ²的正态分布。

如果我们假设过程是潜在的，只能通过一些指标或显变量观察到，那么我们还可以有一个测量模型。最简单的情况是加上测量误差，y_t = x_t + ε_t，其中ε_t也假设服从正态分布N(0, σ_ε²)。当潜在过程x由多个指标测量时，也可以采用因子模型：y_t = Λ x_t + ε_t，其中y_t是观测指标的向量，Λ是因子载荷矩阵。测量误差ε_t现在服从多元正态分布N(0, Σ_ε)。

添加TVP，顾名思义，是假设某些动态参数（而不是在整个过程或观察时间间隔内保持恒定）拥有自己值得研究的轨迹，因此可以用函数形式表示。由于每个参数代表变化模式的一个方面，因此本身也是一个未直接观察到的潜在构念，我们也可以将其视为一个潜在系统变量。

例如，我们可以有一个带有时变AR参数的模型：x_t = μ + β_t (x_t-1- μ) + ξ_t, β_t = f_β(β_t-1, t)。函数f_β可以采用任何函数形式，连续的或突变的，动态的或静态的。这里我们以β作为TVP的示例，但其他参数，如μ，也可以被指定为TVP，这取决于模型的可识别性，即测量值包含足够的信息来唯一识别SSM中设定的所有潜变量。

当目标是在没有先验知识的情况下探索TVP的形状时，可以为f_β选择一个灵活的模型。这里我们考虑将TVP表示为随机游走过程，这已显示出恢复各种数据生成TVP轨迹的潜力，并且当σ_β²= 0时简化为时间不变参数的特殊情况。在模型拟合中，原始过程变量x_t和附加变量β_t在这个例子中被组合成一个新的向量x_t^*，并作为双变量SSM同时进行拟合。方程显示了包含TVP的双变量SSM的潜在过程。该方程也显示了离散时间多元非线性SSM的潜在过程的一般形式。在离散时间SSM中，时间t-1的潜在状态（x__t-1^*）通过某个（可能非线性的）函数f加上噪声ξ_t映射到时间t的潜在状态（x_t^*）。潜在状态然后通过测量模型映射到观测变量。方程中的噪声是影响潜变量并随时间向前传递的过程噪声，而测量模型中的噪声是仅影响观测值且不随时间向前传递的测量噪声。

对于本例的TVP模型，过程噪声具有以下结构：[ξ_xt, ξ_βt]^T ~ N(0, Σ_ξ), Σ_ξ = diag(σ_x², σ_β²)。因此，过程噪声项被假设为服从均值为零、具有对角协方差结构的正态分布。类似地，对于本例的TVP模型，测量模型由给出。在测量模型中，只有潜变量x_t具有非零的因子载荷。TVP仅通过其对x_t的影响间接影响观测值。

状态空间模型估计与异常值检测程序

由于我们将时变参数作为潜变量纳入，我们可以将突然变化视为创新性异常值，因为时变参数的突然变化的影响会通过过程持续下去。为了检测此类创新性异常值，我们利用了拟合SSM时常见的滤波和平滑过程的现成副产品。

当拟合带有测量模型的SSM时，我们面临双重估计问题：即我们需要同时估计过程的潜状态和支配该模型的参数。在我们的方法中，第一部分通过扩展卡尔曼滤波器和固定区间平滑器完成，第二部分通过数值优化对数似然函数完成。

扩展卡尔曼滤波器与对数似然值优化用于参数估计

卡尔曼滤波器最初在信号处理中开发，是一种递归算法，旨在给定噪声观测的情况下估计动态过程的潜在轨迹。它及其各种扩展已广泛应用于状态空间建模框架中。它以一种顺序的、无记忆的方式从第一个到最后一个定义的时间点（通常是观测发生的时间点）运行，并在每个时间点估计潜状态，结合基于之前时间点累积信息的预测和当前时间点观测数据的新信息。传统的卡尔曼滤波器假设动态转移是线性的，过程噪声是正态分布的，因此滤波器可以利用多元高斯分布的清晰统计特性。然而，TVP的存在可能导致原本线性的系统出现非线性动力学（例如，如果TVP在自回归参数β上，通过交互作用产生非线性），因此我们求助于扩展卡尔曼滤波器来处理弱非线性。扩展卡尔曼滤波器使用泰勒级数展开在每个时间点附近对现在的非线性转移进行局部线性化。

大多数滤波算法涉及两个步骤：预测步骤，使用直到时间t-1的所有信息估计时间t的潜状态；更新步骤，使用时间t的观测值细化估计。为了在单个受试者的时间序列上执行扩展卡尔曼滤波器，这两个步骤如下：

步骤1：预测（使用上一次迭代更新步骤的结果）

x_t|t-1^*= E(x_t^*| y_0, ..., y_t-1) = f(x^_t-1|t-1^*)

B_t = ?f(x^_t-1|t-1^*) / ?x_t^*= J_f(x^_t-1|t-1^*)

P_t|t-1 = Cov(x_t^*| y_0, ..., y_t-1) = B_t P_t-1|t-1B_t^T+ Σ_ξ

步骤2：更新（使用预测步骤的结果和观测值y_t）

x^t|t^*= x^t|t-1^*+ K_t v_t

P_t|t = P_t|t-1 - K_t Λ P_t|t-1

K_t = P_t|t-1 Λ^T V_t^-1

v_t = y_t - E(y_t | x^t|t-1^*) = y_t - Λ x^t|t-1^*

V_t = Cov(v_t) = Λ P_t|t-1 Λ^T + Σ_ε

方程中的预测基于SSM将最近的潜状态估计x^_t-1|t-1^*映射到下一个时间点。基于直到时间t-1的所有可用信息，时间t的状态期望值记为x^t|t-1^*。除了将潜状态向前映射到下一个时间点，卡尔曼滤波器还将潜状态的协方差向前映射到下一个时间点。这个预测协方差记为P_t|t-1。为了进行这种映射，扩展卡尔曼滤波器将非线性动态模型f局部线性化为状态转移矩阵B_t，方法是在x^*的最新更新处求f的雅可比矩阵。预测协方差P_t|t-1考虑了（i）从前一个时间点更新的估计的方差（P__t-1|t-1），（ii）使用B_t将该更新方差线性映射到下一个时间点的方式，以及（iii）来自Σ_ξ的模型隐含的过程方差。

在给定直到时间t的所有信息对时间t的潜状态进行预测后，这些预测根据时间t的观测数据更新。这个更新步骤取决于卡尔曼增益K_t、预测误差v_t及其协方差V_t。卡尔曼增益可以看作是新的观测值（受测量误差协方差Σ_ε影响）与当前状态估计（其准确性由P_t|t-1捕捉）之间的相对权重。观测值噪声越大，卡尔曼增益值越低，因此在更新状态预测时赋予新观测值的权重越小。此外，时间t的预测误差v_t影响潜状态更新的方向和幅度。

v_t和V_t都用于通过预测误差分解形式的对数似然函数进行参数优化：log ?(θ) = -1/2 Σ_{j=1}^T [log(2π) + log|V_t_j| + v_t_j^T V_t_j^{-1} v_t_j]，其中θ是模型中所有时间不变参数的集合。这种形式的对数似然函数被用作数值优化的目标函数，以推导θ的最大似然估计。

卡尔曼固定区间平滑器

在获得最大似然估计θ^后，最后一次运行扩展卡尔曼滤波器，然后运行固定区间平滑器，在给定时间不变参数估计的情况下估计潜轨迹。与滤波器使用直到每个时间点t的信息不同，固定区间平滑器以反向递归方式运行，从最后一个时间点开始移动到第一个时间点，并进一步提高了卡尔曼滤波器得到的潜状态估计的准确性。我们根据de Jong和Penzer的定义，在固定区间平滑器中定义了两个新项：u_t和r_t。令p和q分别表示向量x_t^*和y_t的维度。u_t（维度q×1）是除了时间t的观测值外，所有T个观测值未解释的预测误差，而r_t（维度p×1）是基于所有未来预测误差对潜变量的向后调整。可以将r_t视为在时间t累积的扰动。N_t量化了r_t周围的不确定性。

对于每个时间点t，潜向量的平滑估计x^*计算如下：

x_t|T^*= E(x_t^*| y_0, ..., y_T) = x_t|t-1^*+ P_t|t-1 r__t-1

P_t|T = P_t|t-1 - P_t|t-1 N__t-1P_t|t-1

u_t = V_t^-1v_t - K_t^T r_t

r__t-1= Λ^T u_t + B_t^T r_t

L_t = B_t - K_t Λ

N__t-1= Λ^T V_t^-1Λ + L_t^T N_t L_t

卡尔曼滤波器递归地向前工作，创建给定直到时间t的信息的最佳潜状态估计，而固定区间平滑器递归地向后工作，创建给定直到最后一个时间点T的所有信息的进一步改进的潜状态估计。来自卡尔曼滤波器和固定区间平滑器的输出对于检测SSM中的异常值的程序至关重要。

异常值诊断度量

如前所述，时间序列中的异常值可以根据其时间影响分为两种类型。加性异常值是仅影响它们发生的那个时间点，而不影响任何后续时间点的异常值。如果我们用y_kt^*表示观测变量k在时间点t的真实测量值，那么时间t的加性异常值可以表示为y_kt = y_kt^*+ a_kt，其中a_kt是异常值的大小。在心理学研究中，加性异常值可以概念化为异常大的测量误差、数据输入错误或极其短暂的外部影响。

相比之下，创新性异常值代表对底层数据生成过程的冲击，不仅影响它们发生的时间点，还通过过程内的时间依赖性产生的传播影响影响后续时间点。也就是说，如果我们用x_t^#表示没有冲击时时间t的“未污染”过程，那么创新性异常值可以表示为x_t = x_t^#+ w_t，其中w_t是冲击的大小。

在估计潜轨迹时，运行固定区间平滑器得到的副产品可用于推导加性异常值和创新性异常值的检验统计量。在每个时间点t可以选择卡方统计量：ρ_inn,t^*2= r_t^T N_t^{-1} r_t, ρ_add,t^*2= v_t^T V_t^{-1} v_t，其中ρ_inn,t^*2~ χ2(p) 且 ρ_add,t^*2~ χ2(q)。卡方统计量基于广义最小二乘法（GLS），表示与预测误差（创新）和向后调整相关的平方马氏距离。它们检验时间t在观测时间序列的任何维度中不存在异常值的原假设（对于ρ_add,t^*2）或潜在过程（对于ρ_inn,t^*2）。

在许多情况下，也有兴趣识别在每个时间点，哪个特定的变量（潜在的或观测的）可能受到异常值的影响。为此，我们还可以为x_t^*和y_t中的每个变量h分别推导出t统计量，用于创新性异常值和加性异常值。这些t统计量的推导也遵循GLS框架。我们首先为创新性异常值定义一个冲击设计矩阵W_t，为加性异常值定义A_t。这些冲击设计矩阵用于指定要测试哪些变量，其元素为1（被测试）或0（否则）。例如，如果我们有兴趣筛查所有观测和潜变量的加性和创新性异常值，那么W_t和A_t的所有对角元素应设置为1，其他元素设置为0。如果我们只对筛查第一个潜变量的创新性异常值感兴趣，那么只有W_t的第一个对角元素设置为1，W和A中的所有其他元素设置为零。

使用来自平滑器的估计，我们可以推导出以下量：

s_t = A_t^T u_t + W_t^T r_t

Q_t = W_t - K_t A_t

S_t = A_t^T V_t^{-1} A_t + Q_t^T N_t Q_t

其中s_t被称为冲击对比（shock contrast），是平滑过程中调整后的创新和向后调整的组合，捕捉了可归因于加性或创新性异常值的偏差，S_t是其相应的协方差。

对于潜变量或观测变量中的每个维度h（取决于W_t和A_t在对角线索引h处是否为1），可以推导出单变量t统计量：t_ht = s_t(h) / √(S_t(h, h))，其原假设是时间t在维度h上不存在异常值。每个t检验的自由度对于创新性异常值是T - p，对于加性异常值是T - q。

这些理论结果确保了t统计量作为评估潜变量或观测分量中个体水平变化存在的有效检验统计量，并明确了分布假设。由于这些检验是针对每个潜变量维度分别定义的，它们可以直接应用于TVP，TVP现在被作为潜变量处理。这允许在相同的推理框架内检测TVP的突然变化。此外，基于GLS的推导允许这些量直接从单次运行卡尔曼滤波器和平滑器中获得，而无需在备选设定下重新估计模型。在检测创新性异常值方面，t检验在统计功效和I类错误率方面均优于卡方检验。因此，我们专注于使用t统计量来检测TVP中的突然变化。

所提出的SSM方法中的关键变化点检测可以很容易地扩展到几种不同的模型和数据组合。首先，因为检测创新性异常值的过程嵌入在SSM拟合过程中，我们可以很容易地适应具有多元线性测量结构和潜在弱非线性测量结构的模型。其次，这种异常值检测方法最初是为单主体状态空间模型开发的，已被Chow等人扩展到群体背景。在这种拟合SSM到多单元数据的情况下，即使我们拟合一个同质的群体模型，我们仍然能够检测个体特异性的变化点，因为滤波和平滑（以及随之而来的异常值检验统计量）是针对每个单元单独进行的。我们目前的工作通过使该方法适应弱非线性系统，扩展了这条发展路线。

我们通过三个模拟研究证明了该程序的效用和t统计量。模型拟合以及检验统计量的计算是在R中使用dynr包完成的。相关代码已上传至通讯作者的GitHub页面。

模拟研究

本节进行的模拟研究旨在评估我们提出的滤波-平滑方法在各种情景下的性能。我们首先在单主体AR(1)模型（不含任何测量误差，GAM方法未纳入测量误差）下，与先前建立的GAM方法在不同参数设置下进行比较，展示了所提出方法的性能。然后，我们将模拟扩展到三种GAM方法无法处理的情景：（1）具有一个观测变量/指标和测量误差的单主体AR(1)模型；（2）具有一个指标、测量误差和个体特异性变化点的多主体群体AR(1)模型；（3）具有多个指标和因子载荷结构的单主体AR(1)模型。在所有情景中，我们运行了200次蒙特卡洛重复，时间序列总共有T=100个时间点。由于该方法依赖于平滑算法，我们将“成功检测”定义为在真实起始点±T/10范围内，t统计量显示至少有一个显著创新性异常值。相反，我们将“错误检测率”定义为在“成功检测”窗口之外具有统计显著性t统计量的时间点数量除以此类时间点的总数。

模拟研究0：与GAM方法的比较

在比较中，我们研究了突然变化的四种情景，设定点μ和AR参数β各两种：大变化和小变化。突然变化的起始点在观测时间框架的中间。跳跃的值设置为与Albers和Bringmann工作中的值相似。对于μ的条件，大变化意味着变化前后μ的差异是过程噪声标准差的两倍，小变化意味着差异等于标准差。对于β的条件，大变化是从0.1到0.8（差异0.7），小变化是从0.5到0.8（差异0.3）。我们推测AR强度可能会影响滤波和优化过程在估计μ轨迹时的性能，因此我们将设定点条件扩展到包括一系列不同的真实时间不变β值。

我们假设我们对TVP的函数形式没有先验知识，因此我们将TVP设置为遵循离散时间随机游走，以便为算法提供探索参数空间的一些灵活性。为了启动滤波算法，TVP的初始值μ_0是自由估计的，其初始方差固定为零。系统变量x_0的初始值设置为与μ_0相同，其初始方差设置为σ_x²。模型中的所有其他参数都是自由估计的。虽然真实模型中没有测量误差，但测量误差标准差被固定在一个非常小的值，以确保滤波和平滑算法在数值上稳定，并且不会被大测量误差的存在所掩盖。

GAM方法遵循Albers和Bringmann的设定，将每个观测值x_t对其滞后一阶值x__t-1进行回归，将滞后值作为模型中的预测变量以考虑时间依赖性。为了使GAM方法与异常值检测方法具有可比性，我们只为具有真实变化的参数设置样条和变化点，并让另一个参数被估计为时间不变的。也就是说，如果我们期望μ中存在异常值，我们将建立一个以时间作为平滑项的截距的基线模型，然后通过添加一个虚拟编码预测变量（在时间t之前为0，之后为1）为每个时间点t构建一个候选变化点模型，从而在时间t将时间序列分成两块。然后我们将每个候选模型与没有任何变化点的基线模型进行比较。如果候选模型在拟合上（基于AIC或BIC值）比基线模型有改进，则时间点t成为候选变化点。在所有候选变化点中，改进最大的那个被选为第一个检测到的变化点。然后时间序列被这个变化点分成两个序列，并在这两个序列上分别重复相同的迭代过程。该过程重复进行，直到没有候选变化点被识别出来。AIC和BIC的有意义模型改进的阈值被设置为相对宽松的值-5，这表明没有变化点的零模型“相当不可能”是最佳模型。GAM方法中的模型拟合使用R包mgcv进行，使用薄板回归样条和平滑项中的10个基函数。回归样条惩罚的平滑参数使用mgcv包中的限制性最大似然选项进行估计。相关代码已上传至通讯作者的GitHub页面。

检测率如表1所示。我们可以看到，在大多数μ和β的情景中，使用状态空间建模方法观察到了更高的检测率。然而，当真实AR值较大时，状态空间建模方法在检测μ的跳跃方面确实遇到了困难，