Statistik II
Multiple Regression
AG Forschungsmethoden und Kognitive Psychologie, Institut für Psychologie, Universität Bremen
Willkommen bei Stats by Randolph. Hier geht es zur Hauptseite mit weiteren Informationen und Inhalten.
Autor:innen dieser Seite: An den Inhalten dieser Seite haben mitgearbeitet: Markus Janczyk und Valentin Koob. Der Inhalt dieser Seite wird in der Lehre in den Studiengängen Psychologie von der AG Forschungsmethoden und Kognitive Psychologie an der Universität Bremen verwendet, steht aber allen Interessierten zur Verfügung. Rückmeldungen/Fehler/Vorschläge können gesendet werden an randolph@uni-bremen.
Versionshistory:
- v1.0: erste online-gestellte Version (03.3.2024)
7 Multiple Regression: Einführung
In Statistik I haben wir uns bereits mit der einfachen linearen Regression beschäftigt (Teil 6 in Statistik I). “Einfach” meint in diesem Kontext, dass wir nur eine Variable als Prädiktor berücksichtigt haben. Nach einer kurzen Wiederholung stellen wir dann einige Ergänzungen vor, die erst jetzt Sinn ergeben, nachdem wir andere Konzepte und Verfahren bereits kennen.
In der Regel sind messbare Dinge aber nicht nur von einem Umstand abhängig, sondern von mehreren. Um mehrere Einflussgrößen und ihre (kombinierte) Wirkung auf eine Variable zu erfassen, können wir die einfache lineare Regression (mit einem Prädiktor) so erweitern, dass mehrere Prädiktoren im Modell berücksichtigt werden. Wir reden dann von einer multiplen Regression. Weitgehend ist dies eine Verallgemeinerung, aber ein paar Besonderheiten z.B. bezüglich der Interpretation der Koeffizienten kommen hierbei ins Spiel. Wir werden das Grundprinzip zunächst am besonderen Beispiel zweier nicht-korrelierter Prädiktoren einführen und dann im Folgenden den allgemeinen Fall mit \(q\)-vielen Prädiktoren betrachten.
# Pakete die in diesem Teil benutzt werden:
library(car)
library(effects)
library(scatterplot3d)
libary(schoRsch)7.1 Einfache lineare Regression: Kurze Wiederholung und Ergänzungen
7.1.1 Kurze Wiederholung der einfachen, linearen Regression
Die Grundidee der einfachen, linearen Regression ist das Bestimmen einer optimalen Geraden zur Beschreibung empirischer Daten. Die schwarzen Punkte der folgenden Abbildung stellen dabei empirische Werte \((x_i,y_i)\) dar und die roten Punkte sind diejenigen Werte, die die eingezeichnete Gerade für die jeweiligen \(x_i\)-Werte vorhergesagen würde, also \((x_i,\hat{y}_i)\). Als Residuum \(e_i\) wird die vertikale Distanz der empirischen und vorhergesagten Werte bezeichnet, also \(e_i=y_i-\hat{y}_i\).
Die Gerade wird durch \[\hat{Y} = b\cdot X +a \] beschrieben, wobei mit \(\hat{Y}\) die durch die Gleichung vorhergesagten Werte auf dem Kriterium und mit \(X\) die Werte des Prädiktors \(X\) gemeint sind. \(b\) ist dann die Steigung der Geraden und \(a\) ihr Achsenabschnitt. Gesucht werden dann Werte für \(b\) und \(a\), mit denen die Summe der quadrierten Residuen (also der quadrierten Abweichungen der vorhergesagten von den empirischen Werten) minimal wird (“Kriterium der kleinsten Quadrate”), also das Minimum der Funktion \[Q=\sum_{i=1}^n e_i^2=\sum_{i=1}^n (y_i-\hat{y}_i)^2\ .\] Durch Bildung und Nullsetzung der partiellen Ableitungen von \(Q(b,a)\) ergeben sich daraus als optimale Werte für \(b\) und \(a\) (die vollständige Herleitung findet sich hier) \[\begin{equation*} b=\frac{\text{Kov}(X,Y)}{S_X^2}\text{ und } a=M_Y - b\cdot M_X \ . \end{equation*}\] Beide Werte werden als Regressionskoeffizienten bezeichnet; \(b\) dann weiter oft als Regressionsgewicht oder Slope und \(a\) als Achsenabschnitt oder Intercept.
Als Maß der Varianzaufklärung haben wir dann den Anteil der durch die lineare Beziehung aufgeklärten Varianz an der Gesamtvarianz von \(Y\) (vgl. auch Effektstärken in der Varianzanalyse) \[\frac{S^2_{\hat{Y}}}{S_Y^2}=\frac{S^2_{\hat{Y}}}{S^2_{\hat{Y}}+S_E^2}=r^2_{XY}\] hergeleitet. Dies wird als Determinationskoeffizient bezeichnet, der ein Spezialfall des weiter unten eingeführten Maßes \(R^2\) ist.
Als ein Beispiel für die Berechnung einer einfachen, linearen Regression hatten wir folgende Daten verwendet:
Als Regressionsgerade (vgl. auch die nächste Abbildung) ergibt sich: \[Y=-1.461\cdot X+20.051\]
Die Berechnung einer Regression in R erfolgt i.d.R. mit der Funktion
lm(). Zur praktischen Durchführung fassen wir die Variablen
\(X\) und \(Y\) zunächst in einem DataFrame namens
daten zusammen:
daten <- data.frame(X = c(1, 3, 5, 4, 6, 4, 7, 8, 5, 8),
Y = c(19, 18, 14, 13, 15, 12, 9, 10, 9, 7))Dann führen wir die Regression durch und speichern das Ergebnis
inmodell. Eine Zusammenfassung können wir mit der Funktion
summary() anfordern. Eine andere Möglichkeit ist hier die
Funktion S() aus dem Paket car, die mehr
Flexibilität erlaubt. Wir werden in diesem und den folgenden Teilen
immer wieder Funktionen aus diesem Paket benutzen (“car” steht für
companion to applied regression; Fox & Weisberg, 2019):
modell <- lm(Y ~ X, # Berechnung der Regression
data = daten) # Ausgabe mit S() aus dem Paket car
S(modell)## Call: lm(formula = Y ~ X, data = daten)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.051 1.988 10.085 7.97e-06 ***
## X -1.461 0.360 -4.058 0.00364 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 2.412 on 8 degrees of freedom
## Multiple R-squared: 0.6731
## F-statistic: 16.47 on 1 and 8 DF, p-value: 0.003642
## AIC BIC
## 49.76 50.67Auf diesen Grundlagen aufbauend, werden wir in den folgenden Abschnitten nun Aspekte am Beispiel der einfachen, linearen Regression ergänzen, die auch im Fall der multiplen Regression relevant sind, aber Konzepte benötigen, die wir erst nach der Einführung in die Regression in Statistik I kennengelernt haben.
7.1.2 Grafische Überprüfung von Verteilungsannahmen
Bei \(t\)-Tests und Varianzanalysen haben wir Verteilungsannahmen über bestimmte Werte an den Anfang gestellt, um die Verteilung der jeweils hergleiteten Prüfgrößen korrekt bestimmen zu können. Wir werden später sehen, dass dies natürlich auch für Regressionen gemacht wird. Hier wollen wir uns nun zunächst allgemein mit Möglichkeiten beschäftigen, diese Annahmen grafisch zu inspizieren. Dazu verwenden wir als Beispiel eine Stichprobe von \(n=200\) Werten aus einer Standardnormalverteilung:
stichprobe <- round(rnorm(n = 200,
mean = 0,
sd = 1), 2)Eine bereits bekannte Möglichkeit sind Histogramme.
Die folgende Abbildung stellt die gerade generierten Daten als
Histogramm dar. Zusätzlich sind die kleinen Balken an der \(x\)-Achse mit rug()
hinzugefügt worden, deren “Dichte” auch Informationen über die
Häufigkeit des Auftretens der Werte anzeigt. Schließlich wurde noch die
Standardnormalverteilung eingezeichnet, da wir ja wissen, aus welcher
Verteilung die Daten stammen:
hist(stichprobe,
freq = FALSE,
xlim = c(-4, 4),
ylim = c(0, 0.6),
main = "Histogramm")
box() # umrahmt das Diagramm
rug(stichprobe) # fügt die kleinen Balken an der x-Achse hinzu
x <- seq(-4, 4, length.out= 100) # fügt Dichte der...
curve(dnorm(x = x, # ...Standardnormalverteilung hinzu
mean = 0,
sd = 1),
add = TRUE,
col = "black",
lwd = 2) 
Histogramme sind oft hilfreich, aber die Wahl der Kategorienbreite
kann auch Information verzerren. Da wir zudem oft nicht wissen, aus
welcher Verteilung Daten stammen, sondern wir darüber Information
erhalten wollen, können wir auch die “wahre” Verteilung nicht einfach
einzeichnen. Als Alternative hierfür bieten sich sog. Dichteschätzungen
oder Kerndichteplots (Kernel-Density Plots)
an. Dies ist in der folgenden Abbildung illustriert.
density() ist dabei die Standard-R-Funktion (blaue Linie),
adaptiveKernel() ist eine flexiblere Funktion aus dem Paket
car (rote Linie). In etwa kommen beide Funktionen aber zu
einem ähnlichen Ergebnis, welches eine Normalverteilung nahelegt:
hist(stichprobe,
freq = FALSE,
xlim = c(-4, 4),
ylim = c(0, 0.6),
main = "Histogramm mit Kerndichte")
box() # umrahmt das Diagramm
rug(stichprobe) # fügt die kleinen Balken an der x-Achse hinzu
# fügt zwei Kerndichteschätzungen hinzu
lines(density(stichprobe), # Standardkerndichte
col = "blue",
lwd = 2)
lines(adaptiveKernel(stichprobe), # aus car
col = "red",
lwd = 2)
Ein weiteres wichtiges Werkzeug zur Beurteilung von Verteilungen sind
sog. Quantil-Quantil-Plots oder auch Q-Q
Plots genannt, die einen grafischen Vergleich einer empirischen
Verteilung mit einer angenommenen Referenzverteilung (oft ist dies dann
die Normalverteilung) erlauben. Hier werden auf der \(y\)-Achse sortierte empirische Werte gegen
das jeweils theoretisch erwartete Quantil auf der \(x\)-Achse abgebildet. Wenn die Daten zur
Referenzverteilung passen, sollten die Punkte i.W. auf einer Geraden
liegen, ohne systematische Abweichungen in einer Richtung aufzuweisen.
Die einfachste Möglichkeit hierzu bietet die Funktion
qqnorm(), die als Referenzverteilung die Normalverteilung
annimmt (siehe linker Teil der nächsten Abbildung). Die Funktion
qqPlot() aus dem Paket car ist flexibler und
plottet gleich noch einen Toleranzrahmen um die Gerade in welchem die
Punkte weitgehend liegen sollten (siehe rechter Teil der nächsten
Abbildung). Die Option id = list(n=3) gibt dabei an, ob
bzw. wieviele Werte (hier: 3, Standard: 2) markiert werden sollen, die
von der Funktion als “extrem” identifiziert werden:
par(mfrow = c(1, 2))
qqnorm(stichprobe) # einfacher Standard Q-Q Plot
qqPlot(stichprobe, # Q-Q Plot aus car
id = list(n = 3)) 
## [1] 61 171 14Der Vorteil der Funktion qqPlot() ist, dass bei ihr auch
auch andere Referenzverteilungen angenommen werden können, und nicht nur
die Normalverteilung wie – naheliegenderweise – bei
qqnorm(). Um dies zu demonstrieren ziehen wir im nächsten
Beispiel zunächst Daten aus einer Gleichverteilung bei der, etwas
vereinfacht gesagt, alle Werte zwischen 0 und 10 gleich wahrscheinlich
sind. Im linken Teil der Abbildung wird dann eine Normalverteilung als
Referenz benutzt (was der Standard ist). Hier kann man gut sehen, dass
die Punkte nicht auf einer Geraden liegen, sondern systematisch
abweichen. Im rechten Teil der Abbildung ist dann eine Gleichverteilung
als Referenz angenommen worden und hier zeigt sich die erwartete Gerade.
Im Prinzip kann jede in R implementierte Verteilung als Referenz
verwendet werden:
stichprobe2 <- runif(n = 200, # aus Gleichverteilung ziehen
min = 0,
max = 10)
par(mfrow = c(1, 2))
qqPlot(stichprobe2) # Referenz = Normalverteilung## [1] 88 63qqPlot(stichprobe2,
distribution = "unif") # Referenz = Gleichverteilung
## [1] 88 637.1.3 Annahmen der Regression und deren grafische Überprüfung
Um valide Schlussfolgerungen ziehen zu können, hatten wir bei \(t\)-Tests und der Varianzanalyse jeweils Annahmen gemacht, die für die Daten gelten sollen. Im Fall unabhängiger Stichproben waren diese bspw.:
- die \(J\) Stichproben sind unabhängig voneinander gezogen worden
- die abhängige Variable \(Y\) ist mindestens intervallskaliert
- die abhängige Variable ist in jeder der \(J\) Populationen normalverteilt mit gleicher Varianz \(\sigma^2\) (Varianzhomogenität), wobei wir aber verschiedene Erwartungswerte \(\mu_j\) zulassen: \[Y_j\sim N(\mu_j,\sigma^2)\]
Ganz ähnliche Annahmen haben wir bisher im Fall einfacher linearer Regression stillschweigend gemacht, aber nie explizit eingeführt:
- das Kriterium \(Y\) ist mindestens intervallskaliert
- der Zusammenhang von \(Y\) und \(X\) ist linear
- die \(n\) Beobachtungen sind unabhängig voneinander
- für die Residuen werden darüber hinaus folgende Annahmen gemacht:
- bedingte Varianzhomogenität: für jeden Wert von \(X\) haben die Residuen die gleiche Varianz
- der Erwartungswert der Residuen ist 0
- bedingte Normalverteilung (wichtig für Signifikanztests und Konfidenzintervalle): für jeden Wert von \(X\) sind die Residuen normalverteilt \[\epsilon_m \sim N(0,\sigma^2_{\epsilon})\]
Diese Annahmen müssen teilweise natürlich im Rahmen der Datenerhebung
sichergestellt werden, teilweise können wir aber gerade die
Verteilungsannahmen auch grafisch überprüfen. Dazu stellt R einige
diagnostische Abbildungen bereit. Zur Illustration simulieren wir uns
zunächst zwei Variablen in einem Datensatz (mit \(n = 100\) Fällen), die eine Korrelation von
\(r = 0.7\) haben. Eine Möglichkeit,
dies zu tun, bietet die Funktion mvrnorm() aus dem Paket
MASS:
library(MASS)
set.seed(1) # um reproduzierbare Daten zu simulieren
mu <- c(0, 0) # Erwartungswerte sollen Null sein
Sigma <- matrix(c(1, 0.7, # Kovarianzmatrix bauen
0.7, 1),
nrow = 2,
ncol = 2,
byrow = TRUE)
Sigma # Kovarianzmatrix anzeigen## [,1] [,2]
## [1,] 1.0 0.7
## [2,] 0.7 1.0rohwerte <- mvrnorm(n = 100, # wieviele Datenpunkte?
mu = mu, # Erwartungswerte
Sigma = Sigma, # Kovarianzmatrix
empirical = TRUE) # wenn TRUE, dann spezifizieren mu und
# Sigma die empirischen Werte, sonst die
# der Population
daten <- data.frame(var1 = rohwerte[,1],
var2 = rohwerte[,2]) # in DataFrame wandeln
cor_out(cor.test(daten$var1,
daten$var2)) # stimmt die Korrelation?## Text
## 1 r(100) = .70, p < .001Der erste Schritt besteht nun in der Visualisierung als Scatterplot,
um eine Idee des Zusammenhangs zu bekommen. Statt mit der einfachen
Funktion plot(), führen wir hier die Funktion
scatterplot() aus dem Paket car ein. Für diese
Funktion wird wieder die Modellsprache verwendet. Die dicke
durchgezogene Linie ist dabei die Regressionsgerade, während die dicke
gestrichelte Linie eine non-parametrische Schätzung der Mittelwerte
bietet und frei von Annahmen wie Linearität o.ä. ist. Die beiden äußeren
Linien sind non-parametrische Schätzungen der Standardabweichungen der
\(Y\)-Werte zu jedem Wert von \(X\). Außen sind zusätzlich marginal
Boxplots zu finden, die Informationen über die univariaten
Verteilungen beider Variablen bieten:
scatterplot(var2 ~ var1, # aus car
data = daten,
xlab = "Variable 1",
ylab = "Variable 2") 
Eine lineare Beziehung beider Variablen ist also gegeben und sie eignen sich für eine lineare Regression:
modell <- lm(var2 ~ var1,
data = daten)
brief(modell) # die ganz kurze Zusammenfassung## (Intercept) var1
## Estimate 1.89e-17 0.7000
## Std. Error 7.18e-02 0.0721
##
## Residual SD = 0.718 on 98 df, R-squared = 0.49Eine Reihe diagnostischer Abbildung erhält man nun schnell, wenn die
Funktion plot() auf das lm()-Objekt angewendet
wird:
par(mfrow = c(2, 2)) # 2 x 2 Abbildungen
plot(modell)
- Das erste Diagramm (
Residuals vs Fitted) zeigt die vorhergesagten Werte auf der \(x\)-Achse und die Residuen auf der \(y\)-Achse. Die rote Linie wird auch “Lowess-Linie” genannt und ist eine non-parametrische Anpassung des Zusammenhangs beider Variablen. Wenn das Modell richtig spezifiziert ist und Linearität vorliegt, sollten die Punkte unsystematisch verteilt sein und die rote Linie möglichst parallel zur \(x\)-Achse verlaufen. - Das zweite Diagramm (
Q-Q Residuals) ist ein Q-Q-Plot der (standardisierten) Residuen zur Überprüfung der Normalverteiltheitsannahme (siehe oben). - Das dritte Diagramm (
Scale-Location) kann zur Überprüfung der Varianzhomogenität herangezogen werden. Liegt diese vor, sollten die Punkte unsystematisch verteilt sein. - Das letzte Diagramm (
Residuals vs Leverage) kann zur Identifikation von Ausreißern bzw. einflussreichen Werten herangezogen werden. Hierauf wird im nächsten Abschnitt genauer eingegangen.
7.1.4 Einflussreiche Datenpunkte, Leverage
Wann hat ein Datenpunkt einen besonderen Einfluss auf eine
Regressionsgerade bzw. die -koeffizienten? Um sich dem zu nähern ist es
zunächst sinnvoll zu verstehen, dass der Einfluss neuer Datenpunkte auf
eine Regressionsgerade nicht immer gleich ist. In der folgenden
Abbildung stellt die schwarze Linie die Regressionsgerade für die
schwarzen Datenpunkte dar. Mit \(R^2=0.88\) ist hier eine relativ große
Varianzaufklärung erreicht worden. Der rote Punkt hat eine Koordinate
auf der \(x\)-Achse, die dem Mittelwert
der \(x\)-Koordinaten der schwarzen
Punkte entspricht (\(M=6\)). In einem
solchen Fall wird die Steigung der resultierenden Regressionsgeraden
nicht verändert, wohl aber der Achsenabschnitt. Der Anteil aufgeklärter
Varianz wird darüber hinaus geringer (\(R^2=0.55\)). Die Hinzunahme des blauen
Punktes hingegen verändert (auch) die Steigung.
Wie stark der Einfluss ist hängt von zweien seiner Merkmale ab:
- Der Position auf der \(x\)-Achse, die oft als “Leverage” bezeichnet wird.
- Dem Residuum, also wie weit der neue Punkt von der Regressionsgerade entfernt ist.
Ein Wert der weit vom Mittelwert auf der \(x\)-Achse entfert liegt, muss aber gar keinen großen Einfluss haben, wenn er direkt auf der Regressionsgeraden liegt, während einer nahe am Mittelwert einen großen Einfluss haben kann, wenn das Residuum groß ist.
Eine Möglichkeit der Messung von Leverage sind sog. “Hat-Values” (auch “Hebelwerte” genannt): \[h_i=\frac{1}{N} + \frac{(x_i-M_X)^2}{\sum_{i=1}^N (x_i-M_X)^2 }\] Das Minimum dieser Werte ergibt sich aus der Anzahl der eingehenden Beobachtungen und es gilt \(h_i\geq \frac{1}{N}\). Zudem summieren sich alle Werte im Falle der einfachen linearen Regression zu 2 auf, also \(\sum_{i=1}^N h_i = 2\). Idealerweise sollten alle \(h_i\) ähnlich sein und einzelne hohe Hebelwerte deuten auf einzelne einflussreiche Datenpunkte hin (was aber nicht sein muss, da ja auch das Residuum beachtet werden muss; vgl. unten: Cook’s Distance).
Die Berechnung in R erfolgt sehr leicht, indem die Funktion
hatvalues() auf ein lm()-Objekt angewendet
wird:
modell <- lm(y ~ x) # lineares Modell berechnen
h <- hatvalues(modell) # hat-values berechnen...
h # ...und ausgeben## 1 2 3 4 5 6 7
## 0.4269481 0.1883117 0.1883117 0.1542208 0.1542208 0.1883117 0.6996753(Gesamt-)Maße des Einflusses kombinieren Leverage und Residuen in ein
Maß. Das bekannteste davon ist Cook’s Distance: \[D_i=\frac{e_i^2}{(q+1)\cdot MS_{error}}\cdot
\left(\frac{h_i}{(1-h_i)^2} \right)\] Hier ist \(q\) die Anzahl der Prädiktoren im
Regressionsmodell; bei einer einfachen linearen Regression also \(q=1\). Mit \(MS_{error}\) ist eine Schätzung der
Residuenvarianz gemeint. Auch hier sind absolute Werte schwer zu
interpretieren, denn Cook’s Distance hängt u.a. ab von der
Stichprobengröße und Anzahl der Prädiktoren. Die Untersuchung der Daten
auf einflussreiche Datenpunkte sollte daher nicht von starren
Ausschlusskriterien isoliert erfolgen, sondern stets eingebettet in die
Gesamtmenge der Daten erfolgen, wofür sich v.a. graphische Abbildungen
anbieten (bspw. mit Hilfe der Funktion influenceIndexPlot,
siehe unten). Dennoch gibt es ein paar grobe Empfehlungen. So ist eine
häufige Empfehlung, dass ein \(D\geq
1\) auf einen einflussreichen Datenpunkt hinweist (aber manchmal
wird auch \(D\approx 0.5\) bereits als
einflussreich aufgefasst). Andere empfehlen hingegen den Grenzwert an
den Daten zu orientieren und schlagen vor einen Wert als einflussreich
anzusehen, wenn \(D_i > \frac{4}{n - q -
1}\) (mit \(n\) als
Stichprobengröße und \(q\) als Anzahl
der Prädiktoren).
Die Berechnung von Cook’s Distance anhand der Formel erfolgt durch
# um leicht an SS_error zu kommen:
anova.modell <- Anova(modell)
# dann: MS_error = SS_error / df_error
MS.e <- anova.modell$`Sum Sq`[2] / anova.modell$`Df`[2]
# Umsetzung der Formel
q = 1
D <- (modell$residuals^2 / ((q+1) * MS.e ) ) * (h / ((1-h)^2) )
D## 1 2 3 4 5 6 7
## 0.02790906 0.15980277 0.18676730 0.03975840 0.20924749 0.01823899 0.36611170oder einfach mit der entsprechenden Funktion:
cooks.distance(modell) # Cook's Distance## 1 2 3 4 5 6 7
## 0.02790906 0.15980277 0.18676730 0.03975840 0.20924749 0.01823899 0.36611170Eine grafische Übersicht über Hat-Values und Cook’s Distance (und
weitere diagnostische Maße) kann mit
influenceIndexPlot(modell) aus dem Paket car
angefordert werden:
Auf der x-Achse sind die Indizes der Datenpunkte aufgetragen.
Beispielweise bezieht sich der Index 1 auf den ersten Wert des
Prädiktors bzw. Kriteriums (quasi y[1] und
x[1]).
Der erste Plot gibt Cook’s Distance eines jeden Datenpunktes an (vgl. hierzu die Ausgabe der Funktion
cooks.distance()).Der zweite Plot zeigt die sog. studentisierten Residuen an (nach dem Pseudonym „Student“ des Statistikers William Sealy Gosset). Hierbei werden die Residuen derart standardisiert, dass sie unter Gültigkeit der Regressionsannahmen einer \(t\)-Verteilung mit \(n-q-2\) Freiheitsgraden folgen.
Der dritte Plot bestimmt für jedes Residuum einen \(p\)-Wert anhand der theoretischen \(t\)-Verteilung. Da hierbei allerdings mehrfach getestet wird, sind die berechneten \(p\)-Werte Bonferroni korrigiert (mit einer Begrenzung auf \(p\leq1\), falls die Bonferroni-Korrektur \(p\)-Werte größer als 1 ergibt). Ist der korrigierte \(p\) statistisch signifikant, so ist der entsprechende Datenpunkt, relativ zur Regressionsgerade gesehen, vermutlich ein Ausreißer.
Der letzte Plot gibt die Hebelwerte aller Datenpunkte an (vgl. hierzu die Ausgabe der Funktion
hatvalues()).
Eine Shiny-App zur interaktiven Visualisierung der Konzepte Leverage und Cook’s Distance sowie der Auswirkung einzelner Ausreißer auf die Regressionsgerade ist hier zu finden.
7.1.5 Konfidenz- und Prädiktionsintervalle
In Teil
13 von Statistik I haben wir uns auch bereits mit
Konfidenzintervallen befasst, wobei wir den Schwerpunkt auf Mittelwerte
und deren Unterschiede gelegt haben. Prinzipiell können aber für alle
Werte, deren Verteilung bekannt ist, Konfidenzintervalle berechnet
werden. Für die Regressionskoeffizienten werden in der
Ausgabe von summary() bzw S() nur
Signifikanztests angegeben. Mit der Funktion confint()
können aber sehr einfach für die Koeffizienten Konfidenzintervalle
angefordert werden:
confint(modell) # standardmäßig 95%-KI...## 2.5 % 97.5 %
## (Intercept) 7.6669411 11.1382537
## x -0.8399515 -0.3418667confint(modell,
level = 0.99) # ...oder mit level = XX einstellen## 0.5 % 99.5 %
## (Intercept) 6.6800953 12.1250995
## x -0.9815501 -0.2002681Wenn die Regressionskoeffizienten bekannt sind, kann die Regressionsgleichung zur Vorhersage neuer Werte benutzt werden. Dies war als eines der wichtigsten Ziele der Regressionsrechnung eingeführt worden. Natürlich ist diese Vorhersage mit einer gewissen Unsicherheit behaftet und auch hier kann ein Intervall angegeben werden, welches die Unsicherheit darstellt. Allerdings wird hier zwischen zwei Intervallen unterschieden, wobei der Unterschied zwischen beiden Intervallen von der unterschiedlichen Berechnung der Standardfehler herrührt:
- Vorhersage der bedingten Erwartungswerte zu jedem Wert von \(X\) (quasi die Regressionsgerade): Hierfür wird ein Konfidenzintervall berechnet, welches – wegen der unterschiedlichen Leverage-Werte – nicht parallel zur Regressionsgeraden verläuft.
- Vorhersage individueller (neuer) Kriteriumswerte: hierfür wird ein Prädiktionsintervall berechnet, welches größer als das Konfidenzintervall ist.
Zur Berechnung (und vor allem auch zur Visualisierung) können wir die
Funktion predict() verwenden. Einfach angewendet auf ein
lm()-Objekt berechnet diese Funktion die vorhergesagten
Werte für die Werte des Prädiktors:
predict(modell)## 1 2 3 4 5 6 7
## 8.811688 7.038961 7.038961 6.448052 5.266234 4.675325 1.720779Wenn der Funktion aber noch ein DataFrame mit einer Variablen
übergeben wird, deren Name dem Prädiktor entspricht, dann werden die
vorhergesagten Werte samt Intervall berechnet. Der Trick hierbei ist,
dass mit der Funktion seq() ein ganzer Wertebereich
abgedeckt wird. Mit dem Argument interval wird dann
angegeben, ob ein Konfidenz- oder ein Prädiktionsintervall berechnet
werden soll:
modell <- lm(y ~ x) # Modell berechnen
newx <- seq(from = 0, # neue x-Werte
to = 15,
by = 0.01)
conf_interval <- predict(modell,
newdata = data.frame(x = newx),
interval = "confidence",
level = 0.95)
pred_interval <- predict(modell,
newdata = data.frame(x = newx),
interval = "prediction",
level = 0.95)Die Funktion erzeugt eine Matrix mit drei Spalten. Die erste Spalte enthält die vorhergesagten Werte, die zweite und dritte Spalte die Werte der unteren und oberen Grenze der Intervalle:
head(conf_interval) # Anzeigen der ersten sechs Zeilen## fit lwr upr
## 1 9.402597 7.666941 11.13825
## 2 9.396688 7.663176 11.13020
## 3 9.390779 7.659409 11.12215
## 4 9.384870 7.655642 11.11410
## 5 9.378961 7.651873 11.10605
## 6 9.373052 7.648104 11.09800Mit dieser Matrix können die Intervalle einfach geplottet werden. Im
oberen Teil der folgenden Abbildung sind die Intervalle als Kurven
eingezeichnet, während sie unten als durchsichtige Flächen dargestellt
sind (was mit der Funktion polygon() und der Farbwahl über
die Funktion rgb() funktioniert). Im linken Teil sind
jeweils Konfidenzintervalle zu sehen, rechts Prädiktionsintervalle:
par(mfrow = c(2, 2))
# Konfidenzintervall
plot(x, y,
xlab = "X",
ylab = "Y",
main = "Konfidenzintervall",
ylim = c(0, 10),
xlim = c(0, 15),
pch = 19)
lines(newx, conf_interval[,1], col = "black", lty = 1)
lines(newx, conf_interval[,2], col = "blue", lty = 2)
lines(newx, conf_interval[,3], col = "blue", lty = 2)
# Prädiktionsintervall
plot(x, y,
xlab = "X",
ylab = "Y",
main = "Prädiktionsintervall",
ylim = c(0, 10),
xlim = c(0, 15),
pch = 19)
lines(newx, pred_interval[,1], col = "black", lty = 1)
lines(newx, pred_interval[,2], col = "blue", lty = 2)
lines(newx, pred_interval[,3], col = "blue", lty = 2)
# Konfidenzintervall
plot(x, y,
xlab = "X",
ylab = "Y",
main = "Konfidenzintervall",
ylim = c(0, 10),
xlim = c(0, 15),
pch = 19)
lines(newx, conf_interval[,1], col = "black", lty = 1)
polygon(c(newx, rev(newx)), c(conf_interval[,2], rev(conf_interval[,3]) ),
col = rgb(0.7,0.7,0.7,0.4),
border = FALSE)
# Prädiktionsintervall
plot(x, y,
xlab = "X",
ylab = "Y",
main = "Prädiktionsintervall",
ylim = c(0, 10),
xlim = c(0, 15),
pch = 19)
lines(newx, pred_interval[,1], col="black", lty=1)
polygon(c(newx, rev(newx)), c(pred_interval[,2], rev(pred_interval[,3]) ),
col = rgb(0.7,0.7,0.7,0.4),
border = FALSE)
7.2 Erste Erweiterung: 2 orthogonale Prädiktoren
Von orthogonalen Variablen redet man, wenn zwei Variablen nicht miteinander korreliert sind. Dies ist quasi ein “idealer Spezialfall”. Wir betrachten hier nun ein Beispiel aus dem Buch von Baguley (2012, S. 426ff) mit drei Variablen: \(X_1\) und \(X_2\) werden unsere Prädiktoren sein, \(Y\) das Kriterium. Zunächst erstellen wir uns ein passendes DataFrame mit den Daten:
daten <- data.frame(VP = c(1:7),
X1 = c(-3, -2, -1, 0, 1, 2, 3),
X2 = c(2, 2, 0, 4, 5, 1, 1),
Y = c(-7.15, 0.02, -1.27, 0.68, 11.44, 9.43, 3.29))
round(cor(daten[,2:4]), 2) # Korrelationsmatrix der Spalten 2-4 in daten## X1 X2 Y
## X1 1.00 0.00 0.76
## X2 0.00 1.00 0.34
## Y 0.76 0.34 1.00Sowohl \(X_1\) als auch \(X_2\) sind mit \(Y\) korreliert, während \(X_1\) und \(X_2\) selber unkorreliert, d.h. orthogonal zueinander, sind. Die übergeordnete Frage ist nun: Wie kann man aus den beiden Prädiktoren \(X_1\) und \(X_2\) das Kriterium \(Y\) möglichst gut vorhersagen? Die Erweiterung des Modells der einfachen linearen Regression \[\hat{Y}=b\cdot X+a\] sieht konsequenterweise so aus: \[\hat{Y}=b_1\cdot X_1+b_2\cdot X_2 + a\] Ganz ähnlich ist tatsächlich auch das weitere Vorgehen: Gesucht werden nun drei Werte \(b_1\), \(b_2\) und \(a\) mit denen die Summe der quadrierten Abweichungen der vorhergesagten von den empirischen Werten (also der Residuen) minimal wird. Um es vorwegzunehmen: Dies ist im Regelfall nicht mehr ganz einfach; aber deshalb schauen wir uns auch zunächst zwei orthogonale Prädiktoren an.
Während bei einer einfachen linearen Regression die Regressionsgerade eine Gerade im zweidimensionalen Raum spezifiziert hat, wird bei der multiplen Regression mit zwei Prädiktoren eine Ebene im dreidimensionalen Raum spezifiziert. Diese kann über einen 3D-Scatterplot visualisiert werden, wobei in der folgenden Abbildungen bereits die “optimale” Ebene eingezeichnet ist:
library(scatterplot3d)
s3d <- scatterplot3d(daten$X1,
daten$X2,
daten$Y,
xlab = "X1",
ylab = "X2",
zlab = "Y",
pch = 16,
color = "red",
highlight.3d = FALSE,
type = "h",
main="3D Scatterplot",
grid = TRUE,
box = FALSE)
s3d$plane3d(modell)
Auch in diesem Fall können Residuen ganz ähnlich aufgefasst werden: Die
Abstände der empirischen Datenpunkte zur Ebene (statt zur Geraden). Wie
erhalten wir nun die gesuchten Werte für die Regressionskoeffizienten?
Die Lösung mit R ist tatsächlich sehr einfach und sie erfordert
tatsächlich nur eine minimale Ergänzung der Modellgleichung, indem
(additiv) ein weiterer Prädiktor hinzugefügt wird:
modell <- lm(Y ~ X1 + X2,
data = daten)
S(modell)## Call: lm(formula = Y ~ X1 + X2, data = daten)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2823 2.7043 -0.104 0.9219
## X1 2.2446 0.8222 2.730 0.0524 .
## X2 1.2277 1.0019 1.225 0.2876
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 4.351 on 4 degrees of freedom
## Multiple R-squared: 0.6912
## F-statistic: 4.478 on 2 and 4 DF, p-value: 0.09533
## AIC BIC
## 44.53 44.32Erkennbar ist, dass es neben dem Intercept (dem Achsenabschnitt \(a=-0.2823\); also demjenigen Wert, der für
\(X1=X2=0\) vorhergesagt wird) nun zwei
weitere Koeffizienten gibt, nämlich \(b_1=2.2446\) und \(b_2=1.2277\). Sind wir nur an den
Koeffizienten und ihren Standardfehlern interessiert, könenn wir auch
die Funktion brief() zur Ausgabe nutzen:
brief(modell)## (Intercept) X1 X2
## Estimate -0.282 2.245 1.23
## Std. Error 2.704 0.822 1.00
##
## Residual SD = 4.35 on 4 df, R-squared = 0.691Kleine Denkaufgabe: Welchen Wert erhalten Sie, wenn
Sie die Zeile estimate durch die Zeile
Std. Error dividieren?
Wie diese Koeffizienten im allgemeinen Fall berechnet werden, ist nicht ganz trivial (in Teil 12 von Statistik II werden wir kurz darauf zurück kommen, wenn wir das “Allgemeine Lineare Modell” in seinen Grundzügen kennenlernen). Im Fall orthogonaler Prädiktoren reduziert sich das Problem allerdings auf zwei einfache lineare Regressionen, da in diesem Fall die Prädiktoren distinkte Teile der Gesamtvariabilität der Daten binden. Dies können wir sehr leicht nachvollziehen, indem wir zwei einfache lineare Regressionen berechnen:
modell.X1 <- lm(Y ~ X1, # X1 als Prädiktor
data = daten)
brief(modell.X1)## (Intercept) X1
## Estimate 2.35 2.245
## Std. Error 1.72 0.862
##
## Residual SD = 4.56 on 5 df, R-squared = 0.575modell.X2 <- lm(Y ~ X2, # X2 als Prädiktor
data = daten)
brief(modell.X2) ## (Intercept) X2
## Estimate -0.282 1.23
## Std. Error 4.093 1.52
##
## Residual SD = 6.58 on 5 df, R-squared = 0.116Den Achsenabschnitt \(a\) können wir in diesem Fall berechnen als \[a=M_Y-(b_1\cdot M_{X_1}+b_2\cdot M_{X_2})\] Bezogen auf die Beispieldaten ergibt sich also:
b1 <- coef(modell.X1)[2] # Extrahieren der beiden Koeffizienten
b2 <- coef(modell.X2)[2]
a <- mean(daten$Y) - ( b1 * mean(daten$X1) + b2 * mean(daten$X2))
a## X1
## -0.2822727Die so berechneten Koeffizienten entsprechen also denen, die zu Beginn ausgegeben wurden, als beide Prädiktoren in das Modell aufgenommen worden sind.
Auch das Maß \(R^2\) kann an dieser Stelle noch einmal anders betrachtet werden. Im Fall der einfachen linearen Regression hatten wir festgestellt, dass \(R^2\) gleich der quadrierten Korrelation zwischen \(X\) und \(Y\) ist. Im Fall mehrerer Prädiktoren berechnen wir zunächst die vorhergesagten Werte auf Basis der Regressionskoeffizienten; hier für den Fall unserer zwei Prädiktoren: \[\hat Y=b_1\cdot X_1+b_2\cdot X_2 + a=2.2446 \cdot X_1 + 1.2277\cdot X_2-0.2823\] Da wir die beiden Koeffizienten gerade bereits extrahiert haben, können wir \(\hat Y\) sehr einfach selbst berechnen und als neue Variable an den DataFrame anhängen:
daten$y_dach <- b1 * daten$X1 + b2 * daten$X2 + aNun berechnen wir die quadrierte Korrelation zwischen \(Y\) und \(\hat Y\):
cor(daten$Y, daten$y_dach)^2## [1] 0.6912412Dies ist genau der Wert \(R^2\), den
auch die Zusammenfassung von modell ergeben kann: \(R^2\) (bzw. genauer \(R\), die Wurzel aus \(R^2\)) ist also selber eine Korrelation,
nämlich die zwischen empirischen und vorhergesagten Werten; der
multiple Korrelationskoeffizient. Dass zwei orthogonale
Prädiktoren jeweils verschiedene Varianzanteile aufklären, sehen wir
auch daran, dass \(R^2\) hier auch die
Summe der quadrierten Korrelationen zwischen \(Y\) und den beiden Prädiktoren einzeln ist:
\[R^2=r_{Y,X_1}^2 + r_{Y,X_2}^2\]
R2 <- cor(daten$Y, daten$X1)^2 + cor(daten$Y, daten$X2)^2
R2## [1] 0.69124127.3 Allgemeiner Ansatz bei intervallskalierten Prädiktoren und Kriterium
Wir haben gerade den Spezialfall (zweier) orthogonaler Prädiktoren betrachtet und einige anschauliche Zusammenhänge erarbeitet. In der Realität ist die Unkorreliertheit der Prädiktoren allerdings eher selten der Fall, wenngleich die grundlegende Herangehensweise auch bei korrelierten Prädiktoren identisch ist. Auch die generellen Annahmen und Diagnostik der Annahmen sind wie im Fall der einfachen linearen Regression (die letztlich nur einen Spezialfall multipler Regression darstellt). Im Folgenden betrachten wir nun den allgemeinen Ansatz und die Besonderheiten, die sich für die Interpretation der Koeffizienten ergeben.
7.3.1 Datenbeispiel
Wir benutzen hier ein (minimal modifiziertes) Datenbeispiel aus dem
Buch von Field, Miles und Field (2012), welches im Datenordner unter dem
Namen Album_Verkaeufe.dat gespeichert ist. Dieser Datensatz
enthält vier Variablen:
daten <- read.table("./Daten/Album_Verkaeufe.dat",
header = TRUE)
head(daten)Die erste Variable Werbung gibt an, wieviel Euro eine
Plattenfirma in die Werbemaßnahmen eines neuen Musikalbums investiert
hat (in “Tausend Euro”-Einheiten). Die zweite Variable
Verkaeufe gibt an angibt, wie oft das Album (in
Tausendereinheiten) verkauft wurde. Die dritte Variable
Frequenz gibt an, wie häufig Stücke vom Album in der Woche
vor der eigentlichen Veröffentlichung im Radio gespielt wurden und die
vierte Variable Attraktivitaet gibt an, wie optisch
ansprechend die Band auf einer Skala von 1-10 eingeschätzt wurde.
Wollen wir nun aus den Werbeausgaben die Verkäufe vorherzusagen, ist dies ein Fall für eine einfache, lineare Regression. Wir berechnen diese Regression daher und visualisieren die Daten und die Regressionsgerade:
modell1 <- lm(Verkaeufe ~ Werbung,
data = daten)
plot(daten$Werbung, # Scatterplot
daten$Verkaeufe,
pch = 19,
xlab = "Werbung",
ylab = "Verkäufe",
cex = 0.5)
abline(modell1, # Regressionsgerade hinzufügen
lwd = 2,
col = "blue")
S(modell1) # Ausgabe der Regression## Call: lm(formula = Verkaeufe ~ Werbung, data = daten)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.341e+02 7.537e+00 17.799 <2e-16 ***
## Werbung 9.612e-02 9.632e-03 9.979 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 65.99 on 198 degrees of freedom
## Multiple R-squared: 0.3346
## F-statistic: 99.59 on 1 and 198 DF, p-value: < 2.2e-16
## AIC BIC
## 2247.38 2257.27Was bedeutet das Ergebnis nun?
- Zunächst ist im Scatterplot erkennbar, dass die Anzahl der Verkäufe tatsächlich mit den investierten Werbeeuros ansteigt. Zwar gibt es auch Alben, die sich fast ohne Werbung gut verkaufen (links oben im Scatterplot), aber der generelle Trend zeigt den Anstieg.
- \(R^2=0.33\) bedeutet, dass etwa 33% der Varianz der Verkäufe durch die Varianz der Werbeeuros erklärt werden können. Das bedeutet aber auch, dass etwa 67% der Varianz unerklärt bleiben, was wiederum die Beteiligung weiterer Prädiktoren nahelegt.
- Kommen wir nun zu den Regressionskoeffizienten:
- Als Achsenabschnitt ergibt sich \(a=134.1\). Das bedeutet, wenn keine Werbeeuros investiert werden (der Prädiktor also den Wert Null annimmt), würden dennoch etwa 134100 Alben verkauft werden (die Maßeinheit war ja in Tausendern!)
- Die Steigung der Geraden wurde mit \(b=0.096\) berechnet. Dies kann interpretiert werden als Änderung des Kriteriums, wenn der Prädiktor um 1 erhöht wird. Wenn also der Prädiktor “Budget für Werbung” um einen Schritt (d.h. hier: um 1000 Euro) erhöht wird, dann steigt das Kriterium “Anzahl Verkäufe” um 0.096 (d.h. hier: 96 verkaufte Alben). Also: Für 1000 Euro mehr an Werbung gibt es “nur” 96 Verkäufe mehr. Allerdings ist der Beitrag dieses Prädiktor “signifikant” in dem Sinne, als dass der Signifikanztests für die Alternativhypothese \(H_1: \beta\neq 0\) entschieden hat. (Dies drückt sich hier auch dadurch aus, dass der \(F\)-Test den gleichen \(p\)-Wert ergibt, wie der Test der Steigung. Darauf kommen wir später nocheinmal zurück.)
Wir werden diesen Datensatz später weiter verwenden.
7.3.2 q-viele Prädiktoren
Generell kann die bisher verwendete Regressionsgleichung auf beliebig viele additiv zusammengefügte Prädiktoren erweitert werden. Im Fall von \(q\) Prädiktoren wäre die Regressionsgleichung dann also: \[\hat{Y}=a + b_1\cdot X_1 + b_2\cdot X_2 + \ldots + b_q\cdot X_q=a+\sum_{p=1}^q b_p X_p\] Hier sind nun zwei Punkte festzuhalten:
- Die Bestimmung der Koeffizienten ist im allgemeinen Fall nicht mehr als Problem einfacher linearer Regressionen aufzufassen (dies gilt nur bei orthogonalen Prädiktoren). Während für R die Lösung dennoch sehr einfach geht (und wir kommen in Teil 12 kurz darauf zurück), ist eine einfache Lösung von Hand nicht möglich.
- Im Fall eines Prädiktors kann die beste Lösung als Gerade dargestellt werden und bei zwei Prädiktoren kann – wie oben gezeigt – die Lösung als Ebene dargestellt werden. Solche grafischen Illustrationen sind natürlich bei mehr als zwei Prädiktoren dann nicht mehr möglich, da wir dann (mindestens) vier-dimensionale Darstellungen benötigen würden.
7.3.3 Interpretation der Regressionsgewichte
Wir greifen nun wieder das Datenbeispiel auf. Es ist auch
naheliegend, dass ein Album umso häufiger verkauft wird, je häufiger es
vor dem Erscheinen gespielt wurde (Variable Frequenz). Wir
berechnen daher zunächst die einfache, lineare Regression mit
Frequenz als Prädiktor:
modell.frequenz <- lm(Verkaeufe ~ Frequenz,
data = daten)
# Scatterplot
plot(daten$Frequenz,
daten$Verkaeufe,
pch = 19,
xlab = "Frequenz",
ylab = "Verkäufe",
cex = 0.5)
abline(modell.frequenz,
lwd = 2,
col = "blue")
S(modell.frequenz)## Call: lm(formula = Verkaeufe ~ Frequenz, data = daten)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 84.8725 11.2670 7.533 1.73e-12 ***
## Frequenz 3.9392 0.3743 10.524 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 64.79 on 198 degrees of freedom
## Multiple R-squared: 0.3587
## F-statistic: 110.7 on 1 and 198 DF, p-value: < 2.2e-16
## AIC BIC
## 2240.01 2249.91Offensichtlich werden also Alben auch umso häufiger verkauft, je öfters sie vor dem Erscheinen bereits gespielt wurden. Nun betrachten wir beide Prädiktoren gemeinsam, stellen aber zunächst fest, dass beide miteinander korreliert sind, sie also nicht orthogonal sind:
cor(daten$Werbung, daten$Frequenz)## [1] 0.1018828Grafisch lassen sich die Daten wieder als 3D-Scatterplot illustrieren
und die bereits eingezeichnete Regressionsebene verdeutlicht den
Zusammenhang der Verkaufszahlen zu beiden Prädiktoren bereits. Die
Verkaufszahlen nehmen mit zunehmenden Werbeeuros aber auch mit der
Frequenz des Abspielens vor dem Erscheinen zu:
Nun berechnen wir die multiple Regression mit den beiden Prädiktoren
Werbung und Frequenz:
modell.2 <- lm(Verkaeufe ~ Werbung + Frequenz,
data = daten)
summary(modell.2)##
## Call:
## lm(formula = Verkaeufe ~ Werbung + Frequenz, data = daten)
##
## Residuals:
## Min 1Q Median 3Q Max
## -112.121 -30.027 3.952 32.072 155.498
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 41.123811 9.330952 4.407 1.72e-05 ***
## Werbung 0.086887 0.007246 11.991 < 2e-16 ***
## Frequenz 3.588789 0.286807 12.513 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 49.38 on 197 degrees of freedom
## Multiple R-squared: 0.6293, Adjusted R-squared: 0.6255
## F-statistic: 167.2 on 2 and 197 DF, p-value: < 2.2e-16Zunächst einmal sehen wir, dass die aufgeklärte Varianz des Kriteriums mit \(R^2\)=0.63%, also etwa 63%, deutlich gegenüber der einfachen, linearen Regression mit Werbungseuros als alleinigem Prädiktor zugenommen hat. Wir werden auf dieses und andere Maße später nocheinmal zu sprechen kommen.
Als Regressionskoeffizienten haben sich nun ergeben:
- Intercept \(a=41.12\)
- Werbungseuros \(b_1=0.09\)
- Frequenz \(b_2=3.59\)
Die Interpretation der Regressionskoeffizienten \(b_1\) und \(b_2\) unterscheidet sich nun aber von der im Fall nur eines Prädiktors: Sie können nun nur als Regressionsgewicht bedingter einfacher linearer Regressionen aufgefasst werden, d.h. einfacher linearer Regressionen wenn der Wert des jeweils anderen Prädiktors konstant gehalten wird. Diese Interpretation generalisiert auch auf den Fall von mehr als zwei Prädiktoren: Jedes Regressionsgewicht gibt dann die Veränderung des Kriteriums in Abhängigkeit von der Veränderung des Prädiktors an, wenn alle anderen Prädiktoren konstant gehalten werden.
Darüber hinaus gibt es noch eine zweite Interpretation, nämlich die
als Regressionsgewichte von Regressionsresiduen. Die beiden
Regressionsresiduen erhalten wir, indem sowohl das Kriterium als auch
der fragliche Prädiktor vom linearen Einfluss aller anderen Prädiktoren
“befreit” werden. Dies passiert, indem beide Größen als Kriterien einer
Regression behandelt werden und dann deren Residuen weiter verwendet
werden. Einen ganz ähnlichen Gedankengang hatten wir bereits bei der
Partialkorrelation, die wir als Korrelation von Regressionsresiduen
aufgefasst hatten. Bezogen auf unseren Fall heißt das für das
Regressionsgewicht des Prädiktors Werbungseuros (\(b_1=0.09\)), das wir sowohl
Verkaeufe als auch Werbung durch
Frequenz vorhersagen und dann die Residuen beider
Regressionen in einem neuen Modell verwenden:
p1 <- lm(Verkaeufe ~ Frequenz,
data = daten)
p2 <- lm(Werbung ~ Frequenz,
data = daten)
m <- lm(p1$residuals ~ p2$residuals)
brief(m)## (Intercept) p2$residuals
## Estimate 1.18e-14 0.08689
## Std. Error 3.48e+00 0.00723
##
## Residual SD = 49.3 on 198 df, R-squared = 0.422Offensichtlich ergibt sich mit \(0.08689\) der gleiche Wert wie für den
Prädiktor Werbung im Rahmen der multiplen Regression. Wegen
der Verwandtschaft zur Partialkorrelation werden die Regressionsgewichte
in der multiplen Regression auch
(Semi-)Partialregressionsgewichte genannt. (Anmerkung:
Bei einer Partialkorrelation hatten wir beide Variablen vom linearen
Einfluss der auszupartialisierenden dritten Variablen befreit und dann
die Korrelation der Residuen betrachtet. Bei einer
Semipartialkorrelation wird nur eine der beiden Variablen vom linearen
Einfluss der auszupartialisierenden dritten Variablen befreit.)
Die Wirkung eines Prädiktors bei Konstanthaltung der anderen
Prädiktoren kann am besten mit Hilfe sog.
Prädiktor-Effekt-Plots anschaulich visualisiert werden.
Hierzu wird das Paket effects zusätzlich benötigt. Die
Grundidee ist dabei, für jeden Prädiktor eine bedingte Regressionsgerade
zu visualisieren, während alle anderen Prädiktoren auf einen festen Wert
gesetzt werden; der Standardwert ist hierbei der Mittelwert aller
Beobachtungen auf den jeweiligen Prädiktoren. Im einfachsten Fall lassen
sich die Plots dann mit der Funktion predictorEffects()
erstellen (mehr Informationen und Details gibt es hier):
plot(predictorEffects(modell.2),
ylim = c(0,400))
7.3.4 unstandardisierte vs. standardisierte Regressionsgewichte
Bisher haben wir sog. unstandardisierte Regressionsgewichte betrachtet, die auf Grundlage der Metrik, in der die Variablen gemessen wurden, berechnet werden. Darüber hinaus lassen sich standardisierte Regressionsgewichte berechnen, die Veränderungen in der Einheit von Standardabweichungen ausdrücken. Wir bezeichnen die standardisierten Regressionsgewichte als \(\beta\)-Gewichte (aber es finden sich auch andere Bezeichnungen wie z.B. \(b_{1z}\)).
Solche standardisierten Regressionsgewichte erhalten wir z.B. indem
alle beteiligten Variablen \(z\)-standardisiert werden. Wir zeigen dies
hier zunächst noch einmal am Beispiel der einfachen, linearen Regression
mit Verkaeufe als Kriterium und Werbung als
Prädiktor:
modell.1.z <- lm(scale(Verkaeufe) ~ scale(Werbung),
data = daten)
brief(modell.1.z)## (Intercept) scale(Werbung)
## Estimate -2.14e-17 0.578
## Std. Error 5.78e-02 0.058
##
## Residual SD = 0.818 on 198 df, R-squared = 0.335Während der Achsenabschnitt in diesem Fall immer Null ist, ergibt sich ein standardisiertes Regressionsgewicht von \(\beta=0.578\). Im Fall der einfachen, linearen Regression entspricht dies der Korrelation beider Originalvariablen:
cor(daten$Werbung, daten$Verkaeufe)## [1] 0.5784877Ganz ähnlich können wir natürlich auch die standardisierten Regressionsgewichte bei einer multiplen Regression berechnen:
modell.2.z <- lm(scale(Verkaeufe) ~ scale(Werbung) + scale(Frequenz),
data = daten)
brief(modell.2.z)## (Intercept) scale(Werbung) scale(Frequenz)
## Estimate -3.24e-18 0.5229 0.5456
## Std. Error 4.33e-02 0.0436 0.0436
##
## Residual SD = 0.612 on 197 df, R-squared = 0.629Auch in diesem Fall ergibt sich als Achsenabschnitt Null, aber die beiden Regressionsgewichte sind nun so zu interpretieren, dass sie eine Veränderung des Kriteriums in der Einheit Standardabweichung angeben, wenn sich der Prädiktor um eine Standardabweichung verändert und die anderen Prädiktoren konstant gehalten werden. Auch wenn es vielleicht auf den ersten Blick naheliegend ist: die standardisierten Regressionsgewichte entsprechen i.d.R. nicht (semi-)partiellen Korrelationen eines Prädiktors und des Kriteriums. Eine solche Interpretation gilt nur im Fall orthogonaler Prädiktoren!
Eine alternative Berechnung der standardisierten Regressionsgewichte erfolgt unter Berücksichtigung der Standardabweichungen des Kriteriums und des jeweiligen Prädiktors als: \[\beta_p = b_p\cdot\frac{S_{X_p}}{S_Y}\] Im Fall des Prädiktors Werbeeuros kann hier z.B. wie folgt vorgegangen werden:
coef(modell.2)[2] * (sd(daten$Werbung) / sd(daten$Verkaeufe))## Werbung
## 0.5228959Ein häufiges Argument für die Verwendung standardisierter Regressionskoeffizienten ist, wenn man den relativen Einfluss mehrerer Prädiktoren vergleichen will, die in verschiedenen Metriken gemessen wurden (siehe aber z.B. Fox & Weisberg, 2019, S. 187).
7.3.5 Güte des Regressionsmodells
Als Maß der Güte eines Regressionsmodells haben wir bereits \(R^2\) kennengelernt. Ein Problem dieses Maßes ist, dass bei Hinzunahme eines weiteren Prädiktors der Anteil der aufgeklärten Varianz zwangsläufig größer werden muss. Daher gibt es Anpassungen und andere Maße, die gleichzeitig die Anzahl der Prädiktoren mit verrechnen, um dies zu berücksichtigen.
Das einfachste Maß, welches gleichzeitig eine Schätzung der
aufgeklärten Varianz auf Populationsebene darstellt, ist das
adjustierte \(R^2\), abgekürzt
oft als \(R^2_{adj}\). Die Funktion
summary() angewendet auf ein lm()-Objekt gibt
dieses Maß bereits mit aus und berechnet es nach folgender Formel: \[R^2_{adj}=1-(1-R^2)\cdot\frac{n-1}{n-q-1}\]
wobei \(n\) die Anzahl der
Beobachtungen und \(q\) die Anzahl der
Prädiktoren ist. Es gibt aber eine ganze Reihe von Vorschlägen gibt, wie
\(R^2_{adj}\) berechnet werden sollte
und die Ergebnisse können sich im Detail dann auch unterscheiden können
(für mehr Details, siehe Yin & Fan, 2001).
7.4 Inferenzstatistik
Bezüglich der inferenzstatistischen Auswertung können mehrere Herangehensweisen unterschieden werden, die im Folgenden einzeln beschrieben werden.
7.4.1 allgemeiner Modelltest
Die Zusammenfassung der Ergebnisse eines Regressionsmodells liefert immer auch einen \(F\)-Test mit, der testet, ob das Modell in seiner Gänze einen signifikanten Beitrag zur Varianzaufklärung leistet. In anderen Worten heißt dies, dass hier die Nullhypothese \(H_0: R^2=0\) getestet wird. Woher der \(F\)-Wert stammt, können wir uns veranschaulichen, wenn wir eine Quadratsummenzerlegung durchführen, genau wie wir es im Rahmen der Varianzanalyse durchgeführt haben. Der Ursprung der drei Quadratsummen können wir anhand der folgenden Abbildung am Beispiel einer einfachen linearen Regression sehen:
- Die totale Quadratsumme ergibt sich aus der Summe der quadrierten Abweichungen der Daten vom Mittelwert des Kriteriums: \[SS_{total}=\sum_{i=1}^n (y_i-M_Y)^2\]
- Die Fehler-Quadratsumme ergibt sich aus der Summe der quadrierten Residuen: \[SS_{Fehler}=\sum_{i=1}^n (y_i-\hat{y}_i)^2=\sum_{i=1}^n e_i^2\]
- Die Modell-Quadratsumme ergibt sich aus der Summe der quadrierten Abweichungen der Regressionsgeraden (also der vorhergesagten Werte) vom Mittelwert des Kriteriums: \[SS_{Modell}=\sum_{i=1}^n (\hat{y}_i-M_Y)^2\]
Durch diese Abbildung wird bereits klar, dass die totale Quadratsumme die Summe der Fehler-Quadratsumme und der Modell-Quadratsumme ist: \[SS_{total}=SS_{Modell}+SS_{Fehler}\] Zur weiteren Illustration berechnen wir die Quadratsummen für die Daten, die der letzten Abbildung zugrunde liegen:
X <- c(1:8)
Y <- c(3, 8, 5, 8, 12, 11, 18, 15)
M.Y <- mean(Y) # Mittelwert von Y
modell <- lm(Y ~ X) # Regression von Y auf X
( SS.total <- sum( (Y-M.Y)^2) )## [1] 176( SS.Fehler <- sum((modell$residuals)^2) )## [1] 31.14286( SS.Modell <- sum( (modell$fitted.values - M.Y)^2) )## [1] 144.8571SS.Fehler + SS.Modell # = SS.total## [1] 176Ganz analog zur Vorgehensweise bei der Varianzanalyse (vgl. Teil 2) können wir nun einen \(F\)-Bruch als Quotienten der durch ihre Freiheitsgrade dividierten Quadratsummen (das sind dann Mittlere Quadratsummen) bilden: \[F=\frac{ \frac{SS_{Modell}}{df_{Modell}} }{ \frac{SS_{Fehler}}{df_{Fehler}} }=\frac{MS_{Modell}}{MS_{Fehler}}\] Die Freiheitsgrade ergeben sich aus der Anzahl \(q\) der Prädiktoren und der Anzahl \(n\) der Beobachtungen: \[df_{Modell}=q\quad\text{ und }\quad df_{Fehler}=n-q-1\] Die entsprechende Zufallsvariable ist unter der Nullhypothese \(F\)-verteilt mit \(q\) Zählerfreiheitsgraden und \(n-q-1\) Nennerfreiheitsgraden: \[\mathbf{F}\overset{H_0}{\sim} F_{q;n-q-1}\] In unserem Beispiel ergibt dies für \(q = 1\) und \(n=8\):
q <- 1
n <- 8
F <- (SS.Modell/q) / (SS.Fehler/(n-q-1))
F## [1] 27.90826Alternativ kann der \(F\)-Bruch auch direkt aus \(R^2\) berechnet werden (auf diese Formel kommen wir in ganz ähnlicher Form gleich noch einmal zurück): \[F=\frac{(n-q-1)\cdot R^2}{q\cdot(1-R^2)}\] Betrachten wir nun zum Vergleich die Zusammenfassung des Ergebnisses der Regressionsanalyse mit den Daten der Abbildung, finden wir den \(F\)-Wert und die Freiheitsgrade sowie den entsprechenden \(p\)-Wert ganz am Ende der Ausgabe wieder:
S(modell)## Call: lm(formula = Y ~ X)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.6429 1.7752 0.925 0.39044
## X 1.8571 0.3515 5.283 0.00186 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 2.278 on 6 degrees of freedom
## Multiple R-squared: 0.8231
## F-statistic: 27.91 on 1 and 6 DF, p-value: 0.00186
## AIC BIC
## 39.58 39.81Alternativ können wir diese Ergebnisse mit der Funktion
anova() auch in Form der aus Teil
2 bekannten ANOVA-Tabelle erhalten:
anova(modell)Wenngleich wir hier das Beispiel der einfachen, linearen Regression genutzt haben, ist das generelle Vorgehen auch im Fall der multiplen Regression völlig identisch.
7.4.2 Tests der Regressionskoeffizienten
In Teil 12 von Statistik I haben wir bereits festgestellt, wie die Steigung \(b\) in der einfachen linearen Regression gegen die Nullhypothese \(H_0:\beta=0\) getestet wird. In diesem Fall liefern der Test von \(b\) und der gerade besprochene \(F\)-Test die gleiche Information. Dies kann man auch daran sehen, dass wieder die Beziehung \(F=t^2\) gilt. (Eine Generalisierung dieser Beziehung werden wir im nächsten Abschnitt kennenlernen.)
Im Fall mehrerer Prädiktoren kann zwar der allgemeine Modelltest
signifikant ausfallen, dennoch muss nicht für jeden
Regressionskoeffizienten an sich die Alternativhypothese gelten. Die
Tests der einzelnen Nullhypothesen laufen wiederum über \(t\)-Tests und finden sich sämtlich in den
Ausgaben der summary()- oder S()-Funktion
wieder.
7.4.3 Modellvergleich
Weiter oben haben wir bereits den \(F\)-Test der Nullhypothese \(H_0:R^2=0\) kennengelernt. Sogenannte Modellvergleiche ermöglichen darüber hinaus zu testen, ob sich durch Hinzunahme (oder Wegnahme) von Prädiktoren eine signifikante Veränderung von \(R^2\) einstellt. Dieses Konzept spielt bei vielen komplexeren linearen Modellen eine wichtige Rolle. Oft wird ein sog. reduziertes Modell mit einem sog. komplexen Modell verglichen, wobei wichtig ist, dass das reduzierte Modell “genested” im komplexen Modell ist, d.h. alle seine Prädiktoren auch im komplexen Modell enthalten sind.
Der \(F_{change}\)-Wert berechnet sich in diesem Fall als \[F_{change}=\frac{(n-q-1)\cdot (R_{komplex}^2-R_{reduziert}^2)}{J\cdot(1-R_{komplex}^2)}\] Hierbei meint \(n\) wieder die Anzahl der Beobachtungen, \(q\) ist die Anzahl der Prädiktoren im komplexen Modell und \(J\) ist die Anzahl der Prädiktoren, die sich zwischen komplexem und reduziertem Modell unterscheiden.
Als Beispiel betrachten wir den Datensatz mit den Albumsverkäufen.
Bisher haben wir nur zwei Prädiktoren genutzt; es gibt aber noch eine
weitere Variable Attraktivitaet im Datensatz. Wir wollen
hier nun wissen, ob sich ein signifikanter Zuwachs in \(R^2\) ergibt, wenn zum reduzierten Modell
mit Werbung als Prädiktor beide weiteren Variablen
gleichzeitig als Prädiktoren hinzugegeben werden, also wenn das komplexe
Modell drei Prädiktoren enthält: \[\begin{equation*}
\begin{aligned}
\text{reduziertes Modell: }&\hat{Y}=a + b_1\cdot \text{Werbung}\\
\text{komplexes Modell: }&\hat{Y}=a + b_1\cdot \text{Werbung}+
b_2\cdot \text{Frequenz}+ b_3\cdot \text{Attraktivitaet}\\
\end{aligned}
\end{equation*}\]
Die (beiden) Modelle die verglichen werden sollen müssen dann
zunächst mit lm() berechnet werden:
daten <- read.table("./Daten/Album_Verkaeufe.dat",
header = TRUE)
modell.1 <- lm(Verkaeufe ~ Werbung,
data = daten)
modell.3 <- lm(Verkaeufe ~ Werbung + Frequenz + Attraktivitaet,
data = daten)Ohne uns die Ergebnisse nun im einzelnen anzuschauen, berechnen wir
direkt den \(F\)-Test, indem wir die
Funktion anova() auf die beiden Modelle anwenden, wobei die
Modelle in aufsteigender Komplexität übergeben werden:
anova(modell.1, modell.3)Es ergibt sich mit \(F_{change}=96.45\) also ein signifikanter Zuwachs in \(R^2\) durch die Hinzunahme der beiden Prädiktoren.
Unterscheiden sich die beiden Modelle nur um einen Prädiktor, entspricht dieser \(F\)-Test dem \(t\)-Test des entsprechenden Prädiktors gegen die Nullhypothese \(H_0:\beta_p=0\) und es gilt wiederum \(F=t^2\).
7.5 Polynomiale Regression
Im Fall nur eines Prädiktors haben wir bisher immer eine Regressionsgeraden \[ \hat{Y}=bX+a \] berechnet. Die multiple Regression erlaubt uns aber noch eine ganz flexible Erweiterung, um nicht-lineare Zusammenhänge zu erfassen, die sich in Daten finden lassen. Als Beispiel schauen wir einen Fall des sog. Anscombe-Quartetts an. Ganz offenkundig herrscht hier ein kurvilinearer Zusammenhang. Dennoch führen wir zunächst eine einfache lineare Regression durch und schauen uns das Ergebnis an (linker Teil der folgenden Abbildung). Immerhin kann diese Regression bereits etwa 67% der Varianz von \(Y\) aufklären; aber eine Gerade scheint in der Tat nicht die beste Beschreibung der Daten zu sein, was insbesondere auch klar wird, wenn man sich die Residuen ansieht (rechter Teil der folgenden Abbildung).
par(mfrow = c(1,2))
# Daten des Beispiels
X <- c(4:14)
Y <- c(3.1, 4.74, 6.13, 7.26, 8.14, 8.77, 9.14, 9.26, 9.13, 8.74, 8.1)
m1 <- lm(Y ~ X) # Regression vo Y auf X und speichern
# Plot linker Teil
plot(X,
Y,
pch = 19,
cex = 0.5,
xlim = c(0, 20),
ylim =c (0, 15),
ylab = "Y")
abline(m1, lty = 2) # Regressionsgerade einzeichnen
# Plot rechter Teil
plot(X,
m1$residuals, # Residuen der Regression von Y auf X
pch = 19,
cex = 0.5,
xlim = c(0, 20),
ylim = c(-5, 5),
ylab = "E")
S(m1) # Ergebnis der Regression## Call: lm(formula = Y ~ X)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.001 1.125 2.667 0.02576 *
## X 0.500 0.118 4.239 0.00218 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 1.237 on 9 degrees of freedom
## Multiple R-squared: 0.6662
## F-statistic: 17.97 on 1 and 9 DF, p-value: 0.002179
## AIC BIC
## 39.69 40.89Zur Beschreibung solcher Daten bietet sich eine quadratische
Funktion, eine sog. Parabel, an. Oft werden Parabeln in ihrer sog.
Scheitelform beschrieben, z.B. als \[
f(x)=(x-5)^2+10\ .
\] Scheitelform heißt diese Darstellung deshalb, weil der
Scheitelpunkt hier direkt abgelesen werden kann. Im Beispiel liegt er
bei den Koordinaten \((5, 10)\). Die
allgemeine Form dieser Gleichung ist allerdings ein Polynom 2. Grades:
\[\hat{Y}=a + b_1\cdot X + b_2 \cdot X^2\
.\] In dieser Form geschrieben, liegt die Gleichung einer
multiplen Regression vor, wobei der erste Prädiktor \(X\) ist und der zweite Prädiktor die
quadrierten Werte von \(X\), also \(X^2\). Da derartige Funktionen Polynome
heißen, wird diese Regression auch polynomiale Regression
genannt. Im Prinzip können auch Polynome höheren Grades benutzt werden,
sodass die allgemeine Form lautet: \[\hat{Y}=a + b_1\cdot X + b_2 \cdot X^2 + \ldots +
b_q\cdot X^q=a+\sum_{p=1}^q b_p X^p\ .\] Die Umsetzung in R
erfolgt wie bei einer multiplen Regression. In unserem Fall wollen wir
neben \(X\) auch \(X^2\) als Prädiktor in das Modell
aufnehmen. Dazu wird X^2 in die Klammer von
I() geschrieben, um deutlich zu machen, dass hier eine
arithmetische Operation gemeint ist und das ^-Zeichen nicht
im Sinne der Modellsprache interpretiert werden soll:
m2 <- lm(Y ~ X + I(X^2))
S(m2)## Call: lm(formula = Y ~ X + I(X^2))
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.9957343 0.0043299 -1385 <2e-16 ***
## X 2.7808392 0.0010401 2674 <2e-16 ***
## I(X^2) -0.1267133 0.0000571 -2219 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 0.001672 on 8 degrees of freedom
## Multiple R-squared: 1
## F-statistic: 7.378e+06 on 2 and 8 DF, p-value: < 2.2e-16
## AIC BIC
## -104.94 -103.35Da \(R^2=1\) ist, scheint das Modell die Daten also komplett zu erfassen und beschreiben zu können. Als Regressionskoeffizienten haben sich ergeben:
- \(a=-5.9957343\)
- \(b_1=2.7808392\)
- \(b_2=-0.1267133\)
Damit können wir nun die entsprechende Funktion zeichnen, die die Datenpunkte bestmöglich beschreibt:
plot(X, Y, # Scatterplot der Datenpunkte
pch = 19,
cex = 0.5,
xlim = c(0, 20),
ylim = c(0, 15),
ylab = "Y")
x.neu <- c(3:15) # x-Werte für die die Funktion gezeichnet werden soll
# dann die neuen y-Werte entsprechend mit den Koeffizienten des Polynoms berechnen...
y.neu <- -5.9957343 + 2.7808392 * x.neu - 0.1267133 * x.neu^2
lines(x.neu, # ...und als gestrichelte Linie einzeichnen
y.neu,
lty = 2)
Eine andere Variante der Formulierung polynomialer Regressionen bietet
die Funktion poly(), die insbesondere bei Polynomen höherer
Ordnung nützlich ist, da die höchste Potenz dem Argument
degree übergeben wird:
m2.poly <- lm(Y ~ poly(X, degree = 2, raw = TRUE))
S(m2.poly) # gleiches Ergebnis wie oben## Call: lm(formula = Y ~ poly(X, degree = 2, raw = TRUE))
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.9957343 0.0043299 -1385 <2e-16 ***
## poly(X, degree = 2, raw = TRUE)1 2.7808392 0.0010401 2674 <2e-16 ***
## poly(X, degree = 2, raw = TRUE)2 -0.1267133 0.0000571 -2219 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 0.001672 on 8 degrees of freedom
## Multiple R-squared: 1
## F-statistic: 7.378e+06 on 2 and 8 DF, p-value: < 2.2e-16
## AIC BIC
## -104.94 -103.35Zum Abschluss betrachten wir noch Daten, die mit einem Polynom
dritten Grades erzeugt wurden und nutzen Modellvergleiche. Hieran können
wir praktisch sehen, wie sukzessive Prädiktoren einem Modell hinzugefügt
werden können und getestet wird, ob das Modell bzw. dessen
Vorhersagekraft dadurch besser geworden ist. Dazu generieren wir
zunächst zufällig zwei Vektoren der “empirischen Daten”. Der Vektor
X beinhaltet die \(x\)-Werte und wird aus einer
Gleichverteilung gezogen; die Werte von Y resultieren
daraus, dass zunächst Werte mit einem entsprechenden Polynom berechnet
werden und anschließend “nach oben oder unten” versetzt werden, indem
ein Wert aus einer Normalverteilung addiert wird:
set.seed(1232) # zur Reproduzierbarkeit (ändern Sie den Wert ruhig!)
a <- 10 # Koeffizienten vorgeben zur Datengenerierung
b1 <- 1
b2 <- 0.5
b3 <- 0.3
X <- runif(n = 50, # 50 x-Werte aus Gleichverteilung ziehen
min = -10,
max = 10)
Y <- a + b1 * X + b2 * X^2 + b3 * X^3 # Berechnung der y-Werte aus Polynom...
Y <- Y + rnorm(length(X),0,20) # ...und zufällige VerschiebungDie folgenden Abbildungen visualisieren zunächst die Daten und die drei Modelle, die im Anschluss dann erstellt und evaluiert werden. Optisch ist vermutlich direkt klar, dass das kubische Modell \[ \hat{Y}=a + b_1\cdot X + b_2 \cdot X^2 + b_3 \cdot X^3 \] am besten zu den Daten passt.
Diesen Eindruck testen wir nun mit entsprechenden polynomialen Regressionen und den jeweiligen Modellvergleichen. Die drei Modelle, die wir miteinander vergleichen wollen, sind die folgenden:
modell.1 <- lm(Y ~ X) # linear
modell.2 <- lm(Y ~ X + I(X^2)) # quadratisch
modell.3 <- lm(Y ~ X + I(X^2) + I(X^3)) # kubischAls erstes lassen wir die jeweiligen \(R^2\) Werte ausgeben, also den Anteil an Kriteriumsvarianz, den ein Modell aufklären kann:
summary(modell.1)$r.squared
## [1] 0.8689258
summary(modell.2)$r.squared
## [1] 0.8754764
summary(modell.3)$r.squared
## [1] 0.9850817Dass die Modelle die Daten jeweils “besser” widerspiegeln hat sich ja
bereits optisch gezeigt, aber es zeigt sich auch in der ansteigenden
aufgeklärten Varianz, die vom quadratischen zum kubischen Modell noch
einmal einen großen Sprung macht. Nun wollen wir wissen, ob die
jeweiligen Zuwächse in \(R^2\) auch
statistisch signifikant sind. Dazu führen wir zwei Modellgleiche vor,
nämlich die Hinzunahme des quadratischen Anteils zum bloßen linearen
Modell und die Hinzunahme des kubischen Anteils zum quadratischen
Modell. Diese beiden Modellvergleiche können wir direkt aufrufen, indem
wir alle drei Modelle in aufsteigender Komplexität an die Funktion
anova() übergeben:
anova(modell.1, modell.2, modell.3)Es kommt also jeweils durch die Hinzunahme eines weiteren Polynoms als Prädiktor zu einer signifikant höheren Varianzaufklärung. Prinzipiell könnten wir natürlich noch das Polynom 4. Ordnung dem Modell hinzufügen:
modell.4 <- lm(Y ~ X + I(X^2) + I(X^3) + I(X^4))
summary(modell.4)$r.squared## [1] 0.985095Tatsächlich ist das nun erreichte \(R^2\) noch einmal minimal gestiegen. Der entsprechende Modellvergleich zeigt allerdings an, dass diese Zunahme statistisch nicht signifikant ist. Aus Gründen der Sparsamkeit sollten wir daher weiter von einem kubischen Modell, d.h. einem kubischen Zusammenhang bei der Variblen ausgehen.
anova(modell.3, modell.4)7.6 Literatur
Baguley, T. (2012). Serious stat: A guide to advanced statistics for the behavioral sciences. Bloomsbury Publishing.
Fox,J. & Weisberg, S. (2019). An R companion to applied regression. Sage.
Field, A., Miles, J. & Field, Z. (2012). Discovering statistics using R. Sage Publications Ltd.
Yin, P. & Fan, X. (2001). Estimating \(R^2\) shrinkage in multiple regression: A comparison of different analytical methods. The Journal of Experimental Education, 69, 203-224.
8 Multiple Regression: Vertiefung
In Teil 7 haben wir die Regression auf den Fall mehr als eines Prädiktors erweitert. Hier beginnen wir nun mit der Erweiterung um Interaktionen von (intervallskaierten) Prädiktoren. Dann zeigen wir, wie auch nominalskalierte Prädiktoren von einer Regression gehandhabt werden können und, dass hierbei eine Äquivalenz zur Varianzanalyse aus Teil 2 vorliegt. Den Abschluss bilden dann besondere Situationen, die bei den Prädiktoren auftreten und potenziell Probleme bereiten können.
# Pakete die in diesem Teil benutzt werden:
library(car)
library(effects)
library(emmeans)
library(schoRsch)
library(interactions)8.1 Interaktionen zweier intervallskalierter Variablen
8.1.1 Hinzufügen eines Interaktionsterms
Das Konzept der Interaktion haben wir im Rahmen der Varianzanalyse bereits ausgiebig kennengelernt: Eine Interaktion meint, dass die Wirkung einer Variablen von der betrachteten Stufe einer anderen Variablen abhängt. Bisher waren aber beide Variablen nominalskaliert (eben die Faktoren bei der Varianzanalyse). Interaktionen können aber auch generalisiert betrachtet werden, wenn zwei (oder mehr) Variablen (mindestens) intervallskaliert sind.
Verändert eine Variable die Wirkung einer anderen Variablen, so wird
diese oft als Moderatorvariable bezeichnet und die
entsprechende Analyse dann als moderierte (multiple)
Regression. Technisch gesehen ist dies aber nichts anderes als
eine multiple Regression, die einen Interaktionsterm berücksichtigt. Wir
betrachten hier daher den Fall zweier Prädiktoren: \[\hat{Y} = a+b_1\cdot X_1+b_2\cdot X_2\] In
varianzanalytischer Terminologie entsprechen die Effekte von \(X_1\) und \(X_2\) den Haupteffekten, die eine additive
Wirkung auf das Kriterium haben. Der Trick der Interaktion ist nun, eine
Variable \(X_3\) additiv hinzuzufügen,
die aber ihrerseits nicht-additive Effekte repräsentiert. Dies gelingt
ganz einfach, indem die Interaktion von \(X_1\) und \(X_2\) als deren Produkt aufgefasst wird,
also: \[\hat{Y} = a+b_1\cdot X_1+b_2\cdot
X_2+b_3\cdot X_3\quad\text{mit } X_3=X_1\cdot X_2\ .\] Warum
bewirkt diese simple Multiplikation genau das, was eine Interaktion
verbal beschrieben darstellt? Um dies zu sehen, setzen wir \(X_3\) in die Gleichung ein und klammern
dann \(X_2\) aus. Nun sehen wir, dass
die Gleichung einer einfachen linearen Regression resultiert, deren
Achsenabschnitt \(a'\) und Steigung
\(b'\) aber vom Wert auf der
Variablen \(X_1\) abhängen. In anderen
Worten, für Personen mit jeweils anderen Werten auf der Variablen \(X_1\) verändert sich die (lineare)
Beziehung von \(Y\) und \(X_2\) (natürlich kann man auch \(X_1\) ausklammern und die Interaktion
anders herum betrachten): \[\begin{equation*}
\begin{aligned}
\hat{Y} &= a+b_1\cdot X_1+b_2\cdot X_2+b_3\cdot X_1\cdot X_2\\
&= \underset{a'}{\underbrace{a+b_1\cdot
X_1}}+ \underset{b'}{\underbrace{(b_2+b_3\cdot X_1)}}\cdot X_2
\end{aligned}
\end{equation*}\] Das Schöne daran ist, dass \(X_3\) das Zusammenspiel von \(X_1\) und \(X_2\) abbildet, aber additiv hinzugefügt
wird. Dadurch ändert sich eigentlich nichts, außer dass es eben eine
dritte Prädiktorvariable gibt. Tatsächlich kann diese Variable einfach
berechnet werden, indem für jede Person die Werte auf den beiden
Ausgangsvariablen multipliziert werden. Nötig ist dies aber nicht, da R
die Interaktion selbst hinzufügt, wenn diese angefordert wird. Dies
geschieht entweder mit Y ~ X1 + X2 + X1:X2 oder in kurz
Y ~ X1 * X2. Einzig sollte darauf geachtet werden, dass
eine Interaktion nur dann Sinn ergibt, wenn auch die Ausgangsvariablen
im Modell enthalten sind. Interaktionen höherer Ordnung werden ganz
analog formuliert.
8.1.2 Datenbeispiel und Durchführung
Um besser darstellen zu können, wie sich eine Interaktion auswirkt,
nutzen wir den Datensatz regression_ia.dat:
daten_regression_ia <- read.table(file = "./Daten/regression_ia.dat",
header = TRUE)
some(daten_regression_ia) # aus car: zeigt zufällige Datenpunkte anIn dieser Datei befinden sich drei Variablen für 100 Personen. Die
Variable Training beschreibt eine Einschätzung, wieviel
Prozent ihrer Zeit eine Person auf ein Sporttraining aufwendet (0-100%).
Auf der Variable Talent hat jede Person auf einer Skala von
0-20 ihr eigenes Talent eingeschätzt. Die Variable Leistung
dokumentiert auf einer Skala von 1-25 die erbrachte Leistung (höhere
Werte sprechen für besser Leistung).
Zunächst ist es naheliegend, dass sowohl Training als
auch Talent positiv mit Leistung
zusammenhängen:
cor_out(cor.test(daten_regression_ia$Training, daten_regression_ia$Leistung))## Text
## 1 r(100) = .49, p < .001cor_out(cor.test(daten_regression_ia$Talent, daten_regression_ia$Leistung))## Text
## 1 r(100) = .75, p < .001Man würde also positive Werte für \(b_1\) und \(b_2\) in einer Regression mit beiden Prädiktoren erwarten:
lm.result <- lm(Leistung ~ Training + Talent,
data = daten_regression_ia)
S(lm.result, brief = TRUE)## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.658602 0.537084 -8.674 9.73e-14 ***
## Training 0.117861 0.006434 18.320 < 2e-16 ***
## Talent 0.804103 0.031444 25.573 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 1.879 on 97 degrees of freedom
## Multiple R-squared: 0.9023
## F-statistic: 447.9 on 2 and 97 DF, p-value: < 2.2e-16
## AIC BIC
## 414.90 425.32Tatsächlich ist es hier aber auch naheliegend, dass beide Variablen interagieren. Training bewirkt vielleicht eine umso bessere Leistung, je mehr Talent eine Person hat. Dazu fügen wir nun einfach die Interaktion mit hinzu:
lm.result.ia <- lm(Leistung ~ Training * Talent,
data = daten_regression_ia)
S(lm.result.ia, brief = TRUE)## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.893141 0.125760 7.102 2.14e-10 ***
## Training 0.010930 0.002044 5.347 6.03e-07 ***
## Talent 0.300874 0.009716 30.968 < 2e-16 ***
## Training:Talent 0.010057 0.000166 60.605 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 0.3015 on 96 degrees of freedom
## Multiple R-squared: 0.9975
## F-statistic: 1.283e+04 on 3 and 96 DF, p-value: < 2.2e-16
## AIC BIC
## 49.88 62.90Tatsächlich ist auch der Beitrag der Interaktion signifikant.
Interpretationen solcher Interaktionen sind aber für sich genommen oft
nicht ganz einfach. Sehr hilfreich ist es hier, sich die Ergebnisse zu
visualisieren, indem man das Kriterium gegen einen Prädiktor plottet und
dies getrennt für ausgewählte Werte des anderen Prädiktors tut. Welche
Werte man hier wählt, hängt auch von inhaltlichen Gesichtspunkten ab;
häufig wird aber z.B. der Mittelwert, Mittelwert minus eine
Standardabweichung und der Mittelwert plus eine Standardabweichung dafür
genommen. Der folgende R-Code nutzt dafür die Funktion
interact_plot():
interact_plot(lm.result.ia, # lm()-Objekt
pred = Training, # x-Achse
modx = Talent, # Moderator mit separaten Geraden
plot.points = TRUE, # Datenpunkte hinzufügen
x.label = "Training", # Bezeichnung x-Achse
y.label = "Leistung") # Bezeichnung y-Achse
In dieser Abbildung kann man sehr schön sehen, dass in Abhängigkeit
von der Ausprägung auf der Variablen Talent sich die
lineare Beziehung von Leistung und Training
verändert.
8.2 Nominalskalierte Prädiktoren
Nominalskalierte (oder kategoriale) Variablen haben Kategorien als Werte, die keine Ordnung beinhalten. Typische Beispiele sind z.B. Geschlecht (weiblich, männlich, divers), Haarfarbe (schwarz, braun, blond, …), Schulabschluss (Mittlere Reife, Abitur, …) und viele mehr. Auch Gruppenzugehörigkeit im Rahmen eines Experimentes kann dazu gehören, z.B. mit drei Stufen wie Experimentalgruppe 1, Experimentalgruppe 2 und Kontrollgruppe. Solche Variablen werden in anderen Kontexten, insbesondere der Varianzanalyse, auch Faktoren genannt.
Spalten mit Zeichenketten (Strings) als Werte werden von R-Befehlen
i.d.R automatisch als Faktoren behandelt. Andere Variablen müssen
zunächst mit as.factor() (oder factor) in
einen Faktor umgewandelt werden. Tatsächlich stellen derartige Variablen
als Prädiktor aber kein Problem für eine Regression dar. Der Trick ist
dabei, dass eine nominalskalierte Variable durch einen Satz anderer
Variablen kodiert wird, mit denen eine Regression wieder umgehen kann.
Allgemein gilt, dass ein Faktor mit \(J\) Stufen durch \(J-1\) neue Variablen kodiert wird, die dann
(anstelle des eigentlichen Faktors) in die Regression eingehen.
8.2.1 Binärer Prädiktor vs. \(t\)-Test
Wir beginnen hier mit dem Spezialfall einer binären (dichotomen) Variablen, also einer, die nur zwei Werte annehmen kann. Dies könnten z.B. Experimental- versus Kontrollgruppe sein und zum Vergleich der abhängigen Variablen dieser beiden Gruppen würde sich ein \(t\)-Test für zwei unabhängige Stichproben anbieten (siehe Teil 12 in Statistik 1). Bei \(J=2\) Stufen gestaltet sich die Kodierung der beiden Stufen noch sehr einfach da wir nur \(J-1=2-1=1\) Variable dafür benötigen. Nun gibt es zwei Möglichkeiten der Kodierung, die als Dummykodierung und als Effektkodierung bezeichnet werden:
- Dummykodierung: 0 vs. 1
- Effektkodierung: -1 vs. 1
Für den entspechenden \(t\)-Test hat dies keinerlei Auswirkungen. Wir wollen aber hier schauen, wie sich eine Regression mit dem binären Prädiktor Gruppenzugehörigkeit verhält. Daher betrachten wir nun ein Beispiel und erstellen einen Datensatz mit 20 Beobachtungen auf einer Variablen \(Y\). Den Prädiktor kodieren wir dann einmal in Dummyform und einmal in Effektform:
daten <- data.frame(X.dummy = rep(c(0, 1), # dummy-kodiert
each = 10),
X.effekt = rep(c(-1, 1), # effekt-kodiert
each = 10),
Y = c(8, 9, 10, 8, 9, 7, 4, 6, 4, 8,
10, 12, 13, 14, 15, 14, 13, 15, 20, 18))
brief(daten) # zeigt die ersten 3 und die letzten 2 Datenzeilen an## 20 x 3 data.frame (15 rows omitted)
## X.dummy X.effekt Y
## [n] [n] [n]
## 1 0 -1 8
## 2 0 -1 9
## 3 0 -1 10
## . . .
## 19 1 1 20
## 20 1 1 18Mit diesen Daten berechnen wir nun einen \(t\)-Test, indem wir die dummy-kodierte Gruppenvariable benutzen:
t_result <- t.test(Y ~ X.dummy,
data = daten,
var.equal = TRUE)
t_result##
## Two Sample t-test
##
## data: Y by X.dummy
## t = -6.3504, df = 18, p-value = 5.543e-06
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## -9.448902 -4.751098
## sample estimates:
## mean in group 0 mean in group 1
## 7.3 14.4Die Mittelwerte der Gruppen sind also \(7.3\) und \(14.4\). Der Grand Mean ist \(\frac{7.3+14.2}{2}=10.85\) und der Unterschied zwischen beiden Mittelwerten ist \(7.1\). Der Betrag des \(t\)-Werts schließlich ist \(|t|=6.3504\).
Nun formulieren wir die Analyse als Regression und verwenden
X.dummy als (dummy-kodierten) Prädiktor. Da uns vor allem
die Koeffizienten interessieren, lassen wir uns auch diese zunächst
gesondert ausgeben:
dummy_result <- lm(Y ~ X.dummy,
data = daten)
coef(dummy_result)## (Intercept) X.dummy
## 7.3 7.1Hier sehen wir mehrere Beziehungen zu den Mittelwerten der Gruppen:
- Der Achsenabschnitt entspricht dem Mittelwert der mit 0 kodierten (Referenz-)Gruppe.
- Die Steigung entspricht der Differenz beider Mittelwerte.
Auch der \(t\)-Wert des Signifikanztests der Steigung entspricht dem, den auch der \(t\)-Test bereits geliefert hat:
summary(dummy_result)##
## Call:
## lm(formula = Y ~ X.dummy, data = daten)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.40 -1.40 0.15 0.95 5.60
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.3000 0.7906 9.234 3.00e-08 ***
## X.dummy 7.1000 1.1180 6.350 5.54e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.5 on 18 degrees of freedom
## Multiple R-squared: 0.6914, Adjusted R-squared: 0.6743
## F-statistic: 40.33 on 1 and 18 DF, p-value: 5.543e-06Diese Regression und der \(t\)-Test
scheinen also die gleiche Information zu beinhalten. Die Interpretation
der Koeffizienten ändert sich ein wenig, wenn wir die Regression mit
X.effekt als (effekt-kodiertem) Prädiktor berechnen:
effekt_result <- lm(Y ~ X.effekt,
data = daten)
coef(effekt_result)## (Intercept) X.effekt
## 10.85 3.55In diesem Fall gilt:
- Der Achsenabschnitt entspricht dem Grand Mean.
- Die Steigung entspricht nun dem Unterschied vom Grand Mean zu den beiden Mittelwerten (also die halbe Differenz beider Mittelwerte). Im Rahmen der Varianzanalyse haben wir diese Abweichungen Effekte genannt.
- Der \(t\)-Wert des Signifikanztests der Steigung ist wieder identisch zum oben durchgeführten \(t\)-Test.
Insgesamt liefern beide Auswertungsstrategien also i.W. die gleichen Informationen und scheinen ineinander überführbar zu sein. Tatsächlich liefert das sog. Allgemeine Lineare Modell (ALM) (siehe Teil 12 in Statistik II) ein übergeordnetes Dach für Regressionen und Varianzanalysen; beide Verfahren sind also nur Spezialfälle des ALM. Nun werden wir die Form des Umgangs mit nominalskalierten Variablen in einer Regression auf den Fall von mehr als zwei Stufen erweitern.
8.2.2 Kodierung bei mehr als zwei Stufen
Wir benutzen hier im Folgenden das Beispiel aus den Teilen
2-3 mit der Frage, ob die Leistung in einem Sprachtest mit der
Anzahl der Trainingstage zunimmt. Wir laden daher den Datensatz und
stellen zunächst fest, dass die Variable Trainingstage
nicht von vornherein ein Faktor ist:
daten <- read.table("./Daten/daten_training.DAT",
header = TRUE)
some(daten) # aus car: zeigt zufällige Datenpunkte anclass(daten$Trainingstage)## [1] "integer"Es handelt sich vielmehr um eine Variable vom Typ
integer, also eine Zahl, und als solche würde sie von R
auch interpretiert werden. Daher müssen wir sie, wenn wir sie als
nominalskalierten Faktor auffassen wollen, zunächst in einen Faktor
umwandeln. Die einfachste Art ist es, den Faktor quasi sich selbst
zuzuweisen und die Umwandlung per factor() (oder
as.factor()) durchzuführen:
daten$Trainingstage <- factor(daten$Trainingstage)
class(daten$Trainingstage)## [1] "factor"levels(daten$Trainingstage) ## [1] "1" "2" "3"Die letzte Zeile zeigt die (hier: drei) Stufen des Faktors an.
Typischerweise werden die Stufen alphanumerisch geordnet. Möchte man die
Faktorstufen anders benennen, so kann man das labels
Argument der Funktion factor nutzen:
daten$Trainingstage <- factor(daten$Trainingstage,
labels = c("1 Tag", "2 Tage", "3 Tage"))
levels(daten$Trainingstage) ## [1] "1 Tag" "2 Tage" "3 Tage"Hierbei wird allerdings auch die Ausprägung der Variable geändert, dem sollte man sich bewusst sein
some(daten) # vgl. hierzu den Datensatz direkt nach dem LadenDerartige Faktoren können nun ohne Weiteres in eine Modellgleichung aufgenommen werden, wobei die \(J\) Stufen dann automatisch durch \(J-1\) neue Variablen kodiert werden (die manchmal Regressoren genannt werden). Weiter oben haben wir die üblichen Dummy- und Effektkodierungen bereits Ende angesprochen, als wir den \(t\)-Test mit einer Regression mit einem binären Prädiktor verglichen haben.
Die Standardeinstellung in R ist die Dummykodierung, wobei diese hier “Treatment”-Kodierung genannt wird. Das Prinzip ist eigentlich recht einfach und kann am Beispiel des Faktors Trainingstage demonstriert werden, indem wir die Kodierung anzeigen lassen:
contrasts(daten$Trainingstage)## 2 Tage 3 Tage
## 1 Tag 0 0
## 2 Tage 1 0
## 3 Tage 0 1Jede der drei Zeilen steht hierbei für eine Stufe des Faktors. Die beiden Spalten zeigen die Werte, welche die Dummyvariablen (also die Rgressoren) für jede Stufe annehmen:
- Die erste Stufe
1 Tagwird als Referenz bzw. Baseline aufgefasst und wenndaten$Trainingstagediesen Wert annimmt, erhalten beide Dummyvariablen den Wert0erhalten. - Die zweite Stufe
2 Tagewird so kodiert, dass die Dummyvariable2 Tageden Wert1erhält,3 Tagehingegen den Wert0. - Die dritte Stufe
3 Tagewird so kodiert, dass die Dummyvariable3 Tageden Wert1erhält,2 Tagehingegen den Wert0.
Möchte man eine andere Gruppe als Referenzgruppe kodieren, so muss
man die Reihenfolge der in R kodierten Faktorstufen ändern. Aktuell ist
die Stufe 1 Tag die erste Ausprägung in der Auflistung der
Faktorstufen.
levels(daten$Trainingstage)## [1] "1 Tag" "2 Tage" "3 Tage"Die Gruppe 1 Tag wird somit bei Aufruf von
contrasts() als Referenzgruppe deklariert. Man hat nun 2
Möglichkeiten dies zu ändern. Entweder ändert man die Reihenfolge der
Faktorstufen mit Hilfe des levels Argument der Funktion
factor (ACHTUNG: levels, nicht
labels! Ansonsten werden die Ausprägungen geändert!)
daten$Trainingstage <- factor(daten$Trainingstage,
levels = c("2 Tage", "1 Tag", "3 Tage"))
levels(daten$Trainingstage)## [1] "2 Tage" "1 Tag" "3 Tage"contrasts(daten$Trainingstage) # 2 Tage nun als Referenzgruppe kodiert## 1 Tag 3 Tage
## 2 Tage 0 0
## 1 Tag 1 0
## 3 Tage 0 1Oder man nutzt die Funktion relevel() und ihr Argument
ref um die Referenzgruppe explizit zu benennen
daten$Trainingstage <- relevel(daten$Trainingstage,
ref = "3 Tage")
levels(daten$Trainingstage)## [1] "3 Tage" "2 Tage" "1 Tag"contrasts(daten$Trainingstage)## 2 Tage 1 Tag
## 3 Tage 0 0
## 2 Tage 1 0
## 1 Tag 0 1Wir ändern nun die Ordnung der Faktorstufen wieder auf die ursprüngliche Form für den nächsten Abschnitt.
daten$Trainingstage <- factor(daten$Trainingstage,
levels = c("1 Tag", "2 Tage", "3 Tage"))
levels(daten$Trainingstage)## [1] "1 Tag" "2 Tage" "3 Tage"8.2.3 Einfaktorielle Varianzanalyse als Regression
Mit dieser Vorarbeit können wir nun eine Regression mit dem Kriterium
Sprachtest und dem Prädiktor Trainingstage
durchführen. Mathematisch wird dieser Prädiktor durch die beiden
dummykodierten Variablen repräsentiert. Wir können dann auch hier den
Effekt des Prädiktors mit den gleichen Methoden visualisieren, die wir
bereits kennengelernt haben (auch wenn wir bereits andere Varianten zur
Genüge kennen):
modell <- lm(Sprachtest ~ Trainingstage, # wie eine Regression
data = daten)
plot(predictorEffects(modell))
Nun schauen wir uns das Ergebnis der Regressionsanalyse an, zur Erinnerung lassen wir uns vorher noch einmal die Mittelwerte ausgeben:
Tapply(Sprachtest ~ Trainingstage,
fun = mean,
data = daten)## 1 Tag 2 Tage 3 Tage
## 13.1 15.7 18.1S(modell, brief = TRUE)## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.1000 0.7198 18.199 < 2e-16 ***
## Trainingstage2 Tage 2.6000 1.0180 2.554 0.0166 *
## Trainingstage3 Tage 5.0000 1.0180 4.912 3.86e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 2.276 on 27 degrees of freedom
## Multiple R-squared: 0.472
## F-statistic: 12.07 on 2 and 27 DF, p-value: 0.0001801
## AIC BIC
## 139.33 144.93Ein besonderer Fokus gilt hier dann den Regressionskoeffizienten und ihrer inhaltlichen Interpretation:
- Der Achsenabschnitt
(Intercept)ist der Mittelwert der Referenzgruppe (“1 Trainingstag”). - Der Koeffizient
Trainingstage2 Tageist der Wert, um den der Mittelwert der Referenzgruppe verändert werden muss, um zum Mittelwert der Gruppe “2 Trainingstage” zu gelangen, also \(13.1+2.6=15.7\). - Der Koeffizient
Trainingstage3 Tageist der Wert, um den der Mittelwert der Referenzgruppe verändert werden muss, um zum Mittelwert der Gruppe “3 Trainingstage” zu gelangen, also \(13.1+5.0=18.1\). - Die Signifikanztests testen hierbei die jeweiligen Mittelwertdifferenzen (2 Tage vs. 1 Tag bzw. 3 Tage vs. 1 Tag) unter Berücksichtigung des gesamten ANOVA-Modells (siehe auch den Abschnitt zu Kontrasten in Teil 3).
Eine ähnliche – aber zur Interpretation der Ausgabe anderer möglicher Kodierungen essentielle – Sichtweise ergibt sich, wenn wir noch einmal die Kontrastmatrix betrachten:
contrasts(daten$Trainingstage)## 2 Tage 3 Tage
## 1 Tag 0 0
## 2 Tage 1 0
## 3 Tage 0 1Nun können wir mit den Regressionskoeffizienten und den
Kontrastkoeffizienten zeilenweise die Mittelwerte der Gruppen berechnen.
Bezeichnen wir mit \(b_1\) und \(b_2\) die Regressionskoeffizienten für die
beiden Variablen 2 Tage bzw. 3 Tage und mit
\(a\) den Achsenabschnitt, dann ergibt
sich: \[\begin{equation*}
\begin{aligned}
M_{\text{1 Tag}} &= a + (0)\cdot b_1 + (0)\cdot b_2 \\
M_{\text{2 Tage}} &= a + (1)\cdot b_1 + (0)\cdot b_2 \\
M_{\text{3 Tage}} &= a + (0)\cdot b_1 + (1)\cdot b_2 \\
\end{aligned}
\end{equation*}\]
Dass sich hier tatsächlich die Gruppenmittelwerte ergeben, kann mit R leicht nachvollzogen werden:
a <- as.numeric(coef(modell)[1]) # as.numeric() entfernt den Namen und hat hier nur
b1 <- as.numeric(coef(modell)[2]) # ästhetische Gründe
b2 <- as.numeric(coef(modell)[3])
# jetzt die Mittelwerte berechnen:
a + (0)*b1 + (0)*b2 # Mittelwert Gruppe 1## [1] 13.1a + (1)*b1 + (0)*b2 # Mittelwert Gruppe 2## [1] 15.7a + (0)*b1 + (1)*b2 # Mittelwert Gruppe 3## [1] 18.1Wie an den beiden \(t\)-Tests in der Zusammenfassung der Ergebnisse ferner zu sehen ist, unterscheiden sich also beide Gruppen “2 Trainingstage” und “3 Trainingstage” signifikant von der Referenzgruppe (“1 Trainingstag”). Diese Art von Auswertung testet also nicht die Nullhypothese der Omnibus-Varianzanalyse, sondern eben gezielt zwei Kontraste, die mit der Effektkodierung formuliert werden.
Allerdings steckt in der Regression genau die Information, die auch
die entsprechende einfaktorielle Varianzanalyse ergeben würde. Um die
Ergebnisse in Form einer “klassischen” Varianzanalyse anzuzeigen, kann
die Funktion anova() auf das lm-Objekt
angewendet werden. Wie wir sehen, sind dies exakt die gleichen
Ergebnisse, die wir in Teil
2 bereits berechnet haben:
anova(modell)Ein Vorteil, eine ANOVA als Regression aufzufassen, ist, dass es
komplexere und verallgemeinerte lineare Modelle gibt, für die aber die
gleichen Grundprinzipien und Möglichkeiten wie hier gelten. Zudem kann
ein lm()-Objekt leicht weiter verwendet werden, z.B. um
beliebige Kontraste zu berechnen. Eine sehr flexible Funktion dazu ist
emmeans() aus dem Paket emmeans. Die
Dummykodierung bisher hat ja i.W. zwei Stufen mit der Referenzstufe
verglichen. Wollen wir bspw. nun alle drei paarweisen Vergleiche testen,
so geht dies z.B. mit:
emmeans(modell, pairwise ~ Trainingstage)$contrast## contrast estimate SE df t.ratio p.value
## 1 Tag - 2 Tage -2.6 1.02 27 -2.554 0.0425
## 1 Tag - 3 Tage -5.0 1.02 27 -4.912 0.0001
## 2 Tage - 3 Tage -2.4 1.02 27 -2.358 0.0647
##
## P value adjustment: tukey method for comparing a family of 3 estimatesDer Zusatz $contrast sorgt dafür, dass nur die
Vergleiche ausgegeben werden. Unter $emmeans würden auch
sogenannte adjustierte Mittelwerte ausgegeben werden, die aber in einer
einfaktoriellen Varianzanalyse den normalen Gruppenmittelwerten
entsprechen (siehe aber später in diesem Teil). Im letzten Teil der
Ausgabe ist die verwendete Art der Anpassung für multiple Vergleiche zu
finden (hier Tukey’s HSD-Methode); emmeans() wählt
automatisch eine adäquate Variante aus.
8.3 (Additive) Modelle mit intervallskalierten und nominalskalierten Variablen
8.3.1 Grundagen
Bisher haben wir Modelle betrachtet, die entweder nur intervallskalierte oder nur nominalskalierte Prädiktoren berücksichtigt hatten. Es ist aber natürlich sehr naheliegend, dass im Rahmen einer Regression auch Prädiktoren verschiedener Skalenniveaus berücksichtigt werden können. Wir beginnen hier mit dem Beispiel, dass beide Prädiktoren additiv verwendet werden, d.h. wir lassen keine Interaktion beider Prädiktoren zu. Diese Konstellation wird manchmal auch als Kovarianzanalyse (engl. Analysis of Covariance; ANCOVA) bezeichnet und die intervallskalierte Variable wird dann Kovariate genannt.
Natürlich lässt sich eine solche Regression auch entsprechend als multiple Regression der Form \[\hat{Y}=a+b_1\cdot X_1+b_2\cdot X_2,\] schreiben, wobei (die Kovariate) \(X_1\) hier eine (mindestens) intervallskalierte Variable ist und (der spätere Faktor) \(X_2\) eine nominalskalierte Variable ist. Dass wir hier vom “späteren Faktor” reden hat den Grund, dass ein Faktor mit \(J\) Stufen ja eigentlich über die \(J-1\) (dummy-)kodierten Variablen in dem Modell berücksichtigt wird.
Rechnerisch haben wir hier also eine Regression vorliegen und die Mathematik ist ganz ähnlich wie bei den vorher behandelten Regression auch. Der manchmal genutzte Name “Kovarianzanalyse” legt nahe, dass hier zwei Dinge zusammengebracht werden: Zum einen verwenden wir einen Faktor (oder natürlich auch mehrere Fakoren) wie in einer normalen Varianzanalye und sind an Haupteffekten und ggf. Interaktionen interessiert; zum anderen verwenden wir eine intervallskalierte Variable im Regressionsmodell. Der Einbezug dieser Kovariaten im Kontext einer eigentlich interessanten Varianzanalyse hat vor allem einen wichtigen Grund, nämlich den einer Fehler-Reduktion innerhalb der Gruppen. In der (einfaktoriellen) Varianzanalyse hatte \(SS_w\) die nicht durch den Faktor bzw. die Faktoren erklärte Varianz quantifiziert. Der \(F\)-Bruch bei der einfaktoriellen ANOVA lautete: \[F = \frac{MS_A}{MS_w} = \frac{\frac{SS_A}{J-1}}{\frac{SS_w}{N-J}}\]
Wenn \(SS_w\) kleiner wird, wird natürlich der \(F\)-Bruch größer (die Power steigt also). Eine Kovariate kann in diesem Sinne genutzt werden um Varianz, die nicht durch den Faktor bzw. die Faktoren erklärt wird, zu binden und dadurch die restliche nicht-erklärte Varianz (\(SS_w\)) zu verringern.
Wir beginnen mit einem Beispiel für diese Situation. Im Anschluss daran, verwenden wir die gleiche inhaltliche Fragestellung, aber haben die Daten anders simuliert, um Grenzen und Probleme in der Interpretation zu verdeutlichen.
8.3.2 Beispiel 1: Faktor und Kovariate unkorreliert
Als Beispiel für die folgenden Ausführungen dient die Datei
ANCOVA_beispiel_1.dat:
daten_ancova <- read.table(file = "./Daten/ANCOVA_beispiel_1.dat",
header = TRUE)
some(daten_ancova) # aus car: zeigt zufällige Datenpunkte anDie Datei enthält die Daten von 40 Versuchspersonen, die entweder
einer Kontroll- oder einer Experimentalgruppe zugeordnet wurden. Dies
ist z.B. ein typischer Fall, wenn ein Medikament (oder eine andere Form
einer Therapie) an der Experimentalgruppe getestet wird, während die
Kontrollgruppe diese Therapie nicht erhält. Die Variable
Test enthält eine Einschätzung der Krankheitssymptome (je
geringer desto besser); auf die Variable Pretest kommen wir
gleich noch zu sprechen. Als erstes berechnen wir die beiden
Gruppenmittelwerte und dann eine Varianzanalyse für Test
mit Gruppenzugehörigkeit als Faktor:
Tapply(Test ~ Gruppe,
fun = mean,
data = daten_ancova)## Experimental Kontrolle
## 47.7 48.1daten_ancova$Gruppe <- as.factor(daten_ancova$Gruppe)
# die nächste Zeile sorgt dafür, dass die Kontrollgruppe als Referenz mit 0 kodiert wird
daten_ancova$Gruppe <- relevel(daten_ancova$Gruppe,
ref = "Kontrolle")
levels(daten_ancova$Gruppe)## [1] "Kontrolle" "Experimental"aov.result <- aov(Test ~ Gruppe,
data = daten_ancova)
summary(aov.result)## Df Sum Sq Mean Sq F value Pr(>F)
## Gruppe 1 1.6 1.6000 2.764 0.105
## Residuals 38 22.0 0.5789Wir sehen, dass der Mittelwert für die Experimentalgruppe (\(M=47.7\)) zwar geringer ausfällt als für
die Kontrollgruppe (\(M=48.1\)), aber
dieser Unterschied erweist sich nicht als statistisch signifikant (\(p=.105\)), das Medikament also eher keine
Wirkung hat. Nun erinnern wir uns, dass bei der Durchführung der Studie
doch auch das Wohlbefinden gemessen wurde, bevor die Personen
in die Gruppen eingeteilt wurden. Diese Werte sind als Variable
Pretest im Datensatz auch vorhanden. Naheliegenderweise
korrelieren die Testergebnisse vor und nach der Durchführung der Studie
und ein Teil der Variabilität der Variable Test wird
bereits aufgeklärt durch die Variabilität auf der Variablen
Pretest. Dieser Aspekt ist aber in der Varianzanalyse nicht
berücksichtigt und wir holen dies nun nach (im Ernstfall haben wir das
selbstverständlich genau so geplant geplant!):
# korrelieren beide Variablen?
cor_out(cor.test(daten_ancova$Test, daten_ancova$Pretest))## Text
## 1 r(40) = .34, p = .034Nun verwenden wir die Kovariate Pretest und den
eigentlich interessierenden Faktor Gruppe als Prädiktoren
in einer Regression:
modell1 <- lm(Test ~ Pretest + Gruppe,
data = daten_ancova)
S(modell1, brief = TRUE)## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.7571 7.3174 3.657 0.00079 ***
## Pretest 0.4286 0.1469 2.917 0.00597 **
## GruppeExperimental -0.5714 0.2276 -2.511 0.01654 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 0.6953 on 37 degrees of freedom
## Multiple R-squared: 0.2421
## F-statistic: 5.911 on 2 and 37 DF, p-value: 0.005922
## AIC BIC
## 89.32 96.08Was wir hier sehen ist vor allem, dass die Hinzunahme der Kovariaten
Pretest Varianz des Kriteriums Test gebunden
hat und nun der Effekt des Faktors Gruppe signifikant
geworden ist. Hier hat also die Berücksichtigung einer Kovariaten
tatsächlich geholfen, die Power einer relevanten Analyse zu erhöhen.
Wichtig ist anzumerken, dass der Output des Modells mit
S() oder summary() Tests über
Einzelkoeffizienten ausgibt (ähnlich wie bei der Berechnung der
einfaktoriellen ANOVA weiter oben; was hier allerdings genau getestet
wird ist anders, das klären wir gleich). Die Anzahl der Koeffizienten
die für einen dummykodierten Faktor ausgegeben wird hängt ab von der
Anzahl der Stufen des Faktors. Besitzen wir im Allgemeinen einen Faktor
mit \(J\) Stufen, so werden \(J-1\) Einzelkoeffizienten für diesen Faktor
ausgegeben und jeder Koeffizient testet nur eine bestimmte Differenz
gegenüber der Referenzgruppe. Oft interessieren wir uns aber dafür, ob
ein Faktor (mit all seinen Stufen) einen Einfluss hat. Hierfür können
wir die Funktion Anova() aus dem Paket car
benutzen.
Anova(modell1)In der Zeile “Gruppe” wird geprüft, ob der Faktor in seiner
Gesamtheit einen signifikanten Einfluss hat oder nicht. Da wir hier eine
kategoriale Variable mit nur 2 Stufen betrachten und damit nur eine
Dummyvariable besitzen, sind der Test über den Einzelkoeffizienten (in
S() bzw. summary()) und der Test des Faktors
in seiner Gesamtheit identisch (ab \(J>2\) gilt dies natürlich nicht
mehr).
Ergänzende Anmerkung: In Teil
5 hatten wir kurz das Thema Quadratsummenzerlegung angesprochen,
welches im Falle der mehrfaktoriellen ANOVA mit ungleicher Zellbesetzung
relevant wird. Tatsächlich ist die Frage danach, wie man die
Quadratsummen eines statistischen Modells zerlegt, genereller Natur und
stellt sich immer dann, wenn die Prädiktoren (egal ob intervallskaliert
oder nominalskaliert) nicht unabhängig sind. Im Zuge eines
kovarianzanalytischen Designs wird häufig die sog. Typ II
Quadratsummenzerlegung empfohlen (welche bei der Funktion
Anova()der Standard ist). Welche Quadratsummenzerlegung
allerdings gewählt wird, muss nach inhaltlichen Gesichtspunkten
entschieden werden, denn statistisch sind alle drei Typen der
Quadratsummenzerlegung legitim. Tatsächlich gibt es unseres Wissens nach
keinen allgemeinen Konsens, wann welche Quadratsummenzerlegung zu wählen
ist, sodass wir dieses Thema hier nicht weiter austreten werden. Es sei
jedoch noch angemerkt, dass die Typ II und Typ III
Quadratsummenzerlegung bei einem Modell ohne Interaktion(en) identisch
sind.
Typ I Quadratsummenzerlegungen können mit Hilfe von anova()
berechnet werden und Typ II und Typ III Quadratsummenzerlegungen mit
Anova():
# Typ I: Varianzaufklärung wird erst soweit es geht Pretest zugeschlagen
# -> Test des Faktors Gruppe nur anhand des Varianzanteils, der über Pretest
# hinaus geht
anova(modell1) # Typ II: Test von Pretest bzw. Gruppe anhand der Reduktion der Varianz-
# aufklärung wenn man die jeweilige Variable in allen höheren
# Modellkomponenten weglässt
Anova(modell1,
type = "II")# Typ III: Test von Pretest bzw. Gruppe anhand der Reduktion der Varianz-
# aufklärung wenn man nur die jeweilige Modellkomponente weglässt
# -> im Falle ohne Interaktion identisch zu Typ II
Anova(modell1,
type = "III")Nun wollen wir klären, was die Koeffizienten der Regression
eigentlich bedeuten, und damit auch, was dieses Modell genau testet.
Hierfür visualisieren wir die Daten zunächst einmal. Eine Möglichkeit
der Darstellung sieht wie folgt aus, wobei wir die Werte auf der
Variablen Pretest mit etwas “Jitter” versehen haben, damit
nicht viele Punkte einfach aufeinander liegen:
# Plotten der empirischen Daten (orange = Kontrollgruppe, blau = Experimentalgruppe)
plot(jitter(Test, 0.3) ~ jitter(Pretest, 0.3),
data = daten_ancova,
col = ifelse(Gruppe == "Kontrolle", "orange", "blue"),
pch = 19,
ylab = "Test",
xlab = "Pretest")
# einzeichnen der Regressionsgeraden für beide Gruppen
coefs <- coef(modell1)
abline(coef = coefs[1:2], col = "orange")
abline(coef = c(coefs[1] + coefs[3], coefs[2]), col = "blue")
legend(x = "bottomleft",
legend = c("Kontrollgruppe","Experimentalgruppe"),
col = c("orange","blue"),
pch = 19,
lty = 1)
#####################
# Ergänzung um zwei Pfeile zur besseren Erklärung
m_1 = mean(daten_ancova$Pretest[daten_ancova$Gruppe == "Kontrolle"])
preds_1 = predict(object = modell1,
data.frame(Pretest = m_1,
Gruppe = c("Experimental", "Kontrolle")))
arrows(x0 = m_1, x1 = m_1, y0 = preds_1[2], y1 = preds_1[1])
Die Koeffizienten der Regression waren:
coef(modell1)## (Intercept) Pretest GruppeExperimental
## 26.7571429 0.4285714 -0.5714286Hier bedeuten nun:
Interceptgibt den Achsenabschnitt für die Kontrollgruppe an (also den Wert, an dem die orangene Gerade die \(Y\)-Achse schneidet). Ein signifikantes Ergebnis, bedeutet, dass der Achsenabschnitt der Kontrollgruppe statistisch signifikant von 0 abweicht.Pretestgibt die Steigung beider Regressionsgeraden an. Ein signifikantes Ergebnis bedeutet, dass die Steigung (in beiden Gruppen) statistisch signifikant von 0 abweicht.GruppeExperimentalgibt den Wert an, den man zur Regressionsgeraden der Kontrollgruppe (Referenzgruppe) addieren muss, um zur Regressionsgeraden der Experimentalgruppe zu gelangen (in der Abbildung illustriert durch einen Pfeil). Ein signifikantes Ergebnis bedeutet, dass die Regressionslinie für die Experimentalgruppe statistisch signifikant von der Regressionslinie der Kontrollgruppe (parallel) verschoben ist.
Zu beachten ist die implizite Bedeutung des Koeffizienten des Faktors. Der Wert \(-0.57\) gibt die vertikale Differenz der beiden Linien bzgl. des Kriteriums an. Da die Steigungen der Regressiongeraden in beiden Gruppen identisch ist, gilt diese Differenz allerdings für jeden beliebigen Prädiktorwert: Bei einem Pretest Wert von \(49\) beträgt die Differenz genauso \(0.57\) wie bei einem Pretest Wert von \(51\). Bei einem solchen Design vergleicht man also die Werte zweier Gruppen, für den Fall, dass diese vergleichbare Prädiktorwerte besäßen.
In diesem Zusammenhang spricht man auch von sog. “adjustierten Mittelwerten”. Die ursprünglichen Mittelwerte ohne Berücksichtigung der Kovariate lauten \(48.1\) und \(47.7\). Die sog. “adjustierten Mittelwerte” geben nun die “Mittelwerte” der Gruppen an, für den Fall, dass beide Gruppen einen vergleichbaren Wert auf der Kovariate (hier dem Pretest) hätten. Lassen Sie sich allerdings nicht verwirren: Statistisch gesehen sind dies keine wirklichen Mittelwerte, sondern die vorhergesagten Kriteriumswerte beider Gruppen, wenn man den Mittelwert der Kovariate betrachtet. Im folgenden Plot ist der Mittelwert der Kovariate als gestrichelte vertikale Linie eingezeichnet und die adjustierten Mittelwerte als Kreise.
# Plotten der empirischen Daten (orange = Kontrollgruppe, blau = Experimentalgruppe)
plot(jitter(Test, 0.3) ~ jitter(Pretest, 0.3),
data = daten_ancova,
col = ifelse(Gruppe == "Kontrolle", "orange", "blue"),
pch = 19,
ylab = "Test",
xlab = "Pretest")
# einzeichnen der Regressionsgeraden für beide Gruppen
coefs <- coef(modell1)
abline(coef = coefs[1:2], col = "orange")
abline(coef = c(coefs[1] + coefs[3], coefs[2]), col = "blue")
m_1 = mean(daten_ancova$Pretest)
abline(v = m_1, lty = 2)
adj_m = predict(object = modell1,
data.frame(Pretest = m_1,
Gruppe = c("Experimental", "Kontrolle")))
points(adj_m ~ rep(m_1,2), col = c("blue", "orange"), cex = 1.2)
legend(x = "bottomleft",
legend = c("Kontrollgruppe","Experimentalgruppe"),
col = c("orange","blue"),
pch = 19,
lty = 1)
Berechnen können wir diese adjustierten Werte mit der Funktion
predict(), wenn wir für die Variable Pretest
deren Mittelwert einsetzen:
predict(object = modell1, # Vorhersage auf Basis welches Modells?
data.frame(Pretest = mean(daten_ancova$Pretest), # für Mittelwert der Kovariate
Gruppe = c("Kontrolle", "Experimental"))) # und beide Gruppen## 1 2
## 48.18571 47.61429Alternativ können die adjustierten Mittelwerte mit der Funktion
emmeans() berechnet werden:
emmeans(modell1, pairwise ~ Gruppe)$emmeans## Gruppe emmean SE df lower.CL upper.CL
## Kontrolle 48.2 0.158 37 47.9 48.5
## Experimental 47.6 0.158 37 47.3 47.9
##
## Confidence level used: 0.95Eine genauere Betrachtung zum Zustandekommen der adjustierten Mittelwerte ist in einer Ergänzung zu Statistik II zu finden.
Derartige Situationen können natürlich auch mit anderen Funktionen
wie z.B. interact_plot() aus dem Paket
interactions erstellt werden (mehr Informationen gibt es hier):
interact_plot(modell1,
pred = Pretest,
modx = Gruppe,
plot.points = TRUE,
jitter = 0.1,
point.shape = TRUE,
interval = TRUE,
int.width = 0.95,
int.type = "confidence",
x.label = "Pretest",
y.label = "Test")
Alternativ können auch mit der Funktion predictorEffects()
Plots erstellt werden (ändern Sie ~Pretest auch zu
~Gruppe, um die resultierenden Unterschiede in der
Darstellung zu sehen). Da wir hier nur eine Steigung in beiden Gruppen
zugelassen haben (also keinen Interaktionsterm zwischen
Gruppe und Pretest zulassen), wird nur eine
Regressionsgerade gezeichnet (quasi beide Geraden gemittelt):
plot(predictorEffects(modell1, ~Pretest),
lines = list(multiline = FALSE))
Eine Besonderheit der hier betrachteten Daten (die wir bisher
verschwiegen hatten) war, dass Kovariate und Faktor nicht
(substanziell) korreliert waren. In anderen Worten unterscheiden
sich beide Gruppen nicht hinsichtlich der Werte auf der Variablen
Pretest:
aov.result.pretest <- aov(Pretest ~ Gruppe,
data = daten_ancova)
summary(aov.result.pretest)## Df Sum Sq Mean Sq F value Pr(>F)
## Gruppe 1 1.6 1.6000 2.714 0.108
## Residuals 38 22.4 0.5895Wir kommen nun zu dem Fall, dass Kovariate und Faktor korreliert sind, d.h. ein Unterschied zwischen beiden Gruppen auf der Kovariaten vorliegt.
8.3.3 Beispiel 2: Faktor und Kovariate korreliert
Mitunter steht geschrieben, dass eine Kovarianzanalyse auch genutzt
werden kann, um bestehenden Gruppenunterschieden entgegen zu wirken.
Würden also z.B. beide Gruppen in unserem Beispiel sich auf der
Kovariaten Pretest unterscheiden (also bereits vor der
eigentlichen Behandlung), dann wären Kovariate und Faktor konfundiert.
Das manchmal vorgebrachte Argument ist nun, dass die Kovarianzanalyse
diese Konfundierung quasi “herausrechnet” und somit ein statistisches
Verfahren zum Umgang mit solchen Situationen bereitstellt. Als Beispiel
nutzen wir nun den Datensatz in der Datei
ANCOVA_beispiel_2.dat, der inhaltlich dem Beispiel aus dem
vorangegangenen Abschnitt entspricht:
daten_ancova <- read.table(file = "./Daten/ANCOVA_beispiel_2.dat",
header = TRUE)
some(daten_ancova) # aus car: zeigt zufällige Datenpunkte anVisualisiert sehen diese Daten dann wie folgt aus:
und eine Varianzanalyse zeigt nun auch den Unterschied auf der Variablen
Pretest an:
aov.result.pretest <- aov(Pretest ~ Gruppe,
data = daten_ancova)
summary(aov.result.pretest)## Df Sum Sq Mean Sq F value Pr(>F)
## Gruppe 1 105.63 105.63 305.2 <2e-16 ***
## Residuals 38 13.15 0.35
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Das heißt, Personen der Experimentalgruppe haben einen höheren Pretest Wert als Personen in der Kontrollgruppe.
Natürlich können wir auch die Mittelwerte berechnen und dann eine Regression durchführen und uns entsprechend die Koeffizienten sowie deren Inferenzstatistik ausgeben lassen:
# Mittelwerte berechnen
Tapply(Test ~ Gruppe,
fun = mean,
data = daten_ancova)## Experimental Kontrolle
## 49.40 48.85daten_ancova$Gruppe <- as.factor(daten_ancova$Gruppe)
daten_ancova$Gruppe <- relevel(daten_ancova$Gruppe,
ref = "Kontrolle")
# Regression
modell2 = lm(Test ~ Pretest + Gruppe,
data = daten_ancova)
S(modell2, brief = TRUE)## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.9650 10.3586 2.603 0.0132 *
## Pretest 0.4373 0.2069 2.113 0.0414 *
## GruppeExperimental -0.8711 0.7132 -1.221 0.2297
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 0.7504 on 37 degrees of freedom
## Multiple R-squared: 0.21
## F-statistic: 4.918 on 2 and 37 DF, p-value: 0.01276
## AIC BIC
## 95.43 102.18Die Interpretation ist genau wie für Beispiel 1, allerdings tritt
hier sogar etwas eher Unangenehmes auf: Die Mittelwerte sind \(M=48.85\) für die Kontrollgruppe und \(M=49.4\) für die Experimentalgruppe (das
Medikament wäre also eher kontraproduktiv). Der Koeffizient
GruppeExperimental besagt aber nun, dass wir einen Wert von
\(-0.8711\) zur Kontrollgruppe addieren
müssen, um zur Experimentalgruppe zu gelangen.
Der Grund ist in den unterschiedlichen Ausgangssituationen zu finden:
Die Experimentalgruppe hat zwar höhere Krankheitswerte nach der
Medikamenteinnahme, besitzt aber auch schon höhere Krankheitssymtome zum
Zeitpunkt, als die Variable Pretest erhoben wurde.
Vergleicht man also die reinen Mittelwerte der Experimental- und
Kontrollgruppe, so zeigt die Experimentalgruppe einen höheren Testwert.
Die Situation verändert sich allerdings, wenn wir die Pretest Werte als
weitere Variable berücksichtigen. Durch die Hinzunahme der Kovariaten
Pretest, können wir die Differenzen beider Gruppen
ausrechnen, für den Fall beide Gruppen hätten vergleichbare
Pretest Werte (quasi: wir betrachten die Differenz der
beiden Gruppen “unter Bereinigung der unterschiedlichen
Pretest-Werte”). Zum Beispiel könnten wir wieder die
“adjustierten Mittelwerte” berechnen, also die Modellvorhersagen, wenn
sich die Personen auf der Kovariaten nicht unterscheiden.
emmeans(modell2, pairwise ~ Gruppe)$emmeans## Gruppe emmean SE df lower.CL upper.CL
## Kontrolle 49.6 0.376 37 48.8 50.3
## Experimental 48.7 0.376 37 47.9 49.5
##
## Confidence level used: 0.95Nun sehen wir, dass in diesem Fall die Experimentalgruppe geringere
Krankheitssymptome aufweist als die Kontrollgruppe. Dabei ist wichtig zu
verstehen, dass dieser Vergleich statistischer Natur ist, unter linearer
Extrapolation der Werte. Denn faktisch liegen keine vergleichbaren
Pretest-Werte vor und wir können nur abschätzen, wie die
Werte der Experimental- bzw. Kontrollgruppe wären, hätten beide
vergleichbare Pretest-Werte. Insofern scheint uns die
Kovarianzanalyse als statistisches Werkzeug zur Eliminierung von
Konfundierungen ohne Wissen darüber was hier genau verglichen wird, eher
mit Vorsicht zu genießen zu sein (auch wenn es darüber in der Statistik
auch Diskussionen gibt).
8.3.4 Beispiel 3: Interaktion zwischen Faktor und Kovariate
Bisher waren wir immer davon ausgegangen, dass die Steigung der Regressionsgeraden in den Gruppen gleich ist. Diese Annahme wird auch als Homogenität der Regressionssteigungen bezeichnet und ist typisch für den Fall der Kovarianzanalyse. Statistisch bedeutet dies, dass wir die Interaktion nicht in das Modell aufgenommen haben.
In diesem letzten Beispiel zeigen wir nun noch, dass diese Annahme
natürlich dazu führen kann, tatsächliche Beziehungen zwischen Kovariate
und Kriterium zu verfehlen. Dazu nutzen wir nun die Daten in der Datei
ANCOVA_beispiel_3.dat, die wiederum die gleiche Situation
wie in den vorherigen Beispielen beschreiben:
daten_ancova <- read.table(file = "./Daten/ANCOVA_beispiel_3.dat",
header = TRUE)
daten_ancova$Gruppe <- as.factor(daten_ancova$Gruppe)
daten_ancova$Gruppe <- relevel(daten_ancova$Gruppe,
ref = "Kontrolle")Wie allerdings schnell im linken Teil der folgenden Abbildung zu sehen ist, liegt in der Experimentalgruppe ein negativer Zusammenhang und in der Kontrollgruppe ein positiver Zusammenhang vor. Beide Zusammenhänge würden durch ein Modell ohne Interaktion (fast) nicht erfasst werden. Lassen wir die Interaktion allerdings zu, dann erhalten wir eine deutlich bessere Schätzung.
Da wir nun separate Regressionssteigungen für beide Gruppen schätzen
(durch die Interaktion, siehe weiter unten), würden über die Funktion
predictorEffects() auch wieder zwei Abbildungen für jeweils
eine der Gruppen resultieren:
plot(predictorEffects(lm.result.ia, ~Pretest),
lines = list(multiline = FALSE))
Nun berechnen wir beide Modell, d.h. einmal ohne und einmal mit Interaktion und führen einen Modellvergleich durch, aus dem deutlich wird, dass die Hinzunahme der Interaktion das Modell verbessert hat.
modell.ohne.ia <- lm(Test ~ Pretest + Gruppe,
data = daten_ancova)
modell.mit.ia <- lm(Test ~ Pretest * Gruppe,
data = daten_ancova)
anova(modell.ohne.ia, modell.mit.ia)Schließlich können wir uns natürlich auch für das Modell mit Interaktion die Koeffizienten und ihre Interpretation anschauen:
coef(modell.mit.ia)## (Intercept) Pretest
## 31.8826816 0.3240223
## GruppeExperimental Pretest:GruppeExperimental
## 43.1173184 -0.8240223Interceptmeint hier den Achsenabschnitt der Kontrollgruppe.Pretestist die die Steigung der Kontrollgruppe.GruppeExperimentalist der Wert, der auf den Achsenabschnitt der Kontrollgruppe addiert werden muss, um zum Achsenabschnitt der Experimentalgruppe zu gelangen (\(31.88 + 43.12 = 75\)).Pretest:GruppeExperimentalist der Wert, der auf den Steigungskoeffizienten der Kontrollgruppe addiert werden muss, um zur Steigung der Experimentalgruppe zu gelangen (\(0.324 + (-0.824) = -0.5\)).
8.4 Spezialfälle bei Prädiktoren und Hinzunahme von Prädiktoren
8.4.1 Multikollinearität und Varianz-Inflations-Faktor
Kollinearität meint im Zuge der multiplen linearen Regression eine perfekte Korrelation zweier Prädiktoren; Multikollinearität meint, wenn mehrere Prädiktoren miteinander perfekt korreliert sind. (Multi-)Kollinearität ist zunächst ein technisches Problem, wenn zwei Prädiktoren vollständig korreliert sind. In diesem Fall sind die Prädiktoren linear abhängig und das relevante Gleichungssystem ist nicht lösbar. Allerdings spricht man auch von (Multi-)Kollinearität, wenn mehrer Prädiktoren zwar nicht perfekt, aber sehr stark korreliert sind, was zu Problemen für der Präzision der Parameterschätzungen führt (v.a. in Kombination mit kleinen Stichproben). Erkennbar wird dies durch einen größeren Standardfehler der Parameter, welche mit steigender (Multi-)Kollinearität größer werden.
Dazu kommt noch, dass es zunehmend schwer wird, den individuellen Beitrag einzelner Prädiktoren zu ermitteln; was durchaus auch im Fokus des Forschungsinteresses stehen kann. Allerdings ist dies nicht nur ein statistisches Problem, sondern eines, welches bei Kenntnis des Forschungsgegenstandes auch gelöst werden kann. Angenommen, jemand ermittelt die Mathematikleistung (als Kriterium) sowie das Alter und die Schulstufe von Kindern. Vermutlich wird es zwischen allen drei Variablen hohe Korrelationen geben. Eine Regression mit Alter und Schulstufe als Prädiktoren kann nun keine Aussage darüber geben, ob Alter oder Schulstufe wichtiger ist; eine inhaltliche Betrachtung würde aber vermutlich darauf hindeuten, dass nicht das bloße Vergehen von Zeit (also das Alter) die wichtige Variable ist, sondern die Schulstufe.
Um das Problem der (Multi-)Kollinearität zu verdeutlichen, betrachten wir eine multiple lineare Regression mit zwei Prädiktoren. Im folgenden linken Plot ist ein Scatterplot gezeigt, bei dem die Prädiktoren stark korreliert sind. Geometrisch zeigt sich dies dadurch, dass die Punkte sich nicht breit verteilen, sondern einen Art “Schlauch” bilden. Mathematisch können wir zwar nun ohne Probleme eine Regressionebene durch die Punktwolke ziehen und somit Koeffizienten für die Regression erhalten. Da sich die Punkte allerdings nicht breit verteilen, ist die Lösung “instabil”. Das heißt relativ kleine Veränderungen der Punktwolke, beispielweise beim ein erneuten Ziehen von Daten kann zu einer deutlichen Veränderung der Regressionsebene und damit ihrer Koeffizienten führen. Hierfür können wir den rechten Plot betrachten, für den Daten erzeugt wurden, die aus der gleichen Population stammen wie die Daten aus dem linken Plot. Durch das Ziehen anderer Werte gibt es allerdings kleinere Unterschiede in den konkreten Ausprägungen. Demgegenüber ist die Auswirkung auf die Regressionsebene erstaunlich groß:
library(scatterplot3d)
library(MASS)
angle <- 220
set.seed(1)
mu <- c(5, 5)
Sigma <- matrix(c(1, 0.98, # beide Prädiktoren stark korreliert
0.98, 1),
nrow = 2,
ncol = 2)
sample <- mvrnorm(15,
mu,
Sigma,
empirical = TRUE)
daten <- data.frame(VP = 1:15,
x_1 = sample[, 1],
x_2 = sample[, 2])
daten$y = 3 + 2 * daten$x_1 + 1 * daten$x_2 + rnorm(15)
par(mfrow = c(1, 2))
lm1 <- lm(y ~ x_1 + x_2,
data = daten)
s3d <- scatterplot3d(
daten$x_1,
daten$x_2,
daten$y,
xlab = "X1",
ylab = "X2",
zlab = "Y",
ylim = c(0, 8),
xlim = c(0, 8),
zlim = c(5, 30),
pch = 16,
color = "red",
highlight.3d = FALSE,
type = "h",
main = "Stichprobe 1",
grid = TRUE,
box = FALSE,
angle = angle
)
s3d$plane3d(lm1)
# nun wird eine zweite Stichprobe aus der gleichen Population gezogen
set.seed(11)
sample <- mvrnorm(15,
mu,
Sigma,
empirical = TRUE)
daten <- data.frame(VP = 1:15,
x_1 = sample[, 1],
x_2 = sample[, 2])
daten$y = 3 + 2 * daten$x_1 + 1 * daten$x_2 + rnorm(15)
lm1 <- lm(y ~ x_1 + x_2,
data = daten)
s3d <- scatterplot3d(
daten$x_1,
daten$x_2,
daten$y,
xlab = "X1",
ylab = "X2",
zlab = "Y",
ylim = c(0, 8),
xlim = c(0, 8),
zlim = c(5, 30),
pch = 16,
color = "red",
highlight.3d = FALSE,
type = "h",
main = "Stichprobe 2",
grid = TRUE,
box = FALSE,
angle = angle
)
s3d$plane3d(lm1)
In anderen (wenig mathematischen) Worten: Bei (Multi-)Kollinearität werden Regressionsebenen durch die empirischen Daten (= die Punktwolke) nur schlecht austariert und werden instabil. Wenn wir also immer und immer wieder eine Stichprobe ziehen, so ändern sich die Regressionsebenen relativ gesehen stark und die einzelnen Regressionskoeffizienten schwanken entsprechend von Stichprobe zu Stichprobe.
Die zwei wahrscheinlich sinnvollsten Maßnahmen gegen (Multi-)Kollinearität sind eine ausreichend große Stichprobe (zur Erhöhung der Power der einzelnen Tests über die Parameter) und das Vermeiden des Einbezugs von korrelierten Prädiktoren. Kommen diese beiden Maßnahmen nicht in Frage, kann ggf. das Regressionsmodell angepasst werden, indem Prädiktoren weggelassen werden oder die korrelierten Prädiktoren zu einem einzigen Prädiktor “verrechnet” werden – auch wenn diese Maßnahmen besonderer Umsicht bedürfen und sicherlich nicht generell zu empfehlen sind.
Kollinearität kann durch eine Korrelationsmatrix entdeckt werden. Hierbei kann Multikollinearität allerdings leicht übersehen werden. Wie groß die geteilte Variabilität eines Prädiktors \(p\) mit den restlichen \(q-1\) Prädiktoren ist, wird berechnet, indem der Prädiktor \(X_p\) als Kriterium einer multiplen Regression mit den verbleibenden Prädiktoren als Prädiktoren verwendet wird. Die gemeinsame Varianz wird dann in Form von \(R^2\) quantifiziert, oft als \(R^2_{X_p}\) bezeichnet und zur Toleranz \(T\) verrechnet: \[T = 1-R^2_{X_p}\] \(T\) wird 1, wenn ein \(X_p\) orthogonal zu den anderen Prädiktoren ist und wird 0, wenn \(X_p\) perfekt aus den anderen Prädiktoren vorhersagbar ist. Als Varianz-Inflations Faktor (VIF) wird bezeichnet: \[VIF_{X_p}=\frac{1}{1-R^2_{X_p}}=\frac{1}{T}\] Dieser wiederum geht in die Berechnung des geschätzten Standardfehlers \(S_{b_p}\) eines Prädiktors \(p\) mit ein:
\[S_{b_p} = \frac{S_E^2}{(n-1)S_p^2} \cdot VIF_{X_p} \]
Eine Interpretation von VIF ist, dass der Standardfehler eines Prädiktors um \(\sqrt{VIF}\) relativ zu seinem theoretischen Minimum bei vollständiger Abwesenheit von Multikollinearität aufgebläht wird. Es gibt aber keine klaren Angaben dazu, ab welchen Werten für \(VIF\) die Multikollinearität problematisch wird. Manchen Empfehlungen zu Folge gilt dies ab \(VIF\geq 10\) und wenn der durchschnittliche Wert für \(VIF\) “deutlich” größer als 1 ist.
Die praktische Berechnung erfolgt sehr einfach mit der Funktion
vif() aus dem Paket car, der einfach ein
lm()-Objekt übergeben wird. Hier zur Illustration einmal
mit den drei Prädiktoren des Beispiels aus Teil 8:
daten <- read.table("./Daten/Album_Verkaeufe.dat",
header = TRUE)
modell <- lm(Verkaeufe ~ Werbung + Frequenz + Attraktivitaet,
data = daten)
vif(modell) # Einzelwerte## Werbung Frequenz Attraktivitaet
## 1.014593 1.042504 1.038455mean(vif(modell)) # Gesamtmittelwert## [1] 1.03185Zur weiteren Exploration der Auswirkungen von (Multi-)Kollinearität finden Sie hier eine ShinyApp.
8.4.2 Suppression
Manchmal redet man bei gewissen Variablen von Suppressorvariablen, die in bestimmten Konstellationen einer Regression auftreten können. Die ist die Situation, dass ein Prädiktor in ein multiples Regressionsmodell aufgenommen wird, der an sich keinen oder nur einen geringen Zusammenhang mit dem Kriterium aufweist, allerdings in Kombination mit anderen Prädiktoren die Modellvorhersage deutlich verbessert. Solch eine Variable wird als Suppressorvariable bezeichnet. Hierbei kann es sogar vorkommen, dass sich das Vorzeichen von Prädiktoren verändert, sobald die Suppressorvariable aufgenommen wird.
Zur Veranschaulichung dieses Phänomens wollen wir den Datensatz
supression.dat betrachten. Modelliert werden soll hierbei
die von Schüler:innen erreichte Punktzahl in Abhängigkeit der Anzahl an
Lernstunden sowie der Ablenkbarkeit.
daten_sup <- read.table(file = "./Daten/supression.dat",
header = TRUE)
some(daten_sup)Zunächst berechnen wir eine Regression von Punktzahl auf
die Lernstunden und erwarten, dass eine höhere Anzahl an
Lernstunden mit einer höheren Punktzahl einhergeht. Tatsächlich zeigt
sich dieser Zusammenhang und die Güte der Regression ist \(R^2 = 0.15\). Statistisch signifikant ist
der Test jedoch nicht.
reg1 <- lm(Punktzahl ~ Lernstunden,
data = daten_sup)
S(reg1, brief = TRUE)## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.0650 2.3541 2.576 0.0190 *
## Lernstunden 0.1913 0.1054 1.815 0.0863 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 4.637 on 18 degrees of freedom
## Multiple R-squared: 0.1546
## F-statistic: 3.292 on 1 and 18 DF, p-value: 0.0863
## AIC BIC
## 122.01 125.00Analog berechnen wir nun eine Regression von Punktzahl
auf Ablenkbarkeit und erwarten, dass eine höhere
Ablenkbarkeit mit einer geringeren Punktzahl einhergeht. Tatsächlich
zeigt sich aber überraschenderweise auch hier ein deskriptiv positiver
Zusammenhang, sodass eine höhere Ablenkbarkeit scheinbar mit einer
höheren Punktzahl einhergeht. Allerdings ist der Zusammenhang nicht
statistisch signifikant und \(R^2 =
0.004\) ist sehr gering (auch die Korrelation beider Variablen
fällt natürlich mit \(0.06\) sehr
gering aus):
reg2 <- lm(Punktzahl ~ Ablenkbarkeit,
data = daten_sup)
S(reg2, brief = TRUE)## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.1536 2.9377 3.116 0.00597 **
## Ablenkbarkeit 0.1059 0.3849 0.275 0.78639
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 5.033 on 18 degrees of freedom
## Multiple R-squared: 0.004186
## F-statistic: 0.07567 on 1 and 18 DF, p-value: 0.7864
## AIC BIC
## 125.29 128.27Betrachtet man allerdings nun beide Variablen gemeinsam, so verbessert sich die Modellgüte substanziell (mit \(R^2 = 0.33\)) und der Effekt der Lernstunden wird signifikant
reg3 = lm(Punktzahl ~ Lernstunden + Ablenkbarkeit,
data = daten_sup)
S(reg3, brief = TRUE)## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.6399 2.4944 3.464 0.00297 **
## Lernstunden 0.4575 0.1608 2.845 0.01118 *
## Ablenkbarkeit -1.1223 0.5409 -2.075 0.05349 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard deviation: 4.262 on 17 degrees of freedom
## Multiple R-squared: 0.3255
## F-statistic: 4.101 on 2 and 17 DF, p-value: 0.0352
## AIC BIC
## 119.50 123.48Auch der Zusammenhang zwischen der Ablenkbarkeit und der
erreichten Punktzahlen ist nun negativ, sodass (unter
Konstanthaltung der Lernstunden) eine höhere
Ablenkbarkeit zu einer geringeren Punktzahl
führt.
Zusammengefasst zeigt dieses Beispiel also, dass
Ablenkbarkeit kaum bzw. sogar schwach positiv mit der
erreichten Punktzahl zusammenhängt, wenn man die Variable isoliert
betrachtet. Auch die Anzahl an Lernstunden hängt alleine
nicht statistisch signifikant (wenn auch deskriptiv) mit der erreichten
Punktzahl zusammen. Betrachtet man allerdings beide
Variablen gemeinsam, so verbessert sich die Modellgüte substanziell
gegenüber den einzelnen Regressionen (vgl. \(R^2 = 0.004\) bzw. \(R^2 = 0.15\) gegenüber \(R^2 = 0.33\)) und auch der Zusammenhang
zwischen Ablenkbarkeit und Punktzahl verändert
sich gegenüber der einfachen lineare Regression.
Eine Variable die sich so verhält wie die Ablenkbarkeit
in diesem Beispiel nennt man einen Suppressor. Wieso sich diese Variable
so auswirken kann, ist oft schlecht nachvollziehbar ohne die tiefere
Mathematik der Regression zu verstehen. Auch verbale Beschreibungen oder
sog. Venn-Diagramme kommen schnell an ihre Grenzen. Eine sinnvolle
Anschauung liefert allerdings die vektorgeometrische Betrachtung der
Regression im Allgemeinen Linearen Modell (siehe z.B. Wickens, 2014).
Zur weiteren Exploration der Auswirkungen von Suppression finden Sie hier
eine ShinyApp.
8.4.3 Methoden der Hinzunahme von Prädiktoren
Gibt es mehrere Prädiktoren, stellt sich die Frage, ob tatsächlich alle Prädiktoren gleichzeitig in ein Modell aufgenommen werden sollen oder, ob sukzessive mehrere Modelle mit zunehmend mehr Prädiktoren Berücksichtigung finden sollen. Im Folgenden werden verschiedene Herangehensweisen beschrieben, ohne diese werten zu wollen: Verschiedene Forschungstraditionen scheinen hier verschiedene Standardvarianten zu nutzen und das Verfahren der Wahl mag auch ein inhaltlich begründetes sein.
- Alle Prädiktoren gemeinsam: Bei dieser Herangehensweise werden alle Prädiktoren gleichzeitig in ein Modell aufgenommen. Die Wahl der Prädiktoren erfolgt selbstverständlich aus inhaltlichen Gründen, aber es wird dann keine Reihenfolge festgelegt (vgl. Punkt 2). Es gibt durchaus Positionen, die dies als die einzig sinnvolle Art ansehen.
- Hierarchische Regression: Hierbei wird auf Basis vorheriger Studien entschieden, welche Prädiktoren zuerst ein Modell aufgenommen werden sollen und welche danach. Zwischen den Modellen kann mit \(F\)-Tests der \(R^2\)-Zuwachs getestet werden. Eine Variante ist es, hierbei mit den Prädiktoren zu beginnen, von denen man weiß, dass diese einen starken Einfluss auf das Kriterium haben.
- Schrittweise Regression: Bei einer schrittweisen Regression wird über die Hinzunahme eines Prädiktors rein nach einem mathematischen Kriterium entschieden. Diese Vorgehensweise wird wegen ihrer Theorielosigkeit skeptisch betrachtet, aber sei der Vollständigkeit halber kurz erwähnt. Bei einer Vorwärtsregression wird mit einem Modell begonnen, welches nur den Achsenabschnitt berücksichtigt und dann sucht der Computer nach demjenigen Prädiktor, der am meisten Variabilität des Kriteriums aufklärt. Erhöht dieser Prädiktor \(R^2\), wird er beibehalten und der nächstbeste Prädiktor wird gesucht und hinzugefügt. Eine Rückwärtsregression beginnt mit dem vollständigsten Modell und entfernt dann Prädiktoren und testet, ob \(R^2\) kleiner wird. Es gibt auch noch eine beidseitige Variante, bei der ein Modell ständig neu bewertet wird und redundante Prädiktoren entfernt werden.
8.5 Literatur
Wickens, T. D. (2014). The geometry of multivariate statistics. Psychology Press.