PCA funkcije u R-u

U R-u postoje funkcije u standardnom paketu za sprovođenje PCA: prcomp() i princomp(). Koristićemo i paket FactoMineR u kome se nalazi funkcija PCA()

irisClean <- na.omit(iris[,-5])
pca1 <- prcomp(irisClean,scale=FALSE)
pca2 <- princomp(irisClean,cor=FALSE, scores=TRUE)
pca3 <- PCA(irisClean, scale.unit = TRUE, ncp = 5, graph = TRUE)

Argumenti za prcomp()

Argumenti za princomp()

Argumenti za PCA()

Na izlazu dobijamo sledeće vrednosti

Jedan paket za vizuelizaciju PCA je factoextra - grafici nam dolaze iz ggplot2 paketa.

library(factoextra)

Pokazaćemo na primeru kako se koristi ovaj paket. Posmatrajmo dataframe decathlon2 iz ovog paketa. U njoj je dato 27 opservacija i 13 varijabli. Na slici vidimo deo ovog dataframe-a.

Primenićemo PCA na prvih 23 opservacija (svetloplave vrste) i prvih 10 varijabli (roze kolone)

data(decathlon2)
decathlon2.active <- decathlon2[1:23, 1:10]
head(decathlon2.active[, 1:6])
          X100m Long.jump Shot.put High.jump X400m X110m.hurdle
SEBRLE    11.04      7.58    14.83      2.07 49.81        14.69
CLAY      10.76      7.40    14.26      1.86 49.37        14.05
BERNARD   11.02      7.23    14.25      1.92 48.93        14.99
YURKOV    11.34      7.09    15.19      2.10 50.42        15.31
ZSIVOCZKY 11.13      7.30    13.48      2.01 48.62        14.17
McMULLEN  10.83      7.31    13.76      2.13 49.91        14.38

Korelaciona analiza

Koristićemo paket corrr

Prikaz korelacija.

library(magrittr)
library(corrr)

korObjekat <- decathlon2.active %>%
  correlate() %>% # Create correlation data frame (cor_df)
  rearrange() %>% # rearrange by correlations
  shave() # Shave off the upper triangle for a clean result

library(kableExtra)

korMatrica <- fashion(korObjekat)
kable(korMatrica,caption = 'Korelaciona matrica') %>% 
  kableExtra::kable_styling(bootstrap_options = c("striped", "hover","bordered"),full_width = F,position = "center")
Korelaciona matrica
rowname X100m X110m.hurdle X400m Pole.vault X1500m Javeline High.jump Shot.put Discus Long.jump
X100m
X110m.hurdle .73
X400m .59 .58
Pole.vault .09 .14 .25
X1500m -.22 -.05 .29 .39
Javeline -.29 -.07 -.13 .23 .09
High.jump -.40 -.25 -.37 -.50 -.26 .22
Shot.put -.45 -.38 -.31 .02 -.05 .48 .53
Discus -.48 -.53 -.36 -.19 .08 .28 .34 .71
Long.jump -.76 -.59 -.51 .02 .22 .37 .34 .44 .46 
rplot(korObjekat)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

decathlon2.active %>%
  correlate() %>%
  network_plot(min_cor = .2)

Provera adekvatnosti podataka za PCA i određivanje broja komponenti

U zavisnosti od vrednosti Kaiser-Meyer-Olkin testa kvalitet podataka za PCA se ocenjuje na sledeći način:

library(psych)
KMO(decathlon2.active)
Kaiser-Meyer-Olkin factor adequacy
Call: KMO(r = decathlon2.active)
Overall MSA =  0.68
MSA for each item = 
       X100m    Long.jump     Shot.put    High.jump        X400m X110m.hurdle 
        0.77         0.88         0.60         0.59         0.78         0.81 
      Discus   Pole.vault     Javeline       X1500m 
        0.64         0.38         0.76         0.40

Bartletov test sferičnosti je jedan od načina da se analizira korelaciona matrica. Ovaj test analizira postojanje statističke značajnosti odnosno da li postoji korelacija barem između nekih varijabli. Treba imati u vidu da sa povećanjem uzorka Bartlett-ov test postaje sve osetljiviji na otkrivanje korelacije.Da bi PCA imala smisla Bartletov test sferičnosti bi trebao da se pokaže statistički značajan za p < .05 (ili neko drugo odabrano p).

barttest <- function(dat) { #dat is your raw data
  R<-cor(dat)
  p<-ncol(dat)
  n<-nrow(dat)
  chi2<- -((n-1)-((2*p)+5)/6 ) * log(det(R)) #this is the formula
  df<-(p*(p-1)/2)
  crit<-qchisq(.95,df) #critical value
  p<-pchisq(chi2,df,lower.tail=F) #pvalue
  cat("Bartlett's test of sphericity: Hikvadrat(",
      df,")=",chi2,", p=",
      round(p,3),sep="" )   
}

barttest(decathlon2.active)
Bartlett's test of sphericity: Hikvadrat(45)=100.3876, p=0
# provera da li se rezultat naše funkcije 
# slaže sa rezultatom ugrađene funkcije u paketu psych
cortest.bartlett(decathlon2.active)
$chisq
[1] 100.3876

$p.value
[1] 4.176075e-06

$df
[1] 45

Hornova paralelna analiza je postupak koji počiva na pretpostavci da treba zadržati samo one komponente čije su sopstvene vrednosti veće od sopstvenih vrednosti koje je moguće dobiti na osnovu slučajnih podataka sa analognim karakteristikama (npr. isti broj varijabli i opservacija).

library(paran)

horn_pa <- paran(decathlon2.active, graph = TRUE)

Using eigendecomposition of correlation matrix.
Computing: 10%  20%  30%  40%  50%  60%  70%  80%  90%  100%


Results of Horn's Parallel Analysis for component retention
300 iterations, using the mean estimate

-------------------------------------------------- 
Component   Adjusted    Unadjusted    Estimated 
            Eigenvalue  Eigenvalue    Bias 
-------------------------------------------------- 
1           2.924747    4.124213      1.199465
2           1.074310    1.838530      0.764220
-------------------------------------------------- 

Adjusted eigenvalues > 1 indicate dimensions to retain.
(2 components retained)

# Procena koliko komponenti treba zadržati
horn_pa$Retained
[1] 2

Sledeća metoda određivanja broja komponenti zasnovana je na unakrsnoj validaciji. Detaljni opisi metoda unakrsne validacije za PCA i faktorsku analizu dati su ispod:

http://alexhwilliams.info/itsneuronalblog/2018/02/26/crossval/

http://mathfaculty.fullerton.edu/sbehseta/cv-for-pca.pdf

library(missMDA)
estimCV <- estim_ncpPCA(decathlon2.active,method.cv = "loo")

# Procena koliko komponenti treba zadržati
estimCV$ncp
# Procena koliko komponenti treba zadržati
estimCV$ncp
[1] 2

Za sledeću metodu odabira broja komponenti, scree test, potrebni su nam rezultati PCA. Pozivamo funkciju prcomp. Funkcijom print dobijamo sopstvene vrednosti i matricu sopstvenih vektora. Funkcijom summary dobijamo standardne devijacije svake komponente - sopstvene vrednosti, proporciju objašnjene disperzije i kumulativne proporcije.

res.pca <- prcomp(decathlon2.active, scale = TRUE)
print(res.pca)
Standard deviations (1, .., p=10):
 [1] 2.0308159 1.3559244 1.1131668 0.9052294 0.8375875 0.6502944 0.5500742
 [8] 0.5238988 0.3939758 0.3492435

Rotation (n x k) = (10 x 10):
                      PC1         PC2         PC3         PC4        PC5
X100m        -0.418859080  0.13230683 -0.27089959  0.03708806 -0.2321476
Long.jump     0.391064807 -0.20713320  0.17117519 -0.12746997  0.2783669
Shot.put      0.361388111 -0.06298590 -0.46497777  0.14191803 -0.2970589
High.jump     0.300413236  0.34309742 -0.29652805  0.15968342  0.4807859
X400m        -0.345478567 -0.21400770 -0.25470839  0.47592968  0.1240569
X110m.hurdle -0.376265119  0.01824645 -0.40325254 -0.01866477  0.2676975
Discus        0.365965721 -0.03662510 -0.15857927  0.43636361 -0.4873988
Pole.vault   -0.106985591 -0.59549862 -0.08449563 -0.37447391 -0.2646712
Javeline      0.210864329 -0.28475723 -0.54270782 -0.36646463  0.2361698
X1500m        0.002106782 -0.57855748  0.19715884  0.49491281  0.3142987
                      PC6         PC7         PC8         PC9        PC10
X100m         0.054398099 -0.16604375 -0.19988005 -0.76924639  0.12718339
Long.jump    -0.051865558 -0.28056361 -0.75850657 -0.13094589  0.08509665
Shot.put     -0.368739186 -0.01797323  0.04649571  0.12129309  0.62263702
High.jump    -0.437716883  0.05118848  0.16111045 -0.28463225 -0.38244596
X400m        -0.075796432  0.52012255 -0.44579641  0.20854176 -0.09784197
X110m.hurdle  0.004048005 -0.67276768 -0.01592804  0.41058421 -0.04475363
Discus        0.305315353 -0.25946615 -0.07550934  0.03391600 -0.49418361
Pole.vault   -0.503563524 -0.01889413  0.06282691 -0.06540692 -0.39288155
Javeline      0.556821016  0.24281145  0.10086127 -0.10268134 -0.01103627
X1500m        0.064663250 -0.20245828  0.37119711 -0.25950868  0.17991689
summary(res.pca)
Importance of components:
                          PC1    PC2    PC3     PC4     PC5     PC6     PC7
Standard deviation     2.0308 1.3559 1.1132 0.90523 0.83759 0.65029 0.55007
Proportion of Variance 0.4124 0.1839 0.1239 0.08194 0.07016 0.04229 0.03026
Cumulative Proportion  0.4124 0.5963 0.7202 0.80213 0.87229 0.91458 0.94483
                           PC8     PC9   PC10
Standard deviation     0.52390 0.39398 0.3492
Proportion of Variance 0.02745 0.01552 0.0122
Cumulative Proportion  0.97228 0.98780 1.0000

Scree test je grafički metod za određivanje broja komponenti za PCA koji koristi linijski grafik. Sopstvene vrednosti ili procenat objašnjene (zadržane) varijanse se ucrtavaju na grafik za prvih n glavnih komponenti opadajućim redosledom. Sa grafika treba identifikovati poslednji osetan pad u vrednostima, tj. „tačku preloma”, ili “lakat”, nakon koga linija postaje relativno ravna. Do te tačke se uključuju komponente koje će dalje biti razmatrane u PCA. Problem je što grafik može imati više laktova. Scree grafik za procenat zadržane varijanse dobijamo funkcijom fviz_eig.

Na kraju, bitno je naglasiti da ima smisla zadržati samo one komponente čije su sopstvene vrednosti >1, jer komponente sa sopstvenom vrednošću manjom od 1, objašnjavaju manje varijanse od jedne originalne promenljive.

fviz_eig(res.pca)

Rezultati PCA

Koristićemo nadalje funkciju PCA

library("FactoMineR")
res.pca <- PCA(decathlon2.active, graph = FALSE,nc=10)
print(res.pca)
**Results for the Principal Component Analysis (PCA)**
The analysis was performed on 23 individuals, described by 10 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"

Objasnićemo kako se računaju ove vrednosti. Navodimo i funkcije za grafike iz paketa factoextra.

Sada ćemo da objasnimo sve pomenute metode. Razdvojićemo ih u one koje se odnose na razmatranje sopstvenih vrednosti, varijabli i opservacija.

Sopstvene vrednosti

library("factoextra")
eig.val <- get_eigenvalue(res.pca)
eig.val
       eigenvalue variance.percent cumulative.variance.percent
Dim.1   4.1242133        41.242133                    41.24213
Dim.2   1.8385309        18.385309                    59.62744
Dim.3   1.2391403        12.391403                    72.01885
Dim.4   0.8194402         8.194402                    80.21325
Dim.5   0.7015528         7.015528                    87.22878
Dim.6   0.4228828         4.228828                    91.45760
Dim.7   0.3025817         3.025817                    94.48342
Dim.8   0.2744700         2.744700                    97.22812
Dim.9   0.1552169         1.552169                    98.78029
Dim.10  0.1219710         1.219710                   100.00000
sum(eig.val[,1])
[1] 10

U prvoj koloni date su sopstvene vrednosti. Vidimo da je suma svih sopstvenih vrednosti 10, a to je ukupna disperzija. Objašnjenje za ovo je da je trag kovarijacione matrice (suma disperzija) jednak sumi sopstvenih vektora. U drugoj koloni dat je procenat objašnjene disperzije i u trećoj procenat objašnjene disperzije kumulativno. Sopstvene vrednosti koje su veće od 1 znače da odgovarajuća glavna komponenta objašnjava više varijabilnosti nego neka od originalnih (standardizovanih) varijabli. Crtamo grafik sopstvenih vrednosti.

fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50))

Varijable

Informacije o varijablama nakon PCA možemo dobiti funkcijom get_pca_var()

var <- get_pca_var(res.pca)
var
Principal Component Analysis Results for variables
 ===================================================
  Name       Description                                    
1 "$coord"   "Coordinates for the variables"                
2 "$cor"     "Correlations between variables and dimensions"
3 "$cos2"    "Cos2 for the variables"                       
4 "$contrib" "contributions of the variables"

Dakle dobijamo

Varijable možemo na 2D grafiku predstaviti ako uzmemo samo prve dve glavne komponente, i nacrtamo koordinate varijabli u tom sistemu.

fviz_pca_var(res.pca, col.var = "black")

Pravac vektora je određen koordinatama varijabli u prve dve komponente. Dužina vektora, udaljenost od koordinatnog početka predstavlja kvalitet reprezentacije varijable na grafiku komponenti, to jest koliko je one dobro opisane sa ove dve komponente. Ta vrednost se računa kao sum vrednosti cos2 za tu varijablu i prve dve glavne komponente. Ako su ovo male vrednosti, onda odgovarajuća varijabla nije dobro predstavljena sa prve dve glavne komponente. Sa ovog grafika se mogu videti veze i između svih varijabli. Pozitivno korelisane varijable su grupisane zajedno, negativno korelisane varijable su usmerene suprotno. Navešćemo još neke grafike za varijable.

library("corrplot")
corrplot(var$cos2, is.corr=FALSE)

Na ovom grafiku vidimo koliko dobro je koja varijabla prikazana kojom glavnom komponentom. Vrednosti iz odgovarajuće matrice su samo prikazane grafički - što je veći i tamniji kružić to je veća vrednost cos2. Možemo nacrtati i barplot sume cos2 za prve dve komponente.

library("corrplot")
# Ukupni cos2 varijabli za Dim.1 and Dim.2
fviz_cos2(res.pca, choice = "var", axes = 1:2)

Kada sumiramo (po glavnim komponentama) vrednosti cos2 za fiksiranu varijablu, zbir je 1. Ovo je posledica Pitagorine teoreme. Ako je neka varijabla savršeno objašnjena sa prve dve komponente, onda je zbir cos2 za tu varijablu i prve dve komponente jednak 1, ta tačka se nalazi na krugu. Ako je nekoj varijabli potrebno više glavnih komponenti da bi je opisale, onda imamo tačke unutar kruga. Na grafiku varijabli možemo dodati argumente funkciji fviz_pca_var tako da različita boja bude dodeljena različitim varijablama u zavisnosti od njene dužine (sum cos2).

Osim toga koliko je varijabla dobro objašnjena komponentama, možemo govoriti i o tome koliki je doprinos varijabli jednoj komponenti. Doprinos se definiše preko cos2. cos2 za neku varijablu i neku komponentu predstavlja koliko je dobro varijabla opisana tom komponentom. Ako fiksiramo komponentu i gledamo vrednosti cos2 po varijablama, primećujemo veće vrednosti - odgovaraju varijablama koje se pružaju približno u pravcu tih komponenti, i manje vrednosti - odgovaraju varijablama koje nisu sličnog pravca. Vrednost cos2*100/(suma cos2 za tu komponentu) nam daje procenat koliko koja varijabla utiče na neku komponentu.

head(var$contrib, 4)
              Dim.1      Dim.2     Dim.3     Dim.4     Dim.5      Dim.6
X100m     17.544293  1.7505098  7.338659 0.1375524  5.389252  0.2959153
Long.jump 15.293168  4.2904162  2.930094 1.6248594  7.748815  0.2690036
Shot.put  13.060137  0.3967224 21.620432 2.0140727  8.824401 13.5968587
High.jump  9.024811 11.7715838  8.792888 2.5498795 23.115504 19.1596070
               Dim.7      Dim.8     Dim.9     Dim.10
X100m     2.75705260  3.9952035 59.174001  1.6175614
Long.jump 7.87159392 57.5332222  1.714683  0.7241439
Shot.put  0.03230371  0.2161851  1.471201 38.7676858
High.jump 0.26202607  2.5956579  8.101552 14.6264909

Slično kao malopre imamo još dva tipa grafika.

corrplot(var$contrib, is.corr=FALSE)

# Contributions of variables to PC1
fviz_contrib(res.pca, choice = "var", axes = 1, top = 10)

# Contributions of variables to PC2
fviz_contrib(res.pca, choice = "var", axes = 2, top = 10)

Ukupni doprinos neke varijable na dve komponente PC1 i PC2 računa se kao \(\text{contrib=[(C1*Eig1)+ (C2 * Eig2)]/(Eig1 + Eig2)]}\), gde su C1 i C2 doprinosi varijable prvoj i drugoj glavnoj komponenti, a Eig1 i Eig2 odgovarajuće sopstvene vrednosti. Crvena linija na grafiku predstavlja očekivani doprinos na dve komponente. U ovom slučaju očekivani doprinos prediktora na prve dve komponente je \(\text{ [(10* Eig1) + (10 * Eig2)]/(Eig1 + Eig2)}\).

Opservacije

Rezultate o opservacijama možemo dobiti na sledeći način. Iste vrednosti koje smo imali i za varijable sada možemo računati i za opservacije. Koordinate su malopre predstavljale koordinate varijabli u novom sistemu, a u novom koordinatnom sistemu možemo prikazati i tačke - opservacije koje su date u starom sistemu. Slično dobijamo i vrednosti cos2 koje sada daju odgovor na pitanje koliko je kojom komponentom dobro opisana koja opservacija. Male vrednosti cos2 za neku komponentu znače da u novom sistemu ta opservacija ima malu vrednost te komponente, ne značajnu. Doprinos neke opservacije glavnoj komponenti nam govori o tome koje su opservacije najvažnije za tu komponentu. Vrednost blizu jedinice znači veliki cos2 za tu opservaciju i mali za ostale.To bi značilo da se ta opservacija nalazi u pravcu posmatrane komponente.

ind <- get_pca_ind(res.pca)
ind
Principal Component Analysis Results for individuals
 ===================================================
  Name       Description                       
1 "$coord"   "Coordinates for the individuals" 
2 "$cos2"    "Cos2 for the individuals"        
3 "$contrib" "contributions of the individuals"
# Coordinates of individuals
head(ind$coord)
               Dim.1      Dim.2      Dim.3       Dim.4       Dim.5      Dim.6
SEBRLE     0.1955047  1.5890567  0.6424912  0.08389652  1.16829387  0.4743235
CLAY       0.8078795  2.4748137 -1.3873827  1.29838232 -0.82498206 -1.3335117
BERNARD   -1.3591340  1.6480950  0.2005584 -1.96409420  0.08419345  0.4096327
YURKOV    -0.8889532 -0.4426067  2.5295843  0.71290837  0.40782264 -0.1051753
ZSIVOCZKY -0.1081216 -2.0688377 -1.3342591 -0.10152796 -0.20145217 -0.9156681
McMULLEN   0.1212195 -1.0139102 -0.8625170  1.34164291  1.62151286 -0.1907660
                Dim.7       Dim.8       Dim.9      Dim.10
SEBRLE     0.21263396 -0.04443733  0.67416271 -0.03346803
CLAY       0.21770395 -0.63111292  0.06147666  0.32428824
BERNARD    0.41557211 -0.71967500 -0.17390589  0.10130825
YURKOV     0.33217594 -0.11758064  0.11198556  0.12238736
ZSIVOCZKY -0.09023977  0.20688886  0.53485969  0.35625335
McMULLEN  -0.48903363 -0.29967708  0.10799705  0.40201453
# Quality of individuals
head(ind$cos2)
                Dim.1      Dim.2       Dim.3       Dim.4        Dim.5
SEBRLE    0.007530179 0.49747323 0.081325232 0.001386688 0.2689026575
CLAY      0.048701249 0.45701660 0.143628117 0.125791741 0.0507850580
BERNARD   0.197199804 0.28996555 0.004294015 0.411819183 0.0007567259
YURKOV    0.096109800 0.02382571 0.778230322 0.061812637 0.0202279796
ZSIVOCZKY 0.001574385 0.57641944 0.239754152 0.001388216 0.0054654972
McMULLEN  0.002175437 0.15219499 0.110137872 0.266486530 0.3892621478
                Dim.6       Dim.7        Dim.8        Dim.9       Dim.10
SEBRLE    0.044324130 0.008907507 0.0003890334 0.0895406730 0.0002206741
CLAY      0.132690734 0.003536548 0.0297208402 0.0002820119 0.0078471026
BERNARD   0.017913117 0.018436341 0.0552910432 0.0032285723 0.0010956493
YURKOV    0.001345355 0.013419804 0.0016814404 0.0015252249 0.0018217256
ZSIVOCZKY 0.112917691 0.001096685 0.0057644781 0.0385270286 0.0170924251
McMULLEN  0.005387699 0.035406160 0.0132956152 0.0017267328 0.0239268142
# Contributions of individuals
head(ind$contrib)
               Dim.1      Dim.2      Dim.3       Dim.4       Dim.5      Dim.6
SEBRLE    0.04029447  5.9714533  1.4483919  0.03734589  8.45894063  2.3131374
CLAY      0.68805664 14.4839248  6.7537381  8.94458283  4.21794385 18.2829287
BERNARD   1.94740183  6.4234107  0.1411345 20.46819433  0.04393073  1.7252074
YURKOV    0.83308415  0.4632733 22.4517396  2.69663605  1.03075263  0.1137311
ZSIVOCZKY 0.01232413 10.1217143  6.2464325  0.05469230  0.25151025  8.6204180
McMULLEN  0.01549089  2.4310854  2.6102794  9.55055888 16.29493304  0.3741572
              Dim.7      Dim.8      Dim.9     Dim.10
SEBRLE    0.6496730 0.03128047 12.7309990 0.03992784
CLAY      0.6810236 6.30944930  0.1058653 3.74867531
BERNARD   2.4815461 8.20446284  0.8471525 0.36585167
YURKOV    1.5854988 0.21900244  0.3512830 0.53393521
ZSIVOCZKY 0.1170107 0.67803410  8.0133253 4.52411199
McMULLEN  3.4364262 1.42260514  0.3267058 5.76101406

Za crtanje grafika koristimo funkciju fviz_pca_ind().

fviz_pca_ind(res.pca)

Kao i na grafiku varijabli i sada možemo koristiti boje da bi se jasnije videle njihove cos2 vrednosti.

fviz_pca_ind(res.pca, col.ind = "cos2", 
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE # Avoid text overlapping (slow if many points)
             )

Za kraj možemo nacrtati na istom grafiku i komponente i opservacije. Ova slika predstavlja projekciju cele situacije na ravan određenu sa prve dve komponente.

fviz_pca_biplot(res.pca, repel = TRUE,
                col.var = "#2E9FDF", # Variables color
                col.ind = "#696969"  # Individuals color
                )

Ručno računanje vrednosti sa funkcijom prcomp()

Prvo računamo koordinate varijabli u novom sistemu. To su samo skalirani sopstveni vektori.

res.pca <- prcomp(decathlon2.active, scale = TRUE)

var_coord_func <- function(loadings, comp.sdev){
  loadings*comp.sdev
} #funkcija za skalarni proizvod, loadings su koordinate sopstvenog vektora

loadings <- res.pca$rotation #matrica sopstvenih vektora, zove se i matrica rotacije
sdev <- res.pca$sdev #standardne devijacije glavnih komponenti
var.coord <- t(apply(loadings, 1, var_coord_func, sdev)) #na vrste matrice primenjujemo funkciju
head(var.coord[, 1:4])
                    PC1         PC2        PC3         PC4
X100m        -0.8506257  0.17939806 -0.3015564  0.03357320
Long.jump     0.7941806 -0.28085695  0.1905465 -0.11538956
Shot.put      0.7339127 -0.08540412 -0.5175978  0.12846837
High.jump     0.6100840  0.46521415 -0.3300852  0.14455012
X400m        -0.7016034 -0.29017826 -0.2835329  0.43082552
X110m.hurdle -0.7641252  0.02474081 -0.4488873 -0.01689589

Sada računamo cos2 vrednosti - to su samo kvadrirane koordinate.

var.cos2 <- var.coord^2
head(var.cos2[, 1:4])
                   PC1          PC2        PC3          PC4
X100m        0.7235641 0.0321836641 0.09093628 0.0011271597
Long.jump    0.6307229 0.0788806285 0.03630798 0.0133147506
Shot.put     0.5386279 0.0072938636 0.26790749 0.0165041211
High.jump    0.3722025 0.2164242070 0.10895622 0.0208947375
X400m        0.4922473 0.0842034209 0.08039091 0.1856106269
X110m.hurdle 0.5838873 0.0006121077 0.20149984 0.0002854712

Za doprinose prvo koristimo funkciju apply kojom samo sumiramo po kolonama vrednosti - tako dobijemo za jednu komponentu zbir svih cos2 vrednosti. Da smo sabirali po vrstama dobili bismo samo jedinice.

comp.cos2 <- apply(var.cos2, 2, sum) #
contrib <- function(var.cos2, comp.cos2){var.cos2*100/comp.cos2}
var.contrib <- t(apply(var.cos2,1, contrib, comp.cos2))
head(var.contrib[, 1:4])
                   PC1        PC2       PC3         PC4
X100m        17.544293  1.7505098  7.338659  0.13755240
Long.jump    15.293168  4.2904162  2.930094  1.62485936
Shot.put     13.060137  0.3967224 21.620432  2.01407269
High.jump     9.024811 11.7715838  8.792888  2.54987951
X400m        11.935544  4.5799296  6.487636 22.65090599
X110m.hurdle 14.157544  0.0332933 16.261261  0.03483735

Sada prelazimo na računanje za opservacije. Koordinate možemo odmah da dobijemo funkcijom, ali računali bismo ih isto kao što smo to radili za varijable - korelacija sa sopstvenim vektorima.

ind.coord <- res.pca$x
head(ind.coord[, 1:4])
                 PC1        PC2        PC3         PC4
SEBRLE     0.1912074 -1.5541282 -0.6283688  0.08205241
CLAY       0.7901217 -2.4204156  1.3568870  1.26984296
BERNARD   -1.3292592 -1.6118687 -0.1961500 -1.92092203
YURKOV    -0.8694134  0.4328779 -2.4739822  0.69723814
ZSIVOCZKY -0.1057450  2.0233632  1.3049312 -0.09929630
McMULLEN   0.1185550  0.9916237  0.8435582  1.31215266

cos2 vrednosti za opservacije

# Cos2 of individuals
#:::::::::::::::::::::::::::::::::
# 1. square of the distance between an individual and the
# PCA center of gravity
center <- res.pca$center
scale<- res.pca$scale #koje su vrednosti kojima smo centrirai i skalirali podatke
getdistance <- function(ind_row, center, scale){
  return(sum(((ind_row-center)/scale)^2))
  } #ova funkcija je namenjena da za jednu opservaciju izračuna kvadrat rastojanja od centra podataka, ako skaliramo ovo je vec 1
d2 <- apply(decathlon2.active,1,getdistance, center, scale) #funkcijom apply primenjujemo na sve vrste (sve opservacije)
# 2. Compute the cos2. The sum of each row is 1
cos2 <- function(ind.coord, d2){return(ind.coord^2/d2)} #za razliku od varijabli, za racunanje cos2 moramo da skaliramo opservaciju prvo
ind.cos2 <- apply(ind.coord, 2, cos2, d2)
head(ind.cos2[, 1:4])
                  PC1        PC2         PC3         PC4
SEBRLE    0.007530179 0.49747323 0.081325232 0.001386688
CLAY      0.048701249 0.45701660 0.143628117 0.125791741
BERNARD   0.197199804 0.28996555 0.004294015 0.411819183
YURKOV    0.096109800 0.02382571 0.778230322 0.061812637
ZSIVOCZKY 0.001574385 0.57641944 0.239754152 0.001388216
McMULLEN  0.002175437 0.15219499 0.110137872 0.266486530

Doprinos

# Contributions of individuals
#:::::::::::::::::::::::::::::::
contrib <- function(ind.coord, comp.sdev, n.ind){
  100*(1/n.ind)*ind.coord^2/comp.sdev^2
}
ind.contrib <- t(apply(ind.coord, 1, contrib, 
                       res.pca$sdev, nrow(ind.coord)))
head(ind.contrib[, 1:4])
                 PC1        PC2        PC3         PC4
SEBRLE    0.03854254  5.7118249  1.3854184  0.03572215
CLAY      0.65814114 13.8541889  6.4600973  8.55568792
BERNARD   1.86273218  6.1441319  0.1349983 19.57827284
YURKOV    0.79686310  0.4431309 21.4755770  2.57939100
ZSIVOCZKY 0.01178829  9.6816398  5.9748485  0.05231437
McMULLEN  0.01481737  2.3253860  2.4967890  9.13531719

3D scatter plot vizualizacija PCA

Uz pomoć biblioteke plotly moguće je vizualizovati PCA varijable ili opservacije u 3D. Prethodno je potrebno malo pripreme podataka. Ovde vizualizujemo sve zapise iz dataframe-a decathlon2.active u prostoru prve tri glavne komponente, s tim što ćemo pojedinačne zapise obojiti prema takmičenju kojem pripadaju.

library(tidyverse)
# library(janitor)

ind <- get_pca_ind(res.pca)
(ind$coord)
                 Dim.1      Dim.2       Dim.3       Dim.4         Dim.5
SEBRLE       0.1912074 -1.5541282 -0.62836882  0.08205241  1.1426139415
CLAY         0.7901217 -2.4204156  1.35688701  1.26984296 -0.8068483724
BERNARD     -1.3292592 -1.6118687 -0.19614996 -1.92092203  0.0823428202
YURKOV      -0.8694134  0.4328779 -2.47398223  0.69723814  0.3988584116
ZSIVOCZKY   -0.1057450  2.0233632  1.30493117 -0.09929630 -0.1970241089
McMULLEN     0.1185550  0.9916237  0.84355824  1.31215266  1.5858708644
MARTINEAU   -2.3923532  1.2849234 -0.89816842  0.37309771 -2.2433515889
HERNU       -1.8910497 -1.1784614 -0.15641037  0.89130068 -0.1267412520
BARRAS      -1.7744575  0.4125321  0.65817750  0.22872866 -0.2338366980
NOOL        -2.7770058  1.5726757  0.60724821 -1.55548081  1.4241839810
BOURGUIGNON -4.4137335 -1.2635770 -0.01003734  0.66675478  0.4191518468
Sebrle       3.4514485 -1.2169193 -1.67816711 -0.80870696 -0.0250530746
Clay         3.3162243 -1.6232908 -0.61840443 -0.31679906  0.5691645854
Karpov       4.0703560  0.7983510  1.01501662  0.31336354 -0.7974259553
Macey        1.8484623  2.0638828 -0.97928455  0.58469073 -0.0002157834
Warners      1.3873514 -0.2819083  1.99969621 -1.01959817 -0.0405401497
Zsivoczky    0.4715533  0.9267436 -1.72815525 -0.18483138  0.4073029909
Hernu        0.2763118  1.1657260  0.17056375 -0.84869401 -0.6894795441
Bernard      1.3672590  1.4780354  0.83137913  0.74531557  0.8598016482
Schwarzl    -0.7102777 -0.6584251  1.04075176 -0.92717510 -0.2887568007
Pogorelov   -0.2143524 -0.8610557  0.29761010  1.35560294 -0.0150531057
Schoenbeck  -0.4953166 -1.3000530  0.10300360 -0.24927712 -0.6452257128
Barras      -0.3158867  0.8193681 -0.86169481 -0.58935985 -0.7797389436
                  Dim.6       Dim.7        Dim.8        Dim.9      Dim.10
SEBRLE      -0.46389755 -0.20796012  0.043460568 -0.659344137  0.03273238
CLAY         1.30420016 -0.21291866  0.617240611 -0.060125359 -0.31716015
BERNARD     -0.40062867 -0.40643754  0.703856040  0.170083313 -0.09908142
YURKOV       0.10286344 -0.32487448  0.114996135 -0.109524039 -0.11969720
ZSIVOCZKY    0.89554111  0.08825624 -0.202341299 -0.523103099 -0.34842265
McMULLEN     0.18657283  0.47828432  0.293089967 -0.105623196 -0.39317797
MARTINEAU   -0.45666350 -0.29975522 -0.291628488 -0.223417655 -0.61640509
HERNU        0.43623496 -0.56609980 -1.529404317  0.006184409  0.55368016
BARRAS       0.09026010  0.21594095  0.682583078 -0.669282042  0.53085420
NOOL         0.49716399 -0.53205687 -0.433385655 -0.115777808 -0.09622142
BOURGUIGNON -0.08200220 -0.59833739  0.563619921  0.525814030  0.05855882
Sebrle      -0.08279306  0.01016177 -0.030585843 -0.847210682  0.21970353
Clay         0.77715960  0.25750851 -0.580638301  0.409776590 -0.61601933
Karpov      -0.32958134 -1.36365568  0.345306381  0.193055107  0.21721852
Macey       -0.19728082 -0.26927772 -0.363219506  0.368260269  0.21249474
Warners     -0.55673300 -0.26739400 -0.109470797  0.180283071  0.24208420
Zsivoczky   -0.11383190  0.03991159  0.538039776  0.585966156 -0.14271715
Hernu       -0.33168404  0.44308686  0.247293566  0.066908586 -0.20868256
Bernard     -0.32806564  0.36357920  0.006165316  0.279488675  0.32067773
Schwarzl    -0.68891640  0.56568604 -0.687053339 -0.008358849 -0.30211546
Pogorelov   -1.59379599  0.78370119 -0.037623661 -0.130531397 -0.03697576
Schoenbeck   0.16172381  0.85752368 -0.255850722  0.564222295  0.29680481
Barras       1.17415412  0.94512710  0.365550568  0.102255763  0.61186706
# kordsDF3$competition <- decathlon2[1:23,]$Competition

Competition <- decathlon2 %>% 
  select(Competition) %>% 
  slice(1:23) %>% 
  pull(Competition)
  


kordsDF3 <- as.data.frame(ind$coord) %>% 
  rownames_to_column(var = "competitor") %>%
  mutate(competitor=str_to_sentence(str_squish(competitor))) %>% 
  mutate(competitor = factor(competitor, levels = sort(unique(competitor)))) %>%
  add_column(Competition=Competition) %>% 
  janitor::clean_names() %>% 
  select(competitor,dim_1, dim_2, dim_3, competition) %>% 
  arrange(competitor)


# kordsDF3$Competition <- Competition

Funkcija htmlwidgets snima grafika kao poseban HTML fajl.

library(plotly)
library(ggsci)

bojeLancet <- pal_lancet()(5)

labele3 <- vector(mode = "character", 3)
eig_val <- get_eigenvalue(res.pca)

for (i in 1:3) {
  labele3[i] <-
    str_c("PC", i, " (", round(eig_val$variance.percent[i], 2), "%)", sep = "")
}


plotliGrafik <- plot_ly(kordsDF3, x = ~dim_1, y = ~dim_2, z = ~dim_3, color = ~competition, colors = bojeLancet,
                   marker = list(size=8,opacity = 0.5,line = list(color='rgba(0, 0, 0, 1)',width = 1.5,opacity=1)), text = ~paste('Competitor:<b>', competitor,'</b><br>Competition:<em>', competition,'</em>')) %>%
  add_markers() %>% 
  layout(scene = list(xaxis = list(title = labele3[1]),
                      yaxis = list(title = labele3[2]),
                      zaxis = list(title = labele3[3])))


htmlwidgets::saveWidget(as_widget(plotliGrafik), here::here("figures", "plotliGrafik.html"))
plotliGrafik

Rotacija komponenti

# devtools::install_github("thomasp85/patchwork")
library(patchwork)
pca_unrotated <- psych::principal(decathlon2.active, rotate="none", nfactors=2, scores=TRUE)

pca_rotated <- psych::principal(decathlon2.active, rotate="varimax", nfactors=2, scores=TRUE)

pca_rotated2 <- psych::principal(decathlon2.active, rotate="oblimin", nfactors=2, scores=TRUE)

unrotateSol <- biplot.psych(pca_unrotated, main="Unrotated solution")

varimaxSol <- biplot.psych(pca_rotated, main="Varimax rotation")

obliminSol <- biplot.psych(pca_rotated2, main="Oblimin rotation")

Prikaz rotacija.