R语言数据分析与挖掘(第六章):主成分分析(2)——案例讲解
R语言数据分析与挖掘(第六章):主成分分析(2)——案例讲解
这一讲通过一个案例讲解主成分分析。
在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)本文分享自 MedBioInfoCloud 微信公众号，前往查看
如有侵权，请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体分享计划 ，欢迎热爱写作的你一起参与！