PowerAnalysis

PowerAnalysis.R

library(parameters)
library(ggplot2)


# Einfacher t-Test ------------------------------

# Poweranalyse in R, für einen t-Test, mit n = 50 und Effekstärke = 0.3
pwr::pwr.t.test(n = 50, d = .3)

# Ergebnis: eine Power von knapp 32% - d.h., wenn wir wiederholt Stichproben
# mit n = 100 (2x50) ziehen, und sich zwei Gruppen bzgl. eines Mittelwertes
# um 0.3 unterscheiden, erwarten wir aufgrund der Varianz in der Stichprobe
# nur in 32% aller Fälle einen signifikanten Unterschied zwischen diesen beiden
# Gruppen. Genau das machen wir im Folgenden mit der Simulation: wir ziehen
# wiederholt (1000x) je 2x50 Fälle, einmal n=50 mit MW 0 und einmal n=50 mit
# MW 0.3 (d.h. der Unterschied bzw. die "Effektstärke" ist 0.3). Dann berechnen
# wir einen t-Test und speichern den p-Wert für jede Wiederholung. Am Ende
# müssten nur ca. 32% aller p-Werte unter 0.05 liegen.



# Wir machen das erstmal nur einmal, und schauen uns die Verteilung an
set.seed(123)
gruppe1 <- data.frame(x = rnorm(50, mean = 0), gruppe = as.factor(1))
gruppe2 <- data.frame(x = rnorm(50, mean = 0.3), gruppe = as.factor(2))

# zusammenfügen beider Datensätze
sim_data <- rbind(gruppe1, gruppe2)
ggplot(sim_data, aes(x = x, group = gruppe, fill = gruppe)) +
  stat_density(alpha = .2, position = position_dodge(width = 0)) +
  geom_vline(xintercept = mean(gruppe1$x), color = see::flat_colors("blue"), alpha = .5) +
  geom_vline(xintercept = mean(gruppe2$x), color = see::flat_colors("red"), alpha = .5) +
  see::scale_fill_flat()

ggplot(sim_data, aes(x = x, group = gruppe, fill = gruppe)) +
  stat_density(alpha = .2, position = position_dodge(width = 0)) +
  geom_vline(xintercept = median(gruppe1$x), color = see::flat_colors("blue"), alpha = .5) +
  geom_vline(xintercept = median(gruppe2$x), color = see::flat_colors("red"), alpha = .5) +
  see::scale_fill_flat()



# Exkurs: Zusammenhang t-Test und lineare Regression ----------------

# Achtung! Der neue "pipe"-Operator: |> setzt R Version 4.1 voraus!
# Bei Schwierigkeiten "|>" durch "%>%" ersetzen...

# mean of x = Intercept
# mean of y = Intercept + Koeffizient gruppe[2]
# t-statistic, df und p-Wert vom t-Test entsprechen den Angaben von "gruppe[2]"
t.test(gruppe1$x, gruppe2$x)
lm(x ~ gruppe, data = sim_data) |> model_parameters(digits = 4)



# Exkurs: Zusammenhang Korrelation und lineare Regression ----------------

# t-statistic, df und p-Wert vom t-Test entsprechen den Angaben von "x2"
cor_data <- data.frame(x1 = gruppe1$x, x2 = gruppe2$x)
cor.test(cor_data$x1, cor_data$x2)
lm(x1 ~ x2, data = cor_data) |> model_parameters(standardize = "refit", digits = 4)



# Nun die Simulation...

# anzahl der Simulationen
nsim <- 1000

# p-Werte von Simulationen
all_p <- vector("numeric", nsim)

set.seed(1234) # damit das Ergebnis bei jedem gleich ist.

# Fälle pro Gruppe
n <- 50

# 1000 Iterationen
for (i in 1:nsim) {
  # Zufallsvariable mit n = 50 und MW = 0 erstellen
  x1 <- rnorm(n)

  # Zufallsvariable mit n = 50 und MW = 0.3 (= Effektstärke) erstellen
  x2 <- rnorm(n, mean = .3)

  # t-Test berechnen und p-Werte speichern
  all_p[i] <- t.test(x1, x2)$p.value
}

# Anteil der p-Werte, die kleiner als 0.05 sind. "all_p < 0.05" gibt TRUE bzw.
# FALSE zurück. "sum()" summiert alle "TRUE" auf, also die Anzahl aller "TRUEs".
# Dies teilen wir durch die Gesamtanzahl (length()) aller p-Werte
sum(all_p < .05) / length(all_p)

# oder, bei logischen Variablen, die nur TRUE/FALSE enthalten, geht auch der MW
logische_werte <- all_p < 0.05
logische_werte

# tada - eine Power von etwa 32%
mean(logische_werte)



# n = 50, Effektstärke = 0.5

pwr::pwr.t.test(n = 50, d = .5)
all_p <- vector("numeric", 1000)

for (i in 1:1000) {
  x1 <- rnorm(50)
  x2 <- rnorm(50, mean = .5)
  all_p[i] <- t.test(x1, x2)$p.value
}

mean(all_p < .05)





# n = 50, Effektstärke = 0.5, Alternativhypothese nicht ungleich, sondern "größer"
# dies ist die aktuelle Standardeinstellung in "G*Power"

n <- 51
d <- .5
pwr::pwr.t.test(n = n, d = d, alternative = "greater")
all_p <- vector("numeric", 1000)

for (i in 1:1000) {
  x1 <- rnorm(n)
  x2 <- rnorm(n, mean = d)
  all_p[i] <- t.test(x1, x2, alternative = "less")$p.value
}

mean(all_p < .05)







# Lineare Regression ------------------------------

library(datawizard)
library(parameters)
library(sjPlot)
library(dplyr)


# Daten einlesen
d2 <- data_read("./Daten_mergen_Anna/German Wings match.sav") |> 
  datawizard::to_numeric(dummy_factors = FALSE)

# alle Werte im Minusbereich als missing setzen
d2 <- convert_to_na(d2, na = -99:-1)

# nur diese vier Variablen für später auswählen, keine fehlenden Werte
d <- d2 |>
  data_select(select = c("b14a", "q100", "q101", "e11")) |>
  na.omit() # alle Fälle mit fehlenden Werten entfernen

# Achtung! Der neue "pipe"-Operator: |> setzt R Version 4.1 voraus!
# Bei Schwierigkeiten "|>" durch "%>%" ersetzen...


# wie sieht das "ursprüngliche" Modell aus, das anhand der realen
# Daten berechnet wurde?
m <- lm(b14a ~ q100 + q101 + e11, data = d)
summary(m)
model_parameters(m)


# Exkurs: Zugriff auf einzelne Werte aus dem Summary der Regressionsanalyse ---------------

# model_parameters() gibt einen data frame zurück:
mp <- model_parameters(m)

# die Werte des Summaries der Regressionsanalyse
# sind als Spalten/Variablen gespeichert
colnames(mp)
view_df(mp)
class(mp)

# da der data frame auch noch das Klassenattribut "parameters_model" hat,
# gibt es eine spezielle print() Methode, die die Ausgabe formatiert:
mp

# als "reiner" data frame sieht es anders aus
as.data.frame(mp)

# Zugriff auf Koeffizienten also wie der Zugriff auf eine "Variable"
# im normalen data frame, mit $ oder [[]]
betas <- model_parameters(m)$Coefficient # oder coef(m), oder mp$Coefficients





# Exkurs: wie simuliert man ein Regressionsmodell? -------------------------

# Regressionsformel: y = b0 + b1 * x1 + b2 * x2 + ... + e

# 100 Fälle, x1 ist kontinuierlich mit Mittelwert 10 und SD 2,
# x2 ist binär, 60% 0 und 40% 1
set.seed(1234)

n <- 100
x1 <- rnorm(n, mean = 10, sd = 2)
x2 <- rbinom(n, 1, prob = c(.4 , .6))

# intercept: Mittelwert von y, wenn alle anderen Kovarianten auf Mittelwert
# (kontinuierlich) bzw. auf "0" (kategorial) gesetzt werden
b0 <- 30

# Koeffizient von x1
b1 <- 5

# Koeffizient von x2
b2 <- 2

# Residual-Standardfehler (error term)
sig <- sqrt(abs(1 - b1^2 - b2^2))  # komische Formel, hab ich irgendwo gefunden
y <- rnorm(n, b0 + b1 * x1 + b2 * x2, sig)

model <- lm(y ~ x1 + x2)

# die Koeffizienten des Modells entsprechen den Angaben oben
# zu b0 bis b2...
model_parameters(model)

# intercept + b1 * MW(x1) ~ 30 + 5 * 10 ~ 80 als MW für y
mean(y)



# Power für einzelne Koeffizienten
# Dies kann über eine Anova-Tabelle einfach berechnet werden

m <- lm(b14a ~ q100 + q101 + e11, data = d)
anova_m <- aov(m)
model_parameters(anova_m, power = TRUE)


# Summary des Modells. Wir kennen nicht nur die Koeffizienten, sondern
# auch die Residual-Standardabweichung
model_parameters(m, summary = TRUE)


# Über Simulation kann man das auch prüfen. Wir wissen in diesem Fall bereits,
# wie die "wahren" Koeffizienten sind und nehmen diese für die Simulation,
# d.h. wir "simulieren" unser Modell, rechnen 1000x das Regressionsmodell
# mit neu simulierten Daten und schauen dann, wie oft jeder Koeffizienz
# signifikant war. Dies entspricht der "Power" eines Koeffizienten".

ps <- list()
betas <- model_parameters(m)$Coefficient

# für diese "Simulation" nehmen wir die bekannten Koeffizienten. Theoretisch
# müsste die Power dann vergleichbar mit der obigen Anova-Tabelle sein. Wenn wir
# die Koeffizienten (d.h. Effektstärken für die einzelnen Variablen) nicht aus
# eigenen Daten kennen, müssten wir Annahmen über die Effektstärken einzelner
# Koeffizienten treffen, wie wir sie in der Power-Analyse auch habe (z.B.
# "mindestens mittlere Effektstärke").
b0 <- betas[1]
b1 <- betas[2]
b2 <- betas[3]
b3 <- betas[4]

# "n" = Anzahl der Fälle. "nobs" = number of observations, also Anzahl
# der gültigen Fälle im Modell. Dies sind für unser Modell oben 473 Fälle
n <- nobs(m)

# Residual-Standardabweichung
sig <- sqrt(abs(1 - b1^2 - b2^2 - b3^2))  # komische Formel, hab ich irgendwo gefunden
sig

# die obige Formel ergibt eine Residual-Standardabweichung von 0.998. Die
# "tatsächliche" Residual-Standardabweichung im Modell beträgt 0.848, ist also
# relativ ähnlich. Wir verwenden daher hier die "echte".

sig <- insight::get_sigma(m)

set.seed(1234)

# wir erstellen 1000x Zufallsvariablen, die die gleichen Merkmale wie die
# "echten" Daten haben, also kontinuierliches Alter, normalverteilt mit MW
# des Alters im Datensatz (mean(d$q101) = 50.1) und der entsprechenden
# SD (sd(d$q101) = 18.1) usw. Damít berechnen wir ein Regressionsmodell und
# extrahieren die p-Werte. Da die Daten jedesmal zufällig neu generiert werden
# (insgesamt 1000x), variieren natürlich auch die Ergebnisse der Regressions-
# analysen und damit auch die p-Werte. Die Anteil der p-Werte unter 0.05 ist
# die Power eines Koeffizienten.

for (i in 1:1000) {
  # Variable geschlecht: Bernoulli-Folge (nur 0 und 1), d.h. 473 mal (=n)
  # die Werte von 0 bis 1 (=size) ziehen. Wahrscheinlichkeit für die Werte 0
  # und 1 (=prob) entsprechend der Verteilung im Datensatz
  geschlecht <- rbinom(n, size = 1, prob = proportions(table(d$q100)))

  # Variable alter: kontinuierlich, normalverteilt, mit MW und SD
  # entsprechend der Verteilung im Datensatz
  alter <- rnorm(n, mean = mean(d$q101), sd = sd(d$q101))

  # Variable eonkommen: kontinuierlich, normalverteilt, mit MW und SD
  # entsprechend der Verteilung im Datensatz
  einkommen <- rnorm(n, mean = mean(d$e11), sd = sd(d$e11))

  # y, also die AV, simulieren. Da wir uns auf bekannte Koeffizienten und
  # Variablenmerkmale beziehen (d.h. wir verwenden die Koeffizienzen aus
  # dem Regressionsmodell, das auf dem richtigen Datensatz basiert, und
  # auch Sigma des Modells), sollte die simulierte AV "y" im ähnlichen Bereich
  # liegen wie die tatsächliche AV, "b14a".
  y <- rnorm(n, b0 + b1 * geschlecht + b2 * alter + b3 * einkommen, sig)

  # Hier die Ergebnisse aus einer Zufallsserie als Beispiel
  # - gar nicht so schlecht!
  #
  # > summary(y)
  #   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  # 0.2461  2.3111  2.8189  2.8390  3.4386  5.1365
  # > summary(d$b14a)
  #  Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  # 1.000   2.000   3.000   2.782   3.000   4.000

  # Modell mit simulierten Daten berechnen
  m <- lm(y ~ geschlecht + alter + einkommen)

  # Alle 4 p-Werte in einer Liste speichern
  # ps[[i]] <- as.vector(coef(summary(m))[, 4])
  ps[[i]] <- parameters::p_value(m)$p
}

# "lapply" - apply function over list - ist ähnliche wie eine "For-Schleife"

# Power für Geschlecht

p_geschlecht <- unlist(lapply(ps, function(i) {
  i[2]
}))
mean(p_geschlecht < .05)



# Power für Alter

p_alter <- sapply(ps, function(i) {
  i[3]
})
mean(p_alter < .05)



# Power für Einkommen

p_einkommen <- sapply(ps, function(i) {
  i[4]
})
mean(p_einkommen < .05)


# nochmal zum Vergleich
m <- lm(b14a~ q100 + q101 + e11, data = d)
anova_m <- aov(m)
model_parameters(anova_m, power = TRUE)



# Zusammengefasst
#
# 1.) Wir müssen etwas über unser Outcome wissen, also auf welcher Skala wollen
# wir dies erfassen? 1-100? 0-10? Dichotom 0/1? Danach müssen wir darüber Vermutungen
# angestellen, wie der Mittelwert bzw. bei dichotom die prozentuale Verteilung
# des Outcomes in der gesamten Stichprobe aussehen könnte. Merken.
#
# 2.) Wir müssen wissen, wie die Merkmale (Geschlecht, Alter, Bildung, ...) in der
# Bevölkerung verteilt sind. Dies wissen wir aus Statistiken (Destatis) oder
# aus anderen Studien. Dieses Wissen nutzen wir, um unsere Variablen zu simulieren
# (mit rnorm, rbinom, ...).
#
# 3.) Wir stellen Vermutungen an, wie große die "Effektstärken" dieser Merkmale
# jeweils sind (also die Regressionskoeffizienten). Z.B. Um wie viel Punkte
# auf der Outcome-Skala sollen sich Männer von Frauen unterscheiden? Dies wäre
# dann der Koeffizient b1 bis bn, also die "Effektstärke". Diese geschätzten
# Effektstärken brauchen wir für jede unabhänige Varuable, die wir unter 2)
# simuliert haben.
#
# 4.) Der Intercept "b0" wäre: mean(outcome) - (b1 * mean(x1) + b2 * mean(x2) ...)
# "outcome" muss 1x mit den Merkmalen wie in 1.) beschrieben simuliert werden.
# x1 bis xn müssen 1x wie in 2.) beschrieben simuliert werden.
# Beispiel: unser Outcome is Lebensqualität, bei der wir ausgehen, dass der
# Mittelwert in unserer Stichprobe etwa 70 wäre. Der Anteil an Männer/Frauen
# in unserer Untersuchungspopulation beträgt 40/60. Wir gehen davon aus, dass
# in der Gruppe, die wir untersuchen, das durchschnittliche Alter 50 ist.
# Der "Effekt" von Geschlecht ist 5, d.h. Frauen haben durchschnittlich eine um
# 10 Punkte höhere Lebensqualität -> b1 = 5
# Der "Effekt" von Alter ist -0.5, d.h. mit jedem Jahr, das eine Person älter ist,
# sinkt die durchschnittliche Lebensqualität um 0.5 -> b2 = -0.5
# daraus können wir mit der genannten Formel oben b0 berechnen
# Das Outcome wird jetzt nochmal neu berechnet (wir haben ja alle "betas"), und
# wir fügen die Residual-Standardabweichung (Sigma) hinzu.

b1 <- 5
b2 <- -0.5

lebensqualitaet <- rnorm(100, mean = 70, sd = 20)
geschlecht <- rbinom(100, size = 1, prob = c(.4, .6))
alter <- rnorm(100, mean = 50, sd = 15)

b0 <- mean(lebensqualitaet) - (b1 * mean(geschlecht) + b2 * mean(alter))

sig <- sqrt(abs(1 - b1^2 - b2^2))
lebensqualitaet <- rnorm(100, b0 + b1 * geschlecht + b2 * alter, sig)

# jetzt rechnen wir ein Modell mit den simulierten Daten
model <- lm(lebensqualitaet ~ geschlecht + alter)
model_parameters(model)

# je nach "Zufallsstichprobe", erhalten wir hierfür signifikante Werte für
# Alter und Geschlecht. Wenn wir das jetzt 1000x wiederholen, haben wir die
# Power für die jeweiligen Variablen berechnet, für die aktuelle Stichprobengröße
# n = 100.

b1 <- 5
b2 <- -0.5

lebensqualitaet <- rnorm(100, mean = 70, sd = 20)
geschlecht <- rbinom(100, size = 1, prob = c(.4, .6))
alter <- rnorm(100, mean = 50, sd = 15)

b0 <- mean(lebensqualitaet) - (b1 * mean(geschlecht) + b2 * mean(alter))
sig <- sqrt(abs(1 - b1^2 - b2^2))

p_werte <- list()
for (i in 1:1000) {
  geschlecht <- rbinom(100, size = 1, prob = c(.4, .6))
  alter <- rnorm(100, mean = 50, sd = 15)
  lebensqualitaet <- rnorm(100, b0 + b1 * geschlecht + b2 * alter, sig)
  model <- lm(lebensqualitaet ~ geschlecht + alter)
  ps[[i]] <- parameters::p_value(model)$p
}

# Power für Geschlecht
p_geschlecht <- unlist(lapply(ps, function(i) {
  i[2]
}))
mean(p_geschlecht < .05)

# Power für Alter
p_alter <- sapply(ps, function(i) {
  i[3]
})
mean(p_alter < .05)









# power-Analyen für Regressions-Modelle

k <- 8 # anzahl der Prädiktoren
n <- 60

pwr::pwr.f2.test(u = 1, v = n - k, f2 = .15, sig.level = 0.05)



pwr::pwr.f2.test(u = 1, f2 = .15, sig.level = 0.05, power = 0.8)
# sample size = u + v + 2