纺锤波：EEG中纺锤波参数分析和检测框架，并应用于睡眠纺锤波
文章来源于微信公众号（茗创科技），欢迎有兴趣的朋友搜索关注。
导读
目的：EEG纺锤波是一种窄带振荡信号，是研究被试状态和神经功能的生物标志物。现有的纺锤波检测方法大多数通过优化专家标签的一致性来选择算法参数。本文提出了一种基于纺锤波属性稳定性选择算法参数的新框架，并阐明了这些属性对几种算法参数选择的依赖性。
方法：为了演示这种方法，本文开发了一种新算法(Spindler)，该算法使用与Gabor原子匹配的追踪来分解信号，并计算参数值的精细网格中每个点的纺锤波。在将特征曲面作为参数的函数进行计算后，Spindler根据特征曲面几何形状的稳定性选择算法参数。
结果：相对于几种常见的有监督和无监督EEG睡眠纺锤波检测方法，Spindler的性能较好。Spindler是一个开源的MATLAB工具箱(https://github.com/VisLab/EEG-Spindles)。除了Spindler之外，该工具箱还提供了几种其他纺锤波检测算法的实现，以及将地面真值与预测相匹配的标准化方法和理解算法参数表面的框架。
意义：本研究表明，基于物理约束而不是标记数据的参数选择可以提供有效的、全自动的、无监督的纺锤波检测。这项研究还揭示了在不考虑纺锤波属性对参数的依赖情况下应用交叉验证的风险。为优化一种性能指标或匹配方法而选择的参数未针对其他参数进行优化。此外，阐明预测指标在算法参数选择方面的稳定性对于这些算法的实际应用至关重要。
1.前言
纺锤波是在一个或多个头皮区域主导EEG（脑电图）信号的窄带活动突发。睡眠纺锤波是11至17Hz频率范围内的振荡活动，是NREM（非快速眼动）第2阶段睡眠的主要指标。睡眠纺锤波具有定义明确的结构，通常前面有一个称为K复合波的信号特征，它由一个短的负偏差，紧接着是一个大的正偏差，然后是一个大的负偏差组成。每分钟纺锤的数量（速率或密度）、纺锤振幅和纺锤持续时间是个体内部的稳定指标，但在个体之间的差异很大。纺锤波指标通常显示不同睡眠周期内和夜间睡眠过程中的特定变化模式。睡眠纺锤波被认为在记忆巩固中发挥作用，一些研究表明睡眠纺锤波密度和IQ得分呈正相关。睡眠纺锤波模式的异常是阿尔茨海默病等几种神经退行性疾病的生物标志物，对睡眠纺锤波的观察可用于多种睡眠障碍的诊断。
手动评估睡眠阶段、睡眠纺锤波和其他睡眠特征的技术已被美国睡眠医学学会标准化。文献中已经提出了数十种睡眠纺锤波检测的自动化方法，其中许多方法涉及使用带通滤波进行频谱分解、傅里叶变换、小波分解或匹配追踪(MP)。最近提出的一些算法应用变换来分离振荡和瞬态信号，以更有效地近似纺锤波。
几项著名的基准测试工作比较了用不同方法对标准化数据集和指标的睡眠纺锤评级。Warby等人描述了对专家、非专家和睡眠纺锤波评估自动化方法的性能进行基准测试的大量工作。他们使用24位专家对2000个25秒周期的睡眠记录（来自110名被试）的共识评级，制定了一个黄金标准评级。他们发现单个专家评级和黄金标准之间一致性的平均F1为0.75。该研究量化了评分者和算法之间的差异，突显了创建一致标准的难度。他们还评估了几种自动化算法，发现F1与黄金标准的一致性范围为0.21到0.52之间。
不幸的是，来自Warby研究的数据无法提供给其他研究人员用于基准测试。许多研究人员使用两个公开可用的档案：梦境睡眠纺锤波数据集（Dreams Sleep Spindle Database）和蒙特利尔睡眠研究档案(MASS) SS2数据集来比较算法。每个数据集中都有两位专家对纺锤波的评级。对于每个数据集，两位专家的评分者间一致性相对较低，许多研究人员使用联合的专家评级作为评估的地面真值。
当前的纺锤波检测算法分为两类：预选算法参数的无监督算法和通过优化与专家评级相关的性能来设置参数的监督算法。相比之下，本文提出了一种混合方法，该方法明确计算纺锤波属性对算法参数的依赖性，并使用生成的纺锤波属性曲面的几何形状来选择一组“稳定”的算法参数。该方法(Spindler)及其评估框架可在MATLAB工具箱(https://github.com/VisLab/EEG-Spindles.git)中实现。如果提供了地面真值信息，Spindler将会为许多常见指标生成性能表面。为了便于比较，该工具箱还包括其他几种纺锤波检测算法。在可能的情况下，该工具箱使用作者提供的带有包装器的代码来合并常见的处理方法。
本文的结构如下。第2节阐述了Spindler算法，定义了五种用于将检测到的纺锤波与地面真值匹配的常用方法，并介绍了用于评估的性能指标。第2节还描述了用于评估的测试数据和用于比较的其他算法。第3节比较了算法性能并显示了纺锤波属性如何依赖于各种算法参数。第4节讨论了Spindler方法的优点和缺点，以及所提出的方法与其他计算纺锤波方法的关系。第5节提供了结论性意见。
2.方法
2.1 Spindler算法
Spindler算法包括三个步骤：在参数网格上进行信号分解和重建，在参数网格上评估纺锤波属性，以及根据这些纺锤波属性曲面的几何形状选择“最佳”参数。如下所述，纺锤波特性很大程度上取决于两个参数：Ns（用于表示信号的Gabor原子/秒数）和 Tb（应用于重建信号的阈值，以检测纺锤波的存在）。Spindler的其他参数可以改变，但不会强烈影响结果。本文的所有结果都使用表1中总结的默认Spindler参数设置。
表1.Spindler参数
2.1.1 使用Gabor原子的匹配追踪(MP)来表示信号
Spindler将原始EEG数据的单个通道作为输入。本文报告的结果使用的是通道Cz。纺锤波带通滤波在信号分解之前使用EEGLAB的pop_eegfiltnew对纺锤波频率范围（本文为11至17 Hz）内的信号进行滤波。受纺锤几何学的启发，Spindler将信号s(t)近似为Gabor原子的加权和：
Gabor原子g由四个参数表示：
Spindler使用表1中指定参数的原子组成的MP字典。N，表示信号的原子总数，是Ns乘以信号长度（以秒为单位）。f是以Hz为单位的原子振荡频率（对于睡眠纺锤波来说，频率在11到17 Hz范围内，步长为0.5 Hz），σ是原子高斯包络宽度的一半（0.0625、0.125或0.25），φ是相对于高斯峰值的正弦振荡相位（假设为0）。时间偏移τ对应于原子在每个时间点的中心位置。
Spindler采用高效的贪婪策略将候选原子投影到信号上，并依次去除投影绝对值最大的原子的贡献。该算法在选择指定的最大原子数/秒(Ns)后停止，而不是在表征达到指定的重建误差时停止。这种贪婪策略确保了具有大量原子的表征是具有较少原子数的表征的超集。因此，MP表征只需计算一次即可在Ns的参数网格上生成信号。一个有效的表征只需保留添加到表征中的每个原子(a，f，σ，τ，φ)的值，该表征介于参数网格上的当前Ns值和下一个Ns值之间。可以快速地重新计算每个选定Ns的信号重建。默认情况下，Spindler使用[0，0.4]中的39个均匀间隔的Ns值。（对于本文检查过的所有示例，当Ns＞0.3时，纺锤速率存在差异）
应该注意的是，更传统的MP实现在整个频率范围内使用不受限制的Gabor原子字典，然后应用后续处理来分离所需的信号。这种方法的难点在于，当信号具有显著的伪影时，很难使用后续处理来消除这些伪影带来的影响。通过对信号进行带通滤波，然后将表征限制为仅可能表示纺锤波的原子，Spindler可以处理较大的噪声信号。
2.1.2. 参数曲面的生成
为了检测在(Ns，Tb)处重建信号中是否存在纺锤波，Spindler首先取重建信号的绝对值，然后除以这些值的第95个百分位，以减少信号异常值的影响。Spindler将第95个百分位以上的所有缩放值设置为1。这种缩放将Tb限制在区间[0，1]内，从而允许对所有数据集使用通用的参数值网格。Spindler形成一个掩码，在重建信号中标记绝对值超过指定阈值Tb的数据点。Spindler然后将最小纺锤波分离标准应用于掩码以连接相邻的点，并应用最小和最大长度标准以消除伪点。纺锤波被检测为掩码中的连续区域，从而形成参数网格上每个(Ns，Tb)的纺锤波事件列表。
2.1.3 参数选择的启发式方法
给定(Ns，Tb)网格上每个点的纺锤事件列表，Spindler会形成三个属性曲面：纺锤速率、纺锤时间段和纺锤长度。如图1所示，两个具有代表性的睡眠数据集的曲面几何形状在数据集之间是相似的，并被用作参数选择的基础。Ns绘制在水平轴上，每条细线对应于表1中Tb的一个固定值，范围从0（深蓝色）到1（深红色），极值Tb=0和Tb=1分别用虚线和蓝绿色线强调。粗实的蓝色线对应于各个参数的中心值，该值是通过对每个Ns值在Tb=1和Tb=0处的各自参数值进行平均来计算的。
图1.两组代表性数据集的纺锤波特性曲面。
Spindler启发式尝试选择算法参数，使参数值的微小变化不会导致纺锤波属性值的大变化。当Ns值很小时，所有纺锤波速率横截面重叠（图1的第一行），表明在此范围内，纺锤波速率与Tb无关，并且与Ns大致呈线性变化。在几何上，这是因为贪婪MP表示按顺序添加原子，因此当表征信号的原子数量有限时，无论阈值如何，该表征只捕获能量最大（相对长）的纺锤波。然而，随着Ns的增加，纺锤波速率对阈值会变得非常敏感。
当一个给定的阈值可以检测到所有可以检测到的纺锤波，并且向表示添加更多原子不会导致更多纺锤波时，就会发生纺锤率饱和。当这种情况发生时，可以观察到纺锤率曲线（图1的第一行）先下降，然后趋于平稳。纺锤率饱和首先出现在Tb=1（蓝色点状曲线）处，最后出现在其他Tb阈值。限制Ns小于第一个饱和点，并将其作为Ns的最小值，Tb=1的纺锤波速率曲线的斜率为零。这种情况意味着向表征中添加额外的原子不会产生更多的纺锤波。我们也要求Ns足够大，使纺锤率跨阈值的标准偏差大于0。当它为零时，纺锤率完全独立于阈值，这表明即使是最突出的纺锤波也没有足够的原子来表示。Ns的合格范围由图1中的粗灰色水平线标记。
为了获得关于阈值的稳定表示，Spindler通过对Tb=0和Tb=1的曲线求平均来计算中心纺锤波分数曲线。然后Spindler找到 Tb∗，其纺锤波分数曲线的阈值仍然在Ns的合格范围内最接近中心分数曲线。最后，Spindler选择Ns∗的值作为在Ns的合格范围内使Tb∗的主轴长度最小的值。正如结果部分所示，对于将检测到的纺锤波与地面真值相匹配的所有方法来说，(Ns，Tb)平面中没有一个点可以优化性能。Spindler方法的总体策略不是在与专家评级一致性的意义上寻找最优值，而是在参数网格上选择一个中心点，以便参数值的微小变化不会导致属性值的大变化。合格点的范围是专门选择的，以“远离边缘”。
2.2 将地面真值与检测到的纺锤波相匹配
大多数标准性能指标依赖于混淆矩阵（TP、FN、FP和TN）的生成，以匹配检测到的纺锤波与地面真值。当检测结果与地面真值相匹配时，就会出现真阳性(TP)。当地面真值显示了纺锤波，但检测结果显示没有检测到纺锤波时，就会出现假阴性(FN)；而误报(FP)表示，检测到的纺锤波不在真值中时。当检测结果和地面真值均未显示纺锤波时，则表示真阴性(TN)。不幸的是，没有标准的方法来计算纺锤波的混淆矩阵，文献中出现了许多不同的方法，每种方法都有其优缺点。确定两组纺锤波之间评级一致性的方法会显著影响性能结果。
本研究已经实现了五种确定一致性的方法，基于：epoch counts、simple hits、intersections、onsets和times。该工具箱报告了这五种一致性方法(Count，Hit，Intersect，Onset和Time)的结果，列出了一系列从混淆矩阵派生出来的性能指标。
真阴性(TN)的计算是有问题的，而且Time是唯一提供直接计算TN的匹配方法。本文在分析中强调不依赖于真阴性的指标，但为了完整起见，本研究使用Devuyst等人提出的默认TN方法，该方法假设所有纺锤波的长度都是固定的（睡眠纺锤波为1秒），并移除由信号长度的真阳性、假阳性和假阴性表示的时间。然后，默认方法是用剩余时间除以假定的纺锤波长度来计算真阴性。
2.2.1 Count方法
Count方法将在指定长度的段(本文为30秒)中检测到的纺锤波数制成表格。与多个段重叠的纺锤波仅计入第一个重叠间隔。每段中的真阳性(TP)数是检测到的纺锤波数和该区间内真纺锤波数的最小值。如果真正的纺锤波多于检测到的纺锤波，则差值就是该段中的假阴性(FN)数。相反，如果在段中检测到的纺锤波多于真正的纺锤波，则差值是该段中假阳性(FP)的数量。TN使用默认方法计算。整体混淆矩阵为每一段中计算的混淆矩阵的总和。
2.2.2 Hit方法
如果真正纺锤波和检测到的纺锤波之间存在任何重叠，则Hit方法检测TP。如果给定检测到的纺锤波与真正的纺锤波不重叠，则Hit检测FP，如果真正纺锤波不重叠所有检测到的纺锤波，则检测FN。Spindler查找真正纺锤波和检测到的纺锤波的互补间隔（开始和结束），并将互补间隔中的重叠计数为TN。然而，研究者提醒用户不要使用基于TN的性能指标，尤其是对于Hit方法。对于非TN指标，Hit匹配往往比其他方法更乐观地估计性能，而对于TN指标则不太乐观。此外，没有纺锤波的长周期被计为单个TN，而不是由间隔长度进行加权。Hit方法的另一个挑战是，将整个数据集标记为单个纺锤波（TP）得到的精度和召回值为1。Hit性能对参数选择相对不敏感，但如果有太多非常长的事件，则应谨慎使用。
Hit方法与Nonclercq等人使用的方法有一些相似之处。这些作者将信号分解为0.5秒的重叠窗口，间隔125毫秒，并根据窗口的频率幅度特性将每个窗口分类为纺锤波或非纺锤波。他们根据特定窗口是否具有与真正或检测到的纺锤波匹配的标签对窗口进行分类。Nonclercq算法将没有标签的窗口归类为真阴性。因为Nonclercq算法依赖于知道真阳性的频率-幅值的分布，所以它不适合一般实现。
2.2.3 Intersect方法
Intersect方法类似于Warby 等人提出的用于睡眠纺锤波评估的并集/交集规则算法。Intersect方法计算表示两个比较事件（真正和检测到）的时间间隔的交集和并集之间的比率。如果此比率大于相交公差（默认为0.2秒），则相应的纺锤波是潜在匹配项。 Intersect方法选择具有最大合格率的事件对作为匹配，将匹配归类为真阳性，并从进一步考虑中删除事件对。Intersect将不匹配的真实事件归类为假阴性，将不匹配的检测到的事件归类为假阳性。Warby等人没有提供用于计算真阴性(TN)的伪代码。本文使用默认方法来计算Intersect方法的TN。
2.2.4 Onset方法
Onset方法确定在指定的Onset公差范围内，真正和检测到的纺锤波的起始时间是否一致。当检测到的纺锤波在真正纺锤波开始的Onset公差内，就会出现真阳性。当检测到的纺锤波不在任何真正纺锤波的Onset公差内，就会出现假阳性。当真正的纺锤波在其Onset公差内没有检测到纺锤波时，就会出现假阴性。本研究中的Onset是使用默认方法来判断真阴性的。与Intersect方法一样，Onset会从进一步考虑中删除任何匹配的事件。该方法经常用于评估Onset公差为0.5s的睡眠纺锤波。 研究者发现性能结果对Onset公差极为敏感，0.5s给出的结果似乎过于乐观。Spindler默认使用0.3s作为Onset公差。
2.2.5 Time方法
Lawhern等人引入的用于评估alpha纺锤波的Time方法，该方法使用每个状态所花费的实际时间长度，而不是将匹配视为事件。考虑到难以确定纺锤波的准确开始时间，Time允许通过时间公差在任意方向上扩展事件以获得更好的匹配。默认情况下，Time使用0.2s作为时间公差。Time将检测到的纺锤波不与任何真纺锤波重叠的时段归类为假阳性时间，将真纺锤波不与任何检测到的纺锤波重叠的时段归类为假阴性时间。任何未标记为TN、FP或FN的周期都会累积为Null协议或真阴性时间。Time方法不区分单个重叠和多个重叠，而是根据真正纺锤波和检测到的纺锤波重叠的时间量来计算真阳性。
2.3 性能指标
最常见的性能指标来自混淆矩阵（TP、FP；FN、TN）。正如Warby等人所指出的，在这种情况下，真阴性(TN)的计算是不可靠的，因此不应使用基于TN的指标。本文强调的指标如表2所示，但Spindler基于混淆矩阵计算了一组全面的指标，并且可以将这些指标作为Ns和Tb的函数来绘制（图2）。
表2.一些不基于TN的性能指标。
图2.在(Ns，Tb)平面上，F1(左列)和FDR(右列)在5种不同匹配策略(Count；Hit；Intersect；Onset；Time)上的性能。
2.4 用于比较的算法
为了便于比较，Spindler工具箱还包括三个无监督睡眠纺锤波检测算法（CWT-a7、CWT-a8和SEM-a6）和两个监督睡眠纺锤波的检测算法（Spinky和McSleep）。此外，Spindler工具箱还包括有监督和无监督的alpha纺锤波算法，这里不讨论。监督算法将部分数据集的预测与专家评级相匹配，以选择参数，然后计算其余数据的性能。对于本文中的监督预测，研究者交替使用数据的前半部分和后半部分作为各自的训练集，然后对结果进行平均。Spindler也可以用作监督算法（表示为Spindler-S），在训练集上使用地面真值来确定最佳(Ns，Tb)。
2.4.1 连续小波概率估计（CWT-a7 和 CWT-a8）——无监督。
Tsanas和Clifford提出的连续小波概率估计方法，采用131个Morlet尺度和5.4-40.5Hz范围内的伪频率的连续Morlet小波变换对信号进行分解。24个小波基落在11-17Hz的睡眠纺锤波频率范围内。在每个采样时间点，该算法根据该时刻前十个小波系数中有多少在纺锤波频率范围内具有Morlet尺度来估计纺锤波的概率。该算法使用10个最近邻对瞬时概率估计进行平滑处理，并通过应用一些控制纺锤波概率和持续时间的规则将邻近点组合成纺锤波。一个规则示例：如果两个间隔的平均纺锤波概率＞0.7，每个间隔的长度至少为0.1s，且间隔小于0.3s，则将这两个间隔合并为一个纺锤波。Tsanas和Clifford提供了两个版本的算法（用CWTa7和CWT-a8表示），使用不同的方法来平滑瞬时概率函数。这些实现将频率范围作为参数，但对所有阈值、纺锤波持续时间标准和纺锤波分离标准进行硬编码。这些阈值和Morlet小波参数严重依赖于100Hz的信号采样率。Spindler工具箱包装器在执行前将信号重新采样到100Hz。从调用返回时，Spindler工具箱包装器会组合相邻事件、匹配地面真值并计算性能指标。
2.4.2. 信号包络建模 (SEM-a6)——无监督。
Wendt等人提出了一种睡眠纺锤波检测算法，该算法从两个EEG通道中获取输入：中央通道（如C3）和枕部通道（如O1）。该算法在11-16Hz频段进行滤波后，使用来自中央通道的信号包络的局部极值检测候选睡眠纺锤波。如果纺锤波频率为13Hz，并且在定义候选纺锤波的区间内，中央通道的带通功率小于枕骨通道中的功率，则该算法通过移除候选纺锤波来进一步校正“阿尔法效应”。该算法假设输入信号对应于NREM睡眠。本文使用Warby等人提供的算法。该算法本质上依赖于特定的频率滤波器，并且不易适用于检测其他类型的纺锤波。但Spindler工具箱为用户实现该算法提供了一个包装器，用于组合相邻事件、匹配地面真值和计算性能指标。
2.4.3 Spinky Q可调小波估计——监督
Spinky提供了一种睡眠纺锤波检测的监督方法，该方法使用可调Q因子小波变换分解信号。Q是一个硬编码常数，选择它的目的是使小波具有最小周期数，对应于所需频带11-17Hz在0.5秒间隔中的振荡次数。在小波分解之后，Spinky使用Count匹配方法以30秒的间隔确定一个阈值作为Sensitivity-FDR的最大值。Spinky根据信号计算此阈值的最小和最大允许值。Spindler使用这些最大值和最小值之间的值网格（默认间隔10个单位）进行参数分析。Spindler工具箱为纺锤波检测的原始Spinky MATLAB实现提供了一个包装器。请注意，原始Spinky工具箱中提供的K-complex检测和方便手动评级的python用户界面未包含在Spindler工具箱中。
2.4.4 McSleep——监督
McSleep及其单通道前身DETOKS是一种非线性方法，将EEG建模为稀疏瞬态、低频信号和稀疏振荡信号的总和。这些方法需要四个正则化参数（λ0、λ1、λ2、μ）、一个阈值（c）以及一些块大小和迭代大小参数，可以根据信号采样频率固定λ0和λ1。使用具有Count匹配的F1指标优化了其余参数。McSleep工具箱不提供参数优化代码。本研究的分析分别使用Dreams和MASS-SS2数据集中的λ0、λ1和μ参数。Spindler使用λ2和c的网格来生成McSleep的参数曲线。默认情况下，λ2的网格为[20，50]，间距为1，c的网格为[0.5，3]，间距为0.1。本文采用3个通道（中央额叶、中央、中央枕叶）进行分析。
2.5. 用于评估的数据
本文使用两个公开可用的数据集合来评估Spindler算法，这些数据集合是在知情同意的情况下获得的，并且不包含个人身份信息。
2.5.1梦境睡眠纺锤波数据库
本文从公开的梦境睡眠纺锤波数据库(Dreams Sleep Spindlesdatabase)中选取了30分钟的片段。第二位专家的评级仅包括纺锤波开始时间。根据Devuyst等人，本研究假设第二位专家的纺锤波长度为1秒。这些片段没有说明特定的睡眠阶段，也没有进行预处理。在每种情况下，专家都标注了CZ位置对应的中央通道。
2.5.2. 蒙特利尔睡眠档案数据库(MASS)
蒙特利尔睡眠档案队列1子集SS2包括19名健康被试的整晚记录，采样率为256Hz。本研究使用了这19名被试的记录，其中15名由两位专家评级，4名（被试4、8、15和16）仅由专家1评级。整晚的记录也同样有睡眠阶段的评级。本研究提取了每个记录中第2阶段睡眠的最长片段用于分析。
结果
3.1 Spindler作为一种无监督算法
图2显示了具有代表性的睡眠纺锤波数据集（MASS SS2被试1中最长的第2阶段）在两种不同指标和几种匹配方法上的Spindler性能曲线。从图2中可以明显看出几个特征。所有匹配方法的性能曲线具有相同的一般形状，但Onset方法的性能比其他方法差得多，并且对阈值的精确选择非常敏感。请注意，Spindler不一定选择指定性能指标的全局最优值，而是尝试选择一个相对稳定的值。性能结果显示出明显的权衡——对应于最高F1值的参数不会导致低FDR值。Count匹配方法通常比其他方法更乐观。Hit、Intersect和Time方法给出了更多的中心值，这些值在Ns的候选范围内保持相对接近并且相对独立于Tb（由图 2 中的粗灰色水平线标记）。另一个有趣的特征是性能依赖于地面真值的选择。专家1对纺锤波的评级比专家2更保守。因此，专家评级联合的表现与专家2的相似，专家1的FDR更高。
本研究使用图形和统计技术将Spindler性能与其他无监督算法进行比较。补充材料包含带有统计详细信息的表格(stacks.iop.org/JNE/15/066015/mmedia)。图3显示了Spindler使用五种不同的匹配方法（Count，Hit，Intersect，Onset，Time），并在三种常见的无监督算法（CWT-a7、CWT-a8和SEM-a6）的F1（左列）和FDR（右列）上进行性能比较。总体而言，与其他三种无监督算法相比，Spindler具有更高的F1和更低的FDR值。Spindler 与其他算法的配对t 检验证实，对于这两个指标，Spindler的性能明显优于MASS-SS2的其他三种算法，F1 的p值小于10-5，FDR 的p值小于10-14。对于Dreams数据集，Spindler性能与CWT-a7性能没有显著差异，但明显优于其他两种算法。使用其他匹配方法的结果类似。
图3.在Time匹配方法上，Spindler与其他无监督算法的F1(左列)和FDR(右列)性能指标的比较。
为了更好地理解使用这些算法的实际意义，研究者检查了MASS-SS2数据集的三个纺锤波属性：平均纺锤波速率（纺锤波/分钟）、纺锤波时间分数和平均纺锤波长度(s)。这些属性不是独立的，因为纺锤波时间的比例与平均纺锤波长度和平均纺锤波速度的乘积成正比。表3显示了由四种无监督算法以及专家评估者预测的MASS-SS2的平均属性。
表3.MASS-SS2非监督算法和专家评级的均值属性。
专家1的纺锤波率和纺锤波分数明显低于专家2。但是，专家1的纺锤波长度在统计上并不低于专家2，这可能是因为专家2在所选段上的纺锤波长度预测中显示出较大的标准偏差。专家2的纺锤波长度与Spindler、CWT-a7或SEM预测的长度没有显著差异。然而，CWT-a8纺锤波长度预测值明显高于所有其他值。Warby等人的黄金标准数据集专家评级的平均受试者睡眠纺锤波长度为0.74s。他们报告的SEMa6对单个被试平均纺锤波长度的平均值为 0.55s，这与本研究从Dreams数据集中获得的结果一致。
Spindler的纺锤波速率和纺锤波分数(p＜10−3) 明显小于CWT-a7和SEM-a6，但纺锤波长度预测在Spindler、CWT-a7和SEM-a6之间没有显著差异。CWT-8对所有三个属性的估计都显著更高（p＜10-2）。作为比较，Warby 等人报告的黄金标准数据收集的专家共识评级的平均被试睡眠纺锤波率是2.3纺锤波/分钟，一些被试的纺锤波率几乎为零，而其他人则有纺锤波速度在5到10 纺锤波/分钟范围内。这些作者还报告了SEM-a6的被试睡眠纺锤波率，黄金标准聚类平均值约为7个纺锤波/分钟，这与表3和图4中给出的结果一致。Spindler和专家2的平均纺锤率约为5 纺锤波/分钟。CWT-a8的平均纺锤率聚集在12个纺锤波/分钟左右，远远超出正常范围。图4证实了这些结果，该图比较了MASS-SS2 队列中被试的平均纺锤波特性。
图4
3.2. Spindler 作为监督算法 (Spindler-S)
监督算法使用“地面真值”标记的训练数据搜索参数值的最佳组合。此参数搜索还提供了一个机会来检查预测属性相对于参数选择的稳定性，如图5所示，用于最近在文献中报道的两种监督算法：Spinky和McSleep。Spinky（图5的左列）具有单个阈值参数(T)。Spinky 根据数据幅度和频率计算网格范围后使用固定的网格分辨率。图5显示Spinky的平均纺锤波长度相对平坦，但随着T增加超过50，纺锤波速率和纺锤波分数急剧下降。本研究中检查的所有数据集的Spinky参数曲线具有相似的总体形状。
图5
McSleep有大量参数，如方法部分所述。根据既往文献建议，研究者设定MASS-SS2集合的K=256 O=128、μ=0.5、Nit=40、λ0=0.6 和λ1=7。对于给定的集合，这些参数必须凭经验确定，研究者发现结果对λ0和λ1的选择有些敏感，但这里不探讨这项工作的这种依赖性。在(c，λ2)的二维网格上计算属性/性能表面。本文对c在[0.5，3.0]范围内以0.1为增量使用网格，对于λ2在[20，50]范围内以1为增量使用网格。McSleep纺锤波长度与参数选择相对独立，但特别是对于较小的c值，纺锤波速率和纺锤波分数随着λ2的增加而显著下降。
图6显示了F1指标对 Spinky（左列）和 McSleep（右列）参数的依赖性。颜色表示匹配方法。对于所有匹配方法，性能曲线的形状大致相同。图6的性能指标是在整个选定数据段上计算的，而不是将段拆分为训练集和测试集，因此结果与图2相当。根据不同的专家，Spinky和McSleep 具有不同形状的地面真值参数曲线，导致算法预测的纺锤波属性值不同。当专家2的评分可用时，专家评分结果的并集接近专家2的结果。McSleep结果还显示出对参数c的依赖性，这强烈影响其纺锤波率和纺锤波分数预测。
图6
表4使用两个指标（F1和FDR）、两个版本的地面真值（专家1和专家2）和Time匹配方法比较了MASS-SS2集合的不同监督和无监督方法的平均性能。作为参考，表4还显示了每个专家评分使用其他评分者作为真值的表现。表4表明，性能结果很大程度上受哪个专家评估者被认为是地面真值的影响。
表4.MASS-SS2的无监督和有监督纺锤波检测算法的均值属性(SD)。
如图5所示，纺锤波属性（尤其是速率和分数）受到Spinky和McSleep参数选择的强烈影响。这反过来又取决于哪个专家被用作真值。表5比较了两种专家选择的平均属性值。与Spindler或Spinky相比，无论用于确定真值的专家如何，McSleep给出的平均纺锤波长度明显更小。
表5.MASS-SS2的监督算法和专家评级的均值属性(SD)。
图7显示了MASS-SS2数据集按被试细分的结果。对于所有指标和算法，专家2的预测通常高于专家1的预测。有监督Spindler属性预测通常包含来自无监督版本的预测。算法之间的纺锤波分数有很好的一致性。有监督Spindler的长度预测似乎偏高。
图7
讨论
通常很难评估有多少已发布的纺锤波检测算法将在实践中执行。Warby等人的综合研究是了解纺锤波评估变异性来源的重要一步。如Warby等人的综合研究所示，地面真值从来都不是绝对的，专家评估者之间存在重大分歧。报告已知公开可用数据的结果对于重现结果至关重要。然而，除了报告检测到的事件的开始和结束时间之外，研究人员还应该记录较小的片段如何组合成较大的纺锤波，以及如何实现与地面真值的匹配，因为这些选择会极大地影响结果。Spindler工具箱通过为该评估提供通用基础架构，允许在方法之间进行更真实的比较。这项工作的一个关键结果是证明优化性能指标而不考虑预测属性的稳定性可能会导致不切实际的结果。在不考虑这些依赖关系的情况下报告交叉验证结果并不能真实反映性能。
睡眠纺锤波检测方法通常使用Intersect或Onset将检测到的纺锤波与真值相匹配。这两种方法都将匹配视为不同的，并从进一步考虑中删除匹配的事件。这些方法在计算真阴性时假定纺锤波长度是固定的，这是经验数据不支持的假设。研究者同意Warby等人的观点，即由于纺锤波和纺锤波之间的严重不平衡，不应使用依赖于真阴性的性能指标。
本研究专注于平衡精度和召回率的指标F1以及FDR。监督算法的“最佳”参数的选择很大程度上取决于所选择的指标以及其他因素。起始F度量通常低于同一数据集的所有其他度量，并且对起始容忍度的选择非常敏感。Time、Intersect和Hit对Spindler在稳定区域中的F1进行了相当接近的估计。对于Spinky和McSleep，Time比 Intersect稍微乐观一些。Hit指标容易受到不平衡的影响。如果一个算法声明一切都是一个纺锤波，那么精度和召回率及其派生指标可能接近1，而对于其他匹配算法，这些指标几乎为零。Count方法始终高估所有指标和算法的性能。总的来说，研究者认为应该使用 ntersect 或Time来进行性能评估。而且研究者更喜欢Time，因为它给出了TN的内在度量。本文报告的大多数结果都使用Time方法。
本研究检查过的所有算法都有隐藏参数和嵌入的常量值，这些参数是由ad hoc方法选择的。由于研究人员通常使用一小部分测试数据来测试他们的算法，因此临时选择可能不适用于更大的环境。研究者认为，算法充分暴露所有隐藏参数以进行有效评估非常重要。而且评估纺锤波特性对参数选择的依赖性很重要。研究者手动检查了为所有DREAMS和MASS-SS2数据集生成的Spindler参数曲线，并验证参数曲线与图1中显示的曲线相似。算法也可能对伪影敏感，研究者认为要么应该使用标准化的自动化方法删除伪影，要么应该根据未处理的数据评估算法。
文献中已经提出了其他几种使用Gabor原子或Gabor、Fourier 和Dirac原子组合的基于MP的方法。这些算法将信号分解为短窗口（0.5s、16s、20s或30s，具体取决于方法），使用完整的原子字典分解每个窗口中的信号，然后在分解中考虑具有频率的原子和适合纺锤波的尺度。CWT-a7和CWT-a8使用相关策略，使用来自指定字典的Morlet小波使用连续小波变换分解每个时间点的信号。然后，这些算法通过检查每个点的10个小波在纺锤波频率和比例范围内的部分来计算纺锤波的概率。
相比之下，Spindler是一种全局算法，它使用仅由捕获纺锤波的原子组成的字典来分解整个信号。然后Spindler构造一个全局参数曲面来设置特定于数据集的阈值。与大多数其他算法不同，除了纺锤波事件之外，Spindler还生成纺锤波信号的重建。这使得使用这些重建信号研究跨通道的纺锤波同步或频移成为可能。Spindler 还可以处理源信号，例如由独立成分分析产生的信号，以将纺锤波特征映射到网络。
5.结论
本研究对无监督睡眠纺锤波算法（Spindler、SEM-a6、CWT-a7、CWT-a8）的比较表明，Spindler在各种设置和使用未处理的EEG数据中提供了更好的性能。Spindler也可以应用于监督设置(Spindler-S)，并提供与Spinky和McSleep相当的预测性能。这项工作提供了一个文档齐全的工具箱，并且该工具箱实现了Spindler和本文中讨论的其他算法，以促进算法之间的进一步基准测试和比较。Spindler工具箱还包括生成可视化所需的功能。Spindler工具箱可在https://github.com/VisLab/EEG-Spindles上免费获得。
原文：Larocco, J. , Franaszczuk, P. J. , Kerick, S. , & Robbins, K. . (2018). Spindler: a framework for parametric analysis and detection of spindles in eeg with application to sleep spindles. Journal of neural engineering, 15(6), 066015.1-066015.15.
