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

1) Mô tả dữ liệu

Nguồn (data + code R): R. Wehrens. Chemometrics with R: Multivariate Data Analysis in the Natural Science and Life Sciences. Heidelberg, Germany: Springer, 2011.

Dữ liệu về 3 loại rượu (Barbera, Barolo và Grignolino) ở vùng Piedmont, Ý, gồm:

• 14 biến (variables)
• 177 quan trắc (observations)

Rượu Barolo được làm từ nho Nebbiolo, hai loại rượu còn lại được làm từ loại nho có cùng tên với chúng.

Đọc dữ liệu:

wines <- read.csv("wines.csv", header = T, sep = ";") # Đọc và gán dữ liệu từ file wines.csv vào data.frame wines 
head(wines) # Trình bày 10 hàng đầu của data.frame wines

##   alcohol malic.acid  ash ash.alkalinity magnesium tot..phenols flavonoids
## 1   13.20       1.78 2.14           11.2       100         2.65       2.76
## 2   13.16       2.36 2.67           18.6       101         2.80       3.24
## 3   14.37       1.95 2.50           16.8       113         3.85       3.49
## 4   13.24       2.59 2.87           21.0       118         2.80       2.69
## 5   14.20       1.76 2.45           15.2       112         3.27       3.39
## 6   14.39       1.87 2.45           14.6        96         2.50       2.52
##   non.flav..phenols proanth col..int. col..hue OD.ratio proline vintages
## 1              0.26    1.28      4.38     1.05     3.40    1050   Barolo
## 2              0.30    2.81      5.68     1.03     3.17    1185   Barolo
## 3              0.24    2.18      7.80     0.86     3.45    1480   Barolo
## 4              0.39    1.82      4.32     1.04     2.93     735   Barolo
## 5              0.34    1.97      6.75     1.05     2.85    1450   Barolo
## 6              0.30    1.98      5.25     1.02     3.58    1290   Barolo

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

##  [1] "alcohol"           "malic.acid"        "ash"              
##  [4] "ash.alkalinity"    "magnesium"         "tot..phenols"     
##  [7] "flavonoids"        "non.flav..phenols" "proanth"          
## [10] "col..int."         "col..hue"          "OD.ratio"         
## [13] "proline"           "vintages"

Trong 14 biến, biến \(vintages\) là biến nhãn cho biết tên loại rượu của từng quan trắc tương ứng. Ngoại trừ các biến \(col..int\), \(col..hue\), \(OD.ratio\), các biến còn lại đều là các biến mô tả nồng độ (concentration). \(col..int\) = color intensity, \(col..hue\) = color hue và \(OD.ration\) = tỷ số giữa sự hấp thụ (absorbance) tại các bước sóng \(280\) và \(315\) nm.

Thống kê số lượng 3 loại rượu trong mẫu:

table(wines$vintages)

## 
##    Barbera     Barolo Grignolino 
##         48         58         71

2) Thống kê mô tả

Trước tiên, ta tạo một data.frame mới gồm 13 cột đầu tiên của data.frame wine (gồm các biến định lượng) để thực hiện thống kê mô tả (biến thứ 14 là vintages là biến định tính):

dat_wines <- wines[ , 1:13]

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

apply(dat_wines, 2, mean) 

##           alcohol        malic.acid               ash    ash.alkalinity 
##        12.9936723         2.3398870         2.3661582        19.5169492 
##         magnesium      tot..phenols        flavonoids non.flav..phenols 
##        99.5875706         2.2922599         2.0234463         0.3623164 
##           proanth         col..int.          col..hue          OD.ratio 
##         1.5869492         5.0548023         0.9569831         2.6042938 
##           proline 
##       745.0960452

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

apply(dat_wines, 2, var)

##           alcohol        malic.acid               ash    ash.alkalinity 
##      6.541711e-01      1.252865e+00      7.566925e-02      1.112937e+01 
##         magnesium      tot..phenols        flavonoids non.flav..phenols 
##      2.009028e+02      3.924585e-01      9.973170e-01      1.553835e-02 
##           proanth         col..int.          col..hue          OD.ratio 
##      3.266634e-01      5.403051e+00      5.250287e-02      4.971701e-01 
##           proline 
##      9.915196e+04

Dùng đồ thị boxplot:

boxplot(dat_wines)

apply(dat_wines, 2, range)

##      alcohol malic.acid  ash ash.alkalinity magnesium tot..phenols flavonoids
## [1,]   11.03       0.74 1.36           10.6        70         0.98       0.34
## [2,]   14.83       5.80 3.23           30.0       162         3.88       5.08
##      non.flav..phenols proanth col..int. col..hue OD.ratio proline
## [1,]              0.13    0.41      1.28     0.48     1.27     278
## [2,]              0.66    3.58     13.00     1.71     4.00    1680

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:

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

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

sc_wines <- as.data.frame(scale(dat_wines)) 
boxplot(sc_wines)

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

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

##                        alcohol   malic.acid           ash ash.alkalinity
## alcohol             0.65417110   0.09049758  0.0469369158     -0.8185115
## malic.acid          0.09049758   1.25286476  0.0507899043      1.0685076
## ash                 0.04693692   0.05078990  0.0756692476      0.4099291
## ash.alkalinity     -0.81851146   1.06850761  0.4099291217     11.1293702
## magnesium           2.96623909  -0.77817187  1.1194292501     -3.3906972
## tot..phenols        0.14417518  -0.23386224  0.0220882768     -0.6637260
## flavonoids          0.18588386  -0.45754734  0.0313400199     -1.1558031
## non.flav..phenols  -0.01526878   0.04067185  0.0064242906      0.1494548
## proanth             0.05896752  -0.13944694  0.0012705990     -0.3637605
## col..int.           1.03003738   0.65058464  0.1653787819      0.1587988
## col..hue           -0.01396897  -0.14384455 -0.0047387018     -0.2084690
## OD.ratio            0.03274437  -0.28942678  0.0002915896     -0.6308459
## proline           163.26765665 -66.79419363 19.3141210837   -458.9084553
##                      magnesium tot..phenols   flavonoids non.flav..phenols
## alcohol              2.9662391   0.14417518   0.18588386      -0.015268782
## malic.acid          -0.7781719  -0.23386224  -0.45754734       0.040671854
## ash                  1.1194293   0.02208828   0.03134002       0.006424291
## ash.alkalinity      -3.3906972  -0.66372602  -1.15580306       0.149454834
## magnesium          200.9027992   1.84872143   2.64841808      -0.445402863
## tot..phenols         1.8487214   0.39245850   0.54056774      -0.035008105
## flavonoids           2.6484181   0.54056774   0.99731703      -0.066764847
## non.flav..phenols   -0.4454029  -0.03500811  -0.06676485       0.015538354
## proanth              1.8349278   0.21860296   0.37115035      -0.025880961
## col..int.            6.5675033  -0.08213080  -0.40486379       0.040620630
## col..hue             0.1690214   0.06215322   0.12430080      -0.007475017
## OD.ratio             0.4693378   0.30901411   0.55372887      -0.044110003
## proline           1729.6648369  97.81114535 154.45722393     -12.124144196
##                        proanth    col..int.     col..hue      OD.ratio
## alcohol            0.058967517   1.03003738 -0.013968971  0.0327443695
## malic.acid        -0.139446938   0.65058464 -0.143844547 -0.2894267848
## ash                0.001270599   0.16537878 -0.004738702  0.0002915896
## ash.alkalinity    -0.363760497   0.15879883 -0.208469029 -0.6308459168
## magnesium          1.834927773   6.56750329  0.169021379  0.4693377632
## tot..phenols       0.218602956  -0.08213080  0.062153220  0.3090141050
## flavonoids         0.371150347  -0.40486379  0.124300797  0.5537288681
## non.flav..phenols -0.025880961   0.04062063 -0.007475017 -0.0441100026
## proanth            0.326663367  -0.03601878  0.038554493  0.2069052196
## col..int.         -0.036018783   5.40305118 -0.278351336 -0.7141730042
## col..hue           0.038554493  -0.27835134  0.052502869  0.0916705277
## OD.ratio           0.206905220  -0.71417300  0.091670528  0.4971700950
## proline           58.621999230 231.02095816 16.946768683 67.9468011685
##                       proline
## alcohol             163.26766
## malic.acid          -66.79419
## ash                  19.31412
## ash.alkalinity     -458.90846
## magnesium          1729.66484
## tot..phenols         97.81115
## flavonoids          154.45722
## non.flav..phenols   -12.12414
## proanth              58.62200
## col..int.           231.02096
## col..hue             16.94677
## OD.ratio             67.94680
## proline           99151.96231

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

##                       alcohol  malic.acid          ash ash.alkalinity
## alcohol            1.00000000  0.09996298  0.210964395    -0.30334986
## malic.acid         0.09996298  1.00000000  0.164955040     0.28614768
## ash                0.21096440  0.16495504  1.000000000     0.44669776
## ash.alkalinity    -0.30334986  0.28614768  0.446697755     1.00000000
## magnesium          0.25874233 -0.04904903  0.287107111    -0.07170686
## tot..phenols       0.28454303 -0.33351175  0.128175570    -0.31758259
## flavonoids         0.23013326 -0.40932410  0.114083528    -0.34692207
## non.flav..phenols -0.15144545  0.29150054  0.187353998     0.35939511
## proanth            0.12756072 -0.21797499  0.008081623    -0.19077876
## col..int.          0.54788293  0.25005306  0.258642888     0.02047823
## col..hue          -0.07537498 -0.56085396 -0.075181010    -0.27271858
## OD.ratio           0.05741673 -0.36671960  0.001503349    -0.26818563
## proline            0.64106760 -0.18951167  0.222979319    -0.43685781
##                     magnesium tot..phenols flavonoids non.flav..phenols
## alcohol            0.25874233   0.28454303  0.2301333        -0.1514454
## malic.acid        -0.04904903  -0.33351175 -0.4093241         0.2915005
## ash                0.28710711   0.12817557  0.1140835         0.1873540
## ash.alkalinity    -0.07170686  -0.31758259 -0.3469221         0.3593951
## magnesium          1.00000000   0.20820043  0.1871014        -0.2520911
## tot..phenols       0.20820043   1.00000000  0.8640455        -0.4483005
## flavonoids         0.18710135   0.86404554  1.0000000        -0.5363259
## non.flav..phenols -0.25209110  -0.44830051 -0.5363259         1.0000000
## proanth            0.22650394   0.61053272  0.6502540        -0.3632684
## col..int.          0.19933693  -0.05640137 -0.1744106         0.1401924
## col..hue           0.05204238   0.43298739  0.5432075        -0.2617087
## OD.ratio           0.04696129   0.69956639  0.7863720        -0.5018594
## proline            0.38754158   0.49583915  0.4911803        -0.3088858
##                        proanth   col..int.    col..hue     OD.ratio    proline
## alcohol            0.127560717  0.54788293 -0.07537498  0.057416730  0.6410676
## malic.acid        -0.217974990  0.25005306 -0.56085396 -0.366719602 -0.1895117
## ash                0.008081623  0.25864289 -0.07518101  0.001503349  0.2229793
## ash.alkalinity    -0.190778761  0.02047823 -0.27271858 -0.268185630 -0.4368578
## magnesium          0.226503941  0.19933693  0.05204238  0.046961289  0.3875416
## tot..phenols       0.610532721 -0.05640137  0.43298739  0.699566388  0.4958392
## flavonoids         0.650254005 -0.17441056  0.54320753  0.786372013  0.4911803
## non.flav..phenols -0.363268448  0.14019245 -0.26170871 -0.501859420 -0.3088858
## proanth            1.000000000 -0.02711186  0.29439692  0.513415210  0.3257315
## col..int.         -0.027111858  1.00000000 -0.52261546 -0.435743974  0.3156321
## col..hue           0.294396923 -0.52261546  1.00000000  0.567395275  0.2348793
## OD.ratio           0.513415210 -0.43574397  0.56739527  1.000000000  0.3060313
## proline            0.325731494  0.31563211  0.23487929  0.306031309  1.0000000

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

pairs(sc_wines)

Ta có thể nhận thấy rằng hai biến \(tot..phenols\) và \(flavonoids\) có tương quan tuyến tính mạnh (hệ số tương quan bằng \(0.864\)):

plot(sc_wines$tot..phenols, sc_wines$flavonoids, pch = 16, col = "blue")

trong khi biến \(alcohol\) và \(col..hue\) hầu như không có mối tương quan tuyến tính (hệ số tương quan bằng \(-0.0754\)):

plot(sc_wines$alcohol, sc_wines$col..hue, 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 sc_wines.

pca_wines <- princomp(sc_wines)
summary(pca_wines)

## Importance of components:
##                           Comp.1    Comp.2    Comp.3     Comp.4     Comp.5
## Standard deviation     2.1567037 1.5770968 1.2021310 0.95876028 0.92567174
## Proportion of Variance 0.3598307 0.1924128 0.1117946 0.07111109 0.06628744
## Cumulative Proportion  0.3598307 0.5522435 0.6640381 0.73514919 0.80143663
##                            Comp.6     Comp.7     Comp.8     Comp.9    Comp.10
## Standard deviation     0.80075246 0.74085307 0.59055673 0.53623338 0.49539305
## Proportion of Variance 0.04960367 0.04246014 0.02697991 0.02224462 0.01898528
## Cumulative Proportion  0.85104030 0.89350044 0.92048035 0.94272497 0.96171025
##                           Comp.11    Comp.12     Comp.13
## Standard deviation     0.47346227 0.40917666 0.321500291
## Proportion of Variance 0.01734155 0.01295206 0.007996133
## Cumulative Proportion  0.97905180 0.99200387 1.000000000

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

pca_wines$loadings

## 
## Loadings:
##                   Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## alcohol            0.138  0.486  0.209         0.270  0.211         0.401
## malic.acid        -0.246  0.222        -0.533         0.531 -0.434       
## ash                       0.315 -0.624  0.205  0.160  0.155  0.145 -0.171
## ash.alkalinity    -0.237        -0.614                       0.290  0.428
## magnesium          0.135  0.300 -0.136  0.392 -0.706        -0.318 -0.152
## tot..phenols       0.396        -0.145 -0.203  0.133               -0.407
## flavonoids         0.424        -0.149 -0.156                      -0.189
## non.flav..phenols -0.299        -0.169  0.175  0.510 -0.273 -0.588 -0.231
## proanth            0.313        -0.151 -0.392 -0.162 -0.540 -0.362  0.368
## col..int.                 0.528  0.136               -0.414  0.234       
## col..hue           0.300 -0.274         0.419  0.195        -0.237  0.434
## OD.ratio           0.377 -0.165 -0.167 -0.190         0.267              
## proline            0.284  0.370  0.128  0.224  0.167  0.117         0.116
##                   Comp.9 Comp.10 Comp.11 Comp.12 Comp.13
## alcohol            0.488  0.220   0.272   0.242         
## malic.acid               -0.301  -0.119  -0.111         
## ash               -0.315 -0.114   0.485          -0.142 
## ash.alkalinity     0.200  0.123  -0.464                 
## magnesium          0.270  0.105                         
## tot..phenols       0.321 -0.247  -0.324   0.318  -0.464 
## flavonoids               -0.159                   0.832 
## non.flav..phenols  0.185  0.245                   0.114 
## proanth           -0.220          0.250          -0.117 
## col..int.                -0.296          -0.596         
## col..hue           0.118 -0.524          -0.248         
## OD.ratio                  0.517          -0.614  -0.157 
## proline           -0.575  0.194  -0.525                 
## 
##                Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9
## SS loadings     1.000  1.000  1.000  1.000  1.000  1.000  1.000  1.000  1.000
## Proportion Var  0.077  0.077  0.077  0.077  0.077  0.077  0.077  0.077  0.077
## Cumulative Var  0.077  0.154  0.231  0.308  0.385  0.462  0.538  0.615  0.692
##                Comp.10 Comp.11 Comp.12 Comp.13
## SS loadings      1.000   1.000   1.000   1.000
## Proportion Var   0.077   0.077   0.077   0.077
## Cumulative Var   0.769   0.846   0.923   1.000

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

pca_wines$scores[1:10,]

##         Comp.1     Comp.2     Comp.3       Comp.4      Comp.5      Comp.6
##  [1,] 2.223934 -0.3014576  2.0271695  0.281108579  0.25880549  0.92499071
##  [2,] 2.524760  1.0592518 -0.9739613 -0.733645703  0.19804048 -0.55567525
##  [3,] 3.744056  2.7973729  0.1798599 -0.575492236  0.25714173 -0.09982545
##  [4,] 1.017245  0.8858673 -2.0181445  0.431568191 -0.27445613  0.40199816
##  [5,] 3.040574  2.1638681  0.6369402  0.486248271  0.62957197 -0.13044651
##  [6,] 2.451274  1.2036500  0.9854402  0.004664938  1.02718855  0.61172791
##  [7,] 2.055773  1.6358443 -0.1433608  1.196313161 -0.01045482  1.44104903
##  [8,] 2.511320  0.9581190  1.7773376 -0.104420449  0.87123084  0.12082683
##  [9,] 2.760141  0.8221890  0.9861588 -0.373844189  0.43657501 -0.14383730
## [10,] 3.479291  1.3513568  0.4281045 -0.039867610  0.31598814  0.17817536
##            Comp.7     Comp.8      Comp.9     Comp.10      Comp.11      Comp.12
##  [1,] -0.07949880 -1.0235663 -0.31225947  0.13088481  0.152818951 -0.399900961
##  [2,] -0.43112665  0.3346618 -1.17573332  0.00673330  0.274595491 -0.003370376
##  [3,]  0.36389217 -0.6450175  0.06771425  0.37226378 -0.694465299  0.240417778
##  [4,] -0.45343437 -0.4108678  0.33710279 -0.09604793  0.539928345  0.187265337
##  [5,] -0.42010486 -0.3976031 -0.11313320 -0.01996873 -0.388016985  0.379779318
##  [6,] -0.06595238  0.3742163 -0.53262110  0.92254164  0.557313940 -0.188210141
##  [7,] -0.05822750 -0.2276222  0.08138933  0.79090393 -0.187584784 -0.410169540
##  [8,] -0.13395273  0.5056063  0.59911659  0.17695055  0.570545139  0.591536770
##  [9,]  0.86746006 -0.1525266  0.21302021  0.19066911  0.003143058 -0.564267439
## [10,] -0.27231968  1.1947708 -0.48455491 -0.13358646 -0.751897200  0.101901305
##            Comp.13
##  [1,]  0.001895375
##  [2,]  0.021545005
##  [3,] -0.369417990
##  [4,] -0.081588741
##  [5,]  0.144170806
##  [6,] -0.273844055
##  [7,] -0.113143937
##  [8,]  0.139422584
##  [9,] -0.044249196
## [10,]  0.123080932

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

screeplot(pca_wines, type = "lines", main = "Scree plot")

Nhận xét: từ đô thị scree plot, ta thấy rằng có thể giữ lại 4 thành phần chính đầu (giải thích được \(73.51%\) phương sai toàn phần) hoặc 5 thành phần chính đầu (giải thích được \(80.14%\) 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_wines <- pca_wines$sdev^2 
eigenvalues_wines

##    Comp.1    Comp.2    Comp.3    Comp.4    Comp.5    Comp.6    Comp.7    Comp.8 
## 4.6513710 2.4872343 1.4451189 0.9192213 0.8568682 0.6412045 0.5488633 0.3487573 
##    Comp.9   Comp.10   Comp.11   Comp.12   Comp.13 
## 0.2875462 0.2454143 0.2241665 0.1674255 0.1033624

Vẽ các PC scores và loadings:

wines_scores <- pca_wines$scores

# Plot the scores
plot(wines_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(wines_scores[, 1] + 0.2, wines_scores[, 2], label = rownames(sc_wines), col="blue", cex=0.6)

wines_loadings <- pca_wines$loadings

# Plot the loading vector
plot(wines_loadings, xlab = "PC 1 (36%)", ylab = "PC 2 (16.2%)", main = "Loading plot", xlim = c(-0.5,0.8), ylim = c(-0.5,0.8))
abline(h = 0, v = 0, col = "gray")

# Plot the loadings as points
text(wines_loadings[, 1] + 0.1, wines_loadings[, 2] + 0.1, label = rownames(wines_loadings), col="blue", cex = 0.8)

Vẽ đồ thị biplot:

biplot(pca_wines, cex = c(0.6, 1), xlim = c(-0.2, 0.2), ylim = c(-0.2, 0.2), pch = ".",
       xlab = "PC 1 (36%)", ylab = "PC 2 (16.2%)", main = "Biplot/n/n")
abline(v = 0, h = 0, col = 'gray')

Vẽ biplot sử dụng hàm autoplot:

library(ggfortify)

autoplot(pca_wines, loadings = TRUE, loadings.colour = 'blue', loadings.label = TRUE, loadings.label.size = 3)

4) Xác định cụm (clusters) bằng PCA:

• Để xác định chùm (clusters) để phân loại các loại rượu, ta vẽ các scores tương ứng với hai thành phần chính đầu tiên.

autoplot(pca_wines, loadings = TRUE, loadings.colour = 'blue', loadings.label = TRUE, loadings.label.size = 3,
         colour = "vintages", data = wines)

5) Xác định các điểm outliers dùng PCA

Để xác định các điểm outliers, ta vẽ các scores tương ứng với các thành phần chính cuối.

# Plot the scores
plot(wines_scores[,12:13], xlab = "PC 12", ylab = "PC 13", main = "Scores plot")
abline(h = 0, v = 0, col = "gray")

# Plot the scores as points
text(wines_scores[, 12] + 0.07, wines_scores[, 13], label = rownames(sc_wines), col="blue", cex = 0.6)