Lernziele: Am Ende des Kapitels können Sie

den t-Test für Regressionskoeffizienten im linearen Regressionsmodell in R anwenden

die Ergebnisse interpretieren

Simpson’s Paradox erkennen

1 Hypothesentests in der linearen Regression

Im Zusammenhang mit dem linearen Regressionsmodell kann man verschiedene Hypothesen in Bezug auf die Modellparameter testen. In diesem Kapitel konzentrieren wir uns auf den wichtigsten Hypothesentest

H_0:\beta = 0 \text{ versus } H_1:\beta \neq 0 \tag{1.1}

Man testet also, ob die Steigung der Regressionsgeraden gleich 0 ist. Der obige Hypothesentest bedeutet praktisch, dass die erklärende Variable x keinen Einfluss auf das Ergebnis beziehungsweise die erklärte Variable y hat.

Abhängig von der Situation kann die Alternativhypothese H_1:\beta \neq 0 auch durch die Alternative H_1:\beta >0 ersetzen. Oft ist in diesem Fall jedoch der Test

H_0:\beta \leq 0 \text{ versus } H_1:\beta >0 \tag{1.2}

das in der Praxis relevante Problem. Die einseitige Alternative H_1:\beta >0 ist lediglich dann relevant, wenn wir sicher sind, dass die erklärende Variable x keinesfalls einen negativen Einfluss auf y haben kann. Oft hat man jedoch keine solchen Vorabinformationen, und der Test in Gleichung 1.2 ist sinnvoller.

2 Teststatistik

Hier konzentrieren wir uns auf den Test in Gleichung 1.1. Zum Testen der Hypothese H_0:\beta=0 bietet sich zunächst die geschätzte Steigung \hat{\beta} der Regressionsgeraden an. Um nun beurteilen zu können, ob diese Schätzung \hat{\beta} signifikant vom Wert 0 abweicht, standardisieren wir mit Hilfe der Standardabweichung. Dazu ist folgendes Resultat hilfreich:

Theorem 2.1 (Varianzen der Kleinste-Quadrate-Schätzer) Die Kleinste-Quadrate-Schätzer \hat{\alpha} und \hat{\beta} sind normalverteilt und erwartungstreu mit Varianzen

Var(\hat{\alpha}):=\sigma^2 \frac{\sum_{i=1}^n x_i^2}{n \sum_{i=1}^n (x_i-\bar{x})^2} \tag{2.1}

Var(\hat{\beta}):=\sigma^2 \frac{1}{\sum_{i=1}^n (x_i-\bar{x})^2} \tag{2.2} Der Schätzer

s_{y|x}^2:=\frac{1}{n-2}\sum_{i=1}^n (y_i-\hat{\alpha}-\hat{\beta}x_i)^2 \tag{2.3} ist erwartungstreu und (n-2)s_{y|x}^2 / \sigma^2 hat eine \chi_{n-2}^2-Verteilung.

Theorem 2.1 erlaubt uns nun mit Hilfe von Gleichung 2.2 eine Teststatistik zu konstruieren. Die Statistik Z:=\sqrt{\sum_{i=1}^n (x_i-\bar{x})^2}\frac{\hat{\beta}}{\sqrt{\sigma^2}} hat unter der Hypothese H_0:\beta=0 eine N(0,1) Verteilung. Leider können wir Z jedoch noch nicht als Teststatistik benutzen, da wir \sigma^2 nicht kennen. Dieses Problem können wir jedoch lösen, indem wir einen Schätzer für \sigma^2 einsetzen. Ersetzen wir in der obigen Formel für Z also \sigma^2 durch s_{y|x}^2, so erhalten wir folgende Teststatistik:

Definition 2.1 (T-Test) Wir testen die Hypothese H_0:\beta = 0 mit Hilfe der T-Statistik

T:=\sqrt{\sum_{i=1}^n(x_i-\bar{x})^2}\frac{\hat{\beta}}{\sqrt{s_{y|x}^2}}

Diese Teststatistik hat under der Hypothese H_0:\beta=0 eine t_{n-2}-Verteilung.

Testen wir somit entsprechend Gleichung 1.1 zum Niveau \alpha=0.05, so verwerfen wir H_0 wenn T\leq -t_{n-2;0.025} oder T\geq t_{n-2;0.025}. Testen wir entsprechend Gleichung 1.2, so verwerfen wir nur bei großen Werten von T, und können H_0 zum Niveau \alpha=0.05 verwerfen falls T\geq t_{n-2;0.05}.

3 T-Test im linearen Regressionsmodell

Wir verwenden nun einen Datensatz von Gorman, Williams, und Fraser (2014), der Daten zu drei verschiedenen antarktischen Pinguinarten enthält. Diese wurden 2007-2009 von Dr. Kristen Gorman im Rahmen des Palmer Station Long Term Ecological Research Program gesammelt und aufbereitet.

Zunächst installieren und laden wir dazu das R-Paket palmerpenguins.

install.packages("palmerpenguins")

Anschließend laden wir den Datensatz, wobei wir elf Pinguine entfernen die unvollständige Messungen enthalten. Zusätzlich laden wir die Pakete ggplot2 zur Visualisierung und dplyr zum Filtern des Datensatzes.

library(palmerpenguins)
library(dplyr)
library(ggplot2)
theme_set(theme_minimal())
penguins <- na.omit(penguins)
head(penguins)

str(penguins)

tibble [333 × 8] (S3: tbl_df/tbl/data.frame)
 $ species          : Factor w/ 3 levels "Adelie","Chinstrap",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ island           : Factor w/ 3 levels "Biscoe","Dream",..: 3 3 3 3 3 3 3 3 3 3 ...
 $ bill_length_mm   : num [1:333] 39.1 39.5 40.3 36.7 39.3 38.9 39.2 41.1 38.6 34.6 ...
 $ bill_depth_mm    : num [1:333] 18.7 17.4 18 19.3 20.6 17.8 19.6 17.6 21.2 21.1 ...
 $ flipper_length_mm: int [1:333] 181 186 195 193 190 181 195 182 191 198 ...
 $ body_mass_g      : int [1:333] 3750 3800 3250 3450 3650 3625 4675 3200 3800 4400 ...
 $ sex              : Factor w/ 2 levels "female","male": 2 1 1 1 2 1 2 1 2 2 ...
 $ year             : int [1:333] 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...
 - attr(*, "na.action")= 'omit' Named int [1:11] 4 9 10 11 12 48 179 219 257 269 ...
  ..- attr(*, "names")= chr [1:11] "4" "9" "10" "11" ...

Der Datensatz enthält Messungen zu 344 Pinguinen der Spezies Adelie, Chinstrap und und Gentoo. Die 344 Pinguine wurden auf drei verschiedenen Inseln Biscoe, Dream und Torgersen gesichtet, und die erhobenen Variablen sind

Die Schnabellänge bill_length_mm in Millimetern
Die Schnabeldicke bill_depth_mm in Millimetern
Die Flossenlänge flipper_length_mm in Millimetern
Das Gewicht body_mass_g in Gramm
Das Geschlecht sex
Das Jahr year wann der Pinguin gesichtet und vermessen wurde

Abbildung 3.1: Schnabellänge und -dicke von Pinguinen, Grafik von Allison Horst

Abbildung 3.1 veranschaulicht die gemessenen Abschnitte.¹ Zunächst visualisieren wir Korrelationen zwischen den gemessenen Variablen:

library(dplyr)
my_plot <- ggplot(data = penguins,
                         aes(x = bill_length_mm,
                             y = bill_depth_mm)) +
  geom_point() +
  scale_color_manual(values = c("darkorange","purple","cyan4")) +
  labs(title = "Zusammenhang zwischen Schnabeldicke und -länge",
       x = "Schnabellänge (mm)",
       y = "Schnabeldicke (mm)")
my_plot

Abbildung 3.2: Zusammenhang zwischen Schnabeldicke und -länge

Abbildung 3.2 zeigt die Punktwolke die den Zusammenhang zwischen Schnabellänge und -dicke visualisiert. Wir möchten nun formal testen, ob die Schnabellänge einen Einfluss auf die Schnabeldicke hat. Wir könnten natürlich umgekehrt auch testen, ob die Schnabeldicke einen Einfluss auf die Schnabellänge hat. Prinzipiell könnte man an dieser Stelle auch zunächst eine Korrelationsanalyse zwischen den Variablen Schnabeldicke und Schnabellänge durchführen, indem man z.B. einen Korrelationskoeffizienten berechnet. Da wir jedoch statistische Tests illustrieren wollen, nehmen wir der Einfachheit halber an, dass es biologische Gründe gibt, weswegen die Schnabellänge die erklärende Variable sein sollte.

linReg <- lm(bill_depth_mm ~ bill_length_mm, data = penguins)
summary(linReg)


Call:
lm(formula = bill_depth_mm ~ bill_length_mm, data = penguins)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.1548 -1.4291  0.0122  1.3994  4.5004 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    20.78665    0.85417  24.335  < 2e-16 ***
bill_length_mm -0.08233    0.01927  -4.273 2.53e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.92 on 331 degrees of freedom
Multiple R-squared:  0.05227,   Adjusted R-squared:  0.04941 
F-statistic: 18.26 on 1 and 331 DF,  p-value: 2.528e-05

Wir sehen, dass die erklärende Variable Schnabellänge einen signifikanten Einfluss auf die Schnabeldicke hat, wenn wir das Testlevel \alpha=0.05 zu Grunde legen. Wir sehen außerdem, dass \hat{\beta}=-0.08233 ist, also die Schnabellänge einen negativen Einfluss auf die Schnabeldicke zu haben scheint.

Aufgabe 3.1 (Berechung des Schätzers \hat{\beta}) Die R-Funktion lm nimmt uns wie oben gezeigt die Arbeit ab, den Test in Gleichung 1.1 durchzuführen. Wir möchten nun jedoch den Test händisch durchführen und überprüfen, ob das Ergebnis mit dem obigen übereinstimmt. Berechnen Sie dazu zunächst \hat{\beta}:=\frac{\sum_{i=1}^n (x_i-\bar{x})(y_i-\bar{y})}{\sum_{i=1}^n (x_i-\bar{x})^2}.

Tipp zu Aufgabe 3.1 anzeigen

Um \hat{\beta} zu berechnen, können Sie die in Aufgabe 1 angegebene Formel verwenden. Legen Sie sich dazu zunächst Variablen x und y an, etwa so:

x <- penguins$bill_length_mm
y <- penguins$bill_depth_mm

Anschließend verfahren sie analog für \bar{x} und \bar{y} und benutzen dann die Funktion sum um \hat{\beta} zu berechnen.

Aufgabe 3.2 (Berechnung der Teststatistik T) Wir haben in Aufgabe 3.1 nun \hat{\beta} händisch berechnet. Nun wollen wir Gleichung 2.2 verwenden, um den Wert der Teststatistik T zu berechnen. Berechnen Sie dazu zunächst s_{y|x}^2 mit Hilfe der Funktion var und berechnen dann mit Hilfe von Gleichung 2.2 den Wert von T.

Beachten Sie, dass die Schätzer \alpha und \beta gegeben sind durch:

\hat{\alpha} = \bar{y}-\hat{\beta}\bar{x}

\hat{\beta}:=\frac{\sum_{i=1}^n (x_i−\bar{x})(y_i−\bar{y})}{\sum_{i=1}^n (x_i-\bar{x})^2}

wobei \bar{x}:=\sum_{i=1}^n x_i und \bar{y}:=\sum_{i=1}^n y_i.

Tipp zu Aufgabe 3.2 anzeigen

x <- penguins$bill_length_mm
y <- penguins$bill_depth_mm
xbar <- mean(penguins$bill_length_mm)
ybar <- mean(penguins$bill_depth_mm)
hatBeta <- sum((x - xbar) * (y - ybar)) / (sum((x - xbar)^2))
hatBeta

[1] -0.08232675

Wir können die Schätzung s_{y|x}^2 für die Varianz \sigma^2 wie folgt berechnen:

hatAlpha <- ybar -hatBeta*xbar
n <- nrow(penguins)
s_yx_sq <- (1 / (n - 2)) * sum((y - hatAlpha - hatBeta * x)^2)
s_yx_sq

[1] 3.686296

Setzen Sie dieses Ergebnis nun zusammen mit dem Ergebnis aus Aufgabe 1 in die Formel für T ein.

Aufgabe 3.3 (Berechnung des p-Werts) Wir haben in Aufgabe 1 und 2 nun T händisch berechnet. Berechnen Sie nun den p-Wert mit Hilfe der Funktion pt.

4 Simpson’s Paradox

Wir möchten nun dem Zusammenhang von oben näher auf den Grund gehen. Es erscheint paradox, dass Pinguine mit einer längeren Schnabellänge einen weniger dicken Schnabel haben. Man würde erwarten, dass Pinguine mit größeren Schnäbeln generell über größere anatomische Proportionen verfügen und daher auch dickere Schnäbel haben. Der Hypothesentest hat jedoch gezeigt, dass ein umgekehrter Zusammenhang zu gelten scheint. Abbildung 4.1 visualisiert die Ergebnisse der linearen Regression.

my_plot + geom_smooth(method = "lm", se = FALSE)