ZADANIE 7
Kristína Mečiarová
V tomto zadaní budeme opäť pracovať s dátovým súborom
uscrime z knižnice SMS data. Súbor obsahuje 50
pozorovaní a 11 premenných s informáciami o počte výskytov rôznych
trestných činov v rámci konkrétnych regiónov USA, ich rozlohe a
populácii.
Analýza hlavných komponentov
V predchádzajúcej úlohe sme sa zaoberali analýzou hlavných komponentov. Nasledujúca tabuľka zrhrňuje výsledky pre všetkých 9 komponentov. Vidíme, že na pokrytie zhruba 92% variability nám stačí prvých 5 komponentov, pričom prvý komponent vyjadruje približne 50% variability dát. Na zníženie dimenzionality dát a zároveň dostatočné vysvetlenie variability by nám stačilo prvých 5 komponentov namiesto pôvodných 9.
library(DescTools)
library(kableExtra)
library(knitr)
data(uscrime, package = "SMSdata")
data<-uscrime
pcdata<-data[,-c(10,11)]
pcdata<- scale(pcdata, center = TRUE, scale = TRUE)
pc<-prcomp(pcdata)
tab<-summary(pc)$importance
tab %>%
kbl() %>%
kable_styling()| PC1 | PC2 | PC3 | PC4 | PC5 | PC6 | PC7 | PC8 | PC9 | |
|---|---|---|---|---|---|---|---|---|---|
| Standard deviation | 2.112468 | 1.211975 | 1.074463 | 0.8535431 | 0.6728768 | 0.5266941 | 0.4740118 | 0.365127 | 0.3121157 |
| Proportion of Variance | 0.495840 | 0.163210 | 0.128270 | 0.0809500 | 0.0503100 | 0.0308200 | 0.0249700 | 0.014810 | 0.0108200 |
| Cumulative Proportion | 0.495840 | 0.659040 | 0.787320 | 0.8682700 | 0.9185700 | 0.9494000 | 0.9743600 | 0.989180 | 1.0000000 |
Na obrázkoch nižšie, v ktorých využívame 2. a 3. komponent na zobrazenie dát a 3. a 4. komponent (=rôzne “uhly” pohľadu na 9-rozmerné dáta), môžeme pozorovať zhlukovanie dát podľa jednotlivých regiónov. Zdá sa, že zatiaľ čo v južných regiónoch prevládajú vraždy a útoky, v západných sú to krádeže a vlámania a v stredozápadných regiónoch krádeže áut a prepady.
library("ggfortify")
autoplot(pc, data = data, loadings.label = TRUE, colour = 'reg',x=2,y=3, loadings = T)
autoplot(pc, data = data, loadings.label = TRUE, colour = 'reg',x=3,y=4, loadings = T)
Faktorová analýza
Pracujeme s 9-rozmerným náhodným vektorom a v rámci faktorovej analýzy sa pozrieme na vhodný počet faktorov. Keďže pri 6 faktoroch už dostávame záporný počet stupňov voľnosti pri 9 premenných (\(d=\frac{1}{2}(p-k)^2-\frac{1}{2}(p+k)\)), budeme skúmať iba vlastnoti 2,3,4 a 5 faktorov. Nasledujúca tabuľka zachytáva kumulatívny rozptyl a p-hodnoty pre daný počet faktorov.
fa1 <- factanal(pcdata, factors = 2, rotation="varimax")
fa2 <- factanal(pcdata, factors = 3, rotation="varimax")
fa3 <- factanal(pcdata, factors = 4, rotation="varimax")
fa4 <- factanal(pcdata, factors = 5, rotation="varimax")
f1<-c(0.59, fa1$PVAL)
f2<-c(0.69, fa2$PVAL)
f3<-c(0.79, fa3$PVAL)
f4<-c(0.86, fa4$PVAL)
tabul<-rbind(f1,f2,f3,f4)
colnames(tabul)<-c("cum. var", "p-hodnota")
rownames(tabul)<-c("2","3","4","5")
print(tabul)## cum. var p-hodnota
## 2 0.59 0.001115917
## 3 0.69 0.081426919
## 4 0.79 0.292347438
## 5 0.86 0.127794059Vidíme, že zamietame nulovú hypotézu, že sú postačujúce 2 faktory a 3
faktory sú relatívne blízko hladiny významnosti. Na základe najvyššej
p-hodnoty, dostatočne vysokého kumulatívneho rozptylu a vhodnej redukcie
dimenzionality si zvolíme 4 faktory.
V nižšie uvedenej tabuľke pozorujeme, že faktor 1 je vysoko korelovaný s premennými burglary a larceny, faktor 2 s premennými murder a assault, faktor 3 s popu1985 a robbery a faktor 4 s land area.
print(fa3, digits=2, cutoff=.6, sort=TRUE)##
## Call:
## factanal(x = pcdata, factors = 4, rotation = "varimax")
##
## Uniquenesses:
## land area popu 1985 murder rape robbery assault burglary larceny
## 0.00 0.52 0.19 0.27 0.16 0.09 0.10 0.18
## autotheft
## 0.41
##
## Loadings:
## Factor1 Factor2 Factor3 Factor4
## rape
## burglary 0.85
## larceny 0.88
## autotheft
## murder 0.88
## assault 0.87
## popu 1985 0.63
## robbery 0.80
## land area 0.98
##
## Factor1 Factor2 Factor3 Factor4
## SS loadings 2.35 2.01 1.63 1.08
## Proportion Var 0.26 0.22 0.18 0.12
## Cumulative Var 0.26 0.48 0.67 0.79
##
## Test of the hypothesis that 4 factors are sufficient.
## The chi square statistic is 7.32 on 6 degrees of freedom.
## The p-value is 0.292par(mfrow = c(1,2))
load <- fa3$loadings[,1:2]
plot(load,type="n")
text(load,labels=names(uscrime[,-c(10,11)]),cex=.3)
lines(c(-1,1),c(0,0), lty = 3)
lines(c(0,0),c(-1,1), lty = 3)
load <- fa3$loadings[,3:4]
plot(load,type="n")
text(load,labels=names(uscrime[,-c(10,11)]),cex=.3)
lines(c(-1,1),c(0,0), lty = 3)
lines(c(0,0),c(-1,1), lty = 3)
par(mfrow = c(1,1))Na obrázkoch pozorujeme, že premenná rape je kladne a približne rovnako korelovaná so všetkými faktormi a premenná autotheft je pomerne vysoko korelovaná s faktorom 1 a 3. Na základe tejto analýzy by sme mohli v potrebnom modeli uvažovať faktor 1, ktorý by mohol reprezentovať vlámania (burglary, larceny, autotheft), faktor 2, ktorý by reprezentoval napadnutia a vraždy (assault, murder), faktor 3, ako reprezentanta krádeží (robbery, rape, popul1985) a ďalej zvlášť premennú land area, na ktorej sa zdá, že nezávisí intenzita výskytov týchto trestých činov.
rbPal <- colorRampPalette(c('red','blue'))
loadings <- data.frame(fa3$loadings[,1:4])
loadings$colF1 <- rbPal(20)[as.numeric(cut(loadings[,1],breaks = seq(-1,1,length = 20)))]
loadings$nam1 <- paste(row.names(loadings), round(loadings[,1], digits = 3), sep = "\n")
loadings$colF2 <- rbPal(20)[as.numeric(cut(loadings[,2],breaks = seq(-1,1,length = 20)))]
loadings$nam2 <- paste(row.names(loadings), round(loadings[,2], digits = 3), sep = "\n")
loadings$colF3 <- rbPal(20)[as.numeric(cut(loadings[,3],breaks = seq(-1,1,length = 20)))]
loadings$nam3 <- paste(row.names(loadings), round(loadings[,2], digits = 3), sep = "\n")
loadings$colF4 <- rbPal(20)[as.numeric(cut(loadings[,4],breaks = seq(-1,1,length = 20)))]
loadings$nam4 <- paste(row.names(loadings), round(loadings[,2], digits = 3), sep = "\n")
par(mar = c(5, 5, 5, 4) + 0.4)
barplot(abs(loadings[,1]), col = loadings$colF1, ylim = c(0,1), ylab = "Factor 1", names.arg = loadings$nam1, cex.names = 0.5)
barplot(abs(loadings[,2]), col = loadings$colF2, ylim = c(0,1), ylab = "Factor 2", names.arg = loadings$nam2, cex.names = 0.5)
barplot(abs(loadings[,3]), col = loadings$colF3, ylim = c(0,1), ylab = "Factor 3", names.arg = loadings$nam2, cex.names = 0.5)
barplot(abs(loadings[,4]), col = loadings$colF4, ylim = c(0,1), ylab = "Factor 4", names.arg = loadings$nam2, cex.names = 0.5)