1. 引言
帕金森病(PD)作为一种慢性神经退行性疾病,主要影响中老年人群的中枢神经系统,其流行程度仅次于阿尔茨海默病(AD) [1]。近年来,全球范围内帕金森病的患病率显著上升,其导致的残疾人数和死亡人数亦呈现增长趋势。在我国,随着人口老龄化加剧,帕金森病的发病率也呈上升趋势,预计到2030年,我国帕金森病患者数量将占全球总数的一半以上。
帕金森病的确切病因尚未明确,但普遍认为与黑质多巴胺能神经元变性死亡有关[2]。该病起病隐匿、进展缓慢,但致残性强,严重影响患者的生活质量。其运动症状主要包括静止性震颤、肌肉僵硬、运动迟缓等,同时也常伴有非运动症状如认知功能障碍、焦虑抑郁等[3]。因此,早期发现和有效管理PD对于延缓疾病进程、减轻患者痛苦至关重要。
帕金森病综合评分量表(UPDRS)作为评估PD病情的重要工具,能够全面衡量患者的精神、行为、日常活动、运动功能等多个方面,其中运动功能这一项被记作motor-UPDRS,常被单独用来评估患者疾病的严重程度。然而,传统的UPDRS评分方法依赖于医生的经验判断和患者的自我描述,具有一定的主观性。研究表明,PD患者的语音特征在疾病早期即会发生改变,这些细微的语音差异在一定程度上反映了患者的运动功能和神经系统状况。并且随着人工智能和机器学习技术的发展,基于语音数据的PD评估方法受到越来越多学者的关注,成为PD辅助治疗的工具。
因此,本研究旨在利用机器学习算法,结合UCI数据库中远程帕金森语音数据,分别使用高斯过程回归、支持向量回归、随机森林和XGBoost模型对帕金森病患者的motor-UPDRS评分进行预测,用来监测患者帕金森病发生进展情况。通过比较不同预测模型的MAE,RMSE和R2值等指标,选择最优模型为医生提供更为准确、客观的病情评估结果,有助于为PD的早期诊断、病情监测及疗效评估提供新的解决方案。
2. 国内外研究现状
Tsanas [4]等人最早利用Goetz等人的研究记录,通过Intel公司开发的AHTD系统远程收集了42名PD患者的语音数据,提取了16个对PD临床进展有用的语音特征,分别采用线性和非线性回归技术将这些特征映射到UPDRS评分,最终发现CART决策树方法预测误差最小。Asgari [5]等人则采用了不同的语音收集任务,包括持续发声任务、DDK任务和阅读任务,以获取更多的语音信息。随后通过支持向量机的不同核函数回归模型对UPDRS评分量表中的运动评分量表(motor-UPDRS)进行预测。Ömer Eskidere [6]等人首先根据变量之间的关系手动选择了不同组合的特征子集,然后比较了SVM、最小二乘支持向量机(LS-SVM)等机器学习算法对UPDRS数据的预测性能,最终发现,保留全部特征并使用LS-SVM的模型对UPDRS的预测效果最佳。Nilashi [7]等人提出了一种新的混合智能系统,该系统结合了EM算法的高斯混合模型、主成分分析的降维方法和基于自适应神经网络的模糊推理系统(ANFIS)以及支持向量机。通过对实验数据集聚类、降维和比较不同监督预测技术的性能,他们取得了很好的预测效果。一年后,为了解决帕金森病患者增量数据更新问题[8],他们进一步提出了一种结合SOM-NIPALS-ISVR的方法,该方法在减小预测误差的同时大大节省了计算时间。Shahid [9]等人则将临床信息纳入考虑,采用PCA降维方法选择了8个主成分,并结合深度学习技术建立了预测模型,这种结合临床信息和深度学习的方法在预测UPDRS评分时取得了不错的效果,特别是在决定系数上的表现得到了显著提高。
在国内研究方面,Liu [10]等人考虑到数据库中预测motor-UPDRS和total-UPDRS这两个回归任务之间的依赖关系,采用了基于模型过滤的多任务回归方法进行预测。相较于单任务预测方法,他们的方法在预测效果上得到了大幅提高。Ahmed和Zhang等人[11]提出了一种鲁棒智能回归模型,该模型结合了蚁狮优化器(BALO)的二进制版本和基于差异评估的极限学习机。他们首先使用蚁狮优化器选择语音特征,然后利用基于差异评估的极限学习机对UPDRS评分进行连续预测。这种方法在PD评估中展示了良好的鲁棒性和预测性能。
3. 数据来源与描述
3.1. 数据来源
本研究采用了UCI机器学习库中的帕金森病远程监测数据集,这一数据由牛津大学Athanasios Tsanas和Max Little教授团队联合美国医疗中心在英特尔技术支持下收集。数据集通过六个月居家远程监控获得,包含42位早期帕金森病患者(28名男性14名女性)的生物医学语音测量和临床信息,每位患者至少完成20次有效测试。录音数据经特征提取后,结合医生临床评估的线性插值UPDRS评分,形成原始数据集,其中帕金森病患者的临床信息有3个,分别是年龄,性别和测试相对时间,语音信息共有16个,本文中这里我们选择原数据集中的motor-UPDRS即w患者的运动功能评分作为预测变量。表1是数据集具体的变量描述。
Table 1. Parkinson’s disease remote monitoring dataset variables
表1. 帕金森病远程监测数据集变量
作用 |
变量名称 |
变量描述 |
ID |
Subject# |
受试者编号,用于唯一标识个体的整数 |
特征 |
age |
受试者的年龄;整数 |
sex |
受试者的性别;二分类变量:0、男性,1、女性 |
test-time |
受试者自招募入组到进行试验的时间间隔。
整数部分表示自招募以来的天数 |
Jitter (%) |
Jitter类:用于量化基频(Fundamental Frequency)
变化的一系列指标;连续型变量 |
Jitter (Abs) |
Jitter: RAP |
Jitter: PPQ5 |
Jitter: DDP |
续表
|
Shimmer |
Shimmer类:用于量化振幅(Amplitude)变化的
一系列指标;连续型变量 |
Shimmer (dB) |
Shimmer: APQ3 |
Shimmer: APQ5 |
Shimmer: APQ11 |
Shimmer: DDA |
NHR |
用于衡量语音信号中噪声与音调成分比率的
两个重要指标;连续型变量 |
HNR |
RPDE |
一种非线性动态复杂度度量;连续型 |
DFA |
信号分形缩放指数;连续型 |
PPE |
基频变化的非线性测量;连续型 |
目标 |
motor-UPDRS |
由临床医生给出的基于UPDRS运动部分的
线性插值 |
3.2. 数据描述
为了能够更好地把握帕金森病患者数据的整体趋势、分布情况以及发展进程,我们首先对帕金森远程语音数据集进行简单的可视化分析。
Figure 1. Receiver motor-UPDRS distribution map (left) and motor-UPDRS progress map of the top six subjects (right)
图1. 受试者motor-UPDRS分布图(左)与前6名受试者motor-UPDRS进展图(右)
从图1左图中可以看出受试者的motor-UPDRS得分最小值为5.0377,最大值为39.511。motor-UPDRS得分大多落在12~18和25~28这两个区间。为了观察不同受试者之间UPDRS评分的差异和变化趋势,我们以前六个受试者为例绘制图像观察motor-UPDRS评分变化进展,从图2中我们可以看出前六名帕金森病患者在六个月左右的时间内帕金森病进展比较缓慢,但整体呈现上升趋势,这与帕金森病的发病特性保持一致,即进展缓慢且不可逆,同时我们发现不同患者的帕金森病严重程度差异较大。
Figure 2. Pearson correlation coefficient heat map
图2. 皮尔逊相关系数热力图
4. 变量筛选
4.1. 皮尔逊相关系数图
为了观察变量之间的关系,我们首先计算所有变量之间的皮尔逊相关系数并画出热力图。
皮尔逊相关系数的计算公式为:
(1)
从热力图的结果来看,自变量之间如 Jitter类内变量之间的相关系数都在0.8以上,Shimmer类内的变量之间的皮尔逊相关系数也在0.6以上,如果全部将变量引入模型中,可以会出现冗余变量,增加模型的计算量,降低模型的计算效率和预测精度。最后一列显示的结果显示自变量与因变量之间不是简单的线性相关,可能存在更复杂的非线性关系。
4.2. Lasso变量选择
为了确保模型简洁且高效,我们首先采用了Lasso方法进行变量选择,旨在避免引入过多的冗余变量。Lasso变量选择通过施加对回归系数的约束,促使部分系数缩减至零,从而有效筛选出对因变量有显著影响的自变量,实现模型的精简。这种简化不仅有助于我们对数据内在规律进行认知,还显著降低了模型的计算复杂度。
在实际操作中,我们利用了Python软件中的lassoCV包,通过10折交叉验证的方法对数据进行Lasso回归变量选择。我们绘制了变量随不同lambda取值变化的系数压缩图以及均方误差随lambda变化的交叉验证曲线图(如图3)。结果显示当lambda取值为0.0212时,均方误差达到最小值。
Figure 3. Compression plot of LASSO coefficient (left) and variation plot of mean square error (right)
图3. Lasso系数压缩图(左)与均方误差变化图(右)
Figure 4. Coefficient compression with lasso
图4. Lasso系数图
在采用交叉验证得到的最优参数λ建立的Lasso回归模型中,我们得到了各变量的系数图。从图4中可以清晰地观察到,“Jitter (%)”“Jitter: PPQ5”“Shimmer”以及“Shimmer (dB)”这四个变量的系数被最显著地压缩至0,这表明这些变量在模型中的贡献极小,对于预测目标变量motor-UPDRS几乎无影响。为了提升模型的效率和准确性,我们决定在后续的建模过程中将这四个变量从自变量中排除。因此,我们将仅保留剩余的15个变量作为自变量,并基于这些变量建立不同的回归模型对motor-UPDRS进行预测,并比较不同模型的性能。
5. 模型构建
从第四章的分析我们可以看出自变量与因变量之间并不是简单的线性相关,此时传统的线性回归模型对数据的拟合效果并不适合,因此我们考虑高斯过程回归、支持向量回归、随机森林和XGBoost这四个模型作为候选模型,因为其在处理非线性关系时各有优势,我们通过比较常用的评价指标来选取合适的模型,并通过调整参数来获得最佳的模型性能。
5.1. 评价指标
对于motor-UPDRS的预测属于回归问题,我们事先按照4:1的比例划分训练集和测试集,常用以下几个评价指标来衡量模型的预测效果,以选择最优的模型。
1. 均方根误差(Root Mean Squared Error, RMSE)
(2)
MSE是预测值与真实值之间差异的平方的均值的平方根,RMSE的值越小,表示模型的拟合程度越好。
2. 平均绝对误差(Mean Absolute Error, MAE):
(3)
MAE是预测值与真实值之间差异的绝对值的平均值。MAE值越小,表示模型的预测结果与实际值的差异也越小。与RMSE相比,MAE对异常值不敏感。
3. 决定系数(Coefficient of Determination, R2):
(4)
用于衡量模型的拟合优度,它表示因变量的方差中有多少可以被自变量解释。取值范围:从0到1,数值越接近1表示模型对目标变量的拟合程度越好。
5.2. 高斯过程回归
5.2.1. 基本原理
高斯过程回归(GPR)是高斯过程在回归方面的应用。作为一种非参数贝叶斯模型,其基本原理是通过使用高斯过程作为先验分布,结合贝叶斯推断(Bayesian inference)来求解回归问题。GPR中给定输入值和对应的输出值,我们可以假设输入值与输出值之间的关系表示
(5)
其中
,
为噪声因子且服从正态分布
。
由此可以得出观测值
的先验分布为
(6)
对于新的观测点
,观测值
和预测值
的联合先验概率分布为:
(7)
根据贝叶斯定理,可以推导出
的后验分布也是一个高斯过程,通过计算后验分布的条件分布得到:
(8)
其中:
(9)
(10)
5.2.2. 高斯过程回归模型
鉴于不同的核函数适应不同的数据特性和问题需求,合理选择核函数能够显著提升模型的预测性能。常用的核函数有常数核,RBF核,Matern类核和有理二次核。我们基于python中的Gaussian Process Regressor分别使用不同的核函数建立模型,比较不同评价指标值,选出最优模型。
Table 2. GPR model evaluation index results with different kernel functions
表2. 不同核函数的GPR模型评价指标结果
评价指标 核函数 |
MAE |
RMSE |
常数核 |
6.91 |
8.03 |
RBF核 |
5.61 |
8.52 |
Matern核 |
3.67 |
5.77 |
有理二次核 |
3.40 |
4.63 |
从表2中的结果我们可以看出,采用有理二次核函数(Rational Quadratic)建立的模型MAE和RMSE均最小。
5.3. 支持向量回归
5.3.1. 基本原理
支持向量回归(Support Vector Regression, SVR)是一种强大的机器学习算法,是支持向量机在解决回归问题上的应用。其核心思想是在一个高维空间中寻找一个最优的超平面,这个超平面能够最大限度地拟合训练数据点,同时保持与数据点之间的最大间隔。
与传统的回归模型不同,支持向量回归不直接根据预测值
与真实值
之间的绝对误差来计算损失,它通过引入ε-insensitive损失函数,相当于在
周围创建了一个宽度为
的间隔带。这种方法允许预测值在真实值的一个ε范围内波动而不计入损失,这使得模型对噪声和异常值具有较好的鲁棒性。
SVR对于非线性问题有很强的处理能力,我们可以通过引入核函数将原始数据映射到高维空间,使得数据在高维空间内变得线性可分。
此时SVR的解为:
(11)
这里
为核函数。
5.3.2. 支持向量回归模型
我们在python中使用SVR建立回归模型,主要有以下三个参数:① 核函数(Kernel Function):将输入数据映射到高维空间,以便更好地处理非线性问题。② 惩罚参数C:控制模型的复杂度和训练误差之间的平衡。③ gamma参数:对于某些核函数(如RBF核)来说,gamma参数决定了样本点对决策边界的影响程度。
由于径向基函数(Radial Basis Function,也称为高斯核)适用于非线性关系的数据,并能较好地处理大多数问题。这里我们采取RBF核函数建立SVR模型,并采用网格搜索GridSearchCV对其余两个参数进行5折交叉验证选择最优参数组合。我们设置参数搜索范围如表3所示,得到的最优参数组合为C = 10,gamma = 1。
Table 3. SVR parameter settings
表3. SVR参数设置
参数名称 |
搜索范围 |
最优参数 |
C |
0.1, 1, 10, 100 |
10 |
gamma |
0.001, 0.01, 0.1, 1,10 |
1 |
利用最优参数建立SVR模型对motor-UPDRS评分进行预测,得到对应的评价指标结果如表4所示。
Table 4. SVR evaluation index results
表4. SVR评价指标结果
评价指标 |
值 |
RMSE |
4.42 |
MAE |
3.03 |
R2 |
0.69 |
5.4. 随机森林
随机森林算法(Random Forest, RF)是一种基于集成学习中的Bagging思想的机器学习算法。其核心思想是通过构建多棵互不关联的决策树,在随机选择的子样本上选择部分特征独立进行训练,并通过平均或加权平均来集成它们的预测结果进行回归预测[12]。
在随机森林算法中训练每棵决策树时,都会从训练集中通过自主抽样法随机选择一部分样本,并随机选择一部分特征,递归地选择最佳划分特征,将数据集划分为不纯度最小的子集,直到达到提前设置的停止条件为止。这样可以确保每棵决策树都是在不同的样本集和特征子集上进行训练的,相比于单一的决策树模型,随机森林的模型更加多样,从而具有更高的预测精度、更好的泛化能力和更强的鲁棒性。
随机森林的参数设置可以影响模型的复杂度,通过调优参数,我们希望能够找到一个平衡点,使得模型既可以充分地习得训练集的信息,又可以较好的拟合新数据。表5中首先列出了随机森林中几个重要参数及其含义。我们通过python软件下的Random Forest Regressor建立随机森林回归模型,传统的网格搜索需要遍历所有可能的参数组合,这一过程通常需要耗费大量的时间和计算资源。为此我们基于贝叶斯优化的原理建立代理模型和概率模型,在有限的试验次数内更聪明地选择下一个尝试的参数组合,以更高效地逼近全局最优解。这一过程可以通过Optuna基于贝叶斯优化的超参数优化框架库来实现。
Table 5. RF parameter settings
表5. RF参数设置
参数名称 |
参数说明 |
搜索范围 |
最优参数 |
n_estimators |
树的个数 |
(10, 1000) |
850 |
max_depth |
树的最大深度 |
(2, 30) |
22 |
max_features |
决策树随机选择的特征数量 |
“auto”, “sqrt”, “log2” |
auto |
min_samples_split |
拆分节点所需的最小样本数 |
(2, 10) |
2 |
min_samples_leaf |
叶节点所需的最小样本数 |
(1, 10) |
1 |
通过100次参数寻优过程我们得到最优的参数组合,并计算得到表6所示的随机森林模型评价指标结果。
Table 6. RF evaluation index results
表6. RF评价指标结果
评价指标 |
值 |
RMSE |
1.28 |
MAE |
0.59 |
R2 |
0.97 |
5.5. XGBoost
XGBoost是一种基于梯度提升框架的集成学习算法,它是梯度提升决策树(Gradient Boosting Decision Tree,简称GBDT)的一个优化变种[13]。其基本原理是通过集成多个决策树(也称为“弱学习器”)来构建一个强大的预测模型。XGBoost采用了梯度提升算法的思想,通过迭代地添加新的决策树来不断改善模型的预测能力。在每一次迭代中,XGBoost都会根据之前所有树的预测结果与目标值之间的差异(即残差)来训练一棵新的树,并基于这个新树来更新整体的预测。这个过程会一直持续,直到模型在训练数据上的表现达到预定的要求或迭代次数达到预设的最大值。同时为了控制模型的复杂度并防止过拟合,XGBoost还在目标函数中加入了正则化项。目标函数大致可以写为:Objective = L(θ) + Ω(θ),其中L(θ)是误差函数,它体现了模型拟合数据的能力;Ω(θ)是正则化项,用于惩罚复杂模型的参数,解决过拟合问题。
在XGBoost中,模型的超参数主要有树的结构和正则化项等。表7中首先列举了一些XGBoost模型中相对重要的参数以及参数说明,与随机森林模型的建立过程类似,我们在Python中首先使用xgboost库来建立回归模型,设置参数搜索范围,通过optuna框架下的交叉验证和贝叶斯优化等技术,我们可以找到最优的参数组合。
Table 7. XGBoost parameter settings
表7. XGBoost参数设置
参数名称 |
参数说明 |
搜索范围 |
最优参数 |
reg_alpha |
L1正则化的值 |
(1e−2, 10.0) |
0.01 |
reg_lambda |
L2正则化的值 |
(1e−2, 10.0) |
0.2 |
colsample_bytree |
队列采样率 |
(0.5, 1) |
0.9 |
learning_rate |
学习率 |
0.01, 0.05, 0.1, 0.2, 0.3 |
0.1 |
subsample |
样本采样率 |
(0.5, 1) |
0.9 |
max_depth |
最大树深 |
(3, 10) |
9 |
min_child_weight |
最小叶节点的权重 |
(1, 10) |
6 |
表8是最优参数下的模型预测结果:
Table 8. XGBoost evaluation index results
表8. XGBoost评价指标结果
评价指标 |
值 |
RMSE |
1.18 |
MAE |
0.63 |
R² |
0.98 |
5.6. 模型比较与特征重要性分析
我们将建立的四种回归模型的预测结果进行整合,同时与多元线性回归做比较,确定在预测帕金森病(PD)患者motor-UPDRS评分方面最适合的方法,对比结果如图5所示。
Figure 5. Comparison of the results of different regression models
图5. 不同回归模型评价指标结果比较
经过对结果的深入分析,本文所探讨的四种模型在预测帕金森病患者motor-UPDRS评分方面,相较于传统的多元线性回归模型,均展现出了显著的优势。表明这四种模型在处理非线性问题上具有较强的拟合能力。
Figure 6. XGBoost prediction results and true scatter plots
图6. XGBoost预测结果与真实值散点图
Figure 7. XGBoost feature importance diagram
图7. XGBoost特征重要性图
在四种模型中,基于集成学习算法的随机森林和XGBoost模型在总体性能上优于单一模型,即高斯过程回归和支持向量回归。具体来说,XGBoost模型在RMSE、MAE和R2这三个关键评价指标上均表现出色。其中,其MAE值为0.63,虽略高于随机森林的0.58,但两者之间的差距较小。而在RMSE和R2方面,XGBoost的数值分别为1.18和0.98,优于其他模型。从图6也可以看出预测值与真实值的误差 不大,因此我们可以认为在本研究中XGBoost是预测帕金森病患者的motor-UPDRS的最佳选择。
通过图7的特征重要性评分图,我们可以观察到年龄这一特征变量,在帕金森病预测的多个因素中占据首位。鉴于帕金森病作为一种在中老年人群中常见的神经性疾病,其发病概率和严重程度与年龄之间存在着显著的相关性。通常而言,随着年龄的增长,个体罹患帕金森病的风险会逐步上升,且病情可能会随时间的推移而逐渐加重。性别因素位列第二,这提示我们不同性别对帕金森病严重程度可能具有一定程度的影响,这可能是由于生理结构、内分泌差异或遗传因素等多种因素共同作用的结果。此外,在语音特征中,“DFA”“Jitter (DDP)”和“Jitter (Abs)”这三个特征指标排名较为靠前。表明这三种特定的语音信号特征对于预测motor-UPDRS评分具有较为重要的作用。
6. 总结与展望
帕金森病作为一种慢性神经系统疾病,严重影响患者的日常生活质量。因此,准确预测患者的运动功能评分对于制定个性化的治疗方案和评估治疗效果具有重要意义。
本文基于机器学习方法,利用UCI中的帕金森病远程语音数据集对帕金森病患者的运动功能评分(motor-UPDRS)进行了预测研究。首先通过lasso变量选择方法筛选出对预测结果影响较大的15个特征。随后,分别构建了四个不同的机器学习模型:包括高斯过程回归、支持向量回归、随机森林和XGBoost,并对这些模型进行训练调参优化和性能评估。最终结果显示,XGBoost模型在预测帕金森病患者的运动功能评分方面表现最佳,具有较高的预测精度和稳定性。随后本文基于表现较好的XGBoost模型进行了特征重要性分析,揭示了不同特征对motor-UPDRS预测结果的影响程度,为帕金森病的研究提供了有价值的参考。
后续我们还可以采用更多有效的特征选择方法和预测算法,以提高预测精度和泛化能力。同时,我们也可以选择帕金森病的其他相关指标,如认知功能、生活质量等,为患者提供更为全面和个性化的评估与治疗方案。