Zum Datensatz

Für diese Übung brauchst du den Datensatz blood-pressure.csv. Den Datensatz findest du auf Moodle. Der Datensatz enthält folgende Angaben:

  • gender: Geschlecht (female oder male)
  • age: Alter in Jahren
  • height: Körpergrösse in cm
  • weight: Gewicht in kg
  • drug: Medikamenteneinnahme (no oder yes)
  • bp1 und bp2: Blutdruck zu zwei verschiedenen Messzeitpunkten (Prä/Post), in mmHg

Lade den Datensatz herunter und lies ihn in R ein. Prüfe, ob der Datensatz richtig eingelesen wurde (z. B. durch str() oder visuelle Überprüfung).

Korrelation

Aufgaben

  1. Erstelle ein Streudiagramm für die Variablen bp1 und bp2. Ist ein linearer Zusammenhang plausibel? Welcher Korrelationskoeffizient ist in diesem Fall adäquat?
  2. Wie stark ist der Zusammenhang zwischen bp1 und bp2? Tipp: cor(), use = "pairwise.complete.obs"
  3. Wie lautet das 95%-Konfidenzintervall für den Korrelationskoeffizienten nach Pearson? Tipp: cor.test()

Lösungen

Zuerst erstellen wir wieder ein Streudiagramm für die Variablen bp1 und bp2.

library(rio)
df.bp <- import("../data/blood-pressure.csv")
plot(df.bp$bp1, df.bp$bp2)

Ein linearer Zusammenhang scheint plausibel zu sein, weshalb der Korrelationskoeffizient nach Pearson berechnet werden kann.

Die Korrelation können wir mit cor() berechnen:

cor(df.bp$bp1, df.bp$bp2, use = "pairwise.complete.obs")
## [1] 0.9076092
# Wegen fehlender Werte werden nur Beobachtungen mit paarweise vollständigen Daten berücksichtigt.

Der Zusammenhang ist positiv und mit 0.908 sehr stark. Dies ist nicht überraschend, da ja zweimal der Blutdruck an derselben Person gemessen wird. Beachte, dass Personen mit und ohne Medikament nicht getrennt analysiert werden.

Mit cor.test() erhalten wir zudem ein Konfidenzintervall und einen p-Wert (\(H_0: \rho = 0\)):

cor.test(df.bp$bp1, df.bp$bp2, use = "pairwise.complete.obs")
## 
##  Pearson's product-moment correlation
## 
## data:  df.bp$bp1 and df.bp$bp2
## t = 21.072, df = 95, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8646884 0.9373728
## sample estimates:
##       cor 
## 0.9076092
# Wegen fehlender Werte werden nur Beobachtungen mit paarweise vollständigen Daten berücksichtigt.

Bestimmtheitsmass

Aufgabe

  1. Schätze noch einmal das lineare Modell bp2 ~ bp1.

  2. Welchen Wert hat das Bestimmtheitsmass? Wie ist das Bestimmtheitsmass zu interpretieren?

  3. Welcher Zusammenhang besteht zwischen dem Bestimmtheitsmass (R-squared) und dem Korrelationskoeffizienten nach Pearson?

Lösung

lm.bp <- lm(bp2 ~ bp1, df.bp)
summary(lm.bp)
## 
## Call:
## lm(formula = bp2 ~ bp1, data = df.bp)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.0245  -3.3893  -0.0881   3.9283  16.8956 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -10.04294    6.17738  -1.626    0.107    
## bp1           1.06352    0.05047  21.072   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.419 on 95 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.8238, Adjusted R-squared:  0.8219 
## F-statistic:   444 on 1 and 95 DF,  p-value: < 2.2e-16

Das Bestimmtheitsmass ist im Output mit R-squared angegeben und beträgt 0.8238. Bei einer einfachen linearen Regression ist das Bestimmtheitsmass das Quadrat des Korrelationskoeffizienten nach Pearson.

Das Bestimmtheitsmass sagt uns, wie gut das Modell die abhängige Variable erklären kann. In diesem Fall können mit dem Prädiktor bp1 82.38% der Variabilität der abhängigen Variablen bp2 erklärt werden.

Modelldiagnostik

Lineare Regressionsmodelle setzen bestimmte Bedingungen voraus, damit die Resultate korrekt interpretiert werden können. Diese Voraussetzungen gelten insbesondere für die korrekte Schätzung der Standardfehler und somit für die Berechnung von p-Werten und Konfidenzintervallen für Regressionskoeffizienten.

Aufgaben

  1. Was sind die vier wichtigsten Annahmen für die Gültigkeit einer linearen Regression?

  2. Prüfe diese Voraussetzungen mithilfe grafischer Werkzeuge (Plots) für das Modell bp2 ~ bp1.

Lösungen

Nr. Annahme Mathematische Formulierung Bedeutung
1 Fehler mit Erwartungswert null \(\mathbb{E}[\varepsilon_i] = 0\) Das Modell ist im Mittel korrekt spezifiziert.
2 Homoskedastizität \(\text{Var}(\varepsilon_i) = \sigma^2\) Die Fehler haben konstante Varianz (keine Heteroskedastizität).
3 Unabhängigkeit der Fehler \(\varepsilon_i\) unabhängig Die Fehlerterme sind statistisch unabhängig.
4 Normalverteilung der Fehler (optional) \(\varepsilon_i \sim \mathcal{N}(0, \sigma^2)\) Nur für Inferenz (Tests, Konfidenzintervalle) notwendig.
  • Die letzte Annahme (Normalverteilung der Fehler) ist für die Punktschätzung von \(\beta_0\) und \(\beta_1\) (also für die empirische Gerade) nicht nötig, aber für Intervallschätzung und für statistische Tests.
  • Wir schreiben als Annahme für die Fehler oft: \(\epsilon_i \overset{\text{iid}}\sim \mathcal{N}(0, \sigma^2)\), mit iid als unabhängig und gleich verteilt (independent and identically distributed).

Natürlich ist es eine weitere Voraussetzung, dass es sich bei der abhängigen Variablen um eine quantitative Variable handelt und, dass ein linearer Zusammenhang vertretbar ist.

Man kann die Voraussetzungen nicht direkt prüfen, weil wir die wahren Fehler nicht kennen. Stattdessen prüfen wir die Voraussetzungen qualitativ, indem wir die Residuen visualisieren.

Streudiagramm:

plot(df.bp$bp1, df.bp$bp2)

Mithilfe eines Streudiagramms kann man sich die Verteilung der Punkte anschauen. In diesem Fall ist ein linearer Zusammenhang vertretbar (die Punktewolke verteilt sich um eine Linie).

TA-Plot:

plot(lm.bp, which = 1)

Wenn man in R bei der Funktion plot() ein Modell als Objekt verwendet, dann werden verschiedene Grafiken erstellt. Die erste Grafik (darum which = 1) ist ein Residuenplot. Man sieht, dass die Residuen (y-Achse) über den ganzen Bereich der modellierten Werte (x-Achse) gleichmässig und um 0 verteilt sind. Homoskedastizität ist die wichtigste Voraussetzung für ein lineares Regressionsmodell. In diesem Fall gibt es keine Evidenz gegen diese Annahme (der Residuenplot ist praktisch perfekt). Die rote Linie ist eine Hilfslinie. Sie sollte möglichst waagrecht und nahe bei 0 sein.

QQ-Plot:

plot(lm.bp, which = 2)

Wie oben erhält man mit which = 2 einen QQ-Plot für die Residuen. Die Grafik weist nicht darauf hin, dass die Residuen von einer Normalverteilung abweichen. Somit gibt es auch gegen diese Voraussetzung keine Evidenz (sie ist aber bei Weitem nicht so wichtig wie die Homoskedastizität).

Ob die Fehler \(\epsilon_i\) unabhängig sind, lässt sich zum Beispiel anhand des Studiendesigns beurteilen.

Ausblick

Vielleicht hast du dich gefragt, warum es sinnvoll sein soll, das Modell bp2 ~ bp1 zu betrachten, wenn dabei gar keine Aussage über den Effekt des Medikaments gemacht werden kann. Es wäre doch viel interessanter, auch noch die Variable drug in das Modell aufzunehmen. Tatsächlich enthalten lineare Regressionsmodelle häufig mehrere Prädiktorvariablen. Solche Modelle nennt man dann nicht mehr einfache, sondern multiple lineare Regressionen. Als Vorgeschmack folgt hier das obige Modell, ergänzt um die Variable drug.

lm.bp1.drug <- lm(bp2 ~ bp1 + drug, df.bp)
summary(lm.bp1.drug)
## 
## Call:
## lm(formula = bp2 ~ bp1 + drug, data = df.bp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.156  -3.057  -0.206   2.513  13.234 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -10.90899    5.16625  -2.112   0.0374 *  
## bp1           1.09497    0.04247  25.779  < 2e-16 ***
## drugyes      -5.99613    0.92614  -6.474  4.3e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.531 on 94 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.8781, Adjusted R-squared:  0.8755 
## F-statistic: 338.6 on 2 and 94 DF,  p-value: < 2.2e-16

Wie man multiple lineare Regressionsmodelle erstellt und interpretiert, werden wir am nächsten Modultag besprechen.