--- lang: cs format: revealjs: slide-number: true theme: [default, ../theme/presentation_theme.scss, ../theme/presentation_components.scss] embed-resources: true preview-links: false output-file: "presentation" output-ext: "html" highlight-style: github-light include-in-header: "../theme/fonts-include.html" filters: - ../theme/rn-shorthand.lua - ../_extensions/EmilHvitfeldt/roughnotation/rough.lua --- ```{r} #| label: setup #| include: false options(htmltools.dir.version = FALSE) knitr::opts_chunk$set( fig.width = 7, fig.height = 5, fig.align = "center", out.width = "100%" ) library(here) here::i_am("Presentation/presentation.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()) path_materials <- "Presentation/Materials" dir.create( here::here(path_materials), recursive = TRUE, showWarnings = FALSE ) source(here::here("R", "set_r_theme.R")) # Pixel canvases use 150 DPI; match showtext font metrics to the save device. showtext::showtext_opts(dpi = 150) source(here::here("R", "Functions", "format_cz.R")) source(here::here("R", "Functions", "format_cz_math.R")) source(here::here("R", "Functions", "include_local_figure.R")) source(here::here("R", "save_local_figure.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", "compute_se_by_n.R")) source(here::here("R", "Functions", "plot_se_by_n.R")) source(here::here("R", "Functions", "simulate_influence_gallery.R")) source(here::here("R", "Functions", "plot_influence_gallery.R")) source(here::here("R", "Functions", "plot_coverage.R")) source(here::here("R", "Functions", "simulate_band_gallery.R")) source(here::here("R", "Functions", "plot_band_gallery.R")) ``` ```{r} #| label: dynamic-link-generation #| include: false project_name <- basename(here::here()) display_name <- "Biostatistika" github_url <- paste0( "https://github.com/CUNI-NATUR-Biostatistics/", project_name ) slides_url <- paste0( "https://CUNI-NATUR-Biostatistics.github.io/", project_name, "/" ) ``` ## O kolik déle čekáme po delší erupci — a jak jistě to víme? {.course-title .course-title-illustrated .bg-indigo} ::: {.course-title-illustration} ```{r} #| label: zobrazit-titulni-ilustraci #| echo: false #| fig-alt: "Bizon ve slaměném klobouku, svišť v brýlích a krkavec ve fialové pláštěnce pozorují gejzír Old Faithful; drží desky s šedými body a stopky. Na zábradlí sedí výr, vedle telefonu s ikonou gejzíru letí papírové vlaštovky a visí fialové kroužky. Ilustrace vytvořená pomocí AI." #| out-width: "100%" include_local_figure("geyser_titulni_ilustrace.png") ``` ::: ::: {.course-title-mark} ![Logo kurzu Biostatistika](../theme/course-logo-reversed.svg) BIOSTATISTIKA ::: ::: {.course-title-number} L04 ::: ::: {.course-title-formal} Od odhadu k nejistotě ::: ::: {.course-title-link} [cuni-natur-biostatistics.github.io/materialy.html](https://cuni-natur-biostatistics.github.io/materialy.html) ::: ::: {.course-title-caption} Ilustrace vytvořená pomocí AI. ::: ::: {.notes} Titulní ilustrace spojuje postavy z této přednášky: tým bizona, sviště a krkavce, gejzír, hlášení do aplikace GeyserWatch a kroužky. Neobsahuje žádná data ani výsledky. Původ a zadání jsou v `Presentation/Materials/GENERATED_IMAGES.md`. ::: ## Gejzír Old Faithful v číslech ```{r} #| label: nacist-data-pro-otazku #| include: false data_gejzir <- read.csv(file = here::here("data", "old_faithful_2024.csv")) data_gejzir$date <- as.Date(x = data_gejzir$date) data_gejzir$delka_erupce_min <- data_gejzir$duration / 60 ``` :::: {.columns .align-center} ::: {.column width="58%"} ```{r} #| label: zobrazit-fotografii-gejziru #| echo: false #| fig-alt: "Gejzír Old Faithful během erupce: vysoký sloup páry a vody nad šedobílou plošinou, v popředí návštěvníci." include_local_figure("old_faithful_nps_1800.jpg") ``` ::: ::: {.column width="42%"} ::: {.panel-data} [Léto 2024]{.panel-title} - `r format_cz(nrow(data_gejzir))` erupcí v našich datech - erupce trvá přibližně **1,5–5 minut** - čekání na další erupci: obvykle **50–127 minut** - pozorování dobrovolníků ze služby GeyserTimes :::
::: {.takeaway-strip} Napoví délka právě skončené erupce, [jak dlouho budeme čekat]{.text-hgl-orange} na další? ::: ::: :::: ::: {.footer} Foto: [NPS/Jim Peaco](https://www.nps.gov/media/photo/view.htm?id=9C858B86-9CDB-4BB0-8C17-EA567391FAA2), public domain. Rozsahy: National Park Service. ::: ::: {.notes} Otevírám praktickou otázkou návštěvníků. Hlavní otázka z titulní stránky má dvě části: o kolik déle čekáme a jak jistě to víme. Druhou část dnes postupně vybudujeme. ::: {{< include ../pollslive/generated/active.qmd >}} ## Výsledky učení
Po absolvování této přednášky budete schopni: ::: {.incremental} - vysvětlit, proč [jiný soubor pozorování]{.text-hgl-graphite} dává [jiný odhad sklonu]{.text-hgl-indigo_velvet} - najít ve výstupu `summary()` [odhad sklonu]{.text-hgl-indigo_velvet} a jeho [standardní chybu]{.text-hgl-orange}, vysvětlit, co standardní chybu zmenšuje, a rozlišit velikost vztahu od přesnosti odhadu - získat [95% interval spolehlivosti]{.text-hgl-amethyst} sklonu funkcí `confint()` a [pás spolehlivosti]{.text-hgl-amethyst} kolem přímky funkcí `predict()` a interpretovat je bez tvrzení, že obsahují 95 % pozorování nebo dokazují hypotézu - posoudit pomocí [pákového efektu a Cookovy vzdálenosti]{.text-hgl-orange}, zda jedna erupce výrazně mění odhad sklonu - zapsat výsledek jednou větou s [odhadem]{.text-hgl-indigo_velvet}, [SE]{.text-hgl-orange}, [intervalem]{.text-hgl-amethyst} a jednotkami ::: ::: {.notes} Fialová v textu i v grafech znamená odhad modelu a jeho interval, oranžová standardní chybu a šedá naměřená data. Testování hypotéz a p-hodnota přijdou až v následující lekci. ::: # Dají jiné erupce stejnou přímku? ## Každý řádek je jedna erupce ```{r} #| label: ukazat-nacteni-dat #| echo: true #| output-location: column-fragment # Načteme tabulku erupcí data_gejzir <- read.csv(file = "data/old_faithful_2024.csv") # Datum převedeme 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 # Počet erupcí nrow(x = data_gejzir) # Celkový počet chybějících hodnot sum(is.na(x = data_gejzir)) ``` ::: {.footer} [old_faithful_2024.csv](../data/old_faithful_2024.csv): uložte do složky `data` ve svém projektu v R. ::: ::: {.notes} Prediktor je délka erupce v minutách, odezva čekání na další erupci v minutách. Žádný sloupec nemá chybějící hodnoty; podrobnější kontrolu rozsahů najdou studenti ve výukových materiálech. ::: ## Delší erupce, delší čekání? ```{r} #| label: pripravit-graf-vsech-erupci #| include: false plot_data <- 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 = 3 ) + ggplot2::labs( x = "Délka erupce (min)", y = "Čekání na další erupci (min)" ) + theme_biostat(base_size = 20) save_local_figure( plot = plot_data, filename = "vsechna_pozorovani.png", height = 700 ) ``` ```{r} #| label: zobrazit-vsechna-pozorovani #| echo: false #| fig-alt: "Bodový graf všech erupcí: délka erupce na vodorovné ose, čekání na další erupci na svislé. Body tvoří dvě skupiny a celkově stoupají." include_local_figure("vsechna_pozorovani.png") ``` ::: {.takeaway-strip .text-center} Jaký směr má vztah? Kde je nejvíc bodů? Leží body na jediné přímce? ::: ::: {.notes} Krátká diskuze ve dvojicích. Očekávané odpovědi: rostoucí vztah, dvě skupiny erupcí, kolem přímky zůstává rozptýlení. ::: ## Přímka z tohoto léta ```{r} #| label: pripravit-model-vsech-erupci #| include: false mod_gejzir <- lm( formula = cekani_min ~ delka_erupce_min, data = data_gejzir ) odhad_sklonu <- coef(object = mod_gejzir)[["delka_erupce_min"]] plot_model <- plot_data + # 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.6 ) save_local_figure( plot = plot_model, filename = "vsechna_pozorovani_model.png", height = 700 ) ``` ```{r} #| label: zobrazit-model-vsech-pozorovani #| echo: false #| fig-alt: "Stejný bodový graf s fialovou rostoucí přímkou lineárního modelu." include_local_figure("vsechna_pozorovani_model.png") ``` Odhadnutý sklon: [**`r format_cz(odhad_sklonu)` minuty čekání na minutu erupce**]{.text-hgl-indigo_velvet} ::: {.fragment .strip-indigo} **Stejná data + stejný vzorec modelu = vždy stejné odhady.** Funkce `lm()` nic náhodně nelosuje. ::: ::: {.notes} Sklon čteme jako v minulé lekci: dvě erupce, jejichž délky se liší o minutu, se podle modelu v průměru liší čekáním o tento počet minut. Zdůrazním, že model ze stejných dat vždy vrátí stejné číslo; změnit se může jen tehdy, když se změní data. ::: ## Sedí přímka k datům? ```{r} #| label: pripravit-kontrolu-residui #| include: false data_residua <- tibble::tibble( odhadnute_cekani = fitted(object = mod_gejzir), residuum = resid(object = mod_gejzir) ) plot_residua_odhad <- ggplot2::ggplot( data = data_residua, mapping = ggplot2::aes( x = odhadnute_cekani, y = residuum ) ) + ggplot2::labs( x = "Odhadnuté čekání (min)", y = "Residuum (min)" ) + theme_biostat(base_size = 22) + ggplot2::geom_hline( yintercept = 0, linetype = "dashed", color = biostat_cols[["indigo_velvet"]], linewidth = 1 ) + ggplot2::geom_point( color = biostat_cols[["grey_olive"]], alpha = 0.75, size = 3 ) plot_residua_histogram <- ggplot2::ggplot( data = data_residua, mapping = ggplot2::aes(x = residuum) ) + ggplot2::labs( x = "Residuum (min)", y = "Počet erupcí" ) + theme_biostat(base_size = 22) + ggplot2::geom_histogram( bins = 15, fill = biostat_cols[["grey_olive"]], color = biostat_cols[["white"]] ) + ggplot2::geom_vline( xintercept = 0, linetype = "dashed", color = biostat_cols[["indigo_velvet"]], linewidth = 1 ) save_local_figure( plot = plot_residua_odhad, filename = "kontrola_residua_odhad.png", width = 900, height = 700 ) save_local_figure( plot = plot_residua_histogram, filename = "kontrola_residua_histogram.png", width = 900, height = 700 ) ``` :::: {.columns .align-center} ::: {.column width="50%"} ```{r} #| label: zobrazit-residua-odhad #| echo: false #| fig-alt: "Residua proti odhadnutému čekání. Body tvoří dvě skupiny podle krátkých a dlouhých erupcí, obě kolem přerušované nulové čáry, bez oblouku a bez vějíře; několik bodů uprostřed leží výrazně nad nulou." include_local_figure("kontrola_residua_odhad.png") ``` ::: ::: {.column width="50%"} ```{r} #| label: zobrazit-residua-histogram #| echo: false #| fig-alt: "Histogram residuí: většina residuí leží blízko nuly, několik je výrazně kladných." include_local_figure("kontrola_residua_histogram.png") ``` ::: :::: ::: {.fragment} Residua leží kolem nuly bez oblouku a bez vějíře (heteroskedasticity); několik erupcí střední délky leží výrazně nad přímkou. ::: ::: {.fragment .takeaway-strip} [Každý model zkontrolujte, než ho začnete vykládat]{.text-hgl-orange}, i ve vlastních datech: standardní chyba a interval platí jen pro přímku, která data přiměřeně popisuje. ::: ::: {.notes} Krátké opakování z minulé lekce, ne nové učivo. Standardní chyba i interval, ke kterým dnes směřujeme, platí jen tehdy, když přímka data přiměřeně popisuje. Proto residua kontrolujeme vždy, než začneme výsledek interpretovat, i ve vlastních projektech. Stačí `fitted()`, `resid()`, `plot()` a `hist()` z minulé lekce. ::: ## Dostali bychom příště stejnou přímku? ::: {.slide-margin-top-15} ::: {.card-question} ::: {.card-title} Otázka ::: Gejzír pozoruje [někdo jiný v jiných dnech]{.text-hgl-orange}, zaznamená své [body]{.text-hgl-graphite} a fituje stejný model. Jaký [sklon]{.text-hgl-indigo_velvet} dostane? :::
:::: {.columns} ::: {.column width="33%"} ::: {.card-question .card-compact} **A)** Přesně stejný sklon jako my. ::: ::: ::: {.column width="33%"} ::: {.card-question .card-compact} [**B)** Podobný, ale ne stejný sklon.]{.rn-circle-orange} ::: ::: ::: {.column width="33%"} ::: {.card-question .card-compact} **C)** Úplně jiný sklon, nesouvisející s naším. ::: ::: ::::
::: {.takeaway-strip} Zvolte A, B, nebo C a zdůvodněte svou volbu sousedovi. ::: ::: ::: {.notes} Po minutě diskuze zakroužkuji B. Následující slidy s Pepou, Mařenkou a Karlem odpověď ukážou na datech: jiné erupce dají podobnou, ale ne stejnou přímku. ::: ## Pepa pozoroval pouze červen ```{r} #| label: pripravit-pepuv-model #| include: false data_pepa <- data_gejzir[ format(x = data_gejzir$date, format = "%m") == "06", ] mod_pepa <- lm( formula = cekani_min ~ delka_erupce_min, data = data_pepa ) sklon_pepa <- coef(object = mod_pepa)[["delka_erupce_min"]] plot_pepa <- 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 = 3.4 ) + # No band: the confidence band is introduced later in this lesson. ggplot2::geom_abline( intercept = coef(object = mod_pepa)[["(Intercept)"]], slope = sklon_pepa, color = biostat_cols[["indigo_velvet"]], linewidth = 1.6 ) + ggplot2::labs( x = "Délka erupce (min)", y = "Čekání na další erupci (min)", title = "Pepa: červen", subtitle = paste0( "n = ", nrow(data_pepa), ", sklon = ", format_cz(sklon_pepa) ) ) + theme_biostat(base_size = 20) save_local_figure( plot = plot_pepa, filename = "pepa_cerven.png", height = 700 ) ``` ```{r} #| label: zobrazit-pepuv-model #| echo: false #| fig-alt: "Bodový graf Pepových červnových erupcí s fialovou přímkou odhadnutou pouze z nich." include_local_figure("pepa_cerven.png") ``` ::: {.fragment .takeaway-strip} V Pepových datech se dvě erupce, jejichž délky se liší o minutu, v průměru liší čekáním o [**`r format_cz(sklon_pepa)` minuty**]{.text-hgl-indigo_velvet}. ::: ::: {.notes} Pepa ponechal jen červnová pozorování a fitoval stejný model. Jeho výsledek správně popisuje jeho data. ::: ## Co čekáte od Mařenky a Karla? :::: {.columns .align-center} ::: {.column width="44%"} ::: {.panel-olive} [Diskuze se sousedem — 1 min]{.panel-title} Mařenka pozoruje v červenci, Karel v srpnu. Oba použijí stejné proměnné a stejný model jako Pepa. **Budou všechny tři sklony stejné? Proč ano, nebo proč ne?** ::: ::: ::: {.column width="56%"} ```{r} #| label: zobrazit-tym-pozorovatelu #| echo: false #| fig-alt: "Humorná ilustrace tří terénních výzkumníků z Yellowstonu, bizona Pepy, svištice Mařenky a krkavce Karla, z nichž každý drží vlastní desky s několika šedými body. Ilustrace vytvořená pomocí AI." #| out-width: "100%" include_local_figure("geyser_team_photo.png") ``` ::: {.text-size-tiny .text-right} Ilustrace vytvořená pomocí AI. ::: ::: :::: ::: {.notes} Krátká predikce před odhalením tří panelů. ::: ## [Stejný postup]{.text-hgl-indigo_velvet}, tři různé výsledky ```{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 ) ), sklon = purrr::map_dbl( .x = model, .f = ~ coef(.x)[["delka_erupce_min"]] ), popisek = paste0( as.character(vyzkumnik), "\nsklon = ", format_cz(sklon) ) ) |> dplyr::arrange(vyzkumnik) data_mesicni_graf <- data_mesicni |> dplyr::left_join( data_modely_mesicni |> dplyr::select(vyzkumnik, popisek), by = dplyr::join_by(vyzkumnik) ) |> dplyr::mutate( popisek = factor( x = popisek, levels = data_modely_mesicni$popisek ) ) plot_mesice <- 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 = 2.6 ) + # 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.4 ) + ggplot2::facet_wrap( facets = ggplot2::vars(popisek), ncol = 3 ) + ggplot2::labs( x = "Délka erupce (min)", y = "Čekání (min)" ) + theme_biostat(base_size = 20) save_local_figure( plot = plot_mesice, filename = "tri_mesice.png", height = 720 ) ``` ```{r} #| label: zobrazit-tri-mesice #| echo: false #| fig-alt: "Tři bodové grafy vedle sebe pro červen, červenec a srpen, každý s vlastní fialovou přímkou a mírně odlišným sklonem uvedeným v záhlaví." include_local_figure("tri_mesice.png") ``` ::: {.takeaway-strip .fragment} Rozdíly mezi měsíci mohou odrážet [jiné zaznamenané erupce]{.text-hgl-graphite} i skutečnou změnu chování gejzíru. Tři měsíce tyto příčiny neoddělí. ::: ::: {.notes} Měsíce nejsou kontrolovaný experiment. Pro nás je podstatné, že stejný postup dal tři různé odhady. ::: ## Odhad závisí na datech
:::: {.columns} ::: {.column width="48%"} ::: {.panel-result} [Co zůstalo stejné]{.panel-title} - otázka a proměnné - [lineární model]{.text-hgl-indigo_velvet} - funkce `lm()` ::: ::: ::: {.column width="48%"} ::: {.panel-question} [Co se změnilo]{.panel-title} - [pozorované erupce]{.text-hgl-graphite} - [odhadnutá přímka]{.text-hgl-indigo_velvet} a její sklon ::: ::: ::::
::: {.fragment .takeaway-strip} Jednomu souboru pozorování říkáme **výběr** (*sample*). Proměnlivost odhadu mezi možnými výběry je [výběrová variabilita]{.text-hgl-orange} (*sampling variability*). ::: ::: {.notes} Tady poprvé pojmenuji výběr a výběrovou variabilitu, až po zkušenosti se třemi měsíci. Odpovídám tím i na otázku z diskuze: příště bychom pravděpodobně dostali podobnou, ale ne stejnou přímku. ::: # Sto pozorovatelů, sto sklonů ## Skutečnost, výběr, odhad ```{r} #| label: nastavit-hlaseni-geyserwatch #| include: false pocet_hlaseni <- 100L pocet_erupci_v_hlaseni <- 40L ``` [Ve výzkumu]{.text-bold} :::: {.columns .align-center} ::: {.column width="31%"} ::: {.panel-data} [Skutečnost]{.panel-title} Jak se gejzír chová u všech erupcí (cílová populace). Skutečný sklon [nevidíme]{.text-hgl-orange}. ::: ::: ::: {.column width="3%"} **→** ::: ::: {.column width="31%"} ::: {.panel-data} [Výběr]{.panel-title} Erupce, které někdo zaznamenal: našich `r format_cz(nrow(data_gejzir))`. ::: ::: ::: {.column width="3%"} **→** ::: ::: {.column width="31%"} ::: {.panel-result} [Odhad]{.panel-title} Jeden sklon: [`r format_cz(odhad_sklonu)`]{.text-hgl-indigo_velvet}. ::: ::: :::: ::::: {.fragment} [V simulaci]{.text-bold} :::: {.columns .align-center} ::: {.column width="31%"} ::: {.panel-data} [Skutečnost si zvolíme]{.panel-title} Přímka z tohoto léta jako [referenční sklon]{.text-hgl-graphite}. ::: ::: ::: {.column width="3%"} **→** ::: ::: {.column width="31%"} ::: {.panel-data} [`r pocet_hlaseni` výběrů]{.panel-title} Každý ze `r pocet_hlaseni` pozorovatelů: `r pocet_erupci_v_hlaseni` erupcí. ::: ::: ::: {.column width="3%"} **→** ::: ::: {.column width="31%"} ::: {.panel-result} [`r pocet_hlaseni` odhadů]{.panel-title} Sto sklonů. ::: ::: :::: ::::: ::: {.fragment .takeaway-strip} V simulaci [víme, kde je pravda]{.text-hgl-orange}; ve výzkumu máme jeden výběr a skutečnost neznáme. ::: ::: {.notes} Navážu na první lekci: výběr jsou pozorování, která skutečně máme; cílem je říct něco o širším procesu, o skutečnosti. Simulace nám dovolí hrát přírodu: skutečnost si zvolíme sami, a proto můžeme vidět, jak se odhady z různých výběrů od ní liší. Tento rámec budeme potřebovat u standardní chyby i u intervalu spolehlivosti. ::: ## Co udělá [sto dalších pozorovatelů]{.text-hgl-orange}? :::: {.columns .align-center} ::: {.column width="52%"} ::: {.text-bold} Každý pošle do aplikace GeyserWatch [`r pocet_erupci_v_hlaseni` erupcí]{.text-hgl-graphite}, tedy jeden výběr, odhadne [stejný model]{.text-hgl-indigo_velvet} a získá jeden sklon. Jak budou vypadat [všechny sklony společně]{.text-hgl-orange}? ::: ::: {.takeaway-strip} Ve dvojici nakreslete hrubý tvar, který očekáváte po `r pocet_hlaseni` hlášeních. ::: ::: ::: {.column width="48%"} ```{r} #| label: zobrazit-dav-pozorovatelu #| echo: false #| fig-alt: "Humorná ilustrace zástupu návštěvníků a zvířat na chodníku u gejzíru; od každého letí papírová vlaštovka do velkého telefonu s ikonou gejzíru. Ilustrace vytvořená pomocí AI." #| out-width: "100%" include_local_figure("geyserwatch_crowd.png") ``` ::: {.text-size-tiny .text-right} Ilustrace vytvořená pomocí AI. ::: ::: :::: ::: {.panel-data} [Simulace]{.panel-title} Délky erupcí z našich dat; čekání vytvoříme z **referenční přímky** (sklon `r format_cz(odhad_sklonu)`) s odchylkami velkými jako residua. ::: ::: {.notes} Hlášení jsou simulovaná: sto skutečných nezávislých pozorovatelů nemáme. Protože čekání vytváříme z přímky tohoto léta, známe v této ukázce referenční sklon. Skutečný sklon všech erupcí Old Faithful neznáme. ::: ## Sledujte cestu jednoho sklonu ```{r} #| label: vytvorit-hlaseni-a-animaci #| 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 = c(0.95, 0.80, 0.50) ) vec_meze_histogramu <- seq( from = floor(min(data_odhady_hlaseni$odhad_sklonu) * 2) / 2, to = ceiling(max(data_odhady_hlaseni$odhad_sklonu) * 2) / 2, by = 0.25 ) max_cetnost_animace <- max( hist( x = data_odhady_hlaseni$odhad_sklonu, breaks = vec_meze_histogramu, plot = FALSE )$counts ) path_frames <- here::here("Temp", "geyserwatch_frames") unlink(path_frames, recursive = TRUE) dir.create( path_frames, recursive = TRUE, showWarnings = FALSE ) vec_frame_paths <- seq_len(pocet_hlaseni) |> purrr::map_chr( .f = ~ { data_frame_hlaseni <- data_hlaseni |> dplyr::filter(hlaseni == .x) data_frame_odhady <- data_odhady_hlaseni |> dplyr::filter(hlaseni <= .x) data_frame_predchozi_primky <- data_odhady_hlaseni |> dplyr::filter(hlaseni < .x) aktualni_hlaseni <- data_odhady_hlaseni |> dplyr::filter(hlaseni == .x) plot_frame_data <- ggplot2::ggplot( data = data_frame_hlaseni, mapping = ggplot2::aes( x = delka_erupce_min, y = cekani_min ) ) + # No band: simulated report lines in the animation. ggplot2::geom_abline( data = data_frame_predchozi_primky, mapping = ggplot2::aes( intercept = intercept, slope = odhad_sklonu ), inherit.aes = FALSE, color = biostat_cols[["amethyst"]], alpha = 0.10, linewidth = 0.65 ) + ggplot2::geom_point( color = biostat_cols[["grey_olive"]], alpha = 0.75, size = 2.2 ) + ggplot2::geom_abline( intercept = aktualni_hlaseni$intercept, slope = aktualni_hlaseni$odhad_sklonu, color = biostat_cols[["indigo_velvet"]], linewidth = 1.4 ) + ggplot2::coord_cartesian( xlim = range(data_gejzir$delka_erupce_min), ylim = c(45, 135) ) + ggplot2::labs( title = paste0("Hlášení ", .x, "/", pocet_hlaseni), subtitle = paste0( "sklon = ", format_cz(aktualni_hlaseni$odhad_sklonu) ), x = "Délka erupce (min)", y = "Čekání (min)" ) + theme_biostat(base_size = 18) plot_frame_histogram <- ggplot2::ggplot( data = data_frame_odhady, mapping = ggplot2::aes(x = odhad_sklonu) ) + ggplot2::geom_histogram( breaks = vec_meze_histogramu, fill = biostat_cols[["amethyst"]], color = biostat_cols[["white"]] ) + ggplot2::geom_vline( xintercept = aktualni_hlaseni$odhad_sklonu, color = biostat_cols[["indigo_velvet"]], linewidth = 1.4 ) + ggplot2::coord_cartesian( xlim = range(vec_meze_histogramu), ylim = c(0, max_cetnost_animace + 1) ) + ggplot2::labs( title = "Všechny dosavadní sklony", subtitle = "stejný sklon v obou grafech", x = "Sklon (min čekání / min erupce)", y = "Počet hlášení" ) + theme_biostat(base_size = 18) path_frame <- file.path( path_frames, paste0( "frame_", formatC(x = .x, width = 3, flag = "0"), ".png" ) ) save_local_figure( plot = cowplot::plot_grid( plot_frame_data, plot_frame_histogram, ncol = 2 ), filename = basename(path_frame), width = 1600, height = 760, path = path_frame ) path_frame } ) vec_animacni_soubory <- c( rep(vec_frame_paths[[1]], times = 8), vec_frame_paths, rep(vec_frame_paths[[pocet_hlaseni]], times = 16) ) vec_animacni_soubory |> magick::image_read() |> magick::image_animate( fps = 10, loop = 0 ) |> magick::image_write( path = here::here(path_materials, "geyserwatch_sklon_do_histogramu.gif"), format = "gif" ) ``` ```{r} #| label: zobrazit-animaci-hlaseni #| echo: false #| fig-alt: "Animace: vlevo body a fialová přímka jednoho simulovaného hlášení, slabě v pozadí přímky předchozích hlášení; vpravo histogram sklonů, do kterého každé hlášení přidá jeden sklon." include_local_figure("geyserwatch_sklon_do_histogramu.gif") ``` ::: {.strip-indigo} Co se mění z hlášení na hlášení a co zůstává stejné? ::: ::: {.fragment .takeaway-strip} Nový výběr bodů → nová přímka → nový sklon → další díl histogramu. ::: ::: {.notes} Mezi hlášeními se mění vybrané délky erupcí a vytvořená čekání. Stejný zůstává model, referenční přímka i typická velikost odchylek. ::: ## Po `r pocet_hlaseni` hlášeních ```{r} #| label: ulozit-posledni-snimek #| include: false file.copy( from = vec_frame_paths[[pocet_hlaseni]], to = here::here(path_materials, "geyserwatch_100_hlaseni.png"), overwrite = TRUE ) ``` ```{r} #| label: zobrazit-sto-hlaseni #| echo: false #| fig-alt: "Poslední snímek animace: body posledního hlášení se všemi předchozími přímkami v pozadí a úplný histogram sta sklonů." include_local_figure("geyserwatch_100_hlaseni.png") ``` ## Co vidíme — a co ještě nevíme?
:::: {.columns} ::: {.column width="48%"} ::: {.panel-result} [Co vidíme]{.panel-title} - sklony [nejsou stejné]{.text-hgl-orange} - většina leží ve společné oblasti - některé výběry dávají strmější, jiné plošší přímku ::: ::: ::: {.column width="48%"} ::: {.panel-question .fragment} [Co ještě nevíme]{.panel-title} - jak [šířku histogramu]{.text-hgl-orange} popsat jedním číslem - jak ji odhadnout, když máme jen jeden výběr - jaký rozsah sklonů je s daty slučitelný ::: ::: :::: # Jedno číslo pro přesnost ## Jak široce se sklony rozptylují? ```{r} #| label: pripravit-histogram-a-se #| include: false sd_sklonu_hlaseni <- sd(data_odhady_hlaseni$odhad_sklonu) typicka_se_hlaseni <- median(data_odhady_hlaseni$se_sklonu) # The median SE and the SD of the slopes agree to two decimals, so the # slide shows both with three decimals to make clear they differ. plot_histogram_sklonu <- ggplot2::ggplot( data = data_odhady_hlaseni, mapping = ggplot2::aes(x = odhad_sklonu) ) + ggplot2::geom_histogram( breaks = vec_meze_histogramu, fill = biostat_cols[["amethyst"]], color = biostat_cols[["white"]] ) + ggplot2::geom_vline( xintercept = odhad_sklonu, color = biostat_cols[["graphite"]], linewidth = 1 ) + ggplot2::annotate( geom = "label", x = odhad_sklonu, y = max_cetnost_animace + 0.5, label = "referenční sklon", color = biostat_cols[["graphite"]], fill = biostat_cols[["white"]], label.size = 0, size = 6 ) + ggplot2::scale_y_continuous( limits = c(0, max_cetnost_animace + 4), expand = ggplot2::expansion(mult = c(0, 0.02)) ) + ggplot2::labs( x = "Sklon (min čekání / min erupce)", y = "Počet hlášení" ) + theme_biostat(base_size = 20) save_local_figure( plot = plot_histogram_sklonu, filename = "histogram_sklonu.png", width = 1100, height = 800 ) plot_histogram_se <- plot_histogram_sklonu + ggplot2::annotate( geom = "errorbarh", xmin = odhad_sklonu - typicka_se_hlaseni, xmax = odhad_sklonu + typicka_se_hlaseni, y = max_cetnost_animace + 2, height = 0.9, color = biostat_cols[["orange"]], linewidth = 2 ) + ggplot2::annotate( geom = "text", x = odhad_sklonu + typicka_se_hlaseni, y = max_cetnost_animace + 3.2, label = "± 1 SE", hjust = 0, color = biostat_cols[["graphite"]], size = 7 ) save_local_figure( plot = plot_histogram_se, filename = "histogram_a_se.png", width = 1100, height = 800 ) ``` :::: {.columns .align-center} ::: {.column width="58%"} ::: {.r-stack} ```{r} #| label: zobrazit-histogram-sklonu #| echo: false #| fig-alt: "Histogram sta simulovaných sklonů se svislou tmavou čarou s popiskem referenční sklon." include_local_figure("histogram_sklonu.png") ``` ::: {.fragment fragment-index=2} ```{r} #| label: zobrazit-histogram-a-se #| echo: false #| fig-alt: "Stejný histogram s oranžovou úsečkou nad ním, která sahá jednu standardní chybu na každou stranu od referenčního sklonu a pokrývá jeho střední část." include_local_figure("histogram_a_se.png") ``` ::: ::: ::: ::: {.column width="42%"} ::: {.fragment fragment-index=1} ::: {.panel-data} [Ze sta hlášení]{.panel-title} Směrodatná odchylka `r pocet_hlaseni` sklonů: [**`r formatC(x = sd_sklonu_hlaseni, format = "f", digits = 3, decimal.mark = ",")`**]{.text-hgl-orange} ::: ::: ::: {.fragment fragment-index=2} ::: {.panel-result} [Z jediného hlášení]{.panel-title} R vypíše standardní chybu (SE), typicky [**`r formatC(x = typicka_se_hlaseni, format = "f", digits = 3, decimal.mark = ",")`**]{.text-hgl-orange} ::: ::: ::: {.fragment fragment-index=3} ::: {.takeaway-strip} **Standardní chyba (SE)** odhaduje z jediného výběru, jak moc by se sklony typicky lišily mezi výběry stejné velikosti. ::: ::: ::: :::: ::: {.notes} Navážu na první lekci: SD popisuje, jak rozdílná jsou pozorování; SE popisuje, jak přesný je odhad. SD sta sklonů potřebuje všech sto hlášení; SE má každý výzkumník ze svého jediného hlášení. Jednotlivé SE se trefí do šířky histogramu jen přibližně, každá je sama odhadem. Tmavá čára je referenční sklon simulace. Do vzdálenosti jedné SE od ní padnou zhruba dvě třetiny sklonů. Rozsah ± 1 SE ještě není interval spolehlivosti; ten vytvoříme později. Otázku na konci nechám chvíli viset; odpověď dají následující slidy. ::: ## Co dělá sklon stabilnějším? ```{r} #| label: nastavit-srovnani-stability #| include: false pocet_erupci_v_malem_hlaseni <- 10L min_delka_dlouhe_erupce <- 3 ``` ::: {.card-question} ::: {.card-title} Otázka ::: `r pocet_hlaseni` hlášení GeyserWatch, ale vždy se změní jedna věc. Ve které variantě budou mít hlášení [méně rozptýlené sklony]{.text-hgl-orange}? :::
:::: {.columns} ::: {.column width="33%"} ::: {.card-question .card-compact} **Srovnání 1 · body kolem přímky** **A)** body daleko od přímky **B)** body blízko přímky ::: ::: ::: {.column width="33%"} ::: {.card-question .card-compact} **Srovnání 2 · délky erupcí** **A)** délky z celého rozsahu **B)** jen erupce od `r min_delka_dlouhe_erupce` minut ::: ::: ::: {.column width="33%"} ::: {.card-question .card-compact} **Srovnání 3 · počet erupcí** **A)** `r pocet_erupci_v_malem_hlaseni` erupcí v hlášení **B)** `r pocet_erupci_v_hlaseni` erupcí v hlášení ::: ::: ::::
::: {.takeaway-strip} Se sousedem zvolte v každém srovnání A, nebo B. Na odpověď máte minutu. ::: ::: {.notes} Každé srovnání mění jen jednu věc; vše ostatní je jako v GeyserWatch. Hlasování rukou po sloupcích. Odpovědi zatím neprozradím, ukáže je následující slide. ::: ## Tři srovnání, tři odpovědi ```{r} #| label: pripravit-dvojice-stability #| include: false 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 ) # Columns are the three pairs, rows the options A and B; the steadier option # is B, A, B so the answer is not always the same row. data_panely_stability <- tibble::tribble( ~dvojice, ~moznost, ~nastaveni, ~popis, 1L, "A", "scatter_double", "body daleko od přímky", 2L, "A", "geyserwatch", "délky z celého rozsahu", 3L, "A", "eruptions_few", paste0(pocet_erupci_v_malem_hlaseni, " erupcí v hlášení"), 1L, "B", "geyserwatch", "body blízko přímky", 2L, "B", "lengths_long", paste0("jen erupce od ", min_delka_dlouhe_erupce, " min"), 3L, "B", "geyserwatch", paste0(pocet_erupci_v_hlaseni, " erupcí v hlášení") ) |> dplyr::left_join( list_dvojice_stability$data_lines |> dplyr::group_by(nastaveni) |> dplyr::summarise(sd_sklonu = sd(odhad_sklonu)) |> dplyr::ungroup(), by = dplyr::join_by(nastaveni), relationship = "many-to-one", unmatched = "error" ) |> dplyr::mutate( panel = factor( x = paste0(dvojice, moznost, ": ", popis), levels = paste0(dvojice, moznost, ": ", popis) ), text_sd = paste0("SD sklonů ", format_cz(sd_sklonu)) ) data_primky_stability <- data_panely_stability |> 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_panely_stability |> dplyr::left_join( list_dvojice_stability$data_points, by = dplyr::join_by(nastaveni), relationship = "many-to-many" ) # One image per revealed pair; later pairs stay as empty panels so the # layout does not jump between fragments. purrr::walk( .x = 1:3, .f = ~ save_local_figure( plot = plot_stability_pairs( data_points = data_body_stability, data_lines = data_primky_stability, data_panels = data_panely_stability, max_pair = .x, xlim = range(data_gejzir$delka_erupce_min), ylim = range(data_gejzir$cekani_min) + c(-15, 15) ), filename = paste0("dvojice_stability_", .x, ".png"), height = 700 ) ) ``` ::: {.r-stack} ```{r} #| label: zobrazit-dvojice-stability-1 #| echo: false #| fig-alt: "Šest panelů ve dvou řadách. Vyplněné jsou jen panely prvního srovnání: vějíř přímek je při bodech daleko od přímky širší než při bodech blízko přímky. Vpravo dole je v každém vyplněném panelu malý histogram sta sklonů se stejnou osou; u širšího vějíře je širší." include_local_figure("dvojice_stability_1.png") ``` ::: {.fragment} ```{r} #| label: zobrazit-dvojice-stability-2 #| echo: false #| fig-alt: "Stejných šest panelů, nově vyplněné druhé srovnání: při délkách z celého rozsahu je vějíř úzký, při jen dlouhých erupcích jsou přímky krátké a rozbíhají se mnohem víc. Malé histogramy sklonů vpravo dole to ukazují také." include_local_figure("dvojice_stability_2.png") ``` ::: ::: {.fragment} ```{r} #| label: zobrazit-dvojice-stability-3 #| echo: false #| fig-alt: "Všech šest panelů vyplněných; ve třetím srovnání je vějíř hlášení s menším počtem erupcí širší než vějíř hlášení s větším počtem erupcí. Malé histogramy sklonů vpravo dole to ukazují také." include_local_figure("dvojice_stability_3.png") ``` ::: ::: ::: {.fragment .takeaway-strip} Stabilnější sklon: body [blízko přímky]{.text-hgl-orange}, délky z [celého rozsahu]{.text-hgl-orange} (body daleko od sebe drží přímku jako houpačku), [více erupcí]{.text-hgl-orange}. Malý histogram v rohu: užší znamená menší SD sklonů. ::: ::: {.notes} Malý histogram v rohu panelu: sklony ze sta hlášení, ve všech panelech se stejnou osou; užší histogram znamená menší SD. Odpovědi 1B, 2A, 3B. Srovnání 2 je houpačka: body daleko od sebe po vodorovné ose přímku drží pevně, body natlačené k sobě ji nechají snadno naklonit. Všechny tři vlastnosti výzkumník vidí i v jediném souboru dat, a právě z nich R standardní chybu počítá. ::: ## Tři ingredience standardní chyby
$$ \text{SE sklonu} = \frac{\text{typická velikost residua}} {\text{SD délek erupcí} \times \sqrt{\text{počet erupcí} - 1}} $$
:::: {.columns} ::: {.column width="33%"} ::: {.fragment fragment-index=1} ::: {.panel-data} [1 · Čitatel]{.panel-title} Body daleko od přímky → [větší SE]{.text-hgl-orange} ::: ::: ::: ::: {.column width="33%"} ::: {.fragment fragment-index=2} ::: {.panel-data} [2 · Jmenovatel]{.panel-title} Délky z celého rozsahu → [menší SE]{.text-hgl-orange} ::: ::: ::: ::: {.column width="33%"} ::: {.fragment fragment-index=3} ::: {.panel-data} [3 · Jmenovatel]{.panel-title} Více erupcí → [menší SE]{.text-hgl-orange}, ale jen s odmocninou ::: ::: ::: :::: ::: {.notes} Rozptýlení délek je jejich SD z první lekce. „Minus 1“ pochází z toho, jak se počítá SD; u velkého počtu erupcí je to skoro odmocnina z počtu. Jak získat typickou velikost residua, ukáže následující slide. Ve výukových materiálech je stejný vzorec i v obvyklém zápisu se symboly. ::: ## Jak velké je typické residuum? ```{r} #| label: pripravit-rse-ze-sse #| include: false # Displayed values are rounded at every step so the arithmetic on the slide # reproduces the displayed residual standard deviation. rss_zobrazeny <- round(sum(resid(object = mod_gejzir)^2), 1) stupne_volnosti_gejzir <- nrow(data_gejzir) - 2L podil_rss_zobrazeny <- round(rss_zobrazeny / stupne_volnosti_gejzir, 2) rsd_ze_rss_zobrazena <- round(sqrt(podil_rss_zobrazeny), 2) stopifnot( abs(rsd_ze_rss_zobrazena - round(sigma(object = mod_gejzir), 2)) < 0.005 ) ``` ::: {.text-center} Residua jsou svislé odchylky erupcí od přímky. Jejich typickou velikost shrneme podobně jako SD. :::
:::: {.columns} ::: {.column width="33%"} ::: {.panel-data} [1 · Sečteme čtverce]{.panel-title} RSS = [**`r format_cz(rss_zobrazeny, digits = 1)` min²**]{.text-hgl-graphite} Součet čtverců residuí (v minulé lekci SSE): `lm()` volí přímku s nejmenším RSS. ::: ::: ::: {.column width="33%"} ::: {.fragment fragment-index=1} ::: {.panel-data} [2 · Vydělíme n − 2]{.panel-title} `r format_cz(rss_zobrazeny, digits = 1)` / `r stupne_volnosti_gejzir` = **`r format_cz(podil_rss_zobrazeny)` min²** ::: ::: ::: ::: {.column width="33%"} ::: {.fragment fragment-index=2} ::: {.panel-data} [3 · Odmocníme]{.panel-title} √`r format_cz(podil_rss_zobrazeny)` = [**`r format_cz(rsd_ze_rss_zobrazena)` min**]{.text-hgl-orange} ::: ::: ::: :::: ::: {.fragment fragment-index=3} $$ \text{residuální směrodatná odchylka} = \sqrt{\frac{\text{RSS}}{n - 2}} = \sqrt{\frac{`r format_cz_math(rss_zobrazeny, digits = 1, nsmall = 1)`}{`r stupne_volnosti_gejzir`}} = `r format_cz_math(rsd_ze_rss_zobrazena, nsmall = 2)`\ \text{min} $$ ::: ::: {.fragment fragment-index=4 .takeaway-strip} [Residuální směrodatná odchylka]{.text-hgl-orange} (ve výstupu R *Residual standard error*): erupce leží typicky asi `r format_cz(rsd_ze_rss_zobrazena)` min nad přímkou nebo pod ní. ::: ::: {.notes} Odmocněním se vracíme z min² k minutám čekání. n − 2 jsou stupně volnosti; proč se odečítá právě 2, vysvětlíme v následující lekci. R vrátí stejné číslo funkcí sigma() a ve výstupu summary() jej uvidíme jako Residual standard error. ::: ## Spočítejte standardní chybu sami ```{r} #| label: pripravit-vypocet-se #| include: false # Displayed values are rounded so students can redo the arithmetic on a # calculator and reach the displayed result. rezidualni_sd_zobrazena <- round(sigma(object = mod_gejzir), 2) 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 ) ``` :::: {.columns .align-center} ::: {.column width="45%"} ::: {.panel-data} [Všech `r format_cz(nrow(data_gejzir))` erupcí]{.panel-title} - residuální směrodatná odchylka: **`r format_cz(rezidualni_sd_zobrazena)` min** - SD délek erupcí: **`r format_cz(sd_delky_zobrazena, digits = 3)` min** - počet erupcí: **`r nrow(data_gejzir)`** ::: ::: {.panel-olive} [Dvojice — 2 minuty]{.panel-title} $$ \text{SE} = \frac{\text{residuální SD}}{\text{SD délek} \times \sqrt{n - 1}} $$ Kalkulačka v telefonu stačí. ::: ::: ::: {.column width="55%"} ::: {.fragment} $$ \begin{aligned} \text{SE} &= \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}} \\[4pt] &= \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} $$ ::: ::: {.fragment .takeaway-strip} [SE sklonu = `r format_cz(standardni_chyba_zobrazena)` min/min]{.text-hgl-orange} — z jediného souboru dat. ::: ::: :::: ::: {.notes} Typickou velikost residua R vrátí funkcí sigma(), v dalším výstupu ji najdeme jako Residual standard error. SD délek známe z první lekce. Kdo nemá kalkulačku, odhadne odmocninu ze `r nrow(data_gejzir) - 1` zhruba na `r format_cz(odmocnina_n_zobrazena, digits = 1)`. ::: ## Proč má hlášení se `r pocet_erupci_v_hlaseni` erupcemi větší SE? ```{r} #| label: pripravit-se-podle-poctu #| include: false # SE from the formula with the residual SD and spread of lengths of our data; # only the number of eruptions changes. data_krivka_se <- tibble::tibble( pocet = seq(from = 20, to = 4L * nrow(data_gejzir) + 30, by = 2) ) |> dplyr::mutate( se = compute_se_by_n( mod_reference = mod_gejzir, vec_predictor = data_gejzir$delka_erupce_min, vec_n = pocet ) ) data_body_se <- tibble::tibble( pocet = c( pocet_erupci_v_hlaseni, nrow(data_gejzir), 4L * nrow(data_gejzir) ), popis = c("hlášení", "naše data", "4 × naše data") ) |> dplyr::mutate( se = compute_se_by_n( mod_reference = mod_gejzir, vec_predictor = data_gejzir$delka_erupce_min, vec_n = pocet ), text = paste0(popis, ": n = ", pocet, ", SE = ", format_cz(se)) ) se_hlaseni_vzorec <- compute_se_by_n( mod_reference = mod_gejzir, vec_predictor = data_gejzir$delka_erupce_min, vec_n = pocet_erupci_v_hlaseni ) purrr::walk( .x = 2:3, .f = ~ save_local_figure( plot = plot_se_by_n( data_curve = data_krivka_se, data_points = data_body_se, n_points = .x ), filename = paste0("se_podle_poctu_", .x, ".png"), height = 620 ) ) ``` ::: {.r-stack} ```{r} #| label: zobrazit-se-podle-poctu #| echo: false #| fig-alt: "Oranžová křivka standardní chyby sklonu podle vzorce klesá s počtem erupcí, nejdřív rychle, pak pomaleji. Vyznačené body: hlášení a naše data." include_local_figure("se_podle_poctu_2.png") ``` ::: {.fragment} ```{r} #| label: zobrazit-se-podle-poctu-ctyrikrat #| echo: false #| fig-alt: "Stejná křivka s dalším bodem pro čtyřnásobek našich erupcí, kde je standardní chyba poloviční." include_local_figure("se_podle_poctu_3.png") ``` ::: ::: ::: {.fragment .takeaway-strip} Hlášení: vzorec [**`r format_cz(se_hlaseni_vzorec)`**]{.text-hgl-orange}, simulace **`r format_cz(sd_sklonu_hlaseni)`**. [Čtyřikrát víc erupcí]{.text-hgl-orange} → jen [poloviční SE]{.text-hgl-orange}. ::: ::: {.notes} Stejná residua a stejné rozptýlení délek jako v našich datech, mění se jen počet erupcí. Vzorec z jediného souboru dat odhadne šířku histogramu, kterou jsme jinak potřebovali sto hlášení. Před dalším bodem se zeptám: kolik erupcí bychom potřebovali na poloviční SE? Přesnost s dalšími daty roste stále pomaleji. ::: ## Co vypíše `summary()`? :::: {.columns} ::: {.column width="58%" .text-size-smaller} ```{r} #| label: zobrazit-souhrn-modelu #| echo: true # Model ze všech erupcí a jeho souhrn mod_gejzir <- lm( formula = cekani_min ~ delka_erupce_min, data = data_gejzir ) summary(object = mod_gejzir) ``` ::: ::: {.column width="42%" .text-size-smaller} ::: {.panel-olive} [Dvojice — 1 min]{.panel-title} Co z výstupu už znáte? Co vidíte poprvé? ::: ::: {.incremental} - `Call`: vzorec modelu - `Residuals`: pět čísel o residuích - `Coefficients`: tabulka odhadů, řádek pro intercept a pro sklon - `Residual standard error`: [residuální směrodatná odchylka]{.text-hgl-orange} - zbytek: v dalších lekcích ::: ::: :::: ::: {.notes} Je to poprvé, co vidíme celý výstup summary() pro model. Nechám dvojice minutu hledat, co znají: odhady interceptu a sklonu z coef(), residua z minulé lekce. Pak projdeme části shora dolů. t value a Pr(>|t|) patří k testování v další lekci; Multiple R-squared a F-statistic přijdou při porovnávání modelů. Residual standard error není standardní chyba sklonu, i když se tak anglicky jmenuje. Tabulku Coefficients si na dalším slidu přiblížíme. ::: ## Tabulka `Coefficients`: velikost a přesnost ```{r} #| label: pripravit-souhrn-modelu #| include: false standardni_chyba_sklonu <- summary(object = mod_gejzir)$coefficients[ "delka_erupce_min", "Std. Error" ] ``` ::: {.text-size-large} ```{r} #| label: zobrazit-tabulku-koeficientu #| echo: true # Jen tabulka Coefficients ze souhrnu modelu, na 3 platné číslice signif( x = summary(object = mod_gejzir)$coefficients, digits = 3 ) ``` :::
:::: {.columns} ::: {.column width="45%"} ::: {.panel-olive} [Dvojice — 1 min]{.panel-title} Najděte číslo, které popisuje 1. **velikost vztahu**, 2. **přesnost tohoto odhadu**. ::: ::: ::: {.column width="55%"} ::: {.fragment} ::: {.panel-result} [Řádek `delka_erupce_min`]{.panel-title} - `Estimate` = [`r format_cz(odhad_sklonu)` min/min]{.text-hgl-indigo_velvet}: velikost vztahu - `Std. Error` = [`r format_cz(standardni_chyba_sklonu)` min/min]{.text-hgl-orange}: stejné číslo, jaké jste spočítali ::: ::: ::: :::: ::: {.notes} Řádek (Intercept) patří k interceptu; pro naši otázku je podstatný řádek sklonu. Sloupce t value a Pr(>|t|) zatím přeskočíme; vrátíme se k nim v následující lekci. ::: ## Velikost vztahu, nebo přesnost? ::: {.takeaway-strip} Ve dvojici přiřaďte každou otázku k jednomu ze dvou čísel. ::: :::: {.columns} ::: {.column width="50%"} ::: {.card-question} „O kolik minut se v průměru liší čekání po erupcích lišících se o minutu?“ [Velikost vztahu (`Estimate`)]{.fragment .text-hgl-indigo_velvet fragment-index=1} ::: ::: {.card-question} „Jak moc by se odhad lišil mezi výběry?“ [Přesnost odhadu (`Std. Error`)]{.fragment .text-hgl-orange fragment-index=1} ::: ::: ::: {.column width="50%"} ::: {.card-question} „Jak těsně se možné odhady shlukují?“ [Přesnost odhadu (`Std. Error`)]{.fragment .text-hgl-orange fragment-index=1} ::: ::: {.card-question} „Je vztah prudký, nebo pozvolný?“ [Velikost vztahu (`Estimate`)]{.fragment .text-hgl-indigo_velvet fragment-index=1} ::: ::: :::: ::: {.notes} Odpovědi se objeví najednou po krátké diskuzi. Zdůrazním, že velký sklon může být odhadnutý nepřesně a malý velmi přesně. ::: ## O kolik déle čekáme — a jak jistě to víme?
:::: {.columns} ::: {.column width="48%"} ::: {.panel-result} [Už víme]{.panel-title} - velikost vztahu: o [**`r format_cz(odhad_sklonu)` min**]{.text-hgl-indigo_velvet} čekání na minutu erupce - přesnost odhadu: [**SE = `r format_cz(standardni_chyba_sklonu)` min/min**]{.text-hgl-orange} - SE klesá s více erupcemi a s širším rozsahem délek ::: ::: ::: {.column width="48%"} ::: {.fragment} ::: {.panel-question} [Ještě nevíme]{.panel-title} - jaký rozsah sklonů je s daty slučitelný - jak přesně známe celou přímku - zda závěr nestojí na jediné erupci ::: ::: ::: :::: ::: {.notes} Vracím se k otázce z titulní stránky. První část už umíme zodpovědět číslem, druhou jen napůl: SE říká, jak moc by odhad kolísal, ale ještě ne, jaký rozsah sklonů je s daty slučitelný. To jsou tři zbývající části přednášky. ::: # Jaké sklony jsou s daty slučitelné? ## ::: {.slide-margin-top-15} ::: {.card-question} ::: {.card-title} Otázka ::: Každé ze `r pocet_hlaseni` hlášení GeyserWatch si kolem svého sklonu vyznačí rozsah [odhad ± asi 2 SE]{.text-hgl-orange}. Při mnoha opakováních zachytí referenční sklon kolik takových intervalů? :::
:::: {.columns} ::: {.column width="25%"} ::: {.card-question .card-compact} **A)** asi 50 ::: ::: ::: {.column width="25%"} ::: {.card-question .card-compact} **B)** asi 70 ::: ::: ::: {.column width="25%"} ::: {.card-question .card-compact} [**C)** asi 95]{.rn-circle-orange} ::: ::: ::: {.column width="25%"} ::: {.card-question .card-compact} **D)** všech `r pocet_hlaseni` ::: ::: ::::
::: {.takeaway-strip} Zvolte A, B, C, nebo D a zdůvodněte svou volbu sousedovi. ::: ::: ::: {.notes} Navážu na histogram: do vzdálenosti jedné SE od referenčního sklonu padly asi dvě třetiny sklonů. Teď kolem každého sklonu dáme dvě SE na každou stranu. Po diskuzi zakroužkuji C; následující slide to ukáže na stejných hlášeních. ::: ## Sto pozorovatelů, sto 95% intervalů ```{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( vysledek = dplyr::if_else( dolni_mez <= referencni_sklon & horni_mez >= referencni_sklon, "Zachytil referenční sklon", "Minul referenční sklon" ) ) data_pocty_pokrytych <- data_pokryti_intervalu |> dplyr::group_by(hladina) |> dplyr::summarise( pocet = sum(vysledek == "Zachytil referenční sklon"), .groups = "drop" ) pocet_pokrytych_95 <- data_pocty_pokrytych$pocet[data_pocty_pokrytych$hladina == 0.95] rozsah_osy_intervalu <- range( data_pokryti_intervalu$dolni_mez[data_pokryti_intervalu$hladina == 0.95], data_pokryti_intervalu$horni_mez[data_pokryti_intervalu$hladina == 0.95] ) save_local_figure( plot = plot_coverage( data_coverage = data_pokryti_intervalu, level = 0.95, slope_reference = referencni_sklon, xlim = rozsah_osy_intervalu ), filename = "sto_intervalu_pokryti.png", height = 760 ) ``` ::: {.text-center} Stejných `r pocet_hlaseni` hlášení, každé se svým rozsahem odhad ± asi 2 SE (přesně: 95% interval z `confint()`). ::: ```{r} #| label: zobrazit-sto-intervalu #| echo: false #| out-width: "88%" #| fig-alt: "Sto vodorovných intervalů pod sebou se svislou tmavou čarou referenčního sklonu. Téměř všechny jsou fialové a čáru protínají; několik šedých ji míjí." include_local_figure("sto_intervalu_pokryti.png") ``` ::: {.strip-orange} **`r pocet_pokrytych_95` ze `r pocet_hlaseni`** 95% intervalů zachytilo referenční sklon. Ve skutečném výzkumu skutečný sklon neznáme. ::: ::: {.notes} Referenční sklon známe jen proto, že jsme z něj čekání v simulaci vytvořili. Výzkumník s jedním intervalem nepozná, zda je jeho interval mezi těmi, které minuly. ::: ## [95 %]{.text-hgl-orange} patří k [postupu]{.text-hgl-indigo_velvet} :::: {.columns .align-center} ::: {.column width="40%"} ```{r} #| label: zobrazit-hazeni-krouzku #| echo: false #| fig-alt: "Humorná ilustrace bizona, který hází fialové kroužky na tmavý kolík; většina kroužků kolík obepíná, dva leží vedle. Krkavec jako rozhodčí počítá oblázky. Ilustrace vytvořená pomocí AI." #| out-width: "100%" include_local_figure("geyser_ring_toss.png") ``` ::: {.text-size-tiny .text-right} Ilustrace vytvořená pomocí AI. ::: ::: ::: {.column width="60%"} ::: {.panel-result} [95% interval spolehlivosti (*confidence interval*, CI)]{.panel-title} Skutečnost je pevná jako kolík; každý výběr dá jiný interval jako jiný hod kroužkem. Při mnoha srovnatelných opakováních zachytí skutečný sklon přibližně 95 % takto vytvořených intervalů. Tomuto podílu říkáme **pokrytí** (*coverage*). ::: ::: {.panel-warning} [Neříkáme]{.panel-title} - 95 % pozorování je uvnitř - tento interval má 95% šanci obsahovat pravdu - interval dokazuje určitou hodnotu ::: ::: :::: ::: {.notes} Kolík stojí pevně jako skutečný sklon, kroužek se při každém hodu posune jako interval. Většina hodů kolík obepne, některé ne, a z jednoho hodu nepoznáme, jestli patří k těm úspěšným. Interval nepopisuje rozdělení pozorování, ale nejistotu sklonu. Jednomu už vypočtenému intervalu nepřisuzujeme pravděpodobnost; číslo 95 % patří k postupu. Interval také nic nedokazuje, jen ukazuje sklony slučitelné s daty a modelem. ::: ## Jak takový interval vytvořit z jednoho výběru? ```{r} #| label: pripravit-vypocet-intervalu #| 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 slide shows three decimals; students must be able to redo that # arithmetic with the displayed numbers and reach 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) ) ``` ::: {.text-center} Máme [odhad `r format_cz(odhad_sklonu, digits = 3)`]{.text-hgl-indigo_velvet} a [SE `r format_cz(standardni_chyba_sklonu, digits = 3)`]{.text-hgl-orange} (na tři desetinná místa). Kolik standardních chyb přidat na každou stranu? ::: ::: {.fragment} $$ \text{polovina šířky} = \text{kritická hodnota} \times \text{SE} = `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)` $$ ::: ::: {.fragment} $$ \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)` $$ $$ \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)` $$ ::: ::: {.fragment .takeaway-strip} Pro 95% interval je **kritická hodnota** blízko 2; její přesnou hodnotu vysvětlíme v následující lekci. ::: ::: {.notes} Úplný postup slovy, čísly a symboly je ve výukových materiálech. Kritická hodnota závisí na hladině a na počtu pozorování. ::: ## Stejný interval v R
:::: {.columns .align-center} ::: {.column width="58%"} ```{r} #| label: zobrazit-confint #| echo: true # 95% intervaly spolehlivosti koeficientů confint(object = mod_gejzir) ``` ::: ::: {.column width="42%"} ::: {.panel-result} [Řádek `delka_erupce_min`]{.panel-title} [`r format_cz(dolni_interval)` až `r format_cz(horni_interval)`]{.text-hgl-amethyst} minuty čekání na minutu erupce — stejné hranice jako ručně. ::: ::: ::::
::: {.takeaway-strip} Sloupce `2.5 %` a `97.5 %`: 95% interval nechává na každé straně 2,5 %. ::: ## 95% interval spolehlivosti ```{r} #| label: pripravit-interval-95 #| include: false plot_interval_95 <- ggplot2::ggplot() + ggplot2::annotate( geom = "errorbarh", xmin = dolni_interval, xmax = horni_interval, y = 1, height = 0.18, color = biostat_cols[["amethyst"]], linewidth = 2.4 ) + ggplot2::annotate( geom = "point", x = odhad_sklonu, y = 1, color = biostat_cols[["indigo_velvet"]], size = 6 ) + ggplot2::annotate( geom = "text", x = c(dolni_interval, odhad_sklonu, horni_interval), y = c(0.82, 1.18, 0.82), label = c("dolní hranice", "odhad", "horní hranice"), hjust = c(0, 0.5, 1), color = biostat_cols[["graphite"]], size = 7 ) + ggplot2::scale_y_continuous( limits = c(0.65, 1.3), breaks = NULL ) + ggplot2::labs( x = "Sklon (min čekání / min erupce)", y = NULL ) + theme_biostat(base_size = 22) save_local_figure( plot = plot_interval_95, filename = "interval_95.png", height = 520 ) ``` ```{r} #| label: zobrazit-interval-95 #| echo: false #| fig-alt: "Vodorovný fialový interval sklonu s tmavým bodem odhadu uprostřed a popisky dolní hranice, odhadu a horní hranice." include_local_figure("interval_95.png") ``` ::: {.fragment} Pro Old Faithful: [**[`r format_cz(dolni_interval)`; `r format_cz(horni_interval)`]**]{.text-hgl-amethyst} minuty čekání na minutu erupce. ::: ::: {.fragment .takeaway-strip} Interval popisuje nejistotu **sklonu**, ne rozmezí čekání po jedné budoucí erupci — jednotlivá čekání se navíc rozptylují kolem přímky. ::: ## Statistická redakce: opravte titulek
::: {.panel-olive} [Dvojice — 90 sekund]{.panel-title} Redakce napsala: > „Náš interval má 95% šanci obsahovat skutečný sklon.“ Přepište titulek do **jedné přesné věty o postupu**. Nejvýše 15 slov. :::
::: {.fragment} ::: {.panel-result} [Jedna možná oprava]{.panel-title} „Při mnoha opakováních zachytí takto vytvořené intervaly skutečný sklon přibližně v 95 % případů.“ ::: ::: ::: {.notes} Nechám dvě nebo tři dvojice přečíst svou verzi, teprve potom ukážu možnou opravu a porovnám ji s odpovědí C. Správná věta mluví o postupu a opakování, ne o jednom už vypočteném intervalu. ::: ## Co se stane s 80% intervaly? ::: {.slide-margin-top-15} ::: {.card-question} ::: {.card-title} Otázka ::: Stejných `r pocet_hlaseni` hlášení, ale místo 95% intervalů spočítáme 80%. Jaké budou intervaly a kolik jich zachytí skutečný sklon? :::
:::: {.columns} ::: {.column width="50%"} ::: {.card-question .card-compact} [**A)** Užší, a zachytí ho méně z nich.]{.rn-circle-orange} ::: ::: ::: {.column width="50%"} ::: {.card-question .card-compact} **B)** Širší, a zachytí ho více z nich. ::: ::: ::::
::: {.takeaway-strip} Ve dvojici zvolte A, nebo B a připravte jedno vysvětlení. ::: ::: ::: {.notes} Po diskuzi zakroužkuji A. Nižší hladina spolehlivosti používá menší kritickou hodnotu, a proto přidá na každou stranu méně standardních chyb. ::: ## Sto pozorovatelů, nižší hladina spolehlivosti ```{r} #| label: pripravit-sto-intervalu-80-50 #| include: false pocet_pokrytych_80 <- data_pocty_pokrytych$pocet[data_pocty_pokrytych$hladina == 0.80] pocet_pokrytych_50 <- data_pocty_pokrytych$pocet[data_pocty_pokrytych$hladina == 0.50] purrr::walk( .x = c(0.80, 0.50), .f = ~ save_local_figure( plot = plot_coverage( data_coverage = data_pokryti_intervalu, level = .x, slope_reference = referencni_sklon, xlim = rozsah_osy_intervalu ), filename = paste0("sto_intervalu_pokryti_", 100 * .x, ".png"), height = 700 ) ) ``` ::: {.r-stack} ```{r} #| label: zobrazit-sto-intervalu-80 #| echo: false #| out-width: "84%" #| fig-alt: "Stejných sto hlášení s osmdesátiprocentními intervaly: intervaly jsou kratší a více z nich je šedých." include_local_figure("sto_intervalu_pokryti_80.png") ``` ::: {.fragment fragment-index=1} ```{r} #| label: zobrazit-sto-intervalu-50 #| echo: false #| out-width: "84%" #| fig-alt: "Stejných sto hlášení s padesátiprocentními intervaly: intervaly jsou krátké a zhruba polovina z nich referenční čáru míjí." include_local_figure("sto_intervalu_pokryti_50.png") ``` ::: ::: ::: {.takeaway-strip} 80 %: [**`r pocet_pokrytych_80` ze `r pocet_hlaseni`**]{.text-hgl-orange} intervalů zachytilo referenční sklon. [50 %: **`r pocet_pokrytych_50` ze `r pocet_hlaseni`**.]{.fragment fragment-index=1} [Hladina spolehlivosti popisuje úspěšnost postupu.]{.fragment fragment-index=2} ::: ::: {.notes} Stejná hlášení, jen jiná hladina spolehlivosti. Před ukázáním 50 % se zeptám, kolik intervalů asi zachytí tmavou čáru. ::: ## Užší neznamená automaticky přesnější ```{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, 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") ) ) plot_intervaly <- 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.24, color = biostat_cols[["amethyst"]], linewidth = 2 ) + ggplot2::geom_point( color = biostat_cols[["indigo_velvet"]], size = 5 ) + ggplot2::labs( x = "Sklon (min čekání / min erupce)", y = NULL ) + theme_biostat(base_size = 20) save_local_figure( plot = plot_intervaly, filename = "porovnani_intervalu.png", height = 500 ) ``` ```{r} #| label: zobrazit-porovnani-intervalu #| echo: false #| fig-alt: "Čtyři vodorovné intervaly pod sebou se stejným středem; osmdesátiprocentní je nejkratší, devadesátidevítiprocentní nejdelší." include_local_figure("porovnani_intervalu.png") ``` ::: {.box-rule} Intervaly porovnáváme jako přesnější a méně přesné jen při [stejné hladině spolehlivosti]{.text-hgl-amethyst}. ::: ::: {.fragment} Snížením hladiny z 95 % na 80 % získáme [užší interval]{.text-hgl-amethyst}, ale z dat jsme [nezískali víc informace]{.text-hgl-orange}: data i SE zůstaly stejné. ::: # Jak přesně známe celou přímku? ## ::: {.slide-margin-top-15} ::: {.card-question} ::: {.card-title} Otázka ::: Průměrné čekání pro každou délku erupce, celý tento vztah, nazýváme [průměrný trend]{.text-hgl-indigo_velvet}; přímka ho odhaduje. Nakreslíme přímky všech `r pocet_hlaseni` hlášení GeyserWatch do jednoho grafu. Kde se od sebe budou lišit nejméně? :::
:::: {.columns} ::: {.column width="33%"} ::: {.card-question .card-compact} **A)** U krátkých erupcí. ::: ::: ::: {.column width="33%"} ::: {.card-question .card-compact} [**B)** U průměrně dlouhých.]{.rn-circle-orange} ::: ::: ::: {.column width="33%"} ::: {.card-question .card-compact} **C)** U dlouhých erupcí. ::: ::: ::::
::: {.takeaway-strip} Zvolte A, B, nebo C a zdůvodněte svou volbu sousedovi. ::: ::: ::: {.notes} Dosud jsme nejistotu vyjadřovali pro jediné číslo, sklon. Přímka ale dává pro každou délku erupce i odhadnuté průměrné čekání, a to je také odhad. To, jak se průměrné čekání mění s délkou erupce, nazýváme průměrný trend. Po diskuzi zakroužkuji B; následující slide ukáže proč. ::: ## Sto hlášení, sto přímek ```{r} #| label: pripravit-vejir-sta-primek #| include: false # The fan and the band share one x range and the mean eruption length. 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 ) plot_vejir <- 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.7 ) + 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.6 ) + 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 = 20) save_local_figure( plot = plot_vejir, filename = "vejir_sta_primek.png", height = 700 ) ``` ```{r} #| label: zobrazit-vejir-sta-primek #| echo: false #| fig-alt: "Sto tenkých fialových přímek nakreslených přes sebe s tmavou referenční přímkou. 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íř." include_local_figure("vejir_sta_primek.png") ``` ::: {.takeaway-strip} Fialové jsou přímky hlášení, [tmavá je přímka ze všech erupcí]{.text-hgl-graphite}. Přímky se kříží u průměrné délky erupce (čárkovaně) a k okrajům se rozbíhají. ::: ::: {.notes} Přímka z lm() vždy prochází bodem (průměrná délka erupce; průměrné čekání) svých dat, a tyto průměry se mezi hlášeními liší málo. Proto se přímky s různým sklonem kříží blízko středu a nejistota sklonu se projeví tím víc, čím dál jsme od průměrné délky. Hlášení mají jen `r pocet_erupci_v_hlaseni` erupcí, takže vějíř je širší než nejistota ze všech `r format_cz(nrow(data_gejzir))` erupcí, stejně jako SE sklonu v jednom hlášení byla větší. ::: ## Pás pro průměrný trend ```{r} #| label: pripravit-pas-spolehlivosti #| include: false # Band intervals annotated at a short and at the average eruption length delka_kratka <- 2 mat_dve_delky <- predict( object = mod_gejzir, newdata = data.frame(delka_erupce_min = c(delka_kratka, prumerna_delka)), interval = "confidence" ) data_dve_delky <- data.frame( delka_erupce_min = c(delka_kratka, prumerna_delka), lwr = mat_dve_delky[, "lwr"], upr = mat_dve_delky[, "upr"] ) plot_pas <- plot_data + ggplot2::geom_smooth( method = "lm", formula = y ~ x, se = TRUE, color = biostat_cols[["indigo_velvet"]], fill = biostat_cols[["amethyst"]], alpha = 0.25, linewidth = 1.6 ) + ggplot2::geom_vline( xintercept = prumerna_delka, linetype = "dashed", color = biostat_cols[["graphite"]] ) plot_pas_intervaly <- plot_pas + ggplot2::geom_errorbar( data = data_dve_delky, mapping = ggplot2::aes( x = delka_erupce_min, ymin = lwr, ymax = upr ), inherit.aes = FALSE, width = 0.12, linewidth = 1.3, color = biostat_cols[["indigo_velvet"]] ) + ggplot2::geom_label( data = data_dve_delky, mapping = ggplot2::aes( x = delka_erupce_min + 0.1, y = lwr - 3, label = paste0( "šířka ", formatC(x = upr - lwr, format = "f", digits = 1, decimal.mark = ","), " min" ) ), inherit.aes = FALSE, hjust = 0, size = 6, color = biostat_cols[["indigo_velvet"]], fill = scales::alpha(biostat_cols[["white"]], 0.75), label.size = 0 ) save_local_figure( plot = plot_pas_intervaly, filename = "pas_spolehlivosti_intervaly.png", height = 700 ) ``` ```{r} #| label: zobrazit-pas-spolehlivosti-intervaly #| echo: false #| out-width: "88%" #| fig-alt: "Bodový graf všech erupcí s fialovou přímkou a světle fialovým pásem kolem ní. Dvě svislé úsečky ukazují interval pro průměrné čekání: u dvouminutové erupce je širší, u průměrné délky erupce na přerušované čáře užší." include_local_figure("pas_spolehlivosti_intervaly.png") ``` ::: {.takeaway-strip} [Pás spolehlivosti]{.text-hgl-amethyst}: pro každou délku erupce 95% interval pro **průměrné** čekání; hranice vrátí `predict(..., interval = "confidence")`. Nejužší je u průměrné délky erupce (čárkovaně), ne tam, kde je nejvíc bodů. ::: ::: {.notes} Skutečný výzkumník sto hlášení nemá, ale nejistotu přímky odhadne z jednoho souboru dat, podobně jako SE sklonu. Průměrná délka erupce leží v řidší oblasti mezi shluky krátkých a dlouhých erupcí, takže rozhoduje vzdálenost od průměru, ne počet bodů v okolí. Pás je užší než vějíř z minulého slidu, protože vychází ze všech erupcí, ne ze čtyřiceti. Číslo 95 % čteme stejně jako u sklonu: patří k postupu. ::: ## Kdy je pás široký? ```{r} #| label: pripravit-galerii-pasu #| include: false list_galerie_pasu <- simulate_band_gallery() save_local_figure( plot = plot_band_gallery( list_gallery = list_galerie_pasu, base_size = 18 ), filename = "galerie_pasu.png", height = 820 ) ``` ```{r} #| label: zobrazit-galerii-pasu #| echo: false #| out-width: "90%" #| fig-alt: "Čtyři panely simulovaných dat se stejným skutečným trendem (přerušovaná čára). Fialová odhadnutá přímka se světle fialovým pásem: při mnoha bodech s malým rozptylem je pás tenký, při velkém rozptylu a při pěti bodech širší a při bodech jen mezi x = 4,5 a 5,5 se pás mimo data rozevírá přes celý graf." include_local_figure("galerie_pasu.png") ``` ::: {.fragment .takeaway-strip} Pás rozšiřují stejné ingredience jako SE: [velký rozptyl kolem přímky]{.text-hgl-orange}, [málo bodů]{.text-hgl-orange}, [úzký rozsah x]{.text-hgl-orange}. Nejširší je [daleko od dat]{.text-hgl-orange}: tam přímku prodlužujeme naslepo. ::: ::: {.notes} Simulovaná data se stejným skutečným trendem ve všech panelech. Nejdřív se zeptám, ve kterém panelu bude pás u x = 9 nejširší. Panel 4 má body jen mezi 4,5 a 5,5: u nich je pás úzký, ale mimo ně se rychle rozevírá a přímku tam nelze spolehlivě prodlužovat. Naše erupce pokrývají délky od `r format_cz(min(data_gejzir$delka_erupce_min))` do `r format_cz(max(data_gejzir$delka_erupce_min))` min, proto je pás v jejich rozsahu úzký. ::: ## Kolik erupcí leží uvnitř pásu? ```{r} #| label: pripravit-body-v-pasu #| include: false # Band limits at every observed eruption mat_pas_pozorovani <- predict( object = mod_gejzir, interval = "confidence" ) data_gejzir_pas <- data_gejzir |> dplyr::mutate( v_pasu = cekani_min >= mat_pas_pozorovani[, "lwr"] & cekani_min <= mat_pas_pozorovani[, "upr"] ) pocet_v_pasu <- sum(data_gejzir_pas$v_pasu) save_local_figure( plot = plot_pas, filename = "pas_spolehlivosti.png", height = 700 ) plot_pas_uvnitr <- plot_pas + ggplot2::geom_point( data = dplyr::filter(data_gejzir_pas, v_pasu), mapping = ggplot2::aes( x = delka_erupce_min, y = cekani_min ), color = biostat_cols[["orange"]], size = 3.4 ) save_local_figure( plot = plot_pas_uvnitr, filename = "pas_spolehlivosti_uvnitr.png", height = 700 ) ``` ::: {.r-stack} ```{r} #| label: zobrazit-pas-spolehlivosti #| echo: false #| out-width: "82%" #| fig-alt: "Bodový graf všech erupcí s fialovou přímkou a světle fialovým pásem spolehlivosti." include_local_figure("pas_spolehlivosti.png") ``` ::: {.fragment} ```{r} #| label: zobrazit-pas-spolehlivosti-uvnitr #| echo: false #| out-width: "82%" #| fig-alt: "Stejný graf; oranžově jsou zvýrazněné erupce uvnitř pásu. Je jich jen menšina, většina bodů leží mimo pás." include_local_figure("pas_spolehlivosti_uvnitr.png") ``` ::: ::: ::::: {.fragment} ::: {.panel-warning} [Pás ≠ kde leží 95 % erupcí]{.panel-title} Uvnitř leží jen [**`r pocet_v_pasu` ze `r format_cz(nrow(data_gejzir))`**]{.text-hgl-orange} erupcí. Pás ukazuje, kde může ležet průměrný trend. Rozsah pro jedno nové čekání musí být širší. ::: ::::: ::: {.notes} Nejprve nechám studenty odhadnout počet. To, že většina bodů leží mimo pás, samo o sobě není známkou špatného modelu. Jedno nové čekání se od trendu odchýlí o své vlastní rozptýlení, takže rozsah pro něj musí zahrnout dvě nejistoty: polohu trendu a rozptýlení kolem něj. Anglicky se mu říká prediction interval; jeho výpočet zde nepotřebujeme. ::: # Může jedna erupce změnit závěr? ## ```{r} #| label: pripravit-vyrazene-erupce #| include: false # The three eruptions removed during data preparation (GeyserTimes archive) data_vyrazene <- data.frame( 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), oznaceni = c("A", "B", "C") ) data_vyrazene$delka_erupce_min <- data_vyrazene$duration / 60 plot_vyrazene <- plot_pas + ggplot2::geom_point( data = data_vyrazene, color = biostat_cols[["orange"]], size = 5 ) + ggplot2::geom_text( data = data_vyrazene, mapping = ggplot2::aes(label = oznaceni), nudge_y = 6, fontface = "bold", size = 8, color = biostat_cols[["graphite"]] ) save_local_figure( plot = plot_vyrazene, filename = "vyrazene_erupce.png", width = 1200, height = 820 ) ``` ::: {.card-question} ::: {.card-title} Otázka ::: Při přípravě dat jsme [tři erupce vyřadili]{.text-hgl-orange}. Která by po vrácení naklonila přímku nejvíc? :::
:::: {.columns .align-center} ::: {.column width="60%"} ```{r} #| label: zobrazit-vyrazene-erupce #| echo: false #| fig-alt: "Bodový graf erupcí lekce s fialovou přímkou a pásem. Tři oranžové body: A mezi krátkými erupcemi nad přímkou, B mezi dlouhými erupcemi nad přímkou, C úplně vlevo, kratší než všechny ostatní erupce, s dlouhým čekáním." include_local_figure("vyrazene_erupce.png") ``` ::: ::: {.column width="40%"} ::: {.card-question .card-compact} **A)** `r formatC(x = data_vyrazene$delka_erupce_min[1], format = "f", digits = 2, decimal.mark = ",")` min erupce, čekání `r data_vyrazene$cekani_min[1]` min ::: ::: {.card-question .card-compact} **B)** `r formatC(x = data_vyrazene$delka_erupce_min[2], format = "f", digits = 2, decimal.mark = ",")` min erupce, čekání `r data_vyrazene$cekani_min[2]` min ::: ::: {.card-question .card-compact} [**C)** `r formatC(x = data_vyrazene$delka_erupce_min[3], format = "f", digits = 2, decimal.mark = ",")` min erupce, čekání `r data_vyrazene$cekani_min[3]` min]{.rn-circle-orange} ::: ::: :::: ::: {.notes} Hlasování A, B, nebo C se zdůvodněním sousedovi. Po diskuzi zakroužkuji C: leží nejdál od průměrné délky erupce a zároveň daleko od přímky. Erupce byly vyřazeny při přípravě dat; teď se podíváme, jestli to mělo vliv. ::: ## Krajní erupce mají delší páku ```{r} #| label: pripravit-paku-a-cooka #| include: false data_gejzir_120 <- rbind( data_gejzir[, c("delka_erupce_min", "cekani_min")], data_vyrazene[, c("delka_erupce_min", "cekani_min")] ) mod_gejzir_120 <- lm( formula = cekani_min ~ delka_erupce_min, data = data_gejzir_120 ) data_vliv_120 <- data_gejzir_120 |> dplyr::mutate( pakovy_efekt = hatvalues(model = mod_gejzir_120), residuum = resid(object = mod_gejzir_120), cookova_vzdalenost = cooks.distance(model = mod_gejzir_120), oznaceni = c(rep(NA_character_, nrow(data_gejzir)), data_vyrazene$oznaceni) ) plot_pakovy_efekt <- ggplot2::ggplot( data = data_vliv_120, mapping = ggplot2::aes( x = delka_erupce_min, y = pakovy_efekt ) ) + ggplot2::labs( x = "Délka erupce (min)", y = "Pákový efekt" ) + theme_biostat(base_size = 20) + ggplot2::geom_vline( xintercept = mean(data_gejzir_120$delka_erupce_min), linetype = "dashed", color = biostat_cols[["graphite"]] ) + ggplot2::geom_point( color = biostat_cols[["grey_olive"]], size = 3 ) + ggplot2::geom_point( data = dplyr::filter(data_vliv_120, !is.na(oznaceni)), color = biostat_cols[["orange"]], size = 5 ) + ggplot2::geom_text( data = dplyr::filter(data_vliv_120, !is.na(oznaceni)), mapping = ggplot2::aes(label = oznaceni), nudge_y = 0.006, fontface = "bold", size = 8, color = biostat_cols[["graphite"]] ) save_local_figure( plot = plot_pakovy_efekt, filename = "pakovy_efekt.png", width = 1100, height = 820 ) ``` :::: {.columns .align-center} ::: {.column width="55%"} ```{r} #| label: zobrazit-pakovy-efekt #| echo: false #| fig-alt: "Pákový efekt proti délce erupce. Body tvoří oblouk tvaru U s minimem u průměrné délky (přerušovaná čára). Oranžová erupce C úplně vlevo má zdaleka nejvyšší pákový efekt, A je výš než B." include_local_figure("pakovy_efekt.png") ``` ::: ::: {.column width="45%"} ::: {.panel-data} [Pákový efekt (*leverage*)]{.panel-title} Jak daleko leží délka erupce od průměrné délky, poměřeno rozprostřením všech délek. Na čekání nezávisí. ::: ::: {.text-size-smaller} ```{r} #| label: zobrazit-kod-paky #| echo: true #| eval: false # Model ze 117 erupcí a vrácených A, B, C mod_gejzir_120 <- lm( formula = cekani_min ~ delka_erupce_min, data = data_gejzir_120 ) # Pákový efekt každé erupce hatvalues(model = mod_gejzir_120) ``` ::: ::: :::: ::: {.fragment .takeaway-strip} [Krajní erupce]{.text-hgl-orange} je jako dítě na konci houpačky: přímku nakloní nejsnáz. ::: ::: {.notes} Pákový efekt roste se vzdáleností od průměrné délky, stejně jako se rozšiřuje pás spolehlivosti. Je to druhá ingredience SE, rozprostření délek, jen pro každou erupci zvlášť. Erupce C má největší pákový efekt jen proto, že trvala nejkratší dobu. ::: ## Kdy jeden bod pohne přímkou? ```{r} #| label: pripravit-galerii-vlivu #| include: false list_galerie_vlivu <- simulate_influence_gallery() save_local_figure( plot = plot_influence_gallery( list_gallery = list_galerie_vlivu, base_size = 18 ), filename = "galerie_vlivu.png", height = 800 ) ``` ```{r} #| label: zobrazit-galerii-vlivu #| echo: false #| out-width: "90%" #| fig-alt: "Čtyři panely simulovaných dat se stejnými šedými body a jedním přidaným oranžovým bodem. Plná fialová přímka je odhadnutá s bodem, přerušovaná bez něj. Jen v panelu s velkou pákou a velkým residuem se přímky výrazně rozcházejí." include_local_figure("galerie_vlivu.png") ``` ::: {.fragment .takeaway-strip} [Cookova vzdálenost]{.text-hgl-orange}: jak moc by se změnily odhadnuté hodnoty bez bodu. Velká je až při velké páce a velkém residuu současně. ::: ::: {.notes} Simulovaná data: stejné šedé body, v každém panelu jeden přidaný oranžový bod. Plná přímka je odhadnutá s ním, přerušovaná bez něj. V záhlaví je pákový efekt a Cookova vzdálenost přidaného bodu. Nejdřív se zeptám, ve kterém panelu se přímky rozejdou nejvíc. ::: ## Delší páka a velké residuum, velký vliv ```{r} #| label: pripravit-tabulku-vlivu #| include: false table_vliv_abc <- data_vliv_120 |> dplyr::filter(!is.na(oznaceni)) |> dplyr::mutate( Erupce = oznaceni, `Residuum (min)` = format_cz(residuum, digits = 1), `Pákový efekt` = format_cz(pakovy_efekt, digits = 3), `Cookova vzdálenost` = format_cz(cookova_vzdalenost, digits = 3) ) |> dplyr::select( Erupce, `Residuum (min)`, `Pákový efekt`, `Cookova vzdálenost` ) ``` :::: {.columns .align-center} ::: {.column width="40%"} ::: {.panel-data} [Cookova vzdálenost (*Cook's distance*)]{.panel-title} Jak moc by se změnily odhadnuté hodnoty, kdyby erupce chyběla. Velká je, když se dlouhá páka potká s velkým residuem. Takovému pozorování říkáme **vlivné pozorování**. ::: ::: {.text-size-smaller} ```{r} #| label: zobrazit-kod-cooka #| echo: true #| eval: false # Cookova vzdálenost každé erupce cooks.distance(model = mod_gejzir_120) ``` ::: ::: ::: {.column width="60%" .text-size-smaller} ::: {.fragment} ```{r} #| label: zobrazit-tabulku-vlivu #| echo: false knitr::kable( x = table_vliv_abc, align = "crrr" ) ``` ::: ::: :::: ::: {.fragment .takeaway-strip} A a B jsou odlehlé se [stejně velkým residuem]{.text-hgl-orange}; A má delší páku, a proto větší Cookovu vzdálenost. Výrazně vlivná je jen C: spojuje nejdelší páku s největším residuem. ::: ::: {.notes} Erupce s dlouhou pákou, která leží přímo na přímce, by sklon nezměnila. Vliv vzniká až spojením páky s residuem. ::: ## Mění některá z našich erupcí závěr? ```{r} #| label: pripravit-vliv-erupci-lekce #| include: false hranice_cook_lekce <- 4 / nrow(data_gejzir) pocet_nad_hranici_lekce <- sum(cooks.distance(model = mod_gejzir) > hranice_cook_lekce) data_cook_lekce <- tibble::tibble( poradi = seq_len(nrow(data_gejzir)), cookova_vzdalenost = cooks.distance(model = mod_gejzir) ) max_zmena_sklonu_lekce <- max(abs(dfbeta(model = mod_gejzir)[, "delka_erupce_min"])) plot_cook_lekce <- ggplot2::ggplot( data = data_cook_lekce, mapping = ggplot2::aes( x = poradi, y = cookova_vzdalenost ) ) + ggplot2::labs( x = "Pořadí erupce v datech", y = "Cookova vzdálenost" ) + theme_biostat(base_size = 20) + ggplot2::geom_segment( mapping = ggplot2::aes(xend = poradi, yend = 0), color = biostat_cols[["grey_olive"]], linewidth = 1 ) + ggplot2::geom_hline( yintercept = hranice_cook_lekce, linetype = "dashed", color = biostat_cols[["indigo_velvet"]], linewidth = 1 ) + ggplot2::annotate( geom = "text", x = nrow(data_gejzir), y = hranice_cook_lekce, label = "4/n", vjust = -0.6, hjust = 1, size = 7, color = biostat_cols[["indigo_velvet"]] ) save_local_figure( plot = plot_cook_lekce, filename = "cookova_vzdalenost_lekce.png", width = 1100, height = 820 ) ``` :::: {.columns .align-center} ::: {.column width="52%"} ```{r} #| label: zobrazit-cook-lekce #| echo: false #| fig-alt: "Cookova vzdálenost erupcí lekce jako svislé čárky s přerušovanou čarou 4/n. Několik čárek čáru mírně přesahuje, žádná výrazně nevyčnívá." include_local_figure("cookova_vzdalenost_lekce.png") ``` ::: ::: {.column width="48%"} ::: {.panel-data} [Čára 4/n]{.panel-title} 4 děleno počtem erupcí: orientační hodnota, ne hranice pro mazání. V modelu ze `r nrow(data_gejzir)` erupcí ji přesahuje `r pocet_nad_hranici_lekce` z `r nrow(data_gejzir)` erupcí: na ty se podíváme blíž. ::: ::::: {.fragment} ::: {.panel-result} [Vynecháme jednu erupci]{.panel-title} Sklon se změní nejvýš o [**`r format_cz(max_zmena_sklonu_lekce)` min/min**]{.text-hgl-indigo_velvet}, tedy asi o `r format_cz(100 * max_zmena_sklonu_lekce / standardni_chyba_sklonu, digits = 0)` % SE ([`r format_cz(standardni_chyba_sklonu)`]{.text-hgl-orange}). ::: ::::: ::: :::: ::: {.fragment .takeaway-strip} Závěr [nestojí na jediné erupci]{.text-hgl-orange}: vliv jedné erupce porovnáváme se standardní chybou. ::: ::: {.notes} Změnu sklonu po vynechání každé jednotlivé erupce vrací funkce dfbeta(). Erupce nad čarou 4/n tím prověříme všechny najednou. ::: ## Co s vlivným pozorováním? ```{r} #| label: pripravit-model-bez-c #| include: false # The most influential added eruption is C, the last row. index_nejvlivnejsi <- which.max(cooks.distance(model = mod_gejzir_120)) stopifnot( index_nejvlivnejsi == nrow(data_gejzir_120) ) sklon_bez_c <- coef( object = lm( formula = cekani_min ~ delka_erupce_min, data = data_gejzir_120[-index_nejvlivnejsi, ] ) )[["delka_erupce_min"]] sklon_120 <- coef(object = mod_gejzir_120)[["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 = "," ) ``` :::: {.columns .align-center} ::: {.column width="50%"} ::: {.panel-warning} [Prověřit, ne mazat]{.panel-title} - Velká Cookova vzdálenost je důvod podívat se blíž, ne důvod mazat. - O vyřazení rozhoduje znalost dat: chyba záznamu nebo událost mimo proces, který popisujeme. - Když si nejsme jistí, ukážeme výsledek s erupcí i bez ní. ::: ::: ::: {.column width="50%"} ::::: {.fragment} ::: {.panel-data} [Naše tři erupce]{.panel-title} - **A a B** posunou sklon jen o `r text_zmena_dvou_erupci` min/min (z `r format_cz(odhad_sklonu)` na `r format_cz(sklon_bez_c)`). Vyřadili jsme je jen kvůli přehlednosti a po ověření bychom je ponechali. - **C** přidaná k nim posune sklon dál na `r format_cz(sklon_120)` (asi `r format_cz((sklon_bez_c - sklon_120) / standardni_chyba_sklonu, digits = 1)` × SE našeho modelu). Trvala jen `r data_vyrazene$duration[3]` s, tedy pod rozsahem 1,5–5 min podle National Park Service (NPS): důvod k vyřazení máme. ::: ::::: ::: :::: ::: {.notes} Pro C je důvod stejného druhu jako u čekání, kde jsme ponechali jen čekání 50–127 minut uváděné NPS. Erupce A a B navíc posunou sklon jen o `r format_cz(abs(sklon_bez_c - odhad_sklonu))`, mnohem méně, než je SE. ::: # Závěr ## Co už víme o čekání na další erupci? ```{r} #| label: pripravit-ctyri-minuty #| include: false # Mean wait after a 4-minute eruption with its confidence interval vec_ctyri_minuty_text <- formatC( x = predict( object = mod_gejzir, newdata = data.frame(delka_erupce_min = 4), interval = "confidence" )[1, ], format = "f", digits = 2, decimal.mark = "," ) ``` ::: {.box-question} **O kolik déle čekáme po delší erupci — a jak jistě to víme?** :::
:::: {.columns} ::: {.column width="48%"} ::: {.panel-result} [Co už můžeme říct]{.panel-title} Dvě erupce, jejichž délky se liší o minutu, se v průměru liší čekáním o [`r format_cz(odhad_sklonu)` min]{.text-hgl-indigo_velvet}; s daty a modelem jsou slučitelné sklony od [`r format_cz(dolni_interval)` do `r format_cz(horni_interval)`]{.text-hgl-amethyst}. Po čtyřminutové erupci je odhadnuté průměrné čekání [`r vec_ctyri_minuty_text[["fit"]]` min (95% CI [`r vec_ctyri_minuty_text[["lwr"]]`; `r vec_ctyri_minuty_text[["upr"]]`])]{.text-hgl-amethyst}. ::: ::: ::: {.column width="48%"} ::: {.panel-question} [Co zůstává otevřené]{.panel-title} Jak pomocí dat posoudit [konkrétní tvrzení o sklonu]{.text-hgl-orange} — například že délka erupce s čekáním vůbec nesouvisí? ::: ::: :::: ## Napište jednu vědeckou větu :::: {.columns} ::: {.column width="47%"} ::: {.panel-olive} [Dvojice — 3 min]{.panel-title} Věta musí obsahovat: - vztah mezi délkou erupce a čekáním - odhad sklonu - SE a 95% CI - jednotky - opatrný jazyk bez tvrzení o příčině ::: ::: ::: {.column width="49%"} ::: {.panel-data} [Data a výsledek modelu]{.panel-title} - **`r format_cz(nrow(data_gejzir))` erupcí** gejzíru Old Faithful, léto 2024 - délka erupce: `r format_cz(min(data_gejzir$delka_erupce_min))`–`r format_cz(max(data_gejzir$delka_erupce_min))` min - čekání: `r format_cz(min(data_gejzir$cekani_min))`–`r format_cz(max(data_gejzir$cekani_min))` min - 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)`]** - jednotky: min čekání / min erupce ::: ::: :::: ::: {.takeaway-strip .fragment} Vyměňte si větu se sousední dvojicí: jsou v ní proměnné, jednotky a nejistota? ::: ## Jedna možná věta ::: {.slide-margin-top-15} ::: {.text-size-large .text-center} Dvě erupce, jejichž délky se lišily o jednu minutu, se podle lineárního modelu v průměru lišily čekáním na další erupci o [**`r format_cz(odhad_sklonu)`**]{.text-hgl-indigo_velvet} minuty (odhad sklonu = `r format_cz(odhad_sklonu)`, [SE = `r format_cz(standardni_chyba_sklonu)`]{.text-hgl-orange}, [95% CI [`r format_cz(dolni_interval)`; `r format_cz(horni_interval)`]]{.text-hgl-amethyst} minuty čekání na minutu erupce). ::: ::: ## Co ještě omezuje náš závěr? :::: {.columns .align-center} ::: {.column width="46%"} ::: {.text-size-large} Dobrovolníci nezaznamenali [náhodný výběr všech erupcí]{.text-hgl-orange} a blízké erupce mohou být časově propojené. :::
::: {.fragment} Výsledek je opatrný odhad vztahu v zachyceném období — ne neměnný přírodní zákon. Nestojí ale na jediné erupci. ::: ::: ::: {.column width="54%"} ```{r} #| label: zobrazit-nocni-erupci #| echo: false #| fig-alt: "Humorná ilustrace nočního výbuchu gejzíru Old Faithful: bizon a svišť na lavičce spí, stopky a prázdný zápisník jim padají z tlap a erupci sleduje jen výr. Ilustrace vytvořená pomocí AI." #| out-width: "100%" include_local_figure("geyser_night_eruption.png") ``` ::: {.text-size-tiny .text-right} Ilustrace vytvořená pomocí AI. ::: ::: :::: ::: {.notes} Některé erupce, třeba v noci, nikdo nezaznamená. Zaznamenané erupce proto nejsou náhodným výběrem všech erupcí. ::: ## Shrnutí ::: {.slide-margin-top-15} ::: {.incremental} - [Jiný výběr]{.text-hgl-graphite} dává jiný odhad sklonu: výběrová variabilita. - [`Std. Error`]{.text-hgl-orange} je přesnost odhadu: roste s residui, klesá s rozsahem prediktoru a s odmocninou z počtu pozorování. - [95% interval]{.text-hgl-amethyst} a [pás spolehlivosti]{.text-hgl-amethyst}: 95 % patří k postupu, ne k pozorováním. - [Vlivné pozorování]{.text-hgl-orange} = dlouhá páka a velké residuum; prověřit, ne mazat. - Vědecká věta: odhad, SE, interval a jednotky. ::: ::: ::: {.notes} Shrnutí kopíruje výsledky učení. ::: ## Co potřebujete vidět, než napíšete „po minutě erupce navíc čekáme o `r format_cz(odhad_sklonu)` minuty déle“? {.text-center}
::: {.text-size-large .text-bold .text-center} Jaké tři otázky položíte autorovi, který vám ukáže jen tento titulek? :::
::: {.panel-olive} [Řekněte sousedovi jednou větou]{.panel-title} Kterou z otázek byste položili jako první, a proč? ::: ::: {.fragment .takeaway-strip} Jak posoudit konkrétní tvrzení o sklonu? Následující lekce přidá testování hypotéz. ::: ::: {.notes} Otevřená otázka bez jedné správné odpovědi. Očekávané otázky: Jak vypadá graf s pásem spolehlivosti a residua? Kolik erupcí a jak byly zaznamenány? Jaká je SE a interval? Nemění závěr jedna erupce, a vyřadil autor některé erupce? Mění se sklon mezi obdobími? Tvrdí autor příčinu, nebo jen průměrný rozdíl? :::