Aufgabe 1
Lineare Beziehungen zwischen Variablen: Korrelationstest unter
Wir wollen uns ebenfalls die Power für den Korrelationstest ansehen. Dazu müssen wir allerdings korrelierte Variablen generieren. Um das hinzubekommen, müssen wir einige Eigenschaften der Normalverteilung ausnutzen: bspw. dass die Summe zweier normalverteilter Zufallsvariablen wieder normalverteilt ist. Für zwei unabhängige (unkorrelierte) standardnormalverteilte Zufallsvariablen
wieder standardnormalverteilt und um den Korrelationskoeffizienten
N <- 20
set.seed(12345)
X <- rnorm(N)
Z <- rnorm(N)
Y <- 0.5*X + sqrt(1 - 0.5^2)*Z
cor(X, Y) # empirische Korrelation
## [1] 0.579799
sd(X)
## [1] 0.8339354
sd(Y)
## [1] 1.232089
Falls Sie die oben genutzte Formel zur Generierung korrelierter Zufallsvariablen überprüfen wollen, dann setzen Sie doch einmal N = 10^6
(also eine Stichprobe von 1 Mio). Dann sollte die empirische Korrelation sehr nah an der theoretischen liegen. Auch sollten dann die empirischen Standardabweichungen sehr nah an der 1 liegen.
Verwenden Sie für diese Aufgabe stets den Seed 12345 (set.seed(12345)
).
- Betrachten Sie das Modell für eine Stichprobe von
N = 10^6
. Berichten Sie die empirische Korrelation sowie die empirischen Standardabweichungen.
Lösung
N <- 10^6
set.seed(12345)
X <- rnorm(N)
Z <- rnorm(N)
Y <- 0.5*X + sqrt(1 - 0.5^2)*Z
cor(X, Y) # empirische Korrelation
## [1] 0.4994574
sd(X)
## [1] 1.001315
sd(Y)
## [1] 0.9994427
Die Korrelation liegt bei
- Untersuchen Sie die Power des Korrelationstests für eine Korrelation von
und . Führen Sie eine Simulationsstudie durch. Wie groß ist die Power?
Lösung
N <- 20
set.seed(12345)
pcor_H1 <- replicate(n = 10000, expr = {X <- rnorm(N)
Z <- rnorm(N)
Y <- 0.5*X + sqrt(1 - 0.5^2)*Z
cortestH1 <- cor.test(X, Y)
cortestH1$p.value})
mean(pcor_H1 < 0.05) # empirische Power
## [1] 0.6385
Die Power des Korrelationstests für eine Korrelation von 0.5 für
- Stellen Sie die Verteilung der empirischen Korrelationen (für
und ) unter der dar.
Lösung
set.seed(12345)
cors_H1 <- replicate(n = 10000, expr = {X <- rnorm(N)
Z <- rnorm(N)
Y <- 0.5*X + sqrt(1 - 0.5^2)*Z
cor(X, Y)})
summary(cors_H1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.4293 0.3791 0.5080 0.4889 0.6177 0.8997
hist(cors_H1, breaks = 50)

Aufgabe 2
Lineare Beziehungen zwischen Variablen: Korrelationstest unter für ungleiche Varianzen
Wiederholen Sie die Analyse. Verändern Sie diesmal die Varianz der beiden Variablen X
und Y
. X
soll eine Varianz von 9 haben (multiplizieren Sie dazu X
mit 3, nachdem Sie Y
mithilfe von X
und Z
generiert haben), und Y
soll eine Varianz von 0.25 haben (multiplizieren Sie dazu Y
mit 0.5, nachdem Sie Y
mit Hilfe von X
und Z
generiert haben).
- Betrachten Sie das Modell für eine Stichprobe von
N = 10^6
. Berichten Sie die empirische Korrelation sowie die empirischen Standardabweichungen.
Lösung
N <- 10^6
set.seed(12345)
X <- rnorm(N)
Z <- rnorm(N)
Y <- 0.5*X + sqrt(1 - 0.5^2)*Z
X_new <- 3*X
Y_new <- 0.5*Y
cor(X_new, Y_new) # empirische Korrelation
## [1] 0.4994574
sd(X_new)
## [1] 3.003945
sd(Y_new)
## [1] 0.4997214
Die Korrelation liegt bei $\hat{\rho}{X\text{new}Y_\text{new}}=
- Führen Sie eine Simulationsstudie durch (für
und ). Wie verändert sich die Power des Tests durch die veränderten Varianzen?
Lösung
N <- 20
set.seed(12345)
pcor_H1_new <- replicate(n = 10000, expr = {X <- rnorm(N)
Z <- rnorm(N)
Y <- 0.5*X + sqrt(1 - 0.5^2)*Z
X_new <- 3*X
Y_new <- 0.5*Y
cortestH1 <- cor.test(X_new, Y_new)
cortestH1$p.value})
mean(pcor_H1_new < 0.05) # empirische Power
## [1] 0.6385
Die Power des Korrelationstests für eine Korrelation von 0.5 für
Aufgabe 3
Sensitivitätsanalyse für die Korrelation
In den Inhalten zu dieser Sitzung haben Sie neben der Poweranalyse auch die Sensitivitätsanalyse für den WebPower
auch für die Korrelation durchführen.
- Laden Sie zunächst das Paket
WebPower
und schauen Sie sich die Funktion für die Poweranalyse bei einer Korrelation an.
Lösung
library(WebPower)
## Warning: Paket 'WebPower' wurde unter R Version 4.3.1 erstellt
## Lade nötiges Paket: MASS
##
## Attache Paket: 'MASS'
## Das folgende Objekt ist maskiert 'package:dplyr':
##
## select
## Lade nötiges Paket: lme4
## Lade nötiges Paket: Matrix
## Warning: Paket 'Matrix' wurde unter R Version 4.3.2 erstellt
## Lade nötiges Paket: lavaan
## Warning: Paket 'lavaan' wurde unter R Version 4.3.1 erstellt
## This is lavaan 0.6-16
## lavaan is FREE software! Please report any bugs.
##
## Attache Paket: 'lavaan'
## Das folgende Objekt ist maskiert 'package:psych':
##
## cor2cov
## Lade nötiges Paket: parallel
## Lade nötiges Paket: PearsonDS
## Warning: Paket 'PearsonDS' wurde unter R Version 4.3.1 erstellt
?wp.correlation
- Nehmen Sie an, dass Sie eine Gruppe von
Personen untersucht haben. Sie möchten nun wissen, wie groß der Korrelationskoeffizient theoretisch sein müsste, damit eine Power von erreicht werden kann. Das -Fehleriveau soll dabei bei liegen.
Lösung
Für die Sensitivitätsanalyse legen wir in der Funktion wp.correlation
die Stichprobengröße (r
lassen wir hingegen leer:
wp.correlation(n = 50, r = NULL, power = 0.95, alpha = 0.05, alternative = c("two.sided"))
## Power for correlation
##
## n r alpha power
## 50 0.4780569 0.05 0.95
##
## URL: http://psychstat.org/correlation
Es müsste also eine Korrelation von 0.4780569 vorliegen, damit wir mit unserer Stichprobe eine Power von
Aufgabe 4
Type I-Error und Power zu einem -Niveau von des -Test
Wir wollen nun die Power des
Nutzen Sie den Seed 12345 (set.seed(12345)
).
- Führen Sie eine Simulation durch, um das empirische
-Niveau des -Tests zu bestimmen für . Vergleichen Sie das Ergebnis mit dem Ergebnis aus der Sitzung.
Lösung
N <- 20
set.seed(12345)
pt_H0 <- replicate(n = 10000, expr = {X <- rnorm(N)
Y <- rnorm(N)
ttestH1 <- t.test(X, Y, var.equal = TRUE)
ttestH1$p.value})
mean(pt_H0 < 0.001) # empirischer Alpha-Fehler
## [1] 0.0011
Der empirische
- Führen Sie eine Simulation durch, um die empirische Power des
-Tests zu bestimmen für , und . Vergleichen Sie das Ergebnis mit dem Ergebnis aus der Sitzung. Was bedeutet dies für die Wahl der Irrtumswahrscheinlichkeit?
Lösung
set.seed(12345)
pt_H1 <- replicate(n = 10000, expr = {X <- rnorm(N)
Y <- rnorm(N) + 0.5
ttestH1 <- t.test(X, Y, var.equal = TRUE)
ttestH1$p.value})
mean(pt_H1 < 0.001) # empirische Power
## [1] 0.0362
Die empirische Power liegt bei 3.62%. Dieser Wert fällt nun deutlich geringer aus, als die 33.5%, die wir in der Sitzung beobachtet hatten. Dies zeigt nochmal deutlich auf, dass wenn wir unsere Irrtumswahrscheinlichkeit drastisch reduzieren wollen, wir in Kauf nehmen, dass die Power einen Effekt zu finden, wenn dieser da ist, deutlich eingeschränkt wird!
Sie können sich die Power auch für andere Irrtumswahrscheinlichkeiten anschauen, indem Sie die 0.001
ersetzen durch Ihre gewünschte Irrtumswahrscheinlichkeit!
Aufgabe 5
Power-Plots für den -Test
Wir wollen nun die Power des
Nutzen Sie den Seed 12345 (set.seed(12345)
).
- Erstellen Sie einen Power-Plot für die folgenden Effekte
und bei einer Stichprobengröße von . Stellen Sie die Effektgröße auf der x-Achse dar.
Lösung
set.seed(12345)
pt_H1_0 <- replicate(n = 10000, expr = {X <- rnorm(20)
Y <- rnorm(20)
ttestH1 <- t.test(X, Y, var.equal = TRUE)
ttestH1$p.value})
pt_H1_0.25 <- replicate(n = 10000, expr = {X <- rnorm(20)
Y <- rnorm(20) + 0.25
ttestH1 <- t.test(X, Y, var.equal = TRUE)
ttestH1$p.value})
pt_H1_0.5 <- replicate(n = 10000, expr = {X <- rnorm(20)
Y <- rnorm(20) + 0.5
ttestH1 <- t.test(X, Y, var.equal = TRUE)
ttestH1$p.value})
pt_H1_0.75 <- replicate(n = 10000, expr = {X <- rnorm(20)
Y <- rnorm(20) + 0.75
ttestH1 <- t.test(X, Y, var.equal = TRUE)
ttestH1$p.value})
pt_H1_1 <- replicate(n = 10000, expr = {X <- rnorm(20)
Y <- rnorm(20) + 1
ttestH1 <- t.test(X, Y, var.equal = TRUE)
ttestH1$p.value})
pt_H1_1.25 <- replicate(n = 10000, expr = {X <- rnorm(20)
Y <- rnorm(20) + 1.25
ttestH1 <- t.test(X, Y, var.equal = TRUE)
ttestH1$p.value})
t_power_d <- c(mean(pt_H1_0 < 0.05),
mean(pt_H1_0.25 < 0.05),
mean(pt_H1_0.5 < 0.05),
mean(pt_H1_0.75 < 0.05),
mean(pt_H1_1 < 0.05),
mean(pt_H1_1.25 < 0.05))
Ds <- seq(0, 1.25, 0.25)
plot(x = Ds, y = t_power_d, type = "b", main = "Power vs. d")

Dem Plot ist zu entnehmen, dass die Power mit steigender Effektgröße ansteigt.
- Welcher Effekt muss mindestens bestehen, damit die Power bei
liegt?
Lösung
Diesem Plot ist nun zu entnehmen, dass eine Mittelwertsdifferenz von größer 0.8 nötig ist, damit die Power hinreichend groß ist. Außerdem wird in diesem Plot auch ersichtlich, dass wenn die Mittelwertsdifferenz 0 ist, dann sind wir gerade im Fall der
Aufgabe 6
Powervergleich: -Test vs. Wilcoxon-Test
Wir wollen nun die Power des
Nutzen Sie den Seed 12345 (set.seed(12345)
).
- Verwenden Sie das gleiche Setting wie aus der Sitzung und bestimmen Sie die Power des Wilcoxon-Tests für
, und . Vergleichen Sie das Ergebnis mit dem Ergebnis aus der Sitzung.
Lösung
N <- 20
set.seed(12345)
pt_H1_t <- replicate(n = 10000, expr = {X <- rnorm(N)
Y <- rnorm(N) + 0.5
ttestH1 <- t.test(X, Y, var.equal = TRUE)
ttestH1$p.value})
mean(pt_H1_t < 0.05) # empirische Power des t-Tests
## [1] 0.335
set.seed(12345)
pt_H1_W <- replicate(n = 10000, expr = {X <- rnorm(N)
Y <- rnorm(N) + 0.5
wilcoxonH1 <- wilcox.test(X, Y)
wilcoxonH1$p.value})
mean(pt_H1_W < 0.05) # empirische Power des Wilcoxon-Tests
## [1] 0.3198
Die empirische Power des
Bspw. mit solchen Fragen beschäftigen sich Methodiker:innen aus verschiedensten Disziplinen. Wenn Sie sich dafür interessieren, fragen Sie doch gerne in einer der beiden Abteilungen nach!