张春成
2022/10/11阅读：51主题：默认主题
这可能是一种 BIAS
这可能是一种 BIAS
今天与同学讨论，聊到这么一篇文章
Putative rhythms in attentional switching can be explained by aperiodic temporal structure
很多研究发现，在注意力转移过程中，采集到的大脑信号具有“节律性”变化。
而这篇文章试图说明，这种现象可能是由于所采用的置换检验的统计方法破坏了原始信号的时间非周期的连续性而产生的。换句话说，它可能是时间置换检验方法带来的 BIAS，而不是信号本身所具有的属性。
我无意讨论后面的问题，只是单纯地好奇
从纯数据的角度出发，这个 BIAS 是如此严重的吗？
信号生成
实验需要一些数据，但我没有数据，所以我用 perlin 噪声生成一组信号。生成这组信号的参数是相同的，只是随机种子不同，因此可以认为这些信号是“独立同分布的”。同样的信号我生成了 500 次。
这些信号的时域波形如下，另外，我还对它们进行快速傅里叶变换，它们的频谱图也如下。
时域波形图
.png)
.png)
频谱图
.png)
.png)
import scipy
import numpy as np
import pandas as pd
import plotly.express as px
from tqdm.auto import tqdm
from noise import pnoise1
repeat = 500
num = 1000
data = []
for _ in tqdm(range(repeat)):
seed = np.random.randint(500)
octaves = 5
persistence = 0.8
s = np.array([
pnoise1(e, persistence=persistence, base=seed, octaves=octaves)
for e in np.linspace(0, 1, num)
])
s -= np.mean(s)
data.append(s)
data = np.array(data)
print(data.shape)
df = pd.DataFrame()
for j in range(5):
df['signal-{}'.format(j)] = data[j]
px.scatter(df, title='Signals').show()
px.imshow(data, title='Signal').show()
统计检验
经验分布的检验
下面进行统计检验，由于同样的信号我生成了 500 次，因此就相当于得到了 500 个独立同分布样本。我用这些样本当作经验分布，根据经验分布就可以得到某个信号的“概率值”。
当然，这里暗含着一个假设，那就是经验分布的形式是已知的。为了方便起见，在经验数据的支持下，我认为它服从均值和方差已知的正态分布。接下来就可以进行假设检验，
假设某个信号不是来自经验分布
下面就是这个假设为真的概率，即 p 值。可以看到，p 值沿频率轴的分布具有挺奇怪的模式，这表示不同频率下信号的分布不尽相同。并且，这个模式并不是持续不变的，当被测信号不同时，它的 p 值分布曲线也不同。然而，这些 p 值几乎都在 0.05 以上，这表示我们没有足够的证据推翻原假设，也就是说，我们只能相信“它来自经验分布”，并且这个结论几乎在每个频率点都是成立的。
随机置换的检验
接下来，使用时间置换检验方法来进行同样的操作。与上述方法不同的是，它不再需要经验分布，而只依赖于被测信号自身。具体的做法是将被测信号沿时间进行随机打乱，将打乱后的信号作为经验分布的近似，从而实现检验。
它也能够得到 p 值曲线，但这条曲线看上去没有那么多波动，并且，它有一个挺奇怪的问题，就是在低频范围内，它“总是”能够得到显著的差异，即 p 值小于 0.05。事实上，这个值还可以更小甚至相当小。
如果之前“经验分布的检验”结果是正确的话，那么这种随机转换的检验方法显然是有点 BIAS 了。
这可能是由于通过打乱时间点的方法进行的置换检验，将信号退化到了与白噪声相似的功率谱。在此基础上进行的统计检验，无异于检验某一段频率的能量是否显著地大于其他频段。这种检验很难得到他曾经宣称的结论，即某种特殊条件下的脑信号，具有特殊的节律成分。
repeat = 500
for j in tqdm(range(5), 'Permutation test'):
d = data[j]
s = d.copy()
f_list = []
for _ in range(repeat):
np.random.shuffle(s)
f = np.abs(np.fft.fftshift(np.fft.fft(s)))
f_list.append(f)
fs = np.array(f_list)
std = np.std(fs, axis=0)
mean = np.mean(fs, axis=0)
f = np.abs(np.fft.fftshift(np.fft.fft(d)))
z = (f-mean) / std
p = 1 - scipy.special.ndtr(z)
fig = px.scatter(p, title='p value (signal {})'.format(j))
fig.add_hline(y=0.05)
fig.show()
作者介绍