清晨醒来时刻前后是心血管疾病(CVD)事件发生的高峰期,这一点很可能与夜间睡眠结束时交感神经活动的激增有关。本文以 70 位在两年随访期内发生了 CVD 事件和 70 位未发生 CVD 事件的 140 位受试者为研究对象,提出了一种两层模型方法以探究觉醒前的心率变异性(HRV)特征是否有利于两类受试者的区分。在该方法中,第一层采用极端梯度提升算法(XGBoost)构建分类器,通过评估该分类器的特征重要性实现特征筛选;筛选出来的特征作为第二层模型的输入来构建最终的分类器。在第二层模型中,比较了 XGBoost、随机森林和支持向量机三种机器学习算法,以确定基于何种算法建立的模型可以得到最优的分类效果。研究结果显示,基于觉醒前 HRV 特征构建的 XGBoost+XGBoost 模型性能最优,准确率高达 84.3%;在所使用的 HRV 特征中,非线性动力学指标在模型中的重要性优于传统的时域、频域分析指标;其中尺度 1 下的排列熵、尺度 3 下的样本熵较为重要。本研究结果对 CVD 的预防、诊断以及 CVD 风险评估系统的设计有着积极的参考价值。
引言
心血管疾病(cardiovascular disease,CVD)是由心脏和血管异常引起的一系列疾病的总称。目前,我国 CVD 患病率处于持续上升阶段,而因 CVD 致死则占城乡居民总死亡原因的首位[1],因此 CVD 的早期预测对人类的健康和社会发展有着积极的意义[2]。
睡眠,作为一种复杂的生理过程,与机体的各基本系统密切相关,其中就包括心血管系统[3]。心血管系统的稳定性受到众多因素的控制,其中自主神经系统对心率的调节起着至关重要的作用。自主神经系统包括交感神经和副交感神经,交感神经兴奋时,心率会增加;副交感神经兴奋时,心率会降低。自主神经系统活动失衡会增加许多疾病(包括 CVD)的发病率和死亡率[4]。研究表明,在一天的不同时刻,急性 CVD 事件的发生率是不同的,清晨醒来时刻前后是 CVD 事件发生的高峰期[5-9]。虽然导致清晨 CVD 事件增加的原因并未明确,但不少研究认为这一状况可能与夜间睡眠末期交感神经活动的激增有关[10-14]。
自主神经系统使得心血管系统能够对机体内外的易变刺激做出迅速反应,导致逐次心跳的心率呈现出细小的波动,通常被称为心率变异性(heart rate variability,HRV)。HRV 分析方法主要包括时域分析法、频域分析法和非线性动力学分析法[15]。作为一种非侵入、无创的手段,HRV 分析已经广泛用于评估自主神经系统功能和心血管系统健康状况[16-18]。在文献中,也有利用 HRV 分析进行 CVD 风险评估的报道[19-21],但这些研究存在着一些不足。首先,使用的 HRV 数据并非针对睡眠状态。而在本文的前期研究中发现,睡眠 HRV 分析能够有效预测短时的 CVD 事件[22]。其次,采用的分析方法大都基于时域、频域分析,鲜有涉及非线性动力学分析。然而大量的研究表明,心脏动力系统富含非线性成分,以熵分析和幂律分析为代表的非线性动力学分析方法现已广泛用于 HRV 信号的分析中[15]。
因此,鉴于急性 CVD 事件高发于清晨醒来前后的事实、HRV 分析在自主神经功能评估中的重要性以及心脏动力系统的非线性特性,本文综合使用时域、频域、非线性动力学分析方法来挖掘夜间睡眠时的 HRV 信号的特征,并利用机器学习算法构建 CVD 自动预测模型,从而探索睡眠末期的 HRV 信号分析是否可以提供更高的 CVD 预测准确率以及哪些 HRV 指标在 CVD 预测中较为重要,这一结果或对 CVD 的早期预测带来新的突破。
1 研究方法
1.1 研究对象
本文采用的研究数据来源于美国睡眠心脏健康研究项目(sleep heart health study,SHHS)(网址:https://sleepdata.org/datasets/shhs)的开源数据库[23],共采纳 140 名受试者的相关数据,受试者的筛选过程在本文的前期研究中有详细的阐述[22],本文不再赘述。受试者在家中使用便携式多导睡眠监测仪采集多导睡眠图(polysomnography,PSG),其中记录了两个导联的心电图(electrocardiograph,ECG)数据,采样频率为 125 Hz。根据 SHHS 的随访记录,其中的 70 人在 PSG 采集后的两年随访期内发生了 CVD 事件,包括心绞痛、急性心肌梗塞等;另外 70 人在随访期内未报道任何 CVD 事件。本文将发生 CVD 事件的受试者记为 CVD 组,未发生 CVD 事件的受试者记为 non-CVD 组。两组人群在年龄、腰臀比、是否有糖尿病病史以及是否有高血压病史这四个指标上的差异具有统计学意义(P<0.05),而在性别、身高、身体质量指数(body mass index,BMI)、呼吸紊乱指数(apnea hypopnea index,AHI)、呼吸障碍指数(respiratory disturbance index,RDI)、吸烟状况和吸烟频率等指标上的差异没有统计学意义(P>0.05)。
1.2 信号预处理
对于每位受试者,其整夜 ECG 信号经巴特沃斯带通滤波器(0.5~45 Hz)预处理后,根据潘-汤普金斯(Pan-Tompkins,PT)方法检测出所有的 R 波顶点[24],然后计算两个连续 R 波峰之间的间隔,得到 R-R 间期序列。接着,去除 R-R 间期序列中的异位起搏点或伪迹并进行插值处理[25],以获得供后续分析使用的 HRV 信号。
为了对比不同睡眠段 HRV 信号对 CVD 事件的预测能力,如图 1 所示,本文将每位受试者整夜睡眠 HRV 信号划分成 5 个 1 h 的片段,分别定义为睡眠段 1~5,其中睡眠段 5 对应为清晨觉醒前 1 h。若整晚睡眠时长大于 5 h,则每个片段互不重叠,片段之间的分离间期相等;若整晚睡眠时长小于 5 h,则每个片段有重叠,片段之间的重叠间期相等。

1.3 HRV 分析
本文综合考虑了 HRV 的时域、频域和非线性分析方法。时域和非线性分析参数的计算基于每 1 h 的 HRV 数据段进行,而频域参数的计算则基于每 5 min 的 HRV 片段,然后以 1 h 内的 12 个片段的均值作为对应段的特征值。时域参数包括所有 R-R 间期(NN 间期)的标准差(standard deviation of normal to normal,SDNN)、相邻 R-R 间期差值的均方根值(root mean square of successive differences,RMSSD),频域参数包括 0.003~0.4 Hz 频率范围内的总功率(total power,TP)、0.04~0.15 Hz 频率范围内的低频功率(power in low frequency,LF)、0.15~0.4 HZ 频率范围内的高频功率(power in high frequency,HF)以及归一化的高频段能量(HF power in normalized units,HFnorm)[26]。对于非线性动力学分析,本文采用了去趋势波动分析(detrend fluctuation analysis,DFA)和多尺度熵分析。
1.3.1 DFA 分析简介
采用传统方法研究时间序列的长期幂律关系可能会因为噪声或者趋势成分的影响而使得分析结果不可靠[27]。DFA 分析则可以通过去除趋势项将信号与外部环境的噪声分离,从而适用于非平稳时间序列的处理。对于长度为 N 的序列,DFA 算法描述如下:
首先,如式(1)所示,将 N} 零均值处理后进行一阶积分,得到一个新的序列
:
![]() |
其次,给定区间长度 n,将 划分为不重叠的 m 个区间,m 为 N/n 的整数部分。对于划分的每一段序列采用最小二乘法线性拟合出局部趋势
,最后对
去除每个区间的局部趋势并计算出序列的均方根,如式(2)所示:
![]() |
随着 n 的变化呈现幂律分布,其中的斜率表示 DFA 指数,本文记为 DFA_alpha,描述了时间序列的粗糙程度,其值越大表明时间序列越光滑。DFA 是一种能够有效地量化嵌入在非平稳时间序列中的长期相关性的有效方法[28]。
1.3.2 多尺度熵分析
多尺度熵分析引进了时间尺度因子 τ 的概念,通过构造一系列不同粗粒化程度的时间序列,并计算它们的熵值,以获得原始时间序列在不同尺度下的复杂性测度。其具体计算过程包括粗粒化和熵值计算两个部分。
对于 N 点的时间序列 ,粗粒化过程如图 2 所示[29]。以尺度因子 τ 为长度,依次将原时间序列划分为 N/τ 段不重叠的子序列,这些子序列的均值就形成了
在尺度 τ 上的粗粒化序列,记为
。

对于每个粗粒化序列,本文分别采用样本熵(sample entropy,SE)和改进的排列熵(modified permutation entropy, MPE)算法来进行熵分析。SE 算法如式(3)~(6)所示[30]。
将时间序列嵌入到 m 维空间中,该空间中向量如式(3)所示:
![]() |
对于 m 维空间中的向量 Xim,计算它与该空间中其他任一向量 Xim 之间的距离,该距离表示为向量中对应维度上元素距离的最大值,如式(4)所示:
![]() |
然后,对于每一个向量,计算 m 维空间中与之距离小于给定相似容限 r 的向量占比,记为 。对于所有向量,计算
的均值,记为
,如式(5)所示:
![]() |
将嵌入空间的维数增加 1 后重复以上步骤得到 ,由此可计算出 SE 值(符号记为 SE),如式(6)所示:
![]() |
在本文中,嵌入维 m 取值设定为 2,相似容限 r 的取值设定为 0.15×std,std 为原始时间序列的标准差。
具体 PE 算法为[31]:在 m 维的嵌入空间中,对每一个 m 维的向量 Xim 中的元素数值大小按升序重新排列,结果中的每个元素在原始序列中的位序作为一个符号记录下来,则得到由 m 个符号构成的一个排列。以 Pj 来表示每种可能出现的排列的频率,并计算它们的香农熵,即可得到序列的 PE 值(符号记为 PE),如式(7)所示:
![]() |
在 PE 算法中,向量中相等的元素值会被赋予不同的符号。Bian 等[32]提出了 MPE 算法,其中使用相同的符号来表示等值情况。相比 PE 算法,MPE 算法在不同的生理和病理条件下,能够更有效地表征 HRV 信号的复杂性[32]。因此,本文使用 MPE 算法。
本文对每一段 1 h 的 HRV 信号进行了尺度因子从 1 逐一递增至 10 的粗粒化;对于每个粗粒化的序列分别计算了 SE、MPE 值,分别记为 SE1~SE10 和 MPE1~MPE10;同时,还计算了 10 个尺度下 SE、MPE 的均值,分别记为 aver_SE 和 aver_MPE。
1.4 CVD 预测模型构建
本文提出一种两层模型方法来实现 CVD 组和 non-CVD 组的二分类,如图 3 所示。第一层是基于极端梯度提升算法(extreme gradient boosting,XGBoost)的分类器,用于特征筛选。其输入为某段 HRV 信号的时域特征(SDNN、RMSSD)、频域特征(TP、LF、HF、HFnorm)、非线性动力学特征(DFA_alpha、SE1~SE10、MPE1~MPE10、aver_SE、aver_MPE)以及 11 个临床特征(年龄、BMI、身高、腰臀比、吸烟状况、吸烟频率、是否有高血压病史、是否有糖尿病病史、AHI、RDI、性别),共计 40 个特征。使用基于机器学习库的网格搜索及交叉验证方法来寻找最优分类器[33],并对其特征重要性进行排序。在第二层模型中,按照上层模型的特征重要性从高到低依次挑选 2~40 的特征,并构建机器学习模型的输入特征空间,建立 CVD 和 non-CVD 分类器。对于机器学习算法的选择,本文考察了 XGBoost、随机森林(random forest,RF)和支持向量机(support vector machine,SVM)三种算法,来寻找分类性能最优的模型。为了提高机器学习模型的泛化能力,每个模型均采用五折交叉验证来划分训练集和测试集,即将整个数据集随机划分为 5 个大小几乎相同的不相交的子集,使用其中 1 折作为测试数据,其余的 4 折作为训练数据,直到所有折都被依次用作测试数据集。最终对 5 轮预测结果求平均值,从而能更合理地评估模型的泛化能力[34]。

1.5 模型性能评价
典型的二分类问题中存在两个类,定义为“正类”和“负类”。实际为正类并且被预测为正类的样本称为“真正类”(true positive,TP)(符号记为 TP),实际为正类却被预测为负类的样本称为“假负类”(false negative,FN)(符号记为 FN);类似地,实际为负类并且被预测为负类的样本称为“真负类”(true negative,TN)(符号记为 TN),实际为负类却被预测为正类的样本称为“假正类”(false positive,FP)(符号记为 FP)[35-36]。
本文将 CVD 组定义为正类,non-CVD 组定义为负类,并使用六个常用指标评估模型预测性能,即:准确率(accuracy,ACC)(符号记为 ACC),Kappa 系数(符号记为 Kappa),真阴性率(true negative rate,TNR)(符号记为 TNR),真阳性率(true positive rate,TPR)(符号记为 TPR),阳性预测值(positive predictive value,PPV)(符号记为 PPV),F1 值(F1-Score,F1)(符号记为 F1),马修斯相关系数(matthews correlation coefficient,MCC)(符号记为 MCC)。每个指标具体计算方法如式(8)~(14)所示。ACC、Kappa、MCC 常用于综合评价正负类的预测表现,而 TPR、PPV 和 F1 偏向于衡量正类预测表现,TNR 偏向于衡量负类预测表现。
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
2 结果与分析
本文利用提出的两层模型方法探索了基于哪一个睡眠段的 HRV 信号分析能更好地预测 CVD 的问题。如表 1 所示,以五折交叉验证的 ACC 值作为比较标准,列出了利用不同睡眠段 HRV 构建的最优模型的分类性能。使用睡眠段 5,第二层模型采用 XGBoost 算法,使用 13 个输入特征,可以得到正负类综合预测性能最高的模型,其 ACC 值为 84.3%,Kappa 系数为 0.686,MCC 值为 0.689,表示模型的预测值与真实值之间有高度的一致性。表 1 的结果表明,相比其他睡眠时段,觉醒前 HRV 在预测 CVD 事件上具有更高价值。

如图 4 所示,进一步展示了利用不同睡眠段 HRV 构建的最优模型的特征重要性图。对于睡眠段 5,列出了第一层模型中所有特征的重要性排序;而对于睡眠段 1~睡眠段 4,只列出了 ACC 最大情况下筛选出的用于第二层模型的特征。结果显示,无论使用哪一段睡眠数据,年龄与腰臀比都是模型中最重要的特征,这与文献[37]的报道是吻合的;基于非线性动力学分析的 HRV 指标在各个睡眠段,都有入选。对于睡眠段 3,其 ACC 最低(如表 1 所示),入选的特征中 HRV 指标较少。而对于睡眠段 5,在筛选特征个数为 13 时 ACC 数值最大,在这一筛选的特征中 HRV 指标占比较高。这一结果可能意味着模型中 HRV 特征重要性的欠缺会导致较低的 CVD 预测正确率。

虽然对于各个睡眠段,入选第二层的特征有所不同,如图 4 所示,但可以发现,基于非线性动力学分析的 MPE1 指标在每一睡眠段中都会出现,SE3 指标也出现在大多数睡眠段中。对于睡眠段 5,如图 5 所示,组间比较(非配对 t 检验)的结果表明,相对于 non-CVD 组而言,CVD 组的 MPE1 上升,SE3 下降,差异具有统计学意义(P<0.01)。



3 结论
本文提出并采用了一种两层机器学习模型的方法,通过比较利用不同睡眠时期 HRV 特征构建 CVD 预测模型的性能,来探索是否觉醒前的 HRV 更有利于 CVD 事件预测。在正类和负类均衡的情况下,利用觉醒前 HRV 特征构建的预测模型对于 CVD 组和 non-CVD 组的分类准确率最高,达 84.3%。这一结果表明,与其他睡眠段 HRV 相比,觉醒前 HRV 能更准确地预测 CVD 事件。同时,通过对比不同睡眠段 HRV 构建的最优模型的特征重要性排序图,发现 CVD 预测正确率的高低与 HRV 特征重要性占比有关,而 HRV 分析的欠缺可能会导致较低的 CVD 预测正确率。
本研究也表明,夜间睡眠时 HRV 信号的非线性动力学分析也许能更好地捕捉心血管系统潜在的危险。相对于 non-CVD 组,CVD 组在觉醒前 HRV 上,MPE1 和 SE3 都出现了显著的变化。此外,虽然对于各个不同时期的睡眠段,各 HRV 分析指标在模型中的重要性有所不同,但非线性动力学指标 MPE1、SE3 在各睡眠段中均展现出了较高的特征重要性,意味着它们对 CVD 的预测较为重要。本文研究结果进一步证实了心脏动力系统具有非线性特性,在构建 CVD 预测模型时,应综合 HRV 的线性与非线性动力学指标,才能更好地评估潜在的生理过程。
利益冲突声明:本文全体作者均声明不存在利益冲突。
引言
心血管疾病(cardiovascular disease,CVD)是由心脏和血管异常引起的一系列疾病的总称。目前,我国 CVD 患病率处于持续上升阶段,而因 CVD 致死则占城乡居民总死亡原因的首位[1],因此 CVD 的早期预测对人类的健康和社会发展有着积极的意义[2]。
睡眠,作为一种复杂的生理过程,与机体的各基本系统密切相关,其中就包括心血管系统[3]。心血管系统的稳定性受到众多因素的控制,其中自主神经系统对心率的调节起着至关重要的作用。自主神经系统包括交感神经和副交感神经,交感神经兴奋时,心率会增加;副交感神经兴奋时,心率会降低。自主神经系统活动失衡会增加许多疾病(包括 CVD)的发病率和死亡率[4]。研究表明,在一天的不同时刻,急性 CVD 事件的发生率是不同的,清晨醒来时刻前后是 CVD 事件发生的高峰期[5-9]。虽然导致清晨 CVD 事件增加的原因并未明确,但不少研究认为这一状况可能与夜间睡眠末期交感神经活动的激增有关[10-14]。
自主神经系统使得心血管系统能够对机体内外的易变刺激做出迅速反应,导致逐次心跳的心率呈现出细小的波动,通常被称为心率变异性(heart rate variability,HRV)。HRV 分析方法主要包括时域分析法、频域分析法和非线性动力学分析法[15]。作为一种非侵入、无创的手段,HRV 分析已经广泛用于评估自主神经系统功能和心血管系统健康状况[16-18]。在文献中,也有利用 HRV 分析进行 CVD 风险评估的报道[19-21],但这些研究存在着一些不足。首先,使用的 HRV 数据并非针对睡眠状态。而在本文的前期研究中发现,睡眠 HRV 分析能够有效预测短时的 CVD 事件[22]。其次,采用的分析方法大都基于时域、频域分析,鲜有涉及非线性动力学分析。然而大量的研究表明,心脏动力系统富含非线性成分,以熵分析和幂律分析为代表的非线性动力学分析方法现已广泛用于 HRV 信号的分析中[15]。
因此,鉴于急性 CVD 事件高发于清晨醒来前后的事实、HRV 分析在自主神经功能评估中的重要性以及心脏动力系统的非线性特性,本文综合使用时域、频域、非线性动力学分析方法来挖掘夜间睡眠时的 HRV 信号的特征,并利用机器学习算法构建 CVD 自动预测模型,从而探索睡眠末期的 HRV 信号分析是否可以提供更高的 CVD 预测准确率以及哪些 HRV 指标在 CVD 预测中较为重要,这一结果或对 CVD 的早期预测带来新的突破。
1 研究方法
1.1 研究对象
本文采用的研究数据来源于美国睡眠心脏健康研究项目(sleep heart health study,SHHS)(网址:https://sleepdata.org/datasets/shhs)的开源数据库[23],共采纳 140 名受试者的相关数据,受试者的筛选过程在本文的前期研究中有详细的阐述[22],本文不再赘述。受试者在家中使用便携式多导睡眠监测仪采集多导睡眠图(polysomnography,PSG),其中记录了两个导联的心电图(electrocardiograph,ECG)数据,采样频率为 125 Hz。根据 SHHS 的随访记录,其中的 70 人在 PSG 采集后的两年随访期内发生了 CVD 事件,包括心绞痛、急性心肌梗塞等;另外 70 人在随访期内未报道任何 CVD 事件。本文将发生 CVD 事件的受试者记为 CVD 组,未发生 CVD 事件的受试者记为 non-CVD 组。两组人群在年龄、腰臀比、是否有糖尿病病史以及是否有高血压病史这四个指标上的差异具有统计学意义(P<0.05),而在性别、身高、身体质量指数(body mass index,BMI)、呼吸紊乱指数(apnea hypopnea index,AHI)、呼吸障碍指数(respiratory disturbance index,RDI)、吸烟状况和吸烟频率等指标上的差异没有统计学意义(P>0.05)。
1.2 信号预处理
对于每位受试者,其整夜 ECG 信号经巴特沃斯带通滤波器(0.5~45 Hz)预处理后,根据潘-汤普金斯(Pan-Tompkins,PT)方法检测出所有的 R 波顶点[24],然后计算两个连续 R 波峰之间的间隔,得到 R-R 间期序列。接着,去除 R-R 间期序列中的异位起搏点或伪迹并进行插值处理[25],以获得供后续分析使用的 HRV 信号。
为了对比不同睡眠段 HRV 信号对 CVD 事件的预测能力,如图 1 所示,本文将每位受试者整夜睡眠 HRV 信号划分成 5 个 1 h 的片段,分别定义为睡眠段 1~5,其中睡眠段 5 对应为清晨觉醒前 1 h。若整晚睡眠时长大于 5 h,则每个片段互不重叠,片段之间的分离间期相等;若整晚睡眠时长小于 5 h,则每个片段有重叠,片段之间的重叠间期相等。

1.3 HRV 分析
本文综合考虑了 HRV 的时域、频域和非线性分析方法。时域和非线性分析参数的计算基于每 1 h 的 HRV 数据段进行,而频域参数的计算则基于每 5 min 的 HRV 片段,然后以 1 h 内的 12 个片段的均值作为对应段的特征值。时域参数包括所有 R-R 间期(NN 间期)的标准差(standard deviation of normal to normal,SDNN)、相邻 R-R 间期差值的均方根值(root mean square of successive differences,RMSSD),频域参数包括 0.003~0.4 Hz 频率范围内的总功率(total power,TP)、0.04~0.15 Hz 频率范围内的低频功率(power in low frequency,LF)、0.15~0.4 HZ 频率范围内的高频功率(power in high frequency,HF)以及归一化的高频段能量(HF power in normalized units,HFnorm)[26]。对于非线性动力学分析,本文采用了去趋势波动分析(detrend fluctuation analysis,DFA)和多尺度熵分析。
1.3.1 DFA 分析简介
采用传统方法研究时间序列的长期幂律关系可能会因为噪声或者趋势成分的影响而使得分析结果不可靠[27]。DFA 分析则可以通过去除趋势项将信号与外部环境的噪声分离,从而适用于非平稳时间序列的处理。对于长度为 N 的序列,DFA 算法描述如下:
首先,如式(1)所示,将 N} 零均值处理后进行一阶积分,得到一个新的序列
:
![]() |
其次,给定区间长度 n,将 划分为不重叠的 m 个区间,m 为 N/n 的整数部分。对于划分的每一段序列采用最小二乘法线性拟合出局部趋势
,最后对
去除每个区间的局部趋势并计算出序列的均方根,如式(2)所示:
![]() |
随着 n 的变化呈现幂律分布,其中的斜率表示 DFA 指数,本文记为 DFA_alpha,描述了时间序列的粗糙程度,其值越大表明时间序列越光滑。DFA 是一种能够有效地量化嵌入在非平稳时间序列中的长期相关性的有效方法[28]。
1.3.2 多尺度熵分析
多尺度熵分析引进了时间尺度因子 τ 的概念,通过构造一系列不同粗粒化程度的时间序列,并计算它们的熵值,以获得原始时间序列在不同尺度下的复杂性测度。其具体计算过程包括粗粒化和熵值计算两个部分。
对于 N 点的时间序列 ,粗粒化过程如图 2 所示[29]。以尺度因子 τ 为长度,依次将原时间序列划分为 N/τ 段不重叠的子序列,这些子序列的均值就形成了
在尺度 τ 上的粗粒化序列,记为
。

对于每个粗粒化序列,本文分别采用样本熵(sample entropy,SE)和改进的排列熵(modified permutation entropy, MPE)算法来进行熵分析。SE 算法如式(3)~(6)所示[30]。
将时间序列嵌入到 m 维空间中,该空间中向量如式(3)所示:
![]() |
对于 m 维空间中的向量 Xim,计算它与该空间中其他任一向量 Xim 之间的距离,该距离表示为向量中对应维度上元素距离的最大值,如式(4)所示:
![]() |
然后,对于每一个向量,计算 m 维空间中与之距离小于给定相似容限 r 的向量占比,记为 。对于所有向量,计算
的均值,记为
,如式(5)所示:
![]() |
将嵌入空间的维数增加 1 后重复以上步骤得到 ,由此可计算出 SE 值(符号记为 SE),如式(6)所示:
![]() |
在本文中,嵌入维 m 取值设定为 2,相似容限 r 的取值设定为 0.15×std,std 为原始时间序列的标准差。
具体 PE 算法为[31]:在 m 维的嵌入空间中,对每一个 m 维的向量 Xim 中的元素数值大小按升序重新排列,结果中的每个元素在原始序列中的位序作为一个符号记录下来,则得到由 m 个符号构成的一个排列。以 Pj 来表示每种可能出现的排列的频率,并计算它们的香农熵,即可得到序列的 PE 值(符号记为 PE),如式(7)所示:
![]() |
在 PE 算法中,向量中相等的元素值会被赋予不同的符号。Bian 等[32]提出了 MPE 算法,其中使用相同的符号来表示等值情况。相比 PE 算法,MPE 算法在不同的生理和病理条件下,能够更有效地表征 HRV 信号的复杂性[32]。因此,本文使用 MPE 算法。
本文对每一段 1 h 的 HRV 信号进行了尺度因子从 1 逐一递增至 10 的粗粒化;对于每个粗粒化的序列分别计算了 SE、MPE 值,分别记为 SE1~SE10 和 MPE1~MPE10;同时,还计算了 10 个尺度下 SE、MPE 的均值,分别记为 aver_SE 和 aver_MPE。
1.4 CVD 预测模型构建
本文提出一种两层模型方法来实现 CVD 组和 non-CVD 组的二分类,如图 3 所示。第一层是基于极端梯度提升算法(extreme gradient boosting,XGBoost)的分类器,用于特征筛选。其输入为某段 HRV 信号的时域特征(SDNN、RMSSD)、频域特征(TP、LF、HF、HFnorm)、非线性动力学特征(DFA_alpha、SE1~SE10、MPE1~MPE10、aver_SE、aver_MPE)以及 11 个临床特征(年龄、BMI、身高、腰臀比、吸烟状况、吸烟频率、是否有高血压病史、是否有糖尿病病史、AHI、RDI、性别),共计 40 个特征。使用基于机器学习库的网格搜索及交叉验证方法来寻找最优分类器[33],并对其特征重要性进行排序。在第二层模型中,按照上层模型的特征重要性从高到低依次挑选 2~40 的特征,并构建机器学习模型的输入特征空间,建立 CVD 和 non-CVD 分类器。对于机器学习算法的选择,本文考察了 XGBoost、随机森林(random forest,RF)和支持向量机(support vector machine,SVM)三种算法,来寻找分类性能最优的模型。为了提高机器学习模型的泛化能力,每个模型均采用五折交叉验证来划分训练集和测试集,即将整个数据集随机划分为 5 个大小几乎相同的不相交的子集,使用其中 1 折作为测试数据,其余的 4 折作为训练数据,直到所有折都被依次用作测试数据集。最终对 5 轮预测结果求平均值,从而能更合理地评估模型的泛化能力[34]。

1.5 模型性能评价
典型的二分类问题中存在两个类,定义为“正类”和“负类”。实际为正类并且被预测为正类的样本称为“真正类”(true positive,TP)(符号记为 TP),实际为正类却被预测为负类的样本称为“假负类”(false negative,FN)(符号记为 FN);类似地,实际为负类并且被预测为负类的样本称为“真负类”(true negative,TN)(符号记为 TN),实际为负类却被预测为正类的样本称为“假正类”(false positive,FP)(符号记为 FP)[35-36]。
本文将 CVD 组定义为正类,non-CVD 组定义为负类,并使用六个常用指标评估模型预测性能,即:准确率(accuracy,ACC)(符号记为 ACC),Kappa 系数(符号记为 Kappa),真阴性率(true negative rate,TNR)(符号记为 TNR),真阳性率(true positive rate,TPR)(符号记为 TPR),阳性预测值(positive predictive value,PPV)(符号记为 PPV),F1 值(F1-Score,F1)(符号记为 F1),马修斯相关系数(matthews correlation coefficient,MCC)(符号记为 MCC)。每个指标具体计算方法如式(8)~(14)所示。ACC、Kappa、MCC 常用于综合评价正负类的预测表现,而 TPR、PPV 和 F1 偏向于衡量正类预测表现,TNR 偏向于衡量负类预测表现。
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
2 结果与分析
本文利用提出的两层模型方法探索了基于哪一个睡眠段的 HRV 信号分析能更好地预测 CVD 的问题。如表 1 所示,以五折交叉验证的 ACC 值作为比较标准,列出了利用不同睡眠段 HRV 构建的最优模型的分类性能。使用睡眠段 5,第二层模型采用 XGBoost 算法,使用 13 个输入特征,可以得到正负类综合预测性能最高的模型,其 ACC 值为 84.3%,Kappa 系数为 0.686,MCC 值为 0.689,表示模型的预测值与真实值之间有高度的一致性。表 1 的结果表明,相比其他睡眠时段,觉醒前 HRV 在预测 CVD 事件上具有更高价值。

如图 4 所示,进一步展示了利用不同睡眠段 HRV 构建的最优模型的特征重要性图。对于睡眠段 5,列出了第一层模型中所有特征的重要性排序;而对于睡眠段 1~睡眠段 4,只列出了 ACC 最大情况下筛选出的用于第二层模型的特征。结果显示,无论使用哪一段睡眠数据,年龄与腰臀比都是模型中最重要的特征,这与文献[37]的报道是吻合的;基于非线性动力学分析的 HRV 指标在各个睡眠段,都有入选。对于睡眠段 3,其 ACC 最低(如表 1 所示),入选的特征中 HRV 指标较少。而对于睡眠段 5,在筛选特征个数为 13 时 ACC 数值最大,在这一筛选的特征中 HRV 指标占比较高。这一结果可能意味着模型中 HRV 特征重要性的欠缺会导致较低的 CVD 预测正确率。

虽然对于各个睡眠段,入选第二层的特征有所不同,如图 4 所示,但可以发现,基于非线性动力学分析的 MPE1 指标在每一睡眠段中都会出现,SE3 指标也出现在大多数睡眠段中。对于睡眠段 5,如图 5 所示,组间比较(非配对 t 检验)的结果表明,相对于 non-CVD 组而言,CVD 组的 MPE1 上升,SE3 下降,差异具有统计学意义(P<0.01)。



3 结论
本文提出并采用了一种两层机器学习模型的方法,通过比较利用不同睡眠时期 HRV 特征构建 CVD 预测模型的性能,来探索是否觉醒前的 HRV 更有利于 CVD 事件预测。在正类和负类均衡的情况下,利用觉醒前 HRV 特征构建的预测模型对于 CVD 组和 non-CVD 组的分类准确率最高,达 84.3%。这一结果表明,与其他睡眠段 HRV 相比,觉醒前 HRV 能更准确地预测 CVD 事件。同时,通过对比不同睡眠段 HRV 构建的最优模型的特征重要性排序图,发现 CVD 预测正确率的高低与 HRV 特征重要性占比有关,而 HRV 分析的欠缺可能会导致较低的 CVD 预测正确率。
本研究也表明,夜间睡眠时 HRV 信号的非线性动力学分析也许能更好地捕捉心血管系统潜在的危险。相对于 non-CVD 组,CVD 组在觉醒前 HRV 上,MPE1 和 SE3 都出现了显著的变化。此外,虽然对于各个不同时期的睡眠段,各 HRV 分析指标在模型中的重要性有所不同,但非线性动力学指标 MPE1、SE3 在各睡眠段中均展现出了较高的特征重要性,意味着它们对 CVD 的预测较为重要。本文研究结果进一步证实了心脏动力系统具有非线性特性,在构建 CVD 预测模型时,应综合 HRV 的线性与非线性动力学指标,才能更好地评估潜在的生理过程。
利益冲突声明:本文全体作者均声明不存在利益冲突。