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 Veränderung von bp1 nach bp2 (definiert als post-prä, also bp2-bp1). Speichere die Ergebnisse in der Variable bp_change und füge diese dem Datensatz hinzu.

Lösungen

  1. Laden und betrachten des Datensatzes:
library("rio")
df.bp <- import("../data/blood-pressure.csv",stringsAsFactors=TRUE)
## mit stringsAsFactors=TRUE (default ist FALSE) werden character vars. als Faktoren mit ihren Levels deklariert  
str(df.bp)
## 'data.frame':    100 obs. of  9 variables:
##  $ pid   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ gender: Factor w/ 2 levels "female","male": 1 NA 2 1 2 1 2 2 1 2 ...
##  $ drug  : Factor w/ 2 levels "no","yes": 2 1 1 2 1 1 1 2 2 2 ...
##  $ 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            weight     
##  Min.   :  1.00   female:53   no  :50   Min.   :28.90   Min.   :53.00  
##  1st Qu.: 25.75   male  :46   yes :48   1st Qu.:45.40   1st Qu.:63.75  
##  Median : 50.50   NA's  : 1   NA's: 2   Median :51.40   Median :68.50  
##  Mean   : 50.50                         Mean   :51.23   Mean   :69.76  
##  3rd Qu.: 75.25                         3rd Qu.:56.10   3rd Qu.:76.00  
##  Max.   :100.00                         Max.   :79.00   Max.   :91.00  
##                                         NA's   :1                      
##      height           bmi             bp1             bp2       
##  Min.   :154.0   Min.   :17.20   Min.   : 97.0   Min.   : 92.0  
##  1st Qu.:164.0   1st Qu.:23.10   1st Qu.:115.5   1st Qu.:111.0  
##  Median :169.0   Median :24.20   Median :121.0   Median :119.0  
##  Mean   :169.3   Mean   :24.23   Mean   :121.8   Mean   :119.8  
##  3rd Qu.:174.5   3rd Qu.:25.10   3rd Qu.:128.0   3rd Qu.:128.0  
##  Max.   :188.0   Max.   :32.50   Max.   :155.0   Max.   :160.0  
##  NA's   :1       NA's   :1       NA's   :1       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 Blutdruckveränderung:
df.bp$bp_change <- df.bp$bp2 - df.bp$bp1
summary(df.bp$bp_change)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -15.000  -5.500  -2.000  -2.284   2.000  16.000

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_change). Lass dir davon eine Zusammenfassung anzeigen. Interpretiere die Koeffizienten und deren P-Werte.

Lösung

lm.bp <- lm(bp_change ~ drug, data = df.bp)
summary(lm.bp)
## 
## Call:
## lm(formula = bp_change ~ drug, data = df.bp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9.8085 -3.5625  0.1915  2.4375 15.4375 
## 
## 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 Blutdruckveränderung in der Gruppe ohne Medikament. In dieser Gruppe hat sich der Blutdruck also um 0.5625 mmHg erhöht. Der P-Wert ist gross. Somit besteht keine Evidenz gegen \(H_0\). \(\hat{\beta}_1\) ist die Mittelwertsdifferenz der Veränderungen. Die Blutdruckveränderung ist in der Gruppe mit Medikament kleiner, was einer grösseren Reduktion entspricht ( -5.7539894). 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_change erklären kann (siehe R-squared). Um den “Fit” des Modells zu verbessern, also den Anteil der erklärten Varianz von bp_change 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 drug 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 lautet die 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 Blutdruckveränderung für folgende Person:
    • age = 60
    • bmi = 29
    • drug = “yes”

Lösungen

Zuerst wird das entsprechende Modell angepasst

lm.bp.multi <- lm(bp_change ~ drug + bmi + age, data = df.bp)
summary(lm.bp.multi)
## 
## Call:
## lm(formula = bp_change ~ drug + bmi + age, data = df.bp)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.9500  -3.6027  -0.0425   2.8668  14.2940 
## 
## 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_change 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\): Das kleinere Modell ist wahr.

anova(lm.bp, lm.bp.multi)
## Analysis of Variance Table
## 
## Model 1: bp_change ~ drug
## Model 2: bp_change ~ 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_{model2}} = \frac{1794}{91} = 19.71\),

\(MS_{Res}^{model1 - model2} = \frac{SS_{Res}^{model1} - SS_{Res}^{model2}}{df_{model1}-df_{model2}} = \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 drug zu interpretieren?

Wenn die Variable drug die Ausprägung 1 annimmt (was heisst, dass die Person das Medikament nimmt), dann wird die abhängige Variable bp_change um 5.3337722 kleiner, unter Berücksichtigung von bmi und age (negative Steigung, mehr Blutdruckreduktion als die Kontrollgruppe). Inhaltlich heisst das, dass die Blutdruckreduktion mit Medikament um 5.33 mmHg höher ist. Dieser Koeffizient (oder diese Differenz der Veränderungen) unterscheidet sich statistisch signifikant von 0.

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

Wenn die Variable bmi um 1 steigt, dann steigt der Wert der abhängigen Variable bp_change um 0.51 mmHg, unter Berücksichtigung von age und drug (positive Steigung). Inhaltlich heisst das, dass die Blutdruckveränderung mit zunehmendem BMI zunimmt (oder die Reduktion abnimmt). Dieser Koeffizient unterscheidet sich statistisch signifikant von 0.

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

Wenn die Variable age um 1 steigt, dann wird der Wert der abhängigen Variable bp_change um 0.14 kleiner, unter Berücksichtigung von drug und bmi (negative Steigung). Inhaltlich heisst das, dass die Blutdruckveränderung mit zunehmendem Alter sinkt (oder die Reduktion steigt). Dieser Koeffizient unterscheidet sich statistisch signifikant von 0.

Die drei Effekte sind hier noch graphisch dargestellt:

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 verschieden als im einfachen Modell?

Eine der beiden anderen Prädiktorvariablen age und bmi (oder beide) können mit der Variable drug und mit dem Outcome korrelieren und sind damit potentielle Confounder. In der multiplen Regression ist der Effekt von drug der Effekt, der über die Effekte von age und bmi hinausgeht, es wird für age und bmi adjustiert. Der Koeffizient für drug von -0.533 ist also der Koeffizient unter Berücksichtigung des Effekts 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 verteilt 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_change 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.