R语言数据分析与挖掘(第八章):判别分析(3)——费歇尔（Fisher）判别分析
- 2019 年 12 月 13 日
- 筆記
我们之前介绍了判别分析中，因为判别准则的不同，可分为多种判别分析法。常用的有费歇尔（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个样本分类正确。
