近红外光谱脑功能成像(fNIRS):2.实验设计、指标计算与多重比较

一、实验设计的策略与方法

近红外光谱成像（INIRS）作为一种非侵入性脑功能成像技术，为研究大脑活动提供了一种高效、生态效度高的方法。然而，为了充分利用INIRS技术并确保实验结果的准确性和可靠性，研究者必须精心设计实验。这包括合理排布实验条件、选择恰当的刺激呈现方式、决策实验设计类型、精确设定组块时长和事件相关设计的刺激间隔、严格控制可能影响结果的混淆因素、以及合理抽样和分配受试者。

此外，研究者还需深入理解INIRS技术的特性，包括其在血液动力学响应、空间分辨率和信噪比方面的局限性，并在实验设计中采取相应措施以降低这些局限性对研究结果的影响。通过这些综合考虑和策略运用，INIRS实验设计能够为探索大脑功能和认知过程提供坚实基础。

实验条件的合理安排：在自变量和其水平确定后，研究者需要精心安排实验条件的呈现顺序和方式。这不仅包括不同条件的组合，也涉及到同一条件下多次试次的排布，目的是有效激发并测量大脑的神经响应。
刺激呈现方式的策略选择：INIRS实验通常采用两种基本的呈现策略：组块设计和事件相关设计。组块设计通过集中呈现同一类型的刺激，适合于研究需要长时间稳定状态的任务，如观察情绪状态的变化。而事件相关设计则通过随机顺序呈现不同条件的刺激，适合于探索大脑对特定事件的响应，例如在Oddball范式中检测对异常刺激的神经活动。
设计类型的决策依据：选择使用组块设计还是事件相关设计，取决于实验的具体目标和任务特性。如果实验的重点是检测特定条件下的神经活动是否存在显著差异，组块设计可能更为合适；若目标是分析单个事件的神经响应特性，则事件相关设计可能更有优势。
关键参数的精确控制：无论是组块设计还是事件相关设计，都涉及到关键参数的控制。对于组块设计，组块时长是一个关键因素，它影响着血氧响应的累积和恢复过程。事件相关设计中，刺激间隔（ISI）的设置则决定了血氧信号的分离度和可检测性。
混淆因素的有效控制：实验中不可避免地存在一些与自变量无关但可能影响结果的因素。研究者需要通过各种方法，如消除法、恒定法、平衡法或统计控制法，来控制这些混淆因素，以确保实验结果的纯净性。
受试者样本的代表性：受试者的选取和分配对实验结果的普适性至关重要。研究者应采用随机抽样等方法，确保样本能够代表总体，并通过被试间或被试内设计合理分配受试者，以减少个体差异对实验效应的影响。
INIRS技术特性的深入理解：尽管INIRS技术具有高生态效度和能够提供多维度的血红蛋白浓度变化信息，但其在空间分辨率和信噪比方面的局限性也不容忽视。研究者需要在实验设计中考虑到这些技术特性，采取相应措施来提高数据质量和分析的准确性。

二、个体血液动力学响应指标计算

当我们使用fNIRS（功能性近红外光谱成像）的技术来观察大脑活动时，我们想要知道每个人对不同刺激（比如声音、图像等）的反应有多强烈。为了做到这一点，我们需要先建立一个模型，这个模型会猜测fNIRS数据中包含了哪些部分，比如大脑在没有刺激时的正常活动、刺激引起的反应，以及一些干扰信号（比如噪音）。

从个体水平的INIRS（近红外光谱成像）数据中计算血液动力学响应指标，常常使用一般线性模型（GLM）进行数据建模和分析。

一般线性模型（GLM）：
- GLM用于建模INIRS数据，包括基线活动、刺激诱发的血液动力学响应和噪声等成分。
- 公式为 $y=XB+E$ ，其中 y 是观测信号，X 是设计矩阵，B 是模型参数，E 是残差。
- 使用最小二乘法估计模型参数 $\hat{B} = (X^T X)^{-1} X^T Y$ 。
组块设计的血液动力学响应指标计算：
使用哑变量编码实验条件，形成方波形状的任务参考波。接着，通过GLM估计模型参数，得到每种实验条件的血液动力学响应强度。
也就是说，在进行大脑研究时，我们设计了的实验由几个“组块”组成，每个组块都是一种特定的任务或刺激。为了分析大脑对这些任务的反应，我们使用了一种叫做“哑变量”的编码方式，这就像给每个任务贴上不同的标签。

通过这种编码，我们可以创建一个图形，这个图形看起来像一系列的方波，每个方波代表一个任务组块。这个图形就是我们的“任务参考波”，它帮助我们标记出实验中每个任务开始和结束的时间，类似于嵌入式时序中的时钟信号。
事件相关设计的血液动力学响应指标计算：
单次刺激的血液动力学响应通过与血液动力学响应函数（HRF）卷积得到。多个刺激的响应通过线性叠加得到。
通俗的说，我们设计的实验中包含了许多单独的、一次性的刺激（比如突然出现的声音或图像）。我们想要了解大脑对这些单次刺激的反应。

为了做到这一点，我们使用了一个叫做“血液动力学响应函数（HRF）”的工具，这个工具可以帮助我们预测大脑在接收到刺激后血液流动的变化模式。我们将每个单次刺激与这个HRF进行“卷积”，卷积是一种数学操作，可以让我们把每个刺激的反应模式叠加到一起，形成一个总的反应模式。

通过这种方式，我们可以得到每个单次刺激对大脑血液流动的影响，然后将这些影响“线性叠加”，也就是简单地把它们加在一起，来得到所有刺激对大脑的总体影响。这样，我们就能知道大脑对一系列单次刺激的总体反应有多强烈了。
设计矩阵中的其他解释变量：
包括基线、趋势项、噪声等，以提高模型的解释能力和减小残差。
个体实验效应提取：
使用对比向量从模型参数中提取感兴趣的个体效应量。对比向量 c 与模型参数向量 B 的乘积 cB表示特定的个体效应。

通过这些步骤，研究人员可以准确地从INIRS数据中提取出个体对不同类型刺激的血液动力学响应指标，为进一步的神经机制研究提供数据支持。

我们不妨写一段Python代码来模拟INIRS数据，用最小二乘法估计一般线性模型（GLM）的参数。然后，计算预测信号和残差，并绘制观测信号与预测信号以及残差的图形。最后，使用对比向量提取了感兴趣的个体效应量。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import linregress

np.random.seed(0)
n_samples = 100
n_features = 5
y = np.random.randn(n_samples)  # 观测信号
X = np.random.randn(n_samples, n_features)  # 设计矩阵
B = np.random.randn(n_features)  # 模型参数
E = np.random.randn(n_samples)  # 残差

# 使用最小二乘法估计模型参数
B_est, _, _, _ = np.linalg.lstsq(X, y, rcond=None)

y_pred = X.dot(B_est)  # 计算预测值

residuals = y - y_pred  # 计算残差

plt.figure(figsize=(12, 6))
plt.plot(y, label='Observed Signal')
plt.plot(y_pred, label='Predicted Signal')
plt.legend()
plt.title('Observed vs Predicted Signal')
plt.xlabel('Time')
plt.ylabel('Signal')
plt.show()

# 绘制残差
plt.figure(figsize=(12, 6))
plt.plot(residuals, label='Residuals')
plt.legend()
plt.title('Residuals')
plt.xlabel('Time')
plt.ylabel('Residual')
plt.show()

# 提取感兴趣的个体效应量
c = np.array([1, 0, 0, 0, 0])  # 对比向量，假设感兴趣的是第一个效应
cB = c.dot(B_est)

print(f'感兴趣的个体效应量: {cB}')

三、群体统计推断和多重比较

多重比较问题（multiple comparison）在多通道fNIRS研究中很重要。想象你有一篮子苹果，每个苹果代表一个测量通道，而你想要找出哪些苹果是坏的（即哪些通道显示出显著的实验效应）。如果你只检查一个苹果，并且设定一个规则说，如果这个苹果有超过5%的概率是坏的，你就认为它是坏的。这个规则在只检查一个苹果时是合理的，因为你有5%的误差允许范围。

但是，如果你现在要检查篮子里所有的苹果（即所有的测量通道），并且对每个苹果都使用同样的5%的规则，那么随着你检查的苹果数量增加，你最终找到至少一个坏苹果的概率就会大大增加。这就是多重比较问题：当你对多个通道进行检验时，即使所有的通道都没有真正的实验效应，你也很可能因为随机误差而错误地认为某些通道是显著的。

在fNIRS研究中，这意味着如果你对每个通道都使用固定的显著性阈值（比如5%），你最终可能会得到一些假阳性的结果，即错误地认为某些通道在实验中显示出显著效应，而实际上并没有。这就是为什么需要进行多重比较校正，以确保在检查多个通道时，假阳性的概率仍然控制在可接受的范围内。

多重比较校正方法包括FWE（family wise error）校正和FDR（false discovery rate）校正。FWE校正方法如Bonferroni校正和基于随机场理论的校正方法，而FDR校正方法如Benjamini-Hochberg校正方法。

1.FWE（family wise error）校正

FWE（family wise error）校正是一种多重比较校正方法，其目的是控制在一组统计检验中至少出现一个假阳性错误的概率。在fNIRS研究中，由于涉及到多个测量通道，每个通道都可能显示出显著的实验效应，因此需要进行FWE校正来避免假阳性结果的增加。

1.1 FWE校正的原理

FWE校正的基本思想是将所有的通道看作一个“家庭”，每个通道都是这个家庭中的一员。在进行统计检验时，我们设定一个零假设（H0），即“任何家庭成员都未激活”，换句话说就是任何一个通道都没有实验效应的存在。FWE校正的目标是在H0的假设条件下，控制任意（一个或者多个）测量通道出现假激活的概率（即FWE P），从而实现对单通道上的显著性阈值进行校正。

1.2 在fNIRS研究时的操作

设定FWE错误率：首先，研究者需要设定一个FWE错误率，通常为0.05，这意味着我们希望控制在所有通道中至少有一个假阳性错误的概率不超过5%。
计算校正后的显著性阈值：接下来，使用特定的FWE校正方法（如Bonferroni校正）来计算每个通道的校正后的显著性阈值。Bonferroni校正是最简单和直接的FWE校正方法，它假设各个观测通道出现假阳性错误是彼此独立的。校正后的显著性阈值（ $\alpha_{\text{corrected}}$ ）可以通过将原始显著性水平（α）除以通道总数（n）来计算，即 $\alpha_{\text{corrected}} = \frac{\alpha}{n}$ 。
应用校正后的显著性阈值：在实际的统计检验中，对每个通道的p值进行比较时，使用校正后的显著性阈值（ $\alpha_{\text{corrected}}$ ）而不是原始的显著性水平（α）。只有当一个通道的p值小于或等于 $\alpha_{\text{corrected}}$ 时，我们才认为该通道显示出显著的实验效应。
考虑通道间的相关性：在fNIRS研究中，相邻通道的数据可能存在相关性，这会影响FWE校正的效果。因此，实际操作中可能需要考虑通道间的相关性，并采用更为复杂的校正方法，如基于随机场理论的FWE校正方法。

通过这些步骤，FWE校正确保了在多通道fNIRS研究中，即使进行了多次统计检验，假阳性结果的概率仍然控制在预设的水平以下，从而提高了研究结果的可靠性。

2.FDR（false discovery rate）校正

FDR（False Discovery Rate）校正是一种多重比较校正方法，它旨在控制假阳性结果的比例，而不是像FWE校正那样控制至少出现一个假阳性错误的概率。FDR校正特别适用于大规模数据集，如fNIRS研究中的多通道数据，因为它在保持较高统计功效的同时，有效地控制了假阳性结果的比例。

2.1 FDR校正的原理

FDR校正的基本原理是控制假阳性结果在所有被拒绝的零假设中的比例。具体来说，FDR校正的目标是确保在所有被标记为显著的结果中，假阳性结果的比例不超过预设的FDR水平（通常为0.05或0.10）。

2.2 在fNIRS研究时的操作

设定FDR水平：首先，研究者需要设定一个FDR水平，通常为0.05或0.10，这意味着我们希望控制在所有被标记为显著的结果中，假阳性结果的比例不超过5%或10%。
计算校正后的p值：接下来，使用特定的FDR校正方法（如Benjamini-Hochberg (BH) 方法）来计算每个通道的校正后的p值。BH方法是一种常用的FDR校正方法，它按照原始p值从小到大排序，然后对每个p值进行校正。
确定显著性阈值：在BH方法中，首先将所有通道的原始p值从小到大排序，然后从最小的p值开始，逐步检查每个p值是否满足校正后的条件。具体来说，对于第i个通道，如果其原始p值小于或等于 $\left( \frac{i}{n} \right) \times Q$ ，其中n是通道总数，Q是设定的FDR水平，则认为该通道是显著的。
应用校正后的显著性阈值：在实际的统计检验中，对每个通道的p值进行比较时，使用校正后的显著性阈值。只有当一个通道的p值小于或等于校正后的显著性阈值时，我们才认为该通道显示出显著的实验效应。
考虑通道间的相关性：与FWE校正类似，在fNIRS研究中，相邻通道的数据可能存在相关性，这会影响FDR校正的效果。因此，实际操作中可能需要考虑通道间的相关性，并采用更为复杂的校正方法，如基于随机场理论的FDR校正方法。

通过这些步骤，FDR校正确保了在多通道fNIRS研究中，即使进行了多次统计检验，假阳性结果的比例仍然控制在预设的水平以下，从而提高了研究结果的可靠性。

下面是一段代码，代码模拟了fNIRS研究中的多个通道的p值，然后分别进行了FWE和FDR校正。使用statsmodels库中的multipletests函数进行校正，其中method='bonferroni'用于FWE校正，method='fdr_bh'用于FDR校正。最后，代码绘制了原始p值和校正后的p值，并在图中标记了FWE和FDR的显著性阈值。

首先得装库：

pip install statsmodels

详细代码如下：

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.stats.multitest import multipletests

# 模拟数据
np.random.seed(0)
n_channels = 50
p_values = np.random.rand(n_channels)  # 模拟p值

# FWE校正
alpha_fwe = 0.05
corrected_p_values_fwe = multipletests(p_values, alpha=alpha_fwe, method='bonferroni')[1]

# FDR校正
alpha_fdr = 0.05
corrected_p_values_fdr = multipletests(p_values, alpha=alpha_fdr, method='fdr_bh')[1]

# 绘制原始p值和校正后的p值
plt.figure(figsize=(14, 7))

# 原始p值
plt.subplot(1, 2, 1)
plt.scatter(range(n_channels), p_values, color='blue', label='Original p-values')
plt.axhline(y=alpha_fwe, color='red', linestyle='--', label='FWE threshold')
plt.axhline(y=alpha_fdr, color='green', linestyle='--', label='FDR threshold')
plt.title('Original p-values')
plt.xlabel('Channel')
plt.ylabel('p-value')
plt.legend()

# 校正后的p值
plt.subplot(1, 2, 2)
plt.scatter(range(n_channels), corrected_p_values_fwe, color='red', label='FWE corrected p-values')
plt.scatter(range(n_channels), corrected_p_values_fdr, color='green', label='FDR corrected p-values')
plt.axhline(y=alpha_fwe, color='red', linestyle='--', label='FWE threshold')
plt.axhline(y=alpha_fdr, color='green', linestyle='--', label='FDR threshold')
plt.title('Corrected p-values')
plt.xlabel('Channel')
plt.ylabel('Corrected p-value')
plt.legend()

plt.tight_layout()
plt.show()