在社会经济、教育和心理等领域会经常用到主成分分析，其目的是对数据进行降维，用少数几个（综合）指标的信息来反映原有的多个指标的绝大多数信息。
主成分分析的概念
主成分分析（Principal Component Analysis，简称PCA）是将多个指标化为少数几个综合指标的一种统计分析方法，即通过降维技术把多个变量化为几个少数主成分的方法
主成分分析的数学定义是，一个正交化线性变换，把数据变换到一个新的坐标系中，使这个数据的任何投影的第一大方差在第一个坐标（第一主成分）上，第二大方差在第二个坐标（第二主成分）上，以此类推。
(相关资料图)
主成分分析的基本思想是：将原来众多具有一定相关性的指标，重新组合成一组新的相互无关的综合指标来代替原来的指标。具体来说就是，定义p个随机变量X1, X2, ……, Xp，主成分分析就是要把p个指标的问题转变为m个新指标F1, F2, ……, Fm（m<=p），这m个新指标保留了原有变量的主要信息，并且
相互独立（正交）。
为了能更好地理解主成分分析的思想，我们来举个栗子:)
[例]今有14名学生的身高与体重数据，做相关图以显示变量间的关系
生成数据
x1 <- c(147,171,175,159,155,152,158,154,164,168,166,159,164,177) #生成向量x1x2 <- c(32,57,64,41,38,35,44,41,54,57,49,47,46,63) #生成向量x2绘制图形
plot(x1,x2,xlim=c(145,180),ylim=c(25,75))library(shape) #调用R包lines(getellipse(24,3,mid=c(162,48),angle=48),lty=3) #添加椭圆边框lines(c(146,178),c(30,66)); #主成分1lines(c(162,165),c(53,46)); #主成分2求协方差
cov(x1, x2)## [1] 85.86813求相关系数
cor(x1, x2)## [1] 0.9672073主成分分析的目的
主成分分析的主要目的是用较少的变量去解释原资料中的大部分变异（方差），即期望能将手中许多相关性很高的变量转化成互相独立（正交）的变量，并能解释大部分资料之变异的几个新变量，也就是所谓的主成分。
表示两个变量的相关性，在数学上用协方差阵，在统计学上用相关系数阵
X <- data.frame(x1, x2)S <- cov(X); S #协方差阵## x1 x2## x1 77.45604 85.86813## x2 85.86813 101.75824R <- cor(X); R #相关系数阵## x1 x2## x1 1.0000000 0.9672073## x2 0.9672073 1.0000000主成分分析的数学表达
主成分分析的数学表达
从1到p主成分所能承载的信息量应该是递减的
主成分的推导
要解决这个问题就要使用到数学上的谱分解定理：
基于这一定理我们就可以进行主成分的推导了
从公式的结果我们可以看出当a=u1的时候主成分的方差最大，所以：
这样我们就完成了主成分分析的数学推导过程
主成分的性质
在完成了主成分的数学推导之后，我们就来看看主成分的性质
看到这里的第四个性质，虽然有点抽象但也是很重要的，它表达了方差相等的意思，说明我们进行主成分变换之后并没有改变原来数据的信息量。同样抽象的还有第二个式子，它表明了主成分和原来变量之间的相关系数就是主成分的载荷。
主成分分析的R语言实现
在R中，计算主成分分析的函数是princomp()函数。其具体的用法为：princomp(x, cor = FALSE, scores = TRUE, …) 其中的参数意义为：x 数据矩阵或数据框 cor 是否用相关阵（常用），默认为协差阵 scores 是否输出成分得分
用之前生成的身高体重数据框来做
pc <- princomp(X); pc## Call:## princomp(x = X)## ## Standard deviations:## Comp.1 Comp.2 ## 12.795925 1.636326 ## ## 2 variables and 14 observations.结果中就显示了两个主成分的标准差， 但这里用标准差不太好，所以可以转化为方差
pc$sdev^2## Comp.1 Comp.2 ## 163.735703 2.677562在此基础之上我们可以计算主成分的负荷（loadings）（也就是系数） 主成分负荷：
pc$loadings## ## Loadings:## Comp.1 Comp.2## x1 0.656 0.755## x2 0.755 -0.656## ## Comp.1 Comp.2## SS loadings 1.0 1.0## Proportion Var 0.5 0.5## Cumulative Var 0.5 1.0有了主成分负荷（系数）就可以把它带入y=a’x中计算出主成分得分（scores）
pc$scores## Comp.1 Comp.2## [1,] -21.7469609 -1.0753729## [2,] 12.8653792 0.6526072## [3,] 20.7733282 -0.9172151## [4,] -7.0833638 2.0835700## [5,] -11.9712295 1.0305890## [6,] -16.2033944 0.7326289## [7,] -5.4740021 -0.6385533## [8,] -10.3618678 -1.6915342## [9,] 6.0104111 -2.6654363## [10,] 10.8982768 -1.6124553## [11,] 3.5467085 2.1231094## [12,] -2.5532388 -1.8506348## [13,] -0.0297556 2.5801701## [14,] 21.3297090 1.2485274然后绘制一个图
plot(pc$scores, asp = 1);#asp调整x轴和y轴的比例，一定要加；abline(h=0, v=0, lty=3)从图中可以看出，通过主成分的变换后本来相关的x变量变成不相关的主成分——comp.1和comp.2，而且从图中可以看出，数据的信息量主要集中在comp.1上（图中的点主要集中在comp.1上），即第一个主成分之上。
可以进一步用biplot图来表示 首先可以用biplot()函数来绘图，用法：biplot(scores, loadings, …) scores得分 loading载荷
biplot(pc$scores, pc$loadings)abline(h=0, v=0, lty=3)也可以使用ggbiplot()带有ggplot2风格的函数绘图
#library(devtools)#install_github("vqv/ggbipolt")library(ggbiplot)ggbiplot(pc,obs.scale = 1, var.scale = 1, ellipse = TRUE) + scale_y_continuous(limits = c(-20, 20))两张图中x1和x2的位置不同主要是由于第一张图有两个刻度，显然第二张图更能说明主成分的压缩，x1和x2都被压缩到了第一主成分上。
在数学上可以通过计算方差贡献率来进行压缩
summary(pc)可以到，两个变量可以压缩为一个主成分comp.1，这个主成分的方差贡献率高达98%（一般在80%以上就可）
主成分分析的步骤
根据上面的分析，我们可以对主成分分析的一般步骤做一个总结 1、将原始数据标准化，得标准化数据矩阵；2、建立相关系数阵；3、求特征值及特征向量；4、获得主成分
主成分分析函数算法的改进
本文的写作参考了暨南大学王斌会教授的《多元统计分析及R语言建模》一书以及相关的mooc课程，王斌会教授在princomp()函数的基础之上，改进并编写了新的函数，使得能够在计算主成分时忽略掉特征值特别小的主成分，使得压缩更有效率。算法改进主要用到的公式：
改进之后的函数：
msa.pca<-function(X,cor=FALSE,m=2,scores=TRUE,ranks=TRUE,sign=TRUE,plot=TRUE){ if(m<1) return PC=princomp(X,cor=cor) Vi=PC$sdev^2 Vari=data.frame("Variance"=Vi[1:m],"Proportion"=(Vi/sum(Vi))[1:m], "Cumulative"=(cumsum(Vi)/sum(Vi))[1:m]) cat("\n") Loadi=as.matrix(PC$loadings[,1:m]) Compi=as.matrix(PC$scores[,1:m]) if(sign) for (i in 1:m) if(sum(Loadi[,i])<0){ Loadi[,i] = -Loadi[,i] Compi[,i] = -Compi[,i] } pca<-NULL pca$vars=Vari if(m<=1) pca$loadings = data.frame(Comp1=Loadi) else pca$loadings = Loadi; if(scores & !ranks) pca$scores=round(Compi,4) if(scores & plot){ plot(Compi);abline(h=0,v=0,lty=3) text(Compi,row.names(X)) # par(mar=c(4,4,2,3)) # biplot(Compi,Loadi); abline(h=0,v=0,lty=3) # par(mar=c(4,4,1,1)) } if(scores & ranks){ pca$scores=round(Compi,4) Wi=Vi[1:m];Wi Comp=Compi%*%Wi/sum(Wi) Rank=rank(-Comp) pca$ranks=data.frame(Comp=round(Comp,4),Rank=Rank) } pca}简单介绍：该函数用于进行主成分分析（PCA）并输出相关的结果。具体功能如下：
输入参数：X为待分析数据矩阵，cor表示是否进行数据标准化后再进行PCA，m表示保留的主成分数目，scores表示是否输出每个样本在主成分上的投影分数，ranks表示是否输出每个样本在主成分上的排名，sign表示是否将负载荷矩阵和得分矩阵中的负值变为正值，plot表示是否绘制得分矩阵的散点图。
输出结果：返回一个列表，包含以下结果：
vars：每个主成分的方差、方差贡献率和累计方差贡献率。
loadings：每个变量对每个主成分的负载荷。
scores：每个样本在每个主成分上的投影分数。
ranks：每个样本在每个主成分上的排名。
该函数首先调用princomp函数进行PCA，并计算出每个主成分的方差、方差贡献率和累计方差贡献率。然后将负载荷矩阵和得分矩阵的负值变为正值（如果设置了sign参数），并将结果存储在loadings和scores变量中。如果设置了ranks参数，则还会计算每个样本在每个主成分上的排名，并将结果存储在ranks变量中。如果设置了plot参数，则会绘制得分矩阵的散点图。最后将所有结果存储在一个列表中并返回。
用法
可以使用这一个改进算法的函数来进行计算
msa.pca(X)#协方差## $vars## Variance Proportion Cumulative## Comp.1 163.735703 0.98391016 0.9839102## Comp.2 2.677562 0.01608984 1.0000000## ## $loadings## Comp.1 Comp.2## x1 0.6557008 0.7550208## x2 0.7550208 -0.6557008## ## $scores## Comp.1 Comp.2## [1,] -21.7470 -1.0754## [2,] 12.8654 0.6526## [3,] 20.7733 -0.9172## [4,] -7.0834 2.0836## [5,] -11.9712 1.0306## [6,] -16.2034 0.7326## [7,] -5.4740 -0.6386## [8,] -10.3619 -1.6915## [9,] 6.0104 -2.6654## [10,] 10.8983 -1.6125## [11,] 3.5467 2.1231## [12,] -2.5532 -1.8506## [13,] -0.0298 2.5802## [14,] 21.3297 1.2485## ## $ranks## Comp Rank## 1 -21.4144 14## 2 12.6689 3## 3 20.4243 2## 4 -6.9359 10## 5 -11.7620 12## 6 -15.9309 13## 7 -5.3962 9## 8 -10.2224 11## 9 5.8708 5## 10 10.6970 4## 11 3.5238 6## 12 -2.5419 8## 13 0.0122 7## 14 21.0066 1msa.pca(X, cor = T)#相关系数## $vars## Variance Proportion Cumulative## Comp.1 1.96720729 0.98360365 0.9836036## Comp.2 0.03279271 0.01639635 1.0000000## ## $loadings## Comp.1 Comp.2## x1 0.7071068 0.7071068## x2 0.7071068 -0.7071068## ## $scores## Comp.1 Comp.2## [1,] -2.3997 -0.1135## [2,] 1.4199 0.0690## [3,] 2.2626 -0.1067## [4,] -0.7445 0.2323## [5,] -1.2962 0.1170## [6,] -1.7646 0.0851## [7,] -0.6097 -0.0693## [8,] -1.1614 -0.1846## [9,] 0.6180 -0.2964## [10,] 1.1698 -0.1812## [11,] 0.4211 0.2340## [12,] -0.3080 -0.2041## [13,] 0.0361 0.2855## [14,] 2.3566 0.1328## ## $ranks## Comp Rank## 1 -2.3622 14## 2 1.3978 3## 3 2.2238 2## 4 -0.7285 10## 5 -1.2731 12## 6 -1.7343 13## 7 -0.6008 9## 8 -1.1454 11## 9 0.6030 5## 10 1.1476 4## 11 0.4180 6## 12 -0.3063 8## 13 0.0402 7## 14 2.3202 1可以看出改进算法后的函数能够自动绘制图形，输出结果也更加合理
以此我们可以总结出在R语言中进行主成分分析的一般步骤：
案例分析
学会了主成分分析，接下来让我们来进行一个案例分析——对居民消费数据做主成分分析
library(openxlsx)d3.1 <- read.xlsx("mvstats5.xlsx","d3.1",rowNames=T)head(d3.1, n = 10)计算相关系数矩阵
主成分分析法的前提是需要变量具有相关性，因此，我们需要先查看各个变量的相关程度
options(digits = 2)cor(d3.1)#输出相关系数矩阵可视化相关系数
library(corrplot)cor1 <- cor(d3.1)corrplot.mixed(cor1, lower = "ellips", upper = "number")从图中可以看出多个变量之间都存在正相关关系，有些相关关系很强，可以进行后续的主成分分析
首先用R中自带的princomp()函数进行主成分分析
PC <- princomp(d3.1, cor = T)确定主成分个数
通过特征值和特征向量计算方差贡献率和累计方差贡献率，并绘制碎石图
summary(PC)#计算方差贡献率和累计方差贡献率library(ggbiplot)ggscreeplot(PC, type = "pev") + geom_hline(yintercept = 0.1, linetype = "dashed")#绘制碎石图碎石图：可用作主成分分析，为什么叫“碎石”图，意思就是就是一颗石头从上面滚下来，只要取出让石头滚得快的点，取斜率比较大的点，就是该因素的主要因素，主要结合累计贡献率来得出取其中几个因素来作为主要因素，即达到降维的效果。
根据方差贡献率以及碎石图中可知，我们可以把主成分压缩为2维，即保留两个主成分
主成分负荷分析
PC$loadings从负荷矩阵上来看（只用分析前两个），主成分Comp.1在设备、交通、教育、居住、杂项上载荷值较大，可解释为非必须消费主成分；主成分Comp.2在食品、衣着、医疗上载荷值很大，可解释为为反映日常必须消费的主成分。这样通过负荷（系数）的分析我们对压缩的主成分进行了解释。
计算主成分得分
PC$scores[, 1:2]#只计算前两个## Comp.1 Comp.2## 北京 6.122 -1.52## 天津 3.010 -0.54## 河北 -0.888 -0.69## 山西 -1.104 -0.60## 内蒙古 0.533 -1.85## 辽宁 0.094 -0.66## 吉林 -0.327 -1.42## 黑龙江 -1.689 -1.00## 上海 7.085 1.07## 江苏 1.141 0.45## 浙江 3.821 -0.17## 安徽 -1.123 0.35## 福建 1.172 1.38## 江西 -1.669 0.55## 山东 0.481 -0.81## 河南 -1.277 -0.65## 湖北 -1.009 0.12## 湖南 -0.365 -0.20## 广东 4.032 2.48## 广西 -1.627 1.23## 海南 -1.873 2.35## 重庆 0.394 -0.46## 四川 -1.154 0.52## 贵州 -2.014 0.66## 云南 -2.430 0.42## 西藏 -2.720 1.01## 陕西 -0.888 -0.12## 甘肃 -1.325 -0.14## 青海 -1.768 -0.21## 宁夏 -1.317 -0.50## 新疆 -1.318 -1.05绘制主成分得分图
plot(PC$scores)abline(h=0,lty=3);abline(v=0,lty=3)text(PC$scores, labels=rownames(PC$scores))从图中可知，在主成分comp.1上得分最高的五个地区为上海、北京、广东、浙江、天津，且上海、北京的绝对值明显高于其他地区，表明在以设备、交通、教育、居住、杂项为主的消费上，上海、北京的消费水平远远高于其他省份。在主成分comp.2上得分最高的五个地区为广东、海南、福建、广西、上海，可见这部分地区对于食品、衣着、医疗方面的消费较大。
同时也可以用biplot图来可视化得分情况
library(ggbiplot)ggbiplot(PC, choices = 1:2, obs.scale = 1, var.scale = 1, ellipse = T, circle = T, varname.size = 4)+ geom_text(aes(x=PC$scores[,1], y=PC$scores[,2], label=rownames(PC$scores)))这个图中的中心点箭头：表示所有样本在主成分上的平均得分。该箭头通常是指向主成分得分的均值点，它可以用于帮助我们了解样本得分分布的位置和形态。
用改进的函数进行主成分分析
msa.pca(d3.1, cor = T)## $vars## Variance Proportion Cumulative## Comp.1 5.7 0.71 0.71## Comp.2 1.0 0.13 0.84## ## $loadings## Comp.1 Comp.2## 食品 0.35 -0.429## 衣着 0.25 0.677## 设备 0.37 -0.089## 医疗 0.30 0.472## 交通 0.38 -0.324## 教育 0.40 -0.069## 居住 0.37 -0.056## 杂项 0.37 0.118## ## $scores## Comp.1 Comp.2## 北京 6.122 1.52## 天津 3.010 0.54## 河北 -0.887 0.69## 山西 -1.104 0.60## 内蒙古 0.533 1.85## 辽宁 0.094 0.66## 吉林 -0.327 1.42## 黑龙江 -1.689 1.00## 上海 7.085 -1.07## 江苏 1.141 -0.45## 浙江 3.821 0.17## 安徽 -1.123 -0.35## 福建 1.172 -1.38## 江西 -1.669 -0.55## 山东 0.481 0.81## 河南 -1.277 0.65## 湖北 -1.010 -0.12## 湖南 -0.365 0.20## 广东 4.032 -2.48## 广西 -1.627 -1.23## 海南 -1.873 -2.35## 重庆 0.394 0.46## 四川 -1.154 -0.52## 贵州 -2.014 -0.66## 云南 -2.429 -0.42## 西藏 -2.720 -1.01## 陕西 -0.888 0.12## 甘肃 -1.325 0.14## 青海 -1.768 0.21## 宁夏 -1.317 0.50## 新疆 -1.318 1.05## ## $ranks## Comp Rank## 北京 5.419 2## 天津 2.632 5## 河北 -0.646 14## 山西 -0.843 16## 内蒙古 0.734 8## 辽宁 0.180 11## 吉林 -0.059 12## 黑龙江 -1.278 24## 上海 5.838 1## 江苏 0.897 6## 浙江 3.263 3## 安徽 -1.005 20## 福建 0.782 7## 江西 -1.498 26## 山东 0.531 9## 河南 -0.983 19## 湖北 -0.873 17## 湖南 -0.279 13## 广东 3.037 4## 广西 -1.567 27## 海南 -1.946 29## 重庆 0.405 10## 四川 -1.057 22## 贵州 -1.807 28## 云南 -2.122 30## 西藏 -2.459 31## 陕西 -0.734 15## 甘肃 -1.100 23## 青海 -1.466 25## 宁夏 -1.040 21## 新疆 -0.956 18改进之后的函数，能一次输出方差贡献率，负荷，得分以及得分图，最重要的是能够输出综合排名（rank）这是R语言原先的函数没有的。
以上就是主成分分析的所有内容了，建议以后用主成分分析可以直接使用王斌会教授所改进的函数，方便快捷。
