Zum Datensatz

Für diese Übung brauchst du den Datensatz blood-pressure.csv. Den Datensatz findest du auf Moodle. Der Datensatz beinhaltet Informationen zu:

  • Geschlecht
  • Alter
  • Körpergrösse
  • Gewicht
  • ob die Person ein Medikament einnahm (drug)
  • Blutdruck zu zwei verschiedenen Messzeitpunkten (Prä/Post) in mmHg

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

Korrelation

Aufgaben

  1. Erstelle einen Scatterplot 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. Was ist das 95% CI für den Korrelationskoeffizient nach Pearson? Tipp: cor.test()

Lösungen

Zuerst erstellen wir wieder einen Scatterplot 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
# Weil es NAs hat, werden nur Beobachtungen mit paarweis kompletten Daten berücksichtig.

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

Mit cor.test() erhalten wir zudem ein Konfidenzintervall und ein 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
# Weil es NAs hat, werden nur Beobachtungen mit paarweis kompletten Daten berücksichtig.

Bestimmtheitsmass

Aufgabe

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

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

  3. Was ist der der Zusammenhang zwischen R-squared und dem Korrelationskoeffizient 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 kann mit dem Prädiktor bp1 82.38% der Variabilität der abhängigen Variable bp2 erklärt werden.

Modelldiagnostik

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

Aufgaben

  1. Welches sind die wichtigsten drei Voraussetzungen für eine lineare Regression?

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

Lösungen

Die wichtigsten Voraussetzungen sind:

  1. Ein linearer Zusammenhang ist vertretbar.

  2. Die Fehler \(\epsilon\) sind normalverteilt.

  3. Die Varianz \(\sigma_{\epsilon}^2\) ist konstant über den Wertebereich von \(x\).

Natürlich ist es eine weitere Voraussetzung, dass es sich bei der abhängigen Variable um eine quantitative Variable handelt.

Man kann die drei Voraussetzungen nicht direkt prüfen, weil wir die wahren Fehler nicht kennen. Stattdessen prüfen wir qualitativ, in dem die Residuen visualisiert werden.

Scatterplot:

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

Mit Hilfe eines Scatterplots 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 plot() Funktion 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-Achste) gleichmässig und um 0 verteilt sind. Homoskedastizität ist die wichtigste Voraussetzung für ein lineares Regressinsmodell. In diesem Fall ist sie klar erfüllt (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, einfach mit which = 2, erhält man einen QQ-Plot für die Residuen. Die Grafik weist nicht darauf hin, dass die Residuen von einer Normalverteilung abweichen. Somit ist auch diese Voraussetzung erfüllt (sie ist aber bei weitem nicht so wichtig wie die Homoskedastizität).

Ausblick

Vielleicht haben Sie sich gefragt, warum es sinnvoll sein soll, das Modell bp2 ~ bp1 zu betrachten, wenn dabei gar keine Aussage bzgl. des Effekts des Medikaments gemacht werden kann. Es wäre doch viel interessanter, auch noch die Variable drug in das Modell zu nehmen. In der Tat hat man in linearen Regressionsmodellen viel häufiger mehrere Prädiktorvariablen. Solche Modelle nennt man dann nicht mehr einfache, sondern multiple lineare Regressionen. Als Vorgeschmack hier das obige Modell, zusätzlich mit der Variable durg.

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.