9 Lineare Regression in:

Andreas Behr, Ulrich Pötter

Einführung in die Statistik mit R, page 134 - 165

2. Edition 2010, ISBN print: 978-3-8006-3599-3, ISBN online: 978-3-8006-4878-8, https://doi.org/10.15358/9783800648788_134

Series: Vahlens Kurzlehrbücher

Bibliographic information
L R Wohl die wichtigsten Aufgabenstellungen der Statistik ergeben sich aus der Frage nach dem Zusammenhang zwischen mehreren statistischen Variablen. Die lineare Regression ist vermutlich das am häufigsten benutzte statistische Verfahren für diesen Zweck. In diesem Kapitel beschreiben wir seine Grundlagen sowie die Umsetzung in R. 9.1 Grundlagen 9.2 Lineare Modelle in R 9.2.1 Der lm Befehl 9.2.2 Formeln für Regressionen 9.2.3 Modellvergleiche und fehlende Werte 9.2.4 Diagnostik 9.2.5 Ausschluss von Ausreißern 9.2.6 Transformation der abhängigen Variablen 9.2.7 Transformationen der Kovariablen 9.3 Regression in Matrixnotation 9.3.1 Einfache Matrixoperationen in R 9.3.2 Das Regressionsmodell in Matrixnotation 9.3.3 Die Berechnung von Regressionen in R 9.4 Übungsaufgaben . Grundlagen Der Zusammenhang mehrerer statistischer Variablen lässt sich nur schwer mit den Methoden des Kapitels darstellen, zumindest wenn es um mehr als zwei Variable geht. Zwar kann man an Hand der gemeinsamen Verteilungen alle Fragen über den Zusammenhang zwischen Variablen beantworten. Nur sind die gemeinsamen Verteilungen weder einfach zu beschreiben noch einfach zu schätzen oder zu simulieren. Die Grundidee der linearen Regression besteht nun darin, an Stelle der gemeinsamen Verteilung besonders einfache Charakterisierungen der bedingten Verteilung einer Variablen zu betrachten. Insbesondere 5 10 35 83 _W iS o B eh r P öt te r 2A - B g 5 Lineare Regression betrachtet man die bedingten Mittelwerte einer statistischen Variablen Y für alle Werte der bedingenden Variablen X = x ,X = x , . . . ,Xm = xm als lineare Funktion dieser Werte. Für die durch die Werte statistischer Variabler x , x , . . . , xm definierten Gruppen soll für deren Mittelwert approximativ gelten: M(Y X = x ,X = x , . . . ,Xm = xm) b + b x + b x + . . . + bmxm Abbildung 9.1: Lineares stochastisches Modell. Der Einfluss einer Variablen lässt sich dann besonders einfach zusammenfassen: Eine Änderung des Wertes von Xk = xk zu Xk = xk ändert den bedingten Mittelwert um bk(xk xk). Der Einfluss jeder der bedingenden Variablen kann also einfach am zugehörigen Regressionskoeffizienten abgelesen werden. Und ein Vergleich von zwei Gruppen, die durch x , x , . . . , xm bzw. x , x , . . . , xm gekennzeichnet sind, ergibt einen ungefähren Unterschied von b (x x ) + . . . + bm(xm xm) in den Mittelwerten der Gruppen. Im Standardmodell der Statistik lässt sich diese Situation ebenso einfach beschreiben und simulieren. Setzt man für ZufallsvariableX, є, die auf einem gemeinsamenWahrscheinlichkeitsraum definiert sind, Y := b + b X + є (9.1) wobei E(є X = x) = für alle x gelten soll, dann ist offenbar E(Y X = x) = b + b x Wir werden diese Gleichung ein stochastisches lineares Modell nennen. Als nächsten Schritt kann man Daten dieser Art simulieren: Denn auf den ersten Blick muss man nur die Verteilung vonX und die von є angeben und erhält mit entsprechend verteilten Pseudozufallszahlen für alle vorgegebenen Parameter b und b neue Werte von Y, indem einfach die entsprechende Gleichung ausgerechnet wird. Wir können dieses Vorgehen mit den Mitteln des letzten Kapitels direkt in ein R-Programm übersetzen: In stochastischen linearen Modellen betrachtet man die Parameter b und b (und Eigenschaften von є) als gegeben. Man kann sie variieren, aber dann spricht man von verschiedenen stochastischenModellen. DasModell ist linear inX. Das steht imGegensatz zumSprachgebrauch in der Statistik. In statistischen linearen Modellen interessiert man sich für die („unbekannten“) Parameter b und b und nennt ein Modell linear, wenn es linear in diesen Parametern ist. Ob als statistische Variable X oder beliebige Funktionen f (X) erscheinen, ist dabei unerheblich. Das ist offenbar auch die gewünschte Interpretation statistischer Analysen vorgelegter Populationsdaten, in denen man sich für die (approximativen) Unterschiede der Mittelwerte in verschiedenen Subgruppen der Gesamtheit interessiert. . Lineare Modelle in R > set.seed(54) > n < 20 > x < rnorm(n) > epsilon < rnorm(n) > y < 1+2∗x+epsilon > plot(x,y,bty="l") > abline( 1,2) Hier werden also Werte der ZufallsvariablenY erzeugt, die linear von den Werten von X abhängen, zu denen jeweils ein unabhängiger Fehlerterm є addiert wurde. Diese Realisation des linearen Modells ist in Abbildung . dargestellt. Die Linie ist die lineare Funktion y = f (x) = + x im Intervall [ . , ]. Nun müssen Daten und die zugehörigen statistischen Variablen natürlich nicht den Vorgaben eines recht willkürlich gewählten stochastischen linearen Modells ähneln. Aber die Verwendung der linearen Regression setzt das auch gar nicht voraus, weder im Rahmen des stochastischen linearen Modells noch im Rahmen der reinen Datenbeschreibung. Man benutzt einfach die Regressionskoeffizienten, die die bedingte Verteilung von Y gegeben X am besten beschreiben. Dabei soll „am besten“ heißen, dass der durchschnittliche Abstand zwischen den Y(u) und denWerten von b + b X (u) + . . . + bmXm(u) der kleinst mögliche ist. Die lineare Regression ist daher definiert durch die Werte der Regressionskoeffizienten b , b , . . . , bm, die die Werte u (Y(u) b b X (u) b X (u) . . . bmXm(u)) (9.2) minimieren. Diese Werte, die man häufig mit b̂ , b̂ , . . . , b̂m bezeichnet, erzeugen die beste lineare Näherung an die Werte von Y(u), u durch eine Linearkombination der Xi(u), u , i , . . . ,m . Die stochastischen linearen Modelle dienen nur dazu, einen der vielen möglichen Vergleichsmaßstäbe zu konstruieren, die es erlauben, die Eigenschaften des Verfahrens zu beurteilen. Es ist gerade die Fähigkeit der linearen Regression, beliebige Zusammenhänge zumindest imDurchschnitt optimal darzustellen, die seine häufige Verwendung rechtfertigt. . Lineare Modelle in R . . Der lm Befehl In R steht mit der Funktion lm() (wie linearmodel) eine komfortable und flexible Funktion zur Schätzung linearer Modelle bereit. Als Argument ist Die in der Simulation verwandte Annahme der Unabhängigkeit vonX und є ist eine viel stärkere Annahme als die oben benutzte Form, die nur die bedingten Erwartungswerte von є beschränkt, insbesondere also verschiedene Verteilungen von є ebenso wie Abhängigkeiten zwischenX und є zulässt. Lineare Regression der Zusammenhang zwischen unabhängigen Variablen und einer abhängigen Variablen in Form einer Formel anzugeben. Die abhängige Variable Y und die Kombination der unabhängigen Variablen Xi wird in symbolischer Form angegeben. Die einfachste Form einer Formel in R ist y~x1+x2. Die Ergebnisse einer linearen Regression der Variablen y auf die zwei unabhängigen Variablen x1 und x2 erhält man entsprechend durch den Befehl > lm(y~x1+x2) Wir berechnen die beste lineare Regression der Daten, die wir im letzten Abschnitt simuliert haben: > lm(y~x) Call: lm(formula=y~x) Coeicients: (Intercept) x 0.910 2.028 Das Resultat des Befehls lm() ist eine Liste der Ergebnisse einer linearen Regression. Die Liste enthält insbesondere die Koeffizienten b̂ , b̂ , . . . , b̂m, den Aufruf des Befehls lm() (also insbesondere die benutzte Formel) und weitere Elemente, die es erlauben, das vollständige Ergebnis der Regression zu rekonstruieren. Der direkte Aufruf der Funktion lm() liefert nur die spezielle Ausgabedarstellung der Funktion auf dem Bildschirm, also weder alle Elemente der Ergebnisliste noch die üblichen Angaben zu den Ergebnissen, die man vielleicht erwarten würde. Denn in R hat jedes Objekt abhängig von seiner Klassenzugehörigkeit eine Darstellungsform, die gewählt wird, wenn nur der Name des Objektes an R geschickt wird. Die Angabe der Funktion lm() allein liefert daher nur eingeschränkte erste Informationen. Weist man den Ergebnissen aber einen Namen zu, kann man auf alle berechneten Ergebnisse zugreifen. > erg < lm(y~x) > is.list(erg) [1] TRUE > class(erg) [1] "lm" Das Ergebnis des lm() Befehls ist also in der Tat eine Liste der Klasse lm. Insbesondere sind die geschätzten Koeffizienten b̂ , b̂ durch erg$coeicients ansprechbar: > erg$coeicients (Intercept) x 0.9099756 2.0276285 Statt „unabhängige Variable“ werden auch oft die Begriffe „Kovariable“, „Regressor“ oder „Prädiktor“ benutzt. Wir verwenden alle diese Begriffe synonym. . Lineare Modelle in R Sowohl für die Regressionskoeffizienten wie für deren geschätzte Varianz- Kovarianzmatrix gibt es Funktionen, die diese Informationen aus dem Objekt erg extrahieren. Für die Koeffizienten ist dies die Funktion coef() > coef(erg) (Intercept) x 0.9099756 2.0276285 und für die Kovarianzmatrix vcov() > vcov(erg) (Intercept) x (Intercept) 0.05582266 0.00912990 x 0.00912990 0.04262705 Man kann also recht einfach auf Teilergebnisse einer Regression zugreifen, um damit weiterzurechnen. Nun wäre es sehr umständlich, Ergebnisse einer Regression durch Angabe der entsprechenden Elemente der Liste erg einzeln aufzurufen oder nur die Funktionen coef() bzw. vcov() benutzen zu können. Der Befehl summary() liefert einen ersten Überblick über die Ergebnisse, in dem die wichtigsten Maßzahlen in einer übersichtlichen Tabelle zusammengestellt werden: > summary(erg) Call: lm(formula=y~x) Residuals: Min 1Q Median 3Q Max 1.9619 0.4675 0.3422 0.5596 1.7522 Coeicients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.9100 0.2363 3.851 0.00117 ∗∗ x 2.0276 0.2065 9.821 1.18e 08 ∗∗∗ Residual standard error: 1.038 on 18 degrees of freedom Multiple R squared: 0.8427, Adjusted R squared: 0.834 F statistic: 96.45 on 1 and 18 DF, p value: 1.179e 08 Das Ergebnis des summary() Befehls ist ebenfalls eine Liste. Man kann daher einige zusätzliche Angaben aus dieser Liste zu weiteren Berechnungen benutzen: > erg2 < summary(erg) > names(erg2) [1] "call" "terms" "residuals" "coeicients" [5] "aliased" "sigma" "df" "r.squared" [9] "adj.r.squared" "fstatistic" "cov.unscaled" > erg2$sigma ##geschaetzte Varianz der Fehlerterme [1] 1.037952 Namen können wieder abgekürzt werden, solange die Abkürzung eindeutig ist: Lineare Regression > erg2$f #F Statistik value numdf dendf 96.4476 1.0000 18.0000 Besonders interessant ist das Element erg2$coeicients des summary() Objektes. Es enthält nicht nur die geschätzten Koeffizienten, sondern auch die geschätzten Standardfehler, die zugehörige t-Statistik und den beobachteten p-Wert für jede Kovariable. Das Objekt erg2$coeicients ist als Matrix organisiert, wobei die erste Spalte die geschätzten Koeffizienten, die zweite Spalte die Standardfehler, die dritte die t-Werte und die vierte Spalte die beobachteten p-Werte enthält. Will man etwa alle beobachteten Signifikanzniveaus der unabhängigen Variablen erhalten, dann kann man schreiben: > erg2$coef[,4] #abgekürzter Name (Intercept) x 1.169323e 03 1.178836e 08 . . Formeln für Regressionen Die symbolische Form, in der abhängige und unabhängige Variable im Befehl lm() (und in fast allen anderen Versionen von Regressionen, Tabellen etc.) angegeben werden, ist sehr flexibel. Wir illustrieren den Gebrauch an den Angaben im Mikrozensus zu den Miethöhen. Dazu lesen wir zunächst die Daten wieder ein und wählen einige interessierende Variablen aus: > library(foreign) > dat < read.spss("mz02_cf.sav",to.data.frame=T,use.value.labels=F) > hh < dat$ef3∗100+dat$ef4 #Nr des Auswahlbezirks∗100 + #+HHnr im Bezirk als HHId > oo < !duplicated(hh) #Nur eine Angabe/HH Wir betrachten nur die acht Variablen Miethöhe ef462, Wohnungsgröße ef453, Anzahl der Personen in derWohnung ef500 und Anzahl der Erwerbstätigen im Haushalt ef522, sowie die Anzahl der Kinder imHaushalt, die in den Variablen ef528 bis ef531 für Kinder unter Jahren, Kinder im Alter von bis unter , bis unter und bis unter abgelegt sind. Da es in einer Wohnung mehrere Haushalte geben kann, beschränken wir uns auf Wohnungen mit einem Haushalt und schließen auch Gemeinschaftsunterkünfte aus. Außerdem beschränken wir uns zunächst auf die Fälle mit vollständigen Angaben. Die Datenauswahl ist dann: > Wohn < cbind(dat[oo,"ef462"],dat[oo,"ef453"],dat[oo,"ef500"], + dat[oo,"ef522"],dat[oo,"ef22"],dat[oo,"ef528"], + dat[oo,"ef529"],dat[oo,"ef530"],dat[oo,"ef531"]) > Wohn < Wohn[complete.cases(Wohn),] > Wohn < subset(Wohn,Wohn[,1]>0&Wohn[,1]<=1800& + Wohn[,2]>=10&Wohn[,2]<300&Wohn[,3]>0& + Wohn[,4]<9&Wohn[,5]==1& . Lineare Modelle in R + Wohn[,6]<9&Wohn[,7]<9&Wohn[,8]<9& + Wohn[,9]<9 ) > Wohn < data.frame(Miete=Wohn[,1],qm=Wohn[,2], + HHGroesse=Wohn[,3], + ErwPerson=Wohn[,4], + AnzKu3=Wohn[,6],AnzK3b6=Wohn[,7], + AnzK6b10=Wohn[,8], + AnzK10b15=Wohn[,9]) > attach(Wohn) Wir erhalten valide Angaben zu allen Variablen für Haushalten. In Abschnitt . ist schon die gemeinsame Dichte von Miethöhe und Wohnungsgröße dargestellt worden. Mit dem Regressionsansatz können nun aber sehr leicht viel mehr Variable in die Analyse einbezogen werden. Wir wählen zunächst als abhängige Variable die Miethöhe und als unabhängige Variable Wohnungsgröße, Anzahl der Personen im Haushalt und die Zahl der erwerbstätigen Personen imHaushalt. Eine erste Analyse ergibt > erg1 < lm(Miete~qm+HHGroesse+ErwPerson) > summary(erg1) Call: lm(formula=Miete~qm+HHGroesse+ErwPerson) Residuals: Min 1Q Median 3Q Max 552.59 69.02 12.61 52.15 1012.23 Coeicients: Estimate Std. Error t value Pr(>|t|) (Intercept) 41.25132 4.77447 8.640 < 2e 16 ∗∗∗ qm 4.17014 0.07598 54.881 < 2e 16 ∗∗∗ HHGroesse 6.52986 1.81933 3.589 0.000335 ∗∗∗ ErwPerson 23.55631 2.22283 10.597 < 2e 16 ∗∗∗ Residual standard error: 122.1 on 5693 degrees of freedom Multiple R squared: 0.4431, Adjusted R squared: 0.4428 F statistic: 1510 on 3 and 5693 DF, p value: < 2.2e 16 Das Ergebnis ist auf den ersten Blick recht überzeugend. Alle Einflüsse der Kovariablen sind hoch signifikant, das R ist mit % weit besser als in vielen anderen Anwendungen. Aber man sollte vielleicht nicht einfach die Anzahl der Personen betrachten, sondern noch unterscheiden, ob es sich um Kinder oder Erwachsene handelt. Aber wie passt das in die Syntax der Formel? Die Anzahl der Kinder unter Jahren würde man in R durch AnzKu3+AnzK3b6+AnzK6b10+AnzK10b15 ausdrücken. Schreibt man das in den lm() Befehl, dann werden die vier Terme getrennt als unabhängige Variable betrachtet. Um zwischen der Syntax für Formeln und der von Rechenoperationen in R (hier: Addition) zu unterscheiden, kann die Funktion I() (ein Lineare Regression großes I) benutzt werden. Will man getrennt die Anzahl Erwachsener und die Anzahl der Kinder im Haushalt berücksichtigen, ohne zunächst neue Variable zu definieren, dann kann man schreiben: > erg2 < lm(Miete~qm+HHGroesse+I(AnzKu3+AnzK3b6+ + AnzK6b10+AnzK10b15)+ErwPerson) Das Ergebnis (leicht gekürzt) ist Coeicients: Estimate Std. Error t value Pr(>|t|) (Intercept) 45.43474 5.14999 8.822 < 2e 16 ∗∗∗ qm 4.17595 0.07601 54.941 < 2e 16 ∗∗∗ HHGroesse 10.74357 2.66529 4.031 5.63e 05 ∗∗∗ I(AnzKu3 + ... 8.05660 3.72516 2.163 0.0306 ∗ ErwPerson 24.96405 2.31549 10.781 < 2e 16 ∗∗∗ Residual standard error: 122.1 on 5692 degrees of freedom Multiple R squared: 0.4435, Adjusted R squared: 0.4431 F statistic: 1134 on 4 and 5692 DF, p value: < 2.2e 16 Einige weitere Symbole, die im Rahmen einer Formel in lm() (und anderen Regressionsmodellen) eine andere Bedeutung als im Rest von R haben, sind neben + noch ^, :, , ∗, . und %in%. Will man die ursprüngliche Bedeutung dieser Symbole in einer Formel benutzen, muss immer die Form I() benutzt werden. UmMissverständnissen vorzubeugen, ist es ratsam, für Rechenoperationen in Formel-Ausdrücken immer die I() Form zu verwenden, zumal da in einigen Zusatzpaketen Erweiterungen von Formel-Ausdrücken („specials“) definiert sind, die die Bedeutung von Funktionsaufrufen in Formeln verändern können. Die Basiselemente der symbolischen Formeln in R sind+ für das Hinzufügen eines Regressors und für dessen Entfernung. Insbesondere wird ohne weitere Angaben immer eine Konstante (ein Achsenabschnitt) eingeschlossen. Soll aber eine homogene Regressionsfunktion geschätzt werden, kann die automatische Aufnahme durch - auf der rechten Seite des ~ Zeichens unterdrückt werden. Man muss also für eine homogene Regression durch den Nullpunkt von y auf x1 schreiben: y ~ 1+x1. Sind x1 und x2 Kovariable, so werden Interaktionen symbolisch durch x1:x2 angedeutet. Ein Regressionsmodell mit zwei Kovariablen x1 und x2 und deren Interaktion lässt sich also in Rs Formel als y~x1+x2+x1:x2 schreiben. Eine Abkürzung ist (x1+x2)^2. Interaktionen höherer Ordnung einschließlich aller Terme niedrigerer Ordnung erhält man etwa durch (x1+x2+x3)^3. Das liefert neben den Haupteffekten alle Interaktionen der drei Kovariablen einschließlich des Interaktionsterms für alle drei Variablen. Einzelne Interaktionen kann man wieder durch das Symbol ausschließen. So schließt etwa -x1:x3 die Interaktion zwischen x1 und x3 aus. . Lineare Modelle in R Ist man sich über die durch einen Formel-Ausdruck erzeugte Formel nicht sicher (oder möchte diese Information in einem Programm weiter verwenden), dann kann der Befehl terms() benutzt werden. Das erzeugte Objekt enthält weitere Informationen über die Formel, die in den Attributen des Objektes enthalten sind. So liefert > ttt < terms(y ~ (x1+x2+x3)^3 x1:x2) > attr(ttt,"term.labels") #Attribut "term.labels" [1] "x1" "x2" "x3" "x1:x3" "x2:x3" "x1:x2:x3" die symbolische Form aller Terme im Formel-Ausdruck, hier also die drei Haupteffekte x1,x2,x3, die Interaktionen zwischen x3 und x1 bzw. x2 (die Interaktion zwischen x1 und x2 fehlt), sowie die Interaktion aller drei Variablen. . . Modellvergleiche und fehlende Werte Wir haben bei der Analyse der Mikrozensusdaten bisher fehlende Werte durch entsprechende Befehle explizit ausgeschlossen. Man wird immer so vorgehen, wenn entweder von vornhinein sicher ist, welche Variablen man in allen folgenden Analysen verwenden möchte oder wenn man systematische Modellvergleiche und Verfahren der Modellwahl benutzen möchte. Denn Modellvergleiche setzen i.d.R. voraus, dass sich alle betrachtetenModelle auf die gleiche Datenbasis beziehen. Z.B. kann man den zusätzlichen Term im zweiten Modell für die Miethöhe auch durch einen Test zwischen beiden Modellen direkt beurteilen wollen (in unserem Beispiel mit nur einem zusätzlichen Term ohne Interaktionen ist ein solcher Test äquivalent zu dem entsprechenden Test des zusätzlichen Terms in der Regression). Modelle können durch den Befehl anova() verglichen werden, der als Argumente eine Folge von Modellergebnissen gestattet. Ein Vergleich von erg1 und erg2 ergibt sich durch den Aufruf > anova(erg1,erg2) Model 1: Miete ~ qm+HHGroesse+ErwPerson Model 2: Miete ~ qm+HHGroesse+I(AnzKu3+AnzK3b6+ AnzK6b10+AnzK10b15)+ErwPerson Res.Df RSS Df Sum of Sq F Pr(>F) 1 5693 84862030 2 5692 84792350 1 69680 4.6775 0.0306 ∗ der einen einfachen F-Test zwischen den Modellen berechnet. Um eine einheitliche Auswahl von Fällen zu erreichen, die nur gültige Werte für alle interessierenden Variablen enthalten, kann man den Befehl complete.cases() verwenden. Wir haben ihn schon mehrfach benutzt. Der Befehl erwartet als Argumente eine Folge von Vektoren, Matrizen oder Dataframes und gibt einen logischen Vektor zurück, der für diejenigen Zeilen der Argumente den Wert TRUE enthält, in denen keine der Variablen in den Argumenten einen definierten fehlenden Wert NA aufweist. Lineare Regression In den meisten Fällen aber möchte man in verschiedenen Modellen alle Fälle benutzen, die für die gerade betrachteten Variablen Informationen enthalten. Die verschiedenen implementierten Regressionsmodelle und insbesondere der lm() Befehl schließen automatisch alle Fälle aus, die auf einer der im Formel-Ausdruck angegebenen Variablen einen fehlenden Wert ausweisen. Die Information über die tatsächliche Anzahl der Fälle, auf denen eine Analyse beruht, kann man in den summary() Methoden nur indirekt über die Anzahl der Freiheitsgrade (df) bzw. über die in vielen Modellen ausgegebene Anzahl der ausgeschlossenen Fälle erschließen. Diese Information ist aber extrem wichtig für die Interpretation von Ergebnissen. Sie sollte gerade bei der Modellentwicklung immer zuerst gesucht werden! In vielen Fällen ergeben sich extrem gute Ergebnisse, die aber nur auf einem Bruchteil der Daten basieren. Der Befehl lm() hat ein Argument na.action, mit dem die Behandlung fehlender Werte in der Berechnung einer Regression beeinflusst werden kann. Es gibt die Möglichkeiten na.fail, das eine Fehlermeldung produziert, wenn eine der Variablen einen fehlenden Wert (NA) enthält, na.omit, das alle Fälle mit fehlenden Werten ausschließt, sowie na.exclude mit gleichem Verhalten wie na.omit, aber mit dem Unterschied, dass Teilergebnisse (etwa Vektoren von Residuen) so ergänzt werden, dass sie mit den Fällen der Ausgangsdaten übereinstimmen. Diese Version ist fast immer vorzuziehen, wenn Regressionsergebnisse mit verschiedenen Varianten von Kovariablen zumindest informell verglichen werden sollen. Die Voreinstellung dieser Möglichkeiten richtet sich nach den globalen Optionen, die mit options("na.action") abgefragt werden können. Die Voreinstellung ist na.omit, so dass alle Zeilen mit fehlenden Werten in wenigsten einer Variablen ausgeschlossen werden. Will man die Option nicht in allen Aufrufen der Regressionen jeweils explizit ändern, kann man sie auch global durch options(na.action="na.exclude") setzen. . . Diagnostik In fast allen Anwendungen linearer Regressionen kann man sich bei der Beurteilung der Ergebnisse nicht allein auf die Ergebnisse stützen, die der summary() Befehl zur Verfügung stellt. Es bleiben Fragen nach den Gründen für besonders gute oder schlechte Modellanpassung und die Möglichkeit von einflussreichen Ausreißernmuss untersucht werden. Zudem kann ein genauerer Blick auf die Einzelergebnisse der Regression Hinweise für eine Revision des ursprünglichen Modells liefern, die zu deutlich besser interpretierbaren und angemesseneren Ergebnissen führt. In R ist es fast immer hilfreich, sich einen ersten Überblick über Aspekte ungenügenden Modellfits durch den Befehl plot(erg1) zu verschaffen. Hat man das Ergebnis einer linearen Regression Da das Modell erg1 zur Klasse lm gehört, wird durch den plot() Befehl die speziellere Graphikfunktion plot.lm() ausgeführt. Sie ist im Paket stats dokumentiert, nicht im Paket graphics. Der Grund ist, dass jedes Paket eigene Varianten der Befehle plot(), print(), summary() etc. für neue, in diesem Paket definierte Klassen bereitstellen kann. . Lineare Modelle in R (wie erg1 für die Miethöhen imMikrozensus), dann erzeugt der Befehl vier Graphiken mit a) den standardisierten Residuen aller Fälle aufgetragen gegen die durch die Regression vorhergesagten Werte, b) die Wurzeln des Absolutbetrags der Residuen gegen die vorhergesagten Werte, c) eine Q-Q-Graphik der Verteilung der Residuen im Vergleich zu einer Normalverteilung und d) eine Graphik mit den Residuen gegen die Hebelpunkte (Leverages), die allein durch die relative Lage der Kovariablen definiert sind. In Abbildung . sind die Graphiken a) und b) dargestellt (wir haben in den Graphiken viele Punkte ausgeschlossen, um etwas übersichtlichere Druckdarstellungen zu schaffen). (a) (b) Abbildung 9.2: Diagnostische Graphiken. a) Standardisierte Residuen gegen angepasste Werte, b) Wurzel des Absolutbetrags der standardisierten Residuen gegen angepasste Werte. Eingezeichnet ist auch ein lokaler Glätter. Die in beiden Graphiken verwandten standardisierten Residuen werden aus demVektor derResiduen ê(u) := Y(u) b̂ b̂ X (u) . . . b̂mXm(u) durch eine Orientierung an einer speziellen Form des stochastischen Regressionsmodells (9.1) im Rahmen des klassischen Standardmodells der Statistik abgeleitet. Denn wenn man n unabhängige Wiederholungen des Modells (9.1) unterstellt und zusätzlich annimmt, dass alle n Zufallsvariablen є gleiche Varianz unabhängig vomWert der ZufallsvariablenX haben, dann hat auch die Zufallsvariable є̂ := Y b̂ b̂ X . . . b̂mXm eine Varianz. Man muss aber berücksichtigen, dass diese Größen nun nicht mehr unabhängig sein können (ihre Summe ist z.B. immer , denn sie müssen die Bedingung (9.2) erfüllen), noch sind sie identisch verteilt. Man kann dennoch für jede Replikation im klassischen Modell die Varianz von є̂ berechnen und erhält für gegebene statistische Variable X(u) = x V(є̂(u) X(u) = x) = V(є) ( huu) Der zusätzliche Term huu ist gerade der zu u gehörende Hebelpunkt. Er beschreibt den Einfluss einer Realisation der KovariablenX auf die Modellanpassung. Je größer huu ist, desto weiter ist die u-te Realisation vonX von ihrem Erwartungswert entfernt. Desto kleiner aber ist auch die Varianz der zugehörigen Residuen, denn solche weit entfernten Punkte bestimmen wesentlich die Lineare Regression Regressionskoeffizienten b̂ , b̂ , . . .. Das ist ihreHebelwirkung: Je größer huu ist, desto größer ist der Einfluss dieser Elemente u bei der Bestimmung der Lösung der Gleichung (9.2). Der Mittelwert der huu ist immer gerade (m + )/n, wobei m die Anzahl der Kovariablen in der Regression (ohne die Konstante) ist und n die Anzahl der Fälle bezeichnet. Das kann alsHinweis für sehr einflussreiche Kovariablenkonstellationen benutzt werden. Der Maximalwert der huu ist . Wenn tatsächlich huu = für ein u gilt, dann ist die Varianz dieses Residuums , der zugehörige Punkt (Y(u),X (u), . . . ,Xm(u)) liegt also unabhängig von demWert von Y(u) immer auf der Regressionsebene. Verschiebt man also den Wert von Y(u), dann verändert sich die gesamte Regressionsebene so, dass auch der neue Wert auf der Regressionsebene liegt. Die standardisierten Residuen definiert man nun durch ẽ(u) := ê(u)/ σ̂ huu wobei σ̂ ein Schätzer der Varianz der є ist, so dass die ẽ(u) nun annähernd gleiche Varianz haben. Man muss die vorhergehende Formel, die sich nur im Standardmodell der Statistik rechtfertigen lässt, nicht akzeptieren, um die Nützlichkeit der Standardisierung der Residuen einzusehen: Sie berücksichtigt den Einfluss einzelner Werte X(u) auf das Ergebnis einer Regression, ganz unabhängig von spezifischen Verteilungsannahmen über die KovariablenX im stochastischen Modell. Die Graphik . a) zeigt die standardisierten Residuen aufgetragen gegen die angepassten Werte, also die Werte ŷ(u) := b̂ + b̂ X (u) + . . . + b̂mXm(u). Die Graphik deutet bereits auf etliche Probleme mit unserer ersten Regression hin: Es gibt etliche sehr große (standardisierte) Residuen, die einen erheblichen Einfluss auf die Regressionskoeffizienten haben können. Insbesondere im Bereich positiver Residuen sindWerte bis zumMaximumvon . zu beobachten. Zudem wächst die Streuung der Residuen deutlich mit den angepasstenWerten der Regression an. Diesen letzten Effekt kann man einfacher in Graphik . b) beurteilen. Man möchte ja gern einen Eindruck davon gewinnen, wie groß die Varianz der Residuen in der Umgebung verschiedener angepasster Werte ist. Dazu benutzt man eine Transformation, die Wurzel der Absolutbeträge der standardisierten Residuen. Die Graphik enthält auch einen lokalen Durchschnitt der Größen. In der Literatur spricht man oft auch von erwarteten Werten, in Anspielung auf die Modellformulierung E(Y X , . . . ,Xm) = b + b X + . . . + bmXm im Rahmen des klassischen statistischen Modells. Die Sprechweise ist zwar verständlich, aber sie unterstützt die Tendenz, mathematische Modelle mit Zusammenfassungen von Daten zu verwechseln. Deshalb wählen wir zumeist die Bezeichnung „angepasste Werte“. Die naheliegendere Idee, die Quadrate der Residuen zu benutzen, um die lokale Varianz der Residuen zu approximieren, funktioniert nicht besonders gut, weil die Verteilung dieser Größen (unter einem klassischen Modell) sehr schief wäre und man daher keine vernünftige Vorstellung entwickeln kann, wie eine solche Graphik aussehen muss, wenn tatsächlich Daten unter einem mathematischen Modell mit konstanten Varianzen vorliegen würden. . Lineare Modelle in R Er zeigt deutlich, dass die Variabilität der Residuenmit wachsenden angepassten Miethöhen stark steigt. Damit aber bestimmen Kovariablenwerte, die zu großen angepassten Werten führen, die Regressionskoeffizienten systematisch stärker als solche mit geringen angepassten Miethöhen. (a) (b) Abbildung 9.3: Diagnostische Graphiken. a) Q-Q Plot der Quantile der standardisierten Residuen gegen die Quantile einer Normalverteilung, b) Standardisierte Residuen gegenHebelpunkte. Der Quantil-Quantil-Plot der standardisierten Residuen gegen die Quantile der Normalverteilung in Abbildung . a) zeigt für die Daten über Mieten im Mikrozensus eine deutliche Abweichung von einer Normalverteilung (wären die Residuen annähernd normal verteilt, sollten sie in der Nähe der eingezeichneten Geraden liegen). Nun ist die Verteilung der Residuen für die Gültigkeit von Regressionsergebnissen weitgehend irrelevant. Was man an einer solchen Graphik dennoch erkennt, ist die sehr asymmetrische Verteilung der Residuen. Das heißt aber, dass die angepasstenWerte die tatsächliche Miethöhe nur recht ungenau repräsentieren. Es gibt weit gröbere Unterschätzungen der Miethöhe als es Überschätzungen gibt. Zudem ist die Verteilung bei den Unterschätzungen deutlich anders als bei den Überschätzungen. Die zweite Graphik, . b), die die Residuen gegen die Hebelpunkte abträgt, zeigt zudem einige extrem große Hebelpunkte. Der Durchschnittswert der huu ist hier / . , das Maximum der Werte ist aber . , das -fache des Durchschnittswerts. Zudem kommen recht viele relativ große Hebelpunkte vor: Immerhin haben einen Wert über dem Dreifachen des Durchschnitts. Wir werden im nächsten Abschnitt drei Methoden vorstellen, mit denen diese Probleme behandelt werden können. Zunächst aber soll noch auf die Möglichkeiten hingewiesen werden, die R bereitstellt, um auf einzelne diagnostische Größen zugreifen zu können. Man erhält die (unstandardisierten) Residuen aus einem Modell durch residuals(), die Hebelpunkte durch hatvalues(). Die Standardabweichung der Residuen kann man dem summary() Objekt entneh- Lineare Regression men. Um also die standardisierten Residuen zu berechnen (und nicht nur zu malen), kann man schreiben: > huu < hatvalues(erg1) #Hebelpunkte > sig < summary(erg1)$sigma #Standardabweichung > residnorm < residuals(erg1)/sig/sqrt(1 huu) Der Befehl rstandard() liefert das gleiche Ergebnis. Ein weiteres wichtiges Instrument der Modellbeurteilung ist eine Statistik, die ausweist, wie sich Regressionskoeffizienten verändern, wenn einzelne Beobachtungen ausgeschlossen werden. Man nennt diese Größen im Anschluss an die Bezeichnungen im klassischen Buch von Belsley, Kuh und Welsch dfbeta. Der gleichnamige R Befehl dfbeta() ergibt eine Matrix der Dimension n m + , die die Änderungen derm + Koeffizienten (die Spalten der Matrix) für jeden ausgeschlossenen Fall (die Zeilen der Matrix) enthält. Die Werte der dfbeta ermöglichen auch eine andere Interpretation der Standardfehler von Koeffizienten, die traditionell nur im Rahmen des klassischen Modells der Statistik abgeleitet werden. Denn man kann die Variabilität der Regressionskoeffizienten ja auch ganz unabhängig vom Standardmodell als Variabilität der Koeffizienten betrachten, wenn einzelne Beobachtungen ausgeschlossen werden. Diese Sichtweise unterstellt gar kein stochastisches Modell, nicht ein mal die Einbettung in einen beliebigen, möglicherweise komplizierten stochastischen Modellzusammenhang. Probieren wir die Idee, dann ergibt sich: > v0 < vcov(erg1) > dfb < dfbeta(erg1) > v1 < crossprod(dfb) > s0 < sqrt(diag(v0)) > s1 < sqrt(diag(v1));s1 (Intercept) qm HHGroesse ErwPerson 6.636 0.131 2.221 2.279 > (s1 s0)/s0 (Intercept) qm HHGroesse ErwPerson 0.3899 0.7235 0.2208 0.0255 Hierbei haben wir verwandt, dass die Mittelwerte der dfbeta immer sind. Daher entspricht crossprod(dfb) gerade var(dfb)∗(n 1), der Varianz-Kovarianz- David A. Belsley, Edwin Kuh, Roy E. Welsch: Regression Diagnostics. , New York: Wiley. Es gibt auch eine normierte Version dfbetas(). Wir finden sie weniger nützlich, weil die Größen dfbeta sich direkt auf die ja auch inhaltlich interpretierbaren Koeffizientes des Modells beziehen. Weiterhin gibt es eine Variante standardisierter Residuen rstudent(), die die Varianz der Residuen ohne den jeweiligen Fall berechnen. Große Ausreißer führen ja auch zu großen geschätzten Varianzen, die die Beurteilung von Modellproblemen erschweren. Die Funktion rstudent() ist daher besonders in kleinen Datensätzen gut geeignet, Modellabweichungen aufzudecken. . Lineare Modelle in R matrix der dfbeta. Die aus den dfbeta abgeleiteten Standardabweichungen weichen stark von den Standardabweichungen ab, die durch summary(erg1) berechnet werden. Die geschätzte Standardabweichung für den Koeffizienten der Wohnungsgröße ist gar um % größer als in der summary(erg1) Zusammenfassung. Auch das deutet auf einige Probleme mit dieser Variablen hin. . . Ausschluss von Ausreißern Wir suchen zunächst noch einmal Werte mit sehr großen Hebelpunkten. Die Beobachtungen mit den zehn größten Hebelpunkten haben die Werte > o < order(huu,decreasing=T) > Wohn[o[1:10],1:4] Miete qm HHGroesse ErwPerson 3908 1150 260 1 1 997 650 250 1 1 1200 820 240 2 1 805 1000 240 2 2 3000 605 117 9 1 2646 400 200 7 0 533 408 60 7 0 3497 450 200 2 1 2407 600 200 6 1 1267 500 124 8 2 Die großen Hebelpunkte ergeben sich also insbesondere bei sehr großen Wohnungen oder bei sehr großen Haushalten oder bei einer Kombination großer Haushaltsgrößen mit kleinenWohnungen bzw. wenigen Erwerbspersonen. Man könnte nun nach inhaltlichen Kriterien versuchen, einige dieser Werte auszuschließen. Dabei sind die Plausibilität der Angaben und das Interesse an entsprechenden Konstellationen (etwa große Haushalte in kleinenWohnungen) ausschlaggebend, kaum aber statistische Erwägungen. Wir gehen hier aber summarisch vor und demonstrieren nur den statistischen Effekt von Beobachtungen mit großen Hebelpunkten. Dazu schließen wir alle Beobachtungen aus, Die Varianz-Kovarianz-Matrix v1 ist eine Variante der so genannten heteroskedastiekonsistenten (Eicker-Huber-White) Varianzschätzungen. Sie werden u.a. vom Paket sandwich bereitgestellt. Der Befehl vcovHC(erg1) (mit der voreingestellten Variante HC3) ist numerisch identisch mit v1. Nähere Informationen über die Varianten des Pakets sandwich finden sich bei Achim Zeileis: Econometric computing with HC and HAC covariance matrix estimators. Journal of Statistical Software, ( ), , S. – (http://www.jstatsoft.org/v11/i10/) und Achim Zeileis: Object-oriented computation of sandwich estimators. Journal of Statistical Software, ( ), , S. – (http://www.jstatsoft.org/v16/i09/). Eine modellbasierte Rechtfertigung dieser Variante benutzt die Tatsache, dass es sich bei den dfbeta um eine empirische Version der Einflussfunktion handelt. Lineare Regression die mehr als . mal so groß wie der Mittelwert derHebelpunkte sind. Das betrifft Beobachtungen. Entsprechend formal verfahren wir bei großen Residuen. Zunächst betrachten wir große positive Residuen. Die größten Residuen sind > oo < order(residnorm,decreasing=T) > Wohn[oo[1:10],1:4] Miete qm HHGroesse ErwPerson 3055 1700 150 4 2 276 1500 125 1 1 2182 1700 180 4 1 2443 1400 105 2 2 3517 1400 115 4 1 335 1700 180 2 2 2713 1300 100 1 1 1173 1250 120 4 1 3406 1150 96 4 1 3466 1200 110 5 2 Diese Größenordnungen von Mieten für die zugehörigen Wohnungsgrößen sind durchaus möglich, wenn auch relativ groß. Ein anderes Bild ergibt sich bei den sehr kleinen negativen Residuen. Dort erhalten wir > ooo < order(residnorm) > Wohn[ooo[1:10],1:4] Miete qm HHGroesse ErwPerson 2489 100 140 3 2 614 180 150 1 1 5003 126 126 3 3 5209 255 170 1 0 2100 200 150 1 1 997 650 250 1 1 3269 40 100 2 2 3497 450 200 2 1 2646 400 200 7 0 396 125 113 1 1 Diese Angaben implizieren extrem geringe Mieten je qm. Sie mögen auf speziellen Arrangements zwischen Mietern und Vermietern beruhen und durchaus wichtig für die Beurteilung des Wohnungsmarktes sein. Wir gehen Der Ausschluss von Fällen allein aufgrund statistischer Überlegungen und Maßzahlen ist selten zu rechtfertigen. Maßzahlen wie Hebelpunkte, Residuen und dfbeta hängen stark von dem unterstellten Modell ab, die Modellwahl selbst aber ist in sehr vielen Fällen allein durch die Statistikkenntnisse und die Phantasie des Benutzers, die verwandten Statistikprogramme und die Traditionen des Faches bestimmt. Es ist uns wichtig zu betonen, dass ungenügende Werte von Diagnostiken immer auch auf die Möglichkeit verweisen, dass das Ausgangsmodell unangemessen ist oder wesentliche Aspekte der Situation bei der Modellwahl unberücksichtigt geblieben sind. . Lineare Modelle in R aber wieder nur formal vor und schließen Werte mit negativen normierten Residuen kleiner als - . aus. Außerdem schließen wir die Daten mit den größtenHebelpunkten aus. Als Grenzwert haben wir (willkürlich) das . -fache des Durchschnitts derHebelwerte gewählt. Die Datenauswahl ist dann nach dem Ausschluss fehlender Werte durch > Wohn2 < Wohn[huu<10/nrow(Wohn)&residnorm> 2.5,] > detach(Wohn) > attach(Wohn2) gegeben. Insgesamt sind Beobachtungen ausgeschlossen worden. Das neue Regressionsergebnis ist > erg3 < lm(Miete~qm+HHGroesse+ErwPerson) > summary(erg3) Coeicients: Estimate Std. Error t value Pr(>|t|) (Intercept) 34.007 5.077 6.70 2.3e 11 ∗∗∗ qm 4.285 0.087 49.27 < 2e 16 ∗∗∗ HHGroesse 6.542 2.062 3.17 0.0015 ∗∗ ErwPerson 24.859 2.253 11.03 < 2e 16 ∗∗∗ Residual standard error: 113 on 5387 degrees of freedom Multiple R squared: 0.42, Adjusted R squared: 0.419 F statistic: 1.3e+03 on 3 and 5387 DF, p value: <2e 16 . . Transformation der abhängigen Variablen Das Problem der sehr asymmetrischen Verteilung der Residuen und die ungleiche Variabilität der Residuen für verschiedene angepasste Werte kann man oft durch eine passende Transformation der abhängigen Variablen beheben. Aber was ist eine angemessene Transformation? Bei positiven Werten der abhängigen Variablen nimmt man oft den Logarithmus, weil er sicher zu kleineren positiven Residuen und einer symmetrischeren Verteilung der Residuen führt. Aber ist das auch für die Mieten gerechtfertigt? Box und Cox haben vorgeschlagen, die folgende Familie von Transformationen zu betrachten: T(y, λ) := yλ λ λ > log(y) λ = (9.3) George E.P. Box, David R. Cox: An analysis of transformations (with discussion), Journal of the Royal Statistical Society B, , , S. – . Die von Box und Cox vorgeschlagene Methode hat immer wieder zu Diskussionen geführt, denn Regressionen mit verschieden transformierten abhängigen Variablen sind nur schwer vergleichbar. Das gilt um so mehr, wenn die Transformation auf der Basis der Daten gewählt wird. Einen guten Einstieg in die Diskussion bietet der Artikel von David Hinkley und G. Runger:The analysis of transformed data (with discussion), Journal of the American Statistical Association, , , S. – . Lineare Regression λ Abbildung 9.4: Beste Werte für den Transformationsparameter λ in denDaten überMiethöhen. Diese Familie enthält neben dem Logarithmus auch alle positiven Potenzen von y, insbesondere also auch die Wurzeln, die auch oft als Standardtransformation benutzt werden. Man kann dann für jedes vorgelegte λ die zugehörigen optimalen Regressionsparameter bestimmen und dann denjenigen Wert λ̂ wählen, der die quadratischen Abstände (9.2) zwischen den transformierten abhängigen Variablen und der besten zu dieser Transformation passenden Regression minimiert. Abbildung . zeigt die (negativen und durch die Varianz der jeweiligen Residuen normierten) Werte von (9.2), der Anpassung des Modells an die Daten, wobei jeweils für ein vorgegebenes λ eine entsprechende Regression berechnet wird. Die senkrechten Linien geben ein approximatives % Konfidenzintervall für λ an, das also für diese Daten zwischen . und . liegt. Das Paket MASS stellt den Befehl boxcox() bereit, der diese Graphik erzeugt. Wir haben > boxcox(Miete ~ qm+HHGroesse+ErwPerson, + lambda=seq(0.05,0.27,length=99)) benutzt. Die Graphik zeigt, dass die Wahl des Logarithmus sicher zu radikal wäre. Der optimale Wert von λ ist etwa . / . , so dass wir als nächstes eine Regression der . -ten Wurzel der Miethöhe auf unsere drei Variablen betrachten: > traMiete < (Miete^(1/6.5) 1)∗6.5 > erg4 < lm(traMiete ~ qm+HHGroesse+ErwPerson) > summary(erg4) Coeicients: Estimate Std. Error t value Pr(>|t|) (Intercept) 7.012773 0.035098 199.80 <2e 16 ∗∗∗ qm 0.031022 0.000601 51.60 <2e 16 ∗∗∗ HHGroesse 0.036227 0.014254 2.54 0.011 ∗ ErwPerson 0.172660 0.015574 11.09 <2e 16 ∗∗∗ Residual standard error: 0.781 on 5387 degrees of freedom Multiple R squared: 0.444, Adjusted R squared: 0.444 F statistic: 1.44e+03 on 3 and 5387 DF, p value: <2e 16 Wir betrachten noch einmal den Q-Q-Plot für dieses Modell, um zu sehen, ob sich die Probleme asymmetrischer Residualverteilung und großer Ausreißer reduziert haben. Das Ergebnis ist Abbildung . a). Die positiven Residuen sind nun deutlich kleiner und liegen näher an der Geraden. Zudem ist die Verteilung . Lineare Modelle in R nun fast symmetrisch. Zudemwird eine fast konstante Varianzfunktion erreicht. Das zeigt die Graphik . b). (a) (b) Abbildung 9.5: Diagnostische Graphiken nach der Box-Cox Transformation. a) Q-Q Plot der Quantile der standardisierten Residuen gegen die Quantile einer Normalverteilung, b) Wurzel der Beträge der standardisierten Residuen gegen angepasste Werte. Aber wie interpretiert man das Ergebnis? Schließlich ist man an den Miethöhen interessiert, nicht an den . -ten Wurzeln der Mieten. Es ist naheliegend, aber falsch, die Regression mit der . -ten Wurzel einfach durch die . -te Potenz der Regression zurückzutransformieren. Denn wenn der Mittelwert M(Y / . X = x) tatsächlich ungefähr gleich b̂ + b̂ x + . . . + b̂mxm ist, dann ist M(Y X = x) sicher nicht gleich (M(Y / . X = x)) . . Das gilt natürlich für alle Transformationen in der Box-Cox Familie (9.3). Eine einfache Möglichkeit der korrekten Darstellung der Regressionsergebnisse nach der Retransformation in Miethöhen liefert der Durchschnitt der inversen Transformation T- (.) angewandt auf die Regressionsergebnisse an der vorgegebenen Stelle x = (x , . . . , xm), also M̂(Y X = x , . . . ,Xm = xm) = n u U T λ̂ (b̂ + b̂ x + . . . + b̂mxm + ê(u)) In R lässt sich diese Retransformation einfach realisieren. Will man etwa den Einfluss der Wohnungsgröße zwischen und qm für Haushalte mit einer Erwerbsperson und drei Personen im Haushalt zeigen, dann kann man schreiben: > NeuDat < data.frame(qm=20:100,HHGroesse=3,ErwPerson=1) > xb < predict(erg4,NeuDat) Naihua Duan hat das Verfahren näher untersucht und mit anderen Vorschlägen verglichen (Smearing estimate: A nonparametric retransformation method. Journal of the American Statistical Association, , , S. – ). Lineare Regression > e < residuals(erg4) > Werte < sapply(xb,function(xb) mean(((xb+e)/6.5+1)^6.5)) > plot(20:100,Werte,xlab="qm",ylab="Angepasste Miete") Abbildung 9.6: Retransformierte Miethöhen für einen Haushalt mit drei Personen. Dabei enthält der Dataframe NeuDat Zeilen mit den Variablen HHGroesse, qm und ErwPerson, also die gleichen Variablennamen wie die in der Regression benutzten. Der predict() Befehl nimmt als erstes Argument ein Regressionsmodell und als zweites einen Datensatz, der alle unabhängigen Variablen enthalten muss, die im Regressionsmodell benutzt werden. Berechnet werden dann die angepassten Werte b̂ + b̂ x + . . . + b̂mxm für die im zweiten Argument übergebenen Daten. In der dritten Zeile werden nochmals die Residuen berechnet, dann werden für alle Elemente in xb die retransformierten Werte berechnet. Das Ergebnis ist Abbildung . , in der dieMieten fürWohnungen zwischen und qm und keiner Erwerbsperson(gestrichelt), Erwerbsperson (durchgezogen) und Erwerbspersonen (gepunktet) dargestellt sind. . . Transformationen der Kovariablen Nun kann man nicht nur die abhängige Variable transformieren, um eine angemessenere Repräsentation der Daten zu finden. Man kann das auch mit den unabhängigen Variablen tun. Denn das Modell bleibt auch dann linear, wenn man an Stelle der ursprünglichen Variablen (bekannte) Funktionen der Variablen benutzt. Ein Modell der Form b + b g (X ) + . . . + bmgm(Xm) in dem g (.), . . . , gm(.) bekannte Funktionen der Kovariablen sind, nennt man oft ein additivesModell. In unserem Beispiel scheint es zwar plausibel, dass die Wohnungsgröße die Mieten linear bestimmen sollte. Wäre es die einzige Bestimmungsgröße, dann sollte man ja sogar erwarten, dass sich eine Regression durch den Nullpunkt (Konstante = ) ergeben sollte. Der relativ große Wert der Konstanten in den untransformierten Regressionen macht aber deutlich, dass sich die Wohnungspreise nicht so homogen verhalten. Zudem ist der Zusammenhang nach der Transformation zwangsweise nicht mehr linear. Die Transformation erzwingt gerade einen konkaven Verlauf, der aber durch die Daten gar nicht gerechtfertigt sein muss. In der Tat ist der Verlauf in Abbildung . eher unerwartet. Er impliziert eine besonders langsame Steigerung der Miethöhen für kleinere Wohnungen und eine überproportionale Steigerung . Lineare Modelle in R für große Wohnungen. Insbesondere von der Nachfrageseite her würde man aber eher das Gegenteil erwarten. Man könnte als erste Näherung für einen nicht-linearen Einfluss der Wohnungsgröße ein Polynom niedrigen Grades in der Variablen qm verwenden. Man tut das auch oft, muss dann aber mit mehreren Problemen rechnen: Zum einen sind Polynome schon durch wenige Punkte eindeutig bestimmt, so dass einige wenige Datenpunkte ihren gesamten Verlauf bestimmen können. Daher können sie auch nicht an lokale Besonderheiten der Daten an bestimmten Stellen x angepasst werden. Zudem muss man technische Probleme erwarten, weil zumindest über relativ kleine Intervalle die Werte etwa von X(u) und X(u) stark korrelieren werden. Eine Alternative sind Splinefunktionen, die sich stückweise aus Polynomen zusammensetzen. Dabei setzt man die Stücke so zusammen, dass sie an den Stellen, an denen sie zusammentreffen, die gleichen Werte haben und der Übergang so glatt wie möglich ist. Das vermindert den Einfluss einzelner Punkte auf den gesamten Kurvenverlauf. Und es erlaubt die Darstellung beliebiger stetiger Funktionsverläufe. Das Paket splines stellt die wichtigsten Splinefunktionen zur Verfügung. Wir verwenden einen so genannten natürlichen Spline, ein stückweises Polynom dritten Grades, das zusätzlich durch das Verhalten an den Rändern des Definitionsbereichs (im Regressionszusammenhang: demWertebereich von X(u), u ) beschränkt ist. Dort soll es eine lineare Funktion sein. Die Umsetzung in R ist mit dem Paket splines sehr einfach. Der Befehl ns() berechnet eine Splinebasis für einen natürlichen Spline. Die Option df=4 gibt an, dass drei Knoten (die Stellen, an denen die Polynomteile zusammengefügt werden sollen) gewählt werden. Die Stellen der Knoten werden, wenn im Befehl nichts anderes angegeben wird, an den entsprechenden Quantilen der Verteilung von X(u) gesetzt. Je größer man den Parameter df wählt, desto irregulärere Verläufe lassen sich darstellen. Auf der anderen Seite müssen dann auch mehr Parameter geschätzt werden, so dass insbesondere in kleineren Datensätzen df nicht größer als oder gewählt werden sollte. > library(splines) > erg5 < lm(traMiete ~ ns(qm,df=4)+HHGroesse+ErwPerson) > summary(erg5) Estimate Std. Error t value Pr(>|t|) (Intercept) 6.8110 0.1001 68.05 <2e 16 ∗∗∗ Zwar erzwingt die relativ große Konstante für einen großen Bereich an Wohnungsgrößen eine eher fallende Miethöhe je Quadratmeter für größere Wohnungen, der marginale Effekt aber ist konkav. Der Name „Spline“ verweist auf die Konstruktion von Schiffsplanken, bei der elastische Materialien wie Stahl oder Holz an mehreren Stützstellen (den sogenannten Knoten) in eine bestimmte Richtung gebogen werden. Die natürlichen Splines sind die Lösung des mathematischen Problems, bei vorgegebenen Knoten die Gesamtkrümmung der Planke zu minimieren. Wie bei Regressionenmit Potenztermen wie X, X etc. kannman Splines als Linearkombination von transformierten Werten der unabhängigen Variablen schreiben. Lineare Regression ns(qm,df = 4)1 2.2200 0.0950 23.38 <2e 16 ∗∗∗ ns(qm,df = 4)2 2.5946 0.0813 31.91 <2e 16 ∗∗∗ ns(qm,df = 4)3 5.1456 0.2215 23.23 <2e 16 ∗∗∗ ns(qm,df = 4)4 3.3136 0.1156 28.66 <2e 16 ∗∗∗ HHGroesse 0.0364 0.0143 2.55 0.011 ∗ ErwPerson 0.1764 0.0155 11.36 <2e 16 ∗∗∗ Residual standard error: 0.777 on 5384 degrees of freedom Multiple R squared: 0.45, Adjusted R squared: 0.449 F statistic: 733 on 6 and 5384 DF, p value: <2e 16 Das Ergebnis des summary() Befehls ist in diesem Fall aber wenig hilfreich, weil nur die Koeffizienten der Elemente der Splinebasis angegeben werden. Es ist wieder der predict() Befehl, der weiterhilft. > NeuDat < data.frame(qm=20:130,HHGroesse=3,ErwPerson=1) > Traqm < predict(erg5,NeuDat,type="terms",terms=1) > plot(20:130,Traqm,type="l",xlab="qm") Hier benutzen wir den predict() Befehl mit der Option type="terms", um uns die Terme des linearen Modells berechnen zu lassen. Da der Spline der erste Term des Modells ist (was terms(erg5) bestätigt), wählen wir ihn mit der Option terms=1 aus. Das Ergebnis ist in Abbildung . a) wiedergegeben. Auf der transformierten Skala erhalten wir also einen leicht konkaven Verlauf. (a) (b) Abbildung 9.7: Splinefunktion. a) Splinefunktion auf der transformierten Skala. b) Retransformierte Regressionsfunktion mit Splinefunktion (durchgezogene Linie) und mit linearem Term (gepunktete Linie). Der graue Bereich ist ein punktweises % Konfidenzintervall der Regression. Die Graphik . b) zeigt die retransformierte Funktion zusammen mit einem % Konfidenzintervall sowie zum Vergleich die entsprechende Funktion aus demModell mit dem linearen Einfluss der Wohnungsgröße. Um diese Graphik zu erstellen, benutzen wir die predict() Funktion mit der Option interval="conĄdence". Mit dieser Option berechnet predict() neben den an- . Lineare Modelle in R gepassten Werten auch die zugehörigen Konfidenzintervalle und gibt das Ergebnis als Matrix mit drei Spalten zurück. > xb < predict(erg5,NeuDat,interval="conĄdence") > e < residuals(erg5) > Werte < sapply(xb[,1],function(xb)mean(((xb+e)/6.5+1)^6.5)) > WerteU < sapply(xb[,2],function(xb)mean(((xb+e)/6.5+1)^6.5)) > WerteO < sapply(xb[,3],function(xb)mean(((xb+e)/6.5+1)^6.5)) Die Graphik erhält man dann mit > plot(20:130,Werte,type="l",bty="l",xlab="qm") > polygon(c(20:130,130:20),c(WerteU,rev(WerteO)), col=rgb(0.8,0.8,0.8,0.5),border=NA) Es zeigt sich, dass die retransformierten angepassten Werte fast linear in der Wohnungsgröße sind. Sie steigen für kleine Wohnungsgrößen stärker, für große Wohnungen langsamer als in der linearen Spezifikation. Das erscheint uns plausibler als das Ergebnis der linearen Spezifikation. Wir wollen am Schluss dieses Abschnitts noch auf die Struktur der Spline- Objekte eingehen, die etwa durch den Befehl ns() erzeugt werden. Man kann etwa für die Wohnungsgrößen eine Splinebasis durch > sp < ns(qm,df=4) > names(attributes(sp)) [1] "dim" "dimnames" "degree" "knots" [5] "Boundary.knots" "intercept" "class" erzeugen. Die Splinebasis besteht aus transformierten Werten der Variablen qm. Es ist eine Matrix der Dimension n df zusammen mit einigen Attributen, die die Eigenschaften des benutzten Splines angeben. Das ist ganz ähnlich zu der Kodierung eines Polynoms dritten Grades, bei der man die transformierten Variablen X(u), X(u) , und X(u) jeweils als Regressoren bildet, um dann eine einfache lineare Regression mit diesen neuen Kovariablen zu berechnen. Der Unterschied besteht darin, dass die Basis (die transformierten Werte der Variablen qm in unserem Beispiel) so gewählt werden muss, dass sowohl die Randbedingungen (Linearität jenseits der Randpunkte) als auch die Bedingungen über den Zusammenhang der verschiedenen Polynomstücke erhalten bleiben. Der ns() Befehl integriert diese Anforderungen und minimiert die Korrelation zwischen den verschiedenen Komponenten. Die Komponenten des Ergebnisses des ns() Befehls können deshalb einfach als Regressoren in den lm() Befehl einbezogen werden. Allerdings hängt das Ergebnis des ns() Befehls immer von der Verteilung der benutzten Daten ab und kann nicht einfach auf neue Daten übertragen werden. Der Befehl attr() erlaubt es, Informationen über die benutzte Splinebasis zurückzuerhalten. Insbesondere liefert > attr(sp,"Boundary.knots") [1] 10 140 > attr(sp,"knots") 25% 50% 75% 52 64 79 Lineare Regression > attr(sp,"degree") [1] 3 die Position der Grenzknoten, hinter denen der Spline linear sein soll, die Position der inneren Knoten (hier: die Quartile bei df=4) und den Grad der Polynome. Zumeist muss man Splinebasen dann direkt manipulieren, wenn man die Splinebasen konstant halten möchte, um verschiedene Modelle mit unterschiedlichen Variablen und unterschiedlichen Fallzahlen noch sinnvoll vergleichen zu können. Dafür stellt das Paket splines einen eigenen predict() Befehl bereit. Will man etwa die Werte der Splinebasis für die Wohnungsgrößen von bis qm benutzen, aber die gleichen Parameter (Lage der Knoten, Grad der Polynome) benutzen, die sich bei der Berechnung mit allen Daten ergeben haben, dann kann man zunächst explizit die Splinebasis für alle Beobachtungen durch > Splineqm < ns(qm,df=4) erzeugen. Die Klasse dieses neuen Objektes ist > class(Splineqm) [1] "ns" "basis" "matrix" Nun definiert man sich neue Daten und ruft den predict() Befehl auf > Neuqm < 20:22 > neuSplineqm < predict(Splineqm,Neuqm) > neuSplineqm 1 2 3 4 [1,] 0.00639 0.100 0.225 0.125 [2,] 0.00851 0.109 0.246 0.136 [3,] 0.01104 0.118 0.266 0.147 Dies ist wieder ein Beispiel für die Verwendung objektorientierter Methoden in R. predict() ist eine generische Funktion, die auf Grund des Klassenattributes des übergebenen Objektes prüft, ob es dafür eine spezielle Methode gibt. Wenn das der Fall ist, wird die entsprechende Funktion aufgerufen. In diesem Fall gehört das an predict() übergebene Objekt zur Klasse ns. Deshalb wird predict.ns() aufgerufen. Diese Funktion ist in der Hilfe zum Paket splines dokumentiert. . Regression in Matrixnotation R stellt eine Reihe von Matrixoperationen zur Verfügung. Sie bilden die numerische Grundlage für viele statistische Berechnungen: lineare und nicht- Will man sich die Definition der Funktion ansehen, reicht es nicht, einfach ihren Namen (ohne die ()) aufzurufen. Denn die Funktion ist im Namespace des Pakets enthalten. Zugang erhält man durch die Angabe splines:::predict.ns. Alle Funktionen zu einer Klasse erhält man durch methods(class=ns). Man spricht von Methoden in Anlehnung an den allgemeinen Sprachgebrauch in der objektorientierten Programmierung. . Regression in Matrixnotation lineare Regressionen, Faktorenanalysen, Diskriminanzanalysen etc. Die Funktion lm() zur Schätzung linearer Modelle wurde bereits behandelt. In diesem Kapitel soll nun das lineare Regressionsmodell in Matrixnotation behandelt werden. Dabei werden die R-Befehle für Matrizen eingeführt und diskutiert. . . Einfache Matrixoperationen in R Zunächst wollen wir uns wichtige Matrixoperationen anhand einiger Beispiele veranschaulichen. Mit dem Befehl matrix() können Matrizen erzeugt werden. So erzeugte Matrizen haben immer ein Attribut „Dimension“: > x < matrix(1:3) > dim(x) [1] 3 1 Wie ersichtlich, wird ohne Festlegung der Dimension ein Spaltenvektor erzeugt. Ganz analog kann eine Matrix mit vorgegebenen Dimensionen erzeugt werden: > x < matrix(1:6,nrow=2,ncol=3);x [,1] [,2] [,3] [1,] 1 3 5 [2,] 2 4 6 DasMultiplikationszeichen „*“ bewirkt die skalare, elementweiseMultiplikation der Elemente der Matrix: > x∗2 [,1] [,2] [,3] [1,] 2 6 10 [2,] 4 8 12 Um eine Matrixmultiplikation durchzuführen, ist anstelle von ∗ der Befehl %∗% zu verwenden: > y < matrix(c(3,12.4, 0.1),nrow=3) > x%∗%y [1,] 39.7 [2,] 55.0 Was passiert, wenn Rs Vektoren mit Matrizen multipliziert werden sollen? Rs Datentyp vector hat keine Dimension (genauer: kein Attribut dim) und daher ist nicht von vornherein klar, wie (oder ob überhaupt) Rs Vektoren mit Rs Matrizen multipliziert werden können. R verwendet beim binären Operator%∗% zwischenMatrizen und Vektoren (oder Vektoren und Vektoren) Konventionen, die durch die folgenden Beispiele wohl am einfachsten illustriert werden: > u < 1:3 > v < 4:6 > w < 7:8 > x%∗%u Lineare Regression [,1] [1,] 22 [2,] 28 > w%∗%x [,1] [,2] [,3] [1,] 23 53 83 > u%∗%x Error in u %∗% x : non conformable arguments > u%∗%v [,1] [1,] 32 Bei der Matrixmultiplikation von Matrizen mit R Vektoren werden die Vektoren also als entsprechende Matrizen interpretiert, jedenfalls dann, wenn das überhaupt möglich ist. Das Matrixprodukt %∗% zweier Vektoren ist immer das innere Produkt, also i uivi, das dann alsMatrix der Dimension x gespeichert wird. Eine Matrix kann mittels der Funktion t() transponiert werden: > t(matrix(1:6,nrow=3,ncol=2)) [,1] [,2] [,3] [1,] 1 2 3 [2,] 4 5 6 Matrixmultiplikationen können anstelle von „%*%“ auchmit den Funktionen crossprod() bzw. tcrossprod() durchgeführt werden, die insbesondere bei großen Matrizen wesentlich effizienter sind (crossprod() berechnet t(x)%∗%y, tcrossprod x%∗%t(y)): > matrix(1:2,ncol=1)%∗%matrix(3:4,ncol=2) [,1] [,2] [1,] 3 4 [2,] 6 8 > crossprod(matrix(1:2,ncol=2),matrix(3:4,ncol=2)) [,1] [,2] [1,] 3 4 [2,] 6 8 Die Berechnung der Inversen einer Matrix erfolgt durch den Befehl solve(): > x < matrix(1:4,nrow=2) > solve(x) [,1] [,2] [1,] 2 1.5 [2,] 1 0.5 Der Befehl löst somit das Gleichungssystem xa=I nach a auf, wobei I die Einheitsmatrix mit den Werten auf der Diagonalen und den Werten außerhalb der Diagonalen ist. solve(x) ist die Inverse der Matrix x. Dies können wir auch an einem Beispiel überprüfen: . Regression in Matrixnotation > solve(x)%∗%x [,1] [,2] [1,] 1 0 [2,] 0 1 Die Variante solve(x,b) löst das Gleichungssystem xa=b für einen Vektor b. Eine Einheitsmatrix kann mit dem Befehl diag() erzeugt werden: > diag(2) [,1] [,2] [1,] 1 0 [2,] 0 1 Der Befehl diag() ist allerdings kontextabhängig. Wird als Argument eine Matrix übergeben, wird ein Vektor der Diagonalelemente dieser Matrix erzeugt: > diag(matrix(1:4,nrow=2)) [1] 1 4 Wird ein Vektor übergeben, wird eine Diagonalmatrix mit den Vektorelementen auf der Diagonalen erzeugt: > diag(1:2) [,1] [,2] [1,] 1 0 [2,] 0 2 Nachfolgend sollen noch weitere nützliche Rechenregeln betrachtet und mit einfachen Beispielen illustriert werden. Wir können uns z.B. davon überzeugen, dass bei der Matrixmultiplikation das Kommutativgesetz (ab = ba) nicht gilt: > matrix(1:4,nrow=2)%∗%matrix(5:8,nrow=2) [,1] [,2] [1,] 23 31 [2,] 34 46 > matrix(5:8,nrow=2)%∗%matrix(1:4,nrow=2) [,1] [,2] [1,] 19 43 [2,] 22 50 Das Distributivgesetz (a + b) c = ac + bc hingegen gilt: > (matrix(1:4,nrow=2)+matrix(5:8,nrow=2))%∗%matrix(9:12,nrow=2) [,1] [,2] [1,] 154 186 [2,] 192 232 > matrix(1:4,nrow=2)%∗%matrix(9:12,nrow=2)+ + matrix(5:8,nrow=2)%∗%matrix(9:12,nrow=2) [,1] [,2] [1,] 154 186 [2,] 192 232 Ebenso gilt die Transponierungsregel (ab)t = btat: Lineare Regression > t(matrix(1:4,nrow=2)%∗%matrix(5:8,nrow=2)) [,1][,2] [1,] 23 34 [2,] 31 46 > t(matrix(5:8,nrow=2))%∗%t(matrix(1:4,nrow=2)) [,1] [,2] [1,] 23 34 [2,] 31 46 . . Das Regressionsmodell in Matrixnotation Ausgangspunkt ist das lineare Regressionsmodell Y(u) = b + b X (u) + b X (u) + . . . + bmXm(u) + e(u) Werden die n Beobachtungen untereinander geschrieben, dann lässt sich das Modell in Matrixnotation darstellen als: Y = Xb + e Dabei ist Y ein Spaltenvektor der Dimension n , X ist eine Matrix der Dimension n m + , b ist der Spaltenvektor der Regressionskoeffizienten mit der Dimension m + und e ist ein Spaltenvektor der Dimension n . Gesucht ist der Vektor b̂, der am besten zu den Daten passt. Wählt man die quadrierten Abstände zwischen Daten und der Funktion b + b X (u) + b X (u)+ ...+bmXm(u) := Xb als Kriterium, dann soll die Summe der Quadrate der Residuen ê(u) minimiert werden. Die Summe der quadrierten Residuen in Matrixnotation ist ete: ete = (Y Xb)t (Y Xb) = YtY btXtY + btXtXb Diese soll nun durch die Wahl von bminimiert werden. Es wird daher die erste Ableitung nach b gebildet und Null gesetzt: ∂ete ∂b b̂ = XtY + XtXb̂ = NachUmformung resultiert dann die Bestimmungsgleichung für den gesuchten Parametervektor b̂: XtXb̂ = XtY b̂ = XtX XtY Der Schätzwert der Varianz der Störterme ergibt sich in der Matrixnotation als σ̂ = êt ê/(n m) und die geschätzte Varianz-Kovarianz-Matrix des geschätzten Parametervektors b̂ als cov(b̂) = σ̂ XtX .XtX wird alsKreuzproduktmatrix bezeichnet. . Regression in Matrixnotation . . Die Berechnung von Regressionen in R Es soll nun eine lineare Regression in R berechnet werden. Als Beispiel erzeugen wir Daten mit folgendem Prozess: Y = + X + X + e e N ( , ) Die Werte der erklärenden Variablen erzeugen wir als Realisationen einer multivariat normalverteilten Zufallsvariablen. > set.seed(123) > library(MASS) > x12 < mvrnorm(n=5,mu=c(0,0),Sigma=cbind(c(1,0.8),c(0.8,1))) Als erste Spalte der Datenmatrix x setzen wir eine für die Konstante, die beiden erklärenden Variablen werden Spalten und (x12): > x < cbind(1,x12) Den bekannten Parametervektor b benötigen wir ebenfalls zur Datenerzeugung: > b < matrix(1:3,ncol=1) Der deterministische Teil wird additiv von einem Vektor der Störterme überlagert, den wir ebenfalls als Realisation der Normalverteilung erzeugen: > e < rnorm(5) Nun können wir den Vektor y erzeugen: > y < x%∗%b+e Betrachten wir die einzelnen Modellkomponenten. Der zu erklärende Datenvektor enthält nun die folgenden Werte: > t(y) [,1] [,2] [,3] [,4] [,5] [1,] 0.977 0.122 9.200 1.660 1.200 Die Matrix der erklärenden Variablen hat folgende Form: > t(x) [,1] [,2] [,3] [,4] [,5] [1,] 1.0000 1.0000 1.00 1.000 1.0000 [2,] 0.0106 0.0726 1.08 0.150 0.0183 [3,] 1.0741 0.3641 1.88 0.284 0.2636 Nun soll auf Basis der vorliegenden n = Werte der Vektor b geschätzt werden: > bd < solve(t(x)%∗%x)%∗%t(x)%∗%y; t(bd) [,1] [,2] [,3] [1,] 1.198 3.655 2.133 Außer für den Parametervektor könnten wir uns auch noch für weitere Modellinformationen interessieren. Summe der quadrierten Residuen (rss für residual sum of squares): Einige Hinweise auf effizientere numerische Möglichkeiten in R gibt Douglas Bates: Least squares calculations in R. Timing different approaches. R News, , , – . Lineare Regression > res < y x%∗%bd > rss < as.vector(t(res)%∗%res);rss [1] 0.4206 Geschätzte Varianz bzw. Standardabweichung der Residuen: > n < length(y) > m < ncol(x) > sigma2 < rss/(n m);sigma2 [1] 0.2103 > sigma < sigma2^0.5;sigma [1] 0.4586 Kovarianzmatrix der Parameter: > covb < sigma2∗solve(t(x)%∗%x);covb [,1] [,2] [,3] [1,] 0.04911 0.05666 0.01298 [2,] 0.05666 0.60301 0.23085 [3,] 0.01298 0.23085 0.13246 Vektor der Standardabweichungen der Parameter: > sb < sqrt(diag(covb));sb [1] 0.2216 0.7765 0.3640 Vektor der t-Werte: > tb < as.vector(bd/sb);tb [1] 5.408 4.707 5.861 Zweiseitige empirische Signifikanzniveaus (p-values): > pvb < 2∗pt( abs(tb),n m);pvb [1] 0.03253 0.04229 0.02790 R und korrigiertes R : > r2 < as.vector(var(x%∗%bd)/var(y));r2 [1] 0.9935 > r2.adj < 1 (n 1)/(n m)∗(1 r2);r2.adj [1] 0.9870 F-Statistik und Signifikanzniveau der F-Statistik: > f < as.vector(sum((x%∗%bd mean(y))∗∗2)/(m 1)/sigma2);f [1] 152.6 > pvf < 1 pf(f,m 1,n m); pvf [1] 0.00651 Vergleichen wir das Ergebnis mit dem Ergebnis der Prozedur lm(): > summary(lm(y~x12)) Call: lm(formula=y~x12) Residuals: 1 2 3 4 5 0.0769 0.0341 0.0453 0.4073 0.4955 . Übungsaufgaben Coeicients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.198 0.222 5.41 0.033 ∗ x121 3.655 0.777 4.71 0.042 ∗ x122 2.133 0.364 5.86 0.028 ∗ Residual standard error: 0.459 on 2 degrees of freedom Multiple R squared: 0.993, Adjusted R squared: 0.987 F statistic: 153 on 2 and 2 DF, p value: 0.00651 . Übungsaufgaben 1) Erzeugen Sie Realisationen einer standardnormalverteilten Zufallsvariablen (set.seed(1)) und weisen Sie die Realisationen der Variablen x zu. Berechnen Sie Regressionswerte xb durch 1+0.5∗x. Addieren Sie Realisationen (set.seed(3)) einer standardnormalverteilten Zufallsvariablen e und weisen Sie die so entstandenenWerte dem Vektor y zu. Stellen Sie x und y in einem Streudiagramm dar und zeichnen Sie die Funktion f (x) = + . x ein. Benutzen Sie dazu den Befehl abline(). Ermitteln Sie mittels der Funktion lm() die Regressionswerte und tragen Sie auch diese Gerade in Ihr Diagramm ein. (Benutzen Sie wieder den Befehl abline().) 2) Erzeugen Sie einen neuen Vektor v mit Realisationen einer standardnormalverteilten Zufallsvariablen (set.seed( )) und addieren Sie diese Werte zu den xb. Nennen Sie das Ergebnis y2. Tragen Sie auch die neuen Realisationen x,y in Ihr Diagramm ein. Benutzen Sie dazu eine andere Farbe oder andere Symbole. 3) Ermitteln Sie für die neuen Realisationen eine Regressionsgerade und tragen Sie diese zusätzlich in Ihr Diagramm ein. 4) Extrahieren Sie aus dem Objekt der ersten Regression den Vektor der Residuen. 5) Erzeugen Sie die voreingestellten diagnostischen Graphiken für das erste Regressionsmodell. 6) Extrahieren Sie den Wert des R aus dem summary() Objekt des ersten Regressionsmodells. 7) Simulieren Sie Realisationen von jeweils Werten einer standardnormalverteilten Variablen und addieren Sie sie zu den in Aufgabe ) erzeugten Werten von xb. Berechnen Sie für jede der Realisationen eine Regression und speichern Sie die Werte der berechneten Regressionskoeffizienten in einer Matrix (Sie extrahieren die Koeffizienten mit der Funktion coef()). Berechnen Sie die Standardabweichungen der beiden Koeffizienten über alle Realisationen und vergleichen Sie den Lineare Regression Wert mit den im ersten Regressionsmodell geschätzten Standardfehlern der beiden Koeffizienten. 8) Erweitern Sie die Berechnungen in der letzten Aufgabe und speichern Sie auch die geschätzten Standardfehler und Varianzen der Regressionsparameter. Sie können dazu z.B. den Befehl diag(vcov(erg)) verwenden. Berechnen Sie dann den Durchschnitt dieser Werte und vergleichen Sie ihn mit der Standardabweichung und der Varianz, die sich aus der Simulation der letzten Aufgabe ergab. 9) Erzeugen Sie je Realisationen standardnormalverteilter Variabler x1,x2,e und berechnen Sie y := x + exp(sin( x ) + e) Können Sie durch passende Transformationen die funktionale Form der Regressoren wiederentdecken? Was passiert, wenn Sie die Berechnungen mit Realisationen wiederholen? Können Sie sich Funktionen konstruieren, bei denen Ihre Methoden versagen würden? 10) Benutzen Sie die Angaben des Mikrozensus zu Miethöhen, ergänzen Sie aber die Regressionsgleichung um Angaben zum Haushaltseinkommen (ef539). Die Angaben sind gruppiert. Wie würden Sie mit diesen Angaben umgehen, wenn sie als unabhängige Variable in die Regression eingehen sollen? 11) Mit Residuen kannman erstaunliche Dingemachen. Besuchen Sie die Seite http://www4.stat.ncsu.edu/~stefanski/NSF_Supported/ Hidden_Images/stat_res_plots.html für einige Beispiele und entsprechenden R Code (am Schluss der Seite). Lesen Sie dazu die Erläuterungen von Leonard A. Stefanski: Residual (Sur-)Realism. American Statistician, , , S. – (der Artikel ist auch auf der angegebenen Seite erhältlich) und erzeugen Sie sich eine eigene Residuengraphik. 12) Erzeugen Sie eine Matrix X der Dimension mit Realisationen einer standardnormalverteilten Zufallsvariablen und berechnen Sie mittels der Matrixmultiplikation und der Funktion crossprod() XtX sowie mittels der Funktion solve() die Inverse (XtX) . a) Erzeugen Sie eine Variable Y als Summe einer normalverteilten Zufallsvariablen (mean= ,sd= ) und Xb. X sei dabei die um einen Eins-Vektor ergänzte Matrix X und b sei : . b) Ermitteln Sie den Vektor der Regressionsparameter nach der Methode der kleinsten Quadrate mit Hilfe von Matrizenoperationen. c) Ermitteln Sie für die Parameter die Standardfehler, die t-Werte und die p-Werte. d) Ermitteln Sie das Bestimmtheitsmaß.

Chapter Preview

References

Zusammenfassung

Vorteile

- Einführung in die statistische Analyse mit R für Wirtschafts- und Sozialwissenschaftler

- Inklusive hilfreicher Tipps wie "Ansprechende Grafiken mit R gestalten"

Zum Thema

R ist ein Statistikprogramm, das kostenlos über das Internet verbreitet wird und dessen Source Codes frei zugänglich sind.

Aufgrund dieses kostenlosen Angebots gehen immer mehr Dozenten dazu über, neben SPSS auch R zu lehren bzw. SPSS durch R zu ersetzen.

In R steht dem Nutzer die gesamte Bandbreite statistischer Verfahren zur Verfügung. Durch die eigenständige Programmierumgebung ist die Software sehr flexibel und erlaubt notwendige Modifikationen und Erweiterungen verfügbarer Prozeduren.

Zum Werk

Dieses Buch führt leicht verständlich in die statistische Analyse mit R ein. Anhand von Beispielen wird die Umsetzung der wichtigsten Methoden der Statistik, wie sie üblicherweise in den Grundkursen gelehrt werden, mit R vorgestellt.

Das Buch verfolgt entsprechend zwei Ziele:

1. Vorstellung der statistischen Methoden,

2. Benutzung des Werkzeuges R zur Analyse von Daten.

Inhalt

- Grundlagen von R

- Datenbehandlung und graphische Darstellungen mit R

- Datenbeschreibungen (deskriptive Statistik)

- Wahrscheinlichkeitsverteilungen

- Regressionsanalysen

- Optimierungsverfahren

- Simulationen mit R

Neben vielen neuen, wirtschaftsorientierten Beispielen wird nun auch in die Paneldatenanalyse und Stichprobentheorie eingeführt.

Zu den Autoren

Dr. Andreas Behr ist wissenschaftlicher Mitarbeiter am Institut für Statistik und Ökonometrie der Universität Münster.

Dr. Ulrich Pötter ist wissenschaftlicher Mitarbeiter am Institut für Statistik der Universität Bochum.

Zielgruppe

Für Studierende und Dozenten der Wirtschaftswissenschaften im Bachelor an Universitäten und Fachhochschulen.