Zum Datensatz

Für diese Übung brauchst du den Datensatz blood-pressure.csv. Den Datensatz findest du auf Moodle. Der Datensatz enthält folgende Angaben:

  • gender: Geschlecht (female oder male)
  • age: Alter in Jahren
  • height: Körpergrösse in cm
  • weight: Gewicht in kg
  • drug: Medikamenteneinnahme (no oder yes)
  • bp1 und bp2: Blutdruck zu zwei verschiedenen Messzeitpunkten (Prä/Post), in mmHg

Aufgaben

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

  2. Gibt es im Datensatz fehlende Werte (NA)? Tipp: is.na() oder summary()

Lösungen

Den Datensatz kannst du mit der Funktion import() aus dem Paket rio importieren. Achte darauf, dass das Arbeitsverzeichnis richtig eingestellt ist.

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

Fehlende Werte

Mit der Funktion is.na() kann für jede Zelle im Datensatz abgefragt werden, ob es sich um einen NA handelt. Falls ja, resultiert ein TRUE, falls nein ein FALSE. In der Computersprache entspricht ein TRUE dem Wert 1, ein FALSE dem Wert 0. Wenn man mit sum() alle TRUEs zusammenzählt, erhält man die Anzahl fehlender Werte (in diesem Fall 9).

sum(is.na(df.bp))
## [1] 9

Man kann auch die Funktion summary() verwenden; dann wird die Anzahl der fehlenden Werte (NA) pro Variable angezeigt. Achtung: Bei Zeichenkettenvariablen funktioniert das nicht → Ein weiterer Grund, warum man solche Variablen in Faktoren umwandeln sollte.

Ohne Umwandlung:

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

Mit Umwandlung:

df.bp$gender <- factor(df.bp$gender)
df.bp$drug <- factor(df.bp$drug)
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

stringsAsFactors=TRUE beim Einlesen: Wenn man in import() dieses Argument setzt (der Standardwert ist FALSE), muss man Zeichenkettenvariablen nicht mehr explizit als Faktoren deklarieren:

df.test <- import("../data/blood-pressure.csv",stringsAsFactors=TRUE)
str(df.test)
## '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.test)
##       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

Es ist wichtig, sich bewusst zu sein, dass fehlende Werte vorhanden sind. Wir müssen für die Analyse aber nichts weiter unternehmen, weil lm() damit umgehen kann.

Standardmässig gilt: na.action = na.omit. Das bedeutet: Jede Zeile, in der entweder die abhängige Variable oder eine unabhängige Variable fehlende Werte (NA) enthält, wird aus der Analyse ausgeschlossen.

Deskriptive Statistik

Aufgabe

Berechne die folgenden deskriptiven Statistiken für die Variable bp1:

  • Mittelwert
  • Median
  • Varianz
  • Standardabweichung
  • Quartile

Lösung

Beachte: Der Datensatz enthält fehlende Werte (NA) → na.rm = TRUE ist deshalb nötig.

mean(df.bp$bp1, na.rm = TRUE)
## [1] 121.8283
median(df.bp$bp1, na.rm = TRUE)
## [1] 121
var(df.bp$bp1, na.rm = TRUE)
## [1] 119.9804
sd(df.bp$bp1, na.rm = TRUE)
## [1] 10.95356
quantile(df.bp$bp1, na.rm = TRUE, probs = c(0.25, 0.5, 0.75))
##   25%   50%   75% 
## 115.5 121.0 128.0

Grafische Darstellung

Aufgabe

Stelle die Variablen bp1 und bp2 als Boxplot dar (am besten zwei Boxplots in einer Grafik).

Lösung

Für eine übersichtlichere Darstellung wird die x-Achse hier manuell beschriftet. Es wird nicht erwartet, dass du das selbstständig tun kannst.

boxplot(df.bp$bp1, df.bp$bp2, xaxt = "none") # xaxt = none --> keine Beschriftung der x-Achse
axis(1, at = c(1,2), labels = c("bp1", "bp2")) # Manuelle Beschriftung der x-Achse (die beiden Befehle müssen zusammen ausgeführt werden)

Berechnung von Differenzen

Aufgaben

  1. Uns interessiert vor allem die Veränderung des Blutdrucks zwischen den beiden Messzeitpunkten. Berechne für jede Person die Veränderung von der ersten zur zweiten Blutdruckmessung. Erstelle eine neue Variable bp_change (bp2 - bp1) und füge sie dem Datensatz hinzu. Was ist die mittlere Veränderung von bp1 zu bp2?
  2. Was bedeutet in diesem Fall eine positive Zahl?

Lösungen

Die neue Variable bp_change (bp2 - bp1) wird wie folgt erstellt. Den Mittelwert können wir mit mean() berechnen.

df.bp$bp_change <- df.bp$bp2 - df.bp$bp1
mean(df.bp$bp_change, na.rm = TRUE)
## [1] -2.298969

Eine positive Zahl heisst, dass bp2 grösser ist als bp1, d. h., dass die Veränderung positiv ist (und umgekehrt). Demnach hat der Blutdruck im Schnitt abgenommen.

t-Test

Aufgaben

  1. Untersuche, ob sich die mittlere Veränderung (also der Mittelwert der Variable bp_change) von 0 unterscheidet (unabhängig von der Medikamenteneinnahme). Welcher statistische Test ist dafür geeignet?
  2. Führe den richtigen Test durch und interpretiere das Resultat. Tipp: t.test()
  3. Unterscheidet sich bp_change zwischen Personen, die ein Medikament eingenommen haben, und denjenigen, die kein Medikament eingenommen haben? Stelle die Veränderung pro Gruppe zuerst mithilfe zweier Boxplots dar.
  4. Überprüfe, ob sich die mittlere Veränderung zwischen den Gruppen unterscheidet. Welcher Test ist für diese Fragestellung geeignet?
  5. Führe den Test durch. Wie ist der Output zu interpretieren? Was ist der Effekt des Medikaments? Tipp: t.test()

Lösungen

In diesem Fall vergleichen wir nicht zwei Gruppen. Somit ist der Einstichproben-t-Test (oder der gepaarte t-Test, was dasselbe ist) angebracht. Beachte, dass \(H_0: \mu_{bp.change} = 0\).

t.test(df.bp$bp_change, mu = 0)
## 
##  One Sample t-test
## 
## data:  df.bp$bp_change
## t = -4.1656, df = 96, p-value = 6.789e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -3.394473 -1.203465
## sample estimates:
## mean of x 
## -2.298969

Die mittlere Veränderung von -2.3 mmHg ist auf einem Signifikanzniveau von 5% statistisch signifikant (p-Wert < 0.001). Unter \(H_0\) ist es also unwahrscheinlich, diese Veränderung oder eine extremere allein durch Zufall zu finden. Das 95%-Konfidenzintervall reicht von -3.4 bis -1.2. Ob die mittlere Veränderung von 2.3 klinisch relevant ist, sagt uns der t-Test nicht. Bei der Interpretation muss berücksichtigt werden, dass in dieser Stichprobe Personen mit und ohne Medikament nicht separat analysiert wurden. Der Test gibt uns also keine Auskunft über den Effekt des Medikaments.

Im zweiten Teil dieser Aufgabe werden zwei Gruppen (mit und ohne Medikament) verglichen. Daher ist in dieser Situation der Zweistichproben-t-Test angebracht. Es ist eine gute Idee, die Analyse zuerst zu visualisieren, bevor man einen statistischen Test durchführt. In diesem Fall bieten sich Boxplots an:

boxplot(df.bp$bp_change ~ df.bp$drug, xlab = "Medication",ylab="change")

Unten siehst du den Output des Zweistichproben-t-Tests. R macht automatisch eine kleine Korrektur (Welch), weil die Varianzen der beiden Populationen als ungleich angenommen werden. Beachte, dass \(H_0: \mu_{yes} = \mu_{no}\).

t.test(df.bp$bp_change~df.bp$drug)
## 
##  Welch Two Sample t-test
## 
## data:  df.bp$bp_change by df.bp$drug
## t = 6.131, df = 94.869, p-value = 1.992e-08
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  3.894398 7.624309
## sample estimates:
##  mean in group no mean in group yes 
##         0.5510204        -5.2083333

In der Gruppe ohne Medikament hat der Blutdruck um 0.55 zugenommen, während der Blutdruck in der Gruppe mit Medikament um 5.21 abgenommen hat. Der geschätzte mittlere Zwischengruppenunterschied der Veränderung ist demnach -5.76 mmHg (“yes” versus “no” [R rechnet “no”-“yes”]). Diese Mittelwertsdifferenz ist auf einem Signifikanzniveau von 5% statistisch signifikant (p-Wert < 0.001). Sollte also \(H_0\) zutreffen, ist ein solches oder ein noch extremeres Ergebnis sehr unwahrscheinlich. Das angegebene 95%-Konfidenzintervall reicht von -7.62 bis -3.9. Wir vertrauen zu 95% darauf, dass die wahre Blutdruckreduktion durch das Medikament, relativ zur Kontrolle, zwischen diesen Grenzen liegt. Ob ein solcher Effekt klinisch relevant ist, sagt uns der Test nicht.

Anmerkung: Wir gehen hier davon aus, dass die Daten nicht durch andere Einflüsse verzerrt sind.

Lineare Regression mit einer binären Prädiktorvariable

Oben hast du mithilfe des t-Tests untersucht, ob die Veränderung bp2 - bp1 davon abhängig ist, ob eine Person ein Medikament nahm oder nicht. Um diese Frage zu beantworten, kann statt eines t-Tests auch ein lineares Regressionsmodell verwendet werden (was mathematisch äquivalent ist).

Aufgaben

  1. Schätze ein lineares Regressionsmodell mit bp_change als abhängige Variable und drug als Prädiktorvariable: \(Y_i = \beta_0 + \beta_1 x_i + \epsilon_i\), wobei \(Y = bp\_change\) und \(x = drug\). Lass dir eine Zusammenfassung dieses Modells anzeigen. Tipp: lm(), summary().
  2. Wie gross ist die Steigung? Warum sollte dir diese Zahl bekannt vorkommen?
  3. Warum sollte dir der t-Wert von 6.132 bekannt vorkommen? Wie lässt sich dieser anhand des Regressionsoutputs berechnen?
  4. Was bedeutet der Achsenabschnitt im berechneten Modell?

Lösungen

Zuerst passen wir das entsprechende Modell an und schauen uns die Zusammenfassung davon an:

lm.drug <- lm(bp_change ~ drug, data = df.bp)
summary(lm.drug)
## 
## Call:
## lm(formula = bp_change ~ drug, data = df.bp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.792 -3.551 -0.551  2.449 15.449 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.5510     0.6607   0.834    0.406    
## drugyes      -5.7594     0.9392  -6.132 1.98e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.625 on 95 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.2836, Adjusted R-squared:  0.276 
## F-statistic:  37.6 on 1 and 95 DF,  p-value: 1.976e-08

Die Steigung (\(\hat{\beta}_1\)) ist die Mittelwertsdifferenz, welche oben (beim t-Test) schon berechnet wurde. Der Achsenabschnitt ist \(\hat{y}\) (bp_change), wenn \(x\) (drug) 0 ist. Wenn man kein Medikament nimmt (die Variable drug also den Wert 0 annimmt), ist die mittlere Blutdruckveränderung 0.551 (Aufgabe 4.). Dieser Wert ist identisch mit dem, den du mit dem t-Test berechnet hast.

Der t-Wert ist derselbe wie beim t-Test. Der t-Wert berechnet sich aus \(\hat{\beta}_1/se = -5.7594/0.9392\). Ein lineares Regressionsmodell mit einer dichotomen Prädiktorvariable liefert also die gleichen Resultate wie ein t-Test (abgesehen von der Welch-Korrektur). Eine exakte Übereinstimmung ergibt sich bei einem t-Test unter Annahme homogener Varianzen, auch wenn hier bei der angezeigten Anzahl an Dezimalstellen kein Unterschied feststellbar ist.

t.test(df.bp$bp_change~df.bp$drug,var.equal=TRUE)
## 
##  Two Sample t-test
## 
## data:  df.bp$bp_change by df.bp$drug
## t = 6.132, df = 95, p-value = 1.976e-08
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  3.894749 7.623959
## sample estimates:
##  mean in group no mean in group yes 
##         0.5510204        -5.2083333

Anmerkung: Durch die Welch-Korrektur kann es zu kleinen Unterschieden kommen. R macht automatisch eine Welch-Korrektur, falls die Varianzen in den beiden Gruppen nicht identisch sind. Mit var.equal = TRUE als zusätzlichem Argument in der Funktion t.test() kommen t.test() und lm() zum exakt gleichen Ergebnis.

Lineare Regression mit einer kontinuierlichen Prädiktorvariable

Aufgaben

  1. Erstelle ein Streudiagramm für die Variablen bp1 und bp2. Tipp: plot()
  2. Schätze folgendes lineares Regressionsmodell: \(Y_i = \beta_0 + \beta_1 x_i + \epsilon_i\), wobei \(Y = bp2\) und \(x = bp1\). Lass dir eine Zusammenfassung davon ausgeben.
  3. Interpretiere die Regressionskoeffizienten des geschätzten Modells.
  4. Wie lautet das 95%-Konfidenzintervall für \(\beta_1\)? Tipp: confint()
  5. Welchen durchschnittlichen Wert von bp2 schätzt das Modell für eine Person mit bp1 = 140?

Lösungen

Wir visualisieren die Analyse zuerst mit einem Streudiagramm:

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

Anhand der Grafik erwarten wir einen starken Zusammenhang zwischen den beiden Variablen. Wir passen nun das Modell an und schauen uns die Zusammenfassung an:

lm.bp1 <- lm(bp2 ~ bp1, data = df.bp)
summary(lm.bp1)
## 
## 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

Der Achsenabschnitt beträgt -10.043. Das ist der durchschnittlich erwartete Wert von bp2, wenn bp1 = 0 ist. Die inhaltliche Interpretation des Achsenabschnitts, einschliesslich seines Standardfehlers (SE), t-Werts und p-Werts, ist also wenig sinnvoll. Die Steigung beträgt 1.064. Das heisst: Wenn bp1 um 1 steigt, dann steigt bp2 durchschnittlich um 1.064 mmHg. Die Hypothese, dass diese Steigung 0 ist (\(H_0\)), erweist sich als wenig plausibel (p < 0.001).

Eine positive Steigung erscheint in diesem Fall (auf den ersten Blick) wenig sinnvoll. Eigentlich würden wir doch eine negative Steigung erwarten, weil der mittlere Blutdruck abnimmt. Man muss jedoch beachten, dass bei der Regressionsgleichung auch der Achsenabschnitt berücksichtigt werden muss. Wenn man das tut, erhält man für bp2 jeweils einen tieferen Wert als für bp1. Für bp1 = 140 schätzt das Modell z. B. einen Wert von \(-10.043 + 1.064 * 140 = 138.917\). 

Achtung: Dieses Modell sagt nichts über den Effekt des Medikaments aus!

Konfidenzintervalle für die Koeffizienten erhalten wir durch die Funktion confint():

confint(lm.bp1)
##                   2.5 %   97.5 %
## (Intercept) -22.3065973 2.220715
## bp1           0.9633251 1.163722

Wir vertrauen zu 95% darauf, dass sich die wahre Steigung im Intervall [0.963; 1.164] befindet. Das Konfidenzintervall ist konsistent mit dem sehr kleinen p-Wert.

Die geschätzten Parameter aus dem Output können nun in die Regressionsgleichung eingesetzt werden. Welchen durchschnittlichen Wert von bp2 schätzt das Modell für eine Person mit bp1 = 140?

\[ \hat{y_i} = \hat{\beta_0} + \hat{\beta_1} x_i = -10.043 + 1.064 * 140 = 138.917 \]