pacman::p_load(tidyverse, ggplot2, ggthemes, psych, haven, EFAutilities, knitr)4 Hauptkomponentenanalyse (PCA)
4.1 Setup
Wir können die im einführenden Kapitel zur Exploratorischen Datenreduktion bereits bereinigten Daten zur Lebenszufriedenheit einlesen, entweder aus dem data Ordner oder auch direkt von GitHub:
# Aus dem lokalen datenordner
# ls_clean <- read_csv("data/ls_clean.csv")
# Aus GitHub
ls_clean <- read_csv("https://raw.githubusercontent.com/methodenlehre/data/master/statistik_IV/ls_clean.csv")4.2 Parallelanalyse
Zuerst führen wir eine Parallelanalyse durch, um die Anzahl der zu extrahierenden Komponenten zu bestimmen. Weil wir die Parallelanalyse zunächst nur für die PCA berechnen, spezifizieren wir das Argument fa jetzt als "pc" (für principal component). Wenn wir das nicht spezifizieren, wird die Parallelanalyse sowohl für die PCA als auch für die Exploratorische Faktorenanalyse (EFA) zusammen dargestellt werden, was wir hier nicht wollen.
Die Anzahl Iterationen (Zufalls-Samples), die per Default 20 ist, heben wir auf 1000 an, um eine möglichst stabile Lösung zu erhalten.
Als Kriterium für die Anzahl Komponenten wählen wir das Quantil quant = 0.5, also den Median. Das bedeutet, dass die Eigenwerte der tatsächlichen Komponenten, um als substantielle (zu extrahierende) Hauptkomponenten zu gelten, grösser sein müssen als der Eigenwerts-Median der jeweiligen Komponente aus den simulierten Daten (Median von jeweils 1000 simulierten Eigenwerten).
Das unter quant eingestellte Quantil der simulierten Eingenwerteverteilungen bezieht sich allerdings nur auf die im Text ausgegebene Schlussfolgerung, wie viele Komponenten extrahiert werden sollten (im Sinne eines Vergleichs des jeweiligen empirischen Eigenwerts mit diesem Quantil). Die Voreinstellung (default) für diesen Test ist quant = 0.95, ein empirischer Eigenwert wird also nur dann als relevant betrachtet, wenn er den grössten 5 % der simulierten Eigenwerte entspricht.
Die im ausgegeben Plot der Parallelanalyse dargestellte Linie bezieht sich dagegen unabhängig von der quant-Einstellung immer auf den Mittelwert der simulierten Eigenwerteverteilungen. Daher stimmt die ausgegebene Schlussfolgerung nicht in jedem Fall mit der nach dem Plot zu ziehenden Schlussfolgerung überein. Da wir quant = 0.5 verwenden, sollte das für unsere Beispiele aber kaum eine Rolle spielen (da der Median normalerweise sehr nahe am Mittelwert liegt).
parallelAnalyse <- fa.parallel(ls_clean, n.iter = 1000, fa = "pc", quant = 0.5)
Parallel analysis suggests that the number of factors = NA and the number of components = 3 Die Parallelanalyse der PCA zeigt deutlich drei zu extrahierende Komponenten an, da die vierte Komponente aus random Data einen höheren Eigenwert als die im Datenset gefundene vierte Komponente aufweist. Drei Komponenten würden wir wohl auch nach dem Scree-Kriterium extrahieren, da sich nach der dritten Komponente ein deutlicher “Knick” im Verlauf der Eigenwerte zeigt. Nicht zuletzt sollten auch nach dem Kaiser-Kriterium drei Komponenten extrahiert werden (da der Eigenwert der vierten Komponente bereits < 1 ist).
4.3 Unrotierte Lösung mit drei Hauptkomponenten
Zuerst berechnen wir eine unrotierte PCA:
# Unrotierte PCA
principal(ls_clean, nfactors = 3, rotate = "none")Principal Components Analysis
Call: principal(r = ls_clean, nfactors = 3, rotate = "none")
Standardized loadings (pattern matrix) based upon correlation matrix
PC1 PC2 PC3 h2 u2 com
leben1_schule 0.45 0.42 0.12 0.40 0.60 2.1
leben3_schule 0.41 0.60 0.38 0.68 0.32 2.5
leben4_schule 0.50 0.54 0.45 0.75 0.25 2.9
leben2_selbst 0.50 -0.43 0.39 0.58 0.42 2.9
leben6_selbst 0.57 -0.51 0.20 0.63 0.37 2.2
leben5_freunde 0.62 -0.33 0.22 0.54 0.46 1.8
leben7_freunde 0.48 -0.34 0.09 0.36 0.64 1.9
leben8_familie 0.62 0.08 -0.56 0.71 0.29 2.0
leben9_familie 0.60 0.14 -0.63 0.78 0.22 2.1
leben10_familie 0.58 0.05 -0.33 0.44 0.56 1.6
PC1 PC2 PC3
SS loadings 2.90 1.54 1.43
Proportion Var 0.29 0.15 0.14
Cumulative Var 0.29 0.44 0.59
Proportion Explained 0.49 0.26 0.24
Cumulative Proportion 0.49 0.76 1.00
Mean item complexity = 2.2
Test of the hypothesis that 3 components are sufficient.
The root mean square of the residuals (RMSR) is 0.1
with the empirical chi square 211.22 with prob < 5.6e-35
Fit based upon off diagonal values = 0.85Im Output bekommen wir unter SS loadings (Sum of Squared Loadings) zum einen die Eigenwerte der Komponenten (vgl. Scree-Plot Parallelanalyse). Zum anderen die Anteile bzw. die kumulativen Anteile der durch die Komponenten erklärten Varianz, also den Eigenwert dividiert durch die Anzahl Variablen.
In den letzten beiden Zeilen erhält man noch die Anteile bzw. kumulativen Anteile der erklärten Varianz der jeweiligen Komponenten an der insgesamt durch alle Komponenten erklärten Varianz. Z.B. werden durch die drei Komponenten insgesamt 59 % der Gesamtvarianz erklärt (Cumulative Var \(PC3 = 0.59\)). Also ist die Proportion explained für \(PC1 = 0.29/0.59 = 0.49\) und die für \(PC2 = 0.15/0.59 = 0.25\).
Des weiteren wird noch ein Chi-Quadrat-Test zum Test dafür, dass drei Komponenten ausreichen, ausgegeben. Wir betrachten diesen hier nicht, da dieser Test bei der PCA umstritten ist, weil er trotz der rein datenreduzierenden und damit deskriptiven Funktion der PCA einen Test der verbleibenden Residual(ko-)varianz in der Population durchführt. Mit dieser Art von Test beschäftigen wir uns im folgenden Kapitel zur ML-EFA.
Bei der unrotierten Lösung sieht man relativ grosse Querladungen der einzelnen Items. Das ist normal: Die Komponenten wurden so extrahiert, dass die erste Komponente möglichst viel erklärt, d.h. die Ladungen aller Variablen auf der ersten Komponente wurden maximiert (hier z.B. keine Ladung kleiner als 0.41).
Die variable com im Output steht hier für Komplexität der Items. Je höher die Komplexität, desto mehr laden sie auf verschiedenen Komponenten. Die Berechnung dieser statistischen Grösse behandeln wir hier nicht.
h2 steht für Kommunalität, also den Anteil der Varianz einer beobachteten Variablen, der durch die extrahierten Komponenten erklärt wird. u2 steht für Einzigartigkeit (uniqueness), also für den durch die Komponenten unerklärten Varianzanteil der Variablen (u2 = 1 - h2).
Von Hand berechnet:
Um eine etwas grössere Genauigkeit zu erzielen, sind in folgender Tabelle die Ladungen nochmals auf drei Dezimalstellen gerundet.
| PC1 | PC2 | PC3 | |
|---|---|---|---|
| leben1_schule | 0.453 | 0.422 | 0.115 |
| leben3_schule | 0.408 | 0.605 | 0.383 |
| leben4_schule | 0.500 | 0.545 | 0.451 |
| leben2_selbst | 0.501 | -0.429 | 0.387 |
| leben6_selbst | 0.574 | -0.510 | 0.202 |
| leben5_freunde | 0.616 | -0.334 | 0.218 |
| leben7_freunde | 0.483 | -0.338 | 0.093 |
| leben8_familie | 0.625 | 0.084 | -0.563 |
| leben9_familie | 0.604 | 0.140 | -0.630 |
| leben10_familie | 0.576 | 0.050 | -0.330 |
\(~\) Berechnung des Eigenwerts SS loadings der zweiten Komponente:
\(\begin{aligned} Var(H_2) &= \lambda_{12}^2 + \lambda_{22}^2 + \lambda_{32}^2 + \lambda_{42}^2 + \lambda_{52}^2 + \lambda_{62}^2 + \lambda_{72}^2 + \lambda_{82}^2 + \lambda_{92}^2 + \lambda_{102}^2 \\ &= 0.422^2 + 0.605^2 + 0.545^2 + (-0.429)^2 + (-0.51)^2 + (-0.334)^2 + (-0.338)^2 \\ &{~~~~~}+ 0.084^2 + 0.14^2 + 0.05^2 \\ &= 1.54 \end{aligned}\)
\(~\)
Berechnung der Kommunalität h2 des Items leben2_selbst:
\(\begin{aligned} \hat{H}(Y_4) &= \lambda_{41}^2 + \lambda_{42}^2 + \lambda_{43}^2 \\ &= 0.501^2 + (-0.429)^2 + 0.387^2 \\ &= 0.58 \end{aligned}\)
\(~\)
Berechnung der Uniqueness u2 des Items leben8_familie:
\(\begin{aligned} \hat{U}(Y_8) &= 1 - (\lambda_{81}^2 + \lambda_{82}^2 + \lambda_{83}^2) \\ &= 1 - [0.625^2 + 0.084^2 + (-0.563)^2] \\ &= 0.29 \end{aligned}\)
4.4 Orthogonale Rotation mit drei Hauptkomponenten
Jetzt rotieren wir die Komponenten mit der orthogonalen Varimax-Rotation. Dabei werden die Komponenten so rotiert, dass die Items möglichst wenige Querladungen aufweisen, also bestmöglich eine Einfachstruktur darstellen und gleichzeitig die Orthogonalität der Lösung erhalten bleibt (keine Interkorrelationen der Komponenten).
Im Gegensatz zu obliquen Rotationen bleiben hier die Summen der Eigenwerte identisch. Auch wenn in der aktuellen Literatur meist oblique Rotationsmethoden empfohlen werden, sind orthogonale Rotationen wegen ihrer Eigenschaft, dass die Komponenten weiterhin als unabhängige Dimensionen interpretiert werden können, attraktiv (insbesondere wenn sich damit eine Einfachstruktur erreichen lässt).
# Rotierte PCA
pca.3.rotiert <- principal(ls_clean, nfactors = 3, rotate = "varimax")
# Per Default ist das principal-Objekt nicht sortiert. Damit wir die
# Items in der Reihenfolge der Faktorenladungen bekommen, müssen wir
# die Funktion `print.psych()` auf das Objekt anwenden und den
# Parameter `sort = TRUE` setzen.
print.psych(pca.3.rotiert, sort = TRUE)Principal Components Analysis
Call: principal(r = ls_clean, nfactors = 3, rotate = "varimax")
Standardized loadings (pattern matrix) based upon correlation matrix
item RC1 RC3 RC2 h2 u2 com
leben6_selbst 5 0.78 0.13 -0.01 0.63 0.37 1.1
leben2_selbst 4 0.76 -0.05 0.10 0.58 0.42 1.0
leben5_freunde 6 0.70 0.16 0.14 0.54 0.46 1.2
leben7_freunde 7 0.57 0.18 0.02 0.36 0.64 1.2
leben9_familie 9 0.05 0.88 0.07 0.78 0.22 1.0
leben8_familie 8 0.12 0.83 0.08 0.71 0.29 1.1
leben10_familie 10 0.21 0.62 0.14 0.44 0.56 1.3
leben4_schule 3 0.14 0.03 0.85 0.75 0.25 1.1
leben3_schule 2 0.01 0.03 0.82 0.68 0.32 1.0
leben1_schule 1 0.06 0.24 0.58 0.40 0.60 1.4
RC1 RC3 RC2
SS loadings 2.08 1.99 1.80
Proportion Var 0.21 0.20 0.18
Cumulative Var 0.21 0.41 0.59
Proportion Explained 0.35 0.34 0.31
Cumulative Proportion 0.35 0.69 1.00
Mean item complexity = 1.1
Test of the hypothesis that 3 components are sufficient.
The root mean square of the residuals (RMSR) is 0.1
with the empirical chi square 211.22 with prob < 5.6e-35
Fit based upon off diagonal values = 0.85Die Reihenfolge der Darstellung der rotierten Komponenten (RC) ergibt sich aus der Grösse der Eigenwerte/des Varianzanteils nach der Rotation. Jetzt ist der Anteil der erklärten Varianz der drei Komponenten sehr ähnlich, da durch die Rotation das Kriterium der sukzessiv maximalen Varianzaufklärung aufgehoben wurde und die erklärte Varianz mit dem Ziel einer Einfachstruktur möglichst gleichmässig auf die Komponenten verteilt wurde.
Wir können aufgrund folgender Tabelle beispielhaft nochmal einen rotierten Eigenwert von Hand berechnen, diesmal den der ersten Komponente:
| RC1 | RC3 | RC2 | |
|---|---|---|---|
| leben1_schule | 0.057 | 0.243 | 0.578 |
| leben3_schule | 0.014 | 0.031 | 0.823 |
| leben4_schule | 0.139 | 0.027 | 0.854 |
| leben2_selbst | 0.756 | -0.052 | 0.102 |
| leben6_selbst | 0.784 | 0.127 | -0.015 |
| leben5_freunde | 0.701 | 0.163 | 0.143 |
| leben7_freunde | 0.569 | 0.178 | 0.018 |
| leben8_familie | 0.124 | 0.833 | 0.076 |
| leben9_familie | 0.048 | 0.880 | 0.075 |
| leben10_familie | 0.206 | 0.617 | 0.142 |
\(~\)
\(\begin{aligned} Var(H_{1 rotiert}) &= \lambda_{11}^2 + \lambda_{21}^2 + \lambda_{31}^2 + \lambda_{41}^2 + \lambda_{51}^2 + \lambda_{61}^2 + \lambda_{71}^2 + \lambda_{81}^2 + \lambda_{91}^2 + \lambda_{101}^2 \\ &= 0.057^2 + 0.014^2 + 0.139^2 + 0.756^2 + 0.784^2 + 0.701^2 + 0.569^2 + 0.124^2 \\ &{~~~~~}+ 0.048^2 + 0.206^2 \\ &= 2.084 \end{aligned}\)
\(~\)
Daraus folgender Anteil der erklärten Varianz der ersten rotierten Komponente:
Proportion Var \(= \frac{2.084}{10} = 0.2084\)
Die Kommunalität der Items bleibt nach orthogonaler Rotation dagegen unverändert, z.B. für leben2_selbst:
\(\begin{aligned} \hat{H}(Y_4) &= \lambda_{41}^2 + \lambda_{42}^2 + \lambda_{43}^2 \\ &= 0.756^2 + (-0.052)^2 + 0.102^2 \\ &= 0.58 \end{aligned}\)
\(~\)
Tatsächlich gibt es jetzt weniger und deutlich schwächere Querladungen. Entsprechend ist auch die com (Komplexität) der Items gesunken. Jetzt erst zeigt sich eine Einfachstruktur und damit auch, wie die Komponenten inhaltlich zu interpretieren sind: auf der ersten rotierten Komponente (RC1) laden die Items der Inhaltsbereiche Selbst und Freunde, die zweite rotierte Komponente (RC2) ist als Schul-Komponente zu interpretieren und die dritte (RC3) als Familien-Komponente.
Es gibt nur wenige Querladungen. leben10_familie hat eine Ladung von 0.21 auf der Selbst/Freunde-Komponente und leben1_schule hat eine Ladung von 0.24 auf der Familien-Komponente. Alle anderen Ladungen auf nicht zugehörigen Komponenten sind < 0.20.
4.5 Oblique Rotation mit drei Hauptkomponenten
Jetzt benutzen wir als oblique Rotation eine Oblimin-Rotation. Je stärker die Komponenten (tatsächlich) miteinander korrellieren, desto grösser ist der Vorteil von obliquen Rotationen. Diese können anders als orthogonale Rotationen eine Korrelation der latenten Variablen (Komponenten) selbst abbilden und damit bei den Ladungen (die jetzt nicht mehr als Korrelationen, sondern als Partial-Regressionskoeffizienten zu interpretieren sind) eine noch bessere Einfachstruktur erreichen. Falls die Dimensionen tatsächlich orthogonal sein sollten, so wird sich das in Komponentenkorrelationen nahe Null widerspiegeln. Das ist auch der Grund, warum man orthogonale Rotationen nur in seltenen Fällen wirklich benötigt.
# Oblique Rotation
pca.3.rotiert.obli <- principal(ls_clean, nfactors = 3, rotate = "oblimin")
# Sortiert ausgeben lassen:
print.psych(pca.3.rotiert.obli, sort = TRUE)Principal Components Analysis
Call: principal(r = ls_clean, nfactors = 3, rotate = "oblimin")
Standardized loadings (pattern matrix) based upon correlation matrix
item TC1 TC3 TC2 h2 u2 com
leben6_selbst 5 0.79 0.05 -0.08 0.63 0.37 1.0
leben2_selbst 4 0.77 -0.14 0.05 0.58 0.42 1.1
leben5_freunde 6 0.69 0.09 0.08 0.54 0.46 1.1
leben7_freunde 7 0.56 0.12 -0.04 0.36 0.64 1.1
leben9_familie 9 -0.05 0.90 -0.01 0.78 0.22 1.0
leben8_familie 8 0.03 0.84 -0.01 0.71 0.29 1.0
leben10_familie 10 0.13 0.60 0.07 0.44 0.56 1.1
leben4_schule 3 0.06 -0.03 0.86 0.75 0.25 1.0
leben3_schule 2 -0.06 -0.01 0.83 0.68 0.32 1.0
leben1_schule 1 -0.02 0.21 0.56 0.40 0.60 1.3
TC1 TC3 TC2
SS loadings 2.08 2.00 1.79
Proportion Var 0.21 0.20 0.18
Cumulative Var 0.21 0.41 0.59
Proportion Explained 0.35 0.34 0.30
Cumulative Proportion 0.35 0.70 1.00
With component correlations of
TC1 TC3 TC2
TC1 1.00 0.22 0.18
TC3 0.22 1.00 0.17
TC2 0.18 0.17 1.00
Mean item complexity = 1.1
Test of the hypothesis that 3 components are sufficient.
The root mean square of the residuals (RMSR) is 0.1
with the empirical chi square 211.22 with prob < 5.6e-35
Fit based upon off diagonal values = 0.85Gibt es noch substantielle Querladungen? - Das einzige Item, das noch eine gewisse Querladung aufweist, ist leben1_schule. Dieses bezieht sich direkt auf die Schulnoten und zeigt eine Querladung auf der Familien-Komponente (TC3). Auch ist die Ladung dieses Items auf der Schulkomponente (TC2) deutlich niedriger als die der anderen beiden zugehörigen Items. Die Schulnoten scheinen also auch für die Zufriedenheit für den Bereich Familie eine gewisse Bedeutung zu haben. Das ist inhaltlich durchaus nachvollziehbar.
In Bezug auf die Komponentenkorrelationen zeigt sich, dass alle drei Komponenten relativ gleichmässig und gleichzeitig nicht sonderlich stark miteinander korrelieren \((r_{H_1H_2} = 0.18, r_{H_1H_3} = 0.22, r_{H_2H_3} = 0.17)\). Dies kann auf eine übergeordnete Dimension Lebenszufriedenheit hinweisen. Damit werden wir uns demnächst bei der konfirmatorischen Faktorenanalyse (CFA) noch weiter beschäftigen.
Visualisierung:
Wie stark die Items auf Komponenten laden, kann mit factor.plot() dargestellt werden:
factor.plot(pca.3.rotiert.obli, cut = 0.5)
\(~\)
Welche Items “gehören” zu welcher Komponente? Das können wir mit fa.diagram() darstellen:
fa.diagram(pca.3.rotiert.obli)
4.6 Zusammenfassung PCA
Die PCA zeigt also, dass drei Hauptkomponenten insgesamt ca. 60 % der Gesamtvarianz erklären können und dass diese im Vergleich zu den übrigen sieben Komponenten über eine überdurchschnittliche Varianzaufklärung verfügen. Ob dies auch im Sinne einer EFA, die die Korrelationen zwischen den Items vollständig durch Faktoren erkären will, ausreicht, wird sich im nächsten Teil zeigen.
Inhaltlich interessant ist die Kombination der Lebenzufriedenheit mit dem Bereich Selbst und dem Bereich Freunde. Diese beiden Zufriedenheiten scheinen eng miteinander verknüpft zu sein.
4.7 Übung
Für die Übung bearbeiten wir ein Datenset von den 1988 Olympischen Spielen in Seoul. Im Datenset sind die Resultate des olympischen Siebenkampfs (Heptathlon) enthalten. Jede Person hat also insgesamt sieben Resultate. Es handelt sich um Daten von 25 Athletinnen, also um ein sehr kleines n für eine PCA mit 7 Variablen. Das soll uns aber fürs Üben an einem kleinen Beispiel nicht weiter stören. Die Daten sind im Package HSAUR2 gespeichert.
Folgender Code-Chunk lädt erst das Package und dann die Daten heptathlon.
pacman::p_load(HSAUR2)
data("heptathlon", package = "HSAUR2")heptathlon hurdles highjump shot run200m longjump javelin run800m
Joyner-Kersee (USA) 12.69 1.86 15.80 22.56 7.27 45.66 128.51
John (GDR) 12.85 1.80 16.23 23.65 6.71 42.56 126.12
Behmer (GDR) 13.20 1.83 14.20 23.10 6.68 44.54 124.20
Sablovskaite (URS) 13.61 1.80 15.23 23.92 6.25 42.78 132.24
Choubenkova (URS) 13.51 1.74 14.76 23.93 6.32 47.46 127.90
Schulz (GDR) 13.75 1.83 13.50 24.65 6.33 42.82 125.79
Fleming (AUS) 13.38 1.80 12.88 23.59 6.37 40.28 132.54
Greiner (USA) 13.55 1.80 14.13 24.48 6.47 38.00 133.65
Lajbnerova (CZE) 13.63 1.83 14.28 24.86 6.11 42.20 136.05
Bouraga (URS) 13.25 1.77 12.62 23.59 6.28 39.06 134.74
Wijnsma (HOL) 13.75 1.86 13.01 25.03 6.34 37.86 131.49
Dimitrova (BUL) 13.24 1.80 12.88 23.59 6.37 40.28 132.54
Scheider (SWI) 13.85 1.86 11.58 24.87 6.05 47.50 134.93
Braun (FRG) 13.71 1.83 13.16 24.78 6.12 44.58 142.82
Ruotsalainen (FIN) 13.79 1.80 12.32 24.61 6.08 45.44 137.06
Yuping (CHN) 13.93 1.86 14.21 25.00 6.40 38.60 146.67
Hagger (GB) 13.47 1.80 12.75 25.47 6.34 35.76 138.48
Brown (USA) 14.07 1.83 12.69 24.83 6.13 44.34 146.43
Mulliner (GB) 14.39 1.71 12.68 24.92 6.10 37.76 138.02
Hautenauve (BEL) 14.04 1.77 11.81 25.61 5.99 35.68 133.90
Kytola (FIN) 14.31 1.77 11.66 25.69 5.75 39.48 133.35
Geremias (BRA) 14.23 1.71 12.95 25.50 5.50 39.64 144.02
Hui-Ing (TAI) 14.85 1.68 10.00 25.23 5.47 39.14 137.30
Jeong-Mi (KOR) 14.53 1.71 10.83 26.61 5.50 39.26 139.17
Launa (PNG) 16.42 1.50 11.78 26.16 4.88 46.38 163.43
score
Joyner-Kersee (USA) 7291
John (GDR) 6897
Behmer (GDR) 6858
Sablovskaite (URS) 6540
Choubenkova (URS) 6540
Schulz (GDR) 6411
Fleming (AUS) 6351
Greiner (USA) 6297
Lajbnerova (CZE) 6252
Bouraga (URS) 6252
Wijnsma (HOL) 6205
Dimitrova (BUL) 6171
Scheider (SWI) 6137
Braun (FRG) 6109
Ruotsalainen (FIN) 6101
Yuping (CHN) 6087
Hagger (GB) 5975
Brown (USA) 5972
Mulliner (GB) 5746
Hautenauve (BEL) 5734
Kytola (FIN) 5686
Geremias (BRA) 5508
Hui-Ing (TAI) 5290
Jeong-Mi (KOR) 5289
Launa (PNG) 4566# Im Datensatz sind die Namen der Athletinnen in den "rownames" gespeichert.
# rownames sind eine R Eigenschaft (ein Attribut von Datenframes), die wir
# bisher nicht kennengelernt haben und auf die wir auch jetzt nicht näher
# eingehen wollen, da wir es bevorzugen, eine richtige Variable mit den
# Namen zu haben. Daher:
heptathlon <- heptathlon |>
rownames_to_column(var = "name")
# Anschauen:
head(heptathlon) name hurdles highjump shot run200m longjump javelin run800m
1 Joyner-Kersee (USA) 12.69 1.86 15.80 22.56 7.27 45.66 128.51
2 John (GDR) 12.85 1.80 16.23 23.65 6.71 42.56 126.12
3 Behmer (GDR) 13.20 1.83 14.20 23.10 6.68 44.54 124.20
4 Sablovskaite (URS) 13.61 1.80 15.23 23.92 6.25 42.78 132.24
5 Choubenkova (URS) 13.51 1.74 14.76 23.93 6.32 47.46 127.90
6 Schulz (GDR) 13.75 1.83 13.50 24.65 6.33 42.82 125.79
score
1 7291
2 6897
3 6858
4 6540
5 6540
6 6411Die Variable score ist die finale Wertung der Leistung einer Athletin. Diese interessiert uns momentan nicht, wir belassen sie aber für später im Datensatz.
Das Datenframe hat also 7 Variablen, die uns im Moment interessieren:
hurdles: Hürdenlaufhighjump: Hochsprungshot: Kugelstossenrun200m: 200m Lauflongjump: Weitsprungjavelin: Speerwurfrun800m: 800m Lauf
Visualisierung
Wir visualisieren die Daten, um uns einen Überblick zu verschaffen. Dies geht aber nicht so einfach, weil sich die Skalen offensichtlich stark unterscheiden. Zusätzlich sind bei manchen Sportarten kleine Werte erwünscht (run200) und bei anderen grosse Werte (longjump). Folgender Code-Chunk stellt sicher, dass für alle Sportarten grosse Werte wünschenswert sind:
# Wir subtrahieren den Wert von "falsch kodierten" Variablen vom Maximalwert
# dieser Variable:
heptathlon <- heptathlon |>
mutate(
hurdles = max(hurdles) - hurdles,
run200m = max(run200m) - run200m,
run800m = max(run800m) - run800m
)
# Die Werte der veränderten Variablen lassen sich jetzt zwar nicht mehr so
# einfach interpretieren, dafür zeigen sie nun alle in dieselbe Richtung.Jetzt können wir uns die Daten in einem Sternenplot anschauen:
par(mar = c(4, 4, 0.1, 0.1))
heptathlon |>
select(-name, -score) |>
stars(
draw.segments = TRUE,
key.loc = c(-0.7, 5),
nrow = 5,
scale = TRUE,
labels = heptathlon$name
)
Der Sternenplot ist eine simple Methode, um multivariate Daten mit bis zu ungefähr 10 Variablen zu Visualisieren. Der Plot bietet sich besonders an, wenn wir nur wenige Personen haben.
In unserem Fall hat jede Variable ein Segment (und eine Farbe) pro Person. Je grösser dieses Segment ist, desto höher ist der Wert dieser Person auf dieser Variable. Weil die Variablen auf unterschiedlichen Skalen gemessen wurden, setzen wir scale = TRUE. Dies ist bereits der default dieser Funktion und wird hier lediglich als zusätzliche Information explizit erwähnt.
Aufgabe 1
Im Datenframe gibts einen Outlier. Die Werte dieser Athletin unterscheiden sich stark von denen aller anderen Teilnehmerinnen. Entfernen Sie den Outlier basierend auf dem starplot von oben.
Solution
Der Outlier ist Launa. Sie ist sehr gut im Speerwerfen, hat aber sonst sehr kleine Werte.
heptathlon <- heptathlon |>
filter(name != "Launa (PNG)")Aufgabe 2
Das Ziel dieser Übung ist die Reduktion der 7 Scores auf möglichst eine Dimension, um eine informative Rangreihenfolge der Athletinnen zu erstellen, die die Leistung über alle sieben Disziplinen hinweg möglichst optimal widerspiegeln soll. In einem ersten Schritt möchten wir wissen, ob sich diese Daten überhaupt auf weniger bzw. auf eine einzige Dimension reduzieren lassen.
Führen Sie dafür eine Parallelanalyse (für PCA) durch und erklären Sie, wie viele Hauptkomponenten angebracht sind.
Solution
library(psych) # falls nicht mehr geladen
heptathlon |>
select(-name, -score) |>
fa.parallel(fa = "pc", n.iter = 1000, quant = 0.5)
Parallel analysis suggests that the number of factors = NA and the number of components = 1 Die Parallelanalyse empfiehlt eine Lösung mit nur einer Hauptkomponente. Das bedeutet für uns, dass mit einer Hauptkomponente bereits genügend Varianz eingefangen wird, um einen guten Siebenkampf-Score zu berechnen, es aus empirischer Sicht also nicht nötig ist, mehrere Disziplin-Dimensionen zu betrachten (man könnte natürlich auch dann einen Gesamtscore berechnen, wenn sich herausstellen sollte, dass man eigentlich mehrere Komponenten benötigt - dann würde dieser Gesamtscore einfach keine homogene Siebenkampf-Fähigkeit widerspiegeln).
Aufgabe 3
Berechnen Sie eine PCA mit einer Hauptkomponente über die sieben Disziplinen.
Wieviel Prozent Varianz wird mit dieser Hauptkomponente erklärt? (Output)
Berechnen Sie ausserdem den Eigenwert dieser Komponente von Hand und daraus dann nochmals die % erklärte Varianz. (Tipp: die Kommunalitäten h2 ensprechen jetzt den quadrierten Ladungen auf der ersten und einzigen Komponente!)
Solution
# Das letzte Argument könnte wir uns auch sparen, da eine einzige
# Komponente nicht rotiert werden kann
pca_heptathlon <- heptathlon |>
select(-name, -score) |>
principal(factors = 1, rotate = "none")
# Ergebnisobjekt drucken, um Ladungen nach Grösse sortieren zu können
print.psych(pca_heptathlon, sort = TRUE)Principal Components Analysis
Call: principal(r = select(heptathlon, -name, -score), rotate = "none",
factors = 1)
Standardized loadings (pattern matrix) based upon correlation matrix
V PC1 h2 u2 com
longjump 5 0.94 0.88 0.12 1
hurdles 1 0.94 0.88 0.12 1
run200m 4 0.89 0.79 0.21 1
shot 3 0.84 0.70 0.30 1
highjump 2 0.65 0.43 0.57 1
run800m 7 0.63 0.40 0.60 1
javelin 6 0.50 0.25 0.75 1
PC1
SS loadings 4.32
Proportion Var 0.62
Mean item complexity = 1
Test of the hypothesis that 1 component is sufficient.
The root mean square of the residuals (RMSR) is 0.1
with the empirical chi square 9.81 with prob < 0.78
Fit based upon off diagonal values = 0.97Mit der ersten Komponente werden 62% der Varianz erklärt. Dies lässt sich direkt aus dem Output (Proportion Var) ablesen.
Für den Eigenwert (SS Loadings) von Komponente 1 quadrieren wir die Ladung jeder Variablen auf Komponente 1 und summieren diese Werte anschliessend. Oder wir nehmen die Kommunalitäten h2 und summieren diese! Von Hand lässt sich in diesem Fall auch mit R als “Taschenrechner” interpretieren. In der Prüfung müssen Sie das allerdings mit einem normalen Taschenrechner rechnen…
# Erste Möglichkeit:
# Um genauere Werte zu bekommen, extrahieren wir die Ladungen und
# runden sie auf 4 Stellen hinter dem Komma:
loadings <- round(pca_heptathlon$loadings[, 1], 4)
# Das Quadrieren können wir dann auch vektorisiert machen...
loadings_squared <- loadings^2
# Und wenn wir schon dabei sind...
eigenvalue <- sum(loadings_squared)
eigenvalue |> round(2)[1] 4.32# Und daraus dann...
percent_explained <- eigenvalue / 7 * 100
round(percent_explained, 2)[1] 61.77# Zweite Möglichkeit (wirklich von Hand):
loadings # Ladungen ausgeben hurdles highjump shot run200m longjump javelin run800m
0.9365 0.6540 0.8369 0.8881 0.9377 0.5038 0.6298 pctexplained_vonhand <- (0.9365^2 + 0.654^2 + 0.8369^2 +
0.8881^2 + 0.9377^2 + 0.5038^2 + 0.6298^2) / 7 * 100
round(pctexplained_vonhand, 2)[1] 61.77# Dritte Möglichkeit (über die Kommunalitäten h2):
pctexplained_ueberh2 <- sum(pca_heptathlon$communality) / 7 * 100
round(pctexplained_ueberh2, 2)[1] 61.77Der Eigenwert der ersten Komponente ist also 4.32 und sie erklärt 61.77 % der Gesamtvarianz.
Vertiefung: Aufgabe 4
- Extrahieren Sie die
scoresder ersten Hauptkomponente.
Die scores sind die Werte der Teilnehmer projiziert auf die Achse der entsprechenden Hauptkomponente. Mit anderen Worten: die (standardisierte) Summe der mit den Ladungen gewichteten Werten jeder Teilnehmerin auf den sieben Variablen (Disziplinen).
In unserem Beispiel entsprechen die scores der ersten Hauptkomponente den Werten jeder Person auf unserer eindimensionalen Skala für die Fähigkeit im Heptathlon. Das vorher schon erwähnte Applet visualisiert auch genau solche scores (Komponentenwerte).
Wer ist die Gewinnerin des Heptathlons nach unserer Methode?
Schauen Sie ausserdem, wie stark die
scoresmit der offiziellen Score aus dem Datenframeheptathlonzusammenhängt.
Solution
a) Extraktion der Scores
pc1_scores <- pca_heptathlon$scores[, 1]
pc1_scores [1] 2.288003481 1.513916923 1.407268194 0.619493399 0.723043464
[6] 0.460948428 0.458533216 0.304539034 0.183506561 0.251196423
[11] 0.104697557 0.517465086 -0.001449975 -0.052508930 -0.100449355
[16] -0.111817913 -0.317178053 -0.363988434 -0.904582980 -0.879208346
[21] -1.018691636 -1.332495085 -1.876159087 -1.874081972b) Gewinnerin?
Grundsätzlich müssen wir dafür nur die Person (name) mit dem grössten Wert auf pc1_scores identifizieren. Dafür können wir die Funktion which.max() verwenden.
heptathlon$name[which.max(pc1_scores)][1] "Joyner-Kersee (USA)"Jackie Joyner-Kersee aus den USA hat den höchsten Siebenkampf-Score nach der PCA-Methode.
c) Streudiagramm und Korrelation zwischen den PCA-Scores und den Original-Scores
plot(pc1_scores, heptathlon$score)
cor(pc1_scores, heptathlon$score)[1] 0.9931168So wie es scheint haben die Verantwortlichen für die Bewertung der Teilnehmerinnen eine ähnliche Methode angewendet. Auf jeden Fall kommen wir zu sehr ähnlichen Resultaten, was durch die geringen Abweichungen von der Geraden im Streudiagramm sowie durch die Korrelation nahe 1 veranschaulicht wird.
Eine gute Visualisierung und noch etwas mehr Intuition gibts in diesem Applet.