7 Macht und Stichprobengrösse

7.1 Theorie in Kürze

Vor einem Experiment sollte man zumindest überschlagsmässig ausrechnen, ob man den gesuchten Effekt mit den zur Verfügung stehenden Resourcen (d.h. wie viele Stichproben können wir maximal machen) überhaupt ermitteln kann. Bei Experimenten mit Tieren oder Menschen wird das von einer Ethikkommission in der Regel verlangt. Aber auch sonst empfiehlt sich dieses Vorgehen um Zeit- und Ressourcenverschwendungen zu vermeiden.
Nach einem Experiment verwendet man am besten ein Vertrauensintervall um die Ungenauigkeiten in der Schätzung zu kommunizieren.
Einige wenige Tests werden diesbzgl. in R sehr gut unterstützt. Allerdings kommen diese Funktionen schnell an ihre Grenzen.
Mit einer Simulation lässt sich immer die nötige Stichprobengrösse berechnen. Dazu braucht man zwei Zutaten: 1) Man muss Daten von einer konkreten Alternativhypothese simulieren können und 2) man muss diese Daten auswerten können (wie bei der geplanten Auswertung mit den “echten” Daten)
Man kann die Macht / Strichprobengrösse nur mit einer konkreten Alternative (diese kommt vom Anwender mit Kontextwissen) bestimmen. Oft ist das der intellektuell schwierigste Schritt. Wegen der Unsicherheit bei der konkreten Alternative kann die nötige Stichprobe meist nur sehr grob abgeschätzt werden.
Wichtige R Funktionen: , ,

7.2 Vorgefertigte Funktionen: 1-Stichproben t-Test

In R sind der t-Test (1- oder 2-Stichproben), die 1-weg ANVOA und der Proportion-Test mit entsprechenden Funktionen gut unterstützt.

Hier ein Beispiel mit dem 1-Stichproben t-Test. Angenommen, eine Fabrik füllt Wasserflaschen mit 1 Liter ab. Die Genauigkeit der Abfüllung (Standardabweichung) ist 1 ml. Die Abfüllmaschine wurde nun neu eingestellt und es soll geprüft werden, ob die Abfüllmenge hinreichend genau eingestellt ist. Dazu wollen wir \(n\) Flaschen entnehmen und mit einem zweiseitigen 1-Stichproben t-Test auf dem 5% Niveau testen, ob die Abfüllmenge wirklich 1 Liter sein könnte (d.h., das ist die Nullhypothese). Falls die Maschine im Mittel 1 ml zu viel abfüllt (=konkrete Alternative), soll der t-Test mit 80% Wahrscheinlichkeit (=Macht) verwerfen. Wie viele Flaschen müssen wir untersuchen ?

Wir verwenden die Funktion . Das Argument bezeichnet den Unterschied zw. Wert in der Nullhypothese und Wert in der konkreten Alternativhypothese (hier 1 ml). Es kommt bei diesem Test also nur auf den Unterschied zwischen den Hypothesen und nicht auf die absoluten Werte an (das ist z.B. beim Proportion Test nicht so, daher muss man dort sowohl als auch und nicht nur die Differenz angeben). Das Argument bezeichnet die Standardabweichung der Einzelmessungen (hier 1 ml). Die Einheiten sind egal, solange und in den gleichen Einheiten gemessen wurden (Probieren Sie es aus: , also in Einheiten ml, liefert das gleiche wie , also in Einheiten Liter).

power.t.test(delta=1, sd=1, sig.level = 0.05, power=0.8,
             type="one.sample", alternative = "two.sided")

## 
##      One-sample t test power calculation 
## 
##               n = 9.937864
##           delta = 1
##              sd = 1
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided

Wir runden \(n\) immer zur nächst grösseren ganzen Zahl auf. D.h., wir brauchen 10 Flaschen.

Wir können die Funktion auch noch auf andere Arten verwenden (siehe ):

Angenommen, wir können in dem Setting nur 5 Flaschen ziehen. Welche Macht hätten wir dann?

power.t.test(n=5, delta=1, sd=1, sig.level = 0.05,
             type="one.sample", alternative = "two.sided")

## 
##      One-sample t test power calculation 
## 
##               n = 5
##           delta = 1
##              sd = 1
##       sig.level = 0.05
##           power = 0.4013203
##     alternative = two.sided

Die Macht für eine Abweichung von 1 ml ist dann nur 0.4.

Angenommen, wir können 20 Flaschen ziehen und wünschen eine Macht von 80%. Welchen Unterschied zw. Null- und Alternativhypothese () können wir feststellen?

power.t.test(n=20, power = 0.8, sd = 1, sig.level = 0.05,
             type="one.sample", alternative = "two.sided")

## 
##      One-sample t test power calculation 
## 
##               n = 20
##           delta = 0.6604409
##              sd = 1
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided

In diesem Setting könnten wir ein von 0.66 mit einer Macht von 0.8 entdecken.

7.3 Simulation: 1-Stichproben t-Test

7.3.1 Ideale Welt

Die vorgefertigen Funktionen sind zwar sehr komfortabel, haben aber auch Nachteile: Zum einen gelten sie natürlich nur für die genannten Tests (t-Test, 1-weg ANOVA, Proportion Test) und nicht für andere Tests (Binomialtest, Fisher-Test, Lineare Regression, etc.). Zum anderen können sie nur Spezialfälle der genannten Tests behandeln (z.B. wird beim 2-Stichproben t-Test gleiche Gruppengrössen und gleiche Varianzen in beiden Gruppen verlangt).

Mit einer Simulation kann man immer die nötige Stichprobengrösse berechnen.

In einer Welt mit unbeschränkten Ressourcen könnten wir folgendermassen vorgehen, um die Macht für eine Alternative zu bestimmen:

Wiederhole 1000 mal:
Erhebe die Daten in einer Welt, in der die Alternativhypothese gilt (z.B. Stelle die Abfüllmaschine auf Abfüllmenge 1.001, fülle Flaschen ab und ziehe \(n\) Flaschen)
Mache den Test und schaue, ob die Nullhypothese verworfen wird
Wie oft wurde der Test unter den 1000 Wiederholungen verworfen ? Dieser Anteil ist ein Schätzwert für die gesuchte Macht

Natürlich haben wir beschränkte Ressourcen und können so ein Experiment nur einmal (und nicht 1000 mal) machen und zudem können wir die Welt in der Regel nicht in eine Welt gemäss konkreter Alternativhypothese einstellen (es ist z.B. nicht so einfach einen fairen Würfel auf die schnelle so zu fälschen, dass er die 6 mit Wahrscheinlichkeit 1/3 statt 1/6 zeigt.)

Häufig kann uns eine Computersimulation helfen. Vielleicht können wir die Daten ja mit dem Computer simulieren? Bei den Trinkflaschen ist das recht einfach. Wenn wir annehmen, dass die Abfüllmenge normalverteilt mit Erwartungswert \(1.001\) (also gemäss Alternativhypothese) und Standardabweichung \(0.001\) (in Einheiten Liter) ist, können wir mit dem Befehl entsprechende Zufallszahlen generieren:

rnorm(n=10, mean = 1.001, sd = 0.001)

##  [1] 1.001426 1.000705 1.001895 1.001878 1.001822 1.001689 1.001554 1.000938
##  [9] 1.000694 1.000620

Mit R können praktisch für alle erdenklichen Verteilungen Zufallszahlen erzeugt werden. Eine Liste der unterstützen Verteilungen finden Sie mit (und in spezialisierten Paketen gibt es noch mehr Auswahl). Beachten Sie, dass es zu jeder Verteilung vier Funktionen gibt: Immer eine mit den Anfangsbuchstaben (für die Wa.dichte, also pdf), (für die Verteilungsfunktion, also cdf), (für die Quantilsfunktion) und (zum Erzeugen von Zufallszahlen), also am Bsp der Normalverteilung , , und .

Oben haben wir die Macht berechnet, wenn wir 5 Trinkflaschen ziehen würden und dann einen zweiseitigen 1-Stichproben t-Test auf dem 5% machen würden. Das Ergebnis war eine Macht von \(0.4013\). Dieses Ergebnis ist exakt, d.h., wir können mit einer Simulation kein genaueres Ergebnis erzielen. Zur Illustration lösen wir das gleiche Problem aber trotzdem nochmals, diesmal aber mit einer Simulation. Das Ergebnis ist dann weniger genau, aber die Methode der Simulation lässt sich auch auf alle anderen Problemstellungen anwenden.

Zunächst legen wir uns auf eine konkrete Alternative fest. Wir hatten oben: , d.h., die Erwartungswerte unter \(H_0\) und \(H_A\) sind um eine Einheit verschoben. Wir könnten also z.B. konkret festlegen \(H_0:\mu = 1\) und \(H_A:\mu = 1.001\). Die Standardabweichung der Einzelbeobachtungen hatten wir als \(\sigma=0.001\) angenommen (zur Erinnerung: Wir müssen keine Einheiten festlegen, solange und , also \(\sigma\), in der gleichen Einheit angegeben werden). Zudem nehmen wir an, dass die Daten normalverteilt sind.

Jetzt müssen wir überlegen, wie wir 5 Beobachtungen gemäss \(H_A\) simulieren können. Das ist mit R sehr einfach:

dat <- rnorm(n=5, mean=1.001, sd=0.001)

Schliesslich müssen wir überlegen, wie wir mit diesen Daten den t-Test durchführen können und auf dem 5% Niveau entscheiden, ob verworfen wurde:

res <- t.test(x = dat, alternative = "two.sided", alpha = 0.05)
verworfen <- (res$p.value < 0.05)

Jetzt haben wir alle Bausteine zusammen. Die beiden Schritte “Daten gemäss \(H_A\) generieren” und “t-Test anwenden” müssen wir nun sehr oft wiederholen und am Schluss schauen, wie oft verworfen wurde. Dieser Anteil ist eine Schätzung der Macht. Die Genauigkeit der Schätzung können wir mit dem Wurzel-n-Gesetz abschätzen. (Wir verwenden , damit die Simulation exakt reproduziert werden kann.)

set.seed(123)
reps <- 1000
verworfen <- rep(NA, reps)
for (i in 1:reps) {
  dat <- rnorm(n=5, mean=1.001, sd=0.001)
  ## mu ist der Wert unter H_0
  res <- t.test(x = dat, alternative = "two.sided", alpha = 0.05, mu = 1) 
  verworfen[i] <- (res$p.value < 0.05)
}
mean(verworfen) ## geschätzte Macht

## [1] 0.403

sd(verworfen)/sqrt(reps) ## Std.fehler der geschätzten Macht

## [1] 0.01551876

Der geschätzte Wert \(0.403\) ist sehr nahe am wahren Wert (\(0.4013...\)) und der Std.fehler zeigt uns, dass dieses Ergebnis auf +- ein paar Prozent genau ist. Wenn wir es genauer haben wollen, müssen wir die Anzahl Schleifendurchläufe () erhöhen und die Genauigkeit wird gemäss dem Wurzel-n-Gesetz langsam (doppelte Genauigkeit, wenn man vervierfacht) besser.

Häufig möchte man nicht primär die Macht berechnen sondern eine Stichprobengrösse zu einer gegebenen Macht finden. Dazu verpackt man obigen Code einfach in eine Funktion mit Inputvariable und probiert verschiedene Werte von aus, bis man nahe genug an der gewünschten Macht ist.

myPower <- function(n) {
  reps <- 1000
  verworfen <- rep(NA, reps)
  for (i in 1:reps) {
    dat <- rnorm(n=n, mean=1.001, sd=0.001)
    res <- t.test(x = dat, alternative = "two.sided", alpha = 0.05, mu = 1)
    verworfen[i] <- (res$p.value < 0.05)
  }
list(m=mean(verworfen), s=sd(verworfen)/sqrt(reps))
}

Angenommen, wir suchen die nötige Stichprobengrösse um eine Macht von 80% zu erzielen (ich setze wieder einen random seed, damit das Ergebnis genau reproduziert werden kann:

set.seed(42)
myPower(n=5)$m ## zu klein

## [1] 0.403

myPower(n=15)$m ## zu gross

## [1] 0.952

myPower(n=10)$m ## ok

## [1] 0.804

Wie oben, kommen wir auf eine Stichprobengrösse von 10 Flaschen.