---
format:
  revealjs:
    slide-number: true
    theme: [default, ../theme/presentation_theme.scss, ../theme/presentation_components.scss]
    embed-resources: true
    preview-links: false
    output-file: "presentation"
    output-ext: "html"
    highlight-style: github-light
    include-in-header: "../theme/fonts-include.html"
filters:
  - ../theme/rn-shorthand.lua
  - ../_extensions/EmilHvitfeldt/roughnotation/rough.lua
---

```{r}
#| label: setup
#| include: false
#| warning: false
#| message: false
options(htmltools.dir.version = FALSE)
knitr::opts_chunk$set(
  fig.width = 10,
  fig.height = 5.6,
  fig.align = "center",
  out.width = "100%",
  warning = FALSE,
  message = FALSE
)

here::i_am("Presentation/presentation.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())

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

# 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
  }
)

path_materials <- "Presentation/Materials"

# Absolute paths: root.dir is the project root, so knitr must not rewrite
#   the image path relative to the Presentation folder.
include_local_figure <- function(data_source) {
  return(
    knitr::include_graphics(
      path = here::here(path_materials, data_source),
      rel_path = FALSE,
      error = TRUE
    )
  )
}

# Species are categories, not semantic roles: brand categorical palette,
#   deliberately distinct from model purple and residual orange. Shapes
#   repeat the distinction without colour; the written materials and the
#   exercise use the same mapping.
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
)

theme_l09 <- theme_biostat(base_size = 19) +
  ggplot2::theme(
    panel.grid.minor = ggplot2::element_blank(),
    legend.position = "bottom"
  )
```

```{r}
#| label: dynamic link generation
#| include: false
project_name <- basename(here::here())
display_name <- stringr::str_glue("Biostatistika - {project_name}")
github_url <- stringr::str_glue(
  "https://github.com/CUNI-NATUR-Biostatistics/{project_name}"
)
slides_url <- stringr::str_glue(
  "https://CUNI-NATUR-Biostatistics.github.io/{project_name}/"
)
```

## Jak poznáme, že je model dobrý? {.course-title .course-title-illustrated .bg-indigo}

::: {.course-title-mark}
![Logo kurzu Biostatistika](../theme/course-logo-reversed.svg) BIOSTATISTIKA
:::

::: {.course-title-number}
L09
:::

::: {.course-title-formal}
Porovnání kandidátních modelů: R², adjustované R² a AIC
:::

::: {.course-title-link}
[cuni-natur-biostatistics.github.io/materialy.html](https://cuni-natur-biostatistics.github.io/materialy.html)
:::

::: {.course-title-illustration}
```{r}
#| label: zobrazit-titulni-ilustraci
#| echo: false
#| fig-alt: "Liška a veverka u terénního stolku s pásmem pod třemi různými jehličnany; opřené dřevěné destičky ukazují tři uspořádání fialových tyček a sova je z větve zamyšleně posuzuje. Ilustrace vytvořená pomocí AI."
# Generated title illustration; prompt and provenance in GENERATED_IMAGES.md.
include_optional_figure(data_source = "stromy_titulni_ilustrace.png")
```
:::

::: {.course-title-caption}
`r ai_disclosure("stromy_titulni_ilustrace.png")`
:::

::: {.notes}
Titulní ilustrace: liška a veverka měří stromy a sova posuzuje tři uspořádání tyček, tedy tři kandidátní modely téhož vztahu. Žádná destička není označená jako vítěz; to je otázka celé přednášky.
:::

## Co vlastně na stromě měříme?

:::: {.layout-stack .layout-center-y}

::: columns
::: {.column width="58%"}
```{r}
#| label: zobrazit-stromy-mereni
#| echo: false
#| out-width: "100%"
#| fig-alt: "Liška na lesní měřicí stanici omotává pásmo kolem kmene ve výšce prsou a veverka zaměřuje vrchol stromu; vedle sebe stojí tři různé jehličnany."
include_optional_figure(data_source = "stromy_mereni.png")
```

::: {.text-size-tiny .text-right}
`r ai_disclosure("stromy_mereni.png")`
:::
:::
::: {.column width="40%"}
::: {.panel-data}
[U každého jehličnanu známe]{.card-title}

- **druh stromu** — tři severoamerické jehličnany;
- **průměr kmene** ve výšce 1,3 m (cm);
- **výšku stromu** (m).
:::
:::
:::

::: {.fragment}
[Mají stejně silné stromy různých druhů stejnou výšku?]{.text-size-heading .text-bold}
:::

::::

::: {.notes}
Krátce: co by na stromě měřili studenti? Liška omotává pásmo kolem kmene ve výšce prsou a veverka zaměřuje vrchol stromu. Otázku o druzích nechte otevřenou; ilustrace ukazuje jen měření, nikoli výsledek porovnání druhů.
:::

{{< include ../pollslive/generated/active.qmd >}}

## Výsledky učení

Po absolvování této přednášky budete schopni:

::: incremental
- vypočítat [R²]{.text-hgl-indigo_velvet} z celkového součtu čtverců a součtu čtverců residuí a vysvětlit, co popisuje;
- vysvětlit, proč [nejvyšší R²]{.text-hgl-orange} automaticky neznamená nejlepší model;
- rozlišit, jak R², [adjustované R²]{.text-hgl-indigo_velvet} a [AIC]{.text-hgl-indigo_velvet} zacházejí se shodou a složitostí, a používat AIC jen k relativnímu porovnání;
- sestavit malou sadu biologicky zdůvodněných kandidátních modelů se stejnou odezvou a stejnými pozorováními;
- [obhájit volbu modelu]{.text-hgl-orange} pomocí biologické otázky, složitosti, diagnostiky a interpretovatelnosti a pojmenovat její omezení.
:::

# Jedna přímka pro všechny stromy?

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

::: columns
::: {.column width="62%"}
```{r}
#| label: nacist-data-stromu
#| echo: true
data_stromy <-
  read.csv(file = "data/allometrie_stromu.csv")

data_stromy$druh_stromu <-
  factor(x = data_stromy$druh_stromu)

vec_promenne <-
  c("druh_stromu", "prumer_kmene_cm", "vyska_stromu_m")

head(
  x = data_stromy[, vec_promenne],
  n = 4
)
```
:::
::: {.column width="36%" .fragment}

<br>

::: {.panel-data .text-size-smaller}
[`r nrow(data_stromy)` stromů]{.card-title}

```{r}
#| label: pocty-druhu
#| echo: false
#| results: asis
vec_pocty_druhu <- table(data_stromy$druh_stromu)
cat(
  paste0(
    "- *", names(vec_pocty_druhu), "*: ", vec_pocty_druhu, " stromů",
    collapse = "\n"
  )
)
```

- chybějící hodnoty: `r sum(is.na(data_stromy[, vec_promenne]))`
- průměr kmene: `r format_value_cz(min(data_stromy$prumer_kmene_cm), digits = 1)`–`r format_value_cz(max(data_stromy$prumer_kmene_cm), digits = 1)` cm
- výška: `r format_value_cz(min(data_stromy$vyska_stromu_m), digits = 1)`–`r format_value_cz(max(data_stromy$vyska_stromu_m), digits = 1)` m
:::
:::
:::

Jiné stromy než minule: tři druhy jehličnanů a [výška je tentokrát odezva]{.text-hgl-graphite}.

::: {.notes}
V minulé lekci výška stromu pomáhala odhadnout objem dřeva střemchy; byla tedy prediktorem. Tady jde o jiné stromy a výšku chceme popsat.
:::

[Data: J. Marshall, University of Idaho · balíček `lgrdata` 0.1.1 (`allometry`) · CC0]{.footer}

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

```{r}
#| label: zobrazit-stromy
#| echo: false
#| fig-height: 5
#| fig-alt: "Bodový graf výšky stromu proti průměru kmene pro jehličnany tří druhů; výška s průměrem roste."
graf_stromy <-
  ggplot2::ggplot(
    data = data_stromy,
    mapping = ggplot2::aes(
      x = prumer_kmene_cm,
      y = vyska_stromu_m,
      colour = druh_stromu,
      shape = druh_stromu
    )
  ) +
  ggplot2::geom_point(size = 3, alpha = 0.85) +
  ggplot2::scale_colour_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)",
    colour = NULL,
    shape = NULL
  ) +
  theme_l09 +
  ggplot2::theme(legend.text = ggplot2::element_text(face = "italic"))

graf_stromy
```

::: {.panel-question}
**Ve dvojici:** Jaký vzorec v datech převládá? Co by mohlo vadit jedné přímce pro všechny stromy?
:::

## Jedna společná přímka pro všechny druhy

::: columns
::: {.column width="42%"}

<br>

```{r}
#| label: odhadnout-spolecny-model
#| echo: true
mod_spolecny <-
  lm(
    formula = vyska_stromu_m ~
      prumer_kmene_cm,
    data = data_stromy
  )

coef(object = mod_spolecny)
```
:::
::: {.column width="56%"}
```{r}
#| label: zobrazit-spolecnou-primku
#| echo: false
#| fig-width: 7
#| fig-height: 5.4
#| fig-alt: "Body stromů tří druhů a jedna společná fialová přímka modelu."
graf_stromy +
  ggplot2::geom_abline(
    intercept = coef(mod_spolecny)[[1]],
    slope = coef(mod_spolecny)[[2]],
    colour = biostat_cols[["indigo_velvet"]],
    linewidth = 1.4
  ) +
  ggplot2::guides(
    colour = ggplot2::guide_legend(nrow = 2),
    shape = ggplot2::guide_legend(nrow = 2)
  )
```
:::
:::

::: {.fragment}
Sklon: o [`r format_value_cz(coef(mod_spolecny)[[2]], digits = 2)` m vyšší]{.text-hgl-indigo_velvet} odhadovaná výška na každý centimetr průměru — [stejně pro všechny druhy]{.text-hgl-orange}.
:::

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

```{r}
#| label: souhrn-spolecneho-modelu
#| echo: true
#| class-output: text-size-tiny
summary(object = mod_spolecny)
```

[[Známe:]{.fragment} [[volání modelu · residua · koeficienty s nejistotou]{.text-hgl-indigo_velvet}]{.fragment}]{.text-size-smaller}

::: {.notes}
V L08 jsme tyto dva řádky záměrně odložili. Hodnoty zatím nečteme — vrátíme se k nim, až je umíme spočítat.
:::

## Zůstává v residuích rozdíl mezi druhy?

```{r}
#| label: pripravit-residua-spolecneho-modelu
#| include: false
data_residua_spolecny <- data.frame(
  odhad = fitted(object = mod_spolecny),
  residuum = residuals(object = mod_spolecny),
  druh_stromu = data_stromy$druh_stromu
)

data_prumerna_residua <-
  aggregate(
    residuum ~ druh_stromu,
    data = data_residua_spolecny,
    FUN = mean
  )

graf_residua_spolecny <-
  ggplot2::ggplot(
    data = data_residua_spolecny,
    mapping = ggplot2::aes(
      x = odhad,
      y = residuum,
      colour = druh_stromu,
      shape = druh_stromu
    )
  ) +
  ggplot2::geom_hline(
    yintercept = 0,
    linetype = "dashed",
    colour = biostat_cols[["graphite"]]
  ) +
  ggplot2::geom_point(size = 3, alpha = 0.85) +
  ggplot2::scale_colour_manual(values = vec_barvy_druhu) +
  ggplot2::scale_shape_manual(values = vec_symboly_druhu) +
  ggplot2::scale_y_continuous(
    limits = range(data_residua_spolecny$residuum) + c(-1, 1)
  ) +
  ggplot2::labs(
    x = "Odhadnutá výška stromu (m)",
    y = "Residuum (m)",
    colour = NULL,
    shape = NULL
  ) +
  theme_l09 +
  ggplot2::theme(legend.text = ggplot2::element_text(face = "italic")) +
  ggplot2::guides(
    colour = ggplot2::guide_legend(nrow = 2),
    shape = ggplot2::guide_legend(nrow = 2)
  )
```

::: columns
::: {.column width="64%"}
::: {.r-stack}
```{r}
#| label: zobrazit-residua-spolecneho-modelu
#| echo: false
#| fig-width: 8
#| fig-height: 6
#| fig-alt: "Residua společné přímky proti odhadnuté výšce, barevně podle druhu stromu."
graf_residua_spolecny
```

::: {.fragment}
```{r}
#| label: zobrazit-residua-s-prumery-druhu
#| echo: false
#| fig-width: 8
#| fig-height: 6
#| fig-alt: "Stejná residua s vodorovnými čarami průměrného residua každého druhu; Pinus monticola leží převážně nad nulou."
graf_residua_spolecny +
  ggplot2::geom_hline(
    data = data_prumerna_residua,
    mapping = ggplot2::aes(yintercept = residuum, colour = druh_stromu),
    linewidth = 1.4
  )
```
:::
:::
:::
::: {.column width="34%"}

<br>

::: {.panel-question}
Leží residua některého druhu převážně nad nulou, nebo pod ní?
:::

::: {.fragment}
Průměrné residuum *P. monticola*: [`r format_value_cz(data_prumerna_residua$residuum[data_prumerna_residua$druh_stromu == "Pinus monticola"], digits = 1, signed = TRUE)` m]{.text-hgl-orange}

Společná přímka [přehlíží rozdíly mezi druhy]{.text-hgl-orange}.
:::
:::
:::

# Tři představy o růstu stromů

## Jaké jiné přímky by mohly stromy popsat?

::: columns
::: {.column width="62%"}
```{r}
#| label: zobrazit-stromy-pro-nacrt
#| echo: false
#| fig-width: 8
#| fig-height: 6
#| fig-alt: "Stejný bodový graf stromů bez přímek, připravený pro náčrt alternativních modelů."
graf_stromy +
  ggplot2::guides(
    colour = ggplot2::guide_legend(nrow = 2),
    shape = ggplot2::guide_legend(nrow = 2)
  )
```
:::
::: {.column width="36%"}

<br>

::: {.panel-question}
**Ve dvojici:** Načrtněte dvě jiné možnosti než jednu společnou přímku.

Co by každá z nich tvrdila o stromech?
:::

::: {.fragment}
Pomůcka: modely z minulé lekce — [`+`]{.text-hgl-indigo_velvet} a [`*`]{.text-hgl-indigo_velvet}.
:::
:::
:::

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

```{r}
#| label: zobrazit-glyfy-kandidatu
#| echo: false
#| fig-height: 2.2
#| fig-alt: "Tři schémata: jedna přímka, tři rovnoběžné přímky a tři přímky s různými sklony."
vec_nazvy_modelu <- c(
  "Společná přímka",
  "Posun mezi druhy",
  "Různé sklony"
)

data_glyfy <- data.frame(
  model = factor(
    rep(vec_nazvy_modelu, times = c(1, 3, 3)),
    levels = vec_nazvy_modelu
  ),
  skupina = c("a", "a", "b", "c", "a", "b", "c"),
  intercept = c(1, 2, 1, 0, 0.6, 0.3, 0),
  sklon = c(1, 1, 1, 1, 1.5, 1, 0.55)
)

ggplot2::ggplot(data = data_glyfy) +
  ggplot2::geom_abline(
    mapping = ggplot2::aes(intercept = intercept, slope = sklon, group = skupina),
    colour = biostat_cols[["indigo_velvet"]],
    linewidth = 1.6
  ) +
  ggplot2::facet_wrap(facets = ggplot2::vars(model), nrow = 1) +
  ggplot2::coord_cartesian(xlim = c(0, 3), ylim = c(0, 5)) +
  ggplot2::theme_void(base_size = 20) +
  ggplot2::theme(
    strip.text = ggplot2::element_text(face = "bold", margin = ggplot2::margin(b = 6)),
    panel.border = ggplot2::element_rect(colour = biostat_cols[["light_gray"]], fill = NA)
  )
```

::: columns
::: {.column width="32%"}
`vyska ~ prumer`

Všechny druhy sdílejí [stejný intercept i sklon]{.text-hgl-indigo_velvet}.
:::
::: {.column width="32%"}
`vyska ~ prumer + druh`

Při stejném průměru se druhy mohou lišit [výškou]{.text-hgl-indigo_velvet}, sklon sdílejí.
:::
::: {.column width="32%"}
`vyska ~ prumer * druh`

Vztah průměru a výšky se může [mezi druhy lišit]{.text-hgl-indigo_velvet}.
:::
:::

::: {.fragment .takeaway-strip}
Tuto sadu kandidátů určíme [předem]{.text-hgl-orange} — z biologické otázky, ne z výsledků.
:::

## Všechny tři modely odhadneme ze stejných `r nrow(data_stromy)` stromů

:::: {.layout-stack .layout-center-y}

```{r}
#| label: odhadnout-kandidatni-modely
#| echo: true
# Všechny tři kandidáty zapíšeme vedle sebe
mod_spolecny <-
  lm(
    formula = vyska_stromu_m ~ prumer_kmene_cm,
    data = data_stromy
  )

mod_aditivni <-
  lm(
    formula = vyska_stromu_m ~ prumer_kmene_cm + druh_stromu,
    data = data_stromy
  )

mod_interakce <-
  lm(
    formula = vyska_stromu_m ~ prumer_kmene_cm * druh_stromu,
    data = data_stromy
  )
```

::: {.fragment}
::: {.panel-data}
Stejná odezva [`vyska_stromu_m`]{.text-hgl-graphite} · stejných [`r nobs(mod_spolecny)`]{.text-hgl-graphite} stromů v každém modelu
:::
:::

::::

::: {.notes}
Tato podmínka bude důležitá později u AIC; teď ji jen pojmenujeme.
:::

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

```{r}
#| label: pripravit-predikce-kandidatu
#| include: 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 = 100
  ),
  druh_stromu = levels(data_stromy$druh_stromu)
)

list_kandidati <- list(mod_spolecny, mod_aditivni, mod_interakce)

data_predikce <-
  purrr::map2(
    .x = list_kandidati,
    .y = vec_nazvy_modelu,
    .f = ~ data.frame(
      data_mrizka,
      odhad = predict(object = .x, newdata = data_mrizka),
      model = .y
    )
  ) |>
  purrr::list_rbind()

data_predikce$model <- factor(data_predikce$model, levels = vec_nazvy_modelu)

# The common model has one line: draw it once, in model purple
data_predikce_spolecny <-
  data_predikce[
    data_predikce$model == vec_nazvy_modelu[1] &
      data_predikce$druh_stromu == levels(data_stromy$druh_stromu)[1],
  ]

data_body_kandidati <-
  purrr::map2(
    .x = list_kandidati,
    .y = vec_nazvy_modelu,
    .f = ~ data.frame(
      data_stromy[, c("druh_stromu", "prumer_kmene_cm", "vyska_stromu_m")],
      odhad = fitted(object = .x),
      model = .y
    )
  ) |>
  purrr::list_rbind()

data_body_kandidati$model <- factor(data_body_kandidati$model, levels = vec_nazvy_modelu)

graf_kandidati <-
  ggplot2::ggplot(
    data = data_body_kandidati,
    mapping = ggplot2::aes(x = prumer_kmene_cm, y = vyska_stromu_m)
  ) +
  ggplot2::geom_point(
    mapping = ggplot2::aes(colour = druh_stromu, shape = druh_stromu),
    size = 2.2,
    alpha = 0.8
  ) +
  ggplot2::geom_line(
    data = data_predikce[data_predikce$model != vec_nazvy_modelu[1], ],
    mapping = ggplot2::aes(y = odhad, colour = druh_stromu, group = druh_stromu),
    linewidth = 1.2
  ) +
  ggplot2::geom_line(
    data = data_predikce_spolecny,
    mapping = ggplot2::aes(y = odhad),
    colour = biostat_cols[["indigo_velvet"]],
    linewidth = 1.3
  ) +
  ggplot2::facet_wrap(facets = ggplot2::vars(model), nrow = 1) +
  ggplot2::scale_colour_manual(values = vec_barvy_druhu) +
  ggplot2::scale_shape_manual(values = vec_symboly_druhu) +
  ggplot2::labs(
    x = "Průměr kmene (cm)",
    y = "Výška stromu (m)",
    colour = NULL,
    shape = NULL
  ) +
  theme_l09 +
  ggplot2::theme(
    legend.text = ggplot2::element_text(face = "italic"),
    strip.text = ggplot2::element_text(face = "bold")
  )
```

```{r}
#| label: zobrazit-kandidati
#| echo: false
#| fig-height: 4.6
#| fig-alt: "Tři panely se stejnými stromy: jedna společná přímka, tři rovnoběžné přímky a tři přímky s různými sklony."
graf_kandidati
```

::: columns
::: {.column width="32%"}
::: {.card-compact}
**A** · Společná přímka
:::
:::
::: {.column width="32%"}
::: {.card-compact}
**B** · Posun mezi druhy
:::
:::
::: {.column width="32%"}
::: {.card-compact}
**C** · Různé sklony
:::
:::
:::

## Nejblíže procházejí přímky s různými sklony

```{r}
#| label: zobrazit-kandidati-s-residui
#| echo: false
#| fig-height: 4.6
#| fig-alt: "Stejné tři panely se svislými oranžovými úsečkami residuí a součtem čtverců residuí v záhlaví; nejmenší je u modelu Různé sklony."
vec_popisky_rss <- stats::setNames(
  object = paste0(
    vec_nazvy_modelu,
    "\nRSS = ",
    format_value_cz(purrr::map_dbl(list_kandidati, deviance), digits = 0),
    " m²"
  ),
  nm = vec_nazvy_modelu
)

graf_kandidati +
  ggplot2::facet_wrap(
    facets = ggplot2::vars(model),
    nrow = 1,
    labeller = ggplot2::as_labeller(vec_popisky_rss)
  ) +
  ggplot2::geom_segment(
    mapping = ggplot2::aes(
      xend = prumer_kmene_cm,
      y = odhad,
      yend = vyska_stromu_m
    ),
    colour = biostat_cols[["orange"]],
    linewidth = 0.7,
    alpha = 0.9
  )
```

::: {.fragment .text-center}
[Je proto nejlepší?]{.text-size-heading .text-bold .text-hgl-orange}
:::

::: {.notes}
Oranžové úsečky jsou residua v každém panelu; RSS v záhlaví je známý součet jejich čtverců. Rozdíly jsou malé — a právě proto potřebujeme víc než nejmenší RSS.
:::

## Podle čeho poznáme lepší model?

:::: {.layout-stack .layout-center-y}

[Podle čeho byste mezi modely vybírali?]{.text-size-heading .text-bold}

::: incremental
- jak těsně model odpovídá [pozorovaným datům]{.text-hgl-graphite};
- kolik [koeficientů]{.text-hgl-indigo_velvet} potřebuje;
- zda v [residuích]{.text-hgl-orange} nezůstává systematický vzorec;
- zda odpovídá biologické otázce;
- zda jeho výsledek dokážeme smysluplně interpretovat.
:::

::: {.fragment .takeaway-strip}
Žádné jediné číslo neodpoví na všechny tyto otázky.
:::

::::

# Jak změřit shodu modelu s daty?

## Co už víme o součtu čtverců residuí?

:::: {.layout-stack .layout-center-y}

::: columns
::: {.column width="32%"}
::: {.card-data}
[Hledání přímky]{.card-title}

Kandidátní přímky jsme porovnávali [součtem čtverců residuí]{.text-hgl-orange}.
:::
:::
::: {.column width="32%"}
::: {.card-data}
[Zkratka RSS]{.card-title}

Stejnou veličinu jsme označili **RSS** (*residual sum of squares*).
:::
:::
::: {.column width="32%"}
::: {.card-data}
[Skupinové průměry]{.card-title}

Skupinové průměry RSS zmenšily oproti jednomu společnému průměru.
:::
:::
:::

<br>

::: {.panel-question}
Co se stane s RSS, když přímka prochází blíže pozorovaným bodům?
:::

::: {.fragment}
Kratší residua → menší čtverce → [menší RSS]{.text-hgl-orange} = těsnější shoda s daty.
:::

::::

## Je přímka lepší než předpověď průměrem?

```{r}
#| label: pripravit-sest-stromu
#| include: false
vec_radky_sest <- c(1, 11, 22, 32, 43, 58)

data_sest <-
  data_stromy[
    vec_radky_sest,
    c("druh_stromu", "prumer_kmene_cm", "vyska_stromu_m")
  ]

data_sest$strom <- LETTERS[seq_len(nrow(data_sest))]

prumerna_vyska_sest <- mean(data_sest$vyska_stromu_m)

graf_sest_zaklad <-
  ggplot2::ggplot(
    data = data_sest,
    mapping = ggplot2::aes(x = prumer_kmene_cm, y = vyska_stromu_m)
  ) +
  ggplot2::geom_hline(
    yintercept = prumerna_vyska_sest,
    linetype = "dashed",
    colour = biostat_cols[["graphite"]],
    linewidth = 1
  ) +
  ggplot2::annotate(
    geom = "text",
    x = min(data_sest$prumer_kmene_cm),
    y = prumerna_vyska_sest,
    label = paste0("průměr ", format_value_cz(prumerna_vyska_sest, digits = 1), " m"),
    vjust = -0.6,
    hjust = 0,
    colour = biostat_cols[["graphite"]],
    size = 5.5
  ) +
  ggplot2::geom_point(size = 4.5, colour = biostat_cols[["grey_olive"]]) +
  # Letters end just left of each point so they never cover it
  ggplot2::geom_text(
    mapping = ggplot2::aes(label = strom),
    nudge_x = -1.4,
    hjust = 1,
    size = 6,
    fontface = "bold",
    colour = biostat_cols[["graphite"]]
  ) +
  ggplot2::coord_cartesian(
    xlim = range(data_sest$prumer_kmene_cm) + c(-5, 3),
    ylim = range(data_sest$vyska_stromu_m, prumerna_vyska_sest) + c(-3, 3)
  ) +
  ggplot2::labs(x = "Průměr kmene (cm)", y = "Výška stromu (m)") +
  theme_l09
```

::: columns
::: {.column width="62%"}
```{r}
#| label: zobrazit-sest-stromu
#| echo: false
#| fig-width: 8
#| fig-height: 6
#| fig-alt: "Šest skutečných stromů A až F a vodorovná přerušovaná čára jejich průměrné výšky."
graf_sest_zaklad
```
:::
::: {.column width="36%"}

<br>

Šest skutečných stromů z tabulky.

Předpověď bez prediktoru: každému stromu [průměrnou výšku]{.text-hgl-graphite}.

<br>

::: {.panel-question}
Prošla by přímka podle průměru kmene body blíže než průměr?
:::
:::
:::

## Jak velké jsou rozdíly výšek kolem průměru?

```{r}
#| label: spocitat-celkovy-soucet-sest
#| include: false
data_sest$odchylka <- data_sest$vyska_stromu_m - prumerna_vyska_sest
data_sest$ctverec_odchylky <- data_sest$odchylka^2
ss_total_sest <- sum(data_sest$ctverec_odchylky)
```

::: columns
::: {.column width="56%"}
```{r}
#| label: zobrazit-odchylky-sest
#| echo: false
#| fig-width: 7.5
#| fig-height: 6
#| fig-alt: "Stejných šest stromů se svislými úsečkami od průměrné výšky k pozorované výšce."
graf_sest_zaklad +
  ggplot2::geom_segment(
    mapping = ggplot2::aes(
      xend = prumer_kmene_cm,
      y = prumerna_vyska_sest,
      yend = vyska_stromu_m
    ),
    colour = biostat_cols[["graphite"]],
    linewidth = 1.3
  )
```
:::
::: {.column width="42%" .text-size-smaller}
```{r}
#| label: tabulka-odchylky-sest
#| echo: false
knitr::kable(
  x = data.frame(
    Strom = data_sest$strom,
    Odchylka = format_value_cz(data_sest$odchylka, digits = 1),
    `Čtverec` = format_value_cz(data_sest$ctverec_odchylky, digits = 1),
    check.names = FALSE
  ),
  align = "crr"
)
```

::: {.fragment}
Součet čtverců: [`r format_value_cz(ss_total_sest, digits = 1)` m²]{.text-hgl-graphite}
:::
:::
:::

::: {.fragment .takeaway-strip}
Tomuto součtu říkáme [celkový součet čtverců]{.text-hgl-graphite} (*total sum of squares*).
:::

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

```{r}
#| label: spocitat-rss-sest
#| include: false
mod_sest <- lm(vyska_stromu_m ~ prumer_kmene_cm, data = data_sest)
data_sest$odhad <- fitted(object = mod_sest)
data_sest$residuum <- residuals(object = mod_sest)
data_sest$ctverec_residua <- data_sest$residuum^2
ss_res_sest <- sum(data_sest$ctverec_residua)

strom_max_residuum <- data_sest$strom[which.max(data_sest$residuum)]
strom_max_ctverec <- data_sest$strom[which.max(data_sest$ctverec_residua)]
```

::: columns
::: {.column width="56%"}
```{r}
#| label: zobrazit-residua-sest
#| echo: false
#| fig-width: 7.5
#| fig-height: 6
#| fig-alt: "Stejných šest stromů s fialovou odhadnutou přímkou a oranžovými úsečkami residuí."
graf_sest_zaklad +
  ggplot2::geom_abline(
    intercept = coef(mod_sest)[[1]],
    slope = coef(mod_sest)[[2]],
    colour = biostat_cols[["indigo_velvet"]],
    linewidth = 1.4
  ) +
  ggplot2::geom_segment(
    data = data_sest,
    mapping = ggplot2::aes(
      xend = prumer_kmene_cm,
      y = odhad,
      yend = vyska_stromu_m
    ),
    colour = biostat_cols[["orange"]],
    linewidth = 1.6
  )
```
:::
::: {.column width="42%" .text-size-smaller}
```{r}
#| label: tabulka-residua-sest
#| echo: false
knitr::kable(
  x = data.frame(
    Strom = data_sest$strom,
    Residuum = format_value_cz(data_sest$residuum, digits = 1),
    `Čtverec` = format_value_cz(data_sest$ctverec_residua, digits = 1),
    check.names = FALSE
  ),
  align = "crr"
)
```

RSS: [`r format_value_cz(ss_res_sest, digits = 1)` m²]{.text-hgl-orange} · celkový: `r format_value_cz(ss_total_sest, digits = 1)` m²
:::
:::

::: {.fragment}
Který strom leží [nejvýše nad přímkou]{.text-hgl-orange}? A který [nejvíce přispívá k RSS]{.text-hgl-orange}? [Strom `r strom_max_residuum`; k RSS nejvíce `r strom_max_ctverec` (rozhoduje velikost, ne znaménko).]{.fragment}
:::

::: {.notes}
Odpověď: nejvýše nad přímkou leží strom `r strom_max_residuum` (největší kladné residuum), nejvíce k RSS přispívá strom `r strom_max_ctverec` (největší čtverec). Umocnění odstraní znaménko — o příspěvku k RSS rozhoduje velikost residua, ne jeho směr.
:::

## Jakou část rozdílů přímka zachytila?

```{r}
#| label: spocitat-r2-sest
#| include: false
r2_sest <- 1 - ss_res_sest / ss_total_sest
```

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

::: {.panel-question}
Jaká část rozdílů ve výšce [zůstala]{.text-hgl-orange} a jakou část přímka [zachytila]{.text-hgl-indigo_velvet}? Dosaďte RSS `r format_value_cz(ss_res_sest, digits = 1)` m² a celkový součet čtverců `r format_value_cz(ss_total_sest, digits = 1)` m².
:::

::: {.fragment}
$$
\text{zachycený podíl} = 1 - \frac{`r format_value_cz(ss_res_sest, digits = 1, math = TRUE)`}{`r format_value_cz(ss_total_sest, digits = 1, math = TRUE)`}
$$
:::

::: {.fragment}
$$
\text{zachycený podíl} = `r format_value_cz(r2_sest, digits = 3, math = TRUE)`
$$
:::

## Zachycený podíl se jmenuje R²

:::: {.layout-stack .layout-center-y}

::: {.panel-result}
[Koeficient determinace R² (*coefficient of determination*)]{.panel-title}

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

::: incremental
- $SS_{\mathrm{res}}$ — [součet čtverců residuí]{.text-hgl-orange} (RSS) od odhadnuté přímky;
- $SS_{\mathrm{total}}$ — [celkový součet čtverců]{.text-hgl-graphite} odchylek od průměru;
- u šesti stromů: $R^2 = `r format_value_cz(r2_sest, digits = 3, math = TRUE)`$.
:::

::::

## Jak se mění RSS a R², když přímku natáčíme?

```{r}
#| label: pripravit-natoceni-primky
#| include: false
# Every candidate line passes through the centre of the six trees
# (mean diameter, mean height). Slope 0 is the mean-only prediction,
# the fitted slope is the least-squares line from the previous slides.
prumer_x_sest <- mean(data_sest$prumer_kmene_cm)
sklon_odhad_sest <- coef(mod_sest)[["prumer_kmene_cm"]]

spocitat_natoceni <- function(sklon) {
  odhad <- prumerna_vyska_sest + sklon * (data_sest$prumer_kmene_cm - prumer_x_sest)
  rss <- sum((data_sest$vyska_stromu_m - odhad)^2)
  list(
    sklon = sklon,
    odhad = odhad,
    rss = rss,
    r2 = 1 - rss / ss_total_sest
  )
}

stopifnot(
  isTRUE(all.equal(spocitat_natoceni(0)$rss, ss_total_sest)),
  isTRUE(all.equal(spocitat_natoceni(sklon_odhad_sest)$rss, ss_res_sest))
)

nakreslit_natoceni <- function(sklon, velikost_textu = 1) {
  natoceni <- spocitat_natoceni(sklon)
  data_usecky <- data.frame(
    prumer_kmene_cm = data_sest$prumer_kmene_cm,
    vyska_stromu_m = data_sest$vyska_stromu_m,
    odhad = natoceni$odhad
  )

  graf_body <-
    graf_sest_zaklad +
    ggplot2::geom_abline(
      intercept = prumerna_vyska_sest - sklon * prumer_x_sest,
      slope = sklon,
      colour = biostat_cols[["indigo_velvet"]],
      linewidth = 1.4
    ) +
    ggplot2::geom_segment(
      data = data_usecky,
      mapping = ggplot2::aes(
        x = prumer_kmene_cm,
        xend = prumer_kmene_cm,
        y = odhad,
        yend = vyska_stromu_m
      ),
      colour = biostat_cols[["orange"]],
      linewidth = 1.6
    ) +
    ggplot2::labs(
      subtitle = paste0("sklon přímky: ", format_value_cz(sklon, digits = 2), " m/cm")
    )

  data_sloupce <- data.frame(
    soucet = factor(
      c("celkový", "RSS"),
      levels = c("RSS", "celkový")
    ),
    hodnota = c(ss_total_sest, natoceni$rss)
  )

  graf_sloupce <-
    ggplot2::ggplot(
      data = data_sloupce,
      mapping = ggplot2::aes(x = hodnota, y = soucet, fill = soucet)
    ) +
    ggplot2::geom_col(width = 0.6) +
    ggplot2::geom_text(
      mapping = ggplot2::aes(label = format_value_cz(hodnota, digits = 0)),
      hjust = -0.15,
      size = 5.5 * velikost_textu,
      colour = biostat_cols[["graphite"]]
    ) +
    ggplot2::scale_fill_manual(
      values = c(
        "celkový" = biostat_cols[["graphite"]],
        "RSS" = biostat_cols[["orange"]]
      ),
      guide = "none"
    ) +
    ggplot2::scale_x_continuous(limits = c(0, ss_total_sest * 1.35)) +
    ggplot2::labs(
      x = "Součet čtverců (m²)",
      y = NULL,
      title = paste0("R² = ", format_value_cz(natoceni$r2, digits = 3))
    ) +
    theme_l09 +
    ggplot2::theme(
      plot.title = ggplot2::element_text(
        colour = biostat_cols[["indigo_velvet"]],
        face = "bold",
        size = 26 * velikost_textu
      ),
      panel.grid.major.y = ggplot2::element_blank(),
      axis.text.x = ggplot2::element_blank()
    )

  cowplot::plot_grid(
    graf_body,
    graf_sloupce,
    nrow = 1,
    rel_widths = c(1.45, 1)
  )
}

# Smooth back-and-forth sweep with a pause at both ends
vec_krok <- (1 - cos(seq(from = 0, to = pi, length.out = 28))) / 2
vec_sklony_animace <-
  sklon_odhad_sest *
  c(rep(0, 8), vec_krok, rep(1, 8), rev(vec_krok))

adresar_snimku <- file.path(tempdir(), "natoceni_primky")
dir.create(adresar_snimku, showWarnings = FALSE)

vec_soubory_snimku <-
  purrr::imap_chr(
    .x = vec_sklony_animace,
    .f = function(sklon, poradi) {
      soubor <- file.path(adresar_snimku, sprintf("snimek_%03d.png", poradi))
      ggplot2::ggsave(
        filename = soubor,
        plot = nakreslit_natoceni(sklon),
        width = 12,
        height = 5.6,
        dpi = 110,
        bg = biostat_cols[["parchment"]]
      )
      soubor
    }
  )

magick::image_read(vec_soubory_snimku) |>
  magick::image_animate(delay = 8, optimize = TRUE) |>
  magick::image_write(
    path = here::here(path_materials, "natoceni_primky_r2.gif")
  )
```

```{r}
#| label: zobrazit-natoceni-primky
#| echo: false
#| fig-alt: "Animace: šest stromů a fialová přímka, která se natáčí kolem středu dat od vodorovného průměru k odhadnuté přímce a zpět. Oranžové úsečky residuí se zkracují, sloupec RSS se zmenšuje vůči celkovému součtu čtverců a R² roste od 0 k hodnotě odhadnuté přímky."
include_local_figure("natoceni_primky_r2.gif")
```

::: {.panel-question}
Sledujte oranžové úsečky a sloupec RSS: kdy je R² nulové a kdy největší?
:::

::: {.notes}
Všechny přímky procházejí středem dat (průměrný průměr kmene, průměrná výška). Vodorovná přímka je předpověď průměrem, krajní poloha je odhadnutá přímka. Hodnota 1 − RSS / celkový součet čtverců se počítá pro každou polohu; R² modelu je ta pro odhadnutou přímku.
:::

## Stejných šest stromů, čtyři natočení přímky

```{r}
#| label: zobrazit-natoceni-staticky
#| echo: false
#| fig-height: 5.2
#| fig-alt: "Čtyři panely se stejnými šesti stromy: přímka natočená na 0, třetinu, dvě třetiny a celý odhadnutý sklon. Residua se zkracují, RSS klesá a R² roste od 0 k hodnotě odhadnuté přímky."
vec_podily_sklonu <- c(0, 1 / 3, 2 / 3, 1)

list_natoceni_panely <-
  purrr::map(
    .x = vec_podily_sklonu * sklon_odhad_sest,
    .f = function(sklon) {
      natoceni <- spocitat_natoceni(sklon)
      # Lean base for small panels: same points and axes, no text labels
      ggplot2::ggplot(
        data = data_sest,
        mapping = ggplot2::aes(x = prumer_kmene_cm, y = vyska_stromu_m)
      ) +
        ggplot2::geom_hline(
          yintercept = prumerna_vyska_sest,
          linetype = "dashed",
          colour = biostat_cols[["graphite"]],
          linewidth = 0.8
        ) +
        ggplot2::geom_point(size = 3.5, colour = biostat_cols[["grey_olive"]]) +
        ggplot2::coord_cartesian(
          xlim = range(data_sest$prumer_kmene_cm) + c(-3, 3),
          ylim = range(data_sest$vyska_stromu_m, prumerna_vyska_sest) + c(-3, 3)
        ) +
        theme_l09 +
        ggplot2::geom_abline(
          intercept = prumerna_vyska_sest - sklon * prumer_x_sest,
          slope = sklon,
          colour = biostat_cols[["indigo_velvet"]],
          linewidth = 1.2
        ) +
        ggplot2::geom_segment(
          data = data.frame(
            prumer_kmene_cm = data_sest$prumer_kmene_cm,
            vyska_stromu_m = data_sest$vyska_stromu_m,
            odhad = natoceni$odhad
          ),
          mapping = ggplot2::aes(
            x = prumer_kmene_cm,
            xend = prumer_kmene_cm,
            y = odhad,
            yend = vyska_stromu_m
          ),
          colour = biostat_cols[["orange"]],
          linewidth = 1.3
        ) +
        ggplot2::labs(
          x = NULL,
          y = NULL,
          title = paste0("RSS ", format_value_cz(natoceni$rss, digits = 0), " · R² ", format_value_cz(natoceni$r2, digits = 2))
        ) +
        ggplot2::theme(
          plot.title = ggplot2::element_text(size = 15, face = "bold"),
          axis.text = ggplot2::element_text(size = 11)
        )
    }
  )

cowplot::plot_grid(plotlist = list_natoceni_panely, nrow = 1)
```

::: {.fragment}
Vodorovná přímka = předpověď [průměrem]{.text-hgl-graphite}: R² = 0. Odhadnutá přímka: [nejmenší RSS]{.text-hgl-orange}, a proto [největší R²]{.text-hgl-indigo_velvet}.
:::

## Dává `summary()` stejné R² pro všech `r nrow(data_stromy)` stromů?

```{r}
#| label: spocitat-r2-vsechny
#| include: false
ss_total_vsechny <- sum(
  (data_stromy$vyska_stromu_m - mean(data_stromy$vyska_stromu_m))^2
)
ss_res_spolecny <- deviance(object = mod_spolecny)
r2_spolecny <- 1 - ss_res_spolecny / ss_total_vsechny

stopifnot(
  isTRUE(all.equal(r2_spolecny, summary(mod_spolecny)$r.squared))
)
```

:::: {.layout-stack .layout-center-y}

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

::: {.fragment}
```{r}
#| label: r2-ze-summary
#| echo: true
summary(object = mod_spolecny)$r.squared
```
:::

::: {.fragment}
`Multiple R-squared` ve výpisu je [tentýž výpočet]{.text-hgl-indigo_velvet} — nic nového.
:::

::::

##

:::: {.slide-margin-top-15}

::: {.box-question}
[Společná přímka má R² = `r format_value_cz(r2_spolecny, digits = 3)`. Co to znamená?]{.text-size-heading .text-bold}
:::

::: columns
::: {.column width="32%"}
::: {.card-question}
**A**

Model správně předpověděl `r format_value_cz(100 * r2_spolecny, digits = 1)` % stromů.
:::
:::
::: {.column width="32%"}
::: {.card-question}
**B**

[Model zachytil `r format_value_cz(100 * r2_spolecny, digits = 1)` % variability výšky kolem průměru v těchto datech.]{.rn-circle-orange}
:::
:::
::: {.column width="32%"}
::: {.card-question}
**C**

S pravděpodobností `r format_value_cz(r2_spolecny, digits = 3)` je model pravdivý.
:::
:::
:::

::: {.panel-olive}
Přesvědčte souseda, proč jste vybrali svoji možnost.
:::

::: {.fragment}
Obhajitelná je [B]{.text-hgl-indigo_velvet}: R² popisuje podíl zachycené variability v datech, ze kterých byl model odhadnut.
:::

::::

## Co R² popisuje a na co nestačí?

:::: {.layout-stack .layout-center-y}

::: columns
::: {.column width="48%"}
::: {.panel-result}
[R² popisuje]{.panel-title}

- podíl [variability odezvy]{.text-hgl-indigo_velvet} zachycený modelem;
- ve [stejných datech]{.text-hgl-graphite}, na kterých byl model odhadnut.
:::
:::
::: {.column width="48%"}
::: {.panel-question}
[R² neříká]{.panel-title}

- zda model splňuje předpoklady;
- zda je vztah příčinný;
- jak dobře model poslouží pro nové stromy;
- [zda je to nejlepší z našich kandidátů?]{.text-hgl-orange}
:::
:::
:::

::::

# Pomáhá každý koeficient navíc?

## Co udělá R² po přidání druhu a interakce?

:::: {.layout-stack .layout-center-y}

[`prumer` → `prumer + druh` → `prumer * druh`]{.text-size-heading}

[Co se stane s R², když přidáme druh a interakci?]{.text-bold}

::: columns
::: {.column width="32%"}
::: {.card-question}
**A** · R² poroste
:::
:::
::: {.column width="32%"}
::: {.card-question}
**B** · R² může klesnout
:::
:::
::: {.column width="32%"}
::: {.card-question}
**C** · R² zůstane stejné
:::
:::
:::

::::

::: {.notes}
Hlasování rukou. Odpověď na dalším snímku; zde ji neprozrazovat.
:::

## R² po přidání členů neklesne

```{r}
#| label: pripravit-porovnani-r2
#| include: false
data_porovnani <- data.frame(
  model = factor(vec_nazvy_modelu, levels = vec_nazvy_modelu),
  k = purrr::map_int(list_kandidati, ~ length(coef(.x))),
  rss = purrr::map_dbl(list_kandidati, deviance),
  r2 = purrr::map_dbl(list_kandidati, ~ summary(.x)$r.squared)
)
```

::: columns
::: {.column width="56%"}
```{r}
#| label: zobrazit-porovnani-r2
#| echo: false
#| fig-width: 7
#| fig-height: 5.4
#| fig-alt: "R² tří kandidátních modelů roste od společné přímky přes posun mezi druhy k různým sklonům."
ggplot2::ggplot(
  data = data_porovnani,
  mapping = ggplot2::aes(x = model, y = r2, group = 1)
) +
  ggplot2::geom_line(colour = biostat_cols[["indigo_velvet"]], linewidth = 1) +
  ggplot2::geom_point(colour = biostat_cols[["indigo_velvet"]], size = 5) +
  ggplot2::geom_text(
    mapping = ggplot2::aes(label = format_value_cz(r2, digits = 3)),
    # First label sits left of its point so the rising line does not cross it
    hjust = c(1.25, 0.5, 0.5),
    vjust = c(0.5, -1.1, -1.1),
    size = 6,
    colour = biostat_cols[["indigo_velvet"]]
  ) +
  ggplot2::scale_x_discrete(labels = function(x) stringr::str_wrap(x, width = 10)) +
  ggplot2::scale_y_continuous(
    limits = c(min(data_porovnani$r2) - 0.01, max(data_porovnani$r2) + 0.012),
    labels = function(x) format_value_cz(x, digits = 2)
  ) +
  ggplot2::labs(x = NULL, y = "R²") +
  theme_l09
```
:::
::: {.column width="42%"}

<br>

```{r}
#| label: tabulka-porovnani-r2
#| echo: false
knitr::kable(
  x = data.frame(
    Model = data_porovnani$model,
    k = data_porovnani$k,
    RSS = format_value_cz(data_porovnani$rss, digits = 0)
  ),
  align = "lrr"
)
```

::: {.fragment}
<br>

Složitější model umí [zopakovat jednodušší]{.text-hgl-indigo_velvet}, a tak R² nejhůř zůstane stejné; klesnout nemůže.
:::
:::
:::

::: {.fragment .box-warning}
**Vyšší R² není automatické vítězství.** I téměř zbytečný prediktor R² nesníží.
:::

## Které koeficienty jednotlivé modely odhadují?

```{r}
#| label: pripravit-nazvy-koeficientu
#| include: false
list_nazvy_koeficientu <- purrr::map(list_kandidati, ~ names(coef(.x)))

# Non-breaking spaces keep each coefficient name (with the species name)
#   on one line instead of wrapping in the middle
text_koeficienty <- purrr::map_chr(
  list_nazvy_koeficientu,
  ~ paste0("`", gsub(" ", " ", .x, fixed = TRUE), "`", collapse = "<br>")
)
```

::: {.panel-question}
Který seznam patří společné přímce, který posunu mezi druhy a který různým sklonům?
:::

::: columns
::: {.column width="56%"}
::: {.card-indigo .text-size-smaller}
[A[ · Posun mezi druhy · k = `r length(list_nazvy_koeficientu[[2]])`]{.fragment}]{.card-title}

`r text_koeficienty[2]`
:::
:::
::: {.column width="42%"}
::: {.card-indigo .text-size-smaller}
[B[ · Společná přímka · k = `r length(list_nazvy_koeficientu[[1]])`]{.fragment}]{.card-title}

`r text_koeficienty[1]`
:::
:::
:::

::: {.card-indigo .text-size-smaller}
[C[ · Různé sklony · k = `r length(list_nazvy_koeficientu[[3]])`]{.fragment}]{.card-title}

`r text_koeficienty[3]`
:::

## Jak R² opravit za počet koeficientů?

```{r}
#| label: spocitat-adjustovane-r2-aditivni
#| include: false
n_stromu <- nobs(object = mod_aditivni)
k_aditivni <- length(coef(object = mod_aditivni))
r2_aditivni <- summary(object = mod_aditivni)$r.squared
r2_adj_aditivni <- summary(object = mod_aditivni)$adj.r.squared

stopifnot(
  isTRUE(all.equal(
    r2_adj_aditivni,
    1 - (1 - r2_aditivni) * (n_stromu - 1) / (n_stromu - k_aditivni)
  ))
)
```

$$
\text{adjustované R}^2 = 1 - \text{nezachycený podíl} \times \text{úprava za počet stromů a koeficientů}
$$

::: {.fragment}
Posun mezi druhy: $R^2 = `r format_value_cz(r2_aditivni, digits = 3, math = TRUE)`$, [$n = `r n_stromu`$]{.text-hgl-graphite} stromů, [$k = `r k_aditivni`$]{.text-hgl-indigo_velvet} koeficienty

$$
R^2_{\mathrm{adj}} = 1 - \left(1 - `r format_value_cz(r2_aditivni, digits = 3, math = TRUE)`\right) \times \frac{`r n_stromu` - 1}{`r n_stromu` - `r k_aditivni`} = `r format_value_cz(r2_adj_aditivni, digits = 3, math = TRUE)`
$$
:::

::: {.fragment}
**Adjustované R²** (*adjusted R-squared*):

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

::: {.notes}
Ukažte, kde je shoda (1 − R²), velikost dat (n) a složitost (k). S každým koeficientem navíc roste zlomek (n − 1)/(n − k).
:::

## Co udělá adjustované R² po přidání interakce?

:::: {.layout-stack .layout-center-y}

Posun mezi druhy → Různé sklony: R² [`r format_value_cz(data_porovnani$r2[2], digits = 3)` → `r format_value_cz(data_porovnani$r2[3], digits = 3)`]{.text-hgl-indigo_velvet} · *k* `r data_porovnani$k[2]` → `r data_porovnani$k[3]` · *n* = `r n_stromu`

[**A** · adjustované R² vzroste]{.text-size-heading}

[**B** · adjustované R² klesne]{.text-size-heading}

::::

## Adjustované R² po přidání interakce klesne

```{r}
#| label: pripravit-adjustovane-r2-kandidatu
#| include: false
data_porovnani$r2_adj <- purrr::map_dbl(list_kandidati, ~ summary(.x)$adj.r.squared)
```

::: columns
::: {.column width="58%"}
```{r}
#| label: zobrazit-r2-a-adjustovane
#| echo: false
#| fig-width: 7.5
#| fig-height: 5.4
#| fig-alt: "R² tří modelů stále roste, zatímco adjustované R² je nejvyšší u posunu mezi druhy a u různých sklonů klesne."
data_r2_dlouhe <- data.frame(
  model = rep(data_porovnani$model, times = 2),
  mira = factor(
    rep(c("R²", "adjustované R²"), each = nrow(data_porovnani)),
    levels = c("R²", "adjustované R²")
  ),
  hodnota = c(data_porovnani$r2, data_porovnani$r2_adj)
)

ggplot2::ggplot(
  data = data_r2_dlouhe,
  mapping = ggplot2::aes(x = model, y = hodnota, group = mira, linetype = mira)
) +
  ggplot2::geom_line(colour = biostat_cols[["indigo_velvet"]], linewidth = 1) +
  ggplot2::geom_point(colour = biostat_cols[["indigo_velvet"]], size = 4.5) +
  ggplot2::geom_text(
    mapping = ggplot2::aes(
      label = format_value_cz(hodnota, digits = 3),
      vjust = ifelse(mira == "R²", -1, 2)
    ),
    size = 5.5,
    colour = biostat_cols[["indigo_velvet"]]
  ) +
  ggplot2::scale_linetype_manual(values = c("solid", "dashed")) +
  ggplot2::scale_x_discrete(labels = function(x) stringr::str_wrap(x, width = 10)) +
  ggplot2::scale_y_continuous(
    limits = range(data_r2_dlouhe$hodnota) + c(-0.008, 0.01),
    labels = function(x) format_value_cz(x, digits = 2)
  ) +
  ggplot2::labs(x = NULL, y = NULL, linetype = NULL) +
  theme_l09
```
:::
::: {.column width="40%"}

<br>

Dva sklony navíc zlepšily R² jen o [`r format_value_cz(data_porovnani$r2[3] - data_porovnani$r2[2], digits = 3)`]{.text-hgl-indigo_velvet}.

<br>

::: {.fragment .panel-result}
Zlomek $\frac{n-1}{n-k}$ vzrostl víc, než se zmenšil nezachycený podíl → [úprava zvítězila]{.text-hgl-orange}.
:::
:::
:::

## Co už o třech modelech víme?

:::: {.layout-stack .layout-center-y}

::: {.box-question}
Jak poznáme, že je model dobrý?
:::

::: columns
::: {.column width="48%"}
::: {.panel-result}
[Už umíme říct]{.panel-title}

- R² měří [shodu s daty]{.text-hgl-indigo_velvet} a vždy upřednostní nejpružnější model.
- Adjustované R² upřednostní [posun mezi druhy]{.text-hgl-indigo_velvet}.
:::
:::
::: {.column width="48%" .fragment}
::: {.panel-question}
[Zatím nevíme]{.panel-title}

- Kdy se složitost [obecně vyplatí]{.text-hgl-orange}?
- Vypadají residua vybraného modelu v pořádku?
:::
:::
:::

::::

# Kdy se složitost vyplatí?

## Vyplatí se model šitý na míru?

:::: {.layout-stack .layout-center-y}

::: columns
::: {.column width="55%"}
```{r}
#| label: zobrazit-stromy-krejci
#| echo: false
#| out-width: "100%"
#| fig-alt: "Liška krejčí šije jehličnanu oblek s desítkami švů a špendlíků, které kopírují každý hrbol kmene; vedle visí jednoduchý dobře padnoucí kabát a za nimi čeká další strom."
include_optional_figure(data_source = "stromy_krejci.png")
```

::: {.text-size-tiny .text-right}
`r ai_disclosure("stromy_krejci.png")`
:::
:::
::: {.column width="43%"}
::: {.panel-question}
**Ve dvojici:** Posun mezi druhy a různé sklony popisují stromy téměř stejně dobře (R² `r format_value_cz(summary(mod_aditivni)$r.squared, digits = 3)` a `r format_value_cz(summary(mod_interakce)$r.squared, digits = 3)`). První potřebuje [`r length(coef(mod_aditivni))`]{.text-hgl-indigo_velvet} koeficienty, druhý [`r length(coef(mod_interakce))`]{.text-hgl-indigo_velvet}.

Který byste upřednostnili? Co byste chtěli vědět, než se rozhodnete?
:::
:::
:::

::: {.fragment .panel-result}
[Princip úspornosti (*parsimony*)]{.panel-title}

Složitější model upřednostníme jen tehdy, když [dostatečně zlepší shodu]{.text-hgl-orange} s daty.
:::

::::

::: {.notes}
Oblek s desítkami švů sedí jen na tento jeden strom; jednoduchý kabát padne dobře i dalším stromům. Obraz nerozhoduje, který model je lepší; ptá se, zda se šití navíc vyplatí.
:::

## AIC váží shodu i složitost

:::: {.slide-margin-top-15}

::: columns
::: {.column width="52%"}

::: {.panel-result}
[Akaikeho informační kritérium (*Akaike information criterion*, AIC)]{.panel-title}

- lepší shoda AIC [snižuje]{.text-hgl-indigo_velvet};
- každá odhadovaná veličina navíc AIC [zvyšuje]{.text-hgl-orange};
- preferujeme [nižší AIC]{.text-hgl-indigo_velvet}.
:::
:::
::: {.column width="44%"}

```{r}
#| label: aic-spolecneho-modelu
#| echo: true
AIC(object = mod_spolecny)
```

::: {.fragment}
AIC [není podíl]{.text-hgl-orange} a nemá stejnou škálu jako R².
:::
:::
:::

::::

## Model má AIC `r format_value_cz(AIC(mod_spolecny), digits = 2)`. Je to dobrý model?

:::: {.slide-margin-top-25}

::: columns
::: {.column width="32%"}
::: {.card-question}
**A** · Ano, AIC je nízké.
:::
:::
::: {.column width="32%"}
::: {.card-question}
**B** · Ne, AIC je vysoké.
:::
:::
::: {.column width="32%"}
::: {.card-question}
**C** · [Z jedné hodnoty to nepoznáme.]{.rn-circle-orange}
:::
:::
:::

::: {.fragment .box-rule}
AIC je [relativní]{.text-hgl-indigo_velvet}. Porovnáváme jen modely se [stejnou odezvou]{.text-hgl-graphite}, [stejnými pozorováními]{.text-hgl-graphite} a stejným předpokladem o rozdělení residuí.
:::

::::

::: {.notes}
Hlasování a krátké přesvědčování souseda proběhne před odkrytím; kroužek u C a pravidlo přijdou až po něm.
:::

## O kolik je AIC každého kandidáta vyšší než nejnižší?

```{r}
#| label: spocitat-delta-aic
#| include: false
data_porovnani$aic <- purrr::map_dbl(list_kandidati, AIC)
aic_min <- min(data_porovnani$aic)
data_porovnani$delta_aic <- data_porovnani$aic - aic_min
model_min_aic <- as.character(data_porovnani$model[which.min(data_porovnani$aic)])
```

$$
\Delta\text{AIC} = \text{AIC kandidáta} - \text{nejnižší AIC v sadě}
$$

::: {.fragment}

Společná přímka: $`r format_value_cz(data_porovnani$aic[1], digits = 2, math = TRUE)` - `r format_value_cz(aic_min, digits = 2, math = TRUE)` = `r format_value_cz(data_porovnani$delta_aic[1], digits = 2, math = TRUE)`$
:::

::: {.fragment}
```{r}
#| label: tabulka-delta-aic
#| echo: false
knitr::kable(
  x = data.frame(
    Model = data_porovnani$model,
    AIC = format_value_cz(data_porovnani$aic, digits = 2),
    `ΔAIC` = format_value_cz(data_porovnani$delta_aic, digits = 2),
    check.names = FALSE
  ),
  align = "lrr"
)
```
:::

::: {.fragment}
[`r model_min_aic`]{.text-hgl-indigo_velvet} má ΔAIC = 0 · univerzální práh [nepoužíváme]{.text-hgl-orange}
:::

## Ukazují R², adjustované R² a AIC na stejný model?

```{r}
#| label: tabulka-syntezy
#| echo: false
knitr::kable(
  x = data.frame(
    Model = data_porovnani$model,
    k = data_porovnani$k,
    `R² ↑` = format_value_cz(data_porovnani$r2, digits = 3),
    `adj. R² ↑` = format_value_cz(data_porovnani$r2_adj, digits = 3),
    `AIC ↓` = format_value_cz(data_porovnani$aic, digits = 2),
    check.names = FALSE
  ),
  align = "lrrrr"
)
```

<br>

::: {.panel-question}
**Ve dvojici:** Proč R² upřednostní jiný model než adjustované R² a AIC?
:::

::: {.fragment}
R² hodnotí [jen shodu]{.text-hgl-indigo_velvet}. Adjustované R² a AIC [penalizují složitost]{.text-hgl-orange} — dva sklony navíc se nevyplatily. Nejsou to tři hlasy o pravdě.
:::

## Nejlépe podporovaný kandidát nemusí být dobrý model

:::: {.layout-stack .layout-center-y}

::: columns
::: {.column width="50%"}
```{r}
#| label: zobrazit-stromy-porota
#| echo: false
#| out-width: "100%"
#| fig-alt: "Porota lesních zvířat s prázdnými hodnoticími kartičkami posuzuje na pódiu tři kandidáty, každého s jiným uspořádáním rovných tyček; za pootevřenými dveřmi čeká čtvrtý kandidát s prohnutou tyčkou."
include_optional_figure(data_source = "stromy_porota.png")
```

::: {.text-size-tiny .text-right}
`r ai_disclosure("stromy_porota.png")`
:::
:::
::: {.column width="48%"}
::: {.box-warning}
Porovnání řadí [jen modely, které jsme do sady zařadili]{.text-hgl-orange}.

Neříká, že je vítěz obecně dobrý. Model, který jsme vůbec neposuzovali, může stromy popsat lépe.
:::
:::
:::

::: {.fragment}
[Proto vybraný model ještě zkontrolujeme.]{.text-size-heading .text-bold}
:::

::::

::: {.notes}
Porota hodnotí jen kandidáty, které jsme pozvali na pódium. Čtvrtý kandidát za dveřmi s prohnutou tyčkou připomíná, že lepší popis dat může ležet mimo naši sadu.
:::

# Obstojí vybraný model?

## Zůstává v residuích modelu s posunem mezi druhy nějaký vzorec?

```{r}
#| label: pripravit-diagnostiku-aditivniho-modelu
#| include: false
data_diagnostika <- data.frame(
  odhad = fitted(object = mod_aditivni),
  residuum = residuals(object = mod_aditivni),
  druh_stromu = data_stromy$druh_stromu
)

graf_residua_aditivni <-
  ggplot2::ggplot(
    data = data_diagnostika,
    mapping = ggplot2::aes(x = odhad, y = residuum)
  ) +
  ggplot2::geom_hline(yintercept = 0, linetype = "dashed", colour = biostat_cols[["graphite"]]) +
  ggplot2::geom_point(
    mapping = ggplot2::aes(colour = druh_stromu, shape = druh_stromu),
    size = 2.8,
    alpha = 0.85
  ) +
  ggplot2::scale_colour_manual(values = vec_barvy_druhu, guide = "none") +
  ggplot2::scale_shape_manual(values = vec_symboly_druhu, guide = "none") +
  ggplot2::labs(x = "Odhadnutá výška (m)", y = "Residuum (m)", title = "Residua") +
  theme_l09

graf_qq_aditivni <-
  ggplot2::ggplot(data = data_diagnostika, mapping = ggplot2::aes(sample = residuum)) +
  ggplot2::stat_qq(colour = biostat_cols[["grey_olive"]], size = 2.8) +
  ggplot2::stat_qq_line(colour = biostat_cols[["indigo_velvet"]], linewidth = 1.2) +
  ggplot2::labs(x = "Teoretické kvantily", y = "Kvantily residuí", title = "Q–Q graf") +
  theme_l09
```

::: {.r-stack}
```{r}
#| label: zobrazit-diagnostiku-aditivniho-modelu
#| echo: false
#| fig-height: 3.6
#| fig-alt: "Vlevo residua modelu s posunem mezi druhy proti odhadnuté výšce, vpravo Q–Q graf residuí."
cowplot::plot_grid(graf_residua_aditivni, graf_qq_aditivni, nrow = 1)
```

::: {.fragment}
```{r}
#| label: zobrazit-diagnostiku-aditivniho-modelu-s-trendem
#| echo: false
#| fig-height: 3.6
#| fig-alt: "Stejné grafy; v residuích je navíc vyhlazená oranžová křivka, která tvoří oblouk: nízké i nejvyšší odhadnuté výšky mají residua převážně pod nulou, středně vysoké nad nulou."
cowplot::plot_grid(
  graf_residua_aditivni +
    ggplot2::geom_smooth(
      method = "loess",
      formula = y ~ x,
      se = FALSE,
      colour = biostat_cols[["orange"]],
      linewidth = 1.4
    ),
  graf_qq_aditivni,
  nrow = 1
)
```
:::
:::

[**Ve dvojici:** Leží residua v celém rozsahu odhadnutých výšek náhodně kolem nuly?]{.text-size-smaller}

::: {.fragment .text-size-smaller}
Oranžová křivka jen vyhlazuje průběh bodů, aby byl vzor lépe vidět: residua tvoří [oblouk]{.text-hgl-orange} · Q–Q graf [není přímý]{.text-hgl-orange}.
:::

::: {.notes}
Vyhlazená křivka je tu nová (v předchozích lekcích nebyla); stačí říct, že sleduje průměrnou polohu residuí v každém místě grafu. Nízké i nejvyšší stromy leží převážně pod nulou, středně vysoké nad ní. Oblouk je skutečný a týká se všech tří kandidátů (viz residua společné přímky na začátku); druh stromu ho neodstraní. Nechte otázku otevřenou až do závěrečného snímku.
:::

## Které stromy ovlivňují model nejvíc?

```{r}
#| label: cookova-vzdalenost
#| echo: false
#| fig-height: 4.1
#| fig-alt: "Cookovy vzdálenosti všech stromů pro model Posun mezi druhy s přerušovanou orientační čarou 4/n; několik stromů leží nad ní."
data_cook <- data.frame(
  poradi = seq_len(nobs(mod_aditivni)),
  cook = cooks.distance(model = mod_aditivni)
)
hranice_cook <- 4 / nobs(mod_aditivni)
pocet_nad_hranici <- sum(data_cook$cook > hranice_cook)

ggplot2::ggplot(data = data_cook, mapping = ggplot2::aes(x = poradi, y = cook)) +
  ggplot2::geom_segment(
    mapping = ggplot2::aes(xend = poradi, y = 0, yend = cook),
    colour = biostat_cols[["grey_olive"]],
    linewidth = 1.2
  ) +
  ggplot2::geom_hline(
    yintercept = hranice_cook,
    linetype = "dashed",
    colour = biostat_cols[["orange"]],
    linewidth = 1.1
  ) +
  ggplot2::labs(x = "Pořadí stromu v tabulce", y = "Cookova vzdálenost") +
  theme_l09
```

Cookova vzdálenost: o kolik by se odhadnuté hodnoty změnily, kdyby strom chyběl. Nad orientační čarou $4/n$: [`r pocet_nad_hranici`]{.text-hgl-orange} stromy.

::: {.fragment .takeaway-strip}
Důvod [stromy prověřit]{.text-hgl-orange} — ne je automaticky odstranit.
:::

::: {.notes}
Cookova vzdálenost byla v předchozích lekcích jen v doplňujícím textu, proto ji krátce vysvětlete: jak moc by se odhadnuté hodnoty změnily, kdybychom jeden strom vynechali.
:::

## Obhajte model dvěma větami

:::: {.layout-stack .layout-center-y}

::: {.panel-question}
[40 s sami · ve dvojici porovnejte]{.panel-title}

„Ze tří kandidátů vybíráme model ____, protože ____.“

„Omezením tohoto závěru je ____.“
:::

::: {.fragment}
Opřete se aspoň o [dva různé druhy důkazů]{.text-hgl-orange}: shodu, složitost, diagnostiku, biologický smysl.
:::

::::

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

:::: {.layout-stack .layout-center-y}

::: columns
::: {.column width="50%"}
::: {.panel-result}
[Posun mezi druhy]{.panel-title}

- zachytí rozdíly mezi druhy, které společná přímka přehlíží;
- nejvyšší adjustované R² ([`r format_value_cz(max(data_porovnani$r2_adj), digits = 3)`]{.text-hgl-indigo_velvet}) a nejnižší AIC.
:::
:::
::: {.column width="46%" .fragment}
::: {.panel-question}
[Omezení]{.panel-title}

- vybírali jsme jen ze tří kandidátů;
- residua tvoří [oblouk]{.text-hgl-orange}: přímka nevystihuje tvar vztahu;
- `r nobs(mod_aditivni)` stromů, Q–Q graf není přímý;
- pozorovací data: nejde o příčinný efekt druhu.
:::
:::
:::

::: {.fragment .box-result}
Dobrý model nepoznáme podle jediného čísla. Obhajoba propojuje [otázku, shodu, složitost, diagnostiku a interpretovatelnost]{.text-hgl-indigo_velvet}.
:::

::::

## Jak porovnávat modely bez automatického hledání?

::: columns
::: {.column width="56%"}
::: incremental
1. Položte biologickou otázku.
2. Prohlédněte a zobrazte data.
3. Určete malou sadu [zdůvodněných kandidátů]{.text-hgl-indigo_velvet}.
4. Použijte [stejnou odezvu a stejná pozorování]{.text-hgl-graphite}.
5. Porovnejte přímky, residua, R², adjustované R² a AIC.
6. Zkontrolujte předpoklady a vlivná pozorování.
7. Obhajte závěr a [pojmenujte omezení]{.text-hgl-orange}.
:::
:::
::: {.column width="42%" .fragment}

<br>

::: {.box-warning}
Automatické přidávání a odebírání členů vytvoří mnoho modelů, ale [neurčí, které dávají biologický smysl]{.text-hgl-orange}.
:::
:::
:::

## Co si odnést

:::: {.layout-stack .layout-center-y}

::: incremental
- R² = podíl zachycené variability; počítá se z celkového součtu čtverců a RSS.
- Po přidání členu R² [neklesne]{.text-hgl-indigo_velvet}, a proto samo model nevybírá.
- Adjustované R² a AIC [penalizují složitost]{.text-hgl-orange}, každé jinak; AIC je [relativní]{.text-hgl-indigo_velvet}.
- Kandidáty určujeme [předem]{.text-hgl-orange}: stejná odezva, stejná pozorování.
- Obhajoba spojuje [několik druhů důkazů]{.text-hgl-orange} včetně diagnostiky a omezení.
:::

::::

[Data: J. Marshall, University of Idaho · `lgrdata::allometry` · CC0]{.footer}

## Co když žádná přímka nestačí?

:::: {.layout-stack .layout-center-y}

::: columns
::: {.column width="62%"}
```{r}
#| label: zobrazit-stromy-strop
#| echo: false
#| out-width: "100%"
#| fig-alt: "Řada jehličnanů od mladých po staré: jejich výška nejprve roste a pak se ustálí, zatímco obří rovná fialová tyč stoupá dál až do mraků; liška a veverka se na ni udiveně dívají."
include_optional_figure(data_source = "stromy_strop.png")
```

::: {.text-size-tiny .text-right}
`r ai_disclosure("stromy_strop.png")`
:::
:::
::: {.column width="36%"}
[Všechny tři kandidáty jsme kreslili jako přímky — a residua přesto tvoří oblouk.]{.text-size-heading}
:::
:::

::: {.fragment}
[Jak popsat vztah, který není přímkový?]{.text-size-heading .text-bold .text-hgl-orange}
:::

::::

::: {.notes}
Liška a veverka sledují rovnou fialovou tyč, která stoupá dál, i když výška nejstarších stromů už neroste. Jde o myšlenkový obraz věkové řady, nikoli o výsledek porovnání druhů v datech. Otázku nechte otevřenou: odpověď přinese následující lekce o tvaru vztahu a transformacích.
:::