---
title: "Kdy vztah dvou znaků závisí na třetím?"
subtitle: "Interakce číselného a skupinového znaku, dvou skupinových i dvou číselných znaků"
lang: cs
filters:
  - ../theme/semantic-boxes.lua
format:
  html:
    toc: true
    toc-depth: 3
    toc-title: "Obsah"
    number-sections: true
    code-fold: false
    embed-resources: true
    theme: [flatly, ../theme/skripta_theme.scss]
    highlight-style: github
    include-in-header:
      - "../theme/fonts-include.html"
      - "../theme/skripta-favicon.html"
  typst:
    toc: true
    toc-depth: 3
    number-sections: true
    output-file: "skripta_raw.pdf"
    include-in-header:
      - "skripta_theme.typ"
      - "skripta-logo.typ"
    date: last-modified
date: 2026/09/22
date-format: long
date-modified: last-modified
execute:
  message: false
  warning: false
---

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

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

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

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"))
source(here::here("R", "Functions", "build_tree_bands.R"))
source(here::here("R", "Functions", "plot_tree_panels.R"))
use_ragg_brand_fonts()

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

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

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

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

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

## Liší se vztah délky a zadní šířky krunýře mezi samicemi a samci?

Představte si dva kraby se stejně dlouhým krunýřem, jednoho samce a jednu samici. Musí mít stejně širokou zadní část krunýře? A bude případný rozdíl stejný u kratších i delších krabů? To je otázka **vztahu dvou měřených znaků v různých skupinách**, nikoli otázka, které pohlaví je obecně větší.

Proč bychom rozdíl vůbec čekali? Zadní část krunýře leží zhruba nad zadečkem, který má krab podvinutý pod tělem. Samice pod zadečkem nosí oplozená vajíčka, dokud se z nich nevylíhnou larvy. Campbell a Mahon (1974), z jejichž měření budeme vycházet, zjistili, že se samice a samci pobřežních krabů rodu *Leptograpsus* liší poměrem zadní šířky a délky krunýře. Pokud samice s rostoucí velikostí potřebují pro snůšku relativně více místa, mohl by se rozdíl mezi pohlavími s velikostí krabů měnit. Je to biologické očekávání; zda mu odpovídají i naše data, ukážou až měření.

::: {.box-question title="Hlavní otázka této lekce"}
**Je rozdíl v zadní šířce krunýře mezi stejně dlouhou samicí a samcem stejný u menších i větších krabů?**
:::

Zadní šířku budeme odhadovat z délky krunýře a pohlaví. Data jsou pozorovací: samotný rozdíl mezi skupinami nám neřekne, co jej způsobilo.

## Kdy předpovídají dvě přímky stejný sklon?

V minulé lekci jste viděli model se dvěma `r render_glossary_term("prediktor", display = "prediktory")`: číselným znakem a skupinou. Jeho dvě předpovědní přímky byly rovnoběžné. Jak se podle takového modelu změní rozdíl mezi skupinami, když se číselný prediktor zvětší?

::: {.box-question title="Vybavte si minulou lekci"}
Nakreslete dvě rovnoběžné přímky. Vyberte na ose x dvě různé hodnoty. Je svislá vzdálenost mezi přímkami v obou místech stejná?
:::

::: {.box-result collapse="true" title="Možná odpověď"}
Ano. Rovnoběžné přímky mají společný `r render_glossary_term("sklon", display = "sklon")`, takže jejich svislá vzdálenost zůstává stejná.
:::

Zkusme zjistit, zda takové omezení dává smysl i pro kraby.

## Výsledky učení

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

- rozpoznat biologickou otázku, ve které vztah jednoho `r render_glossary_term("prediktor", display = "prediktoru")` k `r render_glossary_term("odezva", display = "odezvě")` závisí na hodnotě nebo úrovni druhého prediktoru;
- fitovat `r render_glossary_term("linearni-model", display = "lineární model")` s `r render_glossary_term("interakce-predikatoru", display = "interakcí")` a zobrazit jeho předpovědi pro kombinace číselných i skupinových prediktorů;
- vyložit koeficienty modelu s interakcí jako podmíněné a porovnávat skupiny nebo různé hodnoty prediktoru pomocí předpovědí při zvolené hodnotě druhého prediktoru;
- vysvětlit interakci dvou skupinových prediktorů jako rozdíl dvou rozdílů;
- omezit závěr na pozorovaný rozsah dat a na souvislost v pozorovacích datech.

## Jeden řádek tabulky představuje jednoho kraba

Použijeme měření krabů z práce Campbell a Mahon (1974). Autoři změřili stejný počet samic a samců dvou barevných forem pobřežního kraba, modré a oranžové, které považují za samostatné druhy rodu *Leptograpsus*. Soubor [krabi.csv](https://cuni-natur-biostatistics.github.io/L08/current/data/krabi.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. Načteme ho a zkontrolujeme názvy a typy sloupců:

```{r}
#| label: nacist-data
#| echo: true
# Načteme tabulku všech změřených krabů
data_krabi_vse <-
  read.csv(file = "data/krabi.csv")

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

Sloupec `sp` označuje barevnou formu (`B` modrá, `O` oranžová), `sex` pohlaví (`F` samice, `M` samec). `CL` je délka krunýře a `RW` šířka jeho zadní části; obě hodnoty jsou v milimetrech. Délka a zadní šířka jsou dvě různé osy téhož krunýře. Ostatní rozměry v této lekci nepotřebujeme.

Pro první otázku vybereme jen kraby **modré formy** a tři potřebné sloupce. Jeden řádek znamená **jednoho změřeného kraba**, ne opakované měření téhož jedince.

```{r}
#| label: vybrat-modrou-formu
#| echo: true
# Vybereme modrou formu a tři znaky potřebné pro otázku
data_krabi <-
  data_krabi_vse[
    data_krabi_vse$sp == "B",
    c("sex", "CL", "RW")
  ]

# Přejmenujeme sloupce podle toho, co znamenají
names(data_krabi) <-
  c("pohlavi", "delka_krunyre_mm", "zadni_sirka_mm")

# Pohlaví převedeme na faktor; první úroveň bude referenční skupinou
data_krabi$pohlavi <-
  factor(
    x = data_krabi$pohlavi,
    levels = c("F", "M"),
    labels = c("samice", "samec")
  )

# Každý řádek tabulky patří jednomu krabovi
head(x = data_krabi, n = 6)
```

::: {.box-data title="Doplňující: jak vznikl soubor krabi.csv" collapse="true"}
Soubor `krabi.csv` je nezměněná kopie tabulky `crabs` z balíčku MASS, který se instaluje spolu s R. Připravený soubor používáme proto, aby všichni pracovali se stejnou verzí dat. Stejnou tabulku si můžete vytvořit sami; kód spusťte v hlavní složce svého projektu:

```{r}
#| label: extra-priprava-krabi
#| echo: true
#| eval: false
# Tabulka je součástí balíčku MASS
data_krabi_zdroj <-
  as.data.frame(x = MASS::crabs)

# Kontrola: 200 krabů, 8 sloupců a žádné chybějící hodnoty
dim(x = data_krabi_zdroj)
sum(is.na(x = data_krabi_zdroj))

# Uložíme tabulku do složky data jako soubor CSV
dir.create(
  path = "data",
  showWarnings = FALSE
)
write.csv(
  x = data_krabi_zdroj,
  file = "data/krabi.csv",
  row.names = FALSE
)
```

Funkce `write.csv()` uloží textové sloupce v uvozovkách, takže se soubor může lišit vzhledem, ne hodnotami: po načtení pomocí `read.csv()` dostanete stejnou tabulku. Výuková kopie byla vytvořena z verze 7.3-66 balíčku MASS.
:::

::: {.box-data title="Původ a dosah dat" collapse="true"}
Campbell, N. A. a R. J. Mahon (1974): A multivariate study of variation in two species of rock crab of the genus *Leptograpsus*. *Australian Journal of Zoology* 22, 417–425. Tabulku šíří balíček MASS (`MASS::crabs`) pod licencí GPL-2 nebo GPL-3; soubor `krabi.csv` je její nezměněná kopie. Autoři vybrali do každé kombinace formy a pohlaví právě `r min(table(data_krabi_vse$sp, data_krabi_vse$sex))` krabů. Popis zdroje nedokládá náhodný `r render_glossary_term("vyber", display = "výběr")` ze všech krabů; naše závěry proto vztahujeme k těmto měřením, ne automaticky k celé populaci.
:::

## Které délky a šířky jsou v datech zastoupené?

Před modelem zkontrolujeme typy `r render_glossary_term("promenna", display = "proměnných")`, chybějící hodnoty a rozsahy. To nám také ukáže, při jakých délkách má smysl porovnávat předpovědi obou pohlaví.

```{r}
#| label: prohlednout-data
#| echo: true
# Ověříme typy sloupců a chybějící hodnoty
str(object = data_krabi)
colSums(x = is.na(data_krabi))

# Rozdělíme kraby podle pohlaví a prohlédneme každou skupinu zvlášť
data_samice <-
  data_krabi[data_krabi$pohlavi == "samice", ]
data_samci <-
  data_krabi[data_krabi$pohlavi == "samec", ]
summary(object = data_samice)
summary(object = data_samci)
```

V souboru je `r nrow(data_krabi)` krabů, z toho `r sum(data_krabi$pohlavi == "samice")` samic a `r sum(data_krabi$pohlavi == "samec")` samců. V těchto třech sloupcích chybí `r sum(is.na(data_krabi))` hodnot. Délka krunýře samic sahá od `r format_value_cz(min(data_samice$delka_krunyre_mm), digits = 1)` do `r format_value_cz(max(data_samice$delka_krunyre_mm), digits = 1)` mm, u samců od `r format_value_cz(min(data_samci$delka_krunyre_mm), digits = 1)` do `r format_value_cz(max(data_samci$delka_krunyre_mm), digits = 1)` mm. Obě skupiny tedy mají společný rozsah pozorovaných délek. To nám později umožní porovnat modelové předpovědi pro stejně dlouhou samici a samce bez extrapolace.

## Jak souvisí délka a zadní šířka krunýře?

Každý bod v následujícím **`r render_glossary_term("bodovy-graf", display = "bodovém grafu (anglicky scatter plot)")`** představuje jednoho kraba. Vodorovná souřadnice je délka krunýře, svislá zadní šířka. Barva a tvar bodu rozlišují pohlaví, takže graf zůstane čitelný i bez barev.

```{r}
#| label: bodovy-graf-krabu
#| echo: true
#| fig-alt: "Bodový graf zadní šířky krunýře proti délce krunýře u modrých krabů; body samic jsou zelené kruhy, body samců růžovofialové trojúhelníky a oba soubory bodů stoupají zleva doprava."
# Barvy a symboly přiřazujeme v celé lekci stejným skupinám
vec_barvy_pohlavi <-
  c(samice = "#009E73", samec = "#CC79A7")
vec_symboly_pohlavi <-
  c(samice = 16, samec = 17)

# Jeden bod představuje jednoho změřeného kraba
plot(
  x = data_krabi$delka_krunyre_mm,
  y = data_krabi$zadni_sirka_mm,
  col = vec_barvy_pohlavi[as.character(data_krabi$pohlavi)],
  pch = vec_symboly_pohlavi[as.character(data_krabi$pohlavi)],
  xlab = "Délka krunýře (mm)",
  ylab = "Zadní šířka krunýře (mm)",
  cex.axis = 1.2,
  cex.lab = 1.2
)
legend(
  x = "topleft",
  legend = names(vec_barvy_pohlavi),
  col = vec_barvy_pohlavi,
  pch = vec_symboly_pohlavi,
  bty = "n"
)
```

::: {.box-question title="Co vidíte v grafu?"}
Když se pohybujete zleva doprava, co se děje se zadní šířkou? Vypadají oba soubory bodů, jako by měly stejný `r render_glossary_term("sklon", display = "sklon")`?
:::

::: {.box-result collapse="true" title="Možná odpověď"}
Body naznačují rostoucí vztah u obou pohlaví. Samice mají při podobné délce obvykle větší zadní šířku, ale samotný pohled na body zatím nestačí k přesnému vyčíslení, zda se tento rozdíl s délkou mění.
:::

## Rovnoběžné přímky předpokládají společný sklon

Začneme známým aditivním modelem. Oba `r render_glossary_term("prediktor", display = "prediktory")` v něm mohou měnit předpovězenou šířku, ale **`r render_glossary_term("sklon", display = "sklon")` vůči délce je pro obě pohlaví společný**.

```{r}
#| label: aditivni-model
#| echo: true
# Aditivní model připouští posun mezi skupinami, ale stejný sklon
mod_rovnobezne <-
  lm(
    formula = zadni_sirka_mm ~ delka_krunyre_mm + pohlavi,
    data = data_krabi
  )
coef(object = mod_rovnobezne)
```

```{r}
#| label: obrazek-rovnobezne-primky
#| echo: false
#| fig-alt: "Body krabů obou pohlaví a dvě rovnoběžné přímky aditivního modelu; přímka samic leží nad přímkou samců a obě vedou jen přes délky pozorované v dané skupině."
# Each line spans only the lengths observed in its own group
data_primky <-
  data_krabi |>
  dplyr::group_by(pohlavi) |>
  dplyr::reframe(
    delka_krunyre_mm = seq(
      from = min(delka_krunyre_mm),
      to = max(delka_krunyre_mm),
      length.out = 100
    )
  )
data_primky$odhad_rovnobezne_mm <-
  predict(object = mod_rovnobezne, newdata = data_primky)

# Shared base plot for all figures of the first example
graf_krabi_zaklad <-
  ggplot2::ggplot(
    data = data_krabi,
    mapping = ggplot2::aes(
      x = delka_krunyre_mm,
      y = zadni_sirka_mm,
      colour = pohlavi,
      shape = pohlavi
    )
  ) +
  ggplot2::geom_point(alpha = 0.85, size = 2) +
  ggplot2::scale_colour_manual(values = vec_barvy_pohlavi) +
  ggplot2::scale_shape_manual(values = vec_symboly_pohlavi) +
  ggplot2::labs(
    x = "Délka krunýře (mm)",
    y = "Zadní šířka krunýře (mm)",
    colour = "Pohlaví",
    shape = "Pohlaví"
  ) +
  theme_biostat(base_size = 14)

graf_krabi_zaklad +
  ggplot2::geom_line(
    data = data_primky,
    mapping = ggplot2::aes(y = odhad_rovnobezne_mm),
    linewidth = 1.1
  )
```

::: {.box-question title="Co předpovídá aditivní model?"}
Co by tento model předpověděl o rozdílu v zadní šířce mezi stejně dlouhými kraby u kratších a delších délek z pozorovaného rozsahu?
:::

::: {.box-result collapse="true" title="Možná odpověď"}
Přímky zůstávají rovnoběžné, takže předpovězený rozdíl je v obou případech stejný. Je to vlastnost modelu, kterou jsme do něj vložili znakem `+`, ne výsledek samostatného testu této vlastnosti.
:::

## Mění se rozdíl mezi samicemi a samci s délkou krunýře?

V `r render_glossary_term("bodovy-graf", display = "bodovém grafu")` se rozdíl mezi skupinami při různých délkách nemusí jevit jako stálý. Na stejné osy proto můžeme orientačně vést přímku každou skupinou zvlášť. Srovnejte jejich směr s rovnoběžnými přímkami předchozího obrázku; body krabů jsou v obou grafech stejné.

```{r}
#| label: obrazek-navrh-ruznych-sklonu
#| echo: false
#| fig-alt: "Stejné body krabů s přímkou vedenou zvlášť každým pohlavím; přímka samic stoupá strměji než přímka samců a vzdálenost mezi nimi se s délkou zvětšuje."
graf_krabi_zaklad +
  ggplot2::geom_smooth(
    method = "lm",
    formula = y ~ x,
    se = FALSE,
    linewidth = 1.1
  )
```

::: {.box-question title="Předpovězte tvar"}
Pokud mají samice při určité společné délce větší předpovězenou zadní šířku a jejich přímka stoupá rychleji, bude rozdíl mezi přímkami u delších krabů menší, stejný, nebo větší?
:::

::: {.box-result collapse="true" title="Možná odpověď"}
Rozdíl by byl větší. Dvě různě strmé přímky totiž mění vzájemnou svislou vzdálenost s délkou.
:::

Když se svislá vzdálenost přímek mění s délkou, závisí vztah délky a šířky na pohlaví. Této závislosti říkáme `r render_glossary_term("interakce-predikatoru", display = "interakce prediktorů")`. Stejnou myšlenku můžeme vyjádřit jedním modelem, který oba `r render_glossary_term("sklon", display = "sklony")` odhadne společně.

## Model s interakcí dovoluje samicím a samcům různé sklony

Předchozí obrázek ukázal, proč zkusit model se dvěma různými `r render_glossary_term("sklon", display = "sklony")`. Délku krunýře v modelu měříme v milimetrech. Zápis `delka_krunyre_mm * pohlavi` dovolí, aby se sklon vůči délce mezi pohlavími lišil. Jednotlivé členy tohoto vzorce rozebereme po konkrétní předpovědi.

```{r}
#| label: model-s-interakci
#| echo: true
# Hvězdička připustí odlišné sklony pro samice a samce
mod_interakce <-
  lm(
    formula = zadni_sirka_mm ~ delka_krunyre_mm * pohlavi,
    data = data_krabi
  )
```

Předpovědní přímky kreslíme jen v rozsahu délek, které byly v dané skupině opravdu pozorovány. Za konci těchto úseků bychom už extrapolovali.

```{r}
#| label: obrazek-interakcni-primky
#| echo: false
#| fig-alt: "Body krabů a dvě nerovnoběžné přímky modelu s interakcí; přímka samic stoupá strměji než přímka samců."
data_primky$odhad_interakce_mm <-
  predict(object = mod_interakce, newdata = data_primky)

graf_interakce <-
  graf_krabi_zaklad +
  ggplot2::geom_line(
    data = data_primky,
    mapping = ggplot2::aes(y = odhad_interakce_mm),
    linewidth = 1.1
  )

graf_interakce
```

Přímka samic zde stoupá strměji než přímka samců. Čáry jsou `r render_glossary_term("odhadnuta-hodnota", display = "odhadnuté hodnoty")` modelu, zatímco body jsou skutečná měření. Svislou vzdálenost přímek vyčíslíme pro jednu společnou délku krunýře.

## Co model předpovídá pro stejně dlouhou samici a samce?

V grafu obou přímek můžeme porovnat jejich výšku nad stejným místem vodorovné osy. Zvolíme délku, která leží uvnitř pozorovaného rozsahu samic i samců a snadno ji najdeme na ose. Při srovnání pohlaví tak držíme délku krunýře stejnou.

```{r}
#| label: predpovedi-pro-stejnou-delku
#| echo: true
# Zvolíme společnou délku uvnitř pozorovaného rozsahu obou pohlaví
delka_porovnani_mm <- 30

# Připravíme stejnou délku pro samici i samce
data_porovnani <-
  data.frame(
    delka_krunyre_mm = rep(x = delka_porovnani_mm, times = 2),
    pohlavi = c("samice", "samec")
  )

# Model odhadne šířku pro obě kombinace
data_porovnani$odhad_sirky_mm <-
  predict(
    object = mod_interakce,
    newdata = data_porovnani
  )
data_porovnani
```

```{r}
#| label: vypocet-rozdilu-predpovedi
#| echo: false
odhad_samice_porovnani <-
  data_porovnani$odhad_sirky_mm[data_porovnani$pohlavi == "samice"]
odhad_samec_porovnani <-
  data_porovnani$odhad_sirky_mm[data_porovnani$pohlavi == "samec"]
rozdil_porovnani <- odhad_samec_porovnani - odhad_samice_porovnani
```

Zvolili jsme **`r delka_porovnani_mm` mm**. Která přímka je nad touto délkou výše? Oba odhady jsou v grafu označeny prázdnými kruhy. Svislá tečkovaná čára ukazuje místo, kde čteme jejich **svislou vzdálenost**; menší body zůstávají naměřenými kraby.

```{r}
#| label: obrazek-dve-predpovedi
#| echo: false
#| fig-alt: "Nerovnoběžné přímky modelu s interakcí a svislá tečkovaná čára při zvolené délce krunýře; dva prázdné kruhy označují odhad samice na horní přímce a odhad samce na dolní přímce."
graf_interakce +
  ggplot2::geom_vline(
    xintercept = delka_porovnani_mm,
    linetype = "dotted",
    colour = biostat_cols[["graphite"]]
  ) +
  ggplot2::geom_point(
    data = data_porovnani,
    mapping = ggplot2::aes(y = odhad_sirky_mm),
    shape = 21,
    fill = biostat_cols[["white"]],
    stroke = 1.7,
    size = 4,
    show.legend = FALSE
  )
```

::: {.box-result title="Rozdíl při zvolené délce"}
Při délce `r delka_porovnani_mm` mm model odhaduje zadní šířku samice `r format_value_cz(odhad_samice_porovnani, digits = 2)` mm a samce `r format_value_cz(odhad_samec_porovnani, digits = 2)` mm. **Rozdíl samec minus samice je `r format_value_cz(rozdil_porovnani, digits = 2, signed = TRUE)` mm.** Samčí předpověď tedy leží o `r format_value_cz(abs(rozdil_porovnani), digits = 2)` mm níže.
:::

Přímky nemají stejný `r render_glossary_term("sklon", display = "sklon")`, a proto se jejich svislá vzdálenost s délkou mění. Koeficienty modelu ukážou, o kolik se tento rozdíl změní při porovnání délek vzdálených o 1 mm.

## Co znamenají koeficienty modelu s interakcí?

Z minulé lekce známe aditivní model se společným `r render_glossary_term("sklon", display = "sklonem")`:

`zadni_sirka_mm ~ delka_krunyre_mm + pohlavi`

Aby se sklon mohl mezi pohlavími lišit, přidáme `r render_glossary_term("interakce-predikatoru", display = "interakční člen")` označený dvojtečkou:

`zadni_sirka_mm ~ delka_krunyre_mm + pohlavi + delka_krunyre_mm:pohlavi`

R nabízí pro celý pravý výraz kratší zápis s hvězdičkou:

`zadni_sirka_mm ~ delka_krunyre_mm * pohlavi`

Obecně tedy platí `x * y = x + y + x:y`. Výpis `summary()` ukáže čtyři koeficienty: `r render_glossary_term("intercept", display = "intercept")` a koeficient pro každý ze tří členů na pravé straně. Jde o tentýž model, který nakreslil obě přímky a vytvořil předpovědi při `r delka_porovnani_mm` mm:

```{r}
#| label: souhrn-interakce
#| echo: true
summary(object = mod_interakce)
```

V tabulce **Coefficients** jsou čtyři řádky. Sloupec **Estimate** obsahuje `r render_glossary_term("odhad", display = "odhady")` v milimetrech nebo v milimetrech šířky na milimetr délky. Referenční skupinou jsou samice, protože jsme je při tvorbě faktoru `pohlavi` uvedli jako první.

```{r}
#| label: odhady-koeficientu
#| echo: false
vec_odhad_koeficientu <- coef(mod_interakce)
sklon_samice <- unname(vec_odhad_koeficientu["delka_krunyre_mm"])
rozdil_sklonu_samec_minus_samice <-
  unname(vec_odhad_koeficientu["delka_krunyre_mm:pohlavisamec"])
sklon_samec <- sklon_samice + rozdil_sklonu_samec_minus_samice
interval_rozdilu_sklonu <-
  confint(mod_interakce)["delka_krunyre_mm:pohlavisamec", ]
prispevek_delky_samice <- delka_porovnani_mm * sklon_samice
prispevek_delky_rozdilu <-
  delka_porovnani_mm * rozdil_sklonu_samec_minus_samice
```

- **`(Intercept)`**, tedy intercept, je odhad zadní šířky **samice s délkou krunýře 0 mm**: `r format_value_cz(vec_odhad_koeficientu[["(Intercept)"]], digits = 3)` mm. Takto malého kraba jsme neměřili; číslo samo o sobě biologicky nevykládáme.
- **`delka_krunyre_mm`** je sklon **samičí** přímky: při délce o 1 mm větší model odhaduje šířku větší o `r format_value_cz(sklon_samice, digits = 3)` mm. Nejde o růst jednoho kraba v čase, ale o porovnání různě dlouhých jedinců.
- **`pohlavisamec`** je rozdíl **samec minus samice při délce 0 mm**: `r format_value_cz(vec_odhad_koeficientu[["pohlavisamec"]], digits = 3, signed = TRUE)` mm. Není to rozdíl mezi pohlavími při každé délce ani odhad pro některého ze změřených krabů.
- **`delka_krunyre_mm:pohlavisamec`** je interakční člen: rozdíl **sklon samce minus sklon samice**, `r format_value_cz(rozdil_sklonu_samec_minus_samice, digits = 3, signed = TRUE)` mm/mm. Záporná hodnota říká, že samčí přímka stoupá pomaleji.

::: {.box-data title="Doplňující: proč s interakcí ponecháváme i oba hlavní členy" collapse="true"}

Hvězdička rozvine model na `delka_krunyre_mm + pohlavi + delka_krunyre_mm:pohlavi`. První dva členy určují referenční přímku a rozdíl mezi skupinami v referenční hodnotě délky; poslední člen popisuje rozdíl sklonů. Tomuto uspořádání se říká hierarchický princip modelu. <!-- TODO: add slug hierarchicky-princip-modelu to pojmy.yaml -->

Kdybychom hlavní člen délky nebo pohlaví bez věcného důvodu odstranili, interakční koeficient by zůstal bez přirozeného základu, vůči kterému se interpretuje. Výjimky existují pro úzce formulované modely, ale v úvodní analýze ponecháváme oba hlavní členy a biologickou otázku čteme z jejich společných předpovědí.
:::

Hodnota 0 mm určuje pouze to, co v tomto zápisu znamenají dva koeficienty. Pro srovnání skutečně měřených délek použijeme oba koeficienty sklonu a dosadíme délku z pozorovaného rozsahu. **Rozdíl `r format_value_cz(rozdil_porovnani, digits = 2, signed = TRUE)` mm jsme už odečetli z předpovědí při `r delka_porovnani_mm` mm.** Následující výpočet nepřináší nový výsledek; ukazuje, jak se stejná hodnota sestaví ze čtyř řádků tabulky **Coefficients**. Odhad šířky samice a rozdíl **samec minus samice** mají tuto stavbu:

$$
\begin{aligned}
\text{odhad šířky samice při zvolené délce}
&= \text{odhad šířky samice při 0 mm} \\
&\quad + \text{zvolená délka}\times\text{samičí sklon} \\
\text{rozdíl samec minus samice při zvolené délce}
&= \text{rozdíl při 0 mm} \\
&\quad + \text{zvolená délka}\times\text{rozdíl sklonů}
\end{aligned}
$$

Při `r delka_porovnani_mm` mm dosadíme odhady z tabulky **Coefficients**:

$$
\begin{aligned}
\text{odhad šířky samice}
&= `r format_value_cz(vec_odhad_koeficientu[["(Intercept)"]], digits = 4, math = TRUE)`
 + `r delka_porovnani_mm`\cdot `r format_value_cz(sklon_samice, digits = 4, math = TRUE)`
 \approx `r format_value_cz(odhad_samice_porovnani, digits = 2, math = TRUE)`\,\text{mm} \\
\text{rozdíl samec minus samice}
&= `r format_value_cz(vec_odhad_koeficientu[["pohlavisamec"]], digits = 4, math = TRUE)`
 + `r delka_porovnani_mm`\cdot (`r format_value_cz(rozdil_sklonu_samec_minus_samice, digits = 4, signed = TRUE, math = TRUE)`)
 \approx `r format_value_cz(rozdil_porovnani, digits = 2, signed = TRUE, math = TRUE)`\,\text{mm}
\end{aligned}
$$

Ve výpočtu používáme nezaokrouhlené koeficienty modelu; zobrazená čtyři desetinná místa jsou jen pro čitelnost. Po vynásobení se výpočet zkrátí na dva součty:

$$
\begin{aligned}
\text{odhad šířky samice}
&\approx `r format_value_cz(vec_odhad_koeficientu[["(Intercept)"]], digits = 4, math = TRUE)`
 + `r format_value_cz(prispevek_delky_samice, digits = 4, math = TRUE)`
 \approx `r format_value_cz(odhad_samice_porovnani, digits = 2, math = TRUE)`\,\text{mm} \\
\text{rozdíl samec minus samice}
&\approx `r format_value_cz(vec_odhad_koeficientu[["pohlavisamec"]], digits = 4, math = TRUE)`
 + (`r format_value_cz(prispevek_delky_rozdilu, digits = 4, signed = TRUE, math = TRUE)`)
 \approx `r format_value_cz(rozdil_porovnani, digits = 2, signed = TRUE, math = TRUE)`\,\text{mm}
\end{aligned}
$$

Stejný postup lze zapsat stručně pro libovolnou délku:

$$
\begin{aligned}
\hat{y}_{\mathrm{samice}}(x) &= \beta_0 + \beta_1 x \\
\hat{y}_{\mathrm{samec}}(x)-\hat{y}_{\mathrm{samice}}(x) &= \beta_2 + \beta_3 x
\end{aligned}
$$

Zde $x$ znamená zvolenou délku krunýře v mm a $\hat{y}$ modelový odhad zadní šířky v mm. $\beta_0$ je řádek `(Intercept)` pro samici při 0 mm, $\beta_1$ její sklon, $\beta_2$ řádek `pohlavisamec` pro rozdíl při 0 mm a $\beta_3$ interakční řádek pro rozdíl sklonů.

Odhad šířky samce získáme přičtením druhého výsledku k prvnímu: `r format_value_cz(odhad_samice_porovnani, digits = 2)` + (`r format_value_cz(rozdil_porovnani, digits = 2, signed = TRUE)`) = `r format_value_cz(odhad_samec_porovnani, digits = 2)` mm. To souhlasí s řádkem samce v tabulce předpovědí. Pro biologické porovnání čteme **předpovědi při konkrétní délce** a jejich rozdíl.

Samčí sklon získáme sečtením samičího sklonu a interakčního koeficientu. Vyjde `r format_value_cz(sklon_samec, digits = 3)` mm/mm oproti `r format_value_cz(sklon_samice, digits = 3)` mm/mm u samic. Rozdíl **samec minus samice** se při porovnání o 1 mm delších krabů změní o `r format_value_cz(rozdil_sklonu_samec_minus_samice, digits = 3, signed = TRUE)` mm. Pro tento rozdíl sklonů získáme z téhož modelu 95% `r render_glossary_term("interval-spolehlivosti", display = "interval spolehlivosti")`:

```{r}
#| label: interval-rozdilu-sklonu
#| echo: true
# Interval patří rozdílu sklonů samce a samice
confint(
  object = mod_interakce,
  parm = "delka_krunyre_mm:pohlavisamec"
)
```

Odhadovaný rozdíl sklonů je `r format_value_cz(rozdil_sklonu_samec_minus_samice, digits = 3)` mm/mm; interval sahá od `r format_value_cz(interval_rozdilu_sklonu[[1]], digits = 3)` do `r format_value_cz(interval_rozdilu_sklonu[[2]], digits = 3)` mm/mm. Popisuje `r render_glossary_term("nejistota", display = "nejistotu")` odhadu při předpokladech modelu, nikoli rozsah šířek jednotlivých krabů.

::: {.box-warning title="Koeficient pohlaví neplatí při každé délce"}
Samotný koeficient `pohlavisamec` popisuje rozdíl při **0 mm**, tedy mimo pozorované délky. Při biologickém porovnání vždy uveďte **konkrétní délku** a porovnejte dvě modelové předpovědi při této délce.
:::

::: {.box-data title="Doplňující: co dalšího ukazuje summary()?" collapse="true"}
Sloupec **Std. Error** popisuje nejistotu jednotlivých koeficientů. **t value** a **Pr(>|t|)** slouží k testování, zda by příslušný *podmíněný* koeficient mohl být nulový při předpokladech modelu. Hvězdičky u `r render_glossary_term("p-hodnota", display = "p-hodnot")` nejsou velikostí biologického rozdílu: například test u `(Intercept)` se ptá na nulovou šířku samice při délce 0 mm, což zde není zajímavá biologická otázka.

Část **Residuals** stručně shrnuje odchylky jednotlivých krabů od přímek. **Residual standard error** je `r render_glossary_term("rezidualni-smerodatna-odchylka", display = "reziduální směrodatná odchylka")`; popisuje typickou velikost těchto odchylek v milimetrech. Údaje **Multiple R-squared**, **Adjusted R-squared** a **F-statistic** zatím nepotřebujeme k výkladu interakce; k `r render_glossary_term("porovnani-modelu", display = "porovnávání")` biologicky zdůvodněných modelů se vrátíme v následující lekci.
:::

::: {.box-data title="Doplňující: modelové průměry při zvolené délce pomocí emmeans" collapse="true"}
Funkci `emmeans()` znáte z předchozích lekcí. Pracuje se stejným odhadnutým modelem jako naše tabulka předpovědí; nastavíme v ní délku krunýře na `r delka_porovnani_mm` mm.

```{r}
#| label: extra-emmeans-sirka
#| echo: true
# Odhadneme průměrnou zadní šířku obou pohlaví při zvolené délce
prumery_sirka <-
  emmeans::emmeans(
    object = mod_interakce,
    specs = ~ pohlavi | delka_krunyre_mm,
    at = list(delka_krunyre_mm = delka_porovnani_mm)
  )
prumery_sirka
```

Sloupec `emmean` obsahuje stejné dvě bodové předpovědi jako naše tabulka z `predict()`. Sloupce `lower.CL` a `upper.CL` vymezují modelový interval pro **`r render_glossary_term("prumer", display = "průměrnou")` šířku skupiny**, ne rozsah šířek jednotlivých krabů. Chceme-li interval přímo pro rozdíl **samec minus samice**, vypočteme jej zvlášť:

```{r}
#| label: extra-emmeans-rozdily-sirky
#| echo: true
# Odečteme samičí odhad od samčího odhadu
rozdily_sirka <-
  emmeans::contrast(
    object = prumery_sirka,
    method = "revpairwise"
  )
confint(object = rozdily_sirka)
```

Odhad rozdílu odpovídá dříve vypočteným `r format_value_cz(rozdil_porovnani, digits = 2, signed = TRUE)` mm. **Interval pro rozdíl vzniká z obou předpovědí a jejich společné nejistoty**; nelze jej spolehlivě vyčíst jen z překryvu dvou intervalů pro jednotlivé průměry. Způsob výběru krabů stále omezuje, na koho lze výsledek přenést.
:::

## Jak dobře přímky popisují změřené kraby?

Přímky procházejí oběma mračny bodů, ale jednotlivé šířky se od nich odchylují. Pro každého kraba je `r render_glossary_term("residuum", display = "residuum")` rozdíl mezi naměřenou a `r render_glossary_term("odhadnuta-hodnota", display = "odhadnutou šířkou")`. Graf residuí proti odhadům pomůže odhalit zakřivení nebo změnu `r render_glossary_term("rozptyl", display = "rozptylu")`, které by dvě přímky nezachytily.

```{r}
#| label: diagnostika-interakcniho-modelu
#| echo: true
#| fig-alt: "Residua modelu s interakcí proti odhadnuté zadní šířce; body samic a samců leží nad i pod vodorovnou přerušovanou nulou bez oblouku a s podobným rozptylem v celém rozsahu odhadů."
# Spočteme rozdíl mezi měřením a předpovědí pro každého kraba
data_krabi$odhad_sirky_mm <-
  fitted(object = mod_interakce)
data_krabi$residuum_mm <-
  data_krabi$zadni_sirka_mm - data_krabi$odhad_sirky_mm

# Nula znamená dokonalou shodu měření a předpovědi
plot(
  x = data_krabi$odhad_sirky_mm,
  y = data_krabi$residuum_mm,
  col = vec_barvy_pohlavi[as.character(data_krabi$pohlavi)],
  pch = vec_symboly_pohlavi[as.character(data_krabi$pohlavi)],
  xlab = "Odhadnutá zadní šířka (mm)",
  ylab = "Residuum (mm)",
  cex.axis = 1.2,
  cex.lab = 1.2
)
abline(h = 0, lty = 2)
```

Hledáme soustavné zakřivení nebo nápadně větší rozptyl na jedné straně: tím bychom poznali, že dvě přímky některou část vztahu nevystihují. V tomto grafu nic takového nevidíme. Residua leží nad nulou i pod ní v celém rozsahu odhadů, netvoří oblouk a jejich rozptyl je u malých i velkých odhadů podobný. Vraťte se také k `r render_glossary_term("bodovy-graf", display = "bodovému grafu")` a prohlédněte kraby na koncích rozsahu délek. V těchto místech je méně sousedních měření, a jednotlivý krajní bod může `r render_glossary_term("sklon", display = "sklon")` přímky více ovlivnit. Takový bod je důvod k ověření měření a citlivosti závěru, ne k automatickému odstranění.

`r render_glossary_term("diagnostika-modelu", display = "Diagnostika modelu")` sama nezaručí, že je vztah v populaci přesně lineární, že výběr krabů byl reprezentativní, ani že pohlaví něco způsobilo. Předpovědi mimo pozorované délky jsme nekontrolovali.

::: {.box-data title="Doplňující: interakci můžeme číst jen tam, kde máme data" collapse="true"}

Rozdíl sklonů je odhadnut z oblastí, kde máme měření obou pohlaví. Kdyby samice byly pozorovány jen mezi kratšími délkami a samci jen mezi delšími délkami bez společného rozsahu, model by sice dokázal nakreslit dvě přímky, ale jejich srovnání při stejné délce by bylo založeno hlavně na extrapolaci.

Před interpretací interakce proto kontrolujeme nejen počet bodů, ale také **společný rozsah hodnot** a zastoupení všech důležitých kombinací. Řídká nebo prázdná část grafu znamená, že modelová předpověď je méně podepřená skutečnými srovnáními. Více parametrů nemůže vytvořit informaci v kombinaci, kterou jsme nepozorovali.
:::

::: {.box-result title="Co ukázal první příklad"}
U těchto změřených krabů modré formy roste odhadnutá zadní šířka s délkou u obou pohlaví. Sklon je u samic přibližně `r format_value_cz(sklon_samice, digits = 3)` mm šířky na 1 mm délky a u samců `r format_value_cz(sklon_samec, digits = 3)` mm/mm. Při `r delka_porovnani_mm` mm je modelový rozdíl **samec minus samice** `r format_value_cz(rozdil_porovnani, digits = 2, signed = TRUE)` mm. Samčí sklon je o `r format_value_cz(abs(rozdil_sklonu_samec_minus_samice), digits = 3)` mm/mm menší, takže u krabů s větší délkou model předpovídá ještě větší náskok samic. To odpovídá biologickému očekávání z úvodu, že samice s rostoucí velikostí mají relativně širší zadní část krunýře. Jde však o souvislost v těchto měřeních, nikoli o důkaz, že ji způsobuje pohlaví nebo nošení snůšky.
:::

## Liší se rozdíl v délce krunýře mezi pohlavími u obou barevných forem?

Vztah zadní šířky k délce jsme sledovali pouze u modré formy. Soubor obsahuje také oranžovou barevnou formu. U mnoha živočichů se samice a samci liší velikostí a tento rozdíl nemusí být stejný u příbuzných druhů nebo forem. Můžeme se proto zeptat na jiný rozměr týchž krabů: **je rozdíl v délce krunýře mezi samci a samicemi stejný u obou forem?** `r render_glossary_term("odezva", display = "Odezvou")` je tentokrát délka krunýře v milimetrech. Barevná forma i pohlaví jsou `r render_glossary_term("kategoricka-promenna", display = "kategoriální proměnné")` neboli skupinové znaky; pro tento příklad nepotřebujeme číselný `r render_glossary_term("prediktor", display = "prediktor")`.

::: {.box-warning title="Co tato data o velikosti krabů neřeknou"}
Autoři vybrali do každé ze čtyř skupin právě `r min(table(data_krabi_vse$sp, data_krabi_vse$sex))` krabů. `r render_glossary_term("prumer", display = "Průměrné")` délky proto závisejí také na tom, které kraby sesbírali a změřili, nejen na tom, jak velcí krabi dané formy a pohlaví v přírodě obvykle jsou. Příklad nám ukáže, jak číst model se dvěma skupinovými znaky; tvrzení o velikosti krabů v přírodě by vyžadovalo náhodný `r render_glossary_term("vyber", display = "výběr")`.
:::

```{r}
#| label: pripravit-obe-formy
#| echo: true
# Ze stejné tabulky vybereme novou odezvu i oba skupinové znaky
data_formy <-
  data_krabi_vse[, c("sp", "sex", "CL")]
names(data_formy) <-
  c("forma", "pohlavi", "delka_krunyre_mm")

# Obě skupinové proměnné převedeme na faktory s referenční úrovní
data_formy$forma <-
  factor(
    x = data_formy$forma,
    levels = c("B", "O"),
    labels = c("modrá", "oranžová")
  )
data_formy$pohlavi <-
  factor(
    x = data_formy$pohlavi,
    levels = c("F", "M"),
    labels = c("samice", "samec")
  )

# Počty ve čtyřech skupinách a chybějící délky
table(data_formy$forma, data_formy$pohlavi)
colSums(x = is.na(data_formy))
summary(object = data_formy$delka_krunyre_mm)
```

V souboru je `r nrow(data_formy)` krabů. V každé ze čtyř kombinací je `r min(table(data_formy$forma, data_formy$pohlavi))` měření; v délce krunýře chybí `r sum(is.na(data_formy$delka_krunyre_mm))` hodnot. Pozorované délky v celém souboru sahají od `r format_value_cz(min(data_formy$delka_krunyre_mm), digits = 1)` do `r format_value_cz(max(data_formy$delka_krunyre_mm), digits = 1)` mm. Jeden řádek stále představuje jednoho kraba. Názvy forem označují dvě skupiny v datovém souboru, nikoli zásah, který bychom krabům přidělili.

## Jaké délky krunýře vidíme ve čtyřech skupinách?

Podívejte se na jednotlivé naměřené délky; každý bod představuje jednoho kraba. Body samic a samců jsou uvnitř každé barevné formy mírně posunuté do stran, aby se méně překrývaly; barvy a tvary bodů pro pohlaví zůstávají stejné jako v předchozích grafech. Modrý a oranžový kroužek u popisku připomíná barevnou formu.

```{r}
#| label: zobrazit-ctyri-skupiny
#| echo: false
#| fig-alt: "Délky krunýře jednotlivých krabů ve čtyřech skupinách: samice a samci modré formy vlevo, samice a samci oranžové formy vpravo; ve všech skupinách se délky výrazně překrývají."
# Position on the x axis keeps the form and separates the sexes within it
posun_pohlavi <- c(samice = -0.13, samec = 0.13)

# Blue and orange circles name the two colour forms; the sexes keep their
#   own colours, so the two groupings cannot be confused
vec_popisky_forem <-
  c("modrá" = "🔵 modrá forma", "oranžová" = "🟠 oranžová forma")
data_formy$x_skupina <-
  as.numeric(data_formy$forma) +
  posun_pohlavi[as.character(data_formy$pohlavi)]

graf_ctyri_skupiny <-
  ggplot2::ggplot(
    data = data_formy,
    mapping = ggplot2::aes(
      x = x_skupina,
      y = delka_krunyre_mm,
      colour = pohlavi,
      shape = pohlavi
    )
  ) +
  ggplot2::geom_point(
    position = ggplot2::position_jitter(
      width = 0.07,
      height = 0,
      seed = 900723
    ),
    alpha = 0.8,
    size = 1.9
  ) +
  ggplot2::scale_x_continuous(
    breaks = seq_along(levels(data_formy$forma)),
    labels = vec_popisky_forem[levels(data_formy$forma)],
    limits = c(0.65, 2.35)
  ) +
  ggplot2::scale_colour_manual(values = vec_barvy_pohlavi) +
  ggplot2::scale_shape_manual(values = vec_symboly_pohlavi) +
  ggplot2::labs(
    x = NULL,
    y = "Délka krunýře (mm)",
    colour = "Pohlaví",
    shape = "Pohlaví"
  ) +
  theme_biostat(base_size = 14)

graf_ctyri_skupiny
```

```{r}
#| label: prumery-ctyr-skupin
#| echo: false
data_prumery_skupin <-
  data_formy |>
  dplyr::group_by(forma, pohlavi) |>
  dplyr::summarise(
    prumer_mm = mean(delka_krunyre_mm),
    .groups = "drop"
  )
data_prumery_skupin$x_skupina <-
  as.numeric(data_prumery_skupin$forma) +
  posun_pohlavi[as.character(data_prumery_skupin$pohlavi)]

mat_prumery_formy <-
  xtabs(
    formula = prumer_mm ~ forma + pohlavi,
    data = data_prumery_skupin
  )

knitr::kable(
  x = data.frame(
    forma = rownames(mat_prumery_formy),
    samice = round(unname(mat_prumery_formy[, "samice"]), digits = 2),
    samec = round(unname(mat_prumery_formy[, "samec"]), digits = 2)
  ),
  col.names = c("Barevná forma", "Samice: průměr (mm)", "Samec: průměr (mm)"),
  caption = "Pozorované průměrné délky krunýře ve čtyřech skupinách.",
  row.names = FALSE
)

rozdil_pohlavi_modra <-
  unname(mat_prumery_formy["modrá", "samec"] - mat_prumery_formy["modrá", "samice"])
rozdil_pohlavi_oranzova <-
  unname(mat_prumery_formy["oranžová", "samec"] - mat_prumery_formy["oranžová", "samice"])
rozdil_rozdilu_formy <- rozdil_pohlavi_oranzova - rozdil_pohlavi_modra
```

::: {.box-question title="Porovnejte obě formy"}
Ve které formě mají samci větší `r render_glossary_term("prumer", display = "průměrnou")` délku než samice? Platí to i ve druhé formě?
:::

::: {.box-result collapse="true" title="Možná odpověď"}
U modré formy činí rozdíl **samec minus samice** `r format_value_cz(rozdil_pohlavi_modra, digits = 2, signed = TRUE)` mm. U oranžové formy je to `r format_value_cz(rozdil_pohlavi_oranzova, digits = 2, signed = TRUE)` mm, takže rozdíl má opačný směr.
:::

Průměry popisují skupiny, zatímco jednotlivé body v obrázku připomínají velkou proměnlivost délek mezi kraby uvnitř každé skupiny.

## Jak model zachytí změnu rozdílu mezi pohlavími podle formy?

Vrátíme se ke **stejnému obrázku se všemi jednotlivými kraby**. Body ani osy neměníme. Přidáme pouze velké kosočtverce pro čtyři `r render_glossary_term("prumer", display = "průměry")` a v každé barevné formě spojíme průměr samic s průměrem samců. Sklon spojnice zde nevyjadřuje spojitou předpověď; pouze zvýrazňuje velikost a směr rozdílu mezi dvěma kategoriemi.

```{r}
#| label: obrazek-rozdilu-mezi-formami
#| echo: false
#| fig-alt: "Stejné body čtyř skupin doplněné o čtyři velké kosočtverce s průměry; u modré formy spojnice stoupá od samic k samcům, u oranžové formy mírně klesá."
graf_ctyri_skupiny +
  ggplot2::geom_line(
    data = data_prumery_skupin,
    mapping = ggplot2::aes(
      x = x_skupina,
      y = prumer_mm,
      group = forma
    ),
    inherit.aes = FALSE,
    linewidth = 1,
    colour = biostat_cols[["graphite"]]
  ) +
  ggplot2::geom_point(
    data = data_prumery_skupin,
    mapping = ggplot2::aes(
      x = x_skupina,
      y = prumer_mm,
      fill = pohlavi
    ),
    inherit.aes = FALSE,
    shape = 23,
    size = 4.5,
    colour = biostat_cols[["white"]],
    show.legend = FALSE
  ) +
  ggplot2::scale_fill_manual(values = vec_barvy_pohlavi)
```

Na modré formě leží průměr samců nad průměrem samic; na oranžové pod ním. Zajímá nás, **o kolik se rozdíl samec minus samice změní mezi formami**. Tento rozdíl dvou rozdílů je druhým typem `r render_glossary_term("interakce-predikatoru", display = "interakce")`:

$$
\begin{aligned}
\text{změna rozdílu mezi pohlavími}
&= \text{průměr délky oranžových samců} \\
&\quad - \text{průměr délky oranžových samic} \\
&\quad - \text{průměr délky modrých samců} \\
&\quad + \text{průměr délky modrých samic}
\end{aligned}
$$

Vezmeme čtyři průměry z tabulky a odečteme vždy dvě hodnoty ze stejného řádku:

$$
\begin{aligned}
\text{změna rozdílu mezi pohlavími}
&= (`r format_value_cz(mat_prumery_formy["oranžová", "samec"], digits = 3, math = TRUE)` - `r format_value_cz(mat_prumery_formy["oranžová", "samice"], digits = 3, math = TRUE)`) - (`r format_value_cz(mat_prumery_formy["modrá", "samec"], digits = 3, math = TRUE)` - `r format_value_cz(mat_prumery_formy["modrá", "samice"], digits = 3, math = TRUE)`) \\
&= `r format_value_cz(rozdil_pohlavi_oranzova, digits = 3, math = TRUE)` - `r format_value_cz(rozdil_pohlavi_modra, digits = 3, math = TRUE)` \\
&= `r format_value_cz(rozdil_rozdilu_formy, digits = 3, math = TRUE)`\,\text{mm}
\end{aligned}
$$

Stručně lze stejné čtyři průměry označit $\bar y_{O,M}$, $\bar y_{O,F}$, $\bar y_{B,M}$ a $\bar y_{B,F}$; $O$ a $B$ označují oranžovou a modrou formu, $M$ a $F$ samce a samici a pruh nad $y$ průměr délky. Pak platí $\Delta = (\bar y_{O,M} - \bar y_{O,F}) - (\bar y_{B,M} - \bar y_{B,F})$, kde $\Delta$ označuje změnu rozdílu mezi pohlavími při přechodu od modré k oranžové formě.

Jediným modelem můžeme tyto čtyři skupinové průměry odhadnout společně. Hvězdička zahrne obě skupiny i jejich interakci:

```{r}
#| label: model-dvou-skupinovych-znaku
#| echo: true
# Referenční kombinací je modrá forma a samice
mod_formy <-
  lm(
    formula = delka_krunyre_mm ~ forma * pohlavi,
    data = data_formy
  )
summary(object = mod_formy)
```

V tabulce **Coefficients** se opět soustřeďte na sloupec **Estimate**. Řádky nejsou čtyři nezávislé průměry: jeden je referenční průměr a zbylé tři popisují, jak se od něj dostaneme k ostatním skupinám.

## Co znamenají čtyři koeficienty modelu se dvěma skupinovými znaky?

Protože referenčními kategoriemi jsou **modrá forma** a **samice**, každý koeficient má přesně určené srovnání. `r render_glossary_term("odhad", display = "Odhady")` vytáhneme ze stejného modelu, který vytvořil výpis `summary()`.

```{r}
#| label: koeficienty-dvou-skupinovych-znaku
#| echo: false
vec_koeficienty_formy <- coef(mod_formy)
vec_interval_formy <- confint(mod_formy)["formaoranžová:pohlavisamec", ]
```

- **`(Intercept)`** je odhad `r render_glossary_term("prumer", display = "průměrné")` délky **modré samice**: `r format_value_cz(vec_koeficienty_formy[["(Intercept)"]], digits = 3)` mm. Najdete jej v tabulce průměrů v řádku modré formy a sloupci samic.
- **`formaoranžová`** je rozdíl **oranžová minus modrá mezi samicemi**: `r format_value_cz(vec_koeficienty_formy[["formaoranžová"]], digits = 3, signed = TRUE)` mm. Po přičtení k `r render_glossary_term("intercept", display = "interceptu")` dostaneme průměr oranžových samic; tento koeficient nepopisuje rozdíl forem u obou pohlaví stejně.
- **`pohlavisamec`** je rozdíl **samec minus samice v modré formě**: `r format_value_cz(vec_koeficienty_formy[["pohlavisamec"]], digits = 3, signed = TRUE)` mm. Souhlasí se svislou vzdáleností obou průměrů u modré formy.
- **`formaoranžová:pohlavisamec`** je změna rozdílu samec minus samice při přechodu **od modré k oranžové formě**: `r format_value_cz(vec_koeficienty_formy[["formaoranžová:pohlavisamec"]], digits = 3, signed = TRUE)` mm. Jde přesně o výše vypočtený rozdíl dvou rozdílů, nikoli o další samostatnou skupinu.

```{r}
#| label: nejistota-rozdilu-rozdilu
#| echo: true
# Interval patří k rozdílu dvou rozdílů, ne k jedné skupině
confint(
  object = mod_formy,
  parm = "formaoranžová:pohlavisamec"
)
```

Modelový 95% `r render_glossary_term("interval-spolehlivosti", display = "interval spolehlivosti")` pro změnu rozdílu sahá od `r format_value_cz(vec_interval_formy[[1]], digits = 2)` do `r format_value_cz(vec_interval_formy[[2]], digits = 2)` mm. `r if (vec_interval_formy[[2]] < 0) "Celý leží pod nulou: při předpokladech modelu data podporují menší rozdíl samec minus samice u oranžové než u modré formy." else if (vec_interval_formy[[1]] > 0) "Celý leží nad nulou: při předpokladech modelu data podporují větší rozdíl samec minus samice u oranžové než u modré formy." else "Zahrnuje také nulovou změnu, takže směr rozdílu mezi formami zůstává podle tohoto modelu nejistý."` Protože autoři vybrali do každé skupiny pevný počet krabů, nelze tento interval automaticky zobecnit na všechny kraby.

Jak model sestaví předpověď pro oranžového samce? K odhadu modré samice přidá rozdíl forem mezi samicemi, rozdíl pohlaví u modré formy a změnu tohoto rozdílu u oranžové formy:

$$
\begin{aligned}
\text{odhad délky oranžového samce}
&= \text{odhad délky modré samice} \\
&\quad + \text{rozdíl forem mezi samicemi} \\
&\quad + \text{rozdíl pohlaví u modré formy} \\
&\quad + \text{změna rozdílu pohlaví u oranžové formy}
\end{aligned}
$$

Dosazené hodnoty jsou právě čtyři řádky výpisu modelu:

$$
\begin{aligned}
\text{odhad délky oranžového samce}
&= `r format_value_cz(vec_koeficienty_formy[["(Intercept)"]], digits = 3, math = TRUE)` + `r format_value_cz(vec_koeficienty_formy[["formaoranžová"]], digits = 3, math = TRUE)` + `r format_value_cz(vec_koeficienty_formy[["pohlavisamec"]], digits = 3, math = TRUE)` + (`r format_value_cz(vec_koeficienty_formy[["formaoranžová:pohlavisamec"]], digits = 3, math = TRUE)`) \\
&= `r format_value_cz(sum(vec_koeficienty_formy[c("(Intercept)", "formaoranžová", "pohlavisamec")]), digits = 3, math = TRUE)` + (`r format_value_cz(vec_koeficienty_formy[["formaoranžová:pohlavisamec"]], digits = 3, math = TRUE)`) \\
&= `r format_value_cz(sum(vec_koeficienty_formy), digits = 3, math = TRUE)`\,\text{mm}
\end{aligned}
$$

V obecném zkráceném zápisu je to $\widehat y_{O,M}=\beta_0+\beta_O+\beta_M+\beta_{OM}$. Symbol $\widehat y_{O,M}$ značí modelový odhad délky oranžového samce, $\beta_0$ referenční průměr modré samice, $\beta_O$ rozdíl forem u samic, $\beta_M$ rozdíl pohlaví u modré formy a $\beta_{OM}$ změnu tohoto rozdílu u oranžové formy.

```{r}
#| label: predpovedi-ctyr-skupin
#| echo: true
# Vytvoříme všechny čtyři kombinace a pro každou spočteme odhad
data_kombinace <-
  expand.grid(
    forma = levels(data_formy$forma),
    pohlavi = levels(data_formy$pohlavi)
  )
data_kombinace$delka_odhad_mm <-
  predict(
    object = mod_formy,
    newdata = data_kombinace
  )
data_kombinace
```

```{r}
#| label: tabulka-predpovedi-ctyr-skupin
#| echo: false
data_predikce_formy <-
  data.frame(
    forma = levels(data_formy$forma),
    samice = data_kombinace$delka_odhad_mm[data_kombinace$pohlavi == "samice"],
    samec = data_kombinace$delka_odhad_mm[data_kombinace$pohlavi == "samec"]
  )

knitr::kable(
  x = data.frame(
    forma = data_predikce_formy$forma,
    samice = round(data_predikce_formy$samice, digits = 2),
    samec = round(data_predikce_formy$samec, digits = 2)
  ),
  col.names = c("Barevná forma", "Samice: odhad (mm)", "Samec: odhad (mm)"),
  caption = "Čtyři odhady z jediného modelu; srovnání pohlaví probíhá uvnitř řádku formy.",
  row.names = FALSE
)

rozdil_oranzova_predikce <-
  data_predikce_formy$samec[data_predikce_formy$forma == "oranžová"] -
  data_predikce_formy$samice[data_predikce_formy$forma == "oranžová"]
```

V řádku oranžové formy je rozdíl samec minus samice `r format_value_cz(rozdil_oranzova_predikce, digits = 2, signed = TRUE)` mm. Získáme jej také sečtením koeficientů `pohlavisamec` a `formaoranžová:pohlavisamec`: `r format_value_cz(vec_koeficienty_formy[["pohlavisamec"]], digits = 3)` + (`r format_value_cz(vec_koeficienty_formy[["formaoranžová:pohlavisamec"]], digits = 3)`) = `r format_value_cz(sum(vec_koeficienty_formy[c("pohlavisamec", "formaoranžová:pohlavisamec")]), digits = 3)` mm. **Samotný koeficient pohlaví tedy platí pro modrou formu; u oranžové se k němu přidá interakce.**

::: {.box-data title="Doplňující: čtyři skupinové předpovědi pomocí emmeans" collapse="true"}
Stejně jako v předchozím příkladu můžeme z již odhadnutého modelu získat předpovědi pomocí `emmeans()`. Tentokrát žádnou číselnou hodnotu nenastavujeme: žádáme o **obě pohlaví uvnitř každé barevné formy**, tedy o čtyři skutečně zastoupené kombinace.

```{r}
#| label: extra-emmeans-ctyri-skupiny
#| echo: true
# Ze stejného modelu získáme odhad pro každou kombinaci formy a pohlaví
prumery_formy <-
  emmeans::emmeans(
    object = mod_formy,
    specs = ~ pohlavi | forma
  )
prumery_formy
```

Sloupec `emmean` ukazuje stejné délky jako tabulka předpovědí z `predict()`. Protože model obsahuje obě skupiny i jejich interakci, tyto čtyři odhady jsou zde také rovny pozorovaným průměrům čtyř skupin. Intervaly `lower.CL` a `upper.CL` se vztahují ke **skupinovým průměrům**. Rozdíl **samec minus samice** pro každou formu vypočteme přímo:

```{r}
#| label: extra-emmeans-rozdily-formy
#| echo: true
# Opačné pořadí dvojice dává rozdíl samec minus samice
rozdily_formy <-
  emmeans::contrast(
    object = prumery_formy,
    method = "revpairwise"
  )
confint(object = rozdily_formy)
```

Řádek modré formy odpovídá rozdílu `r format_value_cz(rozdil_pohlavi_modra, digits = 2, signed = TRUE)` mm a řádek oranžové `r format_value_cz(rozdil_pohlavi_oranzova, digits = 2, signed = TRUE)` mm. Interval u každého řádku se týká **rozdílu pohlaví v dané formě**; interval pro změnu tohoto rozdílu mezi formami jsme získali výše z interakčního koeficientu. Tyto tři otázky mají odlišné intervaly, i když vycházejí z téhož modelu.
:::

## Co o jednotlivých krabech čtyři průměry neříkají?

Model čtyř skupin přesně vystihuje jejich `r render_glossary_term("prumer", display = "průměry")`. To však neznamená, že stejně dobře vystihuje každého kraba: jednotlivé délky se od skupinových průměrů výrazně liší. Zkontrolujeme proto `r render_glossary_term("residuum", display = "residua")` uvnitř skupin; jejich nulový průměr je zde vlastností odhadu skupinových průměrů, nikoli důkazem bezchybného modelu. V houslovém grafu (anglicky *violin plot*) širší část naznačuje, kde se odchylky častěji vyskytují, a každý bod ukazuje skutečné měření jednoho kraba.

<!-- TODO: add slug houslovy-graf to pojmy.yaml -->

```{r}
#| label: residua-ctyr-skupin
#| echo: false
#| fig-alt: "Houslové grafy residuí délky krunýře pro samice a samce v modré a oranžové formě; ve všech čtyřech skupinách leží body nad i pod přerušovanou nulou a rozpětí residuí je velké."
data_formy$residuum_mm <- residuals(object = mod_formy)

ggplot2::ggplot(
  data = data_formy,
  mapping = ggplot2::aes(
    x = pohlavi,
    y = residuum_mm,
    fill = pohlavi
  )
) +
  ggplot2::geom_violin(
    alpha = 0.35,
    width = 0.75,
    colour = biostat_cols[["grey_olive"]]
  ) +
  ggplot2::geom_point(
    position = ggplot2::position_jitter(
      width = 0.11,
      height = 0,
      seed = 900723
    ),
    alpha = 0.65,
    size = 1.3,
    colour = biostat_cols[["graphite"]]
  ) +
  ggplot2::geom_hline(
    yintercept = 0,
    linetype = "dashed",
    colour = biostat_cols[["graphite"]]
  ) +
  ggplot2::facet_wrap(
    facets = ggplot2::vars(forma),
    labeller = ggplot2::labeller(forma = vec_popisky_forem)
  ) +
  ggplot2::scale_fill_manual(values = vec_barvy_pohlavi, guide = "none") +
  ggplot2::labs(
    x = "Pohlaví",
    y = "Residuum: krab − průměr skupiny (mm)"
  ) +
  theme_biostat(base_size = 14)
```

Ve všech skupinách jsou délky na obou stranách odhadu. Rozpětí a krajní hodnoty se mezi skupinami mohou lišit; graf proto nečtěte jako důkaz, že čtyři průměry popisují jednotlivce přesně. Bez popisu způsobu výběru krabů také nelze rozdíly automaticky zobecnit na všechny kraby ani z nich vyvozovat účinek pohlaví nebo barevné formy.

## Přibývá objem dřeva s průměrem kmene stejně rychle u nízkých i vysokých stromů?

V obou příkladech s kraby byl aspoň jeden `r render_glossary_term("prediktor", display = "prediktor")` skupinou. Interakce však může spojovat i dva `r render_glossary_term("numericka-promenna", display = "číselné")` prediktory. Lesníci odhadují, kolik dřeva je v kmeni, z jeho průměru a z výšky stromu. Podíváme se na měření pokácených stromů střemchy pozdní (*Prunus serotina*), severoamerického stromu, který se v českých lesích šíří jako invazní druh.

::: {.box-question title="Co čekáte?"}
Představte si dvojici nízkých stromů stejné výšky, z nichž jeden má kmen o 1 cm tlustší, a stejnou dvojici vysokých stromů. O kolik více dřeva má tlustší strom v každé dvojici? Bude tento rozdíl u nízkých i vysokých stromů stejný?
:::

::: {.box-result collapse="true" title="Možná odpověď"}
Pravděpodobně ne. Centimetr průměru navíc znamená silnější vrstvu dřeva po celé délce kmene, takže delší kmen jí obsahuje více. Toto očekávání vychází z tvaru kmene, ne z fyziologie stromu. Pokud toto očekávání platí, `r render_glossary_term("sklon", display = "sklon")` vztahu objemu a průměru bude záviset na výšce.
:::

Soubor [stromy.csv](https://cuni-natur-biostatistics.github.io/L08/current/data/stromy.csv) je součástí materiálů lekce; uložte jej do stejné složky `data` jako soubor s kraby. Hodnoty jsou v původních amerických jednotkách: průměr kmene v palcích, výška ve stopách a objem dřeva v krychlových stopách. Proto je hned převedeme:

```{r}
#| label: nacist-stromy
#| echo: true
# Načteme tabulku pokácených stromů
data_stromy <-
  read.csv(file = "data/stromy.csv")

# Převedeme jednotky na cm, m a m³
data_stromy$prumer_cm <-
  data_stromy$Girth * 2.54 # 1 palec = 2,54 cm
data_stromy$vyska_m <-
  data_stromy$Height * 0.3048 # 1 stopa = 0,3048 m
data_stromy$objem_m3 <-
  data_stromy$Volume * 0.0283168 # 1 krychlová stopa = 0,0283168 m³

# Zkontrolujeme typy sloupců, chybějící hodnoty a rozsahy
str(object = data_stromy)
colSums(x = is.na(data_stromy))
summary(object = data_stromy[, c("prumer_cm", "vyska_m", "objem_m3")])
```

```{r}
#| label: popis-stromu
#| echo: false
korelace_prumer_vyska <-
  cor(x = data_stromy$prumer_cm, y = data_stromy$vyska_m)
```

Jeden řádek představuje **jeden pokácený strom**. Sloupec `Girth` (anglicky obvod) je podle dokumentace ve skutečnosti průměr kmene měřený asi 1,37 m nad zemí. V souboru je `r nrow(data_stromy)` stromů a chybí `r sum(is.na(data_stromy))` hodnot. Průměr kmene sahá od `r format_value_cz(min(data_stromy$prumer_cm), digits = 1)` do `r format_value_cz(max(data_stromy$prumer_cm), digits = 1)` cm, výška od `r format_value_cz(min(data_stromy$vyska_m), digits = 1)` do `r format_value_cz(max(data_stromy$vyska_m), digits = 1)` m a objem dřeva od `r format_value_cz(min(data_stromy$objem_m3), digits = 2)` do `r format_value_cz(max(data_stromy$objem_m3), digits = 2)` m³. Vyšší stromy mívají i tlustší kmen (`r render_glossary_term("korelace", display = "korelace")` průměru a výšky je `r format_value_cz(korelace_prumer_vyska, digits = 2)`); oba prediktory tedy nenesou úplně nezávislou informaci, podobně jako dvě teploty v minulé lekci.

Stromů je málo a dokumentace neuvádí, kde a jak byly vybrány. Závěry proto vztahujeme k těmto `r nrow(data_stromy)` stromům a k pozorovaným rozsahům průměru a výšky.

::: {.box-data title="Doplňující: odkud pocházejí data o stromech" collapse="true"}
Data pocházejí z knihy Meyer, H. A. (1953): *Forest Mensuration*, Penns Valley Publishers, převzaté do Ryan, T. A., Joiner, B. L. a Ryan, B. F. (1976): *The Minitab Student Handbook*, Duxbury Press. Tabulku šíří základní balíček R `datasets` pod licencí GPL-2 nebo GPL-3; soubor `stromy.csv` je její nezměněná kopie. Stejnou tabulku si můžete vytvořit sami:

```{r}
#| label: extra-priprava-stromy
#| echo: true
#| eval: false
# Tabulka je součástí základní instalace R
data_stromy_zdroj <-
  as.data.frame(x = datasets::trees)

# Kontrola: 31 stromů, 3 sloupce a žádné chybějící hodnoty
dim(x = data_stromy_zdroj)
sum(is.na(x = data_stromy_zdroj))

# Uložíme tabulku do složky data jako soubor CSV
write.csv(
  x = data_stromy_zdroj,
  file = "data/stromy.csv",
  row.names = FALSE
)
```
:::

## Vystihne společný sklon nízké i vysoké stromy?

U krabů jsme kreslili jednu přímku pro každé pohlaví. Výška je však číslo a žádné skupiny nemá. Abychom viděli, jak vypadá vztah objemu a průměru u různě vysokých stromů, rozdělíme stromy **jen pro kreslení** do čtyř výškových pásem se zhruba stejným počtem stromů. Každý panel ukazuje stromy jednoho pásma. Do modelu pásma nevstupují: model dál pracuje s výškou jako s číslem.

Začneme známým aditivním modelem z minulé lekce, který má dva číselné `r render_glossary_term("prediktor", display = "prediktory")` spojené znakem `+`:

```{r}
#| label: aditivni-model-stromy
#| echo: true
# Aditivní model: jeden společný sklon průměru pro všechny výšky
mod_aditivni <-
  lm(
    formula = objem_m3 ~ prumer_cm + vyska_m,
    data = data_stromy
  )
coef(object = mod_aditivni)
```

```{r}
#| label: pasma-stromu
#| echo: false
res_pasma_stromy <- build_tree_bands(data_stromy = data_stromy)
data_pasma_souhrn <- res_pasma_stromy$summary
sklon_aditivni <- unname(coef(mod_aditivni)["prumer_cm"])
```

Přímka v každém panelu je předpověď modelu pro **střední výšku daného pásma** (`r render_glossary_term("median", display = "medián")`: `r paste(format_value_cz(data_pasma_souhrn$vyska_median_m, digits = 1), collapse = "; ")` m) a vede jen přes průměry kmenů pozorované v tomto pásmu. V pásmech je `r paste(data_pasma_souhrn$n, collapse = ", ")` stromů.

```{r}
#| label: obrazek-pasma-aditivni
#| echo: false
#| fig-width: 9
#| fig-height: 3.8
#| fig-alt: "Čtyři panely vedle sebe, jeden pro každé výškové pásmo stromů od nejnižšího vlevo po nejvyšší vpravo; v každém je objem dřeva proti průměru kmene a fialová přímka aditivního modelu se stejným sklonem ve všech panelech. V nejnižším pásmu leží tenké kmeny nad přímkou a tlusté pod ní, v nejvyšším pásmu naopak."
plot_tree_panels(
  res_pasma = res_pasma_stromy,
  mod = mod_aditivni,
  base_size = 13
)
```

::: {.box-question title="Co v panelech vidíte?"}
Mají přímky ve všech panelech stejný `r render_glossary_term("sklon", display = "sklon")`? Leží body nejnižších a nejvyšších stromů kolem svých přímek stejně?
:::

::: {.box-result collapse="true" title="Možná odpověď"}
Přímky jsou rovnoběžné: aditivní model dává všem výškám společný sklon `r format_value_cz(sklon_aditivni, digits = 3)` m³ na 1 cm průměru. V nejnižším pásmu však tenké kmeny leží nad přímkou a tlusté pod ní, takže body stoupají pozvolněji než přímka. V nejvyšším pásmu je to naopak a body stoupají strměji. Společný sklon tedy zřejmě nevystihuje nízké ani vysoké stromy.
:::

## Model s interakcí dovolí, aby se sklon průměru měnil s výškou

Stejně jako u krabů zapíšeme `r render_glossary_term("interakce-predikatoru", display = "interakci")` hvězdičkou. Rozdíl je jen v tom, že druhý `r render_glossary_term("prediktor", display = "prediktor")` je číslo: interakční člen `prumer_cm:vyska_m` dovolí, aby se **`r render_glossary_term("sklon", display = "sklon")` průměru plynule měnil s výškou stromu**.

```{r}
#| label: model-stromy
#| echo: true
# Hvězdička dovolí, aby sklon průměru závisel na výšce
mod_stromy <-
  lm(
    formula = objem_m3 ~ prumer_cm * vyska_m,
    data = data_stromy
  )
coef(object = mod_stromy)
```

Model má opět čtyři koeficienty. Jejich význam rozebereme v další části; nejprve se podíváme, jaké přímky z nich vzniknou.

Na stejné panely nyní nakreslíme přímky modelu s interakcí. Šedá přerušovaná čára opakuje v každém panelu přímku nejnižšího pásma, aby bylo vidět, jak se sklon s výškou mění.

```{r}
#| label: obrazek-pasma-interakce
#| echo: false
#| fig-width: 9
#| fig-height: 3.8
#| fig-alt: "Stejné čtyři panely stromů s fialovými přímkami modelu s interakcí a šedou přerušovanou přímkou nejnižšího pásma v každém panelu; od nejnižšího k nejvyššímu pásmu je fialová přímka stále strmější a vzdaluje se od šedé."
plot_tree_panels(
  res_pasma = res_pasma_stromy,
  mod = mod_stromy,
  flag_reference = TRUE,
  base_size = 13
)
```

```{r}
#| label: sklony-stromy
#| echo: false
vec_koeficienty_stromy <- coef(mod_stromy)
koef_prumer <- unname(vec_koeficienty_stromy["prumer_cm"])
koef_vyska <- unname(vec_koeficienty_stromy["vyska_m"])
koef_interakce_stromy <- unname(vec_koeficienty_stromy["prumer_cm:vyska_m"])

# Teaching input: three heights inside the observed range
vec_vysky_sklon_m <- c(21, 23, 25)
vec_sklony_prumeru <- koef_prumer + koef_interakce_stromy * vec_vysky_sklon_m
```

Jak strmá je přímka při určité výšce? Sklon průměru složíme ze dvou koeficientů z výpisu `coef()`, z řádku `prumer_cm` a z interakčního řádku `prumer_cm:vyska_m`:

$$
\begin{aligned}
\text{sklon průměru při zvolené výšce}
&= \text{koeficient průměru} \\
&\quad + \text{interakční koeficient} \times \text{zvolená výška}
\end{aligned}
$$

Při výšce `r vec_vysky_sklon_m[[1]]` m dosadíme:

$$
\text{sklon průměru} = `r format_value_cz(koef_prumer, digits = 4, math = TRUE)` + `r format_value_cz(koef_interakce_stromy, digits = 5, math = TRUE)` \cdot `r vec_vysky_sklon_m[[1]]` \approx `r format_value_cz(vec_sklony_prumeru[[1]], digits = 3, math = TRUE)`\,\text{m}^3 \text{ na 1 cm}
$$

Stručně $\beta_1 + \beta_3 v$, kde $v$ je zvolená výška v metrech. Koeficienty číslujeme v pořadí výpisu: $\beta_0$ je intercept, $\beta_1$ řádek `prumer_cm`, $\beta_2$ řádek `vyska_m` a $\beta_3$ řádek `prumer_cm:vyska_m`. Stejně dostaneme při výškách `r paste(vec_vysky_sklon_m, collapse = ", ")` m sklony `r paste(format_value_cz(vec_sklony_prumeru, digits = 3), collapse = "; ")` m³ na 1 cm. S každým metrem výšky navíc přibude ke sklonu `r format_value_cz(koef_interakce_stromy, digits = 4)` m³ na 1 cm; to je hodnota interakčního koeficientu.

Stejné sklony ukážou i přímé předpovědi. Zvolíme dva průměry a dvě výšky, které leží uvnitř pozorovaných dat, a pro každou kombinaci spočteme odhad objemu:

```{r}
#| label: predpovedi-stromy
#| echo: true
# Zvolené průměry a výšky leží uvnitř pozorovaných rozsahů
prumery_porovnani_cm <- c(28, 36)
vysky_porovnani_m <- c(21, 25)

# Každá kombinace průměru a výšky dostane jeden řádek
data_porovnani_stromy <-
  expand.grid(
    prumer_cm = prumery_porovnani_cm,
    vyska_m = vysky_porovnani_m
  )
data_porovnani_stromy$objem_odhad_m3 <-
  predict(
    object = mod_stromy,
    newdata = data_porovnani_stromy
  )
data_porovnani_stromy
```

```{r}
#| label: rozdily-predpovedi-stromy
#| echo: false
vec_odhad_stromy <- data_porovnani_stromy$objem_odhad_m3
prirustek_nizky <-
  vec_odhad_stromy[data_porovnani_stromy$vyska_m == vysky_porovnani_m[[1]]] |>
  diff()
prirustek_vysoky <-
  vec_odhad_stromy[data_porovnani_stromy$vyska_m == vysky_porovnani_m[[2]]] |>
  diff()
krok_prumeru_cm <- diff(prumery_porovnani_cm)
```

::: {.box-result title="Stejný rozdíl v průměru, různý rozdíl v objemu"}
Kmen o průměru `r prumery_porovnani_cm[[2]]` cm místo `r prumery_porovnani_cm[[1]]` cm znamená podle modelu u stromu vysokého `r vysky_porovnani_m[[1]]` m o `r format_value_cz(prirustek_nizky, digits = 2)` m³ dřeva více, u stromu vysokého `r vysky_porovnani_m[[2]]` m o `r format_value_cz(prirustek_vysoky, digits = 2)` m³ více. Po vydělení `r krok_prumeru_cm` cm dostaneme právě sklony průměru při těchto výškách.
:::

## Co znamená koeficient průměru, když se sklon mění s výškou?

Vraťme se ke čtyřem koeficientům z výpisu `coef(object = mod_stromy)`:

- **`(Intercept)`** je `r render_glossary_term("odhad", display = "odhad")` objemu stromu s průměrem 0 cm a výškou 0 m. Takový strom neexistuje; číslo samo o sobě nevykládáme.
- **`prumer_cm`** je `r render_glossary_term("sklon", display = "sklon")` průměru **při výšce 0 m**: `r format_value_cz(koef_prumer, digits = 3, signed = TRUE)` m³ na 1 cm. Záporná hodnota neznamená, že by tlustší kmen měl méně dřeva; platí pro výšku, kterou žádný změřený strom nemá.
- **`vyska_m`** je obdobně sklon výšky **při průměru 0 cm**: `r format_value_cz(koef_vyska, digits = 3, signed = TRUE)` m³ na 1 m. Ani tento koeficient sám nevykládáme.
- **`prumer_cm:vyska_m`** je `r render_glossary_term("interakce-predikatoru", display = "interakční člen")`: o `r format_value_cz(koef_interakce_stromy, digits = 4)` m³ na 1 cm se zvětší sklon průměru, když je strom o 1 m vyšší.

::: {.box-warning title="Koeficient průměru neplatí při každé výšce"}
U dvou číselných `r render_glossary_term("prediktor", display = "prediktorů")` leží hodnota 0 často daleko mimo data. Hlavní koeficient pak popisuje situaci, kterou jsme nepozorovali. Sklon jednoho prediktoru proto vždy uvádíme **při zvolené hodnotě druhého** uvnitř pozorovaného rozsahu.
:::

```{r}
#| label: sklony-vysky-stromy
#| echo: false
vec_sklony_vysky <- koef_vyska + koef_interakce_stromy * prumery_porovnani_cm
interval_interakce_stromy <-
  confint(mod_stromy)["prumer_cm:vyska_m", ]
```

Interakci lze číst i obráceně. Stejný koeficient říká, o kolik se zvětší sklon výšky (m³ dřeva na 1 m výšky), když je kmen o 1 cm tlustší. U stromu s průměrem `r prumery_porovnani_cm[[1]]` cm přibude podle modelu na každý metr výšky `r format_value_cz(vec_sklony_vysky[[1]], digits = 3)` m³ dřeva, u stromu s průměrem `r prumery_porovnani_cm[[2]]` cm `r format_value_cz(vec_sklony_vysky[[2]], digits = 3)` m³. Obě věty popisují tutéž vlastnost modelu.

Jak přesně interakční koeficient známe, ukáže 95% `r render_glossary_term("interval-spolehlivosti", display = "interval spolehlivosti")`:

```{r}
#| label: interval-interakce-stromy
#| echo: true
# Interval patří změně sklonu průměru na 1 m výšky
confint(
  object = mod_stromy,
  parm = "prumer_cm:vyska_m"
)
```

Interval pro změnu sklonu průměru (m³ na 1 cm průměru) při zvýšení stromu o 1 m sahá od `r format_value_cz(interval_interakce_stromy[[1]], digits = 4)` do `r format_value_cz(interval_interakce_stromy[[2]], digits = 4)`. `r if (interval_interakce_stromy[[1]] > 0) "Celý leží nad nulou: při předpokladech modelu tedy data podporují strmější vztah objemu a průměru u vyšších stromů." else if (interval_interakce_stromy[[2]] < 0) "Celý leží pod nulou: při předpokladech modelu tedy data podporují pozvolnější vztah objemu a průměru u vyšších stromů." else "Zahrnuje i nulu, takže směr změny sklonu s výškou zůstává podle tohoto modelu nejistý."` Platí to pro `r nrow(data_stromy)` změřených stromů a v pozorovaném rozsahu průměrů a výšek.

V tomto krátkém příkladu jsme nekontrolovali `r render_glossary_term("residuum", display = "residua")` ani vliv jednotlivých stromů, jako jsme to udělali u krabů; `r nrow(data_stromy)` stromů je pro takovou kontrolu malý soubor. Objem kmene navíc roste zhruba s druhou mocninou průměru, takže část rozdílu mezi panely by mohl popsat i zakřivený vztah. Přímky zde proto čteme jako jednoduchý popis v pozorovaném rozsahu, ne jako jediný možný model.

## Kdy musíme při výkladu jednoho prediktoru uvést hodnotu druhého?

První příklad vyjádřil `r render_glossary_term("interakce-predikatoru", display = "interakci")` jako **změnu rozdílu v šířce mezi pohlavími při různých délkách**: dva `r render_glossary_term("sklon", display = "sklony")` nemusí být stejné. Druhý příklad ji vyjádřil jako **změnu rozdílu v délce mezi pohlavími u dvou barevných forem**: dvě srovnání samce se samicí nemusí dát stejný výsledek. Třetí příklad ji vyjádřil jako **plynulou změnu sklonu průměru s výškou stromu**: s každým metrem výšky je přímka o něco strmější. Ve všech třech případech musíme při výkladu jednoho `r render_glossary_term("prediktor", display = "prediktoru")` uvést hodnotu nebo úroveň druhého: ptáme se **„při jaké délce?“**, **„v jaké barevné formě?“** nebo **„při jaké výšce?“**.

::: {.box-data title="Doplňující: interakce závisí na měřítku odezvy" collapse="true"}

V příkladech s kraby popisujeme `r render_glossary_term("odezva", display = "odezvu")` v milimetrech, takže interakce znamená změnu **rozdílu** v milimetrech. Jiná biologická otázka se může ptát na poměr: například zda jsou samci v obou formách o stejných 10 % delší než samice. Stejný procentní poměr může při různých `r render_glossary_term("prumer", display = "průměrech")` vytvářet různé rozdíly v milimetrech.

```{r}
#| label: stromy-logaritmy
#| echo: false
mod_stromy_log <-
  lm(
    formula = log(objem_m3) ~ log(prumer_cm) * log(vyska_m),
    data = data_stromy
  )
interval_interakce_log <-
  confint(mod_stromy_log)["log(prumer_cm):log(vyska_m)", ]
```

Poměry a součiny se na logaritmickém měřítku mění v rozdíly a součty. Příklad se stromy to ukazuje na skutečných datech. Objem kmene je zhruba úměrný druhé mocnině průměru krát výšce, tedy *součinu*. Když objem, průměr i výšku nahradíme jejich logaritmy (`r render_glossary_term("transformace", display = "transformace")`, ke které se vrátíme v lekci o tvaru rozdělení), změní se tento součin v součet a model bez interakce by mohl stačit. Interakční koeficient takového modelu má 95% `r render_glossary_term("interval-spolehlivosti", display = "interval spolehlivosti")` od `r format_value_cz(interval_interakce_log[[1]], digits = 2)` do `r format_value_cz(interval_interakce_log[[2]], digits = 2)`. Interval je široký a zahrnuje nulu, takže na logaritmickém měřítku data interakci nedokládají, ačkoli v krychlových metrech ji ukazují zřetelně.

Tvrzení „interakce existuje“ proto není úplné bez uvedení měřítka modelu. V lekci o binárních a početních odezvách se setkáme s modely, které používají jiné měřítko než přímé milimetry; interakční koeficient se pak nejprve interpretuje na tomto modelovém měřítku a potom pomocí předpovědí v biologicky srozumitelných jednotkách.
:::

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

Vraťme se k hlavní otázce: **Je rozdíl v zadní šířce krunýře mezi stejně dlouhou samicí a samcem stejný u menších i větších krabů?**

U změřených krabů modré formy není. Zadní šířka roste s délkou u obou pohlaví, ale u samic rychleji: o `r format_value_cz(sklon_samice, digits = 3)` mm na 1 mm délky, u samců o `r format_value_cz(sklon_samec, digits = 3)` mm. Rozdíl samec minus samice se proto mění o `r format_value_cz(rozdil_sklonu_samec_minus_samice, digits = 3, signed = TRUE)` mm na každý milimetr délky (95% `r render_glossary_term("interval-spolehlivosti", display = "interval spolehlivosti")` `r format_value_cz(interval_rozdilu_sklonu[[1]], digits = 3)` až `r format_value_cz(interval_rozdilu_sklonu[[2]], digits = 3)` mm/mm). Při délce `r delka_porovnani_mm` mm má samec podle modelu zadní šířku o `r format_value_cz(abs(rozdil_porovnani), digits = 2)` mm menší než stejně dlouhá samice. Graf residuí neukázal zakřivení, které by dvě přímky přehlédly.

Výsledek odpovídá biologickému očekávání z úvodu, že samice mají s rostoucí velikostí relativně širší zadní část krunýře. Data však jsou pozorovací a není doloženo, že krabi byli vybráni náhodně: popisují souvislost u těchto `r nrow(data_krabi)` krabů a v pozorovaném rozsahu délek, nikoli příčinu. Druhý příklad ukázal, že stejná logika platí pro dvě skupinové proměnné: rozdíl v délce mezi pohlavími byl u oranžové formy o `r format_value_cz(abs(rozdil_rozdilu_formy), digits = 2)` mm menší než u modré, ale `r render_glossary_term("prumer", display = "průměrné")` délky závisejí i na tom, které kraby autoři vybrali.

Třetí příklad přenesl stejnou myšlenku na dva číselné `r render_glossary_term("prediktor", display = "prediktory")`. U `r nrow(data_stromy)` změřených stromů byl vztah objemu a průměru tím strmější, čím byl strom vyšší: při `r vysky_porovnani_m[[1]]` m přidal kmen o `r krok_prumeru_cm` cm tlustší podle modelu `r format_value_cz(prirustek_nizky, digits = 2)` m³ dřeva, při `r vysky_porovnani_m[[2]]` m `r format_value_cz(prirustek_vysoky, digits = 2)` m³. Na původním měřítku objemu to odpovídá očekávání z tvaru kmene; platí to jen v pozorovaném rozsahu průměrů a výšek.

## Co bude následovat dál

Aditivní model a model s interakcí jsou dva kandidátní popisy téhož vztahu. Model s interakcí je pružnější, ale má více parametrů a jeho koeficienty se vykládají složitěji. V následující lekci se budeme ptát, jak mezi několika předem biologicky zdůvodněnými modely vybrat ten, který data popisuje dostatečně dobře bez zbytečné složitosti.

## Shrnutí

- `r render_glossary_term("interakce-predikatoru", display = "Interakce prediktorů")` znamená, že vztah jednoho `r render_glossary_term("prediktor", display = "prediktoru")` k `r render_glossary_term("odezva", display = "odezvě")` závisí na hodnotě nebo úrovni druhého; potřebu interakce má podpořit biologická otázka i graf.
- `x * y` je zkratka pro `x + y + x:y`; interakční člen u číselného a skupinového prediktoru popisuje rozdíl `r render_glossary_term("sklon", display = "sklonů")`.
- Hlavní koeficienty v modelu s interakcí jsou podmíněné: platí při nulové hodnotě druhého prediktoru nebo v referenční skupině. Skupiny proto porovnáváme pomocí předpovědí při zvolené hodnotě.
- U dvou skupinových prediktorů je interakce rozdíl dvou rozdílů.
- U dvou číselných prediktorů interakce znamená, že se sklon jednoho plynule mění s hodnotou druhého; sklon proto uvádíme při zvolených hodnotách uvnitř dat.
- Závěr omezujeme na pozorovaný rozsah dat, na způsob jejich získání a na souvislost, nikoli příčinu.

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

Uveďte biologickou otázku, u které byste čekali, že vztah dvou znaků závisí na skupině nebo na hodnotě jiného znaku. Jak by vypadal graf?

::: {.box-result collapse="true" title="Možné odpovědi"}
Například vztah hmotnosti a délky těla u ryb z teplé a studené vody: v teplejší vodě mohou ryby s rostoucí délkou přibírat rychleji, takže přímky by se rozbíhaly. Nebo vztah dávky hnojiva a výnosu u dvou odrůd plodiny, z nichž jedna na hnojivo reaguje silněji. V obou případech by `r render_glossary_term("bodovy-graf", display = "bodový graf")` ukázal dvě skupiny bodů s různým sklonem a aditivní model s rovnoběžnými přímkami by takový rozdíl nezachytil.

Místo skupiny může vztah měnit i další číselný znak: například růst semenáčků s teplotou může být strmější ve vlhčí půdě. Graf by pak měl několik panelů pro různou vlhkost půdy a přímka růstu proti teplotě by byla s rostoucí vlhkostí strmější.
:::

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

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

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