R包做主成分分析(含实例)

注:内容非原创,来自对网站http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials/的翻译。

1 什么是主成分分析?

Principal component analysis is used to extract the important information from a multivariate data table and to express this information as a set of few new variables called principal components

主成分分析(Principal component analysis,PCA)属于多元统计分析中的一种,适用于将具有一定相关性的多个变量化为少数几个相互无关的综合指标(主成分),以方便可视化。

比如我这里有两个地区采集到的多个A虫样本(individuals),测量了31个形态指标(即变量variables),我想知道两个地区的A虫是否存在形态上的差异,就用我的数据做一下主成分分析。

真正理解主成分分析原理需要线性代数的知识,这里只作简单解释:

如下图A,数据本来在xy坐标系中显示出来,通过找到数据展现出最大变异的PC1,再把数据投影到PC1上,我们实现了降维,把在二维坐标系中的数据降到了一维。PC2轴是第二重要的方向,它和PC1正交,也就是说二者之间不相关。

The amount of variance retained by each principal component is measured by the so-called eigenvalue.

每个主成分保留的变异又通过特征值来表示。

在数据集中的变量彼此高度相关时,主成分分析尤其有用。变量相关就说明数据存在冗余。

上图来感受下什么是高度冗余(两个变量高度相关):

什么是较低的冗余(两个变量相关不明显):

2 R中如何实现主成分分析

我们通过R包"FactoMineR"和"factoextra"来实现主成分分析和可视化。

下面我会通过自己的一个实例来展示。

#安装和载入所需要的R包

>install.packages(c("FactoMineR","factoextra","ggplot2"))

>library("FactoMineR")

>library("ggplot2")

>library("factoextra")

“FactoMineR”是用来实现数据分析的,“factoextra”是用来实现基于“ggplot2”的可视化的

首先展示一下我的数据:

前两列是样本编号(样本,individuals)和每个样本的采样地点(sites),不参与数据分析。后面31列对应的是31个变量(观测值或者变量,variables)

#从剪切板读取数据,前提是在excel里面已经复制数据集

>collina=read.table("clipboard",header=TRUE)

#显示数据框的前几行,确保正确读取数据

>head(collina)

#进行主成成分分析,第1、2列不参与数据分析(第一列是巢编号,第二列是采样地点)

>collina.pca <- PCA(collina[,-c(1,2)],graph = FALSE)

一般在数据分析之前要进行标准化处理(standardization),尤其是当变量的测量单位不同的时候(如千米,厘米等)。标准化的目的是为了让变量可比,总的来说就是标准差为1,平均值为0。

进行标准化时,数据可进行如下转化:

用R自带的scale()可以进行标准化,但FactoMineR里的PCA()已经包括了标准化的步骤,因此不需要额外的标准化处理

3 对主成分分析的解读和可视化

3.1 主成分的特征值(eigenvalues)和方差(variance)

前面说过,特征值衡量的是主成分保留的变异量,第一主成分保留的变异量最大,依次递减。第一主成分也与数据集方差最大的方向相同。通过特征值,我们可以决定主成分的个数。

#提取主成分(PCs)的特征值和保留方差的百分比(信息),输入:

>eig.val <- get_eigenvalue(collina.pca)

>eig.val

出来的结果如下图:

第一列是每个主成分的特征值。

第二列是每个特征值解释的变异的百分比,计算方法是某主成分特征值/特征值的加和。

第三列是累计的变异的百分比

在本例中,第一、二主成分能解释48.9%的变异。

我们在前面说过特征值可以用来决定保留的主成分的个数,那么怎么来确定呢?

方法1 选择特征值大于1 的主成分保留

在数据已经进行标准化处理的前提下,主成分特征值大于1,意味着主成分比标准化后的原始数据集的任一变量解释的方差要多。根据这个标准,本例中可以保留前7个主成分。

方法二 根据累计方差百分比

一般认为主成分的累计方差百分比达到80%以上,就可以了。根据这个标准,本例中可以保留前7个主成分。

当然80%这个值是可变的,如果你认为70%就够了也是可以的。

方法3 做陡坡图来确定

陡坡图就是将主成分按照第二列其,特征值解释的变异的百分比的大小顺序来排列的柱状图。当趋势降低到已经变得差不多小,没有太大减小的趋势时,认为前面的主成分就是可以保留的。

#做陡坡图,确定主成分的个数

>fviz_eig(collina.pca, addlabels = TRUE,ylim = c(0, 32))

结果如下图:

根据陡坡图,本例中,可以保留前4个主成分。

3.2获得和可视化变量结果

3.2.1 获得主成分分析中变量的结果

#获得主成分分析中变量的结果

>var <- get_pca_var(collina.pca)

>var

这个命令可以帮助我们得到关于变量的主成分分析的结果,帮助我们进行变量的作图,具体帮助如下:

 var$coord:coordinates of variables to create a scatter plot

var$coord提供用变量做散点图的坐标

var$cos2: represents the quality ofrepresentation for variables on the factor map. It’s calculated as the squaredcoordinates: var.cos2 = var.coord * var.coord.

var$cos2代表的是变量在因子图中代表的变异的质量

var$contrib: contains the contributions(in percentage) of the variables to the principal components. The contributionof a variable (var) to a given principal component is (in percentage) :(var.cos2 * 100) / (total cos2 of the component).

var$contrib代表的是变量对于主成分的贡献(百分比表示)。

#查看前4个变量对主成分的贡献

>head(var$contrib, 4)

# 变量对主成分1(PC1)的贡献作图,显示贡献前20的变量

fviz_contrib(collina.pca, choice ="var", axes = 1, top = 20)

# 变量对主成分1(PC1)的贡献作图,显示贡献前20的变量

fviz_contrib(collina.pca, choice ="var", axes = 2, top = 20)

# 变量对主成分1、2(PC1和PC2)的贡献作图,显示贡献前26的变量

fviz_contrib(res.pca, choice ="var", axes = 1:2, top = 26)

其中红色的虚线代表的是期望的平均贡献,如果变量的超过了这条线,认为是对主成分有重要贡献的。

3.2.2 获得主成分分析中样本(individuals)的结果

>ind <- get_pca_ind(res.pca)

>ind

具体含义和变量的结果类似。

要查看各项结果输入:

# Coordinates of individuals

head(ind$coord)

# Quality of individuals

head(ind$cos2)

# Contributions of individuals

head(ind$contrib)

3.2.3 biplot作图

The representation of variables differs from the plot of the observations: The observations are represented by their projections, but the variables are represented by their correlations.

主成分分析的结果主要通过biplot图来呈现,横坐标是主成分1,纵坐标是主成分2,变量通过与主成分的相关关系表现在图中,样本点通过该样本多个变量的观测值在主成分上的投影表现在图中。输入:

>fviz_pca_biplot(collina.pca,

                #individuals

                geom.ind = "point",#只显示样本点,不包括文本信息,即样本编号

                fill.ind = collina$sites,#按sites进行分组

                pointshape = 21, pointsize = "contrib",#点的尺寸大小表示样本的贡献

                palette = "jco",#配色选择“jco”杂志的调色板

                addEllipses = TRUE,#在相同组的样本点为画出密度椭圆

                #variables

                col.var = "cos2",#变量的颜色表示其在因子图中代表的变异的质量,即cos2

                gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),repel=TRUE,

                legend.title = list(fill = "sites"))

出来的PCA分析的biplot图如下:

我们可以看出来,变量的cos2值的大小通过不同的颜色梯度表示出来,每个样本点对主成分的贡献通过点的大小表现出来。两个采样地的样本点也进行了分组。

但是,由于我自身数据的原因,两个采样点的样本分不太开。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,542评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,596评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 158,021评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,682评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,792评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,985评论 1 291
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,107评论 3 410
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,845评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,299评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,612评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,747评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,441评论 4 333
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,072评论 3 317
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,828评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,069评论 1 267
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,545评论 2 362
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,658评论 2 350

推荐阅读更多精彩内容