Lineare Regression mit Hilfe von lm() und Funktionen zum Extrahieren von Modelleigenschaften

Einordnung des Artikels

• Einführung in die Informatik
  • Einführung in die Programmiersprache R
    • Statistik
      • Lineare Regression mit Hilfe von lm() und Funktionen zum Extrahieren von Modelleigenschaften

In Demonstration zur Methode der kleinsten Quadrate und der Eigenschaften der Regressionsgerade und Regressionsanalyse: Die Varianzzerlegung, das Bestimmtheitsmaß und der Residualplot wurden zwei Beispiele zur Regressionsanalyse angegeben, einmal für stark und einmal für schwach korrelierte Messdaten. Die stark korrelierten Testdaten werden hier ebenfalls eingesetzt; in Letzterem der genannten Artikel werden die Testdaten kurz zusammengefasst.

Die Begriffe und Formeln aus der Statistik, die in den genannten Artikeln eingesetzt werden, werden hier nicht erläutert sondern vorausgesetzt.

Einführung

In Demonstration zur Methode der kleinsten Quadrate und der Eigenschaften der Regressionsgerade und Regressionsanalyse: Die Varianzzerlegung, das Bestimmtheitsmaß und der Residualplot wird mit Hilfe einfacher Testdaten gezeigt, wie man Steigung und y-Abschnitt der Regressionsgerade, weitere Zufallsvariablen, nämlich Regressionswert und Residuum, und das Bestimmtheitsmaß berechnet. Dort wurden weder Herleitungen der Formeln gegebenen (meist nur angedeutet) noch Quelltexte gezeigt, wie man die entsprechenden Größen berechnet.

Hier werden diese Testdaten verwendet, um zu zeigen, wie man mit Hilfe der R-Funktion lm() aus dem Basis-Paket stats die soeben genannten Größen berechnen und einsetzen kann. Dabei wird an einigen Stellen angedeutet, dass die Funktion lm() (und mit ihr verwandte Funktionen) Antworten auf Fragen geben kann, die weit darüber hinausgehen. Ausführlich besprochen werden aber nur die statistischen Konzepte aus den beiden genannten Artikeln.

Der Name der Funktion lm() sowie der Name der Klasse "lm" ihres Rückgabewertes leiten sich von linear model ab.

Vorbereitung: Erzeugen von Testdaten

Das folgende Skript zeigt eine einfache Möglichkeit, wie Testdaten für eine Regressionsanalyse erzeugt werden können. Ausführlich wird diese Methode in Demonstration zur Methode der kleinsten Quadrate und der Eigenschaften der Regressionsgerade beschrieben, hier erfolgt nur eine kurze Beschreibung. Es werden zuerst Messpunkte erzeugt, die auf der Gerade y = m*x liegen, anschließend werden die y-Koordinaten durch einen Zufallsgenerator leicht variiert.

Zeile 1: Damit die Regressionsanalyse reproduzierbar ist, wird der Zufallsgenerator mit einer festen Zahl initialisiert.

Zeile 2: Die Steigung m der Ausgangsgerade wird mit 2 festgelegt. Die Steigung der Regressionsgerade wird dann einen ähnlichen Wert besitzen.

Zeile 3: Es werden 20 x-Werte in äquidistantem Abstand festgelegt.

Zeile 4: Die Funktion jitter() sorgt dafür, dass die y-Werte leicht von m*xs abweichen; dazu wird von jitter() eine Gleichverteilung wie in runif(n, -amount, amount) verwendet.

set.seed(seed = 42)
m <- 2
xs <- (1:20)
ys <- jitter(m*xs, amount = 2)

Diese Testdaten werden unter den Namen xs und ys in allen folgenden Skripten verwendet und nicht nochmals erklärt.

Die wichtigsten Eingabewerte von lm()

Sämtliche Eingabewerte von lm() werden hier nicht besprochen; sie sind dem folgenden Listing zu entnehmen:

lm(formula, data, subset, weights, na.action,
    method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
    singular.ok = TRUE, contrasts = NULL, offset, ...)

Die Eingabewerte formula und data

Der wichtigste Eingabewert der Funktion lm() ist formula als eine symbolische Beschreibung des Modells, für das die Regressionsgerade (und weitere Eigenschaften) berechnet werden soll. Diese symbolische Beschreibung kann auf zwei Arten realisiert werden:

1. Die y- und x-Werte der Messdaten werden unmittelbar mit dem Operator ~ verknüpft.
2. Die Messdaten werden in einem Dataframe abgespeichert; der Operator ~ verknüpft dann die Spalten des Dataframes.

In den beiden folgenden Unterabschnitten werden diese zwei Möglichkeiten gezeigt.

Festlegen des Modells allein mit formula

Das folgende Skript realisiert die erste Möglichkeit. Es wirkt sehr lang; zum Verständnis des Eingabewertes formula sind vorerst nur die Zeilen 1 und 7 nötig, ansonsten werden übliche Funktionen zum Erzeugen eines geeigneten Plots aufgerufen (siehe Abbildung 1).

mod <- lm(formula = ys ~ xs)

plot(x = xs, y = ys, col = "red",
     xlim = c(0, 20),
     ylim = c(min(0, min(ys)), max(ys)),
     xlab = "x", ylab = "y")
abline(mod, col = "blue")
grid()
# Koord.achsen:
abline(h = 0, col = "darkgray")
abline(v = 0, col = "darkgray")

Zur Erklärung:

Zeile 1: In Zeile 1 müssen drei Dinge erklärt werden:

1. Was bedeutet die Symbolik ys ~ xs im Argument formula ?
2. Warum dürfen ys und xs in ys ~ xs nicht vertauscht werden?
3. Welches Ergebnis liefert der Aufruf von lm() und wie kann man dieses Objekt einsetzen?

Zu 1.: Mit xs und ys werden die oben erzeugten Testdaten eingesetzt. Sie werden in lm() im Argument formula durch das Symbol ~ miteinander verknüpft. Dadurch entsteht ein Objekt der Klasse "formula" . Da die Funktion lm() ein lineares Modell erzeugt, verbirgt sich hinter ys ~ xs die Anweisung, dass die gegebenen x- und y-Werte durch eine lineare Funktion y = a*x + b verknüpft werden sollen – und falls die Messpunkte nicht auf einer eindeutigen Gerade liegen, wird mit der Methode der kleinsten Quadrate die Regressionsgerade bestimmt.

Zu 2.: Da die Methode der kleinsten Quadrate die x- und y-Werte nicht gleichberechtigt behandelt, muss die Reihenfolge in ys ~ xs beachtet werden. Die Fehlerfunktion Q(a, b) der Methode der kleinsten Quadrate nimmt die x-Werte als fehlerfrei an und berechnet die Summe der quadratischen Abweichungen der y-Werte von der Geraden y = ax + b. Entsprechend ist ys ~ xs wie die funktionale Abhängigkeit y = y(x) zu lesen, in der x die unabhängigen und y die abhängigen Werte sind.

Zu 3.: Das Ergebnis des Funktionsaufrufs von lm() wird in der Variable mod abgespeichert (dabei soll der Name mod nur andeuten, dass ein Modell entsteht). Genauer handelt es sich um ein Objekt der Klasse "lm" , das weiter unten näher untersucht wird. Salopp ist darin das gesamte lineare Modell enthalten, das die Regressionsanalyse liefert; auf die Eigenschaften des Modells kann durch weitere Funktionsaufrufe zugegriffen werden. Vorerst wird das Objekt mod nur in Zeile 7 verwendet.

Zeile 3 bis 11: Die Graphik-Anweisungen werden hier nicht im Detail erläutert. Die Messpunkte werden rot als Punkte dargestellt; dies geschieht durch den Aufruf von plot(), wo auch weitere Einstellungen für das Diagramm vorgenommen werden. Mit Hilfe von abline() werden die beiden Koordinatenachsen (Zeile 10 und 11) und die Regressionsgerade (Zeile 7) dem Plot hinzugefügt.

Zeile 7: Die Funktion abline() wird meist mit den Argumenten a und b eingesetzt, die für den y-Abschnitt a und die Steigung b einer Gerade stehen (Achhtung: also vertauscht gegenüber der Verwendung y = ax + b, mit der meist die Geraden bezeichnet werden). Man kann anstelle von a und b aber auch ein Argument reg setzen, das ein Regressions-Objekt erhält, das wiederum eine Funktion coef() besitzen muss, um die Regressionskoeffizienten zu extrahieren. Wie weiter unten gezeigt wird, besitzt das von lm() erzeugte Objekt genau diese Methode. Damit sorgt Zeile 7 dafür, dass die Regressionsgerade in blauer Farbe dem Plot hinzugefügt wird.

Abbildung 1 zeigt den Plot, den obiges Skript erzeugt.

Abbildung 1: Darstellung der Testdaten (rote Punkte) und der von lm() berechneten Regressionsgerade (blau). Der Plot wurde mit den Anweisungen aus dem letzten Skript erstellt.

Ablegen der Daten in einem Dataframe

Das folgende Skript realisiert die zweite Möglichkeit, wie formula eingesetzt wird, wenn die Daten in einem Dataframe (Eingabewert data) abgelegt sind. Dies sollte man als den typischen Fall betrachten, denn oben wurden die Testdaten eigens erzeugt, um damit die Regressionsanalyse durchzuführen. Typischerweise liegen die auszuwertenden Daten bereits in einem Dataframe vor und sie werden von dort an die Regressionsanalyse weitergereicht.

# xs, ys wie oben

df.xy <- data.frame(x = xs, y = ys)
str(df.xy)
# 'data.frame': 20 obs. of  2 variables:
# $ x: int  1 2 3 4 5 6 7 8 9 10 ...
# $ y: num  3.66 5.75 5.14 9.32 10.57 ...

mod <- lm(formula = y ~ x, data = df.xy)

plot(x = df.xy$x, y = df.xy$y, col = "red",
     xlim = c(0, 20),
     ylim = c(min(0, min(ys)), max(ys)),
     xlab = "x", ylab = "y",
     main = "Messpunkte und Regressionsgerade")
abline(reg = mod, col = "blue")
grid()
# Koord.achsen:
abline(h = 0, col = "darkgray")
abline(v = 0, col = "darkgray")

Das Skript ist lediglich eine leichte Abwandlung gegenüber oben; erklärt werden nur die relevanten Unterschiede:

Zeile 3: Es wird ein Dataframe df.xy mit zwei Spalten gebildet; die Spalten besitzen die Namen x beziehungsweise y und beinhalten die Testdaten xs und ys aus dem vorherigen Skript.

Zeile 4 bis 7: Mit str(df.xy) wird das Dataframe untersucht.

Zeile 9: Wie oben wird das von lm() erzeugte Objekt abgespeichert. Neu ist der Aufruf von lm():

1. Mit dem Argument data = df.xy wird das Dataframe angesprochen.
2. Die zu verknüpfenden Variablen werden mit formula = y ~ x festgelegt; damit werden die Spalten des Dataframes df.xy bezeichnet, wodurch wieder die bekannten Testdaten geladen werden.

Zeile 11: Für den Plot werden die x- und y-Werte aus dem Dataframe df.xy geladen.

Der mit dem Skript erzeugte Plot ist identisch mit Abbildung 1.

Auswahl einer Teilmenge mit subset

Oft sollen nicht alle Daten zur Auswertung verwendet werden – etwa wenn "Ausreißer" entfernt werden. Damit man die entsprechende Teilmenge nur beim Aufruf von lm() angeben muss und nicht ein eigenes (bereinigtes) Dataframe erzeugen muss, kann die Teilmengenbildung mit dem Argument subset erfolgen.

Im folgenden Skript wird das Dataframe df.xy aus dem letzten Skript verwendet, aber durch die Angabe von subset im Aufruf von lm() werden nur bestimmte Zeilen des Dataframes df.xy ausgewählt (Zeile 3). Dazu erinnere man sich daran, wie die Testdaten erzeugt wurden: man ist von xs <- (1:20) und ys <- 2 * xs ausgegangen, anschließend wurden die Werte ys durch einen Zufallsgenerator leicht variiert. Jetzt werden in subset diejenigen Indizes bestimmt, für die der y-Wert um weniger als 0.5 von 2*x abweicht

Welche Zeilen dies sind, kann man durch die explizite Berechnung und Ausgabe der zugehörigen Indizes leicht feststellen (Zeile 5).

# df.xy wie oben

mod <- lm(formula = y ~ x, data = df.xy, subset = which( abs(df.xy$y - 2 * df.xy$x) < 0.5) )

which( abs(df.xy$y - 2 * df.xy$x) < 0.5)
# [1]  6 11 15 19 20

Abbildung 2 zeigt den zugehörigen Plot, in dem die 5 Messpunkte, die die Bedingung im Argument subset aus Zeile 3 erfüllen, rot gekennzeichnet sind. Beim Vergleich mit Abbildung 1 fällt auf, dass die Regressionsgerade jetzt näher an der Gerade y = 2*x liegt. Wie man die Regressionskoeffizienten aus dem Modell mod extrahieren kann, wird im nächsten Abschnitt gezeigt.

Abbildung 2: Werden für die Regressionsanalyse nur die mit subset charakterisierten Messpunkte ausgewählt (hier die rot gekennzeichneten Punkte), entsteht eine leicht veränderte Regressionsgerade (leicht erkennbar, wenn man die y-Abschnitte in Abbildung 1 und 2 vergleicht).

Der Rückgabewert von lm()

Die Klasse lm

Der Rückgabewert der Funktion lm() ist ein Objekt der Klasse "lm" ; dieses Objekt wird hier kurz vorgestellt, aber man sollte beim Umgang mit der Funktion lm() zwei Dinge beachten:

1. Im Rückgabewert von lm() sind alle Informationen enthalten, die man bei der Berechnung einer Regressionsgerade erwartet. Allerdings ist das Objekt sehr kompliziert aufgebaut, so dass der Umgang damit sehr fehleranfällig ist.
2. Um auf die Informationen zuzugreifen, die man bei der Regressionsanalyse am Häufigsten benötigt, gibt es Hilfsfunktionen, die eine einfache Struktur besitzen: Sie erhalten das Objekt der Klasse "lm" als Eingabewert und liefern als Rückgabewert meist einen Vektor als Rückgabewert, mit dem man leichter weiterrechnen kann als mit dem Objekt der Klasse "lm" .

Das folgende Skript bildet wie oben durch den Aufruf von lm() das Modell mod für die gezeigten Testdaten. Das Objekt mod wird kurz untersucht:

# mod wie in den Beispielen oben:
mod <- lm(formula = y ~ x, data = df.xy)

is.list(mod)
# [1] TRUE
class(mod)
# [1] "lm"

Man erkennt, dass das Objekt mod eine Liste ist, für die sogar das Klassenattribut gesetzt ist, nämlich "lm" .

Übersicht mit names()

Um einen schnellen Überblick zu gewinnen, welche Informationen in dem Objekt mod aus dem letzten Skript enthalten sind, kann man sich mit der Funktion names() anzeigen lassen, welche Attribute das Objekt besitzt, die über einen Namen angesprochen werden können:

mod <- lm(formula = y ~ x, data = df.xy)

names(mod)
# [1] "coefficients"  "residuals"     "effects"       "rank"          "fitted.values" "assign"        "qr"           
# [8] "df.residual"   "xlevels"       "call"          "terms"         "model"

Man erkennt, dass 12 Namen gesetzt sind, von denen nur "coefficients" , "residuals" und "fitted.values" besprochen werden; sie reichen aus, um die Plots zu erzeugen, die in Demonstration zur Methode der kleinsten Quadrate und der Eigenschaften der Regressionsgerade gezeigt wurden.

Möchte man das Objekt mod genauer untersuchen, muss man dessen Struktur anzeigen (Zeile 4):

# mod wie oben:
mod <- lm(formula = y ~ x, data = df.xy)

str(mod)
# List of 12
# $ coefficients : Named num [1:2] 0.827 1.964
# ..- attr(*, "names")= chr [1:2] "(Intercept)" "x"
# $ residuals    : Named num [1:20] 0.8678 0.9926 -1.5755 0.6374 -0.0818 ...
# ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
# $ effects      : Named num [1:20] -95.939 50.655 -1.882 0.357 -0.337 ...
# ..- attr(*, "names")= chr [1:20] "(Intercept)" "x" "" "" ...
# $ rank         : int 2
# $ fitted.values: Named num [1:20] 2.79 4.76 6.72 8.68 10.65 ...
# ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
# $ assign       : int [1:2] 0 1
# $ qr           :List of 5
# ..$ qr   : num [1:20, 1:2] -4.472 0.224 0.224 0.224 0.224 ...
# .. ..- attr(*, "dimnames")=List of 2
# .. .. ..$ : chr [1:20] "1" "2" "3" "4" ...
# .. .. ..$ : chr [1:2] "(Intercept)" "x"
# .. ..- attr(*, "assign")= int [1:2] 0 1
# ..$ qraux: num [1:2] 1.22 1.26
# ..$ pivot: int [1:2] 1 2
# ..$ tol  : num 1e-07
# ..$ rank : int 2
# ..- attr(*, "class")= chr "qr"
# $ df.residual  : int 18
# $ xlevels      : Named list()
# $ call         : language lm(formula = y ~ x, data = df.xy)
# $ terms        :Classes 'terms', 'formula'  language y ~ x
# .. ..- attr(*, "variables")= language list(y, x)
# .. ..- attr(*, "factors")= int [1:2, 1] 0 1
# .. .. ..- attr(*, "dimnames")=List of 2
# .. .. .. ..$ : chr [1:2] "y" "x"
# .. .. .. ..$ : chr "x"
# .. ..- attr(*, "term.labels")= chr "x"
# .. ..- attr(*, "order")= int 1
# .. ..- attr(*, "intercept")= int 1
# .. ..- attr(*, "response")= int 1
# .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
#   .. ..- attr(*, "predvars")= language list(y, x)
# .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
# .. .. ..- attr(*, "names")= chr [1:2] "y" "x"
# $ model        :'data.frame': 20 obs. of  2 variables:
#   ..$ y: num [1:20] 3.66 5.75 5.14 9.32 10.57 ...
# ..$ x: int [1:20] 1 2 3 4 5 6 7 8 9 10 ...
# ..- attr(*, "terms")=Classes 'terms', 'formula'  language y ~ x
# .. .. ..- attr(*, "variables")= language list(y, x)
# .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
# .. .. .. ..- attr(*, "dimnames")=List of 2
# .. .. .. .. ..$ : chr [1:2] "y" "x"
# .. .. .. .. ..$ : chr "x"
# .. .. ..- attr(*, "term.labels")= chr "x"
# .. .. ..- attr(*, "order")= int 1
# .. .. ..- attr(*, "intercept")= int 1
# .. .. ..- attr(*, "response")= int 1
# .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
#   .. .. ..- attr(*, "predvars")= language list(y, x)
# .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
# .. .. .. ..- attr(*, "names")= chr [1:2] "y" "x"
# - attr(*, "class")= chr "lm"

Man kann hierin die Komponenten identifizieren, die oben mit names() angezeigt wurden; wie schon gesagt, ist es nicht nötig, das von lm() erzeugte Objekt mod zu verarbeiten. Wenn das Modell einmal erzeugt wurde, kann man mit ihm die bereitgestellten Hilfsfunktionen aufrufen, die in den folgenden Abschnitten erklärt werden.

Die Berechnung der Regressionskoeffizienten mit coefficients()

Die Funktion coefficients() befindet sich wie lm() im Basis-Paket stats. Anstelle des Namens coefficients() kann auch die Abkürzung coef() verwendet werden.

Im folgenden Skript wird gezeigt, wie man mit Hilfe der Funktion coefficients() die beiden Regressionskoeffizienten erzeugen und einsetzen kann:

mod <- lm(formula = y ~ x, data = df.xy)

cf <- coefficients(object = mod)

str(cf)
# Named num [1:2] 0.827 1.964
# - attr(*, "names")= chr [1:2] "(Intercept)" "x"

cf
# (Intercept)           x 
# 0.8270583   1.9643359 

# Plot der Regressionsgerade:
plot(x = df.xy$x, y = df.xy$y,
     col = "red", type = "p", 
     xlim = c(-1, 21), ylim = c(-1, 41),
     xlab = "x", ylab = "y",
     main = paste0( "Messpunkte und Regr.gerade ax + b, mit a = ", format(cf[2], digits = 4),
                   ", b = ", format(cf[1], digits = 4) ),
     frame.plot = TRUE)
grid()
abline(a = cf[1], b = cf[2], col = "blue")
# Koord.achsen:
abline(h = 0, col = "darkgray")
abline(v = 0, col = "darkgray")

Zeile 1: Mit lm() wird wieder aus den Testdaten das Objekt mod erzeugt.

Zeile 3: An die Funktion coefficients() wird das Objekt mod übergeben; das Ergebnis des Funktionsaufrufs wird in cf abgespeichert.

Zeile 5 bis 7: Man erkennt, dass cf ein Vektor der Länge 2 ist; die Namen der Komponenten zeigen ihre Bedeutung:

• die erste Komponente ist der y-Abschnitt,
• die zweite Komponente ist die Steigung der Regressionsgerade.

Hier kann vielleicht eine Merkregel hilfreich sein, da die Reihenfolge der Koeffizienten vielleicht ungewohnt erscheint:

• Üblicherweise werden Polynome nach Potenzen von x sortiert, so dass man die Regressionsgerade in der Form y = ax + b angibt, also zuerst die Steigung a und zuletzt den y-Abschnitt b.
• Der Rückgabewert von coefficients() vertauscht die Reihenfolge der Koeffizienten. Sie werden so angegeben, wie man sie zum Zeichnen einsetzt: Man muss in ein Diagramm zuerst den y-Abschnitt b einzeichnen und kann dann die Steigung a eintragen. (Nebenbei: auch die Funktion abline() verwendet y-Abschnitt und Steigung in dieser Reihenfolge.)

Zeile 9 bis 11: Die Ausgabe von cf zeigt die Namen und Werte der Regressionskoeffizienten.

Zeile 13 bis 25: Mit den Anweisungen wird der Plot in Abbildung 3 erzeugt. Der Zugriff auf die Regressionskoeffizienten erfolgt dabei in den Zeilen 18, 19 und 22 für die Überschrift des Diagramms sowie für das Eintragen der Regressionsgerade mit Hilfe von abline().

Abbildung 3: Darstellung der Testdaten (rot) und der Regressionsgerade (blau) unter Verwendung der mit coefficients() berechneten Regressionskoeffizienten.

Die Berechnung der Regressionswerte mit fitted.values()

Projiziert man die Messpunkte der Testdaten entlang der y-Richtung auf die Regressionsgerade, so erhält man die Regressionswerte W; genauer die y-Koordinaten der projizierten Messpunkte sind die Regressionswerte:

W = aX + b,

wobei für a und b die Regressionskoeffizienten eingesetzt werden.

Den Vektor der Regressionswerte kann man mit Hilfe der Funktion fitted.values() berechnen. (Die Kurzversion lautet hier fitted().) Das folgende Skript zeigt diese Berechnung und den typischen Einsatz der Regressionswerte für Abbildung 4:

# mod und cf wie oben

regr.values <- fitted.values(object = mod)

str(regr.values)
# Named num [1:20] 2.79 4.76 6.72 8.68 10.65 ...
# - attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...

par(mfrow = c(1, 2))

plot(x = df.xy$x, y = df.xy$y,
     col = "red", type = "p", 
     xlim = c(-1, 21), ylim = c(-1, 41),
     xlab = "x", ylab = "y",
     main = paste0( "Messpunkte und Regr.gerade ax + b, mit a = ", format(cf[2], digits = 4),
                    ", b = ", format(cf[1], digits = 4) ),
     frame.plot = TRUE)
grid()
abline(a = cf[1], b = cf[2], col = "blue")
# Koord.achsen:
abline(h = 0, col = "darkgray")
abline(v = 0, col = "darkgray")

plot(x = df.xy$x, y = regr.values,
     col = "orange", type = "p", 
     lwd = 2,
     xlim = c(-1, 21), ylim = c(-1, 41),
     xlab = "x", ylab = "y",
     main = paste0( "Regressionswerte und Regr.gerade ax + b, mit a = ", format(cf[2], digits = 4),
                    ", b = ", format(cf[1], digits = 4) ),
     frame.plot = TRUE)
grid()
abline(a = cf[1], b = cf[2], col = "blue")
# Koord.achsen:
abline(h = 0, col = "darkgray")
abline(v = 0, col = "darkgray")

Zeile 1: Die Objekte mod und cf werden wie im letzten Skript erzeugt.

Zeile 3: Die Funktion fitted.values() erhält wieder als Eingabewert das Modell mod; das Ergebnis des Funktionsaufrufs wird in regr.values abgespeichert.

Zeile 5 bis 7: Die Struktur von regr.values zeigt, dass es sich um einen Vektor der Länge 20 handelt – also genau so viele Komponenten wie die x- beziehungsweise y-Werte besitzen. Als Namen der Komponenten sind die Nummern 1, ..., 20 gesetzt.

Zeile 9: Es werden zwei Plots nebeneinander erzeugt. Links nochmals der Plot aus Abbildung 3 (zum besseren Vergleich der Testdaten und der Regressionswerte). Rechts werden anstelle der Testdaten die Regressionswerte eingetragen (orange). Sie liegen exakt auf der Regressionsgerade.

Zeile 11 bis 22: Nochmals der Quelltext, mit dem Abbildung 3 erzeugt wurde.

Zeile 24 bis 36: Um anstelle der Testdaten die Regressionswerte zu zeichnen, muss innerhalb von plot() gesetzt werden: y = regr.values (Zeile 24). Ansonsten werden lediglich col und lwd verändert (Zeile 25 und 26).

Abbildung 4: Links: Darstellung der Testdaten (rot) und der Regressionsgerade (blau) unter Verwendung der mit coefficients() berechneten Regressionskoeffizienten (wie in Abbildung 3). Rechts: Anstelle der Testdaten werden für die y-Koordinaten die die Regressionswerte verwendet, die mit fitted.values() berechnet und in regr.values abgespeichert wurden (orange).

Die Berechnung des Residuums mit residuals()

Bildet man die Differenz aus den y-Werten und den Regressionswerten W, erhält man das Residuum U:

U = Y - W.

Man kann jetzt:

1. Entweder die Differenz der tatsächlichen y-Koordinaten zu den Regressionswerten als Fehlerbalken in ein Diagramm wie Abbildung 4 eintragen; dies geschieht unten in Abbildung 5 links. (Dabei sind die Messpunkte rot und die Regressionswerte grün gekennzeichnet.)
2. Oder man kann nur die Werte von U gegen die Werte von X aufzutragen – die Regressionsgerade wird also vollständig aus den Daten herausgerechnet; dies geschieht in Abbildung 5 rechts im sogenannten Residualplot.

Da die Testdaten sehr stark korreliert sind, sind einige der Fehlerbalken in Abbildung 5 links kaum zu erkennen. Im Residualplot wird dann die Streuung der Messpunkte um die Regressionsgerade besser sichtbar.

Berechnet wird das Residuum mit Hilfe der Funktion residuals(). Wie man die Funktion aufruft und mit ihrem Ergebnis Abbildung 5 erzeugt, zeigt das folgende Skript:

# mod, cf und regr.values wie oben

res <- residuals(object = mod)

str(res)
# Named num [1:20] 0.8678 0.9926 -1.5755 0.6374 -0.0818 ...
# - attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...

par(mfrow = c(1, 2))

# Links: Regr.gerade mit Fehlerbalken (berechnet ohne Residuum)
plot(x = df.xy$x, y = regr.values,
     col = "green", type = "p", 
     lwd = 2,
     xlim = c(-1, 21), ylim = c(-1, 41),
     xlab = "x", ylab = "y",
     main = paste0( "Fehlerbalken und Regr.gerade ax + b, mit a = ", format(cf[2], digits = 4),
                    ", b = ", format(cf[1], digits = 4) ),
     frame.plot = TRUE)
grid()
abline(a = cf[1], b = cf[2], col = "blue")
# Messpunkte:
points(x = df.xy$x, y = df.xy$y, col = "red")
# Koord.achsen:
abline(h = 0, col = "darkgray")
abline(v = 0, col = "darkgray")
# x0: x, y0: Regressionswert, x1: x-Wert, y1: y-Wert
lapply(X = df.xy$x, FUN = function(i){segments(x0 = df.xy$x[i], y0 = regr.values[i],
                                          x1 = df.xy$x[i], y1 = df.xy$y[i],
                                          col = "red", lwd = 3)})

# Rechts: Residualplot:
plot(x = df.xy$x, y = res,
     col = "red", type = "h",
     lwd = 3,
     xlab = "x", ylab = "y",
     main = "Residuum",
     frame.plot = TRUE)
grid()
#abline(a = cf[1], b = cf[2], col = "blue")
# Koord.achsen:
abline(h = 0, col = "darkgray")
abline(v = 0, col = "darkgray")

Zeile 1: Man benötigt wieder die Objekte mod (das von lm() erzeugte Modell), regr.values (die Regressionswerte) und cf (der Vektor mit den beiden Regressionskoeffizienten); sie werden wie in den Skripten oben erzeugt (hier nicht gezeigt).

Zeile 3: An die Funktion residuals() wird lediglich das Modell mod übergeben; das Ergebnis des Funktionsaufrufs wird in der Variable res abgespeichert.

Zeile 5 bis 7: Die Ausgabe der Struktur von res zeigt, dass es sich um einen Vektor mit 20 Komponenten handelt.

Zeile 12 bis 30: Der Plot links in Abbildung 5 wird noch ohne die Verwendung des Residuums erzeugt. Die roten Punkte sind die Messpunkte, die grünen Punkte die Regressionswerte und die Fehlerbalken werden mit Hilfe der Funktion segments() erzeugt, wozu die Differenz der y-Koordinaten der Messpunkte und Regressionswerte verwendet wird (dies geschieht innerhalb lapply(), um über alle x-Werte zu iterieren, Zeile 28 bis 30).

Zeile 33 bis 43: Für den Residualplot (in Abbildung 5 rechts) werden in plot() für das Argument y die Werte des Residuums, also der Vektor res, eingesetzt. Alle anderen Eigenschaften des Plots sind unabhängig von res.

Abbildung 5: Links: Blau ist die Regressionsgerade eingezeichnet; die Messpunkte (rot) und die Regressionswerte (grün) werden durch Fehlerbalken (rot) verbunden. Die Länge der Fehlerbalken ist zugleich das Residuum. Rechts: Trägt man nur das Residuum gegen die x-Werte auf, erhält man den sogenannten Residualplot. An ihm lässt sich leichter ablesen, wie stark die Messpunkte um die Regressionsgerade streuen.

Zusammenfassung der Auswertung der Testdaten

In Abbildung 6 sind oben nochmals die Definition von Regressionswert und Residuum gezeigt. In der Tabelle darunter sind in der ersten und zweiten Spalte die x- und y-Werte und in der dritten und vierte Spalte Regressionswert W und Residuum U für die in den Abbildungen 1 bis 5 verwendeten Testdaten zu sehen.

Abbildung 6: Tabellarische Darstellung der Testdaten (Zufallsvariablen X und Y), der Regressionswerte W und des Residuums U.

Die Funktion plot.lm()

Sämtliche Plots, die bisher vorgestellt wurden, verwendeten die einfache Funktion plot() aus dem Basis-Paket graphics. Im Aufruf von plot() müssen je ein Vektor für die x- beziehungsweise y-Werte angegeben werden.

Die Funktion plot() ist generisch implementiert und es gibt auch eine Implementierung für die Klasse "lm" ; diese Version befindet sich im Paket stats.

Wird im Argument x ein Objekt der Klasse "lm" an plot() übergeben, so wird der Funktionsaufruf an plot.lm() weitergereicht. Die Plots, die man dann erstellen kann, verwenden Konzepte der Statistik, die hier nicht besprochen wurden. Daher werden hier auch keine entsprechenden Beispiele angegeben. Nähere Informationen kann man aus der R-Dokumentation entnehmen.

Die – zahlreichen – Eingabewerte von plot.lm(), die dem Kenner einen Eindruck vermitteln, welche Diagramme sich erstellen lassen, sind im folgenden Skript gezeigt:

plot(x, which = c(1,2,3,5),
        caption = list("Residuals vs Fitted", "Normal Q-Q",
        "Scale-Location", "Cook's distance",
        "Residuals vs Leverage",
        expression("Cook's dist vs Leverage " * h[ii] / (1 - h[ii]))),
        panel = if(add.smooth) function(x, y, ...)
        panel.smooth(x, y, iter=iter.smooth, ...) else points,
        sub.caption = NULL, main = "",
        ask = prod(par("mfcol")) < length(which) && dev.interactive(),
        ...,
        id.n = 3, labels.id = names(residuals(x)), cex.id = 0.75,
        qqline = TRUE, cook.levels = c(0.5, 1.0),
        add.smooth = getOption("add.smooth"),
        iter.smooth = if(isGlm && binomialLike) 0 else 3,
        label.pos = c(4,2),
        cex.caption = 1, cex.oma.main = 1.25)

Die Funktion summary.lm()

Die Funktion summary() wird meist verwendet, um für Vektoren eine kleine statistische Auswertung bereitzustellen. Das folgende Skript zeigt ein typisches Beispiel: Das Residuum wird mit summary() untersucht (Zeile 6); die Ergebnisse lassen sich leicht mit dem Plot in Abbildung 5 rechts oder der Tabelle in Abbildung 6 nachvollziehen.

mod <- lm(formula = y ~ x, data = df.xy)
# mod und cf wie oben

res <- residuals(object = mod)

summary(object = res)
#     Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# -2.0031 -0.5535  0.1247  0.0000  0.6950  1.6921

Ähnlich wie für plot() gibt es auch für summary() eine Implementierung für die Klasse "lm" im Paket stats, also eine Funktion summary.lm(), an die ein Aufruf von summary() weitergereicht wird, wenn der Eingabewert ein Objekt der Klasse "lm" ist. Aber auch hier gilt, dass die nach der Ausführung von summary.lm() bereitgestellten Informationen weit über die hier besprochenen Konzepte hinausgehen.

summary(object = mod)
# Call:
#   lm(formula = y ~ x, data = df.xy)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -2.0031 -0.5535  0.1247  0.6950  1.6921 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.82706    0.51048    1.62    0.123    
# x            1.96434    0.04261   46.10   <2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 1.099 on 18 degrees of freedom
# Multiple R-squared:  0.9916,  Adjusted R-squared:  0.9911 
# F-statistic:  2125 on 1 and 18 DF,  p-value: < 2.2e-16

Zumindest einige der Größen sind inzwischen bestens bekannt:

• Der Aufruf, mit dem das Modell erzeugt wurde (Zeile 2 und 3).
• Die Auswertung des Residuums (Zeile 5 bis 7).
• Die Regressionskoeffizienten (Zeile 9 bis 12).
• Das Bestimmtheitsmaß, hier als R-squared bezeichnet, das angibt welcher Anteil der Varianz von Y durch das Regressionsmodell erklärt werden kann (Zeile 17).

Das Objekt, das von summary.lm() zurückgegeben wird, ist – wie man im folgenden Skript sehen kann – eine Liste mit 11 Komponenten:

s <- summary(object = mod)
str(s)
# List of 11
# $ call         : language lm(formula = y ~ x, data = df.xy)
# $ terms        :Classes 'terms', 'formula'  language y ~ x
# .. ..- attr(*, "variables")= language list(y, x)
# .. ..- attr(*, "factors")= int [1:2, 1] 0 1
# .. .. ..- attr(*, "dimnames")=List of 2
# .. .. .. ..$ : chr [1:2] "y" "x"
# .. .. .. ..$ : chr "x"
# .. ..- attr(*, "term.labels")= chr "x"
# .. ..- attr(*, "order")= int 1
# .. ..- attr(*, "intercept")= int 1
# .. ..- attr(*, "response")= int 1
# .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
#   .. ..- attr(*, "predvars")= language list(y, x)
# .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
# .. .. ..- attr(*, "names")= chr [1:2] "y" "x"
# $ residuals    : Named num [1:20] 0.8678 0.9926 -1.5755 0.6374 -0.0818 ...
# ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ...
# $ coefficients : num [1:2, 1:4] 0.8271 1.9643 0.5105 0.0426 1.6201 ...
# ..- attr(*, "dimnames")=List of 2
# .. ..$ : chr [1:2] "(Intercept)" "x"
# .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
# $ aliased      : Named logi [1:2] FALSE FALSE
# ..- attr(*, "names")= chr [1:2] "(Intercept)" "x"
# $ sigma        : num 1.1
# $ df           : int [1:3] 2 18 2
# $ r.squared    : num 0.992
# $ adj.r.squared: num 0.991
# $ fstatistic   : Named num [1:3] 2125 1 18
# ..- attr(*, "names")= chr [1:3] "value" "numdf" "dendf"
# $ cov.unscaled : num [1:2, 1:2] 0.2158 -0.0158 -0.0158 0.0015
# ..- attr(*, "dimnames")=List of 2
# .. ..$ : chr [1:2] "(Intercept)" "x"
# .. ..$ : chr [1:2] "(Intercept)" "x"
# - attr(*, "class")= chr "summary.lm"

Um sich einen schnellen Überblick über diese Liste zu verschaffen, ist es wieder einfacher, die Namen der Attribute abzufragen:

s <- summary(object = mod)
names(s)
# [1] "call"          "terms"         "residuals"     "coefficients"  "aliased"       "sigma"         "df"           
# [8] "r.squared"     "adj.r.squared" "fstatistic"    "cov.unscaled"

Ausführliche Informationen zu der Bedeutung der Komponenten dieser Liste s findet man in der Dokumentation zu summary.lm().