library(tidyverse)
library(afex)
library(emmeans)
library(patchwork)
#for parallel backend
library(doParallel)
library(parallel)
library(foreach)11 Powerschätzungen für neue Experimente
11.1 Wozu Powerschätzungen?
Bei der Planung einer Studie ist es wichtig zu überlegen, wie viele Daten man benötigt, um die Schlussfolgerungen ziehen zu können, die man - gegeben die Hypothesen - gern ziehen möchte.
Das Problem ist: Daten enthalten immer auch Rauschen (noise, Fehler). Je größer das Rauschen, desto zufälliger schwanken die erhobenen Werte - vor allem, wenn man wenige Versuchspersonen und/oder wenige Trials erhebt. Das liegt daran, dass einerseits von Person zu Person (oder Trial zu Trial) hohe Variabilität besteht; andererseits werden Extremwerte nur unvollständig weggemittelt, wenn wenig Daten erhoben werden. Am Ende weiß man dann nicht: ist der beobachtete Effekt “echt” oder nur “noise”?
Dieses Problem geht man mit einer Poweranalyse an. Mit ihr plant man die Studie so, dass die erwarteten Effekte mit einer ausreichend hohen Wahrscheinlichkeit auch entdeckt werden, wenn es sie tatsächlich gibt. Umgekehrt kann man sich dann (ausreichend) sicher sein, dass ein Effekt wahrscheinlich wirklich nicht existiert, wenn die Studie keinen Effekt zeigt.
Wovon hängt die Power meines Experiments ab?
Die Power einer Untersuchung hängt von einer Vielzahl von Parametern ab.
Die wichtigsten sind:
| Größe der zu erwartenden Effekte | je größer der Effekt, desto besser (d.h. mit weniger Datenpunkten) lässt er sich entdecken, d.h. desto eher ist ein Unterschied in den Daten auch statistisch signifikant. | ||
| Anzahl Versuchspersonen | je mehr Versuchspersonen man erhebt, desto zuverlässiger wird die Messung. | |
| Reliabilität der Messung | Hohe Reliabilität bedeutet, dass man an zwei Messterminen einen sehr ähnlichen Wert erhält; dies ist der Fall, wenn das Messinstrument zuverlässig funktioniert. In typischen Experimenten der Allgemeinen und Biopsychologie spielt hierbei z.B. die Anzahl der Trials eine Rolle. Je mehr Trials man erhebt, desto zuverlässiger (= reliabler) wird der Mittelwert erfasst. |
|
| Das statistische Verfahren | verschiedene Verfahren haben unterschiedliche zugrunde liegende Prinzipien. Entsprechend sind verschiedene Verfahren ggf. auch unterschiedlich sensitiv. | |
11.1.1 Was will man bei einer Powerschätzung erreichen?
11.1.1.1 Powerschätzung: Anzahl der VP
Die üblichste Form der Powerschätzung bezieht sich auf die Anzahl der Versuchspersonen.
Hierbei wird eine Effektgröße angenommen (z.B. Cohens d), das statistische Verfahren festgelegt (z.B. T-Test für unabhängige Stichproben) und die gewünschte Wahrscheinlichkeit für einen Fehler zweiter Art (z.B. ein beta von 0.80) festgelegt.
11.1.1.2 Powerschätzung: Effektgröße
Alle diese Parameter sind wichtig für das Ergebnis der Poweranalyse. Der schwierigste Parameter ist häufig die Effektgröße. Sie sagt aus, wie groß der Unterschied zwischen den Experimentalbedingungen sein wird.
Dies ist natürlich schwierig zu sagen, denn man hat das Experiment ja noch nicht gemacht. Man muss die Effektgröße also schätzen. Wenn möglich, sollten Effektgrößen aus der Literatur geschätzt werden.
Hier gilt: je besser die Datenlage, desto verlässlicher die Schätzung. Am meisten kann man Werten aus einer Meta-Analyse trauen. Bei Einzelexperimenten können die Effekte über- oder unterschätzt sein. In beiden Fällen ist aber wichtig, dass die Effekte aus der Literatur von möglichst ähnlichen Experimenten/Paradigmen stammen wie das eigene, geplante Experiment. Die beste Metaanalyse hilft nichts, wenn man etwas ganz anderes plant.
11.2 Analytische vs. simulationsbasierte Powerschätzungen
Für einfache Designs kann die Anzahl an VP analytisch ermittelt werden. Bei komplexeren Designs bieten sich eher Simulationen an, da sie flexibler sind. Wegen der hohen Flexibilität befürworten wir den Einsatz von simulationsbasierten Verfahren.
Vorteile von analytischen und simulationsbasierten Powerschätzungen
11.2.1 Powerschätzung: G*Power
Für einfache Designs kann die Anzahl an VP analytisch ermittelt werden. Das bedeutet, dass es eine mathematische Formel gibt, mit der man die Power berechnen kann. Dazu gibt es allerlei Programme. Das bekannteste ist G*Power.
- G*Power und ähnliche Programme haben allerdings Einschränkungen. Für viele Designs gibt es derzeit keine analytischen Lösungen. Außerdem berechnen die meisten Programme die Power nur für die traditionellen Tests wie t-Test, Anova und Korrelation. Neuere statistische Verfahren wie Linear Mixed Models sind nicht implementiert.
11.2.2 Powerschätzung: Simulationen
Simulationen stellen eine deutlich flexiblere Alternative zu den analytischen Programmen dar.
Bei einer Simulation zieht man Zufallsstichproben, d.h. man tut so, als hätte man die Daten erhoben, in Wirklichkeit hat man sie aber nur erfunden. Dabei legt man der Datensimulation die Effektgröße zugrunde, die man auch bei der analytischen Powerschätzung verwendet hat.
Dadurch “erfindet” man Daten, die einerseits so sind, wie man es sich mit dem Experiment zu finden erhofft; andererseits sind die Daten verrauscht, da sie aus Zufallsverteilungen gezogen werden. Dieses Daten-Erfinden macht man viele Male. Jedes Mal berechnet man das statistische Verfahren und prüft, ob der Effekt signifikant geworden ist. Dieser Schritt ist der wichtigste: da wir die Daten ja so simulieren, dass wir wissen, dass der Effekt vorhanden ist, sollte das statistische Verfahren ihn eigentlich auch immer finden. Immer, wenn die simulierten Daten aber zu verrauscht waren, zeigt das statistische Verfahren fälschlicherweise keine Signifkanz an. Nun können wir die Simulation schrittweise verändern. Bspw. können wir verschiedene Anzahlen von Versuchspersonen simulieren.
Wir probieren also durch stures, massiertes Probieren auf der Basis von simulierten Datensätzen einfach aus, wie viele Versuchspersonen man testen müsste, um einen Effekt zu finden, der so groß ist, wie per Effektgröße in unserer Simulation angegeben. Je mehr Versuchspersonen wir simulieren, desto häufiger wird der statistische Test den Effekt auch entdecken. Wir verwenden dann so viele Versuchspersonen, dass die Anzahl der Nichtentdeckungen unserem gewünschten beta-Fehler (Fehler zweiter Art, oben im Beispiel war das 80%) entspricht.
Dabei muss man im Blick haben: das Ergebnis gilt nur, wenn der echte Effekt so groß ist wie in der Simulation angenommen. Um dem entgegenzutreten, kann man die Simulationen auch über eine ganze Reihe von Effektgrößen erstellen. Dann hat man einen Überblick, wie groß die Power ist, je nachdem wie groß der (unbekannte) Effekt in Wirklichkeit ist.
Vorteil der Simulation: Sehr hohen Flexibilität. Eine Simulation kann mit jeder Art von Design, Hypothese und statistischem Verfahren durchgeführt werden. Der Nachteil ist, dass hierfür extra an die geplante Studie angepasster Code geschrieben werden muss. Hier zeigen wir, eine Step-by-Step-Anleitung, die sich leicht an spezifische Studien anpassen lässt.
11.3 Packages und Daten
Zunächst benötigen wir wieder eine Reihe an Packages. Die meisten sind aus den anderen Kapiteln bekannt. Da man die Statistik viele Male rechnen muss, lohnt es sich, mehrere Prozessoren des Computers parallel einzusetzen. Dies erledigen die packages unter “parallel backend”.
11.4 Ziel, Design, Statistik und Hypothesen definieren
Nehmen wir an, wir wollen die in Kapitel 8 beschriebenen Interaktion von “posture” und “factor_B”, die mit der rm-Anova gezeigt wurde, in einem neuen Experiment replizieren. Wir wollen das gleiche grundsätzliche Design verwenden und also prüfen, ob es eine signifikante Interaktion aus posture und factor_B gibt. Diese Interaktion wurde schon in den Kapitel 8 gezeigt: sie kommt dadurch zustande, dass es bei “Level_2” aus factor_B einen größeren crossing-Effekt (Unterschied zwischen RTs bei crossed und bei uncrossed) gibt als bei “Level_1” aus factor_B.
Wir wollen uns bzgl. des Ergebnisses des Experiments sehr sicher sein und daher eine Power von 0.95 erreichen. Typische Power-Werte sind 0.8, 0.9 und 0.95; wie beim alpha-Fehler gibt es dabei aber kein richtig oder falsch. Bei jeder Studie muss man abwägen, wie gravierend der Fehler zweiter Art - der Effekt existiert, aber ich habe ihn nicht entdeckt - ist und die Power entsprechend wählen.
Wir wollen wissen, wie viele VP wir messen sollten.
11.5 Parameter aus existierenden Daten schätze und ähnliche Daten simulieren
Unser Ziel ist jetzt, Parameter aus Datensatzes zu berechnen, um dann möglichst ähnliche Daten simulieren zu können. Mit Parametern sind hier Mittelwerte, Varianzen und Interkorrelationen zwischen Bedingungen gemeint. So können wir für die Powerschätzung dann zum Beispiel mal nur 10 Personen ziehen und überprüfen, ob die statistischen Verfahren einen Effekt (von dem wir ja wisssen, dass er da ist) entdecken. Bevor wir das mit dem uns schon gut bekannten Datemsatz aus den letzten Kapiteln tun werden, tauchen wir aber kurz noch mal ein bisschen tiefer in Parameterberechnung und Simulation von Daten ab. Im Folgekapitel demonstrieren wir zwei Funktionen, die die Simulation von Daten sehr erleichtern.
11.5.1 Hilfsfunktionen für die Generierung von neuen Daten
Um neue Daten zu erzeigen, die möglichst eine ähnliche Struktur haben wie die empirisch beobachteten Daten, sind eine Reihe an Umformungen und Berechnungen nötig. Diese müssen nicht im Detail verstanden werden. Wir benutzen hier zwei Hilfsfunktionen:
- computeSampleParams berechnet die Parameter, die zur Schätzung nötig sind
- simulateNewData verewndet diese Parameter, um neue Daten zu erzeugen
11.5.1.1 computeSampleParams Funktion
Hier ein paar Infos zur Anwendung:
- Funktion erwartet Daten im langen Format
- Der Datensatz kann mehrere Bedingungen haben, muss aber so spezifiziert sein, dass es nur einen Wert pro Kombination gibt. Das heißt, wenn man mehrere Trials hat und Daten auch auf Trialbasis schätzen möchte, muss es eine Variable geben, die die Trials benennt
- Es müssen der Funktion die Bedeutungen von Variablen benannt werden:
- id.var: wie heißt die Variable, die den VP-Code enthält?
- condition.vars: alle Variablen, die Bedingungen enthalten (character-Vektor)
- value.var: welche Variable enthält (im langen Format) die Werte?
computeSampleParams <- function(data, id.var, condition.vars, value.var) {
#reduce data to the relevant variables
data <- data %>%
select(all_of(c(id.var, condition.vars, value.var)))
#get into wide format
data.wide <- data %>%
pivot_wider(names_from = all_of(condition.vars),
values_from = all_of(value.var), names_sep = "_")
#compute matrices
param.means <- apply(data.wide[, c(-1)], 2, FUN = function(x) mean(x, na.rm=T)) #mean
param.sds <- apply(data.wide[, c(-1)], 2, FUN = function(x) sd(x, na.rm=T)) #sds
param.cor <- cor(data.wide[, c(-1)], use = "pairwise.complete.obs") #correlations
param.cor <- psych::cor.smooth(param.cor) #fix corr mat
return(list(param.means=param.means,
param.sds=param.sds,
param.cor=param.cor))
}11.5.1.2 simulateNewData Funktion
Auch hier ein paar Infos zur Anwendung:
- die Funktion erwartet Matrizen mit Mittelwerten, Standardabweichungen und Korrelationen. Quasi genau, wie durch die obige Funktion auf der Basis von empirischen Daten berechnet. Ein Hinweis: man kann auch andere Werte einsetzen. Dazu mehr später!
- Es ist ganz wichtig, dass die Matrizen, die man eingibt alle nach einem Schema benannt sind, in dem in einem konsistenten Stil die Variablenausprägungen durhc einen Unterstrich getrennt sind.
- Der Funktion muss man auch sagen, die Variablen heißen, genau wie oben. Das ist in dem Fall nur dazu da, dass die Variablen dann “schön” benannt werden können.
- Die Funktion kann auch über das Argument increment.means die Mittelwerte anpassen. Das brauchen wir aber erst später
simulateNewData <- function(param.means, param.sds, param.cor, n.subjects, id.var, condition.vars, value.var, increment.means=1) {
#optional: increment the mean values
if (increment.means!=1) {
param.means <- mean(param.means) + (param.means - mean(param.means)) * increment.means
}
#compute covariance matrix
param.cov <- psych::cor2cov(rho = param.cor, sigma = param.sds)
#simulate the data
data.simulated.wide <- MASS::mvrnorm(n=n.subjects,
mu=param.means,
Sigma = param.cov)
data.simulated.wide <- as.data.frame(data.simulated.wide)
data.simulated.wide[, id.var] <- paste("id", sprintf("%03.f", 1:n.subjects), sep="")
#transform back to long format
data.simulated.long <- data.simulated.wide %>%
pivot_longer(cols = 1:(ncol(data.simulated.wide)-1),
values_to = {{value.var}}) %>%
separate(col = name, into = {{condition.vars}}, sep = "_")
return(data.simulated.long)
}11.5.1.3 Testen der beiden Hilfsfunktionen
Wer etwas genauer verstehen möchte, wie die Funktionen funktionieren, kann sich gerne den Kasten ansehen. Hier wird der “Nerd-Gain” ein bisschen hoch gedreht. Wer sich davon abgeschreckt fühlt, kann den Abschnitt aber auch einfach überspringen…
Demo der Funktionen an “Toy Data”
Hier wollen wir aus Demonstrationszwecken mal vollkommene Fantasiedaten erzeugen.
ToyData <- expand.grid(id=paste("id", 1:30, sep=""), posture=c("crossed", "uncrossed"), time=c("pre", "post"))
ToyData$RT <- rnorm(n = nrow(ToyData), mean=300, sd=25)
#einen Effekt verstecken: crossed post ist wie crossed pre, aber größer
ToyData <- ToyData %>%
pivot_wider(names_from = c(posture, time), values_from = RT) %>%
mutate(crossed_post = crossed_pre + rnorm(n=n(), mean=100, sd=15)) %>% #hier erhöhen wir crossed post
pivot_longer(cols = -id, names_to = c("posture", "time"), names_sep="_", values_to = "RT")
head(ToyData)# A tibble: 6 × 4
id posture time RT
<fct> <chr> <chr> <dbl>
1 id1 crossed pre 319.
2 id1 uncrossed pre 309.
3 id1 crossed post 423.
4 id1 uncrossed post 315.
5 id2 crossed pre 362.
6 id2 uncrossed pre 308.Jetzt werden mit der Funktion computeSampleParams die Parameter des Datensatzes berechnet.
ToyData.params <- computeSampleParams(data = ToyData, id.var = "id", condition.vars = c("posture", "time"), value.var = "RT")
ToyData.params$param.means
crossed_pre uncrossed_pre crossed_post uncrossed_post
307.6486 307.2642 408.0441 294.7963
$param.sds
crossed_pre uncrossed_pre crossed_post uncrossed_post
22.63433 23.14959 27.26097 25.03804
$param.cor
crossed_pre uncrossed_pre crossed_post uncrossed_post
crossed_pre 1.00000000 0.3628359 0.789688173 0.041264533
uncrossed_pre 0.36283592 1.0000000 0.319531241 -0.015670603
crossed_post 0.78968817 0.3195312 1.000000000 0.001949086
uncrossed_post 0.04126453 -0.0156706 0.001949086 1.000000000Hier wird es ein bisschen nerdy! Aber wir können an der obigen Matrix an Mittelwerten und Korrelationen sehen, dass wir ein bisschen an den Daten gespielt haben. Der Mittelwert von crossed post ist höher und crossed pre und post sind deutlich höher korreliert als die anderen Werte.
Nächster Schritt: jetzt generieren wir neue Daten ausgehend von den Parametern.
NewToyData <- simulateNewData(param.means = ToyData.params$param.means,
param.sds = ToyData.params$param.sds,
param.cor = ToyData.params$param.cor,
n.subjects = 500,
id.var = "id",
condition.vars = c("posture", "time"),
value.var = "RT")Jetzt können wir graphisch die neu erzeugten Werte mit den “wirklichen” vergleichen.
p.density.origToyData <- ToyData %>%
ggplot(aes(x=RT, fill=time)) +
geom_density(alpha=0.5) +
facet_wrap(~posture) +
ggtitle("original simulated data")
p.density.NewToyData <- NewToyData %>%
ggplot(aes(x=RT, fill=time)) +
geom_density(alpha=0.5) +
facet_wrap(~posture) +
ggtitle("reproduced simulated data")
p.density.origToyData + p.density.NewToyData
Sie können an den Abbildungen oben erkennen, dass sich die Daten recht ähnlich sehen. Das mag jetzt zunächst nicht so beeindruckend wirken, aber denken Sie darüber nach, dass ja die Reproduktion auf einer viel kleineren Anzahl an Informationen beruht, weil ja lediglich die Mittelwerte, Streuungen und Korrelationen bekannt waren.
11.6 Parameter aus dem bekannten Datensatz schätzen und neue Daten
So, jetzt sind wir gewappnet, um echte Daten zu vewenden.
Hier greifen wir auf den schon oft behandelten Datensatz zurück.
Achtung - Dieses Vorgehen hat den Nachteil, dass wir nicht wissen, wie groß der Messfehler bei der letzten Durchführung war. Besser wäre, man hätte mehrere Studien oder eine Metaanalyse; haben wir aber nicht.
Die Simulation erfolgt auf der Ebene der aggregierten Mittelwerte. Das heißt, dass wir nicht die einzelnen Trials simulieren, die eine Vp durchlaufen hat, sondern direkt die Mittelwerte jeder Bedingung.
load("RTL_beispieldaten_clean.RData")Zunächst werden wieder nur die Daten betrachtet, in denen die VPn richtig geantwortet haben. Wir aggregieren die Daten so, dass mittlere RTs für die Kombinationen aus participant, posture und factor_B übrig bleiben.
dcSim <- dc %>%
filter(respClean==1) %>% # nur korrekte Antworten verwenden
group_by(participant, factor_B, posture) %>% # gruppieren auf dem für uns relevanten Level
summarise(meanRT=mean(rtClean)) %>%
ungroup()`summarise()` has grouped output by 'participant', 'factor_B'. You can override
using the `.groups` argument.dcSim# A tibble: 72 × 4
participant factor_B posture meanRT
<fct> <fct> <fct> <dbl>
1 B006b-08 Level_1 uncr 638.
2 B006b-08 Level_1 crossed 687.
3 B006b-08 Level_2 uncr 611.
4 B006b-08 Level_2 crossed 736.
5 B006b-09 Level_1 uncr 565.
6 B006b-09 Level_1 crossed 634.
7 B006b-09 Level_2 uncr 581.
8 B006b-09 Level_2 crossed 648.
9 B006b-10 Level_1 uncr 772.
10 B006b-10 Level_1 crossed 1026.
# ℹ 62 more rowsHier noch ein Blick auf die Verteilung der Daten.
p.density.empirical <- dcSim %>%
ggplot(aes(x=meanRT, fill=factor_B)) +
geom_density(alpha=0.5) +
facet_wrap(~posture)
p.density.empirical
Unser Design ist ein repeated measurement, d.h. die Versuchspersonen haben mehrere Bedingungen durchlaufen. Der Vorteil von repeated measurements ist ja, dass die Daten innerhalb einer Versuchsperson über die Bedingungen hinweg ähnlicher sind als über die Versuchspersonen hinweg. In anderen Worten: jemand, die in einer Bedingung schnell ist, ist auch in anderen Bedingungen eher schnell; jemand, der in einer Bedingung langsam ist, ist auch in anderen Bedingungen eher langsam.
Anders ausgedrückt korrelieren die Werte innerhalb einer Versuchsperson. Diese Korrelation muss bei der Simulation berücksichtigt werden, denn natürlich wollen wir in unseren simulierten Daten den Effekt der Versuchsperson auch abbilden. Er ist zentral für die Berechnung der Signifikanz.
Dazu können wir direkt die Funktion computeSampleParams verwenden. Wir müssen allerdings aus “Level_1” und “Level_2” (Stufen von “factor_B”) noch den Unterstrich entfernen, damit die Funktion nicht verwirrt wird.
dcSim <- dcSim %>%
mutate(factor_B=gsub("_", "", factor_B))
dcSim.Params <- computeSampleParams(data = dcSim, id.var = "participant",
condition.vars = c("factor_B", "posture"),
value.var = "meanRT")
dcSim.Params$param.means
Level1_uncr Level1_crossed Level2_uncr Level2_crossed
729.1838 840.1628 716.1749 935.3435
$param.sds
Level1_uncr Level1_crossed Level2_uncr Level2_crossed
144.87466 170.39388 98.79463 180.87108
$param.cor
Level1_uncr Level1_crossed Level2_uncr Level2_crossed
Level1_uncr 1.0000000 0.7923623 0.8492734 0.5111338
Level1_crossed 0.7923623 1.0000000 0.6864440 0.8796122
Level2_uncr 0.8492734 0.6864440 1.0000000 0.5186879
Level2_crossed 0.5111338 0.8796122 0.5186879 1.0000000Analog zu oben sieht man auch hier wieder die Matrixen aus Mittelwerten, s.d. und Korrelationen. Sie sehen hier, dass alle Bedingungen korreliert sind. Das ist aus dem oben besprochenen Grund so, dass sich Daten innerhalb der gleichen Personen ähneln.
11.7 Neue Daten mit gleicher Parameterstruktur generieren
11.7.1 Parameter für die simulierten Daten
Zunächst legen wir ein paar Parameter fest, z.B. wie viele “Versuchspersonen” generiert werden sollen. Wir erzeugen einen großen Pool, aus dem wir dann später immer wieder ziehen werden.
set.seed(1111) #das macht das Beispiel reproduzierbar
nVPn = 1000 #Anzahl von simulierten Versuchspersonen11.7.2 Korrelierte Daten erzeugen
In diesem Schritt werden Daten erzeugt, die so gut es geht der oben gewünschten Struktur an Mittelwerten, Standardabweichungen und Interkorrelationen entsprechen. Auch dazu können wir auf die Funktion von oben, diesmal die Funktion simulateNewData zurückgreifen. Die Parameter haben wir ja bereits im letzten Abshcnitt berechnet und setzen wir hier direkt ein.
Der simulierte Datensatz sieht dann so aus:
dSimulated <- simulateNewData(param.means = dcSim.Params$param.means,
param.sds = dcSim.Params$param.sds,
param.cor = dcSim.Params$param.cor,
n.subjects = nVPn,
id.var = "participant",
condition.vars = c("factor_B", "posture"),
value.var = "meanRT")
head(dSimulated)# A tibble: 6 × 4
participant factor_B posture meanRT
<chr> <chr> <chr> <dbl>
1 id001 Level1 uncr 681.
2 id001 Level1 crossed 850.
3 id001 Level2 uncr 721.
4 id001 Level2 crossed 998.
5 id002 Level1 uncr 711.
6 id002 Level1 crossed 577.11.7.3 Unplausible Werte entfernen
Unmögliche Werte
Bei Simulationen kommt es manchmal vor, dass Werte generiert werden, die praktsich unmöglich sind (zum Beispiel sehr kleine oder sogar negative Reaktionszeiten). Solche Fälle sollte man entfernen.
Wir prüfen daher erstmal, ob es irgendwelche RTs gibt, die kleiner sind als 20 ms.
dSimulated %>%
summarise(sum(meanRT < 20)) %>%
as.numeric[1] 0Das ist nicht der Fall. Daher ist alles in Ordnung.
Man sollte sich der Gründlichkeit halber auch die simulierten Daten im Vergleich zu den originalen Daten anschauen und sie vergleichen. Wie das geht, sehen Sie in der Box.
Überprüfung der simulierten Daten
Um zu sehen, ob die Mittelwerte, Standardabweichungen und Interkorrelationen der simulierten Daten wirklich so sind wie bei den Ausgangsdaten, können wir die Parameter der simulierten Daten berechnen.
Wir wenden also wieder genau die gleichen Funktionen wie oben an.
dSimulated.Params <- computeSampleParams(data = dSimulated, id.var = "participant",
condition.vars = c("factor_B", "posture"),
value.var = "meanRT")Zum Vergleich “echte Daten vs. simulierte Daten” plotten wir beides zusammen.
Params.empirical.simulated <- rbind(
data.frame(param="mean", empirical=as.numeric(dcSim.Params$param.means),
simulated=as.numeric(dSimulated.Params$param.means)),
data.frame(param="sd", empirical=as.numeric(dcSim.Params$param.sds),
simulated=as.numeric(dSimulated.Params$param.sds)),
data.frame(param="cor", empirical=as.numeric(dcSim.Params$param.cor[upper.tri(dcSim.Params$param.cor)]),
simulated=as.numeric(dSimulated.Params$param.cor[upper.tri(dSimulated.Params$param.cor)])))
Params.empirical.simulated %>%
ggplot(aes(x=empirical, y=simulated)) +
geom_point() +
geom_abline(slope=1) +
facet_wrap(~param, scales="free")
Es ist zu sehen, dass die Werte sehr ähnlich sind. Die simulierten Daten haben also sehr ähnliche Parameter wie die empirischen Daten.
Nun werden wir noch die Verteilungen als Density-Plots ansehen.
p.density.simulated <- dSimulated %>%
mutate(posture=factor(posture, levels=c("uncr", "crossed")),
factor_B=factor(factor_B, levels=c("Level1", "Level2"), labels=c("Level_1", "Level_2"))) %>%
ggplot(aes(x=meanRT, fill=factor_B)) +
geom_density(alpha=0.5) +
facet_wrap(~posture)
p.density.empirical + ggtitle("empirical") + xlim(c(0, 1600)) +
p.density.simulated + ggtitle("simulated") + xlim(c(0, 1600)) +
plot_layout(ncol=1, guides = "collect") 
Der Vergleich der Abbildung verdeutlicht, dass die genrierten Daten im Großen und Ganzen ganz gut den empirischen Daten entsprichen. Man beobachtet aber auch Unterschiede. Die empirischen Daten haben eine gewisse Tendenz zu einer Bimodalität (zwei Hügel), was die simulierten Daten nicht abbilden.
Warum ist das so?
- Für die Simulation der Daten müssen wir festlegen, aus welcher Grundgesamtheit die künstlichen Daten gezogen werden sollen. Dazu haben wir angenommen, dass die Daten normalverteilt sind. Dass wir dies nun in den erzeugten Daten auch sehen, ist also per Definition so.
Das ist etwas, was man im Hinterkopf behalten sollte. Daten so zu simulieren, wie sie im echten Leben sind, ist ein sehr komplexes Problem. Je ungewöhlicher die Verteilung von empirischen Daten ist, desto schwerer wird es, sie mit den vereinfachten Annahmen einer normalverteilten Population zu simulieren. Unsere Simulation ist nichts anderes als ein Modell des Prozesses, der die experimentellen Daten erzeugt hat.
- Jedes Modell ist nur eine Annäherung an die Wirklichkeit. Ob wir mit der Ähnlichkeit der echten und der erzeugten Daten zufrieden sind, müssen wir selbst entscheiden. Falls nein, müssen wir darüber nachdenken, ob wir eine andere Verteilung zur Erzeugung der erfundenen Daten verwenden können/sollten.
Wie ist das nun in unserem Fall hier? Wir sollten berücksichtigen, dass die empirischen Daten auf nur wenigen Versuchspersonen beruhen.
Die Bimodalität, die wir im echten Datensatz sehen, könnte also ein Zufallsergebnis sein.
Tatsächlich: der Hügel am oberen Ende bei Level_1 uncr beruht zum Beispiel nur auf einem einzigen Wert. Hätten wir mehr Daten erhoben, würden wir möglicherweise eine Normalverteilung beobachtet haben.
Das macht noch einmal deutlich, dass eine Simulation, die sich auf vorhandene Daten beruft, sehr von der Qualität dieser Daten abhängt, die man als Ausgangspunkt nimmt. Mehr Daten führen zu einer akkurateren Schätzung.
11.8 Statistisches Modell definieren
Als nächstes definieren wir das statistische Modell, das wir verwenden, um die Hypothesen zu prüfen. Wir verwenden zunächst eine Versuchspersonenanzahl, die wir selbst definieren. Mit dem Code kann man also schon einmal “rumspielen” und schauen, wie sich die Ergebnisse verändern, wenn man den “nVP”-Parameter verändert.
nVP <- 50 # Anzahl der aus dem Datensatz gezogenen VersuchspersonenJetzt werden aus unserem generierten Datensatz mit 1000 Vpn zufällig 50 Personen gezogen und das statistische Modell gerechnet.
subjects <- unique(dSimulated$participant) #Vektor, der alle VP-Labels enthält
subjects.now <- sample(subjects, size = nVP, replace = F) #Zufallsziehung aus dem Versuchspersonen-Vektor; replace = F bedeutet, dass jede Vp nur einmal gezogen werden darf - ohne Zurücklegen also
# Zufällig gezogene Daten (Auschnitt mit den oben ausgewählten "Personen")
data.powerEst.now <- dSimulated %>%
filter(participant %in% subjects.now) %>%
mutate(posture=factor(posture, levels=c("crossed", "uncr")))
# statistisches Modell definieren: repeated-measures Anova
simRTs.rmanova <- aov_ez(id = "participant", dv = "meanRT", within = c("factor_B", "posture"), data = data.powerEst.now)
# Modell Ergebnisse ausgeben
simRTs.rmanovaAnova Table (Type 3 tests)
Response: meanRT
Effect df MSE F ges p.value
1 factor_B 1, 49 4933.93 9.21 ** .011 .004
2 posture 1, 49 14738.11 104.13 *** .267 <.001
3 factor_B:posture 1, 49 1751.01 88.21 *** .035 <.001
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1Man sieht nun an der Anova, dass die factor_B:posture-Interaktion signifikant ist.
Wir plotten die Ergebnisse erneut, um sie besser nachvollziehen zu können.
as.data.frame(emmeans(simRTs.rmanova, ~factor_B:posture)) %>%
ggplot(aes(x=factor_B, y=emmean, color=posture)) +
geom_point(position = position_dodge(width=0.5), size=3) +
geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL),
position = position_dodge(width=0.5), width=0.2) +
ylab(label = "reaction time [ms]")
11.9 Die Hypothese(n) definieren
Als nächstes formalisieren wir die kritischen Hypothesen (die man im Zusammenhang mit der Abbildung oben gut nachvollziehen kann). Wir sehen die Annahmen als bestätigt, wenn drei Dinge zutreffen:
- Es gibt eine signifikante Interaktion aus posture und factor_B
- Bei posture=crossed ist die RT bei factor_B=Level_2 signifikant höher als factor_B=Level_1
- Die Differenz aus crossed-uncrossed ist bei factor_B=Level_2 signifikant größer als bei factor_B=Level_1
Wir werden diese drei Hypothesen nun mit dem emmeans-Paket als post-hoc-Tests prüfen.
Der erste Test ist einfach, denn man kann den p-Wert einfach dem Anova-Table entnehmen und prüfen, ob er kleiner als 0.05 ist.
p_val_Hyp1 <- as.data.frame(simRTs.rmanova$anova_table) %>%
rownames_to_column() %>%
filter(rowname=="factor_B:posture") %>%
select("Pr(>F)") %>%
as.numeric()
Hyp1_confirmed <- c(p_val_Hyp1 < 0.05)
Hyp1_confirmed[1] TRUEDie Hypothese 1 ist also bestätigt.
Für den zweiten Test müssen wir zwei Dinge betrachten: erstens muss der p-Wert des post-hoc Vergleichs < 0.05 sein; zweitens muss die RT-Differenz Level_1 - Level_2 negativ sein.
simEmm.factorBXposture <- emmeans(simRTs.rmanova, ~factor_B|posture)
simEmm.factorBXposture.pairwise <- pairs(simEmm.factorBXposture, by = NULL, adjust="none")
diff_val_Hyp2 <- as.data.frame(simEmm.factorBXposture.pairwise) %>%
filter(contrast=="Level1 crossed - Level2 crossed") %>%
select("estimate") %>%
as.numeric()
p_val_Hyp2 <- as.data.frame(simEmm.factorBXposture.pairwise) %>%
filter(contrast=="Level1 crossed - Level2 crossed") %>%
select("p.value") %>%
as.numeric()
Hyp2_confirmed <- (diff_val_Hyp2 < 0) & (p_val_Hyp2 < 0.05)
Hyp2_confirmed[1] TRUEHypothese 2 ist also auch bestätigt.
Hypothese 3 ist etwas komplexer, weil sie sich auf einen “Kontrast von einem Kontrast” bezieht: wir testen zunächst die Differenzen zwischen crossed und uncrossed für Level_1 und Level_2 aus factor_B. Dann testen wir, ob die eine dieser Differenzen kleiner ist als die andere. Zum Glück geht das mit dem emmeans-Paket ganz einfach. Auch hier gilt: es muss der p-Wert < 0.05 sein. Und die Differenz (crossed - uncr Level_1) - (crossed - uncr Level_2) muss negativ sein!
simEmm.postureXfactorB <- emmeans(simRTs.rmanova, ~posture|factor_B)
simEmm.postureXfactorB.pairwise <- pairs(pairs(simEmm.postureXfactorB), by = NULL, adjust="none")
diff_val_Hyp3 <- as.data.frame(simEmm.postureXfactorB.pairwise) %>%
filter(contrast=="(crossed - uncr Level1) - (crossed - uncr Level2)") %>%
select("estimate") %>%
as.numeric()
p_val_Hyp3 <- as.data.frame(simEmm.postureXfactorB.pairwise) %>%
filter(contrast=="(crossed - uncr Level1) - (crossed - uncr Level2)") %>%
select("p.value") %>%
as.numeric()
Hyp3_confirmed <- (diff_val_Hyp3 < 0) & (p_val_Hyp3 < 0.05)
Hyp3_confirmed[1] TRUEHypothese 3 ist also auch bestätigt.
Auf diese Weise können wir prüfen, ob alle Hypothesen bestätigt sind. Dies gilt nun also für die eine Stichprobe von 50 Personen, deren Daten wir mit einer bestimmten Effektgröße etc. simuliert haben.
all(Hyp1_confirmed, Hyp2_confirmed, Hyp3_confirmed)[1] TRUE11.10 Eine Simulation starten
11.10.1 Simulation als Funktion
Als nächstes kombinieren wir die Schritte von oben:
- Daten zufällig auswählen
- Statistischess Modell aufstellen
- Hypothesen testen
Dazu nutzen wir eine Funktion.
Die Funktion bekommt als Inputparameter die Range (von - bis) VP, die aus dem (generierten) Datensatz zufällig gezogen werden sollen, den Datensatz selbst und ein Label für die Iteration.
Der Label dient als Identifier für die Iteration, er hat aber an sich keinen Einfluss auf die Ergebnisse.
simulateOneIteration <- function(range.subjects, simdata, iteration=1) {
simulation.results.iteration <- data.frame() #container
subjects <- unique(simdata$participant)
for (k in range.subjects) {
#Daten selegieren
subjects.now <- sample(subjects, size = k, replace = F)
data.powerEst.now <- simdata %>%
filter(participant %in% subjects.now) %>%
mutate(posture=factor(posture, levels=c("crossed", "uncr")))
# statistisches Modell definieren: repeated measures Anova
simRTs.rmanova <- aov_ez(id = "participant", dv = "meanRT",
within = c("factor_B", "posture"),
data = data.powerEst.now)
# die 3 hypothesen testen
# Hypothese 1
p_val_Hyp1 <- as.data.frame(simRTs.rmanova$anova_table) %>%
rownames_to_column() %>%
filter(rowname=="factor_B:posture") %>%
select("Pr(>F)") %>%
as.numeric()
Hyp1_confirmed <- c(p_val_Hyp1 < 0.005)
# Hypothese 2
simEmm.factorBXposture <- emmeans(simRTs.rmanova, ~factor_B|posture)
simEmm.factorBXposture.pairwise <- pairs(simEmm.factorBXposture, by = NULL, adjust="none")
diff_val_Hyp2 <- as.data.frame(simEmm.factorBXposture.pairwise) %>%
filter(contrast=="Level1 crossed - Level2 crossed") %>%
select("estimate") %>%
as.numeric()
p_val_Hyp2 <- as.data.frame(simEmm.factorBXposture.pairwise) %>%
filter(contrast=="Level1 crossed - Level2 crossed") %>%
select("p.value") %>%
as.numeric()
Hyp2_confirmed <- (diff_val_Hyp2 < 0) & (p_val_Hyp2 < 0.05)
# Hypothese 3
simEmm.postureXfactorB <- emmeans(simRTs.rmanova, ~posture|factor_B)
simEmm.postureXfactorB.pairwise <- pairs(pairs(simEmm.postureXfactorB), by = NULL, adjust="none")
diff_val_Hyp3 <- as.data.frame(simEmm.postureXfactorB.pairwise) %>%
filter(contrast=="(crossed - uncr Level1) - (crossed - uncr Level2)") %>%
select("estimate") %>%
as.numeric()
p_val_Hyp3 <- as.data.frame(simEmm.postureXfactorB.pairwise) %>%
filter(contrast=="(crossed - uncr Level1) - (crossed - uncr Level2)") %>%
select("p.value") %>%
as.numeric()
Hyp3_confirmed <- (diff_val_Hyp3 < 0) & (p_val_Hyp3 < 0.05)
#alle Hypothesen erfüllt?
allHyp_confirmed <- all(Hyp1_confirmed, Hyp2_confirmed, Hyp3_confirmed)
#Eregbnisse in dem container speichern
simulation.results.iteration <- rbind(simulation.results.iteration,
data.frame(iteration=iteration,
n.subjects=k,
Hyp1_confirmed=Hyp1_confirmed,
Hyp2_confirmed=Hyp2_confirmed,
Hyp3_confirmed=Hyp3_confirmed,
allHyp_confirmed=allHyp_confirmed))
}
return(simulation.results.iteration)
}Jetzt probieren wir die Funktion aus, um den Output zu verstehen.
range.subjects <- c(5:10)
simulateOneIteration(range.subjects=range.subjects, simdata=dSimulated) iteration n.subjects Hyp1_confirmed Hyp2_confirmed Hyp3_confirmed
1 1 5 FALSE TRUE TRUE
2 1 6 FALSE FALSE TRUE
3 1 7 FALSE FALSE FALSE
4 1 8 TRUE TRUE TRUE
5 1 9 FALSE TRUE TRUE
6 1 10 TRUE TRUE TRUE
allHyp_confirmed
1 FALSE
2 FALSE
3 FALSE
4 TRUE
5 FALSE
6 TRUEIn der letzten Spalte kann man nun also sehen, ob in dieser Iteration alle Hypothesen erfüllt wurden. Aber: das ist nur eine Iteration, das heißt nur eine einzige Zufallsziehung von “n.subjects”. Was wir eigentlich wissen wollen, ist, wie oft die Hypothesen bestätigt werden, wenn wir ganz oft zufällig “n.subjects” aus dem generierten Datensatz ziehen.
11.10.2 Simulation in einer Schleife
Als nächstes passiert also genau das: wir ziehen immer wieder einen Zufallsdatensatz und testen die Hypothesen.
Damit das einigermaßen effizient passiert, verwenden wir einen parallelen “for”-Loop. Das bedeutet, dass der Computer mehrere Prozessoren gleichzeitig verwendet, um zu simulieren. Das verkürzt die Berechnungsdauer um ein Vielfaches.
Zunächst müssen wir auch einige Parameter festlegen:
- Wie viele Personen
- Wie viele Iterationen
- Den Dateinamen in den die Simulationen gespeichert werden
Es gibt hier auch eine Option: “rerunSimulations”. Wenn diese auf FALSE gesetzt wird, werden die Simulationen nicht erneut berechnet. Das macht dann Sinn, wenn man nur noch mit den Ergebnissen rechnen will. Sie werden dann aus der Datei geladen.
# Simulationen rechnen?
rerunSimulations <- F
# Datensatz
simdata=dSimulated
# Parameters für die Simulation
n.iterations = 250
range.subjects = c(3:20)
sim.filename = "Simulations/Power_estimation_simulations.RData"
# XXX simulate all iterations using parallel backend
#time counter
startTime <- Sys.time()
#set up cluster for parallel computing
cl <- parallel::makeCluster(detectCores()) # activate the cluster for the foreach library
doParallel::registerDoParallel(cl)
if (rerunSimulations==T) {
nModels <- n.iterations * length(range.subjects)
print(sprintf("the total number of models to compute is %s", nModels))
# Loop vorbereiten
simulation.results.table <-
foreach(i=1:n.iterations,
.combine = "rbind",
.packages = c("afex","emmeans","dplyr", "tibble")) %dopar%
{
# Infos über Fortschritt werden in eine Textdatei geschrieben
sink("Simulations/sim_iterations.txt", append=T) # Sink-File öffnen und schreiben
cat(paste("iteration", i))
cat("\n")
sink() #end output
# Eine Iteration wird an die Funktion übergeben
simulateOneIteration(range.subjects = range.subjects, simdata = simdata, iteration = i)
}
endTime <- Sys.time()
# Ergebnisse speichern
save(simulation.results.table, startTime, endTime,
file = sim.filename)
print("DONE")
}
if (rerunSimulations==FALSE) {
load(sim.filename)
}
# stop cluster
parallel::stopCluster(cl)Wie lange hat das gedauert?
print("The simulation has taken") [1] "The simulation has taken"endTime - startTimeTime difference of 2.170169 mins11.11 Ergebnisse der Simulation
Als nächstes können wir berechnen, bei welchem Anteil der Iterationen die Hypothesen erfüllt wurden. Dies entspricht dann der Power. Die Ergebnisse werden grafisch dargestellt.
simulation.results.table.summary <- simulation.results.table %>%
group_by(n.subjects) %>%
summarise(power=mean(allHyp_confirmed)*100) %>%
ungroup()
Wie man in der Abbildung sehen kann, werden die 95%-Power (rote Line) ab etwa 15 VPn erreicht.
11.12 Poweranalysen mit selbst-definierten Parametern aus einer Excel-Datei
In vielen Fällen liegt vielleicht kein Datensatz vor, aus dem die Parameter (Mittelwerte, SD, Korrelationen) berecnet werden können. Oder Sie wollen vielleicht mit Absicht von den Parametern abweichen.
Hier geben wir auch dafür noch ein Beispiel und präsentieren auch einen etwas anderen Workflow, bei dem wir die nötigen Parameter in eine Exceltabelle schreiben. Diese Tabelle lesen wir dann in R ein und verfahren ab da, wie zuvor.
Mit diesem Workflow sollte es Ihnen einfach fallen, eine Powerschätzung selbst zu machen. Sie müssen dazu nur die Exceltabelle und die Hypothesen auf Ihren Fall anpassen.
11.12.1 Struktur der Excel-Tablle
Es ist dabei sehr wichtig, dass die Exceltabelle nach einem bestimmten Muster aufgebaut ist. Lesen Sie dazu unbedingt die Box durch.
Wie die Parameter-Exceltabelle aufgebaut sein muss
Die Exceltabelle muss die Information enthalten, die auch in den oben verwendeten Parameter-Datensätzen enthalten war, d.h., Mittelwerte, Standardabweichungen und Interkorrelationen.
Zur Erinnerung, die oben verwendeten Parameter sahen so aus:
dcSim.Params$param.means
Level1_uncr Level1_crossed Level2_uncr Level2_crossed
729.1838 840.1628 716.1749 935.3435
$param.sds
Level1_uncr Level1_crossed Level2_uncr Level2_crossed
144.87466 170.39388 98.79463 180.87108
$param.cor
Level1_uncr Level1_crossed Level2_uncr Level2_crossed
Level1_uncr 1.0000000 0.7923623 0.8492734 0.5111338
Level1_crossed 0.7923623 1.0000000 0.6864440 0.8796122
Level2_uncr 0.8492734 0.6864440 1.0000000 0.5186879
Level2_crossed 0.5111338 0.8796122 0.5186879 1.0000000Das wird dann in der drei Tabellenblätter der Exceltabelle realisiert. Das erste Blatt enthält die Mittelwerte, das zweite die Standardabweichungen und das dritte die Korrelationen.
Dabei ist es wichtig, dass die Variablenstufen hier mit einem Unterstrich voneinander getrennt sind und die Benennung 100% konsistent ist.
Die Exceltabelle müsste also in etwa so aussehen. Der Einfachheit halber sind in der Tabelle mal gerundet die Werte übernommen, die wir vorher hatten.
Beachten Sie, dass man in Excel üblicherweise die Werte mit einem Komma eingibt, nicht einem Punkt, wie in R.
| Blatt 1: Mittelwerte | Blatt 2: SDs | Blatt 3: Korrelationen |
|---|---|---|
 |
 |
 |
Sie sehen also:
die Paraemter sind auf drei Blätter der gleichen xlsx-Datei verteilt (das sieht man an den Tabs unten in der Tabelle)
die Benenung ist konsistent: Faktor1Stufe_Faktor2Stufe
Zudem fällt auf, dass der untere Teil der Korrelationsmatrix freigelassen wurde. Das bietet sich an, weil die Matrix sowieso redundant ist, da sie ja an der Diagonale (mit den 1-Weerten) gespiegelt ist.
11.12.2 Tabelle in R einlesen
Wie in Kapitel 3 beschrieben, können wir die Exceltabelle einfach mit dem Paket “rio” einlesen. Welches Blatt wir einlesen, wird in rio über den Parameter “which” definiert.
library(rio)
customSim.Params.means <- import("Simulations/customSimParams.xlsx", which = 1)
customSim.Params.sds <- import("Simulations/customSimParams.xlsx", which = 2)
customSim.Params.cor <- import("Simulations/customSimParams.xlsx", which = 3)
#turn into named numeric vectors
customSim.Params.Vars <- names(customSim.Params.means)
customSim.Params.means <- as.numeric(customSim.Params.means)
customSim.Params.sds <- as.numeric(customSim.Params.sds)
names(customSim.Params.means) <- customSim.Params.Vars
names(customSim.Params.sds) <- customSim.Params.Vars
#the correlation matrix needs some reshaping
customSim.Params.cor <- as.matrix(customSim.Params.cor[, -1]) #remove first column
rownames(customSim.Params.cor) <- colnames(customSim.Params.cor)
#fill up the part of the matrix below the diagonal ´
customSim.Params.cor <- forceSymmetric(customSim.Params.cor, uplo = "U")Jetzt können wir uns die Matrixen noch einmal ansehen.
customSim.Params.means Level1_uncr Level1_crossed Level2_uncr Level2_crossed
729.18 840.16 716.17 935.34 customSim.Params.sds Level1_uncr Level1_crossed Level2_uncr Level2_crossed
144.87 170.39 98.79 180.87 customSim.Params.cor4 x 4 Matrix of class "dsyMatrix"
Level1_uncr Level1_crossed Level2_uncr Level2_crossed
Level1_uncr 1.00 0.79 0.85 0.51
Level1_crossed 0.79 1.00 0.67 0.88
Level2_uncr 0.85 0.67 1.00 0.52
Level2_crossed 0.51 0.88 0.52 1.0011.12.3 Wiederholung der Poweranalyse mit den eingelsenen Daten
Der Rest ist eigentlich wie oben, außer, dass wir die eingelesenen Parameter verwenden. Zunächst erzeugen wir wieder korrelierte Daten.
Es müssen auch hier wieder die gewünschten Variablennamen (“factor_B” und “posture”) angegeben werden.
nVPn = 1000
dCustomSimulated <- simulateNewData(param.means = customSim.Params.means,
param.sds = customSim.Params.sds,
param.cor = customSim.Params.cor,
n.subjects = nVPn,
id.var = "participant",
condition.vars = c("factor_B", "posture"),
value.var = "meanRT")Wie zuvor, schauen wir uns die Werte mal graphisch an.
p.density.simulated <- dCustomSimulated %>%
mutate(posture=factor(posture, levels=c("uncr", "crossed")),
factor_B=factor(factor_B, levels=c("Level1", "Level2"), labels=c("Level_1", "Level_2"))) %>%
ggplot(aes(x=meanRT, fill=factor_B)) +
geom_density(alpha=0.5) +
facet_wrap(~posture)
p.density.empirical + ggtitle("empirical") + xlim(c(0, 1600)) +
p.density.simulated + ggtitle("simulated") + xlim(c(0, 1600)) +
plot_layout(ncol=1, guides = "collect") 
Da wir die gleichen Faktoren haben wie oben und mal hypothetisch die gleichen Hypothesen prüfen, können wir einfach wieder die Funktion von oben verewenden, um die Poweranalyse durchzuführen.
range.subjects <- c(5:10)
simulateOneIteration(range.subjects=range.subjects, simdata=dCustomSimulated) iteration n.subjects Hyp1_confirmed Hyp2_confirmed Hyp3_confirmed
1 1 5 FALSE FALSE FALSE
2 1 6 FALSE TRUE TRUE
3 1 7 TRUE TRUE TRUE
4 1 8 TRUE TRUE TRUE
5 1 9 TRUE TRUE TRUE
6 1 10 FALSE TRUE TRUE
allHyp_confirmed
1 FALSE
2 FALSE
3 TRUE
4 TRUE
5 TRUE
6 FALSEDa das funktioniert, funktioniert auch die Schleife, in der die Poweranalyse dann durchgeführt wurde. Das einzige Wesentliche, was wir ändern, ist der Datensatz: dCustomSimulated
Da wir auch schon wissen, dass Stichprobengrößen unter 10 irrelevant sind, passen wir auch das range.subjects-Argument an. Das definiert, was für Stichprobengrößen verwendet werden.
# Simulationen rechnen?
rerunSimulations <- F
# Datensatz
simdata=dCustomSimulated
# Parameters für die Simulation
n.iterations = 250
range.subjects = 10:20
sim.filename = "Simulations/Power_estimation_simulations_custom.RData"
# XXX simulate all iterations using parallel backend
#time counter
startTime <- Sys.time()
#set up cluster for parallel computing
cl <- parallel::makeCluster(detectCores()) # activate the cluster for the foreach library
doParallel::registerDoParallel(cl)
if (rerunSimulations==T) {
nModels <- n.iterations * length(range.subjects)
print(sprintf("the total number of models to compute is %s", nModels))
# Loop vorbereiten
simulation.results.table <-
foreach(i=1:n.iterations,
.combine = "rbind",
.packages = c("afex","emmeans","dplyr", "tibble")) %dopar%
{
# Infos über Fortschritt werden in eine Textdatei geschrieben
sink("Simulations/sim_iterations.txt", append=T) # Sink-File öffnen und schreiben
cat(paste("iteration", i))
cat("\n")
sink() #end output
# Eine Iteration wird an die Funktion übergeben
simulateOneIteration(range.subjects = range.subjects, simdata = simdata, iteration = i)
}
endTime <- Sys.time()
# Ergebnisse speichern
save(simulation.results.table, startTime, endTime,
file = sim.filename)
print("DONE")
}
if (rerunSimulations==FALSE) {
load(sim.filename)
}
# stop cluster
parallel::stopCluster(cl)Analog zu oben lönnen wir wieder die Ergebnisse plotten.
simulation.results.table.summary <- simulation.results.table %>%
group_by(n.subjects) %>%
summarise(power=mean(allHyp_confirmed)*100) %>%
ungroup()
Sie sind herzlich eingeladen, mal ein bisschen herumzuprobieren, indem Sie Werte in der Exceltabelle verändern. Was passiert, wenn man die Mittelwertsunterschiede größer oder kleiner macht oder die Streuungen verändert? Probieren Sie es doch einfach mal aus!
Das ist auch ein gute Überleitung zum nächsten Abschnitt über die Sensitivitätsanalyse. Denn diese fragt genau systematisch danach: wie klein kann ich die Effekte machen und sie trotzdem in der Stichprobe mit dem gewählten Design entdecken?
Ein bisschen Vorsicht ist geboten bei der Anpassung der Korrelationsmatrix. Besonders, wenn man hohe Korrelationen und eine weite Streung an Korrelationsen einträgt, kann es sein, dass es einfach nicht möglich ist, entsprechende Daten zu simulieren. Das liegt daran, dass bestimmte Korrelationsmatrizen einfach nicht möglich sind. Überlegen Sie mal: Sie hätten Variablen A, B und C. A soll positiv mit B und C korrelieren. Sie tragen aber auch ein, dass B und C untereinander negativ korrelieren sollen. Wie soll das gehen? R wird in solchen Fällen eine Fehlermeldung geben. Überlegen Sie also gut und verwenden Sie realistische Parameter.
11.13 Sensitivitätsanalyse
Manchmal hat man bereits Daten. Vielleicht sind dort signifikante Effekte drin, vielleicht auch nicht. Da die Daten vorliegen, kann man die Effektgröße ermitteln.
| Signifikante Effekte |
|
| Keine Signifikanten Effekte |
|
Und zuletzt kann man alles kombinieren und fragen:
Mit der Varianz, die in meinen Daten vorhanden ist, welche Power hatte ich da für welche Effektgröße?
Denken Sie nochmal daran, dass die Power sich ergibt aus der Anzahl der Vpn, der Effektgröße und der Varianz. Nun liegt die Varianz schon fest, da die Daten schon vorhanden sind. Entsprechend können wir für jede angenommene Effektgröße ermitteln, wie groß die Power in unserer Studie (mit N Vpn und der beobachteten Varianz) war.
11.13.1 Inkrementierung von Daten
Der Clou bei der Sensitivitätsanalyse liegt in der schrittweisen (oder “inkrementellen”) Verkleinerung (oder wenn man möchte auch Vergrößerung) des Effektes.
Weil das häufig wiederholt wird, gibt es auch hierzu eine Hilfsfunktion, . Die Funktion ist sehr simpel. Sie verändert einfach die Mittelwertsparameter, die verwendet werden, um Daten zu simulieren. Erinnern Sie sich an den letzten Abschnitt, in dem es um das Einlesen von Parametern aus einer Exceltabelle ging, in der man die Parameter setzen kann wie man möchte? Die Funktion macht genau das: sie verändert diese Parameter, genauer, die Größe der Mittelwertsabstände.
Ein wichtiges Argument, das hier gesetzt werden muss, ist “increment.means”. Das Inkrement drückt aus, wie groß der Effekt sein soll, also, wie stark er verändert werden soll. Wenn der Wert 1 ist, würde man nichts verändern. Wenn man ihn auf 0.5 setzt, wäre er noch halb so groß, bei einem Wert von 2 doppelt so groß und bei 0 abwesend. Für Beispiele und Abbildungen, siehe die Box unten.
Wie funktioniert die Inkrementierung?
Das Argument “increment.means” ist ein Skalierungsfaktor für die Mittelwertsparameter. Bei einem Wert von 1 bleiben die Mittelwerte gleich, bei Werten kleiner als 1 rücken die Verteilungen näher zusammen und bei Werten über 1, rücken sie auseinander. Hier eine Demo.
Die Parameter, die wir in der Stichprobe gefunden hatten, haben wir oben mit der Funktion “simulateSampleParams” berechnet und “dcSim.Params” genannt.
dSimIncremented <- data.frame() #empty container
for (i in c(0, 0.1, 0.5, 1, 1.5, 3)){
dSimIncremented.now <- simulateNewData(param.means = dcSim.Params$param.means,
param.sds = dcSim.Params$param.sds,
param.cor = dcSim.Params$param.cor,
n.subjects = nVPn,
id.var = "participant",
condition.vars = c("factor_B", "posture"),
value.var = "meanRT",
increment.means = i)
dSimIncremented.now$increment <- i
dSimIncremented <- rbind(dSimIncremented, dSimIncremented.now)
}
#plot for different increments
dSimIncremented %>%
filter(factor_B=="Level1") %>%
ggplot(aes(x=meanRT, color=posture)) +
geom_density(alpha=0.5) +
facet_wrap(~increment) +
scale_color_manual(values=c("seagreen", "orange")) +
ggtitle("Data for different increments")
In der Abbildung sieht man, wie bei einem Inkrement von 0 die Bedingungen genau aufeinander liegen, also gar kein Effekt mehr da ist. Bei einem Inkrement von 1 sind die Effekte theoretisch wie sie tatsächlich waren und bei Werten größer als 1 werden sie übertrieben. Der Darstellung halber ist nur eine Stufe von factor_B dargestellt.
Man kann es auch noch anders darstellen, um den Verlauf der Versuchspersonen über die Bedingungen klarer zu machen:
#plot for different increments
dSimIncremented %>%
filter(participant %in% paste0("id", sprintf("%03d", 1:20))) %>%
mutate(factor_B=factor(factor_B, levels=c("Level1", "Level2"),
labels=c("L1", "L2")),
posture=factor(posture, levels=c("uncr", "crossed"),
labels=c("ucr", "cr")),
condition=paste(factor_B, posture, sep="")) %>%
ggplot(aes(x=condition, y=meanRT)) +
geom_point() +
geom_line(aes(group=participant)) +
stat_summary(fun=mean, geom="point", size=2, color="red") +
facet_wrap(~increment) +
ggtitle("Simulated data for different increments")
11.13.2 Schrittweise Inkrementierung der Daten für die Sensitivitätsanalyse
Die folgenden Schritte sind ähnlich wie bei der Poweranalyse in dem Sinne, dass wir prüfen, ob wir diese Hypothesen bei den simulierten Daten bestätigen oder nicht. Der Unterschied ist, dass wir nicht die Stichprobengröße reduzieren, sondern die Effektgröße. Wir “inkrementieren” uns also von der originalen Effektgröße in kleinen Schritten nach unten und schauen, ab wann die Statistik nicht mehr signifikant wird.
Für die Implementierung können wir die zuvor verwendete Funktion “simulateOneIteration” verwenden. Wir stellen aber einfach die VP-Anzahl, die aus dem Datensatz gezogen wird auf 18, was der tatsächlichen Stichprobengröße entspricht.
Demo, wie die Funktion auf die inkrementierten Daten angewendet werden kann
Wir verwenden auch gleich mal ein Inkrement von 0.25:
set.seed(1111) #make reproducible
dSimIncremented.now <- simulateNewData(param.means = dcSim.Params$param.means,
param.sds = dcSim.Params$param.sds,
param.cor = dcSim.Params$param.cor,
n.subjects = nVPn,
id.var = "participant",
condition.vars = c("factor_B", "posture"),
value.var = "meanRT",
increment.means = 0.25)Mit dem Datensatz kann man die ANOVA rechnen (beachte aber, dass er 1000 VP hat):
aov_ez(id = "participant", dv = "meanRT",
within = c("factor_B", "posture"),
data = dSimIncremented.now)Anova Table (Type 3 tests)
Response: meanRT
Effect df MSE F ges p.value
1 factor_B 1, 999 5412.77 18.62 *** .001 <.001
2 posture 1, 999 15332.06 131.43 *** .022 <.001
3 factor_B:posture 1, 999 1597.83 121.46 *** .002 <.001
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1Und natürlich den oben definierten Test, ob sich die Hypothesen bestätigen.
simulateOneIteration(range.subjects=18, simdata=dSimIncremented.now) iteration n.subjects Hyp1_confirmed Hyp2_confirmed Hyp3_confirmed
1 1 18 FALSE FALSE FALSE
allHyp_confirmed
1 FALSEMan sieht, dass wir (in dieser Zeihung) mit 18 VP bei einer um 75% reduzierten Effektgröße die Hypothesen nicht mehr bestätigen können. Das ist wenig überraschend, denn eine Verringerung um 75% ist sehr viel.
11.13.2.1 Increments einsetzen in einer Schleife
Wieder verwenden wir eine Schleife… Sie ist ähnlich aufgebaut wie oben, aber:
- wir verwenden immer 18 VP, d.h., range.subjects ist festgesetzt auf 18
- wir verwenden logischerweise verschiedene Inkremente, d.h. es ist noch ein weiter Parameter in der Schleife implementiert.
# Simulationen rechnen?
rerunSimulations <- F
# Datensatz
simdata.params=dSimulated.Params
# Parameters für die Simulation
increments <- seq(1, 0, length.out=101) #1 to zero in 0.01 increments
n.iterations = 250
range.subjects = c(18)
sim.filename = "Simulations/Sensitivity_estimation_simulations.RData"
# XXX simulate all iterations using parallel backend
#time counter
startTime <- Sys.time()
#set up cluster for parallel computing
cl <- parallel::makeCluster(detectCores()) # activate the cluster for the foreach library
doParallel::registerDoParallel(cl)
if (rerunSimulations==T) {
nModels <- n.iterations * length(range.subjects) * length(increments)
print(sprintf("the total number of models to compute is %s", nModels))
simulation.results.table <- data.frame()
for (inc in increments) {
simdata <- simulateNewData(param.means = simdata.params$param.means,
param.sds = simdata.params$param.sds,
param.cor = simdata.params$param.cor,
n.subjects = nVPn,
id.var = "participant",
condition.vars = c("factor_B", "posture"),
value.var = "meanRT",
increment.means = inc)
# Loop vorbereiten
simulation.results.table.iteration <-
foreach(i=1:n.iterations,
.combine = "rbind",
.packages = c("afex","emmeans","dplyr", "tibble")) %dopar%
{
# Infos über Fortschritt werden in eine Textdatei geschrieben
sink("Simulations/sim_iterations.txt", append=T) # Sink-File öffnen und schreiben
cat(paste("iteration", i))
cat("\n")
sink() #end output
# Eine Iteration wird an die Funktion übergeben
simulateOneIteration(range.subjects = range.subjects, simdata = simdata, iteration = i)
}
simulation.results.table.iteration$increment <- inc
#attach to container
simulation.results.table <- rbind(simulation.results.table, simulation.results.table.iteration)
}
endTime <- Sys.time()
# Ergebnisse speichern
save(simulation.results.table, startTime, endTime,
file = sim.filename)
print("DONE")
}
if (rerunSimulations==FALSE) {
load(sim.filename)
}
# stop cluster
parallel::stopCluster(cl)11.13.2.2 Ergebnisse visualisieren
Das funktioniert genau wie zuvor.
simulation.results.table.summary <- simulation.results.table %>%
group_by(increment) %>%
summarise(power=mean(allHyp_confirmed)*100) %>%
ungroup()
Man kann in der Abbildung nun sehen, dass die Effekten nicht viel kleiner sein dürften als sie waren, um bei der verwendeten Stichprobengröße und Testdesign die Hypothesen bestätigen zu konnen. Wenn z.B. der Effekt nur 75% der beobachteten Größe wäre, hätte man auch nur noch eine Power von etwa 75%.
11.13.3 Sensitivitäts-Analyse für verschiedene Versuchspersonen-Anzahl
Zuletzt wollen wir noch fragen: Geht eigentlich auch Variation der Effektgröße und Stichprobengröße auf einmal? Sehr einfach sogar. Man muss nur in der Schleife oben wieder die range.subjects hochsetzen. Allerdings braucht man eine andere Visualisierung.
Hier ein Beispiel. Wir verwenden aber weniger Inkrements, da die Anzahl an Simulationen sonst sehr groß wird.
# Simulationen rechnen?
rerunSimulations <- F
# Datensatz
simdata.params=dSimulated.Params
# Parameters für die Simulation
increments <- seq(1, 0, length.out=11) #1 to zero in 0.1 increments
n.iterations = 250
range.subjects = c(10:25)
sim.filename = "Simulations/Sensitivity_estimation_simulations2.RData"
# XXX simulate all iterations using parallel backend
#time counter
startTime <- Sys.time()
#set up cluster for parallel computing
cl <- parallel::makeCluster(detectCores()) # activate the cluster for the foreach library
doParallel::registerDoParallel(cl)
if (rerunSimulations==T) {
nModels <- n.iterations * length(range.subjects) * length(increments)
print(sprintf("the total number of models to compute is %s", nModels))
simulation.results.table <- data.frame()
for (inc in increments) {
simdata <- simulateNewData(param.means = simdata.params$param.means,
param.sds = simdata.params$param.sds,
param.cor = simdata.params$param.cor,
n.subjects = nVPn,
id.var = "participant",
condition.vars = c("factor_B", "posture"),
value.var = "meanRT",
increment.means = inc)
# Loop vorbereiten
simulation.results.table.iteration <-
foreach(i=1:n.iterations,
.combine = "rbind",
.packages = c("afex","emmeans","dplyr", "tibble")) %dopar%
{
# Infos über Fortschritt werden in eine Textdatei geschrieben
sink("Simulations/sim_iterations.txt", append=T) # Sink-File öffnen und schreiben
cat(paste("iteration", i))
cat("\n")
sink() #end output
# Eine Iteration wird an die Funktion übergeben
simulateOneIteration(range.subjects = range.subjects, simdata = simdata, iteration = i)
}
simulation.results.table.iteration$increment <- inc
#attach to container
simulation.results.table <- rbind(simulation.results.table, simulation.results.table.iteration)
}
endTime <- Sys.time()
# Ergebnisse speichern
save(simulation.results.table, startTime, endTime,
file = sim.filename)
print("DONE")
}
if (rerunSimulations==FALSE) {
load(sim.filename)
}
# stop cluster
parallel::stopCluster(cl)11.13.3.1 Visualisierung der Ergebnisse in einer Heatmap
Weil wir hier zwei Dimensionen auf einmal haben, bietet es sich an, eine sogenannte Heatmap zu verwenden.
simulation.results.table.summary <- simulation.results.table %>%
group_by(n.subjects, increment) %>%
summarise(power=mean(allHyp_confirmed)*100) %>%
ungroup()`summarise()` has grouped output by 'n.subjects'. You can override using the
`.groups` argument.ggplot(data = simulation.results.table.summary, aes(x=n.subjects,
y=as.factor(increment), fill=power)) +
geom_tile() +
geom_text(aes(label = sprintf("%2.0f", power)), color="grey20", size=5) +
scale_fill_gradient2(low = "black", mid = "red", high="yellow", midpoint = 50) +
scale_x_continuous(breaks=range.subjects) +
xlab("number of participants") + ylab("proportion of effect size") 
Diese Abbildung enthält eine Menge Information. Wenn ich zum Beispiel eine Power von 80% anstrebe und annehme, dass der Effekt etwas kleiner, sagen wir 75%, sein könnte, dann sollte ich zumindest 20 VP messen. Noch wichtiger: die Abbildung zeigt, dass viele mögliche Designs absolult zwecklos wären. Mit einer Studie mit wenigen Versuchspersonen (z.B. 10, 11, 12) ist wenig zu erreichen, weil sie nicht genug Power hätten, besonders wenn der Effekt auch nur minimal kleiner als der beobachtete ist.