两个正态随机变量均为小样本时,t—检验法可用来检验它们的数学期望是否有显著差异
在Matlab中执行命令
[h,sig,ci]=ttest2(x,y,0.01,0);
h为一个布尔值1或0只,h=0表示不可以拒绝假设,h=1表示可以拒绝假设,sig为假设成功的概率,ci为x与y均值差的置信区间。
把h=1对应的基因筛选出来作为下一步采用主成分分析进行筛选的输入数据。
最终得到具有显著性差异的基因36个。
4.2、主成分分析法,取出对致病起最终作用的基因
由t检验得出的与患病有关的基因还是过多,因此我们需要进行进一步的筛选。
利用原样本基因变量之间的相关关系,用较少的新变量代替原来较多的变量,并使这些少数变量尽可能多的保留原来较多的变量所反应的信息,这样问题就简单化了。
4.2.1主成分分析法原理:
假定有n个样本,每个样本共有p个变量,构成一个n×p阶的数据矩阵,
记原变量指标为x1,x2,…,xp,设它们降维处理后的综合指标,即新变量为z1,z2,z3,…,zm(m≤p),则
系数lij的确定原则:
①zi与zj(i≠j;i,j=1,2,…,m)相互无关;
②z1是x1,x2,…,xP的一切线性组合中方差最大者,z2是与z1不相关的x1,x2,…,xP的所有线性组合中方差最大者;zm是与z1,z2,……,zm-1都不相关的x1,x2,…xP,的所有线性组合中方差最大者。
新变量指标z1,z2,…,zm分别称为原变量指标x1,x2,…,xP的第1,第2,…,第m主成分。
从以上的分析可以看出,主成分分析的实质就是确定原来变量xj(j=1,2,…,p)在诸主成分zi(i=1,2,…,m)上的荷载lij(i=1,2,…,m;j=1,2,…,p)。
从数学上可以证明,它们分别是相关矩阵m个较大的特征值所对应的特征向量。
4.2.2主成分分析的计算步骤
1、计算相关系数矩阵
rij(i,j=1,2,…,p)为原变量xi与xj的相关系数,rij=rji,其计算公式为
解特征方程 ,常用雅可比法(Jacobi)求出特征值,并使其按大小顺序排列;
分别求出对应于特征值的特征向量,要求=1,即
其中表示向量的第j个分量。
3、计算主成分贡献率及累计贡献率
贡献率:
累计贡献率:
一般取累计贡献率达85%-95%的特征值,所对应的第1、第2、…、第m(m≤p)个主成分。
4、计算主成分载荷
5、各主成分得分
4.2.3主成分分析法在matlab的操作(附录4)。
4.2.4主成分分析法结果集处理
主成分分析法得出影响健康指标8个主要基因和其贡献率如下:
基因序号
2
3
4
5
9
12
15
22
贡献率
0.271463
0.165531
0.129312
0.085919
0.079815
0.061331
0.047829
0.036578
通过各个基因的贡献率计算该8个基因各自对指标影响的权重得:
基因序号
2
3
4
5
9
12
15
22
权重
0.309261
0.18858
0.147317
0.097882
0.090929
0.069871
0.054489
0.041671
最后得出健康指标表达式:
Y=0.3093*X1+0.1886*X2+0.1473*X3+0.0979*X4+0.0909*X5+0.0699*X6+0.0545*X7+0.0417*X8。
通过该表达式计算各人员的Y值得:
1
2
3
4
5
6
7
8
0.0233
0.0369
0.0202
0.0137
0.0189
0.0085
0.0210
0.0232
9
10
11
12
13
14
15
16
0.0104
0.0158
0.0127
0.0120
0.0708
0.0304
0.0135
0.0110
17
18
19
20
21
22
23
24
0.0466
0.0152
0.0236
0.0702
0.0298
0.0289
0.0285
0.0349
25
26
27
28
29
30
31
32
0.0289
0.0372
0.0232
0.0155
0.0074
0.0276
0.0253
0.0168
33
34
35
36
37
38
39
40
0.0161
0.0236
0.0316
0.0231
0.0394
0.0137
0.0342
0.0123
41
42
43
44
45
46
47
48
0.0468
0.0724
0.0108
0.0119
0.0205
0.0148
0.0166
0.0157
49
50
51
52
53
54
55
56
0.0333
0.0231
0.0102
0.0203
0.0209
0.0104
0.0182
0.0105
57
58
59
60
0.0420
0.0134
0.0658
0.0202
通过分析1~20号患者和21~40号健康人员的得分情况将指标分为0~0.02、0.02~0.04、0.04~0.06、0.06~0.08、0.08以上五个段,通过统计各个分段的健康人员和地贫患者的人数得下表:
0~0.02
0.02~0.04
0.04~0.06
0.06~0.08
0.08以上
患者
10
7
1
2
0
正常
6
14
0
0
0
进而得出各分段的患病率:
分段
0~0.02
0.02~0.04
0.04~0.06
0.06~0.08
0.08以上
发病概率
0.63
0.33
1.00
1.00
0.00
通过各分段的发病概率得各人员的患病概率如下:
人员序号
41
42
43
44
45
46
发病概率
1
1
0.63
0.63
0.33
0.63
47
48
49
50
51
52
53
0.63
0.63
0.33
0.33
0.63
0.33
0.33
54
55
56
57
58
59
60
0.63
0.63
0.63
1
0.63
1
0.33
5.模型优缺点
优点:
首先该模型采用主成分分析法可以从大量的指标变量中筛选出少数几综合指标变量,能够最大限度反映出原数据的信息,并且指标之间彼此有一定的相关性,是问题不失真的大大得到简化,该方法能够得出个综合指标,贡献率、累计贡献率、和综合得分,这些数据易于分析。
并且模型能够通过输入某个样品的若干评价健康的主要基因,能够判断该样品患病的概率有多大。
模型整洁,使用起来简单
缺点:
在数据中间处理过程中数据存在舍入误差,稍有影响最终结果,在分段划分区间带有一定的主观性,该模型离通过输入若干评价健康的主要基因百分之百准确判断所求样品是患病。
6.参考文献
1.《概率论与数理统计》、盛骤、高等教育出版社、2008
2.《数学建模与数学实验》、赵静、高等教育出版社、2008
3.《MATLAB基础及其应用教程》、周开利、北京大学出版社、2007
7.附录
附录1:
通过T检验筛选后的基因序号:
基因2
基因3
基因4
基因5
基因9
基因12
基因15
基因22
基因28
基因31
基因33
基因34
基因41
基因42
基因47
基因52
基因54
基因56
基因59
基因62
基因75
基因76
基因78
基因79
基因81
基因84
基因85
基因87
基因89
基因92
基因94
基因96
基因97
基因99
基因101
基因108
附录2:
模型求解主要代码:
2.1建立cwstd.m用用总和标准化法标准化矩阵
%源程序%cwstd.m,用总和标准化法标准化矩阵
functionstd=cwstd(vector)
cwsum=sum(vector,1);%对列求和
[a,b]=size(vector);%矩阵大小,a为行数,b为列数
fori=1:
a
forj=1:
b
std(i,j)=vector(i,j)/cwsum(j);
end
end
3.2建立solve.m文件解决问题
x=xlsread('C:
\Users\dengxinyu\Desktop\桂林理工大学第九届数学建模竞赛暨2013年全国数学建模竞赛选拔赛试题下载\基因链.xls','sheet2');
y=xlsread('C:
\Users\dengxinyu\Desktop\桂林理工大学第九届数学建模竞赛暨2013年全国数学建模竞赛选拔赛试题下载\基因链.xls','sheet3');
z=xlsread('C:
\Users\dengxinyu\Desktop\桂林理工大学第九届数学建模竞赛暨2013年全国数学建模竞赛选拔赛试题下载\基因链.xls','sheet4');
[h,sig,ci]=ttest2(x,y,0.01,0);
x1=[];
y1=[];
z1=[];%x1,y1,z1用于保存t检验存在差异的基因
j=[];%h2用于保存t检验存在差异的基因的序号
h2=[];%h2用于保存t检验存在差异的基因及其序号
fori=1:
110
ifh(i)==1
x1=[x1,x(:
i)];
y1=[y1,y(:
i)];
z1=[z1,z(:
i)];
j=[j,i];
end
end
h2=[j;x1];
v1=cwstd(x1);
v2=cwstd(y1);
v3=cwstd(z1);
vector=v2;
fprintf('相关系数矩阵:
\n')
std=CORRCOEF(vector)%计算相关系数矩阵
fprintf('特征向量(vec)及特征值(val):
\n')
[vec,val]=eig(std)%求特征值(val)及特征向量(vec)
newval=diag(val);
[y,i]=sort(newval);%对特征根进行排序,y为排序结果,i为索引
fprintf('特征根排序:
\n')
forz=1:
length(y)
newy(z)=y(length(y)+1-z);
end
fprintf('%g\n',newy)
rate=y/sum(y);
fprintf('\n贡献率:
\n')
newrate=newy/sum(newy)
sumrate=0;
%newi=[];
fork=length(y):
-1:
1
sumrate=sumrate+rate(k);
newi(length(y)+1-k)=i(k);
ifsumrate>0.85break;
end
end%记下累积贡献率大85%的特征值的序号放入newi中
fprintf('主成分数:
%g\n\n',length(newi));
fprintf('主成分载荷:
\n')
forp=1:
length(newi)
forq=1:
length(y)
result(q,p)=sqrt(newval(newi(p)))*vec(q,newi(p));
end
end%计算载荷
disp(result)
h3=[];%h3用于记录主成分的基因及序号
newrate1=newrate(1:
p);%newate1用于记录主成分的贡献率
fori=1:
p
h3=[h3,h2(:
newi(i))];
end
ratio=[];
fori=1:
p
ratio=[ratio,newrate1(i)/newrate1(p)];%ratio表示系数
end
t=sum(ratio);
fori=1:
p
ratio(i)=ratio(i)/t;
end
score1=[];score2=[];score3=[];
v11=v1(:
1:
p);v21=v2(:
1:
p);v31=v3(:
1:
p);
fori=1:
20
score11=0;score21=0;score31=0;
form=1:
p-6
score11=score11+ratio(m)*v11(i,m);
score21=score21+ratio(m)*v21(i,m);
score31=score31+ratio(m)*v31(i,m);
end
score1=[score1,score11];score2=[score2,score21];score3=[score3,score31];
end
score=[score1;score2;score3];%score的第一行、第二行、第三行用于分别记录健康人员、患者、待审查人员的分数