install.packages("palmerpenguins")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.
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_mmin MillimeternDie Schnabeldicke
bill_depth_mmin MillimeternDie Flossenlänge
flipper_length_mmin MillimeternDas Gewicht
body_mass_gin GrammDas Geschlecht
sexDas Jahr
yearwann der Pinguin gesichtet und vermessen wurde
Abbildung 3.1 veranschaulicht die gemessenen Abschnitte.1 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 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-05Wir 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}.
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.
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)
Wir verfeinern Abbildung 3.2 daher nun wie folgt:
my_plot2 <- ggplot(data = penguins,
aes(x = bill_length_mm,
y = bill_depth_mm,
group = species)) +
geom_point(aes(color = species,
shape = species),
size = 3,
alpha = 0.8) +
scale_color_manual(values = c("darkorange","purple","cyan4")) +
labs(title = "Zusammenhang zwischen Schnabeldicke und -länge nach Spezies",
subtitle = "Spezies: Adelie, Chinstrap und Gentoo Pinguine",
x = "Schnabellänge (mm)",
y = "Schnabeldicke (mm)",
color = "Penguin species",
shape = "Penguin species") +
theme(legend.position = c(0.85, 0.15),
plot.title.position = "plot",
plot.caption = element_text(hjust = 0, face= "italic"),
plot.caption.position = "plot")
my_plot2
Wir sehen nun, dass innerhalb jeder Spezies (Adelie, Chinstrap oder Gentoo) ein positiver Zusammenhang zwischen der Schnabellänge und -dicke vorzuliegen scheint. Die farbigen Punktwolken zeigen, dass mit zunehmender Schnabellänge in Millimetern auch tendenziell die Schnabeldicke zunimmt.
Aufgabe 4.1 (Lineare Regression für die Adelie-Pinguine) Im letzten Abschnitt haben Sie gelernt, wie man mit Hilfe der Funktion lm lineare Regressionsmodelle berechnen und Hypothesentests durchführen kann. Führen Sie den Test in Gleichung 1.1 für die Spezies der Adelie-Pinguine durch, welche der dunkelorangenen Punktwolke in Abbildung 4.2 entsprechen.
Aufgabe 4.2 (Lineare Regression für die Chinstrap-Pinguine) Im letzten Abschnitt haben Sie gelernt, wie man mit Hilfe der Funktion lm lineare Regressionsmodelle berechnen und Hypothesentests durchführen kann. Führen Sie den Test in Gleichung 1.1 für die Spezies der Chinstrap-Pinguine durch, welche der lila Punktwolke in Abbildung 4.2 entsprechen. Verfahren Sie dann analog für die Gentoo-Pinguine, welche der cyanen Punktwolke in Abbildung 4.2 entsprechen.
Falls Sie einen Tipp benötigen, schauen Sie sich die Hilfestellung zu Arbeitsauftrag 4 noch einmal an.
Die Aufgaben 4 und 5 haben verdeutlicht, dass in jeder Spezies ein positiver Zusammenhang zwishen der Schnabeldicke und -länge existiert. Wir können dies wie folgt in Abbildung 4.2 einbauen:
my_plot2 + geom_smooth(method = "lm", se = FALSE, aes(color = species))
Anmerkung (Zum Ausprobieren - Konfidenzintervalle für \hat{\beta}). Ändern Sie im obigen Code einmal se=FALSE auf se=TRUE. Sie erhalten dann zusätzlich eine Visualisierung der Unsicherheit über den Parameter \beta, indem das Konfidenzintervall für \beta verwendet wird um die Bandbreite an linearen Regressionsgeraden darzustellen, die aus dem 95% Konfidenzintervall für \beta resultieren.
5 Erklärung
Das obige Phänomen heißt auch Simpson’s Paradox und wurde in der Literatur vielfach diskutiert (Pearl 2009). In einfachen Worten beschreibt es die Tatsache, dass sich die Richtung eines (linearen) Zusammenhangs umkehren kann, wenn man klassierte Daten ohne die Klassenzugehörigkeiten analysiert (Pearl, Glymour, und Jewell 2016). In unserem Beispiel war der Zusammenhang zwischen Schnabeldicke und -länge ohne die Information der Pinguin-Spezies zunächst signifikant und \hat{\beta} negativ. Unter Berücksichtigung der Klassen, hier der Pinguin-Spezies, wurde das Paradox aufgelöst und in jeder Spezies konnte ein positiver Zusammenhang nachgewiesen werden.
6 Zusammenfassung
In diesem Kapitel haben Sie gelernt wie man mit Hilfe von R lineare Regressionsmodelle analysieren und Hypothesentests für den Regressionskoeffizienten \beta durchführen kann. Zudem haben Sie Simpson’s Paradox kennengelernt, welches häufig bei linearen Regressionsmodellen auftritt, hier anhand des realen Datensatzes der Palmer Pinguine von Gorman, Williams, und Fraser (2014).
Lösungen der Aufgaben
Lizenz
Diese Lerneinheit “Statistische Hypothesentests” wurde von Riko Kelter, Alexander Schurr und Susanne Spies unter Mithilfe von Annika Hirth an der Universität Siegen entwickelt. Es ist lizenziert unter der CC-BY-SA 4.0 Lizenz und ist verfügbar auf ORCA.nrw.
Literatur
Gorman, Kristen B., Tony D. Williams, und William R. Fraser. 2014. „Ecological Sexual Dimorphism and Environmental Variability within a Community of Antarctic Penguins (Genus Pygoscelis)“. PLOS ONE 9 (3): e90081. https://doi.org/10.1371/JOURNAL.PONE.0090081.
Pearl, Judea. 2009. Causality: Models, reasoning, and inference, second edition. New York: Cambridge University Press. https://doi.org/10.1017/CBO9780511803161.
Pearl, Judea, Madelyn Glymour, und Nicholas P. Jewell. 2016. Causal Inference in Statistics: A Primer. Chichester, UK: Wiley.
Fußnoten
Grafik von Allison Horst, frei verfügbar unter https://allisonhorst.github.io/palmerpenguins/articles/art.html↩︎