---
title: "Co nám statistický test dovoluje tvrdit?"
subtitle: "Od velikosti efektu k opatrnému vědeckému závěru"
lang: cs
filters:
  - ../theme/semantic-boxes.lua
format:
  html:
    toc: true
    toc-depth: 3
    toc-title: "Obsah"
    number-sections: true
    code-fold: false
    code-summary: "Zobrazit kód v R"
    embed-resources: true
    theme: [flatly, ../theme/skripta_theme.scss]
    highlight-style: github
    include-in-header:
      - "../theme/fonts-include.html"
      - "../theme/skripta-favicon.html"
  typst:
    toc: true
    toc-depth: 3
    number-sections: true
    output-file: "skripta_raw.pdf"
    include-in-header:
      - "skripta_theme.typ"
      - "skripta-logo.typ"
    date: last-modified
date: 2026/07/30
date-format: long
date-modified: last-modified
---

```{r}
#| label: nastaveni
#| include: false
if (!requireNamespace("here", quietly = TRUE)) install.packages("here")
if (!requireNamespace("readr", quietly = TRUE)) install.packages("readr")
if (!requireNamespace("dplyr", quietly = TRUE)) install.packages("dplyr")
if (!requireNamespace("tidyr", quietly = TRUE)) install.packages("tidyr")
if (!requireNamespace("purrr", quietly = TRUE)) install.packages("purrr")
if (!requireNamespace("ggplot2", quietly = TRUE)) install.packages("ggplot2")
if (!requireNamespace("maps", quietly = TRUE)) install.packages("maps")
if (!requireNamespace("glossary", quietly = TRUE)) install.packages("glossary")

here::i_am("Learning_materials/skripta.qmd")

# Visible code reads data/pie_crab_site_means.csv relative to the project root,
#   exactly as students do in their own R project.
knitr::opts_knit$set(root.dir = here::here())

# Hidden figures use Czech decimal commas on their axes. Visible chunks keep
#   R's default output so that it matches what students see in R.
knitr::knit_hooks$set(
  desetinna_carka = function(before) {
    if (before) {
      options(OutDec = ",")
    } else {
      options(OutDec = ".")
    }
  }
)
knitr::opts_hooks$set(
  echo = function(options) {
    if (isFALSE(options$echo)) {
      options$desetinna_carka <- TRUE
    }
    options
  }
)

source(here::here("R", "set_r_theme.R"))
source(here::here("R", "Functions", "format_cz.R"))
source(here::here("R", "Functions", "format_cz_math.R"))
source(here::here("R", "Functions", "simulovat_odhaleni_sklonu.R"))

local_pojmy <- here::here("..", "slovnik", "pojmy.yaml")

if (file.exists(local_pojmy)) {
  glossary::glossary_path(local_pojmy)
} else {
  remote_pojmy <- tempfile(fileext = ".yaml")
  tryCatch(
    {
      download.file(
        url = paste0(
          "https://raw.githubusercontent.com/",
          "CUNI-NATUR-Biostatistics/slovnik/main/pojmy.yaml"
        ),
        destfile = remote_pojmy,
        quiet = TRUE
      )
      glossary::glossary_path(remote_pojmy)
    },
    error = function(e) {
      message("Glossary YAML unavailable: ", conditionMessage(e))
    }
  )
}

local_glossary_helper <- here::here(
  "..",
  "slovnik",
  "R",
  "Functions",
  "render_glossary_term.R"
)

if (file.exists(local_glossary_helper)) {
  source(local_glossary_helper)
} else {
  tryCatch(
    {
      source(
        paste0(
          "https://raw.githubusercontent.com/",
          "CUNI-NATUR-Biostatistics/slovnik/main/",
          "R/Functions/render_glossary_term.R"
        )
      )
    },
    error = function(e) {
      source(here::here("R", "Functions", "render_glossary_term.R"))
    }
  )
}
```

```{r}
#| label: nastaveni-slovnik
#| echo: false
#| message: false
#| results: asis
if (knitr::is_html_output()) {
  glossary::glossary_style(
    color = biostat_cols[["indigo_velvet"]],
    text_decoration = "underline",
    def_bg = "#FFFFFF",
    def_color = "#2E2E2E"
  )
  cat(
    "<style>",
    "a.glossary:hover .def { display: inline-block; }",
    "a.glossary:hover .def::after {",
    "  content: ' '; position: absolute; top: 100%; left: 50%;",
    "  margin-left: -5px; border-width: 5px; border-style: solid;",
    "  border-color: #FFFFFF transparent transparent transparent;",
    "}",
    "a.glossary:active .def { display: none; }",
    "a.glossary:active .def::after { display: none; }",
    "</style>",
    "<script>",
    "document.addEventListener('DOMContentLoaded', function() {",
    "  document.querySelectorAll('a.glossary').forEach(function(a) {",
    "    a.addEventListener('click', function(e) { e.preventDefault(); });",
    "  });",
    "});",
    "</script>",
    sep = "\n"
  )
}
```

## Úvod

V předchozí lekci jsme k `r render_glossary_term("odhad", display = "odhadu")` `r render_glossary_term("sklon", display = "sklonu")` přidali jeho `r render_glossary_term("nejistota", display = "nejistotu")`: `r render_glossary_term("standardni-chyba", display = "standardní chybu")` a `r render_glossary_term("interval-spolehlivosti", display = "interval spolehlivosti")`. Ve výpisu `summary()` ale zůstaly dva sloupce, které jsme zatím přeskočili: `t value` a `Pr(>|t|)`. Patří k otázce, kterou si klade každá empirická věda: co můžeme na základě dat tvrdit a jaké pozorování by nás přimělo změnit názor?

### Měsíc ze sýra?

Představte si, že někdo pronese naprosto vážně:

::: {.box-data title="Tvrzení, které budeme zkoumat"}
**Měsíc je celý ze sýra.**
:::

Nejspíš se zasmějete. Přesto je na této větě něco užitečného: můžeme přemýšlet, jak bychom ji ověřovali. Mohli bychom zkoumat měsíční horniny, jejich složení nebo fyzikální vlastnosti.

Tvrzení je tedy **testovatelné**. To ale neznamená, že je **věrohodné**. Všechno, co už o Měsíci víme, mluví silně proti němu.

::: {.box-question title="Zkuste si — Čtyři velmi rozdílná tvrzení"}
Nejprve si odpovězte sami. Potom rozbalte možné odpovědi a porovnejte s nimi svou úvahu:

1. Dá se věta zkoumat pomocí pozorování?
2. Je podle současných znalostí věrohodná?
3. Jaké pozorování by vás přimělo změnit názor?

| Tvrzení |
|---|
| Měsíc je celý ze sýra |
| Všechny labutě jsou bílé |
| Holub stěhovavý vyhynul |
| Tento pták je nejkrásnější |
:::

::: {.box-result collapse="true" title="Možné odpovědi"}

| Tvrzení | Co lze pozorovat? | Co by silně změnilo závěr? |
|---|---|---|
| Měsíc je celý ze sýra | složení měsíčních hornin | důvěryhodný vzorek s vlastnostmi sýra |
| Všechny labutě jsou bílé | barvu nalezených labutí | jediná spolehlivě určená ne-bílá labuť |
| Holub stěhovavý vyhynul | přítomnost živých jedinců | jediný spolehlivě doložený živý jedinec |
| Tento pták je nejkrásnější | preference různých lidí | nejde o stejný druh empirického tvrzení |
:::

U vyhynulého holuba je navíc vidět důležitá nesymetrie. Nález živého jedince by závěr změnil okamžitě. Když žádného nenajdeme, jistota je menší: možná jsme hledali na špatném místě nebo jsme ho přehlédli.

::: {.box-question title="Rychlá zastávka — Testovatelné neznamená pravděpodobné"}
Student řekne: „Když lze tvrzení o sýrovém Měsíci testovat, musíme ho brát stejně vážně jako geologické vysvětlení.“ Co se v jeho úvaze pokazilo?
:::

::: {.box-result collapse="true" title="Možná odpověď"}
Testovatelnost říká jen to, že tvrzení můžeme porovnat s pozorováním. Kolik váhy mu dáme, závisí na tom, co už o světě víme. Geologické vysvětlení podporuje obrovské množství měření měsíčních hornin, sýrový Měsíc žádné. Obě tvrzení jsou testovatelná, ale nejsou stejně věrohodná.
:::

Věda tedy není stroj, který po jednom pokusu vytiskne **pravda** nebo **nepravda**. Záleží také na kvalitě měření, uspořádání studie, předpokladech modelu, předchozích znalostech a výsledcích dalších studií.

Statistický test bude jen jedním dílkem tohoto většího obrazu.

## Výsledky učení

Po prostudování této lekce budete schopni:

- rozlišit testovatelné vědecké tvrzení od hodnotového výroku a převést biologické tvrzení na `r render_glossary_term("nulova-hypoteza", display = "nulovou")` a alternativní hypotézu o `r render_glossary_term("sklon", display = "sklonu")` modelu;
- popsat, jak velký vztah jsme v datech našli, posoudit jeho biologický význam a říct, jak přesně jsme ho `r render_glossary_term("odhad", display = "odhadli")`;
- vysvětlit, jak z odhadu, jeho `r render_glossary_term("standardni-chyba", display = "standardní chyby")` a `r render_glossary_term("stupne-volnosti", display = "stupňů volnosti")` vzniká t-statistika a její referenční `r render_glossary_term("rozdeleni", display = "rozdělení")`, a správně interpretovat `r render_glossary_term("p-hodnota", display = "p-hodnotu")`;
- rozlišit chybu I. a II. typu v kontextu biologické otázky;
- napsat závěr v pořadí velikost `r render_glossary_term("efekt", display = "efektu")`, `r render_glossary_term("nejistota", display = "nejistota")`, slučitelnost s nulovou hypotézou a omezení studie.

## Každý řádek je jedno slanisko

Přesuňme se od Měsíce k otázce, u které odpověď předem neznáme:

::: {.box-question title="Hlavní biologická otázka"}
**Jsou samci atlantského kraba houslisty *Minuca pugnax* v `r render_glossary_term("prumer", display = "průměru")` větší na severněji položených slaniskách?**
:::

Biologové si už dlouho všímají, že organismy bývají v chladnějších oblastech větší. Tomuto vzoru se říká Bergmannovo pravidlo. Neplatí automaticky pro každý druh a samo nevysvětluje, proč by měl případný vztah vznikat. Platí něco podobného i pro naše kraby?

V létě 2016 vědci navštívili slaniska mezi Floridou a Massachusetts a na každém změřili dospělé samce. Pro hlavní analýzu máme jejich průměry: každý řádek tabulky představuje **jedno slanisko**. Naší `r render_glossary_term("jednotka-pozorovani", display = "jednotkou pozorování")` je proto slanisko, nikoli jednotlivý krab.

Soubor [pie_crab_site_means.csv](../data/pie_crab_site_means.csv) je součástí materiálů. Pro následující kód jej uložte do složky `data` ve svém projektu v R. Načteme ho a zkontrolujeme názvy a typy sloupců:

```{=typst}
#block(breakable: false)[
```

```{r}
#| label: nacist-a-prohlednout-lokality
#| echo: true
# Načteme tabulku s jedním řádkem pro každé slanisko
data_lokality <-
  read.csv(file = "data/pie_crab_site_means.csv")

# Zkontrolujeme názvy a typy sloupců
str(object = data_lokality)
```

```{=typst}
]
```

Tabulka má pět sloupců:

- `kod_lokality` a `nazev_lokality`: zkratka a název slaniska;
- `zemepisna_sirka`: zeměpisná šířka slaniska ve stupních severní šířky — náš `r render_glossary_term("prediktor", display = "prediktor")`;
- `pocet_krabu`: počet změřených dospělých samců;
- `prumerna_sirka_mm`: průměrná šířka krunýře samců v milimetrech — naše `r render_glossary_term("odezva", display = "odezva")`.

Kontrola dat pokračuje rozsahy číselných sloupců a chybějícími hodnotami:

```{=typst}
#block(breakable: false)[
```

```{r}
#| label: zkontrolovat-lokality
#| echo: true
# Vybereme číselné sloupce
vec_ciselne_sloupce <-
  c("zemepisna_sirka", "pocet_krabu", "prumerna_sirka_mm")

# Minimum, kvartily, průměr a maximum číselných sloupců
summary(object = data_lokality[, vec_ciselne_sloupce])

# Počet chybějících hodnot v každém sloupci
colSums(x = is.na(x = data_lokality))
```

```{=typst}
]
```

Tabulka obsahuje `r nrow(data_lokality)` slanisek mezi `r format_cz(min(data_lokality$zemepisna_sirka), digits = 1)` a `r format_cz(max(data_lokality$zemepisna_sirka), digits = 1)` ° s. š. Na každém bylo změřeno `r min(data_lokality$pocet_krabu)` až `r max(data_lokality$pocet_krabu)` samců a průměrná šířka krunýře se pohybuje od `r format_cz(min(data_lokality$prumerna_sirka_mm), digits = 1)` do `r format_cz(max(data_lokality$prumerna_sirka_mm), digits = 1)` mm. `r if (sum(is.na(data_lokality)) == 0) "Žádný sloupec nemá chybějící hodnoty." else paste0("Chybějících hodnot je ", sum(is.na(data_lokality)), ".")` Rozsahy odpovídají pobřeží mezi Floridou a Massachusetts i velikosti dospělých krabů, takže nic nenaznačuje chybu jednotek.

Celá tabulka je tak malá, že ji můžeme zobrazit najednou:

```{r}
#| label: zobrazit-prumery-lokalit
#| echo: false
#| output: true
#| message: false
#| warning: false
#| tbl-cap: "Průměry sledovaných slanisek, které použijeme v hlavním modelu."
data_lokality |>
  dplyr::transmute(
    `Lokalita` = kod_lokality,
    `Šířka (° s. š.)` =
      format_cz(zemepisna_sirka, digits = 1),
    `n` = pocet_krabu,
    `Průměr (mm)` =
      format_cz(prumerna_sirka_mm, digits = 2)
  ) |>
  knitr::kable()
```

Data pocházejí z výukového souboru [`lterdatasampler::pie_crab`](https://lter.github.io/lterdatasampler/articles/pie_crab_vignette.html), zveřejněného pod licencí CC0. Původní data jsou uložena v [Environmental Data Initiative](https://doi.org/10.6073/pasta/4c27d2e778d3325d3830a5142e3839bb) a související studie vyšla v časopise *Ecology and Evolution* ([Johnson et al., 2019](https://doi.org/10.1002/ece3.5883)).

```{r}
#| label: zobrazit-polohu-lokalit
#| echo: false
#| output: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 7.2
#| fig-cap: >
#|   Každý oranžový bod představuje jedno sledované slanisko podél
#|   atlantského pobřeží USA. Souřadnice pocházejí z tabulky 1
#|   související studie.
#| fig-alt: >
#|   Mapa atlantského pobřeží USA od Floridy po Massachusetts
#|   s oranžovými body sledovaných slanisek.
data_souradnice <-
  data.frame(
    kod_lokality = c(
      "GTM", "SI", "NIB", "ZI", "RC", "VCR", "DB",
      "JC", "CT", "NB", "CC", "BC", "PIE"
    ),
    zemepisna_delka = c(
      -81.327, -81.364, -79.205, -77.944, -76.676,
      -75.942, -75.437, -74.342, -71.981, -71.326,
      -70.068, -70.927, -70.849
    )
  )

data_mapa_lokality <-
  merge(
    x = data_lokality,
    y = data_souradnice,
    by = "kod_lokality"
  )

data_mapa_staty <-
  ggplot2::map_data("state") |>
  dplyr::filter(
    region %in% c(
      "florida", "georgia", "south carolina",
      "north carolina", "virginia", "maryland",
      "delaware", "new jersey", "new york",
      "connecticut", "rhode island", "massachusetts"
    )
  )

ggplot2::ggplot() +
  ggplot2::geom_polygon(
    data = data_mapa_staty,
    mapping = ggplot2::aes(
      x = long,
      y = lat,
      group = group
    ),
    fill = biostat_cols[["amethyst"]],
    alpha = 0.16,
    color = biostat_cols[["grey_olive"]],
    linewidth = 0.45
  ) +
  ggplot2::geom_point(
    data = data_mapa_lokality,
    mapping = ggplot2::aes(
      x = zemepisna_delka,
      y = zemepisna_sirka
    ),
    shape = 21,
    fill = biostat_cols[["orange"]],
    color = biostat_cols[["indigo_velvet"]],
    size = 4.4,
    stroke = 0.8
  ) +
  ggplot2::coord_quickmap(
    xlim = c(-82.5, -68.5),
    ylim = c(29, 44),
    expand = FALSE
  ) +
  ggplot2::labs(
    x = "Zeměpisná délka",
    y = "Zeměpisná šířka"
  ) +
  theme_biostat(base_size = 12)
```

```{r}
#| label: nacist-jednotlive-kraby
#| include: false
data_krabi <-
  readr::read_csv(
    file = here::here("data", "pie_crab.csv"),
    show_col_types = FALSE
  )
```

::: {.box-data collapse="true" title="Doplňující: odkud se vzaly průměry slanisek?"}
**Při prvním čtení můžete následující čtyři boxy „Doplňující“ přeskočit. Hlavní krabí příběh pokračuje grafem průměrů a přímkou.**

Původní tabulka obsahuje `r nrow(data_krabi)` jednotlivých krabů. Pro hlavní analýzu jsme je zprůměrovali po lokalitách.

Data jsou přímo součástí balíčku `lterdatasampler`. Následující kód je načte, srozumitelně přejmenuje, vytvoří průměry slanisek a uloží obě výukové tabulky. V dalších částech lekce už načítáme připravené CSV soubory.

```{r}
#| label: ukazat-stazeni-a-pripravu-dat
#| echo: true
#| eval: false
# Balíček stačí nainstalovat jednou:
# install.packages("lterdatasampler")

data("pie_crab", package = "lterdatasampler")

# Přejmenujeme sloupce tak, aby se s nimi v lekci dobře pracovalo.
data_krabi <-
  data.frame(
    datum = as.Date(pie_crab$date),
    zemepisna_sirka = pie_crab$latitude,
    kod_lokality = pie_crab$site,
    nazev_lokality = pie_crab$name,
    sirka_krunyre_mm = pie_crab$size
  )

# Pro každou lokalitu spočítáme průměrnou šířku krunýře.
data_prumery <-
  aggregate(
    sirka_krunyre_mm ~
      kod_lokality + nazev_lokality + zemepisna_sirka,
    data = data_krabi,
    FUN = mean
  )

# Zvlášť spočítáme počet změřených krabů na každé lokalitě.
data_pocty <-
  aggregate(
    sirka_krunyre_mm ~ kod_lokality,
    data = data_krabi,
    FUN = length
  )
```
:::

::: {.box-data collapse="true" title="Doplňující: spojení průměrů a počtů"}
Průměry a počty teď spojíme do jedné tabulky s jedním řádkem pro každé slanisko. Nakonec uložíme jak jednotlivá měření, tak připravené průměry.

```{r}
#| label: ukazat-spojeni-a-ulozeni-dat
#| echo: true
#| eval: false
# Obě souhrnné tabulky spojíme podle kódu lokality.
data_lokality <-
  merge(
    x = data_prumery,
    y = data_pocty,
    by = "kod_lokality",
    suffixes = c("_prumer", "_pocet")
  )

names(data_lokality)[
  names(data_lokality) == "sirka_krunyre_mm_prumer"
] <- "prumerna_sirka_mm"

names(data_lokality)[
  names(data_lokality) == "sirka_krunyre_mm_pocet"
] <- "pocet_krabu"

# Seřadíme lokality od jihu k severu a ponecháme potřebné sloupce.
data_lokality <-
  data_lokality[
    order(data_lokality$zemepisna_sirka),
    c(
      "kod_lokality",
      "nazev_lokality",
      "zemepisna_sirka",
      "pocet_krabu",
      "prumerna_sirka_mm"
    )
  ]
```
:::

::: {.box-data collapse="true" title="Doplňující: uložení výukových tabulek"}
Připravené tabulky uložíme jako CSV. Díky tomu může hlavní část lekce začít jednoduchým `read.csv()` a nemusí přípravu opakovat při každém renderu.

```{r}
#| label: ukazat-ulozeni-dat
#| echo: true
#| eval: false
# Uložíme jednotlivé kraby i připravené průměry lokalit.
write.csv(
  data_krabi,
  file = "data/pie_crab.csv",
  row.names = FALSE
)

write.csv(
  data_lokality,
  file = "data/pie_crab_site_means.csv",
  row.names = FALSE
)
```
:::

::: {.box-data collapse="true" title="Doplňující: jak vypadají všichni krabi?"}
Na obrázku jsou jednotliví krabi světle šedě. Velké body jsou průměry slanisek, které používáme v modelu.

```{r}
#| label: zobrazit-jednotlive-kraby-navic
#| echo: false
#| output: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 5.4
#| fig-cap: >
#|   Jednotliví krabi sdílejí zeměpisnou šířku svého slaniska.
#|   Velké body zobrazují průměry použité v hlavním modelu.
#| fig-alt: >
#|   Bodový graf šířky krunýře jednotlivých krabů proti zeměpisné
#|   šířce; světlé body jsou jedinci, velké tmavé body průměry slanisek.
ggplot2::ggplot() +
  ggplot2::geom_point(
    data = data_krabi,
    mapping = ggplot2::aes(
      x = zemepisna_sirka,
      y = sirka_krunyre_mm
    ),
    color = biostat_cols[["grey_olive"]],
    alpha = 0.18,
    position = ggplot2::position_jitter(
      width = 0.08,
      height = 0
    )
  ) +
  ggplot2::geom_point(
    data = data_lokality,
    mapping = ggplot2::aes(
      x = zemepisna_sirka,
      y = prumerna_sirka_mm
    ),
    color = biostat_cols[["graphite"]],
    size = 3.5
  ) +
  ggplot2::labs(
    x = "Zeměpisná šířka (° s. š.)",
    y = "Šířka krunýře (mm)"
  ) +
  theme_biostat(base_size = 12)
```

V této lekci proto pracujeme se `r nrow(data_lokality)` průměry lokalit. Proč nelze jednotlivé kraby bez dalšího považovat za `r nrow(data_krabi)` nezávislých odpovědí na otázku o zeměpisné šířce, probereme později u `r render_glossary_term("hierarchicka-data", display = "hierarchických dat")`.
:::

```{r}
#| label: spocitat-vysledky-modelu
#| include: false
mod_krabi <-
  lm(
    formula = prumerna_sirka_mm ~ zemepisna_sirka,
    data = data_lokality
  )

table_koeficienty <-
  summary(mod_krabi)$coefficients

intervaly_koeficientu <-
  confint(mod_krabi)

odhad_sklonu <-
  table_koeficienty[
    "zemepisna_sirka",
    "Estimate"
  ]

standardni_chyba_sklonu <-
  table_koeficienty[
    "zemepisna_sirka",
    "Std. Error"
  ]

interval_sklonu <-
  intervaly_koeficientu[
    "zemepisna_sirka",
  ]

pocet_lokalit <- nrow(data_lokality)
```

## Co říká známá přímka?

::: {.box-question title="Než se podíváte — Jaký obrázek by podpořil naši domněnku?"}
Než si hodnoty zobrazíte v grafu, načrtněte si očekávaný obrázek.

- Co bude na vodorovné ose?
- Co bude na svislé ose?
- Jak by měly body vypadat, kdyby byli krabi na severu větší?
- Jaký vzor by naopak mluvil proti této představě?
:::

Teď už můžeme body opravdu zobrazit. Každý bod je `r render_glossary_term("prumer", display = "průměr")` jednoho slaniska.

```{r}
#| label: zobrazit-prumery-a-primku
#| echo: false
#| output: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 5.3
#| fig-cap: >
#|   Průměrná šířka krunýře na sledovaných slaniskách.
#|   Fialová přímka je jednoduchý lineární model.
#| fig-alt: >
#|   Bodový graf průměrné šířky krunýře proti zeměpisné šířce
#|   s popisky lokalit a rostoucí fialovou přímkou modelu.
ggplot2::ggplot(
  data = data_lokality,
  mapping = ggplot2::aes(
    x = zemepisna_sirka,
    y = prumerna_sirka_mm
  )
) +
  ggplot2::geom_point(
    color = biostat_cols[["grey_olive"]],
    size = 3.5
  ) +
  ggplot2::geom_text(
    mapping = ggplot2::aes(label = kod_lokality),
    nudge_x = 0.18,
    nudge_y = 0.2,
    color = biostat_cols[["graphite"]],
    size = 3.1
  ) +
  ggplot2::geom_abline(
    intercept = coef(mod_krabi)[["(Intercept)"]],
    slope = odhad_sklonu,
    color = biostat_cols[["indigo_velvet"]],
    linewidth = 1.1
  ) +
  ggplot2::labs(
    x = "Zeměpisná šířka (° s. š.)",
    y = "Průměrná šířka krunýře (mm)"
  ) +
  theme_biostat(base_size = 12)
```

Body směrem doprava převážně stoupají. Na severnějších lokalitách tedy mají krabi v průměru širší krunýře. Přímku fitujeme stejným příkazem `lm()`, který už znáte:

```{r}
#| label: fitovat-model-prumeru
#| echo: true
#| output: true
# Fitujeme přímku přes průměry slanisek
mod_krabi <-
  lm(
    formula = prumerna_sirka_mm ~ zemepisna_sirka,
    data = data_lokality
  )

# Zobrazíme celý přehled modelu.
summary(object = mod_krabi)

# Zobrazíme 95% intervaly spolehlivosti obou koeficientů.
confint(object = mod_krabi)
```

V první tabulce už z předchozích lekcí umíme najít řádek se `r render_glossary_term("sklon", display = "sklonem")` (`zemepisna_sirka`). Ve sloupci `Estimate` je jeho `r render_glossary_term("odhad", display = "odhad")` `r format_cz(odhad_sklonu, digits = 3)` mm na stupeň. Vedle něj je známá `r render_glossary_term("standardni-chyba", display = "standardní chyba")` `r format_cz(standardni_chyba_sklonu, digits = 4)` mm na stupeň. Příkaz `confint()` pak ukazuje `r render_glossary_term("interval-spolehlivosti", display = "95% interval")` od `r format_cz(interval_sklonu[[1]], digits = 3)` do `r format_cz(interval_sklonu[[2]], digits = 3)` mm na stupeň.

Model odhaduje také `r render_glossary_term("intercept", display = "intercept")`. Ten by zde popisoval šířku krabů na rovníku, daleko mimo sledovaný rozsah, a pro naši biologickou otázku není zajímavý. Sloupce `t value` a `Pr(>|t|)`, zmíněné v úvodu, zatím necháme stranou.

### Dobrá praxe: po každém modelu zkontrolujeme residua

Výpis koeficientů není konec práce. Stejně jako v minulé lekci se hned po fitování modelu podíváme na `r render_glossary_term("residuum", display = "residua")` — rozdíly mezi pozorovanými průměry lokalit a hodnotami na přímce. Tuto kontrolu budeme dělat pokaždé, když fitujeme `r render_glossary_term("linearni-model", display = "lineární model")`.

```{r}
#| label: spocitat-residua-krabi
#| echo: true
#| output: false
# Z modelu získáme residua a hodnoty odhadnuté přímkou.
vec_residua_krabi <- resid(object = mod_krabi)
vec_odhadnute_hodnoty <- fitted(object = mod_krabi)
```

```{r}
#| label: pripravit-diagnostiku-krabi
#| include: false
data_diagnostika <-
  data_lokality |>
  dplyr::mutate(
    odhadnuta_hodnota = vec_odhadnute_hodnoty,
    residuum = vec_residua_krabi
  )
```

První graf ukazuje `r render_glossary_term("rozdeleni", display = "rozdělení")` residuí. Hodnota nula znamená, že bod leží přesně na přímce; kladná a záporná residua leží na opačných stranách přímky.

```{r}
#| label: zobrazit-hustotu-residui
#| echo: false
#| output: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 4.2
#| fig-cap: >
#|   Vyhlazený odhad rozdělení residuí krabího modelu.
#|   Krátké čárky dole označují jednotlivé lokality.
#| fig-alt: >
#|   Fialová křivka hustoty residuí s čárkovanou svislou čarou
#|   v nule a krátkými oranžovými čárkami pro jednotlivé lokality.
ggplot2::ggplot(
  data = data_diagnostika,
  mapping = ggplot2::aes(x = residuum)
) +
  ggplot2::geom_density(
    fill = biostat_cols[["amethyst"]],
    color = biostat_cols[["indigo_velvet"]],
    alpha = 0.55,
    linewidth = 1
  ) +
  ggplot2::geom_rug(
    color = biostat_cols[["orange"]],
    linewidth = 0.8
  ) +
  ggplot2::geom_vline(
    xintercept = 0,
    color = biostat_cols[["graphite"]],
    linetype = "dashed"
  ) +
  ggplot2::labs(
    x = "Residuum (mm)",
    y = "Hustota"
  ) +
  theme_biostat(base_size = 12)
```

Druhý graf vrací každé residuum k `r render_glossary_term("odhadnuta-hodnota", display = "hodnotě odhadnuté modelem")` pro danou lokalitu. Hledáme body rozptýlené kolem nuly bez zřetelného oblouku nebo trychtýře.

```{r}
#| label: zobrazit-residua-proti-odhadum
#| echo: false
#| output: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 5
#| fig-cap: >
#|   Residua proti hodnotám odhadnutým modelem. Popisky umožňují
#|   spojit neobvyklé residuum s konkrétní lokalitou.
#| fig-alt: >
#|   Bodový graf residuí proti odhadnutým hodnotám s čárkovanou
#|   vodorovnou čarou v nule a popisky lokalit.
ggplot2::ggplot(
  data = data_diagnostika,
  mapping = ggplot2::aes(
    x = odhadnuta_hodnota,
    y = residuum
  )
) +
  ggplot2::geom_hline(
    yintercept = 0,
    color = biostat_cols[["graphite"]],
    linetype = "dashed"
  ) +
  ggplot2::geom_point(
    color = biostat_cols[["grey_olive"]],
    size = 3
  ) +
  ggplot2::geom_text(
    mapping = ggplot2::aes(label = kod_lokality),
    nudge_y = 0.18,
    color = biostat_cols[["graphite"]],
    size = 3.2
  ) +
  ggplot2::labs(
    x = "Hodnota odhadnutá modelem (mm)",
    y = "Residuum (mm)"
  ) +
  theme_biostat(base_size = 12)
```

::: {.box-question title="Přečtěte diagnostické grafy"}

- Jsou kladná i záporná residua zastoupena kolem nuly?
- Vidíte v grafu residuí proti odhadnutým hodnotám oblouk, trychtýř nebo jiný systematický vzor?
- Které lokality mají největší residua?
- Co můžeme z tak malého počtu bodů posoudit jen velmi opatrně?
:::

Grafy neukazují nápadný systematický vzor, ale `r pocet_lokalit` lokalit je málo pro sebevědomou `r render_glossary_term("diagnostika-modelu", display = "diagnostiku")` tvaru rozdělení. Kontrola tedy neodhalila zřejmý problém; stále však musíme pamatovat na malý počet lokalit a jednoduchost modelu.

Stejný modelový výsledek budeme dále číst podmíněně: všechny následující závěry o velikosti vztahu a `r render_glossary_term("nejistota", display = "nejistotě")` předpokládají, že zvolená přímka, chování residuí a nezávislost lokalit jsou pro naši otázku přiměřené. Diagnostika nemůže dokázat správnost modelu, ale pomáhá odhalit zřetelný rozpor dříve, než začneme interpretovat další čísla z modelu.

Sklon `r format_cz(odhad_sklonu, digits = 3)` mm na jeden stupeň říká:

::: {.box-result title="Význam sklonu v našem modelu"}
Podle fitované přímky je na lokalitě o 1° severněji průměrná šířka krunýře větší přibližně o `r format_cz(odhad_sklonu, digits = 3)` mm.
:::

To je odpověď na otázku **„o kolik?“**. Takovému číselnému popisu změny říkáme `r render_glossary_term("efekt", display = "efekt")`; v tomto modelu ho nejlépe vyjadřuje sklon v milimetrech na stupeň. V dalších krocích se budeme ptát, jak přesný odhad je a jak neobvyklý by byl, kdyby ve skutečnosti žádný vztah neexistoval.

### Co znamená sklon v milimetrech na stupeň?

```{r}
#| label: prepocitat-efekt
#| include: false
rozdil_stupnu <- 5
rozsah_sirky <-
  diff(range(data_lokality$zemepisna_sirka))

efekt_zvoleny_rozdil <- odhad_sklonu * rozdil_stupnu
efekt_cely_rozsah <- odhad_sklonu * rozsah_sirky
```

Změna o `r format_cz(odhad_sklonu, digits = 2)` mm na jeden stupeň se těžko představuje. Co by takový sklon znamenal při cestě o `r rozdil_stupnu` stupňů na sever? A co mezi Floridou a Massachusetts?

- na 1°: `r format_cz(odhad_sklonu, digits = 3)` mm;
- na `r rozdil_stupnu`°: `r format_cz(efekt_zvoleny_rozdil, digits = 2)` mm;
- přes celý pozorovaný rozsah `r format_cz(rozsah_sirky, digits = 1)`°: `r format_cz(efekt_cely_rozsah, digits = 2)` mm.

```{r}
#| label: zobrazit-tri-meritka-efektu
#| echo: false
#| output: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 4.8
#| fig-cap: >
#|   Tentýž sklon přepočtený na 1°, na zvolený větší rozdíl
#|   a na celý rozsah sledovaných zeměpisných šířek.
#| fig-alt: >
#|   Tři fialové sloupce odhadnuté změny průměrné šířky krunýře
#|   pro rozdíl 1°, zvolený větší rozdíl a celý pozorovaný rozsah.
vec_popisky_meritek <-
  c(
    "1°",
    paste0(rozdil_stupnu, "°"),
    "Celý pozorovaný rozsah"
  )

data_meritka <-
  data.frame(
    meritko = factor(
      vec_popisky_meritek,
      levels = vec_popisky_meritek
    ),
    rozdil_sirky = c(
      1,
      rozdil_stupnu,
      rozsah_sirky
    )
  ) |>
  dplyr::mutate(
    zmena_mm = odhad_sklonu * rozdil_sirky
  )

ggplot2::ggplot(
  data = data_meritka,
  mapping = ggplot2::aes(
    x = meritko,
    y = zmena_mm
  )
) +
  ggplot2::geom_col(
    fill = biostat_cols[["amethyst"]],
    width = 0.68
  ) +
  ggplot2::geom_text(
    mapping = ggplot2::aes(
      label = paste0(
        format_cz(zmena_mm, digits = 2),
        " mm"
      )
    ),
    vjust = -0.4,
    color = biostat_cols[["graphite"]],
    size = 3.8
  ) +
  ggplot2::scale_y_continuous(
    limits = c(
      0,
      max(data_meritka$zmena_mm) * 1.18
    )
  ) +
  ggplot2::labs(
    x = "Rozdíl zeměpisné šířky",
    y = "Odhadnutá změna průměrné šířky krunýře"
  ) +
  theme_biostat(base_size = 12)
```

Průměrné šířky krunýře se na slaniskách pohybují přibližně mezi `r format_cz(min(data_lokality$prumerna_sirka_mm), digits = 1)` a `r format_cz(max(data_lokality$prumerna_sirka_mm), digits = 1)` mm. Rozdíl kolem `r format_cz(efekt_cely_rozsah, digits = 2)` mm přes celé sledované pobřeží proto stojí za biologickou pozornost.

::: {.box-question title="Zkuste si — Je to biologicky hodně?"}

1. Které měřítko se vám interpretuje nejsnáze: 1°, `r rozdil_stupnu`°, nebo celý rozsah?
2. Je změna přibližně `r format_cz(efekt_cely_rozsah, digits = 2)` mm velká vzhledem k šířce krabů v tabulce?
3. Jakou další biologickou informaci byste potřebovali, abyste význam rozdílu posoudili lépe?
4. Mohla by studie najít velmi přesně změnu, která je biologicky zanedbatelná?
:::

::: {.box-result collapse="true" title="Možné odpovědi"}
Přes celé sledované pobřeží odpovídá sklon změně asi `r format_cz(efekt_cely_rozsah, digits = 2)` mm, tedy zhruba `r format_cz(100 * efekt_cely_rozsah / mean(data_lokality$prumerna_sirka_mm), digits = 0)` % průměrné šířky krunýře. Takový rozdíl lze v terénu změřit a může souviset například s plodností nebo soubojem samců. K lepšímu posouzení bychom potřebovali vědět, jak velké rozdíly ve velikosti mají pro kraby biologické důsledky. A ano: velmi přesná studie, třeba s tisíci lokalitami, může přesně odhadnout i změnu, která je biologicky zanedbatelná.
:::

## Jeden odhad nestačí

`r render_glossary_term("sklon", display = "Sklon")` `r format_cz(odhad_sklonu, digits = 3)` mm na stupeň je náš nejlepší bodový `r render_glossary_term("odhad", display = "odhad")`. Naše slaniska ale nejsou jediná, která jsme mohli navštívit. Jiný `r render_glossary_term("vyber", display = "výběr")` by nejspíš dal trochu jiný sklon.

Právě proto jsme se v minulé lekci učili číst `r render_glossary_term("standardni-chyba", display = "standardní chybu")` a `r render_glossary_term("interval-spolehlivosti", display = "interval spolehlivosti")`.

U krabího sklonu máme:

- odhad `r format_cz(odhad_sklonu, digits = 3)` mm na stupeň;
- standardní chybu `r format_cz(standardni_chyba_sklonu, digits = 4)` mm na stupeň;
- 95% interval od `r format_cz(interval_sklonu[[1]], digits = 3)` do `r format_cz(interval_sklonu[[2]], digits = 3)` mm na stupeň.

```{r}
#| label: zobrazit-odhad-a-interval
#| echo: false
#| output: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 4.2
#| fig-cap: >
#|   Odhad sklonu a jeho 95% interval spolehlivosti.
#| fig-alt: >
#|   Oranžový bod odhadu sklonu s fialovým vodorovným intervalem
#|   spolehlivosti, který celý leží v kladných hodnotách.
data_interval <-
  data.frame(
    vysledek = "Sklon velikosti k zeměpisné šířce",
    odhad = odhad_sklonu,
    dolni = interval_sklonu[[1]],
    horni = interval_sklonu[[2]]
  )

ggplot2::ggplot(
  data = data_interval,
  mapping = ggplot2::aes(
    x = odhad,
    y = vysledek
  )
) +
  ggplot2::geom_errorbar(
    mapping = ggplot2::aes(
      xmin = dolni,
      xmax = horni
    ),
    orientation = "y",
    width = 0.2,
    color = biostat_cols[["indigo_velvet"]],
    linewidth = 1.4
  ) +
  ggplot2::geom_point(
    color = biostat_cols[["orange"]],
    size = 4
  ) +
  ggplot2::labs(
    x = "Změna průměrné šířky krunýře (mm na 1° severněji)",
    y = NULL
  ) +
  theme_biostat(base_size = 12)
```

::: {.box-question title="Zastavte se u obrázku — Co všechno ještě připouští interval?"}

- Kde leží bodový odhad?
- Jaké malé i velké změny připouští celý interval?
- Jsou všechny biologicky stejně zajímavé?
- Co by pro naše sebevědomí znamenal mnohem širší interval?
:::

Všechny hodnoty v intervalu popisují růst velikosti směrem na sever. Liší se v tom, jak rychlý by tento růst mohl být.

## Co kdyby ve skutečnosti žádný trend nebyl?

Z modelu už známe velikost vztahu i `r render_glossary_term("nejistota", display = "nejistotu")` `r render_glossary_term("odhad", display = "odhadu")`. K našemu výsledku teď postavíme jednu přesnou soupeřící možnost.

Představte si skeptickou kolegyni. Podívá se na graf a řekne:

> Možná mezi zeměpisnou šířkou a `r render_glossary_term("prumer", display = "průměrnou velikostí")` krabů v populaci žádný lineární vztah není. Těchto `r pocet_lokalit` bodů náhodou vyšlo rostoucích.

Jak by takový svět vypadal v našem modelu? Kdyby měla přímka `r render_glossary_term("sklon", display = "sklon")` **nula**, průměrná velikost krabů by se při cestě na sever neměnila. To ale neznamená, že by všechny lokality měly úplně stejný průměr.

I ve světě bez severojižního trendu by mezi lokalitami zůstávaly rozdíly. Mohly by vznikat například:

- náhodným `r render_glossary_term("vyber", display = "výběrem")` jiných krabů;
- nepřesností měření;
- místními podmínkami, které v modelu nemáme, třeba potravou nebo teplotou;
- dalšími biologickými rozdíly mezi slanisky.

Budeme tedy skládat simulovaná data ze dvou částí: **systematický sklon nastavíme na nulu**, ale kolem vodorovné přímky ponecháme **náhodnou `r render_glossary_term("variabilita", display = "variabilitu")`**.

Tuto druhou část už známe z `r render_glossary_term("diagnostika-modelu", display = "diagnostiky modelu")`. Rozdíl mezi hodnotou lokality a hodnotou na přímce jsme nazvali `r render_glossary_term("residuum", display = "residuum")`. V simulaci ke každé lokalitě vytvoříme nové náhodné residuum a přičteme ho k vodorovné přímce. Residua budou v průměru nulová, ale jednotlivě mohou být kladná i záporná. Jejich typickou velikost nastavíme podle `r render_glossary_term("rezidualni-smerodatna-odchylka", display = "reziduální směrodatné odchylky")` našeho krabího modelu, přibližně `r format_cz(sigma(mod_krabi), digits = 2)` mm.

Právě tato simulovaná residua způsobí, že přímka fitovaná jednomu náhodnému souboru nemusí vyjít přesně vodorovná — stejně jako skutečná residua způsobují, že pozorované body neleží přesně na fitované přímce.

Porovnáváme tedy dvě možnosti:

| Možnost | Co říká o krabech |
|---|---|
| Nulová možnost | Průměrná velikost se se zeměpisnou šířkou lineárně nemění. Sklon je 0 mm na stupeň. |
| Alternativní možnost | Průměrná velikost se se zeměpisnou šířkou lineárně mění. Sklon není 0 mm na stupeň. |

Nula se zde nevzala z obecného statistického pravidla. Patří k naší konkrétní otázce: v `r render_glossary_term("linearni-model", display = "lineárním modelu")` znamená sklon 0 přesně **žádnou změnu průměrné velikosti při cestě na sever**.

Tutéž dvojici můžeme zapsat zkráceně. Symbol $\beta_{\text{zem. šířka}}$ označuje sklon vztahu k zeměpisné šířce v širší populaci lokalit, který se snažíme odhadnout pomocí našich `r pocet_lokalit` slanisek:

$$
H_0:\ \beta_{\text{zem. šířka}}=0
\qquad\text{a}\qquad
H_A:\ \beta_{\text{zem. šířka}}\ne0
$$

`r render_glossary_term("nulova-hypoteza", display = "Nulová hypotéza")` $H_0$ je právě nulová možnost z tabulky. Alternativní hypotéza $H_A$ připouští růst i pokles. <!-- TODO: add slug alternativni-hypoteza to pojmy.yaml -->

Nulová hypotéza tedy **neříká**, že každý soubor `r pocet_lokalit` slanisek musí dát odhad sklonu přesně nula. Říká, že dlouhodobý systematický trend je nulový. Náhodná variabilita kolem něj může v jednotlivém souboru vytvořit mírně rostoucí i mírně klesající přímku.

Tento nulový svět teď prozkoumáme simulací. Ponecháme stejné zeměpisné šířky jako v našich datech. Pro každou z nich vygenerujeme průměrnou velikost jako **stejnou základní úroveň plus nové náhodné residuum**. V simulaci tedy není žádný systematický růst k severu, ale lokality se stále liší. Každému novému souboru potom fitujeme přímku.

::: {.box-question title="Předpovězte výsledek simulace"}

- Budou všechny odhadnuté sklony přesně nulové?
- Kolem jaké hodnoty se budou shlukovat?
- Mohou občas vzniknout i výrazně kladné nebo záporné sklony?
- Kam by asi patřil náš sklon `r format_cz(odhad_sklonu, digits = 3)` mm na stupeň?
:::

Zkusme si taková data vyrobit. Všechny lokality dostanou stejnou očekávanou velikost, takže vložený sklon je nula. Ke každé lokalitě ale přidáme nové residuum podobně velké jako residua kolem skutečné krabí přímky. Začneme dvaceti takovými náhodnými studiemi.

```{r}
#| label: simulovat-sklony-pri-nule
#| include: false
set.seed(900723)

pocet_simulaci <- 5000
prumer_nuloveho_modelu <-
  mean(data_lokality$prumerna_sirka_mm)
rezidualni_sd <- sigma(mod_krabi)

data_nulove_ukazky <-
  purrr::map(
    .x = seq_len(20),
    .f = ~ {
      data.frame(
        simulace = .x,
        zemepisna_sirka = data_lokality$zemepisna_sirka,
        simulovana_sirka =
          prumer_nuloveho_modelu +
          stats::rnorm(
            n = pocet_lokalit,
            mean = 0,
            sd = rezidualni_sd
        )
      )
    }
  ) |>
  purrr::list_rbind()

vec_nulove_sklony <-
  purrr::map_dbl(
    .x = seq_len(pocet_simulaci),
    .f = ~ {
      data_nulovy_vyber <-
        data_lokality |>
        dplyr::mutate(
          simulovana_sirka =
            prumer_nuloveho_modelu +
            stats::rnorm(
              n = pocet_lokalit,
              mean = 0,
              sd = rezidualni_sd
            )
        )

      mod_nulovy <-
        lm(
          formula = simulovana_sirka ~ zemepisna_sirka,
          data = data_nulovy_vyber
        )

      coef(mod_nulovy)[["zemepisna_sirka"]]
    }
  )

data_nulove_sklony <-
  data.frame(
    sklon = vec_nulove_sklony
  )
```

```{r}
#| label: zobrazit-dvacet-vyberu-pri-nule
#| echo: false
#| output: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 8
#| fig-cap: >
#|   Dvacet náhodných studií vytvořených ve světě s nulovým
#|   sklonem. Od vodorovné přímky jsme ke každé lokalitě
#|   přičetli nové náhodné residuum. Každá studie obsahuje
#|   stejný počet bodů jako naše data a vlastní fitovanou přímku.
#| fig-alt: >
#|   Dvacet malých bodových grafů simulovaných studií bez trendu;
#|   fialové fitované přímky v nich mírně stoupají i klesají.
ggplot2::ggplot(
  data = data_nulove_ukazky,
  mapping = ggplot2::aes(
    x = zemepisna_sirka,
    y = simulovana_sirka
  )
) +
  ggplot2::geom_point(
    color = biostat_cols[["grey_olive"]],
    size = 1.1,
    alpha = 0.75
  ) +
  ggplot2::geom_smooth(
    method = "lm",
    formula = y ~ x,
    se = FALSE,
    color = biostat_cols[["indigo_velvet"]],
    linewidth = 0.7
  ) +
  ggplot2::facet_wrap(
    facets = ggplot2::vars(simulace),
    ncol = 5
  ) +
  ggplot2::labs(
    x = "Zeměpisná šířka",
    y = "Simulovaný průměr (mm)"
  ) +
  theme_biostat(base_size = 9) +
  ggplot2::theme(
    strip.text = ggplot2::element_text(size = 7),
    axis.text = ggplot2::element_text(size = 6)
  )
```

Přestože jsme žádný trend nevložili, některé přímky stoupají a jiné klesají. Náhoda tedy dokáže vytvořit nenulový odhad. Stejnou studii zopakujeme `r format(pocet_simulaci, big.mark = " ")`krát a z každé si ponecháme jen odhadnutý sklon.

```{r}
#| label: zobrazit-sklony-pri-nule
#| echo: false
#| output: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 5
#| fig-cap: !expr paste0("Sklony z ", format(pocet_simulaci, big.mark = " "), " simulovaných souborů, ve kterých byl skutečný sklon nulový. Šedá čárkovaná čára označuje nulu, oranžová pozorovaný sklon.")
#| fig-alt: >
#|   Šedý histogram simulovaných sklonů soustředěný kolem nuly;
#|   oranžová čára pozorovaného sklonu leží daleko vpravo.
ggplot2::ggplot(
  data = data_nulove_sklony,
  mapping = ggplot2::aes(x = sklon)
) +
  ggplot2::geom_histogram(
    bins = 55,
    fill = biostat_cols[["grey_olive"]],
    color = "white",
    linewidth = 0.2
  ) +
  ggplot2::geom_vline(
    xintercept = 0,
    color = biostat_cols[["graphite"]],
    linetype = "dashed",
    linewidth = 0.9
  ) +
  ggplot2::geom_vline(
    xintercept = odhad_sklonu,
    color = biostat_cols[["orange"]],
    linewidth = 1.1
  ) +
  ggplot2::annotate(
    geom = "text",
    x = odhad_sklonu,
    y = Inf,
    label = "náš sklon",
    vjust = 1.4,
    hjust = 1.05,
    color = biostat_cols[["orange"]],
    size = 3.6
  ) +
  ggplot2::labs(
    x = "Odhadnutý sklon (mm na 1° severněji)",
    y = "Počet simulací"
  ) +
  ggplot2::coord_cartesian(
    xlim = c(
      min(data_nulove_sklony$sklon),
      odhad_sklonu + 0.07
    )
  ) +
  theme_biostat(base_size = 12)
```

Většina sklonů leží blízko nuly. Náš pozorovaný sklon je daleko vpravo. Abychom mohli říct **jak daleko**, potřebujeme vzdálenost vyjádřit vzhledem k nejistotě odhadu.

## Jak daleko jsme od vodorovné přímky?

Nejprve si sepíšeme známá čísla:

```{r}
#| label: tabulka-pro-vypocet-t
#| echo: false
#| output: true
t_pozorovane <-
  table_koeficienty[
    "zemepisna_sirka",
    "t value"
  ]

data.frame(
  `Co známe` = c(
    "Pozorovaný sklon (Estimate)",
    "Sklon vodorovné přímky (nulová hodnota)",
    "Rozdíl od nuly (Estimate − nulová hodnota)",
    "Standardní chyba sklonu (Std. Error)"
  ),
  `Hodnota` = c(
    paste0(
      format_cz(odhad_sklonu, digits = 3),
      " mm na 1°"
    ),
    "0 mm na 1°",
    paste0(
      format_cz(odhad_sklonu, digits = 3),
      " mm na 1°"
    ),
    paste0(
      format_cz(standardni_chyba_sklonu, digits = 4),
      " mm na 1°"
    )
  ),
  check.names = FALSE
) |>
  knitr::kable()
```

Otázka zní: **Kolikrát se `r render_glossary_term("standardni-chyba", display = "standardní chyba")` vejde do vzdálenosti mezi naším `r render_glossary_term("sklon", display = "sklonem")` a nulou?**

$$
\frac{
\text{pozorovaný sklon}-\text{nulový sklon}
}{
\text{standardní chyba sklonu}
}
$$

Dosadíme:

$$
\frac{
`r format_cz_math(odhad_sklonu, digits = 3)`-0
}{
`r format_cz_math(standardni_chyba_sklonu, digits = 4)`
}
\approx
`r format_cz_math(t_pozorovane, digits = 2)`.
$$

Náš sklon tedy leží přibližně `r format_cz(t_pozorovane, digits = 2)` standardních chyb od nuly.

Statistici tuto vzdálenost označují písmenem $t$. Stejný výpočet se obecně zapisuje:

$$
t=
\frac{
\hat{\beta}_{\text{šířka}}-0
}{
\operatorname{SE}(\hat{\beta}_{\text{šířka}})
}.
$$

Symbol $\hat{\beta}_{\text{šířka}}$ označuje `r render_glossary_term("odhad", display = "odhadnutý")` sklon; stříška odlišuje odhad z našich dat od neznámého populačního sklonu $\beta_{\text{zem. šířka}}$. Nula je hodnota sklonu podle `r render_glossary_term("nulova-hypoteza", display = "nulové hypotézy")` a $\operatorname{SE}(\hat{\beta}_{\text{šířka}})$ standardní chyba odhadu. Výslednému číslu se říká t-statistika. <!-- TODO: add slug t-statistika to pojmy.yaml -->

::: {.box-question title="Zkuste si — Spočítejte t a řekněte, co znamená"}

1. Odečtěte nulový sklon od pozorovaného.
2. Vydělte rozdíl standardní chybou.
3. Doplňte větu: „Pozorovaný sklon leží přibližně …“
4. Vysvětlete, proč t-statistika není velikostí biologického `r render_glossary_term("efekt", display = "efektu")`.
:::

T-statistika je **standardizovaná vzdálenost od jedné konkrétní hodnoty**. Velké $|t|$ může vzniknout velkým efektem, vysokou přesností, nebo obojím.

## Malá odbočka: kolik volnosti zbývá?

Stejná hodnota $t$ není stejně neobvyklá ve studii se čtyřmi body a ve studii se čtyřmi tisíci body. Při malém množství informace `r render_glossary_term("odhad", display = "odhadujeme")` šum méně stabilně a musíme být opatrnější.

Tuto „zbývající informaci“ popisují `r render_glossary_term("stupne-volnosti", display = "stupně volnosti")`.

### Sudoku: kdy poslední volba přestane být volbou

Statistické stupně volnosti si zkusíme ilustrovat pomocí Sudoku. Sudoku je hlavolam s mřížkou 9 × 9 polí. Mřížka je rozdělena na devět menších čtverců 3 × 3. Číslice 1 až 9 se musí objevit právě jednou v každém řádku, sloupci i malém čtverci.

Podívejme se jen na jeden čtverec 3 × 3. Vlevo známe tři číslice a pro prázdná pole ještě zbývá mnoho možností. Vpravo už jsou umístěny číslice 1 až 8. Poslední pole nemá na výběr: musí v něm být 9.

```{r}
#| label: zobrazit-sudoku-volnost
#| echo: false
#| output: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 4.2
#| fig-cap: >
#|   V neúplném čtverci zbývá mnoho možností. Po umístění osmi
#|   různých číslic je poslední hodnota vynucená.
#| fig-alt: >
#|   Dva čtverce Sudoku 3 × 3: vlevo jsou jen tři číslice,
#|   vpravo osm číslic a oranžové pole s otazníkem pro vynucenou devítku.
data_sudoku <-
  tidyr::expand_grid(
    varianta = c(
      "Málo čísel: mnoho možností",
      "Osm čísel: poslední je vynucené"
    ),
    sloupec = seq_len(3),
    radek = seq_len(3)
  ) |>
  dplyr::arrange(
    varianta,
    radek,
    sloupec
  ) |>
  dplyr::group_by(varianta) |>
  dplyr::mutate(
    pozice = dplyr::row_number(),
    hodnota = dplyr::case_when(
      varianta == "Málo čísel: mnoho možností" &
        pozice == 1 ~ "1",
      varianta == "Málo čísel: mnoho možností" &
        pozice == 5 ~ "4",
      varianta == "Málo čísel: mnoho možností" &
        pozice == 9 ~ "7",
      varianta == "Osm čísel: poslední je vynucené" &
        pozice < 9 ~ as.character(pozice),
      varianta == "Osm čísel: poslední je vynucené" &
        pozice == 9 ~ "?",
      TRUE ~ ""
    ),
    vynucene = hodnota == "?"
  ) |>
  dplyr::ungroup()

ggplot2::ggplot(
  data = data_sudoku,
  mapping = ggplot2::aes(
    x = sloupec,
    y = -radek
  )
) +
  ggplot2::geom_tile(
    mapping = ggplot2::aes(fill = vynucene),
    color = biostat_cols[["graphite"]],
    linewidth = 0.8
  ) +
  ggplot2::geom_text(
    mapping = ggplot2::aes(label = hodnota),
    color = biostat_cols[["indigo_velvet"]],
    size = 8,
    fontface = "bold"
  ) +
  ggplot2::facet_wrap(
    facets = ggplot2::vars(varianta),
    nrow = 1
  ) +
  ggplot2::scale_fill_manual(
    values = c(
      "FALSE" = "white",
      "TRUE" = biostat_cols[["orange"]]
    ),
    guide = "none"
  ) +
  ggplot2::coord_equal() +
  ggplot2::scale_x_continuous(
    breaks = NULL,
    expand = c(0, 0)
  ) +
  ggplot2::scale_y_continuous(
    breaks = NULL,
    expand = c(0, 0)
  ) +
  ggplot2::labs(
    x = NULL,
    y = NULL
  ) +
  theme_biostat(base_size = 12) +
  ggplot2::theme(
    panel.grid = ggplot2::element_blank(),
    axis.text = ggplot2::element_blank(),
    axis.ticks = ggplot2::element_blank()
  )
```

Statistické stupně volnosti jsou trochu podobné této situaci v Sudoku: říkají nám, **kolik volnosti k volbě ještě zbývá**, když už platí určitá omezení. Sudoku tuto myšlenku dobře zviditelní, i když jeho pravidla nejsou statistickým výpočtem.

### Čtyři čísla s průměrem 10

```{r}
#| label: pripravit-priklad-pevneho-prumeru
#| include: false
pocet_cisel_df <- 4
pevny_prumer_df <- 10
cilovy_soucet_df <- pocet_cisel_df * pevny_prumer_df

prvni_cislo_df <- 7
druhe_cislo_df <- 12
treti_cislo_df <- 9
vynucene_cislo_df <-
  cilovy_soucet_df -
  prvni_cislo_df -
  druhe_cislo_df -
  treti_cislo_df

df_pevny_prumer <- pocet_cisel_df - 1
```

Požadovaný součet získáme jako **počet čísel krát jejich pevný průměr**. Chceme `r pocet_cisel_df` čísla s `r render_glossary_term("prumer", display = "průměrem")` `r pevny_prumer_df`, takže jejich součet musí být `r pocet_cisel_df` krát `r pevny_prumer_df`, tedy `r cilovy_soucet_df`.

$$
\text{požadovaný součet}
=
\text{počet čísel}
\times
\text{pevný průměr}.
$$

Zvolme tři čísla. Poslední už musí doplnit součet do `r cilovy_soucet_df`:

$$
`r prvni_cislo_df`+`r druhe_cislo_df`+`r treti_cislo_df`+x=`r cilovy_soucet_df`
$$

$$
x=`r cilovy_soucet_df`-`r prvni_cislo_df`-`r druhe_cislo_df`-`r treti_cislo_df`=`r vynucene_cislo_df`.
$$

Co když známe jen první dvě čísla?

$$
`r prvni_cislo_df`+`r druhe_cislo_df`+x+y=`r cilovy_soucet_df`
$$

Teď neznáme ani $x$, ani $y$. Z jediné rovnice je neumíme oba určit: $x$ si ještě můžeme zvolit a tím je $y$ vynuceno tak, aby součet zůstal `r cilovy_soucet_df`. Ve srovnání s prvním příkladem tedy zbývá **o jednu volbu více**.

Nejde nám o vyřešení rovnice. Sledujeme, kolik hodnot ještě můžeme zvolit nezávisle a kolik jich už musí splnit omezení pevného součtu.

::: {.box-question title="Zkuste vlastní čtveřici"}
Napište tři libovolná čísla a sestavte rovnici pro čtvrté tak, aby jejich průměr byl `r pevny_prumer_df`. Kolik hodnot jste mohli zvolit nezávisle? Která hodnota už byla vynucená?
:::

Pro čtyři čísla s pevným průměrem platí:

$$
\text{stupně volnosti}
=
\text{počet čísel}
-
\text{jedno omezení}
$$

$$
df=`r pocet_cisel_df`-1=`r df_pevny_prumer`.
$$

První tři čísla mohou být zvolena volně. Jakmile známe je a požadovaný průměr, čtvrté musí doplnit součet do `r cilovy_soucet_df`.

Obecně pro $n$ čísel s jedním pevným průměrem platí:

$$
df=n-1,
$$

kde $n$ je počet čísel a jednička představuje jedno omezení jejich součtu.

### Zpět ke krabům

```{r}
#| label: spocitat-stupne-volnosti
#| include: false
stupne_volnosti <- df.residual(mod_krabi)
```

Krabí model má `r pocet_lokalit` průměrů slanisek. To je náš počet `r render_glossary_term("pozorovani", display = "pozorování")`. Z těchto pozorování model odhaduje dvě hodnoty:

- `r render_glossary_term("intercept", display = "intercept")`;
- `r render_glossary_term("sklon", display = "sklon")`.

Po odhadu těchto dvou koeficientů zbývá:

$$
\text{reziduální stupně volnosti}
=
\text{počet pozorování (zde počet lokalit)}
-
\text{počet odhadnutých koeficientů}
$$

$$
df=`r pocet_lokalit`-2=`r stupne_volnosti`.
$$

Stejný vztah můžeme pro `r render_glossary_term("linearni-model", display = "lineární model")` zapsat obecně:

$$
df_{\mathrm{res}}=n-p,
$$

kde $n$ je počet pozorování — u krabů počet lokalit — a $p$ je počet koeficientů odhadnutých z těchto pozorování. Zde tedy $n=`r pocet_lokalit`$ a $p=2$.

::: {.box-question title="Spojte tři obrazy — Kde přesně analogie končí?"}
Doplňte tři věty:

1. Sudoku ukazuje, že …
2. Čtyři čísla s pevným průměrem ukazují přesněji, že …
3. V regresi odhad interceptu a sklonu znamená, že …

:::

::: {.box-result collapse="true" title="Možné odpovědi"}
1. Sudoku ukazuje, že po splnění pravidel může poslední volba přestat být volbou.
2. Čtyři čísla s pevným průměrem ukazují přesněji, že jedno omezení součtu spotřebuje jeden stupeň volnosti.
3. V regresi odhad interceptu a sklonu znamená, že dva stupně volnosti spotřebujeme na tyto dva koeficienty a pro `r render_glossary_term("residuum", display = "residua")` jich zbývá $n-2$.

Analogie končí tam, kde Sudoku vynucuje přesnou číslici. V regresi residua vynucená nejsou; stupně volnosti jen říkají, kolik nezávislé informace o jejich rozptýlení zbývá.
:::

Teď si spojme počet pozorování s tím, kolik volnosti po fitování přímky zbývá. Přímka pokaždé odhaduje dva koeficienty, intercept a sklon:

| Počet nezávislých pozorování | Odhadnuté koeficienty | Zbývající reziduální df |
|---:|---:|---:|
| 3 | 2 | 1 |
| `r pocet_lokalit` | 2 | `r stupne_volnosti` |
| 10 002 | 2 | 10 000 |

Sudoku ukázalo volnost, která zbývá při postupném doplňování hodnot. Příklad s pevným průměrem ukázal stejnou myšlenku přesnou rovnicí. V regresi spotřebujeme část informace na odhad interceptu a sklonu. Když přidáme další nezávislá pozorování, ale stále odhadujeme jen tyto dva koeficienty, zbývá více reziduálních stupňů volnosti.

### Proč jsme s malým počtem lokalit opatrnější?

Pro vyhodnocení t-statistiky používáme `r render_glossary_term("rozdeleni", display = "t-rozdělení")`. Křivka ukazuje, které hodnoty $t$ bychom při opakování studie čekali častěji a které méně často. **Celá plocha pod křivkou představuje všechny možné výsledky dohromady.**

Vyberme si dobře viditelnou hranici: hodnoty alespoň dvě `r render_glossary_term("standardni-chyba", display = "standardní chyby")` od nuly, tedy $t\leq-2$ nebo $t\geq2$. Fialová plocha ukáže, jak velká část možných výsledků leží za těmito hranicemi.

::: {.box-data title="Pouze ilustrační příklad — ještě ne krabí výsledek"}
Hodnotu $|t|=2$ jsme zvolili jen proto, aby byly vybarvené plochy dobře vidět. **Není to naše krabí t-statistika ani výsledek krabího modelu.** Obrázek zatím slouží pouze k porovnání ploch při různých stupních volnosti. Co přesně tato plocha znamená, pojmenujeme v následující části.
:::

```{r}
#| label: pripravit-a-zobrazit-t-krivky
#| echo: false
#| output: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 7
#| fig-cap: >
#|   Ilustrační hranice t = -2 a t = 2 nejsou výsledkem krabího
#|   modelu. Ukazují pouze, že při menším df zůstává za stejnými
#|   hranicemi větší část celkové plochy.
#| fig-alt: >
#|   Tři t-křivky pod sebou pro různé stupně volnosti s fialově
#|   vybarvenými chvosty za hranicemi minus 2 a plus 2 a popiskem plochy.
data_t_krivky <-
  tidyr::expand_grid(
    t = seq(
      from = -5,
      to = 5,
      length.out = 1000
    ),
    df = c(
      1,
      stupne_volnosti,
      10000
    )
  ) |>
  dplyr::mutate(
    hustota = dt(
      x = t,
      df = df
    ),
    krajni_oblast = abs(t) >= 2,
    popisek = factor(
      paste0("df = ", df),
      levels = c(
        "df = 1",
        paste0("df = ", stupne_volnosti),
        "df = 10000"
      )
    )
  )

data_t_plochy <-
  data.frame(
    df = c(1, stupne_volnosti, 10000)
  ) |>
  dplyr::mutate(
    popisek = factor(
      paste0("df = ", df),
      levels = c(
        "df = 1",
        paste0("df = ", stupne_volnosti),
        "df = 10000"
      )
    ),
    podil = 2 * pt(
      q = -2,
      df = df
    ),
    text = paste0(
      "plocha = ",
      formatC(
        podil,
        format = "f",
        digits = 3,
        decimal.mark = ","
      )
    )
  )

ggplot2::ggplot(
  data = data_t_krivky,
  mapping = ggplot2::aes(
    x = t,
    y = hustota
  )
) +
  ggplot2::geom_area(
    data = data_t_krivky |>
      dplyr::filter(krajni_oblast),
    mapping = ggplot2::aes(group = t < 0),
    fill = biostat_cols[["amethyst"]],
    alpha = 0.8
  ) +
  ggplot2::geom_line(
    color = biostat_cols[["graphite"]],
    linewidth = 0.9
  ) +
  ggplot2::geom_vline(
    xintercept = c(-2, 2),
    color = biostat_cols[["graphite"]],
    linetype = "dashed",
    linewidth = 0.8
  ) +
  ggplot2::geom_text(
    data = data_t_plochy,
    mapping = ggplot2::aes(
      x = 0,
      y = 0.035,
      label = text
    ),
    color = biostat_cols[["indigo_velvet"]],
    fontface = "bold",
    size = 3.8
  ) +
  ggplot2::facet_wrap(
    facets = ggplot2::vars(popisek),
    ncol = 1
  ) +
  ggplot2::labs(
    x = "Hodnota t",
    y = "Hustota"
  ) +
  theme_biostat(base_size = 12)
```

::: {.box-question title="Přečtěte plochu pod křivkou"}

- Co představuje celá plocha pod jednou křivkou?
- Co představují dvě fialové části?
- U které křivky zabírají fialové části největší podíl celku?
- Co se s fialovou plochou mění, když roste počet stupňů volnosti?
:::

Při malém $df$ zůstává dál od nuly větší část celkové plochy. Pro tento tvar mají statistici běžný název: říkají, že taková křivka má **těžší chvosty**. Velké absolutní hodnoty $t$ proto nejsou při malém $df$ tak neobvyklé. Krabí model má `r stupne_volnosti` stupňů volnosti, takže v následující části použijeme tvar prostřední křivky — už ale s naší skutečnou hodnotou $t$, nikoli s ilustrační hranicí 2.

## Jak neobvyklé by bylo naše t v nulovém světě?

Vraťme se ke krabům. `r render_glossary_term("nulova-hypoteza", display = "Nulová hypotéza")` popisuje svět, ve kterém je `r render_glossary_term("sklon", display = "sklon")` opravdu nulový. T-křivka s `r render_glossary_term("stupne-volnosti", display = "stupni volnosti")` `r stupne_volnosti` ukazuje, jaké hodnoty $t$ by v takovém světě vznikaly při opakování studie.

Pozorovaný sklon leží přibližně `r format_cz(t_pozorovane, digits = 2)` `r render_glossary_term("standardni-chyba", display = "standardních chyb")` od nuly. Odpovídající t-statistika je tedy $t\approx `r format_cz_math(t_pozorovane, digits = 2)`$. Teď se ptáme, **jaká část nulového `r render_glossary_term("rozdeleni", display = "referenčního rozdělení")` leží alespoň tak daleko od nuly jako náš výsledek**.

Celá plocha pod t-křivkou představuje všechny možné hodnoty $t$ v nulovém světě a má velikost 1. Hledaná část této plochy proto musí být číslo mezi 0 a 1.

| Co hledáme? | V krabím modelu |
|---|---|
| Hranice stejně vzdálené od nuly | $t=-`r format_cz_math(abs(t_pozorovane), digits = 2)`$ a $t=+`r format_cz_math(abs(t_pozorovane), digits = 2)`$ |
| Výsledky ještě vzdálenější | vše nalevo od záporné hranice a napravo od kladné hranice |
| Proč oba směry? | alternativa připouští kladný i záporný vztah |

Obě krajní oblasti se nazývají **chvosty rozdělení**. Protože naše otázka připouští kladný i záporný sklon, budeme sčítat plochu v obou chvostech.

```{r}
#| label: pripravit-p-hodnotu
#| include: false
p_hodnota <-
  table_koeficienty[
    "zemepisna_sirka",
    "Pr(>|t|)"
  ]

data_p_krivka <-
  data.frame(
    t = seq(
      from = -7,
      to = 7,
      length.out = 5000
    )
  ) |>
  dplyr::mutate(
    hustota = dt(
      x = t,
      df = stupne_volnosti
    ),
    chvost = abs(t) >= abs(t_pozorovane),
    strana = dplyr::if_else(t < 0, "levá", "pravá")
  )
```

```{r}
#| label: zobrazit-p-hodnotu
#| echo: false
#| output: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 4.8
#| fig-cap: >
#|   Oranžové pásy označují hodnoty t alespoň tak vzdálené od
#|   nuly jako pozorované t. Hledaný podíl tvoří plocha pod
#|   křivkou uvnitř obou pásů.
#| fig-alt: >
#|   T-křivka s oranžovými čárkovanými hranicemi v minus a plus
#|   pozorovaném t a světle fialovými pásy za nimi; plocha v pásech je velmi malá.
ggplot2::ggplot(
  data = data_p_krivka,
  mapping = ggplot2::aes(
    x = t,
    y = hustota
  )
) +
  ggplot2::annotate(
    geom = "rect",
    xmin = c(
      -7,
      abs(t_pozorovane)
    ),
    xmax = c(
      -abs(t_pozorovane),
      7
    ),
    ymin = -Inf,
    ymax = Inf,
    fill = biostat_cols[["orange"]],
    alpha = 0.12
  ) +
  ggplot2::geom_area(
    data = data_p_krivka |>
      dplyr::filter(chvost),
    mapping = ggplot2::aes(group = strana),
    fill = biostat_cols[["orange"]],
    alpha = 0.9
  ) +
  ggplot2::geom_line(
    color = biostat_cols[["graphite"]],
    linewidth = 0.9
  ) +
  ggplot2::geom_vline(
    xintercept = c(
      -abs(t_pozorovane),
      abs(t_pozorovane)
    ),
    color = biostat_cols[["orange"]],
    linetype = "dashed",
    linewidth = 0.9
  ) +
  ggplot2::annotate(
    geom = "text",
    x = c(
      -abs(t_pozorovane),
      abs(t_pozorovane)
    ),
    y = 0.035,
    label = c(
      "levý chvost",
      "pravý chvost"
    ),
    color = biostat_cols[["indigo_velvet"]],
    size = 3.4
  ) +
  ggplot2::labs(
    x = "Hodnota t při nulovém sklonu",
    y = "Hustota"
  ) +
  theme_biostat(base_size = 12)
```

::: {.box-question title="Najděte obě hledané plochy"}
Na obrázku najděte:

1. pozorovanou kladnou t-statistiku;
2. stejně vzdálenou zápornou hodnotu;
3. všechny hodnoty alespoň tak extrémní vpravo;
4. všechny hodnoty alespoň tak extrémní vlevo.

Proč by vybarvení jen pravého chvostu neodpovídalo otázce, která připouští růst i pokles?
:::

Sečteme plochu v levém a pravém chvostu. **Tomuto součtu se ve statistice říká `r render_glossary_term("p-hodnota", display = "p-hodnota")`.**

::: {.box-rule title="P-hodnota v jedné větě"}
Pokud by nulová hypotéza a předpoklady modelu platily, p-hodnota je pravděpodobnost, že získáme t-statistiku alespoň tak vzdálenou od nuly, jako je ta naše. Je to plocha pod referenční křivkou, a proto nabývá hodnot **od 0 do 1**.
:::

Co z p-hodnoty můžeme přímo přečíst?

- **Vždy leží mezi 0 a 1.** Nula by znamenala žádnou plochu, jednička celou plochu pod křivkou.
- **Je podmíněná nulovým světem a modelem.** Nejde o pravděpodobnost, že je nulová hypotéza pravdivá.
- **U oboustranné otázky zahrnuje oba chvosty.** Počítáme kladné i záporné výsledky alespoň tak vzdálené od nuly.
- **Čím menší je p-hodnota, tím méně jsou data slučitelná s přesně nulovým sklonem.** Takto krajní výsledek by v nulovém světě vznikal méně často.

Pro krabí model je tento součet přibližně `r format_cz(p_hodnota, digits = 6)`. Kdyby byl skutečný sklon nulový a model správně popisoval vznik dat, takto krajní nebo ještě krajnější t-statistika by se objevila přibližně `r format_cz(p_hodnota * 100000, digits = 1)`krát ze 100 000 opakování.

### Tři studenti čtou stejný výsledek

> **Anna:** „P-hodnota `r format_cz(p_hodnota, digits = 6)` znamená, že pravděpodobnost nulové hypotézy je `r format_cz(p_hodnota, digits = 6)`.“
>
> **Boris:** „Tak malá p-hodnota znamená, že biologický `r render_glossary_term("efekt", display = "efekt")` musí být obrovský.“
>
> **Cyril:** „Kdyby byl skutečný sklon nulový a platil náš model, takto krajní výsledek by byl velmi vzácný.“

::: {.box-question title="Kdo výsledek čte správně?"}
Který student je nejblíže správné interpretaci? U zbývajících dvou přesně pojmenujte chybu. Co musí i správná věta přiznat o předpokladech modelu?
:::

Cyril jako jediný odpovídá na otázku, kterou p-hodnota skutečně klade.

::: {.box-warning title="Co p-hodnota neříká"}
Následující výklady jsou nesprávné. P-hodnota není:

- ~~pravděpodobnost, že $H_0$ je pravdivá~~;
- ~~pravděpodobnost, že alternativní hypotéza je pravdivá~~;
- ~~velikost biologického efektu~~;
- ~~pravděpodobnost, že výsledek vznikl „náhodou“~~;
- ~~pravděpodobnost úspěšného zopakování studie~~;
- ~~kontrola správnosti všech předpokladů~~;
- ~~konečný verdikt o Bergmannově pravidle~~.
:::

::: {.box-data collapse="true" title="Doplňující: jedna předem položená otázka a mnoho vyzkoušených otázek"}

Naše p-hodnota odpovídá jedné předem pojmenované otázce o sklonu vůči zeměpisné šířce. Kdybychom však ze stejných dat bez předem daného plánu položili mnoho souvisejících otázek, zvýšila by se šance, že alespoň jedna analýza vytvoří nápadně malou p-hodnotu i bez přesvědčivého biologického signálu.

Proto k interpretaci patří také informace, kolik otázek výsledku předcházelo a zda byly zvoleny před pohledem na výsledky. Předem formulovaná biologická otázka omezuje hledání výsledku, který pouze vypadá zajímavě. V následující lekci uvidíme, jak se stejné riziko objeví, když z jedněch dat položíme několik souvisejících otázek.
:::

## Když musíme udělat rozhodnutí

```{r}
#| label: zvolit-hladinu-vyznamnosti
#| include: false
alfa <- 0.05
```

`r render_glossary_term("p-hodnota", display = "P-hodnota")` popisuje, nakolik jsou data slučitelná s `r render_glossary_term("nulova-hypoteza", display = "nulovou hypotézou")`. V některých situacích ale musíme udělat také **rozhodnutí**: nulovou hypotézu zamítnout, nebo ji nezamítnout. K tomu potřebujeme pravidlo zvolené ještě před pohledem na výsledek.

Představme si mnoho studií, ve kterých nulová hypotéza opravdu platí. Jak často jsme ochotni připustit, že ji naše pravidlo přesto omylem zamítne? Hodnota 5 % například znamená, že při mnoha opakovaných studiích se stejným pravidlem připouštíme pětiprocentní dlouhodobé riziko falešného zamítnutí pravdivé nulové hypotézy.

Tato předem zvolená hranice se nazývá `r render_glossary_term("hladina-vyznamnosti", display = "hladina významnosti")` a označuje se řeckým písmenem $\alpha$ (alfa). Často se používá:

- 5 %, tedy $\alpha=`r format_cz_math(alfa)`$;
- v přísněji nastavených analýzách například 1 %, tedy $\alpha=0{,}01$.

Ani 5 %, ani 1 % není přírodní zákon. Hranice má vycházet z pravidel oboru a z následků nesprávného rozhodnutí. Dva možné způsoby, jak se při takovém rozhodování splést, pojmenujeme v následující části.

### Kritické hodnoty

```{r}
#| label: pripravit-kriticke-hodnoty
#| include: false
kriticka_hodnota <-
  qt(
    p = 1 - alfa / 2,
    df = stupne_volnosti
  )
```

Plochu pod t-křivkou ani její přesné hranice není snadné vypočítat ručně. Proto statistické učebnice před běžným používáním počítačů obsahovaly několik stran předem vypočtených tabulek. Z nich se vyhledávala `r render_glossary_term("kriticka-hodnota", display = "kritická hodnota")` pro zvolené `r render_glossary_term("rozdeleni", display = "rozdělení")`, hladinu významnosti a `r render_glossary_term("stupne-volnosti", display = "stupně volnosti")`.

Tabulka tedy nepočítala přesnou plochu znovu pro každý model. Umožnila jednodušší postup: najít správnou hranici a porovnat s ní pozorovanou testovou statistiku. Pokud $|t|$ leželo za kritickou hodnotou, odpovídalo to p-hodnotě nejvýše $\alpha$ a nulová hypotéza se zamítla. Podrobnější tabulky také dovolovaly p-hodnotu zařadit mezi dvě uvedené hranice, ale obvykle nedaly tolik desetinných míst jako dnešní software.

Nebyla to jedna hranice pro všechny modely. Bylo potřeba zvolit:

- referenční rozdělení — zde t-rozdělení;
- oboustrannou hladinu významnosti — zde $\alpha=`r format_cz_math(alfa)`$;
- správný řádek podle stupňů volnosti — zde $df=`r stupne_volnosti`$.

Malá část takové tabulky by vypadala následovně. Hodnoty říkají, jak daleko od nuly začíná při oboustranné hladině 5 % krajních pět procent t-rozdělení.

| Reziduální stupně volnosti | Kritická hodnota $|t|$ pro $\alpha=`r format_cz_math(alfa)`$ |
|---:|---:|
| **`r stupne_volnosti` — krabí model** | **`r format_cz(kriticka_hodnota, digits = 3)`** |
| 10 000 | `r format_cz(qt(p = 1 - alfa / 2, df = 10000), digits = 3)` |

V řádku pro `r stupne_volnosti` stupňů volnosti najdeme hranice přibližně $-`r format_cz_math(kriticka_hodnota, digits = 3)`$ a $+`r format_cz_math(kriticka_hodnota, digits = 3)`$. Pozorované krabí $|t|=`r format_cz_math(abs(t_pozorovane), digits = 2)`$ leží daleko za nimi.

Teď stejné hranice přidáme přímo do známého t-rozdělení. Fialová vyznačuje předem zvolené kritické oblasti; oranžová ukazuje výsledek našeho modelu.

```{r}
#| label: zobrazit-kriticke-hodnoty-v-t-rozdeleni
#| echo: false
#| output: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 4.8
#| fig-cap: >
#|   Kritické hranice a pozorovaná t-statistika v nulovém
#|   t-rozdělení. Fialové chvosty mají dohromady plochu alfa;
#|   oranžové hranice vymezují oboustrannou p-hodnotu.
#| fig-alt: !expr paste0("Nulové t-rozdělení s fialovými kritickými hranicemi minus a plus ", format_cz(kriticka_hodnota, digits = 3), ", oranžovou pozorovanou hodnotou t rovnou plus ", format_cz(t_pozorovane, digits = 3), " a čárkovanou oranžovou zrcadlovou hranicí.")
kriticke_t_graf <- kriticka_hodnota

data_kriticke_t <-
  tibble::tibble(
    t = seq(
      from = -7,
      to = 7,
      length.out = 5000
    )
  ) |>
  dplyr::mutate(
    hustota = dt(
      x = t,
      df = stupne_volnosti
    ),
    kriticky_chvost = abs(t) >= kriticke_t_graf,
    strana = dplyr::if_else(t < 0, "levá", "pravá")
  )

ggplot2::ggplot(
  data = data_kriticke_t,
  mapping = ggplot2::aes(
    x = t,
    y = hustota
  )
) +
  ggplot2::geom_area(
    data = data_kriticke_t |>
      dplyr::filter(kriticky_chvost),
    mapping = ggplot2::aes(group = strana),
    fill = biostat_cols[["amethyst"]],
    alpha = 0.42
  ) +
  ggplot2::geom_line(
    color = biostat_cols[["graphite"]],
    linewidth = 0.9
  ) +
  ggplot2::geom_vline(
    xintercept = c(-kriticke_t_graf, kriticke_t_graf),
    color = biostat_cols[["amethyst"]],
    linetype = 2,
    linewidth = 0.9
  ) +
  ggplot2::geom_vline(
    xintercept = t_pozorovane,
    color = biostat_cols[["orange"]],
    linewidth = 1.1
  ) +
  ggplot2::geom_vline(
    xintercept = -abs(t_pozorovane),
    color = biostat_cols[["orange"]],
    linetype = 2,
    linewidth = 1.1
  ) +
  ggplot2::annotate(
    geom = "label",
    x = kriticke_t_graf,
    y = 0.31,
    label = paste0(
      "kritické hranice = ±",
      format_cz(kriticke_t_graf, digits = 3)
    ),
    hjust = -0.05,
    color = biostat_cols[["amethyst"]],
    fill = "white",
    linewidth = 0.2,
    size = 3.2
  ) +
  ggplot2::annotate(
    geom = "label",
    x = t_pozorovane,
    y = 0.12,
    label = paste0(
      "pozorované t = +",
      format_cz(t_pozorovane, digits = 3)
    ),
    hjust = 1.05,
    color = biostat_cols[["orange"]],
    fill = "white",
    linewidth = 0.2,
    size = 3.2
  ) +
  ggplot2::annotate(
    geom = "label",
    x = -abs(t_pozorovane),
    y = 0.12,
    label = paste0(
      "zrcadlová hranice = -",
      format_cz(abs(t_pozorovane), digits = 2)
    ),
    hjust = -0.05,
    color = biostat_cols[["orange"]],
    fill = "white",
    linewidth = 0.2,
    size = 3.2
  ) +
  ggplot2::labs(
    x = "t-statistika v nulovém světě",
    y = "Hustota"
  ) +
  theme_biostat(base_size = 12)
```

Fialové chvosty za kritickými hranicemi mají dohromady plochu $\alpha=`r format_cz_math(alfa)`$. Skutečně pozorovaná t-statistika je kladná; záporná oranžová hranice je pouze její zrcadlový protějšek pro oboustrannou otázku. P-hodnota je plocha pod t-křivkou za oběma oranžovými hranicemi. Protože $|t|=`r format_cz_math(abs(t_pozorovane), digits = 2)`$ leží za kritickou hodnotou `r format_cz(kriticka_hodnota, digits = 3)`, vedou porovnání $|t|$ s kritickou hodnotou a porovnání p-hodnoty s $\alpha$ ke stejnému rozhodnutí.

V R stejnou hodnotu vyhledá funkce `qt()`. Funkce `pt()` pak spočítá plochu za pozorovaným t, tedy p-hodnotu. Následující kód používá přímo model, který jsme už fitovali:

```{r}
#| label: porovnat-kritickou-hodnotu
#| echo: true
#| eval: true
# Hranici spočítáme pro alfa = 0,05 a df našeho modelu.
kriticke_t <-
  qt(
    p = 1 - 0.05 / 2,
    df = df.residual(object = mod_krabi)
  )

# t-statistiku přečteme z tabulky koeficientů modelu.
tabulka_koeficientu <- coef(object = summary(object = mod_krabi))
t_krabi <- tabulka_koeficientu["zemepisna_sirka", "t value"]

# Obě krajní plochy t-rozdělení sečteme.
p_krabi <-
  2 * pt(
    q = -abs(x = t_krabi),
    df = df.residual(object = mod_krabi)
  )

c(
  kriticka_hodnota = kriticke_t,
  t_krabi = t_krabi,
  p_hodnota = p_krabi
)
abs(x = t_krabi) > kriticke_t
```

Při běžné analýze se o ruční vyhledávání kritických hodnot ani o výpočet pomocí `qt()` a `pt()` nemusíme starat. R už test koeficientu provedl při fitování modelu. Podívejme se znovu na známý výpis:

```{r}
#| label: znovu-zobrazit-souhrn-pro-p-hodnotu
#| echo: true
#| eval: true
# Znovu zobrazíme celý souhrn stejného modelu.
summary(object = mod_krabi)
```

P-hodnotu pro naši biologickou otázku najdeme v tabulce `Coefficients`:

- vybereme řádek `zemepisna_sirka`, protože právě ten obsahuje `r render_glossary_term("sklon", display = "sklon")` vztahu k zeměpisné šířce;
- ve sloupci `t value` vidíme pozorovanou t-statistiku `r format_cz(t_pozorovane, digits = 3)`;
- ve sloupci `Pr(>|t|)` vidíme oboustrannou p-hodnotu `r format_cz(p_hodnota, digits = 6)`;
- počet reziduálních stupňů volnosti `r stupne_volnosti` najdeme níže ve výpisu na řádku `Residual standard error`.

To je hodnota, kterou budeme při běžné práci číst. Tabulka kritických hodnot, porovnání s $|t|$ a porovnání p-hodnoty s $\alpha$ jsou jen tři zápisy stejného rozhodovacího pravidla. Počítač princip porovnání nezměnil; stejné hranice a plochy však vypočítá za nás.

::: {.box-data collapse="true" title="Doplňující: proč se mluví o t-testu"}
Test, který R provedl pro sklon, patří mezi takzvané `r render_glossary_term("t-test", display = "t-testy")`: odhad vydělíme jeho standardní chybou a výslednou t-statistiku porovnáme s t-rozdělením. V učebnicích a článcích se pod názvem t-test nejčastěji objevuje porovnání průměrů dvou skupin, například kontroly a zásahu. K tomuto t-testu se můžeme vrátit později v kurzu. Uvidíme, že jde o stejný lineární model, jen s prediktorem, který označuje příslušnost ke skupině.
:::

Rozhodovací pravidlo pak zní:

| Výsledek porovnání | Běžné rozhodnutí a označení |
|---|---|
| p-hodnota je nejvýše $\alpha$ | **zamítneme nulovou hypotézu**; výsledek označíme jako **statisticky významný** |
| p-hodnota je větší než $\alpha$ | **nulovou hypotézu nezamítneme**; výsledek není na zvolené hladině statisticky významný |

V těchto materiálech používáme český termín **statisticky významný**. <!-- TODO: add slug statisticky-vyznamny to pojmy.yaml --> V článcích často uvidíte také synonymum **statisticky signifikantní** (anglicky *statistically significant*). P-hodnota je v obou případech menší než předem zvolená rozhodovací hranice.

::: {.box-warning title="Statisticky významný neznamená biologicky důležitý"}
Slovo „významný“ zde **neříká**, že je `r render_glossary_term("efekt", display = "efekt")` velký, biologicky důležitý nebo určitě pravdivý. O biologickém významu rozhodujeme podle velikosti efektu v jeho jednotkách a podle biologického kontextu, nikoli podle samotného porovnání p-hodnoty s hranicí.
:::

U krabů je p-hodnota `r format_cz(p_hodnota, digits = 6)`. Při $\alpha=`r format_cz_math(alfa)`$ proto nulovou hypotézu o přesně nulovém sklonu zamítáme a výsledek nazveme statisticky významným. Odpovídající `r render_glossary_term("interval-spolehlivosti", display = "95% interval spolehlivosti")` nulu neprotíná, takže vede ke stejnému rozhodnutí. Neříkáme tím, že jsme alternativní hypotézu dokázali; říkáme, že přesně nulový sklon je za předpokladů modelu s výsledkem velmi obtížně slučitelný.

### Stejnou otázku můžeme položit intervalu

Teď už víme, proč porovnáváme právě 95% interval s hladinou $\alpha=`r format_cz_math(alfa)`$. U oboustranného testu patří k sobě:

- kritické hranice vymezující prostředních 95 % příslušného referenčního rozdělení;
- 95% interval spolehlivosti vytvořený pomocí stejných hranic;
- zbývajících 5 % rozdělených mezi oba chvosty;
- rozhodovací hranice $\alpha=`r format_cz_math(alfa)`$.

Stejnou otázku proto můžeme položit bez dalšího výpočtu p-hodnoty: **obsahuje 95% interval spolehlivosti nulový sklon?**

```{r}
#| label: zobrazit-interval-proti-nule
#| echo: false
#| output: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 3.4
#| fig-cap: >
#|   Odhad sklonu a jeho 95% interval spolehlivosti. Svislá čára
#|   označuje nulový sklon, tedy žádnou lineární změnu k severu.
#| fig-alt: >
#|   Fialový interval spolehlivosti sklonu s oranžovým bodem odhadu;
#|   čárkovaná svislá čára nulového sklonu leží vlevo mimo interval.
data_interval_proti_nule <-
  data.frame(
    odhad = odhad_sklonu,
    dolni_mez = interval_sklonu[[1]],
    horni_mez = interval_sklonu[[2]],
    radek = 1
  )

ggplot2::ggplot(
  data = data_interval_proti_nule,
  mapping = ggplot2::aes(
    x = odhad,
    y = radek
  )
) +
  ggplot2::geom_vline(
    xintercept = 0,
    color = biostat_cols[["graphite"]],
    linetype = "dashed",
    linewidth = 0.9
  ) +
  ggplot2::geom_segment(
    mapping = ggplot2::aes(
      x = dolni_mez,
      xend = horni_mez,
      yend = radek
    ),
    color = biostat_cols[["indigo_velvet"]],
    linewidth = 2
  ) +
  ggplot2::geom_point(
    color = biostat_cols[["orange"]],
    size = 4
  ) +
  ggplot2::scale_x_continuous(
    limits = c(-0.1, 0.75)
  ) +
  ggplot2::scale_y_continuous(
    breaks = NULL
  ) +
  ggplot2::labs(
    x = "Sklon (mm na 1° severněji)",
    y = NULL
  ) +
  theme_biostat(base_size = 12)
```

Krabí interval od `r format_cz(interval_sklonu[[1]], digits = 3)` do `r format_cz(interval_sklonu[[2]], digits = 3)` mm na stupeň nulu neprotíná. Při oboustranné hladině $\alpha=`r format_cz_math(alfa)`$ proto vede ke stejnému rozhodnutí jako p-hodnota: nulový sklon zamítneme. Kdyby 95% interval nulu protínal, p-hodnota by při stejném modelu nebyla menší než 0,05 a nulovou hypotézu bychom nezamítli.

Interval nám přitom ukazuje víc než samotné rozhodnutí. Vidíme z něj také rozsah kladných velikostí vztahu, které zůstávají slučitelné s daty a modelem.

**Nezamítnout nulovou hypotézu naopak neznamená dokázat, že efekt je přesně nula.**

::: {.box-data collapse="true" title="Doplňující: nula nemusí být biologicky nejdůležitější hranice"}

Nulový sklon je matematicky jednoduchá referenční hodnota, ale biologická otázka může mít důležitější hranici. Výzkumníci mohou například předem určit, jak velký kladný sklon šířky krunýře na jeden stupeň je ještě biologicky zanedbatelný a od jaké velikosti už by měnil jejich biologický závěr. Konkrétní hranice musí vycházet z biologických znalostí a měla by být zdůvodněna před čtením výsledku.

Interval spolehlivosti můžeme porovnat nejen s nulou, ale také s takovou předem zvolenou hranicí. Může například celý ležet pod hranicí biologicky důležitého kladného sklonu, i když obsahuje nulu; v takovém případě mezi hodnotami slučitelnými s daty a modelem není velký kladný efekt, ale interval neříká automaticky totéž o opačném směru. Jindy může interval připouštět nulový i biologicky významný kladný efekt. Formální testy ekvivalence patří do pokročilejšího kurzu; zde stačí číst celý rozsah efektů slučitelných s daty.
:::

## Dva způsoby, jak se splést

Rozhodovací pravidlo už máme: `r render_glossary_term("p-hodnota", display = "p-hodnotu")` porovnáme s předem zvolenou `r render_glossary_term("hladina-vyznamnosti", display = "hladinou významnosti")`. Rozhodnutí ale není totéž co jistota. Může odpovídat skutečnosti, nebo se můžeme splést dvěma různými způsoby.

### Rychlá paměťová pomůcka

U těhotenského testu si představme nulový stav „těhotenství není přítomno“:

- **chyba I. typu neboli falešně pozitivní výsledek:** test hlásí těhotenství v očividně nemožné situaci;
- **chyba II. typu neboli falešně negativní výsledek:** test nehlásí těhotenství v očividně pokročilé fázi.

::: {#fig-tehotensky-test-chyby layout-ncol="2"}
![Falešně pozitivní výsledek: zjevně netěhotný muž překvapeně drží pozitivní těhotenský test.](images/chyba_i_falesne_pozitivni.jpg)

![Falešně negativní výsledek: zjevně těhotná žena překvapeně drží negativní těhotenský test.](images/chyba_ii_falesne_negativni.jpg)

Dvě úmyslně absurdní situace jako paměťová pomůcka: vlevo chyba I. typu, vpravo chyba II. typu.
:::

Je to úmyslně absurdní zkratka. Skutečná definice vždy závisí na tom, jak jsme určili nulový stav.

- **Chyba I. typu:** zamítneme `r render_glossary_term("nulova-hypoteza", display = "nulovou hypotézu")`, přestože je pravdivá. Jde o falešně pozitivní závěr neboli falešný poplach. Předem zvolená hladina významnosti $\alpha$ určuje dlouhodobé riziko právě této chyby.
- **Chyba II. typu:** nulovou hypotézu nezamítneme, přestože ve skutečnosti neplatí a existuje skutečný `r render_glossary_term("efekt", display = "efekt")`.

### Biologický příklad: invazní druh

Představme si průzkum, jehož nulová hypotéza zní:

$$
H_0:\ \text{invazní druh v oblasti není přítomen}.
$$

| Skutečnost | Průzkum hlásí „nepřítomen“ | Průzkum hlásí „přítomen“ |
|---|---|---|
| Druh není přítomen | správné rozhodnutí | **chyba I. typu:** zbytečný zásah |
| Druh je přítomen | **chyba II. typu:** přehlédnutá invaze | správné rozhodnutí |

::: {.box-question title="Rozhodněte — Která chyba je horší?"}

1. Zařaďte oba těhotenské příklady jako chybu I. nebo II. typu.
2. Totéž udělejte pro invazní druh.
3. Která chyba může být v ochraně přírody závažnější?
4. Jak odpověď závisí na ceně zásahu, rychlosti šíření a vratnosti škod?
:::

::: {.box-result collapse="true" title="Možné odpovědi"}
Pozitivní test u muže je chyba I. typu, negativní test v pokročilém těhotenství chyba II. typu. U invazního druhu je zbytečný zásah chybou I. typu a přehlédnutá invaze chybou II. typu. Přehlédnutá invaze bývá závažnější, když se druh rychle šíří a škody jsou nevratné; levný a vratný zásah naopak snese více falešných poplachů.
:::

Žádný typ chyby není vždy horší. Záleží na následcích.

U krabího modelu by chyba I. typu znamenala, že tvrdíme vztah, přestože je populační `r render_glossary_term("sklon", display = "sklon")` nulový. Chyba II. typu by znamenala, že skutečný vztah neodhalíme.

### Když skutečný vztah přehlédneme

```{r}
#| label: spocitat-ilustracni-silu
#| include: false
pocet_simulaci_sila <- 5000
skutecny_sklon_sila <- 0.25
alfa_sila <- alfa
stredni_velikost <- mean(data_lokality$prumerna_sirka_mm)

# Stejné náhodné residuum použijeme pro oba scénáře.
set.seed(900723)
vec_odhalen_vztah <-
  simulovat_odhaleni_sklonu(
    vec_zemepisna_sirka = data_lokality$zemepisna_sirka,
    stredni_velikost = stredni_velikost,
    skutecny_sklon = skutecny_sklon_sila,
    rezidualni_sd = rezidualni_sd,
    pocet_simulaci = pocet_simulaci_sila,
    alfa = alfa_sila
  )

pocet_odhalenych <- sum(vec_odhalen_vztah)
odhad_sily <- mean(vec_odhalen_vztah)

# Stejný výpočet zopakujeme pro sklon odhadnutý z krabích dat.
set.seed(900723)
vec_odhalen_pozorovany_sklon <-
  simulovat_odhaleni_sklonu(
    vec_zemepisna_sirka = data_lokality$zemepisna_sirka,
    stredni_velikost = stredni_velikost,
    skutecny_sklon = odhad_sklonu,
    rezidualni_sd = rezidualni_sd,
    pocet_simulaci = pocet_simulaci_sila,
    alfa = alfa_sila
  )
pocet_odhalenych_pozorovany_sklon <-
  sum(vec_odhalen_pozorovany_sklon)
odhad_sily_pozorovany_sklon <-
  mean(vec_odhalen_pozorovany_sklon)
```

Chyba II. typu popisuje jednotlivý omyl. Zbývá ale praktická otázka: **Jak dobře je naše studie připravena skutečný vztah vůbec odhalit?** Této schopnosti říkáme **síla testu**. <!-- TODO: add slug sila-testu to pojmy.yaml -->

V tomto kontextu znamená „efekt“ vždy **předem zvolenou skutečnou velikost efektu** (*effect size*) v širší populaci — například skutečný sklon `r format_cz(skutecny_sklon_sila)` mm na stupeň. Nejde o pravděpodobnost, že efekt existuje, ani o dodatečné přejmenování sklonu, který jsme právě pozorovali.

Kdybychom opakovali stejně navrženou studii ve světě s touto konkrétní skutečnou velikostí efektu, síla testu by byla podíl opakování, ve kterých bychom nulovou hypotézu správně zamítli. Vyšší síle obvykle pomáhá:

- větší skutečná velikost efektu;
- více nezávislých lokalit;
- menší `r render_glossary_term("variabilita", display = "variabilita")` a tím menší `r render_glossary_term("standardni-chyba", display = "standardní chyba")`;
- méně přísné rozhodovací pravidlo — za cenu vyššího rizika chyby I. typu.

Nízká síla znamená, že častěji nastane chyba II. typu: reálný vztah přehlédneme. Proto **nezamítnutí nulové hypotézy není důkazem nulového efektu**, zvlášť u malé nebo velmi variabilní studie.

### Kolik studií vztah odhalí?

Sílu testu můžeme `r render_glossary_term("odhad", display = "odhadnout")` simulací. Musíme si zvolit **konkrétní skutečný efekt**, který chceme umět odhalit. Zvolíme vztah `r format_cz(skutecny_sklon_sila)` mm na jeden stupeň, tedy asi `r format_cz(100 * skutecny_sklon_sila / odhad_sklonu, digits = 0)` % sklonu pozorovaného u krabů. Ptáme se tak, zda by stejná studie dokázala zachytit i slabší, ale předem zvolený vztah. Výsledek porovnáme také s naším pozorovaným sklonem.

Vytvoříme `r format(pocet_simulaci_sila, big.mark = " ")` hypotetických studií. Každá bude mít stejné zeměpisné šířky jako naše data. Do každé vložíme skutečný sklon `r format_cz(skutecny_sklon_sila)` mm na stupeň a přidáme náhodná `r render_glossary_term("residuum", display = "residua")` podobně velká jako v našem modelu. Potom vždy fitujeme přímku a zaznamenáme, zda při $\alpha=`r format_cz_math(alfa_sila)`$ zamítla nulový sklon.

| Co jsme v příkladu nastavili | Hodnota |
|---|---:|
| Skutečný sklon, který chceme odhalit | `r format_cz(skutecny_sklon_sila)` mm na 1° |
| Počet lokalit v každé studii | `r pocet_lokalit` |
| Velikost náhodných residuí | přibližně `r format_cz(rezidualni_sd, digits = 2)` mm |
| Rozhodovací hranice | $\alpha=`r format_cz_math(alfa_sila)`$ |
| Studie, které nulový sklon zamítly | `r format(pocet_odhalenych, big.mark = " ")` z `r format(pocet_simulaci_sila, big.mark = " ")` |

Sílu v tomto **ilustračním scénáři** odhadneme jako podíl:

$$
\text{odhad síly testu}
=
\frac{\text{počet studií, které zamítly nulový sklon}}
{\text{počet všech simulovaných studií}}.
$$

$$
\frac{`r pocet_odhalenych`}{`r pocet_simulaci_sila`}
= `r format_cz_math(odhad_sily, digits = 3)`
\approx `r format_cz_math(100 * odhad_sily, digits = 1)`\,\%.
$$

Zkráceně bychom stejný výpočet mohli zapsat jako

$$
\widehat{\text{síla}}=\frac{m}{N},
$$

kde $m$ je počet simulovaných studií, které nulový sklon zamítly, a $N$ je počet všech simulovaných studií.

To znamená: kdyby byl skutečný sklon `r format_cz(skutecny_sklon_sila)` mm na stupeň a všechny ostatní podmínky odpovídaly této simulaci, přibližně `r format_cz(100 * odhad_sily, digits = 1)` % stejně navržených studií by při $\alpha=`r format_cz_math(alfa)`$ nulový sklon zamítlo. Zbylé studie by tento skutečný vztah přehlédly — dopustily by se chyby II. typu.

::: {.box-result title="Druhý zvolený scénář: větší skutečný sklon"}
Zvolme ještě druhý scénář, ve kterém je skutečný sklon stejně velký jako náš odhad, tedy `r format_cz(odhad_sklonu, digits = 3)` mm na stupeň. Nulový sklon pak zamítne `r format(pocet_odhalenych_pozorovany_sklon, big.mark = " ")` z `r format(pocet_simulaci_sila, big.mark = " ")` studií a síla vychází přibližně `r format_cz(100 * odhad_sily_pozorovany_sklon, digits = 2)` %. Je vyšší než u sklonu `r format_cz(skutecny_sklon_sila)` mm na stupeň, protože větší skutečný efekt je při stejném počtu lokalit a stejné variabilitě snazší odhalit.

**Ani toto číslo není síla naší již dokončené studie.** Pozorovaný sklon je jen odhad a skutečný sklon může být menší i větší. Sílu proto počítáme pro předem zvolené scénáře, nikoli zpětně z výsledku, který jsme už viděli.
:::

::: {.box-question title="Změňte scénář v hlavě"}
Co by se podle vás stalo se silou, kdyby byl skutečný sklon menší? Co kdybychom měli více nezávislých lokalit nebo menší residua?
:::

::: {.box-result collapse="true" title="Možná odpověď"}
Menší skutečný sklon by se odhaloval hůř, takže síla by klesla. Více nezávislých lokalit nebo menší residua zmenší standardní chybu sklonu, a síla proto vzroste.
:::

Těchto `r format_cz(100 * odhad_sily, digits = 1)` % **není pravděpodobnost, že je náš závěr pravdivý**, ani dodatečné hodnocení pozorované p-hodnoty. Je to vlastnost předem popsaného hypotetického scénáře. Formální plánování síly a potřebné velikosti `r render_glossary_term("vyber", display = "výběru")` necháme na pozdější kurz.

## Co tedy můžeme říct o krabech?

Vraťme se k původní otázce.

**Dokázali jsme, že severněji žijící krabi jsou větší?**

Ne. Máme silný a dobře popsaný výsledek pro `r pocet_lokalit` sledovaných slanisek a jednoduchý model. To není totéž jako definitivní důkaz obecného biologického pravidla.

::: {.box-question title="Která výsledková věta vám řekne více?"}
Porovnejte dva závěry:

> **A:** „$p<0{,}05$, proto jsou severní krabi větší.“
>
> **B:** „`r render_glossary_term("odhad", display = "Odhad")` `r render_glossary_term("sklon", display = "sklonu")` říká, že `r render_glossary_term("prumer", display = "průměrná šířka krunýře")` rostla přibližně o `r format_cz(odhad_sklonu, digits = 3)` mm na stupeň; `r render_glossary_term("interval-spolehlivosti", display = "interval")` od `r format_cz(interval_sklonu[[1]], digits = 3)` do `r format_cz(interval_sklonu[[2]], digits = 3)` mm ukazuje, jak velký vztah je s daty slučitelný.“

Která věta vám o krabech řekne více? Co je špatně ve větě A? Co ještě chybí i ve větě B?
:::

::: {.box-result collapse="true" title="Možná odpověď"}
Více řekne věta B: obsahuje velikost vztahu v jednotkách i jeho `r render_glossary_term("nejistota", display = "nejistotu")`. Věta A zaměňuje statistickou významnost za důkaz a vůbec neříká, o kolik jsou krabi větší. Větě B ještě chybí biologický význam změny, srovnání s nulovým sklonem a omezení studie.
:::

Úplný závěr můžeme poskládat z následujících odpovědí:

| Nejprve se ptáme | Odpověď v krabím příběhu |
|---|---|
| Jak velký je vztah? | `r format_cz(odhad_sklonu, digits = 3)` mm na 1°; `r format_cz(efekt_zvoleny_rozdil, digits = 2)` mm na `r rozdil_stupnu`° |
| Je biologicky zajímavý? | Přes sledovaný rozsah odpovídá asi `r format_cz(efekt_cely_rozsah, digits = 2)` mm |
| Jak přesný je odhad? | 95% CI [`r format_cz(interval_sklonu[[1]], digits = 3)`; `r format_cz(interval_sklonu[[2]], digits = 3)`] mm na 1° |
| Jak dobře se data shodují s nulovým sklonem? | $t\approx `r format_cz_math(t_pozorovane, digits = 2)`$, $df=`r stupne_volnosti`$, $p\approx `r format_cz_math(p_hodnota, digits = 6)`$ |
| Jaké je rozhodnutí při $\alpha=`r format_cz_math(alfa)`$? | výsledek je statisticky významný; `r render_glossary_term("nulova-hypoteza", display = "nulovou hypotézu")` zamítáme |
| Co ještě může závěr omezovat? | `r render_glossary_term("vyber", display = "výběr")` lokalit, prostorová závislost, agregace a jednoduchost modelu |

::: {.box-question title="Vaše verze — Vysvětlete výsledek biologovi"}
Napište pět až šest vět tak, aby čtenář nejprve pochopil, **co se s velikostí krabů mění a o kolik**. Potom přidejte nejistotu, srovnání s nulovým sklonem a nakonec omezení toho, co lze z těchto slanisek zobecnit.
:::

Jedna možná odpověď:

::: {.box-result title="Příklad úplného výsledku"}
Na `r pocet_lokalit` sledovaných slaniskách atlantského pobřeží se průměrná šířka krunýře samců zvyšovala přibližně o `r format_cz(odhad_sklonu, digits = 3)` mm na každý stupeň zeměpisné šířky, tedy asi o `r format_cz(efekt_zvoleny_rozdil, digits = 2)` mm na `r rozdil_stupnu`°. Přes celý sledovaný rozsah odpovídá fitovaná přímka změně přibližně `r format_cz(efekt_cely_rozsah, digits = 2)` mm, která je vzhledem k velikosti krabů biologicky zajímavá. 95% interval spolehlivosti sahal od `r format_cz(interval_sklonu[[1]], digits = 3)` do `r format_cz(interval_sklonu[[2]], digits = 3)` mm na stupeň. Při hladině $\alpha=`r format_cz_math(alfa)`$ byl výsledek statisticky významný; při předpokladech jednoduchého `r render_glossary_term("linearni-model", display = "lineárního modelu")` byla data obtížně slučitelná s přesně nulovým sklonem ($t\approx `r format_cz_math(t_pozorovane, digits = 2)`$, $df=`r stupne_volnosti`$, $p\approx `r format_cz_math(p_hodnota, digits = 6)`$). Výsledek podporuje kladný vztah mezi zeměpisnou šířkou a průměrnou velikostí samců na sledovaných slaniskách, ale sám neprokazuje, že zeměpisná šířka velikost krabů způsobuje.
:::

```{r}
#| label: overit-vliv-jednotlivych-lokalit
#| include: false
vec_sklony_bez_lokality <-
  purrr::map_dbl(
    .x = seq_len(pocet_lokalit),
    .f = ~ {
      mod_bez_lokality <-
        lm(
          formula = prumerna_sirka_mm ~ zemepisna_sirka,
          data = data_lokality[-.x, ]
        )
      coef(mod_bez_lokality)[["zemepisna_sirka"]]
    }
  )
```

### Co se do jedné přímky nevešlo

`r render_glossary_term("diagnostika-modelu", display = "Diagnostické grafy")` jsme prohlédli hned po fitování modelu. Ani dobrý diagnostický obraz ale neřeší všechna omezení studie a zvoleného modelu:

- **Jsou všechna slaniska stejně dobře změřena?** Ne úplně. Počty krabů se mírně liší, takže průměry nemají stejnou přesnost. Jednoduchá nevážená regrese jim přesto dává stejnou váhu.
- **Co se stalo s rozdíly mezi jednotlivými kraby?** Zprůměrováním jsme je z hlavní analýzy odstranili. Pro naši otázku to bylo užitečné, ale část informace jsme ztratili.
- **Nestojí přímka hlavně na několika lokalitách?** Když postupně vynecháme vždy jednu lokalitu, sklon se nejvíce změní po vynechání lokalit `r paste(data_lokality$kod_lokality[order(abs(vec_sklony_bez_lokality - odhad_sklonu), decreasing = TRUE)[1:2]], collapse = " a ")`. Všechny sklony však zůstávají kladné a leží přibližně mezi `r format_cz(min(vec_sklony_bez_lokality), digits = 3)` a `r format_cz(max(vec_sklony_bez_lokality), digits = 3)` mm na stupeň.
- **Mohou si blízká slaniska být podobnější?** Ano. Geografické lokality nemusí být zcela nezávislé a přímka není úplný prostorový ekologický model.
- **Dokazuje vztah příčinu?** Ne. Studie je observační. Zeměpisná šířka může souviset s teplotou, délkou sezóny i dalšími podmínkami.

Kvůli těmto omezením nemusíme výsledek zahodit. Musíme ale být opatrní, když ho zobecňujeme na jiná místa, populace nebo biologické mechanismy.

## Měsíc pořád není ze sýra

`r render_glossary_term("p-hodnota", display = "P-hodnota")` nám neřekne pravděpodobnost, že je Měsíc ze sýra. Stejně tak nám neřekne pravděpodobnost, že je Bergmannovo pravidlo pravdivé.

Odpovídá na mnohem užší otázku: **Jak neobvyklá by byla naše testová statistika, kdyby platilo přesně zvolené nulové tvrzení a předpoklady modelu?**

Vědecký závěr vzniká spojením více částí:

- velikosti `r render_glossary_term("efekt", display = "efektu")`;
- biologického významu;
- `r render_glossary_term("nejistota", display = "nejistoty")`;
- slučitelnosti s konkrétním nulovým tvrzením;
- uspořádání studie;
- předpokladů a možných chyb;
- dosavadních znalostí a dalších studií.

::: {.box-rule title="Hlavní myšlenka lekce"}
**Věda nezískává pravdu z $p<0{,}05$. Staví opatrné a opravitelné závěry z efektu, nejistoty, kvality studie a souvislostí.**
:::

## Co bude následovat dál

Zeměpisná šířka byla číselný `r render_glossary_term("prediktor", display = "prediktor")`, a proto jsme vztah popisovali `r render_glossary_term("sklon", display = "sklonem")`. V následující lekci nebude prediktorem číslo, ale příslušnost ke skupině, například kontrola a zásah. Místo sklonu se budeme ptát na rozdíl mezi dvěma `r render_glossary_term("prumer", display = "průměry")`.

Způsob uvažování zůstane stejný: nejprve velikost a biologický význam rozdílu, potom `r render_glossary_term("nejistota", display = "nejistota")`, srovnání s přesnou nulou a nakonec omezení studie.

## Shrnutí

- Testovatelnost tvrzení není totéž jako jeho věrohodnost.
- Nulový `r render_glossary_term("sklon", display = "sklon")` znamená žádný **systematický lineární trend**, nikoli žádnou náhodnou `r render_glossary_term("variabilita", display = "variabilitu")` mezi lokalitami.
- T-statistika měří vzdálenost od nulového sklonu ve `r render_glossary_term("standardni-chyba", display = "standardních chybách")`.
- `r render_glossary_term("stupne-volnosti", display = "Stupně volnosti")` určují, kterou referenční t-křivku použijeme a kolik plochy zůstává daleko od nuly.
- `r render_glossary_term("p-hodnota", display = "P-hodnota")` je podmíněná plocha v chvostech referenčního `r render_glossary_term("rozdeleni", display = "rozdělení")`, nikoli pravděpodobnost pravdivosti hypotézy.
- „Statisticky významný“ je výsledek rozhodovacího pravidla; neznamená automaticky velký nebo biologicky důležitý `r render_glossary_term("efekt", display = "efekt")`.
- Zamítnutí i nezamítnutí `r render_glossary_term("nulova-hypoteza", display = "nulové hypotézy")` může být chybné a závažnost obou chyb závisí na jejich důsledcích.
- Síla testu popisuje schopnost odhalit předem zvolenou skutečnou velikost efektu; nízká síla zvyšuje riziko, že reálný vztah přehlédneme.
- Dobrá vědecká odpověď spojuje velikost vztahu, `r render_glossary_term("nejistota", display = "nejistotu")`, srovnání s přesnou nulou, `r render_glossary_term("diagnostika-modelu", display = "diagnostiku modelu")` a omezení studie.

### Závěrečná otázka

Kolega vám ukáže jen titulek „p < 0,05: severní krabi jsou větší“. Jaké tři otázky mu položíte, než titulku uvěříte?

::: {.box-result collapse="true" title="Možné otázky"}
O kolik milimetrů na stupeň se krabi liší a je to biologicky hodně? Jak široký je `r render_glossary_term("interval-spolehlivosti", display = "interval spolehlivosti")`? Kolik slanisek studie měla a jak byla vybrána? Sedí přímka k datům a jsou slaniska nezávislá? A tvrdí titulek jen souvislost, nebo i příčinu?
:::

## Slovníček pojmů {.unnumbered}

Pojmy použité v těchto materiálech a jejich definice. Úplný průběžně aktualizovaný slovníček je dostupný také [online](https://cuni-natur-biostatistics.github.io/slovnik/).

```{r}
#| label: slovnicek
#| echo: false
#| results: asis
if (knitr::is_html_output()) {
  glossary::glossary_table()
} else {
  cat(
    paste0(
      "Interaktivní slovníček je dostupný online: ",
      "https://cuni-natur-biostatistics.github.io/slovnik/"
    )
  )
}
```