R语言数据分析与挖掘(第六章):主成分分析(2)——案例讲解

  • 2019 年 12 月 16 日
  • 筆記

这一讲通过一个案例讲解主成分分析。

在R中用于完成主成分分析的函数是princomp()，该函数有2种调用方式：

1.公式形式

基本语法为：

princomp(formula, data = NULL, subset, na.action, ...)

参数介绍：

formula：指定用于主成分的公示对象，类似于回归分析和方差分析中的公式对象；

data：指定用于主成分分析的数据对象，一般为数据框；

subset：指定可选的向量，表示选择的样本子集；

na.action：一个函数，指定缺失数据的处理方法，若为NULL，则使用函数na.omit()删除缺失数据。

2. 矩阵形式

基本语法为：

princomp(x, cor = FALSE, scores = TRUE, covmat = NULL,           subset = rep_len(TRUE, nrow(as.matrix(x))), fix_sign = TRUE, ...)

参数介绍：

X:指定用于主成分分析的数据对象，一般为数值矩阵的数据框:

Cor:逻辑值，指定用于主成分分析中采用的矩阵形式(相关矩阵或协方差阵),为TRUE表示用样本的相关矩阵做主成分分析，为TALSE表示用样本的协方差阵做主成分分析,默认值为FALSE;

Scores:逻辑值，指定是否计算各主成分的分量，即是否计算样本的主成分得分，默认值TRUE;

Covmati:指定协方差阵，或者为cov.wt()提供的协方差列表，当数据不用参数x提供时，可由协方差阵提供;

Subset:指定可选向量，表示选择的样本子集。

函数princomp()的返回值为-个列表，包括:

sdev表示各主成分的标准差;

loadings表示载荷矩阵;

center表示各指标的样本均值;

scale表示各指标的样本标准差;

n.obs表示观测样本的个数；

scores表示主成分得分(当参数scores= TRUE时提供)。

此外，也可以利用其他函数来提取主成分分析的结果:

　　函数summary()可用于提取主成分的信息;

　　函数loadings()可用于提取载荷矩阵;

　　函数predict()可用于计算主成分得分;

　　函数screeplot()可用于绘制主成分的碎石图;

　　函数biplot()可用于绘制数据关于主成分的散点图和原坐标在主成分下的方向。

综合案例:longley数据集的变量降维及回归

下面利用第5章提到的longley数据集进行操作演练，该数据集收集了1947年至1962年16年的宏观经济数据，包含7个变量，分别为物价指数平减后的国民生产总值(GNP.deflator)、 国民生产总值(GNP)、(Amed.Frees)人口总量(Popultin)、年份(Year) 和就业人数(Emplyed)。失业人数(Unemployed)、 军队人数（Amed.Forces），人口总(Population）、年份(Year)和就业人数（Employed）。

> data(longley)  > summary(longley)    GNP.deflator         GNP          Unemployed     Armed.Forces   Min.   : 83.00   Min.   :234.3   Min.   :187.0   Min.   :145.6   1st Qu.: 94.53   1st Qu.:317.9   1st Qu.:234.8   1st Qu.:229.8   Median :100.60   Median :381.4   Median :314.4   Median :271.8   Mean   :101.68   Mean   :387.7   Mean   :319.3   Mean   :260.7   3rd Qu.:111.25   3rd Qu.:454.1   3rd Qu.:384.2   3rd Qu.:306.1   Max.   :116.90   Max.   :554.9   Max.   :480.6   Max.   :359.4     Population         Year         Employed   Min.   :107.6   Min.   :1947   Min.   :60.17   1st Qu.:111.8   1st Qu.:1951   1st Qu.:62.71   Median :116.8   Median :1954   Median :65.50   Mean   :117.4   Mean   :1954   Mean   :65.32   3rd Qu.:122.3   3rd Qu.:1958   3rd Qu.:68.29   Max.   :130.1   Max.   :1962   Max.   :70.55  > corr<-cor(longley)  > library(corrplot)  corrplot 0.84 loaded  > corrplot(corr)

该数据的相关矩阵如上图，可以很直观的看出，数据存在多重共线性问题，下面计算方差膨胀因子(VIF)来判断个变量的共线性问题的严重性：

> library(car)  > vif(lm(Employed~.,data=longley))  GNP.deflator          GNP   Unemployed Armed.Forces   Population     135.53244   1788.51348     33.61889      3.58893    399.15102          Year     758.98060

一般认为：

当0<VIF<10，不存在多重共线性；

当10<=VIF<100,存在较强的多重共线性；

当VIF>= 100，多重共线性非常严重。

上诉输出结果显示，变量GNP.deflator 、GNP、Population和Year存在严重的多重共线性。 为了解决多重共线性问题，下面对数据集进行主成分分析，去掉响应变量后进行分析：

> (pr1<-princomp(longley[,-7], cor = TRUE))  Call:  princomp(x = longley[, -7], cor = TRUE)    Standard deviations:      Comp.1     Comp.2     Comp.3     Comp.4     Comp.5     Comp.6  2.14554820 1.08413122 0.45102702 0.12218125 0.05051797 0.01940897     6  variables and  16 observations.  > summary(pr1)  Importance of components:                            Comp.1    Comp.2     Comp.3      Comp.4  Standard deviation     2.1455482 1.0841312 0.45102702 0.122181253  Proportion of Variance 0.7672295 0.1958901 0.03390423 0.002488043  Cumulative Proportion  0.7672295 0.9631196 0.99702383 0.999511871                               Comp.5       Comp.6  Standard deviation     0.0505179747 1.940897e-02  Proportion of Variance 0.0004253443 6.278469e-05  Cumulative Proportion  0.9999372153 1.000000e+00

输出结果中：Standard deviation表示主成分的标准差，即特征值的开方；Proportion of Variance表示每个主成分的贡献率；Cumulative Proportion表示主成分的累计贡献率，通常情况下，参考该参数的值进行主成分的选择，一般累计贡献率超过85%的主成分，在本案例中，前2个主成分的累计贡献率已经达到96%，即前2个主成分能解释原始变量的96%的信息，故应该选择前2个主成分，剔除后面4个主成分，达到降维的目的。

下面利用函数loadings()输出载荷矩阵：

> loadings(pr1)    Loadings:               Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6  GNP.deflator  0.462         0.149  0.793  0.338  0.135  GNP           0.462         0.278 -0.122 -0.150 -0.818  Unemployed    0.321 -0.596 -0.728               -0.107  Armed.Forces  0.202  0.798 -0.562  Population    0.462         0.196 -0.590  0.549  0.312  Year          0.465         0.128        -0.750  0.450                   Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6  SS loadings     1.000  1.000  1.000  1.000  1.000  1.000  Proportion Var  0.167  0.167  0.167  0.167  0.167  0.167  Cumulative Var  0.167  0.333  0.500  0.667  0.833  1.000

可以得到|成分!原的变量的线性关系:

第1主成分：

F1=0.462*CNPdnator-0.462*GNP-0.321*Unemploycd-0.202*Amed.Forces-0.462*Ppulation-0.465*Year

第2主成分

Fn=-0.596*Unemployed+0.798* Armed.Forces

下面对选样的主成分进行命名，主成分命名的原则应注意以下几点:

由于F1中包含不同变量的信息，在进行命名时需费综合考虑不同变量的信息；需要注意线性关系中系数的符号和绝对值大小。

基于上述原则，第一主成分对应系数的符号相同，均为负号，除变最Unemployed和Amed.Forces以外，其他原始变量的系数绝对值均为0.5左有，故第1主成分反映了经济的萧条程度，经济越萧条，F1的取值越大，各原始变量的取值越小，反之，经济越不萧条，F1的取值越小，各原始变量的取值越大。

第二主成分要与原始变量Unemployed和Amed.Forces有关,且变Unemployed对应的系数为负，Amed.Forces对应的系数为正，系数的绝对值均大于0.5,故第二主成分反映了军备程度，军备程度越高，F2的取值越大，失业人数越小，且军队人数越多。反之，军备程度越低，F2的取值越小，失业人数越多，军队人数越少。

　　尽管上述分析对两个主成分进行了命名，但命名结果还是不能很好地解释原版始变量的信息，因为主成分命名本身就具有较人的难度，我们能做到的就是尽量个面地涵盖原始变量的信息。

下面是利用函数predict()输出主成分的得分。

> predict(pr1)            Comp.1     Comp.2      Comp.3       Comp.4       Comp.5  1947 -3.59294203 -0.7761110  0.31804507 -0.169624216  0.009085721  1948 -3.10924187 -0.8768853  0.66329350  0.130051755  0.063565251  1949 -2.42014988 -1.5905016 -0.50961596 -0.009106066  0.005934717  1950 -2.16257276 -1.3181781 -0.11494311 -0.063265342 -0.063873551  1951 -1.48540767  1.2763227 -0.03004947  0.100660376  0.053969833  1952 -1.04339031  1.9851409 -0.16650384  0.047914580  0.038252420  1953 -0.72546382  1.9732212  0.06933769 -0.073217347  0.022425386  1954  0.03355589  0.6124972 -1.07307182 -0.067293112 -0.002331281  1955  0.10277563  0.7162362 -0.10076710 -0.104426569 -0.102049470  1956  0.46417097  0.5658113  0.30255566  0.018137849 -0.086508832  1957  0.98638447  0.4435320  0.45984009  0.123243732 -0.024470311  1958  1.87669233 -0.8914800 -0.69963411  0.193195372  0.022381351  1959  2.00361238 -0.3992485  0.27468475  0.148640632 -0.037890240  1960  2.43855615 -0.5154723  0.37766452  0.063619308 -0.016767727  1961  3.17896368 -1.0224220 -0.20859150 -0.070338720  0.058280241  1962  3.45445683 -0.1824626  0.43775561 -0.268192230  0.059996493             Comp.6  1947  0.002663801  1948  0.012373058  1949  0.005226048  1950 -0.014125087  1951 -0.044081763  1952 -0.013169060  1953  0.032810005  1954  0.017241531  1955 -0.019544921  1956  0.014604365  1957  0.028044545  1958  0.008369752  1959 -0.024302116  1960  0.004500429  1961 -0.001372125  1962 -0.009238462  > sort(predict(pr1)[,1],decreasing = T)         1962        1961        1960        1959        1958        1957   3.45445683  3.17896368  2.43855615  2.00361238  1.87669233  0.98638447         1956        1955        1954        1953        1952        1951   0.46417097  0.10277563  0.03355589 -0.72546382 -1.04339031 -1.48540767         1950        1949        1948        1947  -2.16257276 -2.42014988 -3.10924187 -3.59294203  > sort(predict(pr1)[,2],decreasing = T)        1952       1953       1951       1955       1954       1956   1.9851409  1.9732212  1.2763227  0.7162362  0.6124972  0.5658113        1957       1962       1959       1960       1947       1948   0.4435320 -0.1824626 -0.3992485 -0.5154723 -0.7761110 -0.8768853        1958       1961       1950       1949  -0.8914800 -1.0224220 -1.3181781 -1.5905016

从第一主成分看，1947、1948和1949年的得分较高，说明那几年的经济较为萧条，得分较低的是1960、1961和1962年，说明那几年的经济形势较好。

从第二主成分来看，1952、1953和1951年得分较高，说明那几年的军备情况较好，1961、1950和1949年的得分较低，说明那几年的军备情况较差。

下面利用函数是绘制主成分碎石图：

screeplot(pr1,type="lines",main="")

绘图形势选择折线图，碎石图展现的是各主成分没有解释道原始数据的信息情况，折线图越靠近x轴，说明主成分解释原始变量的信息越多。可以看出选择前2个主成分是合理的。

下面利用函数biplot()绘制双坐标图，默认情况下，该函数绘制第一、二主成分样本散点图和原始坐标在第一、二主成分下的方向。

biplot(pr1)

主成分回归 主成分回归可以参考前面讲的回归分析看代码：

pre1<-predict(pr1)  longley$F1<-pre1[,1]  longley$F2<-pre1[,2]  longley$F3<-pre1[,3]  lm.pr1<-lm(Employed~F1+F2+F3,data=longley)  summary(lm.pr1)  beta<-coef(lm.pr1)  A<-loadings(pr1)[,1:3]  x.mean<-pr1$center  x.sd<-pr1$scale  coef<-A%*%beta[2:4]/x.sd  beta0<-beta[1]-x.mean%*%coef  c(beta0,coef)

VirMach 便宜 VPS

QNews

QNews