Example Crime
Beispiel: Crime Data (Everitt, 2010)
Beispiel crime data [Everitt (2010, p. 196)][everitt]
wir sehen:
Hauptkomponentenanalyse kann auf Basis der Rohwertmatrix oder auf Basis der Varianz-Kovarianz-Matrix berechnet werden.
individuelle Berechnung von Komponenten-Scores sind möglich auf Basis der Factor-Loadings, hier gezeigt für den State NY
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.3 ✓ purrr 0.3.4
## ✓ tibble 3.1.0 ✓ dplyr 1.0.5
## ✓ tidyr 1.1.3 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/be/crime.txt")
head(dd)## state murder rape robbery assault burglary theft vehicules
## 1 ME 2.0 14.8 28 102 803 2347 164
## 2 NH 2.2 21.5 24 92 755 2208 228
## 3 VT 2.0 21.8 22 103 949 2697 181
## 4 MA 3.6 29.7 193 331 1071 2189 906
## 5 RI 3.5 21.4 119 192 1294 2568 705
## 6 CT 4.6 23.8 192 205 1198 2758 447# get means of types of crimes
dd %>% dplyr::select(murder:vehicules) %>% apply(., 2, mean)## murder rape robbery assault burglary theft vehicules
## 7.25098 34.21765 154.09804 283.35294 1207.07843 2941.96078 393.84314# we construct a pseudo variable m_crime
# which is sort of an average level of crimes for each state
dd$m_crime <- dd %>% dplyr::select(murder:vehicules) %>% scale() %>% rowMeans()
# covariance matrix
( var_cov <- dd %>% dplyr::select(murder:vehicules) %>% var() ) ## murder rape robbery assault burglary theft
## murder 23.20215 40.59768 533.7049 558.2056 1179.872 1327.224
## rape 40.59768 212.31228 1064.0942 1423.6716 4433.301 7060.497
## robbery 533.70490 1064.09424 18993.3702 15134.7247 32058.152 42034.064
## assault 558.20565 1423.67165 15134.7247 22004.3129 44452.372 57988.894
## burglary 1179.87192 4433.30059 32058.1522 44452.3718 177912.834 246021.443
## theft 1327.22404 7060.49671 42034.0639 57988.8941 246021.443 582812.838
## vehicules 616.80016 1855.30482 24216.2757 21156.8565 54570.013 65869.474
## vehicules
## murder 616.8002
## rape 1855.3048
## robbery 24216.2757
## assault 21156.8565
## burglary 54570.0125
## theft 65869.4737
## vehicules 50007.3749# pca based on raw data
# you may notice small differences in sd compared to pca based on covariance matrix
pca.h <- dd %>% dplyr::select(murder:vehicules) %>% princomp()
print(summary(pca.h),digits=3,loadings=TRUE)## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4
## Standard deviation 841.1698854 290.1020494 174.1630176 96.85347460
## Proportion of Variance 0.8471205 0.1007580 0.0363153 0.01123074
## Cumulative Proportion 0.8471205 0.9478785 0.9841938 0.99542455
## Comp.5 Comp.6 Comp.7
## Standard deviation 61.101673321 9.1116054841 2.293672e+00
## Proportion of Variance 0.004469758 0.0000993957 6.298548e-06
## Cumulative Proportion 0.999894306 0.9999937015 1.000000e+00
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## murder 0.998
## rape 0.998
## robbery 0.276 0.361 0.396 0.794
## assault 0.106 0.278 0.106 0.777 -0.544
## burglary 0.427 0.630 -0.629 -0.135
## theft 0.885 -0.435 0.163
## vehicules 0.128 0.510 0.661 -0.469 -0.256# pca based on covariance matrix
pca.cov <- princomp(var_cov)
print(summary(pca.cov),digits=3,loadings=TRUE)## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4
## Standard deviation 2.144662e+05 2.568336e+04 1.059763e+04 3.296795e+03
## Proportion of Variance 9.832221e-01 1.410060e-02 2.400774e-03 2.323367e-04
## Cumulative Proportion 9.832221e-01 9.973227e-01 9.997235e-01 9.999558e-01
## Comp.5 Comp.6 Comp.7
## Standard deviation 1.437771e+03 2.143431e+01 6.972908e-06
## Proportion of Variance 4.418893e-05 9.820934e-09 1.039350e-21
## Cumulative Proportion 1.000000e+00 1.000000e+00 1.000000e+00
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## murder 1.000
## rape -0.999
## robbery 0.214 0.367 0.421 0.799
## assault 0.242 0.799 -0.535
## burglary 0.396 0.738 -0.518 -0.147
## theft 0.908 -0.403 0.118
## vehicules 0.434 0.758 -0.403 -0.258# compute New Yorks component 1 manually
dd.crimes <- dd %>% dplyr::select(murder:vehicules)
# multiply difference from mean by factor loading
sum((dd.crimes[dd$state=="NY",] - apply(dd.crimes, 2, mean)) * pca.h$loadings)## [1] 65.20597# check against scores provided by princomp for NY
princomp(dd.crimes, scores=T)$scores[which(dd$state=="NY"),1]## Comp.1
## 65.20597# princomp(dd.crimes, scores=T)$scores[7,] all Component scores for NY
# add factor scores to data.frame
dd <- cbind(dd, princomp(dd.crimes, scores=T)$scores)
# a visualization of the parallel analysis
crimes.parallel <- psych::fa.parallel(dd.crimes, fa="pc")
## Parallel analysis suggests that the number of factors = NA and the number of components = 1# a possible visualization of the components
dd %>% ggplot(aes(x=pca.h$scores[,1], y=pca.h$scores[,2], color=m_crime)) + geom_point()
‘Loadings’ zeigt die Gewichte, mit denen die Komponenten-Werte gebildet werden können. ‘Proportion Var’ ist der Anteil an erklärter Varianz, den die jeweilige Komponente erklärt. ‘Cumulative Var’ ist der kumulierte Anteil an der Varianz, der durch alle höheren Komponenten inkl. der aktuellen erkärt wird.
Ein Komponenten-Wert ist eine Linearkombination aus den Rohwerten der Variablen und den Gewichten (loadings) für eine Komponente. Konkret wird oben die Abweichungen einer Beobachtung (state) vom Gruppenmittelwert multipliziert mit dem zugehörigen loading. Dies passiert für alle Variablen in der Komponente und wird aufaddiert.
Die Komponenten-Werte können im Quell-Data-Frame mit abgespeichert werden. Dies ist u. U. sinnvoll, wenn man die Anzahl der Variablen verringern möchte und mit den Komponenten weiterrechnen will.
Die Parallelanalyse empfiehlt die nur eine Komponente.
Die Grafik zeigt, dass die Komponenten etwas mit unserem mittleren Kriminalitätsniveau m_crime zu tun haben.