主成分分析和因子分析

#包载入
library(corrplot)
library(psych)
library(GPArotation)
library(nFactors)
library(gplots)
library(RColorBrewer)
1
2
3
4
5
6
7

主成分分析

主成分分析（PCA）是对针对大量相关变量提取获得很少的一组不相关的变量，这些无关变量也成为主成分变量。

数据探索

#本部分引入消费者品牌感知的问卷数据集作为数据pca <- read.csv("http://r-marketing.r-forge.r-project.org/data/rintro-chapter8.csv")summary(pca)
#可以看到，该数据除了一个字符变量外，其他都是1~10的数字变量
head(pca)
1
2
3
4
5
6
7
8

library(corrplot)
#检查两两变量之间的相关关系
corrplot(cor(pca[,1:9]), order = "FPC")
#"FPC" for the first principal component order.
1
2
3
4
5

可以大致发现聚成了三个类别，分别是latest/trendy/fun，serious/leader/perform，bargain/rebuy/value三类，而这也是接下来的分析所要验证的。

提取主成分

#数据表度化
pca.sc <- pca
pca.sc[,1:9] <- scale(pca.sc[,1:9])
#提取主成分
pca.pc <- prcomp(pca.sc[,1:9])
summary(pca.pc)
#判断成分数量
plot(pca.pc, type = "l")
1
2
3
4
5
6
7
8

对于主成分数量，scree图可看到在3类之后，每个主成分解释的方差增值减少。

主成分得分获取

#psych包也有不错的PCA分析输出
library(psych)
fa.parallel(pca.sc[,1:9],fa = "pc") #principal()需要事先知道大致多少成分，该包的碎石图输出
pca.psy <- principal(pca.sc[,1:9], nfactors = 3,rotate = "none")
round(unclass(pca.psy$weights),2) #获取主成分得分
1
2
3
4
5

此处的获得的各个主成分的构成系数，可以导出：PC1 = 0.14*perform + 0.12 leader + latest(-0.21)….

品牌感知图

与此同时，主成分分析运用的另一个重要方面是，通过双标图（biplot）考察不同类别（品牌）之间的关系可视化

#对主成分中的前两个成分映射到二维，但因为直接投射数据主体会面临散点太多，可视度差的问题
biplot(pca.pc)
#所以，可对其类别（品牌）进行映射
pca.mean <- aggregate(.~ brand, pca.sc, mean)
rownames(pca.mean) <- pca.mean[,1]
pca.mean.pc <- prcomp(pca.mean[,-1], scale = T)
summary(pca.mean.pc)
biplot(pca.mean.pc) #品牌感知图
1
2
3
4
5
6
7
8
9

从这张图上就可以考到各个分类（品牌）的调性和位置。

可进一步考察不同品牌间的差异：

pca.mean["a",] - pca.mean["j",]
1
2

解释性因子分析

因子分析（EFA）是用来发现一组变量潜在结构的方法，主要针对无法观测的因子变量，提取获得可观测的显变量。

EFA因子数量确定

# 碎石图和特征值确定因子数目
library(nFactors)
# 多种碎石图方案
nScree(pca.sc[,1:9])#  noc naf nparallel nkaiser
#  1   3   2         3       3# >1 的特征值数量
eigen(cor(pca.sc[,1:9]))#$values
#[1] 2.9792956 2.0965517 1.0792549 0.7272110 0.6375459 0.5348432 0.3901044
#[8] 0.3120464 0.2431469
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

从以上结果可知，因子数在2和3之间选择。

#使用psych包提取公因子
#公因子数量为2的方案
fa(pca.sc[,1:9], nfactors = 2, rotate = "none",fm = "ml") #因子化方法选择最大似然法ml
1
2
3

2因子方案的只解释了44%的方差。

#公因子数量为3的方案
fa(pca.sc[,1:9], nfactors = 3, rotate = "none",fm = "ml") 
1
2

上升到了57%的方差解释度。由此可认为3因子方案更胜一筹。

EFA旋转

两种方法进行旋转

#使用psych的fa来进行正交旋转
#正交旋转
fa.vaf <- fa(pca.sc[,1:9], nfactors = 3, rotate = "varimax",fm = "ml")  
fa.vaf
#斜交旋转
#使用GPArotation的oblimin斜交旋转，factanal的结果相对简练。
library(GPArotation)
fa.ob <- factanal(pca.sc[,1:9], factors = 3, rotation = "oblimin") 
fa.ob
#获得未列出的因子结构矩阵,通过因子模型矩阵*因子关联矩阵获得
fsm <- function(oblique) {
if (class(oblique)[2]=="fa" & is.null(oblique$rotmat)) {warning("Object doesn't look like oblique EFA")
} else {P <- unclass(oblique$loading)F <- P %*% oblique$rotmatreturn(F)
} }
fsm(fa.ob)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

正交旋转

斜交旋转

正交旋转将人为的强制因子间不相关，并且关注的结果集中在因子结构矩阵（correlations of the variables with the factors）。斜交旋转则关注三个矩阵：
* 因子模式矩阵（pattern matrix），结果中的（loading）列出了各个变量和因子变量的标准化回归系数。
* 因子关联矩阵（Factor Correlations matrix），因子间的关联性考察。
* 因子结构矩阵（structure matrix），即因子载荷矩阵，衡量变量与因子间的相关系数。

#相应的替代方案
#正交旋转
fa.va <- factanal(pca.sc[,1:9], factors = 3, rotation = "varimax") 
#斜交旋转
fa.promax <- fa(pca.sc[,1:9], nfactors = 3, rotate = "promax",fm = "ml") fsmfa <- function(oblique) {
if (class(oblique)[2]=="fa" & is.null(oblique$Phi)) {warning("Object doesn't look like oblique EFA")
} else {P <- unclass(oblique$loading)F <- P %*% oblique$Phireturn(F)
} }fsmfa(fa.promax)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

EFA旋转结果的可视化

使用路线图展示潜因子和单独因子之间的关系。

fa.diagram(fa.promax,simple = T,digits = 2)
1

使用热图更直观的展示潜因子和变量之间的关系。

library(gplots)
library(RColorBrewer)
heatmap.2(fa.ob$loadings,col=brewer.pal(9, "GnBu"), trace="none", key=FALSE, dend="none",Colv=FALSE, cexCol = 1.5,main="\n\n\nFactor loadings for brand adjectives")
1
2
3
4
5
6

不同品牌因子得分

#获得各个品牌的三个因子均值
fa.ob <- factanal(pca.sc[,1:9], factors = 3, rotation = "oblimin",scores = "Bartlett") 
fa.score <- data.frame(fa.ob$scores)
fa.score$brand <- pca.sc$brand
fa.score.mean <- aggregate(.~ brand, fa.score, mean)
fa.score.mean
#根据因子均值做热图
rownames(fa.score.mean) <- fa.score.mean[, 1]           # brand names
fa.score.mean <- fa.score.mean[, -1]heatmap.2(as.matrix(fa.score.mean),col=brewer.pal(9, "GnBu"), trace="none", key=FALSE, dend="none",cexCol=1.2, main="\n\n\n\nMean factor score by brand")
1
2
3
4
5
6
7
8
9
10
11
12
13
14