判别分析概述
判别分析是判断个体所属类别的一种多元统计分析方法。它在医学领域有着广泛的应用,主要有疾病诊断、疾病预测和病因学分析。例如,根据病人的症状、生化指标判断病人得的是什么疾病,根据病人症状的严重程度或者指标的高低预测病人的预后等等。比如,高血压、高血糖、动脉硬化程度这些都是脑血管疾病的患病危险因素;那么如果知道了人体的这些指标,并对这些数据进行分析,就可以对尚未明确诊断的人是否发生脑血管疾病进行预测;对于很可能是脑血管疾病的人就可以事先给予预防,或者在入院后尽快得到救治,提高诊疗有效率。
判别分析也属于对事物现象进行分类的统计分析方法,它和聚类分析不同的地方在于:聚类分析(后面会讲)事先并不知道分型情况,而判别分析需要事先知道分型情况,已知的分型数据又叫训练数据。判别分析需要事先得到一些已经明确知道诊断结果的训练数据,利用这些数据建立判别准则,然后依据准则对未知类别的预测值进行判别。如果是对于分类不明的数据,可以先用聚类分析对这组数据进行分类,然后再用判别分析对新建立的类别进行判断。
在判别分析中,因为判别准则的不同,可分为多种判别分析法。常用的有费歇尔(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 数据集进行操作演练,由于该数据集中鸢尾花的种类有三种,下面将原始数据分为训练集和测试集,分别包含随机抽取的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%。
本文分享自 MedBioInfoCloud 微信公众号,前往查看
如有侵权,请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ,欢迎热爱写作的你一起参与!