Zum Datensatz

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

Codebuch:

  • gender: Geschlecht
  • age: Alter
  • height: Körpergrösse
  • weight: Gewicht
  • bmi: Body-Mass-Index in kg/m²
  • drug: Medikamenteneinnahme (no = nein, yes = ja)
  • bp1 und bp2: Blutdruck in mmHg vor (bp1) und nach (bp2) dem Beobachtungszeitraum

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 minus 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

Insgesamt haben 5 Personen mindestens einen fehlenden Wert; 95 vollständige Beobachtungen bleiben übrig. Eine Analyse vollständiger Fälle kann bei systematisch fehlenden Werten verzerrte Schätzungen liefern.

  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. Das ist 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 des Tests auf einen Achsenabschnitt von 0 beträgt ungefähr 0.41. Die Daten liefern somit keine ausreichende Evidenz gegen die Nullhypothese, dass die mittlere Blutdruckveränderung in der Gruppe ohne Medikament 0 ist. \(\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 Wert unterscheidet sich statistisch signifikant von 0.

Einfache vs. multiple lineare Regression

Das einfache Modell erklärt 27.91% der beobachteten Varianz von bp_change (Bestimmtheitsmass \(R^2\)). Mit zusätzlichen Prädiktoren erhalten wir eine multiple lineare Regression. Das unadjustierte \(R^2\) kann beim Hinzufügen von Prädiktoren auf denselben Beobachtungen nicht abnehmen; ein Anstieg allein belegt daher keine bessere Vorhersage für neue Beobachtungen.

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 einfachen 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 0.2791 auf 0.3632; das multiple Modell erklärt damit 36.32% der beobachteten Varianz. Das adjustierte \(R^2\) berücksichtigt zusätzlich die Anzahl geschätzter Parameter und steigt hier von ungefähr 0.2713 auf 0.3422. Auch dieser Anstieg belegt allein keine bessere Vorhersage für neue Beobachtungen.

Mit einem partiellen F-Test vergleichen wir die beiden verschachtelten Modelle, die auf denselben Beobachtungen angepasst wurden. Die Nullhypothese lautet \(H_0: \beta_{bmi} = \beta_{age} = 0\); der Koeffizient von drug bleibt dabei frei.

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 Residuenquadratsumme (RSS) ist 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_{Res}^{model2} = \frac{RSS_{model2}}{df_{Res,model2}} = \frac{1794}{91} = 19.71\),

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

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

wobei \(MS\) für mittlere Quadratsumme steht. Im Nenner steht die mittlere Residuenquadratsumme des grösseren Modells, im Zähler die Abnahme der Residuenquadratsumme pro zusätzlichem Parameter.

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

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

Für Personen mit drug = "yes" ist die geschätzte mittlere Blutdruckveränderung um 5.3337722 mmHg kleiner als für Personen mit drug = "no" bei gleichem BMI und Alter. 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?

Ein um 1 kg/m² höherer BMI ist bei gleichem Alter und gleicher Medikamenteneinnahme mit einer um 0.51 mmHg höheren geschätzten mittleren Blutdruckveränderung assoziiert. 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?

Ein um ein Jahr höheres Alter ist bei gleichem BMI und gleicher Medikamenteneinnahme mit einer um 0.14 mmHg niedrigeren geschätzten mittleren Blutdruckveränderung assoziiert. 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 adjustierten Zusammenhänge sind hier grafisch dargestellt:

Was ist der p-Wert des globalen F-Tests?

Der p-Wert des globalen F-Tests steht unten im Output und beträgt ungefähr \(5.68 \times 10^{-9}\). Der Test prüft \(H_0: \beta_{drug} = \beta_{bmi} = \beta_{age} = 0\) gegen die Alternative, dass mindestens einer dieser Koeffizienten nicht 0 ist. Er vergleicht das multiple Modell mit einem Modell, das nur einen Achsenabschnitt enthält. Unter \(H_0\) erklären die Prädiktoren keine Varianz in der Population; das geschätzte \(R^2\) in der Stichprobe muss deshalb nicht 0 sein.

Warum ist der Regressionskoeffizient für drug im multiplen Modell anders als im einfachen Modell?

Das einfache Modell schätzt den unadjustierten Unterschied der mittleren Blutdruckveränderung zwischen den Medikamentengruppen. Das multiple Modell schätzt diesen Unterschied bei gleichem Alter und BMI; es wird für age und bmi adjustiert. Sind diese Prädiktoren sowohl mit drug als auch mit bp_change assoziiert, kann sich der Koeffizient durch die Adjustierung ändern. Der adjustierte Koeffizient für drug beträgt hier ungefähr −5.33 mmHg. Ob Alter oder BMI Confounder sind und ob der Koeffizient kausal interpretiert werden kann, hängt zusätzlich vom Studiendesign und den kausalen Beziehungen ab; Korrelationen allein reichen dafür nicht aus.

Schätzung der erwarteten Blutdruckveränderung für eine Person mit folgenden Merkmalen:

  • 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

Wir betrachten den Tukey-Anscombe-Plot (Residuen gegen vorhergesagte Werte) und den QQ-Plot der Residuen:

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

Der TA-Plot zeigt keine auffällige Veränderung der Residuenstreuung, und der QQ-Plot zeigt keine ausgeprägten Abweichungen von der Normalverteilung. Diese Befunde beweisen die Annahmen jedoch nicht. Die geglättete rote Linie im TA-Plot weist eine leichte Krümmung auf, die auf eine unzureichend modellierte Form des bedingten Mittelwerts hindeuten könnte. Aus diesem Plot allein lässt sich nicht bestimmen, ob age oder bmi dafür verantwortlich ist. Quadratische Terme könnten bei entsprechender fachlicher Begründung geprüft werden; das Modell bliebe linear in seinen Koeffizienten. Die Unabhängigkeit der Fehler muss anhand des Studiendesigns beurteilt werden. Zusätzlich sollten einflussreiche Beobachtungen untersucht werden.