Forscher:innen der Psychologie oder anderer Natur-, Sozial- und Geisteswissenschaften interessieren sich häufig dafür, wie sich Daten auf einige wenige entscheidende Faktoren herunterbrechen lassen, welche ein theoretisches Erklärungsmodell für die Variation in einem Datensatz liefern. Die Annahme ist hierbei, dass die beobachtbaren Messungen eine Linearkombination (also eine Summe) aus einem systematischen (wahren) und einem unsystematischen (Fehler-)Anteil bilden. Die dahinterliegenden Faktoren sind nicht messbare (latente) Variablen, auf welche, unter gewissen Annahmen, nur anhand der Kovariation zwischen den beobachtbaren Items geschlossen werden kann. Durch diese Zusammenhänge zwischen den Messungen können schließlich Hypothesen für die latenten Variablen untersucht werden. Ein theoriegenerierendes Verfahren, das hierzu häufig verwendet wird, ist die exploratorische Faktorenanalyse (im Folgenden EFA, engl. Exploratory Factor Analysis, vgl. Eid, Gollwitzer & Schmitt, 2017, Kapitel 25. Außerdem können Sie sich Brandt, 2020, Kapitel 23 genauer ansehen, wenn Sie weitere Informationen, bzw. eine zusätzliche Erklärung wünschen).
Wir wollen die EFA zur Auswertung von Beziehungen zwischen Variablen in R
näher kennenlernen. Die EFA ist manchen Gesichtspunkten, sowie im Output in R
, recht verwandt mit der Hauptkomponentenanalyse (PCA) aus dem vergangenen Semester. Allerdings konnten wir bei der PCA kein Erklärungsmodell aufstellen (“wir konnten den Hauptkomponenten nicht so einfach eine inhaltliche Bedeutung zuschreiben”) und auch keine Messfehler (unsystematische Fehleranteile) mitmodellieren. Zudem waren die Hauptkomponenten der PCA Linearkombinationen (also Zusammensetzungen) aus den beobachteten Variablen, bei der EFA hingegen sind die Messungen (beobachteten Variablen) Linearkombinationen aus systematischen (latenten) Variablen sowie Messfehlern.
Bevor wir mit den Analysen beginnen können, laden wir zunächst alle Pakete, welche wir im Folgenden benötigen werden.
library(corrplot) # Korrelationsmatrix grafisch darstellen
library(psych) # EFA durchführen
library(GPArotation) # EFA Lösung rotieren
Datensatz
Wir wollen uns die Faktorenstruktur der Big-5 eines entsprechenden Fragebogens ansehen. Der Originaldatensatz ist ein Onlinedatensatz, wird seit 2012 erfasst und ist auf openpsychometrics.org als .zip downloadbar. Bisher haben über 19700 Proband:innen aus der ganzen Welt daran teilgenommen. Zu jeder der fünf Facetten gibt es 10 Fragen. Der Fragebogen ist auf personality-testing.info einzusehen. Um das Ganze etwas übersichtlicher zu gestalten, betrachten wir einen gekürzten Datensatz. Wir kennen diesen Datensatz, allerdings in anderer Zusammensetzung, bereits aus den Übungen zum vergangenen Semester zum Themenblock MANOVA, weswegen wir hier das Kürzel “EFA” angehängt haben. Im Datensatz Big5_EFA.rda befinden sich 15 Items aus dem Big-5 Persönlichkeitsfragebogen. Hier werden von diesen 10 Items jeweils die ersten drei verwendet. Der Itemwortlaut der verwendeten Items ist wie folgt:
Itemwortlaut | |
Item Nr. | Item |
---|---|
E1 | I am the life of the party. |
E2 | I don't talk a lot. |
E3 | I feel comfortable around people. |
N1 | I get stressed out easily. |
N2 | I am relaxed most of the time. |
N3 | I worry about things. |
A1 | I feel little concern for others. |
A2 | I am interested in people. |
A3 | I insult people. |
C1 | I am always prepared. |
C2 | I leave my belongings around. |
C3 | I pay attention to details. |
O1 | I have a rich vocabulary. |
O2 | I have difficulty understanding abstract ideas. |
O3 | I have a vivid imagination. |
Die Kürzung des vollen Datensatzes lässt sich im Appendix A nachvollziehen. Zusätzlich zu den Persönlichkeitsitems wurden demografische Daten, die mögliche Unterschiede zwischen Personen beschreiben, erfasst.
Daten laden
Wir laden zunächst die Daten: entweder lokal von Ihrem Rechner:
load("C:/Users/Musterfrau/Desktop/Big5_EFA.rda")
oder wir laden sie direkt über die Website:
load(url("https://pandar.netlify.app/daten/Big5_EFA.rda"))
Nun sollte in R
-Studio oben rechts in dem Fenster unter der Rubrik “Data” unser Datensatz mit dem Namen “Big5” erscheinen. Der Datensatz heißt also genauso wie der Datensatz für die MANOVA, enthält aber andere Variablen (also Achtung!).
Übersicht über die Daten
head(Big5, n = 10) # gebe die ersten 10 Zeilen aus
## age engnat gender country E1 E2 E3 N1 N2 N3 A1 A2 A3 C1 C2 C3 O1 O2 O3
## 1 53 1 1 US 4 2 5 1 5 2 1 5 1 4 1 5 4 1 3
## 2 46 1 2 US 2 2 3 2 3 4 1 3 3 4 1 3 3 3 3
## 3 14 2 2 PK 5 1 1 5 1 5 5 1 5 4 1 5 4 5 5
## 4 19 2 2 RO 2 5 2 5 4 4 2 5 4 3 3 4 4 3 5
## 5 25 2 2 US 3 1 3 3 3 3 5 5 3 3 1 5 3 1 1
## 6 31 1 2 US 1 5 2 1 5 4 2 2 3 2 5 4 4 2 1
## 7 20 1 2 US 5 1 5 2 4 2 5 5 1 2 4 3 3 1 5
## 8 23 2 1 IN 4 3 5 1 4 4 2 5 1 4 2 5 3 1 5
## 9 39 1 2 US 3 1 5 2 4 5 1 5 1 4 3 5 3 3 5
## 10 18 1 2 US 1 4 2 5 2 5 2 3 1 5 2 4 4 2 5
Wir sehen, dass in den ersten 4 Spalten die demografischen Daten wie etwa Alter (“age”), Englisch als Muttersprache (“engant”, 1=yes, 2=no, 0=missed), Geschlecht (“gender”, 1=Male, 2=Female, 3=Other, 0=missed) und Herkunftsland (“country”, ISO-kodiert, bspw. “DE” = Deutschland, “FR” = Frankreich, “EM” = Vereinigte Arabische Emirate, “US” = Vereinigten Staaten von Amerika) eingetragen wurden. In den darauf folgenden Spalten sind die Items der Extraversion (engl. extraversion, Items: E1, E2, E3), des Neurotizismus (engl. neuroticism, Items: N1, N2, N3), der Verträglichkeit (engl. agreeableness, Items: A1, A2, A3), der Gewissenhaftigkeit (engl. conscientiousness, Items: C1, C2, C3) und der Offenheit für Erfahrungen (engl. openness, Items: O1, O2, O3) eingetragen. Beispielsweise ist die erste Person des Datensatzes ein 53-jähriger Mann, der Englisch als Muttersprache spricht und in den USA lebt.
Da wir uns in der Praxis nur sehr selten in der glücklichen Lage befinden, einen solch riesigen Datensatz zu haben, wollen wir uns innerhalb des Datensatzes auf Subgruppen beschränken: wir wollen uns zunächst nur Daten von Personen aus Frankreich ansehen. Dazu wählen wir nur diejenigen Zeilen aus, in denen country == "FR"
gilt. Das erreichen wir wie folgt: Mit Big5$country
haben wir Zugriff auf die “Country”-Spalte im Datensatz und können mit == "FR"
prüfen, an welchen Stellen hier “FR” steht, also Personen, die in Frankreich leben. dim
gibt die Dimensionen des Datensatzes wieder.
dim(Big5)
## [1] 19711 19
data_France <- Big5[Big5$country == "FR", ]
dim(data_France)
## [1] 129 19
Dem Output sollte zu entnehmen sein, dass data_France
129 Zeilen (also Proband:innnen, die in Frankreich leben) und 19 Spalten (also Variablen) enthält. Für die weiteren Analysen brauchen wir die demografischen Variablen in dem Datensatz der in Frankreich lebenden Teilnehmer:innen nicht mehr. Aus diesem Grund speichern wir den Datensatz noch einmal ohne die ersten 4 Spalten ab. Anschließend stellen wir die Korrelationsmatrix dieser Daten grafisch dar (den Befehl corrplot
aus dem gleichnamigen Paket kennen wir bereits aus der PCA Sitzung des vergangenen Semesters).
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
#### Visualisierte Korrelationsmatrix in dataFR
corrplot(corr = cor(dataFR), # Korrelationsmatrix (Datengrundlage)
method = "color", # zeichne die Ausprägung der Korrelation farblich kodiert
addCoef.col = "black", # schreibe die Korrelationskoeffizienten in schwarz in die Grafik
number.cex = 0.7) # stelle die Schriftgröße der Koeffizienten ein

Auf den ersten Blick scheinen die Items der gleichen Skala (ausgedrückt durch gleiche Buchstaben pro Item) stärker (betragsmäßig höher) miteinander zu korrelieren. Allerdings sind hier sehr viele Korrelationen abgetragen. Wir wollen uns zunächst nur auf Extraversion und Neurotizismus beschränken.
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
head(dataFR2)
## E1 E2 E3 N1 N2 N3
## 17 1 3 2 4 2 3
## 398 3 3 3 4 3 4
## 488 1 5 2 4 4 5
## 545 1 2 1 4 1 5
## 551 1 4 1 5 1 5
## 656 1 4 2 5 1 5
# zum gleichen Ergebnis würde auch Folgendes kommen (besonders von Relevanz,
# wenn wir bspw. nicht die Position sondern nur die Namen der Variablen kennen!):
head(dataFR[, c("E1", "E2", "E3", "N1", "N2", "N3")])
## E1 E2 E3 N1 N2 N3
## 17 1 3 2 4 2 3
## 398 3 3 3 4 3 4
## 488 1 5 2 4 4 5
## 545 1 2 1 4 1 5
## 551 1 4 1 5 1 5
## 656 1 4 2 5 1 5
Wenn wir uns die Korrelationmatrix des gekürzten Datensatzes dataFR2
ansehen…
# Visualisierte Korrelationsmatrix
corrplot(corr = cor(dataFR2), # Korrelationsmatrix (Datengrundlage)
method = "color", # Zeichne die Ausprägung der Korrelation farblich kodiert
addCoef.col = "black", # schreibe die Korrelationskoeffizienten in schwarz in die Grafik
number.cex = 1) # Stelle die Schriftgröße der Koeffizienten ein

…erkennen wir deutlich, dass die Extraversionsitems und die Neurotizismusitems untereinander jeweils stärker zusammenhängen als zwischen den Konstrukten. Dennoch ist der Grafik zu entnehmen, dass die beiden Konstrukte nicht unabhängig voneinander sind (es gibt Beziehungen zwischen Items der beiden Konstrukte).
Ziel: EFA
Unser Ziel ist es, mit den gegebenen Items eine exploratorische Faktorenanalyse durchzuführen. Wir wollen hierbei die Anzahl der Faktoren mittels einer Parallelanalyse bestimmen und anschließend dieses Modell mit dem

Auf unseren Datensatz angepasst, wollen wir folgendes Modell anpassen:

Natürlich erwarten wir, dass insgesamt 2 Faktoren die Daten am besten beschreiben und dass die konstruktkongruenten Items jeweils auf dem gleichen Faktor am stärksten laden. Aber stützen die Daten diese Hypothese?
Wir wollen im Folgenden
- eine Parallelanalyse durchführen, um in Erfahrung zu bringen, wie viele Faktoren sinnvoll zu den Daten passen
- eine Hauptachsenanalyse mit orthogonaler und obliquer Rotation durchführen
- eine exploratorische Maximum-Likelihood-Faktorenanalyse durchführen und die Passung zu den Daten untersuchen
- im Rahmen der exploratorischen Maximum-Likelihood-Faktorenanalyse die Passung zu den Daten im Vergleich zu konkurrierenden Modellen untersuchen.
Parallelanalyse und Auswahl an Faktoren
Zur Auswahl der Anzahl an Faktoren in der EFA kann auf die Eigenwerte zurückgegriffen werden. Diese Eigenwerte entstehen beispielsweise durch Lösen des Eigenwerteproblems und entsprechen den Varianzen der Faktoren. Hier gilt es, nur solche Faktoren zu wählen, die auch große Varianzen haben. Die Parallelanalyse hatten wir im Zusammenhang mit der Hauptkomponentenanalyse kennengelernt. Hier werden vielfach (z.B. 1000 Mal, für Vergleichbarkeit sogar besser mehr!) unabhängige Daten in dem gleichen Format des ursprünglichen Datensatzes gezogen und eine PCA oder EFA durchgeführt. Die entstehenden Eigenwerte werden der Größe nach sortiert und dann über die Wiederholungen gemittelt. So entsteht ein auf die Stichprobe und Anzahl der Variablen genormter, zufälliger, durchschnittlicher Eigenwerteverlauf. Sind Eigenwerte der tatsächlich beobachteten Daten größer als die der Parallelanalyse, so spricht dies für eine/n bedeutsame/n Komponente/Faktor. Weitere Kriterien zur Auswahl von zu extrahierenden Faktoren im Rahmen der PCA waren das Eigenwerte-größer-1 Kriterium (Kaiser-Guttman-Kriterium) sowie der Scree-Test (Elbow-Criterion, Knick im Eigenwerteverlauf). Weitere Informationen zur EFA sowie zu Wiederholungen der PCA und der Auswahlkriterien können beispielsweise in Eid, Gollwitzer und Schmitt (2017) in Kapitel 25 (Seite 919 und folgend) nachgelesen werden.
Der wesentliche Unterschied zwischen einer EFA und einer PCA ist, dass bei der EFA angenommen wird, dass die beobachteten Variablen systematische (wahre) und unsystematische (Fehler-) Anteile enthalten. Es wird somit ein Erklärungsmodell, welches die Variation zwischen den Variablen erzeugt, postuliert. Bei der PCA werden die beobachteten Variablen als messfehlerfrei angenommen. Eine wichtige Folge aus der (Nicht-) Modellierung der Fehler ist, dass in der Regel die Faktorladungen bei der PCA höher ausfallen als bei der EFA. Dies liegt daran, dass bei der PCA die Variablen mit ihren eigenen Messfehlern, aus welchen auch die Hauptkomponenten unter anderem zusammengesetzt sind, korrelieren. Die Faktorladungen/ Komponentenladungen stehen hierbei (im orthogonalen Fall) für die Korrelation zwischen Item und Faktor/ Komponente; im obliquen Fall sind die Faktorladungen als Regressionkoeffizienten zu interpretieren. Wird eine ML-EFA an die Daten angepasst, so wird zusätzlich noch ein Erklärungsmodell basierend auf Verteilungsannahmen (multivariate Normalverteilung der Faktoren und Fehler — und als Konsequenz daraus auch der Items) herangezogen. Bei der PCA sind die Hauptkomponenten lediglich Linearkombinationen aus den beobachteten Variablen ohne jegliche Verteilungsannahmen (die Hauptkomponenten bestehen aus gewichteten Summen der beobachteten Variablen). Entsprechend ist es bei der PCA auch nicht möglich, konkurierende Modelle gegeneinander zu testen. Es bleiben dort nur die deskriptiven und recht subjektiven Auswahlkriterien für die Anzahl der Komponenten.
Mit Hilfe des fa.parallel
-Befehls aus dem psych
-Paket, welchen wir im Rahmen der PCA bereits kennengelernt hatten, lässt sich ganz einfach der Eigenwerteverlauf inklusive Parallelanalyse grafisch darstellen.
fa.parallel(dataFR2)

## Parallel analysis suggests that the number of factors = 2 and the number of components = 2
Ohne weitere Einstellungen wird der Eigenwerteverlauf der PCA und der EFA ausgegeben. Deutlich zu sehen ist, dass die Eigenwerte der PCA größer ausfallen als die der EFA. Dies liegt erneut daran, dass die Faktoren der EFA lediglich die systematischen Anteile der Variablen enthalten, während die Komponenten der PCA Kompositionen sind - also Zusammensetzungen aus den Variablen; inklusive der Messfehler. Wählen wir fa = "fa"
, so wird uns nur der Verlauf der Eigenwerte auf Basis einer EFA aufgeführt.
fa.parallel(dataFR2, fa = "fa")

## Parallel analysis suggests that the number of factors = 2 and the number of components = NA
Die Grafik zeigt drei Eigenwerteverläufe. FA Actual Data ist der Eigenwerteverlauf unseres Datensatzes. FA Simulated Data ist der Eigenwerteverlauf basierend auf den 1000 simulierten Datensätzen. FA Resampled Data ist der Eigenwerteverlauf von Datensätzen, der durch Resampling, also neues Verteilen unseres Datensatzes entsteht (das ist im Grunde Bootstrapping, was im Rahmen der Sitzung zu SEM näher erläutert wird).
Der Parallelanalyse der EFA ist zu entnehmen, dass voraussichtlich 2 Faktoren genügen, um die Variation im Datensatz zu erklären. Auch die Parallelanalyse der PCA (Grafik zuvor) lässt dies vermuten. Des Weiteren sprechen beide Scree-Tests für einen Knick um den 3. Faktor/die 3. Komponente, was auch für eine Dimensionalität von 2 spricht. Zu guter Letzt zeigt auch das Kaiser-Guttman-Kriterium kein anderes Ergebnis. Allerdings ist dieses Kriterium nur sinnvoll auf den Eigenwerteverlauf der PCA anwendbar, weswegen wir es auch nur im Bezug auf den PCA-Eigenwerteverlauf interpretieren.
Orthogonale und oblique Hauptachsenanalyse
Da unsere Auswahlkriterien einstimmig für 2 Faktoren sprechen und dies auch unsere Hypothese war, modellieren wir zunächst eine orthogonale Hauptachsenanalyse. Dazu nutzen wir den fa
(Factor Analysis) Befehl des psych
-Paketes. Mit Hilfe der Argumente nfactors
und rotate
lässt sich die Anzahl an Faktoren sowie die Rotation auswählen (genauso wie bei der pca
Funktion für die PCA!). Wir wollen hier orthogonal varianzmaximierend (also varimax) rotieren.
fa(dataFR2, nfactors = 2, rotate = "varimax")
## Factor Analysis using method = minres
## Call: fa(r = dataFR2, nfactors = 2, rotate = "varimax")
## Standardized loadings (pattern matrix) based upon correlation matrix
## MR1 MR2 h2 u2 com
## E1 0.69 -0.06 0.48 0.52 1.0
## E2 -0.65 0.00 0.42 0.58 1.0
## E3 0.82 -0.19 0.70 0.30 1.1
## N1 -0.01 0.84 0.70 0.30 1.0
## N2 0.10 -0.58 0.35 0.65 1.1
## N3 -0.05 0.59 0.34 0.66 1.0
##
## MR1 MR2
## SS loadings 1.58 1.42
## Proportion Var 0.26 0.24
## Cumulative Var 0.26 0.50
## Proportion Explained 0.53 0.47
## Cumulative Proportion 0.53 1.00
##
## Mean item complexity = 1
## Test of the hypothesis that 2 factors are sufficient.
##
## df null model = 15 with the objective function = 1.47 with Chi Square = 183.82
## df of the model are 4 and the objective function was 0.07
##
## The root mean square of the residuals (RMSR) is 0.04
## The df corrected root mean square of the residuals is 0.07
##
## The harmonic n.obs is 129 with the empirical chi square 4.98 with prob < 0.29
## The total n.obs was 129 with Likelihood Chi Square = 9.06 with prob < 0.06
##
## Tucker Lewis Index of factoring reliability = 0.886
## RMSEA index = 0.099 and the 90 % confidence intervals are 0 0.187
## BIC = -10.38
## Fit based upon off diagonal values = 0.99
## Measures of factor score adequacy
## MR1 MR2
## Correlation of (regression) scores with factors 0.89 0.88
## Multiple R square of scores with factors 0.79 0.77
## Minimum correlation of possible factor scores 0.58 0.55
Im Output ganz oben erkennen wir die Schätzmethode (hier: minres
, also Minimierung der Residuen). Aus diesem Grund heißen die Faktoren in diesem Output auch MR1 und MR2; für Minimale-Residuen-Faktor 1 und 2 (diese Benennung war uns auch bereits bei der PCA aufgefallen, wo diese je nach Rotation entweder PC1 oder RC1, etc., hießen). Die Faktorladungen zu den zugehörigen Faktoren sind unter Standardized loadings (pattern matrix) based upon correlation matrix
zu sehen. h2
steht für die Kommunalität (u2
ist die “uniqueness” (SS loadings
steht für “Sum of Squares loadings”, also die Quadratsumme der Faktorladungen. Diese ist gleich dem Eigenwert: Proportion Var
betitelt den Anteil der Variation, der durch die jeweiligen Faktoren erklärt werden kann. Cumulative Var
kumuliert, also summiert, diese Anteile bis zum jeweiligen Faktor auf ($\text{CumVar}j = \sum{k=1}^j\theta_k = \theta_1+\dots+\theta_jProportion Explained
setzt die Variation, die durch die Faktoren erklärt wird, in Relation zur gesamten erklärten Varianz (d.h. hier summiert sich die erklärte Varianz immer zu 1, während sich die proportionale Varianz nur zu 1 aufsummiert, wenn die gesamte Variation im Datensatz auf die beiden Variablen zurückzuführen ist). Cumulative Proportion
beschreibt das gleiche wie Cumulative Var
, nur bezieht sie sich hier auf die Proportion Explained
. Bei der Interpretation dieser Kennwerte ist zu bedenken, dass bei der EFA angenommen wird, dass die beobachteten Variablen Messfehler enthalten (also die Reliabilität nicht als 1 angenommen werden kann). Folglich ist die Kommunalität
Außerdem werden durch diesen Befehl sehr viele Informationen ausgegeben. Aus diesen Gründen speichern wir uns diese Analyse als ein Objekt ab, welchem wir dann gezielt Informationen mit Hilfe von ...$...
entlocken können. Welche Argumente entlockt werden können, kann beispielsweise mit names
herausgefunden werden. Wir speichern das Objekt unter dem Namen two_factor ab.
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
names(two_factor) # mögliche Informationen
## [1] "residual" "dof" "chi" "nh" "rms"
## [6] "EPVAL" "crms" "EBIC" "ESABIC" "fit"
## [11] "fit.off" "sd" "factors" "complexity" "n.obs"
## [16] "objective" "criteria" "STATISTIC" "PVAL" "Call"
## [21] "null.model" "null.dof" "null.chisq" "TLI" "CFI"
## [26] "RMSEA" "BIC" "SABIC" "r.scores" "R2"
## [31] "valid" "score.cor" "weights" "rotation" "hyperplane"
## [36] "communality" "communalities" "uniquenesses" "values" "e.values"
## [41] "loadings" "model" "fm" "rot.mat" "Structure"
## [46] "method" "scores" "R2.scores" "r" "np.obs"
## [51] "fn" "Vaccounted" "ECV"
Beispielsweise erhalten wir mit $loadings
die Faktorladungsmatrix sowie Informationen über die Eigenwerte. (Wichtig für später: Die richtigen Kommunalitäten werden mit $communality
angefordert.)
two_factor$loadings
##
## Loadings:
## MR1 MR2
## E1 0.692
## E2 -0.646
## E3 0.819 -0.186
## N1 0.837
## N2 -0.580
## N3 0.585
##
## MR1 MR2
## SS loadings 1.578 1.419
## Proportion Var 0.263 0.236
## Cumulative Var 0.263 0.499
Hier ist relativ deutlich die Zuordnung zu den jeweiligen Faktoren zu sehen. Faktor 1 (MR1) entspräche post-hoc interpretiert (die Theorie wird also aus den Daten generiert; es sind auch andere Interpretationsansätze zulässig) der Extraversion, während der zweite Faktor (MR2) dem Neurotizismus entspräche. Indem wir hinter loadings
eckige Klammern mit einem Komma setzten ([,]
), bekommen wir alle Nachkommastellen ohne Runden angzeigt (hätten wir bspw. round
auf die Ladungsmatrix oben angewendet, würden sich auch die Eigenwerte ändern, weswegen diese Ansicht sich nicht wirklich zur genauen Interpretation der Eigenwerte eignet! Hier hatten wir das Argument Vaccounted
bereits im Rahmen der PCA kennengelernt):
two_factor$loadings[,] # ohne seltsames Runden
## MR1 MR2
## E1 0.69189323 -0.058067393
## E2 -0.64577474 0.004122088
## E3 0.81858640 -0.185912523
## N1 -0.01345489 0.837407770
## N2 0.09879742 -0.580348020
## N3 -0.04734674 0.585385819
Die Faktorladungsmatrix wird auch manchmal Mustermatrix genannt. Die Strukturmatrix enthält Informationen über die Korrelation der Items mit den jeweiligen Faktoren. Sie heißt Structure
. Auch hier ist der Zusatz [,]
sinnvoll!
two_factor$Structure[,]
## MR1 MR2
## E1 0.69189323 -0.058067393
## E2 -0.64577474 0.004122088
## E3 0.81858640 -0.185912523
## N1 -0.01345489 0.837407770
## N2 0.09879742 -0.580348020
## N3 -0.04734674 0.585385819
Wir sehen deutlich, dass die Strukturmatrix sich nicht von der Faktorladungsmatrix unterscheidet. Das liegt daran, dass die Faktoren noch als unkorreliert angenommen werden: somit ist die Faktorladungsmatrix gleich der Strukturmatrix. Genau aus diesem Grund hatten wir uns die Strukturmatrix auch nicht im Rahmen der PCA angesehen — dort hatten wir nur orthogonal rotiert.
Da wir nicht davon ausgehen können, dass die Faktoren unkorreliert sind, wollen wir die gleiche Analyse nun für oblique (“oblimin” in R
) rotierte Faktoren durchführen.
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")
Die einzig neue Information können wir unter With factor correlations of
ablesen: die Korrelation zwischen den Faktoren. Im Output der obliquen Rotation ist zu erkennen, dass sich die Kommunalitäten nicht ändern. Wir können also nicht mehr Variation im Datensatz erklären. Die Varianz wird nur umverteilt, wie den veränderten Eigenwerten neben SS loadings
zu entnehmen ist. Hier hat der erste Faktor einen etwas größeren Eigenwert als im orthogonalen Fall (entsprechend ist der 2. Eigenwert kleiner, da nicht mehr Variation erklärt wird):
two_factor$Vaccounted
## MR1 MR2
## SS loadings 1.5780086 1.4186844
## Proportion Var 0.2630014 0.2364474
## Cumulative Var 0.2630014 0.4994488
## Proportion Explained 0.5265833 0.4734167
## Cumulative Proportion 0.5265833 1.0000000
two_factor_oblimin$Vaccounted
## MR1 MR2
## SS loadings 1.6065736 1.3901195
## Proportion Var 0.2677623 0.2316866
## Cumulative Var 0.2677623 0.4994488
## Proportion Explained 0.5361155 0.4638845
## Cumulative Proportion 0.5361155 1.0000000
Schauen wir uns die Ladungsmatrix an, …
two_factor_oblimin$loadings[,] # Ladungsmatrix
## MR1 MR2
## E1 0.70055173 0.03946625
## E2 -0.65623412 -0.08729631
## E3 0.82326401 -0.07142273
## N1 0.02609114 0.84195602
## N2 0.07286715 -0.57082957
## N3 -0.02032830 0.58319228
… so ist post-hoc interpretiert anzunehmen, dass der erste Faktor die Extraversion abbildet und der zweite den Neurotizismus.
Entlocken Sie doch mal dem Objekt two_factor_oblimin
die latente Kovarianzmatrix, also die Kovarianzmatrix der latenten Variablen. Tipp, der griechische Buchstabe in diesem Zusammenhang ist häufig
## MR1 MR2
## MR1 1.0000000 -0.1852246
## MR2 -0.1852246 1.0000000
Als neue Information entnehmen wir der Korrelationsmatrix der Faktoren, dass die beiden Faktoren negativ korreliert sind zu -0.19. Nun wollen wir nachschauen, ob sich tatsächlich die Strukurmatrix im oblique-rotierten Fall von der Faktorladungsmatrix unterscheidet:
two_factor_oblimin$loadings[,]
## MR1 MR2
## E1 0.70055173 0.03946625
## E2 -0.65623412 -0.08729631
## E3 0.82326401 -0.07142273
## N1 0.02609114 0.84195602
## N2 0.07286715 -0.57082957
## N3 -0.02032830 0.58319228
two_factor_oblimin$Structure[,]
## MR1 MR2
## E1 0.6932416 -0.09029319
## E2 -0.6400647 0.03425442
## E3 0.8364933 -0.22391151
## N1 -0.1298599 0.83712330
## N2 0.1785989 -0.58432636
## N3 -0.1283499 0.58695759
Sie sehen, dass sich nun die Strukturmatrix von der Faktorladungsmatrix unterscheidet. Die Unterschiede sind allerdings nicht sehr groß, da die Korrelation zwischen den beiden Faktoren mit -0.19 betragsmäßig nicht sonderlich groß ausfällt. Weitere Informationen und wie die beiden Matrizen ineinander überführbar sind, erfahren Sie im Appendix C. In Appendix D erfahren Sie, wie Sie Kommunalitäten und Eigenwerte im Rahmen der EFA nur mit Hilfe der loadings
und der Faktorkorrelation Phi
bestimmen. Falls Sie sich nun zu Recht wundern, dass Appendix B fehlt, dann sei gesagt, dass dieser alle Analysen am vollständigen Datensatz durchführt und es durch die Abhängigkeiten der Daten in den Appendizes Sinn macht, diesen zwar im Text später zu erwähnen aber inhaltlich vorzuziehen.
Die Frage ist nun, ob unser Modell überhaupt zu den Daten passt.
Exploratorische Maximum-Likelihood-Faktorenanalyse (ML-EFA)
Wir möchten unsere Analysen nun gegen andere konkurrierende Modelle absichern sowie untersuchen, ob unser zweifaktorielles Modell überhaupt zu den Daten passt. Hierzu müssen wir annehmen, dass unsere Daten multivariat normalverteilt sind. Wie man diese Annahme zumindest deskriptiv untersucht, hatten wir im Zusammenhang mit den Voraussetzungen von statistischen Verfahren kennengelernt (Mahalanobisdistanz sollte approximativ
Im unkorrelierten/orthogonalen Fall sähe die Matrix so aus:
Unter der Normalverteilungsannahme brauchen wir nur (Ko-)Varianzen und Mittelwerte, um unsere Variablen vollständig zu beschreiben. In der Schätzung der ML-EFA geht es uns darum, die Parameter in
Um mit Hilfe von fa
eine ML-EFA durchzuführen, muss dem Argument fm
die entsprechende Bezeichnung "ml"
übergeben werden.
two_factor_ML <- fa(dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")
two_factor_ML
## Factor Analysis using method = ml
## Call: fa(r = dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
## ML1 ML2 h2 u2 com
## E1 0.69 0.04 0.47 0.53 1.0
## E2 -0.65 -0.06 0.41 0.59 1.0
## E3 0.83 -0.06 0.71 0.29 1.0
## N1 0.03 0.84 0.70 0.30 1.0
## N2 0.10 -0.57 0.35 0.65 1.1
## N3 -0.01 0.59 0.36 0.64 1.0
##
## ML1 ML2
## SS loadings 1.61 1.39
## Proportion Var 0.27 0.23
## Cumulative Var 0.27 0.50
## Proportion Explained 0.54 0.46
## Cumulative Proportion 0.54 1.00
##
## With factor correlations of
## ML1 ML2
## ML1 1.00 -0.18
## ML2 -0.18 1.00
##
## Mean item complexity = 1
## Test of the hypothesis that 2 factors are sufficient.
##
## df null model = 15 with the objective function = 1.47 with Chi Square = 183.82
## df of the model are 4 and the objective function was 0.07
##
## The root mean square of the residuals (RMSR) is 0.04
## The df corrected root mean square of the residuals is 0.07
##
## The harmonic n.obs is 129 with the empirical chi square 5.6 with prob < 0.23
## The total n.obs was 129 with Likelihood Chi Square = 8.75 with prob < 0.068
##
## Tucker Lewis Index of factoring reliability = 0.893
## RMSEA index = 0.096 and the 90 % confidence intervals are 0 0.184
## BIC = -10.69
## Fit based upon off diagonal values = 0.99
## Measures of factor score adequacy
## ML1 ML2
## Correlation of (regression) scores with factors 0.90 0.88
## Multiple R square of scores with factors 0.80 0.77
## Minimum correlation of possible factor scores 0.61 0.55
Wir sehen, dass diesmal die Schätzmethode “ml” ist. Auch die Faktoren heißen nun ML1 und ML2. Die Faktorladungen im ML-EFA Modell mit obliquer Rotation sehen den Faktorladungen aus unserer vorigen Analyse sehr ähnlich.
Uns interessiert nun die Modellpassung, also inwiefern unsere Daten von unserem behaupteten Modell abweichen. $STATISTIC
und $PVAL
entlocken der Analyse (abgespeichert als Objekt) den $chi
, was den empirisch und nicht likelihoodbasierten
two_factor_ML$STATISTIC # Likelihood basierter Chi²-Wert
## [1] 8.749298
two_factor_ML$PVAL # p-Wert
## [1] 0.06768059
Dem ist zu entnehmen, dass auf dem Signifikanzniveau von 5% die Hypothese auf Passung der Kovarianz unserer Daten mit der modellimplizierten Kovarianz in der Population nicht verworfen wird. Die Daten widersprechen dem zweifaktoriellen Modell nicht (die untersuchte Null-Hypothese ist:
Modellvergleich: ML-EFA
Das einfaktorielle Modell erhalten wir ganz einfach via:
one_factor_ML <- fa(dataFR2, nfactors = 1, rotate = "oblimin", fm = "ml")
one_factor_ML$STATISTIC # Chi²-Wert
## [1] 76.76935
one_factor_ML$PVAL # p-Wert
## [1] 7.063217e-13
Das einfaktorielle Modell scheint nicht zu den Daten zu passen (Mit einer Irrtumswahrscheinlichkeit von 5% ist davon auszugehen, dass in der Population die Differenz zwischen der Populationskovarianzmatrix und der modellimplizierten Kovarianzmatrix, bzw. der daraus folgenden Likelihoods, nicht 0 ist.). Dennoch wollen wir dies genau wissen und vergleichen die beiden Modelle direkt miteinander. Für einen solchen Vergleich ist es notwendig, dass es sich bei den beiden Modellen um geschachtelte Modelle handelt. Das bedeutet, dass ein Modell durch Restriktionen von Modellparametern aus dem anderen Modell erzeugt werden kann. Das einfaktorielle Modell lässt sich aus dem zweifaktoriellen Modell durch die Restriktion gewinnen, dass alle Ladungen auf dem Faktor 0 sind und die Varianz dieses Faktors dementsprechend ebenfalls 0 ist.
Mit Hilfe des anova
-Befehls, welchen wir schon bei einigen anderen Modellvergleichen im Rahmen der Regression, der logistischen Regression sowie der Multi-Level Modelle kennengelernt haben, lässt sich nun das einfaktorielle mit dem zweifaktoriellen Modell vergleichen.
anova(one_factor_ML, two_factor_ML)
## Model 1 = fa(r = dataFR2, nfactors = 1, rotate = "oblimin", fm = "ml")
## Model 2 = fa(r = dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")
df | d.df | chiSq | d.chiSq | PR | test | empirical | d.empirical | test.echi | BIC | d.BIC | RMSEA | rms | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
1 | 9 | NA | 76.77 | NA | NA | NA | 142.43 | NA | NA | 33.03 | NA | 0.24 | 0.19 |
2 | 4 | 5 | 8.75 | 68.02 | 0 | 13.6 | 5.60 | 136.83 | 27.37 | -10.69 | -68.02 | 0.10 | 0.04 |
Zunächst bekommen wir angezeigt, welche Modelle unter welchem Kürzel gegeneinander getestet werden. Die Modellnummer steht im ANOVA-Output entsprechend für das jeweilige Modell. In der Zeile 2
steht also der Output für das 2-faktorielle Modelle und am Ende dieser Zeile steht auch der Modellvergleich.
In der ersten Spalte stehen beispielsweise die Freiheitsgrade der Modelle (df
). Daneben steht die Differenz der Freiheitsgrade (d.df
) und dahinter stehen verschiedene Fit-Maße, bzw. Modellvergleiche. Wir müssen im Bereich des chiSq und nicht bei empirical (eine Näherung des chiSq
oder empirical
) angezeigt wird und anschließend die Differenz (d.chiSq
oder d.empirical
) sowie die Signifikanzentscheidung berichtet werden (PR test
, welche allerdings nicht bei der Näherung angezeigt wird). Außerdem wird noch eine Rubrik test
, bzw. test.echi
angezeigt, welche noch einmal die d.df
(häufig BIC
, sowie dessen Differenz d.BIC
aufgeführt. Dieses Informationskriterium werden wir bei der CFA und deren Modellpassung noch einmal genauer betrachten.
Insgesamt wird die Null-Hypothese, dass beide Modell die Daten gleich gut beschreiben, verworfen. Wir entscheiden uns — Ockhams Rasiermesser folgend (siehe Eid et al., 2017, p. 787) — somit für das Modell mit mehr Parametern - das weniger restriktive Modell - welches die Daten besser beschreibt: hier das zweifaktorielle Modell. Nun ist die Frage, ob wir das Modell noch weiter verbessern können, indem wir drei anstatt zwei Faktoren verwenden, um die Kovariation zwischen den Variablen zu beschreiben.
# Passt auch eines mit 3 Faktor?
three_factor_ML <- fa(dataFR2, nfactors = 3, rotate = "oblimin", fm = "ml")
three_factor_ML$STATISTIC # Chi²-Wert
## [1] 0.0328432
three_factor_ML$PVAL # p-Wert
## [1] NA
Das dreifaktorielle Modell beschreibt die Daten perfekt. Das liegt daran, dass es im dreifaktoriellen Modell genauso viele Parameter gibt, wie es empirische Informationen im Datensatz gibt. Demnach lässt sich die empirische Korrelationsmatrix perfekt durch die modelltheoretische Korrelationsmatrix (diejenige Korrelationsmatrix, die sich ergibt, wenn nur die Beziehungen zwischen den Variablen bestehen, die durch das Modell angenommen werden) darstellen. Ein Test auf Modellpassung ist in diesem Fall nicht möglich und auch nicht nötig (deshalb wird beim $PVAL
nichts bzw. NA
ausgegeben). Nun vergleichen wir die beiden Modelle (es ist sehr sinnvoll, das komplexere Modell rechts hin zuschreiben, da es ansonsten bei manchen Analysemethoden zu negativen
anova(two_factor_ML, three_factor_ML)
## Warning in data.frame(df = unlist(dfs), d.df = c(NA, delta.df), chiSq = unlist(chi), : row names
## were found from a short variable and have been discarded
## Model 1 = fa(r = dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")
## Model 2 = fa(r = dataFR2, nfactors = 3, rotate = "oblimin", fm = "ml")
## Warning in data.frame(df = unlist(dfs), d.df = c(NA, delta.df), chiSq = unlist(chi), : row names
## were found from a short variable and have been discarded
df | d.df | chiSq | d.chiSq | PR | test | empirical | d.empirical | test.echi | BIC | d.BIC | RMSEA | rms | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
1 | 4 | NA | 8.75 | NA | NA | NA | 5.60 | NA | NA | -10.69 | NA | 0.1 | 0.04 |
2 | 0 | 4 | 0.03 | 8.72 | 0.07 | 2.18 | 0.02 | 5.58 | 1.4 | -10.69 | -8.72 | 0.1 | 0.00 |
Der d.df
liegt bei 4 (anova
-Befehl, steht hier das komplexere Modell links), nicht verworfen. Aus diesem Grund entscheiden wir uns - Ockhams Rasiermesser folgend (siehe Eid et al., 2017, p. 787) - für das sparsamere Modell, also jenes, welches weniger Parameter enthält und somit restriktiver ist, hier: das zweifaktorielle Modell. Denn in der Wissenschaft streben wir danach, Modelle möglichst einfach zu halten!
Wenn Sie davon noch nicht genug haben, so können Sie in Appendix B nachlesen, wie eine EFA am gesamten (gekürzten) Datensatz durchgeführt wird. Dort stehen keine neuen Informationen zur Durchführung, Sie sollten es folglich eher als Übung ansehen, falls Sie sich entschließen, Appendix B durchzuarbeiten!
Die Zuordnung, die wir hier gefunden haben, entspringt der spezifischen Stichprobe, die wir untersucht haben. Wenn wir a priori aufgestellte Theorien über die Faktorstruktur prüfen wollen, können wir uns der konfirmatorischen Faktorenanalyse bedienen, die wir in der nächsten Sitzung betrachten.
Appendix A
Kürzen des Datensatzes
Falls Sie die Originaldaten auf openpsychometrics.org als .zip herunterladen wollen, so können Sie diesen auf die hier verwendeten Daten wie folgt kürzen:
data_full <- read.table("BIG5/data.csv", header = T, sep = "\t") # nach entpacken des .zip liegen die Daten in einem Ordner namens Big5
### Entferne leere Zeilen und Zeilen mit Missings aus dem Datensatz
ind <- apply(data_full, 1, FUN = function(x) any(is.na(x))) # erzeuge eine Variable, welche TRUE ist, wenn mindestens ein Eintrag pro Zeile fehlt und ansonsten FALSE anzeigt
data_full <- data_full[!ind, ] # Wähle nur diejenigen Zeilen, in denen unsere Indikatorvariable "ind" NICHT TRUE anzeigt, also wo alle Einträge vorhanden sind
# !ind (Ausrufezeichen vor ind) negiert die Einträge in ind (Prüfe bspw. !FALSE == TRUE, nicht false ist gleich true)
### Shorten Data Set
Big5 <- data_full[, c(2:4,7,7+rep(1:3,5)+sort(rep(seq(0,40,10),3)))]
# Verwende nur 3 Items pro Skala plus einige demografische Items
Big5 <- data.frame(Big5) # Schreibe Datensatz als data.frame
save(list = c("Big5"), file = "Big5.rda")
# Speichere gekürzten Datensatz in .rda file (dem R-internen Datenformat)
## --> Das ist auch der Datensatz, den wir weiter verwendet haben!
Appendix B
ML-EFA für den gesamten (gekürzten) Datensatz
Für den vollen Datensatz mit jeweils drei Items pro Persönlichkeitsfacette, nehmen wir zunächst an, dass es 5 Faktoren gibt. Dies wird hier allerdings nicht durch die Parallelanalyse gestützt. Wir müssen die Funktion fa.parallel
diesmal auf den vollen (gekürzten) Datensatz anwenden; nämlich auf dataFR
.
fa.parallel(x = dataFR,fa = "fa")

## Parallel analysis suggests that the number of factors = 4 and the number of components = NA
Hier scheinen eher 4 Faktoren sinnvoll. Wir prüfen dennoch erstmal unsere inhaltliche Hypothese, dass es 5 Faktoren gibt, mit Hilfe der oblique Rotierten ML-EFA.
five_factor_ML <- fa(dataFR, nfactors = 5, rotate = "oblimin", fm = "ml")
five_factor_ML$STATISTIC
## [1] 44.07717
five_factor_ML$PVAL # Modell wird durch die Daten nicht verworfen
## [1] 0.3031972
Die Daten scheinen unserem Modell mit 5 Faktoren nicht zu widersprechen. Schauen wir uns die Faktorladungen an, um die Faktoren inhaltlich zu interpretieren.
five_factor_ML$loadings # auch nochmal ohne [,] um die Ausblendehilfe von psych als Unterstützung für die Zuordnung zu nutzen
##
## Loadings:
## ML4 ML3 ML1 ML2 ML5
## E1 0.688
## E2 -0.643 0.155
## E3 0.813 0.117
## N1 0.433 -0.340
## N2 0.986
## N3 0.920
## A1 -0.227 0.160 0.203 -0.276
## A2 0.655 0.242
## A3 -0.186 0.159 0.197
## C1 0.170 0.208 0.179 -0.355
## C2 0.730
## C3 0.177 0.140 0.145 -0.278
## O1 0.998
## O2 0.126 0.247 -0.214 -0.106
## O3 0.229 0.163 0.189
##
## ML4 ML3 ML1 ML2 ML5
## SS loadings 2.091 1.328 1.274 1.194 0.869
## Proportion Var 0.139 0.089 0.085 0.080 0.058
## Cumulative Var 0.139 0.228 0.313 0.392 0.450
five_factor_ML$loadings[,] # alle Dezimalstellen anzeigen
## ML4 ML3 ML1 ML2 ML5
## E1 0.688468319 -0.03665453 -0.012193249 0.02787857 0.071260642
## E2 -0.643272367 0.08577303 0.155240994 0.09656023 0.065735002
## E3 0.812764373 -0.05121616 0.116888237 0.01806219 -0.076129365
## N1 -0.022613674 0.43257679 -0.339987402 -0.07415779 0.069005051
## N2 0.017710987 -0.02860957 0.985613730 -0.02537580 0.026025668
## N3 -0.004230572 0.91962125 -0.026913718 -0.00657639 -0.013229825
## A1 -0.227423669 0.15957275 0.203222904 -0.08704521 -0.275728683
## A2 0.655281001 0.24198774 0.006033115 0.01175666 0.092166503
## A3 -0.186178370 0.15918711 0.095110929 0.19652141 0.097407381
## C1 0.056714149 0.16981702 0.207576065 0.17890272 -0.354809510
## C2 -0.015398173 0.01818279 0.084998138 0.01776008 0.729941525
## C3 0.014355504 0.17741965 0.140061799 0.14464044 -0.278431363
## O1 0.002352582 -0.01037494 -0.023817340 0.99812410 -0.001488256
## O2 0.125971083 0.24681060 0.029753625 -0.21358237 -0.106264287
## O3 -0.080478366 0.22862665 0.163239313 0.18923154 0.082298847
Durch die Rotation sind auch hier die Faktoren anders nummeriert. Der erste Faktor ist hier ML4 (dieser Faktor ist der erste in der Liste, da hier der Eigenwerte nach Rotation maximal ist; vor Rotation hatte ML4 den viert größten Eigenwert). Die höchsten Faktorladungen mit diesem Faktor haben die Items
Dies könnte durchaus daran liegen, dass
fa(dataFR, nfactors = 5, rotate = "varimax", fm = "ml")$loadings[,]
## ML4 ML3 ML1 ML2 ML5
## E1 0.679741350 -0.0298612003 0.04717727 -0.06209588 -0.025294029
## E2 -0.616269745 -0.0009785951 0.07502041 0.24020664 -0.043100326
## E3 0.812929628 -0.0734949752 0.16420656 -0.12171911 0.134579147
## N1 -0.060232149 0.5043444919 -0.37548744 0.06060518 -0.063964536
## N2 0.101251528 -0.2329637937 0.95295799 0.10457873 0.106824247
## N3 -0.016430443 0.8479850127 -0.18607020 0.29295427 0.151708177
## A1 -0.213443248 0.1199149173 0.13194526 -0.04987855 0.278471636
## A2 0.644794316 0.2278750098 0.01990198 0.03032774 -0.004249186
## A3 -0.165880085 0.0456542742 0.02042712 0.28908281 -0.018069436
## C1 0.081940401 0.0183155255 0.08750067 0.11679630 0.435079003
## C2 -0.006651558 0.0264955123 0.16832618 0.26419496 -0.667066743
## C3 0.032388740 0.0576353789 0.03621164 0.11108025 0.344258937
## O1 0.056711714 -0.4160253195 -0.23367971 0.84296969 0.231306357
## O2 0.112086056 0.3032964635 0.02602280 -0.14672358 0.100453282
## O3 -0.056937114 0.0975755222 0.08029510 0.29567696 0.021600591
was wahrscheinlich daran liegt, dass die Kovariation zwischen den Faktoren nicht sehr groß ist:
round(five_factor_ML$Phi, 2) # runde auf 2 Nachkommastellen
## ML4 ML3 ML1 ML2 ML5
## ML4 1.00 -0.05 0.14 -0.07 -0.09
## ML3 -0.05 1.00 -0.30 -0.02 -0.02
## ML1 0.14 -0.30 1.00 0.04 0.02
## ML2 -0.07 -0.02 0.04 1.00 0.00
## ML5 -0.09 -0.02 0.02 0.00 1.00
fa(dataFR, nfactors = 5, rotate = "varimax", fm = "ml")$Phi
## NULL
NULL
zeigt hierbei an, dass es das $Phi
-Objekt nicht gibt. Tatsächlich ist die Kovarianzmatrix im orthogonalen Fall die Einheitsmatrix der Dimension R
:
diag(5) # Einheitsmatrix der Dimension 5x5.
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 0
## [2,] 0 1 0 0 0
## [3,] 0 0 1 0 0
## [4,] 0 0 0 1 0
## [5,] 0 0 0 0 1
Modellvergleich: ML-EFA
Wir schauen uns nun die Passung unseres Modells im Vergleich zu einem vier- und einem sechsfaktoriellen Modell an.
four_factor_ML <- fa(dataFR, nfactors = 4, rotate = "oblimin", fm = "ml")
four_factor_ML$STATISTIC
## [1] 73.04596
four_factor_ML$PVAL
## [1] 0.023097
Das vierfaktorielle Modell wird durch die Daten verworfen (
anova(four_factor_ML, five_factor_ML)
## Model 1 = fa(r = dataFR, nfactors = 4, rotate = "oblimin", fm = "ml")
## Model 2 = fa(r = dataFR, nfactors = 5, rotate = "oblimin", fm = "ml")
df | d.df | chiSq | d.chiSq | PR | test | empirical | d.empirical | test.echi | BIC | d.BIC | RMSEA | rms | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
1 | 51 | NA | 73.05 | NA | NA | NA | 105.81 | NA | NA | -174.80 | NA | 0.06 | 0.06 |
2 | 40 | 11 | 44.08 | 28.97 | 0 | 2.63 | 47.94 | 57.87 | 5.26 | -150.32 | -28.97 | 0.03 | 0.04 |
Wir entscheiden uns hier nun für das fünffaktorielle Modell. Nun wollen wir uns das fünffaktorielle Modell noch im Vergleich zum sechsfaktoriellen Modell ansehen.
six_factor_ML <- fa(dataFR, nfactors = 6, rotate = "oblimin", fm = "ml")
six_factor_ML$STATISTIC
## [1] 29.19415
six_factor_ML$PVAL # Modell wird durch die Daten nicht verworfen
## [1] 0.5074152
Dem sechsfaktoriellen Modell widersprechen die Daten genauso wenig, wie dem fünffakoriellen (beide
anova(five_factor_ML, six_factor_ML)
## Model 1 = fa(r = dataFR, nfactors = 5, rotate = "oblimin", fm = "ml")
## Model 2 = fa(r = dataFR, nfactors = 6, rotate = "oblimin", fm = "ml")
df | d.df | chiSq | d.chiSq | PR | test | empirical | d.empirical | test.echi | BIC | d.BIC | RMSEA | rms | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
1 | 40 | NA | 44.08 | NA | NA | NA | 47.94 | NA | NA | -150.32 | NA | 0.03 | 0.04 |
2 | 30 | 10 | 29.19 | 14.88 | 0.14 | 1.49 | 28.63 | 19.31 | 1.93 | -116.60 | -14.88 | 0.00 | 0.03 |
Der
Wir hätten auch mehrere Tests gleichzeitig durchführen können. Allerdings sollten nicht beliebig konkurriende Theorien getestet werden — Stichwort Alpha-Inflation!
anova(four_factor_ML, five_factor_ML, six_factor_ML)
## Model 1 = fa(r = dataFR, nfactors = 4, rotate = "oblimin", fm = "ml")
## Model 2 = fa(r = dataFR, nfactors = 5, rotate = "oblimin", fm = "ml")
## Model 3 = fa(r = dataFR, nfactors = 6, rotate = "oblimin", fm = "ml")
df | d.df | chiSq | d.chiSq | PR | test | empirical | d.empirical | test.echi | BIC | d.BIC | RMSEA | rms | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
1 | 51 | NA | 73.05 | NA | NA | NA | 105.81 | NA | NA | -174.80 | NA | 0.06 | 0.06 |
2 | 40 | 11 | 44.08 | 28.97 | 0.00 | 2.63 | 47.94 | 57.87 | 5.26 | -150.32 | -28.97 | 0.03 | 0.04 |
3 | 30 | 10 | 29.19 | 14.88 | 0.14 | 1.49 | 28.63 | 19.31 | 1.93 | -116.60 | -14.88 | 0.00 | 0.03 |
Appendix C
Faktorladungsmatrix vs. Strukturmatrix
Um die Beziehung zwischen der Faktorladungsmatrix und der Strukturmatrix genauer zu verstehen, schauen wir uns das zweifaktorielle Modell für den (standardisierten) Datensatz dataFR2
genauer an (standardisiert ist hier wichtig, da dies bedeutet, dass die Mittelwerte alle
Dies ist die Messmodellsgleichung, welche die Beziehung zwischen den latenten Variablen und den Messfehlern mit den beobachtbaren Variablen beschreibt. Die Faktorladungsmatrix
Im Abschnitt zur Hauptachsenanalyse hatten wir erkannt, dass der erste Faktor wahrscheinlich der Extraversion und der zweite wahrscheinlich dem Neurotizismus entspricht. Demnach könnten wir
two_factor
) und das oblique-rotiert geschätzte Modell (bereits geschätzt in two_factor_oblimin
) an.
Im orthogonalen Fall ist dies etwas unspannend:
two_factor$loadings[1, 1] # volle Formel für ersten Eintrag in Strukutrmatrix, da Kovarianz der Faktoren = 0
## [1] 0.6918932
two_factor$Structure[1, 1] # erster Eintrag in der Strukturmatrix
## [1] 0.6918932
Offensichtlich sind beide Einträge gleich, was daran liegt, dass die Faktoren als unkorreliert angenommen werden. Nun zum oblique rotierten Fall:
two_factor_oblimin$loadings[1, 1] # erste Faktorladung im obliquen Modell (unterscheidet sich von dem ersten Eintrag der Strukturmatrix)
## [1] 0.7005517
two_factor_oblimin$loadings[1, 1] + two_factor_oblimin$loadings[1, 2]*two_factor_oblimin$Phi[2, 1] # volle Formel für ersten Eintrag in Strukutrmatrix
## [1] 0.6932416
two_factor_oblimin$Structure[1, 1] # erster Eintrag in der Strukturmatrix
## [1] 0.6932416
Hier ist zu sehen, dass sich Faktorladungsmatrix und Strukturmatrix unterscheiden. Der Unterschied ist nicht sehr hoch, da die Korrelation zwischen den beiden Faktoren lediglich bei
R
leicht empirisch überprüfen. Einen Überblick über die Befehle für Matrix-Algebra in R
finden Sie auf der Quick-R Website, auf welche bereits in der Sitzung zu Einführung in lavaan
aufmerksam gemacht wurde. Außerdem ist im Appendix B der Einführungssitzung zu PsyMSc1 bereits eine Einführung in Matrixalgebra gegeben worden. Wir berechnen nun das Matrixprodukt für den oblique rotieren Fall:
two_factor_oblimin$loadings[,] %*% two_factor_oblimin$Phi[,] # Matrixprodukt
## MR1 MR2
## E1 0.6932416 -0.09029319
## E2 -0.6400647 0.03425442
## E3 0.8364933 -0.22391151
## N1 -0.1298599 0.83712330
## N2 0.1785989 -0.58432636
## N3 -0.1283499 0.58695759
two_factor_oblimin$Structure[,] # Strukturmatrix
## MR1 MR2
## E1 0.6932416 -0.09029319
## E2 -0.6400647 0.03425442
## E3 0.8364933 -0.22391151
## N1 -0.1298599 0.83712330
## N2 0.1785989 -0.58432636
## N3 -0.1283499 0.58695759
%*%
signalisiert R
, dass ein Matrixprodukt und keine komponentenweise Mulitplikation durchzuführen ist. [,]
sorgt dafür, dass nur die Matrizen verwendet werden und nicht der zusätzliche Output, der ggf. durch das psych
-Paket mit ausgegeben wird.
Das ganze funktioniert selbstverständlich auch für den fünffaktoriellen oblique rotierten ML-EFA Fall, den wir uns später angesehen haben, als es darum ging, den gesamten (gekürzten) Datensatz mit Hilfe der ML-EFA zu untersuchen. Das zugehörige Objekt, welches das geschätzte Modell enthält, heißt five_factor_ML
:
five_factor_ML$loadings[,] %*% five_factor_ML$Phi[,] # Matrixprodukt
## ML4 ML3 ML1 ML2 ML5
## E1 0.68012693 -0.06717103 0.09828653 -0.01822411 0.009260453
## E2 -0.63791259 0.06552915 0.04403530 0.14413074 0.125484598
## E3 0.83728352 -0.12337423 0.24576500 -0.03084187 -0.146382973
## N1 -0.09189974 0.53663293 -0.47561225 -0.09483930 0.055041769
## N2 0.15694831 -0.32787315 0.99631061 0.01289071 0.045752916
## N3 -0.04934588 0.92837355 -0.30647378 -0.02583979 -0.032437499
## A1 -0.17544108 0.11612973 0.11370804 -0.06663641 -0.254026633
## A2 0.63566748 0.20739547 0.02728767 -0.03702446 0.027811781
## A3 -0.20229561 0.13312175 0.03054772 0.20947982 0.112785431
## C1 0.09812869 0.10806476 0.16374057 0.18026600 -0.359308148
## C2 -0.07172742 -0.02229676 0.09337011 0.02096820 0.732728884
## C3 0.04129630 0.13719615 0.08822234 0.14594730 -0.280624296
## O1 -0.06746981 -0.02332237 0.01902778 0.99723607 -0.003104410
## O2 0.14261965 0.23841159 -0.03792795 -0.22572173 -0.121938771
## O3 -0.08840219 0.17745445 0.09190619 0.19639438 0.088088488
five_factor_ML$Structure[,] # Strukturmatrix
## ML4 ML3 ML1 ML2 ML5
## E1 0.68012693 -0.06717103 0.09828653 -0.01822411 0.009260453
## E2 -0.63791259 0.06552915 0.04403530 0.14413074 0.125484598
## E3 0.83728352 -0.12337423 0.24576500 -0.03084187 -0.146382973
## N1 -0.09189974 0.53663293 -0.47561225 -0.09483930 0.055041769
## N2 0.15694831 -0.32787315 0.99631061 0.01289071 0.045752916
## N3 -0.04934588 0.92837355 -0.30647378 -0.02583979 -0.032437499
## A1 -0.17544108 0.11612973 0.11370804 -0.06663641 -0.254026633
## A2 0.63566748 0.20739547 0.02728767 -0.03702446 0.027811781
## A3 -0.20229561 0.13312175 0.03054772 0.20947982 0.112785431
## C1 0.09812869 0.10806476 0.16374057 0.18026600 -0.359308148
## C2 -0.07172742 -0.02229676 0.09337011 0.02096820 0.732728884
## C3 0.04129630 0.13719615 0.08822234 0.14594730 -0.280624296
## O1 -0.06746981 -0.02332237 0.01902778 0.99723607 -0.003104410
## O2 0.14261965 0.23841159 -0.03792795 -0.22572173 -0.121938771
## O3 -0.08840219 0.17745445 0.09190619 0.19639438 0.088088488
Hier alle Einträge auf Gleichheit zu untersuchen, ist sehr mühsam. Wir können dies viel einfacher mit einer Differenz tun:
five_factor_ML$loadings[,] %*% five_factor_ML$Phi[,] - five_factor_ML$Structure[,]
## ML4 ML3 ML1 ML2 ML5
## E1 0 0 0 0 0
## E2 0 0 0 0 0
## E3 0 0 0 0 0
## N1 0 0 0 0 0
## N2 0 0 0 0 0
## N3 0 0 0 0 0
## A1 0 0 0 0 0
## A2 0 0 0 0 0
## A3 0 0 0 0 0
## C1 0 0 0 0 0
## C2 0 0 0 0 0
## C3 0 0 0 0 0
## O1 0 0 0 0 0
## O2 0 0 0 0 0
## O3 0 0 0 0 0
Da hier nur Nullen herauskommen, scheinen die Ausdrücke identisch zu sein!
Appendix D
Berechnen von Eigenwerten und Kommunalitäten mit Hilfe von und
Im vorigen Abschnitt hatten wir bemerkt, wie leicht die Strukturmatrix aus der Faktorladungsmatrix
two_factor_ML$communality
## E1 E2 E3 N1 N2 N3
## 0.4715292 0.4068085 0.7128628 0.6964745 0.3520527 0.3565427
diag(two_factor_ML$loadings[,] %*% two_factor_ML$Phi[,] %*% t(two_factor_ML$loadings[,]))
## E1 E2 E3 N1 N2 N3
## 0.4715292 0.4068085 0.7128628 0.6964745 0.3520527 0.3565427
diag(two_factor_ML$Structure[,] %*% t(two_factor_ML$loadings[,]))
## E1 E2 E3 N1 N2 N3
## 0.4715292 0.4068085 0.7128628 0.6964745 0.3520527 0.3565427
In allen 3 Fällen kommen die Kommunalitäten heraus. diag
fordert hierbei die Diagonalelemente einer quadratischen Matrix an. Um nun die Eigenwerte zu erhalten, müssen wir das Matrixprodukt nur umdrehen: Faktorladungsmatrix’ Strukturmatrix, bzw. das Transponierungszeichen tauschen (Strukturmatrix’ Faktorladungsmatrix ginge auch!):
two_factor_ML$Vaccounted # Eigenwerte nach Rotation und Extraktion in SS loadings
## ML1 ML2
## SS loadings 1.6053881 1.3908823
## Proportion Var 0.2675647 0.2318137
## Cumulative Var 0.2675647 0.4993784
## Proportion Explained 0.5357955 0.4642045
## Cumulative Proportion 0.5357955 1.0000000
diag(two_factor_ML$Phi[,] %*% t(two_factor_ML$loadings[,]) %*% two_factor_ML$loadings[,])
## ML1 ML2
## 1.605388 1.390882
diag(t(two_factor_ML$Structure[,]) %*% two_factor_ML$loadings[,])
## ML1 ML2
## 1.605388 1.390882
In der ersten Zeile von two_factor_ML$Vaccounted
stehen die Eigenwerte nach Rotation und Extraktion. Diese sind identisch zur Diagonale der beiden Matrixprodukte!
Appendix E
Prüfen der Voraussetzungen mit der Mahalanobisdistanz und Mardia's Test
Auf multivariate Normalverteilung können wir beispeilsweise **deskriptiv** prüfen, indem wir die Mahalanobisdistanz (die Distanz vom gemeinsame Zentroiden; dem Mittelwert über alle Variablen; unter Berücksichtigung der Kovariation im Datensatz) plotten und sie mit einerMahalanobis_Distanz <- mahalanobis(x = dataFR, cov = cov(dataFR), center = colMeans(dataFR)) # Berechnen der Mahalanobisdistanz
hist(Mahalanobis_Distanz, col = "skyblue", border = "blue", freq = F, breaks = 15) # Histogramm
lines(x = seq(0, max(Mahalanobis_Distanz), 0.01), y = dchisq(x = seq(0, max(Mahalanobis_Distanz), 0.01), df = 15), col = "darkblue", lwd = 4) # Einzeichnen der Dichte

Sie können ja mal Einstellungen verändern und sich deren Konsequenz für die Grafik ansehen!
Das Histogramm scheint nicht perfekt zur
Die Funktion mahalanobis
berechnet die Mahalanobisdistanz pro Proband:in. Als Datenargument braucht sie eine Matrix x
. Die Mahalanobisdistanz ist ein Distanzmaß, welches die korrelative Struktur in den Daten berücksichtigt. Wir übergeben daher mit cov = cov(dataFR)
der Funktion mahalanobis
die empirische Kovarianzmatrix unserer Daten (cov(dataFR)
), um diese Struktur mit zuberücksichtigen. Außerdem müssen die Variablen und deren Variation relativ zu einem Zentroiden angegeben werden. Der Zentroid wird dem center
-Argument übergeben. Wir brauchen also für jede Variable deren Mittelwert. Dies machen wir mit colMeans
.
colMeans(dataFR)
## E1 E2 E3 N1 N2 N3 A1 A2 A3 C1
## 2.558140 2.968992 3.217054 3.372093 3.131783 3.852713 2.620155 3.596899 2.286822 3.100775
## C2 C3 O1 O2 O3
## 3.131783 4.000000 3.945736 2.077519 4.240310
Der hist
-Befehl erzeugt schließlich ein Histogramm der Mahalanobisdistanzen. Mit den Argumenten col = "skyblue"
und border = "blue"
setzten wir die Farben des Histogramms fest. Mit freq = F
sagen wir, dass wir nicht die absoluten sondern die relativen Häufigkeiten angezeigt haben wollen (dies brauchen wir um anschließend die Dichte der breaks = 15
beschließen wir, dass insgesamt ca. 15 Balken gezeichnet werden sollen.
Schließlich zeichnen wir mit lines
eine Line, welche als x-Argument x = seq(0, max(Mahalanobis_Distanz), 0.01)
eine Sequenz von Zahlen von 0 bis zur maximalen Mahalanobisdistanz erhält und in 0.01 Schritten wächst. Gegen diese x-Werte zeichnen wir die Dichte der y = dchisq(x = seq(0, max(Mahalanobis_Distanz), 0.01), df = 15)
. col = "darkblue"
und lwd = 4
setzten jeweils die Linienfarbe und Liniendicke fest. Weitere Informationen zu Verteilungen und wie man diese in R
umsetzt, können im R-Wiki zu Verteilungen, in Wikipedia zu Verteilungen und Dichten oder in einer Kurzzusammenfassung auf statmethods nachgelesen werden. Grundlagen hierzu können außerdem in Eid et al. (2017) in Kapitel 7 ab Seite 171 gefunden werden.
Außerdem können wir auch noch Mardia’s Test auf multivariate Normalverteilung verwenden. Diesen gibt es bspw. im R
-Paket MVN
(für Multi-Variate Normal-distribution). In diesem Paket verwenden wir die mvn
-Funktion, um Mardia’s Test zu verwenden:
library(MVN)
mvn(data = dataFR, mvnTest = "mardia")
## $multivariateNormality
## Test Statistic p value Result
## 1 Mardia Skewness 814.236119391994 0.000288124299329733 NO
## 2 Mardia Kurtosis 1.53049602635342 0.125893996300673 YES
## 3 MVN <NA> <NA> NO
##
## $univariateNormality
## Test Variable Statistic p value Normality
## 1 Anderson-Darling E1 4.9942 <0.001 NO
## 2 Anderson-Darling E2 4.0628 <0.001 NO
## 3 Anderson-Darling E3 4.0445 <0.001 NO
## 4 Anderson-Darling N1 5.1846 <0.001 NO
## 5 Anderson-Darling N2 4.5090 <0.001 NO
## 6 Anderson-Darling N3 8.7268 <0.001 NO
## 7 Anderson-Darling A1 5.2696 <0.001 NO
## 8 Anderson-Darling A2 5.8351 <0.001 NO
## 9 Anderson-Darling A3 9.9801 <0.001 NO
## 10 Anderson-Darling C1 4.7848 <0.001 NO
## 11 Anderson-Darling C2 4.0206 <0.001 NO
## 12 Anderson-Darling C3 8.3567 <0.001 NO
## 13 Anderson-Darling O1 7.4239 <0.001 NO
## 14 Anderson-Darling O2 7.7934 <0.001 NO
## 15 Anderson-Darling O3 11.2213 <0.001 NO
##
## $Descriptives
## n Mean Std.Dev Median Min Max 25th 75th Skew Kurtosis
## E1 129 2.558140 1.3044813 3 1 5 1 3 0.32427985 -1.03222800
## E2 129 2.968992 1.2433447 3 1 5 2 4 0.13071545 -0.98604106
## E3 129 3.217054 1.2683727 3 1 5 2 4 -0.20289793 -0.99959701
## N1 129 3.372093 1.3232876 4 1 5 2 4 -0.31730421 -1.14311645
## N2 129 3.131783 1.1550152 3 1 5 2 4 -0.01411112 -0.95870960
## N3 129 3.852713 1.0975881 4 1 5 3 5 -0.87031995 -0.10813842
## A1 129 2.620155 1.3761007 2 1 5 1 4 0.33995849 -1.17089544
## A2 129 3.596899 1.1693974 4 1 5 3 5 -0.37750788 -0.99209594
## A3 129 2.286822 1.3761007 2 1 5 1 4 0.55150074 -1.20186405
## C1 129 3.100775 1.1239228 3 1 5 2 4 -0.16388511 -0.85645458
## C2 129 3.131783 1.2770807 3 1 5 2 4 -0.17800768 -1.01121119
## C3 129 4.000000 1.0231691 4 1 5 3 5 -0.86845894 -0.01507734
## O1 129 3.945736 1.0178875 4 1 5 3 5 -0.90696625 0.43639201
## O2 129 2.077519 1.1012238 2 1 5 1 3 0.75369527 -0.39530983
## O3 129 4.240310 0.9080817 4 1 5 4 5 -1.16861246 0.88182592
Hier bekommen wir einiges an Output. Der erste Block enthält den Test auf multivariate Normalität (unter $multivariateNormality
). Der nächste Unterpunkt enthält Informationen zur univariate Normalität (unter $univariateNormality
) sowie einen Block zu Deskriptivstatistiken (unter $Descriptives
). Wir konzentrieren uns nur auf den multivariaten Test:
## ## $multivariateNormality
## ## Test Statistic p value Result
## ## 1 Mardia Skewness 814.236119391994 0.000288124299329737 NO
## ## 2 Mardia Kurtosis 1.53049602635342 0.125893996300673 YES
## ## 3 MVN <NA> <NA> NO
## ##
In der ersten Spalte steht der Test Mardia Skewness
oder Mardia Kurtosis
, wo mit multivariater Schiefe (Skewness) und Kurtosis (Wölbung) untersucht wird, ob diese von der Normalverteilung abweichen. Statistic
enthält Informationen zu Mardia’s Teststatistik, welche auch einen zugehörigen p value
hat. Unter Result
steht eine Entscheidung, ob die Schiefe oder Wölbung jeweils von einer Normalverteilung stammen könnte. Ganz unten steht MVN
, also eine globale Entscheidung. Hier steht leider NO
, also kann nicht von multivariater Normalverteilung gesprochen werden, denn die Null-Hypothese, dass die Daten die selbe multivariate Schiefe und Wölbung wie die einer Normalverteilung aufzeigen, musste auf dem
Literatur
Brandt H. (2020). Exploratorische Faktorenanalyse (EFA). In Moosbrugger H., Kelava A. (eds) Testtheorie und Fragebogenkonstruktion. Berlin, Heidelberg: Springer. https://doi.org/10.1007/978-3-662-61532-4_23
Eid, M., Gollwitzer, M., & Schmitt, M. (2017). Statistik und Forschungsmethoden (5. Auflage, 1. Auflage: 2010). Weinheim: Beltz.
- Blau hinterlegte Autor:innenangaben führen Sie direkt zur universitätsinternen Ressource.