Zum Datensatz

Für diese Übung brauchst du den Datensatz blood_pressure.csv. Den Datensatz findest du auf Moodle. Der Datensatz beinhaltet Informationen zu:

  • Geschlecht
  • Alter
  • Körpergrösse
  • Gewicht
  • ob die Person ein Medikament einnahm (drug)
  • Blutdruck zu zwei verschiedenen Messzeitpunkten (Prä/Post) in mmHg

Aufgaben

  1. Lade den Datensatz herunter und lese 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’s)? Tipp: is.na() oder summary()

Lösungen

Den Datensatz kannst du mit der import() Funktion aus dem Package rio importieren. Achte darauf, dass die Working Directory richtig eingestellt ist.

library(rio)
data <- import("~/Documents/Git/stats-msc-bfh/Data/blood_pressure.csv")
str(data)
## '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 is.na() Funktion 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 Cumputersprache 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(data))
## [1] 9

Man kann auch die summary() Funktion verwenden, dann werden die NAs pro Variable angezeigt. Achtung: bei character-Variablen funktioniert das nicht –> Ein weiterer Grund, warum man solche Variablen in Faktoren umwandeln sollte.

Ohne Umwandlung:

summary(data)
##       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:

data$gender <- factor(data$gender)
data$drug <- factor(data$drug)
summary(data)
##       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

Deskriptive Statistik

Aufgabe

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

  • Mittelwert
  • Median
  • Varianz
  • Standardabweichung
  • Quartile

Lösung

Beachte: der Datensatz hat fehlende Werte (NA’s) –> na.rm = TRUE ist darum nötig.

mean(data$bp1, na.rm = TRUE) 
## [1] 121.8283
median(data$bp1, na.rm = TRUE)
## [1] 121
var(data$bp1, na.rm = TRUE)
## [1] 119.9804
sd(data$bp1, na.rm = TRUE)
## [1] 10.95356
quantile(data$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

Damit es schöner aussieht, wird hier die x-Achse noch manuell programmiert. Es wird nicht erwartet, dass du das selbstständig tun kannst.

boxplot(data$bp1, data$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)

Berechung von Differenzen

Aufgaben

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

Lösungen

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

data$bp_diff <- data$bp1 - data$bp2
mean(data$bp_diff, na.rm = TRUE)
## [1] 2.298969

Eine positive Zahl heisst, dass bp1 grösser ist als bp2. Demnach hat der Blutdruck abgenommen. Es gibt keine allgemeine Regel, wie man Differenzen berechnet (bp2 - bp1 wäre auch möglich). Wichtig ist, dass die Interpretation korrekt ist.

t-Test

Aufgaben

  1. Untersuche, ob sich die mittlere Differenz (also der Mittelwert der Variable bp_diff) 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_diff zwischen Personen, welche ein Medikament eingenommen haben und denjenigen, welche kein Medikament zu sich nahmen? Stelle die Differenz pro Gruppe zuerst mittels zwei Boxplots dar.
  4. Überprüfe, ob sich die mittlere Differenz 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 One Sample t-test (oder paired t-test, was das Selbe ist) angebracht. Beachte, dass \(H_0: \mu_{bp.diff} = 0\).

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

Die mittlere Differenz (Abnahme) von 2.3 mmHg ist bei einem alpha-Level von 5% statistisch signifikant (p-value < 0.001). Unter \(H_0\) ist es also unwahrscheinlich, diese Differenz oder eine extremere alleine durch Zufall zu finden. Das 95% CI geht von 1.2 bis 3.4. Ob die mittlere Differenz 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 two sample 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(data$bp_diff ~ data$drug, xlab = "Medication")

Unten siehst du den Output des zweistichproben t-Tests. R macht automatisch eine kleine Korrektur (Welch), weil die Varianzen der beiden Gruppen nicht exakt gleich sind. Beachte, dass \(H_0: \mu_{yes} = \mu_{no}\).

t.test(data$bp_diff~data$drug)
## 
##  Welch Two Sample t-test
## 
## data:  data$bp_diff by data$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:
##  -7.624309 -3.894398
## 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. Die Mittelwertsdifferenz ist demnach -5.76 mmHg. Diese Mittelwertsdifferenz ist bei einem alpha-Level von 5% statistisch signifikant (p-value < 0.001). Sollte also \(H_0\) sitmmen, ist ein solches oder ein noch extremeres Ergebnis sehr unwahrscheinlich. Das angegebene 95% CI geht von -7.6 bis -3.9. Wir vertrauen zu 95% darauf, dass die wahre Blutdruckreduktion durch das Medikament zwischen 7.6 und 3.9 mmHg 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 mit Hilfe des t-Tests untersucht, ob die Differenz bp1 - bp2 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 gebraucht werden (was mathematisch äquivalent ist).

Aufgaben

  1. Rechne ein lineares Regressionsmodell mit bp_diff als abhängige Variable und drug als Prädiktorvariable: \(Y_i = \beta_0 + \beta_1 x_i + \epsilon_i\), wobei \(Y = bp\_diff\) und \(x = drug\). Lass dir eine Zusammenfassung dieses Modell 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 das Intercept im berechneten Modell?

Lösungen

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

lm.sex <- lm(bp_diff ~ drug, data = data)
summary(lm.sex)
## 
## Call:
## lm(formula = bp_diff ~ drug, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.449  -2.449   0.551   3.551   9.792 
## 
## 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 (b1) ist in die Mittelwertsdifferenz, welche oben (beim t-Test) schon berechnet wurde. Das Intercept ist der Wert von y (bp_diff), wenn x (drug) 0 ist. Wenn man kein Medikament nimmt (die Variable drug also den Wert 0 annimmt), ist die mittlere Blutdruckdifferenz -0.551. Dieser Wert ist identisch mit dem, den du mit dem t-Test berechnet hast.

Der t-Wert ist der gleiche, wie derjenige beim t-Test. Der t-Wert berechnet sich aus \(b_1/se\), also 5.7594/0.9392. Ein lineares Regressionsmodell mit einer dichotomen Prädiktorvariable liefert also exakt die gleichen Resultate wie ein t-Test (abgesehen von der Welch-Korrektur).

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ätzliches Argument in der t.test() Funktion wird t.test() und lm() zum exakt gleichen Ergebnis kommen.

Lineare Regression mit einer kontinuierlichen Prädiktorvariable

Aufgaben

  1. Erstelle einen Scatterplot 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. Was ist das 95% CI für \(b_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 Scatterplot:

plot(data$bp1, data$bp2)

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

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

Das Intercept beträgt -10.043. Das ist der durchschnittlich erwarterte Wert von bp2, wenn bp1 = 0 ist. Die inhaltliche Interpretation des Intercepts, inkl. dessen SE, t-Wert und p-Wert ist also wenig sinnvoll. Die Steigung beträgt 1.064. Das heisst wenn bp1 um 1 steigt, dann steigt bp2 durchschnitlich um 1.064 mmHg. Die Hypothese, dass diese Steigung 0 ist (\(H_0\)), erweist sich als wenig plausibel (p < 0.001).

Eine positive Steigung macht in diesem Fall (auf den ersten Blick) wenig Sinn. Eigentlich würde wir doch eine negative Steigung erwarten, weil der mittlere Blutdruck abnimmt. Man muss jedoch beachten, dass bei der Regressionsgleichung auch das Intercept 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 sie die wahre Steigung im Intervall [0.963; 1.164] befindet. Das Konfidenzintervall ist konsistent zum sehr kleinen P-Wert.

Die geschätzten Parameter aus dem Output können nun in die Regressionsgelichung 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 \]