PCA: phân tích dữ liệu về rượu (tutle shell)

1) Mô tả dữ liệu

Nguồn (data + code R): TNV. Hồi quy. NXB TH TP.HCM, 2020.

Dữ liệu đo đạc mai rùa:

• 3 biến (variables)
• 24 quan trắc (observations) Tải thư viện:

library(magrittr)
library(dplyr)
library(GGally)

Đọc dữ liệu:

tt <- read.table("shell.txt", header = TRUE)

Xem 6 hàng đầu:

head(tt)

##   length width height
## 1     98    81     38
## 2    103    84     38
## 3    103    86     42
## 4    105    86     42
## 5    109    88     44
## 6    123    92     50

str(tt)

## 'data.frame':    24 obs. of  3 variables:
##  $ length: int  98 103 103 105 109 123 123 133 133 133 ...
##  $ width : int  81 84 86 86 88 92 95 99 102 102 ...
##  $ height: int  38 38 42 42 44 50 46 51 51 51 ...

Xem cấu trúc:

str(tt)

## 'data.frame':    24 obs. of  3 variables:
##  $ length: int  98 103 103 105 109 123 123 133 133 133 ...
##  $ width : int  81 84 86 86 88 92 95 99 102 102 ...
##  $ height: int  38 38 42 42 44 50 46 51 51 51 ...

mối tương quan:

tt %>% select(length, width, height) %>% ggpairs(upper = list(continuous = "points"), lower = list(continuous = "cor"))

names(tt) # Cho biết tên các biến trong data.frame  

## [1] "length" "width"  "height"

Trong 14 biến các biến trên thấy nhiều biến tương quan

2) Thống kê mô tả

Tính trung bình mẫu của các biến:

apply(tt, 2, mean) 

##    length     width    height 
## 136.04167 102.58333  52.04167

Tính phương sai mẫu của các biến:

apply(tt, 2, var)

##    length     width    height 
## 451.51993 171.73188  64.73732

Dùng đồ thị boxplot:

boxplot(tt)

apply(tt, 2, range)

##      length width height
## [1,]     98    81     38
## [2,]    177   132     67

Nhận xét: từ đồ thị boxplot và kết quả tính phạm vi biến thiên của dữ liệu (range), ta thấy rằng đơn vị và phương sai của các biến rất khác biệt nhau. Do đó ta cần thực hiện chuẩn hóa dữ liệu (standardizing data) trước khi xử lý. Ta có thể thực hiện quy tâm (mean-centering) hoặc/và co giãn (scaling) dữ liệu. Các đồ thị boxplot dưới đây sẽ minh họa các bước xử lý trên.

Quy tâm dữ liệu:

tt.n <- scale(tt, scale = FALSE) # Chỉ thực hiện quy tâm, không co giãn 
tt.n <- as.data.frame(tt.n)
boxplot(tt)

Quy tâm và co giãn dữ liệu (chuẩn hóa):

tt.scale <- as.data.frame(scale(tt)) 
boxplot(tt.scale)

Tính ma trận hiệp phương sai (covariance matrix) và ma trận hệ số tương quan (correlation matrix):

tt.cov <- cov(tt) # Ma trận hiệp phương sai tính từ dữ liệu chưa chuẩn hóa 
tt.cov

##          length    width    height
## length 451.5199 270.9746 165.95471
## width  270.9746 171.7319 101.84420
## height 165.9547 101.8442  64.73732

tt.cor <- cor(tt) # Ma trận hệ số tương quan tính từ dữ liệu chưa chuẩn hóa 
tt.cor

##           length     width    height
## length 1.0000000 0.9731162 0.9706748
## width  0.9731162 1.0000000 0.9659029
## height 0.9706748 0.9659029 1.0000000

Vẽ ma trận đồ thị phân tán:

pairs(tt)

Ta có thể nhận thấy rằng các biến có tương quan tuyến tính mạnh

plot(tt$height , tt$width, pch = 16, col = "blue")

3) Phân tích thành phần chính (PCA):

Ta thực hiện PCA từ bộ dữ liệu đã được chuẩn hóa tt.scale.

pca_tt.scale <- princomp(tt.scale)
summary(pca_tt.scale)

## Importance of components:
##                          Comp.1    Comp.2      Comp.3
## Standard deviation     1.678484 0.1814027 0.157434408
## Proportion of Variance 0.979933 0.0114459 0.008621076
## Cumulative Proportion  0.979933 0.9913789 1.000000000

Trình bày các PC loadings:

pca_tt.scale$loadings

## 
## Loadings:
##        Comp.1 Comp.2 Comp.3
## length  0.578  0.137  0.804
## width   0.577  0.628 -0.522
## height  0.577 -0.766 -0.284
## 
##                Comp.1 Comp.2 Comp.3
## SS loadings     1.000  1.000  1.000
## Proportion Var  0.333  0.333  0.333
## Cumulative Var  0.333  0.667  1.000

Trình bày 3 hàng đầu các PC scores:

pca_tt.scale$scores[1:3,]

##         Comp.1      Comp.2      Comp.3
## [1,] -2.992145  0.05696787 -0.08482498
## [2,] -2.723970  0.23302866 -0.01510844
## [3,] -2.349171 -0.05201949 -0.23582973

Xác định số thành phần chính cần giữ lại: vẽ đồ thị scree plot

screeplot(pca_tt.scale, type = "lines", main = "Scree plot")

Nhận xét: từ đô thị scree plot, ta thấy rằng có thể giữ lại 1 thành phần chính đầu (giải thích được \(98.%\) phương sai toàn phần).

Tiêu chuẩn khác: vì PCA được thực hiện trên dữ liệu chuẩn hóa, ta có thể dùng tiêu chuẩn giữ lại các thành phần chính mà có trị riêng (phương sai) lớn hơn 1. Tính các trị riêng:

eigenvalues_pca <- pca_tt.scale$sdev^2 
eigenvalues_pca

##     Comp.1     Comp.2     Comp.3 
## 2.81730746 0.03290695 0.02478559

Vẽ các PC scores và loadings:

tt_scores <- pca_tt.scale$scores

# Plot the scores
plot(tt_scores, xlab = "PC 1 (36%)", ylab = "PC 2 (16.2%)", main = "Scores plot")
abline(h = 0, v = 0, col = "gray")

# Plot the scores as points
text(tt_scores[, 1] + 0.2, tt_scores[, 2], label = rownames(tt_scores), col="blue", cex=0.6)

4) Phân tích bằng faxtoextra

• Tính toán

Kiểm tra ggpair

tt %>% select(length, width, height) %>% ggpairs(upper = list(continuous = "points"), lower = list(continuous = "cor"))

#Phân tích PCA
pca <- prcomp(tt, scale = TRUE)
print(pca)

## Standard deviations (1, .., p=3):
## [1] 1.7145842 0.1853043 0.1608205
## 
## Rotation (n x k) = (3 x 3):
##              PC1        PC2        PC3
## length 0.5781417 -0.1373949 -0.8042853
## width  0.5771970 -0.6278484  0.5221591
## height 0.5767112  0.7661129  0.2836814

Kiểm tra bằng lệnh summary

summary(pca)

## Importance of components:
##                           PC1     PC2     PC3
## Standard deviation     1.7146 0.18530 0.16082
## Proportion of Variance 0.9799 0.01145 0.00862
## Cumulative Proportion  0.9799 0.99138 1.00000

Có thể thấy rằng PC1 giải thích 97.99 % phương sai: Vẽẽ biểu đồ bằng hàm fviz_eig() trong thư viện factoextra

library(factoextra)
fviz_eig(pca)

#Như vậy PC1 giải thích gần 100% phương sai
fviz_pca_ind(pca, col.ind = "cos2", gradient.cols = c("red","blue", "green1"), repel = TRUE)

#Như vậy component 1 giải thích 98%, component 2 giải thích 1.1 %. 

#biểu đồ 3:
fviz_pca_var(pca, col.var = "contrib", gradient.cols = c("red", "blue", "black"), repel = TRUE)

#Phần trăm của contribution của mỗi một component

#Tóm lại PCA giúp cho 3 biến giảm 1 biến, có thể giải thích được % phương sai, giúp cho phát hiện ngoại vi hiệu quả.


######################