---
title: "Vztahy dvou proměnných a vizualizace"
subtitle: "Podkladové čtenářské materiály k přednášce"
lang: cs
filters:
  - ../theme/semantic-boxes.lua
  - glossary-headings.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/02/26
date-format: long
date-modified: last-modified
---

```{r}
#| label: nastaveni
#| include: false
if (!require("here")) install.packages("here")

here::i_am("Learning_materials/skripta.qmd")
knitr::opts_knit$set(root.dir = here::here())

if (!require("tidyverse")) install.packages("tidyverse")
if (!require("ggbeeswarm")) install.packages("ggbeeswarm")
if (!require("palmerpenguins")) install.packages("palmerpenguins")
if (!require("glossary")) install.packages("glossary")
if (!require("cowplot")) install.packages("cowplot")
if (!require("kableExtra")) install.packages("kableExtra")

library(tidyverse)
library(ggbeeswarm)
library(palmerpenguins)
library(glossary)
glossary::glossary_popup("hover")
library(cowplot)
library(kableExtra)

# Load the shared course figure theme before lesson-specific refinements.
source(here::here("R", "set_r_theme.R"))

# Glossary: point to the canonical slovnik pojmy.yaml hosted on GitHub.
# Terms are looked up by slug ID (e.g. "histogram", "median", "boxplot").
# Usage in text: `r render_glossary_term("slug", display = "Czech term")`
# The glossary_style() call applies brand colours to the hover tooltip (HTML only).
local_pojmy <-
  tempfile(fileext = ".yaml")
tryCatch(
  {
    download.file(
      url = "https://raw.githubusercontent.com/CUNI-NATUR-Biostatistics/slovnik/main/pojmy.yaml",
      destfile = local_pojmy,
      quiet = TRUE
    )
    glossary::glossary_path(local_pojmy)
  },
  error = function(e) {
    message("Glossary YAML unavailable: ", conditionMessage(e))
  }
)
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"))
  }
)
source(here::here("R", "Functions", "make_korelace_panel.R"))

```

```{r}
#| label: nastaveni-slovnik
#| echo: false
#| warning: false
#| message: false
#| results: asis
if (knitr::is_html_output()) {
  glossary::glossary_style(
    color = biostat_cols["indigo_velvet"],
    text_decoration = "underline",
    def_bg = biostat_cols["white"],
    def_color = biostat_cols["graphite"]
  )

}
```


```{r}
#| label: nastaveni-vizualizace
#| echo: false
#| message: false
#| warning: false
ggplot2::theme_set(
  theme_biostat() +
    ggplot2::theme(
      panel.grid.minor = ggplot2::element_blank(),
      strip.background = ggplot2::element_rect(
        fill  = biostat_cols["light_gray"],
        color = biostat_cols["grey_olive"]
      ),
      strip.text = ggplot2::element_text(
        face   = "bold",
        colour = biostat_cols["indigo_velvet"]
      )
    )
)

col_data        <- biostat_cols["grey_olive"]
col_highlight   <- biostat_cols["orange"]
col_median_line <- biostat_cols["indigo_velvet"]
# Species colours follow the palmerpenguins dataset's own established palette
# (https://allisonhorst.github.io/palmerpenguins/), kept by decision of
# O. Mottl (2026-10-06) as a documented exception to the _brand categorical
# palette. Species are groups here, not semantic data/model/focus roles.
barvy_druhu <-
  c(Adelie = "darkorange", Chinstrap = "mediumpurple", Gentoo = "cyan4")

```

## Úvod

V minulé lekci jsme se naučili určit, co představuje jeden řádek, rozpoznat typ jedné `r render_glossary_term("promenna", display = "proměnné")` a zvolit pro ni vhodný souhrn a graf. V této lekci položíme druhou proměnnou vedle první. Stejná pravidla stále platí, ale otázka se mění: nehledáme už jen tvar jednoho `r render_glossary_term("rozdeleni", display = "rozdělení")`, nýbrž vztah mezi dvěma znaky.

U živočichů často čekáme, že větší tělo bude souviset s většími orgány nebo končetinami. U tučňáků si proto můžeme položit jednoduchou otázku: mají jedinci s delší ploutví také větší tělesnou hmotnost? Tato otázka je vhodná pro první krok do světa vztahů dvou proměnných, protože nás nutí přemýšlet nejen o tvaru vzoru, ale i o `r render_glossary_term("variabilita", display = "variabilitě")` dat.

::: {.box-question title="Otázka před prvním grafem"}

Co byste čekali ještě před tím, než se podíváme na data?

- **A:** Delší ploutev téměř vždy znamená vyšší hmotnost.
- **B:** Delší ploutev bude s vyšší hmotností souviset, ale ne dokonale.
- **C:** Délka ploutve a hmotnost spolu souviset nebudou.

Než budete číst dál, zkuste si vybrat, která odpověď vám připadá nejpravděpodobnější, a proč.
:::

::: {.box-result collapse="true" title="Odpověď: očekáváme nedokonalý vztah"}

Biologicky nejpravděpodobnější je **možnost B**: délka ploutve bude s tělesnou hmotností souviset, ale ne dokonale. Tělesné rozměry tučňáků jsou vzájemně provázány, ale každý jedinec je jiný — liší se druhem, pohlavím, věkem i podmínkami prostředí. Dokonalá shoda by v reálných biologických datech byla překvapivá.
:::

## Výsledky učení

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

- přečíst `r render_glossary_term("bodovy-graf", display = "bodový graf")` dvou `r render_glossary_term("numericka-promenna", display = "numerických proměnných")`: směr, tvar a rozptýlení bodů
- vysvětlit, proč jsou biologická data rozptýlená i při zřetelném vztahu a co k popisu přidá `r render_glossary_term("korelace", display = "korelace")`
- zvolit vhodný první graf podle typu dvojice `r render_glossary_term("promenna", display = "proměnných")`
- rozlišit popis asociace od tvrzení o její příčině

## Jaká data budeme používat

Dataset Palmer Penguins zachycuje měření dospělých tučňáků z Palmerova souostroví v letech 2007–2009. Jeden řádek odpovídá jednomu změřenému jedinci. Délka ploutve je v milimetrech a tělesná hmotnost v gramech; druh a ostrov označují skupiny. Původní data zveřejnil projekt Palmer Penguins pod CC0.

Soubor [palmer_penguins.csv](../data/palmer_penguins.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čtením vznikne tabulka `data_puvodni`; prázdná pole přečteme jako chybějící hodnoty `NA`. Příkaz `dim()` ukáže počet řádků a sloupců, `head()` první tři řádky vybraných sloupců:


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

```{r}
#| label: nacist-a-prohlednout-puvodni-data
#| echo: true
#| eval: true
data_puvodni <-
  read.csv(
    file = "data/palmer_penguins.csv",
    na.strings = ""
  )
dim(x = data_puvodni)
head(
  x = data_puvodni[c("species", "flipper_length_mm", "body_mass_g")],
  n = 3
)
```

```{=typst}
]
```


Sloupec `species` obsahuje druh, `flipper_length_mm` délku ploutve a `body_mass_g` hmotnost. Zkontrolujeme, zda některé z obou měření chybí. Funkce `is.na()` označí chybějící hodnoty a `sum()` jejich počet:

```{r}
#| label: zkontrolovat-chybejici-parovana-mereni
#| echo: true
#| eval: true
sum(is.na(x = data_puvodni$flipper_length_mm))
sum(is.na(x = data_puvodni$body_mass_g))
```

Bod potřebuje známou délku i známou hmotnost **téhož** jedince. `complete.cases()` označí řádky s oběma hodnotami; hranaté závorky zachovají celé tyto řádky, aby se dvojice nerozpojily. V nové tabulce použijeme české názvy sloupců:


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

```{r}
#| label: pripravit-data-tucnaci
#| echo: true
#| eval: true
vec_uplne <-
  complete.cases(
    data_puvodni$flipper_length_mm,
    data_puvodni$body_mass_g
  )
data_pary <- data_puvodni[vec_uplne, ]
data_tucnaci <-
  data.frame(
    druh_tucnaka = data_pary$species,
    ostrov = data_pary$island,
    delka_ploutve_mm = data_pary$flipper_length_mm,
    hmotnost_g = data_pary$body_mass_g
  )
```

```{=typst}
]
```


Původní soubor má `r nrow(data_puvodni)` řádků. `r nrow(data_puvodni) - nrow(data_tucnaci)` řádky nemají obě měření; po jejich vynechání zůstává `r nrow(data_tucnaci)` úplných dvojic, tedy `r render_glossary_term("pozorovani", display = "pozorování")`, u nichž známe délku ploutve i tělesnou hmotnost. V datech se objevují `r length(unique(data_tucnaci$druh_tucnaka))` druhy (`r paste(sort(unique(data_tucnaci$druh_tucnaka)), collapse = ", ")`) ze `r length(unique(data_tucnaci$ostrov))` ostrovů. Známé délky ploutve se pohybují od `r min(data_tucnaci$delka_ploutve_mm)` do `r max(data_tucnaci$delka_ploutve_mm)` mm a hmotnosti od `r min(data_tucnaci$hmotnost_g)` do `r max(data_tucnaci$hmotnost_g)` g. Kontrola rozsahu pomáhá odhalit zjevnou chybu jednotek nebo nesmyslnou hodnotu ještě před kreslením grafu.

::: {.box-data collapse="true" title="Doplňující: jak vznikl soubor palmer_penguins.csv"}

Měření nasbírala Kristen Gorman s výzkumnou stanicí Palmer Station (program Long Term Ecological Research) na třech ostrovech Palmerova souostroví v Antarktidě v letech 2007–2009 ([Gorman, Williams a Fraser 2014](https://doi.org/10.1371/journal.pone.0090081)). Allison Horst, Alison Hill a Kristen Gorman data zpřístupnily v balíčku R palmerpenguins ([Horst a kol. 2020](https://doi.org/10.5281/zenodo.3960218)) pod licencí CC0, takže je smíme volně používat a sdílet.

Soubor pro tuto lekci je přímým exportem tabulky `penguins` z balíčku palmerpenguins verze 0.1.1. Obsahuje všech `r nrow(data_puvodni)` řádků a `r ncol(data_puvodni)` původních sloupců (`r paste0("\x60", names(data_puvodni), "\x60", collapse = ", ")`). Žádný řádek jsme nevynechali, nic nepřepočítali ani nepřejmenovali. Chybějící hodnota je v souboru prázdné pole, a proto ji při načtení převádíme na `NA` pomocí `na.strings = ""`. Chybí `r sum(is.na(data_puvodni$flipper_length_mm))` délky ploutve, `r sum(is.na(data_puvodni$body_mass_g))` hmotnosti a `r sum(is.na(data_puvodni$sex))` údajů o pohlaví.

Soubor vytvořil skript [`R/prepare_penguin_data.R`](https://github.com/CUNI-NATUR-Biostatistics/L02/blob/main/R/prepare_penguin_data.R), který před uložením ověří verzi balíčku, počet řádků a sloupců i počty chybějících hodnot. Stejnou tabulku si můžete z balíčku uložit i v základním R; potřebujete k tomu nainstalovaný balíček palmerpenguins:

```{r}
#| label: ulozit-tabulku-tucnaku-do-csv
#| echo: true
#| eval: false
data_export <-
  as.data.frame(palmerpenguins::penguins)
write.csv(
  x = data_export,
  file = "data/palmer_penguins.csv",
  row.names = FALSE,
  na = ""
)
```

Pro výuku používáme samostatný soubor CSV místo balíčku proto, že se nezmění s novou verzí balíčku a v praktiku stačí načíst jediný soubor pomocí `read.csv()`.
:::

Na tomto datasetu je užitečné to, že biologický příběh dává smysl i bez složité statistiky. Delší ploutev může souviset s větší tělesnou velikostí, ale zároveň víme, že do výsledku vstupují i další zdroje rozdílů mezi jedinci.

![Druhy tučňáků v datasetu Palmer Penguins. Ilustrace: Allison Horst.](images/lter_penguins.png){width="85%"}

Na obrázku výše jsou `r length(unique(data_tucnaci$druh_tucnaka))` druhy, se kterými budeme pracovat. Už z ilustrace je vidět, že se mezi nimi mohou lišit tělesné proporce. To je dobrá připomínka, že část rozptýlení bodů v `r render_glossary_term("bodovy-graf", display = "bodovém grafu")` může souviset i s tím, že nesledujeme jediný homogenní typ jedince.

Jako první krok se ale na data podíváme co nejjednodušším způsobem: bez rozlišování druhů, jen jako na vztah mezi délkou ploutve a tělesnou hmotností.


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

::: {.box-data collapse="true" title="Doplňující: sto bodů nemusí znamenat sto nezávislých informací"}

V této tabulce představuje každý řádek jiného změřeného tučňáka. V jiných studiích však může několik řádků patřit stejnému jedinci, mláďatům ze stejného hnízda, organismům ze stejné lokality nebo měřením ze stejného dne. Takové body mohou být navzájem podobnější než body z různých jedinců či lokalit.

Počet teček v grafu proto nemusí být totožný s počtem nezávislých informačních jednotek. Graf vztah stále užitečně ukáže, ale `r render_glossary_term("nejistota", display = "nejistotu")` nelze bezpečně počítat, jako by každý řádek přinášel zcela novou informaci. `r render_glossary_term("hierarchicka-data", display = "Hierarchická data")`, `r render_glossary_term("opakovana-mereni", display = "opakovaná měření")` a pseudoreplikaci budeme podrobně řešit v pozdější lekci. <!-- TODO: add slug pseudoreplikace to pojmy.yaml -->
:::

```{=typst}
]
```



## První graf: délka ploutve a tělesná hmotnost

Nejprve si zobrazíme co nejjednodušší `r render_glossary_term("bodovy-graf", display = "bodový graf")`. V tuto chvíli nechceme nic testovat ani odhadovat pomocí modelu. Chceme si všimnout, jak data vypadají.

```{r}
#| label: prvni-bodovy-graf
#| echo: false
#| message: false
#| warning: false
#| fig-cap: "Bodový graf délky ploutve (mm) a tělesné hmotnosti (g) v datasetu Palmer Penguins."
ggplot2::ggplot(
  data    = data_tucnaci,
  mapping = ggplot2::aes(x = delka_ploutve_mm, y = hmotnost_g)
) +
  ggplot2::geom_point(colour = col_data, alpha = 0.7, size = 2) +
  ggplot2::labs(
    x = "Délka ploutve (mm)",
    y = "Tělesná hmotnost (g)"
  )
```

Tento první graf obsahuje `r nrow(data_tucnaci)` tučňáků. Stejná data si později vyzkoušíte číst i pomocí modelu.

Než si podobný graf nakreslíme sami, je užitečné v R rychle zkontrolovat prvních několik hodnot na obou osách:

```{r}
#| label: nahled-na-promenne-pro-prvni-graf
#| echo: true
#| eval: true
head(
  # Prvních šest délek ploutve — hodnoty na ose x
  x = data_tucnaci$delka_ploutve_mm
)

head(
  # Prvních šest hmotností — hodnoty na ose y
  x = data_tucnaci$hmotnost_g
)
```

Jakmile víme, které dvě `r render_glossary_term("promenna", display = "proměnné")` chceme porovnat, můžeme tentýž první bodový graf rychle nakreslit v základním R:

```{r}
#| label: prvni-bodovy-graf-v-zakladnim-r
#| echo: true
#| eval: true
plot(
  # Každý řádek dat vytvoří jeden bod se dvěma souřadnicemi
  x = data_tucnaci$delka_ploutve_mm,
  y = data_tucnaci$hmotnost_g,
  # Popisky připomínají význam a jednotky obou os
  xlab = "Délka ploutve (mm)",
  ylab = "Tělesná hmotnost (g)",
  # Plné, průhledné body omezí překrývání
  pch = 16,
  col = rgb(red = 0.54, green = 0.54, blue = 0.54, alpha = 0.6)
)
```

## Co v grafu vidíme

Než budete číst dál, zkuste se u grafu zastavit u tří krátkých otázek:

- Jakým směrem se oblak bodů zvedá nebo klesá?
- Jak široké je rozptýlení bodů kolem hlavního směru?
- Vidíte body, které se od většiny ostatních odlišují?

::: {.box-result collapse="true" title="Odpověď: vztah je kladný, ale ne dokonalý"}

Celkový směr vztahu je zřetelně rostoucí: tučňáci s delší ploutví bývají těžší. Body ale neleží na jedné přímce; i při podobné délce ploutve se může tělesná hmotnost lišit o stovky gramů. Vztah v datech tedy opravdu vidíme, ale není natolik jednoduchý, aby se dal popsat jedinou „čistou“ čarou bez dalších otázek.
:::

## Jaké `r render_glossary_term("promenna", display = "proměnné")` budeme sledovat

Po prvním pohledu na graf se teď vyplatí pojmenovat, co jsme vlastně zobrazili. Pracujeme se dvěma `r render_glossary_term("numericka-promenna", display = "numerickými proměnnými")`:

- **osa x**: délka ploutve (anglicky *flipper length*) v milimetrech (`flipper_length_mm`)
- **osa y**: tělesná hmotnost (anglicky *body mass*) v gramech (`body_mass_g`)
- doplňkově `r render_glossary_term("kategoricka-promenna", display = "kategoriální proměnná")` druh tučňáka (anglicky *species*; `species`)

Když jen popisujeme, zda dvě měření souvisejí, je přiřazení os libovolné. Kdybychom osy prohodili — délku svisle a hmotnost vodorovně — bodový mrak by vypadal jinak, ale směr i těsnost vztahu by zůstaly naprosto stejné. Pojmy *`r render_glossary_term("prediktor", display = "prediktor")`* a *`r render_glossary_term("odezva", display = "odezva")`*, které popisují asymetrii os (jedna proměnná „vysvětluje“ druhou), patří k regresním modelům; ty si ukážeme v pozdějších lekcích.

```{r}
#| label: flip-axes-demo
#| echo: false
#| message: false
#| warning: false
#| fig-cap: "Stejná data, dvě orientace os. Tvar mraku bodů se změní, ale směr i těsnost vztahu zůstanou stejné."
p_flip_standard <-
  ggplot2::ggplot(
    data    = data_tucnaci,
    mapping = ggplot2::aes(x = delka_ploutve_mm, y = hmotnost_g)
  ) +
  ggplot2::geom_point(colour = col_data, alpha = 0.7, size = 2) +
  ggplot2::labs(
    x        = "Délka ploutve (mm)",
    y        = "Tělesná hmotnost (g)",
    subtitle = "Délka ploutve na ose x"
  )

p_flip_flipped <-
  ggplot2::ggplot(
    data    = data_tucnaci,
    mapping = ggplot2::aes(x = hmotnost_g, y = delka_ploutve_mm)
  ) +
  ggplot2::geom_point(colour = col_data, alpha = 0.7, size = 2) +
  ggplot2::labs(
    x        = "Tělesná hmotnost (g)",
    y        = "Délka ploutve (mm)",
    subtitle = "Hmotnost na ose x"
  )

cowplot::plot_grid(p_flip_standard, p_flip_flipped, ncol = 2)
```

`r render_glossary_term("bodovy-graf", display = "Bodový graf")` nám umožňuje podívat se současně na směr vztahu i na to, jak moc jsou hodnoty rozptýlené kolem hlavního vzoru.

Když ale sledujeme vztah dvou proměnných obecněji, typ grafu vybíráme podle toho, jaké proměnné máme k dispozici:

- **Dvě numerické proměnné:** hledáme směr vztahu a rozptýlení bodů. Dobrou první volbou je bodový graf (anglicky *scatter plot*).
- **Kategoriální a numerická proměnná:** hledáme rozdíly mezi skupinami a překryv `r render_glossary_term("rozdeleni", display = "rozdělení")`. Dobrou první volbou jsou body rozdělené do skupin, případně `r render_glossary_term("boxplot", display = "krabicový graf")` (anglicky *box plot*) s body.
- **Dvě kategoriální proměnné:** sledujeme, kolik `r render_glossary_term("pozorovani", display = "pozorování")` spadá do jednotlivých kombinací. Začít můžeme tabulkou `r render_glossary_term("cetnost", display = "četností")` nebo `r render_glossary_term("sloupcovy-graf", display = "sloupcovým grafem")` (anglicky *bar chart*).

Obecné pravidlo tedy nezní jen „je v datech nějaký vztah?“, ale také „jaký graf vůbec dává pro tuto dvojici proměnných smysl?“ Ke všem třem situacím se vrátíme v závěru lekce. Do té doby zůstaneme u první z nich a budeme dál číst jeden bodový graf dvou numerických proměnných.

## Rozptýlení ve dvou směrech

U jednorozměrných dat sledujeme rozptýlení jen v jednom směru, například jak moc se hodnoty hmotnosti liší kolem `r render_glossary_term("prumer", display = "průměru")`. V `r render_glossary_term("bodovy-graf", display = "bodovém grafu")` ale kolísají obě osy současně. Jeden tučňák může mít o něco delší ploutev než jiný a zároveň se může lišit i hmotností. Výsledkem je mrak bodů, který má šířku i výšku.

Na obrázku níže je stejný vztah zobrazený dvakrát. Vlevo si všímáme, že i při podobné délce ploutve mohou mít jedinci docela odlišnou hmotnost. Vpravo naopak vidíme, že i při podobné hmotnosti se může měnit délka ploutve.

```{r}
#| label: rozptyl-ve-dvou-smerech-graf
#| echo: false
#| message: false
#| warning: false
#| fig-cap: "Stejný bodový graf můžeme číst ve dvou směrech: při podobné délce ploutve se může lišit hmotnost a při podobné hmotnosti se může lišit délka ploutve. Zvýrazněné body leží uvnitř příslušného pásu."
# Same teaching strips as in the presentation.
strip_delka_lwr   <- 188
strip_delka_upr   <- 192
strip_hmotnost_lwr <- 3800L
strip_hmotnost_upr <- 4200L

data_rozptyl <-
  data_tucnaci |>
  dplyr::mutate(
    in_strip_delka    = delka_ploutve_mm >= strip_delka_lwr &
      delka_ploutve_mm <= strip_delka_upr,
    in_strip_hmotnost = hmotnost_g >= strip_hmotnost_lwr &
      hmotnost_g <= strip_hmotnost_upr
  )

p_rozptyl_delka <-
  ggplot2::ggplot(
    data    = data_rozptyl,
    mapping = ggplot2::aes(x = delka_ploutve_mm, y = hmotnost_g)
  ) +
  ggplot2::annotate(
    "rect",
    xmin  = strip_delka_lwr,
    xmax  = strip_delka_upr,
    ymin  = -Inf,
    ymax  = Inf,
    fill  = col_highlight,
    alpha = 0.2
  ) +
  ggplot2::geom_point(
    data   = dplyr::filter(data_rozptyl, !in_strip_delka),
    colour = col_data,
    alpha  = 0.5,
    size   = 1.6
  ) +
  ggplot2::geom_point(
    data   = dplyr::filter(data_rozptyl, in_strip_delka),
    colour = col_highlight,
    alpha  = 0.9,
    size   = 2.5
  ) +
  ggplot2::labs(
    x        = "Délka ploutve (mm)",
    y        = "Tělesná hmotnost (g)",
    subtitle = "Podobná délka ploutve"
  )

p_rozptyl_hmotnost <-
  ggplot2::ggplot(
    data    = data_rozptyl,
    mapping = ggplot2::aes(x = delka_ploutve_mm, y = hmotnost_g)
  ) +
  ggplot2::annotate(
    "rect",
    xmin  = -Inf,
    xmax  = Inf,
    ymin  = strip_hmotnost_lwr,
    ymax  = strip_hmotnost_upr,
    fill  = col_highlight,
    alpha = 0.2
  ) +
  ggplot2::geom_point(
    data   = dplyr::filter(data_rozptyl, !in_strip_hmotnost),
    colour = col_data,
    alpha  = 0.5,
    size   = 1.6
  ) +
  ggplot2::geom_point(
    data   = dplyr::filter(data_rozptyl, in_strip_hmotnost),
    colour = col_highlight,
    alpha  = 0.9,
    size   = 2.5
  ) +
  ggplot2::labs(
    x        = "Délka ploutve (mm)",
    y        = "Tělesná hmotnost (g)",
    subtitle = "Podobná hmotnost"
  )

cowplot::plot_grid(p_rozptyl_delka, p_rozptyl_hmotnost, ncol = 2)
```

To je důležité pro interpretaci. Když vidíme podlouhlý oblak bodů směřující zleva dole doprava nahoru, znamená to, že vedle celkového trendu stále existuje velké rozptýlení bodů. Dva jedinci s téměř stejnou délkou ploutve nemusí mít stejnou hmotnost a dva podobně těžcí jedinci nemusí mít stejně dlouhou ploutev.

### Připomenutí: `r render_glossary_term("rozptyl", display = "rozptyl")` a `r render_glossary_term("smerodatna-odchylka", display = "směrodatná odchylka")` jedné proměnné

```{r}
#| label: maly-vzorek-rovnice
#| echo: false
#| message: false
#| warning: false
# Malý ilustrační vzorek 4 tučňáků druhu Adelie vybraný tak,
# aby byl v datech dobře vidět kladný vztah mezi délkou ploutve a hmotností.
idx_maly <-
  which(data_tucnaci$druh_tucnaka == "Adelie")[c(1, 2, 9, 14)]

maly_vzorek_delka <-
  data_tucnaci$delka_ploutve_mm[idx_maly]

n_maly <-
  length(maly_vzorek_delka)

prumer_maly_delka <-
  mean(maly_vzorek_delka)

odchylky_delka <- maly_vzorek_delka - prumer_maly_delka

ctverce_odchylek <- odchylky_delka^2
soucet_ctvercu <-
  sum(ctverce_odchylek)
rozptyl_maly <- soucet_ctvercu / (n_maly - 1L)

prumer_maly_delka_txt <-
  formatC(prumer_maly_delka, format = "f", digits = 2, decimal.mark = ",")

ctverce_txt <-
  formatC(ctverce_odchylek, format = "f", digits = 4, decimal.mark = ",")
odchylky_delka_txt <-
  formatC(odchylky_delka, format = "f", digits = 2, decimal.mark = ",")

soucet_ctvercu_txt <-
  formatC(soucet_ctvercu, format = "f", digits = 2, decimal.mark = ",")
rozptyl_maly_txt <-
  formatC(rozptyl_maly, format = "f", digits = 1, decimal.mark = ",")
```

Z předchozí lekce připomeňme souhrny rozptýlení jedné `r render_glossary_term("promenna", display = "proměnné")`: **výběrový rozptyl** (anglicky *sample variance*; zkráceně také jen **rozptyl**) a z něj odvozenou výběrovou směrodatnou odchylku (anglicky *sample standard deviation*, zkratka **SD**).

Výpočet výběrového rozptylu slovy: **od každé hodnoty odečteme průměr, odchylky umocníme na druhou, sečteme je a součet vydělíme počtem `r render_glossary_term("pozorovani", display = "pozorování")` minus jedna**. Jako příklad použijeme malý vzorek `r n_maly` tučňáků druhu Adelie (průměr délky ploutve: `r prumer_maly_delka_txt` mm). Tabulka ukazuje délku ploutve každého z nich, odchylku od průměru a čtverec té odchylky:

```{r}
#| label: tabulka-rozptyl
#| echo: false
#| warning: false
#| message: false
knitr::kable(
  data.frame(
    tucnak   = seq_len(n_maly),
    delka    = maly_vzorek_delka,
    odchylka = odchylky_delka,
    ctverec  = ctverce_odchylek
  ),
  col.names = c(
    "Tučňák",
    "Délka ploutve (mm)",
    paste0("Odchylka od průměru (", prumer_maly_delka_txt, " mm)"),
    "Čtverec odchylky (mm²)"
  ),
  align = "c"
)
```

Stejný výpočet zapsaný názvy veličin:

$$
\begin{gathered}
\text{výběrový rozptyl délky ploutve} \\
= \frac{
  \substack{\text{součet čtverců odchylek} \\ \text{délek ploutví od průměru}}
}{\text{počet tučňáků} - \text{jedna}}
\end{gathered}
$$

Po dosazení délek ploutví a jejich průměru z tabulky:

$$\begin{align}
\text{rozptyl}
  &= \frac{\substack{(`r maly_vzorek_delka[1]` - `r prumer_maly_delka_txt`)^2 + (`r maly_vzorek_delka[2]` - `r prumer_maly_delka_txt`)^2 \\ + (`r maly_vzorek_delka[3]` - `r prumer_maly_delka_txt`)^2 + (`r maly_vzorek_delka[4]` - `r prumer_maly_delka_txt`)^2}}{`r n_maly` - 1} \\[6pt]
  &= \frac{`r ctverce_txt[1]` + `r ctverce_txt[2]` + `r ctverce_txt[3]` + `r ctverce_txt[4]`}{`r n_maly - 1`} \\[6pt]
  &= \frac{`r soucet_ctvercu_txt`}{`r n_maly - 1`} \approx `r rozptyl_maly_txt`\text{ mm}^2
\end{align}$$

Směrodatná odchylka vrací rozptýlení do původních jednotek: je odmocninou výběrového rozptylu. Pro stejnou tabulku tedy platí:

$$\text{směrodatná odchylka} = \sqrt{\text{výběrový rozptyl}}$$

$$\sqrt{\frac{`r soucet_ctvercu_txt`}{`r n_maly - 1`}} \approx `r formatC(sqrt(rozptyl_maly), digits = 2, format = "f", decimal.mark = ",")`\text{ mm}$$

V obecném zápisu označíme délku ploutve $x$, její výběrový rozptyl $s_x^2$ a výběrovou směrodatnou odchylku $s_x$. Pro libovolný soubor $n$ pozorování platí:

$$s^2_x = \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})^2, \qquad s_x = \sqrt{s^2_x}$$

Ve vzorci $x_i$ označuje délku ploutve jednoho tučňáka, $\bar{x}$ průměrnou délku ploutve, $n$ počet pozorování, $s_x^2$ výběrový rozptyl a $s_x$ směrodatnou odchylku v původních jednotkách.

V celém souboru tučňáků vychází směrodatná odchylka délky ploutve přibližně `r round(sd(data_tucnaci$delka_ploutve_mm), 1)` mm a tělesné hmotnosti přibližně `r round(sd(data_tucnaci$hmotnost_g))` g. Tato čísla říkají, o kolik milimetrů nebo gramů se typický jedinec odchyluje od průměru celého souboru. Stejné odchylky od průměru použijeme později u kovariance, kde budeme sledovat dvě proměnné zároveň.

## Proč jsou data rozptýlená

`r render_glossary_term("variabilita", display = "Variabilita")` v biologických datech není chyba, kterou bychom měli odstranit. Je to vlastnost reálného světa. Dva tučňáci mohou mít podobně dlouhou ploutev, ale lišit se druhem, pohlavím, věkem nebo celkovou tělesnou stavbou. K rozptýlení bodů přispívá i běžná nepřesnost měření a to, že žádný biologický vztah nebývá dokonale přesný.

Na tomto místě je užitečné rozlišit několik zdrojů variability, které se v grafu míchají dohromady:

- **biologické rozdíly mezi jedinci**: dva podobní jedinci nemusí být stejně těžcí ani stejně stavění
- **rozdíly mezi skupinami**: část variability vzniká tím, že v jednom grafu spojujeme více druhů
- **nepřesnost a omezení měření**: reálná data nikdy nejsou úplně bez chyby
- **další nevysvětlené vlivy**: sledujeme jen malou část toho, co hmotnost organismu ovlivňuje

V našem datasetu se délka ploutve pohybuje přibližně od `r min(data_tucnaci$delka_ploutve_mm)` do `r max(data_tucnaci$delka_ploutve_mm)` mm a tělesná hmotnost od `r min(data_tucnaci$hmotnost_g)` do `r max(data_tucnaci$hmotnost_g)` g. Už z tohoto rozsahu je vidět, že nejde o několik téměř stejných jedinců, ale o biologicky pestrý soubor měření.


## Stejný graf, ale s viditelnými druhy

Jeden ze zdrojů `r render_glossary_term("variabilita", display = "variability")` můžeme v datech přímo ukázat. Vrátíme se ke stejnému `r render_glossary_term("bodovy-graf", display = "bodovému grafu")` a odkryjeme ještě jednu informaci: druh tučňáka. Původní šedý oblak bodů se rozpadne na několik podskupin, které se částečně překrývají, ale zároveň mají vlastní polohu a vlastní rozptýlení.

Biologicky je to důležité: v původním bodovém grafu jsme viděli „jeden“ vztah, ale část tohoto vztahu je tvořena tím, že různé druhy mají různé typické kombinace délky ploutve a hmotnosti. Jinými slovy: část rozptýlení není jen náhodný šum, ale odráží i skutečnou strukturu dat.

```{r}
#| label: bodovy-graf-podle-druhu
#| echo: false
#| message: false
#| warning: false
#| fig-cap: "Stejný bodový graf délky ploutve a tělesné hmotnosti, body obarvené podle druhu tučňáka."
ggplot2::ggplot(
  data    = data_tucnaci,
  mapping = ggplot2::aes(
    x      = delka_ploutve_mm,
    y      = hmotnost_g,
    colour = druh_tucnaka
  )
) +
  ggplot2::geom_point(alpha = 0.7, size = 2) +
  ggplot2::scale_colour_manual(
    values = barvy_druhu,
    name   = "Druh"
  ) +
  ggplot2::labs(
    x = "Délka ploutve (mm)",
    y = "Tělesná hmotnost (g)"
  )
```

Jakmile jsou druhy vidět, je snazší pochopit, proč byl původní oblak bodů tak široký. Část variability vzniká uvnitř druhů, ale část je dána i tím, že v jednom obrázku mícháme více biologicky odlišných skupin. To je přesně okamžik, kdy nám doplňková `r render_glossary_term("kategoricka-promenna", display = "kategoriální proměnná")` pomůže lépe číst vztah dvou `r render_glossary_term("numericka-promenna", display = "numerických proměnných")`.

::: {.box-data collapse="true" title="Doplňující: souhrnný vztah se může lišit od vztahů uvnitř skupin"}

Šedý graf odpovídal na otázku, jak spolu délka ploutve a hmotnost souvisejí **napříč všemi tučňáky dohromady**. Barevný graf navíc dovoluje položit jinou otázku: jak spolu souvisejí **uvnitř každého druhu**. Tyto odpovědi nemusejí být stejné, protože celkový vzor kombinuje rozdíly mezi jedinci uvnitř druhů s rozdíly mezi typickými hodnotami jednotlivých druhů.

V pozdější lekci se k tomuto mísení vrátíme při tématu `r render_glossary_term("konfundovani", display = "konfundování")`. Zde stačí bezpečné pravidlo: kdykoli barvy odhalí skupiny v jiných částech grafu, popište zvlášť celkový vzor a vzory uvnitř skupin. Zatím nepotřebujeme model ani nový výpočet; důležité je rozlišit dvě otázky, které lze položit nad stejnými body.
:::


## Jak se dvě `r render_glossary_term("promenna", display = "proměnné")` pohybují společně?

```{r}
#| label: pripravit-kovarianci-celeho-souboru
#| echo: false
#| message: false
#| warning: false
# Covariance and normalization
kovariance_val <-
  cov(
    data_tucnaci$delka_ploutve_mm,
    data_tucnaci$hmotnost_g
  )

kovariance_txt <-
  formatC(
    kovariance_val,
    format = "f",
    digits = 0,
    decimal.mark = ","
  )
```

`r render_glossary_term("rozptyl", display = "Rozptyl")` jedné proměnné popisuje, jak moc jsou její hodnoty rozptýleny kolem `r render_glossary_term("prumer", display = "průměru")`. U dvou proměnných nás navíc zajímá, zda se jejich odchylky od průměrů pohybují *společně*: mají tučňáci s nadprůměrně dlouhou ploutví také nadprůměrnou hmotnost?

Odpověď nejlépe uvidíme v grafu. Do `r render_glossary_term("bodovy-graf", display = "bodového grafu")` přidáme dvě přerušované čáry — svislou v průměrné délce ploutve (`r format(round(mean(data_tucnaci$delka_ploutve_mm), 1), decimal.mark = ",")` mm) a vodorovnou v průměrné hmotnosti (`r round(mean(data_tucnaci$hmotnost_g))` g). Graf se tím rozdělí na čtyři kvadranty:

```{r}
#| label: kovariance-schema
#| echo: false
#| message: false
#| warning: false
#| fig-cap: "Bodový graf tučňáků s průměrnou délkou ploutve a průměrnou hmotností (přerušované čáry), které dělí graf na čtyři kvadranty."
prumer_delka    <- mean(data_tucnaci$delka_ploutve_mm)
prumer_hmotnost <-
  mean(data_tucnaci$hmotnost_g)

ggplot2::ggplot(
  data    = data_tucnaci,
  mapping = ggplot2::aes(x = delka_ploutve_mm, y = hmotnost_g)
) +
  ggplot2::geom_vline(
    xintercept = prumer_delka,
    linetype   = "dashed",
    colour     = biostat_cols["graphite"],
    linewidth  = 0.8
  ) +
  ggplot2::geom_hline(
    yintercept = prumer_hmotnost,
    linetype   = "dashed",
    colour     = biostat_cols["graphite"],
    linewidth  = 0.8
  ) +
  ggplot2::geom_point(colour = col_data, alpha = 0.6, size = 1.8) +
  ggplot2::labs(
    x = "Délka ploutve (mm)",
    y = "Tělesná hmotnost (g)"
  )
```

::: {.box-question title="Přečtěte si obrázek"}

Ve kterých kvadrantech se shlukuje nejvíce bodů? Co to říká o tom, jakým směrem se délka ploutve a hmotnost pohybují spolu?
:::

::: {.box-result collapse="true" title="Odpověď: body převažují vpravo nahoře a vlevo dole"}

Nejvíce bodů leží vpravo nahoře a vlevo dole. Tučňáci s nadprůměrně dlouhou ploutví bývají i nadprůměrně těžcí a tučňáci s krátkou ploutví bývají lehčí. Odchylky obou proměnných od průměru tedy mají většinou stejné znaménko.
:::

Pro každého tučňáka můžeme spočítat dvě odchylky — délky ploutve od jejího průměru a hmotnosti od jejího průměru — a vynásobit je. Znaménko tohoto součinu závisí na kvadrantu:

- **vpravo nahoře**: obě odchylky kladné — součin kladný
- **vlevo dole**: obě odchylky záporné — součin je stále kladný, protože $(-)\cdot(-) = (+)$
- **vpravo dole**: odchylky s opačnými znaménky — součin záporný
- **vlevo nahoře**: odchylky s opačnými znaménky — součin záporný

Když součiny všech tučňáků sečteme a vydělíme počtem `r render_glossary_term("pozorovani", display = "pozorování")` minus jedna, dostaneme `r render_glossary_term("kovariance", display = "kovarianci")` (anglicky *covariance*). Rozhoduje znaménko i velikost jednotlivých součinů, ne jen počet bodů v kvadrantech: několik velmi vzdálených bodů může přispět víc než mnoho bodů blízko středu.

```{r}
#| label: pripravit-maly-priklad-kovariance
#| echo: false
#| message: false
#| warning: false
maly_vzorek_hmotnost <-
  data_tucnaci$hmotnost_g[idx_maly]

prumer_maly_hmotnost <-
  mean(maly_vzorek_hmotnost)
odchylky_hmotnost <-
  maly_vzorek_hmotnost - prumer_maly_hmotnost
souciny_odchylek <-
  odchylky_delka * odchylky_hmotnost
soucet_soucinu_maly <-
  sum(souciny_odchylek)
kovariance_maly <-
  soucet_soucinu_maly / (n_maly - 1L)
prumer_maly_hmotnost_txt <-
  formatC(prumer_maly_hmotnost, format = "f", digits = 0, decimal.mark = ",")
odchylky_hmotnost_txt <-
  formatC(odchylky_hmotnost, format = "f", digits = 0, decimal.mark = ",")
souciny_txt_maly <-
  formatC(souciny_odchylek, format = "f", digits = 1, decimal.mark = ",")
soucet_soucinu_maly_txt <-
  formatC(soucet_soucinu_maly, format = "f", digits = 0, decimal.mark = ",")
kovariance_maly_txt <-
  formatC(kovariance_maly, format = "f", digits = 0, decimal.mark = ",")
```

Výpočet kovariance slovy: **u každého tučňáka vynásobíme odchylku délky ploutve od jejího průměru odchylkou hmotnosti od jejího průměru, součiny sečteme a vydělíme počtem pozorování minus jedna**. Pro stejný malý vzorek `r n_maly` tučňáků Adelie jako u rozptylu ukazuje tabulka všechny potřebné hodnoty:

```{r}
#| label: tabulka-kovariance
#| echo: false
#| warning: false
#| message: false
knitr::kable(
  data.frame(
    tucnak     = seq_len(n_maly),
    odch_delka = odchylky_delka,
    odch_hmotn = round(odchylky_hmotnost, 0),
    soucin     = souciny_odchylek
  ),
  col.names = c(
    "Tučňák",
    paste0("Odch. délky (mm), prům. ", prumer_maly_delka_txt, " mm"),
    paste0("Odch. hmotnosti (g), prům. ", prumer_maly_hmotnost_txt, " g"),
    "Součin odchylek (mm·g)"
  ),
  align = "c"
)
```

Stejný výpočet zapsaný názvy obou proměnných:

$$
\begin{gathered}
\text{výběrová kovariance délky ploutve a hmotnosti} \\
= \frac{
  \substack{\text{součet součinů odchylek} \\ \text{délky ploutve a hmotnosti}}
}{\text{počet tučňáků} - \text{jedna}}
\end{gathered}
$$

Po dosazení odchylek z tabulky:

$$\begin{align}
\text{kovariance}
  &= \frac{\substack{(`r odchylky_delka_txt[1]`) \cdot (`r odchylky_hmotnost_txt[1]`) + (`r odchylky_delka_txt[2]`) \cdot (`r odchylky_hmotnost_txt[2]`) \\ + (`r odchylky_delka_txt[3]`) \cdot (`r odchylky_hmotnost_txt[3]`) + (`r odchylky_delka_txt[4]`) \cdot (`r odchylky_hmotnost_txt[4]`)}}{`r n_maly` - 1} \\[6pt]
  &= \frac{`r souciny_txt_maly[1]` + `r souciny_txt_maly[2]` + `r souciny_txt_maly[3]` + `r souciny_txt_maly[4]`}{`r n_maly - 1`} \\[6pt]
  &= \frac{`r soucet_soucinu_maly_txt`}{`r n_maly - 1`} \approx `r kovariance_maly_txt`\text{ mm·g}
\end{align}$$

V obecném zápisu označíme délku ploutve $x$, hmotnost $y$ a jejich výběrovou kovarianci $\operatorname{cov}(x,y)$. Pro $n$ pozorování:

$$\text{cov}(x, y) = \frac{1}{n-1}\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})$$

Zde $x_i$ označuje délku ploutve jednoho tučňáka, $y_i$ jeho hmotnost, $\bar{x}$ a $\bar{y}$ příslušné průměry, $n$ počet tučňáků a $\text{cov}(x,y)$ výběrovou kovarianci obou proměnných.

V celém souboru tučňáků převáží kladné součiny nad zápornými, takže kovariance vychází kladná: přibližně `r kovariance_txt` mm·g.

Kovariance ale závisí na jednotkách. Kdybychom přepočítali délku z milimetrů na centimetry a hmotnost ponechali v gramech, kovariance by se desetkrát zmenšila, přestože se vztah mezi tučňáky nijak nezměnil. Abychom vztah popsali číslem nezávislým na jednotkách, přepočítáme kovarianci pomocí `r render_glossary_term("smerodatna-odchylka", display = "směrodatných odchylek")` obou proměnných.

### Od kovariance k Pearsonově `r render_glossary_term("korelace", display = "korelaci")`

```{r}
#| label: pripravit-normalizaci-kovariance
#| echo: false
#| message: false
#| warning: false
ctverce_odchylek_hmotnost <-
  odchylky_hmotnost^2
soucet_ctvercu_hmotnost <-
  sum(ctverce_odchylek_hmotnost)
pearson_maly <-
  cor(x = maly_vzorek_delka, y = maly_vzorek_hmotnost)
ctverce_hmotnost_txt <-
  formatC(ctverce_odchylek_hmotnost, format = "f", digits = 0, decimal.mark = ",")
soucet_ctvercu_hmotnost_txt <-
  formatC(soucet_ctvercu_hmotnost, format = "f", digits = 0, decimal.mark = ",")
pearson_maly_txt <-
  formatC(pearson_maly, format = "f", digits = 2, decimal.mark = ",")
```

Na stejném malém vzorku `r n_maly` tučňáků si ukážeme plný výpočet Pearsonova korelačního koeficientu. Potřebujeme stejné součiny odchylek jako u kovariance a navíc čtverce odchylek obou proměnných:

```{r}
#| label: tabulka-pearson
#| echo: false
#| warning: false
#| message: false
knitr::kable(
  data.frame(
    tucnak            = seq_len(n_maly),
    odch_delka        = odchylky_delka,
    odch_hmotnost     = round(odchylky_hmotnost, 0),
    soucin_odchylek   = souciny_odchylek,
    ctverec_delky     = ctverce_odchylek,
    ctverec_hmotnosti = round(ctverce_odchylek_hmotnost, 0)
  ),
  col.names = c(
    "Tučňák",
    "Odch. délky (mm)",
    "Odch. hmotnosti (g)",
    "Součin odchylek",
    "Čtverec odch. délky",
    "Čtverec odch. hmotnosti"
  ),
  align = "c"
)
```

Výpočet zapsaný názvy veličin:

$$
\begin{gathered}
\text{Pearsonův korelační koeficient} \\
= \frac{
  \substack{\text{součet součinů odchylek} \\ \text{délky ploutve a hmotnosti}}
}{
  \sqrt{\substack{\text{součet čtverců odchylek} \\ \text{délky ploutve}}}
  \cdot \sqrt{\substack{\text{součet čtverců odchylek} \\ \text{hmotnosti}}}
}
\end{gathered}
$$

Součty čtverců získáme ze dvou posledních sloupců tabulky. Pro délku ploutve rozepíšeme všechny čtyři příspěvky:

$$
\begin{gathered}
(`r odchylky_delka_txt[1]`)^2 + (`r odchylky_delka_txt[2]`)^2 + (`r odchylky_delka_txt[3]`)^2 + (`r odchylky_delka_txt[4]`)^2 \\
= `r ctverce_txt[1]` + `r ctverce_txt[2]` + `r ctverce_txt[3]` + `r ctverce_txt[4]` \\
= `r soucet_ctvercu_txt`\text{ mm}^2
\end{gathered}
$$

Stejně postupujeme u hmotnosti:

$$
\begin{gathered}
(`r odchylky_hmotnost_txt[1]`)^2 + (`r odchylky_hmotnost_txt[2]`)^2 + (`r odchylky_hmotnost_txt[3]`)^2 + (`r odchylky_hmotnost_txt[4]`)^2 \\
= `r ctverce_hmotnost_txt[1]` + `r ctverce_hmotnost_txt[2]` + `r ctverce_hmotnost_txt[3]` + `r ctverce_hmotnost_txt[4]` \\
= `r soucet_ctvercu_hmotnost_txt`\text{ g}^2
\end{gathered}
$$

Součet součinů odchylek vydělíme odmocninou součinu těchto dvou součtů:

$$\begin{align}
\text{korelace}
  &= \frac{\substack{(`r odchylky_delka_txt[1]`) \cdot (`r odchylky_hmotnost_txt[1]`) + (`r odchylky_delka_txt[2]`) \cdot (`r odchylky_hmotnost_txt[2]`) \\ + (`r odchylky_delka_txt[3]`) \cdot (`r odchylky_hmotnost_txt[3]`) + (`r odchylky_delka_txt[4]`) \cdot (`r odchylky_hmotnost_txt[4]`)}}{\sqrt{`r soucet_ctvercu_txt`} \cdot \sqrt{`r soucet_ctvercu_hmotnost_txt`}} \\[6pt]
  &= \frac{`r souciny_txt_maly[1]` + `r souciny_txt_maly[2]` + `r souciny_txt_maly[3]` + `r souciny_txt_maly[4]`}{\sqrt{`r soucet_ctvercu_txt` \cdot `r soucet_ctvercu_hmotnost_txt`}} \\[6pt]
  &= \frac{`r soucet_soucinu_maly_txt`}{\sqrt{`r soucet_ctvercu_txt` \cdot `r soucet_ctvercu_hmotnost_txt`}} \approx `r pearson_maly_txt`
\end{align}$$

```{r}
#| label: pripravit-pearsonuv-vysledek
#| echo: false
#| message: false
#| warning: false
pearson_korelace <-
  cor(
    data_tucnaci$delka_ploutve_mm,
    data_tucnaci$hmotnost_g,
    method = "pearson"
  )

pearson_korelace_txt <-
  formatC(
    pearson_korelace,
    format = "f",
    digits = 2,
    decimal.mark = ","
  )

pearson_from_cov_txt <-
  formatC(
    kovariance_val / (
      sd(data_tucnaci$delka_ploutve_mm) *
        sd(data_tucnaci$hmotnost_g)
    ),
    format = "f",
    digits = 2,
    decimal.mark = ","
  )
```

Tímto přepočtem dostaneme Pearsonův korelační koeficient, který leží vždy v intervalu $[-1, +1]$. Je to totéž jako vydělit kovarianci součinem směrodatných odchylek obou proměnných: kovariance i oba rozptyly obsahují stejné dělení počtem pozorování minus jedna, které se ve zlomku vykrátí. Pro čtyři tučňáky z tabulky vyšlo `r pearson_maly_txt`. V celém souboru tučňáků vychází přibližně `r pearson_from_cov_txt` — přesně stejný výsledek jako přímý výpočet (`r pearson_korelace_txt`), jen odvozený jinou cestou.

Pearsonův korelační koeficient označujeme symbolem $r$. Obecně jej můžeme zapsat bez jednotek takto:

$$r = \frac{\sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^{n}(x_i - \bar{x})^2} \sqrt{\sum_{i=1}^{n}(y_i - \bar{y})^2}}$$

Zde $r$ označuje Pearsonův korelační koeficient, $x_i$ a $y_i$ dvojici hodnot jednoho tučňáka, $\bar{x}$ a $\bar{y}$ průměry obou proměnných a $n$ počet pozorování. Čitatel shrnuje společný pohyb odchylek; jmenovatel jej převede na bezrozměrné měřítko.

## `r render_glossary_term("korelace", display = "Korelace")` jako stručný číselný popis

```{r}
#| label: pripravit-spearmanovu-korelaci
#| echo: false
#| message: false
#| warning: false
spearman_korelace <-
  cor(
  data_tucnaci$delka_ploutve_mm,
  data_tucnaci$hmotnost_g,
  method = "spearman"
)

spearman_korelace_txt <-
  formatC(
  spearman_korelace,
  format = "f",
  digits = 2,
  decimal.mark = ","
)

```

Graf je vždy první krok, ale někdy chceme vztah dvou `r render_glossary_term("numericka-promenna", display = "numerických proměnných")` shrnout i číslem. K tomu slouží korelace. Pearsonovo $r$ jsme před chvílí odvodili z `r render_glossary_term("kovariance", display = "kovariance")` — vydělíme kovarianci součinem `r render_glossary_term("smerodatna-odchylka", display = "směrodatných odchylek")` a dostaneme číslo v rozsahu $[-1, +1]$. V datech tučňáků vychází přibližně `r pearson_korelace_txt`.

Druhou běžnou možností je **Spearmanova korelace** (anglicky *Spearman's rank correlation*). Hodnoty obou proměnných nejprve nahradí jejich pořadím (nejkratší ploutev = 1, druhá nejkratší = 2, …) a korelaci spočítá z těchto pořadí. Ptá se tedy, zda vztah stále roste nebo stále klesá, i když body neleží na přímce; takovému vztahu se odborně říká monotónní. U tučňáků vychází Spearmanova korelace přibližně `r spearman_korelace_txt`.

V R si obě čísla spočítáme přímo nad stejnými dvěma sloupci, které jsme předtím použili pro `r render_glossary_term("bodovy-graf", display = "bodový graf")`:

```{r}
#| label: spocitat-korelaci-v-zakladnim-r
#| echo: true
#| eval: true
cor(
  # Pearsonova korelace shrnuje lineární vztah původních hodnot
  x = data_tucnaci$delka_ploutve_mm,
  y = data_tucnaci$hmotnost_g,
  method = "pearson"
)

cor(
  # Spearmanova korelace místo hodnot porovnává jejich pořadí
  x = data_tucnaci$delka_ploutve_mm,
  y = data_tucnaci$hmotnost_g,
  method = "spearman"
)
```

Takto dostaneme dvě stručná čísla, která shrnují tutéž biologickou otázku trochu jiným způsobem: Pearson zdůrazňuje lineární vztah, Spearman pořadí hodnot.

Na schématu níže jsou čtyři typické situace. Kladný vztah používá skutečná měření tučňáků; záporný vztah, zakřivený vztah a příklad bez zřetelné souvislosti jsou zkonstruované výukové ilustrace, nikoli další terénní měření. Všimněte si, jak se mění tvar bodového mraku a kdy se Pearsonova a Spearmanova korelace od sebe vzdálí.

```{r}
#| label: pripravit-ctyri-situace-korelace
#| echo: false
#| message: false
#| warning: false
# Stejný kladný příklad jako v prezentaci: skutečná data tučňáků
x_linear <- data_tucnaci$delka_ploutve_mm
y_linear <- data_tucnaci$hmotnost_g

# Stejný zakřivený příklad jako v prezentaci: zrychlující růst
# s úměrným šumem; pořadí hodnot roste, ale vztah není přímka.
set.seed(900723)
x_curve_n <- 120
x_curve <-
  sort(runif(x_curve_n, min = 1, max = 10))
signal_curve <- exp(x_curve / 2)
y_curve <-
  signal_curve * exp(rnorm(x_curve_n, mean = 0, sd = 0.2))

pearson_linear_txt <-
  formatC(
  cor(x_linear, y_linear, method = "pearson"),
  format = "f",
  digits = 2,
  decimal.mark = ","
)

spearman_linear_txt <-
  formatC(
  cor(x_linear, y_linear, method = "spearman"),
  format = "f",
  digits = 2,
  decimal.mark = ","
)

pearson_curve_txt <-
  formatC(
  cor(x_curve, y_curve, method = "pearson"),
  format = "f",
  digits = 2,
  decimal.mark = ","
)

spearman_curve_txt <-
  formatC(
  cor(x_curve, y_curve, method = "spearman"),
  format = "f",
  digits = 2,
  decimal.mark = ","
)

# Stejný příklad bez zřetelné souvislosti jako v prezentaci
x_none <-
  seq(1, 10, length.out = 12)
y_none <-
  c(3.2, 4.8, 2.7, 5.1, 3.9, 4.4, 2.9, 5.0, 3.5, 4.1, 3.8, 4.6)

pearson_none_txt <-
  formatC(
    cor(x_none, y_none, method = "pearson"),
    format = "f",
    digits = 2,
    decimal.mark = ","
  )

spearman_none_txt <-
  formatC(
    cor(x_none, y_none, method = "spearman"),
    format = "f",
    digits = 2,
    decimal.mark = ","
  )

# Stejný záporný příklad jako v prezentaci: nadmořská výška a teplota
x_neg <-
  seq(1, 10, length.out = 12)
y_neg <-
  c(6.0, 5.8, 5.5, 5.2, 4.7, 4.5, 4.2, 3.9, 3.4, 3.2, 2.8, 2.5)

set.seed(900723)
vec_noise <-
  rnorm(length(x_neg), mean = 0, sd = 1)
data_nadmorska <-
  data.frame(
    nadmorska_vyska_m = round((x_neg + vec_noise) * 90 + 200),
    prumerna_teplota_c = round(y_neg * 2.4 + 5.5, 1)
  )

x_neg <- data_nadmorska$nadmorska_vyska_m
y_neg <- data_nadmorska$prumerna_teplota_c

pearson_neg_txt <-
  formatC(
    cor(x_neg, y_neg, method = "pearson"),
    format = "f",
    digits = 2,
    decimal.mark = ","
  )

spearman_neg_txt <-
  formatC(
    cor(x_neg, y_neg, method = "spearman"),
    format = "f",
    digits = 2,
    decimal.mark = ","
  )
```

```{r}
#| label: korelace-schemata
#| echo: false
#| message: false
#| warning: false
#| fig-cap: "Čtyři typické situace korelace: silný kladný vztah tučňáků, žádná zřetelná souvislost, silný záporný vztah a stále rostoucí, ale zakřivený vztah. U každého panelu jsou hodnoty Pearsonovy a Spearmanovy korelace."

p_linear <-
  make_korelace_panel(
    x          = x_linear,
    y          = y_linear,
    line_layer = ggplot2::geom_smooth(
      method    = "lm",
      formula   = y ~ x,
      se        = FALSE,
      colour    = col_median_line,
      linewidth = 1.2
    ),
    label = paste0("Tučňáci: silná kladná\nPearson ", pearson_linear_txt, "\nSpearman ", spearman_linear_txt),
    x_label = "Délka ploutve (mm)",
    y_label = "Tělesná hmotnost (g)"
  )

p_none <-
  make_korelace_panel(
    x          = x_none,
    y          = y_none,
    line_layer = ggplot2::geom_hline(
      yintercept = mean(y_none),
      linetype   = "dashed",
      colour     = biostat_cols["grey_olive"],
      linewidth  = 0.8
    ),
    label = paste0("Bez zřetelné souvislosti\nPearson ", pearson_none_txt, "\nSpearman ", spearman_none_txt)
  )

p_neg <-
  make_korelace_panel(
    x          = x_neg,
    y          = y_neg,
    line_layer = ggplot2::geom_smooth(
      method    = "lm",
      formula   = y ~ x,
      se        = FALSE,
      colour    = col_median_line,
      linewidth = 1.2
    ),
    label = paste0("Silná záporná\nPearson ", pearson_neg_txt, "\nSpearman ", spearman_neg_txt),
    x_label = "Nadmořská výška (m)",
    y_label = "Průměrná teplota (°C)"
  )

p_curve <-
  make_korelace_panel(
    x          = x_curve,
    y          = y_curve,
    line_layer = ggplot2::geom_smooth(
      method = "lm",
      formula = y ~ x,
      se = FALSE,
      colour = col_median_line,
      linewidth = 1.2
    ),
    label      = paste0("Rostoucí, ale zakřivený\nPearson ", pearson_curve_txt, "\nSpearman ", spearman_curve_txt)
  )

cowplot::plot_grid(p_linear, p_none, p_neg, p_curve, ncol = 2)
```

- **Pearsonova korelace:** pracuje s původními hodnotami, shrnuje, jak těsně body leží podél přímky, a bývá citlivější na `r render_glossary_term("odlehla-pozorovani", display = "odlehlá pozorování")`.
- **Spearmanova korelace:** pracuje s pořadím hodnot, shrnuje, zda vztah stále roste nebo stále klesá, a bývá vůči odlehlým pozorováním odolnější.

Když body tvoří přibližně přímkový oblak bez výrazně vzdálených pozorování, bývají obě korelace podobné — tak je tomu u tučňáků. V zakřiveném panelu ale Spearmanova korelace zůstává téměř 1, zatímco Pearsonova je nižší: pořadí hodnot roste téměř dokonale, ale body neleží na přímce. Rozdíl mezi oběma čísly tak sám o sobě upozorňuje, že vztah není přímkový nebo že jej ovlivňuje několik vzdálených bodů.

V této lekci je důležité hlavně to, že obě čísla mají podpůrnou roli. Nenahrazují graf ani biologickou interpretaci. Pomáhají jen stručně vyjádřit, že delší ploutve bývají spojeny s vyšší hmotností a že tento vztah je v datech zřetelný i po převodu do číselné podoby.

::: {.box-data collapse="true" title="Doplňující: podobná čísla, různé grafy"}

### Stejná čísla, jiné grafy

Čtyři malé soubory známé jako **Anscombův kvartet** mají velmi podobné `r render_glossary_term("prumer", display = "průměry")`, `r render_glossary_term("rozptyl", display = "rozptyly")` i Pearsonovu korelaci. Nejdřív se podívejte na body: co byste o vztahu řekli v každém panelu, kdybyste znali jen společné číslo $r$?

```{r}
#| label: pripravit-anscombuv-kvartet
#| include: false
#| message: false
#| warning: false

data_anscombe <-
  tibble::tibble(
    soubor = rep(LETTERS[1:4], each = 11),
    x = c(anscombe$x1, anscombe$x2, anscombe$x3, anscombe$x4),
    y = c(anscombe$y1, anscombe$y2, anscombe$y3, anscombe$y4)
  )

souhrn_anscombe <-
  data_anscombe |>
  dplyr::group_by(soubor) |>
  dplyr::summarise(
    prumer_x = mean(x),
    prumer_y = mean(y),
    rozptyl_x = var(x),
    rozptyl_y = var(y),
    korelace = cor(x, y),
    .groups = "drop"
  )
```

```{r}
#| label: zobrazit-anscombuv-kvartet
#| echo: false
#| message: false
#| warning: false
#| fig-cap: "Anscombův kvartet: čtyři bodové grafy se skoro stejnými souhrnnými čísly, ale odlišným tvarem vztahu. Data jsou součástí R."
#| fig-alt: "Čtyři panely: přibližně lineární oblak, zakřivený vztah, lineární oblak s odlehlým bodem a téměř svislý shluk s jedním vzdáleným bodem."
#| fig-height: 5

ggplot2::ggplot(
  data = data_anscombe,
  mapping = ggplot2::aes(x = x, y = y)
) +
  ggplot2::geom_point(
    color = col_data,
    size = 2.6
  ) +
  ggplot2::facet_wrap(
    facets = ggplot2::vars(soubor),
    ncol = 2
  ) +
  ggplot2::labs(
    x = "Hodnota x",
    y = "Hodnota y"
  ) +
  theme_biostat(base_size = 12)
```


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

```{r}
#| label: zobrazit-souhrn-anscombova-kvartetu
#| echo: false
#| message: false
#| warning: false

souhrn_anscombe |>
  dplyr::transmute(
    Soubor = soubor,
    `Průměr x` = formatC(prumer_x, digits = 2, format = "f", decimal.mark = ","),
    `Průměr y` = formatC(prumer_y, digits = 2, format = "f", decimal.mark = ","),
    `Rozptyl x` = formatC(rozptyl_x, digits = 2, format = "f", decimal.mark = ","),
    `Rozptyl y` = formatC(rozptyl_y, digits = 2, format = "f", decimal.mark = ","),
    `Pearsonovo r` = formatC(korelace, digits = 2, format = "f", decimal.mark = ",")
  ) |>
  knitr::kable(align = c("l", rep("r", 5)))
```

```{=typst}
]
```


Čísla jsou téměř stejná, ale panely nejsou zaměnitelné: v jednom se vztah zakřivuje, jinde jej výrazně ovlivňuje jediný bod. Proto začínáme grafem a teprve pak čteme korelaci nebo další souhrn.

:::

## `r render_glossary_term("korelace", display = "Korelace")` není kauzalita

V této chvíli už umíme vztah dvou `r render_glossary_term("numericka-promenna", display = "numerických proměnných")` dobře popsat slovně i číslem. Je ale důležité nepřekročit hranici toho, co nám korelace opravdu říká. Korelace popisuje, že se dvě `r render_glossary_term("promenna", display = "proměnné")` mění společně, ale sama o sobě ještě neukazuje, která proměnná je příčinou druhé.

U tučňáků tedy můžeme bezpečně říct, že delší ploutev bývá spojena s vyšší hmotností. Nemůžeme ale jen z korelace uzavřít, že změna délky ploutve přímo způsobuje změnu hmotnosti. Část vztahu může vznikat i tím, že obě proměnné souvisejí s celkovou velikostí těla nebo s druhem.

Modelování umožní vztah popsat v původních jednotkách a zohlednit další sledované proměnné. Příčinný výklad však závisí na designu studie a předpokladech, které samotný graf, korelace ani přidání modelu nezaručí. Také korelační `r render_glossary_term("odhad", display = "odhad")` má `r render_glossary_term("nejistota", display = "nejistotu")`; její vyjádření není vyhrazené modelům a budeme se mu věnovat v další části kurzu.

::: {.box-question title="Co je správná interpretace?"}

- **A:** Kladná korelace znamená, že jsme prokázali příčinný vztah.
- **B:** Kladná korelace znamená, že proměnné spolu souvisejí, ale příčinu musíme řešit opatrněji.
- **C:** Korelace je jen technické číslo bez biologického významu.
:::

::: {.box-result collapse="true" title="Odpověď: korelace sama neurčuje příčinu"}

Správně je **B**. Korelace je biologicky užitečný souhrn společného vzoru v datech, ale sama ještě nerozlišuje přímou příčinu od dalších možných vysvětlení.
:::

## Jak typ `r render_glossary_term("promenna", display = "proměnných")` určuje výběr grafu

Typ grafu nevybíráme podle toho, co je „hezčí“, ale podle toho, jaký vztah mezi proměnnými chceme přečíst. Jinak čteme dvě `r render_glossary_term("numericka-promenna", display = "numerické proměnné")` na osách, jinak porovnáváme skupiny proti číslům a jinak sledujeme složení dvou `r render_glossary_term("kategoricka-promenna", display = "kategoriálních proměnných")`.

Jako rychlá pomůcka pomáhají dvě otázky:

- Jsou obě proměnné numerické?
- Nebo jedna z nich hlavně rozděluje `r render_glossary_term("pozorovani", display = "pozorování")` do skupin?

### Dvě numerické proměnné: délka ploutve a hmotnost

Právě proto jsme v této lekci začali `r render_glossary_term("bodovy-graf", display = "bodovým grafem")` (anglicky *scatter plot*). Když jsou obě proměnné numerické, bývá bodový graf nejlepší první volbou: každý bod nese dvojici hodnot a my můžeme najednou číst směr vztahu, rozptýlení bodů, možné shluky i `r render_glossary_term("odlehla-pozorovani", display = "odlehlá pozorování")`.

```{r}
#| label: dve-numericke-promenne
#| echo: false
#| message: false
#| warning: false
#| fig-cap: "Dvě numerické proměnné: bodový graf je přirozený první krok, protože každý bod nese dvojici hodnot současně."
ggplot2::ggplot(
  data = data_tucnaci,
  mapping = ggplot2::aes(
    x = delka_ploutve_mm,
    y = hmotnost_g
  )
) +
  ggplot2::geom_point(
    alpha = 0.65,
    colour = col_data,
    size = 2
  ) +
  ggplot2::labs(
    x = "Délka ploutve (mm)",
    y = "Tělesná hmotnost (g)"
  )
```

U tučňáků jsme v něm viděli, že delší ploutve bývají spojené s vyšší hmotností, ale také že vztah není dokonale těsný. Právě to je důvod, proč bodový graf nepřeskakovat. Žádné jediné číslo nám neukáže tvar oblaku bodů tak přímo.

Když ale chceme ve stejných datech zvýraznit jinou vlastnost vztahu, můžeme použít i několik dalších variant. Každá z nich ukazuje něco trochu jiného.

::: {.panel-tabset}

#### Hladký trend (*fitted trend line*)

Čára pomáhá rychle přečíst hlavní směr vztahu, ale body pořád zůstávají viditelné.

```{r}
#| label: galerie-dve-numericke-hladky-trend
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_tucnaci,
  mapping = ggplot2::aes(
    x = delka_ploutve_mm,
    y = hmotnost_g
  )
) +
  ggplot2::geom_point(
    alpha = 0.45,
    colour = col_data,
    size = 1.8
  ) +
  ggplot2::geom_smooth(
    method = "lm",
    formula = y ~ x,
    se = FALSE,
    colour = col_median_line,
    linewidth = 1.2
  ) +
  ggplot2::labs(
    x = "Délka ploutve (mm)",
    y = "Tělesná hmotnost (g)",
    title = "Bodový graf s vyznačeným hlavním trendem"
  )
```

#### Hustota v mřížce (*2D `r render_glossary_term("histogram", display = "histogram")`*)

Když je bodů hodně, můžeme místo jednotlivých bodů ukázat, kde se měření hromadí nejvíce.

```{r}
#| label: galerie-dve-numericke-hustota
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_tucnaci,
  mapping = ggplot2::aes(
    x = delka_ploutve_mm,
    y = hmotnost_g
  )
) +
  ggplot2::geom_bin2d(
    bins = 18
  ) +
  ggplot2::scale_fill_gradient(
    low = biostat_cols["light_gray"],
    high = biostat_cols["graphite"],
    name = "Počet\nbodů"
  ) +
  ggplot2::labs(
    x = "Délka ploutve (mm)",
    y = "Tělesná hmotnost (g)",
    title = "Mřížka hustoty místo jednotlivých bodů"
  )
```

#### Vrstevnice hustoty (*2D density*)

Místo počtů v mřížce můžeme ukázat hladké vrstevnice: čím tmavší oblast, tím víc tučňáků má podobnou kombinaci délky ploutve a hmotnosti.

```{r}
#| label: galerie-dve-numericke-vrstevnice
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_tucnaci,
  mapping = ggplot2::aes(
    x = delka_ploutve_mm,
    y = hmotnost_g
  )
) +
  ggplot2::stat_density_2d(
    mapping = ggplot2::aes(fill = ggplot2::after_stat(level)),
    geom = "polygon",
    bins = 7,
    alpha = 0.75,
    colour = "white",
    linewidth = 0.18
  ) +
  ggplot2::scale_fill_gradient(
    low = biostat_cols["light_gray"],
    high = biostat_cols["graphite"],
    guide = "none"
  ) +
  ggplot2::labs(
    x = "Délka ploutve (mm)",
    y = "Tělesná hmotnost (g)",
    title = "Vrstevnice hustoty bodů"
  )
```

#### Po druzích v panelech (*small multiples*)

Rozdělení do panelů pomáhá odlišit, co je společný trend a co je rozdíl mezi skupinami.

```{r}
#| label: galerie-dve-numericke-po-druzich
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_tucnaci,
  mapping = ggplot2::aes(
    x = delka_ploutve_mm,
    y = hmotnost_g,
    colour = druh_tucnaka
  )
) +
  ggplot2::geom_point(
    alpha = 0.65,
    size = 1.7
  ) +
  ggplot2::facet_wrap(
    facets = ggplot2::vars(druh_tucnaka),
    nrow = 1
  ) +
  ggplot2::scale_x_continuous(breaks = c(180, 200, 220)) +
  ggplot2::scale_colour_manual(
    values = barvy_druhu,
    guide = "none"
  ) +
  ggplot2::labs(
    x = "Délka ploutve (mm)",
    y = "Tělesná hmotnost (g)",
    title = "Stejný vztah rozdělený podle druhu"
  )
```

:::

### Kategoriální a numerická proměnná: druh a hmotnost

Jakmile nás začne zajímat přímo vztah mezi druhem a tělesnou hmotností, bodový graf už není nejpřirozenější první volbou. Teď totiž nesrovnáváme dvě číselné osy, ale kategoriální proměnnou určující skupiny a numerickou proměnnou, která se mezi nimi porovnává.

```{r}
#| label: hmotnost-podle-druhu
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data    = data_tucnaci,
  mapping = ggplot2::aes(
    x      = druh_tucnaka,
    y      = hmotnost_g,
    colour = druh_tucnaka,
    fill   = druh_tucnaka
  )
) +
  ggplot2::geom_boxplot(
    alpha         = 0.35,
    outlier.shape = NA
  ) +
  ggbeeswarm::geom_beeswarm(
    alpha = 0.55,
    size  = 1.4
  ) +
  ggplot2::scale_colour_manual(values = barvy_druhu, guide = "none") +
  ggplot2::scale_fill_manual(values = barvy_druhu, guide = "none") +
  ggplot2::labs(
    x = "Druh tučňáka",
    y = "Tělesná hmotnost (g)"
  )
```

Než si přečtete vysvětlení, zkuste si v grafu najít tři věci:

- Který druh leží typicky nejvýš?
- Kde se `r render_glossary_term("rozdeleni", display = "rozdělení")` hmotnosti mezi druhy překrývá a kde už méně?
- Je rozptýlení uvnitř jednotlivých druhů malé, nebo i v rámci jednoho druhu zůstává docela velké?


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

::: {.box-result collapse="true" title="Odpověď: rozdíly mezi druhy i uvnitř druhů"}

Tučňáci druhu Gentoo leží typicky nejvýš. Adelie a Chinstrap se vzájemně překrývají více. I uvnitř každého druhu je viditelné rozptýlení — takže ani znalost druhu sama o sobě nestačí k přesnému určení hmotnosti konkrétního jedince.
:::

```{=typst}
]
```


Na tomto grafu nečteme směr vztahu jako v bodovém grafu. Místo toho sledujeme, zda se rozdělení tělesné hmotnosti mezi druhy překrývá, nebo zda některé skupiny leží typicky výš či níž. I to je vztah dvou proměnných: hodnota jedné proměnné se mění podle úrovně druhé.

Na první pohled se zdá, že tučňáci druhu Gentoo leží typicky výš než Adelie a Chinstrap. Zároveň ale skupiny nejsou úplně oddělené a i uvnitř každého druhu zůstává viditelné rozptýlení. Pro tuto kombinaci typů proměnných je proto často přirozenější `r render_glossary_term("boxplot", display = "krabicový graf")` (anglicky *box plot*) nebo body rozdělené do skupin než klasický bodový graf.

Totéž srovnání ale můžeme nakreslit i jinak podle toho, zda chceme zdůraznit jednotlivé body, šířku skupin nebo jen jejich shrnutí.

::: {.panel-tabset}

#### Samotné body (*beeswarm plot*)

Nejlépe ukazují jednotlivá měření a překryv skupin, ale hůře shrnují jejich střed.

```{r}
#| label: galerie-kategorie-numericka-body
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_tucnaci,
  mapping = ggplot2::aes(
    x = druh_tucnaka,
    y = hmotnost_g,
    colour = druh_tucnaka
  )
) +
  ggbeeswarm::geom_beeswarm(
    alpha = 0.65,
    size = 1.6
  ) +
  ggplot2::scale_colour_manual(
    values = barvy_druhu,
    guide = "none"
  ) +
  ggplot2::labs(
    x = "Druh tučňáka",
    y = "Tělesná hmotnost (g)",
    title = "Skupiny zobrazené jen jako jednotlivé body"
  )
```

#### Krabice bez bodů (*box plot*)

Rychle shrnou polohu a překryv skupin, ale skryjí jemnější tvar rozložení uvnitř každého druhu.

```{r}
#| label: galerie-kategorie-numericka-krabice
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_tucnaci,
  mapping = ggplot2::aes(
    x = druh_tucnaka,
    y = hmotnost_g,
    fill = druh_tucnaka
  )
) +
  ggplot2::geom_boxplot(
    alpha = 0.45,
    width = 0.55,
    outlier.shape = NA
  ) +
  ggplot2::scale_fill_manual(
    values = barvy_druhu,
    guide = "none"
  ) +
  ggplot2::labs(
    x = "Druh tučňáka",
    y = "Tělesná hmotnost (g)",
    title = "Krabicový graf jako stručné shrnutí skupin"
  )
```

#### Šířka rozložení (*violin plot*)

Hladký obrys dává lepší představu o tom, kde je v každé skupině více a méně hodnot.

```{r}
#| label: galerie-kategorie-numericka-sirka
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_tucnaci,
  mapping = ggplot2::aes(
    x = druh_tucnaka,
    y = hmotnost_g,
    fill = druh_tucnaka
  )
) +
  ggplot2::geom_violin(
    alpha = 0.35,
    trim = FALSE,
    colour = NA
  ) +
  ggplot2::geom_boxplot(
    width = 0.14,
    alpha = 0.45,
    outlier.shape = NA,
    colour = biostat_cols["graphite"]
  ) +
  ggplot2::scale_fill_manual(
    values = barvy_druhu,
    guide = "none"
  ) +
  ggplot2::labs(
    x = "Druh tučňáka",
    y = "Tělesná hmotnost (g)",
    title = "Obrys skupin zdůrazňuje jejich vnitřní tvar"
  )
```

#### Počty v intervalech (*stacked histogram*)

Skládaný histogram ukazuje skutečné počty tučňáků v intervalech hmotnosti; barva rozlišuje druh.

```{r}
#| label: galerie-kategorie-numericka-histogram
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_tucnaci,
  mapping = ggplot2::aes(
    x = hmotnost_g,
    fill = druh_tucnaka
  )
) +
  ggplot2::geom_histogram(
    binwidth = 250,
    boundary = 2500,
    alpha = 0.7,
    position = "stack",
    colour = "white",
    linewidth = 0.35
  ) +
  ggplot2::scale_fill_manual(
    values = barvy_druhu,
    name = "Druh tučňáka"
  ) +
  ggplot2::labs(
    x = "Tělesná hmotnost (g)",
    y = "Počet tučňáků",
    title = "Skládaný histogram hmotnosti podle druhu"
  )
```

#### Relativní složení (*stacked density*)

Skládaná hustota neukazuje počty, ale relativní složení: který druh převažuje při dané hmotnosti.

```{r}
#| label: galerie-kategorie-numericka-skladana-hustota
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_tucnaci,
  mapping = ggplot2::aes(
    x = hmotnost_g,
    fill = druh_tucnaka
  )
) +
  ggplot2::geom_density(
    mapping = ggplot2::aes(y = ggplot2::after_stat(count)),
    alpha = 0.8,
    adjust = 1.25,
    position = "fill",
    colour = "white",
    linewidth = 0.35
  ) +
  ggplot2::scale_y_continuous(
    labels = scales::percent_format(accuracy = 1)
  ) +
  ggplot2::scale_fill_manual(
    values = barvy_druhu,
    name = "Druh tučňáka"
  ) +
  ggplot2::labs(
    x = "Tělesná hmotnost (g)",
    y = "Relativní složení",
    title = "Který druh převažuje při dané hmotnosti"
  )
```

:::

Další užitečný případ jsou dvě kategoriální proměnné najednou.

### Druh a ostrov

```{r}
#| label: pripravit-pocty-druhu-na-ostrovech
#| echo: false
#| message: false
#| warning: false
data_ctetnosti_ostrovy <-
  data_tucnaci |>
  dplyr::count(ostrov, druh_tucnaka, name = "pocet") |>
  tidyr::complete(ostrov, druh_tucnaka, fill = list(pocet = 0L)) |>
  dplyr::mutate(
    barva_popisku = dplyr::if_else(
      condition = pocet >= max(pocet) / 2,
      true = "white",
      false = biostat_cols[["graphite"]]
    )
  )

pocty_na_ostrovech <-
  tapply(
    data_ctetnosti_ostrovy$pocet,
    data_ctetnosti_ostrovy$ostrov,
    sum
  )
```

Jinou situaci dostaneme, když obě proměnné rozdělují data do skupin. V datech tučňáků je takovou dvojicí druh a ostrov, na kterém byl jedinec změřen. Teď nečteme číselné hodnoty na ose, ale počty a složení kombinací kategorií.

```{r}
#| label: druh-a-ostrov
#| echo: false
#| message: false
#| warning: false
#| fig-cap: "Složení druhů na třech ostrovech: každý sloupec představuje 100 % tučňáků změřených na daném ostrově."
ggplot2::ggplot(
  data = data_ctetnosti_ostrovy,
  mapping = ggplot2::aes(
    x = ostrov,
    y = pocet,
    fill = druh_tucnaka
  )
) +
  ggplot2::geom_col(
    position = "fill",
    colour = biostat_cols["graphite"]
  ) +
  ggplot2::scale_fill_manual(
    values = barvy_druhu,
    name = "Druh tučňáka"
  ) +
  ggplot2::scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  ggplot2::labs(
    x = "Ostrov",
    y = "Podíl",
    title = "Dvě kategoriální proměnné:\nsložení skupin podle podílů"
  )
```

U takového grafu nás zajímá něco jiného než u boxplotu. Neptáme se, jaká je typická hodnota, ale jaké je zastoupení druhů na jednotlivých ostrovech. Proto se zde hodí proporční skládaný `r render_glossary_term("sloupcovy-graf", display = "sloupcový graf")` (anglicky *proportional stacked bar chart*): každý ostrov má stejnou výšku 100 % a my snadno porovnáváme podíly.

Než si přečtete vysvětlení, zkuste odpovědět:

- Na kterém ostrově jsme zaznamenali jen jeden druh?
- Který druh se v datech vyskytuje na všech třech ostrovech?
- Co z tohoto grafu nevyčtete o tom, kolik tučňáků bylo na jednotlivých ostrovech změřeno?

::: {.box-result collapse="true" title="Odpověď: složení se mezi ostrovy liší, počty graf skrývá"}

Na ostrově Torgersen jsme zaznamenali jen druh Adelie, který jako jediný žije na všech třech ostrovech. Gentoo se v datech objevuje jen na ostrově Biscoe a Chinstrap jen na ostrově Dream. Proporční graf ale skrývá velikost skupin: všechny sloupce mají stejnou výšku, přestože na Biscoe bylo změřeno `r pocty_na_ostrovech[["Biscoe"]]` tučňáků a na Torgersen jen `r pocty_na_ostrovech[["Torgersen"]]`.
:::

Stejná dvojice kategoriálních proměnných se dá zobrazit i jinými způsoby podle toho, zda chceme zdůraznit podíly, absolutní počty nebo přehled celé tabulky najednou.

::: {.panel-tabset}

#### Proporce (*proportional stacked bar chart*)

Když mají všechny sloupce stejnou výšku, snadno porovnáváme složení ostrovů mezi sebou.

```{r}
#| label: galerie-dve-kategorie-proporce
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_ctetnosti_ostrovy,
  mapping = ggplot2::aes(
    x = ostrov,
    y = pocet,
    fill = druh_tucnaka
  )
) +
  ggplot2::geom_col(
    position = "fill",
    colour = biostat_cols["graphite"]
  ) +
  ggplot2::scale_fill_manual(
    values = barvy_druhu,
    name = "Druh tučňáka"
  ) +
  ggplot2::scale_y_continuous(
    labels = scales::percent_format(accuracy = 1)
  ) +
  ggplot2::labs(
    x = "Ostrov",
    y = "Podíl",
    title = "Proporční sloupce"
  )
```

#### Absolutní počty (*grouped bar chart*)

Vedle sebe postavené sloupce se hodí tehdy, když nás zajímá i velikost jednotlivých skupin, nejen jejich složení.

```{r}
#| label: galerie-dve-kategorie-pocty
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_ctetnosti_ostrovy,
  mapping = ggplot2::aes(
    x = ostrov,
    y = pocet,
    fill = druh_tucnaka
  )
) +
  ggplot2::geom_col(
    position = ggplot2::position_dodge(width = 0.75),
    width = 0.68
  ) +
  ggplot2::scale_fill_manual(
    values = barvy_druhu,
    name = "Druh tučňáka"
  ) +
  ggplot2::labs(
    x = "Ostrov",
    y = "Počet pozorování",
    title = "Sloupce vedle sebe"
  )
```

#### Tepelná mapa (*heat map*)

Když chceme rychle přečíst celou tabulku kombinací, je praktické převést ji na barevnou mřížku.

```{r}
#| label: galerie-dve-kategorie-tepelna-mapa
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_ctetnosti_ostrovy,
  mapping = ggplot2::aes(
    x = ostrov,
    y = druh_tucnaka,
    fill = pocet
  )
) +
  ggplot2::geom_tile(
    colour = "white",
    linewidth = 0.8
  ) +
  ggplot2::geom_text(
    mapping = ggplot2::aes(
      label = pocet,
      colour = barva_popisku
    ),
    size = 4.2,
    show.legend = FALSE
  ) +
  ggplot2::scale_fill_gradient(
    low = biostat_cols["light_gray"],
    high = biostat_cols["graphite"],
    name = "Počet"
  ) +
  ggplot2::scale_colour_identity() +
  ggplot2::labs(
    x = "Ostrov",
    y = "Druh tučňáka",
    title = "Tepelná mapa kombinací kategorií"
  )
```

#### Dlaždicový graf (*waffle chart*)

Každý čtvereček představuje jednoho tučňáka. Graf tak ukáže složení i velikost skupin najednou: ostrov s více změřenými tučňáky má více čtverečků.

```{r}
#| label: pripravit-dlazdice-po-jedincich
#| echo: false
#| message: false
#| warning: false
data_waffle_ostrovy <-
  data_ctetnosti_ostrovy |>
  tidyr::uncount(weights = pocet) |>
  dplyr::group_by(ostrov) |>
  dplyr::mutate(
    id_dlazdice = dplyr::row_number(),
    sloupec = ((id_dlazdice - 1) %% 14) + 1,
    radek = (id_dlazdice - 1) %/% 14
  ) |>
  dplyr::ungroup()
```

```{r}
#| label: galerie-dve-kategorie-waffle
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_waffle_ostrovy,
  mapping = ggplot2::aes(
    x = sloupec,
    y = radek,
    fill = druh_tucnaka
  )
) +
  ggplot2::geom_tile(
    colour = "white",
    linewidth = 0.5
  ) +
  ggplot2::facet_wrap(
    facets = ggplot2::vars(ostrov),
    nrow = 1
  ) +
  ggplot2::coord_equal() +
  ggplot2::scale_fill_manual(
    values = barvy_druhu,
    name = "Druh tučňáka"
  ) +
  ggplot2::scale_x_continuous(
    breaks = NULL
  ) +
  ggplot2::scale_y_continuous(
    breaks = NULL
  ) +
  ggplot2::labs(
    x = NULL,
    y = NULL,
    title = stringr::str_wrap("Dlaždicový graf: každý čtvereček je jeden tučňák", width = 42)
  ) +
  ggplot2::theme(
    panel.grid = ggplot2::element_blank(),
    axis.text = ggplot2::element_blank(),
    axis.ticks = ggplot2::element_blank()
  )
```

#### Bublinové počty (*bubble chart*)

Plocha každé bubliny odpovídá počtu tučňáků v kombinaci druhu a ostrova. Graf je efektní, ale porovnávat plochy je pro oko těžší než porovnávat délky sloupců.

```{r}
#| label: galerie-dve-kategorie-bubliny
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = dplyr::filter(data_ctetnosti_ostrovy, pocet > 0),
  mapping = ggplot2::aes(
    x = ostrov,
    y = druh_tucnaka,
    size = pocet,
    fill = druh_tucnaka
  )
) +
  ggplot2::geom_point(
    shape = 21,
    alpha = 0.8,
    colour = "white",
    stroke = 0.9
  ) +
  ggplot2::geom_text(
    mapping = ggplot2::aes(label = pocet),
    colour = "white",
    fontface = "bold",
    size = 3.4
  ) +
  ggplot2::scale_size_area(max_size = 20, guide = "none") +
  ggplot2::scale_fill_manual(values = barvy_druhu, guide = "none") +
  ggplot2::labs(
    x = "Ostrov",
    y = "Druh tučňáka",
    title = "Stejná tabulka jako plocha bublin"
  )
```

:::

Ještě jiný případ je seřazená osa proti numerické proměnné.

### Seřazené kategorie délky ploutve a hmotnost

```{r}
#| label: pripravit-kategorie-pro-vyber-grafu
#| echo: false
#| message: false
#| warning: false
hranice_delky_ploutve <-
  quantile(
    x = data_tucnaci$delka_ploutve_mm,
    probs = c(1 / 3, 2 / 3),
    na.rm = TRUE
  )

data_tucnaci_typy_grafu <-
  data_tucnaci |>
  dplyr::mutate(
    kategorie_delky_ploutve = cut(
      x = delka_ploutve_mm,
      breaks = c(-Inf, hranice_delky_ploutve[[1]], hranice_delky_ploutve[[2]], Inf),
      labels = c("kratší", "střední", "delší"),
      ordered_result = TRUE
    )
  )

# Ordered categories: a neutral light-to-dark ramp encodes order without
# reusing the semantic orange/purple layers or the species palette.
barvy_kategorii_delky <-
  stats::setNames(
    grDevices::colorRampPalette(
      c(biostat_cols[["light_gray"]], biostat_cols[["graphite"]])
    )(4)[2:4],
    c("kratší", "střední", "delší")
  )
```

Kategoriální proměnná může mít i přirozené pořadí. Abychom si to ukázali na stejných datech, rozdělíme délku ploutve do tří seřazených kategorií: kratší, střední a delší (každá obsahuje přibližně třetinu tučňáků). Jde jen o ukázku práce s `r render_glossary_term("ordinalni-promenna", display = "ordinální proměnnou")`; pro `r render_glossary_term("spojita-promenna", display = "spojitou proměnnou")` bychom vztah s hmotností četli raději přímo v bodovém grafu. Stejný typ grafu použijeme i pro přirozeně seřazené biologické škály, například malý–střední–velký.

Kategorii délky ploutve dáme na osu x a hmotnost jako numerickou proměnnou na osu y. V takové chvíli je užitečný opět krabicový graf (anglicky *box plot*), ale tentokrát je zásadní, že pořadí kategorií na ose nese informaci.

```{r}
#| label: hmotnost-podle-kategorie-delky-ploutve
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_tucnaci_typy_grafu,
  mapping = ggplot2::aes(
    x = kategorie_delky_ploutve,
    y = hmotnost_g,
    fill = kategorie_delky_ploutve
  )
) +
  ggplot2::geom_boxplot(
    alpha = 0.5,
    outlier.shape = NA
  ) +
  ggbeeswarm::geom_beeswarm(
    alpha = 0.35,
    size = 1.2,
    colour = biostat_cols["graphite"]
  ) +
  ggplot2::scale_fill_manual(
    values = barvy_kategorii_delky,
    guide = "none"
  ) +
  ggplot2::labs(
    x = "Kategorie délky ploutve",
    y = "Tělesná hmotnost (g)",
    title = stringr::str_wrap("Ordinální + numerická proměnná: záleží i na pořadí kategorií", width = 42)
  )
```

Než si přečtete vysvětlení, zkuste si všimnout tří věcí:

- Posouvá se typická hmotnost od kratších ploutví k delším?
- Kde se skupiny ještě překrývají a kde už méně?
- Co by se ztratilo, kdybychom pořadí kategorií na ose x promíchali?

Tentokrát neporovnáváme nepojmenované skupiny, ale seřazené stupně. Proto má smysl sledovat, zda se s rostoucí kategorií posouvá i typická hmotnost. Kdybychom pořadí kategorií na ose promíchali, velká část sdělení grafu by se ztratila.

I zde ale existuje více rozumných variant. Některé lépe ukazují jednotlivá pozorování, jiné zase co nejstručnější souhrn podle pořadí.

::: {.panel-tabset}

#### Samotné body (*beeswarm plot*)

Když chceme ponechat všechna měření viditelná, stačí je rozložit do stran v jednotlivých kategoriích.

```{r}
#| label: galerie-ordinalni-numericka-body
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_tucnaci_typy_grafu,
  mapping = ggplot2::aes(
    x = kategorie_delky_ploutve,
    y = hmotnost_g,
    colour = kategorie_delky_ploutve
  )
) +
  ggbeeswarm::geom_beeswarm(
    alpha = 0.55,
    size = 1.4
  ) +
  ggplot2::scale_colour_manual(
    values = barvy_kategorii_delky,
    guide = "none"
  ) +
  ggplot2::labs(
    x = "Kategorie délky ploutve",
    y = "Tělesná hmotnost (g)",
    title = "Všechna měření v seřazených kategoriích"
  )
```

#### Spojené středy (*line plot of group summaries*)

Pokud chceme co nejrychleji číst posun od kratších k delším kategoriím, pomůže jednoduché spojení typických hodnot.

```{r}
#| label: pripravit-mediany-serazenych-kategorii
#| echo: false
#| message: false
#| warning: false
data_typicka_hmotnost_delky <-
  data_tucnaci_typy_grafu |>
  dplyr::group_by(kategorie_delky_ploutve) |>
  dplyr::summarise(
    typicka_hmotnost_g = median(hmotnost_g),
    .groups = "drop"
  )
```

```{r}
#| label: galerie-ordinalni-numericka-stredy
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_typicka_hmotnost_delky,
  mapping = ggplot2::aes(
    x = kategorie_delky_ploutve,
    y = typicka_hmotnost_g,
    group = 1
  )
) +
  ggplot2::geom_line(
    colour = col_median_line,
    linewidth = 1.2
  ) +
  ggplot2::geom_point(
    colour = col_median_line,
    size = 3
  ) +
  ggplot2::labs(
    x = "Kategorie délky ploutve",
    y = "Typická hmotnost (g)",
    title = "Stručné shrnutí podle pořadí kategorií"
  )
```

#### Šířka rozložení (*violin plot*)

Hladký obrys pomůže ukázat, kde jsou v každé seřazené skupině data koncentrovaná a kde naopak řídká.

```{r}
#| label: galerie-ordinalni-numericka-sirka
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_tucnaci_typy_grafu,
  mapping = ggplot2::aes(
    x = kategorie_delky_ploutve,
    y = hmotnost_g,
    fill = kategorie_delky_ploutve
  )
) +
  ggplot2::geom_violin(
    alpha = 0.4,
    trim = FALSE,
    colour = NA
  ) +
  ggplot2::geom_boxplot(
    width = 0.14,
    alpha = 0.45,
    outlier.shape = NA,
    colour = biostat_cols["graphite"]
  ) +
  ggplot2::scale_fill_manual(
    values = barvy_kategorii_delky,
    guide = "none"
  ) +
  ggplot2::labs(
    x = "Kategorie délky ploutve",
    y = "Tělesná hmotnost (g)",
    title = "Obrys rozložení v seřazených kategoriích"
  )
```

:::

### Rychlá tabulka: který graf zvolit jako první?

Tabulka zahrnuje také souhrny jedné proměnné z předchozí lekce: histogram pro rozdělení číselných hodnot a sloupcový graf pro `r render_glossary_term("cetnost", display = "četnosti")` kategorií. U seřazených skupin můžeme porovnat `r render_glossary_term("median", display = "mediány")`; jednotlivá měření přitom nesmíme zaměnit za jejich souhrn.


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

```{r}
#| label: tabulka-vyber-grafu
#| echo: false
#| message: false
#| warning: false
#| results: asis
data_prehled_grafu <-
  tibble::tibble(
    co_chceme_zobrazit = c(
      "Jak vypadá rozdělení jedné numerické proměnné",
      "Kolik pozorování padá do jednotlivých skupin",
      "Jak spolu souvisejí dvě numerické proměnné",
      "Jak se liší numerická proměnná mezi skupinami",
      "Jak se kombinují dvě skupinové proměnné",
      "Jak se mění numerická proměnná napříč seřazenými stupni"
    ),
    typ_promennych = c(
      "jedna numerická",
      "jedna kategoriální",
      "numerická × numerická",
      "kategoriální × numerická",
      "kategoriální × kategoriální",
      "ordinální × numerická"
    ),
    dobry_prvni_graf = c(
      "histogram",
      "sloupcový graf",
      "bodový graf",
      "body ve skupinách; boxplot + body",
      "tabulka četností; sloupce počtů nebo podílů",
      "boxplot se seřazenou osou"
    ),
    co_v_nem_hledat = c(
      "tvar rozdělení, rozptýlení, odlehlá pozorování",
      "četnosti, podíly, nerovnováha skupin",
      "směr, tvar, rozptýlení, odlehlá pozorování",
      "rozdíly mezi skupinami, překryv, rozptýlení uvnitř skupin",
      "složení kombinací, počty vs. podíly",
      "posun mediánu, překryv, monotónní trend"
    )
  )

if (knitr::is_html_output()) {
  nahled_histogram <-
    kableExtra::spec_hist(
      x = c(4, 5, 5, 6, 6, 6, 7, 7, 8, 9),
      width = 180,
      height = 54,
      col = col_data,
      border = "white"
    )

  nahled_sloupce <-
    kableExtra::spec_plot(
      x = 1:4,
      y = c(4, 8, 5, 6),
      width = 180,
      height = 54,
      type = "h",
      lwd = 14,
      col = col_data
    )

  nahled_body <-
    kableExtra::spec_plot(
      x = c(1, 2, 2.4, 3, 3.8, 4.2, 5, 5.8, 6.2),
      y = c(1.1, 1.8, 2.4, 2.2, 3.3, 3.1, 3.9, 4.6, 4.9),
      width = 180,
      height = 54,
      type = "p",
      pch = 16,
      cex = 1.4,
      col = col_data
    )

  nahled_boxplot <-
    kableExtra::spec_boxplot(
      x = c(2.1, 2.4, 2.8, 3.1, 3.2, 3.4, 3.9, 4.5, 4.9),
      width = 180,
      height = 54,
      col = "#D9E6EA",
      border = col_data,
      medcol = col_median_line
    )

  nahled_skladany <-
    paste(
      "<svg xmlns='http://www.w3.org/2000/svg' width='180' height='54' viewBox='0 0 180 54'>",
      "<rect x='14' y='8' width='28' height='16' fill='darkorange'/>",
      "<rect x='14' y='24' width='28' height='11' fill='mediumpurple'/>",
      "<rect x='14' y='35' width='28' height='7' fill='cyan4'/>",
      "<rect x='56' y='8' width='28' height='9' fill='darkorange'/>",
      "<rect x='56' y='17' width='28' height='12' fill='mediumpurple'/>",
      "<rect x='56' y='29' width='28' height='13' fill='cyan4'/>",
      "<rect x='98' y='8' width='28' height='6' fill='darkorange'/>",
      "<rect x='98' y='14' width='28' height='10' fill='mediumpurple'/>",
      "<rect x='98' y='24' width='28' height='18' fill='cyan4'/>",
      "<line x1='8' y1='42.5' x2='134' y2='42.5' stroke='#6F6F6F' stroke-width='1.2'/>",
      "</svg>",
      sep = ""
    )

  nahled_ordinalni <-
    paste(
      "<svg xmlns='http://www.w3.org/2000/svg' width='180' height='54' viewBox='0 0 180 54'>",
      "<line x1='20' y1='14' x2='20' y2='44' stroke='#6F6F6F' stroke-width='1.5'/>",
      "<rect x='13' y='28' width='14' height='10' fill='#D9E6EA' stroke='#6F6F6F' stroke-width='1.2'/>",
      "<line x1='13' y1='33' x2='27' y2='33' stroke='#5D2890' stroke-width='1.5'/>",
      "<line x1='70' y1='10' x2='70' y2='44' stroke='#6F6F6F' stroke-width='1.5'/>",
      "<rect x='63' y='23' width='14' height='12' fill='#D9E6EA' stroke='#6F6F6F' stroke-width='1.2'/>",
      "<line x1='63' y1='29' x2='77' y2='29' stroke='#5D2890' stroke-width='1.5'/>",
      "<line x1='120' y1='8' x2='120' y2='44' stroke='#6F6F6F' stroke-width='1.5'/>",
      "<rect x='113' y='18' width='14' height='14' fill='#D9E6EA' stroke='#6F6F6F' stroke-width='1.2'/>",
      "<line x1='113' y1='25' x2='127' y2='25' stroke='#5D2890' stroke-width='1.5'/>",
      "<line x1='10' y1='44.5' x2='132' y2='44.5' stroke='#6F6F6F' stroke-width='1.2'/>",
      "</svg>",
      sep = ""
    )

  data_prehled_grafu |>
    dplyr::mutate(
      nahled = c(
        stringr::str_replace(nahled_histogram[["svg_text"]], "^<\\?xml[^>]*>\\n?", ""),
        stringr::str_replace(nahled_sloupce[["svg_text"]], "^<\\?xml[^>]*>\\n?", ""),
        stringr::str_replace(nahled_body[["svg_text"]], "^<\\?xml[^>]*>\\n?", ""),
        stringr::str_replace(nahled_boxplot[["svg_text"]], "^<\\?xml[^>]*>\\n?", ""),
        nahled_skladany,
        nahled_ordinalni
      )
    ) |>
    dplyr::select(
      `Typ proměnných` = typ_promennych,
      `Co chceme zobrazit` = co_chceme_zobrazit,
      `Dobrý první graf` = dobry_prvni_graf,
      `Co v něm hledat` = co_v_nem_hledat,
      Nahled = nahled
    ) |>
    knitr::kable(
      format = "html",
      escape = FALSE,
      col.names = c("Typ proměnných", "Co chceme zobrazit", "Dobrý první graf", "Co v něm hledat", "Náhled")
    ) |>
    kableExtra::kable_styling(
      full_width = TRUE,
      bootstrap_options = c("striped", "hover", "condensed")
    ) |>
    kableExtra::column_spec(1, width = "11em") |>
    kableExtra::column_spec(2, width = "16em") |>
    kableExtra::column_spec(3, width = "12em") |>
    kableExtra::column_spec(4, width = "15em") |>
    kableExtra::column_spec(5, width = "12em")
} else {
  data_prehled_grafu |>
    dplyr::transmute(
      `Typ proměnných` = typ_promennych,
      `Co chceme zobrazit` = co_chceme_zobrazit,
      `Dobrý první graf` = dobry_prvni_graf,
      `Co v něm hledat` = co_v_nem_hledat
    ) |>
    knitr::kable()
}
```

```{=typst}
]
```


Tahle tabulka není seznam všech možných grafů. Je to bezpečný první krok pro situaci, kdy si nejdřív potřebujeme ujasnit, co vlastně porovnáváme. Až potom má smysl přidávat další vrstvy, jemnější styl nebo formálnější model.

## Co zatím ještě nemůžeme tvrdit

Častý omyl zní: v grafu je vidět vztah, takže jsme prokázali příčinu. Graf ukázal směr a tvar vztahu a `r render_glossary_term("korelace", display = "korelace")` vyjádřila jeho těsnost číslem. Tím však neodlišíme přímý účinek jedné `r render_glossary_term("promenna", display = "proměnné")` od společné příčiny, například druhu nebo celkové velikosti těla. Pro tvrzení o příčině potřebujeme vhodný design studie a obhajitelné předpoklady; samotné přidání modelu příčinu neprokáže.

Správná formulace pro tuto chvíli je skromnější: v datech vidíme vzor, který je biologicky smysluplný a který chceme v dalších lekcích popsat přesněji.

## Co bude následovat dál

`r render_glossary_term("korelace", display = "Korelace")` již shrnula směr a těsnost vztahu číslem bez jednotek. V další lekci navážeme `r render_glossary_term("linearni-model", display = "lineárním modelem")`: vztah popíšeme v původních jednotkách a budeme se ptát, o kolik gramů se typická hmotnost liší při rozdílu jednoho milimetru délky ploutve. `r render_glossary_term("odhad", display = "Odhady")` a jejich `r render_glossary_term("nejistota", display = "nejistotu")` rozvineme v další části kurzu.

Palmer Penguins ale nabízí i další `r render_glossary_term("promenna", display = "proměnné")`, které se budou hodit pro pozdější příklady, například délku a hloubku zobáku. Stejná logika výběru grafu a interpretace vztahu se pak bude dát použít i na jiné biologické otázky než jen na délku ploutve a hmotnost.

![Schéma měření délky a hloubky zobáku v datech Palmer Penguins. Ilustrace: Allison Horst.](images/culmen_depth.png){width="85%"}

## Shrnutí

- Pro dvě `r render_glossary_term("numericka-promenna", display = "numerické proměnné")` je první volbou `r render_glossary_term("bodovy-graf", display = "bodový graf")`; nejprve jej dobře přečteme a teprve potom přidáváme další vrstvy informace.
- V bodovém grafu čteme současně směr vztahu i rozptýlení bodů ve dvou směrech.
- Biologická data jsou rozptýlená i při zřetelném vztahu; doplňková `r render_glossary_term("kategoricka-promenna", display = "kategoriální proměnná")`, například druh, může ukázat, že část této `r render_glossary_term("variabilita", display = "variability")` odráží skutečnou strukturu dat.
- `r render_glossary_term("korelace", display = "Korelace")` je stručný číselný popis vztahu dvou numerických proměnných, ale nenahrazuje graf ani interpretaci.
- Typ grafu vybíráme podle dvojice `r render_glossary_term("promenna", display = "proměnných")`; biologická otázka rozhodne, co má graf zvýraznit.
- Viditelný vzor ani korelace ještě nejsou důkazem příčiny.

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

Někdo vám ukáže jen číslo *r* = `r pearson_korelace_txt` pro délku ploutve a hmotnost tučňáků, bez grafu, a napíše titulek „Delší ploutev, těžší tučňák“. Jaké tři otázky mu položíte, než titulku uvěříte?

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

::: {.box-result collapse="true" title="Možné otázky"}

Jak vypadá bodový graf — leží body podél přímky, nebo tvoří křivku či oddělené shluky? Jsou v datech skupiny, například druhy, které mohou vztah vytvářet nebo měnit? Kolik `r render_glossary_term("pozorovani", display = "pozorování")` data obsahují a jak byla měřena? A tvrdí titulek jen souvislost, nebo i příčinu?
:::

```{=typst}
]
```

## 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
#| warning: false
#| message: 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/"
    )
  )
}
```
