Nur ein wenig lineare Regression

Der tipping Datensatz wird oft analysiert. Das Verhältnis von Trinkgeld (tip) und Rechnungsbetrag (total_bill) steht dabei im Vordergrund einer lineare Regressionsanalyse. So auch hier. Wir wollen die einzelnen Angaben von R dabei in den Fokus rücken und einmal Hinterfragen, was wir bei der Ausgabe von R eigentlich genau sehen, woher es kommt und wie man es interpretieren kann.

Zunächst laden wir dazu die tipping Daten mittels

library(mosaic)
download.file("https://goo.gl/whKjnl", destfile = "tips.csv")
tips <- read.csv2("tips.csv")

in den Arbeitsspeicher.

Eine lineares Modell wird schnell mit

linMod <- lm(tip ~ total_bill, data = tips)

erstellt. Betrachten wir die Zusammenfassung:

summary(linMod)

## 
## Call:
## lm(formula = tip ~ total_bill, data = tips)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1982 -0.5652 -0.0974  0.4863  3.7434 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.920270   0.159735   5.761 2.53e-08 ***
## total_bill  0.105025   0.007365  14.260  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.022 on 242 degrees of freedom
## Multiple R-squared:  0.4566, Adjusted R-squared:  0.4544 
## F-statistic: 203.4 on 1 and 242 DF,  p-value: < 2.2e-16

Die zentrale Frage bei einer linearen Regression ist, finden wir einen linearen Zusammenhang in unserer Stichprobe, den wir auf die Population (als die Grundgesamtheit) übertragen können.

Die Spalte Estimate im Abschnitt Coefficients liefert uns in unser Stichprobe einen möglichen linearen Zusammenhang gemäß

\[\hat{y}_{\text{tip}} = \hat{\beta}_{\text{0}} + \hat{\beta}_{\text{total_bill}} \cdot x_{\text{total_bill}},\]

mit den Regressionskoeffizienten \(\hat{\beta}_0=0.9202696\) und \(\hat{\beta}_{\text{total_bill}}=0.1050245\).

Graphisch ergibt sich damit das Modell wie folgt:

# Statt plotModel(linMod) besser:
mypanel <- function(x, y) {
    # Scatterplot:
    panel.xyplot(x, y, col = "darkgreen") 
    # Regressionsgerade:
    panel.abline(linMod, col = "red", lwd = 1.2, lty = 2)
}
xyplot(
    tip ~ total_bill, data = tips, 
    panel = mypanel,
    main  = "Streudiagramm der Trinkgelder",
    ylab  = "Trinkgeld",
    xlab  = "Rechnungsbetrag",
    key = list(
            space = "bottom", padding.text = 8,
            lines = list(col = c("red"), lty = c(2), lwd = 1.2),
            text = list(c("Regressionsgerade"))
          )
     )

Was hat es mit dem y-Achsenabschnitt \(\hat{\beta}_0\) auf sich?

Ist es etwa eine Art Grundtrinkgeld, mit dem der Kellern rechnen kann, auch wenn der Kunde gar nichts bestellt?

Nun ja, es so etwas in der Art, aber eben ein rein fiktiver Wert, der durch die Konstruktion der Parameter entsteht. Eine (affin-)lineare Gerade geht nun einmal irgendwann durch die y-Achse (wenn sie nicht parallel dazu ist) und es kann passieren, dass eine sinnvolle Interpretation nicht so ohne weiteres möglich ist.

Wir können aber dieses Grundtrinkgeld heraus nehmen und den y-Achsenabschnitt auf Null setzen. Dazu ziehen wir \(\hat{\beta}_0\) einfach von alle Trinkgeldern ab. Wir erhalten quasi nur noch den Trinkgeldzuwach.

beta_0 <- coef(linMod)["(Intercept)"]  # Grundtrinkgeld
tips$delta_tip <- tips$tip - beta_0    # wird abgezogen

Vergleichen wir das alte lineare Modell mit dem neuen Modell (linModDelta):

linModDelta <- lm(delta_tip ~ total_bill, data = tips)
summary(linModDelta)

## 
## Call:
## lm(formula = delta_tip ~ total_bill, data = tips)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1982 -0.5652 -0.0974  0.4863  3.7434 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.549e-15  1.597e-01    0.00        1    
## total_bill   1.050e-01  7.365e-03   14.26   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.022 on 242 degrees of freedom
## Multiple R-squared:  0.4566, Adjusted R-squared:  0.4544 
## F-statistic: 203.4 on 1 and 242 DF,  p-value: < 2.2e-16

In diesem Modell ist der Wert für den y-Achsenabschnitt numerisch gleich 0. – Ja, da mag zwar \(-4.5487837\times 10^{-15}\) stehen, jedoch sind so kleine Werte der jedem Rechner inne wohnenden Ungenauigkeit in der Gleitkomma-Arithmetik geschuldet und ist faktisch gleich 0.

Der Wert für die Steigung lautet weiterhin \(0.1050245\). Das war auch zu erwarten, denn wir haben unsere Regressionsgerade eigentlich nur um \(\hat{\beta}_0\) nach unten verschoben. (Der Fachmann spricht von einer Translation (Parallelverschiebung)¹ um \(-\hat{\beta}_0\).

# Statt plotModel(linModDelta) besser:
mypanel <- function(x, y) {
    # Scatterplot:
    panel.xyplot(x, y, col = "darkgreen") 
    # Regressionsgerade:
    panel.abline(linModDelta, col = "red", lwd = 1.2, lty = 2)
}
xyplot(
    delta_tip ~ total_bill, data=tips, 
    panel = mypanel,
    main  = "Streudiagramm der Delta Trinkgelder",
    ylab  = "Delta Trinkgeld",
    xlab  = "Rechnungsbetrag",
    key = list(
            space="bottom", padding.text=8,
            lines=list(col=c("red"), lty=c(2), lwd=1.2),
            text=list(c("Regressionsgerade")))
     )

Vergleichen wir die beiden Zusammenfassungen, so stellen wir fest das sich mit Ausnahme der [Intercept] Zeile praktisch nichts geändert hat. Das ist kein Wunder, sondern Absicht!

Die Regressionsgerade stellt für unsere Stichprobe die Gerade mit dem geringsten Fehler an den Datenpunkten dar. Mathematisch heißt das folgendes:

An den \(n=244\) Datenpunkten unserer Stichprobe \((x_i, y_i)=(tips\$total\_bill[i], tips\$tip[i])\) [für \((i=1, \dots, n)\)] sind die Residuen, also die Fehlerterme,

\[ \hat{e}_i =\hat{y}_i - y_i = \left[\hat{\beta}_{\text{0}} + \hat{\beta}_{\text{total_bill}} \cdot x_i\right] - y_i \]

durch die verwendete Methode der kleinsten Quadrate² quadratisch minimal. Kurz:

\[ \sum_{i=1}^n (\hat{e}_i)^2 \text{ ist minimal!} \]

Wir können diese Fehlerterme graphisch ansehen um die Varianz der Residuen zu sehen. Dazu ziehen wir von allen Datenpunkten \(y_i\) den geschätzten Wert \(\hat{y}_i\) ab und erstellen ein neues lineares Modell:

beta_total_bill <- coef(linModDelta)["total_bill"]
tips$error_tip <- (tips$tip - beta_0 - beta_total_bill * tips$total_bill)
linModError <- lm(error_tip ~ total_bill, data = tips)
summary(linModError)

## 
## Call:
## lm(formula = error_tip ~ total_bill, data = tips)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1982 -0.5652 -0.0974  0.4863  3.7434 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.900e-15  1.597e-01       0        1
## total_bill  -8.740e-17  7.365e-03       0        1
## 
## Residual standard error: 1.022 on 242 degrees of freedom
## Multiple R-squared:  6.665e-31,  Adjusted R-squared:  -0.004132 
## F-statistic: 1.613e-28 on 1 and 242 DF,  p-value: 1

Also Diagramm sieht es dann so aus:

# Statt plotModel(linModError) besser:
mypanel <- function(x, y) {
    # Scatterplot:
    panel.xyplot(x, y, col = "darkgreen") 
    # Regressionsgerade:
    panel.abline(linModError, col = "red", lwd = 1.2, lty = 2)
}
xyplot(
    error_tip ~ total_bill, data = tips, 
    panel = mypanel,
    main  = "Streudiagramm der Residuen",
    ylab  = "Residuen",
    xlab  = "Rechnungsbetrag",
    key = list(
            space = "bottom", rows = 3, padding.text = 8,
            lines = list(col=c("red"), lty = c(2), lwd = 1.2),
            text = list(c("Regressionsgerade / x-Achse"))
          )
)

Wir können die Graphik im wesentlichen auch einfacher über den Befehl

xyplot(residuals(linMod) ~ fitted(linMod))

erhalten.

Betrachten wir kurz nur die Residuen:

favstats(~residuals(linMod))

##        min         Q1      median        Q3      max          mean       sd   n
##  -3.198225 -0.5651615 -0.09744499 0.4863111 3.743435 -2.022281e-17 1.019943 244
##  missing
##        0

Wir sehe, dass wir in der Zusammenfassung immer genau diese Werte unter dem Abschnitt Residuals gefunden haben. Minimum, das 1. Quantil, der Median, das 3. Quantil und das Maximum stimmen überein.

Der erwartungstreue und unverzerrte Schätzer für den Standardfehler der Residuen, lautet

\[ \begin{align*} SE_{\text{Residuen}} &= \sqrt{\frac{1}{n-2} \cdot \sum_{i=1}^n (\hat{e_i})^2} = \sqrt{\frac{n-1}{n-2} \cdot \frac{1}{n-1} \cdot \sum_{i=1}^n (\hat{e_i})^2} \\ &= \sqrt{\frac{n-1}{n-2}} \cdot \sqrt{\frac{1}{n-1} \cdot \sum_{i=1}^n (\hat{e_i})^2} \\ &= \sqrt{\frac{n-1}{n-2}} \cdot s_{\text{Residuen}} \end{align*} \]

Also finden wir den Wert Residual standard error aus der Zeile

## Residual standard error: 1.022 on 242 degrees of freedom

in dem wir den in den favstats gefundenen Wert für die Standardabweichung entsprechen korrigieren:

\[ SE_{\text{Residuen}} = \sqrt{\frac{n-1}{n-2}} \cdot s_{\text{Residuen}} = \sqrt{\frac{243}{242}} \cdot 1.0199426 = 1.0220477 \]

Der Median der Residuen ist nicht gleich Null, wie der Mittelwert. (Welcher auch hier als numerisch Null interpretiert werden muss!) Es könnte also eine linkssteile, rechtsschiefe Verteilung der Residuen vorliegen. Betrachten wir dazu das Histogramm:

histogram(~residuals(linMod), nint = 19)

Schon beim ersten Blick auf das Histogramm kann an eine Normalverteilung der Residuen nicht mehr so ganz geglaubt werden.

Ein Shapiro-Wilk-Test³ hat als Nullhypothese die Annahme, dass die Daten normalverteilt sind!

shapiro.test(residuals(linMod))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(linMod)
## W = 0.96728, p-value = 2.171e-05

Davon ist nach dem Ergebnis eben sowenig auszugehen, wie nach einem Blick auf das QQ-Normal-Diagramm:

qqnorm(residuals(linMod), col = "darkgreen")

Ein K.O.-Kriterium für gute Prognosen.

Wie gut aber beschreibt unsere Regressionsgerade die Daten?

Als Maß dafür können wir das Bestimmtheitsmaß nehmen.

Ein kurzer Blick auf die Situation, der Mittelwert der Trinkgelder ist

\[ \bar{y} = \frac{1}{n} \cdot \sum_{i=1}^n y_i = 2.9982787. \]

Wir erhalten so folgendes Diagramm:

mypanel <- function(x, y) {
    panel.xyplot(x, y)
    panel.abline(h = mean(y), lwd = 1.2, lty = 2, col = "darkgreen")
    panel.lmline(x, y, col = "red", lwd = 1.2, lty = 2)
}
xyplot(
    tip ~ total_bill, data = tips, 
    panel = mypanel,
    main  = "Streudiagramm der Trinkgelder",
    ylab  = "Trinkgeld",
    xlab  = "Rechnungsbetrag",
    key = list(
            space = "bottom",
            padding.text = 8,
            columns = 2,
            just = c("center", "bottom"),
            lines = list(col = c("darkgreen", "red"), lty = c(2, 2), lwd = 1.2),
            text = list(c(expression(bar(y)), expression(hat(beta)[0]+hat(beta)[total_bill] * x[total_bill]))),
            text = list(c("Mittelwert Trinkgeld", "Regressionsgerade"))
    )
)

Die Varianz \(s^2_{y_i}=1.9144546\) beschreibt die mittlere quadratische Abweichung der Datenpunkte \(y_i\) vom Mittelwert \(\bar{y}\). Diese Varianz lässt sich Zerlegen in einen Anteil, der durch die Regressionsgerade erklärt wird und in einen Anteil, der durch die Regressionsgerade nicht erklärt wird.

\[ s^2_{y_i} = s^2_{\hat{y}_i} + s^2_{\hat{e}_i} \] Dividiert man beider Seiten durch die Varianz \(s^2_{y_i}\), so normiert man den Ausdruck und kann den Faktor \(\frac{1}{n-1}\) (bzw. \(\frac{1}{n}\)) herauskürzen. Es bleibt dann:

\[ 1 = \frac{\sum_{i=1}^n (\bar{y}- \hat{y_i})^2}{\sum_{i=1}^n (\bar{y}-y_i)^2} + \frac{\sum_{i=1}^n (\hat{e_i})^2}{\sum_{i=1}^n (\bar{y}-y_i)^2} \]

Multipliziert man beide Seiten mit \(\sum_{i=1}^n (y_i)^2\), so erhält man:

\[ \sum_{i=1}^n (\bar{y}- y_i)^2 = \sum_{i=1}^n (\bar{y}- \hat{y_i})^2+ \sum_{i=1}^n (\hat{e_i})^2 \]

Zur Vereinfachung nennt man die einzelnen Summen in dem Ausdruck wie folgt:

• Der erste Ausdruck heißt Gesamtvarianz oder total sum of squares oder kurz \(SS_T\), (oder TSS) er ist die Summe der quadrierten Differenzen

\[ SS_T = \sum_{i=1}^n (\bar{y}-y_i)^2 \]

• Der zweite Ausdruck heißt Modellvarianz oder model sum of squares oder kurz \(SS_M\) (oder RSS), er ist die Summe der quadrierten Differenzen aus dem Mittelwert \(\bar{y}\) und der Punkte auf der Regressionsgeraden \(\hat{y}_i\):

\[ SS_M = \sum_{i=1}^n (\bar{y}-\hat{y}_i)^2 \]

• Der dritte Ausdruck heißt Gesamt-Verhersage-Fehler, Fehlersteuung der Regression oder error sum of squares oder kurz \(SS_E\) (oder ESS), er ist die Summe der quadratischen Differenz aus den Datenpunkten \(y_i\) und den Punkten der Regressionsgeraden \(\hat{y}_i\):

\[ SS_E = \sum_{i=1}^n (\hat{y}_i-y_i)^2 = \sum_{i=1}^n (\hat{e}_i)^2 \]

Wir können daher auch kurz

\[ SS_T = SS_M + SS_E \] schreiben und sparen uns die ganzen Summenzeichen.

Die Güte einer Regression wollen wir durch den Anteil der durch das Model erklärten Varianz (also der \(SS_M\)) ausdrücken und stellen daher nach \(SS_M\) um:

\[ SS_M = SS_T - SS_E \] Teilen wir beide Seiten durch \(SS_T\) also der maximalen (weil totalen) Quadratsumme, so erhalten wir: \[ \frac{SS_M}{SS_T} = \frac{SS_T}{SS_T} - \frac{SS_E}{SS_T} = 1 - \frac{SS_E}{SS_T} \]

Den Ausdruck \(\frac{SS_M}{SS_T}\) nennen wir Bestimmtheitsmaß und schreiben dafür \(R^2\). Es ist ein Wert zwischen 0 und 1, der den Anteil der durch das Modell beschriebenen Varianz in Bezug auf die Gesamtvarianz angibt. Kraft Definition ist \(R^2\) im eindimensionalen Fall tatsächlich das Quadrat des (Pearson-)Korrelationskoeffizienten \(r\). (M.a.W.: \(R^2= r^2\).)

In unserer Zusammenfassung des linearen Models findet sich dieser Wert auch. Und zwar unter dem Begriff:

## Multiple R-squared:  0.4566, 

Es gilt ja:

\[ R^2 = 1 - \frac{SS_E}{SS_T} = 1 - \frac{s^2_{\hat{e}_i}}{s^2_{y_i}} = 1 - \frac{1.0402829}{1.9144546} = 0.4566166 \]

Der Wert

## ..., Adjusted R-squared:  0.4544

erklärt sich daraus⁴, dass das Bestimmheitsmaß um so größer wird je größer die Zahl der unabhängigen Variablen wird. Und zwar unabhöngig davon, ob weitere unabhängige Variablen wirklich einen Beitrag zur Erklärungskraft liefern. Daher nutzt man besser das korrigierte Bestimmtheitsmaß (engl.: adjusted R-squared) \(\bar{R}^2\):

\[ \begin{align*} \bar{R}^2 &= 1- (1-R^2) \cdot \frac{n-1}{n-p-1}\\ &= R^2 - (1-R^2) \cdot \frac{p}{n-p-1} \end{align*} \]

Wobei \(p\) die Anzahl der unabhängigen Variablen im Modell darstellt. In unserem Beispiel gilt daher:

\[ \begin{align*} \bar{R}^2 &= 1 - (1-R^2) \cdot \frac{n-1}{n-p-1} \\ &= 1 - (1- 0.4566166) \cdot \frac{244-1}{244- 1- 1} \\ &= 0.4543712 \end{align*} \]

Vorsicht: Das korrigierte Bestimmtheitsmaß ist nicht mehr an das Intervall \([0; 1]\) gebunden! Es kann negative Werte annehmen, ist in der Regel kleiner als das (unkorrigierte) Bestimmtheitsmaß und erreicht die obere Grenze (\(\bar{R}^2=1\)) genau dann, wenn \(R^2 = 1\) ist.

Bei der Gesamtsignifikanz des Modells (auch Overall-F-Test genannt) wird geprüft, ob mindestens eine Variable einen Erklärungsgehalt für das Modell liefert.

Falls diese Hypothese verworfen wird ist somit das Modell nutzlos. Dieser Test lässt sich so interpretieren als würde man die gesamte Güte des Modells, also das \(R^2\) des Modells, testen. Aus diesem Grund wird der F-Test der Gesamtsignifikanz des Modells auch als Anpassungsgüte-Test bezeichnet. Die Nullhypothese des F-Test der Gesamtsignifikanz des Modells sagt aus, dass alle erklärenden Variablen keinen Einfluss auf die abhängige Variable haben. Sowohl die abhängige Variable als auch die unabhängigen Variablen können binär (kategoriell) oder metrisch sein. Der Wald-Test kann dann die Hypothesen testen (ohne Einbezug des Achsenabschnittes):

\[ H_{0}\colon \beta _{1}=\beta _{2}=\ldots =\beta _{k}\;=\;0\Rightarrow R^{2}=0 \] gegen

\[ H_{1}:\beta _{j}\;\neq \;0\;\mathrm {f{\ddot {u}}r\;mindestens\;ein} \;j\in \{1,\ldots ,k\}\Rightarrow R^{2}\neq 0 \]

Die Teststatistik dieses Tests lautet

\[ \begin{aligned} F\;\;{\stackrel {H_{0}}{=}}{\frac {R^{2}}{1-R^{2}}} \cdot {\frac {n-p-1}{p}}\;\;{\stackrel {H_{0}}{\sim }}\;\;F(p,n-p) \end{aligned}. \]

mit \(p\) und \(\displaystyle (n-p-1)\) Freiheitsgraden. Überschreitet der empirische F-Wert einen kritischen F-Wert, der zu einem a priori festgelegten Signifikanzniveau \(\alpha\), so verwirft man die Nullhypothese \(H_{0}\). Das \(R^{2}\) ist dann ausreichend groß und mindestens ein Regressor trägt also vermutlich genügend viel Information zur Erklärung von \(y\) bei. Es ist naheliegend bei hohen F-Werten die Nullhypothese zu verwerfen, da ein hohes Bestimmtheitsmaß zu einem hohen F-Wert führt. Wenn der Wald-Test für eine oder mehrere unabhängige Variablen die Nullhypothese ablehnt, dann kann man davon ausgehen, dass die zugehörigen Parameter ungleich Null sind, so dass die Variable(n) in das Modell mit einbezogen werden sollten.

In unserem Beispiel ist

\[ F={\frac {R^{2}}{1-R^{2}}} \cdot {\frac {n-p-1}{p}} = \frac{0.4566166}{1-0.4566166} \cdot \frac{244-1-1}{1} = 203.3577233 \]

der Wert in der Zeile

## F-statistic: 203.4 on 1 and 242 DF,  p-value: < 2.2e-16

mit Parametern \(p=1\) und \(n-p-1=242\) Freiheitsgraden.

Der p-Wert von (numerisch) 0, liefert also ein hinreichendes Indiz dafür, dass der Rechnungsbetrag einen echten Beitrag liefert.