R中的主成分分析

本文概述

• PCA简介
• 一个简单的PCA
• 绘制PCA
• 解释结果
• ggbiplot的图形参数
• 自定义ggbiplot
• 添加新样品
• 将新样品投影到原始PCA上
• 包起来

主成分分析(PCA)是一种用于探索性数据分析的有用技术, 可让你更好地可视化包含多个变量的数据集中的变化。对于”宽”数据集, 其中每个样本都有许多变量, 这特别有用。在本教程中, 你将发现R中的PCA。

更具体地说, 你将解决以下主题：

• 首先, 你将对PCA进行介绍：你将了解主成分以及它们与特征值和特征向量的关系。
• 然后, 你将尝试使用具有简单易懂的数据集的简单PCA。
• 接下来, 你将使用上一部分的结果来绘制你的第一个PCA-可视化非常重要！
• 你还将看到如何开始解释这些可视化的结果, 以及
• 如何使用ggbiplot软件包设置绘图的图形参数！
• 当然, 你希望可视化尽可能自定义, 这就是为什么你还将介绍一些对图进行其他自定义的方法的原因！
• 你还将看到如何将新样本添加到绘图中, 最终将新样本投影到原始PCA上。
• 包起来

PCA简介

正如你在导言中已经读到的那样, 当你使用”宽”数据集时, PCA特别方便。但是为什么呢？

好吧, 在这种情况下, 存在许多变量, 你将无法轻松地以原始格式绘制数据, 从而很难了解其中的趋势。 PCA使你可以查看数据的整体”形状”, 从而确定哪些样本彼此相似, 哪些样本非常不同。这可以使我们识别出相似的样本组, 并确定哪些变量使一组与另一组不同。

它的基础数学有点复杂, 因此我不会赘述过多, 但是PCA的基本知识如下：你获取具有多个变量的数据集, 并通过将原始变量变成较小的数字来简化该数据集”主要组件”。

但是这些到底是什么？主成分是数据中的基础结构。它们是方差最大的方向, 是数据分布最多的方向。这意味着我们试图找到一条最好的直线, 当它沿着数据投影时, 它可以最好地散布数据。这是第一个主要成分, 即显示数据中最大变化的直线。

PCA是对给定数据集的线性变换, 它具有一定数量的变量(坐标)和一定数量的空间的值。这种线性变换使该数据集适合新的坐标系, 从而在第一个坐标上找到最大的方差, 而每个后续坐标都与最后一个坐标正交, 并且方差较小。这样, 你可以将y个样本上的x个相关变量集转换为相同样本上p个不相关的主成分集。

在许多变量相互关联的地方, 它们都将对同一主要成分做出重大贡献。每个主成分加总数据集中总变化的一定百分比。如果你的初始变量之间具有很强的相关性, 则只需几个主要成分, 你就可以估算出数据集中的大多数复杂性。随着添加更多主要成分, 你将汇总越来越多的原始数据集。添加其他组件可以使你对总数据集的估计更加准确, 但也更加麻烦。

特征值和特征向量

就像生活中的许多事物一样, 特征向量和特征值成对出现：每个特征向量都有对应的特征值。简而言之, 特征向量是一个方向, 例如”垂直”或” 45度”, 而特征值是一个数字, 告诉你该方向上数据有多少变化。因此, 特征值最高的特征向量是第一个主成分。

因此, 等等, 在一个数据集中可能发现更多的特征值和特征向量？

没错！退出的特征值和特征向量的数量等于数据集具有的维数。在上面的示例中, 有2个变量, 因此数据集是二维的。这意味着有两个特征向量和特征值。同样, 你会在三维数据集中找到三对。

我们可以根据这些特征向量和特征值重新构建数据集, 而无需更改基础信息。请注意, 根据一组特征值和特征向量对数据集进行重新框架并不需要更改数据本身, 而只是从不同角度查看数据, 这应该可以更好地表示数据。

既然你已经了解了PCA背后的一些理论, 那么你已经准备好将所有这些理论付诸实践！

一个简单的PCA

在本节中, 你将使用简单易懂的数据集尝试PCA。你将使用R中内置的mtcars数据集。该数据集包含32种汽车模型的数据, 这些数据取自《美国汽车》杂志(1974年《汽车趋势》杂志)。对于每辆汽车, 你有11个要素, 以不同的单位(美国单位)表示, 它们如下：

* mpg：油耗(每加仑英里(美国))：功能更强大, 更重的汽车往往消耗更多的燃油。

* cyl：汽缸数：功率更大的汽车通常具有更多的汽缸

* disp：排量(立方英寸)：发动机气缸的总容积

* hp：总马力：这是汽车产生的功率的量度

* drat：后轴比率：这描述了驱动轴的转动与车轮的转动如何对应。较高的值会降低燃油效率。

* wt：重量(1000磅)：非常不言自明！

* qsec：1/4英里时间：汽车的速度和加速度

* vs：发动机缸体：表示车辆的发动机形状是” V”形还是更常见的直形。

* am：变速箱：这表示汽车的变速箱是自动(0)还是手动(1)。

*齿轮：前进挡的数量：跑车往往具有更多的挡位。

* carb：化油器数量：与更强大的发动机相关

请注意, 所使用的单位会有所不同, 并且占据不同的比例。

计算主要成分

由于PCA最适合使用数值数据, 因此你将排除两个类别变量(vs和am)。你剩下9列32行的矩阵, 将其传递给prcomp()函数, 将输出分配给mtcars.pca。你还将设置两个参数(居中和缩放)为TRUE。然后, 你可以使用summary()窥视你的PCA对象。

mtcars.pca <- prcomp(mtcars[, c(1:7, 10, 11)], center = TRUE, scale. = TRUE)

summary(mtcars.pca)

## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6
## Standard deviation     2.3782 1.4429 0.71008 0.51481 0.42797 0.35184
## Proportion of Variance 0.6284 0.2313 0.05602 0.02945 0.02035 0.01375
## Cumulative Proportion  0.6284 0.8598 0.91581 0.94525 0.96560 0.97936
##                            PC7    PC8     PC9
## Standard deviation     0.32413 0.2419 0.14896
## Proportion of Variance 0.01167 0.0065 0.00247
## Cumulative Proportion  0.99103 0.9975 1.00000

你获得9个主要组件, 称为PC1-9。每一个都解释了数据集中总变化的百分比。这就是说：PC1解释了总方差的63％, 这意味着数据集中将近三分之二的信息(9个变量)可以仅由一个主成分封装。 PC2解释了23％的方差。因此, 通过了解样本相对于PC1和PC2的位置, 你可以非常准确地了解样本相对于其他样本的位置, 因为仅PC1和PC2可以解释86％的方差。

让我们调用str()来查看你的PCA对象。

str(mtcars.pca)

## List of 5
##  $ sdev    : num [1:9] 2.378 1.443 0.71 0.515 0.428 ...
##  $ rotation: num [1:9, 1:9] -0.393 0.403 0.397 0.367 -0.312 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:9] "mpg" "cyl" "disp" "hp" ...
##   .. ..$ : chr [1:9] "PC1" "PC2" "PC3" "PC4" ...
##  $ center  : Named num [1:9] 20.09 6.19 230.72 146.69 3.6 ...
##   ..- attr(*, "names")= chr [1:9] "mpg" "cyl" "disp" "hp" ...
##  $ scale   : Named num [1:9] 6.027 1.786 123.939 68.563 0.535 ...
##   ..- attr(*, "names")= chr [1:9] "mpg" "cyl" "disp" "hp" ...
##  $ x       : num [1:32, 1:9] -0.664 -0.637 -2.3 -0.215 1.587 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
##   .. ..$ : chr [1:9] "PC1" "PC2" "PC3" "PC4" ...
##  - attr(*, "class")= chr "prcomp"

我不会在这里详细描述结果, 但是你的PCA对象包含以下信息：

• 每个主要成分的中心点($ center), 缩放比例($ scale), 标准差(sdev)
• 初始变量和主成分($ rotation)之间的关系(相关性或反相关性等)
• 每个样本的值, 以主成分($ x)表示

绘制PCA

现在是时候绘制你的PCA了。你将制作一个双向图, 其中包括每个样本在PC1和PC2方面的位置, 还将向你展示初始变量如何映射到该位置。你将使用ggbiplot程序包, 该程序包提供了一种用户友好且漂亮的功能来绘制双线图。双线图是一种绘图类型, 可让你直观地看到样本在我们的PCA中如何相互关联(哪些样本相似, 哪些样本不同), 并同时揭示每个变量对每个主要成分的贡献。

在开始之前, 请不要忘记先安装ggbiplot！

library(devtools)
install_github("vqv/ggbiplot")

接下来, 你可以在PCA上调用ggbiplot：

library(ggbiplot)

ggbiplot(mtcars.pca)

轴被视为从中心点开始的箭头。在这里, 你会看到变量hp, cyl和disp都对PC1有所贡献, 这些变量中的值越高, 则将样本移到该图的右侧。这使你可以看到数据点与轴的关系, 但是如果不知道哪个点对应于哪个样本(汽车), 就不是很有帮助。

你将为ggbiplot提供一个参数：让我们给它mtcars的行名作为标签。这将使用相关汽车的名称来命名每个点：

ggbiplot(mtcars.pca, labels=rownames(mtcars))

现在你可以看到哪些汽车彼此相似。例如, 玛莎拉蒂Bora, 法拉利Dino和福特Pantera L都聚集在顶部。这是有道理的, 因为所有这些都是跑车。

你还能如何尝试更好地了解你的数据？

解释结果

也许如果你看看每辆车的来历。你将它们分为三类(类别？)之一, 每类分别用于美国, 日本和欧洲的汽车。你为此信息创建一个列表, 然后将其传递给ggbiplot的groups参数。你还将椭圆参数设置为TRUE, 这将在每个组周围绘制一个椭圆。

mtcars.country <- c(rep("Japan", 3), rep("US", 4), rep("Europe", 7), rep("US", 3), "Europe", rep("Japan", 3), rep("US", 4), rep("Europe", 3), "US", rep("Europe", 3))

ggbiplot(mtcars.pca, ellipse=TRUE, labels=rownames(mtcars), groups=mtcars.country)

现在你会看到一些有趣的东西：美国汽车在右边形成了一个独特的集群。纵观轴线, 你会发现美式汽车的特点是缸, 分配和重量的值很高。另一方面, 日系汽车的特点是平均每加仑(mpg)高。欧洲汽车处于中间位置, 并且不如任何一组更紧密地聚集。

当然, 你有许多可用的主要组件, 每个组件与原始变量的映射方式都不相同。你还可以通过使用choices参数要求ggbiplot绘制这些其他组件。

让我们看一下PC3和PC4：

ggbiplot(mtcars.pca, ellipse=TRUE, choices=c(3, 4), labels=rownames(mtcars), groups=mtcars.country)

你在这里看不到太多, 但这并不奇怪。 PC3和PC4解释了总变化的很小的百分比, 因此如果你发现它们非常有用并且将各组分开或显示出明显的模式, 这将是令人惊讶的。

让我们花点时间回顾一下：使用mtcars数据集执行PCA后, 我们可以看到美国汽车和日本汽车之间的明显分隔, 以及与cyl, disp, wt和mpg密切相关的主要成分。这为我们将来的分析提供了一些线索。如果我们试图建立一个分类模型来识别汽车的起源, 那么这些变量可能会很有用。

ggbiplot的图形参数

你还可以使用其他一些变量来更改双拼。你可以在数据集的中心添加一个圆圈(circle参数)：

ggbiplot(mtcars.pca, ellipse=TRUE, circle=TRUE, labels=rownames(mtcars), groups=mtcars.country)

你还可以缩放样本(obs.scale)和变量(var.scale)：

ggbiplot(mtcars.pca, ellipse=TRUE, obs.scale = 1, var.scale = 1, labels=rownames(mtcars), groups=mtcars.country)

你也可以使用var.axes完全删除箭头。

ggbiplot(mtcars.pca, ellipse=TRUE, obs.scale = 1, var.scale = 1, var.axes=FALSE, labels=rownames(mtcars), groups=mtcars.country)

自定义ggbiplot

由于ggbiplot基于ggplot函数, 因此你可以像使用任何ggplot一样使用相同的图形参数集来更改双曲线。在这里, 你将要：

• 使用scale_colour_manual()指定要用于组的颜色
• 用ggtitle()添加标题
• 指定minimal()主题
• 用theme()移动图例

ggbiplot(mtcars.pca, ellipse=TRUE, obs.scale = 1, var.scale = 1, labels=rownames(mtcars), groups=mtcars.country) +
  scale_colour_manual(name="Origin", values= c("forest green", "red3", "dark blue"))+
  ggtitle("PCA of mtcars dataset")+
  theme_minimal()+
  theme(legend.position = "bottom")

添加新样品

好的, 假设你要向数据集中添加一个新样本。这是一款非常特别的汽车, 其数据与众不同。它超级强大, 具有60缸发动机, 惊人的燃油经济性, 无齿轮且非常轻便。是木星的”太空船”。

你可以将其添加到现有数据集中, 并查看其相对于其他汽车的位置吗？

你将其添加到mtcars, 创建mtcarsplus, 然后重复分析。你可能希望能够看到最喜欢哪个地区的汽车。

spacecar <- c(1000, 60, 50, 500, 0, 0.5, 2.5, 0, 1, 0, 0)

mtcarsplus <- rbind(mtcars, spacecar)
mtcars.countryplus <- c(mtcars.country, "Jupiter")

mtcarsplus.pca <- prcomp(mtcarsplus[, c(1:7, 10, 11)], center = TRUE, scale. = TRUE)

ggbiplot(mtcarsplus.pca, obs.scale = 1, var.scale = 1, ellipse = TRUE, circle = FALSE, var.axes=TRUE, labels=c(rownames(mtcars), "spacecar"), groups=mtcars.countryplus)+
  scale_colour_manual(name="Origin", values= c("forest green", "red3", "violet", "dark blue"))+
  ggtitle("PCA of mtcars dataset, with extra sample added")+
  theme_minimal()+
  theme(legend.position = "bottom")

但这是一个天真的假设！随着样品的加入, PCA的形状发生了巨大变化。当你更仔细地考虑此结果时, 它实际上是很有意义的。在原始数据集中, 某些变量(例如cyl和mpg)之间具有很强的相关性, 这些变量对PC1有所贡献, 从而沿该轴将各个组彼此分开。但是, 当你使用额外的样本执行PCA时, 将不存在相同的相关性, 从而扭曲了整个数据集。在这种情况下, 效果特别强, 因为你的额外样本在多个方面都是极端异常。

如果要查看新样本与初始PCA产生的组的比较, 则需要将其投影到该PCA上。

将新样品投影到原始PCA上

这意味着定义的主要成分与你的航天器样本无关, 然后你可以通过应用PCA产生的转换来计算相对于其他样本的航天器放置位置。你可以认为这是, 你无需获取所有样本的均值并使航天器偏斜此均值, 而是获取其余样本的均值并查看与此相关的航天器。

这意味着你只需相对于PCA中心(mtcars.pca $ center)缩放航天器的值。然后, 将PCA矩阵的旋转应用于航天器样本。然后, 可以将航天器的投影值rbind()到pca $ x矩阵的其余部分, 然后像以前一样将其传递给ggbiplot：

s.sc <- scale(t(spacecar[c(1:7, 10, 11)]), center= mtcars.pca$center)
s.pred <- s.sc %*% mtcars.pca$rotation


mtcars.plusproj.pca <- mtcars.pca
mtcars.plusproj.pca$x <- rbind(mtcars.plusproj.pca$x, s.pred)


ggbiplot(mtcars.plusproj.pca, obs.scale = 1, var.scale = 1, ellipse = TRUE, circle = FALSE, var.axes=TRUE, labels=c(rownames(mtcars), "spacecar"), groups=mtcars.countryplus)+
  scale_colour_manual(name="Origin", values= c("forest green", "red3", "violet", "dark blue"))+
  ggtitle("PCA of mtcars dataset, with extra sample projected")+
  theme_minimal()+
  theme(legend.position = "bottom")

这个结果完全不同。请注意, 所有其他样本都返回其初始位置, 而航天器则放置在中间附近。你的额外样本不再会使整体分布发生偏差, 但是不能将其分配给特定的组。

但是PCA的预测或重新计算哪个更好？

它在某种程度上取决于你要回答的问题。重新计算表明航天器是一个离群值, 投影表明你不能将其放置在现有组中。当通过PCA进行探索性数据分析时, 执行这两种方法通常很有用。在深入研究数据集之前, 这种探索性分析通常是一个很好的起点。你的PCA会告诉你哪些变量将美国汽车与其他汽车区分开来, 而太空车在我们的数据集中是一个离群值。下一步可能是查看这些关系是否适用于其他汽车, 或者查看汽车如何按品牌或类型(跑车, 四轮驱动汽车等)进行分组。

包起来

所以你有它！

你已经了解了PCA的原理, 如何创建双线图, 如何微调该图, 并了解了两种不同的方法将样品添加到PCA分析中。谢谢阅读！

如果你想了解有关R的更多信息, 请参加srcmini的免费R入门课程。