主成分分析和因子分析

#包載入 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

總結

以上是生活随笔為你收集整理的因子分析模型的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。