Herbst im Wallis, 2021

Einleitung

Im Statistikunterricht schlagen einem die Dozierenden traditionellerweise Tonnen von verschiedenen Tests um die Ohren. Im ersten Moment hat man oft das Gefühl, dass diese Tests alle ganz unterschiedlich funktionieren - schliesslich haben sie ja alle ganz andere Namen. Da kann es schnell passieren, dass man die Übersicht verliert (oder eben, vor lauter Bäumen den Wald nicht sieht).

Ziel dieser Übung ist es zu zeigen, dass 95% Vertrauensintervalle, der z-Test, der einstichproben t-Test und der zweistichproben t-Test gar nicht so unterschiedlich sind. Eigentlich wird immer das genau gleiche Prinzip angewendet.

Wir verwenden für diese Übung Daten von 10 Personen, an welchen ein Blutdruckmedikament getestet wurde. Bei den zehn Personen wurden folgende Effekte in Bezug auf den systolischen Blutdruck (in mm Hg) gemessen (positive Werte bedeuten eine Blutdruckerhöhung):

set.seed(1234)
df <- data.frame(BP = round(rnorm(10, -3, 10),3))
library(tidyverse)
df %>% 
  kable() %>% 
   kable_paper("hover", full_width = F)
BP
-15.071
-0.226
7.844
-26.457
1.291
2.061
-8.747
-8.466
-8.645
-11.900

Deskriptive Statistik

Aufgabe

Berechnen sie die folgenden Statistiken:

  • Mittelwert
  • Standardabweichung
  • n

Lösung

  • Mittelwert
attach(df)
m <- round(mean(BP), 3)
m
## [1] -6.83
  • Standardabweichung
s <- round(sd(BP), 3)
s
## [1] 9.96
  • n
n <- length(BP)
n
## [1] 10

Vertrauensintervalle

Aufgabe

Berechnen Sie die obere und untere Schranke des 95% Vertrauensintervalls (95% CI) für die Blutdruckveränderung.

Lösung

Um das 95% CI zu berechnen, brauchen wir zuerst den Standardfehler:

se <- round(s/sqrt(n), 3)

\(se = \frac{s}{\sqrt{n}} = 3.149\)

Ist die Populationsstandardabweichung bekannt, kann man auf der z-Tabelle den z-Wert nachschauen, welcher nur noch von 2.5% der Mittelwerte übertroffen wird. 2.5%, weil wir ein Intervall für die mittleren 95% der Mittelwerte berechnen. Die restlichen 5% teilen sich auf die beiden Seiten auf.

z2 <- round(qnorm(0.975),3)
z1 <- -z2

ggplot(NULL, aes(c(-3,3))) +
  geom_area(stat = "function", fun = dnorm, fill = "#00998a", xlim = c(-3, z1)) +
  geom_area(stat = "function", fun = dnorm, fill = "grey80", xlim = c(z1, z2)) +
  geom_area(stat = "function", fun = dnorm, fill = "#00998a", xlim = c(z2, 3)) +
  labs(x = "z", y = "") +
  scale_y_continuous(breaks = NULL) +
  scale_x_continuous(breaks = c(z1, z2))

ll <- round(m + -1.96 * se, 3)
ul <- round(m + 1.96 * se, 3)

Wir rechnen die z-Transformation rückwärts, um die beiden Schranken zu erhalten:

Untere Schranke: \(\bar{x} -1.96 * 3.149 = -13.004\)

Obere Schranke: \(\bar{x} + 1.96 * 3.149 = -0.66\)

Die zweite Situation ist jene, welche in der Praxis bedeutend häufiger auftritt: die Populationsstandardabweichung ist nicht bekannt. Alles, was sich verändert ist, dass man sich nun nicht an der z-Verteilung sondern an der t-Verteilung mit \(df = n - 1 = 9\) Freiheitsgeraden orientiert. Die Schwellenwerte werden wiederum unten angezeigt:

t2 <- round(qt(0.975,n-1), 3)
t1 <- -t2

ggplot(NULL, aes(c(-3,3))) +
  geom_area(stat = "function", fun = dnorm, fill = "#00998a", xlim = c(-3, t1)) +
  geom_area(stat = "function", fun = dnorm, fill = "grey80", xlim = c(t1, t2)) +
  geom_area(stat = "function", fun = dnorm, fill = "#00998a", xlim = c(t2, 3)) +
  labs(x = paste("t, df = ", n-1, sep = ""), y = "") +
  scale_y_continuous(breaks = NULL) +
  scale_x_continuous(breaks = c(t1, t2))

ll <- round(m + t1 * se, 3)
ul <- round(m + t2 * se, 3)

Wir rechen die “z-Transformation” rückwärts, um die beiden Schranken zu erhalen:

Untere Schranke: \(\bar{x} + -2.262 * 3.149 = -13.955\)

Obere Schranke: \(\bar{x} + 2.262 * 3.149 = 0.291\)

Interpretation

Wir vertrauen zu 95% darauf, dass der wahre Effekt irgendwo zwischen der unteren und der oberen Schranke ist. Wenn man die z-Verteilung braucht, schliesst das 95% CI die 0 (in diesem Fall unsere Null-Hypothese!) nicht mit ein. Wir wissen also schon jetzt, dass der p-Wert beim zweiseitigen z-Test kleiner sein muss als 5%. Beim CI mit der t-Verteilung ist es genau umgekehrt. Dieses schliesst die 0 mit ein. Somit ist beim zweiseitigen t-Test der p-Wert grösser als 5%.

Einstichproben Tests

Aufgabe

Führen Sie einen Test durch, um zu untersuchen, ob sich die mittlere Blutdrucksdifferenz statistisch signifikant von 0 unterscheidet. Denken Sie an die Formulierung der Hypothesen und überlegen Sie sich, ob ein zwei- oder einseitiger Test angebracht ist.

Lösung

In diesem Fall ist ein zweiseitiger Test angebracht. Anhand der Daten ist ersichtlich, dass es Personen gab, bei welchen der Blutdruck gestiegen ist. Es könnte daher sein, dass das Medikament sich negativ auf den Blutdruck auswirkt.

\(H_{0}: \mu_{0} = 0\)
\(H_{A}: \mu_{0} \ne 0\)

Genau wie bei den Vertrauensintervallen gibt es auch hier zwei Möglichkeiten: Beim z-Test (Populationsstandardabweichung bekannt) suchen wir den kritischen Wert auf der z-Tabelle, beim t-Test (Populationsstandardabweichung nicht bekannt) auf der t-Tabelle mit n-1 Freiheitsgeraden. Unsere Teststatistik (auch Prüfgrösse) berechnet sich jedoch genau gleich über die z-Transformation!

\(z = \frac{\mu_{0}-0}{se} = -2.17\)
\(t = \frac{\mu_{0}-0}{se} = -2.17\)

In der Praxis kennen wir wie oben erwähnt die Populationsstandardabweichung sehr selten, weshalb im Normalfall die t-Verteilung zum Zuge kommt. Die kritischen Werte haben wir bei den 95% CIs schon berechnet, diese sind für die z-Verteilung 1.96 und für die t-Verteilung 2.262. Wie oben angekündigt, befinden wir uns unter der z-Verteilung im Verwerfungsbereich. Wir würden also \(H_{0}\) zugunsten von \(H_{A}\) verwerfen. Bei der t-Verteilung mit 9 Freiheitsgeraden ist es, analog zum Vertrauensintervall, umgekehrt. Die Teststatistik fällt in den Nichtverwerfungsbereich. Die Null-Hypothese wird folglich beibehalten.

Mit Statistikprogrammen kann man den genauen, zweiseitigen p-Wert berechnen. Dieser beträgt bei der z-Verteilung 0.03 und bei der t-Verteilung 0.058

Anmerkung

Der t-Test für abhängige Daten ist genau das gleiche wie der einstichproben t-Test. Man berechnet im ersten Schritt die Differenz der Werte vor und nach einer Intervention für jede Person. Mit diesen Differenzen macht man dann einen einstichproben t-Test.

Zweistichproben t-Test

Gepaarte t-Tests (einstichproben t-Tests) werden in der Pflegeforschung häufig angewendet, wenn man Vorher- Nachhervergeliche macht (sogenannte pre-post designs). Ein grosses Problem solcher Studien ist, dass praktisch nie mit Sicherheit gesagt werden kann, ob der Effekt wirklich auf die Intervention zurückzuführen ist. Bei unserem Beispiel mit der Blutdruckbehandlung belibt beispielsweise unklar, wie sich der Blutdruck ohne die Intervetion verhalten hätte. Studien, welche bezüglich Kausalität deutlich besser eingestuft werden können, sind randomisierte, kontrollierte Studien (RCTs). Übertragen auf unser Beispiel, bekommt in einer solchen Studie eine Gruppe das aktive Medikament, während die andere Gruppe nur ein Placebo bekommt. Die beiden Stichproben sind unabhängig von einander. Die Daten sehen Sie unten:

set.seed(1234)
df <- data.frame(Diff_intervention = round(BP,3),
                 Diff_control = round(BP + rnorm(10, 4, 2),3))
df %>% 
  kable() %>% 
  kable_paper("hover", full_width = F)
Diff_intervention Diff_control
-15.071 -13.48
-0.226 4.33
7.844 14.01
-26.457 -27.15
1.291 6.15
2.061 7.07
-8.747 -5.90
-8.466 -5.56
-8.645 -5.77
-11.900 -9.68

Aufgabe 1

Berechnen Sie Mittelwert, Standardabweichung, n und den Standardfehler für jede Gruppe. Wie gross ist der Unterschied zwischen den beiden Gruppen?

Lösung

attach(df)
m1 <- round(mean(Diff_intervention), 3)
m2 <- round(mean(Diff_control), 3)
s1 <- round(sd(Diff_intervention), 3)
s2 <- round(sd(Diff_control), 3)
n1 <- length(Diff_intervention)
n2 <- length(Diff_control)
se1 <- round(s1/sqrt(n1), 3)
se2 <- round(s2/sqrt(n2), 3)
se_diff <- round(sqrt(se1^2 + se2^2), 3)

Interventionsgrupppe: \(\bar{x}\) = -6.832, \(s\) = 9.958, \(n\) = 10, \(se\) = 3.149

Kontrollgruppe: \(\bar{x}\) = -3.598, \(s\) = 11.949, \(n\) = 10, \(se\) = 3.779

Um den Unterschied der Gruppen numerisch auszudrücken, berechnen wir die Mittelwertsdifferenz:

\(\bar{x}_{diff} = \bar{x}_{inttervention} - \bar{x}_{control} = -3.234\)

Anmerkung: Beim einstichproben t-Test ging es nicht um die Mittelwertsdifferenz, sondern und die mittlere Differenz.

Aufgabe 2

Sie haben nun alle Zutaten, um den t-Test für unabhängige Stichproben durchzuführen. Denken Sie an die Hypothesen und überlegen Sie sich, ob ein ein- oder zweiseitiger Test sinnvoller ist.

Lösung

Aus den selben Gründen wie oben, ist ein zweiseitiger Test angebracht. Die Hypothesen könnten wie folgt aussehen:

\(H_{0}: \mu_{0A} = \mu_{0B}\)

\(H_{A}: \mu_{0A} \ne \mu_{0B}\)

Möglich wäre aber auch diese Formulierung:

\(H_{0}: \mu_{0A} - \mu_{0B} = 0\)

\(H_{A}: \mu_{0A} - \mu_{0B} \ne 0\)

Wie bei den vorangegangen Tests müssen wir unsere Teststatistik über die “z-Transformation” berechnen. Dafür brauchen wir den Standardfehler. Beim einstichproben t-Test brauchen wir den Standardfehler für die mittlere Differenz, beim zweistichproben t-Test brauchen wir den Standardfehler für die Mittelwertsdifferenz. Weil wir an dieser Stelle nicht zu sehr ins Detail gehen wollen, brauchen wir einfach diese Formel, um den Standardfehler für die Mittelwertsdifferenz (bei gleicher Gruppengrösse) zu berechnen (aufmerksame Studierende erinnern sich an den guten alten Pythagoras):

\(se_{diff} = \sqrt{se_{intervention}^2 + se_{control}^2}\)

Danach können wir die Teststatistik berechnen:

\(t = \frac{\bar{x}_{intervention} - \bar{x}_{control}}{se_{diff}} = -0.657\)

Bei der t-Verteilung mit \(n_{intervention} + n_{control} - 2\) Freiheitsgeraden schauen wir den kritischen Wert nach. Dieser beträgt für einen zweiseitigen Test 2.101.

t2 <- round(qt(0.975,n1+n2-2), 3)
t1 <- -t2

ggplot(NULL, aes(c(-3,3))) +
  geom_area(stat = "function", fun = dnorm, fill = "#00998a", xlim = c(-3, t1)) +
  geom_area(stat = "function", fun = dnorm, fill = "grey80", xlim = c(t1, t2)) +
  geom_area(stat = "function", fun = dnorm, fill = "#00998a", xlim = c(t2, 3)) +
  labs(x = paste("t, df = ", n1+n2-2, sep = ""), y = "") +
  scale_y_continuous(breaks = NULL) +
  scale_x_continuous(breaks = c(t1, t2))

Die Teststatistik fällt deutlich in den Nichtverwerfungsbereich (grau). Es besteht keine Evidenz gegen die Nullhypthese, weshalb diese beibehalten wird. Der genaue p-Wert beträgt 0.52 (nur mit Software zu berechnen).

So würde der Output einer Statistiksoftware aussehen (wir haben richtig gerechnet!):

t.test(Diff_intervention, Diff_control, paired = FALSE, var.equal = TRUE)
## 
##  Two Sample t-test
## 
## data:  Diff_intervention and Diff_control
## t = -0.7, df = 18, p-value = 0.5
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -13.6   7.1
## sample estimates:
## mean of x mean of y 
##     -6.83     -3.60

Aufgabe 3

Berechne das 95% Vertrauensintervall für die Mittelwertsdifferenz. Alle Zutaten sind bereits vorhanden ;-)

Lösung

Wir haben sowohl den kritischen Wert, sowie den Standardfehler für die Mittelwertsdifferenz berechnet. Wir müssen die Werte also nur noch in die Formel einsetzen:

Untere Schranke: \(\bar{x}_{diff} - t * se_{diff} = -3.234 - 2.101 * 4.919 = -13.569\)

Obere Schranke: \(\bar{x}_{diff} + t * se_{diff} = -3.234 + 2.101 * 4.919 = 7.101\)

Auch diese Werte stimmen mit jenen der Statistiksoftware überein:

t.test(Diff_intervention, Diff_control)
## 
##  Welch Two Sample t-test
## 
## data:  Diff_intervention and Diff_control
## t = -0.7, df = 17, p-value = 0.5
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -13.59   7.12
## sample estimates:
## mean of x mean of y 
##     -6.83     -3.60

Interpretation

Das Vertrauensintervall schliesst die 0 deutlich mit ein (darum der hohe p-Wert). Der wahre Effekt des Medikaments kann also sowohl positiv wie auch negativ sein. In der Forschung ist diese Situation nicht selten. Nämlich, dass bei pre-post Studien angebliche Effekte gefunden werden, welche dann in RCTs nicht bestätigt werden konnten.

Nicht-parametrischer zweistichproben Test

Der Wilcoxon-Test ist eine nicht-parametrische alternative zum einstichproben t-Test. Analog dazu gibt es auch eine nicht-parametrische Variante zum zweistichproben t-Test: der Mann-Whitney-U Test. Wir brauchen für die gleichen Daten wie oben, um diesen durchzuführen:

df %>% 
  kable() %>% 
  kable_paper("hover", full_width = F)
Diff_intervention Diff_control
-15.071 -13.48
-0.226 4.33
7.844 14.01
-26.457 -27.15
1.291 6.15
2.061 7.07
-8.747 -5.90
-8.466 -5.56
-8.645 -5.77
-11.900 -9.68

Aufgabe

Kommen Sie mit dem Mann-Whitney-U Test ebenfalls zum Schluss, dass die Null-Hypothese beibehalten wird?

Lösung

Analog zum zweistichproben t-Test wird ein zweiseitiger Test durchgeführt. Da es bei nicht-parametrischen Tests nicht um die Analyse von Mittelwerten geht, sehen die Hypothesen etwas anders aus:

\(H_{0}:P(X>Y)=P(Y>X)\) (Es besteht eine 50%-Wahrscheinlichkeit dafür, dass ein zufällig gezogener Wert aus X grösser ist als ein zufällig gezogener Mittelwert aus Y (und umgekehrt)

\(H_{A}:P(X>Y)≠P(Y>X)\) (Die Wahrscheinlichkeit ist nicht 50%, dass ein zufällig gezogener Wert aus X grösser ist als ein zufällig gezogener Mittelwert aus Y (und umgekehrt)

Als erstes ordnen wir die Scores und verteilen die Ränge:

df %>% 
  pivot_longer(1:2, names_to = "Gruppe", values_to = "Blutdruckveraenderung") %>% 
  mutate(Rang = rank(Blutdruckveraenderung)) %>% 
  pivot_wider(names_from = Gruppe, values_from = Rang) %>% 
  rename(Rang_intervention = Diff_intervention) %>% 
  rename(Rang_control = Diff_control) %>% 
  kable() %>% 
  kable_paper("hover", full_width = F)
Blutdruckveraenderung Rang_intervention Rang_control
-15.071 3 NA
-13.485 NA 4
-0.226 13 NA
4.329 NA 16
7.844 19 NA
14.013 NA 20
-26.457 2 NA
-27.148 NA 1
1.291 14 NA
6.149 NA 17
2.061 15 NA
7.073 NA 18
-8.747 7 NA
-5.896 NA 10
-8.466 9 NA
-5.559 NA 12
-8.645 8 NA
-5.774 NA 11
-11.900 5 NA
-9.680 NA 6

Die Rangsummen sind nun wie folgt:

df %>% 
  pivot_longer(1:2, names_to = "Gruppe", values_to = "Blutdruckveraenderung") %>% 
  mutate(Rang = rank(Blutdruckveraenderung)) %>% 
  group_by(Gruppe) %>% 
  summarise(Rangumme = sum(Rang)) %>% 
  kable() %>% 
  kable_paper("hover", full_width = F)
Gruppe Rangumme
Diff_control 115
Diff_intervention 95

Die Prüfgrösse wird wie folgt berechnet:

\(U_{1} = n_{1} * n_{2} + \frac{n_{1} * (n_{1} + 1)}{2} - T_{1}\) = 60

\(U_{2} = n_{1} * n_{2} + \frac{n_{2} * (n_{2} + 1)}{2} - T_{2}\) = 40

\(U = 40\)

Der kritische Wert für einen zweiseitigen Test mit \(\alpha\) = 0.05 beträgt 23 (siehe Tabelle). Da unsere Teststatistik deutlich grösser ist, besteht keine Evidenz gegen \(H_{0}\). Die Entscheidung ist also die gleiche, wie beim zweistichproben t-Test.

So könnte der Output einer Statistiksorftware aussehen:

wilcox.test(Diff_intervention, Diff_control, paired = FALSE, correct = FALSE)
## 
##  Wilcoxon rank sum exact test
## 
## data:  Diff_intervention and Diff_control
## W = 40, p-value = 0.5
## alternative hypothesis: true location shift is not equal to 0

Der p-Wert ist sehr ähnlich, wie beim zweistichproben t-Test.