“在现实世界应用中解决广义线性模型中的自相关问题”
解决广义线性模型中自相关问题的实际应用技巧
深入研究数据科学家最常见的噩梦之一
导言
线性回归中最大的问题之一是自相关残差。在这个背景下，本文重新审视线性回归，深入探讨科克伦-奥可特程序作为解决这个问题的方法，并探讨了在fMRI脑激活分析中的真实应用。
广义线性模型（GLM）再访
线性回归可能是任何数据科学家最重要的工具之一。然而，在时间序列的背景下，经常会看到许多误解。因此，让我们花一些时间重新审视这个概念。GLM在时间序列分析中的主要目标是对变量在一系列时间点上的关系建模。其中Y是目标数据，X是特征数据，B和A是要估计的系数，Ɛ是高斯误差。
指标表示数据的时间演变。以更紧凑的形式书写：
作者提供。
参数的估计是通过最小二乘法（OLS）进行的，该方法假设观测值与模型预测值之间的误差或残差是独立且具有相同分布（i.i.d）。
这意味着残差必须是非自相关的，以确保正确估计系数、模型的有效性以及预测的准确性。
自相关
自相关是指时间序列内观测之间的相关性。我们可以理解为每个数据点与序列中滞后数据点的关系。
自相关函数（ACF）用于检测自相关性。这些方法测量数据点与它们滞后值（t = 1,2,…,40）之间的相关性，揭示数据点与前面或后面的值相关的程度。自相关函数图（图1）显示了不同滞后的相关系数，指示自相关的强度和阴影区域上的统计显著性。
如果某些滞后的系数与零明显不同，则表明存在自相关。
残差中的自相关
残差中的自相关表明时间序列中当前和过去的误差之间存在关系或依赖。这种相关模式表明误差不是随机的，可能受到模型中未考虑的因素的影响。例如，自相关可能导致参数估计有偏，特别是在方差方面，影响变量之间关系的理解。这导致从模型中得出无效推论，对变量之间关系的结论具有误导性。此外，还会导致预测效果不佳，这意味着模型未捕捉到正确的信息。
Cochrane–Orcutt过程
Cochrane–Orcutt过程是一种经济计量学中著名的方法，也在多个领域中广泛应用，用于处理时间序列中的自相关问题，通过线性模型对误差项的串行相关性进行建模 [1,2]。我们已经知道这违反了普通最小二乘（OLS）回归的一个假设，即假设误差（残差）是不相关的 [1]。在文章的后面，我们将使用这个过程来消除自相关性，并检查系数的偏倚程度。
Cochrane–Orcutt过程的步骤如下：
- 1. 初始OLS回归：使用普通最小二乘（OLS）进行初始回归分析，估计模型参数。
- 2. 残差计算：计算初始回归的残差。
- 3. 检验自相关性：通过自相关函数图或Durbin-Watson检验等方法检验残差中是否存在自相关性。如果自相关性不显著，无需继续进行以下步骤。
- 4. 变换：通过对因变量和自变量进行差分变换来消除自相关性。这里的思路是使残差更接近无相关性。
- 5. 对变换后的模型进行回归：使用变换后的模型进行新的回归分析，并计算新的残差。
- 6. 检验自相关性：再次检验新的残差中的自相关性。如果仍存在自相关性，则返回第4步，并进一步变换模型，直到残差不再显示显著的自相关性。
最终模型估计：一旦残差不再显示显著的自相关性，使用Cochrane-Orcutt过程得出的最终模型和系数进行推论和得出结论！
现实世界应用：功能磁共振成像（fMRI）分析
功能磁共振成像（fMRI）简介
功能磁共振成像（fMRI）是一种神经成像技术，通过检测血液流量的变化来测量和绘制脑活动。它依赖于神经活动与增加的血流和氧合有关的原理。在fMRI中，当一个脑区变得活跃时，会触发血液动力学响应，导致血氧水平相关（BOLD）信号的变化。 fMRI数据通常由表示不同时间点脑活动的3D图像组成，因此大脑的每个体素都有自己的时间序列（图2）。
常用线性模型（GLM）
GLM模型假设测量到的fMRI信号是不同因素（特征）的线性组合，例如任务信息与称为血氧反应功能（HRF）的神经活动预期响应的混合。为简单起见，我们将忽略HRF的性质，只假设它是一个重要的特征。
为了理解任务对结果BOLD信号（因变量）的影响，我们将使用GLM。这意味着通过与任务信息相关的统计显著系数来检查影响。因此，X1和X2（自变量）是参与者执行任务的信息，通过数据收集与HRF（图3）进行卷积。
应用于真实数据
为了检查这个真实世界的应用，我们将使用由João Sato教授在巴西圣保罗联邦大学收集的数据，这些数据可以在GitHub上获得。独立变量fmri_data包含来自一个体素（单个时间序列）的数据，但我们可以对大脑中的每个体素进行分析。包含任务信息的因变量为cong和incong。这些变量的解释超出了本文的范围。
#读取数据fmri_img = nib.load('/Users/rodrigo/Medium/GLM_Orcutt/Stroop.nii')cong = np.loadtxt('/Users/rodrigo/Medium/GLM_Orcutt/congruent.txt')incong = np.loadtxt('/Users/rodrigo/Medium/GLM_Orcutt/incongruent.txt')#获取每个体素的时间序列fmri_data = fmri_img.get_fdata()#HRF函数HRF = glover(.5)#任务数据与HRF的卷积conv_cong = np.convolve(cong.ravel(), HRF.ravel(), mode='same')conv_incong = np.convolve(incong.ravel(), HRF.ravel(), mode='same')可视化任务信息变量（特征）。
拟合GLM
使用最小二乘法来拟合模型并估计模型参数，我们得到
import statsmodels.api as sm#选择一个体素（时间序列）y = fmri_data[20,30,30]x = np.array([conv_incong, conv_cong]).T#为预测变量添加常数x = sm.add_constant(x)#拟合线性回归模型model = sm.OLS(y,x).fit()#查看模型摘要print(model.summary())params = model.params
可以看到，当P > |t|小于0.05时，系数X1是统计显著的，这可能意味着任务确实对BOLD信号有影响。但在使用这些参数进行推断之前，检查残差（即y减去预测值）在任何滞后阶段都不是自相关的是非常重要的。否则，我们的估计将具有偏差。
检查残差的自相关
如前所述，ACF图是检查时间序列的自相关性的好方法。
从ACF绘图中可以看到在滞后1处存在高自相关。因此，该线性模型存在偏差，解决这个问题非常重要。
Cochrane-Orcutt解决残差的自相关
Cochrane-Orcutt 程序被广泛应用于fMRI数据分析中，以解决这类问题[2]。在这个特定案例中，滞后1的残差自相关是显著的，因此我们可以使用Cochrane–Orcutt公式来处理自回归项AR(1)。
# LAG 0yt = y[2:180]# LAG 1yt1 = y[1:179]# 计算滞后1的相关系数rho= np.corrcoef(yt,yt1)[0,1]# Cochrane-Orcutt公式Y2= yt - rho*yt1X2 = x[2:180,1:] - rho*x[1:179,1:]拟合变换后的模型
对修正后的Cochrane-Orcutt进行再次拟合。
import statsmodels.api as sm#添加常数到预测变量X2 = sm.add_constant(X2)#拟合线性回归模型model = sm.OLS(Y2,X2).fit()#查看模型摘要print(model.summary())params = model.params
现在，系数X1不再有统计意义，否定了任务对BOLD信号的影响的假设。参数标准误差估计有了显著的变化，表明残差自相关对估计的影响很大。
再次检查自相关
这是有意义的，因为有可能证明当存在自相关时，方差总是会有偏差[1]。
现在残差的自相关已经被消除，估计不再有偏差。如果我们忽略了残差中的自相关，我们可能会认为系数是显著的。然而，经过消除自相关后，发现该参数并不显著，从而避免了对任务与信号关联的错误推断。
结论
在通用线性模型的残差中存在自相关可能导致估计偏差、预测效率低和无效的推断。在实际的fMRI数据分析中，应用Cochrane–Orcutt程序能够有效地消除残差的自相关，避免错误结论，确保模型参数的可靠性和推断分析的准确性。
备注
Cochrane-Orcutt法只是解决残差自相关性问题的一种方法。然而，还有其他方法来解决这个问题，比如Hildreth-Lu过程和一阶差分法[1]。
致谢
该项目灵感来自João Ricardo Sato教授。
本文的笔记本可以在这里找到：链接。
参考资料
[1] Applied Regression Modeling, Iain Pardoe. Wileyl. 2023. https://online.stat.psu.edu/stat462/node/189/
[2] Sato JR, Takahashi DY, Cardoso EF, Martin Mda G, Amaro Júnior E, Morettin PA. 在功能连接性识别中应用干预模型于FMRI. Int J Biomed Imaging. 2006;2006:27483. doi:10.1155/IJBI/2006/27483
