CN102930163A - 一种2型糖尿病风险状态判定方法 - Google Patents

Abstract

本发明涉及一种2型糖尿病风险状态判定方法，属于生物医学技术领域。本方法使用EM聚类和C4.5分类提取影响2型糖尿病发病的关键属性，并将人群根据自身特点划分到不同聚类簇；根据关键属性提取结果和Logistic回归方法，结合代谢综合征和高危人群判定标准判定个人或人群风险状态。在获得影响2型糖尿病的关键属性的同时，考虑人群不同特点构造不同风险状态判定模型，提升模型的普适性和实用性，并满足实时处理及移动计算的需求；通过量化分析方法评价关键属性的作用程度，具有复杂度低、准确率高、划分风险状态详细的特点。

Description

一种2型糖尿病风险状态判定方法

技术领域

本发明涉及一种2型糖尿病风险状态判定方法，属于生物医学技术领域。

背景技术

随着社会经济发展，生活方式改变以及人口老龄化等危险因素迅速增加，2型糖尿病及其并发症患病率迅速上升。据2008年调查结果显示，在20岁以上的成人中，年龄标化的糖尿病的患病率为9.7%，糖尿病前期的比例更高达15.5%。因此，早期对糖尿病危险人群采取预防措施，做到未病先防是控制糖尿病患病率的有效途径。糖尿病风险状态判定工具作为人群筛查的重要工具，正越来越受到医学人员的重视。糖尿病风险状态判定主要倾向于对患病风险状态的评估，是指对个人所处的风险状态所作的分析和推断，是风险决策的基础，糖尿病风险状态判定模型能有效并准确预测个人发生疾病的风险状态，可以减轻医学人员的工作量，同时也可以准确的针对高风险人群采取针对性的预防措施。

2型糖尿病风险状态判定需要解决2个基本问题：1.如何提取与2型糖尿病发病紧密相关的关键属性，增强风险判定的普适性；2.如何根据关键属性合理的判定个人发病风险状态。综观现有2型糖尿病风险状态判定方法，通常使用的方法有：

1.在关键属性提取方面：

根据算法原理一般分为过滤法和包装法，过滤法一般包括以下方法：

(1)Relief方法：该方法是根据统计相关性标准提取关键属性，根据特征值的区分能力来评价特征的相关度，即关键属性应该使同类的样本接近，而使不同类的样本之间远离，基本思想是：对样本进行抽样，在根据抽取的样本与同类、不同类的两个最接近样本的差异计算相关度，从而确定每个属性的不同权重。Relief方法选择的关键属性相关性较强，且能够处理离散和连续属性，但该方法不能消除冗余属性，且由于计算样本间距离将产生较大的时间开销，无法满足高维数据关键属性提取的时间性能要求。

(2)主成分分析方法：该方法研究各个属性之间的相关关系，将原来一组具有一定相关性的属性，通过变换成为一组新的属性集合作为关键属性，通过这种变换达到用较少的新属性代替原来较多属性的目的，并且使新属性尽可能多保留原来较多属性反映的信息。但是主成分分析法涉及到特征方程等矩阵运算，算法的时间消耗不能满足高维数据关键属性提取的要求。

(3)粗糙集方法：该方法是在保持属性的分类能力下，不断筛除冗余属性从而获得关键属性集合。粗糙集方法一般根据差别矩阵、属性重要度或JOHNSON约简方法，通过属性依赖程度的不同消除对分类结果影响较低的属性，达到提取关键属性的目的。虽然该方法可以有效地删除不相关属性，但没有考虑噪声数据的影响，并且计算效率低。

(4)信息熵方法：该方法主要用于信息理论中分析信息不确定度，也可用作属性作用度的评价，即关键属性提取。其基本思想是根据计算信息增益等方法划分数据，并从新计算划分后的数据增益，典型方法有ID3和C4.5，但该方法时间复杂度高。

(5)遗传算法：该方法将解表示为以二进制串编码的“染色体”，在执行算法前，给出假设解的“染色体”，然后把这些假设解置于具体问题也即“环境”中，按照一定原则从中选择出较适应环境的编码串模拟生物遗传过程的复制、交叉、变异产生更适应环境的新一代。依此进化，最后就会逐渐收敛到最适应环境的一个编码串上，也即最优解。通过该过程即可提取关键属性。但该方法需要不停迭代计算，算法时间复杂度高，在关键属性提取中很少应用。

包装法：该方法将学习算法作为测试用的黑盒子，利用相关的学习算法对属性子集进行评价，其主要思想是用训练数据和相应的学习算法训练一个分类模型，然后用测试数据来评估这个分类器的分类准确率，通过迭代提取关键属性，同时能够发现比较适合的学习算法以及算法的相关参数设定值。包装法优点是对学习算法的支持度高，缺点是该模型需要耗费大量的时间进行学习和训练，时间复杂度高，效率低，不适用于学习算法经常变动的情况。

2.在风险状态判定方法方面：

(1)芬兰糖尿病风险评分法：芬兰风险评分问卷（FINDRISC）由Lindstrom于2003年提出，是第一个通过队列研究获得的糖尿病风险评估模型，被认为是最具权威性的2型糖尿病风险评估工具，在全球应用广泛，可在无需医务人员的协同帮助下，实现自我评分预测。芬兰风险评分模型采用两组不同随访人群数据作为风险评估研究的数据源，该方法通过对随机抽样人群随访10年，并对关于年龄、体质指数(BMI)、腹围腹臀围比、血压、家族史、饮食和运动习惯等问题评分，每个问题均有固定评分标准，根据评分细则，个人可通过计算每个相应危险因素的危险系数分值，并各项相加即可得到个人风险评分，其分值越高，患糖尿病风险越大。在风险评分模型中，个人风险评分值范围在0~20分，通过临床队列分析，大于等于9分者则需要进一步诊断和检查。经过1987年和1992年这两个队列研究，模型的合理性得到其灵敏性是78%和81%，特异性是77%和76%，阳性预测值是0.13和0.05，AUC为0.85和0.87，结果显示该风险评分模型具有较好的预测效果。

(2)多因素模型法：该方法以是否患有糖尿病为因变量，关键属性为自变量，进行多元回归（Logistic回归或Cox回归），获得回归系数，然后将回归系数转化成相应的风险分值，建立模型，计算个人总风险分值，并根据ROC曲线获得风险分值的判断点。该方法通过比较个人总风险分值与判断点的大小，判定个人处于糖尿病期或糖尿病前期。多因素模型法在不同人群中经过验证，具有一定准确性，但该方法仅能判定人群是否处于糖尿病期或糖尿病前期，并不能判定其他风险状态。

(3)单因素加权分数法：该方法采用Logistic回归进行建模，将各关键属性的OR值转化成相应的风险分值，建立模型并计算总的风险分值，该方法按照总风险分值将糖尿病风险状态划分为5个状态。但该方法未考虑人群特点，使得建立模型没有普适性，不同国家和不同地区的模型也大不相同。

(4)HCI的糖尿病风险评估方法：是由吴海云、潘平等人于2007年提出的，是针对我国成年人糖尿病风险程度评估的模型之一，该方法根据多学科专家小组提出我国成年人糖尿病发病的主要风险因素及其相对风险度，建立的可通过病史及生活问卷资料评估个体糖尿病发病风险的计算方法，应用哈佛癌症风险指数计算公式，根据个人风险因素的相对风险度以及同性别年龄组的相对风险度计算个人患病的相对风险度，该模型可用于评估我国成年人个体糖尿病发病风险，提示个体不同风险因素对其糖尿病发病风险的影响。该方法建立的评估模型方法可以通过大型人群队列研究，经多因素回归分析得出某一疾病的相关预测变量，模型将患病风险程度划分为5个等级，具有一定的借鉴意义，且便于实现，并用于网络的健康管理***及社区预防保健等机构使用，但是该方法并未给出相关准确性的验证，所以该方法在实际应用的价值还有待商榷。

综上所述，对于糖尿病风险状态判定问题，现有方法选取的关键属性均不相同，没有代表性；同时，现有风险状态判定方法没有考虑人群内在特点，建立模型不具有普适性，对于不同人群的判定效果不甚理想。

发明内容

本发明的目的是解决2型糖尿病风险状态判定问题，提出一种基于关键属性提取和logistic回归的2型糖尿病风险状态判定方法，其中关键属性提取采用EM聚类和C4.5分类实现对影响2型糖尿病发病的关键属性提取。

本发明的设计原理为使用EM聚类和C4.5分类提取影响2型糖尿病发病的关键属性，并将人群根据自身特点划分到不同聚类簇；根据关键属性提取结果和Logistic回归方法，结合代谢综合征和高危人群判定标准判定个人或人群风险状态。在获得影响2型糖尿病的关键属性的同时，考虑人群不同特点构造不同风险状态判定模型，提升模型的普适性和实用性，并满足实时处理及移动计算的需求。

本发明的技术方案是通过如下步骤实现的：

步骤1，将N个被评估对象作为样本集S，其中每个对象包含M个影响2型糖尿病发病的关键属性；将关键属性作为列，不同样本对应的属性值作为行，建立样本集S的矩阵表示[s_(a+c)b]；并采用EM聚类方法对样本集S进行聚类，得到k个聚类簇。具体实现方法为：

步骤1.1，首先将N个对象根据是否含有缺失数据分为两个样本集：完全数据集X和缺失数据集Y。

所述完全数据集X为M个关键属性对应的数据均无缺失的对象集合。集合中所有对象关键属性的数据构成矩阵[x_ab]，其中矩阵的行数a表示完全数据样本对象，列数b表示完全数据样本各属性对应的数据，b=1,2,...,M。

所述缺失数据集Y为M个关键属性对应的数据中存在一个或多个缺失的对象集合。集合中所有对象的关键属性数据形成一个矩阵[y_cb]，行数c表示缺失数据样本对象，列数b表示缺失数据样本各属性对应的数据。

步骤1.2，设置聚类簇个数k（k≤N），由于样本集S中任意第i个样本n_i（1≤i≤N），均服从k个聚类簇的混合概率分布：

$n_i~p\left(n_i\right|\mathrm{\Θ})=\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\π}}_{j}p\left({\mathrm{\θ}}_i\right|{\mathrm{\θ}}_j)$

其中Θ＝(π₁，π₂，...，π_k，θ₁，θ₂，...，θ_k)表示k个聚类簇的混合概率分布参数，θ_j表示第j个聚类簇的概率分布参数，π_j表示n_i来自第j个聚类簇的可能性，j＝1,2，...,k，π₁+π₂+...+π_k＝1。

因此，设定一组参数值 ${\mathrm{\Θ}}^0=({\mathrm{\π}}_1^0,{\mathrm{\π}}_2^0,...,{\mathrm{\π}}_k^0,{\mathrm{\θ}}_1^0,{\mathrm{\θ}}_2^0,...,{\mathrm{\θ}}_k^0)$ 作为样本集S的初始混合概率分布参数估计值。

步骤1.3，将步骤1.2给定的初始混合概率分布参数估计值Θ⁰代入缺失数据集Y，获得y_cb的后验分布概率为：

$p\left(y_{\mathrm{cb}}\right|x_{\mathrm{ab}},{\mathrm{\Θ}}^0)=\frac{{\mathrm{\π}}_k^0{p}_{y_{\mathrm{cb}}}\left(x_{\mathrm{ab}}\right|{\mathrm{\θ}}_k^0)}{{\mathrm{\Σ}}_{k=1}^M{\mathrm{\π}}_k^0{p}_k\left(x_{\mathrm{ab}}\right|{\mathrm{\θ}}_k^0)}$ 其中x_ab∈X，y_cb∈Y

由于样本集间各样本的独立性，得到数据集Y的后验分布函数为：

$p\left(Y\right|X,{\mathrm{\Θ}}^0)=\underset{b=1}{\overset{M}{\mathrm{\Π}}}p\left(y_{\mathrm{cb}}\right|x_{\mathrm{ab}},{\mathrm{\Θ}}^0)$

步骤1.4，根据步骤1.3获得的后验分布函数，对完全数据的对数似然函数ln L(Θ|X，Y)关于缺失数据求期望，并将得到M个对应每一列的期望值，用于代替缺失数据集Y中相应列数的各个缺失数据，得到新样本集Y＇。

所述ln L(Θ|X，Y)为完全数据关于缺失数据的对数似然函数，其中

$\ln L\left(\mathrm{\Θ}\right|X,Y)=\ln p(x,y|\mathrm{\Θ})=\underset{b=1}{\overset{M}{\mathrm{\Σ}}}\ln p\left(x_{\mathrm{ab}}\right|y_{\mathrm{cb}})p\left(y_{\mathrm{cb}}\right)$

步骤1.5，根据步骤1.4计算获得的样本集Y＇以及完全数据集X，重新计算样本集S的最大似然参数Q(Θ,Θ⁰)：

Q(Θ,Θ⁰)=∑ln(L(Θ|X，Y)p(Y|X,Θ⁰))

步骤1.6，将Q(Θ,Θ⁰)极大化，得到Θ¹，满足Q(Θ¹，Θ⁰)=max Q(Θ，Θ⁰)，用Θ¹代替Θ⁰，代入步骤1.3。

步骤1.7，执行步骤1.3至步骤1.6的迭代过程α次，直至||Q(Θ^α+1,Θ^α)-Q(Θ^α,Θ^α-1)||<ε，停止迭代。最终得到的Θ^α即为k个聚类簇的混合概率分布参数估计值Θ。其中ε为依据聚类准确度设定的精度值。

步骤1.8，利用混合概率分布参数估计值Θ，分别计算各样本对象属于聚类簇j的后验条件概率密度值，即各对象隶属于每个簇的隶属概率。按照样本与聚类簇的隶属概率最大的原则，把样本集S中的各个样本分入k个聚类簇。

步骤2，对步骤1得到的每个聚类簇分别采用C4.5分类方法进行决策树训练，获得k个决策树。

其中第j个聚类簇的决策树的具体构建过程如下：

步骤2.1，根据中国糖尿病协会提出的2型糖尿病判定准则，将聚类簇j中的对象分为患病和未患病两个类别P和Q，并且聚类簇j分别含有g个属于类别P的样本对象和h个属于类别Q的样本对象。

步骤2.2，计算所有对象分别属于类别P和Q的信息量为：

$\mathrm{Info}\left(j\right)=\mathrm{Info}(P,Q)=-(\frac{g}{g+h}*\log \frac{g}{g+h}+\frac{h}{g+h}*\log \frac{h}{g+h}).$

步骤2.3，在第b列属性数据中选取一个数据值A_β，将本列中属性值大于等于A_β的样本对象划分至类别S₁₁，小于A_β的样本对象划分至类别S₁₂，形成两个子类集合；其中S₁₁中包含e₁个属于类别P的样本对象，f₁个属于类别Q的样本对象；S₁₂中包含e₂个属于类别P的样本对象，f₂个属于类别Q的样本对象。然后计算两个子类的信息量为：

$\mathrm{Info}(A_{\mathrm{\&beta;}},j)=\underset{o=1}{\overset{2}{\mathrm{\&Sigma;}}}\frac{e_o+f_o}{g+h}*\mathrm{Info}(P,Q)$

步骤2.4，将第b列中各个属性值分别作为A_β，按步骤2.3所述的方法获得各属性值相应的信息量，并进一步得到各属性值的信息增益值：

Gain(A_β,j)=Info(j)-Info(A_β,j)

步骤2.5，按照步骤2.3至步骤2.4的方法，对每一个列求取各属性值的信息增益值，并将得到的所有信息增益值中最大值对应的属性作为决策树的根节点；将按照信息增益值中最大值对应的A_β划分的子类S₁₁和S₁₂作为划分下一层节点的数据集。

步骤2.6，根据划分出的数据集S₁₁和S₁₂中各样本对象分别属于类别P或类别Q的个数，按照步骤2.2所述的方法，分别计算数据集S₁₁和S₁₂的信息量Info(S₁₁)和Info(S₁₂)。

步骤2.7，按照步骤2.3至步骤2.6所述的方法，重复进行λ次划分，直到被划分节点的所有样本对象属于同一类别或者所有属性均划分完毕，则停止节点划分，获得（λ+1）层决策树。

步骤3，针对每个属性，根据其在步骤2获得的k个决策树中出现的层数，以及在该层出现的次数和，计算该属性的层数系数。

对于样本集S，第b列属性的层数系数L_b为：

$L_b=\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}\frac{\underset{w=1}{\overset{\mathrm{\&lambda;}+1}{\mathrm{\&Sigma;}}}\left(\frac{1}{2^w}{t}_w\right)}{\underset{w=1}{\overset{\mathrm{\&lambda;}+1}{\mathrm{\&Sigma;}}}{t}_w}$

其中t_w表示在第w层出现了t次。

将层数系数作为衡量各属性影响2型糖尿病发病的作用程度，并选取L_b>δ的B个属性作为影响发病的主要属性集合U(U₁,U₂，...,U_B)，其中1<B＜M，δ为主要属性作用程度的阈值。

步骤4，将样本对象属于类别P记为P_d=1，属于类别Q记为P_d=0，将P_d作为因变量，步骤3获得的主要属性集合U对应的所有数据U_vz作为自变量（z=1,2,...,B，v=1,2,...,(g+h)），分别对步骤1获得的各个聚类簇中的样本对象进行Logistic回归建模，获得各主要属性对应的回归系数β_z，构建k个Logistic回归模型。

所述Logistic回归模型的形式为：

$\ln \left(\frac{P_v}{1-P_v}\right)={\mathrm{\&beta;}}_1{U}_{v1}+{\mathrm{\&beta;}}_2{U}_{v2}+...+{\mathrm{\&beta;}}_B{U}_{\mathrm{vB}}(v=1,...,g+h),$ 其中P_v为第j个聚类簇中第v个样本对象属于类别P的概率。

其中第j个聚类簇的Logistic回归模型的具体构建过程如下：

步骤4.1，获得第v个样本是否属于类别P的条件概率函数为：

$P_v\left(P_d\right|U_{v1},U_{v2},...,U_{\mathrm{vB}})=\frac{\exp \left({\mathrm{\&Sigma;}}_{z=1}^B{\mathrm{\&beta;}}_z{U}_{\mathrm{vz}}\right)}{1+\exp \left({\mathrm{\&Sigma;}}_{z=1}^B{\mathrm{\&beta;}}_z{U}_{\mathrm{vz}}\right)}$

步骤4.2，根据步骤4.1获得各个样本对象的条件概率函数，计算聚类簇j的最大似然函数，并求对数，转换为对数似然函数：

$l(P,U;{\mathrm{\&beta;}}_1,...,{\mathrm{\&beta;}}_B)=\underset{v=1}{\overset{g+h}{\mathrm{\&Sigma;}}}\left(P_d\underset{z=1}{\overset{B}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_z{U}_{\mathrm{vz}}\right)-\underset{v=1}{\overset{g+h}{\mathrm{\&Sigma;}}}\log [1+\underset{z=1}{\overset{B}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_z{U}_{\mathrm{vz}}]$

步骤4.3，对步骤4.2获得的对数似然函数中的每一个β_z求偏导，并令偏导数等于0，获得B个对数似然方程，组成对数似然方程组，求解方程组得到自变量U_vz对应的回归系数估计值β_z，建立Logistic回归模型。

步骤5，根据步骤4获得的Logistic回归模型，计算各聚类簇中的样本对象属于类别P的概率即各样本对象患病的概率，并将各个聚类簇中的样本对象按照年龄段和性别划分R组，计算每个聚类簇的各组样本对象的相对发病概率RR。通过ROC曲线获得人群相对发病概率阈值，然后比较个人相对发病概率RR值与阈值的大小，并结合代谢综合征判定方法，判定未患病的各样本对象的风险状态。将未患病样本对应的风险状态划分为四个等级：无风险、低风险、中风险和高风险。

所述第v个样本对象的相对发病概率RR=样本v的患病概率/基线发病概率。其中患病概率为根据该样本对象对应的Logistic回归模型，计算得到的样本v属于类别P的概率；基线发病概率为将第r（r=1，...,R）组中所有样本对象相应的关键属性集合U的数据均值带入相应Logistic回归模型计算得到的概率值。其中每个聚类簇获得R个基线发病概率值。

所述ROC曲线为将各样本对象相对发病概率RR中的不同值作为阈值，将各样本对象划分为患病P＇和未患病Q＇两类，并与该样本对象在步骤2.1中按照中国糖尿病协会准则判定结果对比，将该样本对象同时属于P和P＇的作为确认患病的对象，将该样本对象同时属于Q和Q＇的作为排除患病的对象。计算所有样本对象的灵敏度和特异度，并将灵敏度为纵坐标，(1-特异度)为横坐标，绘制出的曲线，并选取[灵敏度+(1-特异度)]最大的RR值作为相对发病概率的阈值。

其中，所述灵敏度=确认患病对象个数/属于类别P的对象个数；特异度=排除患病的对象个数/属于类别Q的对象个数。

有益效果

相比于粗糙集方法、遗传算法，本发明采用EM聚类和C4.5分类的关键属性提取技术的时间复杂度低。

相比于Relief方法、主成分分析方法，本发明采用的关键属性提取技术在保证高准确率的同时，通过量化分析方法评价关键属性的作用程度。

相比于芬兰糖尿病风险评分法、多因素模型法，本发明采用的关键属性提取和Logistic回归的风险状态判定方法大多采用介入性属性，具有准确率高、划分风险状态详细的特点，且判定结合代谢综合征等糖尿病前期症状，使得判定风险程度更为准确。

相比于单因素加权分数法、HCI的糖尿病风险评估方法，本发明通过EM聚类将不同人群具有不同特性问题考虑在内，针对不同人群特点提取不同关键属性，并通过医学人员使用证明，具有普适性和实用性。

附图说明

图1为本发明的风险状态判定原理图；

图2为具体实施方式中17946条横截面数据源关键属性提取流程图；

图3为具体实施方式中对17946条横截面数据进行EM聚类，且聚类簇为3时第一个聚类簇构建决策树的结果图；

图4为具体实施方式中相对发病概率ROC曲线图。

具体实施方式

为了更好的说明本发明的目的和优点，下面结合附图和实施例对本发明方法的实施方式做进一步详细说明。

分别以2001年2月至9月从北京市西城区和海淀区科研院所17946人整群抽样横断面数据、59839条横截面数据和2288条中科院01-07年随访回访数据作为输入，设计并部署3项测试：(1)针对17946条横截面数据的影响2型糖尿病发病的主要属性提取测试；(2)针对59839条横截面数据的风险状态判定可行性测试；(3)针对2288条随访数据的风险状态判定有效性测试。

下面将对上述3个测试流程逐一进行说明，所有测试均在同一台计算机上完成，具体配置为：Intel双核CPU（主频1.8G），1G内存，WindowsXP SP3操作***。

在上述1测试中，均使用EM聚类和C4.5分类的关键属性提取方法。关键属性提取方法流程如图1所示。

另外，在上述2、3测试中，使用相同的EM聚类的混合概率分布和Logistic回归模型。风险状态判定方法流程如图2所示。

1.针对17946条横截面数据的关键属性提取方法测试

该数据源来自于2001年2月至9月从北京市西城区和海淀区科研院所17946人整群抽样横断面健康调查，包含101维属性。经过预处理以后，筛选出67维13781条数据构成试验数据源。数据源处理为5部分，全部数据集合、女性数据集合、男性数据集合、有家族史数据集合和无家族史数据集合。EM聚类试验目的通过设置簇的个数使聚类的效果最好，对数似然比最小。该测试把EM聚类的簇个数选定为3和4。然后将EM聚类算法应用到5种不同的数据源上，把对应的数据源都聚成3和4类人群，即共10组聚类试验，得到5*3+5*4=35个人群，数据源划分如图3所示；然后，利用血糖门限值（6.1mmol·L^-1,5.85mmol·L^-1,5.6mmol·L^-1,5.26mmol·L^-1）对聚类以后的35个人群进行状态标定，大于门限值的标定为‘Y’，其余标定为‘N’；对标定以后的数据利用C4.5算法进行分类训练，可以得到35*4=140个分类决策树。最后，通过对决策树信息的统计，得到与II型糖尿病最相关的属性集合。每类数据源提取关键属性具体流程为：

步骤1，设置聚类的簇的个数为3或4。

步骤2，根据设置的聚类个数进行EM聚类。

步骤3，利用血糖门限值（如6.1mmol·L^-1）标定数据，将标定变量作为C4.5的分类标签变量。将聚类后获得的3个或4个数据源进行C4.5分类，获得对应数据源的决策树，获得决策树的形式如图3所示。

步骤4，统计决策树中各层出现的属性的次数，并计算各属性对应的层数系数，获得属性的作用程度。将各属性作用程度排序，提取前9个属性作为关键属性。

2.针对59839条横截面数据的风险状态判定可行性测试

该数据来源于全国各城市累积数据，经过预处理以后，筛选出14维59839条数据构成测试数据源。数据源采用提取的9维主要属性作为输入属性，采用EM聚类将59839条划分到不同的3个聚类簇中。对于每个聚类簇的Logistic回归的构建，设置步进概率的进入概率为0.05，删除概率为0.1，最大迭代次数为50，模型输入因变量为是否患病，自变量为9维关键属性，需计算参数为各属性的参数估计值、标准误差以及Wald检验值。测试3采用同样参数设置，以下不再敖述。

该测试不同聚类簇的Logistic回归模型计算个人相对发病概率，并结合各类综合征判定实现对个人的风险状态判定。通过对该数据源进行判定，并统计每一个风险状态所占比例，体现判定结果在各个风险状态所占比例的合理性以及该判定方法的可行性。具体风险状态判定流程为：

步骤1，根据9维关键属性，设定EM初始聚类个数为3，进行EM聚类。将数据源划分到3个聚类簇中。

步骤2，分别对3个聚类簇中的数据源进行logistic回归训练，将是否患病作为因变量，9维关键属性对应数据作为自变量，获得3个logistic回归模型。

步骤3，根据3个logistic回归模型分别计算对应数据源中个人相对发病概率，并通过ROC曲线选取合适判断点为2.2，ROC曲线图如图4。

步骤4，通过比较个人相对发病概率值与判断点值的大小，结合代谢综合征以及高维人群判定标准，判定个人发病风险状态。

3.针对2288条随访数据的风险状态判定有效性测试

该数据源为01-07年中科院7年随访数据，经过预处理，保留14维关键属性，共2288条。数据源通过与测试2相同判定方法判定个人发病风险状态。EM聚类的混合概率分布与Logistic回归模型均与测试2相同。

该测试分别通过对01年和07年数据源进行风险状态判定，并统计每一个风险状态状态下人群的数目，以及每一个风险状态状态下01年未发病而在07年时患病的人群数目，因为发病风险状态高的人群患病概率会大，所以该测试通过二者的比例关系验证风险状态下发病结果与风险状态程度关系的合理性。

测试结果

对于测试(1)，表1列出了本发明所提的关键属性提取方法中各属性在决策树中出现的次数。

表1  各属性在决策树中出现次数统计表

由上述结果可知，对II型糖尿病发病影响较大的九个重要发病危险因素为：血糖、年龄、高密度脂蛋白、收缩压、舒张压、胆固醇、体质系数、腹围、甘油三酯。很明显，在全部数据和不同性别的实验测试中都明显地起到比其他因素更重要的作用。该结果与临床医学经验证明，所提出的9维属性均为影响2型糖尿病发病的重要关键属性。

对于测试(2)，表2列出本发明所提出的风险状态判定方法的合理性测试，表中分别体现不同风险状态下人群的分布比例，未患病的人群中，随着发病风险状态越高，人数比例越少，各个发病危险状态所占比例比较符合人群特征。其中风险状态判定结果中高危风险下的数据占整个人群的比高于医学常识，主要是由于采集横截面数据源的原因，在该数据源中有一部分是患有代谢综合症的人群，导致这部分发病危险状态的数据所占的比例较高，从总体上证明本发明的合理性和可行性。

表2  不同风险状态人群分布

对于测试(3)，表3列出本发明所提出的风险状态判定方法的有效性测试，表中为每一个风险状态下人群的数目，以及每一个风险状态下01年未发病而在07年时患病的人群数目和比例。

表3  不同风险状态人群7年后发病人数分布

发病风险状态越高的人群，在7年以后的发病率也越高；01年时风险状态越低的人群，在07年发病的比例就越少，反之就越高；在各个人群中，各个风险状态的分布呈递减的趋势，而且每个风险状态的人群比例和实际人群的分布也一致，证明风险状态判定算法是有效的。

本发明在关键属性提取方法中，通过提取的9个关键属性，确立影响2型糖尿病发病的关键属性。通过对不同风险状态人群分布验证以及不同风险状态N年后发病分布验证，证明风险状态判定方法的合理性和有效性。

Claims (5)

1.一种2型糖尿病风险状态判定方法，其特征在于：包括以下步骤：

步骤1，将N个被评估对象作为样本集S，其中每个对象包含M个影响2型糖尿病发病的关键属性；将关键属性作为列，不同样本对应的属性值作为行，建立样本集S的矩阵表示；采用EM聚类方法对样本集S进行聚类，得到k个聚类簇；具体实现方法为：

步骤1.1，将N个对象根据是否含有缺失数据分为完全数据集X和缺失数据集Y；

步骤1.2，设置聚类簇个数k，以及样本集S的初始混合概率分布参数估计值

θ_j表示第j个聚类簇的概率分布参数，π_j表示第i个样本n_i来自第j个聚类簇的可能性，j＝1,2，...,k，π₁+π₂+...+π_k＝1；k≤N，1≤i≤N；

步骤1.3，将步骤1.2给定的初始混合概率分布参数估计值Θ⁰代入缺失数据集Y，获得y_cb的后验分布概率为：

$p\left(y_{\mathrm{cb}}\right|x_{\mathrm{ab}},{\mathrm{\&Theta;}}^0)=\frac{{\mathrm{\&pi;}}_k^0{p}_{y_{\mathrm{cb}}}\left(x_{\mathrm{ab}}\right|{\mathrm{\&theta;}}_k^0)}{{\mathrm{\&Sigma;}}_{k=1}^M{\mathrm{\&pi;}}_k^0{p}_k\left(x_{\mathrm{ab}}\right|{\mathrm{\&theta;}}_k^0)}$

其中x_ab∈X，y_cb∈Y，a表示完全数据样本对象，b表示完全数据样本各属性对应的数据，c表示缺失数据样本对象，b=1,2,...,M；

数据集Y的后验分布函数为：

$p\left(Y\right|X,{\mathrm{\&Theta;}}^0)=\underset{b=1}{\overset{M}{\mathrm{\&Pi;}}}p\left(y_{\mathrm{cb}}\right|x_{\mathrm{ab}},{\mathrm{\&Theta;}}^0);$

步骤1.4，根据步骤1.3获得的后验分布函数，对完全数据的对数似然函数ln L(Θ|X，Y)求期望，得到M个期望值，用其代替缺失数据集Y中相应列数的各个缺失数据，得到新样本集Y＇；

所述 $\ln L\left(\mathrm{\&Theta;}\right|X,Y)=\ln p(x,y|\mathrm{\&Theta;})=\underset{b=1}{\overset{M}{\mathrm{\&Sigma;}}}\ln p\left(x_{\mathrm{ab}}\right|y_{\mathrm{cb}})p\left(y_{\mathrm{cb}}\right);$

步骤1.5，根据步骤1.4计算获得的样本集Y＇以及完全数据集X，计算样本集S的最大似然参数Q(Θ,Θ⁰)：

Q(Θ,Θ⁰)=∑ln(L(Θ|X，Y)p(Y|X,Θ⁰))

步骤1.6，将Q(Θ,Θ⁰)极大化，得到Θ¹，满足Q(Θ¹，Θ⁰)=max Q(Θ，Θ⁰)，用Θ¹代替Θ⁰，代入步骤1.3；

步骤1.7，执行步骤1.3至步骤1.6的迭代过程α次，直至||Q(Θ^α+1,Θ^α)-Q(Θ^α,Θ^α-1)||<ε，停止迭代；得到的Θ^α作为k个聚类簇的混合概率分布参数估计值Θ；

步骤1.8，利用混合概率分布参数估计值Θ，分别计算各样本对象属于聚类簇j的后验条件概率密度值；按照样本与聚类簇的隶属概率最大的原则，把样本集S中的各个样本分入k个聚类簇；

步骤2，对步骤1得到的每个聚类簇分别采用C4.5分类方法进行决策树训练，获得k个决策树；

步骤3，针对每个属性，根据其在步骤2获得的k个决策树中出现的层数，以及在该层出现的次数和，计算该属性的层数系数；

对于样本集S，第b列属性的层数系数L_b为：

$L_b=\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}\frac{\underset{w=1}{\overset{\mathrm{\&lambda;}+1}{\mathrm{\&Sigma;}}}\left(\frac{1}{2^w}{t}_w\right)}{\underset{w=1}{\overset{\mathrm{\&lambda;}+1}{\mathrm{\&Sigma;}}}{t}_w}$

其中t_w表示在第w层出现了t次；

将层数系数作为衡量各属性影响2型糖尿病发病的作用程度，并选取L_b>δ的B个属性作为影响发病的主要属性集合U(U₁,U₂，...,U_B)，其中1<B＜M；

步骤4，将样本对象属于类别P记为P_d=1，属于类别Q记为P_d=0，将P_d作为因变量，步骤3获得的主要属性集合U对应的所有数据U_vz作为自变量，v=1,2,...,(g+h)，z=1,2,...,B，分别对步骤1获得的各个聚类簇中的样本对象进行Logistic回归建模，获得各主要属性对应的回归系数β_z，构建k个Logistic回归模型；

步骤5，根据步骤4获得的Logistic回归模型，计算各聚类簇中的样本对象属于类别P的概率，并将各个聚类簇中的样本对象按照年龄段和性别划分R组，计算每个聚类簇的各组样本对象的相对发病概率RR；通过ROC曲线获得人群相对发病概率阈值，然后比较个人RR值与阈值的大小，并结合代谢综合征判定方法，判定未患病的各样本对象的风险状态。将未患病样本对应的风险状态划分为四个等级：无风险、低风险、中风险和高风险；

所述相对发病概率阈值为[灵敏度+(1-特异度)]最大的RR值；

所述第v个样本对象的相对发病概率RR=样本v的患病概率/基线发病概率；其中患病概率为根据该样本对象对应的Logistic回归模型，计算得到的样本v属于类别P的概率；基线发病概率为将第r组中所有样本对象相应的关键属性集合U的数据均值带入相应Logistic回归模型计算得到的概率值；其中r=1，...,R，每个聚类簇获得R个基线发病概率值。

2.根据权利要求1所述的一种2型糖尿病风险状态判定方法，其特征在于：所述完全数据集X为M个关键属性对应的数据均无缺失的对象集合，其矩阵表示为[x_ab]；所述缺失数据集Y为M个关键属性对应的数据中存在一个或多个缺失的对象集合，其矩阵表示为[y_cb]。

3.根据权利要求1所述的一种2型糖尿病风险状态判定方法，其特征在于：所述步骤2中，第j个聚类簇的决策树的具体构建过程为：

步骤2.1，根据中国糖尿病协会提出的2型糖尿病判定准则，将聚类簇j中的对象分为患病类别P和未患病类别Q，其中P包含g个样本对象，Q包含h个样本对象；

步骤2.2，计算所有对象分别属于类别P和Q的信息量为：

$\mathrm{Info}\left(j\right)=\mathrm{Info}(P,Q)=-(\frac{g}{g+h}*\log \frac{g}{g+h}+\frac{h}{g+h}*\log \frac{h}{g+h});$

步骤2.3，在第b列属性数据中选取一个数据值A_β，将本列中属性值大于等于A_β的样本对象划分至子类S₁₁，小于A_β的样本对象划分至子类S₁₂；其中S₁₁包含e₁个属于类别P的样本对象，f₁个属于类别Q的样本对象；S₁₂包含e₂个属于类别P的样本对象，f₂个属于类别Q的样本对象；两个子类的信息量为：

$\mathrm{Info}(A_{\mathrm{\&beta;}},j)=\underset{o=1}{\overset{2}{\mathrm{\&Sigma;}}}\frac{e_o+f_o}{g+h}*\mathrm{Info}(P,Q);$

步骤2.4，将第b列中各个属性值分别作为A_β，按步骤2.3所述的方法获得各属性值相应的信息量，并求得各属性值的信息增益值：

Gain(A_β,j)=Info(j)-Info(A_β,j)

步骤2.5，按照步骤2.3至步骤2.4的方法，对每一列求取各属性值的信息增益值，并将得到的信息增益值中最大值对应的属性作为决策树的根节点；将按照信息增益值中最大值对应的A_β划分的子类S₁₁和S₁₂作为划分下一层节点的数据集；

步骤2.6，根据划分出的数据集S₁₁和S₁₂中各样本对象分别属于类别P或类别Q的个数，分别计算数据集S₁₁和S₁₂的信息量Info(S₁₁)和Info(S₁₂)；

步骤2.7，按照步骤2.3至步骤2.6所述的方法，进行λ次划分，直到被划分节点的所有样本对象属于同一类别或所有属性均划分完毕，则停止节点划分，获得（λ+1）层决策树。

4.根据权利要求1所述的一种2型糖尿病风险状态判定方法，其特征在于：所述Logistic回归模型的形式为：

$\ln \left(\frac{P_v}{1-P_v}\right)={\mathrm{\&beta;}}_1{U}_{v1}+{\mathrm{\&beta;}}_2{U}_{v2}+...+{\mathrm{\&beta;}}_B{U}_{\mathrm{vB}}(v=1,...,g+h),$ 其中P_v为第j个聚类簇中第v个样本对象属于类别P的概率；

其中第j个聚类簇的Logistic回归模型的具体构建过程为：

步骤4.1，第v个样本属于类别P的条件概率函数为：

$P_v\left(P_d\right|U_{v1},U_{v2},...,U_{\mathrm{vB}})=\frac{\exp \left({\mathrm{\&Sigma;}}_{z=1}^B{\mathrm{\&beta;}}_z{U}_{\mathrm{vz}}\right)}{1+\exp \left({\mathrm{\&Sigma;}}_{z=1}^B{\mathrm{\&beta;}}_z{U}_{\mathrm{vz}}\right)}$

步骤4.2，根据步骤4.1获得各个样本对象的条件概率函数，计算聚类簇j的最大似然函数，并转换为对数似然函数：

$l(P,U;{\mathrm{\&beta;}}_1,...,{\mathrm{\&beta;}}_B)=\underset{v=1}{\overset{g+h}{\mathrm{\&Sigma;}}}\left(P_d\underset{z=1}{\overset{B}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_z{U}_{\mathrm{vz}}\right)-\underset{v=1}{\overset{g+h}{\mathrm{\&Sigma;}}}\log [1+\underset{z=1}{\overset{B}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_z{U}_{\mathrm{vz}}]$

步骤4.3，对步骤4.2获得的对数似然函数中的每一个β_z求偏导；令偏导数等于0，获得B个对数似然方程，组成对数似然方程组；求解方程组得到自变量U_vz对应的回归系数估计值β_z，建立Logistic回归模型。

5.根据权利要求1所述的一种2型糖尿病风险状态判定方法，其特征在于：所述ROC曲线为将各样本对象不同的相对发病概率作为阈值，将各样本对象划分为患病类P＇和未患病类Q＇，并与根据中国糖尿病协会的2型糖尿病判定准则得到的患病类P、未患病类Q对比，将同时属于P和P＇的样本对象作为确认患病对象，将同时属于Q和Q＇的样本对象作为排除患病对象；计算所有样本对象的灵敏度和特异度，并将灵敏度为纵坐标，(1-特异度)为横坐标，绘制出的曲线；其中，灵敏度=确认患病对象个数/属于类别P的对象个数；特异度=排除患病对象个数/属于类别Q的对象个数。

Images