Übung zu Mittelwertsvergleichen und P-Werten

Zum Datensatz

In der Übung zu den Konfidenzintervallen hast du bereits den Datensatz weight-loss.csv kennengelernt. Hier geht es mit den gleichen Daten weiter. In dieser Übung stehen jedoch nicht die CI’s, sondern die P-Werte im Mittelpunkt. Du wirst sehen, dass 95% CI’s und P-Werte viel gemeinsam haben und auf dem selben Grundprinzip beruhen.

Aufgabe

Importiere den Datensatz weight-loss.csv in R und kontrolliere, ob alles richtig funktioniert hat.

Lösung

library(rio)
df.wloss <- import("../data/weight-loss.csv")
str(df.wloss)

## 'data.frame':    22 obs. of  4 variables:
##  $ ID   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ group: chr  "Control" "Control" "Control" "Control" ...
##  $ wl1  : int  4 4 4 3 5 6 6 5 5 3 ...
##  $ wl3  : int  3 3 1 1 2 4 4 1 1 2 ...

Die Variable group ist als character hinterlegt. Es ist eine gute Angewohnheit, solche Variablen in Faktoren umzuwandeln:

df.wloss$group <- factor(df.wloss$group)

Übung 1:

Vergleich eines Mittelwertes mit einem Referenzwert Wert (One Sample t-test):

In gewissen Situationen kann es sein, dass \(H_{0}\) einen festgelegten Wert annimmt (muss natürlich inhaltlich begründet sein). Überprüfe ob sich die mittlere Gewichtsreduktion nach einem Monat (wl1) statistisch signifikant von 0 Unterscheidet (unabhängig von der Gruppenzugehörigkeit). \(H_{0}\) wäre dabei:

\(H_{0}: \mu = 0\)

Aufgaben

Wie gross ist der p-Wert (zuerst manuelle Berechnung)?
Wie gross ist der p-Wert (mit Hilfe der t.test() Funktion)?
Wie ist dieser p-Wert zu interpretieren?
Wie hängen das 95% CI (wird beim t.test() angezeigt) und der p-Wert zusammen?

Lösungen

Wie gross ist der p-Wert (zuerst manuelle Berechnung)?

H0 <- 0
mean.wl1 <- mean(df.wloss$wl1)
sd.wl1 <- sd(df.wloss$wl1)
n.wl1 <- length(df.wloss$wl1)
se.wl1 <- sd.wl1/sqrt(n.wl1)
z <- mean.wl1/se.wl1
p <- 2* (1-pnorm(z))
p

## [1] 0

Anmerkungen:

Der p-Wert ist in Wirklichkeit nicht 0, sondern einfach sehr sehr klein…
Beachte, dass der p-Wert mit 2 multipliziert wird, da wir im Normalfall zweiseitige P-Werte berechnen.

Wie gross ist der p-Wert (mit Hilfe der t.test() Funktion)?

t.test(df.wloss$wl1, mu = 0) # entspricht dem One Sample t-test

## 
##  One Sample t-test
## 
## data:  df.wloss$wl1
## t = 13.323, df = 21, p-value = 1.034e-11
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  4.449685 6.095769
## sample estimates:
## mean of x 
##  5.272727

Hier sehen wir, dass der p-Wert nicht 0 sondern \(1.034 * 10^{-11}\) ist, also 0.00000000001034.
Auch hier wird der zweiseitige (two-tailed) P-Wert berechnet. Für einen einseitigen Test müsste man explizit alternative = "less" oder alternative = "greater" definieren, womit der P-Wert dann halb so gross wäre:

t.test(df.wloss$wl1, mu = 0, alternative = "greater")

## 
##  One Sample t-test
## 
## data:  df.wloss$wl1
## t = 13.323, df = 21, p-value = 5.171e-12
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  4.591715      Inf
## sample estimates:
## mean of x 
##  5.272727

Wie ist dieser P-Wert zu interpretieren?

Die Wahrscheinlichkeit, den beobachteten Wert (oder einen extremeren, also noch weiter weg von \(H_0\)) unter dem \(H_0\)-Modell zu erhalten, beträgt 1.034e-11, also 0.00000000001034.

Anmerkung: Auch hier erhälst du das 95% CI. Dieses entspricht exakt den Werten, welche du bei der Übung zu den Konfidenzintervallen berechnet hast. Die t.test() Funktion ist also eine weitere Möglichkeit, um CI’s zu berechnen.

Wie hängt das 95% CI (wird beim t-Test angezeigt) und der P-Wert zusammen?

Das 95% CI umschliesst den von \(H_{0}\) postulierten Wert 0 deutlich nicht. Man könnte auch sagen, das 95% CI schneidet \(H_{0}\) nicht. Wenn ein 95% CI \(H_{0}\) nicht schneidet, dann ist der P-Werte sicher kleiner als 5%. Je deutlicher das CI \(H_{0}\) nicht schneidet, desto kleiner ist der P-Wert.

Übung 2

In der Übung zu den Konfidenzintervallen hast du das 95% CI für die Mittelwertsdifferenz der beiden Gruppen zum ersten Zeitpunkt berechnet. Nun geht es zusätzlich um den P-Wert für diese Mittelwertsdifferenz (Two sample t-Test).

Aufgaben

Berechne den P-Wert für die Mittelwertsdifferenz zwischen den Gruppen zum Zeitpunk wl1 (manuelle Berechnung unter Verwendung der Z-Verteilung).
Berechne den P-Wert für die Mittelwertsdifferenz zwischen den Gruppen zum Zeitpunk wl1 (mit Hilfe der t.test() Funktion).
Interpretieren den Output des t-Tests.

Lösungen

Berechne den P-Wert für die Mittelwertsdifferenz zwischen den Gruppen zum Zeitpunk wl1 (manuelle Berechnung unter Verwendung der z-Verteilung).

mean.CG <- mean(df.wloss$wl1[df.wloss$group == "Control"])
sd.CG <- sd(df.wloss$wl1[df.wloss$group == "Control"])
n.CG <- length(df.wloss$wl1[df.wloss$group == "Control"])
se.CG <- sd.CG/sqrt(n.CG)

mean.IG <- mean(df.wloss$wl1[df.wloss$group == "DietEx"])
sd.IG <- sd(df.wloss$wl1[df.wloss$group == "DietEx"])
n.IG <- length(df.wloss$wl1[df.wloss$group == "DietEx"])
se.IG <- sd.IG/sqrt(n.IG)

mean.diff <- mean.CG - mean.IG
se.mean.diff <- sqrt(se.CG^2+se.IG^2) # Vorgehen bei homogenen Varianzen, siehe Anmerkung unten.

z <- mean.diff/se.mean.diff

p <- 2*pnorm(z)

p

## [1] 0.02688307

Anmerkung: Weil die Varianzen (und der Stichprobenumfang) der beiden Gruppen nicht gleich ist, wäre die Berechnung noch exakter, wenn zuerst eine gepoolte Standardabweichung berechnet werden würde. Wenn die Varianzen nicht homogen sind, sollte der Welch Test verwendet werden (Standardeinstellung in R).

Berechne den P-Wert für die Mittelwertsdifferenz zwischen den Gruppen zum Zeitpunk wl1 (mit Hilfe der t.test() Funktion).

t.test(wl1 ~ group, data = df.wloss)

## 
##  Welch Two Sample t-test
## 
## data:  wl1 by group
## t = -2.2132, df = 11.94, p-value = 0.04712
## alternative hypothesis: true difference in means between group Control and group DietEx is not equal to 0
## 95 percent confidence interval:
##  -3.37451537 -0.02548463
## sample estimates:
## mean in group Control  mean in group DietEx 
##                   4.5                   6.2

Wie erwartet ist der P-Wert mit der z-Verteilung kleiner als mit der t-Verteilung.

Interpretieren den Output des t-Tests.

Die Mittelwertsdifferenz beträgt 1.7kg [95% CI = 0.03 bis 3.37] zugunsten der Interventionsgruppe. Der P-Wert beträgt 0.047. Das heisst unter der Annahme, dass \(H_0\) korrekt ist, ist die Wahrscheinlichkeit, dass man eine solche oder eine noch extremere Differenz zwischen den Gruppen rein durch Zufall findet, 4.7%. Diese Differenz ist bei einem Signifikanzlevel von 5% statistisch signifikant.

Zusammenfassung

Weil wir den Populationsparameter \(\mu\) (oder \(\mu_1 - \mu_2\)) nicht kennen, nehmen wir die Stichprobenmittelwerte als besten Schätzer. Die Verteilung von \(\mu\) (oder \(\mu_1 - \mu_2\)) schätzen wir mit dem Standardfehler \(SE\). Via Teststatistik kann anhand der kumulativen Wahrscheinlichkeitsfunktion ein P-Wert berechnet werden.
Im Normalfall haben wir eine ungerichtete Alternativhypothese und somit testen wir zweiseitig. Dies hat zur Folge, dass der P-Wert verdoppelt wird (gilt bei der Normal- und t-Verteilung).
Der P-Wert ist die Wahrscheinlichkeit, den beobachteten Wert (oder einen extremeren) unter dem \(H_0\)-Modell zu erhalten.
Wenn ein einzelnes 95% CI \(H_0\) nicht überdeckt, ist der P-Wert kleiner als 5%.
Wenn zwei 95% CI’s zweier Gruppen sich nicht überlappen, ist der P-Wert kleiner als 5%. Achtung, das Umgekehrte ist nicht wahr: Die beiden CI’s können sich “überlappen”, trotzdem überlappt das CI für \(mu_1-mu2\) \(H_0\) nicht (da mehr Daten zur Verfügung).
Um CI’s zu berechnen, eignen sich die Funktionen CI(), group.CI() oder t.test().
P-Werte für Einstichproben-Tests oder Zweistichproben-Tests berechnen wir mit der t.test() Funktion.