Zum Datensatz

Für diese Übung brauchst du den Datensatz “blood-pressure.csv”.

Code-Book:

  • gender: Geschlecht
  • age: Alter
  • height: Körpergrösse
  • weight: Gewicht
  • drug: ob die Person ein Medikament einnahm (drug)
  • bp1 und bp2: Blutdruck zu zwei verschiedenen Messzeitpunkten (Prä/Post)

Aufgaben

  1. Lade und inspiziere den Datensatz.
  2. Entferne alle Beobachtungen mit fehlenden Werten. Tipp: na.omit()
    Wie viele Beobachtungen bleiben übrig?
  3. Berechne für jede Person die Differenz zwischen bp1 und bp2, was der Blutdruckreduktion entspricht. Speichere die Ergebnisse in der Variable bp_red und füge diese dem Datensatz hinzu. Für die Lösungen wurde bp1 - bp2 gerechnet (umgekehrt wäre aber nicht falsch).

Lösungen

  1. Laden und betrachten des Datensatzes:
library("rio")
df.bp <- import("../data/blood-pressure.csv")
str(df.bp)
## 'data.frame':    100 obs. of  9 variables:
##  $ pid   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ gender: chr  "female" NA "male" "female" ...
##  $ drug  : chr  "yes" "no" "no" "yes" ...
##  $ age   : num  NA 56.5 45.9 64 58 52.9 77.8 56.9 49.9 65.3 ...
##  $ weight: int  66 82 65 59 55 75 84 68 66 68 ...
##  $ height: int  173 NA 165 160 174 182 169 177 169 175 ...
##  $ bmi   : num  22.1 NA 23.9 23 18.2 22.6 29.4 21.7 23.1 22.2 ...
##  $ bp1   : int  129 122 109 126 NA 106 128 144 118 132 ...
##  $ bp2   : int  123 122 105 116 139 109 133 133 117 120 ...
summary(df.bp)
##       pid            gender              drug                age       
##  Min.   :  1.00   Length:100         Length:100         Min.   :28.90  
##  1st Qu.: 25.75   Class :character   Class :character   1st Qu.:45.40  
##  Median : 50.50   Mode  :character   Mode  :character   Median :51.40  
##  Mean   : 50.50                                         Mean   :51.23  
##  3rd Qu.: 75.25                                         3rd Qu.:56.10  
##  Max.   :100.00                                         Max.   :79.00  
##                                                         NA's   :1      
##      weight          height           bmi             bp1       
##  Min.   :53.00   Min.   :154.0   Min.   :17.20   Min.   : 97.0  
##  1st Qu.:63.75   1st Qu.:164.0   1st Qu.:23.10   1st Qu.:115.5  
##  Median :68.50   Median :169.0   Median :24.20   Median :121.0  
##  Mean   :69.76   Mean   :169.3   Mean   :24.23   Mean   :121.8  
##  3rd Qu.:76.00   3rd Qu.:174.5   3rd Qu.:25.10   3rd Qu.:128.0  
##  Max.   :91.00   Max.   :188.0   Max.   :32.50   Max.   :155.0  
##                  NA's   :1       NA's   :1       NA's   :1      
##       bp2       
##  Min.   : 92.0  
##  1st Qu.:111.0  
##  Median :119.0  
##  Mean   :119.8  
##  3rd Qu.:128.0  
##  Max.   :160.0  
##  NA's   :2

Um die anschliessenden Analysen etwas zu vereinfachen, verwenden wir nur Beobachtungen mit kompletten Werten. Dieses Vorgehen ist didaktisch und nicht inhaltlich begründet.

df.bp <- na.omit(df.bp)
nrow(df.bp)
## [1] 95

Man sieht, dass insgesamt 5 Personen fehlende Werte hatten.

  1. Berechnung der Blutdruckreduktion:
df.bp$bp_red <- df.bp$bp1 - df.bp$bp2

Einfache lineare Regression

Aufgabe

Berechne zuerst eine einfache lineare Regression. Also ein Regressionsmodell, mit nur einer Prädiktorvariable (drug) und einer abhängigen Variable (bp_red). Lass dir davon eine Zusammenfassung anzeigen. Interpretiere die Koeffizienten und deren P-Werte.

Lösung

lm.bp <- lm(bp_red ~ drug, data = df.bp)
summary(lm.bp)
## 
## Call:
## lm(formula = bp_red ~ drug, data = df.bp)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.4375  -2.4375  -0.1915   3.5625   9.8085 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.5625     0.6745  -0.834    0.406    
## drugyes       5.7540     0.9590   6.000 3.75e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.673 on 93 degrees of freedom
## Multiple R-squared:  0.2791, Adjusted R-squared:  0.2713 
## F-statistic:    36 on 1 and 93 DF,  p-value: 3.754e-08

Das Intercept entspricht der durchschnittlichen Blutdruckreduktion in der Gruppe ohne Medikament. In dieser Gruppe hat sich der Blutdruck also um 0.55 mmHg erhöht. Der P-Wert ist gross. Somit besteht keine Evidenz gegen \(H_0\). \(b_1\) ist die Mittelwertsdifferenz zwischen den beiden Gruppen. Die Blutdruckreduktion ist in der Gruppe mit Medikament also 5.76 mmHg höher. Dieser Werte unterscheidet sich statistisch signifikant von 0.

Einfache vs. multiple lineare Regression

Das Modell oben sagt dir, dass die Variable drug 27.91% der Varianz von bp_red erklären kann (siehe R-squared). Um den “Fit” des Modells zu verbessern, also den Anteil der erklärten Varianz von bp_red zu steigern, können zusätzliche Variablen ins Modell genommen werden. Aus der simple linear regression wird also eine multiple linear regression.

Hier noch einmal die Form eines einfachen und eines multiplen Regressionsmodells:

Einfache lineare Regression:

\[ Y_i = \beta_0 + \beta_1 x_i + \epsilon_i, i = 1, ..., n, \]

Multiple lineare Regression:

\[ Y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2}+ ... + \beta_mx_{im} + \epsilon_i, i = 1, ..., n, \]

Aufgaben

  1. Füge dem einfachen Modell die Variablen bmi und age hinzu. Lass dir eine Zusammenfassung dieses Modells ausgeben.
  2. Was passiert mit R-squared im Vergleich zum einfachen Modell und was bedeutet das?
  3. Ist das multiple Modell signifikant besser als das einfache? Vergleiche die beiden Modelle. Tipp: anova()
  4. Wie ist der Regressionskoeffizient der Prädiktorvariable durg zu interpretieren?
  5. Wie ist der Regressionskoeffizient der Prädiktorvariable bmi zu interpretieren?
  6. Wie ist der Regressionskoeffizient der Prädiktorvariable age zu interpretieren?
  7. Was ist der P-Wert des globalen F-Tests? Wie lautetdie Hypothese, welche damit getestet wird?
  8. Warum ist der Regressionskoeffizient für drug im multiplen Modell ein anderer als im simplen Modell?
  9. Schätze die Blutdruckreduktion für folgende Person:
    • age = 60
    • bmi = 29
    • drug = “yes”

Lösungen

Zuerst wird das entsprechende Modell angepasst

lm.bp.multi <- lm(bp_red ~ drug + bmi + age, data = df.bp)
summary(lm.bp.multi)
## 
## Call:
## lm(formula = bp_red ~ drug + bmi + age, data = df.bp)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.2940  -2.8668   0.0425   3.6027  11.9500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.02228    5.89225   0.852  0.39626    
## drugyes      5.33377    0.92961   5.738 1.25e-07 ***
## bmi         -0.51083    0.23182  -2.204  0.03008 *  
## age          0.13763    0.04769   2.886  0.00487 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.44 on 91 degrees of freedom
## Multiple R-squared:  0.3632, Adjusted R-squared:  0.3422 
## F-statistic:  17.3 on 3 and 91 DF,  p-value: 5.676e-09

\(R^2\) steigt von .2791 auf .3632. Somit hat sich der “Fit” des Modells verbessert. Weil mehr Prädiktorvariablen ins Modell genommen wurden, kann ein grösserer Anteil der Varianz von bp_red erklärt werden (neu 36.32%). Wenn wir Modelle mit mehreren Prädiktoren betrachten, interpretieren wir in der Regel das Adjusted R-squared.

Wir können testen, ob sich die beiden Modelle unterscheiden. \(H_0\): Die Modelle unterscheiden sich nicht.

anova(lm.bp, lm.bp.multi)
## Analysis of Variance Table
## 
## Model 1: bp_red ~ drug
## Model 2: bp_red ~ drug + bmi + age
##   Res.Df    RSS Df Sum of Sq     F   Pr(>F)   
## 1     93 2031.1                               
## 2     91 1794.0  2    237.12 6.014 0.003523 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Die Summe der quadrierten Residuen RSS sind im komplexeren Modell rund \(2031.1 - 1794 = 237.1\) tiefer. Der entsprechende P-Wert ist kleiner als 5% und darum ziehen wir die Schlussfolgerung, dass das komplexe Modell die abhängige Variable statistisch signifikant besser erklärt, als das einfache Modell. Den P-Wert könnte man auch manuell berechnen:

\(MS_{tot}^{model2} = \frac{SS_{tot}^{model2}}{df_{tot}^{model2}} = \frac{1794}{91} = 19.71\),

\(MS_{Res}^{model1 - model2} = \frac{SS_{Res}^{model1} - SS_{Res}^{model2}}{df_{tot}^{model1}-df_{tot}^{mode2}} = \frac{237.1}{2} = 118.56\),

und \(F = \frac{MS_{Res}^{model1 - model2}}{MS_{tot}^{model2}} = \frac{118.56}{19.71} = 6.014\), mit \(df = \{2, 91\}\),

wobei \(MS\) für Mean Squares steht.

pf(6.014, 2, 91, lower.tail = FALSE)
## [1] 0.003523048

Wie ist der Regressionskoeffizient der Prädiktorvariable durg zu interpretieren?

Wenn die Variable drug die Ausprägung 1 annimmt (was heisst, dass die Person das Medikament nimmt), dann steigt die abhängige Variable bp_red um 5.33 unter Berücksichtigung von bmi und age (positive Steigung). Inhaltlich heisst das, dass die Blutdruckreduktion mit Medikament um 5.33 mmHg höher ist. Dieser Koeffizient (oder diese Differenz) unterscheidet sich statistisch signifikant von 0.

Wie ist der Regressionskoeffizient der Prädiktorvariable bmi zu interpretieren?

Wenn die Variable bmi um 1 steigt, dann sinkt der Wert der abhängigen Variable bp_red um 0.51 mmHg, unter Berücksichtigung von age und drug (negative Steigung). Inhaltlich heisst das, dass die Blutdruckreduktion mit zunehmendem BMI abnimmt (wenn dieser Effekt linear konstant ist). Dieser Koeffizient von -0.51 unterscheidet sich statistisch signifikant von 0.

Wie ist der Regressionskoeffizient der Prädiktorvariable age zu interpretieren?

Wenn die Variable age um 1 steigt, dann steigt der Wert der abhängigen Variable bp_red um 0.14, unter Berücksichtigung von drug und bmi (positive Steigung). Inhaltlich heisst das, dass die Blutdruckreduktion mit zunehmendem Alter steigt (wenn dieser Effekt linear konstant ist). Dieser Koeffizient von 0.14 unterscheidet sich statistisch signifikant von 0.

Was ist der P-Wert des gesamten Modells/ globaler F-Test?

Hier muss der P-Wert des F-Tests unten rechts im Output betrachtet werden. Dieser P-Wert testet die Null-Hypothese \(\beta_1 = \beta_2 = ... = \beta_m = 0\). Das ist gleichbedeutend mit \(R^2 = 0\).

Warum ist der Regressionskoeffizient für drug im multiplen Modell ein anderer als im simplen Modell?

Eine der beiden anderen Prädiktorvariablen age und bmi (oder beide) scheint mit der Variable drug zu korrelieren. Weil die multiple Regression als Addition der einzelnen Koeffizienten zu verstehen ist, wird dafür adjustiert. Der Koeffizient für drug von 0.533 für ist also der Koeffizient unter Berücksichtigung der Variablen age und bmi. Dieses “Berücksichtigen” ist in Studien sehr hilfreich, weil man so beispielsweise für den BMI korrigieren kann, falls dieser in den beiden Gruppen nicht gleich ist.

Schätzung der durchschnittlichen Blutdruckreduktion für folgende Person:

  • age = 60
  • bmi = 29
  • drug = “yes”

\[ \hat{y} = \hat{\beta_0} + \hat{\beta_1}drug + \hat{\beta_2}bmi + \hat{\beta_3}age \]

\[ = 5.02228+5.33377+(-0.51083)*29+0.13763*60 = 3.8 \]

In R kann man dies wie folgt tun:

predict(lm.bp.multi,
        newdata = data.frame(
          age = 60, 
          bmi = 29, 
          drug = "yes"
          )
        )
##        1 
## 3.800148

Modelldiagnostik

Aufgabe

Erstelle die nötigen diagnostischen Plots der multiplen linearen Regression und interpretiere diese.

Lösung

Wie immer schauen wir uns den Tukey-Anscombe Plot (Residuenplot) und den QQ-Plot mit den Residuen an:

par(mfrow = c(1,2)) # um zwei Grafiken nebeneinander zu erstellen
plot(lm.bp.multi, which = c(1,2))

Sowohl der TA-Plot sowie der QQ-Plot liefern keine offensichtlichen Hinweise, dass die Annahme bzgl. Homoskedastizität oder Normalverteilung der Fehler verletzt sein könnte. Der Smoother (rote Linie) im TA-Plot weist eine leichte Krümmung, konvex nach unten, auf. Dies könnte darauf hinweisen, dass age und/oder bmi nicht linear mit bp_red zusammenhängt/en, sondern dass z.B. eine quadratische Funktion besser geeignet wäre. Es ist möglich, quadratische Terme in ein Modell einzuschliessen, das Modell bleibt dennoch linear.