---
title: "Jak poznáme, že je model dobrý?"
subtitle: "Podkladové čtenářské materiály k přednášce"
lang: cs
filters:
  - ../theme/semantic-boxes.lua
format:
  html:
    toc: true
    toc-depth: 3
    toc-title: "Obsah"
    number-sections: true
    code-fold: false
    code-summary: "Zobrazit kód v R"
    embed-resources: true
    theme: [flatly, ../theme/skripta_theme.scss]
    highlight-style: github
    include-in-header:
      - "../theme/fonts-include.html"
      - "../theme/skripta-favicon.html"
  typst:
    toc: true
    toc-depth: 3
    number-sections: true
    output-file: "skripta_raw.pdf"
    include-in-header:
      - "skripta_theme.typ"
      - "skripta-logo.typ"
    date: last-modified
date: 2026/09/30
date-format: long
date-modified: last-modified
---

```{r}
#| label: nastaveni
#| include: false
here::i_am("Learning_materials/skripta.qmd")

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

knitr::opts_chunk$set(
  message = FALSE,
  warning = FALSE,
  fig.width = 7,
  fig.height = 4.6
)

# Hidden figures and tables use Czech decimal commas. 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
  }
)

# Network add-ons must not block a local render.
options(timeout = 15)

source(here::here("R", "set_r_theme.R"))
source(here::here("R", "format_value_cz.R"))
source(here::here("R", "Functions", "use_ragg_brand_fonts.R"))
use_ragg_brand_fonts()

# Glossary: use local slovnik/pojmy.yaml when available, otherwise remote.
local_pojmy <- here::here("..", "slovnik", "pojmy.yaml")

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

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

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

```{r}
#| label: nastaveni-slovnik
#| echo: false
#| 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"]]
  )
}
```

::: {.box-question title="Hlavní otázka"}
Jaké důkazy potřebujeme, abychom mohli obhájit, že model je pro naši biologickou otázku dobrý?
:::

## Co už s daty umíme udělat?

Při analýze dat už umíme propojit několik kroků: data prohlédnout, zobrazit, položit biologickou otázku, odhadnout model, interpretovat jeho výstup a zkontrolovat jeho předpoklady.

::: {.box-question title="Sestavte pracovní postup"}
V jakém pořadí byste použili kroky **prohlédnout data**, **zobrazit vztah**, **formulovat model**, **interpretovat výstup** a **zkontrolovat `r render_glossary_term("residuum", display = "residua")`**? Než budete pokračovat, sestavte si vlastní pořadí.
:::

Jedna obhajitelná cesta je položit biologickou otázku, prohlédnout a zobrazit data, formulovat model, interpretovat jeho výstup a zkontrolovat residua. Jednotlivé kroky se mohou opakovat: `r render_glossary_term("diagnostika-modelu", display = "diagnostika")` nás například může vrátit k biologické otázce nebo k formulaci modelu. V této lekci známý postup použijeme na nový soubor dat a rozšíříme jej o `r render_glossary_term("porovnani-modelu", display = "porovnání několika modelů")`.

## Výsledky učení

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

- vypočítat `r render_glossary_term("r2", display = "koeficient determinace R²")` (*coefficient of determination*) z celkového součtu čtverců a `r render_glossary_term("soucet-ctvercu-residui", display = "součtu čtverců residuí")` a vysvětlit, co popisuje;
- vysvětlit, proč nejvyšší R² automaticky neznamená nejlepší model;
- rozlišit, jak R², adjustované R² (*adjusted R-squared*) a `r render_glossary_term("aic", display = "Akaikeho informační kritérium")` (*Akaike information criterion*, AIC) zacházejí se shodou a složitostí, a používat AIC jen k relativnímu `r render_glossary_term("porovnani-modelu", display = "porovnání modelů")`;
- sestavit malou sadu biologicky zdůvodněných kandidátních modelů se stejnou `r render_glossary_term("odezva", display = "odezvou")` a stejnými pozorováními;
- obhájit volbu modelu pomocí biologické otázky, složitosti, `r render_glossary_term("diagnostika-modelu", display = "diagnostiky")` a interpretovatelnosti a pojmenovat její omezení.

## Jeden řádek tabulky představuje jeden strom

Budeme pracovat s měřením tří druhů severoamerických jehličnanů. Každý řádek dat představuje jeden strom. Pro tuto lekci využijeme druh stromu, průměr kmene měřený ve výšce 1,3 m a výšku stromu. V minulé lekci výška stromu pomáhala odhadnout objem dřeva u jiných stromů; tady jde o jiné stromy a výška je naopak `r render_glossary_term("odezva", display = "odezvou")`, kterou chceme popsat.

Soubor [allometrie_stromu.csv](https://cuni-natur-biostatistics.github.io/L09/current/data/allometrie_stromu.csv) je součástí materiálů lekce. Pro následující kód jej uložte do složky `data` ve svém projektu v R.

```{r}
#| label: nacist-data-stromu
#| echo: true
#| eval: true

# Načteme výukovou tabulku stromů
data_stromy <-
  read.csv(file = "data/allometrie_stromu.csv")

# Druh stromu převedeme na kategoriální proměnnou
data_stromy$druh_stromu <-
  factor(x = data_stromy$druh_stromu)
```

Tabulka obsahuje `r nrow(data_stromy)` stromů a `r nlevels(data_stromy$druh_stromu)` druhy. Než začneme modelovat, zkontrolujeme několik řádků, strukturu `r render_glossary_term("promenna", display = "proměnných")`, chybějící hodnoty a zastoupení druhů.

```{r}
#| label: prohlednout-data-stromu
#| echo: true
#| eval: true

# Prohlédneme první řádky proměnných použitých v modelu
head(
  x = data_stromy[
    ,
    c(
      "druh_stromu",
      "prumer_kmene_cm",
      "vyska_stromu_m"
    )
  ],
  n = 6
)

# Ověříme datové typy a počet chybějících hodnot
str(object = data_stromy)
colSums(x = is.na(data_stromy))

# Zkontrolujeme počet stromů jednotlivých druhů
table(
  data_stromy$druh_stromu,
  useNA = "ifany"
)

# Zkontrolujeme rozsahy průměru a výšky
range(
  x = data_stromy$prumer_kmene_cm,
  na.rm = TRUE
)

range(
  x = data_stromy$vyska_stromu_m,
  na.rm = TRUE
)

# Porovnáme rozsahy průměru kmene mezi druhy
tapply(
  X = data_stromy$prumer_kmene_cm,
  INDEX = data_stromy$druh_stromu,
  FUN = range,
  na.rm = TRUE
)
```

Průměry kmenů sahají od `r format_value_cz(min(data_stromy$prumer_kmene_cm), digits = 1)` do `r format_value_cz(max(data_stromy$prumer_kmene_cm), digits = 1)` cm a výšky od `r format_value_cz(min(data_stromy$vyska_stromu_m), digits = 1)` do `r format_value_cz(max(data_stromy$vyska_stromu_m), digits = 1)` m. Rozsahy průměru kmene jednotlivých druhů se překrývají, takže druhy můžeme porovnávat alespoň v části společného rozsahu. Hodnoty jsou pro měřené stromy věrohodné; nevidíme záporné hodnoty ani zjevnou chybu jednotek.

::: {.box-data title="Doplňující: původ dat" collapse="true"}
Data pocházejí z objektu `allometry` v balíčku [`lgrdata`](https://cran.r-project.org/package=lgrdata), verze 0.1.1. Měření jsou připisována Johnu Marshallovi z University of Idaho a byla zveřejněna pod licencí CC0. Pokud je tato verze balíčku `lgrdata` nainstalovaná, původní tabulku lze načíst přímo z něj; není třeba stahovat ani rozbalovat archiv. Následující kód ověří verzi balíčku, načte původní tabulku, přejmenuje proměnné a uloží soubor použitý v této lekci.

```{r}
#| label: reprodukovat-data-stromu
#| echo: true
#| eval: false

# Ověříme verzi zdrojového balíčku
stopifnot(
  as.character(
    x = packageVersion(pkg = "lgrdata")
  ) == "0.1.1"
)

# Načteme původní objekt z balíčku
data(
  list = "allometry",
  package = "lgrdata"
)

data_allometry_source <- allometry

# Přeložíme zkratky druhů a vybereme výukové proměnné
vec_nazvy_druhu <-
  c(
    "PIMO" = "Pinus monticola",
    "PIPO" = "Pinus ponderosa",
    "PSME" = "Pseudotsuga menziesii"
  )

data_allometry_teaching <-
  data.frame(
    druh_stromu = unname(
      obj = vec_nazvy_druhu[
        as.character(x = data_allometry_source$species)
      ]
    ),
    prumer_kmene_cm = data_allometry_source$diameter,
    vyska_stromu_m = data_allometry_source$height,
    plocha_listu_m2 = data_allometry_source$leafarea,
    hmotnost_vetvi_kg = data_allometry_source$branchmass
  )

# Ověříme počet řádků a chybějící hodnoty
stopifnot(
  nrow(data_allometry_teaching) ==
    nrow(data_allometry_source),
  !anyNA(data_allometry_teaching)
)

# Uložíme stejnou tabulku, jakou načítá tato lekce
write.csv(
  x = data_allometry_teaching,
  file = "data/allometrie_stromu.csv",
  row.names = FALSE,
  fileEncoding = "UTF-8"
)
```
:::

## Jak se výška mění s průměrem kmene?

Vztah dvou `r render_glossary_term("numericka-promenna", display = "numerických proměnných")` zobrazíme pomocí `r render_glossary_term("bodovy-graf", display = "bodového grafu (anglicky scatter plot)")`. Druh stromu rozlišuje barva i tvar bodu, ale každý bod stále představuje jeden konkrétní strom.

```{r}
#| label: zobrazit-vztah-vysky-a-prumeru
#| echo: true
#| eval: true
#| fig-cap: "Výška stromu v závislosti na průměru kmene. Barva a tvar bodu rozlišují druhy stromů."
#| fig-alt: "Bodový graf výšky stromu proti průměru kmene u jehličnanů tří druhů. Výška s průměrem kmene roste; body druhu Pinus monticola leží v části grafu výše než body ostatních druhů."

# Barvy a tvary druhů budeme používat stejně ve všech grafech
vec_barvy_druhu <-
  c(
    "Pinus monticola" = "#4E9A8A",
    "Pinus ponderosa" = "#C2613A",
    "Pseudotsuga menziesii" = "#2E6CA4"
  )

vec_symboly_druhu <-
  c(
    "Pinus monticola" = 16,
    "Pinus ponderosa" = 17,
    "Pseudotsuga menziesii" = 15
  )

# Zobrazíme vztah průměru kmene a výšky stromu
plot(
  x = data_stromy$prumer_kmene_cm,
  y = data_stromy$vyska_stromu_m,
  col = vec_barvy_druhu[as.character(data_stromy$druh_stromu)],
  pch = vec_symboly_druhu[as.character(data_stromy$druh_stromu)],
  cex = 1.2,
  cex.axis = 1.2,
  cex.lab = 1.2,
  xlab = "Průměr kmene ve výšce 1,3 m (cm)",
  ylab = "Výška stromu (m)"
)

legend(
  x = "bottomright",
  legend = names(vec_barvy_druhu),
  col = vec_barvy_druhu,
  pch = vec_symboly_druhu,
  text.font = 3,
  bty = "n"
)
```

::: {.box-question title="Co vidíte v datech?"}
Popište nejsilnější společný vzorec a jednu vlastnost grafu, kvůli které by jediná přímka nemusela stačit.
:::

Výška se s rostoucím průměrem kmene většinou zvětšuje. Barvy druhů však v některých částech grafu leží v odlišných výškových pásmech, takže jedna společná přímka může přehlížet rozdíly mezi druhy.

## Jak vztah popíšeme jedním `r render_glossary_term("linearni-model", display = "lineárním modelem")`?

Jako první pokus použijeme lineární model, který výšku všech stromů popíše jednou společnou přímkou. `r render_glossary_term("odezva", display = "Odezvou")` je výška stromu a `r render_glossary_term("prediktor", display = "prediktorem")` je průměr kmene.

```{r}
#| label: odhadnout-spolecny-model
#| echo: true
#| eval: true
#| fig-alt: "Stejný bodový graf výšky proti průměru kmene s jednou fialovou přímkou společného modelu pro všechny tři druhy."

# Odhadneme jednu společnou přímku pro všechny stromy
mod_spolecny <-
  lm(
    formula = vyska_stromu_m ~ prumer_kmene_cm,
    data = data_stromy
  )

# Přidáme odhadnutou přímku do grafu
plot(
  x = data_stromy$prumer_kmene_cm,
  y = data_stromy$vyska_stromu_m,
  col = vec_barvy_druhu[as.character(data_stromy$druh_stromu)],
  pch = vec_symboly_druhu[as.character(data_stromy$druh_stromu)],
  cex = 1.2,
  cex.axis = 1.2,
  cex.lab = 1.2,
  xlab = "Průměr kmene ve výšce 1,3 m (cm)",
  ylab = "Výška stromu (m)"
)

abline(
  reg = mod_spolecny,
  col = "#5D2890",
  lwd = 3
)
```

Tento model záměrně ignoruje druh stromu a popisuje všechny stromy jedinou přímkou.

## Co už umíme přečíst ze `summary()`?

Výstup `summary()` obsahuje části, se kterými už umíme pracovat, i dva dosud nevysvětlené řádky. Prohlédněte si celý výpis a zkuste je rozlišit.

```{r}
#| label: zobrazit-souhrn-spolecneho-modelu
#| echo: true
#| eval: true

# Zobrazíme úplný souhrn odhadnutého modelu
summary(object = mod_spolecny)
```

::: {.box-question title="Známé a nové části výstupu"}
Které části výstupu už dokážete interpretovat? Které dva řádky zatím představují nové otázky?
:::

Už umíme číst volání modelu, souhrn `r render_glossary_term("residuum", display = "residuí")` a tabulku `r render_glossary_term("parametry-modelu", display = "koeficientů")` s jejich `r render_glossary_term("nejistota", display = "nejistotou")`. Novými otázkami zůstává význam řádků **Multiple R-squared** a **Adjusted R-squared**.

## Co model přehlíží?

Z předchozích lekcí už známe `r render_glossary_term("residuum", display = "residua")`: rozdíly mezi pozorovanými a `r render_glossary_term("odhadnuta-hodnota", display = "modelem odhadnutými hodnotami")`. Vhodný model by v nich neměl zanechávat snadno rozpoznatelný systematický vzorec. Zde se proto neptáme, co residuum znamená, ale zda společná přímka přehlíží rozdíl mezi druhy.

```{r}
#| label: zobrazit-residua-podle-druhu
#| echo: true
#| eval: true
#| fig-cap: "Residua společného modelu podle odhadnuté výšky. Barva a tvar bodu rozlišují druh stromu."
#| fig-alt: "Residua společného modelu proti odhadnuté výšce. Body druhu Pinus monticola leží převážně nad přerušovanou nulovou čarou."

# Zobrazíme residua a zachováme barvy a tvary druhů
plot(
  x = fitted(object = mod_spolecny),
  y = residuals(object = mod_spolecny),
  col = vec_barvy_druhu[as.character(data_stromy$druh_stromu)],
  pch = vec_symboly_druhu[as.character(data_stromy$druh_stromu)],
  cex = 1.2,
  cex.axis = 1.2,
  cex.lab = 1.2,
  xlab = "Odhadnutá výška stromu (m)",
  ylab = "Residuum (m)"
)

abline(
  h = 0,
  col = "#2E2E2E",
  lty = 2
)
```

```{r}
#| label: zkontrolovat-normalitu-spolecneho-modelu
#| echo: true
#| eval: true
#| fig-cap: "Q–Q graf residuí společného modelu."
#| fig-alt: "Q–Q graf residuí společného modelu: body leží v prostřední části blízko přímky, na obou koncích se od ní odchylují."

# Porovnáme residua s očekáváním normálního rozdělení
qqnorm(
  y = residuals(object = mod_spolecny),
  pch = 19,
  cex.axis = 1.2,
  cex.lab = 1.2,
  main = "Q–Q graf residuí",
  xlab = "Teoretické kvantily",
  ylab = "Kvantily residuí (m)"
)

qqline(
  y = residuals(object = mod_spolecny),
  col = "#2E2E2E",
  lwd = 2
)
```

::: {.box-question title="Detektivní práce s residui"}
Zůstává v residuích společného modelu rozdíl mezi druhy? Co by takový vzorec znamenal pro jednu společnou přímku?
:::

```{r}
#| label: spocitat-prumerna-residua-druhu
#| echo: false

vec_prumerna_residua_spolecny <-
  tibble::tibble(
    druh_stromu = as.character(data_stromy$druh_stromu),
    residuum = residuals(object = mod_spolecny)
  ) |>
  dplyr::group_by(druh_stromu) |>
  dplyr::summarise(
    prumerne_residuum = mean(residuum),
    .groups = "drop"
  ) |>
  tibble::deframe()
```

Residua druhu *Pinus monticola* leží převážně nad nulou; jejich průměr je `r format_value_cz(vec_prumerna_residua_spolecny[["Pinus monticola"]], digits = 1, signed = TRUE)` m. Společná přímka tedy výšku tohoto druhu soustavně podhodnocuje a výšku ostatních dvou druhů spíše nadhodnocuje (průměrná residua `r format_value_cz(vec_prumerna_residua_spolecny[["Pinus ponderosa"]], digits = 1, signed = TRUE)` a `r format_value_cz(vec_prumerna_residua_spolecny[["Pseudotsuga menziesii"]], digits = 1, signed = TRUE)` m). Takový vzorec v residuích ukazuje, že model přehlíží rozdíl mezi druhy.

## Je jeden model jedinou rozumnou možností?

Společná přímka odpovídá jedné biologické představě: stejně silné stromy mají podle modelu stejnou `r render_glossary_term("odhadnuta-hodnota", display = "odhadnutou výšku")` bez ohledu na druh. Pozorovaný graf a `r render_glossary_term("residuum", display = "residua")` ale připouštějí také rovnoběžné přímky jednotlivých druhů nebo přímky s rozdílnými `r render_glossary_term("sklon", display = "sklony")`.

## Tři modely vyjadřují tři biologické představy

Každý kandidátní model vyjadřuje jinou představu o vztahu průměru kmene, druhu stromu a jeho výšky.

| Kandidát | Vzorec v R | Biologická představa |
|---|---|---|
| Společná přímka | `vyska_stromu_m ~ prumer_kmene_cm` | Všechny druhy sdílejí stejný `r render_glossary_term("intercept", display = "intercept")` i `r render_glossary_term("sklon", display = "sklon")`. |
| Posun mezi druhy | `vyska_stromu_m ~ prumer_kmene_cm + druh_stromu` | Druhy mohou mít odlišnou výšku při stejném průměru, ale sdílejí sklon. |
| Různé sklony | `vyska_stromu_m ~ prumer_kmene_cm * druh_stromu` | Vztah mezi průměrem a výškou se může mezi druhy lišit. |

```{r}
#| label: odhadnout-kandidatni-modely
#| echo: true
#| eval: true

# Model 1: jedna společná přímka
mod_spolecny <-
  lm(
    formula = vyska_stromu_m ~ prumer_kmene_cm,
    data = data_stromy
  )

# Model 2: rovnoběžné přímky druhů
mod_aditivni <-
  lm(
    formula = vyska_stromu_m ~
      prumer_kmene_cm +
      druh_stromu,
    data = data_stromy
  )

# Model 3: druhově specifické sklony
mod_interakce <-
  lm(
    formula = vyska_stromu_m ~
      prumer_kmene_cm *
      druh_stromu,
    data = data_stromy
  )
```

## Který model prochází nejblíže pozorovaným stromům?

```{r}
#| label: pripravit-predikce-kandidatu
#| echo: false
#| output: false
#| message: false
#| warning: false
#| error: false

data_mrizka <-
  expand.grid(
    prumer_kmene_cm = seq(
      from = min(data_stromy$prumer_kmene_cm),
      to = max(data_stromy$prumer_kmene_cm),
      length.out = 120
    ),
    druh_stromu = levels(data_stromy$druh_stromu)
  )

data_predikce_spolecny <-
  data_mrizka |>
  dplyr::mutate(
    vyska_stromu_m = predict(
      object = mod_spolecny,
      newdata = data_mrizka
    ),
    model = "Společná přímka"
  )

data_predikce_aditivni <-
  data_mrizka |>
  dplyr::mutate(
    vyska_stromu_m = predict(
      object = mod_aditivni,
      newdata = data_mrizka
    ),
    model = "Posun mezi druhy"
  )

data_predikce_interakce <-
  data_mrizka |>
  dplyr::mutate(
    vyska_stromu_m = predict(
      object = mod_interakce,
      newdata = data_mrizka
    ),
    model = "Různé sklony"
  )

data_predikce <-
  dplyr::bind_rows(
    data_predikce_spolecny,
    data_predikce_aditivni,
    data_predikce_interakce
  ) |>
  dplyr::mutate(
    model = factor(
      model,
      levels = c(
        "Společná přímka",
        "Posun mezi druhy",
        "Různé sklony"
      )
    )
  )

data_predikce_spolecny_graf <-
  data_predikce |>
  dplyr::filter(
    model == "Společná přímka",
    druh_stromu == levels(data_stromy$druh_stromu)[1]
  )
```

```{r}
#| label: zobrazit-kandidatni-modely
#| echo: false
#| message: false
#| warning: false
#| fig-cap: "Stejná pozorování a tři kandidátní modely. Panely se liší pouze strukturou odhadnutých přímek."
#| fig-alt: "Tři panely se stejnými stromy: vlevo jedna fialová společná přímka, uprostřed tři rovnoběžné přímky druhů, vpravo tři přímky druhů s různými sklony."
#| fig-width: 9
#| fig-height: 4.8

ggplot2::ggplot(
  data = data_stromy,
  mapping = ggplot2::aes(
    x = prumer_kmene_cm,
    y = vyska_stromu_m,
    color = druh_stromu,
    shape = druh_stromu
  )
) +
  ggplot2::geom_point(size = 2.3, alpha = 0.8) +
  ggplot2::geom_line(
    data = data_predikce_spolecny_graf,
    mapping = ggplot2::aes(
      x = prumer_kmene_cm,
      y = vyska_stromu_m
    ),
    inherit.aes = FALSE,
    color = biostat_cols[["indigo_velvet"]],
    linewidth = 1.2
  ) +
  ggplot2::geom_line(
    data = data_predikce |>
      dplyr::filter(model != "Společná přímka"),
    mapping = ggplot2::aes(
      x = prumer_kmene_cm,
      y = vyska_stromu_m,
      color = druh_stromu,
      group = druh_stromu
    ),
    inherit.aes = FALSE,
    linewidth = 1.1
  ) +
  ggplot2::facet_wrap(
    facets = ggplot2::vars(model),
    nrow = 1
  ) +
  ggplot2::scale_color_manual(values = vec_barvy_druhu) +
  ggplot2::scale_shape_manual(values = vec_symboly_druhu) +
  ggplot2::labs(
    x = "Průměr kmene ve výšce 1,3 m (cm)",
    y = "Výška stromu (m)",
    color = NULL,
    shape = NULL
  ) +
  theme_biostat(base_size = 13) +
  ggplot2::theme(
    legend.position = "bottom",
    legend.text = ggplot2::element_text(face = "italic")
  )
```

::: {.box-question title="Porovnejte vzdálenosti bodů od přímek"}
Ve kterém panelu zůstávají mezi pozorovanými body a odhadnutými přímkami nejkratší svislé vzdálenosti?
:::

Podle obrázku procházejí pozorovaným bodům nejblíže přímky v panelu **Různé sklony**. Rozdíly mezi body a přímkami však chceme vyjádřit číslem, ne pouze odhadovat okem.

## Jak poznáme, který model je lepší?

Při `r render_glossary_term("porovnani-modelu", display = "porovnání modelů")` potřebujeme rozlišit několik druhů důkazů:

- jak těsně model odpovídá pozorovaným datům;
- kolik `r render_glossary_term("parametry-modelu", display = "parametrů")` potřebuje;
- zda v `r render_glossary_term("residuum", display = "residuích")` nezůstává systematický vzorec;
- zda odpovídá biologické otázce;
- zda jeho výsledek dokážeme smysluplně interpretovat.

Žádné jediné číslo neodpoví na všechny tyto otázky.

## Co už víme o `r render_glossary_term("residuum", display = "residuích")` a jejich čtvercích?

Když jsme poprvé hledali přímku, porovnávali jsme kandidáty pomocí `r render_glossary_term("soucet-ctvercu-residui", display = "součtu čtverců residuí")`. V další lekci jsme stejnou veličinu označili anglickou zkratkou **RSS** (*residual sum of squares*) a použili ji při výpočtu `r render_glossary_term("rezidualni-smerodatna-odchylka", display = "reziduální směrodatné odchylky")`. Při porovnání jednoho společného `r render_glossary_term("prumer", display = "průměru")` s několika skupinovými průměry jsme také viděli, že složitější model může RSS zmenšit.

RSS proto není nový pojem. V této lekci jej spojíme s druhým součtem čtverců, který popisuje rozdíly mezi pozorovanými výškami ještě před použitím `r render_glossary_term("prediktor", display = "prediktorů")`.

::: {.box-question title="Vzpomeňte si na význam RSS"}
Co se stane se součtem čtverců residuí, když odhadnutá přímka prochází blíže pozorovaným bodům?
:::

Čím jsou residua kratší, tím menší jsou jejich čtverce i jejich součet. **Menší RSS proto znamená těsnější shodu modelu s pozorovanými daty.**

{{< pagebreak >}}

## Je přímka lepší než předpověď `r render_glossary_term("prumer", display = "průměrem")`?

Než posoudíme odhadnutou přímku, potřebujeme jednoduchý výchozí způsob odhadu. Pokud bychom o jednotlivých stromech nevěděli nic dalšího, mohli bychom každému přiřadit průměrnou výšku. Tento průměrový model nám umožní popsat, jak moc se pozorované výšky liší ještě před použitím průměru kmene.

Následující výpočet používá šest skutečných stromů z datové tabulky. Malý výběr umožňuje vidět každý krok výpočtu. Hodnoty kandidátních modelů v dalších tabulkách vycházejí ze všech `r nrow(data_stromy)` stromů.

```{r}
#| label: pripravit-priklad-r2
#| echo: false
#| output: false

vec_radky_priklad_r2 <- c(1, 11, 22, 32, 43, 58)

data_priklad_r2 <-
  data_stromy[
    vec_radky_priklad_r2,
    c(
      "druh_stromu",
      "prumer_kmene_cm",
      "vyska_stromu_m"
    )
  ]

data_priklad_r2$strom <- LETTERS[seq_len(nrow(data_priklad_r2))]
prumerna_vyska_priklad <- mean(data_priklad_r2$vyska_stromu_m)

mod_priklad_r2 <-
  lm(
    formula = vyska_stromu_m ~ prumer_kmene_cm,
    data = data_priklad_r2
  )

data_priklad_r2$odchylka_od_prumeru <-
  data_priklad_r2$vyska_stromu_m -
  prumerna_vyska_priklad

data_priklad_r2$ctverec_odchylky <-
  data_priklad_r2$odchylka_od_prumeru^2

data_priklad_r2$odhadnuta_vyska <-
  fitted(object = mod_priklad_r2)

data_priklad_r2$residuum <-
  residuals(object = mod_priklad_r2)

data_priklad_r2$absolutni_residuum <-
  abs(data_priklad_r2$residuum)

data_priklad_r2$ctverec_residua <-
  data_priklad_r2$residuum^2

ss_total_priklad <-
  sum(data_priklad_r2$ctverec_odchylky)

ss_res_priklad <-
  sum(data_priklad_r2$ctverec_residua)

strom_max_total <-
  data_priklad_r2$strom[
    which.max(data_priklad_r2$ctverec_odchylky)
  ]

index_max_residuum <-
  which.max(data_priklad_r2$residuum)

index_max_absolutni_residuum <-
  which.max(data_priklad_r2$absolutni_residuum)

strom_max_residuum <-
  data_priklad_r2$strom[index_max_residuum]

strom_max_absolutni_residuum <-
  data_priklad_r2$strom[index_max_absolutni_residuum]

# Expanded squares, then the collapsed squared values, for the equations
text_ss_total_priklad <-
  paste0(
    "(",
    format_value_cz(
      x = data_priklad_r2$odchylka_od_prumeru,
      digits = 2,
      math = TRUE
    ),
    ")^2",
    collapse = " + "
  )

text_ss_total_ctverce_priklad <-
  paste(
    format_value_cz(
      x = data_priklad_r2$ctverec_odchylky,
      digits = 2,
      math = TRUE
    ),
    collapse = " + "
  )

text_ss_res_priklad <-
  paste0(
    "(",
    format_value_cz(
      x = data_priklad_r2$residuum,
      digits = 2,
      math = TRUE
    ),
    ")^2",
    collapse = " + "
  )

text_ss_res_ctverce_priklad <-
  paste(
    format_value_cz(
      x = data_priklad_r2$ctverec_residua,
      digits = 2,
      math = TRUE
    ),
    collapse = " + "
  )
```

```{r}
#| label: zobrazit-vyber-pro-r2
#| echo: false
#| results: asis

knitr::kable(
  x = data_priklad_r2[
    ,
    c(
      "strom",
      "prumer_kmene_cm",
      "vyska_stromu_m"
    )
  ],
  digits = 2,
  row.names = FALSE,
  col.names = c(
    "Strom",
    "Průměr kmene (cm)",
    "Výška (m)"
  ),
  caption = "Šest stromů použitých k porovnání průměrového modelu a odhadnuté přímky."
)
```

## Jak velké jsou rozdíly výšek před použitím `r render_glossary_term("prediktor", display = "prediktorů")`?

`r render_glossary_term("prumer", display = "Průměrná výška")` šesti stromů je `r format_value_cz(prumerna_vyska_priklad, digits = 2)` m. Svislé úsečky ukazují odchylku každého `r render_glossary_term("pozorovani", display = "pozorování")` od tohoto průměru.

```{r}
#| label: zobrazit-celkovou-variabilitu
#| echo: false
#| fig-cap: "Odchylky šesti pozorovaných výšek od jejich průměru (přerušovaná čára). Délky úseček jsou před umocněním kladné nebo záporné odchylky."
#| fig-alt: "Šest stromů A až F s pozorovanými výškami a přerušovanou vodorovnou čarou průměrné výšky; svislé úsečky spojují každý bod s průměrem, nejdelší má strom E pod průměrem a strom F nad ním."

plot(
  x = seq_len(nrow(data_priklad_r2)),
  y = data_priklad_r2$vyska_stromu_m,
  pch = 19,
  cex = 1.3,
  cex.axis = 1.2,
  cex.lab = 1.2,
  col = biostat_cols[["grey_olive"]],
  xaxt = "n",
  xlab = "Strom",
  ylab = "Výška stromu (m)",
  ylim = range(
    data_priklad_r2$vyska_stromu_m,
    prumerna_vyska_priklad
  )
)

axis(
  side = 1,
  at = seq_len(nrow(data_priklad_r2)),
  labels = data_priklad_r2$strom,
  cex.axis = 1.2
)

abline(
  h = prumerna_vyska_priklad,
  col = biostat_cols[["graphite"]],
  lty = 2,
  lwd = 2
)

segments(
  x0 = seq_len(nrow(data_priklad_r2)),
  y0 = prumerna_vyska_priklad,
  x1 = seq_len(nrow(data_priklad_r2)),
  y1 = data_priklad_r2$vyska_stromu_m,
  col = biostat_cols[["graphite"]],
  lwd = 3
)
```

Kladné a záporné odchylky by se při prostém součtu vyrušily. Proto je umocníme na druhou a sečteme. Vznikne **celkový součet čtverců** (*total sum of squares*), tedy číselné vyjádření celkové `r render_glossary_term("variabilita", display = "variability")` pozorovaných výšek kolem jejich průměru.

<!-- TODO: add slug celkovy-soucet-ctvercu to pojmy.yaml -->

```{r}
#| label: tabulka-celkove-variability
#| echo: false
#| results: asis

knitr::kable(
  x = data_priklad_r2[
    ,
    c(
      "strom",
      "vyska_stromu_m",
      "odchylka_od_prumeru",
      "ctverec_odchylky"
    )
  ],
  digits = 2,
  row.names = FALSE,
  col.names = c(
    "Strom",
    "Výška",
    "Odchylka od průměru",
    "Čtverec odchylky"
  ),
  caption = "Od pozorovaných výšek k celkovému součtu čtverců."
)
```

Celkový součet čtverců sčítá čtverce odchylek všech šesti stromů:

$$
\text{celkový součet čtverců} =
\text{čtverec odchylky stromu A} + \cdots +
\text{čtverec odchylky stromu F}.
$$

Do vztahu dosadíme odchylky z tabulky, umocníme je a sečteme:

$$
\begin{align}
\text{celkový součet čtverců}
&= `r text_ss_total_priklad` \\
&= `r text_ss_total_ctverce_priklad` \\
&= `r format_value_cz(ss_total_priklad, digits = 2, math = TRUE)`.
\end{align}
$$

Obecně označíme pozorovanou hodnotu $y_i$, průměr všech pozorovaných hodnot $\bar{y}$ a počet pozorování $n$:

$$
SS_{\mathrm{total}} =
\sum_{i=1}^{n}
\left(y_i - \bar{y}\right)^2.
$$

Symbol $SS_{\mathrm{total}}$ označuje celkový součet čtverců. Každý člen $y_i - \bar{y}$ je odchylka jedné pozorované výšky od průměrné výšky. Hodnota `r format_value_cz(ss_total_priklad, digits = 2)` m² popisuje celkovou čtvercovou odchylku šesti výšek od jejich průměru.

::: {.box-question title="Který strom přispívá nejvíce?"}
Najděte v grafu a tabulce strom s největším příspěvkem k celkové variabilitě. Proč o velikosti příspěvku rozhoduje délka úsečky **i její umocnění na druhou**?
:::

Největší příspěvek má strom `r strom_max_total`, protože má největší absolutní odchylku od průměru. Umocnění zachová příspěvky kladných i záporných odchylek a větší odchylky zvýrazní.

## Kolik rozdílů zůstalo po odhadnutí přímky?

Místo `r render_glossary_term("prumer", display = "průměru")` teď použijeme odhadnutou přímku. Každá úsečka vede od `r render_glossary_term("odhadnuta-hodnota", display = "modelem odhadnuté výšky")` k pozorované výšce stejného stromu. Její orientovaná délka je známé `r render_glossary_term("residuum", display = "residuum")`:

$$
\text{residuum} =
\text{pozorovaná výška} -
\text{odhadnutá výška}.
$$

```{r}
#| label: zobrazit-residualni-variabilitu
#| echo: false
#| fig-cap: "Residua stejné šestice stromů. Každá oranžová úsečka spojuje odhadnutou výšku na fialové přímce s pozorovanou výškou."
#| fig-alt: "Šest stromů A až F podle průměru kmene a výšky s fialovou odhadnutou přímkou; oranžové svislé úsečky residuí spojují každý bod s přímkou."

plot(
  x = data_priklad_r2$prumer_kmene_cm,
  y = data_priklad_r2$vyska_stromu_m,
  pch = 19,
  cex = 1.3,
  cex.axis = 1.2,
  cex.lab = 1.2,
  col = biostat_cols[["grey_olive"]],
  xlab = "Průměr kmene (cm)",
  ylab = "Výška stromu (m)",
  xlim = range(data_priklad_r2$prumer_kmene_cm) + c(-6, 3)
)

abline(
  reg = mod_priklad_r2,
  col = biostat_cols[["indigo_velvet"]],
  lwd = 2
)

segments(
  x0 = data_priklad_r2$prumer_kmene_cm,
  y0 = data_priklad_r2$odhadnuta_vyska,
  x1 = data_priklad_r2$prumer_kmene_cm,
  y1 = data_priklad_r2$vyska_stromu_m,
  col = biostat_cols[["orange"]],
  lwd = 3
)

# Letters sit to the left of each point so they do not cover it
text(
  x = data_priklad_r2$prumer_kmene_cm,
  y = data_priklad_r2$vyska_stromu_m,
  labels = data_priklad_r2$strom,
  pos = 2,
  offset = 0.8,
  font = 2,
  cex = 1.1
)
```

Residuum nese znaménko: bod nad přímkou má kladné residuum, bod pod přímkou záporné. Absolutní hodnota residua udává délku svislé úsečky mezi bodem a přímkou. Čtverec residua znaménko odstraní a určuje příspěvek daného stromu k součtu čtverců residuí.

<!-- TODO: add slug absolutni-hodnota to pojmy.yaml -->

Čtverce residuí sečteme. Dostaneme již známý `r render_glossary_term("soucet-ctvercu-residui", display = "součet čtverců residuí")` (*residual sum of squares*, RSS): kolik čtvercové odchylky zůstalo po odhadnutí modelu.

```{r}
#| label: tabulka-residualni-variability
#| echo: false
#| results: asis

knitr::kable(
  x = data_priklad_r2[
    ,
    c(
      "strom",
      "vyska_stromu_m",
      "odhadnuta_vyska",
      "residuum",
      "absolutni_residuum",
      "ctverec_residua"
    )
  ],
  digits = 2,
  row.names = FALSE,
  col.names = c(
    "Strom",
    "Pozorovaná výška",
    "Odhadnutá výška",
    "Residuum",
    "Absolutní hodnota",
    "Čtverec residua"
  ),
  caption = "Od residuí k součtu čtverců residuí (RSS)."
)
```

Součet čtverců residuí sčítá čtverce residuí všech šesti stromů:

$$
\text{součet čtverců residuí} =
\text{čtverec residua stromu A} + \cdots +
\text{čtverec residua stromu F}.
$$

Do vztahu dosadíme residua z tabulky, umocníme je a sečteme:

$$
\begin{align}
\text{součet čtverců residuí}
&= `r text_ss_res_priklad` \\
&= `r text_ss_res_ctverce_priklad` \\
&= `r format_value_cz(ss_res_priklad, digits = 2, math = TRUE)`.
\end{align}
$$

Index $i$ označuje pořadí `r render_glossary_term("pozorovani", display = "pozorování")`. Residuum $i$-tého pozorování označíme jako $e_i$:

$$
SS_{\mathrm{res}} =
\sum_{i=1}^{n} e_i^2.
$$

Symbol $SS_{\mathrm{res}}$ označuje RSS. Jeho hodnota `r format_value_cz(ss_res_priklad, digits = 2)` m² je menší než celkový součet čtverců, protože přímka sleduje část rozdílů mezi výškami lépe než samotný průměr. **Celkový součet čtverců a RSS jsou dvě stavební veličiny, které použijeme společně v následujícím výpočtu.**

::: {.box-question title="Residuum, jeho velikost a čtverec"}
Který strom má číselně nejvyšší residuum? Který strom má residuum s největší absolutní hodnotou, a proto také největší čtverec residua? Jsou to stejné stromy?
:::

Číselně nejvyšší residuum má strom `r strom_max_residuum`: `r format_value_cz(data_priklad_r2$residuum[index_max_residuum], digits = 2)` m. Největší absolutní hodnotu residua má strom `r strom_max_absolutni_residuum`: jeho residuum je `r format_value_cz(data_priklad_r2$residuum[index_max_absolutni_residuum], digits = 2)` m, jeho absolutní hodnota je `r format_value_cz(data_priklad_r2$absolutni_residuum[index_max_absolutni_residuum], digits = 2)` m a čtverec residua je `r format_value_cz(data_priklad_r2$ctverec_residua[index_max_absolutni_residuum], digits = 2)` m². `r if (index_max_residuum == index_max_absolutni_residuum) "Jde o stejný strom." else "Nejde o stejné stromy."` Nejde o chybu: umocnění odstraní znaménko, takže příspěvek k RSS určuje velikost residua bez ohledu na to, zda je kladné, nebo záporné.

## Jakou část `r render_glossary_term("variabilita", display = "variability")` zachytil model?

Máme dvě veličiny: celkový součet čtverců před použitím `r render_glossary_term("prediktor", display = "prediktoru")` a `r render_glossary_term("soucet-ctvercu-residui", display = "součet čtverců residuí (RSS)")` po odhadnutí přímky. Jejich poměr říká, jaká část původních rozdílů zůstala nevysvětlena. Jeho doplněk do jedné proto vyjadřuje **zachycený podíl variability**.

```{r}
#| label: vypocitat-r2-prikladu
#| echo: false
#| output: false

r2_priklad <-
  1 -
  ss_res_priklad /
  ss_total_priklad
```

Zachycený podíl vyjádříme pomocí dvou veličin, které už známe:

$$
\text{zachycený podíl variability} =
1 -
\frac{\text{součet čtverců residuí}}
{\text{celkový součet čtverců}}.
$$

Do vztahu dosadíme všech šest residuí a všech šest odchylek od průměru:

$$
\text{zachycený podíl variability} =
1 -
\frac{`r text_ss_res_priklad`}
{`r text_ss_total_priklad`}.
$$

Po sečtení čtverců dostaneme:

$$
\text{zachycený podíl variability} =
1 -
\frac{`r format_value_cz(ss_res_priklad, digits = 2, math = TRUE)`}
{`r format_value_cz(ss_total_priklad, digits = 2, math = TRUE)`}
=
`r format_value_cz(r2_priklad, digits = 3, math = TRUE)`.
$$

Tento podíl nazýváme `r render_glossary_term("r2", display = "koeficient determinace R²")` (**coefficient of determination, R-squared**). Obecně jej zapíšeme:

$$
R^2 =
1 -
\frac{SS_{\mathrm{res}}}
{SS_{\mathrm{total}}}.
$$

Zkratka $SS$ znamená „součet čtverců“. $SS_{\mathrm{res}}$ je součet čtverců residuí od odhadnutého modelu. $SS_{\mathrm{total}}$ je součet čtverců odchylek pozorované `r render_glossary_term("odezva", display = "odezvy")` od jejího `r render_glossary_term("prumer", display = "průměru")`.

```{r}
#| label: pripravit-cviceni-r2
#| echo: false
#| output: false

ss_total_cviceni <- 500
ss_res_cviceni <- 125

r2_cviceni <-
  1 -
  ss_res_cviceni /
  ss_total_cviceni
```

::: {.box-question title="Vypočítejte a interpretujte"}
Kdyby $SS_{\mathrm{total}} = `r format_value_cz(ss_total_cviceni, digits = 0, math = TRUE)`$ a $SS_{\mathrm{res}} = `r format_value_cz(ss_res_cviceni, digits = 0, math = TRUE)`$, jaké by bylo R²? Napište také jednu větu o tom, co by takový výsledek říkal o odezvě.
:::

R² by bylo $1 - `r format_value_cz(ss_res_cviceni, digits = 0, math = TRUE)`/`r format_value_cz(ss_total_cviceni, digits = 0, math = TRUE)` = `r format_value_cz(r2_cviceni, digits = 2, math = TRUE)`$. Model by tedy v datech zachytil `r format_value_cz(100 * r2_cviceni, digits = 0)` % variability odezvy kolem jejího průměru.

## Dostaneme stejné `r render_glossary_term("r2", display = "R²")` také ze `summary()`?

V malém příkladu jsme mohli rozepsat každý čtverec. Pro všech `r nrow(data_stromy)` stromů použijeme stejný výpočet, ale dosadíme už jen oba součty čtverců.

```{r}
#| label: overit-r2-pro-vsechny-stromy
#| echo: false
#| output: false

prumerna_vyska_vsechny <-
  mean(data_stromy$vyska_stromu_m)

ss_total_vsechny <-
  sum(
    (
      data_stromy$vyska_stromu_m -
      prumerna_vyska_vsechny
    )^2
  )

ss_res_spolecny <-
  deviance(object = mod_spolecny)

r2_spolecny_vypocet <-
  1 -
  ss_res_spolecny /
  ss_total_vsechny

r2_spolecny_summary <-
  summary(object = mod_spolecny)$r.squared

stopifnot(
  isTRUE(
    all.equal(
      target = r2_spolecny_vypocet,
      current = r2_spolecny_summary
    )
  )
)
```

$$
R^2 =
1 -
\frac{`r format_value_cz(ss_res_spolecny, digits = 2, math = TRUE)`}
{`r format_value_cz(ss_total_vsechny, digits = 2, math = TRUE)`}
=
`r format_value_cz(r2_spolecny_vypocet, digits = 3, math = TRUE)`.
$$

Stejnou hodnotu vrátí také `summary(mod_spolecny)$r.squared`: `r format_value_cz(r2_spolecny_summary, digits = 3)`. Výpis tedy nepřináší jinou definici; pouze za nás provede tentýž výpočet.

## Co nám `r render_glossary_term("r2", display = "R²")` říká — a co ne?

Pro všech `r nrow(data_stromy)` stromů je R² společné přímky `r format_value_cz(r2_spolecny_vypocet, digits = 3)`. Model tedy ve **stejných datech, na kterých byl odhadnut**, zachycuje přibližně `r format_value_cz(100 * r2_spolecny_vypocet, digits = 1)` % pozorované `r render_glossary_term("variabilita", display = "variability")` výšky kolem jejího `r render_glossary_term("prumer", display = "průměru")`.

R² však neříká:

- že model správně předpověděl `r format_value_cz(100 * r2_spolecny_vypocet, digits = 1)` % jednotlivých stromů;
- že model splňuje předpoklady;
- že vztah je příčinný;
- že model bude stejně dobrý pro nové stromy.

::: {.box-question title="Která interpretace je obhajitelná?"}
Vyberte a zdůvodněte jednu možnost:

- **A:** R² vyjadřuje podíl správně předpovězených stromů.
- **B:** R² vyjadřuje podíl zachycené variability `r render_glossary_term("odezva", display = "odezvy")` v těchto datech.
- **C:** R² vyjadřuje pravděpodobnost, že je model pravdivý.
:::

Obhajitelná je možnost **B**: R² vyjadřuje podíl zachycené variability odezvy v datech použitých k odhadu modelu.

## Proč `r render_glossary_term("r2", display = "R²")` po přidání členu neklesne?

Stejný výpočet R² můžeme použít pro všechny tři kandidáty na všech `r nrow(data_stromy)` stromech. Kandidátní modely jsou vnořené: složitější model vždy obsahuje všechny možnosti jednoduššího modelu a přidává další `r render_glossary_term("parametry-modelu", display = "koeficienty")`. Při odhadu metodou nejmenších čtverců tedy může složitější model zopakovat jednodušší řešení, nebo `r render_glossary_term("soucet-ctvercu-residui", display = "součet čtverců residuí (RSS)")` ještě zmenšit.

```{r}
#| label: pripravit-porovnani-modelu
#| echo: false
#| output: false

list_modely <-
  list(
    mod_spolecny,
    mod_aditivni,
    mod_interakce
  )

vec_nazvy_modelu <-
  c(
    "Společná přímka",
    "Posun mezi druhy",
    "Různé sklony"
  )

data_porovnani_modelu <-
  data.frame(
    model = vec_nazvy_modelu,
    n = purrr::map_int(
      .x = list_modely,
      .f = nobs
    ),
    pocet_koeficientu = purrr::map_int(
      .x = list_modely,
      .f = function(model) {
        length(coef(object = model))
      }
    ),
    ss_total = rep(
      x = ss_total_vsechny,
      times = length(list_modely)
    ),
    ss_res = purrr::map_dbl(
      .x = list_modely,
      .f = deviance
    ),
    r2 = purrr::map_dbl(
      .x = list_modely,
      .f = function(model) {
        summary(object = model)$r.squared
      }
    )
  )
```

```{r}
#| label: zobrazit-r2-kandidatu
#| echo: false
#| results: asis

knitr::kable(
  x = data_porovnani_modelu[
    ,
    c(
      "model",
      "pocet_koeficientu",
      "ss_res",
      "r2"
    )
  ],
  digits = c(0, 0, 1, 3),
  col.names = c(
    "Kandidátní model",
    "Počet koeficientů",
    "Součet čtverců residuí (RSS)",
    "R²"
  ),
  caption = "Při přidávání koeficientů součet čtverců residuí (RSS) klesá a R² roste."
)
```

::: {.box-warning title="Vyšší R² není automatické vítězství"}
Nejpružnější kandidát získává nejvyšší R² už proto, že má více možností přiblížit se stejným pozorováním. Samotné R² neumí rozhodnout, zda malé zlepšení shody stojí za další koeficienty.
:::

::: {.box-question title="Co udělá zbytečný člen?"}
Představte si, že do modelu přidáme `r render_glossary_term("prediktor", display = "prediktor")`, který téměř nic nepřinese. Může R² klesnout? Dokazuje jeho případné zvýšení, že je nový model užitečnější?
:::

R² nemůže klesnout, protože rozšířený model může vždy zopakovat odhady jednoduššího modelu. Případné zvýšení však samo nedokazuje, že jsou nové koeficienty užitečné.

## Které `r render_glossary_term("parametry-modelu", display = "koeficienty")` jednotlivé modely odhadují?

Počet koeficientů není abstraktní číslo. Každý z nich odpovídá konkrétní části modelu: `r render_glossary_term("intercept", display = "interceptu")`, `r render_glossary_term("sklon", display = "sklonu")` nebo rozdílu mezi druhy.

::: {.box-question title="Přiřaďte koeficienty k modelům"}
Než se podíváte na následující seznam, rozhodněte, který model potřebuje pouze intercept a jeden sklon, který přidává rozdíly druhů a který navíc přidává rozdíly ve sklonech.
:::

```{r}
#| label: zobrazit-koeficienty-kandidatu
#| echo: false
#| results: asis

vec_radky_koeficientu <-
  purrr::map2_chr(
    .x = vec_nazvy_modelu,
    .y = list_modely,
    .f = function(nazev_modelu, model) {
      nazvy_koeficientu <-
        names(coef(object = model))

      n_koeficientu <- length(nazvy_koeficientu)

      # Czech plural: 2-4 koeficienty, 5 and more koeficientů
      slovo_koeficient <-
        dplyr::if_else(
          condition = n_koeficientu %in% 2:4,
          true = "koeficienty",
          false = "koeficientů"
        )

      paste0(
        "- **",
        nazev_modelu,
        " — ",
        n_koeficientu,
        " ",
        slovo_koeficient,
        ":**\n",
        paste0(
          "    - `",
          nazvy_koeficientu,
          "`",
          collapse = "\n"
        )
      )
    }
  )

cat(
  paste(
    vec_radky_koeficientu,
    collapse = "\n"
  ),
  "\n"
)
```

Společná přímka má intercept a jeden sklon. Aditivní model přidává dva rozdíly druhů vůči referenčnímu druhu. `r render_glossary_term("interakce-predikatoru", display = "Interakční model")` navíc přidává dva rozdíly ve sklonech. Právě tyto konkrétní možnosti navíc mohou zlepšit shodu s daty.

## Stačí zlepšení `r render_glossary_term("r2", display = "R²")` na další `r render_glossary_term("parametry-modelu", display = "koeficienty")`?

<!-- TODO: add slug adjustovane-r2 to pojmy.yaml -->

```{r}
#| label: pripravit-adjustovane-r2
#| echo: false
#| output: false

data_porovnani_modelu$adjustovane_r2 <-
  purrr::map_dbl(
    .x = list_modely,
    .f = function(model) {
      summary(object = model)$adj.r.squared
    }
  )

r2_aditivni <-
  summary(object = mod_aditivni)$r.squared

podil_nezachyceny_aditivni <-
  1 - r2_aditivni

faktor_upravy_aditivni <-
  (nobs(object = mod_aditivni) - 1) /
  (
    nobs(object = mod_aditivni) -
    length(coef(object = mod_aditivni))
  )
```

**Adjustované R² (adjusted R-squared)** se ptá, zda je nárůst obyčejného R² dostatečně velký vzhledem ke dvěma dalším údajům:

- $n$ je počet `r render_glossary_term("pozorovani", display = "pozorování")` použitých k odhadu modelu;
- $k$ je počet odhadovaných regresních koeficientů **včetně `r render_glossary_term("intercept", display = "interceptu")`**.

V našich datech mají všechny modely stejné $n = `r nrow(data_stromy)`$, ale postupně odhadují `r paste(data_porovnani_modelu$pocet_koeficientu, collapse = ", ")` koeficientů. Pro aditivní model známe:

- obyčejné R² = `r format_value_cz(summary(mod_aditivni)$r.squared, digits = 3)`;
- počet pozorování $n = `r nobs(mod_aditivni)`$;
- počet regresních koeficientů $k = `r length(coef(mod_aditivni))`$.

Úprava vychází z podílu `r render_glossary_term("variabilita", display = "variability")`, který obyčejné R² nezachytilo. Tento podíl zvětší podle počtu pozorování a koeficientů:

$$
\text{adjustované R²} =
1 -
\left(
\text{podíl variability nezachycený modelem}
\times
\text{úprava za počet pozorování a koeficientů}
\right).
$$

Do vztahu dosadíme hodnoty aditivního modelu:

$$
R^2_{\mathrm{adj}} =
1 -
\left(1 - `r format_value_cz(summary(mod_aditivni)$r.squared, digits = 3, math = TRUE)`\right)
\times
\frac{`r nobs(mod_aditivni)` - 1}
{`r nobs(mod_aditivni)` - `r length(coef(mod_aditivni))`}
=
`r format_value_cz(summary(mod_aditivni)$adj.r.squared, digits = 3, math = TRUE)`.
$$

Po výpočtu obou částí v závorce dostaneme:

$$
R^2_{\mathrm{adj}} =
1 -
\left(
`r format_value_cz(podil_nezachyceny_aditivni, digits = 3, math = TRUE)`
\times
`r format_value_cz(faktor_upravy_aditivni, digits = 3, math = TRUE)`
\right)
=
`r format_value_cz(summary(mod_aditivni)$adj.r.squared, digits = 3, math = TRUE)`.
$$

Obecná podoba stejného výpočtu je:

$$
R^2_{\mathrm{adj}} =
1 -
\left(1 - R^2\right)
\times
\frac{n - 1}{n - k}.
$$

Faktor $\frac{n - 1}{n - k}$ je při přidání koeficientů větší. Pokud nové koeficienty nesníží residuální variabilitu dostatečně, adjustované R² se zmenší.

::: {.box-question title="Označte význam každé části"}
Ve výpočtu výše najděte část popisující shodu modelu, velikost datového souboru a počet odhadovaných koeficientů. Která část se změní, když přidáme `r render_glossary_term("interakce-predikatoru", display = "interakci")`?
:::

Shodu vyjadřuje $1-R^2$, velikost datového souboru představuje $n$ a počet koeficientů představuje $k$. Po přidání interakce se změní R² i $k$; počet pozorování $n$ zůstane stejný.

## Kdy může adjustované R² klesnout?

Porovnáme `r render_glossary_term("r2", display = "obyčejné R²")` po přidání `r render_glossary_term("interakce-predikatoru", display = "interakce")`. Ta přidává další `r render_glossary_term("sklon", display = "sklony")` a má smysl jen tehdy, pokud dostatečně změní `r render_glossary_term("residuum", display = "residua")` vzhledem k počtu nových `r render_glossary_term("parametry-modelu", display = "koeficientů")`.

```{r}
#| label: zobrazit-adjustovane-r2-kandidatu
#| echo: false
#| results: asis

knitr::kable(
  x = data_porovnani_modelu[
    ,
    c(
      "model",
      "pocet_koeficientu",
      "r2",
      "adjustovane_r2"
    )
  ],
  digits = c(0, 0, 3, 3),
  col.names = c(
    "Kandidátní model",
    "Počet koeficientů",
    "R²",
    "Adjustované R²"
  ),
  caption = "Obyčejné R² vždy roste, ale adjustované R² může po přidání málo užitečných koeficientů klesnout."
)
```

Přechod z aditivního modelu k interakčnímu přidá dva sklony. Obyčejné R² se zvýší z `r format_value_cz(summary(mod_aditivni)$r.squared, digits = 3)` na `r format_value_cz(summary(mod_interakce)$r.squared, digits = 3)`, ale adjustované R² klesne z `r format_value_cz(summary(mod_aditivni)$adj.r.squared, digits = 3)` na `r format_value_cz(summary(mod_interakce)$adj.r.squared, digits = 3)`. Zisk ve shodě je pro tuto úpravu příliš malý.

::: {.box-question title="Předpověď před výsledkem"}
Kdyby nový koeficient téměř nezměnil residua, co se stane s obyčejným a adjustovaným R²? Vysvětlete obě odpovědi pomocí jejich výpočtu.
:::

Obyčejné R² se nezmenší, protože RSS se po přidání koeficientu nemůže zvětšit. Adjustované R² však pravděpodobně klesne: téměř stejný nezachycený podíl se násobí větší úpravou za počet koeficientů.

## Jak vyvážit shodu a jednoduchost modelu?

`r render_glossary_term("r2", display = "R²")` hodnotí shodu s pozorovanými daty, ale nezvýhodňuje jednodušší vysvětlení. Při `r render_glossary_term("porovnani-modelu", display = "porovnání modelů")` proto používáme také **princip úspornosti (parsimony)**: složitější model má být upřednostněn jen tehdy, když přinese dostatečné zlepšení shody.

<!-- TODO: add slug princip-uspornosti to pojmy.yaml -->

::: {.box-question title="Kdy stojí složitost za to?"}
Posun mezi druhy a různé sklony popisují stromy téměř stejně dobře (R² `r format_value_cz(summary(object = mod_aditivni)$r.squared, digits = 3)` a `r format_value_cz(summary(object = mod_interakce)$r.squared, digits = 3)`). První model potřebuje `r length(coef(object = mod_aditivni))` `r render_glossary_term("parametry-modelu", display = "koeficienty")`, druhý `r length(coef(object = mod_interakce))`. Který byste předběžně upřednostnili a jakou další informaci byste chtěli znát?
:::

Předběžně bychom upřednostnili jednodušší model. Potřebujeme však posoudit, zda malé zlepšení shody složitějšího modelu vyváží `r length(coef(object = mod_interakce)) - length(coef(object = mod_aditivni))` koeficienty navíc.

## Jak `r render_glossary_term("aic", display = "AIC")` vyjadřuje úspornost modelu?

**Akaikeho informační kritérium (Akaike information criterion, AIC)** je alternativní způsob `r render_glossary_term("porovnani-modelu", display = "porovnání kandidátů")`. V jednom čísle vyvažuje:

- shodu modelu s pozorovanými daty;
- penalizaci za složitost modelu.

Lepší shoda AIC snižuje, zatímco další odhadované veličiny je zvyšují. Proto při srovnání stejných dat preferujeme **nižší AIC**. AIC nevyjadřuje podíl `r render_glossary_term("variabilita", display = "variability")` a nemá stejnou škálu jako R².

## Co nám říkají rozdíly `r render_glossary_term("aic", display = "AIC")` mezi kandidáty?

```{r}
#| label: pripravit-porovnani-aic
#| echo: false
#| output: false

data_porovnani_modelu$aic <-
  purrr::map_dbl(
    .x = list_modely,
    .f = AIC
  )

nejnizsi_aic <-
  min(data_porovnani_modelu$aic)

data_porovnani_modelu$rozdil_od_minima_aic <-
  data_porovnani_modelu$aic -
  nejnizsi_aic
```

Hodnotu AIC jednoho modelu získáme přímo funkcí `AIC()`:

```{r}
#| label: zobrazit-aic-spolecneho-modelu
#| echo: true
#| eval: true

# Vypočítáme AIC společné přímky
AIC(object = mod_spolecny)
```

Samotné AIC `r format_value_cz(AIC(mod_spolecny), digits = 2)` nelze označit za „dobré“ nebo „špatné“. Jde o relativní kritérium: číslo získává význam až ve `r render_glossary_term("porovnani-modelu", display = "srovnání předem určených kandidátů")`.

Porovnávané modely přitom musí používat:

- stejnou `r render_glossary_term("odezva", display = "odezvu")`;
- stejná `r render_glossary_term("pozorovani", display = "pozorování")`;
- stejný způsob pravděpodobnostního popisu `r render_glossary_term("residuum", display = "residuí")` — zde u všech kandidátů předpokládáme residua s `r render_glossary_term("rozdeleni", display = "normálním rozdělením")` a stejnou základní strukturou.

Kandidáty porovnáme tak, že od AIC každého z nich odečteme nejnižší AIC v celé sadě:

$$
\text{rozdíl od nejnižší AIC} =
\text{AIC kandidáta} -
\text{nejnižší AIC v sadě}.
$$

Pro společnou přímku dosadíme její AIC a nejnižší AIC ze tří kandidátů:

$$
\text{rozdíl od nejnižší AIC} =
`r format_value_cz(data_porovnani_modelu$aic[1], digits = 2, math = TRUE)` -
`r format_value_cz(nejnizsi_aic, digits = 2, math = TRUE)`.
$$

Po odečtení dostaneme:

$$
\text{rozdíl od nejnižší AIC} =
`r format_value_cz(data_porovnani_modelu$rozdil_od_minima_aic[1], digits = 2, math = TRUE)`.
$$

Tento rozdíl se často označuje jako **ΔAIC (delta AIC)**.

<!-- TODO: add slug delta-aic to pojmy.yaml -->

Obecně jej zapíšeme:

$$
\Delta AIC_i = AIC_i - AIC_{\min}.
$$

Symbol $i$ označuje konkrétní kandidátní model, $AIC_i$ jeho hodnotu AIC a $AIC_{\min}$ nejnižší AIC v celé kandidátní sadě. Nejlépe podporovaný kandidát má ΔAIC rovné 0; ostatní hodnoty ukazují, o kolik je jejich AIC větší. V této lekci záměrně nepoužíváme univerzální práh, který by rozhodnutí automatizoval.

```{r}
#| label: zobrazit-relativni-aic
#| echo: false
#| results: asis

knitr::kable(
  x = data_porovnani_modelu[
    ,
    c(
      "model",
      "aic",
      "rozdil_od_minima_aic"
    )
  ],
  digits = c(0, 2, 2),
  col.names = c(
    "Kandidátní model",
    "AIC",
    "Rozdíl od nejnižší AIC (ΔAIC)"
  ),
  caption = "AIC a ΔAIC v předem určené kandidátní sadě."
)
```

::: {.box-question title="Proč nestačí jedno číslo?"}
Kolega oznámí: „Můj model má AIC 395, takže je dobrý.“ Jaká informace chybí, abyste mohli jeho tvrzení posoudit?
:::

Chybí AIC dalších modelů odhadnutých pro stejnou odezvu a stejná pozorování. Kandidáty porovnávají rozdíly mezi těmito hodnotami, ne jedna hodnota sama.

::: {.box-warning title="Nejlepší kandidát nemusí být dobrý model"}
Porovnání vybírá nejlépe podporovaný model **pouze z modelů, které jsme zařadili do kandidátní sady**. Neříká, že je tento model obecně dobrý ani že jsme našli nejlepší možný popis dat. Jiný model, který jsme vůbec neposuzovali, může stejná data popsat podstatně lépe.
:::

::: {.box-data title="Doplňující: Z čeho se AIC skládá?" collapse="true"}
AIC sčítá dvě části. První vyjadřuje **nesoulad modelu s daty**: vychází z takzvané věrohodnosti (*likelihood*), tedy z toho, jak pravděpodobná by byla pozorovaná data, kdyby model platil. Čím lépe model data vystihuje, tím je tato část menší. Druhá část je **penalizace za složitost**: AIC přičítá 2 za každou veličinu, kterou model odhaduje.

U modelu z `lm()` se kromě regresních `r render_glossary_term("parametry-modelu", display = "koeficientů")` odhaduje ještě jedna veličina: `r render_glossary_term("smerodatna-odchylka", display = "směrodatná odchylka")` residuí. Společná přímka proto odhaduje `r length(coef(object = mod_spolecny))` koeficienty a jednu směrodatnou odchylku, celkem `r length(coef(object = mod_spolecny)) + 1` veličiny. Ručně AIC počítat nemusíte; následující kód jen ukazuje, že `AIC()` sčítá právě tyto dvě části.

```{r}
#| label: rozlozit-aic
#| echo: true
#| eval: true

# Nesoulad s daty: −2 × logaritmus věrohodnosti modelu
nesoulad_spolecny <-
  -2 * as.numeric(logLik(object = mod_spolecny))

# Počet odhadovaných veličin: koeficienty a směrodatná odchylka residuí
pocet_velicin_spolecny <-
  length(coef(object = mod_spolecny)) + 1

# Součet obou částí
nesoulad_spolecny + 2 * pocet_velicin_spolecny

# Stejnou hodnotu vrátí funkce AIC()
AIC(object = mod_spolecny)
```
:::

## Ukazují `r render_glossary_term("r2", display = "R²")`, adjustované R² a `r render_glossary_term("aic", display = "AIC")` na stejný model?

Každá míra zachází se shodou a složitostí jinak. Jejich společná tabulka proto není hlasováním, ale podkladem pro argument.

```{r}
#| label: zobrazit-uplne-porovnani-modelu
#| echo: false
#| results: asis

knitr::kable(
  x = data_porovnani_modelu[
    ,
    c(
      "model",
      "pocet_koeficientu",
      "r2",
      "adjustovane_r2",
      "aic"
    )
  ],
  digits = c(0, 0, 3, 3, 2),
  col.names = c(
    "Kandidátní model",
    "Počet koeficientů",
    "R²",
    "Adjustované R²",
    "AIC"
  ),
  caption = "Tři různé souhrny stejné, předem definované sady kandidátních modelů."
)
```

::: {.box-question title="Čtěte tabulku jako argument"}
Který model má nejvyšší R²? Který má nejvyšší adjustované R² a nejnižší AIC? Vysvětlete neshodu pomocí toho, jak každá míra zachází se shodou a složitostí, ne jako tři rovnocenné hlasy o tom, který model je „pravdivý“.
:::

`r render_glossary_term("interakce-predikatoru", display = "Interakční model")` má nejvyšší R², protože má nejvíce možností přizpůsobit se datům. Aditivní model má nejvyšší adjustované R² a nejnižší AIC, protože dodatečné `r render_glossary_term("sklon", display = "sklony")` interakčního modelu nepřinesly dostatečné zlepšení shody. Žádná z těchto hodnot neudává pravděpodobnost, že je model pravdivý.

## Stačí dobré skóre k tomu, aby byl model dobrý?

Srovnávací skóre nenahrazuje `r render_glossary_term("diagnostika-modelu", display = "diagnostiku modelu")`. U kandidáta podporovaného adjustovaným R² a AIC kontrolujeme vztah `r render_glossary_term("residuum", display = "residuí")` k `r render_glossary_term("odhadnuta-hodnota", display = "odhadnutým hodnotám")`, `r render_glossary_term("qq-graf", display = "Q–Q graf")` a vliv jednotlivých `r render_glossary_term("pozorovani", display = "pozorování")`.

Do grafu residuí přidáme vyhlazenou křivku, kterou spočítá funkce `lowess()`. Křivka v každém místě grafu sleduje, kde residua v okolí v průměru leží. Nic nového neodhaduje, jen pomáhá vidět vzor, který tvoří body.

```{r}
#| label: zkontrolovat-residua-aditivniho-modelu
#| echo: true
#| eval: true
#| fig-cap: "Residua aditivního modelu podle odhadnuté výšky. Oranžová křivka pouze vyhlazuje průběh bodů, aby byl vzor lépe vidět."
#| fig-alt: "Residua aditivního modelu proti odhadnuté výšce s oranžovou vyhlazenou křivkou. Residua nízkých i nejvyšších odhadnutých výšek leží převážně pod nulou, středně vysokých nad nulou, takže křivka tvoří oblouk."

# Zkontrolujeme, zda residua leží kolem nuly bez vzoru
plot(
  x = fitted(object = mod_aditivni),
  y = residuals(object = mod_aditivni),
  col = vec_barvy_druhu[as.character(data_stromy$druh_stromu)],
  pch = vec_symboly_druhu[as.character(data_stromy$druh_stromu)],
  cex = 1.2,
  cex.axis = 1.2,
  cex.lab = 1.2,
  xlab = "Odhadnutá výška stromu (m)",
  ylab = "Residuum (m)"
)

abline(
  h = 0,
  col = "#2E2E2E",
  lty = 2
)

# Vyhlazená křivka sleduje průměrnou polohu residuí
lines(
  x = lowess(
    x = fitted(object = mod_aditivni),
    y = residuals(object = mod_aditivni)
  ),
  col = "#F3A712",
  lwd = 3
)
```

```{r}
#| label: spocitat-oblouk-residui
#| echo: false

# Teaching input: two fitted heights that split the trees into
#   low, middle and tall groups for describing the arch
vec_hranice_pasem_m <- c(16, 32)

# Mean residual in each band of fitted heights, for all three candidates
data_oblouk_residui <-
  list(
    spolecny = mod_spolecny,
    aditivni = mod_aditivni,
    interakce = mod_interakce
  ) |>
  purrr::imap(
    .f = ~ tibble::tibble(
      model = .y,
      pasmo = as.integer(
        cut(
          x = fitted(object = .x),
          breaks = c(-Inf, vec_hranice_pasem_m, Inf)
        )
      ),
      residuum = residuals(object = .x)
    )
  ) |>
  purrr::list_rbind() |>
  dplyr::group_by(model, pasmo) |>
  dplyr::summarise(
    prumerne_residuum = mean(residuum),
    .groups = "drop"
  )

# Signed band means of each model as text, e.g. "−5,9; +3,1 a −2,1"
vec_text_oblouk <-
  data_oblouk_residui |>
  dplyr::arrange(model, pasmo) |>
  dplyr::mutate(
    hodnota = format_value_cz(
      x = prumerne_residuum,
      digits = 1,
      signed = TRUE
    )
  ) |>
  dplyr::group_by(model) |>
  dplyr::summarise(
    text = stringr::str_glue("{hodnota[1]}; {hodnota[2]} a {hodnota[3]}"),
    .groups = "drop"
  ) |>
  tibble::deframe()

vec_residua_pasem_aditivni <-
  data_oblouk_residui |>
  dplyr::filter(model == "aditivni") |>
  dplyr::arrange(pasmo) |>
  dplyr::pull(prumerne_residuum)
```

Residua neleží kolem nuly náhodně. U stromů s odhadnutou výškou pod `r vec_hranice_pasem_m[[1]]` m je průměrné residuum `r format_value_cz(vec_residua_pasem_aditivni[[1]], digits = 1, signed = TRUE)` m, mezi `r vec_hranice_pasem_m[[1]]` a `r vec_hranice_pasem_m[[2]]` m `r format_value_cz(vec_residua_pasem_aditivni[[2]], digits = 1, signed = TRUE)` m a nad `r vec_hranice_pasem_m[[2]]` m `r format_value_cz(vec_residua_pasem_aditivni[[3]], digits = 1, signed = TRUE)` m. Residua tak tvoří **oblouk**: nízké i nejvyšší stromy model nadhodnocuje, středně vysoké podhodnocuje. Přímka tedy nevystihuje tvar vztahu mezi průměrem kmene a výškou; výška s průměrem kmene nejspíš roste stále pomaleji. Ve stejných pásmech mají stejný oblouk i residua společné přímky (průměrná residua `r vec_text_oblouk[["spolecny"]]` m) a modelu s různými sklony (`r vec_text_oblouk[["interakce"]]` m). Oblouk tedy nevyřeší ani druh stromu, ani rozdílné sklony.

```{r}
#| label: zkontrolovat-normalitu-aditivniho-modelu
#| echo: true
#| eval: true
#| fig-cap: "Q–Q graf residuí aditivního modelu."
#| fig-alt: "Q–Q graf residuí aditivního modelu: body tvoří mírně esovitou křivku a na obou koncích se odchylují od přímky."

# Zkontrolujeme tvar rozdělení residuí
qqnorm(
  y = residuals(object = mod_aditivni),
  pch = 19,
  cex.axis = 1.2,
  cex.lab = 1.2,
  main = "Q–Q graf residuí",
  xlab = "Teoretické kvantily",
  ylab = "Kvantily residuí (m)"
)

qqline(
  y = residuals(object = mod_aditivni),
  col = "#2E2E2E",
  lwd = 2
)
```

Body v Q–Q grafu se na obou koncích od přímky odchylují. Část této odchylky může souviset s obloukem: residua, která zanechává nevhodný tvar modelu, nemusí mít rozdělení předpokládané lineárním modelem.

Zbývá otázka, zda některé pozorování model neovlivňuje neobvykle silně. `r render_glossary_term("cookova-vzdalenost", display = "Cookova vzdálenost")` (*Cook's distance*) shrnuje, jak silně by se odhadnuté hodnoty změnily po vynechání jednoho pozorování. Orientační čára $4/n$ není automatický test; pouze označuje stromy, které si zaslouží bližší kontrolu.

```{r}
#| label: zkontrolovat-vliv-aditivniho-modelu
#| echo: true
#| eval: true
#| fig-cap: "Cookovy vzdálenosti pro aditivní model. Přerušovaná čára zobrazuje orientační hodnotu 4/n."
#| fig-alt: "Svislé čárky Cookových vzdáleností všech stromů v pořadí tabulky s přerušovanou orientační čarou 4/n; několik stromů ji přesahuje."

# Vypočítáme vliv jednotlivých pozorování
vec_cookova_vzdalenost <-
  cooks.distance(model = mod_aditivni)

hranice_cook <-
  4 / nobs(object = mod_aditivni)

# Zobrazíme, která pozorování mají největší vliv
plot(
  x = vec_cookova_vzdalenost,
  type = "h",
  lwd = 2,
  cex.axis = 1.2,
  cex.lab = 1.2,
  xlab = "Pořadí stromu v tabulce",
  ylab = "Cookova vzdálenost"
)

abline(
  h = hranice_cook,
  col = "#F3A712",
  lty = 2,
  lwd = 2
)
```

Nad orientační čarou leží `r sum(vec_cookova_vzdalenost > hranice_cook)` pozorování. To není důvod je odstranit; je to důvod vrátit se k jejich měření, biologickému významu a vlivu na závěr.

::: {.box-question title="Rozhodnutí po kontrole modelu"}
Podporují graf residuí, Q–Q graf a Cookovy vzdálenosti závěr ze srovnávací tabulky, nebo vyžadují opatrnější formulaci? Uveďte jeden důvod pro zvolený model a jedno omezení.
:::

::: {.box-result collapse="true" title="Možná odpověď"}
Srovnávací tabulka podporuje aditivní model jako nejlepšího ze tří kandidátů: zachytí rozdíly mezi druhy a má nejvyšší adjustované R² i nejnižší AIC. Diagnostika ale vyžaduje opatrnější formulaci. Residua tvoří oblouk, takže žádná z přímek nevystihuje tvar vztahu, a Q–Q graf není přímý. Stromy nad čarou $4/n$ je vhodné zkontrolovat, ne automaticky vyřadit.
:::

## Který model lze pro tyto stromy nejlépe obhájit?

Aditivní model s průměrem kmene a druhem stromu představuje pro tato data nejlépe obhajitelný model **ze tří posuzovaných kandidátů**. Oproti jedné společné přímce zachytí rozdíly mezi druhy. `r render_glossary_term("interakce-predikatoru", display = "Interakční model")` sice dosáhne nejvyššího `r render_glossary_term("r2", display = "R²")`, ale jeho dodatečné `r render_glossary_term("sklon", display = "sklony")` nepřinesou dostatečný zisk podle adjustovaného R² ani `r render_glossary_term("aic", display = "AIC")`.

Tento závěr zůstává podmíněný. Neříká, že je aditivní model nejlepší ze všech možných modelů. `r render_glossary_term("residuum", display = "Residua")` všech tří kandidátů tvoří oblouk, takže přímky nevystihují zakřivený tvar vztahu. Data navíc obsahují omezený počet stromů, `r render_glossary_term("qq-graf", display = "Q–Q graf")` není přímý a pozorovací studie sama nedokládá příčinný `r render_glossary_term("efekt", display = "efekt")` druhu.

::: {.box-result title="Odpověď na hlavní otázku"}
Dobrý model nepoznáme podle jediného čísla. Obhajoba musí propojit biologickou otázku, shodu s daty, složitost, `r render_glossary_term("diagnostika-modelu", display = "diagnostiku")`, interpretovatelnost a omezení dat.
:::

{{< pagebreak >}}

## Jak `r render_glossary_term("porovnani-modelu", display = "porovnávat modely")` bez automatického hledání?

1. Položte biologickou otázku.
2. Prohlédněte a zobrazte data.
3. Definujte malou sadu biologicky zdůvodněných kandidátů.
4. Použijte pro všechny kandidáty stejnou `r render_glossary_term("odezva", display = "odezvu")` a stejná `r render_glossary_term("pozorovani", display = "pozorování")`.
5. Porovnejte odhadnuté vztahy, `r render_glossary_term("residuum", display = "residua")`, `r render_glossary_term("r2", display = "R²")`, adjustované R² a `r render_glossary_term("aic", display = "AIC")`.
6. Zkontrolujte předpoklady a vlivná pozorování.
7. Obhajte závěr a pojmenujte jeho omezení.

::: {.box-warning title="Automatické hledání nenahrazuje biologickou otázku"}
Funkce pro automatické přidávání a odebírání členů mohou vytvořit dlouhý seznam modelů, ale samy neurčí, které modely mají biologický smysl. V této lekci proto nepoužíváme krokový výběr ani vyhledávání všech možných kombinací.
:::

## Co bude následovat dál

Všechny tři kandidáty jsme kreslili jako přímky, a jejich `r render_glossary_term("residuum", display = "residua")` přesto tvoří oblouk. Problém tedy není jen ve volbě `r render_glossary_term("prediktor", display = "prediktorů")`: vztah výšky a průměru kmene je `r render_glossary_term("nelinearita", display = "nelineární")`. V následující lekci se naučíme popisovat tvar `r render_glossary_term("rozdeleni", display = "rozdělení")` hodnot i tvar vztahu a uvidíme, jak může vhodná `r render_glossary_term("transformace", display = "transformace")` proměnných pomoci lineárnímu modelu popsat i zakřivený vztah. Kandidátní tvary vztahu budeme stejně jako v této lekci porovnávat až poté, co je biologicky zdůvodníme.

## Shrnutí

- `r render_glossary_term("r2", display = "R²")` je podíl `r render_glossary_term("variabilita", display = "variability")` `r render_glossary_term("odezva", display = "odezvy")`, který model zachytil v datech; počítá se z celkového součtu čtverců a `r render_glossary_term("soucet-ctvercu-residui", display = "součtu čtverců residuí")`.
- R² po přidání členu do vnořeného `r render_glossary_term("linearni-model", display = "lineárního modelu")` nemůže klesnout, a proto samo nevybírá model.
- Adjustované R² a `r render_glossary_term("aic", display = "AIC")` vyvažují shodu a složitost, každé jinak; AIC má smysl jen jako rozdíl mezi modely se stejnou odezvou a stejnými `r render_glossary_term("pozorovani", display = "pozorováními")`.
- Kandidáty definujeme podle biologické otázky předtím, než porovnáme jejich skóre.
- Konečná obhajoba musí zahrnout také `r render_glossary_term("diagnostika-modelu", display = "diagnostiku")`, interpretovatelnost a omezení dat.

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

Vyberte si biologickou otázku ze svého oboru. Jaké dva nebo tři kandidátní modely byste sestavili ještě před analýzou a co by vás přesvědčilo, abyste dali přednost jednoduššímu z nich?

::: {.box-result collapse="true" title="Možná odpověď"}
Například otázka, jak hmotnost ptáčat závisí na jejich věku a na tom, zda hnízdo leží v lese, nebo ve městě. Kandidáty mohou být jedna společná přímka, rovnoběžné přímky pro les a město a přímky s různými sklony, pokud čekáme, že ptáčata v jednom prostředí přibývají rychleji. Jednoduššímu modelu bychom dali přednost, kdyby složitější model zlepšil shodu jen nepatrně (adjustované R² by kleslo a AIC by nebylo nižší), `r render_glossary_term("residuum", display = "residua")` jednoduššího modelu nevykazovala vzor a biologická otázka nevyžadovala rozdílné sklony.
:::

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

Pojmy použité v těchto materiálech a jejich definice:

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