rmd

Beispiel: Crime Data (Everitt, 2010)

Beispiel crime data [Everitt (2010, p. 196)][everitt]

wir sehen:

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.