In diesem Artikel beschreiben wir, wie sich eine Regressionsanalyse in R umsetzen lässt.
Die Daten
Für die Analysen nutzen wir den Datensatz wage2 aus dem Paket wooldrige. Hierzu müssen wir zunächst das Paket wooldrige mittels library(wooldridge) in R laden. Falls das Paket bei Euch nicht geladen wird, so liegt es daran, dass es nicht installiert ist. Mit install.packages("wooldridge") installiert ihr es in R. Der Datensatz wage2 befindeet sich nach dem Laden des Pakets wooldridge im Workspace von R. Der Datensatz wage2 umfassst Daten zu N = 935 Berufstätigen 17 verschiedene Merkmale. Wir werden hier die Variablen lwage das logarithmierte monatliche Einkommen (abhängige Variable) und die unabhängigen Variablen hours (durchschnittliche wöchentliche Arbeitsstunden), IQ (Intelligenzquotient der Befragten), KWW (Ein Test zu arbeitsbezogenen Fähigkeiten), educ (Jahre an Schulbilidung), exper (Jahre an Berufserfahrung), tenure (Jahr bei aktuellem Arbeitgeber) sowie age (Alter in Jahren)verwenden. Für die genaue Beschreibung des Datensatzes ist auf seine Dokumentation in R verwiesen. Ihr erhaltet sie, wenn ihr in R ?wage2 eingebt.
Die Analyse
Wir möchten den lwage durch die unabhängigen Variablen hours, IQ, KWW, educ, exper, tenure und age in einer multplen linearen Regression erklären. Mit der Funktion lm können lineare Regressionen in R berechnet werden. Der Befehl für unser Modell, welches wir als Objekt model abspeichern, lautet:
Mit der Funktion summary angewendet auf unser Modell erhalten wir eine Modellzusammenfassung in Form von p-Werten, Bestimmtheistmaßen und weiteren Kennzahlen. Wir haben uns nicht die Ergebnisse dieser Regression ausgeben lassen, da wir zunächst die Annahmen des Modells überprüfen uund daraufhin das Modell gegebenenfalls anpassen werden.
Zunächst prüfen wir die Linearität. Ein Signifikanztest, welcher auf Linearität in einem Regressionsmodell testet, ist der Rainbow Test. Dieser ist im Paket lmtest implementiert. Die zugehörige Funktion lautet raintest. Mit raintest(model) berechnet R die Ergebnisse.
Obiges Bild zeigt die Ergebnisse des Rainbow Tests. Der Rainbow Test liefert hierbei ein nicht-signifikantes Ergebnis, F(468, 459) = 0,97, p = 0,622. Somit konnte kein Indiz für eine nicht vorliegende Linearität beobachtet werden. Die Linearität kann in diesem Modell somit angenommen werden.
Weiterhin müssen wir unser Modell auf Autokorrelation untersuchen. Dies können wir mittels des Duerbin-Watson-Test machen. Dieser ist im Paket car mit der Funktion durbinWatsonTest implementiert. Der Durbin-Watson-Test liefert hierbei ein signifikantes Ergebnis, DW = 1,78, p = 0,002. Dies spricht zunächst für ein Autokorrelationsproblem. Betrachten wir jedoch die Korrelation der Residuen zum Lag 1, so zeigt sich, dass hier eine schwach Korrelation, r = 0,11 vorliegt. Dieser Wert ist ebenfalls in unterem Output aufgeführt.
Es ist also anzunehmen, dass die Autokorrelation der Residuen eher schwach ist. Dies zeigt sich ebenfalls in der unteren Grafik. Hier sind die Residuen zum Lag 1 gegen die Residuen abgetragen. Wenn überhaupt würde man eine Gerade mit einer vergleichsweise schwachen positiven Steigung durch diese Punktwolke legen, so wie es die geringe Korrelation schon angedeutet hat. Unter diesen Gesichtspunkten ist ein Autokorrelationsproblem in diesem Modell eher auszuschließen, weswegen die Annahme der Unkorrliertheit der Residuen, also keine Autokorrelation, und damit eine Erfüllung der Regressionsvoraussetzung angenommen wird. Das signifikante Ergebnis bei dem Durbin-Watson-Test ist vermutlich auf die hohe Stichprobengröße und die damit einhergenehnde hohe Power zurückzuführen. Sprich: Trotz signifikantem Ergebnis ist die Korrelation marginal und damit praktisch nicht relevant.
Weiterhin prüfen Ihr nun die Annahme der Homoskedastizität. Mittels des Breusch-Pagan-Tests kann dies erfolgen. Der Bresuch-Pagan-Test ist in R im Paket lmtest mit der Funktion bptest implementiert. Wir erhalten das Ergebnis mit bptest(model). Kommender Output zeigt die Resultate. Hierbei liefert der Breusch-Pagan-Test χ2(7) = 44,09, p = 0,000.
Mittels des Befehls plot(fitted(model), residuals(model), xlab = "geschätzte Werte", ylab = "Residuen") erhalten wir einen Residualplot. Dieser ist in kommender Grafik abgebildet. Die Streuung wirkt naehzu gleichmäßig. Hier sind wir datenanalystisch in der selben Situation, wie bei der Prüfung auf Autokorrelation. Das signifikante Ergebnis bei dem Breusch-Pagan-Test ist vermutlich auf eine Hohe Power zurückzuführen. Somit wird also Homoskedastizität festgestellt und eine Verletzung dieser Annahme ausgeschlossen.
In R kann das Modell bezüglich diverser Kriterien weiter untersucht werden. Zum Beispiel kann geprüft werden, ob vergessene Einflussgrößen vorliegen. Dies ist beispielsweise mit der RESET-Test nach Ramsey möglich. In R erhalten Sie ihn mit dem Befehl reset(model). Dieser liefert ein nicht-signifikantes Ergebnis, F(3, 63) = 1,70, p = 0,176. Somit liegen keine vergessenen Regressoren vor.
Nun wird die Annahme der Normalverteilung untersucht. Hierfür müssen die Residuen untersucht werden. Dies gelingt Euch mittels eines Q-Q-PLots. Die Residuen erhalten wir in R über residuals(model). Mit dem Befehl qqnorm(residuals(model), xlab = "Quantile der Normalverteilung", ylab = "Quantile der Residuen", main = NULL) lässt sich dann ein Q-Q-Plot für die Residuen erstellen, bei welchem die Quantile der Reisduen mit denen einer Normalverteilung verglichen werden. Kommende Grafik zeigt uns einen solchen Plot. Mit qqline(residuals(model)) wird eine Hilfslinie in die Grafik eingezeichnet. Es sind hierbei an den Rändern leichte Abweichungen von einer Normalverteilung zu erkennen.
Kommendes Bild zeigt Euch nochmal die genutzte Code-Folge für die Erstellung des Q-Q-Plots.
Abschließend prüfen wir unser Modell noch auf ein etwaiges Multikollinearitätsproblem. Dies gelingt uns mittels der Varianzinflationsfaktoren (VIF). Diese erhalten wir mittels des Befehls vif(model) in R.
Es zeigt sich, dass alle VIF kleiner 10 sind. Somit ist von keiner Multikollinearität zwischen unseren Prädiktoren auszugehen.
Alle Regressionsannahme waren soweit (annähernd) erfüllt. Es müssen also hier keine weiteren datenanalytische Maßnahmen erfolgen. Somit kann nun das Modell untersucht werden. Mit summary(model) erhalten wir in R eine Zusammenfassung der zentralen Regressionsergebnisse.
Es zeigt sich, dass der Effekt der durchschnittlichen Wochenstunden hours auf dem 5% Niveau signifikant negativ ist. Die Effekte der Arbeitsbezogenen Fähigkeiten KWW, der Intelligenzquotient der Befragten IQ, die Jahre an Schulbilidung educ, die Jahre an Arbeitserfahrung exper und die Jahre beim aktuellen Arbeitgeber tenure haben jeweils auf dem 5% Niveau einem signifikant posiitven Effekt auf das logarithmierte Einkommen. Der Effekt des Alters ist positiv und nicht-signifikant. Weiterhin liegt ein adjustiertes Bestimmtheitsmaß von 0,195 vor. Somit werden durch dieses Modell 19,5% der Varianz des logarithmierten Einkommens erklärt.
