--- title: "O kolik déle čekáme po delší erupci — a jak jistě to víme?" subtitle: "Od odhadu k nejistotě: výběrová variabilita a interval spolehlivosti" lang: cs filters: - ../theme/semantic-boxes.lua - glossary-headings.lua format: html: toc: true toc-depth: 3 toc-title: "Obsah" number-sections: true code-fold: false 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/07/21 date-format: long date-modified: last-modified --- ```{r} #| label: nastaveni #| include: false if (!require("here")) install.packages("here") here::i_am("Learning_materials/skripta.qmd") # Visible code reads data/old_faithful_2024.csv relative to the project root, # exactly as students do in their own R project. knitr::opts_knit$set(root.dir = here::here()) # Hidden figures use Czech decimal commas on their axes. Visible chunks keep # R's default output so that it matches what students see in R. knitr::knit_hooks$set( desetinna_carka = function(before) { if (before) { options(OutDec = ",") } else { options(OutDec = ".") } } ) knitr::opts_hooks$set( echo = function(options) { if (isFALSE(options$echo)) { options$desetinna_carka <- TRUE } options } ) if (!require("tidyverse")) install.packages("tidyverse") if (!require("glossary")) install.packages("glossary") library(tidyverse) library(glossary) glossary::glossary_popup("hover") 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", "simulate_reports.R")) source(here::here("R", "Functions", "summarise_reports.R")) source(here::here("R", "Functions", "simulate_stability_pairs.R")) source(here::here("R", "Functions", "plot_stability_pairs.R")) source(here::here("R", "Functions", "simulate_influence_gallery.R")) source(here::here("R", "Functions", "plot_influence_gallery.R")) source(here::here("R", "Functions", "simulate_band_gallery.R")) source(here::here("R", "Functions", "plot_band_gallery.R")) # Glossary: use local slovnik/pojmy.yaml when available, otherwise remote fallback. # Terms are looked up by slug ID (e.g. "histogram", "median", "boxplot"). # Usage in text: `r render_glossary_term("slug", display = "Czech term")` local_pojmy <- here::here("..", "slovnik", "pojmy.yaml") if (file.exists(local_pojmy)) { glossary::glossary_path(local_pojmy) } else { remote_pojmy <- tempfile(fileext = ".yaml") tryCatch( { download.file( url = "https://raw.githubusercontent.com/CUNI-NATUR-Biostatistics/slovnik/main/pojmy.yaml", destfile = remote_pojmy, quiet = TRUE ) glossary::glossary_path(remote_pojmy) }, error = function(e) { message("Glossary YAML unavailable: ", conditionMessage(e)) } ) } tryCatch( { source( paste0( "https://raw.githubusercontent.com/", "CUNI-NATUR-Biostatistics/slovnik/main/", "R/Functions/render_glossary_term.R" ) ) }, error = function(e) { source(here::here("R/Functions/render_glossary_term.R")) } ) ``` ```{r} #| label: nastaveni-slovnik #| echo: false #| message: false #| results: asis if (knitr::is_html_output()) { glossary::glossary_style( color = biostat_cols[["indigo_velvet"]], text_decoration = "underline", def_bg = biostat_cols[["white"]], def_color = biostat_cols[["graphite"]] ) } ``` ## Úvod Gejzír Old Faithful v Yellowstonském národním parku střídá erupce a období klidu. Podle [National Park Service](https://www.nps.gov/places/old-faithful-geyser.htm) trvá erupce přibližně 1,5 až 5 minut a čekání na další erupci se obvykle pohybuje od 50 do 127 minut. Návštěvníci i správa parku se proto ptají, zda délka právě skončené erupce napoví, jak dlouho budeme čekat na další. V předchozí lekci jsme se naučili odhadnout `r render_glossary_term("linearni-model", display = "lineárním modelem")`, o kolik se v průměru liší `r render_glossary_term("odezva", display = "odezva")` dvou `r render_glossary_term("pozorovani", display = "pozorování")`, jejichž `r render_glossary_term("prediktor", display = "prediktor")` se liší o jednu jednotku. `r render_glossary_term("odhad", display = "Odhad")` z jednoho souboru dat ale není celá odpověď. Kdyby gejzír pozoroval někdo jiný v jiných dnech, zaznamenal by jiné erupce a jeho odhad by mohl vyjít jinak. Potřebujeme proto znát nejen velikost `r render_glossary_term("efekt", display = "efektu")`, ale také to, jak přesně jsme ji odhadli. ::: {.box-question title="Hlavní otázka této lekce"} **O kolik déle čekáme po delší erupci — a jak jistě to víme?** ::: ::: {.box-question title="Co čekáte ještě před prvním grafem?"} Představte si, že gejzír pozorují další lidé v jiných dnech a každý z nich spočítá vlastní `r render_glossary_term("sklon", display = "sklon")` přímky. Jaké sklony dostanou? - **A:** Přesně stejný sklon jako my. - **B:** Podobné, ale ne stejné sklony. - **C:** Zcela odlišné sklony, které spolu nesouvisejí. Než budete číst dál, vyberte možnost, která vám připadá nejpravděpodobnější, a zdůvodněte ji. ::: ::: {.box-result collapse="true" title="Odpověď: podobné, ale ne stejné sklony"} Nejpravděpodobnější je **možnost B**. Každý pozorovatel zaznamená jiné erupce. Pokud vztah mezi délkou erupce a čekáním skutečně existuje, budou si odhady podobné; protože však body nikdy neleží přesně na přímce, budou se navzájem trochu lišit. Tato lekce ukazuje, jak velký takový rozdíl bývá a jak jej popsat čísly. ::: ## Výsledky učení Po prostudování této lekce budete schopni: - vysvětlit, proč jiný soubor `r render_glossary_term("pozorovani", display = "pozorování")` dává jiný `r render_glossary_term("odhad", display = "odhad")` `r render_glossary_term("sklon", display = "sklonu")`; - najít ve výstupu `summary()` odhad sklonu a jeho `r render_glossary_term("standardni-chyba", display = "standardní chybu")`, vysvětlit, co standardní chybu zmenšuje, a rozlišit velikost vztahu od přesnosti odhadu; - získat `r render_glossary_term("interval-spolehlivosti", display = "95% interval spolehlivosti")` sklonu funkcí `confint()` a `r render_glossary_term("pas-spolehlivosti", display = "pás spolehlivosti")` kolem přímky funkcí `predict()` a interpretovat je bez tvrzení, že obsahují 95 % pozorování nebo dokazují hypotézu; - posoudit pomocí `r render_glossary_term("pakovy-efekt", display = "pákového efektu")` a `r render_glossary_term("cookova-vzdalenost", display = "Cookovy vzdálenosti")`, zda jedna erupce výrazně mění odhad sklonu; - zapsat výsledek jednou větou s odhadem, standardní chybou, intervalem a jednotkami. ## Každý řádek je jedna erupce Používáme skutečná `r render_glossary_term("pozorovani", display = "pozorování")` erupcí Old Faithful z léta 2024, která zaznamenali dobrovolníci prostřednictvím služby [GeyserTimes](https://geysertimes.org/). Jeden řádek tabulky odpovídá jedné erupci a čekání na erupci následující; erupce je tedy `r render_glossary_term("jednotka-pozorovani", display = "jednotkou pozorování")`. Soubor [old_faithful_2024.csv](../data/old_faithful_2024.csv) je součástí materiálů. Pro následující kód jej uložte do složky `data` ve svém projektu v R. Délka erupce je v souboru v sekundách; pro výklad ji převedeme na minuty, aby obě měření měla stejnou jednotku: ```{r} #| label: data-gejzir-prehled #| echo: true # Načteme tabulku erupcí ze souboru CSV data_gejzir <- read.csv(file = "data/old_faithful_2024.csv") # Datum převedeme z textu na kalendářní datum data_gejzir$date <- as.Date(x = data_gejzir$date) # Délku erupce převedeme ze sekund na minuty data_gejzir$delka_erupce_min <- data_gejzir$duration / 60 # Zkontrolujeme názvy a typy sloupců str(object = data_gejzir) ``` V lekci budeme pracovat se dvěma `r render_glossary_term("numericka-promenna", display = "numerickými proměnnými")`: - `delka_erupce_min`: délka právě pozorované erupce v minutách — náš `r render_glossary_term("prediktor", display = "prediktor")`; - `cekani_min`: čekání do následující erupce v minutách — naše `r render_glossary_term("odezva", display = "odezva")`. Sloupec `date` použijeme k rozdělení pozorování podle měsíců. Ostatní sloupce popisují identifikátor erupce, čas a to, zda erupci zachytila webkamera. Kontrola dat pokračuje rozsahy obou měření, chybějícími hodnotami a počtem erupcí v jednotlivých měsících: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: data-gejzir-kontrola #| echo: true # Minimum, kvartily, průměr a maximum obou měření summary( object = data_gejzir[, c("delka_erupce_min", "cekani_min")] ) # Počet chybějících hodnot v každém sloupci colSums(x = is.na(x = data_gejzir)) # Počet erupcí v jednotlivých měsících (06 = červen) table(format(x = data_gejzir$date, format = "%m")) ``` ```{=typst} ] ``` ```{r} #| label: spocitat-erupce-po-mesicich #| include: false vec_pocty_mesicu <- as.vector(table(format(x = data_gejzir$date, format = "%m"))) ``` Tabulka obsahuje `r format_cz(nrow(data_gejzir))` erupcí od `r format(x = min(data_gejzir$date), format = "%d. %m. %Y")` do `r format(x = max(data_gejzir$date), format = "%d. %m. %Y")`: v červnu, červenci a srpnu postupně `r paste(format_cz(vec_pocty_mesicu), collapse = ", ")`. Erupce trvaly od `r format_cz(min(data_gejzir$delka_erupce_min))` do `r format_cz(max(data_gejzir$delka_erupce_min))` minuty a čekání se pohybovalo od `r format_cz(min(data_gejzir$cekani_min))` do `r format_cz(max(data_gejzir$cekani_min))` minut; žádný sloupec nemá chybějící hodnoty. Oba rozsahy leží uvnitř hodnot uváděných National Park Service, takže nic nenaznačuje chybu jednotek. ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: tabulka-prvnich-erupci #| echo: false #| message: false #| warning: false #| tbl-cap: "Prvních šest erupcí: datum, délka erupce a čekání na další erupci." data_gejzir |> dplyr::transmute( Datum = format(x = date, format = "%d. %m. %Y"), `Délka erupce (min)` = format_cz(delka_erupce_min), `Čekání (min)` = format_cz(cekani_min) ) |> head(n = 6) |> knitr::kable(align = "lrr") ``` ```{=typst} ] ``` ::: {.box-data collapse="true" title="Doplňující: jak vznikl soubor old_faithful_2024.csv"} Erupce Old Faithful zaznamenávají dobrovolníci v databázi [GeyserTimes](https://geysertimes.org/). Tabulku z roku 2024 připravili autoři učebnice ModernDive jako tabulku `old_faithful_2024` v balíčku [moderndive](https://moderndive.github.io/moderndive/reference/old_faithful_2024.html). Pro tuto lekci jsme stejný postup zopakovali přímo z [oficiálního archivu erupcí Old Faithful](https://geysertimes.org/archive/geysers/Old_Faithful_eruptions.tsv.gz): - ponechali jsme pouze přesně zaznamenané hlavní erupce z června, července a srpna 2024; - čekání jsme spočítali jako rozdíl časů dvou po sobě jdoucích erupcí a ponechali intervaly od 50 do 127 minut, což odpovídá [rozsahu uváděnému National Park Service](https://www.nps.gov/places/old-faithful-geyser.htm); delší mezera v dobrovolnických záznamech pravděpodobně znamená, že některá erupce nebyla zaznamenána; - odstranili jsme tři erupce, které vypadaly jako `r render_glossary_term("odlehla-pozorovani", display = "odlehlá pozorování")`, jednu z každého měsíce; jejich identifikátory jsou uvedeny přímo v kódu. Erupce dlouhá jen 41 sekund leží pod rozsahem 1,5 až 5 minut uváděným National Park Service. Zbylé dvě jsme odstranili jen kvůli přehlednějšímu výukovému příkladu; v kapitole „Může jedna erupce změnit sklon?“ je vrátíme zpět a ukážeme, že sklon téměř neovlivňují, takže bychom je po takovém ověření v datech ponechali. Základní postup lekce si můžete vyzkoušet i s klasickou tabulkou `faithful`, která je součástí R. Neobsahuje však data erupcí, takže s ní nelze zopakovat rozdělení podle měsíců: ```{r} #| label: alternativa-faithful-v-r #| echo: true #| eval: false # Přejmenujeme sloupce klasické tabulky podle této lekce data_faithful <- data.frame( delka_erupce_min = faithful$eruptions, cekani_min = faithful$waiting ) ``` Úplný záznam stažení a přípravy dat používá i pokročilejší kroky, které pro lekci není nutné rozumět ani spouštět. Pokud GeyserTimes požádá o ověření v prohlížeči, uložte archiv ručně do složky `data` pod názvem `Old_Faithful_eruptions.tsv.gz`. ```{r} #| label: stahnout-a-pripravit-data-gejziru #| echo: true #| eval: false #| code-fold: true #| code-summary: "Ukázat kód pro stažení a přípravu dat" # Sestavíme adresu oficiálního archivu Old Faithful archive_url <- paste0( "https://geysertimes.org/archive/geysers/", "Old_Faithful_eruptions.tsv.gz" ) # Určíme, kam se stažený soubor uloží v projektu archive_file <- here::here( "data", "Old_Faithful_eruptions.tsv.gz" ) # Soubor stáhneme jen tehdy, když jej ještě nemáme if ( !file.exists(archive_file) ) { download.file( url = archive_url, destfile = archive_file, mode = "wb" ) } # Načteme archiv, pojmenujeme potřebné sloupce a převedeme čas data_archive <- readr::read_tsv( file = gzfile(description = archive_file), na = c("", "NULL"), quote = "", show_col_types = FALSE ) |> dplyr::transmute( eruption_id = as.numeric(eruptionID), primary_id = as.numeric(associated_primaryID), geyser, datetime = lubridate::as_datetime( x = as.numeric(eruption_time_epoch), tz = "America/Denver" ), exact = exact == 1, near_start = ns == 1, in_eruption = ie == 1, questionable = q == 1, webcam = wc == 1, duration = as.numeric(duration_seconds) ) |> # Ponecháme pouze přesně zaznamenané hlavní erupce dplyr::filter( is.na(primary_id) | eruption_id == primary_id, exact, !near_start, !in_eruption, !questionable ) |> dplyr::arrange(datetime) |> # Z časů dvou sousedních erupcí vypočítáme dobu čekání dplyr::mutate( waiting = as.numeric( difftime( time1 = dplyr::lead(datetime), time2 = datetime, units = "mins" ) ) ) # Vybereme tři měsíce a proměnné potřebné pro výuku data_export <- data_archive |> dplyr::filter( datetime >= lubridate::ymd( "2024-06-01", tz = "America/Denver" ), datetime < lubridate::ymd( "2024-09-01", tz = "America/Denver" ), !is.na(duration), dplyr::between(waiting, left = 50, right = 127) ) |> dplyr::transmute( eruption_id, geyser, date = lubridate::as_date(datetime), time = lubridate::hour(datetime) * 100 + lubridate::minute(datetime), cekani_min = waiting, webcam = dplyr::if_else(webcam, "Yes", "No"), duration ) # Odstraníme tři předem prověřená odlehlá pozorování outlier_ids <- c( 1461337, # červen 1468130, # červenec 1471579 # srpen ) data_export <- data_export |> dplyr::filter(!eruption_id %in% outlier_ids) # Ověříme, že příprava dat skončila očekávaným počtem řádků expected_rows <- 117L if ( nrow(data_export) != expected_rows ) { stop("Zdrojový soubor má neočekávaný počet řádků.") } # Uložíme stabilní výukovou tabulku readr::write_csv( x = data_export, file = here::here("data", "old_faithful_2024.csv"), na = "" ) ``` Pro výuku používáme samostatný soubor CSV, protože v praktiku stačí načíst jediný soubor pomocí `read.csv()`. Původ dat a podmínky jejich použití shrnuje také soubor `data/README.md`. ::: ## Delší erupce, delší čekání? ```{r} #| label: vztah-delky-a-cekani #| echo: false #| message: false #| warning: false #| fig-cap: "Každý bod je jedna zaznamenaná erupce: její délka a čekání na další erupci." #| fig-alt: "Bodový graf všech erupcí. Délka erupce je na vodorovné ose, čekání na další erupci na svislé. Body tvoří dvě skupiny, krátké erupce kolem dvou minut a dlouhé kolem čtyř minut, a celkově stoupají." ggplot2::ggplot( data = data_gejzir, mapping = ggplot2::aes( x = delka_erupce_min, y = cekani_min ) ) + ggplot2::geom_point( color = biostat_cols[["grey_olive"]], alpha = 0.65, size = 2.3 ) + ggplot2::labs( x = "Délka erupce (min)", y = "Čekání na další erupci (min)" ) + theme_biostat(base_size = 13) ``` ::: {.box-question title="Co v grafu vidíte?"} 1. Jaký směr má vztah mezi délkou erupce a čekáním? 2. Které části grafu obsahují nejvíce bodů? 3. Leží body na jediné přímce? ::: ::: {.box-result collapse="true" title="Odpověď: rostoucí vztah s rozptýlenými body"} Vztah je rostoucí: po delších erupcích zpravidla čekáme déle. Body tvoří dvě skupiny — krátké erupce kolem dvou minut a dlouhé erupce kolem čtyř minut — a ani v jedné neleží na jediné přímce. Při stejné délce erupce se čekání liší o několik minut. ::: Stejný vztah shrneme `r render_glossary_term("linearni-model", display = "lineárním modelem")`, jak jsme to dělali v předchozí lekci: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: fitovat-model-vsech-erupci #| echo: true # Odhadneme čekání podle délky erupce ze všech erupcí mod_gejzir <- lm( formula = cekani_min ~ delka_erupce_min, data = data_gejzir ) # Zobrazíme odhadnutý intercept a sklon coef(object = mod_gejzir) ``` ```{=typst} ] ``` ```{r} #| label: pripravit-odhad-vsech-erupci #| echo: false odhad_sklonu <- coef(object = mod_gejzir)[["delka_erupce_min"]] ``` ```{r} #| label: graf-primky-vsech-erupci #| echo: false #| message: false #| warning: false #| fig-cap: "Stejné erupce s přímkou odhadnutou ze všech zaznamenaných erupcí." #| fig-alt: "Stejný bodový graf jako výše s fialovou rostoucí přímkou lineárního modelu." ggplot2::ggplot( data = data_gejzir, mapping = ggplot2::aes( x = delka_erupce_min, y = cekani_min ) ) + ggplot2::geom_point( color = biostat_cols[["grey_olive"]], alpha = 0.65, size = 2.3 ) + # No band: the confidence band is introduced later in this lesson. ggplot2::geom_abline( intercept = coef(object = mod_gejzir)[["(Intercept)"]], slope = odhad_sklonu, color = biostat_cols[["indigo_velvet"]], linewidth = 1.1 ) + ggplot2::labs( x = "Délka erupce (min)", y = "Čekání na další erupci (min)" ) + theme_biostat(base_size = 13) ``` Odhadnutý `r render_glossary_term("sklon", display = "sklon")` je `r format_cz(odhad_sklonu)` minuty čekání na minutu erupce: dvě erupce, jejichž délky se liší o jednu minutu, se podle modelu v průměru liší čekáním na další erupci přibližně o `r format_cz(odhad_sklonu)` minuty. Stejná data a stejný vzorec modelu dávají vždy stejné `r render_glossary_term("odhad", display = "odhady")`; funkce `lm()` nic náhodně nelosuje. Odhad se může změnit jen tehdy, když se změní pozorované erupce. ### Sedí přímka k datům? Než budeme odhadu věřit, zkontrolujeme jako v předchozí lekci graf `r render_glossary_term("residuum", display = "residuí")` proti `r render_glossary_term("odhadnuta-hodnota", display = "odhadnutým hodnotám")`: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: graf-rezidui-gejzir #| echo: true #| fig-width: 8 #| fig-height: 5 #| fig-alt: "Odhadnuté čekání proti residuu každé erupce. Body tvoří dvě skupiny kolem nuly; několik erupcí střední délky leží výrazně nad nulovou čarou." plot( # Každý bod spojuje odhadnutou hodnotu s residuem stejné erupce x = fitted(object = mod_gejzir), y = resid(object = mod_gejzir), xlab = "Odhadnuté čekání (min)", ylab = "Residuum (min)", pch = 16, col = "#8A8A8A", cex = 1.3, cex.axis = 1.4, cex.lab = 1.4 ) abline( # Body hledáme rozptýlené kolem nulového residua h = 0, lty = 2, col = "#5D2890" ) ``` ```{=typst} ] ``` Residua tvoří dvě skupiny podle krátkých a dlouhých erupcí a v obou leží kolem nuly bez zřetelného vějíře, tedy bez rozšiřování rozptylu residuí (`r render_glossary_term("heteroskedasticita", display = "heteroskedasticity")`). Několik erupcí střední délky ale leží výrazně nad nulou: po nich se čekalo déle, než přímka odhaduje. Přímka tedy vztah dobře shrnuje, ale ne dokonale. Tuto výhradu si ponecháme, protože i `r render_glossary_term("standardni-chyba", display = "standardní chyba")` a `r render_glossary_term("interval-spolehlivosti", display = "interval spolehlivosti")`, ke kterým lekce směřuje, platí jen za předpokladu, že model data přiměřeně popisuje. Data jsou navíc skutečná, ale dobrovolníci nezaznamenali náhodně všechny erupce: pozorovali v dostupných termínech a blízké erupce mohou být časově propojené. Při interpretaci proto oddělujeme **vztah v zaznamenaných datech** od širšího procesu, který chceme popsat. ## Pepa: jeden měsíc, jeden `r render_glossary_term("odhad", display = "odhad")` Představme si výzkumníka **Pepu**, který gejzír pozoroval pouze v červnu. Ze všech erupcí ponechá jen ty červnové a odhadne stejný model: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: fitovat-pepuv-model #| echo: true # Ponecháme jen řádky, jejichž datum patří do června data_pepa <- data_gejzir[ format(x = data_gejzir$date, format = "%m") == "06", ] # Odhadneme stejný vztah pouze z Pepových dat mod_pepa <- lm( formula = cekani_min ~ delka_erupce_min, data = data_pepa ) # Zobrazíme odhadnutý intercept a sklon coef(object = mod_pepa) ``` ```{=typst} ] ``` ```{r} #| label: pepova-data-a-primka #| echo: false #| message: false #| warning: false #| fig-cap: "Pepova červnová pozorování a přímka odhadnutá pouze z nich." #| fig-alt: "Bodový graf červnových erupcí s fialovou rostoucí přímkou odhadnutou pouze z nich." ggplot2::ggplot( data = data_pepa, mapping = ggplot2::aes( x = delka_erupce_min, y = cekani_min ) ) + ggplot2::geom_point( color = biostat_cols[["grey_olive"]], alpha = 0.75, size = 2.5 ) + # No band: the confidence band is introduced later in this lesson. ggplot2::geom_abline( intercept = coef(object = mod_pepa)[["(Intercept)"]], slope = coef(object = mod_pepa)[["delka_erupce_min"]], color = biostat_cols[["indigo_velvet"]], linewidth = 1.1 ) + ggplot2::labs( x = "Délka erupce (min)", y = "Čekání na další erupci (min)" ) + theme_biostat(base_size = 13) ``` Pepa zaznamenal `r format_cz(nrow(data_pepa))` erupcí. Jeho `r render_glossary_term("odhad", display = "odhad")` `r render_glossary_term("sklon", display = "sklonu")` je `r format_cz(coef(object = mod_pepa)[["delka_erupce_min"]])` minuty čekání na minutu erupce: v jeho datech se dvě erupce, jejichž délky se liší o jednu minutu, v průměru liší čekáním přibližně o `r format_cz(coef(object = mod_pepa)[["delka_erupce_min"]])` minuty. Tento výsledek správně popisuje **Pepova data**. Neznamená však, že by každý další pozorovatel musel získat stejný sklon. ## Mařenka a Karel: jiné měsíce, jiné přímky ```{r} #| label: pripravit-mesicni-modely #| include: false vec_poradi_vyzkumniku <- c( "Pepa – červen", "Mařenka – červenec", "Karel – srpen" ) data_mesicni <- data_gejzir |> dplyr::mutate( mesic = format(x = date, format = "%m"), vyzkumnik = dplyr::case_when( mesic == "06" ~ "Pepa – červen", mesic == "07" ~ "Mařenka – červenec", mesic == "08" ~ "Karel – srpen" ), vyzkumnik = factor( x = vyzkumnik, levels = vec_poradi_vyzkumniku ) ) |> tidyr::drop_na(vyzkumnik) data_modely_mesicni <- data_mesicni |> dplyr::group_by(vyzkumnik) |> tidyr::nest() |> dplyr::ungroup() |> dplyr::mutate( model = purrr::map( .x = data, .f = ~ lm( formula = cekani_min ~ delka_erupce_min, data = .x ) ), pocet_erupci = purrr::map_int( .x = data, .f = ~ nrow(.x) ), intercept = purrr::map_dbl( .x = model, .f = ~ coef(.x)[["(Intercept)"]] ), sklon = purrr::map_dbl( .x = model, .f = ~ coef(.x)[["delka_erupce_min"]] ), popisek = paste0( as.character(vyzkumnik), "\nn = ", pocet_erupci, "; sklon = ", format_cz(sklon) ) ) |> dplyr::arrange(vyzkumnik) |> dplyr::mutate( popisek = factor( x = popisek, levels = popisek ) ) data_mesicni_graf <- data_mesicni |> dplyr::left_join( data_modely_mesicni |> dplyr::select(vyzkumnik, popisek), by = dplyr::join_by(vyzkumnik) ) ``` Dalšími pozorovateli jsou **Mařenka**, která sledovala gejzír v červenci, a **Karel**, který jej sledoval v srpnu. Oba použijí stejné `r render_glossary_term("promenna", display = "proměnné")` a stejný model jako Pepa. ::: {.box-question title="Co čekáte?"} Budou mít Pepa, Mařenka a Karel stejný `r render_glossary_term("sklon", display = "sklon")`? Pokud ne, čím se jejich výpočty liší? ::: ```{r} #| label: porovnat-mesicni-primky #| echo: false #| message: false #| warning: false #| fig-width: 10 #| fig-height: 4.5 #| fig-cap: "Tři lineární modely odhadnuté stejným postupem na erupcích ze tří měsíců. V záhlaví každého panelu je počet erupcí a odhad sklonu." #| fig-alt: "Tři bodové grafy vedle sebe pro červen, červenec a srpen. Každý má vlastní fialovou přímku s mírně odlišným sklonem." ggplot2::ggplot( data = data_mesicni_graf, mapping = ggplot2::aes( x = delka_erupce_min, y = cekani_min ) ) + ggplot2::geom_point( color = biostat_cols[["grey_olive"]], alpha = 0.7, size = 1.9 ) + # No band: the confidence band is introduced later in this lesson. ggplot2::geom_smooth( method = "lm", formula = y ~ x, se = FALSE, color = biostat_cols[["indigo_velvet"]], linewidth = 1 ) + ggplot2::facet_wrap( facets = ggplot2::vars(popisek), ncol = 3 ) + ggplot2::labs( x = "Délka erupce (min)", y = "Čekání na další erupci (min)" ) + theme_biostat(base_size = 15) ``` ```{r} #| label: tabulka-mesicnich-koeficientu #| echo: false #| message: false #| warning: false #| tbl-cap: "Koeficienty tří modelů odhadnutých odděleně po měsících." data_modely_mesicni |> dplyr::transmute( `Výzkumník a měsíc` = as.character(vyzkumnik), `Počet erupcí` = format_cz(pocet_erupci), `Intercept (min)` = format_cz(intercept), `Sklon (min/min)` = format_cz(sklon) ) |> knitr::kable(align = "lrrr") ``` Pepa, Mařenka i Karel použili **stejný postup**: stejné proměnné, stejný model a stejnou funkci `lm()`. Změnila se pouze pozorovaná data — a s nimi přímka i sklon. Rozdíly mezi měsíci mohou zahrnovat náhodnou proměnlivost i skutečné změny chování gejzíru; tři měsíce proto nejsou dokonale kontrolovaný experiment. Jsou ale prvním konkrétním důkazem, že **`r render_glossary_term("odhad", display = "odhad")` není neměnné číslo oddělené od dat, ze kterých vznikl**. Jednomu souboru `r render_glossary_term("pozorovani", display = "pozorování")`, ze kterého počítáme odhad, se odborně říká `r render_glossary_term("vyber", display = "výběr")` (anglicky *sample*). Česky se používá také **výběrový soubor** nebo méně formálně **vzorek**. Proměnlivost odhadu mezi možnými výběry nazýváme `r render_glossary_term("vyberova-variabilita", display = "výběrová variabilita")` (anglicky *sampling variability*): jiné pozorované erupce mohou vytvořit jinou přímku. ## Sto hlášení v aplikaci GeyserWatch ```{r} #| label: nastavit-hlaseni-geyserwatch #| include: false pocet_hlaseni <- 100L pocet_erupci_v_hlaseni <- 40L pocet_ukazkovych_hlaseni <- 15L vec_hladiny_pokryti <- c(0.95, 0.80, 0.50) ``` Tři měsíce nám daly pouze tři přímky. Abychom viděli, jak se `r render_glossary_term("odhad", display = "odhady")` chovají při mnoha opakováních, zavedeme fiktivní aplikaci **GeyserWatch**. Do ní pošle hlášení `r format_cz(pocet_hlaseni)` návštěvníků; každé hlášení obsahuje `r format_cz(pocet_erupci_v_hlaseni)` erupcí. Z každého hlášení odhadneme stejný model a uložíme jeho `r render_glossary_term("sklon", display = "sklon")`. Pomůže oddělit tři pojmy. Rozdíl mezi výběrem a `r render_glossary_term("populace", display = "cílovou populací")` jsme zmínili v první lekci; přibývá k nim odhad: - **Skutečnost** (širší proces, cílová populace) je to, jak se gejzír chová u všech erupcí. Skutečný sklon nikdy přímo nevidíme. - **`r render_glossary_term("vyber", display = "Výběr")`** jsou erupce, které někdo skutečně zaznamenal, například našich `r format_cz(nrow(data_gejzir))`. - **Odhad** je číslo spočítané z výběru, například sklon `r format_cz(odhad_sklonu)`. Ve výzkumu máme jeden výběr a jeden odhad, skutečnost neznáme. V `r render_glossary_term("simulace", display = "simulaci")` si naopak skutečnost zvolíme sami a vytvoříme z ní mnoho výběrů. Proto u simulace známe skutečný (referenční) sklon a můžeme sledovat, jak daleko od ní odhady z jednotlivých výběrů padají. Hlášení jsou **simulovaná**. Délky erupcí v každém hlášení vybereme z našich skutečných erupcí. Čekání k nim vytvoříme podle přímky odhadnuté ze všech erupcí a k přímce přidáme náhodné odchylky, které jsou typicky stejně velké jako `r render_glossary_term("residuum", display = "residua")` kolem ní. Každé hlášení tak přináší nová čekání, jako by šlo o nové erupce. Stejný zůstává model, přímka, ze které čekání vytváříme, i typická velikost odchylek. Sklon této přímky (`r format_cz(odhad_sklonu)` minuty na minutu) je proto v simulaci známý **referenční sklon**; skutečný sklon všech erupcí Old Faithful neznáme. ::: {.box-data collapse="true" title="Doplňující: proč hlášení simulujeme a co z toho plyne"} Sto skutečných nezávislých pozorovatelů nemáme. Simulace vychází z přímky odhadnuté ze všech `r format_cz(nrow(data_gejzir))` erupcí. **V této ukázce proto známe „správný“ sklon** — je to sklon `r format_cz(odhad_sklonu)` minuty na minutu, ze kterého jsme čekání vytvořili. Tato **referenční hodnota** platí jen pro simulaci; skutečný sklon všech erupcí Old Faithful neznáme. Mezi hlášeními se mění vybrané délky erupcí a náhodné odchylky čekání. Stejný zůstává model, referenční přímka i typická velikost odchylek. Proč nepoužíváme přímo výběry ze `r format_cz(nrow(data_gejzir))` skutečných erupcí? Všechna hlášení by pak sdílela tytéž erupce, byla by si podobnější než skutečně nové výběry a jejich sklony by se lišily méně, než odpovídá novým `r render_glossary_term("pozorovani", display = "pozorováním")`. ::: ```{r} #| label: vytvorit-hlaseni-geyserwatch #| include: false set.seed(900723) data_hlaseni <- simulate_reports( data_source = data_gejzir, mod_reference = mod_gejzir, n_reports = pocet_hlaseni, n_eruptions = pocet_erupci_v_hlaseni ) data_odhady_hlaseni <- summarise_reports( data_reports = data_hlaseni, vec_levels = vec_hladiny_pokryti ) data_ukazkova_hlaseni <- data_hlaseni |> dplyr::filter(hlaseni <= pocet_ukazkovych_hlaseni) |> dplyr::left_join( data_odhady_hlaseni |> dplyr::select(hlaseni, odhad_sklonu), by = dplyr::join_by(hlaseni) ) |> dplyr::mutate( popisek = paste0( "Hlášení ", formatC(x = hlaseni, width = 2, flag = "0"), "\nsklon = ", format_cz(odhad_sklonu) ) ) ``` Takto vypadá `r format_cz(pocet_ukazkovych_hlaseni)` ze `r format_cz(pocet_hlaseni)` hlášení: ```{r} #| label: patnact-hlaseni-a-modelu #| echo: false #| message: false #| warning: false #| fig-width: 11 #| fig-height: 7.5 #| fig-cap: !expr paste0("Prvních ", pocet_ukazkovych_hlaseni, " simulovaných hlášení GeyserWatch. Každý panel obsahuje ", pocet_erupci_v_hlaseni, " erupcí, vlastní přímku a vlastní sklon.") #| fig-alt: "Mřížka malých bodových grafů simulovaných hlášení. Každý obsahuje body jednoho hlášení a fialovou přímku; některé přímky jsou strmější, jiné plošší." ggplot2::ggplot( data = data_ukazkova_hlaseni, mapping = ggplot2::aes( x = delka_erupce_min, y = cekani_min ) ) + ggplot2::geom_point( color = biostat_cols[["grey_olive"]], alpha = 0.55, size = 0.9 ) + # No band: simulated report lines; the band is introduced later in this lesson. ggplot2::geom_smooth( method = "lm", formula = y ~ x, se = FALSE, color = biostat_cols[["indigo_velvet"]], linewidth = 0.7 ) + ggplot2::facet_wrap( facets = ggplot2::vars(popisek), ncol = 5 ) + ggplot2::labs( x = "Délka erupce (min)", y = "Čekání na další erupci (min)" ) + theme_biostat(base_size = 12) + ggplot2::theme( strip.text = ggplot2::element_text(size = 10) ) ``` *Které přímky jsou strmější a které plošší? Dokázali byste z jediného panelu poznat, jak moc by se mohl lišit panel další?* Každý panel opakuje stejný postup, ale body, přímky i sklony se mění. Kdybychom měli pouze jedno hlášení, neviděli bychom, jak široká je množina dalších rozumně možných odhadů. Všech `r format_cz(pocet_hlaseni)` hlášení nám dává `r format_cz(pocet_hlaseni)` sklonů a jejich `r render_glossary_term("histogram", display = "histogram")` je ukazuje společně: ```{r} #| label: histogram-sta-sklonu #| echo: false #| message: false #| warning: false #| fig-width: 8 #| fig-height: 4.8 #| fig-cap: !expr paste0("Všech ", pocet_hlaseni, " simulovaných hlášení GeyserWatch vytvořilo ", pocet_hlaseni, " odhadů sklonu.") #| fig-alt: "Histogram sta odhadů sklonu. Většina leží uprostřed, méně odhadů je na obou okrajích." graf_histogram_sklonu <- ggplot2::ggplot( data = data_odhady_hlaseni, mapping = ggplot2::aes(x = odhad_sklonu) ) + ggplot2::geom_histogram( bins = 18, fill = biostat_cols[["amethyst"]], color = biostat_cols[["white"]] ) + ggplot2::labs( x = "Odhad sklonu (min čekání / min erupce)", y = "Počet hlášení" ) + theme_biostat(base_size = 12) graf_histogram_sklonu ``` Většina sklonů leží ve společné oblasti, ale nejsou totožné. Odhady z různých výběrů tvoří vlastní `r render_glossary_term("rozdeleni", display = "rozdělení")`: některé hodnoty jsou časté, jiné vzácné. ## Co měří `r render_glossary_term("standardni-chyba", display = "standardní chyba")`? ```{r} #| label: spocitat-rozptyleni-sklonu #| include: false sd_sklonu_hlaseni <- sd(data_odhady_hlaseni$odhad_sklonu) typicka_se_hlaseni <- median(data_odhady_hlaseni$se_sklonu) vec_kvantily_se_hlaseni <- c(0.10, 0.90) vec_rozsah_se_hlaseni <- quantile( x = data_odhady_hlaseni$se_sklonu, probs = vec_kvantily_se_hlaseni, names = FALSE ) pocet_v_jedne_se <- sum( abs(data_odhady_hlaseni$odhad_sklonu - odhad_sklonu) <= typicka_se_hlaseni ) ``` V první lekci jsme odlišili dvě otázky. `r render_glossary_term("smerodatna-odchylka", display = "Směrodatná odchylka")` (SD) odpovídá na otázku **jak rozdílná jsou `r render_glossary_term("pozorovani", display = "pozorování")`?** Druhá otázka zní **jak přesný je náš `r render_glossary_term("odhad", display = "odhad")`?** `r render_glossary_term("histogram", display = "Histogram")` `r render_glossary_term("sklon", display = "sklonů")` nám na ni dává první odpověď: kdyby se odhady z různých `r render_glossary_term("vyber", display = "výběrů")` téměř nelišily, byl by histogram úzký; čím víc se liší, tím je širší. Šířku histogramu můžeme shrnout stejně jako jakékoli jiné hodnoty, směrodatnou odchylkou. Všech `r format_cz(pocet_hlaseni)` sklonů z GeyserWatch má směrodatnou odchylku **`r format_cz(sd_sklonu_hlaseni)` minuty čekání na minutu erupce**. Skutečný výzkumník ale tolik hlášení nemá — má jediné. Přesto mu R tuto velikost odhadne. Z jediného hlášení vypočítá u sklonu číslo, které ve výstupu `summary()` nazývá `Std. Error` (celý výstup si prohlédneme později). Porovnejme je se směrodatnou odchylkou všech sklonů: ```{r} #| label: tabulka-se-hlaseni #| echo: false #| message: false #| warning: false #| tbl-cap: "Rozptýlení sta sklonů a standardní chyba, kterou R vypočítá z jediného hlášení." tibble::tibble( `Co počítáme` = c( paste0( "Směrodatná odchylka ", pocet_hlaseni, " sklonů (potřebuje všech ", pocet_hlaseni, " hlášení)" ), "Std. Error z jednoho hlášení, typická hodnota (medián)", paste0( "Std. Error z jednoho hlášení, rozmezí u prostředních ", format_cz(100 * diff(vec_kvantily_se_hlaseni), digits = 0), " % hlášení" ) ), Hodnota = c( formatC(x = sd_sklonu_hlaseni, format = "f", digits = 3, decimal.mark = ","), formatC(x = typicka_se_hlaseni, format = "f", digits = 3, decimal.mark = ","), paste0( format_cz(vec_rozsah_se_hlaseni[[1]]), " až ", format_cz(vec_rozsah_se_hlaseni[[2]]) ) ) ) |> knitr::kable(align = "lr") ``` Obě čísla jsou si velmi blízká, i když nejsou stejná. `Std. Error` je `r render_glossary_term("standardni-chyba", display = "standardní chyba")` odhadu sklonu (anglicky *standard error*, zkratka **SE**). **Standardní chyba odhaduje, jak moc by se odhady sklonu typicky lišily mezi výběry stejné velikosti — a odhaduje to z jediného výběru.** Má stejné jednotky jako sklon (minuty čekání na minutu délky erupce). Hodnota z jednoho hlášení se od hlášení k hlášení trochu liší, protože i ona je odhadem. ```{r} #| label: histogram-a-jedna-se #| echo: false #| message: false #| warning: false #| fig-width: 8 #| fig-height: 5.2 #| fig-cap: "Histogram simulovaných sklonů s referenčním sklonem simulace (tmavá čára) a rozsahem plus minus jedna typická standardní chyba (oranžová úsečka)." #| fig-alt: "Histogram sta sklonů. Svislá tmavá čára označuje referenční sklon. Nad histogramem je oranžová úsečka sahající jednu standardní chybu na každou stranu; pokrývá střední část histogramu." max_cetnost_sklonu <- max( ggplot2::ggplot_build(graf_histogram_sklonu)$data[[1]]$count ) data_rozsah_se <- tibble::tibble( dolni = odhad_sklonu - typicka_se_hlaseni, horni = odhad_sklonu + typicka_se_hlaseni, y = max_cetnost_sklonu + 1.5 ) graf_histogram_sklonu + ggplot2::geom_vline( xintercept = odhad_sklonu, color = biostat_cols[["graphite"]], linetype = "dashed", linewidth = 0.9 ) + ggplot2::geom_errorbarh( data = data_rozsah_se, mapping = ggplot2::aes( xmin = dolni, xmax = horni, y = y ), inherit.aes = FALSE, height = 0.7, color = biostat_cols[["orange"]], linewidth = 1.6 ) + ggplot2::annotate( geom = "text", x = odhad_sklonu + typicka_se_hlaseni, y = max_cetnost_sklonu + 1.5, label = " ± 1 SE", hjust = 0, color = biostat_cols[["graphite"]], size = 4 ) + ggplot2::scale_y_continuous( expand = ggplot2::expansion(mult = c(0, 0.18)) ) ``` Do vzdálenosti jedné typické standardní chyby od referenčního sklonu padlo `r format_cz(pocet_v_jedne_se)` ze `r format_cz(pocet_hlaseni)` sklonů, tedy asi `r format_cz(100 * pocet_v_jedne_se / pocet_hlaseni, digits = 0)` %. Rozsah **± jedna standardní chyba** ale ještě není `r render_glossary_term("interval-spolehlivosti", display = "interval spolehlivosti")`. Je to měřítko přesnosti, ze kterého se interval vytváří. *Jak ale může R poznat šířku histogramu z jediného hlášení, když ostatní hlášení nikdy neviděl?* ## Co dělá sklon stabilnějším? ```{r} #| label: pripravit-dvojice-stability #| include: false pocet_erupci_v_malem_hlaseni <- 10L min_delka_dlouhe_erupce <- 3 list_dvojice_stability <- simulate_stability_pairs( data_source = data_gejzir, mod_reference = mod_gejzir, n_reports = pocet_hlaseni, n_eruptions = pocet_erupci_v_hlaseni, n_eruptions_small = pocet_erupci_v_malem_hlaseni, min_length_long = min_delka_dlouhe_erupce ) # Each pair changes one SE ingredient; the GeyserWatch setting appears in all # three pairs. The steadier option is B, A, B so the answer is not one side. data_panely_stability <- tibble::tribble( ~dvojice, ~moznost, ~nastaveni, "1", "A", "scatter_double", "1", "B", "geyserwatch", "2", "A", "geyserwatch", "2", "B", "lengths_long", "3", "A", "eruptions_few", "3", "B", "geyserwatch" ) |> dplyr::mutate( popis = dplyr::case_when( dvojice == "1" & moznost == "A" ~ "body daleko od přímky", dvojice == "2" & moznost == "B" ~ paste0( "jen erupce od ", format_cz(min_delka_dlouhe_erupce), " min" ), dvojice == "3" & moznost == "A" ~ paste0( pocet_erupci_v_malem_hlaseni, " erupcí v hlášení" ), dvojice == "1" ~ "body blízko přímky", dvojice == "2" ~ "délky z celého rozsahu", dvojice == "3" ~ paste0(pocet_erupci_v_hlaseni, " erupcí v hlášení") ) ) data_sd_stability <- list_dvojice_stability$data_lines |> dplyr::group_by(nastaveni) |> dplyr::summarise( sd_sklonu = sd(odhad_sklonu) ) |> dplyr::ungroup() data_popisky_stability <- data_panely_stability |> dplyr::left_join( data_sd_stability, by = dplyr::join_by(nastaveni) ) |> dplyr::mutate( panel = paste0( dvojice, moznost, ": ", popis, "\nSD sklonů = ", format_cz(sd_sklonu) ), panel = factor(x = panel, levels = panel) ) data_primky_stability <- data_popisky_stability |> dplyr::select(nastaveni, panel, dvojice) |> dplyr::left_join( list_dvojice_stability$data_lines, by = dplyr::join_by(nastaveni), relationship = "many-to-many" ) |> dplyr::mutate( y = intercept + odhad_sklonu * x_od, yend = intercept + odhad_sklonu * x_do ) data_body_stability <- data_popisky_stability |> dplyr::select(nastaveni, panel, dvojice) |> dplyr::left_join( list_dvojice_stability$data_points, by = dplyr::join_by(nastaveni), relationship = "many-to-many" ) vec_sd_stability <- data_sd_stability$sd_sklonu |> rlang::set_names(data_sd_stability$nastaveni) ``` Šířka `r render_glossary_term("histogram", display = "histogramu")` `r render_glossary_term("sklon", display = "sklonů")` ale nemusí být vždy stejná. Představme si tři srovnání dvou variant aplikace GeyserWatch. V každém srovnání se obě varianty liší jen jednou věcí: ```{=typst} #block(breakable: false)[ ``` ::: {.box-question title="Co čekáte?"} Ve které variantě, A, nebo B, se budou sklony jednotlivých hlášení lišit méně? | Srovnání | Co se mění | Varianta A | Varianta B | |:--|:--|:--|:--| | 1 | body | daleko od přímky | blízko přímky | | 2 | délky erupcí | z celého rozsahu | jen od `r format_cz(min_delka_dlouhe_erupce)` min | | 3 | počet erupcí | `r pocet_erupci_v_malem_hlaseni` | `r pocet_erupci_v_hlaseni` | ::: ```{=typst} ] ``` ::: {.box-result collapse="true" title="Odpověď: 1B, 2A, 3B"} Stabilnější sklony dávají **body blízko přímky (1B)**, **délky erupcí z celého rozsahu (2A)** a **více erupcí v hlášení (3B)**. Následující graf to ukazuje na simulovaných hlášeních. ::: ```{r} #| label: graf-dvojice-stability #| echo: false #| message: false #| warning: false #| fig-width: 9 #| fig-height: 5.6 #| fig-cap: !expr paste0("Tři srovnání simulovaných variant aplikace, každá varianta se ", pocet_hlaseni, " hlášeními. Fialové přímky patří k jednotlivým hlášením, šedé body ukazují erupce jednoho z nich. V záhlaví panelu je směrodatná odchylka (SD) sklonů všech hlášení, vpravo dole histogram těchto sklonů; jeho vodorovná osa je ve všech panelech stejná, užší histogram tedy znamená menší SD.") #| fig-alt: "Šest panelů ve třech řadách po dvou. V každém panelu je vějíř stovky fialových přímek. Užší vějíř mají panely s body blízko přímky, s délkami erupcí z celého rozsahu a se čtyřiceti erupcemi v hlášení; širší vějíř mají panely s body daleko od přímky, s jen dlouhými erupcemi a s deseti erupcemi v hlášení. Vpravo dole je v každém panelu malý histogram sklonů; u užších vějířů je užší." plot_stability_pairs( data_points = data_body_stability, data_lines = data_primky_stability, data_panels = data_popisky_stability, max_pair = 3, xlim = range(data_gejzir$delka_erupce_min), ylim = range(data_gejzir$cekani_min) + c(-15, 15), base_size = 13, n_col = 2, show_sd_label = FALSE ) + ggplot2::labs(y = "Čekání na další erupci (min)") + ggplot2::theme( strip.text = ggplot2::element_text(size = 11) ) ``` Vždy jsme změnili jen jednu věc a vějíř přímek se rozšířil nebo zúžil: 1. **Body daleko od přímky** dávají širší vějíř. Když se čekání kolem přímky hodně rozptylují, může každé hlášení náhodou zachytit body, které přímku víc nakloní. Při dvojnásobném rozptýlení bodů vychází i `r render_glossary_term("smerodatna-odchylka", display = "směrodatná odchylka")` sklonů zhruba dvojnásobná (`r format_cz(vec_sd_stability[["scatter_double"]])` místo `r format_cz(vec_sd_stability[["geyserwatch"]])`). 2. **Délky erupcí z celého rozsahu** dávají užší vějíř. Přímku si můžeme představit jako houpačku: body daleko od sebe po vodorovné ose ji drží pevně, body natlačené blízko sebe ji nechají snadno naklonit. Když hlášení obsahují jen erupce od `r format_cz(min_delka_dlouhe_erupce)` minut, je směrodatná odchylka sklonů `r format_cz(vec_sd_stability[["lengths_long"]])` místo `r format_cz(vec_sd_stability[["geyserwatch"]])`. 3. **Více erupcí v hlášení** dává užší vějíř. Hlášení s `r pocet_erupci_v_malem_hlaseni` erupcemi má směrodatnou odchylku sklonů `r format_cz(vec_sd_stability[["eruptions_few"]])`, hlášení se `r pocet_erupci_v_hlaseni` erupcemi `r format_cz(vec_sd_stability[["geyserwatch"]])`. Tyto tři vlastnosti výzkumník zná i z jediného souboru dat: vidí, jak daleko leží body od přímky, jak široce jsou rozprostřené délky erupcí a kolik erupcí zaznamenal. Právě z nich R `r render_glossary_term("standardni-chyba", display = "standardní chybu")` počítá. ## Standardní chyba z jediného výběru ```{r} #| label: pripravit-vypocet-standardni-chyby #| include: false # The displayed arithmetic uses rounded values at every step so readers can # repeat it on a calculator; R's unrounded result is reported separately. rezidualni_sd_zobrazena <- round(sigma(object = mod_gejzir), 2) rss_zobrazeny <- round(sum(resid(object = mod_gejzir)^2), 1) podil_rss_zobrazeny <- round(rss_zobrazeny / (nrow(data_gejzir) - 2), 2) stopifnot( abs(round(sqrt(podil_rss_zobrazeny), 2) - rezidualni_sd_zobrazena) < 0.005 ) sd_delky_zobrazena <- round(sd(data_gejzir$delka_erupce_min), 3) odmocnina_n_zobrazena <- round(sqrt(nrow(data_gejzir) - 1), 2) jmenovatel_zobrazeny <- round(sd_delky_zobrazena * odmocnina_n_zobrazena, 2) standardni_chyba_zobrazena <- round(rezidualni_sd_zobrazena / jmenovatel_zobrazeny, 2) standardni_chyba_summary <- summary(object = mod_gejzir)$coefficients[ "delka_erupce_min", "Std. Error" ] stopifnot( abs(standardni_chyba_zobrazena - standardni_chyba_summary) < 0.01 ) # The same formula with the number of eruptions of one GeyserWatch report # and of one small report from the third pair standardni_chyba_hlaseni_vzorec <- sigma(object = mod_gejzir) / (sd(data_gejzir$delka_erupce_min) * sqrt(pocet_erupci_v_hlaseni - 1)) standardni_chyba_male_hlaseni_vzorec <- sigma(object = mod_gejzir) / (sd(data_gejzir$delka_erupce_min) * sqrt(pocet_erupci_v_malem_hlaseni - 1)) ``` `r render_glossary_term("standardni-chyba", display = "Standardní chyba")` `r render_glossary_term("sklon", display = "sklonu")` spojuje právě tři vlastnosti dat z předchozí kapitoly; budeme jim říkat ingredience standardní chyby: 1. jak daleko leží body od přímky — typická velikost `r render_glossary_term("residuum", display = "residua")`; 2. jak široce jsou rozprostřené délky erupcí — jejich `r render_glossary_term("smerodatna-odchylka", display = "směrodatná odchylka")`; 3. kolik erupcí máme. První vlastnost stojí v čitateli: větší residua standardní chybu zvětšují. Zbylé dvě stojí ve jmenovateli: širší rozprostření délek a více erupcí ji zmenšují. $$ \text{standardní chyba sklonu} = \frac{\text{typická velikost residua}} {\text{směrodatná odchylka délek erupcí} \times \sqrt{\text{počet erupcí} - 1}}. $$ Odečtení jedničky pod odmocninou souvisí s tím, jak se počítá směrodatná odchylka (vysvětlíme níže u obvyklého zápisu vzorce). Při desítkách erupcí je to prakticky odmocnina z jejich počtu. ### Typická velikost residua Residua známe z předchozí lekce: jsou to svislé odchylky erupcí od přímky. Jejich typickou velikost shrneme podobně jako směrodatnou odchylku. Residua umocníme a sečteme. Tento `r render_glossary_term("soucet-ctvercu-residui", display = "součet čtverců residuí")` (anglicky *residual sum of squares*, **RSS**) známe z minulé lekce, kde jsme ho značili SSE: funkce `lm()` vybírá přímku právě tak, aby byl co nejmenší. Z něj odhadneme typickou velikost residua. Součet vydělíme počtem erupcí zmenšeným o 2 a odmocníme, abychom se vrátili k minutám čekání: $$ \text{residuální směrodatná odchylka} = \sqrt{\frac{\text{součet čtverců residuí}}{\text{počet erupcí} - 2}}. $$ V R to spočítáme takto: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: spocitat-rezidualni-sd #| echo: true # Residua: svislé odchylky erupcí od přímky (v minutách) vec_residua <- resid(object = mod_gejzir) # Počet erupcí v modelu pocet_erupci <- nrow(x = data_gejzir) # Součet čtverců residuí vydělíme počtem erupcí minus 2 a odmocníme rezidualni_sd <- sqrt(x = sum(vec_residua^2) / (pocet_erupci - 2)) rezidualni_sd ``` ```{=typst} ] ``` Výsledek se nazývá `r render_glossary_term("rezidualni-smerodatna-odchylka", display = "residuální směrodatná odchylka")` (ve výstupu R *Residual standard error*): erupce leží typicky asi `r format_cz(rezidualni_sd_zobrazena)` minuty nad přímkou nebo pod ní. V R ji vrátí funkce `sigma()`. Výraz počet erupcí minus 2 se nazývá `r render_glossary_term("stupne-volnosti", display = "stupně volnosti")` (anglicky *degrees of freedom*); proč se odečítá právě 2, vysvětlíme v následující lekci. Po dosazení zaokrouhlených hodnot: $$ \text{residuální směrodatná odchylka} = \sqrt{\frac{`r format_cz_math(rss_zobrazeny, digits = 1, nsmall = 1)`}{`r nrow(data_gejzir) - 2`}} = \sqrt{`r format_cz_math(podil_rss_zobrazeny, nsmall = 2)`} = `r format_cz_math(rezidualni_sd_zobrazena, nsmall = 2)`\ \text{min}. $$ Obecně zapisujeme: $$ \hat{\sigma} = \sqrt{\frac{\operatorname{RSS}}{n - 2}}. $$ Zde $\hat{\sigma}$ označuje residuální směrodatnou odchylku, $\operatorname{RSS}$ součet čtverců residuí a $n$ počet erupcí. ### Rozptýlení délek a počet erupcí ```{=typst} #block(breakable: false)[ ``` Zbylé dvě ingredience spočítáme funkcemi, které známe z první lekce: ```{r} #| label: spocitat-standardni-chybu-rucne #| echo: true # Jak široce jsou rozprostřené délky erupcí sd_delky <- sd(x = data_gejzir$delka_erupce_min) sd_delky # Standardní chyba sklonu ze tří ingrediencí standardni_chyba_rucne <- rezidualni_sd / (sd_delky * sqrt(x = pocet_erupci - 1)) standardni_chyba_rucne ``` ```{=typst} ] ``` Po dosazení zaokrouhlených hodnot dostaneme: $$ \begin{aligned} \text{standardní chyba sklonu} &= \frac{`r format_cz_math(rezidualni_sd_zobrazena, nsmall = 2)`} {`r format_cz_math(sd_delky_zobrazena, digits = 3, nsmall = 3)` \times \sqrt{`r nrow(data_gejzir)` - 1}} \\ &= \frac{`r format_cz_math(rezidualni_sd_zobrazena, nsmall = 2)`} {`r format_cz_math(sd_delky_zobrazena, digits = 3, nsmall = 3)` \times `r format_cz_math(odmocnina_n_zobrazena, nsmall = 2)`} = \frac{`r format_cz_math(rezidualni_sd_zobrazena, nsmall = 2)`}{`r format_cz_math(jmenovatel_zobrazeny, nsmall = 2)`} = `r format_cz_math(standardni_chyba_zobrazena, nsmall = 2)`. \end{aligned} $$ Standardní chyba sklonu ze všech `r format_cz(nrow(data_gejzir))` erupcí je `r format_cz(standardni_chyba_rucne)` minuty čekání na minutu erupce. Stejné číslo vypisuje i funkce `summary()`. ### Souhlasí vzorec se simulací? Ve vzorci můžeme dosadit i počet erupcí jednoho hlášení GeyserWatch. Kdyby hlášení mělo stejně rozptýlená residua i délky erupcí jako naše data a obsahovalo `r pocet_erupci_v_hlaseni` erupcí, vyšla by standardní chyba `r format_cz(standardni_chyba_hlaseni_vzorec)`. Směrodatná odchylka sta simulovaných sklonů byla `r format_cz(sd_sklonu_hlaseni)`. Vzorec tedy z jediného výběru odhadne šířku `r render_glossary_term("histogram", display = "histogramu")`, který bychom jinak museli vytvořit z mnoha hlášení. Vzorec také vysvětluje, co jsme viděli ve třech srovnáních. Dvojnásobná typická velikost residua zdvojnásobí standardní chybu. Délky erupcí rozprostřené jen v úzkém rozsahu mají malou směrodatnou odchylku, a standardní chyba proto roste. Počet erupcí stojí ve jmenovateli pod odmocninou: **standardní chyba klesá s odmocninou z počtu `r render_glossary_term("pozorovani", display = "pozorování")`**, ne s počtem samotným. Čtyřnásobek erupcí ji proto zmenší jen na polovinu, ne na čtvrtinu. Ve třetím srovnání vyšel rozdíl o něco větší: směrodatná odchylka sklonů `r format_cz(vec_sd_stability[["eruptions_few"]])` u deseti erupcí proti `r format_cz(vec_sd_stability[["geyserwatch"]])` u čtyřiceti. Vzorec se stejným rozprostřením délek dává pro deset erupcí `r format_cz(standardni_chyba_male_hlaseni_vzorec)`. Hlášení s deseti erupcemi ale mívají rozprostření délek velmi různé: někdy všechny délky padnou do jednoho shluku a přímka se pak opírá o úzký rozsah. Vzorec s rozprostřením délek z našich dat proto tento rozdíl trochu podcení; standardní chyba každého jednotlivého hlášení ale vychází z jeho vlastních délek. ### Obvyklý zápis vzorce ```{=typst} #block(breakable: false)[ ``` V učebnicích se tentýž vzorec obvykle zapisuje takto: $$ \operatorname{SE}(\hat{\beta}_1) = \frac{\hat{\sigma}}{\sqrt{S_{xx}}}, \qquad S_{xx} = \sum_{i=1}^{n}(x_i - \bar{x})^2 = (n - 1)\, s_x^2. $$ ```{=typst} ] ``` Zde $\hat{\beta}_1$ označuje odhad sklonu, $\operatorname{SE}(\hat{\beta}_1)$ jeho standardní chybu, $\hat{\sigma}$ residuální směrodatnou odchylku, $n$ počet erupcí, $x_i$ délku $i$-té erupce, $\bar{x}$ průměrnou délku a $s_x$ směrodatnou odchylku délek. $S_{xx}$ je součet čtverců odchylek délek erupcí od jejich průměru. Protože se směrodatná odchylka počítá z tohoto součtu vyděleného $n-1$, platí $S_{xx} = (n-1)\,s_x^2$ a oba zápisy dávají stejné číslo. ## Kde ve výstupu najdeme velikost a přesnost? ```{r} #| label: pripravit-souhrn-modelu #| include: false table_souhrn_modelu <- summary(object = mod_gejzir)$coefficients standardni_chyba_sklonu <- table_souhrn_modelu[ "delka_erupce_min", "Std. Error" ] ``` Funkce `summary()` vypíše pro model ze všech `r format_cz(nrow(data_gejzir))` erupcí `r render_glossary_term("odhad", display = "odhady")` i jejich `r render_glossary_term("standardni-chyba", display = "standardní chyby")` najednou. Celý výstup vypadá tak, jak jej R skutečně vypisuje: ```{r} #| label: precist-odhad-a-standardni-chybu #| echo: true # Souhrn obsahuje odhad sklonu i jeho standardní chybu summary(object = mod_gejzir) ``` ::: {.box-question title="Najděte ve výstupu dvě čísla"} 1. Které číslo popisuje **velikost vztahu** mezi délkou erupce a čekáním? 2. Které číslo popisuje **přesnost tohoto odhadu**? ::: ```{=typst} #block(breakable: false)[ ``` ::: {.box-result collapse="true" title="Odpověď: řádek delka_erupce_min, sloupce Estimate a Std. Error"} V tabulce `Coefficients` najdeme řádek `delka_erupce_min`. Sloupec `Estimate` obsahuje odhad `r render_glossary_term("sklon", display = "sklonu")` `r format_cz(odhad_sklonu)` — velikost vztahu. Sloupec `Std. Error` obsahuje jeho standardní chybu `r format_cz(standardni_chyba_sklonu)` — přesnost odhadu. Obě čísla mají jednotku minuty čekání na minutu erupce. ::: ```{=typst} ] ``` Výstup můžeme číst shora: 1. `Call` připomíná vzorec modelu `cekani_min ~ delka_erupce_min`: čekání odhadujeme z délky předchozí erupce. 2. `Residuals` shrnují `r render_glossary_term("residuum", display = "residua")`, tedy svislé odchylky bodů od přímky z předchozí lekce. Pět čísel (minimum, dolní kvartil, medián, horní kvartil, maximum) čteme stejně jako `summary()` dat v první lekci: residua mají být rozložená kolem nuly zhruba souměrně. 3. V tabulce `Coefficients` obsahuje sloupec `Estimate` stejné odhady, které vrací `coef(mod_gejzir)`. Sloupec `Std. Error` vedle nich ukazuje jejich standardní chyby. 4. Řádek `Residual standard error` je `r render_glossary_term("rezidualni-smerodatna-odchylka", display = "residuální směrodatná odchylka")` `r format_cz(rezidualni_sd_zobrazena)` z předchozí kapitoly, tedy první ingredience standardní chyby. Za ní následuje počet `r render_glossary_term("stupne-volnosti", display = "stupňů volnosti")`. Přestože se tento řádek v angličtině jmenuje *standard error*, není to standardní chyba odhadu sklonu, ale typická velikost residuí. Sloupce `t value` a `Pr(>|t|)` a poslední řádky výstupu zatím přeskočíme. Souvisí s posuzováním konkrétních tvrzení o `r render_glossary_term("parametry-modelu", display = "parametrech modelu")`, kterému se budeme věnovat v následující lekci. `Std. Error` sklonu je přesně číslo, které jsme spočítali ze tří ingrediencí: `r format_cz(standardni_chyba_sklonu)`. Je menší než typická standardní chyba jednoho hlášení (`r format_cz(typicka_se_hlaseni)`), protože model ze všech erupcí vychází z `r format_cz(nrow(data_gejzir))` `r render_glossary_term("pozorovani", display = "pozorování")`, hlášení jen ze `r format_cz(pocet_erupci_v_hlaseni)`. **Velikost vztahu a přesnost odhadu jsou dvě různé informace** — velký sklon může být odhadnutý nepřesně a malý sklon velmi přesně. ::: {.box-result title="Vzor zápisu: odhad a standardní chyba"} „Dvě erupce, jejichž délky se lišily o jednu minutu, se podle `r render_glossary_term("linearni-model", display = "lineárního modelu")` v průměru lišily čekáním na další erupci o `r format_cz(odhad_sklonu)` minuty (odhad sklonu = `r format_cz(odhad_sklonu)`, SE = `r format_cz(standardni_chyba_sklonu)` minuty čekání na minutu erupce).“ Věta nejprve interpretuje vztah v běžném jazyce a až za ní uvádí kompaktní statistický zápis. Samotná dvojice čísel bez proměnných a jednotek by nestačila. ::: ::: {.box-data collapse="true" title="Doplňující: standardní chyba průměru"} ```{r} #| label: pripravit-standardni-chybu-prumeru #| include: false prumer_cekani <- mean(data_gejzir$cekani_min) sd_cekani <- sd(data_gejzir$cekani_min) standardni_chyba_prumeru <- sd_cekani / sqrt(nrow(data_gejzir)) ``` Standardní chybu nemá jen sklon. Má ji každý odhad, například i `r render_glossary_term("prumer", display = "průměr")`. Průměrné čekání na další erupci je v našich datech `r format_cz(prumer_cekani)` minuty. Kdyby někdo zaznamenal jiných `r format_cz(nrow(data_gejzir))` erupcí, vyšel by mu trochu jiný průměr. Jak moc by se takové průměry typicky lišily, odhadne standardní chyba průměru: $$ \text{standardní chyba průměru} = \frac{\text{směrodatná odchylka čekání}}{\sqrt{\text{počet erupcí}}}. $$ Po dosazení: $$ \text{standardní chyba průměru} = \frac{`r format_cz_math(sd_cekani, nsmall = 2)`}{\sqrt{`r nrow(data_gejzir)`}} = `r format_cz_math(standardni_chyba_prumeru, nsmall = 2)`\ \text{min}. $$ Obecně zapisujeme $\operatorname{SE}(\bar{y}) = s / \sqrt{n}$, kde $\bar{y}$ je průměr, $s$ směrodatná odchylka hodnot a $n$ jejich počet. Logika je stejná jako u sklonu: větší rozptýlení jednotlivých čekání standardní chybu zvětšuje a více erupcí ji zmenšuje, opět s odmocninou z jejich počtu. Chybí jen rozprostření délek erupcí, protože průměr žádnou přímku nenaklání. Průměr je ve skutečnosti také lineární model, jen bez `r render_glossary_term("prediktor", display = "prediktoru")`, tzv. `r render_glossary_term("nulovy-model", display = "nulový model")`. Ve vzorci modelu jej zapíšeme jako `~ 1`; jediný odhadnutý koeficient, `(Intercept)`, je pak průměrné čekání: ```{r} #| label: standardni-chyba-prumeru-model #| echo: true # Standardní chyba průměru podle vzorce sd(x = data_gejzir$cekani_min) / sqrt(x = nrow(x = data_gejzir)) # Stejný průměr jako model bez prediktoru mod_prumer <- lm( formula = cekani_min ~ 1, data = data_gejzir ) # Odhad a standardní chyba z tabulky Coefficients summary(object = mod_prumer)$coefficients[, c("Estimate", "Std. Error")] ``` Model bez prediktoru vrací stejný průměr i stejnou standardní chybu jako vzorec. ::: ::: {.box-data collapse="true" title="Doplňující: výpočet standardní chyby sklonu krok za krokem"} ```{r} #| label: pripravit-maly-vypocet-standardni-chyby #| include: false data_maly_se <- data_gejzir |> dplyr::slice_head(n = 5) |> dplyr::select( delka_erupce_min, cekani_min ) mod_maly_se <- lm( formula = cekani_min ~ delka_erupce_min, data = data_maly_se ) pocet_maly_se <- nrow(data_maly_se) standardni_chyba_maly_se <- summary(mod_maly_se)$coefficients[ "delka_erupce_min", "Std. Error" ] # The worked example uses the values displayed in the table (two decimals) # at every step, so readers can repeat the arithmetic exactly. R's own # unrounded result is reported separately below. prumer_delky_maly_se <- round(mean(data_maly_se$delka_erupce_min), 2) data_maly_se <- data_maly_se |> dplyr::mutate( erupce = dplyr::row_number(), delka = round(delka_erupce_min, 2), odhadnute_cekani = round(fitted(mod_maly_se), 2), residuum = round(residuals(mod_maly_se), 2), residuum_ctverec = round(residuum^2, 2), odchylka_delky = delka - prumer_delky_maly_se, odchylka_delky_ctverec = round(odchylka_delky^2, 2) ) soucet_ctvercu_residui_maly_se <- sum(data_maly_se$residuum_ctverec) podil_rss_maly_se <- round(soucet_ctvercu_residui_maly_se / (pocet_maly_se - 2), 2) rezidualni_smerodatna_odchylka_maly_se <- round(sqrt(podil_rss_maly_se), 2) soucet_ctvercu_delky_maly_se <- sum(data_maly_se$odchylka_delky_ctverec) odmocnina_sxx_maly_se <- round(sqrt(soucet_ctvercu_delky_maly_se), 2) standardni_chyba_rucne_maly_se <- round( rezidualni_smerodatna_odchylka_maly_se / odmocnina_sxx_maly_se, 2 ) # Hand rounding may differ from R only in the last displayed digit. stopifnot( abs(standardni_chyba_rucne_maly_se - standardni_chyba_maly_se) < 0.05 ) text_residua_mocniny <- paste0( "(", format_cz_math(data_maly_se$residuum, nsmall = 2), ")^2" ) text_residua_ctverce <- format_cz_math(data_maly_se$residuum_ctverec, nsmall = 2) text_odchylky_mocniny <- paste0( "(", format_cz_math(data_maly_se$odchylka_delky, nsmall = 2), ")^2" ) text_odchylky_ctverce <- format_cz_math(data_maly_se$odchylka_delky_ctverec, nsmall = 2) ``` Vzorec z kapitoly „Standardní chyba z jediného výběru“ můžeme zpomalit tak, aby každý krok šel spočítat na kalkulačce. Použijeme malý soubor prvních `r format_cz(pocet_maly_se)` erupcí. Všechny hodnoty v tabulkách níže jsou zaokrouhlené na dvě desetinná místa a v každém kroku počítáme s nimi. První tabulka ukazuje, jak daleko leží jednotlivé erupce od přímky odhadnuté z těchto pěti erupcí: ```{r} #| label: tabulka-maleho-vypoctu-residua #| echo: false #| message: false #| warning: false #| tbl-cap: "Pět erupcí: naměřené a odhadnuté čekání, residuum a jeho čtverec (čekání, odhady a residua v minutách)." data_maly_se |> dplyr::transmute( `Erupce` = erupce, `Čekání` = format_cz(cekani_min), `Odhadnuté čekání` = format_cz(odhadnute_cekani), `Residuum` = sub("^-", "−", format_cz(residuum)), `Čtverec residua` = format_cz(residuum_ctverec) ) |> knitr::kable(align = "rrrrr") ``` Sloupec `Residuum` ukazuje residuum (anglicky *residual*), tedy svislou odchylku bodu od přímky odhadnuté z těchto pěti erupcí. Slovo „chyba“ by zde bylo méně přesné, protože skutečnou chybu měření neznáme. Standardní chyba sklonu porovnává dvě informace: jak daleko leží body od přímky a jak široce jsou rozprostřené délky erupcí kolem svého průměru. **Krok 1: `r render_glossary_term("soucet-ctvercu-residui", display = "součet čtverců residuí")` (RSS).** Druhé mocniny všech residuí sečteme. $$ \text{součet čtverců residuí} = (\text{residuum erupce 1})^2 + \cdots + (\text{residuum poslední erupce})^2 $$ Dosadíme residua ze sloupce `Residuum`, umocníme je (sloupec `Čtverec residua`) a sečteme: $$ \begin{aligned} \text{součet čtverců residuí} &= `r paste(text_residua_mocniny[1:3], collapse = " + ")` \\ &\quad + `r paste(text_residua_mocniny[4:5], collapse = " + ")` \\ &= `r paste(text_residua_ctverce, collapse = " + ")` \\ &= `r format_cz_math(soucet_ctvercu_residui_maly_se, nsmall = 2)`. \end{aligned} $$ Tentýž výpočet můžeme zapsat kratším vzorcem: $$ \operatorname{RSS} = \sum_{i=1}^{n} e_i^2. $$ Zde $e_i$ označuje residuum $i$-té erupce, $i$ postupně prochází všechna pozorování a $n$ je jejich počet. Značka $\sum$ znamená „sečti následující hodnotu pro všechna pozorování“. **Krok 2: residuální směrodatná odchylka.** Součet čtverců vydělíme počtem erupcí minus 2 a odmocníme, čímž se vrátíme k minutám čekání. $$ \text{residuální směrodatná odchylka} = \sqrt{ \frac{\text{součet čtverců residuí}}{\text{počet erupcí} - 2} }. $$ $$ \begin{aligned} \text{residuální směrodatná odchylka} &= \sqrt{\frac{`r format_cz_math(soucet_ctvercu_residui_maly_se, nsmall = 2)`}{`r pocet_maly_se` - 2}} \\ &= \sqrt{`r format_cz_math(podil_rss_maly_se, nsmall = 2)`} \\ &= `r format_cz_math(rezidualni_smerodatna_odchylka_maly_se, nsmall = 2)`. \end{aligned} $$ Obecně zapisujeme: $$ \hat{\sigma} = \sqrt{\frac{\operatorname{RSS}}{n - 2}}. $$ Zde $\hat{\sigma}$ označuje residuální směrodatnou odchylku, $\operatorname{RSS}$ právě vypočtený součet čtverců residuí a $n$ počet erupcí; $n-2$ jsou stupně volnosti. **Krok 3: součet čtverců odchylek délek erupcí.** Druhou informací je, jak daleko jsou délky erupcí od své průměrné hodnoty `r format_cz(prumer_delky_maly_se)` minuty. V tabulce jsou délky i jejich odchylky v minutách erupce, jejich čtverce v $\mathrm{min}^2$. ```{r} #| label: tabulka-maleho-vypoctu-odchylky-delky #| echo: false #| message: false #| warning: false #| tbl-cap: "Stejných pět erupcí: délka erupce, její odchylka od průměrné délky a čtverec této odchylky (délky a odchylky v minutách erupce, čtverce odchylek v min²)." data_maly_se |> dplyr::transmute( `Erupce` = erupce, `Délka erupce` = format_cz(delka, digits = 2), `Odchylka od průměru` = sub("^-", "−", format_cz(odchylka_delky)), `Čtverec odchylky` = format_cz(odchylka_delky_ctverec) ) |> knitr::kable(align = "rrrr") ``` $$ \begin{aligned} \text{součet čtverců odchylek délek} &= (\text{délka 1} - \text{průměr})^2 + \cdots \\ &\quad + (\text{délka 5} - \text{průměr})^2. \end{aligned} $$ Dosadíme odchylky ze sloupce `Odchylka od průměru` a jejich čtverce ze sloupce `Čtverec odchylky`: $$ \begin{aligned} \text{součet čtverců odchylek délek erupcí} &= `r paste(text_odchylky_mocniny[1:3], collapse = " + ")` \\ &\quad + `r paste(text_odchylky_mocniny[4:5], collapse = " + ")` \\ &= `r paste(text_odchylky_ctverce, collapse = " + ")` \\ &= `r format_cz_math(soucet_ctvercu_delky_maly_se, nsmall = 2)`. \end{aligned} $$ Obecně zapisujeme: $$ S_{xx} = \sum_{i=1}^{n}(x_i - \bar{x})^2. $$ Zde $x_i$ označuje délku $i$-té erupce, $\bar{x}$ průměrnou délku a $S_{xx}$ součet čtverců odchylek délek erupcí od průměru (anglicky *sum of squares of x*). Čím rozmanitější délky erupcí máme a čím víc jich je, tím větší je $S_{xx}$. **Krok 4: standardní chyba sklonu.** Standardní chyba sklonu spojí obě informace: je to podíl residuální směrodatné odchylky a odmocniny ze součtu čtverců odchylek délek erupcí. $$ \text{standardní chyba sklonu} = \frac{\text{residuální směrodatná odchylka}} {\sqrt{\text{součet čtverců odchylek délek erupcí}}}. $$ $$ \begin{aligned} \text{standardní chyba sklonu} &= \frac{`r format_cz_math(rezidualni_smerodatna_odchylka_maly_se, nsmall = 2)`} {\sqrt{`r format_cz_math(soucet_ctvercu_delky_maly_se, nsmall = 2)`}} \\ &= \frac{`r format_cz_math(rezidualni_smerodatna_odchylka_maly_se, nsmall = 2)`} {`r format_cz_math(odmocnina_sxx_maly_se, nsmall = 2)`} \\ &= `r format_cz_math(standardni_chyba_rucne_maly_se, nsmall = 2)`. \end{aligned} $$ R počítá s nezaokrouhlenými hodnotami a pro těchto pět erupcí vypíše `r format_cz(standardni_chyba_maly_se)`; případný rozdíl v poslední číslici způsobuje jen zaokrouhlování. Je to tentýž vzorec $\operatorname{SE}(\hat{\beta}_1) = \hat{\sigma} / \sqrt{S_{xx}}$ jako v kapitole o standardní chybě z jediného výběru, jen dosazený pro pět erupcí. Pět erupcí dává mnohem větší standardní chybu než všech `r format_cz(nrow(data_gejzir))` erupcí (`r format_cz(standardni_chyba_sklonu)`): $S_{xx}$ má jen pět členů. ::: ### Co říká standardní chyba o našem sklonu? ```{r} #| label: pripravit-navrat-k-modelu #| include: false table_koeficienty <- summary(object = mod_gejzir)$coefficients standardni_chyba_interceptu <- table_koeficienty["(Intercept)", "Std. Error"] procento_se_ze_sklonu <- 100 * standardni_chyba_sklonu / odhad_sklonu ``` Pro otázku, o kolik déle čekáme po delší erupci, potřebujeme z tabulky `Coefficients` modelu ze všech `r format_cz(nrow(data_gejzir))` erupcí jen řádek `delka_erupce_min` a jeho první dva sloupce: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: odhad-a-standardni-chyba-sklonu #| echo: true # Odhad sklonu a jeho standardní chyba z tabulky Coefficients summary(object = mod_gejzir)$coefficients[ "delka_erupce_min", c("Estimate", "Std. Error") ] ``` ```{=typst} ] ``` Obě čísla čteme společně a obě mají jednotku minuty čekání na minutu erupce: - **`Estimate` = `r format_cz(odhad_sklonu)`** je velikost vztahu: dvě erupce, jejichž délky se liší o jednu minutu, se podle modelu v průměru liší čekáním na další erupci o `r format_cz(odhad_sklonu)` minuty. - **`Std. Error` = `r format_cz(standardni_chyba_sklonu)`** je přesnost tohoto odhadu. Kdyby někdo za stejných podmínek zaznamenal jiný výběr `r format_cz(nrow(data_gejzir))` erupcí a odhadl stejný model, dostal by jiný odhad sklonu. Takové odhady by typicky kolísaly zhruba o `r format_cz(standardni_chyba_sklonu)` minuty kolem skutečného sklonu — stejně jako kolísaly sklony hlášení GeyserWatch, jen méně, protože jedno hlášení mělo jen `r format_cz(pocet_erupci_v_hlaseni)` erupcí (typická standardní chyba `r format_cz(typicka_se_hlaseni)`). Ve srovnání s velikostí sklonu je standardní chyba malá: odpovídá asi `r format_cz(procento_se_ze_sklonu, digits = 0)` % odhadu. Náš sklon je tedy odhadnutý poměrně přesně, i když `r format_cz(odhad_sklonu)` není jediná hodnota, kterou data připouštějí. Vlastní standardní chybu má i `r render_glossary_term("intercept", display = "intercept")` (`r format_cz(standardni_chyba_interceptu)` minuty) a čte se stejně; pro naši otázku je ale podstatný sklon. ```{=typst} #block(breakable: false)[ ``` ::: {.box-question title="Kolik erupcí na poloviční standardní chybu?"} Kolik erupcí bychom museli zaznamenat, aby standardní chyba sklonu klesla na polovinu, tedy asi na `r format_cz(standardni_chyba_sklonu / 2)`? Předpokládejte, že residua i délky erupcí budou rozptýlené stejně jako dosud. ::: ```{=typst} ] ``` ::: {.box-result collapse="true" title="Odpověď: čtyřikrát víc"} Počet erupcí stojí ve vzorci pod odmocninou. Aby se jmenovatel zdvojnásobil, musí se počet erupcí zhruba zčtyřnásobit: místo `r format_cz(nrow(data_gejzir))` erupcí bychom jich potřebovali asi `r format_cz(4 * nrow(data_gejzir))`. Přesnost se tedy s dalšími daty zlepšuje stále pomaleji. ::: ::: {.box-data collapse="true" title="Doplňující: co standardní chyba nezahrnuje"} Standardní chyba popisuje `r render_glossary_term("vyberova-variabilita", display = "výběrovou variabilitu")` **za předpokladu, že použitý model a způsob zacházení s pozorováními jsou přiměřené**. Sama automaticky neopraví systematickou chybu měření, nereprezentativní `r render_glossary_term("vyber", display = "výběr")`, `r render_glossary_term("opakovana-mereni", display = "opakovaná měření")` vydávaná za nezávislá ani zřetelné zakřivení (`r render_glossary_term("nelinearita", display = "nelinearitu")`) či změnu rozptýlení residuí. Malá SE tedy neznamená, že studie nemůže být systematicky zavádějící. Můžeme velmi přesně odhadnout vztah v nereprezentativním souboru nebo získat příliš optimistickou SE, pokud jsme přehlédli závislost pozorování. Proto se `r render_glossary_term("nejistota", display = "nejistota")` z tabulky modelu čte společně s `r render_glossary_term("diagnostika-modelu", display = "diagnostikou")` z předchozí lekce a s popisem získání dat. Zda odhad nestojí na jediném pozorování, zkontrolujeme v kapitole „Může jedna erupce změnit sklon?“. ::: *Jaký rozsah sklonů je tedy s našimi daty slučitelný?* ## Jak vyjádřit rozsah `r render_glossary_term("nejistota", display = "nejistoty")`? ### Připomenutí `r render_glossary_term("kvantil", display = "kvantilů")` ```{r} #| label: pripravit-kvantily-sklonu #| include: false vec_kvantily_sklonu <- quantile( x = data_odhady_hlaseni$odhad_sklonu, probs = c(0.025, 0.975), names = FALSE ) ``` V první lekci jsme pomocí kvantilů (anglicky *quantiles*) rozdělovali seřazené hodnoty na části. Například 2,5% kvantil odděluje nejnižších 2,5 % hodnot a 97,5% kvantil nejvyšších 2,5 %. Mezi těmito dvěma hranicemi proto leží prostředních 95 % `r render_glossary_term("rozdeleni", display = "rozdělení")`. Na `r render_glossary_term("histogram", display = "histogramu")` sta `r render_glossary_term("sklon", display = "sklonů")` z GeyserWatch vyznačíme právě tyto dva kvantily: ```{r} #| label: ukazat-kvantily-sklonu #| echo: false #| message: false #| warning: false #| fig-width: 8 #| fig-height: 5.2 #| fig-cap: "Mezi svislými čarami leží prostředních 95 % simulovaných sklonů." #| fig-alt: "Histogram sta sklonů s oranžově podbarvenou střední oblastí mezi 2,5% a 97,5% kvantilem." ggplot2::ggplot( data = data_odhady_hlaseni, mapping = ggplot2::aes(x = odhad_sklonu) ) + ggplot2::annotate( geom = "rect", xmin = vec_kvantily_sklonu[[1]], xmax = vec_kvantily_sklonu[[2]], ymin = -Inf, ymax = Inf, fill = biostat_cols[["orange"]], alpha = 0.35 ) + ggplot2::geom_histogram( bins = 18, fill = biostat_cols[["amethyst"]], color = biostat_cols[["white"]] ) + ggplot2::geom_vline( xintercept = vec_kvantily_sklonu, color = biostat_cols[["orange"]], linetype = "dashed", linewidth = 1 ) + ggplot2::labs( x = "Odhad sklonu (min čekání / min erupce)", y = "Počet hlášení" ) + theme_biostat(base_size = 12) ``` Kvantily vždy patří k určitému rozdělení hodnot; zde popisují `r format_cz(pocet_hlaseni)` simulovaných sklonů. Skutečný výzkumník ale tolik hlášení nemá. `r render_glossary_term("interval-spolehlivosti", display = "Interval spolehlivosti")` (anglicky *confidence interval*, **CI**; setkáte se i s označením konfidenční interval) proto vychází z jediného souboru: použije `r render_glossary_term("odhad", display = "odhad")` sklonu, jeho `r render_glossary_term("standardni-chyba", display = "standardní chybu")` a předpoklady modelu, aby vyjádřil rozsah hodnot `r render_glossary_term("efekt", display = "efektu")` slučitelných s daty. ### Od odhadu a standardní chyby k intervalu ```{r} #| label: pripravit-interval-95 #| include: false hladina_hlavni <- 0.95 kriticka_hodnota_95 <- qt( p = 1 - (1 - hladina_hlavni) / 2, df = df.residual(mod_gejzir) ) polovina_intervalu_95 <- kriticka_hodnota_95 * standardni_chyba_sklonu table_interval_95 <- confint( object = mod_gejzir, parm = "delka_erupce_min", level = hladina_hlavni ) dolni_interval <- table_interval_95[1, 1] horni_interval <- table_interval_95[1, 2] stopifnot( isTRUE(all.equal(dolni_interval, odhad_sklonu - polovina_intervalu_95)), isTRUE(all.equal(horni_interval, odhad_sklonu + polovina_intervalu_95)) ) # The worked steps show three decimals; the displayed arithmetic must # reproduce the displayed results. stopifnot( round(round(kriticka_hodnota_95, 3) * round(standardni_chyba_sklonu, 3), 3) == round(polovina_intervalu_95, 3), round(round(odhad_sklonu, 3) - round(polovina_intervalu_95, 3), 3) == round(dolni_interval, 3), round(round(odhad_sklonu, 3) + round(polovina_intervalu_95, 3), 3) == round(horni_interval, 3) ) ``` Interval vytvoříme ze tří čísel: odhadu sklonu, jeho standardní chyby a `r render_glossary_term("kriticka-hodnota", display = "kritické hodnoty")` (anglicky *critical value*). Kritická hodnota určuje, kolik standardních chyb přidáme na každou stranu od odhadu. Závisí na zvolené `r render_glossary_term("hladina-spolehlivosti", display = "hladině spolehlivosti")` (anglicky *confidence level*), nejčastěji 95 %. Co toto číslo přesně znamená, ukážeme na simulaci v oddílu „Jak číst hladinu spolehlivosti?“. Pro 95% interval je kritická hodnota blízko čísla 2; její přesná hodnota závisí také na počtu `r render_glossary_term("pozorovani", display = "pozorování")` a podrobně ji vysvětlíme v následující lekci. #### 1. Vzdálenost od odhadu k hranici intervalu {.unnumbered .unlisted} Vzdálenost každé hranice intervalu od odhadu se nazývá **polovina šířky intervalu** nebo `r render_glossary_term("mez-chyby", display = "mez chyby")` (anglicky *margin of error*). Samotným hranicím se v odborných textech říká také **meze intervalu**, anglicky *confidence limits*. V této lekci používáme slovo **hranice**. $$ \text{polovina šířky intervalu} = \text{kritická hodnota} \times \text{standardní chyba sklonu}. $$ Pro 95% interval je kritická hodnota `r format_cz(kriticka_hodnota_95, digits = 3)`. Mezivýsledky zde uvádíme na tři desetinná místa, aby výpočet přesně navazoval: $$ \begin{aligned} \text{polovina šířky intervalu} &= `r format_cz_math(kriticka_hodnota_95, digits = 3, nsmall = 3)` \times `r format_cz_math(standardni_chyba_sklonu, digits = 3, nsmall = 3)` \\ &= `r format_cz_math(polovina_intervalu_95, digits = 3, nsmall = 3)`. \end{aligned} $$ Obecně zapisujeme: $$ M = t_{1-\alpha/2,\,n-2}\operatorname{SE}(\hat{\beta}_1). $$ Zde $M$ označuje polovinu šířky intervalu, $t_{1-\alpha/2,\,n-2}$ kritickou hodnotu a $\operatorname{SE}(\hat{\beta}_1)$ standardní chybu odhadu sklonu. Písmeno $\alpha$ je doplněk hladiny do jedné: $\alpha = 1 - \text{hladina}$, takže pro 95% interval je $\alpha = `r format_cz_math(1 - hladina_hlavni, nsmall = 2)`$. Písmeno $t$ označuje rozdělení, ze kterého kritická hodnota pochází, a $n-2$ jsou opět `r render_glossary_term("stupne-volnosti", display = "stupně volnosti")`; obojí vysvětlíme v následující lekci. #### 2. Dolní hranice {.unnumbered .unlisted} $$ \text{dolní hranice} = \text{odhad sklonu} - \text{polovina šířky intervalu}. $$ $$ \begin{aligned} \text{dolní hranice} &= `r format_cz_math(odhad_sklonu, digits = 3, nsmall = 3)` - `r format_cz_math(polovina_intervalu_95, digits = 3, nsmall = 3)` \\ &= `r format_cz_math(dolni_interval, digits = 3, nsmall = 3)`. \end{aligned} $$ Obecně zapisujeme $L = \hat{\beta}_1 - M$. Zde $L$ označuje dolní hranici, $\hat{\beta}_1$ odhad sklonu a $M$ právě vypočtenou polovinu šířky. #### 3. Horní hranice {.unnumbered .unlisted} $$ \text{horní hranice} = \text{odhad sklonu} + \text{polovina šířky intervalu}. $$ $$ \begin{aligned} \text{horní hranice} &= `r format_cz_math(odhad_sklonu, digits = 3, nsmall = 3)` + `r format_cz_math(polovina_intervalu_95, digits = 3, nsmall = 3)` \\ &= `r format_cz_math(horni_interval, digits = 3, nsmall = 3)`. \end{aligned} $$ Obecně zapisujeme $U = \hat{\beta}_1 + M$. Zde $U$ označuje horní hranici. #### 4. Celý interval {.unnumbered .unlisted} $$ \text{95\%}\ \text{interval} = \left[\text{dolní hranice};\,\text{horní hranice}\right]. $$ Po dosazení obou vypočtených hranic: $$ \text{95\%}\ \text{interval} = \left[`r format_cz_math(dolni_interval)`\,;\, `r format_cz_math(horni_interval)`\right]. $$ Celý interval můžeme zapsat kratším vzorcem: $$ I_{1-\alpha} = \left[L;U\right] = \hat{\beta}_1 \pm t_{1-\alpha/2,\,n-2} \operatorname{SE}(\hat{\beta}_1). $$ Značka $I_{1-\alpha}$ označuje interval se zvolenou hladinou spolehlivosti $1-\alpha$; pro 95% interval platí $1-\alpha=`r format_cz_math(hladina_hlavni, nsmall = 2)`$. #### Stejný interval v R {.unnumbered .unlisted} V praxi interval nepočítáme ručně. Funkce `confint()` vypíše intervaly pro všechny koeficienty modelu: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: interval-confint #| echo: true # 95% intervaly spolehlivosti interceptu a sklonu confint(object = mod_gejzir) ``` ```{=typst} ] ``` V řádku `delka_erupce_min` najdeme stejné hranice, jaké jsme vypočítali ručně: `r format_cz(dolni_interval)` a `r format_cz(horni_interval)` minuty čekání na minutu erupce. Názvy sloupců `2.5 %` a `97.5 %` připomínají, že 95% interval nechává na každé straně 2,5 %. ### Jak se interval mění s požadovanou hladinou spolehlivosti? ```{r} #| label: pripravit-intervaly-ruznych-hladin #| include: false vec_hladiny_intervalu <- c(0.80, 0.90, 0.95, 0.99) data_intervaly <- vec_hladiny_intervalu |> purrr::map( .f = ~ { table_jeden_interval <- confint( object = mod_gejzir, parm = "delka_erupce_min", level = .x ) tibble::tibble( hladina = .x, kriticka_hodnota = qt( p = 1 - (1 - .x) / 2, df = df.residual(mod_gejzir) ), dolni = table_jeden_interval[1, 1], horni = table_jeden_interval[1, 2] ) } ) |> purrr::list_rbind() |> dplyr::mutate( popisek = factor( x = paste0(format_cz(100 * hladina, digits = 0), "% interval"), levels = paste0(format_cz(100 * vec_hladiny_intervalu, digits = 0), "% interval") ) ) ``` Hodnota 95 % není jediná možnost. Stejným postupem můžeme vytvořit například 80%, 90% nebo 99% interval; v `confint()` stačí nastavit argument `level`, například `level = 0.80`. ::: {.box-question title="Co čekáte?"} Seřaďte 80%, 90%, 95% a 99% interval od nejužšího po nejširší. Musí interval s vyšší hladinou spolehlivosti být širší, nebo užší? ::: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: tabulka-ruznych-intervalu #| echo: false #| message: false #| warning: false #| tbl-cap: "Intervaly spolehlivosti sklonu pro čtyři zvolené hladiny (min čekání / min erupce)." data_intervaly |> dplyr::transmute( `Hladina` = paste0(format_cz(100 * hladina, digits = 0), " %"), `Kritická hodnota` = format_cz(kriticka_hodnota), `Dolní hranice` = format_cz(dolni), `Odhad` = format_cz(odhad_sklonu), `Horní hranice` = format_cz(horni) ) |> knitr::kable(align = "rrrrr") ``` ```{=typst} ] ``` ```{r} #| label: porovnat-ruzne-intervaly #| echo: false #| message: false #| warning: false #| fig-width: 8 #| fig-height: 4.8 #| fig-cap: "Vyšší požadovaná hladina vede k širšímu intervalu kolem stejného odhadu." #| fig-alt: "Čtyři vodorovné intervaly pod sebou se stejným středem. Osmdesátiprocentní je nejkratší, devadesátidevítiprocentní nejdelší." ggplot2::ggplot( data = data_intervaly, mapping = ggplot2::aes( x = odhad_sklonu, y = popisek ) ) + ggplot2::geom_errorbarh( mapping = ggplot2::aes( xmin = dolni, xmax = horni ), height = 0.25, color = biostat_cols[["amethyst"]], linewidth = 1.4 ) + ggplot2::geom_point( color = biostat_cols[["indigo_velvet"]], size = 3.3 ) + ggplot2::labs( x = "Odhad sklonu (min čekání / min erupce)", y = NULL ) + theme_biostat(base_size = 12) ``` Vyšší hladina spolehlivosti používá větší kritickou hodnotu, a proto přidá na každou stranu od odhadu více standardních chyb. Proto je 80% interval nejužší a 99% interval nejširší. Co přesně vyšší hladina spolehlivosti získá, ukáže následující oddíl. ### Jak číst hladinu spolehlivosti? ```{r} #| label: pripravit-pokryti-intervalu #| include: false # The simulation was generated from the full-data line, so its slope is the # known reference value of this demonstration. referencni_sklon <- odhad_sklonu data_pokryti_intervalu <- data_odhady_hlaseni |> dplyr::select(hlaseni, odhad_sklonu, intervaly) |> tidyr::unnest(cols = intervaly) |> dplyr::mutate( pokryva = dolni_mez <= referencni_sklon & horni_mez >= referencni_sklon ) data_popisky_pokryti <- data_pokryti_intervalu |> dplyr::group_by(hladina) |> dplyr::summarise( pocet_pokrytych = sum(pokryva), .groups = "drop" ) |> dplyr::mutate( panel = paste0( format_cz(100 * hladina, digits = 0), "% interval: ", pocet_pokrytych, " ze ", pocet_hlaseni ) ) |> dplyr::arrange(match(x = hladina, table = vec_hladiny_pokryti)) pocet_pokrytych_95 <- data_popisky_pokryti$pocet_pokrytych[data_popisky_pokryti$hladina == 0.95] pocet_pokrytych_80 <- data_popisky_pokryti$pocet_pokrytych[data_popisky_pokryti$hladina == 0.80] pocet_pokrytych_50 <- data_popisky_pokryti$pocet_pokrytych[data_popisky_pokryti$hladina == 0.50] ``` Vraťme se ke všem `r format_cz(pocet_hlaseni)` hlášením GeyserWatch. Každý návštěvník by ze svého hlášení získal **jeden sklon, jednu standardní chybu a jeden interval** — přesně tak, jak to děláme s našimi daty. Protože hlášení jsou simulovaná, známe referenční sklon `r format_cz(odhad_sklonu)`, ze kterého vznikla. Můžeme proto u každého intervalu ověřit, zda referenční sklon zachytil. Následující trojice panelů používá **stejná hlášení** a mění pouze hladinu spolehlivosti. ```{r} #| label: porovnat-pokryti-sta-intervalu #| echo: false #| message: false #| warning: false #| fig-width: 10 #| fig-height: 6.5 #| fig-cap: !expr paste0("Stejných ", pocet_hlaseni, " simulovaných hlášení s 95%, 80% a 50% intervaly. Plné fialové intervaly zachytily referenční sklon (tmavá čára), šedé přerušované jej minuly.") #| fig-alt: "Tři panely, v každém vodorovné intervaly všech hlášení pod sebou a svislá tmavá čára referenčního sklonu. S nižší hladinou spolehlivosti jsou intervaly kratší a více z nich je šedých a přerušovaných, tedy čáru míjí." data_pokryti_graf <- data_pokryti_intervalu |> dplyr::left_join( data_popisky_pokryti, by = dplyr::join_by(hladina) ) |> dplyr::mutate( panel = factor( x = panel, levels = data_popisky_pokryti$panel ) ) ggplot2::ggplot( data = data_pokryti_graf, mapping = ggplot2::aes( x = odhad_sklonu, y = hlaseni, xmin = dolni_mez, xmax = horni_mez, colour = pokryva, linetype = pokryva ) ) + ggplot2::geom_errorbarh( height = 0, linewidth = 0.8, alpha = 0.9 ) + ggplot2::geom_point(size = 0.7) + ggplot2::geom_vline( xintercept = referencni_sklon, colour = biostat_cols[["graphite"]], linewidth = 0.9 ) + ggplot2::facet_wrap( facets = ggplot2::vars(panel), nrow = 1 ) + ggplot2::scale_colour_manual( values = c( "TRUE" = biostat_cols[["indigo_velvet"]], "FALSE" = biostat_cols[["grey_olive"]] ), guide = "none" ) + ggplot2::scale_linetype_manual( values = c( "TRUE" = "solid", "FALSE" = "22" ), guide = "none" ) + ggplot2::labs( x = "Odhad sklonu a interval (min čekání / min erupce)", y = "Hlášení" ) + theme_biostat(base_size = 14) + ggplot2::theme( axis.text.y = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), panel.grid = ggplot2::element_blank(), strip.text = ggplot2::element_text(size = 12) ) ``` Z `r format_cz(pocet_hlaseni)` intervalů zachytilo referenční sklon `r pocet_pokrytych_95` při hladině 95 %, `r pocet_pokrytych_80` při hladině 80 % a `r pocet_pokrytych_50` při hladině 50 %. Při dalším opakování simulace by počty vyšly trochu jinak, ale kolísaly by kolem hodnot odpovídajících hladině. Podíl intervalů, které hledanou hodnotu zachytí, se nazývá `r render_glossary_term("pokryti-intervalu", display = "pokrytí")` (anglicky *coverage*). **Hladina spolehlivosti popisuje, jak často postup při mnoha opakováních zachytí skutečný sklon** — v simulaci referenční sklon, ve skutečném výzkumu neznámý sklon celého procesu. U 95% intervalů jej asi každý dvacátý mine, a výzkumník, který má jediný interval, nepozná, zda je jeho interval mezi nimi. Je důležité rozlišit dvě věty: - U kvantilů histogramu jsme vyznačili oblast, ve které leží prostředních 95 % simulovaných **odhadů sklonu**. - U 95% intervalu spolehlivosti číslo 95 % říká, jaký podíl opakovaně vytvořených **intervalů** dlouhodobě zachytí skutečný sklon. Jeden vypočtený interval tedy neříká, že uvnitř musí ležet 95 ze 100 sklonů. Neříká ani, že skutečný sklon má po výpočtu 95% pravděpodobnost být uvnitř. Číslo 95 % popisuje úspěšnost postupu při mnoha opakováních. Nižší hladina spolehlivosti, například 80 %, vytvoří z týchž dat užší interval. To ale neznamená, že jsme z dat získali přesnější informaci: pouze jsme přijali, že postup bude častěji minout. **Užší interval znamená vyšší přesnost pouze tehdy, porovnáváme-li intervaly se stejnou hladinou spolehlivosti**, například dva 95% intervaly získané z různě informativních dat. ::: {.box-rule title="Proč právě 95 %?"} Hodnota 95 % není dána přírodním zákonem ani našimi daty. Je to ustálená vědecká konvence a kompromis: vyšší hladina spolehlivosti dává širší interval, nižší užší interval, který ale častěji mine. V následující lekci uvidíme, proč se 95 % často objevuje společně s posuzováním statistických tvrzení. Pro tuto lekci je důležité vědět, že **95 % je standardní volba, nikoli jediná správná volba**. ::: ::: {.box-warning title="Co interval spolehlivosti neříká"} - **Neobsahuje prostředních 95 % erupcí ani čekacích dob.** Popisuje nejistotu odhadu sklonu, nikoli rozdělení jednotlivých pozorování. - **Nedokazuje, že některá hodnota je pravdivá nebo nepravdivá.** Ukazuje, které velikosti efektu jsou s daty a modelem slučitelné. - **Neříká, že skutečný sklon má po výpočtu tohoto intervalu 95% pravděpodobnost ležet uvnitř.** Číslo 95 % patří k dlouhodobému postupu vytváření intervalů. ::: ::: {.box-result title="Vzor zápisu: odhad, standardní chyba a interval"} „Dvě erupce, jejichž délky se lišily o jednu minutu, se podle `r render_glossary_term("linearni-model", display = "lineárního modelu")` v průměru lišily čekáním na další erupci o `r format_cz(odhad_sklonu)` minuty (odhad sklonu = `r format_cz(odhad_sklonu)`, SE = `r format_cz(standardni_chyba_sklonu)`, 95% CI [`r format_cz(dolni_interval)`; `r format_cz(horni_interval)`] minuty čekání na minutu erupce).“ Tento zápis obsahuje velikost vztahu i jeho přesnost. Stále musí být doprovázen popisem dat, modelu a jednotek. ::: ## Jak přesně známe celou přímku? Zatím jsme `r render_glossary_term("nejistota", display = "nejistotu")` vyjadřovali pro jediné číslo — `r render_glossary_term("sklon", display = "sklon")`. Přímka ale dává ještě něco dalšího: pro každou délku erupce odhadnuté průměrné čekání. I výška přímky je tedy `r render_glossary_term("odhad", display = "odhad")`, a proto má vlastní nejistotu. To, jak se průměrné čekání mění s délkou erupce, budeme nazývat **`r render_glossary_term("prumerny-trend", display = "průměrný trend")`**; přímka modelu je jeho odhadem. ```{=typst} #block(breakable: false)[ ``` ::: {.box-question title="Co čekáte?"} Kdybychom do jednoho grafu nakreslili přímky všech `r format_cz(pocet_hlaseni)` hlášení GeyserWatch, kde by se od sebe lišily nejméně? - **A:** U krátkých erupcí. - **B:** U průměrně dlouhých erupcí. - **C:** U dlouhých erupcí. ::: ```{=typst} ] ``` ::: {.box-result collapse="true" title="Odpověď: u průměrně dlouhých erupcí"} Správně je **možnost B**. Přímky se nejméně liší u průměrně dlouhých erupcí a směrem ke krátkým i dlouhým erupcím se rozcházejí. Proč, ukáže následující graf. ::: ### Sto hlášení, sto přímek ```{r} #| label: pripravit-vejir-primek #| include: false # Every report line is drawn over the observed range of eruption lengths. vec_rozsah_delky <- range(data_gejzir$delka_erupce_min) prumerna_delka <- mean(data_gejzir$delka_erupce_min) intercept_vsech_erupci <- coef(object = mod_gejzir)[["(Intercept)"]] data_vejir <- data_odhady_hlaseni |> dplyr::transmute( hlaseni, x = vec_rozsah_delky[1], xend = vec_rozsah_delky[2], y = intercept + odhad_sklonu * x, yend = intercept + odhad_sklonu * xend ) # Heights of all report lines at a short and at the average eruption length delka_kratka <- 2 vec_vyska_kratka <- data_odhady_hlaseni$intercept + data_odhady_hlaseni$odhad_sklonu * delka_kratka vec_vyska_prumerna <- data_odhady_hlaseni$intercept + data_odhady_hlaseni$odhad_sklonu * prumerna_delka ``` ```{r} #| label: vejir-sta-primek #| echo: false #| message: false #| warning: false #| fig-width: 8 #| fig-height: 5 #| fig-cap: !expr paste0("Přímky všech ", pocet_hlaseni, " simulovaných hlášení GeyserWatch (fialově), referenční přímka ze všech erupcí (tmavě) a průměrná délka erupce (přerušovaná čára).") #| fig-alt: "Sto tenkých fialových přímek nakreslených přes sebe. Uprostřed, u svislé přerušované čáry průměrné délky erupce, leží těsně u sebe; směrem ke krátkým i dlouhým erupcím se rozbíhají jako vějíř. Tmavá referenční přímka prochází středem vějíře." ggplot2::ggplot(data = data_vejir) + ggplot2::geom_segment( mapping = ggplot2::aes( x = x, xend = xend, y = y, yend = yend, group = hlaseni ), color = biostat_cols[["indigo_velvet"]], alpha = 0.25, linewidth = 0.5 ) + ggplot2::annotate( geom = "segment", x = vec_rozsah_delky[1], xend = vec_rozsah_delky[2], y = intercept_vsech_erupci + odhad_sklonu * vec_rozsah_delky[1], yend = intercept_vsech_erupci + odhad_sklonu * vec_rozsah_delky[2], color = biostat_cols[["graphite"]], linewidth = 1.1 ) + ggplot2::geom_vline( xintercept = prumerna_delka, linetype = "dashed", color = biostat_cols[["graphite"]] ) + ggplot2::labs( x = "Délka erupce (min)", y = "Odhadnuté čekání (min)" ) + theme_biostat(base_size = 13) ``` Každá fialová přímka patří k jednomu hlášení. Přímky mají různé sklony, ale všechny procházejí blízko středu dat. Přímka odhadnutá funkcí `lm()` má totiž tu vlastnost, že vždy prochází bodem, jehož souřadnicemi jsou průměrná délka erupce a průměrné čekání v datech, ze kterých byla odhadnuta. Tyto průměry se mezi hlášeními liší jen málo. Přímky s různým sklonem se proto kříží blízko průměrné délky erupce (`r format_cz(prumerna_delka)` minuty) a směrem k okrajům se rozbíhají jako vějíř. Po erupcích dlouhých `r format_cz(delka_kratka)` minuty odhadují přímky hlášení průměrné čekání od `r formatC(min(vec_vyska_kratka), format = "f", digits = 1, decimal.mark = ",")` do `r formatC(max(vec_vyska_kratka), format = "f", digits = 1, decimal.mark = ",")` minuty, u průměrné délky erupce jen od `r formatC(min(vec_vyska_prumerna), format = "f", digits = 1, decimal.mark = ",")` do `r formatC(max(vec_vyska_prumerna), format = "f", digits = 1, decimal.mark = ",")` minuty. **Nejistota sklonu se projeví tím víc, čím dál jsme od průměrné délky erupce.** Každé hlášení obsahuje jen `r format_cz(pocet_erupci_v_hlaseni)` erupcí. Vějíř je proto širší než nejistota přímky odhadnuté ze všech `r format_cz(nrow(data_gejzir))` erupcí, stejně jako byla standardní chyba sklonu v jednom hlášení větší než ve výstupu pro všechny erupce. ### `r render_glossary_term("pas-spolehlivosti", display = "Pás spolehlivosti")` kolem přímky Skutečný výzkumník sto hlášení nemá. Nejistotu přímky ale umí odhadnout ze svých dat, podobně jako odhadl `r render_glossary_term("standardni-chyba", display = "standardní chybu")` sklonu: kolem přímky nakreslí `r render_glossary_term("pas-spolehlivosti", display = "pás spolehlivosti")`. Kód, kterým pás nakreslíte sami, ukážeme o kousek níže. ```{r} #| label: pas-spolehlivosti-vsech-erupci #| echo: false #| message: false #| warning: false #| fig-width: 8 #| fig-height: 5 #| fig-cap: "Všechny erupce, přímka modelu a její 95% pás spolehlivosti. Přerušovaná čára označuje průměrnou délku erupce." #| fig-alt: "Bodový graf všech erupcí s fialovou rostoucí přímkou a světle fialovým pásem kolem ní. Pás je nejužší u svislé přerušované čáry průměrné délky erupce, která leží v řidší oblasti mezi shluky krátkých a dlouhých erupcí, a k oběma okrajům se rozšiřuje. Většina bodů leží mimo pás." graf_pas_vsech_erupci <- ggplot2::ggplot( data = data_gejzir, mapping = ggplot2::aes( x = delka_erupce_min, y = cekani_min ) ) + ggplot2::geom_point( color = biostat_cols[["grey_olive"]], alpha = 0.65, size = 2.3 ) + ggplot2::geom_smooth( method = "lm", formula = y ~ x, se = TRUE, color = biostat_cols[["indigo_velvet"]], fill = biostat_cols[["amethyst"]], alpha = 0.2, linewidth = 1.1 ) + ggplot2::geom_vline( xintercept = prumerna_delka, linetype = "dashed", color = biostat_cols[["graphite"]] ) + ggplot2::labs( x = "Délka erupce (min)", y = "Čekání na další erupci (min)" ) + theme_biostat(base_size = 13) graf_pas_vsech_erupci ``` Pás spolehlivosti (anglicky *confidence band*) ukazuje pro každou délku erupce 95% `r render_glossary_term("interval-spolehlivosti", display = "interval spolehlivosti")` pro **průměrné čekání** po erupcích této délky. Skládá se tedy z mnoha intervalů, jednoho pro každou délku erupce. Má podobný tvar jako vějíř přímek: je nejužší u průměrné délky erupce a k okrajům se rozšiřuje. Je ale užší než vějíř, protože vychází ze všech `r format_cz(nrow(data_gejzir))` erupcí. Nejužší místo pásu neleží tam, kde je nejvíc erupcí. Průměrná délka erupce (`r format_cz(prumerna_delka)` minuty) leží v řidší oblasti mezi shluky krátkých a dlouhých erupcí, a pás je přesto právě tam nejužší. Rozhoduje vzdálenost od průměrné délky erupce, ne počet bodů v okolí. Jeden soubor dat dává jeden pás, stejně jako jeden interval z `confint()`. Číslo 95 % čteme stejně jako u sklonu: popisuje úspěšnost postupu při mnoha opakováních, ne jistotu jednoho nakresleného pásu. ### Kdy je pás široký? ```{r} #| label: pripravit-galerii-pasu #| include: false # Same settings as the function defaults used in the presentation. pocet_bodu_galerie <- 40L pocet_bodu_galerie_malo <- 5L vec_rozsah_x_galerie <- c(1, 9) vec_rozsah_x_galerie_uzky <- c(4.5, 5.5) x_sirky_pasu <- 9 list_galerie_pasu <- simulate_band_gallery( n_points = pocet_bodu_galerie, n_points_few = pocet_bodu_galerie_malo, x_range_full = vec_rozsah_x_galerie, x_range_narrow = vec_rozsah_x_galerie_uzky, x_width_check = x_sirky_pasu ) vec_sirka_pasu <- rlang::set_names( list_galerie_pasu$data_panels$sirka_pasu, list_galerie_pasu$data_panels$panel_id ) ``` ```{=typst} #block(breakable: false)[ ``` ::: {.box-question title="Co čekáte?"} Představte si čtyři `r render_glossary_term("simulace", display = "simulované")` soubory dat se stejným skutečným trendem. Ve kterém z nich bude pás spolehlivosti u x = `r x_sirky_pasu` nejširší? - **1:** `r pocet_bodu_galerie` bodů blízko přímky, x od `r vec_rozsah_x_galerie[1]` do `r vec_rozsah_x_galerie[2]`. - **2:** `r pocet_bodu_galerie` bodů daleko od přímky, x od `r vec_rozsah_x_galerie[1]` do `r vec_rozsah_x_galerie[2]`. - **3:** jen `r pocet_bodu_galerie_malo` bodů blízko přímky, x od `r vec_rozsah_x_galerie[1]` do `r vec_rozsah_x_galerie[2]`. - **4:** `r pocet_bodu_galerie` bodů blízko přímky, x jen od `r format_cz(vec_rozsah_x_galerie_uzky[1])` do `r format_cz(vec_rozsah_x_galerie_uzky[2])`. ::: ```{=typst} ] ``` ::: {.box-result collapse="true" title="Odpověď: soubor 4"} Nejširší je pás v souboru 4, přestože jeho body leží blízko přímky. Hodnota x = `r x_sirky_pasu` je od nich daleko. Široké jsou ale i pásy v souborech 2 a 3. Všechny čtyři pásy ukazuje následující graf. ::: ```{r} #| label: galerie-pasu #| echo: false #| message: false #| warning: false #| fig-width: 8 #| fig-height: 5.2 #| fig-cap: !expr paste0("Čtyři simulované soubory dat se stejným skutečným trendem (přerušovaná čára). Fialová přímka je odhadnutá z bodů, světle fialový pás je její 95% pás spolehlivosti, nakreslený i mimo rozsah bodů. V záhlaví panelu je šířka pásu u x = ", x_sirky_pasu, ".") #| fig-alt: "Čtyři panely simulovaných dat se stejnou přerušovanou čarou skutečného trendu. Při mnoha bodech blízko přímky je pás tenký, při bodech daleko od přímky a při pěti bodech je širší. Když body leží jen mezi x = 4,5 a 5,5, je pás u nich úzký, ale mimo ně se rozevírá přes celý graf." plot_band_gallery( list_gallery = list_galerie_pasu, base_size = 12 ) ``` Šířku pásu určují stejné tři ingredience jako standardní chybu sklonu. Soubor 1 je má všechny příznivé a jeho pás je u x = `r x_sirky_pasu` široký jen `r format_cz(vec_sirka_pasu[["1"]], digits = 1)`. Body daleko od přímky (soubor 2) ho rozšíří na `r format_cz(vec_sirka_pasu[["2"]], digits = 1)`, pouhých `r pocet_bodu_galerie_malo` bodů (soubor 3) na `r format_cz(vec_sirka_pasu[["3"]], digits = 1)`. Ve všech souborech se navíc pás rozšiřuje se vzdáleností od průměru x, stejně jako vějíř přímek. Soubor 4 ukazuje, co se děje mimo rozsah dat. Mezi x = `r format_cz(vec_rozsah_x_galerie_uzky[1])` a `r format_cz(vec_rozsah_x_galerie_uzky[2])`, kde body leží, je pás úzký. Kdo by ale přímku prodloužil až k x = `r x_sirky_pasu`, dostane pás široký `r format_cz(vec_sirka_pasu[["4"]], digits = 1)`: data tam o poloze přímky skoro nic neříkají. Použití modelu mimo rozsah pozorovaných hodnot se nazývá `r render_glossary_term("extrapolace", display = "extrapolace")`. I široký pás přitom ukazuje jen nejistotu odhadnuté přímky. Nezaručuje, že vztah mimo data vůbec zůstává přímkový. Naše erupce trvaly od `r format_cz(min(data_gejzir$delka_erupce_min))` do `r format_cz(max(data_gejzir$delka_erupce_min))` minuty. Krátkých i dlouhých erupcí máme dost, a proto je náš pás v celém tomto rozsahu úzký. Odhady průměrného čekání ale platí jen pro erupce v tomto rozsahu. ### Pás v číslech: `predict()` Interval pro jednu zvolenou délku erupce spočítáme funkcí `predict()` s argumentem `interval = "confidence"`. Zvolíme krátkou erupci, průměrnou délku a dlouhou erupci: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: predict-tri-delky #| echo: true # Tři délky erupce: krátká, průměrná a dlouhá data_tri_delky <- data.frame( delka_erupce_min = c(2, mean(x = data_gejzir$delka_erupce_min), 4.5) ) # Odhad průměrného čekání a jeho 95% interval spolehlivosti predict( object = mod_gejzir, newdata = data_tri_delky, interval = "confidence" ) ``` ```{=typst} ] ``` ```{r} #| label: pripravit-sirky-pasu #| include: false mat_tri_delky <- predict( object = mod_gejzir, newdata = data_tri_delky, interval = "confidence" ) vec_sirky_pasu <- mat_tri_delky[, "upr"] - mat_tri_delky[, "lwr"] ``` Sloupec `fit` je odhadnuté průměrné čekání — stejné číslo, jaké leží na přímce. Sloupce `lwr` a `upr` (z anglického *lower* a *upper*) jsou dolní a horní hranice intervalu. Po dvouminutové erupci je interval široký `r format_cz(vec_sirky_pasu[[1]], digits = 1)` minuty, u průměrné délky erupce jen `r format_cz(vec_sirky_pasu[[2]], digits = 1)` minuty a po erupci dlouhé 4,5 minuty `r format_cz(vec_sirky_pasu[[3]], digits = 1)` minuty. Prostřední řádek tedy odpovídá nejužšímu místu pásu. Když stejný výpočet zopakujeme pro mnoho délek erupce v celém pozorovaném rozsahu a hranice spojíme čarami, dostaneme celý pás. Délky erupce a hranice pásu spočítáme takto: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: pas-zakladni-r-vypocet #| echo: true # Padesát délek erupce od nejkratší po nejdelší pozorovanou data_mrizka <- data.frame( delka_erupce_min = seq( from = min(data_gejzir$delka_erupce_min), to = max(data_gejzir$delka_erupce_min), length.out = 50 ) ) # Pro každou délku odhadnuté čekání (fit) a hranice intervalu (lwr, upr) mat_pas <- predict( object = mod_gejzir, newdata = data_mrizka, interval = "confidence" ) ``` ```{=typst} ] ``` Body, přímku a obě hranice pásu nakreslíme takto: ```{r} #| label: pas-zakladni-r-graf #| echo: true #| fig-width: 8 #| fig-height: 4.6 #| fig-alt: "Bodový graf všech erupcí s plnou fialovou odhadnutou přímkou a dvěma přerušovanými fialovými čarami nad a pod ní; čáry jsou si nejblíž u průměrné délky erupce a k okrajům se vzdalují." # Body erupcí plot( x = data_gejzir$delka_erupce_min, y = data_gejzir$cekani_min, xlab = "Délka erupce (min)", ylab = "Čekání na další erupci (min)", pch = 16, col = "#8A8A8A", cex = 1.3, cex.axis = 1.2, cex.lab = 1.3 ) # Odhadnutá přímka plnou čarou, hranice pásu přerušovaně lines( x = data_mrizka$delka_erupce_min, y = mat_pas[, "fit"], col = "#5D2890", lwd = 2 ) lines( x = data_mrizka$delka_erupce_min, y = mat_pas[, "lwr"], col = "#5D2890", lty = 2, lwd = 2 ) lines( x = data_mrizka$delka_erupce_min, y = mat_pas[, "upr"], col = "#5D2890", lty = 2, lwd = 2 ) ``` Funkce `lines()` přidá do již nakresleného grafu čáru spojující zadané body. Pás nekreslíme mimo pozorovaný rozsah délek erupce: tam model nemá data, o která by se opřel. ::: {.box-data collapse="true" title="Doplňující: kratší alternativa pomocí geom_smooth()"} V balíčku `{ggplot2}` nakreslí přímku i pás jediná vrstva `geom_smooth(method = "lm")`. Argument `se = TRUE` přidá pás spolehlivosti, `se = FALSE` nakreslí jen přímku. Zkratka se hodí pro rychlý graf, ale hranice pásu jako čísla nevrátí; ty získáme z `predict()`. ```{r} #| label: pas-ggplot2 #| echo: true #| fig-width: 8 #| fig-height: 4 #| fig-alt: "Stejný graf vytvořený balíčkem ggplot2: šedé body, fialová přímka a světle fialový pás spolehlivosti." # Body erupcí, přímka modelu a 95% pás spolehlivosti ggplot2::ggplot( data = data_gejzir, mapping = ggplot2::aes( x = delka_erupce_min, y = cekani_min ) ) + ggplot2::geom_point(colour = "#8A8A8A") + ggplot2::geom_smooth( method = "lm", formula = y ~ x, se = TRUE, colour = "#5D2890", fill = "#86579E" ) + ggplot2::labs( x = "Délka erupce (min)", y = "Čekání na další erupci (min)" ) + ggplot2::theme_bw(base_size = 14) ``` ::: ### Co pás neříká ```{r} #| label: pripravit-body-mimo-pas #| include: false # Band limits at every observed eruption length mat_pas_pozorovani <- predict( object = mod_gejzir, interval = "confidence" ) vec_mimo_pas <- data_gejzir$cekani_min < mat_pas_pozorovani[, "lwr"] | data_gejzir$cekani_min > mat_pas_pozorovani[, "upr"] pocet_mimo_pas <- sum(vec_mimo_pas) # Mean wait after a 4-minute eruption, used in the answer to the main question delka_ctyri_minuty <- 4 mat_ctyri_minuty <- predict( object = mod_gejzir, newdata = data.frame(delka_erupce_min = delka_ctyri_minuty), interval = "confidence" ) vec_ctyri_minuty_text <- formatC( x = mat_ctyri_minuty[1, ], format = "f", digits = 2, decimal.mark = "," ) ``` ```{=typst} #block(breakable: false)[ ``` ::: {.box-warning title="Na co si dát pozor"} - **Pás neobsahuje 95 % erupcí.** Mimo pás leží `r format_cz(100 * pocet_mimo_pas / nrow(data_gejzir), digits = 0)` % erupcí (`r pocet_mimo_pas` ze `r format_cz(nrow(data_gejzir))`). Pás popisuje, kde může ležet průměrný trend, nikoli kde leží jednotlivá čekání. To, že většina bodů leží mimo pás, proto samo o sobě není známkou špatného modelu. - **Pás neříká, jak dlouho budeme čekat po jedné konkrétní příští erupci.** Jedna erupce se od průměrného trendu odchýlí, stejně jako se od přímky odchylují už pozorované body. Rozsah pro jedno nové čekání musí zahrnout dvě nejistoty současně: nejistotu polohy trendu a rozptýlení jednotlivých čekání kolem něj. Takový rozsah se nazývá `r render_glossary_term("predikcni-interval", display = "predikční interval")` (anglicky *prediction interval*) a je širší než pás; jeho výpočet zde nepotřebujeme. - **Pás neověřuje, že přímka je vhodný model.** Platí jen tehdy, když přímka data přiměřeně popisuje (proto kontrolujeme `r render_glossary_term("residuum", display = "residua")`), a jen v rozsahu pozorovaných délek erupce. ::: ```{=typst} ] ``` ::: {.box-result title="Vzor zápisu: průměrné čekání a pás spolehlivosti"} „Po erupcích dlouhých `r format_cz(delka_ctyri_minuty)` minuty model odhaduje průměrné čekání na další erupci `r vec_ctyri_minuty_text[["fit"]]` minuty (95% CI [`r vec_ctyri_minuty_text[["lwr"]]`; `r vec_ctyri_minuty_text[["upr"]]`] minuty).“ K obrázku: „Plná čára ukazuje přímku lineárního modelu, světlý pás je 95% pás spolehlivosti pro průměrné čekání po erupcích dané délky.“ ::: ## Může jedna erupce změnit sklon? `r render_glossary_term("pas-spolehlivosti", display = "Pás spolehlivosti")` se k okrajům rozšiřuje: čím dál od průměrné délky erupce, tím méně přesně známe výšku přímky. Na okrajích rozsahu má přímka nejméně opory, a jedna erupce ji tam proto může snáz naklonit. Zatím jsme pracovali se `r format_cz(nrow(data_gejzir))` erupcemi, ale při přípravě dat jsme tři erupce vyřadili (viz rámeček o vzniku souboru). Vrátíme je zpět, označíme je písmeny A, B a C a podíváme se, co s přímkou udělají: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: vratit-vyrazene-erupce #| echo: true # Tři erupce vyřazené při přípravě dat (hodnoty z archivu GeyserTimes); # v každém sloupci jsou v pořadí erupce A, B a C data_vyrazene <- data.frame( oznaceni = c("A", "B", "C"), eruption_id = c(1461337, 1468130, 1471579), date = as.Date(x = c("2024-06-04", "2024-07-08", "2024-08-01")), cekani_min = c(94, 121, 106), duration = c(121, 256, 41) ) # Délku erupce převedeme ze sekund na minuty data_vyrazene$delka_erupce_min <- data_vyrazene$duration / 60 # Připojíme je za erupce lekce; stačí dva sloupce modelu vec_sloupce_modelu <- c("delka_erupce_min", "cekani_min") data_gejzir_120 <- rbind( data_gejzir[, vec_sloupce_modelu], data_vyrazene[, vec_sloupce_modelu] ) nrow(x = data_gejzir_120) ``` ```{=typst} ] ``` ```{r} #| label: pripravit-vyrazene-erupce #| include: false # Option texts computed from the typed raw values (seconds -> minutes) vec_delky_vyrazene_text <- format_cz(data_vyrazene$delka_erupce_min, digits = 1) ``` ```{r} #| label: graf-vyrazene-erupce #| echo: false #| message: false #| warning: false #| fig-width: 8 #| fig-height: 5 #| fig-cap: "Erupce lekce (šedě) s přímkou a pásem spolehlivosti odhadnutými z nich a tři vyřazené erupce A, B a C (oranžově)." #| fig-alt: "Bodový graf erupcí lekce s fialovou přímkou a světlým pásem. Tři oranžové body jsou označené písmeny: A mezi krátkými erupcemi nad přímkou, B mezi dlouhými erupcemi nad přímkou a C úplně vlevo, kratší než všechny ostatní erupce, s dlouhým čekáním." ggplot2::ggplot( data = data_gejzir, mapping = ggplot2::aes( x = delka_erupce_min, y = cekani_min ) ) + ggplot2::labs( x = "Délka erupce (min)", y = "Čekání na další erupci (min)" ) + theme_biostat(base_size = 13) + ggplot2::geom_point( color = biostat_cols[["grey_olive"]], alpha = 0.65, size = 2.3 ) + ggplot2::geom_smooth( method = "lm", formula = y ~ x, se = TRUE, color = biostat_cols[["indigo_velvet"]], fill = biostat_cols[["amethyst"]], alpha = 0.2, linewidth = 1.1 ) + ggplot2::geom_point( data = data_vyrazene, color = biostat_cols[["orange"]], size = 3.6 ) + ggplot2::geom_text( data = data_vyrazene, mapping = ggplot2::aes(label = oznaceni), nudge_y = 5, fontface = "bold", size = 5, color = biostat_cols[["graphite"]] ) ``` ```{=typst} #block(breakable: false)[ ``` ::: {.box-question title="Co čekáte?"} Která z vyřazených erupcí by přímku po vrácení do dat naklonila nejvíc? - **A:** Erupce dlouhá asi `r vec_delky_vyrazene_text[[1]]` minuty, po které se čekalo `r data_vyrazene$cekani_min[[1]]` minut. - **B:** Erupce dlouhá asi `r vec_delky_vyrazene_text[[2]]` minuty, po které se čekalo `r data_vyrazene$cekani_min[[2]]` minut. - **C:** Erupce dlouhá jen `r data_vyrazene$duration[[3]]` sekund, po které se čekalo `r data_vyrazene$cekani_min[[3]]` minut. ::: ```{=typst} ] ``` ::: {.box-result collapse="true" title="Odpověď: erupce C"} Nejvíc přímku nakloní **erupce C**. Leží nejdál od průměrné délky erupce a zároveň daleko od přímky. Proč právě tato kombinace rozhoduje, ukážou následující dvě kapitoly. ::: Stejný model odhadneme ze všech `r format_cz(nrow(data_gejzir_120))` erupcí: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: fitovat-model-120-erupci #| echo: true # Stejný vzorec modelu, data se 120 erupcemi mod_gejzir_120 <- lm( formula = cekani_min ~ delka_erupce_min, data = data_gejzir_120 ) coef(object = mod_gejzir_120) ``` ```{=typst} ] ``` ```{r} #| label: pripravit-model-120-erupci #| include: false sklon_120 <- coef(object = mod_gejzir_120)[["delka_erupce_min"]] zmena_sklonu_120 <- odhad_sklonu - sklon_120 ``` Tři erupce navíc změnily `r render_glossary_term("sklon", display = "sklon")` z `r format_cz(odhad_sklonu)` na `r format_cz(sklon_120)` minuty čekání na minutu erupce, tedy o `r format_cz(zmena_sklonu_120)`. To je asi `r format_cz(zmena_sklonu_120 / standardni_chyba_sklonu, digits = 1)`násobek `r render_glossary_term("standardni-chyba", display = "standardní chyby")` sklonu (`r format_cz(standardni_chyba_sklonu)`). Tři `r render_glossary_term("pozorovani", display = "pozorování")` z `r format_cz(nrow(data_gejzir_120))` tak posunula `r render_glossary_term("odhad", display = "odhad")` víc, než kolik odpovídá jeho `r render_glossary_term("vyberova-variabilita", display = "výběrové variabilitě")` vyjádřené standardní chybou. ### Krajní erupce mají delší páku Ve srovnání variant aplikace jsme přímku přirovnali k houpačce: body rozprostřené daleko od sebe ji drží pevně. Platí to i obráceně: jedna erupce daleko od průměrné délky má **delší páku** a může přímku naklonit víc než erupce uprostřed. Tuto vlastnost každého pozorování měří `r render_glossary_term("pakovy-efekt", display = "pákový efekt")` (anglicky *leverage*); mohli jste se s ním setkat už v doplňujícím rámečku minulé lekce. V R jej vrátí funkce `hatvalues()`: ```{r} #| label: spocitat-pakovy-efekt #| echo: true # Pákový efekt každé ze 120 erupcí vec_pakovy_efekt <- hatvalues(model = mod_gejzir_120) # Pákový efekt erupcí A, B a C (poslední tři řádky) data.frame( oznaceni = data_vyrazene$oznaceni, pakovy_efekt = tail(x = vec_pakovy_efekt, n = 3) ) # Rozsah pákového efektu všech erupcí range(vec_pakovy_efekt) ``` ```{r} #| label: pripravit-pakovy-efekt #| include: false vec_paka_abc <- tail(x = vec_pakovy_efekt, n = 3) ``` ```{r} #| label: graf-pakovy-efekt #| echo: true #| fig-width: 8 #| fig-height: 4.2 #| fig-alt: "Pákový efekt proti délce erupce. Body tvoří oblouk tvaru písmene U s nejnižšími hodnotami u průměrné délky erupce, označené svislou přerušovanou čarou. Oranžový bod C úplně vlevo má zdaleka nejvyšší pákový efekt; body A a B jsou na oblouku níž." plot( # Pákový efekt závisí jen na délce erupce x = data_gejzir_120$delka_erupce_min, y = vec_pakovy_efekt, # Místo nad nejvyšším bodem pro jeho popisek ylim = c(0, 1.15 * max(vec_pakovy_efekt)), xlab = "Délka erupce (min)", ylab = "Pákový efekt", pch = 16, col = "#8A8A8A", cex = 1.3, cex.axis = 1.4, cex.lab = 1.4 ) # Tři vrácené erupce zvýrazníme oranžově points( x = tail(x = data_gejzir_120$delka_erupce_min, n = 3), y = tail(x = vec_pakovy_efekt, n = 3), pch = 16, col = "#F3A712", cex = 1.8 ) # Písmena A, B a C nad vrácenými erupcemi text( x = tail(x = data_gejzir_120$delka_erupce_min, n = 3), y = tail(x = vec_pakovy_efekt, n = 3), labels = data_vyrazene$oznaceni, pos = 3, font = 2 ) # Průměrná délka erupce abline( v = mean(x = data_gejzir_120$delka_erupce_min), lty = 2, col = "#5D2890" ) ``` Pákový efekt nemá jednotky. Nejmenší je u erupcí blízko průměrné délky (`r format_cz(min(vec_pakovy_efekt), digits = 3)`) a roste na obě strany, stejně jako se rozšiřuje pás spolehlivosti. Vyjadřuje, jak daleko leží délka jedné erupce od průměrné délky, poměřeno rozprostřením všech délek, tedy druhou ingrediencí standardní chyby. Čekání na další erupci ani `r render_glossary_term("residuum", display = "residuum")` na pákový efekt vliv nemají: erupce C má největší pákový efekt (`r format_cz(vec_paka_abc[[3]], digits = 3)`) jen proto, že trvala nejkratší dobu. Erupce A (`r format_cz(vec_paka_abc[[1]], digits = 3)`) leží dál od průměrné délky než erupce B (`r format_cz(vec_paka_abc[[2]], digits = 3)`). ### Cookova vzdálenost: co kdyby erupce chyběla? Dlouhá páka sama o sobě přímku nenakloní. Erupce s velkým pákovým efektem, která leží přímo na přímce, by sklon nezměnila. Vliv vzniká teprve tehdy, když se dlouhá páka potká s velkým residuem. Obojí spojuje do jednoho čísla `r render_glossary_term("cookova-vzdalenost", display = "Cookova vzdálenost")` (anglicky *Cook's distance*). Pro každou erupci vyjadřuje, jak moc by se změnily `r render_glossary_term("odhadnuta-hodnota", display = "odhadnuté hodnoty")` všech erupcí, kdybychom právě tuto erupci vynechali a přímku odhadli znovu. Jak spolu páka a residuum působí, ukazují simulovaná data na obrázku níže. Ve všech čtyřech panelech jsou stejné šedé body a jeden přidaný oranžový bod, který má malé nebo velké residuum a malý nebo velký pákový efekt. Plná přímka je odhadnutá s přidaným bodem, přerušovaná bez něj. ```{r} #| label: graf-galerie-vlivu #| echo: false #| message: false #| warning: false #| fig-width: 9 #| fig-height: 6.4 #| fig-cap: "Simulovaná data: stejné šedé body a jeden přidaný oranžový bod se čtyřmi kombinacemi pákového efektu a residua. Plná přímka je odhadnutá s přidaným bodem, přerušovaná bez něj. V záhlaví je pákový efekt a Cookova vzdálenost přidaného bodu." #| fig-alt: "Čtyři panely simulovaných dat. Přidaný oranžový bod leží buď uprostřed, nebo daleko vpravo, a buď blízko přímky, nebo daleko od ní. Plná a přerušovaná přímka se výrazně rozcházejí jen v panelu, kde je bod daleko vpravo a zároveň daleko od přímky." list_galerie_vlivu <- simulate_influence_gallery() plot_influence_gallery( list_gallery = list_galerie_vlivu, base_size = 12 ) ``` ```{r} #| label: pripravit-galerii-vlivu #| include: false vec_cook_galerie <- list_galerie_vlivu$data_panels$cookova_vzdalenost ``` Bod uprostřed s malým residuem přímku téměř nezmění (Cookova vzdálenost `r format_cz(vec_cook_galerie[[1]], digits = 3)`). Velké residuum uprostřed ji jen posune nahoru (`r format_cz(vec_cook_galerie[[2]])`) a bod s dlouhou pákou, který leží blízko přímky, ji jen potvrdí (`r format_cz(vec_cook_galerie[[3]])`). Teprve když se dlouhá páka potká s velkým residuem, přímka se výrazně natočí (`r format_cz(vec_cook_galerie[[4]], digits = 1)`). ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: spocitat-cookovu-vzdalenost #| echo: true # Cookova vzdálenost každé ze 120 erupcí vec_cookova_vzdalenost <- cooks.distance(model = mod_gejzir_120) # Erupce A, B a C: délka, residuum, pákový efekt a Cookova vzdálenost data.frame( oznaceni = data_vyrazene$oznaceni, delka_erupce_min = tail(x = data_gejzir_120$delka_erupce_min, n = 3), residuum = tail(x = resid(object = mod_gejzir_120), n = 3), pakovy_efekt = tail(x = vec_pakovy_efekt, n = 3), cookova_vzdalenost = tail(x = vec_cookova_vzdalenost, n = 3) ) ``` ```{=typst} ] ``` ```{r} #| label: pripravit-cookovu-vzdalenost #| include: false vec_cook_abc <- tail(x = vec_cookova_vzdalenost, n = 3) vec_residua_abc <- tail(x = resid(object = mod_gejzir_120), n = 3) pocet_nad_hranici_120 <- sum(vec_cookova_vzdalenost > 4 / nrow(data_gejzir_120)) # The text names A and C as the eruptions above 4/n. stopifnot( identical( names(which(vec_cookova_vzdalenost > 4 / nrow(data_gejzir_120))), c("118", "120") ) ) ``` Erupce A a B leží nad přímkou zhruba stejně daleko (residuum `r format_cz(vec_residua_abc[[1]], digits = 1)` a `r format_cz(vec_residua_abc[[2]], digits = 1)` minuty). Erupce A je ale dál od průměrné délky, má delší páku, a proto vyšší Cookovu vzdálenost (`r format_cz(vec_cook_abc[[1]], digits = 3)` proti `r format_cz(vec_cook_abc[[2]], digits = 3)`). Erupce C spojuje největší pákový efekt s největším residuem a její Cookova vzdálenost (`r format_cz(vec_cook_abc[[3]])`) je mnohonásobně vyšší než u všech ostatních. Tři pojmy proto držíme odděleně. `r render_glossary_term("odlehla-pozorovani", display = "Odlehlé pozorování")` má velké residuum, leží daleko od přímky. Pozorování s velkým pákovým efektem leží daleko od průměrné hodnoty prediktoru. **`r render_glossary_term("vlivne-pozorovani", display = "Vlivné pozorování")`** výrazně mění odhad, a to obvykle tehdy, když se obojí potká. Erupce A a B jsou odlehlé, ale vlivná je jen erupce C. Cookovy vzdálenosti všech erupcí zobrazíme v pořadí řádků tabulky. Přerušovaná čára ukazuje orientační hodnotu 4/n, kde n je počet erupcí v modelu: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: graf-cookova-vzdalenost #| echo: true #| fig-width: 8 #| fig-height: 4.2 #| fig-alt: "Cookova vzdálenost podle pořadí erupce jako svislé čárky. Téměř všechny jsou nízké; poslední čárka, erupce C, vyčnívá vysoko nad všechny ostatní i nad přerušovanou orientační čáru 4/n." # Orientační hodnota 4/n hranice_cook <- 4 / nrow(x = data_gejzir_120) plot( # Každá svislá čárka představuje jednu erupci x = vec_cookova_vzdalenost, type = "h", lwd = 2, col = "#8A8A8A", xlab = "Pořadí erupce v datech", ylab = "Cookova vzdálenost", cex.axis = 1.4, cex.lab = 1.4 ) abline( # Orientační čára, nikoli hranice pro mazání dat h = hranice_cook, lty = 2, lwd = 2, col = "#5D2890" ) ``` ```{=typst} ] ``` Vedle erupce C se všechny ostatní čárky i čára 4/n krčí u nuly. Hodnotu 4/n přesahuje Cookova vzdálenost jen u `r pocet_nad_hranici_120` erupcí: kromě C také u erupce A. Hodnota 4/n není přísná hranice. Označuje pozorování, na která se vyplatí podívat blíž. Jak velký vliv erupce C má, zjistíme přímo: přímku odhadneme znovu bez ní. ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: fitovat-model-bez-nejvlivnejsi #| echo: true # Pořadí erupce s největší Cookovou vzdáleností index_nejvlivnejsi <- which.max(x = vec_cookova_vzdalenost) index_nejvlivnejsi # Stejný model bez této jediné erupce mod_gejzir_bez_c <- lm( formula = cekani_min ~ delka_erupce_min, data = data_gejzir_120[-index_nejvlivnejsi, ] ) coef(object = mod_gejzir_bez_c) ``` ```{=typst} ] ``` ```{r} #| label: pripravit-model-bez-c #| include: false sklon_bez_c <- coef(object = mod_gejzir_bez_c)[["delka_erupce_min"]] # formatC() keeps the trailing zero (0,20), which format_cz() would drop. text_zmena_dvou_erupci <- formatC( x = abs(sklon_bez_c - odhad_sklonu), format = "f", digits = 2, decimal.mark = "," ) ``` Bez jediné erupce C se sklon změní z `r format_cz(sklon_120)` na `r format_cz(sklon_bez_c)` minuty čekání na minutu erupce. Tato jedna erupce tedy způsobila většinu rozdílu mezi modelem se `r format_cz(nrow(data_gejzir_120))` a se `r format_cz(nrow(data_gejzir))` erupcemi; zbylé dvě erupce sklon změnily už jen o `r text_zmena_dvou_erupci`. ### Mění některá z našich erupcí závěr? ```{r} #| label: pripravit-vliv-erupci-lekce #| include: false vec_cook_lekce <- cooks.distance(model = mod_gejzir) pocet_nad_hranici_lekce <- sum(vec_cook_lekce > 4 / nrow(data_gejzir)) max_zmena_sklonu_lekce <- max(abs(dfbeta(model = mod_gejzir)[, "delka_erupce_min"])) podil_zmeny_ze_se <- max_zmena_sklonu_lekce / standardni_chyba_sklonu ``` Stejnou kontrolu uděláme pro model, se kterým pracuje celá lekce. Bez erupce C už je čára 4/n v grafu dobře vidět: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: graf-cookova-vzdalenost-lekce #| echo: true #| fig-width: 8 #| fig-height: 4.2 #| fig-alt: "Cookova vzdálenost erupcí lekce jako svislé čárky s přerušovanou orientační čarou 4/n. Několik čárek čáru mírně přesahuje, žádná výrazně nevyčnívá." plot( # Cookova vzdálenost erupcí, se kterými pracuje lekce x = cooks.distance(model = mod_gejzir), type = "h", lwd = 2, col = "#8A8A8A", xlab = "Pořadí erupce v datech", ylab = "Cookova vzdálenost", cex.axis = 1.4, cex.lab = 1.4 ) abline( # Orientační hodnota 4/n pro 117 erupcí h = 4 / nrow(x = data_gejzir), lty = 2, lwd = 2, col = "#5D2890" ) ``` ```{=typst} ] ``` Funkce `dfbeta()` vrátí pro každou erupci, o kolik by se změnil každý koeficient, kdybychom právě tuto erupci vynechali. Tím se podíváme blíž na všechny erupce najednou: ```{=typst} #block(breakable: false)[ ``` ```{r} #| label: vliv-erupci-lekce #| echo: true # Největší Cookova vzdálenost a orientační hodnota 4/n pro 117 erupcí max(cooks.distance(model = mod_gejzir)) 4 / nrow(x = data_gejzir) # Změna sklonu po vynechání každé jednotlivé erupce vec_zmena_sklonu <- dfbeta(model = mod_gejzir)[, "delka_erupce_min"] # Největší změna sklonu bez ohledu na směr max(abs(x = vec_zmena_sklonu)) ``` ```{=typst} ] ``` Hodnotu 4/n překračuje Cookova vzdálenost u `r pocet_nad_hranici_lekce` erupcí. Největší z těchto hodnot je ale `r format_cz(max(vec_cook_lekce), digits = 3)`, tedy mnohem méně než u erupce C. A hlavně: vynechání kterékoli jednotlivé erupce, včetně těchto `r pocet_nad_hranici_lekce`, změní sklon nejvýš o `r format_cz(max_zmena_sklonu_lekce)` minuty čekání na minutu erupce. To je asi `r format_cz(100 * podil_zmeny_ze_se, digits = 0)` % standardní chyby sklonu (`r format_cz(standardni_chyba_sklonu)`). **Žádná jednotlivá erupce našeho souboru závěr nenese.** Vliv jednoho pozorování proto čteme na stejném měřítku jako nejistotu: změnu sklonu po vynechání erupce porovnáme se standardní chybou. ### Co s vlivným pozorováním? ```{=typst} #block(breakable: false)[ ``` ::: {.box-rule title="Vlivné pozorování nejdřív prověříme, nemažeme"} - Velká Cookova vzdálenost není důvod pozorování smazat. Je to důvod podívat se na ně blíž. - O vyřazení rozhoduje znalost dat: chyba záznamu, jiný druh události nebo pozorování mimo proces, který chceme popsat. Samotná vlivnost důvodem není. - Když si nejsme jistí, uvedeme výsledek s pozorováním i bez něj. ::: ```{=typst} ] ``` U erupce C máme důvod, který nezávisí na tom, jak nakloní přímku. Trvala jen `r data_vyrazene$duration[[3]]` sekund, tedy méně, než je rozsah 1,5 až 5 minut uváděný National Park Service. Mohlo jít o chybu záznamu nebo o neobvyklou událost mimo běžné erupce, které chceme popsat. Je to stejný druh důvodu, kvůli kterému jsme už při přípravě dat ponechali jen čekání od 50 do 127 minut. Erupce A a B takový důvod nemají. Odstranili jsme je jen proto, aby výukový příklad vypadal přehledněji. Jak jsme ukázali, tyto dvě erupce sklon téměř neovlivňují, a podle pravidla výše bychom je po takovém ověření v datech ponechali. Závěry lekce se tím nemění: sklon by se posunul jen o `r text_zmena_dvou_erupci`, mnohem méně, než je jeho standardní chyba. ## Co můžeme říct o Old Faithful? Vraťme se k hlavní otázce: **o kolik déle čekáme po delší erupci — a jak jistě to víme?** V zaznamenaných datech od června do srpna 2024 se dvě erupce, jejichž délky se liší o jednu minutu, podle modelu v průměru liší čekáním na další erupci o `r format_cz(odhad_sklonu)` minuty. To je **velikost vztahu**. `r render_glossary_term("interval-spolehlivosti", display = "95% interval spolehlivosti")` sahá od `r format_cz(dolni_interval)` do `r format_cz(horni_interval)` minuty čekání na minutu erupce. To je **přesnost `r render_glossary_term("odhad", display = "odhadu")`**: data a model jsou slučitelné s tímto rozsahem kladných `r render_glossary_term("sklon", display = "sklonů")`, ne pouze s jediným číslem. Interval se vztahuje ke **sklonu**, tedy k průměrnému rozdílu čekání mezi erupcemi různé délky. Pro zvolenou délku erupce můžeme odhadnout i samotné průměrné čekání: po erupcích dlouhých `r format_cz(delka_ctyri_minuty)` minuty je to podle modelu `r vec_ctyri_minuty_text[["fit"]]` minuty (95% CI [`r vec_ctyri_minuty_text[["lwr"]]`; `r vec_ctyri_minuty_text[["upr"]]`] minuty). Ani tento interval ale neříká, jak dlouho budeme čekat po jedné konkrétní budoucí erupci: ta se od průměrného trendu odchýlí o své vlastní rozptýlení, jak ukázal `r render_glossary_term("pas-spolehlivosti", display = "pás spolehlivosti")` a body mimo něj. Odhad přitom nestojí na jediné erupci: vynechání kterékoli z `r format_cz(nrow(data_gejzir))` erupcí změní sklon nejvýš o `r format_cz(max_zmena_sklonu_lekce)`, tedy asi o `r format_cz(100 * podil_zmeny_ze_se, digits = 0)` % jeho standardní chyby. Erupci dlouhou jen `r data_vyrazene$duration[[3]]` sekund jsme vyřadili z důvodu, který nesouvisí s tím, jak nakloní přímku: leží mimo rozsah délek, které chceme popsat. Interpretace stále závisí na použitém modelu a na způsobu získání dat. Dobrovolnická `r render_glossary_term("pozorovani", display = "pozorování")` nejsou náhodný `r render_glossary_term("vyber", display = "výběr")` všech erupcí. Rozdělení podle měsíců ukázalo, že se odhad sklonu mezi obdobími liší; ze tří měsíců ale nepoznáme, zda jde jen o `r render_glossary_term("vyberova-variabilita", display = "výběrovou variabilitu")`, nebo i o změnu chování gejzíru. Výsledek proto čteme jako opatrný odhad vztahu v zachyceném období, nikoli jako neměnný přírodní zákon. ## Co bude následovat dál V této lekci jsme k `r render_glossary_term("odhad", display = "odhadu")` `r render_glossary_term("sklon", display = "sklonu")` přidali jeho `r render_glossary_term("nejistota", display = "nejistotu")`: `r render_glossary_term("standardni-chyba", display = "standardní chybu")`, `r render_glossary_term("interval-spolehlivosti", display = "interval spolehlivosti")` a `r render_glossary_term("pas-spolehlivosti", display = "pás spolehlivosti")` kolem přímky. Ověřili jsme také, že odhad nestojí na jediné erupci. V následující lekci se zeptáme, jak pomocí dat posoudit konkrétní tvrzení o `r render_glossary_term("parametry-modelu", display = "parametru modelu")` — například zda jsou data slučitelná s tím, že délka erupce s čekáním vůbec nesouvisí. Vrátíme se přitom ke sloupcům `t value` a `Pr(>|t|)`, které jsme zatím přeskočili. ## Shrnutí - Jiný výběr `r render_glossary_term("pozorovani", display = "pozorování")` dává jiný `r render_glossary_term("odhad", display = "odhad")` `r render_glossary_term("sklon", display = "sklonu")`; tuto proměnlivost mezi `r render_glossary_term("vyber", display = "výběry")` nazýváme `r render_glossary_term("vyberova-variabilita", display = "výběrová variabilita")`. - Ve výstupu `summary()` popisuje `Estimate` velikost vztahu a `Std. Error` jeho přesnost. `r render_glossary_term("standardni-chyba", display = "Standardní chyba")` odhaduje z jediného výběru, jak moc by se odhady typicky lišily mezi výběry stejné velikosti. U sklonu ji zvětšují velká `r render_glossary_term("residuum", display = "residua")` a zmenšují ji široce rozprostřené hodnoty `r render_glossary_term("prediktor", display = "prediktoru")` a více pozorování; klesá s odmocninou z jejich počtu. - `r render_glossary_term("interval-spolehlivosti", display = "95% interval spolehlivosti")` z `confint()` ukazuje rozsah sklonů slučitelných s daty a modelem. `r render_glossary_term("hladina-spolehlivosti", display = "Hladina spolehlivosti")` 95 % popisuje úspěšnost postupu při opakování, nikoli podíl pozorování uvnitř intervalu ani důkaz hypotézy. - `r render_glossary_term("pas-spolehlivosti", display = "Pás spolehlivosti")` kolem přímky ukazuje `r render_glossary_term("nejistota", display = "nejistotu")` `r render_glossary_term("prumerny-trend", display = "průměrného trendu")` pro každou délku erupce; hranice získáme z `predict()` s `interval = "confidence"`. Je nejužší u průměrné délky erupce a neobsahuje 95 % erupcí. - `r render_glossary_term("pakovy-efekt", display = "Pákový efekt")` (`hatvalues()`) ukazuje, jak daleko leží pozorování od průměru prediktoru; `r render_glossary_term("cookova-vzdalenost", display = "Cookova vzdálenost")` (`cooks.distance()`) spojuje pákový efekt s velikostí residua. Změnu sklonu po vynechání pozorování porovnáváme se standardní chybou. `r render_glossary_term("vlivne-pozorovani", display = "Vlivné pozorování")` prověříme dřív, než o něm rozhodneme; jen kvůli vlivu je nemažeme. - Vědecký zápis spojuje interpretaci v běžném jazyce s odhadem, standardní chybou, intervalem a jednotkami. ### Závěrečná otázka Někdo vám ukáže jen titulek „Po minutě erupce navíc čekáme o `r format_cz(odhad_sklonu)` minuty déle“. Co potřebujete vidět, než titulku uvěříte? ```{=typst} #block(breakable: false)[ ``` ::: {.box-result collapse="true" title="Možné otázky"} Jak vypadá `r render_glossary_term("bodovy-graf", display = "bodový graf")` a graf residuí — sedí přímka k datům? Je v grafu kolem přímky pás spolehlivosti, aby bylo vidět, jak přesně známe průměrný trend? Kolik erupcí autor měl, kdy a jak byly zaznamenány? Jak přesný je odhad — jaká je standardní chyba a 95% interval spolehlivosti? Nemění závěr jedna erupce, a vyřadil autor některé erupce, případně proč? Mění se sklon mezi obdobími? A tvrdí titulek jen průměrný rozdíl mezi erupcemi různé délky, nebo i to, že delší erupce čekání způsobuje? ::: ```{=typst} ] ``` ## Slovníček pojmů {.unnumbered} Pojmy použité v těchto materiálech a jejich definice. Úplný průběžně aktualizovaný slovníček je dostupný také [online](https://cuni-natur-biostatistics.github.io/slovnik/). ```{r} #| label: slovnicek #| echo: false #| results: asis if ( knitr::is_html_output() ) { glossary::glossary_table() } else { cat( paste0( "Interaktivní slovníček je dostupný online: ", "https://cuni-natur-biostatistics.github.io/slovnik/" ) ) } ```