In martscht/PsyMSc1: Übungen für die Seminare im Sommersemester 2020

library(learnr)       # package to generate interactive HTML
library(gradethis)    # grade Code from input in HTML
library(shiny)        # nice appearance in R
library(fontawesome)  # nice fonts
library(corrplot)     # Korrelationsmatrix grafisch darstellen
library(psych)        # EFA durchführen
library(GPArotation)  # EFA Lösung rotieren

data(Big5, package = 'PsyMSc1')
knitr::opts_chunk$set(exercise.checker = gradethis::grade_learnr)

Einleitung zur exploratorischen Faktorenanalyse (EFA)

Forscher der Psychologie oder anderer Natur-, Sozial- und Geisteswissenschaften interessieren sich häufig dafür, dass sich Daten runterbrechen lassen auf einige wenige entscheidende Faktoren, 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 Varibalen untersucht werden. Ein theoriegenerierendes Verfahren, das hierzu häufig verwendet wird, ist die exploratorische Faktorenanalyse (EFA). Wir wollen dieses Verfahren zur Auswertung von Beziehungen zwischen Variablen in R näher kennenlernen.

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

Den gesamten R-Code, der in dieser Sitzung genutzt wird, können Sie r fontawesome::fa("download") hier herunterladen.

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 r fontawesome::fa("download") openpsychometrics.org als .zip downloadbar. Bisher haben über 19700 Probanden aus der ganzen Welt 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. Im Datensatz Big5.rda befinden sich 15 Items aus dem Big-5 Persönlichkeits Fragebogen. Hier werden von diesen 10 Items jeweils die ersten drei verwendet. Der Itemwortlaut der verwendeten Items ist

library(htmlTable)
output <- c(c("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."))
output <- matrix(output, ncol = 2, nrow = 3*5, byrow = T)
htmlTable(output, header =  paste(c("Item Nr.", "Item")),
         caption="Itemwortlaut", align = "cl")

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.

# Datensatz laden 
data("Big5", package = "PsyMSc1") # shortend Big 5 Questionnaire Data Set
head(Big5, n = 10) # gebe die ersten 10 Zeilen aus

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. Versuchen Sie doch mal den Datensatz zu kürzen und anschließend herauszufinden, wie viele Zeilen und Spalten der Datensatz der in Frankreich lebenden Menschen hat.

data_France <- ...
...

# Sie benötigen hierzu den Datensatz Big5 und müssen bei diesem die Zeilen ansprechen, um hier nur diejenigen auszuwählen, in denen Menschen stehen, die in Frankreich leben, bspw. gibt 
Big5[1, ] # die erste Zeile aus

Big5$country == "FR" # ist ein Vektor, der an den Stellen TRUE stehen hat, in welchen die Person in Frankreich lebt
Big5[Big5$country == "FR", ] # zeigt dijenigen Personen, die in Frankreich leben, jetz müssen Sie dieses Objekt nur noch abspeichern und mit dim() die Zeilen und Spalten abfragen!

# dies ist die korrekte Lösung:
data_France <- Big5[Big5$country == "FR", ]
dim(data_France)

# For further use
data_France <- Big5[Big5$country == "FR", ]

# For further use in exercises
data_France <- Big5[Big5$country == "FR", ]

Dem Output sollte zu entnehmen sein, dass data_France r dim(data_France)[1] Zeilen (also Probanden, die in Frankreich leben) und r dim(data_France)[2] Spalten (also Variablen) enthält. Für die weiteren Analysen brauchen wir die demografischen Variablen in dem Datensatz der in Frankreich lebenden Teilnehmer nicht mehr. Aus diesem Grund speichern wir den Datensatz noch einmal ohne die ersten 4 Spalten ab.

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

# For further use
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"

# For further use in exercises
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3

Auf den ersten Blick scheinen die Items der gleichen Skala 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)

Visualisieren Sie doch mal die Korrelationsmatrix in unserem gekürzten Datensatz dataFR2.

# Visualisierte Korrelationsmatrix
corrplot(corr = cor(...), # 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

# 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

Hier ist nun deutlich zu sehen, 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 $\chi^2$-Test (Likelihood-Quotiententest/ Likelihood-Differenzentest/ $\chi^2$-Differenzentest) gegen konkurrierende Modelle vergleichen. Hierbei wollen wir die oblique rotierte und die orthogonal rotierte Lösung vergleichen und hinsichtlich unserer Daten interpretieren. Das Modell, was wir an unsere Daten anpassen wollen, sieht für 6 Variablen ($V_1,\dots,V_6$) im Allgemeinen erst einmal so aus:

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 {#Parallelanalyse_1}

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 (PCA, principal component analysis) kennengelernt. Hier werden vielfach (z.B. 1000 Mal) 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 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 (Ellbow-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 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).

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)

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.

fa.parallel(dataFR2, fa = "fa")

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. Da beim Resampling die Systematiken unseres Datensatzes verloren gehen, sollte das Resampling einen ähnlichen Eigenwerteverlauf wie die unabhängig gezogenen Datensätze erzeugen.

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 {#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. Wir wollen hier orthogonal varianzmaximierend (also varimax) rotieren.

fa(dataFR2, ...)

fa(dataFR2, nfactors = ..., rotate = ...)

# Finale Lösung:
fa(dataFR2, nfactors = 2, rotate = "varimax")

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. 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 ($h^2$), also den Anteil an systematischer Variation, die auf die 2 Faktoren zurückzuführen ist (diese kann ähnlich der Reliabilität interpretiert werden). u2 ist die "uniqueness" ($u^2$), also der unerklärte Anteil. Offensichtlich gilt $u^2 = 1-h^2$ oder $h^2 + u^2 = 1$. Unter den Faktorladungen erhalten wir Informationen über die Faktoren. SS loadings steht für "Sum of Squares loadings", also die Quadratsumme der Faktorladungen. Diese ist gleich dem Eigenwert: $\theta_j = \Sigma_{i=1}^p\lambda_{ij}^2 = \lambda_{1j}^2+\dots+\lambda_{pj}^2$ (Spaltenquadratsumme der Faktorladungen), mit $p=$ Anzahl an Variablen (hier $p=6$). Proportion Var betitelt den Anteil der Variation, der durch die jeweiligen Faktoren erklärt werden kann. Cumulative Var kumuliert, also summiert, diese Anteile auf bis zum jeweiligen Faktor auf ($\text{CumVar}i = \sum{j=1}^i\theta_j = \theta_1+\dots+\theta_i$, also $\text{CumVar}_1=\theta_1$ und $\text{CumVar}_2=\theta_1+\theta_2$). Proportion 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 $h^2$ nicht 1 und wir können nicht unbedingt davon ausgehen, dass die Faktoren die gesamte Variation der Daten erklären.

Da durch diesen Befehl sehr viele Informationen ausgegeben werden, 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. Speichern Sie das Objekt unter dem Namen two_factor ab und wenden die Funktion names darauf an.

two_factor <- ...
...

# Intro into excercise
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3

# Begin grading:
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
res <- names(two_factor)

grade_result(
  fail_if(~ (!identical(.result, res)), 'Haben Sie den korrekten Datensatz angesprochen und die richtigen Argumente verwendet?'),
  pass_if(~ identical(.result, res)),
  correct = 'Sehr gut! Sie haben das Objekt aus der Übung zuvor ["fa(dataFR2, nfactors = 2, rotate = "varimax")"] abgespeichert und können sich nun ansehen, welche Informationen es enthält!',
  incorrect = 'Leider falsch.',
  glue_correct = '{.correct}',
  glue_incorrect = '{.incorrect} {.message}')

# for later use
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")

# For further use in exercises
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")

Beispielsweise erhalten wir mit $loadings die Faktorladungsmatrix sowie Informationen über die Eigenwerte.

two_factor$loadings

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. Die Faktorladungsmatrix wird auch manchmal Mustermatrix genannt. Die Strukturmatrix enthält Informationen über die Korrelation der Items mit den jeweiligen Faktoren. Versuchen Sie diese über ...$... auszugeben.

names(two_factor) # zeigt, welche Argumente möglich sind!

# Intro into excercise
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")

# Begin grading:
res <- two_factor$Structure 

grade_result(
  fail_if(~ (!identical(.result, res) & !identical(.result, names(two_factor))), 'Schauen Sie bspw. mit names() nochmal nach, welche Argumente möglich wären.'),
    fail_if(~ (identical(.result, names(two_factor))), 'names sollte Ihnen helfen das Argument zu finden. Wahrscheinlich haben Sie aus Versehen *Submit Answer* anstatt *Run Code* gedrückt.'),
  pass_if(~ identical(.result, res)),
  correct = 'Sehr gut! $Structure ist das richtige Argument. Der Name verrät hier, dass es sich um die Strukturmatrix handel!',
  incorrect = 'Leider falsch.',
  glue_correct = '{.correct}',
  glue_incorrect = '{.incorrect} {.message}')

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.

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 <- ...
two_factor_oblimin

# Intro into excercise
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")

# Begin grading:
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")
res <- two_factor_oblimin

grade_result(
  fail_if(~ (!identical(.result, res)), 'Haben Sie den korrekten Datensatz angesprochen und die richtigen Argumente verwendet?'),
  pass_if(~ identical(.result, res)),
  correct = 'Sehr gut! Sie haben das Objekt aus der Übung zuvor abgespeichert und können sich nun ansehen, welche Informationen es enthält!',
  incorrect = 'Leider falsch.',
  glue_correct = '{.correct}',
  glue_incorrect = '{.incorrect} {.message}')

# for later use
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")

# For further use in exercises
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
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_oblimin$loadings # Ladungsmatrix

Der Ladungsmatrix ist post-hoc interpretiert auch zu entnehmen, 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 $\Phi$.

names(two_factor_oblimin) # zeigt, welche Argumente möglich sind!

# Intro into excercise
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")

# Begin grading:
res <- two_factor_oblimin$Phi # Korrelationsmatrix der Faktoren

grade_result(
  fail_if(~ (!identical(.result, res) & !identical(.result, names(two_factor_oblimin))), 'Schauen Sie bspw. mit names() nochmal nach, welche Argumente möglich wären.'),
    fail_if(~ (identical(.result, names(two_factor_oblimin))), 'names sollte Ihnen helfen das Argument zu finden. Wahrscheinlich haben Sie aus Versehen *Submit Answer* anstatt *Run Code* gedrückt.'),
  pass_if(~ identical(.result, res)),
  correct = 'Sehr gut! $Phi ist das richtige Argument. Phi wird im Rahmen der Kovarianzanalyse häufig für die Benennung der latenten Kovarianzmatrix  verwendet - so auch hier!',
  incorrect = 'Leider falsch.',
  glue_correct = '{.correct}',
  glue_incorrect = '{.incorrect} {.message}')

Als neue Information entnehmen wir der Korrelationsmatrix der Faktoren, dass die beiden Faktoren negativ korreliert sind zu r round(two_factor_oblimin$Phi[2], 2). Nun wollen wir nachschauen, ob sich tatsächlich die Strukurmatrix im oblique-rotierten Fall von der Faktorladungsmatrix unterscheidet:

names(two_factor_oblimin) # zeigt, welche Argumente möglich sind!

# Intro into excercise
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")

# Begin grading:
res <- two_factor_oblimin$Structure 

grade_result(
  fail_if(~ (!identical(.result, res) & !identical(.result, names(two_factor_oblimin))), 'Schauen Sie bspw. mit names() nochmal nach, welche Argumente möglich wären.'),
    fail_if(~ (identical(.result, names(two_factor_oblimin))), 'names sollte Ihnen helfen das Argument zu finden. Wahrscheinlich haben Sie aus Versehen *Submit Answer* anstatt *Run Code* gedrückt.'),
  pass_if(~ identical(.result, res)),
  correct = 'Sehr gut! $Structure ist das richtige Argument. Der Name verrät hier, dass es sich um die Strukturmatrix handel!',
  incorrect = 'Leider falsch.',
  glue_correct = '{.correct}',
  glue_incorrect = '{.incorrect} {.message}')

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 r round(two_factor_oblimin$Phi[2], 2) betragsmäßig nicht sonderlich groß ausfällt. Weitere Informationen und wie die beiden Matrizen ineinander überführbar sind, erfahren Sie im Appendix B.

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 $\chi^2$-verteilt sein, siehe hierzu im Appendix C nach). Mit Hilfe dieser Verteilungsannahme können wir die Maximum-Likelihood-Schätzmethode nutzen, um die Parameter in unserem Modell zu schätzen. Die Likelihood ist die Wahrscheinlichkeit unserer Daten gegeben das Modell. Sie hängt somit von den beobachteten Daten ab (den Ausprägungen der Personen auf den Variablen), hat die Gestalt unseres Modells (hier: Faktorstruktur mit normalverteilten Faktoren und Fehlern) und wird parametrisiert durch die Koeffizienten (Parameter) in unserem Modell (hier: $\lambda_{..}$ und $\theta_{..}$ im unkorrelierten Modell vor Rotation; Achtung: $\theta_{..}$ hat zwei Indizes und ist nicht mit den Eigenwerten zu verwechseln, welche nur einen Index haben; bspw. steht $\theta_{22}$ für die modellimplizierte Fehlervarianz des zweiten Items). Die durch das Modell implizierte Kovarianz oder Korrelationsmatrix wird mit $\Sigma$ betitelt und setzt sich folgendermaßen zusammen (Die Matrizen $\Lambda$ und $\Theta$ enthalten jeweils die Faktorladungen $\lambda_{..}$ und die Fehlervarianzen $\theta_{..}$):

$$\Sigma := \Lambda \Lambda' + \Theta.$$

Da die Normalveteilung bereits durch die Kovarianzmatrix sowie die Mittelwerte eindeutig zu bestimmen ist, reicht es, die Kovarianzmatrix mit Hilfe unserer Modellparameter so anzupassen, dass sie möglichst nah an der Kovarianzmatrix der Daten liegt; die Likelihood ist hier dann maximal. Die Parameter, die die Likelihood maximieren, werden Maxmimum-Likelihood-Schätzer genannt. Sie sind die Schätzung der Parameter, die am wahrscheinlichsten eine solche Datenkonstellation hervorrufen; gegeben das Modell.

Um mit Hilfe von fa eine ML-EFA durchzuführen, muss dem Argument fm die entsprechende Bezeichnung "ml" übergeben werden. $STATISTIC (nicht $chi, was den empirisch und nicht likelihoodbasierten $\chi^2$-Wert ausgibt; sinnvoll wenn Modellvoraussetzungen nicht erfüllt sind) und $PVAL entlocken der Analyse (abgespeichert als Objekt) den $\chi^2$ Wert und den zugehörigen p-Wert bei r fa(dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")$dof Freiheitsgraden.

### ML
two_factor_ML <- fa(dataFR2, nfactors = 2, rotate = "oblimin", ...)
two_factor_ML

### ML
two_factor_ML <- fa(dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")
two_factor_ML

### ML for later use
two_factor_ML <- fa(dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")

# For further use in exercises
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")
two_factor_ML <- fa(dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")

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 zuvorigen Analyse sehr ähnlich. Uns interessiert nun die Modellpassung. Der Likelihood Chi Square ist der richtige $\chi^2$-Wert. Diesen entlocken wir nun dem two_factor_ML Objekt:

two_factor_ML$... # Likelihood basierter Chi^2-Wert
two_factor_ML$... # p-Wert

two_factor_ML$STATISTIC # Likelihood basierter Chi^2-Wert
two_factor_ML$PVAL # p-Wert

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. Vielleicht reicht auch ein Faktor aus, um die Variation in unserem Datensatz zu beschreiben? Wir wollen unser Modell mit zwei Faktoren gegen eines mit einem und eines mit drei Faktoren absichern.

Modellvergleich: ML-EFA

# Passt auch eines mit 1 Faktor?
one_factor_ML <- fa(dataFR2, ...)
one_factor_ML$... # Chi-Quadratwert 
one_factor_ML$... # p-Wert

# Passt auch eines mit 1 Faktor?
one_factor_ML <- fa(dataFR2, nfactors = 1, rotate = "oblimin", fm = "ml")
one_factor_ML$STATISTIC # Chi-Quadratwert 
one_factor_ML$PVAL # p-Wert

# for later use
one_factor_ML <- fa(dataFR2, nfactors = 1, rotate = "oblimin", fm = "ml")

# For further use in exercises
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")
two_factor_ML <- fa(dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")
one_factor_ML <- fa(dataFR2, nfactors = 1, rotate = "oblimin", fm = "ml")

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 testen die beiden geschachtelten (die Modell lassen sich auseinander gewinnen: das zweifaktorielle Modell lässt sich aus dem restrikiveren, einfaktoriellen Modell durch Freisetzen der Faktorladungen auf dem 2. Faktor sowie durch Freisetzung der Kovarianz zwischen den Faktoren gewinnen).

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.

# Intro into excercise
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")
two_factor_ML <- fa(dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")
one_factor_ML <- fa(dataFR2, nfactors = 1, rotate = "oblimin", fm = "ml")

# Begin grading:
res <- anova(one_factor_ML, two_factor_ML) # 2 Faktoren passen besser zu den Daten!
sol2 <- anova(two_factor_ML, one_factor_ML) # 2 Faktoren passen besser zu den Daten!

grade_result(
  fail_if(~ (!identical(.result, res) & !identical(.result, sol2)), 'Der Befehl anova muss in diesem Zusammenhang zwei Argumente entgegennehmen: die zwei Objekte der ML-EFAs mit jeweils 1 und 2 Faktoren, welche wir zuvor erstellt haben! Achten Sie entsprechend darauf, dass Sie auch wirklich die ML-EFA-Objekte verwenden.'),
    fail_if(~ (identical(.result, sol2)), 'Das ist an sich RICHTIG, allerdings sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell zuerst übergeben werden. Hier: one_factor_ML, da sonst die df negativ sein können und auch der Chi^2-Wert negativ sein könnte!'),
  pass_if(~ identical(.result, res)),
  correct = 'Sehr gut! Dies ist der Modellvergleich des ein- und zweifaktoriellen Modells. Hier sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell zuerst übergeben werden; also: one_factor_ML, da sonst die df negativ sein können und auch der Chi^2-Wert negativ sein könnte!',
  incorrect = 'Leider falsch.',
  glue_correct = '{.correct}',
  glue_incorrect = '{.incorrect} {.message}')

Im Output des anova Befehls müssen wir im Bereich des Chisq-Tests und nicht beim Emp Chisq (eine Nährung des $\chi^2$ Wertes, wenn Annahmen verletzt sind) nachsehen. Der $\chi^2$-Differenzwert liegt bei r round(anova(one_factor_ML, two_factor_ML)[[4]][2], 3) mit einen zugehörigen p-Wert von de facto r round(anova(one_factor_ML, two_factor_ML)[[5]][2], 7) (e-13 bedeutet $10^{-13}$; also eine Zahl, bei welcher die Dezimalstelle um 13 Stellen nach links verschoben wird, z.B. 1.234e-13$=110^{-13}=0.0000000000001234$ - eine sehr kleine Zahl!). Delta Df (häufig $\Delta df$) gibt die Anzahl an Freiheitsgraden an (hier: df = r anova(one_factor_ML, two_factor_ML)[[2]][2]). Somit wird die Null-Hypothese, dass beide Modell die Daten gleich gut beschreiben verworfen. Wir entscheiden uns 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, ...)
three_factor_ML$... # Chi-Quadratwert 
three_factor_ML$... # p-Wert

# Passt auch eines mit 3 Faktor?
three_factor_ML <- fa(dataFR2, nfactors = 3, rotate = "oblimin", fm = "ml")
three_factor_ML$STATISTIC # Chi-Quadratwert 
three_factor_ML$PVAL # p-Wert

# for later use
three_factor_ML <- fa(dataFR2, nfactors = 3, rotate = "oblimin", fm = "ml")

# For further use in exercises
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")
two_factor_ML <- fa(dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")
one_factor_ML <- fa(dataFR2, nfactors = 1, rotate = "oblimin", fm = "ml")
three_factor_ML <- fa(dataFR2, nfactors = 3, rotate = "oblimin", fm = "ml")

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:

# Intro into excercise
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")
two_factor_ML <- fa(dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")
one_factor_ML <- fa(dataFR2, nfactors = 1, rotate = "oblimin", fm = "ml")
three_factor_ML <- fa(dataFR2, nfactors = 3, rotate = "oblimin", fm = "ml")

# Begin grading:
res <- anova(two_factor_ML, three_factor_ML) # 2 Faktoren passen besser zu den Daten!
sol2 <- anova(three_factor_ML, two_factor_ML) # 2 Faktoren passen besser zu den Daten!

grade_result(
  fail_if(~ (!identical(.result, res) & !identical(.result, sol2)), 'Der Befehl anova muss in diesem Zusammenhang zwei Argumente entgegennehmen: die zwei Objekte der ML-EFAs mit jeweils 2 und 3 Faktoren, welche wir zuvor erstellt haben! Achten Sie entsprechend darauf, dass Sie auch wirklich die ML-EFA-Objekte verwenden.'),
    fail_if(~ (identical(.result, sol2)), 'Das ist an sich RICHTIG, allerdings sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell zuerst übergeben werden. Hier: two_factor_ML, da sonst die df negativ sein können und auch der Chi^2-Wert negativ sein könnte!'),
  pass_if(~ identical(.result, res)),
  correct = 'Sehr gut! Dies ist der Modellvergleich des zwei- und dreifaktoriellen Modells. Hier sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell zuerst übergeben werden; also: two_factor_ML, da sonst die df negativ sein können und auch der Chi^2-Wert negativ sein könnte!',
  incorrect = 'Leider falsch.',
  glue_correct = '{.correct}',
  glue_incorrect = '{.incorrect} {.message}')

Der $\chi^2$-Differenzwert liegt hier bei r round(anova(two_factor_ML, three_factor_ML)[[4]][2], 3) mit einen zugehörigen p-Wert von r round(anova(two_factor_ML, three_factor_ML)[[5]][2], 3). Delta Df liegt bei r anova(two_factor_ML, three_factor_ML)[[2]][2] ($\Delta df$ = r anova(two_factor_ML, three_factor_ML)[[2]][2]). Somit wird die Null-Hypothese, dass beide Modell die Daten gleich gut beschreiben, bzw. dass das sparsamere Modell die Daten genauso gut beschreiben kann, wie das komplexere Modell ($H_0:\Sigma_{3-Fakt.} = \Sigma_{2-Fakt.}$; im Gegensatz zum 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.

ML-EFA für den gesamten (gekürtzten) Datensatz {#fivefactorML}

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üzt. Wir müssen die Funktion fa.parallel diesmal auf den vollen (gekürtzten) Datensatz anwenden; namlich auf dataFR.

fa.parallel(x = dataFR)

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.

### ML
five_factor_ML <- ...
five_factor_ML$STATISTIC
five_factor_ML$PVAL # Modell wird durch die Daten nicht verworfen

### ML
five_factor_ML <- fa(dataFR, nfactors = 5, rotate = "oblimin", fm = "ml")
five_factor_ML$STATISTIC
five_factor_ML$PVAL # Modell wird durch die Daten nicht verworfen

### ML
five_factor_ML <- fa(dataFR, nfactors = 5, rotate = "oblimin", fm = "ml")

# For further use in exercises
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")
two_factor_ML <- fa(dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")
one_factor_ML <- fa(dataFR2, nfactors = 1, rotate = "oblimin", fm = "ml")
three_factor_ML <- fa(dataFR2, nfactors = 3, rotate = "oblimin", fm = "ml")
five_factor_ML <- fa(dataFR, nfactors = 5, rotate = "oblimin", fm = "ml")

Die Daten scheinen unserem Modell mit 5 Faktoren nicht zu wiedersprechen. Schauen wir uns die Faktorladungen an, um die Faktoren inhaltlich zu interpretieren.

five_factor_ML$loadings

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 $E_1$, $E_2$, $E_3$ und $A_2$. Somit könnte man diesen am ehesten post-hoc (die Theorie wird also aus den Daten generiert; es sind auch andere Interpretationsansätze zulässig) als Extraversion interpretieren. Allerdings scheinen die Items der Extraversion einiges mit jenen der Verträglichkeit ($A_{...}$) gemeinsam zu haben. Dies könnte mit unter damit zusammen hängen, dass diese beiden Items am ehesten etwas mit sozialer Erwünschtheit zu tun haben. Auf dem Faktor ML3 laden vor allem die Items $N_1$ und $N_3$. Allerdings lädt $N_2$ besonders auf ML1.
Dies könnte durchaus daran liegen, dass $N_1$ ("I get stressed out easily.") und $N_3$ ("I worry about things.") negativ kodiert sind, während $N_2$ ("I am relaxed most of the time.") positiv kodiert ist. Somit scheint ML3 ein Faktor der Sorgen, also des Neurotizismus zu sein, während ML1 eher für einen Faktor der Gelassenheit spricht; beispielsweise laden hier auch positiv Items der Extraversion und Verträglichkeit. Das die Items des Neurotizismus auf unterschiedlichen Faktoren laden und unterschiedliche Vorzeichen aufweisen, kann für Methodeneffekte sprechen (Unterschiede die zustande kommen, da unterschiedliche Methoden, hier: Itemvormulierungen [positiv vs. negativ], verwendet werden.). Auch auf ML2 und ML5 laden jeweils nur ein Item besonders stark: $O_1$ auf ML2 und $C_2$ auf ML5. Insgesamt muss geschlussfolgert werden, dass zwar die fünffaktorielle Struktur durch die Daten nicht verworfen wird, aber dass die oblique rotierte Lösung keine eindeutige Zuordnung der Items aufweist. Allerdings bringt auch eine varimax-rotierte Lösung keine Verbesserung der Interpretierbarkeit, da diese neben der Einfachstruktur in der Faktorladungsmatrix noch die Unkorreliertheit der Faktoren berücksichtigen muss (in der varimax-rotierten Lösung sind dafür die Konstrukte nicht überlappend, was allerdings auch eine strenge Annahme ist):

fa(dataFR, nfactors = 5, rotate = "varimax", fm = "ml")$loadings

was wahrscheinlich daran liegt, dass die Kovariation zwischen den Faktoren nicht sehr groß ist:

round(five_factor_ML$Phi, 2) # runde auf 2 Nachkommastellen
fa(dataFR, nfactors = 5, rotate = "varimax", fm = "ml")$Phi

NULL zeigt hierbei an, dass es das $Phi -Objekt nicht gibt. Tatsächlich ist die Kovarianzmatrix im orthogonalen Fall die Einheitsmatrix der Dimension $5\times5$: $$\begin{pmatrix} 1& 0&0&0&0 \ 0& 1&0&0&0 \ 0& 0&1&0&0\ 0& 0&0&1&0 \ 0& 0&0&0&1 \end{pmatrix}$$ In R:

diag(5) # Einheitsmatrix der Dimension 5x5.

Modellvergleich: ML-EFA

Wir schauen uns nun die Passung unseres Modells im Vergleich zu einem vier- und einem sechsfaktoriellen Modell an.

# Passt auch eines mit 4 Faktor?
four_factor_ML <- ...
four_factor_ML$STATISTIC
four_factor_ML$PVAL 

### ML
four_factor_ML <- fa(dataFR, nfactors = 4, rotate = "oblimin", fm = "ml")
four_factor_ML$STATISTIC
four_factor_ML$PVAL 

### ML
four_factor_ML <- fa(dataFR, nfactors = 4, rotate = "oblimin", fm = "ml")

# For further use in exercises
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")
two_factor_ML <- fa(dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")
one_factor_ML <- fa(dataFR2, nfactors = 1, rotate = "oblimin", fm = "ml")
three_factor_ML <- fa(dataFR2, nfactors = 3, rotate = "oblimin", fm = "ml")
five_factor_ML <- fa(dataFR, nfactors = 5, rotate = "oblimin", fm = "ml")
four_factor_ML <- fa(dataFR, nfactors = 4, rotate = "oblimin", fm = "ml")

Das vierfaktorielle Modell wird durch die Daten verworfen ($p<0.05$). Nun zum Modellvergleich:

# führen Sie eine anova zwischen dem vier- und dem fünffaktoriellen Modell durch

# Intro into excercise
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")
two_factor_ML <- fa(dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")
one_factor_ML <- fa(dataFR2, nfactors = 1, rotate = "oblimin", fm = "ml")
three_factor_ML <- fa(dataFR2, nfactors = 3, rotate = "oblimin", fm = "ml")
five_factor_ML <- fa(dataFR, nfactors = 5, rotate = "oblimin", fm = "ml")
four_factor_ML <- fa(dataFR, nfactors = 4, rotate = "oblimin", fm = "ml")
six_factor_ML <- fa(dataFR, nfactors = 6, rotate = "oblimin", fm = "ml")

# Begin grading:
res <- anova(four_factor_ML, five_factor_ML) 
sol2 <- anova(five_factor_ML, four_factor_ML) 

grade_result(
  fail_if(~ (!identical(.result, res) & !identical(.result, sol2)), 'Der Befehl anova muss in diesem Zusammenhang zwei Argumente entgegennehmen: die zwei Objekte der ML-EFAs mit jeweils 4 und 5 Faktoren, welche wir zuvor erstellt haben! Achten Sie entsprechend darauf, dass Sie auch wirklich die ML-EFA-Objekte verwenden.'),
    fail_if(~ (identical(.result, sol2)), 'Das ist an sich RICHTIG, allerdings sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell zuerst übergeben werden. Hier: four_factor_ML, da sonst die df negativ sein können und auch der Chi^2-Wert negativ sein könnte!'),
  pass_if(~ identical(.result, res)),
  correct = 'Sehr gut! Dies ist der Modellvergleich des vier- und fünffaktoriellen Modells. Hier sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell zuerst übergeben werden; also: four_factor_ML, da sonst die df negativ sein können und auch der Chi^2-Wert negativ sein könnte!',
  incorrect = 'Leider falsch.',
  glue_correct = '{.correct}',
  glue_incorrect = '{.incorrect} {.message}')

Die Frage ist nun, für welches Modell wir uns entscheiden?

question("Für welches Modell entscheiden Sie sich und warum?",
  answer("Modell mit 5 Faktoren, da wir annehmen, dass der Big-5 Fragebogen  aus 5 Persönlichkeitsfaktoren besteht."),
  answer("Modell mit 5 Faktoren, weil es das weniger sparsame Modell ist."),
  answer("Modell mit 4 Faktoren, weil es das sparsamere Modell ist. "),
  answer("Modell mit 5 Faktoren, weil der p-Wert (p < 0.05) anzeigt, dass die Null-Hypothese, dass beide Modell gleich gut zu den Daten passen, verworfen wird. Wir entscheiden uns daher für das weniger restriktive Modell mit mehr Parametern; hier: das Modell mit fünf Faktoren.", correct = TRUE),
  answer("Modell mit 4 Faktoren, weil der p-Wert (p < 0.05) anzeigt, dass die Null-Hypothese, dass beide Modell gleich gut zu den Daten passen, verworfen wird. Wir entscheiden uns daher für das restriktivere Modell mit weniger Parametern; hier: das Modell mit vier Faktoren.")
)

Nun wollen wir uns das fünffaktoirelle Modell noch im Vergleich zum sechsfaktoriellen Modell ansehen.

# Passt auch eines mit 6 Faktor?
six_factor_ML <- ...
six_factor_ML$STATISTIC
six_factor_ML$PVAL # Modell wird durch die Daten nicht verworfen

### ML
six_factor_ML <- fa(dataFR, nfactors = 6, rotate = "oblimin", fm = "ml")
six_factor_ML$STATISTIC
six_factor_ML$PVAL # Modell wird durch die Daten nicht verworfen

### ML
six_factor_ML <- fa(dataFR, nfactors = 6, rotate = "oblimin", fm = "ml")

# For further use in exercises
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")
two_factor_ML <- fa(dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")
one_factor_ML <- fa(dataFR2, nfactors = 1, rotate = "oblimin", fm = "ml")
three_factor_ML <- fa(dataFR2, nfactors = 3, rotate = "oblimin", fm = "ml")
five_factor_ML <- fa(dataFR, nfactors = 5, rotate = "oblimin", fm = "ml")
four_factor_ML <- fa(dataFR, nfactors = 4, rotate = "oblimin", fm = "ml")
six_factor_ML <- fa(dataFR, nfactors = 6, rotate = "oblimin", fm = "ml")

Dem sechsfaktoriellen Modell widersprechen die Daten genauso wenig, wie dem Fünffakoriellen (beide $p>0.05$). Dies war zu erwarten, da wir durch Hinzunahme des sechsten Faktors die Komplexität des Modell erhöht haben, wodurch sich das Modell stärker der konkreten Datenlage annähern kann. Wie in den Folien zur PCA und EFA gesehen, bedeuten mehr Faktoren eine detailgetreuere Abbildung der ursprünglichen Datenlage (siehe auch Eid et al., 2017, Kapitel 25).

# führen Sie eine anova zwischen dem fünf- und dem sechsfaktoriellen Modell durch

# Intro into excercise
data_France <- Big5[Big5$country == "FR", ]
dataFR <- data_France[, -c(1:4)] # entferne demografische Daten und speichere als "dataFR"
dataFR2 <- dataFR[,1:6] # Zunächst wählen wir die ersten 6 Items: E1 bis E3 und N1 bis N3
two_factor <- fa(dataFR2, nfactors = 2, rotate = "varimax")
two_factor_oblimin <- fa(dataFR2, nfactors = 2, rotate = "oblimin")
two_factor_ML <- fa(dataFR2, nfactors = 2, rotate = "oblimin", fm = "ml")
one_factor_ML <- fa(dataFR2, nfactors = 1, rotate = "oblimin", fm = "ml")
three_factor_ML <- fa(dataFR2, nfactors = 3, rotate = "oblimin", fm = "ml")
five_factor_ML <- fa(dataFR, nfactors = 5, rotate = "oblimin", fm = "ml")
four_factor_ML <- fa(dataFR, nfactors = 4, rotate = "oblimin", fm = "ml")
six_factor_ML <- fa(dataFR, nfactors = 6, rotate = "oblimin", fm = "ml")

# Begin grading:
res <- anova(five_factor_ML, six_factor_ML) 
sol2 <- anova(six_factor_ML, five_factor_ML) 

grade_result(
  fail_if(~ (!identical(.result, res) & !identical(.result, sol2)), 'Der Befehl anova muss in diesem Zusammenhang zwei Argumente entgegennehmen: die zwei Objekte der ML-EFAs mit jeweils 5 und 6 Faktoren, welche wir zuvor erstellt haben! Achten Sie entsprechend darauf, dass Sie auch wirklich die ML-EFA-Objekte verwenden.'),
    fail_if(~ (identical(.result, sol2)), 'Das ist an sich RICHTIG, allerdings sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell zuerst übergeben werden. Hier: five_factor_ML, da sonst die df negativ sein können und auch der Chi^2-Wert negativ sein könnte!'),
  pass_if(~ identical(.result, res)),
  correct = 'Sehr gut! Dies ist der Modellvergleich des fünf- und sechsfaktoriellen Modells. Hier sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell zuerst übergeben werden; also: five_factor_ML, da sonst die df negativ sein können und auch der Chi^2-Wert negativ sein könnte!',
  incorrect = 'Leider falsch.',
  glue_correct = '{.correct}',
  glue_incorrect = '{.incorrect} {.message}')

Der $\chi^2$-Differenzwert liegt hier bei r round(anova(five_factor_ML, six_factor_ML)[[4]][2], 3) mit einen zugehörigen p-Wert von r ifelse(anova(five_factor_ML, six_factor_ML)[[5]][2] < 0.01, "< 0.01", paste(round(anova(five_factor_ML, six_factor_ML)[[5]][2], 4))) mit $\Delta df$ = r anova(five_factor_ML, six_factor_ML)[[2]][2]. Mit diesem Test wird geprüft, ob das sparsamere Modell die Daten schlechter abbildet. Die Nullhypothese ist also, dass das sparsamere Modell die Daten genauso gut beschreiben kann, wie das komplexere Modell ($H_0:\Sigma_{5-Fakt.} = \Sigma_{4-Fakt.}$). Da in diesem Fall der p-Wert größer als $.05$ ist, wird diese Nullhypothese nicht verworfen und wir entscheiden uns - Ockhams Rasiermesser folgend (siehe Eid et al., 2017, p. 787) - für das sparsamere Modell.

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 {#AppendixA}

Kürzen des Datensatzes

Falls Sie die Originaldaten auf r fontawesome::fa("download") 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 {#AppendixB}

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 $0$ sind und wir somit diese ignorieren können): $$\begin{pmatrix}E_1\E_2\E_3\N_1\N_2\N_3 \end{pmatrix} = \begin{pmatrix} \lambda_{11} & \lambda_{12}\ \lambda_{21} & \lambda_{22}\ \lambda_{31} & \lambda_{32}\ \lambda_{41} & \lambda_{42}\ \lambda_{51} & \lambda_{52}\ \lambda_{61} & \lambda_{62} \end{pmatrix} \begin{pmatrix}\xi_1\\xi_2 \end{pmatrix} + \begin{pmatrix}\varepsilon_{E_1}\\varepsilon_{E_2}\\varepsilon_{E_3}\\varepsilon_{N_1}\\varepsilon_{N_2}\\varepsilon_{N_3} \end{pmatrix}$$

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 $\xi_1=\xi_\text{Extraversion}$ und $\xi_2=\xi_\text{Neurotizismus}$ nennen. Bennen wir nun die Faktorladungsmatrix als $\Lambda$ und die Korrelationsmatrix der latenten Variablen $\Phi$ (die Diagonaleinträge sind $1$). Die Kovarianz zwischen dem ersten Extraversionitem und dem Extraversionfaktor ist folgendermaßen zu berechnen (wir rechnen hier mit Kovarianzen, da dies im Allgemeinen deutlich einfacher ist als mit Korrelationen zu rechen. Außerdem sind hier alle Variablen standardisiert und somit sind Korrelation und Kovarianz identisch; über die Rechenregeln und die Beziehungen zwischen Korrelation und Kovarianz können sie in Eid, et al. (2017) S. 195-196 und folgend und S.570-571 und folgend nachlesen): $$\begin{align}\mathbf{C}ov[E_1, &\xi_1]\&= \mathbf{C}ov[\lambda_{11}\xi_1 + \lambda_{12}\xi_2+\varepsilon_{E_1}, \xi_1] = \lambda_{11}\mathbf{C}ov[\xi_1, \xi_1] + \lambda_{12}\mathbf{C}ov[\xi_2, \xi_1] +\mathbf{C}ov[\varepsilon_{E_1}, \xi_1] = \lambda_{11}\phi_{11} + \lambda_{12}\phi_{21}\&= \lambda_{11} + \lambda_{12}\phi_{21}\end{align}$$ $\mathbf{C}ov[\varepsilon_{E_1}, \xi_1]=0$ gilt, da die Fehler als unabhängig von allen weiteren Variablen im Modell angenommen werden. Außerdem sind $\phi_{11}=1$ und $\phi_{21}=\phi_{12}$ die Varianz von $\xi_1$ und die Kovarianz/Korrelation zwischen $\xi_1$ und $\xi_2$ und entsprechend Einträge von $\Phi$. Aus dieser Rechnung folgt, dass der erste Eintrag in der Strukturmatrix (an der Stelle 1. Zeile, 1. Spalte) $\lambda_{11} + \lambda_{12}\phi_{21}$ sein muss. Hier ist zu erkennen, dass falls die Korrelation zwischen den latenten Variablen als $0$ angenommen wird (im orthogonalen Fall gilt dann $\phi_{21}=\phi_{12}=0$), dann ist die Strukurmatrix gleich der Faktorladungsmatrix und der erste Eintrag lautet $\lambda_{11}$. Wir schauen uns dies empirisch für das orthogonale Modell (bereits geschätzt in 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
two_factor$Structure[1, 1] # erster Eintrag in der Strukturmatrix

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)
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
two_factor_oblimin$Structure[1, 1] # erster Eintrag in der Strukturmatrix

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 $\hat{\phi}{21}$=r round(two_factor_oblimin$Phi[2, 1], 4) liegt und somit $\hat{\lambda}{12}\hat{\phi}{21}=$ r round(two_factor_oblimin$loadings[1, 2]*two_factor_oblimin$Phi[2, 1], 4) keine große Veränderung zu $\hat{\lambda}{11}$ mit sich bringt. In Matrixschreibweise lässt sich die Strukturmatrix unkompliziert bestimmen. Sie wird durch folgenden Ausdruck berechnet: $$\Lambda\Phi$$ Dies können wir in 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 ersten Sitzung aufmerksam gemacht wurde. Wir berechnen nun das Matrixprodukt für den oblique rotieren Fall:

two_factor_oblimin$loadings[1:6, 1:2] %*% two_factor_oblimin$Phi[1:2, 1:2] # Matrixprodukt
two_factor_oblimin$Structure[1:6, 1:2] # Strukturmatrix

%*% signalisiert R, dass ein Matrixprodukt und keine komponentenweise Mulitplikation durchzuführen ist. Wir haben hier in die Matrizen rein indiziert, damit auch alle Elemente verwendet und angezeigt werden.

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[1:15,1:5] %*% five_factor_ML$Phi[1:5,1:5] # Matrixprodukt
five_factor_ML$Structure[1:15,1:5] # Strukturmatrix

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[1:15,1:5] %*% five_factor_ML$Phi[1:5,1:5] - five_factor_ML$Structure[1:15,1:5]

Da hier nur Nullen herauskommen, scheinen die Ausdrücke identisch zu sein!

Appendix C {#AppendixC}

Prüfen der Voraussetzungen mit der Mahalanobisdistanz

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 einer $\chi^2$-Verteilung vergleichen; wobei $df=p$ und $p=$ Anzahl an Variablen (hier $df=p=15$).

Mahalanobis_Distanz <- mahalanobis(x = dataFR, cov = cov(dataFR), center = apply(X = dataFR, MARGIN = 2, FUN = mean)) # 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 $\chi^2$ Verteilung zu passen. Allerdings sind die Abweichungen auch nicht enorm. Wir verwerfen auf Basis des Histogramms die Normalverteilungsannahme nicht, sollten die Ergebnisse aber trotzdem unter Vorbehalt interpretiert werden.

Die Funktion mahalanobis berechnet die Mahalanobisdistanz pro Proband. 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 berechnen wir mit mean und der apply-Funktion (Man könnte dies natürlich auch "zu Fuß" für jede Variable einzeln machen, die Mittelwert in einem Vektor abspeichern und diesen hier angeben). Die Funktion apply führt auf die Datenmatrix, welche dem Argument X übergeben wird, entweder über die Zeilen MARGIN = 1 oder über die Spalten MARGIN = 2 (hier gewählt) die Funktion aus, welche im Argument FUN angegeben wird. So wird mit FUN = mean der Mittelwert über alle Variablen in dataFR berechnet.

apply(X = dataFR, MARGIN = 2, FUN = mean)

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 $\chi^2$-Verteilung einzuzeichnen). Mit 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 $\chi^2(df=15)$-Verteilung ein: 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.

Literatur

Eid, M., Gollwitzer, M., & Schmitt, M. (2017). Statistik und Forschungsmethoden (5. Auflage, 1. Auflage: 2010). Weinheim: Beltz.

• Blau hinterlegte Autorenangaben führen Sie direkt zur universitätsinternen Ressource.

martscht/PsyMSc1 documentation built on July 11, 2020, 3:45 a.m.