R语言数据分析与挖掘(第八章):判别分析(1)——距离判别法

  • 2019 年 12 月 13 日
  • 筆記

判别分析概述

判别分析是判断个体所属类别的一种多元统计分析方法。它在医学领域有着广泛的应用,主要有疾病诊断、疾病预测和病因学分析。例如,根据病人的症状、生化指标判断病人得的是什么疾病,根据病人症状的严重程度或者指标的高低预测病人的预后等等。比如,高血压、高血糖、动脉硬化程度这些都是脑血管疾病的患病危险因素;那么如果知道了人体的这些指标,并对这些数据进行分析,就可以对尚未明确诊断的人是否发生脑血管疾病进行预测;对于很可能是脑血管疾病的人就可以事先给予预防,或者在入院后尽快得到救治,提高诊疗有效率。

判别分析也属于对事物现象进行分类的统计分析方法,它和聚类分析不同的地方在于:聚类分析(后面会讲)事先并不知道分型情况,而判别分析需要事先知道分型情况,已知的分型数据又叫训练数据。判别分析需要事先得到一些已经明确知道诊断结果的训练数据,利用这些数据建立判别准则,然后依据准则对未知类别的预测值进行判别。如果是对于分类不明的数据,可以先用聚类分析对这组数据进行分类,然后再用判别分析对新建立的类别进行判断。

 在判别分析中,因为判别准则的不同,可分为多种判别分析法。常用的有费歇尔(Fisher)判别分析、贝叶斯(Bayes)判别分析和距离判别分析。我们这里先介绍距离判别法。

距离判别的基本思想是样品X离哪个总体的距离最近,就判断X属于哪个总体。在判别法中根据不同的功能需求,会经常用到dist()、mahalanobis()和wmd()这三个函数。

函数介绍

1

dist()函数

进行距离判别的第一步是计算"距离",最常用的函数为dist(),该函数按照指定方法计算数据矩阵之间的距离,书写格式为:

dist(x, method = "euclidean", diag = FALSE, upper = FALSE, p = 2)

参数介绍:

X:指定用于计算距离的数据对象,可以为一个矩阵、数据框或者“dist" 类的对象;

Method:指定测度距离的方法,"euclidean" 表示欧氏距离,"maximum" 表示切比雪夫距离,"manhattan"表示绝对值距离,"canberra"表示 Lance距离,"binary"表示定性变量距离,"minkowski"表示明科夫斯基距离,使用时要指定p值,默认值为"euclidean";

Diag:逻辑值,为TRUE时表示输出距离矩阵的对角线,默认值为FALSE;

Upper:逻辑值,为TRUE时表示输出距离矩阵的上三角部分,默认值为FALSE;

p:指定明科夫斯基距离的权数,当参数method设置为"minkowski"时, 此参数需要进行设定,默认值为2。

2

mahalanobis()函数

不难发现,函数dist()不能用于计算马氏距离,下面介绍一个专门用于计算马氏距离的函数: mahalanobis(), 其基本书写格式为:

mahalanobis(x, center, cov, inverted = FALSE, ...)

参数介绍:

x:指定用于计算距离的数据对象,p维的数据向量或矩阵;

Center:指定分布的均值,即总体均值;

Cov:指定分布的协力差,即总体协方差,一般用样本的协方差进行评估;

inverted:逻辑值,若为TRUE,则表示参数cov应该包括协方差的逆。

3

wmd()函数

上述介绍的两个函数均返回距离值,而不能直接判别,下面介绍一个可直接用于判别的函数: wmd(), 该函数存在于WMDR包中,可用于实现加权马氏距离的判别,它利用函数mahalanobis()计算出马氏距离,然后进行判别分析,最终返回包含结果和准确度的表单,其基本书写格式为:

wmd(TrnX,TrnG,Tweight = NUL, TstX = NULL, var.equal = F)

参数介绍:

TmX:指定训练集的数据对象,可以为矩阵或数据框;

TrnG:一个因子类的向量,用于指定已知的训练样本的分类;

Tweight:指定权重,若没有进行指定,则软件默认使用主成分分析中的相应贡献率作为权重,默认值为NUll,表示不进行加权,采用传统的马氏距离判别法;

TstX: 指定测试集的数据对象,可以为向量、矩阵或数据框,若为向量,则将被识别为单个案例的行向量,默认值为NULL,表示直接对训练集进行判别;

Var.equal: 指定不同类别之间是否具有相同的协方差,默认值为F。

需要注意的是,函数wmd()中训练集的样本量与测试集的样本量相等,否则R语言会报错。

综合案例:基于距离判别的iris数据集分类

下面利用iris 数据集进行操作演练,由于该数据集中鸢尾花的种类有三种,下面将原始数据分为训练集和测试集,分别包含随机抽取的100个样本,具体操作如下:

> set.seed(1234)  > sa<-sample(1:150,100)  > sa    [1]  18  93  91  92 126 149   2  34  95  73  98  76  40 127 138 114  39   [18]  36  25  31  42 136  21   6  28 102  66 113 125 137  55  32 133  60   [35]  22  88  23  30 112  90  61  71 143  67  35  53 109  50 132  78   8   [52] 131 104  49  15  48  47  70  17 101 148   4 146 144 128 110  26  43   [69]   5  46  10 140  87  85   7 134  29 121  24 142  65  33  80  37  13   [86]  59 122  20  68 120  62  54 100  58 141  74 147 111 145  38  > dtrain<-iris[sa,1:5]  > head(dtrain)      Sepal.Length Sepal.Width Petal.Length Petal.Width    Species  18           5.1         3.5          1.4         0.3     setosa  93           5.8         2.6          4.0         1.2 versicolor  91           5.5         2.6          4.4         1.2 versicolor  92           6.1         3.0          4.6         1.4 versicolor  126          7.2         3.2          6.0         1.8  virginica  149          6.2         3.4          5.4         2.3  virginica  > dtest<-iris[-sa,1:4]  > head(dtest)     Sepal.Length Sepal.Width Petal.Length Petal.Width  1           5.1         3.5          1.4         0.2  3           4.7         3.2          1.3         0.2  9           4.4         2.9          1.4         0.2  11          5.4         3.7          1.5         0.2  12          4.8         3.4          1.6         0.2  14          4.3         3.0          1.1         0.1

上面代码表示:首先利用函数set.seed()释放随机种子,然后利用函数sample()从1到150中随机抽取100个数据,将抽取的数据对应到数据中观测样本的编号,得到训练数据集,剩下的数据集即为测试集。下面对训练集中的观测样本安装变量Species分类:

> d1<-subset(dtrain,Species=="setosa");dim(d1)  [1] 38  5  > d2<-subset(dtrain,Species=="versicolor");dim(d2)  [1] 29  5  > d3<-subset(dtrain,Species=="virginica");dim(d3)  [1] 33  5

上述输出结果显示:在训练集中种类为setosa的鸢尾花共有38条记录,种类为versicolor()和virginica的鸢尾花分别有29和33条记录。

下面利用函数mahalanobis()计算马氏距离:

> ma1<-mahalanobis(dtest,colMeans(d1[,1:4]),cov(d1[,1:4]))  > ma2<-mahalanobis(dtest,colMeans(d2[,1:4]),cov(d2[,1:4]))  > ma3<-mahalanobis(dtest,colMeans(d3[,1:4]),cov(d3[,1:4]))  > (distance<-cbind(ma1,ma2,ma3,iris[-sa,5]))               ma1         ma2        ma3  1      0.3590061 135.3154437 271.444570 1  3      1.5544105 112.3495585 242.943354 1  9      3.4278559  86.6318539 203.172929 1  11     1.8993357 151.9381072 290.189235 1  12     1.9068175 117.1927197 235.675433 1  14     9.0397389 110.4697730 241.640736 1  16     9.8088108 211.1672458 359.672511 1  19     6.3669561 147.5440882 282.272697 1  27     4.0876272 104.1682597 226.561182 1  41     1.8477792 132.0256229 270.980300 1  44    17.8963186 102.9131186 222.375954 1  45    12.0006760 125.8813043 238.336324 1  51   519.2519291   5.8057277  26.778995 2  52   483.9369049   2.7504623  21.480586 2  56   440.2499733   3.2297851  14.484601 2  57   550.1907724   3.7024849  16.423962 2  63   314.6573331   6.6213187  26.256011 2  64   505.8581500   1.9220704  11.529449 2  69   509.6105640  14.5675407  11.236568 2  72   347.7106464   2.1919434  28.310069 2  75   406.4106931   1.7786523  25.078831 2  77   548.0360755   3.6213575  14.421208 2  79   485.3833731   0.7923886  12.021558 2  81   278.4175024   1.7677099  25.346488 2  82   246.8996580   2.8295498  30.772019 2  83   308.7761180   1.1340186  28.096945 2  84   653.7060829   8.5350447   3.495719 2  86   506.4900210   6.8295042  22.639212 2  89   358.5380439   3.0143070  26.291627 2  94   189.8000871   4.8847665  37.271789 2  96   357.0214367   5.0103361  27.493277 2  97   378.1639201   1.3599246  21.536322 2  99   169.1113735  11.7826462  52.221092 2  103 1041.4600884  22.9217779   1.827084 3  105 1041.9268906  28.4050506   2.251143 3  106 1286.1243510  36.4572358  10.796066 3  107  560.0828547  18.2689904  10.068365 3  108 1082.7402091  29.1658678   8.484063 3  115  951.7079875  57.3638916   6.100539 3  116  943.8604134  31.8419254   1.768798 3  117  814.6300092  10.5481905   1.706876 3  118 1336.8526085  31.3734151  11.353519 3  119 1489.6049927  69.4397271  29.160649 3  123 1302.1372419  44.5459977  16.435043 3  124  666.3834809  11.6169643   2.831568 3  129  948.9251799  26.6719664   1.741693 3  130  864.8827533  15.9814974   6.508066 3  135  747.0655076  29.1338728  11.601745 3  139  635.8783895   7.2968399   4.865241 3  150  708.1090088   9.5614468   3.787757 3

上述代码表示:分别对训练集计算三种类别的马氏距离,其中函数colMeans()表示按列计算均值;训练集中每一个观测样本分别对应三个马氏距离,然后利用函数cbind()将三个马氏距离值与原始数据集中测试样本对应的分类合并在一起,输出结果如上所示。对于测试集中的每一个观测样本而言,三个马氏距离中最小的那一个所对应的类别即为测试样本属于的类别,如第一条记录中,第一个马氏距离的值明显小于另外两个,故第一条记录应归为第一类,即该鸢尾花属于setosa类,与原始数据集中的分类一致,说明分类正确。

上述分类过程没有直接得到结果,需要比较每一个测试样本对应的三个马氏距离进行分类,如需了解该分类的效率,还要进行后续操作,分类过程较为烦琐,下面利用函数wmd()直接得到分类结果:

install.packages("WMDR")# 对高版本的R已经不适用  library(WMDR)  dta<-iris[,1:4]  species<-gl(3,50)  wmd(dta,species)  wmd(dta,species,diag(rep(0.25,4)))#设定权重,增加准确率。

函数wmd()的输出结果中,第一部分表示对150个观测值进行分类的结果(由于函数中没有指定测试集和训练集,故软件默认训练集和测试集均为同一个); "num of wrong judgement"表示判别错误的样本编号,可以看到有8个样品错判; "samples divided to"和"samples actually belongs to" 分别表示上述样品的错误判别归类和真实类别;"percent of right judgement" 表示判别结果的正确率,由(150-8) /150得到正确率为94.7%。