In martscht/PsyBSc7: Übungen für die Praktika im Sommersemester 2020
source('setup.R')
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(car)
library(MASS)
library(ggplot2)
library(lm.beta) # erforderlich für standardiserte Gewichte
data("Schulleistungen", package = "PsyBSc7")
knitr::opts_chunk$set(exercise.checker = gradethis::grade_learnr)
Vorbereitung
Übungs-Datensatz (wie beim letzten Termin)
Eine Stichprobe von 100 Schülerinnen und Schülern hat einen Lese- und einen Mathematiktest, sowie zusätzlich einen allgemeinen Intelligenztest, bearbeitet. Im Datensatz enthalten ist zudem das Geschlecht (female, 0=m, 1=w). Wir laden als erstes die nötigen R-Pakete sowie den Datensatz.
library(car)
library(MASS)
library(ggplot2)
library(lm.beta) # erforderlich für standardiserte Regressionsgewichte
# Datensatz laden
data("Schulleistungen", package = "PsyBSc7")
Sagen Sie mit Hilfe eines Regressionsmodells die Leseleistung durch das Geschlecht und durch die Intelligenz vorher (im Datensatz Schulleistungen) und speichern sie diese Modell unter dem Namen mod ab. Geben Sie das Objekt mod aus.
# Modell zur Vorhersage der Leseleistung
mod <- ...
...
# Verwenden Sie den Befehl
lm(...~...)
# Finale Lösung:
mod <- lm(reading ~ female + IQ, data = Schulleistungen)
mod
# For further use in exercises
mod <- lm(reading ~ female + IQ, data = Schulleistungen)
Im Output sehen wir die Parameterschätzungen unseres Regressionsmodells:
$$\text{Reading}_i=b_0+b_1\text{Female}_i+b_2\text{IQ}_i+\varepsilon_i,$$
für $i=1,\dots,100=:n$. Wir wollen uns die Ergebnisse unserer Regressionsanalyse noch detaillierter anschauen. Dazu können wir wieder die summary Funktion anwenden. Um auch die standardisierten Ergebnisse zu erhalten, verwenden wir die Funktion lm.beta (lm steht hier für lineares Modell und beta für die standardisierten Koeffizienten. Achtung: Häufig werden allerdings auch unstandardisierten Regressionkoeffizienten als $\beta$s bezeichnet, sodass darauf stehts zu achten ist - siehe bspw. Appendix A).
# Verwenden Sie den Befehl
lm.beta(...) # um das geschätzte Modell zu standardisieren
# Verwenden Sie den Befehl
summary(...) # um dem standardisierten Modell die Zusammenfassungsstatistiken zu entlocken
#Finale Lösung
summary(lm.beta(mod))
# For further use
mod <- lm(reading ~ female + IQ, data = Schulleistungen)
# for later use
stdsum <- summary(lm.beta(mod))
An der standardisierten Lösung fällt auf, dass das Interzept $0$ ist. Dies macht völlig Sinn, da beim Standardisieren die Mittelwerte aller Variablen auf $0$ und die Standardabweichungen auf $1$ gesetzt werden. So auch beim Kriterium, der Leseleistung. Damit beträgt der Mittelwert des Kriteriums und der Mittelwert aller Prädiktore jeweils $0$. Das Interzept gibt an, welchen Wert das Kriterium hat, wenn alle Prädiktoren den Wert $0$ aufweisen. Somit muss das Interzept hier genau $0$ betragen. Den Koeffizienten entnehmen wir, dass die Prädiktoren female und IQ signifikante Anteile am Kriterium erklären (Mit einer Irrtumswahrscheinlichkeit von 5% ist der Regressionkoeffizent und damit die lineare Beziehung zwischen dem jeweiligen Prädiktor (z.B. IQ) und der Leseleistung in der Population nicht 0). Der standardisierte Koeffizient des IQ ist hierbei deutlich größer ($\hat{\beta}2=$ r round(stdsum$coefficients[3,2], 2)) als vom Geschlecht ($\hat{\beta}_1=$ r round(stdsum$coefficients[2,2], 2)). Damit ist ersichtlich, dass der IQ mehr Variation am Kriterium Leseleisstug erklärt. Achtung: Diese Information ist den unstandardisierten Koeffizienten nicht zu entnehmen, da ein direkter Vergleich der unstandardisierten Koeffizienten nicht möglich ist. Der unstandardisierte Effekt des Geschlechts auf die Leseleistung beträgt $\hat{b}_1=$ r round(stdsum$coefficients[2,1], 2) und der Effekt des IQ lediglich $\hat{b}_2=$ r round(stdsum$coefficients[3,1], 2). Dies liegt natürlich daran, dass der IQ eine wesentlich größere Streuung ($\sigma{IQ}^2$=r round(var(Schulleistungen$IQ), 2)) aufweist als die dummykodierte Variable female ($\sigma_{female}^2$=r round(var(Schulleistungen$female), 2)). Insgesamt werden r round(stdsum$r.squared*100, 2)% der Variation der Leseleistung durch die beiden Prädiktoren erklärt.
Gegenstand der folgenden Sitzung ist die Prüfung von Voraussetzungen der Regressionsanalyse:
- Linearität (Wiederholung)
- Homoskedastizität (Wiederholung)
- Normalverteilung der Residuen (Wiederholung)
- Multikollinearität
- Identifikation von Ausreißern und einflussreichen Datenpunkten
Im Folgenden werden wir mit den unstandardisierten Modell weiterarbeiten, welches wir im Objekt mod gespeichert haben.
Linearität
Eine grafische Prüfung der partiellen Linearität zwischen den einzelnen Prädiktorvariablen und dem Kriterium kann durch partielle Regressionsplots (engl. Partialplots) erfolgen. Dafür sagen wir in einem Zwischenschritt einen einzelnen Prädiktor durch alle anderen Prädiktoren im Modell vorher und speichern die Residuen, die sich aus dieser Regression ergeben. Diese kennzeichnen den eigenständigen Anteil, den ein Prädiktor nicht mit den anderen Prädiktoren gemein hat. Dann werden die Residuen aus dieser Vorhersage gegen die Residuen des Kriteriums bei Vorhersage durch den jeweiligen Prädiktor dargestellt. Diese Grafiken können Hinweise auf systematische nicht-lineare Zusammenhänge geben, die in der Modellspezifikation nicht berücksichtigt wurden. Die zugehörige R-Funktion des car Pakets heißt avPlots und braucht als Argument lediglich das spezifizierte Regressionsmodell mod.
avPlots(model = ..., pch = 16, lwd = 4)
# Finale Lösung:
avPlots(model = mod, pch = 16, lwd = 4)
Mit Hilfe der Argumente pch=16 und lwd=4 werden die Darstellung der Punkte (ausgefüllt anstatt leer) sowie die Dicke der Linie (vierfache Dicke) manipuliert. Den Achsenbeschriftungen ist zu entnehmen, dass auf der Y-Achse jeweils reading | others dargestellt ist. Die vertikale Linie | steht hierbei für den mathematischen Ausdruck gegeben. Others steht hierbei für alle weiteren (anderen) Prädiktoren im Modell. Dies bedeutet, dass es sich hierbei um die Residuen aus der Regression von reading auf alle anderen Prädiktoren handelt. Bei den unabhängigen Variablen (UV, female, IQ) steht UV | others also jeweils für die jeweilige UV gegeben der anderen UVs im Modell. Somit beschreiben die beiden Plots jeweils die Beziehungen, die die UVs über die anderen UVs im Modell hinaus mit dem Kriterium (AV, abhängige Variable) haben.
Verteilung der Residuen
Homoskedastizität
Die Varianz der Residuen sollte unabhängig von den Ausprägungen der Prädiktoren sein. Dies wird i.d.R. grafisch geprüft, indem die Residuen $e_i$ gegen die vorhergesagten Werte $\hat{y}_i$ geplottet werden: der sogenannte Residuenplot (engl.: residual plot). In diesem Streudiagramm sollten die Residuen gleichmäßig über $\hat{y}_i$ streuen und keine systematischen Trends (linear, quadratisch, auffächernd, o.ä.) erkennbar sein. Prüfung nicht-linearer Zusammenhänge: Die Funktion residualPlots des Pakets car erzeugt separate Streudiagramme für die Residuen in Abhängigkeit von jedem einzelnen Prädiktor $x_j$ und von den vorhergesagten Werten $\hat{y}_i$ ("Fitted Values"); als Input braucht sie das Modell mod. Zusätzlich wird für jeden Plot ein quadratischer Trend eingezeichnet und auf Signifikanz getestet, wodurch eine zusätzliche Prüfung auf nicht-lineare Effekte erfolgt. Sind diese Test nicht signifikant, ist davon auszugehen, dass diese Effekte nicht vorliegen und die Voraussetzungen nicht verletzt sind.
# Grafisch (+ Test auf Nicht-Linearität)
residualPlots(model = ..., pch = 16)
# Finale Lösung:
residualPlots(model = mod, pch = 16)
Prüfung linearer Zusammenhänge: Die Funktion ncvTest (als Input braucht sie mod) prüft, ob die Varianz der Residuen bedeutsam linear (!) mit den vorhergesagten Werten zusammenhängt. Wird dieser Test signifikant, ist die Annahme der Homoskedaszidität verletzt (Wikipedia: Breusch-Pagan). Wird er nicht signifikant, kann dennoch eine Verletzung vorliegen, z.B. ein nicht-linearer Zusammenhang.
# Test For Non-Constant Error Variance
ncvTest(model = ...)
# Finale Lösung:
ncvTest(mod)
Normalverteilung
Voraussetzung für die Signifikanztests im Kontext der linearen Regression ist die Normalverteilung der Residuen. Auch diese Annahme wird i.d.R. grafisch geprüft. Hierfür bietet sich zum einen ein Histogramm der Residuen an, zum anderen ein Q-Q-Diagramm (oder "Quantile-Quantile-Plot"). Zusätzlich zur grafischen Darstellung erlaubt z.B. der Shapiro-Wilk-Test auf Normalität (shapiro.test) einen Signifkanztest der Nullhypothese, dass die Residuen normalverteilt sind. Auch der Kolmogorov-Smirnov (ks.test) Test, welcher eine deskriptive mit einer theoretischen Verteilung vergleicht, wäre denkbar.
Wir beginnen mit der Vorbereitung der Daten. Wir wollen die studentisierten Residuen grafisch mit dem R-Paket ggplot2darstellen. Der Befehl studres aus dem MASS Paket speichert die Residuen aus einem lm Objekt, also aus dem definierten Modell mod, und studentisiert diese. Das Studentisieren der Residuen bezeichnet eine Art der Standardisierung, sodass anschließend der Mittelwert 0 und die Varianz 1 ist (somit lassen sich solche Plots immer gleich interpretieren). Speichern Sie die Residuen aus mod als res ab und erzeugen Sie daraus einen data.frame, den Sie unter dem Namen df_res speichern. Diesen wollen wir später in den Grafiken mit ggplot2 verwenden.
# Verwenden Sie die Befehle
studres(...)
# Verwenden Sie die Befehle
res <- studres(...)
# Verwenden Sie die Befehle
res <- studres(...)
# und
df_res <- data.frame(...)
# Verwenden Sie die Befehle
res <- studres(...)
# und
df_res <- data.frame(...)
# geben Sie das Objekt des Data.frames zurück (sehen Sie sich dieses an)
df_res
# Intro into excercise
mod <- lm(reading ~ female + IQ, data = Schulleistungen)
# Begin grading:
res <- studres(mod) # Studentisierte Residuen als Objekt speichern
df_res <- data.frame(res) # als Data.Frame für ggplot
res <- df_res
grade_result(
fail_if(~ (!identical(.result, res)), 'Haben Sie die richtigen Input-Argumente verwendet? Sie wollen die studentisierten Resiuden abspeichern - als data.frame. Geben Sie diese Objekt aus (indem Sie den Namen des Objekts ohne Zuordnung als letztes schreiben).'),
pass_if(~ identical(.result, res)),
correct = 'Sehr gut! Sie haben erfolgreich die Residuen als data.frame abgespeichert.',
incorrect = 'Leider falsch.',
glue_correct = '{.correct}',
glue_incorrect = '{.incorrect} {.message}')
res <- studres(mod) # Studentisierte Residuen als Objekt speichern
df_res <- data.frame(res) # als Data.Frame für ggplot
# For further use in exercises
mod <- lm(reading ~ female + IQ, data = Schulleistungen)
res <- studres(mod) # Studentisierte Residuen als Objekt speichern
df_res <- data.frame(res) # als Data.Frame für ggplot
Im folgenden Block sehen wir den Code für ein Histogramm in ggplot2-Notation. Hier sind einige Zusatzeinstellungen gewählt, die das Histogramm optisch aufbereiten. Vervollständigen Sie den Code mit Hilfe der Tipps.
# Grafisch: Histogramm mit Normalverteilungskurve
ggplot(data = ???, aes(x = ??)) +
geom_histogram(aes(y =..density..),
bins = 15, # Wie viele Balken sollen gezeichnet werden?
colour = "blue", # Welche Farbe sollen die Linien der Balken haben?
fill = "skyblue") + # Wie sollen die Balken gefüllt sein?
stat_function(fun = dnorm, args = list(mean = mean(???), sd = sd(???)), col = "darkblue") + # Füge die Normalverteilungsdiche hinzu und nutze den empirischen Mittelwert und die empirische Standardabweichung, wähle dunkelblau als Linienfarbe
labs(title = "Histogramm der Residuen mit Normalverteilungsdichte", x = "Residuen") # Füge eigenen Titel und Achsenbeschriftung hinzu
# an
data # wird in ggplot immer der data.frame mit den Daten übergeben
# der Dataframe enthält in unserem Fall eine Variable, welche die Residuen enthält. Den Namen haben wir selbst vergeben. Diesen können Sit mit
colnames()
# in Erfahrung bringen (welcher die Namen der Spalten ausgibt)
# Auf die x-Achse in einem Histogramm wird immer die Variable selbst abgetragen. Diese müssen sie an
aes(x = ...) # übergeben
# mit
stat_function(fun = dnorm, args = list(...)) # übergeben wir eine statistische Funktion,
# welche dazu geplottet werden soll.
# Hier ist dies die Dichte der Normalverteilung "dnorm". Diese Funktion wird normalerweise mit
dnorm(x = ..., mean = ..., sd = ...) # angesprochen. Wir können mit dieser die Dichte der Normalverteilung
# am Punkt x berechnen. Um welche Normalverteilung es sich handelt, geben wir mit "mean" und "sd" an.
# Die Normalverteilung ist eindeutig durch Mittelwert und Standardabweichung definiert.
# mit
stat_function(fun = dnorm, args = list(...)) # übergeben wir eine statistische Funktion,
# welche dazu geplottet werden soll.
# Hier ist dies die Dichte der Normalverteilung "dnorm". Diese Funktion wird normalerweise mit
dnorm(x = ..., mean = ..., sd = ...) # angesprochen. Wir können mit dieser die Dichte der Normalverteilung
# am Punkt x berechnen. Um welche Normalverteilung es sich handelt, geben wir mit "mean" und "sd" an.
# Die Normalverteilung ist eindeutig durch Mittelwert und Standardabweichung definiert.
# Die Argumente, die dnorm normalerweise erhält, können wir in args als Liste in stat_function abgeben
# Hier brauchen wir allerdings "x" nicht, da das übrigbleibenden (wenn wir mean und sd angeben) Argument
# für den Plot manipuliert wird (von ggplot2). Somit erhalten wir für die Ausprägungen entlang der
# x-Achse die Dichte der Normalverteilung (und zwar der Residuen!).
# Finale Lösung:
ggplot(data = df_res, aes(x = res)) +
geom_histogram(aes(y =..density..),
bins = 15, # Wie viele Balken sollen gezeichnet werden?
colour = "blue", # Welche Farbe sollen die Linien der Balken haben?
fill = "skyblue") + # Wie sollen die Balken gefüllt sein?
stat_function(fun = dnorm, args = list(mean = mean(res), sd = sd(res)), col = "darkblue") + # Füge die Normalverteilungsdiche "dnorm" hinzu und nutze den empirischen Mittelwert und die empirische Standardabweichung "args = list(mean = mean(res), sd = sd(res))", wähle dunkelblau als Linienfarbe
labs(title = "Histogramm der Residuen mit Normalverteilungsdichte", x = "Residuen") # Füge eigenen Titel und Achsenbeschriftung hinzu
Nutzen wir nur die Defaulteinstellung des Histogramms (bis auf bins = 15 - für die Vergleichbarkeit der beiden Grafiken), sieht es so aus:
# hier nochmal nur das Nötigste:
ggplot(data = df_res, aes(x = res)) +
geom_histogram(aes(y =..density..), bins = 15) +
stat_function(fun = dnorm, args = list(mean = mean(res), sd = sd(res)))
Ähnliche Informationen sind aus dem Q-Q-Plot zu entnehmen. Dafür kann der Befehl qqplot aus dem car Paket genutzt werden. Hier kann die Verteilung, gegen die geprüft werden soll, direkt mit distribution festgelegt werden kann.
# Grafisch: Q-Q-Diagramm mit der car Funktion qqPlot
qqPlot(..., pch = 16, distribution = ...)
qqplot() # brauch als Argument das Regressionmodell "mod".
# Die Normalverteilung heißt in R "norm".
# Finale Lösung:
qqPlot(mod, pch = 16, distribution = "norm")
Nun wollen wir die Voraussetzung zusätzlich mittels statistischer Tests überprüfen. Diese testen im Grunde genommen die Null-Hypothese: $H_0:$ "Normalverteilung liegt vor". (Wir wissen allerdings, dass wir Hypothesen nur verwerfen können und nicht wirklich annehmen, weshalb diese Tests immer nur unter Vorbehalt interpretiert werden sollten - Aussagen wie etwa: "$H_0$ nicht verworfen; daraus folgt, dass in der Population Normalverteilung herrscht", sind nicht zulässig. Wird die $H_0$ nicht verworfen, nehmen wir nur weiterhin an, dass Normalverteilung vorliegt, ein Beweis ist dies jedoch nicht).
Die Funktion shapiro.test und ks.test benötigen die Variable, die auf Normalverteilung geprüft werden soll, als Argument; hier also res.
# Test auf Abwweichung von der Normalverteilung mit dem Shapiro Test
shapiro.test(...)
# Finale Lösung:
shapiro.test(res)
# Test auf Abwweichung von der Normalverteilung mit dem Kolmogorov-Smirnov Test
ks.test(..., "...", mean(...), sd(...))
# Die Argumente in
ks.test # sind wie folgt angeordnet:
# (Variable, "Verteilung (Wahrscheinlichkeitswerte)", Argument1_der_Verteilung, Argument2_der_Verteilung)
# Variable = res
# "Verteilung (Wahrscheinlichkeitswerte)" = "pnorm", die Verteilungswerte der Normalverteilung mit den Argumenten q=Quantil, sowie Mittelwert und Standardabweichung
# Argument1_der_Verteilung = Mittelwert von res
# Argument2_der_Verteilung = Standardabweichung von res
# Finale Lösung:
ks.test(res, "pnorm", mean(res), sd(res))
Für unser Modell zeigen alle Ergebnisse übereinstimmend, dass die Residuen normalverteilt sind (bzw. lehnen diese Hypothese nicht ab: Shapiro-Wilk Test und Kolmogorov-Smirnov Tests gegen die Normalverteilung). Durch die Angabe von distribution = "norm als Argument für die Funktion qqPlot des car-Pakets erzeugt den Vergleich der vorliegenden Variable gegen eine Normalverteilung (auch der Vergleich zu anderen Verteilungen wäre möglich). Dem Kolmogorov-Smirnov Test gegen die Normalverteilung ks.test muss zunächst als erstes Argument res, die Variable, deren Verteilung untersucht werden soll, übergeben werden. "pnorm" gibt an, dass die empirische Verteilung unserer Daten res mit der Normalverteilung verglichen werden soll (wird hier eine andere Verteilung gewählt, so kann gegen diese geprüft werden; auch können zwei empirisiche Verteilungen auf Gleichheit untersucht werden). Anschließend geben wir die verteilungsspezifischen Parameter; hier mean(res) (der Mittelwert) sowie sd(res) (die Standardabweichung) an, da die Normalverteilung genau durch diese beiden Parameter bestimmt ist (hier könnten auch theoretisch angenommene Parameter gewählt werden, z.B. 0 und 1 für eine Mittelwert von 0 und einer Varianz von 1 und es kann geprüft werden, ob die beobachteten Daten dieser Verteilung widersprechen).
Wäre der Shapiro-Wilk-Normality Tests oder der Kolmogorov-Smirnov Tests gegen die Normalverteilung signifikant auf dem 5%-Niveau, würde dies dafür sprechen, dass "mit einer Irrtumswahrscheinlichkeit von 5% die Null-Hypothese auf Gleichheit der Normalverteilung mit unserer empirischen Verteilung verworfen werden muss, da in der Population die Abweichung der empirischen zur theoretischen Normalverteilung mit der Irrtumswahrscheinlichkeit von 5% ungleich 0 ist".
Schiefe Verteilungen, welche sich im Histogramm oder im Q-Q-Plot wiederspiegeln, sowie zu signifikanten Shapiro-Wilk-Normality Tests oder Kolmogorov-Smirnov Tests gegen die Normalverteilung führen können, könnten Indizien für nichtlineare (z.B. quadratische) Beziehungen im Datensatz sein. Weitere Analysen wären hierfür nötig (siehe hierzu PsyBSc7::Sitzung_8()).
Multikollinearität
Multikollinearität ist ein potenzielles Problem der multiplen Regressionsanalyse und liegt vor, wenn zwei oder mehrere Prädiktoren hoch miteinander korrelieren. Hohe Multikollinearität
- schränkt die mögliche multiple Korrelation ein, da die Prädiktoren redundant sind und überlappende Varianzanteile in $y$ erklären.
- erschwert die Identifikation von bedeutsamen Prädiktoren, da deren Effekte untereinander konfundiert sind (die Effekte können schwer voneinander getrennt werden).
- bewirkt eine Erhöung der Standardfehler der Regressionkoeffizienten (der Standardfehler ist die Standardabweichung zu der Varianz der Regressionskoeffizienten bei wiederholter Stichprobenziehung und Schätzung). Dies bedeutet, dass die Schätzungen der Regressionsparameter instabil, und damit weniger verlässlich, werden.
Für weiterführende Informationen zur Instabilität und Standardfehlern siehe Appendix A.
Multikollinearität kann durch Inspektion der bivariaten Zusammenhänge (Korrelationsmatrix) der Prädiktoren $x_j$ untersucht werden, dies kann aber nicht alle Formen von Multikollinearität aufdecken. Darüber hinaus ist die Berechung der sogennanten Toleranz und des Varianzinflationsfaktors (VIF) für jeden Prädiktor möglich. Hierfür wird für jeden Prädiktor $x_j$ der Varianzanteil $R_j^2$ berechnet, der durch Vorhersage von $x_j$ durch alle anderen Prädiktoren in der Regression erklärt wird. Toleranz und VIF sind wie folgt definiert:
- $T_j = 1-R^2_j = \frac{1}{VIF_j}$
- $VIF = \frac{1}{1-R^2_j} = \frac{1}{T_j}$
Offensichtlich genügt eine der beiden Statistiken, da sie vollständig ineinander überführbar und damit redundant sind. Empfehlungen als Grenzwert für Kollinearitätsprobleme sind z. B. $VIF_j>10$ ($T_j<0.1$) (siehe Eid, Gollwitzer, & Schmitt, 2017, S. 712 und folgend). Die Varianzinflationsfaktoren der Prädiktoren im Modell können mit der Funktion vif des car-Paktes bestimmt werden, der Toleranzwert als Kehrwert des VIFs.
# Korrelation der Prädiktoren
cor(Schulleistungen$female, Schulleistungen$IQ)
Die Prädiktoren sind nur schwach negativ korreliert. Wir schauen uns trotzdem den VIF und die Toleranz an. Dazu übergeben wir wieder das definierte Regressionsmodell an vif.
# Varianzinflationsfaktoren
vif(...)
# Finale Lösung:
vif(mod)
Berechnen Sie die Toleranz.
# Verwenden Sie den Befehl
vif(...) # von zuvor und berechnen Sie den Kehrwert
# Intro into excercise
mod <- lm(reading ~ female + IQ, data = Schulleistungen)
# Begin grading:
res <- 1/vif(mod)
grade_result(
fail_if(~ (!identical(.result, res)), 'Haben Sie die richtigen Input-Argumente verwendet? Sie wollen die Toleranz als Kehrwert des VIF bestimmen.'),
pass_if(~ identical(.result, res)),
correct = 'Sehr gut! Sie haben erfolgreich die Toleranz aus dem VIF berechnet.',
incorrect = 'Leider falsch.',
glue_correct = '{.correct}',
glue_incorrect = '{.incorrect} {.message}')
Für unser Modell wird ersichtlich, dass die Prädiktoren praktisch unkorreliert sind und dementsprechend kein Multikollinearitätsproblem vorliegt. Unabhängigkeit folgt hieraus allerdings nicht, da nichtlineare Beziehungen zwischen den Variablen bestehen könnten, die durch diese Indizes nicht abgebildet werden.
Identifikation von Ausreißern und einflussreichen Datenpunkten
Die Plausibilität unserer Daten ist enorm wichtig. Aus diesem Grund sollten Aureißer oder einflussreiche Datenpunkte analysiert werden. Diese können bspw. durch Eingabefehler entstehen (Alter von 211 Jahren anstatt 21) oder es sind seltene Fälle (hochintelligentes Kind in einer Normstichprobe), welche so in natürlicher Weise (aber mit sehr geringer Häufigkeit) auftreten können. Es muss dann entschieden werden, ob Ausreißer die Repräsentativität der Stichprobe gefährden und ob diese besser ausgeschlossen werden sollten.
Hebelwerte $h_j$ erlauben die Identifikation von Ausreißern aus der gemeinsamen Verteilung der unabhängigen Variablen, d.h. einzelne Fälle, die weit entfernt vom Mittelwert der gemeinsamen Verteilung der unabhängigen Variablen liegen und somit einen starken Einfluss auf die Regressionsgewichte haben können. Diese werden mit der Funktion hatvalues ermittelt. Kriterien zur Beurteilung der Hebelwerte variieren, so werden von Eid et al. (2017, S. 707 und folgend) Grenzen von $2\cdot k / n$ für große und $3\cdot k / n$ für kleine Stichproben vorgeschlagen, in den Vorlesungsfolien werden Werte von $4/n$ als auffällig eingestuft (hierbei ist $k$ die Anzahl an Prädiktoren und $n$ die Anzahl der Beobachtungen). Alternativ zu einem festen Cut-Off-Kriterium kann die Verteilung der Hebelwerte inspiziert und diejenigen Werte kritisch betrachtet werden, die aus der Verteilung ausreißen. Die Funktion hatvalues erzeugt die Hebelwerte aus einem Regression-Objekt. Wir wollen diese als Histogramm darstellen.
# Hebelwerte
n <- length(residuals(...))
h <- hatvalues(...) # Hebelwerte
df_h <- data.frame(h) # als Data.Frame für ggplot
ggplot(data = df_h, aes(x = h)) +
geom_histogram(aes(y =..density..), bins = 15)+
geom_vline(xintercept = 4/n, col = "red") # Cut-off bei 4/n
# Finale Lösung
n <- length(residuals(mod))
h <- hatvalues(mod) # Hebelwerte
df_h <- data.frame(h) # als Data.Frame für ggplot
ggplot(data = df_h, aes(x = h)) +
geom_histogram(aes(y =..density..), bins = 15)+
geom_vline(xintercept = 4/n, col = "red") # Cut-off bei 4/n
Cook's Distanz $CD_i$ gibt eine Schätzung, wie stark sich die Regressionsgewichte verändern, wenn eine Person $i$ aus dem Datensatz entfernt wird. Fälle, derem Elimination zu einer deutlichen Veränderung der Ergebnisse führen würden, sollten kritisch geprüft werden. Als einfache Daumenregel gilt, dass $CD_i>1$ auf einen einflussreichen Datenpunkt hinweist. Cook's Distanz kann mit der Funktion cooks.distance ermittelt werden.
# Cooks Distanz
CD <- cooks.distance(...) # Cooks Distanz
df_CD <- data.frame(CD) # als Data.Frame für ggplot
ggplot(data = df_CD, aes(x = CD)) +
geom_histogram(aes(y =..density..), bins = 15)+
geom_vline(xintercept = 1, col = "red") # Cut-Off bei 1
# Cooks Distanz
CD <- cooks.distance(mod) # Cooks Distanz
df_CD <- data.frame(CD) # als Data.Frame für ggplot
ggplot(data = df_CD, aes(x = CD)) +
geom_histogram(aes(y =..density..), bins = 15)+
geom_vline(xintercept = 1, col = "red") # Cut-Off bei 1
Die Funktion influencePlot des car-Paktes erzeugt ein "Blasendiagramm" zur simultanen grafischen Darstellung von Hebelwerten (auf der x-Achse), studentisierten Residuen (auf der y-Achse) und Cooks Distanz (als Größe der Blasen). Vertikale Bezugslinien markieren das Doppelte und Dreifache des durchschnittlichen Hebelwertes, horizontale Bezugslinien die Werte -2, 0 und 2 auf der Skala der studentisierten Residuen. Fälle, die nach einem der drei Kriterien als Ausreißer identifiziert werden, werden im Streudiagramm durch ihre Zeilennummer gekennzeichnet. Diese Zeilennummern können verwendet werden, um sich die Daten der auffälligen Fälle anzeigen zu lassen. Sie werden durch InfPlot ausgegeben werden. Auf diese kann durch as.numeric(row.names(InfPlot)) zugegriffen werden.
# Blasendiagramm mit Hebelwerten, studentisierten Residuen und Cooks Distanz
# In "IDs" werden die Zeilennummern der auffälligen Fälle gespeichert,
# welche gleichzeitig als Zahlen im Blasendiagramm ausgegeben werden
InfPlot <- influencePlot(mod)
IDs <- as.numeric(row.names(InfPlot))
# Werte der identifizierten Fälle
InfPlot
InfPlot <- influencePlot(mod)
IDs <- as.numeric(row.names(InfPlot))
Schauen wir uns die möglichen Ausreißer an und standardisieren die Ergebnisse für eine bessere Interpretierbarkeit.
# Rohdaten der auffälligen Fälle (gerundet für bessere Übersichtlichkeit)
round(Schulleistungen[IDs,],2)
# z-Standardisierte Werte der auffälligen Fälle
round(scale(Schulleistungen)[IDs,],2)
Die Funktion scale z-standardisiert den Datensatz, mit Hilfe von [IDs,], werden die entsprechenden Zeilen der Ausreißer aus dem Datensatz ausgegeben und anschließend auf 2 Nachkommastellen gerundet. Mit Hilfe der z-standardisieren Ergebnisse lassen sich Ausreißer hinsichtlich ihrer Ausprägungen einordnen:
Interpretation
Was ist an den fünf identifizierten Fällen konkret auffällig?
- Fall 6: recht niedrige Lesekompetenz bei gleichzeitig durchschnittlichem bis überdurchschnittlichem IQ
- Fall 9: Sehr hohe Werte in IQ, Lesen & Mathe
- Fall 33: Unterdurchschnittliche Leseleistung "trotz" eher durchschnittlicher Intelligenz
- Fall 80: Sehr niedriger IQ, gleichzeitig durchschnittliche bis überdurchschnittliche Lesekompetenz
- Fall 99: Sehr niedrige Werte in IQ, Lesekompetenz und Mathematik
Die Entscheidung, ob Ausreißer oder auffällige Datenpunkte aus Analysen ausgeschlossen werden, ist schwierig und kann nicht pauschal beantwortet werden. Im vorliegenden Fall wäre z.B. zu überlegen, ob die Intelligenztestwerte der Fälle 80 und 99, die im Bereich von Lernbehinderung oder sogar geistiger Behinderung liegen, in einer Stichprobe von Schülern/innen aus Regelschulen als glaubwürdige Messungen interpretiert werden können oder als Hinweise auf mangelndes Commitment bei der Beantwortung.
Aufgaben
Die Aufgaben zu dieser Sitzung finden Sie im Appendix B.
Appendix A {#AppendixA}
Multikollinearität und Standardfehler
Im Folgenden stehen $\beta$s für unstandardisierte Regressionskoeffizienten.
Für eine einfache Regressionsgleichung mit $$Y_i=\beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + \varepsilon_i$$
kann die selbe Gleichung auch in Matrixnotation formuliert werden $$Y = X\beta + \varepsilon.$$ Hier ist $X$ die Systemmatrix, welche die Zeilenvektoren $X_i=(1, x_{i1}, x_{i2})$ enthält. Die Standardfehler, welche die Streuung der Parameter $\beta:=(\beta_0,\beta_1,\beta_2)$ beschreiben, lassen sich wie folgt ermitteln. Wir bestimmen zunächst die Matrix $I$ wie folgt
$$I:=(X'X)^{-1}\hat{\sigma}^2_e,$$
wobei $\hat{\sigma}^2_e$ die Residualvarianz unserer Regressionanalyse ist (also der nicht-erklärte Anteil an der Varianz von $Y$). Aus der Matrix $I$ erhalten wir die Standardfehler sehr einfach: Sie stehen im Quadrat auf der Diagonale. Das heißt, die Standardfehler sind $SE(\beta)=\sqrt{\text{diag}(I)}$ (Wir nehmen mit $\text{diag}$ die Diagonalelemente aus $I$ und ziehen aus diesen jeweils die Wurzel: der erste Eintrag ist der $SE$ von $\beta_0$; also $SE(\beta_0)=\sqrt{I_{11}}$; der zweite von $\beta_1$;$SE(\beta_1)=\sqrt{I_{22}}$; usw.). Was hat das nun mit der Kollinearität zu tun? Wir wissen, dass in $X'X$ die Information über die Kovariation im Datensatz steckt (dafür muss nur noch durch die Stichprobengröße geteilt werden und das Vektorprodukt der Mittelwerte abgezogen werden; damit wir eine Zentrierung um die Mittelwert sowie eine Normierung an der Stichprobengröße vorgenommen). Beispielsweise lässt sich die empirische Kovarianzmatrix $S$ zweier Variablen $z_1$ und $z_2$ sehr einfach bestimmen mit $Z:=(z_1, z_2)$:
$$ S=\frac{1}{n}Z'Z - \begin{pmatrix}\overline{z}_1\\overline{z}_2 \end{pmatrix}\begin{pmatrix}\overline{z}_1&\overline{z}_2 \end{pmatrix}.$$
Weitere Informationen hierzu (Kovarianzmatrix und Standardfehler) sind im Appendix B (sowie auch in einigen Kapiteln von Eid et al. (2017) Unterpunkt 5.2-5.4 bzw. ab Seite 1058) nachzulesen.
Insgesamt bedeutet dies, dass die Standardfehler von der Inversen der Kovarianzmatrix unserer Daten sowie von der Residualvarianz abhängen. Sie sind also groß, wenn die Residualvarianz groß ist (damit ist die Vorhersage von $Y$ schlecht) oder wenn die Inverse der Kovarianzmatrix groß ist (also wenn die Variablen stark redundant sind und somit hoch mit einander korrelieren). Nehmen wir dazu der Einfachheit halber an, dass $\hat{\sigma}_e^2=1$ (es geht hier nur um eine numerische Präsentation der Effekte, nicht um ein sinnvolles Modell) sowie $n = 100$ (Stichprobengröße). Zusätzlich gehen wir von zentrierten Variablen (Mittelwert von 0) aus. Dann lässt sich aus $X'X$ durch Division durch $100$ die Kovarianzmatrix der Variablen bestimmen. Wir gucken uns drei Fälle an mit
Fall 1: $X'X=\begin{pmatrix}100&0&0\0&100&0\0&0&100\end{pmatrix}.\qquad$
Fall 2: $X'X=\begin{pmatrix}100&0&0\0&100&99\0&99&100\end{pmatrix} \qquad$ und
Fall 3: $X'X=\begin{pmatrix}100&0&0\0&100&100\0&100&100\end{pmatrix}$.
Hierbei ist zu beachten, dass $X$ die Systemmatrix ist, welche auch die $1$ des Interzepts enthält. Natürlich ist eine Variable von einer Konstanten unabhängig, weswegen die erste Zeile und Spalte von $X'X$ jeweils der Vektor $(100, 0, 0)$ ist. Die zugehörigen Korrelationsmatrizen können durch Divison durch 100 berechnet werden (da wir zentrierte Variablen haben, die Stichprobengröße gleich 100 ist und die Varianzen der Variablen gerade 100 sind!). Wir betrachten nur die Minormatrizen, aus welchen die 1. Zeile und die 1. Spalte gestrichen wurden. Diese teilen wir durch 100 und erhalten die Korrelationsmatrix der Variablen:
Fall 1: $\Sigma_1=\begin{pmatrix}1&0\0&1\end{pmatrix}.\qquad$
Fall 2: $\Sigma_2=\begin{pmatrix}1&.99\.99&1\end{pmatrix} \qquad$ und
Fall 3: $\Sigma_3=\begin{pmatrix}1&1\1&1\end{pmatrix}$.
Im Fall 1 sind die zwei Variablen unkorreliert. Die Inverse ist leicht zu bilden.
XX_1 <- matrix(c(100,0,0,
0,100,0,
0,0,100),3,3)
XX_1 # Die Matrix X'X im Fall 1
I_1 <- solve(XX_1)*1 # I (*1 wegen Residualvarianz = 1)
I_1
sqrt(diag(I_1)) # Wurzel aus den Diagonalelementen der Inverse = SE, wenn sigma_e^2=1
Die Standardfehler sind nicht sehr groß: alle liegen bei $0.1$.
Im Fall 2 sind die zwei Variablen fast perfekt (zu $.99$) korreliert - es liegt hohe Multikollinearität vor. Die Inverse ist noch zu bilden. Die Standardfehler sind deutlich erhöht im Vergleich zu Fall 1.
XX_2 <- matrix(c(100,0,0,
0,100,99,
0,99,100),3,3)
XX_2 # Die Matrix X'X im Fall 2
I_2 <- solve(XX_2)*1 # I (*1 wegen Residualvarianz = 1)
I_2
sqrt(diag(I_2)) # SEs im Fall 2
sqrt(diag(I_1)) # SEs im Fall 1
Die Standardfehler des Fall 2 sind sehr groß im Vergleich zu Fall 1 (mehr als sieben Mal so groß); nur der Standardfehler des Interzept bleibt gleich. Die Determinante von $X'X$ in Fall 2 liegt deutlich näher an $0$ im Vergleich zu Fall 1; hier: $10^6$.
det(XX_2) # Determinante Fall 2
det(XX_1) # Determinante Fall 1
Im Fall 3 sind die zwei Variablen perfekt korreliert - es liegt perfekte Multikollinearität vor. Die Inverse kann nicht gebildet werden (da $\text{det}(X'X) = 0$). Die Standardfehler können nicht berechnet werden. Eine Fehlermeldung wird ausgegeben.
XX_3 <- matrix(c(100,0,0,
0,100,100,
0,100,100),3,3)
XX_3 # Die Matrix X'X im Fall 3
det(XX_3) # Determinante on X'X im Fall 3
I_3 <- solve(XX_3)*1 # I (*1 wegen Residualvarianz = 1)
I_3
sqrt(diag(I_3)) # Wurzel aus den Diagonalelementen der Inverse = SE, wenn sigma_e^2=1
# hier wird eine Fehlermeldung ausgegeben, wodurch der Code nicht ausführbar ist und I_3 nicht gebildet werden kann:
# Error in solve.default(XX_3) :
# Lapack routine dgesv: system is exactly singular: U[2,2] = 0
Der VIF bzw. die Toleranz quantifizieren die Korrelation zwischen den beiden Variablen. Der VIF wäre in diesen Analysen im 1. Fall für beide Variablen 1, im 2. Fall für beide Variabeln 50.25 und im 3. Fall nicht berechenbar (bzw. $\infty$). Entsprechend wäre die Toleranz im 1. Fall 1 und 1, im 2. Fall 0.02 und 0.02 sowie im 3. Fall 0 und 0.
Dieser Exkurs zeigt, wie sich die Multikolinearität auf die Standardfehler und damit auf die Präzision der Parameterschätzung auswirkt. Inhaltlich bedeutet dies, dass die Prädiktoren redundant sind und nicht klar ist, welchem Prädiktor die Effekte zugeschrieben werden können.
Die Matrix $I$ ist im Zusammenhang mit der Maximum-Likelihood-Schätzung die Inverse der Fischer-Information und enthält die Informationen der Kovariationen der Parameterschätzer (diese Informationen enthält sie hier im Übrigen auch!).
Appendix B {#AppendixB}
Aufgaben
Datensatz (wie beim letzten Termin)
Der Beispieldatensatz enthält Daten zu Lesekompetenz aus der deutschen Stichprobe der PISA-Erhebung in Deutschland 2009.
library(car)
library(MASS)
library(lm.beta) # erforderlich für standardiserte Gewichte
# Datensatz laden
data("PISA2009", package = "PsyBSc7")
Im Einzelnen enthält der Datensatz die folgenden Variablen:
- Klassenstufe (
Grade)
- Alter in Jahren (
Age)
- Geschlecht (
Female, 0=m / 1=w)
- Lesekompetenz (
Reading)
- Lesefreude (
JoyRead)
- Lernzeit in Minuten für Deutsch (
LearnMins)
- Sozialstatus (
HISEI, "highest international socio-economic index of occupational status")
- Fragebogen-Score für kulturelle Besitztümer zu Hause (
CultPoss, z. B. klassische Literatur, Kunstwerke)
- Anzahl Bücher zu Hause (
Books, eigentlich eine ordinale Ratingskala)
- Anzahl Fernseher zu Hause (
TVs)
- Anzahl Computer zu Hause (
Computers)
- Anzahl Autos zu Hause (
Cars)
- Migrationshintergrund (
MigHintergrund, 0=beide Eltern in D geboren, 1=min. 1 Elternteil im Ausland geboren)
- Bildungsabschluss des Vaters (
FatherEdu, International Standard Classification of Education)
- Bildungsabschluss der Mutter (
MotherEdu, International Standard Classification of Education)
Vorhersage von Lesekompetenz mit individuellen Merkmalen der Schüler/innen
In der letzten Sitzung hatten wir herausgefunden, dass der Sozialstatus (HISEI), der Bildungsabschluss der Mutter (MotherEdu) und die Zahl der Bücher zu Hause (Books) bedeutsame Prädiktoren für die Lesekompetenz der Schüler/innen sind. Berechnen Sie dieses Modell erneut.
# Berechnung des Modells und Ausgabe der Ergebnisse
m1 <- lm(Reading ~ HISEI + MotherEdu + Books, data = PISA2009)
summary(lm.beta(m1))
Aufgabe 1. Prüfen Sie die Annahme der Homoskedastizität
a) durch Inspektion der Residuenplots und
b) den Breusch–Pagan Test auf nicht-konstante Varianz.
Aufgabe 2. Prüfen Sie die Annahme der Linearität
a) durch Inspektion der partiellen Regressionsplots und
b) durch den Test für quadratische Trends, der mit der Funktion residualPlots erzeugt wurde.
Aufgabe 3. Püfen Sie die Normalverteilung der Residuen
a) grafisch (Histogramm, Q-Q-Plot) und
b) über den Shapiro-Wilks-Test sowie den Kolmogorov-Smirnov-Test.
Aufgabe 4. Prüfen Sie die Multikollinearität der Prädiktoren im Modell
a) durch Inspektion der bivariaten Korrelationen und
b) durch Berechnung von VIF und Toleranz.
Aufgabe 5. Ausreißer und einflussreiche Datenpunkte
a) Inspizieren Sie die Verteilung der Hebelwerte $h_i$.
b) Inspizieren Sie die Verteilung von Cooks Distanz $CD_i$.
c) Sichten Sie die Rohdaten und standardisierten Werte auffälliger Fälle (verwenden Sie die Funktion influencePlot).
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.
$$ Diese Sitzung basiert zum Teil auf Unterlagen von Prof. Johannes Hartig aus dem SoSe 2019.*
martscht/PsyBSc7 documentation built on Sept. 1, 2020, 10:50 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
source('setup.R') 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(car) library(MASS) library(ggplot2) library(lm.beta) # erforderlich für standardiserte Gewichte data("Schulleistungen", package = "PsyBSc7") knitr::opts_chunk$set(exercise.checker = gradethis::grade_learnr)
Vorbereitung
Übungs-Datensatz (wie beim letzten Termin)
Eine Stichprobe von 100 Schülerinnen und Schülern hat einen Lese- und einen Mathematiktest, sowie zusätzlich einen allgemeinen Intelligenztest, bearbeitet. Im Datensatz enthalten ist zudem das Geschlecht (female, 0=m, 1=w). Wir laden als erstes die nötigen R-Pakete sowie den Datensatz.
library(car) library(MASS) library(ggplot2) library(lm.beta) # erforderlich für standardiserte Regressionsgewichte # Datensatz laden data("Schulleistungen", package = "PsyBSc7")
Sagen Sie mit Hilfe eines Regressionsmodells die Leseleistung durch das Geschlecht und durch die Intelligenz vorher (im Datensatz Schulleistungen) und speichern sie diese Modell unter dem Namen mod ab. Geben Sie das Objekt mod aus.
# Modell zur Vorhersage der Leseleistung mod <- ... ...
# Verwenden Sie den Befehl lm(...~...)
# Finale Lösung: mod <- lm(reading ~ female + IQ, data = Schulleistungen) mod
# For further use in exercises mod <- lm(reading ~ female + IQ, data = Schulleistungen)
Im Output sehen wir die Parameterschätzungen unseres Regressionsmodells:
$$\text{Reading}_i=b_0+b_1\text{Female}_i+b_2\text{IQ}_i+\varepsilon_i,$$
für $i=1,\dots,100=:n$. Wir wollen uns die Ergebnisse unserer Regressionsanalyse noch detaillierter anschauen. Dazu können wir wieder die summary Funktion anwenden. Um auch die standardisierten Ergebnisse zu erhalten, verwenden wir die Funktion lm.beta (lm steht hier für lineares Modell und beta für die standardisierten Koeffizienten. Achtung: Häufig werden allerdings auch unstandardisierten Regressionkoeffizienten als $\beta$s bezeichnet, sodass darauf stehts zu achten ist - siehe bspw. Appendix A).
# Verwenden Sie den Befehl lm.beta(...) # um das geschätzte Modell zu standardisieren
# Verwenden Sie den Befehl summary(...) # um dem standardisierten Modell die Zusammenfassungsstatistiken zu entlocken
#Finale Lösung summary(lm.beta(mod))
# For further use mod <- lm(reading ~ female + IQ, data = Schulleistungen)
# for later use stdsum <- summary(lm.beta(mod))
An der standardisierten Lösung fällt auf, dass das Interzept $0$ ist. Dies macht völlig Sinn, da beim Standardisieren die Mittelwerte aller Variablen auf $0$ und die Standardabweichungen auf $1$ gesetzt werden. So auch beim Kriterium, der Leseleistung. Damit beträgt der Mittelwert des Kriteriums und der Mittelwert aller Prädiktore jeweils $0$. Das Interzept gibt an, welchen Wert das Kriterium hat, wenn alle Prädiktoren den Wert $0$ aufweisen. Somit muss das Interzept hier genau $0$ betragen. Den Koeffizienten entnehmen wir, dass die Prädiktoren female und IQ signifikante Anteile am Kriterium erklären (Mit einer Irrtumswahrscheinlichkeit von 5% ist der Regressionkoeffizent und damit die lineare Beziehung zwischen dem jeweiligen Prädiktor (z.B. IQ) und der Leseleistung in der Population nicht 0). Der standardisierte Koeffizient des IQ ist hierbei deutlich größer ($\hat{\beta}2=$ r round(stdsum$coefficients[3,2], 2)) als vom Geschlecht ($\hat{\beta}_1=$ r round(stdsum$coefficients[2,2], 2)). Damit ist ersichtlich, dass der IQ mehr Variation am Kriterium Leseleisstug erklärt. Achtung: Diese Information ist den unstandardisierten Koeffizienten nicht zu entnehmen, da ein direkter Vergleich der unstandardisierten Koeffizienten nicht möglich ist. Der unstandardisierte Effekt des Geschlechts auf die Leseleistung beträgt $\hat{b}_1=$ r round(stdsum$coefficients[2,1], 2) und der Effekt des IQ lediglich $\hat{b}_2=$ r round(stdsum$coefficients[3,1], 2). Dies liegt natürlich daran, dass der IQ eine wesentlich größere Streuung ($\sigma{IQ}^2$=r round(var(Schulleistungen$IQ), 2)) aufweist als die dummykodierte Variable female ($\sigma_{female}^2$=r round(var(Schulleistungen$female), 2)). Insgesamt werden r round(stdsum$r.squared*100, 2)% der Variation der Leseleistung durch die beiden Prädiktoren erklärt.
Gegenstand der folgenden Sitzung ist die Prüfung von Voraussetzungen der Regressionsanalyse:
- Linearität (Wiederholung)
- Homoskedastizität (Wiederholung)
- Normalverteilung der Residuen (Wiederholung)
- Multikollinearität
- Identifikation von Ausreißern und einflussreichen Datenpunkten
Im Folgenden werden wir mit den unstandardisierten Modell weiterarbeiten, welches wir im Objekt mod gespeichert haben.
Linearität
Eine grafische Prüfung der partiellen Linearität zwischen den einzelnen Prädiktorvariablen und dem Kriterium kann durch partielle Regressionsplots (engl. Partialplots) erfolgen. Dafür sagen wir in einem Zwischenschritt einen einzelnen Prädiktor durch alle anderen Prädiktoren im Modell vorher und speichern die Residuen, die sich aus dieser Regression ergeben. Diese kennzeichnen den eigenständigen Anteil, den ein Prädiktor nicht mit den anderen Prädiktoren gemein hat. Dann werden die Residuen aus dieser Vorhersage gegen die Residuen des Kriteriums bei Vorhersage durch den jeweiligen Prädiktor dargestellt. Diese Grafiken können Hinweise auf systematische nicht-lineare Zusammenhänge geben, die in der Modellspezifikation nicht berücksichtigt wurden. Die zugehörige R-Funktion des car Pakets heißt avPlots und braucht als Argument lediglich das spezifizierte Regressionsmodell mod.
avPlots(model = ..., pch = 16, lwd = 4)
# Finale Lösung: avPlots(model = mod, pch = 16, lwd = 4)
Mit Hilfe der Argumente pch=16 und lwd=4 werden die Darstellung der Punkte (ausgefüllt anstatt leer) sowie die Dicke der Linie (vierfache Dicke) manipuliert. Den Achsenbeschriftungen ist zu entnehmen, dass auf der Y-Achse jeweils reading | others dargestellt ist. Die vertikale Linie | steht hierbei für den mathematischen Ausdruck gegeben. Others steht hierbei für alle weiteren (anderen) Prädiktoren im Modell. Dies bedeutet, dass es sich hierbei um die Residuen aus der Regression von reading auf alle anderen Prädiktoren handelt. Bei den unabhängigen Variablen (UV, female, IQ) steht UV | others also jeweils für die jeweilige UV gegeben der anderen UVs im Modell. Somit beschreiben die beiden Plots jeweils die Beziehungen, die die UVs über die anderen UVs im Modell hinaus mit dem Kriterium (AV, abhängige Variable) haben.
Verteilung der Residuen
Homoskedastizität
Die Varianz der Residuen sollte unabhängig von den Ausprägungen der Prädiktoren sein. Dies wird i.d.R. grafisch geprüft, indem die Residuen $e_i$ gegen die vorhergesagten Werte $\hat{y}_i$ geplottet werden: der sogenannte Residuenplot (engl.: residual plot). In diesem Streudiagramm sollten die Residuen gleichmäßig über $\hat{y}_i$ streuen und keine systematischen Trends (linear, quadratisch, auffächernd, o.ä.) erkennbar sein. Prüfung nicht-linearer Zusammenhänge: Die Funktion residualPlots des Pakets car erzeugt separate Streudiagramme für die Residuen in Abhängigkeit von jedem einzelnen Prädiktor $x_j$ und von den vorhergesagten Werten $\hat{y}_i$ ("Fitted Values"); als Input braucht sie das Modell mod. Zusätzlich wird für jeden Plot ein quadratischer Trend eingezeichnet und auf Signifikanz getestet, wodurch eine zusätzliche Prüfung auf nicht-lineare Effekte erfolgt. Sind diese Test nicht signifikant, ist davon auszugehen, dass diese Effekte nicht vorliegen und die Voraussetzungen nicht verletzt sind.
# Grafisch (+ Test auf Nicht-Linearität) residualPlots(model = ..., pch = 16)
# Finale Lösung: residualPlots(model = mod, pch = 16)
Prüfung linearer Zusammenhänge: Die Funktion ncvTest (als Input braucht sie mod) prüft, ob die Varianz der Residuen bedeutsam linear (!) mit den vorhergesagten Werten zusammenhängt. Wird dieser Test signifikant, ist die Annahme der Homoskedaszidität verletzt (Wikipedia: Breusch-Pagan). Wird er nicht signifikant, kann dennoch eine Verletzung vorliegen, z.B. ein nicht-linearer Zusammenhang.
# Test For Non-Constant Error Variance ncvTest(model = ...)
# Finale Lösung: ncvTest(mod)
Normalverteilung
Voraussetzung für die Signifikanztests im Kontext der linearen Regression ist die Normalverteilung der Residuen. Auch diese Annahme wird i.d.R. grafisch geprüft. Hierfür bietet sich zum einen ein Histogramm der Residuen an, zum anderen ein Q-Q-Diagramm (oder "Quantile-Quantile-Plot"). Zusätzlich zur grafischen Darstellung erlaubt z.B. der Shapiro-Wilk-Test auf Normalität (shapiro.test) einen Signifkanztest der Nullhypothese, dass die Residuen normalverteilt sind. Auch der Kolmogorov-Smirnov (ks.test) Test, welcher eine deskriptive mit einer theoretischen Verteilung vergleicht, wäre denkbar.
Wir beginnen mit der Vorbereitung der Daten. Wir wollen die studentisierten Residuen grafisch mit dem R-Paket ggplot2darstellen. Der Befehl studres aus dem MASS Paket speichert die Residuen aus einem lm Objekt, also aus dem definierten Modell mod, und studentisiert diese. Das Studentisieren der Residuen bezeichnet eine Art der Standardisierung, sodass anschließend der Mittelwert 0 und die Varianz 1 ist (somit lassen sich solche Plots immer gleich interpretieren). Speichern Sie die Residuen aus mod als res ab und erzeugen Sie daraus einen data.frame, den Sie unter dem Namen df_res speichern. Diesen wollen wir später in den Grafiken mit ggplot2 verwenden.
# Verwenden Sie die Befehle studres(...)
# Verwenden Sie die Befehle res <- studres(...)
# Verwenden Sie die Befehle res <- studres(...) # und df_res <- data.frame(...)
# Verwenden Sie die Befehle res <- studres(...) # und df_res <- data.frame(...) # geben Sie das Objekt des Data.frames zurück (sehen Sie sich dieses an) df_res
# Intro into excercise mod <- lm(reading ~ female + IQ, data = Schulleistungen) # Begin grading: res <- studres(mod) # Studentisierte Residuen als Objekt speichern df_res <- data.frame(res) # als Data.Frame für ggplot res <- df_res grade_result( fail_if(~ (!identical(.result, res)), 'Haben Sie die richtigen Input-Argumente verwendet? Sie wollen die studentisierten Resiuden abspeichern - als data.frame. Geben Sie diese Objekt aus (indem Sie den Namen des Objekts ohne Zuordnung als letztes schreiben).'), pass_if(~ identical(.result, res)), correct = 'Sehr gut! Sie haben erfolgreich die Residuen als data.frame abgespeichert.', incorrect = 'Leider falsch.', glue_correct = '{.correct}', glue_incorrect = '{.incorrect} {.message}')
res <- studres(mod) # Studentisierte Residuen als Objekt speichern df_res <- data.frame(res) # als Data.Frame für ggplot
# For further use in exercises mod <- lm(reading ~ female + IQ, data = Schulleistungen) res <- studres(mod) # Studentisierte Residuen als Objekt speichern df_res <- data.frame(res) # als Data.Frame für ggplot
Im folgenden Block sehen wir den Code für ein Histogramm in ggplot2-Notation. Hier sind einige Zusatzeinstellungen gewählt, die das Histogramm optisch aufbereiten. Vervollständigen Sie den Code mit Hilfe der Tipps.
# Grafisch: Histogramm mit Normalverteilungskurve ggplot(data = ???, aes(x = ??)) + geom_histogram(aes(y =..density..), bins = 15, # Wie viele Balken sollen gezeichnet werden? colour = "blue", # Welche Farbe sollen die Linien der Balken haben? fill = "skyblue") + # Wie sollen die Balken gefüllt sein? stat_function(fun = dnorm, args = list(mean = mean(???), sd = sd(???)), col = "darkblue") + # Füge die Normalverteilungsdiche hinzu und nutze den empirischen Mittelwert und die empirische Standardabweichung, wähle dunkelblau als Linienfarbe labs(title = "Histogramm der Residuen mit Normalverteilungsdichte", x = "Residuen") # Füge eigenen Titel und Achsenbeschriftung hinzu
# an data # wird in ggplot immer der data.frame mit den Daten übergeben
# der Dataframe enthält in unserem Fall eine Variable, welche die Residuen enthält. Den Namen haben wir selbst vergeben. Diesen können Sit mit colnames() # in Erfahrung bringen (welcher die Namen der Spalten ausgibt)
# Auf die x-Achse in einem Histogramm wird immer die Variable selbst abgetragen. Diese müssen sie an aes(x = ...) # übergeben
# mit stat_function(fun = dnorm, args = list(...)) # übergeben wir eine statistische Funktion, # welche dazu geplottet werden soll. # Hier ist dies die Dichte der Normalverteilung "dnorm". Diese Funktion wird normalerweise mit dnorm(x = ..., mean = ..., sd = ...) # angesprochen. Wir können mit dieser die Dichte der Normalverteilung # am Punkt x berechnen. Um welche Normalverteilung es sich handelt, geben wir mit "mean" und "sd" an. # Die Normalverteilung ist eindeutig durch Mittelwert und Standardabweichung definiert.
# mit stat_function(fun = dnorm, args = list(...)) # übergeben wir eine statistische Funktion, # welche dazu geplottet werden soll. # Hier ist dies die Dichte der Normalverteilung "dnorm". Diese Funktion wird normalerweise mit dnorm(x = ..., mean = ..., sd = ...) # angesprochen. Wir können mit dieser die Dichte der Normalverteilung # am Punkt x berechnen. Um welche Normalverteilung es sich handelt, geben wir mit "mean" und "sd" an. # Die Normalverteilung ist eindeutig durch Mittelwert und Standardabweichung definiert. # Die Argumente, die dnorm normalerweise erhält, können wir in args als Liste in stat_function abgeben # Hier brauchen wir allerdings "x" nicht, da das übrigbleibenden (wenn wir mean und sd angeben) Argument # für den Plot manipuliert wird (von ggplot2). Somit erhalten wir für die Ausprägungen entlang der # x-Achse die Dichte der Normalverteilung (und zwar der Residuen!).
# Finale Lösung: ggplot(data = df_res, aes(x = res)) + geom_histogram(aes(y =..density..), bins = 15, # Wie viele Balken sollen gezeichnet werden? colour = "blue", # Welche Farbe sollen die Linien der Balken haben? fill = "skyblue") + # Wie sollen die Balken gefüllt sein? stat_function(fun = dnorm, args = list(mean = mean(res), sd = sd(res)), col = "darkblue") + # Füge die Normalverteilungsdiche "dnorm" hinzu und nutze den empirischen Mittelwert und die empirische Standardabweichung "args = list(mean = mean(res), sd = sd(res))", wähle dunkelblau als Linienfarbe labs(title = "Histogramm der Residuen mit Normalverteilungsdichte", x = "Residuen") # Füge eigenen Titel und Achsenbeschriftung hinzu
Nutzen wir nur die Defaulteinstellung des Histogramms (bis auf bins = 15 - für die Vergleichbarkeit der beiden Grafiken), sieht es so aus:
# hier nochmal nur das Nötigste: ggplot(data = df_res, aes(x = res)) + geom_histogram(aes(y =..density..), bins = 15) + stat_function(fun = dnorm, args = list(mean = mean(res), sd = sd(res)))
Ähnliche Informationen sind aus dem Q-Q-Plot zu entnehmen. Dafür kann der Befehl qqplot aus dem car Paket genutzt werden. Hier kann die Verteilung, gegen die geprüft werden soll, direkt mit distribution festgelegt werden kann.
# Grafisch: Q-Q-Diagramm mit der car Funktion qqPlot qqPlot(..., pch = 16, distribution = ...)
qqplot() # brauch als Argument das Regressionmodell "mod".
# Die Normalverteilung heißt in R "norm".
# Finale Lösung: qqPlot(mod, pch = 16, distribution = "norm")
Nun wollen wir die Voraussetzung zusätzlich mittels statistischer Tests überprüfen. Diese testen im Grunde genommen die Null-Hypothese: $H_0:$ "Normalverteilung liegt vor". (Wir wissen allerdings, dass wir Hypothesen nur verwerfen können und nicht wirklich annehmen, weshalb diese Tests immer nur unter Vorbehalt interpretiert werden sollten - Aussagen wie etwa: "$H_0$ nicht verworfen; daraus folgt, dass in der Population Normalverteilung herrscht", sind nicht zulässig. Wird die $H_0$ nicht verworfen, nehmen wir nur weiterhin an, dass Normalverteilung vorliegt, ein Beweis ist dies jedoch nicht).
Die Funktion shapiro.test und ks.test benötigen die Variable, die auf Normalverteilung geprüft werden soll, als Argument; hier also res.
# Test auf Abwweichung von der Normalverteilung mit dem Shapiro Test shapiro.test(...)
# Finale Lösung: shapiro.test(res)
# Test auf Abwweichung von der Normalverteilung mit dem Kolmogorov-Smirnov Test ks.test(..., "...", mean(...), sd(...))
# Die Argumente in ks.test # sind wie folgt angeordnet: # (Variable, "Verteilung (Wahrscheinlichkeitswerte)", Argument1_der_Verteilung, Argument2_der_Verteilung)
# Variable = res # "Verteilung (Wahrscheinlichkeitswerte)" = "pnorm", die Verteilungswerte der Normalverteilung mit den Argumenten q=Quantil, sowie Mittelwert und Standardabweichung # Argument1_der_Verteilung = Mittelwert von res # Argument2_der_Verteilung = Standardabweichung von res
# Finale Lösung: ks.test(res, "pnorm", mean(res), sd(res))
Für unser Modell zeigen alle Ergebnisse übereinstimmend, dass die Residuen normalverteilt sind (bzw. lehnen diese Hypothese nicht ab: Shapiro-Wilk Test und Kolmogorov-Smirnov Tests gegen die Normalverteilung). Durch die Angabe von distribution = "norm als Argument für die Funktion qqPlot des car-Pakets erzeugt den Vergleich der vorliegenden Variable gegen eine Normalverteilung (auch der Vergleich zu anderen Verteilungen wäre möglich). Dem Kolmogorov-Smirnov Test gegen die Normalverteilung ks.test muss zunächst als erstes Argument res, die Variable, deren Verteilung untersucht werden soll, übergeben werden. "pnorm" gibt an, dass die empirische Verteilung unserer Daten res mit der Normalverteilung verglichen werden soll (wird hier eine andere Verteilung gewählt, so kann gegen diese geprüft werden; auch können zwei empirisiche Verteilungen auf Gleichheit untersucht werden). Anschließend geben wir die verteilungsspezifischen Parameter; hier mean(res) (der Mittelwert) sowie sd(res) (die Standardabweichung) an, da die Normalverteilung genau durch diese beiden Parameter bestimmt ist (hier könnten auch theoretisch angenommene Parameter gewählt werden, z.B. 0 und 1 für eine Mittelwert von 0 und einer Varianz von 1 und es kann geprüft werden, ob die beobachteten Daten dieser Verteilung widersprechen).
Wäre der Shapiro-Wilk-Normality Tests oder der Kolmogorov-Smirnov Tests gegen die Normalverteilung signifikant auf dem 5%-Niveau, würde dies dafür sprechen, dass "mit einer Irrtumswahrscheinlichkeit von 5% die Null-Hypothese auf Gleichheit der Normalverteilung mit unserer empirischen Verteilung verworfen werden muss, da in der Population die Abweichung der empirischen zur theoretischen Normalverteilung mit der Irrtumswahrscheinlichkeit von 5% ungleich 0 ist".
Schiefe Verteilungen, welche sich im Histogramm oder im Q-Q-Plot wiederspiegeln, sowie zu signifikanten Shapiro-Wilk-Normality Tests oder Kolmogorov-Smirnov Tests gegen die Normalverteilung führen können, könnten Indizien für nichtlineare (z.B. quadratische) Beziehungen im Datensatz sein. Weitere Analysen wären hierfür nötig (siehe hierzu PsyBSc7::Sitzung_8()).
Multikollinearität
Multikollinearität ist ein potenzielles Problem der multiplen Regressionsanalyse und liegt vor, wenn zwei oder mehrere Prädiktoren hoch miteinander korrelieren. Hohe Multikollinearität
- schränkt die mögliche multiple Korrelation ein, da die Prädiktoren redundant sind und überlappende Varianzanteile in $y$ erklären.
- erschwert die Identifikation von bedeutsamen Prädiktoren, da deren Effekte untereinander konfundiert sind (die Effekte können schwer voneinander getrennt werden).
- bewirkt eine Erhöung der Standardfehler der Regressionkoeffizienten (der Standardfehler ist die Standardabweichung zu der Varianz der Regressionskoeffizienten bei wiederholter Stichprobenziehung und Schätzung). Dies bedeutet, dass die Schätzungen der Regressionsparameter instabil, und damit weniger verlässlich, werden. Für weiterführende Informationen zur Instabilität und Standardfehlern siehe Appendix A.
Multikollinearität kann durch Inspektion der bivariaten Zusammenhänge (Korrelationsmatrix) der Prädiktoren $x_j$ untersucht werden, dies kann aber nicht alle Formen von Multikollinearität aufdecken. Darüber hinaus ist die Berechung der sogennanten Toleranz und des Varianzinflationsfaktors (VIF) für jeden Prädiktor möglich. Hierfür wird für jeden Prädiktor $x_j$ der Varianzanteil $R_j^2$ berechnet, der durch Vorhersage von $x_j$ durch alle anderen Prädiktoren in der Regression erklärt wird. Toleranz und VIF sind wie folgt definiert:
- $T_j = 1-R^2_j = \frac{1}{VIF_j}$
- $VIF = \frac{1}{1-R^2_j} = \frac{1}{T_j}$
Offensichtlich genügt eine der beiden Statistiken, da sie vollständig ineinander überführbar und damit redundant sind. Empfehlungen als Grenzwert für Kollinearitätsprobleme sind z. B. $VIF_j>10$ ($T_j<0.1$) (siehe Eid, Gollwitzer, & Schmitt, 2017, S. 712 und folgend). Die Varianzinflationsfaktoren der Prädiktoren im Modell können mit der Funktion vif des car-Paktes bestimmt werden, der Toleranzwert als Kehrwert des VIFs.
# Korrelation der Prädiktoren cor(Schulleistungen$female, Schulleistungen$IQ)
Die Prädiktoren sind nur schwach negativ korreliert. Wir schauen uns trotzdem den VIF und die Toleranz an. Dazu übergeben wir wieder das definierte Regressionsmodell an vif.
# Varianzinflationsfaktoren vif(...)
# Finale Lösung: vif(mod)
Berechnen Sie die Toleranz.
# Verwenden Sie den Befehl vif(...) # von zuvor und berechnen Sie den Kehrwert
# Intro into excercise mod <- lm(reading ~ female + IQ, data = Schulleistungen) # Begin grading: res <- 1/vif(mod) grade_result( fail_if(~ (!identical(.result, res)), 'Haben Sie die richtigen Input-Argumente verwendet? Sie wollen die Toleranz als Kehrwert des VIF bestimmen.'), pass_if(~ identical(.result, res)), correct = 'Sehr gut! Sie haben erfolgreich die Toleranz aus dem VIF berechnet.', incorrect = 'Leider falsch.', glue_correct = '{.correct}', glue_incorrect = '{.incorrect} {.message}')
Für unser Modell wird ersichtlich, dass die Prädiktoren praktisch unkorreliert sind und dementsprechend kein Multikollinearitätsproblem vorliegt. Unabhängigkeit folgt hieraus allerdings nicht, da nichtlineare Beziehungen zwischen den Variablen bestehen könnten, die durch diese Indizes nicht abgebildet werden.
Identifikation von Ausreißern und einflussreichen Datenpunkten
Die Plausibilität unserer Daten ist enorm wichtig. Aus diesem Grund sollten Aureißer oder einflussreiche Datenpunkte analysiert werden. Diese können bspw. durch Eingabefehler entstehen (Alter von 211 Jahren anstatt 21) oder es sind seltene Fälle (hochintelligentes Kind in einer Normstichprobe), welche so in natürlicher Weise (aber mit sehr geringer Häufigkeit) auftreten können. Es muss dann entschieden werden, ob Ausreißer die Repräsentativität der Stichprobe gefährden und ob diese besser ausgeschlossen werden sollten.
Hebelwerte $h_j$ erlauben die Identifikation von Ausreißern aus der gemeinsamen Verteilung der unabhängigen Variablen, d.h. einzelne Fälle, die weit entfernt vom Mittelwert der gemeinsamen Verteilung der unabhängigen Variablen liegen und somit einen starken Einfluss auf die Regressionsgewichte haben können. Diese werden mit der Funktion hatvalues ermittelt. Kriterien zur Beurteilung der Hebelwerte variieren, so werden von Eid et al. (2017, S. 707 und folgend) Grenzen von $2\cdot k / n$ für große und $3\cdot k / n$ für kleine Stichproben vorgeschlagen, in den Vorlesungsfolien werden Werte von $4/n$ als auffällig eingestuft (hierbei ist $k$ die Anzahl an Prädiktoren und $n$ die Anzahl der Beobachtungen). Alternativ zu einem festen Cut-Off-Kriterium kann die Verteilung der Hebelwerte inspiziert und diejenigen Werte kritisch betrachtet werden, die aus der Verteilung ausreißen. Die Funktion hatvalues erzeugt die Hebelwerte aus einem Regression-Objekt. Wir wollen diese als Histogramm darstellen.
# Hebelwerte n <- length(residuals(...)) h <- hatvalues(...) # Hebelwerte df_h <- data.frame(h) # als Data.Frame für ggplot ggplot(data = df_h, aes(x = h)) + geom_histogram(aes(y =..density..), bins = 15)+ geom_vline(xintercept = 4/n, col = "red") # Cut-off bei 4/n
# Finale Lösung n <- length(residuals(mod)) h <- hatvalues(mod) # Hebelwerte df_h <- data.frame(h) # als Data.Frame für ggplot ggplot(data = df_h, aes(x = h)) + geom_histogram(aes(y =..density..), bins = 15)+ geom_vline(xintercept = 4/n, col = "red") # Cut-off bei 4/n
Cook's Distanz $CD_i$ gibt eine Schätzung, wie stark sich die Regressionsgewichte verändern, wenn eine Person $i$ aus dem Datensatz entfernt wird. Fälle, derem Elimination zu einer deutlichen Veränderung der Ergebnisse führen würden, sollten kritisch geprüft werden. Als einfache Daumenregel gilt, dass $CD_i>1$ auf einen einflussreichen Datenpunkt hinweist. Cook's Distanz kann mit der Funktion cooks.distance ermittelt werden.
# Cooks Distanz CD <- cooks.distance(...) # Cooks Distanz df_CD <- data.frame(CD) # als Data.Frame für ggplot ggplot(data = df_CD, aes(x = CD)) + geom_histogram(aes(y =..density..), bins = 15)+ geom_vline(xintercept = 1, col = "red") # Cut-Off bei 1
# Cooks Distanz CD <- cooks.distance(mod) # Cooks Distanz df_CD <- data.frame(CD) # als Data.Frame für ggplot ggplot(data = df_CD, aes(x = CD)) + geom_histogram(aes(y =..density..), bins = 15)+ geom_vline(xintercept = 1, col = "red") # Cut-Off bei 1
Die Funktion influencePlot des car-Paktes erzeugt ein "Blasendiagramm" zur simultanen grafischen Darstellung von Hebelwerten (auf der x-Achse), studentisierten Residuen (auf der y-Achse) und Cooks Distanz (als Größe der Blasen). Vertikale Bezugslinien markieren das Doppelte und Dreifache des durchschnittlichen Hebelwertes, horizontale Bezugslinien die Werte -2, 0 und 2 auf der Skala der studentisierten Residuen. Fälle, die nach einem der drei Kriterien als Ausreißer identifiziert werden, werden im Streudiagramm durch ihre Zeilennummer gekennzeichnet. Diese Zeilennummern können verwendet werden, um sich die Daten der auffälligen Fälle anzeigen zu lassen. Sie werden durch InfPlot ausgegeben werden. Auf diese kann durch as.numeric(row.names(InfPlot)) zugegriffen werden.
# Blasendiagramm mit Hebelwerten, studentisierten Residuen und Cooks Distanz # In "IDs" werden die Zeilennummern der auffälligen Fälle gespeichert, # welche gleichzeitig als Zahlen im Blasendiagramm ausgegeben werden InfPlot <- influencePlot(mod) IDs <- as.numeric(row.names(InfPlot)) # Werte der identifizierten Fälle InfPlot
InfPlot <- influencePlot(mod) IDs <- as.numeric(row.names(InfPlot))
Schauen wir uns die möglichen Ausreißer an und standardisieren die Ergebnisse für eine bessere Interpretierbarkeit.
# Rohdaten der auffälligen Fälle (gerundet für bessere Übersichtlichkeit) round(Schulleistungen[IDs,],2) # z-Standardisierte Werte der auffälligen Fälle round(scale(Schulleistungen)[IDs,],2)
Die Funktion scale z-standardisiert den Datensatz, mit Hilfe von [IDs,], werden die entsprechenden Zeilen der Ausreißer aus dem Datensatz ausgegeben und anschließend auf 2 Nachkommastellen gerundet. Mit Hilfe der z-standardisieren Ergebnisse lassen sich Ausreißer hinsichtlich ihrer Ausprägungen einordnen:
Interpretation
Was ist an den fünf identifizierten Fällen konkret auffällig?
- Fall 6: recht niedrige Lesekompetenz bei gleichzeitig durchschnittlichem bis überdurchschnittlichem IQ
- Fall 9: Sehr hohe Werte in IQ, Lesen & Mathe
- Fall 33: Unterdurchschnittliche Leseleistung "trotz" eher durchschnittlicher Intelligenz
- Fall 80: Sehr niedriger IQ, gleichzeitig durchschnittliche bis überdurchschnittliche Lesekompetenz
- Fall 99: Sehr niedrige Werte in IQ, Lesekompetenz und Mathematik
Die Entscheidung, ob Ausreißer oder auffällige Datenpunkte aus Analysen ausgeschlossen werden, ist schwierig und kann nicht pauschal beantwortet werden. Im vorliegenden Fall wäre z.B. zu überlegen, ob die Intelligenztestwerte der Fälle 80 und 99, die im Bereich von Lernbehinderung oder sogar geistiger Behinderung liegen, in einer Stichprobe von Schülern/innen aus Regelschulen als glaubwürdige Messungen interpretiert werden können oder als Hinweise auf mangelndes Commitment bei der Beantwortung.
Aufgaben
Die Aufgaben zu dieser Sitzung finden Sie im Appendix B.
Appendix A {#AppendixA}
Multikollinearität und Standardfehler
Im Folgenden stehen $\beta$s für unstandardisierte Regressionskoeffizienten.
Für eine einfache Regressionsgleichung mit $$Y_i=\beta_0 + \beta_1x_{i1} + \beta_2x_{i2} + \varepsilon_i$$ kann die selbe Gleichung auch in Matrixnotation formuliert werden $$Y = X\beta + \varepsilon.$$ Hier ist $X$ die Systemmatrix, welche die Zeilenvektoren $X_i=(1, x_{i1}, x_{i2})$ enthält. Die Standardfehler, welche die Streuung der Parameter $\beta:=(\beta_0,\beta_1,\beta_2)$ beschreiben, lassen sich wie folgt ermitteln. Wir bestimmen zunächst die Matrix $I$ wie folgt $$I:=(X'X)^{-1}\hat{\sigma}^2_e,$$ wobei $\hat{\sigma}^2_e$ die Residualvarianz unserer Regressionanalyse ist (also der nicht-erklärte Anteil an der Varianz von $Y$). Aus der Matrix $I$ erhalten wir die Standardfehler sehr einfach: Sie stehen im Quadrat auf der Diagonale. Das heißt, die Standardfehler sind $SE(\beta)=\sqrt{\text{diag}(I)}$ (Wir nehmen mit $\text{diag}$ die Diagonalelemente aus $I$ und ziehen aus diesen jeweils die Wurzel: der erste Eintrag ist der $SE$ von $\beta_0$; also $SE(\beta_0)=\sqrt{I_{11}}$; der zweite von $\beta_1$;$SE(\beta_1)=\sqrt{I_{22}}$; usw.). Was hat das nun mit der Kollinearität zu tun? Wir wissen, dass in $X'X$ die Information über die Kovariation im Datensatz steckt (dafür muss nur noch durch die Stichprobengröße geteilt werden und das Vektorprodukt der Mittelwerte abgezogen werden; damit wir eine Zentrierung um die Mittelwert sowie eine Normierung an der Stichprobengröße vorgenommen). Beispielsweise lässt sich die empirische Kovarianzmatrix $S$ zweier Variablen $z_1$ und $z_2$ sehr einfach bestimmen mit $Z:=(z_1, z_2)$: $$ S=\frac{1}{n}Z'Z - \begin{pmatrix}\overline{z}_1\\overline{z}_2 \end{pmatrix}\begin{pmatrix}\overline{z}_1&\overline{z}_2 \end{pmatrix}.$$ Weitere Informationen hierzu (Kovarianzmatrix und Standardfehler) sind im Appendix B (sowie auch in einigen Kapiteln von Eid et al. (2017) Unterpunkt 5.2-5.4 bzw. ab Seite 1058) nachzulesen.
Insgesamt bedeutet dies, dass die Standardfehler von der Inversen der Kovarianzmatrix unserer Daten sowie von der Residualvarianz abhängen. Sie sind also groß, wenn die Residualvarianz groß ist (damit ist die Vorhersage von $Y$ schlecht) oder wenn die Inverse der Kovarianzmatrix groß ist (also wenn die Variablen stark redundant sind und somit hoch mit einander korrelieren). Nehmen wir dazu der Einfachheit halber an, dass $\hat{\sigma}_e^2=1$ (es geht hier nur um eine numerische Präsentation der Effekte, nicht um ein sinnvolles Modell) sowie $n = 100$ (Stichprobengröße). Zusätzlich gehen wir von zentrierten Variablen (Mittelwert von 0) aus. Dann lässt sich aus $X'X$ durch Division durch $100$ die Kovarianzmatrix der Variablen bestimmen. Wir gucken uns drei Fälle an mit
Fall 1: $X'X=\begin{pmatrix}100&0&0\0&100&0\0&0&100\end{pmatrix}.\qquad$ Fall 2: $X'X=\begin{pmatrix}100&0&0\0&100&99\0&99&100\end{pmatrix} \qquad$ und Fall 3: $X'X=\begin{pmatrix}100&0&0\0&100&100\0&100&100\end{pmatrix}$.
Hierbei ist zu beachten, dass $X$ die Systemmatrix ist, welche auch die $1$ des Interzepts enthält. Natürlich ist eine Variable von einer Konstanten unabhängig, weswegen die erste Zeile und Spalte von $X'X$ jeweils der Vektor $(100, 0, 0)$ ist. Die zugehörigen Korrelationsmatrizen können durch Divison durch 100 berechnet werden (da wir zentrierte Variablen haben, die Stichprobengröße gleich 100 ist und die Varianzen der Variablen gerade 100 sind!). Wir betrachten nur die Minormatrizen, aus welchen die 1. Zeile und die 1. Spalte gestrichen wurden. Diese teilen wir durch 100 und erhalten die Korrelationsmatrix der Variablen:
Fall 1: $\Sigma_1=\begin{pmatrix}1&0\0&1\end{pmatrix}.\qquad$ Fall 2: $\Sigma_2=\begin{pmatrix}1&.99\.99&1\end{pmatrix} \qquad$ und Fall 3: $\Sigma_3=\begin{pmatrix}1&1\1&1\end{pmatrix}$.
Im Fall 1 sind die zwei Variablen unkorreliert. Die Inverse ist leicht zu bilden.
XX_1 <- matrix(c(100,0,0, 0,100,0, 0,0,100),3,3) XX_1 # Die Matrix X'X im Fall 1 I_1 <- solve(XX_1)*1 # I (*1 wegen Residualvarianz = 1) I_1 sqrt(diag(I_1)) # Wurzel aus den Diagonalelementen der Inverse = SE, wenn sigma_e^2=1
Die Standardfehler sind nicht sehr groß: alle liegen bei $0.1$.
Im Fall 2 sind die zwei Variablen fast perfekt (zu $.99$) korreliert - es liegt hohe Multikollinearität vor. Die Inverse ist noch zu bilden. Die Standardfehler sind deutlich erhöht im Vergleich zu Fall 1.
XX_2 <- matrix(c(100,0,0, 0,100,99, 0,99,100),3,3) XX_2 # Die Matrix X'X im Fall 2 I_2 <- solve(XX_2)*1 # I (*1 wegen Residualvarianz = 1) I_2 sqrt(diag(I_2)) # SEs im Fall 2 sqrt(diag(I_1)) # SEs im Fall 1
Die Standardfehler des Fall 2 sind sehr groß im Vergleich zu Fall 1 (mehr als sieben Mal so groß); nur der Standardfehler des Interzept bleibt gleich. Die Determinante von $X'X$ in Fall 2 liegt deutlich näher an $0$ im Vergleich zu Fall 1; hier: $10^6$.
det(XX_2) # Determinante Fall 2 det(XX_1) # Determinante Fall 1
Im Fall 3 sind die zwei Variablen perfekt korreliert - es liegt perfekte Multikollinearität vor. Die Inverse kann nicht gebildet werden (da $\text{det}(X'X) = 0$). Die Standardfehler können nicht berechnet werden. Eine Fehlermeldung wird ausgegeben.
XX_3 <- matrix(c(100,0,0, 0,100,100, 0,100,100),3,3) XX_3 # Die Matrix X'X im Fall 3 det(XX_3) # Determinante on X'X im Fall 3
I_3 <- solve(XX_3)*1 # I (*1 wegen Residualvarianz = 1) I_3 sqrt(diag(I_3)) # Wurzel aus den Diagonalelementen der Inverse = SE, wenn sigma_e^2=1 # hier wird eine Fehlermeldung ausgegeben, wodurch der Code nicht ausführbar ist und I_3 nicht gebildet werden kann: # Error in solve.default(XX_3) : # Lapack routine dgesv: system is exactly singular: U[2,2] = 0
Der VIF bzw. die Toleranz quantifizieren die Korrelation zwischen den beiden Variablen. Der VIF wäre in diesen Analysen im 1. Fall für beide Variablen 1, im 2. Fall für beide Variabeln 50.25 und im 3. Fall nicht berechenbar (bzw. $\infty$). Entsprechend wäre die Toleranz im 1. Fall 1 und 1, im 2. Fall 0.02 und 0.02 sowie im 3. Fall 0 und 0.
Dieser Exkurs zeigt, wie sich die Multikolinearität auf die Standardfehler und damit auf die Präzision der Parameterschätzung auswirkt. Inhaltlich bedeutet dies, dass die Prädiktoren redundant sind und nicht klar ist, welchem Prädiktor die Effekte zugeschrieben werden können.
Die Matrix $I$ ist im Zusammenhang mit der Maximum-Likelihood-Schätzung die Inverse der Fischer-Information und enthält die Informationen der Kovariationen der Parameterschätzer (diese Informationen enthält sie hier im Übrigen auch!).
Appendix B {#AppendixB}
Aufgaben
Datensatz (wie beim letzten Termin)
Der Beispieldatensatz enthält Daten zu Lesekompetenz aus der deutschen Stichprobe der PISA-Erhebung in Deutschland 2009.
library(car) library(MASS) library(lm.beta) # erforderlich für standardiserte Gewichte # Datensatz laden data("PISA2009", package = "PsyBSc7")
Im Einzelnen enthält der Datensatz die folgenden Variablen:
- Klassenstufe (
Grade) - Alter in Jahren (
Age) - Geschlecht (
Female, 0=m / 1=w) - Lesekompetenz (
Reading) - Lesefreude (
JoyRead) - Lernzeit in Minuten für Deutsch (
LearnMins) - Sozialstatus (
HISEI, "highest international socio-economic index of occupational status") - Fragebogen-Score für kulturelle Besitztümer zu Hause (
CultPoss, z. B. klassische Literatur, Kunstwerke) - Anzahl Bücher zu Hause (
Books, eigentlich eine ordinale Ratingskala) - Anzahl Fernseher zu Hause (
TVs) - Anzahl Computer zu Hause (
Computers) - Anzahl Autos zu Hause (
Cars) - Migrationshintergrund (
MigHintergrund, 0=beide Eltern in D geboren, 1=min. 1 Elternteil im Ausland geboren) - Bildungsabschluss des Vaters (
FatherEdu, International Standard Classification of Education) - Bildungsabschluss der Mutter (
MotherEdu, International Standard Classification of Education)
Vorhersage von Lesekompetenz mit individuellen Merkmalen der Schüler/innen
In der letzten Sitzung hatten wir herausgefunden, dass der Sozialstatus (HISEI), der Bildungsabschluss der Mutter (MotherEdu) und die Zahl der Bücher zu Hause (Books) bedeutsame Prädiktoren für die Lesekompetenz der Schüler/innen sind. Berechnen Sie dieses Modell erneut.
# Berechnung des Modells und Ausgabe der Ergebnisse m1 <- lm(Reading ~ HISEI + MotherEdu + Books, data = PISA2009) summary(lm.beta(m1))
Aufgabe 1. Prüfen Sie die Annahme der Homoskedastizität
a) durch Inspektion der Residuenplots und
b) den Breusch–Pagan Test auf nicht-konstante Varianz.
Aufgabe 2. Prüfen Sie die Annahme der Linearität
a) durch Inspektion der partiellen Regressionsplots und
b) durch den Test für quadratische Trends, der mit der Funktion residualPlots erzeugt wurde.
Aufgabe 3. Püfen Sie die Normalverteilung der Residuen
a) grafisch (Histogramm, Q-Q-Plot) und
b) über den Shapiro-Wilks-Test sowie den Kolmogorov-Smirnov-Test.
Aufgabe 4. Prüfen Sie die Multikollinearität der Prädiktoren im Modell
a) durch Inspektion der bivariaten Korrelationen und
b) durch Berechnung von VIF und Toleranz.
Aufgabe 5. Ausreißer und einflussreiche Datenpunkte
a) Inspizieren Sie die Verteilung der Hebelwerte $h_i$.
b) Inspizieren Sie die Verteilung von Cooks Distanz $CD_i$.
c) Sichten Sie die Rohdaten und standardisierten Werte auffälliger Fälle (verwenden Sie die Funktion influencePlot).
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.
$$ Diese Sitzung basiert zum Teil auf Unterlagen von Prof. Johannes Hartig aus dem SoSe 2019.*
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.