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.

Lösungen

Laden und Betrachten des Datensatzes:

library("rio")
d.bp <- import("../../data/blood-pressure.csv")
str(d.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(d.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

Beachte, dass es fehlende Werte gibt!

ANOVA vs lineare Modelle

Was ist eigentlich der Unterschied zwischen einem linearen Modell und einer ANOVA (Analysis Of Variance)? Mathematisch sind die beiden “Methoden” äquivalent. Was sich unterscheidet sind die R-Codes und die R-Outputs. Grundsätzlich kann man nichts falsch machen wenn man immer die lm() Funktion verwendet.

Wenn man Mittelwerte von mehr als zwei Gruppen vergleicht (bei zwei nehmen wir ja den t-Test), dann sprechen wir typischerweise von ANOVA. In diesem Fall ist die unabhängige Variable die Gruppenzugehörigkeit, also nominal skaliert. Bei einer ANOVA wird die unabhängige Variable oft auch als Faktor bezeichnet. Ist die unabhängige Variable quantitativ skaliert, spricht man eher von linearen Modellen. In den folgenden Aufgaben werden die Gemeinsamkeiten von ANOVA und multipler Regression noch einmal verdeutlicht.

Aufgaben

  1. Erstelle eine neue Variable bmi_cat. Für Personen, welche einen BMI < 18.5 haben, soll diese den Ausprägungsgrad Underweight annehmen. Die weiteren Ausprägungen sind Normal (18.5 \(\le\) BMI <25) und Overweight (BMI \(\ge\) 25). Füge die Variable bmi_cat dem Datensatz hinzu.
  2. Ist der mittlere Blutdruck (bp1) in allen drei Gruppen der Gleiche (\(H_0\)) oder unterscheidet sich mindestens einer der Mittelwerte (\(H_A\)) ? Stelle den BMI der drei Gruppen mittels Boxplots dar.
  3. Berechne den Einfluss von bmi_cat auf bp1 mit lm(). Lass dir von R eine Zusammenfassung ausgeben. Lass dir ebenfalls die Koeffizienten und deren 95% CIs anzeigen.
  4. Berechne den Einfluss von bmi_cat auf bp1 mit aov(). Lass dir von R eine Zusammenfassung ausgeben. Lass dir ebenfalls die Koeffizienten und deren 95% CIs anzeigen.
  5. Erläutere, inwiefern sich die Outputs unterscheiden.

Lösungen

Zuerst wird die Variable bmi_cat erstellt (auch hier gäbe es wieder viele verschiedene Möglichkeiten).

d.bp$bmi_cat[d.bp$bmi<=18.5] <- "Underweight"
d.bp$bmi_cat[d.bp$bmi>18.5 & d.bp$bmi<25] <- "Normal"
d.bp$bmi_cat[d.bp$bmi>=25] <- "Overweight"
d.bp$bmi_cat <- factor(d.bp$bmi_cat)
table(d.bp$bmi_cat)
## 
##      Normal  Overweight Underweight 
##          66          31           2

Wir sehen, dass die Kategorie Underweight nur zweimal vorkommt. Zudem sehen wir, dass bei einer Person die Angabe zum BMI fehlt.

Boxplots:

boxplot(bp1~bmi_cat, data = d.bp)

Der Boxplot macht für die Gruppe Underweight keinen Sinn. Zur Erinnerung: Die Linie repräsentiert nicht den Mittelwert sondern den Median.

Anpassen des linearen Modells:

lm.bp <- lm(bp1 ~ bmi_cat, data = d.bp)
summary(lm.bp)
## 
## Call:
## lm(formula = bp1 ~ bmi_cat, data = d.bp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.323  -6.879  -0.161   5.677  36.121 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         118.879      1.252  94.944  < 2e-16 ***
## bmi_catOverweight     9.444      2.215   4.264 4.74e-05 ***
## bmi_catUnderweight   -3.879     10.249  -0.378    0.706    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.17 on 95 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.164,  Adjusted R-squared:  0.1464 
## F-statistic: 9.318 on 2 and 95 DF,  p-value: 0.0002018

Der globale F-Test ist mit einem P-Wert von 0.0002018 statistisch signifikant (F-Statistik: 9.318, df = 2, 95). In der Gruppe Normal beträgt der mittlere Blutdruck 118.879. In der Gruppe Overweight ist er 9.444 höher. Diese Differenz ist statistisch signifikant. Der Blutdruck in der Gruppe Underweight ist 3.879 tiefer (Minuszeichen beachten). Diese Differenz ist statistisch nicht signifikant (p = 0.706, beachte: n in dieser Gruppe ist tief…). Anhand der Variable bmi_cat lässt sich 16.4% der Varianz von bp1 erklären.

95% CIs:

confint(lm.bp)
##                         2.5 %    97.5 %
## (Intercept)        116.393057 121.36452
## bmi_catOverweight    5.046766  13.84082
## bmi_catUnderweight -24.225371  16.46780

Anpassung des Modells mit aov():

aov.bp <- aov(bp1 ~ bmi_cat, data = d.bp)
summary(aov.bp)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## bmi_cat      2   1928   964.1   9.318 0.000202 ***
## Residuals   95   9830   103.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2 observations deleted due to missingness

Der P-Wert für das gesamte Modell lässt sich in der Zeile für bmi_cat ablesen. F-Wert (inkl. den Freiheitsgraden) und P-Wert sind genau gleich, wie im linearen Modell oben (kleiner Unterschied wegen dem Runden). Hier ist ersichtlich, dass der lm()-Funktion-Output etwas mehr Fleisch am Knochen hat. Man kann diesen aber auch aus einem aov-Objekt herauskitzeln:

summary.lm(aov.bp)
## 
## Call:
## aov(formula = bp1 ~ bmi_cat, data = d.bp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.323  -6.879  -0.161   5.677  36.121 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         118.879      1.252  94.944  < 2e-16 ***
## bmi_catOverweight     9.444      2.215   4.264 4.74e-05 ***
## bmi_catUnderweight   -3.879     10.249  -0.378    0.706    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.17 on 95 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.164,  Adjusted R-squared:  0.1464 
## F-statistic: 9.318 on 2 and 95 DF,  p-value: 0.0002018
confint(aov.bp)
##                         2.5 %    97.5 %
## (Intercept)        116.393057 121.36452
## bmi_catOverweight    5.046766  13.84082
## bmi_catUnderweight -24.225371  16.46780

Fazit: Ob ANOVA oder multiple lineare Regression, ist Geschmackssache und hängt primär davon ab, welche Informationen man braucht. Möchte man einzelne Faktoren testen (siehe weiter unten), bietet der Output von aov() gewisse Vorteile. Stehen einzelne Koeffizienten im Zentrum, dass liefert die lm() Funktion den besseren Output.

Sum of Squares, Mean Squares

Betrachte für die folgenden Aufgaben den aov-Output des Modells oben (bp1 ~ bmi_cat)

Aufgaben

  1. Die Mean Squares (MS oder Mean Sq) des Modells und der Residuen berechnen sich jeweils aus \(SS/df\). Wie lauten die \(MS_{model}\) und die \(MS_{Residuals}\) und was sind deren entsprechende Freiheitsgrade?
  2. Rechne den F-Wert manuell nach und bestimme die entsprechenden Freiheitsgrade.
  3. Mit der Funktion pf(F, df1, df2, lower.tail = FALSE) kannst du dir den P-Wert von R berechnen lasssen. Wie gross ist dieser?
  4. Die Sum of Squares (SS) des Modells und die SS der Residuen ergeben die totalen SS:
SS_tot <- 1928 + 9830
SS_tot
## [1] 11758

Die \(SS_{total}\) entsprechen der gesamten Variabilität der abhängigen Variable (bp1). Die \(SS_{model}\) entsprechen der Variabilität, welche durch das Modell erklärt werden kann. Der Rest ist “Error”: \(SS_{residuals}\). Wir können berechnen, wie viel der gesamten Variabilität durch das Modell erklärt werden kann:

\[\frac{SS_{model}}{SS_{total}}\]

Berechne dieses Ergebnis und interpretiere es.

Lösungen

summary(aov.bp)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## bmi_cat      2   1928   964.1   9.318 0.000202 ***
## Residuals   95   9830   103.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2 observations deleted due to missingness

Die Mean Squares berechnen sich wie folgt:

\(MS_{model} = \frac{SumSq_{bmi.cat}}{df} = \frac{1928}{2} = 964\), mit \(df = 2\)

\(MS_{residuals} = \frac{SumSq_{residuals}}{df} = \frac{9830}{95} = 103.5\), mit \(df = 95\)

Mit den Mean Squares kann der F-Wert und damit der P-Wert berechnet werden:

\(F-Wert = MS_{model}/MS_{residuals}\) mit den Freiheitsgraden \(df_{model}\) und \(df_{residuals}\),

\(\frac{MS_{model}}{MS_{residuals}} = \frac{964}{103.5} = 9.31401\), mit \(df_{model} = 2\) und \(df_{residuals} = 95\).

pf(9.31401, 2, 95, lower.tail = FALSE)
## [1] 0.0002024477

Die Summe aus \(SS_{model}\) und \(SS_{residuals}\) ergibt die \(SS_{total}\)

SS_tot <- 1928 + 9830
SS_tot
## [1] 11758

Die \(SS_{total}\) entsprechen der gesamten Variabilität der abhängigen Variable (bp1). Die \(SS_{model}\) entspricht der Variabilität, welche durch das Modell erklärt werden kann. Der Rest ist “Error” (\(SS_{residuals}\)). Wir können berechnen, wie viel der gesamten Variabilität durch das Modell erklärt werden kann:

\(SS_{model}/SS_{total} = 1928/11758 = 0.164\)

1928 / SS_tot
## [1] 0.1639735

Das bedeutet, das Modell kann 16.4% der gesamten Variabilität der abhängigen Variable erklären. Der Wert ist nichts anderes als R-squared (vgl. Output zum multiplen linearen Modell).

ANOVA mit mehr als einem Faktor

Es ist möglich, mehr als eine qualitative Prädiktorvariable in das Modell einzuschliessen. Man spricht dann nicht mehr von einer One-Way ANOVA (1 Faktor) sondern von Two-Way-ANOVA (2 Faktoren) bzw. von factorial designs (mehrere Faktoren).

Aufgabe

  1. Passe das Modell bp1 ~ bmi_cat + gender an. Einmal mit lm() und einmal mit aov().

  2. Interpretiere beide Outputs. Was sind Gemeinsamkeiten, was sind Unterschiede?

  3. Google (oder frage Chat GPT), wie die Daten nach BMI-Gruppe und nach Geschlecht darstellen kannst.

Lösungen

Anpassen des Modells:

aov.bp.2 <- aov(bp1 ~ bmi_cat + gender, data = d.bp)
lm.bp.2 <- lm(bp1 ~ bmi_cat + gender, data = d.bp)
summary(aov.bp.2)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## bmi_cat      2   1928   964.1   9.991 0.000116 ***
## gender       1    759   758.7   7.862 0.006134 ** 
## Residuals   94   9071    96.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2 observations deleted due to missingness
summary(lm.bp.2)
## 
## Call:
## lm(formula = bp1 ~ bmi_cat + gender, data = d.bp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -21.553  -5.953   0.047   4.958  31.691 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         116.953      1.391  84.096  < 2e-16 ***
## bmi_catOverweight     6.244      2.424   2.575  0.01157 *  
## bmi_catUnderweight   -1.953      9.921  -0.197  0.84440    
## gendermale            6.356      2.267   2.804  0.00613 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.824 on 94 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.2285, Adjusted R-squared:  0.2039 
## F-statistic: 9.281 on 3 and 94 DF,  p-value: 1.951e-05

Wir sehen wiederum die Gemeinsamkeiten und Unterschiede, welche oben schon besprochen wurden. Was nun anders ist, ist der P-Wert für bmi_cat. Im linearen Modell erhalten wir jeweils die P-Werte (und 95% CIs via confint()) für die Mittelwertsdifferenz zwischen dem Referenzlevel (Normal) und den anderen Levels (unter Berücksichtigung der BMI-Kategorie). Im aov-Output hingegen erhalten wir den P-Wert für den Effekt von bmi_cat insgesamt und nicht für einzelne Vergleiche. Man nennt diese Gesamteffekte eine Faktors main effects. Es kommt daher auf die Forschungsfrage an, welcher Output geeigneter ist.

Bei gender sind die P-Werte der beiden Modelle wieder identisch. Warum? Immer dann, wenn ein Prädiktor nur 1 Freiheitsgrad hat, sind die P-Werte der beiden Outputs gleich. gender hat nur zwei Levels und weil \(df = \#levels - 1\), eben nur ein Freiheitsgrad.

Technische Anmerkung: Falls du statt bp1 ~ bmi_cat + gender, bp1 ~ gender + bmi_cat angepasst hast, bekommst du mit aov() andere P-Werte. Bei aov() ist die Reihenfolge der Prädiktoren wichtig. F-Tests in aov() und anova() sind sequentiell (Reihenfolge ist wichtig), t-Tests in summary.lm() bzw. lm() sind marginal (Reihenfolge spielt keine Rolle).

Darstellung der Daten

Einmal mehr gibt es verschiedene Möglichkeiten. Eine ist, das viel gelobte Package ggplot2 zu verwenden. Dieses Packages mag Faktoren, weshalb diese zuerst noch erstellt werden.

library(tidyverse)
d.bp$gender <- factor(d.bp$gender)
d.bp <- na.omit(d.bp) # NAs rausnhemen

library(ggplot2)
ggplot(d.bp, aes(x = bmi_cat, y = bp1, fill = gender)) +
  geom_boxplot()