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)
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.
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).
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
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.
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.
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.
# 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.
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.
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.
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!
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!
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.
Eid, M., Gollwitzer, M., & Schmitt, M. (2017). Statistik und Forschungsmethoden (5. Auflage, 1. Auflage: 2010). Weinheim: Beltz.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.