基于Lasso-Logistic回归的帕金森疾病声学特征诊断研究

Diagnostic Study of Acoustic Features in Parkinson’s Disease Based on Lasso-Logistic Regression

摘要: 本文基于埃斯特雷马杜拉(西班牙)帕金森病区域协会的80名欧洲受试者声学特征数据，结合lasso回归提出两阶段变量选择法对44个声学特征因子筛选，最后得到6个显著的声学特征因子：Gender、Shim_loc、MFCC3、HNR35、PPE、GNE。将上述因素通过多因素logistic回归构建患PD疾病风险的列线图模型，并从多个角度验证该模型的有效性和校准性。结果表明，早期PD患者基底神经节运动调节功能异常，声学数据中MFCC3和HNR35数值偏低，PPE、GNE和Shim_loc数值偏高，临床表现为发音时声带振动的最大频率降低，声音低沉，进一步说明所构建的列线图模型可以根据不同的声学特征较好地诊断研究对象患有PD疾病的风险高低。今后声学特征有望成为早期PD诊断的重要生物标记物，为疾病远程筛查提供辅助手段。

Abstract: In this paper, based on the acoustic profile data of 80 European subjects from the Regional Association of Parkinson’s Disease in Extremadura (Spain), a two-stage variable selection method was proposed to screen 44 acoustic profile factors in conjunction with lasso regression, and finally six significant acoustic profile factors were obtained: gender, Shim_loc, MFCC3, HNR35, PPE, and GNE. The above factors were combined to construct a column plot model of the risk of developing PD disease by multifactorial logistic regression to construct a column-line graph model of the risk of developing PD disease, and validate the validity and calibration of the model from multiple perspectives. The results showed that early PD patients had abnormal motor regulation of the basal ganglia, low values of MFCC3 and HNR35, high values of PPE, GNE and Shim_loc in the acoustic data, and reduced maximal frequency of vocal fold vibration during articulation and muffled voices, which indicated that the constructed columnar plot model could diagnose the risk of PD in the subjects according to the different acoustic features. In the future, acoustic features are expected to become important biomarkers for early PD diagnosis and provide an aid for remote screening of the disease.

1. 引言

帕金森病(PD)是一种常见的神经退行性疾病，临床表现的特征是静止性震颤，肌强直，运动迟缓，姿势步态障碍等运动症状。PD的运动控制机理十分复杂，由于没有明确的诊断测试，诊断任务十分困难，特别是在运动症状出现的早期阶段因不严重而难以诊断。国外有学者研究发现，在PD的早期阶段，90%患者可能会出现不同程度的语言障碍。PD患者的语音障碍主要是构音异常和发声减弱，表现为语速减慢、语言停顿异常、持续发声障碍、重音异常、单一高音、强度见多、发声与呼吸同步性降低等等。合并语音障碍的PD患者更易出现不良的生理和心理状况，进而严重影响生活质量。

目前的临床工作中，神经科医师接诊PD患者时，重心多集中在肢体运动症状，例如：震颤、僵硬、动作迟缓等，往往对语音障碍的关注度较低。此外，由于人类听觉能力的局限性，在PD早期仅仅凭借听觉无法灵敏捕捉到患者语音特征的异常变化，使得语音障碍很难被感知。因此，根据PD患者早期阶段的声学特征进行分析，对疾病的诊断非常必要，对病情评估及康复治疗也有重要意义。本文基于Naranjo等人研究的80名欧洲受试者声学特征数据 [1] ，结合lasso回归提出两阶段变量选择法对显著声学特征因子进行筛选，将筛选后的因子使用多因素logistic回归建立的预测模型构建列线图，同时验证其有效性，利用有效的声学特征为早期帕金森疾病诊断提供可能的理论依据和科学手段。

2. 方法

2.1. 一般资料

本文引用Naranjo等人研究的80名欧洲受试者声学特征数据进行研究 [1] 。Naranjo等通过考虑重复录音，进行了PD受试者与健康个体的区分实验 [1] 。共有80名50岁以上的受试者参与了这项研究。其中40例健康：男性22例(55%)和女性18例(45%)，其中40例受帕金森病影响：男性27例(67.5%)，女性13例(32.5%)。对照组的平均年龄(±标准差)为66.38 ± 8.38岁，帕金森病患者的平均年龄为69.58 ± 7.82岁。参与这项研究的帕金森病患者是埃斯特雷马杜拉(西班牙)帕金森病区域协会的成员。该研究方案得到了埃斯特雷马杜拉大学生物伦理委员会的批准。所有受试者都签署了知情同意书。每个语音记录被处理以提供44个声学特征，即每个语音记录一个44维向量。根据提取的特征是否有相关的表述，将其分成6组：性别、音高局部扰动测量、振幅局部扰动测量、噪声特征、频谱包络测量(基于梅尔频率外频系数的0至12阶频谱测量及其导数)和非线性特征，详情见表1。

Table 1. Key indicators of the acoustic characterisation dataset

表1. 声学特征数据集主要指标

2.2. 两阶段变量选择法

通过观察声学特征可知，6组特征集的其中4组：音高局部扰动测量c₂、振幅局部扰动测量c₃、噪声特征c₄、频谱包络测量c₅之间的变量存在高度相关。以相对抖动和相对平均扰动为例，二者的皮尔逊相关系数高达0.99，不同组别的各变量之间也存在高度相关的情况(详见附件)。为了避免各变量间存在的多重共线性问题，在变量筛选的第一阶段从4组中选择一个代表性特征是合理的。6组特征集的44个特征变量，经过第一阶段筛选将变成10个特征量。在第二阶段，将具有统计学意义的变量通过LASSO方法进行筛选，选择系数不为0的特征变量。

在第一阶段中，本文引入Naranjo等人的筛选方法 [2] ，利用两变量间简单相关系数 $\rho \left(x_k,x_l\right)$ 定义二者的不相似度，其表达式见公式(1)，并定义每个变量相对于其组中其他变量的差异度量为 $\delta \left(x_k\right)$ ，见公式(2)。然后，变量选择标准在于找到每组中差异最小的变量 $x_{kg}=\arg {\min }_{l\in c_g}\delta \left(x_k\right)$ 。这为每个组提供了一个表示变量。显然，当每组只有一个变量时，这个准则的特殊应用提供了唯一的变量作为代表变量。

$d\left(x_k,x_l\right)=1-\left|\rho \left(x_k,x_l\right)\right|$ ，其中 $\rho \left(x_k,x_l\right)=\frac{{\displaystyle {\sum }_{i=1}^n\left(x_{ik}-{\bar{x}}_k\right)\left(x_{il}-{\bar{x}}_l\right)}}{\sqrt{{\displaystyle {\sum }_{i=1}^n{\left(x_{ik}-{\bar{x}}_k\right)}^2}}\sqrt{{\displaystyle {\sum }_{i=1}^n{\left(x_{il}-{\bar{x}}_l\right)}^2}}}$ (1)

$\delta \left(x_k\right)={\displaystyle {\sum }_{l\in c_g}d\left(x_k,x_l\right)}$ (2)

在第二阶段中，考虑到数据集中存在“性别”计数资料和其他计量资料，利用R软件对计量资料进行正态性检验，符合正态分布的计量资料则以x ± s表示，采用独立样本t检验比较各指标在健康组与患病组之间的分布差异，非正态分布的计量值则用中位数(四分位间距)表示，采用Mann-Whitney U检验进行组间比较；计数资料以例数或百分比表示，组间比较采用χ2检验。再根据单因素分析结果，将有统计学意义的变量通过LASSO方法进行筛选，选择系数不为0的特征变量 [3] 。Lasso回归模型 [4] 相较于多元线性回归模型，在计算损失函数时采用了改良后的普通最小二乘估计，通过最小二乘估计改进特征间的共线性影响，其关键点在于在线性回归的损失函数后加一个L1正则化项，如式(3)所示：

${\min }_{\omega }\left({\Vert y-X\omega \Vert }_2^2+\lambda {\Vert \omega \Vert }_1\right),\lambda >0$ (3)

式中，X为输入的特征矩阵；y为输出矩阵； $\omega$ 为模型的参数向量； $\lambda$ 为惩罚系数。Lasso回归方程的解见式(4)，其中I为单位矩阵。

$\omega =\arg {\min }_{\omega }\left({\Vert y-X\omega \Vert }_2^2+\lambda {\Vert \omega \Vert }_1\right)={\left(X^{\text{T}}X\right)}^{-1}\left(X^{\text{T}}y-0.5\lambda I\right)$ (4)

采用CV准则选择合适的惩罚系数，其中训练集可表示为 $\left\{\left(x_1,y_1\right),\cdots ,\left(x_{t-1},y_{t-1}\right),\left(x_{t+1},y_{t+1}\right),\cdots ,\left(x_n,y_n\right)\right\}$ ，测试集 $\left(x_t,y_t\right)$ 在训练集上取拟合lasso回归模型，在测试集中求得y的拟合值为 $\stackrel{⌣}{y}$ ，计算均方误差 $MSE={\left(y_t-{\stackrel{⌣}{y}}_t\right)}^2$ ，一直计算n次，得到n个均方误差的均值 $CV\left(\lambda \right)$ ，以此来选择惩罚系数 $\lambda$ 。

$CV(\lambda )=\frac{1}{n}{\displaystyle {\sum }_{j=1}^{n}MSE,j}$ (5)

2.3. 建立多因素Logistic回归列线图模型

通过前文声学特征筛选后，通过多因素logistic回归分析将筛选出的显著因子纳入模型，利用R软件绘制列线图来预测患PD疾病的风险。最后，通过评估受试者声学特征(ROC)曲线下的面积，Calibration曲线的校准度来评价列线图预测模型的预测能力和校准性。以P < 0.05为差异有统计学意义。

Logistic回归分析属于非线性回归，它是研究因变量为二项分类或多项分类结果与某些影响因素之间关系的一种多重回归分析方法 [5] 。在疾病的病因学研究中，经常需要分析疾病的发生与各危险因素之间的定量关系。Logistic回归模型较好地解决了因变量为二分类变量无法满足线性回归的基本假设条件问题，是医学研究，特别是流行病学病因研究中最常用的分析方法之一。在本题中，设因变量y是二分类变量(取值为1代表PD患者，0代表健康者)，不同的声学特征 $x_1,x_2,\cdots$ 可以作为生物标志物反映受试者患病与否从而影响y取值。假设某受试者具有m个自变量的声学特征，即在m个自变量作用下为PD患者发生的条件概率 $p=p\left(y=1|x_1,x_2,\cdots ,x_m\right)$ ，则logistic回归模型可表示为公式(6)，其中 ${\beta }_0$ 为常数项， ${\beta }_1,\cdots ,{\beta }_m$ 为偏回归系数。

$p=\frac{\exp \left({\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_m{x}_m\right)}{1+\exp \left({\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_m{x}_m\right)}$ (6)

对(6)式进行logit变换 $\text{logit}\left(p\right)=\ln p/\left(1-p\right)$ ，即logistic回归模型可以表示成如下的线性形式。很好地将问题转化为回归的问题。对于所谓的logit变换，就是通过对p进行变换之后再次纳入回归模型，得到的模型记为“logistic回归模型”。

$\ln \frac{p}{1-p}={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_m{x}_m$ (7)

在临床医学的研究者，往往将logistic回归模型得到的结果以列线图的形式呈现出来。列线图(Alignment Diagram)，又称诺莫图(Nomogram图)，用来把多因素回归分析结果(logistic回归和cox回归)用图形方式表现出来，将多个预测指标进行整合，然后采用带有刻度的线段，按照一定的比例绘制在同一平面上，从而用以表达预测模型中各个变量之间的相互关系。根据模型中各个影响因素对结局变量的贡献程度(回归系数的大小)，给每个影响因素的每个取值水平进行赋分，然后再将各个评分相加得到总评分，最后通过总评分与结局事件发生概率之间的函数转换关系，从而计算出该个体结局事件的预测值。

3. 实证分析

3.1. 声学特征筛选结果

利用R软件分别计算4组特征集：音高局部扰动测量c₂、振幅局部扰动测量c₃、噪声特征c₄、频谱包络测量c₅之间的各变量的皮尔逊相关系数，并绘制相关热力图，发现不同组别的各变量之间存在高度相关的情况(详见附录)。利用前文定义的差异度量选取各组别中差异度最小的变量分别为：Jitter_rel、Shim_loc、HNR35、MFCC3和Delta8。再结合其余2组所含有的声学特征指标，第一阶段共筛选出10个具有代表性的声学特征：Gender、Jitter_rel、Shim_loc、HNR35、MFCC3、Delta8、RPDE、DFA、PPE和GNE。

利用R软件对患病组和健康组不同指标计量资料进行正态性检验，发现均符合正态分布。为此，利用两独立样本t检验比较组间差异；对于“Gender”的计数资料，采用χ2检验比较组间差异，得到结果见表2。由表中数据可知，第一阶段选取的10个代表性特征在两组之间均存在显著性差异，说明所选的10个声学特征在统计意义上能反映出PD患者与健康患者语音功能上的变化，进而可作为生物标志物去诊断帕金森疾病。其中RPDE和Shim_loc等指标均值在患病组和健康组之间虽没有明显差别，P值却小于0.05，在统计学中具有较为明显的差异。

Table 2. Comparison of clinical data between the diseased and healthy groups

表2. 患病组与健康组的临床资料比较

将患病组与健康组之间差异具有统计学意义的10项声学特征进行lassou回归，利用R软件绘制各指标回归系数路径图和交叉验证图(见图1)，选取最优 $\lambda =0.0733473$ ，此时保留6个回归系数非0的特征指标：Gender、Shim_loc、MFCC3、HNR35、PPE和GNE。

Figure 1. Path diagram of Lasso regression coefficients (left) with cross-validation (right)

图1. Lasso回归系数路径图(左)与交叉验证图(右)

3.2. Logistic回归列线图进行PD诊断

利用R软件将选取的6个声学特征变量(Gender、Shim_loc、MFCC3、HNR35、PPE、GNE)纳入多因素logistic回归分析，并构建制研究对象是否患有帕金森疾病的预测列线图，见图2。每个声学特征变量可在评分轴分别得到一个具体分值，再将6个评分值相加得到总评分，在总评分轴上找到相应得分位置，即可得到患者是否患有PD疾病的概率。如研究对象的性别为女，Shim_loc为0.08、MFCC3为1.1、HNR35 < 75、PPE ≥ 0.6、GNE < 0.93，则总评分大致为138分，对应的总评分风险轴上的预测概率为0.96，因此该研究对象被判定为PD疾病的高危人群 [6] 。通过建立logistic回归列线图预测模型，即可根据研究对象的发声强弱、语速、语言停顿等声学特征对患病概率进行预测，从而诊断研究对象是否患PD疾病。

Figure 2. Diagnostic line diagram of PD diseases

图2. PD疾病诊断列线图

对该模型采用Bootstrap法进行内部验证，重复抽样次数为1000次，同时绘制出Calibration曲线，见图3。校正曲线和理想曲线贴合良好，说明该模型具备较好的预测能力。利用ROC曲线分析该模型预测患PD风险的效率，AUC为0.923125 (95% CI: 0.8912~0.9612)，灵敏度为0.8900，特异度为0.6808，见图4。依据Decision曲线(见图4)分析，当阈概率在20%以上时，使用该模型诊断是否患PD疾病的净收益率明显大于使用全部声学特征指标进行诊断的净收益率，且该模型净收益率是在大于0，说明该模型效果良好，对患病风险评估具有临床意义。

Figure 3. Calibration curve

图3. Calibration曲线

Figure 4. ROC curve (left) and Decision curve (right)

图4. ROC曲线(左)与Decision曲线(右)

4. 讨论

在PD患者中，声带僵硬和弯曲引起声带质量和张力的变化。因此，主体不能适当地保持持续发声，表现出不稳定的基频，即高抖动 [7] 。这与得到的结果相匹配，表明Shim_loc是PD的一个重要特征。此外，声带弯曲还会导致声门关闭不全，使得未发声的空气泄漏，产生多余的噪声，因此预期的HNR值会降低 [7] 。这与我们的实验结果一致，表明GNE、MFCC3和HNR35也是PD的特征指标。最后，在Jan Rusz等人的研究中 [8] ，患有睡眠障碍的PD患者发声响度降低，声音嘶哑，进一步证实了PPE在PD诊断中的价值。综上所述，我们通过多因素logistic回归构建了一个预测PD风险的列线图模型，该模型基于6大特征：Gender、Shim_loc、MFCC3、HNR35、PPE和GNE。这些特征可以有效地用于诊断PD的风险，有助于早期诊断和患者的个体化治疗。除此之外，对于特征变量的筛选，本文引用了一种两阶段变量选择和分类的方法。对于特征量较大的实验数据，第一阶段减少变量数量，第二阶段采用基于lasso的变量选择和分类方法，从而避免多重共线性，为今后变量筛选提供参考。

然而，目前国内针对帕金森病(PD)患者的语音声学分析研究仍然不够充分，缺乏标准的语音测试模式。此外，语言任务的多样性以及性别等因素都可能对声学特征的分析结果产生影响。因此，本研究未能充分考虑受试者内部的声学差异。另外，PD患者在进行持续发音时存在一定的困难，这可能导致一些重复的发音无法用于后续的特征提取过程。总体而言，该研究具有一定的局限性。如何在允许受试者进行重复实验的前提下，准确有效地提取和处理患病者的声学数据特征或许会成为声学特征诊断疾病的一个热门方向。随着对PD语言研究的不断深入，声学特征有望成为早期PD诊断的重要生物标记物，为疾病的远程筛查提供有力支持。

致谢

感谢老师和同学一路上的指导与帮助。

附录

1) 各组特征集不相关性 $d\left(x_k,x_l\right)$

2) 各组特征集热力图