-
Notifications
You must be signed in to change notification settings - Fork 1
Regressionsanalyse mit dem Commander
Dieses File soll die generelle Vorgehensweise bei einer Regressionsanalyse in R wiedergeben. Es geht nicht darum das beste Modell aus den Daten zufinden, sondern ein aus der Theorie abgeleitetes Modell anhand von empirisch gemessen Daten zu überprüfen. Das Modell aus der Theorie soll für uns sein: Benzinverbrauch lässt sich linear über das Gewicht und die Pferdestärke vorrausagen.
Es gibt zwei Möglichkeiten die relevanten Pakete zu laden. Direkt über den library Befehl:
# Zum Laden des Paketes *car*
library(car)
oder über die oberere Menüleiste unter packages und loadpackages. Dann erscheint folgendes Menü:
Wir benötigen nur die beiden Pakete QuantPsyc und car
Der R-Commander ist ein weiteres Paket. Der Paketname lautet Rcmdr. Er kann durch die beiden oben genannten Möglichkeiten geöffnet werden.
# R-Commander
library(Rcmdr)
Es gibt viele Möglichkeiten Datenmatrixen in den R-Commander zu laden.Es hängt davon ab, wo und wie die Daten vorhanden sind. Daten können z.B. als CSV-Format vorliegen. Mit Excel lassen sich Daten in diesem Format abspeichern. Diese Daten können unter Datenmanagment unter Importiere Dateien gefunden werden.  Darüber gelagt man zu einem Menü. Dort muss man angeben, ob man die Daten auf dem lokalen Systemgespreicht hat oder man sie direkt von einer Internetseite ziehen möchte. Desweiteren muss man das Format der Daten definieren und Angaben über das Datentreenfeldzeichen und das Dezimaltrennfeldzeichen machen. Auch kann man den Namen der Datenmatrix festlegen.  Hat man dies alles getan öffnet sich ein Explorerfenster. Dort wählt man dann die entsprechende Datei aus.
Es können aber auch Daten geladen werden, welche bereits in R vorinstalliert sind. Dazu geht man zu Daten in Paketen und kann die Daten mit dem Befehl Lese Datenmatrix aus geladenen Paketen bekommen.
Wir benötigen die Daten mtcars aus dem Paket datasets.
Der RCode lautet dazu:
# Datenmatrix mtcars
data(mtcars, package="datasets")
Bevor man anfängt das Regressionsmodell zu erstellen. Kann man sich bereits anschauen, ob sich ein Zusammenhang in den Daten graphisch zeigt. Auch kann man so die Daten auf Homoskedastizität testen und schon einen ersten Blick auf Ausreißer werden. Dazu zeichnet man für jede der Variablen (UVs), welche zur Vorhersage im Modell gelten sollen, ein Streudiagramm mit der vorherzusagenden Varibale (AV). In unserem Fall sind die UVs im Datensatz: hp (Horsepower also PS) und wt (weight, steht für Gewicht). Dazu muss man im Menü unter Grafik auf Streudiagramme gehen.  Dann erscheint folgendes Menü:  In diesem Menü setzen wir als X-Achse hp und als y-Achse mpg.
Dann klicken dort weiter auf Optionen und haben viele Möglichkeiten. 
Wir setzen als Plot-Option ein Häckchen bei Kleinste-Quadrat Linie und Beschriften die Achsen des Graphen. Der Commander generiert uns daraufhin folgenden Befehl und dieser kreiert uns dann die nachfolgende Graphik.
# Streudiagramm Pferdestärke zu Verbrauch
scatterplot(mpg~hp, reg.line=lm, smooth=FALSE, spread=FALSE,
id.method='mahal', id.n = 2, boxplots=FALSE, span=0.5, xlab="Pferdestärle",
ylab="Meilen pro Gallone", main="Pferdestärke zum Verbrauch", data=mtcars)
In der Graphik kann man erkennen, dass es einen linearen Zusammenhang zugeben scheint. Man sieht aber auch, dass es zwei Ausreißer aus den Daten gibt. Diese beiden Automodelle nennt uns R in der Grafik direkt. ( Toyota Corolla und Maserati Bora). Zusätzlich gibt uns R die Spaltennummer der Daten an. (Toyota 20 und Maserati 31). Diese werden wir uns später nochmal anschauen. Nun erstellen wir dieselbe Grafik mit dem Gewicht als UV. Wir ändern die Einstellungen in dem eben besprochenden Menü und erhalten folgenden Befehl und folgenden Grafik:
# Streudiagramm Gewicht zu Verbrauch
scatterplot(mpg~wt, reg.line=lm, smooth=FALSE, spread=FALSE,
id.method='mahal', id.n = 2, boxplots=FALSE, span=0.5, xlab="Gewicht",
ylab="Meilen pro Gallone", main="Gewicht zum Verbrauch", data=mtcars)
Auch hier scheint es einen Zusammenhang zugeben. Es fallen auch wieder zwei Ausreißer auf.
Zur Erstellung unseres Regressionsmodells klicken wir im Commander in das Menü Statistik/Regressionsmodelle/LineareRegression.  Dann erscheint folgendes Fester: 
Das Modell nennen wir Verbrauch. Als Abhängige Variable wählen wir mpg und als Unabhängige Variablen wt und hp. Der Commander generiert uns daraufhin folgenden Befehl:
# Erstellung des Models
Verbrauch <- lm(mpg~hp+wt, data=mtcars)
summary(Verbrauch)
Die Befehlszeile lässt sich wie folgt lesen: Das Modell heißt Verbrach und es wird durch den LM Befehl generiert. Die AV mpg steht vor dem geschwungenen Minus. Die zur Vorhersage genutzten Variablen stehen dahinter. Die Daten stammen aus dem Datenset mtcars. Der Summary-Befehl gibt uns die Ausgabe relevanten Kennwerte des Modells.
# Zusammenfassung der Kennwerte
summary(Verbrauch)
Diese Ausgabe gibt uns nun den F-test für das Gesamtmodell an. Dieser gibt uns einen p-Wert von 9.109e-12 an. D.h. er ist kleiner als 0,05. Die Varianzaufklärung wird uns durch das adjustierte Rquadrat gegeben und beträgt 0,8148. Die beiden Koeefizenten hp und wt scheinen in unserem Modell eine Vorhersagekraft zuhaben, da ihr p-Wert ebenfalls signifikant geworden ist. Die einzelnen Koeffizenten sind aber nicht standartisiert und lassen sich so nicht vergleichen. Dazu brauchen wir die standartisiereten Regressionskoeffizenten. Diese finden sich nicht im R-Commander sondern können nur durch einen Befehl im Paket QuantPsyc erstellt werden. Hierzu laden wir erst dieses Paket und führen den Befehl (lm.beta(Model)) aus.
# laden des Paketes
library(QuantPsyc)
# stadartisierte Koeffizenten
lm.beta(Verbrauch)
Nun können wir erkennen, dass der standartisierte Betakoeffizent des Gewichts größer ist als der der PS Anzahl.
Um die Teststärke in Gpower für das Model zu berechnen, müssen wir in Gpower verschiedene Einstellungen tätigen. Diese kann man alle im folgenden Screenshot erkennen. Wichtig ist, dass wir die Effektstärke in f quadrat umrechnen. Dies geht auf dem mittleren linken Knopf determine. Dann öffnet sich rechts ein Fenster, wo wir unserer Adjustiertes R Quadrat einsetzen. Wir halten dann einen fquadratwert, welcher sich automatisch übertragen lässt.Nun tragen wir die anderen Parameter ein: Alpha 0,05, N= 32. Dadurch erhalten wir eine Teststärke von 1. Wir haben also ausreichend Power. Siehe auch Bild: 
Zuerst berechnen wir die standardisierten und studentisierten Residuen.
# Studentisierte Residuen:
rstudent(Verbrauch)
# standartisierte Residuen:
rstandard(Verbrauch)
Wir können jetzt die standartisierten und studentisierten Residuen betrachten und uns die Abweichungen genauer anschauen. Es fällt auf das drei Werte relativ hoch liegen. Wir können an dieser Stelle die Werte der studentisierten Residuen direkt in unserer Datenmatrix abspeichern. Dazu gehen wir im Commander auf Modelle/Add observation statistics to data
  Hier können wir nun die Statistiken auswählen, die uns R mit in die Datamatrix editieren soll. Jetzt werden schon alle ausgeführt, wir werden uns aber erstmal auf die Residuen beziehen. Die generierten Befehle dafür lauten:
# Fitted sind die Vorhergesagten Werte
mtcars$fitted.Verbrauch <- fitted(Verbrauch)
# Das sind die Residuen
mtcars$residuals.Verbrauch <- residuals(Verbrauch)
# Studentisierte Residuen
mtcars$rstudent.Verbrauch <- rstudent(Verbrauch)
# hat values
mtcars$hatvalues.Verbrauch <- hatvalues(Verbrauch)
# cooks.distance
mtcars$cooks.distance.Verbrauch <- cooks.distance(Verbrauch)
mtcars$obsNumber <- 1:nrow(mtcars)
Was bei dem Befehl auffält ist, dass der R-Commander uns nur die studentisierten und nicht die standartisierten Residuen gibt. Die studentisierten sind bei kleineren Stichproben bessere Schätzungen. Falls wir uns aber auch noch die standartisierten einspeichern wollten, können wir uns den Befehl nach dem selben Schema selbst erstellen.
# Der Befehl ist genauso, wie oben. Nur die Funktion und der Name sind von rstudent in rstandard geändert
mtcars$rstandard.Verbrauch <- rstandard(Verbrauch)
Wir können uns die Werte in der Datenmatrix anschauen.  Bei den Studentisierten Residuen sollten nur 5% größer als 2 (bzw. kleiner als -2) sein und nur 1% größer als 2,5 ( bzw. kleiner als 2,5) Um dies zu testen könnten wir rausfinden, wie viele studentisierte Residuen in diesen Bereich fallen. Das würde in unserem Fall noch mit Auszählen gehen, aber R bietet uns dort Hilfe an. Zuerst lassen wir uns alle Fälle ausgeben, die zwischen -2 und 2 liegen. Dies geht mit dem folgenden Befehl:
# Der Befehl gibt uns an, welche Werte zwischen 2 und -2 liegen
mtcars$rstudent.Verbrauch > 2| mtcars$rstudent.Verbrauch < -2
Wir erhalten als Output TRUE oder False. True bedeutet, dass der Wert größer beziehungsweise kleiner ist. Wenn wir die Funktion ebenfalls abspeichern können wir, über die Summenfunktion direkt die Anzahl erfahren.
# Zuordnung einer Funktion für die Datenmatrix
mtcars$groß.studentisierteresiduen2<-mtcars$rstudent.Verbrauch > 2| mtcars$rstudent.Verbrauch < -2
# Summenbefehl
sum(mtcars$groß.studentisierteresiduen2)
Wir haben also 3 Residuen, die größer sind als 2. Um zu schauen, ob sie auch größer sind als 2,5 führen wir den gleiche Befehl für den Bereich zwischen -2,5 und 2,5 aus.
# Zuordnung einer Funktion für die Datenmatrix, Funktionsname hat sich geändert.
mtcars$groß.studentisierteresiduen25<-mtcars$rstudent.Verbrauch > 2.5| mtcars$rstudent.Verbrauch < -2.5
# Summenbefehl
sum(mtcars$groß.studentisierteresiduen2)
Es fallen immer noch 3 Werte aus dem Rahmen. Diese werden wir weiterhin beobachten.
Wir können die studentisierten Residuen auch graphisch betrachten. Dazu haben wir mehrere Möglichkeiten. Wir können die Studentisierten Residuen in einem Streudiagramm gegen die vorhergesagten Werte plotten. Dies geht auf dieselbe Art und Weise, wie oben. ![Streudiagrammfittedvsresidualoptionen] (C:/Users/Matthias/Desktop/R Daten/Markdownprojekt/Screenshots/Streudiagrammfittedvsresidualoptionen.png) ![Streudiagrammfittedvsresidualoptionen2] (C:/Users/Matthias/Desktop/R Daten/Markdownprojekt/Screenshots/Streudiagrammfittedvsresidualoptionen2.png)
scatterplot(rstudent.Verbrauch~fitted.Verbrauch, reg.line=lm, smooth=FALSE,
spread=FALSE, id.method='mahal', id.n = 2, boxplots=FALSE, span=0.5,
xlab="vorhergesagte Werte", ylab="studentisierte Residuen",
main="studentisierte Residuen vs vorhergesagte Werte", data=mtcars)
Auch könnte man sich Histogramm für die Daten zeichnen lassen:   
# Das ist der Befehl fürs Histogramm
with(mtcars, Hist(rstudent.Verbrauch, scale="percent", breaks=10,
col="darkgray"))
Eine andere Möglichkeit diesen Graphen zu erzeugen findet man unter Modelle/Grafiken/Grundlegende diagnostische Grafiken
oldpar <- par(oma=c(0,0,3,0), mfrow=c(2,2))
plot(Verbrauch)
par(oldpar)
Dafür verwenden wir den Breusch-Pagan Test. Er befindet sich unter Modelle/Numerical diganostics:   Es wird folgender Befehl generiert:
# Die Funktion lädt sich automatisch ein neues Paket
library(lmtest, pos=4)
#Das hier ist er eigentliche Befehl:
bptest(mpg ~ hp + wt, varformula = ~ fitted.values(Verbrauch),
studentize=TRUE, data=mtcars)
Der Test wird signifikant, d.h. es liegt keine Homoskeadastität vor Jetzt müssten wir eigentlich Maßnahmen ergreifen, um dem entgegen zu wirken.
Die hat values haben wir bereits in R erzeugt. Wir finden sie in unserer Datenmatrix.    Der Einfluss von einer hat value ist zu groß, wenn hii > 2(k+1) / n. Wir haben Unabhängige Variablen und 32 Datenpunkte. Daraus folgt:
# Man kann direkt in R den Wert berechnen
2*(2+1)/32
Wir könne mit Hilfe von R wieder prüfen, wie viele Werte drüber liegen dürfen. Dies geht wieder, indem wir eine neue Funktion bilden.
# Hier erstellen wir wieder eine Funktion
mtcars$groß.hatvalues<-mtcars$hatvalues.Verbrauch > .1875
# Hier lassen wir uns die Summe rausgeben
sum(mtcars$groß.hatvalues)
Es gibt also 3 Datenpunkte, die die kritische hatvalue von 0.1875 überschreiten. Wir können uns das ganze auch grapfisch mit Hilfe eines Histogramms anschauen. ( Wie man zum Histogrammcode kommt, kann man oben nochmal nachschauen)
with(mtcars, Hist(hatvalues.Verbrauch, scale="frequency", breaks=10,
col="darkgray"))
Man kann erkennen, das ein Wert über eine sehr große hat value verfügt.Man müsste diesen Datenpunkt ausschließen.
Cooks distance wird uns wieder direkt in R berechnet und uns in die Datenmatrix eingeführt. Uns interessieren Werte, welche über 1 liegen. Werte über 4 müssen ausgeschlossen werden. Wir bilden wieder die übliche Funktion, um die Werte zu identifizieren.
# Hier erstellen wir wieder eine Funktion
mtcars$groß.cdistance1<-mtcars$cooks.distance.Verbrauch > 1
# Hier lassen wir uns die Summe rausgeben
sum(mtcars$groß.cdistance1)
# Hier erstellen wir wieder eine Funktion
mtcars$groß.cdistance4<-mtcars$cooks.distance.Verbrauch > 4
# Hier lassen wir uns die Summe rausgeben
sum(mtcars$groß.cdistance4)
Die Anzahl der Werte, die über 1 liegen ist 0. Deshalb musste die Anzahl der Werte, welche über 4 liegen auch 0 sein. Wir können uns dies auch nochmal graphisch mit Hilfe eines Histogramms anschauen.
with(mtcars, Hist(cooks.distance.Verbrauch, scale="frequency",
breaks="Sturges", col="darkgray"))
Man kann erkennen, dass die einzelnen Werte noch sehr weit von der 1 entfernt sind.
Wir testen nun noch die Unabhängigkeit der Fehler. Dazu verwenden wir den Durbin-Watson Test. Diesen findet man unter Modelle/Numerical Diagnostic   Wir wählen Rho = 0 aus, da wir die Hypothese testen, ob eine Autokorrelation vorliegt. Es wird folgender Befehl generiert:
dwtest(mpg ~ hp + wt, alternative="two.sided", data=mtcars)
Die Fehler sind in unseren Daten also nicht unabhängig voneinander.
Um die Multikollinearität zu analysieren berechnen wir die Variance-Inflation-Faktoren (VIFs) für unseren beiden Unabhängigen Variablen. R berechnet uns diese direkt. Modelle/numerical diagnostics/Variance-Inflation-Factor 
vif(Verbrauch)
Die Werte liegen nicht über 10. Wir können sie miteinnehmen.