Parameter

Wir gehen von einem linearen Einfluss des Starrens auf mu aus.

Effekt skalieren

effect_size <- -1 # 1 Sekunde weniger Zeit

Bisher zeigt der Term b * x nur die Veränderung von einer Einheit an. Wie groß diese ist, müssen wir uns überlegen, um darauf basierend die nötige Stichprobengröße zu bestimmen, s.u..

Um den Effekt in Sekunden zu berechnen, muss x mit dem erwarteten Effekt multipliziert werden: Bei einem erwarteten Effekt von 5 Sekunden: b * x * 5

Plausibilität von a

# mittlere Überquerungsgeschwindigkeit ist 'mean' und variiert um 'sd'
sample_a <- rnorm( N , mean=12 , sd=.4 )

ggplot() +
    aes(x=sample_a) +
    stat_density(
        adjust = .5,
        fill = random_color(),
        alpha = .5) +
    geom_vline(xintercept = mean(sample_a), color = random_color()) +
    theme_minimal()

Frage 3: Sind die Werte für a plausibel, und was bedeutet die Linie?

Plausibilität von b

# mittlerer Effekt 'meanlog', variiert um 'sd_log'
sample_b <- rlnorm( n=N , meanlog=0 , sdlog=.02 )

ggplot() +
    stat_density(
        aes(x=sample_b),
        adjust = .5,
        fill = random_color(),
        alpha = .5) +
    geom_vline(xintercept = mean(sample_b), color = random_color())+
    theme_minimal()

Frage 4: Sind die Werte für b plausibel?

Plausibilität von mu

# Unabhängige Variable
# x = starren beschleunigt in Sekunden
x_0 = 0  # kein starren: kein Effekt
x_1 = 1 # starren: eine Einheit schneller, als sonst

# x_bar = mean(c(x_0, x_1))
# x_bar

mu_1 <- sample_a + sample_b * (x_1 * effect_size) # den Effekt in Sekunden umrechnen
mu_0 <- sample_a + sample_b * (x_0 * effect_size)

avg_mu_sim <- tibble(sample_a, sample_b, mu_1, mu_0)

ggplot(data=avg_mu_sim) +
    stat_density(aes(x=mu_0), fill = random_color(), alpha = .5) +
    geom_vline(xintercept = mean(mu_0), color = random_color()) +
    stat_density(aes(x=mu_1), fill = random_color(), alpha = .5) +
    geom_vline(xintercept = mean(mu_1), color = random_color())+
    theme_minimal()

Plausibilität von sigma

# wird erst später benötigt
sample_sigma <- runif( n=N , min=0 , max=2) # wird erst später relevant

ggplot() +
    stat_density(
        aes(x=sample_sigma),
        adjust = .5,
        fill = random_color(),
        alpha = .5) +
    geom_vline(xintercept = mean(sample_sigma), color = random_color())+
    theme_minimal()

folgt

Parameter für die Simulation

# Anzahl der Simulationen
nsims <- 10000 # hat nur was mit der anzunehmenden Genauigkeit der Power zutun

# 1. Stichprobengröße
nsubjects <- 40

# 2. Signifikanzniveau
type_I_error <- .01 # alpha

# 4. Effektgröße
effect_size <- effect_size # oben schon festgelegt

• H₁: Menschen benötigen -1 Sekunden weniger Zeit die Straße zu überqueren, wenn sie angestarrt werden), wir legen fest, dass wir in 1% der Fälle H₁ fälschlicherweise akzeptieren und wir legen fest, dass wir 40 Versuchspersonen erheben.
• nsims: Für ein Experiment liegt diese quasi bei 1.
  • Daher wollen wir natürlich eine möglichst hohe Power und strenges Signifikanzniveau haben. Wir möchten aber natürlich nicht so viele Versuchspersonen erheben und dabei muss die Effektstärke auch noch realistisch sein.

Simulation

Durch die Simulation können wir nun rumprobieren und entscheiden, was das für unsere Stichprobengröße bedeutet.

# Ein paar Dinge, die wir später benötigen
df_sim_data <- tibble() # Einen Datensatz für alle simulierten Daten
p_values <- tibble() # Einen Datensatz für alle p-Werte
starring <- c(x_0, x_1) # Eine Manipulation, s.o.

for (sim in 1:nsims){

    time_crossing <-
        sapply(
          starring ,
          function(stare)

            rnorm(

                # Wie viele Versuchspersonen sollen gezogen werden?
                n = nsubjects/2,
                # Das lineare Model (s.o.), welches die Zeit abhängig vom Starren vorhersagt
                mean = sample_a + sample_b * stare * effect_size ,
                # Die Varianz, wie unterschiedlich lange Versuchspersonen brauchen
                sd = sample_sigma ) )

    sim_data <-
        tibble(
          sim_n = sim,
          no_stare=time_crossing[,1],
          stare=time_crossing[,2])

    # t-test rechnen
    sim_ttest <-
      t.test(
        sim_data$stare,
        sim_data$no_stare,
        conf.level = 1 - type_I_error, # zu 1 - Fehlerrate wird H_1 korrekt angenommen
        var.equal = TRUE, # per Definition (s.o.)
        paired = FALSE) # keine Messwiederholung)

    # # oder als lineares model
    # sim_data_lm <-  gather(
    #   sim_data, key = "stare", value = "crossing_time", stare:no_stare)
    #
    # sim_lm <-
    #   summary(
    #     lm(crossing_time ~ 1 + stare, data = sim_data_lm))

    # p-werte sammeln
    sim_p_value <-
        tibble(
          sim_n =sim,
          # p_value = sim_lm$coefficients["starestare", "Pr(>|t|)"],
          p_value = sim_ttest$p.value, # nur den p-Wert behalten
          t_value = sim_ttest$statistic,
          critical_t_value = qt(1-type_I_error/2, sim_ttest$parameter))

    df_sim_data <- rbind(df_sim_data, sim_data)
    p_values <- rbind(p_values, sim_p_value)

    }

Ergebnis

In dem Model legen wir einen Unterschied zwischen angestarrt und nicht angestarrt auf -1 Sekunden fest, H₁ ist in unserem Model wahr.

• Aufgrund von (zufälligen) Schwankungen (bspw. durch Messfehler), wird der kritische p-Wert (0.01), nicht immer unterboten.
  • Genauso, wie der p-Wert auch manchmal unterboten wird, obwohl H₀ wahr ist, Fehler 1. Art!
  • Die Wahrscheinlichkeit versuchen wir anhand der Simulation herauszufinden.

p_values$below_alpha <- ifelse(p_values$p_value <= type_I_error, 1, 0)

type_II_error <- sum(p_values$below_alpha)/nrow(p_values)

ggplot(data = p_values, aes(x = p_value)) +
    stat_density(
          adjust = .5,
          fill = random_color(),
          alpha = .5) + # anderes alpha!
    geom_vline(
        xintercept = type_I_error,
        color = random_color())+
    theme_minimal()

# # 1000 simulierte Studien
# df_sim_data <-
#     df_sim_data %>%
#     group_by(sim_n) %>%
#     summarise(
#         mean_stare = mean(stare),
#         # sd_stare = sd(stare),
#         mean_no_stare = mean(no_stare),
#         # sd_no_stare = sd(no_stare),
#         effect_size = mean_stare - mean_no_stare) %>%
#     cbind(p_values)
#
# DT::datatable(df_sim_data)