我们之前介绍了判别分析中,因为判别准则的不同,可分为多种判别分析法。常用的有费歇尔(Fisher)判别分析、贝叶斯(Bayes)判别分析和距离判别分析。在上2篇文章中(判别分析——距离判别法和贝叶斯(Bayes)判别分析)介绍了距离判别分析和贝叶斯判别,本文将介绍贝费歇尔(Fisher)判别分析。
Fisher判别,又称典则判别(canonical discriminant),适用于两类和多类判别。我们将结合两类判别的问题,来介绍一下Fisher判别的原理。
已知有A类和B类两类观察对象,A类有a例,B类有b例,分别记录了X1,X2,……Xm个观察指标,我们称这m个观察指标为判别指标或变量。
Fisher判别法就是找到一个线性组合:
使得综合指标Z在A类的均数与B类的均数的差异尽可能大,而两类的类内综合指标的变异(S2A+S2B)尽可能小,也就是类间差异尽可能大,类内变异尽可能小,即使
达到最大,此时综合指标的公式便称为Fisher判别函数,C1,C2,……,Cm即为判别系数。
建立判别函数后,我们逐例计算出综合指标Zi,求得A类的均数、B类的均数及总均数,按照下式计算判别界值:
如果A类均值大于B类的话,最终的判别规则如下:
?收集22例肝硬化患者的3个指标(腹水量X1,肝长径X2,肝短径X3)中心化、标准化后的资料,其中早期患者A类12例,晚期患者B类10例,如果让我们做Fisher判别:
Step1找到一个类间差异尽可能大,类内变异尽可能小判别函数,各系数通过合并协方差阵代入解方程可得,即Z=-0.070X1+0.225X2-0.318X3;
Step2 逐例计算综合指标Zi,计算出A类、B类的均数和总均数分别为1.428,-1.722,-0.004;
Step3 确定界值,进行两类判别Zc=(1.428+1.722)/2=-0.147,那么-0.147即为界值,Z值高于-0.147即判别为A类,低于则判别为B类。
Step4 判别效果评价,一般要求判别函数的误判概率小于0.1或0.2才有应用价值。本例有4例错判,那么误判概率为4/22=18.2%。
表 22例患者3项指标观察结果
在R语言中,用与进行Fisher判别的最常用函数为lda(),该函数在包MASS中,有2种调用方式。
公式形式:
lda(formula, data, ..., subset, na.action)
Formula:指定用于Fisher判别的公示对象;
Data:指定用于Fisher判别的数据对象,一般为数据框,且优先采用公式中指定的变量数据;
Subset:指定可选向量,表示选择的样本子集;
Na.action:一个函数,指定缺失数据的处理方法,若为NULL,则使用函数
na.omit()删除缺失数据。
矩阵形式:
lda(x, grouping, prior = proportions, tol = 1.0e-4,
method, CV = FALSE, nu, ...)
数据框形式:
lda(x, grouping, ..., subset, na.action)
参数介绍:
X:指定用于Fisher判别的数据对象,可以为矩阵、数据框和包含解释变量的矩阵;
Grouping:因子向量,用于指定样本属于哪一类;
Subset:指定可选向量,表示选择的样本子集;
Na.action:一个函数,指定缺失数据的处理方法,若为NULL,则使用函数na. omit()删除缺失数据;
Prior:指定各个类别的先验概率,默认值为已有训练样本的计算结果;
tol:控制精度,用于判断矩阵是否奇异;
Method:字符串,用于指定估计方法,“mle”表示极大似然估计,“mve”表示使用cov. mve进行估计,“t”表示基于t分布的稳健估计;
CV:逻辑值,若为TRUE,则表示返回值中将包括舍一法的交叉验证结果;
Nu;当参数method设定为“t”时,此处需要设定t分布的自由度。
函数lda()的返回值包括:调用方式、先验慨率、每一类样本的均值和线性判别系数( Fisher判别属于线性判别)。此外,还可以利用函数predict()输出Fisher判别的结果。
下面以iris数据集进行操作演练,首先对数据集中的分类变量进行数据转换,将莺尾花的三个类别分别用1,2,3替代:
> library(MASS)
> data(iris)
> diris<- data.frame(rbind(iris3[,,1], iris3[,,2], iris3[,,3]),species =rep(c(1,2,3), rep(50,3)))
> head(diris)
Sepal.Length Petal.Length Petal.Width Species log.slength species
1 5.1 1.4 0.2 setosa 1.629241 1
2 4.9 1.4 0.2 setosa 1.589235 1
3 4.7 1.3 0.2 setosa 1.547563 1
4 4.6 1.5 0.2 setosa 1.526056 1
5 5.0 1.4 0.2 setosa 1.609438 1
6 5.4 1.7 0.4 setosa 1.686399 1
下面提取训练集和测试集,由于函数lda()中要求训练集的样本与测试集的样本量相等,故此处的训练集和测试集均为75,具体操作如下:
> sa<-sample(1:150, 75)
> sa
[1] 48 46 24 6 32 118 76 131 119 7 64 37 43 70 25 103 27 35
[19] 143 106 72 84 40 79 42 63 129 60 30 93 9 139 85 123 18 58
[37] 57 145 20 95 96 5 133 2 26 75 33 53 147 115 13 89 107 144
[55] 120 50 142 141 122 61 124 126 100 150 17 77 104 132 14 73 82 130
[73] 11 78 110
> table(diris$species[sa])
1 2 3
26 23 26
上面结果显示:测试数据集中莺尾花类别分别为1,2,3的样本量分别为26,23和26,利用函数lda()进行Fisher()判别分析的代码如下:
> z <- lda(species ~ ., diris, prior = c(1,1,1)/3, subset = sa)
> z
Call:
lda(species ~ ., data = diris, prior = c(1, 1, 1)/3, subset = sa)
Prior probabilities of groups:
1 2 3
0.3333333 0.3333333 0.3333333
Group means:
Sepal.L. Sepal.W. Petal.L. Petal.W.
1 4.900 3.352174 1.430435 0.2173913
2 6.024 2.800000 4.344000 1.3560000
3 6.600 2.911111 5.607407 2.0333333
Coefficients of linear discriminants:
LD1 LD2
Sepal.L. 1.357925 -0.2184939
Sepal.W. 1.352581 -1.7520073
Petal.L. -2.580115 1.2345873
Petal.W. -2.928385 -3.1793546
Proportion of trace:
LD1 LD2
0.9958 0.0042
> pre<-predict(z, diris[-sa, ])
> pre$class
[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 2 2 2 2 2 2 2 2 2
[37] 3 2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 2 3 3 3 3
[73] 3 3 3
Levels: 1 2 3
> head(pre$posterior,n=20)
1 2 3
3 1 2.090293e-22 4.088528e-45
7 1 9.587920e-21 1.818849e-42
9 1 1.669810e-17 1.720576e-38
13 1 9.791804e-22 1.457286e-44
15 1 2.446685e-35 1.190135e-62
16 1 1.892658e-31 1.336548e-56
17 1 9.352663e-29 4.043571e-53
19 1 1.966823e-26 2.742911e-50
21 1 3.093880e-23 2.702536e-46
24 1 3.094668e-17 2.374637e-37
25 1 2.669539e-17 5.117169e-38
26 1 2.024177e-19 3.391083e-41
27 1 1.246972e-19 8.594884e-41
28 1 5.185051e-25 1.245723e-48
29 1 1.561341e-25 2.107706e-49
32 1 3.367160e-23 8.900937e-46
34 1 1.709891e-32 1.436627e-58
37 1 3.808029e-29 2.273465e-54
39 1 2.580772e-19 6.413240e-41
40 1 9.124537e-24 5.889644e-47
> head(pre$x,n=20)
LD1 LD2
3 8.061759 0.04341695
7 7.645632 -0.47961183
9 6.990596 0.75802603
13 7.961862 0.81336319
15 10.895552 -1.72199091
16 9.941080 -2.66643913
17 9.373435 -1.97180479
19 8.906347 -1.05038186
21 8.250776 0.03390468
24 6.829625 -0.67915279
25 6.919998 0.41191849
26 7.424586 0.69864642
27 7.379941 -0.63802740
28 8.630472 -0.34451473
29 8.753226 -0.29277273
32 8.181122 -0.84888370
34 10.242668 -1.75992672
37 9.553873 -0.65698037
39 7.383865 0.45936658
40 8.359422 -0.14746461
> class<-pre$class
> diris$species[-sa]
[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 2 2 2 2 2 2 2 2 2
[37] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
[73] 3 3 3
> sum(class==diris$species[-sa])
[1] 72
上面结果中,Call表示调用方法;Prior probabilities of groups表示先验概率;Group means表示每一类样本的均值;Coefficients of linear discriminants表示线性判别系数;Proportion of trace表示比例值。
然后利用class()函数将判别结果展示出来,该函数输出结果为1个列表,其中class、posterior和x三个子列表,分别表示分类结果,后验概率,由于输出量较大,故后面只展示20行记录。
在测试数据集的75条记录中,72个样本分类正确。
本文分享自 MedBioInfoCloud 微信公众号,前往查看
如有侵权,请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ,欢迎热爱写作的你一起参与!