Eine Regression in R ist vielleicht etwas ungewohnt, dafür liefert diese in kürzester Zeit Regressionen für jedes nur erdenkliche Modell und gibt mit nur wenigen Befehlen Statistiken zu den Residuen aus.
Für die Regression laden wir die Daten von Anscombe:
> attach(anscombe) # die Daten in den Suchpfad aufnehmen > anscombe x1 x2 x3 x4 y1 y2 y3 y4 1 10 10 10 8 8.04 9.14 7.46 6.58 2 8 8 8 8 6.95 8.14 6.77 5.76 3 13 13 13 8 7.58 8.74 12.74 7.71 4 9 9 9 8 8.81 8.77 7.11 8.84 ... > par.backup <- par(mfrow=c(2,2)) > par.backup $mfrow [1] 1 1 > for(i in 1:4) { + x <- get( paste("x", i, sep="")) + y <- get( paste("x", i, sep="")) + plot(x, y, main=paste("Plot x", i, " gegen y", i, sep="" )) + }
Der R-Code lädt die Anscombe Daten, welche 4 verschiedene Messreihen enthalten. Als nächstes wird die Grafik in 4 Subgrafiken unterteilt und die alten Einstellungen in par.backup gespeichert. Darauf folgt eine gezählte Schleife, welche mit for(){} aufgerufen wird. Innerhalb der runden Klammern steht dabei die Variable, die gezählt werden soll. In diesem Beispiel wird die Variable i von 1 auf 4 gezählt. Innerhalb der Schleife werden aus dem Datensatz die zusammengehörenden Variablen x1, x2, x3, x4 und y1, y2, y3, y4 in jedem Durchgang in den Variablen x und y gespeichert und abschließend mit der Funktion plot() als Scatterplot dargestellt. Die Funktion paste(Argument1, Argument2,..., sep=) setzt dabei die verschiedenen Argumente mit dem in sep= angegebenen Trennzeichen zusammen.
R besitzt die Möglichkeit, jedes Modell gegen Daten zu fitten. Alles was dafür nötig ist, ist das Modell in einer für R verständlichen Form zu beschreiben. Und das ist auch die einzige Hürde, die es zu meistern gilt! Ein Modell wird dabei mit einer Formel beschrieben, für die es eine eigene Notation gibt. So existieren verschiedene Operatoren, die man innerhalb der Formel verwenden kann:
Ein paar Beispiele sollen das verdeutlichen: Sie wollen folgende Formel als Modell in R darstellen: \[ y = \alpha_0 +\alpha_1x_1 + \alpha_2x_2 + \epsilon \] Die abhängige Variable ist hierbei \( y \) und die unabhängigen Variablen sind \( x_1 \) und \( x_2\). Die Verschiebung des Modells an der y-Achse ist in der Standardeinstellung enthalten und \( \epsilon \) entspricht den Residuen zwischen Modell und Datensatz. Das Modell in Formelnotation in R wäre folgendes:
y ~ x1 + x2
Ist eine Verschiebung auf der y-Achse zum Beispiel nicht erwünscht, so sähe die Formel: \[ y = \alpha_1x_1 + \alpha_2x_2 + \epsilon \] als Modell in R folgender Maßen aus:
y ~ x1 +x2 -1
Die Funktion in R für lineare Regression lautet \verb+lm()+ Die Abbildung zeigt, dass es sich im Plot x1 gegen y1 wahrscheinlich um einen linearen Zusammenhang handelt. Eine lineare Regression nach der Formel: \[ y = \alpha_0 + \alpha_1x + \epsilon \] entspricht dem Modell \verb+y~x+ in R. Folgender Code erzeugt eine lineare Regression:
> reg.lin <- lm(y1~x1, data=anscombe) # lineare Regression > reg.lin Call: lm(formula = y1 ~ x1, data = anscombe) Coefficients: (Intercept) x1 3.0001 0.5001 > summary(reg.lin) Call: lm(formula = y1 ~ x1, data = anscombe) Residuals: Min 1Q Median 3Q Max -1.92127 -0.45577 -0.04136 0.70941 1.83882 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.0001 1.1247 2.667 0.02573 * x1 0.5001 0.1179 4.241 0.00217 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.237 on 9 degrees of freedom Multiple R-squared: 0.6665, Adjusted R-squared: 0.6295 F-statistic: 17.99 on 1 and 9 DF, p-value: 0.002170 > res <- residuals(reg.lin) > res 1 2 3 4 5 6 0.03900000 -0.05081818 -1.92127273 1.30909091 -0.17109091 -0.04136364 7 8 9 10 11 1.23936364 -0.74045455 1.83881818 -1.68072727 0.17945455 > reg.lin.2 <- lm(y1~x1 -1, data=anscombe) > summary(reg.lin.2) Call: lm(formula = y1 ~ x1 - 1, data = anscombe) Residuals: Min 1Q Median 3Q Max -2.7784 -0.5962 0.5756 1.4586 2.4592 Coefficients: Estimate Std. Error t value Pr(>|t|) x1 0.79680 0.04962 16.06 1.81e-08 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.57 on 10 degrees of freedom Multiple R-squared: 0.9627, Adjusted R-squared: 0.9589 F-statistic: 257.9 on 1 and 10 DF, p-value: 1.812e-08 > plot(x1, y1, xlim=c(0,15), ylim=c(0,15)) > abline(reg.lin, col="blue") > abline(reg.lin.2, col="red")
Sie können außerdem noch mit der Funktion plot() verschiedene Grafiken zur Residualanalyse ausgeben lassen:
> par(mfrow=c(2,2)) > plot(reg.lin)
Residual vs Fitted Der klassische Residualplot trägt die Residuen gegen die angepassten Werte auf. Idealerweise sollten die Residuen dabei Null sein, was der gestrichelten Linie entspricht. Die tatsächlichen Residuen sind mit Hilfe der roten Linie eingezeichnet.
Mit Hilfe der Funktion predict() lassen sich aus dem Modell weitere Werte vorhersagen. Dabei benötigt die Funktion als 2. Argument ein data.frame mit den zu verwendenden neuen Werten für die Vorhersage.
> new <- data.frame(x1=seq(-20,20,0.5)) > reg.pred.p <- predict(reg.lin, new, interval="prediction") > reg.pred.c <- predict(reg.lin, new, interval="confidence") > matplot(new$x1, cbind(reg.pred.p, reg.pred.c[,-1]), + lty=c(1,2,2,3,3), col=c(1,2,2,3,3), type="l", ylab="predicted y") > points(x1,y1) > legend(5,-5, legend=c("predicted y", "pred. interval", "confidence + interval"), lty=c(1,2,3), col=c(1,2,3))
Um in R nichtlineare Modelle zu fitten, muss man die Funktion nls() verwenden, welche für nonlinear least squares steht. In der vorherigen Abbildung erkennt man, dass die Daten \( y_2 \) \( x_2 \) einem quadratischen Modell folgen.
> reg.quad <- nls(y2 ~ a*x2^2 + b*x2 + c, data=anscombe, + start=list(a=1, b=1, c=1)) > summary(reg.quad) Formula: y2 ~ a * x2^2 + b * x2 + c Parameters: Estimate Std. Error t value Pr(>|t|) a -0.1267133 0.0000571 -2219 <2e-16 *** b 2.7808392 0.0010401 2674 <2e-16 *** c -5.9957343 0.0043299 -1385 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.001672 on 8 degrees of freedom Number of iterations to convergence: 1 Achieved convergence tolerance: 6.097e-08 >fitted(reg.quad) [1] 9.141329 8.141329 8.740629 8.768042 9.261189 8.100210 6.127622 3.100210 [9] 9.127622 7.261189 4.740629 attr(,"label") [1] "Fitted values"
Im Aufruf von nls() wird als erstes das entsprechende Modell definiert (hier eine klassische Parabel \( y=ax^2+bx+c\) ), gefolgt von den entsprechenden Daten. Zusätzlich benötigt die Funktion eine Liste mit Startwerten für die einzelnen Koeffizienten, welche mit start= übergeben wird. Die Funktion summary() zeigt eine Zusammenfassung der Statistik des Fits an und mit Hilfe der Funktion fitted() lassen sich die mit Hilfe des Modells angepassten Werte ausgeben.
> new <- data.frame(x2=seq(3,15,0.05)) > lm2.pred <- predict(lm2,new) > plot(new$x2, lm2.pred, type="l", col=2) > points(x2, y2)
R besitzt vielfältige Möglichkeiten verschiedene Grafiken zu erstellen. Dabei unterscheidet man ganz allgemein zwischen \emph{high-level} und \emph{low-level}-Funktionen. Dabei erzeugen die high-level Funktionen komplette Grafiken mit Achsenbeschriftung etc. Mit low-level-Funktionen können dagegen verschiedene Elemente zu einer Grafik hinzugefügt werden (z.Bsp. Pfeile, Legenden etc.) Die Grafiken werden dabei auf einem so genannten Device ausgegeben, wo verschiedene Devices zur Auswahl stehen:
Als optionales Argument kann man den Funktionen für die verschiedenen Devices auch einen Dateinamen übergeben: png("Grafik1.png") Außerdem sind Angaben über die Größe und Auflösung möglich. (siehe Hilfe) Nach dem die Grafik fertig erstellt wurde, muss das Grafikdevice noch mit dev.off() geschlossen werden, damit die Bilddateien von anderen Programmen gelesen werden können!
Die verschiedenen Grafikfunktionen besitzen eine Menge zusätzlicher Parameter, mit denen sich das Aussehen der Grafik beeinflussen lässt. So kann man die Achsenbeschriftung, den Titel und die Farben einstellen. Zu den einzelnen Möglichkeiten muss man die Hilfe lesen, hier jedoch ein paar sehr allgemeine Parameter, die meistens unterstützt werden:
par() ist eine sehr mächtige Grafikfunktion, mit der sich alle Einstellungen rund um Grafiken in R ändern lassen. Von den Einstellungen zu den Rändern bis hin zur Unterteilung in Untergrafiken. Die Möglichkeiten sind vielseitig und würden hier den Rahmen sprengen, deshalb sei auf die Hilfe verwiesen: ?par Die häufigste Änderung ist dabei das erzeugen von Untergrafiken. Das Beispiel erzeugt einen Plot mit 4 Untergrafiken:
> par(mfrow = c(2,2)) > for (i in 1:4){ + curve(x^i, from= -3, to=3, lty=i) + }
Der Befehl zum Unterteilen des Grafikfensters in Untergrafiken lautet mfrow= und ist ein Argument der Funktion par(). Dabei wird in einem Vektor die Anzahl der Spalten und Zeilen übergeben. Es bietet sich an, die alten Einstellungen in einem Objekt zu speichern, um sie zurück setzen zu können. Dies ist vor allem nützlich, wenn man mit verschiedenen Argumenten experimentiert:
> par.back <- par(mfrow = c(2,2)) # speichern der alten Einstellungen > par(par.back) # alte Einstellungen # wiederherstellen
> png() > set.seed(10) > x <- rnorm(1000) > hist(x, freq=FALSE, xlim=c(-4,4), ylim=c(0,0.5), ylab="Dichte", main="Dichte einer Standardnormalverteilung", col="grey") > curve(dnorm, from=-4, to=4, add=TRUE, lty=3, lwd=2, col="red") >legend(2, 0.5, legend=c("empirisch", "theoretisch"), + col=c("grey","red"), lwd=2, lty=c(1,3)) > arrows(-2.1, dnorm(-2.1)+0.1, -2.1, dnorm(-2.1), length=0.15, code=2,lwd=2) > text(-2.1,dnorm(-2.1)+0.11, labels=signif(dnorm(-2.1),3)) > text(2, 0.4, adj=0, cex=1, + expression(f(x)==frac(1,sigma*sqrt(2*pi)) ~~ + e^{frac(-(x-mu)^2,2*sigma^2)})) > dev.off()