R语言数据分析与挖掘(第八章):判别分析(3)——费歇尔（Fisher）判别分析

DoubleHelix

发布于 2019-12-13 10:11:22

8.2K0

发布于 2019-12-13 10:11:22

文章被收录于专栏：生物信息云

概述

我们之前介绍了判别分析中，因为判别准则的不同，可分为多种判别分析法。常用的有费歇尔（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判别的结果。

案例：基于Fisher 判别的iris 数据集分类

下面以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个样本分类正确。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2019-10-24，如有侵权请联系 cloudcommunity@tencent.com 删除

机器学习