Lernziele in R

Simulation zur \(\alpha\)-Fehler-Kumulierung

\(\alpha\)-Fehlerkumulierung in mehrfaktoriellen Designs

Im Beitrag alpha-Inflation (Teil 1) haben wir die Grundprinzipien zum kumulierten \(\alpha\)-Fehler simuliert. Hier wollen wir die Probleme für mehrfaktorielle Designs erweitern. Zunächst definieren wir für alle Simulationen:

n = 100 # Teilstichproben mit 100 Personen
k = 1000 # 1000 Durchläufe

Einfaktorielle Vergleiche

Einführend betrachten wir (noch mal) die einfaktorielle ANOVA an einem Beispiel eines einfaktoriellen Designs mit 4 experimentellen Bedingungen an. Das bedeutet, wir untersuchen die Mittelwertsunterschiede von vier verschiedenen Gruppen (z. B. den vier Himmelsrichtungen) hinsichtlich eines Merkmales (z. B. “Stimmung”):

	Norden	Osten	Süden	Westen
Stimmung	μ_A	μ_B	μ_C	μ_D

Wir vergleichen vier Gruppen jeweils miteinander, haben also 6 Paar-Vergleiche. Laut Formel erwarten wir für den kumulierten \(\alpha\)-Fehler:

\[{\displaystyle {\bar {\alpha }}=1-\left(1-\alpha _{}\right)^{6}} = 0.26\]

Berechnen wir zunächst diese 6 einzelnen \(t\)-Tests (alle unabhängig voneinander, was wir bereits im ersten Teil diskutiert haben):

set.seed(202302)
p_values_6t <- replicate(
  n = k, expr = {
    ps <- c(
      t.test(x = rnorm(n), y = rnorm(n))$p.value,
      t.test(x = rnorm(n), y = rnorm(n))$p.value,
      t.test(x = rnorm(n), y = rnorm(n))$p.value,
      t.test(x = rnorm(n), y = rnorm(n))$p.value,
      t.test(x = rnorm(n), y = rnorm(n))$p.value,
      t.test(x = rnorm(n), y = rnorm(n))$p.value
      )
    return(ps)
    }
  )
mean(
  apply(p_values_6t, MARGIN = 2, FUN = function(x){any(x <= .05)})
)

[1] 0.242

Wir haben nur eine minimale Abweichung und Simulation \(\bar{\alpha}_{sim} = .24\) und \(\bar{\alpha}_{formal} = .26\). Wichtig ist an dieser Stelle, dass die Diskrepanz zu \(\alpha = .05\) sehr deutlich ist.

Berechnen wir nun stattdessen eine ANOVA mittels aov:

##Leerer data.frame für die p-Werte erstellen
df_aov1x4 <- c()

for(i in 1:k){
  data.aov <- data.frame(Himmelsrichtung = rep(c("Norden", "Süden", "Osten", "Westen")),
                          Stimmung = rnorm(n * 4))
  model = aov(Stimmung ~ Himmelsrichtung, data = data.aov)
  df_aov1x4[i] = summary(model)[[1]][1,5] 
  }

Wir haben nun das Ergebnis von \(k = 1.000\) ANOVAs und 1.000 \(p\)-Werte basierend auf 4 Gruppen mit je \(n_{Norden} = n_{Süden} = n_{Osten} = n_{Westen} = 100\) Stimmungswerten.

# hier nur die ersten 10 Werte
round(df_aov1x4[1:10], digits = 4)

 [1] 0.2346 0.4128 0.7521 0.3394 0.1112 0.6336 0.6788 0.5230 0.9653 0.0504

Wie viele dieser 1.000 \(p\)-Werte sind signifikant?

mean(df_aov1x4 <= .05)

[1] 0.053

Wir sehen: \(\bar{\alpha}_{sim} = .053\) signifikante Ergebnisse. Die Berechnung durch eine ANOVA bringt also statistische Vorteile!

Zweifaktorielle Vergleiche

Was passiert, wenn wir Interaktionen in der ANOVA ins Spiel bringen? Betrachten wir eine ANOVA mit 2-Faktoren (Temperatur und Regen):

		Temperatur
Regen	warm	mittel	kalt
ja	μ_jw	μ_jm	μ_jk	μ_j
nein	μ_nw	μ_nm	μ_nk	μ_n
	μ_w	μ_m	μ_k	μ

Die Mittelwerte der jeweiligen Gruppen μ_jw bis μ_nk fassen wieder die Stimmung der einzelnen Gruppen zusammen. Für einfache Mittelwertsvergleiche, würden wir hier bereits 14 \(t\)-Tests benötigen. Laut Formel für den kumulierten α-Fehler würde das bedeuten:

\[{\displaystyle {\bar {\alpha }}=1-\left(1-\alpha _{}\right)^{14}} = 0.51\] Für die ANOVA ergeben sich in dem zweifaktoriellen Design 3 Signifikanztests.

Haupteffekt: Regen \[H_0: μ~_{j}~ = μ~_{n}\]
Haupteffekt: Temperatur \[H_0: μ~_{w}~ = μ~_{m} = μ~_{k}\]
Interaktion: Regen x Temepratur (gilt für alle Zellen, hier ist nur eine exemplarisch herausgegriffen!) \[H_0: μ_{jw}~ - μ_{j} - μ_{w} + μ = 0\]

In 50 % der Fälle begehen wir einen α-Fehler! Wie sieht es in einer Simulation mit der ANOVA aus?

df_p_aov2x3 <- data.frame() #Leeren dataframe für die p-Werte erstellen
for(i in 1:k){
  data.test <- data.frame(Regen = rep(c("Ja", "Nein"), each = 100),
                          Temperatur = rep(c("warm", "mittel", "kalt"), each = 100, times = 2),
                          Stimmung = rnorm(n*6))
  model = aov(Stimmung ~ Regen * Temperatur, data = data.test)
  output = summary(model)[[1]][1:3,5] 
  df_p_aov2x3 = rbind(df_p_aov2x3, output) 
}
names(df_p_aov2x3) = c("Regen_p", "Wetter_p", "Regen:Wetter_p")
df_p_aov2x3 <- t(df_p_aov2x3) #transformieren

Wir haben nun das Ergebnis von k = 3.000 ANOVAs (1000 Durchläufe für alle 3 Effekte). Das ergibt eine Tabelle df_p_aov2x3, mit 3 Zeilen und 1.000 zugehörigen \(p\)-Werten je Zeile:

round(df_p_aov2x3[, 1:9], digits = 4)

                 [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]
Regen_p        0.2318 0.4392 0.2546 0.9828 0.0480 0.7228 0.3619 0.0111 0.1544
Wetter_p       0.6673 0.6570 0.9334 0.2696 0.3276 0.6770 0.4071 0.9927 0.5130
Regen:Wetter_p 0.9181 0.0786 0.5629 0.9428 0.3533 0.8162 0.6714 0.7375 0.4302

Wie viele dieser 1.000 Analysen besitzen \(p\)-Werte \(\leq0.05\) ?

mean(apply(df_p_aov2x3, MARGIN = 2, FUN = function(x){any(x <= .05)}))

[1] 0.14

Offensichtlich taucht auch hier ein kumulierter \(\alpha\)-Fehler auf. Das hatten wir so nicht erwartet. Sollte nicht eigentlich der Vorteil darin liegen, den kumulierten \(\alpha\)-Fehler zu vermeiden? Zunächst einmal sind die simulierten \(\bar{\alpha}_{sim} = .14\) deutlich kleiner als die mit \(t\)-Tests zu erwartenden \(\bar{\alpha}_{formal} = .51\). Auch für die mehrfaktorielle ANOVA gilt, dass die Wahrscheinlichkeit des \(\alpha\)-Fehlers mit steigender Zahl an Tests steigt!