Power analýza

1 Balíčky

Pro power analýzu budeme používat dva nové balíčky: WebPower a paramtest.

library(WebPower)
library(paramtest)
library(tidyverse)

2 Statistická síla

Statistická síla je vlastně podmíněná pravděpodobnost toho, že pokud alternativní hypotéza platí (čili nulová hypotéza neplatí), a v populaci tedy existuje takový efekt jako očekáváme, podaří se nám zamítnout nulovou hypotézu. Je doplňkem k chybě II. typu (\(\beta\)) a proto se značí 1 – \(\beta\). Vzrůstá s velikostí vzorku, velikostí účinku v populaci a liberálnější hladinou \(\alpha\). Je důležitá především proto, že studie s malou statistickou silou nejsou příliš užitečné. Malá statistická síla znamená, že je velmi nepravděpodobné, že se nám podaří zamítnou nulovou hypotézu, i když v populaci existuje nezanedbatelný účinek. V důsledku toho jsou výsledky “underpowered” výzkumu málo informativní: nevíme, jestli nesignifikantní výsledek je dán jen nedostatečnou statistickou silou, nebo tím, že zkoumaný efekt je v populaci zanedbatelný (blízký nule). Pokud máme naopak velký vzorek, a tedy vysokou statistickou sílu, je velmi nepravděpodobné, že výsledek statistického testu bude nesignifikantní, pokud v populaci existuje nezanedbatelný účinek. Pozn.: “Účinkem” může být míněno mnoho věcí, např. rozdíl mezi průměry, podíl vysvětleného rozptylu, velikost korelace, regresního koeficientu apod.

V zásadě při power analýze operujeme se čtyřmi “veličinami”:

1. hladinou alfa, což je předem stanovená pravděpodobnost chyby I. typu;
2. velikostí účinku;
3. velikostí vzorku;
4. statistickou silou.

Platí přitom, že pokud známe tři z těchto “veličin”, můžeme dopočíst hodnotu zbýbající. Proto existuje více typů power analýzy.

1. Apriorní power analysis. To je ta nejdůležitější. Na základě stanovené hladiny \(\alpha\), očekávané velikosti účinku a požadované statistické síly chceme dopředu stanovit potřebnou velikost vzorku.
2. Post-hoc power analysis. Zde chceme dodatečně zjistit pozorovanou statistickou sílu na základě velikosti vzorku, hladiny \(\alpha\) a pozorované velikosti účinku. Obvykle to není moc užitečné, protože pozorovaný účinek se může více či méně lišit od “skutečného” (populačního) účinku.
3. Sensitivity power analysis. Na základě hladiny \(\alpha\), velisti vzorku a požadované statistické síly chceme vypočíst potřebný účinek. Toto není moc používané, ale může se hodit např. tehdy, když víme, že můžeme získat maximálně pouze určitý počet respondentů, nemáme tušení o velikosti účinku, ale víme, s jakou hladinou \(\alpha\) chceme pracovat a jaké statistické síly chceme dosáhnout. Na základě toho pak můžeme vypočíst požadovaný účinek a posoudit, jestli je se jedná o realistickou hodnotu (jestli tak velký očekávaný efekt v populaci je realisitký).
4. Criterion power analysis. Zde “dopočítáváme” hladinu \(\alpha\) na základě stanovené velikosti vzorku, očekávané velikosti účinku a požadované statistické síly.
5. Compromise power analysis. Zde jako vstupní hodnoty vkládáme očekávanou velikost účinku, poměr mezi pravděpodobností chyby I. a II. typu (\(\alpha\)) a \(\beta\)) a velikost vzorku. Zjistíme, s jakou hladinou alfa máme pracovat (a jaké je implikovaná statistická síla), abychom zachovali stanovený poměr mezi pravděpodobnostmi obou typů chyb. Hodí se, když víme, jak “nákladná” by byla chyba I. a II. typu a na základě toho chceme stanovit hladinu \(\alpha\). Pokud by např. chyba I. typu obnášela čtyřnásobně vyšší náklady než chyba II. typu, stanovili bychom poměr \(\beta\)/\(\alpha\) na hodnotu 4 či více.

3 WebPower

Balíček WebPower je takovou obdobou programu GPower v tom smyslu, že umožňuje odhad statistické síly pro základní, ale nejčastěji používané statistické testy/modely. Má mnoho funkcí v závislosti na statistickém testu/modelu, pro který chceme provést power analýzu, a proto se podíváme jen na některé z nich. Přehled všech funckí můžete najít zde.

3.1 Korelace

Nejdříve si vyzkoušíme odhad statistické síly pro korelace Dejme tomu, že zkoumáme vztah mezi stresem a well-beingem Na základě předchozích výzkumů očekáváme korelaci r = 0.3. Můžeme si položit otázku: kdybychom získali měřili stres a well-being u 50 osob, jakou statistickou sílu bychom měli při alfa = 0.05?

Můžeme k tomu použít funkci wp.correlation(). Jeden z klíčových argumentů (power, n, r, alpha) musíme nechat prázdný, protože jsme neuvedli power, vypočte se právě ona

wp.correlation(n = 50, r = 0.3, alpha = .05)

># Power for correlation
># 
>#      n   r alpha     power
>#     50 0.3  0.05 0.5728731
># 
># URL: http://psychstat.org/correlation

Defaultně se počítá s oboustranným testem, ale můžeme použít argument alternative pro stanovení “směru” alternativní hypotézy Také můžeme změnit nulovou hypotézu pomocí argumentu rho0 (defalt je rho0 = 0)

Kdybychom chtěli vypočíst statistckou sílu pro různé velikosti vzorku, můžeme postupovat takto. Funkce balíčku webpower() umožňují, aby argument n mohl zahrnovat více hodnot, takže n může být číselný vektor s více hodnotami. Zde vidíme, že při očekávané korelaci 0,30 v populaci a hladině alfa = 0,05 by velikost vzorku 10, 20 ani 30 nestačila k dosažení dostatečné statistické síly.

wp.correlation(n = c(10, 20, 30), # Jeden z argumentů může zahrnovat více hodnot
               r = 0.3, alpha = .05)

># Power for correlation
># 
>#      n   r alpha     power
>#     10 0.3  0.05 0.1377062
>#     20 0.3  0.05 0.2582591
>#     30 0.3  0.05 0.3731596
># 
># URL: http://psychstat.org/correlation

Když uvedeme power, ale nikoli n, vypočte se potřebná velikost vzorku.

wp.correlation(power = .8, r = 0.3, alpha = .05)

># Power for correlation
># 
>#            n   r alpha power
>#     83.94932 0.3  0.05   0.8
># 
># URL: http://psychstat.org/correlation

Často si nemusíme být uplně jisti ohledně očekávaného efektu. V případě korelací tedy můžeme vypočíst potřebnou velikost vzorku pro různé hodnoty očekávané korelace. Můžeme si pomoci funkcí seq(), která dovede vygenerovat řadu čísel od – do po stanovených krocích.

# Řada čísel od 0.1 do 0.3 (včetně) po krocích velikosti 0.01
rs <- seq(0.1, 0.3, by = 0.01)
rs

>#  [1] 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23 0.24
># [16] 0.25 0.26 0.27 0.28 0.29 0.30

Tak jsme si vytvořili vektor s hodnotami korelací od 0,1 do 0,3. Ten teď iterativně vložíme do funkce wp.correlation() pomocí map() a výsledného bojektu si pak extrahujeme všechny relevantní statistiky do jednoho dataframu pomocí map_df()

rs %>% 
  map(~wp.correlation(power = .8, 
                      r = .x, 
                      alpha = .05)) %>% 
  map_df(~.x[c("r","alpha","power", "n")])

># # A tibble: 21 × 4
>#        r alpha power     n
>#    <dbl> <dbl> <dbl> <dbl>
>#  1  0.1   0.05   0.8  782.
>#  2  0.11  0.05   0.8  645.
>#  3  0.12  0.05   0.8  542.
>#  4  0.13  0.05   0.8  461.
>#  5  0.14  0.05   0.8  397.
>#  6  0.15  0.05   0.8  346.
>#  7  0.16  0.05   0.8  303.
>#  8  0.17  0.05   0.8  268.
>#  9  0.18  0.05   0.8  239.
># 10  0.19  0.05   0.8  214.
># # … with 11 more rows

Co když chceme zjistit power pro různé velikosti vzorku, různé r a různou alfa? Vytvoříme si list relevantních hodnot

input_list <- list(
  n = seq(10, 100, by = 1), # velkost vzorku od 10 do 100
  r = seq(0.1, 0.6, by = 0.1), # velikost korelace od 0.1 po 0.6
  alpha = c(0.05, 0.01) # hladina alfa
)
input_list

># $n
>#  [1]  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28
># [20]  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47
># [39]  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65  66
># [58]  67  68  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85
># [77]  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100
># 
># $r
># [1] 0.1 0.2 0.3 0.4 0.5 0.6
># 
># $alpha
># [1] 0.05 0.01

Pak použijeme funkcícross_df() a získáme všechny možné kombinace n, r a alpha jako řádky dataframu

input_df <- cross_df(input_list)
input_df

># # A tibble: 1,092 × 3
>#        n     r alpha
>#    <dbl> <dbl> <dbl>
>#  1    10   0.1  0.05
>#  2    11   0.1  0.05
>#  3    12   0.1  0.05
>#  4    13   0.1  0.05
>#  5    14   0.1  0.05
>#  6    15   0.1  0.05
>#  7    16   0.1  0.05
>#  8    17   0.1  0.05
>#  9    18   0.1  0.05
># 10    19   0.1  0.05
># # … with 1,082 more rows

Když máme tento dataframe, který má sloupce pojmenovány přesně podle argumentů funkce, kterou chceme iterativně uplatnit, můžeme jej vložit do funkce pmap() společně s funkcí, do které chceme tyto argumenty iterativně (po řádcích) vložit, a to je wp.correlation()

pwr <- input_df %>% 
  pmap(wp.correlation)

Každý prvek pwr je list, který obsahuje tyto prvky:

str(pwr[[1]])

># List of 7
>#  $ n          : num 10
>#  $ r          : num 0.1
>#  $ alpha      : num 0.05
>#  $ power      : num 0.0575
>#  $ alternative: chr "two.sided"
>#  $ method     : chr "Power for correlation"
>#  $ url        : chr "http://psychstat.org/correlation"
>#  - attr(*, "class")= chr "webpower"

Pokud bychom některé z nich chtěli extrahovat do jednoho dataframu, můžeme použít map_df()

pwr <- pwr %>% 
  map_df(~.x[c("r", "alpha", "n", "power")])
pwr

># # A tibble: 1,092 × 4
>#        r alpha     n  power
>#    <dbl> <dbl> <dbl>  <dbl>
>#  1   0.1  0.05    10 0.0575
>#  2   0.1  0.05    11 0.0590
>#  3   0.1  0.05    12 0.0605
>#  4   0.1  0.05    13 0.0619
>#  5   0.1  0.05    14 0.0632
>#  6   0.1  0.05    15 0.0645
>#  7   0.1  0.05    16 0.0658
>#  8   0.1  0.05    17 0.0670
>#  9   0.1  0.05    18 0.0682
># 10   0.1  0.05    19 0.0695
># # … with 1,082 more rows

Výsledný dataframe pak můžeme použít pro tvorbu grafu power curves. Graf níže např. ukazuje, jak roste statistická síla (na ose Y) v závislosti na velikosti vzorku (na ose X) pro různou velikost účinku v populaci (různou korelaci v populaci, vyznačena barevně) a alfa (0.01 a 0.05).

pwr %>% 
  ggplot(aes(n, power, color = as.factor(r))) +
  geom_line(size = 1) +
  facet_wrap(~alpha, labeller = "label_both") +
  geom_hline(yintercept = .8, linetype = 2) +
  labs(color = "r") + 
  theme(legend.position = "top")

Abychom si usnadnili práci, můžeme si vytvoři vlastní funkci, která budou stačit dva argumenty (plus argument dot-dot-dot pro doplnění jakýkoli dalších argumentů do fn, které jsou konstantní a nepotřebujeme je tedy umisťovat do input_list):

1. input: objekt se vstupními argumenty, které chceme vložit do funkce fn). Pokud je input není typu data.frame, vytvoří se z něj dataframe se všemi kombinacemi prvků,
2. fn: některá z funkcí balíčku WebPower.
3. dot-dot-dot: pro dodatečné argument(y), které bychom chtěli také vložit do fn, ale nechce se nám je stanovovat v rámci input_list, protože jsou konstatní.

pwr_analysis <- function(input, fn, ...) {
  
  if (!is.data.frame(input)) {
    input <- input %>% 
      purrr::cross_df()
  }
  
  out_list <- pmap(input, fn, ...)
  
  out_df <- out_list %>% 
    purrr::map_df(~.x[names(.x)])
  
  out_df
}

S touto funkcí pak stačí vytvořit list se vstupními argumenty. Pokud chceme vypočíst power pro statistiský test korelačního koeficientu v závislosti na velikost vzorku, velikosti korelace a hladiny alfa, vytvořili bychom si tento list a jako fn uvedli funkci wp.correlation

pwr_output <- pwr_analysis(
  list(n = seq(10, 100, by = 1), # velkost vzorku od 10 do 100
       r = seq(0.1, 0.6, by = 0.1), # velikost korelace od 0.1 po 0.6
       alpha = c(0.05, 0.01)), # hladina alfa
  fn = wp.correlation
)
pwr_output

># # A tibble: 1,092 × 7
>#        n     r alpha  power alternative method                url               
>#    <dbl> <dbl> <dbl>  <dbl> <chr>       <chr>                 <chr>             
>#  1    10   0.1  0.05 0.0575 two.sided   Power for correlation http://psychstat.…
>#  2    11   0.1  0.05 0.0590 two.sided   Power for correlation http://psychstat.…
>#  3    12   0.1  0.05 0.0605 two.sided   Power for correlation http://psychstat.…
>#  4    13   0.1  0.05 0.0619 two.sided   Power for correlation http://psychstat.…
>#  5    14   0.1  0.05 0.0632 two.sided   Power for correlation http://psychstat.…
>#  6    15   0.1  0.05 0.0645 two.sided   Power for correlation http://psychstat.…
>#  7    16   0.1  0.05 0.0658 two.sided   Power for correlation http://psychstat.…
>#  8    17   0.1  0.05 0.0670 two.sided   Power for correlation http://psychstat.…
>#  9    18   0.1  0.05 0.0682 two.sided   Power for correlation http://psychstat.…
># 10    19   0.1  0.05 0.0695 two.sided   Power for correlation http://psychstat.…
># # … with 1,082 more rows

Podobně si pro usnadnění můžeme vytvořit funkci pro vytváření grafu zachycujícího power curves. Tato funkce bude mít tyto tři povinné argumenty:

• data: použitý vstupní dataset,
• x: proměnná na ose X (defaultně "n" jako velikost vzorku),
• y: proměnná na ose Y (defaultně "power" čili statistická síla).

A tyto nepovinné argumenty:

• color: proměnná mapovaná na barvu.
• facet_wrap: proměnná pro vytvoření fazet.
• facet_grid: dvě proměnné (řádky, sloupce) pro vytvoření fazet.
• hline: zobrazit horizontální čáru se zadaným průsečíkem v ose Y.
• vline: zobrazit vertikální čáru se zadaným průsečíkem v ose X.

pwr_plot  <- function(data, x = "n", y = "power", 
                      color = NULL, facet_wrap = NULL, facet_grid = NULL,
                      hline = NULL, vline = NULL) {
  
  # Převést sloupce color na factor, pokud ještě není factor
  if (!is.null(color) && !is.factor(data[[color]])) {
    data[[color]] <- factor(data[[color]])
  }
  
  # Vytvořit základní kostru grafu
  p <- data %>% 
    ggplot(aes_string(x = x, y = y)) +
    theme(legend.position = "top") +
    expand_limits(y = 0)
  
  # Přidat křivku. Pokud je color vyplněno, přidat křivky různé barvy
  if (is.null(color)) {
    p <- p + geom_line(aes(group = 1))
  } else {
    p <- p + 
      geom_line(aes_string(color = color))
  }
  
  # Přidat horizontální nebo vertikální čáru
  if (!is.null(hline)) {
    p <- p + geom_hline(yintercept = hline,
                        linetype = 2)
  }
  
  if (!is.null(vline)) {
    p <- p + geom_hline(yintercept = vline,
                        linetype = 2)
  }
  
  # Přidat fazety, pokud jsou vyplněny
  if (!is.null(facet_wrap)) {
    p <- p + facet_wrap(vars(.data[[facet_wrap]]), 
                    labeller = "label_both")
  } else if (!is.null(facet_grid)) {
    p <- p + facet_grid(rows = vars(.data[[facet_grid[1]]]),
                    cols = vars(.data[[facet_grid[2]]]),
                    labeller = "label_both")
  }
  
  return(p)
}

A rovnou si tuto funkci můžeme vyzkoušet:

pwr_plot(pwr_output,
         x = "n", y = "power",
         color = "r", facet_wrap = "alpha")

Cvičení

Spočítejte statistickou sílu pro oboustranné testy tří korelačních koeficientů r = 0,1, 0,2 a 0,3 při hladině alfa = 0,05 a velikosti vzorku od 30 do 1000 a vytvořte vhodné grafy ukazující, jak statistická síla těchto tří korelačních koeficientů závisí na velikosti vzorku. Klidně k tomu použijte funkce pwr_analysis() a pwr_plot().

4 t-testy

4.1 Párový t-test

Dejme tomu, že chceme testovat efekt nějaké intervence; máme skupinu studentů a měříme je před a po intervenci (čili párová data). Očekáváme velikost účinku d = 0,5. Kolik účastníků potřebujeme k dosažení statistické síly 0,8? Můžeme k tomu použít funkci wp.t()

wp.t(d = .5, power = 0.8, alpha = 0.05, type = "paired")

># Paired t-test
># 
>#            n   d alpha power
>#     33.36713 0.5  0.05   0.8
># 
># NOTE: n is number of *pairs*
># URL: http://psychstat.org/ttest

Pro dosažení 80% statistické síly ba nám stačilo získat 34 osob.

Cvičení

Spočítejte, jako velikost vzorku bychom potřebovali pro dosažení aspoň 80% statistické síly u oboustranného párového t-testu s hladinou alfa = 0,05, ale pro různé velikostí účinku Cohenovo d = 0,1 až 0,8 (po krocích o velikosti 0,01). Vhodně výsledky vizualizujte.

4.2 Nezávislý t-test

Dejme tomu, že místo vnitrosubjektového použijeme mezisubjektový desing, takže máme dvě nezávislé skupiny, jedna bude vystavena placebo-intervenci a druhá skutečné intervenci. Očekáváme menší velikost účinku (d = 0,3) a chceme mít vyvážený počet osob v obou skupinách

wp.t(d = .3, power = 0.8, alpha = 0.05, type = 'two.sample')

># Two-sample t-test
># 
>#            n   d alpha power
>#     175.3847 0.3  0.05   0.8
># 
># NOTE: n is number in *each* group
># URL: http://psychstat.org/ttest

Co kdyby ale poměr skupin nebyl vyvážený. Např. kdybychom chtěli poskytnout skutečnou intervenci co nejvíce lidem, ale věděli bychom, že celková velikost vzorku, kterou můžeme získat, je 500 osob. Stále bychom ale chtěli mít sílu testu aspoň 0,8. Jak moc mohou být skupiny nevyvážené?

Začneme tím, že se vytvoříme tibble se různými kombinacemi velikostí skupin, které dohromady dávají 500. Sloupce musíme pojmenovat podle názvů argumentů ve funkci wp.t(), které chceme měnit (pro nápovědu s přehledem argumentů se nebojte použít ?wp.t() nebo help(wp.t))

input <- tibble(
  n1 = 2:498,
  n2 = 500 - n1,
)

Ten pak vložíme do funkce pwr_analysis() a přidáme jako další argumenty velikost účinku d, hladinu alfa a typ t-testu ("two.sample.2n" znamená dvouvýběrový t-tests s různě velkými skupinami). Statistickou sílu jsme neuvedli jako argument, protože tu nám to má dopočíst.

output <- input %>% 
  pwr_analysis(fn = wp.t, 
               d = .3, alpha = 0.05, type = "two.sample.2n")
output

># # A tibble: 497 × 9
>#       n1    n2     d alpha  power alternative note                  method url  
>#    <int> <dbl> <dbl> <dbl>  <dbl> <chr>       <chr>                 <chr>  <chr>
>#  1     2   498   0.3  0.05 0.0707 two.sided   n1 and n2 are number… Unbal… http…
>#  2     3   497   0.3  0.05 0.0811 two.sided   n1 and n2 are number… Unbal… http…
>#  3     4   496   0.3  0.05 0.0916 two.sided   n1 and n2 are number… Unbal… http…
>#  4     5   495   0.3  0.05 0.102  two.sided   n1 and n2 are number… Unbal… http…
>#  5     6   494   0.3  0.05 0.113  two.sided   n1 and n2 are number… Unbal… http…
>#  6     7   493   0.3  0.05 0.123  two.sided   n1 and n2 are number… Unbal… http…
>#  7     8   492   0.3  0.05 0.134  two.sided   n1 and n2 are number… Unbal… http…
>#  8     9   491   0.3  0.05 0.145  two.sided   n1 and n2 are number… Unbal… http…
>#  9    10   490   0.3  0.05 0.155  two.sided   n1 and n2 are number… Unbal… http…
># 10    11   489   0.3  0.05 0.166  two.sided   n1 and n2 are number… Unbal… http…
># # … with 487 more rows

Nakonec si vytvoříme graf s výsledky.

output %>% 
  pwr_plot(x = "n1", y = "power",
           hline = 0.8)

Jak můžeme vidět, pro dosažení alespoň 80% statistické síly při celkové velikost vzorku 500 osob by maximální nepoměr mezi velikostmi skupin mohl být přibližně 125 : 375. Pokud bychom chtěli přesné číslo:

output %>% 
  filter(power >= .80) %>% 
  slice_min(n1)

># # A tibble: 1 × 9
>#      n1    n2     d alpha power alternative note                    method url  
>#   <int> <dbl> <dbl> <dbl> <dbl> <chr>       <chr>                   <chr>  <chr>
># 1   114   386   0.3  0.05 0.802 two.sided   n1 and n2 are number i… Unbal… http…

Cvičení

Spočítejte, jaké statistické síly pro oboustranný nezávislý t-test s hladinou alfa = 0,05 bychom dosáhli při velikosti vzorku od 5 do 500 osob v každé skupině (předpokládejte početně vyrovnané skupiny) a pro velikosti účinku Cohenovo d = 0,1, 0,2, 0,3 atd. až 0,8. Výsledky vhodně vizualizujte.

5 ANOVA

5.1 One-Way ANOVA

Dejme tomu, že chceme srovnat čtyři skupiny a očekáváme velikost účinku \(R^2 = 0,2\). Kdybychom mohli získat pouze pro každou skupinu pouze 10 účastníků (tj. dohromady 40 osob), jakou statistickou sílu by měl F-test ověřující shodu průměrů skupin.

Nejdříve vypočteme Cohenovo f, protože to je velikost účinku, která se používá v rámci power analýzy pro analýzu rozptylu. Je to vlastně odmocnina z poměru mezi přírůstkem vysvětleného rozptylu (po zahrnutí daného efektu do modelu) vůči zbývajícímu nevysvětlenému rozptylu. K jejímu výpočtu z R² si můžeme rovnou vytvořit vlastní funkci: (r2_full - r2_null je onen přírůstek).

r2_to_f <- function(r2_full, # Vysvětlený rozptyl po přidání efektu/efektů
                    r2_null = 0, # Původní vysvětlený rozptyl
                    squared = FALSE) { # Chceme čtverec Cohenova F?
  
  if (squared == TRUE) {
    (r2_full - r2_null) / (1 - r2_full)
  } else {
    sqrt((r2_full - r2_null) / (1 - r2_full))
  }
}

cohen_f <- r2_to_f(r2_full = .2)
cohen_f

># [1] 0.5

Pak všechny vstupní parametry vložíme do funkce wp.anova():

wp.anova(k = 4, # Počet skupin
         n = 40, # Celková velikost vzorku
         alpha = .05, # Hladina alfa
         f = cohen_f) # Cohenovo f

># Power for One-way ANOVA
># 
>#     k  n   f alpha     power
>#     4 40 0.5  0.05 0.7096896
># 
># NOTE: n is the total sample size (overall)
># URL: http://psychstat.org/anova

Jak můžeme vidět, při těchto vstupních parametrech bychom dosáhli přibližně 71% statistické síly.

Také můžeme vypočíst minimální velikost účinku pro požadovanou statistickou sílu, velikost vzorku a hladinu alfa:

wp.anova(f = NULL, k = 4, n = 100, power = 0.8, alpha = 0.05)

># Power for One-way ANOVA
># 
>#     k   n         f alpha power
>#     4 100 0.3369881  0.05   0.8
># 
># NOTE: n is the total sample size (overall)
># URL: http://psychstat.org/anova

Protože podíl vysvětleného roztptylu je zřejmě srozumitelnější než Cohenovo f, můžeme si vytvořit funkci pro převod f na R²

f_to_r2 <- function(f) {
  f^2 / (1 + f^2)
}

A použít ji na předchozí výsledek

f <- wp.anova(f = NULL, k = 4, n = 100, power = 0.8, alpha = 0.05)$f
f_to_r2(f)

># [1] 0.10198

Kdychom tedy srovnávali čtyři skupiny po 25 lidech (celkový vzorek by tedy zahrnoval 100 osob), pracovali s hladinou alfa = 0,05 a chtěli dosáhnout alespoň 80% statistické síly, musela by velikost účinku v populaci být minimálně okolo R² = 0,10.

Stejně tak bychom mohli chtít zjistit požadovanou velikost vzorku pro danou velikost účinku, alfu a počet skupin:

wp.anova(f = cohen_f, k = 4, n = NULL, power = 0.8, alpha = 0.05)

># Power for One-way ANOVA
># 
>#     k        n   f alpha power
>#     4 47.70445 0.5  0.05   0.8
># 
># NOTE: n is the total sample size (overall)
># URL: http://psychstat.org/anova

Jak můžeme vidět, v tomto případě by to bylo 48 osob (když zaokrouhlíme n nahoru na celé číslo)

Funkce wp.kanova() je zobecněním pro více nezávislých proměnných (k-way ANOVA). Má tyto argumenty (ten z nich, který chceme dopočítat, necháváme opět prázdný):

• n: Velikost vzorku.
• power: Požadovaná statistická síla.
• ndf: stupně volnosti pro daná efekt. Pro hlavní efekty to bude počet úrovní daného faktoru mínus jedna. Pro interakční efekty to bude násobek. Kdyby měl např. jeden faktor tři úrovně a druhý faktor čtyři úrovně a chtěli bychom odhadnout interakci mezi nimi, stupně volnosti pro tento interakční efekt by byly (3 – 1) * (4 – 1) = 6.
• ng: Celkový počet skupin (kombinací všech úrovní faktorů).
• f: velikost účinku Cohenovo f.
• alpha: Hladina alfa, se kterou chceme pracovat.

Dejme tomu, že chceme udělat power analýzu pro interakční efekt dvou faktorů. Jeden z nich má tři úrovně a druhý z nich čtyři, takže stupně volnosti pro tento interakční efekt činí df = 6 a celkový počet skupin je \(3\times 4 = 12\). Na základě předchozímu výzkumu očekáváme velikost účinku parciální éta na druhou \(\eta_p = 0.3\), což je poměr mezi rozptylem vysvětleným daným efektem vůči nevysvětlenému rozptylu v modelu bez daného efektu. Opět si tuto velikost účinku převedeme na Cohenovo f pomocí jednoduché funkce a vložíme výsledek do wp.kanova() v argumentu f.

etap2_to_f <- function(etap2) {
  sqrt(etap2 / (1 - etap2))
}

f <- etap2_to_f(0.3)
f 

># [1] 0.6546537

wp.kanova(n = NULL,  
          power = .8, # požadovaná statistická síla
          ndf = 6, # stupně volnosti pro daný efekt
          ng = 12, # celkový počet skupin (násek počtu úrovní jednotlivých faktorů)
          f = cohen_f, # velikost účinku
          alpha = .05) # Hladina alfa

># Multiple way ANOVA analysis
># 
>#            n ndf      ddf   f ng alpha power
>#     61.70124   6 49.70124 0.5 12  0.05   0.8
># 
># NOTE: Sample size is the total sample size
># URL: http://psychstat.org/kanova

Jak můžeme vidět, při zadaných vstupních parametrech bychom potřebovali vzorek o velikosti 62 osob (reálně o něco více, třeba \(6\times 12 = 72\), protože 62 osob by nešlo rovnoměrně rozdělit do 12 skupin)

5.2 Lineární regrese

Výzkumník se domnívá, že průměr známek ze střední školy a výsledek TSP vysvětluje 25 % rozptylu v průměru známek na vysoké škole. Pokud má vzorek 50 osob, jaká je velikost účinku pro tento efekt (pracujeme-li s hladinou alfa = 0.01).

Pro odpoveď na tuto otázu můžeme použít funkci wp.regression(). Jako velikost účinku do ní ale musíme zadat čtverec Cohenova f, ale už máme funkci pro převod R² na Cohenovo f, argumentem squared si jen řekneme o umocnění na druhou.

f2 <- r2_to_f(0.25, squared = TRUE)
f2

># [1] 0.3333333

A pak požadované parametry vložíme do funkce wp.regression()

wp.regression(n = 50, # velikost vzorku
              p1 = 2, # počet prediktorů
              f2 = f2, # Cohenovo f2
              alpha = 0.01) # Hladina alfa

># Power for multiple regression
># 
>#      n p1 p2        f2 alpha    power
>#     50  2  0 0.3333333  0.01 0.840395
># 
># URL: http://psychstat.org/regression

Jiný výzkumník se domnívá, že hodnocení z přijímacího pohovoru vysvětluje dalších 5 % rozptylu známek na VŠ. Jakou velikost vzorku potřebujeme získat pro dosažení 80% statistické síly při hladině alfa = 0.01? Opět si nejdříve vypočteme, jakmu Cohenovu f² daný přírůstek vysvětleného rozptylu odpovídá.

f2 <- r2_to_f(0.3, # Celkový vysvětlený rozptyl
              r2_null = 0.25, # Původní vysvětlený rozptyl
              squared = TRUE) # Chceme čtverec Cohenova f
f2

># [1] 0.07142857

A poté všechny parametry vložíme do funkce wp.regression():

wp.regression(n = NULL, 
              power = .8, # požadovaná statistická síla
              p1 = 3, # Celkový počet prediktorů
              p2 = 2, # Počet prediktorů v původním modelu
              f2 = f2, # Velikost účinku Cohenovo f2
              alpha = 0.01) # Hladina alfa

># Power for multiple regression
># 
>#           n p1 p2         f2 alpha power
>#     166.894  3  2 0.07142857  0.01   0.8
># 
># URL: http://psychstat.org/regression

Jak vidíme, pro dosažení 80% statistické síly je v tomto případě zapotřebí vzorek 167 osob.

6 Balíček paramtest

Nejflexibilnější přístup představuje power analýza pomocí simulace (ale zároveň nejsložitější pro naprogramování). Obecně balíček paramtest umožňuje simulovat výběrové rozdělení jakýchkoli statistik potřebujeme.

Začněme relativně jednuše: výpočtem statistické síly pro nezávislý t-test. Zajímá nás statistická síla při počtu osob 20 a 40 v první a druhé skupině, velikosti účinku Cohenovo d = 0,5 a hladině alfa = 0,05.

Samozřejmě bychom mohli použít přesné, analytické řešení bez simulace např. pomocí wp.t() funkce z WebPower.

wp.t(n1 = 20, n2 = 40, d = 0.5, alpha = .05, type = "two.sample.2n")

># Unbalanced two-sample t-test
># 
>#     n1 n2   d alpha     power
>#     20 40 0.5  0.05 0.4347675
># 
># NOTE: n1 and n2 are number in *each* group
># URL: http://psychstat.org/ttest2n

Ale tento přístup má omezení, protože např. předpokládá shodný rozptyl skupin.

6.1 Nezávislý t-test a funkce run_test()

Balíček paramtest nejprve vyžaduje napsat si vlastní funkci pro generování dat, provedení daného statistického testu a uložení výsledků

První argument by měl být vžy simNum, ten slouží ke stanovení počtu simulací. Další argumenty pak upřesňují, jaká data budeme generovat, a odvíjejí se od toho, které funkce k tomu použijeme a jaký test chceme použít.

Obecně musí vytvořená funkce mít následující strukturu

1. Kód sloužící ke generování dat.
2. Kód sloužící k odhadu modelu (provedení statistického testu) a uložení tohoto modelu (fit-objektu.
3. Kód sloužící k extrakci potřebných informací (statistik) z modelu (fit-objektu. Obvykle se jedná o klíčové koeficienty modelu, testovés statistiky, odhady standardních chyb nebo p-hodnoty.
4. Stanovení výstupních vektoru. Musí jím být atomický numerický vektoru, jehož jsou pojmenovány.

Protože chceme provést t-test, musíme simulovat data pro dvě skupiny. Využijeme k tomu funkci rnorm(), která náhodně generuje n-hodnot z normálního rozdělení s průměrem m a směrodatnou odchylkou sd, proto toto budou jedny ze vstupních argumentů. Další argumenty jsou doplňkové, ale mohou se hodit, pokud bycho chtěli měnit typ nezávislého t-testu pomocí argumentu var.equal a hladinu alfa pomocí argumentu alpha.

t_fun <- function(simNum, # argument sloužící ke stanovení počtu simulací
                  n1, n2, m1, m2, sd1, sd2, # argumenty sloužící pro simulaci vzorků z populace
                  var.equal = FALSE, # argument používaný funkcí t.test
                  alpha = .05) { # arugment pro stanovení hladiny alfa
  
  x1 <- rnorm(n1, mean = m1, sd = sd1) # simulace dat pro první skupinu
  x2 <- rnorm(n2, mean = m2, sd = sd2) # simulace dat pro druhou skupinu
  
  t <- t.test(x1, x2, var.equal = var.equal)  # provedení t-testu
  
  stat <- t$statistic # uložení testové statistiky t
  p <- t$p.value # uložení p-hodnoty
  
  return(c(t = stat, p = p, sig = (p < alpha))) # output funkce
  # výstupní hodnoty funkce t_fun: named vektor s 
  # testovou statistikou, p-hodnotou a logickou hodnotou informující o tom,
  # zda bychom zamítli nulovou hypotézu
}

Pro provedení simulace samotné pak použijeme funkci run_test(). Nejprve uvedeme, kterou funkci využijeme ke generaci dat, provedení statistického testu a uložení potřebných statistik (tj. funkci t_fun), kterou jsme si vytvořili. Pak uvedem počet simulací (ideálně aspoň 1000) pomocí argumentu n.iter, formát výstupu ("data.frame" nebo "list") a nakonec všechny argumenty, které budou využity funkcí t_fun. Pro použití více jader procesoru a urychlení simulace můžeme použít argument parallel = "snow" a dále pomocí ncpus specifikovat počet jader (nejlépe ne všechny, alespoň jedno nechat volné pro ostatní aplikace, jinak se vám mohou sekat ostatní programy)

pwr_ttest <- run_test(
  t_fun, # použitá funkce
  n.iter = 2000, # počet simulací (kolikrát je funkce spuštěna)
  output = "data.frame", # formát výstupu: list nebo data.frame
  parallel = "snow", # Použití více jader procoru
  ncpus = 6, # Kolik jader použít
  n1 = 20, n2 = 40, m1 = 0, m2 = 0.5, sd1 = 0.5, sd2 = 1.5, # argumenty vložené
  var.equal = FALSE, alpha = 0.05) # do funkce t_fun

># Running 2,000 tests...

Výsledky simulace extrahujeme pomocí funkce results() Každý řádek je výsledek jedné simulace

res <- results(pwr_ttest) %>% 
  as_tibble()
res

># # A tibble: 2,000 × 4
>#     iter    t.t        p   sig
>#    <int>  <dbl>    <dbl> <dbl>
>#  1     1 -1.90  0.0631       0
>#  2     2 -0.602 0.550        0
>#  3     3 -1.83  0.0725       0
>#  4     4 -2.05  0.0458       1
>#  5     5 -0.654 0.516        0
>#  6     6 -3.59  0.000735     1
>#  7     7 -2.74  0.00808      1
>#  8     8 -0.463 0.645        0
>#  9     9 -1.80  0.0767       0
># 10    10 -0.189 0.850        0
># # … with 1,990 more rows

Odhad power pak vypočteme jednoduše jako podíl signifikantních výsledků

res_summary <- res %>%
  summarise(power = mean(sig),
            sum_sig = sum(sig))
res_summary

># # A tibble: 1 × 2
>#   power sum_sig
>#   <dbl>   <dbl>
># 1 0.454     907

Kdybychom chtěli vypočíst interval spolehlivosti pro tento odhad, můžeme použít funkci binom.test():

binom.test(x = res_summary$sum_sig, # počet "úspěchů" (signifikantních výsledků) 
           n = nrow(res)) # počet pokusů (celkový počet simulací)

># 
>#  Exact binomial test
># 
># data:  res_summary$sum_sig and nrow(res)
># number of successes = 907, number of trials = 2000, p-value = 3.48e-05
># alternative hypothesis: true probability of success is not equal to 0.5
># 95 percent confidence interval:
>#  0.4315137 0.4756229
># sample estimates:
># probability of success 
>#                 0.4535

6.2 Různé vstupní hodnoty a funkce grid_search()

V předchozím případě jsem použili funkci run_test(), protože jsme jsme použili jenom jedny vstupní hodnoty. Pomocí funkce grid_search jich ale můžeme použít několik, například různé velikosti první (n1 = c(10, 20, 30)) a druhé (n2 = c(40, 50, 60))) skupiny:

pwr_ttest <- grid_search(
  t_fun, # použitá funkce
  params = list(n1 = c(10, 20, 30), n2 = c(40, 50, 60)), # různé velikosti vzorku
  n.iter = 2000, # počet simulací
  output= 'data.frame', # výstupní formát 
  parallel = 'snow', # paralelní zpracování (využití více jader procesoru)
  ncpus = 6, # počet jader procesoru k využití
  m1 = 0, m2 = 0.5, sd1 = 0.5, sd2 = 1.5) # ostatní vstupní argumenty

># Running 18,000 tests...

Výsledky si opět můžeme uložit do objektu res, spočítat si podíl signifikantních výsledků v závislosti na velikosti skupin a vytvořit si graf pomocí funkce pwr_plot()

# Uložení výsledků
res <- results(pwr_ttest)

# Shrnutí výsledků
res_summary <- res %>%
  group_by(n1.test, n2.test) %>% 
  summarise(power = mean(sig))

># `summarise()` has grouped output by 'n1.test'. You can override using the
># `.groups` argument.

res_summary

># # A tibble: 9 × 3
># # Groups:   n1.test [3]
>#   n1.test n2.test power
>#     <dbl>   <dbl> <dbl>
># 1      10      40 0.399
># 2      10      50 0.456
># 3      10      60 0.516
># 4      20      40 0.478
># 5      20      50 0.532
># 6      20      60 0.628
># 7      30      40 0.486
># 8      30      50 0.548
># 9      30      60 0.634

# Graf s výsledky
res_summary %>% 
  pwr_plot(x = "n1.test", y = "power", color = "n2.test")

6.2.1 Cvičení

Zkuste generovat data za předpokladu, že nulová hypotéza platí (obě populace mají shodný průměr rovný např. m1 = m2 = 0). První populce má ale mnohem vyšší směrodatnou odchylku (např. sd1 = 5) než druhá (sd2 = 1). Zkuste, jestli bude očekávaný podíl signifikantních výsledků odpovídat pozorovanému (simulovanému), i když použijeme t-test předpokládající shodu rozptylů (var.equal = TRUE). Protože nulová hypotéza platí, očekávali bychom, že pouze 5% p-hodnot bude menších než 0,05. Zkuste to ověřit pro různý poměr velikostí skupin n1/n2: 10/90, 25/75 a 50/50.

6.3 Simulace na základě korelační matice

Data můžeme simulovat i na základě očekávané populační matice korelací (nebo kovariancí plus vektoru s průměry). Obvykle bychom k jejímu odhadu použili vzorek z předchozí studie, (nejlépe dotatečně velký a reprezentativní).

Dejme tomu, že máme pouze tři proměnné a tohle jsou očekávané korelace mezi nimi (na diagonále jsou jedničky, protože korelace proměnné se sebou samou je vždy 1):

cormat <- matrix(
  c(1.0, 0.3, 0.6,
    0.3, 1.0, 0.4,
    0.6, 0.4, 1.0),
  byrow = TRUE, nrow = 3,
  dimnames = list(
    c("x", "y", "z"),
    c("x", "y", "z"))
)
cormat

>#     x   y   z
># x 1.0 0.3 0.6
># y 0.3 1.0 0.4
># z 0.6 0.4 1.0

Pak můžeme na základě této korelační matice vygenerovat multivariačně normálních data. Pro zjednodušení pracujeme se standardizovanými proměnnými (čili sigma = 1 a mí = 0). Funkce mvrnorm() z balíčku MASS předpokládá multivariačně normální rozdělení, které je definováno variačně-kovariační maticí plus vektorem s průměry. Protože pracujeme s korelacemi, ale korelace jsou v podstatě kovariance mezi standardizovanými proměnnými, můžeme do funkce mvrnorm() zadat korelační matici a nastavit průměry proměnných na nulu.

Takto to vypadá, když si vygenerujeme např. 100 hodnot a sestavíme lineárně regresní model, kde x a z predikují y

df <- MASS::mvrnorm(100, 
                    mu = rep(0, times = 3), # Průměry proměnných
                    Sigma = cormat) %>%  # Variačně-kovariační matice
  as_tibble()

df

># # A tibble: 100 × 3
>#          x      y      z
>#      <dbl>  <dbl>  <dbl>
>#  1  0.326   1.79   0.811
>#  2 -0.0151  2.23  -0.201
>#  3 -1.92    0.439 -1.28 
>#  4  0.863  -0.932  0.427
>#  5  0.125  -0.208 -1.23 
>#  6 -1.07   -0.925 -0.337
>#  7 -0.238   1.68   1.26 
>#  8  0.105  -0.184 -0.574
>#  9 -0.152  -0.555  0.244
># 10  1.28   -0.181  1.45 
># # … with 90 more rows

# Korelace
df %>% 
  cor()

>#           x         y         z
># x 1.0000000 0.3656976 0.6197690
># y 0.3656976 1.0000000 0.4899691
># z 0.6197690 0.4899691 1.0000000

# Protože pracujeme se standardizovanými proměnnými,
# průsečík bude z definice nula
fit <- lm(y ~ 0 + x + z, data = df)
summary(fit)

># 
># Call:
># lm(formula = y ~ 0 + x + z, data = df)
># 
># Residuals:
>#     Min      1Q  Median      3Q     Max 
># -2.6743 -0.7955 -0.1335  0.5646  2.3310 
># 
># Coefficients:
>#   Estimate Std. Error t value Pr(>|t|)    
># x  0.09061    0.11596   0.781 0.436459    
># z  0.48784    0.12165   4.010 0.000118 ***
># ---
># Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
># 
># Residual standard error: 1.021 on 98 degrees of freedom
># Multiple R-squared:   0.25,  Adjusted R-squared:  0.2347 
># F-statistic: 16.33 on 2 and 98 DF,  p-value: 7.572e-07

Kdybychom chtěli simulovat power pro prediktory v lineárně regresním modelu, museli bychom si opět vytvořit vlastní funkci:

lm_fun <- function(simNum, # musí obsahovat argument simNum pro stanovení počtu simulací
                    n, # velikost vzorku
                    cov_mat, # kovariační nebo korelační matice 
                    mean_vector = NULL,# Vektor s průměry
                    formula) { # rovnice lineární regrese 

  # Pokud je mean vektor prázdný, předpokládat standardizované proměnné
  # s průměrem nula
  if (is.null(mean_vector)) {
    mu <- rep(0, ncol(cov_mat))
  } else {
    mu <- mean_vector
  }
  
  # Nejdříve si generujeme data z multivariačně normálního rozdělení
  data <- MASS::mvrnorm(n = n, mu = mu, Sigma = cov_mat)
  
  # Fitujeme lineárně regresní model
  fit <- lm(formula, data = as.data.frame(data))
  
  # Extrahujeme koeficienty
  # Funkce k extrakci koeficientů do jednoho atomického vektoru
  get_coefs <- function(fit) {
    out <- fit |> 
      summary() |> 
      coef() |> 
      as.data.frame() |> 
      setNames(c("b", "se", "t", "p"))
    
    row_names <- rownames(out)
    row_names[row_names == "(Intercept)"] <- "int"
    
    out <- out |> 
      dplyr::mutate(predictor = row_names) |> 
      tidyr::pivot_wider(names_from = predictor, values_from = -predictor)
    
    nms <- colnames(out)
    
    out <- as.numeric(out)
    names(out) <- nms
    
    out
  }
  
  # Výstupní hodnoty
  return(get_coefs(fit))
}

Drobný problém je, že pokud chceme použít více jader procesoru (paralelní výpočty) v následné simulace, musíme ve funkci výše používat pouze základní funkce z base R, anebo pokud se nejedná o funkce z base R, musíme se na ně explicitně odkazovat pomocí nazev_balicku::funkce. Proto je výše také použit nativní pipe operátor |> místo %>%.

Pokud ale funkci správně vytvoříme, můžeme pak pomocí grid_search() provést simulaci pro různé velikosti vzorku, např. 50, 100 a 200 osob.

pwr_lm <- grid_search(lm_fun, 
                      params = list(n = c(50, 100, 200)), # Velikosti vzorku
                      n.iter = 1000, # počet simulací
                      output = 'data.frame', # Formát výstupu
                      cov_mat = cormat, # kovariační matice
                      formula = y ~ 0 + x + z, # regresní rovnice
                      parallel = 'snow', # Paralelní výpočty 
                      ncpus = 6) # POčet použitých jader procesoru

># Running 3,000 tests...

Pak opět můžeme shrnout výsledky:

res <- results(pwr_lm) %>% 
  as_tibble()

res %>% 
  group_by(n.test) %>% 
  summarise(
    power_x = mean(p_x < 0.05), # Podíl p-hodnot menších než 0,05
    power_z = mean(p_z < 0.05)
  )

># # A tibble: 3 × 3
>#   n.test power_x power_z
>#    <dbl>   <dbl>   <dbl>
># 1     50   0.08    0.543
># 2    100   0.12    0.83 
># 3    200   0.215   0.994

Jak můžeme vidět, pro dosažení dostatečné statistické síly k identifikaci efektu prediktoru x by nám ani 200 osob nestačilo, zatímco pro identifikaci efektu prediktoru z by nám mělo stačit 100 až 200 osob.

6.3.1 Cvičení

Chceme pomocí lineární regrese zkoumat efekt veku (vekr), školního prospěchu (gpa), monitorování rodiči (monit), důvěry vůči vrstevníkům (duv_v) a optimismu na životní spokojenost (ziv_sp). Předpokládáme, že tyto proměnné mají přibližně multivariačně normální rozdělení a že variance a kovariance mezi tětito proměnnými jsou takovéto:

cov_mat <- structure(c(0.23, 0.02,-0.07, 0.10, 0.00, 0.81, 
                       0.02, 3.96, 0.56,-0.30, 0.18,-0.01, 
                      -0.07, 0.56, 0.56,-0.10, 0.01,-0.14, 
                       0.10,-0.30,-0.10, 0.32,-0.01, 0.30, 
                       0.00, 0.18, 0.01,-0.01, 0.26, 0.07, 
                       0.81,-0.01,-0.14, 0.30, 0.07, 9.71),
          dim = c(6L, 6L), 
          dimnames = list(c("ziv_sp", "vekr", "gpa", "monit", 
                            "duv_v", "optim"), 
                          c("ziv_sp", "vekr", "gpa", "monit", 
                            "duv_v", "optim")))
cov_mat

>#        ziv_sp  vekr   gpa monit duv_v optim
># ziv_sp   0.23  0.02 -0.07  0.10  0.00  0.81
># vekr     0.02  3.96  0.56 -0.30  0.18 -0.01
># gpa     -0.07  0.56  0.56 -0.10  0.01 -0.14
># monit    0.10 -0.30 -0.10  0.32 -0.01  0.30
># duv_v    0.00  0.18  0.01 -0.01  0.26  0.07
># optim    0.81 -0.01 -0.14  0.30  0.07  9.71

A průměry proměnných takovéto:

mean_vec <- c(ziv_sp = 2.9, vekr = 14.11, gpa = 2.2, monit = 2.8, 
              duv_v = 2.7, optim = 16.32)
mean_vec

># ziv_sp   vekr    gpa  monit  duv_v  optim 
>#   2.90  14.11   2.20   2.80   2.70  16.32

Využijte funkce lm_fun() a grid_search() a odhadněte statistickou sílu testů jednotlivých regresních koeficientů pro velikosti vzorku: 30, 100, 300, 500 a 1000 osob.