library(tidyverse)
library(broom)
library(kableExtra) 
library(ggthemes)
library(ggsci)
library(scales)

d <- read.table(file = 'a1-a15.txt', header = T, sep = "", stringsAsFactors = F, dec = ',') %>% 
    as_tibble

Til að meta fjölda þátttakenda notum við úrtaksmeðaltal og -staðalfrávik á muni í munnmagni milli 20 min og 0 min frá pilot stúdíunni. Við gerum líka ráð fyrir að þið viljið framkvæma þrjú t-próf, eftir 5 min, 10 min og 20 min. Þess vegna framkvæmum við Bonferroni leiðréttingu sem útvíkkar öryggisbilin til að niðurstöður skekkist ekki vegna fjölda tilgátuprófana.

Gefum okkur að næsta tilraun muni hafa sama meðaltal og sama staðalfrávik. Þá gefa einhliða neðri 95% mörk okkur að meðaltalið liggi að öllu jöfnu á bilinu \([\hat\mu, \infty]\) þar sem \(\hat\mu\) er y-gildið á myndritinu fyrir neðan fyrir hvern þátttakendafjölda. Bæst væri ef þátttakendur væru fleiri en 40.

my_fun <- function(n) {
    0.01482 - qt(1 - 0.05 / 3, n - 1) * 0.03792 / sqrt(n)
}

data_frame(n = seq(20, 60)) %>% 
    mutate(lower = my_fun(n)) %>% 
    ggplot(aes(n, lower)) +
    geom_hline(yintercept = 0, lty = 2, alpha = 0.3) +
    geom_line() +
    theme_bw() +
    labs(x = "Fjöldi þátttakenda", y = "Neðri mörk öryggisbils",
         title = "Neðri 95% öryggismörk í einhliðahliða t-prófi með bonferroni leiðréttingu") +
    ggsave("t_conf.png", width = 10, height = 6)

Næsta tilraun átti samt að vera nákvæmari svo við getum skoðað þátttakendafjölda með mismunandi staðalfrávikum. Við sjáum að það hefur töluverð áhrif á nákvæmnina. Besta leiðin til að fá sem nákvæmast mat er því að minnka staðalfrávikið í mælingum. Myndritið bendir til að við viljum >40 þátttakendur til að eiga góðar líkur á að hafa 0 utan öryggismarka.

Við skulum þó alltaf hafa í huga að þetta séu vænt neðri mörk og því ekki víst að 0 verði utan öryggismarka þó svo að við höfum nægan þátttakendafjölda.

my_fun <- function(n, s) {
    0.01482 - qt(1 - 0.05 / 3, n - 1) * s / sqrt(n)
}

crossing(n = seq(15, 60),
         s = seq(0.03, 0.04, 0.002)) %>% 
    mutate(lower = map2_dbl(n, s, my_fun)) %>% 
    ggplot(aes(n, lower, col = factor(s))) +
    geom_hline(yintercept = 0, lty = 2, alpha = 0.3) +
    geom_line() +
    theme_bw() +
    scale_color_jama() +
    scale_x_continuous(breaks = scales::pretty_breaks(8)) +
    labs(x = "Fjöldi þátttakenda", 
         y = "Vænt neðri mörk öryggisbils",
         col = "Staðalfrávik í \nnýjum mælingum",
         title = "Neðri 95% öryggismörk í einhliðahliða t-prófi með bonferroni leiðréttingu",
         subtitle = "Fyrri staðalfrávik var uþb 0.038") +
    ggsave("t_conf_more_precision.png", width = 10, height = 6)

Reiknum nú powerið fyrir prófin. Power segir okkur hlutfallslega hversu oft við munum geta hafnað núlltilgátunni að magnið í munni sé 0 eftir 20 min. Myndin segir það sama og hinar:

Ef við viljum eiga góðar líkur á að fá 0 utan öryggismarka þarf >40 þátttakendur. Oft er miðað við 80% power í vísindalegum rannsóknum og samsvarar það í kringum 41 þátttakanda ef staðalfrávikið í mælingum verður það sama og í pilot stúdíunni

my_power <- function(n, s) {
    1 - pt(qt(1 - 0.05, n - 1) - (0.01482 - 0) / (s / sqrt(n)), n - 1)
}

crossing(n = seq(15, 60),
         s = seq(0.03, 0.04, 0.002)) %>% 
    mutate(power = map2_dbl(n, s, my_power)) %>% 
    ggplot(aes(n, power, col = factor(s))) +
    geom_line() +
    theme_bw() +
    scale_color_jama() +
    scale_x_continuous(breaks = scales::pretty_breaks(8)) +
    scale_y_continuous(labels = percent) +
    labs(x = "Fjöldi þátttakenda", y = "Power",
         col = "Staðalfrávik í \nnýjum mælingum",
         title = "Líkur á að hafna u = 0 fyrir mismunandi n og staðalfrávik") +
    ggsave("est_power.png", width = 10, height = 6)