Übungsblatt 2

Nicht zusammengesetzte Hypothesentests, Metaanalyse

  1. Das folgende R-Skript enthält eine Simulation von N = 1000 voneinander unabhängigen Hypothesentests mit nicht zusammengesetzten Hypothesen.

    library(pwr)

N <- 1000 # Anzahl Hypothesentests
rho <- 0.5 # Basisrate
alpha <- 0.05 # Signifikanzniveau der einzelnen Hypothesentests
delta <- 0.3 # Wahre Effektstärke, falls ein Effekt vorliegt
n <- 175 # Stichprobengröße pro Gruppe (t-Test, unabhängige Stichproben, ungerichtet)
power <- pwr.t.test(n = n, d = delta, sig.level = alpha)$power # Power der einzelnen Hypothesentests

set.seed(1) # Seed garantiert immer die gleichen Zufallszahlen für die Simulation 

# Simuliere N wahre Hypothesen mit Basisrate rho:
# 0: H0 ist war, 1: H1 ist wahr
H_Wahrheit <- sample(c(rep(0, time = N * rho), rep(1, times = N * (1 - rho))))

## Simulation ##
# Simuliere Testentscheidungen
H_Entscheidung <- rep(NA, N) # Erstelle leeren Vektor der Länge N
for (j in 1:N) {
  if (H_Wahrheit[j] == 0) {
    H_Entscheidung[j] <- rbinom(1, size = 1, prob = alpha)
    }
  if (H_Wahrheit[j] == 1) {
    H_Entscheidung[j] <- rbinom(1, size = 1, prob = power)
  }
  }
####

# Vergleiche wahre Hypothesen mit Testentscheidungen
Tabelle <- table(H_Wahrheit, H_Entscheidung)

    1. Bonus: Beschreiben Sie kurz, was im gegebenen R-Skript im Abschnitt ## Simulation ## passiert.

      Für jeden der N Hypothesentests wird mithilfe des „For Loops” eine Testentscheidung getroffen. Im ersten Durchgang des For Loops hat die Variable j den Wert 1, im zweiten Durchgang den Wert 2, …, im letzten Durchgang den Wert N, also for (j in 1:N). Falls für den Hypothesentest j die H0j wahr ist, also if (H_Wahrheit[j] == 0), wird eine Bernoulli-verteilte Zufallsvariable mit Parameter \(\pi = \alpha\) gezogen. Falls für den Hypothesentest j die H1j wahr ist, also if (H_Wahrheit[j] == 1), wird eine Bernoulli-verteilte Zufallsvariable mit Parameter \(\pi = 1 - \beta\) gezogen. Der Vektor H_Entscheidung enthält damit nach dem letzten Durchgang des For Loops an einer Stelle j eine 0, falls im Hypothesentest j eine Entscheidung für die Nullhypothese H0j getroffen wurde und eine 1, falls eine Entscheidung für die Alternativhypothese H1j getroffen wurde.

    2. Wie viele Fehler 1. Art werden in der Simulation beobachtet.

      Ergebnis der Simulation:

      Tabelle
                H_Entscheidung
H_Wahrheit   0   1
         0 474  26
         1 111 389

      Es werden 26 Fehler 1. Art beobachtet.

    3. Wie hoch ist in der Simulation der beobachtete Anteil der Fehler 1. Art an allen Tests, in denen die Nullhypothese gilt?

      Der beobachtete Anteil der Fehler 1. Art an allen Tests, in denen die Nullhypothese gilt beträgt:

      26 / (474 + 26) = 0.052

    4. Berechnen Sie für die vorgegebene Situation die Wahrscheinlichkeit für einen Fehler erster Art. Verwenden Sie bei der Berechnung ein Signifikanzniveau von \(\alpha = 0.05\) für die einzelnen Hypothesentests.

      Wir wollen die Wahrscheinlichkeit für einen Fehler erster Art berechnen, also die Wahrscheinlichkeit bei den N = 1000 voneinander unabhängigen Hypothesentests mindestens einen Fehler 1. Art zu begehen, wenn wir dabei wissen, dass in der Hälfte der Tests die Nullhypothese richtig ist. \[\alpha^{*} = 1 - (1 - \alpha)^{N} = 1 - (1 - 0.05)^{500} \approx 1\]

      Es kommt also bei 1000 Hypothesentests mit einem Signifikanzniveau von 0.05, falls für die Hälfte der Tests die Nullhypothese gilt, mit fast 100 prozentiger Sicherheit zu mindestens einem Fehler 1. Art.

    5. Wie hoch ist in der Simulation der beobachtete Anteil der falschen Entscheidungen für die Alternativhypothese an allen Entscheidungen für die Alternativhypothese.

      Der beobachtete Anteil der falschen Entscheidungen für die Alternativhypothese an allen Entscheidungen für die Alternativhypothese beträgt:

      26 / (26 + 389) = 0.063

    6. Berechnen Sie für die vorgegebene Situation die False Discovery Rate. Verwenden Sie bei der Berechnung die in der Simulation festgelegten Werte für \(\alpha,\rho,1 - \beta\) .

      \[FDR = \frac{\alpha \cdot \rho}{\alpha \cdot \rho + (1 - \beta) \cdot (1 - \rho)} = \frac{0.05 \cdot 0.5}{0.05 \cdot 0.5 + 0.8 \cdot (1 - 0.5)} = 0.059\]

      Obwohl die FWER sehr hoch ist, ist die FDR mit einem Wert von 0.059 immer noch vertretbar klein. Im Gegensatz zur FWER ist die FDR hier eine für die Praxis sehr sinnvolle Größe, da man sich für die einzelnen Tests interessiert. Ist die FDR klein, ist es plausibel, dass im Falle eines signifikanten Testergebnisses auch tatsächlich die Alternativhypothese wahr ist.

    7. Kontrollieren Sie mit der Bonferroni-Methode die Family-Wise Error Rate auf den Wert 0.05 und wiederholen Sie die Simulation (Hinweis: Achten Sie darauf, im Skript auch die Power neu zu berechnen). Wie hoch ist der beobachtete Anteil der Entscheidungen für die Alternativhypothese an allen Tests, in denen die Alternativhypothese gilt.

      Unter Verwendung des Bonferroni-korrigierten Signifikanzniveaus von

      alpha_new <- 0.05 / N
alpha_new # Signifikanzniveau der einzelnen Hypothesentests
      [1] 5e-05

      ergibt sich:

      ### Neue Simulation der Testentscheidungen: Parallel zum Code der Aufgabenstellung, nur mit korrigiertem alpha.

N <- 1000 # Anzahl Hypothesentests
rho <- 0.5 # Basisrate
alpha <- alpha_new
delta <- 0.3 # Wahre Effektstärke, falls ein Effekt vorliegt
n <- 175 # Stichprobengröße pro Gruppe (t-Test, unabhängige Stichproben, ungerichtet)
power <- pwr.t.test(n = n, d = delta, sig.level = alpha)$power # Power der einzelnen Hypothesentests

set.seed(1) # Seed garantiert immer die gleichen Zufallszahlen für die Simulation 
# Simuliere N wahre Hypothesen mit Basisrate rho:
# 0: H0 ist war, 1: H1 ist wahr
H_Wahrheit <- sample(c(rep(0, time = N * rho), rep(1, times = N * (1 - rho))))

## Simulation ##
# Simuliere Testentscheidungen
H_Entscheidung <- rep(NA, N) # Erstelle leeren Vektor der Länge N
for (j in 1:N) {
 if (H_Wahrheit[j] == 0) {
   H_Entscheidung[j] <- rbinom(1, size = 1, prob = alpha)
   }
 if (H_Wahrheit[j] == 1) {
  H_Entscheidung[j] <- rbinom(1, size = 1, prob = power)

  }}
####

# Vergleiche wahre Hypothesen mit Testentscheidungen
Tabelle <- table(H_Wahrheit, H_Entscheidung)
Tabelle
                H_Entscheidung
H_Wahrheit   0   1
         0 500   0
         1 445  55

      Der beobachtete Anteil der Entscheidungen für die H1 an allen Tests in denen die H1 gilt beträgt 55 / (445 + 55) = 0.110.

      Würde man die FWER kontrollieren tritt zwar mit hoher Sicherheit kein Fehler 1. Art mehr auf und auch die FDR ist nahe 0. Allerdings ist dann die Power so gering, dass fast keine der wahren Alternativhypothesen tatsächlich erkannt wird. Dies wäre in der Praxis keine sinnvolle Situation.

  2. Sie vermuten, dass Personen mit Narzissmus im Mittel eine höhere (stetige) Depressionsschwere aufweisen als Personen ohne Narzissmus. Dabei können Sie auf N = 3 bereits veröffentlichte Studien zurückgreifen. Die folgende Tabelle zeigt die Ergebnisse der drei Studien, wobei t die Realisation der Teststatistik (t-Test für unabhängige Stichproben mit gerichteten Hypothesen) bezeichnet, p den dazugehörigen p-Wert und n1 sowie n2 die jeweiligen Stichprobengrößen der einzelnen Studien.

    \(t\) \(p\) \(n_1\) \(n_2\)
    Studie 1 1.717 0.045 40 40
    Studie 2 3.200 0.001 20 20
    Studie 3 - 0.842 0.800 200 200

    1. Stellen Sie die statistischen Hypothesen für den Hypothesentest der Metaanalyse auf.

      \[H_0: \delta \leq 0\]

      \[H_1: \delta > 0\]

    2. Berechnen Sie für alle drei Studien die einzelnen Schätzwerte \(d_{j}\) für Cohen's \(\delta\).

      \[ \begin{aligned} d_{1} &= t_{1} \cdot \sqrt{\frac{n_{11} + n_{21}}{n_{11} \cdot n_{21}}} = 1.717 \cdot \sqrt{\frac{40 + 40}{40 \cdot 40}} = 0.383 \\ d_{2} &= 3.200 \cdot \sqrt{\frac{20 + 20}{20 \cdot 20}} = 1.012 \\ d_{3} &= - 0.842 \cdot \sqrt{\frac{200 + 200}{200 \cdot 200}} = - 0.084 \end{aligned} \]

    3. Berechnen Sie den metaanalytischen Schätzwert \(d\) für Cohen's \(\delta\).

      \[\begin{aligned} {\widehat{\sigma}}_{1\ wert}^{2} &= \frac{n_{11} + n_{21}}{n_{11} \cdot n_{21}} + \frac{d_{1}^{2}}{2 \cdot \left( n_{11} + n_{21} \right)} = \frac{40 + 40}{40 \cdot 40} + \frac{{0.383}^{2}}{2 \cdot (40 + 40)} = 0.051 \\ {\widehat{\sigma}}_{2\ wert}^{2} &= \frac{20 + 20}{20 \cdot 20} + \frac{{1.012}^{2}}{2 \cdot (20 + 20)} = 0.113\\ {\widehat{\sigma}}_{3\ wert}^{2} &= \frac{200 + 200}{200 \cdot 200} + \frac{{- 0.084}^{2}}{2 \cdot (200 + 200)} = 0.010 \\ \\ d &= \frac{1}{\sum_{j = 1}^{N}w_{j}}\sum_{j = 1}^{N}{w_{j} \cdot d_{j}} = \\ &= \frac{1}{\frac{1}{0.051} + \frac{1}{0.113} + \frac{1}{0.010}} \cdot \left( \frac{1}{0.051} \cdot 0.383 + \frac{1}{0.113} \cdot 1.012 + \frac{1}{0.010} \cdot - 0.084 \right) = \\ &= 0.063 \end{aligned}\]

    4. Geben Sie den kritischen Bereich für den metaanalytischen Hypothesentest bei einem Signifikanzniveau von 0.005 an.

      Die Teststatistik des metaanalytischen Hypothesentests ist unter der H0 approximativ standardnormalverteilt:

      qnorm(0.995, 0, 1)
      [1] 2.575829

      \[K_{T} = \lbrack + 2.576;\ + \infty\lbrack\]

    5. Treffen Sie die Testentscheidung und interpretieren Sie das Ergebnis.

      \[\begin{aligned} t &= \left( {\widehat{\theta}}_{wert} - \theta_{0} \right)\sqrt{\sum_{j = 1}^{N}w_{j}} = 0.063\sqrt{\frac{1}{0.051} + \frac{1}{0.113} + \frac{1}{0.010}} = 0.714 \\ t &= 0.714\ \notin \ K_{T} = \lbrack + 2.576;\ + \infty\rbrack \end{aligned}\]

      Entscheidung für die Nullhypothese, dass Personen mit Narzissmus im Mittel keine höhere Depressionsschwere aufweisen als Personen ohne Narzissmus.

  3. In fünf verschiedenen Studien wurde jeweils überprüft, ob sich Personen mit Angststörung und Personen mit Depression in ihrer mittleren negativen Selbstbewertung unterscheiden. Laden Sie die folgenden Datensätze herunter:
    , , , ,
    Verwenden Sie für alle Teilaufgaben ein Signifikanzniveau von 0.005.

    Datensätze einlesen 
    daten1 <- read.csv2('Studie1.csv')
daten2 <- read.csv2('Studie2.csv')
daten3 <- read.csv2('Studie3.csv')
daten4 <- read.csv2('Studie4.csv')
daten5 <- read.csv2('Studie5.csv')

    1. Berechnen Sie jeweils die Power der einzelnen Hypothesentests unter der Annahme, dass ein Effekt von \(\delta = 0.3\) vorliegt.

      library(pwr)
## Stichprobengrößen pro Gruppe
n1 <- nrow(daten1)
n2 <- nrow(daten2)
n3 <- nrow(daten3)
n4 <- nrow(daten4)
n5 <- nrow(daten5)

## Power einzelne Studien
p1 <- pwr.t.test(n = n1, d = 0.3, sig.level = 0.005)$power
p2 <- pwr.t.test(n = n2, d = 0.3, sig.level = 0.005)$power
p3 <- pwr.t.test(n = n3, d = 0.3, sig.level = 0.005)$power
p4 <- pwr.t.test(n = n4, d = 0.3, sig.level = 0.005)$power
p5 <- pwr.t.test(n = n5, d = 0.3, sig.level = 0.005)$power

p1
      [1] 0.06696447
      p2
      [1] 0.05617617
      p3
      [1] 0.04421589
      p4
      [1] 0.06256211
      p5
      [1] 0.427441

    2. Wie hoch ist unter dieser Annahme die Wahrscheinlichkeit dafür, dass alle einzelnen Hypothesentests signifikant werden?

      ## Wahrscheinlichkeit dafür, dass alle einzelnen Hypothesentests
## signifikant werden
p1 * p2 * p3 * p4 * p5
      [1] 4.447976e-06

    3. Berechnen Sie die Power des Hypothesentests unter der Annahme, dass ein Effekt von \(\delta = 0.3\) vorliegt, für den Fall, dass Sie die Daten aus allen fünf Studien zusammenfügen würden.

      ## Power zusammengefasste Studien
pwr.t.test(n = n1 + n2 + n3 + n4 + n5, d = 0.3, sig.level = 0.005)
      
     Two-sample t test power calculation 

              n = 297
              d = 0.3
      sig.level = 0.005
          power = 0.7985948
    alternative = two.sided

NOTE: n is number in *each* group

    4. Führen Sie die Hypothesentests in allen fünf Studien getrennt durch und treffen Sie jeweils eine Testentscheidung.

      ## t-Test Studie 1
t.test(daten1$depression, daten1$angst, var.equal = TRUE)
      
 Two Sample t-test

data:  daten1$depression and daten1$angst
t = 2.9677, df = 78, p-value = 0.003982
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.218807 1.110645
sample estimates:
  mean of x   mean of y 
 0.58613300 -0.07859292

      Der mithilfe von R ermittelte \(p\)-Wert beträgt \(p = 0.0040\). Da \(p = 0.0040 < \alpha = 0.005\) entscheiden wir uns für die Alternativhypothese.

      ## t-Test Studie 2
t.test(daten2$depression, daten2$angst, var.equal = TRUE)
      
 Two Sample t-test

data:  daten2$depression and daten2$angst
t = 4.4292, df = 68, p-value = 3.523e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.5132888 1.3550106
sample estimates:
mean of x mean of y 
1.0388052 0.1046555

      Der mithilfe von R ermittelte \(p\)-Wert beträgt \(p < 0.0001\). Da \(p < 0.0001 < \alpha = 0.005\) entscheiden wir uns für die Alternativhypothese.

      ## t-Test Studie 3
t.test(daten3$depression, daten3$angst, var.equal = TRUE)
      
 Two Sample t-test

data:  daten3$depression and daten3$angst
t = -0.17864, df = 56, p-value = 0.8589
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.6599655  0.5518950
sample estimates:
mean of x mean of y 
0.1319587 0.1859939

      Der mithilfe von R ermittelte \(p\)-Wert beträgt \(p = 0.8589\). Da \(p = 0.8589 > \alpha = 0.005\) entscheiden wir uns für die Nullhypothese.

      ## t-Test Studie 4
t.test(daten4$depression, daten4$angst, var.equal = TRUE)
      
 Two Sample t-test

data:  daten4$depression and daten4$angst
t = -0.79885, df = 74, p-value = 0.4269
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.6169705  0.2638369
sample estimates:
mean of x mean of y 
0.0563887 0.2329555

      Der mithilfe von R ermittelte \(p\)-Wert beträgt \(p = 0.4269\). Da \(p = 0.4269 > \alpha = 0.005\) entscheiden wir uns für die Nullhypothese

      ## t-Test Studie 5
t.test(daten5$depression, daten5$angst, var.equal = TRUE)
      
 Two Sample t-test

data:  daten5$depression and daten5$angst
t = -0.21165, df = 308, p-value = 0.8325
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.2572783  0.2073056
sample estimates:
  mean of x   mean of y 
-0.06593888 -0.04095250

      Der mithilfe von R ermittelte \(p\)-Wert beträgt \(p = 0.8325\). Da \(p = 0.8325 > \alpha = 0.005\) entscheiden wir uns für die Nullhypothese.

    5. Fügen Sie alle Datensätze mithilfe der Funktion rbind() zusammen, führen Sie in dieser Gesamtstichprobe den Hypothesentest durch und berechnen Sie mithilfe der Funktion cohen.d() aus dem effsize package den Punktschätzwert und ein 95%-Konfidenzintervall für \(\delta\).

      Hinweise 
      daten <- rbind(daten1, daten2, daten3, daten4, daten5)
cohen.d(Variable1, Variable2) # Befehl ist im Paket effsize
      ## t-Test und KI für delta zusammengefasste Studien
daten <- rbind(daten1, daten2, daten3, daten4, daten5)
t.test(daten$depression, daten$angst, var.equal = TRUE)
      
 Two Sample t-test

data:  daten$depression and daten$angst
t = 1.8415, df = 592, p-value = 0.06605
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.01055955  0.32796580
sample estimates:
mean of x mean of y 
0.1870456 0.0283425 
      library(effsize)
cohen.d(daten$depression, daten$angst)
      
Cohen's d

d estimate: 0.1511116 (negligible)
95 percent confidence interval:
      lower       upper 
-0.01028429  0.31250758

    6. Berechnen Sie die Punktschätzwerte für \(\delta\) getrennt für die einzelnen Studien, berechnen Sie auf der Basis dieser Schätzwerte ein kombiniertes metaanalytisches 95%-KI für \(\delta\) und führen Sie den kombinierten metaanalytischen Hypothesentest durch.

      ## Metaanalyse alle Studien

## Cohens d einzelne Studien
d1 <- cohen.d(daten1$depression, daten1$angst)$estimate
d2 <- cohen.d(daten2$depression, daten2$angst)$estimate
d3 <- cohen.d(daten3$depression, daten3$angst)$estimate
d4 <- cohen.d(daten4$depression, daten4$angst)$estimate
d5 <- cohen.d(daten5$depression, daten5$angst)$estimate
ds <- c(d1, d2, d3, d4, d5)

## Gewichte einzelne Studien
w1 <- 1 / (((2*n1)/(n1^2)) + ((d1^2) / (4*n1)))
w2 <- 1 / (((2*n2)/(n2^2)) + ((d2^2) / (4*n2)))
w3 <- 1 / (((2*n3)/(n3^2)) + ((d3^2) / (4*n3)))
w4 <- 1 / (((2*n4)/(n4^2)) + ((d4^2) / (4*n4)))
w5 <- 1 / (((2*n5)/(n5^2)) + ((d5^2) / (4*n5)))
ws <- c(w1, w2, w3, w4, w5)

## Kombinierter Schätzwert für delta
d <- sum(ds * ws) / sum(ws)
d
      [1] 0.1571459
      ## KI für delta
d - 1.96/sqrt(sum(ws))
      [1] -0.005501724
      d + 1.96/sqrt(sum(ws))
      [1] 0.3197936
      ## Hypothesentest für delta
t_statistik <- d * sqrt(sum(ws))
t_statistik
      [1] 1.893701
      2 * pnorm(-t_statistik)
      [1] 0.05826471

      Die metaanalytischen Methoden ergeben einen Schätzwert von \(\hat{\delta}_{Wert} = 0.16\) und ein Konfidenzintervall von \(KI_{\delta} = [-0.01, 0.32]\). Es gibt zwar leichte Abweichungen zu den aus allen Datensätzen zusammen berechneten Werten oben, jedoch sind diese vernachlässigbar. Genauso wie oben würden wir den metaanalytischen Hypothesentest bei einem \(\alpha = 0.05\) als nicht signifikant einschätzen.

      Für den (üblichen) Fall, dass nicht alle Daten in ihrer Rohform für eine Metaanalyse vorliegen, sind die hier verwendeten metaanalytischen Methoden eine gute Alternative um einen Effekt zu schätzen bzw. zu testen.

    7. Führen Sie die Berechnungen aus Teilaufgabe f nur auf Basis der Studien durch, in denen in Teilaufgabe d ein signifikantes Ergebnis herauskam.

      ## Metaanalyse nur mit Studien 1 und 2
ds12 <- c(d1, d2)
ws12 <- c(w1, w2)
d12 <- sum(ds12 * ws12) / sum(ws12)
d12 - 1.96/sqrt(sum(ws12))
      [1] 0.5057774
      d12 + 1.96/sqrt(sum(ws12))
      [1] 1.17505
      t12 <- d12 * sqrt(sum(ws12))
t12
      [1] 4.922392
      2 * pnorm(-t12)
      [1] 8.549259e-07

      Wir sehen also, dass in einem Fall der durch Publikationsbias auch tatsächlich auftreten könnte (nur signifikante Studien werden verwendet) die Metaanalyse zu einer fehlerhaften Einschätzung des Effekts kommen würde.

      Für die Generalisierbarkeit der Aussagen einer Metaanalyse ist es also entscheidend, welche Studien darin aufgenommen werden. Eine systematische Vorauswahl, wie beispielsweise nur signifikante Studien zu verwenden, für dazu, dass metaanalytische Schätzwerte und Hypothesentests im Sinne genau dieser Vorauswahl verzerrt werden.