一文了解R语言数据分析 —-主成分分析

专题介绍：R是一种广泛用于数据分析和统计计算的强大语言，于上世纪90年代开始发展起来。得益于全世界众多爱好者的无尽努力，大家继而开发出了一种基于R但优于R基本文本编辑器的R Studio（用户的界面体验更好）。也正是由于全世界越来越多的数据科学社区和用户对R包的慷慨贡献，让R语言在全球范围内越来越流行。其中一些R包，例如MASS，SparkR， ggplot2，使数据操作,可视化和计算功能越来越强大。R是用于统计分析、绘图的语言和操作环境。R是属于GNU系统的一个自由、免费、源代码开放的软件，它是一个用于统计计算和统计制图的优秀工具。R作为一种统计分析软件，是集统计分析与图形显示于一体的。它可以运行于UNIX、Windows和Macintosh的操作系统上，而且嵌入了一个非常方便实用的帮助系统，相比于其他统计分析软件，R的学术性开发比较早，适合生物学和医学等学术学科的科研人员使用。

(文末有赠书福利）

在数据挖掘项目中，经常会遇到的情况是有很多特征可以用，这是一件好事，但是有的时候数据中存在很多冗余情况，也就是说数据存在相关性或者共线性。在这种情况下对于分析带来了很多麻烦。不必要的特征太多会造成模型的过于复杂，共线性相关性会造成模型的不稳定，即数据微小的变化会造成模型结果很大的变化。主成分分析是解决这种问题的一个工具。

一、概述

主成分分析简称PCA，PCA是一个很好的预处理工具，用于挑选最相关的变量线性组合，并在我们的预测模型中使用它们。PCA帮助我们找出解释数据中变化最大的变量。PCA也是解决数据间存在共线性问题的一种方式。

PCA的唯一缺点是它以无监督的方式生成主成分，也就是与数据的标签没有任何关系。如果最终的目的是构建模型的话，选出来的主成分不一定对于标签或者说因变量有很好的预测能力。在本章节会介绍PCA的一些基本概念，并且使用R语言进行实现。

1. 维度相关的问题

有两个因素会使得数据的维度变得非常高。首先，具有不同级别的分类变量，其次是数据集中的冗余变量。当我们在建模中包含分类变量时，我们使用虚拟变量，并且分类变量中的级别数越多，我们创建的虚拟变量就越多，这会增加维度。冗余变量代表数据集中的信息不够精简，多个特征描述了同一信息，例如，数据集中有三个特征，收益，成本，利润率，这三个特征就存在冗余，因为第三个特征是通过前两个特征生成的，因此可以删除掉冗余特征。

多重共线性（共线性）可以定义为回归模型中的两个或更多个变量高度相关的现象，这意味着，模型数据中存在一定的冗余。

变量之间存在冗余变量和多重共线性可以以多种方式影响变量的分析，其中一些在下面例子：

• 破坏参数估计的稳定性。

• 混淆模型解释。

• 增加过度拟合的风险。

• 增加计算时间。

输入数据的微小变化（如删除或添加单个变量）会导致模型发生较大变化，甚至导致参数符号更改。此外，多重共线性的存在增加了系数估计的方差或标准误差，使其对微小变化敏感，从而导致解释困难。

当我们讨论数据集的最佳回归模型时，这意味着模型会考虑预测变量，其中每个预测变量与目标/因变量相关，但相互之间几乎没有关联。这种模型被称为低噪声模型，具有很强的鲁棒性。如果变量本身之间的相关性很高，并且我们在训练数据集上获得了很高的预测能力，那么我们可能无法实现在测试数据集上相同。

2.检测多重共线性

在了解了多重共线性如何影响我们的分析和模型的预测能力之后，学习如何检测数据中多重共线性的存在非常重要。

以下方法可用于检测多重共线性：

• R平方值很大，但没有β重量统计上显着的，即整体模型的F检验是显着的，但个体系数估计的t检验不是。

• 变量对之间的相关性很大。

• 差异通货膨胀因素。

3. 方差膨胀因子

在检测多重共线性时提到的一点是，变量对之间的相关性很大，但有时只关注预测变量对之间的相关性可能是一个限制因素。

成对相关可能很小，但三个或更多变量之间存在线性相关。例如，X1，X2和X3之间的关系由线性方程定义，例如：X3 = 3X1 + 4X2 +Error。对于这种情况，我们使用方差膨胀因子。

VIF衡量估计回归系数的方差，可以通过模型中预测变量之间相关性的存在来夸大。如果VIF的值> 10，则表明变量之间存在一些严重的相关性。

我们不会讨论VIF的数学推导。 VIF的公式为：

预测变量“k”的VIF是通过对剩余预测变量上的第k个变量进行回归来计算R平方值而获得的。

多重共线性与高纬度的补救措施包括如下：

• 检查其中一个变量是否重复。

• 删除冗余变量。

• 通过收集更多数据来增加样本量。

• 平均中心预测变量。

• 预测变量的标准化 – 如果均值中心没有影响。

• 主成分分析，岭回归，偏最小二乘回归。

在本文档中，我们将借助示例讨论主成分分析。

二、主成分分析详解

PCA也称为特征缩减或特征提取PCA使用数据集中已存在的许多功能创建新功能。在PCA中，不是直接在解释变量上回归因变量，而是将解释变量的主成分用作回归量。人们通常只使用所有主成分的子集进行回归，从而使PCR成为一种正则化程序。通常选择具有较高方差的主要成分作为回归量。

PCA的一个主要用途在于克服多重共线性问题。PCA可以通过在回归步骤中排除一些低方差主成分来恰当地处理这种情况。此外，通过仅对所有主成分的子集进行回归，PCR可以通过显着降低表征基础模型的参数的有效数量来减少尺寸。

1. 主成分分析的定义

假设，最初在您的模型中有p个预测变量，每个都有n条记录。所以，数据矩阵如下所示：

2. 主成分分析的简单原理

假设在最后一步中生成的个新变量是它背后有两个条件，

• 新变量将是不相关的（或正交的）。

• 最大的变化将沿着，第二大变化将沿着，第三大变化将沿着，依此类推。

关于主成分的原理其实很简单，其本质上是数据的线形转换，首先使用用表示$X_1,X_2, , 的方差 – 协方差矩阵。

3. 主成分分析的算法

三、使用R语言进行主成分分析

在R中有很多方法可以实现主成分分析，基础的prcomp函数，psych包中的principal函数都可以构建主成分分析。在进行主成分分析之前，需要判断是需要选取多少个主成分。判断的准则有很多，包括：

• 根据经验或者理论知识进行判断主成分个数。

• 根据解释变量方差的积累值的阈值来判断保留的主成分个数。

• 通过变量间的关系矩阵来判断需要保留的主成分个数。

• 基于特征值进行判断，根据Kaiser-Harris准则建议保存特质大于1的主成分。

• 使用Cattell 碎石图进行判断，碎石图绘制了特征值与主成分数量，这类图可以清晰的展示特征值与主成分数之间的关系，图形变化最大之前的主成分都可以保留。

• 平行分析法，其原理是模拟一个与原数据集相同大小的矩阵来判断提取的特征值，若真实的某个特征值大于随机数据矩阵的平均特征值，则可以保留。

利用psych包中的fa.parallel函数可以对，特征值大于1，碎石检验，平行分析三种准则进行评价。所使用到的数据集是psych包中自带的一个数据集，Thurstone。使用fa.parallel构建碎石图的参数是需要进行主成分分析的数据集。碎石图的结果如图5.1所示。

library(psych)
fa.parallel(Thurstone)
## Warning in fa.parallel(Thurstone): It seems as if you are using a
## correlation matrix, but have not specified the number of cases. The number
## of subjects is arbitrarily set to be 100
## Parallel analysis suggests that the number of factors =  3  and the number of components =  1

图1碎石图的结果给出了三种准则的评判结果，从结果中可以发现，特征值大于1的特征包含三个；碎石图的曲线从第一个特征到第二个特征之间有急剧的变化，因此，选取第一个主成分；根据平行分析法，有一个特征大于随机矩阵的特征，总而言之，选取一个主成分最合适。

1. 主成分分析的实现

使用prcomp函数和psych包的principal函数都可以实现主成分分析。对于prcomp而言，传入需要进行主成分分析的原始数据即可得到主成分分析的结果：

com <- prcomp(Thurstone)
com
## Standard deviations (1, .., p=9):
## [1] 4.068067e-01 3.717632e-01 1.747973e-01 1.614655e-01 1.328383e-01
## [6] 1.137641e-01 8.417453e-02 5.966690e-02 5.671327e-17
## 
## Rotation (n x k) = (9 x 9):
##                           PC1         PC2        PC3         PC4
## Sentences         -0.51809588  0.01904674  0.2219406 -0.10347521
## Vocabulary        -0.50023933  0.09586015  0.1948562 -0.03929555
## Sent.Completion   -0.51043041  0.05906835  0.1922287 -0.05061761
## First.Letters      0.17889726  0.41142018  0.2781580  0.28375796
## Four.Letter.Words  0.22201272  0.34066448  0.3832537  0.47143714
## Suffixes           0.06316668  0.48794084 -0.5349420 -0.38650203
## Letter.Series      0.12799298 -0.46305936  0.1213871  0.04342013
## Pedigrees         -0.18269477 -0.36539660 -0.4557516  0.57758413
## Letter.Group       0.29283881 -0.34002238  0.3780991 -0.44445752
##                           PC5          PC6         PC7         PC8
## Sentences          0.06248519 -0.111471957  0.40722485 -0.68008332
## Vocabulary         0.09218721 -0.016188107  0.42092377  0.71770206
## Sent.Completion   -0.10617941  0.099931122 -0.76715103  0.04517829
## First.Letters      0.05821963  0.752293096  0.13448197 -0.06117109
## Four.Letter.Words -0.04884194 -0.625099183 -0.02856458  0.03809150
## Suffixes           0.12836566 -0.114136797  0.06024721  0.04481668
## Letter.Series      0.77851211  0.004042136 -0.05517789  0.05033184
## Pedigrees         -0.32964974  0.079329046  0.14775329  0.03158129
## Letter.Group      -0.48920364  0.035988647  0.14513339  0.09806774
##                           PC9
## Sentences         0.163013221
## Vocabulary        0.004977171
## Sent.Completion   0.290922025
## First.Letters     0.222883485
## Four.Letter.Words 0.264802257
## Suffixes          0.535941908
## Letter.Series     0.375364298
## Pedigrees         0.392434091
## Letter.Group      0.432241146

从结果中，得到了所有的主成分，然后根据之前的判断的主成分个数，选取相应的主成分，即可得到所要的结果。

另外一种方式是使用principal进行主成分分析，其可以通过原始数据矩阵或者相关系数矩阵进行主成分分析，其主要的参数如下：

• r：原始数据矩阵或者相关系数矩阵。

• nfactors：选取的主成分个数，默认为1。

然后，使用第二种方式进行主成分分析：

com <- principal(r = Thurstone,nfactors = 1)
com
## Principal Components Analysis
## Call: principal(r = Thurstone, nfactors = 1)
## Standardized loadings (pattern matrix) based upon correlation matrix
##                    PC1   h2   u2 com
## Sentences         0.82 0.67 0.33   1
## Vocabulary        0.84 0.70 0.30   1
## Sent.Completion   0.80 0.65 0.35   1
## First.Letters     0.72 0.52 0.48   1
## Four.Letter.Words 0.71 0.51 0.49   1
## Suffixes          0.67 0.45 0.55   1
## Letter.Series     0.67 0.45 0.55   1
## Pedigrees         0.70 0.50 0.50   1
## Letter.Group      0.64 0.41 0.59   1
## 
##                 PC1
## SS loadings    4.85
## Proportion Var 0.54
## 
## Mean item complexity =  1
## Test of the hypothesis that 1 component is sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.12 
## 
## Fit based upon off diagonal values = 0.94

接下来分析结果，PC1栏是成分荷载（componentloadings），表示的是观测变量与主成分的相关系数，如果不止提取一个主成分，那么会有PC2栏，或者PC3栏目。成分荷载可以用来解释成分的含义，从结果中可以看出，PC1与所有的变量都非常相关，因此，这个主成分是可以用来进行一般性评价的纬度。h2是公因子方差，表示主成分对于每一个变量的解释程度。u2表示唯一性，指的是变量无法被主成分解释的部分。SSloadings 包含了与主成分相关联的特征值。Proportion Var表示每个主成分对整个数据集的解释程度。从结果中可以看出，第一个主成分解释了所有变量54%的方差。

2. 主成分分析案例

在这里，这个例子将使用USArrests数据，该数据包含四个连续变量和50条记录。首先创建一个散点图矩阵来检查变量之间的关系，首先加载USArrests数据集，glimpse用于查看数据的结构，apply函数的使用方式类似于for循环，只是其效率更高，apply(USArrests,2,mean)表示计算每一列的均值，apply的第一个参数需要处理的数据集，第二个参数用于指定是对行进行处理还是对列进行处理，2表示对列进行处理，第三个参数用于指定使用函数进行处理。

data("USArrests")
head(USArrests,3)
##         Murder Assault UrbanPop Rape
## Alabama   13.2     236       58 21.2
## Alaska    10.0     263       48 44.5
## Arizona    8.1     294       80 31.0
library(dplyr)
 #Explore the structure of the data
usarrests <- tbl_df(USArrests)
glimpse(usarrests)
## Observations: 50
## Variables: 4
## $ Murder   <dbl> 13.2, 10.0, 8.1, 8.8, 9.0, 7.9, 3.3, 5.9, 15.4, 17.4, 5…
## $ Assault  <int> 236, 263, 294, 190, 276, 204, 110, 238, 335, 211, 46, 1…
## $ UrbanPop <int> 58, 48, 80, 50, 91, 78, 77, 72, 80, 60, 83, 54, 83, 65,…
## $ Rape     <dbl> 21.2, 44.5, 31.0, 19.5, 40.6, 38.7, 11.1, 15.8, 31.9, 2…
apply(USArrests,2,mean)
##   Murder  Assault UrbanPop     Rape 
##    7.788  170.760   65.540   21.232
apply(USArrests,2,var) #
##     Murder    Assault   UrbanPop       Rape 
##   18.97047 6945.16571  209.51878   87.72916

每个变量的方差都有很大差异。在PCA中，平均值不起作用，但方差在定义PC中起主要

normalize <- function(x)
{
  return((x - mean(x))/sd(x))
}

# create normalized variables
usarrests2 <- usarrests %>%
  apply(2, normalize) %>%
  as.data.frame()

接下来对数据进行探索性分析，分析特征的分布以及特征之间的一个相关性，这里使用的是GGally的ggpairs函数，其可以绘制出特征的分布图，特征之间的散点图，以及特征的相关系数。ggpairs的第一个参数是数据集，第二个参数是columns用于指定选择哪些变量进行绘图。特征相关性结果如图2所示。

library(GGally)
## Loading required package: ggplot2
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
## 
## Attaching package: 'GGally'
## The following object is masked from 'package:dplyr':
## 
##     nasa
ggpairs(data = usarrests2, 
        columns = 1:4,progress = F)

从图2看出，Murder与Assault，Assault与Rape，Murder与Rape，UrbanPop与Rape都是高度相关的。因此，数据存在比较严重的共线性问题。然后判断主成分个数，结果的碎石图如图3所示。

fa.parallel(usarrests2)
## Parallel analysis suggests that the number of factors =  2  and the number of components =  1

从碎石图来看，可以选取一个主成分。

然后分别使用进行主成分分析：

pca.out<-prcomp(usarrests2)
pca.out
## Standard deviations (1, .., p=4):
## [1] 1.5748783 0.9948694 0.5971291 0.4164494
## 
## Rotation (n x k) = (4 x 4):
##                 PC1        PC2        PC3         PC4
## Murder   -0.5358995  0.4181809 -0.3412327  0.64922780
## Assault  -0.5831836  0.1879856 -0.2681484 -0.74340748
## UrbanPop -0.2781909 -0.8728062 -0.3780158  0.13387773
## Rape     -0.5434321 -0.1673186  0.8177779  0.08902432
summary(pca.out)
## Importance of components:
##                           PC1    PC2     PC3     PC4
## Standard deviation     1.5749 0.9949 0.59713 0.41645
## Proportion of Variance 0.6201 0.2474 0.08914 0.04336
## Cumulative Proportion  0.6201 0.8675 0.95664 1.00000