--- title: "Když prediktor není číslo" subtitle: "Od rozdílu průměrů k analýze rozptylu" lang: cs filters: - ../theme/semantic-boxes.lua format: html: toc: true toc-depth: 3 toc-title: "Obsah" number-sections: true code-fold: false code-summary: "Zobrazit kód v R" embed-resources: true theme: [flatly, ../theme/skripta_theme.scss] highlight-style: github include-in-header: - "../theme/fonts-include.html" - "../theme/skripta-favicon.html" typst: toc: true toc-depth: 3 number-sections: true output-file: "skripta_raw.pdf" include-in-header: - "skripta_theme.typ" - "skripta-logo.typ" date: last-modified date: 2026/08/06 date-format: long date-modified: last-modified --- ```{r} #| label: nastaveni #| include: false if (!requireNamespace("here", quietly = TRUE)) install.packages("here") if (!requireNamespace("readr", quietly = TRUE)) install.packages("readr") if (!requireNamespace("dplyr", quietly = TRUE)) install.packages("dplyr") if (!requireNamespace("tidyr", quietly = TRUE)) install.packages("tidyr") if (!requireNamespace("ggplot2", quietly = TRUE)) install.packages("ggplot2") if (!requireNamespace("emmeans", quietly = TRUE)) install.packages("emmeans") if (!requireNamespace("ggbeeswarm", quietly = TRUE)) { install.packages("ggbeeswarm") } if (!requireNamespace("glossary", quietly = TRUE)) { install.packages("glossary") } here::i_am("Learning_materials/skripta.qmd") # Visible code reads data/palmer_penguins.csv relative to the project root, # exactly as students do in their own R project. knitr::opts_knit$set(root.dir = here::here()) # Hidden figures use Czech decimal commas on their axes. Visible chunks keep # R's default output so that it matches what students see in R. knitr::knit_hooks$set( desetinna_carka = function(before) { if (before) { options(OutDec = ",") } else { options(OutDec = ".") } } ) knitr::opts_hooks$set( echo = function(options) { if (isFALSE(options$echo)) { options$desetinna_carka <- TRUE } options } ) source(here::here("R", "set_r_theme.R")) source(here::here("R", "Functions", "format_cz.R")) source(here::here("R", "Functions", "format_cz_math.R")) source(here::here("R", "Functions", "format_p_cz.R")) # Species colours follow the palmerpenguins dataset's own established palette # (https://allisonhorst.github.io/palmerpenguins/), kept by decision of # O. Mottl (2026-10-06, L02) as a documented exception to the _brand # categorical palette. Species are groups here, not semantic roles. barvy_druhu <- c( "Adelie" = "darkorange", "Chinstrap" = "mediumpurple", "Gentoo" = "cyan4" ) popisky_druhu <- c( "Adelie" = "tučňák kroužkový (Adélie)", "Chinstrap" = "tučňák uzdičkový (Chinstrap)", "Gentoo" = "tučňák oslí (Gentoo)" ) 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 #| message: false #| results: asis if (knitr::is_html_output()) { glossary::glossary_style( color = biostat_cols[["indigo_velvet"]], text_decoration = "underline", def_bg = "#FFFFFF", def_color = "#2E2E2E" ) cat( "", "", sep = "\n" ) } ``` ## Úvod ::: {.box-question title="Hlavní biologická otázka"} **Liší se `r render_glossary_term("prumer", display = "průměrná")` hmotnost těla mezi druhy tučňáků?** ::: V předchozí lekci jsme převáděli biologické tvrzení na `r render_glossary_term("nulova-hypoteza", display = "nulovou hypotézu")` o koeficientu známého modelu. `r render_glossary_term("prediktor", display = "Prediktorem")` byla `r render_glossary_term("numericka-promenna", display = "numerická proměnná")` (zeměpisná šířka) a otázka zněla, zda je `r render_glossary_term("sklon", display = "sklon")` slučitelný s nulou. V této lekci bude prediktorem `r render_glossary_term("kategoricka-promenna", display = "kategoriální proměnná")`: druh tučňáka. Způsob uvažování zůstává stejný. Výsledek začíná jednotlivými `r render_glossary_term("pozorovani", display = "pozorováními")`, pokračuje rozdílem v gramech a jeho `r render_glossary_term("nejistota", display = "nejistotou")` a testová otázka jej pouze doplňuje. ## Výsledky učení Po prostudování této lekce dokážete: - fitovat `lm(y ~ skupina)` pro dvě i více skupin a přečíst jeho koeficienty; - interpretovat `r render_glossary_term("intercept", display = "intercept")` jako `r render_glossary_term("prumer", display = "průměr")` referenční skupiny a ostatní koeficienty jako rozdíly vůči této skupině; - vysvětlit `r render_glossary_term("t-test", display = "dvouvýběrový t-test")` a `r render_glossary_term("anova", display = "analýzu rozptylu (ANOVA)")` jako otázky položené `r render_glossary_term("linearni-model", display = "lineárnímu modelu")`; - rozlišit celkový F-test od `r render_glossary_term("post-hoc-test", display = "post-hoc porovnání")` dvojic a rozhodnout, kdy párová porovnání dávají smysl; - napsat závěr v biologicky smysluplných jednotkách: rozdíl, jeho `r render_glossary_term("nejistota", display = "nejistota")`, test a omezení studie. ## Data: jeden řádek je jeden změřený tučňák ```{r} #| label: nacist-data-tucnaku #| echo: false #| output: false #| message: false #| warning: false #| error: false data_penguins_raw <- readr::read_csv( file = here::here("data", "palmer_penguins.csv"), na = c("", "NA"), show_col_types = FALSE ) count_all <- nrow(data_penguins_raw) count_missing_mass <- data_penguins_raw |> dplyr::summarise( count = sum(is.na(hmotnost_tela_g)) ) |> dplyr::pull(count) data_penguins <- data_penguins_raw |> dplyr::filter(!is.na(hmotnost_tela_g)) |> dplyr::mutate( druh = factor( druh, levels = c("Adelie", "Chinstrap", "Gentoo") ) ) count_complete <- nrow(data_penguins) data_penguins_two <- data_penguins |> dplyr::filter(druh %in% c("Adelie", "Gentoo")) |> dplyr::mutate( druh = droplevels(druh) ) ``` Data obsahují tělesná měření dospělých tučňáků tří druhů z Palmerova souostroví v Antarktidě. Jeden řádek představuje jednoho změřeného jedince. Naše `r render_glossary_term("odezva", display = "odezva")` je hmotnost těla v gramech a `r render_glossary_term("kategoricka-promenna", display = "kategoriálním prediktorem")` je druh. ![Druhy tučňáků v datasetu Palmer Penguins a barvy, které pro ně budeme používat. Ilustrace: Allison Horst.](images/lter_penguins.png){fig-alt="Kreslené ilustrace tří druhů tučňáků: tučňák uzdičkový, tučňák oslí a tučňák kroužkový." width="85%"} Stejně jako v lekci o grafech a typech `r render_glossary_term("promenna", display = "proměnných")` budeme tučňáka kroužkového (Adélie) zobrazovat oranžově, tučňáka uzdičkového (Chinstrap) fialově a tučňáka oslího (Gentoo) modrozeleně. Tabulku načteme a známými funkcemi zjistíme, jaké hodnoty a datové typy obsahuje. `str()` ukáže strukturu tabulky a typ každé proměnné, `summary()` základní souhrny sloupců a `colSums(is.na())` počet chybějících hodnot v každém sloupci. ```{r} #| label: pripravit-data-viditelne #| echo: true #| eval: true # Načteme připravenou tabulku s jedním řádkem na tučňáka data_tucnaci <- read.csv( file = "data/palmer_penguins.csv", na.strings = c("", "NA") ) # Zkontrolujeme strukturu a datové typy str(data_tucnaci) # Prohlédneme si souhrny hodnot summary(data_tucnaci) # Spočítáme chybějící hodnoty v každém sloupci colSums(x = is.na(x = data_tucnaci)) ``` Tabulka má `r count_all` řádků. Sloupec `druh` je zatím text (`chr`) a u `r count_missing_mass` tučňáků chybí hmotnost. Hmotnost se pohybuje od `r format_cz(min(data_penguins$hmotnost_tela_g), 0)` do `r format_cz(max(data_penguins$hmotnost_tela_g), 0)` g, což jsou u dospělých tučňáků těchto druhů věrohodné hodnoty. Model hmotnosti nemůže použít řádek, ve kterém hmotnost chybí: bez naměřené `r render_glossary_term("odezva", display = "odezvy")` nelze spočítat `r render_glossary_term("residuum", display = "residuum")`. Ponecháme proto jen tučňáky se změřenou hmotností. Ostatní chybějící položky (například pohlaví) zatím nevadí, protože je v modelu nepoužíváme. Sloupec `druh` zároveň převedeme na **faktor**, tedy kategoriální proměnnou s pevně daným pořadím úrovní. Toto pořadí bude později důležité pro čtení koeficientů. ```{r} #| label: odstranit-chybejici-hmotnosti-viditelne #| echo: true #| eval: true # Ponecháme jen tučňáky se změřenou hmotností data_tucnaci <- data_tucnaci[ !is.na(data_tucnaci$hmotnost_tela_g), ] # Druh převedeme na faktor s pevným pořadím úrovní data_tucnaci$druh <- factor( x = data_tucnaci$druh, levels = c("Adelie", "Chinstrap", "Gentoo") ) # Ověříme počet ponechaných řádků a úrovně druhu nrow(x = data_tucnaci) levels(x = data_tucnaci$druh) ``` ```{r} #| label: ukazat-prvni-radky #| echo: false #| output: true #| message: false #| warning: false #| error: false data_penguins_raw |> dplyr::select( druh, ostrov, hmotnost_tela_g, pohlavi, rok ) |> dplyr::slice(1:6) |> knitr::kable( col.names = c( "Druh", "Ostrov", "Hmotnost těla (g)", "Pohlaví", "Rok" ), caption = "Prvních šest řádků připravené tabulky." ) ``` {{< pagebreak >}} ::: {.box-data title="Doplňující: získání a příprava dat" collapse="true"} Pro průchod lekcí stačí připravený soubor `data/palmer_penguins.csv`. Ten záměrně zachovává všech `r count_all` řádků, aby byly chybějící hodnoty při prvním prohlédnutí viditelné. Pokud si ale chcete celý postup zopakovat, následující kód provede celý sled kroků: **získá data → vybere a přejmenuje sloupce → připraví je pro model hmotnosti → uloží vlastní pracovní kopii**. Kód vynechá řádky bez hmotnosti hned po stažení, takže vlastní kopie bude mít jen `r count_complete` řádků. Druh v ní rovnou uloží jako faktor se stejným pořadím úrovní jako výše. Tento blok můžete přeskočit; další příklady používají tabulku `data_tucnaci` připravenou výše. ```{r} #| label: stahnout-a-pripravit-puvodni-data #| echo: true #| eval: false # Tento příkaz stačí spustit pouze jednou # install.packages("palmerpenguins") # Získáme data a prohlédneme chybějící hmotnosti data_puvodni <- palmerpenguins::penguins sum(is.na(data_puvodni$body_mass_g)) data_puvodni <- data_puvodni[!is.na(data_puvodni$body_mass_g), ] # Vybereme, přejmenujeme a upravíme sloupce pro lekci data_tucnaci <- data.frame( druh = factor( x = data_puvodni$species, levels = c("Adelie", "Chinstrap", "Gentoo") ), ostrov = as.character(data_puvodni$island), delka_zobaku_mm = data_puvodni$bill_length_mm, hloubka_zobaku_mm = data_puvodni$bill_depth_mm, delka_kridla_mm = data_puvodni$flipper_length_mm, hmotnost_tela_g = data_puvodni$body_mass_g, pohlavi = as.character(data_puvodni$sex), rok = data_puvodni$year ) # Uložíme si zpracovanou tabulku pro další práci write.csv( x = data_tucnaci, file = "data/palmer_penguins_vlastni_kopie.csv", row.names = FALSE, na = "" ) ``` ::: Balíček `palmerpenguins` obsahuje kopii původních tabulek z Environmental Data Initiative. [Dokumentace balíčku palmerpenguins](https://allisonhorst.github.io/palmerpenguins/articles/download.html) uvádí také cestu k neupraveným archivním souborům. Měření pocházejí z výzkumu tučňáků rodu *Pygoscelis* v letech 2007–2009. Data shromáždila Kristen Gorman ve spolupráci s Palmer Station Long Term Ecological Research a byla zveřejněna pod licencí CC0. Připravený soubor a jeho přesná provenance jsou popsány v `data/README.md`; biologickou studii uvádějí Gorman, Williams a Fraser (2014), DOI: [10.1371/journal.pone.0090081](https://doi.org/10.1371/journal.pone.0090081). ## Hmotnosti tučňáků kroužkových a oslích Pro první srovnání ponecháme jen dva druhy: tučňáka kroužkového (Adélie) a tučňáka oslího (Gentoo). V prvním grafu jsou jejich hmotnosti smíchané dohromady. Každý bod představuje jedno `r render_glossary_term("pozorovani", display = "pozorování")`, tedy jednoho změřeného jedince. ```{r} #| label: pripravit-souhrn-smisene-hmotnosti #| echo: false #| output: false #| message: false #| warning: false #| error: false data_summary_pooled <- data_penguins_two |> dplyr::summarise( n = dplyr::n(), prumer_g = mean(hmotnost_tela_g), median_g = median(hmotnost_tela_g), minimum_g = min(hmotnost_tela_g), maximum_g = max(hmotnost_tela_g), iqr_g = IQR(hmotnost_tela_g) ) ``` ```{r} #| label: zobrazit-smisene-hmotnosti #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Hmotnosti tučňáků dvou druhů předtím, než použijeme informaci o druhu. Všechny body jsou záměrně šedé." #| fig-alt: "Šedé body hmotností tučňáků kroužkových a oslích v jednom společném sloupci s obrysem rozdělení; hodnoty sahají přibližně od tří do šesti kilogramů." #| fig-width: 8 #| fig-height: 4.6 ggplot2::ggplot( data = data_penguins_two, mapping = ggplot2::aes( x = "všichni jedinci", y = hmotnost_tela_g ) ) + ggplot2::geom_violin( fill = biostat_cols[["grey_olive"]], alpha = 0.15, linewidth = 0.8, show.legend = FALSE ) + ggbeeswarm::geom_quasirandom( width = 0.25, alpha = 0.55, size = 1.7, color = biostat_cols[["grey_olive"]] ) + ggplot2::labs( x = NULL, y = "Hmotnost těla (g)" ) ``` Bez další `r render_glossary_term("promenna", display = "proměnné")` umíme popsat jen tento společný soubor. V lekci o popisu jedné proměnné jsme k tomu používali počet pozorování, minimum, maximum, `r render_glossary_term("prumer", display = "průměr")`, `r render_glossary_term("median", display = "medián")` a `r render_glossary_term("iqr", display = "mezikvartilové rozpětí (IQR)")`. Stejné souhrny nyní spočítáme v R: ```{r} #| label: souhrn-smisene-hmotnosti-viditelne #| echo: true #| eval: true # Pro první pohled vybereme dva druhy, ale zatím je nesrovnáváme data_dva_druhy <- data_tucnaci[ data_tucnaci$druh %in% c("Adelie", "Gentoo"), ] data_dva_druhy$druh <- droplevels(data_dva_druhy$druh) # Připomeneme si známé souhrny pro všechny hmotnosti dohromady c( pocet = length(data_dva_druhy$hmotnost_tela_g), prumer = mean(data_dva_druhy$hmotnost_tela_g), median = median(data_dva_druhy$hmotnost_tela_g), minimum = min(data_dva_druhy$hmotnost_tela_g), maximum = max(data_dva_druhy$hmotnost_tela_g), IQR = IQR(data_dva_druhy$hmotnost_tela_g) ) ``` Společný průměr `r format_cz(data_summary_pooled$prumer_g, 1)` g shrnuje střed šedého oblaku, ale zakrývá možnost, že se v něm mísí odlišné skupiny. Tuto možnost prověří proměnná `druh`. Proměnná `druh` je `r render_glossary_term("kategoricka-promenna", display = "kategoriální")`: její hodnoty nejsou měření na číselné škále, ale názvy kategorií. Funkce `table()` ukáže, kolik řádků (jedinců) patří do každé kategorie: ```{r} #| label: zobrazit-pocty-druhu #| echo: true #| eval: true table(data_tucnaci$druh) ``` V celé tabulce jsou tři druhy. Pro první model ponecháváme pouze Adélie a Gentoo, které jsme vybrali do objektu `data_dva_druhy`. ::: {.box-question title="Co by tvrdila nulová hypotéza?"} V předchozí lekci jsme testovali `r render_glossary_term("nulova-hypoteza", display = "nulovou hypotézu")` o koeficientu modelu. Zkuste vyslovit nulovou hypotézu pro průměrnou hmotnost Adélie a Gentoo. Jak by podle ní vypadal graf, ve kterém body obarvíme podle druhu? ::: ::: {.box-result collapse="true" title="Možná odpověď"} Nulová hypotéza tvrdí, že průměrná hmotnost obou druhů je v širší populaci stejná: rozdíl průměrů je 0 g. Body obou druhů by pak ležely kolem stejné výšky. Jednotlivci by se stále lišili, ale žádný druh by systematicky neležel výše. ::: ```{r} #| label: zobrazit-jednotlive-hmotnosti-podle-druhu #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Stejná měření po přidání kategoriální proměnné druh. Barva i vodorovná poloha nyní ukazují příslušnost ke skupině." #| fig-alt: "Hmotnosti rozdělené do dvou sloupců podle druhu: oranžové body tučňáků kroužkových leží převážně níže než modrozelené body tučňáků oslích, obě skupiny se částečně překrývají." #| fig-width: 8 #| fig-height: 5 ggplot2::ggplot( data = data_penguins_two, mapping = ggplot2::aes( x = druh, y = hmotnost_tela_g, color = druh ) ) + ggplot2::geom_violin( mapping = ggplot2::aes( fill = druh ), alpha = 0.15, linewidth = 0.8, show.legend = FALSE ) + ggbeeswarm::geom_quasirandom( width = 0.22, alpha = 0.65, size = 1.8, show.legend = FALSE ) + ggplot2::scale_color_manual(values = barvy_druhu) + ggplot2::scale_x_discrete(labels = popisky_druhu) + ggplot2::labs( x = NULL, y = "Hmotnost těla (g)" ) ``` Po rozdělení vidíme silný posun: body tučňáků oslích leží převážně výše. Zároveň se hodnoty uvnitř každého druhu liší a oba soubory se částečně překrývají. Druh nám tedy pomohl odhalit strukturu, kterou šedý graf skrýval. ### Průměr a `r render_glossary_term("smerodatna-odchylka", display = "směrodatná odchylka")` každého druhu ```{r} #| label: spocitat-souhrny-dvou-druhu #| echo: false #| output: false #| message: false #| warning: false #| error: false data_summary_two <- data_penguins_two |> dplyr::group_by(druh) |> dplyr::summarise( n = dplyr::n(), prumer_g = mean(hmotnost_tela_g), sd_g = sd(hmotnost_tela_g), .groups = "drop" ) |> dplyr::mutate( x = c(1, 2) ) mean_adelie <- data_summary_two |> dplyr::filter(druh == "Adelie") |> dplyr::pull(prumer_g) mean_gentoo <- data_summary_two |> dplyr::filter(druh == "Gentoo") |> dplyr::pull(prumer_g) difference_gentoo_adelie <- mean_gentoo - mean_adelie ``` ```{r} #| label: ukazat-souhrny-dvou-druhu #| echo: false #| output: true #| message: false #| warning: false #| error: false data_summary_two |> dplyr::transmute( Druh = popisky_druhu[as.character(druh)], `Počet jedinců` = n, `Průměr (g)` = round(prumer_g, digits = 1), `SD (g)` = round(sd_g, digits = 1) ) |> knitr::kable( caption = paste0( "Počet pozorování, průměr a směrodatná odchylka ", "hmotnosti ve dvou skupinách." ) ) ``` Průměr shrnuje polohu jedné skupiny, ale neříká, jak se liší jednotlivci uvnitř ní. Proto vedle něj zobrazíme `r render_glossary_term("smerodatna-odchylka", display = "směrodatnou odchylku")` (SD) i všechny body. Svislá úsečka kolem průměru sahá od průměru minus jedna SD po průměr plus jedna SD. ```{r} #| label: zobrazit-prumery-dvou-druhu #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Body jsou jednotliví tučňáci, silné vodorovné úsečky jsou průměry a tenké svislé úsečky ukazují průměr ± 1 SD." #| fig-alt: "Body hmotností obou druhů s vodorovnou úsečkou průměru a svislou úsečkou průměr plus minus jedna směrodatná odchylka; průměr tučňáků oslích leží výrazně výše." #| fig-width: 8 #| fig-height: 5 ggplot2::ggplot( data = data_penguins_two, mapping = ggplot2::aes( x = druh, y = hmotnost_tela_g, color = druh ) ) + ggplot2::geom_violin( mapping = ggplot2::aes( fill = druh ), alpha = 0.15, linewidth = 0.8, show.legend = FALSE ) + ggbeeswarm::geom_quasirandom( width = 0.22, alpha = 0.45, size = 1.7, show.legend = FALSE ) + ggplot2::geom_errorbar( data = data_summary_two, mapping = ggplot2::aes( x = druh, ymin = prumer_g - sd_g, ymax = prumer_g + sd_g, color = druh ), width = 0.16, linewidth = 0.9, inherit.aes = FALSE, show.legend = FALSE ) + ggplot2::geom_segment( data = data_summary_two, mapping = ggplot2::aes( x = x - 0.28, xend = x + 0.28, y = prumer_g, yend = prumer_g, color = druh ), linewidth = 2.2, inherit.aes = FALSE, show.legend = FALSE ) + ggplot2::scale_color_manual(values = barvy_druhu) + ggplot2::scale_x_discrete(labels = popisky_druhu) + ggplot2::labs( x = NULL, y = "Hmotnost těla (g)" ) ``` Průměr tučňáků oslích je o `r format_cz(difference_gentoo_adelie, 1)` g vyšší než průměr tučňáků kroužkových. To je hlavní věcný rozdíl, který bude náš model popisovat. ## Rozdíl `r render_glossary_term("prumer", display = "průměrů")` Dva průměry se mohou lišit o stejnou hodnotu, ale síla informace závisí také na `r render_glossary_term("variabilita", display = "variabilitě")` uvnitř skupin. Následující graf proto skrývá jednotlivé body a ponechává pouze tři součásti: dva průměry, jejich rozdíl a `r render_glossary_term("smerodatna-odchylka", display = "SD")` každé skupiny. ```{r} #| label: zobrazit-rozdil-prumeru-a-sd #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Silné vodorovné úsečky jsou průměry, svislé úsečky ukazují ± 1 SD a závorka vpravo měří rozdíl mezi průměry. Celkový rozsah dat zde záměrně nezobrazujeme." #| fig-alt: "Dva skupinové průměry hmotnosti s úsečkami plus minus jedna SD a svislou šipkou vpravo, která měří rozdíl mezi průměry v gramech." #| fig-width: 8 #| fig-height: 5 ggplot2::ggplot( data = data_summary_two, mapping = ggplot2::aes( x = x, y = prumer_g, color = druh ) ) + ggplot2::geom_errorbar( mapping = ggplot2::aes( ymin = prumer_g - sd_g, ymax = prumer_g + sd_g ), width = 0.14, linewidth = 1 ) + ggplot2::geom_segment( mapping = ggplot2::aes( x = x - 0.26, xend = x + 0.26, y = prumer_g, yend = prumer_g ), linewidth = 2.3 ) + ggplot2::annotate( geom = "segment", x = 2.52, xend = 2.52, y = mean_adelie, yend = mean_gentoo, color = biostat_cols[["graphite"]], linewidth = 1, arrow = grid::arrow(ends = "both", length = grid::unit(0.12, "cm")) ) + ggplot2::annotate( geom = "label", x = 2.57, y = mean(c(mean_adelie, mean_gentoo)), label = paste0(format_cz(difference_gentoo_adelie, 0), " g"), hjust = 0, color = biostat_cols[["graphite"]], fill = "white", linewidth = 0.2, size = 4.2 ) + ggplot2::scale_color_manual(values = barvy_druhu) + ggplot2::scale_x_continuous( breaks = c(1, 2), labels = popisky_druhu[c("Adelie", "Gentoo")], limits = c(0.65, 2.95) ) + ggplot2::labs( x = NULL, y = "Hmotnost těla (g)", color = "Druh" ) + ggplot2::theme(legend.position = "none") ``` Svislá závorka měří systematický rozdíl, který chceme vysvětlit druhem. Úsečky SD připomínají, že každý průměr zastupuje celý oblak různých jedinců. Rozdíl průměrů proto nikdy nebudeme číst bez informace o variabilitě kolem nich. ### Stejný rozdíl, jiná variabilita ```{r} #| label: pripravit-myslenkovy-pokus-variabilita #| echo: false #| output: false #| message: false #| warning: false #| error: false # Teaching inputs of a constructed thought experiment, identical to the # generator of fixed_gap_varying_sd.png in the presentation. vec_pattern <- as.numeric(scale(seq(-1, 1, length.out = 10))) prumer_skupiny_a <- 3500 rozdil_myslenkovy_pokus <- 1000 vec_sd_scenaru <- c(250, 500, 1000) # Abstract groups A and B use the _brand categorical palette (teal and # terracotta), not the semantic model/focus colours. barvy_skupin_ab <- c("A" = "#4E9A8A", "B" = "#C2613A") vec_popisky_scenaru <- paste0("SD = ", format_cz(vec_sd_scenaru, 0), " g") data_variability_demo <- tidyr::crossing( sd_g = vec_sd_scenaru, skupina = c("A", "B"), index = seq_along(vec_pattern) ) |> dplyr::mutate( prumer_g = ifelse( skupina == "A", prumer_skupiny_a, prumer_skupiny_a + rozdil_myslenkovy_pokus ), hmotnost_g = prumer_g + sd_g * vec_pattern[index], scenar = factor( paste0("SD = ", format_cz(sd_g, 0), " g"), levels = vec_popisky_scenaru ) ) ``` Následující obrázek není další `r render_glossary_term("vyber", display = "výběr")` tučňáků. Je to záměrně zjednodušený myšlenkový pokus se dvěma vymyšlenými skupinami A a B po deseti jedincích. Ve všech třech panelech se průměry skupin liší o stejných `r format_cz(rozdil_myslenkovy_pokus, 0)` g; panely se liší jen tím, jak daleko jsou jedinci rozptýleni kolem průměru své skupiny. ```{r} #| label: zobrazit-myslenkovy-pokus-variabilita #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Myšlenkový pokus se dvěma vymyšlenými skupinami A a B: rozdíl mezi průměry (vodorovné úsečky) je ve všech třech panelech stejný, směrodatná odchylka uvnitř skupin zleva doprava roste." #| fig-alt: "Tři panely se dvěma skupinami bodů A a B. Průměry skupin se ve všech panelech liší stejně; vlevo jsou body těsně u průměrů a skupiny se nepřekrývají, vpravo jsou body široce rozptýlené a skupiny se výrazně překrývají." #| fig-width: 10.5 #| fig-height: 4.6 ggplot2::ggplot( data = data_variability_demo, mapping = ggplot2::aes( x = skupina, y = hmotnost_g, color = skupina ) ) + ggbeeswarm::geom_quasirandom( width = 0.12, size = 2.2, alpha = 0.75, show.legend = FALSE ) + ggplot2::stat_summary( fun = mean, geom = "crossbar", width = 0.55, linewidth = 0.8, show.legend = FALSE ) + ggplot2::facet_wrap( facets = ggplot2::vars(scenar), nrow = 1 ) + ggplot2::scale_color_manual(values = barvy_skupin_ab) + ggplot2::scale_y_continuous( limits = c(1800, 6200), breaks = seq(2000, 6000, 1000) ) + ggplot2::labs( x = NULL, y = "Ilustrační hmotnost (g)" ) ``` ::: {.box-question title="Kde je rozdíl nejpřesvědčivější?"} Rozdíl průměrů je ve všech panelech stejný. Ve kterém panelu bychom nejsnáze odlišili systematický rozdíl mezi skupinami od běžných rozdílů mezi jedinci? A ve kterém nejhůře? ::: ::: {.box-result collapse="true" title="Možná odpověď"} Nejsnáze v levém panelu: skupiny se téměř nepřekrývají a rozdíl průměrů je mnohem větší než typická odchylka jedince od průměru jeho skupiny. Nejhůře v pravém panelu: rozdíl průměrů je zde stejně velký jako směrodatná odchylka uvnitř skupin a skupiny se výrazně překrývají. Stejný rozdíl průměrů tak může být přesvědčivý i sotva rozeznatelný. ::: **Rozdíl průměrů tedy nečteme izolovaně. Posuzujeme jej vzhledem k variabilitě, která zůstává kolem skupinových průměrů.** Z tohoto porovnání vychází také název metody `r render_glossary_term("anova", display = "analýza rozptylu")`, ke které se v této lekci dostaneme: porovnává variabilitu zachycenou rozdílnými skupinovými průměry s residuální variabilitou uvnitř skupin. Rozdíl hmotností zatím pouze popisujeme. Abychom k němu mohli připojit nejistotu a test, **fitujeme `r render_glossary_term("linearni-model", display = "lineární model")`**: `r render_glossary_term("odezva", display = "odezvou")` bude hmotnost těla a `r render_glossary_term("prediktor", display = "prediktorem")` druh. Před vložením `r render_glossary_term("kategoricka-promenna", display = "kategoriální proměnné")` do známé formule si připomeneme, jak model pracoval s `r render_glossary_term("numericka-promenna", display = "numerickým prediktorem")`. ## Stejný `r render_glossary_term("linearni-model", display = "lineární model")`, jiný `r render_glossary_term("prediktor", display = "prediktor")` ### Připomenutí: přímka pro délku ploutve a hmotnost V lekci o vztazích dvou proměnných jsme u tučňáků kreslili bodový graf délky ploutve a tělesné hmotnosti; v lekci o lineárním modelu jsme podobný vztah popsali přímkou. Přímku si připomeneme právě na tučňácích: délka ploutve (sloupec `delka_kridla_mm`) bude prediktorem a hmotnost těla `r render_glossary_term("odezva", display = "odezvou")`. ```{r} #| label: pripravit-pripomenuti-ciselneho-modelu #| echo: false #| output: false #| message: false #| warning: false #| error: false data_flipper_model <- data_penguins |> dplyr::filter(!is.na(delka_kridla_mm)) mod_flipper <- lm( formula = hmotnost_tela_g ~ delka_kridla_mm, data = data_flipper_model ) coef_flipper <- coef(mod_flipper) intercept_flipper <- unname(coef_flipper[["(Intercept)"]]) slope_flipper <- unname(coef_flipper[["delka_kridla_mm"]]) data_flipper_model <- data_flipper_model |> dplyr::mutate( fitted_flipper = fitted(mod_flipper), residual_flipper = resid(mod_flipper) ) # One clearly visible residual in the middle of the flipper range. data_flipper_residual <- data_flipper_model |> dplyr::filter(dplyr::between(delka_kridla_mm, 195, 205)) |> dplyr::slice_max(abs(residual_flipper), n = 1, with_ties = FALSE) # Teaching input: the slope is shown for a 10 mm longer flipper. x_slope_start <- 190 x_slope_end <- x_slope_start + 10 y_slope_start <- intercept_flipper + slope_flipper * x_slope_start y_slope_end <- intercept_flipper + slope_flipper * x_slope_end x_flipper_min <- min(data_flipper_model$delka_kridla_mm) x_flipper_max <- max(data_flipper_model$delka_kridla_mm) ``` Každý bod je jeden tučňák. Délka ploutve určuje jeho vodorovnou polohu a fialová přímka ukazuje, jakou hmotnost model pro tuto délku `r render_glossary_term("odhadnuta-hodnota", display = "odhaduje")`. ```{r} #| label: zobrazit-ciselny-model-primka #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Pozorované délky ploutve a hmotnosti tučňáků s odhadnutou přímkou. Každá hodnota numerického prediktoru vybírá jedno místo na přímce." #| fig-alt: "Bodový graf: šedé body tučňáků stoupají zleva doprava od krátkých ploutví a nízké hmotnosti k dlouhým ploutvím a vysoké hmotnosti; fialová přímka prochází středem oblaku." #| fig-width: 9 #| fig-height: 3.6 ggplot2::ggplot( data = data_flipper_model, mapping = ggplot2::aes( x = delka_kridla_mm, y = hmotnost_tela_g ) ) + ggplot2::geom_point( alpha = 0.55, size = 1.7, color = biostat_cols[["grey_olive"]] ) + ggplot2::annotate( geom = "segment", x = x_flipper_min, xend = x_flipper_max, y = intercept_flipper + slope_flipper * x_flipper_min, yend = intercept_flipper + slope_flipper * x_flipper_max, color = biostat_cols[["indigo_velvet"]], linewidth = 1.2 ) + ggplot2::labs( x = "Délka ploutve (mm)", y = "Hmotnost těla (g)" ) ``` {{< pagebreak >}} ### Kde je `r render_glossary_term("intercept", display = "intercept")`? V lekci o lineárním modelu jsme intercept a sklon psali jako $a$ a $b$. V lekcích o nejistotě a testování hypotéz jsme pro stejné koeficienty používali řecké písmeno $\beta$ (beta): $\beta_0$ je intercept a $\beta_1$ sklon v širší populaci. Stříška označuje odhad z našich dat, takže $\hat{\beta}_0$ je totéž co dřívější $a$ a $\hat{\beta}_1$ totéž co $b$. Intercept $\hat{\beta}_0$ je odhadnutá hmotnost tučňáka s ploutví dlouhou 0 mm. Nula leží daleko mimo pozorované délky ploutví (`r x_flipper_min`–`r x_flipper_max` mm), proto plnou přímku prodloužíme čárkovaně až k nule. Vychází $\hat{\beta}_0 = `r format_cz_math(intercept_flipper, 0)`$ g, tedy záporná hmotnost. Intercept umíme v rovnici vypočítat, ale tady nemá biologický smysl: je to extrapolace daleko za naše data. ```{r} #| label: zobrazit-ciselny-model-intercept #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Intercept je výška přímky při hodnotě prediktoru 0. Čárkovaná část upozorňuje, že mezi nulou a nejkratší pozorovanou ploutví nemáme žádná data." #| fig-alt: "Stejný bodový graf s osou délky ploutve od nuly. Fialová přímka je čárkovaně prodloužena doleva až k nule, kde fialový bod označuje intercept hluboko pod nulovou hmotností." #| fig-width: 9 #| fig-height: 4.4 ggplot2::ggplot( data = data_flipper_model, mapping = ggplot2::aes( x = delka_kridla_mm, y = hmotnost_tela_g ) ) + ggplot2::geom_hline( yintercept = 0, color = biostat_cols[["graphite"]], linewidth = 0.4 ) + ggplot2::geom_point( alpha = 0.45, size = 1.5, color = biostat_cols[["grey_olive"]] ) + ggplot2::annotate( geom = "segment", x = 0, xend = x_flipper_min, y = intercept_flipper, yend = intercept_flipper + slope_flipper * x_flipper_min, color = biostat_cols[["indigo_velvet"]], linetype = 2, linewidth = 0.9 ) + ggplot2::annotate( geom = "segment", x = x_flipper_min, xend = x_flipper_max, y = intercept_flipper + slope_flipper * x_flipper_min, yend = intercept_flipper + slope_flipper * x_flipper_max, color = biostat_cols[["indigo_velvet"]], linewidth = 1.2 ) + ggplot2::annotate( geom = "point", x = 0, y = intercept_flipper, color = biostat_cols[["indigo_velvet"]], size = 3.4 ) + ggplot2::annotate( geom = "label", x = 8, y = intercept_flipper, label = paste0("β̂₀ = ", format_cz(intercept_flipper, 0), " g"), hjust = 0, color = biostat_cols[["graphite"]], fill = "white", linewidth = 0.2, size = 3.6 ) + ggplot2::scale_x_continuous( breaks = seq(0, 250, 50), limits = c(0, 240) ) + ggplot2::labs( x = "Délka ploutve (mm)", y = "Hmotnost těla (g)" ) ``` {{< pagebreak >}} ### Kde je `r render_glossary_term("sklon", display = "sklon")`? Sklon $\hat{\beta}_1$ je změna odhadnuté hmotnosti při prodloužení ploutve o 1 mm. Jeden milimetr by byl v grafu téměř neviditelný, proto zvýrazníme prodloužení o `r x_slope_end - x_slope_start` mm. Odpovídající změna výšky přímky je $10 \cdot \hat{\beta}_1 = `r format_cz_math(10 * slope_flipper, 0)`$ g. ```{r} #| label: zobrazit-ciselny-model-sklon #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Sklon převádí změnu prediktoru na změnu odhadnuté hodnoty. Kvůli čitelnosti je zobrazena změna při prodloužení ploutve o 10 mm." #| fig-alt: "Bodový graf s fialovou přímkou; tmavá vodorovná šipka vyznačuje prodloužení ploutve o deset milimetrů a fialová svislá šipka odpovídající nárůst odhadnuté hmotnosti." #| fig-width: 9 #| fig-height: 4.6 ggplot2::ggplot( data = data_flipper_model, mapping = ggplot2::aes( x = delka_kridla_mm, y = hmotnost_tela_g ) ) + ggplot2::geom_point( alpha = 0.35, size = 1.6, color = biostat_cols[["grey_olive"]] ) + ggplot2::annotate( geom = "segment", x = x_flipper_min, xend = x_flipper_max, y = intercept_flipper + slope_flipper * x_flipper_min, yend = intercept_flipper + slope_flipper * x_flipper_max, color = biostat_cols[["indigo_velvet"]], linewidth = 1.2 ) + ggplot2::annotate( geom = "segment", x = x_slope_start, xend = x_slope_end, y = y_slope_start, yend = y_slope_start, color = biostat_cols[["graphite"]], linewidth = 1, arrow = grid::arrow(ends = "both", length = grid::unit(0.10, "cm")) ) + ggplot2::annotate( geom = "segment", x = x_slope_end, xend = x_slope_end, y = y_slope_start, yend = y_slope_end, color = biostat_cols[["indigo_velvet"]], linewidth = 1.1, arrow = grid::arrow(ends = "both", length = grid::unit(0.10, "cm")) ) + ggplot2::annotate( geom = "label", x = mean(c(x_slope_start, x_slope_end)), y = y_slope_start - 260, label = "změna x = 10 mm", color = biostat_cols[["graphite"]], fill = "white", linewidth = 0.2, size = 3.5 ) + ggplot2::annotate( geom = "label", x = x_slope_end + 1, y = mean(c(y_slope_start, y_slope_end)), label = paste0("10 · β̂₁ = ", format_cz(10 * slope_flipper, 0), " g"), hjust = 0, color = biostat_cols[["graphite"]], fill = "white", linewidth = 0.2, size = 3.5 ) + ggplot2::labs( x = "Délka ploutve (mm)", y = "Hmotnost těla (g)" ) ``` {{< pagebreak >}} ### Kde je `r render_glossary_term("residuum", display = "residuum")`? Residuum je svislá vzdálenost naměřené hodnoty od hodnoty na přímce. U zvýrazněného tučňáka s ploutví dlouhou `r data_flipper_residual$delka_kridla_mm` mm porovnáme naměřenou hmotnost `r format_cz(data_flipper_residual$hmotnost_tela_g, 0)` g s odhadem modelu `r format_cz(data_flipper_residual$fitted_flipper, 0)` g. Residuum je jejich rozdíl, `r format_cz(data_flipper_residual$residual_flipper, 0)` g. ```{r} #| label: zobrazit-ciselny-model-residuum #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Residuum jednoho pozorování je svislá vzdálenost mezi naměřenou hodnotou (plný oranžový bod) a odhadnutou hodnotou na přímce (prázdný fialový bod)." #| fig-alt: "Bodový graf s fialovou přímkou; u jednoho tučňáka spojuje oranžová svislá úsečka jeho naměřenou hmotnost s odhadnutou hmotností na přímce." #| fig-width: 9 #| fig-height: 4.6 ggplot2::ggplot( data = data_flipper_model, mapping = ggplot2::aes( x = delka_kridla_mm, y = hmotnost_tela_g ) ) + ggplot2::geom_point( alpha = 0.28, size = 1.5, color = biostat_cols[["grey_olive"]] ) + ggplot2::annotate( geom = "segment", x = x_flipper_min, xend = x_flipper_max, y = intercept_flipper + slope_flipper * x_flipper_min, yend = intercept_flipper + slope_flipper * x_flipper_max, color = biostat_cols[["indigo_velvet"]], linewidth = 1.2 ) + ggplot2::geom_segment( data = data_flipper_residual, mapping = ggplot2::aes( x = delka_kridla_mm, y = hmotnost_tela_g, xend = delka_kridla_mm, yend = fitted_flipper ), color = biostat_cols[["orange"]], linewidth = 1.3, inherit.aes = FALSE ) + ggplot2::geom_point( data = data_flipper_residual, mapping = ggplot2::aes( x = delka_kridla_mm, y = hmotnost_tela_g ), color = biostat_cols[["orange"]], size = 3.2, inherit.aes = FALSE ) + ggplot2::geom_point( data = data_flipper_residual, mapping = ggplot2::aes( x = delka_kridla_mm, y = fitted_flipper ), shape = 21, fill = "white", color = biostat_cols[["indigo_velvet"]], stroke = 1.1, size = 3.2, inherit.aes = FALSE ) + ggplot2::annotate( geom = "label", x = data_flipper_residual$delka_kridla_mm + 1.5, y = mean(c( data_flipper_residual$hmotnost_tela_g, data_flipper_residual$fitted_flipper )), label = paste0( "e = ", format_cz(data_flipper_residual$residual_flipper, 0), " g" ), hjust = 0, color = biostat_cols[["graphite"]], fill = "white", linewidth = 0.2, size = 3.5 ) + ggplot2::labs( x = "Délka ploutve (mm)", y = "Hmotnost těla (g)" ) ``` {{< pagebreak >}} ### Stejné koeficienty ve výstupu R Funkce `lm()` vrátí stejné dva koeficienty, které jsme v grafech hledali odděleně. ```{r} #| label: fitovat-model-ploutve-viditelne #| echo: true #| eval: true # Odhadneme přímku pro hmotnost podle délky ploutve mod_ploutev <- lm( formula = hmotnost_tela_g ~ delka_kridla_mm, data = data_tucnaci ) # Zobrazíme odhadnutý intercept a sklon coef(object = mod_ploutev) ``` Řádek `(Intercept)` je **intercept**, tedy výška přímky při délce ploutve 0 mm. Řádek `delka_kridla_mm` je **sklon**, tedy změna odhadnuté hmotnosti při prodloužení ploutve o 1 mm. Slovy: $$ \text{odhadnutá hmotnost} = \text{intercept} + \text{sklon} \cdot \text{délka ploutve}. $$ Pro tučňáka s ploutví dlouhou `r x_slope_start` mm: $$ \text{odhadnutá hmotnost} = `r format_cz_math(intercept_flipper, 1)` + `r format_cz_math(slope_flipper, 2)` \cdot `r x_slope_start` = `r format_cz_math(intercept_flipper + slope_flipper * x_slope_start, 1)`\ \text{g}. $$ Obecně, kde $x_i$ je délka ploutve tučňáka $i$, $\hat{y}_i$ jeho odhadnutá hmotnost, $\hat{\beta}_0$ intercept a $\hat{\beta}_1$ sklon: $$ \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i. $$ Residuum dál měří svislou vzdálenost konkrétního `r render_glossary_term("pozorovani", display = "pozorování")` od této přímky. U numerického prediktoru jsme tedy naměřenou hodnotu rozdělili na hodnotu odhadnutou modelem a residuum: $$ \text{naměřená hodnota} = \text{odhadnutá hodnota} + \text{residuum}, $$ $$ y_i = \hat{y}_i + e_i = \hat{\beta}_0 + \hat{\beta}_1 x_i + e_i. $$ Číslo $x_i$ určovalo polohu pozorování na vodorovné ose a výraz $\hat{\beta}_0 + \hat{\beta}_1 x_i$ vybral odpovídající výšku na přímce. Prediktor druh ale obsahuje názvy, ne čísla. Aby mohl vstoupit do stejného výpočtu, model si z něj vytvoří číselnou proměnnou. ## `r render_glossary_term("linearni-model", display = "Lineární model")` s `r render_glossary_term("kategoricka-promenna", display = "kategoriální proměnnou")` Vraťme se k hlavní biologické otázce: **liší se `r render_glossary_term("prumer", display = "průměrná hmotnost")` Adélie a Gentoo?** `r render_glossary_term("odezva", display = "Odezvou")` je stále hmotnost těla v gramech. Místo numerické délky ploutve ale modelu předáváme název druhu. Lineární model potřebuje způsob, jak dvě kategorie vložit do rovnice. Pro dva druhy stačí `r render_glossary_term("binarni-promenna", display = "binární proměnná")` **„je Gentoo?“**: hodnota 0 znamená „ne“, tedy Adélie, a hodnota 1 znamená „ano“, tedy Gentoo. Protože takovou proměnnou zapínáme a vypínáme, budeme jí říkat **přepínač**. Nejde o nové biologické měření ani pořadí druhů. Je to pouze číselný zápis příslušnosti ke skupině. Tabulka jej ukazuje jako nový sloupec vedle dvou skutečných `r render_glossary_term("pozorovani", display = "pozorování")`. ```{r} #| label: pripravit-prepinac-druhu #| echo: false #| output: false #| message: false #| warning: false #| error: false data_switch_examples <- data_penguins_two |> dplyr::group_by(druh) |> dplyr::slice(1) |> dplyr::ungroup() |> dplyr::mutate( je_gentoo = ifelse(druh == "Gentoo", 1L, 0L) ) mod_two_species <- lm( formula = hmotnost_tela_g ~ druh, data = data_penguins_two ) coef_two <- coef(mod_two_species) summary_two <- coef(summary(mod_two_species)) confint_two <- confint(mod_two_species) intercept_two <- unname(coef_two[["(Intercept)"]]) difference_two <- unname(coef_two[["druhGentoo"]]) se_difference_two <- summary_two["druhGentoo", "Std. Error"] t_difference_two <- summary_two["druhGentoo", "t value"] p_difference_two <- summary_two["druhGentoo", "Pr(>|t|)"] df_two <- df.residual(mod_two_species) data_binary_plot <- data_penguins_two |> dplyr::mutate( je_gentoo = ifelse(druh == "Gentoo", 1L, 0L), fitted_g = fitted(mod_two_species) ) |> dplyr::group_by(druh) |> dplyr::arrange(hmotnost_tela_g, .by_group = TRUE) |> dplyr::mutate( x_plot = je_gentoo + ( ( sin(seq_len(dplyr::n()) * 12.9898) * 43758.5453 ) %% 1 - 0.5 ) * 0.24 ) |> dplyr::ungroup() ``` ```{r} #| label: ukazat-prepinac-druhu #| echo: false #| output: true #| message: false #| warning: false #| error: false data_switch_examples |> dplyr::transmute( Druh = popisky_druhu[as.character(druh)], `Naměřená hmotnost (g)` = hmotnost_tela_g, `je Gentoo?` = je_gentoo ) |> knitr::kable( caption = paste0( "Dvě skutečná pozorování a číselný přepínač ", "pro jejich skupinu." ) ) ``` {{< pagebreak >}} V následujícím grafu tvoří stejná proměnná vodorovnou osu. Model stále dostává přesně hodnotu 0 nebo 1. Jednotlivé body v každé skupině pouze nepatrně posuneme doleva a doprava, aby se navzájem nepřekrývaly; tento zobrazovací posun nevytváří nové hodnoty `r render_glossary_term("prediktor", display = "prediktoru")`. Vodorovné úsečky jsou dvě `r render_glossary_term("odhadnuta-hodnota", display = "odhadnuté hodnoty")` a tečkované svislice jsou vzdálenosti jednotlivých měření od těchto odhadů. ```{r} #| label: zobrazit-binarni-prediktor #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-alt: "Body hmotností nad hodnotami 0 (Adélie) a 1 (Gentoo) s vodorovnou úsečkou průměru každého druhu a tečkovanými svislicemi od bodů k průměrům; šipka vpravo měří rozdíl průměrů." #| fig-cap: "Kategorie druhu zapsaná binární proměnnou 0/1. Model přiřadí všem jedincům s hodnotou 0 průměr Adélie a všem jedincům s hodnotou 1 průměr Gentoo. Body jsou jen kvůli překryvu mírně posunuty do stran; tečkované úsečky jsou residua." #| fig-width: 9 #| fig-height: 5.8 ggplot2::ggplot( data = data_binary_plot, mapping = ggplot2::aes( x = x_plot, y = hmotnost_tela_g, color = druh ) ) + ggplot2::geom_segment( mapping = ggplot2::aes( xend = x_plot, yend = fitted_g ), color = biostat_cols[["grey_olive"]], linetype = 3, alpha = 0.45, linewidth = 0.45 ) + ggplot2::geom_point(alpha = 0.55, size = 1.6) + ggplot2::annotate( geom = "segment", x = -0.18, xend = 0.18, y = intercept_two, yend = intercept_two, color = barvy_druhu[["Adelie"]], linewidth = 2.1 ) + ggplot2::annotate( geom = "segment", x = 0.82, xend = 1.18, y = intercept_two + difference_two, yend = intercept_two + difference_two, color = barvy_druhu[["Gentoo"]], linewidth = 2.1 ) + ggplot2::annotate( geom = "segment", x = 1.27, xend = 1.27, y = intercept_two, yend = intercept_two + difference_two, color = biostat_cols[["indigo_velvet"]], linewidth = 1, arrow = grid::arrow(ends = "both", length = grid::unit(0.11, "cm")) ) + ggplot2::annotate( geom = "label", x = 1.30, y = intercept_two + difference_two / 2, label = paste0("β̂₁ = ", format_cz(difference_two, 0), " g"), hjust = 0, color = biostat_cols[["graphite"]], fill = "white", linewidth = 0.2, size = 4 ) + ggplot2::scale_color_manual(values = barvy_druhu) + ggplot2::scale_x_continuous( breaks = c(0, 1), labels = c("0\nAdélie: ne", "1\nGentoo: ano"), limits = c(-0.25, 1.62) ) + ggplot2::labs( x = "Binární proměnná „je Gentoo?“", y = "Hmotnost těla (g)", color = "Druh" ) + ggplot2::theme(legend.position = "none") ``` Nula a jednička nejsou další biologická měření a neříkají, že jeden druh je „vyšší“ nebo „lepší“. Jsou jen jazykem, kterým modelu sdělíme příslušnost ke skupině. Stejně jako u numerického prediktoru zůstává odhadnutá hodnota ve tvaru `r render_glossary_term("intercept", display = "intercept")` plus koeficient krát hodnota prediktoru. ### Rovnice se dvěma hodnotami přepínače **Intercept** je v tomto modelu průměrná hmotnost referenčního druhu Adélie. **Druhý koeficient** je rozdíl mezi průměrem Gentoo a průměrem Adélie. U numerického prediktoru byl druhý koeficient `r render_glossary_term("sklon", display = "sklon")`; zde má nový význam skupinového rozdílu. Slovy: $$ \begin{aligned} \text{odhadnutá hmotnost těla} &= \text{průměrná hmotnost druhu Adélie} \\ &\quad + \text{rozdíl Gentoo vůči Adélie} \cdot \text{„je Gentoo?“}. \end{aligned} $$ Pro Adélie má přepínač hodnotu 0: $$ \begin{aligned} \text{odhadnutá hmotnost Adélie} &= `r format_cz_math(intercept_two, 1)` + `r format_cz_math(difference_two, 1)` \cdot 0 \\ &= `r format_cz_math(intercept_two, 1)`\ \text{g}. \end{aligned} $$ Pro Gentoo má přepínač hodnotu 1: $$ \begin{aligned} \text{odhadnutá hmotnost Gentoo} &= `r format_cz_math(intercept_two, 1)` + `r format_cz_math(difference_two, 1)` \cdot 1 \\ &= `r format_cz_math(intercept_two + difference_two, 1)`\ \text{g}. \end{aligned} $$ Obecně zapíšeme stejný vztah se symboly, které už známe. Symbol $x_i$ je nyní přepínač „je Gentoo?“ pro tučňáka $i$ (0 nebo 1), $\hat{\beta}_0$ odhadnutý průměr Adélie a $\hat{\beta}_1$ odhadnutý rozdíl Gentoo minus Adélie: $$ \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_i, $$ $$ y_i = \hat{\beta}_0 + \hat{\beta}_1 x_i + e_i. $$ Zápis se od numerického prediktoru nezměnil. Změnil se význam $x_i$ a tím i význam koeficientu $\hat{\beta}_1$: už nepopisuje změnu při zvýšení prediktoru o jednu měrnou jednotku, ale rozdíl jedné skupiny vůči referenční skupině. ::: {.box-warning title="Nula a jednička nejsou pořadí druhů"} Kdybychom vyměnili referenční druh, změnil by se přepínač, intercept i znaménko koeficientu. Skupinové průměry a biologický obsah dat by však zůstaly stejné. Kódování je jazyk modelu, ne vlastnost tučňáků. ::: ## Dvouskupinový model v R V R zůstává funkce `lm()` stejná. Do pravé strany formule vložíme `druh`. R ale musí rozhodnout, vůči kterému druhu bude zapisovat rozdíl. Proto zvolí jednu úroveň faktoru jako **referenční skupinu**. Ve výchozím kódování používá R jako referenci první úroveň faktoru. Při přípravě dat jsme pořadí nastavili na `Adelie`, `Chinstrap`, `Gentoo`; po ponechání dvou druhů proto zůstane první Adélie. Volba reference je nutná, protože informace „je Adélie?“ a „je Gentoo?“ jsou u dvou druhů navzájem dokonale propojené: známe-li jednu odpověď, druhou už dopočítáme. Model proto použije jeden společný základ – `r render_glossary_term("intercept", display = "intercept")` reference – a jeden rozdíl druhé skupiny vůči němu. ```{r} #| label: fitovat-model-dvou-druhu-viditelne #| echo: true #| eval: true # Ponecháme dva druhy a odstraníme nevyužitou třetí úroveň data_dva_druhy <- data_tucnaci[ data_tucnaci$druh %in% c("Adelie", "Gentoo"), ] data_dva_druhy$druh <- droplevels(data_dva_druhy$druh) # Zkontrolujeme pořadí úrovní; první je reference levels(data_dva_druhy$druh) # Fitujeme známý lineární model s kategoriálním prediktorem mod_dva_druhy <- lm( formula = hmotnost_tela_g ~ druh, data = data_dva_druhy ) # Prohlédneme si celý známý výstup modelu summary(mod_dva_druhy) ``` {{< pagebreak >}} ```{r} #| label: vysvetlit-koeficienty-dvou-druhu #| echo: false #| output: true #| message: false #| warning: false #| error: false tibble::tibble( `Řádek výstupu` = c("(Intercept)", "druhGentoo"), `Odhad (g)` = round( c(intercept_two, difference_two), digits = 1 ), `Biologický význam` = c( "průměrná hmotnost Adélie", "průměr Gentoo minus průměr Adélie" ), `Kde číslo vidíme v grafu` = c( "vzdálenost průměrové úsečky Adélie od nuly", "svislá vzdálenost mezi dvěma průměry" ) ) |> knitr::kable( caption = "Dvouskupinové koeficienty a jejich význam." ) ``` Tabulka překládá řádky výstupu do slov. V následujícím grafu jsou stejná dvě čísla umístěna přímo k datům: intercept je výška `r render_glossary_term("prumer", display = "průměru")` Adélie od nuly a koeficient `druhGentoo` je vzdálenost od průměru Adélie k průměru Gentoo. ```{r} #| label: zobrazit-koeficienty-dvou-druhu-v-grafu #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Dva koeficienty dvouskupinového modelu přímo v datech. Intercept β̂₀ je vzdálenost průměrové úsečky Adélie od nuly; koeficient β̂₁ je rozdíl Gentoo minus Adélie." #| fig-alt: "Body obou druhů na ose hmotnosti začínající nulou; šedá šipka od nuly k průměru Adélie označuje intercept a fialová šipka mezi průměry Adélie a Gentoo označuje rozdílový koeficient." #| fig-width: 9 #| fig-height: 5.6 ggplot2::ggplot( data = data_binary_plot, mapping = ggplot2::aes( x = x_plot, y = hmotnost_tela_g, color = druh ) ) + ggplot2::geom_point(alpha = 0.28, size = 1.5) + ggplot2::annotate( geom = "segment", x = -0.12, xend = 0.12, y = intercept_two, yend = intercept_two, color = barvy_druhu[["Adelie"]], linewidth = 1.8 ) + ggplot2::annotate( geom = "segment", x = 0.88, xend = 1.12, y = intercept_two + difference_two, yend = intercept_two + difference_two, color = barvy_druhu[["Gentoo"]], linewidth = 1.8 ) + ggplot2::annotate( geom = "segment", x = -0.23, xend = -0.23, y = 0, yend = intercept_two, color = biostat_cols[["graphite"]], linewidth = 1, arrow = grid::arrow(ends = "both", length = grid::unit(0.10, "cm")) ) + ggplot2::annotate( geom = "label", x = -0.20, y = intercept_two / 2, label = paste0( "β̂₀ = ", format( round(intercept_two, 1), big.mark = " ", decimal.mark = ",", nsmall = 1 ), " g" ), hjust = 0, color = biostat_cols[["graphite"]], fill = "white", linewidth = 0.2, size = 3.6 ) + ggplot2::annotate( geom = "segment", x = 0.52, xend = 0.52, y = intercept_two, yend = intercept_two + difference_two, color = biostat_cols[["indigo_velvet"]], linewidth = 1.1, arrow = grid::arrow(ends = "both", length = grid::unit(0.10, "cm")) ) + ggplot2::annotate( geom = "label", x = 0.55, y = intercept_two + difference_two / 2, label = paste0( "β̂₁ = ", format( round(difference_two, 1), big.mark = " ", decimal.mark = ",", nsmall = 1 ), " g" ), hjust = 0, color = biostat_cols[["indigo_velvet"]], fill = "white", linewidth = 0.2, size = 3.6 ) + ggplot2::annotate( geom = "label", x = 1, y = max( data_binary_plot$hmotnost_tela_g[data_binary_plot$je_gentoo == 1] ) + 180, label = paste0( "β̂₀ + β̂₁ = ", format( round(intercept_two + difference_two, 1), big.mark = " ", decimal.mark = ",", nsmall = 1 ), " g" ), color = biostat_cols[["graphite"]], fill = "white", linewidth = 0.2, size = 3.6 ) + ggplot2::scale_color_manual(values = barvy_druhu) + ggplot2::scale_x_continuous( breaks = c(0, 1), labels = c("Adélie\n0", "Gentoo\n1"), limits = c(-0.34, 1.28) ) + ggplot2::scale_y_continuous( limits = c(0, max(data_binary_plot$hmotnost_tela_g) + 400), expand = ggplot2::expansion(mult = c(0, 0.02)) ) + ggplot2::labs( x = "Binární proměnná „je Gentoo?“", y = "Hmotnost těla (g)", color = "Druh" ) + ggplot2::theme(legend.position = "none") ``` Intercept `r format_cz(intercept_two, 1)` g je průměr referenční skupiny Adélie. Koeficient `druhGentoo` `r format_cz(difference_two, 1)` g je rozdíl **Gentoo minus Adélie**. Součet obou koeficientů dává `r format_cz(intercept_two + difference_two, 1)` g, tedy průměr Gentoo. ::: {.box-data title="Doplňující: co se změní, když zvolíme jinou referenci?" collapse="true"} Referenci můžeme změnit funkcí `relevel()`. Když jako základ zvolíme Gentoo, intercept se stane průměrem Gentoo a rozdílový koeficient bude Adélie minus Gentoo. `r render_glossary_term("odhadnuta-hodnota", display = "Odhadnuté skupinové průměry")`, `r render_glossary_term("residuum", display = "residua")` i biologický obsah modelu zůstanou stejné; změní se pouze způsob jejich zápisu. ```{r} #| label: zmenit-referencni-druh #| echo: true #| eval: true data_reference_gentoo <- data_dva_druhy data_reference_gentoo$druh <- relevel( x = data_reference_gentoo$druh, ref = "Gentoo" ) levels(data_reference_gentoo$druh) mod_reference_gentoo <- lm( formula = hmotnost_tela_g ~ druh, data = data_reference_gentoo ) coef(mod_reference_gentoo) ``` Při referenci Gentoo je intercept `r format_cz(intercept_two + difference_two, 1)` g a koeficient `druhAdelie` je `r format_cz(-difference_two, 1)` g. Součet opět vrátí průměr Adélie `r format_cz(intercept_two, 1)` g. ::: ::: {.box-rule title="Ve studii s kontrolou nastavte referenci záměrně"} Pokud má studie skutečnou kontrolní skupinu, je obvykle nejčitelnější nastavit ji jako referenci. Intercept potom vyjadřuje průměr kontroly a každý léčebný koeficient odpovídá přímo otázce **ošetření minus kontrola**. ```r data_studie$skupina <- relevel( x = data_studie$skupina, ref = "kontrola" ) levels(data_studie$skupina) ``` Tato volba pouze zpřehlední koeficienty. Nevytváří randomizaci, kontrolní podmínku ani kauzální důkaz — ty závisejí na designu studie. Změna reference také nemění odhadnuté skupinové průměry, residua ani celkový test modelu. ::: ::: {.box-warning title="Co netestuje řádek interceptu"} `r render_glossary_term("p-hodnota", display = "P-hodnota")` v řádku interceptu porovnává průměr Adélie s nulovou hmotností. Nulová hmotnost není naše biologická otázka. Otázku rozdílu druhů ve dvouskupinovém modelu najdeme v řádku `druhGentoo`. ::: ## `r render_glossary_term("odhadnuta-hodnota", display = "Odhadnuté hodnoty")` a `r render_glossary_term("residuum", display = "residua")` nezmizely V modelu s číselným `r render_glossary_term("prediktor", display = "prediktorem")` měla každá délka ploutve vlastní hodnotu na přímce. V modelu se dvěma skupinami existují jen dvě odhadnuté hodnoty: `r render_glossary_term("prumer", display = "průměr")` Adélie a průměr Gentoo. Každý jedinec stejného druhu dostane stejnou odhadnutou hmotnost, ale vlastní residuum. ```{r} #| label: spocitat-odhadnute-hodnoty-residua #| echo: false #| output: false #| message: false #| warning: false #| error: false data_two_model <- data_binary_plot |> dplyr::transmute( druh, hmotnost_tela_g, odhadnuta_hmotnost_g = fitted_g, residuum_g = hmotnost_tela_g - fitted_g, x_plot = x_plot + 1 ) data_one_penguin <- data_two_model |> dplyr::filter(druh == "Adelie") |> dplyr::slice_max(abs(residuum_g), n = 1, with_ties = FALSE) one_measured <- data_one_penguin |> dplyr::pull(hmotnost_tela_g) one_fitted <- data_one_penguin |> dplyr::pull(odhadnuta_hmotnost_g) one_residual <- data_one_penguin |> dplyr::pull(residuum_g) ``` ```{r} #| label: zobrazit-odhadnute-hodnoty-residua #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Každý jedinec sdílí odhadnutou hodnotu svého druhu. Body jsou kvůli překryvu mírně posunuty do stran. Tenké svislé úsečky spojují skutečné body s odhadem a jsou tedy residui; jedno residuum je zvýrazněno." #| fig-alt: "Body obou druhů s průměrovou úsečkou každého druhu; šedé tečkované svislice spojují body s průměrem a jedna tmavá svislice u tučňáka kroužkového zvýrazňuje velké residuum." #| fig-width: 8 #| fig-height: 5.5 ggplot2::ggplot( data = data_two_model, mapping = ggplot2::aes( x = x_plot, y = hmotnost_tela_g, color = druh ) ) + ggplot2::geom_segment( mapping = ggplot2::aes( xend = x_plot, yend = odhadnuta_hmotnost_g ), color = biostat_cols[["grey_olive"]], alpha = 0.42, linewidth = 0.45, linetype = 3 ) + ggplot2::geom_point( alpha = 0.45, size = 1.4 ) + ggplot2::geom_segment( data = data_one_penguin, mapping = ggplot2::aes( x = x_plot, xend = x_plot, y = hmotnost_tela_g, yend = odhadnuta_hmotnost_g ), color = biostat_cols[["graphite"]], linewidth = 1.2, inherit.aes = FALSE ) + ggplot2::geom_point( data = data_one_penguin, mapping = ggplot2::aes( x = x_plot, y = hmotnost_tela_g ), color = biostat_cols[["graphite"]], size = 3, inherit.aes = FALSE ) + ggplot2::geom_segment( data = data_summary_two, mapping = ggplot2::aes( x = x - 0.28, xend = x + 0.28, y = prumer_g, yend = prumer_g, color = druh ), linewidth = 2.2, inherit.aes = FALSE ) + ggplot2::scale_color_manual(values = barvy_druhu) + ggplot2::scale_x_continuous( breaks = c(1, 2), labels = popisky_druhu[c("Adelie", "Gentoo")] ) + ggplot2::labs( x = NULL, y = "Hmotnost těla (g)", color = "Druh" ) + ggplot2::theme(legend.position = "none") ``` Pro zvýrazněného tučňáka kroužkového (jeho residuum je ze všech Adélie největší) platí slovy: $$ \text{residuum} = \text{naměřená hmotnost} - \text{odhadnutá hmotnost}. $$ Po dosazení čísel ze stejného řádku dat a z výšky průměrové úsečky: $$ \begin{aligned} \text{residuum} &= `r format_cz_math(one_measured, 1)` - `r format_cz_math(one_fitted, 1)` \\ &= `r format_cz_math(one_residual, 1)`\ \text{g}. \end{aligned} $$ Obecný zápis používá stejné symboly jako u přímky: $$ e_i = y_i - \hat{y}_i. $$ {{< pagebreak >}} ```{r} #| label: ziskat-odhadnute-hodnoty-residua-viditelne #| echo: true #| eval: true # Z modelu získáme odhadnutou hmotnost každého jedince data_dva_druhy$odhadnuta_hmotnost_g <- fitted(mod_dva_druhy) # Residuum je měření minus odhadnutá hodnota data_dva_druhy$residuum_g <- resid(mod_dva_druhy) # Zkontrolujeme první řádky obou nových sloupců head( data_dva_druhy[ , c( "druh", "hmotnost_tela_g", "odhadnuta_hmotnost_g", "residuum_g" ) ] ) ``` ::: {.box-result title="Co zůstalo stejné"} `r render_glossary_term("linearni-model", display = "Lineární model")` stále rozděluje každou naměřenou hodnotu na systematickou část a residuum. U číselného prediktoru byla systematickou částí přímka; u `r render_glossary_term("kategoricka-promenna", display = "kategoriálního prediktoru")` jsou jí skupinové průměry. ::: {{< pagebreak >}} ## Jak velký a přesný je rozdíl dvou druhů? Z dlouhého výstupu `summary()` vytáhneme jen dva `r render_glossary_term("odhad", display = "odhadnuté koeficienty")`. Funkce `coef()` je užitečná právě ve chvíli, kdy víme, na které hodnoty se chceme zaměřit: ```{r} #| label: pripomenout-koeficienty-dvou-druhu #| echo: true #| eval: true coef(mod_dva_druhy) ``` Koeficient `druhGentoo` odhaduje rozdíl `r format_cz(difference_two, 1)` g. Podobně jako v předchozích dvou lekcích odhad doplníme `r render_glossary_term("standardni-chyba", display = "standardní chybou")` a `r render_glossary_term("interval-spolehlivosti", display = "95% intervalem spolehlivosti")`. Interval nevzniká jako další nezávislé číslo. Připomeňme si z lekce o nejistotě, jak se vypočítá ze stejného odhadu a standardní chyby. Slovy: $$ 95\%\ \text{interval} = \text{odhad rozdílu} \pm \text{kritická hodnota} \cdot \text{standardní chyba rozdílu}. $$ Po dosazení odhadu, kritické hodnoty (vysvětlené níže) a standardní chyby z řádku `druhGentoo`: $$ 95\%\ \text{interval} = `r format_cz_math(difference_two, 1)` \pm `r format_cz_math(qt(0.975, df = df_two), 3)` \cdot `r format_cz_math(se_difference_two, 2)` = `r format_cz_math(confint_two["druhGentoo", 1], 1)`\ \text{až}\ `r format_cz_math(confint_two["druhGentoo", 2], 1)`\ \text{g}. $$ Se symboly, kde $\hat{\beta}_1$ je odhadnutý rozdíl Gentoo minus Adélie a $\operatorname{SE}(\hat{\beta}_1)$ jeho standardní chyba: $$ 95\%\ \text{interval} = \hat{\beta}_1 \pm \underbrace{t_{0{,}975,\,df}}_{\text{kritická hodnota}} \cdot \operatorname{SE}(\hat{\beta}_1). $$ Pro 95% interval necháváme 2,5 % `r render_glossary_term("rozdeleni", display = "t-rozdělení")` v každém krajním chvostu. `r render_glossary_term("kriticka-hodnota", display = "Kritická hodnota")` proto závisí na zvolené `r render_glossary_term("hladina-spolehlivosti", display = "hladině 95 %")` a na `r df_two` residuálních `r render_glossary_term("stupne-volnosti", display = "stupních volnosti")`; zde je `r format_cz(qt(0.975, df = df_two), 3)`. {{< pagebreak >}} ```{r} #| label: vypocitat-interval-dvou-druhu-viditelne #| echo: true #| eval: true # Vezmeme odhad a standardní chybu z řádku druhGentoo tabulka_koeficientu <- coef(summary(mod_dva_druhy)) odhad_rozdilu <- tabulka_koeficientu["druhGentoo", "Estimate"] standardni_chyba <- tabulka_koeficientu["druhGentoo", "Std. Error"] # Kritická hodnota pro 95% interval a residuální stupně volnosti kriticka_hodnota <- qt( p = 0.975, df = df.residual(mod_dva_druhy) ) dolni_mez <- odhad_rozdilu - kriticka_hodnota * standardni_chyba horni_mez <- odhad_rozdilu + kriticka_hodnota * standardni_chyba c(dolni_mez = dolni_mez, horni_mez = horni_mez) # Stejný výpočet pohodlně provede confint() confint(mod_dva_druhy)["druhGentoo", ] ``` ```{r} #| label: pripravit-interval-dvou-druhu #| echo: false #| output: false #| message: false #| warning: false #| error: false data_difference_interval <- tibble::tibble( porovnani = "Gentoo minus Adélie", odhad = difference_two, dolni = confint_two["druhGentoo", 1], horni = confint_two["druhGentoo", 2] ) ``` ```{r} #| label: zobrazit-interval-dvou-druhu #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Odhad rozdílu průměrné hmotnosti Gentoo minus Adélie a jeho 95% interval spolehlivosti. Nula znamená stejné průměry." #| fig-alt: "Fialový bod odhadu rozdílu s oranžovou vodorovnou úsečkou intervalu spolehlivosti; celý interval leží daleko vpravo od čárkované nuly." #| fig-width: 8 #| fig-height: 3.2 ggplot2::ggplot( data = data_difference_interval, mapping = ggplot2::aes( x = odhad, y = porovnani ) ) + ggplot2::geom_vline( xintercept = 0, linetype = 2, color = biostat_cols[["grey_olive"]] ) + ggplot2::geom_errorbar( mapping = ggplot2::aes( xmin = dolni, xmax = horni ), orientation = "y", width = 0.15, linewidth = 1.1, color = biostat_cols[["orange"]] ) + ggplot2::geom_point( size = 3.5, color = biostat_cols[["indigo_velvet"]] ) + ggplot2::labs( x = "Rozdíl průměrné hmotnosti (g)", y = NULL ) ``` Interval sahá od `r format_cz(confint_two["druhGentoo", 1], 1)` do `r format_cz(confint_two["druhGentoo", 2], 1)` g. Celý leží nad nulou, ale jeho hlavním sdělením není pouze „nula uvnitř / nula vně“. Ukazuje také, jak velký rozdíl je s daty a modelem slučitelný. ### Je rozdíl slučitelný s nulou? `r render_glossary_term("nulova-hypoteza", display = "Nulová hypotéza")` a alternativa mluví o rozdílu `r render_glossary_term("prumer", display = "průměrných")` hmotností obou druhů v širší populaci, ze které změření tučňáci pocházejí: $$ H_0: \text{průměr Gentoo} - \text{průměr Adélie} = 0, $$ $$ H_A: \text{průměr Gentoo} - \text{průměr Adélie} \neq 0. $$ Se symboly jde o tvrzení o populačním koeficientu: $H_0: \beta_1 = 0$ a $H_A: \beta_1 \neq 0$. Pro výpočet si znovu zobrazíme celý výstup modelu. V řádku `druhGentoo` stojí vedle sebe všechny tři hodnoty, které potřebujeme: odhad rozdílu, jeho standardní chyba a vypočtená t-statistika. {{< pagebreak >}} ```{r} #| label: pripomenout-summary-pro-t-test #| echo: true #| eval: true summary(mod_dva_druhy) ``` V řádku `druhGentoo` tedy čteme `Estimate` = `r format_cz(difference_two, 2)`, `Std. Error` = `r format_cz(se_difference_two, 2)` a `t value` = `r format_cz(t_difference_two, 2)`. Následující rovnice pouze ukazuje, jak R ze dvou prvních sloupců získalo třetí. Z předchozí lekce známe vztah mezi odhadem, nulovou hodnotou a jeho `r render_glossary_term("nejistota", display = "nejistotou")`. Slovy: $$ \text{t-statistika} = \frac{ \text{odhadnutý rozdíl}-\text{nulový rozdíl} }{ \text{standardní chyba rozdílu} }. $$ Po dosazení čísel z řádku `druhGentoo`: $$ \begin{aligned} t &= \frac{ `r format_cz_math(difference_two, 1)`-0 }{ `r format_cz_math(se_difference_two, 1)` } \\ &= `r format_cz_math(t_difference_two, 2)`. \end{aligned} $$ Obecně: $$ t = \frac{\hat{\beta}_1-0}{\operatorname{SE}(\hat{\beta}_1)}. $$ Stejně jako v předchozí lekci můžeme výsledek umístit do t-rozdělení očekávaného při nulové hypotéze. Musíme ale rozlišit dvě různá čísla: - **kritické hranice** jsou zde $-`r format_cz_math(qt(0.975, df = df_two), 3)`$ a $+`r format_cz_math(qt(0.975, df = df_two), 3)`$. Oddělují prostředních 95 % nulového t-rozdělení od krajních 5 %: v každém chvostu zůstává 2,5 %. Stejnou kladnou kritickou hodnotu jsme použili už při výpočtu 95% intervalu spolehlivosti; - **pozorovaná t-statistika** je $+`r format_cz_math(t_difference_two, 2)`$. Tu jsme právě vypočítali z řádku `druhGentoo` v `summary()`. Kritické hranice tedy **nejsou** výsledkem našeho testu. Vymezují předem zvolené rozhodovací oblasti pro $\alpha=0{,}05$. Náš výsledek $t=+`r format_cz_math(t_difference_two, 2)`$ leží mnohem dále od nuly. Protože testujeme oboustrannou otázku, sledujeme také stejně vzdálenou zrcadlovou hranici $-`r format_cz_math(abs(t_difference_two), 2)`$. `r render_glossary_term("p-hodnota", display = "P-hodnota")` je součet ploch nulového t-rozdělení za těmito dvěma hranicemi, nikoli plocha za kritickými hodnotami. Následující graf proto používá dva panely: první zvětší kritické hranice a oblasti o celkové ploše $\alpha=0{,}05$; druhý ukáže skutečně pozorované $t$ a jeho zrcadlovou hranici na celé ose. ```{r} #| label: pripravit-t-rozdeleni-dvou-druhu #| echo: false #| output: false #| message: false #| warning: false #| error: false t_critical <- qt(0.975, df = df_two) data_t_distribution <- dplyr::bind_rows( tibble::tibble( panel = "kritické hranice a alfa", t = seq(-6, 6, length.out = 1601) ), tibble::tibble( panel = "pozorované t a zrcadlová hranice", t = seq(-25, 25, length.out = 2401) ) ) |> dplyr::mutate( density = dt(t, df = df_two), tail = dplyr::case_when( t <= -t_critical ~ "levý chvost", t >= t_critical ~ "pravý chvost", TRUE ~ NA_character_ ), panel = factor( panel, levels = c( "kritické hranice a alfa", "pozorované t a zrcadlová hranice" ) ) ) data_t_critical <- tidyr::expand_grid( panel = levels(data_t_distribution$panel), t = c(-t_critical, t_critical) ) |> dplyr::mutate( panel = factor(panel, levels = levels(data_t_distribution$panel)) ) data_t_observed <- tibble::tibble( panel = factor( "pozorované t a zrcadlová hranice", levels = levels(data_t_distribution$panel) ), t = t_difference_two ) data_t_mirror <- tibble::tibble( panel = factor( "pozorované t a zrcadlová hranice", levels = levels(data_t_distribution$panel) ), t = -abs(t_difference_two) ) data_t_p_regions <- tibble::tibble( panel = factor( rep("pozorované t a zrcadlová hranice", 2), levels = levels(data_t_distribution$panel) ), xmin = c(-25, abs(t_difference_two)), xmax = c(-abs(t_difference_two), 25) ) ``` ```{r} #| label: zobrazit-t-rozdeleni-dvou-druhu #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: !expr paste0("T-rozdělení při nulové hypotéze. Levý panel ukazuje kritické hranice a fialové chvosty o celkové ploše alfa = 0,05. V pravém panelu fialové čáry tyto hranice připomínají, zatímco oranžové pásy označují oblasti za pozorovaným t = ", format_cz(t_difference_two, 2), " a zrcadlovou hranicí −", format_cz(abs(t_difference_two), 2), ". P-hodnota je plocha pod t-křivkou uvnitř těchto oranžových oblastí.") #| fig-alt: "Dva panely t-rozdělení. Vlevo zvonovitá křivka s fialově vybarvenými krajními chvosty za kritickými hranicemi. Vpravo stejná křivka na široké ose; oranžová čára pozorovaného t a jeho zrcadlová hranice leží tak daleko od nuly, že plocha za nimi není vidět." #| fig-width: 10 #| fig-height: 4.6 ggplot2::ggplot( data = data_t_distribution, mapping = ggplot2::aes(x = t, y = density) ) + ggplot2::geom_area( data = data_t_distribution |> dplyr::filter( panel == "kritické hranice a alfa", !is.na(tail) ), mapping = ggplot2::aes(group = tail), fill = biostat_cols[["amethyst"]], alpha = 0.42 ) + ggplot2::geom_rect( data = data_t_p_regions, mapping = ggplot2::aes( xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf ), fill = biostat_cols[["orange"]], alpha = 0.12, inherit.aes = FALSE ) + ggplot2::geom_line( color = biostat_cols[["graphite"]], linewidth = 1 ) + ggplot2::geom_vline( data = data_t_critical, mapping = ggplot2::aes(xintercept = t), color = biostat_cols[["amethyst"]], linetype = 2, linewidth = 0.9, inherit.aes = FALSE ) + ggplot2::geom_vline( data = data_t_observed, mapping = ggplot2::aes(xintercept = t), color = biostat_cols[["orange"]], linewidth = 1.1, inherit.aes = FALSE ) + ggplot2::geom_vline( data = data_t_mirror, mapping = ggplot2::aes(xintercept = t), color = biostat_cols[["orange"]], linetype = 2, linewidth = 1.1, inherit.aes = FALSE ) + ggplot2::geom_label( data = tibble::tibble( panel = factor( "kritické hranice a alfa", levels = levels(data_t_distribution$panel) ), t = t_critical, density = 0.31, label = paste0("kritické hranice = ±", format_cz(t_critical, 3)) ), mapping = ggplot2::aes(x = t, y = density, label = label), hjust = -0.05, color = biostat_cols[["amethyst"]], fill = "white", linewidth = 0.2, size = 3.2, inherit.aes = FALSE ) + ggplot2::geom_label( data = tibble::tibble( panel = factor( "pozorované t a zrcadlová hranice", levels = levels(data_t_distribution$panel) ), t = t_difference_two, density = 0.12, label = paste0("pozorované t = +", format_cz(t_difference_two, 2)) ), mapping = ggplot2::aes(x = t, y = density, label = label), hjust = 1.05, color = biostat_cols[["orange"]], fill = "white", linewidth = 0.2, size = 3.2, inherit.aes = FALSE ) + ggplot2::geom_label( data = tibble::tibble( panel = factor( "pozorované t a zrcadlová hranice", levels = levels(data_t_distribution$panel) ), t = -abs(t_difference_two), density = 0.12, label = paste0( "zrcadlová hranice = -", format_cz(abs(t_difference_two), 2) ) ), mapping = ggplot2::aes(x = t, y = density, label = label), hjust = -0.05, color = biostat_cols[["orange"]], fill = "white", linewidth = 0.2, size = 3.2, inherit.aes = FALSE ) + ggplot2::facet_wrap( facets = ggplot2::vars(panel), nrow = 1, scales = "free_x" ) + ggplot2::labs( x = "t-statistika očekávaná při nulové hypotéze", y = "Hustota" ) + ggplot2::theme( legend.position = "none", strip.text = ggplot2::element_text(face = "bold") ) ``` Pozorovaný rozdíl leží `r format_cz(t_difference_two, 2)` standardních chyb nad nulou, takže skutečně pozorovaná t-statistika je $+`r format_cz_math(t_difference_two, 2)`$. Záporná hodnota $-`r format_cz_math(abs(t_difference_two), 2)`$ nebyla pozorována; je to zrcadlová hranice potřebná pro oboustrannou otázku. Při `r df_two` residuálních stupních volnosti je součet ploch za oběma hranicemi, tedy p-hodnota, `r format_p_cz(p_difference_two)`. Za předpokladů modelu jsou data velmi obtížně slučitelná s přesně nulovým rozdílem mezi těmito dvěma druhy. Slovo „za předpokladů“ je důležité: residua tohoto modelu kontrolujeme stejně jako u přímky. Jak to vypadá u kategoriálního prediktoru, ukážeme níže na modelu se třemi druhy; stejné grafy platí i pro dva druhy. ### Historická poznámka: `r render_glossary_term("t-test", display = "dvouvýběrový t-test")` v literatuře Statistické metody se historicky vyučovaly jako soubor samostatně pojmenovaných testů pro různé situace: jeden test pro dva průměry, jiný pro více průměrů, další pro vztah dvou numerických proměnných. Takové názvy byly praktické v době ručních výpočtů a tabulek, a proto zůstávají v učebnicích, článcích i nabídkách statistických programů dodnes. Mnoho těchto klasických postupů ale můžeme zapsat jako speciální případy jednoho regresního modelu. Dvouskupinový t-test se společným `r render_glossary_term("rozptyl", display = "rozptylem")` používá za stejných předpokladů tentýž odhad rozdílu, standardní chybu, t-statistiku a p-hodnotu jako koeficient `druhGentoo` v našem dvouskupinovém `r render_glossary_term("linearni-model", display = "lineárním modelu")`. Nemusíme se tedy učit nový způsob uvažování: biologickou otázku, odhad, nejistotu i test stále čteme z jednoho modelového rámce. Název **dvouvýběrový t-test** však budeme potkávat v další literatuře. Následující příklad proto ukazuje, jak známý modelový výsledek přeložit do historického názvu testu a zpět. Nejde o další nesouvisející analýzu. Aby znaménko odpovídalo našemu koeficientu, odečteme Adélie od Gentoo. ::: {.box-result title="Proč tento historický název uvádíme"} Když v knize uvidíte dvouvýběrový t-test se společným rozptylem, můžete si jej představit jako lineární model `y ~ skupina` se dvěma úrovněmi. Při čtení výsledku se znovu ptejte na rozdíl průměrů, jeho standardní chybu, t-statistiku a p-hodnotu. ::: {{< pagebreak >}} ```{r} #| label: porovnat-model-a-t-test #| echo: true #| eval: true # Připravíme hmotnosti obou skupin ve směru Gentoo minus Adélie vec_gentoo <- data_dva_druhy$hmotnost_tela_g[ data_dva_druhy$druh == "Gentoo" ] vec_adelie <- data_dva_druhy$hmotnost_tela_g[ data_dva_druhy$druh == "Adelie" ] # Test se společným rozptylem odpovídá dvouskupinovému lm() t.test( x = vec_gentoo, y = vec_adelie, var.equal = TRUE ) ``` ```{r} #| label: ukazat-shodu-modelu-t-testu #| echo: false #| output: true #| message: false #| warning: false #| error: false res_t_equal <- t.test( x = data_penguins_two |> dplyr::filter(druh == "Gentoo") |> dplyr::pull(hmotnost_tela_g), y = data_penguins_two |> dplyr::filter(druh == "Adelie") |> dplyr::pull(hmotnost_tela_g), var.equal = TRUE ) tibble::tibble( Postup = c( "koeficient v lm()", "t-test se společným rozptylem" ), `Rozdíl Gentoo − Adélie (g)` = round( c( difference_two, unname(res_t_equal$estimate[[1]] - res_t_equal$estimate[[2]]) ), digits = 1 ), `t-statistika` = round( c(t_difference_two, unname(res_t_equal$statistic)), digits = 2 ), df = round( c(df_two, unname(res_t_equal$parameter)), digits = 0 ) ) |> knitr::kable( caption = paste0( "Stejná dvouskupinová otázka zapsaná jako koeficient ", "lineárního modelu a jako t-test." ) ) ``` ::: {.box-warning title="Výchozí t.test() v R je Welchův test"} Bez argumentu `var.equal = TRUE` používá R Welchův t-test, který nevyžaduje společný rozptyl skupin a může mít jiné stupně volnosti i mírně jiný interval. Není proto pravda, že každý výstup `t.test()` musí být číselně totožný s obyčejným `lm()`. Totožnost v tabulce platí pro stejnou otázku a stejné předpoklady. ::: ## Tři druhy, stejná formule Do dat se vrací tučňák uzdičkový (Chinstrap). Do `lm()` nemusíme přidávat nový druh ručně: faktor `druh` už obsahuje všechny tři úrovně. ```{r} #| label: fitovat-model-tri-druhy-skryte #| echo: false #| output: false #| message: false #| warning: false #| error: false mod_three_species <- lm( formula = hmotnost_tela_g ~ druh, data = data_penguins ) coef_three <- coef(mod_three_species) summary_three <- coef(summary(mod_three_species)) confint_three <- confint(mod_three_species) data_summary_three <- data_penguins |> dplyr::group_by(druh) |> dplyr::summarise( n = dplyr::n(), prumer_g = mean(hmotnost_tela_g), sd_g = sd(hmotnost_tela_g), .groups = "drop" ) |> dplyr::mutate( x = c(1, 2, 3) ) mean_chinstrap <- data_summary_three |> dplyr::filter(druh == "Chinstrap") |> dplyr::pull(prumer_g) difference_chinstrap_adelie <- mean_chinstrap - mean_adelie ``` ```{r} #| label: zobrazit-tri-druhy-prumery #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Tři druhy s průměrem ± 1 SD. Svislé závorky ukazují oba rozdíly, které model přímo vyjadřuje vůči referenčnímu druhu Adélie." #| fig-alt: "Body hmotností tří druhů s průměry a úsečkami plus minus jedna SD; průměry Adélie a Chinstrap leží téměř stejně vysoko, průměr Gentoo výrazně výše. Dvě šipky vyznačují rozdíly vůči Adélie." #| fig-width: 9 #| fig-height: 5.5 ggplot2::ggplot( data = data_penguins, mapping = ggplot2::aes( x = druh, y = hmotnost_tela_g, color = druh ) ) + ggbeeswarm::geom_quasirandom( width = 0.22, alpha = 0.42, size = 1.55, show.legend = FALSE ) + ggplot2::geom_errorbar( data = data_summary_three, mapping = ggplot2::aes( x = druh, ymin = prumer_g - sd_g, ymax = prumer_g + sd_g, color = druh ), width = 0.14, linewidth = 0.9, inherit.aes = FALSE, show.legend = FALSE ) + ggplot2::geom_segment( data = data_summary_three, mapping = ggplot2::aes( x = x - 0.28, xend = x + 0.28, y = prumer_g, yend = prumer_g, color = druh ), linewidth = 2.2, inherit.aes = FALSE ) + ggplot2::annotate( geom = "segment", x = 1.55, xend = 1.55, y = mean_adelie, yend = mean_chinstrap, color = biostat_cols[["graphite"]], linewidth = 1.1, arrow = grid::arrow(ends = "both", length = grid::unit(0.07, "cm")) ) + ggplot2::annotate( geom = "label", x = 1.58, y = mean(c(mean_adelie, mean_chinstrap)) + 115, label = paste0( "Chinstrap − Adélie\n", format_cz(difference_chinstrap_adelie, 0), " g" ), hjust = 0, color = biostat_cols[["graphite"]], fill = "white", linewidth = 0.2, size = 3.5 ) + ggplot2::annotate( geom = "segment", x = 2.62, xend = 2.62, y = mean_adelie, yend = mean_gentoo, color = biostat_cols[["graphite"]], linewidth = 1.1, arrow = grid::arrow(ends = "both", length = grid::unit(0.10, "cm")) ) + ggplot2::annotate( geom = "label", x = 2.66, y = mean(c(mean_adelie, mean_gentoo)), label = paste0( "Gentoo − Adélie\n", format_cz(difference_gentoo_adelie, 0), " g" ), hjust = 0, color = biostat_cols[["graphite"]], fill = "white", linewidth = 0.2, size = 3.5 ) + ggplot2::scale_color_manual(values = barvy_druhu) + ggplot2::scale_x_discrete(labels = popisky_druhu) + ggplot2::labs( x = NULL, y = "Hmotnost těla (g)" ) + ggplot2::theme(legend.position = "none") ``` `r render_glossary_term("prumer", display = "Průměry")` Adélie a Chinstrap leží blízko sebe, zatímco průměr Gentoo je mnohem výše. To je přesně situace, ve které potřebujeme odlišit tři otázky: jak se jednotlivé skupiny liší od reference, zda jsou všechny tři průměry celkově stejné a které konkrétní dvojice chceme porovnat. {{< pagebreak >}} ```{r} #| label: fitovat-model-tri-druhy-viditelne #| echo: true #| eval: true # Fitujeme stejnou formuli na všech třech druzích mod_tri_druhy <- lm( formula = hmotnost_tela_g ~ druh, data = data_tucnaci ) # Koeficienty stále čteme vzhledem k referenčnímu druhu Adélie coef(mod_tri_druhy) confint(mod_tri_druhy) ``` ### Tři skupiny potřebují dva přepínače Adélie zůstává referencí. Model proto vytvoří přepínač pro Chinstrap a druhý pro Gentoo. ```{r} #| label: ukazat-dva-prepinace #| echo: false #| output: true #| message: false #| warning: false #| error: false tibble::tibble( Druh = popisky_druhu[c("Adelie", "Chinstrap", "Gentoo")], `je Chinstrap?` = c(0L, 1L, 0L), `je Gentoo?` = c(0L, 0L, 1L), `Odhadnutá hmotnost` = c( "β̂₀", "β̂₀ + β̂₁", "β̂₀ + β̂₂" ) ) |> knitr::kable( caption = paste0( "Kódování tří skupin při referenčním druhu Adélie. ", "Přepínače nejsou biologické veličiny." ) ) ``` Slovy: $$ \begin{aligned} \text{odhadnutá hmotnost těla} &= \text{průměr Adélie} \\ &\quad + \text{rozdíl Chinstrap vůči Adélie} \cdot \text{„je Chinstrap?“} \\ &\quad + \text{rozdíl Gentoo vůči Adélie} \cdot \text{„je Gentoo?“}. \end{aligned} $$ Po dosazení skutečných koeficientů a obou přepínačů z tabulky skončí každý řádek právě na průměru daného druhu: $$ \begin{aligned} \widehat{\text{hmotnost}}_{\text{Adélie}} &= `r format_cz_math(coef_three[["(Intercept)"]], 1)` + `r format_cz_math(coef_three[["druhChinstrap"]], 1)`\cdot 0 + `r format_cz_math(coef_three[["druhGentoo"]], 1)`\cdot 0 = `r format_cz_math(coef_three[["(Intercept)"]], 1)`\ \text{g}, \\ \widehat{\text{hmotnost}}_{\text{Chinstrap}} &= `r format_cz_math(coef_three[["(Intercept)"]], 1)` + `r format_cz_math(coef_three[["druhChinstrap"]], 1)`\cdot 1 + `r format_cz_math(coef_three[["druhGentoo"]], 1)`\cdot 0 = `r format_cz_math(coef_three[["(Intercept)"]] + coef_three[["druhChinstrap"]], 1)`\ \text{g}, \\ \widehat{\text{hmotnost}}_{\text{Gentoo}} &= `r format_cz_math(coef_three[["(Intercept)"]], 1)` + `r format_cz_math(coef_three[["druhChinstrap"]], 1)`\cdot 0 + `r format_cz_math(coef_three[["druhGentoo"]], 1)`\cdot 1 = `r format_cz_math(coef_three[["(Intercept)"]] + coef_three[["druhGentoo"]], 1)`\ \text{g}. \end{aligned} $$ Obecně: symbol $x_{\text{Chinstrap},i}$ je číselná odpověď na otázku „je `r render_glossary_term("pozorovani", display = "pozorování")` $i$ Chinstrap?“ a $x_{\text{Gentoo},i}$ na otázku „je pozorování $i$ Gentoo?“. Oba nabývají hodnoty 0 nebo 1. Koeficient $\hat{\beta}_0$ je odhadnutý průměr Adélie, $\hat{\beta}_1$ rozdíl Chinstrap minus Adélie a $\hat{\beta}_2$ rozdíl Gentoo minus Adélie: $$ \hat{y}_i = \hat{\beta}_0 + \hat{\beta}_1 x_{\text{Chinstrap},i} + \hat{\beta}_2 x_{\text{Gentoo},i}. $$ ```{r} #| label: vysvetlit-koeficienty-tri-druhy #| echo: false #| output: true #| message: false #| warning: false #| error: false data_coef_three <- tibble::tibble( radek = c("(Intercept)", "druhChinstrap", "druhGentoo"), odhad_g = unname(coef_three), vyznam = c( "průměr Adélie", "Chinstrap minus Adélie", "Gentoo minus Adélie" ) ) data_coef_three |> dplyr::transmute( `Řádek výstupu` = radek, `Odhad (g)` = round(odhad_g, digits = 1), `Význam při referenci Adélie` = vyznam ) |> knitr::kable( caption = "Koeficienty třískupinového modelu." ) ``` ::: {.box-warning title="Koeficienty neobsahují všechny dvojice"} Tabulka s referencí Adélie přímo odpovídá na otázky Chinstrap minus Adélie a Gentoo minus Adélie. Rozdíl Gentoo minus Chinstrap v ní není samostatný řádek, i když jej lze ze stejného modelu odvodit. ::: ## Popisuje model data rozumně? Než třískupinový model použijeme k testům, zkontrolujeme stejně jako u přímky `r render_glossary_term("diagnostika-modelu", display = "diagnostiku modelu")`. U `r render_glossary_term("kategoricka-promenna", display = "kategoriálního prediktoru")` má graf `r render_glossary_term("residuum", display = "residuí")` proti `r render_glossary_term("odhadnuta-hodnota", display = "odhadnutým hodnotám")` zvláštní, ale očekávaný tvar: model odhaduje jen tři skupinové `r render_glossary_term("prumer", display = "průměry")`, takže body tvoří svislé pásy, jeden pro každý druh. Průměry Adélie a Chinstrap se liší jen o `r format_cz(difference_chinstrap_adelie, 0)` g, proto jejich pásy na vodorovné ose téměř splývají a v grafu vidíme zdánlivě jen dva. ::: {.box-warning title="Co od diagnostiky potřebujeme"} Test předpokládá nezávislá `r render_glossary_term("pozorovani", display = "pozorování")`, přibližně srovnatelnou `r render_glossary_term("variabilita", display = "variabilitu")` residuí ve skupinách a přibližně normální tvar residuí. Graf residuí proti odhadnutým hodnotám pomáhá posoudit srovnatelnost `r render_glossary_term("rozptyl", display = "rozptylu")`; `r render_glossary_term("qq-graf", display = "Q–Q graf")` pomáhá posoudit tvar `r render_glossary_term("rozdeleni", display = "rozdělení")`. Nezávislost však z těchto grafů nepoznáme — musíme ji zdůvodnit způsobem, jakým byla pozorování získána. ::: ::: {.box-data collapse="true" title="Doplňující: co když se skupiny liší variabilitou?"} Rozdílný rozptyl residuí neznamená, že rozdíl průměrů zmizel. Znamená, že běžný model se společným residuálním rozptylem nemusí správně vyjádřit nejistotu tohoto rozdílu. Nápovědou může být výrazně širší pás residuí u jednoho druhu nebo skupinové rozdělení s velmi odlišnou šířkou. Pro dvě skupiny už jsme viděli Welchův t-test, který společný rozptyl nevyžaduje. Pro více skupin existuje Welchova ANOVA a v modelovém rámci také další způsoby výpočtu odolnější nejistoty. Volba řešení závisí na otázce, velikosti skupin a tvaru dat; diagnostický graf proto nemá sloužit jako automatický příkaz k transformaci nebo mazání pozorování. ::: ```{r} #| label: pripravit-diagnostiku-tri-druhy #| echo: false #| output: false #| message: false #| warning: false #| error: false data_diagnostics <- data_penguins |> dplyr::mutate( odhadnuta_hmotnost_g = fitted(mod_three_species), residuum_g = resid(mod_three_species) ) ``` ```{r} #| label: zobrazit-residua-proti-odhadum #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Residua proti odhadnutým hodnotám. Svislé pásy odpovídají skupinovým průměrům a jsou vlastností tohoto modelu, ne samy o sobě závadou. Pásy Adélie a Chinstrap téměř splývají, protože jejich průměry jsou skoro stejné." #| fig-alt: "Residua tučňáků proti odhadnuté hmotnosti: oranžové body Adélie a fialové body Chinstrap tvoří společný svislý pás vlevo, modrozelené body Gentoo pás vpravo; oba pásy jsou rozptýleny kolem nulové čáry podobně široce." #| fig-width: 8 #| fig-height: 5 ggplot2::ggplot( data = data_diagnostics, mapping = ggplot2::aes( x = odhadnuta_hmotnost_g, y = residuum_g, color = druh ) ) + ggplot2::geom_hline( yintercept = 0, linetype = 2, color = biostat_cols[["grey_olive"]] ) + ggplot2::geom_point( alpha = 0.58, size = 1.7, position = ggplot2::position_jitter(width = 22, height = 0) ) + ggplot2::scale_color_manual( values = barvy_druhu, labels = popisky_druhu ) + ggplot2::labs( x = "Odhadnutá hmotnost (g)", y = "Residuum (g)", color = "Druh" ) ``` Ve všech pásech jsou residua rozptýlena kolem nuly. Rozsah není dokonale totožný, ale žádná skupina nevytváří dramaticky odlišný vějíř a žádný jednotlivý bod zjevně neurčuje celý výsledek. Stejné dva grafy nakreslíte v základním R. Graf residuí proti odhadnutým hodnotám vznikne funkcí `plot()`, vodorovnou čáru v nule přidá `abline()`. Q–Q graf nakreslí `qqnorm()` a referenční přímku `qqline()`. Na rozdíl od grafů výše nebudou body obarvené podle druhu. ```{r} #| label: diagnostika-tri-druhy-viditelne #| echo: true #| eval: false # Residua proti hodnotám odhadnutým modelem plot( x = fitted(object = mod_tri_druhy), y = resid(object = mod_tri_druhy), xlab = "Odhadnutá hmotnost (g)", ylab = "Residuum (g)" ) # Přerušovaná čára v nule abline(h = 0, lty = 2) # Q–Q graf residuí s referenční přímkou qqnorm(y = resid(object = mod_tri_druhy)) qqline(y = resid(object = mod_tri_druhy)) ``` ```{r} #| label: zobrazit-qq-graf-residui #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Q–Q graf residuí. Body převážně sledují referenční přímku; menší odchylky na koncích jsou u reálných dat běžné." #| fig-alt: "Q–Q graf: šedé body residuí leží téměř na tmavé referenční přímce, jen na obou koncích se od ní mírně vzdalují." #| fig-width: 7 #| fig-height: 5 ggplot2::ggplot( data = data_diagnostics, mapping = ggplot2::aes( sample = residuum_g ) ) + ggplot2::stat_qq( color = biostat_cols[["grey_olive"]], alpha = 0.55, size = 1.6 ) + ggplot2::stat_qq_line( color = biostat_cols[["graphite"]], linewidth = 0.9 ) + ggplot2::labs( x = "Teoretické kvantily normálního rozdělení", y = "Kvantily residuí" ) ``` Q–Q graf porovnává seřazená residua s hodnotami očekávanými při přibližně normálním tvaru. Body většinou sledují přímku. Model proto pro srovnání skupin nevykazuje problém, který by převzal hlavní biologický příběh, a můžeme jej použít k testům. To neznamená, že jsou předpoklady „dokázané“; znamená to, že grafy neukazují zásadní rozpor. ## Jedna celková otázka: `r render_glossary_term("anova", display = "analýza rozptylu")` U dvou skupin stačil jediný koeficient: Gentoo minus Adélie. U tří skupin už model obsahuje dva rozdílové koeficienty a jejich samostatné řádky neodpovídají přímo na jednu společnou otázku. Celková otázka zní, zda přidání `r render_glossary_term("promenna", display = "proměnné")` `druh` zlepší popis dat jako celek: $$ H_0: \text{všechny tři průměrné hmotnosti jsou stejné}. $$ Alternativa neříká, která dvojice se liší. Říká pouze, že alespoň jeden `r render_glossary_term("prumer", display = "průměr")` není stejný jako ostatní. ### Jeden společný průměr proti třem skupinovým průměrům `r render_glossary_term("nulova-hypoteza", display = "Nulová hypotéza")` odpovídá jednoduššímu modelu bez `r render_glossary_term("prediktor", display = "prediktoru")`: všem tučňákům přiřadí jeden společný průměr. Model s druhem je složitější a přiřadí každému jedinci průměr jeho druhu. Otázku analýzy rozptylu proto můžeme přeformulovat: **zkrátí tři skupinové průměry `r render_glossary_term("residuum", display = "residua")` natolik, aby přidaná složitost modelu stála za to?** ```{r} #| label: pripravit-porovnani-jednoho-a-tri-prumeru #| echo: false #| output: false #| message: false #| warning: false #| error: false mod_one_mean <- lm( formula = hmotnost_tela_g ~ 1, data = data_penguins ) common_mean_g <- unname(coef(mod_one_mean)[["(Intercept)"]]) data_plot_positions <- data_penguins |> dplyr::group_by(druh) |> dplyr::arrange(hmotnost_tela_g, .by_group = TRUE) |> dplyr::mutate( x_plot = as.numeric(druh) + ( ( sin(seq_len(dplyr::n()) * 12.9898) * 43758.5453 ) %% 1 - 0.5 ) * 0.26 ) |> dplyr::ungroup() data_model_comparison <- dplyr::bind_rows( data_plot_positions |> dplyr::mutate( model = "jeden společný průměr", fitted_g = common_mean_g ), data_plot_positions |> dplyr::mutate( model = "tři skupinové průměry", fitted_g = predict( mod_three_species, newdata = data_plot_positions ) ) ) |> dplyr::mutate( residual_g = hmotnost_tela_g - fitted_g, model = factor( model, levels = c( "jeden společný průměr", "tři skupinové průměry" ) ) ) data_common_line <- tibble::tibble( model = factor( "jeden společný průměr", levels = levels(data_model_comparison$model) ), x = 0.68, xend = 3.32, fitted_g = common_mean_g ) data_group_lines <- data_summary_three |> dplyr::transmute( model = factor( "tři skupinové průměry", levels = levels(data_model_comparison$model) ), druh, x_start = x - 0.32, x_end = x + 0.32, fitted_g = prumer_g ) sse_one_mean <- sum(resid(mod_one_mean)^2) sse_three_means <- sum(resid(mod_three_species)^2) sse_reduction <- sse_one_mean - sse_three_means data_group_shifts <- data_summary_three |> dplyr::transmute( druh, x = x + 0.30, common_mean_g = common_mean_g, group_mean_g = prumer_g, model = factor( "jeden společný průměr", levels = levels(data_model_comparison$model) ) ) data_sse_labels <- tibble::tibble( model = factor( c("jeden společný průměr", "tři skupinové průměry"), levels = levels(data_model_comparison$model) ), x = 1.05, y = max(data_penguins$hmotnost_tela_g) + 160, label = paste0( "součet čtverců residuí = ", format( round(c(sse_one_mean, sse_three_means)), big.mark = " ", scientific = FALSE ), " g²" ) ) data_ss_parts <- tibble::tibble( cast = factor( c("zachyceno rozdíly druhů", "zůstává uvnitř druhů"), levels = c("zůstává uvnitř druhů", "zachyceno rozdíly druhů") ), ss = c(sse_reduction, sse_three_means), label = paste0( c("druh: ", "residua: "), format( round(c(sse_reduction, sse_three_means)), big.mark = " ", scientific = FALSE ), " g²" ) ) ``` ```{r} #| label: zobrazit-jeden-proti-trem-prumerum #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Stejná pozorování popsaná jedním společným průměrem a třemi skupinovými průměry. Body jsou kvůli překryvu mírně posunuty do stran. Tečkované svislice jsou všechna residua; šipky v levém panelu ukazují posun od společného průměru k průměrům druhů." #| fig-alt: "Dva panely se stejnými body tří druhů. Vlevo jedna společná vodorovná čára a dlouhé tečkované svislice residuí, vpravo tři skupinové průměry a kratší residua." #| fig-width: 10 #| fig-height: 5.8 ggplot2::ggplot( data = data_model_comparison, mapping = ggplot2::aes( x = x_plot, y = hmotnost_tela_g, color = druh ) ) + ggplot2::geom_segment( mapping = ggplot2::aes( xend = x_plot, yend = fitted_g ), color = biostat_cols[["grey_olive"]], linetype = 3, alpha = 0.38, linewidth = 0.38 ) + ggplot2::geom_point( alpha = 0.38, size = 1.25 ) + ggplot2::geom_segment( data = data_common_line, mapping = ggplot2::aes( x = x, xend = xend, y = fitted_g, yend = fitted_g ), color = biostat_cols[["graphite"]], linewidth = 1.8, inherit.aes = FALSE, show.legend = FALSE ) + ggplot2::geom_segment( data = data_group_lines, mapping = ggplot2::aes( x = x_start, xend = x_end, y = fitted_g, yend = fitted_g, color = druh ), linewidth = 1.8, inherit.aes = FALSE, show.legend = FALSE ) + ggplot2::geom_segment( data = data_group_shifts, mapping = ggplot2::aes( x = x, xend = x, y = common_mean_g, yend = group_mean_g, color = druh ), linewidth = 1.1, arrow = grid::arrow(ends = "both", length = grid::unit(0.09, "cm")), inherit.aes = FALSE, show.legend = FALSE ) + ggplot2::geom_label( data = data_sse_labels, mapping = ggplot2::aes( x = x, y = y, label = label ), hjust = 0, size = 3.3, color = biostat_cols[["graphite"]], fill = "white", linewidth = 0.2, inherit.aes = FALSE ) + ggplot2::facet_wrap( facets = ggplot2::vars(model), nrow = 1 ) + ggplot2::scale_color_manual(values = barvy_druhu) + ggplot2::scale_x_continuous( breaks = c(1, 2, 3), labels = c("Adélie", "Chinstrap", "Gentoo") ) + ggplot2::labs( x = NULL, y = "Hmotnost těla (g)", color = "Druh" ) + ggplot2::theme(legend.position = "none") ``` V levém panelu jsou residua vzdálenosti od jednoho společného průměru. V pravém panelu jsou vzdálenosti od průměru vlastního druhu. `r render_glossary_term("soucet-ctvercu-residui", display = "Součet čtverců residuí")` klesne z `r format_cz(sse_one_mean, 0)` na `r format_cz(sse_three_means, 0)` g². Rozdílné skupinové průměry tedy vysvětlí část `r render_glossary_term("variabilita", display = "variability")`, která v jednodušším modelu zůstávala v residuích. Proč residua umocňujeme? Kdybychom pouze sečetli kladné a záporné vzdálenosti od průměru, navzájem by se vyrušily. Čtverce jsou vždy nezáporné a velkým odchylkám dávají větší váhu. Stejnou myšlenku součtu čtverců residuí už známe z fitování přímky. Pokles součtu čtverců můžeme zobrazit jako rozdělení původní chyby jednoduššího modelu na dvě části: $$ \underbrace{\text{SSE jednoho průměru}}_{`r format_cz_math(sse_one_mean, 0)`} = \underbrace{\text{SS zachycený druhem}}_{`r format_cz_math(sse_reduction, 0)`} + \underbrace{\text{SSE tří průměrů}}_{`r format_cz_math(sse_three_means, 0)`}. $$ Obecně: $\mathrm{SSE}_{\text{jeden průměr}} = \mathrm{SS}_{\text{druh}} + \mathrm{SSE}_{\text{tři průměry}}$, kde SSE je součet čtverců residuí daného modelu a $\mathrm{SS}_{\text{druh}}$ část, kterou zachytily rozdíly průměrů druhů. ```{r} #| label: zobrazit-rozklad-souctu-ctvercu #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Původní součet čtverců residuí se po přidání druhu rozdělí na část zachycenou rozdíly skupinových průměrů a část, která stále zůstává uvnitř druhů." #| fig-alt: "Vodorovný sloupec součtu čtverců rozdělený na fialovou část zachycenou rozdíly druhů a šedou část, která zůstává uvnitř druhů; fialová část je větší." #| fig-width: 9 #| fig-height: 3.2 ggplot2::ggplot( data = data_ss_parts, mapping = ggplot2::aes( x = "jeden společný průměr", y = ss, fill = cast ) ) + ggplot2::geom_col(width = 0.55) + ggplot2::geom_label( mapping = ggplot2::aes(label = label), position = ggplot2::position_stack(vjust = 0.5), color = biostat_cols[["graphite"]], fill = "white", linewidth = 0.2, label.padding = grid::unit(0.15, "lines"), size = 3.5 ) + ggplot2::coord_flip() + ggplot2::scale_y_continuous( labels = function(x) paste(format_cz(x / 1e6, 0), "mil.") ) + ggplot2::scale_fill_manual( values = c( "zůstává uvnitř druhů" = biostat_cols[["grey_olive"]], "zachyceno rozdíly druhů" = biostat_cols[["amethyst"]] ) ) + ggplot2::labs( x = NULL, y = "Součet čtverců (g²)", fill = NULL ) + ggplot2::theme(legend.position = "bottom") ``` ### Co porovnává `r render_glossary_term("f-statistika", display = "F-statistika")` Samotný pokles součtu čtverců nestačí. Složitější model téměř vždy residua alespoň trochu zkrátí, protože dostal více koeficientů. Potřebujeme proto porovnat dvě veličiny na férovém měřítku: 1. kolik variability zachytil každý přidaný rozdíl mezi skupinami; 2. jak velký čtverec residua typicky zůstává uvnitř skupin. F-statistika je poměr právě těchto dvou množství. Protože porovnává dvě variability, jmenuje se metoda **analýza rozptylu**, i když biologická otázka mluví o průměrech. Slovy: $$ \text{F-statistika} = \frac{ \text{variabilita připadající na rozdíly mezi skupinami} }{ \text{residuální variabilita uvnitř skupin} }. $$ ```{r} #| label: pripravit-anova-vysledek #| echo: false #| output: false #| message: false #| warning: false #| error: false table_anova <- anova(mod_three_species) f_anova <- table_anova["druh", "F value"] p_anova <- table_anova["druh", "Pr(>F)"] df_group <- table_anova["druh", "Df"] df_residual <- table_anova["Residuals", "Df"] ms_group <- table_anova["druh", "Mean Sq"] ms_residual <- table_anova["Residuals", "Mean Sq"] data_f_distribution <- dplyr::bind_rows( tibble::tibble( panel = "nulové F-rozdělení pod lupou", f = seq(0, 8, length.out = 1601) ), tibble::tibble( panel = "pozorované F na celé ose", f = seq(0, ceiling(f_anova * 1.03), length.out = 2401) ) ) |> dplyr::mutate( density = df(f, df1 = df_group, df2 = df_residual), panel = factor( panel, levels = c( "nulové F-rozdělení pod lupou", "pozorované F na celé ose" ) ) ) data_f_observed <- tibble::tibble( panel = factor( "pozorované F na celé ose", levels = levels(data_f_distribution$panel) ), f = f_anova ) ``` Čísla už jsme viděli ve dvou konkrétních grafech. V grafu **jeden společný průměr proti třem skupinovým průměrům** byl součet čtverců residuí jednoho průměru `r format_cz(sse_one_mean, 0)` g² a součet čtverců residuí tří průměrů `r format_cz(sse_three_means, 0)` g². Jejich rozdíl `r format_cz(sse_reduction, 0)` g² tvoří fialovou část grafu rozkladu součtu čtverců – část zachycenou druhem. Součty čtverců převedeme na srovnatelné **průměrné čtverce** dělením příslušnými `r render_glossary_term("stupne-volnosti", display = "stupni volnosti")`. Model přidal `r df_group` rozdílové koeficienty; uvnitř skupin zbývá `r df_residual` residuálních stupňů volnosti: $$ \begin{aligned} \text{průměrný čtverec mezi skupinami} &= \frac{ `r format_cz_math(sse_one_mean, 0)` - `r format_cz_math(sse_three_means, 0)` }{ \text{`r df_group`} } = \frac{ `r format_cz_math(sse_reduction, 0)` }{ \text{`r df_group`} } = `r format_cz_math(ms_group, 0)`, \\ \text{průměrný čtverec uvnitř skupin} &= \frac{ `r format_cz_math(sse_three_means, 0)` }{ \text{`r df_residual`} } = `r format_cz_math(ms_residual, 0)`. \end{aligned} $$ {{< pagebreak >}} Kdyby druh nepomáhal, zlepšení připadající na jeden nový koeficient by bylo podobně velké jako běžný čtverec residua a $F$ by bylo kolem 1. Čím je čitatel větší vůči jmenovateli, tím více skupinové rozdíly vyčnívají nad variabilitu jedinců uvnitř skupin. V našich datech jejich poměr dává: $$ F = \frac{ `r format_cz_math(ms_group, 0)` }{ `r format_cz_math(ms_residual, 0)` } = `r format_cz_math(f_anova, 2)`. $$ Obecně: $$ F = \frac{MS_{\text{mezi skupinami}}}{MS_{\text{uvnitř skupin}}}. $$ Symbol $MS$ označuje průměrný čtverec (*mean square*): součet čtverců dělený příslušnými stupni volnosti. F-statistika tak spojuje všechny předchozí kroky: rozdíly průměrů mění `r render_glossary_term("odhadnuta-hodnota", display = "odhadnuté hodnoty")`, tím zkracují residua, z residuí počítáme součty čtverců a jejich férově přepočtený poměr testuje společnou nulovou hypotézu. ```{r} #| label: zobrazit-anova-vysledek #| echo: true #| eval: true # Položíme celkovou otázku stejnému třískupinovému modelu anova(mod_tri_druhy) ``` {{< pagebreak >}} ### Jak číst tabulku `anova()` řádek po řádku Tabulka je kompaktní, ale každý údaj odpovídá kroku, který jsme už viděli v rovnici nebo grafu: 1. **Řádek `druh`** popisuje část variability zachycenou přidáním tří skupinových průměrů místo jednoho společného průměru. **Řádek `Residuals`** popisuje variabilitu, která zůstává kolem skupinových průměrů. 2. **`Df` pro `druh` je `r df_group`.** Tři skupinové průměry vyjádříme jedním referenčním průměrem a dvěma rozdíly, takže druh přidává $`r nlevels(data_penguins$druh)`-1=`r df_group`$ stupně volnosti. **`Df` pro residua je `r df_residual`.** Máme `r nobs(mod_three_species)` `r render_glossary_term("pozorovani", display = "pozorování")` a odhadujeme tři skupinové průměry, proto zbývá $`r nobs(mod_three_species)`-`r nlevels(data_penguins$druh)`=`r df_residual`$ residuálních stupňů volnosti. 3. **`Sum Sq` pro `druh` je `r format_cz(sse_reduction, 0)` g².** Je to přesně rozdíl mezi SSE jednoho a tří průměrů a fialová část předchozího rozkladového grafu. **`Sum Sq` pro `Residuals` je `r format_cz(sse_three_means, 0)` g²**, tedy šedá část, která zůstává uvnitř druhů. 4. **`Mean Sq`** vznikne vydělením každého `Sum Sq` jeho `Df`: `r format_cz(ms_group, 0)` g² mezi skupinami a `r format_cz(ms_residual, 0)` g² uvnitř skupin. 5. **`F value` je `r format_cz(f_anova, 2)`.** Je to poměr těchto dvou průměrných čtverců. Hodnota mnohem větší než 1 znamená, že rozdíly skupinových průměrů jsou velké vůči běžné variabilitě jedinců kolem těchto průměrů. 6. **`Pr(>F)` je `r render_glossary_term("p-hodnota", display = "p-hodnota")`.** Nejde o další nezávislý výpočet z dat; vzniká umístěním pozorovaného $F$ do nulového `r render_glossary_term("rozdeleni", display = "F-rozdělení")` s `r df_group` stupni volnosti pro čitatele a `r df_residual` pro jmenovatele. ### Jak se z F-statistiky stane p-hodnota Pokud by všechny tři populační průměry byly stejné a platily předpoklady modelu, opakované `r render_glossary_term("vyber", display = "výběry")` by vytvářely hodnoty z F-rozdělení $F(`r df_group`, `r df_residual`)$. P-hodnota je plocha napravo od naší pozorované hodnoty `r format_cz(f_anova, 2)`: $$ p = P\!\left( F_{`r df_group`,`r df_residual`} \ge `r format_cz_math(f_anova, 2)` \mid H_0 \right). $$ ```{r} #| label: vypocitat-p-hodnotu-z-f #| echo: true #| eval: true # F a stupně volnosti vezmeme z tabulky anova() tabulka_anova <- anova(object = mod_tri_druhy) f_druh <- tabulka_anova["druh", "F value"] df_druh <- tabulka_anova["druh", "Df"] df_residua <- tabulka_anova["Residuals", "Df"] # Plocha F-rozdělení napravo od pozorovaného F pf( q = f_druh, df1 = df_druh, df2 = df_residua, lower.tail = FALSE ) ``` Argument `lower.tail = FALSE` říká R, že chceme pravý chvost: pravděpodobnost hodnot stejně velkých nebo větších než pozorované $F$. U našeho modelu je tato plocha tak malá, že ji v grafu prakticky neuvidíme. ```{r} #| label: zobrazit-f-rozdeleni #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "F-rozdělení očekávané při nulové hypotéze. Levý panel ukazuje běžné hodnoty blízko jedné; pravý panel umisťuje pozorovanou F-statistiku na skutečnou osu. P-hodnota je plocha napravo od oranžové čáry a je zde příliš malá, aby byla viditelná." #| fig-alt: "Dva panely F-rozdělení: vlevo křivka s většinou hodnot blízko jedné, vpravo stejná křivka na široké ose a oranžová čára pozorovaného F daleko vpravo, kde je plocha pod křivkou neviditelná." #| fig-width: 10 #| fig-height: 4.4 ggplot2::ggplot( data = data_f_distribution, mapping = ggplot2::aes(x = f, y = density) ) + ggplot2::geom_line( color = biostat_cols[["graphite"]], linewidth = 1 ) + ggplot2::geom_vline( data = data_f_observed, mapping = ggplot2::aes(xintercept = f), color = biostat_cols[["orange"]], linewidth = 1.1, inherit.aes = FALSE ) + ggplot2::geom_label( data = tibble::tibble( panel = factor( "pozorované F na celé ose", levels = levels(data_f_distribution$panel) ), f = f_anova, density = 0.20, label = paste0("pozorované F = ", format_cz(f_anova, 2)) ), mapping = ggplot2::aes(x = f, y = density, label = label), hjust = 1.05, color = biostat_cols[["orange"]], fill = "white", linewidth = 0.2, size = 3.2, inherit.aes = FALSE ) + ggplot2::facet_wrap( facets = ggplot2::vars(panel), nrow = 1, scales = "free_x" ) + ggplot2::labs( x = "F-statistika očekávaná při nulové hypotéze", y = "Hustota" ) + ggplot2::theme( strip.text = ggplot2::element_text(face = "bold") ) ``` Pro tučňáky tedy vychází $F(`r df_group`, `r df_residual`)=`r format_cz_math(f_anova, 2)`$ a p-hodnota je `r format_p_cz(p_anova)`. Za předpokladů modelu jsou data velmi obtížně slučitelná s hypotézou, že všechny tři druhy mají stejnou průměrnou hmotnost. ::: {.box-warning title="Celkový výsledek neříká, že se liší každá dvojice"} Malá p-hodnota analýzy rozptylu podporuje tvrzení, že alespoň jeden průměr se liší. Neurčuje však, které dvojice se liší, ani zda jsou všechny rozdíly biologicky důležité. ::: ## `r render_glossary_term("post-hoc-test", display = "Post-hoc")`: které dvojice se liší? F-test odpověděl na celkovou otázku: data jsou obtížně slučitelná s hypotézou tří stejných `r render_glossary_term("prumer", display = "průměrů")`. Neříká ale, které druhy se liší. Na to se ptají navazující otázky o konkrétních dvojicích. U tří druhů existují tři páry: Chinstrap minus Adélie, Gentoo minus Adélie a Gentoo minus Chinstrap. Kdybychom každý pár testovali jako izolovanou novou studii, s rostoucím počtem testů by rostla šance na alespoň jeden náhodně extrémní výsledek. Tukeyho post-hoc metoda proto vyhodnotí celou rodinu párových otázek společně. Získáme `r render_glossary_term("interval-spolehlivosti", display = "simultánní intervaly")` a `r render_glossary_term("p-hodnota", display = "p-hodnoty")` upravené pro všechna tři porovnání. ::: {.box-rule title="Dvě různé otázky položené stejnému modelu"} - **Celkový F-test** se ptá, zda jsou **všechny skupinové průměry stejné**. Odpovídá jedním číslem pro všechny skupiny. - **Post-hoc porovnání** se ptají, **které konkrétní dvojice** se liší a o kolik gramů. Tukeyho úprava sama zohlední, že se díváme na více dvojic současně, a nepotřebuje proto předchozí „významný“ F-test jako vstupenku. Párová porovnání dávají smysl, když máme **tři a více skupin** a biologická otázka se týká konkrétních dvojic. U dvou skupin je jediná dvojice a odpověď už dává koeficient modelu. Když jsou párové rozdíly malé a jejich intervaly obsahují nulu, poctivý závěr zní, že data rozdíly neodlišují od nuly. Problémem by bylo vybrat si dodatečně jednu „zajímavou“ dvojici a uvést její neupravenou p-hodnotu, jako by šlo o jedinou předem položenou otázku. ::: ::: {.box-data collapse="true" title="Doplňující: proč jsou Tukeyho intervaly širší"} Obyčejný pětadevadesátiprocentní interval chrání dlouhodobé pokrytí jedné předem vybrané veličiny. Zde ale chceme současně číst tři rozdíly. Kdyby měl každý neopravený interval pětiprocentní riziko minutí, šance, že alespoň jeden z celé sady mine svou skutečnou hodnotu, by byla větší než pět procent. Tukeyho metoda proto nastaví jednotlivé intervaly opatrněji tak, aby přibližně devadesát pět procent opakovaných **celých sad tří intervalů** pokrylo všechny tři skutečné rozdíly současně. Za tuto společnou ochranu obvykle platíme širšími intervaly a většími upravenými p-hodnotami než při izolovaném čtení jediné dvojice. ::: ::: {.box-data collapse="true" title="Doplňující: plánované kontrasty a post-hoc porovnání nejsou stejná otázka"} V této lekci klademe celkovou otázku a k ní systematicky prohlížíme všechny tři dvojice. To je post-hoc postup. Výzkumník však někdy ještě před sběrem dat určí jeden biologicky motivovaný kontrast, například „potápěčský druh proti průměru dvou pobřežních druhů“. Taková otázka není dodatečným hledáním nejzajímavější dvojice. Předem plánovaný kontrast může být věcně přesnější než všechna párová srovnání a stejně jako Tukeyho porovnání nemusí čekat na výsledek celkového F-testu. Stále však musíme transparentně uvést všechny plánované otázky a při více kontrastech zvážit jejich společnou chybovost. Výpočet kontrastů patří do pokročilejší práce s modely. ::: Funkce `emmeans()` pochází z balíčku `emmeans`; pokud jej nemáte, nainstalujte jej jednou příkazem `install.packages("emmeans")` v Console. Pracuje přímo s již fitovaným objektem `mod_tri_druhy`. **Nefitujeme nový model.** Pouze z téhož `r render_glossary_term("linearni-model", display = "lineárního modelu")` získáme průměrnou hmotnost odhadnutou pro každý druh. V tomto jednoduchém modelu s jediným `r render_glossary_term("kategoricka-promenna", display = "kategoriálním prediktorem")` jsou tyto hodnoty stejné jako skupinové průměry, které už známe z grafu. ```{r} #| label: odhadnout-prumery-emmeans-viditelne #| echo: true #| eval: true # Ze stejného modelu získáme odhadnutý průměr každého druhu prumery_druhu <- emmeans::emmeans( object = mod_tri_druhy, specs = ~ druh ) prumery_druhu ``` Výstup objektu `prumery_druhu` ukazuje všechny tři průměry odhadnuté modelem. Sloupec `emmean` obsahuje průměrnou hmotnost: `r format_cz(as.data.frame(prumery_druhu)$emmean[1], 1)` g pro Adélie, `r format_cz(as.data.frame(prumery_druhu)$emmean[2], 1)` g pro Chinstrap a `r format_cz(as.data.frame(prumery_druhu)$emmean[3], 1)` g pro Gentoo. Sloupec `SE` je `r render_glossary_term("standardni-chyba", display = "standardní chyba")` každého průměru, `df` udává residuální `r render_glossary_term("stupne-volnosti", display = "stupně volnosti")` a `lower.CL` s `upper.CL` vymezují 95% interval pro každý skupinový průměr zvlášť. Tyto intervaly popisují `r render_glossary_term("nejistota", display = "nejistotu")` jednotlivých průměrů; samy ještě nejsou intervaly jejich vzájemných rozdílů. Ty vzniknou až porovnáním odhadnutých průměrů mezi sebou. Metoda `revpairwise` zachová směr, který jsme používali v grafech: Chinstrap minus Adélie, Gentoo minus Adélie a Gentoo minus Chinstrap. Argument `adjust = "tukey"` říká, že intervaly a p-hodnoty mají zohlednit celou rodinu tří porovnání. {{< pagebreak >}} ```{r} #| label: post-hoc-emmeans-viditelne #| echo: true #| eval: true # Ze stejných odhadnutých průměrů vytvoříme tři párové rozdíly post_hoc_druhy <- emmeans::contrast( object = prumery_druhu, method = "revpairwise", adjust = "tukey" ) summary( object = post_hoc_druhy, # První TRUE zobrazí intervaly, druhé t-testy a p-hodnoty infer = c(TRUE, TRUE) ) ``` Ve výstupu sledujeme řádky jednotlivých dvojic. Sloupec `estimate` je odhadnutý rozdíl průměrů, `lower.CL` a `upper.CL` jsou meze simultánního intervalu a `p.value` je p-hodnota upravená pro společnou rodinu porovnání. Poznámky pod výstupem potvrzují, že Tukeyho korekce byla použita pro intervaly i p-hodnoty. Přehledněji si stejné hodnoty uspořádáme do tabulky. ```{r} #| label: spocitat-post-hoc-porovnani #| echo: false #| output: false #| message: false #| warning: false #| error: false data_tukey <- summary( object = post_hoc_druhy, infer = c(TRUE, TRUE) ) |> as.data.frame() |> dplyr::transmute( porovnani = factor( contrast, levels = c( "Chinstrap - Adelie", "Gentoo - Adelie", "Gentoo - Chinstrap" ), labels = c( "Chinstrap minus Adélie", "Gentoo minus Adélie", "Gentoo minus Chinstrap" ) ), odhad = estimate, dolni = lower.CL, horni = upper.CL, p_adjusted = p.value ) ``` ```{r} #| label: zobrazit-post-hoc-tabulku #| echo: false #| output: true #| message: false #| warning: false #| error: false data_tukey |> dplyr::transmute( Porovnání = porovnani, `Rozdíl (g)` = round(odhad, digits = 1), `Dolní mez (g)` = round(dolni, digits = 1), `Horní mez (g)` = round(horni, digits = 1), `Upravená p-hodnota` = format_p_cz(p_adjusted) ) |> knitr::kable( caption = paste0( "Tukeyho post-hoc rozdíly, simultánní intervaly a p-hodnoty ", "upravené pro tři porovnání." ) ) ``` Tabulka zachovává hlavní čísla z výstupu `summary(post_hoc_druhy)`. Graf je ukazuje na jedné ose: bod ukáže odhad rozdílu, vodorovná úsečka jeho simultánní interval a svislá čára nulu. ```{r} #| label: zobrazit-post-hoc-intervaly #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Tři Tukeyho post-hoc rozdíly a simultánní 95% intervaly. Nula znamená stejné průměry dané dvojice." #| fig-alt: "Tři fialové body rozdílů s oranžovými úsečkami intervalů: interval Chinstrap minus Adélie protíná čárkovanou nulu, oba rozdíly s Gentoo leží daleko vpravo od nuly." #| fig-width: 9 #| fig-height: 4.5 ggplot2::ggplot( data = data_tukey, mapping = ggplot2::aes( x = odhad, y = porovnani ) ) + ggplot2::geom_vline( xintercept = 0, linetype = 2, color = biostat_cols[["grey_olive"]] ) + ggplot2::geom_errorbar( mapping = ggplot2::aes( xmin = dolni, xmax = horni ), orientation = "y", width = 0.14, linewidth = 1, color = biostat_cols[["orange"]] ) + ggplot2::geom_point( size = 3.2, color = biostat_cols[["indigo_velvet"]] ) + ggplot2::labs( x = "Rozdíl průměrné hmotnosti (g)", y = NULL ) ``` Rozdíl Chinstrap minus Adélie je `r format_cz(data_tukey$odhad[1], 1)` g a simultánní interval od `r format_cz(data_tukey$dolni[1], 1)` do `r format_cz(data_tukey$horni[1], 1)` g protíná nulu. Data tedy tento malý rozdíl neodlišují od nuly. Oba rozdíly zahrnující Gentoo jsou přibližně `r format_cz(data_tukey$odhad[3], 0)` až `r format_cz(data_tukey$odhad[2], 0)` g a jejich intervaly leží celé nad nulou. ::: {.box-result title="Celková a párová otázka dávají různé informace"} `r render_glossary_term("anova", display = "Analýza rozptylu")` ukázala, že všechny tři průměry nejsou stejné. Párová porovnání upřesnila proč: Gentoo se zřetelně liší od obou dalších druhů, zatímco Adélie a Chinstrap mají v tomto `r render_glossary_term("vyber", display = "výběru")` velmi podobné průměrné hmotnosti. ::: ## Co jednoduchý model neoddělil Rozdíly mezi druhy jsou v těchto datech **asociace**. Jedinci nebyli náhodně přiřazeni ke druhům a druhy nebyly na ostrovech zastoupeny stejně. ```{r} #| label: ukazat-druhy-a-ostrovy #| echo: false #| output: true #| message: false #| warning: false #| error: false data_penguins |> dplyr::count(druh, ostrov) |> tidyr::pivot_wider( names_from = ostrov, values_from = n, values_fill = 0 ) |> dplyr::mutate( Druh = popisky_druhu[as.character(druh)] ) |> dplyr::select(Druh, Biscoe, Dream, Torgersen) |> knitr::kable( caption = "Počet jedinců jednotlivých druhů na třech ostrovech." ) ``` Všichni tučňáci uzdičkoví v tomto `r render_glossary_term("vyber", display = "výběru")` pocházejí z ostrova Dream a všichni tučňáci oslí z ostrova Biscoe. Adélie byli měřeni na všech třech ostrovech. Jednoduchý model proto neumí oddělit druh od prostředí ostrova. Pohlaví také souvisí s hmotností. Následující graf tuto strukturu pouze zviditelňuje; stále nefitujeme model s více `r render_glossary_term("prediktor", display = "prediktory")`. ```{r} #| label: zobrazit-pohlavi-uvnitr-druhu #| echo: false #| output: true #| message: false #| warning: false #| error: false #| fig-cap: "Hmotnosti tří druhů rozlišené podle zaznamenaného pohlaví. Šedé body mají pohlaví nezaznamenané." #| fig-alt: "Body hmotností tří druhů rozlišené barvou a tvarem podle pohlaví; uvnitř každého druhu leží samci převážně výše než samice." #| fig-width: 9 #| fig-height: 5.5 data_sex_plot <- data_penguins |> dplyr::mutate( pohlavi_plot = dplyr::case_when( pohlavi == "female" ~ "samice", pohlavi == "male" ~ "samec", TRUE ~ "nezaznamenáno" ) ) ggplot2::ggplot( data = data_sex_plot, mapping = ggplot2::aes( x = druh, y = hmotnost_tela_g, shape = pohlavi_plot, color = pohlavi_plot ) ) + ggbeeswarm::geom_quasirandom( width = 0.24, alpha = 0.58, size = 1.6 ) + ggplot2::scale_color_manual( values = c( # Sex is a grouping variable: _brand categorical palette. "samice" = "#C2613A", "samec" = "#2E6CA4", "nezaznamenáno" = biostat_cols[["grey_olive"]] ) ) + ggplot2::scale_x_discrete(labels = popisky_druhu) + ggplot2::labs( x = NULL, y = "Hmotnost těla (g)", shape = "Pohlaví", color = "Pohlaví" ) ``` ## Co tedy můžeme říct o tučňácích? Správný závěr nezačíná `r render_glossary_term("p-hodnota", display = "p-hodnotou")`. Jeho jádrem je velikost rozdílů v gramech a jejich `r render_glossary_term("nejistota", display = "nejistota")`; testová informace a omezení studie je doplňují. ::: {.box-question title="Sestavte závěr z pěti částí"} Napište krátký odstavec v pořadí: biologická otázka → rozdíly v gramech → `r render_glossary_term("interval-spolehlivosti", display = "intervaly")` → celkový a párový test → omezení zobecnění a kauzality. ::: ::: {.box-result collapse="true" title="Příklad opatrného závěru"} Ptali jsme se, zda se liší průměrná hmotnost tří druhů tučňáků. V souboru `r count_complete` dospělých tučňáků z Palmerova souostroví byli tučňáci oslí (Gentoo) v průměru přibližně o `r format_cz(data_tukey$odhad[2], 0)` g těžší než tučňáci kroužkoví (Adélie) a o `r format_cz(data_tukey$odhad[3], 0)` g těžší než tučňáci uzdičkoví (Chinstrap); simultánní 95% intervaly obou rozdílů ležely celé nad nulou (od `r format_cz(data_tukey$dolni[2], 0)` do `r format_cz(data_tukey$horni[2], 0)` g a od `r format_cz(data_tukey$dolni[3], 0)` do `r format_cz(data_tukey$horni[3], 0)` g). Průměry Adélie a Chinstrap se lišily jen o `r format_cz(data_tukey$odhad[1], 0)` g a jejich interval od `r format_cz(data_tukey$dolni[1], 0)` do `r format_cz(data_tukey$horni[1], 0)` g zahrnoval nulu, takže data tento malý rozdíl neodlišují od nuly. Celková analýza rozptylu byla s hypotézou tří stejných průměrů obtížně slučitelná ($F(`r df_group`, `r df_residual`) = `r format_cz_math(f_anova, 2)`$; p-hodnota `r format_p_cz(p_anova)`). Jde o observační data: druh je v souboru propojený s ostrovem a hmotnost souvisí také s pohlavím, proto rozdíly nelze vykládat jako přímý účinek druhu a závěr platí pro sledované kolonie. ::: ## Co bude následovat dál Model s jediným kategoriálním prediktorem odpovídá na jednu otázku: jak se liší průměry druhů. Neumí ale oddělit druh od ostrova nebo pohlaví. Změnil by se odhad rozdílu mezi druhy, kdybychom porovnávali jedince stejného pohlaví nebo současně zohlednili ostrov? Na to potřebujeme model s více `r render_glossary_term("prediktor", display = "prediktory")`, kterým se zabývá následující lekce. ## Shrnutí - `r render_glossary_term("kategoricka-promenna", display = "Kategoriální prediktor")` nemění `r render_glossary_term("linearni-model", display = "lineární model")` v úplně jinou metodu; mění způsob, jakým model vytváří `r render_glossary_term("odhadnuta-hodnota", display = "odhadnuté hodnoty")`. - Ve dvouskupinovém modelu je `r render_glossary_term("intercept", display = "intercept")` `r render_glossary_term("prumer", display = "průměr")` referenční skupiny a druhý koeficient rozdíl druhé skupiny vůči referenci. - Odhadnuté hodnoty jsou skupinové průměry a `r render_glossary_term("residuum", display = "residua")` jsou vzdálenosti jednotlivých `r render_glossary_term("pozorovani", display = "pozorování")` od jejich skupinového průměru. - `r render_glossary_term("t-test", display = "Dvouvýběrový t-test")` se společným `r render_glossary_term("rozptyl", display = "rozptylem")` je test koeficientu dvouskupinového lineárního modelu. - Rozdíl mezi průměry posuzujeme vzhledem k `r render_glossary_term("variabilita", display = "variabilitě")` uvnitř skupin. Na tomto porovnání stojí `r render_glossary_term("anova", display = "analýza rozptylu")`: celkový F-test říká, zda jsou všechny průměry slučitelné se stejnou hodnotou, ale ne, která dvojice se liší. - `r render_glossary_term("post-hoc-test", display = "Post-hoc porovnání")` odpovídají na navazující otázky o konkrétních dvojicích a s Tukeyho úpravou zohledňují, že zkoumáme více dvojic současně. - Biologický závěr začíná rozdíly v gramech a jejich `r render_glossary_term("nejistota", display = "nejistotou")` a končí předpoklady a omezeními studie. ### Jak lekce navazuje na předchozí | Dřívější myšlenka | Jak se objevila v této lekci | |---|---| | Typ proměnné určuje vhodný popis | Numerickou hmotnost jsme vysvětlovali kategoriálním druhem. | | Graf předchází modelu | Jednotlivci, překryv a variabilita byli vidět dříve než průměry a koeficienty. | | $y_i=\hat{y}_i+e_i$ | Odhadnuté hodnoty se změnily z bodů na přímce na skupinové průměry; residua zůstala vzdálenostmi od odhadu. | | Efekt a jeho nejistota jsou dvě informace | Rozdíly v gramech jsme doplnili intervaly spolehlivosti. | | Test je otázka o parametru | Dvouvýběrový t-test se stal testem koeficientu a analýza rozptylu celkovou otázkou o všech skupinových průměrech. | | Závěr není jen p-hodnota | Výsledek spojil velikost rozdílů, intervaly, testy, diagnostiku a omezení studie. | ### Závěrečná otázka ::: {.box-question title="Co potřebujete vědět, než napíšete „druh určuje hmotnost tučňáka“?"} Představte si, že vám někdo ukáže jen tuto větu a malou p-hodnotu analýzy rozptylu. Jaké tři otázky mu položíte? ::: ::: {.box-result collapse="true" title="Možná odpověď"} O kolik gramů se druhy liší a jak široké jsou intervaly těchto rozdílů? Které konkrétní dvojice se liší, a které ne (Adélie a Chinstrap jsou si velmi podobné)? Sedí model k datům, tedy jak vypadají residua? A je možné rozdíl připsat druhu, když jsou druhy v souboru propojené s ostrovy a hmotnost souvisí i s pohlavím? Slovo „určuje“ navíc naznačuje příčinu, kterou observační data sama neprokážou. ::: ## 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/" ) ) } ```