---
title: "Jakou délku korunního lístku odhadneme z jeho šířky?"
subtitle: "Lineární model s jedním numerickým prediktorem"
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
    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/15
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")

# Visible code reads data/kosatce.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
  }
)

if (!require("tidyverse")) install.packages("tidyverse")
if (!require("ggbeeswarm")) install.packages("ggbeeswarm")
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("glossary")) install.packages("glossary")

library(tidyverse)
library(ggbeeswarm)
library(ggplot2)
library(glossary)
glossary::glossary_popup("hover")

# kableExtra's table helper uses jQuery; Quarto does not load it by default.
if (knitr::is_html_output()) {
  knitr::knit_meta_add(list(rmarkdown::html_dependency_jquery()))
}

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"))

# Glossary: use local slovnik/pojmy.yaml when available, otherwise remote fallback.
# 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 <- 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 = "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))
    }
  )
}
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 = biostat_cols[["white"]],
    def_color = biostat_cols[["graphite"]]
  )

}
```

## Úvod

V minulé lekci jsme popisovali vztah délky ploutve a hmotnosti tučňáků pomocí `r render_glossary_term("bodovy-graf", display = "bodového grafu")` a `r render_glossary_term("korelace", display = "korelace")`. Korelace řekla, jak těsně dvě měření souvisejí, ale neodpověděla na otázku „o kolik“. V této lekci tuto otázku položíme korunním lístkům kosatců: obě měření ponecháme v původních centimetrech a pomocí `r render_glossary_term("linearni-model", display = "lineárního modelu")` odhadneme délku lístku z jeho šířky.

::: {.box-question title="Hlavní otázka této lekce"}

**Mají širší korunní lístky kosatců také větší délku, a pokud ano, jak velká je tato změna?**
:::

Biologická otázka určuje role měření. Korelace se po prohození os nezmění; model pro odhad délky z šířky však řeší jinou otázku než model pro odhad šířky z délky. Délka korunního lístku je proto `r render_glossary_term("odezva", display = "odezva")` a šířka `r render_glossary_term("prediktor", display = "prediktor")`.

::: {.box-question title="Co čekáte ještě před prvním grafem?"}

- **A:** Širší lístky budou téměř vždy delší, a to u všech květů o stejný kus.
- **B:** Širší lístky budou v průměru delší, ale jednotlivé květy se budou od tohoto průměrného vztahu lišit.
- **C:** Šířka a délka korunního lístku spolu souviset nebudou.

Než budete číst dál, vyberte možnost, která vám připadá nejpravděpodobnější, a zdůvodněte ji.
:::

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

Biologicky nejpravděpodobnější je **možnost B**. Rozměry částí jednoho květu spolu obvykle souvisejí, ale květy se liší druhem, věkem i podmínkami, ve kterých rostly. Model proto popíše průměrný vztah a u jednotlivých květů zůstanou odchylky.
:::

## Výsledky učení

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

- převést biologickou otázku na vztah mezi `r render_glossary_term("odezva", display = "odezvou")` a `r render_glossary_term("prediktor", display = "prediktorem")`
- fitovat model `lm(odezva ~ prediktor)` a interpretovat směr i velikost `r render_glossary_term("efekt", display = "efektu")`
- vysvětlit `r render_glossary_term("residuum", display = "residuum")` jako rozdíl mezi naměřenou a `r render_glossary_term("odhadnuta-hodnota", display = "odhadnutou hodnotou")`
- číst graf residuí a rozpoznat v něm zbývající vzor

## Každý květ má dvě měření

Pro tuto lekci používáme reálná měření kosatců, která sesbíral botanik Edgar Anderson a později je publikoval Ronald Fisher. Jeden řádek tabulky odpovídá jednomu květu. Šířka a délka korunního lístku jsou `r render_glossary_term("numericka-promenna", display = "numerické proměnné")` v centimetrech; druh kosatce je `r render_glossary_term("kategoricka-promenna", display = "kategoriální proměnná")`.

V prvním modelu použijeme záměrně **jen jeden numerický** `r render_glossary_term("prediktor", display = "prediktor")` a **druh zatím do modelu nepřidáme**. Na konci se k tomuto rozhodnutí vrátíme při diagnostice.

::: {.box-data collapse="true" title="Doplňující: proč kosatce?"}

- měří biologicky srozumitelné veličiny ve stejných jednotkách
- obsahují přirozenou `r render_glossary_term("variabilita", display = "variabilitu")` mezi květy
- hodnoty šířky začínají blízko nuly, takže `r render_glossary_term("intercept", display = "intercept")` můžeme ukázat bez stlačení dat do malé části grafu
- umožňují plynulý přechod od `r render_glossary_term("bodovy-graf", display = "bodového grafu")` k `r render_glossary_term("linearni-model", display = "lineárnímu modelu")`

:::

Soubor [kosatce.csv](../data/kosatce.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_kosatce`; funkce `str()` ukáže názvy, typy a první hodnoty sloupců:

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

```{r}
#| label: data-kosatce-prehled
#| echo: true
# Načteme tabulku kosatců ze souboru CSV
data_kosatce <-
  read.csv(file = "data/kosatce.csv")

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

```{=typst}
]
```

Tabulka má čtyři sloupce:

- `sirka_listku`: šířka korunního lístku v centimetrech (`r render_glossary_term("prediktor", display = "prediktor")`)
- `delka_listku`: délka korunního lístku v centimetrech (`r render_glossary_term("odezva", display = "odezva")`)
- `druh`: druh kosatce, který zatím ponecháme mimo model
- `delka_kalisniho_listku`: délka kališního lístku v centimetrech; v této lekci ji nepoužijeme

Dále zkontrolujeme rozsahy hodnot, chybějící hodnoty a počet květů jednotlivých druhů:

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

```{r}
#| label: data-kosatce-kontrola
#| echo: true
# Minimum, kvartily, průměr a maximum numerických sloupců
summary(
  object = data_kosatce[, c("sirka_listku", "delka_listku")]
)

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

# Počet květů jednotlivých druhů
table(data_kosatce$druh)
```

```{=typst}
]
```

V datech je `r format_cz(nrow(data_kosatce))` květů ze `r format_cz(length(unique(data_kosatce$druh)))` druhů. Šířky se pohybují od `r format_cz(min(data_kosatce$sirka_listku))` do `r format_cz(max(data_kosatce$sirka_listku))` cm a délky od `r format_cz(min(data_kosatce$delka_listku))` do `r format_cz(max(data_kosatce$delka_listku))` cm; žádný sloupec nemá chybějící hodnoty. Jsou to malé kladné rozměry odpovídající květům, takže nic nenaznačuje chybu jednotek ani nesmyslnou hodnotu.

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

Měření korunních a kališních lístků tří druhů kosatců (*Iris setosa*, *I. versicolor* a *I. virginica*) sesbíral Edgar Anderson (Anderson 1935, *The irises of the Gaspé Peninsula*, Bulletin of the American Iris Society 59: 2–5). Ronald Fisher je použil v článku o rozlišování druhů pomocí několika měření ([Fisher 1936](https://doi.org/10.1111/j.1469-1809.1936.tb02137.x)). Data jsou součástí základní instalace R jako tabulka `iris` v balíčku `datasets`; podrobnosti uvádí [dokumentace datasetu](https://stat.ethz.ch/R-manual/R-patched/library/datasets/html/iris.html).

Soubor pro tuto lekci obsahuje všech `r format_cz(nrow(data_kosatce))` řádků tabulky `iris`. Žádný řádek ani naměřenou hodnotu jsme nezměnili; vybrali jsme jen čtyři sloupce a dali jim české názvy: `sirka_listku` (původně `Petal.Width`), `delka_listku` (`Petal.Length`), `druh` (`Species`) a `delka_kalisniho_listku` (`Sepal.Length`). Protože je R šířeno pod licencí GNU GPL verze 2 nebo 3, šíříme za stejných podmínek i tento výběr.

Soubor vytvořil skript [`R/prepare_l03_data.R`](https://github.com/CUNI-NATUR-Biostatistics/L03/blob/main/R/prepare_l03_data.R), který před uložením ověří počet řádků a absenci chybějících hodnot. Stejnou tabulku si můžete vytvořit i sami v základním R:

```{r}
#| label: ulozit-tabulku-kosatcu-do-csv
#| echo: true
#| eval: false
data_export <-
  data.frame(
    sirka_listku = iris$Petal.Width,
    delka_listku = iris$Petal.Length,
    druh = as.character(iris$Species),
    delka_kalisniho_listku = iris$Sepal.Length
  )
write.csv(
  x = data_export,
  file = "data/kosatce.csv",
  row.names = FALSE
)
```

Pro výuku používáme samostatný soubor CSV, protože v praktiku stačí načíst jediný soubor pomocí `read.csv()`.
:::

### Šířka a délka jednotlivých květů

Než budeme mluvit o modelu, podívejme se na jednotlivé květy.

Z předchozích lekcí už víme dvě základní věci:

- obě měřené `r render_glossary_term("promenna", display = "proměnné")` jsou numerické, tedy měřené na číselné škále
- každou proměnnou můžeme shrnout samostatně, ještě než se podíváme na jejich vztah

**Porovnejte rozsahy obou měření.** Stačí jejich samostatná shrnutí k určení toho, zda má širší lístek také větší délku? Doplňující grafy ukazují každou proměnnou samostatně; vztah vyžaduje zachování dvojic ze stejného květu.

::: {.box-data collapse="true" title="Doplňující: zopakování každé proměnné samostatně"}

::: {.panel-tabset}
## Šířka korunních lístků

**Jaké šířky jsou pro tyto květy běžné a jaké krajní?** Porovnejte minimum, `r render_glossary_term("median", display = "medián")` a maximum se zobrazenými měřeními.

```{r}
#| label: shrnuti-sirka-listku
#| echo: true
# Zkontrolujeme polohu a rozptýlení šířek korunních lístků
summary(object = data_kosatce$sirka_listku)
```

```{r}
#| label: graf-sirka-listku
#| fig-alt: "Houslový a krabicový graf šířek korunních lístků v cm."
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_kosatce,
  mapping = ggplot2::aes(x = "Šířka", y = sirka_listku)
) +
  ggplot2::geom_violin(
    trim = TRUE,
    fill = biostat_cols[["grey_olive"]],
    alpha = 0.25,
    color = biostat_cols[["graphite"]]
  ) +
  ggplot2::geom_boxplot(
    width = 0.12,
    outlier.shape = NA,
    fill = biostat_cols[["white"]],
    color = biostat_cols[["graphite"]]
  ) +
  ggplot2::labs(x = NULL, y = "Šířka korunního lístku (cm)") +
  theme_biostat()
```

## Délka korunních lístků

**Mají délky stejný rozsah jako šířky?** Obě veličiny měříme v cm; z rozdílných rozsahů samotných ještě nepoznáme vztah v jednotlivých květech.

```{r}
#| label: shrnuti-delka-listku
#| echo: true
# Zkontrolujeme polohu a rozptýlení délek korunních lístků
summary(object = data_kosatce$delka_listku)
```

```{r}
#| label: graf-delka-listku
#| fig-alt: "Houslový a krabicový graf délek korunních lístků v cm."
#| echo: false
#| message: false
#| warning: false
ggplot2::ggplot(
  data = data_kosatce,
  mapping = ggplot2::aes(x = "Délka", y = delka_listku)
) +
  ggplot2::geom_violin(
    trim = TRUE,
    fill = biostat_cols[["grey_olive"]],
    alpha = 0.25,
    color = biostat_cols[["graphite"]]
  ) +
  ggplot2::geom_boxplot(
    width = 0.12,
    outlier.shape = NA,
    fill = biostat_cols[["white"]],
    color = biostat_cols[["graphite"]]
  ) +
  ggplot2::labs(x = NULL, y = "Délka korunního lístku (cm)") +
  theme_biostat()
```
:::
:::

Jaký typ grafu dává největší smysl, když máme dvě `r render_glossary_term("numericka-promenna", display = "numerické proměnné")` měřené na stejném květu?

::: {.box-result collapse="true" title="Odpověď: pro dvě měření zachováme jejich dvojice"}

Dobrá volba je `r render_glossary_term("bodovy-graf", display = "bodový graf")`, protože chceme vidět jednotlivé květy a vztah mezi dvěma měřeními.
:::

### Jak spolu souvisejí šířka a délka?

**Prohlédněte si body:** roste délka s šířkou, je průběh přibližně přímý a leží všechny květy stejně blízko trendu? Svůj popis porovnejte s následujícími vodítky.

V našem případě je `r render_glossary_term("prediktor", display = "prediktorem")` **šířka korunního lístku** a podle ní vysvětlujeme rozdíly v **délce korunního lístku**.

*Co z grafu sledujte zejména:*

- **směr vztahu:** zda se s rostoucí šířkou obvykle zvětšuje i délka
- **neobvyklé květy:** jestli jsou v datech květy, které se chovají výrazně jinak než ostatní
- **rozptýlení kolem trendu:** jak daleko leží jednotlivé hodnoty od hlavního trendu


::: {.panel-tabset}

## Jednotlivé květy v `r render_glossary_term("bodovy-graf", display = "bodovém grafu")`

```{r}
#| label: jednoduse-bodovy-graf
#| fig-alt: "Každý bod je jeden květ, vodorovně šířka a svisle délka v cm; širší lístky jsou obvykle delší."
#| echo: true
#| fig-width: 8
#| fig-height: 5
plot(
  # Každý květ vytvoří jeden bod se dvěma naměřenými souřadnicemi
  x = data_kosatce$sirka_listku,
  y = data_kosatce$delka_listku,
  # Popisky vysvětlují význam a jednotky obou os
  xlab = "Šířka korunního lístku (cm)",
  ylab = "Délka korunního lístku (cm)",
  main = "Bodový graf šířky a délky korunního lístku",
  pch = 16,
  col = "#8A8A8A",
  cex = 1.4,
  cex.axis = 1.5,
  cex.lab = 1.5,
  cex.main = 1.5
)
```

## Graf pomocí {ggplot2}

```{r}
#| label: jednoduse-bodovy-graf-ggplot2
#| fig-alt: "Stejné dvojice šířek a délek kosatců jako v základním bodovém grafu."
#| echo: true
#| fig-width: 8
#| fig-height: 5

ggplot2::ggplot(
  # Data určují řádky; aes() přiřadí sloupce k osám
  data = data_kosatce,
  mapping = ggplot2::aes(x = sirka_listku, y = delka_listku)
) +
  # Jeden řádek dat se zobrazí jako jeden bod
  ggplot2::geom_point(
    size = 2,
    alpha = 0.8,
    color = "#8A8A8A"
  ) +
  # Popíšeme osy i celý graf
  ggplot2::labs(
    x = "Šířka korunního lístku (cm)",
    y = "Délka korunního lístku (cm)",
    title = "Bodový graf šířky a délky korunního lístku"
  ) +
  ggplot2::theme_bw(base_size = 18)
```

:::

### Která přímka vystihuje měření?

**Vyberte z grafu jednu přímku a zdůvodněte svou volbu.** Které vzdálenosti bodů od přímky by mělo pravidlo pro výběr zohlednit?

Pro stejná data existuje **více kandidátních přímek** a každá vede k **trochu jiné interpretaci**. Každá přímka má přitom jiné `r render_glossary_term("parametry-modelu", display = "parametry")` (viz níže).

```{r}
#| label: mnoho-primek-jednim-grafem
#| fig-alt: "Stejná měření kosatců se třemi různě umístěnými přímkami A, B a C."
#| echo: false
#| message: false
#| warning: false
# Every candidate is a model line, so all use the model colour; line type
#   and a direct label distinguish the alternatives.
kandidati_cary <-
  tibble::tibble(
    nazev = c("A", "B", "C"),
    a = c(0.5, 1.1, 1.7),
    b = c(1.8, 2.2, 2.6),
    typ_cary = c("dotted", "solid", "dashed")
  )

popisky_car <-
  tibble::tibble(
    nazev = kandidati_cary$nazev,
    x = rep(max(data_kosatce$sirka_listku) + 0.05, 3),
    y = kandidati_cary$a + kandidati_cary$b * (max(data_kosatce$sirka_listku) + 0.05)
  )

ggplot2::ggplot(
  data = data_kosatce,
  mapping = ggplot2::aes(x = sirka_listku, y = delka_listku)
) +
  ggplot2::geom_point(size = 2, alpha = 0.8, color = biostat_cols[["grey_olive"]]) +
  ggplot2::geom_abline(
    data = kandidati_cary,
    mapping = ggplot2::aes(intercept = a, slope = b, linetype = nazev),
    color = biostat_cols[["indigo_velvet"]],
    linewidth = 1,
    show.legend = FALSE
  ) +
  ggplot2::geom_label(
    data = popisky_car,
    mapping = ggplot2::aes(x = x, y = y, label = nazev),
    hjust = 0,
    color = biostat_cols[["indigo_velvet"]],
    fill = biostat_cols[["white"]],
    fontface = "bold",
    show.legend = FALSE
  ) +
  ggplot2::scale_linetype_manual(
    values = setNames(kandidati_cary$typ_cary, kandidati_cary$nazev)
  ) +
  ggplot2::labs(
    x = "Šířka korunního lístku (cm)",
    y = "Délka korunního lístku (cm)",
    title = "Mnoho možných přímek přes stejné body"
  ) +
  ggplot2::coord_cartesian(
    xlim = c(0, max(data_kosatce$sirka_listku) + 0.3),
    ylim = c(0, max(popisky_car$y) + 0.5)
  ) +
  theme_biostat()
```

Než rozhodneme, která přímka vystihuje data nejlépe, potřebujeme umět libovolnou přímku jednoznačně popsat.

## Jak popíšeme polohu přímky?

Každou přímku popíšeme dvěma `r render_glossary_term("parametry-modelu", display = "parametry")`:

- `a`: `r render_glossary_term("intercept", display = "intercept")` (anglicky *intercept*; česky také **průsečík**)
- `b`: `r render_glossary_term("sklon", display = "sklon")` (anglicky *slope*; česky také **regresní koeficient**)

Intuitivně:

- `a` **posouvá přímku nahoru nebo dolů** (mění základní úroveň `r render_glossary_term("odhadnuta-hodnota", display = "odhadnuté délky")`)
- `b` **mění strmost přímky**, tedy jak rychle roste odhadnutá délka s rostoucí šířkou

**Před prohlédnutím grafů odhadněte:** který parametr mění výšku přímky a který rozdíl mezi jejími odhady při dvou šířkách? V prvním grafu se mění pouze intercept, v druhém pouze sklon.

::: {.panel-tabset}
## Co se změní, když přímku posuneme?

```{r}
#| label: intercept-jednoduse
#| fig-alt: "Tři rovnoběžné přímky se liší průsečíkem se svislou osou, sklon zůstává stejný."
#| echo: false
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 4

data_intercepty <-
  data.frame(
    intercept = kandidati_cary$a,
    sklon = kandidati_cary$b[kandidati_cary$nazev == "B"],
    popisek = paste0("a = ", format_cz(kandidati_cary$a)),
    typ_cary = c("dotted", "solid", "dashed")
  )

ggplot2::ggplot() +
  ggplot2::geom_point(
    data = data_kosatce,
    mapping = ggplot2::aes(x = sirka_listku, y = delka_listku),
    size = 1.8,
    alpha = 0.6,
    color = biostat_cols[["grey_olive"]]
  ) +
  ggplot2::geom_abline(
    data = data_intercepty,
    mapping = ggplot2::aes(
      intercept = intercept,
      slope = sklon,
      linetype = popisek
    ),
    color = biostat_cols[["indigo_velvet"]],
    linewidth = 1,
    show.legend = FALSE
  ) +
  ggplot2::geom_text(
    data = data_intercepty,
    mapping = ggplot2::aes(
      x = 0,
      y = intercept,
      label = format_cz(intercept)
    ),
    hjust = 1.3,
    size = 5,
    fontface = "bold",
    color = biostat_cols[["indigo_velvet"]],
    show.legend = FALSE
  ) +
  ggplot2::scale_linetype_manual(
    values = setNames(
      data_intercepty$typ_cary,
      data_intercepty$popisek
    )
  ) +
  ggplot2::scale_x_continuous(
    limits = c(0, max(data_kosatce$sirka_listku) + 0.2),
    expand = ggplot2::expansion(mult = c(0, 0))
  ) +
  ggplot2::scale_y_continuous(
    breaks = data_intercepty$intercept,
    labels = NULL,
    limits = c(0, max(data_kosatce$delka_listku) + 0.5),
    expand = ggplot2::expansion(mult = c(0, 0))
  ) +
  ggplot2::labs(
    x = "Šířka korunního lístku (cm)",
    y = "Délka korunního lístku (cm)"
  ) +
  ggplot2::coord_cartesian(clip = "off") +
  theme_biostat() +
  ggplot2::theme(
    axis.title.y = ggplot2::element_text(
      margin = ggplot2::margin(r = 38)
    ),
    plot.margin = ggplot2::margin(5.5, 5.5, 5.5, 14)
  ) +
  ggplot2::guides(linetype = "none")
```


## Co se změní, když přímku nakloníme?

```{r}
#| label: sklon-jednoduse
#| fig-alt: "Tři přímky se společným interceptem a různými sklony."
#| echo: false
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 4
data_sklony <-
  data.frame(
    intercept = kandidati_cary$a[kandidati_cary$nazev == "B"],
    sklon = kandidati_cary$b,
    popisek = paste0("b = ", format_cz(kandidati_cary$b)),
    typ_cary = c("dotted", "solid", "dashed")
  )

ggplot2::ggplot() +
  ggplot2::geom_point(
    data = data_kosatce,
    mapping = ggplot2::aes(x = sirka_listku, y = delka_listku),
    size = 1.8,
    alpha = 0.6,
    color = biostat_cols[["grey_olive"]]
  ) +
  ggplot2::geom_abline(
    data = data_sklony,
    mapping = ggplot2::aes(
      intercept = intercept,
      slope = sklon,
      linetype = popisek
    ),
    color = biostat_cols[["indigo_velvet"]],
    linewidth = 1,
    show.legend = FALSE
  ) +
  ggplot2::geom_label(
    data = data_sklony,
    mapping = ggplot2::aes(
      x = max(data_kosatce$sirka_listku),
      y = intercept + sklon * max(data_kosatce$sirka_listku),
      label = popisek
    ),
    hjust = 1,
    color = biostat_cols[["indigo_velvet"]],
    fill = biostat_cols[["white"]],
    size = 4
  ) +
  ggplot2::scale_linetype_manual(
    values = setNames(data_sklony$typ_cary, data_sklony$popisek)
  ) +
  ggplot2::labs(
    x = "Šířka korunního lístku (cm)",
    y = "Délka korunního lístku (cm)"
  ) +
  ggplot2::coord_cartesian(
    xlim = c(0, max(data_kosatce$sirka_listku) + 0.2),
    ylim = c(0, 8)
  ) +
  theme_biostat()
```

:::

Obrázek níže ukazuje čtyři kombinace obou parametrů nad stejnými květy. V řádcích se mění intercept `a`, ve sloupcích sklon `b`.

```{r}
#| label: kombinace-a-b
#| fig-alt: "Čtyři kombinace interceptu a sklonu nad stejnými daty květů."
#| echo: false
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 6
data_kombinace_ab <-
  tidyr::crossing(
    a = kandidati_cary$a[c(1, 3)],
    b = kandidati_cary$b[c(1, 3)]
  )

ggplot2::ggplot(
  data = data_kosatce,
  mapping = ggplot2::aes(x = sirka_listku, y = delka_listku)
) +
  ggplot2::geom_point(
    size = 1.3,
    alpha = 0.55,
    color = biostat_cols[["grey_olive"]]
  ) +
  ggplot2::geom_abline(
    data = data_kombinace_ab,
    ggplot2::aes(intercept = a, slope = b),
    color = biostat_cols[["indigo_velvet"]],
    linewidth = 1
  ) +
  ggplot2::facet_grid(
    rows = ggplot2::vars(a),
    cols = ggplot2::vars(b),
    labeller = ggplot2::labeller(
      a = function(x) paste0("a = ", chartr(".", ",", x)),
      b = function(x) paste0("b = ", chartr(".", ",", x))
    )
  ) +
  ggplot2::labs(
    x = "Šířka korunního lístku (cm)",
    y = "Délka korunního lístku (cm)"
  ) +
  ggplot2::coord_cartesian(
    xlim = c(0, max(data_kosatce$sirka_listku) + 0.2),
    ylim = c(0, max(data_kosatce$delka_listku) + 0.5)
  ) +
  theme_biostat()
```

Každou kandidátní přímku tedy popíšou dvě čísla. Jak ale poznáme, která z nich nejlépe vystihuje vztah mezi šířkou a délkou korunního lístku?

::: {.box-rule collapse="true" title="Jak daleko leží naměřené délky od přímky?"}

Potřebujeme kritérium, které využije vzdálenosti naměřených délek od přímky při dané šířce. Proto potřebujeme `r render_glossary_term("odhadnuta-hodnota", display = "odhadnuté délky")` a `r render_glossary_term("residuum", display = "residua")`.
:::

## Čtyři květy a jedna zvolená přímka

```{r}
#| label: parametry-kandidatni-primky
#| echo: false
#| message: false
#| warning: false
a0 <- kandidati_cary$a[kandidati_cary$nazev == "B"]
b0 <- kandidati_cary$b[kandidati_cary$nazev == "B"]
```

### Čtyři květy napříč prostorem dat

Než porovnáme kandidátní přímky ve všech datech, ukážeme princip na čtyřech květech rozprostřených napříč grafem.

::: {.box-data title="Čtyři skutečné květy pro ruční výpočet"}

Vybrané řádky zachovávají dvojice měření stejného květu. Čísla řádků určují příklad; nejsou součástí statistické metody.

```{r}
#| label: vyber-4-bodu
#| echo: true
# Vybereme čtyři květy a ponecháme jen dvě modelované proměnné
data_ctyri_kvety <-
  data_kosatce[
    c(3, 82, 117, 101),
    c("sirka_listku", "delka_listku")
  ]

# Zobrazíme malou tabulku, se kterou budeme počítat ručně
data_ctyri_kvety
```

:::

::: {.box-data title="Konkrétní kandidátní přímka"}
Jako konkrétní příklad použijeme přímku **(B)** s `r render_glossary_term("intercept", display = "interceptem")` a = **`r format_cz(a0)`** a `r render_glossary_term("sklon", display = "sklonem")` b = **`r format_cz(b0)`**.
:::

`r render_glossary_term("residuum", display = "Residuum")` je **naměřená délka minus odhadnutá délka** v cm: svislý rozdíl mezi bodem a přímkou při stejné šířce. Úsečka znázorňuje jeho velikost, poloha bodu vůči přímce určuje znaménko. Pokud je **kladné**, model délku podhodnocuje; pokud **záporné**, délku nadhodnocuje.

Residua přitom **očekáváme i u užitečného modelu**. Jednotlivé organismy se v důsledku přirozené biologické `r render_glossary_term("variabilita", display = "variability")` liší, takže žádná jednoduchá přímka nepopíše každý květ přesně. **Samotná existence residuí proto neznamená, že je model špatný.**

V grafu níže jsou residua vidět přímo jako oranžové úsečky mezi očíslovanými šedými body a fialovou přímkou B.

```{r}
#| label: rezidua-graficky
#| fig-alt: "Čtyři očíslované květy a přímka B; oranžové svislé úsečky spojují naměřené a odhadnuté délky, květy 2 a 3 leží nad přímkou, květy 1 a 4 pod ní."
#| echo: false
#| message: false
#| warning: false
#| fig-width: 7.5
#| fig-height: 4.2

vec_odhadnuta_delka <-
  a0 + b0 * data_ctyri_kvety$sirka_listku

ggplot2::ggplot(
  data = tibble::tibble(
    kvet = seq_len(nrow(data_ctyri_kvety)),
    sirka_listku = data_ctyri_kvety$sirka_listku,
    delka_listku = data_ctyri_kvety$delka_listku,
    odhadnuta_delka = vec_odhadnuta_delka
  ),
  mapping = ggplot2::aes(x = sirka_listku, y = delka_listku)
) +
  ggplot2::geom_abline(
    intercept = a0,
    slope = b0,
    color = biostat_cols[["indigo_velvet"]],
    linewidth = 1
  ) +
  ggplot2::geom_segment(
    mapping = ggplot2::aes(xend = sirka_listku, yend = odhadnuta_delka),
    color = biostat_cols[["orange"]],
    linewidth = 1.4
  ) +
  ggplot2::geom_point(
    size = 3,
    color = biostat_cols[["grey_olive"]]
  ) +
  ggplot2::geom_text(
    mapping = ggplot2::aes(label = kvet),
    nudge_x = -0.08,
    fontface = "bold",
    color = biostat_cols[["graphite"]]
  ) +
  ggplot2::labs(
    x = "Šířka korunního lístku (cm)",
    y = "Délka korunního lístku (cm)"
  ) +
  theme_biostat()
```

### Jaký rozdíl zůstane u jednoho květu?

**Leží čtvrtý květ (vpravo nahoře) nad, nebo pod přímkou B?** Z grafu odhadněte znaménko rozdílu; výpočtem jej ověříme.

Pro jeden konkrétní květ použijeme tento postup:

1. Vezmeme **šířku** vybraného korunního lístku.
2. Dosadíme ji do vztahu **odhadnutá délka = intercept + sklon × šířka**.
3. Vypočítáme odhadnutou délku tohoto korunního lístku podle zvolené přímky.

Odhadnutou délku vyjadřuje tento vztah mezi veličinami:

$$
\begin{aligned}
\text{odhadnutá délka čtvrtého korunního lístku}
&=
\text{intercept} \\
&\quad +
\text{sklon}
\cdot
\text{šířka čtvrtého korunního lístku}.
\end{aligned}
$$

Šířka čtvrtého květu `r format_cz(data_ctyri_kvety$sirka_listku[4])` cm určuje **vodorovnou polohu bodu 4** v grafu výše. Intercept `r format_cz(a0)` a sklon `r format_cz(b0)` popisují **fialovou přímku**. Po dosazení:

$$
\begin{aligned}
\text{odhadnutá délka}_4
&= \text{intercept} + \text{sklon} \cdot \text{šířka}_4 \\
&= `r format_cz_math(a0)` + `r format_cz_math(b0)` \cdot `r format_cz_math(data_ctyri_kvety$sirka_listku[4])` \\
&= `r format_cz_math(a0)` + `r format_cz_math(round(b0 * data_ctyri_kvety$sirka_listku[4], 2))` \\
&= `r format_cz_math(round(vec_odhadnuta_delka[4], 2))`\ \text{cm}
\end{aligned}
$$

Výsledek `r format_cz(round(vec_odhadnuta_delka[4], 2))` cm je v grafu výše **výška fialové přímky** při šířce `r format_cz(data_ctyri_kvety$sirka_listku[4])` cm, tedy **horní konec čtvrté oranžové úsečky**. Naměřená délka `r format_cz(data_ctyri_kvety$delka_listku[4])` cm je naopak **výška bodu 4**. Residuum je rozdíl těchto dvou výšek:

$$
\begin{aligned}
\text{residuum čtvrtého korunního lístku}
&=
\text{naměřená délka čtvrtého korunního lístku} \\
&\quad -
\text{odhadnutá délka čtvrtého korunního lístku}.
\end{aligned}
$$

S hodnotami odečtenými z bodu a přímky:

$$
\begin{aligned}
\text{residuum}_4
&= \text{naměřená délka}_4 - \text{odhadnutá délka}_4 \\
&= `r format_cz_math(data_ctyri_kvety$delka_listku[4])` - `r format_cz_math(round(vec_odhadnuta_delka[4], 2))` \\
&= `r format_cz_math(round(data_ctyri_kvety$delka_listku[4] - vec_odhadnuta_delka[4], 2))`\ \text{cm}
\end{aligned}
$$

Hodnota `r format_cz(round(data_ctyri_kvety$delka_listku[4] - vec_odhadnuta_delka[4], 2))` cm odpovídá **čtvrté oranžové úsečce**. Samotná úsečka má délku `r format_cz(abs(round(data_ctyri_kvety$delka_listku[4] - vec_odhadnuta_delka[4], 2)))` cm; **záporné znaménko residua** říká, že bod 4 leží pod fialovou přímkou.

```{r}
#| label: odhadnuta-delka-jednoho-kvetu
#| echo: true
# Intercept (cm) a sklon (cm/cm) přímky B
a0 <- 1.1
b0 <- 2.2

# Dosadíme šířku čtvrtého květu do rovnice zvolené přímky
vec_odhadnuta_delka_4 <-
  a0 + b0 * data_ctyri_kvety$sirka_listku[4]

vec_odhadnuta_delka_4
```

```{r}
#| label: kontrola-shody-primky-b
#| include: false
stopifnot(
  identical(a0, kandidati_cary$a[kandidati_cary$nazev == "B"]),
  identical(b0, kandidati_cary$b[kandidati_cary$nazev == "B"])
)
```

### Naměřená délka, odhad a residuum čtyř květů

Stejný výpočet nyní použijeme pro všechny čtyři květy. Získáme tak:

- odhadnuté délky jednotlivých květů
- rozdíly mezi naměřenou a odhadnutou délkou (residua)

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

```{r}
#| label: odhadnute-delky-ctyr-kvetu
#| echo: true
vec_odhadnuta_delka <-
  a0 + b0 * data_ctyri_kvety$sirka_listku

# Residuum je naměřená hodnota minus odhad přímky
vec_residuum <-
  data_ctyri_kvety$delka_listku - vec_odhadnuta_delka

table_pred <-
  data.frame(
    kvet = seq_len(length.out = nrow(x = data_ctyri_kvety)),
    sirka = data_ctyri_kvety$sirka_listku,
    namerena_delka = data_ctyri_kvety$delka_listku,
    odhadnuta_delka = round(x = vec_odhadnuta_delka, digits = 2),
    residuum = round(x = vec_residuum, digits = 2)
  )

# Zobrazíme měření, odhad i residuum vedle sebe
table_pred
```

```{=typst}
]
```

*Čtení tabulky krok za krokem (všechny hodnoty jsou v cm):*

- sloupec `kvet`: pořadí květu, stejné jako čísla bodů v grafu
- sloupec `sirka`: šířka korunního lístku konkrétního květu
- sloupec `namerena_delka`: skutečně naměřená délka korunního lístku
- sloupec `odhadnuta_delka`: délka odhadnutá pomocí přímky B
- sloupec `residuum`: rozdíl mezi naměřenou a odhadnutou délkou

Květy 2 a 3 mají kladné residuum (leží nad přímkou), květy 1 a 4 záporné (leží pod ní).

Stejný výpočet můžeme pro každý květ zapsat pomocí symbolů. Index $i$ označuje pořadí květu. Naměřenou hodnotu označíme $y_i$, šířku $x_i$, odhadnutou hodnotu $\hat{y}_i$, intercept $a$, sklon $b$ a **residuum** (anglicky *residual*) $e_i$.

$$
\hat{y}_i = a + b \cdot x_i
$$

$$
e_i = y_i - \hat{y}_i
$$

V prvním vzorci je $\hat{y}_i$ odhadnutá délka jednoho lístku, $x_i$ jeho naměřená šířka, $a$ intercept a $b$ sklon. Ve druhém vzorci je $y_i$ skutečně naměřená délka a $e_i$ rozdíl mezi měřením a odhadem.

### Jak porovnat vzdálenosti všech květů od přímky?

Nabízí se residua jednoduše sečíst. Pro čtyři květy a přímku B ale vyjde:

$$
`r format_cz_math(vec_residuum[1], nsmall = 2)` + `r format_cz_math(vec_residuum[2], nsmall = 2)` + `r format_cz_math(vec_residuum[3], nsmall = 2)` + (`r format_cz_math(vec_residuum[4], nsmall = 2)`) = `r format_cz_math(sum(vec_residuum))`\ \text{cm}.
$$

Kladná a záporná residua se vyruší, přestože přímka žádným ze čtyř bodů neprochází. Prostý součet proto nepozná, jak daleko body od přímky leží.

`r render_glossary_term("soucet-ctvercu-residui", display = "Součet čtverců residuí")` (anglicky *sum of squared errors* nebo *residual sum of squares*, zkratka **SSE**) tento problém řeší: každé residuum nejprve umocní na druhou. Říká, jak moc jsou body daleko od modelové přímky. **Čím menší SSE, tím lépe přímka vystihuje data.**

Proč čtverce? Protože:

- **odstraní znaménko** (kladná a záporná residua se neodečtou)
- **více penalizují velké odchylky**
- **dávají jedno srovnatelné číslo** pro každou kandidátní přímku

Výpočet slovy pro tuto přímku:

$$
\text{součet čtverců residuí}
=
\left(\text{residuum prvního květu}\right)^2
+
\cdots
+
\left(\text{residuum posledního květu}\right)^2.
$$

S residui z tabulky výše:

$$
\begin{aligned}
\operatorname{SSE}
&= (`r format_cz_math(round(vec_residuum[1], 2), nsmall = 2)`)^2 + (`r format_cz_math(round(vec_residuum[2], 2), nsmall = 2)`)^2 \\
&\quad + (`r format_cz_math(round(vec_residuum[3], 2), nsmall = 2)`)^2 + (`r format_cz_math(round(vec_residuum[4], 2), nsmall = 2)`)^2 \\
&= `r format_cz_math(vec_residuum[1]^2, digits = 4)` + `r format_cz_math(vec_residuum[2]^2, digits = 4)` \\
&\quad + `r format_cz_math(vec_residuum[3]^2, digits = 4)` + `r format_cz_math(vec_residuum[4]^2, digits = 4)` \\
&= `r format_cz_math(sum(vec_residuum^2), digits = 4)`\ \text{cm}^2.
\end{aligned}
$$

Obecně pro všechna `r render_glossary_term("pozorovani", display = "pozorování")`:

$$
SSE = \sum_{i=1}^{n} e_i^2
$$

Zde $e_i$ označuje residuum jednoho květu, $e_i^2$ jeho čtverec, $n$ počet květů a znak $\sum$ znamená, že čtverce všech residuí sečteme.

Residua mají jednotku cm, jejich čtverce a SSE jednotku cm². Porovnávat SSE smíme pro stejné květy a stejné jednotky.

Stejný výpočet nyní provedeme pro každou kandidátní přímku a pro všech `r format_cz(nrow(data_kosatce))` květů. Výsledkem je jediné srovnatelné číslo pro každou přímku:

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

```{r}
#| label: sse-kandidatnich-primek
#| echo: false
#| output: true
#| message: false
#| warning: false
kandidati <-
  tibble::tibble(
    model = kandidati_cary$nazev,
    a = kandidati_cary$a,
    b = kandidati_cary$b
  ) |>
  dplyr::rowwise() |>
  dplyr::mutate(
    SSE = sum(
      (
        data_kosatce$delka_listku -
          (.data$a + .data$b * data_kosatce$sirka_listku)
      )^2
    )
  ) |>
  dplyr::ungroup()

knitr::kable(
  kandidati,
  digits = 2,
  col.names = c("Přímka", "Intercept (cm)", "Sklon (cm/cm)", "SSE (cm²)"),
  format.args = list(decimal.mark = ","),
  caption = "Porovnání kandidátních přímek ve všech datech"
)
```

```{=typst}
]
```

Která z nabízených přímek má nejmenší SSE? Vysvětlete, proč její SSE nelze přímo porovnat s předchozím součtem pro čtyři květy.

::: {.box-result collapse="true" title="Odpověď: nejlepší přímka má nejmenší SSE"}

Z nabízených přímek je nejlepší B, protože má nejmenší SSE. Tabulka používá všech `r format_cz(nrow(data_kosatce))` květů, zatímco ruční příklad pouze čtyři; jednotlivé součty proto hodnotí jiné množství měření. `r render_glossary_term("linearni-model", display = "Lineární model")` ale není omezen na tyto tři možnosti: hledá takový intercept a sklon, které SSE minimalizují mezi všemi možnými přímkami.
:::

### Naměřená délka se skládá z odhadu a residua

Naměřenou délku čtvrtého květu můžeme rozdělit na délku odhadnutou přímkou a residuum:

$$
\begin{aligned}
\text{naměřená délka čtvrtého korunního lístku}
&=
\text{odhadnutá délka čtvrtého korunního lístku} \\
&\quad +
\text{residuum čtvrtého korunního lístku}.
\end{aligned}
$$

S hodnotami čtvrtého květu:

$$
\begin{aligned}
\text{naměřená délka}_4
&= \text{odhadnutá délka}_4 + \text{residuum}_4 \\
&= `r format_cz_math(round(vec_odhadnuta_delka[4], 2))` + (`r format_cz_math(round(vec_residuum[4], 2), nsmall = 2)`) \\
&= `r format_cz_math(data_ctyri_kvety$delka_listku[4])`\ \text{cm}
\end{aligned}
$$

Pomocí již zavedených obecných symbolů:

$$
y_i = \hat{y}_i + e_i = a + b \cdot x_i + e_i
$$

Interpretace v biologickém příběhu:

- $y_i$ je naměřená délka jednoho korunního lístku a $x_i$ jeho naměřená šířka
- $\hat{y}_i$ je délka odhadnutá přímkou pro tuto šířku
- `a` (intercept) je odhad délky korunního lístku při nulové šířce; nula leží těsně za rozsahem pozorovaných dat
- `b` (sklon) říká, o kolik se v průměru liší odhadnutá délka dvou korunních lístků, jejichž šířky se liší o 1 cm
- $e_i$ je nevysvětlená část variability (residuum); index $i$ označuje pořadí květu

{{< pagebreak >}}

### O kolik se liší odhadnutá délka?

Sklon zvolené přímky spočítáme ze dvou bodů, které na této přímce leží; nestačí dvě libovolná pozorování z dat. Použijeme šířky prvního a čtvrtého květu:

```{r}
#| label: dva-body-na-primce
#| echo: false
#| output: true
#| message: false
#| warning: false
x_sklon <- data_ctyri_kvety$sirka_listku[c(1, 4)]
odhadnuta_delka_sklon <- a0 + b0 * x_sklon
b_sklon <-
  (odhadnuta_delka_sklon[2] - odhadnuta_delka_sklon[1]) /
    (x_sklon[2] - x_sklon[1])

knitr::kable(
  tibble::tibble(
    bod = c("při šířce květu 1", "při šířce květu 4"),
    sirka = x_sklon,
    odhadnuta_delka = round(odhadnuta_delka_sklon, 2)
  ),
  col.names = c("Bod na přímce", "Šířka (cm)", "Odhadnutá délka (cm)"),
  format.args = list(decimal.mark = ","),
  caption = "Dva body ležící na zvolené přímce"
)
```

Sklon přímky vyjadřuje vztah mezi dvěma změnami:

$$
\text{sklon přímky}
=
\frac{
  \text{změna odhadnuté délky korunního lístku}
}{
  \text{změna šířky korunního lístku}
}.
$$

Se souřadnicemi dvou bodů na přímce:

$$
\text{sklon přímky}
=
\frac{`r format_cz_math(round(odhadnuta_delka_sklon[2], 2))` - `r format_cz_math(round(odhadnuta_delka_sklon[1], 2))`}
{`r format_cz_math(x_sklon[2])` - `r format_cz_math(x_sklon[1])`}
=
\frac{`r format_cz_math(round(odhadnuta_delka_sklon[2] - odhadnuta_delka_sklon[1], 2))`}
{`r format_cz_math(round(x_sklon[2] - x_sklon[1], 2))`}
=
`r format_cz_math(round(b_sklon, 2))`
$$

Obecně:

$$
b = \frac{\Delta \hat{y}}{\Delta x}
$$

Znak $\Delta$ znamená změnu neboli rozdíl dvou hodnot. $\Delta \hat{y}$ je rozdíl odhadnutých délek a $\Delta x$ rozdíl šířek dvou bodů na přímce. Jednotka sklonu je cm délky na cm šířky (cm/cm); obě veličiny jsou délky, ale jejich biologické role zůstávají odlišné.


## Přímka pro délku korunního lístku v R

Funkce `lm()` fituje `r render_glossary_term("linearni-model", display = "lineární model")`. Její zápis si přečteme jako obecnou gramatiku:

```r
# Obecný zápis pojmenuje role proměnných a dat
lm(
  formula = odezva ~ prediktor,
  data = data
)
```

V tomto zápisu jsou `r render_glossary_term("odezva", display = "odezva")`, `r render_glossary_term("prediktor", display = "prediktor")` a `data` zástupné názvy. Znak `~` můžeme číst jako „odezvu vysvětlujeme pomocí prediktoru“.

Pro naši biologickou otázku dosadíme:

- odezva: `delka_listku`
- prediktor: `sirka_listku`
- data: `data_kosatce`

Konkrétní volání odděluje dvě operace:

- **`lm()` odhadne `r render_glossary_term("parametry-modelu", display = "parametry modelu")` z dat**
- **`coef(object = mod_listky)` vrátí dva parametry**, které už znáte z kandidátních přímek

```{r}
#| label: fit-jednoducheho-modelu
#| echo: true
# Odhadneme vztah mezi šířkou a délkou korunního lístku
mod_listky <-
  lm(
    formula = delka_listku ~ sirka_listku,
    data = data_kosatce
  )

# Zobrazíme odhadnutý intercept a sklon
coef(object = mod_listky)
```

Co vidíme za parametry (koeficienty) modelu:

- `r render_glossary_term("intercept", display = "Intercept")` odpovídá parametru `a`
- odhad u `sirka_listku` odpovídá `r render_glossary_term("sklon", display = "sklonu")` `b`

**Přečtěte odhad u `sirka_listku` v biologických jednotkách.** Jak se podle modelu liší `r render_glossary_term("odhadnuta-hodnota", display = "odhadnuté délky")` dvou květů, jejichž šířky se liší o 1 cm? Platí tento rozdíl přesně pro každou dvojici květů?

::: {.box-result collapse="true" title="Odpověď na biologickou otázku"}

Odhadnutý sklon je `r format_cz(round(coef(object = mod_listky)[["sirka_listku"]], 2))` cm/cm. Když porovnáme dva květy, jejichž korunní lístky se liší šířkou o 1 cm, model pro květ se širším lístkem odhaduje v průměru o `r format_cz(round(coef(object = mod_listky)[["sirka_listku"]], 2))` cm delší korunní lístek.

Pro každou konkrétní dvojici květů tento rozdíl přesně neplatí. Jednotlivé květy mají `r render_glossary_term("residuum", display = "residua")`, takže rozdíl jejich naměřených délek může být větší i menší. Sklon popisuje průměrný vztah, ne pravidlo pro každý květ.
:::

Odhadnutý intercept je `r format_cz(round(coef(object = mod_listky)[["(Intercept)"]], 2))` cm. Popisuje odhadnutou délku při nulové šířce, která leží těsně mimo rozsah pozorovaných dat. V tomto biologickém příběhu je proto důležitější sklon.

Model tak odpovídá na **směr a velikost vztahu**. Ani dobře vypadající lineární vztah ale automaticky neznamená **kauzální vztah**: z modelu zatím čteme pouze **asociaci mezi `r render_glossary_term("promenna", display = "proměnnými")`**.

### Přímka odhadnutá ze všech květů

Odhadnuté délky pro všechny květy vypočítáme pomocí `predict()`. Je to stejný výpočet jako ruční dosazení u čtyř květů, jen pro celý dataset:

```{r}
#| label: odhadnute-delky-vsech-kvetu
#| echo: true
vec_odhadnuta_delka <-
  predict(
    object = mod_listky,
    newdata = data_kosatce
  )

# Zkontrolujeme prvních šest odhadnutých délek
head(x = vec_odhadnuta_delka)
```

Graf naměřených hodnot s přímkou odhadnutou modelem ukazuje, že přímka opravdu odpovídá číslům z modelu:

::: {.panel-tabset}

### Základní graf v R

```{r}
#| label: graf-modelu-base-r
#| echo: true
#| fig-width: 8
#| fig-height: 5
#| fig-alt: "Všechna měření kosatců s přímkou ze dvou koeficientů modelu lm."
# Naměřené dvojice zůstávají body ve svých původních jednotkách
plot(
  x = data_kosatce$sirka_listku,
  y = data_kosatce$delka_listku,
  xlab = "Šířka korunního lístku (cm)",
  ylab = "Délka korunního lístku (cm)",
  pch = 16,
  col = "#8A8A8A",
  cex.axis = 1.4,
  cex.lab = 1.4
)

# Přímku určuje intercept a sklon odhadnutého modelu
abline(
  a = coef(object = mod_listky)[["(Intercept)"]],
  b = coef(object = mod_listky)[["sirka_listku"]],
  col = "#5D2890",
  lwd = 2
)
```

### Graf pomocí {ggplot2}

```{r}
#| label: graf-s-predict
#| fig-alt: "Všechna měření kosatců s fialovou přímkou odhadnutou pomocí lm."
#| echo: true
#| message: false
#| warning: false
data_model <-
  data.frame(
    # Původní naměřené hodnoty
    sirka_listku = data_kosatce$sirka_listku,
    delka_listku = data_kosatce$delka_listku,
    # Délky odhadnuté stejným modelem
    odhadnuta_delka = predict(
      object = mod_listky,
      newdata = data_kosatce
    )
  )

# Body jsou naměřené květy; čára spojuje odhady modelu
ggplot2::ggplot(
  data = data_model,
  mapping = ggplot2::aes(x = sirka_listku, y = delka_listku)
) +
  ggplot2::geom_point(
    size = 2,
    alpha = 0.75,
    color = "#8A8A8A"
  ) +
  ggplot2::geom_line(
    ggplot2::aes(y = odhadnuta_delka),
    color = "#5D2890",
    linewidth = 1
  ) +
  ggplot2::labs(
    x = "Šířka korunního lístku (cm)",
    y = "Délka korunního lístku (cm)"
  ) +
  ggplot2::theme_bw(base_size = 18)
```

:::

::: {.box-data collapse="true" title="Doplňující: kratší alternativa pomocí geom_smooth()"}

Balíček `{ggplot2}` umí lineární model odhadnout a nakreslit přímo pomocí `geom_smooth(method = "lm")`. Tuto zkratku používejte hlavně pro **první rychlou vizualizaci vztahu**. Vrstva sice model interně odhadne, ale nevytvoří snadno dostupný objekt modelu, ze kterého bychom mohli získat koeficienty, residua nebo další výsledky. Pro vlastní analýzu proto model fitujeme zvlášť pomocí `lm()` a uložíme jej do objektu. Argument `se = FALSE` zde vypíná zobrazení pásu `r render_glossary_term("nejistota", display = "nejistoty")`, ke kterému se vrátíme v další lekci.

```{r}
#| label: graf-modelu-geom-smooth
#| fig-alt: "Stejná modelová přímka vytvořená zkratkou geom_smooth."
#| echo: true
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 5

ggplot2::ggplot(
  # Použijeme původní měření šířky a délky
  data = data_kosatce,
  mapping = ggplot2::aes(x = sirka_listku, y = delka_listku)
) +
  ggplot2::geom_point(
    size = 2,
    alpha = 0.75,
    color = "#8A8A8A"
  ) +
  ggplot2::geom_smooth(
    # Funkce odhadne lineární model a rovnou nakreslí jeho přímku
    method = "lm",
    se = FALSE,
    color = "#5D2890"
  ) +
  ggplot2::labs(
    x = "Šířka korunního lístku (cm)",
    y = "Délka korunního lístku (cm)"
  ) +
  ggplot2::theme_bw(base_size = 18)
```
:::

## Co zůstává mezi měřením a odhadem?

`r render_glossary_term("residuum", display = "Residua")` jsou uložena přímo v objektu modelu. Získáme je funkcí `resid()` a uložíme do vektoru, se kterým můžeme dále pracovat:

```{r}
#| label: ziskat-rezidua-z-modelu
#| echo: true
vec_rezidua <-
  # resid() vrátí residuum pro každý květ použitý při fitování
  resid(object = mod_listky)

# Zkontrolujeme prvních šest residuí
head(x = vec_rezidua)
```

**Kladné residuum** znamená, že naměřená délka leží nad odhadem modelu. **Záporné residuum** znamená, že naměřená délka leží pod odhadem modelu.

### Jak jsou residua rozdělená?

`r render_glossary_term("histogram", display = "Histogram")` už znáte z předchozích lekcí. **Která znaménka jsou častá a jak daleko residua sahají od nuly?** Tentokrát na vodorovné ose nejsou původní délky lístků, ale rozdíly mezi naměřenou a `r render_glossary_term("odhadnuta-hodnota", display = "odhadnutou délkou")`:

```{r}
#| label: histogram-rezidui
#| fig-alt: "Histogram residuí v cm s přerušovanou čarou v nule."
#| echo: true
#| fig-width: 8
#| fig-height: 5
hist(
  # Histogram ukazuje rozdělení residuí
  x = vec_rezidua,
  breaks = 12,
  main = "Histogram residuí",
  xlab = "Residuum (cm)",
  ylab = "Četnost",
  col = "#D9D9D9",
  border = "white",
  cex.axis = 1.4,
  cex.lab = 1.4,
  cex.main = 1.4
)

abline(
  # Svislá čára označí nulové residuum
  v = 0,
  lty = 2,
  col = "#5D2890"
)
```

Přerušovaná čára označuje **nulové residuum**, tedy přesnou shodu naměřené a odhadnuté délky. Histogram ukazuje, **která residua jsou častá** a zda některé květy leží od odhadu neobvykle daleko. Neříká ale, **u kterých odhadnutých délek** tato residua vznikají.

### Mění se residua podél odhadnuté délky?

Proto ke každému květu spojíme odhadnutou délku z `fitted()` s jeho residuem z `resid()`:

```{r}
#| label: data-pro-graf-rezidui
#| echo: true
vec_odhadnuta_delka <-
  # fitted() vrátí odhadnutou délku každého květu
  fitted(object = mod_listky)

data_rezidua <-
  data.frame(
    odhadnuta_delka = vec_odhadnuta_delka,
    reziduum = vec_rezidua
  )

# Zkontrolujeme, že každý odhad má odpovídající residuum
head(x = data_rezidua)
```

Stejný graf ukazují dvě záložky: v základním R a pomocí balíčku `{ggplot2}`.

::: {.panel-tabset}

### Základní graf v R

```{r}
#| label: graf-rezidui-base-r
#| fig-alt: "Odhadnutá délka proti residuu každého květu; přerušovaná čára značí nulové residuum."
#| echo: true
#| fig-width: 8
#| fig-height: 5
plot(
  # Každý bod spojuje odhad modelu s residuem stejného květu
  x = data_rezidua$odhadnuta_delka,
  y = data_rezidua$reziduum,
  xlab = "Odhadnutá délka (cm)",
  ylab = "Residuum (cm)",
  main = "Odhadnuté délky vs. residua",
  pch = 16,
  col = "#8A8A8A",
  cex = 1.3,
  cex.axis = 1.4,
  cex.lab = 1.4,
  cex.main = 1.4
)

abline(
  # V dobře padnoucím modelu hledáme body rozptýlené kolem nuly
  h = 0,
  lty = 2,
  col = "#5D2890"
)
```

### Graf pomocí {ggplot2}

```{r}
#| label: graf-rezidui-ggplot2
#| fig-alt: "Stejné odhadnuté délky a residua v grafu ggplot2."
#| echo: true
#| fig-width: 8
#| fig-height: 5

ggplot2::ggplot(
  # Stejná data jako v předchozím grafu v základním R
  data = data_rezidua,
  mapping = ggplot2::aes(x = odhadnuta_delka, y = reziduum)
) +
  ggplot2::geom_hline(
    # Referenční čára ukazuje nulové residuum
    yintercept = 0,
    linetype = 2,
    color = "#5D2890"
  ) +
  ggplot2::geom_point(
    size = 2,
    alpha = 0.8,
    color = "#8A8A8A"
  ) +
  ggplot2::labs(
    x = "Odhadnutá délka (cm)",
    y = "Residuum (cm)",
    title = "Odhadnuté délky vs. residua"
  ) +
  ggplot2::theme_bw(base_size = 18)
```

:::

*Jak tento graf číst:*

- **každý bod** představuje jeden květ
- **vodorovná poloha** ukazuje délku odhadnutou modelem
- **svislá poloha** ukazuje rozdíl mezi naměřenou a odhadnutou délkou
- **bod na nulové čáře** má naměřenou délku přesně rovnou odhadu
- **bod nad nulou** model podhodnotil, **bod pod nulou** model nadhodnotil

Pro první kontrolu stačí jednoduchá otázka: *tvoří body kolem nuly neuspořádaný mrak, nebo v nich zůstává viditelný vzor?* **Viditelný vzor** znamená, že model v datech něco systematicky nezachytil.

**Najděte jednu oblast grafu, kde residua leží převážně nad nebo pod nulou.** Co to říká o odhadovaných délkách v této oblasti? Směr odchylky pojmenujte bez návrhu nového modelu; k biologickému vysvětlení se vrátíme s informací o druhu.

## Co musí být rozumné, aby přímka odpovídala na otázku?

`r render_glossary_term("linearni-model", display = "Lineární model")` je zjednodušený popis dat. Nemusí být dokonale pravdivý, ale jeho zjednodušení nesmí být v přímém rozporu s otázkou a daty. U jednoduchého modelu sledujeme čtyři odlišné oblasti:

1. **Tvar vztahu:** průměrná délka by se v pozorovaném rozsahu šířek měla měnit přibližně po přímce. Zřetelný oblouk v `r render_glossary_term("residuum", display = "residuích")` by upozornil, že lineární tvar nestačí.
2. **Rozptýlení kolem přímky:** svislý `r render_glossary_term("rozptyl", display = "rozptyl")` residuí by se neměl systematicky zvětšovat nebo zmenšovat. Trychtýřovitý tvar by znamenal, že model zachází se všemi částmi vztahu jako se stejně proměnlivými, přestože data tak nevypadají.
3. **Tvar residuí:** aby pozdější odhady `r render_glossary_term("nejistota", display = "nejistoty")` byly přiměřené, nesmí tvar residuí ukazovat silný rozpor s použitým modelem. Několik drobných odchylek od ideálního tvaru je v reálných datech běžných.
4. **Nezávislost pozorování:** zda jsou `r render_glossary_term("pozorovani", display = "pozorování")` nezávislá, vyplývá především ze způsobu získání dat. Graf residuí sám nepozná, že několik měření pochází ze stejného jedince, lokality nebo časové řady.

### Čtyři otázky, čtyři různé stopy v grafu

Následující schémata nejsou další čtyři datasety ani automatický test. Jsou to čtyři záměrně zjednodušené obrazy, které pomáhají pojmenovat, na jaký druh problému se díváme.

```{r}
#| label: pripravit-schematicke-vzory-residui
#| include: false

pocet_schematickych_bodu <- 40
poradi_schematickych_bodu <- seq_len(pocet_schematickych_bodu)
osa_schematickych_odhadu <- seq(
  from = -2,
  to = 2,
  length.out = pocet_schematickych_bodu
)
schematicky_sum <-
  0.30 * sin(poradi_schematickych_bodu * sqrt(2)) +
  0.16 * cos(poradi_schematickych_bodu * sqrt(3))

nazvy_schematickych_vzoru <- c(
  "Bez zjevného vzoru\nžádný nápadný rozpor",
  "Oblouk\nlineární tvar nestačí",
  "Trychtýř\nrozptýlení se mění",
  "Běhy v pořadí\nmožná závislost"
)

data_schematicke_vzory <-
  dplyr::bind_rows(
    tibble::tibble(
      vzor = nazvy_schematickych_vzoru[1],
      osa = osa_schematickych_odhadu,
      residuum = schematicky_sum
    ),
    tibble::tibble(
      vzor = nazvy_schematickych_vzoru[2],
      osa = osa_schematickych_odhadu,
      residuum =
        0.42 * (osa_schematickych_odhadu^2 - mean(osa_schematickych_odhadu^2)) +
        0.10 * sin(poradi_schematickych_bodu * sqrt(2))
    ),
    tibble::tibble(
      vzor = nazvy_schematickych_vzoru[3],
      osa = osa_schematickych_odhadu,
      residuum =
        (0.08 + 0.28 * (osa_schematickych_odhadu - min(osa_schematickych_odhadu))) *
        sin(poradi_schematickych_bodu * sqrt(5))
    ),
    tibble::tibble(
      vzor = nazvy_schematickych_vzoru[4],
      osa = poradi_schematickych_bodu,
      residuum =
        rep(c(-0.55, 0.55, -0.38, 0.38), each = 10) +
        0.08 * sin(poradi_schematickych_bodu)
    )
  ) |>
  dplyr::mutate(
    vzor = factor(vzor, levels = nazvy_schematickych_vzoru)
  )
```

```{r}
#| label: schematicke-vzory-residui
#| echo: false
#| fig-width: 8
#| fig-height: 6.2
#| fig-cap: "Schématické stopy v grafech residuí. Vodorovná čárkovaná čára označuje nulové residuum; spodní pravý panel jako jediný používá na vodorovné ose pořadí měření."
#| fig-alt: "Čtyři schematické bodové grafy residuí: neuspořádaný mrak kolem nuly, oblouk, trychtýř s rostoucím rozptýlením a dlouhé běhy bodů nad a pod nulou v pořadí měření."

ggplot2::ggplot(
  data = data_schematicke_vzory,
  mapping = ggplot2::aes(x = osa, y = residuum)
) +
  ggplot2::geom_hline(
    yintercept = 0,
    linetype = 2,
    color = biostat_cols[["indigo_velvet"]]
  ) +
  ggplot2::geom_point(
    size = 1.8,
    alpha = 0.85,
    color = biostat_cols[["grey_olive"]]
  ) +
  ggplot2::facet_wrap(
    facets = ggplot2::vars(vzor),
    ncol = 2,
    scales = "free_x"
  ) +
  ggplot2::labs(
    x = "Odhadnutá hodnota; vpravo dole pořadí měření",
    y = "Residuum",
    title = "Stejný typ grafu může ukázat různé stopy"
  ) +
  theme_biostat()
```

Vlevo nahoře nevidíme nápadný rozpor, ale ani tento obraz správnost modelu nedokazuje. Oblouk upozorňuje na tvar vztahu, trychtýř na měnící se rozptýlení. Dlouhé běhy nad a pod nulou se ukážou teprve tehdy, když zachováme pořadí měření; nezávislost proto posuzujeme společně s informací o tom, jak data vznikla.

### `r render_glossary_term("odezva", display = "Odezva")` a residua mohou mít jiný tvar

U kosatců můžeme porovnat dvě skutečná `r render_glossary_term("rozdeleni", display = "rozdělení")`. Levý `r render_glossary_term("histogram", display = "histogram")` popisuje naměřené délky korunních lístků. Pravý histogram popisuje odchylky těchto délek od `r render_glossary_term("odhadnuta-hodnota", display = "hodnot odhadnutých")` přímkou.

```{r}
#| label: pripravit-porovnani-odezvy-a-residui
#| include: false

data_tvar_odezvy_a_residui <-
  dplyr::bind_rows(
    tibble::tibble(
      co_zobrazujeme = "Naměřená délka korunního lístku",
      hodnota_cm = data_kosatce$delka_listku
    ),
    tibble::tibble(
      co_zobrazujeme = "Residuum kolem přímky",
      hodnota_cm = data_rezidua$reziduum
    )
  )
```

```{r}
#| label: porovnani-tvaru-odezvy-a-residui
#| echo: false
#| fig-width: 8
#| fig-height: 3.8
#| fig-cap: "Naměřená odezva a residua stejného modelu odpovídají na jiné otázky, a proto jejich rozdělení nemusí mít stejný tvar. Vodorovné osy panelů mají odlišný rozsah."
#| fig-alt: "Dva histogramy dat kosatců. Vlevo je rozdělení naměřené délky korunního lístku, vpravo rozdělení residuí kolem lineárního modelu."

ggplot2::ggplot(
  data = data_tvar_odezvy_a_residui,
  mapping = ggplot2::aes(x = hodnota_cm, fill = co_zobrazujeme)
) +
  ggplot2::geom_histogram(bins = 15, color = "white", linewidth = 0.3) +
  ggplot2::facet_wrap(
    facets = ggplot2::vars(co_zobrazujeme),
    nrow = 1,
    scales = "free_x"
  ) +
  ggplot2::scale_fill_manual(
    values = c(
      "Naměřená délka korunního lístku" = biostat_cols[["grey_olive"]],
      "Residuum kolem přímky" = biostat_cols[["orange"]]
    ),
    guide = "none"
  ) +
  ggplot2::labs(
    x = "Hodnota (cm)",
    y = "Počet květů",
    title = "Původní měření není totéž co odchylka od modelu"
  ) +
  theme_biostat()
```

Histogram odezvy se ptá, jaké délky jsme naměřili. Histogram residuí se ptá, jak se skutečné délky liší od přímky při dané šířce. Ani jeden histogram sám neposoudí linearitu, měnící se rozptýlení nebo nezávislost; k tomu potřebujeme předchozí grafy a znalost uspořádání studie.

`r render_glossary_term("diagnostika-modelu", display = "Diagnostika")` proto není zkouška, kterou model jednou provždy „udělá“ nebo „neudělá“. Spojuje grafy, znalost uspořádání studie a otázku, zda případný problém může změnit biologický závěr. V této lekci problémy pouze rozpoznáváme; jejich řešení přijdou postupně v dalších lekcích.

::: {.box-data collapse="true" title="Doplňující: předpoklady se týkají residuí kolem modelu"}

Lineární model nevyžaduje, aby samotná šířka nebo délka korunních lístků měla dokonale normální rozdělení. Důležité je, jak se chovají odchylky délky od hodnot očekávaných při dané šířce. Proto diagnostikujeme residua kolem přímky, nikoli pouze histogram původní odezvy.

Ani residua nemusí vytvořit dokonalou zvonovitou křivku. U malého souboru lze tvar posoudit jen hrubě; u velkého souboru zase i nepatrná odchylka může vypadat „statisticky zjistitelně“. Smyslem kontroly je hledat odchylku natolik výraznou, že zpochybňuje použitý model nebo jeho nejistotu, nikoli získat mechanické razítko správnosti.
:::

::: {.box-data collapse="true" title="Doplňující: Q–Q graf residuí"}

Tvar histogramu residuí závisí také na tom, kolik sloupců zvolíme. Podrobnější pohled na tvar residuí nabízí **`r render_glossary_term("qq-graf", display = "Q–Q graf")`** (z anglického *quantile–quantile*). Využívá `r render_glossary_term("kvantil", display = "kvantily")` – obecný název pro percentily a kvartily z první lekce, tedy hodnoty, pod nimiž leží určitý podíl seřazených dat. Residua seřadíme od nejmenšího po největší a každé porovnáme s hodnotou, kterou bychom na stejném místě pořadí čekali u souměrného zvonovitého (normálního) rozdělení.

::: {.panel-tabset}

#### Základní graf v R

```{r}
#| label: qq-graf-rezidui-base-r
#| fig-alt: "Seřazená residua proti kvantilům normálního rozdělení s kvartilovou referenční přímkou."
#| echo: true
#| fig-width: 8
#| fig-height: 5
qqnorm(
  # Seřazená residua porovnáme s kvantily normálního rozdělení
  y = vec_rezidua,
  main = "Q–Q graf residuí",
  xlab = "Teoretické kvantily normálního rozdělení",
  ylab = "Residuum (cm)",
  pch = 16,
  col = "#8A8A8A",
  cex = 1.3,
  cex.axis = 1.4,
  cex.lab = 1.4,
  cex.main = 1.4
)

qqline(
  # Referenční přímka prochází prvním a třetím kvartilem residuí
  y = vec_rezidua,
  lty = 2,
  lwd = 2,
  col = "#5D2890"
)
```

#### Graf pomocí {ggplot2}

```{r}
#| label: qq-graf-rezidui-ggplot2
#| fig-alt: "Stejný Q–Q graf residuí se zvýrazněnou kvartilovou přímkou."
#| echo: true
#| fig-width: 8
#| fig-height: 5

ggplot2::ggplot(
  # Stejná residua jako v předchozím grafu v základním R
  data = data_rezidua,
  mapping = ggplot2::aes(sample = reziduum)
) +
  ggplot2::stat_qq_line(
    # Referenční přímka prochází prvním a třetím kvartilem residuí
    linetype = 2,
    color = "#5D2890"
  ) +
  ggplot2::stat_qq(
    size = 2,
    alpha = 0.8,
    color = "#8A8A8A"
  ) +
  ggplot2::labs(
    x = "Teoretické kvantily normálního rozdělení",
    y = "Residuum (cm)",
    title = "Q–Q graf residuí"
  ) +
  ggplot2::theme_bw(base_size = 18)
```

:::

*Jak tento graf číst:*

- **každý bod** je residuum jednoho květu; **svislá osa** ukazuje jeho skutečnou hodnotu v cm
- **vodorovná osa** ukazuje, kde by residuum se stejným pořadím leželo u normálního rozdělení, měřeno ve `r render_glossary_term("smerodatna-odchylka", display = "směrodatných odchylkách")` od středu
- **přerušovaná přímka** ukazuje, kde by body ležely při dokonale zvonovitém tvaru
- **body blízko přímky:** tvar residuí se od zvonovitého rozdělení výrazně neliší
- **oba konce se odklánějí opačnými směry** (vlevo dolů, vpravo nahoru): velkých kladných i záporných residuí je víc, než by odpovídalo zvonovité křivce
- **body tvoří oblouk:** residua jsou nesouměrná, na jedné straně nuly sahají dál než na druhé

U kosatců leží prostřední residua téměř přesně na přímce. Na obou koncích se body mírně odklánějí opačnými směry: několik květů má o něco větší kladná i záporná residua, než by odpovídalo zvonovité křivce. Odchylka se týká hlavně nejkrajnějších bodů a je mírná; sama o sobě proto použití lineárního modelu nezpochybňuje.

Q–Q graf není test, který model „schválí“. I residua vzniklá přesně ze zvonovitého rozdělení se na koncích grafu trochu vlní, zvlášť u malého počtu pozorování. Hledáme proto výraznou a souvislou odchylku, ne každý bod mimo přímku.
:::

{{< pagebreak >}}

::: {.box-data collapse="true" title="Doplňující: odlehlý, krajní a vlivný bod"}

`r render_glossary_term("odlehla-pozorovani", display = "Odlehlé pozorování")` má neobvykle velkou odchylku od toho, co jsme čekali. **Krajní bod** má neobvyklou hodnotu `r render_glossary_term("prediktor", display = "prediktoru")`, například mimo většinu pozorovaných šířek. **Vlivný bod** je takový, jehož přítomnost výrazně mění odhadnutý `r render_glossary_term("intercept", display = "intercept")` nebo `r render_glossary_term("sklon", display = "sklon")`.

Jeden květ může splňovat jednu, dvě nebo všechny tři vlastnosti. Žádná z nich sama neznamená chybu ani povolení řádek smazat. Nejdříve ověříme měření a biologický kontext, potom zkoumáme citlivost závěru. Jak „krajnost“ a vliv jednotlivých květů změřit, ukazují dva následující doplňující rámečky.
:::

::: {.box-data collapse="true" title="Doplňující: pákový efekt"}

**`r render_glossary_term("pakovy-efekt", display = "Pákový efekt")`** (anglicky *leverage*) měří, jak neobvyklou hodnotu prediktoru má jedno pozorování. Přímku si můžeme představit jako páku: květ, jehož šířka leží daleko od průměrné šířky, ji může naklonit snáz než květ uprostřed rozsahu. V R pákový efekt získáme funkcí `hatvalues()`:

```{r}
#| label: pakovy-efekt-vypocet
#| echo: true
vec_pakovy_efekt <-
  # hatvalues() vrátí pákový efekt pro každý květ
  hatvalues(model = mod_listky)

# Spojíme měření, residua a pákový efekt jednotlivých květů
data_vliv_kvetu <-
  data.frame(
    sirka_listku = data_kosatce$sirka_listku,
    delka_listku = data_kosatce$delka_listku,
    reziduum = vec_rezidua,
    pakovy_efekt = vec_pakovy_efekt
  )

# Zkontrolujeme nejmenší a největší pákový efekt
range(vec_pakovy_efekt)
```

Pákový efekt nemá jednotky. S jediným prediktorem závisí pouze na šířce lístku: čím dál je šířka od průměrné šířky, tím je hodnota větší, a to na obě strany. Nejmenší hodnotu proto mají květy se šířkou blízko průměru. Naměřená délka ani residuum na pákový efekt vliv nemají. V následujícím grafu proto vyneseme obě vlastnosti z předchozího rámečku proti sobě: krajnost na vodorovnou osu a odlehlost na svislou.

::: {.panel-tabset}

#### Základní graf v R

```{r}
#| label: pakovy-efekt-base-r
#| fig-alt: "Pákový efekt proti residuu; velká šířka a velké residuum představují jiné vlastnosti květu."
#| echo: true
#| fig-width: 8
#| fig-height: 5
plot(
  # Vodorovně krajnost šířky, svisle odchylka od přímky
  x = data_vliv_kvetu$pakovy_efekt,
  y = data_vliv_kvetu$reziduum,
  xlab = "Pákový efekt",
  ylab = "Residuum (cm)",
  main = "Pákový efekt vs. residua",
  pch = 16,
  col = "#8A8A8A",
  cex = 1.3,
  cex.axis = 1.4,
  cex.lab = 1.4,
  cex.main = 1.4
)

abline(
  # Nulové residuum: naměřená délka je rovna odhadu
  h = 0,
  lty = 2,
  lwd = 2,
  col = "#5D2890"
)
```

#### Graf pomocí {ggplot2}

```{r}
#| label: pakovy-efekt-ggplot2
#| fig-alt: "Stejná srovnání pákového efektu a residuí v grafu ggplot2."
#| echo: true
#| fig-width: 8
#| fig-height: 5

ggplot2::ggplot(
  # Stejná data jako v předchozím grafu v základním R
  data = data_vliv_kvetu,
  mapping = ggplot2::aes(x = pakovy_efekt, y = reziduum)
) +
  ggplot2::geom_hline(
    # Nulové residuum: naměřená délka je rovna odhadu
    yintercept = 0,
    linetype = 2,
    color = "#5D2890"
  ) +
  ggplot2::geom_point(
    size = 2,
    alpha = 0.8,
    color = "#8A8A8A"
  ) +
  ggplot2::labs(
    x = "Pákový efekt",
    y = "Residuum (cm)",
    title = "Pákový efekt vs. residua"
  ) +
  ggplot2::theme_bw(base_size = 18)
```

:::

*Jak tento graf číst:*

- **vpravo** leží květy s nejvíce neobvyklou šířkou lístku
- **daleko od nulové čáry** leží květy, jejichž délka se nejvíce liší od odhadu
- body tvoří **svislé sloupce**: mnoho květů má stejnou šířku a stejný pákový efekt mají i květy stejně daleko od průměrné šířky na opačných stranách – velmi úzké i velmi široké lístky

Úplně vpravo leží květy s nejširšími lístky (`r format_cz(max(data_vliv_kvetu$sirka_listku))` cm). Všechny leží pod nulovou čarou: model jejich délku nadhodnocuje. Květ, který je zároveň vpravo a daleko od nuly, může přímku naklonit nejvíce. Právě tuto kombinaci shrnuje `r render_glossary_term("cookova-vzdalenost", display = "Cookova vzdálenost")`.
:::

::: {.box-data collapse="true" title="Doplňující: Cookova vzdálenost"}

**Cookova vzdálenost** (anglicky *Cook's distance*) spojuje obě vlastnosti z předchozího rámečku do jednoho čísla. Pro každý květ vyjadřuje, jak moc by se změnily odhadnuté délky všech květů, kdybychom právě tento květ vynechali a přímku odhadli znovu. Velkou hodnotu mívá hlavně květ, který má velké residuum a zároveň velký pákový efekt.

```{r}
#| label: cookova-vzdalenost-vypocet
#| echo: true
vec_cookova_vzdalenost <-
  # cooks.distance() vrátí Cookovu vzdálenost pro každý květ
  cooks.distance(model = mod_listky)

# Doplníme Cookovu vzdálenost a pořadí květu v datech
data_vliv_kvetu$cookova_vzdalenost <- vec_cookova_vzdalenost
data_vliv_kvetu$poradi_kvetu <-
  seq_len(length.out = nrow(x = data_vliv_kvetu))

# Orientační hodnota 4/n, kde n je počet květů v modelu
hranice_cook <- 4 / nrow(x = data_vliv_kvetu)
hranice_cook
```

::: {.panel-tabset}

#### Základní graf v R

```{r}
#| label: cookova-vzdalenost-base-r
#| fig-alt: "Cookova vzdálenost podle pořadí květu s orientační čarou 4/n."
#| echo: true
#| fig-width: 8
#| fig-height: 5
plot(
  # Každá svislá čárka představuje jeden květ
  x = data_vliv_kvetu$poradi_kvetu,
  y = data_vliv_kvetu$cookova_vzdalenost,
  type = "h",
  lwd = 2,
  col = "#8A8A8A",
  xlab = "Pořadí květu v datech",
  ylab = "Cookova vzdálenost",
  main = "Cookova vzdálenost květů",
  cex.axis = 1.4,
  cex.lab = 1.4,
  cex.main = 1.4
)

abline(
  # Orientační čára, nikoli hranice pro mazání dat
  h = hranice_cook,
  lty = 2,
  lwd = 2,
  col = "#5D2890"
)
```

#### Graf pomocí {ggplot2}

```{r}
#| label: cookova-vzdalenost-ggplot2
#| fig-alt: "Stejné Cookovy vzdálenosti a orientační čára v grafu ggplot2."
#| echo: true
#| fig-width: 8
#| fig-height: 5

ggplot2::ggplot(
  # Stejná data jako v předchozím grafu v základním R
  data = data_vliv_kvetu,
  mapping = ggplot2::aes(x = poradi_kvetu, y = cookova_vzdalenost)
) +
  ggplot2::geom_segment(
    # Každá svislá čárka vede od nuly k hodnotě jednoho květu
    mapping = ggplot2::aes(xend = poradi_kvetu, yend = 0),
    linewidth = 0.8,
    color = "#8A8A8A"
  ) +
  ggplot2::geom_hline(
    # Orientační čára, nikoli hranice pro mazání dat
    yintercept = hranice_cook,
    linetype = 2,
    color = "#5D2890"
  ) +
  ggplot2::labs(
    x = "Pořadí květu v datech",
    y = "Cookova vzdálenost",
    title = "Cookova vzdálenost květů"
  ) +
  ggplot2::theme_bw(base_size = 18)
```

:::

Nad orientační čarou leží `r format_cz(sum(vec_cookova_vzdalenost > hranice_cook))` z `r format_cz(length(vec_cookova_vzdalenost))` květů. Hodnota 4/n není test ani pravidlo pro mazání dat; pouze upozorňuje, které květy si zaslouží bližší pohled. Seřadíme je od největší Cookovy vzdálenosti:

```{r}
#| label: cookova-vzdalenost-tabulka
#| echo: true
data_vliv_serazene <-
  data_vliv_kvetu[
    # Seřadíme květy od největší Cookovy vzdálenosti
    order(
      data_vliv_kvetu$cookova_vzdalenost,
      decreasing = TRUE
    ),
  ]

# Vybereme sloupce, které chceme u květů porovnat
vec_sloupce_tabulky <-
  c(
    "sirka_listku",
    "delka_listku",
    "reziduum",
    "pakovy_efekt",
    "cookova_vzdalenost"
  )

# Zobrazíme pět nejvlivnějších květů se zaokrouhlenými hodnotami;
# čísla vlevo jsou pořadí květů v datech
round(
  x = head(x = data_vliv_serazene[, vec_sloupce_tabulky], n = 5),
  digits = 3
)
```

Nejvlivnější květ má široký lístek (`r format_cz(data_vliv_serazene$sirka_listku[1])` cm) a jeho naměřená délka leží `r format_cz(abs(round(data_vliv_serazene$reziduum[1], 2)))` cm pod odhadem přímky. Ke kterým druhům patří květy nad orientační čarou?

```{r}
#| label: cookova-vzdalenost-druhy
#| echo: true
# Spočítáme druhy květů s Cookovou vzdáleností nad orientační čarou
table(data_kosatce$druh[vec_cookova_vzdalenost > hranice_cook])
```

Všechny tyto květy patří druhu *virginica*, který má nejširší lístky, a tedy i velký pákový efekt. Model `r format_cz(sum(vec_rezidua[vec_cookova_vzdalenost > hranice_cook] < 0))` z nich nadhodnocuje a `r format_cz(sum(vec_rezidua[vec_cookova_vzdalenost > hranice_cook] > 0))` naopak výrazně podhodnocuje. Proč jedna společná přímka pro tři druhy nestačí, ukazuje kapitola *Které druhy model nadhodnocuje a podhodnocuje?*.

Jak moc nejvlivnější květ přímkou skutečně hýbe, zjistíme přímo: odhadneme stejný model bez něj a porovnáme koeficienty.

```{r}
#| label: citlivost-bez-vlivneho-kvetu
#| echo: true
# Najdeme řádek květu s největší Cookovou vzdáleností
radek_nejvlivnejsi <- which.max(x = vec_cookova_vzdalenost)

mod_bez_kvetu <-
  lm(
    # Stejný model odhadneme bez jediného květu
    formula = delka_listku ~ sirka_listku,
    data = data_kosatce[-radek_nejvlivnejsi, ]
  )

# Porovnáme intercept a sklon obou modelů
rbind(
  vsechny_kvety = coef(object = mod_listky),
  bez_nejvlivnejsiho = coef(object = mod_bez_kvetu)
)
```

Sklon se změní z `r format_cz(round(coef(object = mod_listky)[["sirka_listku"]], 2))` na `r format_cz(round(coef(mod_bez_kvetu)[["sirka_listku"]], 2))` cm délky na 1 cm šířky. Ani nejvlivnější květ tedy přímku podstatně neposune a závěr o vztahu šířky a délky se nemění. Ani vysoká Cookova vzdálenost ale sama o sobě není důvod řádek smazat. Je to důvod ověřit měření, zvážit biologický kontext a porovnat závěr s květem a bez něj – tak jako zde.
:::

## Které druhy model nadhodnocuje a podhodnocuje?

Ve všech dosavadních grafech vypadaly květy stejně a druh kosatce v modelu není. Barvu teď použijeme pouze jako pomůcku, abychom zjistili, zda se jednotlivé druhy v datech chovají podobně.

### Druhy v původních měřeních

Stejné barvy použijeme pro každý druh v obou následujících grafech:

```{r}
#| label: barvy-druhu-kosatcu
#| include: false
# Species are groups, not semantic data/model/focus roles: use the _brand
#   categorical palette (clr_skupiny) instead of the semantic colours.
barvy_druhu <-
  c(
    setosa = "#4E9A8A",
    versicolor = "#C2613A",
    virginica = "#2E6CA4"
  )
```

```{r}
#| label: puvodni-data-podle-druhu
#| fig-alt: "Šířka a délka kosatců s odlišenými druhy, které tvoří tři skupiny."
#| echo: false
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 5

ggplot2::ggplot(
  data = data_kosatce,
  mapping = ggplot2::aes(
    x = sirka_listku,
    y = delka_listku,
    color = druh
  )
) +
  ggplot2::geom_point(size = 2.3, alpha = 0.8) +
  ggplot2::scale_color_manual(values = barvy_druhu) +
  ggplot2::labs(
    x = "Šířka korunního lístku (cm)",
    y = "Délka korunního lístku (cm)",
    color = "Druh"
  ) +
  theme_biostat()
```

**Porovnejte polohu tří barevných skupin.** Shrnuje společná přímka pouze změny uvnitř jednoho druhu, nebo také rozdíly mezi druhy?

Barvy ukazují, že **tři druhy netvoří jeden společný bodový mrak**. Jedna přímka vedená přes všechny květy proto shrnuje nejen vztah šířky a délky, ale také **rozdíly mezi druhy**.

### Druhy v grafu `r render_glossary_term("odhadnuta-hodnota", display = "odhadnutých délek")` a `r render_glossary_term("residuum", display = "residuí")`

Ke každému residuu doplníme druh stejného květu:

```{r}
#| label: pridat-druh-k-reziduim
#| include: false
data_rezidua_druhy <-
  data_rezidua |>
  dplyr::mutate(
    druh = data_kosatce$druh
  )
```

```{r}
#| label: graf-rezidui-podle-druhu
#| fig-alt: "Residua s odlišenými druhy, jejichž skupiny leží vůči nule různě."
#| echo: false
#| message: false
#| warning: false
#| fig-width: 8
#| fig-height: 5

ggplot2::ggplot(
  data = data_rezidua_druhy,
  mapping = ggplot2::aes(
    x = odhadnuta_delka,
    y = reziduum,
    color = druh
  )
) +
  ggplot2::geom_hline(
    yintercept = 0,
    linetype = 2,
    color = biostat_cols[["grey_olive"]]
  ) +
  ggplot2::geom_point(size = 2.3, alpha = 0.8) +
  ggplot2::scale_color_manual(values = barvy_druhu) +
  ggplot2::labs(
    x = "Odhadnutá délka (cm)",
    y = "Residuum (cm)",
    color = "Druh"
  ) +
  theme_biostat()
```

**Které barevné skupiny leží častěji nad nulou a které pod ní?** Pro každý druh rozhodněte, zda model jeho délky spíše podhodnocuje, nebo nadhodnocuje. Použijte stejnou konvenci naměřená minus odhadnutá jako u jednoho květu.

```{r}
#| label: podil-rezidui-nad-nulou
#| echo: false
vec_podil_nad_nulou <-
  tapply(
    X = data_rezidua_druhy$reziduum > 0,
    INDEX = data_rezidua_druhy$druh,
    FUN = mean
  )
```

::: {.box-result collapse="true" title="Odpověď: model jednotlivé druhy odhaduje různě"}

U *versicolor* leží nad nulou `r format_cz(100 * vec_podil_nad_nulou[["versicolor"]], digits = 0)` % residuí, takže model jejich délky spíše **podhodnocuje**. U *setosa* leží pod nulou `r format_cz(100 * (1 - vec_podil_nad_nulou[["setosa"]]), digits = 0)` % residuí a model jejich délky spíše **nadhodnocuje**. U *virginica* leží pod nulou `r format_cz(100 * (1 - vec_podil_nad_nulou[["virginica"]]), digits = 0)` % residuí; nadhodnocení se týká hlavně květů s nejširšími lístky vpravo v grafu.
:::

Po obarvení vidíme, že residua jednotlivých druhů **neleží kolem nuly úplně stejně**. Model některé druhy častěji **nadhodnocuje** a jiné častěji **podhodnocuje**. Informace o druhu tedy zůstala v části dat, kterou jednoduchý model nevysvětlil.

V této lekci stačí toto omezení rozpoznat a popsat. Druh jako další `r render_glossary_term("prediktor", display = "prediktor")` přidáme do modelu až v některé z navazujících lekcí.

## Co přímka popisuje a co zůstává otevřené?

Jednoduchý `r render_glossary_term("linearni-model", display = "lineární model")` je výborný začátek, ale ne univerzální řešení. V této lekci jsme odhalili jeden důvod k opatrnosti: v `r render_glossary_term("residuum", display = "residuích")` zůstaly rozdíly mezi druhy kosatců.

Z grafů již umíme rozpoznat varovné vzory:

- zakřivený průběh může ukazovat na `r render_glossary_term("nelinearita", display = "nelinearitu")`
- postupné rozšiřování nebo zužování bodového mraku souvisí s pojmem `r render_glossary_term("heteroskedasticita", display = "heteroskedasticita")`
- jednotlivé body daleko od ostatních mohou být `r render_glossary_term("odlehla-pozorovani", display = "odlehlá pozorování")`
- barevně oddělené skupiny mohou upozornit na důležitý `r render_glossary_term("prediktor", display = "prediktor")`, který v modelu chybí

*V této lekci stačí vzor graficky rozpoznat a popsat. Formální testy ani nápravu modelu zatím nevyžadujeme.* Později podle biologické otázky zvážíme další prediktor, `r render_glossary_term("transformace", display = "transformaci")`, jiný tvar vztahu nebo jiný model.

**Hlavní princip zůstává stejný:** nejprve omezení modelu transparentně pojmenujeme a teprve potom volíme vhodné rozšíření.

## Co bude následovat dál

V následující lekci ponecháme stejný modelový rámec a přidáme první inferenční vrstvu: `r render_glossary_term("vyberova-variabilita", display = "výběrovou variabilitu")` a `r render_glossary_term("interval-spolehlivosti", display = "intervaly spolehlivosti")` jako vyjádření `r render_glossary_term("nejistota", display = "nejistoty")` odhadu. Testování hypotéz a `r render_glossary_term("p-hodnota", display = "p-hodnota")` přijdou až v pozdější lekci.

## Shrnutí

- Biologickou otázku převádíme na vztah mezi `r render_glossary_term("odezva", display = "odezvou")` a `r render_glossary_term("prediktor", display = "prediktorem")`; u kosatců je délka korunního lístku odezva a šířka prediktor.
- Funkce `lm()` najde přímku s nejmenším `r render_glossary_term("soucet-ctvercu-residui", display = "součtem čtverců residuí")`. Její `r render_glossary_term("sklon", display = "sklon")` udává, o kolik se v průměru liší odezva dvou `r render_glossary_term("pozorovani", display = "pozorování")`, jejichž prediktor se liší o jednu jednotku.
- `r render_glossary_term("residuum", display = "Residuum")` je rozdíl mezi naměřenou a `r render_glossary_term("odhadnuta-hodnota", display = "odhadnutou hodnotou")`; residua nejsou „odpad“, ale informace o tom, co model nevysvětluje.
- V grafu residuí hledáme neuspořádaný mrak kolem nuly; zbývající vzor, například oddělené skupiny druhů, upozorňuje na omezení modelu.

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

Někdo vám ukáže jen odhadnutý sklon `r format_cz(round(coef(object = mod_listky)[["sirka_listku"]], 2))` cm/cm pro šířku a délku korunních lístků kosatců, bez grafu, a napíše titulek „O centimetr širší lístek, o `r format_cz(round(coef(object = mod_listky)[["sirka_listku"]], 2))` cm delší lístek“. 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á `r render_glossary_term("bodovy-graf", display = "bodový graf")` a graf residuí — tvoří residua neuspořádaný mrak kolem nuly, nebo v nich zůstává vzor? Jsou v datech skupiny, například druhy, které vztah vytvářejí nebo mění? Kolik pozorování data obsahují a jak byla měřena? Jak přesný je odhad sklonu? A tvrdí titulek jen průměrný rozdíl mezi květy, nebo i to, že širší lístek délku způsobuje?
:::

```{=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
#| 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/"
    )
  )
}
```