PCA作图⾥⾯的箭头是⼲嘛⽤的?
作图的⽬的是希望在图⾥⾯发现问题或者解释问题,当然更本质⼀点就是你想解决什么问题?
前⼏天做了⼀个PCA的图,图是画出来了,但是问题有很多,⽐如说主成分是是啥意思,图⾥⾯的箭头有什么含义?为了不做⽆意义的重复,所以写⼀篇⽂章尝试做⼀个解释。
数据线 英文我们以R语⾔⾃带的数据集iris作为例⼦来演⽰。
data(iris)
iris翻译成中⽂就是鸢(yuan, 第⼀声)尾花(如下图), 我建议你在R语⾔⾥⽤?iris了解更多这个数据集的出处。
iris
先让我们⽤head简单看下这个数据集的前⾯⼏⾏。其中前⾯四列是⼀朵鸢尾花的⼀些形态特征衡量值,⾃⾏翻译各个单词的含义。最后⼀列是这朵花的所属物种,根据分类,这些花来⾃于"tosa, versicolor,virginica"
> head(iris)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.
2 tosa
2 4.9 3.0 1.4 0.2 tosa
3 4.7 3.2 1.3 0.2 tosahighlands
4 4.6 3.1 1.
5 0.2 tosa
5 5.0 3.
6 1.4 0.2 tosa
6 5.4 3.9 1.
7 0.4 tosa
在平常⼈眼⾥, 看到这些数值可能没有任何想法,也不知道这到底是什么品种,但是对于⼀个研究鸢尾花的⼈,这些数值可能会在他们的脑海⾥重构出⼀朵鸢尾花,然后迅速判断出这是什么品种。我和阅读⽂章的⼈都⼀样,也不知道如何区分鸢尾花的品种,如果想去研究这种规律,应该怎么做呢?
⼈类的学习过程过程⼤多数多看多提出⼀些模型规律,然后基于这个规律去判断,如果错了改正模型,直到找到⼀个好⽤的标准为⽌。因此,后⾯我们也就是在当前的数据集下多试试,看看能不能提出⼀些模型。
可视化是⾮常强⼤的⼯具,让我们先看看能不能⽤“ Sepal.Length”和“ Sepal.Width”进⾏区分
钢铁侠英文
library(ggplot2)
ggplot(iris,aes(x=Sepal.Length, y=Sepal.Width)) + geom_point(aes(colour=Species))
supports
fig1
从上图中,我们发现根据萼⽚(pal)的宽度和长度似乎能够区分出"tosa"和"versicolor","virginica",但是"versicolor","virginica"不太好分。
让我们再试试花瓣(petal)的长度和宽度。从下图,我们可以发现在这两个属性下,不同品种的鸢尾花似乎有了明显的分界线。
ggplot(iris,aes(x=Petal.Length, y=Petal.Width)) + geom_point(aes(colour=Species))
image
我们可以建⽴⼀个模型就是
当花瓣长度⼩于2时就⼀定是tosa
当花瓣的宽度在0.75~1.75之间,花瓣的长度在3~5,则是versicolor,
当花瓣的宽度⼤于1.75时,花瓣的长度⼤于5,则是virginica
这个模型依旧还不完善,存在⼀些点难以区分。不过我们还可以以花瓣长度和萼⽚长度,花瓣宽度和萼⽚宽度等你能想到的组合进⾏作图,说不定能找到更好的模型。更好的情况是,如果⼈类能够想象出⼀个四维的空间,在这个空间⾥⾯同时展现这四个变量,或许就能更加直观的找到规律,可惜,我们⽬前还没有这种能⼒
集合名词
有没有⼀种⽅法能够把这种多维数据降维,也就是抓住数据集的主要⽭盾呢?当然是有,这个⽅法就是100多年前,由卡尔·⽪尔森提出的主成分分析,通过将原来众多具有⼀定相关性的变量,重新组合成⼀组新的互相⽆关的综合变量。⽐如说我们可以将花瓣长度和宽度组合成花瓣的⾯积,当然实际计算可能更加复杂,你需要看。
反正PCA最终的⽬的就是让原本的数据在降维后⽅差最⼤,也就是能把数据分开。我们⼀般⼈记住这⼀点,然后会⽤软件,能解读结果就好。
R语⾔能做主成分分析功能很多,⽐如说prcomp,princomp,principal. 这⾥⽤到是prcomp,它使⽤的是奇异值分解⽽不是特征值分解,⾄于这两个有啥区别,其实这不是我关⼼的,当然你想深⼊了解到话,建议阅读。
代码如下:
halogen
prcomp会返回⼀个列表,⾥⾯包含如下内容,你可以⽤sults)查看
sdev: 对应每个主成分的标准偏差,值越⼤说明解释变异越⼤。
沪江考研
butterflykissrotation: 变量的负荷(loading)矩阵,⽤于衡量⼀个变量在各个主成分重要性。
x: 原矩阵经线性变化后的矩阵(下图不显⽰,可以⽤sults$x提取)
PCA results
虽然有现成的⼯具提取sults⾥的数据进⾏作图,但这样⽆助于理解,我们需要⾃⼰⼿动提取结果进⾏作图。
iris_rotate_df <- as.data.sults$x[,c("PC1","PC2")])
iris_rotate_df$Species <- iris$Species
p <- ggplot(iris_rotate_df, aes(x=PC1, y=PC2)) + geom_point(aes(colour=Species))
print(p)
PCA plot
这个结果虽然也⽆法完美地区分出"versicolor"和"virginica",但是⽐我们盲⽬的寻找好多了。这⾥的主成分是原来的各个变量线性组合结果。
当然我们还有⼀个问题,各个变量在不同的主成分⾥的权重是多少?这个就要看负荷矩阵了,也就是Rotation部分。
在PCA中,我们将协⽅差矩阵分成了标量部分(特征值)和有向部分(特征向量),通过计算特征向量和开⽅后特征值的乘积,就称之为负荷(loading)
我们发现在PC1⾥⾯Petal.Length的绝对值最⼤,⽽在PC2⾥⾯,Sepal.Width的绝对值最⼤,于是我们就可以从原来的数据集中提取这两个变量进⾏作图了。
ggplot(iris,aes(x=Sepal.Width, y=Petal.Length)) + geom_point(aes(colour=Species))
下⾯的作图直接画图的结果,是不是感觉和上⾯的PCA图很像。当我⼿动将其旋转后得到下⾯的右图,你会发现这和PCA的结果⼏乎⼀摸⼀样。
loading
那么我们如何在图上展⽰各个变量在各个PC的占⽐呢?这个就需要画⼏个箭头了,
但是p + annotate("gment", x=0,xend=0.35,y=0,yend=-0.65,arrow=arrow()) +
annotate("text", x=0.35,y=-0.68, label="Sepal.Length") +
annotate("gment", x=0,xend=-0.08,y=0,yend=-0.73,arrow=arrow()) +
annotate("text", x=-0.08,y=-0.75, label="Sepal.Width")
下图中"Septal.Width"朝下,Petal.Length朝右,如果你⼿动旋转⼀下,就是差不多是以"Septal.With"为X轴,"Petal.Length"为Y轴的结果
深圳英孚教育
了。 所以箭头的朝向不重要,重点是长度。
arrow
当然也有现成的包来完成上图, 不过知道原理后的你应该会解读结果了吧
library(ggord)
sults, iris$Species, size=1.25)