PCA and FA commands
PCA
1. prcomp() (stats)
prcomp(x, scale = FALSE)
princomp(x, cor = FALSE, scores = TRUE)
Arguments for princomp()
x: A numeric matrix or data frame.
cor: A logical value. If TRUE, then data will be centred and also scaled before the
analysis.
scores: A logical value. If TRUE, then coordinates on each principal component
are calculated.
1. prcomp() (stats)
2. princomp() (stats)
3. PCA() (FactoMineR)
4. dudi.pca() (ade4)
5. acp() (amap)
> library(FactoMineR) #Author DataFlair
> pca <- PCA(mtcars[,c(1:7,10,11)], scale. = TRUE)
> summary(pca)
> pca$eig #Author DataFlair
> pca$var$coord #DataFlair
> library(devtools)
> install_github("vqv/ggbiplot")
> library(ggbiplot) #AuthorDataFlair
> ggbiplot(pca)
Pricipal Components Analysis
# entering raw data and extracting PCs
# from the correlation matrix
fit <- princomp(mydata, cor=TRUE)
summary(fit) # print variance accounted for
loadings(fit) # pc loadings
plot(fit,type="lines") # scree plot
fit$scores # the principal components
biplot(fit)
Use cor=FALSE to base the principal components on the covariance matrix. Use
the covmat= option to enter a correlation or covariance matrix directly. If entering a
covariance matrix, include the option n.obs=.
The principal( ) function in the psych package can be used to extract and rotate
principal components.
Varimax Rotated Principal Components
# retaining 5 components
library(psych)
fit <- principal(mydata, nfactors=5, rotate="varimax")
fit # print results
mydata can be a raw data matrix or a covariance matrix. Pairwise deletion of missing data is
used. rotate can "none", "varimax", "quatimax", "promax", "oblimin", "simplimax", or "cluster"
Thực hành phân tích PCA trong R
Cài các packages sử dụng trong PCA: install.packages(c("FactoMineR",
"factoextra"))
Gọi các packages
library("FactoMineR")
library("factoextra")
## Loading required package: ggplot2
## Welcome! Want to learn more? See two factoextra-related books at
https://goo.gl/ve3WBa
Đọc dữ liệu vào R
Sử dụng data set Iris
data(iris)
head(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
Tách các biến định lượng thành một data set mới
df <- iris[, 1:4]
head(df)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 5.1 3.5 1.4 0.2
## 2 4.9 3.0 1.4 0.2
## 3 4.7 3.2 1.3 0.2
## 4 4.6 3.1 1.5 0.2
## 5 5.0 3.6 1.4 0.2
## 6 5.4 3.9 1.7 0.4
Sử dụng hàm PCA trong FactoMineR package
PCA(X, scale.unit = TRUE, ncp = 5, graph = TRUE)
“scale.unit = TRUE”: Một giá trị logic, dữ liệu sẽ được chia tỉ lệ thành đơn vị
phương sai trước khi phân tích. Việc tiêu chuẩn hóa theo cùng một thang
đo này tránh cho một số biến trở nên thống trị chỉ vì các đơn vị đo lường lớn
của chúng. Điều này giúp cho các biến có thể so sánh được.
library("FactoMineR")
pca <- PCA(df, graph = FALSE)
print(pca)
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 150 individuals, described by 4 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"
Trích xuất giá trị eigenvalue và các phương sai của các PC bằng factoextra
package
library(factoextra)
eig.val <- get_eigenvalue(pca)
eig.val
## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 2.91849782 72.9624454 72.96245
## Dim.2 0.91403047 22.8507618 95.81321
## Dim.3 0.14675688 3.6689219 99.48213
## Dim.4 0.02071484 0.5178709 100.00000
Tỉ lệ của các biến thiên được thể hiện bằng một giá trị eigenvalue ở cột thứ
2. VD: Dim.1 có eigen value là 2.918, tương ứng với tỉ lệ % phương sai là
72.96 (= 2.918/4)
Các giá trị eigen được sử dụng để xác định số lượng các thành phần chính
cần giữ lại sau PCA (Kaiser 1961):
Nếu giá trị eigen > 1 nói lên rằng các thành phần chính (PCs) chiếm nhi ều
phương sai hơn so với một trong các biến ban đầu. Đây thường được dùng
như một điểm giới hạn để xác định các PC được giữ lại.*
Chúng ta có thể giới hạn số lượng thành phần chính mà số đó chiếm một
phần nhất định của tổng phương sai (VD > 70%).
Trực quan tỉ lệ các thành phần chính bằng biểu đồ scree với các
hàm fviz_eig() hoặc fviz_screeplot() trong packages factoextra
fviz_eig(pca, addlabels = TRUE, ylim = c(0, 100))
Trực quan hoá kết quả PCA bằng biểu đồ với mỗi nhóm gắn với các ký hiệu và
màu sắc khác nhau.
fviz_pca_ind(pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = iris$Species, # color by groups
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
Phân tích PCA với các gói lệnh có sẵn trong R
Tính các thành phần chính PCA
myPr <- prcomp(iris[, -5])
myPr <- prcomp(iris[, -5], scale = TRUE)
myPr # Kiểm tra kết quả
## Standard deviations (1, .., p=4):
## [1] 1.7083611 0.9560494 0.3830886 0.1439265
##
## Rotation (n x k) = (4 x 4):
## PC1 PC2 PC3 PC4
## Sepal.Length 0.5210659 -0.37741762 0.7195664 0.2612863
## Sepal.Width -0.2693474 -0.92329566 -0.2443818 -0.1235096
## Petal.Length 0.5804131 -0.02449161 -0.1421264 -0.8014492
## Petal.Width 0.5648565 -0.06694199 -0.6342727 0.5235971
Thành phần chính đầu tiên có tương quan thuận với chiều dài đài hoa, chi ều
dài cánh hoa và chiều rộng cánh hoa (Ba biến này có mối tương quan cao
trong phân tích biểu đồ nhiệt phân cụm).
Chiều rộng đài hoa là biến số gần như giống nhau giữa ba loài với độ lệch
chuẩn nhỏ.
PC2 chủ yếu được xác định bởi chiều rộng lá đài (Sepal.Width), ít hơn bởi
chiều dài lá đài.
Lưu ý rằng “scale = TRUE” trong lệnh trên có nghĩa là dữ liệu được chuẩn hóa
trước khi phân tích PCA, do đó mỗi biến có đơn vị phương sai.
plot(myPr, ylim = c(0,4)) # Biểu đồ phương sai mỗi thành phần chính thu được
plot(myPr, type = "l") # Biểu đồ phương sai mỗi thành phần chính thu được
biplot(myPr)
biplot(myPr, scale = 0)
Tách điểm các thành phần chính (PC score)
str(myPr) # Kiểm tra cấu trúc đối tượng, liệt kê tất cả các thành phần
## List of 5
## $ sdev : num [1:4] 1.708 0.956 0.383 0.144
## $ rotation: num [1:4, 1:4] 0.521 -0.269 0.58 0.565 -0.377 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:4] "Sepal.Length" "Sepal.Width" "Petal.Length"
"Petal.Width"
## .. ..$ : chr [1:4] "PC1" "PC2" "PC3" "PC4"
## $ center : Named num [1:4] 5.84 3.06 3.76 1.2
## ..- attr(*, "names")= chr [1:4] "Sepal.Length" "Sepal.Width"
"Petal.Length" "Petal.Width"
## $ scale : Named num [1:4] 0.828 0.436 1.765 0.762
## ..- attr(*, "names")= chr [1:4] "Sepal.Length" "Sepal.Width"
"Petal.Length" "Petal.Width"
## $ x : num [1:150, 1:4] -2.26 -2.07 -2.36 -2.29 -2.38 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:4] "PC1" "PC2" "PC3" "PC4"
## - attr(*, "class")= chr "prcomp"
myPr$x # Giá trị toạ độ mới cho mỗi observation
## PC1 PC2 PC3 PC4
## [1,] -2.25714118 -0.478423832 0.127279624 0.024087508
## [2,] -2.07401302 0.671882687 0.233825517 0.102662845
## [3,] -2.35633511 0.340766425 -0.044053900 0.028282305
## [4,] -2.29170679 0.595399863 -0.090985297 -0.065735340
## [5,] -2.38186270 -0.644675659 -0.015685647 -0.035802870
## [6,] -2.06870061 -1.484205297 -0.026878250 0.006586116
## [7,] -2.43586845 -0.047485118 -0.334350297 -0.036652767
## [8,] -2.22539189 -0.222403002 0.088399352 -0.024529919
iris2 <- cbind(iris, myPr$x) # Gộp dữ liệu cũ với dữ liệu toạ độ mới
head(iris2)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species PC1
## 1 5.1 3.5 1.4 0.2 setosa -2.257141
## 2 4.9 3.0 1.4 0.2 setosa -2.074013
## 3 4.7 3.2 1.3 0.2 setosa -2.356335
## 4 4.6 3.1 1.5 0.2 setosa -2.291707
## 5 5.0 3.6 1.4 0.2 setosa -2.381863
## 6 5.4 3.9 1.7 0.4 setosa -2.068701
## PC2 PC3 PC4
## 1 -0.4784238 0.12727962 0.024087508
## 2 0.6718827 0.23382552 0.102662845
## 3 0.3407664 -0.04405390 0.028282305
## 4 0.5953999 -0.09098530 -0.065735340
## 5 -0.6446757 -0.01568565 -0.035802870
## 6 -1.4842053 -0.02687825 0.006586116
Vẽ biểu đồ PCA
library(ggplot2) # install.packages("ggplot2")
ggplot(iris2, aes(PC1, PC2, col = Species, fill = Species)) +
stat_ellipse(geom = "polygon", col = "black", alpha = 0.5) +
geom_point(shape = 21, col = "black") # Biểu diễn bằng điểm
Tuỳ biến
percentVar <- round(100 * summary(myPr)$importance[2, 1:2], 0) # Tính %
các phương sai
ggplot(iris2, aes(PC1, PC2, color = Species, shape = Species)) + # Vẽ biểu
đồ bằng ggplot2
geom_point(size = 2) +
xlab(paste0("PC1: ", percentVar[1], "% variance")) + # x
label
ylab(paste0("PC2: ", percentVar[2], "% variance")) + # y
label
ggtitle("Principal component analysis (PCA)") + # Tiêu
đề
theme(aspect.ratio = 1)
Kết quả trên là một phép chiếu của dữ liệu 4 chiều của data set iris trên
không gian 2 chiều bằng cách sử dụng hai thành phần chính đầu tiên.
Chúng ta có thể thấy rằng chỉ riêng thành phần chính đầu tiên cũng hữu ích
trong việc phân biệt ba loài.
Chúng ta có thể thấy: Nếu PC1 <-1, thì là Iris setosa. Nếu PC1> 1,5 thì là Iris
virginica. Nếu -1 <PC1 <1, thì Iris versicolor.
Tài liệu tham khảo
http://www.sthda.com/english/articles/31-principal-component-methods-in-r-
practical-guide/112-pca-principal-component-analysis-essentials/
EFA
factanal() function of the build-
in stats package
Evaluate data
Descriptive Statistics
Correlation matrix
KMO Test for the adequacy of data
Bartlett’s test for sphericity
R-functions
For an overview of related R-functions used by Radiant to conduct factor analysis see Multivariate >
Factor.
The key functions used in the pre_factor tool are cor from the stats package, eigen from base,
and cortest.bartlett and KMO from the psych package.
RUNNING FA
Cronbach’s alpha
QUEST <- data.frame(
Q1=c(1,5,2,3,4,2,3,4,3,2),
Q2=c(2,4,1,2,4,1,2,5,2,1),
Q3=c(2,5,1,3,3,2,2,4,2,2))
#install.packages("psych")
library(psych)
alpha(QUEST)
Determine Number of Factors to
Extract
ev <- eigen(cor(mydata)) # get eigenvalues
ev$values
fit <- factanal(mydata, Nfacs, rotation="varimax")
print(fit, digit=2, cutoff=0.5, sort=TRUE
Exploratory Factor Analysis
The factanal( ) function produces maximum likelihood factor analysis.
# Maximum Likelihood Factor Analysis
# entering raw data and extracting 3 factors,
# with varimax rotation
fit <- factanal(mydata, 3, rotation="varimax")
print(fit, digits=2, cutoff=.3, sort=TRUE)
# plot factor 1 by factor 2
load <- fit$loadings[,1:2]
plot(load,type="n") # set up plot
text(load,labels=names(mydata),cex=.7) # add variable names
The rotation= options include "varimax", "promax", and "none". Add the
option scores="regression" or "Bartlett" to produce factor scores. Use
the covmat= option to enter a correlation or covariance matrix directly. If entering a
covariance matrix, include the option n.obs=.
The factor.pa( ) function in the psych package offers a number of factor analysis
related functions, including principal axis factoring.
CFA
install.packages(“lavaan”)
library(lavaan)
Then we define the model by specifying the relationship between items and
factors:
path <- ‘
f1 =~ E1 + E2 + E3 + E4 + E5 + E6 + E7 + E8 + E9 + E10
f2 =~ N1 + N2 + N3 + N4 + N5 + N6 + N7 + N8 + N9 + N10 f3 =~ A1 + A2 + A3 + A4 +
A5 + A6 + A7 + A8 + A9 + A10 f4 =~ C1 + C2 + C3 + C4 + C5 + C6 + C7 + C8 + C9 +
C10 f5 =~ O1 + O2 + O3 + O4 + O5 + O6 + O7 + O8 + O9 + O10
fit the model and output the results:
model <- cfa(path, data= pcafit)
summary(model, fit.measures=TRUE)
Criteria:
Model fit:
Chi square test (CMIN), CMIN/df <2 (or 3)
the CFI and TLI cut-off scores should be above 0.9.
The value of RMSEA and SRMR is less than 0.05 (but 0.08 is acceptable
All of the estimate coefficients loadings are significant
