地贫患者的基因筛查问题
摘 要
地中海贫血(简称“地贫”)是全球广为流行、危害极为严重的遗传性溶血性疾病，全世界至少有3.45亿人携带地中海贫血的致病基因。医学上通过大人群的基因筛查来预防地贫患儿的出生。
本文应用统计学原理，对病人以及健康人的110个基因进行分析，采用Fisher判别模型建立判别标准和多元统计模型spss 软件进行筛选。
问题一，利用费希尔模型判别待测者是否患有地贫，以编号1~20地贫患者的样本，编号21~40健康人员的样本，分别作为模版建立模型，用mathlab软件求解得到待测组的患病者编号41~60个是待筛查人员的样本。
问题二，为确定“地贫”样本与“健康”样本在基因链上的区别。以及癌症样本中是否有子类。我们用1~20数据为标准化并确立相关系数矩阵，求出相关矩阵的特征值和特征向量，然后通过前m 个主成 分的累计贡献率满足 来确定贡献率矩阵，从而得出各种基因的权值，又利用初始特征值需大于 1，再运用逐步剔除法得出关键基因
关键字：地贫患者的基因 Fisher判别 筛查 相关系数矩阵
1 问题重述
化验指标能够协助医生诊断。人们到医院就诊时，诊断就诊人员是否患肾炎时通常要化验人体内各种元素含量。表B.1是确诊病例的化验结果，其中1－30号病例是已经确诊为肾炎病人的化验结果；31－60号病例是已经确定为健康人的结果。表B.2是就诊人员的化验结果。
1. 根据表B.1中的数据，提出一种或多种简便的判别方法，判别属于患者或健康人的方法，并检验你提出方法的正确性。
2. 按照1提出的方法，判断表B.2中的30名就诊人员的化验结果进行判别，判定他(她)们是肾炎病人还是健康人。
3. 能否根据表B.1的数据特征，确定哪些指标是影响人们患肾炎的关键或主要因素，以便减少化验的指标。
4. 根据3的结果，重复2的工作。
5. 对2和4的结果作进一步的分析。
2 问题分析
问题解决的关键是如何正确判断正常人与患者之间的差异，利用所给数据，可以选择用医学统计方法[1]中的判别分析法[2]进行分析。从题目给出的表中可以得出以下信息：
1) 表中分别给出正常人与患者各30组数据，每组数据各包含7种元素(Zn、Cu、Fe、Ca、Mg、K、Na)在人体中的含量。通过对这些数据进行分析，可以从中找出数据差异，根据判别法确定判别标准。利用所得判别标准，与就诊人员的化验结果比较可以判别其是否患有肾炎；
2) 采用某医院化验的30名健康人和30名肾炎患者体内的7种微量元素的含量作为基础数据，从中选取20名健康人和20名肾炎患者的数据作为样本，进行采样分析，建立判定模型，并利用剩余的 10 名健康人和 10 名患者对判定模型进行正确率检验，即可得到判定模型的正确率。
3) 利用不同的判别法，用Matlab软件编程对数据进行处理，得到相应判别标准。将就诊人员数据与判别值比较，即可对就诊人员是否患病做出明确判断。这样就解决了第二问；
4) 在第一问基础上利用EXCEL将各元素比重进行作图分析，即可直观地确定各种元素的权重，从而筛选出关键因素。这样就解决了第三问；
5) 在第四问中将第三问得出的关键元素对应的数据的在判别程序中保留，其余元素的数据删除。再利用Matlab软件对就诊人员数据进行处理，判断其是否患病。
6) 在第五问中，将第二问与第四问所得结果进行直观比较，观察其差异进行详细分析。
3 模型假设
1) 表中所给数据全部正确；
2) 除表中所给元素外，其它元素对肾炎的影响甚小；
3) 外界条件对肾炎患者影响不计。
4 符号说明
注：约定0表示正常人，1表示患者，例如x0表示正常人数据。
A：患者矩阵
B：正常人矩阵
S1：患者离差矩阵
S2：正常人离差矩阵
：正常人数据代表，即均值
：患者数据代表，即均值
：组A代表(平均值的判别值)
：组B代表(平均值的判别值)
：正常人数据协方差矩阵
：患者数据协方差矩阵
：患者数据马氏距离定义
：正常人样本马氏距离定义
：Fisher判别模型临界值
：马氏距离判别模型判别函数
5 模型建立
该问题的关键是如何判别一个人是否患有肾炎，属于分类问题，可以采用统计方法中的判别分析方法进行处理。题目只需区分病体和健康体，所以可以采用二类群体判别的方法。可选取判别模型中的Fisher判别模型、马氏距离判别模型和BP神经网络分类模型进行研究。
5.1 不等协差阵的两总体Fisher判别模型[2][3]
5.1.1 基本思想
从两个总体中抽取具有p个指标的样品观测数据，借助方差分析的思想造一个判别函数或称判别式：
y = c1x1 + c2x2 + … + cpxp
其中系数c1，c2，…，cp 确定的原则是使两组间的区别最大，而使每个组内部的离差最小。有了判别式后，对于一个新样品，将它的p个指标值代入判别式中求出y值，然后与判别临界值(或称分界点后面给出)进行比较，就可以判别它应属于哪一个总体。
5.1.2 建立过程
将患者与正常人数据分别表示成矩阵如下：
算出两组数据的代表，即平均值：
(正常人均值) (患者均值)
其中，j = 1，…，7。
做新的矩阵A1，B1及两组数据的离差矩阵S1，S2如下：
则有：
最优判别函数的系数C1，C2，…，C7为下述方程组的解：
即：
可以写出判别函数：
(S-1)
算出两组数据平均值的判别值：
则临界值为：
故若将待测样本值代入S-1式得出相应y值，可得出如下判别结果：
1) 时，若，则判别该对象患病，属于组A；若，则判别该对象正常，属于组B；
2) 时，若，则判别该对象正常，属于组B；若，则判别该对象患病，属于组A.
5.2 马氏距离判别模型[2]
5.2.1 基本思想
首先根据已知分类的数据，分别计算各类的中心即分组(类)的均值。判别准则是对任给的一次观测，若它与第i类的中心距离最近，就认为它来自第i类。
5.2.2 建立过程
设，，，分别为A、B的均值向量和协方差阵。距离定义采用马氏距离，即：
根据问题要求，将G0对应于正常人数据组B，将G1对应于患者数据组A。
首先计算X到A、B两个总体的距离，分别记为和，按距离最近准则判别归类，则可写成：
记，则有：
然后比较和的大小，按距离最近准则判别归类。
① 若正常人与患者数据的协方差相同,由得到判别函数为：
其中，判别标准是
② 若正常人与患者数据的协方差不同，即，则定义判别函数：
5.3 BP神经网络分类模型[1][5]
5.3.1 基本思想
基于误差反向传播(Back propagation)算法的多层前馈网络(Multiple - layer feed forward network ,简记为 BP 神经网络)，是目前应用最成功和广泛的人工神经网络(如下图所示)，它基于信息的正向传播和误差的反向回馈来训练整个 BP 神经网络，使得输出层误差的平方和达到最小。
5.3.2 建立过程
为避免变量各指标数值间差异太悬殊，而导致小数值被大数值所淹没，首先利用premnmx函数对各样本数据进行归一化处理，归一化处理后其值分布在[-1,1]区间内。然后将表B1的正常人和患者数据作为训练样本对网络进行学习训练，网络训练结束后，输入待判样本数据，即可通过网络获得判别结果。
6 模型检验
6.1 利用Fisher判别法的模型检验 [3]
代入数据可得出下列结果：
= -0.0108，= -0.1515，= -0.0812。
进行判别函数的检验，检验其是否有效：
计算统计量，其中，p为向量维数，p = 7：
给定显著水平a=0.05，从F分布表中查出：
所作出的判别函数有效
6.2 三种判别法的正确率检验
采用某医院化验的30名健康人和30名肾炎患者体内的7种微量元素的含量作为基础数据，从中选取20名健康人和20名肾炎患者的数据作为样本，进行采样分析，分别建立了马氏距离判定模型、 费歇尔( Fisher)判定模型和BP 神经网络分类模型 3 种判定模型，并利用剩余的 10 名健康人和 10 名患者对判定模型进行检验，得到 3 种判定模型的正确率如下：
实际结果
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
正确率
Fisher判别程序输出结果
1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1
95%
马氏距离判别程序输出结果
1 1 0 1 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0
85%
BP神经网络分类模型判别结果
实际输出值
0.9991 0.9999 0.5070 0.9159 0.9747 0.9996 0.0761 1.0000 0.9998 0.7341 0.0001 0.0032 0.0023 0.0055 0.0009 0.0008 0.0011 0.0016 0.0021 0.0118
95%
判定值
1 1 1 1 1 1 0 1 1 1 0 0 0 0 0 0 0 0 0 0
注：对于BP模型而言，输出向量是以0、1为标准对样本进行分类，趋近于0分为“0”类，趋近于1的划分为“1”类
由上表可见，正确率都达到85%以上，故建立模型正确有效。
7 模型求解
根据已建立的模型，利用Matlab软件进行编程求解，对B2中30名就诊人员化验结果判断结果如下：
① 马氏距离判别模型判别结果：
11010 10000 11000 10000 00101 00000
② Fisher判别模型判别结果：
11011 10110 11101 10010 00101 00000
③ BP神经网络分类模型判别结果(数据已进行分类处理)：
11111 00110 11100 00000 00101 01000
其中，BP神经网络分类模型判别结果误差分析图如下：
8 模型结果分析
8.1 对问题三解决
数学期望表示样本的平均值，通过正常人与患者各元素含量的期望差的绝对值与正常人的期望之比大小表示各元素对总体结果的影响程度。比值越大，表明该元素对诊断结果的决定性作用越大。将各比值进行归一化处理得到图表如下：
通过图表观察发现，首先由于K的比例值比较突出，即可选择单独对K利用程序进行对样本测试，记得测试方案A；其次，Ca的比值也相对较高，也选择可对Ca和K进行测试，即得测试方案B；最后，Mg和Fe相当，故可选择Fe、Ca、Mg、K为主要元素进行测试，可得测试方案C。
8.2 对问题四的解决
根据问题三所得测试方案，利用M软件求解如下：
Fisher判别模型判别结果
方案A
11010 00100 11000 00000 00110 00000
方案B
11111 11110 11100 10010 00111 01000
方案C
11111 11110 11100 10010 00111 01000
马氏距离判别模型判别结果
方案A
11110 11100 11011 11011 01111 01010
方案B
11011 11110 11100 10000 00100 00000
方案C
11011 10010 11000 10000 00100 00000
BP神经网络分类模型判别结果
方案A
11110 00100 11000 00000 00110 01010
方案B
11010 00110 11100 00000 00110 01000
方案C
11111 10110 11100 00000 00110 01000
8.3 对问题五的解决
对于A、B、C三种测试方案利在不同判别模型下得到的差异数据个数和相对误差如下表：
模型名称
指 标
方案A
方案B
方案C
Fisher判别模型
差异数据个数
9
5
5
相对误差
30%
16.7%
16.7%
马氏距离判别模型
差异数据个数
12
6
3
相对误差
40%
20%
10%
BP神经网络分类模型
差异数据个数
8
4
3
相对误差
26.7%
13.8%
10%
注：1、综合方案将7种元素全部考虑
2、相应模型下各方案的相对误差=相应模型下各方案结果相对综合方案结果的差异数据个数30
从表格中可以看出，在各自模型下，方案A的相对误差都比较明显，而方案C的相对误差最小，故最优应选用方案C，方案B其次，而方案A应舍弃。
从上述模型检验中可以看出，相对于用马氏距离判别模型，费歇尔(Fisher)判定模型和BP 神经网络模型判定是否患病的正确率较高。
9 模型的评价与应用
模型以判别法为核心，根据实验数据构建出判别标准，较准确地解决了题目所给问题。很好的完成了对待测就诊人员的分类(即诊断其是否患肾炎)。在解决第一问中，马氏距离判别模型、费希尔(Fisher)判别模型和BP神经网络分类模型三个模型对数据中所有相关元素的进行了综合的考虑，得到相对较好的结果。为了减少化验指标，利用关键因素，所得的简化指标模型也能相对较好地对就诊人员数据进行分类。但是马氏距离判别法和费希尔判别法处理半结构化和非结构化问题摆脱不了评判过程中的随机性和评判者主观上的不确定性，存在误判现象，为减少误判，提高判别精度我们提出向BP神经网络判别模型改进。BP神经网络判定模型，它基于信息的正向传播和误差的反向回馈来训练整个BP神经网络，使得输出层误差的平方和达到最小。BP神经网络有模拟推理和自动学习性，更接近人脑的自组织和并行处理能力，使其具有一定的智能性，在建模判定方面的能力优于其它方法。
在实际生活中，处理一些分类及判别问题，都可以采用我们的模型进行分析。
10 参考文献
【1】 赵国贞、周圣武，《关于肾炎诊断的统计方法研究》，《云南大学学报(自然科学版)》，，30(S2)：186~191；
【2】 刘焕彬等，《数学模型与实验》，北京，科学出版社，；
【3】 周玉元、周铁军，《DNA序列分类的Fisher判别法》，《湖南农业大学学报(自然科学版)》，第29卷第3期，10月；
【4】 王正林、刘 明，《精通MATLABA 7》，北京，电子工业出版社，.7；
【5】 唐万梅，《基于MATLAB的神经网络判别模型》，《数学的实践与认识》，第36卷第2期，2月。
11 附录
Matlab求解程序如下：
马氏距离判别模型主程序：
clc
clear
A=[166 15.8 24.5700 112 179 513 ;
185 15.7 31.5701 125 184 427 ;
193 9.80 25.9541 163 128 642 ;
159 14.2 39.7896 99.2 239 726 ;
226 16.2 23.8606 152 70.3 218 ;
171 9.29 9.29307 187 45.5 257 ;
201 13.3 26.6551 101 49.4 141 ;
147 14.5 30.0659 102 154 680 ;
172 8.85 7.86551 75.7 98.4 318 ;
156 11.5 32.5639 107 103 552 ;
132 15.9 17.7578 92.4 1314 1372 ;
182 11.3 11.3767 111 264 672 ;
186 9.26 37.1958 233 73.0 347 ;
162 8.23 27.1625 108 62.4 465 ;
150 6.63 21.0627 140 179 639 ;
159 10.7 11.7612 190 98.5 390 ;
117 16.1 7.04988 95.5 136 572 ;
181 10.1 4.041437 184 101 542 ;
146 20.7 23.81232 128 150 1092 ;
42.3 10.3 9.70629 93.7 439 888 ;
28.2 12.4 53.1370 44.1 454 852 ;
154 13.8 53.3621 105 160 723 ;
179 12.2 17.91139 150 45.2 218 ;
13.5 3.36 16.8135 32.6 51.6 182 ;
175 5.84 24.9807 123 55.6 126 ;
113 15.8 47.3626 53.6 168 627 ;
50.5 11.6 6.30608 58.9 58.9 139 ;
78.6 14.6 9.70421 70.8 133 464 ;
90.0 3.27 8.17622 52.3 770 852 ;
178 28.8 32.4992 112 70.2 169 ;];
B=[213 19.1 36.22220249 40.0 168 ;
170 13.929.81285226 47.9 330 ;
162 13.219.81521166 36.2 133 ;
203 13.090.81544162 98.90 394 ;
167 13.114.12278212 46.3 134 ;
164 12.918.62993197 36.3 94.5 ;
167 15.027.02056260 64.6 237 ;
158 14.437.01025101 44.6 72.5 ;
133 22.831.01633401 180 899 ;
156 135322 67471090 228 810 ;
169 8.00308 106899.1 53.0 289 ;
247 17.38.652554241 77.9 373 ;
166 8.1062.81233252 134 649 ;
209 6.4386.92157288 74.0 219 ;
182 6.4961.73870432 143 367 ;
235 15.623.41806166 68.8 188 ;
173 19.117.02497295 65.8 287 ;
151 19.764.22031403 182 874 ;
191 65.435.05361392 137 688 ;
223 24.486.03603353 97.7 479 ;
221 20.1155 3172368 150 739 ;
217 25.028.22343373 110 494 ;
164 22.235.52212281 153 549 ;
173 8.9936.01624216 103 257 ;
202 18.617.73785225 31.0 67.3 ;
182 17.324.83073246 50.7 109 ;
211 24.017.03836428 73.5 351 ;
246 21.593.22112354 71.7 195 ;
164 16.138.02135152 64.3 240 ;
179 21.035.01560226 47.9 330 ;];
X=[58.2 5.4229.7323 138 179 513 ;
106 1.8740.5542 177 184 427 ;
152 0.8012.51332176 128 646 ;
85.5 1.703.99503 62.3 238 762.6;
144 0.7015.1547 79.7 71.0 218.5;
85.7 1.094.2 790 170 45.8 257.9;
144 0.309.11417 552 49.5 141.5;
170 4.169.32943 260 155 680.8;
176 0.5727.3318 133 99.4 318.8;
192 7.0632.91969343 103 553 ;
188 8.2822.61208231 1314 1372 ;
153 5.8734.8328 163 264 672.5;
143 2.8415.7265 123 73.0 347.5;
213 19.136.22220249 62.0 465.8;
192 20.123.81606156 40.0 168 ;
171 10.530.5672 145 47.0 330.5;
162 13.219.81521166 36.2 133 ;
203 13.090.81544162 98.9 394.5;
164 20.128.91062161 47.3 134.5;
167 13.114.12278212 36.5 96.5 ;
164 12.918.62993197 65.5 237.8;
167 15.027.02056260 44.8 72.0 ;
158 14.437.01025101 180 899.5;
133 22.831.31633401 228 289 ;
169 8.030.8106899.1 53.0 817 ;
247 17.38.652554241 77.5 373.5;
185 3.9031.31211190 134 649.8;
209 6.4386.92157288 74.0 219.8;
182 6.4961.73870432 143 367.5;
235 15.623.41806166 68.9 188 ;];
Ca=cov(A);
Cb=cov(B);
CA=inv(Ca);
CB=inv(Cb);
Ua=sum(A)./30;
Ub=sum(B)./30;
for i=1:30
Xa=X(i,:)-Ua;
DA(i)=Xa*CA*Xa;
Xb=X(i,:)-Ub;
DB(i)=Xb*CB*Xb;
w(i)=DA(i)-DB(i);
end
for i=1:30
if w(i)>0
x(i)=0;
else x(i)=1;
end
end
x
Fisher判别模型主程序：
clc
clear
A=[…]; %A数组赋值与马氏模型中A的赋值情况相同%
B=[…]; %B数组赋值与马氏模型中B的赋值情况相同%
X=[…]; %X数组赋值与马氏模型中X的赋值情况相同%
Ua=sum(A)./30;
Ub=sum(B)./30;
d=(Ua-Ub);
EA=[Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;Ua;]
EB=[Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;Ub;];
A0=(A-EA);
B0=(B-EB);
S1=A0*A;
S2=B0*B;
S=S1+S2;
S0=inv(S);
C=S0*d;
m=A*C;
Ya=sum(m)/30;
n=B*C;
Yb=sum(n)/30;
Ye=(Ya+Yb)/2;
for i=1:30
Y(i)=X(i,:)*C;
if Y(i)>Ye
x(i)=1;
else x(i)=0;
end
end
x
BP神经网络分类模型程序：
clear
clc
A=[…]; %A数组赋值与马氏模型中A的赋值情况相同%
B=[…]; %B数组赋值与马氏模型中B的赋值情况相同%
X=[…]; %X数组赋值与马氏模型中X的赋值情况相同%
d=[A;B];
q=premnmx(d);
X=premnmx(X);
X=X;
p=q;
t=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
net=newff(minmax(p),[3 1],{tansig logsig},traingdx);
net.trainParam.show=50;
net.trainParam.lr=0.1;
net.trainParam.lr_inc=1.05;
net.trainParam.epochs=1000;
net.trainParam.goal=0.01;
net=train(net,p,t);
y=sim(net,p);
w=sim(net,X)
- 14 -
展开阅读全文
如果觉得《数模-化验结果判别及matlab程序 数模-化验结果判别及matlab程序.doc》对你有帮助，请点赞、收藏，并留下你的观点哦！
