############################################################################### ## Inhalt: ## Kapitel 7.1 (Regressionsmodelle) des Kurses von Ruckdeschel & Kohl ############################################################################### ########################################################### ## Lesen Sie bitte LineareModelle.pdf und Abschnitt 7.1.1 des Skripts von ## Ruckdeschel und Kohl! ########################################################### ########################################################### ## 1. Einfache lineare Regression ## Datensatz aus dem Paket ISwR (Introductory Statistics with R) von P. Dalgaard. ########################################################### library(ISwR) ## Wir laden die Daten und fügen Sie zum Suchpfad hinzu data(thuesen) ?thuesen attach(thuesen) plot(short.velocity ~ blood.glucose, ylab = "fasting blood glucose [mmol/l]", xlab = "mean circumferential shortening velocity [%/s]", main = "Ventricular shortening velocity - data", pch = 20) ####################################### ## Einfache Lineare Regression y = a + b*x + eps ## Der KQ-Schätzer für Lineare Modelle ist in der Funktion "lm" implementiert. fit1 <- lm(short.velocity ~ blood.glucose) ## Das Ergebnis ist von der S3-Klasse "lm" class(fit1) ## und enthält viele Informationen str(fit1) ## Es gibt eine Vielzahl von Methoden für Objekte der Klasse "lm" ## ein erstes, wichtiges Beispiel ist die Funktion "summary" summary(fit1) ## Mit dem F-Test wird das Modell gegen das Nullmodell (nur Achsenabschnitt) ## getestet. Im vorliegenden Fall der einfachen linearen Regression, ist ## der t-Test für den Parameter zur Variable "blood.glucose" äquivalent ## zu diesem F-Test. Daher auch die identischen p-Werte. ## Im vorliegenden Beispiel erreichen ist "blood.glucose" zum 5% Nivau ## signifikant und analog auch das Modell im Vergleich zum Nullmodell. ## R^2 ist der Anteil der Varianz, der durch das Modell erklärt ist. ## adjusted R^2: berücksichtigt zusätzlich die Anzahl der Parameter ## im Modell. ## Man kann R^2 nämlich einfach dadurch steigern, indem man die Anzahl der ## Parameter im Modell erhöht. ## vgl. auch http://en.wikipedia.org/wiki/Coefficient_of_determination ## und insbesondere http://www.forecastingprinciples.com/rulesforcheaters.html ## ("Wie erreiche ich ein hohes R^2?"). ## Im vorliegenden Beispiel ist R^2 mit 0.1737 bzw. das adjustierte R^2 ## mit 0.1343 relativ klein. ####################################### ## Einfache Lineare Regression: Zusammenhang zwischen R^2 und Korrelationskoeffizient cor.test(short.velocity, blood.glucose, use = "pairwise.complete.obs") cor(short.velocity, blood.glucose, use = "pairwise.complete.obs") sqrt(summary(fit1)$r.squared) ## In diesem Fall ist R^2 gerade das Quadrat des Korrelationskoeffizienten! ####################################### ## Konfidenzintervalle pred.frame <- data.frame(blood.glucose = 4:20) ## Auch die Funktion "predict" besitzt eine Methode für die S3-Klasse "lm" ## Mit "predict" können Vorhersagen berechnet werden. predict(fit1, newdata = pred.frame) ## "predict" kann aber auch dazu verwendet werden, um Konfidenzintervalle zu ## berechnen. ## Konfidenzintervall für die geschätzte Regressionsgerade fit1.c <- predict(fit1, pred.frame, interval = "confidence") ## Konfidenzintervall für die Vorhersage fit1.p <- predict(fit1, pred.frame, interval = "prediction") ## Plot von Daten, Regressionsgerade und Konfidenzintervallen plot(short.velocity ~ blood.glucose, ylab = "fasting blood glucose [mmol/l]", xlab = "mean circumferential shortening velocity [%/s]", main = "Ventricular shortening velocity - regression", ylim = range(fit1.p), pch = 19) matlines(pred.frame$blood.glucose, fit1.c, col = c(1, "darkred", "darkred"), lwd = 2, lty = rep(1, 3)) matlines(pred.frame$blood.glucose, fit1.p, col = c(1, "orange", "orange"), lwd = 2, lty = rep(1, 3)) legend("topleft", legend = c("KQ-Schaetzer", "confidence", "prediction"), fill = c(1, "darkred", "orange")) ####################################### ## Aufräumen: Löschen des Datensatzes vom Suchpfad detach(thuesen) ########################################################### ## 2. Beispiel: Regressionseffekt (regression to mediocrity) ## Beim wiederholten Testen eines Zusammenhangs beobachtet man oft das folgende ## Phänomen ## - Die Spitzengruppe wird im Durchschnitt (relativ) schlechter. ## - Die Schlusslichtgruppe wird im Durchschnitt (relativ) besser. ## Dieser sog. Regressionseffekt (regression to mediocrity) liegt an der ## Zwetschgenform der Datenwolke und kann bei der Regression des zweiten (Y) ## auf das erste Resultat (X) beobachtet werden. ## Man sollte sich durch dieses Phänomen, welches ein rein statistischen Effekt ## zur Ursache hat, nicht verleiten lassen und auf eine "wirkliche" Ursache schließen! ####################################### ## Beschrieben wurde dieses Phänomen zum ersten Mal von Francis Galton im Jahr 1877. ## Er untersuchte die Größe von Vätern und Söhnen. Der Datensatz findet sich im ## R Paket UsingR von J. Verzani library(UsingR) ## vgl. Beispiel in ?father.son data(father.son) attach(father.son) ## ähnlich zum Cover von Statistics von Freedman, Pisani, Purves und Adhikari plot(sheight ~ fheight, pch = 20, ylab = "Größe des Sohns", xlab = "Größe des Vaters") abline(a = 0, b = 1, lty = 2, lwd = 2) abline(lm(sheight ~ fheight), lty=1, lwd=2) ## kleine Väter haben also (relativ gesehen) größere Söhne, während große Väter ## (relativ gesehen) kleinere Söhne haben. ## Erklärung: ## Wir betrachten eine Normalverteilte Zufallsgröße mit Mittelwert 100 und ## Standardabweichung 10. ## Wir führen eine erste Messung durch und erhalten den Wert 120. ## Dann ist es bei der zweiten Messung viel wahrscheinlicher einen Wert < 120 ## zu erhalten als einen Wert > 120 curve(dnorm(x, mean = 100, sd = 10), from = 65, to = 135, n = 501) abline(v = 120, lty = 2) pnorm(120, mean = 100, sd = 10) pnorm(120, mean = 100, sd = 10, lower.tail = FALSE) ## Aufräumen detach(father.son) ########################################################### ## 3. Simpson-Paradoxon ## Die Kenntnis bzw. Unkenntnis eines kann zu völlig gegensätzlichen Ergebnissen ## führen. ## Studiendauer x <- c(12, 11, 11, 10, 12, 10, 11, 9) ## Jahreseinkommen nach Abschluss des Studiums y <- c(38, 38, 34, 40, 33, 30, 30, 35) plot(x, y, pch = 20, xlab = "Studiendauer in Semester", ylab = "Jahreseinkommen in TSD-Euro") abline(lm(y ~ x), lwd = 2) ## Also lieber länger studieren!? ## Nein!!! plot(x, y, pch = 20, xlab = "Studiendauer in Semester", ylab = "Jahreseinkommen in TSD-Euro", main = "Simpson-Paradoxon") abline(lm(y ~ x), lwd = 2) abline(lm(y[1:4] ~ x[1:4]), col = "blue", lwd = 2) abline(lm(y[5:8] ~ x[5:8]), col = "red", lwd = 2) legend("topright", legend=c("gesamt", "Note <=2", "Note > 2"), fill=c("black", "blue", "red")) ## Aufräumen rm(x, y) ########################################################### ## 4. Diagnostik ########################################################### ########################################################### ## Wir definieren nun die für die Regressionsdiagnostik sehr wichtigen Residuen. ## Gegeben seien Beobachtungen x_1, ..., x_n und y_1, ..., y_n sowie eine ## Regressionsfunktion f(x). Dann heißen y_i - f(x_i) (i = 1, ..., n) die ## Residuen. ########################################################### ########################################################### ## Bei den Residuen handelt es sich um den vertikalen Abstand zwischen dem ## Punkt und der Regressionsgerade (Vorhersage) unter Berücksichtigung des ## Vorzeichens. Das Vorzeigen sagt uns, ob der Punkt ober- (+) oder unterhalb (-) ## der Geraden liegt. ########################################################### attach(thuesen) plot(short.velocity ~ blood.glucose, ylab = "fasting blood glucose [mmol/l]", xlab = "mean circumferential shortening velocity [%/s]", main = "Ventricular shortening velocity - residuals", pch = 20) abline(fit1, lwd = 2) segments(blood.glucose[-16], fitted(fit1), blood.glucose[-16], short.velocity[-16]) ## Sind die Voraussetzungen für das Lineare Modell erfüllt, so sollten die ## Residuen die Eigenschaften der stochastischen Fehler im Linearen Modell ## besitzen. ## 1. Homoskedastizität (Variabilität von y unabhängig von x) ## 2. Identisch verteilt ## 3. Unkorreliert bzw. stochastisch unabhängig ## 4. Normalverteilt (für Testtheorie) ####################################### ## Im Idealfall sollte man also das folgende homogene Bild erhalten, bei dem ## die Residuen keinerlei Struktur zeigen. plot(rnorm(1000), pch = 20, ylab = "Residuen", main = "Modell korrekt") abline(h = c(-2, 2), lty = 2) ## Wir haben bewusst viele Zufallszahlen erzeugt, damit die homogene Struktur ## besser sichtbar wird. In der Praxis hat man es jedoch häufig mit deutlich ## weniger Daten zu tun. par(mfrow = c(3, 3)) for(i in 1:9) plot(rnorm(30), pch = 20, ylab = "Residuen", main = "Modell korrekt") ## Wir plotten nun Residuen für den Fall, dass Heteroskedastizität vorliegt plot(1:1000, sqrt(1:1000)*rnorm(1000), ylab = "Residuen", main = "Heteroskedastizität", pch = 20) ## Wir betrachten dieselbe Situation jetzt aber mit weniger Daten par(mfrow = c(3, 3)) for(i in 1:9) plot(1:30, sqrt(1:30)*rnorm(30), ylab = "Residuen", main = "Heteroskedastizität", pch = 20) ## Jetzt kann es schon deutlich schwieriger sein, die Heteroskedastizität zu erkennen. ## Um diese bei wenigen Beobachtungen rel. sicher erkennen zu können, muss sie ## deutlicher ausgeprägt sein. par(mfrow = c(3, 3)) for(i in 1:9) plot(1:30, 1:30*rnorm(30), ylab = "Residuen", main = "Heteroskedastizität", pch = 20) ####################################### ## Ein weiter Punkt, der anhand der Residuen untersucht werden kann, ist, ob das ## Modell nichtlineare Komponenten enthält plot(1:1000, cos((1:1000)*pi/500) + rnorm(1000), ylab = "Residuen", main = "Nichtlinearität", pch = 20) ## Auch dies wieder bei weniger Beobachtungen schwieriger zu erkennen. par(mfrow = c(3, 3)) for(i in 1:9) plot(1:30, cos((1:30)*pi/15) + rnorm(30), ylab = "Residuen", main = "Nichtlinearität", pch = 20) ## Um diese sicher zu erkennen, muss sie bei weniger Beobachtungen stärker ## ausgeprägt sein. par(mfrow = c(3, 3)) for(i in 1:9) plot(1:30, 2*cos((1:30)*pi/15) + rnorm(30), ylab = "Residuen", main = "Nichtlinearität", pch = 20) ####################################### ## Die (Auto-)Korrelation der Residuen läßt sich graphisch mit Hilfe der Funktion ## acf untersuchen acf(rnorm(50), lwd = 2, main = "Korrelation der Residuen") ## Liegen korrelierte Daten vor, so erhalten wir x <- rnorm(500) y <- rnorm(500) r1 <- 0.4 ## vorgegebener Korrelationskoeffizient a1 <- sqrt((1-r1)/(1+r1)) x1 <- x + a1*y y1 <- x - a1*y X <- as.vector(t(cbind(x1, y1))) acf(X, lwd = 2, main = "Korrelation der Residuen") ## Bei weniger Daten wird es wieder deutlich schwieriger par(mfrow = c(3, 3)) for(i in 1:9){ x <- rnorm(15) y <- rnorm(15) r1 <- 0.4 ## vorgegebener Korrelationskoeffizient a1 <- sqrt((1-r1)/(1+r1)) x1 <- x + a1*y y1 <- x - a1*y X <- as.vector(t(cbind(x1, y1))) acf(X, lwd = 2, main = "Korrelation der Residuen") } ## Die Korrelation muss stärker ausgeprägt sein par(mfrow = c(3, 3)) for(i in 1:9){ x <- rnorm(15) y <- rnorm(15) r1 <- 0.8 ## vorgegebener Korrelationskoeffizient a1 <- sqrt((1-r1)/(1+r1)) x1 <- x + a1*y y1 <- x - a1*y X <- as.vector(t(cbind(x1, y1))) acf(X, lwd = 2, main = "Korrelation der Residuen") } ####################################### ## Normalverteilung der Residuen kann man graphisch mittels dem ## qq-Plot überprüfen. X <- rt(1000, df = 10) qqnorm(X) qqline(X) ## Bei weniger Beobachtungen fällt die Entscheidung wieder deutlich schwerer par(mfrow = c(3, 3)) for(i in 1:9){ X <- rt(30, df = 10) qqnorm(X) qqline(X) } ## Wir verstärken die Nichtnormalität par(mfrow = c(3, 3)) for(i in 1:9){ X <- rt(30, df = 2) qqnorm(X) qqline(X) } ####################################### ## Wir sehen uns nun die entsprechenden Residuenplots für den Blutzucker-Datensatz an. plot(fit1$resid, pch = 20, ylab = "Residuen", main = "Ventricular shortening velocity - Residuen") acf(fit1$resid, lwd = 2, main = "Korrelation der Residuen") library(PASWR) ntester(fit1$resid) ## Es sind keine deutlichen Auffälligkeiten zu erkennen. Am auffälligsten erscheint ## noch der qq-Plot. ####################################### ## Natürlich sind auch in R bereits derartige diagnostische Plots eingebaut. ## Es gibt hierfür eine Methode der Funktion "plot" für die Klasse "lm". ## Wir verwenden weiter die Blutzucker-Daten plot(fit1, which = 1) ## Hat man multivariate Regressoren, stellt sich die Frage, was man bei den ## Residuenplots auf der x-Achse abträgt. Üblicherweise wählt man die ## gefitteten Werte. ## Im vorliegenden Fall gibt es keine deutlichen Auffälligkeiten (Strukturen). ## Der zweite Plot ist ein qq-Plot der standardisierten Residuen plot(fit1, which = 2) ## Im dritten Plot ist die Wurzel des Absolutbetrags der Residuen gegen die ## gefitteten Werte abgetragen. plot(fit1, which = 3) ## Im vierten Plot ist der sog. Cook's Abstand abgetragen. Dieser Plot dient ## dazu potentielle Ausreißer zu identifizieren. plot(fit1, which = 4) ## Auch die Plots 5 und 6 dienen dazu Ausreißer zu identifizieren sowie sog. ## Hebelpunkte; d.h., Punkte die eine starken Einfluss auf die geschätzte ## Regression haben und somit die Schätzung "aushebeln" können. plot(fit1, which = 5) plot(fit1, which = 6) ########################################################### ## Es gibt natürlich auch eine Vielzahl von Tests, mit denen sich die ## Voraussetzungen an das Lineare Modell überprüfen lassen. Diese Test ## erkennen in der Regel aber nur bestimmte, spezielle Abweichungen, weshalb ## aus meiner Sicht die Residuenplots den Tests vorzuziehen sind. ## Beispiele solcher Tests sind: ## 1. Heteroskedastizität: ## - studentized Breusch-Pagan Test ## - Goldfeld-Quandt Test ## - Harrison-McCabe Test ## 2. Linearität: ## - Harvey-Collier Test ## - Rainbow Test ## 3. Korrelation: ## - Breusch-Godfrey Test ## - Durbin-Watson Test ## 4. Normalverteilungstests: siehe TestenUndSchätzen.R ########################################################### ########################################################### ## Robuste Regressionschätzer ## Am Beispiel das Blutzucker-Datensatzes ########################################################### require(MASS) require(robustbase) require(RColorBrewer) palette(brewer.pal(5, "Dark2")) ## M-Schätzer fit.M <- rlm(short.velocity ~ blood.glucose) ## MM-Schätzer fit.MM <- rlm(short.velocity ~ blood.glucose, method = "MM") ## LTS-Schätzer fit.LTS <- ltsReg(short.velocity ~ blood.glucose) ## AL-Schätzer ## regressor distribution: design measure require(RobRex) K <- DiscreteMVDistribution(cbind(1,blood.glucose)) ## ALc-Schätzer system.time(IC1 <- rgsOptIC.ALc(r = 0.5, K = K, theta = fit.M$coeff, scale = fit.M$s)) est1 <- oneStepEstimator(cbind(1,as.matrix(thuesen)), IC1, c(fit.M$coeff, fit.M$s)) ## AL-Schätzer system.time(IC2 <- rgsOptIC.AL(r = 0.5, K = K, theta = fit.LTS$coeff, scale = fit.LTS$s)) est2 <- oneStepEstimator(cbind(1,as.matrix(thuesen)), IC2, c(fit.LTS$coeff, fit.LTS$s)) ## Plot der verschiedenen Regressionsschätzer plot(short.velocity ~ blood.glucose, ylab = "fasting blood glucose [mmol/l]", xlab = "mean circumferential shortening velocity [%/s]", main = "Ventricular shortening velocity", pch = 19) abline(fit1, lwd = 2, col = "black") abline(fit.M, lwd = 2, col = 1) abline(fit.MM, lwd = 2, col = 3) abline(fit.LTS, lwd = 2, col = 2) lines(c(1, c(blood.glucose,25)), est1[1] + est1[2]*c(1,c(blood.glucose,25)), col = 4, lwd = 2) lines(c(1, c(blood.glucose,25)), est2[1] + est2[2]*c(1,c(blood.glucose,25)), col = 5, lwd = 2) legend("topleft", legend = c("KQ", "M", "MM", "LTS", "AL1", "AL2"), fill = c("black", 1, 3, 2, 4, 5), ncol = 2) ########################################################### ## Robuste Regressionschätzer ## Am Beispiel: Jährliche Anzahl von Anrufen ########################################################### data(phones) attach(phones) plot(calls~year, pch = 20) fit2 <- lm(calls ~ year) ## M-Schätzer fit2.M <- rlm(calls ~ year, maxit = 50) ## MM-Schätzer fit2.MM <- rlm(calls ~ year, maxit = 50, method = "MM") ## LTS-Schätzer fit2.LTS <- ltsReg(calls ~ year) ## AL-Schätzer ## regressor distribution: design measure K <- DiscreteMVDistribution(cbind(1,year)) ## ALc-Schätzer system.time(IC1 <- rgsOptIC.ALc(r = 0.5, K = K, theta = fit2.M$coeff, scale = fit2.M$s)) est1 <- oneStepEstimator(cbind(1,year,calls), IC1, c(fit2.M$coeff, fit2.M$s)) ## AL-Schätzer system.time(IC2 <- rgsOptIC.AL(r = 0.5, K = K, theta = fit2.LTS$coeff, scale = fit2.LTS$s)) est2 <- oneStepEstimator(cbind(1,year,calls), IC2, c(fit2.LTS$coeff, fit2.LTS$s)) ## Plot der verschiedenen Regressionsschätzer plot(calls ~ year, ylab = "phone calls [Mio.]", xlab = "year", main = "Belgium Phone Calls 1950-1973", pch = 20) abline(fit2, lwd = 2, col = "black") abline(fit2.M, lwd = 2, col = 1) abline(fit2.MM, lwd = 2, col = 3) abline(fit2.LTS, lwd = 2, col = 2) lines(c(1, c(year,75)), est1[1] + est1[2]*c(1,c(year,75)), col = 4, lwd = 2) lines(c(1, c(year,75)), est2[1] + est2[2]*c(1,c(year,75)), col = 5, lwd = 2) legend("topleft", legend = c("KQ", "M", "MM", "LTS", "AL1", "AL2"), fill = c("black", 1, 3, 2, 4, 5), ncol = 2)