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 diese Überlegungen auf mehrfaktorielle Designs erweitern. Zunächst definieren wir für alle Simulationen:

n = 100 # Gruppengröße (100 Personen pro Gruppe)
k = 1000 # 1000 Durchläufe

Einfaktorielle Vergleiche

Einführend betrachten wir (noch einmal) eine einfaktorielle ANOVA mit vier experimentellen Bedingungen. Das bedeutet, wir untersuchen die Mittelwertsunterschiede von vier verschiedenen Gruppen (z. B. den vier Himmelsrichtungen) hinsichtlich eines Merkmales z. B. “Stimmung” als AV):

	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

In der Simulation erhalten wir \(\bar{\alpha}_{sim} = .24\) während der formale Wert \(\bar{\alpha}_{formal} = .26\) beträgt. Wichtig ist an dieser Stelle, dass die Diskrepanz zu \(\alpha = 0.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"), each = n),
                          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.4727 0.9056 0.2458 0.1973 0.0770 0.2627 0.4914 0.7543 0.0532 0.4600

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

mean(df_aov1x4 <= .05)

[1] 0.053

Wir sehen: \(\bar{\alpha}_{sim} = 0.053\) signifikante Ergebnisse. Die ANOVA vermeidet hier die starke \(\alpha\)-Inflation, die bei mehreren unabhängigen \(t\)-Tests auftreten würde.

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 15 \(t\)-Tests benötigen. Laut Formel für den kumulierten \(\alpha\)-Fehler würde das bedeuten:

\[{\displaystyle {\bar {\alpha }}=1-\left(1-\alpha _{}\right)^{15}} = 0.54\] 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 Temperatur (gilt für alle Zellen, hier ist nur eine exemplarisch herausgegriffen!) \[H_0: μ_{jw}~ - μ_{j} - μ_{w} + μ = 0\]

In 54 % der Fälle (einzelner Mittelwertsvergleiche) begehen wir einen α-Fehler! Wie sieht es in einer Simulation mit der ANOVA aus?

df_p_aov2x3 <- data.frame() # Leerer 
for(i in 1:k){
  data.test <- data.frame(
  Regen = rep(c("Ja", "Nein"), each = n * 3),
  Temperatur = rep(rep(c("warm", "mittel", "kalt"), each = n), 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", "Temperatur_p", "Regen:Temperatur_p")

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 Spalten und 1.000 zugehörigen \(p\)-Werten je Spalte:

# Wir zeigen nur die ersten 10 Zeilen
round(df_p_aov2x3[1:10, 1:3], digits = 4)

   Regen_p Temperatur_p Regen:Temperatur_p
1   0.9182       0.6673             0.4514
2   0.0529       0.6570             0.3788
3   0.6264       0.9334             0.3313
4   0.7552       0.2696             0.9895
5   0.5200       0.3276             0.0618
6   0.9597       0.6770             0.7674
7   0.6767       0.4071             0.4834
8   0.1322       0.9927             0.0902
9   0.3393       0.5130             0.2464
10  0.3422       0.0447             0.1759

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

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

[1] 0.142

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} = .54\). Auch bei mehrfaktoriellen ANOVA-Designs steigt die Wahrscheinlichkeit für einen \(\alpha\)-Fehler, je mehr Effekte (Haupteffekte und Interaktionen) gleichzeitig getestet werden.